diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/Process_Library.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/Process_Library.F90 index 7f538648f..e23bb48b4 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/Process_Library.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/Process_Library.F90 @@ -1537,21 +1537,21 @@ subroutine partition_dblgss( fQi, & ! IN ! Compute square roots of some variables so we don't have to do it again - if (w_sec > 0.0) then + if (w_sec > w_thresh*w_thresh) then sqrtw2 = sqrt(w_sec) Skew_w = w3var / (sqrtw2*sqrtw2*sqrtw2) else sqrtw2 = w_thresh Skew_w = 0. endif - if (thlsec > 0.0) then + if (thlsec > 1e-6) then sqrtthl = sqrt(thlsec) skew_thl = hl3 / sqrtthl**3 else sqrtthl = 1e-3 skew_thl = 0. endif - if (qwsec > 0.0) then + if (qwsec > 1e-8*total_water*total_water) then sqrtqt = sqrt(qwsec) skew_qw = qt3/sqrtqt**3 else @@ -1700,7 +1700,7 @@ 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.lt.1e-8) then ! if no residual skewness + if (aterm .lt. 1e-3 .or. aterm.gt.0.499 .or. abs(Skew_qw).lt.1e-8) then ! if no residual skewness qw1_1 = total_water qw1_2 = total_water qw2_1 = qwsec @@ -1851,7 +1851,6 @@ subroutine partition_dblgss( fQi, & ! IN ! this is qs evaluated at Tl qs1 = om1 * (0.622*esval1_1/max(esval1_1,pval-0.378*esval1_1)) & + (1.-om1) * (0.622*esval2_1/max(esval2_1,pval-0.378*esval2_1)) -! qs1 = GEOS_QSAT( Tl1_1, pval*0.01 ) beta1 = (lstarn1*lstarn1*onebrvcp) / (Tl1_1*Tl1_1) @@ -1862,18 +1861,10 @@ subroutine partition_dblgss( fQi, & ! 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 esval1_2 = MAPL_EQsat(Tl1_2) esval2_2 = MAPL_EQsat(Tl1_2,OverIce=.TRUE.) om2 = max(0.,min(1.,1.-fQi)) !max(0.,min(1.,a_bg*(Tl1_2-tbgmin))) lstarn2 = lcond + (1.-om2)*lfus -! 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)) @@ -1945,7 +1936,7 @@ subroutine partition_dblgss( fQi, & ! IN ! finally, compute the SGS cloud fraction - diag_frac = aterm*C1 + onema*C2 + cld_sgs = aterm*C1 + onema*C2 ! om1 = max(0.,min(1.,(Tl1_1-tbgmin)*a_bg)) ! om2 = max(0.,min(1.,(Tl1_2-tbgmin)*a_bg)) @@ -1970,21 +1961,13 @@ subroutine partition_dblgss( fQi, & ! IN ! + tkesbdiss(i,j,k) * (dtn/cp) ! tke dissipative heating ! Update moisture fields - - -! qc = diag_ql + diag_qi -! qi = diag_qi -! qwv = total_water - diag_qn - cld_sgs = diag_frac - - if (sqrtqt>0.0 .AND. sqrtw2>0.0) then - rwqt = (1.-0.5)*wqwsec/(sqrtqt*sqrtw2) -! rwqt = (wqwsec)/(sqrtqt*sqrtw2) + if (sqrtqt>1e-4*total_water .AND. sqrtw2>w_thresh) then + rwqt = 0.5*wqwsec/(sqrtqt*sqrtw2) ! rwqt = max(-1.,min(1.,pdf_rwqt)) else rwqt = 0.0 end if - if (sqrtthl>0.0 .AND. sqrtw2>0.0) then + if (sqrtthl>1e-3 .AND. sqrtw2>w_thresh) then rwthl = wthlsec/(sqrtthl*sqrtw2) ! rwthl = max(-1.,min(1.,pdf_rwth)) else @@ -2007,8 +1990,7 @@ subroutine partition_dblgss( fQi, & ! IN wthv_sec = wthlsec + wrk*wqwsec & + (fac_cond-bastoeps)*wqls & + (fac_sub-bastoeps) *wqis - -! + ((lstarn1/cp)-thv(i,j,k))*0.5*(wqp_sec(i,j,kd)+wqp_sec(i,j,ku)) +! + ((lstarn1/cp)-thv(i,j,k))*0.5*(wqp_sec(i,j,kd)+wqp_sec(i,j,ku)) end subroutine partition_dblgss diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSturbulence_GridComp/GEOS_TurbulenceGridComp.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSturbulence_GridComp/GEOS_TurbulenceGridComp.F90 index efb0a9d5c..daedadf94 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSturbulence_GridComp/GEOS_TurbulenceGridComp.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSturbulence_GridComp/GEOS_TurbulenceGridComp.F90 @@ -1487,6 +1487,34 @@ subroutine SetServices ( GC, RC ) RC=STATUS ) VERIFY_(STATUS) + call MAPL_AddExportSpec(GC, & + LONG_NAME = 'trade_inversion_base_pressure', & + UNITS = 'Pa', & + SHORT_NAME = 'TRINVBS', & + DIMS = MAPL_DimsHorzOnly, & + VLOCATION = MAPL_VLocationNone, & + RC=STATUS ) + VERIFY_(STATUS) + + call MAPL_AddExportSpec(GC, & + LONG_NAME = 'trade_inversion_temperature_jump', & + UNITS = 'K', & + SHORT_NAME = 'TRINVDELT', & + DIMS = MAPL_DimsHorzOnly, & + VLOCATION = MAPL_VLocationNone, & + RC=STATUS ) + VERIFY_(STATUS) + + call MAPL_AddExportSpec(GC, & + LONG_NAME = 'trade_inversion_frequency', & + UNITS = '1', & + SHORT_NAME = 'TRINVFRQ', & + DIMS = MAPL_DimsHorzOnly, & + VLOCATION = MAPL_VLocationNone, & + RC=STATUS ) + VERIFY_(STATUS) + + call MAPL_AddExportSpec(GC, & LONG_NAME = 'Buoyancy_jump_across_inversion', & UNITS = 'm s-2', & @@ -1971,22 +1999,6 @@ 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', & @@ -2178,34 +2190,13 @@ subroutine SetServices ( GC, RC ) RC=STATUS ) VERIFY_(STATUS) - - call MAPL_AddInternalSpec(GC, & - LONG_NAME = 'sensitivity_of_tendency_to_surface_value_for_s', & - SHORT_NAME = 'DKSS', & - UNITS = 's-1', & - DIMS = MAPL_DimsHorzVert, & - VLOCATION = MAPL_VLocationCenter, & - RESTART = MAPL_RestartSkip, & - RC=STATUS ) - VERIFY_(STATUS) - - call MAPL_AddInternalSpec(GC, & - LONG_NAME = 'sensitivity_of_tendency_to_surface_value_for_q', & - SHORT_NAME = 'DKQQ', & - UNITS = 's-1', & - DIMS = MAPL_DimsHorzVert, & - VLOCATION = MAPL_VLocationCenter, & - RESTART = MAPL_RestartSkip, & - RC=STATUS ) - VERIFY_(STATUS) - call MAPL_AddInternalSpec(GC, & - LONG_NAME = 'sensitivity_of_tendency_to_surface_value_for_u', & - SHORT_NAME = 'DKUU', & - UNITS = 's-1', & + LONG_NAME = 'SHOC_TKE_dissipation_rate', & + SHORT_NAME = 'TKEDISS', & + UNITS = 'J kg-1 s-1', & DIMS = MAPL_DimsHorzVert, & VLOCATION = MAPL_VLocationCenter, & - RESTART = MAPL_RestartSkip, & + RESTART = MAPL_RestartSkip, & RC=STATUS ) VERIFY_(STATUS) @@ -2616,7 +2607,6 @@ subroutine RUN1 ( GC, IMPORT, EXPORT, CLOCK, RC ) real, dimension(:,:,:), pointer :: AKSS, BKSS, CKSS, YS real, dimension(:,:,:), pointer :: AKQQ, BKQQ, CKQQ, YQV,YQL,YQI real, dimension(:,:,:), pointer :: AKUU, BKUU, CKUU, YU,YV - real, dimension(:,:,:), pointer :: DKSS, DKQQ, DKUU ! SHOC-related variables integer :: DO_SHOC, SCM_SL @@ -2792,12 +2782,6 @@ subroutine RUN1 ( GC, IMPORT, EXPORT, CLOCK, RC ) ! edmf variables ! - call MAPL_GetPointer(INTERNAL, DKSS, 'DKSS', RC=STATUS) - VERIFY_(STATUS) - call MAPL_GetPointer(INTERNAL, DKQQ, 'DKQQ', RC=STATUS) - VERIFY_(STATUS) - call MAPL_GetPointer(INTERNAL, DKUU, 'DKUU', RC=STATUS) - VERIFY_(STATUS) ! a,b,c and rhs for s call MAPL_GetPointer(INTERNAL, AKSS, 'AKSS', RC=STATUS) VERIFY_(STATUS) @@ -2965,11 +2949,12 @@ subroutine REFRESH(IM,JM,LM,RC) real, dimension(IM,JM,1:LM-1) :: TVE, RDZ real, dimension(IM,JM,LM) :: THV, TV, Z, DMI, PLO, QL, QI, QA, TSM, USM, VSM real, dimension(IM,JM,0:LM) :: ZL0 + real, dimension(IM,JM) :: drycblh integer, dimension(IM,JM) :: SMTH_LEV ! real, dimension(:,:,:), pointer :: MFQTSRC, MFTHSRC, MFW, MFAREA real, dimension(:,:,:), pointer :: EKH, EKM, KHLS, KMLS, KHRAD, KHSFC - real, dimension(:,: ), pointer :: BSTAR, USTAR, PPBL, WERAD, WESFC,VSCRAD,KERAD,DBUOY,ZSML,ZCLD,ZRADML,FRLAND + real, dimension(:,: ), pointer :: BSTAR, USTAR, PPBL, WERAD, WESFC,VSCRAD,KERAD,DBUOY,ZSML,ZCLD,ZRADML,FRLAND,TRINVBS,TRINVFRQ,TRINVDELT real, dimension(:,: ), pointer :: TCZPBL => null() real, dimension(:,: ), pointer :: ZPBL2 => null() real, dimension(:,: ), pointer :: ZPBL10P => null() @@ -2991,10 +2976,10 @@ subroutine REFRESH(IM,JM,LM,RC) real, dimension(:,:,:), pointer :: AKQODT, CKQODT real, dimension(:,:,:), pointer :: AKVODT, CKVODT - real, dimension(:,:,:), pointer :: LSHOC,BRUNTSHOC,BRUNTDRY, BRUNTEDGE,ISOTROPY, & + real, dimension(:,:,:), pointer :: LSHOC,BRUNTSHOC,ISOTROPY, & LSHOC1,LSHOC2,LSHOC3, & SHOCPRNUM,& - TKEBUOY,TKESHEAR,TKEDISS,TKETRANS, & + TKEBUOY,TKESHEAR,TKEDISS,TKEDISSx,TKETRANS, & SL2, SL3, W2, W3, WQT, WSL, SLQT, W3CANUTO, QT2DIAG,SL2DIAG,SLQTDIAG real, dimension(:,:), pointer :: LMIX, edmf_depth @@ -3101,7 +3086,7 @@ subroutine REFRESH(IM,JM,LM,RC) real, dimension( IM, JM, LM ) :: QPL,QPI integer :: DO_SHOC, DOPROGQT2, DOCANUTO real :: SL2TUNE, QT2TUNE, SLQT2TUNE, & - QT3_TSCALE, AFRC_TSCALE + SKEW_TGEN, SKEW_TDIS real :: PDFSHAPE real :: lambdadiss @@ -3249,25 +3234,25 @@ 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.08, RC=STATUS); VERIFY_(STATUS) + call MAPL_GetResource (MAPL, SHOCPARAMS%LAMBDA, trim(COMP_NAME)//"_SHC_LAMBDA:", default=0.25, 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) 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=10.0, RC=STATUS); VERIFY_(STATUS) + call MAPL_GetResource (MAPL, SHOCPARAMS%LENFAC1, trim(COMP_NAME)//"_SHC_LENFAC1:", default=8.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=3.0, RC=STATUS); VERIFY_(STATUS) + call MAPL_GetResource (MAPL, SHOCPARAMS%LENFAC3, trim(COMP_NAME)//"_SHC_LENFAC3:", default=1.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 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 = 5.0 , RC=STATUS); VERIFY_(STATUS) + call MAPL_GetResource (MAPL, QT2TUNE, 'QT2TUNE:', DEFAULT = 9.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 = 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, SKEW_TDIS, 'SKEW_TDIS:', DEFAULT = 1600.0, RC=STATUS); VERIFY_(STATUS) + call MAPL_GetResource (MAPL, SKEW_TGEN, 'SKEW_TGEN:', DEFAULT = 900.0, RC=STATUS); VERIFY_(STATUS) call MAPL_GetResource (MAPL, DOCANUTO, 'DOCANUTO:', DEFAULT = 0, RC=STATUS); VERIFY_(STATUS) ! Get pointers from export state... @@ -3335,6 +3320,12 @@ subroutine REFRESH(IM,JM,LM,RC) VERIFY_(STATUS) call MAPL_GetPointer(EXPORT, DBUOY, 'DBUOY', RC=STATUS) VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT, TRINVBS, 'TRINVBS', RC=STATUS) + VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT, TRINVDELT, 'TRINVDELT', RC=STATUS) + VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT, TRINVFRQ, 'TRINVFRQ', RC=STATUS) + VERIFY_(STATUS) call MAPL_GetPointer(EXPORT, VSCRAD, 'VSCRAD', RC=STATUS) VERIFY_(STATUS) call MAPL_GetPointer(EXPORT, VSCsfc, 'VSCSFC', RC=STATUS) @@ -3487,7 +3478,7 @@ subroutine REFRESH(IM,JM,LM,RC) VERIFY_(STATUS) !========== SHOC =========== - call MAPL_GetPointer(EXPORT, TKEDISS, 'TKEDISS', RC=STATUS) + call MAPL_GetPointer(EXPORT, TKEDISSx, 'TKEDISS', RC=STATUS) VERIFY_(STATUS) call MAPL_GetPointer(EXPORT, TKEBUOY, 'TKEBUOY', RC=STATUS) VERIFY_(STATUS) @@ -3509,13 +3500,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) + call MAPL_GetPointer(INTERNAL, TKEDISS, 'TKEDISS' , RC=STATUS); VERIFY_(STATUS) ! Initialize some arrays LWCRT = RADLW - RADLWC @@ -3644,21 +3632,21 @@ 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) + call MAPL_GetResource (MAPL, MFPARAMS%ENTUFAC, "EDMF_ENTUFAC:", default=2.0, 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=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.25, RC=STATUS) - call MAPL_GetResource (MAPL, MFPARAMS%ENT0LTS, "EDMF_ENT0LTS:", default=1.2, RC=STATUS) + call MAPL_GetResource (MAPL, MFPARAMS%ENT0, "EDMF_ENT0:", default=0.4, RC=STATUS) + call MAPL_GetResource (MAPL, MFPARAMS%ENT0LTS, "EDMF_ENT0LTS:", default=0.8, RC=STATUS) ! L0 if ET==1 call MAPL_GetResource (MAPL, MFPARAMS%L0, "EDMF_L0:", default=100., RC=STATUS) ! L0fac if ET==2 @@ -3685,6 +3673,9 @@ subroutine REFRESH(IM,JM,LM,RC) ! call MAPL_GetResource (MAPL, EDMF_RHO_QB, "EDMF_RHO_QB:", default=0.5, RC=STATUS) ! call MAPL_GetResource (MAPL, C_KH_MF, "C_KH_MF:", default=0., RC=STATUS) ! call MAPL_GetResource (MAPL, NumUpQ, "EDMF_NumUpQ:", default=1, RC=STATUS) + call MAPL_GetResource (MAPL, MFPARAMS%TREFF, "EDMF_TREFF:", default=100., RC=STATUS) + else + call MAPL_GetResource (MAPL, MFPARAMS%TREFF, "EDMF_TREFF:", default=0., RC=STATUS) end if call MAPL_GetResource(MAPL, SCM_SL, 'SCM_SL:', DEFAULT=0 ) @@ -3740,14 +3731,21 @@ subroutine REFRESH(IM,JM,LM,RC) cu => cu_scm ct => ct_scm cq => ct_scm - ustar_scm = 0.3 ! sqrt(CU*UU/RHOS) +! ustar_scm = sqrt(CU*sqrt(U(:,:,LM)**2+V(:,:,LM)**2)/RHOE(:,:,LM)) + ustar_scm = 0.25 !sqrt(CU*U(:,:,LM)/RHOE(:,:,LM)) +! print *,'ustar=',ustar_scm,' cu=',cu + bstar_scm = 0.002 ! 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)) +! bstar_scm = (MAPL_GRAV/(RHOE(:,:,LM)*ustar_scm)) * & +! (SH/THV(:,:,LM) + MAPL_VIREPS*EVAP) ustar => ustar_scm sh => sh_scm evap => evap_scm + print *,'bstar=',bstar_scm,' ustar=',ustar_scm + call MAPL_TimerOff(MAPL,"---SURFACE") end if @@ -3774,6 +3772,7 @@ subroutine REFRESH(IM,JM,LM,RC) qvsrc = 0.0 qlsrc = 0.0 + IF(DOMF /= 0) then call RUN_EDMF(1, IM, 1, JM, 1, LM, DT, & @@ -3789,14 +3788,13 @@ subroutine REFRESH(IM,JM,LM,RC) T, & THL, & THV, & - QT, & Q, & - QL, & - QI, & + QLTOT, & + QITOT, & SH, & EVAP, & - FRLAND, & - ZPBL, & + FRLAND, & + ZPBL, & ! MFTHSRC, MFQTSRC, MFW, MFAREA, & ! CLASP inputs !== Outputs for trisolver == ae3, & @@ -3816,8 +3814,6 @@ subroutine REFRESH(IM,JM,LM,RC) mfqt3, & mfsl3, & mfwqt, & -! mfqt2, & -! mfsl2, & mfslqt, & mfwsl, & !== Outputs for SHOC == @@ -3849,14 +3845,12 @@ subroutine REFRESH(IM,JM,LM,RC) 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(mfaw)) mfaw = aw3 !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 @@ -3867,6 +3861,15 @@ subroutine REFRESH(IM,JM,LM,RC) 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)) + do i = 1,IM + do j = 1,JM + k = LM + do while (edmfdrya(i,j,k).gt.0.01 .and. edmfmoista(i,j,k).lt.1e-4 .and. k.gt.1) + k = k-1 + end do + drycblh(i,j) = ZL0(i,j,k) + end do + end do ELSE ! if there is no mass-flux ae3 = 1.0 @@ -3901,8 +3904,6 @@ subroutine REFRESH(IM,JM,LM,RC) 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 -! 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 @@ -3912,7 +3913,8 @@ subroutine REFRESH(IM,JM,LM,RC) if (associated(edmf_wsl)) edmf_wsl = mfwsl if (associated(edmf_tke)) edmf_tke = mftke if (associated(EDMF_FRC)) EDMF_FRC = 0. - + + drycblh = 0. ENDIF call MAPL_TimerOff(MAPL,"---MASSFLUX") @@ -3953,15 +3955,15 @@ subroutine REFRESH(IM,JM,LM,RC) WTHV2(:,:,1:LM), & BUOYF(:,:,1:LM), & MFTKE(:,:,0:LM), & - ZPBL(:,:), & + DRYCBLH(:,:), & !== Input-Outputs == TKESHOC(:,:,1:LM), & TKH(:,:,1:LM), & !== Outputs == KM(:,:,1:LM), & ISOTROPY(:,:,1:LM), & - !== Diagnostics == ! not used elsewhere TKEDISS, & + !== Diagnostics == ! not used elsewhere TKEBUOY, & TKESHEAR, & LSHOC, & @@ -3975,6 +3977,7 @@ subroutine REFRESH(IM,JM,LM,RC) !== Tuning params == SHOCPARAMS ) + if (associated(TKEDISSx)) TKEDISSx = TKEDISS KH(:,:,1:LM) = TKH(:,:,1:LM) call MAPL_TimerOff (MAPL,name="---SHOC" ,RC=STATUS) @@ -4450,7 +4453,7 @@ subroutine REFRESH(IM,JM,LM,RC) end if ! TKE ! Update the higher order moments required for the ADG PDF - if ( (PDFSHAPE.eq.5) .AND. (DO_SHOC /= 0) ) then + if ( (PDFSHAPE.ge.5) ) then SL = T + (MAPL_GRAV*Z - MAPL_ALHL*QLTOT - MAPL_ALHS*QITOT)/MAPL_CP call update_moments(IM, JM, LM, DT, & SH, & ! in @@ -4493,8 +4496,8 @@ subroutine REFRESH(IM,JM,LM,RC) sl2tune, & qt2tune, & slqt2tune, & - qt3_tscale, & - afrc_tscale, & + skew_tgen, & + skew_tdis, & docanuto ) end if @@ -4744,6 +4747,47 @@ subroutine REFRESH(IM,JM,LM,RC) end if ! SBITOP, SBIFRQ + ! Trade inversion base height + if (associated(TRINVBS)) then + TRINVBS = MAPL_UNDEF + TRINVDELT = MAPL_UNDEF + TRINVFRQ = 0. + do I = 1,IM + do J = 1,JM + K = LM + + do while (PLO(I,J,K).gt.95000.) + K = K-1 + end do + do L = K,1,-1 ! K is first level above 950mb + if (PLO(I,J,L).lt.60000.) exit + + if (T(I,J,L-1).ge.T(I,J,L)) then ! if next level is warmer... + LTOP = L ! L is index of minimum T so far + do while (T(I,J,LTOP).ge.T(I,J,L)) ! find depth of warm layer + LTOP = LTOP-1 + end do + LTOP = LTOP+1 ! LTOP is index of highest level inside warm layer + + if ( MAXVAL(T(I,J,LTOP:L))-T(I,J,L).ge.0.5 .or. & + (MAXVAL(T(I,J,LTOP:L))-T(I,J,L).gt.0.01 .and. PLO(I,J,L)-PLO(I,J,LTOP)>2500.) ) then + + ! only save if DELTA-T exceeds any previous inversion + if ( TRINVFRQ(I,J).eq.0. .or. & + (TRINVFRQ(I,J).ne.0. .and. MAXVAL(T(I,J,LTOP:L))-T(I,J,L).gt.TRINVDELT(I,J)) ) then + TRINVBS(I,J) = PLO(I,J,L) + TRINVDELT(I,J) = MAXVAL(T(I,J,LTOP:L))-T(I,J,L) + TRINVFRQ(I,J) = 1. + end if + + end if + end if ! next level warmer + + end do ! L vert loop + + end do + end do + end if SELECT CASE(PBLHT_OPTION) @@ -4815,23 +4859,25 @@ subroutine REFRESH(IM,JM,LM,RC) PPBL = MAX(PPBL,PLO(:,:,KPBLMIN)) end if + RHOAW3=RHOE*AW3 ! mass flux (positive) + ! Second difference coefficients for scalars; RDZ is RHO/DZ, DMI is (G DT)/DP ! --------------------------------------------------------------------------- - CKS(:,:,1:LM-1) = -KH(:,:,1:LM-1) * RDZ(:,:,1:LM-1) + CKS(:,:,1:LM-1) = -(KH(:,:,1:LM-1)+MFPARAMS%TREFF*RHOAW3(:,:,1:LM-1)) * RDZ(:,:,1:LM-1) AKS(:,:,1 ) = 0.0 AKS(:,:,2:LM ) = CKS(:,:,1:LM-1) * DMI(:,:,2:LM ) CKS(:,:,1:LM-1) = CKS(:,:,1:LM-1) * DMI(:,:,1:LM-1) CKS(:,:, LM ) = -CT * DMI(:,:, LM ) - ! Fill KH at level LM+1 with CT * RDZ for diagnostic output + ! Fill KH at level LM+1 with CT / RDZ for diagnostic output ! --------------------------------------------------------- KH(:,:,LM) = CT * Z(:,:,LM)*((MAPL_RGAS * TV(:,:,LM))/PLE(:,:,LM)) TKH = KH - ! Water vapor can differ at the surface - !-------------------------------------- + ! Water vapor exchange coefficient can differ at the surface + !----------------------------------------------------------- AKQ = AKS CKQ = CKS @@ -4850,10 +4896,10 @@ subroutine REFRESH(IM,JM,LM,RC) CKV(:,:, LM ) = - CU * DMI(:,:, LM ) EKV(:,:, LM ) = MAPL_GRAV * CU - ! Fill KM at level LM with CU * RDZ for diagnostic output + ! Fill KM at level LM with CU / RDZ for diagnostic output ! ------------------------------------------------------- - KM(:,:,LM) = CU * (PLE(:,:,LM)/(MAPL_RGAS * TV(:,:,LM))) / Z(:,:,LM) + KM(:,:,LM) = CU * Z(:,:,LM)*(MAPL_RGAS * TV(:,:,LM)) / PLE(:,:,LM) !Z/rho ! Setup the tridiagonal matrix ! ---------------------------- @@ -5063,10 +5109,6 @@ subroutine REFRESH(IM,JM,LM,RC) call VTRISOLVESURF(BKQ,CKQ,DKQ) call VTRISOLVESURF(BKV,CKV,DKV) - call VTRISOLVESURF(BKSS,CKSS,DKSS) - call VTRISOLVESURF(BKQQ,CKQQ,DKQQ) - call VTRISOLVESURF(BKUU,CKUU,DKUU) - call MAPL_TimerOff(MAPL,"---DECOMP") if(ALLOC_TCZPBL) deallocate(TCZPBL) @@ -5124,6 +5166,7 @@ subroutine DIFFUSE(IM,JM,LM,RC) integer :: KM, K,L logical :: FRIENDLY logical :: WEIGHTED + logical :: OPT ! selects algorithm in VTRISOLVE real, dimension(IM,JM,LM) :: DP real(kind=MAPL_R8), dimension(IM,JM,LM) :: SX @@ -5320,6 +5363,8 @@ subroutine DIFFUSE(IM,JM,LM,RC) ! Pick the right exchange coefficients !------------------------------------- +OPT = .TRUE. + if ( (trim(name) /= 'S' ) .and. (trim(name) /= 'Q' ) .and. & (trim(name) /= 'QLLS') .and. (trim(name) /= 'QILS') .and. & (trim(name) /= 'U' ) .and. (trim(name) /= 'V' )) then @@ -5348,32 +5393,34 @@ subroutine DIFFUSE(IM,JM,LM,RC) elseif (trim(name) =='S') then CX => CT - DX => DKSS + DX => DKS AK => AKSS; BK => BKSS; CK => CKSS SX=S+YS elseif (trim(name)=='Q') then CX => CQ - DX => DKQQ + DX => DKQ AK => AKQQ; BK => BKQQ; CK => CKQQ SX=S+YQV elseif (trim(name)=='QLLS') then CX => CQ - DX => DKQQ + DX => DKQ AK => AKQQ; BK => BKQQ; CK => CKQQ SX=S+YQL +! OPT = .FALSE. elseif (trim(name)=='QILS') then CX => CQ - DX => DKQQ + DX => DKQ AK => AKQQ; BK => BKQQ; CK => CKQQ SX=S+YQI +! OPT = .FALSE. elseif (trim(name)=='U') then CX => CU - DX => DKUU + DX => DKV AK => AKUU; BK => BKUU; CK => CKUU SX=S+YU elseif (trim(name)=='V') then CX => CU - DX => DKUU + DX => DKV AK => AKUU; BK => BKUU; CK => CKUU SX=S+YV end if @@ -5382,7 +5429,7 @@ subroutine DIFFUSE(IM,JM,LM,RC) ! Solve for semi-implicit changes. This modifies SX ! ------------------------------------------------- - call VTRISOLVE(AK,BK,CK,SX,SG) + call VTRISOLVE(AK,BK,CK,SX,SG,OPT) ! Compute the surface fluxes !--------------------------- @@ -5631,11 +5678,11 @@ subroutine UPDATE(IM,JM,LM,LATS,RC) real, dimension(:,: ), pointer :: DSG, SF, SDF, SRFDIS real, dimension(:,: ), pointer :: HGTLM5, LM50M real, dimension(:,: ), pointer :: KETRB, KESRF, KETOP, KEINT - real, dimension(:,:,:), pointer :: DKS, DKV, DKQ, DKSS, DKUU, DKQQ, DKX, EKV, FKV + real, dimension(:,:,:), pointer :: DKS, DKV, DKQ, DKX, EKV, FKV real, dimension(:,:,:), pointer :: DPDTTRB real, dimension(:,:,:), pointer :: QTFLXTRB, SLFLXTRB, WSL, WQT, MFWSL, & MFWQT, TKH, UFLXTRB, VFLXTRB, QTX, SLX, & - SLFLXMF, QTFLXMF, MFAW + SLFLXMF, QTFLXMF, MFAW, TKEDISS integer :: KM, K, L, I, J logical :: FRIENDLY @@ -5664,7 +5711,7 @@ subroutine UPDATE(IM,JM,LM,LATS,RC) real :: SHVC_ALPHA, SHVC_EFFECT, SHVC_SCALING logical :: DO_SHVC logical :: ALLOC_TMP - integer :: KS + integer :: KS, DO_SHOC real :: HGT_SURFACE, WGTSUM ! For idealized SCM surface layer @@ -5681,6 +5728,7 @@ subroutine UPDATE(IM,JM,LM,LATS,RC) ALLOC_TMP = .FALSE. call MAPL_GetPointer(INTERNAL, TKH , 'TKH' , RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(INTERNAL, TKEDISS, 'TKEDISS' , RC=STATUS); VERIFY_(STATUS) call MAPL_GetPointer(EXPORT, QTX , 'QT' , ALLOC=.TRUE., RC=STATUS); VERIFY_(STATUS) call MAPL_GetPointer(EXPORT, SLX , 'SL' , ALLOC=.TRUE., RC=STATUS); VERIFY_(STATUS) @@ -5724,6 +5772,8 @@ subroutine UPDATE(IM,JM,LM,LATS,RC) call MAPL_GetResource( MAPL, HGT_SURFACE, 'HGT_SURFACE:', default=HGT_SURFACE, RC=STATUS ) VERIFY_(STATUS) + call MAPL_GetResource (MAPL, DO_SHOC, trim(COMP_NAME)//"_DO_SHOC:", default=0, RC=STATUS); VERIFY_(STATUS) + ! SHVC Resource parameters. SHVC_EFFECT can be set to zero to turn-off SHVC. ! SHVC_EFFECT = 1. is the tuned value for 2 degree horizontal resolution. ! It should be set to a lower number at higher resolution. @@ -5765,12 +5815,7 @@ subroutine UPDATE(IM,JM,LM,LATS,RC) VERIFY_(STATUS) call MAPL_GetPointer(INTERNAL, DKQ, 'DKQ', RC=STATUS) VERIFY_(STATUS) - call MAPL_GetPointer(INTERNAL, DKQQ, 'DKQQ', RC=STATUS) - VERIFY_(STATUS) - call MAPL_GetPointer(INTERNAL, DKSS, 'DKSS', RC=STATUS) - VERIFY_(STATUS) - call MAPL_GetPointer(INTERNAL, DKUU, 'DKUU', RC=STATUS) - VERIFY_(STATUS) + call MAPL_GetPointer(INTERNAL, EKV, 'EKV', RC=STATUS) VERIFY_(STATUS) call MAPL_GetPointer(INTERNAL, FKV, 'FKV', RC=STATUS) @@ -5966,6 +6011,10 @@ subroutine UPDATE(IM,JM,LM,LATS,RC) ! Loop over all quantities to be diffused. !----------------------------------------- + if (associated(INTDIS) .and. DO_SHOC /= 0) then + INTDIS(:,:,1:LM) = -1.*TKEDISS(:,:,1:LM)*DP(:,:,1:LM)/MAPL_CP + endif + TRACERS: do KS=1,KM K = KK(KS) @@ -6029,13 +6078,13 @@ subroutine UPDATE(IM,JM,LM,LATS,RC) RETURN_(ESMF_FAILURE) end if if( trim(NAME)=='QV' ) then - DKX => DKQQ + DKX => DKQ end if if( trim(NAME)=='S') then - DKX => DKSS + DKX => DKS end if if( trim(NAME)=='U' .or. trim(NAME)=='V' ) then - DKX => DKUU + DKX => DKV end if ! Update diffused quantity @@ -6054,10 +6103,12 @@ subroutine UPDATE(IM,JM,LM,LATS,RC) if( TYPE=='U' ) then if(associated(INTDIS)) then - DF(:,:,1:LM-1) = (0.5/MAPL_CP)*EKV(:,:,1:LM-1)*(SX(:,:,1:LM-1)-SX(:,:,2:LM))**2 ! Shear - DF(:,:, LM ) = 0.0 ! no shear at the surface, surface friction added later - INTDIS(:,:,1:LM-1) = INTDIS(:,:,1:LM-1) + DF - INTDIS(:,:,2:LM ) = INTDIS(:,:,2:LM ) + DF + if (DO_SHOC==0) then + DF(:,:,1:LM-1) = (0.5/MAPL_CP)*EKV(:,:,1:LM-1)*(SX(:,:,1:LM-1)-SX(:,:,2:LM))**2 ! Shear + DF(:,:, LM ) = 0.0 ! no shear at the surface, surface friction added later + INTDIS(:,:,1:LM-1) = INTDIS(:,:,1:LM-1) + DF + INTDIS(:,:,2:LM ) = INTDIS(:,:,2:LM ) + DF + endif ! Add surface dissipation to lower levels do J=1,JM do I=1,IM @@ -6266,6 +6317,7 @@ subroutine UPDATE(IM,JM,LM,LATS,RC) if (associated(QTFLXTRB)) QTFLXTRB = tmp3d + QTFLXMF if (associated(WQT)) WQT = 0.5*( tmp3d(:,:,1:LM)+tmp3d(:,:,0:LM-1) + QTFLXMF(:,:,1:LM)+QTFLXMF(:,:,0:LM-1) ) end if + if (associated(SLFLXTRB).or.associated(WSL)) then tmp3d(:,:,1:LM-1) = (SL(:,:,1:LM-1)-SL(:,:,2:LM))/(ZLO(:,:,1:LM-1)-ZLO(:,:,2:LM)) tmp3d(:,:,1:LM-1) = -1.*TKH(:,:,1:LM-1)*tmp3d(:,:,1:LM-1) @@ -6813,13 +6865,14 @@ end subroutine VTRISOLVESURF ! !INTERFACE: - subroutine VTRISOLVE ( A,B,C,Y,YG ) + subroutine VTRISOLVE ( A,B,C,Y,YG,OPT ) ! !ARGUMENTS: real, dimension(:,:,:), intent(IN ) :: A, B, C real(kind=MAPL_R8), dimension(:,:,:), intent(INOUT) :: Y real, dimension(:,:), intent(IN) :: YG + logical, intent(IN) :: OPT ! !DESCRIPTION: Solves tridiagonal system that has been LU decomposed ! $LU x = f$. This is done by first solving $L g = f$ for $g$, and @@ -6864,9 +6917,10 @@ subroutine VTRISOLVE ( A,B,C,Y,YG ) if(size(YG)>0) then Y(:,:,LM) = (Y(:,:,LM) - C(:,:,LM) * YG )*B(:,:,LM) - else + else if (OPT) then Y(:,:,LM) = Y(:,:,LM)*B(:,:,LM-1)/(B(:,:,LM-1) - A(:,:,LM)*(1.0+C(:,:,LM-1)*B(:,:,LM-1) )) - ! Y(:,:,LM) = Y(:,:,LM)*B(:,:,LM)/( 1.0+C(:,:,LM)*B(:,:,LM) ) ! Alternate formulation + else + Y(:,:,LM) = Y(:,:,LM)*B(:,:,LM)/( 1.0+C(:,:,LM)*B(:,:,LM) ) ! Alternate formulation endif do L = LM-1,1,-1 diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSturbulence_GridComp/edmf.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSturbulence_GridComp/edmf.F90 index 1b60b8cca..13bb5b5b2 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSturbulence_GridComp/edmf.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSturbulence_GridComp/edmf.F90 @@ -55,6 +55,7 @@ module edmf_mod real :: MFLIMFAC real :: ICE_RAMP real :: PRCPCRIT + real :: TREFF endtype EDMFPARAMS_TYPE type (EDMFPARAMS_TYPE) :: MFPARAMS @@ -74,7 +75,6 @@ SUBROUTINE RUN_EDMF(its,ite, jts,jte, kts,kte, dt, & ! Inputs t3, & thl3, & thv3, & - qt3, & qv3, & ql3, & qi3, & @@ -137,7 +137,6 @@ SUBROUTINE RUN_EDMF(its,ite, jts,jte, kts,kte, dt, & ! Inputs V3, & T3, & THL3, & - QT3, & THV3, & QV3, & QL3, & @@ -361,7 +360,7 @@ SUBROUTINE RUN_EDMF(its,ite, jts,jte, kts,kte, dt, & ! Inputs end if end do lts = lts - thv3(IH,JH,kte) - L0 = L0/( 1.0 + (mfparams%ent0lts/mfparams%ent0-1.)*(0.5+0.5*tanh(0.3*(lts-19.))) ) + L0 = L0/( 1.0 + (mfparams%ent0lts/mfparams%ent0-1.)*(0.5+0.5*tanh(0.3*(lts-18.5))) ) end if else ! if mfparams%ET not 2 L0 = mfparams%L0 @@ -376,7 +375,6 @@ SUBROUTINE RUN_EDMF(its,ite, jts,jte, kts,kte, dt, & ! Inputs v(k)=v3(IH,JH,kte-k+kts) thl(k)=thl3(IH,JH,kte-k+kts) thv(k)=thv3(IH,JH,kte-k+kts) - qt(k)=qt3(IH,JH,kte-k+kts) qv(k)=qv3(IH,JH,kte-k+kts) ql(k)=ql3(IH,JH,kte-k+kts) qi(k)=qi3(IH,JH,kte-k+kts) @@ -385,7 +383,6 @@ SUBROUTINE RUN_EDMF(its,ite, jts,jte, kts,kte, dt, & ! Inputs ui(k) = 0.5*( u3(IH,JH,kte-k+kts) + u3(IH,JH,kte-k+kts-1) ) vi(k) = 0.5*( v3(IH,JH,kte-k+kts) + v3(IH,JH,kte-k+kts-1) ) thli(k) = 0.5*( thl3(IH,JH,kte-k+kts) + thl3(IH,JH,kte-k+kts-1) ) - qti(k) = 0.5*( qt3(IH,JH,kte-k+kts) + qt3(IH,JH,kte-k+kts-1) ) qvi(k) = 0.5*( qv3(IH,JH,kte-k+kts) + qv3(IH,JH,kte-k+kts-1) ) qli(k) = 0.5*( ql3(IH,JH,kte-k+kts) + ql3(IH,JH,kte-k+kts-1) ) qii(k) = 0.5*( qi3(IH,JH,kte-k+kts) + qi3(IH,JH,kte-k+kts-1) ) @@ -393,7 +390,6 @@ SUBROUTINE RUN_EDMF(its,ite, jts,jte, kts,kte, dt, & ! Inputs ui(k) = u3(IH,JH,kte-k+kts-1) vi(k) = v3(IH,JH,kte-k+kts-1) thli(k) = thl3(IH,JH,kte-k+kts-1) - qti(k) = qt3(IH,JH,kte-k+kts-1) qvi(k) = qv3(IH,JH,kte-k+kts-1) qli(k) = ql3(IH,JH,kte-k+kts-1) qii(k) = qi3(IH,JH,kte-k+kts-1) @@ -403,17 +399,17 @@ SUBROUTINE RUN_EDMF(its,ite, jts,jte, kts,kte, dt, & ! Inputs ui(kte) = u(kte) vi(kte) = v(kte) thli(kte) = thl(kte) - qti(kte) = qt(kte) qvi(kte) = qv(kte) qli(kte) = ql(kte) qii(kte) = qi(kte) ui(kts-1) = u(kts) vi(kts-1) = v(kts) thli(kts-1) = thl(kts) ! approximate - qti(kts-1) = qt(kts) qvi(kts-1) = qv(kts) qli(kts-1) = ql(kts) qii(kts-1) = qi(kts) + qt = qv+ql+qi + qti = qvi+qli+qii DO k=kts-1,kte rhoe(k) = rhoe3(IH,JH,kte-k+kts-1) @@ -618,7 +614,7 @@ SUBROUTINE RUN_EDMF(its,ite, jts,jte, kts,kte, dt, & ! Inputs Wn2=UPW(K-1,I)**2+2.*MFPARAMS%WA*B*(ZW(k)-ZW(k-1)) ELSE EntW=exp(-2.*WP*(ZW(k)-ZW(k-1))) - Wn2=EntW*UPW(k-1,I)**2+MFPARAMS%WA*B/WP*(1.-EntW) + Wn2=EntW*UPW(k-1,I)**2+(1.-EntW)*MFPARAMS%WA*B/WP END IF @@ -703,7 +699,7 @@ SUBROUTINE RUN_EDMF(its,ite, jts,jte, kts,kte, dt, & ! Inputs K = KTS tmp = 0. tmp2 = 0. - DO WHILE (ZW(K).lt.70. .and. K.lt.KTE) + DO WHILE (ZW(K).lt.100. .and. K.lt.KTE) tmp = tmp + 0.5*SUM(UPA(K,:)*UPW(K,:)*UPW(K,:)) tmp2 = tmp2 + TKE3(IH,JH,KTE-K+KTS) ! UPW(K,:) = UPW(K,:)*exp(-(100.-ZW(K))**2/1e4) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSturbulence_GridComp/shoc.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSturbulence_GridComp/shoc.F90 index 92cb6eac1..9299604c9 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSturbulence_GridComp/shoc.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSturbulence_GridComp/shoc.F90 @@ -370,15 +370,14 @@ subroutine tke_shoc() ! TKE boyancy production term. wthv_sec (buoyancy flux) is calculated in Moist GridComp. - wrk = 0.5 * (tkh(i,j,ku)+tkh(i,j,kd)) - if (shocparams%BUOYOPT==2) then 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) + wrk = 0.5 * (tkh(i,j,ku)*brunt_edge(i,j,ku)+tkh(i,j,kd)*brunt_edge(i,j,kd)) + a_prod_bu = -1.*wrk + (ggr / thv(i,j,k))*wthv_mf(i,j,k) end if - buoy_sgs = brunt(i,j,k) + buoy_sgs = 0.5*(brunt_edge(i,j,ku)+brunt_edge(i,j,kd)) ! buoy_sgs = - a_prod_bu / (wrk + 0.0001) ! tkh is eddy thermal diffussivity !Compute $c_k$ (variable Cee) for the TKE dissipation term following Eq. 11 in Deardorff (1980) @@ -390,7 +389,7 @@ subroutine tke_shoc() Cee = Cek* (pt19 + pt51*smix/grd) - wrk = 0.5 * wrk * (prnum(i,j,ku) + prnum(i,j,kd)) + wrk = 0.25 * (tkh(i,j,ku)+tkh(i,j,kd)) * (prnum(i,j,ku) + prnum(i,j,kd)) if (nx.eq.1) then a_prod_sh = min(min(tkhmax,wrk)*def2(i,j,k),0.0001) ! TKE shear production term else @@ -442,16 +441,13 @@ subroutine tke_shoc() do j=1,ny do i=1,nx ! Calculate "return-to-isotropy" eddy dissipation time scale, see Eq. 8 in BK13 - if (brunt_edge(i,j,k) <= bruntmin) then + if (brunt_edge(i,j,k) <= 1e-5 .or. zl(i,j,k).lt.0.5*zpbl(i,j)) then isotropy(i,j,k) = max(30.,min(max_eddy_dissipation_time_scale,0.5*(tscale1(i,j,k)+tscale1(i,j,k-1)))) else wrk = 0.5*(tscale1(i,j,k)+tscale1(i,j,k-1)) isotropy(i,j,k) = max(30.,min(max_eddy_dissipation_time_scale,wrk/(1.0+lambda*brunt_edge(i,j,k)*wrk*wrk))) ! isotropy(i,j,k) = max(30.,min(max_eddy_dissipation_time_scale,wrk/(1.0+lambda*0.5*(brunt(i,j,k)+brunt(i,j,k-1))*wrk*wrk))) endif -! 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 (tke(i,j,k).lt.2e-4) isotropy(i,j,k) = 30. wrk1 = ck / prnum(i,j,k) @@ -494,7 +490,7 @@ subroutine calc_numbers() DU = (U(:,:,1:nzm-1) - U(:,:,2:nzm))**2 + & ! shear on edges (V(:,:,1:nzm-1) - V(:,:,2:nzm))**2 - DU = MAX( SQRT(DU) / adzi(:,:,1:nzm-1), 0.005 ) + DU = MIN( MAX( SQRT(DU) / adzi(:,:,1:nzm-1), 0.005 ), 0.025 ) RI = 0.0 RI(:,:,2:nz-1) = ggr*( (THV(:,:,2:nzm) - THV(:,:,1:nzm-1)) / adzi(:,:,1:nzm-1) ) & @@ -654,8 +650,8 @@ subroutine eddy_length() ! Calculate the measure of PBL depth, Eq. 11 in BK13 (Is this really PBL depth?) ! cldbasek(:,:) = 1 - do j=1,ny - do i=1,nx +! do j=1,ny +! do i=1,nx ! do k=1,nzm ! if (zl(i,j,k).gt.3000. .or. cld_sgs(i,j,k).gt.0.01) exit @@ -690,11 +686,11 @@ subroutine eddy_length() ! cldbasek(i,j) = cldbasek(i,j) + 1 ! end do - kk = 1 - do while (zl(i,j,kk) .lt. 3000. .or. kk.eq.nzm) - kk = kk + 1 - end do - lts(i,j) = thv(i,j,kk) - thv(i,j,1) +! kk = 1 +! do while (zl(i,j,kk) .lt. 3000. .or. kk.eq.nzm) +! kk = kk + 1 +! end do +! lts(i,j) = thv(i,j,kk) - thv(i,j,1) ! Alternate cloud base calculation ! tep = tabs(i,j,1) @@ -707,8 +703,8 @@ subroutine eddy_length() ! end do ! zcb(i,j) = max(200.,zl(i,j,kk-1)) !kk-1 is highest level *before* condensation ! if (nx.eq.1) print *,'zcb=',zcb(i,j) - enddo - enddo +! enddo +! enddo !Calculate length scale outside of cloud, Eq. 10 in BK13 (Eq. 4.12 in Pete's dissertation) @@ -802,22 +798,27 @@ subroutine eddy_length() + (bbb*fac_cond-tabs(i,j,k))*(qpl(i,j,k)-qpl(i,j,k-1)) & + (bbb*fac_sub -tabs(i,j,k))*(qpi(i,j,k)-qpi(i,j,k-1)) ) end if + end do + + end do + end do ! Reduction of mixing length in the stable regions (where B.-V. freq. > 0) is required. ! Here we find regions of Brunt-Vaisalla freq. > 0 for later use. - if (brunt(i,j,k) >= bruntmin) then - brunt2(i,j,k) = brunt(i,j,k) - else - brunt2(i,j,k) = bruntmin + brunt_edge(:,:,1) = brunt_edge(:,:,2) + brunt_edge(:,:,nz) = brunt_edge(:,:,nzm) + do i=1,nx + do j=1,ny + do k=1,nzm + brunt2(i,j,k) = (1.-cld_sgs(i,j,k))*0.5*(brunt_edge(i,j,k-1)+brunt_edge(i,j,k)) + cld_sgs(i,j,k)*min(brunt_edge(i,j,k-1),brunt_edge(i,j,k)) + if (brunt2(i,j,k) < 1e-5 .or. zl(i,j,k).lt.0.5*zpbl(i,j)) then + brunt2(i,j,k) = 1e-10 endif - end do - end do end do - brunt_edge(:,:,1) = brunt_edge(:,:,2) - brunt_edge(:,:,nz) = brunt_edge(:,:,nzm) + brunt2(:,:,1) = brunt2(:,:,2) brunt2(:,:,nzm) = brunt2(:,:,nzm-1) @@ -889,17 +890,16 @@ subroutine eddy_length() ! smixt1(i,j,k) = sqrt(400.*tkes*vonk*zl(i,j,k))*shocparams%LENFAC1 ! original SHOC, includes TKE ! Turbulent length scale - 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(brunt2(i,j,k))) + smixt2(i,j,k) = sqrt(l_par(i,j,k)*400.*tkes)*shocparams%LENFAC2 + ! Stability length scale, reduced influence below 300m + smixt3(i,j,k) = max(0.05,tkes)*shocparams%LENFAC3/(sqrt(brunt2(i,j,k))) !=== Combine component length scales === if (shocparams%LENOPT .eq. 1) then ! JPL blending approach (w/SHOC length scales) - wrk1 = 2./(1./smixt2(i,j,k)+1./smixt3(i,j,k)) + wrk1 = sqrt(3./(1./smixt2(i,j,k)**2+1./smixt3(i,j,k)**2)) if (zl(i,j,k).lt.300.) then - smixt(i,j,k) = wrk1 + (smixt1(i,j,k)-wrk1)*exp(-(zl(i,j,k)/100.)**2) + smixt(i,j,k) = wrk1 + (smixt1(i,j,k)-wrk1)*exp(-(zl(i,j,k)/60.)) else smixt(i,j,k) = wrk1 end if @@ -919,13 +919,11 @@ subroutine eddy_length() end if ! Enforce minimum and maximum length scales - wrk = 0.1*min(200.,adzl(i,j,k)) ! Minimum 0.1 of local dz (up to 200 m) - if (zl(i,j,k) .lt. 5000.) then - smixt(i,j,k) = max(wrk, smixt(i,j,k)) - else if (zl(i,j,k) .lt. 9500) then ! Between 5-10 km the max length scale reduces with height - smixt(i,j,k) = max(wrk, min(max_eddy_length_scale*(1e4-zl(i,j,k))/5e3,smixt(i,j,k))) - else - smixt(i,j,k) = max(wrk, min(max_eddy_length_scale*0.1,smixt(i,j,k))) + wrk = 20. !0.5*min(100.,adzl(i,j,k)) ! Minimum 0.1 of local dz (up to 200 m) + if (zl(i,j,k) .lt. 2000.) then + smixt(i,j,k) = max(wrk, smixt(i,j,k)) + else if (zl(i,j,k).gt.zpbl(i,j)) then ! if above 2 km and dry CBL top, cap length scale + smixt(i,j,k) = max(wrk, min(200.,smixt(i,j,k))) end if end do end do @@ -1033,8 +1031,8 @@ subroutine update_moments( IM, JM, LM, & ! in hl2tune, & qt2tune, & hlqt2tune, & - qt3_tscale, & - afrc_tscale,& + skew_tgen, & + skew_tdis, & docanuto ) @@ -1077,8 +1075,8 @@ subroutine update_moments( IM, JM, LM, & ! in real, intent(in ) :: HL2TUNE, & ! tuning parameters HLQT2TUNE, & QT2TUNE, & - QT3_TSCALE, & - AFRC_TSCALE + SKEW_TGEN, & + SKEW_TDIS integer, intent(in ) :: DOPROGQT2, & ! prognostic QT2 switch DOCANUTO @@ -1138,7 +1136,7 @@ subroutine update_moments( IM, JM, LM, & ! in ! Second moment of total water mixing ratio. Eq 3 in BK13 qtgrad(:,:,k) = wrk2 / (ZL(:,:,k)-ZL(:,:,k+1)) - qt2_edge(:,:,k) = (KH(:,:,k)*qtgrad(:,:,k)-MFWQT(:,:,k)-WQT_DC(:,:,k))*qtgrad(:,:,k) ! gradient production + qt2_edge(:,:,k) = (KH(:,:,k)*qtgrad(:,:,k)-MFWQT(:,:,k)-0.*WQT_DC(:,:,k))*qtgrad(:,:,k) ! gradient production qt2_edge_nomf(:,:,k) = (KH(:,:,k)*qtgrad(:,:,k))*qtgrad(:,:,k) ! gradient production ! Covariance of total water mixing ratio and liquid/ice water static energy. Eq 5 in BK13 @@ -1157,16 +1155,6 @@ subroutine update_moments( IM, JM, LM, & ! in qtgrad(:,:,0) = qtgrad(:,:,1) - ! Update PDF_A - if (AFRC_TSCALE.gt.0.) then - pdf_a = (pdf_a+mffrc)/(1.+DT/AFRC_TSCALE) - else - pdf_a = pdf_a/(1.-DT/AFRC_TSCALE) - end if - where (mffrc.gt.pdf_a) - pdf_a = mffrc - end where - pdf_a = min(0.5,max(0.,pdf_a)) do k=1,LM @@ -1192,7 +1180,7 @@ subroutine update_moments( IM, JM, LM, & ! in else onemmf = 1.0 - MFFRC(:,:,k) - w2(:,:,k) = onemmf*0.667*TKE(:,:,k) + MFW2(:,:,k) + w2(:,:,k) = onemmf*0.667*TKE(:,:,k) !+ MFW2(:,:,k) ! hl2(:,:,k) = onemmf*0.5*( hl2_edge(:,:,kd) + hl2_edge(:,:,ku) ) !+ MFHL2(:,:,k) hl2(:,:,k) = 0.5*( hl2_edge(:,:,kd) + hl2_edge(:,:,ku) ) @@ -1202,11 +1190,11 @@ subroutine update_moments( IM, JM, LM, & ! in if (DOPROGQT2 /= 0) then wrk3 = QT2TUNE*1.5e-4 ! dissipation qt2(:,:,k) = (qt2(:,:,k)+DT*wrk1) / (1. + DT*wrk3) - qt2diag(:,:,k) = QT2TUNE*ISOTROPY(:,:,k)*0.5*(qt2_edge(:,:,kd)+qt2_edge(:,:,ku)) + qt2diag(:,:,k) = QT2TUNE*ISOTROPY(:,:,k)*0.5*(qt2_edge_nomf(:,:,kd)+qt2_edge_nomf(:,:,ku)) else ! qt2(:,:,k) = QT2TUNE*ISOTROPY(:,:,k)*wrk1 + MFQT2(:,:,k) qt2(:,:,k) = QT2TUNE*ISOTROPY(:,:,k)*wrk1 - qt2diag(:,:,k) = QT2TUNE*ISOTROPY(:,:,k)*0.5*(qt2_edge_nomf(:,:,kd)+qt2_edge_nomf(:,:,ku)) + qt2diag(:,:,k) = 1.0*ISOTROPY(:,:,k)*0.5*(qt2_edge_nomf(:,:,kd)+qt2_edge_nomf(:,:,ku)) end if hlqt(:,:,k) = onemmf*0.5*( hlqt_edge(:,:,kd) + hlqt_edge(:,:,ku) ) + MFHLQT(:,:,k) @@ -1216,12 +1204,12 @@ subroutine update_moments( IM, JM, LM, & ! in end if - 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)) + 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, 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) + ! Restrict QT variance, 2-25% of total water. + qt2(:,:,k) = max(min(qt2(:,:,k),(0.25*QT(:,:,k))**2),(0.02*QT(:,:,k))**2) + qt2diag(:,:,k) = max(min(qt2diag(:,:,k),(0.25*QT(:,:,k))**2),(0.02*QT(:,:,k))**2) hl2(:,:,k) = max(min(hl2(:,:,k),HL2MAX),HL2MIN) hl2diag(:,:,k) = max(min(hl2diag(:,:,k),HL2MAX),HL2MIN) @@ -1232,8 +1220,20 @@ subroutine update_moments( IM, JM, LM, & ! in end do + ! Update PDF_A + if (SKEW_TDIS.gt.0.) then +! pdf_a = (pdf_a+mffrc+2.*0.5*(cnv_mfc(:,:,1:LM)+cnv_mfc(:,:,0:LM-1)))/(1.+DT/AFRC_TSCALE) + pdf_a = (pdf_a+mffrc*DT/SKEW_TGEN)/(1.+DT/SKEW_TDIS) + else + pdf_a = pdf_a/(1.-DT/SKEW_TDIS) + end if + where (mffrc.gt.pdf_a) + pdf_a = mffrc + end where + pdf_a = min(0.5,max(0.,pdf_a)) + if (DOCANUTO==0) then - qt3 = ( qt3 + max(MFQT3,0.) ) / ( 1. + DT/QT3_TSCALE ) + qt3 = ( qt3 + max(MFQT3,0.)*DT/SKEW_TGEN ) / ( 1. + DT/SKEW_TDIS ) hl3 = MFHL3 w3 = MFW3 else @@ -1273,10 +1273,6 @@ subroutine update_moments( IM, JM, LM, & ! in ! brunt = (bet(i,j,k)/thedz)*(thv(i,j,kc)-thv(i,j,kb)) -! if (abs(thedz).le.1e-10) thedz = sign(1e-10,thedz) -! if (abs(thedz).eq.1e-10) print *,'thedz' -! if (abs(thedz2).le.1e-10) thedz2 = sign(1e-10,thedz2) -! if (abs(thedz2).eq.1e-10) print *,'thedz2' thedz = 1. / thedz thedz2 = 1. / thedz2 @@ -1288,8 +1284,8 @@ subroutine update_moments( IM, JM, LM, & ! in avew = 0.5*(0.667*TKE(i,j,k)+0.667*TKE(i,j,kb)) if (abs(avew).ge.1e10) avew = sign(1e10,avew) -! if (abs(avew).eq.1e10) print *,'avew' - cond = 1.2*sqrt(max(1.0e-16,2.*avew*avew*avew)) + + cond = 1.2*sqrt(max(1.0e-20,2.*avew*avew*avew)) wrk1b = bet2*iso wrk2b = thedz2*wrk1b*wrk1b*iso wrk3b = hl2diag(i,j,kc) - hl2diag(i,j,kb) @@ -1313,20 +1309,20 @@ subroutine update_moments( IM, JM, LM, & ! in ! Compute the "omega" terms, see Eq. 6 in C01 (B.6 in Pete's dissertation) dum = 1.-a5*buoy_sgs2 - if (abs(dum).le.1e-16) dum = sign(1e-16,dum) -! if (abs(dum).eq.1e-16) print *,'1.-a5*buoy_sgs2' + if (abs(dum).le.1e-20) dum = sign(1e-20,dum) + omega0 = a4 / dum omega1 = omega0 / (c+c) omega2 = omega1*f3+(5./4.)*omega0*f4 ! Compute the X0, Y0, X1, Y1 terms, see Eq. 5 a-b in C01 (B.5 in Pete's dissertation) dum = 1.-(a1+a3)*buoy_sgs2 - if (abs(dum).le.1e-16) dum = sign(1e-16,dum) -! if (abs(dum).eq.1e-16) print *,'1.-(a1+a3)*buoy_sgs2' + if (abs(dum).le.1e-20) dum = sign(1e-20,dum) + wrk1b = 1.0 / dum dum = 1.-a3*buoy_sgs2 - if (abs(dum).le.1e-16) dum = sign(1e-16,dum) -! if (abs(dum).eq.1e-16) print *,'1.-a3*buoy_sgs2' + if (abs(dum).le.1e-20) dum = sign(1e-20,dum) + wrk2b = 1.0 / dum X0 = wrk1b * (a2*buoy_sgs2*(1.-a3*buoy_sgs2)) Y0 = wrk2b * (2.*a2*buoy_sgs2*X0) @@ -1342,8 +1338,8 @@ subroutine update_moments( IM, JM, LM, & ! in ! than the estimate - limit w3. dum = c-1.2*X0+AA0 - if (abs(dum).le.1e-16) dum = sign(1e-16,dum) -! if (abs(dum).eq.1e-16) print *,'c-1.2*X0+AA0=',dum + if (abs(dum).le.1e-20) dum = sign(1e-20,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) @@ -1357,11 +1353,11 @@ subroutine update_moments( IM, JM, LM, & ! in w3can(i,j,LM) = w3can(i,j,LM-1) enddo enddo - w3 = w3can +! 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 diff --git a/GEOSagcm_GridComp/GEOSsuperdyn_GridComp/GEOSdatmodyn_GridComp/GEOS_DatmoDynGridComp.F90 b/GEOSagcm_GridComp/GEOSsuperdyn_GridComp/GEOSdatmodyn_GridComp/GEOS_DatmoDynGridComp.F90 index 6df0afe2f..dfed49e33 100644 --- a/GEOSagcm_GridComp/GEOSsuperdyn_GridComp/GEOSdatmodyn_GridComp/GEOS_DatmoDynGridComp.F90 +++ b/GEOSagcm_GridComp/GEOSsuperdyn_GridComp/GEOSdatmodyn_GridComp/GEOS_DatmoDynGridComp.F90 @@ -1629,7 +1629,7 @@ subroutine RUN ( GC, IMPORT, EXPORT, CLOCK, RC ) ! whenever datmodyn starts using gocart, this will need to be a real value -- ! see fvdycore for example if(associated(DUMMYAREA)) then - DUMMYAREA=1.0 + DUMMYAREA=1e10 end if if(associated(DUMMYDXC)) then