diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/ConvPar_GF2020.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/ConvPar_GF2020.F90 index d2e0e6f36..12e6a3b6e 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/ConvPar_GF2020.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/ConvPar_GF2020.F90 @@ -11079,14 +11079,13 @@ subroutine get_liq_ice_number_conc(itf,ktf,its,ite, kts,kte,ierr,ktop& outnliq(i,k) = max(0.0, make_DropletNumber(tqliq, nwfa (i,k))/rho(i,k)) enddo - !-- convert in tendencies - outnice = outnice * dtinv ! unit [1/s] - outnliq = outnliq * dtinv ! unit [1/s] - !--- for update - ! nwfa =nwfa + outnliq*dtime - ! nifa =nifa + outnice*dtime - + !-- convert in tendencies + outnice = outnice * dtinv ! unit [1/s] + outnliq = outnliq * dtinv ! unit [1/s] enddo + !--- for update + ! nwfa =nwfa + outnliq*dtime + ! nifa =nifa + outnice*dtime end subroutine get_liq_ice_number_conc !DSM { diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/Process_Library.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/Process_Library.F90 index e59bb0848..2310c2424 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/Process_Library.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/Process_Library.F90 @@ -253,7 +253,7 @@ module GEOSmoist_Process_Library public :: find_cldtop, find_cldbase, gw_prof public :: make_IceNumber, make_DropletNumber, make_RainNumber public :: dissipative_ke_heating - public :: pdffrac, pdfcondensate, partition_dblgss + public :: pdffrac, pdfcondensate, precalc_dblgss, partition_dblgss, partition_dblgss2 public :: SIGMA_DX, SIGMA_EXP public :: CNV_FRACTION_MIN, CNV_FRACTION_MAX public :: SH_MD_DP, DBZ_VAR_INTERCP, DBZ_LIQUID_SKIN, LIQ_RADII_PARAM, ICE_RADII_PARAM @@ -1532,6 +1532,7 @@ subroutine partition_dblgss( fQi, & ! IN tabs, & qwv, & qc, & ! OUT + qsat, & ! qi, & omega, & ! IN zl, & @@ -1586,6 +1587,7 @@ subroutine partition_dblgss( fQi, & ! IN ! real, intent(in ) :: DT ! timestep [s] real, intent(in ) :: tabs ! absolute temperature [K] real, intent(in ) :: qwv ! specific humidity [kg kg-1] + real, intent(in ) :: qsat real, intent( out) :: qc ! liquid+ice condensate [kg kg-1] real, intent(in ) :: omega ! resolved pressure velocity real, intent(in ) :: zl ! layer heights [m] @@ -1677,10 +1679,10 @@ subroutine partition_dblgss( fQi, & ! IN w_first = - rog * omega * thv / pval ! Initialize cloud variables to zero - diag_qn = 0.0 - diag_frac = 0.0 - diag_ql = 0.0 - diag_qi = 0.0 +! diag_qn = 0.0 +! diag_frac = 0.0 +! diag_ql = 0.0 +! diag_qi = 0.0 pkap = (pval/100000.0) ** kapa @@ -1708,7 +1710,16 @@ subroutine partition_dblgss( fQi, & ! IN skew_qw = 0. endif -! Find parameters of the double Gaussian PDF of vertical velocity +! qsat = GEOS_QSAT( tabs, pval ) +! if (skew_qw.lt.1e-3 .and. total_water+3.*sqrtqt.lt.qsat) then + if (qt3.lt.1e-12 .and. total_water+4.*sqrtqt.lt.(1.-0.14*sqrtthl)*qsat) then + wqls = 0. + wqis = 0. + cld_sgs = 0. + qc = 0. + else + + ! Find parameters of the double Gaussian PDF of vertical velocity ! aterm = pdf_a @@ -2130,6 +2141,8 @@ subroutine partition_dblgss( fQi, & ! IN 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) + end if ! small RH conditional + ! diagnostic buoyancy flux. Includes effects from liquid water, ice ! condensate, liquid & ice precipitation wrk = epsv * thv @@ -2143,6 +2156,152 @@ subroutine partition_dblgss( fQi, & ! IN end subroutine partition_dblgss + + subroutine precalc_dblgss(a, wqt, wth, tbar, t2bar, t3bar, qt2bar, qt3bar, & + w2bar, w3bar, beta, rwqt, rwth, rtqt, t1, t2, & + sigt1, sigt2, qt1, qt2, sigqt1, sigqt2, w1, w2, & + sigw1, sigw2) + ! Input arguments + real(4), intent(in) :: a, wqt, wth, tbar, t2bar, t3bar, qt2bar, qt3bar, w2bar, w3bar + + ! Output arguments + real(4), intent(out) :: beta, rwqt, rwth, rtqt, t1, t2, sigt1, sigt2 + real(4), intent(out) :: qt1, qt2, sigqt1, sigqt2, w1, w2, sigw1, sigw2 + + ! Local variables + real(4) :: tmp, fac, qt3bar_local + + ! Calculate beta + beta = (MAPL_RGAS/MAPL_RVAP) * (MAPL_ALHL/(MAPL_RGAS*tbar)) * (MAPL_ALHL/(MAPL_CP*tbar)) + + ! If some skewness is present + if (a > 0.001 .and. a /= 0.5 .and. & + (qt3bar > 0.0 .or. t3bar > 0.0 .or. w3bar > 0.0)) then + + ! Empirical adjustment to ensure realizability + qt3bar_local = min(qt3bar, (sqrt(qt2bar) * (0.6 + 1.5 * exp(-10.0 * a)))**3) + + tmp = 1.0 - (1.0 - a) - a**3 / (1.0 - a)**2 + + if (tmp /= 0.0) then + t2 = sign(1.0, t3bar) * abs(t3bar / tmp)**(1.0/3.0) + w2 = sign(1.0, w3bar) * abs(w3bar / tmp)**(1.0/3.0) + qt2 = abs(qt3bar_local / tmp)**(1.0/3.0) + else + t2 = sign(1.0, t3bar) * abs(t3bar)**(1.0/3.0) + w2 = sign(1.0, w3bar) * abs(w3bar)**(1.0/3.0) + qt2 = abs(qt3bar_local)**(1.0/3.0) + end if + + t1 = a * t2 / (a - 1.0) + w1 = a * w2 / (a - 1.0) + + fac = sqrt(a / (1.0 - a)) + tmp = max(0.0, t2bar - (1.0 - a) * t1**2 - a * t2**2) + sigt1 = sqrt(tmp / (1.0 - a + a * fac)) + sigt2 = fac * sigt1 + tmp = max(0.0, w2bar - (1.0 - a) * w1**2 - a * w2**2) + sigw1 = sqrt(tmp / (1.0 - a + a * fac)) + sigw2 = fac * sigw1 + + else + + t1 = 0.0 + t2 = 0.0 + sigt1 = sqrt(t2bar) + sigt2 = sigt1 + w1 = 0.0 + w2 = 0.0 + sigw1 = sqrt(w2bar) + sigw2 = sigw1 + + end if + + if (a > 0.001 .and. qt3bar > 0.0 .and. a /= 0.5) then + + qt1 = a * qt2 / (1.0 - a) + + fac = sqrt(a / (1.0 - a)) + tmp = max(0.0, qt2bar - (1.0 - a) * qt1**2 - a * qt2**2) + sigqt1 = sqrt(tmp / (1.0 - a + a * fac)) + sigqt2 = fac * sigqt1 + + else + + qt1 = 0.0 + qt2 = 0.0 + sigqt1 = sqrt(qt2bar) + sigqt2 = sigqt1 + + end if + + ! Set correlation coefficients + rwqt = 0. + rwth = 0. + rtqt = -0.2 + + end subroutine precalc_dblgss + + subroutine partition_dblgss2(exner, a, beta, rwqt, rwth, rtqt, t1pert, t2pert, & + sigt1, sigt2, qtbar, qt1, qt2, sigqt1, sigqt2, & + w1, w2, sigw1, sigw2, qsat, dqsat, qc, cld, wqc) + ! Input arguments + real(4), intent(in) :: exner, a, beta, rwqt, rwth, rtqt, t1pert, t2pert + real(4), intent(in) :: sigt1, sigt2, qtbar, qt1, qt2, sigqt1, sigqt2 + real(4), intent(in) :: w1, w2, sigw1, sigw2, qsat, dqsat + + ! Output arguments + real(4), intent(out) :: qc, cld, wqc + + ! Local variables + real(4) :: qsat1, qsat2, cq, ct, s1, s2, sigs1, sigs2 + real(4) :: cld1, cld2, qc1, qc2 + + real(4), parameter :: sqrt2 = 1.4142135 + real(4), parameter :: sqrt2pi = 2.5066282 + + ! Calculate saturation mixing ratio for plume 1 + qsat1 = qsat + dqsat * t1pert + + cq = 1.0 / (1.0 + beta * qsat1) + ct = (1.0 + beta * qtbar) * exner * & + (MAPL_CP/MAPL_ALHL) * beta * qsat * cq**2 + + s1 = qtbar + qt1 - qsat1 * (1.0 + beta * (qtbar + qt1)) * cq + + sigs1 = max(1e-6*qtbar,sqrt(ct**2 * sigt1**2 + cq**2 * sigqt1**2 - & + 2.0 * ct * cq * sigt1 * sigqt1 * rtqt)) + + cld1 = 0.5 * (1.0 + erf(s1 / (sqrt2 * sigs1))) + + qc1 = max(0.0, s1 * cld1 + (sigs1 / sqrt2pi) * exp(-0.5 * (s1/sigs1)**2)) + + if (a > 0.001) then ! If skewed (not single gaussian) + qsat2 = qsat + dqsat * t2pert + s2 = qtbar + qt2 - qsat2 * (1.0 + beta * (qtbar + qt2)) / & + (1.0 + beta * qsat2) + sigs2 = max(1e-6*qtbar,sqrt(ct**2 * sigt2**2 + cq**2 * sigqt2**2 - & + 2.0 * ct * cq * sigt2 * sigqt2 * rtqt)) + cld2 = 0.5 * (1.0 + erf(s2 / (sqrt2 * sigs2))) + qc2 = max(0.0, s2 * cld2 + (sigs2 / sqrt2pi) * exp(-0.5 * (s2/sigs2)**2)) + else + cld2 = cld1 + qc2 = qc1 + end if + + qc = qc1 * (1.0 - a) + qc2 * a + + wqc = (1.0 - a) * w1 * qc1 + a * w2 * qc2 +! wqc = (1.0 - a) * (w1 * qc1 + cld1 * (cq * sigw1 * sigqt1 * rwqt - & +! ct * sigw1 * sigt1 * rwth)) + & +! a * (w2 * qc2 + cld2 * (cq * sigw2 * sigqt2 * rwqt - & +! ct * sigw2 * sigt2 * rwth)) + + cld = (1.0 - a) * cld1 + a * cld2 + + end subroutine partition_dblgss2 + + subroutine hystpdf( & DT , & ALPHA , & @@ -2167,9 +2326,9 @@ subroutine hystpdf( & QT2 , & HLQT , & W3 , & - W2 , & - MFQT3 , & - MFHL3 , & + W2bar , & + QT3 , & + HL3 , & PDF_A , & ! can remove these after development PDFITERS , & #ifdef PDFDIAG @@ -2190,7 +2349,7 @@ subroutine hystpdf( & PDF_RWQT, & #endif WTHV2, & - WQL, & + WQC, & needs_preexisting, & USE_BERGERON, & SC_ICE ) @@ -2199,14 +2358,14 @@ subroutine hystpdf( & integer, intent(in) :: PDFSHAPE real, intent(inout) :: TE,QV,QLLS,QILS,CLLS,QLCN,QICN,CLCN,PDF_A real, intent(in) :: NL,NI,CNVFRC,SRF_TYPE - real, intent(in) :: WHL,WQT,HL2,QT2,HLQT,W3,W2,MFQT3,MFHL3 + real, intent(in) :: WHL,WQT,HL2,QT2,HLQT,W3,W2bar,QT3,HL3 #ifdef PDFDIAG real, intent(out) :: PDF_SIGW1, PDF_SIGW2, PDF_W1, PDF_W2, & PDF_SIGHL1, PDF_SIGHL2, PDF_HL1, PDF_HL2, & PDF_SIGQT1, PDF_SIGQT2, PDF_QT1, PDF_QT2, & PDF_RHLQT, PDF_RWHL, PDF_RWQT #endif - real, intent(out) :: WTHV2, WQL + real, intent(out) :: WTHV2, WQC real, intent(out) :: PDFITERS logical, intent(in) :: needs_preexisting, USE_BERGERON real, optional , intent(in) :: SC_ICE @@ -2223,8 +2382,10 @@ subroutine hystpdf( & real :: QCx, QC, fQi, QCi, qsnx real :: dQICN, dQLCN, dQILS, dQLLS, Nfac, NLv, NIv - real :: tmpARR + real :: tmpARR, tmp, tmp2 real :: alhxbcp, DQCALL + + real :: exner, thv, bastoeps, beta, rwqt, rwhl, rhlqt, t1, t2, sigt1, sigt2, q1, q2, w1, w2, sigw1, sigw2 ! internal scalars integer :: N, nmax @@ -2243,8 +2404,34 @@ subroutine hystpdf( & DQS = GEOS_DQSAT( TEn, PL, QSAT=QSx ) QVn = ( QV - QSx*CLCN )*tmpARR - QT = QCn + QVn !Total LS water after microphysics - + QT = max(0.,QCn + QVn) !Total LS water after microphysics + ! - this can be negative because QV does not account + ! for saturation inside convective cloud fraction + + if (pdfshape .eq. 6) then + exner = (pl / 1e3)**MAPL_KAPPA + call precalc_dblgss( pdf_a, wqt, whl, TEn, hl2, hl3, qt2, qt3, w2bar, w3, & ! inputs + beta, rwqt, rwhl, rhlqt, t1, t2, sigt1, sigt2, & ! outputs + q1, q2, sigmaqt1, sigmaqt2, w1, w2, sigw1, sigw2) +#ifdef PDFDIAG + PDF_SIGW1 = sigw1 + PDF_SIGW2 = sigw2 + PDF_W1 = w1 + PDF_W2 = w2 + PDF_SIGHL1 = sigt1 + PDF_SIGHL2 = sigt2 + PDF_HL1 = te + gravbcp*ZL - alhxbcp*QCn + t1 + PDF_HL2 = te + gravbcp*ZL - alhxbcp*QCn + t2 + PDF_SIGQT1 = sigmaqt1 + PDF_SIGQT2 = sigmaqt2 + PDF_QT1 = qt + q1 + PDF_QT2 = qt + q2 + PDF_RHLQT = rhlqt + PDF_RWHL = rwhl + PDF_RWQT = rwqt +#endif + end if + nmax = 20 do n=1,nmax @@ -2283,11 +2470,14 @@ subroutine hystpdf( & fQi = ice_fraction( TEn, CNVFRC,SRF_TYPE ) alhxbcp = (1.0-fQi)*alhlbcp + fQi*alhsbcp HL = TEn + gravbcp*ZL - alhxbcp*QCn - + tmp = max( QT2, WQT*WQT/max(W2bar,1e-4) ) + tmp2 = max( HL2, WHL*WHL/max(W2bar,1e-4) ) + call partition_dblgss(fQi, & TEn, & QVn, & QCn, & + QSn, & 0.0, & ! assume OMEGA=0 ZL, & PL*100., & @@ -2295,13 +2485,13 @@ subroutine hystpdf( & HL, & WHL, & WQT, & - HL2, & - QT2, & + tmp2, & + tmp, & HLQT, & W3, & - W2, & - MFQT3, & - MFHL3, & + W2bar, & + QT3, & + HL3, & PDF_A, & #ifdef PDFDIAG PDF_SIGW1, & @@ -2321,9 +2511,27 @@ subroutine hystpdf( & PDF_RWQT, & #endif WTHV2, & - WQL, & + WQC, & CFn) - endif + + elseif (PDFSHAPE.eq.6) then + + if (qt+q2+2.*sigmaqt2.gt.qsn) then + ! Update the liquid water static energy + fQi = ice_fraction( TEn, CNVFRC,SRF_TYPE ) + alhxbcp = (1.0-fQi)*alhlbcp + fQi*alhsbcp + HL = TEn + gravbcp*ZL - alhxbcp*QCn + + call partition_dblgss2( exner, pdf_a, beta, rwqt, rwhl, rhlqt, t1, t2, & + sigt1, sigt2, qt, q1, q2, sigmaqt1, sigmaqt2, & + w1, w2, sigw1, sigw2, QSN, DQS, QCn, CFn, WQC ) + else + QCn = 0. + CFn = 0. + WQC = 0. + end if + + endif IF(USE_BERGERON) THEN DQCALL = QCn - QCp @@ -2359,8 +2567,8 @@ subroutine hystpdf( & ! This next line needs correcting - need proper d(del qc)/dT derivative for triangular ! for now, just use relaxation of 1/2 of top-hat. QCn = QCp + 0.5*(QCn-QCp)/(1.-(CFn*(ALPHA-1.)-(QCn/QSn))*DQS*alhxbcp) - elseif(PDFSHAPE.eq.5) then - QCn = QCp + 0.5*(QCn-QCp) + elseif(PDFSHAPE.ge.5) then + QCn = QCp + 0.75*(QCn-QCp) endif QVn = QVp - (QCn - QCp) @@ -2372,6 +2580,14 @@ subroutine hystpdf( & enddo ! qsat iteration + if (PDFSHAPE.eq.6) then ! Double Gaussian buoyancy flux calculation + thv = (TEn/exner) * (1.+QVn*MAPL_H2OMW/MAPL_AIRMW) + bastoeps = (MAPL_RVAP/MAPL_RGAS) * thv + + wthv2 = whl + (MAPL_H2OMW/MAPL_AIRMW)*thv*wqt & + + (alhxbcp-bastoeps)*wqc + end if + ! Now take {\em New} condensate and partition into ice and liquid ! large-scale diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/cloudnew.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/cloudnew.F90 index 3f2006297..5be88bfd8 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/cloudnew.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/cloudnew.F90 @@ -1996,6 +1996,7 @@ subroutine hystpdf_new( & TEn, & QVn, & QCn, & + qsn, & 0.0, & ! assume OMEGA=0 ZL, & PL*100., & diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSturbulence_GridComp/GEOS_TurbulenceGridComp.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSturbulence_GridComp/GEOS_TurbulenceGridComp.F90 index 5540383df..4c8c76d62 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSturbulence_GridComp/GEOS_TurbulenceGridComp.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSturbulence_GridComp/GEOS_TurbulenceGridComp.F90 @@ -282,7 +282,7 @@ subroutine SetServices ( GC, RC ) ! !IMPORT STATE: call MAPL_AddImportSpec(GC, & - LONG_NAME = 'surface geopotential height', & + LONG_NAME = 'surface geopotential', & UNITS = 'm+2 s-2', & SHORT_NAME = 'PHIS', & DIMS = MAPL_DimsHorzOnly, & @@ -480,46 +480,6 @@ subroutine SetServices ( GC, RC ) RC=STATUS ) VERIFY_(STATUS) -! call MAPL_AddImportSpec(GC, & -! SHORT_NAME = 'MFTHSRC', & -! LONG_NAME = 'mass_flux_source_temperature_perturbation', & -! UNITS = 'K', & -! DIMS = MAPL_DimsHorzVert, & -! VLOCATION = MAPL_VLocationCenter, & -! RESTART = MAPL_RestartSkip, & -! RC=STATUS ) -! VERIFY_(STATUS) - -! call MAPL_AddImportSpec(GC, & -! SHORT_NAME = 'MFQTSRC', & -! LONG_NAME = 'mass_flux_source_humidity_perturbation', & -! UNITS = 'kg kg-1', & -! DIMS = MAPL_DimsHorzVert, & -! VLOCATION = MAPL_VLocationCenter, & -! RESTART = MAPL_RestartSkip, & -! RC=STATUS ) -! VERIFY_(STATUS) - -! call MAPL_AddImportSpec(GC, & -! SHORT_NAME = 'MFW', & -! LONG_NAME = 'mass_flux_initial_vertical_velocity', & -! UNITS = 'm s-1', & -! DIMS = MAPL_DimsHorzVert, & -! VLOCATION = MAPL_VLocationCenter, & -! RESTART = MAPL_RestartSkip, & -! RC=STATUS ) -! VERIFY_(STATUS) - -! call MAPL_AddImportSpec(GC, & -! SHORT_NAME = 'MFAREA', & -! LONG_NAME = 'mass_flux_area_fraction', & -! UNITS = '1', & -! DIMS = MAPL_DimsHorzVert, & -! VLOCATION = MAPL_VLocationCenter, & -! RESTART = MAPL_RestartSkip, & -! RC=STATUS ) -! VERIFY_(STATUS) - call MAPL_AddImportSpec(GC, & SHORT_NAME = 'Z0', & LONG_NAME = 'surface_roughness', & @@ -923,24 +883,6 @@ subroutine SetServices ( GC, RC ) RC=STATUS ) VERIFY_(STATUS) -! call MAPL_AddExportSpec(GC, & -! LONG_NAME = 'EDMF_updraft_contribution_to_total_water_variance', & -! UNITS = 'kg2 kg-2', & -! SHORT_NAME = 'EDMF_QT2' , & -! DIMS = MAPL_DimsHorzVert, & -! VLOCATION = MAPL_VLocationCenter, & -! RC=STATUS ) -! VERIFY_(STATUS) - -! call MAPL_AddExportSpec(GC, & -! LONG_NAME = 'Liquid_static_energy_variance_diagnosed_from_updrafts', & -! UNITS = 'K2', & -! SHORT_NAME = 'EDMF_SL2' , & -! DIMS = MAPL_DimsHorzVert, & -! VLOCATION = MAPL_VLocationCenter, & -! RC=STATUS ) -! VERIFY_(STATUS) - call MAPL_AddExportSpec(GC, & LONG_NAME = 'Liquid_static_energy_flux_from_updrafts', & UNITS = 'K s-1', & @@ -1023,32 +965,32 @@ subroutine SetServices ( GC, RC ) VERIFY_(STATUS) - call MAPL_AddExportSpec(GC, & - LONG_NAME = 'Diagnostic_liquid_water_static_energy_variance', & - UNITS = 'K2' , & - SHORT_NAME = 'SL2DIAG' , & - DIMS = MAPL_DimsHorzVert, & - VLOCATION = MAPL_VLocationCenter, & - RC=STATUS ) - VERIFY_(STATUS) +! call MAPL_AddExportSpec(GC, & +! LONG_NAME = 'Diagnostic_liquid_water_static_energy_variance', & +! UNITS = 'K2' , & +! SHORT_NAME = 'SL2DIAG' , & +! DIMS = MAPL_DimsHorzVert, & +! VLOCATION = MAPL_VLocationCenter, & +! RC=STATUS ) +! VERIFY_(STATUS) - call MAPL_AddExportSpec(GC, & - LONG_NAME = 'Diagnostic_total_water_variance', & - UNITS = 'kg2 kg-2' , & - SHORT_NAME = 'QT2DIAG' , & - DIMS = MAPL_DimsHorzVert, & - VLOCATION = MAPL_VLocationCenter, & - RC=STATUS ) - VERIFY_(STATUS) +! call MAPL_AddExportSpec(GC, & +! LONG_NAME = 'Diagnostic_total_water_variance', & +! UNITS = 'kg2 kg-2' , & +! SHORT_NAME = 'QT2DIAG' , & +! DIMS = MAPL_DimsHorzVert, & +! VLOCATION = MAPL_VLocationCenter, & +! RC=STATUS ) +! VERIFY_(STATUS) - call MAPL_AddExportSpec(GC, & - LONG_NAME = 'Diagnostic_liquid_static_energy_total_water_covariance',& - UNITS = 'K kg kg-1' , & - SHORT_NAME = 'SLQTDIAG' , & - DIMS = MAPL_DimsHorzVert, & - VLOCATION = MAPL_VLocationCenter, & - RC=STATUS ) - VERIFY_(STATUS) +! call MAPL_AddExportSpec(GC, & +! LONG_NAME = 'Diagnostic_liquid_static_energy_total_water_covariance',& +! UNITS = 'K kg kg-1' , & +! SHORT_NAME = 'SLQTDIAG' , & +! DIMS = MAPL_DimsHorzVert, & +! VLOCATION = MAPL_VLocationCenter, & +! RC=STATUS ) +! VERIFY_(STATUS) call MAPL_AddExportSpec(GC, & LONG_NAME = 'Third_moment_of_liquid_water_static_energy', & @@ -1068,14 +1010,14 @@ subroutine SetServices ( GC, RC ) RC=STATUS ) VERIFY_(STATUS) - call MAPL_AddExportSpec(GC, & - LONG_NAME = 'Third_moment_of_vertical_velocity_Canuto_estimate', & - UNITS = 'm3 s-3', & - SHORT_NAME = 'W3CANUTO' , & - DIMS = MAPL_DimsHorzVert, & - VLOCATION = MAPL_VLocationCenter, & - RC=STATUS ) - VERIFY_(STATUS) +! call MAPL_AddExportSpec(GC, & +! LONG_NAME = 'Third_moment_of_vertical_velocity_Canuto_estimate', & +! UNITS = 'm3 s-3', & +! SHORT_NAME = 'W3CANUTO' , & +! DIMS = MAPL_DimsHorzVert, & +! VLOCATION = MAPL_VLocationCenter, & +! RC=STATUS ) +! VERIFY_(STATUS) call MAPL_AddExportSpec(GC, & LONG_NAME = 'Vertical_velocity_variance', & @@ -1167,33 +1109,6 @@ subroutine SetServices ( GC, RC ) RC=STATUS ) VERIFY_(STATUS) - call MAPL_AddExportSpec(GC, & - SHORT_NAME = 'SLFLXMF', & - LONG_NAME = 'liquid_water_static_energy_flux_by_MF', & - UNITS = 'K m s-1', & - DIMS = MAPL_DimsHorzVert, & - VLOCATION = MAPL_VLocationEdge, & - RC=STATUS ) - VERIFY_(STATUS) - - call MAPL_AddExportSpec(GC, & - SHORT_NAME = 'QTFLXMF', & - LONG_NAME = 'total_water_flux_by_MF', & - UNITS = 'kg kg-1 m s-1', & - DIMS = MAPL_DimsHorzVert, & - VLOCATION = MAPL_VLocationEdge, & - RC=STATUS ) - VERIFY_(STATUS) - - call MAPL_AddExportSpec(GC, & - SHORT_NAME = 'MFAW', & - LONG_NAME = 'EDMF_kinematic_mass_flux', & - UNITS = 'm s-1', & - DIMS = MAPL_DimsHorzVert, & - VLOCATION = MAPL_VLocationEdge, & - RC=STATUS ) - VERIFY_(STATUS) - call MAPL_AddExportSpec(GC, & SHORT_NAME = 'TRI', & LONG_NAME = 'diffusion_tendencies', & @@ -1224,41 +1139,41 @@ subroutine SetServices ( GC, RC ) RC=STATUS ) VERIFY_(STATUS) - call MAPL_AddExportSpec(GC, & - LONG_NAME = 'air_temperature', & - UNITS = 'K', & - SHORT_NAME = 'T', & - DIMS = MAPL_DimsHorzVert, & - VLOCATION = MAPL_VLocationCenter, & - RC=STATUS ) - VERIFY_(STATUS) +! call MAPL_AddExportSpec(GC, & +! LONG_NAME = 'air_temperature', & +! UNITS = 'K', & +! SHORT_NAME = 'T', & +! DIMS = MAPL_DimsHorzVert, & +! VLOCATION = MAPL_VLocationCenter, & +! RC=STATUS ) +! VERIFY_(STATUS) - call MAPL_AddExportSpec(GC, & - LONG_NAME = 'eastward_wind', & - UNITS = 'm s-1', & - SHORT_NAME = 'U', & - DIMS = MAPL_DimsHorzVert, & - VLOCATION = MAPL_VLocationCenter, & - RC=STATUS ) - VERIFY_(STATUS) +! call MAPL_AddExportSpec(GC, & +! LONG_NAME = 'eastward_wind', & +! UNITS = 'm s-1', & +! SHORT_NAME = 'U', & +! DIMS = MAPL_DimsHorzVert, & +! VLOCATION = MAPL_VLocationCenter, & +! RC=STATUS ) +! VERIFY_(STATUS) - call MAPL_AddExportSpec(GC, & - LONG_NAME = 'northward_wind', & - UNITS = 'm s-1', & - SHORT_NAME = 'V', & - DIMS = MAPL_DimsHorzVert, & - VLOCATION = MAPL_VLocationCenter, & - RC=STATUS ) - VERIFY_(STATUS) +! call MAPL_AddExportSpec(GC, & +! LONG_NAME = 'northward_wind', & +! UNITS = 'm s-1', & +! SHORT_NAME = 'V', & +! DIMS = MAPL_DimsHorzVert, & +! VLOCATION = MAPL_VLocationCenter, & +! RC=STATUS ) +! VERIFY_(STATUS) - call MAPL_AddExportSpec(GC, & - LONG_NAME = 'specific_humidity', & - UNITS = 'kg kg-1', & - SHORT_NAME = 'QV', & - DIMS = MAPL_DimsHorzVert, & - VLOCATION = MAPL_VLocationCenter, & - RC=STATUS ) - VERIFY_(STATUS) +! call MAPL_AddExportSpec(GC, & +! LONG_NAME = 'specific_humidity', & +! UNITS = 'kg kg-1', & +! SHORT_NAME = 'QV', & +! DIMS = MAPL_DimsHorzVert, & +! VLOCATION = MAPL_VLocationCenter, & +! RC=STATUS ) +! VERIFY_(STATUS) call MAPL_AddExportSpec(GC, & LONG_NAME = 'total_momentum_diffusivity', & @@ -1279,7 +1194,7 @@ subroutine SetServices ( GC, RC ) VERIFY_(STATUS) call MAPL_AddExportSpec(GC, & - LONG_NAME = 'Richardson_number_from_Louis', & + LONG_NAME = 'Richardson_number', & UNITS = '1', & SHORT_NAME = 'RI', & DIMS = MAPL_DimsHorzVert, & @@ -1360,7 +1275,7 @@ subroutine SetServices ( GC, RC ) VERIFY_(STATUS) call MAPL_AddExportSpec(GC, & - LONG_NAME = 'Blackadar_length_scale_for_heat', & + LONG_NAME = 'Blackadar_length_scale_for_scalars', & UNITS = 'm', & SHORT_NAME = 'ALH', & DIMS = MAPL_DimsHorzVert, & @@ -1998,14 +1913,6 @@ subroutine SetServices ( GC, RC ) VLOCATION = MAPL_VLocationCenter, RC=STATUS ) VERIFY_(STATUS) - call MAPL_AddExportSpec(GC, & - SHORT_NAME = 'LMIX', & - LONG_NAME = 'mixed_layer_depth_from_SHOC', & - UNITS = 'm', & - DIMS = MAPL_DimsHorzOnly, & - VLOCATION = MAPL_VLocationNone, RC=STATUS ) - VERIFY_(STATUS) - call MAPL_AddExportSpec(GC, & SHORT_NAME = 'LSHOC1', & LONG_NAME = 'dissipation_length_term1_from_SHOC', & @@ -2510,6 +2417,33 @@ subroutine SetServices ( GC, RC ) VLOCATION = MAPL_VLocationCenter, RC=STATUS ) VERIFY_(STATUS) + call MAPL_AddInternalSpec(GC, & + SHORT_NAME = 'MFAW', & + LONG_NAME = 'updraft_normalized_mass_flux', & + UNITS = 'm s-1', & + DEFAULT = 0.0, & + DIMS = MAPL_DimsHorzVert, & + VLOCATION = MAPL_VLocationEdge, RC=STATUS ) + VERIFY_(STATUS) + + call MAPL_AddInternalSpec(GC, & + SHORT_NAME = 'SLFLXMF', & + LONG_NAME = 'liquid_water_static_energy_flux_by_updrafts',& + UNITS = 'K m s-1', & + DEFAULT = 0.0, & + DIMS = MAPL_DimsHorzVert, & + VLOCATION = MAPL_VLocationEdge, RC=STATUS ) + VERIFY_(STATUS) + + call MAPL_AddInternalSpec(GC, & + SHORT_NAME = 'QTFLXMF', & + LONG_NAME = 'total_water_flux_by_updrafts', & + UNITS = 'kg kg-1 m s-1', & + DEFAULT = 0.0, & + DIMS = MAPL_DimsHorzVert, & + VLOCATION = MAPL_VLocationEdge, RC=STATUS ) + VERIFY_(STATUS) + !EOS ! Set the Profiling timers @@ -2649,8 +2583,9 @@ subroutine RUN1 ( GC, IMPORT, EXPORT, CLOCK, RC ) ! SHOC-related variables integer :: DO_SHOC, SCM_SL - real, dimension(:,:,:), pointer :: TKESHOC,TKH,QT2,QT3,WTHV2,WQT_DC,PDF_A - + real, dimension(:,:,:), pointer :: TKESHOC,TKH,QT2,QT3,WTHV2,WQT_DC, & + PDF_A,MFAW,QTFLXMF,SLFLXMF + real, dimension(:,:), pointer :: EVAP, SH ! Idealized SCM surface layer variables @@ -2816,7 +2751,12 @@ subroutine RUN1 ( GC, IMPORT, EXPORT, CLOCK, RC ) VERIFY_(STATUS) call MAPL_GetPointer(INTERNAL, PDF_A, 'PDF_A', RC=STATUS) VERIFY_(STATUS) - + call MAPL_GetPointer(INTERNAL, MFAW, 'MFAW', RC=STATUS) + VERIFY_(STATUS) + call MAPL_GetPointer(INTERNAL, QTFLXMF, 'QTFLXMF', RC=STATUS) + VERIFY_(STATUS) + call MAPL_GetPointer(INTERNAL, SLFLXMF, 'SLFLXMF', RC=STATUS) + VERIFY_(STATUS) ! ! edmf variables ! @@ -2991,7 +2931,6 @@ subroutine REFRESH(IM,JM,LM,RC) 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 :: Z0, Z0H, EIS real, dimension(:,: ), pointer :: BSTAR, USTAR, PPBL, WERAD, WESFC,VSCRAD,KERAD,DBUOY,ZSML,ZCLD,ZRADML,FRLAND,TRINVBS,TRINVFRQ,TRINVDELT @@ -3020,8 +2959,8 @@ subroutine REFRESH(IM,JM,LM,RC) LSHOC1,LSHOC2,LSHOC3, & SHOCPRNUM,& TKEBUOY,TKESHEAR,TKEDISS,TKEDISSx, & - SL2, SL3, W2, W3, WSL, SLQT, W3CANUTO, QT2DIAG,SL2DIAG,SLQTDIAG - real, dimension(:,:), pointer :: LMIX, edmf_depth + SL2, SL3, W2, W3, WSL, SLQT !, W3CANUTO, QT2DIAG,SL2DIAG,SLQTDIAG + real, dimension(:,:), pointer :: edmf_depth ! EDMF variables real, dimension(:,:,:), pointer :: edmf_dry_a,edmf_moist_a,edmf_frc, edmf_dry_w,edmf_moist_w, & @@ -3033,12 +2972,10 @@ subroutine REFRESH(IM,JM,LM,RC) edmf_w2, & !edmf_qt2, edmf_sl2, & edmf_w3, edmf_wqt, edmf_slqt, & edmf_wsl, edmf_qt3, edmf_sl3, & - edmf_entx, edmf_tke, slflxmf, & - qtflxmf, mfaw, edmf_dqrdt, edmf_dqsdt, & - ssrcmf,qvsrcmf,qlsrcmf,qisrcmf + edmf_entx, edmf_tke, & + edmf_dqrdt, edmf_dqsdt real, dimension(IM,JM,0:LM) :: ae3,aw3,aws3,awqv3,awql3,awqi3,awu3,awv3 - real, dimension(IM,JM,1:LM) :: ssrc,qvsrc,qlsrc,qisrc real, dimension(IM,JM) :: zpbl_test @@ -3083,7 +3020,7 @@ subroutine REFRESH(IM,JM,LM,RC) ! all fluxes from surface layer theory ! else: use prescribed surface exchange coefficients real :: SCM_SH ! prescribed surface sensible heat flux (Wm-1) (for SCM_SL_FLUX == 1) - real :: SCM_EVAP ! prescribed surface latent heat flux (Wm-1) (for SCM_SL_FLUX == 1) + real :: SCM_LH ! prescribed surface latent heat flux (Wm-1) (for SCM_SL_FLUX == 1) real :: SCM_Z0 ! surface roughness length (m) real :: SCM_ZETA ! Monin-Obkhov length scale (m) (for SCM_SL_FLUX == 3) real :: SCM_RH_SURF ! Surface relative humidity @@ -3096,7 +3033,7 @@ subroutine REFRESH(IM,JM,LM,RC) real, dimension(:,:), pointer :: SHOBS, LHOBS ! mass-flux constants/parameters - integer :: DOMF, NumUp, DOCLASP + integer :: DOMF, NumUp real :: L0,L0fac real, dimension(IM,JM) :: L02 @@ -3117,7 +3054,7 @@ subroutine REFRESH(IM,JM,LM,RC) real, dimension(im,jm,0:lm) :: RHOE, RHOAW3, edmf_mf, mfwsl, mfwqt, mftke real, dimension(im,jm,lm) :: buoyf, mfw2, mfw3, mfqt3, & mfsl3, mfqt2, mfsl2, & - mfslqt, edmf_ent !mfwhl, edmf_ent + mfslqt, edmf_ent real :: a1,a2 real, dimension(IM,JM,LM) :: dum3d,tmp3d,WVP @@ -3128,9 +3065,9 @@ subroutine REFRESH(IM,JM,LM,RC) ! variables associated with SHOC real, dimension( IM, JM, LM ) :: QPL,QPI - integer :: DO_SHOC, DOPROGQT2, DOCANUTO + integer :: DO_SHOC, DOPROGQT2 real :: SL2TUNE, QT2TUNE, SLQT2TUNE, & - SKEW_TGEN, SKEW_TDIS + SKEW_TGEN, SKEW_TDIS, FREE_ATM_QT2 real :: PDFSHAPE real :: lambdadiss @@ -3197,11 +3134,6 @@ subroutine REFRESH(IM,JM,LM,RC) call MAPL_GetPointer(IMPORT, Z0H, 'Z0H', RC=STATUS); VERIFY_(STATUS) call MAPL_GetPointer(IMPORT, EIS, 'EIS', RC=STATUS); VERIFY_(STATUS) - ! Imports for CLASP heterogeneity coupling in EDMF -! call MAPL_GetPointer(IMPORT, MFTHSRC, 'MFTHSRC',RC=STATUS); VERIFY_(STATUS) -! call MAPL_GetPointer(IMPORT, MFQTSRC, 'MFQTSRC',RC=STATUS); VERIFY_(STATUS) -! call MAPL_GetPointer(IMPORT, MFW, 'MFW' ,RC=STATUS); VERIFY_(STATUS) -! call MAPL_GetPointer(IMPORT, MFAREA, 'MFAREA' ,RC=STATUS); VERIFY_(STATUS) ! Get turbulence parameters from configuration !--------------------------------------------- @@ -3302,32 +3234,34 @@ 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.25, RC=STATUS); VERIFY_(STATUS) + call MAPL_GetResource (MAPL, SHOCPARAMS%PRNUM, trim(COMP_NAME)//"_SHC_PRNUM:", default=-0.9, RC=STATUS); VERIFY_(STATUS) + call MAPL_GetResource (MAPL, SHOCPARAMS%LAMBDA, trim(COMP_NAME)//"_SHC_LAMBDA:", default=0.5, 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=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=1.0, RC=STATUS); VERIFY_(STATUS) + call MAPL_GetResource (MAPL, SHOCPARAMS%LENFAC1, trim(COMP_NAME)//"_SHC_LENFAC1:", default=8., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetResource (MAPL, SHOCPARAMS%LENFAC2, trim(COMP_NAME)//"_SHC_LENFAC2:", default=2., RC=STATUS); VERIFY_(STATUS) + call MAPL_GetResource (MAPL, SHOCPARAMS%LENFAC3, trim(COMP_NAME)//"_SHC_LENFAC3:", default=1., RC=STATUS); VERIFY_(STATUS) call MAPL_GetResource (MAPL, SHOCPARAMS%BUOYOPT, trim(COMP_NAME)//"_SHC_BUOY_OPTION:", default=2, RC=STATUS); VERIFY_(STATUS) + call MAPL_GetResource (MAPL, PDFSHAPE, 'PDFSHAPE:', DEFAULT = 6.0 , RC=STATUS); VERIFY_(STATUS) + else + call MAPL_GetResource (MAPL, PDFSHAPE, 'PDFSHAPE:', DEFAULT = 1.0 , 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 = 1.0 , RC=STATUS); VERIFY_(STATUS) call MAPL_GetResource (MAPL, SLQT2TUNE, 'SLQT2TUNE:', DEFAULT = 7.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_TDIS, 'SKEW_TDIS:', DEFAULT = 900.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) + call MAPL_GetResource (MAPL, FREE_ATM_QT2, 'FREE_ATM_QT2:', DEFAULT = 0.05, RC=STATUS); VERIFY_(STATUS) ! Get pointers from export state... !----------------------------------- - PDFALLOC = (PDFSHAPE.eq.5) + PDFALLOC = (PDFSHAPE.ge.5) call MAPL_GetPointer(EXPORT, KH, 'KH', ALLOC=.TRUE., RC=STATUS) VERIFY_(STATUS) @@ -3473,8 +3407,8 @@ subroutine REFRESH(IM,JM,LM,RC) VERIFY_(STATUS) call MAPL_GetPointer(EXPORT, w3, 'W3', ALLOC=PDFALLOC, RC=STATUS) VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT, w3canuto,'W3CANUTO', ALLOC=PDFALLOC, RC=STATUS) - VERIFY_(STATUS) +! call MAPL_GetPointer(EXPORT, w3canuto,'W3CANUTO', ALLOC=PDFALLOC, RC=STATUS) +! VERIFY_(STATUS) call MAPL_GetPointer(EXPORT, w2, 'W2', ALLOC=PDFALLOC, RC=STATUS) VERIFY_(STATUS) call MAPL_GetPointer(EXPORT, sl3, 'SL3', ALLOC=PDFALLOC, RC=STATUS) @@ -3485,12 +3419,12 @@ subroutine REFRESH(IM,JM,LM,RC) ! VERIFY_(STATUS) call MAPL_GetPointer(EXPORT, wsl, 'WSL', ALLOC=PDFALLOC, RC=STATUS) VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT, qt2diag, 'QT2DIAG', ALLOC=PDFALLOC, RC=STATUS) - VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT, sl2diag, 'SL2DIAG', ALLOC=PDFALLOC, RC=STATUS) - VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT, slqtdiag, 'SLQTDIAG', ALLOC=PDFALLOC, RC=STATUS) - VERIFY_(STATUS) +! call MAPL_GetPointer(EXPORT, qt2diag, 'QT2DIAG', ALLOC=PDFALLOC, RC=STATUS) +! VERIFY_(STATUS) +! call MAPL_GetPointer(EXPORT, sl2diag, 'SL2DIAG', ALLOC=PDFALLOC, RC=STATUS) +! VERIFY_(STATUS) +! call MAPL_GetPointer(EXPORT, slqtdiag, 'SLQTDIAG', ALLOC=PDFALLOC, RC=STATUS) +! VERIFY_(STATUS) call MAPL_GetPointer(EXPORT, edmf_wqt, 'EDMF_WQT', ALLOC=PDFALLOC, RC=STATUS) VERIFY_(STATUS) call MAPL_GetPointer(EXPORT, edmf_wsl, 'EDMF_WSL', ALLOC=PDFALLOC, RC=STATUS) @@ -3499,14 +3433,6 @@ subroutine REFRESH(IM,JM,LM,RC) VERIFY_(STATUS) call MAPL_GetPointer(EXPORT, edmf_mfx, 'EDMF_MF', ALLOC=PDFALLOC, 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, qisrcmf, 'QISRCMF', 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) @@ -3539,12 +3465,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, mfaw, 'MFAW', RC=STATUS) - VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT, slflxmf, 'SLFLXMF', RC=STATUS) - VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT, qtflxmf, 'QTFLXMF', RC=STATUS) - VERIFY_(STATUS) !========== SHOC =========== call MAPL_GetPointer(EXPORT, TKEDISSx, 'TKEDISS', RC=STATUS) @@ -3559,8 +3479,6 @@ subroutine REFRESH(IM,JM,LM,RC) VERIFY_(STATUS) call MAPL_GetPointer(EXPORT, LSHOC1, 'LSHOC1', RC=STATUS) VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT, LMIX, 'LMIX', RC=STATUS) - VERIFY_(STATUS) call MAPL_GetPointer(EXPORT, LSHOC2, 'LSHOC2', RC=STATUS) VERIFY_(STATUS) call MAPL_GetPointer(EXPORT, LSHOC3, 'LSHOC3', RC=STATUS) @@ -3695,20 +3613,19 @@ subroutine REFRESH(IM,JM,LM,RC) ! get updraft constants call MAPL_GetResource (MAPL, DOMF, "EDMF_DOMF:", default=0, RC=STATUS) - MFPARAMS%DOTRACERS = .false. if ( DOMF /= 0 ) then ! number of updrafts call MAPL_GetResource (MAPL, MFPARAMS%NUP, "EDMF_NUMUP:", default=10, RC=STATUS) - call MAPL_GetResource (MAPL, MFPARAMS%DOTRACERS, "EDMF_DOTRACERS:", default=.false., RC=STATUS) - + call MAPL_GetResource (MAPL, MFPARAMS%DOTRACERS, "EDMF_DOTRACERS:", default=.true., 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%PWMAX, "EDMF_PWMAX:", default=3., 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) + call MAPL_GetResource (MAPL, MFPARAMS%WC, "EDMF_WC:", default=0.016, RC=STATUS) + call MAPL_GetResource (MAPL, MFPARAMS%WCTHRESH, "EDMF_WCTHRESH:", default=11., 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) @@ -3716,16 +3633,12 @@ subroutine REFRESH(IM,JM,LM,RC) ! 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.4, RC=STATUS) - call MAPL_GetResource (MAPL, MFPARAMS%ENT0LTS, "EDMF_ENT0LTS:", default=0.8, RC=STATUS) + call MAPL_GetResource (MAPL, MFPARAMS%ENT0, "EDMF_ENT0:", default=0.35, RC=STATUS) ! L0 if ET==1 call MAPL_GetResource (MAPL, MFPARAMS%L0, "EDMF_L0:", default=100., RC=STATUS) ! L0fac if ET==2 call MAPL_GetResource (MAPL, MFPARAMS%L0fac, "EDMF_L0FAC:", default=10., RC=STATUS) - call MAPL_GetResource (MAPL, MFPARAMS%MFLIMFAC, "EDMF_MFLIMFAC:", default=2.5, RC=STATUS) - ! factor to multiply the eddy-diffusivity with - call MAPL_GetResource (MAPL, MFPARAMS%EDfac, "EDMF_EDFAC:", default=1., RC=STATUS) - call MAPL_GetResource (MAPL, MFPARAMS%DOCLASP, "EDMF_DOCLASP:", default=0, RC=STATUS) + call MAPL_GetResource (MAPL, MFPARAMS%MFLIMFAC, "EDMF_MFLIMFAC:", default=2.0, RC=STATUS) call MAPL_GetResource (MAPL, MFPARAMS%ICE_RAMP, "EDMF_ICE_RAMP:", default=-40.0, RC=STATUS ) call MAPL_GetResource (MAPL, MFPARAMS%ENTRAIN, "EDMF_ENTRAIN:", default=0, RC=STATUS) call MAPL_GetResource (MAPL, MFPARAMS%STOCHFRAC, "EDMF_STOCHASTIC:", default=0.5, RC=STATUS) @@ -3733,20 +3646,10 @@ subroutine REFRESH(IM,JM,LM,RC) call MAPL_GetResource (MAPL, MFPARAMS%IMPLICIT, "EDMF_IMPLICIT:", default=1, RC=STATUS) call MAPL_GetResource (MAPL, MFPARAMS%PRCPCRIT, "EDMF_PRCPCRIT:", default=-1., RC=STATUS) call MAPL_GetResource (MAPL, MFPARAMS%UPABUOYDEP,"EDMF_UPABUOYDEP:", default=1, RC=STATUS) - - ! Future options -! call MAPL_GetResource (MAPL, EDMF_THERMAL_PLUME, "EDMF_THERMAL_PLUME:", default=0, RC=STATUS) -! 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_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) -! 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) + MFPARAMS%TREFF = 0. + MFPARAMS%DOTRACERS = .false. end if call MAPL_GetResource(MAPL, SCM_SL, 'SCM_SL:', DEFAULT=0 ) @@ -3758,7 +3661,7 @@ subroutine REFRESH(IM,JM,LM,RC) call MAPL_GetResource(MAPL, SCM_SL_FLUX, 'SCM_SL_FLUX:', DEFAULT=0 ) call MAPL_GetResource(MAPL, SCM_SH, 'SCM_SH:', DEFAULT=0. ) - call MAPL_GetResource(MAPL, SCM_EVAP, 'SCM_EVAP:', DEFAULT=0. ) + call MAPL_GetResource(MAPL, SCM_LH, 'SCM_LH:', DEFAULT=0. ) call MAPL_GetResource(MAPL, SCM_Z0, 'SCM_Z0:', DEFAULT=1.E-4 ) call MAPL_GetResource(MAPL, SCM_RH_SURF, 'SCM_RH_SURF:', DEFAULT=0.98 ) call MAPL_GetResource(MAPL, SCM_TSURF, 'SCM_TSURF:', DEFAULT=298.76 ) ! S6 @@ -3778,7 +3681,7 @@ subroutine REFRESH(IM,JM,LM,RC) if ( SCM_SL_FLUX == 1 ) then sh_scm(:,:) = scm_sh - evap_scm(:,:) = scm_evap/MAPL_ALHL + evap_scm(:,:) = scm_lh/MAPL_ALHL elseif ( SCM_SL_FLUX == 2 ) then sh_scm(:,:) = shobs evap_scm(:,:) = lhobs/MAPL_ALHL @@ -3814,10 +3717,8 @@ subroutine REFRESH(IM,JM,LM,RC) bstar => bstar_scm ustar => ustar_scm - - print *,'bstar=',bstar_scm,' ustar=',ustar_scm - +! print *,'bstar=',bstar_scm,' ustar=',ustar_scm call MAPL_TimerOff(MAPL,"---SURFACE") end if @@ -3830,7 +3731,7 @@ subroutine REFRESH(IM,JM,LM,RC) !=============================================================== call MAPL_TimerOn(MAPL,"---MASSFLUX") -! Initialize EDMF output variables needed for update_moments +! Initialize EDMF output variables mfsl2 = 0.0 mfslqt = 0.0 mfqt2 = 0.0 @@ -3841,11 +3742,15 @@ subroutine REFRESH(IM,JM,LM,RC) mfwqt = 0.0 mfwsl = 0.0 mftke = 0.0 - ssrc = 0.0 - qvsrc = 0.0 - qlsrc = 0.0 - qisrc = 0.0 - + YS = 0.0 + YQV = 0.0 + YQL = 0.0 + YQI = 0.0 + YU = 0.0 + YV = 0.0 + awqv3 = 0.0 + awql3 = 0.0 + awqi3 = 0.0 IF(DOMF /= 0) then @@ -3869,7 +3774,6 @@ subroutine REFRESH(IM,JM,LM,RC) EVAP, & FRLAND, & ZPBL, & -! MFTHSRC, MFQTSRC, MFW, MFAREA, & ! CLASP inputs !== Outputs for trisolver == ae3, & aw3, & @@ -3879,10 +3783,12 @@ subroutine REFRESH(IM,JM,LM,RC) awqi3, & awu3, & awv3, & - ssrc, & - qvsrc, & - qlsrc, & - qisrc, & + YS, & + YQV, & + YQL, & + YQI, & + YU, & + YV, & !== Outputs for ADG PDF == mfw2, & mfw3, & @@ -3915,18 +3821,15 @@ subroutine REFRESH(IM,JM,LM,RC) EDMF_PLUMES_THL, & EDMF_PLUMES_QT ) + mfaw = aw3 + slflxmf = (aws3-awql3*mapl_alhl-awqi3*mapl_alhs)/mapl_cp + qtflxmf = awqv3+awql3+awqi3 + !=== 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 = 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(qisrcmf)) qisrcmf = qisrc if (associated(edmf_w2)) edmf_w2 = mfw2 if (associated(edmf_w3)) edmf_w3 = mfw3 if (associated(edmf_qt3)) edmf_qt3 = mfqt3 @@ -3938,13 +3841,13 @@ subroutine REFRESH(IM,JM,LM,RC) 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 + 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 @@ -3957,10 +3860,13 @@ subroutine REFRESH(IM,JM,LM,RC) awu3 = 0.0 awv3 = 0.0 buoyf = 0.0 + mfaw = 0.0 + slflxmf = 0.0 + qtflxmf = 0.0 if (associated(edmf_dry_a)) edmf_dry_a = 0.0 if (associated(edmf_moist_a)) edmf_moist_a = 0.0 -! if (associated(edmf_dry_w)) edmf_dry_w = MAPL_UNDEF + if (associated(edmf_dry_w)) edmf_dry_w = MAPL_UNDEF if (associated(edmf_moist_w)) edmf_moist_w = MAPL_UNDEF if (associated(edmf_dry_qt)) edmf_dry_qt = MAPL_UNDEF if (associated(edmf_moist_qt)) edmf_moist_qt = MAPL_UNDEF @@ -3974,13 +3880,6 @@ subroutine REFRESH(IM,JM,LM,RC) if (associated(edmf_buoyf)) edmf_buoyf = 0.0 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(qisrcmf)) qisrcmf = 0.0 - if (associated(qvsrcmf)) qvsrcmf = 0.0 - if (associated(slflxmf)) slflxmf = 0.0 - if (associated(qtflxmf)) qtflxmf = 0.0 if (associated(edmf_w2)) edmf_w2 = mfw2 if (associated(edmf_w3)) edmf_w3 = mfw3 if (associated(edmf_qt3)) edmf_qt3 = mfqt3 @@ -3994,16 +3893,14 @@ subroutine REFRESH(IM,JM,LM,RC) drycblh = 0. ENDIF + + + call MAPL_TimerOff(MAPL,"---MASSFLUX") !!!================================================================= !!!=========================== SHOC ============================== -!!!================================================================= -! Description -! -! -! !!!================================================================= if ( DO_SHOC /= 0 ) then @@ -4016,6 +3913,7 @@ subroutine REFRESH(IM,JM,LM,RC) call RUN_SHOC( IM, JM, LM, LM+1, DT, & !== Inputs == + SH(:,:), & PLO(:,:,1:LM), & ZL0(:,:,0:LM), & Z(:,:,1:LM), & @@ -4044,7 +3942,6 @@ subroutine REFRESH(IM,JM,LM,RC) TKEBUOY, & TKESHEAR, & LSHOC, & - LMIX, & LSHOC1, & LSHOC2, & LSHOC3, & @@ -4557,6 +4454,8 @@ subroutine REFRESH(IM,JM,LM,RC) call update_moments(IM, JM, LM, DT, & SH, & ! in EVAP, & + AREA, & + ZPBL, & Z, & ZLE, & KH, & @@ -4566,8 +4465,6 @@ subroutine REFRESH(IM,JM,LM,RC) QT, & SL, & EDMF_FRC, & -! edmf_mf(:,:,1:LM)/rhoe(:,:,1:LM), & -! MFQT2, & MFQT3, & MFSL3, & MFW2, & @@ -4583,19 +4480,14 @@ subroutine REFRESH(IM,JM,LM,RC) sl3, & w2, & w3, & - w3canuto, & - wsl, & slqt, & - qt2diag, & - sl2diag, & - slqtdiag, & doprogqt2, & ! tuning parameters sl2tune, & qt2tune, & slqt2tune, & skew_tgen, & skew_tdis, & - docanuto ) + free_atm_qt2 ) end if @@ -4916,7 +4808,7 @@ subroutine REFRESH(IM,JM,LM,RC) ZPBL = MIN(ZPBL,Z(:,:,KPBLMIN)) KPBL = MAX(KPBL,float(KPBLMIN)) - ! Calc KPBL using surface turbulence, for use in shallow scheme + ! Calc KPBL using surface turbulence, for use in UW shallow scheme if (associated(KPBL_SC)) then KPBL_SC = MAPL_UNDEF do I = 1, IM @@ -5068,26 +4960,7 @@ subroutine REFRESH(IM,JM,LM,RC) ! Y-s ... these are rhs - mean value - surface flux ! (these are added in the diffuse and vrtisolve) - - -! -! 2:LM -> 1:LM-1, 1:LM-1 -> 0:LM-2 -! - 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) + 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) + QISRC(:,:,1:LM-1) ) - - YU(:,:,1:LM-1) = DMI(:,:,1:LM-1)*( RHOE(:,:,1:LM-1)*AWU3(:,:,1:LM-1) - RHOE(:,:,0:LM-2)*AWU3(:,:,0:LM-2) ) - YV(:,:,1:LM-1) = DMI(:,:,1:LM-1)*( RHOE(:,:,1:LM-1)*AWV3(:,:,1:LM-1) - RHOE(:,:,0:LM-2)*AWV3(:,:,0:LM-2) ) - + ! Add prescribed surface fluxes if ( SCM_SL /= 0 .and. (SCM_SL_FLUX == 1 .or. SCM_SL_FLUX == 2) ) then YS(:,:,LM) = YS(:,:,LM) + DMI(:,:,LM)*SH(:,:)!/RHOE(:,:,LM) @@ -5174,29 +5047,6 @@ subroutine REFRESH(IM,JM,LM,RC) AKV = AKX BKV = BKX - ! - ! LU decomposition for the mass-flux variables - ! - AKX=AKSS - BKX=BKSS - call VTRILU(AKX,BKX,CKSS) - BKSS=BKX - AKSS=AKX - - AKX=AKQQ - BKX=BKQQ - call VTRILU(AKX,BKX,CKQQ) - BKQQ=BKX - AKQQ=AKX - - AKX=AKUU - BKX=BKUU - call VTRILU(AKX,BKX,CKUU) - BKUU=BKX - AKUU=AKX - - - ! Get the sensitivity of solution to a unit ! change in the surface value. B and C are ! not modified. @@ -5206,6 +5056,27 @@ subroutine REFRESH(IM,JM,LM,RC) call VTRISOLVESURF(BKQ,CKQ,DKQ) call VTRISOLVESURF(BKV,CKV,DKV) + ! + ! LU decomposition for the mass-flux variables + ! + AKX=AKSS + BKX=BKSS + call VTRILU(AKX,BKX,CKSS) + BKSS=BKX + AKSS=AKX + + AKX=AKQQ + BKX=BKQQ + call VTRILU(AKX,BKX,CKQQ) + BKQQ=BKX + AKQQ=AKX + + AKX=AKUU + BKX=BKUU + call VTRILU(AKX,BKX,CKUU) + BKUU=BKX + AKUU=AKX + call MAPL_TimerOff(MAPL,"---DECOMP") if(ALLOC_TCZPBL) deallocate(TCZPBL) @@ -5274,7 +5145,7 @@ subroutine DIFFUSE(IM,JM,LM,RC) ! Parameters for idealized SCM surface layer integer :: SCM_SL, SCM_SL_FLUX - real :: SCM_SH, SCM_EVAP + real :: SCM_SH, SCM_LH ! EDMF transport real :: EntExp, EntDyn @@ -5308,7 +5179,7 @@ subroutine DIFFUSE(IM,JM,LM,RC) VERIFY_(STATUS) call MAPL_GetResource(MAPL, SCM_SH, 'SCM_SH:', default=0., RC=STATUS) VERIFY_(STATUS) - call MAPL_GetResource(MAPL, SCM_EVAP, 'SCM_EVAP:', default=0., RC=STATUS) + call MAPL_GetResource(MAPL, SCM_LH, 'SCM_LH:', default=0., RC=STATUS) VERIFY_(STATUS) CU => cu_scm @@ -5468,6 +5339,8 @@ subroutine DIFFUSE(IM,JM,LM,RC) end if end if + + ! Pick the right exchange coefficients !------------------------------------- @@ -5541,13 +5414,13 @@ subroutine DIFFUSE(IM,JM,LM,RC) DX => DKQ AK => AKQQ; BK => BKQQ; CK => CKQQ SX=S+YQL -! OPT = .FALSE. + if (DO_SHOC/=0.) OPT = .FALSE. elseif (trim(name)=='QILS') then CX => CQ DX => DKQ AK => AKQQ; BK => BKQQ; CK => CKQQ SX=S+YQI -! OPT = .FALSE. + if (DO_SHOC/=0.) OPT = .FALSE. elseif (trim(name)=='U') then CX => CU DX => DKV @@ -5574,7 +5447,7 @@ subroutine DIFFUSE(IM,JM,LM,RC) if ( trim(name) == 'S' ) then SF(:,:) = scm_sh elseif ( trim(name) == 'Q' ) then - SF(:,:) = scm_evap/mapl_alhl + SF(:,:) = scm_lh/mapl_alhl end if else if ( SCM_SL /= 0 .and. SCM_SL_FLUX ==2 ) then if ( trim(name) == 'S' ) then @@ -5815,8 +5688,8 @@ subroutine UPDATE(IM,JM,LM,LATS,RC) real, dimension(:,: ), pointer :: KETRB, KESRF, KETOP, KEINT 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, & + real, dimension(:,:,:), pointer :: QTFLXTRB, SLFLXTRB, WSL, WQT, & + TKH, UFLXTRB, VFLXTRB, QTX, SLX, & SLFLXMF, QTFLXMF, MFAW, TKEDISS integer :: KM, K, L, I, J @@ -5873,10 +5746,10 @@ subroutine UPDATE(IM,JM,LM,LATS,RC) call MAPL_GetPointer(EXPORT, VFLXTRB , 'VFLXTRB' , RC=STATUS); VERIFY_(STATUS) ! MF contribution, used to calculate TRB fluxes above - call MAPL_GetPointer(EXPORT, SLFLXMF , 'SLFLXMF' , ALLOC=.TRUE., RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT, QTFLXMF , 'QTFLXMF' , ALLOC=.TRUE., RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT, MFAW , 'MFAW' , ALLOC=.TRUE., RC=STATUS); VERIFY_(STATUS) - + call MAPL_GetPointer(INTERNAL, SLFLXMF , 'SLFLXMF' , RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(INTERNAL, QTFLXMF , 'QTFLXMF' , RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(INTERNAL, MFAW , 'MFAW' , RC=STATUS); VERIFY_(STATUS) + ! Used in update_moments for ADG PDF (requires all of above) call MAPL_GetPointer(EXPORT, WSL, 'WSL' , ALLOC=.TRUE., RC=STATUS); VERIFY_(STATUS) call MAPL_GetPointer(EXPORT, WQT, 'WQT' , ALLOC=.TRUE., RC=STATUS); VERIFY_(STATUS) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSturbulence_GridComp/edmf.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSturbulence_GridComp/edmf.F90 index eae0b4a6c..41ff47479 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSturbulence_GridComp/edmf.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSturbulence_GridComp/edmf.F90 @@ -8,7 +8,8 @@ module edmf_mod use MAPL_ConstantsMod, only: mapl_epsilon, mapl_grav, mapl_cp, & mapl_alhl, mapl_p00, mapl_vireps, & - mapl_alhs, mapl_kappa, mapl_rgas + mapl_alhs, mapl_alhf, mapl_kappa, & + mapl_pi, mapl_celsius_to_kelvin use MAPL_Mod, only: mapl_undef @@ -16,31 +17,19 @@ module edmf_mod implicit none -real, parameter :: & - WSTARmin = 1.e-3, & - zpblmin = 100., & - onethird = 1./3., & - r = 2. - - type EDMFPARAMS_TYPE +type EDMFPARAMS_TYPE logical :: DOTRACERS integer :: DISCRETE integer :: IMPLICIT integer :: ENTRAIN - integer :: DOCLASP integer :: NUP - integer :: THERMAL_PLUME - integer :: TEST - integer :: DEBUG integer :: ET integer :: UPABUOYDEP real :: L0 real :: L0fac real :: STOCHFRAC real :: ENTUFAC - real :: EDFAC real :: ENT0 - real :: ENT0LTS real :: ALPHATH real :: ALPHAQT real :: ALPHAW @@ -48,88 +37,85 @@ module edmf_mod real :: PWMIN real :: WA real :: WB - real :: AU0 - real :: CTH1 - real :: CTH2 - real :: RH0_QB - real :: C_KH_MF + real :: WC + real :: WCTHRESH real :: MFLIMFAC real :: ICE_RAMP real :: PRCPCRIT real :: TREFF - endtype EDMFPARAMS_TYPE - type (EDMFPARAMS_TYPE) :: MFPARAMS +endtype EDMFPARAMS_TYPE +type (EDMFPARAMS_TYPE) :: MFPARAMS public run_edmf, mfparams contains - -SUBROUTINE RUN_EDMF(its,ite, jts,jte, kts,kte, dt, & ! Inputs - phis, & - zlo3, & - zw3, & - pw3, & - rhoe3, & - tke3, & - u3, & - v3, & - t3, & - thl3, & - thv3, & - qv3, & - ql3, & - qi3, & - wthl2, & - wqt2, & - frland, & - pblh2, & - ! CLASP inputs for surface heterogeneity - ! mfsrcthl, mfsrcqt, mfw, mfarea, & - ! outputs - variables needed for solver - ae3, & - aw3, & - aws3, & - awqv3, & - awql3, & - awqi3, & - awu3, & - awv3, & - ssrc3, & - qvsrc3, & - qlsrc3, & - qisrc3, & - ! Outputs required for SHOC and ADG PDF - mfw2, & - mfw3, & - mfqt3, & - mfhl3, & - mfwqt, & - mfhlqt, & - mfwhl, & - mftke, & - buoyf, & - edmfmf, & - dry_a3, & - moist_a3, & - dqrdt, & - dqsdt, & - ! Diagnostic outputs - updraft properties - dry_w3, & - moist_w3, & - dry_qt3, & - moist_qt3, & - dry_thl3, & - moist_thl3, & - dry_u3, & - moist_u3, & - dry_v3, & - moist_v3, & - moist_qc3, & - entx, & - mfdepth, & - edmf_plumes_w, & - edmf_plumes_thl, & - edmf_plumes_qt ) + !==== Inputs =================================================== +SUBROUTINE RUN_EDMF(its,ite, jts,jte, kts,kte, dt, & ! Index limits and timestep (s) + phis, & ! Surface geopotential (m2 s-2) + zlo3, & ! Surface-relative heights (m) + zw3, & ! Surface-relative edge heights (m) + pw3, & ! Edge pressures (Pa) + rhoe3, & ! Edge densities (kg m-3) + tke3, & ! Turbulent kinetic energy (m2 s-2) + u3, & ! U wind component (m s-1) + v3, & ! V wind component (m s-1) + t3, & ! Temperature (K) + thl3, & ! Liquid water potential temperature (K) + thv3, & ! Virtual potential temperature (K) + qv3, & ! Specific humidity (kg kg-1) + ql3, & ! Specific large-scale cloud liquid (kg kg-1) + qi3, & ! Specific large-scale cloud ice (kg kg-1) + wthl2, & ! Surface sensible heat flux (W m-2) + wqt2, & ! Surface evaporation (kg m-2 s-1) + frland, & ! Land fraction + pblh2, & ! PBL height (m) + !==== Outputs - variables needed for solver ==================== + ae3, & ! Environmental area fraction + aw3, & ! Area-weighted updraft vertical velocity (m s-1) + aws3, & ! Dry static energy flux (J m s-1 kg-1) + awqv3, & ! Specific humidity flux (kg kg-1 m s-1) + awql3, & ! Specific liquid flux (kg kg-1 m s-1) + awqi3, & ! Specific ice flux (kg kg-1 m s-1) + awu3, & ! Kinematic U momentum flux (m2 s-2) + awv3, & ! Kinematic V momentum flux (m2 s-2) + YS, & ! Dry static energy increment for trisolver (J kg-1) + YQV, & ! Specific humidity increment for trisolver (kg kg-1) + YQL, & ! Liquid water increment for trisolver (kg kg-1) + YQI, & ! Ice water increment (kg kg-1) + YU, & ! U wind increment (m s-1) + YV, & ! V wind increment (m s-1) + !==== Outputs required for SHOC and ADG PDF ==================== + mfw2, & ! Area-weighted vertical velocity squared (m2 s-2) + mfw3, & ! Area-weighted vertical velocity cubed (m3 s-3) + mfqt3, & ! Area-weighted updraft total water anomaly cubed (kg3 kg-3) + mfhl3, & ! Area-weighted updraft liquid static energy anomaly cubed (K3) + mfwqt, & ! Kinematic total water flux on midlevels (kg kg-1 m s-1) + mfhlqt, & ! Updraft total water-temperature anomaly covariance (kg kg-1 K) + mfwhl, & ! Kinematic static energy flux on midlevels (J kg-1 m s-1) + mftke, & ! Updraft kinetic energy (m2 s-2) + buoyf, & ! Updraft buoyancy flux (K m s-1) + edmfmf, & ! Updraft mass flux (kg m s-1) + dry_a3, & ! Dry updraft area fraction + moist_a3, & ! Cloudy updraft area fraction + dqrdt, & ! Liquid precipitation tendency + dqsdt, & ! Frozen precipitation tendency + !==== Diagnostic outputs - updraft properties ================== + dry_w3, & ! Dry updraft mean vertical velocity (m s-1) + moist_w3, & ! Cloudy updraft mean vertical velocity (m s-1) + dry_qt3, & ! Dry updraft mean total water (kg kg-1) + moist_qt3, & ! Cloudy updraft mean total water (kg kg-1) + dry_thl3, & ! Dry updraft liquid water potential temperature (K) + moist_thl3, & ! Cloudy updraft liquid water potential temperature (K) + dry_u3, & ! Dry updraft mean u wind (m s-1) + moist_u3, & ! Cloudy updraft mean u wind (m s-1) + dry_v3, & ! Dry updraft mean v wind (m s-1) + moist_v3, & ! Cloudy updraft mean v wind (m s-1) + moist_qc3, & ! Cloudy updraft mean total condensate (kg kg-1) + entx, & ! Mean fractional lateral mixing rate + mfdepth, & ! Diagnostic updraft depth used in lateral mixing (m-1) + edmf_plumes_w, & ! Individual updraft vertical velocities (m s-1) + edmf_plumes_thl, & ! Individual updraft liquid water potential temperatures (K) + edmf_plumes_qt ) ! Individual updraft total waters (kg kg-1) INTEGER, INTENT(IN) :: ITS,ITE,JTS,JTE,KTS,KTE @@ -154,9 +140,6 @@ SUBROUTINE RUN_EDMF(its,ite, jts,jte, kts,kte, dt, & ! Inputs FRLAND, & PHIS - ! CLASP inputs - REAL,DIMENSION(ITS:ITE,JTS:JTE,KTS:KTE) :: mfsrcqt,mfsrcthl,mfw,mfarea - ! Required outputs REAL,DIMENSION(ITS:ITE,JTS:JTE,KTS-1:KTE), INTENT(OUT) :: dry_a3, & moist_a3, & @@ -173,18 +156,22 @@ SUBROUTINE RUN_EDMF(its,ite, jts,jte, kts,kte, dt, & ! Inputs mfwqt, & mftke - REAL,DIMENSION(ITS:ITE,JTS:JTE,KTS:KTE), INTENT(OUT) :: buoyf,mfw2,mfw3,mfqt3,mfhl3,& - mfhlqt,dqrdt,dqsdt,ssrc3,qvsrc3,qlsrc3,qisrc3 + REAL,DIMENSION(ITS:ITE,JTS:JTE,KTS:KTE), INTENT(OUT) :: buoyf, mfw2, mfw3, & + mfqt3, mfhl3, mfhlqt, & + dqrdt, dqsdt + REAL,DIMENSION(ITS:ITE,JTS:JTE,KTS:KTE), INTENT(INOUT) :: YS, YQV, YQL, & + YQI, YU, YV ! Diagnostic outputs - REAL, DIMENSION(:,:), POINTER :: mfdepth - REAL, DIMENSION(:,:,:), POINTER :: dry_w3, moist_w3, & - dry_qt3, moist_qt3, & - dry_thl3, moist_thl3, & - dry_u3, moist_u3, & - dry_v3, moist_v3, & - moist_qc3, entx + REAL, DIMENSION(:,:), POINTER :: mfdepth + + REAL, DIMENSION(:,:,:), POINTER :: dry_w3, moist_w3, & + dry_qt3, moist_qt3, & + dry_thl3, moist_thl3, & + dry_u3, moist_u3, & + dry_v3, moist_v3, & + moist_qc3, entx REAL, DIMENSION(:,:,:,:), POINTER :: EDMF_PLUMES_W, & EDMF_PLUMES_THL, & @@ -194,54 +181,53 @@ SUBROUTINE RUN_EDMF(its,ite, jts,jte, kts,kte, dt, & ! Inputs !============= Local variables ============= ! updraft properties - REAL,DIMENSION(KTS-1:KTE,1:MFPARAMS%NUP) :: UPW,UPTHL,UPQT,UPQL,UPQI, & - UPA,UPU,UPV,UPTHV + REAL,DIMENSION(KTS-1:KTE,1:MFPARAMS%NUP) :: UPW, UPTHL, UPQT, & + UPQL, UPQI, UPA, & + UPU, UPV, UPTHV ! entrainment variables - REAl,DIMENSION(KTS:KTE,1:MFPARAMS%NUP) :: ENT,ENTf + REAl,DIMENSION(KTS:KTE,1:MFPARAMS%NUP) :: ENT, ENTf INTEGER,DIMENSION(KTS:KTE,1:MFPARAMS%NUP) :: ENTi - INTEGER :: K,I,IH,JH,NUP2 - REAL :: wthv,wstar,qstar,thstar,sigmaW,sigmaQT,sigmaTH,z0, & - wmin,wmax,wlv,wtv,wp - REAL :: B,QTn,THLn,THVn,QCn,QP,Un,Vn,Wn2,EntEXP,EntEXPU,EntW,wf - -! internal flipped variables (GEOS5) - REAL,DIMENSION(KTS:KTE) :: U,V,THL,QT,THV,QV,QL,QI,ZLO,QR,QS - REAL,DIMENSION(KTS-1:KTE) :: ZW,P,THLI,QTI + REAL,DIMENSION(ITS:ITE,JTS:JTE,KTS:KTE) :: tmp3d, dmi + + INTEGER :: K,KTMP,I,IH,JH,NUP2 + REAL :: wthv,wstar,qstar,thstar, & + sigmaW,sigmaQT,sigmaTH, & + wmin,wmax,wlv,wtv,wp, & + B,QTn,THLn,THVn,QCn,QP, & + Un,Vn,Wn2,EntEXP,EntEXPU,& + EntW,wf, WTHL, WQT, PBLH + + ! internal flipped variables (GEOS) + REAL,DIMENSION(KTS:KTE) :: U,V,THL,QT,THV,QV,QL,QI,ZLO,QR,QS + REAL,DIMENSION(KTS-1:KTE) :: ZW,P,THLI,QTI REAL,DIMENSION(KTS-1:KTE) :: UI, VI, QVI, QLI, QII -! internal surface cont - REAL :: UST,WTHL,WQT,PBLH + ! Updraft diagnostics REAL,DIMENSION(KTS-1:KTE) :: dry_a, moist_a,dry_w,moist_w, & dry_qt,moist_qt,dry_thl,moist_thl, & dry_u,moist_u,dry_v,moist_v, moist_qc + REAL,DIMENSION(KTS-1:KTE) :: s_aw,s_aws,s_awqv,s_awql,s_awqi,s_awu,s_awv REAL,DIMENSION(KTS:KTE) :: s_buoyf REAL,DIMENSION(KTS-1:KTE) :: s_aw2,s_aw3,s_aqt3,s_ahl3,s_aqt2, & s_ahlqt,s_awqt,s_ahl2,s_awhl,qte -! exner function - REAL,DIMENSION(KTS:KTE) :: exf,dp,pmid - REAL,DIMENSION(KTS-1:KTE) :: exfh,tmp1d - REAL,DIMENSION(KTS-1:KTE) :: rhoe + REAL,DIMENSION(KTS:KTE) :: exf,dp,pmid,wcfac + REAL,DIMENSION(KTS-1:KTE) :: exfh,rhoe - REAL :: L0,ztop,tmp,tmp2,ltm,MFsrf,QTsrfF,THVsrfF,mft,mfthvt,mf,factor,lts - INTEGER, DIMENSION(2) :: seedmf,the_seed + ! temporary/dummy variables + REAL :: tmp, tmp2 + + REAL :: L0,ztop,ltm,QTsrfF,THVsrfF,mft,mfthvt,mf,factor + INTEGER, DIMENSION(2) :: the_seed LOGICAL :: calc_avg_diag ! min values to avoid singularities - REAL,PARAMETER :: & - WSTARmin=1.e-3, & - PBLHmin=100. - - ! temporary, set - mfsrcthl = -999. - mfsrcqt = -999. - mfw = -999. - mfarea = -999. - - tmp1d(:) = 1e-3 + REAL, PARAMETER :: & + WSTARmin = 1.e-3, & + PBLHmin = 100. ! If any average diagnostics requested, perform required calculations, ! otherwise skip for efficiency @@ -277,10 +263,6 @@ SUBROUTINE RUN_EDMF(its,ite, jts,jte, kts,kte, dt, & ! Inputs awqi3 =0. awu3 =0. awv3 =0. - ssrc3 =0. - qvsrc3=0. - qlsrc3=0. - qisrc3=0. buoyf =0. mfw2 =0. mfw3 =0. @@ -293,13 +275,11 @@ SUBROUTINE RUN_EDMF(its,ite, jts,jte, kts,kte, dt, & ! Inputs edmfmf=0. dqrdt =0. dqsdt =0. - ! mfqt2 =0. - ! mfhl2 =0. if (associated(entx)) entx = mapl_undef - ! this is the environmental area - by default 1. - ae3=MFPARAMS%EDfac + ! Initialize the environmental area. Updraft area will be subtracted below. + ae3=1. DO IH=ITS,ITE ! loop over the horizontal dimensions @@ -312,17 +292,22 @@ SUBROUTINE RUN_EDMF(its,ite, jts,jte, kts,kte, dt, & ! Inputs pblh=max(pblh,pblhmin) wthv=wthl+mapl_epsilon*thv3(IH,JH,kte)*wqt - ! if CLASP enabled: mass flux is input - ! if CLASP disabled: mass-flux if positive surface buoyancy flux and - ! TKE at 2nd model level above threshold - 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 + ! calc average TKE below 100 m + tmp = 0. + tmp2 = 0. + k = kte + do while (zlo3(IH,JH,k)<100. .and. k>1) + tmp2 = tmp2 + (zw3(IH,JH,k-1)-zw3(IH,JH,k)) + tmp = tmp+tke3(IH,JH,k)*(zw3(IH,JH,k-1)-zw3(IH,JH,k)) + k = k-1 + end do + tmp = tmp/tmp2 ! avg TKE - if (MFPARAMS%doclasp/=0) then - nup2 = count(mfsrcthl(IH,JH,1:MFPARAMS%NUP)>=-2.0) - else - nup2 = MFPARAMS%NUP - end if + ! Activate mass-flux only if positive surface buoyancy flux + ! and mean TKE below 100m exceeds threshold (for stability) + IF (wthv > 0.0 .and. tmp>0.05 .and. phis(IH,JH).lt.3e4) then + + nup2 = MFPARAMS%NUP UPW=0. UPTHL=0. @@ -339,27 +324,16 @@ SUBROUTINE RUN_EDMF(its,ite, jts,jte, kts,kte, dt, & ! Inputs ! Estimate scale height for entrainment calculation if (mfparams%ET == 2 ) then - pmid = 0.5*(pw3(IH,JH,kts-1:kte-1)+pw3(IH,JH,kts:kte)) - call calc_mf_depth(kts,kte,t3(IH,JH,:),zlo3(IH,JH,:)-zw3(IH,JH,kte),qv3(IH,JH,:),pmid,ztop,wthv,wqt) - L0 = max(min(ztop,2500.),500.) / mfparams%L0fac - if (associated(mfdepth)) mfdepth(IH,JH) = ztop - - ! Reduce L0 over ocean where LTS is large to encourage StCu formation - lts = 0.0 - if (FRLAND(IH,JH)<0.5) then - do k = kte-1,kts+1,-1 - if (zlo3(IH,JH,k)-zw3(IH,JH,kte).gt.3000.0) then - lts = thv3(IH,JH,k+1) - exit - 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-18.5))) ) - end if + pmid = 0.5*(pw3(IH,JH,kts-1:kte-1)+pw3(IH,JH,kts:kte)) + call calc_mf_depth(kts,kte,t3(IH,JH,:),zlo3(IH,JH,:)-zw3(IH,JH,kte),qv3(IH,JH,:),pmid,ztop,wthv,wqt) + L0 = max(min(ztop,2500.),500.) / mfparams%L0fac + if (associated(mfdepth)) mfdepth(IH,JH) = ztop else ! if mfparams%ET not 2 - L0 = mfparams%L0 + L0 = mfparams%L0 end if + if (ztop.gt.100.) then + ! ! flipping variables ! @@ -471,12 +445,27 @@ SUBROUTINE RUN_EDMF(its,ite, jts,jte, kts,kte, dt, & ! Inputs sigmaQT=MFPARAMS%AlphaQT*qstar sigmaTH=MFPARAMS%AlphaTH*thstar - if (MFPARAMS%doclasp/= 0) then - wmin=2.*sigmaW - wmax=2.*sigmaW - else - wmin=sigmaW*MFPARAMS%pwmin - wmax=sigmaW*MFPARAMS%pwmax + wmin=sigmaW*MFPARAMS%pwmin + wmax=sigmaW*MFPARAMS%pwmax + + ! Identify inversions below 1.5km, calculate stability in overlying 1km to define + ! a dynamic pressure deceleration factor in the updraft w equation below. + wcfac = 0. + tmp = 0. + k = kts+1 + do while (zlo(k).lt.1500.) + if ( t3(IH,JH,kte-k).gt.t3(IH,JH,kte-k+1) ) then + tmp = thv(k) ! THV at inversion + exit + end if + k = k+1 + end do + if (tmp.ne.0.) then + ktmp = k + do while (zlo(ktmp).lt.zlo(k)+1e3) + ktmp = ktmp+1 + end do + wcfac(1:k) = min(10.,max(0.,thv(ktmp)-thv(k)-MFPARAMS%WCTHRESH))*exp(-(zlo(k)-zlo(1:k))/100. ) end if ! define surface conditions @@ -485,28 +474,18 @@ SUBROUTINE RUN_EDMF(its,ite, jts,jte, kts,kte, dt, & ! Inputs wlv=wmin+(wmax-wmin)/(real(NUP2))*(real(i)-1.) wtv=wmin+(wmax-wmin)/(real(NUP2))*real(i) - if (MFPARAMS%doclasp/=0) then - UPW(kts-1,I) = MFW(IH,JH,I) - UPA(kts-1,I)=MFAREA(IH,JH,I) !0.5*(ERF(3.0/sqrt(2.))-ERF(1.0/sqrt(2.)))/real(NUP) ! assume equal size for now + UPW(kts-1,I)=min(0.5*(wlv+wtv), 5.) + if (MFPARAMS%UPABUOYDEP/=0) then + UPA(kts-1,I)=(0.5+0.5*TANH((wthv-0.02)/0.09))*(0.5*ERF(wtv/(sqrt(2.)*sigmaW))-0.5*ERF(wlv/(sqrt(2.)*sigmaW))) else - UPW(kts-1,I)=min(0.5*(wlv+wtv), 5.) - if (MFPARAMS%UPABUOYDEP/=0) then - UPA(kts-1,I)=(0.5+0.5*TANH((wthv-0.02)/0.09))*(0.5*ERF(wtv/(sqrt(2.)*sigmaW))-0.5*ERF(wlv/(sqrt(2.)*sigmaW))) - else - UPA(kts-1,I)=(0.5*ERF(wtv/(sqrt(2.)*sigmaW))-0.5*ERF(wlv/(sqrt(2.)*sigmaW))) - end if + UPA(kts-1,I)=(0.5*ERF(wtv/(sqrt(2.)*sigmaW))-0.5*ERF(wlv/(sqrt(2.)*sigmaW))) end if UPU(kts-1,I)=U(kts) UPV(kts-1,I)=V(kts) - if (MFPARAMS%doclasp/=0) then ! if CLASP, use tile-based perturbations - UPQT(kts-1,I)=QT(kts)+MFSRCQT(IH,JH,I) - UPTHV(kts-1,I)=THV(kts)+MFSRCTHL(IH,JH,I) - else - UPQT(kts-1,I)=QT(kts)+0.32*UPW(kts-1,I)*sigmaQT/sigmaW - UPTHV(kts-1,I)=THV(kts)+0.58*UPW(kts-1,I)*sigmaTH/sigmaW - end if + UPQT(kts-1,I)=QT(kts)+0.32*UPW(kts-1,I)*sigmaQT/sigmaW + UPTHV(kts-1,I)=THV(kts)+0.58*UPW(kts-1,I)*sigmaTH/sigmaW ENDDO ! NUP @@ -562,11 +541,11 @@ SUBROUTINE RUN_EDMF(its,ite, jts,jte, kts,kte, dt, & ! Inputs DO I=1,NUP2 ! loop over updrafts if (UPW(K-1,I).gt.0.) then - if (MFPARAMS%ENTRAIN==3) then ! dynamic entrainment rates - ENT(K,I) = MFPARAMS%ENT0*max(1e-4,B)/max(0.1,UPW(K,I)**2) - elseif (MFPARAMS%ENTRAIN==4) then - ENT(K,I) = ENT(K,I)*(1.+3.0*sqrt(TKE3(IH,JH,I))/max(0.1,UPW(K-1,I))) - end if +! if (MFPARAMS%ENTRAIN==3) then ! dynamic entrainment rates +! ENT(K,I) = MFPARAMS%ENT0*max(1e-4,B)/max(0.1,UPW(K,I)**2) +! elseif (MFPARAMS%ENTRAIN==4) then +! ENT(K,I) = ENT(K,I)*(1.+3.0*sqrt(TKE3(IH,JH,I))/max(0.1,UPW(K-1,I))) +! end if EntExp = exp(-ENT(K,I)*(ZW(k)-ZW(k-1))) EntExpU = exp(-ENT(K,I)*(ZW(k)-ZW(k-1))*MFPARAMS%EntUFac) @@ -592,7 +571,8 @@ SUBROUTINE RUN_EDMF(its,ite, jts,jte, kts,kte, dt, & ! Inputs ! vertical velocity B=mapl_grav*(0.5*(THVn+UPTHV(k-1,I))/THV(k)-1.) - WP=MFPARAMS%WB*ENT(K,I) + ! represent deceleration from dynamic pressure approaching inversion + WP=MFPARAMS%WB*ENT(K,I)+MFPARAMS%WC*wcfac(k) IF (WP==0.) THEN Wn2=UPW(K-1,I)**2+2.*MFPARAMS%WA*B*(ZW(k)-ZW(k-1)) ELSE @@ -611,20 +591,6 @@ 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))) ! qlflx divergence - 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 - tmp2 = max(0.,UPA(K,I)*(RHOE(K)*UPW(K,I)*UPQI(K,I)-RHOE(K-1)*UPW(K-1,I)*UPQI(K-1,I))) ! qiflx divergence - qisrc3(IH,JH,KTE+KTS-K) = qisrc3(IH,JH,KTE+KTS-K) + tmp2 - tmp = max(0.,UPA(K,I)*(RHOE(K-1)*UPW(K-1,I)*UPQL(K-1,I)-RHOE(K)*UPW(K,I)*UPQL(K,I))) ! qlflx convergence - ! if ql convergence, assume ice came from ql, with remainder from qv - qlsrc3(IH,JH,KTE+KTS-K) = qlsrc3(IH,JH,KTE+KTS-K) - min(tmp,tmp2) - qvsrc3(IH,JH,KTE+KTS-K) = qvsrc3(IH,JH,KTE+KTS-K) - (tmp2-min(tmp,tmp2)) - ssrc3(IH,JH,KTE+KTS-K) = ssrc3(IH,JH,KTE+KTS-K) + tmp2*MAPL_ALHS ELSE UPW(K,I) = 0. UPA(K,I) = 0. @@ -663,30 +629,24 @@ SUBROUTINE RUN_EDMF(its,ite, jts,jte, kts,kte, dt, & ! Inputs UPA = factor*UPA QR = factor*QR QS = factor*QS - ssrc3 = factor*ssrc3 - qvsrc3 = factor*qvsrc3 - qlsrc3 = factor*qlsrc3 - qisrc3 = factor*qisrc3 ! Rescale UPA if MF TKE more than half of prognostic TKE near surface ! Prevents instability due to MF without KH K = KTS tmp = 0. tmp2 = 0. - 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) + factor = 1. + DO WHILE (ZW(K).lt.200. .and. SUM(UPW(K-1,:)).gt.0.) + tmp = 0.25*SUM(UPA(K,:)*UPW(K,:)*UPW(K,:)+UPA(K-1,:)*UPW(K-1,:)*UPW(K-1,:)) + tmp2 = TKE3(IH,JH,KTE-K+KTS) + factor = min(factor,0.5*tmp2/tmp) K = K+1 END DO - if (tmp.gt.0.5*tmp2) then - factor = 0.5*tmp2/tmp + if (factor.lt.1.) then +! factor = 0.5*tmp2/tmp UPA = factor*UPA QR = factor*QR QS = factor*QS - ssrc3 = factor*ssrc3 - qvsrc3 = factor*qvsrc3 - qlsrc3 = factor*qlsrc3 - qisrc3 = factor*qisrc3 end if DO k=KTS,KTE @@ -865,7 +825,8 @@ SUBROUTINE RUN_EDMF(its,ite, jts,jte, kts,kte, dt, & ! Inputs if (associated(moist_qc3)) moist_qc3(IH,JH,KTS-1:KTE) = moist_qc(KTE:KTS-1:-1) - ! Note values were initialized to zero above + ! Note values were initialized to zero above. + ! Ending loop at interface above lowest layer. DO K=KTS-1,KTE-1 ! outputs - variables needed for solver aw3(IH,JH,K) = s_aw(KTE+KTS-K-1) @@ -875,14 +836,17 @@ SUBROUTINE RUN_EDMF(its,ite, jts,jte, kts,kte, dt, & ! Inputs awqi3(IH,JH,K) = s_awqi(KTE+KTS-K-1) awu3(IH,JH,K) = s_awu(KTE+KTS-K-1) awv3(IH,JH,K) = s_awv(KTE+KTS-K-1) - ae3(IH,JH,K) = (1.-dry_a(KTE+KTS-K-1)-moist_a(KTE+KTS-K-1))*MFPARAMS%EDfac + ae3(IH,JH,K) = (1.-dry_a(KTE+KTS-K-1)-moist_a(KTE+KTS-K-1)) mfwhl(IH,JH,K) = s_awhl(KTE+KTS-K-1) mfwqt(IH,JH,K) = s_awqt(KTE+KTS-K-1) ENDDO - mfwhl(IH,JH,KTE) = s_awhl(KTS-1) - mfwqt(IH,JH,KTE) = s_awqt(KTS-1) +! mfwhl(IH,JH,KTE) = s_awhl(KTS-1) +! mfwqt(IH,JH,KTE) = s_awqt(KTS-1) + + s_awqv(KTS-1) = 0. + s_awql(KTS-1) = 0. + s_awqi(KTS-1) = 0. - ! buoyancy is defined on full levels DO k=kts,kte @@ -897,8 +861,7 @@ SUBROUTINE RUN_EDMF(its,ite, jts,jte, kts,kte, dt, & ! Inputs 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 - ENDDO - + ENDDO where (UPA.eq.0.) UPW = MAPL_UNDEF @@ -909,19 +872,67 @@ SUBROUTINE RUN_EDMF(its,ite, jts,jte, kts,kte, dt, & ! Inputs if (associated(EDMF_PLUMES_THL)) EDMF_PLUMES_THL(IH,JH,KTS-1:KTE,:) = upthl(KTE:KTS-1:-1,:) if (associated(EDMF_PLUMES_QT)) EDMF_PLUMES_QT(IH,JH,KTS-1:KTE,:) = upqt(KTE:KTS-1:-1,:) - + end if ! IF ( mfdepth>100m) END IF ! IF ( wthv > 0.0 ) - ENDDO ! JH loop over horizontal area + ENDDO ! JH loop over horizontal area ENDDO ! IH + ! Calculate explicit part of state update for trisolver + ! Note vertical index is flipped, as in TurbGridComp. + dmi = mapl_grav * dt / ( pw3(:,:,kts:kte)-pw3(:,:,kts-1:kte-1) ) + YS(:,:,kte) = -rhoe3(:,:,kte-1) * aws3(:,:,kte-1) + YQV(:,:,kte) = -rhoe3(:,:,kte-1) * awqv3(:,:,kte-1) + YQL(:,:,kte) = -rhoe3(:,:,kte-1) * awql3(:,:,kte-1) + YQI(:,:,kte) = -rhoe3(:,:,kte-1) * awqi3(:,:,kte-1) + YU(:,:,kte) = -rhoe3(:,:,kte-1) * awu3(:,:,kte-1) + YV(:,:,kte) = -rhoe3(:,:,kte-1) * awv3(:,:,kte-1) + + YS(:,:,kts:kte-1) = ( rhoe3(:,:,kts:kte-1) *aws3(:,:,kts:kte-1) & + -rhoe3(:,:,kts-1:kte-2)*aws3(:,:,kts-1:kte-2) ) + YQV(:,:,kts:kte-1) = ( rhoe3(:,:,kts:kte-1) *awqv3(:,:,kts:kte-1) & + -rhoe3(:,:,kts-1:kte-2)*awqv3(:,:,kts-1:kte-2) ) + YQL(:,:,kts:kte-1) = ( rhoe3(:,:,kts:kte-1) *awql3(:,:,kts:kte-1) & + -rhoe3(:,:,kts-1:kte-2)*awql3(:,:,kts-1:kte-2) ) + YQI(:,:,kts:kte-1) = ( rhoe3(:,:,kts:kte-1) *awqi3(:,:,kts:kte-1) & + -rhoe3(:,:,kts-1:kte-2)*awqi3(:,:,kts-1:kte-2) ) + YU(:,:,kts:kte-1) = ( rhoe3(:,:,kts:kte-1) *awu3(:,:,kts:kte-1) & + -rhoe3(:,:,kts-1:kte-2)*awu3(:,:,kts-1:kte-2) ) + YV(:,:,kts:kte-1) = ( rhoe3(:,:,kts:kte-1) *awv3(:,:,kts:kte-1) & + -rhoe3(:,:,kts-1:kte-2)*awv3(:,:,kts-1:kte-2) ) + + ! Deal with implied condensation. + ! Where QI flux diverges and QL flux converges, assume that QI came from QL + where (YQI.lt.0. .and. YQL.gt.0.) + tmp3d = min(YQL,-YQI) + YQL = YQL - tmp3d + YQI = YQI + tmp3d + YS = YS + tmp3d*MAPL_ALHF + end where + where (YQI.lt.0.) ! where WQI diverges and no WQL + YQV = YQV - YQI + YS = YS - YQI*MAPL_ALHS + YQI = 0. + end where + where (YQL.lt.0.) ! where WQL diverges, assume condensation occurred + YQV = YQV + YQL + YS = YS - YQL*MAPL_ALHL + YQL = 0. + end where + + YS = dmi * YS + YQV = dmi * YQV + YQL = dmi * YQL + YQI = dmi * YQI + YU = dmi * YU + YV = dmi * YV + END SUBROUTINE run_edmf subroutine calc_mf_depth(kts,kte,t,z,q,p,ztop,wthv,wqt) integer, intent(in ) :: kts, kte -! real, intent(in ) :: ent0 real, intent(in ), dimension(kts:kte) :: t, z, q, p real, intent(in ) :: wthv, wqt real, intent( out) :: ztop @@ -929,16 +940,13 @@ subroutine calc_mf_depth(kts,kte,t,z,q,p,ztop,wthv,wqt) real :: tep,z1,z2,t1,t2,qp,pp,qsp,dqp,dqsp,wstar,qstar,thstar,sigmaQT,sigmaTH integer :: k - wstar=max(0.1,(mapl_grav/300.*wthv*1e3)**(1./3.)) ! convective velocity scale + wstar=max(0.1,(mapl_grav*wthv*1e3/t(kte))**(1./3.)) ! convective velocity scale qstar=max(0.,wqt)/wstar thstar=max(0.,wthv)/wstar -! sigmaW=MFPARAMS%AlphaW*wstar sigmaQT=2.0*qstar sigmaTH=2.0*thstar -! print *,'sigQT=',sigmaQT,' sigTH=',sigmaTH,' wstar=',wstar - tep = t(kte)+max(0.1,sigmaTH) ! parcel values qp = q(kte)+sigmaQT @@ -953,12 +961,9 @@ 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.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. ) dqp = max( qp - qsp, 0. )/(1.+(MAPL_ALHL/MAPL_CP)*dqsp ) @@ -1018,7 +1023,6 @@ subroutine condensation_edmf(QT,THL,P,THV,QC,wf,ice_ramp) QS=geos_qsat(T,P,pascals=.true.,ramp=ice_ramp) QC=max(QT-QS,0.) THV=(THL+get_alhl(T,ice_ramp)/mapl_cp*QC/EXN)*(1.+MAPL_VIREPS*(QT-QC)-QC) -!THV=(THL+get_alhl(T,ice_ramp)/mapl_cp*QC/EXN)*(1.+(mapl_epsilon)*(QT-QC)-QC) wf=water_f(T,ice_ramp) end subroutine condensation_edmf @@ -1099,7 +1103,7 @@ function water_f(T,iceramp) Tmax=0. Tmin=Tmax-abs(iceramp) - Tw=T-273.16 + Tw=T-mapl_celsius_to_kelvin ! water fraction IF (Tw>Tmax) THEN @@ -1122,7 +1126,7 @@ subroutine Poisson(istart,iend,jstart,jend,mu,POI,seed) integer,dimension(2), intent(in) :: seed integer :: seed_len -integer :: IH,JH,idum,p +integer :: IH,JH,p integer,allocatable :: theseed(:) call random_seed(SIZE=seed_len) @@ -1138,7 +1142,6 @@ subroutine Poisson(istart,iend,jstart,jend,mu,POI,seed) do ih=istart,iend do jh=jstart,jend -! poi(IH,JH)=poidev(mu(IH,JH),idum) poi(IH,JH)=poidev(mu(IH,JH)) enddo enddo @@ -1150,8 +1153,8 @@ end subroutine Poisson ! FUNCTION poidev(xm,idum) FUNCTION poidev(xm) INTEGER idum - REAL poidev,xm,PI - PARAMETER (PI=3.141592654) + REAL poidev,xm!,PI +! PARAMETER (PI=3.141592654) !CU USES gammln,ran1 REAL alxm,em,g,oldm,sq,t,y SAVE alxm,g,oldm,sq @@ -1175,7 +1178,7 @@ FUNCTION poidev(xm) g=xm*alxm-gammln(xm+1.) endif !1 y=tan(PI*ran1(idum)) -1 y=tan(PI*ran1()) +1 y=tan(MAPL_PI*ran1()) em=sq*y+xm if (em.lt.0.) goto 1 em=int(em) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSturbulence_GridComp/shoc.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSturbulence_GridComp/shoc.F90 index 64c696d3d..61732c343 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSturbulence_GridComp/shoc.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSturbulence_GridComp/shoc.F90 @@ -47,17 +47,17 @@ module shoc contains subroutine run_shoc( nx, ny, nzm, nz, dtn, & ! in - prsl_inv, phii_inv, phil_inv, & ! in + shf, prsl_inv, phii_inv, phil_inv, & ! in u_inv, v_inv, & ! in omega_inv, tabs_inv, qwv_inv, & ! in qi_inv, qc_inv, qpi_inv, & ! in qpl_inv, cld_sgs_inv, wthv_sec_inv, & ! in - wthv_mf_inv, tke_mf, zpbl, & ! in + wthv_mf_inv, tke_mf, dryzpbl, & ! in tke_inv, tkh_inv, & ! inout tkm_inv, isotropy_inv, & ! out tkesbdiss_inv, tkesbbuoy_inv, & ! out tkesbshear_inv, & ! out - smixt_inv, lmix_out, smixt1_inv, & ! out + smixt_inv, smixt1_inv, & ! out smixt2_inv,smixt3_inv, & ! out bruntmst_inv, ri_inv, prnum_inv, & ! out shocparams ) @@ -83,6 +83,7 @@ subroutine run_shoc( nx, ny, nzm, nz, dtn, & ! in integer, intent(in) :: nz ! Number of layer interfaces (= nzm + 1) type(shocparams_type), intent(in) :: shocparams real, intent(in ) :: dtn ! Physics time step, s + real, intent(in ) :: shf(nx,ny) ! Surface sensible heat flux, W/m2 real, intent(in ) :: prsl_inv (nx,ny,nzm) ! mean layer presure real, intent(in ) :: phii_inv (nx,ny,nz ) ! interface geopotential height real, intent(in ) :: phil_inv (nx,ny,nzm) ! layer geopotential height @@ -100,7 +101,7 @@ subroutine run_shoc( nx, ny, nzm, nz, dtn, & ! in real, intent(in ) :: qpi_inv (nx,ny,nzm) ! snow mixing ratio, kg/kg real, intent(in ) :: cld_sgs_inv(nx,ny,nzm) ! sgs cloud fraction real, intent(in ) :: tke_mf (nx,ny,nz) ! MF vertical velocity on edges, m/s - real, intent(in ) :: zpbl (nx,ny) ! PBLH diagnosed in TurbGridComp + real, intent(in ) :: dryzpbl (nx,ny) ! Dry mass flux depth real, intent(inout) :: tke_inv (nx,ny,nzm) ! turbulent kinetic energy. m**2/s**2 real, intent(inout) :: tkh_inv (nx,ny,nzm) ! eddy scalar diffusivity real, intent( out) :: tkm_inv (nx,ny,nzm) ! eddy momentum diffusivity @@ -111,7 +112,6 @@ subroutine run_shoc( nx, ny, nzm, nz, dtn, & ! in real, dimension(:,:,:), pointer :: tkesbshear_inv ! shear production real, dimension(:,:,:), pointer :: smixt_inv ! dissipation length scale - real, dimension(:,:), pointer :: lmix_out ! mixed layer depth real, dimension(:,:,:), pointer :: smixt1_inv ! length scale, term 1 real, dimension(:,:,:), pointer :: smixt2_inv ! length scale, term 2 real, dimension(:,:,:), pointer :: smixt3_inv ! length scale, term 3 @@ -124,11 +124,12 @@ subroutine run_shoc( nx, ny, nzm, nz, dtn, & ! in ! SHOC tunable parameters real :: lambda - real, parameter :: min_tke = 1e-4 ! Minumum TKE value, m**2/s**2 + real, parameter :: min_tke = 1e-3 ! Minumum TKE value, m**2/s**2 real, parameter :: max_tke = 10. ! Maximum TKE value, m**2/s**2 -! Maximum turbulent eddy length scale, m +! Limits on eddy length scale, m real, parameter :: max_eddy_length_scale = 2000. + real, parameter :: min_eddy_length_scale = 40. ! Maximum "return-to-isotropy" time scale, s real, parameter :: max_eddy_dissipation_time_scale = 2000. @@ -185,8 +186,6 @@ subroutine run_shoc( nx, ny, nzm, nz, dtn, & ! in real, dimension(nx,ny,nzm) :: total_water, brunt2, def2, thv, l_par real, dimension(nx,ny,nz) :: brunt_edge, brunt_dry, brunt_cld - real, dimension(nx,ny) :: l_inf, l_mix, zcb, lts!, l_par!, denom, numer, cldarr - real lstarn, depth, omn, betdz, betdze, bbb, & term, qsatt, dqsat, thedz, conv_var, & tkes, pval, pkap, thlsec, qwsec, & @@ -298,7 +297,6 @@ subroutine run_shoc( nx, ny, nzm, nz, dtn, & ! in if (associated(tkesbbuoy_inv)) tkesbbuoy_inv(:,:,1:nzm) = tkesbbuoy(:,:,nzm:1:-1) if (associated(tkesbshear_inv)) tkesbshear_inv(:,:,1:nzm) = tkesbshear(:,:,nzm:1:-1) - if (associated(lmix_out)) lmix_out(:,:) = l_mix if (associated(smixt_inv)) smixt_inv(:,:,1:nzm) = smixt(:,:,nzm:1:-1) if (associated(smixt1_inv)) smixt1_inv(:,:,1:nzm) = smixt1(:,:,nzm:1:-1) if (associated(smixt2_inv)) smixt2_inv(:,:,1:nzm) = smixt2(:,:,nzm:1:-1) @@ -322,7 +320,7 @@ subroutine tke_shoc() real grd,betdz,betdze,Cek,Cee,lstarn, bbb, omn, omp,qsatt,dqsat, smix, & buoy_sgs,a_prod_sh,a_prod_bu,a_diss,a_prod_bu_debug, buoy_sgs_debug, & - wrk, wrk1, wtke, wtk2, rdtn, tke_env + wrk, wrk1, wtke, wtk2, rdtn, tke_env, lambda_zfac integer i,j,k,ku,kd,itr real, dimension(nx,ny,nzm) :: tscale1 @@ -396,7 +394,7 @@ 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)) + wrk2 = min_tke*(1.+9.*exp(-zl(i,j,k)/100.))+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(wrk2, wtke), max_tke) a_diss = wrk*sqrt(wtke) ! Coefficient in the TKE dissipation term @@ -431,11 +429,13 @@ subroutine tke_shoc() do i=1,nx ! Calculate "return-to-isotropy" eddy dissipation time scale, see Eq. 8 in BK13 ! ignore stability dependence within the lower CBL, to prevent occasional - if (brunt_edge(i,j,k) <= 1e-5 .or. zl(i,j,k).lt.0.7*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)))) + !wrk = 2./(1./tscale1(i,j,k)+1./tscale1(i,j,k-1)) + wrk = 0.5*(tscale1(i,j,k)+tscale1(i,j,k-1)) + lambda_zfac = 0.5+0.5*tanh((zl(i,j,k)-0.75*dryzpbl(i,j)-100.)/100) + if (brunt_edge(i,j,k) <= 1e-5) then ! .or. zl(i,j,k).lt.0.7*dryzpbl(i,j)) then + isotropy(i,j,k) = max(30.,min(max_eddy_dissipation_time_scale,wrk)) 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*lambda_zfac*brunt_edge(i,j,k)*wrk*wrk))) endif if (tke(i,j,k).lt.2e-4) isotropy(i,j,k) = 30. @@ -650,9 +650,9 @@ subroutine eddy_length() + (bbb*fac_cond - (1.+fac_cond*dqsat)*tabs(i,j,k))*(qpl(i,j,kc)-qpl(i,j,kb)) & + (bbb*fac_sub - (1.+fac_sub*dqsat)*tabs(i,j,k))*(qpi(i,j,kc)-qpi(i,j,kb)) ) + if (k.gt.1) then bbb = 0.5*(bbb + (1. + epsv*qsatt-wrk-qpl(i,j,k-1)-qpi(i,j,k-1) & + 1.61*tabs(i,j,k-1)*dqsat) / (1.+lstarn*dqsat) ) - if (k.gt.1) then brunt_edge(i,j,k) = 0.5*(cld_sgs(i,j,k)+cld_sgs(i,j,k-1))*betdz*(bbb*(hl(i,j,k)-hl(i,j,k-1)) & + (bbb*lstarn - (1.+lstarn*dqsat)*tabs(i,j,k)) & * (total_water(i,j,k)-total_water(i,j,k-1)) & @@ -671,8 +671,8 @@ subroutine eddy_length() + (bbb*fac_sub -tabs(i,j,k))*(qpi(i,j,kc)-qpi(i,j,kb)) ) ! Calc outside-of-cloud on edges - bbb = 0.5*(bbb + 1. + epsv*qv(i,j,k-1) - qpl(i,j,k-1) - qpi(i,j,k-1)) if (k.gt.1) then + bbb = 0.5*(bbb + 1. + epsv*qv(i,j,k-1) - qpl(i,j,k-1) - qpi(i,j,k-1)) brunt_edge(i,j,k) = brunt_edge(i,j,k) + (1.-0.5*(cld_sgs(i,j,k)+cld_sgs(i,j,k-1)))*betdz*( bbb*(hl(i,j,k)-hl(i,j,k-1)) & + epsv*tabs(i,j,k)*(total_water(i,j,k)-total_water(i,j,k-1)) & + (bbb*fac_cond-tabs(i,j,k))*(qpl(i,j,k)-qpl(i,j,k-1)) & @@ -682,7 +682,7 @@ subroutine eddy_length() ! 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) < 1e-5 .or. zl(i,j,k).lt.0.7*zpbl(i,j)) then + if (brunt(i,j,k) < 1e-5 .or. zl(i,j,k).lt.0.75*dryzpbl(i,j)) then brunt2(i,j,k) = bruntmin else brunt2(i,j,k) = brunt(i,j,K) @@ -698,37 +698,48 @@ subroutine eddy_length() brunt2(:,:,nzm) = brunt2(:,:,nzm-1) -!=========== Length scale calculations =========== - do k=1,nzm + !=========== Length scale calculations =========== + ! Based on Eq. 10 in BK13 (Eq. 4.12 in Pete's dissertation) + do k=1,nzm-1 do j=1,ny do i=1,nx - tkes = sqrt(tke(i,j,k)) - -! Calculate turbulent length scale in the boundary layer. -! See Eq. 10 in BK13 (Eq. 4.12 in Pete's dissertation) - - !---------------------------------- - ! calculate parcel mixing length - !---------------------------------- - kk = k - wrk = thv(i,j,k)+0.2 ! 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 ! 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 ( shocparams%LENOPT .lt. 4 ) then ! SHOC-MF length scale + tkes = sqrt(tke(i,j,k)) + + !---------------------------------- + ! calculate parcel mixing length + !---------------------------------- + kk = k + wrk = thv(i,j,k)+0.2 ! upward T perturbation + do while (kk.lt.nzm .and. wrk .gt. thv(i,j,kk) ) + kk = kk+1 + end do + kk = kk-1 + + if (abs(thv(i,j,kk+1)-thv(i,j,kk)).gt.0.01) then + 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))) + else + l_par(i,j,k) = zl(i,j,kk) + end if + + kk = k + wrk = thv(i,j,k)-0.2 ! downward T perturbation + do while (kk .gt. 1 .and. wrk .lt. thv(i,j,kk)) + kk = kk-1 + end do + if ( kk.eq.1 .and. wrk.lt.thv(i,j,1) ) kk = kk-1 + kk = kk+1 + if (abs(thv(i,j,kk+1)-thv(i,j,kk)).gt.0.01) then + 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+1)-zl(i,j,kk))/(thv(i,j,kk+1)-thv(i,j,kk))) + else + l_par(i,j,k) = l_par(i,j,k) - zl(i,j,kk) + end if + l_par(i,j,k) = max(min(l_par(i,j,k),1500.),25.) + + + if ( shocparams%LENOPT .lt. 4 ) then ! Surface length scale smixt1(i,j,k) = vonk*zl(i,j,k)*shocparams%LENFAC1 @@ -738,8 +749,16 @@ 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.05,tkes)*shocparams%LENFAC3/(sqrt(brunt2(i,j,k))) - + ! Reduce sensitivity to stability within dry CBL or SBL, + ! but retain full sensitivity in free atmosphere. This is + ! a 'kludge' to increase TKE in stable BLs while suppressing + ! it in cumulus layers. + if ( zl(i,j,k).lt.0.75*dryzpbl(i,j) .or. zl(i,j,k).lt.500. ) then + smixt3(i,j,k) = max(0.05,tkes)*4.*shocparams%LENFAC3/(sqrt(brunt2(i,j,k))) + else + smixt3(i,j,k) = max(0.05,tkes)*shocparams%LENFAC3/(sqrt(brunt2(i,j,k))) + end if + !=== Combine component length scales === if (shocparams%LENOPT .eq. 1) then ! JPL blending approach (w/SHOC length scales) wrk1 = SQRT(3./(1./smixt2(i,j,k)**2+1./smixt3(i,j,k)**2)) @@ -749,36 +768,34 @@ subroutine eddy_length() smixt(i,j,k) = wrk1 end if else if (shocparams%LENOPT .eq. 2) then ! Harmonic mean - smixt(i,j,k) = min(max_eddy_length_scale, 3./(1./smixt1(i,j,k)+1./smixt2(i,j,k)+1./smixt3(i,j,k)) ) + smixt(i,j,k) = 3./(1./smixt1(i,j,k)+1./smixt2(i,j,k)+1./smixt3(i,j,k)) else if (shocparams%LENOPT .eq. 3) then ! SHOC classic approach - smixt(i,j,k) = min(max_eddy_length_scale, SQRT(3.)/SQRT(1./smixt1(i,j,k)**2+1./smixt2(i,j,k)**2+1./smixt3(i,j,k)**2) ) + smixt(i,j,k) = SQRT(3.)/SQRT(1./smixt1(i,j,k)**2+1./smixt2(i,j,k)**2+1./smixt3(i,j,k)**2) end if else if (shocparams%LENOPT .eq. 4) then ! JPL Length scale (Suselj et al 2012) 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*zpbl(i,j)))) + smixt(i,j,k) = 3.3*shocparams%LENFAC1*(wrk1 + (vonk*zl(i,j,k)-wrk1)*exp(-zl(i,j,k)/(0.1*dryzpbl(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) end if - ! Enforce minimum and maximum length scales - wrk = 40. !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 + ! Enforce minimum and maximum length scales, and reduce final + ! length scale exponentially with height above the dry CBL. + smixt(i,j,k) = min(max_eddy_length_scale,max(min_eddy_length_scale,smixt(i,j,k))) + if (zl(i,j,k).gt.dryzpbl(i,j)) smixt(i,j,k) = smixt(i,j,k)*(0.025+0.975*exp(-(zl(i,j,k)-dryzpbl(i,j))/3000.)) end do end do end do - - + smixt(:,:,nzm) = smixt(:,:,nzm-1) + smixt1(:,:,nzm) = smixt1(:,:,nzm-1) + smixt2(:,:,nzm) = smixt2(:,:,nzm-1) + smixt3(:,:,nzm) = smixt3(:,:,nzm-1) end subroutine eddy_length - end subroutine run_shoc @@ -786,6 +803,8 @@ subroutine update_moments( IM, JM, LM, & ! in DT, & ! in SH, & ! in EVAP, & ! in + AREA, & ! in + ZPBL, & ! in ZL, & ! in ZLE, & ! in KH, & ! in @@ -810,26 +829,22 @@ subroutine update_moments( IM, JM, LM, & ! in hl3, & ! out w2, & ! out w3, & ! out - w3can, & ! out -! wqt, & ! out - whl, & ! out hlqt, & ! out - qt2diag, & ! out - hl2diag, & ! out - hlqtdiag, & ! out doprogqt2, & ! tuning parameters hl2tune, & qt2tune, & hlqt2tune, & skew_tgen, & skew_tdis, & - docanuto ) + free_atm_qt2 ) integer, intent(in ) :: IM, JM, LM ! dimensions real, intent(in ) :: DT ! timestep [s] real, intent(in ) :: SH (IM,JM) ! surface sensible heat flux real, intent(in ) :: EVAP (IM,JM) ! surface evaporation + real, intent(in ) :: AREA (IM,JM) ! grid box area [m2] + real, intent(in ) :: ZPBL (IM,JM) ! grid box area [m2] real, intent(in ) :: ZL (IM,JM,LM) ! heights [m] real, intent(in ) :: ZLE (IM,JM,0:LM) ! edge heights [m] real, intent(in ) :: KH (IM,JM,0:LM) ! diffusivity @@ -854,56 +869,28 @@ subroutine update_moments( IM, JM, LM, & ! in real, intent( out) :: hl3 (IM,JM,LM) ! third moment static energy real, intent( out) :: w2 (IM,JM,LM) ! vertical velocity variance real, intent( out) :: w3 (IM,JM,LM) ! third moment vertical velocity - real, intent( out) :: w3can(IM,JM,LM) ! third moment vertical velocity -! real, intent( out) :: wqt (IM,JM,LM) ! vertical flux of total water - real, intent( out) :: whl (IM,JM,LM) ! vertical flux of liquid water static energy real, intent( out) :: hlqt (IM,JM,LM) ! total water, static energy covariance - real, intent( out) :: qt2diag(IM,JM,LM) - real, intent( out) :: hl2diag(IM,JM,LM) - real, intent( out) :: hlqtdiag(IM,JM,LM) real, intent(in ) :: HL2TUNE, & ! tuning parameters HLQT2TUNE, & QT2TUNE, & SKEW_TGEN, & - SKEW_TDIS + SKEW_TDIS, & + FREE_ATM_QT2 - integer, intent(in ) :: DOPROGQT2, & ! prognostic QT2 switch - DOCANUTO + integer, intent(in ) :: DOPROGQT2 ! prognostic QT2 switch real, parameter :: HL2MIN = 0.0005 - real, parameter :: HL2MAX = 2.0 + real, parameter :: HL2MAX = 0.25 ! Local variables integer :: k, kd, ku real, dimension(IM,JM) :: wrk1, wrk2, wrk3 real, dimension(IM,JM) :: sm, onemmf real, dimension(IM,JM,0:LM) :: qt2prod_edge, & - qt2prod_edge_nomf, & hl2_edge, & - hl2_edge_nomf, & - wqt_edge, & - whl_edge, & hlqt_edge,& qtgrad - real, dimension(IM,JM,LM) :: adzl, bet, whl_can -!======= Canuto variables - integer i, j, kb, kc, km1 - real bet2, f0, f1, f2, f3, f4, f5, iso, isosqr, & - omega0, omega1, omega2, X0, Y0, X1, Y1, AA0, AA1, buoy_sgs2, & - thedz, thedz2, cond, wrk, avew, wrk1b, wrk2b, wrk3b, dum -! See Eq. 7 in C01 (B.7 in Pete's dissertation) - real, parameter :: c=7.0, a0=0.52/(c*c*(c-2.)), a1=0.87/(c*c), & - a2=0.5/c, a3=0.6/(c*(c-2.)), a4=2.4/(3.*c+5.), & - a5=0.6/(c*(3.*c+5)) -!======== - - bet = 9.806/300. - - qt2diag = 0. - hl2diag = 0. - hlqtdiag = 0. - whl_edge(:,:,LM) = SH(:,:)/cp ! used only for Canuto below ! Initial calculations on edges @@ -915,19 +902,19 @@ subroutine update_moments( IM, JM, LM, & ! in ! SGS vertical flux liquid/ice water static energy. Eq 1 in BK13 wrk1 = HL(:,:,k) - HL(:,:,k+1) - whl_edge(:,:,k) = - wrk3 * wrk1 ! Second moment of liquid/ice water static energy. Eq 4 in BK13 - hl2_edge_nomf(:,:,k) = HL2TUNE * sm * wrk1 * wrk1 - hl2_edge(:,:,k) = HL2TUNE * 0.5*ISOTROPY(:,:,k) * & - (wrk3*wrk1-MFWHL(:,:,k)) * wrk1/(ZL(:,:,k)-ZL(:,:,k+1)) + hl2_edge(:,:,k) = HL2TUNE * 0.5*ISOTROPY(:,:,k) * (wrk3*wrk1-MFWHL(:,:,k)) & + * min(0.01,wrk1/(ZL(:,:,k)-ZL(:,:,k+1))) ! limit gradient to 1K/100m ! Total water gradient - qtgrad(:,:,k) = (QT(:,:,k) - QT(:,:,k+1)) / (ZL(:,:,k)-ZL(:,:,k+1)) + ! Here we limit very large negative gradients to 20%/100m + ! to avoid excessive QT2 production in stratocumulus layers. + wrk2 = (QT(:,:,k) - QT(:,:,k+1)) + qtgrad(:,:,k) = max(-0.002*qt(:,:,k+1),wrk2 / (ZL(:,:,k)-ZL(:,:,k+1))) ! Mean gradient production of total water variance, with and without MF contribution qt2prod_edge(:,:,k) = (KH(:,:,k)*qtgrad(:,:,k)-MFWQT(:,:,k)-0.*WQT_DC(:,:,k))*qtgrad(:,:,k) - qt2prod_edge_nomf(:,:,k) = (KH(:,:,k)*qtgrad(:,:,k))*qtgrad(:,:,k) ! Covariance of total water mixing ratio and liquid/ice water static energy. Eq 5 in BK13 hlqt_edge(:,:,k) = HLQT2TUNE * sm * wrk1 * wrk2 @@ -935,12 +922,13 @@ subroutine update_moments( IM, JM, LM, & ! in ! set lower boundary conditions hl2_edge(:,:,LM) = hl2_edge(:,:,LM-1) - hl2_edge_nomf(:,:,LM) = hl2_edge_nomf(:,:,LM-1) qt2prod_edge(:,:,LM) = qt2prod_edge(:,:,LM-1) - qt2prod_edge_nomf(:,:,LM) = qt2prod_edge_nomf(:,:,LM-1) hlqt_edge(:,:,LM) = hlqt_edge(:,:,LM-1) + ! value to relax QT2 above PBL, based on grid area + wrk3 = FREE_ATM_QT2*SQRT(SQRT(AREA/1.e10)) + ! Full level calculations do k=1,LM kd = k-1 @@ -949,41 +937,36 @@ subroutine update_moments( IM, JM, LM, & ! in onemmf = 1.0 - MFFRC(:,:,k) - w2(:,:,k) = onemmf*0.667*TKE(:,:,k) - + w2(:,:,k) = 1./(1./(0.667*TKE(:,:,k))+100./zl(:,:,k)) + hl2(:,:,k) = 0.5*( hl2_edge(:,:,kd) + hl2_edge(:,:,ku) ) - hl2diag(:,:,k) = 0.5*( hl2_edge_nomf(:,:,kd) + hl2_edge_nomf(:,:,ku) ) wrk1 = 0.5*(qt2prod_edge(:,:,kd)+qt2prod_edge(:,:,ku)) if (DOPROGQT2 /= 0) then -! wrk3 = QT2TUNE*1.5e-4 ! dissipation - qt2(:,:,k) = (qt2(:,:,k)+wrk1*DT) / (1. + DT/SKEW_TGEN) - qt2diag(:,:,k) = QT2TUNE*ISOTROPY(:,:,k)*0.5*(qt2prod_edge_nomf(:,:,kd)+qt2prod_edge_nomf(:,:,ku)) + ! Above the PBL, relax prognostic QT std dev to a specified fraction + ! of the mean total water. + where (zl(:,:,k).gt.zpbl+500. .and. mffrc(:,:,k).lt.1e-4) + wrk2 = (wrk3*qt(:,:,k))**2 + elsewhere + wrk2 = 0. + end where + qt2(:,:,k) = (qt2(:,:,k)+(wrk1+wrk2/SKEW_TDIS)*DT) / (1. + DT/SKEW_TDIS) else -! qt2(:,:,k) = QT2TUNE*ISOTROPY(:,:,k)*wrk1 qt2(:,:,k) = QT2TUNE*SKEW_TGEN*wrk1 - qt2diag(:,:,k) = 1.0*ISOTROPY(:,:,k)*0.5*(qt2prod_edge_nomf(:,:,kd)+qt2prod_edge_nomf(:,:,ku)) end if hlqt(:,:,k) = onemmf*0.5*( hlqt_edge(:,:,kd) + hlqt_edge(:,:,ku) ) + MFHLQT(:,:,k) - hlqtdiag(:,:,k) = 0.5*( hlqt_edge(:,:,kd) + hlqt_edge(:,:,ku) ) - - whl(:,:,k) = onemmf*0.5*( whl_edge(:,:,kd) + whl_edge(:,:,ku)) + 0.5*( mfwhl(:,:,kd) + mfwhl(:,:,ku) ) - whl_can(:,:,k) = whl(:,:,k) ! 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) - + ! Restrict HL variance hl2(:,:,k) = max(min(hl2(:,:,k),HL2MAX),HL2MIN) - hl2diag(:,:,k) = max(min(hl2diag(:,:,k),HL2MAX),HL2MIN) ! Ensure realizibility hlqt(:,:,k) = sign( min( abs(hlqt(:,:,k)), sqrt(hl2(:,:,k)*qt2(:,:,k)) ), hlqt(:,:,k) ) - hlqtdiag(:,:,k) = sign( min( abs(hlqtdiag(:,:,k)), sqrt(hl2diag(:,:,k)*qt2diag(:,:,k)) ), hlqtdiag(:,:,k) ) end do - + ! Update PDF_A and third moments if (DOPROGQT2 /= 0) then if (SKEW_TDIS.gt.0.) then @@ -1001,139 +984,9 @@ subroutine update_moments( IM, JM, LM, & ! in end if pdf_a = min(0.5,max(0.,pdf_a)) - if (DOCANUTO==0) then - qt3 = ( qt3 + max(MFQT3,0.)*DT/SKEW_TGEN ) / ( 1. + DT/SKEW_TDIS ) - hl3 = MFHL3 + qt3 = min(( qt3 + max(MFQT3,0.)*DT/SKEW_TGEN ) / ( 1. + DT/SKEW_TDIS ),10.*qt2**1.5) + hl3 = max(MFHL3,-10.*hl2**1.5) w3 = MFW3 - else - - -!============ Canuto 2001 estimate of third moments ============== -! This code was retained for diagnostic, comparative purpose only -!================================================================= - do k=1,LM - km1 = k - 1 - do j=1,JM - do i=1,IM - adzl(i,j,k) = (ZLE(i,j,km1) - ZLE(i,j,k)) ! level thickness - enddo - end do - end do - - do k=2,LM - - kb = k-1 - kc = k+1 - - do j=1,JM - do i=1,IM - - if(k == 1) then - kb = 1 - kc = 2 - thedz = adzl(i,j,kc) - thedz2 = thedz - elseif(k == LM) then - kb = LM-1 - kc = LM - thedz = adzl(i,j,k) - thedz2 = thedz - else - thedz = adzl(i,j,k) - thedz2 = adzl(i,j,k)+adzl(i,j,kb) - endif - - thedz = 1. / thedz - thedz2 = 1. / thedz2 - - iso = isotropy(i,j,k) !0.5*(isotropy(i,j,k)+isotropy(i,j,kb)) - isosqr = iso*iso ! Two-level average of "return-to-isotropy" time scale squared - buoy_sgs2 = isosqr*0.5*(brunt(i,j,k)+brunt(i,j,kb)) - bet2 = ggr/hl(i,j,k) !0.5*(bet(i,j,k)+bet(i,j,kb)) !Two-level average of BV frequency squared - -! Compute functions f0-f5, see Eq, 8 in C01 (B.8 in Pete's dissertation) - avew = w2(i,j,k) !0.5*(0.667*TKE(i,j,k)+0.667*TKE(i,j,kb)) -! if (abs(avew).ge.1e10) avew = sign(1e10,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) - wrk3b = hl2diag(i,j,kb) - hl2diag(i,j,kc) - - f0 = wrk2b * wrk1b * whl_can(i,j,k) * wrk3b - -! wrk = whl_can(i,j,kc) - whl_can(i,j,kb) - wrk = whl_can(i,j,kb) - whl_can(i,j,kc) - - f1 = wrk2b * (wrk*whl_can(i,j,k) + 0.5*w2(i,j,k)*wrk3b) - - wrk1b = bet2*isosqr - f2 = thedz*wrk1b*whl_can(i,j,k)*0.667*(TKE(i,j,kb)-TKE(i,j,k)) & - + (thedz2+thedz2)*bet2*isosqr*w2(i,j,k)*wrk - - f3 = thedz2*wrk1b*wrk*w2(i,j,k) + thedz*bet2*isosqr*(whl_can(i,j,k)*(tke(i,j,kb)-tke(i,j,k))) - - wrk1b = thedz*iso*w2(i,j,k) -! f4 = wrk1b*(0.667*TKE(i,j,kb)-0.667*TKE(i,j,k) + tke(i,j,kb)-tke(i,j,k)) - f4 = wrk1b*(w2(i,j,kb)-w2(i,j,k) + tke(i,j,kb)-tke(i,j,k)) - - f5 = wrk1b*0.667*(TKE(i,j,kb)-TKE(i,j,k)) - -! 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-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-20) dum = sign(1e-20,dum) - - wrk1b = 1.0 / dum - dum = 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) - X1 = wrk1b * (a0*f0+a1*f1+a2*(1.-a3*buoy_sgs2)*f2) - Y1 = wrk2b * (2.*a2*(buoy_sgs2*X1+(a0/a1)*f0+f1)) - -! Compute the A0, A1 terms, see Eq. 5d in C01 (B.5 in Pete's dissertation) - AA0 = omega0*X0 + omega1*Y0 - AA1 = omega0*X1 + omega1*Y1 + omega2 - -! Finally, we have the third moment of w, see Eq. 4c in C01 (B.4 in Pete's dissertation) -! cond is an estimate of third moment from second oment - If the third moment is larger -! than the estimate - limit w3. - - dum = c-1.2*X0+AA0 - 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)) - w3can(i,j,k) = (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. - - end do - end do - end do - do j=1,JM - do i=1,IM - 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 - - end if ! DOCANUTO conditional end subroutine update_moments