Skip to content

Commit

Permalink
Merge pull request NCAR#982 from grantfirl/ufs-dev-PR19
Browse files Browse the repository at this point in the history 
UFS-dev PR#19
  • Loading branch information
@grantfirl
grantfirl authored Nov 21, 2022
2 parents 6374c88 + c65aa90 commit 312f63e
Show file tree
Hide file tree
Showing 5 changed files with 100 additions and 34 deletions.
13 changes: 10 additions & 3 deletions physics/GFS_rrtmg_pre.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -73,7 +73,8 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, lextop, ltp, &
use surface_perturbation, only: cdfnor,ppfbet

! For Thompson MP
use module_mp_thompson, only: calc_effectRad, Nt_c, &
use module_mp_thompson, only: calc_effectRad, &
Nt_c_l, Nt_c_o, &
re_qc_min, re_qc_max, &
re_qi_min, re_qi_max, &
re_qs_min, re_qs_max
Expand Down Expand Up @@ -245,6 +246,7 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, lextop, ltp, &
real (kind=kind_phys) :: alpha0,beta0,m,s,cldtmp,tmp_wt,cdfz
real (kind=kind_phys) :: max_relh
integer :: iflag
integer :: islmsk

integer :: ids, ide, jds, jde, kds, kde, &
ims, ime, jms, jme, kms, kme, &
Expand Down Expand Up @@ -748,7 +750,11 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, lextop, ltp, &
qc_mp (i,k) = tracer1(i,k,ntcw)/(1.-qvs)
qi_mp (i,k) = tracer1(i,k,ntiw)/(1.-qvs)
qs_mp (i,k) = tracer1(i,k,ntsw)/(1.-qvs)
nc_mp (i,k) = nt_c*orho(i,k)
if(nint(slmsk(i)) == 1) then
nc_mp (i,k) = Nt_c_l*orho(i,k)
else
nc_mp (i,k) = Nt_c_o*orho(i,k)
endif
ni_mp (i,k) = tracer1(i,k,ntinc)/(1.-qvs)
enddo
enddo
Expand Down Expand Up @@ -877,13 +883,14 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, lextop, ltp, &
end do
!> - Call Thompson's subroutine calc_effectRad() to compute effective radii
do i=1,im
islmsk = nint(slmsk(i))
! Effective radii [m] are now intent(out), bounds applied in calc_effectRad
!tgs: progclduni has different limits for ice radii (10.0-150.0) than
! calc_effectRad (4.99-125.0 for WRFv3.8.1; 2.49-125.0 for WRFv4+)
! it will raise the low limit from 5 to 10, but the high limit will remain 125.
call calc_effectRad (tlyr(i,:), plyr(i,:)*100., qv_mp(i,:), qc_mp(i,:), &
nc_mp(i,:), qi_mp(i,:), ni_mp(i,:), qs_mp(i,:), &
effrl(i,:), effri(i,:), effrs(i,:), 1, lm )
effrl(i,:), effri(i,:), effrs(i,:), islmsk, 1, lm )
! Scale Thompson's effective radii from meter to micron
do k=1,lm
effrl(i,k) = MAX(re_qc_min, MIN(effrl(i,k), re_qc_max))*1.e6
Expand Down
20 changes: 14 additions & 6 deletions physics/GFS_rrtmgp_cloud_mp.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,8 +12,8 @@ module GFS_rrtmgp_cloud_mp
use rrtmgp_lw_cloud_optics, only: &
radliq_lwr => radliq_lwrLW, radliq_upr => radliq_uprLW,&
radice_lwr => radice_lwrLW, radice_upr => radice_uprLW
use module_mp_thompson, only: calc_effectRad, Nt_c, re_qc_min, re_qc_max, re_qi_min, &
re_qi_max, re_qs_min, re_qs_max
use module_mp_thompson, only: calc_effectRad, Nt_c_l, Nt_c_o, re_qc_min, re_qc_max, &
re_qi_min, re_qi_max, re_qs_min, re_qs_max
use module_mp_thompson_make_number_concentrations, only: make_IceNumber, &
make_DropletNumber, make_RainNumber

Expand Down Expand Up @@ -254,7 +254,7 @@ subroutine GFS_rrtmgp_cloud_mp_run(nCol, nLev, nTracers, ncnd, i_cldliq, i_cldic
! Update particle size using modified mixing-ratios from Thompson.
call cmp_reff_Thompson(nLev, nCol, i_cldliq, i_cldice, i_cldsnow, i_cldice_nc, &
i_cldliq_nc, i_twa, q_lay, p_lay, t_lay, tracer, con_eps, con_rd, ltaerosol,&
mraerosol, effrin_cldliq, effrin_cldice, effrin_cldsnow)
mraerosol, lsmask, effrin_cldliq, effrin_cldice, effrin_cldsnow)
cld_reliq = effrin_cldliq
cld_reice = effrin_cldice
cld_resnow = effrin_cldsnow
Expand Down Expand Up @@ -820,7 +820,7 @@ function cld_frac_XuRandall(p_lay, qs_lay, relhum, cld_mr, alpha)
!! \section cmp_reff_Thompson_gen General Algorithm
subroutine cmp_reff_Thompson(nLev, nCol, i_cldliq, i_cldice, i_cldsnow, i_cldice_nc, &
i_cldliq_nc, i_twa, q_lay, p_lay, t_lay, tracer, con_eps, con_rd, ltaerosol, &
mraerosol, effrin_cldliq, effrin_cldice, effrin_cldsnow)
mraerosol, lsmask, effrin_cldliq, effrin_cldice, effrin_cldsnow)
implicit none

! Inputs
Expand All @@ -830,6 +830,7 @@ subroutine cmp_reff_Thompson(nLev, nCol, i_cldliq, i_cldice, i_cldsnow, i_cldice
real(kind_phys), intent(in) :: con_eps,con_rd
real(kind_phys), dimension(:,:),intent(in) :: q_lay, p_lay, t_lay
real(kind_phys), dimension(:,:,:),intent(in) :: tracer
real(kind_phys), dimension(:), intent(in) :: lsmask

! Outputs
real(kind_phys), dimension(:,:), intent(inout) :: effrin_cldliq, effrin_cldice, &
Expand All @@ -840,6 +841,7 @@ subroutine cmp_reff_Thompson(nLev, nCol, i_cldliq, i_cldice, i_cldsnow, i_cldice
real(kind_phys) :: rho, orho
real(kind_phys),dimension(nCol,nLev) :: qv_mp, qc_mp, qi_mp, qs_mp, ni_mp, nc_mp, &
nwfa, re_cloud, re_ice, re_snow
integer :: ilsmask

! Prepare cloud mixing-ratios and number concentrations for calc_effectRa
do iLay = 1, nLev
Expand All @@ -863,7 +865,11 @@ subroutine cmp_reff_Thompson(nLev, nCol, i_cldliq, i_cldice, i_cldsnow, i_cldice
nc_mp(iCol,iLay) = make_DropletNumber(qc_mp(iCol,iLay)*rho, nwfa(iCol,iLay)*rho) * orho
endif
else
nc_mp(iCol,iLay) = nt_c*orho
if (nint(lsmask(iCol)) == 1) then !land
nc_mp(iCol,iLay) = nt_c_l*orho
else
nc_mp(iCol,iLay) = nt_c_o*orho
endif
endif
if (qi_mp(iCol,iLay) > 1.e-12 .and. ni_mp(iCol,iLay) < 100.) then
ni_mp(iCol,iLay) = make_IceNumber(qi_mp(iCol,iLay)*rho, t_lay(iCol,iLay)) * orho
Expand All @@ -873,9 +879,11 @@ subroutine cmp_reff_Thompson(nLev, nCol, i_cldliq, i_cldice, i_cldsnow, i_cldice

! Compute effective radii for liquid/ice/snow.
do iCol=1,nCol
ilsmask = nint(lsmask(iCol))
call calc_effectRad (t_lay(iCol,:), p_lay(iCol,:), qv_mp(iCol,:), qc_mp(iCol,:), &
nc_mp(iCol,:), qi_mp(iCol,:), ni_mp(iCol,:), qs_mp(iCol,:), &
re_cloud(iCol,:), re_ice(iCol,:), re_snow(iCol,:), 1, nLev )
re_cloud(iCol,:), re_ice(iCol,:), re_snow(iCol,:), ilsmask, &
1, nLev )
do iLay = 1, nLev
re_cloud(iCol,iLay) = MAX(re_qc_min, MIN(re_cloud(iCol,iLay), re_qc_max))
re_ice(iCol,iLay) = MAX(re_qi_min, MIN(re_ice(iCol,iLay), re_qi_max))
Expand Down
82 changes: 62 additions & 20 deletions physics/module_mp_thompson.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -92,7 +92,9 @@ MODULE module_mp_thompson
!.. scheme. In 2-moment cloud water, Nt_c represents a maximum of
!.. droplet concentration and nu_c is also variable depending on local
!.. droplet number concentration.
REAL, PARAMETER :: Nt_c = 100.E6
!REAL, PARAMETER :: Nt_c = 100.E6
REAL, PARAMETER :: Nt_c_o = 50.E6
REAL, PARAMETER :: Nt_c_l = 100.E6
REAL, PARAMETER, PRIVATE:: Nt_c_max = 1999.E6

!..Declaration of constants for assumed CCN/IN aerosols when none in
Expand All @@ -108,7 +110,7 @@ MODULE module_mp_thompson
REAL, PARAMETER, PRIVATE:: mu_r = 0.0
REAL, PARAMETER, PRIVATE:: mu_g = 0.0
REAL, PARAMETER, PRIVATE:: mu_i = 0.0
REAL, PRIVATE:: mu_c
REAL, PRIVATE:: mu_c_o, mu_c_l

!..Sum of two gamma distrib for snow (Field et al. 2005).
!.. N(D) = M2**4/M3**3 * [Kap0*exp(-M2*Lam0*D/M3)
Expand Down Expand Up @@ -150,7 +152,6 @@ MODULE module_mp_thompson
REAL, PARAMETER, PRIVATE:: fv_s = 100.0
REAL, PARAMETER, PRIVATE:: av_g = 442.0
REAL, PARAMETER, PRIVATE:: bv_g = 0.89
REAL, PARAMETER, PRIVATE:: av_i = 1493.9
REAL, PARAMETER, PRIVATE:: bv_i = 1.0
REAL, PARAMETER, PRIVATE:: av_c = 0.316946E8
REAL, PARAMETER, PRIVATE:: bv_c = 2.0
Expand Down Expand Up @@ -534,7 +535,8 @@ SUBROUTINE thompson_init(is_aerosol_aware_in, &
!.. disp=SQRT((mu+2)/(mu+1) - 1) so mu varies from 15 for Maritime
!.. to 2 for really dirty air. This not used in 2-moment cloud water
!.. scheme and nu_c used instead and varies from 2 to 15 (integer-only).
mu_c = MIN(15., (1000.E6/Nt_c + 2.))
mu_c_l = MIN(15., (1000.E6/Nt_c_l + 2.))
mu_c_o = MIN(15., (1000.E6/Nt_c_o + 2.))

!> - Compute Schmidt number to one-third used numerous times
Sc3 = Sc**(1./3.)
Expand Down Expand Up @@ -889,7 +891,7 @@ SUBROUTINE thompson_init(is_aerosol_aware_in, &

if (mpirank==mpiroot) write (*,*)'creating microphysics lookup tables ... '
if (mpirank==mpiroot) write (*,'(a, f5.2, a, f5.2, a, f5.2, a, f5.2)') &
' using: mu_c=',mu_c,' mu_i=',mu_i,' mu_r=',mu_r,' mu_g=',mu_g
' using: mu_c_o=',mu_c_o,' mu_i=',mu_i,' mu_r=',mu_r,' mu_g=',mu_g

!> - Call table_ccnact() to read a static file containing CCN activation of aerosols. The
!! data were created from a parcel model by Feingold & Heymsfield with
Expand Down Expand Up @@ -982,7 +984,7 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, &
nwfa, nifa, nwfa2d, nifa2d, &
tt, th, pii, &
p, w, dz, dt_in, dt_inner, &
sedi_semi, decfl, &
sedi_semi, decfl, lsm, &
RAINNC, RAINNCV, &
SNOWNC, SNOWNCV, &
ICENC, ICENCV, &
Expand Down Expand Up @@ -1037,6 +1039,7 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, &
REAL, DIMENSION(ims:ime, kms:kme, jms:jme), OPTIONAL, INTENT(INOUT):: &
nc, nwfa, nifa
REAL, DIMENSION(ims:ime, jms:jme), OPTIONAL, INTENT(IN):: nwfa2d, nifa2d
INTEGER, DIMENSION(ims:ime, jms:jme), INTENT(IN):: lsm
REAL, DIMENSION(ims:ime, kms:kme, jms:jme), OPTIONAL, INTENT(INOUT):: &
re_cloud, re_ice, re_snow
REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT):: pfils, pflls
Expand Down Expand Up @@ -1117,6 +1120,7 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, &
REAL, DIMENSION(its:ite, jts:jte):: pcp_ra, pcp_sn, pcp_gr, pcp_ic
REAL:: dt, pptrain, pptsnow, pptgraul, pptice
REAL:: qc_max, qr_max, qs_max, qi_max, qg_max, ni_max, nr_max
INTEGER:: lsml
REAL:: rand1, rand2, rand3, rand_pert_max
INTEGER:: i, j, k, m
INTEGER:: imax_qc,imax_qr,imax_qi,imax_qs,imax_qg,imax_ni,imax_nr
Expand Down Expand Up @@ -1419,8 +1423,13 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, &
nifa1d(k) = nifa(i,k,j)
enddo
else
lsml = lsm(i,j)
do k = kts, kte
nc1d(k) = Nt_c/rho(k)
if(lsml == 1) then
nc1d(k) = Nt_c_l/rho(k)
else
nc1d(k) = Nt_c_o/rho(k)
endif
nwfa1d(k) = 11.1E6
nifa1d(k) = naIN1*0.01
enddo
Expand All @@ -1429,7 +1438,7 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, &
!> - Call mp_thompson()
call mp_thompson(qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, &
nr1d, nc1d, nwfa1d, nifa1d, t1d, p1d, w1d, dz1d, &
pptrain, pptsnow, pptgraul, pptice, &
lsml, pptrain, pptsnow, pptgraul, pptice, &
#if ( WRF_CHEM == 1 )
rainprod1d, evapprod1d, &
#endif
Expand Down Expand Up @@ -1698,7 +1707,7 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, &
enddo
!> - Call calc_effectrad()
call calc_effectRad (t1d, p1d, qv1d, qc1d, nc1d, qi1d, ni1d, qs1d, &
re_qc1d, re_qi1d, re_qs1d, kts, kte)
re_qc1d, re_qi1d, re_qs1d, lsml, kts, kte)
do k = kts, kte
re_cloud(i,k,j) = MAX(re_qc_min, MIN(re_qc1d(k), re_qc_max))
re_ice(i,k,j) = MAX(re_qi_min, MIN(re_qi1d(k), re_qi_max))
Expand Down Expand Up @@ -1841,7 +1850,7 @@ END SUBROUTINE thompson_finalize
!> @{
subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, &
nr1d, nc1d, nwfa1d, nifa1d, t1d, p1d, w1d, dzq, &
pptrain, pptsnow, pptgraul, pptice, &
lsml, pptrain, pptsnow, pptgraul, pptice, &
#if ( WRF_CHEM == 1 )
rainprod, evapprod, &
#endif
Expand Down Expand Up @@ -1879,6 +1888,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, &
REAL, DIMENSION(kts:kte), INTENT(IN):: p1d, w1d, dzq
REAL, INTENT(INOUT):: pptrain, pptsnow, pptgraul, pptice
REAL, INTENT(IN):: dt
INTEGER, INTENT(IN):: lsml
REAL, INTENT(IN):: rand1, rand2, rand3
! Extended diagnostics, most arrays only allocated if ext_diag is true
LOGICAL, INTENT(IN) :: ext_diag
Expand Down Expand Up @@ -1982,6 +1992,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, &
REAL:: Ef_ra, Ef_sa, Ef_ga
REAL:: dtsave, odts, odt, odzq, hgt_agl, SR
REAL:: xslw1, ygra1, zans1, eva_factor
REAL:: av_i
INTEGER:: i, k, k2, n, nn, nstep, k_0, kbot, IT, iexfrq
INTEGER, DIMENSION(5):: ksed1
INTEGER:: nir, nis, nig, nii, nic, niin
Expand All @@ -2006,6 +2017,8 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, &
odt = 1./dt
odts = 1./dtsave
iexfrq = 1
! Transition value of coefficient matching at crossover from cloud ice to snow
av_i = av_s * D0s ** (bv_s - bv_i)

!+---+-----------------------------------------------------------------+
!> - Initialize Source/sink terms. First 2 chars: "pr" represents source/sink of
Expand Down Expand Up @@ -2210,7 +2223,13 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, &
endif
nc(k) = MIN( DBLE(Nt_c_max), ccg(1,nu_c)*ocg2(nu_c)*rc(k) &
/ am_r*lamc**bm_r)
if (.NOT. (is_aerosol_aware .or. merra2_aerosol_aware)) nc(k) = Nt_c
if (.NOT. (is_aerosol_aware .or. merra2_aerosol_aware)) then
if (lsml == 1) then
nc(k) = Nt_c_l
else
nc(k) = Nt_c_o
endif
endif
else
qc1d(k) = 0.0
nc1d(k) = 0.0
Expand Down Expand Up @@ -2919,13 +2938,13 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, &

!> - Deposition nucleation of dust/mineral from DeMott et al (2010)
!! we may need to relax the temperature and ssati constraints.
if ( (ssati(k).ge. 0.25) .or. (ssatw(k).gt. eps &
if ( (ssati(k).ge. 0.15) .or. (ssatw(k).gt. eps &
.and. temp(k).lt.253.15) ) then
if (dustyIce .AND. (is_aerosol_aware .or. merra2_aerosol_aware)) then
xnc = iceDeMott(tempc,qv(k),qvs(k),qvsi(k),rho(k),nifa(k))
xnc = xnc*(1.0 + 50.*rand3)
else
xnc = MIN(250.E3, TNO*EXP(ATO*(T_0-temp(k))))
xnc = MIN(1000.E3, TNO*EXP(ATO*(T_0-temp(k))))
endif
xni = ni(k) + (pni_rfz(k)+pni_wfz(k))*dtsave
pni_inu(k) = 0.5*(xnc-xni + abs(xnc-xni))*odts
Expand Down Expand Up @@ -3273,7 +3292,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, &
lami = cie(2)/5.E-6
xni = MIN(4999.D3, cig(1)*oig2*xri/am_i*lami**bm_i)
niten(k) = (xni-ni1d(k)*rho(k))*odts*orho
elseif (xDi.gt. 300.E-6) then
elseif (xDi.gt. 300.E-6) then
lami = cie(2)/300.E-6
xni = cig(1)*oig2*xri/am_i*lami**bm_i
niten(k) = (xni-ni1d(k)*rho(k))*odts*orho
Expand Down Expand Up @@ -3389,7 +3408,13 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, &
if ((qc1d(k) + qcten(k)*DT) .gt. R1) then
rc(k) = (qc1d(k) + qcten(k)*DT)*rho(k)
nc(k) = MAX(2., MIN((nc1d(k)+ncten(k)*DT)*rho(k), Nt_c_max))
if (.NOT. (is_aerosol_aware .or. merra2_aerosol_aware)) nc(k) = Nt_c
if (.NOT. (is_aerosol_aware .or. merra2_aerosol_aware)) then
if(lsml == 1) then
nc(k) = Nt_c_l
else
nc(k) = Nt_c_o
endif
endif
L_qc(k) = .true.
else
rc(k) = R1
Expand Down Expand Up @@ -3560,7 +3585,11 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, &
if (is_aerosol_aware .or. merra2_aerosol_aware) then
xnc = MAX(2., activ_ncloud(temp(k), w1d(k)+rand3, nwfa(k)))
else
xnc = Nt_c
if(lsml == 1) then
xnc = Nt_c_l
else
xnc = Nt_c_o
endif
endif
pnc_wcd(k) = 0.5*(xnc-nc(k) + abs(xnc-nc(k)))*odts*orho

Expand Down Expand Up @@ -3630,7 +3659,13 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, &
rc(k) = MAX(R1, (qc1d(k) + DT*qcten(k))*rho(k))
if (rc(k).eq.R1) L_qc(k) = .false.
nc(k) = MAX(2., MIN((nc1d(k)+ncten(k)*DT)*rho(k), Nt_c_max))
if (.NOT. (is_aerosol_aware .or. merra2_aerosol_aware)) nc(k) = Nt_c
if (.NOT. (is_aerosol_aware .or. merra2_aerosol_aware)) then
if(lsml == 1) then
nc(k) = Nt_c_l
else
nc(k) = Nt_c_o
endif
endif
qv(k) = MAX(1.E-10, qv1d(k) + DT*qvten(k))
temp(k) = t1d(k) + DT*tten(k)
rho(k) = 0.622*pres(k)/(R*temp(k)*(qv(k)+0.622))
Expand Down Expand Up @@ -4235,7 +4270,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, &
xDi = (bm_i + mu_i + 1.) * ilami
if (xDi.lt. 5.E-6) then
lami = cie(2)/5.E-6
elseif (xDi.gt. 300.E-6) then
elseif (xDi.gt. 300.E-6) then
lami = cie(2)/300.E-6
endif
ni1d(k) = MIN(cig(1)*oig2*qi1d(k)/am_i*lami**bm_i, &
Expand Down Expand Up @@ -5749,7 +5784,7 @@ END FUNCTION delta_p
!! distribution, not the second part, which is the larger sizes.

subroutine calc_effectRad (t1d, p1d, qv1d, qc1d, nc1d, qi1d, ni1d, qs1d, &
& re_qc1d, re_qi1d, re_qs1d, kts, kte)
& re_qc1d, re_qi1d, re_qs1d, lsml, kts, kte)

IMPLICIT NONE

Expand All @@ -5766,6 +5801,7 @@ subroutine calc_effectRad (t1d, p1d, qv1d, qc1d, nc1d, qi1d, ni1d, qs1d, &
DOUBLE PRECISION:: lamc, lami
LOGICAL:: has_qc, has_qi, has_qs
INTEGER:: inu_c
INTEGER:: lsml
real, dimension(15), parameter:: g_ratio = (/24,60,120,210,336, &
& 504,720,990,1320,1716,2184,2730,3360,4080,4896/)

Expand All @@ -5781,7 +5817,13 @@ subroutine calc_effectRad (t1d, p1d, qv1d, qc1d, nc1d, qi1d, ni1d, qs1d, &
rho(k) = 0.622*p1d(k)/(R*t1d(k)*(qv1d(k)+0.622))
rc(k) = MAX(R1, qc1d(k)*rho(k))
nc(k) = MAX(2., MIN(nc1d(k)*rho(k), Nt_c_max))
if (.NOT. (is_aerosol_aware .or. merra2_aerosol_aware)) nc(k) = Nt_c
if (.NOT. (is_aerosol_aware .or. merra2_aerosol_aware)) then
if( lsml == 1) then
nc(k) = Nt_c_l
else
nc(k) = Nt_c_o
endif
endif
if (rc(k).gt.R1 .and. nc(k).gt.R2) has_qc = .true.
ri(k) = MAX(R1, qi1d(k)*rho(k))
ni(k) = MAX(R2, ni1d(k)*rho(k))
Expand Down
Loading

0 comments on commit 312f63e

Please sign in to comment.