Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Auto PR - develop → MAPL-v3 - pull in mapz implementation for T/PT #79

Merged
merged 4 commits into from
Sep 29, 2023
Merged

Auto PR - develop → MAPL-v3 - pull in mapz implementation for T/PT #79

Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
0bf3168
Merge pull request #72 from GEOS-ESM/geos/develop
mathomp4 Jun 8, 2023
28872fd
Merge branch 'geos/main' into geos/develop
mathomp4 Sep 28, 2023
4501bfb
pull in mapz implementation for T/PT
sdrabenh Sep 29, 2023
59ecf88
Merge pull request #78 from GEOS-ESM/feature/sdrabenh/hwt_spring_exp
mathomp4 Sep 29, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
152 changes: 76 additions & 76 deletions model/fv_mapz.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -212,7 +212,7 @@ subroutine Lagrangian_to_Eulerian(last_step, consv, ps, pe, delp, pkz, pk, &
real(kind=8), dimension(is:ie,js:je):: te_2d, zsum0, zsum1
real, dimension(is:ie,js:je):: dpln
real, dimension(is:ie,km) :: q2, dp2, w2
real, dimension(is:ie,km+1):: pe1, pe2, pk1, pk2, pn2, phis
real, dimension(is:ie,km+1):: pe1, pe2, pk1, pk2, pn1, pn2, phis
real, dimension(is:ie+1,km+1):: pe0, pe3
real, dimension(is:ie):: gz, cvm
real(kind=8):: tesum, zsum, dtmp
Expand Down Expand Up @@ -336,7 +336,7 @@ subroutine Lagrangian_to_Eulerian(last_step, consv, ps, pe, delp, pkz, pk, &
!$OMP ak,bk,nq,isd,ied,jsd,jed,kord_tr,fill, adiabatic, &
!$OMP hs,w,ws,kord_wz,rrg,kord_mt,consv,remap_option,gmao_remap) &
!$OMP private(gz,cvm,bkh,dp2, &
!$OMP pe0,pe1,pe2,pe3,pk1,pk2,pn2,phis,q2,w2,dpln,dlnp)
!$OMP pe0,pe1,pe2,pe3,pk1,pk2,pn1,pn2,phis,q2,w2,dpln,dlnp)
do 1000 j=js,je+1

do k=1,km+1
Expand Down Expand Up @@ -478,6 +478,7 @@ subroutine Lagrangian_to_Eulerian(last_step, consv, ps, pe, delp, pkz, pk, &
do k=1,km+1
do i=is,ie
pk1(i,k) = pk(i,j,k)
pn1(i,k) = peln(i,k,j)
enddo
enddo

Expand All @@ -495,44 +496,31 @@ subroutine Lagrangian_to_Eulerian(last_step, consv, ps, pe, delp, pkz, pk, &
enddo
enddo

if (remap_t) then
if (remap_te) then
!----------------------------------
! map T in log P
! map TE in log P
!----------------------------------
if ( gmao_remap > 0 ) then
call map1_gmao (km, pe1, pt, &
km, pe2, pt, &
is, ie, j, isd, ied, jsd, jed, akap, gmao_remap, T_VAR=3, conserv=.false.)
else
call map_scalar(km, peln(is,1,j), pt, gz, &
km, pn2, pt, &
is, ie, j, isd, ied, jsd, jed, 1, abs(kord_tm), t_min)
endif
elseif (remap_pt) then
!------------------------------------
! map PT in P^KAPPA
!------------------------------------
if ( gmao_remap > 0 ) then
call map1_gmao (km, pe1, pt, &
km, pe2, pt, &
is, ie, j, isd, ied, jsd, jed, akap, gmao_remap, T_VAR=2, conserv=.false.)
call map1_gmao (km, pe1, te, &
km, pe2, te, &
is, ie, j, isd, ied, jsd, jed, akap, gmao_remap, P_MAP=1, conserv=.true.)
else
call map1_ppm (km, pe1, pt, gz, &
km, pe2, pt, &
is, ie, j, isd, ied, jsd, jed, 1, abs(kord_tm))
call map_scalar(km, pn1, te, gz, &
km, pn2, te, &
is, ie, j, isd, ied, jsd, jed, 1, abs(kord_tm), te_min)
endif
elseif (remap_te) then
else
!----------------------------------
! map TE in log P
! map T or PT in log P
!----------------------------------
if ( gmao_remap > 0 ) then
call map1_gmao (km, pe1, te, &
km, pe2, te, &
is, ie, j, isd, ied, jsd, jed, akap, gmao_remap, T_VAR=1, conserv=.true.)
call map1_gmao (km, pe1, pt, &
km, pe2, pt, &
is, ie, j, isd, ied, jsd, jed, akap, gmao_remap, P_MAP=1, conserv=.false.)
else
call map_scalar(km, peln(is,1,j), te, gz, &
km, pn2, te, &
is, ie, j, isd, ied, jsd, jed, 1, abs(kord_tm), te_min)
call map_scalar(km, pn1, pt, gz, &
km, pn2, pt, &
is, ie, j, isd, ied, jsd, jed, 1, abs(kord_tm), t_min)
endif
endif

Expand Down Expand Up @@ -570,7 +558,7 @@ subroutine Lagrangian_to_Eulerian(last_step, consv, ps, pe, delp, pkz, pk, &
enddo
call map1_ppm (km, pe1, delz, gz, &
km, pe2, delz, &
is, ie, j, isd, ied, jsd, jed, 1, abs(kord_tm))
is, ie, j, isd, ied, jsd, jed, 1, abs(kord_wz))
do k=1,km
do i=is,ie
delz(i,j,k) = -delz(i,j,k)*dp2(i,k)
Expand Down Expand Up @@ -1476,6 +1464,7 @@ subroutine map_scalar( km, pe1, q1, qs, &
real q4(4,i1:i2,km)
real pl, pr, qsum, dp, esl
integer i, k, l, m, k0
integer LM1,LP0,LP1

do k=1,km
do i=i1,i2
Expand All @@ -1485,15 +1474,45 @@ subroutine map_scalar( km, pe1, q1, qs, &
enddo

! Compute vertical subgrid distribution
if ( kord >7 ) then
if ( kord > 7 ) then
call scalar_profile( qs, q4, dp1, km, i1, i2, iv, kord, q_min )
else
call ppm_profile( q4, dp1, km, i1, i2, iv, kord )
endif

do i=i1,i2
! Interpolate field onto target Pressures
! ---------------------------------------
do i=i1,i2
k0 = 1
do 555 k=1,kn
LM1 = 1
LP0 = 1
do while( LP0.le.km )
if (pe1(i,LP0).lt.pe2(i,k)) then
LP0 = LP0+1
else
exit
endif
enddo
LM1 = max(LP0-1,1)
LP0 = min(LP0, km)
! Extrapolate Linearly above first model level
! ----------------------------------------------------
if( LM1.eq.1 .and. LP0.eq.1 ) then
q2(i,j,k) = q1(i,j,1) + ( q1(i,j,2)-q1(i,j,1) )*( pe2(i,k)-pe1(i,1) ) &
/( pe1(i,2)-pe1(i,1) )
! Extrapolate Linearly below last model level
! ---------------------------------------------------
else if( LM1.eq.km .and. LP0.eq.km ) then
q2(i,j,k) = q1(i,j,km) + ( q1(i,j,km)-q1(i,j,km-1) )*( pe2(i,k )-pe1(i,km ) ) &
/( pe1(i,km)-pe1(i,km-1) )
! Interpolate Linearly between levels 1 => 2 and km-1 => km
! -----------------------------------------------------------------
else if( LM1.eq.1 .or. LP0.eq.km ) then
q2(i,j,k) = q1(i,j,LP0) + ( q1(i,j,LM1)-q1(i,j,LP0) )*( pe2(i,k )-pe1(i,LP0) ) &
/( pe1(i,LM1)-pe1(i,LP0) )
else

do l=k0,km
! locate the top edge: pe2(i,k)
if( pe2(i,k) >= pe1(i,l) .and. pe2(i,k) <= pe1(i,l+1) ) then
Expand Down Expand Up @@ -1529,6 +1548,8 @@ subroutine map_scalar( km, pe1, q1, qs, &
endif
enddo
123 q2(i,j,k) = qsum / ( pe2(i,k+1) - pe2(i,k) )

endif
555 continue
enddo

Expand Down Expand Up @@ -2005,6 +2026,7 @@ subroutine scalar_profile(qs, a4, delp, km, i1, i2, iv, kord, qmin)
return
endif
!----- Perfectly linear scheme --------------------------------

!------------------
! Apply constraints
!------------------
Expand Down Expand Up @@ -2285,7 +2307,7 @@ subroutine scalar_profile(qs, a4, delp, km, i1, i2, iv, kord, qmin)
do i=i1,i2
a4(4,i,k) = 3.*(2.*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k)))
enddo
else ! kord = 11, 13
elseif ( abs(kord)<16 ) then ! kord = 11, 13
do i=i1,i2
if ( ext5(i,k) .and. (ext5(i,k-1).or.ext5(i,k+1).or.a4(1,i,k)<qmin) ) then
! Noisy region:
Expand Down Expand Up @@ -3787,15 +3809,15 @@ end subroutine moist_cp
! !INTERFACE:
subroutine map1_gmao( km, pe1, q1, &
kn, pe2, q2, i1, i2, &
j, ibeg, iend, jbeg, jend, akap, gmao_remap, T_VAR, conserv)
j, ibeg, iend, jbeg, jend, akap, gmao_remap, P_MAP, conserv)
implicit none

! !INPUT PARAMETERS:
integer, intent(in) :: i1 ! Starting longitude
integer, intent(in) :: i2 ! Finishing longitude
real, intent(in) :: akap
integer, intent(in) :: T_VAR ! Thermodynamic variable to remap
! 1:TE 2:T 3:PT
integer, intent(in) :: P_MAP ! Thermodynamic variable to remap
! 1:TE-logP 2:PT-PK 3:T-logP
integer, intent(in) :: gmao_remap ! 3:cubic 2:quadratic 1:linear
logical, intent(in) :: conserv
integer, intent(in) :: j ! Current latitude
Expand Down Expand Up @@ -3844,66 +3866,47 @@ subroutine map1_gmao( km, pe1, q1, &
! Initialization
! --------------

select case (T_VAR)
select case (P_MAP)
case(1)
! Total Energy Remapping in Log(P)
! Remapping in Log(P)
do k=1,km
qx(:,k) = q1(i1:i2,j,k)
logpl1(:,k) = log( r2*(pe1(:,k)+pe1(:,k+1)) )
enddo
do k=1,kn
logpl2(:,k) = log( r2*(pe2(:,k)+pe2(:,k+1)) )
enddo

do k=1,km-1
dlogp1(:,k) = logpl1(:,k+1)-logpl1(:,k)
enddo

case(2)
! Temperature Remapping in Log(P)
do k=1,km
qx(:,k) = q1(i1:i2,j,k)
logpl1(:,k) = log( r2*(pe1(:,k)+pe1(:,k+1)) )
enddo
do k=1,kn
logpl2(:,k) = log( r2*(pe2(:,k)+pe2(:,k+1)) )
enddo

do k=1,km-1
dlogp1(:,k) = logpl1(:,k+1)-logpl1(:,k)
enddo

case(3)
! Potential Temperature Remapping in P^KAPPA
! Remapping in P^KAPPA
do k=1,km
qx(:,k) = q1(i1:i2,j,k)
logpl1(:,k) = exp( akap*log( r2*(pe1(:,k)+pe1(:,k+1))) )
enddo
do k=1,kn
logpl2(:,k) = exp( akap*log( r2*(pe2(:,k)+pe2(:,k+1))) )
enddo

do k=1,km-1
dlogp1(:,k) = logpl1(:,k+1)-logpl1(:,k)
enddo

end select

if (conserv) then
! Compute vertical integral of Input TE
! -------------------------------------
! Compute vertical integral of Input field
! ----------------------------------------
vsum1(:) = r0
do i=i1,i2
do k=1,km
vsum1(i) = vsum1(i) + qx(i,k)*( pe1(i,k+1)-pe1(i,k) )
enddo
vsum1(i) = vsum1(i) / ( pe1(i,km+1)-pe1(i,1) )
enddo

endif

! Interpolate TE onto target Pressures
! ------------------------------------
! Interpolate field onto target Pressures
! ---------------------------------------
do i=i1,i2
do k=1,kn
LM1 = 1
Expand All @@ -3918,24 +3921,24 @@ subroutine map1_gmao( km, pe1, q1, &
LM1 = max(LP0-1,1)
LP0 = min(LP0, km)

! Extrapolate Linearly in LogP above first model level
! Extrapolate Linearly above first model level
! ----------------------------------------------------
if( LM1.eq.1 .and. LP0.eq.1 ) then
q2(i,j,k) = qx(i,1) + ( qx(i,2)-qx(i,1) )*( logpl2(i,k)-logpl1(i,1) ) &
/( logpl1(i,2)-logpl1(i,1) )

! Extrapolate Linearly in LogP below last model level
! Extrapolate Linearly below last model level
! ---------------------------------------------------
else if( LM1.eq.km .and. LP0.eq.km ) then
q2(i,j,k) = qx(i,km) + ( qx(i,km)-qx(i,km-1) )*( logpl2(i,k )-logpl1(i,km ) ) &
/( logpl1(i,km)-logpl1(i,km-1) )

! Interpolate Linearly in LogP between levels 1 => 2 and km-1 => km
! Interpolate Linearly between levels 1 => 2 and km-1 => km
! -----------------------------------------------------------------
else if( LM1.eq.1 .or. LP0.eq.km ) then
q2(i,j,k) = qx(i,LP0) + ( qx(i,LM1)-qx(i,LP0) )*( logpl2(i,k )-logpl1(i,LP0) ) &
/( logpl1(i,LM1)-logpl1(i,LP0) )
! Interpolate Cubicly in LogP between other model levels
! Interpolate Cubicly between other model levels
! ------------------------------------------------------
else
LP1 = LP0+1
Expand Down Expand Up @@ -3979,27 +3982,24 @@ subroutine map1_gmao( km, pe1, q1, &

enddo
enddo
if (conserv) then

! Compute vertical integral of Output TE
! --------------------------------------
if (conserv) then
! Compute vertical integral of Output field
! -----------------------------------------
vsum2(:) = r0
do i=i1,i2
do k=1,kn
vsum2(i) = vsum2(i) + q2(i,j,k)*( pe2(i,k+1)-pe2(i,k) )
enddo
vsum2(i) = vsum2(i) / ( pe2(i,kn+1)-pe2(i,1) )
enddo

! Adjust Final TE to conserve
! ---------------------------
! Adjust Final field to conserve
! ------------------------------
do i=i1,i2
do k=1,kn
q2(i,j,k) = q2(i,j,k) + vsum1(i)-vsum2(i)
! q2(i,j,k) = q2(i,j,k) * vsum1(i)/vsum2(i)
enddo
enddo

endif

return
Expand Down
Loading