!------------------------------------------------------------------------------ ! Global Modeling and Assimilation Office (GMAO) ! ! Goddard Earth Observing System (GEOS) ! ! MAPL Component ! !------------------------------------------------------------------------------ ! #include "MAPL_ErrLog.h" ! !> !### MODULE: `linearVerticalInterpolation_mod` ! ! Author: GMAO SI-Team ! ! The module `linearVerticalInterpolation_mod` contains routines for ! linear vertical interpolation in the pressure to the KAPPA formulation. ! module linearVerticalInterpolation_mod use ESMF use MAPL_BaseMod use MAPL_Constants, only: MAPL_KAPPA, MAPL_RGAS, MAPL_CP, MAPL_GRAV use MAPL_ExceptionHandling use, intrinsic :: iso_fortran_env, only: REAL64 ! implicit none ! private ! public :: vertInterpolation_pressKappa ! !------------------------------------------------------------------------------ CONTAINS !------------------------------------------------------------------------------ !> ! 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. ! 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 type(ESMF_Info) :: infoh ! !------------------------------------------------------------------------------ ! get dimensions, allocate call ESMF_FieldGet(fModel,grid=grid,rc=status) _VERIFY(STATUS) call ESMF_InfoGetFromHost(fModel,infoh,rc=status) _VERIFY(STATUS) call ESMF_InfoGet(infoh,'UNITS',units,rc=status) _VERIFY(STATUS) call ESMF_InfoGet(infoh,'LONG_NAME',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_InfoGetFromHost(grid,infoh,rc=status) _VERIFY(STATUS) call ESMF_InfoGet(infoh,key='GridAK',values=ak,rc=status) _VERIFY(STATUS) call ESMF_InfoGet(infoh,key='GridBK',values=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 !EOC !------------------------------------------------------------------------------ !> ! Main vertical interpolation formula: ! !$$ !\alpha_{L} = \frac{ \alpha_{k-1}P_{L}^{\kappa} - \alpha_{k}P_{L}^{\kappa} + ! \alpha_{k}P_{k-1}^{\kappa} - \alpha_{k-1}P_{k}^{\kappa} } ! { P_{k-1}^{\kappa} - P_{k}^{\kappa} } !$$ ! ! In this particular case, we have: ! !$$ ! val = \frac{ var2 \times ple^{KAP} - var1 \times ple^{KAP} + ! var1 \times plev2^{KAP} - var2 \times plev1^{KAP} } ! { plev2^{KAP} - plevs^{KAP} } !$$ FUNCTION interpVal (var1, var2, ple, plev1, plev2, KAP) result(val) ! real :: var1 !! variable (pressure level) value at level k real :: var2 !! variable (pressure level) value at level k-1 real :: ple !! model pressure at current level real :: plev1 !! pressure value at level k real :: plev2 !! pressure value at level k-1 real :: KAP !! coefficient KAPPA real :: val !! Returned value ! !------------------------------------------------------------------------------ val = ( var1*(ple**KAP - plev2**KAP) - & var2*(ple**KAP - plev1**KAP) ) / (plev1**KAP - plev2**KAP) END FUNCTION interpVal ! !------------------------------------------------------------------------------ !> ! For each horizontal grid point, the subroutine `searchVerticalLocations` ! searches the last pressure level that is below lowest model levels and ! the first pressure level that is above the last model level. ! This allows to speed up the calculations while doing interpolations. ! subroutine searchVerticalLocations(pl_mod, plevs, kbeg, kend, & im, jm, lmmod, lmpres) ! integer, intent(IN) :: im !! number of longitude grid points integer, intent(IN) :: jm !! number of latitude grid points integer, intent(IN) :: lmmod !! number of model vertical levels integer, intent(IN) :: lmpres !! number of pressure vertical levels real, intent(IN) :: pl_mod(im,jm,lmmod) !! 3D mid pressure on model levels real, intent(IN) :: plevs(lmpres) !! pressure levels integer, intent(OUT) :: kbeg(im,jm) !! starting location integer, intent(OUT) :: kend(im,jm) !! ending location ! integer :: i, j, L, k ! !------------------------------------------------------------------------------ !----------------------------------------------------------------------- ! Find the starting location. ! The kbeg in each column is the index of model level data which has the ! smaller pressure (highest level) that is trapped by pressure levels. !----------------------------------------------------------------------- kbeg(:,:) = 1 DO j = 1, jm DO i = 1, im model_pressure_beg : DO L = 2, lmmod !pressure_search_beg : DO k = lmpres, 2, -1 pressure_search_beg : DO k = 2, lmpres IF ( (pl_mod(i,j,L) .LE. plevs(k )) .AND. & (pl_mod(i,j,L) .GT. plevs(k-1))) THEN kbeg(i,j) = L EXIT model_pressure_beg ELSE IF (pl_mod(i,j,L) .LT. plevs(k-1)) THEN EXIT pressure_search_beg END IF ENDDO pressure_search_beg ENDDO model_pressure_beg ENDDO ENDDO !------------------------------------------------------------------------ ! Find the ending location. ! The kend in each column is the index of model level data which has the ! largest pressure (lowest level) that is trapped by pressure levels. !------------------------------------------------------------------------ kend(:,:) = lmmod DO j = 1, jm DO i = 1, im model_pressure_end : DO L = lmmod, 2, -1 !pressure_search_end : DO k = 2, lmpres pressure_search_end : DO k = lmpres, 2, -1 IF ( (pl_mod(i,j,L) .LE. plevs(k )) .AND. & (pl_mod(i,j,L) .GT. plevs(k-1))) THEN kend(i,j) = L EXIT model_pressure_end !ELSE IF (pl_mod(i,j,L) .LT. plevs(k-1)) THEN ELSE IF (pl_mod(i,j,L) .GT. plevs(k)) THEN EXIT pressure_search_end END IF ENDDO pressure_search_end ENDDO model_pressure_end ENDDO ENDDO return end subroutine searchVerticalLocations ! !------------------------------------------------------------------------------ end module linearVerticalInterpolation_mod