subroutine lcmod(pepx,tcepx,rhepx,nepx,upx,zpx,radx, + cex,alphax,iicex,icweightx,iprtx, + vvavg,vmflux) c c Version 1.0.1 c c This is a program for a 1-dimensional Lagragian cloud model. c c Input variables: c pepx(nepx): Pressure levels (hPa) of environmental sounding c ordered from lowest to highest pressure c tcepx(nepx): Temperature (deg C) at pepx levels c rhepx(nepx): Relative humidity (%) at pepx levels c nepx: No. of pepx levels c upx: Initial parcel vertical velocity (m/s) c zpx: Initial parcel height (m) c radx: Inital parcel radius (m) c cex: Entrainment coefficient (n.d.) c alphax: Precipitation removal coefficient (1/sec) c iicex: Flag for including (=1) or neglecting (=0) ice phase c icweightx: Flag for including (=1) or neglecity condensate weight c on buoyancy c iprtx: Flag for printing output (set to 0 to supress all printing c by lcmod routine except fatal errors. c c Output: c vvavg: Vertically averaged vertical velocity c vmflux: Vertically averaged convective mass flux per unit area. c c Modified: 13 Jan 2016 (v1.0.1) to add error flag to septoz routine c so lsdiag.z processing continues when sounding can not be c interpolated. c c Modified: 31 Jan 2023 (KM) for NCO specifications c c Arrays for input variables dimension pepx(nepx),tcepx(nepx),rhepx(nepx) c c Arrays for saturation vapor pressure functions parameter (mxt=500) dimension evt(0:mxt+1), eit(0:mxt+1) dimension tvp(0:mxt+1),tvpc(0:mxt+1) dimension et(0:mxt+1),et1(0:mxt+1),et2(0:mxt+1) c c Vapor pressure variables common /vpr/ tvpmx,tvpmn,dt,tvp,tvpc,evt,eit,et,et1,et2 common /vpi/ nvp c c Reference state variables common /ref/ tref,paref,eref,rhoaref,rhocref,rhovref,daref c c Physical and numerical constants common /pncon/ pi,g,ra,cpa,cva,rv,cpv,cvv,tck,tf,dtf,ce,cd,alpha c c Arrays for enivornmental sounding versus P parameter (mel=800) dimension tkep(mel),rhep(mel),ppep(mel),tvkep(mel) common /envts/ tkep,rhep,ppep,tvkep,nep c c Arrays for environmental sounding versus z parameter (mz=1000) dimension zze(mz),tze(mz),tvze(mz),pze(mz),rhze(mz) dimension wze(mz),wsze(mz),eze(mz),esze(mz),sze(mz) common /envz/ zze,tze,tvze,pze,rhze,wze,wsze,eze,esze,sze,nez c c Arrays for evenly spaced output grid parameter (mto=10000,mzo=1000) dimension tto(mto),zto(mto),pto(mto),rhoto(mto),uto(mto),rto(mto) dimension tpto(mto),teto(mto),tvpto(mto),tveto(mto) dimension wvpto(mto),wspto(mto),wveto(mto),wcpto(mto) common /tov/ tto,zto,pto,rhoto,uto,rto,tpto,teto,tvpto,tveto, + wvpto,wspto,wveto,wcpto c dimension tzo(mto),zzo(mto),pzo(mto),rhozo(mto),uzo(mto),rzo(mto) dimension tpzo(mto),tezo(mto),tvpzo(mto),tvezo(mto) dimension wvpzo(mto),wspzo(mto),wvezo(mto),wcpzo(mto) common /zov/ tzo,zzo,pzo,rhozo,uzo,rzo,tpzo,tezo,tvpzo,tvezo, + wvpzo,wspzo,wvezo,wcpzo c c Physical factor flags common /pflag/ iice,icweight c c Unit numbers data lulog,ludat /6,10/ c c ++ Parameter specification ++ c c Initialize time (t0), parcel height (m), vertical velocity (m/s) c and radius (m) tsec = 0.0 zp = zpx up = upx rp = radx c c Specify the time step (sec) and the number of time steps dtsec = 5.0 nts = 500 c c Set iupos=1 to stop model integration when vertical velocity c becomes negative iupos=1 c c Specify the time step interval for printing output ntsp = 1 c c Specify the time step label interval for printing output ntsl= 50 c c Set iice=1 to include the ice phase iice = iicex c c Set icweight=1 to include the effect of condensate weight on buoyancy icweight = icweightx c c Specify the height increment max height (m) for the evenly c spaced output and calculation grid dzo = 100.0 dzmax = 15000.0 c c ++ Model set-up ++ c Initialize physical and numerical constants call pncint c c Over-ride alpha and ce with the values from the routine c calling arguments ce = cex alpha = alphax c c Initialize saturation vapor pressure tables call svptab c c ++ Environmental sounding initialization and unit conversion ++ nep = nepx do i=1,nep ppep(i) = 100.0*pepx(i) tkep(i) = tcepx(i) + tck rhep(i) = rhepx(i) enddo c c Interpolate environmental sounding to function of z call septoz(iprtx,dzmax,ierrs) c if (ierrs .ne. 0) then vmflux = 9999. vvavg = 9999. return endif c c Initialize all time accumulation array variables to zero do i=1,mto tto(i) = 0.0 zto(i) = 0.0 pto(i) = 0.0 rhoto(i) = 0.0 uto(i) = 0.0 rto(i) = 0.0 tpto(i) = 0.0 teto(i) = 0.0 tvpto(i) = 0.0 tveto(i) = 0.0 wvpto(i) = 0.0 wspto(i) = 0.0 wveto(i) = 0.0 wcpto(i) = 0.0 enddo c c ** Start parcel initialization *** c c Get environmental p,t,tv,w and s at tsec=0 call tdenv(zp,pe,te,tve,wve,se,lulog) wce = 0.0 we = wve + wce c c Initialize parcel pressure (pp0), water mixing ratio (wp0), c and entropy mixing ratio (sp0) pp = pe wp = wve sp = se c c Specify first guess parcel T and water vapor pressure c for TD diagnostic routine tfg1 = te tfg2 = te pvfg= wp*pp/(wp + ra/rv) c c Diagnose parcel temp, water vapor and condensate mixing ratios, c water vapor pressure and total parcel density call tddiag(pp,sp,wp,tfg1,tfg2,pvfg, t1,t2,tp,wvp,wcp,pvp,rhop) call tvcal1(tp,pp,wvp,tvpar) c c Save t1 and t2 for the next time step tfg1 = t1 tfg2 = t2 c c Calculate saturation mixing ratio of the parcel call escal(tp,es,des,d2es,rl) wsp = ra*es/(rv*(pp-es)) c c Calculate initial parcel mass pmass = rhop*(4.0/3.0)*pi*(rp**3) c c Print initial parcel values ilab = 1 if (iprtx .ne. 0) then call pprint(lulog,ilab,tsec,zp,pp,rhop,up,pmass,rp, + t1,t2,tp,te,tvpar,tve,wvp,wsp,wve,wcp,wce,sp,se) endif c c Save initial parcel and environmental variables nsave = 1 call vsave(nsave,tsec,zp,pp,rhop,up,rp,tp,te,tvpar,tve, + wvp,wsp,wve,wcp) c c ** End Parcel Initialization *** c c ** Begin time integration section ** c c Calculate time step factors for Adams-Bashforth time integration dtth = 1.5*dtsec dtoh = 0.5*dtsec c do i=1,nts tsec = tsec + dtsec c c Calculate time tendency of prognostic variables call ttcal(tvpar,tve,wvp,wcp,we,up,rp,pmass,sp,se, + upt1,pmasst1,wpt1,spt1,zpt1) c c write(6,830) i,upt1,pmasst1,wpt1,spt1,zpt1 c 830 format('i,ut,pmt,wt,st,st ',i4,5(e11.4)) c if (i .eq. 1) then c Forward time step up = up + dtsec*upt1 pmass = pmass + dtsec*pmasst1 wp = wp + dtsec*wpt1 sp = sp + dtsec*spt1 zp = zp + dtsec*zpt1 else c Adams-Bashforth time step up = up + dtth*upt1 - dtoh*upt0 pmass = pmass + dtth*pmasst1 - dtoh*pmasst0 wp = wp + dtth*wpt1 - dtoh*wpt0 sp = sp + dtth*spt1 - dtoh*spt0 zp = zp + dtth*zpt1 - dtoh*zpt0 endif c c Check for negative vertical velocity if (iupos .eq. 1) then if (up .lt. 0.0) then if (iprtx .ne. 0) then write(lulog,860) zp 860 format(/,' Vertical velocity negative, ', + ' Max parcel height reached ',f7.0,' m') endif exit endif endif c if (zp .ge. zze(nez)) then if (iprtx .ne. 0) then write(lulog,861) zp/1000.0 861 format(/,' Parcel height above model top, z=',f6.2,'km') endif exit endif c c Save tendencies for the next time step upt0 = upt1 pmasst0 = pmasst1 wpt0 = wpt1 spt0 = spt1 zpt0 = zpt1 c Update diagnostic variables c c Get environmental p,t,tv,w and s at new z call tdenv(zp,pe,te,tve,wve,se,lulog) wce = 0.0 we = wve + wce c c Set parcel pressure equal to environment pressure pp = pe c c Specify first guess parcel T and water vapor pressure c for TD diagnostic routine tfg = tp pvfg= wp*pp/(wp + ra/rv) c c Diagnose parcel temp, water vapor and condensate mixing ratios, c water vapor pressure and total parcel density call tddiag(pp,sp,wp,tfg1,tfg2,pvfg, t1,t2,tp,wvp,wcp,pvp,rhop) call tvcal1(tp,pp,wvp,tvpar) c c Save t1 and t2 for the next time step tfg1 = t1 tfg2 = t2 c c write(6,831) i,t1,t2,tp,wvp,wcp,pvp,rhop c 831 format('i,t1,t2,tp,wvp,wcp,pvp,rhop ',i3,1x,7(e11.4,1x)) c c Calculate saturation mixing ratio of the parcel call escal(tp,es,des,d2es,rl) wsp = ra*es/(rv*(pp-es)) c c Diagnose parcel radius from the mass rp = ( 3.0*pmass/(4.0*pi*rhop) )**(1.0/3.0) c if (mod(i,ntsp) .eq. 0) then if (mod(i,ntsl) .eq. 0) then ilab = 1 else ilab = 0 endif c if (iprtx .ne. 0) then call pprint(lulog,ilab,tsec,zp,pp,rhop,up,pmass,rp, + t1,t2,tp,te,tvpar,tve,wvp,wsp,wve,wcp,wce,sp,se) endif c endif c c Save initial parcel and environmental variables nsave = i+1 call vsave(nsave,tsec,zp,pp,rhop,up,rp,tp,te,tvpar,tve, + wvp,wsp,wve,wcp) c enddo c c Write output variables on evenly spaced vertical grid call vezdat(ludat,nts+1,dzo,dzmax,nzo,iprtx) c c Calculate average vertical mass flux call avmf(nzo,zzo,uzo,rhozo,tpzo,tezo,vvavg,vmflux,cape,cin) c if (iprtx .ne. 0) then write(lulog,310) vvavg,vmflux,cape,cin 310 format('vvavg=',f6.2,' vmflux=',f6.2, ' cape=',f7.1, + ' cin=',f6.1) endif c return end subroutine pncint c This routine initializes needed physical and c numerical constants. MKS units are used c c Reference state variables common /ref/ tref,paref,eref,rhoaref,rhocref,rhovref,daref c c Physical and numerical constants common /pncon/ pi,g,ra,cpa,cva,rv,cpv,cvv,tck,tf,dtf,ce,cd,alpha c c Physical factor flags common /pflag/ iice,icweight c c PI pi = 3.14159265 c c Gravitational constant g = 9.81 c c Dry air properties ra = 287.01 cpa= 1005.7 cva= 718.7 c c Water vapor properties rv = 462.5 cpv= 1875.0 cvv= 1412.5 c c Conversion from Centigrade to Kelvin tck = 273.15 c c Temperature parameters for liquid to ice transition zone if (iice .eq. 1) then tf = -7.5+tck dtf= 2.5 else tf = -200.0+tck dtf= 2.5 endif c c Entrainment coefficient (nd) ce = 0.1 c c Drag coefficient (nd) cd = 0.0 c c Precipitation coefficient (1/sec) alpha = 1.0/600.0 c return end subroutine svptab c This routine calculates tables of the saturation vapor c pressure (hPa) as a function of temperature in Kelvin (tvp). c The vapor pressure over water and ice (evt, eit) are first c calculated using polynomials from Flatau (1993) and c then for a combined water-ice mixture (et) using the c method from Ooyama (1990). The 1st and 2nd derivative of c et wrt tvp are then calculated using finite difference c formulas. c c Reference state variables common /ref/ tref,paref,eref,rhoaref,rhocref,rhovref,daref c c Physical and numerical constants common /pncon/ pi,g,ra,cpa,cva,rv,cpv,cvv,tck,tf,dtf,ce,cd,alpha c c Arrays for saturation vapor pressure functions parameter (mxt=500) dimension evt(0:mxt+1),eit(0:mxt+1) dimension tvp(0:mxt+1),tvpc(0:mxt+1) dimension et(0:mxt+1),et1(0:mxt+1),et2(0:mxt+1) c common /vpr/ tvpmx,tvpmn,dt,tvp,tvpc,evt,eit,et,et1,et2 common /vpi/ nvp c c Latent heat array dimension rlt(0:mxt+1) common /dvp/ rlt c c Arrays for polynomial coefficients dimension a(0:6),b(0:6) c c Arrays for temperatures weights for water-air mixing dimension ww(0:mxt+1),wi(0:mxt+1),dwwdt(0:mxt+1) dimension enum(0:mxt+1) c c Array for E smoothing dimension evf(0:mxt+1) c c Specify the start, ending and increment of c temperature (K) for the vp tables c tvpmn = -120.0+tck tvpmx = 50.0+tck dt = 0.5 dt2 = 2.0*dt dts = dt*dt nvp = 1 + ifix( 1.0e-4 + (tvpmx-tvpmn)/dt ) c c Set nfilt to number of times to smooth et(n) nfilt=1 c do n=0,nvp+1 tvp(n) = tvpmn + dt*float(n-1) tvpc(n)= tvp(n)-tck enddo c c Calculate water and ice weights as function of T do n=0,nvp+1 tn = (tvp(n)-tf)/dtf ww(n) = 0.5*(1.0+tanh(tn)) wi(n) = 1.0-ww(n) c dwwdt(n) = (1.0-(tanh(tn))**2)/(2.0*dtf) enddo c c Specify polynomial coefficients for calculation c of evt and eit. a(0) = 6.11176750 a(1) = 0.443986062 a(2) = 0.143053301e-01 a(3) = 0.265027242e-03 a(4) = 0.302246994e-05 a(5) = 0.203886313e-07 a(6) = 0.638780966e-10 c b(0) = 6.10952665 b(1) = 0.501948366 b(2) = 0.186288989e-01 b(3) = 0.403488906e-03 b(4) = 0.539797852e-05 b(5) = 0.420713632e-07 b(6) = 0.147271071e-09 c c Speficy reference state T and P tref = 273.15 paref = 1.0e+5 c c Calculate ev and ei at the reference state evref = 0.0 eiref = 0.0 trefc = tref - tck c do k=0,6 ttk = trefc**k evref = evref + a(k)*ttk eiref = eiref + b(k)*ttk enddo c evref = 100.0*evref eiref = 100.0*eiref c c Calculate evt and eit from polynomial series c c Specify minimum T value (C) for utilizing tables. c Asympototic exponential forms are used for colder T. c Also, specify the width of the blending zone (tblend). ttmin = -50.0 tblend= 1.0 c do n=0,nvp+1 tt = (tvpc(n)) c evt(n) = 0.0 eit(n) = 0.0 c c Asymptotic form evt2 = 1459.9*exp(.1083*tt) c eit2 = 1464.2*exp(.1182*tt) c tti = -50.0-tt c if (tti .ge. 0.0) then c eit2 = 3.942*exp(-.106*(tti)**1.09) c else c eit2 = 3.942*exp( .106*(-tti)**1.09) c endif c ttk = tt + tck eit2 = 3.942*exp(6127.6*(1.0/223.15 - 1.0/ttk)) c c General form evt1 = 0.0 eit1 = 0.0 do k=0,6 ttk = tt**k evt1 = evt1 + a(k)*ttk eit1 = eit1 + b(k)*ttk enddo c c Convert to hPa evt1 = 100.0*evt1 eit1 = 100.0*eit1 c w1 = 0.5*(1.0 + tanh( (tt-ttmin)/tblend )) w2 = 1.0 - w1 c c write(6,*) tt,w1,w2,evt1,evt2,eit1,eit2 c evt(n) = w1*evt1 + w2*evt2 eit(n) = w1*eit1 + w2*eit2 enddo c c Initialize water-ice vp with water vp do n=0,nvp+1 et(n) = evt(n) enddo c c Specify index of temperature array to begin c integration for water-ice mixture vp calculation nts = nvp-10 c ts = tvp(nts) evs = evt(nts) eis = eit(nts) c c Analytic part of e(t) calculation do n=nts,0,-1 et(n) = alog(evref) + (ts/tvp(n))*alog(evs/evref) + + ww(n)*alog(evt(n)/evref) + - ww(nts)*(ts/tvp(n))*alog(evs/evref) + + wi(n)*alog(eit(n)/eiref) + - wi(nts)*(ts/tvp(n))*alog(eis/eiref) enddo c c Calculate the integral term in the e(t) calculation etm = 0.0 enum(nts) = 0.0 do n=nts-1,0,-1 dtm = tvp(n)-tvp(n+1) tm1 = tvp(n+1)*dwwdt(n+1)*(alog(eit(n+1)/eiref) + -alog(evt(n+1)/evref)) tm2 = tvp(n )*dwwdt(n )*(alog(eit(n )/eiref) + -alog(evt(n )/evref)) etm = dtm*0.5*(tm1+tm2) enum(n) = enum(n+1) + etm enddo c c Add integral term to analytic term do n=nts,0,-1 et(n) = et(n) + enum(n)/tvp(n) enddo c c Take exponential to find evt do n=nts,0,-1 et(n) = exp(et(n)) enddo c if (nfilt .gt. 0) then do k=1,nfilt evf( 0) = et( 0) evf(nts) = et(nts) c do n=1,nts-1 evf(n) = 0.25*et(n-1) + 0.5*et(n) + 0.25*et(n+1) enddo c do n=1,nts-1 et(n) = evf(n) enddo enddo c endif c et1(0) = 0.0 et2(0) = 0.0 et1(nvp+1) = 0.0 et2(nvp+1) = 0.0 c c Calculate first and second derivative of et do n=1,nvp et1(n) = (et(n+1)-et(n-1))/dt2 et2(n) = (et(n+1)+et(n-1)-2.0*et(n))/dts enddo c c Calculate latent heat of water-ice mixture c from C-C equation rlt(0) = 0.0 rlt(nvp+1) = 0.0 do n=1,nvp rlt(n) = et1(n)*rv*tvp(n)*tvp(n)/et(n) enddo c iprint = 0 if (iprint .eq. 1) then write(6,801) 801 format(' T Wv Ev Ei E ', + ' dE/dT d2E/dT2 L') c do n=1,nvp write(6,800) n,tvpc(n),ww(n),evt(n),eit(n), + et(n),et1(n),et2(n),rlt(n) 800 format(1x,i4,f6.1,1x,f4.2,1x,e10.4,1x,e10.4,4(1x,e10.4)) enddo endif c c Specify the rest of the reference state variables eref = evref rhoaref = paref/(ra*tref) rhocref = 0.0 rhovref = eref/(rv*tref) c call escal(tref,es,des,d2es,rl) daref = rl/tref c return end subroutine tvcal(t,p,rh,tv) c This routine calculates the virtual temperature tv (K) c given the temperature t (K), pressure p (Pa), and relative c humidity rh (%) c c Physical and numerical constants common /pncon/ pi,g,ra,cpa,cva,rv,cpv,cvv,tck,tf,dtf,ce,cd,alpha c call escal(t,es,des,d2es,rl) ws = ra*es/(rv*(p-es)) w = ws*rh/100.0 c1 = (rv-ra)/ra v1 = w/(1.0+w) tv = t*(1.0 + c1*v1) c return end subroutine tvcal1(t,p,w,tv) c This routine calculates the virtual temperature tv (K) c given the temperature t (K), pressure p (Pa), and c water vapor mixing ratio wv (kg/kg) c c Physical and numerical constants common /pncon/ pi,g,ra,cpa,cva,rv,cpv,cvv,tck,tf,dtf,ce,cd,alpha c c1 = (rv-ra)/ra v1 = w/(1.0+w) tv = t*(1.0 + c1*v1) c return end subroutine escal(t,es,des,d2es,rl) c This routine calculates the saturation vapor c pressure (es in Pa) from the temperature (t in K) by c interpolating from a pre-calculated table. c c The first (des) and second (d2es) derivatives of es wrt to temperature c and the latent heat (rl) are also calcualted. c c Arrays for saturation vapor pressure functions parameter (mxt=500) dimension evt(0:mxt+1), eit(0:mxt+1) dimension tvp(0:mxt+1),tvpc(0:mxt+1) dimension et(0:mxt+1),et1(0:mxt+1),et2(0:mxt+1) c c Latent heat array common /dvp/ rlt dimension rlt(0:mxt+1) c c Physical and numerical constants common /pncon/ pi,g,ra,cpa,cva,rv,cpv,cvv,tck,tf,dtf,ce,cd,alpha c c Vapor pressure variables common /vpr/ tvpmx,tvpmn,dt,tvp,tvpc,evt,eit,et,et1,et2 common /vpi/ nvp c if (t .ge. tvpmx) then es = et(0) des = et1(0) d2es = et2(0) rl = rlt(0) elseif (t .le. tvpmn) then es = et(nvp) des = et1(nvp) d2es= et2(nvp) rl = rlt(nvp) else c Find index closest to but less than tt id1 = ifix( (t-tvpmn)/dt ) id2 = id1 + 1 c e1 = et(id1) e2 = et(id2) f1 = et1(id1) f2 = et1(id2) g1 = et2(id1) g2 = et2(id2) r1 = rlt(id1) r2 = rlt(id2) c w2 = (t-tvp(id1))/dt w1 = 1.0 - w2 c es = w1*e1 + w2*e2 des = w1*f1 + w2*f2 d2es = w1*g1 + w2*g2 rl = w1*r1 + w2*r2 endif c return end subroutine scal(p,t,w, sa,sm,s,wc) c This routine calculates the dry air and moisture mass-specific c entropies (sa and sm) and the total air entropy mixing ratio (s) c given the total pressure p, temperature t and total water mixing ratio (w). c The mixing ratio of the condensed water is also calculated if c the air is saturated. MKS units are used for all variables. c c Physical and numerical constants common /pncon/ pi,g,ra,cpa,cva,rv,cpv,cvv,tck,tf,dtf,ce,cd,alpha c c Reference state variables common /ref/ tref,paref,eref,rhoaref,rhocref,rhovref,daref c c Calculate the saturation vapor pressure (e) call escal(t,es,des,d2es,rl) c c Calculate the vapor pressure from the total mixing ratio pv = w*p/(w + ra/rv) c if (pv .le. es) then c The air is unsaturated wc = 0.0 pa = p - pv rhoa = pa/(ra*t) rhov = pv/(rv*t) rhom = rhov c sa = cva*alog(t/tref) - ra*alog(rhoa/rhoaref) sm = cvv*alog(t/tref) - rv*alog(rhov/rhovref) + daref else c The air is saturated pa = p - es rhoa = pa/(ra*t) rhov = es/(rv*t) rhom = w*rhoa c ws = rhov/rhoa wc = w - ws c sa = cva*alog(t/tref) - ra*alog(rhoa/rhoaref) sm = cvv*alog(t/tref) - rv*alog(rhov/rhovref) + daref - + (rl/t)*(1.0-(rhov/rhom)) c endif c s = (sa*rhoa + sm*rhom)/rhoa c c write(6,*) t,p,pv,rhoa,rhom return end subroutine scal1(p,t,w, sa1,sa2,sm1,sm2,s1,s2,wc) c This routine calculates the dry air and moisture mass-specific c entropies (sa and sm) and the total air entropy mixing ratio (s) c given the total pressure p, temperature t and total water mixing ratio (w). c The mixing ratio of the condensed water is also calculated if c the air is saturated. MKS units are used for all variables. c c In this version both sm1,sm2 and s1,s2 are returned (the unsaturated c and saturated versions of sm and s) for testing purposes. c c Physical and numerical constants common /pncon/ pi,g,ra,cpa,cva,rv,cpv,cvv,tck,tf,dtf,ce,cd,alpha c c Reference state variables common /ref/ tref,paref,eref,rhoaref,rhocref,rhovref,daref c c Calculate the saturation vapor pressure (e) call escal(t,es,des,d2es,rl) c c Calculate the vapor pressure from the total mixing ratio pv = w*p/(w + ra/rv) c c Unsaturated case wc = 0.0 pa = p - pv rhoa = pa/(ra*t) rhov = pv/(rv*t) rhom = rhov c sa1 = cva*alog(t/tref) - ra*alog(rhoa/rhoaref) sm1 = cvv*alog(t/tref) - rv*alog(rhov/rhovref) + daref s1 = (sa1*rhoa + sm1*rhom)/rhoa c c Saturated case pa = p - es rhoa = pa/(ra*t) rhov = es/(rv*t) rhom = w*rhoa c ws = rhov/rhoa wc = w - ws c sa2 = cva*alog(t/tref) - ra*alog(rhoa/rhoaref) sm2 = cvv*alog(t/tref) - rv*alog(rhov/rhovref) + daref - + (rl/t)*(1.0-(rhov/rhom)) s2 = (sa2*rhoa + sm2*rhom)/rhoa c return end subroutine septoz(iprtx,dzmax,ierr) c This routine converts the environmental sounding c from a function of P to a function of z, and calcualtes c needed environmental moisture variables. c c Arrays for enivornmental sounding versus P parameter (mel=800) dimension tkep(mel),rhep(mel),ppep(mel),tvkep(mel) common /envts/ tkep,rhep,ppep,tvkep,nep c c Arrays for environmental sounding versus z parameter (mz=1000) dimension zze(mz),tze(mz),tvze(mz),pze(mz),rhze(mz) dimension wze(mz),wsze(mz),eze(mz),esze(mz),sze(mz) common /envz/ zze,tze,tvze,pze,rhze,wze,wsze,eze,esze,sze,nez c dimension rlze(mz) c c Physical and numerical constants common /pncon/ pi,g,ra,cpa,cva,rv,cpv,cvv,tck,tf,dtf,ce,cd,alpha c c Local variables dimension zp(mel),saze(mz),smze(mz) c c Specify min and max z values (meters) zb = 0.0 c zt = 16.0e+3 zt = dzmax + 500.0 dz = 0.5e+3 nez = 1 + ifix(0.0001 + (zt-zb)/dz) c do k=1,nez zze(k) = zb + dz*float(k-1) enddo c c Calculate Tv at pressure levels do n=1,nep call tvcal(tkep(n),ppep(n),rhep(n),tvkep(n)) enddo c c Calculate the height of the environmental pressure levels zp(nep) = 0.0 do n=nep-1,1,-1 p2 = ppep(n ) p1 = ppep(n+1 ) tv2 = tvkep(n ) tv1 = tvkep(n+1) if (tv1 .eq. tv2) tv1 = tv1 + 0.001 c zp(n) = zp(n+1) + (ra/g)*(tv2-tv1)*alog(p1/p2)/alog(tv2/tv1) enddo c if (iprtx .ne. 0) then do n=1,nep write(6,301) ppep(n),zp(n),tkep(n)-tck,tvkep(n)-tck,rhep(n) 301 format(1x,'P=',f7.0,' Z=',f7.0,' T=',f6.1,' Tv=',f6.1, + ' RH',f6.1) enddo endif c c Calculate T and P as a function of height do k=1,nez c c Find the environmental pressure levels above and below the c current height level do n=1,nep-1 if (zze(k) .lt. zp(n) .and. zze(k) .ge. zp(n+1)) then tv1 = tvkep(n+1) tv2 = tvkep(n ) t1 = tkep(n+1) t2 = tkep(n ) p1 = ppep(n+1) p2 = ppep(n ) z1 = zp(n+1) z2 = zp(n ) r1 = rhep(n+1) r2 = rhep(n ) c if (tv1 .eq. tv2) tv1 = tv1 + 0.001 c exit !found pressure levels, leave inner loop endif if (n .eq. nep-1) then c Reached end of loop without finding pressure levels write(6,900) zze(k) 900 format(' z=',f8.1,' out of z domain in routine septoz') c stop ierr=1 return endif enddo c gammav= (tv2-tv1)/(z2-z1) gammat= (t2 - t1)/(z2-z1) gammar= (r2 - r1)/(z2-z1) aa = -g/(ra*gammav) c tvze(k) = tv1 + gammav*(zze(k)-z1) tze(k) = t1 + gammat*(zze(k)-z1) rhze(k) = r1 + gammar*(zze(k)-z1) pze(k) = p1*(tvze(k)/tv1)**aa c call escal(tze(k),esze(k),des,d2es,rl) wsze(k) = ra*esze(k)/(rv*(pze(k)-esze(k))) wze(k) = wsze(k)*rhze(k)/100.0 eze(k) = wze(k)*pze(k)/(wze(k) + ra/rv) rlze(k) = rl c call scal(pze(k),tze(k),wze(k),sa,sm,sze(k),wc) saze(k) = sa smze(k) = sm c enddo c if (iprtx .ne. 0) then write(6,299) 299 format(' Z T Tv P RH ws w e', + ' Sa Sm S L') do k=nez,1,-1 write(6,300) zze(k)/1000.0,tze(k)-tck,tvze(k)-tck, + pze(k)/100.0,rhze(k),wsze(k)*1000.,wze(k)*1000., + eze(k)/100.0,saze(k),smze(k),sze(k),rlze(k) c 300 format( 'Z=',f6.1,' T= ',f6.1,' Tv=',f6.1, c + ' P=',f6.1,' RH=',f6.1,' ws=',f6.2,1x, c + ' w=',f6.2,' e=',f6.2,' sa=',f6.1,' sm=',f6.1, c + ' s=',f6.1,' L=',e10.3) 300 format(5(f6.1,1x),3(f6.2,1x),3(f7.1,1x),e10.3) enddo endif c ierr=0 return end subroutine tdenv(z,pe,te,tve,we,se,lulog) c This routine calculates the thermodynamic properties c of the parcel environment at level z. It is assumed that c the parcel environment is unsaturated. c c Reference state variables common /ref/ tref,paref,eref,rhoaref,rhocref,rhovref,daref c c Physical and numerical constants common /pncon/ pi,g,ra,cpa,cva,rv,cpv,cvv,tck,tf,dtf,ce,cd,alpha c c Arrays for environmental sounding versus z parameter (mz=1000) dimension zze(mz),tze(mz),tvze(mz),pze(mz),rhze(mz) dimension wze(mz),wsze(mz),eze(mz),esze(mz),sze(mz) common /envz/ zze,tze,tvze,pze,rhze,wze,wsze,eze,esze,sze,nez c c Find height level closest to z do k=1,nez-1 if (z .ge. zze(k) .and. z .le. zze(k+1)) then p1 = pze(k) p2 = pze(k+1) t1 = tze(k) t2 = tze(k+1) tv1 = tvze(k) tv2 = tvze(k+1) r1 = rhze(k) r2 = rhze(k+1) e1 = eze(k) e2 = eze(k+1) w1 = wze(k) w2 = wze(k+1) z1 = zze(k) z2 = zze(k+1) se1 = sze(k) se2 = sze(k+1) c exit endif if (k .eq. nez-1) then c Reached the end of the loop without finding heights write(lulog,900) z 900 format(/,'z outside of model domain, z=',f8.1) stop endif enddo c a1 = (z2-z)/(z2-z1) a2 = (z-z1)/(z2-z1) c pe = a1*p1 + a2*p2 te = a1*t1 + a2*t2 tve= a1*tv1+ a2*tv2 we = a1*w1 + a2*w2 se = a1*se1+ a2*se2 c return end subroutine tddiag(p,s,w,tfg1,tfg2,pvfg, t1,t2,t,wv,wc,pv,rho) c The routine diagnoses the temperature t, water vapor mixing c ratio (wv), water vapor condensate (wc) mixing ration, water vapor pressure c (pv) and total air density (rho), given the total pressure (p), entropy c mixing ratio (s), and total water mixing ratio (w). c c First guess temperatures for t1 and t2 (tfg1 and tfg2) and c water vapor pressure (pvfg) need to be provided. c c Physical and numerical constants common /pncon/ pi,g,ra,cpa,cva,rv,cpv,cvv,tck,tf,dtf,ce,cd,alpha c c Specify the number of iterations nit = 10 c c Find T1 wv1 = w wc1 = 0.0 c pv1 = wv1*p/(wv1 + ra/rv) pa = p - pv1 c t1 = tfg1 do i=1,nit rhoa = pa/(ra*t1) call t1froms(s,w,rhoa, t1) enddo c c Find T2 c t2 = tfg2 do i=1,nit call escal(t2,es,des,d2es,rl) pv2 = es pa = p - pv2 rhoa = pa/(ra*t2) wv2 = es/(rv*t2*rhoa) wc2 = w - wv2 tstart = t2 c call t2froms(s,w,rhoa,tstart, t2) c write(6,*) 'i98,tstart,t2 ',i,tstart,t2 enddo c if (t1 .ge. t2) then c if (t1 .ge. 0) then t = t1 wv = wv1 wc = wc1 pv = pv1 else t = t2 wv = wv2 wc = wc2 pv = pv2 endif c c Calculate final value of total air density pa = p - pv rhoa = pa/(ra*t) rhov = wv*rhoa rhoc = wc*rhoa rho = rhoa + rhov + rhoc c c call scal1(p,t1,w, sa11,sa12,sm11,sm12,s11,s12,wc1) c call scal1(p,t2,w, sa21,sa22,sm21,sm22,s21,s22,wc2) c c write(6,800) t1,t2,rhoa,rhov,rhoc,rho,pv c 800 format('t1,t2,rhoa,rhov,rhoc,rho,pv: ', c + f7.2,1x,f7.2,1x,4(e11.4,1x),1x,f6.1) c c write(6,801) sa11,sm11,s11,sa12,sm12,s12,wc1, c + sa21,sm21,s21,sa22,sm22,s22,wc2 c 801 format('sa11,sm11,s11,sa12,sm12,s12,wc1: ',7(e11.4,1x),/, c + 'sa21,sm21,s21,sa22,sm22,s22,wc2: ',7(e11.4,1x)) c return end subroutine t1froms(s,w,rhoa, t1) c This routine calculates the temperature t1 from the c entropy density s, total water mixing ratio w and dry c air density rhoa assuming unsaturated air, following c Ooyama (1990). c c Reference state variables common /ref/ tref,paref,eref,rhoaref,rhocref,rhovref,daref c c Physical and numerical constants common /pncon/ pi,g,ra,cpa,cva,rv,cpv,cvv,tck,tf,dtf,ce,cd,alpha c eta = w*rhoa c c Calculate t1 from the analytic formula fnum = s + ra*alog(rhoa/rhoaref) + w*rv*alog(eta/rhovref) - + w*daref fden = cva + w*cvv t1 = tref*exp(fnum/fden) c return end subroutine t2froms(s,w,rhoa,tstart, t2) c This routine calculates the temperature t from the c entropy density s, total water mixing ratio w and dry c air density rhoa using an iterative procedure. The input c variable tstart is the first guess for the iteration. c The formulation follows Ooyama (1990). c c Reference state variables common /ref/ tref,paref,eref,rhoaref,rhocref,rhovref,daref c c Physical and numerical constants common /pncon/ pi,g,ra,cpa,cva,rv,cpv,cvv,tck,tf,dtf,ce,cd,alpha c eta = w*rhoa c c Calculate t2 using an iterative procedure t2 = tstart nit = 10 do i=1,nit t = t2 sa = cva*alog(t/tref) - ra*alog(rhoa/rhoaref) c call escal(t,es,des,d2es,rl) etas = es/(rv*t) c sm1 = cvv*alog(t/tref) sm2 = -rv*alog(etas/rhovref) sm3 = daref sm4 = (etas-eta)*(rl/t)/eta sm = sm1 + sm2 + sm3 + sm4 c dsadt = cva/t c dsm1 = cpv/t dsm2 = (rv/es)*des*((t/es)*des-2.0) dsm3 = d2es*((1.0/eta) - (rv*t)/es) dsmdt = dsm1 + dsm2 + dsm3 c c siter = sa + w*sm c write(6,750) i,t,sa,sm,siter,s c 750 format('i,t,sa,sm,siter,s ',i2,1x,f6.2,4(e11.4)) c f = sa + w*sm - s dfdt = dsadt + w*dsmdt c c write(6,748) t,etas,eta,rl,f,dfdt c 748 format(' t,etas,eta,rl,f,dfdt: ',3x,6(e12.5)) c write(6,749) i,t,sm1,sm2,sm3,sm4,sm c 749 format('i,t,sm1,2,3,4,sm: ',i2,1x,6(e12.5)) t2 = t - f/dfdt enddo c return end subroutine pprint(lulog,ilab,tsec,zp,pp,rhop,up,pmass,rp, + t1,t2,tp,te,tvp,tve,wvp,wsp,wve,wcp,wce,sp,se) c c This routine prints the parcel and environmental variables c at a single time c common /pncon/ pi,g,ra,cpa,cva,rv,cpv,cvv,tck,tf,dtf,ce,cd,alpha c c Unit conversions for printing ppw = pp/100.0 c t1w = t1 -tck t2w = t2 -tck tpw = tp -tck tew = te -tck tvpw = tvp-tck tvew = tve-tck c wvpw = wvp*1000.0 wvew = wve*1000.0 wcpw = wcp*1000.0 wcew = wce*1000.0 c wspw = wsp*1000.0 c if (ilab .eq. 1) then write(lulog,100) 100 format(/,' t z P rhop u M Rp', + ' T1 T2 Tp Te Tvp Tve', + ' wvp wsp wve wcp wce Sp Se') endif c write(lulog,300) tsec,zp,ppw,rhop,up,pmass,rp, + t1w,t2w,tpw,tew,tvpw,tvew, + wvpw,wspw,wvew,wcpw,wcew,sp,se 300 format(f6.1,1x,f6.0,1x,f6.1,1x,f6.2,1x,f6.1,1x,e10.3,1x,f6.1,1x, + 6(f6.1,1x), + 5(f6.1,1x),2(f6.1,1x)) c return end subroutine ttcal(tvp,tve,wvp,wcp,we,up,rp,pmass,sp,se, + upt,pmasst,wpt,spt,zpt) c c This routine calculates the time tendency of the 5 prognostic variables c of the cloud model: c c upt: Vertical velocity tendency c pmasst: Parcel total mass tendency c wpt: Parcel total water mixing ratio c dpt: Parcel entropy mixing ratio tendency c zpt: Parcel height tendency c c Physical and numerical constants common /pncon/ pi,g,ra,cpa,cva,rv,cpv,cvv,tck,tf,dtf,ce,cd,alpha c c Physical factor flags common /pflag/ iice,icweight c c ++ Preliminary calculations if (rp .le. 0.0) then ri = 0.0 else ri = 1.0/rp endif c etfac = abs(up)*ce*ri c c ++ upt calculation f1 = (tvp-tve)/tve c if (icweight .eq. 1) then f2 = wcp/(1.0+wvp) else f2 = 0.0 endif c upt = g*( (f1-f2)/(1.0+f2) ) - (ce+cd)*up*up*ri c c ++ pmasst calculation pmasst = etfac*pmass c c ++ wpt calculation wp = wvp+wcp wpt = etfac*(we-wp) - alpha*wcp c c ++ spt calculation spt = etfac*(se-sp) c c ++ zpt calculation zpt = up c c write(6,*) 'upt,pmasst,wpt,spt,zpt',upt,pmasst,wpt,spt,zpt return end subroutine vsave(nsave,tsec,zp,pp,rhop,up,rp,tp,te,tvp,tve, + wvp,wsp,wve,wcp) c This routine saves the primary parcel and environmental variables c as a function of time for later interpolation to an evenly spaced c vertical grid. c c Arrays for evenly spaced output grid parameter (mto=10000,mzo=1000) dimension tto(mto),zto(mto),pto(mto),rhoto(mto),uto(mto),rto(mto) dimension tpto(mto),teto(mto),tvpto(mto),tveto(mto) dimension wvpto(mto),wspto(mto),wveto(mto),wcpto(mto) common /tov/ tto,zto,pto,rhoto,uto,rto,tpto,teto,tvpto,tveto, + wvpto,wspto,wveto,wcpto c n = nsave c tto(n) = tsec zto(n) = zp pto(n) = pp rhoto(n) = rhop uto(n) = up rto(n) = rp tpto(n) = tp teto(n) = te tvpto(n) = tvp tveto(n) = tve wvpto(n) = wvp wspto(n) = wsp wveto(n) = wve wcpto(n) = wcp c return end subroutine vezdat(ludat,ntmax,dzo,dzmax,nzo,iprtx) c This routine prints the primary variables as a function of z c on an evenly spaced grid. c c Physical and numerical constants common /pncon/ pi,g,ra,cpa,cva,rv,cpv,cvv,tck,tf,dtf,ce,cd,alpha c c Arrays for evenly spaced output grid parameter (mto=10000,mzo=1000) dimension tto(mto),zto(mto),pto(mto),rhoto(mto),uto(mto),rto(mto) dimension tpto(mto),teto(mto),tvpto(mto),tveto(mto) dimension wvpto(mto),wspto(mto),wveto(mto),wcpto(mto) common /tov/ tto,zto,pto,rhoto,uto,rto,tpto,teto,tvpto,tveto, + wvpto,wspto,wveto,wcpto c dimension tzo(mto),zzo(mto),pzo(mto),rhozo(mto),uzo(mto),rzo(mto) dimension tpzo(mto),tezo(mto),tvpzo(mto),tvezo(mto) dimension wvpzo(mto),wspzo(mto),wvezo(mto),wcpzo(mto) common /zov/ tzo,zzo,pzo,rhozo,uzo,rzo,tpzo,tezo,tvpzo,tvezo, + wvpzo,wspzo,wvezo,wcpzo c c Local variables dimension tdezo(mto),tdpzo(mto) c c Calculate output z points nzo = 1 + ifix(0.001 + dzmax/dzo) do k=1,nzo zzo(k) = dzo*float(k-1) enddo c c Initialize output variables to missing vmiss = 0.0 do k=1,nzo tzo(k) = vmiss pzo(k) = vmiss rhozo(k) = vmiss uzo(k) = vmiss rzo(k) = vmiss tpzo(k) = vmiss tezo(k) = vmiss tvpzo(k) = vmiss tvezo(k) = vmiss wvpzo(k) = vmiss wspzo(k) = vmiss wvezo(k) = vmiss wcpzo(k) = vmiss enddo c c Interpolate time values to z values do k=1,nzo z = zzo(k) c call ttoz(uto,zto,ntmax,z,fz,ierr) if (ierr .eq. 0) uzo(k) = fz c call ttoz(pto,zto,ntmax,z,fz,ierr) if (ierr .eq. 0) pzo(k) = fz c call ttoz(tpto,zto,ntmax,z,fz,ierr) if (ierr .eq. 0) tpzo(k) = fz - tck c call ttoz(teto,zto,ntmax,z,fz,ierr) if (ierr .eq. 0) tezo(k) = fz - tck c call ttoz(tvpto,zto,ntmax,z,fz,ierr) if (ierr .eq. 0) tvpzo(k) = fz - tck c call ttoz(tveto,zto,ntmax,z,fz,ierr) if (ierr .eq. 0) tvezo(k) = fz - tck c call ttoz(rto,zto,ntmax,z,fz,ierr) if (ierr .eq. 0) rzo(k) = fz c call ttoz(rhoto,zto,ntmax,z,fz,ierr) if (ierr .eq. 0) rhozo(k) = fz c call ttoz(wvpto,zto,ntmax,z,fz,ierr) if (ierr .eq. 0) wvpzo(k) = fz c call ttoz(wveto,zto,ntmax,z,fz,ierr) if (ierr .eq. 0) wvezo(k) = fz c call ttoz(wcpto,zto,ntmax,z,fz,ierr) if (ierr .eq. 0) wcpzo(k) = fz c enddo c c Calculate dewpoint temperature of environment do k=1,nzo ttmp = tezo(k)+tck ptmp = pzo(k) wtmp = wvezo(k) c call tdcal(ptmp,ttmp,wtmp,tdtmp) c tdezo(k) = tdtmp-273.15 enddo c c Calculate dewpoint temperature of parcel do k=1,nzo ttmp = tpzo(k)+tck ptmp = pzo(k) wtmp = wvpzo(k) c call tdcal(ptmp,ttmp,wtmp,tdtmp) c tdpzo(k) = tdtmp-273.15 enddo c if (iprtx .ne. 0) then c Open the output file open(file='lcmod.dat',unit=ludat,form='formatted', + status='replace') c c Write output file write(ludat,100) 100 format(' z P RHOp u Rp ', + ' Tp Te Tvp Tve Tdp Tde wvp wve', + ' wcp') c do k=nzo,1,-1 write(ludat,110) zzo(k)/1000.0,pzo(k)/100.0,rhozo(k),uzo(k), + rzo(k),tpzo(k),tezo(k),tvpzo(k),tvezo(k), + tdpzo(k),tdezo(k),wvpzo(k)*1000.0, + wvezo(k)*1000.0,wcpzo(k)*1000.0 110 format(f6.2,1x,f6.1,1x,f6.3,1x,8(f6.1,1x),3(f6.2,1x)) enddo c close(ludat) endif c return end subroutine ttoz(ft,zt,ntmax,z,fz,ierr) c This routine finds a function value (fz) at the given value of c z by linearly interpolating between the points in the function ft. c The two values to interpolate from are determined from the c z as a function of t (zt). c c If z is outside of the range of zt, c fz is set to zero and ierr=1. c dimension ft(ntmax),zt(ntmax) c c Search for the interval of ft in which z lies. do k=1,ntmax-1 if (z .ge. zt(k) .and. z .le. zt(k+1)) then k1 = k k2 = k+1 exit endif if (k .eq. ntmax-1) then c If you get to here the z value was outside the range of zt. ierr = 1 return endif enddo c ierr = 0 c z1 = zt(k1) z2 = zt(k2) f1 = ft(k1) f2 = ft(k2) if (z1 .eq. z2) then fz = f1 else w1 = (z2-z)/(z2-z1) w2 = (z-z1)/(z2-z1) fz = w1*f1 + w2*f2 endif c return end subroutine avmf(nzo,zzo,uzo,rhozo,tvp,tve,vvavg,vmflux,cape,cin) c This routine calculates the vertically averaged vertical mass flux. c The average vertical velocity (vvavg), CAPE and CIN are also calculated. c dimension zzo(nzo),uzo(nzo),rhozo(nzo) dimension tvp(nzo),tve(nzo) c common /pncon/ pi,g,ra,cpa,cva,rv,cpv,cvv,tck,tf,dtf,ce,cd,alpha c vmflux = 0.0 vvavg = 0.0 do i=1,nzo vmflux = vmflux + rhozo(i)*uzo(i) vvavg = vvavg + uzo(i) enddo c if (nzo .ne. 0) then vmflux = vmflux/float(nzo) vvavg = vvavg/float(nzo) endif c cape = 0.0 cin = 0.0 dz = zzo(2)-zzo(1) c icape=0 icin =0 c do i=2,nzo bterm = g*(tvp(i)-tve(i))/(tve(i)+tck) c if (bterm .gt. 0.0 .and. icape .eq. 0) icape = 1 if (bterm .lt. 0.0 .and. icin .eq. 0) icin = 1 c if (icape .eq. 1) then cape = cape + dz*bterm endif c if (icin .eq. 1) then cin = cin + dz*abs(bterm) endif c if (bterm .lt. 0.0 .and. icape .eq. 1) icape= -1 if (bterm .gt. 0.0 .and. icin .eq. 1) icin = -1 c c write(6,*) 'z,dt: ',zzo(i),tvp(i)-tve(i) enddo c return end subroutine tdcal(p,t,w,td) c This routine calculated the dew point temperature (td in K) c given the pressure (p in Pa), temperature (t in K) and c mixing ratio (w in kg/kg). c c Physical and numerical constants common /pncon/ pi,g,ra,cpa,cva,rv,cpv,cvv,tck,tf,dtf,ce,cd,alpha c dt = 0.05 tmin = -100.0+tck td = tmin c do i=-1,2000 t1 = t-dt*float(i) t2 = t1-dt c if (t1 .le. tmin) return c call escal(t1,es,des,d2es,rl) ws1 = ra*es/(rv*(p-es)) c call escal(t2,es,des,d2es,rl) ws2 = ra*es/(rv*(p-es)) c if (ws1 .ge. w .and. ws2 .lt. w) then td = t1 return endif enddo c return end