vertInterpolation_pressKappa Subroutine

public subroutine vertInterpolation_pressKappa(fModel, fPres, ps, plevs, undef, PHIS, phis_units, rc)

The subroutine vertInterpolation_pressKappa performs a linear vertical interpolation from pressure levels to model levels in the pressure kappa formulation. It is assumed that we have a top to bottom ordering of the vertical levels.

Arguments

Type IntentOptional Attributes Name
type(ESMF_Field), intent(inout) :: fModel
type(ESMF_Field), intent(inout) :: fPres
type(ESMF_Field), intent(inout) :: ps
real, intent(in) :: plevs(:)

pressure levels
real, intent(in) :: undef

undefined value
type(ESMF_Field), intent(inout), optional :: PHIS

Surface Geopotential Heights
character(len=*), intent(in), optional :: phis_units

units of phis
integer, intent(out), optional :: rc

Calls

proc~~vertinterpolation_presskappa~~CallsGraph proc~vertinterpolation_presskappa vertInterpolation_pressKappa ESMF_AttributeGet ESMF_AttributeGet proc~vertinterpolation_presskappa->ESMF_AttributeGet ESMF_UtilStringLowerCase ESMF_UtilStringLowerCase proc~vertinterpolation_presskappa->ESMF_UtilStringLowerCase esmf_fieldget esmf_fieldget proc~vertinterpolation_presskappa->esmf_fieldget interface~mapl_assert MAPL_Assert proc~vertinterpolation_presskappa->interface~mapl_assert proc~mapl_gridget MAPL_GridGet proc~vertinterpolation_presskappa->proc~mapl_gridget proc~mapl_return MAPL_Return proc~vertinterpolation_presskappa->proc~mapl_return proc~mapl_verify MAPL_Verify proc~vertinterpolation_presskappa->proc~mapl_verify proc~mapl_gridget->ESMF_AttributeGet proc~mapl_gridget->proc~mapl_return proc~mapl_gridget->proc~mapl_verify ESMF_DistGridGet ESMF_DistGridGet proc~mapl_gridget->ESMF_DistGridGet ESMF_GridGet ESMF_GridGet proc~mapl_gridget->ESMF_GridGet proc~mapl_distgridget MAPL_DistGridGet proc~mapl_gridget->proc~mapl_distgridget proc~mapl_getimsjms MAPL_GetImsJms proc~mapl_gridget->proc~mapl_getimsjms proc~mapl_gridhasde MAPL_GridHasDE proc~mapl_gridget->proc~mapl_gridhasde at at proc~mapl_return->at insert insert proc~mapl_return->insert proc~mapl_throw_exception MAPL_throw_exception proc~mapl_return->proc~mapl_throw_exception proc~mapl_verify->proc~mapl_throw_exception proc~mapl_distgridget->proc~mapl_verify proc~mapl_distgridget->ESMF_DistGridGet proc~mapl_getimsjms->interface~mapl_assert proc~mapl_getimsjms->proc~mapl_return proc~mapl_getimsjms->proc~mapl_verify interface~mapl_sort MAPL_Sort proc~mapl_getimsjms->interface~mapl_sort proc~mapl_gridhasde->proc~mapl_return proc~mapl_gridhasde->proc~mapl_verify proc~mapl_gridhasde->ESMF_DistGridGet proc~mapl_gridhasde->ESMF_GridGet ESMF_DELayoutGet ESMF_DELayoutGet proc~mapl_gridhasde->ESMF_DELayoutGet

Source Code

      subroutine vertInterpolation_pressKappa (fModel, fPres, ps, plevs, &
                                  undef, phis, phis_units, rc)
!
      type(ESMF_Field),           intent(INOUT) :: fModel  
      type(ESMF_Field),           intent(INOUT) :: fPres  
      type(ESMF_Field),           intent(INOUT) :: ps  
      real,                       intent(IN)    :: plevs(:)   !! pressure levels
      real,                       intent(IN)    :: undef      !! undefined value
      type(ESMF_Field), optional, intent(INOUT) :: PHIS       !! Surface Geopotential Heights
      character(len=*), optional, intent(in)    :: phis_units !! units of phis
      integer, optional,          intent(OUT)   :: rc
!
      integer :: im, jm, lmMod, lmPres
      integer :: i, j, L, k, k2, k3, dims(3)
      integer, allocatable :: kbeg(:,:),kend(:,:)
      real    :: vdef, vsurf, TH
      real,allocatable :: pl_mod(:,:,:),ple_mod(:,:,:)
      character(len=1) :: vartype
      real(REAL64), allocatable :: ak(:),bk(:)
      integer :: status
      real    :: gfactor
      type(ESMF_Grid) :: grid
      real, pointer   :: vMod(:,:,:), vPres(:,:,:), vPS(:,:), vPHIS(:,:)
      character(ESMF_MAXSTR) :: vname, units
!
!------------------------------------------------------------------------------
      ! get dimensions, allocate
      call ESMF_FieldGet(fModel,grid=grid,rc=status)
      _VERIFY(STATUS)
      call ESMF_AttributeGet(fModel,name='UNITS',value=units,rc=status)
      _VERIFY(STATUS)
      call ESMF_AttributeGet(fModel,name='LONG_NAME',value=vname,rc=status)
      _VERIFY(STATUS)
      vname = ESMF_UtilStringLowerCase(vname,rc=status)
      call MAPL_GridGet(grid, localCellCountPerDim=dims,rc=status)
      _VERIFY(STATUS)
      im = dims(1)
      jm = dims(2)
      lmMod = dims(3)
      allocate(ak(lmMod+1),stat=status)
      _VERIFY(STATUS)
      allocate(bk(lmMod+1),stat=status)
      _VERIFY(STATUS)
      allocate(pl_mod(im,jm,lmMod),stat=status)
      _VERIFY(STATUS)
      allocate(ple_mod(im,jm,lmMod+1),stat=status)
      _VERIFY(STATUS)
      allocate(kbeg(im,jm),stat=status)
      _VERIFY(STATUS)
      allocate(kend(im,jm),stat=status)
      _VERIFY(STATUS)

      call ESMF_FieldGet(fPres,grid=grid,rc=status)
      _VERIFY(STATUS)
      call MAPL_GridGet(grid, localCellCountPerDim=dims,rc=status)
      _VERIFY(STATUS)
      lmPres = dims(3)

      ! given PS, get AK, BK from the grid to get ple
      call ESMF_FieldGet(PS,0,farrayPtr=vPS,rc=status)
      _VERIFY(STATUS)
      call ESMF_FieldGet(PS,grid=grid,rc=status)
      _VERIFY(STATUS)
      call ESMF_AttributeGet(grid,name="GridAK",valuelist=ak,rc=status)
      _VERIFY(STATUS)
      call ESMF_AttributeGet(grid,name="GridBK",valuelist=bk,rc=status)
      _VERIFY(STATUS)
      do i=1,lmmod+1
         ple_mod(:,:,i)=ak(i)+bk(i)*vPS(:,:)
      enddo
      do i=1,lmmod
         pl_mod(:,:,i)=(ple_mod(:,:,i)+ple_mod(:,:,i+1))/2.0
      enddo
      
      call ESMF_FieldGet(fModel,0,farrayPtr=vMod,rc=status)
      _VERIFY(STATUS)
      call ESMF_FieldGet(fPres,0,farrayPtr=vPres,rc=status)
      _VERIFY(STATUS)

      !---------------------------------------------------------------
      ! For each grid point (i,j), determine the range of model levels
      ! (kbeg to kend) that will be used for the main interpolation
      ! formula.
      !---------------------------------------------------------------
      call searchVerticalLocations(pl_mod, plevs, kbeg, kend, im, jm, lmmod, lmpres)

      !-------------------------------
      ! Compute trapped interpolations
      !-------------------------------
    
      DO j = 1, jm
         DO i = 1, im
            DO L = kbeg(i,j), kend(i,j)
               SEARCH_LEVELS : DO k = 2, lmpres
                  IF (plevs(k)   .GE. pl_mod(i,j,L) .AND. &
                      plevs(k-1) .LT. pl_mod(i,j,L) ) THEN
                     IF ( (vpres(i,j,k) .NE. undef) .AND. (vpres(i,j,k-1) .NE. undef) ) THEN
                        ! Apply the pressure to the kappa linear interpolation
                        vmod(i,j,L) = interpVal(vpres(i,j,k), vpres(i,j,k-1), pl_mod(i,j,L), &
                                                plevs(k), plevs(k-1), MAPL_KAPPA)
                     ELSEIF ( (vpres(i,j,k) .NE. undef) .AND. (vpres(i,j,k-1) .EQ. undef) ) THEN
                        vmod(i,j,L) = vpres(i,j,k)
                     ELSEIF ( (vpres(i,j,k) .EQ. undef) .AND. (vpres(i,j,k-1) .NE. undef) ) THEN
                        vmod(i,j,L) = vpres(i,j,k-1)
                     ELSEIF ( (vpres(i,j,k) .EQ. undef) .AND. (vpres(i,j,k-1) .EQ. undef) ) THEN
                        ! Search for the first non undefined value
                        vdef = undef
                        ! First: search by going toward the surface
                        SEARCH_DEF3: DO k3 = k+1,lmpres
                           IF (vpres(i,j,k3) .NE. undef) THEN
                              vdef = vpres(i,j,k3)
                              EXIT SEARCH_DEF3
                           END IF
                        END DO SEARCH_DEF3

                        IF (vdef .EQ. undef) THEN
                        ! Second: search by going toward the top
                           SEARCH_DEF4: DO k3 = k-2, 1, -1
                              IF (vpres(i,j,k3) .NE. undef) THEN
                                 vdef = vpres(i,j,k3)
                                 EXIT SEARCH_DEF4
                              END IF
                           END DO SEARCH_DEF4
                        END IF
                        if (vdef == undef) then
                           _RETURN(ESMF_FAILURE)
                        end if
                        vmod(i,j,L) = vdef
                     END IF
                     EXIT SEARCH_LEVELS
                  END IF
               ENDDO SEARCH_LEVELS
            ENDDO
         ENDDO
      ENDDO

      !--------------------------------------
      ! Do extrapolation on levels below kbeg
      ! (toward the top of the atmosphere)
      !--------------------------------------
      DO j = 1, jm
         DO i = 1, im
            vdef = vpres(i,j,1)
            IF (vdef .EQ. undef) THEN
               ! Search for the first non undefined value
               SEARCH_DEF1: DO k = 2, lmpres
                  IF (vpres(i,j,k) .NE. undef) THEN
                     vdef = vpres(i,j,k)
                     EXIT SEARCH_DEF1
                  END IF
               END DO SEARCH_DEF1
            END IF
            DO L = 1, kbeg(i,j) -1
               vmod(i,j,L) = vdef
            ENDDO
         ENDDO
      ENDDO
      
      !--------------------------------------
      ! Do extrapolation on levels above kend
      ! (toward the surface) get var type
      !--------------------------------------
      ! = 'W' for winds ( U & V)
      ! = 'T' for temperature
      ! = 'H' for height                        
      ! = 'O' for any other variable
      if (index(vname,"wind")/=0) then
         varType = 'W'
      else if (index(vname,"temperature")/=0) then
         varType = 'T'
      else if (index(vname,"height")/=0) then
         varType = 'H'
         _ASSERT(present(PHIS),'needs informative message')
         call ESMF_FieldGet(PHIS,0,farrayPtr=vPHIS,rc=status)
         _VERIFY(STATUS)
         _ASSERT(present(phis_units),'needs informative message')
         if (trim(units)=="m" .and. trim(phis_units)=="m+2 s-2") then
            gfactor = MAPL_GRAV
         else
            gfactor = 1.0
         end if   
      else
         varType = 'O'
      end if                  
      DO j = 1, jm
         DO i = 1, im
            vdef = vpres(i,j,lmpres)
            IF (vdef .EQ. undef) THEN
               ! Search for the first non undefined value
               SEARCH_DEF2: DO k2 = lmpres-1, 1, -1
                  IF (vpres(i,j,k2) .NE. undef) THEN
                     vdef = vpres(i,j,k2)
                     EXIT SEARCH_DEF2
                  END IF
               END DO SEARCH_DEF2
               DO L = kend(i,j)+1, lmmod
                  IF (varType .EQ. 'W') THEN
                     vsurf = 0.0
                     vmod(i,j,L) = interpVal(vsurf, vdef, pl_mod(i,j,L), &
                                                  plevs(lmpres), plevs(k2), MAPL_KAPPA)
                     !vmod(i,j,L) = interpVal(vdef, vsurf, pl_mod(i,j,L), &
                     !                        plevs(k2), plevs(lmpres),  MAPL_KAPPA)
                  ELSE IF (varType .EQ. 'T') THEN
                     ! Compute the potential temperature at pressure level k2
                     TH = vdef / ( (plevs(k2)*1.00E-05)**(MAPL_RGAS/MAPL_CP) )
                     ! Reverse back to the temperature at model level L
                     vmod(i,j,L) = TH * ( (pl_mod(i,j,L)*1.00E-05)**(MAPL_RGAS/MAPL_CP) )
                  ELSE IF (varType .EQ. 'H') THEN
                     vsurf = vPHIS(i,j)/gfactor
                     vmod(i,j,L) = interpVal(vsurf, vdef, pl_mod(i,j,L), &
                                                  plevs(lmpres), plevs(k2), MAPL_KAPPA)
                  ELSE
                     vmod(i,j,L) = vdef
                  END IF
               ENDDO
            ELSE
               DO L = kend(i,j)+1, lmmod
                  vmod(i,j,L) = vpres(i,j,lmpres)
               ENDDO
            END IF
         ENDDO
      ENDDO
      
      deallocate(ak)
      deallocate(bk)
      deallocate(pl_mod)
      deallocate(ple_mod)
      deallocate(kbeg)
      deallocate(kend)

      _RETURN(ESMF_SUCCESS)

      end subroutine vertInterpolation_pressKappa