Program DIAGNOS
      parameter (l = 21,m = 13,np = 7 ,k = 1)
      real ug(l,m,np),vg(l,m,np),tg(l,m,np)
      real u(l,m,k),v(l,m,k),t(l,m,k)
      real zonav(m,k),areav (k),devzon(l,m,k),devaer(m,k)
      data gleft,gright,glat1,glat2/-25.,25.,-15.,15./
      data p ,del/500.,2.5/
c
c  Open input files .
c

      open (20,file='uv21.dat  ',status='old')
      open (30,file='temp21.dat',status='old')
c
c  Read from multilevel input data .
c
  15  format(10f8.2)
c
      do 13100 ip = 1, np
c
         read (20,15) ((ug(i,j,ip),i=1,l),j=1,m)
         read (20,15) ((vg(i,j,ip),i=1,l),j=1,m)
         read (30,15) ((tg(i,j,ip),i=1,l),j=1,m)
c         
13100 continue
c
c  Seletc data for 500 mb level for this example .
c  Data in multilevel files are arranged as : 
c  1000-850-700-500-300-200-100 mb
c
      do 13102 j  = 1 , m
      do 13102 i  = 1 , l
c      
         u(i,j,k) =  ug (i,j,4)  
         v(i,j,k) =  vg (i,j,4)  
         t(i,j,k) =  tg (i,j,4)
c           
13102 continue
c
c  Compute the zonal and area average of the
c  temperature .
c 
c changed l,m,k to parameter statement inside routine averag1 (glenn)
       call AVERAG1 (t,zonav,areav,gleft,gright,glat1,glat2,del)
c      call AVERAG1 (t,zonav,areav,l,m,k,gleft,gright,glat1,glat2,del)
c
c  Compute the deviation from the zonal average .
c
c dropped l,m,k to use parameter statement
      call DEVIAT  (t,zonav,areav,devzon,devaer) 
c
c  Compute the eddy available potential energy
c  and the zonal available potential energy .
c
c dropped l,m,k to use parameter statement
      call APEZANDE (t,p,az,as,gleft,gright,glat1,glat2,del)
c
c  Compute the eddy kinetic  energy and the 
c  zonal kinetic energy .Conversion of wind 
c  speed inot cgs units is done first .
c
      do 13104 j  = 1,  m
      do 13104 i  = 1,  l
c      
         u(i,j,1) = u(i,j,k)* 100.
         v(i,j,1) = v(i,j,k)* 100.
c      
13104 continue

                                                           
c dropped l,m,k as told above
      call KEZANDE (u,v,ekz,eke,gleft,gright,glat1,glat2,del)
c      
1000  format(21f8.2)
      stop
      end
c
      subroutine DEVIAT (a,zonav,aerav,devzon,devaer) 
c
c  This subroutine computes the departure of 3-D field 'a' 
c  from the zonal average and the departure of the 
c  zonal average of 'a' from the area average.
c
c  Inputs to this routine :
c
c  a        : the 3-d input array of size l x m x kk
c  zonav    : the zonal average of the array a. 
c             zonav is an array of size  m x kk.
c  aerav    : the area average of the array a. aerav is 
c             an array of size kk.
c  l        : the array size in the zonal direction.
c  m        : the array size in the meridional direction.
c  kk       : the array size in the vertical direction.
c
c  Outputs of the routine :
c
c  devzon   : the departure of the array a from the zonal average.
c             devzon is also a 3-d array of size l x m x kk.
c  devaer   : the departure of the zonal average of a from the
c             area average of a. devaer is a 2-d array of size m x kk.
c  NOTE     : all variables in this routine should be in cgs units.
c
      parameter (l = 21,m = 13,np = 7 ,kk = 1)
c
      dimension a(l,m,kk), zonav(m,kk), aerav(kk),
     &          devzon(l,m,kk), devaer(m,kk)
c 
      do 13100  k = 1, kk
      do 13102  j = 1, m 
      do 13104  i = 1, l
c
c  Compute term devzon .
c
         devzon (i,j,k) = a (i,j,k) - zonav (j,k)

13104 continue
c
c  Compute term devaer .
c
         devaer (j,k)   = zonav (j,k) - aerav (k)
13102 continue 
13100 continue
      return
      end
      subroutine AVERAG1 (a,zonav,aerav,
     &            gleft,gright,glat1,glat2,del)
c
c  Given a 3-d field 'a', this subroutine computes the
c  zonal average of 'a' as well as its area average .
c
c  Inputs to the routine :
c
c  a        : the input array of size l x m x kk
c  l        : the array size in the zonal direction.
c  m        : the array size in the meridional direction.
c  kk       : the array size in the vertical direction.
c  gleft    : the western most longitude. (in degrees)
c  gright   : the eastern most longitude. (in degrees)
c  glat1    : the southern most latitude. (in degrees)
c  glat2    : the northern most latitude. (in degrees)
c  del      : the grid spacing in the zonal and the
c             meridional directions (in degrees).
c
c  Outputs of the routine :
c
c  zonav    : the zonal average of the array a. zonav is
c             an array of size m x kk.
c  aerav    : the area average of the array a. areav is
c             an array of size kk.
c  note     : all variables in this routine should be
c             in cgs units.
c
      parameter (l = 21,m = 13,np = 7 ,kk = 1)
c
      dimension a(l,m,kk), zonav(m,kk), aerav(kk)
      pi          = 4.0*atan (1.0)
      r           = 180./pi
c
c  Compute the zonal average
c
      do 13110  k = 1, kk
      do 13111  j = 1, m
         b        = 0.0
      do 13112 i  = 1, l
cb         if (i.eq.1.or.i.eq.l) b = b+del*a(i,j,k)/2.
         b        = b + del*a(i,j,k)
13112 continue
      zonav(j,k)  = b/(gright-gleft)
13111 continue
c
c  Compute the area average
c
      c           = 0.0
      theta       = glat1/r
      do 13113  j = 1, m
         c        = c + cos(theta)*(del/r)*zonav(j,k)
         theta    = theta + del/r
13113 continue
      aerav(k)    = c/(sin(glat2/r) - sin(glat1/r))
13110 continue
      return
      end
      subroutine APEZANDE (temp,p,az,ae,gleft,gright,glat1,
     &                     glat2,del)
c
c  This subroutine computes the eddy available potential
c  and the zonal available potential energy 
c
c  Inputs to this routine :
c
c  temp   : the 3-d temperature field (l x m x kk ).
c  p      : the kk number of pressure levels.
c  gleft  : the western most longitude. (in degrees)
c  gright : the eastern most longitude. (in degrees)
c  glat1  : the southern most latitude. (in degrees)
c  glat2  : the northern most latitude. (in degrees)
c  del    : the grid spacing in the zonal and the 
c           meridional directions (in degrees).
c  l      : the array size in the zonal direction.
c  m      : the array size in the meridional direction.
c  kk     : the array size in the vertical direction.
c
c  Output to this routine :
c
c  az     : the zonal available potential energy.
c  ae     : the eddy  available potential energy.
c
c  NOTE   : all variables in this routine should be in cgs 
c           units except pressure which is in mb.
c
      parameter (l = 21,m = 13,np = 7 ,kk = 1)
c      
      dimension temp(l,m,kk)  , static(kk)  , tza(m,kk),taa(kk)
      dimension devtza(l,m,kk), devtaa(m,kk), dum(l,m,kk)
      dimension p(kk),aape(kk), zake(m,kk)  , avaer(kk)

      pi          = 4.0*atan(1.0)
      r           = pi/180.
c
c  Compute averages and deviations for the temperature field.
c
      call AVERAG1 (temp,tza,taa,gleft,gright,glat1,glat2,del)
      call DEVIAT  (temp,tza,taa,devtza,devtaa) 
c
c  Compute the eddy available potential energy in g/sec**2=erg/cm*2
c  dum represents the following :
c  dum  = (tprime**2)

      do 13110  k = 1, kk
      do 13112  j = 1, m 
      do 13114  i = 1, l
       dum(i,j,k) = devtza(i,j,k)*devtza(i,j,k) 
13114 continue
13112 continue 
13110 continue
c
c  aape is the area average of dum 
c
      call AVERAG1 (dum,zake,aape,gleft,gright,glat1,glat2,del)
c
c  Compute the static stability parameter.
c
      call STASTB  (taa,p,static)
c
c  Divide aape by (2.0*static stability)
c
      do 13116  k = 1, kk
         aape(k)  = (aape(k)/(2.0*static(k)))
13116  continue
c
c  Integrate vertically the term aape to obtain the 
c  eddy available potential energy.

      call INTEGR (aape,ae) 
c
c  ae is the eddy available potential energy in units of
c  g/sec**2 = erg/cm**2  when units are corrected .
c
c  Display results
c
      write(6,1004) ae 
1004  format(2x, 'The eddy available potential energy is ', e10.3,
     &   ' erg/cm**2')
c
c  Compute the zonal available potential energy in g/sec**2=erg/cm*2
c
c  Compute 'tstar' term. the tstar term 
c  is denoted by zake in this routine.
c
c  zake = (zonal average of temp.- the area average of temp.)**2.
c
      do 13118  k = 1, kk
      do 13118  j = 1, m
        zake(j,k) = tza(j,k) - taa(k)
        zake(j,k) = zake(j,k)*zake(j,k)
13118 continue
c
c  Compute the area average avaer of the zake term.
c
      call AVMERID (zake,avaer,glat1,glat2,del)
c
c  Divide avaer by (2.0*static stability)
c
      do 13120  k = 1, kk
13120    avaer(k) = (avaer(k)/(2.0*static(k)))
c
c  Integrate vertically the term avaer to obtain the zonal
c  available potential energy.
c
      call INTEGR (avaer,az) 
c
c  az is the zonal available potential energy in units of
c  g/sec**2 = erg/cm**2  when units are corrected.
c
c  Display  results
c
      write(6,1006) az 
1006  format(2x,'The zonal available potential energy is ', e10.3,
     &   ' erg/cm**2')
      return
      end
      subroutine KEZANDE (uspd, vspd, ekz, eke, gleft, gright,
     &                    glat1,glat2,del)
c
c  This subroutine calculates the zonal kinetic energy and
c  the eddy kinetic energy.
c
c  Inputs to this routine :
c
c  uspd   : the x-component of the wind. a 3-d array of size l x m x kk
c  vspd   : the y-component of the wind. a 3-d array of size l x m x kk
c  gleft  : the wester most longitude.  (degrees)
c  gright : the easter most longitude.  (degrees)
c  glat1  : the southern most latitude. (in degrees)
c  glat2  : the northern most latitude. (in degrees)
c  del    : the grid spacing in the zonal and the meridional directions.
c           (in degrees).
c  l      : the array size in the zonal direction.
c  m      : the array size in the meridional direction.
c  kk     : the array size in the vertical direction.
c
c  Output of this routine :
c
c  ekz    : the zonal kinetic energy   
c  eke    : the eddy  kinetic energy
c
c  Note   : all variables in this routine should be in cgs units.
c  g      : is the acceleration due to gravity at the surface.
c
      parameter (l = 21,m = 13,np = 7 ,kk = 1)
c      
      dimension uspd  (l,m,kk), vspd(l,m,kk), devuza(l,m,kk)
      dimension devvza(l,m,kk), dum (l,m,kk), zake(m,kk),aake(kk)
      dimension h(m,kk), uza(m,kk), vza(m,kk), aad(kk), uaa(kk)
      dimension devuaa(m,kk), devvaa(m,kk), vaa(kk)
c
      g            = 980.6
c
c  Compute the zonal and the area average of uspd.
c
      call AVERAG1 (uspd,uza,uaa,gleft,gright,glat1,glat2,del)
c
c  Compute the deviation of uspd from the zonal average uza, this is
c  represented by the term devuza. Furthermore calculate the deviation
c  of the zonal average uza from the area average uaa, this term is
c  represented by devuaa.
c
      call DEVIAT (uspd,uza,uaa,devuza,devuaa) 
c
c  Compute the zonal average and the area average of vspd
c
      call AVERAG1 (vspd,vza,vaa,gleft,gright,glat1,glat2,del)
c
c  Compute the deviation of vspd from the zonal average vza, this is
c  represented by the term devvza. Compute the deviation of the zonal
c  average vza from the area average vaa, this term is devvaa.
c
      call DEVIAT (vspd,vza,vaa,devvza,devvaa) 
c
c  Compute the eddy kinetic energy .
c  dum represents the following term.
c  dum  =  (uprime)**2 + (vprime)**2
c
      do 13110 k   = 1, kk
      do 13112 j   = 1, m 
      do 13114 i   = 1, l
         dum(i,j,k)= devuza(i,j,k)*devuza(i,j,k) +
     &               devvza(i,j,k)*devvza(i,j,k)
13114 continue
13112 continue 
13110 continue
c
c  Compute the zonal average and the area average of dum.
c  zake is the zonal average of dum.
c  aake is the area average of dum.
c
      call AVERAG1 (dum,zake,aake,gleft,gright,glat1,glat2,del)
c
c  Divide the area average by 2.0 *g.
c
      do 13116 k   = 1, kk
13116    aake(k)   = aake(k)/(2.0*g)  
c
c  Array aake contains  the area average .
c  Integrate aake vertically.
c
      call INTEGR (aake,eke)
c
c  eke is the eddy kinetic energy 
c
      write(6,1000) eke 
1000  format(2x,'The eddy kinetic energy is ', e10.3, '  erg/cm**2') 
c
c  Compute the zonal kinetic energy.
c
      do 13118  k = 1, kk
      do 13120  j = 1, m 
         h(j,k)   = uza(j,k)*uza(j,k) + vza(j,k)*vza(j,k) 
13120 continue 
13118 continue
c
c  Compute the meridional average of h.
c
      call AVMERID (h,aad,glat1,glat2,del) 
c
c  Divide the area average by 2.0*g.
c
      do 13122  k = 1, kk 
         aad(k)   = aad(k)/(2.0*g)
13122 continue 
c
c  Integrate ,aad,vertically.
c
      call INTEGR (aad, ekz ) 
c
c  Display ekz 
c
      write(6,1002) ekz
1002  format(2x,'The zonal kinetic energy is ', e10.3, ' erg/cm**2')
      return
      end
      subroutine AVMERID (f,avaer,glat1,glat2,del) 
c
c  This subroutine computes the meridional mean of a 2-d array f.
c  f is a 2-d array of dimension (m,kk).
c
c  Inputs to this routine :
c
c  f        : the 2-d array of size m x kk.
c  glat1    : the southern most latitude (in degrees).
c  glat2    : the northern most latitude (in degrees).
c  del      : the grid spacing in the zonal and the 
c             meridional directions (in degrees).
c  m        : the array size in the meridional direction.
c  kk       : the array size in the vertical   direction.
c
c  Output of this routine :

c  avaer    : the meridional average of the field f. This is a 
c             1-d array of size kk.
c
c  NOTE     : all variables in this routine should be in cgs units.
c
      parameter (l = 21,m = 13,np = 7 ,kk = 1)
c      
      dimension f(m,kk),avaer(kk)
      pi          = 4.0*atan(1.0)
      r           = 180./pi
      do 13100  k = 1, kk 
         theta    = glat1/r
         g        = 0.0
      do 13102  j = 1, m 
         g        = g + cos(theta)*(del/r)*f(j,k)
         theta    = theta+(del/r) 
13102 continue 
      avaer(k)    = g/(sin(glat2/r)-sin(glat1/r))
13100 continue 
      return
      end
      subroutine INTEGR(a,sum)
c
c  This subroutine performs vertical integration .The vertical 
c  coordinate is pressure (in mb). The pressure levels are separated
c  by 100 mb. In the case the pressure levels are uneven ,one has to
c  make necessary changes before integration.
c
c  Inputs to this routine :
c  a        : the field to be integrated .a is a 1-d array of size kk.
c  kk       : is the dimension of the array a. this is also the number of
c             vertical levels.
c
c  Output of this routine :
c
c  sum      : represents the vertical integral of the field a. 
c  NOTE     : all variables in this routine should be in cgs units
c             except pressure, which has units of mb.
c
      parameter (l = 21,m = 13,np = 7 ,kk = 1)
c
      dimension a(kk)
      sum         = 0.0
      do 13100  i = 1, kk
         if (i.eq.1) sum = sum + 0.5*a(i)
         if (i.gt.1.and.i.lt.kk) sum = sum + a(i)
         if (i.eq.kk) sum = sum + 0.5*a(i)
13100  continue
c
c  Multiply sum by 1000. to convert mb to cgs units.
c
      sum         = 1000.*sum
      return
      end
      subroutine STASTB (tmpaa,p,static)
c
c  This routine calculates the static stability parameter.
c  The static stability parameter is defined as follows :
c
c  static(k) = gravity*(rdgas*tempareaav(k)/cp*press(k)
c              -d(tempareaav(k))/d(p(k))
c
c  Inputs to this routine :
c
c  tmpaa    : the area average of the temperature field. this is a 1-d 
c             array of size kk.
c  p        : the kk number of pressure levels (bottom to top).
c  kk       : the number of vertical, levels.
c
c  Output of this routine :
c
c  static   : the static stability parameter. this is a 1-d array of size kk.
c
c  other variable used    :
c
c  dt       : this represents the term d(tmpaa)/d(p).this is 
c             also 1-d array of size kk.
c  r        : the universal gas constant
c  g        : the acceleration due to gravity at the surface of earth.
c  cp       : the specific heat at constant pressure.
c
c  NOTE     : all variables in this routine should be in cgs units except
c             pressure, which is in mb.
c
      parameter (l = 21,m = 13,np = 7 ,kk = 1)
c
      dimension tmpaa(kk),p(kk),static(kk),dt(kk)
c
c  static stability in units of (k**2/cm)
c  using    :
c  g        = 980.616  cm/sec**2 
c  r        = 2.8704e6 erg/(g*k)
c cp        = 1.004e7  erg/(g*k)
      cp          = 1.004e7
      g           = 980.616
      r           = 2.8704e6
      t1          = g*r/cp
      do 13110  k = 1, kk 
c      
         if (k.eq.1)   dt(k) = (tmpaa(k+1)-tmpaa(k))/(p(k+1)-p(k))
         if (k.eq.kk-1)dt(k) = (tmpaa(k)-tmpaa(k-1))/(p(k) - p(k-1)) 
         if (k.eq.1) go to 13121 
         if (k.eq.kk-1) go to 13121
c          
         dt(k)    = (tmpaa(k+1)-tmpaa(k-1))/(p(k+1)-p(k-1))
c
c  The if statements are used to keep dt range within the region of
c  data values . Multiplication by 1000.0 changes pressure from mb.
c  to cgs units.

13121 static(k)   = ((t1*tmpaa(k))/(1000.0*p(k))) + (g/1000.0)*dt(k)
13110 continue
      return
      end