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