MODULE module_mp_gsfcgce USE module_wrf_error USE module_mp_radar REAL, PRIVATE :: rd1, rd2, al, cp REAL, PRIVATE :: c38, c358, c610, c149, & c879, c172, c409, c76, & c218, c580, c141 REAL, PRIVATE :: ag, bg, as, bs, & aw, bw, bgh, bgq, & bsh, bsq, bwh, bwq REAL, PRIVATE :: tnw, tns, tng, & roqs, roqg, roqr REAL, PRIVATE :: zrc, zgc, zsc, & vrc, vgc, vsc REAL, PRIVATE :: alv, alf, als, t0, t00, & avc, afc, asc, rn1, bnd1, & rn2, bnd2, rn3, rn4, rn5, & rn6, rn7, rn8, rn9, rn10, & rn101,rn10a, rn11,rn11a, rn12 REAL, PRIVATE :: rn14, rn15,rn15a, rn16, rn17, & rn17a,rn17b,rn17c, rn18, rn18a, & rn19,rn19a,rn19b, rn20, rn20a, & rn20b, bnd3, rn21, rn22, rn23, & rn23a,rn23b, rn25,rn30a, rn30b, & rn30c, rn31, beta, rn32 REAL, PRIVATE, DIMENSION( 31 ) :: rn12a, rn12b, rn13, rn25a REAL, PRIVATE :: rn10b, rn10c, rnn191, rnn192, rn30, & rnn30a, rn33, rn331, rn332 REAL, PRIVATE, DIMENSION( 31 ) :: aa1, aa2 DATA aa1/.7939e-7, .7841e-6, .3369e-5, .4336e-5, .5285e-5, & .3728e-5, .1852e-5, .2991e-6, .4248e-6, .7434e-6, & .1812e-5, .4394e-5, .9145e-5, .1725e-4, .3348e-4, & .1725e-4, .9175e-5, .4412e-5, .2252e-5, .9115e-6, & .4876e-6, .3473e-6, .4758e-6, .6306e-6, .8573e-6, & .7868e-6, .7192e-6, .6513e-6, .5956e-6, .5333e-6, & .4834e-6/ DATA aa2/.4006, .4831, .5320, .5307, .5319, & .5249, .4888, .3894, .4047, .4318, & .4771, .5183, .5463, .5651, .5813, & .5655, .5478, .5203, .4906, .4447, & .4126, .3960, .4149, .4320, .4506, & .4483, .4460, .4433, .4413, .4382, & .4361/ REAL , PARAMETER :: xnor = 8.0e6 REAL , PARAMETER :: xnos = 1.6e7 REAL , PARAMETER :: xnoh = 2.0e5 REAL , PARAMETER :: xnog = 4.0e6 REAL , PARAMETER :: rhohail = 917. REAL , PARAMETER :: rhograul = 400. CONTAINS SUBROUTINE gsfcgce( th, & qv, ql, & qr, qi, & qs, & rho, pii, p, dt_in, z, & ht, dz8w, grav, & rhowater, rhosnow, & itimestep, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte, & rainnc, rainncv, & snownc, snowncv, sr, & graupelnc, graupelncv, & refl_10cm, diagflag, do_radar_ref, & f_qg, qg, & ihail, ice2 & ) IMPLICIT NONE INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde , & ims,ime, jms,jme, kms,kme , & its,ite, jts,jte, kts,kte INTEGER, INTENT(IN ) :: itimestep, ihail, ice2 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), & INTENT(INOUT) :: & th, & qv, & ql, & qr, & qi, & qs, & qg REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), & INTENT(IN ) :: & rho, & pii, & p, & dz8w, & z REAL, DIMENSION( ims:ime , jms:jme ), & INTENT(INOUT) :: rainnc, & rainncv, & snownc, & snowncv, & sr, & graupelnc, & graupelncv REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT):: & refl_10cm LOGICAL, OPTIONAL, INTENT(IN) :: diagflag INTEGER, OPTIONAL, INTENT(IN) :: do_radar_ref REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN) :: ht REAL, INTENT(IN ) :: dt_in, & grav, & rhowater, & rhosnow LOGICAL, INTENT(IN), OPTIONAL :: F_QG INTEGER :: itaobraun, istatmin, new_ice_sat, id INTEGER :: i, j, k INTEGER :: iskip, ih, icount, ibud, i24h REAL :: hour REAL , PARAMETER :: cmin=1.e-20 REAL :: dth, dqv, dqrest, dqall, dqall1, rhotot, a1, a2 LOGICAL :: flag_qg INTEGER:: NCALL=0 IF (NCALL .EQ. 0) THEN xam_r = 3.14159*rhowater/6. xbm_r = 3. xmu_r = 0. xam_s = 3.14159*rhosnow/6. xbm_s = 3. xmu_s = 0. if (ihail .eq. 1) then xam_g = 3.14159*rhohail/6. else xam_g = 3.14159*rhograul/6. endif xbm_g = 3. xmu_g = 0. call radar_init NCALL = 1 ENDIF itaobraun = 1 i24h=nint(86400./dt_in) if (mod(itimestep,i24h).eq.1) then write(6,*) 'ihail=',ihail,' ice2=',ice2 if (ice2.eq.0) then write(6,*) 'Running 3-ice scheme in GSFCGCE with' if (ihail.eq.0) then write(6,*) ' ice, snow and graupel' else if (ihail.eq.1) then write(6,*) ' ice, snow and hail' else write(6,*) 'ihail has to be either 1 or 0' stop endif else if (ice2.eq.1) then write(6,*) 'Running 2-ice scheme in GSFCGCE with' write(6,*) ' ice and snow' else if (ice2.eq.2) then write(6,*) 'Running 2-ice scheme in GSFCGCE with' write(6,*) ' ice and graupel' else if (ice2.eq.3) then write(6,*) 'Running warm rain only scheme in GSFCGCE without any ice' else write(6,*) 'gsfcgce_2ice in namelist.input has to be 0, 1, 2, or 3' stop endif endif new_ice_sat = 2 istatmin = 180 id = 0 ibud = 0 call fall_flux(dt_in, qr, qi, qs, qg, p, & rho, z, dz8w, ht, rainnc, & rainncv, grav,itimestep, & rhowater, rhosnow, & snownc, snowncv, sr, & graupelnc, graupelncv, & ihail, ice2, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte ) call consat_s (ihail, itaobraun) iskip = 1 if (iskip.eq.0) then call negcor(qv,rho,dz8w,ims,ime,jms,jme,kms,kme, & itimestep,1, & its,ite,jts,jte,kts,kte) call negcor(ql,rho,dz8w,ims,ime,jms,jme,kms,kme, & itimestep,2, & its,ite,jts,jte,kts,kte) call negcor(qr,rho,dz8w,ims,ime,jms,jme,kms,kme, & itimestep,3, & its,ite,jts,jte,kts,kte) call negcor(qi,rho,dz8w,ims,ime,jms,jme,kms,kme, & itimestep,4, & its,ite,jts,jte,kts,kte) call negcor(qs,rho,dz8w,ims,ime,jms,jme,kms,kme, & itimestep,5, & its,ite,jts,jte,kts,kte) call negcor(qg,rho,dz8w,ims,ime,jms,jme,kms,kme, & itimestep,6, & its,ite,jts,jte,kts,kte) endif call SATICEL_S( dt_in, IHAIL, itaobraun, ICE2, istatmin, & new_ice_sat, id, & th, qv, ql, qr, & qi, qs, qg, & rho, pii, p, itimestep, & refl_10cm, diagflag, do_radar_ref, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte & ) END SUBROUTINE gsfcgce SUBROUTINE fall_flux ( dt, qr, qi, qs, qg, p, & rho, z, dz8w, topo, rainnc, & rainncv, grav, itimestep, & rhowater, rhosnow, & snownc, snowncv, sr, & graupelnc, graupelncv, & ihail, ice2, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte ) IMPLICIT NONE INTEGER, INTENT(IN ) :: ihail, ice2, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte INTEGER, INTENT(IN ) :: itimestep REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), & INTENT(INOUT) :: qr, qi, qs, qg REAL, DIMENSION( ims:ime , jms:jme ), & INTENT(INOUT) :: rainnc, rainncv, & snownc, snowncv, sr, & graupelnc, graupelncv REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), & INTENT(IN ) :: rho, z, dz8w, p REAL, INTENT(IN ) :: dt, grav, rhowater, rhosnow REAL, DIMENSION( ims:ime , jms:jme ), & INTENT(IN ) :: topo REAL, DIMENSION( kts:kte ) :: sqrhoz REAL :: tmp1, term0 REAL :: pptrain, pptsnow, & pptgraul, pptice REAL, DIMENSION( kts:kte ) :: qrz, qiz, qsz, qgz, & zz, dzw, prez, rhoz, & orhoz INTEGER :: k, i, j REAL, DIMENSION( kts:kte ) :: vtr, vts, vtg, vti REAL :: dtb, pi, consta, constc, gambp4, & gamdp4, gam4pt5, gam4bbar REAL , PARAMETER :: & constb = 0.8, constd = 0.11, o6 = 1./6., & cdrag = 0.6 REAL , PARAMETER :: abar = 19.3, bbar = 0.37, & p0 = 1.0e5 REAL , PARAMETER :: rhoe_s = 1.29 INTEGER :: min_q, max_q REAL :: t_del_tv, del_tv, flux, fluxin, fluxout ,tmpqrz LOGICAL :: notlast dtb=dt pi=acos(-1.) consta=2115.0*0.01**(1-constb) constc=78.63*0.01**(1-constd) gambp4=ggamma(constb+4.) gamdp4=ggamma(constd+4.) gam4pt5=ggamma(4.5) gam4bbar=ggamma(4.+bbar) j_loop: do j = jts, jte i_loop: do i = its, ite pptrain = 0. pptsnow = 0. pptgraul = 0. pptice = 0. do k = kts, kte qrz(k)=qr(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/rhoz(k)) zz(k)=z(i,k,j) dzw(k)=dz8w(i,k,j) enddo DO k = kts, kte qiz(k)=qi(i,k,j) ENDDO DO k = kts, kte qsz(k)=qs(i,k,j) ENDDO IF (ice2 .eq. 0) THEN DO k = kts, kte qgz(k)=qg(i,k,j) ENDDO ELSE DO k = kts, kte qgz(k)=0. ENDDO ENDIF 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/rhoz(k)/qrz(k)) tmp1=sqrt(tmp1) vtr(k)=consta*gambp4*sqrhoz(k)/tmp1**constb vtr(k)=vtr(k)/6. if (k .eq. 1) then del_tv=amin1(del_tv,0.9*(zz(k)-topo(i,j))/vtr(k)) else del_tv=amin1(del_tv,0.9*(zz(k)-zz(k-1))/vtr(k)) endif else vtr(k)=0. endif enddo if (max_q .ge. min_q) then 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=rhoz(k)*vtr(k)*qrz(k) flux=(fluxin-fluxout)/rhoz(k)/dzw(k) qrz(k)=qrz(k)+del_tv*flux qrz(k)=amax1(0.,qrz(k)) qr(i,k,j)=qrz(k) 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/rhoz(min_q-1)/dzw(min_q-1) qr(i,min_q-1,j)=qrz(min_q-1) endif else notlast=.false. endif ENDDO 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/rhoz(k)/qsz(k)) tmp1=sqrt(tmp1) vts(k)=constc*gamdp4*sqrhoz(k)/tmp1**constd vts(k)=vts(k)/6. if (k .eq. 1) then del_tv=amin1(del_tv,0.9*(zz(k)-topo(i,j))/vts(k)) else del_tv=amin1(del_tv,0.9*(zz(k)-zz(k-1))/vts(k)) endif else vts(k)=0. endif enddo if (max_q .ge. min_q) then 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=rhoz(k)*vts(k)*qsz(k) flux=(fluxin-fluxout)/rhoz(k)/dzw(k) qsz(k)=qsz(k)+del_tv*flux qsz(k)=amax1(0.,qsz(k)) qs(i,k,j)=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/rhoz(min_q-1)/dzw(min_q-1) qs(i,min_q-1,j)=qsz(min_q-1) endif else notlast=.false. endif ENDDO if (ice2.eq.0) then 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 if (ihail .eq. 1) then min_q=min0(min_q,k) max_q=max0(max_q,k) tmp1=sqrt(pi*rhohail*xnoh/rhoz(k)/qgz(k)) tmp1=sqrt(tmp1) term0=sqrt(4.*grav*rhohail/3./rhoz(k)/cdrag) vtg(k)=gam4pt5*term0*sqrt(1./tmp1) vtg(k)=vtg(k)/6. if (k .eq. 1) then del_tv=amin1(del_tv,0.9*(zz(k)-topo(i,j))/vtg(k)) else del_tv=amin1(del_tv,0.9*(zz(k)-zz(k-1))/vtg(k)) endif else min_q=min0(min_q,k) max_q=max0(max_q,k) tmp1=sqrt(pi*rhograul*xnog/rhoz(k)/qgz(k)) tmp1=sqrt(tmp1) tmp1=tmp1**bbar tmp1=1./tmp1 term0=abar*gam4bbar/6. vtg(k)=term0*tmp1*(p0/prez(k))**0.4 if (k .eq. 1) then del_tv=amin1(del_tv,0.9*(zz(k)-topo(i,j))/vtg(k)) else del_tv=amin1(del_tv,0.9*(zz(k)-zz(k-1))/vtg(k)) endif endif else vtg(k)=0. endif enddo if (max_q .ge. min_q) then 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=rhoz(k)*vtg(k)*qgz(k) flux=(fluxin-fluxout)/rhoz(k)/dzw(k) qgz(k)=qgz(k)+del_tv*flux qgz(k)=amax1(0.,qgz(k)) qg(i,k,j)=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/rhoz(min_q-1)/dzw(min_q-1) qg(i,min_q-1,j)=qgz(min_q-1) endif else notlast=.false. endif ENDDO ENDIF 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) vti(k)= 3.29 * (rhoz(k)* qiz(k))** 0.16 if (k .eq. 1) then del_tv=amin1(del_tv,0.9*(zz(k)-topo(i,j))/vti(k)) else del_tv=amin1(del_tv,0.9*(zz(k)-zz(k-1))/vti(k)) endif else vti(k)=0. endif enddo if (max_q .ge. min_q) then 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=rhoz(k)*vti(k)*qiz(k) flux=(fluxin-fluxout)/rhoz(k)/dzw(k) qiz(k)=qiz(k)+del_tv*flux qiz(k)=amax1(0.,qiz(k)) qi(i,k,j)=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/rhoz(min_q-1)/dzw(min_q-1) qi(i,min_q-1,j)=qiz(min_q-1) endif else notlast=.false. endif ENDDO snowncv(i,j) = pptsnow snownc(i,j) = snownc(i,j) + pptsnow graupelncv(i,j) = pptgraul graupelnc(i,j) = graupelnc(i,j) + pptgraul RAINNCV(i,j) = pptrain + pptsnow + pptgraul + pptice RAINNC(i,j) = RAINNC(i,j) + pptrain + pptsnow + pptgraul + pptice sr(i,j) = 0. if (RAINNCV(i,j) .gt. 0.) sr(i,j) = (pptsnow + pptgraul + pptice) / RAINNCV(i,j) ENDDO i_loop ENDDO j_loop END SUBROUTINE fall_flux 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_gsfcgce: INPUT TO GAMMA FUNCTION TOO LARGE, X=',E12.5) WRITE(wrf_err_message,100)X CALL wrf_error_fatal3("",868,& 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 SUBROUTINE negcor ( X, rho, dz8w, & ims,ime, jms,jme, kms,kme, & itimestep, ics, & its,ite, jts,jte, kts,kte ) REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), & INTENT(INOUT) :: X REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), & INTENT(IN ) :: rho, dz8w integer, INTENT(IN ) :: itimestep, ics REAL :: A0, A1, A2 A1=0. A2=0. do k=kts,kte do j=jts,jte do i=its,ite A1=A1+max(X(i,k,j), 0.)*rho(i,k,j)*dz8w(i,k,j) A2=A2+max(-X(i,k,j), 0.)*rho(i,k,j)*dz8w(i,k,j) enddo enddo enddo A0=0.0 if (A1.NE.0.0.and.A1.GT.A2) then A0=(A1-A2)/A1 if (mod(itimestep,540).eq.0) then if (ics.eq.1) then write(61,*) 'kms=',kms,' kme=',kme,' kts=',kts,' kte=',kte write(61,*) 'jms=',jms,' jme=',jme,' jts=',jts,' jte=',jte write(61,*) 'ims=',ims,' ime=',ime,' its=',its,' ite=',ite endif if (ics.eq.1) then write(61,*) 'qv timestep=',itimestep write(61,*) ' A1=',A1,' A2=',A2,' A0=',A0 else if (ics.eq.2) then write(61,*) 'ql timestep=',itimestep write(61,*) ' A1=',A1,' A2=',A2,' A0=',A0 else if (ics.eq.3) then write(61,*) 'qr timestep=',itimestep write(61,*) ' A1=',A1,' A2=',A2,' A0=',A0 else if (ics.eq.4) then write(61,*) 'qi timestep=',itimestep write(61,*) ' A1=',A1,' A2=',A2,' A0=',A0 else if (ics.eq.5) then write(61,*) 'qs timestep=',itimestep write(61,*) ' A1=',A1,' A2=',A2,' A0=',A0 else if (ics.eq.6) then write(61,*) 'qg timestep=',itimestep write(61,*) ' A1=',A1,' A2=',A2,' A0=',A0 else write(61,*) 'wrong cloud specieis number' endif endif do k=kts,kte do j=jts,jte do i=its,ite X(i,k,j)=A0*AMAX1(X(i,k,j), 0.0) enddo enddo enddo endif END SUBROUTINE negcor SUBROUTINE consat_s (ihail,itaobraun) integer :: itaobraun real :: cn0 al = 2.5e10 cp = 1.004e7 rd1 = 1.e-3 rd2 = 2.2 cpi=4.*atan(1.) cpi2=cpi*cpi grvt=980. cd1=6.e-1 cd2=4.*grvt/(3.*cd1) tca=2.43e3 dwv=.226 dva=1.718e-4 amw=18.016 ars=8.314e7 scv=2.2904487 t0=273.16 t00=238.16 alv=2.5e10 alf=3.336e9 als=2.8336e10 avc=alv/cp afc=alf/cp asc=als/cp rw=4.615e6 cw=4.187e7 ci=2.093e7 c76=7.66 c358=35.86 c172=17.26939 c409=4098.026 c218=21.87456 c580=5807.695 c610=6.1078e3 c149=1.496286e-5 c879=8.794142 c141=1.4144354e7 if(ihail .eq. 1) then roqg=.9 ag=sqrt(cd2*roqg) bg=.5 tng=.002 else roqg=.4 ag=351.2 bg=.37 tng=.04 endif tns=.16 roqs=.1 as=78.63 bs=.11 aw=2115. bw=.8 roqr=1. tnw=.08 bgh=.5*bg bsh=.5*bs bwh=.5*bw bgq=.25*bg bsq=.25*bs bwq=.25*bw ga3b = gammagce(3.+bw) ga4b = gammagce(4.+bw) ga6b = gammagce(6.+bw) ga5bh = gammagce((5.+bw)/2.) ga3g = gammagce(3.+bg) ga4g = gammagce(4.+bg) ga5gh = gammagce((5.+bg)/2.) ga3d = gammagce(3.+bs) ga4d = gammagce(4.+bs) ga5dh = gammagce((5.+bs)/2.) ac1=aw ac2=ag ac3=as bc1=bw cc1=as dc1=bs zrc=(cpi*roqr*tnw)**0.25 zsc=(cpi*roqs*tns)**0.25 zgc=(cpi*roqg*tng)**0.25 vrc=aw*ga4b/(6.*zrc**bw) vsc=as*ga4d/(6.*zsc**bs) vgc=ag*ga4g/(6.*zgc**bg) rn1=9.4e-15 bnd1=6.e-4 rn2=1.e-3 bnd2=2.0e-3 rn3=.25*cpi*tns*cc1*ga3d esw=1. rn4=.25*cpi*esw*tns*cc1*ga3d eri=.1 rn5=.25*cpi*eri*tnw*ac1*ga3b ami=1./(24.*6.e-9) rn6=cpi2*eri*tnw*ac1*roqr*ga6b*ami esr=.5 rn7=cpi2*esr*tnw*tns*roqs esr=1. rn8=cpi2*esr*tnw*tns*roqr rn9=cpi2*tns*tng*roqs rn10=2.*cpi*tns rn101=.31*ga5dh*sqrt(cc1) rn10a=als*als/rw rn10b=alv/tca rn10c=ars/(dwv*amw) rn11=2.*cpi*tns/alf rn11a=cw/alf ami50=3.84e-6 ami40=3.08e-8 eiw=1. ui50=100. ri50=2.*5.e-3 cmn=1.05e-15 rn12=cpi*eiw*ui50*ri50**2 do 10 k=1,31 y1=1.-aa2(k) rn13(k)=aa1(k)*y1/(ami50**y1-ami40**y1) rn12a(k)=rn13(k)/ami50 rn12b(k)=aa1(k)*ami50**aa2(k) rn25a(k)=aa1(k)*cmn**aa2(k) 10 continue egw=1. rn14=.25*cpi*egw*tng*ga3g*ag egi=.1 rn15=.25*cpi*egi*tng*ga3g*ag egi=1. rn15a=.25*cpi*egi*tng*ga3g*ag egr=1. rn16=cpi2*egr*tng*tnw*roqr rn17=2.*cpi*tng rn17a=.31*ga5gh*sqrt(ag) rn17b=cw-ci rn17c=cw apri=.66 bpri=1.e-4 bpri=0.5*bpri rn18=20.*cpi2*bpri*tnw*roqr rn18a=apri rn19=2.*cpi*tng/alf rn19a=.31*ga5gh*sqrt(ag) rn19b=cw/alf rnn191=.78 rnn192=.31*ga5gh*sqrt(ac2/dva) rn20=2.*cpi*tng rn20a=als*als/rw rn20b=.31*ga5gh*sqrt(ag) bnd3=2.e-3 rn21=1.e3*1.569e-12/0.15 erw=1. rn22=.25*cpi*erw*ac1*tnw*ga3b rn23=2.*cpi*tnw rn23a=.31*ga5bh*sqrt(ac1) rn23b=alv*alv/rw if (itaobraun .eq. 0) then cn0=1.e-8 beta=-.6 elseif (itaobraun .eq. 1) then cn0=1.e-6 beta=-.46 endif rn25=cn0 rn30a=alv*als*amw/(tca*ars) rn30b=alv/tca rn30c=ars/(dwv*amw) rn31=1.e-17 rn32=4.*51.545e-4 rn30=2.*cpi*tng rnn30a=alv*alv*amw/(tca*ars) rn33=4.*tns rn331=.65 rn332=.44*sqrt(ac3/dva)*ga5dh return END SUBROUTINE consat_s SUBROUTINE saticel_s (dt, ihail, itaobraun, ice2, istatmin, & new_ice_sat, id, & ptwrf, qvwrf, qlwrf, qrwrf, & qiwrf, qswrf, qgwrf, & rho_mks, pi_mks, p0_mks,itimestep, & refl_10cm, diagflag, do_radar_ref, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte & ) IMPLICIT NONE integer, parameter :: nt=2880, nt2=2*nt integer :: itaobraun,ice2,ihail,new_ice_sat,id,istatmin integer :: itimestep real :: tairccri, cn0, dt integer ids,ide,jds,jde,kds,kde integer ims,ime,jms,jme,kms,kme integer its,ite,jts,jte,kts,kte integer i,j,k, kp real :: a0 ,a1 ,a2 ,afcp ,alvr ,ami100 ,ami40 ,ami50 ,ascp ,avcp ,betah & ,bg3 ,bgh5 ,bs3 ,bs6 ,bsh5 ,bw3 ,bw6 ,bwh5 ,cmin ,cmin1 ,cmin2 ,cp409 & ,cp580 ,cs580 ,cv409 ,d2t ,del ,dwvp ,ee1 ,ee2 ,f00 ,f2 ,f3 ,ft ,fv0 ,fvs & ,pi0 ,pir ,pr0 ,qb0 ,r00 ,r0s ,r101f ,r10ar ,r10t ,r11at ,r11rt ,r12r ,r14f & ,r14r ,r15af ,r15ar ,r15f ,r15r ,r16r ,r17aq ,r17as ,r17r ,r18r ,r19aq ,r19as & ,r19bt ,r19rt ,r20bq ,r20bs ,r20t ,r22f ,r23af ,r23br ,r23t ,r25a ,r25rt ,r2ice & ,r31r ,r32rt ,r3f ,r4f ,r5f ,r6f ,r7r ,r8r ,r9r ,r_nci ,rft ,rijl2 ,rp0 ,rr0 & ,rrq ,rrs ,rt0 ,scc ,sccc ,sddd ,see ,seee ,sfff ,smmm ,ssss ,tb0 ,temp ,ucog & ,ucor ,ucos ,uwet ,vgcf ,vgcr ,vrcf ,vscf ,zgr ,zrr ,zsr real, dimension (its:ite,jts:jte,kts:kte) :: fv real, dimension (its:ite,jts:jte,kts:kte) :: dpt, dqv real, dimension (its:ite,jts:jte,kts:kte) :: qcl, qrn, & qci, qcs, qcg real, dimension (ims:ime, kms:kme, jms:jme) :: ptwrf, qvwrf real, dimension (ims:ime, kms:kme, jms:jme) :: qlwrf, qrwrf, & qiwrf, qswrf, qgwrf real, dimension (ims:ime, kms:kme, jms:jme) :: rho_mks real, dimension (ims:ime, kms:kme, jms:jme) :: pi_mks real, dimension (ims:ime, kms:kme, jms:jme) :: p0_mks real, dimension (its:ite,jts:jte) :: & vg, zg, & ps, pg, & prn, psn, & pwacs, wgacr, & pidep, pint, & qsi, ssi, & esi, esw, & qsw, pr, & ssw, pihom, & pidw, pimlt, & psaut, qracs, & psaci, psacw, & qsacw, praci, & pmlts, pmltg, & asss, y1, y2 real, dimension (its:ite,jts:jte) :: & praut, pracw, & psfw, psfi, & dgacs, dgacw, & dgaci, dgacr, & pgacs, wgacs, & qgacw, wgaci, & qgacr, pgwet, & pgaut, pracs, & psacr, qsacr, & pgfr, psmlt, & pgmlt, psdep, & pgdep, piacr, & y5, scv, & tca, dwv, & egs, y3, & y4, ddb real, dimension (its:ite,jts:jte) :: & pt, qv, & qc, qr, & qi, qs, & qg, tair, & tairc, rtair, & dep, dd, & dd1, qvs, & dm, rq, & rsub1, col, & cnd, ern, & dlt1, dlt2, & dlt3, dlt4, & zr, vr, & zs, vs, & pssub, & pgsub, dda real, dimension (its:ite,jts:jte,kts:kte) :: rho real, dimension (kts:kte) :: & tb, qb, rho1, & ta, qa, ta1, qa1, & coef, z1, z2, z3, & am, am1, ub, vb, & wb, ub1, vb1, rrho, & rrho1, wbx real, dimension (its:ite,jts:jte,kts:kte) :: p0, pi, f0 real, dimension (kts:kte) :: & fd, fe, & st, sv, & sq, sc, & se, sqa real, dimension (kts:kte) :: & srro, qrro, sqc, sqr, & sqi, sqs, sqg, stqc, & stqr, stqi, stqs, stqg real, dimension (nt) :: & tqc, tqr, tqi, tqs, tqg real, dimension (ims:ime,jms:jme) :: & y0, ts0, qss0 integer, dimension (its:ite,jts:jte) :: it integer, dimension (its:ite,jts:jte, 4) :: ics integer :: i24h integer :: iwarm real :: r2is, r2ig REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT):: refl_10cm REAL, DIMENSION(kts:kte):: qv1d, t1d, p1d, qr1d, qs1d, qg1d, dBZ LOGICAL, OPTIONAL, INTENT(IN) :: diagflag INTEGER, OPTIONAL, INTENT(IN) :: do_radar_ref do k=kts,kte do j=jts,jte do i=its,ite rho(i,j,k)=rho_mks(i,k,j)*0.001 p0(i,j,k)=p0_mks(i,k,j)*10.0 pi(i,j,k)=pi_mks(i,k,j) dpt(i,j,k)=ptwrf(i,k,j) dqv(i,j,k)=qvwrf(i,k,j) qcl(i,j,k)=qlwrf(i,k,j) qrn(i,j,k)=qrwrf(i,k,j) qci(i,j,k)=qiwrf(i,k,j) qcs(i,j,k)=qswrf(i,k,j) qcg(i,j,k)=qgwrf(i,k,j) enddo enddo enddo do k=kts,kte do j=jts,jte do i=its,ite fv(i,j,k)=sqrt(rho(i,j,2)/rho(i,j,k)) enddo enddo enddo d2t=dt itaobraun=1 if ( itaobraun.eq.0 ) then cn0=1.e-8 elseif ( itaobraun.eq.1 ) then cn0=1.e-6 endif r2ig=1. r2is=1. if (ice2 .eq. 1) then r2ig=0. r2is=1. endif if (ice2 .eq. 2) then r2ig=1. r2is=0. endif iwarm = 0 if (ice2 .eq. 3 ) iwarm = 1 cmin=1.e-19 cmin1=1.e-20 cmin2=1.e-12 ucor=3071.29/tnw**0.75 ucos=687.97*roqs**0.25/tns**0.75 ucog=687.97*roqg**0.25/tng**0.75 uwet=4.464**0.95 rijl2 = 1. / (ide-ids) / (jde-jds) do j=jts,jte do i=its,ite it(i,j)=1 enddo enddo f2=rd1*d2t f3=rd2*d2t ft=dt/d2t rft=rijl2*ft a0=.5*istatmin*rijl2 rt0=1./(t0-t00) bw3=bw+3. bs3=bs+3. bg3=bg+3. bsh5=2.5+bsh bgh5=2.5+bgh bwh5=2.5+bwh bw6=bw+6. bs6=bs+6. betah=.5*beta r10t=rn10*d2t r11at=rn11a*d2t r19bt=rn19b*d2t r20t=-rn20*d2t r23t=-rn23*d2t r25a=rn25 ami50=3.76e-8 ami100=1.51e-7 ami40=2.41e-8 do 1000 k=kts,kte kp=k+1 tb0=0. qb0=0. do 2000 j=jts,jte do 2000 i=its,ite rp0=3.799052e3/p0(i,j,k) pi0=pi(i,j,k) pir=1./(pi(i,j,k)) pr0=1./p0(i,j,k) r00=rho(i,j,k) r0s=sqrt(rho(i,j,k)) rr0=1./rho(i,j,k) rrs=sqrt(rr0) rrq=sqrt(rrs) f0(i,j,k)=al/cp/pi(i,j,k) f00=f0(i,j,k) fv0=fv(i,j,k) fvs=sqrt(fv(i,j,k)) zrr=1.e5*zrc*rrq zsr=1.e5*zsc*rrq zgr=1.e5*zgc*rrq cp409=c409*pi0 cv409=c409*avc cp580=c580*pi0 cs580=c580*asc alvr=r00*alv afcp=afc*pir avcp=avc*pir ascp=asc*pir vrcf=vrc*fv0 vscf=vsc*fv0 vgcf=vgc*fv0 vgcr=vgc*rrs dwvp=c879*pr0 r3f=rn3*fv0 r4f=rn4*fv0 r5f=rn5*fv0 r6f=rn6*fv0 r7r=rn7*rr0 r8r=rn8*rr0 r9r=rn9*rr0 r101f=rn101*fvs r10ar=rn10a*r00 r11rt=rn11*rr0*d2t r12r=rn12*r00 r14r=rn14*rrs r14f=rn14*fv0 r15r=rn15*rrs r15ar=rn15a*rrs r15f=rn15*fv0 r15af=rn15a*fv0 r16r=rn16*rr0 r17r=rn17*rr0 r17aq=rn17a*rrq r17as=rn17a*fvs r18r=rn18*rr0 r19rt=rn19*rr0*d2t r19aq=rn19a*rrq r19as=rn19a*fvs r20bq=rn20b*rrq r20bs=rn20b*fvs r22f=rn22*fv0 r23af=rn23a*fvs r23br=rn23b*r00 r25rt=rn25*rr0*d2t r31r=rn31*rr0 r32rt=rn32*d2t*rrs pt(i,j)=dpt(i,j,k) qv(i,j)=dqv(i,j,k) qc(i,j)=qcl(i,j,k) qr(i,j)=qrn(i,j,k) qi(i,j)=qci(i,j,k) qs(i,j)=qcs(i,j,k) qg(i,j)=qcg(i,j,k) if (qc(i,j) .le. cmin1) qc(i,j)=0.0 if (qr(i,j) .le. cmin1) qr(i,j)=0.0 if (qi(i,j) .le. cmin1) qi(i,j)=0.0 if (qs(i,j) .le. cmin1) qs(i,j)=0.0 if (qg(i,j) .le. cmin1) qg(i,j)=0.0 tair(i,j)=(pt(i,j)+tb0)*pi0 tairc(i,j)=tair(i,j)-t0 zr(i,j)=zrr zs(i,j)=zsr zg(i,j)=zgr vr(i,j)=0.0 vs(i,j)=0.0 vg(i,j)=0.0 IF (IWARM .EQ. 1) THEN qi(i,j)=0.0 qs(i,j)=0.0 qg(i,j)=0.0 dep(i,j)=0. pint(i,j)=0. psdep(i,j)=0. pgdep(i,j)=0. dd1(i,j)=0. pgsub(i,j)=0. psmlt(i,j)=0. pgmlt(i,j)=0. pimlt(i,j)=0. psacw(i,j)=0. piacr(i,j)=0. psfw(i,j)=0. pgfr(i,j)=0. dgacw(i,j)=0. dgacr(i,j)=0. psacr(i,j)=0. wgacr(i,j)=0. pihom(i,j)=0. pidw(i,j)=0. if (qr(i,j) .gt. cmin1) then dd(i,j)=r00*qr(i,j) y1(i,j)=dd(i,j)**.25 zr(i,j)=zrc/y1(i,j) vr(i,j)=max(vrcf*dd(i,j)**bwq, 0.) endif pracw(i,j)=0. praut(i,j)=0.0 pracw(i,j)=r22f*qc(i,j)/zr(i,j)**bw3 y1(i,j)=qc(i,j)-bnd3 if (y1(i,j).gt.0.0) then praut(i,j)=r00*y1(i,j)*y1(i,j)/(1.2e-4+rn21/y1(i,j)) endif Y1(I,J)=QC(I,J)/D2T PRAUT(I,J)=MIN(Y1(I,J), PRAUT(I,J)) PRACW(I,J)=MIN(Y1(I,J), PRACW(I,J)) Y1(I,J)=(PRAUT(I,J)+PRACW(I,J))*D2T if (qc(i,j) .lt. y1(i,j) .and. y1(i,j) .ge. cmin2) then y2(i,j)=qc(i,j)/(y1(i,j)+cmin2) praut(i,j)=praut(i,j)*y2(i,j) pracw(i,j)=pracw(i,j)*y2(i,j) qc(i,j)=0.0 else qc(i,j)=qc(i,j)-y1(i,j) endif PR(I,J)=(PRAUT(I,J)+PRACW(I,J))*D2T QR(I,J)=QR(I,J)+PR(I,J) cnd(i,j)=0.0 tair(i,j)=(pt(i,j)+tb0)*pi0 y1(i,j)=1./(tair(i,j)-c358) qsw(i,j)=rp0*exp(c172-c409*y1(i,j)) dd(i,j)=cp409*y1(i,j)*y1(i,j) dm(i,j)=qv(i,j)+qb0-qsw(i,j) cnd(i,j)=dm(i,j)/(1.+avcp*dd(i,j)*qsw(i,j)) cnd(i,j)=max(-qc(i,j), cnd(i,j)) pt(i,j)=pt(i,j)+avcp*cnd(i,j) qv(i,j)=qv(i,j)-cnd(i,j) qc(i,j)=qc(i,j)+cnd(i,j) ern(i,j)=0.0 if(qr(i,j).gt.0.0) then tair(i,j)=(pt(i,j)+tb0)*pi0 rtair(i,j)=1./(tair(i,j)-c358) qsw(i,j)=rp0*exp(c172-c409*rtair(i,j)) ssw(i,j)=(qv(i,j)+qb0)/qsw(i,j)-1.0 dm(i,j)=qv(i,j)+qb0-qsw(i,j) rsub1(i,j)=cv409*qsw(i,j)*rtair(i,j)*rtair(i,j) dd1(i,j)=max(-dm(i,j)/(1.+rsub1(i,j)), 0.0) y1(i,j)=.78/zr(i,j)**2+r23af*scv(i,j)/zr(i,j)**bwh5 y2(i,j)=r23br/(tca(i,j)*tair(i,j)**2)+1./(dwv(i,j) & *qsw(i,j)) ern(i,j)=r23t*ssw(i,j)*y1(i,j)/y2(i,j) ern(i,j)=min(dd1(i,j),qr(i,j),max(ern(i,j),0.)) pt(i,j)=pt(i,j)-avcp*ern(i,j) qv(i,j)=qv(i,j)+ern(i,j) qr(i,j)=qr(i,j)-ern(i,j) endif ELSE if (qr(i,j) .gt. cmin1) then dd(i,j)=r00*qr(i,j) y1(i,j)=dd(i,j)**.25 zr(i,j)=zrc/y1(i,j) vr(i,j)=max(vrcf*dd(i,j)**bwq, 0.) endif if (qs(i,j) .gt. cmin1) then dd(i,j)=r00*qs(i,j) y1(i,j)=dd(i,j)**.25 zs(i,j)=zsc/y1(i,j) vs(i,j)=max(vscf*dd(i,j)**bsq, 0.) endif if (qg(i,j) .gt. cmin1) then dd(i,j)=r00*qg(i,j) y1(i,j)=dd(i,j)**.25 zg(i,j)=zgc/y1(i,j) if(ihail .eq. 1) then vg(i,j)=max(vgcr*dd(i,j)**bgq, 0.) else vg(i,j)=max(vgcf*dd(i,j)**bgq, 0.) endif endif if (qr(i,j) .le. cmin2) vr(i,j)=0.0 if (qs(i,j) .le. cmin2) vs(i,j)=0.0 if (qg(i,j) .le. cmin2) vg(i,j)=0.0 y1(i,j)=c149*tair(i,j)**1.5/(tair(i,j)+120.) dwv(i,j)=dwvp*tair(i,j)**1.81 tca(i,j)=c141*y1(i,j) scv(i,j)=1./((rr0*y1(i,j))**.1666667*dwv(i,j)**.3333333) psaut(i,j)=0.0 psaci(i,j)=0.0 praci(i,j)=0.0 piacr(i,j)=0.0 psacw(i,j)=0.0 qsacw(i,j)=0.0 dd(i,j)=1./zs(i,j)**bs3 if (tair(i,j).lt.t0) then esi(i,j)=exp(.025*tairc(i,j)) psaut(i,j)=r2is*max(rn1*esi(i,j)*(qi(i,j)-bnd1) ,0.0) psaci(i,j)=r2is*r3f*esi(i,j)*qi(i,j)*dd(i,j) psacw(i,j)=r2is*0.5*r4f*qc(i,j)*dd(i,j) praci(i,j)=r2is*r5f*qi(i,j)/zr(i,j)**bw3 piacr(i,j)=r2is*r6f*qi(i,j)*(zr(i,j)**(-bw6)) else qsacw(i,j)=r2is*r4f*qc(i,j)*dd(i,j) endif pracw(i,j)=r22f*qc(i,j)/zr(i,j)**bw3 praut(i,j)=0.0 y1(i,j)=qc(i,j)-bnd3 if (y1(i,j).gt.0.0) then praut(i,j)=r00*y1(i,j)*y1(i,j)/(1.2e-4+rn21/y1(i,j)) endif psfw(i,j)=0.0 psfi(i,j)=0.0 pidep(i,j)=0.0 if(tair(i,j).lt.t0.and.qi(i,j).gt.cmin) then y1(i,j)=max( min(tairc(i,j), -1.), -31.) it(i,j)=int(abs(y1(i,j))) y1(i,j)=rn12a(it(i,j)) y2(i,j)=rn12b(it(i,j)) psfw(i,j)=r2is* & max(d2t*y1(i,j)*(y2(i,j)+r12r*qc(i,j))*qi(i,j),0.0) rtair(i,j)=1./(tair(i,j)-c76) y2(i,j)=exp(c218-c580*rtair(i,j)) qsi(i,j)=rp0*y2(i,j) esi(i,j)=c610*y2(i,j) ssi(i,j)=(qv(i,j)+qb0)/qsi(i,j)-1. r_nci=min(1.e-6*exp(-.46*tairc(i,j)),1.) dm(i,j)=max( (qv(i,j)+qb0-qsi(i,j)), 0.) rsub1(i,j)=cs580*qsi(i,j)*rtair(i,j)*rtair(i,j) y3(i,j)=1./tair(i,j) dd(i,j)=y3(i,j)*(rn30a*y3(i,j)-rn30b)+rn30c*tair(i,j)/esi(i,j) y1(i,j)=206.18*ssi(i,j)/dd(i,j) pidep(i,j)=y1(i,j)*sqrt(r_nci*qi(i,j)/r00) dep(i,j)=dm(i,j)/(1.+rsub1(i,j))/d2t if(dm(i,j).gt.cmin2) then a2=1. if(pidep(i,j).gt.dep(i,j).and.pidep(i,j).gt.cmin2) then a2=dep(i,j)/pidep(i,j) pidep(i,j)=dep(i,j) endif psfi(i,j)=r2is*a2*.5*qi(i,j)*y1(i,j)/(sqrt(ami100) & -sqrt(ami40)) elseif(dm(i,j).lt.-cmin2) then pidep(i,j)=0. psfi(i,j)=0. else pidep(i,j)=0. psfi(i,j)=0. endif endif if(qc(i,j)+qr(i,j).lt.1.e-4) then ee1=.01 else ee1=1. endif ee2=0.09 egs(i,j)=ee1*exp(ee2*tairc(i,j)) if (tair(i,j).ge.t0) egs(i,j)=1.0 y1(i,j)=abs(vg(i,j)-vs(i,j)) y2(i,j)=zs(i,j)*zg(i,j) y3(i,j)=5./y2(i,j) y4(i,j)=.08*y3(i,j)*y3(i,j) y5(i,j)=.05*y3(i,j)*y4(i,j) dd(i,j)=y1(i,j)*(y3(i,j)/zs(i,j)**5+y4(i,j)/zs(i,j)**3 & +y5(i,j)/zs(i,j)) pgacs(i,j)=r2ig*r2is*r9r*egs(i,j)*dd(i,j) if (ihail.eq.1) then dgacs(i,j)=pgacs(i,j) else dgacs(i,j)=0. endif wgacs(i,j)=r2ig*r2is*r9r*dd(i,j) y1(i,j)=1./zg(i,j)**bg3 if(ihail .eq. 1) then dgacw(i,j)=r2ig*max(r14r*qc(i,j)*y1(i,j), 0.0) else dgacw(i,j)=r2ig*max(r14f*qc(i,j)*y1(i,j), 0.0) endif qgacw(i,j)=dgacw(i,j) y1(i,j)=abs(vg(i,j)-vr(i,j)) y2(i,j)=zr(i,j)*zg(i,j) y3(i,j)=5./y2(i,j) y4(i,j)=.08*y3(i,j)*y3(i,j) y5(i,j)=.05*y3(i,j)*y4(i,j) dd(i,j)=r16r*y1(i,j)*(y3(i,j)/zr(i,j)**5+y4(i,j)/zr(i,j)**3 & +y5(i,j)/zr(i,j)) dgacr(i,j)=r2ig*max(dd(i,j), 0.0) qgacr(i,j)=dgacr(i,j) if (tair(i,j).ge.t0) then dgacs(i,j)=0.0 wgacs(i,j)=0.0 dgacw(i,j)=0.0 dgacr(i,j)=0.0 else pgacs(i,j)=0.0 qgacw(i,j)=0.0 qgacr(i,j)=0.0 endif dgaci(i,j)=0.0 wgaci(i,j)=0.0 pgwet(i,j)=0.0 if (tair(i,j).lt.t0) then y1(i,j)=qi(i,j)/zg(i,j)**bg3 if (ihail.eq.1) then dgaci(i,j)=r2ig*r15r*y1(i,j) wgaci(i,j)=r2ig*r15ar*y1(i,j) else dgaci(i,j)=0. wgaci(i,j)=r2ig*r15af*y1(i,j) endif if (tairc(i,j).ge.-50.) then if (alf+rn17c*tairc(i,j) .eq. 0.) then write(91,*) itimestep, i,j,k, alf, rn17c, tairc(i,j) endif y1(i,j)=1./(alf+rn17c*tairc(i,j)) if (ihail.eq.1) then y3(i,j)=.78/zg(i,j)**2+r17aq*scv(i,j)/zg(i,j)**bgh5 else y3(i,j)=.78/zg(i,j)**2+r17as*scv(i,j)/zg(i,j)**bgh5 endif y4(i,j)=alvr*dwv(i,j)*(rp0-(qv(i,j)+qb0)) & -tca(i,j)*tairc(i,j) dd(i,j)=y1(i,j)*(r17r*y4(i,j)*y3(i,j) & +(wgaci(i,j)+wgacs(i,j))*(alf+rn17b*tairc(i,j))) pgwet(i,j)=r2ig*max(dd(i,j), 0.0) endif endif y1(i,j)=qc(i,j)/d2t psacw(i,j)=min(y1(i,j), psacw(i,j)) praut(i,j)=min(y1(i,j), praut(i,j)) pracw(i,j)=min(y1(i,j), pracw(i,j)) psfw(i,j)= min(y1(i,j), psfw(i,j)) dgacw(i,j)=min(y1(i,j), dgacw(i,j)) qsacw(i,j)=min(y1(i,j), qsacw(i,j)) qgacw(i,j)=min(y1(i,j), qgacw(i,j)) y1(i,j)=(psacw(i,j)+praut(i,j)+pracw(i,j)+psfw(i,j) & +dgacw(i,j)+qsacw(i,j)+qgacw(i,j))*d2t qc(i,j)=qc(i,j)-y1(i,j) if (qc(i,j) .lt. 0.0) then a1=1. if (y1(i,j) .ne. 0.0) a1=qc(i,j)/y1(i,j)+1. psacw(i,j)=psacw(i,j)*a1 praut(i,j)=praut(i,j)*a1 pracw(i,j)=pracw(i,j)*a1 psfw(i,j)=psfw(i,j)*a1 dgacw(i,j)=dgacw(i,j)*a1 qsacw(i,j)=qsacw(i,j)*a1 qgacw(i,j)=qgacw(i,j)*a1 qc(i,j)=0.0 endif wgacr(i,j)=pgwet(i,j)-dgacw(i,j)-wgaci(i,j)-wgacs(i,j) y2(i,j)=dgacw(i,j)+dgaci(i,j)+dgacr(i,j)+dgacs(i,j) if (pgwet(i,j).ge.y2(i,j)) then wgacr(i,j)=0.0 wgaci(i,j)=0.0 wgacs(i,j)=0.0 else dgacr(i,j)=0.0 dgaci(i,j)=0.0 dgacs(i,j)=0.0 endif y1(i,j)=qi(i,j)/d2t psaut(i,j)=min(y1(i,j), psaut(i,j)) psaci(i,j)=min(y1(i,j), psaci(i,j)) praci(i,j)=min(y1(i,j), praci(i,j)) psfi(i,j)= min(y1(i,j), psfi(i,j)) dgaci(i,j)=min(y1(i,j), dgaci(i,j)) wgaci(i,j)=min(y1(i,j), wgaci(i,j)) y2(i,j)=(psaut(i,j)+psaci(i,j)+praci(i,j)+psfi(i,j) & +dgaci(i,j)+wgaci(i,j))*d2t qi(i,j)=qi(i,j)-y2(i,j)+pidep(i,j)*d2t if (qi(i,j).lt.0.0) then a2=1. if (y2(i,j) .ne. 0.0) a2=qi(i,j)/y2(i,j)+1. psaut(i,j)=psaut(i,j)*a2 psaci(i,j)=psaci(i,j)*a2 praci(i,j)=praci(i,j)*a2 psfi(i,j)=psfi(i,j)*a2 dgaci(i,j)=dgaci(i,j)*a2 wgaci(i,j)=wgaci(i,j)*a2 qi(i,j)=0.0 endif dlt3(i,j)=0.0 dlt2(i,j)=0.0 if (tair(i,j).lt.t0) then if (qr(i,j).lt.1.e-4) then dlt3(i,j)=1.0 dlt2(i,j)=1.0 endif if (qs(i,j).ge.1.e-4) then dlt2(i,j)=0.0 endif endif if (ice2 .eq. 1) then dlt3(i,j)=1.0 dlt2(i,j)=1.0 endif pr(i,j)=(qsacw(i,j)+praut(i,j)+pracw(i,j)+qgacw(i,j))*d2t ps(i,j)=(psaut(i,j)+psaci(i,j)+psacw(i,j)+psfw(i,j) & +psfi(i,j)+dlt3(i,j)*praci(i,j))*d2t pg(i,j)=((1.-dlt3(i,j))*praci(i,j)+dgaci(i,j)+wgaci(i,j) & +dgacw(i,j))*d2t y1(i,j)=abs(vr(i,j)-vs(i,j)) y2(i,j)=zr(i,j)*zs(i,j) y3(i,j)=5./y2(i,j) y4(i,j)=.08*y3(i,j)*y3(i,j) y5(i,j)=.05*y3(i,j)*y4(i,j) pracs(i,j)=r2ig*r2is*r7r*y1(i,j)*(y3(i,j)/zs(i,j)**5 & +y4(i,j)/zs(i,j)**3+y5(i,j)/zs(i,j)) psacr(i,j)=r2is*r8r*y1(i,j)*(y3(i,j)/zr(i,j)**5 & +y4(i,j)/zr(i,j)**3+y5(i,j)/zr(i,j)) qsacr(i,j)=psacr(i,j) if (tair(i,j).ge.t0) then pracs(i,j)=0.0 psacr(i,j)=0.0 else qsacr(i,j)=0.0 endif pgaut(i,j)=0.0 pgfr(i,j)=0.0 if (tair(i,j) .lt. t0) then y2(i,j)=exp(rn18a*(t0-tair(i,j))) temp = 1./zr(i,j) temp = temp*temp*temp*temp*temp*temp*temp pgfr(i,j)=r2ig*max(r18r*(y2(i,j)-1.)* & temp, 0.0) endif y1(i,j)=qr(i,j)/d2t y2(i,j)=-qg(i,j)/d2t piacr(i,j)=min(y1(i,j), piacr(i,j)) dgacr(i,j)=min(y1(i,j), dgacr(i,j)) wgacr(i,j)=min(y1(i,j), wgacr(i,j)) wgacr(i,j)=max(y2(i,j), wgacr(i,j)) psacr(i,j)=min(y1(i,j), psacr(i,j)) pgfr(i,j)= min(y1(i,j), pgfr(i,j)) del=0. if(wgacr(i,j) .lt. 0.) del=1. y1(i,j)=(piacr(i,j)+dgacr(i,j)+(1.-del)*wgacr(i,j) & +psacr(i,j)+pgfr(i,j))*d2t qr(i,j)=qr(i,j)+pr(i,j)-y1(i,j)-del*wgacr(i,j)*d2t if (qr(i,j) .lt. 0.0) then a1=1. if(y1(i,j) .ne. 0.) a1=qr(i,j)/y1(i,j)+1. piacr(i,j)=piacr(i,j)*a1 dgacr(i,j)=dgacr(i,j)*a1 if (wgacr(i,j).gt.0.) wgacr(i,j)=wgacr(i,j)*a1 pgfr(i,j)=pgfr(i,j)*a1 psacr(i,j)=psacr(i,j)*a1 qr(i,j)=0.0 endif prn(i,j)=d2t*((1.-dlt3(i,j))*piacr(i,j)+dgacr(i,j) & +wgacr(i,j)+(1.-dlt2(i,j))*psacr(i,j)+pgfr(i,j)) ps(i,j)=ps(i,j)+d2t*(dlt3(i,j)*piacr(i,j) & +dlt2(i,j)*psacr(i,j)) pracs(i,j)=(1.-dlt2(i,j))*pracs(i,j) y1(i,j)=qs(i,j)/d2t pgacs(i,j)=min(y1(i,j), pgacs(i,j)) dgacs(i,j)=min(y1(i,j), dgacs(i,j)) wgacs(i,j)=min(y1(i,j), wgacs(i,j)) pgaut(i,j)=min(y1(i,j), pgaut(i,j)) pracs(i,j)=min(y1(i,j), pracs(i,j)) psn(i,j)=d2t*(pgacs(i,j)+dgacs(i,j)+wgacs(i,j) & +pgaut(i,j)+pracs(i,j)) qs(i,j)=qs(i,j)+ps(i,j)-psn(i,j) if (qs(i,j).lt.0.0) then a2=1. if (psn(i,j) .ne. 0.0) a2=qs(i,j)/psn(i,j)+1. pgacs(i,j)=pgacs(i,j)*a2 dgacs(i,j)=dgacs(i,j)*a2 wgacs(i,j)=wgacs(i,j)*a2 pgaut(i,j)=pgaut(i,j)*a2 pracs(i,j)=pracs(i,j)*a2 psn(i,j)=psn(i,j)*a2 qs(i,j)=0.0 endif y2(i,j)=d2t*(psacw(i,j)+psfw(i,j)+dgacw(i,j)+piacr(i,j) & +dgacr(i,j)+wgacr(i,j)+psacr(i,j)+pgfr(i,j)) pt(i,j)=pt(i,j)+afcp*y2(i,j) qg(i,j)=qg(i,j)+pg(i,j)+prn(i,j)+psn(i,j) psmlt(i,j)=0.0 pgmlt(i,j)=0.0 tair(i,j)=(pt(i,j)+tb0)*pi0 if (tair(i,j).ge.t0) then tairc(i,j)=tair(i,j)-t0 y1(i,j)=tca(i,j)*tairc(i,j)-alvr*dwv(i,j) & *(rp0-(qv(i,j)+qb0)) y2(i,j)=.78/zs(i,j)**2+r101f*scv(i,j)/zs(i,j)**bsh5 dd(i,j)=r11rt*y1(i,j)*y2(i,j)+r11at*tairc(i,j) & *(qsacw(i,j)+qsacr(i,j)) psmlt(i,j)=r2is*max(0.0, min(dd(i,j), qs(i,j))) if(ihail.eq.1) then y3(i,j)=.78/zg(i,j)**2+r19aq*scv(i,j)/zg(i,j)**bgh5 else y3(i,j)=.78/zg(i,j)**2+r19as*scv(i,j)/zg(i,j)**bgh5 endif dd1(i,j)=r19rt*y1(i,j)*y3(i,j)+r19bt*tairc(i,j) & *(qgacw(i,j)+qgacr(i,j)) pgmlt(i,j)=r2ig*max(0.0, min(dd1(i,j), qg(i,j))) pt(i,j)=pt(i,j)-afcp*(psmlt(i,j)+pgmlt(i,j)) qr(i,j)=qr(i,j)+psmlt(i,j)+pgmlt(i,j) qs(i,j)=qs(i,j)-psmlt(i,j) qg(i,j)=qg(i,j)-pgmlt(i,j) endif if (qc(i,j).le.cmin1) qc(i,j)=0.0 if (qi(i,j).le.cmin1) qi(i,j)=0.0 tair(i,j)=(pt(i,j)+tb0)*pi0 if(tair(i,j).le.t00) then pihom(i,j)=qc(i,j) else pihom(i,j)=0.0 endif if(tair(i,j).ge.t0) then pimlt(i,j)=qi(i,j) else pimlt(i,j)=0.0 endif pidw(i,j)=0.0 if (tair(i,j).lt.t0 .and. tair(i,j).gt.t00) then tairc(i,j)=tair(i,j)-t0 y1(i,j)=max( min(tairc(i,j), -1.), -31.) it(i,j)=int(abs(y1(i,j))) y2(i,j)=aa1(it(i,j)) y3(i,j)=aa2(it(i,j)) y4(i,j)=exp(abs(beta*tairc(i,j))) y5(i,j)=(r00*qi(i,j)/(r25a*y4(i,j)))**y3(i,j) pidw(i,j)=min(r25rt*y2(i,j)*y4(i,j)*y5(i,j), qc(i,j)) endif y1(i,j)=pihom(i,j)-pimlt(i,j)+pidw(i,j) pt(i,j)=pt(i,j)+afcp*y1(i,j)+ascp*(pidep(i,j))*d2t qv(i,j)=qv(i,j)-(pidep(i,j))*d2t qc(i,j)=qc(i,j)-y1(i,j) qi(i,j)=qi(i,j)+y1(i,j) pint(i,j)=0.0 if ( itaobraun.eq.1 ) then tair(i,j)=(pt(i,j)+tb0)*pi0 if (tair(i,j) .lt. t0) then if (qi(i,j) .le. cmin2) qi(i,j)=0. tairc(i,j)=tair(i,j)-t0 rtair(i,j)=1./(tair(i,j)-c76) y2(i,j)=exp(c218-c580*rtair(i,j)) qsi(i,j)=rp0*y2(i,j) esi(i,j)=c610*y2(i,j) ssi(i,j)=(qv(i,j)+qb0)/qsi(i,j)-1. ami50=3.76e-8 y1(i,j)=1./tair(i,j) tairccri=tairc(i,j) if(tairccri.le.-30.) tairccri=-30. y2(i,j)=exp(betah*tairccri) y3(i,j)=sqrt(qi(i,j)) dd(i,j)=y1(i,j)*(rn10a*y1(i,j)-rn10b)+rn10c*tair(i,j) & /esi(i,j) pidep(i,j)=max(r32rt*ssi(i,j)*y2(i,j)*y3(i,j)/dd(i,j), 0.e0) r_nci=min(cn0*exp(beta*tairc(i,j)),1.) dd(i,j)=max(1.e-9*r_nci/r00-qi(i,j)*1.e-9/ami50, 0.) dm(i,j)=max( (qv(i,j)+qb0-qsi(i,j)), 0.0) rsub1(i,j)=cs580*qsi(i,j)*rtair(i,j)*rtair(i,j) dep(i,j)=dm(i,j)/(1.+rsub1(i,j)) pint(i,j)=max(min(dd(i,j), dm(i,j)), 0.) pint(i,j)=min(pint(i,j)+pidep(i,j), dep(i,j)) if (pint(i,j) .le. cmin2) pint(i,j)=0. pt(i,j)=pt(i,j)+ascp*pint(i,j) qv(i,j)=qv(i,j)-pint(i,j) qi(i,j)=qi(i,j)+pint(i,j) endif endif if ( itaobraun.eq.0 ) then tair(i,j)=(pt(i,j)+tb0)*pi0 if (tair(i,j) .lt. t0) then if (qi(i,j) .le. cmin1) qi(i,j)=0. tairc(i,j)=tair(i,j)-t0 dd(i,j)=r31r*exp(beta*tairc(i,j)) rtair(i,j)=1./(tair(i,j)-c76) y2(i,j)=exp(c218-c580*rtair(i,j)) qsi(i,j)=rp0*y2(i,j) esi(i,j)=c610*y2(i,j) ssi(i,j)=(qv(i,j)+qb0)/qsi(i,j)-1. dm(i,j)=max( (qv(i,j)+qb0-qsi(i,j)), 0.) rsub1(i,j)=cs580*qsi(i,j)*rtair(i,j)*rtair(i,j) dep(i,j)=dm(i,j)/(1.+rsub1(i,j)) pint(i,j)=max(min(dd(i,j), dm(i,j)), 0.) y1(i,j)=1./tair(i,j) y2(i,j)=exp(betah*tairc(i,j)) y3(i,j)=sqrt(qi(i,j)) dd(i,j)=y1(i,j)*(rn10a*y1(i,j)-rn10b) & +rn10c*tair(i,j)/esi(i,j) pidep(i,j)=max(r32rt*ssi(i,j)*y2(i,j)*y3(i,j)/dd(i,j), 0.) pint(i,j)=pint(i,j)+pidep(i,j) pint(i,j)=min(pint(i,j),dep(i,j)) pt(i,j)=pt(i,j)+ascp*pint(i,j) qv(i,j)=qv(i,j)-pint(i,j) qi(i,j)=qi(i,j)+pint(i,j) endif endif if (new_ice_sat .eq. 0) then tair(i,j)=(pt(i,j)+tb0)*pi0 cnd(i,j)=rt0*(tair(i,j)-t00) dep(i,j)=rt0*(t0-tair(i,j)) y1(i,j)=1./(tair(i,j)-c358) y2(i,j)=1./(tair(i,j)-c76) qsw(i,j)=rp0*exp(c172-c409*y1(i,j)) qsi(i,j)=rp0*exp(c218-c580*y2(i,j)) dd(i,j)=cp409*y1(i,j)*y1(i,j) dd1(i,j)=cp580*y2(i,j)*y2(i,j) if (qc(i,j).le.cmin) qc(i,j)=cmin if (qi(i,j).le.cmin) qi(i,j)=cmin if (tair(i,j).ge.t0) then dep(i,j)=0.0 cnd(i,j)=1. qi(i,j)=0.0 endif if (tair(i,j).lt.t00) then cnd(i,j)=0.0 dep(i,j)=1. qc(i,j)=0.0 endif y5(i,j)=avcp*cnd(i,j)+ascp*dep(i,j) y1(i,j)=qc(i,j)*qsw(i,j)/(qc(i,j)+qi(i,j)) y2(i,j)=qi(i,j)*qsi(i,j)/(qc(i,j)+qi(i,j)) y4(i,j)=dd(i,j)*y1(i,j)+dd1(i,j)*y2(i,j) qvs(i,j)=y1(i,j)+y2(i,j) rsub1(i,j)=(qv(i,j)+qb0-qvs(i,j))/(1.+y4(i,j)*y5(i,j)) cnd(i,j)=cnd(i,j)*rsub1(i,j) dep(i,j)=dep(i,j)*rsub1(i,j) if (qc(i,j).le.cmin) qc(i,j)=0. if (qi(i,j).le.cmin) qi(i,j)=0. cnd(i,j)=max(-qc(i,j),cnd(i,j)) dep(i,j)=max(-qi(i,j),dep(i,j)) pt(i,j)=pt(i,j)+avcp*cnd(i,j)+ascp*dep(i,j) qv(i,j)=qv(i,j)-cnd(i,j)-dep(i,j) qc(i,j)=qc(i,j)+cnd(i,j) qi(i,j)=qi(i,j)+dep(i,j) endif if (new_ice_sat .eq. 1) then tair(i,j)=(pt(i,j)+tb0)*pi0 cnd(i,j)=rt0*(tair(i,j)-t00) dep(i,j)=rt0*(t0-tair(i,j)) y1(i,j)=1./(tair(i,j)-c358) y2(i,j)=1./(tair(i,j)-c76) qsw(i,j)=rp0*exp(c172-c409*y1(i,j)) qsi(i,j)=rp0*exp(c218-c580*y2(i,j)) dd(i,j)=cp409*y1(i,j)*y1(i,j) dd1(i,j)=cp580*y2(i,j)*y2(i,j) y5(i,j)=avcp*cnd(i,j)+ascp*dep(i,j) y1(i,j)=rt0*(tair(i,j)-t00)*qsw(i,j) y2(i,j)=rt0*(t0-tair(i,j))*qsi(i,j) if (tair(i,j).ge.t0) then dep(i,j)=0.0 cnd(i,j)=1. y2(i,j)=0. y1(i,j)=qsw(i,j) endif if (tair(i,j).lt.t00) then cnd(i,j)=0.0 dep(i,j)=1. y2(i,j)=qsi(i,j) y1(i,j)=0. endif y4(i,j)=dd(i,j)*y1(i,j)+dd1(i,j)*y2(i,j) qvs(i,j)=y1(i,j)+y2(i,j) rsub1(i,j)=(qv(i,j)+qb0-qvs(i,j))/(1.+y4(i,j)*y5(i,j)) cnd(i,j)=cnd(i,j)*rsub1(i,j) dep(i,j)=dep(i,j)*rsub1(i,j) cnd(i,j)=max(-qc(i,j),cnd(i,j)) dep(i,j)=max(-qi(i,j),dep(i,j)) pt(i,j)=pt(i,j)+avcp*cnd(i,j)+ascp*dep(i,j) qv(i,j)=qv(i,j)-cnd(i,j)-dep(i,j) qc(i,j)=qc(i,j)+cnd(i,j) qi(i,j)=qi(i,j)+dep(i,j) endif if (new_ice_sat .eq. 2) then dep(i,j)=0.0 cnd(i,j)=0.0 tair(i,j)=(pt(i,j)+tb0)*pi0 if (tair(i,j) .ge. 253.16) then y1(i,j)=1./(tair(i,j)-c358) qsw(i,j)=rp0*exp(c172-c409*y1(i,j)) dd(i,j)=cp409*y1(i,j)*y1(i,j) dm(i,j)=qv(i,j)+qb0-qsw(i,j) cnd(i,j)=dm(i,j)/(1.+avcp*dd(i,j)*qsw(i,j)) cnd(i,j)=max(-qc(i,j), cnd(i,j)) pt(i,j)=pt(i,j)+avcp*cnd(i,j) qv(i,j)=qv(i,j)-cnd(i,j) qc(i,j)=qc(i,j)+cnd(i,j) endif if (tair(i,j) .le. 258.16) then y2(i,j)=1./(tair(i,j)-c76) qsi(i,j)=rp0*exp(c218-c580*y2(i,j)) dd1(i,j)=cp580*y2(i,j)*y2(i,j) dep(i,j)=(qv(i,j)+qb0-qsi(i,j))/(1.+ascp*dd1(i,j)*qsi(i,j)) dep(i,j)=max(-qi(i,j),dep(i,j)) pt(i,j)=pt(i,j)+ascp*dep(i,j) qv(i,j)=qv(i,j)-dep(i,j) qi(i,j)=qi(i,j)+dep(i,j) endif endif psdep(i,j)=0.0 pgdep(i,j)=0.0 pssub(i,j)=0.0 pgsub(i,j)=0.0 tair(i,j)=(pt(i,j)+tb0)*pi0 if(tair(i,j).lt.t0) then if(qs(i,j).lt.cmin1) qs(i,j)=0.0 if(qg(i,j).lt.cmin1) qg(i,j)=0.0 rtair(i,j)=1./(tair(i,j)-c76) qsi(i,j)=rp0*exp(c218-c580*rtair(i,j)) ssi(i,j)=(qv(i,j)+qb0)/qsi(i,j)-1. y1(i,j)=r10ar/(tca(i,j)*tair(i,j)**2)+1./(dwv(i,j) & *qsi(i,j)) y2(i,j)=.78/zs(i,j)**2+r101f*scv(i,j)/zs(i,j)**bsh5 psdep(i,j)=r10t*ssi(i,j)*y2(i,j)/y1(i,j) pssub(i,j)=psdep(i,j) psdep(i,j)=r2is*max(psdep(i,j), 0.) pssub(i,j)=r2is*max(-qs(i,j), min(pssub(i,j), 0.)) if(ihail.eq.1) then y2(i,j)=.78/zg(i,j)**2+r20bq*scv(i,j)/zg(i,j)**bgh5 else y2(i,j)=.78/zg(i,j)**2+r20bs*scv(i,j)/zg(i,j)**bgh5 endif pgsub(i,j)=r2ig*r20t*ssi(i,j)*y2(i,j)/y1(i,j) dm(i,j)=qv(i,j)+qb0-qsi(i,j) rsub1(i,j)=cs580*qsi(i,j)*rtair(i,j)*rtair(i,j) y1(i,j)=dm(i,j)/(1.+rsub1(i,j)) psdep(i,j)=r2is*min(psdep(i,j),max(y1(i,j),0.)) y2(i,j)=min(y1(i,j),0.) pssub(i,j)=r2is*max(pssub(i,j),y2(i,j)) dd(i,j)=max((-y2(i,j)-qs(i,j)), 0.) pgsub(i,j)=r2ig*min(dd(i,j), qg(i,j), max(pgsub(i,j),0.)) if(qc(i,j)+qi(i,j).gt.1.e-5) then dlt1(i,j)=1. else dlt1(i,j)=0. endif psdep(i,j)=dlt1(i,j)*psdep(i,j) pssub(i,j)=(1.-dlt1(i,j))*pssub(i,j) pgsub(i,j)=(1.-dlt1(i,j))*pgsub(i,j) pt(i,j)=pt(i,j)+ascp*(psdep(i,j)+pssub(i,j)-pgsub(i,j)) qv(i,j)=qv(i,j)+pgsub(i,j)-pssub(i,j)-psdep(i,j) qs(i,j)=qs(i,j)+psdep(i,j)+pssub(i,j) qg(i,j)=qg(i,j)-pgsub(i,j) endif ern(i,j)=0.0 if(qr(i,j).gt.0.0) then tair(i,j)=(pt(i,j)+tb0)*pi0 rtair(i,j)=1./(tair(i,j)-c358) qsw(i,j)=rp0*exp(c172-c409*rtair(i,j)) ssw(i,j)=(qv(i,j)+qb0)/qsw(i,j)-1.0 dm(i,j)=qv(i,j)+qb0-qsw(i,j) rsub1(i,j)=cv409*qsw(i,j)*rtair(i,j)*rtair(i,j) dd1(i,j)=max(-dm(i,j)/(1.+rsub1(i,j)), 0.0) y1(i,j)=.78/zr(i,j)**2+r23af*scv(i,j)/zr(i,j)**bwh5 y2(i,j)=r23br/(tca(i,j)*tair(i,j)**2)+1./(dwv(i,j) & *qsw(i,j)) ern(i,j)=r23t*ssw(i,j)*y1(i,j)/y2(i,j) ern(i,j)=min(dd1(i,j),qr(i,j),max(ern(i,j),0.)) pt(i,j)=pt(i,j)-avcp*ern(i,j) qv(i,j)=qv(i,j)+ern(i,j) qr(i,j)=qr(i,j)-ern(i,j) endif ENDIF if (qc(i,j) .le. cmin1) qc(i,j)=0. if (qr(i,j) .le. cmin1) qr(i,j)=0. if (qi(i,j) .le. cmin1) qi(i,j)=0. if (qs(i,j) .le. cmin1) qs(i,j)=0. if (qg(i,j) .le. cmin1) qg(i,j)=0. dpt(i,j,k)=pt(i,j) dqv(i,j,k)=qv(i,j) qcl(i,j,k)=qc(i,j) qrn(i,j,k)=qr(i,j) qci(i,j,k)=qi(i,j) qcs(i,j,k)=qs(i,j) qcg(i,j,k)=qg(i,j) scc=0. see=0. dd(i,j)=max(-cnd(i,j), 0.) cnd(i,j)=max(cnd(i,j), 0.) dd1(i,j)=max(-dep(i,j), 0.)+pidep(i,j)*d2t dep(i,j)=max(dep(i,j), 0.) sccc=cnd(i,j) seee=dd(i,j)+ern(i,j) sddd=dep(i,j)+pint(i,j)+psdep(i,j)+pgdep(i,j) ssss=dd1(i,j) + pgsub(i,j) smmm=psmlt(i,j)+pgmlt(i,j)+pimlt(i,j) sfff=d2t*(psacw(i,j)+piacr(i,j)+psfw(i,j) & +pgfr(i,j)+dgacw(i,j)+dgacr(i,j)+psacr(i,j) & +wgacr(i,j))+pihom(i,j)+pidw(i,j) 2000 continue 1000 continue do k=kts,kte do j=jts,jte do i=its,ite ptwrf(i,k,j) = dpt(i,j,k) qvwrf(i,k,j) = dqv(i,j,k) qlwrf(i,k,j) = qcl(i,j,k) qrwrf(i,k,j) = qrn(i,j,k) qiwrf(i,k,j) = qci(i,j,k) qswrf(i,k,j) = qcs(i,j,k) qgwrf(i,k,j) = qcg(i,j,k) enddo enddo enddo IF ( PRESENT (diagflag) ) THEN if (diagflag .and. do_radar_ref == 1) then do j=jts,jte do i=its,ite DO K=kts,kte t1d(k)=ptwrf(i,k,j)*pi_mks(i,k,j) p1d(k)=p0_mks(i,k,j) qv1d(k)=qvwrf(i,k,j) qr1d(k)=qrwrf(i,k,j) ENDDO if (ice2.eq.0) then DO K=kts,kte qs1d(k)=qswrf(i,k,j) qg1d(k)=qgwrf(i,k,j) ENDDO elseif (ice2.eq.1) then DO K=kts,kte qs1d(k)=qswrf(i,k,j) ENDDO elseif (ice2.eq.2) then DO K=kts,kte qs1d(k)=0. qg1d(k)=qgwrf(i,k,j) ENDDO elseif (ice2.eq.3) then DO K=kts,kte qs1d(k)=0. qg1d(k)=0. ENDDO endif call refl10cm_gsfc (qv1d, qr1d, qs1d, qg1d, & t1d, p1d, dBZ, kts, kte, i, j, ihail) do k = kts, kte refl_10cm(i,k,j) = MAX(-35., dBZ(k)) enddo enddo enddo endif ENDIF return END SUBROUTINE saticel_s real function gammagce (xx) real*8 cof(6),stp,half,one,fpf,x,tmp,ser data cof,stp / 76.18009173,-86.50532033,24.01409822, & -1.231739516,.120858003e-2,-.536382e-5, 2.50662827465 / data half,one,fpf / .5, 1., 5.5 / x=xx-one tmp=x+fpf tmp=(x+half)*log(tmp)-tmp ser=one do j=1,6 x=x+one ser=ser+cof(j)/x enddo gammln=tmp+log(stp*ser) gammagce=exp(gammln) return END FUNCTION gammagce subroutine refl10cm_gsfc (qv1d, qr1d, qs1d, qg1d, & t1d, p1d, dBZ, kts, kte, ii, jj, ihail) IMPLICIT NONE INTEGER, INTENT(IN):: kts, kte, ii, jj, ihail REAL, DIMENSION(kts:kte), INTENT(IN):: & qv1d, qr1d, qs1d, qg1d, t1d, p1d REAL, DIMENSION(kts:kte), INTENT(INOUT):: dBZ REAL, DIMENSION(kts:kte):: temp, pres, qv, rho REAL, DIMENSION(kts:kte):: rr, rs, rg 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. REAL, PARAMETER:: PIx=3.1415926536 do k = kts, kte dBZ(k) = -35.0 enddo do k = kts, kte temp(k) = t1d(k) 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) N0_r(k) = xnor lamr = (xam_r*xcrg(3)*N0_r(k)/rr(k))**(1./xcre(1)) ilamr(k) = 1./lamr L_qr(k) = .true. else rr(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) = xnos 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) if (ihail.eq.1) then N0_g(k) = xnoh else N0_g(k) = xnog endif 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 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 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/PIx)*(6.0/PIx) & * (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/PIx)*(6.0/PIx) & * (xam_g/900.0)*(xam_g/900.0) & * N0_g(k)*xcgg(4)*ilamg(k)**xcge(4) enddo if (melti .and. k_0.ge.kts+1) then do k = k_0-1, kts, -1 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 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_gsfc END MODULE module_mp_gsfcgce