PROGRAM arpsenscv,110
!
!##################################################################
!##################################################################
!######                                                      ######
!######                  PROGRAM ARPSENSCV                   ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!
!
!-----------------------------------------------------------------------
!  PURPOSE:
!  Read in a series of ARPS history files, generate 2-D pruducts data
!   for display. The output may be in the same grid as the input
!   data, or can be different from it.
!  The output data is on a qni x qnj grid with its southwest corner
!   at qswcorn(lat degree),qswcorw(lon degree) specified in enscv.inc,
!   where qnecorn and qnecorw are also set. However, the mapproj,
!   trulon, trulat(1),trulat(2) parameters should be specified in this
!   program. Please check the NEWMAP subroutine.
!
!  Author :    Dingchen Hou
!  History: Apr. 30 1998, developmed from the framework of ARPSPLT.
!        Sep. 15, 1999, modified to include perturbations in soil
!          variables.
!        Feb 2002 (F. KONG), modified to read nx,ny,nz directly from
!          input file; add 'iaccu' into output_data namelist; & fix bugs
!-----------------------------------------------------------------------
!
!  DATA ARRAYS READ IN:
!
!    x        x-coordinate of grid points in physical/comp. space (m)
!    y        y-coordinate of grid points in physical/comp. space (m)
!    z        z-coordinate of grid points in computational space (km)
!    zp       z-coordinate of grid points in computational space (m)
!
!    uprt     x-component of perturbation velocity (m/s)
!    vprt     y-component of perturbation velocity (m/s)
!    wprt     vertical component of perturbation velocity in Cartesian
!             coordinates (m/s).
!
!    ptprt    perturbation potential temperature (K)
!    pprt     perturbation pressure (Pascal)
!
!    qvprt    perturbation water vapor mixing ratio (kg/kg)
!    qc       Cloud water mixing ratio (kg/kg)
!    qr       Rainwater mixing ratio (kg/kg)
!    qi       Cloud ice mixing ratio (kg/kg)
!    qs       Snow mixing ratio (kg/kg)
!    qh       Hail mixing ratio (kg/kg)
!    tke      Turbulent Kinetic Energy ((m/s)**2)
!    kmh      Horizontal turb. mixing coef. for momentum ( m**2/s )
!    kmv      Vertical turb. mixing coef. for momentum ( m**2/s )
!
!    ubar     Base state x-velocity component (m/s)
!    vbar     Base state y-velocity component (m/s)
!    wbar     Base state z-velocity component (m/s)
!    ptbar    Base state potential temperature (K)
!    pbar     Base state pressure (Pascal)
!    rhobar   Base state air density (kg/m**3)
!    qvbar    Base state water vapor mixing ratio (kg/kg)
!
!    soiltyp  Soil type
!    stypfrct Soil type fraction
!    vegtyp   Vegetation type
!    lai      Leaf Area Index
!    roufns   Surface roughness
!    veg      Vegetation fraction
!
!    tsfc     soil skin temperature (K)
!    tsoil    Deep soil temperature (K)
!    wetsfc   Surface soil moisture
!    wetdp    Deep soil moisture
!    wetcanp  Canopy water amount
!    raing    Grid supersaturation rain
!    rainc    Cumulus convective rain(mm)
!    raing    Total rain (rainc+raing)(mm)
!
!    psl      Sea level pressure (mb)
!
!  CALCULATED DATA ARRAYS:
!
!    u        x-component of velocity (m/s)
!    v        y-component of velocity (m/s)
!    w        z-component of velocity (m/s)
!    pt       potential temperature (K)
!    qv       water vapor mixing ratio (kg/kg)
!    td       dew-point temperature (C)
!    cape1    CAPE  (J/kg)
!    cin1     CIN   (J/kg)
!    thet     theta_E (K)
!    heli     helicity (m2/s2)
!    srlfl    storm-relative low-level flow (0-2km AGL)
!    srmfl    storm-relative mid-level flow (2-9km AGL)
!    shr37    7km - 3km wind shear
!    ustrm    Estimated storm motion (Bob Johns)
!    vstrm    Estimated storm motion (Bob Johns)
!    capst    CAPE strength
!    blcon    boundary layer convergence
!    ct       convective temperature
!
!  WORK ARRAYS:
!
!    tem1     Temporary work array.
!    tem2     Temporary work array.
!    tem3     Temporary work array.
!    tem4     Temporary work array.
!    tem5     Temporary work array.
!    tem6     Temporary work array.
!    tem7     Temporary work array.
!    tem8     Temporary work array.
!    tem9     Temporary work array.
!
!  (These arrays are defined and used locally (i.e. inside this
!   subroutine), they may also be passed into routines called by
!   this one. Exiting the call to this subroutine, these temporary
!   work arrays may be used for other purposes, and therefore their
!   contents may be overwritten. Please examine the usage of work
!   arrays before you make any change to the code.)
!
!-----------------------------------------------------------------------
!
!  Arrays for plots on constant pressure levels
!
!-----------------------------------------------------------------------
!
!   hgt      z-coor of scalar point in physical space (10m)
!   t700     Temperature (K) on 700mb pressure grids

!   zps3d    negative log pressure(Pascal) at ARPS grid points
!   algpzc   -log(pressure) at scalar grid points
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
 
!-----------------------------------------------------------------------
!  Include files:
!-----------------------------------------------------------------------

  INCLUDE 'indtflg.inc'
  INCLUDE 'globcst.inc'
  INCLUDE 'grid.inc'
  INCLUDE 'phycst.inc'
  INCLUDE 'enscv.inc'
!
!-----------------------------------------------------------------------
!
!  Some universal constants
!
!-----------------------------------------------------------------------
!
  REAL*4        kappa,gamma,ex1,ex2,t00,p00,mbtopa
  PARAMETER     (kappa=287.053/cp,                                      &
                gamma=6.5,              & ! 6.5 K/km
                ex1=0.1903643,          & ! R*gamma/g
                ex2=5.2558774,          & ! g/R/gamma
                mbtopa=100.)

!-----------------------------------------------------------------------
!
!  Define the dimensions nx, ny, nz
!
!-----------------------------------------------------------------------
  INTEGER :: nx, ny, nz

!-----------------------------------------------------------------------
!
!  Arrays to be read in:
!
!-----------------------------------------------------------------------
!
  REAL, ALLOCATABLE :: x     (:)      ! The x-coord. of the physical and
                                      ! computational grid. Defined at u-point.
  REAL, ALLOCATABLE :: y     (:)      ! The y-coord. of the physical and
                                      ! computational grid. Defined at v-point.
  REAL, ALLOCATABLE :: z     (:)      ! The z-coord. of the computational grid.
                                      ! Defined at w-point on the staggered grid.
  REAL, ALLOCATABLE :: zp    (:,:,:)  ! The physical height coordinate defined at
                                      ! w-point of the staggered grid.

  REAL, ALLOCATABLE :: u     (:,:,:)  ! Total u-velocity (m/s)
  REAL, ALLOCATABLE :: v     (:,:,:)  ! Total v-velocity (m/s)
  REAL, ALLOCATABLE :: w     (:,:,:)  ! Total w-velocity (m/s)
  REAL, ALLOCATABLE :: ptprt (:,:,:)  ! Perturbation potential temperature
                                      ! from that of base state atmosphere (K)
  REAL, ALLOCATABLE :: pprt  (:,:,:)  ! Perturbation pressure from that
                                      ! of base state atmosphere (Pascal)
  REAL, ALLOCATABLE :: qv    (:,:,:)  ! Water vapor specific humidity (kg/kg)
  REAL, ALLOCATABLE :: qc    (:,:,:)  ! Cloud water mixing ratio (kg/kg)
  REAL, ALLOCATABLE :: qr    (:,:,:)  ! Rain water mixing ratio (kg/kg)
  REAL, ALLOCATABLE :: qi    (:,:,:)  ! Cloud ice mixing ratio (kg/kg)
  REAL, ALLOCATABLE :: qs    (:,:,:)  ! Snow mixing ratio (kg/kg)
  REAL, ALLOCATABLE :: qh    (:,:,:)  ! Hail mixing ratio (kg/kg)
  REAL, ALLOCATABLE :: tke   (:,:,:)  ! Turbulent Kinetic Energy ((m/s)**2)
  REAL, ALLOCATABLE :: kmh   (:,:,:)  ! Horizontal turb. mixing coef. for
                                      ! momentum. ( m**2/s )
  REAL, ALLOCATABLE :: kmv   (:,:,:)  ! Vertical turb. mixing coef. for
                                      ! momentum. ( m**2/s )
  REAL, ALLOCATABLE :: ubar  (:,:,:)  ! Base state u-velocity (m/s)
  REAL, ALLOCATABLE :: vbar  (:,:,:)  ! Base state v-velocity (m/s)
  REAL, ALLOCATABLE :: wbar  (:,:,:)  ! Base state w-velocity (m/s)
  REAL, ALLOCATABLE :: ptbar (:,:,:)  ! Base state potential temperature (K)
  REAL, ALLOCATABLE :: pbar  (:,:,:)  ! Base state pressure (Pascal)
  REAL, ALLOCATABLE :: rhobar(:,:,:)  ! Base state density rhobar
  REAL, ALLOCATABLE :: qvbar (:,:,:)  ! Base state water vapor specific humidity
                                      ! (kg/kg)
  INTEGER :: nstyps                   ! Number of soil type
!  PARAMETER ( nstyps = 4 )

  INTEGER, ALLOCATABLE :: soiltyp(:,:,:) ! Soil type
  REAL, ALLOCATABLE :: stypfrct(:,:,:)   ! Soil type fraction
  INTEGER, ALLOCATABLE :: vegtyp (:,:)   ! Vegetation type
  REAL, ALLOCATABLE :: lai    (:,:)      ! Leaf Area Index
  REAL, ALLOCATABLE :: roufns (:,:)      ! Surface roughness
  REAL, ALLOCATABLE :: veg    (:,:)      ! Vegetation fraction

  REAL, ALLOCATABLE :: tsfc   (:,:,:)    ! Soil Skin Temperature (K)
  REAL, ALLOCATABLE :: tsoil  (:,:,:)    ! Deep soil temperature (K)
  REAL, ALLOCATABLE :: wetsfc (:,:,:)    ! Surface soil moisture
  REAL, ALLOCATABLE :: wetdp  (:,:,:)    ! Deep soil moisture
  REAL, ALLOCATABLE :: wetcanp(:,:,:)    ! Canopy water amount
  REAL, ALLOCATABLE :: snowdpth(:,:)            ! Snow depth (m)

  REAL, ALLOCATABLE :: raing(:,:)        ! Grid supersaturation rain
  REAL, ALLOCATABLE :: rainc(:,:)        ! Cumulus convective rain
  REAL, ALLOCATABLE :: raint(:,:)        ! Total rain (rainc+raing)

  REAL, ALLOCATABLE :: prcrate(:,:,:)    ! precipitation rate (kg/(m**2*s))
                                         ! prcrate(1,1,1) = total precip. rate
                                         ! prcrate(1,1,2) = grid scale precip. rate
                                         ! prcrate(1,1,3) = cumulus precip. rate
                                         ! prcrate(1,1,4) = microphysics precip. rate

  REAL, ALLOCATABLE :: radfrc(:,:,:)     ! Radiation forcing (K/s)
  REAL, ALLOCATABLE :: radsw (:,:)       ! Solar radiation reaching the surface
  REAL, ALLOCATABLE :: rnflx (:,:)       ! Net radiation flux absorbed by surface

  REAL, ALLOCATABLE :: usflx (:,:)       ! Surface flux of u-momentum (kg/(m*s**2))
  REAL, ALLOCATABLE :: vsflx (:,:)       ! Surface flux of v-momentum (kg/(m*s**2))
  REAL, ALLOCATABLE :: ptsflx(:,:)       ! Surface heat flux (K*kg/(m*s**2))
  REAL, ALLOCATABLE :: qvsflx(:,:)       ! Surface moisture flux (kg/(m**2*s))
!
!-----------------------------------------------------------------------
!
!  Arrays derived from the read-in arrays
!
!-----------------------------------------------------------------------
!
  REAL, ALLOCATABLE :: uprt  (:,:,:)  ! Perturbation u-velocity (m/s)
  REAL, ALLOCATABLE :: vprt  (:,:,:)  ! Perturbation v-velocity (m/s)
  REAL, ALLOCATABLE :: wprt  (:,:,:)  ! Perturbation w-velocity (m/s)
  REAL, ALLOCATABLE :: pt    (:,:,:)  ! Total poten
  REAL, ALLOCATABLE :: qvprt (:,:,:)

  REAL, ALLOCATABLE :: xc   (:,:,:)   ! x-coor of sacalar point (km)
  REAL, ALLOCATABLE :: yc   (:,:,:)   ! y-coor of sacalar point (km)
  REAL, ALLOCATABLE :: zc   (:,:,:)   ! z-coor of sacalar point in computational
                                      ! space (km)
  REAL, ALLOCATABLE :: zpc  (:,:,:)   ! z-coor of sacalar point in physical
                                      ! space (km)

  REAL, ALLOCATABLE :: hterain(:,:)   ! The height of the terrain.

!-----------------------------------------------------------------------
!
!  Arrays for plots on constant pressure levels
!
!-----------------------------------------------------------------------
!
  REAL, ALLOCATABLE :: hgt  (:,:,:)   ! z-coor of sacalar point in physical
                                      ! space (10m)
  REAL, ALLOCATABLE :: t700  (:,:)    ! Temperature (K) on 700mb pressure grids
  REAL, ALLOCATABLE :: qs700  (:,:)   ! Temperature (K) on 700mb pressure grids

  REAL, ALLOCATABLE :: algpzc(:,:,:)  ! -log(pressure) at scalar grid points

!-----------------------------------------------------------------------
!
!  Array for CAPE , CIN
!
!-----------------------------------------------------------------------
!
  REAL, ALLOCATABLE :: wrk1(:),wrk2(:),wrk3(:),wrk4(:),wrk5(:),wrk6(:)
  REAL, ALLOCATABLE :: wrk7(:),wrk8(:),wrk9(:),wrk10(:),wrk11(:),wrk12(:)

  REAL, ALLOCATABLE :: lcl(:,:)    ! The lifting condensation level
  REAL, ALLOCATABLE :: lfc(:,:)    ! Level of free convection
  REAL, ALLOCATABLE :: el(:,:)     ! Equilibrium level
  REAL, ALLOCATABLE :: twdf(:,:)   ! Max wet-bulb potential temperature difference
  REAL, ALLOCATABLE :: mbe(:,:)    ! Moist Convective Potential Energy
                        !   (MCAPE, includes water loading)

  REAL, ALLOCATABLE :: li1(:,:)    ! Lifted Index (K)
  REAL, ALLOCATABLE :: cape1(:,:)  ! CAPE  (J/kg)
  REAL, ALLOCATABLE :: cin1(:,:)   ! CIN   (J/kg)
  REAL, ALLOCATABLE :: thet(:,:)   ! theta_E (K)
  REAL, ALLOCATABLE :: heli(:,:)   ! helicity
  REAL, ALLOCATABLE :: brn1(:,:)   ! Bulk Richardson Number (Weisman and Klemp)
  REAL, ALLOCATABLE :: brnu(:,:)   ! Shear parameter of BRN, "U"
  REAL, ALLOCATABLE :: srlfl(:,:)  ! storm-relative low-level flow (0-2km AGL)
  REAL, ALLOCATABLE :: srmfl(:,:)  ! storm-relative mid-level flow (2-9km AGL)
  REAL, ALLOCATABLE :: shr37(:,:)  ! 7km - 3km wind shear
  REAL, ALLOCATABLE :: ustrm(:,:)  ! Estimated storm motion (Bob Johns)
  REAL, ALLOCATABLE :: vstrm(:,:)  ! Estimated storm motion (Bob Johns)

  REAL, ALLOCATABLE :: capst(:,:)  ! cap strength
  REAL, ALLOCATABLE :: blcon(:,:)  ! boundary llayer convergence

  REAL, ALLOCATABLE :: ct(:,:)     ! convective temperature

  REAL :: oe                       ! function for caculate theta_e
  REAL :: alttostpr

  REAL, ALLOCATABLE :: sinlat(:,:) ! Sin of latitude at each grid point
  SAVE sinlat
  REAL :: lat1,lon1,lat2,lon2
  REAL :: pw1d,totw
  REAL :: pi
!
!-----------------------------------------------------------------------
! OUTPUT arrays
!-----------------------------------------------------------------------
  INTEGER :: fileun
  REAL :: resltn
  INTEGER :: houroday, minohour, dayoweek, dateomon, month1, year1
  REAL, ALLOCATABLE :: mslp(:,:)
  REAL, ALLOCATABLE :: hgt250(:,:), vort250(:,:), uwind250(:,:), vwind250(:,:)
  REAL, ALLOCATABLE :: hgt500(:,:), vort500(:,:), uwind500(:,:), vwind500(:,:)
  REAL, ALLOCATABLE :: hgt850(:,:), vort850(:,:), uwind850(:,:), vwind850(:,:)
  REAL, ALLOCATABLE :: temp850(:,:)
  REAL, ALLOCATABLE :: thck(:,:), rh700(:,:), omega700(:,:),temp2m(:,:), dewp2m(:,:)
  REAL, ALLOCATABLE :: accppt(:,:), convppt(:,:)
  REAL, ALLOCATABLE :: tpptold(:,:), cpptold(:,:)
  REAL, ALLOCATABLE :: sreh(:,:),li(:,:),cape(:,:),brn(:,:),cin(:,:),llmc(:,:),pw(:,:)
  REAL, ALLOCATABLE :: cmpref(:,:)
  INTEGER :: ni,nj
  REAL, ALLOCATABLE :: f2(:,:),x2(:,:),y2(:,:)
  REAL, ALLOCATABLE :: x1(:,:),y1(:,:)
!
!-----------------------------------------------------------------------
!  Output Grid
!-----------------------------------------------------------------------
!
  INTEGER :: max2d_out2
!  PARAMETER (MAX2D_OUT2=max((nx-3)*(ny-3),MAX2D_OUT))
  PARAMETER (max2d_out2=max2d_out)

  REAL :: mslpn(max2d_out2)
  REAL :: hgt250n(max2d_out2), vort250n(max2d_out2),                    &
       uwind250n(max2d_out2), vwind250n(max2d_out2)
  REAL :: hgt500n(max2d_out2), vort500n(max2d_out2),                    &
       uwind500n(max2d_out2), vwind500n(max2d_out2)
  REAL :: hgt850n(max2d_out2), vort850n(max2d_out2),                    &
       uwind850n(max2d_out2), vwind850n(max2d_out2)
  REAL :: temp850n(max2d_out2)
  REAL :: thckn(max2d_out2), rh700n(max2d_out2), omega700n(max2d_out2), &
       temp2mn(max2d_out2), dewp2mn(max2d_out2)
  REAL :: accpptn(max2d_out2), convpptn(max2d_out2)
  REAL :: srehn(max2d_out2), lin(max2d_out2), capen(max2d_out2),        &
       brnn(max2d_out2),                                                &
       cinn(max2d_out2), llmcn(max2d_out2), pwn(max2d_out2)
  REAL :: cmprefn(max2d_out2)
  REAL :: xn(max2d_out2),yn(max2d_out2)
  REAL :: x1n(max2d_out2),y1n(max2d_out2)
  REAL :: latn(max2d_out2),lonn(max2d_out2)
!
!-----------------------------------------------------------------------
!
!  Temporary work arrays for general use
!
!-----------------------------------------------------------------------
!
  REAL, ALLOCATABLE :: tem1(:,:,:)
  REAL, ALLOCATABLE :: tem2(:,:,:)
  REAL, ALLOCATABLE :: tem3(:,:,:)
  REAL, ALLOCATABLE :: tem4(:,:,:)
  REAL, ALLOCATABLE :: tem5(:,:,:)
  REAL, ALLOCATABLE :: tem6(:,:,:)

  INTEGER, ALLOCATABLE :: item(:,:)
  REAL, ALLOCATABLE :: var2(:,:)
  REAL, ALLOCATABLE :: pt2m(:,:),qv2m(:,:)
  REAL, ALLOCATABLE :: vtwo1(:,:),vtwo2(:,:),vtwo3(:,:)
  REAL, ALLOCATABLE :: vtwo4(:,:),vtwo5(:,:),vtwo6(:,:)
  REAL, ALLOCATABLE :: qvsfc   (:,:,:)    ! Soil Skin effective qv
  REAL, ALLOCATABLE :: psf(:,:),psl(:,:)

  INTEGER :: ijq

  INTEGER, ALLOCATABLE :: item1(:,:),item2(:,:),item3(:,:),item4(:,:),       &
          item5(:,:),item6(:,:),item7(:,:)
  REAL, ALLOCATABLE :: rtem1(:,:),rtem2(:,:),rtem3(:,:),rtem4(:,:),          &
       rtem5(:,:),rtem6(:,:),rtem7(:,:)
!
!-----------------------------------------------------------------------
!
!  Work arrays to be used in interpolation subroutines
!
!-----------------------------------------------------------------------
!
  REAL, ALLOCATABLE :: b1(:,:)
  REAL, ALLOCATABLE :: t1(:,:),t2(:,:),t3(:,:),t4(:,:)
!
!-----------------------------------------------------------------------
!
!  Common blocks for plotting control parameters
!
!-----------------------------------------------------------------------
  REAL :: z01                      ! altitude (meters) of a horizontal
                                   ! slice so specified
  COMMON /sliceh/z01
!-----------------------------------------------------------------------
!
!  Misc. local variables
!
!-----------------------------------------------------------------------
  REAL :: swx,ctrx,swy,ctry,tem
  REAL :: ctrx1, ctry1
  REAL :: iorig
  REAL :: rho
!
  INTEGER,ALLOCATABLE :: slice(:,:)
  REAL,ALLOCATABLE :: xi1(:), yi1(:),xi2(:),yi2(:),zi1(:),ppil(:)
  CHARACTER (LEN=80) :: label
  CHARACTER (LEN=80) :: filename
  CHARACTER (LEN=80) :: grdbasfn
  CHARACTER (LEN=80) :: hisfile(25)
  CHARACTER (LEN=80) :: outfile(25)
  INTEGER :: hinfmt,nchin
  INTEGER :: nslice(7), indxslic
  
  INTEGER :: ireturn
  INTEGER :: mode,lengbf,lenfil
  INTEGER :: i,j,k,ibgn,iend,jbgn,jend,kbgn,kend,ist,jst,kst,iob
  INTEGER :: ibgn1, iend1
  INTEGER :: nxpic,nypic,islice,jslice,kslice,layout
  REAL :: ctime,angl,aspect
!  real xr,yr,zr,x1,x2,y1,y2,z1,z2,xlimit,ylimit
  REAL :: zmax,zmin
  REAL :: factor,pmb,drot
  INTEGER :: imax,jmax,kmax,imin,jmin,kmin
  REAL :: wmin, wmax, xorold, yorold, tk, tdk, psta
  REAL :: xbgn,xend,ybgn,yend,zbgn,zend
  LOGICAL :: fexist
  INTEGER :: onvf
  INTEGER :: nf
  INTEGER :: lenstag, linstit, lmodel
!
!-----------------------------------------------------------------------
  INTEGER :: intrpl   ! flag indicating whether to interpolate output

  REAL :: truelat(2)
  REAL :: trltnew(2),trlnnew,dxnew
  INTEGER :: mptjnew
  INTEGER :: nhisfile
  INTEGER :: iout(25),icape,iplot,iterr,iaccu
  NAMELIST /history_data/ hinfmt, nhisfile, grdbasfn,hisfile
  NAMELIST /output_data/ outfile,iout,icape,iplot,iterr,iaccu
  NAMELIST /output_grid/ mptjnew,trltnew,trlnnew,dxnew,                 &
      qni,qnj,qswcorn,qswcorw,qnecorn,qnecorw,intrpl,                   &
      enstag,instit,model
  INTEGER :: istatus
!
!-----------------------------------------------------------------------
!
!  Function f_qvsat and inline directives for Cray PVP
!
!-----------------------------------------------------------------------
!
  REAL :: f_qvsat

!fpp$ expand (f_qvsat)
!dir$ inline always f_qvsat

!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

  WRITE(6,'(/ 16(/5x,a)//)')                                            &
      '###############################################################', &
      '###############################################################', &
      '####                                                       ####', &
      '####                 Welcome to ARPSENSCV                  ####', &
      '####                                                       ####', &
      '####       A graphic analysis program for model ARPS 4.5   ####', &
      '####                                                       ####', &
      '####               Program version 4.5                     ####', &
      '####                                                       ####', &
      '####            The graphic plotting is based              ####', &
      '####              on graphic package ZXPLOT                ####', &
      '####               by Ming Xue CAPS/OU                     ####', &
      '####                                                       ####', &
      '###############################################################', &
      '###############################################################'

!-----------------------------------------------------------------------
!
!  Read grid/base name, then get the dimensions
!
!-----------------------------------------------------------------------

  READ(5,history_data,END=100)
  WRITE(6,'(/a/)')'Namelist block history_data sucessfully read in.'

  IF ( hinfmt == 9 ) GO TO 10

  lengbf = len_trim(grdbasfn)
  WRITE(6,'(a,a/)')' The grid/base name is ', grdbasfn(1:lengbf)

  DO nf=1,nhisfile
    lenfil = len_trim(hisfile(nf))
    WRITE(6,'(1x,a,a)') 'Input values for hisfile: ',trim(hisfile(nf))
  END DO

  GO TO 10
  100   WRITE(6,'(a)')                                                  &
         'Error reading NAMELIST block history_data. STOP'
         STOP
  10 CONTINUE
!
  CALL get_dims_from_data(hinfmt,trim(grdbasfn),                    &
       nx,ny,nz,nstyps, ireturn)

  IF( ireturn /= 0 ) THEN
    PRINT*,'Problem occured when trying to get dimensions from data.'
    PRINT*,'Program stopped.'
    STOP
  END IF

  ni=nx-1   ! F. KONG add to fix a bug for zero ni,nj
  nj=ny-1

  WRITE(6,'(3(a,i5))') 'nx =',nx,', ny=',ny,', nz=',nz

  print*,'nstyps =', nstyps
  
!-----------------------------------------------------------------------
!
! Allocate nx, ny, nz dependent arrays and initiliaze to zero
!
!-----------------------------------------------------------------------
  ALLOCATE(  x=0
  ALLOCATE(  y=0
  ALLOCATE(  z=0
  ALLOCATE(zp(nx,ny,nz),STAT=istatus)
  zp=0
  ALLOCATE(  u=0
  ALLOCATE(  v=0
  ALLOCATE(  w=0
  ALLOCATE(ptprt(nx,ny,nz),STAT=istatus)
  ptprt=0
  ALLOCATE(pprt(nx,ny,nz),STAT=istatus)
  pprt=0
  ALLOCATE(qv(nx,ny,nz),STAT=istatus)
  qv=0
  ALLOCATE(qc(nx,ny,nz),STAT=istatus)
  qc=0
  ALLOCATE(qr(nx,ny,nz),STAT=istatus)
  qr=0
  ALLOCATE(qi(nx,ny,nz),STAT=istatus)
  qi=0
  ALLOCATE(qs(nx,ny,nz),STAT=istatus)
  qs=0
  ALLOCATE(qh(nx,ny,nz),STAT=istatus)
  qh=0
  ALLOCATE(tke(nx,ny,nz),STAT=istatus)
  tke=0
  ALLOCATE(kmh(nx,ny,nz),STAT=istatus)
  kmh=0
  ALLOCATE(kmv(nx,ny,nz),STAT=istatus)
  kmv=0
  ALLOCATE(ubar(nx,ny,nz),STAT=istatus)
  ubar=0
  ALLOCATE(vbar(nx,ny,nz),STAT=istatus)
  vbar=0
  ALLOCATE(wbar(nx,ny,nz),STAT=istatus)
  wbar=0
  ALLOCATE(ptbar(nx,ny,nz),STAT=istatus)
  ptbar=0
  ALLOCATE(pbar(nx,ny,nz),STAT=istatus)
  pbar=0
  ALLOCATE(rhobar(nx,ny,nz),STAT=istatus)
  rhobar=0
  ALLOCATE(qvbar(nx,ny,nz),STAT=istatus)
  qvbar=0
  ALLOCATE(soiltyp(nx,ny,nstyps),STAT=istatus)
  soiltyp=0
  ALLOCATE(stypfrct(nx,ny,nstyps),STAT=istatus)
  stypfrct=0
  ALLOCATE(vegtyp(nx,ny),STAT=istatus)
  vegtyp=0
  ALLOCATE(lai(nx,ny),STAT=istatus)
  lai=0
  ALLOCATE(roufns(nx,ny),STAT=istatus)
  roufns=0
  ALLOCATE(veg(nx,ny),STAT=istatus)
  veg=0
  ALLOCATE(tsfc(nx,ny,0:nstyps),STAT=istatus)
  tsfc=0
  ALLOCATE(tsoil(nx,ny,0:nstyps),STAT=istatus)
  tsoil=0
  ALLOCATE(wetsfc(nx,ny,0:nstyps),STAT=istatus)
  wetsfc=0
  ALLOCATE(wetdp(nx,ny,0:nstyps),STAT=istatus)
  wetdp=0
  ALLOCATE(wetcanp(nx,ny,0:nstyps),STAT=istatus)
  wetcanp=0
  ALLOCATE(snowdpth(nx,ny),STAT=istatus)
  snowdpth=0
  ALLOCATE(raing(nx,ny),STAT=istatus)
  raing=0
  ALLOCATE(rainc(nx,ny),STAT=istatus)
  rainc=0
  ALLOCATE(raint(nx,ny),STAT=istatus)
  raint=0
  ALLOCATE(prcrate(nx,ny,4),STAT=istatus)
  prcrate=0
  ALLOCATE(radfrc(nx,ny,nz),STAT=istatus)
  radfrc=0
  ALLOCATE(radsw(nx,ny),STAT=istatus)
  radsw=0
  ALLOCATE(rnflx(nx,ny),STAT=istatus)
  rnflx=0
  ALLOCATE(usflx(nx,ny),STAT=istatus)
  usflx=0
  ALLOCATE(vsflx(nx,ny),STAT=istatus)
  vsflx=0
  ALLOCATE(ptsflx(nx,ny),STAT=istatus)
  ptsflx=0
  ALLOCATE(qvsflx(nx,ny),STAT=istatus)
  qvsflx=0
  ALLOCATE(uprt(nx,ny,nz),STAT=istatus)
  uprt=0
  ALLOCATE(vprt(nx,ny,nz),STAT=istatus)
  vprt=0
  ALLOCATE(wprt(nx,ny,nz),STAT=istatus)
  wprt=0
  ALLOCATE(pt(nx,ny,nz),STAT=istatus)
  pt=0
  ALLOCATE(qvprt(nx,ny,nz),STAT=istatus)
  qvprt=0
  ALLOCATE(xc(nx,ny,nz),STAT=istatus)
  xc=0
  ALLOCATE(yc(nx,ny,nz),STAT=istatus)
  yc=0
  ALLOCATE(zc(nx,ny,nz),STAT=istatus)
  zc=0
  ALLOCATE(zpc(nx,ny,nz),STAT=istatus)
  zpc=0
  ALLOCATE(hterain(nx,ny),STAT=istatus)
  hterain=0
  ALLOCATE(hgt(nx,ny,nz),STAT=istatus)
  hgt=0
  ALLOCATE(t700(nx,ny),STAT=istatus)
  t700=0
  ALLOCATE(qs700(nx,ny),STAT=istatus)
  qs700=0
  ALLOCATE(algpzc(nx,ny,nz),STAT=istatus)
  algpzc=0
  ALLOCATE(wrk1(nz),STAT=istatus)
  wrk1=0
  ALLOCATE(wrk2(nz),STAT=istatus)
  wrk2=0
  ALLOCATE(wrk3(nz),STAT=istatus)
  wrk3=0
  ALLOCATE(wrk4(nz),STAT=istatus)
  wrk4=0
  ALLOCATE(wrk5(nz),STAT=istatus)
  wrk5=0
  ALLOCATE(wrk6(nz),STAT=istatus)
  wrk6=0
  ALLOCATE(wrk7(nz),STAT=istatus)
  wrk7=0
  ALLOCATE(wrk8(nz),STAT=istatus)
  wrk8=0
  ALLOCATE(wrk9(nz),STAT=istatus)
  wrk9=0
  ALLOCATE(wrk10(nz),STAT=istatus)
  wrk10=0
  ALLOCATE(wrk11(nz),STAT=istatus)
  wrk11=0
  ALLOCATE(wrk12(nz),STAT=istatus)
  wrk12=0
  ALLOCATE(lcl(nx,ny),STAT=istatus)
  lcl=0
  ALLOCATE(lfc(nx,ny),STAT=istatus)
  lfc=0
  ALLOCATE(el(nx,ny),STAT=istatus)
  el=0
  ALLOCATE(twdf(nx,ny),STAT=istatus)
  twdf=0
  ALLOCATE(mbe(nx,ny),STAT=istatus)
  mbe=0
  ALLOCATE(li1(nx,ny),STAT=istatus)
  li1=0
  ALLOCATE(cape1(nx,ny),STAT=istatus)
  cape1=0
  ALLOCATE(cin1(nx,ny),STAT=istatus)
  cin1=0
  ALLOCATE(thet(nx,ny),STAT=istatus)
  thet=0
  ALLOCATE(heli(nx,ny),STAT=istatus)
  heli=0
  ALLOCATE(brn1(nx,ny),STAT=istatus)
  brn1=0
  ALLOCATE(brnu(nx,ny),STAT=istatus)
  brnu=0
  ALLOCATE(srlfl(nx,ny),STAT=istatus)
  srlfl=0
  ALLOCATE(srmfl(nx,ny),STAT=istatus)
  srmfl=0
  ALLOCATE(shr37(nx,ny),STAT=istatus)
  shr37=0
  ALLOCATE(ustrm(nx,ny),STAT=istatus)
  ustrm=0
  ALLOCATE(vstrm(nx,ny),STAT=istatus)
  vstrm=0
  ALLOCATE(capst(nx,ny),STAT=istatus)
  capst=0
  ALLOCATE(blcon(nx,ny),STAT=istatus)
  blcon=0
  ALLOCATE(ct(nx,ny),STAT=istatus)
  ct=0
  ALLOCATE(sinlat(nx,ny),STAT=istatus)
  sinlat=0
  ALLOCATE(mslp(nx-1,ny-1),STAT=istatus)
  mslp=0
  ALLOCATE(hgt250(nx-1,ny-1),STAT=istatus)
  hgt250=0
  ALLOCATE(vort250(nx-1,ny-1),STAT=istatus)
  vort250=0
  ALLOCATE(uwind250(nx-1,ny-1),STAT=istatus)
  uwind250=0
  ALLOCATE(vwind250(nx-1,ny-1),STAT=istatus)
  vwind250=0
  ALLOCATE(hgt500(nx-1,ny-1),STAT=istatus)
  hgt500=0
  ALLOCATE(vort500(nx-1,ny-1),STAT=istatus)
  vort500=0
  ALLOCATE(uwind500(nx-1,ny-1),STAT=istatus)
  uwind500=0
  ALLOCATE(vwind500(nx-1,ny-1),STAT=istatus)
  vwind500=0
  ALLOCATE(hgt850(nx-1,ny-1),STAT=istatus)
  hgt850=0
  ALLOCATE(vort850(nx-1,ny-1),STAT=istatus)
  vort850=0
  ALLOCATE(uwind850(nx-1,ny-1),STAT=istatus)
  uwind850=0
  ALLOCATE(vwind850(nx-1,ny-1),STAT=istatus)
  vwind850=0
  ALLOCATE(temp850(nx-1,ny-1),STAT=istatus)
  temp850=0
  ALLOCATE(thck(nx-1,ny-1),STAT=istatus)
  thck=0
  ALLOCATE(rh700(nx-1,ny-1),STAT=istatus)
  rh700=0
  ALLOCATE(omega700(nx-1,ny-1),STAT=istatus)
  omega700=0
  ALLOCATE(temp2m(nx-1,ny-1),STAT=istatus)
  temp2m=0
  ALLOCATE(dewp2m(nx-1,ny-1),STAT=istatus)
  dewp2m=0
  ALLOCATE(accppt(nx-1,ny-1),STAT=istatus)
  accppt=0
  ALLOCATE(convppt(nx-1,ny-1),STAT=istatus)
  convppt=0
  ALLOCATE(tpptold(nx-1,ny-1),STAT=istatus)
  tpptold=0
  ALLOCATE(cpptold(nx-1,ny-1),STAT=istatus)
  cpptold=0
  ALLOCATE(sreh(nx-1,ny-1),STAT=istatus)
  sreh=0
  ALLOCATE(li(nx-1,ny-1),STAT=istatus)
  li=0
  ALLOCATE(cape(nx-1,ny-1),STAT=istatus)
  cape=0
  ALLOCATE(brn(nx-1,ny-1),STAT=istatus)
  brn=0
  ALLOCATE(cin(nx-1,ny-1),STAT=istatus)
  cin=0
  ALLOCATE(llmc(nx-1,ny-1),STAT=istatus)
  llmc=0
  ALLOCATE(pw(nx-1,ny-1),STAT=istatus)
  pw=0
  ALLOCATE(cmpref(nx-1,ny-1),STAT=istatus)
  cmpref=0
  ALLOCATE(f2(nx-1,ny-1),STAT=istatus)
  f2=0
  ALLOCATE(x2(nx-1,ny-1),STAT=istatus)
  x2=0
  ALLOCATE(y2(nx-1,ny-1),STAT=istatus)
  y2=0
  ALLOCATE(x1(nx-1,ny-1),STAT=istatus)
  x1=0
  ALLOCATE(y1(nx-1,ny-1),STAT=istatus)
  y1=0
  ALLOCATE(tem1(nx,ny,nz),STAT=istatus)
  tem1=0
  ALLOCATE(tem2(nx,ny,nz),STAT=istatus)
  tem2=0
  ALLOCATE(tem3(nx,ny,nz),STAT=istatus)
  tem3=0
  ALLOCATE(tem4(nx,ny,nz),STAT=istatus)
  tem4=0
  ALLOCATE(tem5(nx,ny,nz),STAT=istatus)
  tem5=0
  ALLOCATE(tem6(nx,ny,nz),STAT=istatus)
  tem6=0
  ALLOCATE(item(nx,ny),STAT=istatus)
  ALLOCATE(var2(nx,ny),STAT=istatus)
  var2=0
  ALLOCATE(pt2m(nx,ny),STAT=istatus)
  pt2m=0
  ALLOCATE(qv2m(nx,ny),STAT=istatus)
  qv2m=0
  ALLOCATE(vtwo1(nx,ny),STAT=istatus)
  vtwo1=0
  ALLOCATE(vtwo2(nx,ny),STAT=istatus)
  vtwo2=0
  ALLOCATE(vtwo3(nx,ny),STAT=istatus)
  vtwo3=0
  ALLOCATE(vtwo4(nx,ny),STAT=istatus)
  vtwo4=0
  ALLOCATE(vtwo5(nx,ny),STAT=istatus)
  vtwo5=0
  ALLOCATE(vtwo6(nx,ny),STAT=istatus)
  vtwo6=0
  ALLOCATE(qvsfc(nx,ny,0:nstyps),STAT=istatus)
  qvsfc=0
  ALLOCATE(psf(nx,ny),STAT=istatus)
  psf=0
  ALLOCATE(psl(nx,ny),STAT=istatus)
  psl=0
  ALLOCATE(item1(nx,ny),STAT=istatus)
  item1=0
  ALLOCATE(item2(nx,ny),STAT=istatus)
  item2=0
  ALLOCATE(item3(nx,ny),STAT=istatus)
  item3=0
  ALLOCATE(item4(nx,ny),STAT=istatus)
  item4=0
  ALLOCATE(item5(nx,ny),STAT=istatus)
  item5=0
  ALLOCATE(item6(nx,ny),STAT=istatus)
  item6=0
  ALLOCATE(item7(nx,ny),STAT=istatus)
  item7=0
  ALLOCATE(rtem1(nx,ny),STAT=istatus)
  rtem1=0
  ALLOCATE(rtem2(nx,ny),STAT=istatus)
  rtem2=0
  ALLOCATE(rtem3(nx,ny),STAT=istatus)
  rtem3=0
  ALLOCATE(rtem4(nx,ny),STAT=istatus)
  rtem4=0
  ALLOCATE(rtem5(nx,ny),STAT=istatus)
  rtem5=0
  ALLOCATE(rtem6(nx,ny),STAT=istatus)
  rtem6=0
  ALLOCATE(rtem7(nx,ny),STAT=istatus)
  rtem7=0
  ALLOCATE(b1(nx,ny),STAT=istatus)
  b1=0
  ALLOCATE(t1(nx,ny),STAT=istatus)
  t1=0
  ALLOCATE(t2(nx,ny),STAT=istatus)
  t2=0
  ALLOCATE(t3(nx,ny),STAT=istatus)
  t3=0
  ALLOCATE(t4(nx,ny),STAT=istatus)
  t4=0
  ALLOCATE(slice(nx+ny+nz,3),STAT=istatus)
  ALLOCATE(xi1(nx+ny),STAT=istatus)
  ALLOCATE(yi1(nx+ny),STAT=istatus)
  ALLOCATE(xi2(nx+ny),STAT=istatus)
  ALLOCATE(yi2(nx+ny),STAT=istatus)
  ALLOCATE(zi1(nz),STAT=istatus)
  ALLOCATE(ppil(nz),STAT=istatus)

!-----------------------------------------------------------------------
!
!  default output grid values (SAMEX settings):
!
!-----------------------------------------------------------------------
!
  intrpl = 1
  enstag = 'ARPSENS'
  instit = 'CAPS_OU'
  model = 'ARPS5.0.0'
  qswcorn = 27.0844
  qswcorw = -110.068
  qnecorn = 45.1552
  qnecorw = -68.10911
  qni = 117
  qnj = 81
  mptjnew=2
  trltnew(1)=39.0
  trltnew(2)=39.0
  trlnnew=-100.0
  dxnew=30.0

!  Read remaining namelist

  READ(5,output_data,END=200)
  WRITE(6,'(/a/)')'Namelist block output_data sucessfully read in.'
  GO TO 20
  200   WRITE(6,'(a)')                                                  &
         'Error reading NAMELIST block output_data.  ',                 &
         'Default values used.'
  20 CONTINUE
  DO nf=1,nhisfile
    WRITE(6,'(1x,a,a)') 'Input values for outfile: ',trim(outfile(nf))
  END DO

  READ(5,output_grid,END=250)
  WRITE(6,'(/a/)')'Namelist block output_grid sucessfully read in.'
  GO TO 30
  250   WRITE(6,'(a)')                                                  &
         'Error reading NAMELIST block output_grid.  ',                 &
         'Default values used.'
  30 CONTINUE

  lenstag = len_trim(enstag)
  linstit = len_trim(instit)
  lmodel = len_trim(model)
  WRITE (6,*) 'Descriptor enstag: ', enstag(1:lenstag)
  WRITE (6,*) 'Descriptor instit: ', instit(1:linstit)
  WRITE (6,*) 'Descriptor model: ', model(1:lmodel)
  WRITE(6,'(1x,a,3I5)')                                                 &
       'Grid variables intrpl,qni,qnj:', intrpl,qni,qnj
  WRITE(6,'(1x,a,4F11.5)')                                              &
       'Grid variables qswcorn, qswcorw, qnecorn, qnecorw:',            &
       qswcorn, qswcorw, qnecorn, qnecorw
  WRITE(6,'(1x,a,I5,4F10.3)')                                           &
       'Grid variables mptjnew,trltnew,trlnnew,dxnew:',                 &
                       mptjnew,trltnew,trlnnew,dxnew

  IF (intrpl == 1 .AND. qni*qnj > max2d_out2) THEN
    WRITE (6,'(a)')                                                     &
        'Error: qni*nqj > MAX2D_OUT (in enscv.inc).  Stopping.'
    STOP
  END IF
  IF (intrpl == 0 .AND. (nx-3)*(ny-3) > max2d_out2) THEN
    WRITE (6,'(a)')                                                     &
        'Error: (nx-3)*(ny-3) > MAX2D_OUT (in enscv.inc).  Stopping.'
    STOP
  END IF
  IF (intrpl == 0 .AND. (qni /= (nx-3).OR.qnj /= (ny-3) ) ) THEN
    WRITE (6,'(a,a)')                                                   &
        ' WARNING: qni =/= nx-3 and/or qnj=/=ny-3. ',                   &
       'correction is made TO make qni=nx-3 AND qnj=ny-3 '
    qni=nx-3
    qnj=ny-3
  END IF

  IF (iplot >= 1.OR.iterr > 0) CALL xdevic
  DO nf=1, nhisfile

    filename=hisfile(nf)
    lenfil = len_trim(filename)
    WRITE(6,'(a,a/)')                                                   &
        ' The file name of data set is ', filename(1:lenfil)

    15    CONTINUE ! also continue to read another time recode
                   ! from GrADS file

!-----------------------------------------------------------------------
!
!  Read all input data arrays
!
!-----------------------------------------------------------------------
!
    CALL dtaread(nx,ny,nz, nstyps,                                      &
                hinfmt, nchin,grdbasfn(1:lengbf),lengbf,                &
                filename(1:lenfil),lenfil,ctime,                         &
                x,y,z,zp, uprt ,vprt ,wprt ,ptprt, pprt ,               &
                qvprt, qc, qr, qi, qs, qh, tke,kmh,kmv,                 &
                ubar, vbar, wbar, ptbar, pbar, rhobar, qvbar,           &
                soiltyp,stypfrct,vegtyp,lai,roufns,veg,                 &
                tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth,               &
         raing,rainc,prcrate,                                           &
         radfrc,radsw,rnflx,                                            &
         usflx,vsflx,ptsflx,qvsflx,                                     &
                ireturn, tem1,tem2, tem3)
    year1=year
    month1=month
    dateomon=day
    houroday=hour
    minohour=ctime
    PRINT *,year1,month1,dateomon,houroday,ctime
!
!-----------------------------------------------------------------------
!
!  ireturn = 0 for a successful read
!  For hinfmt=9, i.e. the GraDs format data, ireturn is used as a
!  flag indicating if there is any data at more time level to be read.
!
!-----------------------------------------------------------------------
!
    IF( ireturn /= 0 .AND. hinfmt /= 9 ) GO TO 7000 !unsuccessful read
!
!-----------------------------------------------------------------------
!
!  write(6,'(/5x,a,i4)')
!    :     'mapproj in the data was ',mapproj

!  write(6,'(/5x,a,f10.3,a)')
!    :     'trulat1 in the data was ',trulat1,' degree North'

!  write(6,'(/5x,a,f10.3,a)')
!    :     'trulat2 in the data was ',trulat2,' degree North'

!  write(6,'(/5x,a,f10.3,a)')
!    :     'trulon in the data was ',trulon,' degree East'

!  write(6,'(/5x,a,e15.8)')
!    :     'sclfct in the data was ',sclfct
!
!-----------------------------------------------------------------------
!   Convert the rainfal arries and store the old ones.
    DO j=1,ny-1
      DO i=1,nx-1
        raint(i,j)=raing(i,j)+rainc(i,j)
!        IF (nf == 1) THEN     ! don't know what's the purpose (F.KONG)
!          raint(i,j)=0.0E0
!          rainc(i,j)=0.0E0
!        END IF
      END DO
    END DO
!      print *,raint(35,35),rainc(35,35),'rain'
    CALL ary2dcv(nx,ny,raint,ni,nj,accppt)
    CALL ary2dcv(nx,ny,rainc,ni,nj,convppt)
!      print *,accppt(35,35),convppt(35,35),'rain'

    IF (iaccu == 1) THEN
! convert accppt from 0h-now accul. rain to (Tnf-1->Tnf) accul. rain
! (by subtracting tpptold) and store the original value in tpptold
      IF (nf == 1) THEN  ! to clear the 1st time level values
      CALL ary2dcv(nx,ny,raint,ni,nj,tpptold)
      CALL ary2dcv(nx,ny,rainc,ni,nj,cpptold)
      END IF
    CALL pptsto(ni,nj,accppt,tpptold)
    CALL pptsto(ni,nj,convppt,cpptold)
!      print *,accppt(35,35),convppt(35,35),'rain'
    END IF

    IF (iout(nf) /= 1) CYCLE
!
!
!-----------------------------------------------------------------------
!
!  Set coordinates at the grid volume center
!
!-----------------------------------------------------------------------
!
    DO k=1,nz
      DO j=1,ny
        DO i=1,nx-1
          xc(i,j,k) = (x(i)+x(i+1))*0.5 /1000.0
        END DO
      END DO
    END DO

    DO k=1,nz
      DO j=1,ny-1
        DO i=1,nx
          yc(i,j,k) = (y(j)+y(j+1))*0.5 /1000.0
        END DO
      END DO
    END DO

    DO k=1,nz-1
      DO j=1,ny
        DO i=1,nx
          zc(i,j,k) = (z(k)+z(k+1))*0.5 /1000.0
          zpc(i,j,k)= (zp(i,j,k)+zp(i,j,k+1))*0.5 /1000.0
          hgt(i,j,k)=1000.0*zpc(i,j,k)
        END DO
      END DO
    END DO
    DO j=1,ny-1
      DO i=1,nx-1
        x2(i,j) = (x(i)+x(i+1))*0.5
        y2(i,j) = (y(j)+y(j+1))*0.5
        x1(i,j)=x2(i,j)/1000.0
        y1(i,j)=y2(i,j)/1000.0
      END DO
    END DO
!-----------------------------------------------------------------------
!  Find the lat/log of the points in the output grid, LatN, LonN
!  if interpolation is necessary. Set the flag for INTERPOLATION
!  Note that in NEWMAP, setmaps and setorig are used.
    IF (intrpl == 1) THEN
      CALL newmap(qni,qnj,latn,lonn,dxnew,                              &
                  qswcorn,qswcorw,trltnew,trlnnew,mptjnew)
      intrpl=1
    END IF
!-----------------------------------------------------------------------
! Set the map parameters and put the origin the same as used in defining
! xc and yc, x and y etc. The origin is actually point (2,2)
    dx = xc(3,2,2)-xc(2,2,2)
    dy = yc(2,3,2)-yc(2,2,2)
    truelat(1) = trulat1
    truelat(2) = trulat2

    CALL setmapr( mapproj, 1.0 , truelat , trulon)
                               ! set up parameters for map projection

    CALL lltoxy( 1,1, ctrlat,ctrlon, ctrx, ctry )
!     print *,ctrlat,ctrlon,ctrx,ctry
!     print *,dx,dy,nx,ny,xc(1,2,2),xc(2,2,2),yc(2,1,2),yc(2,2,2)

    swx = ctrx- (nx-3)*dx*0.5*1000.0
    swy = ctry- (ny-3)*dy*0.5*1000.0

!    print *, swx,swy,'swx,swy of the original data's physical domain'
    CALL setorig( 1, swx, swy)
!-----------------------------------------------------------------------
!  Calculate the Coriolis para. f=2wsin(lat) at the data's scalar points
    pi=ATAN(1.0)*4.0
    DO j=1,ny-1
      DO i=1,nx-1
        CALL xytoll( 1,1, x2(i,j),y2(i,j),lat1,lon1)
        f2(i,j)=SIN(lat1*pi/180.0)*2*omega
      END DO
    END DO
    PRINT *,f2(1,1),f2(nx-1,ny-1), 'f values at data map corners'
!-----------------------------------------------------------------------
!   Find the SW and NE corners of the outputgrid if they are not set
!   (and intrpl=0)
!   The second most SW or NE scalar point are chosen as these points
!   and so qni=nx-3 and qnj=nx-3 should be satisfied
!  *** for direct assignment instead of interpolation,set the new - grid
!   corners and the flag for NO INTERPOLATION.also calculate latN,lonN
    IF (intrpl == 0) THEN
      CALL xytoll(1,1,x2(2,2),y2(2,2),qswcorn,qswcorw)
      CALL xytoll(1,1,x2(nx-2,ny-2),y2(nx-2,ny-2),qnecorn,qnecorw)
      intrpl=0
      dxnew=dx
      trltnew(1)=truelat(1)
      trltnew(2)=truelat(2)
      trlnnew=trulon
      DO j=1,qnj
        DO i=1,qni
          CALL xytoll(1,1,x2(i+1,j+1),y2(i+1,j+1),lat1,lon1)
          ijq =  i  + qni*(j-1)
          latn(ijq)=lat1
          lonn(ijq)=lon1
        END DO
      END DO
    END IF
!
!-----------------------------------------------------------------------
!   Find the lat/lon of the SW and NE corners of the data grid
!   and print out for comparison with the output grid's corners
    CALL xytoll( 1,1, x2(1,1),y2(1,1),lat1,lon1)
    CALL xytoll( 1,1, x2(nx-1,ny-1),y2(nx-1,ny-1),lat2,lon2)
    PRINT *,lat1,lon1,lat2,lon2,' Data Map corners'
    PRINT *,qswcorn,qswcorw,qnecorn,qnecorw,'Output Map corners'
!
!-----------------------------------------------------------------------
!  Calculate the following
!  1. the x,y coordinates of the points in the output grid relatvie to
!  its SW corner,unit being km, used for plotting and CLLMC,x1N + y1N.
!  2. the x,y coordinates of the points in the output grid in the  old
!  grid,unit being m, used for interpolation,xN + yN. comparable to
!  x2,y2.
    DO j=1,qnj
      DO i=1,qni
        ijq =  i + qni*(j-1)
        x1n(ijq)=(i-1)*dxnew
        y1n(ijq)=(j-1)*dxnew
        CALL lltoxy( 1,1, latn(ijq),lonn(ijq), xn(ijq), yn(ijq) )
        IF (i == 1.AND.j == 1.OR.i == qni.AND.j == qnj) THEN
          PRINT *,'SW/NE of output grid, lat,lon,x and y in input grid'
          PRINT *,i,j,latn(ijq),lonn(ijq),xn(ijq),yn(ijq)
        END IF
      END DO
    END DO
!----------------------------------------------------------------------
!
    DO j=1,nj
      DO i=1,ni
        sreh(i,j)=0.0
        li(i,j)=0.0
        cape(i,j)=0.0
        brn(i,j)=0.0
        cin(i,j)=0.0
        llmc(i,j)=0.0
        pw(i,j)=0.0
      END DO
    END DO

!   combine prt and bar
    DO k=1,nz
      DO j=1,ny
        DO i=1,nx
          u(i,j,k)=uprt(i,j,k)+ubar(i,j,k)
          v(i,j,k)=vprt(i,j,k)+vbar(i,j,k)
          w(i,j,k)=wprt(i,j,k)+wbar(i,j,k)
          qv(i,j,k)=MAX(0.0,qvprt(i,j,k)+qvbar(i,j,k))
          pt(i,j,k)=ptprt(i,j,k)+ptbar(i,j,k)
          tem1(i,j,k)=pprt(i,j,k)+pbar(i,j,k)
          tem2(i,j,k)=0.0
          tem3(i,j,k)=0.0
          tem4(i,j,k)=0.0
        END DO
      END DO
    END DO
! convert u,v,w to scalar points.
    DO k=1,nz-1
      DO j=1,ny-1
        DO i=1,nx-1
          tem2(i,j,k)=0.5*(u(i,j,k)+u(i+1,j,k))
          tem3(i,j,k)=0.5*(v(i,j,k)+v(i,j+1,k))
          tem4(i,j,k)=0.5*(w(i,j,k)+w(i,j,k+1))
        END DO
      END DO
    END DO
    u=0
    v=0
    w=0
    DO k=1,nz-1
      DO j=1,ny-1
        DO i=1,nx-1
          u(i,j,k)=tem2(i,j,k)
          v(i,j,k)=tem3(i,j,k)
          w(i,j,k)=tem4(i,j,k)
        END DO
      END DO
    END DO
!
!----------------------------------------------------------------------
!
!    Calculate temperature (K) at ARPS grid points
!  Calculate dew-point temperature td using enhanced Teten's formula.
!
!-----------------------------------------------------------------------
!
    CALL temper(nx,ny,nz,pt, pprt ,pbar,tem2)
    CALL getdew(nx,ny,nz, 1,nx-1,1,ny-1,1,nz-1, tem1, tem2, qv,tem3)
    CALL getqvs(nx,ny,nz, 1,nx-1,1,ny-1,1,nz-1, tem1,tem2,tem4)
    ice=1
    CALL reflec(nx,ny,nz,rhobar,qr,qs,qh,tem5)
!
!-----------------------------------------------------------------------
!
!  Set negative logrithm pressure (Pascal) coordinates
!   for scalar and w grid points
!
!-----------------------------------------------------------------------
!
    DO k=1,nz-1
      DO j=1,ny-1
        DO i=1,nx-1
          algpzc(i,j,k) = -ALOG(tem1(i,j,k))
        END DO
      END DO
    END DO
!-----------------------------------------------------------------------
! calculate surface pressure (-log(ps)=psf),pmsl=psl,and moist. con. llmc
    CALL psfsl(nx,ny,tem2(1,1,2),qv(1,1,2),tem1(1,1,2),                 &
               hgt(1,1,2),zp(1,1,2),psf,psl)
    CALL cllmc(nx,ny,nz,u,v,qv,psf,zp,rhobar,algpzc,llmc,               &
                     x2,y2, vtwo1,vtwo2)
    CALL edgfill(llmc,1,nx-1,2,nx-2,1,ny-1,2,ny-2,                      &
               1,1,1,1)

!   calculate precipitable water and composite reflectivity
    DO j=1,ny-1
      DO i=1,nx-1
        cmpref(i,j)=0.0
        pw1d=0.0
        DO k=2,nz-1
!     totw=qv(i,j,k)+qc(i,j,k)+qr(i,j,k)
!  :       +qi(i,j,k)+qs(i,j,k)+qh(i,j,k)
          totw=qv(i,j,k)
          rho=tem1(i,j,k)/(rd*tem2(i,j,k)*(1.0+0.61*qv(i,j,k)))
          pw1d=pw1d+totw*(zp(i,j,k+1)-zp(i,j,k))*rho
          cmpref(i,j)=MAX(cmpref(i,j),tem5(i,j,k))
        END DO
        pw(i,j)=pw1d*1000.0
      END DO
    END DO

! calculate pt2m and qv2m by calling PTQV2M
    CALL ptqv2m(nx,ny,nz, nstyps,zp,                                    &
                u ,v ,w ,ptprt, pprt, qvprt,                            &
                ptbar, pbar, rhobar, qvbar,                             &
                soiltyp,stypfrct,vegtyp,lai,roufns,veg,                 &
                tsfc,wetsfc,wetcanp,                                    &
                pt2m,qv2m, vtwo1,vtwo2,                                 &
                qvsfc,vtwo3,                                            &
                vtwo4,vtwo5,vtwo6)
!-----------------------------------------------------------------------
!    Generate other 2-D data sets.
    z01=-ALOG(25000.0)
!      CALL secthrza(nx,ny,nz,hgt,algpzc,var2)
    CALL pintp(nx,ny,nz,hgt,algpzc,var2,3,z01,                          &
               tem2,qv,hgt(1,1,2),zp(1,1,2),psf)
    CALL ary2dcv(nx,ny,var2,ni,nj,hgt250)
!      CALL secthrza(nx,ny,nz,u,algpzc,var2)
    CALL pintp(nx,ny,nz,u,algpzc,var2,1,z01,                            &
               tem2,qv,hgt(1,1,2),zp(1,1,2),psf)
    CALL ary2dcv(nx,ny,var2,ni,nj,uwind250)
!      CALL secthrza(nx,ny,nz,v,algpzc,var2)
    CALL pintp(nx,ny,nz,v,algpzc,var2,1,z01,                            &
               tem2,qv,hgt(1,1,2),zp(1,1,2),psf)
    CALL ary2dcv(nx,ny,var2,ni,nj,vwind250)
    CALL vrtcalc(uwind250,vwind250,ni,nj,vort250,f2,x2,y2)
    z01=-ALOG(50000.0)
!      CALL secthrza(nx,ny,nz,hgt,algpzc,var2)
    CALL pintp(nx,ny,nz,hgt,algpzc,var2,3,z01,                          &
               tem2,qv,hgt(1,1,2),zp(1,1,2),psf)
    CALL ary2dcv(nx,ny,var2,ni,nj,hgt500)
!      CALL secthrz (nx,ny,nz,u,algpzc,var2)
    CALL pintp(nx,ny,nz,u,algpzc,var2,1,z01,                            &
               tem2,qv,hgt(1,1,2),zp(1,1,2),psf)
    CALL ary2dcv(nx,ny,var2,ni,nj,uwind500)
!      CALL secthrza(nx,ny,nz,v,algpzc,var2)
    CALL pintp(nx,ny,nz,v,algpzc,var2,1,z01,                            &
               tem2,qv,hgt(1,1,2),zp(1,1,2),psf)
    CALL ary2dcv(nx,ny,var2,ni,nj,vwind500)
    CALL vrtcalc(uwind500,vwind500,ni,nj,vort500,f2,x2,y2)
    z01=-ALOG(85000.0)
!      CALL secthrza(nx,ny,nz,hgt,algpzc,var2)
    CALL pintp(nx,ny,nz,hgt,algpzc,var2,3,z01,                          &
               tem2,qv,hgt(1,1,2),zp(1,1,2),psf)
    CALL ary2dcv(nx,ny,var2,ni,nj,hgt850)
!      CALL secthrza(nx,ny,nz,u,algpzc,var2)
    CALL pintp(nx,ny,nz,u,algpzc,var2,1,z01,                            &
               tem2,qv,hgt(1,1,2),zp(1,1,2),psf)
    CALL ary2dcv(nx,ny,var2,ni,nj,uwind850)
!      CALL secthrza(nx,ny,nz,v,algpzc,var2)
    CALL pintp(nx,ny,nz,v,algpzc,var2,1,z01,                            &
               tem2,qv,hgt(1,1,2),zp(1,1,2),psf)
    CALL ary2dcv(nx,ny,var2,ni,nj,vwind850)
    CALL vrtcalc(uwind850,vwind850,ni,nj,vort850,f2,x2,y2)
!      CALL secthrza(nx,ny,nz,tem2,algpzc,var2)
    CALL pintp(nx,ny,nz,tem2,algpzc,var2,2,z01,                         &
               tem2,qv,hgt(1,1,2),zp(1,1,2),psf)
    CALL ary2dcv(nx,ny,var2,ni,nj,temp850)
    z01=-ALOG(100000.0)
!      CALL secthrza(nx,ny,nz,hgt,algpzc,var2)
    CALL pintp(nx,ny,nz,hgt,algpzc,var2,3,z01,                          &
               tem2,qv,hgt(1,1,2),zp(1,1,2),psf)
    CALL ary2dcv(nx,ny,var2,ni,nj,thck)
    z01=-ALOG(70000.0)
!      CALL secthrza(nx,ny,nz,tem2,algpzc,t700)
    CALL pintp(nx,ny,nz,tem2,algpzc,t700,2,z01,                         &
               tem2,qv,hgt(1,1,2),zp(1,1,2),psf)
!      CALL secthrza(nx,ny,nz,tem4,algpzc,qs700)
    CALL pintp(nx,ny,nz,tem4,algpzc,qs700,1,z01,                        &
               tem2,qv,hgt(1,1,2),zp(1,1,2),psf)
!      CALL secthrza(nx,ny,nz,qv,algpzc,var2)
    CALL pintp(nx,ny,nz,qv,algpzc,var2,1,z01,                           &
               tem2,qv,hgt(1,1,2),zp(1,1,2),psf)
    CALL ary2dcv(nx,ny,var2,ni,nj,rh700)
!      CALL secthrza(nx,ny,nz,w,algpzc,var2)
    CALL pintp(nx,ny,nz,w,algpzc,var2,1,z01,                            &
               tem2,qv,hgt(1,1,2),zp(1,1,2),psf)
    CALL ary2dcv(nx,ny,var2,ni,nj,omega700)
    DO j=1,ny-1
      DO i=1,nx-1
        mslp(i,j) = psl(i,j)
        thck(i,j) = hgt500(i,j)-thck(i,j)
        rh700(i,j)= rh700(i,j)/qs700(i,j)*100.0
        temp850(i,j)=temp850(i,j)-273.15
        omega700(i,j)=omega700(i,j)*(-700000.0*g/287.053)/t700(i,j)
      END DO
    END DO

    IF (minohour < 1800) THEN
! temp2m, dewp2m scheme 1, direct interpolation from grid values
      z01=2.0
      CALL secthrza(nx,ny,nz,tem3,zc,var2)
      CALL ary2dcv(nx,ny,var2,ni,nj,dewp2m)
      z01=2.0
      CALL secthrza(nx,ny,nz,tem2,zc,var2)
      CALL ary2dcv(nx,ny,var2,ni,nj,temp2m)
    ELSE
!temp2m dewp2m scheme 2,boundary layer physics (PTQV2M's results are used)
      z01=2.0
      CALL secthrza(nx,ny,nz,tem1,zc,vtwo1)   !vtwo1, pressure
      DO j=1,ny-1
        DO i=1,nx-1
          vtwo2(i,j)=pt2m(i,j)*(vtwo1(i,j)/100000.0)**rddcp
        END DO
      END DO
      CALL ary2dcv(nx,ny,vtwo2,ni,nj,temp2m)
      CALL getdew(nx,ny,1,1,nx-1,1,ny-1,1,1, vtwo1, vtwo2, qv2m,vtwo3)
                                      !vtwo3,dew point temperature
      CALL ary2dcv(nx,ny,vtwo3,ni,nj,dewp2m)
!*** for the initial time,dewp at the level k=2 is used as dewp2m
!  CALL ARY2DCV(nx,ny,tem3(1,1,2),ni,nj,dewp2m)
    END IF

!  unit conversion
    DO j=1,ny-1
      DO i=1,nx-1
        temp2m(i,j)=temp2m(i,j)-273.15
        dewp2m(i,j)=dewp2m(i,j)-273.15
        IF (dewp2m(i,j) > temp2m(i,j)) dewp2m(i,j)=temp2m(i,j)
      END DO
    END DO
!----------------------------------------------------------------------
!
!    Calculate sea level pressure (mb)
!    Reduction method: Benjamin and Miller: 1990, MWR, vol.118, No.10,
!                   Page: 2100-2101
!
!-----------------------------------------------------------------------
!
!    do 700 j=1,ny-1
!    do 700 i=1,nx-1
!     p00 = 0.01*(pbar(i,j,2)+pprt(i,j,2))
!     if(p00.le.700.0) then
!       t00=tem2(i,j,2)
!     else
!       t00 = t700(i,j)*(p00/700.0)**ex1
!     end if
!     mslp(i,j) = p00*(1.0+gamma*zpc(i,j,2)/t00)**ex2
!     mslp(i,j) = psl(i,j)*0.01
!     thck(i,j) = hgt500(i,j)-8.0*(mslp(i,j)-1000.0)
!     thck(i,j) = hgt500(i,j)-thck(i,j)
    700     CONTINUE

!-----------------------------------------------------------------------
!
!  Calculate CAPE and CIN
!
!-----------------------------------------------------------------------
!
    DO k=1,nz-1
      DO j=1,ny-1
        DO i=1,nx-1
          tem4(i,j,k)=qv(i,j,k)/(1.-qv(i,j,k))
        END DO
      END DO
    END DO
    CALL edgfill(u,1,nx,1,nx-1,1,ny,1,ny-1,                             &
               1,nz,1,nz-1)
    CALL edgfill(v,1,nx,1,nx-1,1,ny,1,ny-1,                             &
               1,nz,1,nz-1)
    CALL edgfill(tem4,1,nx,1,nx-1,1,ny,1,ny-1,                          &
               1,nz,1,nz-1)
    CALL edgfill(tem1,1,nx,1,nx-1,1,ny,1,ny-1,                          &
               1,nz,1,nz-1)
    CALL edgfill(tem2,1,nx,1,nx-1,1,ny,1,ny-1,                          &
               1,nz,1,nz-1)
    CALL edgfill(tem3,1,nx,1,nx-1,1,ny,1,ny-1,                          &
               1,nz,1,nz-1)
    IF (icape == 2) THEN
      PRINT *, ' calling cape_ncep'
      CALL cape_ncep(nx,ny,nz,                                          &
                     tem2,qv,tem1,cape1,cin1,li1,                       &
                     tem3,tem4,tem5,tem6,                               &
                     vtwo1,vtwo2,vtwo3,vtwo4,vtwo5,item,                &
          item1,item2,item3,item4,item5,item6,item7,   & ! work arrays
      rtem1,rtem2,rtem3,rtem4,rtem5,rtem6,rtem7)   ! work arrays
      PRINT *, ' back from cape_ncep'
    ELSE IF (icape == 1) THEN
      PRINT *, ' calling arps_be'
      CALL arps_be(nx,ny,nz,                                            &
             tem1,zpc,tem2,tem4,                                        &
             lcl,lfc,el,twdf,li1,cape1,mbe,cin1,capst,                  &
             wrk1,wrk2,wrk3,wrk4,wrk5,wrk6,                             &
             wrk7,wrk8,wrk9,wrk10,wrk11,wrk12,t1)
      PRINT *, ' back from arps_be'
    ELSE
      PRINT *,'icape option is incorrect, chose 1 or 2'
      STOP
    END IF
    CALL ary2dcv(nx,ny,cape1,ni,nj,cape)
    CALL ary2dcv(nx,ny,cin1,ni,nj,cin)
    CALL ary2dcv(nx,ny,li1,ni,nj,li)
!
!-----------------------------------------------------------------------
!
!  Calculate helicity and other shear-related paramaters
!
!-----------------------------------------------------------------------
!
    PRINT *, ' calling calcshr'
    CALL xytomf(nx,ny,x,y,b1)
    CALL calcshr(nx,ny,nz,x,y,zp,b1,                                    &
        tem1,tem2,u,v,cape1,                                            &
        shr37,ustrm,vstrm,srlfl,srmfl,heli,brn1,brnu,blcon,             &
        t1,t2,t3,t4,tem4)
    PRINT *, ' back from calcshr'
    CALL ary2dcv(nx,ny,heli,ni,nj,sreh)
    CALL ary2dcv(nx,ny,brn1,ni,nj,brn)
!-----------------------------------------------------------------------
!   Interpolate the 2_d variables
!-----------------------------------------------------------------------

    IF (iterr == 1) THEN
      DO i=1,nx-1
        DO j=1,ny-1
          cmpref(i,j)=zp(i,j,2)
        END DO
      END DO
    END IF

    CALL      d2intrpa(x2,y2,nx-1,ny-1,                                 &
                       xn,yn,qni,qnj,intrpl,                            &
                       mslp,mslpn,                                      &
                       hgt250,hgt250n,                                  &
                       vort250,vort250n,                                &
                       uwind250,uwind250n,                              &
                       vwind250,vwind250n,                              &
                       hgt500,hgt500n,                                  &
                       vort500,vort500n,                                &
                       uwind500,uwind500n,                              &
                       vwind500,vwind500n,                              &
                       hgt850,hgt850n,                                  &
                       vort850,vort850n,                                &
                       uwind850,uwind850n,                              &
                       vwind850,vwind850n,                              &
                       temp850,temp850n,                                &
                       thck,thckn,                                      &
                       rh700,rh700n,                                    &
                       omega700,omega700n,                              &
                       temp2m,temp2mn,                                  &
                       dewp2m,dewp2mn,                                  &
                       accppt,accpptn,                                  &
                       convppt,convpptn,                                &
                       sreh,srehn,                                      &
                       li,lin,                                          &
                       cape,capen,                                      &
                       brn,brnn,                                        &
                       cin,cinn,                                        &
                       llmc,llmcn,                                      &
                       pw,pwn,                                          &
                       cmpref,cmprefn)
!-----------------------------------------------------------------------
!   Find the information about the time and date of the datafile
!   and add the lead time (minohour) to it, find the day of week
!  read (filename(1:26),'(2x,I4,I2,I2,I2,8x,I6)')
!    :        year1, month1,dateomon,houroday,minohour
!  print *, filename(1:26),'DATE VARIABLES'
    PRINT *, year1,month1,dateomon,houroday,minohour
    CALL calender(year1,month1,dateomon,houroday,minohour,dayoweek)
!  print *, year1, month1,dateomon,houroday,minohour,dayoweek
!  print *, 'year1, month1,dateomon,houroday,minohour,dayoweek'
!   Write out the output 2-D fields
    WRITE (0,*) 'About to call EnsWrtAll'
    filename=outfile(nf)
    lenfil = len_trim(filename)
    WRITE(6,'(a,a/)')                                                   &
        ' The file name of output set is ', filename(1:lenfil)
    CALL enswrtall (filename(1:lenfil), 2, dxnew,                       &
                    houroday,minohour,dayoweek,dateomon,month1,year1,   &
                    mslpn, hgt250n, vort250n, uwind250n, vwind250n,     &
                    hgt500n, vort500n, uwind500n, vwind500n,            &
                    hgt850n, vort850n, uwind850n, vwind850n,            &
                    temp850n, thckn, rh700n, omega700n,                 &
                    temp2mn, dewp2mn, accpptn, convpptn,                &
                    srehn, lin, capen, brnn, cinn, llmcn,pwn,cmprefn)

!-----------------------------------------------------------------------
!  Plotting of the relavent parameters if iplot=1
!  Plotting of the terrain contour if iterr=1
    IF (iterr == 1) THEN
      CALL a2plt(cmprefn,x1n,y1n,qni,qnj,dxnew,100.0,'HTER')
    ELSE IF (iplot == 1) THEN
      CALL a2plt(accpptn,x1n,y1n,qni,qnj,dxnew,1.0,'TPPT')
      CALL a2plt(convpptn,x1n,y1n,qni,qnj,dxnew,1.0,'CPPT')
    ELSE IF (iplot == 2) THEN
      CALL a2plt(mslpn,x1n,y1n,qni,qnj,dxnew,200.0,'MSLP')
      CALL a2plt(hgt250n,x1n,y1n,qni,qnj,dxnew,50.0,'H250')
      CALL a2plt(hgt500n,x1n,y1n,qni,qnj,dxnew,20.0,'H500')
      CALL a2plt(hgt850n,x1n,y1n,qni,qnj,dxnew,20.0,'H850')
      CALL a2plt(vort250n,x1n,y1n,qni,qnj,dxnew,0.00002,'V250')
      CALL a2plt(vort500n,x1n,y1n,qni,qnj,dxnew,0.00002,'V500')
      CALL a2plt(vort850n,x1n,y1n,qni,qnj,dxnew,0.00002,'V850')
      CALL a2plt(temp850n,x1n,y1n,qni,qnj,dxnew,2.0,'T850')
      CALL a2plt(thckn,x1n,y1n,qni,qnj,dxnew,20.0,'THCK')
      CALL a2plt(rh700n,x1n,y1n,qni,qnj,dxnew,10.0,'RH70')
      CALL a2plt(omega700n,x1n,y1n,qni,qnj,dxnew,2.0,'OM70')
      CALL a2plt(temp2mn,x1n,y1n,qni,qnj,dxnew,2.0,'TC2M')
      CALL a2plt(dewp2mn,x1n,y1n,qni,qnj,dxnew,2.0,'TD2M')
      CALL a2plt(accpptn,x1n,y1n,qni,qnj,dxnew,2.0,'TPPT')
      CALL a2plt(convpptn,x1n,y1n,qni,qnj,dxnew,2.0,'CPPT')
      CALL a2plt(srehn,x1n,y1n,qni,qnj,dxnew,50.0,'SREH')
      CALL a2plt(lin,x1n,y1n,qni,qnj,dxnew,2.0,'LIDX')
      CALL a2plt(capen,x1n,y1n,qni,qnj,dxnew,500.0,'CAPE')
      CALL a2plt(brnn,x1n,y1n,qni,qnj,dxnew,10.0,'BRN ')
      CALL a2plt(cinn,x1n,y1n,qni,qnj,dxnew,50.0,'CIN ')
      CALL a2plt(llmcn,x1n,y1n,qni,qnj,dxnew,0.0001,'LLMC')
      CALL a2plt(pwn,x1n,y1n,qni,qnj,dxnew,2000.0,' PW ')
      CALL a2plt(cmprefn,x1n,y1n,qni,qnj,dxnew,2000.0,'CPRF')
    END IF

    7000  CONTINUE

  END DO
  WRITE(6,'(a)') 'Program completed.'
  IF (iplot >= 1) CALL xgrend
  STOP
END PROGRAM arpsenscv


SUBROUTINE a2plt(a,x,y,nx,ny,dx,zinc,title) 26
  IMPLICIT NONE
  INTEGER :: nx,ny,ncl,mode
  INTEGER :: i,j
  REAL :: a(nx,ny),x(nx,ny),y(nx,ny),dx
  INTEGER :: iw(212,152)
  REAL :: cl(100),zmax,zmin,zinc
  CHARACTER (LEN=4) :: title
  CALL xpspac(0.1,0.9,0.05,0.15+80.0/116.0*0.8)
!  call XMAP(0.0,2400.0,0.0,2080.0)
!  call XMAP(0.0,3480.0,0.0,2400.0)
  CALL xmap(0.0,(nx-1)*dx,0.0,(ny-1)*dx)
  CALL xaxsca(0.0,(nx-1)*dx,dx,0.0,(ny-1)*dx,dx)
!  call XAXSCA(0.0,3480.0,30.0,0.0,2400.0,30.0)
!  call XAXSCA(0.0,2400.0,32.0,0.0,2080.0,32.0)
  CALL xbordr

  ncl=100
  zmax=-1.0E-10
  zmin=+1.0E+10
  DO j=2,ny-1
    DO i=2,nx-1
      IF (a(i,j) > zmax) zmax=a(i,j)
      IF (a(i,j) < zmin) zmin=a(i,j)
    END DO
  END DO

  IF (ABS(zmax-zmin) < 1.0E-8) GO TO 200
  cl(1)=INT(zmin/zinc)*zinc
  cl(2)=cl(1)+zinc
  mode=1
  CALL xconta(  zmax=cl(ncl)
  zmin=cl(1)
  zinc=cl(MIN(2,ncl))-cl(1)
  200 CONTINUE

  PRINT *,zmax,zmin,zinc
  CALL xclimt(zmax,zmin,zinc,0.0)
  CALL xchsiz(80.00)
  CALL xcharc(1500.,-100.,title)
  CALL xframe
  RETURN
END SUBROUTINE a2plt
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc     C


SUBROUTINE pptsto(ni,nj,a,b) 2
  INTEGER :: ni,nj,i,j
  REAL :: a(ni,nj),b(ni,nj),c
  DO j=1,nj
    DO i=1,ni
      c=a(i,j)-b(i,j)
      b(i,j)=a(i,j)
      a(i,j)=c
    END DO
  END DO
  RETURN
END SUBROUTINE pptsto
!


SUBROUTINE ary2dcv(nx,ny,a,ni,nj,b) 26
  IMPLICIT NONE
  INTEGER :: nx,ny,ni,nj
  REAL :: a(nx,ny),b(ni,nj)
  INTEGER :: i,j
!  print *,a(35,35),a(78,68)
  DO j=1,nj
    DO i=1,ni
      b(i,j)=a(i,j)
    END DO
  END DO
  RETURN
END SUBROUTINE ary2dcv


SUBROUTINE vrtcalc(u,v,ni,nj,b,f,x,y) 3,1
  IMPLICIT NONE
  INTEGER :: ni,nj
  REAL :: u(ni,nj),v(ni,nj),b(ni,nj),f(ni,nj),x(ni,nj),y(ni,nj)
  INTEGER :: i,j
  DO j=1,nj
    DO i=1,ni
      b(i,j)=0.0
    END DO
  END DO
  DO j=2,nj-1
    DO i=2,ni-1
      b(i,j)=(v(i+1,j)-v(i-1,j))/(x(i+1,j)-x(i-1,j))                    &
            -(u(i,j+1)-u(i,j-1))/(y(i,j+1)-y(i,j-1))                    &
            +f(i,j)
    END DO
  END DO
  CALL edgfill(b,1,ni,2,ni-1,1,nj,2,nj-1,1,1,1,1)
  RETURN
END SUBROUTINE vrtcalc


SUBROUTINE enswrthd (filenm, fileun,                                    & 1,3
           enstag, instit, model, resltn, resun,                        &
           houroday, minohour, dayoweek, dateomon,                      &
           month, year)
!
!======================================================================-
!             &&&&    G E N E R A L    I N F O R M A T I O N    &&&&
!======================================================================-
!
!  PURPOSE: This output subroutine produces the header for SAMEX ASCII
!           fields for conversion to images.
!
!  AUTHOR:  Henry Neeman, Univ. of Oklahoma, Project SAMEX
!           (hneeman@tornado.caps.ou.edu)
!
!  HISTORY:
!        1998/03/18  First written [HJN]
!
!  INPUT ARGUMENTS:
!        filenm:   string for filename to output data to
!        fileun:   integer for Fortran file unit
!        enstag:   tag for ensemble (e.g., 'SAMEX')
!        instit:   string for name of institution
!                    (e.g., 'CAPS', 'NCEP', 'NCAR', 'AFWA', 'NSSL')
!        model:    string for name of model
!                    (e.g., 'ARPS4.3.3', 'MM5V2')
!          Note: the model name should not contain any blanks or tabs.
!        resltn:   real for resolution (e.g., 30 for 30 km resolution)
!        resun:    string for name of resolution units (e.g., 'km')
!        houroday: integer for hour of the day (e.g., 20)
!        minohour: integer for minute of the hour (e.g., 30)
!        dayoweek: integer for day the week
!                    1=Sun, 2=Mon, 3=Tue, 4=Wed, 5=Thu, 6=Fri, 7=Sat
!        dateomon: integer for date of the month
!        month:    integer for month (1=Jan, 2=Feb, etc.)
!        year:     integer for year (e.g., 1998)
!
!  OUTPUT:
!        ASCII file header for 2D data field
!
!  I/O:
!        Output data to a Fortran unformatted ASCII file.
!
!  SPECIAL REQUIREMENTS:
!        <none>
!
!  OTHER INFORMATION:
!        <none>
!
!======================================================================-
!                    %%%%    D E C L A R A T I O N S    %%%%
!======================================================================-
!
  IMPLICIT NONE
! - - - - - - - - -  Include files - - - - - - - - - - - - - - - - - - -
!
! - - - - - - - - -  Constant declarations - - - - - - - - - - - - - - -
!
! - - - - - - - - -  Argument declarations - - - - - - - - - - - - - - -
  CHARACTER (LEN=*) :: filenm
  INTEGER :: fileun
  CHARACTER (LEN=*) :: enstag, instit, model
  REAL :: resltn
  CHARACTER (LEN=*) :: resun
  INTEGER :: houroday, minohour, dayoweek, dateomon, month, year
  INTEGER :: lenstag, linstit, lmodel
! - - - - - - - - -  Global/External declarations  - - - - - - - - - - -
!
! - - - - - - - - -  Local declarations  - - - - - - - - - - - - - - - -
  INTEGER :: iost
! - - - - - - - - -  Data statements - - - - - - - - - - - - - - - - - -
!
! - - - - - - - - -  Statement functions - - - - - - - - - - - - - - - -
!
!======================================================================-
!                 @@@@    E X E C U T A B L E    C O D E    @@@@
!======================================================================-
!
  100 FORMAT (a)
  110 FORMAT (a, ' ', a, ' ', a, ' ', g10.3, ' ', a)
  120 FORMAT (i2, i3, i2, i3, i3, i5)
  lenstag = 0
  CALL strlnth( enstag, lenstag)
  linstit = 0
  CALL strlnth( instit, linstit)
  lmodel = 0
  CALL strlnth( model, lmodel)
  OPEN (UNIT=fileun,IOSTAT=iost,ERR=150,FILE=filenm)
  WRITE (fileun,*) enstag(1:lenstag)
  WRITE (fileun,*) instit(1:linstit), ' ', model(1:lmodel),             &
                   ' ', resltn, ' ', resun
  WRITE (fileun,*) houroday, ' ', minohour, ' ', dayoweek, ' ',         &
                   dateomon, ' ', month, ' ', year
  RETURN
  150 WRITE (0,160) filenm
  RETURN
  160 FORMAT ('Error: cannot open file ', a,                            &
              ' for writing ensemble header!')
END SUBROUTINE enswrthd

!
!
!                   ########################################
!                   ########################################
!                   ########################################
!                   ########                        ########
!                   ########        ENSEMBLE        ########
!                   ########                        ########
!                   ########################################
!                   ########################################
!                   ########################################
!
!


SUBROUTINE enswrtq (filenm, fileun,                                     & 29
           qtag, qunits,                                                &
           qswcorn, qswcorw, qnecorn, qnecorw,                          &
           qstag,                                                       &
           qni, qnj, qfield,                                            &
           islast)
!
!======================================================================-
!             &&&&    G E N E R A L    I N F O R M A T I O N    &&&&
!======================================================================-
!
!  PURPOSE: This output subroutine produces each SAMEX ASCII quantity
!           field for conversion to images.
!
!  AUTHOR:  Henry Neeman, Univ. of Oklahoma, Project SAMEX
!           (hneeman@tornado.caps.ou.edu)
!
!  HISTORY:
!        1998/03/18  First written [HJN]
!
!  INPUT ARGUMENTS:
!        filenm: string for filename to output data to
!        fileun:  integer for Fortran file unit
!        qtag: string for quantity tag (e.g., 'MSLP')
!        qunits: string for quantity units (e.g., 'km')
!        qswcorn: real for quantity southwest corner north coordinate
!        qswcore: real for quantity southwest corner west  coordinate
!        qnecorn: real for quantity northeast corner north coordinate
!        qnecore: real for quantity northeast corner west  coordinate
!        qstag: string for quantity staggering (for now, CENTER)
!        qni: integer for quantity number of cells along the major axis
!        qnj: integer for quantity number of cells along the minor axis
!        qfield:  real array for the quantity field data
!        islast:  logical for the quantity being the last in the file
!
!  OUTPUT:
!        ASCII file for 2D data field
!
!  I/O:
!        Output data to a Fortran unformatted ASCII file.
!
!  SPECIAL REQUIREMENTS:
!        <none>
!
!  OTHER INFORMATION:
!        <none>
!
!======================================================================-
!                    %%%%    D E C L A R A T I O N S    %%%%
!======================================================================-
!
  IMPLICIT NONE
! - - - - - - - - -  Include files - - - - - - - - - - - - - - - - - - -
!
! - - - - - - - - -  Constant declarations - - - - - - - - - - - - - - -
!
! - - - - - - - - -  Argument declarations - - - - - - - - - - - - - - -
  CHARACTER (LEN=*) :: filenm
  INTEGER :: fileun
  CHARACTER (LEN=*) :: qtag, qunits
  REAL :: qswcorn, qswcorw, qnecorn, qnecorw
  CHARACTER (LEN=*) :: qstag
  INTEGER :: qni, qnj
  REAL :: qfield(qni,qnj)
  LOGICAL :: islast
! - - - - - - - - -  Global/External declarations  - - - - - - - - - - -
!
! - - - - - - - - -  Local declarations  - - - - - - - - - - - - - - - -
  INTEGER :: iost
  INTEGER :: i, j
! - - - - - - - - -  Data statements - - - - - - - - - - - - - - - - - -

! - - - - - - - - -  Statement functions - - - - - - - - - - - - - - - -

!======================================================================-
!                 @@@@    E X E C U T A B L E    C O D E    @@@@
!======================================================================-
!
  200 FORMAT (a, ' ', a)
  210 FORMAT (g10.5, ' N ', g10.5, ' W ', g10.5, ' N ', g10.5, ' W ')
  220 FORMAT (a)
  230 FORMAT (i5, ' ', i5)
  240 FORMAT (e16.8, e16.8, e16.8, e16.8, e16.8)
  WRITE (fileun,*) qtag, ' ', qunits
  WRITE (fileun,*) qswcorn, ' N ', qswcorw, ' W ',                      &
                   qnecorn, ' N ', qnecorw, ' W'
  WRITE (fileun,*) qstag
  WRITE (fileun,*) qni, ' ', qnj
  DO j = 1, qnj
    WRITE (fileun,240) (qfield(i,j),i=1,qni)
  END DO
  IF (islast) THEN
    CLOSE (UNIT=fileun,IOSTAT=iost,ERR=250,STATUS='KEEP')
  ELSE
    WRITE (fileun,*) ' '
  END IF
  RETURN
  250 WRITE (0,260) filenm
  RETURN
  260 FORMAT ('Error: cannot close file ', a,                           &
              ' after writing ensemble data!')
END SUBROUTINE enswrtq
!
!
!                   ########################################
!                   ########################################
!                   ########################################
!                   ########                        ########
!                   ########        ENSEMBLE        ########
!                   ########                        ########
!                   ########################################
!                   ########################################
!                   ########################################
!
!


SUBROUTINE enswrtall (filenm, fileun, resltn,                           & 1,30
           houroday, minohour, dayoweek, dateomon,                      &
           month, year,                                                 &
           mslp, hgt250, vort250, uwind250, vwind250,                   &
           hgt500, vort500, uwind500, vwind500,                         &
           hgt850, vort850, uwind850, vwind850,                         &
           temp850, thck, rh700, omega700,                              &
           temp2m, dewp2m, accppt, convppt,                             &
           sreh, li, cape, brn, cin, llmc,pw,cmpref)
!
!======================================================================-
!             &&&&    G E N E R A L    I N F O R M A T I O N    &&&&
!======================================================================-
!
!  PURPOSE: This output subroutine produces the header for SAMEX ASCII
!           fields for conversion to images.
!
!  AUTHOR:  Henry Neeman, Univ. of Oklahoma, Project SAMEX
!           (hneeman@tornado.caps.ou.edu)
!
!  HISTORY:
!        1998/03/18  First written [HJN]
!
!  INPUT ARGUMENTS:
!        filenm:   string for filename to output data to
!        fileun:   integer for Fortran file unit
!          Note: the model name should not contain any blanks or tabs.
!        resltn:   real for resolution (e.g., 30 for 30 km resolution)
!        houroday: integer for hour of the day (e.g., 20)
!        minohour: integer for minute of the hour (e.g., 30)
!        dayoweek: integer for day the week
!                    1=Sun, 2=Mon, 3=Tue, 4=Wed, 5=Thu, 6=Fri, 7=Sat
!        dateomon: integer for date of the month
!        month:    integer for month (1=Jan, 2=Feb, etc.)
!        year:     integer for year (e.g., 1998)
!        mslp:     real array for mean sea level pressure in pascals
!        hgt250:   real array for 250 mb height in m
!        vort250:  real array for 250 mb absolute vorticity in 1/s
!        uwing250: real array for 250 mb u wind velocity in m/s
!        vwing250: real array for 250 mb v wind velocity in m/s
!        hgt500:   real array for 500 mb height in m
!        vort500:  real array for 500 mb absolute vorticity in 1/s
!        uwing500: real array for 500 mb u wind velocity in m/s
!        vwing500: real array for 500 mb v wind velocity in m/s
!        hgt850:   real array for 850 mb height in m
!        vort850:  real array for 850 mb absolute vorticity in 1/s
!        uwing850: real array for 850 mb u wind velocity in m/s
!        vwing850: real array for 850 mb v wind velocity in m/s
!        temp850:  real array for 850 mb temperature in degrees C
!        thck:     real array for 1000-500 mb thickness in m
!        rh700:    real array for 700 mb relative humidity in %
!        omega700: real array for 700 mb omega vertical velocity in -microb/s
!        temp2m:   real array for 2 m temperature in degrees C
!        dewp2m:   real array for 2 m dewpoint in degrees C
!        accppt:   real array for accumulated precipitation in mm
!        convppt:  real array for convective precipitation in mm
!        sreh:     real array for storm-relative environmental helicity
!                    in m/s^2
!        li:       real array for lifted index in degrees C
!        cape:     real array for convective available potential energy in J/kg
!        brn:      real array for bulk Richardson number (unitless)
!        cin:      real array for convective inhibition in J/kg
!        llmc:     real array for low-level moisture convergence in g/kg/s
!        pw:       real array for precipitable water in g/m^2
!
!  OUTPUT:
!        ASCII file header for 2D data field
!
!  I/O:
!        Output data to a Fortran unformatted ASCII file.
!
!  SPECIAL REQUIREMENTS:
!        <none>
!
!  OTHER INFORMATION:
!        <none>
!
!======================================================================-
!                    %%%%    D E C L A R A T I O N S    %%%%
!======================================================================-
!
  IMPLICIT NONE
! - - - - - - - - -  Include files - - - - - - - - - - - - - - - - - - -
  INCLUDE 'enscv.inc'
! - - - - - - - - -  Constant declarations - - - - - - - - - - - - - - -
!
! - - - - - - - - -  Argument declarations - - - - - - - - - - - - - - -
  CHARACTER (LEN=*) :: filenm
  INTEGER :: fileun
  REAL :: resltn
  INTEGER :: houroday, minohour, dayoweek, dateomon, month, year
  REAL :: mslp(qni,qnj)
  REAL :: hgt250(qni,qnj), vort250(qni,qnj),                            &
       uwind250(qni,qnj), vwind250(qni,qnj)
  REAL :: hgt500(qni,qnj), vort500(qni,qnj),                            &
       uwind500(qni,qnj), vwind500(qni,qnj)
  REAL :: hgt850(qni,qnj), vort850(qni,qnj),                            &
       uwind850(qni,qnj), vwind850(qni,qnj)
  REAL :: temp850(qni,qnj)
  REAL :: thck(qni,qnj), rh700(qni,qnj), omega700(qni,qnj),             &
       temp2m(qni,qnj), dewp2m(qni,qnj)
  REAL :: accppt(qni,qnj), convppt(qni,qnj)
  REAL :: sreh(qni,qnj), li(qni,qnj), cape(qni,qnj), brn(qni,qnj),      &
       cin(qni,qnj), llmc(qni,qnj), pw(qni,qnj)
  REAL :: cmpref(qni,qnj)
! - - - - - - - - -  Global/External declarations  - - - - - - - - - - -
!
! - - - - - - - - -  Local declarations  - - - - - - - - - - - - - - - -
!
! - - - - - - - - -  Data statements - - - - - - - - - - - - - - - - - -
!
! - - - - - - - - -  Statement functions - - - - - - - - - - - - - - - -
!
!======================================================================-
!                 @@@@    E X E C U T A B L E    C O D E    @@@@
!======================================================================-
!
  PRINT *, 'starting (inside) EnsWrt'
  CALL enswrthd (filenm, fileun,                                        &
                 enstag, instit, model, resltn, resun,                  &
                 houroday, minohour, dayoweek, dateomon,                &
                 month, year)
  CALL enswrtq (filenm, fileun,                                         &
                mslptag, mslpunits,                                     &
                qswcorn, qswcorw, qnecorn, qnecorw,                     &
                qstag,                                                  &
                qni, qnj, mslp,                                         &
               .false.)
  CALL enswrtq (filenm, fileun,                                         &
                hgt250tag, hgt250units,                                 &
                qswcorn, qswcorw, qnecorn, qnecorw,                     &
                qstag,                                                  &
                qni, qnj, hgt250,                                       &
                .false.)
  CALL enswrtq (filenm, fileun,                                         &
                vort250tag, vort250units,                               &
                qswcorn, qswcorw, qnecorn, qnecorw,                     &
                qstag,                                                  &
                qni, qnj, vort250,                                      &
                .false.)
  CALL enswrtq (filenm, fileun,                                         &
                uwind250tag, uwind250units,                             &
                qswcorn, qswcorw, qnecorn, qnecorw,                     &
                qstag,                                                  &
                qni, qnj, uwind250,                                     &
                .false.)
  CALL enswrtq (filenm, fileun,                                         &
                vwind250tag, vwind250units,                             &
                qswcorn, qswcorw, qnecorn, qnecorw,                     &
                qstag,                                                  &
                qni, qnj, vwind250,                                     &
                .false.)
  CALL enswrtq (filenm, fileun,                                         &
                hgt500tag, hgt500units,                                 &
                qswcorn, qswcorw, qnecorn, qnecorw,                     &
                qstag,                                                  &
                qni, qnj, hgt500,                                       &
                .false.)
  CALL enswrtq (filenm, fileun,                                         &
                vort500tag, vort500units,                               &
                qswcorn, qswcorw, qnecorn, qnecorw,                     &
                qstag,                                                  &
                qni, qnj, vort500,                                      &
                .false.)
  CALL enswrtq (filenm, fileun,                                         &
                uwind500tag, uwind500units,                             &
                qswcorn, qswcorw, qnecorn, qnecorw,                     &
                qstag,                                                  &
                qni, qnj, uwind500,                                     &
                .false.)
  CALL enswrtq (filenm, fileun,                                         &
                vwind500tag, vwind500units,                             &
                qswcorn, qswcorw, qnecorn, qnecorw,                     &
                qstag,                                                  &
                qni, qnj, vwind500,                                     &
                .false.)
  CALL enswrtq (filenm, fileun,                                         &
                hgt850tag, hgt850units,                                 &
                qswcorn, qswcorw, qnecorn, qnecorw,                     &
                qstag,                                                  &
                qni, qnj, hgt850,                                       &
                .false.)
  CALL enswrtq (filenm, fileun,                                         &
                vort850tag, vort850units,                               &
                qswcorn, qswcorw, qnecorn, qnecorw,                     &
                qstag,                                                  &
                qni, qnj, vort850,                                      &
                .false.)
  CALL enswrtq (filenm, fileun,                                         &
                uwind850tag, uwind850units,                             &
                qswcorn, qswcorw, qnecorn, qnecorw,                     &
                qstag,                                                  &
                qni, qnj, uwind850,                                     &
                .false.)
  CALL enswrtq (filenm, fileun,                                         &
                vwind850tag, vwind850units,                             &
                qswcorn, qswcorw, qnecorn, qnecorw,                     &
                qstag,                                                  &
                qni, qnj, vwind850,                                     &
                .false.)
  CALL enswrtq (filenm, fileun,                                         &
                temp850tag, temp850units,                               &
                qswcorn, qswcorw, qnecorn, qnecorw,                     &
                qstag,                                                  &
                qni, qnj, temp850,                                      &
                .false.)
  CALL enswrtq (filenm, fileun,                                         &
                thcktag, thckunits,                                     &
                qswcorn, qswcorw, qnecorn, qnecorw,                     &
                qstag,                                                  &
                qni, qnj, thck,                                         &
                .false.)
  CALL enswrtq (filenm, fileun,                                         &
                rh700tag, rh700units,                                   &
                qswcorn, qswcorw, qnecorn, qnecorw,                     &
                qstag,                                                  &
                qni, qnj, rh700,                                        &
                .false.)
  CALL enswrtq (filenm, fileun,                                         &
                omega700tag, omega700units,                             &
                qswcorn, qswcorw, qnecorn, qnecorw,                     &
                qstag,                                                  &
                qni, qnj, omega700,                                     &
                .false.)
  CALL enswrtq (filenm, fileun,                                         &
                temp2mtag, temp2munits,                                 &
                qswcorn, qswcorw, qnecorn, qnecorw,                     &
                qstag,                                                  &
                qni, qnj, temp2m,                                       &
                .false.)
  CALL enswrtq (filenm, fileun,                                         &
                dewp2mtag, dewp2munits,                                 &
                qswcorn, qswcorw, qnecorn, qnecorw,                     &
                qstag,                                                  &
                qni, qnj, dewp2m,                                       &
                .false.)
  CALL enswrtq (filenm, fileun,                                         &
                accppttag, accpptunits,                                 &
                qswcorn, qswcorw, qnecorn, qnecorw,                     &
                qstag,                                                  &
                qni, qnj, accppt,                                       &
                .false.)
  CALL enswrtq (filenm, fileun,                                         &
                convppttag, convpptunits,                               &
                qswcorn, qswcorw, qnecorn, qnecorw,                     &
                qstag,                                                  &
                qni, qnj, convppt,                                      &
                .false.)
  CALL enswrtq (filenm, fileun,                                         &
                srehtag, srehunits,                                     &
                qswcorn, qswcorw, qnecorn, qnecorw,                     &
                qstag,                                                  &
                qni, qnj, sreh,                                         &
                .false.)
  CALL enswrtq (filenm, fileun,                                         &
                litag, liunits,                                         &
                qswcorn, qswcorw, qnecorn, qnecorw,                     &
                qstag,                                                  &
                qni, qnj, li,                                           &
                .false.)
  CALL enswrtq (filenm, fileun,                                         &
                capetag, capeunits,                                     &
                qswcorn, qswcorw, qnecorn, qnecorw,                     &
                qstag,                                                  &
                qni, qnj, cape,                                         &
                .false.)
  CALL enswrtq (filenm, fileun,                                         &
                brntag, brnunits,                                       &
                qswcorn, qswcorw, qnecorn, qnecorw,                     &
                qstag,                                                  &
                qni, qnj, brn,                                          &
                .false.)
  CALL enswrtq (filenm, fileun,                                         &
                cintag, cinunits,                                       &
                qswcorn, qswcorw, qnecorn, qnecorw,                     &
                qstag,                                                  &
                qni, qnj, cin,                                          &
                .false.)
  CALL enswrtq (filenm, fileun,                                         &
                llmctag, llmcunits,                                     &
                qswcorn, qswcorw, qnecorn, qnecorw,                     &
                qstag,                                                  &
                qni, qnj, llmc,                                         &
                .false.)
  CALL enswrtq (filenm, fileun,                                         &
                pwtag, pwunits,                                         &
                qswcorn, qswcorw, qnecorn, qnecorw,                     &
                qstag,                                                  &
                qni, qnj, pw,                                           &
                .false.)
  CALL enswrtq (filenm, fileun,                                         &
                crftag, crfunits,                                       &
                qswcorn, qswcorw, qnecorn, qnecorw,                     &
                qstag,                                                  &
                qni, qnj, cmpref,                                       &
                .true.)
  RETURN
END SUBROUTINE enswrtall
!
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE SECTHRZA                   ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE secthrza(nx,ny,nz,s,z,ss1) 3
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Interpolate 3-D data to a given horizontal level.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Ming Xue & Hao Jin
!  12/18/92.
!
!  MODIFICATION HISTORY:
!
!  9/10/94 (Weygandt & Y. Lu)
!  Cleaned up documentation.
!
!-----------------------------------------------------------------------
!

!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
!  INPUT :
!
!    nx       Number of grid points in the x-direction (east/west)
!    ny       Number of grid points in the y-direction (north/south)
!    nz       Number of grid points in the vertical
!
!    s        3-dimensional array of data to contour
!    s1       2-dimensional array of data to contour
!
!    x        x coordinate of grid points in physical/comp. space (m)
!    y        y coordinate of grid points in physical/comp. space (m)
!    z        z coordinate of grid points in physical space (m)
!
!    z01      value of x for first i grid point to plot
!
!  OUTPUT:
!    ss1      interpolated 3-D data to a given horizontal level
!
!-----------------------------------------------------------------------
!
!  Parameters of output
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  INTEGER :: nx,ny,nz          ! Number of grid points in 3 directions

  REAL :: s(nx,ny,nz)          ! 3-dimensional array of data to contour
  REAL :: z(nx,ny,nz)          ! z coordinate of grid points
                               ! in physical space (m)
  REAL :: ss1(nx,ny)           ! interpolated 3-D data to a
                               ! given horizontal level
!
!-----------------------------------------------------------------------
!
!  Common blocks for plotting control parameters
!
!-----------------------------------------------------------------------
!
  REAL :: z01                      ! the given height of the slice
  COMMON /sliceh/z01

!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------
!
  INTEGER :: i,j,k
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'globcst.inc'
!
!-----------------------------------------------------------------------
!
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!

!
!-----------------------------------------------------------------------
!
!  Find index for interpolation
!
!-----------------------------------------------------------------------
!
  DO j=1,ny-1
    DO i=1,nx-1
      IF(z01 <= z(i,j,1)) GO TO 11
      IF(z01 >= z(i,j,nz-1)) GO TO 12
      DO k=2,nz-2
        IF(z01 >= z(i,j,k).AND.z01 < z(i,j,k+1)) GO TO 15
      END DO

      11    k=1
      GO TO 15
      12    k=nz-1
      GO TO 15

      15    ss1(i,j)=s(i,j,k)+(s(i,j,k+1)-s(i,j,k))*                    &
                     (z01-z(i,j,k))/(z(i,j,k+1)-z(i,j,k))

!-----------------------------------------------------------------------
!
!  If the data point is below the ground level, set the
!  data value to the missing value.
!
!-----------------------------------------------------------------------

!  IF( z01.lt. z(i,j,2) ) ss1(i,j) = -9999.0

    END DO
  END DO

  RETURN
END SUBROUTINE secthrza


SUBROUTINE ptqv2m(nx,ny,nz, nstyps,zp,                                  & 1,5
           u ,v ,w ,ptprt, pprt, qvprt,                                 &
           ptbar, pbar, rhobar, qvbar,                                  &
           soiltyp,stypfrct,vegtyp,lai,roufns,veg,                      &
           tsfc,wetsfc,wetcanp,                                         &
           tem1,tem2, tem3,tem4,                                        &
           qvsfc,ptsfc,                                                 &
           vke2,ptke2,qvke2)

  IMPLICIT NONE

  INCLUDE 'phycst.inc'
  INCLUDE 'globcst.inc'
  INCLUDE 'sfcphycst.inc'

! input
  INTEGER :: nx,ny,nz,nstyps
  REAL :: zp(nx,ny,nz)
  REAL :: u(nx,ny,nz),v(nx,ny,nz),w(nx,ny,nz)
  REAL :: ptprt(nx,ny,nz),pprt(nx,ny,nz),qvprt(nx,ny,nz)
  REAL :: ptbar(nx,ny,nz),pbar(nx,ny,nz),qvbar(nx,ny,nz)
  REAL :: rhobar(nx,ny,nz)
  INTEGER :: soiltyp(nx,ny,nstyps)     ! Soil type
  REAL :: stypfrct(nx,ny,nstyps)    ! Soil type fraction
  INTEGER :: vegtyp (nx,ny) ! Vegetation type
  REAL :: lai    (nx,ny) ! Leaf Area Index
  REAL :: roufns (nx,ny) ! Surface roughness
  REAL :: veg    (nx,ny) ! Vegetation fraction

  REAL :: tsfc   (nx,ny,0:nstyps)    ! Soil Skin Temperature (K)
  REAL :: wetsfc (nx,ny,0:nstyps)    ! Surface soil moisture
  REAL :: wetcanp(nx,ny,0:nstyps)    ! Canopy water amount

! output
  REAL :: tem1(nx,ny),tem2(nx,ny)

!  temperary arrays
  REAL :: tem3(nx,ny),tem4(nx,ny)
  REAL :: ptsfc(nx,ny), qvsfc(nx,ny,0:nstyps)
  REAL :: vke2(nx,ny),ptke2(nx,ny),qvke2(nx,ny)

! temperary variables
  INTEGER :: i,j,k,is
  REAL :: z1,z1drou,z1droup
  REAL :: usuf,vsuf,wsuf,tsuf
  REAL :: psfc,qvsatts, rhgs, pterm, pi
  PARAMETER (pi=3.141592654)
  REAL :: gumove,gvmove
  DATA gumove /0.0/
  DATA gvmove /0.0/

! field moisture capacity
  REAL :: wfc(13)
  DATA wfc / .135, .150, .195, .255, .240, .255,                        &
             .322, .325, .310, .370, .367, 1.0E-25, 1.00/

  REAL :: f_qvsat

!fpp$ expand (f_qvsat)
!dir$ inline always f_qvsat

  PRINT *, 'starting PTQV2M'

  IF (nstyp <= 1) THEN
! for nstyp=1, the DTAREAD subroutine returns tsfc etc. with  (i,j,0)
! filled with data and (i,j,1) not assigned. This loop fills this Gap.
    DO j = 1, ny
      DO i = 1, nx
        tsfc(i,j,1)=tsfc(i,j,0)
        wetsfc(i,j,1)=wetsfc(i,j,0)
        wetcanp(i,j,1)=wetcanp(i,j,0)
        stypfrct(i,j,1)=1.0
      END DO
    END DO
  END IF

  DO j = 1, ny
    DO i = 1, nx
      qvsfc(i,j,0) = 0.0
      tem1(i,j)=0.0
      tem2(i,j)=0.0
      tem3(i,j)=0.0
      tem4(i,j)=0.0
    END DO
  END DO

!    deduce qvsfc(effctive qv at surface) ,adapted from INISFC
  PRINT *,'NSTYPS/NSTYP',nstyps,nstyp
  DO is=1,nstyp
    DO j = 1, ny-1
      DO i = 1, nx-1
        psfc = .5 * ( pbar(i,j,1) + pprt(i,j,1)                         &
                    + pbar(i,j,2) + pprt(i,j,2) )
        qvsatts = f_qvsat( psfc, tsfc(i,j,is) )
        IF ( soiltyp(i,j,is) == 12 .OR. soiltyp(i,j,is) == 13) THEN
                                                     ! Ice and water
          rhgs = 1.0
        ELSE
          pterm = .5 + SIGN( .5, wetsfc(i,j,is) - wfc(soiltyp(i,j,is)) )
          rhgs = pterm + (1.-pterm)                                     &
                       * 0.5 * ( 1. - COS( wetsfc(i,j,is) * pi          &
                                         / wfc(soiltyp(i,j,is)) ) )
        END IF
        qvsfc(i,j,is) = rhgs * qvsatts
      END DO
    END DO
  END DO

  DO is=1,nstyp
    DO j = 1, ny-1
      DO i = 1, nx-1
        qvsfc(i,j,0) = qvsfc(i,j,0) + qvsfc(i,j,is)*stypfrct(i,j,is)
      END DO
    END DO
  END DO

!  calculate wind speed,pt and qv at k=2,  pt at surface
  DO j = 1, ny-1
    DO i = 1, nx-1
      usuf=gumove+0.5*(u(i,j,2)+u(i+1,j,2))
      vsuf=gvmove+0.5*(v(i,j,2)+v(i,j+1,2))
      wsuf=       0.5*(w(i,j,2)+w(i,j,3))
      vke2(i,j)=MAX( SQRT(usuf*usuf+vsuf*vsuf+wsuf*wsuf),vsfcmin)
      ptke2(i,j)=ptprt(i,j,2)+ptbar(i,j,2)
      qvke2(i,j)=qvprt(i,j,2)+qvbar(i,j,2)
      psfc=0.5*(pprt(i,j,1)+pbar(i,j,1)                                 &
               +pprt(i,j,2)+pbar(i,j,2))
      ptsfc(i,j)=tsfc(i,j,0)*(100000.0/psfc)**rddcp
    END DO
  END DO

!  calculate C_u and C_pt for nuetral condition stored in tem1 and
!  tem2, respectively
  DO j = 1, ny-1
    DO i = 1, nx-1
      z1=0.5*(zp(i,j,3)-zp(i,j,2))
      z1=MIN(z1,z1limit)
      z1drou=z1/roufns(i,j)
      z1droup=(z1+roufns(i,j))/roufns(i,j)
      IF (z1droup < 0.0) THEN
        PRINT *,i,j,z1droup,z1,roufns(i,j)
        STOP
      END IF
      tem1(i,j)=kv/LOG(z1droup)
      tem2(i,j)=tem1(i,j)/prantl0l
      IF (soiltyp(i,j,1) == 13.AND.stypfrct(i,j,1) >= 0.5) THEN
        tem1(i,j)=SQRT((0.4+0.079*vke2(i,j))*1.e-3)
        tem2(i,j)=1.14E-3/tem1(i,j)
      END IF
    END DO
  END DO

!  calculate C_pt according to land/water and stability
  CALL cptca(nx,ny,nz,nx-1,ny-1,zp,roufns,vke2(1,1),                    &
            ptsfc(1,1),ptke2(1,1),tem2(1,1),tem3(1,1),1)
  CALL cptcwtra(nx,ny,nz,nx-1,ny-1,zp,soiltyp,vke2(1,1),                &
            ptsfc(1,1),ptke2(1,1),tem2(1,1),tem3(1,1),1)
  CALL cptca(nx,ny,nz,nx-1,ny-1,zp,roufns,vke2(1,1),                    &
            ptsfc(1,1),ptke2(1,1),tem2(1,1),tem4(1,1),2)
  CALL cptcwtra(nx,ny,nz,nx-1,ny-1,zp,soiltyp,vke2(1,1),                &
            ptsfc(1,1),ptke2(1,1),tem2(1,1),tem4(1,1),2)
  DO j = 1, ny-1
    DO i = 1, nx-1
      tem1(i,j)=ptsfc(i,j)+tem3(i,j)/tem4(i,j)*(ptke2(i,j)-ptsfc(i,j))
      tem2(i,j)=qvsfc(i,j,0)+tem3(i,j)/tem4(i,j)                        &
                                      *(qvke2(i,j)-qvsfc(i,j,0))
!  tem2(i,j)=qvke2(i,j)
    END DO
  END DO

  PRINT *, 'finishing PTQV2M'
  RETURN
END SUBROUTINE ptqv2m

!
!##################################################################
!##################################################################
!######                                                      ######
!######                   SUBROUTINE CPTC                    ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE cptca(nx,ny,nz,iend,jend,zp,roufns,wspd,ptsfc,pt1,           & 2
           c_ptneu,c_pt,i1)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Compute c_pt (a nondimensional temperature scale)
!
!  The quantity c_pt is used by the subroutine SFCFLX to obtain
!  surface fluxes for both the unstable and stable cases.
!
!-----------------------------------------------------------------------
!
!  AUTHORS: V. Wong, X. Song and N. Lin
!  9/10/1993
!
!  For the unstable case (Bulk Richardson number bulkri < 0 ), the
!  formulation is based on the paper "On the Analytical Solutions of
!  Flux-Profile relationships for the Atmospheric Surface Layer" by
!  D.W. Byun in J. of Applied Meteorlogy, July 1990, pp. 652-657.
!
!  For stable case, the formulation is based on the paper "A Short
!  History of the Operational PBL - Parameterization at ECMWF" by
!  J.F.Louis, M. Tiedtke and J.F. Geleyn in "Workshop on Planetary
!  Boundary Layer Parameterization", a publication by European Centre
!  for Medium Range Weather Forecasts, 25-27 Nov. 1981.
!
!  MODIFICATION HISTORY:
!
!  9/04/1994 (K. Brewster, V. Wong and X. Song)
!  Facelift.
!  2/27/95 (V. Wong and X. Song)
!
!  02/07/96 (V.Wong and X.Song)
!  Set an upper limiter, z1limit, for depth of the surface layer z1.
!
!  05/01/97 (V. Wong and X. Tan)
!  Changed the computation of temperature scale
!
!  05/29/97 (V. Wong and X. Tan)
!  Modified the formulation considering the height of the surface
!  layer z1 may equal zero.
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    nx       Number of grid points in the x-direction (east/west)
!    ny       Number of grid points in the y-direction (north/south)
!    nz       Number of grid points in the vertical
!
!    iend     i-index where evaluation ends.
!    jend     j-index where evaluation ends.
!
!    zp       The physical height coordinate defined at w-point of
!             staggered grid.
!    wspd     Surface wind speed (m/s), defined as
!             sqrt(usuf*usuf + vsuf*vsuf + wsuf*wsuf)
!    roufns   Surface roughness
!
!    ptsfc    Potential temperature at the ground level (K)
!    pt1      Potential temperature at the 1st scalar point above
!             ground level, (K)
!
!    c_ptneu  Temperature scale (K) at neutral state, defined by
!             surface heat flux / friction velocity
!
!  OUTPUT:
!
!    c_pt     Temperature scale (K), defined by
!             surface heat flux / friction velocity
!
!-----------------------------------------------------------------------
!

!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  INTEGER :: nx,ny,nz          ! The number grid points in 3 directions

  INTEGER :: iend              ! i-index where evaluation ends
  INTEGER :: jend              ! j-index where evaluation ends

  REAL :: zp    (nx,ny,nz)     ! The physical height coordinate
                               ! defined at w-point of staggered grid.
  REAL :: roufns(nx,ny)        ! Surface roughness length

  REAL :: wspd  (nx,ny)        ! Surface wind speed (m/s)

  REAL :: ptsfc(nx,ny)         ! Potential temperature at the
                               ! ground level (K)
  REAL :: pt1   (nx,ny)        ! Potential temperature at the
                               ! 1st scalar
                               ! point above ground level, (K)
  REAL :: c_ptneu(nx,ny)       ! Normalized temperature scale (K)
                               ! at neutral state

  REAL :: c_pt  (nx,ny)        ! Temperature scale (K)
!
!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------
!
  INTEGER :: i, j

  REAL :: z1                  ! The height of 1st scalar point above
                              ! ground level (m)
  REAL :: dz0                 ! z1-roufns, where roufns is defined as
                              ! surface roughness length
  REAL :: bulkri              ! Bulk Richardson number

  REAL :: stabp               ! Monin-Obukhov STABility Parameter
                              ! (zeta)

  REAL :: y1,y0,psih          ! Intermediate parameters needed
  REAL :: z1drou,qb3pb2
  REAL :: c7,c8
  REAL :: sb,qb,pb,thetab,tb  ! During  computations
  REAL :: a,b,c,d
  REAL :: tempan,sqrtqb

  REAL :: zt                  ! Thermal roughness length
  REAL :: dzt,ztdrou

  REAL :: z1droup,ztdroup
  REAL :: dz00
  INTEGER :: i1
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'phycst.inc'
  INCLUDE 'globcst.inc'
  INCLUDE 'sfcphycst.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Following Byun (1990).
!
!-----------------------------------------------------------------------
!
  PRINT *,'starting PTCP'
  DO j=1,jend
    DO i=1,iend

      dz00  = 0.5*(zp(i,j,3)-zp(i,j,2))
      dz00  = MIN(dz00,z1limit)
      dz00  = dz00 - roufns(i,j)
      bulkri = (g/ptsfc(i,j))*(pt1(i,j)-ptsfc(i,j))*dz00/               &
               (wspd(i,j)*wspd(i,j))

      IF (i1 == 1) THEN
        z1  = 0.5*(zp(i,j,3)-zp(i,j,2))
      ELSE
        z1=2.0
      END IF

      z1  = MIN(z1,z1limit)

      zt = ztz0*roufns(i,j)
      dzt = z1-zt
      ztdrou = z1/zt
      ztdroup = (z1+zt)/zt

      dz0 = z1-roufns(i,j)
      z1drou = z1/roufns(i,j)
      z1droup = (z1+roufns(i,j))/roufns(i,j)


      IF (bulkri <= 0.0) THEN
!
!-----------------------------------------------------------------------
!
!  Unstable case: See equations (28)-(34) in Byun (1990).
!
!-----------------------------------------------------------------------
!
        bulkri = MAX (bulkri,-10.0)

        sb =bulkri/prantl0l

        qb=oned9*(c1l+c2l*sb*sb)
        pb=oned54*(c3l+c4l*sb*sb)

        qb3pb2=qb**3-pb*pb
        c7 = (z1*dzt*LOG(z1droup)*LOG(z1droup))/(dz0*dz0*LOG(ztdroup))

        IF( qb3pb2 >= 0.0 ) THEN

          sqrtqb = SQRT(qb)
          tempan = MAX( -1.0, MIN( 1.0, pb/(sqrtqb**3) ) )

          thetab=ACOS(tempan)
          stabp =c7*(-2.0*SQRT(qb)*COS(thetab/3.0)+c5l)

        ELSE

          tb    =(SQRT(-qb3pb2)+ABS(pb))**oned3
          stabp =c7*(-(tb+qb/tb)+c5l)

        END IF
!
!-----------------------------------------------------------------------
!
!  According to equation (15) in Byun (1990).
!
!-----------------------------------------------------------------------
!
        c8=gammaml*stabp
        y1=SQRT(1.0 - c8)
        y0=SQRT(1.0 - c8/ztdrou)

        psih=2.0*LOG((y1+1.0)/(y0+1.0))
!
!-----------------------------------------------------------------------
!
!  Compute c_pt via equation (11) in Byun (1981).
!
!-----------------------------------------------------------------------
!
        c_pt(i,j)=kv / (prantl0l*(LOG(ztdroup)-psih))

      ELSE
!
!-----------------------------------------------------------------------
!
!  Stable case: See Louis et al (1981).
!
!-----------------------------------------------------------------------
!
        a=kv/LOG(ztdroup)
        b=5.0
        d=5.0
        c=SQRT(1.0+d*bulkri)

        c_pt(i,j) = SQRT(1.0+2.0*b*bulkri/c)
        c_pt(i,j) = a*c_pt(i,j)/(prantl0l*(1.0+3.0*b*bulkri*c))

      END IF

    END DO
  END DO

  RETURN
END SUBROUTINE cptca
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE CPTCWTR                    ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE cptcwtra(nx,ny,nz,iend,jend,zp,soiltyp,wspd,                 & 2
           ptsfc,pt1,c_ptwtrneu,c_ptwtr,i1)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Compute c_ptwtr (the product of ustar and ptstar) at sea case.
!  Note:  a temperature scale defined as surface heat flux divided
!  by the friction velocity.
!
!  The quantity c_ptwtr is used by the subroutine SFCFLX to obtain
!  surface fluxes for both the unstable and stable cases.
!
!-----------------------------------------------------------------------
!
!  AUTHORS: V. Wong and X. Song
!  8/04/1994
!
!  MODIFICATION HISTORY:
!
!  2/27/95 (V.W. and X.S.)
!
!  1/12/96 (V.W. and X.S.)
!  Changed the calculation related to zo over the sea.
!  Added kvwtr to denote the Von Karman constant over the sea;
!  Set a lower limiter for zo, zolimit, and an upper limiter for z1, z1limit.
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    nx       Number of grid points in the x-direction (east/west)
!    ny       Number of grid points in the y-direction (north/south)
!    nz       Number of grid points in the vertical
!
!    iend     i-index where evaluation ends.
!    jend     j-index where evaluation ends.
!
!    zp       The physical height coordinate defined at w-point of
!             staggered grid.
!    wspd     Surface wind speed (m/s), defined as
!             sqrt(usuf*usuf + vsuf*vsuf + wsuf*wsuf)
!    ptsfc    Potential temperature at the ground level (K)
!    pt1      Potential temperature at the 1st scalar point above
!             ground level, (K)
!
!
!  OUTPUT:
!
!    c_ptwtr  The product of ustar and ptstar. ptstar is temperature
!             scale (K), defined by surface heat flux / friction velocity
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  INTEGER :: nx,ny,nz          ! The number grid points in 3 directions

  INTEGER :: iend              ! i-index where evaluation ends
  INTEGER :: jend              ! j-index where evaluation ends

  REAL :: zp    (nx,ny,nz)     ! The physical height coordinate defined at
                               ! w-point of staggered grid.
  INTEGER :: soiltyp(nx,ny)    ! Soil type at each point.

  REAL :: wspd  (nx,ny)        ! Surface wind speed (m/s)

  REAL :: ptsfc(nx,ny)         ! Potential temperature at the ground level
                               ! (K)
  REAL :: pt1   (nx,ny)        ! Potential temperature at the 1st scalar
                               ! point above ground level, (K)

  REAL :: c_ptwtrneu(nx,ny)       ! Product of ustar and ptstar
  REAL :: c_ptwtr(nx,ny)       ! Product of ustar and ptstar
!
!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------
!
  INTEGER :: i, j

  REAL :: z1                   ! The height of 1st scalar point above
                               ! ground level (m)
  REAL :: zo                   ! Defined as surface momentum roughness
                               ! length
  REAL :: zt                   ! Defined as surface heat transfer
                               ! roughness length
  REAL :: dzo                  ! z1-zo
  REAL :: dzt                  ! z1-zt
  REAL :: z1dzo                ! z1/zo
  REAL :: z1dzt                ! z1/zt
  REAL :: xcdn                 ! Cdn (sea)
  REAL :: xcdh                 ! Hot-wired value of Cdh (sea)

  REAL :: bulkri               ! Bulk Richardson number
  REAL :: stabp                ! Monin-Obukhov STABility Parameter
                               ! (zeta)

  REAL :: y1,y0,psih           ! Intermediate parameters needed
  REAL :: qb3pb2
  REAL :: c7,c8
  REAL :: sb,qb,pb,thetab,tb   ! during  computations
  REAL :: a,b,c,d
  REAL :: tempan,sqrtqb
  REAL :: z10,zo0,zt0,dzo0,dzt0
  REAL :: i1
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'phycst.inc'
  INCLUDE 'globcst.inc'
  INCLUDE 'sfcphycst.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  PRINT *,'starting PTCPWTR'
  xcdh = 1.1E-3

  DO j=1,jend
    DO i=1,iend
      IF ( soiltyp(i,j) == 13) THEN
!  IF (soiltyp(i,j,1).eq.13.and.stypfrct(i,j,1).ge.0.5) THEN

!  calculate bulk Ri with the lowest level above the ground
        xcdn = (0.4+0.07*wspd(i,j)) * 1.e-3

        z10  = 0.5*(zp(i,j,3)-zp(i,j,2))
        z10  = MIN(z10,z1limit)
!
! calculate zo and zt
!
        zo0 =0.032*xcdn*wspd(i,j)*wspd(i,j)/9.8
        zo0 = MAX(zo0,zolimit)
        zt0 = z10 * EXP( -kv*SQRT(xcdn)/(prantl0w*xcdh) )

        dzo0 = z10-zo0
        dzt0 = z10-zt0

        bulkri =(g/ptsfc(i,j))*(pt1(i,j)-ptsfc(i,j))*dzo0*dzo0/         &
            (dzt0*wspd(i,j)*wspd(i,j))

!    fined Ri   bulk calcuation

        xcdn = (0.4+0.07*wspd(i,j)) * 1.e-3

        IF (i1 == 1) THEN
          z1  = 0.5*(zp(i,j,3)-zp(i,j,2))
        ELSE
          z1=2.0
        END IF
        z1  = MIN(z1,z1limit)
!
! calculate zo and zt
!
        zo =0.032*xcdn*wspd(i,j)*wspd(i,j)/9.8
        zo = MAX(zo,zolimit)
        zt = z1 * EXP( -kv*SQRT(xcdn)/(prantl0w*xcdh) )

        dzo = z1-zo
        z1dzo = z1/zo
        dzt = z1-zt
        z1dzt = z1/zt

        IF (bulkri <= 0.0) THEN
!
!-----------------------------------------------------------------------
!
!  Unstable case: A modified formulation, which is similar to
!  equations (28)-(34) in Byun (1990).
!
!-----------------------------------------------------------------------
!
          bulkri = MAX (bulkri,-10.0)

          sb =bulkri/prantl0w

          qb=oned9*(c1w+c2w*sb*sb)
          pb=oned54*(c3w+c4w*sb*sb)

          qb3pb2=qb**3-pb*pb
          c7 = (z1*dzt/(dzo*dzo))*(ALOG(z1dzo)*ALOG(z1dzo)/ALOG(z1dzt))

          IF( qb3pb2 >= 0.0 ) THEN

            sqrtqb = SQRT(qb)
            tempan = MAX( -1.0, MIN( 1.0, pb/(sqrtqb**3) ) )

            thetab=ACOS(tempan)
            stabp =c7*(-2.0*SQRT(qb)*COS(thetab/3.0)+c5w)

          ELSE

            tb    =(SQRT(-qb3pb2)+ABS(pb))**oned3
            stabp =c7*(-(tb+qb/tb)+c5w)

          END IF
!
!-----------------------------------------------------------------------
!
!  According to a modified equation, which is similar to equation (14)
!  in Byun (1990).
!
!-----------------------------------------------------------------------
!
          c8=gammamw*stabp
          y1=SQRT(1.0 - c8)
          y0=SQRT(1.0 - c8/z1dzt)

          psih=2.0*ALOG((y1+1.0)/(y0+1.0))
!
!-----------------------------------------------------------------------
!
!  Compute ptstar via equation (11) in Byun (1981).
!
!-----------------------------------------------------------------------
!
          c_ptwtr(i,j)=kv / (prantl0w*(ALOG(z1dzt)-psih))

        ELSE
!
!-----------------------------------------------------------------------
!
!  Stable case: With the modified formulation in Louis et al (1981).
!
!-----------------------------------------------------------------------
!
          a=kv*kv/(prantl0w*LOG(z1dzo)*LOG(z1dzt))
          b=5.0
          d=5.0
          c=SQRT(1.0+d*bulkri)

          c_ptwtr(i,j) = SQRT(1.0+2.0*b*bulkri/c)
          c_ptwtr(i,j) = a*c_ptwtr(i,j)/(prantl0l*(1.0+3.0*b*bulkri*c))

        END IF
      END IF

    END DO
  END DO

  RETURN
END SUBROUTINE cptcwtra


SUBROUTINE psfsl(nx,ny,t2,qv2,p2,z2,z1,psf,psl) 1
  IMPLICIT NONE
  INTEGER :: nx,ny
  REAL :: t2(nx,ny),qv2(nx,ny),p2(nx,ny)
  REAL :: z2(nx,ny),z1(nx,ny)
  REAL :: psf(nx,ny),psl(nx,ny) !output

  REAL :: gammam,rd,rv,tsh,zsh,rog,fvirt
  REAL :: tvsf,tvsl,tv2
  INTEGER :: i,j

  gammam=-6.5E-3
  rd=2.8705E2
  rv=4.615E2
  rog=rd/9.8
  fvirt=rv/rd-1.0
  zsh=75.0
  tsh=290.66

  DO j=1,ny-1
    DO i=1,nx-1
      tv2=t2(i,j)*(1.+fvirt*qv2(i,j))
      tvsf=tv2-gammam*(z2(i,j)-z1(i,j))
      psf(i,j)=p2(i,j)*EXP(2.*(z2(i,j)-z1(i,j))/(rog*(tvsf+tv2)))
      IF (z1(i,j) > zsh) THEN
        tvsl=tvsf-gammam*z1(i,j)
        IF (tvsl > tsh) THEN
          IF (tvsf > tsh) THEN
            tvsl=tsh-0.005*(tvsf-tsh)**2
          ELSE
            tvsl=tsh
          END IF
        END IF
      ELSE
        tvsl=tvsf
      END IF
      psl(i,j)=psf(i,j)*EXP(2.0*z1(i,j)/(rog*(tvsf+tvsl)))
      psf(i,j)=-ALOG(psf(i,j))
    END DO
  END DO
  RETURN
END SUBROUTINE psfsl
!


SUBROUTINE pintp(nx,ny,nz,s3d,nlgp3d,s2d,INDEX,nlgp01,                  & 15
           t3d,qv3d,z2,zsf,nlgpsf)
  IMPLICIT NONE
  INTEGER :: nx,ny,nz
  INTEGER :: i,j,k
  REAL :: s3d(nx,ny,nz),nlgp3d(nx,ny,nz),s2d(nx,ny)
  INTEGER :: INDEX
  REAL :: nlgp01
  REAL :: t3d(nx,ny,nz),qv3d(nx,ny,nz)
  REAL :: z2(nx,ny),zsf(nx,ny)
  REAL :: nlgpsf(nx,ny)

  REAL :: gammam,rd,rv,tsh,zsh,rog,fvirt,gammas
  REAL :: tvsf,tvsl,tv2
  REAL :: tvkm,part

  gammam=-6.5E-3
  rd=2.8705E2
  rv=4.615E2
  rog=rd/9.8
  fvirt=rv/rd-1.0
  zsh=75.0
  tsh=290.66

!  print *,'starting PINTP'
  DO j=1,ny-1
    DO i=1,nx-1
      IF(nlgp01 <= nlgp3d(i,j,2)) GO TO 11
      IF(nlgp01 >= nlgp3d(i,j,nz-1)) GO TO 12
      DO k=2,nz-2
        IF(nlgp01 >= nlgp3d(i,j,k).AND.nlgp01 < nlgp3d(i,j,k+1)) GO TO 15
      END DO

      11    k=2
      GO TO 15
      12    k=nz-2
      GO TO 15

      15    s2d(i,j)=s3d(i,j,k)+(s3d(i,j,k+1)-s3d(i,j,k))*              &
                   (nlgp01-nlgp3d(i,j,k))/(nlgp3d(i,j,k+1)-nlgp3d(i,j,k))
!  if (i.eq.1.and.j.eq.1) then
!  print *,i,j,k,nlgp01
!  print *,s2d(i,j),s3d(i,j,k),s3d(i,j,k+1),s3d(i,j,k),
!    :       nlgp01,nlgp3d(i,j,k),nlgp3d(i,j,k+1),nlgp3d(i,j,k)
!  stop
!  endif

      IF (nlgp01 <= nlgpsf(i,j)) THEN
        IF (INDEX == 1.OR.INDEX == 2) THEN
          s2d(i,j)=s3d(i,j,2)
        ELSE IF (INDEX == 3) THEN
          tvsf=t3d(i,j,2)*(1.+fvirt*qv3d(i,j,2))-gammam*(z2(i,j)-zsf(i,j))
          IF (zsf(i,j) > zsh) THEN
            tvsl=tvsf-gammam*zsf(i,j)
            IF (tvsl > tsh) THEN
              IF (tvsf > tsh) THEN
                tvsl=tsh-0.005*(tvsf-tsh)**2
              ELSE
                tvsl=tsh
              END IF
            END IF
            gammas=(tvsf-tvsl)/zsf(i,j)
          ELSE
            gammas=0.0
          END IF
          part=rog*(nlgp01-nlgpsf(i,j))
          s2d(i,j)=zsf(i,j)+tvsf*part/(1-0.5*gammas*part)
        END IF
      ELSE IF (nlgp01 > nlgpsf(i,j).AND.nlgp01 <= nlgp3d(i,j,2)) THEN
        IF (INDEX == 1) THEN
          s2d(i,j)=s3d(i,j,2)
        ELSE IF (INDEX == 2) THEN
          s2d(i,j)=s2d(i,j)
        ELSE IF (INDEX == 3) THEN
          s2d(i,j)=s3d(i,j,2)+(zsf(i,j)-s3d(i,j,2))*                    &
                 (nlgp01-nlgp3d(i,j,2))/(nlgpsf(i,j)-nlgp3d(i,j,2))
        END IF
      ELSE IF (nlgp01 > nlgp3d(i,j,nz-1)) THEN
        IF (INDEX == 1) THEN
          s2d(i,j)=s3d(i,j,nz-1)
        ELSE IF (INDEX == 2) THEN
          s2d(i,j)=s3d(i,j,nz-1)
        ELSE IF (INDEX == 3) THEN
          tvkm=t3d(i,j,nz-1)*(1.0+fvirt*qv3d(i,j,nz-1))
          s2d(i,j)=s3d(i,j,nz-1)+tvkm*rog*(nlgp01-nlgp3d(i,j,nz-1))
        END IF
      END IF
    END DO
  END DO

  RETURN
END SUBROUTINE pintp


SUBROUTINE cllmc(nx,ny,nz,u,v,qv,algpsf,zp,rhobar,algpzs,llmc,          & 1
           x,y, tem1,tem2)
  IMPLICIT NONE
  INTEGER :: nx,ny,nz
  REAL :: u(nx,ny,nz),v(nx,ny,nz),qv(nx,ny,nz)
  REAL :: algpzs(nx,ny,nz),zp(nx,ny,nz),rhobar(nx,ny,nz)
  REAL :: algpsf(nx,ny)
  REAL :: llmc(nx-1,ny-1)
  REAL :: tem1(nx,ny),tem2(nx,ny)
  REAL :: x(nx,ny),y(nx,ny)
  REAL :: ubar,vbar,qbar,mbar,dz
  REAL :: algpst,pst,psf
  INTEGER :: i,j,k,kt
  PRINT *,'starting CLLMC'
  DO j=1,ny-1
    DO i=1,nx-1
      llmc(i,j)=0
    END DO
  END DO
  DO j=1,ny-1
    DO i=1,nx-1
      psf=EXP(-algpsf(i,j))
      pst=psf-5000.0
      algpst=-ALOG(pst)
      ubar=0.0
      vbar=0.0
      qbar=0.0
      mbar=0.0
      DO k=2,nz-1
        IF (algpst < algpzs(i,j,k+1)) THEN
          kt=k
          GO TO 300
        END IF
      END DO
      PRINT *,'error, check CLLMC'
      STOP
      300   DO k=2,kt
        dz=zp(i,j,k+1)-zp(i,j,k)
        IF (k == kt) THEN
          dz=dz*(algpst-0.5*(algpzs(i,j,k)+algpzs(i,j,k-1)))            &
               /(0.5*(algpzs(i,j,k+1)-algpzs(i,j,k-1)))
        END IF
        ubar=ubar+u(i,j,k)*rhobar(i,j,k)*dz
        vbar=vbar+v(i,j,k)*rhobar(i,j,k)*dz
        qbar=qbar+qv(i,j,k)*rhobar(i,j,k)*dz
        mbar=mbar+rhobar(i,j,k)*dz
      END DO
      ubar=ubar/mbar
      vbar=vbar/mbar
      qbar=qbar/mbar
      tem1(i,j)=ubar*qbar
      tem2(i,j)=vbar*qbar
    END DO
  END DO

  DO j=2,ny-2
    DO i=2,nx-2
      llmc(i,j)=-1000.0*((tem1(i+1,j)-tem1(i-1,j))/(x(i+1,j)-x(i-1,j))  &
                        +(tem2(i,j+1)-tem2(i,j-1))/(y(i,j+1)-y(i,j-1)))
    END DO
  END DO

  PRINT *,'finishing CLLMC'
  RETURN
END SUBROUTINE cllmc


SUBROUTINE cape_ncep(nx,ny,nz,tk1,qv1,p1,cape,cins,li,                  & 1,1
           tk,qv,p,pint,plcl,peql,p1d,t1d,qv1d,lmh,                     &
           item1,item2,item3,item4,item5,item6,item7,   & ! work arrays
           rtem1,rtem2,rtem3,rtem4,rtem5,rtem6,rtem7)   ! work arrays
!
  IMPLICIT NONE
  INTEGER :: nx,ny,nz
  REAL :: tk1(nx,ny,nz),qv1(nx,ny,nz),p1(nx,ny,nz)
  REAL :: cape(nx,ny),cins(nx,ny),li(nx,ny)
  REAL :: tk(nx,ny,nz),qv(nx,ny,nz),p(nx,ny,nz),pint(nx,ny,nz)
  REAL :: plcl(nx,ny),peql(nx,ny)
  REAL :: p1d(nx,ny),t1d(nx,ny),qv1d(nx,ny)
  INTEGER :: lmh(nx,ny)

  INTEGER :: i,j,l,lk,lmhij,itype

  INTEGER :: item1(nx,ny),item2(nx,ny),item3(nx,ny),item4(nx,ny),       &
          item5(nx,ny),item6(nx,ny),item7(nx,ny)
  REAL :: rtem1(nx,ny),rtem2(nx,ny),rtem3(nx,ny),rtem4(nx,ny),          &
       rtem5(nx,ny),rtem6(nx,ny),rtem7(nx,ny)

!----------------------------------------------------------------------
!***********************************************************************
!--------------PREPARATIONS---------------------------------------------
  DO j=1,ny
    DO i=1,nx
      lmh(i,j)=nz-2
      lmhij=nz-2
      pint(i,j,1)=2500.
      DO l=1,lmhij
        lk=lmhij-l+2
        p(i,j,l)=p1(i,j,lk)
        tk(i,j,l)=tk1(i,j,lk)
        qv(i,j,l)=qv1(i,j,lk)
        pint(i,j,l+1)=(p1(i,j,lk)+p1(i,j,lk-1))*0.5
      END DO
    END DO
  END DO
!
  itype=1

  CALL calcape(itype,p1d,t1d,qv1d,                                      &
                     tk,qv,p,pint,lmh,nx,ny,nz-2,nz-1,                  &
                     cape,cins,li,plcl,peql,                            &
      item1,item2,item3,item4,item5,item6,item7,   & ! work arrays
  rtem1,rtem2,rtem3,rtem4,rtem5,rtem6,rtem7)   ! work arrays

!
!   WRITE(6,*) ' CAPE= ',CAPE(I,J)
!   WRITE(6,*) ' CIN= ',CINS(I,J)
!   WRITE(6,*) ' LCL= ',PLCL(I,J)
!   WRITE(6,*) ' EQ LVL= ',PEQL(I,J)

  RETURN
END SUBROUTINE cape_ncep


SUBROUTINE calcape(itype,p1d,t1d,q1d,                                   & 1,4
           t,q,p,pint,lmh,im,jm,lm,lp1,                                 &
           cape,cins,li,plcl,peql,                                      &
           ieql,iptb,ithtb,ilres,jlres,ihres,jhres,   & ! work arrays
           lcl,tpar,pp,qq,psp,thesp,ape)              ! work arrays
!$$$  SUBPROGRAM DOCUMENTATION BLOCK
!             .      .    .
! SUBPROGRAM:    CALCAPE     COMPUTES CAPE AND CINS
!   PRGRMMR: TREADON         ORG: W/NMC2     DATE: 93-02-10
!
! ABSTRACT:
!
!  THIS ROUTINE COMPUTES CAPE AND CINS GIVEN TEMPERATURE,
!  PRESSURE, AND SPECIFIC HUMIDTY.  IN "STORM AND CLOUD
!  DYNAMICS" (1989, ACADEMIC PRESS) COTTON AND ANTHES DEFINE
!  CAPE (EQUATION 9.16, P501) AS
!
!               EL
!      CAPE =  SUM G * LN(THETAP/THETAA) DZ
!              LCL
!
!  WHERE,
!   EL    = EQUILIBRIUM LEVEL,
!  LCL    = LIFTING CONDENSTATION LEVEL,
!    G    = GRAVITATIONAL ACCELERATION,
!  THETAP = LIFTED PARCEL POTENTIAL TEMPERATURE,
!  THETAA = AMBIENT POTENTIAL TEMPERATURE.
!
!  NOTE THAT THE INTEGRAND LN(THETAP/THETAA) APPROXIMATELY
!  EQUALS (THETAP-THETAA)/THETAA.  THIS RATIO IS OFTEN USED
!  IN THE DEFINITION OF CAPE/CINS.
!
!  TWO TYPES OF CAPE/CINS CAN BE COMPUTED BY THIS ROUTINE.  THE
!  SUMMATION PROCESS IS THE SAME FOR BOTH CASES.  WHAT DIFFERS

!  IS THE DEFINITION OF THE PARCEL TO LIFT.  FOR ITYPE=1 THE
!  PARCEL WITH THE WARMEST THETA-E IN A DPBND PASCAL LAYER ABOVE
!  THE MODEL SURFACE IS LIFTED.  THE ARRAYS P1D, T1D, AND Q1D
!  ARE NOT USED.  FOR ITYPE=2 THE ARRAYS P1D, T1D, AND Q1D
!  DEFINE THE PARCEL TO LIFT IN EACH COLUMN.  BOTH TYPES OF
!  CAPE/CINS MAY BE COMPUTED IN A SINGLE EXECUTION OF THE POST
!  PROCESSOR.
!
!  THIS ALGORITHM PROCEEDS AS FOLLOWS.
!  FOR EACH COLUMN,
!     (1)  INITIALIZE RUNNING CAPE AND CINS SUM TO 0.0
!     (2)  COMPUTE TEMPERATURE AND PRESSURE AT THE LCL USING
!          LOOK UP TABLE (PTBL).  USE EITHER PARCEL THAT GIVES
!          MAX THETAE IN LOWEST DPBND ABOVE GROUND (ITYPE=1)
!          OR GIVEN PARCEL FROM T1D,Q1D,...(ITYPE=2).
!     (3)  COMPUTE THE TEMP OF A PARCEL LIFTED FROM THE LCL.
!          WE KNOW THAT THE PARCEL'S
!          EQUIVALENT POTENTIAL TEMPERATURE (THESP) REMAINS
!          CONSTANT THROUGH THIS PROCESS.  WE CAN
!          COMPUTE TPAR USING THIS KNOWLEDGE USING LOOK
!          UP TABLE (SUBROUTINE TTBLEX).
!     (4)  FIND THE EQUILIBRIUM LEVEL.  THIS IS DEFINED AS THE
!          HIGHEST POSITIVELY BUOYANT LAYER.
!          (IF THERE IS NO POSITIVELY BUOYANT LAYER, CAPE/CINS
!           WILL BE ZERO)
!     (5)  COMPUTE CAPE/CINS.
!          (A) COMPUTE THETAP.  WE KNOW TPAR AND P.
!          (B) COMPUTE THETAA.  WE KNOW T AND P.
!     (6)  ADD G*(THETAP-THETAA)*DZ TO THE RUNNING CAPE OR CINS SUM.
!          (A) IF THETAP > THETAA, ADD TO THE CAPE SUM.
!          (B) IF THETAP < THETAA, ADD TO THE CINS SUM.
!     (7)  ARE WE AT EQUILIBRIUM LEVEL?
!          (A) IF YES, STOP THE SUMMATION.
!          (B) IF NO, CONTIUNUE THE SUMMATION.
!     (8)  ENFORCE LIMITS ON CAPE AND CINS (I.E. NO NEGATIVE CAPE)
!
! PROGRAM HISTORY LOG:
!   93-02-10  RUSS TREADON
!   93-06-19  RUSS TREADON - GENERALIZED ROUTINE TO ALLOW FOR
!                         TYPE 2 CAPE/CINS CALCULATIONS.
!   94-09-23  MIKE BALDWIN - MODIFIED TO USE LOOK UP TABLES
!                         INSTEAD OF COMPLICATED EQUATIONS.
!   94-10-13  MIKE BALDWIN - MODIFIED TO CONTINUE CAPE/CINS CALC
!                         UP TO AT HIGHEST BUOYANT LAYER.
!   98-03-25  MIKE BALDWIN - STAND ALONE VERSION FOR SAMEX
!
! USAGE:    CALL CALCAPE(ITYPE,P1D,T1D,Q1D,
!    &                   T,Q,P,PINT,LMH,IM,JM,LM,
!    &                   CAPE,CINS)
!   INPUT ARGUMENT LIST:
!  ITYPE,P1D,T1D,Q1D,T,Q,P,PINT,LMH,IM,JM,LM
!
!    CONTROL VARIABLES:
!    ITYPE    - INTEGER FLAG SPECIFYING HOW PARCEL TO LIFT IS
!               IDENTIFIED.  SEE COMMENTS ABOVE.
!    IM - 1ST HORIZONTAL DIMENSION OF ARRAY
!    JM - 2ND HORIZONTAL DIMENSION OF ARRAY
!    LM - VERTICAL DIMENSION OF ARRAY
!    LMH - ARRAY THAT CONTAINS THE NUMBER OF ABOVE GROUND LEVELS
!          AT EACH GRID POINT    DIMENSIONED (IM,JM)
!
!    3-D ARRAYS OF ATMOSPHERIC CONDITIONS:
!    MID-LAYER VARIABLES:
!    T - TEMP (K)         DIMENSIONED TO (IM,JM,LM)
!    Q - SPEC HUM (kg/kg) DIMENSIONED TO (IM,JM,LM)

!    P - PRESSURE (Pa)    DIMENSIONED TO (IM,JM,LM)
!
!    LAYER-INTERFACE VARIABLES:
!    PINT - PRESSURE (Pa) DIMENSIONED TO (IM,JM,LM+1)
!
!    2-D ARRAYS IF SENDING IN A SPECIFIC PARCEL TO LIFT (ITYPE=2):
!    P1D(IM,JM) - ARRAY OF PRESSURE (Pa) OF PARCELS TO LIFT.
!    T1D(IM,JM) - ARRAY OF TEMPERATURE (K) OF PARCELS TO LIFT.
!    Q1D(IM,JM) - ARRAY OF SPECIFIC HUMIDITY (kg/kg) OF PARCELS TO LIFT.
!
!   OUTPUT ARGUMENT LIST:
!  CAPE(IM,JM) - CONVECTIVE AVAILABLE POTENTIAL ENERGY (J/KG)
!  CINS(IM,JM) - CONVECTIVE INHIBITION (J/KG)
!  PLCL(IM,JM) - PRESSURE OF LCL (Pa)
!  PEQL(IM,JM) - PRESSURE OF EQ LEVEL (Pa)
!
!   OUTPUT FILES:
!  STDOUT  - RUN TIME STANDARD OUT.
!
!   SUBPROGRAMS CALLED:
!  UTILITIES:
!    TTBLEQ   - LOOKUP TABLE ROUTINE TO GET T FROM THETAE AND P
!    TTBLE1   - LOOKUP TABLE ROUTINE TO GET T FROM THETAE AND P
!    TTBLEX   - LOOKUP TABLE ROUTINE TO GET T FROM THETAE AND P
!
!
!   ATTRIBUTES:
!  LANGUAGE: FORTRAN
!  MACHINE : CRAY Y-MP
!$$$
!
!
!
!  INCLUDE/SET PARAMETERS.  CONSTANTS ARE FROM BOLTON (MWR, 1980).
!----------------------------------------------------------------------
!  PARAMETER (IM=IMM,JM=JMM,LM=LMM,LP1=LM+1)
  PARAMETER                                                             &
      (h10e5=100000.e0                                                  &
      , epsq=2.e-12                                                     &
      , g=9.8E0,cp=1004.6E0,capa=0.28589641E0,rog=287.04/g)
!----------------------------------------------------------------------
  PARAMETER                                                             &
      (itb=076,jtb=134,itbq=152,jtbq=440)
  PARAMETER (dpbnd=70.e2,                                               &
      elivw=2.72E6,elocp=elivw/cp)
!
!  DECLARE VARIABLES.
!
  INTEGER :: ieql(im,jm)
  INTEGER :: lcl(im,jm)
!
  REAL :: p1d(im,jm),t1d(im,jm),q1d(im,jm)
  REAL :: cape(im,jm),cins(im,jm),li(im,jm)
  REAL :: plcl(im,jm),peql(im,jm)
  REAL :: tpar(im,jm,lm)

  DIMENSION                                                             &
      p  (im,jm,lm),pint  (im,jm,lm+1)                                  &
      ,t(im,jm,lm),q(im,jm,lm)                                          &
      ,lmh  (im,jm)
  DIMENSION                                                             &
      qs0 (jtb),sqs (jtb),the0  (itb),sthe  (itb)                       &
      ,                    the0q(itbq),stheq(itbq)                      &
      ,ptbl  (itb,jtb),ttbl  (jtb,itb),ttblq(jtbq,itbq)
!
  DIMENSION                                                             &
      iptb  (im,jm),ithtb (im,jm)                                       &
      ,pp    (im,jm)                                                    &
      ,qq    (im,jm)                                                    &
      ,psp   (im,jm),thesp (im,jm)                                      &
      ,ilres (im*jm),jlres (im*jm),ihres (im*jm),jhres (im*jm)          &
      ,ape   (im,jm,lm)
!
!
  DATA itabl/0/,thl/210./,plq/70000./,ptt/2500./
!
  SAVE ptbl,ttbl,rdq,rdth,rdp,rdthe,pl,thl,qs0,sqs,sthe,the0,           &
      ttblq,rdpq,rdtheq,plq,stheq,the0q
!-----------------------------------------------------------------------
!
!   SET UP LOOKUP TABLE
!
  IF (itabl == 0) THEN
    CALL table1(ptbl,ttbl,ptt,                                          &
        rdq,rdth,rdp,rdthe,pl,thl,qs0,sqs,sthe,the0)
    CALL tableq(ttblq,rdpq,rdtheq,plq,thl,stheq,the0q)
    itabl=1
  END IF
  imjm=im*jm
!-----------------------------------------------------------------------
!
!
!**************************************************************
!  START CALCAPE HERE.

!
!
!  COMPUTE CAPE/CINS
!
!     WHICH IS: THE SUM FROM THE LCL TO THE EQ LEVEL OF
!          G * (LN(THETAP) - LN(THETAA)) * DZ
!
!          (POSITIVE AREA FOR CAPE, NEGATIVE FOR CINS)
!
!     WHERE:
!          THETAP IS THE PARCEL THETA
!          THETAA IS THE AMBIENT THETA
!          DZ IS THE THICKNESS OF THE LAYER
!
!      USING LCL AS LEVEL DIRECTLY BELOW SATURATION POINT
!      AND EQ LEVEL IS THE HIGHEST POSITIVELY BUOYANT LEVEL.
!
!      IEQL = EQ LEVEL
!
!  INITIALIZE CAPE AND CINS ARRAYS
!
  DO j = 1,jm
    DO i = 1,im
      cape(i,j) = 0.0
      cins(i,j) = 0.0
      thesp(i,j)= 0.0
      ieql(i,j) = lm+1
    END DO
  END DO
  DO l=1,lm
    DO j = 1,jm
      DO i = 1,im
        tpar(i,j,l)= 0.0
        IF (l <= lmh(i,j)) THEN
          apests=p(i,j,l)
          ape(i,j,l)=(h10e5/apests)**capa
          IF(        END IF
      END DO
    END DO
  END DO
!
!  NOTE THAT FOR TYPE 1 CAPE/CINS ARRAYS P1D, T1D, Q1D
!  ARE DUMMY ARRAYS.
!
!-------FOR ITYPE=1-----------------------------------------------------
!---------FIND MAXIMUM THETA E LAYER IN LOWEST DPBND ABOVE GROUND-------
!-------FOR ITYPE=2-----------------------------------------------------
!---------FIND THETA E LAYER OF GIVEN T1D, Q1D, P1D---------------------
!--------------TRIAL MAXIMUM BUOYANCY LEVEL VARIABLES-------------------
  DO kb=1,lm
    IF (itype == 1.OR.kb == 1) THEN
      DO j=1,jm
        DO i=1,im
          lmhk   =lmh(i,j)
          psfck  =p(i,j,lmhk)
          pkl = p(i,j,kb)
          IF (itype == 2.OR.(pkl >= psfck-dpbnd.AND.pkl <= psfck)) THEN
            IF (itype == 1) THEN
              tbtk   =t(i,j,kb)
              qbtk   =q(i,j,kb)
              apebtk =ape(i,j,kb)
!DHOU added the following statement to calculate the temperature of
! the parcil if it diabatically goes to 500hpa. This is stored in li.
              li(i,j)=tbtk*(50000.0/p(i,j,kb))**capa
            ELSE
              pkl    =p1d(i,j)
              tbtk   =t1d(i,j)
              qbtk   =q1d(i,j)
              apebtk =(h10e5/pkl)**capa
            END IF
!--------------SCALING POTENTIAL TEMPERATURE & TABLE INDEX--------------
            tthbtk =tbtk*apebtk
            tthk   =(tthbtk-thl)*rdth
            qq  (i,j)=tthk-AINT(tthk)
            ittbk  =INT(tthk)+1
!--------------KEEPING INDICES WITHIN THE TABLE-------------------------
            IF(ittbk < 1)    THEN
              ittbk  =1
              qq  (i,j)=0.0
            END IF
            IF(ittbk >= jtb)  THEN
              ittbk  =jtb-1
              qq  (i,j)=0.0
            END IF
!--------------BASE AND SCALING FACTOR FOR SPEC. HUMIDITY---------------
            bqs00k=qs0(ittbk)
            sqs00k=sqs(ittbk)
            bqs10k=qs0(ittbk+1)
            sqs10k=sqs(ittbk+1)
!--------------SCALING SPEC. HUMIDITY & TABLE INDEX---------------------
            bqk    =(bqs10k-bqs00k)*qq(i,j)+bqs00k
            sqk    =(sqs10k-sqs00k)*qq(i,j)+sqs00k
            tqk    =(qbtk-bqk)/sqk*rdq
            pp  (i,j)=tqk-AINT(tqk)
            iq     =INT(tqk)+1
!--------------KEEPING INDICES WITHIN THE TABLE-------------------------
            IF(iq < 1)    THEN
              iq     =1
              pp  (i,j)=0.0
            END IF

            IF(iq >= itb)  THEN
              iq     =itb-1
              pp  (i,j)=0.0
            END IF
!--------------SATURATION PRESSURE AT FOUR SURROUNDING TABLE PTS.-------
            it=ittbk
            p00k=ptbl(iq  ,it  )
            p10k=ptbl(iq+1,it  )
            p01k=ptbl(iq  ,it+1)
            p11k=ptbl(iq+1,it+1)
!--------------SATURATION POINT VARIABLES AT THE BOTTOM-----------------
            tpspk=p00k+(p10k-p00k)*pp(i,j)+(p01k-p00k)*qq(i,j)          &
                +(p00k-p10k-p01k+p11k)*pp(i,j)*qq(i,j)
            apespk=(h10e5/tpspk)**capa
            tthesk=tthbtk*EXP(elocp*qbtk*apespk/tthbtk)
!--------------CHECK FOR MAXIMUM THETA E--------------------------------
            IF(tthesk > thesp(i,j))    THEN
              psp  (i,j)=tpspk
              thesp(i,j)=tthesk
            END IF
          END IF
        END DO
      END DO
    END IF
  END DO
!-----CHOOSE LAYER DIRECTLY BELOW PSP AS LCL AND------------------------
!-----ENSURE THAT THE LCL IS ABOVE GROUND.------------------------------
!-------(IN SOME RARE CASES FOR ITYPE=2, IT IS NOT)---------------------
  DO l=1,lm
    DO j=1,jm
      DO i=1,im
        IF (l < lmh(i,j)) THEN
          pkl = p(i,j,l)
          IF (pkl < psp(i,j)) THEN
            lcl(i,j)=l+1
            plcl(i,j)=p(i,j,l+1)
          END IF
        END IF
      END DO
    END DO
  END DO
!-----------------------------------------------------------------------
!---------FIND TEMP OF PARCEL LIFTED ALONG MOIST ADIABAT (TPAR)---------
!-----------------------------------------------------------------------
  DO ivi=1,lm
    l=lp1-ivi
!--------------SCALING PRESSURE & TT TABLE INDEX------------------------
    knuml=0
    knumh=0
    DO j=1,jm
      DO i=1,im
        pkl = p(i,j,l)
        IF(l <= lcl(i,j)) THEN
          IF(pkl < plq)THEN
            knuml=knuml+1
            ilres(knuml)=i
            jlres(knuml)=j
          ELSE
            knumh=knumh+1
            ihres(knumh)=i
            jhres(knumh)=j
          END IF
        END IF
      END DO
    END DO
!***
!***  COMPUTE PARCEL TEMPERATURE ALONG MOIST ADIABAT FOR PRESSURE<PLQ
!**
    IF(knuml > 0)THEN
      CALL ttblex(tpar(1,1,l),ttbl,im,jm,imjm,itb,jtb                   &
          ,             knuml,ilres,jlres                               &
          ,             p(1,1,l),pl,qq,pp,rdp,the0,sthe                 &
          ,             rdthe,thesp,iptb,ithtb)
    END IF

!***
!***  COMPUTE PARCEL TEMPERATURE ALONG MOIST ADIABAT FOR PRESSURE>PLQ
!**
    IF(knumh > 0)THEN
      CALL ttblex(tpar(1,1,l),ttblq,im,jm,imjm,itbq,jtbq                &
          ,             knumh,ihres,jhres                               &
          ,             p(1,1,l),plq,qq,pp,rdpq,the0q,stheq             &
          ,             rdtheq,thesp,iptb,ithtb)
    END IF

!------------SEARCH FOR EQ LEVEL----------------------------------------
    DO n=1,knumh
      i=ihres(n)
      j=jhres(n)
      IF(tpar(i,j,l) > t(i,j,l)) THEN
        ieql(i,j)=l
        peql(i,j)=p(i,j,l)
      END IF
    END DO
    DO n=1,knuml
      i=ilres(n)
      j=jlres(n)
      IF(tpar(i,j,l) > t(i,j,l)) THEN
        ieql(i,j)=l

        peql(i,j)=p(i,j,l)
      END IF
    END DO
!-----------------------------------------------------------------------
  END DO
!-----------------------------------------------------------------------
!  DHOU added this section to calculate LI (Lifted Index=T500-Tp500)
  DO j=1,jm
    DO i=1,im
      DO k=1,lm
        IF (p(i,j,k+1) >= 50000.0.AND.p(i,j,k) < 50000.0) THEN
          t500=t(i,j,k+1)-(ALOG(              *(t(i,j,k+1)-t(i,j,k))/(ALOG(          IF (psp(i,j) < 50000.0) THEN
            tp500=li(i,j)
          ELSE
            tp500=tpar(i,j,k+1)-(ALOG(                *(tpar(i,j,k+1)-tpar(i,j,k))/(ALOG(          END IF
          li(i,j)=t500-tp500
          CYCLE
        END IF
      END DO
    END DO
  END DO
!------------COMPUTE CAPE AND CINS--------------------------------------
  DO j=1,jm
    DO i=1,im
      lclk=lcl(i,j)
      ieqk=ieql(i,j)
      DO l=ieqk,lclk
        presk=p(i,j,l)
        dp=pint(i,j,l+1)-pint(i,j,l)
        dzkl=t(i,j,l)*(q(i,j,l)*0.608+1.0)*rog*dp/presk
        thetap=tpar(i,j,l)*(h10e5/presk)**capa
        thetaa=t(i,j,l)*(h10e5/presk)**capa
        IF (thetap < thetaa)                                            &
            cins(i,j)=cins(i,j)+g*(ALOG(thetap)-ALOG(thetaa))*dzkl
        IF (thetap > thetaa)                                            &
            cape(i,j)=cape(i,j)+g*(ALOG(thetap)-ALOG(thetaa))*dzkl
      END DO
    END DO
  END DO
!
!  ENFORCE LOWER LIMIT OF 0.0 ON CAPE AND UPPER
!  LIMIT OF 0.0 ON CINS.
!
  DO j = 1,jm
    DO i = 1,im
      cape(i,j) = AMAX1(0.0,cape(i,j))
      cins(i,j) = AMIN1(cins(i,j),0.0)
    END DO
  END DO
!
!  END OF ROUTINE.
!
  RETURN
END SUBROUTINE calcape


SUBROUTINE table1(ptbl,ttbl,pt                                          & 1,2
           ,                rdq,rdth,rdp,rdthe,pl,thl,qs0,sqs,sthe,the0)
!  ******************************************************************
!  *                                                                *
!  *    GENERATE VALUES FOR LOOK-UP TABLES USED IN CONVECTION       *
!  *                                                                *
!  ******************************************************************
  PARAMETER                                                             &
      (itb=076,jtb=134)
  PARAMETER                                                             &
      (thh=365.,ph=105000.                                              &
      , pq0=379.90516,a1=610.78,a2=17.2693882,a3=273.16,a4=35.86        &
      , r=287.04,cp=1004.6,eliwv=2.683E6,eps=1.e-9)
  DIMENSION                                                             &
      ptbl  (itb,jtb),ttbl  (jtb,itb),qsold (jtb),pold  (jtb)           &
      , qs0   (jtb),sqs   (jtb),qsnew (jtb)                             &
      , y2p   (jtb),app   (jtb),aqp   (jtb),pnew  (jtb)                 &
      , told  (jtb),theold(jtb),the0  (itb),sthe  (itb)                 &
      , y2t   (jtb),thenew(jtb),apt   (jtb),aqt   (jtb),tnew  (jtb)
!--------------COARSE LOOK-UP TABLE FOR SATURATION POINT----------------
  kthm=jtb
  kpm=itb
  kthm1=kthm-1
  kpm1=kpm-1
!
  pl=pt
!
  dth=(thh-thl)/REAL(kthm-1)
  dp =(ph -pl )/REAL(kpm -1)
!
  rdth=1./dth
  rdp=1./dp
  rdq=kpm-1
!
  th=thl-dth
!-----------------------------------------------------------------------
  DO kth=1,kthm
    th=th+dth
    p=pl-dp
    DO kp=1,kpm
      p=p+dp
      ape=(100000./p)**(r/cp)
      qsold(kp)=pq0/p*EXP(a2*(th-a3*ape)/(th-a4*ape))
      pold(kp)=p
    END DO
!
    qs0k=qsold(1)
    sqsk=qsold(kpm)-qsold(1)
    qsold(1  )=0.
    qsold(kpm)=1.
!
    DO kp=2,kpm1
      qsold(kp)=(qsold(kp)-qs0k)/sqsk
!
      IF((qsold(kp)-qsold(kp-1)) < eps) qsold(kp)=qsold(kp-1)+eps
!
    END DO
!
    qs0(kth)=qs0k
    sqs(kth)=sqsk
!-----------------------------------------------------------------------
    qsnew(1  )=0.

    qsnew(kpm)=1.
    dqs=1./REAL(kpm-1)
!
    DO kp=2,kpm1
      qsnew(kp)=qsnew(kp-1)+dqs
    END DO
!
    y2p(1   )=0.
    y2p(kpm )=0.
!
    CALL spline(jtb,kpm,qsold,pold,y2p,kpm,qsnew,pnew,app,aqp)
!
    DO kp=1,kpm
      ptbl(kp,kth)=pnew(kp)
    END DO
!-----------------------------------------------------------------------
  END DO
!--------------COARSE LOOK-UP TABLE FOR T(P) FROM CONSTANT THE----------
  p=pl-dp
  DO kp=1,kpm
    p=p+dp
    th=thl-dth
    DO kth=1,kthm
      th=th+dth
      ape=(100000./p)**(r/cp)
      qs=pq0/p*EXP(a2*(th-a3*ape)/(th-a4*ape))
      told(kth)=th/ape
      theold(kth)=th*EXP(eliwv*qs/(cp*told(kth)))
    END DO
!
    the0k=theold(1)
    sthek=theold(kthm)-theold(1)
    theold(1   )=0.
    theold(kthm)=1.
!
    DO kth=2,kthm1
      theold(kth)=(theold(kth)-the0k)/sthek
!
      IF((theold(kth)-theold(kth-1)) < eps) theold(kth)=theold(kth-1)  +  eps
!
    END DO
!
    the0(kp)=the0k
    sthe(kp)=sthek
!-----------------------------------------------------------------------
    thenew(1  )=0.
    thenew(kthm)=1.
    dthe=1./REAL(kthm-1)
    rdthe=1./dthe
!
    DO kth=2,kthm1
      thenew(kth)=thenew(kth-1)+dthe
    END DO
!
    y2t(1   )=0.
    y2t(kthm)=0.
!
    CALL spline(jtb,kthm,theold,told,y2t,kthm,thenew,tnew,apt,aqt)
!
    DO kth=1,kthm
      ttbl(kth,kp)=tnew(kth)
    END DO
!-----------------------------------------------------------------------
  END DO
!
  RETURN
END SUBROUTINE table1
!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&


SUBROUTINE tableq(ttblq                                                 & 1,1
           ,                 rdp,rdthe,pl,thl,sthe,the0)
!  ******************************************************************
!  *                                                                *
!  *        GENERATE VALUES FOR FINER LOOK-UP TABLES USED           *
!  *                       IN CONVECTION                            *
!  *                                                                *
!  ******************************************************************
  PARAMETER                                                             &
      (itb=152,jtb=440)
  PARAMETER                                                             &
      (thh=325.,ph=105000.                                              &
      , pq0=379.90516,a1=610.78,a2=17.2693882,a3=273.16,a4=35.86        &
      , r=287.04,cp=1004.6,eliwv=2.683E6,eps=1.e-9)
  DIMENSION                                                             &
      ttblq (jtb,itb)                                                   &
      , told  (jtb),theold(jtb),the0  (itb),sthe  (itb)                 &
      , y2t   (jtb),thenew(jtb),apt   (jtb),aqt   (jtb),tnew  (jtb)
!--------------COARSE LOOK-UP TABLE FOR SATURATION POINT----------------
  kthm=jtb
  kpm=itb
  kthm1=kthm-1
  kpm1=kpm-1
!
  dth=(thh-thl)/REAL(kthm-1)
  dp =(ph -pl )/REAL(kpm -1)
!
  rdp=1./dp
  th=thl-dth
!--------------COARSE LOOK-UP TABLE FOR T(P) FROM CONSTANT THE----------
  p=pl-dp
  DO kp=1,kpm
    p=p+dp
    th=thl-dth
    DO kth=1,kthm
      th=th+dth
      ape=(100000./p)**(r/cp)
      qs=pq0/p*EXP(a2*(th-a3*ape)/(th-a4*ape))
      told(kth)=th/ape

      theold(kth)=th*EXP(eliwv*qs/(cp*told(kth)))
    END DO
!
    the0k=theold(1)
    sthek=theold(kthm)-theold(1)
    theold(1   )=0.
    theold(kthm)=1.
!
    DO kth=2,kthm1
      theold(kth)=(theold(kth)-the0k)/sthek
!
      IF((theold(kth)-theold(kth-1)) < eps) theold(kth)=theold(kth-1)  +  eps
!
    END DO
!
    the0(kp)=the0k
    sthe(kp)=sthek
!-----------------------------------------------------------------------
    thenew(1  )=0.
    thenew(kthm)=1.
    dthe=1./REAL(kthm-1)
    rdthe=1./dthe
!
    DO kth=2,kthm1
      thenew(kth)=thenew(kth-1)+dthe
    END DO
!
    y2t(1   )=0.
    y2t(kthm)=0.
!
    CALL spline(jtb,kthm,theold,told,y2t,kthm,thenew,tnew,apt,aqt)
!
    DO kth=1,kthm
      ttblq(kth,kp)=tnew(kth)
    END DO
!-----------------------------------------------------------------------
  END DO
!
  RETURN
END SUBROUTINE tableq


SUBROUTINE ttblex(tref,ttbl,im,jm,imxjm,itb,jtb                         & 4
           ,                 knum,iarr,jarr                             &
           ,                 pijl,pl,qq,pp,rdp,the0                     &
           ,                 sthe,rdthe,thesp,iptb,ithtb)
!  ******************************************************************
!  *                                                                *
!  *         EXTRACT TEMPERATURE OF THE MOIST ADIABAT FROM          *
!  *                    THE APPROPRIATE TTBL                        *
!  *                                                                *
!  ******************************************************************
!----------------------------------------------------------------------
!  PARAMETER(IM=1,JM=1,LM=50,IMXJM=IM*JM)
  DIMENSION                                                             &
      tref(im,jm),ttbl(jtb,itb),iarr(imxjm),jarr(imxjm)                 &
      ,pijl(im,jm),qq(im,jm),pp(im,jm),the0(itb)                        &
      ,sthe(itb),thesp(im,jm),iptb(im,jm),ithtb(im,jm)
!-----------------------------------------------------------------------
  DO kk=1,knum
!--------------SCALING PRESSURE & TT TABLE INDEX------------------------
    i=iarr(kk)
    j=jarr(kk)
    pk=pijl(i,j)
    tpk=(pk-pl)*rdp
    qq(i,j)=tpk-AINT(tpk)
    iptb(i,j)=INT(tpk)+1
!--------------KEEPING INDICES WITHIN THE TABLE-------------------------
    IF(iptb(i,j) < 1)THEN
      iptb(i,j)=1
      qq(i,j)=0.
    END IF
    IF(iptb(i,j) >= itb)THEN
      iptb(i,j)=itb-1
      qq(i,j)=0.
    END IF
!--------------BASE AND SCALING FACTOR FOR THE--------------------------
    iptbk=iptb(i,j)
    bthe00k=the0(iptbk)
    sthe00k=sthe(iptbk)
    bthe10k=the0(iptbk+1)
    sthe10k=sthe(iptbk+1)
!--------------SCALING THE & TT TABLE INDEX-----------------------------
    bthk=(bthe10k-bthe00k)*qq(i,j)+bthe00k
    sthk=(sthe10k-sthe00k)*qq(i,j)+sthe00k
    tthk=(thesp(i,j)-bthk)/sthk*rdthe
    pp(i,j)=tthk-AINT(tthk)
    ithtb(i,j)=INT(tthk)+1
!--------------KEEPING INDICES WITHIN THE TABLE-------------------------
    IF(ithtb(i,j) < 1)THEN
      ithtb(i,j)=1
      pp(i,j)=0.
    END IF
    IF(ithtb(i,j) >= jtb)THEN
      ithtb(i,j)=jtb-1

      pp(i,j)=0.
    END IF
!--------------TEMPERATURE AT FOUR SURROUNDING TT TABLE PTS.------------
    ith=ithtb(i,j)
    ip=iptb(i,j)
    t00k=ttbl(ith  ,ip  )
    t10k=ttbl(ith+1,ip  )
    t01k=ttbl(ith  ,ip+1)
    t11k=ttbl(ith+1,ip+1)
!--------------PARCEL TEMPERATURE-------------------------------------
    tref(i,j)=(t00k+(t10k-t00k)*pp(i,j)+(t01k-t00k)*qq(i,j)             &
             +(t00k-t10k-t01k+t11k)*pp(i,j)*qq(i,j))
  END DO
!
  RETURN
END SUBROUTINE ttblex
!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&


SUBROUTINE spline(jtb,nold,xold,yold,y2,nnew,xnew,ynew,p,q) 6
!  ******************************************************************
!  *                                                                *
!  *  THIS IS A ONE-DIMENSIONAL CUBIC SPLINE FITTING ROUTINE        *
!  *  PROGRAMED FOR A SMALL SCALAR MACHINE.                         *
!  *                                                                *
!  *  PROGRAMER Z. JANJIC, YUGOSLAV FED. HYDROMET. INST., BEOGRAD  *
!  *                                                                *
!  *                                                                *
!  *                                                                *
!  *  NOLD - NUMBER OF GIVEN VALUES OF THE FUNCTION.  MUST BE GE 3. *
!  *  XOLD - LOCATIONS OF THE POINTS AT WHICH THE VALUES OF THE     *
!  *         FUNCTION ARE GIVEN.  MUST BE IN ASCENDING ORDER.       *
!  *  YOLD - THE GIVEN VALUES OF THE FUNCTION AT THE POINTS XOLD.   *
!  *  Y2   - THE SECOND DERIVATIVES AT THE POINTS XOLD.  IF NATURAL *
!  *         SPLINE IS FITTED Y2(1)=0. AND Y2(NOLD)=0. MUST BE      *
!  *         SPECIFIED.                                             *
!  *  NNEW - NUMBER OF VALUES OF THE FUNCTION TO BE CALCULATED.     *
!  *  XNEW - LOCATIONS OF THE POINTS AT WHICH THE VALUES OF THE     *
!  *         FUNCTION ARE CALCULATED.  XNEW(K) MUST BE GE XOLD(1)   *
!  *         AND LE XOLD(NOLD).                                     *
!  *  YNEW - THE VALUES OF THE FUNCTION TO BE CALCULATED.           *
!  *  P, Q - AUXILIARY VECTORS OF THE LENGTH NOLD-2.                *
!  *                                                                *
!  ******************************************************************
  DIMENSION                                                             &
      xold(jtb),yold(jtb),y2(jtb),p(jtb),q(jtb)                         &
      ,xnew(jtb),ynew(jtb)
!-----------------------------------------------------------------------
  noldm1=nold-1
!
  dxl=xold(2)-xold(1)
  dxr=xold(3)-xold(2)
  dydxl=(yold(2)-yold(1))/dxl
  dydxr=(yold(3)-yold(2))/dxr
  rtdxc=.5/(dxl+dxr)
!
  p(1)= rtdxc*(6.*(dydxr-dydxl)-dxl*y2(1))
  q(1)=-rtdxc*dxr
!
  IF(nold == 3) GO TO 700
!-----------------------------------------------------------------------
  k=3
!
  100   dxl=dxr
  dydxl=dydxr
  dxr=xold(k+1)-xold(k)
  dydxr=(yold(k+1)-yold(k))/dxr
  dxc=dxl+dxr

  den=1./(dxl*q(k-2)+dxc+dxc)
!
  p(k-1)= den*(6.*(dydxr-dydxl)-dxl*p(k-2))
  q(k-1)=-den*dxr
!
  k=k+1
  IF(k < nold) GO TO 100
!-----------------------------------------------------------------------
  700   k=noldm1
!
  200   y2(k)=p(k-1)+q(k-1)*y2(k+1)
!
  k=k-1
  IF(k > 1) GO TO 200
!-----------------------------------------------------------------------
  k1=1
!
  300   xk=xnew(k1)
!
  DO k2=2,nold
    IF(xold(k2) <= xk) CYCLE
    kold=k2-1
    GO TO 450
  END DO
  ynew(k1)=yold(nold)
  GO TO 600
!
  450   IF(k1 == 1)   GO TO 500
  IF(k == kold) GO TO 550
!
  500   k=kold
!
  y2k=y2(k)
  y2kp1=y2(k+1)
  dx=xold(k+1)-xold(k)
  rdx=1./dx
!
!VVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVV
!  WRITE(6,5000) K,Y2K,Y2KP1,DX,RDX,YOLD(K),YOLD(K+1)
!5000 FORMAT(' K=',I4,' Y2K=',E12.4,' Y2KP1=',E12.4,' DX=',E12.4,' RDX='
!    2,E12.4,' YOK=',E12.4,' YOP1=',E12.4)
!AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
  ak=.1666667*rdx*(y2kp1-y2k)
  bk=.5*y2k
  ck=rdx*(yold(k+1)-yold(k))-.1666667*dx*(y2kp1+y2k+y2k)
!
  550   x=xk-xold(k)
  xsq=x*x
!
  ynew(k1)=ak*xsq*x+bk*xsq+ck*x+yold(k)
!
  600 CONTINUE

  k1=k1+1
  IF(k1 <= nnew) GO TO 300
!-----------------------------------------------------------------------
  RETURN
END SUBROUTINE spline


SUBROUTINE newmap(nqi,nqj,latn,lonn,dx,                                 & 1,4
           qswcorn,qswcorw,trulat,trulon,mapproj)
  IMPLICIT NONE
  INTEGER :: nqi,nqj,i,j
  REAL :: latn(nqi,nqj),lonn(nqi,nqj)
  REAL :: qswcorn,qswcorw,dx
  REAL :: trulat(2),trulon
  INTEGER :: mapproj
!
  REAL :: x,y,x1,y1
  REAL :: lat,lon
!
  CALL setmapr(mapproj,1.0,trulat,trulon)
  lat=qswcorn
  lon=qswcorw
  CALL lltoxy(1,1,lat,lon,x,y)
  CALL setorig(1,x,y)
  DO j=1,nqj
    DO i=1,nqi
      x1=(i-1)*dx*1000.0
      y1=(j-1)*dx*1000.0
      CALL xytoll(1,1,x1,y1,lat,lon)
      latn(i,j)=lat
      lonn(i,j)=lon
    END DO
  END DO
  RETURN
END SUBROUTINE newmap


SUBROUTINE d2intrpa(xold,yold,nxold,nyold,                              & 1
           xnew,ynew,nxnew,nynew,intrpl,                                &
           aold1,anew1,                                                 &
           aold2,anew2,                                                 &
           aold3,anew3,                                                 &
           aold4,anew4,                                                 &
           aold5,anew5,                                                 &
           aold6,anew6,                                                 &
           aold7,anew7,                                                 &
           aold8,anew8,                                                 &
           aold9,anew9,                                                 &
           aold10,anew10,                                               &
           aold11,anew11,                                               &
           aold12,anew12,                                               &
           aold13,anew13,                                               &
           aold14,anew14,                                               &
           aold15,anew15,                                               &
           aold16,anew16,                                               &
           aold17,anew17,                                               &
           aold18,anew18,                                               &
           aold19,anew19,                                               &
           aold20,anew20,                                               &
           aold21,anew21,                                               &
           aold22,anew22,                                               &
           aold23,anew23,                                               &
           aold24,anew24,                                               &
           aold25,anew25,                                               &
           aold26,anew26,                                               &
           aold27,anew27,                                               &
           aold28,anew28,                                               &
           aold29,anew29)
  IMPLICIT NONE
  INTEGER :: nxold,nyold, nxnew,nynew,intrpl
  REAL :: xold(nxold,nyold),xnew(nxnew,nynew)
  REAL :: yold(nxold,nyold),ynew(nxnew,nynew)
  REAL :: aold1(nxold,nyold),anew1(nxnew,nynew)
  REAL :: aold2(nxold,nyold),anew2(nxnew,nynew)
  REAL :: aold3(nxold,nyold),anew3(nxnew,nynew)
  REAL :: aold4(nxold,nyold),anew4(nxnew,nynew)
  REAL :: aold5(nxold,nyold),anew5(nxnew,nynew)
  REAL :: aold6(nxold,nyold),anew6(nxnew,nynew)
  REAL :: aold7(nxold,nyold),anew7(nxnew,nynew)
  REAL :: aold8(nxold,nyold),anew8(nxnew,nynew)
  REAL :: aold9(nxold,nyold),anew9(nxnew,nynew)
  REAL :: aold10(nxold,nyold),anew10(nxnew,nynew)
  REAL :: aold11(nxold,nyold),anew11(nxnew,nynew)
  REAL :: aold12(nxold,nyold),anew12(nxnew,nynew)
  REAL :: aold13(nxold,nyold),anew13(nxnew,nynew)
  REAL :: aold14(nxold,nyold),anew14(nxnew,nynew)
  REAL :: aold15(nxold,nyold),anew15(nxnew,nynew)
  REAL :: aold16(nxold,nyold),anew16(nxnew,nynew)
  REAL :: aold17(nxold,nyold),anew17(nxnew,nynew)
  REAL :: aold18(nxold,nyold),anew18(nxnew,nynew)
  REAL :: aold19(nxold,nyold),anew19(nxnew,nynew)
  REAL :: aold20(nxold,nyold),anew20(nxnew,nynew)
  REAL :: aold21(nxold,nyold),anew21(nxnew,nynew)
  REAL :: aold22(nxold,nyold),anew22(nxnew,nynew)
  REAL :: aold23(nxold,nyold),anew23(nxnew,nynew)
  REAL :: aold24(nxold,nyold),anew24(nxnew,nynew)
  REAL :: aold25(nxold,nyold),anew25(nxnew,nynew)
  REAL :: aold26(nxold,nyold),anew26(nxnew,nynew)
  REAL :: aold27(nxold,nyold),anew27(nxnew,nynew)
  REAL :: aold28(nxold,nyold),anew28(nxnew,nynew)
  REAL :: aold29(nxold,nyold),anew29(nxnew,nynew)
  INTEGER :: i,j,i1,i2,j1,j2
  INTEGER :: io,jo
  REAL :: xs,ys,s1,s2,s3,s4,sgrid

  PRINT *,'starting D2INTRPA'

  PRINT *,nxold,nyold
  PRINT *,nxnew,nynew
  PRINT *,xold(1,1),yold(1,1),xnew(1,1),ynew(1,1)
  PRINT *,xold(1,nyold),yold(1,nyold),xnew(1,nynew),ynew(1,nynew)
  PRINT *,xold(nxold,nyold),yold(nxold,nyold),                          &
          xnew(nxnew,nynew),ynew(nxnew,nynew)
  PRINT *,xold(nxold,1),yold(nxold,1),xnew(nxnew,1),ynew(nxnew,1)

  DO j=1,nynew
    DO i=1,nxnew

! Direct assignment in the case of intrpl=0 (No intepolation)
      IF (intrpl == 0) THEN
        io=i+1
        jo=j+1
        anew1(i,j)=aold1(io,jo)
        anew2(i,j)=aold2(io,jo)
        anew3(i,j)=aold3(io,jo)
        anew4(i,j)=aold4(io,jo)
        anew5(i,j)=aold5(io,jo)
        anew6(i,j)=aold6(io,jo)
        anew7(i,j)=aold7(io,jo)
        anew8(i,j)=aold8(io,jo)
        anew9(i,j)=aold9(io,jo)
        anew10(i,j)=aold10(io,jo)
        anew11(i,j)=aold11(io,jo)
        anew12(i,j)=aold12(io,jo)
        anew13(i,j)=aold13(io,jo)
        anew14(i,j)=aold14(io,jo)
        anew15(i,j)=aold15(io,jo)
        anew16(i,j)=aold16(io,jo)
        anew17(i,j)=aold17(io,jo)
        anew18(i,j)=aold18(io,jo)
        anew19(i,j)=aold19(io,jo)
        anew20(i,j)=aold20(io,jo)
        anew21(i,j)=aold21(io,jo)
        anew22(i,j)=aold22(io,jo)
        anew23(i,j)=aold23(io,jo)
        anew24(i,j)=aold24(io,jo)
        anew25(i,j)=aold25(io,jo)
        anew26(i,j)=aold26(io,jo)
        anew27(i,j)=aold27(io,jo)
        anew28(i,j)=aold28(io,jo)
        anew29(i,j)=aold29(io,jo)
        CYCLE
      END IF

!  Interpolation, first finding the old points surroding the new point
      xs=xnew(i,j)
      ys=ynew(i,j)

      IF (xs < xold(1,1)) THEN
        i1=1
        i2=2
      ELSE IF (xs >= xold(nxold,1)) THEN
        i1=nxold-1
        i2=nxold
      ELSE
        DO io=1,nxold-1
          IF (xs >= xold(io,1).AND.xs < xold(io+1,1)) THEN
            i1=io
            i2=io+1
            EXIT
          END IF
        END DO
        205     CONTINUE
      END IF

      IF (ys < yold(1,1)) THEN
        j1=1
        j2=2
      ELSE IF (ys >= yold(1,nyold)) THEN
        j1=nyold-1
        j2=nyold
      ELSE
        DO jo=1,nyold-1
          IF (ys >= yold(1,jo).AND.ys < yold(1,jo+1)) THEN
            j1=jo
            j2=jo+1
            EXIT
          END IF
        END DO
        305     CONTINUE
      END IF

      DO jo=1,nyold
        DO io=1,nxold
! dirtect assignment for a new point, very close to an old point.
          IF (ABS(xs-xold(io,jo)) < 1.0E2.AND. ABS(ys-yold(io,jo)) < 1.0E2) THEN
            anew1(i,j)=aold1(io,jo)
            anew2(i,j)=aold2(io,jo)
            anew3(i,j)=aold3(io,jo)
            anew4(i,j)=aold4(io,jo)
            anew5(i,j)=aold5(io,jo)
            anew6(i,j)=aold6(io,jo)
            anew7(i,j)=aold7(io,jo)
            anew8(i,j)=aold8(io,jo)
            anew9(i,j)=aold9(io,jo)
            anew10(i,j)=aold10(io,jo)
            anew11(i,j)=aold11(io,jo)
            anew12(i,j)=aold12(io,jo)
            anew13(i,j)=aold13(io,jo)
            anew14(i,j)=aold14(io,jo)
            anew15(i,j)=aold15(io,jo)
            anew16(i,j)=aold16(io,jo)
            anew17(i,j)=aold17(io,jo)
            anew18(i,j)=aold18(io,jo)
            anew19(i,j)=aold19(io,jo)
            anew20(i,j)=aold20(io,jo)
            anew21(i,j)=aold21(io,jo)
            anew22(i,j)=aold22(io,jo)
            anew23(i,j)=aold23(io,jo)
            anew24(i,j)=aold24(io,jo)
            anew25(i,j)=aold25(io,jo)
            anew26(i,j)=aold26(io,jo)
            anew27(i,j)=aold27(io,jo)
            anew28(i,j)=aold28(io,jo)
            anew29(i,j)=aold29(io,jo)
            GO TO 405
          END IF
        END DO
      END DO

!   interpolation
      s1=SQRT((xs-xold(i1,j1))**2+(ys-yold(i1,j1))**2)
      s2=SQRT((xs-xold(i2,j1))**2+(ys-yold(i2,j1))**2)
      s3=SQRT((xs-xold(i2,j2))**2+(ys-yold(i2,j2))**2)
      s4=SQRT((xs-xold(i1,j2))**2+(ys-yold(i1,j2))**2)
      s1=1.0/s1
      s2=1.0/s2
      s3=1.0/s3
      s4=1.0/s4
      sgrid=1.0/(s1+s2+s3+s4)
      anew1(i,j)=(aold1(i1,j1)*s1+aold1(i2,j1)*s2                       &
                +aold1(i2,j2)*s3+aold1(i1,j2)*s4)*sgrid
      anew2(i,j)=(aold2(i1,j1)*s1+aold2(i2,j1)*s2                       &
                +aold2(i2,j2)*s3+aold2(i1,j2)*s4)*sgrid
      anew3(i,j)=(aold3(i1,j1)*s1+aold3(i2,j1)*s2                       &
                +aold3(i2,j2)*s3+aold3(i1,j2)*s4)*sgrid
      anew4(i,j)=(aold4(i1,j1)*s1+aold4(i2,j1)*s2                       &
                +aold4(i2,j2)*s3+aold4(i1,j2)*s4)*sgrid
      anew5(i,j)=(aold5(i1,j1)*s1+aold5(i2,j1)*s2                       &
                +aold5(i2,j2)*s3+aold5(i1,j2)*s4)*sgrid
      anew6(i,j)=(aold6(i1,j1)*s1+aold6(i2,j1)*s2                       &
                +aold6(i2,j2)*s3+aold6(i1,j2)*s4)*sgrid
      anew7(i,j)=(aold7(i1,j1)*s1+aold7(i2,j1)*s2                       &
                +aold7(i2,j2)*s3+aold7(i1,j2)*s4)*sgrid
      anew8(i,j)=(aold8(i1,j1)*s1+aold8(i2,j1)*s2                       &
                +aold8(i2,j2)*s3+aold8(i1,j2)*s4)*sgrid
      anew9(i,j)=(aold9(i1,j1)*s1+aold9(i2,j1)*s2                       &
                +aold9(i2,j2)*s3+aold9(i1,j2)*s4)*sgrid
      anew10(i,j)=(aold10(i1,j1)*s1+aold10(i2,j1)*s2                    &
                +aold10(i2,j2)*s3+aold10(i1,j2)*s4)*sgrid
      anew11(i,j)=(aold11(i1,j1)*s1+aold11(i2,j1)*s2                    &
                +aold11(i2,j2)*s3+aold11(i1,j2)*s4)*sgrid
      anew12(i,j)=(aold12(i1,j1)*s1+aold12(i2,j1)*s2                    &
                +aold12(i2,j2)*s3+aold12(i1,j2)*s4)*sgrid
      anew13(i,j)=(aold13(i1,j1)*s1+aold13(i2,j1)*s2                    &
                +aold13(i2,j2)*s3+aold13(i1,j2)*s4)*sgrid
      anew14(i,j)=(aold14(i1,j1)*s1+aold14(i2,j1)*s2                    &
                +aold14(i2,j2)*s3+aold14(i1,j2)*s4)*sgrid
      anew15(i,j)=(aold15(i1,j1)*s1+aold15(i2,j1)*s2                    &
                +aold15(i2,j2)*s3+aold15(i1,j2)*s4)*sgrid
      anew16(i,j)=(aold16(i1,j1)*s1+aold16(i2,j1)*s2                    &
                +aold16(i2,j2)*s3+aold16(i1,j2)*s4)*sgrid
      anew17(i,j)=(aold17(i1,j1)*s1+aold17(i2,j1)*s2                    &
                +aold17(i2,j2)*s3+aold17(i1,j2)*s4)*sgrid
      anew18(i,j)=(aold18(i1,j1)*s1+aold18(i2,j1)*s2                    &
                +aold18(i2,j2)*s3+aold18(i1,j2)*s4)*sgrid
      anew19(i,j)=(aold19(i1,j1)*s1+aold19(i2,j1)*s2                    &
                +aold19(i2,j2)*s3+aold19(i1,j2)*s4)*sgrid
      anew20(i,j)=(aold20(i1,j1)*s1+aold20(i2,j1)*s2                    &
                +aold20(i2,j2)*s3+aold20(i1,j2)*s4)*sgrid
      anew21(i,j)=(aold21(i1,j1)*s1+aold21(i2,j1)*s2                    &
                +aold21(i2,j2)*s3+aold21(i1,j2)*s4)*sgrid
      anew22(i,j)=(aold22(i1,j1)*s1+aold22(i2,j1)*s2                    &
                +aold22(i2,j2)*s3+aold22(i1,j2)*s4)*sgrid
      anew23(i,j)=(aold23(i1,j1)*s1+aold23(i2,j1)*s2                    &
                +aold23(i2,j2)*s3+aold23(i1,j2)*s4)*sgrid
      anew24(i,j)=(aold24(i1,j1)*s1+aold24(i2,j1)*s2                    &
                +aold24(i2,j2)*s3+aold24(i1,j2)*s4)*sgrid
      anew25(i,j)=(aold25(i1,j1)*s1+aold25(i2,j1)*s2                    &
                +aold25(i2,j2)*s3+aold25(i1,j2)*s4)*sgrid
      anew26(i,j)=(aold26(i1,j1)*s1+aold26(i2,j1)*s2                    &
                +aold26(i2,j2)*s3+aold26(i1,j2)*s4)*sgrid
      anew27(i,j)=(aold27(i1,j1)*s1+aold27(i2,j1)*s2                    &
                +aold27(i2,j2)*s3+aold27(i1,j2)*s4)*sgrid
      anew28(i,j)=(aold28(i1,j1)*s1+aold28(i2,j1)*s2                    &
                +aold28(i2,j2)*s3+aold28(i1,j2)*s4)*sgrid
      anew29(i,j)=(aold29(i1,j1)*s1+aold29(i2,j1)*s2                    &
                +aold29(i2,j2)*s3+aold29(i1,j2)*s4)*sgrid

!    if (i.eq.1.and.j.eq.1.or.i.eq.nxnew.and.j.eq.nynew) then
!      print *,i1,j1,i2,j2
!       print *,anew1(i,j),aold1(i1,j1),aold1(i2,j1),
!    :             aold1(i2,j2),aold1(i1,j2)
!       print *,anew2(i,j),aold2(i1,j1),aold2(i2,j1),
!    :             aold2(i2,j2),aold2(i1,j2)
!    endif
      405     CONTINUE
    END DO
  END DO
  PRINT *,'finishing D2INTRPA'

  RETURN
END SUBROUTINE d2intrpa
!
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE FILZERO                    ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE filzero( a, n) 4

!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!    Fill vector a with zeros.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: St. Paul, Dead Sea Scrolls
!
!  MODIFICATIONS:
!    6/09/92  Added full documentation (K. Brewster)
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!
  INTEGER :: n
  REAL :: a(n)
!
!-----------------------------------------------------------------------
!
!  Misc. local variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: i
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  DO i=1,n
    a(i)=0.0
  END DO

  RETURN
END SUBROUTINE filzero
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE IFILZERO                   ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE ifilzero( ia, n )

!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!    Fill vector a with zeros.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: St. Paul, Dead Sea Scrolls
!
!  MODIFICATIONS:
!    6/09/92  Added full documentation (K. Brewster)
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!
  INTEGER :: n
  INTEGER :: ia(n)
!
!-----------------------------------------------------------------------
!
!  Misc. local variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: i
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  DO i=1,n
    ia(i)=0
  END DO

  RETURN
END SUBROUTINE ifilzero
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE TEMPER                     ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE temper ( nx,ny,nz,theta, ppert, pbar, t ) 1

!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!    Using a version of Poisson's formula, calculate temperature.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Joe Bradley
!    12/05/91
!
!  MODIFICATIONS:
!    Modified by Ming Xue so that arrays are only defined at
!             one time level.
!    6/09/92  Added full documentation and phycst include file for
!             rddcp=Rd/Cp  (K. Brewster)
!
!-----------------------------------------------------------------------
!
!  INPUT:
!    nx       Number of grid points in the x-direction (east/west)
!    ny       Number of grid points in the y-direction (north/south)
!    nz       Number of grid points in the vertical
!
!    theta    Potential temperature (degrees Kelvin)
!    ppert    Perturbation pressure (Pascals)
!    pbar     Base state pressure (Pascals)
!
!  OUTPUT:
!
!    t        Temperature (degrees Kelvin)
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!
  INTEGER :: nx,ny,nz
!
  REAL :: theta(nx,ny,nz)      ! potential temperature (degrees Kelvin)
  REAL :: ppert(nx,ny,nz)      ! perturbation pressure (Pascals)
  REAL :: pbar (nx,ny,nz)      ! base state pressure (Pascals)
!
  REAL :: t    (nx,ny,nz)      ! temperature (degrees Kelvin)
!
!-----------------------------------------------------------------------
!
!  Include file
!
!-----------------------------------------------------------------------
!
  INCLUDE 'phycst.inc'
!
!-----------------------------------------------------------------------
!
!  Misc. local variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: i,j,k
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
!-----------------------------------------------------------------------
!
!  Calculate the temperature using Poisson's formula.
!
!-----------------------------------------------------------------------
!
  DO k=1,nz-1
    DO j=1,ny-1
      DO i=1,nx-1

        t(i,j,k) = theta(i,j,k) *                                       &
             (((ppert(i,j,k) + pbar(i,j,k)) / p0) ** rddcp)

      END DO
    END DO
  END DO

  RETURN
END SUBROUTINE temper


SUBROUTINE calender(year,mon,day,hour,ds,wkd) 1
! find the calender date for forecast initiated at hourZ,dayth of mon, Year
! with a lead time of ds seconds. the day of week is also found.
! effective only for date later than 1st Jan 1998, not exceeding 28 Feb. 2100
  IMPLICIT NONE
  INTEGER :: year,mon,day,hour,MIN,sec,ds,wkd
  INTEGER :: mon1,dday
  INTEGER :: y0,wkd0 !on Jan. 1st,the year of y0,the weekday index is wkd0
  INTEGER :: i,dindex
  y0=1998
  wkd0=5
! cross minute hour and day
  sec=0
  MIN=0
  sec=sec+ds
  MIN=MIN+(sec-MOD(sec,60))/60
  sec=MOD(sec,60)
  hour=hour+(MIN-MOD(MIN,60))/60
  MIN=MOD(MIN,60)
  dday=day+(hour-MOD(hour,24))/24
  hour=MOD(hour,24)
  ds=MIN
! cross the month
  IF  ( (mon == 1.OR.mon == 3.OR.mon == 5.OR.mon == 7                   &
        .OR.mon == 8.OR.mon == 10.OR.mon == 12).AND.dday > 31) THEN
    dday=dday-31
    mon1=mon+1
  ELSE IF ( (mon == 4.OR.mon == 6.OR.mon == 9.OR.mon == 11)             &
            .AND.dday > 30) THEN
    dday=dday-30
    mon1=mon+1
  ELSE IF (mon == 2) THEN
    IF (MOD(year,4) == 0.AND.dday > 29) THEN
      dday=dday-29
      mon1=mon+1
    ELSE IF (MOD(year,4) /= 0.AND.dday > 28) THEN
      dday=dday-28
      mon1=mon+1
    ELSE
      mon1=mon
    END IF
  ELSE
    mon1=mon
  END IF
  day=dday
  mon=mon1
!*  cross the year
  IF (mon > 12) THEN
    mon=mon-12
    year=year+1
  END IF
!*  day of week
  dindex=-1
  DO i=y0,year-1
    dindex=dindex+365
    IF (MOD(i,4) == 0) dindex=dindex+1
  END DO
  DO i=1,mon-1
    IF (i == 1.OR.i == 3.OR.i == 5.OR.i == 7                            &
          .OR.i == 8.OR.i == 10.OR.i == 12) THEN
      dindex=dindex+31
    ELSE IF (i == 4.OR.i == 6.OR.i == 9.OR.i == 11) THEN
      dindex=dindex+30
    ELSE IF (i == 2.AND.MOD(year,4) == 0) THEN
      dindex=dindex+29
    ELSE IF (i == 2.AND.MOD(year,4) /= 0) THEN
      dindex=dindex+28
    END IF
  END DO
  dindex=dindex+day
  wkd=wkd0+dindex
  PRINT *,wkd
  wkd=MOD(wkd,7)
  IF (wkd == 0) THEN
    wkd=7
  END IF
  RETURN
END SUBROUTINE calender