!WRF:MODEL_LAYER:PHYSICS ! MODULE module_mp_lin USE module_wrf_error ! REAL , PARAMETER, PRIVATE :: RH = 1.0 ! REAL , PARAMETER, PRIVATE :: episp0 = 0.622*611.21 REAL , PARAMETER, PRIVATE :: xnor = 8.0e6 REAL , PARAMETER, PRIVATE :: xnos = 3.0e6 ! Lin ! REAL , PARAMETER, PRIVATE :: xnog = 4.0e4 ! REAL , PARAMETER, PRIVATE :: rhograul = 917. ! Hobbs REAL , PARAMETER, PRIVATE :: xnog = 4.0e6 REAL , PARAMETER, PRIVATE :: rhograul = 400. ! REAL , PARAMETER, PRIVATE :: & qi0 = 1.0e-3, ql0 = 7.0e-4, qs0 = 6.0E-4, & xmi50 = 4.8e-10, xmi40 = 2.46e-10, & constb = 0.8, constd = 0.25, & o6 = 1./6., cdrag = 0.6, & avisc = 1.49628e-6, adiffwv = 8.7602e-5, & axka = 1.4132e3, di50 = 1.0e-4, xmi = 4.19e-13, & cw = 4.187e3, vf1s = 0.78, vf2s = 0.31, & xni0 = 1.0e-2, xmnin = 1.05e-18, bni = 0.5, & ci = 2.093e3 CONTAINS !------------------------------------------------------------------- ! Lin et al., 1983, JAM, 1065-1092, and ! Rutledge and Hobbs, 1984, JAS, 2949-2972 ! May 2009 - Changes and corrections from P. Blossey (U. Washington) !------------------------------------------------------------------- SUBROUTINE lin_et_al(th & ,qv, ql, qr & ,qi, qs & ,rho, pii, p & ,dt_in & ,z,ht, dz8w & ,grav, cp, Rair, rvapor & ,XLS, XLV, XLF, rhowater, rhosnow & ,EP2,SVP1,SVP2,SVP3,SVPT0 & , RAINNC, RAINNCV & , SNOWNC, SNOWNCV & , GRAUPELNC, GRAUPELNCV, SR & ,ids,ide, jds,jde, kds,kde & ,ims,ime, jms,jme, kms,kme & ,its,ite, jts,jte, kts,kte & ! Optional ,qlsink, precr, preci, precs, precg & , F_QG,F_QNDROP & , qg, qndrop & ) !------------------------------------------------------------------- IMPLICIT NONE !------------------------------------------------------------------- ! ! Shuhua 12/17/99 ! INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde , & ims,ime, jms,jme, kms,kme , & its,ite, jts,jte, kts,kte REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), & INTENT(INOUT) :: & th, & qv, & ql, & qr ! REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), & INTENT(IN ) :: & rho, & pii, & p, & dz8w REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), & INTENT(IN ) :: z REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN) :: ht REAL, INTENT(IN ) :: dt_in, & grav, & Rair, & rvapor, & cp, & XLS, & XLV, & XLF, & rhowater, & rhosnow REAL, INTENT(IN ) :: EP2,SVP1,SVP2,SVP3,SVPT0 REAL, DIMENSION( ims:ime , jms:jme ), & INTENT(INOUT) :: RAINNC, & RAINNCV, & SR ! Optional REAL, DIMENSION( ims:ime , jms:jme ), & OPTIONAL, & INTENT(INOUT) :: SNOWNC, & SNOWNCV, & GRAUPELNC, & GRAUPELNCV REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), & OPTIONAL, & INTENT(INOUT) :: & qi, & qs, & qg, & qndrop REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), & OPTIONAL, INTENT(OUT ) :: & qlsink, & ! cloud water conversion to rain (/s) precr, & ! rain precipitation rate at all levels (kg/m2/s) preci, & ! ice precipitation rate at all levels (kg/m2/s) precs, & ! snow precipitation rate at all levels (kg/m2/s) precg ! graupel precipitation rate at all levels (kg/m2/s) LOGICAL, INTENT(IN), OPTIONAL :: F_QG, F_QNDROP ! LOCAL VAR INTEGER :: min_q, max_q REAL, DIMENSION( its:ite , jts:jte ) & :: rain, snow, graupel,ice REAL, DIMENSION( kts:kte ) :: qvz, qlz, qrz, & qiz, qsz, qgz, & thz, & tothz, rhoz, & orhoz, sqrhoz, & prez, zz, & precrz, preciz, precsz, precgz, & qndropz, & dzw, preclw LOGICAL :: flag_qg, flag_qndrop ! REAL :: dt, pptrain, pptsnow, pptgraul, rhoe_s, & gindex, pptice real :: qndropconst INTEGER :: i,j,k ! flag_qg = .false. flag_qndrop = .false. IF ( PRESENT ( f_qg ) ) flag_qg = f_qg IF ( PRESENT ( f_qndrop ) ) flag_qndrop = f_qndrop ! dt=dt_in rhoe_s=1.29 qndropconst=100.e6 !sg gindex=1.0 IF (.not.flag_qg) gindex=0. j_loop: DO j = jts, jte i_loop: DO i = its, ite ! !- write data from 3-D to 1-D ! DO k = kts, kte qvz(k)=qv(i,k,j) qlz(k)=ql(i,k,j) qrz(k)=qr(i,k,j) thz(k)=th(i,k,j) ! rhoz(k)=rho(i,k,j) orhoz(k)=1./rhoz(k) prez(k)=p(i,k,j) sqrhoz(k)=sqrt(rhoe_s*orhoz(k)) tothz(k)=pii(i,k,j) zz(k)=z(i,k,j) dzw(k)=dz8w(i,k,j) END DO IF (flag_qndrop .AND. PRESENT( qndrop )) THEN DO k = kts, kte qndropz(k)=qndrop(i,k,j) ENDDO ELSE DO k = kts, kte qndropz(k)=qndropconst ENDDO ENDIF DO k = kts, kte qiz(k)=qi(i,k,j) qsz(k)=qs(i,k,j) ENDDO IF ( flag_qg .AND. PRESENT( qg ) ) THEN DO k = kts, kte qgz(k)=qg(i,k,j) ENDDO ELSE DO k = kts, kte qgz(k)=0. ENDDO ENDIF ! pptrain=0. pptsnow=0. pptgraul=0. pptice=0. CALL clphy1d( dt, qvz, qlz, qrz, qiz, qsz, qgz, & qndropz,flag_qndrop, & thz, tothz, rhoz, orhoz, sqrhoz, & prez, zz, dzw, ht(I,J), preclw, & precrz, preciz, precsz, precgz, & pptrain, pptsnow, pptgraul, pptice, & grav, cp, Rair, rvapor, gindex, & XLS, XLV, XLF, rhowater, rhosnow, & EP2,SVP1,SVP2,SVP3,SVPT0, & kts, kte, i, j ) ! ! Precipitation from cloud microphysics -- only for one time step ! ! unit is transferred from m to mm ! rain(i,j)=pptrain snow(i,j)=pptsnow graupel(i,j)=pptgraul ice(i,j)=pptice sr(i,j)=(pptice+pptsnow+pptgraul)/(pptice+pptsnow+pptgraul+pptrain+1.e-12) ! RAINNCV(i,j)= pptrain + pptsnow + pptgraul + pptice RAINNC(i,j)=RAINNC(i,j) + pptrain + pptsnow + pptgraul + pptice IF(PRESENT(SNOWNCV))SNOWNCV(i,j)= pptsnow + pptice IF(PRESENT(SNOWNC))SNOWNC(i,j)=SNOWNC(i,j) + pptsnow + pptice IF(PRESENT(GRAUPELNCV))GRAUPELNCV(i,j)= pptgraul IF(PRESENT(GRAUPELNC))GRAUPELNC(i,j)=GRAUPELNC(i,j) + pptgraul ! !- update data from 1-D back to 3-D ! ! IF ( present(qlsink) .and. present(precr) ) THEN !sg beg DO k = kts, kte if(ql(i,k,j)>1.e-20) then qlsink(i,k,j)=-preclw(k)/ql(i,k,j) else qlsink(i,k,j)=0. endif precr(i,k,j)=precrz(k) END DO END IF !sg end DO k = kts, kte qv(i,k,j)=qvz(k) ql(i,k,j)=qlz(k) qr(i,k,j)=qrz(k) th(i,k,j)=thz(k) END DO ! IF ( flag_qndrop .AND. PRESENT( qndrop ) ) THEN !sg beg DO k = kts, kte qndrop(i,k,j)=qndropz(k) ENDDO END IF !sg end DO k = kts, kte qi(i,k,j)=qiz(k) qs(i,k,j)=qsz(k) ENDDO IF ( present(preci) ) THEN !sg beg DO k = kts, kte preci(i,k,j)=preciz(k) ENDDO END IF IF ( present(precs) ) THEN DO k = kts, kte precs(i,k,j)=precsz(k) ENDDO END IF !sg end IF ( flag_qg .AND. PRESENT( qg ) ) THEN DO k = kts, kte qg(i,k,j)=qgz(k) ENDDO IF ( present(precg) ) THEN !sg beg DO k = kts, kte precg(i,k,j)=precgz(k) ENDDO !sg end END IF ELSE !sg beg IF ( present(precg) ) precg(i,:,j)=0. !sg end ENDIF ! ENDDO i_loop ENDDO j_loop END SUBROUTINE lin_et_al !----------------------------------------------------------------------- SUBROUTINE clphy1d(dt, qvz, qlz, qrz, qiz, qsz, qgz, & qndropz,flag_qndrop, & thz, tothz, rho, orho, sqrho, & prez, zz, dzw, zsfc, preclw, & precrz, preciz, precsz, precgz, & pptrain, pptsnow, pptgraul, & pptice, grav, cp, Rair, rvapor, gindex, & XLS, XLV, XLF, rhowater, rhosnow, & EP2,SVP1,SVP2,SVP3,SVPT0, & kts, kte, i, j ) !----------------------------------------------------------------------- IMPLICIT NONE !----------------------------------------------------------------------- ! This program handles the vertical 1-D cloud micphysics !----------------------------------------------------------------------- ! avisc: constant in empirical formular for dynamic viscosity of air ! =1.49628e-6 [kg/m/s] = 1.49628e-5 [g/cm/s] ! adiffwv: constant in empirical formular for diffusivity of water ! vapor in air ! = 8.7602e-5 [kgm/s3] = 8.7602 [gcm/s3] ! axka: constant in empirical formular for thermal conductivity of air ! = 1.4132e3 [m2/s2/K] = 1.4132e7 [cm2/s2/K] ! qi0: mixing ratio threshold for cloud ice aggregation [kg/kg] ! xmi50: mass of a 50 micron ice crystal ! = 4.8e-10 [kg] =4.8e-7 [g] ! xmi40: mass of a 40 micron ice crystal ! = 2.46e-10 [kg] = 2.46e-7 [g] ! di50: diameter of a 50 micro (radius) ice crystal ! =1.0e-4 [m] ! xmi: mass of one cloud ice crystal ! =4.19e-13 [kg] = 4.19e-10 [g] ! oxmi=1.0/xmi ! ! xni0=1.0e-2 [m-3] The value given in Lin et al. is wrong.(see ! Hsie et al.(1980) and Rutledge and Hobbs(1983) ) ! bni=0.5 [K-1] ! xmnin: mass of a natural ice nucleus ! = 1.05e-18 [kg] = 1.05e-15 [g] This values is suggested by ! Hsie et al. (1980) ! = 1.0e-12 [kg] suggested by Rutlegde and Hobbs (1983) ! rhowater: density of water=1.0 g/cm3=1000.0 kg/m3 ! consta: constant in empirical formular for terminal ! velocity of raindrop ! =2115.0 [cm**(1-b)/s] = 2115.0*0.01**(1-b) [m**(1-b)/s] ! constb: constant in empirical formular for terminal ! velocity of raindrop ! =0.8 ! xnor: intercept parameter of the raindrop size distribution ! = 0.08 cm-4 = 8.0e6 m-4 ! araut: time sacle for autoconversion of cloud water to raindrops ! =1.0e-3 [s-1] ! ql0: mixing ratio threshold for cloud watercoalescence [kg/kg] ! vf1r: ventilation factors for rain =0.78 ! vf2r: ventilation factors for rain =0.31 ! rhosnow: density of snow=0.1 g/cm3=100.0 kg/m3 ! constc: constant in empirical formular for terminal ! velocity of snow ! =152.93 [cm**(1-d)/s] = 152.93*0.01**(1-d) [m**(1-d)/s] ! constd: constant in empirical formular for terminal ! velocity of snow ! =0.25 ! xnos: intercept parameter of the snowflake size distribution ! vf1s: ventilation factors for snow =0.78 ! vf2s: ventilation factors for snow =0.31 ! !---------------------------------------------------------------------- INTEGER, INTENT(IN ) :: kts, kte, i, j REAL, DIMENSION( kts:kte ), & INTENT(INOUT) :: qvz, qlz, qrz, qiz, qsz, & qndropz, & qgz, thz REAL, DIMENSION( kts:kte ), & INTENT(IN ) :: tothz, rho, orho, sqrho, & prez, zz, dzw REAL, INTENT(IN ) :: dt, grav, cp, Rair, rvapor, & XLS, XLV, XLF, rhowater, & rhosnow,EP2,SVP1,SVP2,SVP3,SVPT0 REAL, DIMENSION( kts:kte ), INTENT(OUT) :: preclw, & precrz, preciz, precsz, precgz REAL, INTENT(INOUT) :: pptrain, pptsnow, pptgraul, pptice REAL, INTENT(IN ) :: zsfc logical, intent(in) :: flag_qndrop !sg ! local vars REAL :: obp4, bp3, bp5, bp6, odp4, & dp3, dp5, dp5o2 ! temperary vars REAL :: tmp, tmp0, tmp1, tmp2,tmp3, & tmp4,delta2,delta3, delta4, & tmpa,tmpb,tmpc,tmpd,alpha1, & qic, abi,abr, abg, odtberg, & vti50,eiw,eri,esi,esr, esw, & erw,delrs,term0,term1,araut, & constg2, vf1r, vf2r,alpha2, & Ap, Bp, egw, egi, egs, egr, & constg, gdelta4, g1sdelt4, & factor, tmp_r, tmp_s,tmp_g, & qlpqi, rsat, a1, a2, xnin INTEGER :: k ! REAL, DIMENSION( kts:kte ) :: oprez, tem, temcc, theiz, qswz, & qsiz, qvoqswz, qvoqsiz, qvzodt, & qlzodt, qizodt, qszodt, qrzodt, & qgzodt REAL, DIMENSION( kts:kte ) :: psnow, psaut, psfw, psfi, praci, & piacr, psaci, psacw, psdep, pssub, & pracs, psacr, psmlt, psmltevp, & prain, praut, pracw, prevp, pvapor, & pclw, pladj, pcli, pimlt, pihom, & pidw, piadj, pgraupel, pgaut, & pgfr, pgacw, pgaci, pgacr, pgacs, & pgacip,pgacrp,pgacsp,pgwet, pdry, & pgsub, pgdep, pgmlt, pgmltevp, & qschg, qgchg ! REAL, DIMENSION( kts:kte ) :: qvsbar, rs0, viscmu, visc, diffwv, & schmidt, xka REAL, DIMENSION( kts:kte ) :: vtr, vts, vtg, & vtrold, vtsold, vtgold, vtiold, & xlambdar, xlambdas, xlambdag, & olambdar, olambdas, olambdag REAL :: episp0k, dtb, odtb, pi, pio4, & pio6, oxLf, xLvocp, xLfocp, consta, & constc, ocdrag, gambp4, gamdp4, & gam4pt5, Cpor, oxmi, gambp3, gamdp3,& gambp6, gam3pt5, gam2pt75, gambp5o2,& gamdp5o2, cwoxlf, ocp, xni50, es ! REAL :: qvmin=1.e-20 REAL :: gindex REAL :: temc1,save1,save2,xni50mx ! for terminal velocity flux INTEGER :: min_q, max_q REAL :: t_del_tv, del_tv, flux, fluxin, fluxout ,tmpqrz LOGICAL :: notlast ! REAL :: tmp_tem, tmp_temcc !bloss ! real :: liqconc, dis, beta, kappa, p0, xc, capn,rhocgs !sg: begin ! liqconc = liquid water content (g cm^-3) ! capn = droplet number concentration (# cm^-3) ! dis = relative dispersion (dimensionless) between 0.2 and 1. ! Written by Yangang Liu based on Liu et al., GRL 32, 2005. ! Autoconversion rate p = p0 * (threshold function) ! p0 = "base" autoconversion rate (g cm^-3 s^-1) ! kappa = constant in Long kernel = [kappa2 * (3/(4*pi*rhow))^3] in Liu papers ! beta = Condensation rate constant = (beta6)^6 in Liu papers ! xc = Normalized critical mass ! *********************************************************** if(flag_qndrop)then dis = 0.5 ! droplet dispersion, set to 0.5 per SG 8-Nov-2006 ! Give empirical constants kappa=1.1d10 ! Calculate Condensation rate constant beta = (1.0d0+3.0d0*dis**2)*(1.0d0+4.0d0*dis**2)* & (1.0d0+5.0d0*dis**2)/((1.0d0+dis**2)*(1.0d0+2.0d0*dis**2)) endif !sg: end dtb=dt odtb=1./dtb pi=acos(-1.) pio4=acos(-1.)/4. pio6=acos(-1.)/6. ocp=1./cp oxLf=1./xLf xLvocp=xLv/cp xLfocp=xLf/cp consta=2115.0*0.01**(1-constb) constc=152.93*0.01**(1-constd) ocdrag=1./Cdrag ! episp0k=RH*episp0 episp0k=RH*ep2*1000.*svp1 ! gambp4=ggamma(constb+4.) gamdp4=ggamma(constd+4.) gam4pt5=ggamma(4.5) Cpor=cp/Rair oxmi=1.0/xmi gambp3=ggamma(constb+3.) gamdp3=ggamma(constd+3.) gambp6=ggamma(constb+6) gam3pt5=ggamma(3.5) gam2pt75=ggamma(2.75) gambp5o2=ggamma((constb+5.)/2.) gamdp5o2=ggamma((constd+5.)/2.) cwoxlf=cw/xlf delta2=0. delta3=0. delta4=0. ! !----------------------------------------------------------------------- ! oprez 1./prez ( prez : pressure) ! qsw saturated mixing ratio on water surface ! qsi saturated mixing ratio on ice surface ! episp0k RH*e*saturated pressure at 273.15 K ! qvoqsw qv/qsw ! qvoqsi qv/qsi ! qvzodt qv/dt ! qlzodt ql/dt ! qizodt qi/dt ! qszodt qs/dt ! qrzodt qr/dt ! qgzodt qg/dt ! ! temcc temperature in dregee C ! obp4=1.0/(constb+4.0) bp3=constb+3.0 bp5=constb+5.0 bp6=constb+6.0 odp4=1.0/(constd+4.0) dp3=constd+3.0 dp5=constd+5.0 dp5o2=0.5*(constd+5.0) ! do k=kts,kte oprez(k)=1./prez(k) enddo do k=kts,kte qlz(k)=amax1( 0.0,qlz(k) ) qiz(k)=amax1( 0.0,qiz(k) ) qvz(k)=amax1( qvmin,qvz(k) ) qsz(k)=amax1( 0.0,qsz(k) ) qrz(k)=amax1( 0.0,qrz(k) ) qgz(k)=amax1( 0.0,qgz(k) ) qndropz(k)=amax1( 0.0,qndropz(k) ) !sg ! tem(k)=thz(k)*tothz(k) temcc(k)=tem(k)-273.15 ! ! qswz(k)=episp0k*oprez(k)* & ! exp( svp2*temcc(k)/(tem(k)-svp3) ) es=1000.*svp1*exp( svp2*temcc(k)/(tem(k)-svp3) ) qswz(k)=ep2*es/(prez(k)-es) if (tem(k) .lt. 273.15 ) then ! qsiz(k)=episp0k*oprez(k)* & ! exp( 21.8745584*(tem(k)-273.16)/(tem(k)-7.66) ) es=1000.*svp1*exp( 21.8745584*(tem(k)-273.16)/(tem(k)-7.66) ) qsiz(k)=ep2*es/(prez(k)-es) if (temcc(k) .lt. -40.0) qswz(k)=qsiz(k) else qsiz(k)=qswz(k) endif ! qvoqswz(k)=qvz(k)/qswz(k) qvoqsiz(k)=qvz(k)/qsiz(k) theiz(k)=thz(k)+(xlvocp*qvz(k)-xlfocp*qiz(k))/tothz(k) enddo ! ! !----------------------------------------------------------------------- ! In this simple stable cloud parameterization scheme, only five ! forms of water substance (water vapor, cloud water, cloud ice, ! rain and snow are considered. The prognostic variables are total ! water (qp),cloud water (ql), and cloud ice (qi). Rain and snow are ! diagnosed following Nagata and Ogura, 1991, MWR, 1309-1337. Eq (A7). ! the micro physics are based on (1) Hsie et al.,1980, JAM, 950-977 ; ! (2) Lin et al., 1983, JAM, 1065-1092 ; (3) Rutledge and Hobbs, 1983, ! JAS, 1185-1206 ; (4) Rutledge and Hobbs, 1984, JAS, 2949-2972. !----------------------------------------------------------------------- ! ! rhowater: density of water=1.0 g/cm3=1000.0 kg/m3 ! rhosnow: density of snow=0.1 g/cm3=100.0 kg/m3 ! xnor: intercept parameter of the raindrop size distribution ! = 0.08 cm-4 = 8.0e6 m-4 ! xnos: intercept parameter of the snowflake size distribution ! = 0.03 cm-4 = 3.0e6 m-4 ! xnog: intercept parameter of the graupel size distribution ! = 4.0e-4 cm-4 = 4.0e4 m-4 ! consta: constant in empirical formular for terminal ! velocity of raindrop ! =2115.0 [cm**(1-b)/s] = 2115.0*0.01**(1-b) [m**(1-b)/s] ! constb: constant in empirical formular for terminal ! velocity of raindrop ! =0.8 ! constc: constant in empirical formular for terminal ! velocity of snow ! =152.93 [cm**(1-d)/s] = 152.93*0.01**(1-d) [m**(1-d)/s] ! constd: constant in empirical formular for terminal ! velocity of snow ! =0.25 ! avisc: constant in empirical formular for dynamic viscosity of air ! =1.49628e-6 [kg/m/s] = 1.49628e-5 [g/cm/s] ! adiffwv: constant in empirical formular for diffusivity of water ! vapor in air ! = 8.7602e-5 [kgm/s3] = 8.7602 [gcm/s3] ! axka: constant in empirical formular for thermal conductivity of air ! = 1.4132e3 [m2/s2/K] = 1.4132e7 [cm2/s2/K] ! qi0: mixing ratio threshold for cloud ice aggregation [kg/kg] ! = 1.0e-3 g/g = 1.0e-3 kg/gk ! ql0: mixing ratio threshold for cloud watercoalescence [kg/kg] ! = 2.0e-3 g/g = 2.0e-3 kg/gk ! qs0: mixing ratio threshold for snow aggregation ! = 6.0e-4 g/g = 6.0e-4 kg/gk ! xmi50: mass of a 50 micron ice crystal ! = 4.8e-10 [kg] =4.8e-7 [g] ! xmi40: mass of a 40 micron ice crystal ! = 2.46e-10 [kg] = 2.46e-7 [g] ! di50: diameter of a 50 micro (radius) ice crystal ! =1.0e-4 [m] ! xmi: mass of one cloud ice crystal ! =4.19e-13 [kg] = 4.19e-10 [g] ! oxmi=1.0/xmi ! ! if gindex=1.0 include graupel ! if gindex=0. no graupel ! ! do k=kts,kte psnow(k)=0.0 psaut(k)=0.0 psfw(k)=0.0 psfi(k)=0.0 praci(k)=0.0 piacr(k)=0.0 psaci(k)=0.0 psacw(k)=0.0 psdep(k)=0.0 pssub(k)=0.0 pracs(k)=0.0 psacr(k)=0.0 psmlt(k)=0.0 psmltevp(k)=0.0 ! prain(k)=0.0 praut(k)=0.0 pracw(k)=0.0 prevp(k)=0.0 ! pvapor(k)=0.0 ! pclw(k)=0.0 preclw(k)=0.0 !sg pladj(k)=0.0 ! pcli(k)=0.0 pimlt(k)=0.0 pihom(k)=0.0 pidw(k)=0.0 piadj(k)=0.0 enddo ! !!! graupel ! do k=kts,kte pgraupel(k)=0.0 pgaut(k)=0.0 pgfr(k)=0.0 pgacw(k)=0.0 pgaci(k)=0.0 pgacr(k)=0.0 pgacs(k)=0.0 pgacip(k)=0.0 pgacrP(k)=0.0 pgacsp(k)=0.0 pgwet(k)=0.0 pdry(k)=0.0 pgsub(k)=0.0 pgdep(k)=0.0 pgmlt(k)=0.0 pgmltevp(k)=0.0 qschg(k)=0. qgchg(k)=0. enddo ! ! ! Set rs0=episp0*oprez(k) ! episp0=e*saturated pressure at 273.15 K ! e = 0.622 ! DO k=kts,kte rs0(k)=ep2*1000.*svp1/(prez(k)-1000.*svp1) END DO ! !*********************************************************************** ! Calculate precipitation fluxes due to terminal velocities. !*********************************************************************** ! !- Calculate termianl velocity (vt?) of precipitation q?z !- Find maximum vt? to determine the small delta t ! !-- rain ! t_del_tv=0. del_tv=dtb notlast=.true. DO while (notlast) ! min_q=kte max_q=kts-1 ! do k=kts,kte-1 if (qrz(k) .gt. 1.0e-8) then min_q=min0(min_q,k) max_q=max0(max_q,k) tmp1=sqrt(pi*rhowater*xnor/rho(k)/qrz(k)) tmp1=sqrt(tmp1) vtrold(k)=o6*consta*gambp4*sqrho(k)/tmp1**constb if (k .eq. 1) then del_tv=amin1(del_tv,0.9*(zz(k)-zsfc)/vtrold(k)) else del_tv=amin1(del_tv,0.9*(zz(k)-zz(k-1))/vtrold(k)) endif else vtrold(k)=0. endif enddo if (max_q .ge. min_q) then ! !- Check if the summation of the small delta t >= big delta t ! (t_del_tv) (del_tv) (dtb) t_del_tv=t_del_tv+del_tv ! if ( t_del_tv .ge. dtb ) then notlast=.false. del_tv=dtb+del_tv-t_del_tv endif ! fluxin=0. do k=max_q,min_q,-1 fluxout=rho(k)*vtrold(k)*qrz(k) flux=(fluxin-fluxout)/rho(k)/dzw(k) tmpqrz=qrz(k) qrz(k)=qrz(k)+del_tv*flux fluxin=fluxout enddo if (min_q .eq. 1) then pptrain=pptrain+fluxin*del_tv else qrz(min_q-1)=qrz(min_q-1)+del_tv* & fluxin/rho(min_q-1)/dzw(min_q-1) endif ! else notlast=.false. endif ENDDO ! !-- snow ! t_del_tv=0. del_tv=dtb notlast=.true. DO while (notlast) ! min_q=kte max_q=kts-1 ! do k=kts,kte-1 if (qsz(k) .gt. 1.0e-8) then min_q=min0(min_q,k) max_q=max0(max_q,k) tmp1=sqrt(pi*rhosnow*xnos/rho(k)/qsz(k)) tmp1=sqrt(tmp1) vtsold(k)=o6*constc*gamdp4*sqrho(k)/tmp1**constd if (k .eq. 1) then del_tv=amin1(del_tv,0.9*(zz(k)-zsfc)/vtsold(k)) else del_tv=amin1(del_tv,0.9*(zz(k)-zz(k-1))/vtsold(k)) endif else vtsold(k)=0. endif enddo if (max_q .ge. min_q) then ! ! !- Check if the summation of the small delta t >= big delta t ! (t_del_tv) (del_tv) (dtb) t_del_tv=t_del_tv+del_tv if ( t_del_tv .ge. dtb ) then notlast=.false. del_tv=dtb+del_tv-t_del_tv endif ! fluxin=0. do k=max_q,min_q,-1 fluxout=rho(k)*vtsold(k)*qsz(k) flux=(fluxin-fluxout)/rho(k)/dzw(k) qsz(k)=qsz(k)+del_tv*flux qsz(k)=amax1(0.,qsz(k)) fluxin=fluxout enddo if (min_q .eq. 1) then pptsnow=pptsnow+fluxin*del_tv else qsz(min_q-1)=qsz(min_q-1)+del_tv* & fluxin/rho(min_q-1)/dzw(min_q-1) endif ! else notlast=.false. endif ENDDO ! !-- grauupel ! t_del_tv=0. del_tv=dtb notlast=.true. ! DO while (notlast) ! min_q=kte max_q=kts-1 ! do k=kts,kte-1 if (qgz(k) .gt. 1.0e-8) then min_q=min0(min_q,k) max_q=max0(max_q,k) tmp1=sqrt(pi*rhograul*xnog/rho(k)/qgz(k)) tmp1=sqrt(tmp1) term0=sqrt(4.*grav*rhograul*0.33334/rho(k)/cdrag) vtgold(k)=o6*gam4pt5*term0*sqrt(1./tmp1) if (k .eq. 1) then del_tv=amin1(del_tv,0.9*(zz(k)-zsfc)/vtgold(k)) else del_tv=amin1(del_tv,0.9*(zz(k)-zz(k-1))/vtgold(k)) endif else vtgold(k)=0. endif enddo if (max_q .ge. min_q) then ! ! !- Check if the summation of the small delta t >= big delta t ! (t_del_tv) (del_tv) (dtb) t_del_tv=t_del_tv+del_tv if ( t_del_tv .ge. dtb ) then notlast=.false. del_tv=dtb+del_tv-t_del_tv endif ! fluxin=0. do k=max_q,min_q,-1 fluxout=rho(k)*vtgold(k)*qgz(k) flux=(fluxin-fluxout)/rho(k)/dzw(k) qgz(k)=qgz(k)+del_tv*flux qgz(k)=amax1(0.,qgz(k)) fluxin=fluxout enddo if (min_q .eq. 1) then pptgraul=pptgraul+fluxin*del_tv else qgz(min_q-1)=qgz(min_q-1)+del_tv* & fluxin/rho(min_q-1)/dzw(min_q-1) endif ! else notlast=.false. endif ! ENDDO ! !-- cloud ice (03/21/02) follow Vaughan T.J. Phillips at GFDL ! t_del_tv=0. del_tv=dtb notlast=.true. ! DO while (notlast) ! min_q=kte max_q=kts-1 ! do k=kts,kte-1 if (qiz(k) .gt. 1.0e-8) then min_q=min0(min_q,k) max_q=max0(max_q,k) vtiold(k)= 3.29 * (rho(k)* qiz(k))** 0.16 ! Heymsfield and Donner if (k .eq. 1) then del_tv=amin1(del_tv,0.9*(zz(k)-zsfc)/vtiold(k)) else del_tv=amin1(del_tv,0.9*(zz(k)-zz(k-1))/vtiold(k)) endif else vtiold(k)=0. endif enddo if (max_q .ge. min_q) then ! ! !- Check if the summation of the small delta t >= big delta t ! (t_del_tv) (del_tv) (dtb) t_del_tv=t_del_tv+del_tv if ( t_del_tv .ge. dtb ) then notlast=.false. del_tv=dtb+del_tv-t_del_tv endif fluxin=0. do k=max_q,min_q,-1 fluxout=rho(k)*vtiold(k)*qiz(k) flux=(fluxin-fluxout)/rho(k)/dzw(k) qiz(k)=qiz(k)+del_tv*flux qiz(k)=amax1(0.,qiz(k)) fluxin=fluxout enddo if (min_q .eq. 1) then pptice=pptice+fluxin*del_tv else qiz(min_q-1)=qiz(min_q-1)+del_tv* & fluxin/rho(min_q-1)/dzw(min_q-1) endif ! else notlast=.false. endif ! ENDDO do k=kts,kte-1 !sg beg precrz(k)=rho(k)*vtrold(k)*qrz(k) preciz(k)=rho(k)*vtiold(k)*qiz(k) precsz(k)=rho(k)*vtsold(k)*qsz(k) precgz(k)=rho(k)*vtgold(k)*qgz(k) enddo !sg end precrz(kte)=0. !wig - top level never set for vtXold vars preciz(kte)=0. !wig precsz(kte)=0. !wig precgz(kte)=0. !wig ! Microphysics processes ! DO 2000 k=kts,kte ! qvzodt(k)=amax1( 0.0,odtb*qvz(k) ) qlzodt(k)=amax1( 0.0,odtb*qlz(k) ) qizodt(k)=amax1( 0.0,odtb*qiz(k) ) qszodt(k)=amax1( 0.0,odtb*qsz(k) ) qrzodt(k)=amax1( 0.0,odtb*qrz(k) ) qgzodt(k)=amax1( 0.0,odtb*qgz(k) ) !*********************************************************************** !***** diagnose mixing ratios (qrz,qsz), terminal ***** !***** velocities (vtr,vts), and slope parameters in size ***** !***** distribution(xlambdar,xlambdas) of rain and snow ***** !***** follows Nagata and Ogura, 1991, MWR, 1309-1337. Eq (A7) ***** !*********************************************************************** ! !**** assuming no cloud water can exist in the top two levels due to !**** radiation consideration ! !! if !! unsaturated, !! no cloud water, rain, ice, snow and graupel !! then !! skip these processes and jump to line 2000 ! ! tmp=qiz(k)+qlz(k)+qsz(k)+qrz(k)+qgz(k)*gindex if( qvz(k)+qlz(k)+qiz(k) .lt. qsiz(k) & .and. tmp .eq. 0.0 ) go to 2000 !! calculate terminal velocity of rain ! if (qrz(k) .gt. 1.0e-8) then tmp1=sqrt(pi*rhowater*xnor*orho(k)/qrz(k)) xlambdar(k)=sqrt(tmp1) olambdar(k)=1.0/xlambdar(k) vtrold(k)=o6*consta*gambp4*sqrho(k)*olambdar(k)**constb else vtrold(k)=0. olambdar(k)=0. endif ! ! if (qrz(k) .gt. 1.0e-12) then if (qrz(k) .gt. 1.0e-8) then tmp1=sqrt(pi*rhowater*xnor*orho(k)/qrz(k)) xlambdar(k)=sqrt(tmp1) olambdar(k)=1.0/xlambdar(k) vtr(k)=o6*consta*gambp4*sqrho(k)*olambdar(k)**constb else vtr(k)=0. olambdar(k)=0. endif ! !! calculate terminal velocity of snow ! if (qsz(k) .gt. 1.0e-8) then tmp1=sqrt(pi*rhosnow*xnos*orho(k)/qsz(k)) xlambdas(k)=sqrt(tmp1) olambdas(k)=1.0/xlambdas(k) vtsold(k)=o6*constc*gamdp4*sqrho(k)*olambdas(k)**constd else vtsold(k)=0. olambdas(k)=0. endif ! ! if (qsz(k) .gt. 1.0e-12) then if (qsz(k) .gt. 1.0e-8) then tmp1=sqrt(pi*rhosnow*xnos*orho(k)/qsz(k)) xlambdas(k)=sqrt(tmp1) olambdas(k)=1.0/xlambdas(k) vts(k)=o6*constc*gamdp4*sqrho(k)*olambdas(k)**constd else vts(k)=0. olambdas(k)=0. endif ! !! calculate terminal velocity of graupel ! if (qgz(k) .gt. 1.0e-8) then tmp1=sqrt( pi*rhograul*xnog*orho(k)/qgz(k)) xlambdag(k)=sqrt(tmp1) olambdag(k)=1.0/xlambdag(k) term0=sqrt(4.*grav*rhograul*0.33334*orho(k)*ocdrag) vtgold(k)=o6*gam4pt5*term0*sqrt(olambdag(k)) else vtgold(k)=0. olambdag(k)=0. endif ! ! if (qgz(k) .gt. 1.0e-12) then if (qgz(k) .gt. 1.0e-8) then tmp1=sqrt( pi*rhograul*xnog*orho(k)/qgz(k)) xlambdag(k)=sqrt(tmp1) olambdag(k)=1.0/xlambdag(k) term0=sqrt(4.*grav*rhograul*0.33334*orho(k)*ocdrag) vtg(k)=o6*gam4pt5*term0*sqrt(olambdag(k)) else vtg(k)=0. olambdag(k)=0. endif ! !*********************************************************************** !***** compute viscosity,difusivity,thermal conductivity, and ****** !***** Schmidt number ****** !*********************************************************************** !c------------------------------------------------------------------ !c viscmu: dynamic viscosity of air kg/m/s !c visc: kinematic viscosity of air = viscmu/rho (m2/s) !c avisc=1.49628e-6 kg/m/s=1.49628e-5 g/cm/s !c viscmu=1.718e-5 kg/m/s in RH !c diffwv: Diffusivity of water vapor in air !c adiffwv = 8.7602e-5 (8.794e-5 in MM5) kgm/s3 !c = 8.7602 (8.794 in MM5) gcm/s3 !c diffwv(k)=2.26e-5 m2/s !c schmidt: Schmidt number=visc/diffwv !c xka: thermal conductivity of air J/m/s/K (Kgm/s3/K) !c xka(k)=2.43e-2 J/m/s/K in RH !c axka=1.4132e3 (1.414e3 in MM5) m2/s2/k = 1.4132e7 cm2/s2/k !c------------------------------------------------------------------ viscmu(k)=avisc*tem(k)**1.5/(tem(k)+120.0) visc(k)=viscmu(k)*orho(k) diffwv(k)=adiffwv*tem(k)**1.81*oprez(k) schmidt(k)=visc(k)/diffwv(k) xka(k)=axka*viscmu(k) if (tem(k) .lt. 273.15) then ! !*********************************************************************** !********* snow production processes for T < 0 C ********** !*********************************************************************** !c !c (1) ICE CRYSTAL AGGREGATION TO SNOW (Psaut): Lin (21) !c! psaut=alpha1*(qi-qi0) !c! alpha1=1.0e-3*exp(0.025*(T-T0)) !c ! alpha1=1.0e-3*exp( 0.025*temcc(k) ) alpha1=1.0e-3*exp( 0.025*temcc(k) ) ! if(temcc(k) .lt. -20.0) then tmp1=-7.6+4.0*exp( -0.2443e-3*(abs(temcc(k))-20)**2.455 ) qic=1.0e-3*exp(tmp1)*orho(k) else qic=qi0 end if !testing ! tmp1=amax1( 0.0,alpha1*(qiz(k)-qic) ) ! psaut(k)=amin1( tmp1,qizodt(k) ) tmp1=odtb*(qiz(k)-qic)*(1.0-exp(-alpha1*dtb)) psaut(k)=amax1( 0.0,tmp1 ) !c !c (2) BERGERON PROCESS TRANSFER OF CLOUD WATER TO SNOW (Psfw) !c this process only considered when -31 C < T < 0 C !c Lin (33) and Hsie (17) !c !c! !c! parama1 and parama2 functions must be user supplied !c! ! testing if( qlz(k) .gt. 1.0e-10 ) then temc1=amax1(-30.99,temcc(k)) ! print*,'temc1',temc1,qlz(k) a1=parama1( temc1 ) a2=parama2( temc1 ) tmp1=1.0-a2 !! change unit from cgs to mks a1=a1*0.001**tmp1 !c! dtberg is the time needed for a crystal to grow from 40 to 50 um !c ! odtberg=1.0/dtberg odtberg=(a1*tmp1)/(xmi50**tmp1-xmi40**tmp1) ! !c! compute terminal velocity of a 50 micron ice cystal ! vti50=constc*di50**constd*sqrho(k) ! eiw=1.0 save1=a1*xmi50**a2 save2=0.25*pi*eiw*rho(k)*di50*di50*vti50 ! tmp2=( save1 + save2*qlz(k) ) ! !! maximum number of 50 micron crystals limited by the amount !! of supercool water ! xni50mx=qlzodt(k)/tmp2 ! !! number of 50 micron crystals produced ! ! xni50=qiz(k)*( 1.0-exp(-dtb*odtberg) )/xmi50 xni50=amin1(xni50,xni50mx) ! tmp3=odtb*tmp2/save2*( 1.0-exp(-save2*xni50*dtb) ) psfw(k)=amin1( tmp3,qlzodt(k) ) !testing ! psfw(k)=0. !0915 if( temcc(k).gt.-30.99 ) then !0915 a1=parama1( temcc(k) ) !0915 a2=parama2( temcc(k) ) !0915 tmp1=1.0-a2 !! change unit from cgs to mks !0915 a1=a1*0.001**tmp1 !c! dtberg is the time needed for a crystal to grow from 40 to 50 um !c! odtberg=1.0/dtberg !0915 odtberg=(a1*tmp1)/(xmi50**tmp1-xmi40**tmp1) !c! number of 50 micron crystals produced !0915 xni50=qiz(k)*dtb*odtberg/xmi50 !c! need to calculate the terminal velocity of a 50 micron !c! ice cystal !0915 vti50=constc*di50**constd*sqrho(k) !0915 eiw=1.0 !0915 tmp2=xni50*( a1*xmi50**a2 + & !0915 0.25*qlz(k)*pi*eiw*rho(k)*di50*di50*vti50 ) !0915 psfw(k)=amin1( tmp2,qlzodt(k) ) !0915 psfw(k)=0. !c !c (3) REDUCTION OF CLOUD ICE BY BERGERON PROCESS (Psfi): Lin (34) !c this process only considered when -31 C < T < 0 C !c tmp1=xni50*xmi50-psfw(k) psfi(k)=amin1(tmp1,qizodt(k)) ! testing ! psfi(k)=0. end if ! !0915 tmp1=qiz(k)*odtberg !0915 psfi(k)=amin1(tmp1,qizodt(k)) ! testing !0915 psfi(k)=0. !0915 end if ! if(qrz(k) .le. 0.0) go to 1000 ! ! Processes (4) and (5) only need when qrz > 0.0 ! !c !c (4) CLOUD ICE ACCRETION BY RAIN (Praci): Lin (25) !c may produce snow or graupel !c eri=1.0 !0915 tmp1=qiz(k)*pio4*eri*xnor*consta*sqrho(k) !0915 tmp2=tmp1*gambp3*olambdar(k)**bp3 !0915 praci(k)=amin1( tmp2,qizodt(k) ) save1=pio4*eri*xnor*consta*sqrho(k) tmp1=save1*gambp3*olambdar(k)**bp3 praci(k)=qizodt(k)*( 1.0-exp(-tmp1*dtb) ) !c !c (5) RAIN ACCRETION BY CLOUD ICE (Piacr): Lin (26) !c !0915 tmp2=tmp1*rho(k)*pio6*rhowater*gambp6*oxmi* & !0915 olambdar(k)**bp6 !0915 piacr(k)=amin1( tmp2,qrzodt(k) ) tmp2=qiz(k)*save1*rho(k)*pio6*rhowater*gambp6*oxmi* & olambdar(k)**bp6 piacr(k)=amin1( tmp2,qrzodt(k) ) ! 1000 continue ! if(qsz(k) .le. 0.0) go to 1200 ! ! Compute the following processes only when qsz > 0.0 ! !c !c (6) ICE CRYSTAL ACCRETION BY SNOW (Psaci): Lin (22) !c esi=exp( 0.025*temcc(k) ) save1=pio4*xnos*constc*gamdp3*sqrho(k)* & olambdas(k)**dp3 tmp1=esi*save1 psaci(k)=qizodt(k)*( 1.0-exp(-tmp1*dtb) ) !0915 tmp1=pio4*xnos*constc*gamdp3*sqrho(k)* & !0915 olambdas(k)**dp3 !0915 tmp2=qiz(k)*esi*tmp1 !0915 psaci(k)=amin1( tmp2,qizodt(k) ) !c !c (7) CLOUD WATER ACCRETION BY SNOW (Psacw): Lin (24) !c esw=1.0 tmp1=esw*save1 psacw(k)=qlzodt(K)*( 1.0-exp(-tmp1*dtb) ) !0915 tmp2=qlz(k)*esw*tmp1 !0915 psacw(k)=amin1( tmp2,qlzodt(k) ) !c !c (8) DEPOSITION/SUBLIMATION OF SNOW (Psdep/Pssub): Lin (31) !c includes consideration of ventilation effect !c !c abi=2*pi*(Si-1)/rho/(A"+B") !c tmpa=rvapor*xka(k)*tem(k)*tem(k) tmpb=xls*xls*rho(k)*qsiz(k)*diffwv(k) tmpc=tmpa*qsiz(k)*diffwv(k) abi=2.0*pi*(qvoqsiz(k)-1.0)*tmpc/(tmpa+tmpb) ! !c vf1s,vf2s=ventilation factors for snow !c vf1s=0.78,vf2s=0.31 in LIN ! tmp1=constc*sqrho(k)*olambdas(k)**dp5/visc(k) tmp2=abi*xnos*( vf1s*olambdas(k)*olambdas(k)+ & vf2s*schmidt(k)**0.33334*gamdp5o2*sqrt(tmp1) ) tmp3=odtb*( qvz(k)-qsiz(k) ) ! !bloss: BEGIN !the old implementation would give the wrong results if olambdas(k) ==0 !which can lead to positive pssub, i.e. tmp2=0 but tmp3> 0 ! if( tmp2 .le. 0.0) then if( tmp3 .le. 0.0) then tmp2=amax1( tmp2,tmp3) pssub(k)=amin1(0.,amax1( tmp2,-qszodt(k) )) !bloss: END psdep(k)=0.0 else psdep(k)=amin1( tmp2,tmp3 ) pssub(k)=0.0 end if !0915 psdep(k)=amax1(0.0,tmp2) !0915 pssub(k)=amin1(0.0,tmp2) !0915 pssub(k)=amax1( pssub(k),-qszodt(k) ) ! if(qrz(k) .le. 0.0) go to 1200 ! ! Compute processes (9) and (10) only when qsz > 0.0 and qrz > 0.0 ! !c !c (9) ACCRETION OF SNOW BY RAIN (Pracs): Lin (27) !c esr=1.0 tmpa=olambdar(k)*olambdar(k) tmpb=olambdas(k)*olambdas(k) tmpc=olambdar(k)*olambdas(k) tmp1=pi*pi*esr*xnor*xnos*abs( vtr(k)-vts(k) )*orho(k) tmp2=tmpb*tmpb*olambdar(k)*(5.0*tmpb+2.0*tmpc+0.5*tmpa) tmp3=tmp1*rhosnow*tmp2 pracs(k)=amin1( tmp3,qszodt(k) ) !c !c (10) ACCRETION OF RAIN BY SNOW (Psacr): Lin (28) !c tmp3=tmpa*tmpa*olambdas(k)*(5.0*tmpa+2.0*tmpc+0.5*tmpb) tmp4=tmp1*rhowater*tmp3 psacr(k)=amin1( tmp4,qrzodt(k) ) ! 1200 continue ! else ! !*********************************************************************** !********* snow production processes for T > 0 C ********** !*********************************************************************** ! if (qsz(k) .le. 0.0) go to 1400 !c !c (1) CLOUD WATER ACCRETION BY SNOW (Psacw): Lin (24) !c esw=1.0 tmp1=esw*pio4*xnos*constc*gamdp3*sqrho(k)* & olambdas(k)**dp3 psacw(k)=qlzodt(k)*( 1.0-exp(-tmp1*dtb) ) !0915 tmp1=pio4*xnos*constc*gamdp3*sqrho(k)* & !0915 olambdas(k)**dp3 !0915 tmp2=qlz(k)*esw*tmp1 !0915 psacw(k)=amin1( tmp2,qlzodt(k) ) !c !c (2) ACCRETION OF RAIN BY SNOW (Psacr): Lin (28) !c esr=1.0 tmpa=olambdar(k)*olambdar(k) tmpb=olambdas(k)*olambdas(k) tmpc=olambdar(k)*olambdas(k) tmp1=pi*pi*esr*xnor*xnos*abs( vtr(k)-vts(k) )*orho(k) tmp2=tmpa*tmpa*olambdas(k)*(5.0*tmpa+2.0*tmpc+0.5*tmpb) tmp3=tmp1*rhowater*tmp2 psacr(k)=amin1( tmp3,qrzodt(k) ) !c !c (3) MELTING OF SNOW (Psmlt): Lin (32) !c Psmlt is negative value ! delrs=rs0(k)-qvz(k) term1=2.0*pi*orho(k)*( xlv*diffwv(k)*rho(k)*delrs- & xka(k)*temcc(k) ) tmp1=constc*sqrho(k)*olambdas(k)**dp5/visc(k) tmp2=xnos*( vf1s*olambdas(k)*olambdas(k)+ & vf2s*schmidt(k)**0.33334*gamdp5o2*sqrt(tmp1) ) tmp3=term1*oxlf*tmp2-cwoxlf*temcc(k)*( psacw(k)+psacr(k) ) tmp4=amin1(0.0,tmp3) psmlt(k)=amax1( tmp4,-qszodt(k) ) !c !c (4) EVAPORATION OF MELTING SNOW (Psmltevp): HR (A27) !c but use Lin et al. coefficience !c Psmltevp is a negative value !c tmpa=rvapor*xka(k)*tem(k)*tem(k) tmpb=xlv*xlv*rho(k)*qswz(k)*diffwv(k) tmpc=tmpa*qswz(k)*diffwv(k) tmpd=amin1( 0.0,(qvoqswz(k)-0.90)*qswz(k)*odtb ) ! abr=2.0*pi*(qvoqswz(k)-1.0)*tmpc/(tmpa+tmpb) abr=2.0*pi*(qvoqswz(k)-0.90)*tmpc/(tmpa+tmpb) ! !**** allow evaporation to occur when RH less than 90% !**** here not using 100% because the evaporation cooling !**** of temperature is not taking into account yet; hence, !**** the qsw value is a little bit larger. This will avoid !**** evaporation can generate cloud. ! !c vf1s,vf2s=ventilation factors for snow !c vf1s=0.78,vf2s=0.31 in LIN ! tmp1=constc*sqrho(k)*olambdas(k)**dp5/visc(k) tmp2=abr*xnos*( vf1s*olambdas(k)*olambdas(k)+ & vf2s*schmidt(k)**0.33334*gamdp5o2*sqrt(tmp1) ) tmp3=amin1(0.0,tmp2) tmp3=amax1( tmp3,tmpd ) psmltevp(k)=amax1( tmp3,-qszodt(k) ) 1400 continue ! end if !*********************************************************************** !********* rain production processes ********** !*********************************************************************** ! !c !c (1) AUTOCONVERSION OF RAIN (Praut): RH !sg: begin if(flag_qndrop)then if( qndropz(k) >= 1. ) then ! Liu et al. autoconversion scheme rhocgs=rho(k)*1.e-3 liqconc=rhocgs*qlz(k) ! (kg/kg) to (g/cm3) capn=1.0e-3*rhocgs*qndropz(k) ! (#/kg) to (#/cm3) ! rate function if(liqconc.gt.1.e-10)then p0=(kappa*beta/capn)*(liqconc*liqconc*liqconc) xc=9.7d-17*capn*sqrt(capn)/(liqconc*liqconc) ! Calculate autoconversion rate (g/g/s) if(xc.lt.10.)then praut(k)=(p0/rhocgs) * ( 0.5d0*(xc*xc+2*xc+2.0d0)* & (1.0d0+xc)*exp(-2.0d0*xc) ) endif endif endif else !sg: end !c araut=afa*rho !c afa=0.001 Rate coefficient for autoconvergence !c !c araut=1.0e-3 !c araut=0.001 !testing ! tmp1=amax1( 0.0,araut*(qlz(k)-ql0) ) ! praut(k)=amin1( tmp1,qlzodt(k) ) tmp1=odtb*(qlz(k)-ql0)*( 1.0-exp(-araut*dtb) ) praut(k)=amax1( 0.0,tmp1 ) endif !sg !c !c (2) ACCRETION OF CLOUD WATER BY RAIN (Pracw): Lin (51) !c erw=1.0 ! tmp1=qlz(k)*pio4*erw*xnor*consta*sqrho(k) ! tmp2=tmp1*gambp3*olambdar(k)**bp3 ! pracw(k)=amin1( tmp2,qlzodt(k) ) tmp1=pio4*erw*xnor*consta*sqrho(k)* & gambp3*olambdar(k)**bp3 pracw(k)=qlzodt(k)*( 1.0-exp(-tmp1*dtb) ) !c !c (3) EVAPORATION OF RAIN (Prevp): Lin (52) !c Prevp is negative value !c !c Sw=qvoqsw : saturation ratio !c tmpa=rvapor*xka(k)*tem(k)*tem(k) tmpb=xlv*xlv*rho(k)*qswz(k)*diffwv(k) tmpc=tmpa*qswz(k)*diffwv(k) !bloss: BEGIN ! tmpd=amin1(0.0,(qvoqswz(k)-0.90)*qswz(k)*odtb) ! set max allowed evaporation to 90% of the amount that ! would induce saturation wrt liquid in a single step. tmpd = qswz(k)*xlv/(rvapor*tem(k)**2) ! d(qsat_liq)/dT tmpd = min( 0., 0.9*odtb*(qvz(k) + qlz(k) - qswz(k)) & / (1. + xlvocp * tmpd) ) abr=2.0*pi*(qvoqswz(k)-1.0)*tmpc/(tmpa+tmpb) ! abr=2.0*pi*(qvoqswz(k)-0.90)*tmpc/(tmpa+tmpb) !bloss: END ! !c vf1r,vf2r=ventilation factors for rain !c vf1r=0.78,vf2r=0.31 in RH, LIN and MM5 ! vf1r=0.78 vf2r=0.31 tmp1=consta*sqrho(k)*olambdar(k)**bp5/visc(k) tmp2=abr*xnor*( vf1r*olambdar(k)*olambdar(k)+ & vf2r*schmidt(k)**0.33334*gambp5o2*sqrt(tmp1) ) tmp3=amin1( 0.0,tmp2 ) tmp3=amax1( tmp3,tmpd ) prevp(k)=amax1( tmp3,-qrzodt(k) ) ! ! if(iout .gt. 0) write(20,*)'tmp1,tmp2,tmp3=',tmp1,tmp2,tmp3 ! if(iout .gt. 0) write(20,*)'qlz,qiz,qrz=',qlz(k),qiz(k),qrz(k) ! if(iout .gt. 0) write(20,*)'tem,qsz,qvz=',tem(k),qsz(k),qvz(k) ! if (gindex .eq. 0.) goto 900 ! if (tem(k) .lt. 273.15) then ! ! !-- graupel !*********************************************************************** !********* graupel production processes for T < 0 C ********** !*********************************************************************** !c !c (1) AUTOCONVERSION OF SNOW TO FORM GRAUPEL (Pgaut): Lin (37) !c pgaut=alpha2*(qsz-qs0) !c qs0=6.0E-4 !c alpha2=1.0e-3*exp(0.09*temcc(k)) Lin (38) ! alpha2=1.0e-3*exp(0.09*temcc(k)) ! ! testing ! tmp1=alpha2*(qsz(k)-qs0) ! tmp1=amax1(0.0,tmp1) ! pgaut(k)=amin1( tmp1,qszodt(k) ) tmp1=odtb*(qsz(k)-qs0)*(1.0-exp(-alpha2*dtb)) pgaut(k)=amax1( 0.0,tmp1 ) !c !c (2) FREEZING OF RAIN TO FORM GRAUPEL (Pgfr): Lin (45) !c positive value !c Constant in Bigg freezing Aplume=Ap=0.66 /k !c Constant in raindrop freezing equ. Bplume=Bp=100./m/m/m/s ! if (qrz(k) .gt. 1.e-8 ) then Bp=100. Ap=0.66 tmp1=olambdar(k)*olambdar(k)*olambdar(k) tmp2=20.*pi*pi*Bp*xnor*rhowater*orho(k)* & (exp(-Ap*temcc(k))-1.0)*tmp1*tmp1*olambdar(k) Pgfr(k)=amin1( tmp2,qrzodt(k) ) else Pgfr(k)=0 endif !c !c if (qgz(k) = 0.0) skip the other step below about graupel !c if (qgz(k) .eq. 0.0) goto 4000 !c !c Comparing Pgwet(wet process) and Pdry(dry process), !c we will pick up the small one. !c !c --------------- !c | dry processes | !c --------------- !c !c (3) ACCRETION OF CLOUD WATER BY GRAUPEL (Pgacw): Lin (40) !c egw=1.0 !c Cdrag=0.6 drag coefficients for hairstone !c constg=sqrt(4.*grav*rhograul*0.33334*orho(k)/Cdrag) !c egw=1.0 constg=sqrt(4.*grav*rhograul*0.33334*orho(k)*oCdrag) tmp1=pio4*xnog*gam3pt5*constg*olambdag(k)**3.5 tmp2=qlz(k)*egw*tmp1 Pgacw(k)=amin1( tmp2,qlzodt(k) ) !c !c (4) ACCRETION OF ICE CRYSTAL BY GRAUPEL (Pgaci): Lin (41) !c egi=1. for wet growth !c egi=0.1 for dry growth !c egi=0.1 tmp2=qiz(k)*egi*tmp1 pgaci(k)=amin1( tmp2,qizodt(k) ) !c !c (5) ACCRETION OF SNOW BY GRAUPEL (Pgacs) : Lin (29) !c Compute processes (6) only when qsz > 0.0 and qgz > 0.0 !c egs=exp(0.09*temcc(k)) tmpa=olambdas(k)*olambdas(k) tmpb=olambdag(k)*olambdag(k) tmpc=olambdas(k)*olambdag(k) tmp1=pi*pi*xnos*xnog*abs( vts(k)-vtg(k) )*orho(k) tmp2=tmpa*tmpa*olambdag(k)*(5.0*tmpa+2.0*tmpc+0.5*tmpb) tmp3=tmp1*egs*rhosnow*tmp2 Pgacs(k)=amin1( tmp3,qszodt(k) ) !c !c (6) ACCRETION OF RAIN BY GRAUPEL (Pgacr): Lin (42) !c Compute processes (6) only when qrz > 0.0 and qgz > 0.0 !c egr=1. !c egr=1. tmpa=olambdar(k)*olambdar(k) tmpb=olambdag(k)*olambdag(k) tmpc=olambdar(k)*olambdag(k) tmp1=pi*pi*xnor*xnog*abs( vtr(k)-vtg(k) )*orho(k) tmp2=tmpa*tmpa*olambdag(k)*(5.0*tmpa+2.0*tmpc+0.5*tmpb) tmp3=tmp1*egr*rhowater*tmp2 pgacr(k)=amin1( tmp3,qrzodt(k) ) !c !c (7) Calculate total dry process effect Pdry(k) !c Pdry(k)=Pgacw(k)+pgaci(k)+Pgacs(k)+pgacr(k) !c --------------- !c | wet processes | !c --------------- !c !c (3) ACCRETION OF ICE CRYSTAL BY GRAUPEL (Pgacip): Lin (41) !c egi=1. for wet growth !c egi=0.1 for dry growth !c tmp2=10.*pgaci(k) pgacip(k)=amin1( tmp2,qizodt(k) ) !c !c (4) ACCRETION OF SNOW BY GRAUPEL ((Pgacsp) : Lin (29) !c Compute processes (6) only when qsz > 0.0 and qgz > 0.0 !c egs=exp(0.09*(tem(k)-273.15)) when T < 273.15 k !c tmp3=Pgacs(k)*1.0/egs Pgacsp(k)=amin1( tmp3,qszodt(k) ) !c !c (5) WET GROWTH OF GRAUPEL (Pgwet) : Lin (43) !c may involve Pgacs or Pgaci and !c must include PPgacw or Pgacr, or both. !c ( The amount of Pgacw which is not able !c to freeze is shed to rain. ) IF(temcc(k).gt.-40.)THEN term0=constg*olambdag(k)**5.5/visc(k) !c !c vf1s,vf2s=ventilation factors for graupel !c vf1s=0.78,vf2s=0.31 in LIN !c Cdrag=0.6 drag coefficient for hairstone !c constg2=vf1s*olambdag(k)*olambdag(k)+ !c vf2s*schmidt(k)**0.33334*gam2pt75*sqrt(term0) delrs=rs0(k)-qvz(k) tmp0=1./(xlf+cw*temcc(k)) tmp1=2.*pi*xnog*(rho(k)*xlv*diffwv(k)*delrs-xka(k)* & temcc(k))*orho(k)*tmp0 constg2=vf1s*olambdag(k)*olambdag(k)+ & vf2s*schmidt(k)**0.33334*gam2pt75*sqrt(term0) tmp3=tmp1*constg2+(Pgacip(k)+Pgacsp(k))* & (1-Ci*temcc(k)*tmp0) tmp3=amax1(0.0,tmp3) Pgwet(k)=amin1(tmp3,qlzodt(k)+qszodt(k)+qizodt(k) ) !bloss !c !c Comparing Pgwet(wet process) and Pdry(dry process), !c we will apply the small one. !c if dry processes then delta4=1.0 !c if wet processes then delta4=0.0 ! if ( Pdry(k) .lt. Pgwet(k) ) then delta4=1.0 else delta4=0.0 endif ELSE delta4=1.0 ENDIF !c !c !c (6) Pgacrp(k)=Pgwet(k)-Pgacw(k)-Pgacip(k)-Pgacsp(k) !c if Pgacrp(k) > 0. then some of the rain is frozen to hail !c if Pgacrp(k) < 0. then some of the cloud water collected !c by the hail is unable to freeze and is !c shed as rain. !c Pgacrp(k)=Pgwet(k)-Pgacw(k)-Pgacip(k)-Pgacsp(k) !c !c (8) DEPOSITION/SUBLIMATION OF GRAUPEL (Pgdep/Pgsub): Lin (46) !c includes ventilation effect !c constg=sqrt(4.*grav*rhograul*0.33334*orho(k)/Cdrag) !c constg2=vf1s*olambdag(k)*olambdag(k)+ !c vf2s*schmidt(k)**0.33334*gam2pt75*constg !c !c abg=2*pi*(Si-1)/rho/(A"+B") !c tmpa=rvapor*xka(k)*tem(k)*tem(k) tmpb=xls*xls*rho(k)*qsiz(k)*diffwv(k) tmpc=tmpa*qsiz(k)*diffwv(k) abg=2.0*pi*(qvoqsiz(k)-1.0)*tmpc/(tmpa+tmpb) !c !c vf1s,vf2s=ventilation factors for graupel !c vf1s=0.78,vf2s=0.31 in LIN !c Cdrag=0.6 drag coefficient for hairstone !c term0=constg*olambdag(k)**5.5/visc(k) constg2=vf1s*olambdag(k)*olambdag(k)+ & vf2s*schmidt(k)**0.33334*gam2pt75*sqrt(term0) tmp2=abg*xnog*constg2 pgdep(k)=amax1(0.0,tmp2) pgsub(k)=amin1(0.0,tmp2) pgsub(k)=amax1( pgsub(k),-qgzodt(k) ) 4000 continue else ! !*********************************************************************** !********* graupel production processes for T > 0 C ********** !*********************************************************************** ! !c !c (1) ACCRETION OF CLOUD WATER BY GRAUPEL (Pgacw): Lin (40) !c egw=1.0 !c Cdrag=0.6 drag coefficients for hairstone !c constg=sqrt(4.*grav*rhograul*0.33334*orho(k)/Cdrag) egw=1.0 constg=sqrt(4.*grav*rhograul*0.33334*orho(k)*oCdrag) tmp1=pio4*xnog*gam3pt5*constg*olambdag(k)**3.5 tmp2=qlz(k)*egw*tmp1 Pgacw(k)=amin1( tmp2,qlzodt(k) ) !c !c (2) ACCRETION OF RAIN BY GRAUPEL (Pgacr): Lin (42) !c Compute processes (5) only when qrz > 0.0 and qgz > 0.0 !c egr=1. !c egr=1. tmpa=olambdar(k)*olambdar(k) tmpb=olambdag(k)*olambdag(k) tmpc=olambdar(k)*olambdag(k) tmp1=pi*pi*xnor*xnog*abs( vtr(k)-vtg(k) )*orho(k) tmp2=tmpa*tmpa*olambdag(k)*(5.0*tmpa+2.0*tmpc+0.5*tmpb) tmp3=tmp1*egr*rhowater*tmp2 pgacr(k)=amin1( tmp3,qrzodt(k) ) !c !c (3) GRAUPEL MELTING TO FORM RAIN (Pgmlt): Lin (47) !c Pgmlt is negative value !c constg=sqrt(4.*grav*rhograul*0.33334*orho(k)/Cdrag) !c constg2=vf1s*olambdag(k)*olambdag(k)+ !c vf2s*schmidt(k)**0.33334*gam2pt75*constg !c Cdrag=0.6 drag coefficients for hairstone ! delrs=rs0(k)-qvz(k) term1=2.0*pi*orho(k)*( xlv*diffwv(k)*rho(k)*delrs- & xka(k)*temcc(k) ) term0=sqrt(4.*grav*rhograul*0.33334*orho(k)*ocdrag) & *olambdag(k)**5.5/visc(k) constg2=vf1s*olambdag(k)*olambdag(k)+ & vf2s*schmidt(k)**0.33334*gam2pt75*sqrt(term0) tmp2=xnog*constg2 tmp3=term1*oxlf*tmp2-cwoxlf*temcc(k)*( pgacw(k)+pgacr(k) ) tmp4=amin1(0.0,tmp3) pgmlt(k)=amax1( tmp4,-qgzodt(k) ) !c !c (4) EVAPORATION OF MELTING GRAUPEL (Pgmltevp) : HR (A19) !c but use Lin et al. coefficience !c Pgmltevp is a negative value !c abg=2.0*pi*(qvoqsiz(k)-1.0)*tmpc/(tmpa+tmpb) !c tmpa=rvapor*xka(k)*tem(k)*tem(k) tmpb=xlv*xlv*rho(k)*qswz(k)*diffwv(k) tmpc=tmpa*qswz(k)*diffwv(k) tmpd=amin1( 0.0,(qvoqswz(k)-0.90)*qswz(k)*odtb ) !c !c abg=2*pi*(Si-1)/rho/(A"+B") !c abg=2.0*pi*(qvoqswz(k)-0.90)*tmpc/(tmpa+tmpb) ! !**** allow evaporation to occur when RH less than 90% !**** here not using 100% because the evaporation cooling !**** of temperature is not taking into account yet; hence, !**** the qgw value is a little bit larger. This will avoid !**** evaporation can generate cloud. ! !c vf1s,vf2s=ventilation factors for snow !c vf1s=0.78,vf2s=0.31 in LIN !c constg=sqrt(4.*grav*rhograul*0.33334*orho(k)/Cdrag) !c constg2=vf1s*olambdag(k)*olambdag(k)+ !c vf2s*schmidt(k)**0.33334*gam2pt75*constg ! tmp2=abg*xnog*constg2 tmp3=amin1(0.0,tmp2) tmp3=amax1( tmp3,tmpd ) pgmltevp(k)=amax1( tmp3,-qgzodt(k) ) !c !c (5) ACCRETION OF SNOW BY GRAUPEL (Pgacs) : Lin (29) !c Compute processes (3) only when qsz > 0.0 and qgz > 0.0 !c egs=1.0 !c egs=1. tmpa=olambdas(k)*olambdas(k) tmpb=olambdag(k)*olambdag(k) tmpc=olambdas(k)*olambdag(k) tmp1=pi*pi*xnos*xnog*abs( vts(k)-vtg(k) )*orho(k) tmp2=tmpa*tmpa*olambdag(k)*(5.0*tmpa+2.0*tmpc+0.5*tmpb) tmp3=tmp1*egs*rhosnow*tmp2 Pgacs(k)=amin1( tmp3,qszodt(k) ) endif ! 900 continue !cc !c !c********************************************************************** !c***** combine all processes together and avoid negative ***** !c***** water substances !*********************************************************************** !c if ( temcc(k) .lt. 0.0) then !,delta4,1.-delta4 !c !c gdelta4=gindex*delta4 !c g1sdelt4=gindex*(1.-delta4) !c gdelta4=gindex*delta4 g1sdelt4=gindex*(1.-delta4) !c !c combined water vapor depletions !c !cc graupel tmp=psdep(k)+pgdep(k)*gindex if ( tmp .gt. qvzodt(k) ) then factor=qvzodt(k)/tmp psdep(k)=psdep(k)*factor pgdep(k)=pgdep(k)*factor*gindex end if !c !c combined cloud water depletions !c tmp=praut(k)+psacw(k)+psfw(k)+pracw(k)+gindex*Pgacw(k) if ( tmp .gt. qlzodt(k) ) then factor=qlzodt(k)/tmp praut(k)=praut(k)*factor psacw(k)=psacw(k)*factor psfw(k)=psfw(k)*factor pracw(k)=pracw(k)*factor !cc graupel Pgacw(k)=Pgacw(k)*factor*gindex end if !c !c combined cloud ice depletions !c tmp=psaut(k)+psaci(k)+praci(k)+psfi(k)+Pgaci(k)*gdelta4 & +Pgacip(k)*g1sdelt4 if (tmp .gt. qizodt(k) ) then factor=qizodt(k)/tmp psaut(k)=psaut(k)*factor psaci(k)=psaci(k)*factor praci(k)=praci(k)*factor psfi(k)=psfi(k)*factor !cc graupel Pgaci(k)=Pgaci(k)*factor*gdelta4 Pgacip(k)=Pgacip(k)*factor*g1sdelt4 endif !c !c combined all rain processes !c tmp_r=piacr(k)+psacr(k)-prevp(k)-praut(k)-pracw(k) & +Pgfr(k)*gindex+Pgacr(k)*gdelta4 & +Pgacrp(k)*g1sdelt4 if (tmp_r .gt. qrzodt(k) ) then factor=qrzodt(k)/tmp_r piacr(k)=piacr(k)*factor psacr(k)=psacr(k)*factor prevp(k)=prevp(k)*factor !cc graupel Pgfr(k)=Pgfr(k)*factor*gindex Pgacr(k)=Pgacr(k)*factor*gdelta4 Pgacrp(k)=Pgacrp(k)*factor*g1sdelt4 endif !c !c if qrz < 1.0E-4 and qsz < 1.0E-4 then delta2=1. !c (all Pracs and Psacr become to snow) !c if qrz >= 1.0E-4 or qsz >= 1.0E-4 then delta2=0. !c (all Pracs and Psacr become to graupel) !c if (qrz(k) .lt. 1.0E-4 .and. qsz(k) .lt. 1.0E-4) then delta2=1.0 else delta2=0.0 endif ! !cc graupel !c !c if qrz(k) < 1.0e-4 then delta3=1. means praci(k) --> qs !c piacr(k) --> qs !c if qrz(k) > 1.0e-4 then delta3=0. means praci(k) --> qg !c piacr(k) --> qg : Lin (20) if (qrz(k) .lt. 1.0e-4) then delta3=1.0 else delta3=0.0 endif ! !c !c if gindex = 0.(no graupel) then delta2=1.0 !c delta3=1.0 !c if (gindex .eq. 0.) then delta2=1.0 delta3=1.0 endif ! !c !c combined all snow processes !c tmp_s=-pssub(k)-(psaut(k)+psaci(k)+psacw(k)+psfw(k)+ & psfi(k)+praci(k)*delta3+piacr(k)*delta3+ & psdep(k))+Pgaut(k)*gindex+Pgacs(k)*gdelta4+ & Pgacsp(k)*g1sdelt4+Pracs(k)*(1.-delta2)- & Psacr(k)*delta2 if ( tmp_s .gt. qszodt(k) ) then factor=qszodt(k)/tmp_s pssub(k)=pssub(k)*factor Pracs(k)=Pracs(k)*factor !cc graupel Pgaut(k)=Pgaut(k)*factor*gindex Pgacs(k)=Pgacs(k)*factor*gdelta4 Pgacsp(k)=Pgacsp(k)*factor*g1sdelt4 endif !cc graupel ! ! if (gindex .eq. 0.) goto 998 !c !c combined all graupel processes !c if(delta4.lt.0.5) then !Re-define pgwet to account for limiting of pgacrp, ! pgacip, pgacw and pgacsp above pgwet(k) = pgacrp(k) + pgacw(k) + pgacip(k) + pgacsp(k) end if tmp_g=-pgaut(k)-pgfr(k)-Pgacw(k)*delta4-Pgaci(k)*delta4 & -Pgacr(k)*delta4-Pgacs(k)*delta4 & -pgwet(k)*(1.-delta4)-pgsub(k)-pgdep(k) & -psacr(k)*(1-delta2)-Pracs(k)*(1-delta2) & -praci(k)*(1-delta3)-piacr(k)*(1-delta3) if (tmp_g .gt. qgzodt(k)) then factor=qgzodt(k)/tmp_g pgsub(k)=pgsub(k)*factor endif 998 continue !c !c calculate new water substances, thetae, tem, and qvsbar !c !cc graupel pvapor(k)=-pssub(k)-psdep(k)-prevp(k)-pgsub(k)*gindex & -pgdep(k)*gindex qvz(k)=amax1( qvmin,qvz(k)+dtb*pvapor(k) ) pclw(k)=-praut(k)-pracw(k)-psacw(k)-psfw(k)-pgacw(k)*gindex if(flag_qndrop)then if( qlz(k) > 1e-20 ) & qndropz(k)=amax1( 0.0,qndropz(k)+dtb*pclw(k)*qndropz(k)/qlz(k) ) !sg endif qlz(k)=amax1( 0.0,qlz(k)+dtb*pclw(k) ) pcli(k)=-psaut(k)-psfi(k)-psaci(k)-praci(k)-pgaci(k)*gdelta4 & -Pgacip(k)*g1sdelt4 qiz(k)=amax1( 0.0,qiz(k)+dtb*pcli(k) ) tmp_r=piacr(k)+psacr(k)-prevp(k)-praut(k)-pracw(k) & +Pgfr(k)*gindex+Pgacr(k)*gdelta4 & +Pgacrp(k)*g1sdelt4 232 format(i2,1x,6(f9.3,1x)) prain(k)=-tmp_r qrz(k)=amax1( 0.0,qrz(k)+dtb*prain(k) ) tmp_s=-pssub(k)-(psaut(k)+psaci(k)+psacw(k)+psfw(k)+ & psfi(k)+praci(k)*delta3+piacr(k)*delta3+ & psdep(k))+Pgaut(k)*gindex+Pgacs(k)*gdelta4+ & Pgacsp(k)*g1sdelt4+Pracs(k)*(1.-delta2)- & Psacr(k)*delta2 psnow(k)=-tmp_s qsz(k)=amax1( 0.0,qsz(k)+dtb*psnow(k) ) qschg(k)=qschg(k)+psnow(k) qschg(k)=psnow(k) !cc graupel tmp_g=-pgaut(k)-pgfr(k)-Pgacw(k)*delta4-Pgaci(k)*delta4 & -Pgacr(k)*delta4-Pgacs(k)*delta4 & -pgwet(k)*(1.-delta4)-pgsub(k)-pgdep(k) & -psacr(k)*(1-delta2)-Pracs(k)*(1-delta2) & -praci(k)*(1-delta3)-piacr(k)*(1-delta3) 252 format(i2,1x,6(f12.9,1x)) 262 format(i2,1x,7(f12.9,1x)) pgraupel(k)=-tmp_g pgraupel(k)=pgraupel(k)*gindex qgz(k)=amax1( 0.0,qgz(k)+dtb*pgraupel(k)) ! qgchg(k)=qgchg(k)+pgraupel(k) qgchg(k)=pgraupel(k) qgz(k)=qgz(k)*gindex tmp=ocp/tothz(k)*xLf*(qschg(k)+qgchg(k)) theiz(k)=theiz(k)+dtb*tmp thz(k)=theiz(k)-(xLvocp*qvz(k)-xLfocp*qiz(k))/tothz(k) tem(k)=thz(k)*tothz(k) temcc(k)=tem(k)-273.15 !bloss: Moves computation of saturation mixing ratio below ! if( temcc(k) .lt. -40.0 ) qswz(k)=qsiz(k) ! qlpqi=qlz(k)+qiz(k) ! if ( qlpqi .eq. 0.0 ) then ! qvsbar(k)=qsiz(k) ! else ! qvsbar(k)=( qiz(k)*qsiz(k)+qlz(k)*qswz(k) )/qlpqi ! endif ! else !c !c combined cloud water depletions !c tmp=praut(k)+psacw(k)+pracw(k)+pgacw(k)*gindex if ( tmp .gt. qlzodt(k) ) then factor=qlzodt(k)/tmp praut(k)=praut(k)*factor psacw(k)=psacw(k)*factor pracw(k)=pracw(k)*factor !cc graupel pgacw(k)=pgacw(k)*factor*gindex end if !c !c combined all snow processes !c tmp_s=-(psmlt(k)+psmltevp(k))+Pgacs(k)*gindex if (tmp_s .gt. qszodt(k) ) then factor=qszodt(k)/tmp_s psmlt(k)=psmlt(k)*factor psmltevp(k)=psmltevp(k)*factor !cc graupel Pgacs(k)=Pgacs(k)*factor*gindex endif !c !c !cc graupel !c ! if (gindex .eq. 0.) goto 997 !c !c combined all graupel processes !c tmp_g=-pgmlt(k)-pgacs(k)-pgmltevp(k) if (tmp_g .gt. qgzodt(k)) then factor=qgzodt(k)/tmp_g pgmltevp(k)=pgmltevp(k)*factor pgmlt(k)=pgmlt(k)*factor endif !c 997 continue !c !c combined all rain processes !c tmp_r=-prevp(k)-(praut(k)+pracw(k)+psacw(k)-psmlt(k)) & +pgmlt(k)*gindex-pgacw(k)*gindex if (tmp_r .gt. qrzodt(k) ) then factor=qrzodt(k)/tmp_r prevp(k)=prevp(k)*factor endif !c !c !c calculate new water substances and thetae !c pvapor(k)=-psmltevp(k)-prevp(k)-pgmltevp(k)*gindex qvz(k)=amax1( qvmin,qvz(k)+dtb*pvapor(k)) pclw(k)=-praut(k)-pracw(k)-psacw(k)-pgacw(k)*gindex if(flag_qndrop)then if( qlz(k) > 1e-20 ) & qndropz(k)=amax1( 0.0,qndropz(k)+dtb*pclw(k)*qndropz(k)/qlz(k) ) !sg endif qlz(k)=amax1( 0.0,qlz(k)+dtb*pclw(k) ) pcli(k)=0.0 qiz(k)=amax1( 0.0,qiz(k)+dtb*pcli(k) ) tmp_r=-prevp(k)-(praut(k)+pracw(k)+psacw(k)-psmlt(k)) & +pgmlt(k)*gindex-pgacw(k)*gindex 242 format(i2,1x,7(f9.6,1x)) prain(k)=-tmp_r tmpqrz=qrz(k) qrz(k)=amax1( 0.0,qrz(k)+dtb*prain(k) ) tmp_s=-(psmlt(k)+psmltevp(k))+Pgacs(k)*gindex psnow(k)=-tmp_s qsz(k)=amax1( 0.0,qsz(k)+dtb*psnow(k) ) ! qschg(k)=qschg(k)+psnow(k) qschg(k)=psnow(k) !cc graupel tmp_g=-pgmlt(k)-pgacs(k)-pgmltevp(k) ! write(*,272)k,pgmlt(k),pgacs(k),pgmltevp(k), 272 format(i2,1x,3(f12.9,1x)) pgraupel(k)=-tmp_g*gindex qgz(k)=amax1( 0.0,qgz(k)+dtb*pgraupel(k)) ! qgchg(k)=qgchg(k)+pgraupel(k) qgchg(k)=pgraupel(k) qgz(k)=qgz(k)*gindex ! tmp=ocp/tothz(k)*xLf*(qschg(k)+qgchg(k)) theiz(k)=theiz(k)+dtb*tmp thz(k)=theiz(k)-(xLvocp*qvz(k)-xLfocp*qiz(k))/tothz(k) tem(k)=thz(k)*tothz(k) temcc(k)=tem(k)-273.15 ! qswz(k)=episp0k*oprez(k)* & ! exp( svp2*temcc(k)/(tem(k)-svp3) ) !bloss: Moved computation of saturation mixing ratio below ! es=1000.*svp1*exp( svp2*temcc(k)/(tem(k)-svp3) ) ! qswz(k)=ep2*es/(prez(k)-es) ! qsiz(k)=qswz(k) ! qvsbar(k)=qswz(k) ! end if preclw(k)=pclw(k) !sg ! !*********************************************************************** !********** saturation adjustment ********** !*********************************************************************** !bloss: BEGIN ! ! compute saturation specific humidity at the temperature ! which would be experienced if all cloud liquid and ice ! were to become vapor. This will make for a consistent ! check of saturation. Previously, qvsbar was computed ! without accounting for the change in temperature due to ! evaporation/sublimation, and as a result would ! incorrectly identify some points as subsaturated. tmp_tem = tem(k) ! updated temperature from above. if(qlz(k)+qiz(k).gt. 0.) then ! account for latent heat if all qlz and qiz were converted to vapor tmp_tem = tem(k) - xlvocp*qlz(k) - (xlvocp+xlfocp)*qiz(k) end if ! compute temperature in celsius tmp_temcc = tmp_tem - 273.15 ! estimate lower bound on saturation vapor pressure ! (ice if <0C, liquid aboce) if (tmp_temcc .lt. 0.) then es=1000.*svp1*exp( 21.8745584*(tmp_tem-273.16)/(tmp_tem-7.66) ) !ice else es=1000.*svp1*exp( svp2*tmp_temcc/(tmp_tem-svp3) ) !liquid end if ! compute lower bound on saturation mixing ratio. qvsbar(k)=ep2*es/(prez(k)-es) ! !bloss: END ! ! allow supersaturation exits linearly from 0% at 500 mb to 50% ! above 300 mb ! 5.0e-5=1.0/(500mb-300mb) ! rsat=1.0+0.5*(50000.0-prez(k))*5.0e-5 rsat=amax1(1.0,rsat) rsat=amin1(1.5,rsat) rsat=1.0 if( qvz(k)+qlz(k)+qiz(k) .lt. rsat*qvsbar(k) ) then !c !c unsaturated !c qvz(k)=qvz(k)+qlz(k)+qiz(k) qlz(k)=0.0 qiz(k)=0.0 thz(k)=theiz(k)-(xLvocp*qvz(k)-xLfocp*qiz(k))/tothz(k) tem(k)=thz(k)*tothz(k) temcc(k)=tem(k)-273.15 go to 1800 ! else !c !c saturated !c ! pladj(k)=qlz(k) piadj(k)=qiz(k) ! CALL satadj(qvz, qlz, qiz, prez, theiz, thz, tothz, kts, kte, & k, xLvocp, xLfocp, episp0k, EP2,SVP1,SVP2,SVP3,SVPT0 ) ! pladj(k)=odtb*(qlz(k)-pladj(k)) piadj(k)=odtb*(qiz(k)-piadj(k)) ! pclw(k)=pclw(k)+pladj(k) pcli(k)=pcli(k)+piadj(k) pvapor(k)=pvapor(k)-( pladj(k)+piadj(k) ) ! thz(k)=theiz(k)-(xLvocp*qvz(k)-xLfocp*qiz(k))/tothz(k) tem(k)=thz(k)*tothz(k) temcc(k)=tem(k)-273.15 ! qswz(k)=episp0k*oprez(k)* & ! exp( svp2*temcc(k)/(tem(k)-svp3) ) es=1000.*svp1*exp( svp2*temcc(k)/(tem(k)-svp3) ) qswz(k)=ep2*es/(prez(k)-es) if (tem(k) .lt. 273.15 ) then ! qsiz(k)=episp0k*oprez(k)* & ! exp( 21.8745584*(tem(k)-273.16)/(tem(k)-7.66) ) es=1000.*svp1*exp( 21.8745584*(tem(k)-273.16)/(tem(k)-7.66) ) qsiz(k)=ep2*es/(prez(k)-es) if (temcc(k) .lt. -40.0) qswz(k)=qsiz(k) else qsiz(k)=qswz(k) endif qlpqi=qlz(k)+qiz(k) if ( qlpqi .eq. 0.0 ) then qvsbar(k)=qsiz(k) else qvsbar(k)=( qiz(k)*qsiz(k)+qlz(k)*qswz(k) )/qlpqi endif end if ! !*********************************************************************** !***** melting and freezing of cloud ice and cloud water ***** !*********************************************************************** qlpqi=qlz(k)+qiz(k) if(qlpqi .le. 0.0) go to 1800 ! !c !c (1) HOMOGENEOUS NUCLEATION WHEN T< -40 C (Pihom) !c if(temcc(k) .lt. -40.0) pihom(k)=qlz(k)*odtb !c !c (2) MELTING OF ICE CRYSTAL WHEN T> 0 C (Pimlt) !c if(temcc(k) .gt. 0.0) pimlt(k)=qiz(k)*odtb !c !c (3) PRODUCTION OF CLOUD ICE BY BERGERON PROCESS (Pidw): Hsie (p957) !c this process only considered when -31 C < T < 0 C !c if(temcc(k) .lt. 0.0 .and. temcc(k) .gt. -31.0) then !c! !c! parama1 and parama2 functions must be user supplied !c! a1=parama1( temcc(k) ) a2=parama2( temcc(k) ) !! change unit from cgs to mks a1=a1*0.001**(1.0-a2) xnin=xni0*exp(-bni*temcc(k)) pidw(k)=xnin*orho(k)*(a1*xmnin**a2) end if ! pcli(k)=pcli(k)+pihom(k)-pimlt(k)+pidw(k) pclw(k)=pclw(k)-pihom(k)+pimlt(k)-pidw(k) qlz(k)=amax1( 0.0,qlz(k)+dtb*(-pihom(k)+pimlt(k)-pidw(k)) ) qiz(k)=amax1( 0.0,qiz(k)+dtb*(pihom(k)-pimlt(k)+pidw(k)) ) ! CALL satadj(qvz, qlz, qiz, prez, theiz, thz, tothz, kts, kte, & k, xLvocp, xLfocp, episp0k ,EP2,SVP1,SVP2,SVP3,SVPT0) thz(k)=theiz(k)-(xLvocp*qvz(k)-xLfocp*qiz(k))/tothz(k) tem(k)=thz(k)*tothz(k) temcc(k)=tem(k)-273.15 ! qswz(k)=episp0k*oprez(k)* & ! exp( svp2*temcc(k)/(tem(k)-svp3) ) es=1000.*svp1*exp( svp2*temcc(k)/(tem(k)-svp3) ) qswz(k)=ep2*es/(prez(k)-es) if (tem(k) .lt. 273.15 ) then ! qsiz(k)=episp0k*oprez(k)* & ! exp( 21.8745584*(tem(k)-273.16)/(tem(k)-7.66) ) es=1000.*svp1*exp( 21.8745584*(tem(k)-273.16)/(tem(k)-7.66) ) qsiz(k)=ep2*es/(prez(k)-es) if (temcc(k) .lt. -40.0) qswz(k)=qsiz(k) else qsiz(k)=qswz(k) endif qlpqi=qlz(k)+qiz(k) if ( qlpqi .eq. 0.0 ) then qvsbar(k)=qsiz(k) else qvsbar(k)=( qiz(k)*qsiz(k)+qlz(k)*qswz(k) )/qlpqi endif 1800 continue ! !*********************************************************************** !********** integrate the productions of rain and snow ********** !*********************************************************************** !c 2000 continue !--------------------------------------------------------------------- ! !*********************************************************************** !****** Write terms in cloud physics to time series dataset ***** !*********************************************************************** ! ! open(unit=24,form='formatted',status='new', ! & file='cloud.dat') !9030 format(10e12.6) ! write(24,*)'tmp' ! write(24,9030) (tem(k),k=kts+1,kte) ! write(24,*)'qiz' ! write(24,9030) (qiz(k),k=kts+1,kte) ! write(24,*)'qsz' ! write(24,9030) (qsz(k),k=kts+1,kte) ! write(24,*)'qrz' ! write(24,9030) (qrz(k),k=kts+1,kte) ! write(24,*)'qgz' ! write(24,9030) (qgz(k),k=kts+1,kte) ! write(24,*)'qvoqsw' ! write(24,9030) (qvoqswz(k),k=kts+1,kte) ! write(24,*)'qvoqsi' ! write(24,9030) (qvoqsiz(k),k=kts+1,kte) ! write(24,*)'vtr' ! write(24,9030) (vtr(k),k=kts+1,kte) ! write(24,*)'vts' ! write(24,9030) (vts(k),k=kts+1,kte) ! write(24,*)'vtg' ! write(24,9030) (vtg(k),k=kts+1,kte) ! write(24,*)'pclw' ! write(24,9030) (pclw(k),k=kts+1,kte) ! write(24,*)'pvapor' ! write(24,9030) (pvapor(k),k=kts+1,kte) ! write(24,*)'pcli' ! write(24,9030) (pcli(k),k=kts+1,kte) ! write(24,*)'pimlt' ! write(24,9030) (pimlt(k),k=kts+1,kte) ! write(24,*)'pihom' ! write(24,9030) (pihom(k),k=kts+1,kte) ! write(24,*)'pidw' ! write(24,9030) (pidw(k),k=kts+1,kte) ! write(24,*)'prain' ! write(24,9030) (prain(k),k=kts+1,kte) ! write(24,*)'praut' ! write(24,9030) (praut(k),k=kts+1,kte) ! write(24,*)'pracw' ! write(24,9030) (pracw(k),k=kts+1,kte) ! write(24,*)'prevp' ! write(24,9030) (prevp(k),k=kts+1,kte) ! write(24,*)'psnow' ! write(24,9030) (psnow(k),k=kts+1,kte) ! write(24,*)'psaut' ! write(24,9030) (psaut(k),k=kts+1,kte) ! write(24,*)'psfw' ! write(24,9030) (psfw(k),k=kts+1,kte) ! write(24,*)'psfi' ! write(24,9030) (psfi(k),k=kts+1,kte) ! write(24,*)'praci' ! write(24,9030) (praci(k),k=kts+1,kte) ! write(24,*)'piacr' ! write(24,9030) (piacr(k),k=kts+1,kte) ! write(24,*)'psaci' ! write(24,9030) (psaci(k),k=kts+1,kte) ! write(24,*)'psacw' ! write(24,9030) (psacw(k),k=kts+1,kte) ! write(24,*)'psdep' ! write(24,9030) (psdep(k),k=kts+1,kte) ! write(24,*)'pssub' ! write(24,9030) (pssub(k),k=kts+1,kte) ! write(24,*)'pracs' ! write(24,9030) (pracs(k),k=kts+1,kte) ! write(24,*)'psacr' ! write(24,9030) (psacr(k),k=kts+1,kte) ! write(24,*)'psmlt' ! write(24,9030) (psmlt(k),k=kts+1,kte) ! write(24,*)'psmltevp' ! write(24,9030) (psmltevp(k),k=kts+1,kte) ! write(24,*)'pladj' ! write(24,9030) (pladj(k),k=kts+1,kte) ! write(24,*)'piadj' ! write(24,9030) (piadj(k),k=kts+1,kte) ! write(24,*)'pgraupel' ! write(24,9030) (pgraupel(k),k=kts+1,kte) ! write(24,*)'pgaut' ! write(24,9030) (pgaut(k),k=kts+1,kte) ! write(24,*)'pgfr' ! write(24,9030) (pgfr(k),k=kts+1,kte) ! write(24,*)'pgacw' ! write(24,9030) (pgacw(k),k=kts+1,kte) ! write(24,*)'pgaci' ! write(24,9030) (pgaci(k),k=kts+1,kte) ! write(24,*)'pgacr' ! write(24,9030) (pgacr(k),k=kts+1,kte) ! write(24,*)'pgacs' ! write(24,9030) (pgacs(k),k=kts+1,kte) ! write(24,*)'pgacip' ! write(24,9030) (pgacip(k),k=kts+1,kte) ! write(24,*)'pgacrP' ! write(24,9030) (pgacrP(k),k=kts+1,kte) ! write(24,*)'pgacsp' ! write(24,9030) (pgacsp(k),k=kts+1,kte) ! write(24,*)'pgwet' ! write(24,9030) (pgwet(k),k=kts+1,kte) ! write(24,*)'pdry' ! write(24,9030) (pdry(k),k=kts+1,kte) ! write(24,*)'pgsub' ! write(24,9030) (pgsub(k),k=kts+1,kte) ! write(24,*)'pgdep' ! write(24,9030) (pgdep(k),k=kts+1,kte) ! write(24,*)'pgmlt' ! write(24,9030) (pgmlt(k),k=kts+1,kte) ! write(24,*)'pgmltevp' ! write(24,9030) (pgmltevp(k),k=kts+1,kte) !**** below if qv < qvmin then qv=qvmin, ql=0.0, and qi=0.0 ! do k=kts+1,kte if ( qvz(k) .lt. qvmin ) then qlz(k)=0.0 qiz(k)=0.0 qvz(k)=amax1( qvmin,qvz(k)+qlz(k)+qiz(k) ) end if enddo ! END SUBROUTINE clphy1d !--------------------------------------------------------------------- ! SATURATED ADJUSTMENT !--------------------------------------------------------------------- SUBROUTINE satadj(qvz, qlz, qiz, prez, theiz, thz, tothz, & kts, kte, k, xLvocp, xLfocp, episp0k, EP2,SVP1,SVP2,SVP3,SVPT0) !--------------------------------------------------------------------- IMPLICIT NONE !--------------------------------------------------------------------- ! This program use Newton's method for finding saturated temperature ! and saturation mixing ratio. ! ! In this saturation adjustment scheme we assume ! (1) the saturation mixing ratio is the mass weighted average of ! saturation values over liquid water (qsw), and ice (qsi) ! following Lord et al., 1984 and Tao, 1989 ! ! (2) the percentage of cloud liquid and cloud ice will ! be fixed during the saturation calculation !--------------------------------------------------------------------- ! INTEGER, INTENT(IN ) :: kts, kte, k REAL, DIMENSION( kts:kte ), & INTENT(INOUT) :: qvz, qlz, qiz ! REAL, DIMENSION( kts:kte ), & INTENT(IN ) :: prez, theiz, tothz REAL, INTENT(IN ) :: xLvocp, xLfocp, episp0k REAL, INTENT(IN ) :: EP2,SVP1,SVP2,SVP3,SVPT0 ! LOCAL VARS INTEGER :: n REAL, DIMENSION( kts:kte ) :: thz, tem, temcc, qsiz, & qswz, qvsbar REAL :: qsat, qlpqi, ratql, t0, t1, tmp1, ratqi, tsat, absft, & denom1, denom2, dqvsbar, ftsat, dftsat, qpz, & gindex, es REAL :: tem_noliqice, qsat_noliqice !bloss ! !--------------------------------------------------------------------- thz(k)=theiz(k)-(xLvocp*qvz(k)-xLfocp*qiz(k))/tothz(k) tem(k)=tothz(k)*thz(k) !bloss: BEGIN tem_noliqice = tem(k) - xlvocp*qlz(k) - (xLvocp+xLfocp)*qiz(k) if (tem_noliqice .gt. 273.15) then es=1000.*svp1*exp( svp2*(tem_noliqice-svpt0)/(tem_noliqice-svp3) ) qsat_noliqice=ep2*es/(prez(k)-es) else qsat_noliqice=episp0k/prez(k)* & exp( 21.8745584*(tem_noliqice-273.15)/(tem_noliqice-7.66) ) end if !bloss: END qpz=qvz(k)+qlz(k)+qiz(k) if (qpz .lt. qsat_noliqice) then qvz(k)=qpz qiz(k)=0.0 qlz(k)=0.0 go to 400 end if qlpqi=qlz(k)+qiz(k) if( qlpqi .ge. 1.0e-5) then ratql=qlz(k)/qlpqi ratqi=qiz(k)/qlpqi else t0=273.15 ! t1=233.15 t1=248.15 tmp1=( t0-tem(k) )/(t0-t1) tmp1=amin1(1.0,tmp1) tmp1=amax1(0.0,tmp1) ratqi=tmp1 ratql=1.0-tmp1 end if ! ! !-- saturation mixing ratios over water and ice !-- at the outset we will follow Bolton 1980 MWR for !-- the water and Murray JAS 1967 for the ice ! !-- dqvsbar=d(qvsbar)/dT !-- ftsat=F(Tsat) !-- dftsat=d(F(T))/dT ! ! First guess of tsat tsat=tem(k) absft=1.0 ! do 200 n=1,20 denom1=1.0/(tsat-svp3) denom2=1.0/(tsat-7.66) ! qswz(k)=episp0k/prez(k)* & ! exp( svp2*denom1*(tsat-273.15) ) es=1000.*svp1*exp( svp2*denom1*(tsat-svpt0) ) qswz(k)=ep2*es/(prez(k)-es) if (tem(k) .lt. 273.15) then ! qsiz(k)=episp0k/prez(k)* & ! exp( 21.8745584*denom2*(tsat-273.15) ) es=1000.*svp1*exp( 21.8745584*denom2*(tsat-273.15) ) qsiz(k)=ep2*es/(prez(k)-es) if (tem(k) .lt. 233.15) qswz(k)=qsiz(k) else qsiz(k)=qswz(k) endif qvsbar(k)=ratql*qswz(k)+ratqi*qsiz(k) ! ! if( absft .lt. 0.01 .and. n .gt. 3 ) go to 300 if( absft .lt. 0.01 ) go to 300 ! dqvsbar=ratql*qswz(k)*svp2*243.5*denom1*denom1+ & ratqi*qsiz(k)*21.8745584*265.5*denom2*denom2 ftsat=tsat+(xlvocp+ratqi*xlfocp)*qvsbar(k)- & tothz(k)*theiz(k)-xlfocp*ratqi*(qvz(k)+qlz(k)+qiz(k)) dftsat=1.0+(xlvocp+ratqi*xlfocp)*dqvsbar tsat=tsat-ftsat/dftsat absft=abs(ftsat) 200 continue 9020 format(1x,'point can not converge, absft,n=',e12.5,i5) ! 300 continue if( qpz .gt. qvsbar(k) ) then qvz(k)=qvsbar(k) qiz(k)=ratqi*( qpz-qvz(k) ) qlz(k)=ratql*( qpz-qvz(k) ) else qvz(k)=qpz qiz(k)=0.0 qlz(k)=0.0 end if 400 continue END SUBROUTINE satadj !---------------------------------------------------------------- REAL FUNCTION parama1(temp) !---------------------------------------------------------------- IMPLICIT NONE !---------------------------------------------------------------- ! This program calculate the parameter for crystal growth rate ! in Bergeron process !---------------------------------------------------------------- REAL, INTENT (IN ) :: temp REAL, DIMENSION(32) :: a1 INTEGER :: i1, i1p1 REAL :: ratio data a1/0.100e-10,0.7939e-7,0.7841e-6,0.3369e-5,0.4336e-5, & 0.5285e-5,0.3728e-5,0.1852e-5,0.2991e-6,0.4248e-6, & 0.7434e-6,0.1812e-5,0.4394e-5,0.9145e-5,0.1725e-4, & 0.3348e-4,0.1725e-4,0.9175e-5,0.4412e-5,0.2252e-5, & 0.9115e-6,0.4876e-6,0.3473e-6,0.4758e-6,0.6306e-6, & 0.8573e-6,0.7868e-6,0.7192e-6,0.6513e-6,0.5956e-6, & 0.5333e-6,0.4834e-6/ i1=int(-temp)+1 i1p1=i1+1 ratio=-(temp)-float(i1-1) parama1=a1(i1)+ratio*( a1(i1p1)-a1(i1) ) END FUNCTION parama1 !---------------------------------------------------------------- REAL FUNCTION parama2(temp) !---------------------------------------------------------------- IMPLICIT NONE !---------------------------------------------------------------- ! This program calculate the parameter for crystal growth rate ! in Bergeron process !---------------------------------------------------------------- REAL, INTENT (IN ) :: temp REAL, DIMENSION(32) :: a2 INTEGER :: i1, i1p1 REAL :: ratio data a2/0.0100,0.4006,0.4831,0.5320,0.5307,0.5319,0.5249, & 0.4888,0.3849,0.4047,0.4318,0.4771,0.5183,0.5463, & 0.5651,0.5813,0.5655,0.5478,0.5203,0.4906,0.4447, & 0.4126,0.3960,0.4149,0.4320,0.4506,0.4483,0.4460, & 0.4433,0.4413,0.4382,0.4361/ i1=int(-temp)+1 i1p1=i1+1 ratio=-(temp)-float(i1-1) parama2=a2(i1)+ratio*( a2(i1p1)-a2(i1) ) END FUNCTION parama2 !---------------------------------------------------------------- REAL FUNCTION ggamma(X) !---------------------------------------------------------------- IMPLICIT NONE !---------------------------------------------------------------- REAL, INTENT(IN ) :: x REAL, DIMENSION(8) :: B INTEGER ::j, K1 REAL ::PF, G1TO2 ,TEMP DATA B/-.577191652,.988205891,-.897056937,.918206857, & -.756704078,.482199394,-.193527818,.035868343/ PF=1. TEMP=X DO 10 J=1,200 IF (TEMP .LE. 2) GO TO 20 TEMP=TEMP-1. 10 PF=PF*TEMP 100 FORMAT(//,5X,'module_mp_lin: INPUT TO GAMMA FUNCTION TOO LARGE, X=',E12.5) WRITE(wrf_err_message,100)X CALL wrf_error_fatal(wrf_err_message) 20 G1TO2=1. TEMP=TEMP - 1. DO 30 K1=1,8 30 G1TO2=G1TO2 + B(K1)*TEMP**K1 ggamma=PF*G1TO2 END FUNCTION ggamma !---------------------------------------------------------------- END MODULE module_mp_lin