Skip to content

Commit

Permalink
Merge pull request #879 from GEOS-ESM/feature/narnold1/shocmf_updates
Browse files Browse the repository at this point in the history 
SHOC-EDMF updates
  • Loading branch information
@sdrabenh
sdrabenh authored Jan 24, 2024
2 parents 5ce0644 + 971b1e7 commit d1830cd
Show file tree
Hide file tree
Showing 4 changed files with 276 additions and 172 deletions.
44 changes: 22 additions & 22 deletions GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/Process_Library.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1784,18 +1784,18 @@ subroutine partition_dblgss( dt, & ! IN

! Partition based on temperature for the first plume

IF (Tl1_1 >= tbgmax) THEN
esval1_1 = MAPL_EQsat(Tl1_1)
lstarn1 = lcond
ELSE IF (Tl1_1 < tbgmin) THEN
esval1_1 = MAPL_EQsat(Tl1_1,OverIce=.TRUE.)
lstarn1 = lsub
ELSE
! IF (Tl1_1 >= tbgmax) THEN
! esval1_1 = MAPL_EQsat(Tl1_1)
! lstarn1 = lcond
! ELSE IF (Tl1_1 < tbgmin) THEN
! esval1_1 = MAPL_EQsat(Tl1_1,OverIce=.TRUE.)
! lstarn1 = lsub
! ELSE
esval1_1 = MAPL_EQsat(Tl1_1)
esval2_1 = MAPL_EQsat(Tl1_1,OverIce=.TRUE.)
om1 = 1.-fQi !max(0.,min(1.,a_bg*(Tl1_1-tbgmin))) ! may be inconsistent with hystpdf ice fraction
om1 = max(0.,min(1.,1.-fQi)) !max(0.,min(1.,a_bg*(Tl1_1-tbgmin))) ! may be inconsistent with hystpdf ice fraction
lstarn1 = lcond + (1.-om1)*lfus
ENDIF
! ENDIF

! this is qs evaluated at Tl
qs1 = om1 * (0.622*esval1_1/max(esval1_1,pval-0.378*esval1_1)) &
Expand All @@ -1811,18 +1811,18 @@ subroutine partition_dblgss( dt, & ! IN
beta2 = beta1
ELSE

IF (Tl1_2 < tbgmin) THEN
esval1_2 = MAPL_EQsat(Tl1_2,OverIce=.TRUE.)
lstarn2 = lsub
ELSE IF (Tl1_2 >= tbgmax) THEN
esval1_2 = MAPL_EQsat(Tl1_2)
lstarn2 = lcond
ELSE
! IF (Tl1_2 < tbgmin) THEN
! esval1_2 = MAPL_EQsat(Tl1_2,OverIce=.TRUE.)
! lstarn2 = lsub
! ELSE IF (Tl1_2 >= tbgmax) THEN
! esval1_2 = MAPL_EQsat(Tl1_2)
! lstarn2 = lcond
! ELSE
esval1_2 = MAPL_EQsat(Tl1_2)
esval2_2 = MAPL_EQsat(Tl1_2,OverIce=.TRUE.)
om2 = 1.-fQi !max(0.,min(1.,a_bg*(Tl1_2-tbgmin)))
om2 = max(0.,min(1.,1.-fQi)) !max(0.,min(1.,a_bg*(Tl1_2-tbgmin)))
lstarn2 = lcond + (1.-om2)*lfus
ENDIF
! ENDIF

qs2 = om2 * (0.622*esval1_2/max(esval1_2,pval-0.378*esval1_2)) &
+ (1.-om2) * (0.622*esval2_2/max(esval2_2,pval-0.378*esval2_2))
Expand Down Expand Up @@ -1902,11 +1902,11 @@ subroutine partition_dblgss( dt, & ! IN
qn1 = min(qn1,qw1_1)
qn2 = min(qn2,qw1_2)

! ql1 = qn1*om1
! ql2 = qn2*om2
ql1 = qn1*om1
ql2 = qn2*om2

! qi1 = qn1 - ql1
! qi2 = qn2 - ql2
qi1 = qn1 - ql1
qi2 = qn2 - ql2

qc = min(max(0.0, aterm*qn1 + onema*qn2), total_water)
! diag_ql = min(max(0.0, aterm*ql1 + onema*ql2), diag_qn)
Expand Down
69 changes: 51 additions & 18 deletions GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSturbulence_GridComp/GEOS_TurbulenceGridComp.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -628,6 +628,24 @@ subroutine SetServices ( GC, RC )
! mass-flux export states
!

call MAPL_AddExportSpec(GC, &
LONG_NAME = 'EDMF_rain_tendency', &
UNITS = 'kg kg-1 s-1', &
SHORT_NAME = 'EDMF_DQRDT', &
DIMS = MAPL_DimsHorzVert, &
VLOCATION = MAPL_VLocationCenter, &
RC=STATUS )
VERIFY_(STATUS)

call MAPL_AddExportSpec(GC, &
LONG_NAME = 'EDMF_snow_tendency', &
UNITS = 'kg kg-1 s-1', &
SHORT_NAME = 'EDMF_DQSDT', &
DIMS = MAPL_DimsHorzVert, &
VLOCATION = MAPL_VLocationCenter, &
RC=STATUS )
VERIFY_(STATUS)

call MAPL_AddExportSpec(GC, &
LONG_NAME = 'Vertical_velocity_of_individual_EDMF_plumes', &
UNITS = 'm s-1', &
Expand Down Expand Up @@ -2829,7 +2847,7 @@ subroutine REFRESH(IM,JM,LM,RC)
edmf_w3, edmf_wqt, edmf_slqt, &
edmf_wsl, edmf_qt3, edmf_sl3, &
edmf_entx, edmf_tke, slflxmf, &
qtflxmf, mfaw
qtflxmf, mfaw, edmf_dqrdt, edmf_dqsdt

real, dimension(IM,JM,0:LM) :: ae3,aw3,aws3,awqv3,awql3,awqi3,awu3,awv3

Expand Down Expand Up @@ -3053,19 +3071,19 @@ subroutine REFRESH(IM,JM,LM,RC)
call MAPL_GetResource (MAPL, SHOCPARAMS%CEFAC, trim(COMP_NAME)//"_SHC_CEFAC:", default=1.0, RC=STATUS); VERIFY_(STATUS)
call MAPL_GetResource (MAPL, SHOCPARAMS%CESFAC, trim(COMP_NAME)//"_SHC_CESFAC:", default=4., RC=STATUS); VERIFY_(STATUS)
call MAPL_GetResource (MAPL, SHOCPARAMS%LENOPT, trim(COMP_NAME)//"_SHC_LENOPT:", default=3, RC=STATUS); VERIFY_(STATUS)
call MAPL_GetResource (MAPL, SHOCPARAMS%LENFAC1, trim(COMP_NAME)//"_SHC_LENFAC1:", default=4.0, RC=STATUS); VERIFY_(STATUS)
call MAPL_GetResource (MAPL, SHOCPARAMS%LENFAC2, trim(COMP_NAME)//"_SHC_LENFAC2:", default=1.0, RC=STATUS); VERIFY_(STATUS)
call MAPL_GetResource (MAPL, SHOCPARAMS%LENFAC3, trim(COMP_NAME)//"_SHC_LENFAC3:", default=2.0, RC=STATUS); VERIFY_(STATUS)
call MAPL_GetResource (MAPL, SHOCPARAMS%LENFAC1, trim(COMP_NAME)//"_SHC_LENFAC1:", default=10.0, RC=STATUS); VERIFY_(STATUS)
call MAPL_GetResource (MAPL, SHOCPARAMS%LENFAC2, trim(COMP_NAME)//"_SHC_LENFAC2:", default=2.0, RC=STATUS); VERIFY_(STATUS)
call MAPL_GetResource (MAPL, SHOCPARAMS%LENFAC3, trim(COMP_NAME)//"_SHC_LENFAC3:", default=1.5, RC=STATUS); VERIFY_(STATUS)
call MAPL_GetResource (MAPL, SHOCPARAMS%BUOYOPT, trim(COMP_NAME)//"_SHC_BUOY_OPTION:", default=2, RC=STATUS); VERIFY_(STATUS)
end if

call MAPL_GetResource (MAPL, PDFSHAPE, 'PDFSHAPE:', DEFAULT = 1.0 , RC=STATUS); VERIFY_(STATUS)
call MAPL_GetResource (MAPL, DOPROGQT2, 'DOPROGQT2:', DEFAULT = 1 , RC=STATUS); VERIFY_(STATUS)
call MAPL_GetResource (MAPL, SL2TUNE, 'SL2TUNE:', DEFAULT = 3.0 , RC=STATUS); VERIFY_(STATUS)
call MAPL_GetResource (MAPL, QT2TUNE, 'QT2TUNE:', DEFAULT = 6.0 , RC=STATUS); VERIFY_(STATUS)
call MAPL_GetResource (MAPL, SLQT2TUNE, 'SLQT2TUNE:', DEFAULT = 6.0 , RC=STATUS); VERIFY_(STATUS)
call MAPL_GetResource (MAPL, QT3_TSCALE, 'QT3_TSCALE:', DEFAULT = 2400.0, RC=STATUS); VERIFY_(STATUS)
call MAPL_GetResource (MAPL, AFRC_TSCALE,'AFRC_TSCALE:',DEFAULT = 1800.0, RC=STATUS); VERIFY_(STATUS)
call MAPL_GetResource (MAPL, SL2TUNE, 'SL2TUNE:', DEFAULT = 4.0 , RC=STATUS); VERIFY_(STATUS)
call MAPL_GetResource (MAPL, QT2TUNE, 'QT2TUNE:', DEFAULT = 7.0 , RC=STATUS); VERIFY_(STATUS)
call MAPL_GetResource (MAPL, SLQT2TUNE, 'SLQT2TUNE:', DEFAULT = 7.0 , RC=STATUS); VERIFY_(STATUS)
call MAPL_GetResource (MAPL, QT3_TSCALE, 'QT3_TSCALE:', DEFAULT = 1400.0, RC=STATUS); VERIFY_(STATUS)
call MAPL_GetResource (MAPL, AFRC_TSCALE,'AFRC_TSCALE:',DEFAULT = 1400.0, RC=STATUS); VERIFY_(STATUS)
call MAPL_GetResource (MAPL, DOCANUTO, 'DOCANUTO:', DEFAULT = 0, RC=STATUS); VERIFY_(STATUS)

! Get pointers from export state...
Expand Down Expand Up @@ -3183,6 +3201,16 @@ subroutine REFRESH(IM,JM,LM,RC)
VERIFY_(STATUS)
call MAPL_GetPointer(EXPORT, ZLES, 'ZLES', RC=STATUS)
VERIFY_(STATUS)
call MAPL_GetPointer(EXPORT, EDMF_PLUMES_W, 'EDMF_PLUMES_W', RC=STATUS)
VERIFY_(STATUS)
call MAPL_GetPointer(EXPORT, EDMF_PLUMES_QT, 'EDMF_PLUMES_QT', RC=STATUS)
VERIFY_(STATUS)
call MAPL_GetPointer(EXPORT, EDMF_PLUMES_THL, 'EDMF_PLUMES_THL', RC=STATUS)
VERIFY_(STATUS)
call MAPL_GetPointer(EXPORT, edmf_dqrdt, 'EDMF_DQRDT', ALLOC=.true., RC=STATUS)
VERIFY_(STATUS)
call MAPL_GetPointer(EXPORT, edmf_dqsdt, 'EDMF_DQSDT', ALLOC=.true., RC=STATUS)
VERIFY_(STATUS)
call MAPL_GetPointer(EXPORT, edmf_buoyf, 'EDMF_BUOYF', RC=STATUS)
VERIFY_(STATUS)
! call MAPL_GetPointer(EXPORT, edmf_sl2, 'EDMF_SL2', RC=STATUS)
Expand Down Expand Up @@ -3431,10 +3459,12 @@ subroutine REFRESH(IM,JM,LM,RC)
! number of updrafts
call MAPL_GetResource (MAPL, MFPARAMS%NUP, "EDMF_NUMUP:", default=10, RC=STATUS)
! boundaries for the updraft area (min/max sigma of w pdf)
call MAPL_GetResource (MAPL, MFPARAMS%PWMIN, "EDMF_PWMIN:", default=1.3, RC=STATUS)
call MAPL_GetResource (MAPL, MFPARAMS%PWMIN, "EDMF_PWMIN:", default=1., RC=STATUS)
call MAPL_GetResource (MAPL, MFPARAMS%PWMAX, "EDMF_PWMAX:", default=3., RC=STATUS)
!
call MAPL_GetResource (MAPL, MFPARAMS%ENTWFAC, "EDMF_ENTWFAC:", default=0.333, RC=STATUS)
call MAPL_GetResource (MAPL, MFPARAMS%ENTUFAC, "EDMF_ENTUFAC:", default=1.2, RC=STATUS)
call MAPL_GetResource (MAPL, MFPARAMS%WA, "EDMF_WA:", default=1.0, RC=STATUS)
call MAPL_GetResource (MAPL, MFPARAMS%WB, "EDMF_WB:", default=1.5, RC=STATUS)
! coefficients for surface forcing, appropriate for L137
call MAPL_GetResource (MAPL, MFPARAMS%AlphaW, "EDMF_ALPHAW:", default=0.05, RC=STATUS)
call MAPL_GetResource (MAPL, MFPARAMS%AlphaQT, "EDMF_ALPHAQT:", default=0.5, RC=STATUS)
Expand All @@ -3443,20 +3473,22 @@ subroutine REFRESH(IM,JM,LM,RC)
call MAPL_GetResource (MAPL, MFPARAMS%ET, "EDMF_ET:", default=2, RC=STATUS)
! constant entrainment rate
call MAPL_GetResource (MAPL, MFPARAMS%ENT0, "EDMF_ENT0:", default=0.25, RC=STATUS)
call MAPL_GetResource (MAPL, MFPARAMS%ENT0LTS, "EDMF_ENT0LTS:", default=0.75, RC=STATUS)
call MAPL_GetResource (MAPL, MFPARAMS%ENT0LTS, "EDMF_ENT0LTS:", default=1.5, RC=STATUS)
! L0 if ET==1
call MAPL_GetResource (MAPL, MFPARAMS%L0, "EDMF_L0:", default=100., RC=STATUS)
! L0fac if ET==2
call MAPL_GetResource (MAPL, MFPARAMS%L0fac, "EDMF_L0FAC:", default=10., RC=STATUS)
call MAPL_GetResource (MAPL, MFPARAMS%MFLIMFAC, "EDMF_MFLIMFAC:", default=3.5, RC=STATUS)
call MAPL_GetResource (MAPL, MFPARAMS%MFLIMFAC, "EDMF_MFLIMFAC:", default=2.5, RC=STATUS)
! factor to multiply the eddy-diffusivity with
call MAPL_GetResource (MAPL, MFPARAMS%EDfac, "EDMF_EDFAC:", default=1., RC=STATUS)
call MAPL_GetResource (MAPL, MFPARAMS%DOCLASP, "EDMF_DOCLASP:", default=0, RC=STATUS)
call MAPL_GetResource (MAPL, MFPARAMS%ICE_RAMP, "EDMF_ICE_RAMP:", default=-40.0, RC=STATUS )
call MAPL_GetResource (MAPL, MFPARAMS%ENTRAIN, "EDMF_ENTRAIN:", default=0, RC=STATUS)
call MAPL_GetResource (MAPL, MFPARAMS%STOCHFRAC, "EDMF_STOCHASTIC:", default=0.5, RC=STATUS)
call MAPL_GetResource (MAPL, MFPARAMS%DISCRETE, "EDMF_DISCRETE_TYPE:", default=0, RC=STATUS)
call MAPL_GetResource (MAPL, MFPARAMS%DISCRETE, "EDMF_DISCRETE_TYPE:", default=1, RC=STATUS)
call MAPL_GetResource (MAPL, MFPARAMS%IMPLICIT, "EDMF_IMPLICIT:", default=1, RC=STATUS)
call MAPL_GetResource (MAPL, MFPARAMS%PRCPCRIT, "EDMF_PRCPCRIT:", default=-1., RC=STATUS)
call MAPL_GetResource (MAPL, MFPARAMS%UPABUOYDEP,"EDMF_UPABUOYDEP:", default=1, RC=STATUS)

! Future options
! call MAPL_GetResource (MAPL, EDMF_THERMAL_PLUME, "EDMF_THERMAL_PLUME:", default=0, RC=STATUS)
Expand Down Expand Up @@ -3529,7 +3561,7 @@ subroutine REFRESH(IM,JM,LM,RC)
! bstar_scm = (MAPL_GRAV/(RHOS*sqrt(CM*max(UU,1.e-30)/RHOS))) * &
! (CT*(TH-TA-(MAPL_GRAV/MAPL_CP)*DZ)/TA + MAPL_VIREPS*CQ*(QH-QA))

ustar => ustar_scm
! ustar => ustar_scm
sh => sh_scm
evap => evap_scm

Expand Down Expand Up @@ -3604,6 +3636,7 @@ subroutine REFRESH(IM,JM,LM,RC)
buoyf, &
edmf_mf, & ! needed for ADG PDF
edmfdrya, edmfmoista, & ! outputs for ADG PDF
edmf_dqrdt, edmf_dqsdt, & ! output for micro
!== Diagnostics, not used elsewhere ==
edmf_dry_w, &
edmf_moist_w, &
Expand Down Expand Up @@ -5991,7 +6024,7 @@ subroutine UPDATE(IM,JM,LM,LATS,RC)
tmp3d(:,:,1:LM-1) = -1.*TKH(:,:,1:LM-1)*tmp3d(:,:,1:LM-1)
tmp3d(:,:,LM) = tmp3d(:,:,LM-1)
tmp3d(:,:,0) = 0.0
if (associated(QTFLXMF)) then
if (associated(QTFLXMF).and.MFPARAMS%IMPLICIT.eq.1) then
QTFLXMF(:,:,1:LM-1) = QTFLXMF(:,:,1:LM-1)-MFAW(:,:,1:LM-1)*QT(:,:,1:LM-1)
QTFLXMF(:,:,LM) = QTFLXMF(:,:,LM-1)
QTFLXMF(:,:,0) = 0.
Expand All @@ -6004,9 +6037,9 @@ subroutine UPDATE(IM,JM,LM,LATS,RC)
tmp3d(:,:,1:LM-1) = -1.*TKH(:,:,1:LM-1)*tmp3d(:,:,1:LM-1)
tmp3d(:,:,LM) = tmp3d(:,:,LM-1)
tmp3d(:,:,0) = 0.0
if (associated(SLFLXMF)) then
if (associated(SLFLXMF).and.MFPARAMS%IMPLICIT.eq.1) then
SLFLXMF(:,:,1:LM-1) = SLFLXMF(:,:,1:LM-1)-MFAW(:,:,1:LM-1)*SL(:,:,1:LM-1)/MAPL_CP
SLFLXMF(:,:,LM) = 0.
SLFLXMF(:,:,LM) = SLFLXMF(:,:,LM-1)
SLFLXMF(:,:,0) = 0.
end if
if (associated(SLFLXTRB)) SLFLXTRB = tmp3d/MAPL_CP + SLFLXMF
Expand Down
Loading

0 comments on commit d1830cd

Please sign in to comment.