Skip to content

Commit

Permalink
Addressing comments from Dustin
Browse files Browse the repository at this point in the history
  • Loading branch information
@joeolson42
joeolson42 committed Oct 29, 2023
1 parent 872f488 commit 277bd48
Showing 1 changed file with 23 additions and 30 deletions.
53 changes: 23 additions & 30 deletions physics/module_bl_mynn.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3046,7 +3046,6 @@ SUBROUTINE mym_turbulence ( &
! q-variance (pdq), and covariance (pdc)
pdk(k) = elq*( sm(k)*gm(k) &
& +sh(k)*gh(k)+gamv ) + &
!fix? & TKEprodTD(k)
& 0.5*TKEprodTD(k) ! xmchen
pdt(k) = elh*( sh(k)*dtl(k)+gamt )*dtl(k)
pdq(k) = elh*( sh(k)*dqw(k)+gamq )*dqw(k)
Expand Down Expand Up @@ -3091,7 +3090,6 @@ SUBROUTINE mym_turbulence ( &
!qBUOY1D(k) = elq*(sh(k)*(-dTdz*grav/thl(k)) + gamv) !! ORIGINAL CODE

!! Buoyncy term takes the TKEprodTD(k) production now
!fix? qBUOY1D(k) = elq*(sh(k)*gh(k)+gamv)+TKEprodTD(k) !staggered
qBUOY1D(k) = elq*(sh(k)*gh(k)+gamv)+0.5*TKEprodTD(k) ! xmchen

!!!Dissipation Term (now it evaluated in mym_predict)
Expand Down Expand Up @@ -3801,7 +3799,7 @@ SUBROUTINE mym_condensation (kts,kte, &

!Diagnostic statistical scheme of Chaboureau and Bechtold (2002), JAS
!but with use of higher-order moments to estimate sigma
pblh2=MAX(10.,pblh1)
pblh2=MAX(10._kind_phys,pblh1)
zagl = 0.
dzm1 = 0.
DO k = kts,kte-1
Expand All @@ -3811,7 +3809,7 @@ SUBROUTINE mym_condensation (kts,kte, &
t = th(k)*exner(k)
xl = xl_blend(t) ! obtain latent heat
qsat_tk= qsat_blend(t, p(k)) ! saturation water vapor mixing ratio at tk and p
rh(k) = MAX(MIN(1.00,qw(k)/MAX(1.E-10,qsat_tk)),0.001)
rh(k) = MAX(MIN(1.0_kind_phys,qw(k)/MAX(1.E-10,qsat_tk)),0.001_kind_phys)

!dqw/dT: Clausius-Clapeyron
dqsl = qsat_tk*ep_2*xlv/( r_d*t**2 )
Expand Down Expand Up @@ -3840,7 +3838,7 @@ SUBROUTINE mym_condensation (kts,kte, &
!sgm(k) = max( sgm(k), qsat_tk*0.035 )

!introduce vertical grid spacing dependence on min sgm
wt = max(500. - max(dz(k)-100.,0.0), 0.0)/500. !=1 for dz < 100 m, =0 for dz > 600 m
wt = max(500. - max(dz(k)-100.,0.0), 0.0_kind_phys)/500. !=1 for dz < 100 m, =0 for dz > 600 m
sgm(k) = sgm(k) + sgm(k)*0.2*(1.0-wt) !inflate sgm for coarse dz

!allow min sgm to vary with dz and z.
Expand All @@ -3855,23 +3853,23 @@ SUBROUTINE mym_condensation (kts,kte, &
rh_hack= rh(k)
!ensure adequate RH & q1 when qi is at least 1e-9
if (qi(k)>1.e-9) then
rh_hack =min(1.0, rhcrit + 0.07*(9.0 + log10(qi(k))))
rh_hack =min(1.0_kind_phys, rhcrit + 0.07*(9.0 + log10(qi(k))))
rh(k) =max(rh(k), rh_hack)
!add rh-based q1
q1_rh =-3. + 3.*(rh_hack-rhcrit)/(1.-rhcrit)
q1(k) =max(q1_rh, q1(k) )
endif
!ensure adequate RH & q1 when qc is at least 1e-6
if (qc(k)>1.e-6) then
rh_hack =min(1.0, rhcrit + 0.09*(6.0 + log10(qc(k))))
rh_hack =min(1.0_kind_phys, rhcrit + 0.09*(6.0 + log10(qc(k))))
rh(k) =max(rh(k), rh_hack)
!add rh-based q1
q1_rh =-3. + 3.*(rh_hack-rhcrit)/(1.-rhcrit)
q1(k) =max(q1_rh, q1(k) )
endif
!ensure adequate RH & q1 when qs is at least 1e-7 (above the PBLH)
if (qs(k)>1.e-8 .and. zagl .gt. pblh2) then
rh_hack =min(1.0, rhcrit + 0.07*(8.0 + log10(qs(k))))
rh_hack =min(1.0_kind_phys, rhcrit + 0.07*(8.0 + log10(qs(k))))
rh(k) =max(rh(k), rh_hack)
!add rh-based q1
q1_rh =-3. + 3.*(rh_hack-rhcrit)/(1.-rhcrit)
Expand Down Expand Up @@ -3953,10 +3951,10 @@ SUBROUTINE mym_condensation (kts,kte, &
elseif (q1k .ge. -2.5 .and. q1k .lt. -1.7) then
Fng = 3.0 + exp(-3.8*(q1k+1.7))
else
Fng = min(23.9 + exp(-1.6*(q1k+2.5)), 60.)
Fng = min(23.9 + exp(-1.6*(q1k+2.5)), 60._kind_phys)
endif

cfmax = min(cldfra_bl1D(k), 0.6)
cfmax = min(cldfra_bl1D(k), 0.6_kind_phys)
!Further limit the cf going into vt & vq near the surface
zsl = min(max(25., 0.1*pblh2), 100.)
wt = min(zagl/zsl, 1.0) !=0 at z=0 m, =1 above ekman layer
Expand Down Expand Up @@ -4921,9 +4919,6 @@ SUBROUTINE mynn_tendencies(kts,kte,i, &
sqi2(k) = 0.0 ! if sqw2 > qsat
sqc2(k) = 0.0
ENDIF
!dqv(k) = (sqv2(k) - sqv(k))/delt
!dqc(k) = (sqc2(k) - sqc(k))/delt
!dqi(k) = (sqi2(k) - sqi(k))/delt
ENDDO
ENDIF

Expand All @@ -4932,7 +4927,7 @@ SUBROUTINE mynn_tendencies(kts,kte,i, &
! WATER VAPOR TENDENCY
!=====================
DO k=kts,kte
Dqv(k)=(sqv2(k)/(1.-sqv2(k)) - qv(k))/delt
Dqv(k)=(sqv2(k) - sqv(k))/delt
!if (sqv2(k) < 0.0)print*,"neg qv:",sqv2(k),k
ENDDO

Expand All @@ -4943,7 +4938,7 @@ SUBROUTINE mynn_tendencies(kts,kte,i, &
!print*,"FLAG_QC:",FLAG_QC
IF (FLAG_QC) THEN
DO k=kts,kte
Dqc(k)=(sqc2(k)/(1.-sqv2(k)) - qc(k))/delt
Dqc(k)=(sqc2(k) - sqc(k))/delt
!if (sqc2(k) < 0.0)print*,"neg qc:",sqc2(k),k
ENDDO
ELSE
Expand Down Expand Up @@ -4971,7 +4966,7 @@ SUBROUTINE mynn_tendencies(kts,kte,i, &
!===================
IF (FLAG_QI) THEN
DO k=kts,kte
Dqi(k)=(sqi2(k)/(1.-sqv2(k)) - qi(k))/delt
Dqi(k)=(sqi2(k) - sqi(k))/delt
!if (sqi2(k) < 0.0)print*,"neg qi:",sqi2(k),k
ENDDO
ELSE
Expand All @@ -4985,7 +4980,7 @@ SUBROUTINE mynn_tendencies(kts,kte,i, &
!===================
IF (.false.) THEN !disabled
DO k=kts,kte
Dqs(k)=(sqs2(k)/(1.-sqs2(k)) - qs(k))/delt
Dqs(k)=(sqs2(k) - sqs(k))/delt
ENDDO
ELSE
DO k=kts,kte
Expand All @@ -5009,10 +5004,11 @@ SUBROUTINE mynn_tendencies(kts,kte,i, &
ELSE !-MIX CLOUD SPECIES?
!CLOUDS ARE NOT NIXED (when bl_mynn_cloudmix == 0)
DO k=kts,kte
Dqc(k)=0.
Dqc(k) =0.
Dqnc(k)=0.
Dqi(k)=0.
Dqi(k) =0.
Dqni(k)=0.
Dqs(k) =0.
ENDDO
ENDIF

Expand Down Expand Up @@ -6007,30 +6003,27 @@ SUBROUTINE DMP_mf( &
! Criteria (1)
maxwidth = min(dx*dcut, lmax)
!Criteria (2)
maxwidth = min(maxwidth, 1.1*PBLH)
maxwidth = min(maxwidth, 1.1_kind_phys*PBLH)
! Criteria (3)
if ((landsea-1.5) .lt. 0) then !land
maxwidth = MIN(maxwidth, 0.5*cloud_base)
maxwidth = MIN(maxwidth, 0.5_kind_phys*cloud_base)
else !water
maxwidth = MIN(maxwidth, 0.9*cloud_base)
maxwidth = MIN(maxwidth, 0.9_kind_phys*cloud_base)
endif
! Criteria (4)
wspd_pbl=SQRT(MAX(u(kts)**2 + v(kts)**2, 0.01))
wspd_pbl=SQRT(MAX(u(kts)**2 + v(kts)**2, 0.01_kind_phys))
!Note: area fraction (acfac) is modified below
! Criteria (5) - only a function of flt (not fltv)
if ((landsea-1.5).LT.0) then !land
!width_flx = MAX(MIN(1000.*(0.6*tanh((flt - 0.050)/0.03) + .5),1000.), 0.)
!width_flx = MAX(MIN(1000.*(0.6*tanh((flt - 0.040)/0.03) + .5),1000.), 0.)
width_flx = MAX(MIN(1000.*(0.6*tanh((fltv - 0.040)/0.04) + .5),1000.), 0.)
width_flx = MAX(MIN(1000.*(0.6*tanh((fltv - 0.040)/0.04) + .5),1000._kind_phys), 0._kind_phys)
else !water
!width_flx = MAX(MIN(1000.*(0.6*tanh((flt - 0.003)/0.01) + .5),1000.), 0.)
width_flx = MAX(MIN(1000.*(0.6*tanh((fltv - 0.007)/0.02) + .5),1000.), 0.)
width_flx = MAX(MIN(1000.*(0.6*tanh((fltv - 0.007)/0.02) + .5),1000._kind_phys), 0._kind_phys)
endif
maxwidth = MIN(maxwidth, width_flx)
minwidth = lmin
!allow min plume size to increase in large flux conditions (eddy diffusivity should be
!large enough to handle the representation of small plumes).
if (maxwidth .ge. (lmax - 1.0) .and. fltv .gt. 0.2)minwidth = lmin + dlmin*min((fltv-0.2)/0.3, 1.)
if (maxwidth .ge. (lmax - 1.0) .and. fltv .gt. 0.2)minwidth = lmin + dlmin*min((fltv-0.2)/0.3, 1._kind_phys)

if (maxwidth .le. minwidth) then ! deactivate MF component
nup2 = 0
Expand All @@ -6048,7 +6041,7 @@ SUBROUTINE DMP_mf( &
! Find coef C for number size density N
cn = 0.
d =-1.9 !set d to value suggested by Neggers 2015 (JAMES).
dl = (maxwidth - minwidth)/real(nup-1)
dl = (maxwidth - minwidth)/real(nup-1,kind=kind_phys)
do i=1,NUP
! diameter of plume
l = minwidth + dl*real(i-1)
Expand Down

0 comments on commit 277bd48

Please sign in to comment.