c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE ENERGY                     ######
c     ######                                                      ######
c     ######                Copyright (c) 1992                    ######
c     ######  Center for the Analysis and Prediction of Storms    ######
c     ######    University of Oklahoma.  All rights reserved.     ######
c     ######                                                      ######
c     ##################################################################
c     ##################################################################
c
      SUBROUTINE ENERGY(nx,ny,nz,
     :           u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh,rhobar,
     :           x,y,z,zp,hterain, j3, 
     :           rhostr,urho,vrho,wrho,tem4) 
c
c#######################################################################
c
c     PURPOSE:
c
c     Compute the density weighted domain average of grid scale kinetic 
c     energy, momentum, potential temperature and potential temperature 
c     variance. 
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c     11/13/91.   
c
c     MODIFICATION HISTORY:
c
c     6/06/92 (M. Xue)
c     Added full documentation.
C
c#######################################################################
c
c     INPUT:
c
c       nx       Number of grid points in the x-direction (east/west)
c       ny       Number of grid points in the y-direction (north/south)
c       nz       Number of grid points in the vertical
c
c       u        x component of velocity (m/s)
c       v        y component of velocity (m/s)
c       w        Vertical component of Cartesian velocity (m/s)
c       ptprt    Perturbation potential temperature (K)
c       pprt     Perturbation pressure (Pascal)
c       qv       Water vapor specific humidity (kg/kg)
c       qc       Cloud water mixing ratio (kg/kg)
c       qr       Rainwater mixing ratio (kg/kg)
c       qi       Cloud ice mixing ratio (kg/kg)
c       qs       Snow mixing ratio (kg/kg)
c       qh       Hail mixing ratio (kg/kg)
c
c       rhobar   Base state density (kg/m**3)
c
c       x        x coordinate of grid points in physical/comp. space (m)
c       y        y coordinate of grid points in physical/comp. space (m)
c       z        z coordinate of grid points in computational space (m)
c       zp       Vertical coordinate of grid points in physical space (m)
c       hterain  The height of terrain (m)
c
c       j3       Coordinate transform Jacobian  d(zp)/dz
c
c     OUTPUT:
c
c       None.
c
c     WORK ARRAY:
c
c       rhostr   j3*rhobar, a work array.
c       urho     rhobar*u, a work array.
c       vrho     rhobar*v, a work array.
c       wrho     rhobar*w, a work array.
c       tem4     work array.
c
c#######################################################################
c

c
c#######################################################################
c
c     Variable Declarations.
c
c#######################################################################
c
      implicit none 
c
      integer nt           ! The no. of t-levels of t-dependent arrays.
      integer tpast        ! Index of time level for the past time.
      integer tpresent     ! Index of time level for the present time.
      integer tfuture      ! Index of time level for the future time.

      parameter (nt=3, tpast=1, tpresent=2, tfuture=3)

      integer nx,ny,nz          ! Number of grid points in 3 directions
c
      real u     (nx,ny,nz,nt)  ! Total u-velocity (m/s)
      real v     (nx,ny,nz,nt)  ! Total v-velocity (m/s)
      real w     (nx,ny,nz,nt)  ! Total w-velocity (m/s)
      real ptprt (nx,ny,nz,nt)  ! Perturbation potential temperature (K)
      real pprt  (nx,ny,nz,nt)  ! Perturbation pressure (Pascal)
c
      real qv    (nx,ny,nz,nt)  ! Water vapor specific humidity (kg/kg)
      real qc    (nx,ny,nz,nt)  ! Cloud water mixing ratio (kg/kg)
      real qr    (nx,ny,nz,nt)  ! Rain water mixing ratio (kg/kg)
      real qi    (nx,ny,nz,nt)  ! Cloud ice mixing ratio (kg/kg)
      real qs    (nx,ny,nz,nt)  ! Snow mixing ratio (kg/kg)
      real qh    (nx,ny,nz,nt)  ! Hail mixing ratio (kg/kg)
c
      real rhobar(nx,ny,nz)     ! Base state air density (kg/m**3)

      real x     (nx)           ! The x-coord. of the physical and
                                ! computational grid. Defined at u-point.
      real y     (ny)           ! The y-coord. of the physical and
                                ! computational grid. Defined at v-point.
      real z     (nz)           ! The z-coord. of the computational grid.
                                ! Defined at w-point on the staggered grid.
      real zp    (nx,ny,nz)     ! The physical height coordinate defined at
                                ! w-point of the staggered grid.

      real hterain(nx,ny)       ! The height of the terrain.

      real j3    (nx,ny,nz)     ! Coordinate transform Jacobian defined as
                                ! d( zp )/d( z )
c
      real rhostr(nx,ny,nz)     ! Temporary work array
      real urho  (nx,ny,nz)     ! Temporary work array
      real vrho  (nx,ny,nz)     ! Temporary work array
      real wrho  (nx,ny,nz)     ! Temporary work array
      real tem4  (nx,ny,nz)     ! Temporary work array
c
c#######################################################################
c
c     Misc. local variables: 
c
c#######################################################################
c
      integer ncalls
      data ncalls /0/
      save ncalls

      integer istat
      integer tlevel            ! Time level at which data are printed.
      integer i,j,k

      character engfn*80        ! File name of the energy etc results
      integer lengfn            ! String length of the file name

      real tmass     ! Total mass of the air in model domain (kg)
      real keu       ! Contribution of u to the total kinetic energy
      real kev       ! Contribution of v to the total kinetic energy
      real kew       ! Contribution of w to the total kinetic energy
      real ke        ! Domain average kinetic energy per unit mass (m/s)**2
      real ptvari    ! Domain average variance in ptprt per unit mass (K**2) 

      real momntu    ! Domain average x-component of momentum 
                     ! per unit mass (m/s)
      real momntv    ! Domain average y-component of momentum 
                     ! per unit mass (m/s)
      real momntw    ! Domain average z-component of momentum 
                     ! per unit mass (m/s)
      real ptavg     ! Domain average potential temperature perturbation (K)
 
      real dxdz05,dydz05,dxdy05
c
c#######################################################################
c
c     Include files: 
c
c#######################################################################
c
      include 'phycst.inc'
      include 'globcst.inc'
c
c#######################################################################
c
c     Routines called:
c
c#######################################################################
c
      external aamult
      external uvwrho
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c

c
C#######################################################################
c
c     Dump restart data files every nrstout number of time steps. 
c
C#######################################################################
c            
       tlevel=tpresent
c
c#######################################################################
c
c     Calculate urho=j3*rhobar*u, vrho=j3*rhobar*v, wrho=j3*rhobar*w
c
c#######################################################################
c
      CALL aamult(rhobar,j3,nx,ny,nz,1,nx-1,1,ny-1,1,nz-1, rhostr)

      CALL uvwrho(nx,ny,nz,
     :            u(1,1,1,tlevel),v(1,1,1,tlevel),w(1,1,1,tlevel),
     :            rhostr, urho,vrho,wrho, tem4)

c
c#######################################################################
c
c     Calculate the total air mass inside the phyiscal boundaries 
c     based on base state air density.
c
c#######################################################################
c
      tmass = 0.0

      DO 50 k=2,nz-2
      DO 50 j=2,ny-2
      DO 50 i=2,nx-2
        tmass = tmass + rhostr(i,j,k)
     :                *(x(i+1)-x(i))*(y(j+1)-y(j))*(z(k+1)-z(k))
50    CONTINUE
c
c#######################################################################
c
c     Contribution of u to the total kinetic energy: 
c
c#######################################################################
c

      keu = 0.0

      DO 100 k=2,nz-2
      DO 100 j=2,ny-2

        dydz05 = (y(j+1)-y(j))*(z(k+1)-z(k))*0.5

        keu = keu + 0.5*urho(  2, j,k)*u(  2, j,k,tlevel)
     :                *dydz05*(x(3)-x(1))

        keu = keu + 0.5*urho(nx-1,j,k)*u(nx-1,j,k,tlevel)
     :                *dydz05*(x(nx)-x(nx-2))

        DO 110 i=3,nx-2
          keu = keu + urho(i,j,k)*u(i,j,k,tlevel)
     :                *dydz05*(x(i+1)-x(i-1))
110     CONTINUE

100   CONTINUE

c
c#######################################################################
c
c     Contribution of v to the total kinetic energy: 
c
c#######################################################################
c

      kev = 0.0

      DO 200 k=2,nz-2
      DO 200 i=2,nx-2

        dxdz05 = (x(i+1)-x(i))*(z(k+1)-z(k))*0.5

        kev = kev + 0.5*vrho(i, 2,  k)*v(i, 2,  k,tlevel)
     :                *dxdz05*(y(3)-y(1))

        kev = kev + 0.5*vrho(i,ny-1,k)*v(i,ny-1,k,tlevel)
     :                *dxdz05*(y(ny)-y(ny-2))

        DO 210 j=3,ny-2
          kev = kev + vrho(i,j,k)*v(i,j,k,tlevel)
     :                *dxdz05*(y(j+1)-y(j-1))
210     CONTINUE

200   CONTINUE
c
c#######################################################################
c
c     Contribution of w to the total kinetic energy: 
c
c#######################################################################
c
      kew = 0.0

      DO 300 i=2,nx-2
      DO 300 j=2,ny-2

        dxdy05 = (x(i+1)-x(i))*(y(j+1)-y(j))*0.5

        kew = kew + 0.5*wrho(i,j, 2  )*w(i,j, 2  ,tlevel)
     :                *dxdy05*(z(3)-z(1))
        kew = kew + 0.5*wrho(i,j,nz-1)*w(i,j,nz-1,tlevel)
     :                *dxdy05*(z(nz)-z(nz-2))

        DO 310 k=3,nz-2
          kew = kew + wrho(i,j,k)*w(i,j,k,tlevel)
     :                *dxdy05*(z(k+1)-z(k-1))
310     CONTINUE

300   CONTINUE
c
c#######################################################################
c
c     The domain average kinetic energy per unit mass: 
c
c#######################################################################
c
      ke = 0.5*(keu+kev+kew)/tmass
c
c#######################################################################
c
c     The domain average potential temperature variance: 
c
c#######################################################################
c
      ptvari = 0.0

      DO 400 k=2,nz-2
      DO 400 j=2,ny-2
      DO 400 i=2,nx-2

        ptvari = ptvari + rhostr(i,j,k)*ptprt(i,j,k,tlevel)**2
     :                  *(x(i+1)-x(i))*(y(j+1)-y(j))*(z(k+1)-z(k))

400   CONTINUE

      ptvari = ptvari/tmass

c
c#######################################################################
c
c     Calcualte the domain average momentum in x direction. 
c
c#######################################################################
c

      momntu = 0.0

      DO 600 k=2,nz-2
      DO 600 j=2,ny-2

        dydz05 = (y(j+1)-y(j))*(z(k+1)-z(k))*0.5

        momntu = momntu + 0.5*urho(  2, j,k)*dydz05*(x(3)-x(1))

        momntu = momntu + 0.5*urho(nx-1,j,k)*dydz05*(x(nx)-x(nx-2))

        DO 610 i=3,nx-2
          momntu = momntu + urho(i,j,k)*dydz05*(x(i+1)-x(i-1))
610     CONTINUE

600   CONTINUE

      momntu = momntu/tmass
c
c#######################################################################
c
c     Calcualte the domain average momentum in y direction. 
c
c#######################################################################
c

      momntv = 0.0

      DO 700 k=2,nz-2
      DO 700 i=2,nx-2

        dxdz05 = (x(i+1)-x(i))*(z(k+1)-z(k))*0.5

        momntv = momntv + 0.5*vrho(i, 2,  k)*dxdz05*(y(3)-y(1))

        momntv = momntv + 0.5*vrho(i,ny-1,k)*dxdz05*(y(ny)-y(ny-2))

        DO 710 j=3,ny-2
          momntv = momntv + vrho(i,j,k)*dxdz05*(y(j+1)-y(j-1))
710     CONTINUE

700   CONTINUE

      momntv = momntv/tmass
c
c#######################################################################
c
c     Calcualte the domain average momentum in z direction. 
c
c#######################################################################
c
      momntw = 0.0

      DO 800 i=2,nx-2
      DO 800 j=2,ny-2

        dxdy05 = (x(i+1)-x(i))*(y(j+1)-y(j))*0.5

        momntw = momntw + 0.5*wrho(i,j, 2  )*dxdy05*(z(3)-z(1))
        momntw = momntw + 0.5*wrho(i,j,nz-1)*dxdy05*(z(nz)-z(nz-2))

        DO 810 k=3,nz-2
          momntw = momntw + wrho(i,j,k)*dxdy05*(z(k+1)-z(k-1))
810     CONTINUE

800   CONTINUE

      momntw = momntw/tmass

c
c#######################################################################
c
c     The mass-weighted domain average of potential temperature 
c     perturbation (K) (ptprt per unit mass) 
c
c#######################################################################
c
      ptavg = 0.0

      DO 900 k=2,nz-2
      DO 900 j=2,ny-2
      DO 900 i=2,nx-2

        ptavg = ptavg + rhostr(i,j,k)*ptprt(i,j,k,tlevel)
     :                  *(x(i+1)-x(i))*(y(j+1)-y(j))*(z(k+1)-z(k))

900   CONTINUE

      ptavg = ptavg/tmass

c
c#######################################################################
c
c     Write the results into a file 
c
c#######################################################################
c
      IF(ncalls.eq.0) THEN

        IF( dirname.ne.' ' ) THEN

          engfn  = dirname(1:ldirnam)//'/'//runname(1:6)//'.eng'
          lengfn = 10 + ldirnam + 1

        ELSE

          engfn  = runname(1:6)//'.eng'
          lengfn = 10

        ENDIF

        OPEN(ncheng,form='formatted',status='unknown',
     :              file=engfn(1:lengfn),iostat=istat)

        IF( istat.ne.0) THEN

          write(6,'(/a,i2,/a/)')
     :    ' Error occured when opening file '//runname(1:6)//'.eng'//
     :    ' using FORTRAN unit ',ncheng,' Program stopped in ENERGY.'
          STOP

        ENDIF

        write(ncheng,'(a)') ''''//runname//''''

        write(ncheng,'(t4,a,t15,a,t30,a,t45,a,t60,a,t75,a,t90,a)')
     :       '''time','    KE','  u-momentum','  v-momentum',
     :         '  w-momentum','    ptavg','    ptvari'''

      ENDIF

      write(ncheng,'(f9.3,6f15.7)')
     :      curtim,ke,momntu,momntv,momntw,ptvari,ptavg

      ncalls = 1

      RETURN
      END
