#if ( RWORDSIZE == 4 ) # define VREC vsrec # define VSQRT vssqrt #else # define VREC vrec # define VSQRT vsqrt #endif !------------------------------------------------------------------------------- module module_mp_wdm6 !------------------------------------------------------------------------------- ! use module_utility, only: WRFU_Clock, WRFU_Alarm use module_domain, only: HISTORY_ALARM, Is_alarm_tstep use module_mp_radar ! ! !------------------------------------------------------------------------------- ! ! constant : ! ! dtcldcr - maximum time step for minor loops ! n0r - intercept parameter rain ! n0s - intercept parameter snow ! n0g - intercept parameter graupel ! n0smax - maximum n0s (t=-90C unlimited) ! alpha - 0.122 exponen factor for n0s ! avtr, bvtr - a constant for terminal velocity of rain ! avts, bvts - a constant for terminal velocity of snow ! avtg, bvtg - a constant for terminal velocity of graupel ! lamdarmax - limited maximum value for slope parameter of rain ! lamdarmin - limited minimum value for slope parameter of rain ! lamdasmax - limited maximum value for slope parameter of snow ! lamdagmax - limited maximum value for slope parameter of graupel ! r0 - 8 microm in contrast to 10 micro m ! peaut - collection efficiency ! xncr - maritime cloud in contrast to 3.e8 in tc80 ! xmyu - the dynamic viscosity kg m-1 s-1 ! dicon - constant for the cloud-ice diamter ! dimax - limited maximum value for the cloud-ice diamter ! pfrz1, pfrz2 - constant in Biggs freezing ! qcrmin - minimum values for qr and qs ! ncmin - minimum value for nc ! nrmin - minimum value for nr ! eacrc - now/cloud-water collection efficiency ! qs0 - threshold amount for aggretion to occur ! satmax - maximum saturation value for CCN activation ! 1.008 for maritime /1.0048 for conti ! actk - parameter for the CCN activation ! actr - radius of activated CCN drops ! ncrk1, ncrk2 - Long's collection kernel coefficient ! di100, di600, di2000 - parameter related with accretion and collection ! of cloud drops ! di82 - dimater related with raindrops evaporation ! di15 - auto conversion takes place beyond this diameter ! !------------------------------------------------------------------------------- real, parameter, private :: dtcldcr = 120. real, parameter, private :: n0r = 8.e6 real, parameter, private :: n0s = 2.e6 ! real, parameter, private :: n0g = 4.e6 real, parameter, private :: n0smax = 1.e11 real, parameter, private :: dens = 100.0 real, parameter, private :: alpha = .12 real, parameter, private :: avtr = 841.9 real, parameter, private :: bvtr = 0.8 real, parameter, private :: avts = 11.72 real, parameter, private :: bvts = .41 ! real, parameter, private :: avtg = 330. ! real, parameter, private :: bvtg = 0.8 real, parameter, private :: lamdacmax = 5.e5 real, parameter, private :: lamdacmin = 2.e4 real, parameter, private :: lamdarmax = 5.e4 real, parameter, private :: lamdarmin = 2.e3 real, parameter, private :: lamdasmax = 1.e5 ! real, parameter, private :: lamdagmax = 6.e4 real, parameter, private :: r0 = .8e-5 real, parameter, private :: peaut = .55 real, parameter, private :: xncr = 3.e8 real, parameter, private :: xncr0 = 5.e7 real, parameter, private :: xncr1 = 5.e8 real, parameter, private :: xmyu = 1.718e-5 real, parameter, private :: dicon = 11.9 real, parameter, private :: dimax = 500.e-6 real, parameter, private :: pfrz1 = 100. real, parameter, private :: pfrz2 = 0.66 real, parameter, private :: qcrmin = 1.e-9 real, parameter, private :: ncmin = 1.e1 real, parameter, private :: nrmin = 1.e-2 real, parameter, private :: eacrc = 1.0 real, parameter, private :: qs0 = 6.e-4 real, parameter, private :: satmax = 1.0048 real, parameter, private :: actk = 0.6 real, parameter, private :: actr = 1.5 real, parameter, private :: ncrk1 = 3.03e3 real, parameter, private :: ncrk2 = 2.59e15 real, parameter, private :: di100 = 1.e-4 real, parameter, private :: di600 = 6.e-4 real, parameter, private :: di2000 = 2000.e-6 real, parameter, private :: di82 = 82.e-6 real, parameter, private :: di15 = 15.e-6 real, save :: & qc0,qc1,qck1,pidnc,bvtr1,bvtr2,bvtr3,bvtr4,bvtr5, & bvtr6,bvtr7, bvtr2o5,bvtr3o5, & g1pbr,g2pbr,g3pbr,g4pbr,g5pbr,g6pbr,g7pbr, & g5pbro2,g7pbro2,pi, & pvtr,pvtrn,eacrr,pacrr,pidn0r,pidnr, & precr1,precr2,xmmax,roqimax,bvts1,bvts2, & bvts3,bvts4,g1pbs,g3pbs,g4pbs,g5pbso2, & pvts,pacrs,precs1,precs2,pidn0s,xlv1,pacrc, & bvtg1,bvtg2,bvtg3,bvtg4,g1pbg,g3pbg,g4pbg, & g5pbgo2,pvtg,pacrg,precg1,precg2,pidn0g, & n0g,avtg,bvtg,deng,lamdagmax, & rslopecmax,rslopec2max,rslopec3max, & rslopermax,rslopesmax,rslopegmax, & rsloperbmax,rslopesbmax,rslopegbmax, & rsloper2max,rslopes2max,rslopeg2max, & rsloper3max,rslopes3max,rslopeg3max ! contains !=================================================================== ! subroutine wdm6(th, q, qc, qr, qi, qs, qg, & nn, nc, nr, & den, pii, p, delz, & delt,g, cpd, cpv, ccn0, rd, rv, t0c, & ep1, ep2, qmin, & xls, xlv0, xlf0, den0, denr, & cliq,cice,psat, & xland, & rain, rainncv, & snow, snowncv, & sr, & refl_10cm, diagflag, do_radar_ref, & graupel, graupelncv, & itimestep, & has_reqc, has_reqi, has_reqs, & ! for radiation re_cloud, re_ice, re_snow, & ! for radiation ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte & ) !------------------------------------------------------------------- implicit none !------------------------------------------------------------------- integer , intent(in ) :: itimestep integer , intent(in ) :: & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte real , intent(in ) :: delt, g, ccn0, & rd, rv, t0c, & cpd, cpv, & den0, qmin, & ep1, ep2, xls, & xlv0, xlf0, & cliq, cice, & psat, denr real, dimension(ims:ime,jms:jme) , intent(in ) :: xland real, dimension(ims:ime,kms:kme,jms:jme) , intent(in ) :: den real, dimension(ims:ime,kms:kme,jms:jme) , intent(in ) :: pii real, dimension(ims:ime,kms:kme,jms:jme) , intent(in ) :: p real, dimension(ims:ime,kms:kme,jms:jme) , intent(in ) :: delz real, dimension(ims:ime,kms:kme,jms:jme) , intent(inout) :: th real, dimension(ims:ime,kms:kme,jms:jme) , intent(inout) :: q real, dimension(ims:ime,kms:kme,jms:jme) , intent(inout) :: qc real, dimension(ims:ime,kms:kme,jms:jme) , intent(inout) :: qi real, dimension(ims:ime,kms:kme,jms:jme) , intent(inout) :: qr real, dimension(ims:ime,kms:kme,jms:jme) , intent(inout) :: qs real, dimension(ims:ime,kms:kme,jms:jme) , intent(inout) :: qg real, dimension(ims:ime,kms:kme,jms:jme) , intent(inout) :: nn real, dimension(ims:ime,kms:kme,jms:jme) , intent(inout) :: nc real, dimension(ims:ime,kms:kme,jms:jme) , intent(inout) :: nr real, dimension(ims:ime,jms:jme) , intent(inout) :: rain real, dimension(ims:ime,jms:jme) , intent(inout) :: rainncv real, dimension(ims:ime,jms:jme) , intent(inout) :: sr real, dimension(ims:ime,jms:jme), optional, intent(inout) :: snow real, dimension(ims:ime,jms:jme), optional, intent(inout) :: snowncv real, dimension(ims:ime,jms:jme), optional, intent(inout) :: graupel real, dimension(ims:ime,jms:jme), optional, intent(inout) :: graupelncv ! ! for radiation connecting ! integer , intent(in ) :: has_reqc integer , intent(in ) :: has_reqi integer , intent(in ) :: has_reqs real, dimension(ims:ime,kms:kme,jms:jme), intent(inout) :: re_cloud real, dimension(ims:ime,kms:kme,jms:jme), intent(inout) :: re_ice real, dimension(ims:ime,kms:kme,jms:jme), intent(inout) :: re_snow ! ! for reflectivity ! real, dimension(ims:ime,kms:kme,jms:jme), intent(inout) :: refl_10cm logical , optional, intent(in ) :: diagflag integer , optional, intent(in ) :: do_radar_ref ! ! local variables ! real, dimension(its:ite,kts:kte) :: t real, dimension(its:ite,kts:kte,2) :: qci real, dimension(its:ite,kts:kte,3) :: qrs, ncr integer :: i,j,k ! ! for reflectivity ! real, dimension(kts:kte) :: qv1d real, dimension(kts:kte) :: t1d real, dimension(kts:kte) :: p1d real, dimension(kts:kte) :: qr1d real, dimension(kts:kte) :: nr1d real, dimension(kts:kte) :: qs1d real, dimension(kts:kte) :: qg1d real, dimension(kts:kte) :: dBZ ! ! to calculate effective radius for radiation ! real, dimension(kts:kte) :: qc1d, nc1d, den1d real, dimension(kts:kte) :: qi1d real, dimension(kts:kte) :: re_qc, re_qi, re_qs ! if (itimestep .eq. 1) then do j = jms,jme do k = kms,kme do i = ims,ime nn(i,k,j) = ccn0 enddo enddo enddo endif ! do j = jts,jte do k = kts,kte do i = its,ite t(i,k) = th(i,k,j)*pii(i,k,j) qci(i,k,1) = qc(i,k,j) qci(i,k,2) = qi(i,k,j) qrs(i,k,1) = qr(i,k,j) qrs(i,k,2) = qs(i,k,j) qrs(i,k,3) = qg(i,k,j) ncr(i,k,1) = nn(i,k,j) ncr(i,k,2) = nc(i,k,j) ncr(i,k,3) = nr(i,k,j) enddo enddo ! ! Sending array starting locations of optional variables may cause ! troubles, so we explicitly change the call. ! call wdm62D(t, q(ims,kms,j), qci, qrs, ncr & ,den(ims,kms,j) & ,p(ims,kms,j), delz(ims,kms,j) & ,delt,g, cpd, cpv, ccn0, rd, rv, t0c & ,ep1, ep2, qmin & ,xls, xlv0, xlf0, den0, denr & ,cliq,cice,psat & ,j & ,xland(ims,j) & ,rain(ims,j),rainncv(ims,j) & ,sr(ims,j) & ,ids,ide, jds,jde, kds,kde & ,ims,ime, jms,jme, kms,kme & ,its,ite, jts,jte, kts,kte & ,snow(ims,j),snowncv(ims,j) & ,graupel(ims,j),graupelncv(ims,j) & ) ! do k = kts,kte do i = its,ite th(i,k,j) = t(i,k)/pii(i,k,j) qc(i,k,j) = qci(i,k,1) qi(i,k,j) = qci(i,k,2) qr(i,k,j) = qrs(i,k,1) qs(i,k,j) = qrs(i,k,2) qg(i,k,j) = qrs(i,k,3) nn(i,k,j) = ncr(i,k,1) nc(i,k,j) = ncr(i,k,2) nr(i,k,j) = ncr(i,k,3) enddo enddo ! if ( present (diagflag) ) then if (diagflag .and. do_radar_ref == 1) then do i = its,ite do k = kts,kte t1d(k) = th(i,k,j)*pii(i,k,j) p1d(k) = p(i,k,j) qv1d(k) = q(i,k,j) qr1d(k) = qr(i,k,j) nr1d(k) = nr(i,k,j) qs1d(k) = qs(i,k,j) qg1d(k) = qg(i,k,j) enddo call refl10cm_wdm6 (qv1d, qr1d, nr1d, qs1d, qg1d, & t1d, p1d, dBZ, kts, kte, i, j) do k = kts, kte refl_10cm(i,k,j) = max(-35., dBZ(k)) enddo enddo endif endif ! ! calculate effective radius of cloud, ice, and snow ! if (has_reqc.ne.0 .and. has_reqi.ne.0 .and. has_reqs.ne.0) then do i = its,ite do k = kts,kte re_qc(k) = 2.51e-6 re_qi(k) = 10.01e-6 re_qs(k) = 25.e-6 ! t1d(k) = th(i,k,j)*pii(i,k,j) den1d(k)= den(i,k,j) qc1d(k) = qc(i,k,j) qi1d(k) = qi(i,k,j) qs1d(k) = qs(i,k,j) nc1d(k) = nc(i,k,j) enddo ! call effectRad_wdm6(t1d, qc1d, nc1d, qi1d, qs1d, den1d, & qmin, t0c, re_qc, re_qi, re_qs, & kts, kte, i, j) ! do k = kts,kte re_cloud(i,k,j) = max(2.51e-6, min(re_qc(k), 50.e-6)) re_ice(i,k,j) = max(10.01e-6, min(re_qi(k), 125.e-6)) re_snow(i,k,j) = max(25.e-6, min(re_qs(k), 999.e-6)) enddo enddo endif enddo ! end subroutine wdm6 ! !------------------------------------------------------------------------------ ! subroutine wdm62D(t, q, qci, qrs, ncr, den, p, delz & ,delt,g, cpd, cpv, ccn0, rd, rv, t0c & ,ep1, ep2, qmin & ,xls, xlv0, xlf0, den0, denr & ,cliq,cice,psat & ,lat & ,slmsk & ,rain,rainncv & ,sr & ,ids,ide, jds,jde, kds,kde & ,ims,ime, jms,jme, kms,kme & ,its,ite, jts,jte, kts,kte & ,snow,snowncv & ,graupel,graupelncv & ) !------------------------------------------------------------------------------ implicit none !------------------------------------------------------------------------------- ! ! This code is a WRF double-moment 6-class GRAUPEL phase ! microphyiscs scheme (wdm6). The WDM microphysics scheme predicts ! number concentrations for warm rain species including clouds and ! rain. cloud condensation nuclei (ccn) is also predicted. ! The cold rain species including ice, snow, graupel follow the ! WRF single-moment 6-class microphysics (wsm6, Hong and Lim 2006) ! in which theoretical background for WSM ice phase microphysics is ! based on Hong et al. (2004). A new mixed-phase terminal velocity ! for precipitating ice is introduced in WSM6 (Dudhia et al. 2008). ! The WDM scheme is described in Lim and Hong (2010). ! ! wdm6 cloud scheme ! ! Coded by Kyo-Sun Lim and Song-You Hong (Yonsei Univ.) Fall 2008 ! ! Implemented by Kyo-Sun Lim and Jimy Dudhia (NCAR) Winter 2008 ! ! further modifications : ! semi-lagrangian sedimentation (JH,2010),hong, aug 2009 ! ==> higher accuracy and efficient at lower resolutions ! reflectivity computation from greg thompson, lim, jun 2011 ! ==> only dignostic, but with removal of too large drops ! add hail option from afwa, aug 2014 ! ==> switch graupel or hail by changing no, den, fall vel. ! effetive radius of hydrometeors, bae from kiaps, jan 2015 ! ==> consistency in solar insolation of rrtmg radiation ! bug fix in melting terms, bae from kiaps, nov 2015 ! ==> density of air is divided, which has not been ! bug fix in sedimentation of rain, bae from kiaps, mar 2017 ! change autoconversion rate to equation, bae from kiaps, oct 2017 ! nccn > 100 cm-3, bae from kiaps, oct 2017 ! bug fix hydrometeros update before condensation process of rain, ! bae from kiaps, oct2017 ! ! history log : ! ! Reference) Lim and Hong (LH, 2010) Mon. Wea. Rev. ! Juang and Hong (JH, 2010) Mon. Wea. Rev. ! Hong, Dudhia, Chen (HDC, 2004) Mon. Wea. Rev. ! Hong and Lim (HL, 2006) J. Korean Meteor. Soc. ! Cohard and Pinty (CP, 2000) Quart. J. Roy. Meteor. Soc. ! Khairoutdinov and Kogan (KK, 2000) Mon. Wea. Rev. ! Dudhia, Hong and Lim (DHL, 2008) J. Meteor. Soc. Japan ! ! Lin, Farley, Orville (LFO, 1983) J. Appl. Meteor. ! Rutledge, Hobbs (RH83, 1983) J. Atmos. Sci. ! Rutledge, Hobbs (RH84, 1984) J. Atmos. Sci. ! !------------------------------------------------------------------------------- ! ! input : ! deltim - timestep ! g, cpd, cpv, t0c, den0, - constant ! rd, rv, ep1, ep2, qmin, ! xls, xlv0, xlf0, denr, ! cliq, cice, psat ! ids, ide, jds, jde, kds, kde - dimension ! ims, ime, jms, jme, kms, kme ! its, ite, jts, jte, kts, kte ! ncloud - number of hydrometeor ! p - pressure, Pa ! delz - depth of model layer, m ! * with chemistry ! slmsk - land/sea mask, 0: sea, 1: land ! so4 - mixing ratio of sulfate ! ss - mixing ratio of sea salt ! oc - mixing ratio of organic carbon ! ! inout : ! t1 - temperautre ! q1 - specific humidity ! q2 - mixing ratio of cloud water, rain, ice, ! snow, and graupel ! qc, qr, qi, qs, qg ! n2 - number concentration of ccn, cloud droplet, ! and raindrop ! nn, nc, nr ! rain, rainncv - precipitation ! sr - ratio of snow to rain ! ! all units are in m.k.s. and source/sink terms in kg kg-1 s-1. !------------------------------------------------------------------------------- integer , intent(in ) :: & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte integer , intent(in ) :: lat real , intent(in ) :: delt real , intent(in ) :: g, cpd, cpv, t0c, & den0, rd, rv, & ep1, ep2, qmin, & xls, xlv0, xlf0, & cliq, cice, psat, & denr real , intent(in ) :: ccn0 real, dimension(ims:ime) , intent(in ) :: slmsk real, dimension(ims:ime,kms:kme) , intent(in ) :: p real, dimension(ims:ime,kms:kme) , intent(in ) :: delz real, dimension(ims:ime,kms:kme) , intent(in ) :: den real, dimension(its:ite,kts:kte) , intent(inout) :: t real, dimension(its:ite,kts:kte,2) , intent(inout) :: qci real, dimension(its:ite,kts:kte,3) , intent(inout) :: qrs real, dimension(its:ite,kts:kte,3) , intent(inout) :: ncr real, dimension(ims:ime,kms:kme) , intent(inout) :: q real, dimension(ims:ime) , intent(inout) :: rain real, dimension(ims:ime) , intent(inout) :: rainncv real, dimension(ims:ime) , intent(inout) :: sr real, dimension(ims:ime), optional , intent(inout) :: snow real, dimension(ims:ime), optional , intent(inout) :: snowncv real, dimension(ims:ime), optional , intent(inout) :: graupel real, dimension(ims:ime), optional , intent(inout) :: graupelncv ! ! local variables ! real, dimension(its:ite,kts:kte) :: dend real, dimension(its:ite,kts:kte) :: qcr real, dimension(its:ite,kts:kte,3) :: rh real, dimension(its:ite,kts:kte,3) :: qs real, dimension(its:ite,kts:kte,3) :: rslope, rslope2, rslope3, rslopeb real, dimension(its:ite,kts:kte,3) :: falk, fall real, dimension(its:ite,kts:kte,3) :: work1 real, dimension(its:ite,kts:kte,3) :: qrs_tmp real, dimension(its:ite,kts:kte,2) :: avedia real, dimension(its:ite,kts:kte) :: rslopec, rslopec2,rslopec3 real, dimension(its:ite,kts:kte) :: workn, falln, falkn real, dimension(its:ite,kts:kte) :: worka, workr real, dimension(its:ite,kts:kte) :: den_tmp, delz_tmp, ncr_tmp real, dimension(its:ite,kts:kte) :: lamdr_tmp real, dimension(its:ite,kts:kte) :: lamdc_tmp real, dimension(its:ite,kts:kte) :: falkc, work1c, work2c, fallc real, dimension(its:ite,kts:kte) :: dqr,dnr real, dimension(its:ite,kts:kte) :: pcact, prevp, psdep, pgdep, praut, & psaut, pgaut, pracw, psacw, pgacw, & pgacr, pgacs, psaci, pgmlt, praci, & piacr, pracs, psacr, pgaci, pseml, & pgeml, paacw, pigen, pidep, pcond, & psmlt, psevp, pgevp real, dimension(its:ite,kts:kte) :: nraut, nracw, ncevp, nccol, nrcol, & nsacw, ngacw, niacr, nsacr, ngacr, & naacw, nseml, ngeml, ncact real, dimension(its:ite,kts:kte) :: xl, cpm, work2, denfac, xni, n0sfac, & qsum real, dimension(its:ite,kts:kte) :: denqrs1, denqrs2, denqrs3 real, dimension(its:ite,kts:kte) :: denqr1, denncr3, denqci real, dimension(its:ite) :: delqrs1, delqrs2, delqrs3, delncr3 real, dimension(its:ite) :: delqi real, dimension(its:ite) :: tstepsnow, tstepgraup real :: gfac, sfac ! ! top level for computing of cloud microphysical process ! integer :: kte_in ! ! variables for optimization ! real, dimension(its:ite) :: tvec1 integer, dimension(its:ite) :: mnstep, numndt integer, dimension(its:ite) :: mstep, numdt logical, dimension(its:ite) :: flgcld real :: temp real :: vt2ave real :: holdc, holdci real :: cpmcal, xlcal, lamdac, diffus, & viscos, xka, venfac, conden, diffac, & x, y, z, a, b, c, d, e, & ndt, qdt, holdrr, holdrs, holdrg, & supcol, supcolt, pvt, coeres, supsat, & dtcld, xmi, eacrs, satdt, qimax, & diameter, xni0, roqi0, & fallsum, fallsum_qsi, fallsum_qg, & vt2i, vt2r, vt2s, vt2g, acrfac, & egs, egi, coecol, & xlwork2, factor, source, value, & nfrzdtr, nfrzdtc, & taucon, lencon, lenconcr, & xlf, pfrzdtc, pfrzdtr, supice, & alpha2, delta2, delta3 ! integer :: i, j, k, mstepmax, & iprt, latd, lond, loop, loops, ifsat, & n, idim, kdim ! ! Temporaries used for inlining fpvs function ! real :: dldti, xb, xai, tr, xbi, xa, hvap, & cvap, hsub, dldt, ttp !------------------------------------------------------------------------------- ! ! compute internal functions ! cpmcal(x) = cpd*(1.-max(x,qmin))+max(x,qmin)*cpv xlcal(x) = xlv0-xlv1*(x-t0c) ! ! size distributions: (x=mixing ratio, y=air density): ! valid for mixing ratio > 1.e-9 kg/kg. ! ! optimizatin : A**B => exp(log(A)*(B)) ! lamdac(x,y,z)= exp(log(((pidnc*z)/(x*y)))*((.33333333))) ! ! diffus: diffusion coefficient of the water vapor ! viscos: kinematic viscosity(m2s-1) ! diffus(x,y) = 8.794e-5 * exp(log(x)*(1.81)) / y ! 8.794e-5*x**1.81/y viscos(x,y) = 1.496e-6 * (x*sqrt(x)) /(x+120.)/y ! 1.496e-6*x**1.5/(x+120.)/y xka(x,y) = 1.414e3*viscos(x,y)*y diffac(a,b,c,d,e) = d*a*a/(xka(c,d)*rv*c*c)+1./(e*diffus(c,b)) venfac(a,b,c) = exp(log((viscos(b,c)/diffus(b,a)))*((.3333333))) & /sqrt(viscos(b,c))*sqrt(sqrt(den0/c)) conden(a,b,c,d,e) = (max(b,qmin)-c)/(1.+d*d/(rv*e)*c/(a*a)) ! idim = ite-its+1 kdim = kte-kts+1 ! kte_in = kte ! ! paddint 0 for negative values generated by dynamics ! do k = kts, kte do i = its, ite qci(i,k,1) = max(qci(i,k,1),0.0) qrs(i,k,1) = max(qrs(i,k,1),0.0) qci(i,k,2) = max(qci(i,k,2),0.0) qrs(i,k,2) = max(qrs(i,k,2),0.0) qrs(i,k,3) = max(qrs(i,k,3),0.0) ncr(i,k,1) = min(max(ncr(i,k,1),1.e8),2.e10) ncr(i,k,2) = max(ncr(i,k,2),0.0) ncr(i,k,3) = max(ncr(i,k,3),0.0) enddo enddo ! imsi do k = kts,kte do i = its,ite dend(i,k) = den(i,k) enddo enddo ! ! latent heat for phase changes and heat capacity. neglect the ! changes during microphysical process calculation ! emanuel(1994) ! do k = kts,kte_in do i = its,ite cpm(i,k) = cpmcal(q(i,k)) xl(i,k) = xlcal(t(i,k)) enddo enddo ! qcr(:,:) = 0.0 do i = its,ite if(slmsk(i).eq.2) then ! water qcr(i,:) = qc0 else qcr(i,:) = qc1 endif enddo ! do k = kts,kte do i = its,ite delz_tmp(i,k) = delz(i,k) den_tmp(i,k) = den(i,k) enddo enddo ! ! initialize the surface rain, snow, graupel ! do i = its,ite rainncv(i) = 0. if(present(snowncv) .and. present(snow)) snowncv(i) = 0. if(present (graupelncv) .and. present(graupel)) graupelncv(i) = 0. sr(i) = 0. ! ! new local array to catch step snow and graupel ! tstepsnow(i) = 0. tstepgraup(i) = 0. enddo ! ! compute the minor time steps. ! loops = max(nint(delt/dtcldcr),1) dtcld = delt/loops if(delt.le.dtcldcr) dtcld = delt ! do loop = 1,loops ! ! initialize the large scale variables ! do i = its,ite mstep(i) = 1 mnstep(i) = 1 flgcld(i) = .true. enddo ! do k = kts, kte call VREC( tvec1(its), den(its,k), ite-its+1) do i = its, ite tvec1(i) = tvec1(i)*den0 enddo call VSQRT( denfac(its,k), tvec1(its), ite-its+1) enddo ! ! Inline expansion for fpvs ! qs(i,k,1) = fpvs(t(i,k),0,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c) ! qs(i,k,2) = fpvs(t(i,k),1,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c) ! hsub = xls hvap = xlv0 cvap = cpv ttp = t0c+0.01 dldt = cvap-cliq xa = -dldt/rv xb = xa+hvap/(rv*ttp) dldti = cvap-cice xai = -dldti/rv xbi = xai+hsub/(rv*ttp) ! do k = kts,kte_in do i = its,ite tr = ttp/t(i,k) qs(i,k,1) = psat*exp(log(tr)*(xa))*exp(xb*(1.-tr)) qs(i,k,1) = min(qs(i,k,1),0.99*p(i,k)) qs(i,k,1) = ep2*qs(i,k,1)/(p(i,k)-qs(i,k,1)) qs(i,k,1) = max(qs(i,k,1),qmin) rh(i,k,1) = max(q(i,k)/qs(i,k,1),qmin) tr=ttp/t(i,k) if(t(i,k).lt.ttp) then qs(i,k,2) = psat*exp(log(tr)*(xai))*exp(xbi*(1.-tr)) else qs(i,k,2) = psat*exp(log(tr)*(xa))*exp(xb*(1.-tr)) endif qs(i,k,2) = min(qs(i,k,2),0.99*p(i,k)) qs(i,k,2) = ep2*qs(i,k,2)/(p(i,k)-qs(i,k,2)) qs(i,k,2) = max(qs(i,k,2),qmin) rh(i,k,2) = max(q(i,k)/qs(i,k,2),qmin) enddo enddo ! !---------------------------------------------------------------- ! initialize the variables for microphysical physics ! ! do k = kts, kte do i = its, ite prevp(i,k) = 0. psdep(i,k) = 0. pgdep(i,k) = 0. praut(i,k) = 0. psaut(i,k) = 0. pgaut(i,k) = 0. pracw(i,k) = 0. praci(i,k) = 0. piacr(i,k) = 0. psaci(i,k) = 0. psacw(i,k) = 0. pracs(i,k) = 0. psacr(i,k) = 0. pgacw(i,k) = 0. paacw(i,k) = 0. pgaci(i,k) = 0. pgacr(i,k) = 0. pgacs(i,k) = 0. pigen(i,k) = 0. pidep(i,k) = 0. pcond(i,k) = 0. psmlt(i,k) = 0. pgmlt(i,k) = 0. pseml(i,k) = 0. pgeml(i,k) = 0. psevp(i,k) = 0. pgevp(i,k) = 0. pcact(i,k) = 0. falk(i,k,1) = 0. falk(i,k,2) = 0. falk(i,k,3) = 0. fall(i,k,1) = 0. fall(i,k,2) = 0. fall(i,k,3) = 0. fallc(i,k) = 0. falkc(i,k) = 0. falln(i,k) =0. falkn(i,k) =0. xni(i,k) = 1.e3 nsacw(i,k) = 0. ngacw(i,k) = 0. naacw(i,k) = 0. niacr(i,k) = 0. nsacr(i,k) = 0. ngacr(i,k) = 0. nseml(i,k) = 0. ngeml(i,k) = 0. nracw(i,k) = 0. nccol(i,k) = 0. nrcol(i,k) = 0. ncact(i,k) = 0. nraut(i,k) = 0. ncevp(i,k) = 0. enddo enddo ! do k = kts, kte do i = its, ite if(qci(i,k,1).le.qmin .or. ncr(i,k,2).le.ncmin ) then rslopec(i,k) = rslopecmax rslopec2(i,k) = rslopec2max rslopec3(i,k) = rslopec3max else rslopec(i,k) = 1./lamdac(qci(i,k,1),den(i,k),ncr(i,k,2)) rslopec2(i,k) = rslopec(i,k)*rslopec(i,k) rslopec3(i,k) = rslopec2(i,k)*rslopec(i,k) endif !------------------------------------------------------------- ! Ni: ice crystal number concentraiton [HDC 5c] !------------------------------------------------------------- temp = (den(i,k)*max(qci(i,k,2),qmin)) temp = sqrt(sqrt(temp*temp*temp)) xni(i,k) = min(max(5.38e7*temp,1.e3),1.e6) enddo enddo !---------------------------------------------------------------- ! compute the fallout term: ! first, vertical terminal velosity for minor loops !---------------------------------------------------------------- do k = kts, kte do i = its, ite qrs_tmp(i,k,1) = qrs(i,k,1) qrs_tmp(i,k,2) = qrs(i,k,2) qrs_tmp(i,k,3) = qrs(i,k,3) ncr_tmp(i,k) = ncr(i,k,3) enddo enddo call slope_wdm6(qrs_tmp,ncr_tmp,den_tmp,denfac,t,rslope,rslopeb,rslope2, & rslope3,work1,workn,its,ite,kts,kte) ! ! vt update for qr and nr ! mstepmax = 1 numdt = 1 do k = kte, kts, -1 do i = its, ite work1(i,k,1) = work1(i,k,1)/delz(i,k) workn(i,k) = workn(i,k)/delz(i,k) numdt(i) = max(nint(max(work1(i,k,1),workn(i,k))*dtcld+.5),1) if(numdt(i).ge.mstep(i)) mstep(i) = numdt(i) enddo enddo do i = its, ite if(mstepmax.le.mstep(i)) mstepmax = mstep(i) enddo ! do n = 1, mstepmax k = kte do i = its, ite if(n.le.mstep(i)) then falk(i,k,1) = dend(i,k)*qrs(i,k,1)*work1(i,k,1)/mstep(i) falkn(i,k) = ncr(i,k,3)*workn(i,k)/mstep(i) fall(i,k,1) = fall(i,k,1)+falk(i,k,1) falln(i,k) = falln(i,k)+falkn(i,k) qrs(i,k,1) = max(qrs(i,k,1)-falk(i,k,1)*dtcld/dend(i,k),0.) ncr(i,k,3) = max(ncr(i,k,3)-falkn(i,k)*dtcld,0.) endif enddo ! do k = kte_in-1,kts,-1 do i = its,ite if(n.le.mstep(i)) then falk(i,k,1) = dend(i,k)*qrs(i,k,1)*work1(i,k,1)/mstep(i) falkn(i,k) = ncr(i,k,3)*workn(i,k)/mstep(i) fall(i,k,1) = fall(i,k,1)+falk(i,k,1) falln(i,k) = falln(i,k)+falkn(i,k) dqr(i,k) = min(falk(i,k,1)*dtcld/dend(i,k),qrs(i,k,1)) dqr(i,k+1) = min(falk(i,k+1,1)*delz(i,k+1)/delz(i,k) & *dtcld/dend(i,k),qrs(i,k+1,1)) dnr(i,k) = min(falkn(i,k)*dtcld,ncr(i,k,3)) dnr(i,k+1) = min(falkn(i,k+1)*delz(i,k+1)/delz(i,k)*dtcld, & ncr(i,k+1,3)) qrs(i,k,1) = max(qrs(i,k,1)-dqr(i,k)+dqr(i,k+1),0.) ncr(i,k,3) = max(ncr(i,k,3)-dnr(i,k)+dnr(i,k+1),0.) endif enddo enddo do k = kts, kte do i = its, ite qrs_tmp(i,k,1) = qrs(i,k,1) ncr_tmp(i,k) = ncr(i,k,3) enddo enddo call slope_rain(qrs_tmp,ncr_tmp,den_tmp,denfac,t,rslope,rslopeb,rslope2, & rslope3,work1,workn,its,ite,kts,kte) do k = kte, kts, -1 do i = its, ite work1(i,k,1) = work1(i,k,1)/delz(i,k) workn(i,k) = workn(i,k)/delz(i,k) enddo enddo enddo ! ! for semi ! do k = kte, kts, -1 do i = its, ite qsum(i,k) = max( (qrs(i,k,2)+qrs(i,k,3)), 1.E-15) if(qsum(i,k) .gt. 1.e-15 ) then worka(i,k) = (work1(i,k,2)*qrs(i,k,2) + work1(i,k,3)*qrs(i,k,3)) & /qsum(i,k) else worka(i,k) = 0. endif denqrs2(i,k) = den(i,k)*qrs(i,k,2) denqrs3(i,k) = den(i,k)*qrs(i,k,3) enddo enddo call nislfv_rain_plm6(idim,kdim,den_tmp,denfac,t,delz_tmp,worka, & denqrs2,denqrs3,delqrs2,delqrs3,dtcld,1,1) do k = kts, kte do i = its, ite qrs(i,k,2) = max(denqrs2(i,k)/den(i,k),0.) qrs(i,k,3) = max(denqrs3(i,k)/den(i,k),0.) fall(i,k,2) = denqrs2(i,k)*worka(i,k)/delz(i,k) fall(i,k,3) = denqrs3(i,k)*worka(i,k)/delz(i,k) enddo enddo do i = its, ite fall(i,1,2) = delqrs2(i)/delz(i,1)/dtcld fall(i,1,3) = delqrs3(i)/delz(i,1)/dtcld enddo do k = kts, kte do i = its, ite qrs_tmp(i,k,1) = qrs(i,k,1) qrs_tmp(i,k,2) = qrs(i,k,2) qrs_tmp(i,k,3) = qrs(i,k,3) ncr_tmp(i,k) = ncr(i,k,3) enddo enddo call slope_wdm6(qrs_tmp,ncr_tmp,den_tmp,denfac,t,rslope,rslopeb,rslope2, & rslope3,work1,workn,its,ite,kts,kte) ! do k = kte, kts, -1 do i = its, ite supcol = t0c-t(i,k) n0sfac(i,k) = max(min(exp(alpha*supcol),n0smax/n0s),1.) if(t(i,k).gt.t0c) then !--------------------------------------------------------------- ! psmlt: melting of snow [HL A33] [RH83 A25] ! (T>T0: QS->QR) !--------------------------------------------------------------- xlf = xlf0 work2(i,k) = venfac(p(i,k),t(i,k),den(i,k)) if(qrs(i,k,2).gt.0.) then coeres = rslope2(i,k,2)*sqrt(rslope(i,k,2)*rslopeb(i,k,2)) psmlt(i,k) = xka(t(i,k),den(i,k))/xlf*(t0c-t(i,k))*pi/2. & *n0sfac(i,k)*(precs1*rslope2(i,k,2) & +precs2*work2(i,k)*coeres)/den(i,k) psmlt(i,k) = min(max(psmlt(i,k)*dtcld/mstep(i),-qrs(i,k,2) & /mstep(i)),0.) qrs(i,k,2) = qrs(i,k,2) + psmlt(i,k) qrs(i,k,1) = qrs(i,k,1) - psmlt(i,k) !------------------------------------------------------------------- ! nsmlt: melting of snow [LH A27] ! (T>T0: ->NR) !------------------------------------------------------------------- if(qrs(i,k,2).gt.qcrmin) then sfac = rslope(i,k,2)*n0s*n0sfac(i,k)/qrs(i,k,2) ncr(i,k,3) = ncr(i,k,3) - sfac*psmlt(i,k) endif t(i,k) = t(i,k) + xlf/cpm(i,k)*psmlt(i,k) endif !--------------------------------------------------------------- ! pgmlt: melting of graupel [HL A23] [LFO 47] ! (T>T0: QG->QR) !--------------------------------------------------------------- if(qrs(i,k,3).gt.0.) then coeres = rslope2(i,k,3)*sqrt(rslope(i,k,3)*rslopeb(i,k,3)) pgmlt(i,k) = xka(t(i,k),den(i,k))/xlf*(t0c-t(i,k))*(precg1 & *rslope2(i,k,3) + precg2*work2(i,k)*coeres) & /den(i,k) pgmlt(i,k) = min(max(pgmlt(i,k)*dtcld/mstep(i), & -qrs(i,k,3)/mstep(i)),0.) qrs(i,k,3) = qrs(i,k,3) + pgmlt(i,k) qrs(i,k,1) = qrs(i,k,1) - pgmlt(i,k) !------------------------------------------------------------------- ! ngmlt: melting of graupel [LH A28] ! (T>T0: ->NR) !------------------------------------------------------------------- if(qrs(i,k,3).gt.qcrmin) then gfac = rslope(i,k,3)*n0g/qrs(i,k,3) ncr(i,k,3) = ncr(i,k,3) - gfac*pgmlt(i,k) endif t(i,k) = t(i,k) + xlf/cpm(i,k)*pgmlt(i,k) endif endif enddo enddo !--------------------------------------------------------------- ! Vice [ms-1] : fallout of ice crystal [HDC 5a] !--------------------------------------------------------------- do k = kte, kts, -1 do i = its, ite if(qci(i,k,2).le.0.) then work1c(i,k) = 0. else xmi = den(i,k)*qci(i,k,2)/xni(i,k) diameter = max(min(dicon * sqrt(xmi),dimax), 1.e-25) work1c(i,k) = 1.49e4*exp(log(diameter)*(1.31)) endif enddo enddo ! ! forward semi-laglangian scheme (JH), PCM (piecewise constant), (linear) ! do k = kte, kts, -1 do i = its, ite denqci(i,k) = den(i,k)*qci(i,k,2) enddo enddo call nislfv_rain_plmr(idim,kdim,den_tmp,denfac,t,delz_tmp,work1c,denqci,denqci, & delqi,dtcld,1,0,0) do k = kts, kte do i = its, ite qci(i,k,2) = max(denqci(i,k)/den(i,k),0.) enddo enddo do i = its, ite fallc(i,1) = delqi(i)/delz(i,1)/dtcld enddo ! !---------------------------------------------------------------- ! rain (unit is mm/sec;kgm-2s-1: /1000*delt ===> m)==> mm for wrf ! do i = its, ite fallsum = fall(i,kts,1)+fall(i,kts,2)+fall(i,kts,3)+fallc(i,kts) fallsum_qsi = fall(i,kts,2)+fallc(i,kts) fallsum_qg = fall(i,kts,3) if(fallsum.gt.0.) then rainncv(i) = fallsum*delz(i,kts)/denr*dtcld*1000. + rainncv(i) rain(i) = fallsum*delz(i,kts)/denr*dtcld*1000. + rain(i) endif If(fallsum_qsi.gt.0.) then tstepsnow(i) = fallsum_qsi*delz(i,kts)/denr*dtcld*1000. + tstepsnow(i) IF( PRESENT (snowncv) .AND. PRESENT (snow)) THEN snowncv(i) = fallsum_qsi*delz(i,kts)/denr*dtcld*1000. + snowncv(i) snow(i) = fallsum_qsi*delz(i,kts)/denr*dtcld*1000. + snow(i) ENDIF ENDIF IF(fallsum_qg.gt.0.) then tstepgraup(i) = fallsum_qg*delz(i,kts)/denr*dtcld*1000. & + tstepgraup(i) IF( PRESENT (graupelncv) .AND. PRESENT (graupel)) THEN graupelncv(i) = fallsum_qg*delz(i,kts)/denr*dtcld*1000. & + graupelncv(i) graupel(i) = fallsum_qg*delz(i,kts)/denr*dtcld*1000. + graupel(i) ENDIF ENDIF IF ( PRESENT (snowncv)) THEN if(fallsum.gt.0.)sr(i)=(snowncv(i) + graupelncv(i))/(rainncv(i)+1.e-12) ELSE if(fallsum.gt.0.)sr(i)=(tstepsnow(i) + tstepgraup(i))/(rainncv(i)+1.e-12) ENDIF enddo ! !--------------------------------------------------------------- ! pimlt: instantaneous melting of cloud ice [HL A47] [RH83 A28] ! (T>T0: QI->QC) !--------------------------------------------------------------- do k = kts, kte do i = its, ite supcol = t0c-t(i,k) xlf = xls-xl(i,k) if(supcol.lt.0.) xlf = xlf0 if(supcol.lt.0 .and. qci(i,k,2).gt.0.) then qci(i,k,1) = qci(i,k,1) + qci(i,k,2) !--------------------------------------------------------------- ! nimlt: instantaneous melting of cloud ice [LH A18] ! (T>T0: ->NC) !-------------------------------------------------------------- ncr(i,k,2) = ncr(i,k,2) + xni(i,k) t(i,k) = t(i,k) - xlf/cpm(i,k)*qci(i,k,2) qci(i,k,2) = 0. endif !--------------------------------------------------------------- ! pihmf: homogeneous of cloud water below -40c [HL A45] ! (T<-40C: QC->QI) !--------------------------------------------------------------- if(supcol.gt.40. .and. qci(i,k,1).gt.0.) then qci(i,k,2) = qci(i,k,2) + qci(i,k,1) !--------------------------------------------------------------- ! nihmf: homogeneous of cloud water below -40c [LH A17] ! (T<-40C: NC->) !--------------------------------------------------------------- if(ncr(i,k,2).gt.0.) ncr(i,k,2) = 0. t(i,k) = t(i,k) + xlf/cpm(i,k)*qci(i,k,1) qci(i,k,1) = 0. endif !--------------------------------------------------------------- ! pihtf: heterogeneous of cloud water [HL A44] ! (T0>T>-40C: QC->QI) !--------------------------------------------------------------- if(supcol.gt.0. .and. qci(i,k,1).gt.qmin) then supcolt=min(supcol,70.) pfrzdtc = min(pi*pi*pfrz1*(exp(pfrz2*supcolt)-1.)*denr/den(i,k) & *ncr(i,k,2)*rslopec3(i,k)*rslopec3(i,k)/18.*dtcld & ,qci(i,k,1)) !--------------------------------------------------------------- ! nihtf: heterogeneous of cloud water [LH A16] ! (T0>T>-40C: NC->) !--------------------------------------------------------------- if(ncr(i,k,2).gt.ncmin) then nfrzdtc = min(pi*pfrz1*(exp(pfrz2*supcolt)-1.)*ncr(i,k,2) & *rslopec3(i,k)/6.*dtcld,ncr(i,k,2)) ncr(i,k,2) = ncr(i,k,2) - nfrzdtc endif qci(i,k,2) = qci(i,k,2) + pfrzdtc t(i,k) = t(i,k) + xlf/cpm(i,k)*pfrzdtc qci(i,k,1) = qci(i,k,1)-pfrzdtc endif !--------------------------------------------------------------- ! pgfrz: freezing of rain water [HL A20] [LFO 45] ! (TQG) !--------------------------------------------------------------- if(supcol.gt.0. .and. qrs(i,k,1).gt.0.) then supcolt=min(supcol,70.) pfrzdtr = min(140.*(pi*pi)*pfrz1*ncr(i,k,3)*denr/den(i,k) & *(exp(pfrz2*supcolt)-1.)*rslope3(i,k,1)*rslope3(i,k,1) & *dtcld,qrs(i,k,1)) !--------------------------------------------------------------- ! ngfrz: freezing of rain water [LH A26] ! (T ) !--------------------------------------------------------------- if(ncr(i,k,3).gt.nrmin) then nfrzdtr = min(4.*pi*pfrz1*ncr(i,k,3)*(exp(pfrz2*supcolt)-1.) & *rslope3(i,k,1)*dtcld, ncr(i,k,3)) ncr(i,k,3) = ncr(i,k,3) - nfrzdtr endif qrs(i,k,3) = qrs(i,k,3) + pfrzdtr t(i,k) = t(i,k) + xlf/cpm(i,k)*pfrzdtr qrs(i,k,1) = qrs(i,k,1) - pfrzdtr endif enddo enddo ! do k = kts, kte do i = its, ite ncr(i,k,2) = max(ncr(i,k,2),0.0) ncr(i,k,3) = max(ncr(i,k,3),0.0) enddo enddo ! !---------------------------------------------------------------- ! update the slope parameters for microphysics computation ! do k = kts, kte do i = its, ite qrs_tmp(i,k,1) = qrs(i,k,1) qrs_tmp(i,k,2) = qrs(i,k,2) qrs_tmp(i,k,3) = qrs(i,k,3) ncr_tmp(i,k) = ncr(i,k,3) enddo enddo call slope_wdm6(qrs_tmp,ncr_tmp,den_tmp,denfac,t,rslope,rslopeb,rslope2, & rslope3,work1,workn,its,ite,kts,kte) do k = kts, kte do i = its, ite !----------------------------------------------------------------- ! compute the mean-volume drop diameter [LH A10] ! for raindrop distribution !----------------------------------------------------------------- avedia(i,k,2) = rslope(i,k,1)*((24.)**(.3333333)) ! if(qci(i,k,1).le.qmin .or. ncr(i,k,2).le.ncmin) then rslopec(i,k) = rslopecmax rslopec2(i,k) = rslopec2max rslopec3(i,k) = rslopec3max else rslopec(i,k) = 1./lamdac(qci(i,k,1),den(i,k),ncr(i,k,2)) rslopec2(i,k) = rslopec(i,k)*rslopec(i,k) rslopec3(i,k) = rslopec2(i,k)*rslopec(i,k) endif !-------------------------------------------------------------------- ! compute the mean-volume drop diameter [LH A7] ! for cloud-droplet distribution !-------------------------------------------------------------------- avedia(i,k,1) = rslopec(i,k) enddo enddo ! do k = kts, kte do i = its, ite work1(i,k,1) = diffac(xl(i,k),p(i,k),t(i,k),den(i,k),qs(i,k,1)) work1(i,k,2) = diffac(xls,p(i,k),t(i,k),den(i,k),qs(i,k,2)) work2(i,k) = venfac(p(i,k),t(i,k),den(i,k)) enddo enddo ! !=============================================================== ! ! warm rain processes ! ! - follows the double-moment processes in Lim and Hong ! !=============================================================== ! do k = kts, kte do i = its, ite supsat = max(q(i,k),qmin)-qs(i,k,1) satdt = supsat/dtcld !--------------------------------------------------------------- ! praut: auto conversion rate from cloud to rain [LH 9] [CP 17] ! (QC->QR) !-------------------------------------------------------------- lencon = 2.7e-2*den(i,k)*qci(i,k,1)*(1.e20/16.*rslopec2(i,k) & *rslopec2(i,k)-0.4) lenconcr = max(1.2*lencon, qcrmin) if(qci(i,k,1).gt.qcr(i,k)) then praut(i,k) = qck1*qci(i,k,1)**(7./3.)*ncr(i,k,2)**(-1./3.) praut(i,k) = min(praut(i,k),qci(i,k,1)/dtcld) !---------------------------------------------------------------------- ! nraut: auto conversion rate from cloud to rain [LH A6] [CP 18 & 19] ! (NC->NR) !------------------------------------------------------------------------ nraut(i,k) = 3.5e9*den(i,k)*praut(i,k) if(qrs(i,k,1).gt.lenconcr) & nraut(i,k) = ncr(i,k,3)/qrs(i,k,1)*praut(i,k) nraut(i,k) = min(nraut(i,k),ncr(i,k,2)/dtcld) endif !--------------------------------------------------------------- ! pracw: accretion of cloud water by rain [LH 10] [CP 22 & 23] ! (QC->QR) ! nracw: accretion of cloud water by rain [LH A9] ! (NC->) !--------------------------------------------------------------- if(qrs(i,k,1).ge.lenconcr) then if(avedia(i,k,2).ge.di100) then nracw(i,k) = min(ncrk1*ncr(i,k,2)*ncr(i,k,3)*(rslopec3(i,k) & + 24.*rslope3(i,k,1)),ncr(i,k,2)/dtcld) pracw(i,k) = min(pi/6.*(denr/den(i,k))*ncrk1*ncr(i,k,2) & *ncr(i,k,3)*rslopec3(i,k)*(2.*rslopec3(i,k) & + 24.*rslope3(i,k,1)),qci(i,k,1)/dtcld) else nracw(i,k) = min(ncrk2*ncr(i,k,2)*ncr(i,k,3)*(2.*rslopec3(i,k) & *rslopec3(i,k)+5040.*rslope3(i,k,1) & *rslope3(i,k,1)),ncr(i,k,2)/dtcld) pracw(i,k) = min(pi/6.*(denr/den(i,k))*ncrk2*ncr(i,k,2) & *ncr(i,k,3)*rslopec3(i,k)*(6.*rslopec3(i,k) & *rslopec3(i,k)+5040.*rslope3(i,k,1)*rslope3(i,k,1)) & ,qci(i,k,1)/dtcld) endif endif !---------------------------------------------------------------- ! nccol: self collection of cloud water [LH A8] [CP 24 & 25] ! (NC->) !---------------------------------------------------------------- if(avedia(i,k,1).ge.di100) then nccol(i,k) = ncrk1*ncr(i,k,2)*ncr(i,k,2)*rslopec3(i,k) else nccol(i,k) = 2.*ncrk2*ncr(i,k,2)*ncr(i,k,2)*rslopec3(i,k) & *rslopec3(i,k) endif !---------------------------------------------------------------- ! nrcol: self collection of rain-drops and break-up [LH A21] [CP 24 & 25] ! (NR->) !---------------------------------------------------------------- if(qrs(i,k,1).ge.lenconcr) then if(avedia(i,k,2).lt.di100) then nrcol(i,k) = 5040.*ncrk2*ncr(i,k,3)*ncr(i,k,3)*rslope3(i,k,1) & *rslope3(i,k,1) elseif(avedia(i,k,2).ge.di100 .and. avedia(i,k,2).lt.di600) then nrcol(i,k) = 24.*ncrk1*ncr(i,k,3)*ncr(i,k,3)*rslope3(i,k,1) elseif(avedia(i,k,2).ge.di600 .and. avedia(i,k,2).lt.di2000) then coecol = -2.5e3*(avedia(i,k,2)-di600) nrcol(i,k) = 24.*exp(coecol)*ncrk1*ncr(i,k,3)*ncr(i,k,3) & *rslope3(i,k,1) else nrcol(i,k) = 0. endif endif !--------------------------------------------------------------- ! prevp: evaporation/condensation rate of rain [HL A41] ! (QV->QR or QR->QV) !--------------------------------------------------------------- if(qrs(i,k,1).gt.0.) then coeres = rslope(i,k,1)*sqrt(rslope(i,k,1)*rslopeb(i,k,1)) prevp(i,k) = (rh(i,k,1)-1.)*ncr(i,k,3)*(precr1*rslope(i,k,1) & + precr2*work2(i,k)*coeres)/work1(i,k,1) if(prevp(i,k).lt.0.) then prevp(i,k) = max(prevp(i,k),-qrs(i,k,1)/dtcld) prevp(i,k) = max(prevp(i,k),satdt/2) !---------------------------------------------------------------- ! Nrevp: evaporation/condensation rate of rain [LH A14] ! (NR->NCCN) !---------------------------------------------------------------- if(prevp(i,k).eq.-qrs(i,k,1)/dtcld) then ncr(i,k,1) = ncr(i,k,1)+ncr(i,k,3) ncr(i,k,3) = 0. endif else ! prevp(i,k) = min(prevp(i,k),satdt/2) endif endif enddo enddo ! !=============================================================== ! ! cold rain processes ! ! - follows the revised ice microphysics processes in HDC ! - the processes same as in RH83 and RH84 and LFO behave ! following ice crystal hapits defined in HDC, inclduing ! intercept parameter for snow (n0s), ice crystal number ! concentration (ni), ice nuclei number concentration ! (n0i), ice diameter (d) ! !=============================================================== ! do k = kts, kte do i = its, ite supcol = t0c-t(i,k) n0sfac(i,k) = max(min(exp(alpha*supcol),n0smax/n0s),1.) supsat = max(q(i,k),qmin)-qs(i,k,2) satdt = supsat/dtcld ifsat = 0 !------------------------------------------------------------- ! Ni: ice crystal number concentraiton [HDC 5c] !------------------------------------------------------------- ! xni(i,k) = min(max(5.38e7*(den(i,k) & ! *max(qci(i,k,2),qmin))**0.75,1.e3),1.e6) temp = (den(i,k)*max(qci(i,k,2),qmin)) temp = sqrt(sqrt(temp*temp*temp)) xni(i,k) = min(max(5.38e7*temp,1.e3),1.e6) eacrs = exp(0.07*(-supcol)) ! xmi = den(i,k)*qci(i,k,2)/xni(i,k) diameter = min(dicon * sqrt(xmi),dimax) vt2i = 1.49e4*diameter**1.31 vt2r=pvtr*rslopeb(i,k,1)*denfac(i,k) vt2s=pvts*rslopeb(i,k,2)*denfac(i,k) vt2g=pvtg*rslopeb(i,k,3)*denfac(i,k) qsum(i,k) = max((qrs(i,k,2)+qrs(i,k,3)),1.e-15) if(qsum(i,k) .gt. 1.e-15) then vt2ave=(vt2s*qrs(i,k,2)+vt2g*qrs(i,k,3))/(qsum(i,k)) else vt2ave=0. endif if(supcol.gt.0. .and. qci(i,k,2).gt.qmin) then if(qrs(i,k,1).gt.qcrmin) then !------------------------------------------------------------- ! praci: Accretion of cloud ice by rain [HL A15] [LFO 25] ! (TQR) !------------------------------------------------------------- acrfac = 6.*rslope2(i,k,1)+4.*diameter*rslope(i,k,1) + diameter**2 praci(i,k) = pi*qci(i,k,2)*ncr(i,k,3)*abs(vt2r-vt2i)*acrfac/4. praci(i,k) = min(praci(i,k),qci(i,k,2)/dtcld) !------------------------------------------------------------- ! piacr: Accretion of rain by cloud ice [HL A19] [LFO 26] ! (TQS or QR->QG) !------------------------------------------------------------- piacr(i,k) = pi*pi*avtr*ncr(i,k,3)*denr*xni(i,k)*denfac(i,k) & *g7pbr*rslope3(i,k,1)*rslope2(i,k,1)*rslopeb(i,k,1) & /24./den(i,k) piacr(i,k) = min(piacr(i,k),qrs(i,k,1)/dtcld) endif !------------------------------------------------------------- ! niacr: Accretion of rain by cloud ice [LH A25] ! (T) !------------------------------------------------------------- if(ncr(i,k,3).gt.nrmin) then niacr(i,k) = pi*avtr*ncr(i,k,3)*xni(i,k)*denfac(i,k)*g4pbr & *rslope2(i,k,1)*rslopeb(i,k,1)/4. niacr(i,k) = min(niacr(i,k),ncr(i,k,3)/dtcld) endif !------------------------------------------------------------- ! psaci: Accretion of cloud ice by snow [HDC 10] ! (TQS) !------------------------------------------------------------- if(qrs(i,k,2).gt.qcrmin) then acrfac = 2.*rslope3(i,k,2)+2.*diameter*rslope2(i,k,2) & + diameter**2*rslope(i,k,2) psaci(i,k) = pi*qci(i,k,2)*eacrs*n0s*n0sfac(i,k) & *abs(vt2ave-vt2i)*acrfac/4. psaci(i,k) = min(psaci(i,k),qci(i,k,2)/dtcld) endif !------------------------------------------------------------- ! pgaci: Accretion of cloud ice by graupel [HL A17] [LFO 41] ! (TQG) !------------------------------------------------------------- if(qrs(i,k,3).gt.qcrmin) then egi = exp(0.07*(-supcol)) acrfac = 2.*rslope3(i,k,3)+2.*diameter*rslope2(i,k,3) & + diameter**2*rslope(i,k,3) pgaci(i,k) = pi*egi*qci(i,k,2)*n0g*abs(vt2ave-vt2i)*acrfac/4. pgaci(i,k) = min(pgaci(i,k),qci(i,k,2)/dtcld) endif endif !------------------------------------------------------------- ! psacw: Accretion of cloud water by snow [HL A7] [LFO 24] ! (TQS, and T>=T0: QC->QR) !------------------------------------------------------------- if(qrs(i,k,2).gt.qcrmin .and. qci(i,k,1).gt.qmin) then psacw(i,k) = min(pacrc*n0sfac(i,k)*rslope3(i,k,2)*rslopeb(i,k,2) & *qci(i,k,1)*denfac(i,k),qci(i,k,1)/dtcld) endif !------------------------------------------------------------- ! nsacw: Accretion of cloud water by snow [LH A12] ! (NC ->) !------------------------------------------------------------- if(qrs(i,k,2).gt.qcrmin .and. ncr(i,k,2).gt.ncmin) then nsacw(i,k) = min(pacrc*n0sfac(i,k)*rslope3(i,k,2)*rslopeb(i,k,2) & *ncr(i,k,2)*denfac(i,k),ncr(i,k,2)/dtcld) endif !------------------------------------------------------------- ! pgacw: Accretion of cloud water by graupel [HL A6] [LFO 40] ! (TQG, and T>=T0: QC->QR) !------------------------------------------------------------- if(qrs(i,k,3).gt.qcrmin .and. qci(i,k,1).gt.qmin) then pgacw(i,k) = min(pacrg*rslope3(i,k,3)*rslopeb(i,k,3)*qci(i,k,1) & *denfac(i,k),qci(i,k,1)/dtcld) endif !------------------------------------------------------------- ! ngacw: Accretion of cloud water by graupel [LH A13] ! (NC-> !------------------------------------------------------------- if(qrs(i,k,3).gt.qcrmin .and. ncr(i,k,2).gt.ncmin) then ngacw(i,k) = min(pacrg*rslope3(i,k,3)*rslopeb(i,k,3)*ncr(i,k,2) & *denfac(i,k),ncr(i,k,2)/dtcld) endif !------------------------------------------------------------- ! paacw: Accretion of cloud water by averaged snow/graupel ! (TQG or QS, and T>=T0: QC->QR) !------------------------------------------------------------- if(qsum(i,k) .gt. 1.e-15 ) then paacw(i,k) = (qrs(i,k,2)*psacw(i,k)+qrs(i,k,3)*pgacw(i,k))/(qsum(i,k)) !------------------------------------------------------------- ! naacw: Accretion of cloud water by averaged snow/graupel ! (Nc->) !------------------------------------------------------------- naacw(i,k) = (qrs(i,k,2)*nsacw(i,k)+qrs(i,k,3)*ngacw(i,k))/(qsum(i,k)) endif !------------------------------------------------------------- ! pracs: Accretion of snow by rain [HL A11] [LFO 27] ! (TQG) !------------------------------------------------------------- if(qrs(i,k,2).gt.qcrmin .and. qrs(i,k,1).gt.qcrmin) then if(supcol.gt.0) then acrfac = 5.*rslope3(i,k,2)*rslope3(i,k,2) & + 4.*rslope3(i,k,2)*rslope2(i,k,2)*rslope(i,k,1) & + 1.5*rslope2(i,k,2)*rslope2(i,k,2)*rslope2(i,k,1) pracs(i,k) = pi*pi*ncr(i,k,3)*n0s*n0sfac(i,k)*abs(vt2r-vt2ave) & *(dens/den(i,k))*acrfac pracs(i,k) = min(pracs(i,k),qrs(i,k,2)/dtcld) endif !------------------------------------------------------------- ! psacr: Accretion of rain by snow [HL A10] [LFO 28] ! (TQS or QR->QG) (T>=T0: enhance melting of snow) !------------------------------------------------------------- acrfac = 30.*rslope3(i,k,1)*rslope2(i,k,1)*rslope(i,k,2) & +10.*rslope2(i,k,1)*rslope2(i,k,1)*rslope2(i,k,2) & + 2.*rslope3(i,k,1)*rslope3(i,k,2) psacr(i,k) = pi*pi*ncr(i,k,3)*n0s*n0sfac(i,k)*abs(vt2ave-vt2r) & *(denr/den(i,k))*acrfac psacr(i,k) = min(psacr(i,k),qrs(i,k,1)/dtcld) endif if(qrs(i,k,2).gt.qcrmin .and. ncr(i,k,3).gt.nrmin) then !------------------------------------------------------------- ! nsacr: Accretion of rain by snow [LH A23] ! (T) !------------------------------------------------------------- acrfac = 1.5*rslope2(i,k,1)*rslope(i,k,2) & + 1.0*rslope(i,k,1)*rslope2(i,k,2)+.5*rslope3(i,k,2) nsacr(i,k) = pi*ncr(i,k,3)*n0s*n0sfac(i,k)*abs(vt2ave-vt2r) & *acrfac nsacr(i,k) = min(nsacr(i,k),ncr(i,k,3)/dtcld) endif !------------------------------------------------------------- ! pgacr: Accretion of rain by graupel [HL A12] [LFO 42] ! (TQG) (T>=T0: enhance melting of graupel) !------------------------------------------------------------- if(qrs(i,k,3).gt.qcrmin .and. qrs(i,k,1).gt.qcrmin) then acrfac = 30.*rslope3(i,k,1)*rslope2(i,k,1)*rslope(i,k,3) & +10.*rslope2(i,k,1)*rslope2(i,k,1)*rslope2(i,k,3) & + 2.*rslope3(i,k,1)*rslope3(i,k,3) pgacr(i,k) = pi*pi*ncr(i,k,3)*n0g*abs(vt2ave-vt2r)*(denr/den(i,k)) & *acrfac pgacr(i,k) = min(pgacr(i,k),qrs(i,k,1)/dtcld) endif !------------------------------------------------------------- ! ngacr: Accretion of rain by graupel [LH A24] ! (T) !------------------------------------------------------------- if(qrs(i,k,3).gt.qcrmin .and. ncr(i,k,3).gt.nrmin) then acrfac = 1.5*rslope2(i,k,1)*rslope(i,k,3) & + 1.0*rslope(i,k,1)*rslope2(i,k,3) + .5*rslope3(i,k,3) ngacr(i,k) = pi*ncr(i,k,3)*n0g*abs(vt2ave-vt2r)*acrfac ngacr(i,k) = min(ngacr(i,k),ncr(i,k,3)/dtcld) endif ! !------------------------------------------------------------- ! pgacs: Accretion of snow by graupel [HL A13] [LFO 29] ! (QS->QG) : This process is eliminated in V3.0 with the ! new combined snow/graupel fall speeds !------------------------------------------------------------- if(qrs(i,k,3).gt.qcrmin .and. qrs(i,k,2).gt.qcrmin) then pgacs(i,k) = 0. endif if(supcol.le.0) then xlf = xlf0 !------------------------------------------------------------- ! pseml: Enhanced melting of snow by accretion of water [HL A34] ! (T>=T0: QS->QR) !------------------------------------------------------------- if(qrs(i,k,2).gt.0.) & pseml(i,k) = min(max(cliq*supcol*(paacw(i,k)+psacr(i,k)) & /xlf,-qrs(i,k,2)/dtcld),0.) !-------------------------------------------------------------- ! nseml: Enhanced melting of snow by accretion of water [LH A29] ! (T>=T0: ->NR) !-------------------------------------------------------------- if (qrs(i,k,2).gt.qcrmin) then sfac = rslope(i,k,2)*n0s*n0sfac(i,k)/qrs(i,k,2) nseml(i,k) = -sfac*pseml(i,k) endif !------------------------------------------------------------- ! pgeml: Enhanced melting of graupel by accretion of water [HL A24] [RH84 A21-A22] ! (T>=T0: QG->QR) !------------------------------------------------------------- if(qrs(i,k,3).gt.0.) & pgeml(i,k) = min(max(cliq*supcol*(paacw(i,k)+pgacr(i,k))/xlf & ,-qrs(i,k,3)/dtcld),0.) !-------------------------------------------------------------- ! ngeml: Enhanced melting of graupel by accretion of water [LH A30] ! (T>=T0: -> NR) !-------------------------------------------------------------- if (qrs(i,k,3).gt.qcrmin) then gfac = rslope(i,k,3)*n0g/qrs(i,k,3) ngeml(i,k) = -gfac*pgeml(i,k) endif endif if(supcol.gt.0) then !------------------------------------------------------------- ! pidep: Deposition/Sublimation rate of ice [HDC 9] ! (TQI or QI->QV) !------------------------------------------------------------- if(qci(i,k,2).gt.0. .and. ifsat.ne.1) then pidep(i,k) = 4.*diameter*xni(i,k)*(rh(i,k,2)-1.)/work1(i,k,2) supice = satdt-prevp(i,k) if(pidep(i,k).lt.0.) then pidep(i,k) = max(max(pidep(i,k),satdt/2),supice) pidep(i,k) = max(pidep(i,k),-qci(i,k,2)/dtcld) else pidep(i,k) = min(min(pidep(i,k),satdt/2),supice) endif if(abs(prevp(i,k)+pidep(i,k)).ge.abs(satdt)) ifsat = 1 endif !------------------------------------------------------------- ! psdep: deposition/sublimation rate of snow [HDC 14] ! (TQS or QS->QV) !------------------------------------------------------------- if(qrs(i,k,2).gt.0. .and. ifsat.ne.1) then coeres = rslope2(i,k,2)*sqrt(rslope(i,k,2)*rslopeb(i,k,2)) psdep(i,k) = (rh(i,k,2)-1.)*n0sfac(i,k)*(precs1*rslope2(i,k,2) & + precs2*work2(i,k)*coeres)/work1(i,k,2) supice = satdt-prevp(i,k)-pidep(i,k) if(psdep(i,k).lt.0.) then psdep(i,k) = max(psdep(i,k),-qrs(i,k,2)/dtcld) psdep(i,k) = max(max(psdep(i,k),satdt/2),supice) else psdep(i,k) = min(min(psdep(i,k),satdt/2),supice) endif if(abs(prevp(i,k)+pidep(i,k)+psdep(i,k)).ge.abs(satdt)) ifsat = 1 endif !------------------------------------------------------------- ! pgdep: deposition/sublimation rate of graupel [HL A21] [LFO 46] ! (TQG or QG->QV) !------------------------------------------------------------- if(qrs(i,k,3).gt.0. .and. ifsat.ne.1) then coeres = rslope2(i,k,3)*sqrt(rslope(i,k,3)*rslopeb(i,k,3)) pgdep(i,k) = (rh(i,k,2)-1.)*(precg1*rslope2(i,k,3) & + precg2*work2(i,k)*coeres)/work1(i,k,2) supice = satdt-prevp(i,k)-pidep(i,k)-psdep(i,k) if(pgdep(i,k).lt.0.) then pgdep(i,k) = max(pgdep(i,k),-qrs(i,k,3)/dtcld) pgdep(i,k) = max(max(pgdep(i,k),satdt/2),supice) else pgdep(i,k) = min(min(pgdep(i,k),satdt/2),supice) endif if(abs(prevp(i,k)+pidep(i,k)+psdep(i,k)+pgdep(i,k)).ge. & abs(satdt)) ifsat = 1 endif !------------------------------------------------------------- ! pigen: generation(nucleation) of ice from vapor [HL 50] [HDC 7-8] ! (TQI) !------------------------------------------------------------- if(supsat.gt.0. .and. ifsat.ne.1) then supice = satdt-prevp(i,k)-pidep(i,k)-psdep(i,k)-pgdep(i,k) xni0 = 1.e3*exp(0.1*supcol) roqi0 = 4.92e-11*xni0**1.33 pigen(i,k) = max(0.,(roqi0/den(i,k)-max(qci(i,k,2),0.))/dtcld) pigen(i,k) = min(min(pigen(i,k),satdt),supice) endif ! !------------------------------------------------------------- ! psaut: conversion(aggregation) of ice to snow [HDC 12] ! (TQS) !------------------------------------------------------------- if(qci(i,k,2).gt.0.) then qimax = roqimax/den(i,k) psaut(i,k) = max(0.,(qci(i,k,2)-qimax)/dtcld) endif ! !------------------------------------------------------------- ! pgaut: conversion(aggregation) of snow to graupel [HL A4] [LFO 37] ! (TQG) !------------------------------------------------------------- if(qrs(i,k,2).gt.0.) then alpha2 = 1.e-3*exp(0.09*(-supcol)) pgaut(i,k) = min(max(0.,alpha2*(qrs(i,k,2)-qs0)),qrs(i,k,2)/dtcld) endif endif ! !------------------------------------------------------------- ! psevp: Evaporation of melting snow [HL A35] [RH83 A27] ! (T>=T0: QS->QV) !------------------------------------------------------------- if(supcol.lt.0.) then if(qrs(i,k,2).gt.0. .and. rh(i,k,1).lt.1.) then coeres = rslope2(i,k,2)*sqrt(rslope(i,k,2)*rslopeb(i,k,2)) psevp(i,k) = (rh(i,k,1)-1.)*n0sfac(i,k)*(precs1*rslope2(i,k,2) & +precs2*work2(i,k)*coeres)/work1(i,k,1) psevp(i,k) = min(max(psevp(i,k),-qrs(i,k,2)/dtcld),0.) endif !------------------------------------------------------------- ! pgevp: Evaporation of melting graupel [HL A25] [RH84 A19] ! (T>=T0: QG->QV) !------------------------------------------------------------- if(qrs(i,k,3).gt.0. .and. rh(i,k,1).lt.1.) then coeres = rslope2(i,k,3)*sqrt(rslope(i,k,3)*rslopeb(i,k,3)) pgevp(i,k) = (rh(i,k,1)-1.)*(precg1*rslope2(i,k,3) & + precg2*work2(i,k)*coeres)/work1(i,k,1) pgevp(i,k) = min(max(pgevp(i,k),-qrs(i,k,3)/dtcld),0.) endif endif enddo enddo ! ! !---------------------------------------------------------------- ! check mass conservation of generation terms and feedback to the ! large scale ! do k = kts, kte do i = its, ite ! delta2=0. delta3=0. if(qrs(i,k,1).lt.1.e-4 .and. qrs(i,k,2).lt.1.e-4) delta2=1. if(qrs(i,k,1).lt.1.e-4) delta3=1. if(t(i,k).le.t0c) then ! ! cloud water ! value = max(qmin,qci(i,k,1)) source = (praut(i,k)+pracw(i,k)+paacw(i,k)+paacw(i,k)) & *dtcld if (source.gt.value) then factor = value/source praut(i,k) = praut(i,k)*factor pracw(i,k) = pracw(i,k)*factor paacw(i,k) = paacw(i,k)*factor endif ! ! cloud ice ! value = max(qmin,qci(i,k,2)) source = (psaut(i,k)-pigen(i,k)-pidep(i,k)+praci(i,k)+psaci(i,k) & +pgaci(i,k))*dtcld if (source.gt.value) then factor = value/source psaut(i,k) = psaut(i,k)*factor pigen(i,k) = pigen(i,k)*factor pidep(i,k) = pidep(i,k)*factor praci(i,k) = praci(i,k)*factor psaci(i,k) = psaci(i,k)*factor pgaci(i,k) = pgaci(i,k)*factor endif ! ! rain ! value = max(qmin,qrs(i,k,1)) source = (-praut(i,k)-prevp(i,k)-pracw(i,k)+piacr(i,k) & +psacr(i,k)+pgacr(i,k))*dtcld if (source.gt.value) then factor = value/source praut(i,k) = praut(i,k)*factor prevp(i,k) = prevp(i,k)*factor pracw(i,k) = pracw(i,k)*factor piacr(i,k) = piacr(i,k)*factor psacr(i,k) = psacr(i,k)*factor pgacr(i,k) = pgacr(i,k)*factor endif ! ! snow ! value = max(qmin,qrs(i,k,2)) source = -(psdep(i,k)+psaut(i,k)-pgaut(i,k)+paacw(i,k) & +piacr(i,k)*delta3+praci(i,k)*delta3 & -pracs(i,k)*(1.-delta2)+psacr(i,k)*delta2 & +psaci(i,k)-pgacs(i,k) )*dtcld if (source.gt.value) then factor = value/source psdep(i,k) = psdep(i,k)*factor psaut(i,k) = psaut(i,k)*factor pgaut(i,k) = pgaut(i,k)*factor paacw(i,k) = paacw(i,k)*factor piacr(i,k) = piacr(i,k)*factor praci(i,k) = praci(i,k)*factor psaci(i,k) = psaci(i,k)*factor pracs(i,k) = pracs(i,k)*factor psacr(i,k) = psacr(i,k)*factor pgacs(i,k) = pgacs(i,k)*factor endif ! ! graupel ! value = max(qmin,qrs(i,k,3)) source = -(pgdep(i,k)+pgaut(i,k) & +piacr(i,k)*(1.-delta3)+praci(i,k)*(1.-delta3) & +psacr(i,k)*(1.-delta2)+pracs(i,k)*(1.-delta2) & +pgaci(i,k)+paacw(i,k)+pgacr(i,k)+pgacs(i,k))*dtcld if (source.gt.value) then factor = value/source pgdep(i,k) = pgdep(i,k)*factor pgaut(i,k) = pgaut(i,k)*factor piacr(i,k) = piacr(i,k)*factor praci(i,k) = praci(i,k)*factor psacr(i,k) = psacr(i,k)*factor pracs(i,k) = pracs(i,k)*factor paacw(i,k) = paacw(i,k)*factor pgaci(i,k) = pgaci(i,k)*factor pgacr(i,k) = pgacr(i,k)*factor pgacs(i,k) = pgacs(i,k)*factor endif ! ! cloud ! value = max(ncmin,ncr(i,k,2)) source = (nraut(i,k)+nccol(i,k)+nracw(i,k) & +naacw(i,k)+naacw(i,k))*dtcld if (source.gt.value) then factor = value/source nraut(i,k) = nraut(i,k)*factor nccol(i,k) = nccol(i,k)*factor nracw(i,k) = nracw(i,k)*factor naacw(i,k) = naacw(i,k)*factor endif ! ! rain ! value = max(nrmin,ncr(i,k,3)) source = (-nraut(i,k)+nrcol(i,k)+niacr(i,k)+nsacr(i,k)+ngacr(i,k) & )*dtcld if (source.gt.value) then factor = value/source nraut(i,k) = nraut(i,k)*factor nrcol(i,k) = nrcol(i,k)*factor niacr(i,k) = niacr(i,k)*factor nsacr(i,k) = nsacr(i,k)*factor ngacr(i,k) = ngacr(i,k)*factor endif ! work2(i,k)=-(prevp(i,k)+psdep(i,k)+pgdep(i,k)+pigen(i,k)+pidep(i,k)) ! update q(i,k) = q(i,k)+work2(i,k)*dtcld qci(i,k,1) = max(qci(i,k,1)-(praut(i,k)+pracw(i,k) & +paacw(i,k)+paacw(i,k))*dtcld,0.) qrs(i,k,1) = max(qrs(i,k,1)+(praut(i,k)+pracw(i,k) & +prevp(i,k)-piacr(i,k)-pgacr(i,k) & -psacr(i,k))*dtcld,0.) qci(i,k,2) = max(qci(i,k,2)-(psaut(i,k)+praci(i,k) & +psaci(i,k)+pgaci(i,k)-pigen(i,k)-pidep(i,k)) & *dtcld,0.) qrs(i,k,2) = max(qrs(i,k,2)+(psdep(i,k)+psaut(i,k)+paacw(i,k) & -pgaut(i,k)+piacr(i,k)*delta3 & +praci(i,k)*delta3+psaci(i,k)-pgacs(i,k) & -pracs(i,k)*(1.-delta2)+psacr(i,k)*delta2) & *dtcld,0.) qrs(i,k,3) = max(qrs(i,k,3)+(pgdep(i,k)+pgaut(i,k) & +piacr(i,k)*(1.-delta3) & +praci(i,k)*(1.-delta3)+psacr(i,k)*(1.-delta2) & +pracs(i,k)*(1.-delta2)+pgaci(i,k)+paacw(i,k) & +pgacr(i,k)+pgacs(i,k))*dtcld,0.) ncr(i,k,2) = max(ncr(i,k,2)+(-nraut(i,k)-nccol(i,k)-nracw(i,k) & -naacw(i,k)-naacw(i,k))*dtcld,0.) ncr(i,k,3) = max(ncr(i,k,3)+(nraut(i,k)-nrcol(i,k)-niacr(i,k) & -nsacr(i,k)-ngacr(i,k))*dtcld,0.) xlf = xls-xl(i,k) xlwork2 = -xls*(psdep(i,k)+pgdep(i,k)+pidep(i,k)+pigen(i,k)) & -xl(i,k)*prevp(i,k)-xlf*(piacr(i,k)+paacw(i,k) & +paacw(i,k)+pgacr(i,k)+psacr(i,k)) t(i,k) = t(i,k)-xlwork2/cpm(i,k)*dtcld else ! ! cloud water ! value = max(qmin,qci(i,k,1)) source= (praut(i,k)+pracw(i,k)+paacw(i,k)+paacw(i,k)) & *dtcld if (source.gt.value) then factor = value/source praut(i,k) = praut(i,k)*factor pracw(i,k) = pracw(i,k)*factor paacw(i,k) = paacw(i,k)*factor endif ! ! rain ! value = max(qmin,qrs(i,k,1)) source = (-paacw(i,k)-praut(i,k)+pseml(i,k)+pgeml(i,k) & -pracw(i,k)-paacw(i,k)-prevp(i,k))*dtcld if (source.gt.value) then factor = value/source praut(i,k) = praut(i,k)*factor prevp(i,k) = prevp(i,k)*factor pracw(i,k) = pracw(i,k)*factor paacw(i,k) = paacw(i,k)*factor pseml(i,k) = pseml(i,k)*factor pgeml(i,k) = pgeml(i,k)*factor endif ! ! snow ! value = max(qcrmin,qrs(i,k,2)) source=(pgacs(i,k)-pseml(i,k)-psevp(i,k))*dtcld if (source.gt.value) then factor = value/source pgacs(i,k) = pgacs(i,k)*factor psevp(i,k) = psevp(i,k)*factor pseml(i,k) = pseml(i,k)*factor endif ! ! graupel ! value = max(qcrmin,qrs(i,k,3)) source=-(pgacs(i,k)+pgevp(i,k)+pgeml(i,k))*dtcld if (source.gt.value) then factor = value/source pgacs(i,k) = pgacs(i,k)*factor pgevp(i,k) = pgevp(i,k)*factor pgeml(i,k) = pgeml(i,k)*factor endif ! ! cloud ! value = max(ncmin,ncr(i,k,2)) source = (+nraut(i,k)+nccol(i,k)+nracw(i,k)+naacw(i,k) & +naacw(i,k))*dtcld if (source.gt.value) then factor = value/source nraut(i,k) = nraut(i,k)*factor nccol(i,k) = nccol(i,k)*factor nracw(i,k) = nracw(i,k)*factor naacw(i,k) = naacw(i,k)*factor endif ! ! rain ! value = max(nrmin,ncr(i,k,3)) source = (-nraut(i,k)+nrcol(i,k)-nseml(i,k)-ngeml(i,k) & )*dtcld if (source.gt.value) then factor = value/source nraut(i,k) = nraut(i,k)*factor nrcol(i,k) = nrcol(i,k)*factor nseml(i,k) = nseml(i,k)*factor ngeml(i,k) = ngeml(i,k)*factor endif ! work2(i,k)=-(prevp(i,k)+psevp(i,k)+pgevp(i,k)) ! update q(i,k) = q(i,k)+work2(i,k)*dtcld qci(i,k,1) = max(qci(i,k,1)-(praut(i,k)+pracw(i,k) & +paacw(i,k)+paacw(i,k))*dtcld,0.) qrs(i,k,1) = max(qrs(i,k,1)+(praut(i,k)+pracw(i,k) & +prevp(i,k)+paacw(i,k)+paacw(i,k)-pseml(i,k) & -pgeml(i,k))*dtcld,0.) qrs(i,k,2) = max(qrs(i,k,2)+(psevp(i,k)-pgacs(i,k) & +pseml(i,k))*dtcld,0.) qrs(i,k,3) = max(qrs(i,k,3)+(pgacs(i,k)+pgevp(i,k) & +pgeml(i,k))*dtcld,0.) ncr(i,k,2) = max(ncr(i,k,2)+(-nraut(i,k)-nccol(i,k)-nracw(i,k) & -naacw(i,k)-naacw(i,k))*dtcld,0.) ncr(i,k,3) = max(ncr(i,k,3)+(nraut(i,k)-nrcol(i,k)+nseml(i,k) & +ngeml(i,k))*dtcld,0.) xlf = xls-xl(i,k) xlwork2 = -xl(i,k)*(prevp(i,k)+psevp(i,k)+pgevp(i,k)) & -xlf*(pseml(i,k)+pgeml(i,k)) t(i,k) = t(i,k)-xlwork2/cpm(i,k)*dtcld endif enddo enddo ! ! Inline expansion for fpvs ! qs(i,k,1) = fpvs(t(i,k),0,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c) ! qs(i,k,2) = fpvs(t(i,k),1,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c) hsub = xls hvap = xlv0 cvap = cpv ttp=t0c+0.01 dldt=cvap-cliq xa=-dldt/rv xb=xa+hvap/(rv*ttp) dldti=cvap-cice xai=-dldti/rv xbi=xai+hsub/(rv*ttp) do k = kts, kte do i = its, ite tr=ttp/t(i,k) qs(i,k,1)=psat*exp(log(tr)*(xa))*exp(xb*(1.-tr)) qs(i,k,1) = min(qs(i,k,1),0.99*p(i,k)) qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k) - qs(i,k,1)) qs(i,k,1) = max(qs(i,k,1),qmin) tr=ttp/t(i,k) if(t(i,k).lt.ttp) then qs(i,k,2)=psat*exp(log(tr)*(xai))*exp(xbi*(1.-tr)) else qs(i,k,2)=psat*exp(log(tr)*(xa))*exp(xb*(1.-tr)) endif qs(i,k,2) = min(qs(i,k,2),0.99*p(i,k)) qs(i,k,2) = ep2 * qs(i,k,2) / (p(i,k) - qs(i,k,2)) qs(i,k,2) = max(qs(i,k,2),qmin) rh(i,k,1) = max(q(i,k) / qs(i,k,1),qmin) enddo enddo ! do k = kts,kte_in do i = its,ite qrs_tmp(i,k,1) = qrs(i,k,1) qrs_tmp(i,k,2) = qrs(i,k,2) qrs_tmp(i,k,3) = qrs(i,k,3) ncr_tmp(i,k) = ncr(i,k,3) enddo enddo ! call slope_wdm6(qrs_tmp,ncr_tmp,den_tmp,denfac,t,rslope,rslopeb,rslope2, & rslope3,work1,workn,its,ite,kts,kte) do k = kts, kte do i = its, ite !----------------------------------------------------------------- ! re-compute the mean-volume drop diameter [LH A10] ! for raindrop distribution !----------------------------------------------------------------- avedia(i,k,2) = rslope(i,k,1)*((24.)**(.3333333)) !---------------------------------------------------------------- ! Nrevp_s: evaporation/condensation rate of rain [LH A14] ! (NR->NC) !---------------------------------------------------------------- if(avedia(i,k,2).le.di82) then ncr(i,k,2) = ncr(i,k,2)+ncr(i,k,3) ncr(i,k,3) = 0. !---------------------------------------------------------------- ! Prevp_s: evaporation/condensation rate of rain [LH A15] [KK 23] ! (QR->QC) !---------------------------------------------------------------- qci(i,k,1) = qci(i,k,1)+qrs(i,k,1) qrs(i,k,1) = 0. endif enddo enddo ! do k = kts, kte do i = its, ite !--------------------------------------------------------------- ! rate of change of cloud drop concentration due to CCN activation ! pcact: QV -> QC [LH 8] [KK 14] ! ncact: NCCN -> NC [LH A2] [KK 12] !--------------------------------------------------------------- if(rh(i,k,1).gt.1.) then ncact(i,k) = max(0.,((ncr(i,k,1)+ncr(i,k,2)) & *min(1.,(rh(i,k,1)/satmax)**actk) - ncr(i,k,2)))/dtcld ncact(i,k) =min(ncact(i,k),max(ncr(i,k,1),0.)/dtcld) pcact(i,k) = min(4.*pi*denr*(actr*1.E-6)**3*ncact(i,k)/ & (3.*den(i,k)),max(q(i,k),0.)/dtcld) q(i,k) = max(q(i,k)-pcact(i,k)*dtcld,0.) qci(i,k,1) = max(qci(i,k,1)+pcact(i,k)*dtcld,0.) ncr(i,k,1) = max(ncr(i,k,1)-ncact(i,k)*dtcld,0.) ncr(i,k,2) = max(ncr(i,k,2)+ncact(i,k)*dtcld,0.) t(i,k) = t(i,k)+pcact(i,k)*xl(i,k)/cpm(i,k)*dtcld endif !--------------------------------------------------------------- ! pcond:condensational/evaporational rate of cloud water [HL A46] [RH83 A6] ! if there exists additional water vapor condensated/if ! evaporation of cloud water is not enough to remove subsaturation ! (QV->QC or QC->QV) !--------------------------------------------------------------- tr=ttp/t(i,k) qs(i,k,1)=psat*exp(log(tr)*(xa))*exp(xb*(1.-tr)) qs(i,k,1) = min(qs(i,k,1),0.99*p(i,k)) qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k) - qs(i,k,1)) qs(i,k,1) = max(qs(i,k,1),qmin) work1(i,k,1) = conden(t(i,k),q(i,k),qs(i,k,1),xl(i,k),cpm(i,k)) work2(i,k) = qci(i,k,1)+work1(i,k,1) pcond(i,k) = min(max(work1(i,k,1)/dtcld,0.),max(q(i,k),0.)/dtcld) if(qci(i,k,1).gt.0. .and. work1(i,k,1).lt.0.) & pcond(i,k) = max(work1(i,k,1),-qci(i,k,1))/dtcld !---------------------------------------------------------------- ! ncevp: evpration of Cloud number concentration [LH A3] ! (NC->NCCN) !---------------------------------------------------------------- if(pcond(i,k).eq.-qci(i,k,1)/dtcld) then ncr(i,k,2) = 0. ncr(i,k,1) = ncr(i,k,1)+ncr(i,k,2) endif ! q(i,k) = q(i,k)-pcond(i,k)*dtcld qci(i,k,1) = max(qci(i,k,1)+pcond(i,k)*dtcld,0.) t(i,k) = t(i,k)+pcond(i,k)*xl(i,k)/cpm(i,k)*dtcld enddo enddo ! !---------------------------------------------------------------- ! padding for small values ! do k = kts, kte do i = its, ite if(qci(i,k,1).le.qmin) qci(i,k,1) = 0.0 if(qci(i,k,2).le.qmin) qci(i,k,2) = 0.0 if(qrs(i,k,1).ge.qcrmin .and. ncr(i,k,3) .ge. nrmin) then lamdr_tmp(i,k) = exp(log(((pidnr*ncr(i,k,3)) & /(den(i,k)*qrs(i,k,1))))*((.33333333))) if(lamdr_tmp(i,k) .le. lamdarmin) then lamdr_tmp(i,k) = lamdarmin ncr(i,k,3) = den(i,k)*qrs(i,k,1)*lamdr_tmp(i,k)**3/pidnr elseif(lamdr_tmp(i,k) .ge. lamdarmax) then lamdr_tmp(i,k) = lamdarmax ncr(i,k,3) = den(i,k)*qrs(i,k,1)*lamdr_tmp(i,k)**3/pidnr endif endif if(qci(i,k,1).ge.qmin .and. ncr(i,k,2) .ge. ncmin ) then lamdc_tmp(i,k) = exp(log(((pidnc*ncr(i,k,2)) & /(den(i,k)*qci(i,k,1))))*((.33333333))) if(lamdc_tmp(i,k) .le. lamdacmin) then lamdc_tmp(i,k) = lamdacmin ncr(i,k,2) = den(i,k)*qci(i,k,1)*lamdc_tmp(i,k)**3/pidnc elseif(lamdc_tmp(i,k) .ge. lamdacmax) then lamdc_tmp(i,k) = lamdacmax ncr(i,k,2) = den(i,k)*qci(i,k,1)*lamdc_tmp(i,k)**3/pidnc endif endif enddo enddo enddo ! big loops ! end subroutine wdm62d !------------------------------------------------------------------------------ real function rgmma(x) !------------------------------------------------------------------------------ implicit none !----------------------------------------------------------------------------- ! rgmma function: use infinite product form ! real, parameter :: euler = 0.577215664901532 real :: x, y integer :: i !------------------------------------------------------------------------------ if(x.eq.1.) then rgmma=0. else rgmma = x*exp(euler*x) do i = 1,10000 y = float(i) rgmma = rgmma*(1.000+x/y)*exp(-x/y) enddo rgmma = 1./rgmma endif ! end function rgmma ! !-------------------------------------------------------------------------- real function fpvs(t,ice,rd,rv,cvap,cliq,cice,hvap,hsub,psat,t0c) !-------------------------------------------------------------------------- implicit none !-------------------------------------------------------------------------- REAL t,rd,rv,cvap,cliq,cice,hvap,hsub,psat,t0c,dldt,xa,xb,dldti, & xai,xbi,ttp,tr INTEGER ice ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ttp=t0c+0.01 dldt=cvap-cliq xa=-dldt/rv xb=xa+hvap/(rv*ttp) dldti=cvap-cice xai=-dldti/rv xbi=xai+hsub/(rv*ttp) tr=ttp/t if(t.lt.ttp .and. ice.eq.1) then fpvs=psat*(tr**xai)*exp(xbi*(1.-tr)) else fpvs=psat*(tr**xa)*exp(xb*(1.-tr)) endif ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - END FUNCTION fpvs !------------------------------------------------------------------- SUBROUTINE wdm6init(den0,denr,dens,cl,cpv, ccn0, hail_opt, allowed_to_read) ! RAS !------------------------------------------------------------------- IMPLICIT NONE !------------------------------------------------------------------- !.... constants which may not be tunable REAL, INTENT(IN) :: den0,denr,dens,cl,cpv,ccn0 INTEGER, INTENT(IN) :: hail_opt ! RAS LOGICAL, INTENT(IN) :: allowed_to_read ! RAS13.1 define graupel parameters as graupel-like or hail-like, ! depending on namelist option IF (hail_opt .eq. 1) THEN !Hail! n0g = 4.e4 deng = 700. avtg = 285.0 bvtg = 0.8 lamdagmax = 2.e4 ELSE !Graupel! n0g = 4.e6 deng = 500 avtg = 330.0 bvtg = 0.8 lamdagmax = 6.e4 ENDIF ! pi = 4.*atan(1.) xlv1 = cl-cpv ! qc0 = 4./3.*pi*denr*r0**3.*xncr0/den0 qc1 = 4./3.*pi*denr*r0**3.*xncr1/den0 qck1 = .104*9.8*peaut/(xncr*denr)**(1./3.)/xmyu*den0**(4./3.) ! 7.03 pidnc = pi*denr/6. ! bvtr1 = 1.+bvtr bvtr2 = 2.+bvtr bvtr3 = 3.+bvtr bvtr4 = 4.+bvtr bvtr5 = 5.+bvtr bvtr6 = 6.+bvtr bvtr7 = 7.+bvtr bvtr2o5 = 2.5+.5*bvtr bvtr3o5 = 3.5+.5*bvtr g1pbr = rgmma(bvtr1) g2pbr = rgmma(bvtr2) g3pbr = rgmma(bvtr3) g4pbr = rgmma(bvtr4) ! 17.837825 g5pbr = rgmma(bvtr5) g6pbr = rgmma(bvtr6) g7pbr = rgmma(bvtr7) g5pbro2 = rgmma(bvtr2o5) g7pbro2 = rgmma(bvtr3o5) pvtr = avtr*g5pbr/24. pvtrn = avtr*g2pbr eacrr = 1.0 pacrr = pi*n0r*avtr*g3pbr*.25*eacrr precr1 = 2.*pi*1.56 precr2 = 2.*pi*.31*avtr**.5*g7pbro2 pidn0r = pi*denr*n0r pidnr = 4.*pi*denr ! xmmax = (dimax/dicon)**2 roqimax = 2.08e22*dimax**8 ! bvts1 = 1.+bvts bvts2 = 2.5+.5*bvts bvts3 = 3.+bvts bvts4 = 4.+bvts g1pbs = rgmma(bvts1) !.8875 g3pbs = rgmma(bvts3) g4pbs = rgmma(bvts4) ! 12.0786 g5pbso2 = rgmma(bvts2) pvts = avts*g4pbs/6. pacrs = pi*n0s*avts*g3pbs*.25 precs1 = 4.*n0s*.65 precs2 = 4.*n0s*.44*avts**.5*g5pbso2 pidn0s = pi*dens*n0s ! pacrc = pi*n0s*avts*g3pbs*.25*eacrc ! bvtg1 = 1.+bvtg bvtg2 = 2.5+.5*bvtg bvtg3 = 3.+bvtg bvtg4 = 4.+bvtg g1pbg = rgmma(bvtg1) g3pbg = rgmma(bvtg3) g4pbg = rgmma(bvtg4) g5pbgo2 = rgmma(bvtg2) pacrg = pi*n0g*avtg*g3pbg*.25 pvtg = avtg*g4pbg/6. precg1 = 2.*pi*n0g*.78 precg2 = 2.*pi*n0g*.31*avtg**.5*g5pbgo2 pidn0g = pi*deng*n0g ! rslopecmax = 1./lamdacmax rslopermax = 1./lamdarmax rslopesmax = 1./lamdasmax rslopegmax = 1./lamdagmax rsloperbmax = rslopermax ** bvtr rslopesbmax = rslopesmax ** bvts rslopegbmax = rslopegmax ** bvtg rslopec2max = rslopecmax * rslopecmax rsloper2max = rslopermax * rslopermax rslopes2max = rslopesmax * rslopesmax rslopeg2max = rslopegmax * rslopegmax rslopec3max = rslopec2max * rslopecmax rsloper3max = rsloper2max * rslopermax rslopes3max = rslopes2max * rslopesmax rslopeg3max = rslopeg2max * rslopegmax !+---+-----------------------------------------------------------------+ !..Set these variables needed for computing radar reflectivity. These !.. get used within radar_init to create other variables used in the !.. radar module. xam_r = PI*denr/6. xbm_r = 3. xmu_r = 1. xam_s = PI*dens/6. xbm_s = 3. xmu_s = 0. xam_g = PI*deng/6. xbm_g = 3. xmu_g = 0. call radar_init !+---+-----------------------------------------------------------------+ ! END SUBROUTINE wdm6init !------------------------------------------------------------------------------ subroutine slope_wdm6(qrs,ncr,den,denfac,t,rslope,rslopeb,rslope2,rslope3, & vt,vtn,its,ite,kts,kte) IMPLICIT NONE INTEGER :: its,ite, jts,jte, kts,kte REAL, DIMENSION( its:ite , kts:kte,3) :: & qrs, & rslope, & rslopeb, & rslope2, & rslope3, & vt REAL, DIMENSION( its:ite , kts:kte) :: & ncr, & vtn, & den, & denfac, & t REAL, PARAMETER :: t0c = 273.15 REAL, DIMENSION( its:ite , kts:kte ) :: & n0sfac REAL :: lamdar, lamdas, lamdag, x, y, z, supcol integer :: i, j, k !---------------------------------------------------------------- ! size distributions: (x=mixing ratio, y=air density): ! valid for mixing ratio > 1.e-9 kg/kg. ! ! Optimizatin : A**B => exp(log(A)*(B)) lamdar(x,y,z)= exp(log(((pidnr*z)/(x*y)))*((.33333333))) lamdas(x,y,z)= sqrt(sqrt(pidn0s*z/(x*y))) ! (pidn0s*z/(x*y))**.25 lamdag(x,y)= sqrt(sqrt(pidn0g/(x*y))) ! (pidn0g/(x*y))**.25 ! do k = kts, kte do i = its, ite supcol = t0c-t(i,k) !--------------------------------------------------------------- ! n0s: Intercept parameter for snow [m-4] [HDC 6] !--------------------------------------------------------------- n0sfac(i,k) = max(min(exp(alpha*supcol),n0smax/n0s),1.) if(qrs(i,k,1).le.qcrmin .or. ncr(i,k).le.nrmin ) then rslope(i,k,1) = rslopermax rslopeb(i,k,1) = rsloperbmax rslope2(i,k,1) = rsloper2max rslope3(i,k,1) = rsloper3max else rslope(i,k,1) = min(1./lamdar(qrs(i,k,1),den(i,k),ncr(i,k)),1.e-3) rslopeb(i,k,1) = rslope(i,k,1)**bvtr rslope2(i,k,1) = rslope(i,k,1)*rslope(i,k,1) rslope3(i,k,1) = rslope2(i,k,1)*rslope(i,k,1) endif if(qrs(i,k,2).le.qcrmin) then rslope(i,k,2) = rslopesmax rslopeb(i,k,2) = rslopesbmax rslope2(i,k,2) = rslopes2max rslope3(i,k,2) = rslopes3max else rslope(i,k,2) = 1./lamdas(qrs(i,k,2),den(i,k),n0sfac(i,k)) rslopeb(i,k,2) = rslope(i,k,2)**bvts rslope2(i,k,2) = rslope(i,k,2)*rslope(i,k,2) rslope3(i,k,2) = rslope2(i,k,2)*rslope(i,k,2) endif if(qrs(i,k,3).le.qcrmin) then rslope(i,k,3) = rslopegmax rslopeb(i,k,3) = rslopegbmax rslope2(i,k,3) = rslopeg2max rslope3(i,k,3) = rslopeg3max else rslope(i,k,3) = 1./lamdag(qrs(i,k,3),den(i,k)) rslopeb(i,k,3) = rslope(i,k,3)**bvtg rslope2(i,k,3) = rslope(i,k,3)*rslope(i,k,3) rslope3(i,k,3) = rslope2(i,k,3)*rslope(i,k,3) endif vt(i,k,1) = pvtr*rslopeb(i,k,1)*denfac(i,k) vt(i,k,2) = pvts*rslopeb(i,k,2)*denfac(i,k) vt(i,k,3) = pvtg*rslopeb(i,k,3)*denfac(i,k) vtn(i,k) = pvtrn*rslopeb(i,k,1)*denfac(i,k) if(qrs(i,k,1).le.0.0) vt(i,k,1) = 0.0 if(qrs(i,k,2).le.0.0) vt(i,k,2) = 0.0 if(qrs(i,k,3).le.0.0) vt(i,k,3) = 0.0 if(ncr(i,k).le.0.0) vtn(i,k) = 0.0 enddo enddo END subroutine slope_wdm6 !----------------------------------------------------------------------------- subroutine slope_rain(qrs,ncr,den,denfac,t,rslope,rslopeb,rslope2,rslope3, & vt,vtn,its,ite,kts,kte) IMPLICIT NONE INTEGER :: its,ite, jts,jte, kts,kte REAL, DIMENSION( its:ite , kts:kte) :: & qrs, & ncr, & rslope, & rslopeb, & rslope2, & rslope3, & vt, & vtn, & den, & denfac, & t REAL, PARAMETER :: t0c = 273.15 REAL, DIMENSION( its:ite , kts:kte ) :: & n0sfac REAL :: lamdar, x, y, z, supcol integer :: i, j, k !---------------------------------------------------------------- ! size distributions: (x=mixing ratio, y=air density): ! valid for mixing ratio > 1.e-9 kg/kg. lamdar(x,y,z)= exp(log(((pidnr*z)/(x*y)))*((.33333333))) ! do k = kts, kte do i = its, ite if(qrs(i,k).le.qcrmin .or. ncr(i,k).le.nrmin) then rslope(i,k) = rslopermax rslopeb(i,k) = rsloperbmax rslope2(i,k) = rsloper2max rslope3(i,k) = rsloper3max else rslope(i,k) = min(1./lamdar(qrs(i,k),den(i,k),ncr(i,k)),1.e-3) rslopeb(i,k) = rslope(i,k)**bvtr rslope2(i,k) = rslope(i,k)*rslope(i,k) rslope3(i,k) = rslope2(i,k)*rslope(i,k) endif vt(i,k) = pvtr*rslopeb(i,k)*denfac(i,k) vtn(i,k) = pvtrn*rslopeb(i,k)*denfac(i,k) if(qrs(i,k).le.0.0) vt(i,k) = 0.0 if(ncr(i,k).le.0.0) vtn(i,k) = 0.0 enddo enddo END subroutine slope_rain !------------------------------------------------------------------------------ subroutine slope_snow(qrs,den,denfac,t,rslope,rslopeb,rslope2,rslope3, & vt,its,ite,kts,kte) IMPLICIT NONE INTEGER :: its,ite, jts,jte, kts,kte REAL, DIMENSION( its:ite , kts:kte) :: & qrs, & rslope, & rslopeb, & rslope2, & rslope3, & vt, & den, & denfac, & t REAL, PARAMETER :: t0c = 273.15 REAL, DIMENSION( its:ite , kts:kte ) :: & n0sfac REAL :: lamdas, x, y, z, supcol integer :: i, j, k !---------------------------------------------------------------- ! size distributions: (x=mixing ratio, y=air density): ! valid for mixing ratio > 1.e-9 kg/kg. lamdas(x,y,z)= sqrt(sqrt(pidn0s*z/(x*y))) ! (pidn0s*z/(x*y))**.25 ! do k = kts, kte do i = its, ite supcol = t0c-t(i,k) !--------------------------------------------------------------- ! n0s: Intercept parameter for snow [m-4] [HDC 6] !--------------------------------------------------------------- n0sfac(i,k) = max(min(exp(alpha*supcol),n0smax/n0s),1.) if(qrs(i,k).le.qcrmin)then rslope(i,k) = rslopesmax rslopeb(i,k) = rslopesbmax rslope2(i,k) = rslopes2max rslope3(i,k) = rslopes3max else rslope(i,k) = 1./lamdas(qrs(i,k),den(i,k),n0sfac(i,k)) rslopeb(i,k) = rslope(i,k)**bvts rslope2(i,k) = rslope(i,k)*rslope(i,k) rslope3(i,k) = rslope2(i,k)*rslope(i,k) endif vt(i,k) = pvts*rslopeb(i,k)*denfac(i,k) if(qrs(i,k).le.0.0) vt(i,k) = 0.0 enddo enddo END subroutine slope_snow !---------------------------------------------------------------------------------- subroutine slope_graup(qrs,den,denfac,t,rslope,rslopeb,rslope2,rslope3, & vt,its,ite,kts,kte) IMPLICIT NONE INTEGER :: its,ite, jts,jte, kts,kte REAL, DIMENSION( its:ite , kts:kte) :: & qrs, & rslope, & rslopeb, & rslope2, & rslope3, & vt, & den, & denfac, & t REAL, PARAMETER :: t0c = 273.15 REAL, DIMENSION( its:ite , kts:kte ) :: & n0sfac REAL :: lamdag, x, y, z, supcol integer :: i, j, k !---------------------------------------------------------------- ! size distributions: (x=mixing ratio, y=air density): ! valid for mixing ratio > 1.e-9 kg/kg. lamdag(x,y)= sqrt(sqrt(pidn0g/(x*y))) ! (pidn0g/(x*y))**.25 ! do k = kts, kte do i = its, ite !--------------------------------------------------------------- ! n0s: Intercept parameter for snow [m-4] [HDC 6] !--------------------------------------------------------------- if(qrs(i,k).le.qcrmin)then rslope(i,k) = rslopegmax rslopeb(i,k) = rslopegbmax rslope2(i,k) = rslopeg2max rslope3(i,k) = rslopeg3max else rslope(i,k) = 1./lamdag(qrs(i,k),den(i,k)) rslopeb(i,k) = rslope(i,k)**bvtg rslope2(i,k) = rslope(i,k)*rslope(i,k) rslope3(i,k) = rslope2(i,k)*rslope(i,k) endif vt(i,k) = pvtg*rslopeb(i,k)*denfac(i,k) if(qrs(i,k).le.0.0) vt(i,k) = 0.0 enddo enddo END subroutine slope_graup !--------------------------------------------------------------------------------- !------------------------------------------------------------------- SUBROUTINE nislfv_rain_plmr(im,km,denl,denfacl,tkl,dzl,wwl,rql,rnl,precip,dt,id,iter,rid) !------------------------------------------------------------------- ! ! for non-iteration semi-Lagrangain forward advection for cloud ! with mass conservation and positive definite advection ! 2nd order interpolation with monotonic piecewise linear method ! this routine is under assumption of decfl < 1 for semi_Lagrangian ! ! dzl depth of model layer in meter ! wwl terminal velocity at model layer m/s ! rql cloud density*mixing ration ! precip precipitation ! dt time step ! id kind of precip: 0 test case; 1 raindrop ! iter how many time to guess mean terminal velocity: 0 pure forward. ! 0 : use departure wind for advection ! 1 : use mean wind for advection ! > 1 : use mean wind after iter-1 iterations ! rid : 1 for number 0 for mixing ratio ! ! author: hann-ming henry juang ! implemented by song-you hong ! implicit none integer im,km,id real dt real dzl(im,km),wwl(im,km),rql(im,km),rnl(im,km),precip(im) real denl(im,km),denfacl(im,km),tkl(im,km) ! integer i,k,n,m,kk,kb,kt,iter,rid real tl,tl2,qql,dql,qqd real th,th2,qqh,dqh real zsum,qsum,dim,dip,c1,con1,fa1,fa2 real allold, allnew, zz, dzamin, cflmax, decfl real dz(km), ww(km), qq(km), nr(km), wd(km), wa(km), wa2(km), was(km) real den(km), denfac(km), tk(km) real wi(km+1), zi(km+1), za(km+1) real qn(km), qr(km),tmp(km),tmp1(km),tmp2(km),tmp3(km) real dza(km+1), qa(km+1), qmi(km+1), qpi(km+1) ! precip(:) = 0.0 ! i_loop : do i=1,im ! ----------------------------------- dz(:) = dzl(i,:) qq(:) = rql(i,:) nr(:) = rnl(i,:) if(rid .eq. 1) nr(:) = rnl(i,:)/denl(i,:) ww(:) = wwl(i,:) den(:) = denl(i,:) denfac(:) = denfacl(i,:) tk(:) = tkl(i,:) ! skip for no precipitation for all layers allold = 0.0 do k=1,km allold = allold + qq(k) enddo if(allold.le.0.0) then cycle i_loop endif ! ! compute interface values zi(1)=0.0 do k=1,km zi(k+1) = zi(k)+dz(k) enddo ! ! save departure wind wd(:) = ww(:) n=1 100 continue ! plm is 2nd order, we can use 2nd order wi or 3rd order wi ! 2nd order interpolation to get wi wi(1) = ww(1) wi(km+1) = ww(km) do k=2,km wi(k) = (ww(k)*dz(k-1)+ww(k-1)*dz(k))/(dz(k-1)+dz(k)) enddo ! 3rd order interpolation to get wi fa1 = 9./16. fa2 = 1./16. wi(1) = ww(1) wi(2) = 0.5*(ww(2)+ww(1)) do k=3,km-1 wi(k) = fa1*(ww(k)+ww(k-1))-fa2*(ww(k+1)+ww(k-2)) 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 decfl = (wi(k+1)-wi(k))*dt/dz(k) if( decfl .gt. con1 ) then wi(k) = wi(k+1) - con1*dz(k)/dt endif enddo ! compute arrival point do k=1,km+1 za(k) = zi(k) - wi(k)*dt enddo ! do k=1,km dza(k) = za(k+1)-za(k) enddo dza(km+1) = zi(km+1) - za(km+1) ! ! computer deformation at arrival point do k=1,km qa(k) = qq(k)*dz(k)/dza(k) qr(k) = qa(k)/den(k) if(rid .eq. 1) qr(k) = qa(K) enddo qa(km+1) = 0.0 ! call maxmin(km,1,qa,' arrival points ') ! ! compute arrival terminal velocity, and estimate mean terminal velocity ! then back to use mean terminal velocity if( n.le.iter ) then if(rid.eq.1) then call slope_rain(nr,qr,den,denfac,tk,tmp,tmp1,tmp2,tmp3,wa,wa2,1,1,1,km) else call slope_rain(qr,nr,den,denfac,tk,tmp,tmp1,tmp2,tmp3,wa,wa2,1,1,1,km) endif if(rid.eq.1) wa(:) = wa2(:) if( n.ge.2 ) wa(1:km)=0.5*(wa(1:km)+was(1:km)) do k=1,km !#ifdef DEBUG ! print*,' slope_wsm3 ',qr(k)*1000.,den(k),denfac(k),tk(k),tmp(k),tmp1(k),tmp2(k),ww(k),wa(k) !#endif ! mean wind is average of departure and new arrival winds ww(k) = 0.5* ( wd(k)+wa(k) ) enddo was(:) = wa(:) n=n+1 go to 100 endif ! ! estimate values at arrival cell interface with monotone do k=2,km dip=(qa(k+1)-qa(k))/(dza(k+1)+dza(k)) dim=(qa(k)-qa(k-1))/(dza(k-1)+dza(k)) if( dip*dim.le.0.0 ) then qmi(k)=qa(k) qpi(k)=qa(k) else qpi(k)=qa(k)+0.5*(dip+dim)*dza(k) qmi(k)=2.0*qa(k)-qpi(k) if( qpi(k).lt.0.0 .or. qmi(k).lt.0.0 ) then qpi(k) = qa(k) qmi(k) = qa(k) endif endif enddo qpi(1)=qa(1) 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 kt=1 intp : do k=1,km kb=max(kb-1,1) kt=max(kt-1,1) ! find kb and kt if( zi(k).ge.za(km+1) ) then exit intp else find_kb : do kk=kb,km if( zi(k).le.za(kk+1) ) then kb = kk exit find_kb else cycle find_kb endif enddo find_kb find_kt : do kk=kt,km if( zi(k+1).le.za(kk) ) then kt = kk exit find_kt else cycle find_kt endif enddo find_kt kt = kt - 1 ! compute q with piecewise constant method if( kt.eq.kb ) then tl=(zi(k)-za(kb))/dza(kb) th=(zi(k+1)-za(kb))/dza(kb) tl2=tl*tl th2=th*th qqd=0.5*(qpi(kb)-qmi(kb)) qqh=qqd*th2+qmi(kb)*th qql=qqd*tl2+qmi(kb)*tl qn(k) = (qqh-qql)/(th-tl) else if( kt.gt.kb ) then tl=(zi(k)-za(kb))/dza(kb) tl2=tl*tl qqd=0.5*(qpi(kb)-qmi(kb)) qql=qqd*tl2+qmi(kb)*tl dql = qa(kb)-qql zsum = (1.-tl)*dza(kb) qsum = dql*dza(kb) if( kt-kb.gt.1 ) then do m=kb+1,kt-1 zsum = zsum + dza(m) qsum = qsum + qa(m) * dza(m) enddo endif th=(zi(k+1)-za(kt))/dza(kt) th2=th*th qqd=0.5*(qpi(kt)-qmi(kt)) dqh=qqd*th2+qmi(kt)*th zsum = zsum + th*dza(kt) qsum = qsum + dqh*dza(kt) qn(k) = qsum/zsum endif cycle intp endif ! enddo intp ! ! rain out sum_precip: do k=1,km if( za(k).lt.0.0 .and. za(k+1).lt.0.0 ) then precip(i) = precip(i) + qa(k)*dza(k) cycle sum_precip else if ( za(k).lt.0.0 .and. za(k+1).ge.0.0 ) then precip(i) = precip(i) + qa(k)*(0.0-za(k)) exit sum_precip endif exit sum_precip enddo sum_precip ! ! replace the new values rql(i,:) = qn(:) ! ! ---------------------------------- enddo i_loop ! END SUBROUTINE nislfv_rain_plmr !------------------------------------------------------------------- SUBROUTINE nislfv_rain_plm6(im,km,denl,denfacl,tkl,dzl,wwl,rql,rql2, precip1, precip2,dt,id,iter) !------------------------------------------------------------------- ! ! for non-iteration semi-Lagrangain forward advection for cloud ! with mass conservation and positive definite advection ! 2nd order interpolation with monotonic piecewise linear method ! this routine is under assumption of decfl < 1 for semi_Lagrangian ! ! dzl depth of model layer in meter ! wwl terminal velocity at model layer m/s ! rql cloud density*mixing ration ! precip precipitation ! dt time step ! id kind of precip: 0 test case; 1 raindrop ! iter how many time to guess mean terminal velocity: 0 pure forward. ! 0 : use departure wind for advection ! 1 : use mean wind for advection ! > 1 : use mean wind after iter-1 iterations ! ! author: hann-ming henry juang ! implemented by song-you hong ! implicit none integer im,km,id real dt real dzl(im,km),wwl(im,km),rql(im,km),rql2(im,km),precip(im),precip1(im),precip2(im) real denl(im,km),denfacl(im,km),tkl(im,km) ! integer i,k,n,m,kk,kb,kt,iter,ist real tl,tl2,qql,dql,qqd real th,th2,qqh,dqh real zsum,qsum,dim,dip,c1,con1,fa1,fa2 real allold, allnew, zz, dzamin, cflmax, decfl real dz(km), ww(km), qq(km), qq2(km), wd(km), wa(km), wa2(km), was(km) real den(km), denfac(km), tk(km) real wi(km+1), zi(km+1), za(km+1) real qn(km), qr(km),qr2(km),tmp(km),tmp1(km),tmp2(km),tmp3(km) real dza(km+1), qa(km+1), qa2(km+1),qmi(km+1), qpi(km+1) ! precip(:) = 0.0 precip1(:) = 0.0 precip2(:) = 0.0 ! i_loop : do i=1,im ! ----------------------------------- dz(:) = dzl(i,:) qq(:) = rql(i,:) qq2(:) = rql2(i,:) ww(:) = wwl(i,:) den(:) = denl(i,:) denfac(:) = denfacl(i,:) tk(:) = tkl(i,:) ! skip for no precipitation for all layers allold = 0.0 do k=1,km allold = allold + qq(k) + qq2(k) enddo if(allold.le.0.0) then cycle i_loop endif ! ! compute interface values zi(1)=0.0 do k=1,km zi(k+1) = zi(k)+dz(k) enddo ! ! save departure wind wd(:) = ww(:) n=1 100 continue ! plm is 2nd order, we can use 2nd order wi or 3rd order wi ! 2nd order interpolation to get wi wi(1) = ww(1) wi(km+1) = ww(km) do k=2,km wi(k) = (ww(k)*dz(k-1)+ww(k-1)*dz(k))/(dz(k-1)+dz(k)) enddo ! 3rd order interpolation to get wi fa1 = 9./16. fa2 = 1./16. wi(1) = ww(1) wi(2) = 0.5*(ww(2)+ww(1)) do k=3,km-1 wi(k) = fa1*(ww(k)+ww(k-1))-fa2*(ww(k+1)+ww(k-2)) 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 decfl = (wi(k+1)-wi(k))*dt/dz(k) if( decfl .gt. con1 ) then wi(k) = wi(k+1) - con1*dz(k)/dt endif enddo ! compute arrival point do k=1,km+1 za(k) = zi(k) - wi(k)*dt enddo ! do k=1,km dza(k) = za(k+1)-za(k) enddo dza(km+1) = zi(km+1) - za(km+1) ! ! computer deformation at arrival point do k=1,km qa(k) = qq(k)*dz(k)/dza(k) qa2(k) = qq2(k)*dz(k)/dza(k) qr(k) = qa(k)/den(k) qr2(k) = qa2(k)/den(k) enddo qa(km+1) = 0.0 qa2(km+1) = 0.0 ! call maxmin(km,1,qa,' arrival points ') ! ! compute arrival terminal velocity, and estimate mean terminal velocity ! then back to use mean terminal velocity if( n.le.iter ) then call slope_snow(qr,den,denfac,tk,tmp,tmp1,tmp2,tmp3,wa,1,1,1,km) call slope_graup(qr2,den,denfac,tk,tmp,tmp1,tmp2,tmp3,wa2,1,1,1,km) do k = 1, km tmp(k) = max((qr(k)+qr2(k)), 1.E-15) IF ( tmp(k) .gt. 1.e-15 ) THEN wa(k) = (wa(k)*qr(k) + wa2(k)*qr2(k))/tmp(k) ELSE wa(k) = 0. ENDIF enddo if( n.ge.2 ) wa(1:km)=0.5*(wa(1:km)+was(1:km)) do k=1,km !#ifdef DEBUG ! print*,' slope_wsm3 ',qr(k)*1000.,den(k),denfac(k),tk(k),tmp(k),tmp1(k),tmp2(k), & ! ww(k),wa(k) !#endif ! mean wind is average of departure and new arrival winds ww(k) = 0.5* ( wd(k)+wa(k) ) enddo was(:) = wa(:) n=n+1 go to 100 endif ist_loop : do ist = 1, 2 if (ist.eq.2) then qa(:) = qa2(:) endif ! precip(i) = 0. ! ! estimate values at arrival cell interface with monotone do k=2,km dip=(qa(k+1)-qa(k))/(dza(k+1)+dza(k)) dim=(qa(k)-qa(k-1))/(dza(k-1)+dza(k)) if( dip*dim.le.0.0 ) then qmi(k)=qa(k) qpi(k)=qa(k) else qpi(k)=qa(k)+0.5*(dip+dim)*dza(k) qmi(k)=2.0*qa(k)-qpi(k) if( qpi(k).lt.0.0 .or. qmi(k).lt.0.0 ) then qpi(k) = qa(k) qmi(k) = qa(k) endif endif enddo qpi(1)=qa(1) 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 kt=1 intp : do k=1,km kb=max(kb-1,1) kt=max(kt-1,1) ! find kb and kt if( zi(k).ge.za(km+1) ) then exit intp else find_kb : do kk=kb,km if( zi(k).le.za(kk+1) ) then kb = kk exit find_kb else cycle find_kb endif enddo find_kb find_kt : do kk=kt,km if( zi(k+1).le.za(kk) ) then kt = kk exit find_kt else cycle find_kt endif enddo find_kt kt = kt - 1 ! compute q with piecewise constant method if( kt.eq.kb ) then tl=(zi(k)-za(kb))/dza(kb) th=(zi(k+1)-za(kb))/dza(kb) tl2=tl*tl th2=th*th qqd=0.5*(qpi(kb)-qmi(kb)) qqh=qqd*th2+qmi(kb)*th qql=qqd*tl2+qmi(kb)*tl qn(k) = (qqh-qql)/(th-tl) else if( kt.gt.kb ) then tl=(zi(k)-za(kb))/dza(kb) tl2=tl*tl qqd=0.5*(qpi(kb)-qmi(kb)) qql=qqd*tl2+qmi(kb)*tl dql = qa(kb)-qql zsum = (1.-tl)*dza(kb) qsum = dql*dza(kb) if( kt-kb.gt.1 ) then do m=kb+1,kt-1 zsum = zsum + dza(m) qsum = qsum + qa(m) * dza(m) enddo endif th=(zi(k+1)-za(kt))/dza(kt) th2=th*th qqd=0.5*(qpi(kt)-qmi(kt)) dqh=qqd*th2+qmi(kt)*th zsum = zsum + th*dza(kt) qsum = qsum + dqh*dza(kt) qn(k) = qsum/zsum endif cycle intp endif ! enddo intp ! ! rain out sum_precip: do k=1,km if( za(k).lt.0.0 .and. za(k+1).lt.0.0 ) then precip(i) = precip(i) + qa(k)*dza(k) cycle sum_precip else if ( za(k).lt.0.0 .and. za(k+1).ge.0.0 ) then precip(i) = precip(i) + qa(k)*(0.0-za(k)) exit sum_precip endif exit sum_precip enddo sum_precip ! ! replace the new values if(ist.eq.1) then rql(i,:) = qn(:) precip1(i) = precip(i) else rql2(i,:) = qn(:) precip2(i) = precip(i) endif enddo ist_loop ! ! ---------------------------------- enddo i_loop ! END SUBROUTINE nislfv_rain_plm6 !+---+-----------------------------------------------------------------+ subroutine refl10cm_wdm6 (qv1d, qr1d, nr1d, qs1d, qg1d, & t1d, p1d, dBZ, kts, kte, ii, jj) IMPLICIT NONE !..Sub arguments INTEGER, INTENT(IN):: kts, kte, ii, jj REAL, DIMENSION(kts:kte), INTENT(IN):: & qv1d, qr1d, nr1d, qs1d, qg1d, t1d, p1d REAL, DIMENSION(kts:kte), INTENT(INOUT):: dBZ !..Local variables REAL, DIMENSION(kts:kte):: temp, pres, qv, rho REAL, DIMENSION(kts:kte):: rr, nr, rs, rg REAL:: temp_C DOUBLE PRECISION, DIMENSION(kts:kte):: ilamr, ilams, ilamg DOUBLE PRECISION, DIMENSION(kts:kte):: N0_r, N0_s, N0_g DOUBLE PRECISION:: lamr, lams, lamg LOGICAL, DIMENSION(kts:kte):: L_qr, L_qs, L_qg REAL, DIMENSION(kts:kte):: ze_rain, ze_snow, ze_graupel DOUBLE PRECISION:: fmelt_s, fmelt_g INTEGER:: i, k, k_0, kbot, n LOGICAL:: melti DOUBLE PRECISION:: cback, x, eta, f_d REAL, PARAMETER:: R=287. !+---+ do k = kts, kte dBZ(k) = -35.0 enddo !+---+-----------------------------------------------------------------+ !..Put column of data into local arrays. !+---+-----------------------------------------------------------------+ do k = kts, kte temp(k) = t1d(k) temp_C = min(-0.001, temp(K)-273.15) qv(k) = MAX(1.E-10, qv1d(k)) pres(k) = p1d(k) rho(k) = 0.622*pres(k)/(R*temp(k)*(qv(k)+0.622)) if (qr1d(k) .gt. 1.E-9) then rr(k) = qr1d(k)*rho(k) nr(k) = nr1d(k)*rho(k) lamr = (xam_r*xcrg(3)*xorg2*nr(k)/rr(k))**xobmr ilamr(k) = 1./lamr N0_r(k) = nr(k)*xorg2*lamr**xcre(2) L_qr(k) = .true. else rr(k) = 1.E-12 nr(k) = 1.E-12 L_qr(k) = .false. endif if (qs1d(k) .gt. 1.E-9) then rs(k) = qs1d(k)*rho(k) N0_s(k) = min(n0smax, n0s*exp(-alpha*temp_C)) lams = (xam_s*xcsg(3)*N0_s(k)/rs(k))**(1./xcse(1)) ilams(k) = 1./lams L_qs(k) = .true. else rs(k) = 1.E-12 L_qs(k) = .false. endif if (qg1d(k) .gt. 1.E-9) then rg(k) = qg1d(k)*rho(k) N0_g(k) = n0g lamg = (xam_g*xcgg(3)*N0_g(k)/rg(k))**(1./xcge(1)) ilamg(k) = 1./lamg L_qg(k) = .true. else rg(k) = 1.E-12 L_qg(k) = .false. endif enddo !+---+-----------------------------------------------------------------+ !..Locate K-level of start of melting (k_0 is level above). !+---+-----------------------------------------------------------------+ melti = .false. k_0 = kts do k = kte-1, kts, -1 if ( (temp(k).gt.273.15) .and. L_qr(k) & .and. (L_qs(k+1).or.L_qg(k+1)) ) then k_0 = MAX(k+1, k_0) melti=.true. goto 195 endif enddo 195 continue !+---+-----------------------------------------------------------------+ !..Assume Rayleigh approximation at 10 cm wavelength. Rain (all temps) !.. and non-water-coated snow and graupel when below freezing are !.. simple. Integrations of m(D)*m(D)*N(D)*dD. !+---+-----------------------------------------------------------------+ do k = kts, kte ze_rain(k) = 1.e-22 ze_snow(k) = 1.e-22 ze_graupel(k) = 1.e-22 if (L_qr(k)) ze_rain(k) = N0_r(k)*xcrg(4)*ilamr(k)**xcre(4) if (L_qs(k)) ze_snow(k) = (0.176/0.93) * (6.0/PI)*(6.0/PI) & * (xam_s/900.0)*(xam_s/900.0) & * N0_s(k)*xcsg(4)*ilams(k)**xcse(4) if (L_qg(k)) ze_graupel(k) = (0.176/0.93) * (6.0/PI)*(6.0/PI) & * (xam_g/900.0)*(xam_g/900.0) & * N0_g(k)*xcgg(4)*ilamg(k)**xcge(4) enddo !+---+-----------------------------------------------------------------+ !..Special case of melting ice (snow/graupel) particles. Assume the !.. ice is surrounded by the liquid water. Fraction of meltwater is !.. extremely simple based on amount found above the melting level. !.. Uses code from Uli Blahak (rayleigh_soak_wetgraupel and supporting !.. routines). !+---+-----------------------------------------------------------------+ if (melti .and. k_0.ge.kts+1) then do k = k_0-1, kts, -1 !..Reflectivity contributed by melting snow if (L_qs(k) .and. L_qs(k_0) ) then fmelt_s = MAX(0.005d0, MIN(1.0d0-rs(k)/rs(k_0), 0.99d0)) eta = 0.d0 lams = 1./ilams(k) do n = 1, nrbins x = xam_s * xxDs(n)**xbm_s call rayleigh_soak_wetgraupel (x,DBLE(xocms),DBLE(xobms), & fmelt_s, melt_outside_s, m_w_0, m_i_0, lamda_radar, & CBACK, mixingrulestring_s, matrixstring_s, & inclusionstring_s, hoststring_s, & hostmatrixstring_s, hostinclusionstring_s) f_d = N0_s(k)*xxDs(n)**xmu_s * DEXP(-lams*xxDs(n)) eta = eta + f_d * CBACK * simpson(n) * xdts(n) enddo ze_snow(k) = SNGL(lamda4 / (pi5 * K_w) * eta) endif !..Reflectivity contributed by melting graupel if (L_qg(k) .and. L_qg(k_0) ) then fmelt_g = MAX(0.005d0, MIN(1.0d0-rg(k)/rg(k_0), 0.99d0)) eta = 0.d0 lamg = 1./ilamg(k) do n = 1, nrbins x = xam_g * xxDg(n)**xbm_g call rayleigh_soak_wetgraupel (x,DBLE(xocmg),DBLE(xobmg), & fmelt_g, melt_outside_g, m_w_0, m_i_0, lamda_radar, & CBACK, mixingrulestring_g, matrixstring_g, & inclusionstring_g, hoststring_g, & hostmatrixstring_g, hostinclusionstring_g) f_d = N0_g(k)*xxDg(n)**xmu_g * DEXP(-lamg*xxDg(n)) eta = eta + f_d * CBACK * simpson(n) * xdtg(n) enddo ze_graupel(k) = SNGL(lamda4 / (pi5 * K_w) * eta) endif enddo endif do k = kte, kts, -1 dBZ(k) = 10.*log10((ze_rain(k)+ze_snow(k)+ze_graupel(k))*1.d18) enddo end subroutine refl10cm_wdm6 !+---+-----------------------------------------------------------------+ !----------------------------------------------------------------------- subroutine effectRad_wdm6 (t, qc, nc, qi, qs, rho, qmin, t0c, & re_qc, re_qi, re_qs, kts, kte, ii, jj) !----------------------------------------------------------------------- ! Compute radiation effective radii of cloud water, ice, and snow for ! double-moment microphysics.. ! These are entirely consistent with microphysics assumptions, not ! constant or otherwise ad hoc as is internal to most radiation ! schemes. ! Coded and implemented by Soo Ya Bae, KIAPS, January 2015. !----------------------------------------------------------------------- implicit none !..Sub arguments integer, intent(in) :: kts, kte, ii, jj real, intent(in) :: qmin real, intent(in) :: t0c real, dimension( kts:kte ), intent(in):: t real, dimension( kts:kte ), intent(in):: qc real, dimension( kts:kte ), intent(in):: nc real, dimension( kts:kte ), intent(in):: qi real, dimension( kts:kte ), intent(in):: qs real, dimension( kts:kte ), intent(in):: rho real, dimension( kts:kte ), intent(inout):: re_qc real, dimension( kts:kte ), intent(inout):: re_qi real, dimension( kts:kte ), intent(inout):: re_qs !..Local variables integer:: i,k integer :: inu_c real, dimension( kts:kte ):: ni real, dimension( kts:kte ):: rqc real, dimension( kts:kte ):: rnc real, dimension( kts:kte ):: rqi real, dimension( kts:kte ):: rni real, dimension( kts:kte ):: rqs real :: cdm2 real :: temp real :: supcol, n0sfac, lamdas real :: diai ! diameter of ice in m double precision :: lamc logical :: has_qc, has_qi, has_qs !..Minimum microphys values real, parameter :: R1 = 1.E-12 real, parameter :: R2 = 1.E-6 real, parameter :: pi = 3.1415926536 real, parameter :: bm_r = 3.0 real, parameter :: obmr = 1.0/bm_r real, parameter :: cdm = 5./3. !----------------------------------------------------------------------- has_qc = .false. has_qi = .false. has_qs = .false. cdm2 = rgmma(cdm) do k=kts,kte ! for cloud rqc(k) = max(R1, qc(k)*rho(k)) rnc(k) = max(R2, nc(k)*rho(k)) if (rqc(k).gt.R1 .and. rnc(k).gt.R2) has_qc = .true. ! for ice rqi(k) = max(R1, qi(k)*rho(k)) temp = (rho(k)*max(qi(k),qmin)) temp = sqrt(sqrt(temp*temp*temp)) ni(k) = min(max(5.38e7*temp,1.e3),1.e6) rni(k)= max(R2, ni(k)*rho(k)) if (rqi(k).gt.R1 .and. rni(k).gt.R2) has_qi = .true. ! for snow rqs(k) = max(R1, qs(k)*rho(k)) if (rqs(k).gt.R1) has_qs = .true. enddo if (has_qc) then do k=kts,kte if (rqc(k).le.R1 .or. rnc(k).le.R2) CYCLE lamc = 2.*cdm2*(pidnc*nc(k)/rqc(k))**obmr re_qc(k) = max(2.51e-6,min(sngl(1.0d0/lamc),50.e-6)) enddo endif if (has_qi) then do k=kts,kte if (rqi(k).le.R1 .or. rni(k).le.R2) CYCLE diai = 11.9*sqrt(rqi(k)/ni(k)) re_qi(k) = max(10.01e-6,min(0.75*0.163*diai,125.e-6)) enddo endif if (has_qs) then do k=kts,kte if (rqs(k).le.R1) CYCLE supcol = t0c-t(k) n0sfac = max(min(exp(alpha*supcol),n0smax/n0s),1.) lamdas = sqrt(sqrt(pidn0s*n0sfac/rqs(k))) re_qs(k) = max(25.e-6,min(0.5*(1./lamdas),999.e-6)) enddo endif end subroutine effectRad_wdm6 !----------------------------------------------------------------------- END MODULE module_mp_wdm6