diff --git a/physics/module_mp_thompson.F90 b/physics/module_mp_thompson.F90 index b71c8f1f7..fc5a0a2fb 100644 --- a/physics/module_mp_thompson.F90 +++ b/physics/module_mp_thompson.F90 @@ -980,7 +980,7 @@ END SUBROUTINE thompson_init !> @{ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & nwfa, nifa, nwfa2d, nifa2d, & - aero_ind_fdb, tt, th, pii, & + tt, th, pii, & p, w, dz, dt_in, dt_inner, & sedi_semi, decfl, & RAINNC, RAINNCV, & @@ -994,7 +994,7 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & vt_dbz_wt, first_time_step, & re_cloud, re_ice, re_snow, & has_reqc, has_reqi, has_reqs, & - rand_perturb_on, & + aero_ind_fdb, rand_perturb_on, & kme_stoch, & rand_pert, spp_prt_list, spp_var_list, & spp_stddev_cutoff, n_var_spp, & @@ -1037,7 +1037,6 @@ 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 - LOGICAL, OPTIONAL, INTENT(IN):: aero_ind_fdb 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 @@ -1071,6 +1070,7 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & LOGICAL, INTENT (IN) :: reset_dBZ ! Extended diagnostics, array pointers only associated if ext_diag flag is .true. LOGICAL, INTENT (IN) :: ext_diag + LOGICAL, OPTIONAL, INTENT(IN):: aero_ind_fdb REAL, DIMENSION(:,:,:), INTENT(INOUT):: & !vts1, txri, txrc, & prw_vcdc, & @@ -1483,10 +1483,15 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & !.. Changed 13 May 2013 to fake emissions in which nwfa2d is aerosol !.. number tendency (number per kg per second). if (is_aerosol_aware) then - if ( .not. aero_ind_fdb) then - nwfa1d(kts) = nwfa1d(kts) + nwfa2d(i,j)*dt - nifa1d(kts) = nifa1d(kts) + nifa2d(i,j)*dt - endif + if ( PRESENT (aero_ind_fdb) ) then + if ( .not. aero_ind_fdb) then + nwfa1d(kts) = nwfa1d(kts) + nwfa2d(i,j)*dt + nifa1d(kts) = nifa1d(kts) + nifa2d(i,j)*dt + endif + else + nwfa1d(kts) = nwfa1d(kts) + nwfa2d(i,j)*dt + nifa1d(kts) = nifa1d(kts) + nifa2d(i,j)*dt + end if do k = kts, kte nc(i,k,j) = nc1d(k) @@ -2220,7 +2225,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & ni(k) = MAX(R2, ni1d(k)*rho(k)) if (ni(k).le. R2) then lami = cie(2)/5.E-6 - ni(k) = MIN(499.D3, cig(1)*oig2*ri(k)/am_i*lami**bm_i) + ni(k) = MIN(4999.D3, cig(1)*oig2*ri(k)/am_i*lami**bm_i) endif L_qi(k) = .true. lami = (am_i*cig(2)*oig1*ni(k)/ri(k))**obmi @@ -2228,7 +2233,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 - ni(k) = MIN(499.D3, cig(1)*oig2*ri(k)/am_i*lami**bm_i) + ni(k) = MIN(4999.D3, cig(1)*oig2*ri(k)/am_i*lami**bm_i) elseif (xDi.gt. 300.E-6) then lami = cie(2)/300.E-6 ni(k) = cig(1)*oig2*ri(k)/am_i*lami**bm_i @@ -2434,7 +2439,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & !+---+-----------------------------------------------------------------+ do k = kte, kts, -1 ygra1 = alog10(max(1.E-9, rg(k))) - zans1 = 3.0 + 2./7.*(ygra1+8.) + rand1 + zans1 = 3.4 + 2./7.*(ygra1+8.) + rand1 N0_exp = 10.**(zans1) N0_exp = MAX(DBLE(gonv_min), MIN(N0_exp, DBLE(gonv_max))) lam_exp = (N0_exp*am_g*cgg(1)/rg(k))**oge1 @@ -2462,12 +2467,9 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & do k = kts, kte !> - Rain self-collection follows Seifert, 1994 and drop break-up -!! follows Verlinde and Cotton, 1993. RAIN2M +!! follows Verlinde and Cotton, 1993. Updated after Saleeby et al 2022. RAIN2M if (L_qr(k) .and. mvd_r(k).gt. D0r) then -!-GT Ef_rr = 1.0 -!-GT if (mvd_r(k) .gt. 1500.0E-6) then - Ef_rr = 1.0 - EXP(2300.0*(mvd_r(k)-1950.0E-6)) -!-GT endif + Ef_rr = MAX(-0.1, 1.0 - EXP(2300.0*(mvd_r(k)-1950.0E-6))) pnr_rcr(k) = Ef_rr * 2.0*nr(k)*rr(k) endif @@ -2933,7 +2935,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & !> - Freezing of aqueous aerosols based on Koop et al (2001, Nature) xni = smo0(k)+ni(k) + (pni_rfz(k)+pni_wfz(k)+pni_inu(k))*dtsave - if ((is_aerosol_aware .or. merra2_aerosol_aware) .AND. homogIce .AND. (xni.le.499.E3) & + if ((is_aerosol_aware .or. merra2_aerosol_aware) .AND. homogIce .AND. (xni.le.4999.E3) & & .AND.(temp(k).lt.238).AND.(ssati(k).ge.0.4) ) then xnc = iceKoop(temp(k),qv(k),qvs(k),nwfa(k), dtsave) pni_iha(k) = xnc*odts @@ -3269,7 +3271,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 - xni = MIN(499.D3, cig(1)*oig2*xri/am_i*lami**bm_i) + 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 lami = cie(2)/300.E-6 @@ -3280,8 +3282,8 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & niten(k) = -ni1d(k)*odts endif xni=MAX(0.,(ni1d(k) + niten(k)*dtsave)*rho(k)) - if (xni.gt.499.E3) & - niten(k) = (499.E3-ni1d(k)*rho(k))*odts*orho + if (xni.gt.4999.E3) & + niten(k) = (4999.E3-ni1d(k)*rho(k))*odts*orho !> - Rain tendency qrten(k) = qrten(k) + (prr_wau(k) + prr_rcw(k) & @@ -3510,7 +3512,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & !+---+-----------------------------------------------------------------+ do k = kte, kts, -1 ygra1 = alog10(max(1.E-9, rg(k))) - zans1 = 3.0 + 2./7.*(ygra1+8.) + rand1 + zans1 = 3.4 + 2./7.*(ygra1+8.) + rand1 N0_exp = 10.**(zans1) N0_exp = MAX(DBLE(gonv_min), MIN(N0_exp, DBLE(gonv_max))) lam_exp = (N0_exp*am_g*cgg(1)/rg(k))**oge1 @@ -3959,7 +3961,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & pfll1(k) = pfll1(k) + sed_r(k)*DT*onstep(1) enddo - if (rr(kts).gt.R1*10.) & + if (rr(kts).gt.R1*1000.) & pptrain = pptrain + sed_r(kts)*DT*onstep(1) enddo else !if(.not. sedi_semi) @@ -4053,7 +4055,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & pfil1(k) = pfil1(k) + sed_i(k)*DT*onstep(2) enddo - if (ri(kts).gt.R1*10.) & + if (ri(kts).gt.R1*1000.) & pptice = pptice + sed_i(kts)*DT*onstep(2) enddo endif @@ -4082,7 +4084,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & pfil1(k) = pfil1(k) + sed_s(k)*DT*onstep(3) enddo - if (rs(kts).gt.R1*10.) & + if (rs(kts).gt.R1*1000.) & pptsnow = pptsnow + sed_s(kts)*DT*onstep(3) enddo endif @@ -4112,7 +4114,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & pfil1(k) = pfil1(k) + sed_g(k)*DT*onstep(4) enddo - if (rg(kts).gt.R1*10.) & + if (rg(kts).gt.R1*1000.) & pptgraul = pptgraul + sed_g(kts)*DT*onstep(4) enddo else ! if(.not. sedi_semi) then @@ -4137,7 +4139,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & vtg = 0. if (rg(k).gt. R1) then ygra1 = alog10(max(1.E-9, rg(k))) - zans1 = 3.0 + 2./7.*(ygra1+8.) + rand1 + zans1 = 3.4 + 2./7.*(ygra1+8.) + rand1 N0_exp = 10.**(zans1) N0_exp = MAX(DBLE(gonv_min), MIN(N0_exp, DBLE(gonv_max))) lam_exp = (N0_exp*am_g*cgg(1)/rg(k))**oge1 @@ -4237,7 +4239,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & lami = cie(2)/300.E-6 endif ni1d(k) = MIN(cig(1)*oig2*qi1d(k)/am_i*lami**bm_i, & - 499.D3/rho(k)) + 4999.D3/rho(k)) endif qr1d(k) = qr1d(k) + qrten(k)*DT nr1d(k) = MAX(R2/rho(k), nr1d(k) + nrten(k)*DT) @@ -5702,7 +5704,7 @@ end FUNCTION iceKoop !+---+-----------------------------------------------------------------+ !>\ingroup aathompson -!! Helper routine for Phillips et al (2008) ice nucleation. Trude +!! Helper routine for Phillips et al (2008) ice nucleation. REAL FUNCTION delta_p (yy, y1, y2, aa, bb) IMPLICIT NONE @@ -5745,6 +5747,7 @@ END FUNCTION delta_p !! schemes. Since only the smallest snowflakes should impact !! radiation, compute from first portion of complicated Field number !! 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) @@ -5860,6 +5863,7 @@ end subroutine calc_effectRad !! library of routines. The meltwater fraction is simply the amount !! of frozen species remaining from what initially existed at the !! melting level interface. + subroutine calc_refl10cm (qv1d, qc1d, qr1d, nr1d, qs1d, qg1d, & t1d, p1d, dBZ, rand1, kts, kte, ii, jj, melti, & vt_dBZ, first_time_step) @@ -5892,7 +5896,7 @@ subroutine calc_refl10cm (qv1d, qc1d, qr1d, nr1d, qs1d, qg1d, & REAL, DIMENSION(kts:kte):: ze_rain, ze_snow, ze_graupel DOUBLE PRECISION:: N0_exp, N0_min, lam_exp, lamr, lamg - REAL:: a_, b_, loga_, tc0 + REAL:: a_, b_, loga_, tc0, SR DOUBLE PRECISION:: fmelt_s, fmelt_g INTEGER:: i, k, k_0, kbot, n @@ -5915,7 +5919,7 @@ subroutine calc_refl10cm (qv1d, qc1d, qr1d, nr1d, qs1d, qg1d, & else do_vt_dBZ = .false. allow_wet_snow = .true. - allow_wet_graupel = .true. + allow_wet_graupel = .false. endif do k = kts, kte @@ -6031,7 +6035,7 @@ subroutine calc_refl10cm (qv1d, qc1d, qr1d, nr1d, qs1d, qg1d, & if (ANY(L_qg .eqv. .true.)) then do k = kte, kts, -1 ygra1 = alog10(max(1.E-9, rg(k))) - zans1 = 3.0 + 2./7.*(ygra1+8.) + rand1 + zans1 = 3.4 + 2./7.*(ygra1+8.) + rand1 N0_exp = 10.**(zans1) N0_exp = MAX(DBLE(gonv_min), MIN(N0_exp, DBLE(gonv_max))) lam_exp = (N0_exp*am_g*cgg(1)/rg(k))**oge1 @@ -6085,7 +6089,8 @@ subroutine calc_refl10cm (qv1d, qc1d, qr1d, nr1d, qs1d, qg1d, & !..Reflectivity contributed by melting snow if (allow_wet_snow .and. L_qs(k) .and. L_qs(k_0) ) then - fmelt_s = MAX(0.05d0, MIN(1.0d0-rs(k)/rs(k_0), 0.99d0)) + SR = MAX(0.01, MIN(1.0 - rs(k)/(rs(k) + rr(k)), 0.99)) + fmelt_s = DBLE(SR*SR) eta = 0.d0 oM3 = 1./smoc(k) M0 = (smob(k)*oM3) @@ -6108,7 +6113,8 @@ subroutine calc_refl10cm (qv1d, qc1d, qr1d, nr1d, qs1d, qg1d, & !..Reflectivity contributed by melting graupel if (allow_wet_graupel .and. L_qg(k) .and. L_qg(k_0) ) then - fmelt_g = MAX(0.05d0, MIN(1.0d0-rg(k)/rg(k_0), 0.99d0)) + SR = MAX(0.01, MIN(1.0 - rg(k)/(rg(k) + rr(k)), 0.99)) + fmelt_g = DBLE(SR*SR) eta = 0.d0 lamg = 1./ilamg(k) do n = 1, nrbins @@ -6214,7 +6220,7 @@ SUBROUTINE semi_lagrange_sedim(km,dzl,wwl,rql,precip,pfsan,dt,R1) real zsum,qsum,dim,dip,con1,fa1,fa2 real allold, decfl real dz(km), ww(km), qq(km) - real wi(km+1), zi(km+1), za(km+2) !hmhj + real wi(km+1), zi(km+1), za(km+2) real qn(km) real dza(km+1), qa(km+1), qmi(km+1), qpi(km+1) real net_flx(km) @@ -6265,12 +6271,12 @@ SUBROUTINE semi_lagrange_sedim(km,dzl,wwl,rql,precip,pfsan,dt,R1) enddo wi(km) = 0.5*(ww(km)+ww(km-1)) wi(km+1) = ww(km) -! + ! terminate of top of raingroup do k=2,km if( ww(k).eq.0.0 ) wi(k)=ww(k-1) enddo -! + ! diffusivity of wi con1 = 0.05 do k=km,1,-1 @@ -6283,18 +6289,18 @@ SUBROUTINE semi_lagrange_sedim(km,dzl,wwl,rql,precip,pfsan,dt,R1) do k=1,km+1 za(k) = zi(k) - wi(k)*dt enddo - za(km+2) = zi(km+1) !hmhj -! - do k=1,km+1 !hmhj + za(km+2) = zi(km+1) + + do k=1,km+1 dza(k) = za(k+1)-za(k) enddo -! + ! computer deformation at arrival point do k=1,km qa(k) = qq(k)*dz(k)/dza(k) enddo qa(km+1) = 0.0 -! + ! estimate values at arrival cell interface with monotone do k=2,km dip=(qa(k+1)-qa(k))/(dza(k+1)+dza(k)) @@ -6315,7 +6321,7 @@ SUBROUTINE semi_lagrange_sedim(km,dzl,wwl,rql,precip,pfsan,dt,R1) qmi(1)=qa(1) qmi(km+1)=qa(km+1) qpi(km+1)=qa(km+1) -! + ! interpolation to regular point qn = 0.0 kb=1 @@ -6335,7 +6341,7 @@ SUBROUTINE semi_lagrange_sedim(km,dzl,wwl,rql,precip,pfsan,dt,R1) cycle find_kb endif enddo find_kb - find_kt : do kk=kt,km+2 !hmhj + find_kt : do kk=kt,km+2 if( zi(k+1).le.za(kk) ) then kt = kk exit find_kt @@ -6378,24 +6384,21 @@ SUBROUTINE semi_lagrange_sedim(km,dzl,wwl,rql,precip,pfsan,dt,R1) endif cycle intp endif -! + enddo intp -! + ! rain out sum_precip: do k=1,km if( za(k).lt.0.0 .and. za(k+1).le.0.0 ) then -!hmhj precip = precip + qa(k)*dza(k) net_flx(k) = qa(k)*dza(k) cycle sum_precip else if ( za(k).lt.0.0 .and. za(k+1).gt.0.0 ) then -!hmhj -!hmhj precip(i) = precip(i) + qa(k)*(0.0-za(k)) - th = (0.0-za(k))/dza(k) !hmhj - th2 = th*th !hmhj - qqd = 0.5*(qpi(k)-qmi(k)) !hmhj - qqh = qqd*th2+qmi(k)*th !hmhj - precip = precip + qqh*dza(k) !hmhj + th = (0.0-za(k))/dza(k) + th2 = th*th + qqd = 0.5*(qpi(k)-qmi(k)) + qqh = qqd*th2+qmi(k)*th + precip = precip + qqh*dza(k) net_flx(k) = qqh*dza(k) exit sum_precip endif @@ -6413,11 +6416,9 @@ SUBROUTINE semi_lagrange_sedim(km,dzl,wwl,rql,precip,pfsan,dt,R1) ! ! replace the new values rql(:) = max(qn(:),R1) -! -! ---------------------------------- -! + END SUBROUTINE semi_lagrange_sedim -!+---+-----------------------------------------------------------------+ + !+---+-----------------------------------------------------------------+ !+---+-----------------------------------------------------------------+ END MODULE module_mp_thompson diff --git a/physics/mp_thompson.F90 b/physics/mp_thompson.F90 index 80762a5f1..ffe1a03d6 100644 --- a/physics/mp_thompson.F90 +++ b/physics/mp_thompson.F90 @@ -686,7 +686,6 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & if (is_aerosol_aware .or. merra2_aerosol_aware) then call mp_gt_driver(qv=qv, qc=qc, qr=qr, qi=qi, qs=qs, qg=qg, ni=ni, nr=nr, & nc=nc, nwfa=nwfa, nifa=nifa, nwfa2d=nwfa2d, nifa2d=nifa2d, & - aero_ind_fdb=aero_ind_fdb, & tt=tgrs, p=prsl, w=w, dz=dz, dt_in=dtstep, dt_inner=dt_inner, & sedi_semi=sedi_semi, decfl=decfl, & rainnc=rain_mp, rainncv=delta_rain_mp, & @@ -696,7 +695,8 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & refl_10cm=refl_10cm, & diagflag=diagflag, do_radar_ref=do_radar_ref_mp, & has_reqc=has_reqc, has_reqi=has_reqi, has_reqs=has_reqs, & - rand_perturb_on=spp_mp_opt, kme_stoch=kme_stoch, & + aero_ind_fdb=aero_ind_fdb, rand_perturb_on=spp_mp_opt, & + kme_stoch=kme_stoch, & rand_pert=spp_wts_mp, spp_var_list=spp_var_list, & spp_prt_list=spp_prt_list, n_var_spp=n_var_spp, & spp_stddev_cutoff=spp_stddev_cutoff, & @@ -887,6 +887,29 @@ subroutine get_niwfa(aerfld, nifa, nwfa, ncol, nlev) ! mass. NIFA is mainly summarized over five dust bins and NWFA over the ! other 10 bins. The parameters besides each bins are carefully tuned ! for a good performance of the scheme. + ! + ! The fields for the last index of the aerfld array + ! are specified as below. + ! 1: dust bin 1, 0.1 to 1.0 micrometers + ! 2: dust bin 2, 1.0 to 1.8 micrometers + ! 3: dust bin 3, 1.8 to 3.0 micrometers + ! 4: dust bin 4, 3.0 to 6.0 micrometers + ! 5: dust bin 5, 6.0 to 10.0 micrometers + ! 6: sea salt bin 1, 0.03 to 0.1 micrometers + ! 7: sea salt bin 2, 0.1 to 0.5 micrometers + ! 8: sea salt bin 3, 0.5 to 1.5 micrometers + ! 9: sea salt bin 4, 1.5 to 5.0 micrometers + ! 10: sea salt bin 5, 5.0 to 10.0 micrometers + ! 11: Sulfate, 0.35 (mean) micrometers + ! 15: water-friendly organic carbon, 0.35 (mean) micrometers + ! + ! Bin densities are as follows: + ! 1: dust bin 1: 2500 kg/m2 + ! 2-5: dust bin 2-5: 2650 kg/m2 + ! 6-10: sea salt bins 6-10: 2200 kg/m2 + ! 11: sulfate: 1700 kg/m2 + ! 15: organic carbon: 1800 kg/m2 + implicit none integer, intent(in)::ncol, nlev real (kind=kind_phys), dimension(:,:,:), intent(in) :: aerfld diff --git a/physics/radiation_clouds.f b/physics/radiation_clouds.f index edfa94439..f7e017849 100644 --- a/physics/radiation_clouds.f +++ b/physics/radiation_clouds.f @@ -2487,7 +2487,7 @@ subroutine progcld_thompson & do i = 1, IX cwp(i,k) = max(0.0, clw(i,k,ntcw) * dz(i,k)*1.E6) crp(i,k) = 0.0 - snow_mass_factor = 0.99 + snow_mass_factor = 0.90 cip(i,k) = max(0.0, (clw(i,k,ntiw) & & + (1.0-snow_mass_factor)*clw(i,k,ntsw))*dz(i,k)*1.E6) if (re_snow(i,k) .gt. snow_max_radius)then @@ -3378,7 +3378,7 @@ SUBROUTINE cal_cldfra3(CLDFRA, qv, qc, qi, qs, dz, & REAL:: RH_00L, RH_00O, RH_00 REAL:: entrmnt=0.5 INTEGER:: k - REAL:: TC, qvsi, qvsw, RHUM, delz + REAL:: TC, qvsi, qvsw, RHUM, delz, var_temp REAL, DIMENSION(kts:kte):: qvs, rh, rhoa integer:: ndebug = 0 @@ -3438,7 +3438,8 @@ SUBROUTINE cal_cldfra3(CLDFRA, qv, qc, qi, qs, dz, & CLDFRA(K) = 1.0 elseif (((qc(k)+qi(k)).gt.1.E-10) .and. & & ((qc(k)+qi(k)).lt.1.E-6)) then - CLDFRA(K) = MIN(0.99, 0.1*(11.0 + log10(qc(k)+qi(k)))) + var_temp = MIN(0.99, 0.1*(11.0 + log10(qc(k)+qi(k)))) + CLDFRA(K) = var_temp*var_temp else IF ((XLAND-1.5).GT.0.) THEN !--- Ocean @@ -3448,24 +3449,27 @@ SUBROUTINE cal_cldfra3(CLDFRA, qv, qc, qi, qs, dz, & ENDIF tc = MAX(-80.0, t(k) - 273.15) - if (tc .lt. -12.0) RH_00 = RH_00L + if (tc .lt. -24.0) RH_00 = RH_00L if (tc .gt. 20.0) then CLDFRA(K) = 0.0 - elseif (tc .ge. -12.0) then + elseif (tc .ge. -24.0) then RHUM = MIN(rh(k), 1.0) - CLDFRA(K) = MAX(0., 1.0-SQRT((1.001-RHUM)/(1.001-RH_00))) + var_temp = SQRT(SQRT((1.001-RHUM)/(1.001-RH_00))) + CLDFRA(K) = MAX(0., 1.0-var_temp) else if (max_relh.gt.1.12 .or. (.NOT.(modify_qvapor)) ) then !..For HRRR model, the following look OK. RHUM = MIN(rh(k), 1.45) - RH_00 = RH_00 + (1.45-RH_00)*(-12.0-tc)/(-12.0+85.) - CLDFRA(K) = MAX(0.,1.0-SQRT((1.46-RHUM)/(1.46-RH_00))) + RH_00 = RH_00 + (1.45-RH_00)*(-24.0-tc)/(-24.0+85.) + var_temp = SQRT(SQRT((1.46-RHUM)/(1.46-RH_00))) + CLDFRA(K) = MAX(0., 1.0-var_temp) else !..but for the GFS model, RH is way lower. RHUM = MIN(rh(k), 1.05) - RH_00 = RH_00 + (1.05-RH_00)*(-12.0-tc)/(-12.0+85.) - CLDFRA(K) = MAX(0.,1.0-SQRT((1.06-RHUM)/(1.06-RH_00))) + RH_00 = RH_00 + (1.05-RH_00)*(-24.0-tc)/(-24.0+85.) + var_temp = SQRT(SQRT((1.06-RHUM)/(1.06-RH_00))) + CLDFRA(K) = MAX(0., 1.0-var_temp) endif endif if (CLDFRA(K).gt.0.) CLDFRA(K)=MAX(0.01,MIN(CLDFRA(K),0.99)) @@ -3784,6 +3788,8 @@ SUBROUTINE adjust_cloudFinal(cfr, qc, qi, Rho,dz, kts,kte) END SUBROUTINE adjust_cloudFinal +!+---+-----------------------------------------------------------------+ + subroutine cloud_fraction_XuRandall & & ( IX, NLAY, plyr, clwf, rhly, qstl, & ! --- inputs & cldtot ) & ! --- outputs @@ -3808,7 +3814,6 @@ subroutine cloud_fraction_XuRandall & do k = 1, NLAY do i = 1, IX clwt = 1.0e-6 * (plyr(i,k)*0.001) -! clwt = 2.0e-6 * (plyr(i,k)*0.001) if (clwf(i,k) > clwt) then @@ -3818,8 +3823,6 @@ subroutine cloud_fraction_XuRandall & tem1 = min(max(sqrt(sqrt(onemrh*qstl(i,k))),0.0001),1.0) tem1 = 2000.0 / tem1 -! tem1 = 1000.0 / tem1 - value = max( min( tem1*(clwf(i,k)-clwm), 50.0 ), 0.0 ) tem2 = sqrt( sqrt(rhly(i,k)) ) @@ -3830,6 +3833,8 @@ subroutine cloud_fraction_XuRandall & end subroutine cloud_fraction_XuRandall +!+---+-----------------------------------------------------------------+ + subroutine cloud_fraction_mass_flx_1 & & ( IX, NLAY, lmfdeep2, xrc3, plyr, clwf, rhly, qstl, & ! --- inputs & cldtot ) & ! --- outputs @@ -3879,6 +3884,8 @@ subroutine cloud_fraction_mass_flx_1 & end subroutine cloud_fraction_mass_flx_1 +!+---+-----------------------------------------------------------------+ + subroutine cloud_fraction_mass_flx_2 & & ( IX, NLAY, lmfdeep2, xrc3, plyr, clwf, rhly, qstl, & ! --- inputs & cldtot ) & ! --- outputs