From 6b559960f651c57fc6d55e08af07ff1fe25dd5d2 Mon Sep 17 00:00:00 2001 From: Nathan Arnold Date: Tue, 30 Jan 2024 22:30:33 -0500 Subject: [PATCH 1/8] Intermediate commit. New SHOC L2 based on parcel, adjusted parameters, adjusted mfdepth calc. --- .../GEOSmoist_GridComp/GEOS_MoistGridComp.F90 | 2 +- .../GEOSmoist_GridComp/Process_Library.F90 | 29 +++- .../GEOSmoist_GridComp/cloudnew.F90 | 3 +- .../GEOS_TurbulenceGridComp.F90 | 142 +++++++++++++++--- .../GEOSturbulence_GridComp/edmf.F90 | 43 +++++- .../GEOSturbulence_GridComp/shoc.F90 | 100 +++++++----- 6 files changed, 248 insertions(+), 71 deletions(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/GEOS_MoistGridComp.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/GEOS_MoistGridComp.F90 index 790cf4174..8d5dc7195 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/GEOS_MoistGridComp.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/GEOS_MoistGridComp.F90 @@ -666,7 +666,7 @@ subroutine SetServices ( GC, RC ) call MAPL_AddImportSpec(GC, & SHORT_NAME = 'RADSW', & - LONG_NAME = 'air_temperature_tendency_due_to_longwave', & + LONG_NAME = 'air_temperature_tendency_due_to_shortwave', & UNITS = 'K s-1', & DIMS = MAPL_DimsHorzVert, & VLOCATION = MAPL_VLocationCenter, & diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/Process_Library.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/Process_Library.F90 index 363c2bd9d..05d575644 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/Process_Library.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/Process_Library.F90 @@ -1,6 +1,7 @@ ! $Id$ #include "MAPL_Generic.h" +!#define PDFDIAG 1 !============================================================================= !BOP @@ -1290,8 +1291,7 @@ end subroutine pdfcondensate ! References in the comments in this code are given to ! ! the Appendix A of Pete Bogenschutz's dissertation. ! !------------------------------------------------------------------------! - subroutine partition_dblgss( dt, & ! IN - fQi, & + subroutine partition_dblgss( fQi, & ! IN tabs, & qwv, & qc, & ! OUT @@ -1346,7 +1346,7 @@ subroutine partition_dblgss( dt, & ! IN use MAPL_SatVaporMod, only: MAPL_EQsat real, intent(in ) :: fQi ! ice fraction - real, intent(in ) :: DT ! timestep [s] +! real, intent(in ) :: DT ! timestep [s] real, intent(in ) :: tabs ! absolute temperature [K] real, intent(in ) :: qwv ! specific humidity [kg kg-1] real, intent( out) :: qc ! liquid+ice condensate [kg kg-1] @@ -1905,13 +1905,31 @@ subroutine partition_dblgss( dt, & ! IN rwthl = 0.0 end if + if (abs(rwqt).gt.1e3) print *,'rwqt>1e3' + if (abs(rwthl).gt.1e3) print *,'rwthl>1e3' + if (abs(C1).gt.1e3) print *,'C1>1e3' + if (abs(C2).gt.1e3) print *,'C2>1e3' + if (abs(w1_1).gt.1e3) print *,'w1_1>1e3' + if (abs(w1_2).gt.1e3) print *,'w1_2>1e3' + if (abs(qw2_2).gt.1e3) print *,'qw2_2>1e3' + if (abs(qw2_1).gt.1e3) print *,'qw2_1>1e3' + if (abs(cqt1).gt.1e3) print *,'cqt1>1e3' + if (abs(cqt2).gt.1e3) print *,'cqt2>1e3' + if (abs(cthl1).gt.1e3) print *,'cthl1>1e3' + if (abs(cthl2).gt.1e3) print *,'cthl2>1e3' + if (abs(thl2_1).gt.1e3) print *,'thl2_1>1e3' + if (abs(thl2_2).gt.1e3) print *,'thl2_2>1e3' + wql1 = C1*(cqt1*sqrt(w2_1)*sqrt(qw2_1)*rwqt-cthl1*sqrt(w2_1)*sqrt(thl2_1)*rwthl) wql2 = C2*(cqt2*sqrt(w2_2)*sqrt(qw2_2)*rwqt-cthl2*sqrt(w2_2)*sqrt(thl2_2)*rwthl) - + if (abs(wql1).gt.1e3) print *,'wql1>1e3' + if (abs(wql2).gt.1e3) print *,'wql2>1e3' ! Compute the liquid water flux wqls = aterm * ((w1_1-w_first)*ql1+wql1) + onema * ((w1_2-w_first)*ql2+wql2) wqis = aterm * ((w1_1-w_first)*qi1) + onema * ((w1_2-w_first)*qi2) + if (isnan(wqls)) print *,'wqls is nan at p=',pval + if (wqls>1e3) print *,'wqls >1e3 at p=',pval ! diagnostic buoyancy flux. Includes effects from liquid water, ice ! condensate, liquid & ice precipitation @@ -2072,8 +2090,7 @@ subroutine hystpdf( & alhxbcp = (1.0-fQi)*alhlbcp + fQi*alhsbcp HL = TEn + gravbcp*ZL - alhxbcp*QCn - call partition_dblgss(DT/nmax, & - fQi, & + call partition_dblgss(fQi, & TEn, & QVn, & QCn, & diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/cloudnew.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/cloudnew.F90 index 7b684e1fc..3f2006297 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/cloudnew.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/cloudnew.F90 @@ -1992,8 +1992,7 @@ subroutine hystpdf_new( & ALHX = (1.0-fQi)*MAPL_ALHL + fQi*MAPL_ALHS HL = TEn + (mapl_grav/mapl_cp)*ZL - (ALHX/MAPL_CP)*QCn - call partition_dblgss(DT/nmax, & - fQi, & + call partition_dblgss(fQi, & TEn, & QVn, & QCn, & diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSturbulence_GridComp/GEOS_TurbulenceGridComp.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSturbulence_GridComp/GEOS_TurbulenceGridComp.F90 index 02a2be257..72817571d 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSturbulence_GridComp/GEOS_TurbulenceGridComp.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSturbulence_GridComp/GEOS_TurbulenceGridComp.F90 @@ -1020,6 +1020,16 @@ subroutine SetServices ( GC, RC ) RC=STATUS ) VERIFY_(STATUS) + + call MAPL_AddExportSpec(GC, & + LONG_NAME = 'EDMF_plume_depth_for_entrainment', & + UNITS = 'm', & + SHORT_NAME = 'EDMF_DEPTH2', & + DIMS = MAPL_DimsHorzOnly, & + VLOCATION = MAPL_VLocationNone, & + RC=STATUS ) + VERIFY_(STATUS) + call MAPL_AddExportSpec(GC, & LONG_NAME = 'EDMF_mass_flux', & UNITS = 'kg m s-1', & @@ -1029,6 +1039,33 @@ subroutine SetServices ( GC, RC ) RC=STATUS ) VERIFY_(STATUS) + call MAPL_AddExportSpec(GC, & + LONG_NAME = 'EDMF_dry_static_energy_source_term', & + UNITS = 'J kg-1 s-1', & + SHORT_NAME = 'SSRCMF', & + DIMS = MAPL_DimsHorzVert, & + VLOCATION = MAPL_VLocationCenter, & + RC=STATUS ) + VERIFY_(STATUS) + + call MAPL_AddExportSpec(GC, & + LONG_NAME = 'EDMF_specific_humidity_source_term', & + UNITS = 'kg kg-1 s-1', & + SHORT_NAME = 'QVSRCMF', & + DIMS = MAPL_DimsHorzVert, & + VLOCATION = MAPL_VLocationCenter, & + RC=STATUS ) + VERIFY_(STATUS) + + call MAPL_AddExportSpec(GC, & + LONG_NAME = 'EDMF_liquid_water_source_term', & + UNITS = 'kg kg-1 s-1', & + SHORT_NAME = 'QLSRCMF', & + DIMS = MAPL_DimsHorzVert, & + VLOCATION = MAPL_VLocationCenter, & + RC=STATUS ) + VERIFY_(STATUS) + call MAPL_AddExportSpec(GC, & SHORT_NAME = 'SLFLXMF', & LONG_NAME = 'liquid_water_static_energy_flux_by_MF', & @@ -1872,6 +1909,22 @@ subroutine SetServices ( GC, RC ) VLOCATION = MAPL_VLocationCenter, RC=STATUS ) VERIFY_(STATUS) + call MAPL_AddExportSpec(GC, & + SHORT_NAME = 'BRUNTDRY', & + LONG_NAME = 'Brunt_Vaisala_frequency_from_SHOC', & + UNITS = 's-1', & + DIMS = MAPL_DimsHorzVert, & + VLOCATION = MAPL_VLocationCenter, RC=STATUS ) + VERIFY_(STATUS) + + call MAPL_AddExportSpec(GC, & + SHORT_NAME = 'BRUNTEDGE', & + LONG_NAME = 'Brunt_Vaisala_frequency_from_SHOC', & + UNITS = 's-1', & + DIMS = MAPL_DimsHorzVert, & + VLOCATION = MAPL_VLocationEdge, RC=STATUS ) + VERIFY_(STATUS) + call MAPL_AddExportSpec(GC, & LONG_NAME = 'edge_height_above_surface', & SHORT_NAME = 'ZLES', & @@ -2829,12 +2882,12 @@ subroutine REFRESH(IM,JM,LM,RC) real, dimension(:,:,:), pointer :: AKQODT, CKQODT real, dimension(:,:,:), pointer :: AKVODT, CKVODT - real, dimension(:,:,:), pointer :: LSHOC,BRUNTSHOC,ISOTROPY, & + real, dimension(:,:,:), pointer :: LSHOC,BRUNTSHOC,BRUNTDRY, BRUNTEDGE,ISOTROPY, & LSHOC1,LSHOC2,LSHOC3, & SHOCPRNUM,& TKEBUOY,TKESHEAR,TKEDISS,TKETRANS, & SL2, SL3, W2, W3, WQT, WSL, SLQT, W3CANUTO, QT2DIAG,SL2DIAG,SLQTDIAG - real, dimension(:,:), pointer :: LMIX, edmf_depth + real, dimension(:,:), pointer :: LMIX, edmf_depth, edmf_depth2 ! EDMF variables real, dimension(:,:,:), pointer :: edmf_dry_a,edmf_moist_a,edmf_frc, edmf_dry_w,edmf_moist_w, & @@ -2847,9 +2900,11 @@ 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, edmf_dqrdt, edmf_dqsdt + qtflxmf, mfaw, edmf_dqrdt, edmf_dqsdt, & + ssrcmf,qvsrcmf,qlsrcmf real, dimension(IM,JM,0:LM) :: ae3,aw3,aws3,awqv3,awql3,awqi3,awu3,awv3 + real, dimension(IM,JM,1:LM) :: ssrc,qvsrc,qlsrc real, dimension(IM,JM) :: zpbl_test @@ -3065,7 +3120,7 @@ subroutine REFRESH(IM,JM,LM,RC) call MAPL_GetResource (MAPL, DO_SHOC, trim(COMP_NAME)//"_DO_SHOC:", default=0, RC=STATUS); VERIFY_(STATUS) if (DO_SHOC /= 0) then call MAPL_GetResource (MAPL, SHOCPARAMS%PRNUM, trim(COMP_NAME)//"_SHC_PRNUM:", default=-1.0, RC=STATUS); VERIFY_(STATUS) - call MAPL_GetResource (MAPL, SHOCPARAMS%LAMBDA, trim(COMP_NAME)//"_SHC_LAMBDA:", default=0.04, RC=STATUS); VERIFY_(STATUS) + call MAPL_GetResource (MAPL, SHOCPARAMS%LAMBDA, trim(COMP_NAME)//"_SHC_LAMBDA:", default=0.08, RC=STATUS); VERIFY_(STATUS) call MAPL_GetResource (MAPL, SHOCPARAMS%TSCALE, trim(COMP_NAME)//"_SHC_TSCALE:", default=400., RC=STATUS); VERIFY_(STATUS) call MAPL_GetResource (MAPL, SHOCPARAMS%CKVAL, trim(COMP_NAME)//"_SHC_CK:", default=0.1, RC=STATUS); VERIFY_(STATUS) call MAPL_GetResource (MAPL, SHOCPARAMS%CEFAC, trim(COMP_NAME)//"_SHC_CEFAC:", default=1.0, RC=STATUS); VERIFY_(STATUS) @@ -3073,7 +3128,7 @@ subroutine REFRESH(IM,JM,LM,RC) 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=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%LENFAC3, trim(COMP_NAME)//"_SHC_LENFAC3:", default=3.0, RC=STATUS); VERIFY_(STATUS) call MAPL_GetResource (MAPL, SHOCPARAMS%BUOYOPT, trim(COMP_NAME)//"_SHC_BUOY_OPTION:", default=2, RC=STATUS); VERIFY_(STATUS) end if @@ -3257,6 +3312,12 @@ subroutine REFRESH(IM,JM,LM,RC) VERIFY_(STATUS) call MAPL_GetPointer(EXPORT, edmf_mfx, 'EDMF_MF', RC=STATUS) VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT, ssrcmf, 'SSRCMF', RC=STATUS) + VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT, qvsrcmf, 'QVSRCMF', RC=STATUS) + VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT, qlsrcmf, 'QLSRCMF', RC=STATUS) + VERIFY_(STATUS) call MAPL_GetPointer(EXPORT, edmf_dry_a, 'EDMF_DRY_A', RC=STATUS) VERIFY_(STATUS) call MAPL_GetPointer(EXPORT, edmf_moist_a, 'EDMF_MOIST_A', RC=STATUS) @@ -3289,6 +3350,8 @@ subroutine REFRESH(IM,JM,LM,RC) VERIFY_(STATUS) call MAPL_GetPointer(EXPORT, edmf_depth, 'EDMF_DEPTH', RC=STATUS) VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT, edmf_depth2, 'EDMF_DEPTH2', RC=STATUS) + VERIFY_(STATUS) call MAPL_GetPointer(EXPORT, mfaw, 'MFAW', RC=STATUS) VERIFY_(STATUS) call MAPL_GetPointer(EXPORT, slflxmf, 'SLFLXMF', RC=STATUS) @@ -3319,6 +3382,10 @@ subroutine REFRESH(IM,JM,LM,RC) VERIFY_(STATUS) call MAPL_GetPointer(EXPORT, BRUNTSHOC, 'BRUNTSHOC', ALLOC=PDFALLOC, RC=STATUS) VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT, BRUNTDRY, 'BRUNTDRY', RC=STATUS) + VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT, BRUNTEDGE, 'BRUNTEDGE', RC=STATUS) + VERIFY_(STATUS) call MAPL_GetPointer(EXPORT, SHOCPRNUM,'SHOCPRNUM', RC=STATUS) VERIFY_(STATUS) @@ -3446,6 +3513,33 @@ subroutine REFRESH(IM,JM,LM,RC) call MAPL_TimerOff(MAPL,"---PRELIMS") + ! find thv-based inversion depth + if (associated(EDMF_DEPTH2)) then + tmp3d(:,:,1:LM-1) = (thv(:,:,2:LM)-thv(:,:,1:LM-1)) / (Z(:,:,2:LM)-Z(:,:,1:LM-1)) + do J = 1,JM + do I = 1,IM + temparray = 0. + K = LM + DO WHILE(Z(I,J,K).lt.3e3) + K = K-1 + END DO +! print *,'3km hgt=',Z(I,J,K) + ! at this point, K is index of full level above 3km height + locmax = maxloc(tmp3d(I,J,K:),dim=1) ! returns edge index of max thv gradient + K = K-1+locmax +! print *,'thv max hgt=',Z(I,J,K) + temparray(K:LM-1) = (tmp3d(I,J,K-1:LM-2)-tmp3d(I,J,K:LM-1)) / (Z(I,J,K-1:LM-2)-Z(I,J,K:LM-1)) ! 2nd deriv of THV on full levs +! print *,'2nd deriv=',temparray(K:LM-1) + DO WHILE( K.lt.LM ) + K = K+1 + if (temparray(K).gt.temparray(K+1)) exit + END DO +! print *,'inversion hgt=',Z(I,J,K) + + EDMF_DEPTH2(I,J) = Z(I,J,K) + end do + end do + end if ! Calculate liquid water potential temperature (THL) and total water (QT) EXF=T/TH @@ -3462,17 +3556,17 @@ subroutine REFRESH(IM,JM,LM,RC) 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%ENTUFAC, "EDMF_ENTUFAC:", default=1.2, RC=STATUS) + call MAPL_GetResource (MAPL, MFPARAMS%ENTUFAC, "EDMF_ENTUFAC:", default=1.6, 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) - call MAPL_GetResource (MAPL, MFPARAMS%AlphaTH, "EDMF_ALPHATH:", default=2.89, RC=STATUS) + call MAPL_GetResource (MAPL, MFPARAMS%AlphaQT, "EDMF_ALPHAQT:", default=1.0, RC=STATUS) + call MAPL_GetResource (MAPL, MFPARAMS%AlphaTH, "EDMF_ALPHATH:", default=1.0, RC=STATUS) ! Entrainment rate options 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%ENT0, "EDMF_ENT0:", default=0.3, 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) @@ -3557,11 +3651,11 @@ subroutine REFRESH(IM,JM,LM,RC) cu => cu_scm ct => ct_scm cq => ct_scm -! ustar_scm = 0.3! sqrt(CM*UU/RHOS) + ustar_scm = 0.3 ! sqrt(CU*UU/RHOS) ! 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 @@ -3587,6 +3681,9 @@ subroutine REFRESH(IM,JM,LM,RC) mfwqt = 0.0 mfwsl = 0.0 mftke = 0.0 + ssrc = 0.0 + qvsrc = 0.0 + qlsrc = 0.0 IF(DOMF /= 0) then @@ -3621,6 +3718,9 @@ subroutine REFRESH(IM,JM,LM,RC) awqi3, & awu3, & awv3, & + ssrc, & + qvsrc, & + qlsrc, & !== Outputs for ADG PDF == mfw2, & mfw3, & @@ -3662,6 +3762,9 @@ subroutine REFRESH(IM,JM,LM,RC) if (associated(mfaw)) mfaw = edmf_mf/rhoe if (associated(slflxmf)) slflxmf = (aws3-awql3*mapl_alhl-awqi3*mapl_alhs)/mapl_cp if (associated(qtflxmf)) qtflxmf = awqv3+awql3+awqi3 + if (associated(ssrcmf)) ssrcmf = ssrc + if (associated(qvsrcmf)) qvsrcmf = qvsrc + if (associated(qlsrcmf)) qlsrcmf = qlsrc ! if (associated(edmf_sl2)) edmf_sl2 = mfsl2 ! if (associated(edmf_qt2)) edmf_qt2 = mfqt2 if (associated(edmf_w2)) edmf_w2 = mfw2 @@ -3703,6 +3806,9 @@ subroutine REFRESH(IM,JM,LM,RC) if (associated(edmf_entx)) edmf_entx = MAPL_UNDEF if (associated(edmf_mfx)) edmf_mfx = 0.0 if (associated(mfaw)) mfaw = 0.0 + if (associated(ssrcmf)) ssrcmf = 0.0 + if (associated(qlsrcmf)) qlsrcmf = 0.0 + if (associated(qvsrcmf)) qvsrcmf = 0.0 if (associated(slflxmf)) slflxmf = 0.0 if (associated(qtflxmf)) qtflxmf = 0.0 ! if (associated(edmf_sl2)) edmf_sl2 = mfsl2 @@ -3774,6 +3880,8 @@ subroutine REFRESH(IM,JM,LM,RC) LSHOC2, & LSHOC3, & BRUNTSHOC, & + BRUNTDRY, & + BRUNTEDGE, & RI, & SHOCPRNUM, & !== Tuning params == @@ -4728,16 +4836,16 @@ subroutine REFRESH(IM,JM,LM,RC) ! ! 2:LM -> 1:LM-1, 1:LM-1 -> 0:LM-2 ! - YS(:,:,LM) = -DMI(:,:,LM)*RHOE(:,:,LM-1)*AWS3(:,:,LM-1) - YQV(:,:,LM) = -DMI(:,:,LM)*RHOE(:,:,LM-1)*AWQV3(:,:,LM-1) - YQL(:,:,LM) = -DMI(:,:,LM)*RHOE(:,:,LM-1)*AWQL3(:,:,LM-1) + YS(:,:,LM) = -DMI(:,:,LM)*( RHOE(:,:,LM-1)*AWS3(:,:,LM-1) + SSRC(:,:,LM) ) + YQV(:,:,LM) = -DMI(:,:,LM)*( RHOE(:,:,LM-1)*AWQV3(:,:,LM-1) + QVSRC(:,:,LM) ) + YQL(:,:,LM) = -DMI(:,:,LM)*( RHOE(:,:,LM-1)*AWQL3(:,:,LM-1) + QLSRC(:,:,LM) ) YQI(:,:,LM) = -DMI(:,:,LM)*RHOE(:,:,LM-1)*AWQI3(:,:,LM-1) YU(:,:,LM) = -DMI(:,:,LM)*RHOE(:,:,LM-1)*AWU3(:,:,LM-1) YV(:,:,LM) = -DMI(:,:,LM)*RHOE(:,:,LM-1)*AWV3(:,:,LM-1) - YS(:,:,1:LM-1) = DMI(:,:,1:LM-1)*( RHOE(:,:,1:LM-1)*AWS3(:,:,1:LM-1) - RHOE(:,:,0:LM-2)*AWS3(:,:,0:LM-2) ) - YQV(:,:,1:LM-1) = DMI(:,:,1:LM-1)*( RHOE(:,:,1:LM-1)*AWQV3(:,:,1:LM-1) - RHOE(:,:,0:LM-2)*AWQV3(:,:,0:LM-2) ) - YQL(:,:,1:LM-1) = DMI(:,:,1:LM-1)*( RHOE(:,:,1:LM-1)*AWQL3(:,:,1:LM-1) - RHOE(:,:,0:LM-2)*AWQL3(:,:,0:LM-2) ) + YS(:,:,1:LM-1) = DMI(:,:,1:LM-1)*( RHOE(:,:,1:LM-1)*AWS3(:,:,1:LM-1) - RHOE(:,:,0:LM-2)*AWS3(:,:,0:LM-2) + SSRC(:,:,1:LM-1) ) + YQV(:,:,1:LM-1) = DMI(:,:,1:LM-1)*( RHOE(:,:,1:LM-1)*AWQV3(:,:,1:LM-1) - RHOE(:,:,0:LM-2)*AWQV3(:,:,0:LM-2) + QVSRC(:,:,1:LM-1) ) + YQL(:,:,1:LM-1) = DMI(:,:,1:LM-1)*( RHOE(:,:,1:LM-1)*AWQL3(:,:,1:LM-1) - RHOE(:,:,0:LM-2)*AWQL3(:,:,0:LM-2) + QLSRC(:,:,1:LM-1) ) YQI(:,:,1:LM-1) = DMI(:,:,1:LM-1)*( RHOE(:,:,1:LM-1)*AWQI3(:,:,1:LM-1) - RHOE(:,:,0:LM-2)*AWQI3(:,:,0:LM-2) ) YU(:,:,1:LM-1) = DMI(:,:,1:LM-1)*( RHOE(:,:,1:LM-1)*AWU3(:,:,1:LM-1) - RHOE(:,:,0:LM-2)*AWU3(:,:,0:LM-2) ) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSturbulence_GridComp/edmf.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSturbulence_GridComp/edmf.F90 index 5063baa04..e1d96d8e1 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSturbulence_GridComp/edmf.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSturbulence_GridComp/edmf.F90 @@ -93,6 +93,9 @@ SUBROUTINE RUN_EDMF(its,ite, jts,jte, kts,kte, dt, & ! Inputs awqi3, & awu3, & awv3, & + ssrc3, & + qvsrc3, & + qlsrc3, & ! Outputs required for SHOC and ADG PDF mfw2, & mfw3, & @@ -170,7 +173,7 @@ SUBROUTINE RUN_EDMF(its,ite, jts,jte, kts,kte, dt, & ! Inputs mftke REAL,DIMENSION(ITS:ITE,JTS:JTE,KTS:KTE), INTENT(OUT) :: buoyf,mfw2,mfw3,mfqt3,mfhl3,&!mfqt2,mfhl2,& - mfhlqt,dqrdt,dqsdt + mfhlqt,dqrdt,dqsdt,ssrc3,qvsrc3,qlsrc3 ! Diagnostic outputs @@ -278,6 +281,9 @@ SUBROUTINE RUN_EDMF(its,ite, jts,jte, kts,kte, dt, & ! Inputs awqi3 =0. awu3 =0. awv3 =0. + ssrc3 =0. + qvsrc3=0. + qlsrc3=0. buoyf =0. mfw2 =0. mfw3 =0. @@ -615,6 +621,7 @@ SUBROUTINE RUN_EDMF(its,ite, jts,jte, kts,kte, dt, & ! Inputs Wn2=EntW*UPW(k-1,I)**2+MFPARAMS%WA*B/WP*(1.-EntW) END IF + IF (Wn2>0.) THEN UPW(K,I)=sqrt(Wn2) UPTHV(K,I)=THVn @@ -625,6 +632,13 @@ SUBROUTINE RUN_EDMF(its,ite, jts,jte, kts,kte, dt, & ! Inputs UPU(K,I)=Un UPV(K,I)=Vn UPA(K,I)=UPA(K-1,I) + + ! trisolver source terms due to condensation. assumes that condensation is responsible + ! for any increase in condesate flux with height. ignores lateral mixing. + tmp = max(0.,UPA(K,I)*(RHOE(K)*UPW(K,I)*UPQL(K,I)-RHOE(K-1)*UPW(K-1,I)*UPQL(K-1,I))) + qlsrc3(IH,JH,KTE+KTS-K) = qlsrc3(IH,JH,KTE+KTS-K) + tmp + qvsrc3(IH,JH,KTE+KTS-K) = qvsrc3(IH,JH,KTE+KTS-K) - tmp + ssrc3(IH,JH,KTE+KTS-K) = ssrc3(IH,JH,KTE+KTS-K) + tmp*MAPL_ALHL ELSE UPW(K,I) = 0. UPA(K,I) = 0. @@ -706,6 +720,21 @@ SUBROUTINE RUN_EDMF(its,ite, jts,jte, kts,kte, dt, & ! Inputs ENDDO DQRDT(IH,JH,KTS:KTE) = QR(KTE:KTS:-1)/DT DQSDT(IH,JH,KTS:KTE) = QS(KTE:KTS:-1)/DT + + ! + ! do tracer transport with 'bulk' plume + ! + +! DO n=1,NTR + ! Find tracer flux profile +! DO k=KTS-1,KTE +! trflx( +! END DO + + ! Add tracer tendency directly to bundle + +! END DO + ! ! writing updraft properties for output ! all variables, except Areas are now multipled by the area @@ -826,7 +855,7 @@ SUBROUTINE RUN_EDMF(its,ite, jts,jte, kts,kte, dt, & ! Inputs + mapl_alhl*( UPQL(K,i) - QLI(K) ) & + mapl_alhs*( UPQI(K,I) - QII(K) ) end if - ltm=exfh(k)*(UPTHL(K,i)-THLI(K)) !+mapl_grav*zw(k)/mapl_cp + ltm=exfh(k)*(UPTHL(K,i)-THLI(K)) s_aws(k) = s_aws(K)+UPA(K,i)*UPW(K,i)*tmp ! for trisolver s_ahl2(k) = s_ahl2(K)+UPA(K,i)*ltm*ltm s_ahl3(k) = s_ahl3(K)+UPA(K,i)*ltm*ltm*ltm @@ -942,7 +971,7 @@ subroutine calc_mf_depth(kts,kte,t,z,q,p,ztop,wthv,wqt) thstar=max(0.,wthv)/wstar ! sigmaW=MFPARAMS%AlphaW*wstar - sigmaQT=1.0*qstar + sigmaQT=2.0*qstar sigmaTH=2.0*thstar ! print *,'sigQT=',sigmaQT,' sigTH=',sigmaTH,' wstar=',wstar @@ -961,8 +990,10 @@ subroutine calc_mf_depth(kts,kte,t,z,q,p,ztop,wthv,wqt) tep = tep - MAPL_GRAV*( z2-z1 )/MAPL_CP - qp = qp + (0.25/300.)*(z2-z1)*(q(k)-qp) - tep = tep + (0.25/300.)*(z2-z1)*(t(k)-tep) +! qp = qp + (0.25/300.)*(z2-z1)*(q(k)-qp) +! tep = tep + (0.25/300.)*(z2-z1)*(t(k)-tep) + qp = qp + (0.7/1000.)*(z2-z1)*(q(k)-qp) ! assume fractional entrainment rate of 0.75/km + tep = tep + (0.7/1000.)*(z2-z1)*(t(k)-tep) ! print *,'mfdepth: tep=',tep,' pp=',pp dqsp = GEOS_DQSAT(tep , pp , qsat=qsp, pascals=.true. ) @@ -973,7 +1004,7 @@ subroutine calc_mf_depth(kts,kte,t,z,q,p,ztop,wthv,wqt) ! compare Tv env vs parcel - if ( t2*(1.+MAPL_VIREPS*q(k)) .ge. tep*(1.+MAPL_VIREPS*qp)+0.1 ) then + if ( t2*(1.+MAPL_VIREPS*q(k)) .ge. tep*(1.+MAPL_VIREPS*qp)+0.2 ) then ztop = 0.5*(z2+z1) exit end if diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSturbulence_GridComp/shoc.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSturbulence_GridComp/shoc.F90 index 3f3e1d0bd..bb552ea93 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSturbulence_GridComp/shoc.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSturbulence_GridComp/shoc.F90 @@ -60,7 +60,7 @@ subroutine run_shoc( nx, ny, nzm, nz, dtn, & ! in smixt_inv, lmix_out, smixt1_inv, & ! out smixt2_inv,smixt3_inv, & ! out ! bruntmst_inv, bruntedg_inv, ri_inv, prnum_inv, & ! out - bruntmst_inv, ri_inv, prnum_inv, & ! out + bruntmst_inv, bruntdry_inv, bruntedg_inv, ri_inv, prnum_inv, & ! out shocparams ) @@ -118,7 +118,8 @@ subroutine run_shoc( nx, ny, nzm, nz, dtn, & ! in real, dimension(:,:,:), pointer :: smixt2_inv ! length scale, term 2 real, dimension(:,:,:), pointer :: smixt3_inv ! length scale, term 3 real, dimension(:,:,:), pointer :: bruntmst_inv ! moist Brunt vaisala frequency -! real, dimension(:,:,:), pointer :: bruntedg_inv ! Brunt vaisala frequency on edges + real, dimension(:,:,:), pointer :: bruntdry_inv ! Brunt vaisala frequency on edges + real, dimension(:,:,:), pointer :: bruntedg_inv ! Brunt vaisala frequency on edges real, dimension(:,:,:), pointer :: ri_inv real, dimension(:,:,:), pointer :: prnum_inv @@ -182,8 +183,8 @@ subroutine run_shoc( nx, ny, nzm, nz, dtn, & ! in real isotropy (nx,ny,nzm) ! "Return-to-isotropy" eddy dissipation time scale, s real brunt (nx,ny,nzm) ! Moist Brunt-Vaisalla frequency, s^-1 - real, dimension(nx,ny,nzm) :: total_water, brunt2, def2, thv, brunt_smooth - real, dimension(nx,ny,nz) :: brunt_edge + real, dimension(nx,ny,nzm) :: total_water, brunt2, def2, thv, brunt_smooth, l_par + real, dimension(nx,ny,nz) :: brunt_edge, brunt_dry real, dimension(nx,ny) :: l_inf, l_mix, zcb, lts!, l_par!, denom, numer, cldarr @@ -192,7 +193,7 @@ subroutine run_shoc( nx, ny, nzm, nz, dtn, & ! in tkes, pval, pkap, thlsec, qwsec, & qwthlsec, wqwsec, wthlsec, dum, sm, & prespot, wrk, wrk1, wrk2, wrk3, & - tkeavg, dtqw, dtqi, l_par + tkeavg, dtqw, dtqi !, l_par integer i,j,k,km1,ku,kd,ka,kb,kinv,strt,fnsh,cnvl integer, dimension(nx,ny) :: cldbasek @@ -254,7 +255,7 @@ subroutine run_shoc( nx, ny, nzm, nz, dtn, & ! in prespot = (100000.0*wrk) ** kapa ! Exner function bet(i,j,k) = ggr/(tabs(i,j,k)*prespot) ! Moorthi thv(i,j,k) = thv(i,j,k)*prespot ! Moorthi -! +! ! Lapse rate * height = reference temperature gamaz(i,j,k) = gocp * zl(i,j,k) @@ -264,6 +265,7 @@ subroutine run_shoc( nx, ny, nzm, nz, dtn, & ! in enddo enddo enddo +! print *,'thv=',thv(i,j,1:20) ! Define vertical grid increments for later use in the vertical differentiation @@ -306,7 +308,8 @@ subroutine run_shoc( nx, ny, nzm, nz, dtn, & ! in if (associated(smixt2_inv)) smixt2_inv(:,:,1:nzm) = smixt2(:,:,nzm:1:-1) if (associated(smixt3_inv)) smixt3_inv(:,:,1:nzm) = smixt3(:,:,nzm:1:-1) -! if (associated(bruntedg_inv)) bruntedg_inv(:,:,0:nzm) = brunt_edge(:,:,nz:1:-1) + if (associated(bruntedg_inv)) bruntedg_inv(:,:,0:nzm) = brunt_edge(:,:,nz:1:-1) + if (associated(bruntdry_inv)) bruntdry_inv(:,:,0:nzm) = l_par(:,:,nz:1:-1) if (associated(bruntmst_inv)) bruntmst_inv(:,:,1:nzm) = brunt(:,:,nzm:1:-1) if (associated(prnum_inv)) prnum_inv(:,:,0:nz-1) = prnum(:,:,nz:1:-1) if (associated(ri_inv)) ri_inv(:,:,0:nz-1) = ri(:,:,nz:1:-1) @@ -371,6 +374,7 @@ subroutine tke_shoc() wrk = 0.5 * (tkh(i,j,ku)+tkh(i,j,kd)) if (shocparams%BUOYOPT==2) then +! if (zl(i,j,k).lt.1000.) wthv_sec(i,j,k) = max(0.,wthv_sec(i,j,k)) ! prevent decoupling a_prod_bu = (ggr / thv(i,j,k)) * wthv_sec(i,j,k) else a_prod_bu = -1.*wrk*brunt(i,j,k) + (ggr / thv(i,j,k))*wthv_mf(i,j,k) @@ -404,11 +408,11 @@ subroutine tke_shoc() do itr=1,nitr ! iterate for implicit solution wtke = min(max(min_tke, wtke), max_tke) a_diss = wrk*sqrt(wtke) ! Coefficient in the TKE dissipation term - if (a_diss.ne.-1.) then +! if (a_diss.ne.-1.) then wtke = wrk1 / (1.+a_diss) - else - wtke = wrk1 / (1.01+a_diss) - end if +! else +! wtke = wrk1 / (1.01+a_diss) +! end if wtke = tkef1*wtke + tkef2*wtk2 ! tkef1+tkef2 = 1.0 wtk2 = wtke @@ -449,6 +453,7 @@ subroutine tke_shoc() ! if (k.ge.cldbasek(i,j)) then ! isotropy(i,j,k) = min(200.+(0.5+0.5*tanh(0.3*(lts(i,j)-19.)))*(max_eddy_dissipation_time_scale-200.),isotropy(i,j,k)) ! end if +! if (zl(i,j,k).gt.500.) isotropy(i,j,k) = 200. if (tke(i,j,k).lt.2e-4) isotropy(i,j,k) = 30. wrk1 = ck / prnum(i,j,k) @@ -625,6 +630,7 @@ subroutine eddy_length() ! brunt_edge(i,j,k) = (2.*ggr/(thv(i,j,k)+thv(i,j,kb)))*(thv(i,j,k)-thv(i,j,kb))/adzi(i,j,k) !g/thv/dz *(thv-thv) + brunt_dry(i,j,k) = (thv(i,j,kc)-thv(i,j,kb))*betdz ! Reinitialize the mixing length related arrays to zero smixt(i,j,k) = 1.0 ! shoc_mod module variable smixt @@ -675,8 +681,9 @@ subroutine eddy_length() kk = kk+1 end do dum = (thv(i,j,kk-1)-thv(i,j,kk-2)) + if (abs(dum) .gt. 1e-3) then - l_mix(i,j) = max(zl(i,j,kk-1)+(thv(i,j,3)+0.4-thv(i,j,kk-1))*(zl(i,j,kk-1)-zl(i,j,kk-2))/dum,100.) + l_mix(i,j) = max(zl(i,j,kk-1)+0.*(thv(i,j,3)+0.4-thv(i,j,kk-1))*(zl(i,j,kk-1)-zl(i,j,kk-2))/dum,100.) else l_mix(i,j) = max(zl(i,j,kk-1),100.) end if @@ -815,6 +822,18 @@ subroutine eddy_length() brunt_edge(:,:,1) = brunt_edge(:,:,2) brunt_edge(:,:,nz) = brunt_edge(:,:,nzm) + do j=1,ny + do i=1,nx + kk = 1 + l_mix(i,j) = 0. + do while (l_mix(i,j).lt.1e-5 .and. kk.lt.nzm) + l_mix(i,j) = l_mix(i,j) + brunt(i,j,kk)*adzl(i,j,kk) + kk = kk+1 + end do + l_mix(i,j) = zl(i,j,kk) + end do + end do +! print *,'lmix=',l_mix(i,j) ! brunt_dry = max( bruntmin, brunt_dry ) @@ -829,6 +848,7 @@ subroutine eddy_length() do j=1,ny do i=1,nx + ! if (k.eq.1) print *,'l_mix=',l_mix(i,j) tkes = sqrt(tke(i,j,k)) @@ -840,24 +860,24 @@ subroutine eddy_length() !---------------------------------- ! calculate parcel mixing length !---------------------------------- -! kk = k -! tep = thv(i,j,k)+max(0.,wthv_sec(i,j,k)) / (sqrt(2.)*tkes) ! upward T perturbation -! do while (tep .gt. thv(i,j,kk) .and. kk.lt.nzm) -! kk = kk+1 -! end do -! ktop = kk -! l_par(i,j) = zl(i,j,ktop) !+ (thv(i,j,ktop-1)-tep)* & -!! ((zl(i,j,ktop)-zl(i,j,ktop-1)))/(thv(i,j,ktop)-thv(i,j,ktop-1)) -! kk = k -! tep = thv(i,j,k)-max(0.,wthv_sec(i,j,k)) / (sqrt(2.)*tkes) ! downward T perturbation -! do while (tep .lt. thv(i,j,kk) .and. kk .gt. 1) -! kk = kk-1 -! end do -! l_par(i,j) = l_par(i,j) - zl(i,j,kk) !- (thv(i,j,kk+1)-tep)* & -!! ((zl(i,j,kk)-zl(i,j,kk+1)))/(thv(i,j,kk)-thv(i,j,kk+1)) -! l_par(i,j) = max(min(l_par(i,j),1500.),25.) + kk = k + wrk = thv(i,j,k)+0.2 !max(0.001,3.*wthv_sec(i,j,k)) / max(0.05,sqrt(0.667)*tkes) ! upward T perturbation + do while (wrk .gt. thv(i,j,kk+1) .and. kk.lt.nzm) + kk = kk+1 + end do + l_par(i,j,k) = zl(i,j,kk) + max(0.,(wrk-thv(i,j,kk))* & + (zl(i,j,kk+1)-zl(i,j,kk)) / (thv(i,j,kk+1)-thv(i,j,kk))) + kk = k + wrk = thv(i,j,k)-0.2 !max(0.001,3.*wthv_sec(i,j,k)) / max(0.05,sqrt(0.667)*tkes) ! downward T perturbation + do while (wrk .lt. thv(i,j,kk-1) .and. kk .gt. 1) + kk = kk-1 + end do + l_par(i,j,k) = l_par(i,j,k) - zl(i,j,kk) + max(0.,(thv(i,j,kk)-wrk)* & + (zl(i,j,kk)-zl(i,j,kk-1))/(thv(i,j,kk)-thv(i,j,kk-1))) + l_par(i,j,k) = max(min(l_par(i,j,k),1500.),25.) +! if (zl(i,j,k).lt.600.) print *,'z=',zl(i,j,k),' thv=',thv(i,j,k),' thv pert=',max(0.001,3.*wthv_sec(i,j,k)) / max(0.05,sqrt(0.667)*tkes),' l_par=',l_par(i,j,k) + ! if (brunt_smooth(i,j,k).gt.1e-5) l_par(i,j) = max(25.,l_par(i,j)/2.) -! wrk2 = (tscale*tkes*0.1*l_par(i,j)) !---------------------------------- ! calculate 'TKE' mixing length @@ -884,12 +904,14 @@ subroutine eddy_length() ! smixt1(i,j,k) = sqrt(400.*tkes*vonk*zl(i,j,k))*shocparams%LENFAC1 ! Turbulent length scale - if (zl(i,j,k).lt.l_mix(i,j)) then +! if (zl(i,j,k).lt.zl(i,j,cldbasek(i,j))) then +! if (zl(i,j,k).lt.l_mix(i,j)) then ! smixt2(i,j,k) = sqrt(0.1*zpbl(i,j)*400.*tkes)*shocparams%LENFAC2 - smixt2(i,j,k) = sqrt(min(1000.,l_mix(i,j))*400.*tkes)*(shocparams%LENFAC2) - else - smixt2(i,j,k) = 400.*tkes*shocparams%LENFAC2 - end if +! smixt2(i,j,k) = sqrt(min(1000.,l_mix(i,j))*400.*tkes)*(shocparams%LENFAC2) + smixt2(i,j,k) = sqrt(l_par(i,j,k)*400.*tkes)*(shocparams%LENFAC2) +! else +! smixt2(i,j,k) = 400.*tkes*shocparams%LENFAC2 +! end if ! Stability length scale smixt3(i,j,k) = max(0.1,tkes)*shocparams%LENFAC3/(sqrt(brunt_smooth(i,j,k))) @@ -1302,9 +1324,9 @@ subroutine update_moments( IM, JM, LM, & ! in wrk1b = bet2*isosqr f2 = thedz*wrk1b*whl_can(i,j,k)*0.667*(TKE(i,j,k)-TKE(i,j,kb)) & - + (thedz2+thedz2)*bet(i,j,k)*isosqr*wrk + + (thedz2+thedz2)*bet(i,j,k)*isosqr*avew*wrk - f3 = thedz2*wrk1b*wrk + thedz*bet2*isosqr*(whl_can(i,j,k)*(tke(i,j,k)-tke(i,j,kb))) + f3 = thedz2*wrk1b*wrk*avew + thedz*bet2*isosqr*(whl_can(i,j,k)*(tke(i,j,k)-tke(i,j,kb))) wrk1b = thedz*iso*avew f4 = wrk1b*(0.667*TKE(i,j,k)-0.667*TKE(i,j,kb) + tke(i,j,k)-tke(i,j,kb)) @@ -1345,7 +1367,6 @@ subroutine update_moments( IM, JM, LM, & ! in if (abs(dum).le.1e-16) dum = sign(1e-16,dum) ! if (abs(dum).eq.1e-16) print *,'c-1.2*X0+AA0=',dum w3can(i,j,k) = max(-cond, min(cond, (AA1-1.2*X1-1.5*f5)/dum)) - ! Implemetation of the C01 approach in this subroutine is nearly complete ! (the missing part are Eqs. 5c and 5e which are very simple) ! therefore it's easy to diagnose other third order moments obtained in C01 using this code. @@ -1358,10 +1379,11 @@ subroutine update_moments( IM, JM, LM, & ! in w3can(i,j,LM) = w3can(i,j,LM-1) enddo enddo + w3 = w3can !! skew_w = w3 / w2**1.5 -! qt3 = 1.2*w3*(qt2/w2)**1.5 -! hl3 = w3 * (hl2 / w2)**1.5 + qt3 = 1.2*w3*(qt2/w2)**1.5 + hl3 = w3 * (hl2 / w2)**1.5 end if ! DOCANUTO conditional From f5a2138279aa17c34b63a8e0710ba5b28ee45741 Mon Sep 17 00:00:00 2001 From: Nathan Arnold Date: Thu, 1 Feb 2024 10:04:41 -0500 Subject: [PATCH 2/8] MFQT3 and MFHL3 set to 0 below MF cloud base. no QT3 diffusion. PWMIN to 1.2. total TKE in KH calc, and use EDMF_TKE as min limit. --- .../GEOS_PhysicsGridComp.F90 | 12 ++++----- .../GEOSturbulence_GridComp/edmf.F90 | 9 +++++-- .../GEOSturbulence_GridComp/shoc.F90 | 25 ++++++------------- 3 files changed, 20 insertions(+), 26 deletions(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOS_PhysicsGridComp.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOS_PhysicsGridComp.F90 index e7d1336b9..2bcaba087 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOS_PhysicsGridComp.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOS_PhysicsGridComp.F90 @@ -1698,12 +1698,12 @@ subroutine Initialize ( GC, IMPORT, EXPORT, CLOCK, RC ) call MAPL_FieldBundleAdd (BUNDLE, FIELD, RC=STATUS ) VERIFY_(STATUS) - call ESMF_StateGet (GEX(TURBL), 'QT3' , FIELD, RC=STATUS ) - VERIFY_(STATUS) - call ESMF_AttributeSet(FIELD, NAME="DiffuseLike" ,VALUE="Q", RC=STATUS ) - VERIFY_(STATUS) - call MAPL_FieldBundleAdd (BUNDLE, FIELD, RC=STATUS ) - VERIFY_(STATUS) +! call ESMF_StateGet (GEX(TURBL), 'QT3' , FIELD, RC=STATUS ) +! VERIFY_(STATUS) +! call ESMF_AttributeSet(FIELD, NAME="DiffuseLike" ,VALUE="Q", RC=STATUS ) +! VERIFY_(STATUS) +! call MAPL_FieldBundleAdd (BUNDLE, FIELD, RC=STATUS ) +! VERIFY_(STATUS) ! Add Friendlies from Physics call MAPL_GridCompGetFriendlies(GC, "TURBULENCE", BUNDLE, RC=STATUS ) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSturbulence_GridComp/edmf.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSturbulence_GridComp/edmf.F90 index e1d96d8e1..f6df882d7 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSturbulence_GridComp/edmf.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSturbulence_GridComp/edmf.F90 @@ -929,9 +929,14 @@ SUBROUTINE RUN_EDMF(its,ite, jts,jte, kts,kte, dt, & ! Inputs buoyf(IH,JH,K) = s_buoyf(KTE+KTS-K) ! can be used in SHOC mfw2(IH,JH,K) = 0.5*(s_aw2(KTE+KTS-K-1)+s_aw2(KTE+KTS-K)) mfw3(IH,JH,K) = 0.5*(s_aw3(KTE+KTS-K-1)+s_aw3(KTE+KTS-K)) - mfqt3(IH,JH,K) = 0.5*(s_aqt3(KTE+KTS-K-1)+s_aqt3(KTE+KTS-K)) - mfhl3(IH,JH,K) = 0.5*(s_ahl3(KTE+KTS-K-1)+s_ahl3(KTE+KTS-K)) mfhlqt(IH,JH,K) = 0.5*(s_ahlqt(KTE+KTS-K-1)+s_ahlqt(KTE+KTS-K)) + if (SUM(moist_a(KTS-1:KTE+KTS-K)).le.1e-4) then + mfqt3(IH,JH,K) = 0. + mfhl3(IH,JH,K) = 0. + else + mfqt3(IH,JH,K) = 0.5*(s_aqt3(KTE+KTS-K-1)+s_aqt3(KTE+KTS-K)) + mfhl3(IH,JH,K) = 0.5*(s_ahl3(KTE+KTS-K-1)+s_ahl3(KTE+KTS-K)) + end if ! mfhl2(IH,JH,K)=0.5*(s_ahl2(KTE+KTS-K-1)+s_ahl2(KTE+KTS-K)) ! no longer needed ! mfqt2(IH,JH,K)=0.5*(s_aqt2(KTE+KTS-K-1)+s_aqt2(KTE+KTS-K)) ! no longer needed ENDDO diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSturbulence_GridComp/shoc.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSturbulence_GridComp/shoc.F90 index bb552ea93..44d2d330f 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSturbulence_GridComp/shoc.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSturbulence_GridComp/shoc.F90 @@ -374,7 +374,6 @@ subroutine tke_shoc() wrk = 0.5 * (tkh(i,j,ku)+tkh(i,j,kd)) if (shocparams%BUOYOPT==2) then -! if (zl(i,j,k).lt.1000.) wthv_sec(i,j,k) = max(0.,wthv_sec(i,j,k)) ! prevent decoupling a_prod_bu = (ggr / thv(i,j,k)) * wthv_sec(i,j,k) else a_prod_bu = -1.*wrk*brunt(i,j,k) + (ggr / thv(i,j,k))*wthv_mf(i,j,k) @@ -405,8 +404,9 @@ subroutine tke_shoc() wrk = (dtn*Cee)/smixt(i,j,k) wrk1 = wtke + dtn*(a_prod_sh+a_prod_bu) + wrk2 = min_tke+0.5*(tke_mf(i,j,nz-k+1)+tke_mf(i,j,nz-k)) do itr=1,nitr ! iterate for implicit solution - wtke = min(max(min_tke, wtke), max_tke) + wtke = min(max(wrk2, wtke), max_tke) a_diss = wrk*sqrt(wtke) ! Coefficient in the TKE dissipation term ! if (a_diss.ne.-1.) then wtke = wrk1 / (1.+a_diss) @@ -418,7 +418,7 @@ subroutine tke_shoc() enddo - tke(i,j,k) = min(max(min_tke, wtke), max_tke) + tke(i,j,k) = min(max(wrk2, wtke), max_tke) tscale1(i,j,k) = (dtn+dtn) / a_diss ! See Eq 8 in BK13 (note typo, flipped num/denom) @@ -453,7 +453,6 @@ subroutine tke_shoc() ! if (k.ge.cldbasek(i,j)) then ! isotropy(i,j,k) = min(200.+(0.5+0.5*tanh(0.3*(lts(i,j)-19.)))*(max_eddy_dissipation_time_scale-200.),isotropy(i,j,k)) ! end if -! if (zl(i,j,k).gt.500.) isotropy(i,j,k) = 200. if (tke(i,j,k).lt.2e-4) isotropy(i,j,k) = 30. wrk1 = ck / prnum(i,j,k) @@ -461,9 +460,9 @@ subroutine tke_shoc() ! tkh(i,j,k) = 0.5*( smixt(i,j,k)*sqrt(tke(i,j,k)) & ! alternate form ! + smixt(i,j,k-1)*sqrt(tke(i,j,k-1)) ) - tke_env = max(min_tke,0.5*(tke(i,j,k)+tke(i,j,k-1))-tke_mf(i,j,nz-k+1)) - tkh(i,j,k) = wrk1*isotropy(i,j,k) & - * 1.*(tke_env) ! remove MF TKE +! tke_env = max(min_tke,0.5*(tke(i,j,k)+tke(i,j,k-1))-0.*tke_mf(i,j,nz-k+1)) + tkh(i,j,k) = wrk1*isotropy(i,j,k)*0.5*(tke(i,j,k)+tke(i,j,k-1)) ! & +! *(tke_env) ! remove MF TKE tkh(i,j,k) = min(tkh(i,j,k),tkhmax) end do ! i end do ! j @@ -877,7 +876,6 @@ subroutine eddy_length() l_par(i,j,k) = max(min(l_par(i,j,k),1500.),25.) ! if (zl(i,j,k).lt.600.) print *,'z=',zl(i,j,k),' thv=',thv(i,j,k),' thv pert=',max(0.001,3.*wthv_sec(i,j,k)) / max(0.05,sqrt(0.667)*tkes),' l_par=',l_par(i,j,k) -! if (brunt_smooth(i,j,k).gt.1e-5) l_par(i,j) = max(25.,l_par(i,j)/2.) !---------------------------------- ! calculate 'TKE' mixing length @@ -899,19 +897,11 @@ subroutine eddy_length() if ( shocparams%LENOPT .lt. 4 ) then ! SHOC-MF length scale ! Surface length scale -! smixt1(i,j,k) = vonk*zl(i,j,k)*exp(MIN(1000.,zl(i,j,k))**2/4e4)*shocparams%LENFAC1 smixt1(i,j,k) = vonk*zl(i,j,k)*shocparams%LENFAC1 -! smixt1(i,j,k) = sqrt(400.*tkes*vonk*zl(i,j,k))*shocparams%LENFAC1 +! smixt1(i,j,k) = sqrt(400.*tkes*vonk*zl(i,j,k))*shocparams%LENFAC1 ! original SHOC, includes TKE ! Turbulent length scale -! if (zl(i,j,k).lt.zl(i,j,cldbasek(i,j))) then -! if (zl(i,j,k).lt.l_mix(i,j)) then -! smixt2(i,j,k) = sqrt(0.1*zpbl(i,j)*400.*tkes)*shocparams%LENFAC2 -! smixt2(i,j,k) = sqrt(min(1000.,l_mix(i,j))*400.*tkes)*(shocparams%LENFAC2) smixt2(i,j,k) = sqrt(l_par(i,j,k)*400.*tkes)*(shocparams%LENFAC2) -! else -! smixt2(i,j,k) = 400.*tkes*shocparams%LENFAC2 -! end if ! Stability length scale smixt3(i,j,k) = max(0.1,tkes)*shocparams%LENFAC3/(sqrt(brunt_smooth(i,j,k))) @@ -1257,7 +1247,6 @@ subroutine update_moments( IM, JM, LM, & ! in qt3 = ( qt3 + max(MFQT3,0.) ) / ( 1. + DT/QT3_TSCALE ) hl3 = MFHL3 w3 = MFW3 -! w3 = (1.-MFFRC)*0.5*((Sqrt(w2)-MFAW/(1.-MFFRC))**3 + (-Sqrt(w2)-MFAW/(1.-MFFRC))**3) + MFW3 else ! pre-define adzl, From 71745acef077803972b7e7a0a9720968962d2bcd Mon Sep 17 00:00:00 2001 From: Nathan Arnold Date: Tue, 6 Feb 2024 09:53:56 -0500 Subject: [PATCH 3/8] EDMF parameter updates, added PDFDIAG ifdef to micro interfaces --- .../GEOS_BACM_1M_InterfaceMod.F90 | 31 +++++++++++++ .../GEOS_GFDL_1M_InterfaceMod.F90 | 43 +++++++++++++++++++ .../GEOS_MGB2_2M_InterfaceMod.F90 | 42 +++++++++++++++++- .../GEOS_THOM_1M_InterfaceMod.F90 | 43 +++++++++++++++++++ .../GEOSmoist_GridComp/Process_Library.F90 | 36 ++++++++-------- .../GEOS_TurbulenceGridComp.F90 | 6 +-- 6 files changed, 179 insertions(+), 22 deletions(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/GEOS_BACM_1M_InterfaceMod.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/GEOS_BACM_1M_InterfaceMod.F90 index 6551c126b..73894bc52 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/GEOS_BACM_1M_InterfaceMod.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/GEOS_BACM_1M_InterfaceMod.F90 @@ -1,6 +1,7 @@ ! $Id$ #include "MAPL_Generic.h" +!#define PDFDIAG 1 !============================================================================= !BOP @@ -348,6 +349,12 @@ subroutine BACM_1M_Run (GC, IMPORT, EXPORT, CLOCK, RC) real, pointer, dimension(:,: ) :: LS_SNR, CN_SNR, AN_SNR, SC_SNR real, pointer, dimension(:,:,:) :: PTR3D real, pointer, dimension(:,: ) :: PTR2D +#ifdef PDFDIAG + real, pointer, dimension(:,:,:) :: PDF_W1, PDF_W2, PDF_SIGW1, PDF_SIGW2, & + PDF_QT1, PDF_QT2, PDF_SIGQT1, PDF_SIGQT2, & + PDF_TH1, PDF_TH2, PDF_SIGTH1, PDF_SIGTH2, & + PDF_RQTTH, PDF_RWTH, PDF_RWQT +#endif call ESMF_GridCompGet( GC, CONFIG=CF, RC=STATUS ) VERIFY_(STATUS) @@ -591,6 +598,24 @@ subroutine BACM_1M_Run (GC, IMPORT, EXPORT, CLOCK, RC) call MAPL_GetPointer(EXPORT, WTHV2, 'WTHV2' , ALLOC=.TRUE., RC=STATUS); VERIFY_(STATUS) call MAPL_GetPointer(EXPORT, WQL, 'WQL' , ALLOC=.TRUE., RC=STATUS); VERIFY_(STATUS) +#ifdef PDFDIAG + call MAPL_GetPointer(EXPORT, PDF_W1, 'PDF_W1' , ALLOC=.TRUE., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT, PDF_W2, 'PDF_W2' , ALLOC=.TRUE., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT, PDF_SIGW1, 'PDF_SIGW1' , ALLOC=.TRUE., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT, PDF_SIGW2, 'PDF_SIGW2' , ALLOC=.TRUE., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT, PDF_QT1, 'PDF_QT1' , ALLOC=.TRUE., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT, PDF_QT2, 'PDF_QT2' , ALLOC=.TRUE., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT, PDF_SIGQT1, 'PDF_SIGQT1' , ALLOC=.TRUE., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT, PDF_SIGQT2, 'PDF_SIGQT2' , ALLOC=.TRUE., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT, PDF_TH1, 'PDF_TH1' , ALLOC=.TRUE., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT, PDF_TH2, 'PDF_TH2' , ALLOC=.TRUE., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT, PDF_SIGTH1, 'PDF_SIGTH1' , ALLOC=.TRUE., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT, PDF_SIGTH2, 'PDF_SIGTH2' , ALLOC=.TRUE., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT, PDF_RQTTH, 'PDF_RQTTH' , ALLOC=.TRUE., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT, PDF_RWTH, 'PDF_RWTH' , ALLOC=.TRUE., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT, PDF_RWQT, 'PDF_RWQT' , ALLOC=.TRUE., RC=STATUS); VERIFY_(STATUS) +#endif + call PROGNO_CLOUD ( & IM*JM, LM , & DT_MOIST , & @@ -689,6 +714,12 @@ subroutine BACM_1M_Run (GC, IMPORT, EXPORT, CLOCK, RC) VFALLRN_AN , VFALLRN_LS ,VFALLRN_CN ,VFALLRN_SC , & PDF_A, PDFITERS, & DQVDT_macro, DQLDT_macro, DQIDT_macro, DQADT_macro, & +#ifdef PDFDIAG + PDF_SIGW1, PDF_SIGW2, PDF_W1, PDF_W2, & + PDF_SIGTH1, PDF_SIGTH2, PDF_TH1, PDF_TH2, & + PDF_SIGQT1, PDF_SIGQT2, PDF_QT1, PDF_QT2, & + PDF_RQTTH, PDF_RWTH, PDF_RWQT, & +#endif WTHV2, WQL, & NACTL, & NACTI, & diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/GEOS_GFDL_1M_InterfaceMod.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/GEOS_GFDL_1M_InterfaceMod.F90 index 50c5c0c1f..8b5140755 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/GEOS_GFDL_1M_InterfaceMod.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/GEOS_GFDL_1M_InterfaceMod.F90 @@ -1,6 +1,7 @@ ! $Id$ #include "MAPL_Generic.h" +!#define PDFDIAG 1 !============================================================================= !BOP @@ -342,6 +343,12 @@ subroutine GFDL_1M_Run (GC, IMPORT, EXPORT, CLOCK, RC) real, pointer, dimension(:,:) :: DBZ_MAX, DBZ_1KM, DBZ_TOP, DBZ_M10C real, pointer, dimension(:,:,:) :: PTR3D real, pointer, dimension(:,: ) :: PTR2D +#ifdef PDFDIAG + real, pointer, dimension(:,:,:) :: PDF_W1, PDF_W2, PDF_SIGW1, PDF_SIGW2, & + PDF_QT1, PDF_QT2, PDF_SIGQT1, PDF_SIGQT2, & + PDF_TH1, PDF_TH2, PDF_SIGTH1, PDF_SIGTH2, & + PDF_RQTTH, PDF_RWTH, PDF_RWQT +#endif ! Local variables real :: facEIS @@ -539,6 +546,25 @@ subroutine GFDL_1M_Run (GC, IMPORT, EXPORT, CLOCK, RC) DQSDT_macro=QSNOW DQGDT_macro=QGRAUPEL +#ifdef PDFDIAG + call MAPL_GetPointer(EXPORT, PDF_W1, 'PDF_W1' , ALLOC=.TRUE., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT, PDF_W2, 'PDF_W2' , ALLOC=.TRUE., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT, PDF_SIGW1, 'PDF_SIGW1' , ALLOC=.TRUE., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT, PDF_SIGW2, 'PDF_SIGW2' , ALLOC=.TRUE., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT, PDF_QT1, 'PDF_QT1' , ALLOC=.TRUE., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT, PDF_QT2, 'PDF_QT2' , ALLOC=.TRUE., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT, PDF_SIGQT1, 'PDF_SIGQT1' , ALLOC=.TRUE., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT, PDF_SIGQT2, 'PDF_SIGQT2' , ALLOC=.TRUE., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT, PDF_TH1, 'PDF_TH1' , ALLOC=.TRUE., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT, PDF_TH2, 'PDF_TH2' , ALLOC=.TRUE., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT, PDF_SIGTH1, 'PDF_SIGTH1' , ALLOC=.TRUE., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT, PDF_SIGTH2, 'PDF_SIGTH2' , ALLOC=.TRUE., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT, PDF_RQTTH, 'PDF_RQTTH' , ALLOC=.TRUE., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT, PDF_RWTH, 'PDF_RWTH' , ALLOC=.TRUE., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT, PDF_RWQT, 'PDF_RWQT' , ALLOC=.TRUE., RC=STATUS); VERIFY_(STATUS) +#endif + + ! Include shallow precip condensates if present call MAPL_GetPointer(EXPORT, PTR3D, 'SHLW_PRC3', ALLOC=.TRUE., RC=STATUS); VERIFY_(STATUS) if (associated(PTR3D)) then @@ -612,6 +638,23 @@ subroutine GFDL_1M_Run (GC, IMPORT, EXPORT, CLOCK, RC) SL3(I,J,L) , & PDF_A(I,J,L) , & PDFITERS(I,J,L), & +#ifdef PDFDIAG + PDF_SIGW1(I,J,L), & + PDF_SIGW2(I,J,L), & + PDF_W1(I,J,L), & + PDF_W2(I,J,L), & + PDF_SIGTH1(I,J,L), & + PDF_SIGTH2(I,J,L), & + PDF_TH1(I,J,L), & + PDF_TH2(I,J,L), & + PDF_SIGQT1(I,J,L), & + PDF_SIGQT2(I,J,L), & + PDF_QT1(I,J,L), & + PDF_QT2(I,J,L), & + PDF_RQTTH(I,J,L), & + PDF_RWTH(I,J,L), & + PDF_RWQT(I,J,L), & +#endif WTHV2(I,J,L) , & WQL(I,J,L) , & .false. , & diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/GEOS_MGB2_2M_InterfaceMod.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/GEOS_MGB2_2M_InterfaceMod.F90 index faa49a587..7d2070ebf 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/GEOS_MGB2_2M_InterfaceMod.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/GEOS_MGB2_2M_InterfaceMod.F90 @@ -1,6 +1,7 @@ ! $Id$ #include "MAPL_Generic.h" +!#define PDFDIAG 1 !============================================================================= !BOP @@ -507,7 +508,12 @@ subroutine MGB2_2M_Run (GC, IMPORT, EXPORT, CLOCK, RC) real, pointer, dimension(:,:,:) :: RHCRIT3D real, pointer, dimension(:,:,:) :: PTR3D real, pointer, dimension(:,: ) :: PTR2D - +#ifdef PDFDIAG + real, pointer, dimension(:,:,:) :: PDF_W1, PDF_W2, PDF_SIGW1, PDF_SIGW2, & + PDF_QT1, PDF_QT2, PDF_SIGQT1, PDF_SIGQT2, & + PDF_TH1, PDF_TH2, PDF_SIGTH1, PDF_SIGTH2, & + PDF_RQTTH, PDF_RWTH, PDF_RWQT +#endif !2m real, pointer, dimension(:,:,:) :: SC_ICE, CDNC_NUC, INC_NUC, PFRZ, & @@ -1413,6 +1419,23 @@ subroutine MGB2_2M_Run (GC, IMPORT, EXPORT, CLOCK, RC) DQSDT_macro=QSNOW DQGDT_macro=QGRAUPEL +#ifdef PDFDIAG + call MAPL_GetPointer(EXPORT, PDF_W1, 'PDF_W1' , ALLOC=.TRUE., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT, PDF_W2, 'PDF_W2' , ALLOC=.TRUE., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT, PDF_SIGW1, 'PDF_SIGW1' , ALLOC=.TRUE., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT, PDF_SIGW2, 'PDF_SIGW2' , ALLOC=.TRUE., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT, PDF_QT1, 'PDF_QT1' , ALLOC=.TRUE., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT, PDF_QT2, 'PDF_QT2' , ALLOC=.TRUE., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT, PDF_SIGQT1, 'PDF_SIGQT1' , ALLOC=.TRUE., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT, PDF_SIGQT2, 'PDF_SIGQT2' , ALLOC=.TRUE., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT, PDF_TH1, 'PDF_TH1' , ALLOC=.TRUE., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT, PDF_TH2, 'PDF_TH2' , ALLOC=.TRUE., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT, PDF_SIGTH1, 'PDF_SIGTH1' , ALLOC=.TRUE., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT, PDF_SIGTH2, 'PDF_SIGTH2' , ALLOC=.TRUE., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT, PDF_RQTTH, 'PDF_RQTTH' , ALLOC=.TRUE., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT, PDF_RWTH, 'PDF_RWTH' , ALLOC=.TRUE., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT, PDF_RWQT, 'PDF_RWQT' , ALLOC=.TRUE., RC=STATUS); VERIFY_(STATUS) +#endif ! Include shallow precip condensates if present call MAPL_GetPointer(EXPORT, PTR3D, 'SHLW_PRC3', RC=STATUS); VERIFY_(STATUS) @@ -1473,6 +1496,23 @@ subroutine MGB2_2M_Run (GC, IMPORT, EXPORT, CLOCK, RC) SL3(I,J,L) , & PDF_A(I,J,L) , & PDFITERS(I,J,L), & +#ifdef PDFDIAG + PDF_SIGW1(I,J,L), & + PDF_SIGW2(I,J,L), & + PDF_W1(I,J,L), & + PDF_W2(I,J,L), & + PDF_SIGTH1(I,J,L), & + PDF_SIGTH2(I,J,L), & + PDF_TH1(I,J,L), & + PDF_TH2(I,J,L), & + PDF_SIGQT1(I,J,L), & + PDF_SIGQT2(I,J,L), & + PDF_QT1(I,J,L), & + PDF_QT2(I,J,L), & + PDF_RQTTH(I,J,L), & + PDF_RWTH(I,J,L), & + PDF_RWQT(I,J,L), & +#endif WTHV2(I,J,L) , & WQL(I,J,L) , & .false. , & diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/GEOS_THOM_1M_InterfaceMod.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/GEOS_THOM_1M_InterfaceMod.F90 index 460f4e392..9ff2269a0 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/GEOS_THOM_1M_InterfaceMod.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/GEOS_THOM_1M_InterfaceMod.F90 @@ -1,6 +1,7 @@ ! $Id$ #include "MAPL_Generic.h" +!#define PDFDIAG 1 !============================================================================= !BOP @@ -353,6 +354,13 @@ subroutine THOM_1M_Run (GC, IMPORT, EXPORT, CLOCK, RC) real, pointer, dimension(:,:,:) :: DBZ3D real, pointer, dimension(:,:,:) :: PTR3D real, pointer, dimension(:,: ) :: PTR2D +#ifdef PDFDIAG + real, pointer, dimension(:,:,:) :: PDF_W1, PDF_W2, PDF_SIGW1, PDF_SIGW2, & + PDF_QT1, PDF_QT2, PDF_SIGQT1, PDF_SIGQT2, & + PDF_TH1, PDF_TH2, PDF_SIGTH1, PDF_SIGTH2, & + PDF_RQTTH, PDF_RWTH, PDF_RWQT +#endif + ! Thompson Pointers for inputs real, dimension(:,:,:), allocatable, target :: inputs real, dimension(:,:,:), pointer :: qv => null() @@ -679,6 +687,24 @@ subroutine THOM_1M_Run (GC, IMPORT, EXPORT, CLOCK, RC) DQSDT_macro=QSNOW DQGDT_macro=QGRAUPEL +#ifdef PDFDIAG + call MAPL_GetPointer(EXPORT, PDF_W1, 'PDF_W1' , ALLOC=.TRUE., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT, PDF_W2, 'PDF_W2' , ALLOC=.TRUE., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT, PDF_SIGW1, 'PDF_SIGW1' , ALLOC=.TRUE., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT, PDF_SIGW2, 'PDF_SIGW2' , ALLOC=.TRUE., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT, PDF_QT1, 'PDF_QT1' , ALLOC=.TRUE., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT, PDF_QT2, 'PDF_QT2' , ALLOC=.TRUE., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT, PDF_SIGQT1, 'PDF_SIGQT1' , ALLOC=.TRUE., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT, PDF_SIGQT2, 'PDF_SIGQT2' , ALLOC=.TRUE., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT, PDF_TH1, 'PDF_TH1' , ALLOC=.TRUE., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT, PDF_TH2, 'PDF_TH2' , ALLOC=.TRUE., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT, PDF_SIGTH1, 'PDF_SIGTH1' , ALLOC=.TRUE., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT, PDF_SIGTH2, 'PDF_SIGTH2' , ALLOC=.TRUE., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT, PDF_RQTTH, 'PDF_RQTTH' , ALLOC=.TRUE., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT, PDF_RWTH, 'PDF_RWTH' , ALLOC=.TRUE., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT, PDF_RWQT, 'PDF_RWQT' , ALLOC=.TRUE., RC=STATUS); VERIFY_(STATUS) +#endif + ! Include shallow precip condensates if present call MAPL_GetPointer(EXPORT, PTR3D, 'SHLW_PRC3', ALLOC=.TRUE., RC=STATUS); VERIFY_(STATUS) if (associated(PTR3D)) then @@ -751,6 +777,23 @@ subroutine THOM_1M_Run (GC, IMPORT, EXPORT, CLOCK, RC) SL3(I,J,L) , & PDF_A(I,J,L) , & PDFITERS(I,J,L), & +#ifdef PDFDIAG + PDF_SIGW1(I,J,L), & + PDF_SIGW2(I,J,L), & + PDF_W1(I,J,L), & + PDF_W2(I,J,L), & + PDF_SIGTH1(I,J,L), & + PDF_SIGTH2(I,J,L), & + PDF_TH1(I,J,L), & + PDF_TH2(I,J,L), & + PDF_SIGQT1(I,J,L), & + PDF_SIGQT2(I,J,L), & + PDF_QT1(I,J,L), & + PDF_QT2(I,J,L), & + PDF_RQTTH(I,J,L), & + PDF_RWTH(I,J,L), & + PDF_RWQT(I,J,L), & +#endif WTHV2(I,J,L) , & WQL(I,J,L) , & .false. , & diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/Process_Library.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/Process_Library.F90 index 05d575644..96f1c8e1e 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/Process_Library.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/Process_Library.F90 @@ -1905,31 +1905,31 @@ subroutine partition_dblgss( fQi, & ! IN rwthl = 0.0 end if - if (abs(rwqt).gt.1e3) print *,'rwqt>1e3' - if (abs(rwthl).gt.1e3) print *,'rwthl>1e3' - if (abs(C1).gt.1e3) print *,'C1>1e3' - if (abs(C2).gt.1e3) print *,'C2>1e3' - if (abs(w1_1).gt.1e3) print *,'w1_1>1e3' - if (abs(w1_2).gt.1e3) print *,'w1_2>1e3' - if (abs(qw2_2).gt.1e3) print *,'qw2_2>1e3' - if (abs(qw2_1).gt.1e3) print *,'qw2_1>1e3' - if (abs(cqt1).gt.1e3) print *,'cqt1>1e3' - if (abs(cqt2).gt.1e3) print *,'cqt2>1e3' - if (abs(cthl1).gt.1e3) print *,'cthl1>1e3' - if (abs(cthl2).gt.1e3) print *,'cthl2>1e3' - if (abs(thl2_1).gt.1e3) print *,'thl2_1>1e3' - if (abs(thl2_2).gt.1e3) print *,'thl2_2>1e3' + if (isnan(rwqt)) print *,'rwqt>1e3' + if (isnan(rwthl)) print *,'rwthl>1e3' + if (isnan(C1)) print *,'C1>1e3' + if (isnan(C2)) print *,'C2>1e3' + if (isnan(w1_1)) print *,'w1_1>1e3' + if (isnan(w1_2)) print *,'w1_2>1e3' + if (isnan(qw2_2)) print *,'qw2_2>1e3' + if (isnan(qw2_1)) print *,'qw2_1>1e3' + if (isnan(cqt1)) print *,'cqt1>1e3' + if (isnan(cqt2)) print *,'cqt2>1e3' + if (isnan(cthl1)) print *,'cthl1>1e3' + if (isnan(cthl2)) print *,'cthl2>1e3' + if (isnan(thl2_1)) print *,'thl2_1>1e3' + if (isnan(thl2_2)) print *,'thl2_2>1e3' wql1 = C1*(cqt1*sqrt(w2_1)*sqrt(qw2_1)*rwqt-cthl1*sqrt(w2_1)*sqrt(thl2_1)*rwthl) wql2 = C2*(cqt2*sqrt(w2_2)*sqrt(qw2_2)*rwqt-cthl2*sqrt(w2_2)*sqrt(thl2_2)*rwthl) - if (abs(wql1).gt.1e3) print *,'wql1>1e3' - if (abs(wql2).gt.1e3) print *,'wql2>1e3' + if (isnan(wql1)) print *,'wql1>1e3' + if (isnan(wql2)) print *,'wql2>1e3' ! Compute the liquid water flux wqls = aterm * ((w1_1-w_first)*ql1+wql1) + onema * ((w1_2-w_first)*ql2+wql2) wqis = aterm * ((w1_1-w_first)*qi1) + onema * ((w1_2-w_first)*qi2) - if (isnan(wqls)) print *,'wqls is nan at p=',pval - if (wqls>1e3) print *,'wqls >1e3 at p=',pval + if (isnan(wqls)) wqls = 0. + if (isnan(wqis)) wqis = 0. ! diagnostic buoyancy flux. Includes effects from liquid water, ice ! condensate, liquid & ice precipitation diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSturbulence_GridComp/GEOS_TurbulenceGridComp.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSturbulence_GridComp/GEOS_TurbulenceGridComp.F90 index 72817571d..db4c79ff5 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSturbulence_GridComp/GEOS_TurbulenceGridComp.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSturbulence_GridComp/GEOS_TurbulenceGridComp.F90 @@ -3553,7 +3553,7 @@ 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., RC=STATUS) + call MAPL_GetResource (MAPL, MFPARAMS%PWMIN, "EDMF_PWMIN:", default=1.2, RC=STATUS) call MAPL_GetResource (MAPL, MFPARAMS%PWMAX, "EDMF_PWMAX:", default=3., RC=STATUS) ! call MAPL_GetResource (MAPL, MFPARAMS%ENTUFAC, "EDMF_ENTUFAC:", default=1.6, RC=STATUS) @@ -3562,12 +3562,12 @@ subroutine REFRESH(IM,JM,LM,RC) ! 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=1.0, RC=STATUS) - call MAPL_GetResource (MAPL, MFPARAMS%AlphaTH, "EDMF_ALPHATH:", default=1.0, RC=STATUS) + call MAPL_GetResource (MAPL, MFPARAMS%AlphaTH, "EDMF_ALPHATH:", default=1.0, RC=STATUS) ! Entrainment rate options 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.3, RC=STATUS) - call MAPL_GetResource (MAPL, MFPARAMS%ENT0LTS, "EDMF_ENT0LTS:", default=1.5, RC=STATUS) + call MAPL_GetResource (MAPL, MFPARAMS%ENT0LTS, "EDMF_ENT0LTS:", default=1.2, RC=STATUS) ! L0 if ET==1 call MAPL_GetResource (MAPL, MFPARAMS%L0, "EDMF_L0:", default=100., RC=STATUS) ! L0fac if ET==2 From de0867f8c2f167f22694d528a9e27e6dc41c208d Mon Sep 17 00:00:00 2001 From: Nathan Arnold Date: Fri, 9 Feb 2024 22:06:04 -0500 Subject: [PATCH 4/8] Stability fixes in ADG PDF, EDMF parameter changes, qt2 min increased --- .../GEOSmoist_GridComp/Process_Library.F90 | 46 +++++++++---------- .../GEOS_TurbulenceGridComp.F90 | 16 +++++-- .../GEOSturbulence_GridComp/edmf.F90 | 2 +- .../GEOSturbulence_GridComp/shoc.F90 | 10 +++- 4 files changed, 44 insertions(+), 30 deletions(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/Process_Library.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/Process_Library.F90 index 96f1c8e1e..f90ad6718 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/Process_Library.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/Process_Library.F90 @@ -1460,7 +1460,7 @@ subroutine partition_dblgss( fQi, & ! IN sqrtthl = sqrt(thlsec) skew_thl = hl3 / sqrtthl**3 else - sqrtthl = 0.0 + sqrtthl = 1e-3 skew_thl = 0. endif if (qwsec > 0.0) then @@ -1496,7 +1496,7 @@ subroutine partition_dblgss( fQi, & ! IN ! pdf_a = max(0.,min(0.5,mffrc)) ! end if - if (pdf_a>1e-3) then + if (pdf_a>1e-3 .and. pdf_a<0.5) then aterm = pdf_a else aterm = 0.5 @@ -1612,21 +1612,26 @@ subroutine partition_dblgss( fQi, & ! IN IF (qwsec <= rt_tol*rt_tol .or. abs(w1_2-w1_1) <= w_thresh) THEN ! if no active updrafts - if (aterm .lt. 1e-3 .or. aterm.gt.0.499 .or. Skew_qw.eq.0.) then ! if no residual skewness + if (aterm .lt. 1e-3 .or. aterm.gt.0.499 .or. Skew_qw.lt.1e-8) then ! if no residual skewness qw1_1 = total_water qw1_2 = total_water qw2_1 = qwsec qw2_2 = qwsec sqrtqw2_1 = sqrt(qw2_1) sqrtqw2_2 = sqrt(qw2_2) + else ! qw1_1 = total_water ! qw1_2 = total_water ! qw2_1 = qwsec ! qw2_2 = qwsec - wrk1 = min(10.,skew_qw*sqrtqt**3) ! third moment qt +! wrk1 = max(min(10.,skew_qw*sqrtqt**3) ! third moment qt + wrk1 = qt3 qw1_1 = total_water + (wrk1/(2.*aterm-aterm**3/onema**2))**(1./3.) qw1_2 = (total_water -aterm*qw1_1)/onema + if (isnan(qw1_1)) print *,'qw1_1 nan L1636, qt=',total_water,' wrk1=',wrk1,' onema=',onema,' aterm=',aterm + if (isnan(qw1_2)) print *,'qw1_2 nan L1636' + qw2_1 = qwsec - min(0.5*qwsec,max(0.,(aterm/onema)*(qw1_1-total_water)**2)) qw2_2 = qw2_1 sqrtqw2_1 = sqrt(qw2_1) @@ -1637,7 +1642,7 @@ subroutine partition_dblgss( fQi, & ! IN ! corrtest2 = max(-1.0,min(1.0,wqtntrgs/(sqrtw2*sqrtqt))) corrtest2 = max(-1.0,min(1.0,0.5*wqwsec/(sqrtw2*sqrtqt))) - + qw1_1 = - corrtest2 / w1_2 ! A.7 qw1_2 = - corrtest2 / w1_1 ! A.8 @@ -1905,31 +1910,12 @@ subroutine partition_dblgss( fQi, & ! IN rwthl = 0.0 end if - if (isnan(rwqt)) print *,'rwqt>1e3' - if (isnan(rwthl)) print *,'rwthl>1e3' - if (isnan(C1)) print *,'C1>1e3' - if (isnan(C2)) print *,'C2>1e3' - if (isnan(w1_1)) print *,'w1_1>1e3' - if (isnan(w1_2)) print *,'w1_2>1e3' - if (isnan(qw2_2)) print *,'qw2_2>1e3' - if (isnan(qw2_1)) print *,'qw2_1>1e3' - if (isnan(cqt1)) print *,'cqt1>1e3' - if (isnan(cqt2)) print *,'cqt2>1e3' - if (isnan(cthl1)) print *,'cthl1>1e3' - if (isnan(cthl2)) print *,'cthl2>1e3' - if (isnan(thl2_1)) print *,'thl2_1>1e3' - if (isnan(thl2_2)) print *,'thl2_2>1e3' - wql1 = C1*(cqt1*sqrt(w2_1)*sqrt(qw2_1)*rwqt-cthl1*sqrt(w2_1)*sqrt(thl2_1)*rwthl) wql2 = C2*(cqt2*sqrt(w2_2)*sqrt(qw2_2)*rwqt-cthl2*sqrt(w2_2)*sqrt(thl2_2)*rwthl) - if (isnan(wql1)) print *,'wql1>1e3' - if (isnan(wql2)) print *,'wql2>1e3' ! Compute the liquid water flux wqls = aterm * ((w1_1-w_first)*ql1+wql1) + onema * ((w1_2-w_first)*ql2+wql2) wqis = aterm * ((w1_1-w_first)*qi1) + onema * ((w1_2-w_first)*qi2) - if (isnan(wqls)) wqls = 0. - if (isnan(wqis)) wqis = 0. ! diagnostic buoyancy flux. Includes effects from liquid water, ice ! condensate, liquid & ice precipitation @@ -2130,6 +2116,18 @@ subroutine hystpdf( & WQL, & CFn) endif + if (isnan(CFn)) then + print *,'CFn is nan! PL=',PL,' TEP=',TEP,' QVp=',QVp,' QT2=',QT2 + CFn = CFp + end if + if (isnan(WQL)) then + print *,'WQL is nan! PL=',PL,' TEP=',TEP,' QVp=',QVp,' QT2=',QT2 + end if + if (isnan(QCn)) then + print *,'QCn is nan! PL=',PL,' TEP=',TEP,' QVp=',QVp,' QT2=',QT2 + QCn = QCp + end if + IF(USE_BERGERON) THEN DQCALL = QCn - QCp diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSturbulence_GridComp/GEOS_TurbulenceGridComp.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSturbulence_GridComp/GEOS_TurbulenceGridComp.F90 index db4c79ff5..8f241cd5b 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSturbulence_GridComp/GEOS_TurbulenceGridComp.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSturbulence_GridComp/GEOS_TurbulenceGridComp.F90 @@ -3553,7 +3553,7 @@ 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.2, 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%ENTUFAC, "EDMF_ENTUFAC:", default=1.6, RC=STATUS) @@ -3562,11 +3562,11 @@ subroutine REFRESH(IM,JM,LM,RC) ! 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=1.0, RC=STATUS) - call MAPL_GetResource (MAPL, MFPARAMS%AlphaTH, "EDMF_ALPHATH:", default=1.0, RC=STATUS) + call MAPL_GetResource (MAPL, MFPARAMS%AlphaTH, "EDMF_ALPHATH:", default=1.0, RC=STATUS) ! Entrainment rate options 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.3, RC=STATUS) + call MAPL_GetResource (MAPL, MFPARAMS%ENT0, "EDMF_ENT0:", default=0.25, RC=STATUS) call MAPL_GetResource (MAPL, MFPARAMS%ENT0LTS, "EDMF_ENT0LTS:", default=1.2, RC=STATUS) ! L0 if ET==1 call MAPL_GetResource (MAPL, MFPARAMS%L0, "EDMF_L0:", default=100., RC=STATUS) @@ -4409,6 +4409,16 @@ subroutine REFRESH(IM,JM,LM,RC) qt3_tscale, & afrc_tscale, & docanuto ) +! do I = 1,IM +! do J = 1,JM +! if (PHIS(I,J).gt.1e3) then +! do L = 1,LM +! if (Z(I,J,L).lt.10e3) qt2(I,J,L) = max( qt2(I,J,L), (0.05*qt(I,J,L))**2 ) +! end do +! end if +! end do +! end do + end if KPBLMIN = count(PREF < 50000.) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSturbulence_GridComp/edmf.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSturbulence_GridComp/edmf.F90 index f6df882d7..1b60b8cca 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSturbulence_GridComp/edmf.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSturbulence_GridComp/edmf.F90 @@ -320,7 +320,7 @@ SUBROUTINE RUN_EDMF(its,ite, jts,jte, kts,kte, dt, & ! Inputs ! if CLASP disabled: mass-flux if positive surface buoyancy flux and ! TKE at 2nd model level above threshold ! IF ( (wthv > 0.0 .and. TKE3(IH,JH,kte-1)>0.01 .and. MFPARAMS%doclasp==0 .and. phis(IH,JH).lt.2e4) & - IF ( (wthv > 0.0 .and. MFPARAMS%doclasp==0 .and. phis(IH,JH).lt.2e4) & + IF ( (wthv > 0.0 .and. MFPARAMS%doclasp==0 .and. phis(IH,JH).lt.3e4) & .or. (any(mfsrcthl(IH,JH,1:MFPARAMS%NUP) >= -2.0) .and. MFPARAMS%doclasp/=0)) then ! print *,'wthv=',wthv,' wqt=',wqt,' wthl=',wthl,' edmfdepth=',edmfdepth(IH,JH) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSturbulence_GridComp/shoc.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSturbulence_GridComp/shoc.F90 index 44d2d330f..9058ca99c 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSturbulence_GridComp/shoc.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSturbulence_GridComp/shoc.F90 @@ -1232,9 +1232,15 @@ subroutine update_moments( IM, JM, LM, & ! in whl_can(:,:,k) = onemmf*0.5*( whl_edge(:,:,kd) + whl_edge(:,:,ku) ) !+ mfwhl(:,:,kd) + mfwhl(:,:,ku)) ! Restrict QT variance, 1-25% of total water. - qt2(:,:,k) = max(min(qt2(:,:,k),(0.25*QT(:,:,k))**2),(0.01*QT(:,:,k))**2) +! where (ZL(:,:,k).lt.1.6e4) + qt2(:,:,k) = max(min(qt2(:,:,k),(0.25*QT(:,:,k))**2),(0.03*QT(:,:,k))**2) + qt2diag(:,:,k) = max(min(qt2diag(:,:,k),(0.25*QT(:,:,k))**2),(0.03*QT(:,:,k))**2) +! elsewhere +! qt2(:,:,k) = max(min(qt2(:,:,k),(0.25*QT(:,:,k))**2),(0.01*QT(:,:,k))**2) +! qt2diag(:,:,k) = max(min(qt2diag(:,:,k),(0.25*QT(:,:,k))**2),(0.01*QT(:,:,k))**2) +! end where + hl2(:,:,k) = max(min(hl2(:,:,k),HL2MAX),HL2MIN) - qt2diag(:,:,k) = max(min(qt2diag(:,:,k),(0.25*QT(:,:,k))**2),(0.01*QT(:,:,k))**2) hl2diag(:,:,k) = max(min(hl2diag(:,:,k),HL2MAX),HL2MIN) ! Ensure realizibility From 6b6aeb879ada6aaae3513fc440032d2496e0df42 Mon Sep 17 00:00:00 2001 From: Nathan Arnold Date: Fri, 9 Feb 2024 22:17:32 -0500 Subject: [PATCH 5/8] Removed EDMF_DEPTH2 diagnostic --- .../GEOS_TurbulenceGridComp.F90 | 42 +------------------ 1 file changed, 1 insertion(+), 41 deletions(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSturbulence_GridComp/GEOS_TurbulenceGridComp.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSturbulence_GridComp/GEOS_TurbulenceGridComp.F90 index ec650af35..26db2f5ef 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSturbulence_GridComp/GEOS_TurbulenceGridComp.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSturbulence_GridComp/GEOS_TurbulenceGridComp.F90 @@ -1060,16 +1060,6 @@ subroutine SetServices ( GC, RC ) RC=STATUS ) VERIFY_(STATUS) - - call MAPL_AddExportSpec(GC, & - LONG_NAME = 'EDMF_plume_depth_for_entrainment', & - UNITS = 'm', & - SHORT_NAME = 'EDMF_DEPTH2', & - DIMS = MAPL_DimsHorzOnly, & - VLOCATION = MAPL_VLocationNone, & - RC=STATUS ) - VERIFY_(STATUS) - call MAPL_AddExportSpec(GC, & LONG_NAME = 'EDMF_mass_flux', & UNITS = 'kg m s-1', & @@ -2974,7 +2964,7 @@ subroutine REFRESH(IM,JM,LM,RC) SHOCPRNUM,& TKEBUOY,TKESHEAR,TKEDISS,TKETRANS, & SL2, SL3, W2, W3, WQT, WSL, SLQT, W3CANUTO, QT2DIAG,SL2DIAG,SLQTDIAG - real, dimension(:,:), pointer :: LMIX, edmf_depth, edmf_depth2 + real, dimension(:,:), pointer :: LMIX, edmf_depth ! EDMF variables real, dimension(:,:,:), pointer :: edmf_dry_a,edmf_moist_a,edmf_frc, edmf_dry_w,edmf_moist_w, & @@ -3437,8 +3427,6 @@ subroutine REFRESH(IM,JM,LM,RC) VERIFY_(STATUS) call MAPL_GetPointer(EXPORT, edmf_depth, 'EDMF_DEPTH', RC=STATUS) VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT, edmf_depth2, 'EDMF_DEPTH2', RC=STATUS) - VERIFY_(STATUS) call MAPL_GetPointer(EXPORT, mfaw, 'MFAW', RC=STATUS) VERIFY_(STATUS) call MAPL_GetPointer(EXPORT, slflxmf, 'SLFLXMF', RC=STATUS) @@ -3600,34 +3588,6 @@ subroutine REFRESH(IM,JM,LM,RC) call MAPL_TimerOff(MAPL,"---PRELIMS") - ! find thv-based inversion depth - if (associated(EDMF_DEPTH2)) then - tmp3d(:,:,1:LM-1) = (thv(:,:,2:LM)-thv(:,:,1:LM-1)) / (Z(:,:,2:LM)-Z(:,:,1:LM-1)) - do J = 1,JM - do I = 1,IM - temparray = 0. - K = LM - DO WHILE(Z(I,J,K).lt.3e3) - K = K-1 - END DO -! print *,'3km hgt=',Z(I,J,K) - ! at this point, K is index of full level above 3km height - locmax = maxloc(tmp3d(I,J,K:),dim=1) ! returns edge index of max thv gradient - K = K-1+locmax -! print *,'thv max hgt=',Z(I,J,K) - temparray(K:LM-1) = (tmp3d(I,J,K-1:LM-2)-tmp3d(I,J,K:LM-1)) / (Z(I,J,K-1:LM-2)-Z(I,J,K:LM-1)) ! 2nd deriv of THV on full levs -! print *,'2nd deriv=',temparray(K:LM-1) - DO WHILE( K.lt.LM ) - K = K+1 - if (temparray(K).gt.temparray(K+1)) exit - END DO -! print *,'inversion hgt=',Z(I,J,K) - - EDMF_DEPTH2(I,J) = Z(I,J,K) - end do - end do - end if - ! Calculate liquid water potential temperature (THL) and total water (QT) EXF=T/TH THL=TH-(MAPL_ALHL*QL+MAPL_ALHS*QI)/(MAPL_CP*EXF) From fa0fb9636619aed10ec64f4055706934b7df4873 Mon Sep 17 00:00:00 2001 From: Nathan Arnold Date: Mon, 12 Feb 2024 10:47:58 -0500 Subject: [PATCH 6/8] Cleanup and parameter changes --- .../GEOSmoist_GridComp/Process_Library.F90 | 19 ---- .../GEOS_TurbulenceGridComp.F90 | 76 ++++++--------- .../GEOSturbulence_GridComp/shoc.F90 | 97 ++++++++----------- 3 files changed, 72 insertions(+), 120 deletions(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/Process_Library.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/Process_Library.F90 index afe3a0e6c..6a6b82424 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/Process_Library.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/Process_Library.F90 @@ -1669,8 +1669,6 @@ subroutine partition_dblgss( fQi, & ! IN wrk1 = qt3 qw1_1 = total_water + (wrk1/(2.*aterm-aterm**3/onema**2))**(1./3.) qw1_2 = (total_water -aterm*qw1_1)/onema - if (isnan(qw1_1)) print *,'qw1_1 nan L1636, qt=',total_water,' wrk1=',wrk1,' onema=',onema,' aterm=',aterm - if (isnan(qw1_2)) print *,'qw1_2 nan L1636' qw2_1 = qwsec - min(0.5*qwsec,max(0.,(aterm/onema)*(qw1_1-total_water)**2)) qw2_2 = qw2_1 @@ -1914,11 +1912,6 @@ subroutine partition_dblgss( fQi, & ! IN 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) -! diag_qi = diag_qn - diag_ql - -!!! temporary -! if (abs(qc-diag_qn)>0.001) print *,'SHOC: t=',tabs,' s1=',s1,' qn1=',qn1,' qs1=',qs1,' qt1=',qw1_1 ! Update temperature variable based on diagnosed cloud properties @@ -2156,18 +2149,6 @@ subroutine hystpdf( & WQL, & CFn) endif - if (isnan(CFn)) then - print *,'CFn is nan! PL=',PL,' TEP=',TEP,' QVp=',QVp,' QT2=',QT2 - CFn = CFp - end if - if (isnan(WQL)) then - print *,'WQL is nan! PL=',PL,' TEP=',TEP,' QVp=',QVp,' QT2=',QT2 - end if - if (isnan(QCn)) then - print *,'QCn is nan! PL=',PL,' TEP=',TEP,' QVp=',QVp,' QT2=',QT2 - QCn = QCp - end if - IF(USE_BERGERON) THEN DQCALL = QCn - QCp diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSturbulence_GridComp/GEOS_TurbulenceGridComp.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSturbulence_GridComp/GEOS_TurbulenceGridComp.F90 index 26db2f5ef..0687658e9 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSturbulence_GridComp/GEOS_TurbulenceGridComp.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSturbulence_GridComp/GEOS_TurbulenceGridComp.F90 @@ -1070,20 +1070,20 @@ subroutine SetServices ( GC, RC ) VERIFY_(STATUS) call MAPL_AddExportSpec(GC, & - LONG_NAME = 'EDMF_dry_static_energy_source_term', & - UNITS = 'J kg-1 s-1', & - SHORT_NAME = 'SSRCMF', & + LONG_NAME = 'EDMF_dry_static_energy_source_term', & + UNITS = 'J kg-1 s-1', & + SHORT_NAME = 'SSRCMF', & DIMS = MAPL_DimsHorzVert, & - VLOCATION = MAPL_VLocationCenter, & + VLOCATION = MAPL_VLocationCenter, & RC=STATUS ) VERIFY_(STATUS) call MAPL_AddExportSpec(GC, & - LONG_NAME = 'EDMF_specific_humidity_source_term', & + LONG_NAME = 'EDMF_specific_humidity_source_term', & UNITS = 'kg kg-1 s-1', & SHORT_NAME = 'QVSRCMF', & DIMS = MAPL_DimsHorzVert, & - VLOCATION = MAPL_VLocationCenter, & + VLOCATION = MAPL_VLocationCenter, & RC=STATUS ) VERIFY_(STATUS) @@ -1092,7 +1092,7 @@ subroutine SetServices ( GC, RC ) UNITS = 'kg kg-1 s-1', & SHORT_NAME = 'QLSRCMF', & DIMS = MAPL_DimsHorzVert, & - VLOCATION = MAPL_VLocationCenter, & + VLOCATION = MAPL_VLocationCenter, & RC=STATUS ) VERIFY_(STATUS) @@ -3212,10 +3212,10 @@ subroutine REFRESH(IM,JM,LM,RC) 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 = 4.0 , RC=STATUS); VERIFY_(STATUS) - call MAPL_GetResource (MAPL, QT2TUNE, 'QT2TUNE:', DEFAULT = 7.0 , RC=STATUS); VERIFY_(STATUS) + call MAPL_GetResource (MAPL, QT2TUNE, 'QT2TUNE:', DEFAULT = 5.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, QT3_TSCALE, 'QT3_TSCALE:', DEFAULT = 1600.0, RC=STATUS); VERIFY_(STATUS) + call MAPL_GetResource (MAPL, AFRC_TSCALE,'AFRC_TSCALE:',DEFAULT = 1600.0, RC=STATUS); VERIFY_(STATUS) call MAPL_GetResource (MAPL, DOCANUTO, 'DOCANUTO:', DEFAULT = 0, RC=STATUS); VERIFY_(STATUS) ! Get pointers from export state... @@ -3635,8 +3635,6 @@ subroutine REFRESH(IM,JM,LM,RC) ! call MAPL_GetResource (MAPL, EDMF_THERMAL_PLUME, "EDMF_THERMAL_PLUME:", default=0, RC=STATUS) ! call MAPL_GetResource (MAPL, EDMF_TEST, "EDMF_TEST:" , default=0, RC=STATUS) ! call MAPL_GetResource (MAPL, EDMF_DEBUG, "EDMF_DEBUG:", default=0, RC=STATUS) -! call MAPL_GetResource (MAPL, EDMF_WA, "EDMF_WA:", default=1., RC=STATUS) -! call MAPL_GetResource (MAPL, EDMF_WB, "EDMF_WB:", default=1.5, RC=STATUS) ! call MAPL_GetResource (MAPL, EDMF_AU0, "EDMF_AU0:", default=0.14, RC=STATUS) ! call MAPL_GetResource (MAPL, EDMF_CTH1, "EDMF_CTH1:", default=7.2, RC=STATUS) ! call MAPL_GetResource (MAPL, EDMF_CTH2, "EDMF_CTH2:", default=1.1, RC=STATUS) @@ -3734,7 +3732,7 @@ subroutine REFRESH(IM,JM,LM,RC) IF(DOMF /= 0) then - call RUN_EDMF(1, IM, 1, JM, 1, LM, DT, & + call RUN_EDMF(1, IM, 1, JM, 1, LM, DT, & !== Inputs == PHIS, & Z, & @@ -3802,28 +3800,29 @@ subroutine REFRESH(IM,JM,LM,RC) EDMF_PLUMES_THL, & EDMF_PLUMES_QT ) - if (associated(edmf_dry_a)) edmf_dry_a = edmfdrya + !=== Fill Exports === + if (associated(edmf_dry_a)) edmf_dry_a = edmfdrya if (associated(edmf_moist_a)) edmf_moist_a = edmfmoista - if (associated(edmf_buoyf)) edmf_buoyf = buoyf - if (associated(edmf_mfx)) edmf_mfx = edmf_mf - if (associated(mfaw)) mfaw = edmf_mf/rhoe - if (associated(slflxmf)) slflxmf = (aws3-awql3*mapl_alhl-awqi3*mapl_alhs)/mapl_cp - if (associated(qtflxmf)) qtflxmf = awqv3+awql3+awqi3 - if (associated(ssrcmf)) ssrcmf = ssrc - if (associated(qvsrcmf)) qvsrcmf = qvsrc - if (associated(qlsrcmf)) qlsrcmf = qlsrc + if (associated(edmf_buoyf)) edmf_buoyf = buoyf + if (associated(edmf_mfx)) edmf_mfx = edmf_mf + if (associated(mfaw)) mfaw = edmf_mf/rhoe + if (associated(slflxmf)) slflxmf = (aws3-awql3*mapl_alhl-awqi3*mapl_alhs)/mapl_cp + if (associated(qtflxmf)) qtflxmf = awqv3+awql3+awqi3 + if (associated(ssrcmf)) ssrcmf = ssrc + if (associated(qvsrcmf)) qvsrcmf = qvsrc + if (associated(qlsrcmf)) qlsrcmf = qlsrc ! if (associated(edmf_sl2)) edmf_sl2 = mfsl2 ! if (associated(edmf_qt2)) edmf_qt2 = mfqt2 - if (associated(edmf_w2)) edmf_w2 = mfw2 - if (associated(edmf_w3)) edmf_w3 = mfw3 - if (associated(edmf_qt3)) edmf_qt3 = mfqt3 - if (associated(edmf_sl3)) edmf_sl3 = mfsl3 - if (associated(edmf_wqt)) edmf_wqt = mfwqt - if (associated(edmf_slqt)) edmf_slqt = mfslqt - if (associated(edmf_wsl)) edmf_wsl = mfwsl - if (associated(edmf_tke)) edmf_tke = mftke - if (associated(EDMF_FRC)) EDMF_FRC = 0.5*(edmfdrya(:,:,0:LM-1)+edmfdrya(:,:,1:LM) & - + edmfmoista(:,:,0:LM-1)+edmfmoista(:,:,1:LM)) + if (associated(edmf_w2)) edmf_w2 = mfw2 + if (associated(edmf_w3)) edmf_w3 = mfw3 + if (associated(edmf_qt3)) edmf_qt3 = mfqt3 + if (associated(edmf_sl3)) edmf_sl3 = mfsl3 + if (associated(edmf_wqt)) edmf_wqt = mfwqt + if (associated(edmf_slqt)) edmf_slqt = mfslqt + if (associated(edmf_wsl)) edmf_wsl = mfwsl + if (associated(edmf_tke)) edmf_tke = mftke + if (associated(EDMF_FRC)) EDMF_FRC = 0.5*(edmfdrya(:,:,0:LM-1)+edmfdrya(:,:,1:LM) & + + edmfmoista(:,:,0:LM-1)+edmfmoista(:,:,1:LM)) ELSE ! if there is no mass-flux ae3 = 1.0 @@ -3892,7 +3891,7 @@ subroutine REFRESH(IM,JM,LM,RC) call MAPL_TimerOn (MAPL,name="---SHOC" ,RC=STATUS) VERIFY_(STATUS) - call RUN_SHOC( IM, JM, LM, LM+1, DT, & + call RUN_SHOC( IM, JM, LM, LM+1, DT, & !== Inputs == PLO(:,:,1:LM), & ZL0(:,:,0:LM), & @@ -3927,8 +3926,6 @@ subroutine REFRESH(IM,JM,LM,RC) LSHOC2, & LSHOC3, & BRUNTSHOC, & - BRUNTDRY, & - BRUNTEDGE, & RI, & SHOCPRNUM, & !== Tuning params == @@ -4456,15 +4453,6 @@ subroutine REFRESH(IM,JM,LM,RC) qt3_tscale, & afrc_tscale, & docanuto ) -! do I = 1,IM -! do J = 1,JM -! if (PHIS(I,J).gt.1e3) then -! do L = 1,LM -! if (Z(I,J,L).lt.10e3) qt2(I,J,L) = max( qt2(I,J,L), (0.05*qt(I,J,L))**2 ) -! end do -! end if -! end do -! end do end if diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSturbulence_GridComp/shoc.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSturbulence_GridComp/shoc.F90 index 9058ca99c..524791332 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSturbulence_GridComp/shoc.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSturbulence_GridComp/shoc.F90 @@ -59,8 +59,8 @@ subroutine run_shoc( nx, ny, nzm, nz, dtn, & ! in tkesbshear_inv, & ! out smixt_inv, lmix_out, smixt1_inv, & ! out smixt2_inv,smixt3_inv, & ! out -! bruntmst_inv, bruntedg_inv, ri_inv, prnum_inv, & ! out - bruntmst_inv, bruntdry_inv, bruntedg_inv, ri_inv, prnum_inv, & ! out + bruntmst_inv, ri_inv, prnum_inv, & ! out +! bruntmst_inv, bruntdry_inv, bruntedg_inv, ri_inv, prnum_inv, & ! out shocparams ) @@ -118,8 +118,8 @@ subroutine run_shoc( nx, ny, nzm, nz, dtn, & ! in real, dimension(:,:,:), pointer :: smixt2_inv ! length scale, term 2 real, dimension(:,:,:), pointer :: smixt3_inv ! length scale, term 3 real, dimension(:,:,:), pointer :: bruntmst_inv ! moist Brunt vaisala frequency - real, dimension(:,:,:), pointer :: bruntdry_inv ! Brunt vaisala frequency on edges - real, dimension(:,:,:), pointer :: bruntedg_inv ! Brunt vaisala frequency on edges +! real, dimension(:,:,:), pointer :: bruntdry_inv ! Brunt vaisala frequency on edges +! real, dimension(:,:,:), pointer :: bruntedg_inv ! Brunt vaisala frequency on edges real, dimension(:,:,:), pointer :: ri_inv real, dimension(:,:,:), pointer :: prnum_inv @@ -183,8 +183,8 @@ subroutine run_shoc( nx, ny, nzm, nz, dtn, & ! in real isotropy (nx,ny,nzm) ! "Return-to-isotropy" eddy dissipation time scale, s real brunt (nx,ny,nzm) ! Moist Brunt-Vaisalla frequency, s^-1 - real, dimension(nx,ny,nzm) :: total_water, brunt2, def2, thv, brunt_smooth, l_par - real, dimension(nx,ny,nz) :: brunt_edge, brunt_dry + real, dimension(nx,ny,nzm) :: total_water, brunt2, def2, thv, l_par +! real, dimension(nx,ny,nz) :: brunt_edge, brunt_dry real, dimension(nx,ny) :: l_inf, l_mix, zcb, lts!, l_par!, denom, numer, cldarr @@ -196,7 +196,7 @@ subroutine run_shoc( nx, ny, nzm, nz, dtn, & ! in tkeavg, dtqw, dtqi !, l_par integer i,j,k,km1,ku,kd,ka,kb,kinv,strt,fnsh,cnvl - integer, dimension(nx,ny) :: cldbasek +! integer, dimension(nx,ny) :: cldbasek real, parameter :: bruntmin = 1e-7 real, parameter :: vonk = 0.4 @@ -265,7 +265,6 @@ subroutine run_shoc( nx, ny, nzm, nz, dtn, & ! in enddo enddo enddo -! print *,'thv=',thv(i,j,1:20) ! Define vertical grid increments for later use in the vertical differentiation @@ -308,8 +307,8 @@ subroutine run_shoc( nx, ny, nzm, nz, dtn, & ! in if (associated(smixt2_inv)) smixt2_inv(:,:,1:nzm) = smixt2(:,:,nzm:1:-1) if (associated(smixt3_inv)) smixt3_inv(:,:,1:nzm) = smixt3(:,:,nzm:1:-1) - if (associated(bruntedg_inv)) bruntedg_inv(:,:,0:nzm) = brunt_edge(:,:,nz:1:-1) - if (associated(bruntdry_inv)) bruntdry_inv(:,:,0:nzm) = l_par(:,:,nz:1:-1) +! if (associated(bruntedg_inv)) bruntedg_inv(:,:,0:nzm) = brunt_edge(:,:,nz:1:-1) +! if (associated(bruntdry_inv)) bruntdry_inv(:,:,0:nzm) = brunt_dry(:,:,nz:1:-1) if (associated(bruntmst_inv)) bruntmst_inv(:,:,1:nzm) = brunt(:,:,nzm:1:-1) if (associated(prnum_inv)) prnum_inv(:,:,0:nz-1) = prnum(:,:,nz:1:-1) if (associated(ri_inv)) ri_inv(:,:,0:nz-1) = ri(:,:,nz:1:-1) @@ -654,7 +653,7 @@ subroutine eddy_length() ! brunt_edge(:,:,nz) = brunt_edge(:,:,nzm) ! Calculate the measure of PBL depth, Eq. 11 in BK13 (Is this really PBL depth?) - cldbasek(:,:) = 1 +! cldbasek(:,:) = 1 do j=1,ny do i=1,nx @@ -675,22 +674,21 @@ subroutine eddy_length() ! Identify mixed layer top as level where THV exceeds THV(3) + 0.4 K ! Interpolate for final height based on gradient ! Ignore single isolated levels - kk = 4 - do while (thv(i,j,3)+0.4 .gt. thv(i,j,kk) .or. thv(i,j,3)+0.4 .gt. thv(i,j,kk+1)) - kk = kk+1 - end do - dum = (thv(i,j,kk-1)-thv(i,j,kk-2)) +! kk = 4 +! do while (thv(i,j,3)+0.4 .gt. thv(i,j,kk) .or. thv(i,j,3)+0.4 .gt. thv(i,j,kk+1)) +! kk = kk+1 +! end do +! dum = (thv(i,j,kk-1)-thv(i,j,kk-2)) - if (abs(dum) .gt. 1e-3) then - l_mix(i,j) = max(zl(i,j,kk-1)+0.*(thv(i,j,3)+0.4-thv(i,j,kk-1))*(zl(i,j,kk-1)-zl(i,j,kk-2))/dum,100.) - else - l_mix(i,j) = max(zl(i,j,kk-1),100.) - end if +! if (abs(dum) .gt. 1e-3) then +! l_mix(i,j) = max(zl(i,j,kk-1)+0.*(thv(i,j,3)+0.4-thv(i,j,kk-1))*(zl(i,j,kk-1)-zl(i,j,kk-2))/dum,100.) +! else +! l_mix(i,j) = max(zl(i,j,kk-1),100.) +! end if - do while ((zl(i,j,cldbasek(i,j)).lt.300.) .or. (cld_sgs(i,j,cldbasek(i,j)).lt.0.001 .and. cldbasek(i,j).lt.nzm)) - cldbasek(i,j) = cldbasek(i,j) + 1 - end do -! print *,'cldbase=',zl(i,j,cldbasek(i,j)) +! do while ((zl(i,j,cldbasek(i,j)).lt.300.) .or. (cld_sgs(i,j,cldbasek(i,j)).lt.0.001 .and. cldbasek(i,j).lt.nzm)) +! cldbasek(i,j) = cldbasek(i,j) + 1 +! end do kk = 1 do while (zl(i,j,kk) .lt. 3000. .or. kk.eq.nzm) @@ -820,36 +818,27 @@ subroutine eddy_length() end do brunt_edge(:,:,1) = brunt_edge(:,:,2) brunt_edge(:,:,nz) = brunt_edge(:,:,nzm) - - do j=1,ny - do i=1,nx - kk = 1 - l_mix(i,j) = 0. - do while (l_mix(i,j).lt.1e-5 .and. kk.lt.nzm) - l_mix(i,j) = l_mix(i,j) + brunt(i,j,kk)*adzl(i,j,kk) - kk = kk+1 - end do - l_mix(i,j) = zl(i,j,kk) - end do - end do -! print *,'lmix=',l_mix(i,j) + brunt2(:,:,1) = brunt2(:,:,2) + brunt2(:,:,nzm) = brunt2(:,:,nzm-1) + +! do j=1,ny +! do i=1,nx +! kk = 1 +! l_mix(i,j) = 0. +! do while (l_mix(i,j).lt.1e-5 .and. kk.lt.nzm) +! l_mix(i,j) = l_mix(i,j) + brunt(i,j,kk)*adzl(i,j,kk) +! kk = kk+1 +! end do +! l_mix(i,j) = zl(i,j,kk) +! end do +! end do ! brunt_dry = max( bruntmin, brunt_dry ) - brunt_smooth = brunt - brunt_smooth(:,:,nzm) = brunt(:,:,nzm-1) - brunt_smooth(:,:,1) = brunt(:,:,1) - brunt_smooth = max( bruntmin, brunt_smooth ) - - - do k=1,nzm do j=1,ny do i=1,nx - -! if (k.eq.1) print *,'l_mix=',l_mix(i,j) - tkes = sqrt(tke(i,j,k)) ! Calculate turbulent length scale in the boundary layer. @@ -874,7 +863,6 @@ subroutine eddy_length() l_par(i,j,k) = l_par(i,j,k) - zl(i,j,kk) + max(0.,(thv(i,j,kk)-wrk)* & (zl(i,j,kk)-zl(i,j,kk-1))/(thv(i,j,kk)-thv(i,j,kk-1))) l_par(i,j,k) = max(min(l_par(i,j,k),1500.),25.) -! if (zl(i,j,k).lt.600.) print *,'z=',zl(i,j,k),' thv=',thv(i,j,k),' thv pert=',max(0.001,3.*wthv_sec(i,j,k)) / max(0.05,sqrt(0.667)*tkes),' l_par=',l_par(i,j,k) !---------------------------------- @@ -904,7 +892,7 @@ subroutine eddy_length() smixt2(i,j,k) = sqrt(l_par(i,j,k)*400.*tkes)*(shocparams%LENFAC2) ! Stability length scale - smixt3(i,j,k) = max(0.1,tkes)*shocparams%LENFAC3/(sqrt(brunt_smooth(i,j,k))) + smixt3(i,j,k) = max(0.1,tkes)*shocparams%LENFAC3/(sqrt(brunt2(i,j,k))) !=== Combine component length scales === @@ -924,7 +912,7 @@ subroutine eddy_length() wrk2 = 1.0/(400.*tkes) wrk3 = sqrt(brunt2(i,j,k))/(0.7*tkes) wrk1 = 1.0/(wrk2+wrk3) - smixt(i,j,k) = 3.3*shocparams%LENFAC1*(wrk1 + (vonk*zl(i,j,k)-wrk1)*exp(-zl(i,j,k)/(0.1*l_mix(i,j)))) + smixt(i,j,k) = 3.3*shocparams%LENFAC1*(wrk1 + (vonk*zl(i,j,k)-wrk1)*exp(-zl(i,j,k)/(0.1*l_par(i,j)))) smixt1(i,j,k) = 3.3*shocparams%LENFAC1/wrk2 smixt2(i,j,k) = 3.3*shocparams%LENFAC1/wrk3 smixt3(i,j,k) = 3.3*shocparams%LENFAC1*vonk*zl(i,j,k) @@ -1231,14 +1219,9 @@ subroutine update_moments( IM, JM, LM, & ! in whl(:,:,k) = onemmf*0.5*( whl_edge(:,:,kd) + whl_edge(:,:,ku) ) !+ MFWHL(:,:,k) whl_can(:,:,k) = onemmf*0.5*( whl_edge(:,:,kd) + whl_edge(:,:,ku) ) !+ mfwhl(:,:,kd) + mfwhl(:,:,ku)) - ! Restrict QT variance, 1-25% of total water. -! where (ZL(:,:,k).lt.1.6e4) + ! Restrict QT variance, 3-25% of total water. qt2(:,:,k) = max(min(qt2(:,:,k),(0.25*QT(:,:,k))**2),(0.03*QT(:,:,k))**2) qt2diag(:,:,k) = max(min(qt2diag(:,:,k),(0.25*QT(:,:,k))**2),(0.03*QT(:,:,k))**2) -! elsewhere -! qt2(:,:,k) = max(min(qt2(:,:,k),(0.25*QT(:,:,k))**2),(0.01*QT(:,:,k))**2) -! qt2diag(:,:,k) = max(min(qt2diag(:,:,k),(0.25*QT(:,:,k))**2),(0.01*QT(:,:,k))**2) -! end where hl2(:,:,k) = max(min(hl2(:,:,k),HL2MAX),HL2MIN) hl2diag(:,:,k) = max(min(hl2diag(:,:,k),HL2MAX),HL2MIN) From 5fb0f1e6e38978fe2e3ac8c1fe0418ec01225d0a Mon Sep 17 00:00:00 2001 From: Nathan Arnold Date: Mon, 12 Feb 2024 17:38:48 -0500 Subject: [PATCH 7/8] Fix for compile --- .../GEOSphysics_GridComp/GEOSturbulence_GridComp/shoc.F90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSturbulence_GridComp/shoc.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSturbulence_GridComp/shoc.F90 index 524791332..92cb6eac1 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSturbulence_GridComp/shoc.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSturbulence_GridComp/shoc.F90 @@ -184,7 +184,7 @@ subroutine run_shoc( nx, ny, nzm, nz, dtn, & ! in real brunt (nx,ny,nzm) ! Moist Brunt-Vaisalla frequency, s^-1 real, dimension(nx,ny,nzm) :: total_water, brunt2, def2, thv, l_par -! real, dimension(nx,ny,nz) :: brunt_edge, brunt_dry + real, dimension(nx,ny,nz) :: brunt_edge !, brunt_dry real, dimension(nx,ny) :: l_inf, l_mix, zcb, lts!, l_par!, denom, numer, cldarr @@ -628,7 +628,7 @@ subroutine eddy_length() ! brunt_edge(i,j,k) = (2.*ggr/(thv(i,j,k)+thv(i,j,kb)))*(thv(i,j,k)-thv(i,j,kb))/adzi(i,j,k) !g/thv/dz *(thv-thv) - brunt_dry(i,j,k) = (thv(i,j,kc)-thv(i,j,kb))*betdz +! brunt_dry(i,j,k) = (thv(i,j,kc)-thv(i,j,kb))*betdz ! Reinitialize the mixing length related arrays to zero smixt(i,j,k) = 1.0 ! shoc_mod module variable smixt @@ -912,7 +912,7 @@ subroutine eddy_length() wrk2 = 1.0/(400.*tkes) wrk3 = sqrt(brunt2(i,j,k))/(0.7*tkes) wrk1 = 1.0/(wrk2+wrk3) - smixt(i,j,k) = 3.3*shocparams%LENFAC1*(wrk1 + (vonk*zl(i,j,k)-wrk1)*exp(-zl(i,j,k)/(0.1*l_par(i,j)))) + smixt(i,j,k) = 3.3*shocparams%LENFAC1*(wrk1 + (vonk*zl(i,j,k)-wrk1)*exp(-zl(i,j,k)/(0.1*zpbl(i,j)))) smixt1(i,j,k) = 3.3*shocparams%LENFAC1/wrk2 smixt2(i,j,k) = 3.3*shocparams%LENFAC1/wrk3 smixt3(i,j,k) = 3.3*shocparams%LENFAC1*vonk*zl(i,j,k) From 768a88c78f753f0c398eed21dc4d01240adfe4ce Mon Sep 17 00:00:00 2001 From: Nathan Arnold Date: Tue, 13 Feb 2024 16:52:20 -0500 Subject: [PATCH 8/8] Removed longname change in MoistGC to make 0-diff --- .../GEOSmoist_GridComp/GEOS_MoistGridComp.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/GEOS_MoistGridComp.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/GEOS_MoistGridComp.F90 index b6a0d5005..46c83c19d 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/GEOS_MoistGridComp.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/GEOS_MoistGridComp.F90 @@ -666,7 +666,7 @@ subroutine SetServices ( GC, RC ) call MAPL_AddImportSpec(GC, & SHORT_NAME = 'RADSW', & - LONG_NAME = 'air_temperature_tendency_due_to_shortwave', & + LONG_NAME = 'air_temperature_tendency_due_to_longwave', & UNITS = 'K s-1', & DIMS = MAPL_DimsHorzVert, & VLOCATION = MAPL_VLocationCenter, &