program KUO c c This program computes cloud related heat source (q1), c moisture sink(q2), and rainfall rate for a vertical c column over a given grid point. c This column model requires input of temperature, c dew point temperature and vertical velocity at n c levels in the vertical .It also requires the 700 mb c relative vorticity at that point . The moistening c parameter,b and the meso-scale convergence, eta are c computed before the main subroutine cvheat is called . c c c Definitions of variables c c dp : pressure increment (mb ) c tp : temperature (K ) c td : dew point temperature (K ) c w : vertical p-velocity (mb/s ) c wbar : average vertical velocity c vort : vorticity at 700 mb (sec-1) c pt : pressure at temperature levels (mb ) c parameter (n = 10) real pt(n),t(n),td(n),w(n) ,q1(n), q2(n) c c Input data. c data pt/1000.,900.,800.,700.,600.,500.,400.,300.,200.,100./ data t/301.46,296.45,290.51,285.24,279.28,270.57,258.77, + 243.39,222.21,189.53/ data td/298.79,294.64,286.67,275.86,268.12,260.73,251.75, + 241.81,220.63,187.96/ data w/-0.16250e-02,-0.17760e-02,-0.18789e-02,-0.20905e-02, + -0.21773e-02,-0.18447e-02,-0.13510e-02,-0.10983e-02, + -0.91707e-03,-0.30293e-03/ c open(2,file ='prof.dat',status='unknown') c vort = 0.81442e-05 dp = 100. c c Define coefficients of multiple regression for c heating and moistening parameters c a = 0.158e05 b = 0.304e03 c = 0.476 p = 0.107e05 q = 0.107e03 r = 0.870 c c Compute the average vertical velocity c wbar = 0. do 7100 k = 1, n wbar = wbar + w(k)/float(n) 7100 continue c c Compute b and eta c rnumer = a*vort+b*wbar+c denom = (a+p)*vort+(b+q)*wbar+(c+r) yb = rnumer/denom eta = denom-1. c c Reset parameter b to [0,1] in case it is out of range. c this can happen when running the code over lagre domain. c the coeffiecients given here were obtained using GATE c domain only . c if (yb.gt.1.0) yb = 1.0 if (yb.lt.0.0) yb = 0.0 c c Call the main routine to compute the heating ,moistening c and rainfall rates. c call CVHEAT (pt,t,td,w,yb,eta,n,dp,q1,q2,rrr) c c Display outpts for the column c do 7102 k = 1, n - 1 print *,q1(k),q2(k),rrr write (2,1000)q1(k),q2(k),rrr 7102 continue 1000 format (2x,3e12.5) stop end subroutine CVHEAT (pt,t,td,w,yb,eta,n,dp,q1,q2,rr) c c This subroutine computes the heating and moistening c contributions due to deep convection ala kuo (1965). c The kuo cloud is defined by its structure function c for the vertical distribution of heating (tc-te) c and moistening (qc-qe) c c c Definitions of some variables c c ts : temperature of a local adiabat c q : specific humidity c qs : saturation specific humidity c w : vertical p-velocity c xqt : moisture supply required to raise the large scale c temperature to cloud temperature c xqq : moisture supply required to increase moisture to c : saturation at the cloud temperature c yat : partitioning factor for potential temperature c yaq : partitioning factor for specific humidity c eta : mesoscale convergence parameter c yb : moistening parameter c cdt : cloud time scale (1200s) c q1 : apparent heat source (deg/day) c q2 : apparent moisture sink (deg/day) c conv : large scale available moisture c lclb : cloud base c lclt : cloud top c qadv : vertical advection of specific humidity c tadv : vertical advection of potential temperature c real pt(13),t(13),td(13),w(13),tcnvt(13) real q(13),qs(13),rh(13),ts(13),pconst(13) real lconv(13),theta(13),cconv(13),qsadv(13) real qadv(13),tadv(13),q1(13),q2(13),ptf(13) c c Define constants . c data rho,grav,eps1,day/1.25,9.81,1.e-02,86400./ data cp,yl,cdt,rcp,eps/1004.5,2.5e06,1200.,0.286,0.622/ n1 = n-1 gasr = 2000./7. cpinv = 1 ./cp cdtinv = 1./cdt fact = dp*100./grav c c factor for vertical advection of theta at level 1 based c on the gradient between level 1 and level 2 . c cabcd = 1.0 c c Critical relative humidity for convection at cloud base. c rhcrit = .80 c c Critical relative humidity factor for stable heating. c rhc = 1. c c Constants for heating calculation c cd1 = 6.11*0.622 cd2 = 25.22 cd3 = 273.16 cd4 = 5.33 cd5 = cd2*cd3 cd6 = eps*yl/gasr cd7 = cp/yl c c Initialize the heating terms to zero. c do 7110 k = 1, n1 q1(k) = 0. 7110 q2(k) = 0. c do 7111 k = 1, n c c Conversion factor from theta to temperature . c tcnvt(k) = (pt(k)/1000.)**rcp c c Coefficient of diabatic heating. c pconst(k) = cpinv/tcnvt(k) theta(k) = t(k)/tcnvt(k) 7111 continue c c Compute q, qs and rh . c do 7112 k = 1, n call QFRMTD (pt(k),t(k),t(k),qs(k)) call QFRMTD (pt(k),td(k),td(k),q(k)) rh(k) = q(k)/qs(k) 7112 continue c c Check stability. c c lconv.....1 unstable c lconv.....0 stable c lconv....-1 downward motion c do 7113 k = 2, n if (w(k)) 7532,7531,7531 7532 xx = w(k)/dp qsadv(k) = (qs(k-1)-qs(k))*xx qadv(k) = (q(k-1)-q(k))*xx tadv(k) = (theta(k-1)-theta(k))*xx x = pconst(k)*yl*qsadv(k)+tadv(k) if (x.ge.0) lconv(k) = 0 if (x.lt.0) lconv(k) = 1 go to 7113 7531 lconv(k) = -1 7113 continue x = w(1)/dp qsadv(1) = x*(qs(1)-qs(2))*0.5*cabcd c c Find cloud base. c lclb = 100 do 7114 k = 1, n1 if (lconv(k).ne.1) go to 7114 plcl = pt(k) if (rh(k).lt.rhcrit)go to 7114 tlcl = t(k) qlcl = q(k) qslcl = qs(k) rhlcl = qlcl/qslcl if (rhlcl.lt.rhcrit) go to 7114 lclb = k go to 7541 7114 continue 7541 continue if (lclb.ge.n1) go to 400 c c Compute local moist adiabat. c do 7115 k = 1, n ptf(k) = pt(k)**(-rcp) 7115 continue do 7116 k = 1, n1 cconv(k) = yl*cpinv*rhc*ptf(k)*cd1/pt(k) ts(k) = t(k) 7116 continue lclbp = lclb+1 do 7117 k = lclbp, n1 cx1 = cconv(k) cx2 = -ts(k-1)*ptf(k-1)-yl*cpinv*ptf(k)*qs(k-1) c c ts is obtained through an iterative process c using the Newton Raphson's method c assuming a constant lapse rate as first guess. c ts(k) = ts(k-1)-dp*0.06 gts = ts(k) do 7118 kk= 1, 30 xx = 1./gts x = xx*cd3 ft1 = cx1*exp(cd2*(1.-x))*x**cd4 ft2 = gts*ptf(k)+cx2 ft = ft1+ft2 dftdt = ft1*xx*(cd5/gts-cd4)+ptf(k) gts = gts-ft/dftdt if (abs(gts-ts(k)).lt.eps1) go to 7581 ts(k) = gts 7118 continue print 5001 stop 7581 continue ts(k) = gts x = cd3/ts(k) qs(k) = cd1*exp(cd2*(1.-x))*x**cd4/pt(k)*rhc 7117 continue c c Find cloud top. c lclt = lclb do 7119 k = lclbp, n1 if (ts(k).lt.t(k)) go to 7598 lclt = k 7119 continue 7598 if (lclt.le.lclb+1) go to 400 c c Compute available moisture supply c conv = 0. delp = (pt(lclb)-pt(lclt))*100. do 7120 k = lclb,lclt x = -qadv(k) if (x.le.0.) go to 7120 lconv(k) = 999 conv = conv+x 7120 continue if (conv.le.0.) go to 99 convi = conv convm = conv*eta c c Compute the required moistening and heating for cloud. c xqt = 0. xqq = 0. xqcool = 0. xxx = cd7*cdtinv do 7121 k = lclb, lclt x = cd7*tcnvt(k)*tadv(k) if (k.eq.lclb) x = 0. dtemp = ts(k)-t(k) dq = qs(k)-q(k) xy = xxx*dtemp+x if (xy.lt.0.) go to 7591 xqt = xqt+xy 7591 continue xqq = dq+xqq xqcool = xqcool-x 7121 continue if (xqcool.ge.0.) go to 400 xqq = xqq*cdtinv*fact if (xqt) 7592,7592,7593 7592 continue yat = 0. if (xqq.le.0.) go to 400 yaq = (conv*(1.+eta)*yb)/xqq go to 7594 7593 continue if (xqq) 7595,7595,7596 7595 continue yaq = 0. yat = (conv*(1.+eta)*(1.-yb))/xqt go to 7594 7596 continue c c Compute the fractional area of a grid square c for both heating and moistening . c yat = (conv*(1.+eta)*(1.-yb))/xqt yaq = (conv*(1.+eta)*yb)/xqq 7594 continue x1 = yaq*cdtinv do 7122 k = lclb,l clt dtemp = ts(k)-t(k) x = yat*dtemp*cdtinv/tcnvt(k) dq = qs(k)-q(k) q1(k) = (x+tadv(k)*yat)*day*tcnvt(k) if (k.eq.lclb) q1(k) = 0. q2(k) = -(x1*dq+qadv(k))*day/cd7 7122 continue go to 300 400 continue do 7123 k = lclb, n1 if (lconv(k).eq.999) lconv(k) = 0 7123 continue 300 continue go to 1001 99 print 1002 1001 continue c c Sum q1 for rainfall rate c q1sum = 0. do 7124 k = 1, n 7124 q1sum = q1sum + q1(k) factr = (1000./981.)*(cp/yl)*dp*10. rr = q1sum*factr return 120 format(1x,'lconv',2x,13i4) 190 format(1x,'pt',1x,f6.0,2x,'q1',2x,e12.5,2x,'q2',2x,e12.5) 1002 format(10x,/,'integrated moisture convergence less + than zero') 5001 format(10x,/,'no convergence on moist adiabat') end subroutine QFRMTD (p,t,td,q) c c This subroutine returns specific humidity values from an c input of temperature, dew point values and pressure. c c input parameter : c p : pressure in millibar c t : temperature in degrees kelvin c td : dew point temperature in degrees kelvin c output parameter : c q : specific humidity in gm/gm c a = 17.26 b = 35.86 if (td.le.263.00) a = 21.87 if (td.le.263.00) b = 7.66 e = 6.11*exp(a*(td-273.16)/(td-b)) q = 0.622*e/(p-0.378*e) return end