!
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE INITIAL                    ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!
!EMK BMJ

SUBROUTINE initial(mptr,nx,ny,nz,nxndg,nyndg,nzndg,nstyps,exbcbufsz,    & 3,9
           u,v,w,wcont,ptprt,pprt,qv,qc,qr,qi,qs,qh,tke,                &
           udteb,udtwb,udtnb,udtsb,vdteb,vdtwb,vdtnb,vdtsb,             &
           wdteb,wdtwb,wdtnb,wdtsb,pdteb,pdtwb,pdtnb,pdtsb,             &
           sdteb,sdtwb,sdtnb,sdtsb,                                     &
           ubar,vbar,ptbar,pbar,ptbari,pbari,                           &
           rhostr,rhostri,qvbar,ppi,csndsq,                             &
           x,y,z,zp,hterain, mapfct,                                    &
           j1,j2,j3,aj3x,aj3y,aj3z,j3inv,                               &
           trigs1,trigs2,ifax1,ifax2,                                   &
           wsave1,wsave2,vwork1,vwork2,                                 &
           sinlat, kmh,kmv,rprntl,                                      &
           soiltyp,stypfrct,vegtyp,lai,roufns,veg,                      &
           tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth,ptsfc,qvsfc,        &
           ptcumsrc,qcumsrc,w0avg,nca,kfraincv,                         &
           cldefi,xland,bmjraincv,                                      &
           raing,rainc,prcrate,exbcbuf,bcrlx,radfrc,radsw,              &
           rnflx,usflx,vsflx,ptsflx,qvsflx,                             &
           uincr,vincr,wincr,pincr,ptincr,qvincr,                       &
           qcincr,qrincr,qiincr,qsincr,qhincr,                          &
           temxy1,tem1,tem2,tem3,tem4,tem5,tem6,tem7,                   &
           tem8,tem9,tem10,tem11,tem12,tem13,                           &
           tem14,tem15,tem16,tem17,tem18,tem19,                         &
           tem20,tem21,tem22,tem23,tem24,tem25,tem26,                   &
           tem1_0,tem2_0,tem3_0)
!EMK BMJ

!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Initialize the model parameters and variables, including base state
!  variables, dependent variables and grid structure.
!
!  This is the main driver for model initializations.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Ming Xue
!  11/5/92.
!
!  MODIFICATION HISTORY:
!
!  5/02/92 (M. Xue)
!  Added full documentation.
!
!  5/03/92 (M. Xue)
!  Further documentation.
!
!  9/14/1992 (M. Xue)
!  Different surface drag coefficients defined for momentum, heat and
!  moisture.
!  Three options included for the Coriolis force calculations.
!
!  2/12/94 (Yuhe Liu)
!  Surface data and variables added for surface energy budget model.
!
!  6/10/94 (M. Xue &AS)
!  Added call to initpwr.
!
!  02/07/1995 (Yuhe Liu)
!  Added a new 2-D permanent array, veg(nx,ny), to the argument list
!
!  08/30/1995 (Yuhe Liu)
!  Moved the initialization of external boundary arrays from the
!  main program to this subroutine.
!
!  10/31/95   (D. Weber)
!  Added trigs1,trigs2,ifax1,ifax2 for use in the upper w-p
!  radiation condition.
!
!  1/22/96 (Donghai Wang & Yuhe Liu)
!  Added the map projection factor to ARPS governing equations.
!
!  07/22/97 (D. Weber)
!  Added wsave1,wsave2,vwork1,vwork2 for use in the even fft version
!  of the upper w-p radiation condition (fftopt=2).
!
!  08/01/97 (Zonghui Huo)
!  Added Kain-fritsch cumulus parameterization scheme.
!
!  11/05/97 (D. Weber)
!  Added temxy5 array for use in the bottom boundary condition for
!  pertrubation pressure (hydrostatic).
!
!  11/06/97 (D. Weber)
!  Added three additional levels to the mapfct array.  The three
!  levels (4,5,6) represent the inverse of the first three in order.
!  The inverse map factors are computed to improve efficiency.
!
!  4/15/1998 (Donghai Wang)
!  Added the source terms to the right hand terms of the qc,qr,qi,qs
!  equations due to the K-F cumulus parameterization.
!
!  4/15/1998 (Donghai Wang)
!  Added the running average vertical velocity (array w0avg)
!  for the K-F cumulus parameterization scheme.
!
!  9/15/1998 (D. Weber)
!  Added ptbari, pbari (inverse) for use in optimizing the code.
!
!  8/31/1998 (K. Brewster)
!  Added call to ININUDGE to initialize nudging terms.
!
!  11/18/1998 (K. Brewster)
!  Changed pibar to ppi (full pi) and moved initialization.
!
!  12/09/1998 (Donghai Wang)
!  Added the snow cover.
!
!  07/10/2001 (K. Brewster)
!  Added increment arrays to argument list and to call to ininudge. 
!
!  13 March 2002 (Eric Kemp)
!  Added arrays for WRF BMJ cumulus scheme.
!
!  April 2002 (Fanyou Kong)
!  Added WRF new Kain-Fritsch (April 2002 version: KF_ETA) scheme
!  initialization (lookup table)
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    mptr     Grid identifier.
!    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
!    nxndg    Number of x grid points for nudging (1 or nx)
!    nyndg    Number of y grid points for nudging (1 or ny)
!    nzndg    Number of z grid points for nudging (1 or nz)
!
!  OUTPUT:
!
!    u        x-component of velocity at all time levels (m/s).
!    v        y-component of velocity at all time levels (m/s).
!    w        z-component of velocity at all time levels (m/s).
!    wcont    Contravariant vertical velocity (m/s)
!    ptprt    Perturbation potential temperature at all time levels (K)
!    pprt     Perturbation pressure at all time levels (Pascal)
!    qv       Water vapor specific humidity at all time levels (kg/kg)
!    qc       Cloud water mixing ratio at all time levels (kg/kg)
!    qr       Rainwater mixing ratio at all time levels (kg/kg)
!    qi       Cloud ice mixing ratio at all time levels (kg/kg)
!    qs       Snow mixing ratio at all time levels (kg/kg)
!    qh       Hail mixing ratio at all time levels (kg/kg)
!    tke      Turbulent Kinetic Energy ((m/s)**2)
!
!    udteb    Time tendency of u field at east boundary (m/s**2)
!    udtwb    Time tendency of u field at west boundary (m/s**2)
!    udtnb    Time tendency of u field at north boundary (m/s**2)
!    udtsb    Time tendency of u field at south boundary (m/s**2)
!
!    vdteb    Time tendency of v field at east boundary (m/s**2)
!    vdtwb    Time tendency of v field at west boundary (m/s**2)
!    vdtnb    Time tendency of v field at north boundary (m/s**2)
!    vdtsb    Time tendency of v field at south boundary (m/s**2)
!
!    wdteb    Time tendency of w field at east boundary (m/s**2)
!    wdtwb    Time tendency of w field at west boundary (m/s**2)
!    wdtnb    Time tendency of w field at north boundary (m/s**2)
!    wdtsb    Time tendency of w field at south boundary (m/s**2)
!
!    pdteb    Time tendency of pprt field at east boundary (PASCAL/s)
!    pdtwb    Time tendency of pprt field at west boundary (PASCAL/s)
!    pdtnb    Time tendency of pprt field at north boundary (PASCAL/s)
!    pdtsb    Time tendency of pprt field at south boundary (PASCAL/s)
!
!    sdteb    Time tendency of a scalar field at east boundary (m/s**2)
!    sdtwb    Time tendency of a scalar field at west boundary (m/s**2)
!    sdtnb    Time tendency of a scalar field at north boundary (m/s**2)
!    sdtsb    Time tendency of a scalar field at south boundary (m/s**2)
!
!    ubar     Base state x-velocity component (m/s)
!    vbar     Base state y-velocity component (m/s)
!    ptbar    Base state potential temperature (K)
!    pbar     Base state pressure (Pascal)
!    ptbari   Inverse Base state potential temperature (K)
!    pbari    Inverse Base state pressure (Pascal)
!    rhostr   Base state density rhobar times j3 (kg/m**3)
!    rhostri  Inverse base state density rhobar times j3 (kg/m**3)
!    qvbar    Base state water vapor specific humidity (kg/kg).
!    ppi      Exner function
!    csndsq   Sound wave speed squared.
!
!    x        x-coordinate of grid points in computational space (m)
!    y        y-coordinate of grid points in computational space (m)
!    z        z-coordinate of grid points in computational space (m)
!    zp       Vertical coordinate of grid points in physical space (m)
!    hterain  Terrain height (m)
!    mapfct   Map factors at scalar, u and v points
!
!    j1       Coordinate transformation Jacobian -d(zp)/dx
!    j2       Coordinate transformation Jacobian -d(zp)/dy
!    j3       Coordinate transformation Jacobian  d(zp)/dz
!    aj3x     Avgx of the coordinate transformation Jacobian  d(zp)/dz
!    aj3y     Avgy of the coordinate transformation Jacobian  d(zp)/dz
!    aj3z     Avgz of the coordinate transformation Jacobian  d(zp)/dz
!    j3inv    Inverse of the coordinate transformation j3
!
!    trigs1   Array containing pre-computed trig function for fftopt=1.
!    trigs2   Array containing pre-computed trig function for fftopt=1.
!    ifax1    Array containing the factors of nx for fftopt=1.
!    ifax2    Array containing the factors of ny for fftopt=1.
!
!    vwork1   2-D work array for fftopt=2 option.
!    vwork2   2-D work array for fftopt=2 option.
!    wsave1   Work array for fftopt=2 option.
!    wsave2   Work array for fftopt=2 option.
!
!    sinlat   Sin of latitude at each grid point
!
!    kmh      Horizontal turb. mixing coef. for momentum ( m**2/s )
!    kmv      Vertical turb. mixing coef. for momentum ( m**2/s )
!    rprntl   Reciprocal of Prandtl number
!
!    soiltyp  Soil type
!    stypfrct  Soil type fraction
!    vegtyp   Vegetation type
!    lai      Leaf Area Index
!    roufns   Surface roughness
!    veg      Vegetation fraction
!
!    tsfc     Skin temperature at the ground or ocean surface (K)
!    tsoil    Deep soil temperature (K) (in deep 1 m layer)
!    wetsfc   Surface soil moisture
!    wetdp    Deep soil moisture
!    wetcanp  Canopy water amount
!    ptsfc    Ground surface potential temperature (K)
!    qvsfc    Effective S.H. at sfc.
!
!    ptcumsrc Source term in pt-equation due to cumulus parameterization
!    qcumsrc  Source term in water equations due to cumulus parameterization
!
!    nca      Counter for CAPE release in the Kain-Fritsch scheme
!    kfraincv K-F convective precipitation rate
!    cldefi   BMJ cloud efficiency
!    xland    BMJ land/sea mask
!    bmjraincv BMJ convective precipitation amount
!
!    radfrc   Radiation forcing (K/s)
!    radsw    Solar radiation reaching the surface
!    rnflx    Net absorbed radiation by the surface
!
!    raing    Grid supersaturation rain
!    rainc    Cumulus convective rain
!
!    usflx    Surface flux of u-momentum (kg/(m*s**2))
!    vsflx    Surface flux of v-momentum (kg/(m*s**2))
!    ptsflx   Surface heat flux (K*kg/(m**2 * s ))
!    qvsflx   Surface moisture flux of (kg/(m**2 * s))
!
!    temxy1   Temporary work array
!
!    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.
!    tem10    Temporary work array.
!    tem11    Temporary work array.
!    tem12    Temporary work array.
!    tem13    Temporary work array.
!    tem14    Temporary work array.
!    tem15    Temporary work array.
!    tem16    Temporary work array.
!    tem17    Temporary work array.
!    tem18    Temporary work array.
!    tem19    Temporary work array.
!    tem20    Temporary work array.
!    tem21    Temporary work array.
!    tem22    Temporary work array.
!    tem23    Temporary work array.
!    tem24    Temporary work array.
!    tem25    Temporary work array.
!    tem26    Temporary work array.
!
!    tem1_0   Temporary work array.
!    tem2_0   Temporary work array.
!    tem3_0   Temporary work array.
!
!-----------------------------------------------------------------------
!

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

  INTEGER :: mptr              ! Grid identifier

  INCLUDE 'timelvls.inc'

  INTEGER :: nx,ny,nz          ! The number of grid points in 3 directions
  INTEGER :: nxndg,nyndg,nzndg ! The number of grid points in 3 directions

  REAL :: u     (nx,ny,nz,nt)  ! Total u-velocity (m/s).

  REAL :: v     (nx,ny,nz,nt)  ! Total v-velocity (m/s).
  REAL :: w     (nx,ny,nz,nt)  ! Total w-velocity (m/s).
  REAL :: wcont (nx,ny,nz)     ! Contravariant vertical velocity (m/s)
  REAL :: ptprt (nx,ny,nz,nt)  ! Perturbation potential temperature
                               ! from that of base state atmosphere (Kelvin).
  REAL :: pprt  (nx,ny,nz,nt)  ! Perturbation pressure from that
                               ! of base state atmosphere (Pascal).
  REAL :: qv    (nx,ny,nz,nt)  ! Water vapor specific humidity (kg/kg).
  REAL :: qc    (nx,ny,nz,nt)  ! Cloud water mixing ratio (kg/kg).
  REAL :: qr    (nx,ny,nz,nt)  ! Rain water mixing ratio (kg/kg).
  REAL :: qi    (nx,ny,nz,nt)  ! Cloud ice mixing ratio (kg/kg).
  REAL :: qs    (nx,ny,nz,nt)  ! Snow mixing ratio (kg/kg).
  REAL :: qh    (nx,ny,nz,nt)  ! Hail mixing ratio (kg/kg).
  REAL :: tke   (nx,ny,nz,nt)  ! Turbulent Kinetic Energy ((m/s)**2)

  REAL :: udteb (ny,nz)        ! T-tendency of u at e-boundary (m/s**2)
  REAL :: udtwb (ny,nz)        ! T-tendency of u at w-boundary (m/s**2)
  REAL :: udtnb (nx,nz)        ! T-tendency of u at n-boundary (m/s**2)
  REAL :: udtsb (nx,nz)        ! T-tendency of u at s-boundary (m/s**2)

  REAL :: vdteb (ny,nz)        ! T-tendency of v at e-boundary (m/s**2)
  REAL :: vdtwb (ny,nz)        ! T-tendency of v at w-boundary (m/s**2)
  REAL :: vdtnb (nx,nz)        ! T-tendency of v at n-boundary (m/s**2)
  REAL :: vdtsb (nx,nz)        ! T-tendency of v at s-boundary (m/s**2)

  REAL :: wdteb (ny,nz)        ! T-tendency of w at e-boundary (m/s**2)
  REAL :: wdtwb (ny,nz)        ! T-tendency of w at w-boundary (m/s**2)
  REAL :: wdtnb (nx,nz)        ! T-tendency of w at n-boundary (m/s**2)
  REAL :: wdtsb (nx,nz)        ! T-tendency of w at s-boundary (m/s**2)

  REAL :: pdteb (ny,nz)        ! T-tendency of pprt at e-boundary (PASCAL/s)
  REAL :: pdtwb (ny,nz)        ! T-tendency of pprt at w-boundary (PASCAL/s)
  REAL :: pdtnb (nx,nz)        ! T-tendency of pprt at n-boundary (PASCAL/s)
  REAL :: pdtsb (nx,nz)        ! T-tendency of pprt at s-boundary (PASCAL/s)

  REAL :: sdteb (ny,nz)        ! T-tendency of w at e-boundary (m/s**2)
  REAL :: sdtwb (ny,nz)        ! T-tendency of w at w-boundary (m/s**2)
  REAL :: sdtnb (nx,nz)        ! T-tendency of w at n-boundary (m/s**2)
  REAL :: sdtsb (nx,nz)        ! T-tendency of w at s-boundary (m/s**2)

  REAL :: ubar  (nx,ny,nz)     ! Base state u-velocity (m/s).
  REAL :: vbar  (nx,ny,nz)     ! Base state v-velocity (m/s).
  REAL :: ptbar (nx,ny,nz)     ! Base state potential temperature (K)
  REAL :: pbar  (nx,ny,nz)     ! Base state pressure (Pascal).
  REAL :: ptbari(nx,ny,nz)     ! Inverse Base state pot. temperature (K)
  REAL :: pbari (nx,ny,nz)     ! Inverse Base state pressure (Pascal).
  REAL :: rhostr(nx,ny,nz)     ! Base state density rhobar times j3.
  REAL :: rhostri(nx,ny,nz)    ! Inverse base state density rhobar times j3.
  REAL :: qvbar (nx,ny,nz)     ! Base state water vapor specific humidity
                               ! (kg/kg).
  REAL :: ppi   (nx,ny,nz)     ! Exner function.
  REAL :: csndsq(nx,ny,nz)     ! Sound wave speed squared.

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

  REAL :: mapfct(nx,ny,8)      ! Map factors at scalar, u and v points

  REAL :: j1    (nx,ny,nz)     ! Coordinate transformation Jacobian -d(zp)/dx.
  REAL :: j2    (nx,ny,nz)     ! Coordinate transformation Jacobian -d(zp)/dy.
  REAL :: j3    (nx,ny,nz)     ! Coordinate transformation Jacobian  d(zp)/dz.
  REAL :: aj3x  (nx,ny,nz)     ! Coordinate transformation Jacobian defined
                               ! as d( zp )/d( z ) AVERAGED IN THE X-DIR.
  REAL :: aj3y  (nx,ny,nz)     ! Coordinate transformation Jacobian defined
                               ! as d( zp )/d( z ) AVERAGED IN THE Y-DIR.
  REAL :: aj3z  (nx,ny,nz)     ! Coordinate transformation Jacobian defined
                               ! as d( zp )/d( z ) AVERAGED IN THE Z-DIR.
  REAL :: j3inv (nx,ny,nz)     ! Inverse of J3
  REAL :: trigs1(3*(nx-1)/2+1) ! Array containing pre-computed trig
                               ! function for fftopt=1.
  REAL :: trigs2(3*(ny-1)/2+1) ! Array containing pre-computed trig
                               ! function for fftopt=1.
  INTEGER :: ifax1(13)         ! Array containing the factors of nx for
                               ! fftopt=1.
  INTEGER :: ifax2(13)         ! Array containing the factors of ny for
                               ! fftopt=1.

  REAL :: vwork1 (nx+1,ny+1)   ! 2-D work array for fftopt=2.
  REAL :: vwork2 (ny,nx+1)     ! 2-D work array for fftopt=2.
  REAL :: wsave1 (3*(ny-1)+15) ! Work array for fftopt=2.
  REAL :: wsave2 (3*(nx-1)+15) ! Work array for fftopt=2.

  REAL :: sinlat(nx,ny)        ! Sin of latitude at each grid point

  REAL :: kmh   (nx,ny,nz)     ! Horizontal turb. mixing coef. for
                               ! momentum. ( m**2/s )
  REAL :: kmv   (nx,ny,nz)     ! Vertical turb. mixing coef. for
                               ! momentum. ( m**2/s )
  REAL :: rprntl(nx,ny,nz)     ! Reciprocal of Prandtl number


  INTEGER :: nstyps                  ! Number of soil type
  INTEGER :: soiltyp(nx,ny,nstyps)   ! Soil types at grids
  REAL :: stypfrct(nx,ny,nstyps)  ! Fraction of soil types
  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)    ! Temperature at ground (K)
                                     ! (in top 1 cm layer)
  REAL :: tsoil  (nx,ny,0:nstyps)    ! Deep soil temperature (K)
                                     ! (in deep 1 m layer)
  REAL :: wetsfc (nx,ny,0:nstyps)    ! Surface soil moisture
                                     ! in the top 1 cm layer
  REAL :: wetdp  (nx,ny,0:nstyps)    ! Deep soil moisture
                                     ! in the deep 1 m layer
  REAL :: wetcanp(nx,ny,0:nstyps)    ! Canopy water amount
  REAL :: snowdpth(nx,ny)            ! Snow depth (m)

  REAL :: qvsfc  (nx,ny,0:nstyps)    ! Effective qv at sfc.
  REAL :: ptsfc  (nx,ny)             ! Ground surface potential
                                     ! temperature (K)

  REAL :: ptcumsrc(nx,ny,nz)   ! Source term in pt-equation due
                               ! to cumulus parameterization
  REAL :: qcumsrc(nx,ny,nz,5)  ! Source term in water equations due
                               ! to cumulus parameterization:
                               ! qcumsrc(1,1,1,1) for qv equation
                               ! qcumsrc(1,1,1,2) for qc equation
                               ! qcumsrc(1,1,1,3) for qr equation
                               ! qcumsrc(1,1,1,4) for qi equation
                               ! qcumsrc(1,1,1,5) for qs equation
  REAL :: w0avg(nx,ny,nz)      ! a closing running average vertical
                               ! velocity in 10min for K-F scheme
  REAL :: kfraincv(nx,ny)      ! K-F convective rainfall (cm)
  INTEGER :: nca(nx,ny)        ! K-F counter for CAPE release

!EMK BMJ   
  REAL,INTENT(INOUT) :: cldefi(nx,ny)      ! BMJ cloud efficiency
  REAL,INTENT(INOUT) :: xland(nx,ny)       ! BMJ land mask
                                           !   (1.0 = land, 2.0 = sea)
  REAL,INTENT(INOUT) :: bmjraincv(nx,ny)   ! BMJ convective rainfall (cm)
!EMK END

  REAL :: radfrc(nx,ny,nz)     ! Radiation forcing (K/s)
  REAL :: radsw(nx,ny)         ! Solar radiation reacing the surface
  REAL :: rnflx(nx,ny)         ! Net absorbed radiation by the surface

  REAL :: raing(nx,ny)         ! Grid supersaturation rain
  REAL :: rainc(nx,ny)         ! Cumulus convective rain
  REAL :: prcrate(nx,ny,4)     ! precipitation rates (kg/(m**2*s))
                               ! prcrate(1,1,1) = total precipitation rate
                               ! prcrate(1,1,2) = grid scale precip. rate
                               ! prcrate(1,1,3) = cumulus precip. rate
                               ! prcrate(1,1,4) = microphysics precip. rate

  REAL :: usflx (nx,ny)        ! Surface flux of u-momentum (kg/(m*s**2))
  REAL :: vsflx (nx,ny)        ! Surface flux of v-momentum (kg/(m*s**2))
  REAL :: ptsflx(nx,ny)        ! Surface heat flux (K*kg/(m*s**2))
  REAL :: qvsflx(nx,ny)        ! Surface moisture flux (kg/(m**2*s))

  INTEGER :: exbcbufsz         ! EXBC buffer size
  REAL :: exbcbuf( exbcbufsz ) ! EXBC buffer array
  REAL :: bcrlx(nx,ny)         ! EXBC relaxation coefficients

  REAL :: uincr(nxndg,nyndg,nzndg)      ! Analysis increment for u
  REAL :: vincr(nxndg,nyndg,nzndg)      ! Analysis increment for v
  REAL :: wincr(nxndg,nyndg,nzndg)      ! Analysis increment for w
  REAL :: pincr(nxndg,nyndg,nzndg)      ! Analysis increment for p
  REAL :: ptincr(nxndg,nyndg,nzndg)     ! Analysis increment for pt
  REAL :: qvincr(nxndg,nyndg,nzndg)     ! Analysis increment for qv
  REAL :: qcincr(nxndg,nyndg,nzndg)     ! Analysis increment for qc
  REAL :: qrincr(nxndg,nyndg,nzndg)     ! Analysis increment for qr
  REAL :: qiincr(nxndg,nyndg,nzndg)     ! Analysis increment for qi
  REAL :: qsincr(nxndg,nyndg,nzndg)     ! Analysis increment for qs
  REAL :: qhincr(nxndg,nyndg,nzndg)     ! Analysis increment for qh

  REAL :: temxy1(nx,ny)        ! Temporary work array

  REAL :: tem1  (nx,ny,nz)     ! Temporary work array.
  REAL :: tem2  (nx,ny,nz)     ! Temporary work array.
  REAL :: tem3  (nx,ny,nz)     ! Temporary work array.
  REAL :: tem4  (nx,ny,nz)     ! Temporary work array.
  REAL :: tem5  (nx,ny,nz)     ! Temporary work array.
  REAL :: tem6  (nx,ny,nz)     ! Temporary work array.
  REAL :: tem7  (nx,ny,nz)     ! Temporary work array.
  REAL :: tem8  (nx,ny,nz)     ! Temporary work array.
  REAL :: tem9  (nx,ny,nz)     ! Temporary work array.
  REAL :: tem10 (nx,ny,nz)     ! Temporary work array.
  REAL :: tem11 (nx,ny,nz)     ! Temporary work array.
  REAL :: tem12 (nx,ny,nz)     ! Temporary work array.
  REAL :: tem13 (nx,ny,nz)     ! Temporary work array.
  REAL :: tem14 (nx,ny,nz)     ! Temporary work array.
  REAL :: tem15 (nx,ny,nz)     ! Temporary work array.
  REAL :: tem16 (nx,ny,nz)     ! Temporary work array.
  REAL :: tem17 (nx,ny,nz)     ! Temporary work array.
  REAL :: tem18 (nx,ny,nz)     ! Temporary work array.
  REAL :: tem19 (nx,ny,nz)     ! Temporary work array.
  REAL :: tem20 (nx,ny,nz)     ! Temporary work array.
  REAL :: tem21 (nx,ny,nz)     ! Temporary work array.
  REAL :: tem22 (nx,ny,nz)     ! Temporary work array.
  REAL :: tem23 (nx,ny,nz)     ! Temporary work array.
  REAL :: tem24 (nx,ny,nz)     ! Temporary work array.
  REAL :: tem25 (nx,ny,nz)     ! Temporary work array.
  REAL :: tem26 (nx,ny,nz)     ! Temporary work array.

  REAL :: tem1_0(0:nx,0:ny,0:nz)     ! Temporary work array.
  REAL :: tem2_0(0:nx,0:ny,0:nz)     ! Temporary work array.
  REAL :: tem3_0(0:nx,0:ny,0:nz)     ! Temporary work array.
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'globcst.inc'
  INCLUDE 'bndry.inc'
  INCLUDE 'indtflg.inc'
  INCLUDE 'phycst.inc'
  INCLUDE 'exbc.inc'
  INCLUDE 'nudging.inc'
!
!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------
!
  INTEGER :: i,j
  INTEGER :: ireturn
  REAL :: tem

!EMK BMJ   
  LOGICAL :: restart
  
  INTEGER,PARAMETER :: vegwaterflag = 14
  INTEGER,PARAMETER :: xland_waterflag = 2
  INTEGER,PARAMETER :: xland_landflag = 1
!EMK END

!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  mgrid = mptr

  grdin = 0
  basin = 0
  varin = 0
  mstin = 0
  rainin= 0
  prcin = 0
  icein = 0
  trbin = 0
  sfcin = 0
  landin= 0
  radin = 0
  flxin = 0
!wdt update - init0 no longer necessary (arrays set to 0 after allocation)
!
!-----------------------------------------------------------------------
!
!  INITialize model array VARiables.
!
!-----------------------------------------------------------------------
!
!EMK BMJ
  CALL initgrdvar(nx,ny,nz,nt,nstyps,exbcbufsz,                         &
                  x,y,z,zp,hterain,mapfct,                              &
                  j1,j2,j3,aj3x,aj3y,aj3z,j3inv,                        &
                  u,v,w,wcont,ptprt,pprt,qv,qc,qr,qi,qs,qh,tke,         &
                  udteb, udtwb, vdtnb, vdtsb,                           &
                  pdteb,pdtwb ,pdtnb ,pdtsb,                            &
                  trigs1,trigs2,ifax1,ifax2,                            &
                  wsave1,wsave2,vwork1,vwork2,                          &
                  ubar,vbar,ptbar,pbar,ptbari,pbari,                    &
                  rhostr,rhostri,qvbar,ppi,csndsq,                      &
                  soiltyp,stypfrct,vegtyp,lai,roufns,veg,               &
                  tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth,qvsfc,       &
                  ptcumsrc,qcumsrc,w0avg,nca,kfraincv,                  &
                  cldefi,xland,bmjraincv,                               &
                  raing,rainc,prcrate,exbcbuf,                          &
                  radfrc,radsw,rnflx,                                   &
                  usflx,vsflx,ptsflx,qvsflx,                            &
                  tem1,tem2,tem3,tem4,tem5,                             &
                  tem6,tem7,tem8,tem9)

  DO j=1,ny-1
    DO i=1,nx-1
      tem = 0.5 * ( pprt(i,j,1,2)+pbar(i,j,1)                           &
                  + pprt(i,j,2,2)+pbar(i,j,2) )
      ptsfc(i,j)=tsfc(i,j,0)*(p0/tem)**rddcp
    END DO
  END DO

!
!-----------------------------------------------------------------------
!
!  Calculate the sin of the lattitude of each grid point, to be used
!  in the calculation of latitude-dependent Coriolis parameters.
!
!-----------------------------------------------------------------------
!
  CALL gtsinlat(nx,ny,x,y, sinlat, tem1,tem2, tem3)
!
!-----------------------------------------------------------------------
!
!  Initialize arrays that store the lookup table data.
!
!-----------------------------------------------------------------------
!

  CALL initlktb

!
!-----------------------------------------------------------------------
!
!  Initialize the external boundary data array.
!
!-----------------------------------------------------------------------
!
  IF( lbcopt == 2 .AND. mptr == 1 ) THEN

    ireturn = 0

    CALL extbdtini(nx,ny,nz,                                            &
                   u,v,w,ptprt,pprt,                                    &
                   qv,qc,qr,qi,qs,qh,ptbar,pbar,                        &
                   exbcbuf(nu0exb),  exbcbuf(nv0exb),                   &
                   exbcbuf(nw0exb),  exbcbuf(npt0exb),                  &
                   exbcbuf(npr0exb), exbcbuf(nqv0exb),                  &
                   exbcbuf(nqc0exb), exbcbuf(nqr0exb),                  &
                   exbcbuf(nqi0exb), exbcbuf(nqs0exb),                  &
                   exbcbuf(nqh0exb), exbcbuf(nudtexb),                  &
                   exbcbuf(nvdtexb), exbcbuf(nwdtexb),                  &
                   exbcbuf(nptdtexb),exbcbuf(nprdtexb),                 &
                   exbcbuf(nqvdtexb),exbcbuf(nqcdtexb),                 &
                   exbcbuf(nqrdtexb),exbcbuf(nqidtexb),                 &
                   exbcbuf(nqsdtexb),exbcbuf(nqhdtexb),                 &
                   bcrlx,                                               &
                   tem1,tem2,tem3,tem4,tem5,tem6,                       &
                   tem7,tem8,tem9,tem10,tem11,ireturn)

    IF( ireturn == 1 ) THEN
      WRITE (6,'(a/a)')                                                 &
          'Can not find the external boundary data. Dump the',          &
          'history file and restart file and then STOP the model.'
      CALL arpsstop('arpsstop called from initial with ext boundary file',1)
    ELSE IF( ireturn == 2 ) THEN
      WRITE (6,'(a/a)')                                                 &
          'Can not open the external boundary data. Dump the history',  &
          'file and restart file and then STOP the model.'
      CALL arpsstop('arpsstop called from initial with opening ext      &
            & boundary file ',1)
    ELSE IF( ireturn == 3 ) THEN
      WRITE (6,'(a/a)')                                                 &
          'Read errors in the external boundary data file. Dump the',   &
          'history file and restart file and then STOP the model.'
      CALL arpsstop('arpsstop called from initial with reading ext      &
            & boundary file ',1)
    ELSE IF( ireturn /= 0 ) THEN
      WRITE (6,'(a/a)')                                                 &
          'Other errors in getting the external boundary data. Dump the', &
          'history file and restart file and then STOP the model.'
      CALL arpsstop('arpsstop called from initial with the ext      &
            & boundary file ',1)
    END IF

  END IF

!
!-----------------------------------------------------------------------
!
!  Initialize nudging assimilation variables
!
!-----------------------------------------------------------------------
!
  IF( nudgopt > 0 )                                                     &
    CALL ininudge(nxndg,nyndg,nzndg,                                    &
                  uincr,vincr,wincr,pincr,ptincr,qvincr,                &
                  qcincr,qrincr,qiincr,qsincr,qhincr,ireturn)

!-----------------------------------------------------------------------
!
! Initialize KF_ETA arrays and look-up tables.
!
!-----------------------------------------------------------------------

!KFY KF_ETA
  IF (cnvctopt == 5) THEN

!...Now initialize KF_ETA look-up tables

    IF (initopt == 2) THEN
      restart = .TRUE.
    ELSE
      restart = .FALSE.
    END IF

    CALL interface_wrf_kfinit(nx,ny,nz,nca,restart)
  END IF ! IF (cnvctopt == 5) THEN

!KFY KF_ETA END

!-----------------------------------------------------------------------
!
! Initialize BMJ arrays and look-up tables.
!
!-----------------------------------------------------------------------

!EMK BMJ
  IF (cnvctopt == 4) THEN
    DO j = 1,ny-1
      DO i = 1,nx-1
        IF (vegtyp(i,j) == vegwaterflag) THEN
          xland(i,j) = xland_waterflag
        ELSE
          xland(i,j) = xland_landflag
        END IF
      END DO ! DO i = 1,nx
    END DO ! DO j = 1,ny

!...Now initialize BMJ look-up tables

    IF (initopt == 2) THEN
      restart = .TRUE.
    ELSE
      restart = .FALSE.
    END IF
    
    CALL interface_wrf_bmjinit(nx,ny,nz,cldefi,restart)
  END IF ! IF (cnvctopt == 4) THEN
          
!EMK END

!
!-----------------------------------------------------------------------
!
!  End of model initialization.
!
!-----------------------------------------------------------------------
!
  RETURN
END SUBROUTINE initial
!
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE INIGRD                     ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE inigrd(nx,ny,nz,                                             & 3,18
           x,y,z,zp,hterain,mapfct,j1,j2,j3,zp1d,dzp1d,tem1)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Initialize the model grid variables.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Ming Xue
!  3/17/1991.
!
!  MODIFICATION HISTORY:
!
!  5/02/92 (M. Xue)
!  Added full documentation.
!
!  5/03/92 (M. Xue)
!  Further documentation.
!
!  4/10/93 (M. Xue & Hao Jin)
!  Add the terrain.
!
!  9/10/94 (Weygandt & Y. Lu)
!  Cleaned up documentation.
!
!  12/22/94 (Yuhe Liu)
!  Changed variable tuning, which was hard wired inside this
!  subroutine, to strhtune, which is input from namelist input file.
!
!  1/22/96 (Donghai Wang & Yuhe Liu)
!  Added the map projection factor to ARPS governing equations.
!
!  11/06/97 (D. Weber)
!  Added three additional levels to the mapfct array.  The three
!  levels (4,5,6) represent the inverse of the first three in order.
!  The inverse map factors are computed to improve efficiency.
!
!  2000/04/11 (Gene Bassett)
!  Only update the terrain data with a call to BCS2D for ternopt=1.
!
!-----------------------------------------------------------------------
!
!  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
!
!  OUTPUT:
!
!    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 (m)
!    zp       Vertical coordinate of grid points in physical space
!             (m)
!    hterain  Terrain height (m)
!    mapfct   Map factors at scalar, u and v points
!
!    j1       Coordinate transformation Jacobian -d(zp)/dx
!    j2       Coordinate transformation Jacobian -d(zp)/dy
!    j3       Coordinate transformation Jacobian  d(zp)/dz
!
!  WORK ARRAY:
!
!    zp1d     Temporary work array.
!    dzp1d    Temporary work array.
!    tem1     Temporary work array.
!
!-----------------------------------------------------------------------
!

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

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

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

  REAL :: j1(nx,ny,nz)    ! Coordinate transformation Jacobian
                          ! -d(zp)/dx.
  REAL :: j2(nx,ny,nz)    ! Coordinate transformation Jacobian
                          ! -d(zp)/dy.
  REAL :: j3(nx,ny,nz)    ! Coordinate transformation Jacobian
                          ! d(zp)/dz.

  REAL :: hterain(nx,ny)  ! Terrain height.

  REAL :: mapfct(nx,ny,8)      ! Map factors at scalar, u and v points

  REAL :: zp1d (nz)       ! Temporary array
  REAL :: dzp1d(nz)       ! Temporary array
  REAL :: tem1(nx,ny,nz)  ! Temporary array
!
!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------
!
  INTEGER :: i,j,k
  REAL :: xs, ys, zs, radnd, pi2
  REAL :: zflat1,z1,z2

  REAL :: zpmax

  INTEGER :: mptag
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'globcst.inc'
  INCLUDE 'grid.inc'          ! Grid & map parameters.
  INCLUDE 'phycst.inc'
  INCLUDE 'bndry.inc'
  INCLUDE 'mp.inc'            ! Message passing parameters.

!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
!-----------------------------------------------------------------------
!
!  Define a uniform model grid.
!
!-----------------------------------------------------------------------
!
  CALL setgrd ( nx,ny, x,y )
!
!-----------------------------------------------------------------------
!
!  Set map factor at scalar, u, and v points
!
!-----------------------------------------------------------------------
!
  IF ( mpfctopt /= 0 ) THEN

    DO j=1,ny-1
      DO i=1,nx-1
        xs = 0.5*(x(i)+x(i+1))
        ys = 0.5*(y(j)+y(j+1))
        CALL xytomf( 1,1,xs,ys,mapfct(i,j,1) )
        IF(maptest == 1)THEN   ! set up a symmetric field...
          mapfct(i,j,1) = 1.0 + ABS((xs-0.5*(x(1)+x(nx)))               &
                                   /(x(nx)-x(1))                        &
                                   +(ys-0.5*(y(1)+y(ny)))               &
                                   /(y(ny)-y(1)))
        END IF
        mapfct(i,j,4) = 1.0/mapfct(i,j,1)
        mapfct(i,j,7) = mapfct(i,j,1)*mapfct(i,j,1)
        mapfct(i,j,8) = 0.25*mapfct(i,j,1)  ! for use in sovlwpim
                                            ! and wcontra...
      END DO
    END DO

    DO j=1,ny-1
      DO i=1,nx
        ys = 0.5*(y(j)+y(j+1))
        CALL xytomf( 1,1,x(i),ys,mapfct(i,j,2) )
        IF(maptest == 1)THEN   ! set up a symmetric field...
          mapfct(i,j,2) = 1.0 + ABS((x(i)-0.5*(x(1)+x(nx)))             &
                                   /(x(nx)-x(1))                        &
                                   +(ys-0.5*(y(1)+y(ny)))               &
                                   /(y(ny)-y(1)))
        END IF
        mapfct(i,j,5) = 1.0/mapfct(i,j,2)
      END DO
    END DO

    DO j=1,ny
      DO i=1,nx-1
        xs = 0.5*(x(i)+x(i+1))
        CALL xytomf( 1,1,xs,y(j),mapfct(i,j,3) )
        IF(maptest == 1)THEN   ! set up a symmetric field...
          mapfct(i,j,3) = 1.0 + ABS((xs-0.5*(x(1)+x(nx)))               &
                                   /(x(nx)-x(1))                        &
                                   +(y(j)-0.5*(y(1)+y(ny)))             &
                                   /(y(ny)-y(1)))
        END IF
        mapfct(i,j,6) = 1.0/mapfct(i,j,3)
      END DO
    END DO

  ELSE

    DO k=1,7
      DO j=1,ny
        DO i=1,nx
          mapfct(i,j,k) = 1.0
          mapfct(i,j,8) = 0.25
        END DO
      END DO
    END DO

  END IF
!
!-----------------------------------------------------------------------
!
!  Print map factor at scalar, u, and v points for the purpose of
!  debug
!
!-----------------------------------------------------------------------
!
!   CALL getunit( mpunit )
!   CALL asnctl ('NEWLOCAL', 1, ierr)
!   CALL asnfile('mapfactor.data', '-F f77 -N ieee', ierr)
!
!   OPEN (mpunit,file='mapfactor.data',form='unformatted')
!
!   CALL edgfill(mapfct(1,1,1),1,nx,1,nx-1, 1,ny,1,ny-1, 1,1,1,1)
!   write(mpunit) ((mapfct(i,j,1),i=1,nx),j=1,ny)
!
!   CALL edgfill(mapfct(1,1,2),1,nx,1,nx,   1,ny,1,ny-1, 1,1,1,1)
!   write(mpunit) ((mapfct(i,j,2),i=1,nx),j=1,ny)
!
!   CALL edgfill(mapfct(1,1,3),1,nx,1,nx-1, 1,ny,1,ny,   1,1,1,1)
!   write(mpunit) ((mapfct(i,j,3),i=1,nx),j=1,ny)
!
!   CLOSE( mpunit )
!   CALL retunit( mpunit )
!
!-----------------------------------------------------------------------
!
!  End of debug print.
!
!-----------------------------------------------------------------------
!
  DO k=1,nz
    z(k) = zrefsfc + (k-2) * dz
  END DO
!
!-----------------------------------------------------------------------
!
!  Specify the terrain
!
!-----------------------------------------------------------------------
!
  IF( ternopt == 0 ) THEN     ! No terrain, the ground is flat

    DO j=1,ny-1
      DO i=1,nx-1
        hterain(i,j) = zrefsfc
      END DO
    END DO

  ELSE IF( ternopt == 1 ) THEN  ! Bell-shaped mountain
!
!-----------------------------------------------------------------------
!
!  Define the bell-shaped mountain
!
!-----------------------------------------------------------------------
!
    pi2 = 1.5707963267949

    DO j=1,ny-1
      DO i=1,nx-1
        xs = (x(i)+x(i+1))*0.5
        ys = (y(j)+y(j+1))*0.5
        zs = z(2)

        IF( mntwidy < 0.0 .OR. runmod == 2 ) THEN ! 2-d terrain in x-z plane.

          radnd = 1.0+((xs-mntctrx)/mntwidx)**2

        ELSE IF( mntwidx < 0.0 .OR. runmod == 3 ) THEN ! 2-d terrain in y-z plane.

          radnd = 1.0+((ys-mntctry)/mntwidy)**2

        ELSE                             ! 3-d terrain

          radnd = 1.0+((xs-mntctrx)/mntwidx)**2                         &
                 +((ys-mntctry)/mntwidy)**2
        END IF

        hterain(i,j) = zrefsfc + hmount/radnd

      END DO
    END DO
!
!-----------------------------------------------------------------------
!
!  Make sure that the terrain satisfies the specified boundary
!  conditions.
!
!-----------------------------------------------------------------------
!
!call test_dump (hterain,"XXXinit_hterain",nx,ny,1,0,0)
    IF (mp_opt > 0) THEN
      CALL acct_interrupt(mp_acct)
      CALL mpsend1dew(hterain,nx,ny,ebc,wbc,0,mptag,tem1)
      CALL mprecv1dew(hterain,nx,ny,ebc,wbc,0,mptag,tem1)
      CALL mpsend1dns(hterain,nx,ny,nbc,sbc,0,mptag,tem1)
      CALL mprecv1dns(hterain,nx,ny,nbc,sbc,0,mptag,tem1)
    END IF
!call test_dump (hterain,"XXXBinit_hterain",nx,ny,1,0,1)
    CALL acct_interrupt(bc_acct)
    CALL bcs2d(nx,ny,hterain, ebc,wbc,nbc,sbc)
    CALL acct_stop_inter
!call test_dump (hterain,"XXXAinit_hterain",nx,ny,1,0,1)

  ELSE IF( ternopt == 2 ) THEN          ! Read from terrain data base

!
!-----------------------------------------------------------------------
!
!    Read in the terrain data.
!
!-----------------------------------------------------------------------
!
!  blocking inserted for ordering i/o for message passing
    DO i=0,nprocs-1,max_fopen
      IF(myproc >= i.AND.myproc <= i+max_fopen-1)THEN

        CALL readtrn( nx,ny,dx,dy, terndta,                             &
               mapproj,trulat1,trulat2,trulon,sclfct,ctrlat,ctrlon,     &
               hterain )

      END IF
      IF (mp_opt > 0) CALL mpbarrier
    END DO

  END IF
!
!-----------------------------------------------------------------------
!
!  Set up a stretched vertical grid.
!
!  For strhopt=1, function y = a+b*x**3 is used to specify dz as a
!                              function of k.
!  For strhopt=2, function y = c + a*tanh(b*x) is used to specify dz
!                              as a function of k.
!
!-----------------------------------------------------------------------
!
  IF ( strhopt == 0 ) THEN

    DO k=1,nz
      zp1d(k) = z(k)
    END DO

  ELSE IF ( strhopt == 1 .OR.strhopt == 2 ) THEN

    z1 = zrefsfc + MAX(0.0, MIN(dlayer1, z(nz-2)-zrefsfc ))
    z2 = z1      + MAX(0.0, MIN(dlayer2, z(nz-1)-z1      ))

    IF( dlayer1 >= (nz-3)*dzmin ) THEN
      WRITE(6,'(/1x,a,f13.3,/a,f13.3,a,a)')                             &
          'Can not setup a vertical grid with uniform dz=',dzmin,       &
          ' over the depth of ',dlayer1,' please specify a smaller ',   &
          'value of dlayer1. Program stopped INIGRD.'
      CALL arpsstop('arpsstop called from INIGRD with ther vertical grid ' &
            ,1)
    END IF

    CALL strhgrd(nz,strhopt,zrefsfc,z1,z2,z(nz-1),                      &
                 dzmin,strhtune, zp1d,dzp1d)

  ELSE

    WRITE(6,'(1x,a,i3,a/)')                                             &
        'Invalid vertical grid stretching option, strhopt was ',strhopt, &
        '. Program stopped in INIGRD.'
      CALL arpsstop('arpsstop called from INIGRD with stretching ' ,1)

  END IF
!
!-----------------------------------------------------------------------
!
!  Physical height of computational grid defined as
!
!  Zp=(z-zrefsfc)*(Zm-hterain)/(Zm-zrefsfc)+hterain for z=<Zm.
!  ZP=z for z>Zm
!
!  where Zm the height at which the grid levels becomes flat.
!  Hm < Zm =< Ztop, hm is the height of mountain and Ztop the height
!  of model top.
!
!-----------------------------------------------------------------------
!
  DO k=nz-1,2,-1
    IF(zp1d(k) <= zflat) THEN
      zflat1 = zp1d(k)
      EXIT
    END IF
  END DO
  zflat1=MAX(MIN(
  DO k=2,nz-1

    IF(zp1d(k) > zflat1) THEN
      DO j=1,ny-1
        DO i=1,nx-1
          zp(i,j,k)=zp1d(k)
        END DO
      END DO
    ELSE
      DO j=1,ny-1
        DO i=1,nx-1
          zp(i,j,k)=(zp1d(k)-zrefsfc)*(zflat1-hterain(i,j))             &
                     /(zflat1-zrefsfc)+hterain(i,j)
        END DO
      END DO
    END IF

  END DO

  DO j=1,ny-1
    DO i=1,nx-1
      zp(i,j,2)=hterain(i,j)
      zp(i,j,1)=2.0*zp(i,j,2)-zp(i,j,3)
      zp(i,j,nz)=2.0*zp(i,j,nz-1)-zp(i,j,nz-2)
    END DO
  END DO
!
!-----------------------------------------------------------------------
!
!  Calculate transformation Jacobians J1, J2 and J3.
!
!-----------------------------------------------------------------------
!
  CALL jacob(nx,ny,nz,x,y,z,zp,j1,j2,j3)

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


SUBROUTINE strhgrd(nz,strhopt,z0,z1,z2,ztop,dzmin,strhtune, z,dzk) 3,1
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  To construct a vertically stretched grid.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Ming Xue
!  10/17/94.
!
!  MODIFICATION HISTORY:
!
!  05/11/95 (Jinxing Zong and MX)
!
!  A bug fix for the case of nonzero zrefsfc, the reference height of
!  the surface. Results not affected for zrefsfc=0 (default value).
!
!-----------------------------------------------------------------------
!
!
!  INPUT:
!
!    nz       The vertical dimension of ARPS grid.
!
!
!  OUTPUT:
!
!    z        Array containing the height of veritical grid levels.
!    dzk      Array containing the grid spacing between vertical levels
!
!-----------------------------------------------------------------------

  IMPLICIT NONE

  INTEGER :: nz
  INTEGER :: strhopt
  REAL :: z0
  REAL :: z1
  REAL :: z2
  REAL :: ztop
  REAL :: dzmin
  REAL :: strhtune

  REAL :: z   (nz)
  REAL :: dzk (nz)

  REAL :: rnzh,dzm
  REAL :: a,b,c,hnew,zkmid,dzu
  INTEGER :: nzh,nzl,k
  REAL :: dz

  IF( (z1-z0) == (nz-3)*dzmin.AND.(z1-z0) == (ztop-z0) ) THEN

    dz = (ztop-z0)/(nz-3)
    DO k=1,nz-1
      dzk(k)= dz
    END DO
    DO k=1,nz
      z(k)=z0 + (k-2) * dz
    END DO

    WRITE(6,'(/1x,a,f13.3,/a,f13.3)')                                   &
        'Layer 1 depth was as deep as the entire domain. i',            &
        'A uniform vertical grid is assumed with dz=',dz,               &
        ' over the model depth of ',ztop-z0

    RETURN

  END IF

  IF(z1 < z0) z1 = z0
  IF(z2 > ztop) z2 = ztop

  nzl = (z1-z0)/dzmin

  IF( (z1-z0) >= (nz-3)*dzmin ) THEN
    WRITE(6,'(/1x,a,f13.3,/a,f13.3,a,a)')                               &
        'Can not setup a vertical grid with uniform dz=',dzmin,         &
        ' over the depth of ',z1-z0,' please specify a smaller',        &
        ' value of dlayer1 '
      CALL arpsstop('arpsstop called from STRHGRD with stretching ' ,1)
  END IF

  IF( z2 >= ztop ) THEN
    dzm = (ztop-z0-nzl*dzmin)/(nz-3-nzl)
!    print*, nzl*dzmin + (nz-3-nzl)*dzm
    nzh = 0
    dzu = 2*dzm - dzmin
  ELSE
    a = 2*(nz-3-nzl)
    b = 2*z0-ztop-z2-(nz-3-3*nzl)*dzmin
    c = dzmin*(z2-z0-nzl*dzmin)
    dzm = (-b + SQRT(b*b-4*a*c) )/(2*a)

    rnzh = (ztop-z2)/(2*dzm-dzmin)
    nzh = INT(rnzh)

    hnew = nzl*dzmin + nzh*(2*dzm-dzmin) +                              &
          (nz-3-nzl-nzh)*dzm + z0

    IF( nzh /= 0 ) THEN
      dzu = (2*dzm-dzmin) + (ztop-hnew)/nzh
    ELSE IF( nz-3-nzl-nzh /= 0 ) THEN
      dzm = dzm + (ztop-hnew)/(nz-3-nzl-nzh)
      dzu = (2*dzm-dzmin)
    END IF

  END IF

  DO k=1,nzl+1
    dzk(k)=dzmin
  END DO


  IF( strhopt == 1 ) THEN

    a   = (dzm-dzmin)
    DO k=nzl+2,nz-2-nzh
      dzk(k)= dzm+a*                                                    &
          ((2.0*FLOAT(k-nzl-2)/FLOAT(nz-4-nzh-nzl)-1.0) )**3
    END DO


  ELSE

    zkmid=0.5*FLOAT( nz-nzh+nzl)

    IF( nzl+2-zkmid == 0.0 ) THEN
      b = 0.0
    ELSE
      b= strhtune* 2.0/(nzl+2-zkmid)
    END IF

    a=(dzmin-dzm)/TANH( strhtune* 2.0)

    DO k=nzl+2,nz-2-nzh
      dzk(k)=dzm + a*TANH(b*(FLOAT(k)-zkmid))
    END DO

  END IF

  DO k=nz-2-nzh+1, nz-2
    dzk(k)= dzu
  END DO

  dzk(nz-1) = dzk(nz-2)
  dzk(nz  ) = dzk(nz-1)


  z(2) = z0
  DO k=3,nz-1
    z(k) = z(k-1)+dzk(k-1)
  END DO

  z(1) =z(2)-dzk(1)
  z(nz-1)=ztop
  z(nz)=z(nz-1)+dzk(nz-1)

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


SUBROUTINE jacob(nx,ny,nz,x,y,z,zp, j1,j2,j3) 6,10
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Calculate transformation Jacobians J1, J2 and J3.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Ming Xue
!  3/17/1991.
!
!  MODIFICATION HISTORY:
!
!  4/10/93 (M. Xue & Hao Jin)
!  Add the terrain.
!
!  9/2/94 (M. Xue)
!  Loop 710 that resets j2 on north and south boundary to
!  zero gradient deleted. It shouldn't have been there.
!
!  9/10/94 (Weygandt & Y. Lu)
!  Cleaned up documentation.
!
!-----------------------------------------------------------------------
!
!  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
!
!  OUTPUT:
!
!    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 (m)
!    zp       Vertical coordinate of grid points in physical space
!             (m)
!
!    j1       Coordinate transformation Jacobian -d(zp)/dx
!    j2       Coordinate transformation Jacobian -d(zp)/dy
!    j3       Coordinate transformation Jacobian  d(zp)/dz
!
!-----------------------------------------------------------------------
!

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

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

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

  REAL :: j1(nx,ny,nz)    ! Coordinate transformation Jacobian
                          ! -d(zp)/dx.
  REAL :: j2(nx,ny,nz)    ! Coordinate transformation Jacobian
                          ! -d(zp)/dy.
  REAL :: j3(nx,ny,nz)    ! Coordinate transformation Jacobian
                          ! d(zp)/dz.
!
!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------
!
  INTEGER :: i,j,k

  INTEGER :: mptag
  REAL, allocatable :: mp_tem(:)
  INTEGER :: istat
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'bndry.inc'
  INCLUDE 'globcst.inc'
  INCLUDE 'mp.inc'            ! Message passing parameters.
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  IF (mp_opt > 0) THEN
    ALLOCATE (mp_tem(MAX(nx,ny)*nz),stat=istat)
    IF (istat /= 0) THEN
      CALL arpsstop ("arpsstop called from JACOB: ERROR allocating mp_tem" &
           ,1)
    END IF
  END IF
!
!-----------------------------------------------------------------------
!
!  Calculate transformation Jacobian J1 defined as -del(zp)/del(x)
!
!-----------------------------------------------------------------------
!
  DO k=1,nz
    DO j=1,ny-1
      DO i=2,nx-1
        j1(i,j,k) = -2 * (zp(i,j,k)-zp(i-1,j,k)) / (x(i+1)-x(i-1))
      END DO
    END DO
  END DO

!
!-----------------------------------------------------------------------
!
!  X - boundary conditions of j1
!
!-----------------------------------------------------------------------
!
!call test_dump (zp,"XXXinit_zp",nx,ny,nz,0,1)
!call test_dump (x,"XXXinit_x",nx,1,1,1,1)
!call test_dump (j1,"XXXinit_j1",nx,ny,nz,1,0)
  IF (mp_opt > 0) THEN
    CALL acct_interrupt(mp_acct)
    CALL mpsend2dew(j1,nx,ny,nz,ebc,wbc,1,mptag,mp_tem)
    CALL mprecv2dew(j1,nx,ny,nz,ebc,wbc,1,mptag,mp_tem)
  END IF

!call test_dump (j1,"XXXBinit_j1",nx,ny,nz,1,1)
  CALL acct_interrupt(bc_acct)

  IF(wbc == 1) THEN            ! Rigid wall boundary condition

    DO k=1,nz
      DO j=1,ny-1
        j1( 1,j,k)=j1( 3 ,j,k)
      END DO
    END DO

  ELSE IF(wbc == 2) THEN        ! Periodic boundary condition.

    DO k=1,nz
      DO j=1,ny-1
        j1( 1,j,k)=j1(nx-2,j,k)
      END DO
    END DO

  ELSE IF(wbc /= 0) THEN

    DO k=1,nz
      DO j=1,ny-1
        j1( 1,j,k)=j1( 2 ,j,k)
      END DO
    END DO

  END IF

  IF(ebc == 1) THEN             ! Rigid wall boundary condition

    DO k=1,nz
      DO j=1,ny-1
        j1(nx,j,k)=j1(nx-2,j,k)
      END DO
    END DO

  ELSE IF(ebc == 2) THEN         ! Periodic boundary condition.

    DO k=1,nz
      DO j=1,ny-1
        j1(nx,j,k)=j1( 3 ,j,k)
      END DO
    END DO

  ELSE IF(ebc /= 0) THEN

    DO k=1,nz
      DO j=1,ny-1
        j1(nx,j,k)=j1(nx-1,j,k)
      END DO
    END DO

  END IF
!call test_dump (j1,"XXXCinit_j1",nx,ny,nz,1,1)

  DO k=1,nz
    DO i=1,nx
      j1(i,ny,k) = j1(i,ny-1,k)
    END DO
  END DO

!call test_dump (j1,"XXXAinit_j1",nx,ny,nz,1,1)
  CALL acct_stop_inter
!
!-----------------------------------------------------------------------
!
!  Calculate transformation Jacobian J2 defined as -del(zp)/del(y)
!
!-----------------------------------------------------------------------
!
  DO k=1,nz
    DO i=1,nx-1
      DO j=2,ny-1
        j2(i,j,k) = -2 * (zp(i,j,k)-zp(i,j-1,k)) / (y(j+1)-y(j-1))
      END DO
    END DO
  END DO
!
!-----------------------------------------------------------------------
!
!  Y - boundary conditions of j2
!
!-----------------------------------------------------------------------
!
!call test_dump (j2,"XXXinit_j2",nx,ny,nz,2,0)
  IF (mp_opt > 0) THEN
    CALL acct_interrupt(mp_acct)
    CALL mpsend2dns(j2,nx,ny,nz,nbc,sbc,2,mptag,mp_tem)
    CALL mprecv2dns(j2,nx,ny,nz,nbc,sbc,2,mptag,mp_tem)
  END IF

  CALL acct_interrupt(bc_acct)

  IF(sbc == 1) THEN            ! Rigid wall boundary condition

    DO k=1,nz
      DO i=1,nx-1
        j2(i, 1,k)=j2(i, 3 ,k)
      END DO
    END DO

  ELSE IF(sbc == 2) THEN        ! Periodic boundary condition.

    IF (mp_opt == 0) THEN
      DO k=1,nz
      DO i=1,nx-1
        j2(i, 1,k)=j2(i,ny-2,k)
      END DO
      END DO
    END IF

  ELSE IF(sbc /= 0) THEN

    DO k=1,nz
      DO i=1,nx-1
        j2(i, 1,k)=j2(i, 2 ,k)
      END DO
    END DO

  END IF

  IF(nbc == 1) THEN             ! Rigid wall boundary condition

    DO k=1,nz
      DO i=1,nx-1
        j2(i,ny,k)=j2(i,ny-2,k)
      END DO
    END DO

  ELSE IF(nbc == 2) THEN         ! Periodic boundary condition.

    IF (mp_opt == 0) THEN
      DO k=1,nz
      DO i=1,nx-1
        j2(i,ny,k)=j2(i, 3 ,k)
      END DO
      END DO
    END IF

  ELSE IF(nbc /= 0) THEN

    DO k=1,nz
      DO i=1,nx-1
        j2(i,ny,k)=j2(i,ny-1,k)
      END DO
    END DO

  END IF

  DO k=1,nz
    DO j=1,ny
      j2(nx,j,k) = j2(nx-1,j,k)
    END DO
  END DO

!call test_dump (j2,"XXXAinit_j2",nx,ny,nz,2,1)
  CALL acct_stop_inter
!
!-----------------------------------------------------------------------
!
!  Calculate transformation Jacobian J3 defined as del(zp)/del(z)
!
!-----------------------------------------------------------------------
!
  DO k=1,nz-1
    DO j=1,ny-1
      DO i=1,nx-1
        j3(i,j,k) = (zp(i,j,k+1)-zp(i,j,k))/(z(k+1)-z(k))
      END DO
    END DO
  END DO

  DO k=1,nz-1
    DO j=1,ny-1
      j3(nx,j,k) = j3(nx-1,j,k)
    END DO
  END DO

  DO k=1,nz-1
    DO i=1,nx
      j3(i,ny,k) = j3(i,ny-1,k)
    END DO
  END DO

  DO j=1,ny
    DO i=1,nx
      j3(i,j,nz) = j3(i,j,nz-1)
    END DO
  END DO

  IF (mp_opt > 0) THEN
    DEALLOCATE (mp_tem,stat=istat)
  END IF

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

!EMK BMJ

SUBROUTINE initgrdvar(nx,ny,nz,nts,nstyps,exbcbufsz,                    & 2,46
           x,y,z,zp,hterain,mapfct,                                     &
           j1,j2,j3,aj3x,aj3y,aj3z,j3inv,                               &
           u,v,w,wcont,ptprt,pprt,qv,qc,qr,qi,qs,qh,tke,                &
           udteb, udtwb, vdtnb, vdtsb,                                  &
           pdteb,pdtwb ,pdtnb ,pdtsb,                                   &
           trigs1,trigs2,ifax1,ifax2,                                   &
           wsave1,wsave2,vwork1,vwork2,                                 &
           ubar,vbar,ptbar,pbar,ptbari,pbari,                           &
           rhostr,rhostri,qvbar,ppi,csndsq,                             &
           soiltyp,stypfrct,vegtyp,lai,roufns,veg,                      &
           tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth,qvsfc,              &
           ptcumsrc,qcumsrc,w0avg,nca,kfraincv,                         &
           cldefi,xland,bmjraincv,                                      &
           raing,rainc,prcrate,exbcbuf,                                 &
           radfrc,radsw,rnflx,                                          &
           usflx,vsflx,ptsflx,qvsflx,                                   &
           tem1,tem2,tem3,tem4,tem5,                                    &
           tem6,tem7,tem8,tem9)
!EMK END
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Initialize the model array variables.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Ming Xue
!  3/17/1992.
!
!  MODIFICATION HISTORY:
!
!  5/02/92 (M. Xue)
!  Added full documentation.
!
!  5/03/92 (M. Xue)
!  Further documentation.
!
!  10/7/1992 (M. Xue)
!  Added call to subroutine extinit, the option three of
!  initialization.
!
!  9/10/94 (Weygandt & Y. Lu)
!  Cleaned up documentation.
!
!  02/07/1995 (Yuhe Liu)
!  Added a new 2-D permanent array, veg(nx,ny), to the argument list
!
!  10/31/95  (D. Weber)
!  Added trigs1,trigs2,ifax1,ifax2 for use in the fft code for the
!  upper radiation condition.
!
!  1/22/96 (Donghai Wang & Yuhe Liu)
!  Added the map projection factor to ARPS governing equations.
!
!  07/22/97 (D. Weber)
!  Added wsave1,wsave2,vwork1,vwork2 for use in the even fft version
!  of the upper w-p radiation condition (fftopt=2).
!
!  08/01/97 (Zonghui Huo)
!  Added Kain-fritsch cumulus parameterization scheme.
!
!  11/06/97 (D. Weber)
!  Added three additional levels to the mapfct array.  The three
!  levels (4,5,6) represent the inverse of the first three in order.
!  The inverse map factors are computed to improve efficiency.
!
!  12/05/97 (K. Brewster)
!  Added argument, nt, so that routines that do not require more
!  than one time level can be initialized using less memory.
!
!  4/15/1998 (Donghai Wang)
!  Added the source terms to the right hand terms of the qc,qr,qi,qs
!  equations due to the K-F cumulus parameterization.
!
!  4/15/1998 (Donghai Wang)
!  Added the running average vertical velocity (array w0avg)
!  for the K-F cumulus parameterization scheme.
!
!  11/18/98 (Keith Brewster)
!  Changed pibar to ppi (full pi).
!
!  12/09/1998 (Donghai Wang)
!  Added the snow cover.
!
!  11/06/2001 (Yunheng Wang)
!  Added mpupdatei calling for ict/icb to solve radiation forcing
!  difference for MPI run.
!
!  2002/02/28 (Gene Bassett)
!  Replaced mpupdatei for ict/icb to a call to mpmax0.
!
!  13 March 2002 (Eric Kemp)
!  Added arrays for WRF BMJ cumulus scheme.
!
!  04/10/2002 (Yunheng Wang)
!  Subsituted mpupdatei calls for mpmax0 calls again
!  because mpmax0 calls were not correct.
!
!-----------------------------------------------------------------------
!
!
!  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
!    nts      Number of time levels to be initialized.
!
!    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 (m)
!    zp       Vertical coordinate of grid points in physical space
!             (m)
!    hterain  The height of terrain (m)
!
!    mapfct   Map factors at scalar, u and v points
!
!    j1       Coordinate transformation Jacobian -d(zp)/dx
!    j2       Coordinate transformation Jacobian -d(zp)/dy
!    j3       Coordinate transformation Jacobian  d(zp)/dz
!    aj3x     Avgx of the coordinate transformation Jacobian  d(zp)/dz
!    aj3y     Avgy of the coordinate transformation Jacobian  d(zp)/dz
!    aj3z     Avgz of the coordinate transformation Jacobian  d(zp)/dz
!    trigs1   Array containing pre-computed trig function for fftopt=1.
!    trigs1   Array containing pre-computed trig function for fftopt=1.
!    ifax1    Array containing the factors of nx for fftopt=1.
!    ifax2    Array containing the factors of ny for fftopt=1.
!
!    vwork1   2-D work array for fftopt=2.
!    vwork2   2-D work array for fftopt=2.
!    wsave1   Work array for fftopt=2.
!    wsave2   Work array for fftopt=2.
!
!  OUTPUT:
!
!    u        x component of velocity at all time levels (m/s)
!    v        y component of velocity at all time levels (m/s)
!    w        Vertical component of Cartesian velocity
!             at all time levels (m/s)
!    wcont    Contravariant vertical velocity (m/s)
!    ptprt    Perturbation potential temperature at all time levels
!             (K)
!    pprt     Perturbation pressure at all time levels (Pascal)
!    qv       Water vapor specific humidity at all time levels
!             (kg/kg)
!    qc       Cloud water mixing ratio at all time levels (kg/kg)
!    qr       Rainwater mixing ratio at all time levels (kg/kg)
!    qi       Cloud ice mixing ratio at all time levels (kg/kg)
!    qs       Snow mixing ratio at all time levels (kg/kg)
!    qh       Hail mixing ratio at all time levels (kg/kg)
!    tke      Turbulent Kinetic Energy ((m/s)**2)
!
!    udteb    Time tendency of u field at east boundary (m/s**2)
!    udtwb    Time tendency of u field at west boundary (m/s**2)
!
!    vdtnb    Time tendency of v field at north boundary (m/s**2)
!    vdtsb    Time tendency of v field at south boundary (m/s**2)
!
!    pdteb    Time tendency of pprt field at east boundary (PASCAL/s)
!    pdtwb    Time tendency of pprt field at west boundary (PASCAL/s)
!    pdtnb    Time tendency of pprt field at north boundary
!             (PASCAL/s)
!    pdtsb    Time tendency of pprt field at south boundary
!             (PASCAL/s)
!
!    ubar     Base state zonal velocity component (m/s)
!    vbar     Base state meridional velocity component (m/s)
!    ptbar    Base state potential temperature (K)
!    pbar     Base state pressure (Pascal)
!    ptbari   Inverse Base state potential temperature (K)
!    pbari    Inverse Base state pressure (Pascal)
!    rhostr   Base state density rhobar times j3 (kg/m**3)
!    rhostri  Inverse base state density rhobar times j3 (kg/m**3)
!    qvbar    Base state water vapor specific humidity (kg/kg)
!    ppi      Exner function.
!    csndsq   Sound wave speed squared.
!
!    soiltyp  Soil type
!    vegtyp   Vegetation type
!    lai      Leaf Area Index
!    roufns   Surface roughness
!    veg      Vegetation fraction
!
!    tsfc     Skin temperature at the ground or ocean surface (K)
!    tsoil    Deep soil temperature (K) (in deep 1 m layer)
!    wetsfc   Surface soil moisture
!    wetdp    Deep soil moisture
!    wetcanp  Canopy water amount
!    snowdpth Snow depth (m)
!    qvsfc    Effective S.H. at sfc.
!
!    ptcumsrc Source term in pt-equation due to cumulus parameterization
!    qcumsrc Source term in water equations due to cumulus parameterization
!    kfraincv K-F convective rainfall (cm)
!    nca      K-F counter for CAPE release
!    cldefi   BMJ cloud efficiency
!    xland    BMJ land/sea mask
!    bmjraincv BMJ convective rainfall (cm)
!
!    radfrc   Radiation forcing (K)
!    radsw    Solar radiation reaching the surface
!    rnflx    Net absorbed radiation by the surface
!
!    raing    Grid supersaturation rain
!    rainc    Cumulus convective rain
!
!  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.
!
!-----------------------------------------------------------------------
!

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

  INTEGER :: nts               ! Number of time levels to be initialized.
  INTEGER :: tpast             ! Index of time level for the past time.
  INTEGER :: tpresent          ! Index of time level for the present
                               ! time.
  INTEGER :: tfuture           ! Index of time level for the future
                               ! time.
  INTEGER :: nx,ny,nz          ! The number of grid points in 3
                               ! directions

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

  REAL :: mapfct(nx,ny,8)      ! Map factors at scalar, u and v points

  REAL :: j1    (nx,ny,nz)     ! Coordinate transformation Jacobian
                               ! defined as - d( zp )/d( x ).
  REAL :: j2    (nx,ny,nz)     ! Coordinate transformation Jacobian
                               ! defined as - d( zp )/d( y ).
  REAL :: j3    (nx,ny,nz)     ! Coordinate transformation Jacobian
                               ! defined as d( zp )/d( z ).
  REAL :: aj3x  (nx,ny,nz)     ! Coordinate transformation Jacobian defined
                               ! as d( zp )/d( z ) AVERAGED IN THE X-DIR.
  REAL :: aj3y  (nx,ny,nz)     ! Coordinate transformation Jacobian defined
                               ! as d( zp )/d( z ) AVERAGED IN THE Y-DIR.
  REAL :: aj3z  (nx,ny,nz)     ! Coordinate transformation Jacobian defined
                               ! as d( zp )/d( z ) AVERAGED IN THE Z-DIR.
  REAL :: j3inv (nx,ny,nz)     ! Inverse of j3

  REAL :: trigs1(3*(nx-1)/2+1) ! Array containing pre-computed trig
                               ! function for fftopt=1.
  REAL :: trigs2(3*(ny-1)/2+1) ! Array containing pre-computed trig
                               ! function for fftopt=1.
  INTEGER :: ifax1(13)         ! Array containing the factors of nx for
                               ! fftopt=1.
  INTEGER :: ifax2(13)         ! Array containing the factors of ny for
                               ! fftopt=1.

  REAL :: vwork1 (nx+1,ny+1)   ! 2-D work array for fftopt=2.
  REAL :: vwork2 (ny,nx+1)     ! 2-D work array for fftopt=2.
  REAL :: wsave1 (3*(ny-1)+15) ! Work array for fftopt=2.
  REAL :: wsave2 (3*(nx-1)+15) ! Work array for fftopt=2.

  REAL :: u     (nx,ny,nz,nts) ! Total u-velocity (m/s)
  REAL :: v     (nx,ny,nz,nts) ! Total v-velocity (m/s)
  REAL :: w     (nx,ny,nz,nts) ! Total w-velocity (m/s)
  REAL :: wcont (nx,ny,nz)     ! Contravariant vertical velocity (m/s)
  REAL :: ptprt (nx,ny,nz,nts) ! Perturbation potential temperature (K)
  REAL :: pprt  (nx,ny,nz,nts) ! Perturbation pressure (Pascal)

  REAL :: qv    (nx,ny,nz,nts) ! Water vapor specific humidity (kg/kg)
  REAL :: qc    (nx,ny,nz,nts) ! Cloud water mixing ratio (kg/kg)
  REAL :: qr    (nx,ny,nz,nts) ! Rain water mixing ratio (kg/kg)
  REAL :: qi    (nx,ny,nz,nts) ! Cloud ice mixing ratio (kg/kg)
  REAL :: qs    (nx,ny,nz,nts) ! Snow mixing ratio (kg/kg)
  REAL :: qh    (nx,ny,nz,nts) ! Hail mixing ratio (kg/kg)
  REAL :: tke   (nx,ny,nz,nts) ! Turbulent Kinetic Energy ((m/s)**2)

  REAL :: udteb (ny,nz)        ! T-tendency of u at e-boundary (m/s**2)
  REAL :: udtwb (ny,nz)        ! T-tendency of u at w-boundary (m/s**2)

  REAL :: vdtnb (nx,nz)        ! T-tendency of v at n-boundary (m/s**2)
  REAL :: vdtsb (nx,nz)        ! T-tendency of v at s-boundary (m/s**2)

  REAL :: pdteb (ny,nz)        ! T-tendency of pprt at e-boundary
                               ! (PASCAL/s)
  REAL :: pdtwb (ny,nz)        ! T-tendency of pprt at w-boundary
                               ! (PASCAL/s)
  REAL :: pdtnb (nx,nz)        ! T-tendency of pprt at n-boundary
                               ! (PASCAL/s)
  REAL :: pdtsb (nx,nz)        ! T-tendency of pprt at s-boundary
                               ! (PASCAL/s)

  REAL :: ubar  (nx,ny,nz)     ! Base state u-velocity (m/s)
  REAL :: vbar  (nx,ny,nz)     ! Base state v-velocity (m/s)
  REAL :: ptbar (nx,ny,nz)     ! Base state potential temperature (K)
  REAL :: pbar  (nx,ny,nz)     ! Base state pressure (Pascal).
  REAL :: ptbari(nx,ny,nz)     ! Inverse Base state pot. temperature (K)
  REAL :: pbari (nx,ny,nz)     ! Inverse Base state pressure (Pascal).
  REAL :: rhostr(nx,ny,nz)     ! Base state density rhobar times j3.
  REAL :: rhostri(nx,ny,nz)    ! Inv. base state density rhobar times j3.
  REAL :: qvbar (nx,ny,nz)     ! Base state water vapor specific
                               ! humidity(kg/kg)
  REAL :: ppi   (nx,ny,nz)     ! Exner function.
  REAL :: csndsq(nx,ny,nz)     ! Sound wave speed squared.

  INTEGER :: nstyps                  ! Number of soil types
  INTEGER :: soiltyp(nx,ny,nstyps)   ! Soil types at grid points
  REAL :: stypfrct(nx,ny,nstyps)  ! Fraction of soil types
  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)    ! Temperature at ground (K)
                                     ! (in top 1 cm layer)
  REAL :: qvsfc  (nx,ny,0:nstyps)    ! Effective qv at sfc.
  REAL :: tsoil  (nx,ny,0:nstyps)    ! Deep soil temperature (K)
                                     ! (in deep 1 m layer)
  REAL :: wetsfc (nx,ny,0:nstyps)    ! Surface soil moisture in the top
                                     ! 1 cm layer
  REAL :: wetdp  (nx,ny,0:nstyps)    ! Deep soil moisture in the deep
                                     ! 1 m layer
  REAL :: wetcanp(nx,ny,0:nstyps)    ! Canopy water amount
  REAL :: snowdpth(nx,ny)            ! Snow depth (m)

  REAL :: ptcumsrc(nx,ny,nz)   ! Source term in pt-equation due
                               ! to cumulus parameterization
  REAL :: qcumsrc(nx,ny,nz,5)  ! Source term in water equations due
                               ! to cumulus parameterization:
                               ! qcumsrc(1,1,1,1) for qv equation
                               ! qcumsrc(1,1,1,2) for qc equation
                               ! qcumsrc(1,1,1,3) for qr equation
                               ! qcumsrc(1,1,1,4) for qi equation
                               ! qcumsrc(1,1,1,5) for qs equation
  REAL :: w0avg(nx,ny,nz)      ! a closing running average vertical
                               ! velocity in 10min for K-F scheme
  REAL :: kfraincv(nx,ny)      ! K-F convective rainfall (cm)
  INTEGER :: nca(nx,ny)        ! K-F counter for CAPE release

!EMK BMJ
  REAL,INTENT(INOUT) :: cldefi(nx,ny)      ! BMJ cloud efficiency
  REAL,INTENT(INOUT) :: xland(nx,ny)       ! BMJ land mask
                                           !   (1.0 = land, 2.0 = sea)
  REAL,INTENT(INOUT) :: bmjraincv(nx,ny)   ! BMJ convective rainfall (cm)
!EMK END

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

  REAL :: usflx (nx,ny)        ! Surface flux of u-momentum (kg/(m*s**2))
  REAL :: vsflx (nx,ny)        ! Surface flux of v-momentum (kg/(m*s**2))
  REAL :: ptsflx(nx,ny)        ! Surface heat flux (K*kg/(m*s**2))
  REAL :: qvsflx(nx,ny)        ! Surface moisture flux (kg/(m**2*s))

  REAL :: raing(nx,ny)         ! Grid supersaturation rain
  REAL :: rainc(nx,ny)         ! Cumulus convective rain
  REAL :: prcrate(nx,ny,4)     ! precipitation rate (kg/(m**2*s))
                               ! prcrate(1,1,1) = total precipitation rate
                               ! prcrate(1,1,2) = grid scale precip. rate
                               ! prcrate(1,1,3) = cumulus precip. rate
                               ! prcrate(1,1,4) = microphysics precip. rate

  INTEGER :: exbcbufsz         ! EXBC buffer size
  REAL :: exbcbuf( exbcbufsz ) ! EXBC buffer array

  REAL :: tem1  (nx,ny,nz)     ! Temporary work array.
  REAL :: tem2  (nx,ny,nz)     ! Temporary work array.
  REAL :: tem3  (nx,ny,nz)     ! Temporary work array.
  REAL :: tem4  (nx,ny,nz)     ! Temporary work array.
  REAL :: tem5  (nx,ny,nz)     ! Temporary work array.
  REAL :: tem6  (nx,ny,nz)     ! Temporary work array.
  REAL :: tem7  (nx,ny,nz)     ! Temporary work array.
  REAL :: tem8  (nx,ny,nz)     ! Temporary work array.
  REAL :: tem9  (nx,ny,nz)     ! Temporary work array.
!
!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------
!
  INTEGER :: i,j,k
  REAL :: temp

  REAL :: alatpro(2)
  REAL :: sclf
  REAL :: dxscl             ! Model x-direction grid spacing
                            ! normalized by the map scale
                            ! dxscl=dx/sclf
  REAL :: dyscl             ! Model y-direction grid spacing
                            ! normalized by the map scale
                            ! dyscl=dy/sclf
  REAL :: xs,ys, swx,swy, ctrx, ctry
  REAL zpmax

  INTEGER :: mptag
!wdt update
  REAL :: rmin,rmax
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'globcst.inc'
  INCLUDE 'grid.inc'          ! Grid & map parameters.
  INCLUDE 'bndry.inc'
  INCLUDE 'phycst.inc'
  INCLUDE 'mp.inc'            ! Message passing parameters.
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
!-----------------------------------------------------------------------
!
!  If initopt = 1, initialize the model fields using intialization
!                  routines. Typically from an analytical
!                  definition of initial perturbations.
!
!  If initopt = 2, initialize the model fields from a restart file
!                  produced by a previous model run.
!
!  If initopt = 3, initialize the model fields by reading in an
!                  external data file.
!
!-----------------------------------------------------------------------
!
  IF( nts > 1 ) THEN
    tpast=1
    tpresent=2
    tfuture=3
  ELSE
    tpast=1
    tpresent=1
    tfuture=1
  END IF

  IF( initopt == 1 ) THEN
!
!-----------------------------------------------------------------------
!
!  Initialization of model GRiD definition arrays.
!
!-----------------------------------------------------------------------
!
    CALL inigrd(nx,ny,nz,                                               &
                x,y,z,zp,hterain,mapfct,j1,j2,j3,tem1,tem2,tem3)
!
!-----------------------------------------------------------------------
!
!  Define the base state atmospheric variables.
!
!-----------------------------------------------------------------------
!
    IF ((inibasopt == 1) .AND. (max_fopen < nprocs)) THEN
      CALL wrtcomment("ERROR: for message passing version with "//      &
                    "inibasopt=1, max_fopen (in arps.input)",1)
      CALL arpsstop ("arpsstop called from initgrdvar due to nproc_x-y  &
                    & nproc_x*nproc_y (in arps.input).",1)
    END IF

!  blocking inserted for ordering i/o for message passing
    DO i=0,nprocs-1,max_fopen
      IF(myproc >= i.AND.myproc <= i+max_fopen-1)THEN
        CALL inibase(nx,ny,nz, ubar,vbar,ptbar,pbar,ptbari,pbari,       &
                   rhostr,rhostri,qvbar,                                &
                   x,y,z,zp,j3, tem1,tem2,tem3,tem4)
      END IF
      IF (mp_opt > 0) CALL mpbarrier
    END DO
!
!-----------------------------------------------------------------------
!
!  Initialize time dependent model variables.
!
!-----------------------------------------------------------------------
!
    CALL initdvr(nx,ny,nz,nts,                                          &
                 ubar,vbar,ptbar,pbar,rhostr,qvbar,                     &
                 x,y,z,zp,hterain, j1,j2,j3,                            &
                 u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh,                    &
                 ptcumsrc,qcumsrc,raing,rainc,prcrate,tem1)
!
!-----------------------------------------------------------------------
!
!  Set the time tendencies of u, v and pprt on the lateral boundaries
!  to zero for the initial time step
!
!  These tendencies will be used by lateral boundary condition options
!  4 and 5 only.
!
!-----------------------------------------------------------------------
!
    DO k=1,nz
      DO j=1,ny
        udteb(j,k) = 0.0
        udtwb(j,k) = 0.0
        pdteb(j,k) = 0.0
        pdtwb(j,k) = 0.0
      END DO
    END DO

    DO k=1,nz
      DO i=1,nx
        vdtnb(i,k) = 0.0
        vdtsb(i,k) = 0.0
        pdtnb(i,k) = 0.0
        pdtsb(i,k) = 0.0
      END DO
    END DO

  ELSE IF( initopt == 2) THEN   ! restart run
!
!-----------------------------------------------------------------------
!
!  Read in restart data from a restart file to initialize fields
!  u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh at time tpast and tpresent,
!  the base state variables ubar,vbar,ptbar,pbar,rhostr,qvbar,
!  and the time tendencies of variables at the lateral boundries.
!
!  Fields at tfuture are set to equal to those at tpresent.
!
!  This subroutine also sets the value of tstart to the restart
!  data time. The value from input file in over-written.
!
!-----------------------------------------------------------------------
!
!EMK BMJ
    CALL rstin(nx,ny,nz,nts, nstyps, exbcbufsz,                         &
               u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh,tke,                  &
               udteb, udtwb, vdtnb, vdtsb,                              &
               pdteb, pdtwb, pdtnb, pdtsb,                              &
               ubar,vbar,ptbar,pbar,rhostr,qvbar,                       &
               x,y,z,zp,hterain,mapfct,j1,j2,j3,                        &
               soiltyp,stypfrct,vegtyp,lai,roufns,veg,                  &
               tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth,qvsfc,          &
               ptcumsrc,qcumsrc,w0avg,nca,kfraincv,                     &
               cldefi,xland,bmjraincv,                                  &
               radfrc,radsw,rnflx,                                      &
               raing,rainc,prcrate,exbcbuf,tem1,tem2)
!EMK END
!
!-----------------------------------------------------------------------
!
!  Set map projection parameters which were not stored in restart
!  data file.
!
!-----------------------------------------------------------------------
!
    alatpro(1) = trulat1
    alatpro(2) = trulat2

    IF( sclfct /= 1.0) THEN
      sclf  = 1.0/sclfct
      dxscl = dx*sclf
      dyscl = dy*sclf
    ELSE
      sclf  = 1.0
      dxscl = dx
      dyscl = dy
    END IF

    CALL setmapr( mapproj,sclf,alatpro,trulon )

    CALL lltoxy( 1,1, ctrlat,ctrlon, ctrx, ctry )

!       swx = ctrx - (float(nx-3)/2.) * dxscl
!       swy = ctry - (float(ny-3)/2.) * dyscl
    swx = ctrx - (FLOAT(nproc_x*(nx-3))/2.) * dxscl
    swy = ctry - (FLOAT(nproc_y*(ny-3))/2.) * dyscl

    CALL setorig( 1, swx, swy)
                               ! set up the model origin to the coord.

    CALL setcornerll(nx,ny, x, y)

  ELSE IF( initopt == 3) THEN  ! External data input.
!
!-----------------------------------------------------------------------
!
!  Read in externally supplied initial fields.  These fields include
!  u, v, w, ptprt, pprt, qv, qc, qr, qi, qs, and qh at time level
!  tpresent, and the base state variables ubar, vbar, ptbar, pbar,
!  rhostr,qvbar.
!
!  Fields at tpast and tfuture are set to their values at tpresent.
!
!-----------------------------------------------------------------------
!
    CALL extinit(nx,ny,nz,nts,nstyps,                                   &
                 u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh,tke,                &
                 ubar,vbar,ptbar,pbar,rhostr,qvbar,                     &
                 x,y,z,zp,hterain,j1,j2,j3,                             &
                 soiltyp,stypfrct,vegtyp,lai,roufns,veg,                &
                 tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth,qvsfc,        &
                 ptcumsrc,qcumsrc,raing,rainc,prcrate,                  &
                 radfrc,radsw,rnflx,                                    &
                 usflx,vsflx,ptsflx,qvsflx,                             &
                 tem1,tem2,tem3,tem4,tem5,tem6,tem7,tem8,tem9)
!
!-----------------------------------------------------------------------
!
!  Set map projection parameters which were not stored in history
!  data file.
!
!-----------------------------------------------------------------------
!
    alatpro(1) = trulat1
    alatpro(2) = trulat2

    IF( sclfct /= 1.0) THEN
      sclf  = 1.0/sclfct
      dxscl = dx*sclf
      dyscl = dy*sclf
    ELSE
      sclf  = 1.0
      dxscl = dx
      dyscl = dy
    END IF

    CALL setmapr( mapproj,sclf,alatpro,trulon )

    CALL lltoxy( 1,1, ctrlat,ctrlon, ctrx, ctry )

!       swx = ctrx - (float(nx-3)/2.) * dxscl
!       swy = ctry - (float(ny-3)/2.) * dyscl
    swx = ctrx - (FLOAT(nproc_x*(nx-3))/2.) * dxscl
    swy = ctry - (FLOAT(nproc_y*(ny-3))/2.) * dyscl

    CALL setorig( 1, swx, swy)
                               ! set up the model origin to the coord.

    IF ( mpfctopt /= 0 ) THEN

      DO j=1,ny-1
        DO i=1,nx-1
          xs = 0.5*(x(i)+x(i+1))
          ys = 0.5*(y(j)+y(j+1))
          CALL xytomf( 1,1,xs,ys,mapfct(i,j,1) )
          mapfct(i,j,4) = 1.0/mapfct(i,j,1)
          mapfct(i,j,7) = mapfct(i,j,1)*mapfct(i,j,1)
          mapfct(i,j,8) = 0.25*mapfct(i,j,1)
        END DO
      END DO

      DO j=1,ny-1
        DO i=1,nx
          ys = 0.5*(y(j)+y(j+1))
          CALL xytomf( 1,1,x(i),ys,mapfct(i,j,2) )
          mapfct(i,j,5) = 1.0/mapfct(i,j,2)
        END DO
      END DO

      DO j=1,ny
        DO i=1,nx-1
          xs = 0.5*(x(i)+x(i+1))
          CALL xytomf( 1,1,xs,y(j),mapfct(i,j,3) )
          mapfct(i,j,6) = 1.0/mapfct(i,j,3)
        END DO
      END DO

    ELSE

      DO k=1,7
        DO j=1,ny
          DO i=1,nx
            mapfct(i,j,k) = 1.0
            mapfct(i,j,8) = 0.25
          END DO
        END DO
      END DO

    END IF

    CALL setcornerll(nx,ny, x, y)
!
!-----------------------------------------------------------------------
!
!  Set the time tendencies of u, v and pprt on the lateral boundaries
!  to zero for the initial time step
!
!  These tendencies will be used by lateral boundary condition options
!  4 and 5 only.
!
!-----------------------------------------------------------------------
!
    DO k=1,nz
      DO j=1,ny
        udteb(j,k) = 0.0
        udtwb(j,k) = 0.0
        pdteb(j,k) = 0.0
        pdtwb(j,k) = 0.0
      END DO
    END DO

    DO k=1,nz
      DO i=1,nx
        vdtnb(i,k) = 0.0
        vdtsb(i,k) = 0.0
        pdtnb(i,k) = 0.0
        pdtsb(i,k) = 0.0
      END DO
    END DO

  ELSE

    WRITE(6,'(1x,a,i3,a)')                                              &
         'Invalid option for intialization, INITOPT was ',initopt,      &
         ' Job stopped.'

      CALL arpsstop ("arpsstop called from initgrdvar due to initopt",1)

  END IF
!
!-----------------------------------------------------------------------
!
!  Define the reversed vertical indeces of height cldh2m and cldm2l
!  which represent the levels that separate high, middle, and low
!  clouds in the computation of the solar radiation transfer
!
!-----------------------------------------------------------------------
!
  ict = nz-2
  icb = nz-2

  DO k=nz-2,2,-1
    tem1(1,1,k) = (zp(1,1,k)-zp(1,1,2))*(zflat-zrefsfc)                 &
                 /(zflat-zp(1,1,2))+zrefsfc
  END DO

  DO k=nz-2,2,-1
    IF ( tem1(1,1,k) <= cldh2m) THEN
      ict = k
      EXIT
    END IF
  END DO

! for bit-for-bit MP accuracy:
  CALL mpupdatei(ict, 1)

  DO k=nz-2,2,-1
    IF ( tem1(1,1,k) <= cldm2l) THEN
      icb = k
      EXIT
    END IF
  END DO

! for bit-for-bit MP accuracy:
  CALL mpupdatei(icb, 1)

  DO k=1,nz-1
    DO j=1,ny-1
      DO i=1,nx-1
        j3inv(i,j,k)=1.0/j3(i,j,k)
      END DO
    END DO
  END DO

  DO k=1,nz-1
    DO j=1,ny-1
      DO i=2,nx-1
        aj3x(i,j,k)=0.5*(j3(i,j,k)+j3(i-1,j,k))
      END DO
    END DO
  END DO

!call test_dump (aj3x,"XXXinit_aj3x",nx,ny,nz,1,0)
  IF (mp_opt > 0) THEN
    CALL acct_interrupt(mp_acct)
    CALL mpsend2dew(aj3x,nx,ny,nz,ebc,wbc,1,mptag,tem2)
    CALL mprecv2dew(aj3x,nx,ny,nz,ebc,wbc,1,mptag,tem2)
!    CALL mpsend2dns(aj3x,nx,ny,nz,1,mptag,tem2)
!    CALL mprecv2dns(aj3x,nx,ny,nz,1,mptag,tem2)
  END IF
!call test_dump (aj3x,"XXXBinit_aj3x",nx,ny,nz,1,1)
  CALL acct_interrupt(bc_acct)
  CALL bcsu(nx,ny,nz,1,ny-1,1,nz-1,ebc,wbc,aj3x)
  CALL acct_stop_inter
!call test_dump (aj3x,"XXXAinit_aj3x",nx,ny,nz,1,1)

  DO k=1,nz-1
    DO j=2,ny-1
      DO i=1,nx-1
        aj3y(i,j,k)=0.5*(j3(i,j,k)+j3(i,j-1,k))
      END DO
    END DO
  END DO

!call test_dump (aj3y,"XXXinit_aj3y",nx,ny,nz,2,0)
  IF (mp_opt > 0) THEN
    CALL acct_interrupt(mp_acct)
!    CALL mpsend2dew(aj3y,nx,ny,nz,2,mptag,tem2)
!    CALL mprecv2dew(aj3y,nx,ny,nz,2,mptag,tem2)
    CALL mpsend2dns(aj3y,nx,ny,nz,nbc,sbc,2,mptag,tem2)
    CALL mprecv2dns(aj3y,nx,ny,nz,nbc,sbc,2,mptag,tem2)
  END IF
  CALL acct_interrupt(bc_acct)
  CALL bcsv(nx,ny,nz,1,nx-1,1,nz-1,nbc,sbc,aj3y)
  CALL acct_stop_inter
!call test_dump (aj3y,"XXXAinit_aj3y",nx,ny,nz,2,1)

  DO k=2,nz-1
    DO j=1,ny-1
      DO i=1,nx-1
        aj3z(i,j,k)=0.5*(j3(i,j,k)+j3(i,j,k-1))
      END DO
    END DO
  END DO

  CALL acct_interrupt(bc_acct)
  CALL bcsw(nx,ny,nz,1,nx-1,1,ny-1,tbc,bbc,aj3z)
  CALL acct_stop_inter
!call test_dump (aj3z,"XXXAinit_aj3z",nx,ny,nz,3,1)
!
!-----------------------------------------------------------------------
!
!  Calculate the trigs1,trigs2,ifax1,ifax2 arrays by calling set99
!
!                              OR
!
!  calculate wsave1,wsave2 by calling vcosi.
!
!-----------------------------------------------------------------------
!
  DO i=1,13
    ifax1(i) = 0
    ifax2(i) = 0
  END DO

  DO i=1,3*(nx-1)/2+1
    trigs1(i) = 0
  END DO

  DO j=1,3*(ny-1)/2+1
    trigs2(j) = 0
  END DO

  DO j=1,ny+1
    DO i=1,nx+1
      vwork1(i,j) = 0.0
    END DO
  END DO

  DO i=1,nx+1
    DO j=1,ny
      vwork2(j,i) = 0.0
    END DO
  END DO

  DO i=1,3*(nx-1)+15
    wsave2(i) = 0.0
  END DO

  DO i=1,3*(ny-1)+15
    wsave1(i) = 0.0
  END DO

  IF ( tbc == 4 ) THEN  ! set up the fft work arrays for use in the
                        ! upper radiation boundary condition.
    IF ( fftopt == 1 ) THEN     ! set up periodic work arrays...
      IF ( ny == 4 ) THEN       ! set up trigs in x direction only
        CALL set99(trigs1,ifax1,nx-1)    ! NOTE: nx must be ODD!!!!
                                         ! and of special character...
                                         ! see fft99f.f for details....
      ELSE IF ( nx == 4 ) THEN     ! set up trigs in y direction only

        CALL set99(trigs2,ifax2,ny-1)    ! NOTE: ny must be ODD!!!!
                                         ! and of special character...
                                         ! see fft99f.f for details....
      ELSE    ! set up for 2-d transform...

        CALL set99(trigs1,ifax1,nx-1)    ! NOTE: nx must be ODD!!!!
                                         ! and of special character...
                                         ! see fft99f.f for details....
        CALL set99(trigs2,ifax2,ny-1)    ! NOTE: ny must be ODD!!!!
                                         ! and of special character...
                                         ! see fft99f.f for details....
      END IF

    ELSE IF ( fftopt == 2 ) THEN   ! set up the cos fft arrays...

      IF(ny == 4)THEN   ! set up function in x direction only...

        CALL vcosti(nx-1,wsave2)         ! nx should be even.

      ELSE IF(nx == 4)THEN ! set up function in y direction only...

        CALL vcosti(ny-1,wsave1)         ! ny should be even.

      ELSE   ! set up functions for 2-d transform...

        CALL vcosti(ny-1,wsave1)         ! ny should be even.
        CALL vcosti(nx-1,wsave2)         ! nx should be even.

      END IF  ! end of run type if block...

    END IF   ! end of fftopt if block.....

  END IF
!
!-----------------------------------------------------------------------
!
!  Find the lowest model layer (index rayklow) that is entirely or
!  partially contained in the specified Rayleigh damping (sponge)
!  layers.
!
!  The Rayleigh damping is then applied only to layers with
!  k greater than or equal to rayklow.
!
!-----------------------------------------------------------------------
!
  rayklow = nz-1

  DO k=nz-1,2,-1

    zpmax = zp(1,1,k)

    DO j=1,ny-1
      DO i=1,nx-1
        zpmax = MAX( zp(i,j,k), zpmax )
      END DO
    END DO

! for bit-for-bit accuracy with MP version:
    rmin = zpmax
    call mpmax0(zpmax,rmin)

    IF( zpmax < zbrdmp ) THEN
      rayklow = MAX(2, k+1)
      EXIT
    END IF

  END DO
!
!-----------------------------------------------------------------------
!
!  Calculate Exner function and store in ppi
!
!-----------------------------------------------------------------------
!
  CALL setppi(nx,ny,nz,nts,tpresent,pprt,pbar,ppi)
!
!-----------------------------------------------------------------------
!
!  Calculate and store the sound wave speed squared in csndsq.
!
!-----------------------------------------------------------------------
!
  IF(csopt == 1) THEN       ! Original definition of sound speed.
    DO k=1,nz-1
      DO j=1,ny-1
        DO i=1,nx-1
          csndsq(i,j,k)= cpdcv*pbar(i,j,k)*j3(i,j,k)/rhostr(i,j,k)
        END DO
      END DO
    END DO

  ELSE IF(csopt == 2) THEN   ! Original sound speed multiplied
                             ! by a factor
    temp = csfactr**2*cpdcv
    DO k=1,nz-1
      DO j=1,ny-1
        DO i=1,nx-1
          csndsq(i,j,k)= temp * pbar(i,j,k)*j3(i,j,k)/rhostr(i,j,k)
        END DO
      END DO
    END DO
  ELSE                      ! Specified constant sound speed.
    DO k=1,nz-1
      DO j=1,ny-1
        DO i=1,nx-1
          csndsq(i,j,k)= csound * csound
        END DO
      END DO
    END DO
  END IF
!
!-----------------------------------------------------------------------
!
!  Fill the edges of base state arrays that are otherwise undefined.
!  This is done for safty reason only.
!
!-----------------------------------------------------------------------
!
  DO k=1,nz-1
    DO j=1,ny
      vbar  (nx,j,k) = vbar  (nx-1,j,k)
    END DO
  END DO

  DO k=1,nz-1
    DO i=1,nx
      ubar  (i,ny,k) = ubar  (i,ny-1,k)
    END DO
  END DO

  DO i=1,nx
    DO j=1,ny
      ubar  (i,j,nz) = ubar  (i,j,nz-1)
      vbar  (i,j,nz) = vbar  (i,j,nz-1)
    END DO
  END DO

  DO k=1,nz-1
    DO j=1,ny-1
      ptbar (nx,j,k) = ptbar (nx-1,j,k)
      pbar  (nx,j,k) = pbar  (nx-1,j,k)
      ppi   (nx,j,k) = ppi   (nx-1,j,k)
      qvbar (nx,j,k) = qvbar (nx-1,j,k)
      csndsq(nx,j,k) = csndsq(nx-1,j,k)
      rhostr(nx,j,k) = rhostr(nx-1,j,k)
    END DO
  END DO

  DO k=1,nz-1
    DO i=1,nx
      ptbar (i,ny,k) = ptbar (i,ny-1,k)
      pbar  (i,ny,k) = pbar  (i,ny-1,k)
      ppi   (i,ny,k) = ppi   (i,ny-1,k)
      qvbar (i,ny,k) = qvbar (i,ny-1,k)
      csndsq(i,ny,k) = csndsq(i,ny-1,k)
      rhostr(i,ny,k) = rhostr(i,ny-1,k)
    END DO
  END DO

  DO i=1,nx
    DO j=1,ny
      ptbar (i,j,nz) = ptbar (i,j,nz-1)
      pbar  (i,j,nz) = pbar  (i,j,nz-1)
      ppi   (i,j,nz) = ppi   (i,j,nz-1)
      qvbar (i,j,nz) = qvbar (i,j,nz-1)
      csndsq(i,j,nz) = csndsq(i,j,nz-1)
      rhostr(i,j,nz) = rhostr(i,j,nz-1)
    END DO
  END DO

  DO k=1,nz-1
    DO j=1,ny-1
      DO i=1,nx-1
        ptbari(i,j,k)   = 1.0/ptbar(i,j,k)
        pbari(i,j,k)    = 1.0/pbar(i,j,k)
        rhostri(i,j,k)  = 1.0/rhostr(i,j,k)
      END DO
    END DO
  END DO

  IF ( sfcphy > 0 ) THEN

    CALL initsfc(nx,ny,nz,nstyps,                                       &
                 pbar,pprt(1,1,1,1),                                    &
                 ptbar,ptprt(1,1,1,1),                                  &
                 qvbar,qv(1,1,1,1),                                     &
                 soiltyp,stypfrct, vegtyp, lai,roufns,veg,tem1,         &
                 tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth,qvsfc)

  END IF

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


SUBROUTINE initdvr(nx,ny,nz,nts,                                        & 1,9
           ubar,vbar,ptbar,pbar,rhostr,qvbar,                           &
           x,y,z,zp,hterain, j1,j2,j3,                                  &
           u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh,                          &
           ptcumsrc,qcumsrc,raing,rainc,prcrate,tem1)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Initialize the model time dependent variables.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Ming Xue
!  10/10/1991.
!
!  MODIFICATION HISTORY:
!
!  5/25/92 (M. Xue)
!  Added full documentation.
!
!  5/03/92 (M. Xue)
!  Further documentation.
!
!  4/10/93 (M. Xue & Hao Jin)
!  Add the terrain.
!
!  9/10/94 (Weygandt & Y. Lu)
!  Cleaned up documentation.
!
!  01/28/95 (G. Bassett)
!  Added pt0opt=5 (soup can shaped perturbation to ptprt).
!
!-----------------------------------------------------------------------
!
!  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
!    nts      Number of time levels to be initialized.
!
!    ubar     Base state x-velocity component (m/s)
!    vbar     Base state y-velocity component (m/s)
!    ptbar    Base state potential temperature (K)
!    pbar     Base state pressure (Pascal)
!    rhostr   Base state density rhobar times j3 (kg/m**3)
!    qvbar    Base state water vapor specific humidity (kg/kg).
!
!    x        x-coordinate of grid points in computational space (m)
!    y        y-coordinate of grid points in computational space (m)
!    z        z-coordinate of grid points in computational space (m)
!    zp       Vertical coordinate of grid points in physical space(m)
!    hterain  Terrain height (m)
!
!    j1       Coordinate transformation Jacobian -d(zp)/dx
!    j2       Coordinate transformation Jacobian -d(zp)/dy
!    j3       Coordinate transformation Jacobian  d(zp)/dz
!
!  OUTPUT:
!
!    u        x-component of velocity at all time levels (m/s).
!    v        y-component of velocity at all time levels (m/s).
!    w        z-component of velocity at all time levels (m/s).
!    ptprt    Perturbation potential temperature at all time levels
!             (K)
!    pprt     Perturbation pressure at all time levels (Pascal)
!    qv       Water vapor specific humidity at all time levels
!             (kg/kg)
!    qc       Cloud water mixing ratio at all time levels (kg/kg)
!    qr       Rainwater mixing ratio at all time levels (kg/kg)
!    qi       Cloud ice mixing ratio at all time levels (kg/kg)
!    qs       Snow mixing ratio at all time levels (kg/kg)
!    qh       Hail mixing ratio at all time levels (kg/kg)
!
!    ptcumsrc Source term in pt-equation due to cumulus parameterization
!    qcumsrc Source term in water equations due to cumulus parameterization
!    raing    Grid supersaturation rain
!    rainc    Cumulus convective rain
!    prcrate  Precipitation ratesrain
!
!  WORK ARRAYS:
!
!    tem1     Temporary work array.
!
!
!-----------------------------------------------------------------------
!

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


  INTEGER :: nx,ny,nz          ! The number of grid points in 3
                               ! directions
  INTEGER :: nts               ! Number of time levels to be initialized.
  INTEGER :: tpast             ! Index of time level for the past time.
  INTEGER :: tpresent          ! Index of time level for the present
                               ! time.
  INTEGER :: tfuture           ! Index of time level for the future
                               ! time.

  REAL :: ubar  (nx,ny,nz)     ! Base state u-velocity (m/s)
  REAL :: vbar  (nx,ny,nz)     ! Base state v-velocity (m/s)
  REAL :: ptbar (nx,ny,nz)     ! Base state potential temperature (K)
  REAL :: pbar  (nx,ny,nz)     ! Base state pressure (Pascal).
  REAL :: rhostr(nx,ny,nz)     ! Base state density rhobar times j3.
  REAL :: qvbar (nx,ny,nz)     ! Base state water vapor specific
                               ! humidity(kg/kg)

  REAL :: x     (nx)           ! x-coord. of the physical and compu-
                               ! tational grid. Defined at u-point.
  REAL :: y     (ny)           ! y-coord. of the physical and compu-
                               ! tational grid. Defined at v-point.
  REAL :: z     (nz)           ! z-coord. of the computational grid.
                               ! Defined at w-point on the staggered
                               ! grid.
  REAL :: zp    (nx,ny,nz)     ! Physical height coordinate defined at
                               ! w-point of the staggered grid.
  REAL :: hterain(nx,ny)       ! Terrain height (m).

  REAL :: j1    (nx,ny,nz)     ! Coordinate transformation Jacobian
                               ! -d(zp)/dx.
  REAL :: j2    (nx,ny,nz)     ! Coordinate transformation Jacobian
                               ! -d(zp)/dy.
  REAL :: j3    (nx,ny,nz)     ! Coordinate transformation Jacobian
                               ! d(zp)/dz.

  REAL :: u     (nx,ny,nz,nts) ! Total u-velocity (m/s)
  REAL :: v     (nx,ny,nz,nts) ! Total v-velocity (m/s)
  REAL :: w     (nx,ny,nz,nts) ! Total w-velocity (m/s)
  REAL :: ptprt (nx,ny,nz,nts) ! Perturbation potential temperature (K)
  REAL :: pprt  (nx,ny,nz,nts) ! Perturbation pressure (Pascal)

  REAL :: qv    (nx,ny,nz,nts) ! Water vapor specific humidity (kg/kg)
  REAL :: qc    (nx,ny,nz,nts) ! Cloud water mixing ratio (kg/kg)
  REAL :: qr    (nx,ny,nz,nts) ! Rain water mixing ratio (kg/kg)
  REAL :: qi    (nx,ny,nz,nts) ! Cloud ice mixing ratio (kg/kg)
  REAL :: qs    (nx,ny,nz,nts) ! Snow mixing ratio (kg/kg)
  REAL :: qh    (nx,ny,nz,nts) ! Hail mixing ratio (kg/kg)

  REAL :: ptcumsrc(nx,ny,nz)   ! Source term in pt-equation due
                               ! to cumulus parameterization
  REAL :: qcumsrc(nx,ny,nz,5)  ! Source term in water equations due
                               ! to cumulus parameterization:
                               ! qcumsrc(1,1,1,1) for qv equation
                               ! qcumsrc(1,1,1,2) for qc equation
                               ! qcumsrc(1,1,1,3) for qr equation
                               ! qcumsrc(1,1,1,4) for qi equation
                               ! qcumsrc(1,1,1,5) for qs equation
  REAL :: raing (nx,ny)        ! Grid supersaturation rain
  REAL :: rainc (nx,ny)        ! Cumulus convective rain
  REAL :: prcrate(nx,ny,4)     ! 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 :: tem1(nx,ny,nz)       ! Temporary work array

!
!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------
!
  REAL :: xs, ys, zs
  REAL :: us, vs, ws, rhobar
  REAL :: radnd , pi,pi2,pi4
  INTEGER :: i,j,k, n
  INTEGER :: iseed,ibgn,iend,jbgn,jend,kbgn,kend
  INTEGER :: ebc1,wbc1,nbc1,sbc1

  REAL :: amplitud
  REAL :: knumx,lnumy,mnumz
  REAL :: lnthx,lnthy,lnthz
  REAL :: lambda,lambdah,lambda2

  INTEGER :: mptag
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'globcst.inc'
  INCLUDE 'grid.inc'          ! Grid & map parameters.
  INCLUDE 'phycst.inc'
  INCLUDE 'bndry.inc'
  INCLUDE 'mp.inc'            ! Message passing parameters.
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  IF( nts > 1 ) THEN
    tpast=1
    tpresent=2
    tfuture=3
  ELSE
    tpast=1
    tpresent=1
    tfuture=1
  END IF
!
!-----------------------------------------------------------------------
!
!  Specify the initial potential temperature perturbation.
!
!-----------------------------------------------------------------------
!

  IF( pt0opt == 1 .OR. pt0opt == 6 ) THEN  ! Bubble shaped perturbation

!
!-----------------------------------------------------------------------
!
!  Define a potential temperature perturbation for a bubble-shaped
!  disturbance.
!
!-----------------------------------------------------------------------
!

    pi2 = 1.5707963267949

    IF( ptpert0 == 0.0) THEN

      DO k= 1,nz-1
        DO j= 1,ny-1
          DO i= 1,nx-1
            ptprt(i,j,k,1) = 0.0
          END DO
        END DO
      END DO


    ELSE

      DO k= 1,nz-1
        DO j= 1,ny-1
          DO i= 1,nx-1

!wdt multiple bubbles for timing tests:
            !xs = (x(i)+x(i+1))*0.5
            !ys = (y(j)+y(j+1))*0.5
            xs = (x(i)+x(i+1))*0.5-x(1)
            ys = (y(j)+y(j+1))*0.5-y(1)
            zs = (zp(i,j,k)+zp(i,j,k+1))*0.5

            IF( pt0rady < 0.0 .OR. runmod == 2 ) THEN
                                         ! 2-d bubble in x-z plane.

              radnd = SQRT( ((xs-pt0ctrx)/pt0radx)**2                   &
                         + ((zs-pt0ctrz)/pt0radz)**2 )

            ELSE IF( pt0radx < 0.0 .OR. runmod == 3 ) THEN
                                         ! 2-d bubble in y-z plane.

              radnd = SQRT( ((ys-pt0ctry)/pt0rady)**2                   &
                         + ((zs-pt0ctrz)/pt0radz)**2 )

            ELSE                         ! 3-d bubble

              radnd = SQRT( ((xs-pt0ctrx)/pt0radx)**2 +                 &
                           ((ys-pt0ctry)/pt0rady)**2 +                  &
                           ((zs-pt0ctrz)/pt0radz)**2 )
            END IF

            IF(radnd >= 1.0) THEN
              ptprt(i,j,k,1) = 0.0
            ELSE
              ptprt(i,j,k,1) = ptpert0*(COS(pi2*radnd )**2)
            END IF

          END DO
        END DO
      END DO

      IF(pt0opt == 6) THEN ! Perturbation speficied in T'.

        DO k= 1,nz-1
          DO j= 1,ny-1
            DO i= 1,nx-1
              ptprt(i,j,k,1) = ptprt(i,j,k,1)*(p0/pbar(i,j,k))**rddcp
            END DO
          END DO
        END DO

      END IF

    END IF

  ELSE IF( pt0opt == 2 .OR. pt0opt == 3 ) THEN
                                          ! Random field
!
!-----------------------------------------------------------------------
!
!  Define a potential temperature perturbation by a random function.
!  This ensures that the horizontal average of the perturbation in a
!  specified domain is zero.
!
!  Fill an array with zeros.
!
!-----------------------------------------------------------------------
!
    DO k=1,nz-1
      DO j=1,ny-1
        DO i=1,nx-1
          ptprt(i,j,k,1) = 0.0
        END DO
      END DO
    END DO
!
!-----------------------------------------------------------------------
!
!    The following parameters define the portion of domain to be
!    initialized with a random potential temperture perturbation.
!    Users can modify them to fit their needs.
!
!-----------------------------------------------------------------------
!
    ibgn = 1
    iend = nx-1
    jbgn = 1
    jend = ny-1
    kbgn = 1
    kend = nz-1

    iseed = -100

    DO k= kbgn,kend

      CALL ranary(nx,ny,ibgn,iend,jbgn,jend,                            &
                  iseed,ptpert0,ptprt(1,1,k,1))

    END DO

    IF( pt0opt == 3 ) THEN ! Symmetric random perturbation
!
!-----------------------------------------------------------------------
!
!  Create a random perturbation field symmetric about both the x and y
!  axes.
!
!-----------------------------------------------------------------------

      DO k=1,nz-1
        DO j=1,ny-1
          DO i=1,nx/2
            ptprt(i,j,k,1) = ptprt(nx-i,j,k,1)
          END DO
        END DO
      END DO

      DO k=1,nz-1
        DO i=1,nx-1
          DO j=1,ny/2
            ptprt(i,j,k,1) = ptprt(i,ny-j,k,1)
          END DO
        END DO
      END DO

      IF( nx == ny ) THEN
        DO k=1,nz-1
          DO j=1,ny-1
            DO i=1,nx-1
              tem1(i,j,k) = (ptprt(i,j,k,1)+ptprt(j,i,k,1))*0.5
            END DO
          END DO
        END DO

        DO k=1,nz-1
          DO i=1,nx-1
            DO j=1,ny-1
              ptprt(i,j,k,1) = tem1(i,j,k)
            END DO
          END DO
        END DO
      END IF

    END IF

  ELSE IF( pt0opt == 4 ) THEN       ! Bubble given in Skamarock and
                                    ! Klemp (1994)
    pi2 = 1.5707963267949

    DO k=1,nz-1
      DO i=1,nx-1
        DO j=1,ny-1
          xs = (x(i)+x(i+1))*0.5
          zs = (zp(i,j,k)+zp(i,j,k+1))*0.5
          ptprt(i,j,k,1) = ptpert0                                      &
              *SIN(pi2*2*zs/((nz-3)*dz))/(1+((xs-pt0ctrx)/pt0radx)**2)
        END DO
      END DO
    END DO

  ELSE IF( pt0opt == 5 ) THEN       ! Soup can shaped perturbation

    DO k= 1,nz-1
      DO j= 1,ny-1
        DO i= 1,nx-1

          xs = (x(i)+x(i+1))*0.5
          ys = (y(j)+y(j+1))*0.5
          zs = (zp(i,j,k)+zp(i,j,k+1))*0.5
          IF ( ABS(zs-pt0ctrz)/pt0radz < 1 ) THEN

            IF( pt0rady < 0.0 .OR. runmod == 2 ) THEN
                                         ! 2-d bubble in x-z plane.

              radnd = ABS(xs-pt0ctrx)/pt0radx

            ELSE IF( pt0radx < 0.0 .OR. runmod == 3 ) THEN
                                         ! 2-d bubble in y-z plane.

              radnd = ABS(ys-pt0ctry)/pt0rady

            ELSE                         ! 3-d bubble

              radnd = SQRT( ((xs-pt0ctrx)/pt0radx)**2 +                 &
                            ((ys-pt0ctry)/pt0rady)**2 )

            END IF

            IF(radnd >= 1.0) THEN
              ptprt(i,j,k,1) = 0.0
            ELSE
              ptprt(i,j,k,1) = ptpert0
            END IF

          ELSE

            ptprt(i,j,k,1) = 0.0

          END IF

        END DO
      END DO
    END DO

  END IF

  ebc1=0
  wbc1=0
  sbc1=0
  nbc1=0

  IF( ebc == 1 .OR.ebc == 2 .OR. ebc == 3 )  ebc1=ebc
  IF( wbc == 1 .OR.wbc == 2 .OR. wbc == 3 )  wbc1=wbc
  IF( sbc == 1 .OR.sbc == 2 .OR. sbc == 3 )  sbc1=sbc
  IF( nbc == 1 .OR.nbc == 2 .OR. nbc == 3 )  nbc1=nbc

!call test_dump (ptprt,"XXXinit_ptprt",nx,ny,nz,0,0)
  IF (mp_opt > 0) THEN
    CALL acct_interrupt(mp_acct)
    CALL mpsend2dew(ptprt,nx,ny,nz,ebc1,wbc1,0,mptag,tem1)
    CALL mprecv2dew(ptprt,nx,ny,nz,ebc1,wbc1,0,mptag,tem1)
    CALL mpsend2dns(ptprt,nx,ny,nz,nbc1,sbc1,0,mptag,tem1)
    CALL mprecv2dns(ptprt,nx,ny,nz,nbc1,sbc1,0,mptag,tem1)
  END IF
  CALL acct_interrupt(bc_acct)
  CALL bcsclr(nx,ny,nz,dtbig,                                           &
              ptprt(1,1,1,1),ptprt(1,1,1,1),                            &
              ptprt(1,1,1,1),tem1,tem1,tem1,tem1,                       &
              ebc1,wbc1,nbc1,sbc1,tbc,bbc)
  CALL acct_stop_inter
!call test_dump (ptprt,"XXXAinit_ptprt",nx,ny,nz,0,1)

  DO n=1,nts
    DO k=1,nz-1
      DO j=1,ny-1
        DO i=1,nx-1
          ptprt(i,j,k,n)=ptprt(i,j,k,1)
          pprt(i,j,k,n)=0.0
        END DO
      END DO
    END DO
  END DO

  DO n=1,nts
    DO k=1,nz-1
      DO j=1,ny-1
        DO i=1,nx
          u(i,j,k,n)=ubar(i,j,k)
        END DO
      END DO
    END DO
  END DO

  DO n=1,nts
    DO k=1,nz-1
      DO j=1,ny
        DO i=1,nx-1
          v(i,j,k,n)=vbar(i,j,k)
        END DO
      END DO
    END DO
  END DO

  DO n=1,nts
    DO k=1,nz
      DO j=1,ny-1
        DO i=1,nx-1
          w(i,j,k,n)=0.0
        END DO
      END DO
    END DO
  END DO


  DO n=1,nts
    DO k=1,nz-1
      DO j=1,ny-1
        DO i=1,nx-1
          qv(i,j,k,n)= qvbar(i,j,k)
        END DO
      END DO
    END DO
  END DO

  DO n=1,nts
    DO k=1,nz-1
      DO j=1,ny-1
        DO i=1,nx-1
          qc(i,j,k,n)=0.0
        END DO
      END DO
    END DO
  END DO

  DO n=1,nts
    DO k=1,nz-1
      DO j=1,ny-1
        DO i=1,nx-1
          qr(i,j,k,n)= 0.0
        END DO
      END DO
    END DO
  END DO

  DO n=1,nts
    DO k=1,nz-1
      DO j=1,ny-1
        DO i=1,nx-1
          qi(i,j,k,n)= 0.0
        END DO
      END DO
    END DO
  END DO

  DO n=1,nts
    DO k=1,nz-1
      DO j=1,ny-1
        DO i=1,nx-1
          qs(i,j,k,n)= 0.0
        END DO
      END DO
    END DO
  END DO

  DO n=1,nts
    DO k=1,nz-1
      DO j=1,ny-1
        DO i=1,nx-1
          qh(i,j,k,n)= 0.0
        END DO
      END DO
    END DO
  END DO

  DO k=1,nz-1
    DO j=1,ny-1
      DO i=1,nx-1
        ptcumsrc(i,j,k)= 0.0
        qcumsrc(i,j,k,1)= 0.0
        qcumsrc(i,j,k,2)= 0.0
        qcumsrc(i,j,k,3)= 0.0
        qcumsrc(i,j,k,4)= 0.0
        qcumsrc(i,j,k,5)= 0.0
      END DO
    END DO
  END DO

  DO j=1,ny-1
    DO i=1,nx-1
      raing(i,j)= 0.0
      rainc(i,j)= 0.0
      prcrate(i,j,1)= 0.0
      prcrate(i,j,2)= 0.0
      prcrate(i,j,3)= 0.0
      prcrate(i,j,4)= 0.0
    END DO
  END DO

!
!-----------------------------------------------------------------------
!
!  The following setup is to overwrite the total u, v, w, pprt,
!  ptprt, and qv for Beltrami flow initial conditions.
!
!-----------------------------------------------------------------------
!
  IF ( pt0opt == 7 ) THEN

    pi = 3.1415926535898
    pi2 = 2*pi
    pi4 = 4*pi

    amplitud = 2.0                      ! amplitude A=2 m/s
    tmixopt  = 1                        ! constant viscosity option
    tmixcst  = 1.0                      ! viscosity=1 m**2/s
    tmixvert = 1.0                      ! compute all mixing terms

    wbc = 2                             ! reset boundary conditions
    ebc = 2                             ! to periodical condition
    nbc = 2
    sbc = 2
    tbc = 2
    bbc = 2

    lnthx = dx*(nx-3)                   ! length in x
    lnthy = dy*(ny-3)                   ! length in y
    lnthz = dz*(nz-3)                   ! length in z

    knumx = pi4/lnthx                   ! wave number in x-dir
    lnumy = pi2/lnthy                   ! wave number in y-dir
    mnumz = pi2/lnthz                   ! wave number in z-dir

    lambdah = knumx*knumx + lnumy*lnumy
    lambda2 = lambdah + mnumz*mnumz
    lambda  = SQRT( lambda2 )

!    print *, ' amplitude = ',amplitud
    PRINT *, ' lnthx   = ',lnthx,                                       &
             ' lnthy   = ',lnthy,                                       &
             ' lnthz   = ',lnthz
    PRINT *, ' knumx   = ',knumx,                                       &
             ' lnumy   = ',lnumy,                                       &
             ' mnumz   = ',mnumz
    PRINT *, ' lambda1 = ',lambdah,                                     &
             ' lambda2 = ',lambda2,                                     &
             ' lambda  = ',lambda

    DO k=1,nz-1
      DO j=1,ny-1
        DO i=1,nx
          ys = 0.5*(y(j+1)+y(j))
          zs = 0.5*(z(k+1)+z(k))
          u(i,j,k,1) = -amplitud/lambdah                                &
                     *(lambda*lnumy                                     &
                      *COS(knumx*x(i))*SIN(lnumy*ys)*SIN(mnumz*zs)      &
                      +mnumz*knumx                                      &
                      *SIN(knumx*x(i))*COS(lnumy*ys)*COS(mnumz*zs))
        END DO
      END DO
    END DO

    DO k=1,nz-1
      DO j=1,ny
        DO i=1,nx-1
          xs = 0.5*(x(i+1)+x(i))
          zs = 0.5*(z(k+1)+z(k))
          v(i,j,k,1)= amplitud/lambdah                                  &
                    * (lambda*knumx                                     &
                      *SIN(knumx*xs)*COS(lnumy*y(j))*SIN(mnumz*zs)      &
                      -mnumz*lnumy                                      &
                      *COS(knumx*xs)*SIN(lnumy*y(j))*COS(mnumz*zs))
        END DO
      END DO
    END DO

    DO k=1,nz
      DO j=1,ny-1
        DO i=1,nx-1
          xs = 0.5*(x(i+1)+x(i))
          ys = 0.5*(y(j+1)+y(j))
          w(i,j,k,1)=amplitud                                           &
                    *COS(knumx*xs)*COS(lnumy*ys)*SIN(mnumz*z(k))
        END DO
      END DO
    END DO

    DO k=1,nz
      DO j=1,ny-1
        DO i=1,nx-1
          zs = 0.5*(z(k+1)+z(k))

          us = 0.5*(u(i+1,j,k,1)+u(i,j,k,1))
          vs = 0.5*(v(i,j+1,k,1)+v(i,j,k,1))
          ws = 0.5*(w(i,j,k+1,1)+w(i,j,k,1))

          rhobar = rhostr(i,j,k)/j3(i,j,k)
          pprt(i,j,k,1) = p0-rhobar*(0.5*(us*us+vs*vs+ws*ws)+g*zs)      &
                        - pbar(i,j,k)
        END DO
      END DO
    END DO

    DO n=1,nts
      DO k=1,nz
        DO j=1,ny
          DO i=1,nx
            u    (i,j,k,n) = u(i,j,k,1)
            v    (i,j,k,n) = v(i,j,k,1)
            w    (i,j,k,n) = w(i,j,k,1)
            pprt (i,j,k,n) = pprt(i,j,k,1)
            ptprt(i,j,k,n) = 0.0
            qv   (i,j,k,n) = 0.0
          END DO
        END DO
      END DO
    END DO

  END IF

  RETURN
END SUBROUTINE initdvr

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


SUBROUTINE extinit(nx,ny,nz,nts, nstyps,                                & 1,36
           u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh,tke,                      &
           ubar,vbar,ptbar,pbar,rhostr,qvbar,                           &
           x,y,z,zp,hterain,j1,j2,j3,                                   &
           soiltyp,stypfrct,vegtyp,lai,roufns,veg,                      &
           tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth,qvsfc,              &
           ptcumsrc,qcumsrc,raing,rainc,prcrate,                        &
           radfrc,radsw,rnflx,                                          &
           usflx,vsflx,ptsflx,qvsflx,                                   &
           uprt,vprt,qvprt,kmh,kmv,wbar,                                &
           tem6,tem7,tem8)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Read in the initial fields from an externally supplied initial
!  data file.
!
!  These fields include u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh
!  at time level tpresent, and the base state variables
!  ubar,vbar,ptbar,pbar,rhostr,qvbar.
!
!  Fields at tpast and tfuture are set to their values at tpresent.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Ming Xue
!  10/7/1992
!
!  MODIFICATION HISTORY:
!
!  3/25/94 (G. Bassett, M. Xue)
!  Modified to read in regular binary history dumps.
!
!  9/10/94 (Weygandt & Y. Lu)
!  Cleaned up documentation.
!
!  01/14/1995 (M. Xue)
!  Corrected where jacob was called. It should be called
!  before do loop 90.
!
!  03/29/1997 (M. Xue)
!  Modification made so that when values of mapproj,trulat1,trulat2,
!  trulon,ctrlat,ctrlon in the input data do not match those in
!  input file, the values in the data are used instead of those
!  in input, as was done before. Warning messages will be printed.
!
!  12/09/1998 (Donghai Wang)
!  Added the snow cover.
!
!  03/19/2002 (Keith Brewster)
!  Corrected print statement related to mis-match in data and user times.
!
!-----------------------------------------------------------------------
!
!  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
!    nts      Number of time levels to be initialized.
!
!  OUTPUT:
!
!    u        x component of velocity at times tpast and tpresent
!             (m/s)
!    v        y component of velocity at times tpast and tpresent
!             (m/s)
!    w        Vertical component of Cartesian velocity at times
!             tpast and tpresent (m/s)
!    ptprt    Perturbation potential temperature at times tpast and
!             tpresent (K)
!    pprt     Perturbation pressure at times tpast and tpresent
!             (Pascal)
!
!    qv       Water vapor specific humidity at times tpast and
!             tpresent (kg/kg)
!    qc       Cloud water mixing ratio at times tpast and tpresent
!             (kg/kg)
!    qr       Rainwater mixing ratio at times tpast and tpresent
!             (kg/kg)
!    qi       Cloud ice mixing ratio at times tpast and tpresent
!             (kg/kg)
!    qs       Snow mixing ratio at times tpast and tpresent (kg/kg)
!    qh       Hail mixing ratio at times tpast and tpresent (kg/kg)
!    tke      Turbulence kinetic energy (m**2/s)
!
!    ubar     Base state zonal velocity component (m/s)
!    vbar     Base state meridional velocity component (m/s)
!    ptbar    Base state potential temperature (K)
!    pbar     Base state pressure (Pascal)
!    rhostr   Base state density rhobar times j3 (kg/m**3)
!    qvbar    Base state water vapor specific humidity (kg/kg)
!
!    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 (m)
!    zp       Vertical coordinate of grid points in physical space(m)
!    hterain  The height of terrain (m)
!
!    j1       Coordinate transformation Jacobian -d(zp)/dx
!    j2       Coordinate transformation Jacobian -d(zp)/dy
!    j3       Coordinate transformation Jacobian  d(zp)/dz
!
!    soiltyp  Soil type
!    vegtyp   Vegetation type
!    lai      Leaf Area Index
!    roufns   Surface roughness
!    veg      Vegetation fraction
!
!    tsfc     Skin temperature at the ground or ocean surface (K)
!    tsoil    Deep soil temperature (K) (in deep 1 m layer)
!    wetsfc   Surface soil moisture
!    wetdp    Deep soil moisture
!    wetcanp  Canopy water amount
!    snowdpth Snow depth (m)
!    qvsfc    Effective S.H. at sfc.
!
!    ptcumsrc Source term in pt-equation due to cumulus parameterization
!    qcumsrc Source term in water equations due to cumulus parameterization
!    raing    Grid supersaturation rain
!    rainc    Cumulus convective rain
!    prcrate  Precipitation rates
!
!    radfrc   Radiation forcing (K/s)
!    radsw    Solar radiation reaching the surface
!    rnflx    Net radiation flux absorbed by surface
!
!    usflx    Surface flux of u-momentum (kg/(m*s**2))
!    vsflx    Surface flux of v-momentum (kg/(m*s**2))
!    ptsflx   Surface heat flux (K*kg/(m**2 * s ))
!    qvsflx   Surface moisture flux of (kg/(m**2 * s))
!
!    tstart   The time when the time integration starts, which is set
!             to the time of the restart data
!
!  WORK ARRAYS:
!
!    tem6     Temporary work array.
!    tem7     Temporary work array.
!    tem8     Temporary work array.
!
!-----------------------------------------------------------------------
!

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

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

  INTEGER :: nts               ! Number of time levels to be initialized.
  INTEGER :: tpast             ! Index of time level for the past time.
  INTEGER :: tpresent          ! Index of time level for the present
                               ! time.
  INTEGER :: tfuture           ! Index of time level for the future
                               ! time.

  REAL :: u     (nx,ny,nz,nts)  ! Total u-velocity (m/s)
  REAL :: v     (nx,ny,nz,nts) ! Total v-velocity (m/s)
  REAL :: w     (nx,ny,nz,nts) ! Total w-velocity (m/s)
  REAL :: ptprt (nx,ny,nz,nts) ! Perturbation potential temperature (K)
  REAL :: pprt  (nx,ny,nz,nts) ! Perturbation pressure (Pascal)

  REAL :: qv    (nx,ny,nz,nts) ! Water vapor specific humidity (kg/kg)
  REAL :: qc    (nx,ny,nz,nts) ! Cloud water mixing ratio (kg/kg)
  REAL :: qr    (nx,ny,nz,nts) ! Rain water mixing ratio (kg/kg)
  REAL :: qi    (nx,ny,nz,nts) ! Cloud ice mixing ratio (kg/kg)
  REAL :: qs    (nx,ny,nz,nts) ! Snow mixing ratio (kg/kg)
  REAL :: qh    (nx,ny,nz,nts) ! Hail mixing ratio (kg/kg)
  REAL :: tke   (nx,ny,nz,nts) ! Turbulence kinetic energy (m**2/s)

  REAL :: ubar  (nx,ny,nz)     ! Base state u-velocity (m/s)
  REAL :: vbar  (nx,ny,nz)     ! Base state v-velocity (m/s)
  REAL :: ptbar (nx,ny,nz)     ! Base state potential temperature (K)
  REAL :: pbar  (nx,ny,nz)     ! Base state pressure (Pascal)
  REAL :: rhostr(nx,ny,nz)     ! Base state density rhobar times j3.
  REAL :: qvbar (nx,ny,nz)     ! Base state water vapor specific
                               ! humidity(kg/kg)

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

  REAL :: hterain(nx,ny)       ! Terrain height (m).

  REAL :: j1    (nx,ny,nz)     ! Coordinate transformation Jacobian
                               ! -d(zp)/dx.
  REAL :: j2    (nx,ny,nz)     ! Coordinate transformation Jacobian
                               ! -d(zp)/dy.
  REAL :: j3    (nx,ny,nz)     ! Coordinate transformation Jacobian
                               ! d(zp)/dz.

  INTEGER :: nstyps                  ! Number of soil types
  INTEGER :: soiltyp(nx,ny,nstyps)   ! Soil types at grid points
  REAL :: stypfrct(nx,ny,nstyps)  ! Fraction of soil types
  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) ! Temperature at ground (K)
                                  ! (in top 1 cm layer)
  REAL :: qvsfc  (nx,ny,0:nstyps) ! Effective qv at sfc.
  REAL :: tsoil  (nx,ny,0:nstyps) ! Deep soil temperature (K)
                                  ! (in deep 1 m layer)
  REAL :: wetsfc (nx,ny,0:nstyps) ! Surface soil moisture in the top
                                  ! 1 cm layer
  REAL :: wetdp  (nx,ny,0:nstyps) ! Deep soil moisture in the deep
                                  ! 1 m layer
  REAL :: wetcanp(nx,ny,0:nstyps) ! Canopy water amount
  REAL :: snowdpth(nx,ny)         ! Snow depth (m)


  REAL :: ptcumsrc(nx,ny,nz)   ! Source term in pt-equation due
                               ! to cumulus parameterization
  REAL :: qcumsrc(nx,ny,nz,5)  ! Source term in water equations due
                               ! to cumulus parameterization:
                               ! qcumsrc(1,1,1,1) for qv equation
                               ! qcumsrc(1,1,1,2) for qc equation
                               ! qcumsrc(1,1,1,3) for qr equation
                               ! qcumsrc(1,1,1,4) for qi equation
                               ! qcumsrc(1,1,1,5) for qs equation
  REAL :: raing(nx,ny)         ! Grid supersaturation rain
  REAL :: rainc(nx,ny)         ! Cumulus convective rain
  REAL :: prcrate(nx,ny,4)     ! 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 :: radfrc(nx,ny,nz)     ! Radiation forcing (K/s)
  REAL :: radsw (nx,ny)        ! Solar radiation reaching the surface
  REAL :: rnflx (nx,ny)        ! Net radiation flux absorbed by surface

  REAL :: usflx (nx,ny)        ! Surface flux of u-momentum (kg/(m*s**2))
  REAL :: vsflx (nx,ny)        ! Surface flux of v-momentum (kg/(m*s**2))
  REAL :: ptsflx(nx,ny)        ! Surface heat flux (K*kg/(m*s**2))
  REAL :: qvsflx(nx,ny)        ! Surface moisture flux (kg/(m**2*s))

  REAL :: uprt (nx,ny,nz)      ! Temporary array
  REAL :: vprt (nx,ny,nz)      ! Temporary array
  REAL :: qvprt(nx,ny,nz)      ! Temporary array
  REAL :: kmh   (nx,ny,nz)     ! Horizontal turb. mixing coef. for
                               ! momentum. ( m**2/s )
  REAL :: kmv   (nx,ny,nz)     ! Vertical turb. mixing coef. for
                               ! momentum. ( m**2/s )
  REAL :: wbar (nx,ny,nz)      ! Temporary array
  REAL :: tem6 (nx,ny,nz)      ! Temporary array
  REAL :: tem7 (nx,ny,nz)      ! Temporary array
  REAL :: tem8 (nx,ny,nz)      ! Temporary array
!
!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------
!
  INTEGER :: lengbf, lenfil
  INTEGER :: i, j, k, ireturn, tim
  REAL :: amax,amin
  CHARACTER (LEN=80) :: runnamesv
  INTEGER :: nocmnt_old          ! Number of comment lines
  CHARACTER (LEN=80 ) :: cmnt_old(50)  ! String of comments on this model run
  INTEGER :: nchin
  REAL :: time_tmp
  REAL :: tstopsv,thisdmpsv,mapprojsv,latitudsv,ctrlatsv,ctrlonsv
  INTEGER :: monthsv,daysv,yearsv,hoursv,minutesv,secondsv
  REAL :: trulat1sv,trulat2sv,trulonsv
  REAL :: dxsv,dysv,strhoptsv,dzsv,dzminsv,zrefsfcsv,dlayer1sv,dlayer2sv,  &
          strhtunesv,zflatsv
  REAL :: p0inv,eps,tvbar

  INTEGER :: abstsec0,abstsec1
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'globcst.inc'
  INCLUDE 'grid.inc'          ! Grid & map parameters.
  INCLUDE 'phycst.inc'
  INCLUDE 'mp.inc'            ! Message passing parameters.
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  p0inv=1./p0

  IF( nts > 1 ) THEN
    tpast=1
    tpresent=2
    tfuture=3
  ELSE
    tpast=1
    tpresent=1
    tfuture=1
  END IF
!
!-----------------------------------------------------------------------
!
!  Read in the initial fields from the input data file.
!
!-----------------------------------------------------------------------
!
  lengbf = 80
  CALL strlnth( inigbf, lengbf )

  lenfil = 80
  CALL strlnth( inifile, lenfil )

  IF (mp_opt > 0) THEN
    runnamesv = inigbf
    WRITE(inigbf,'(a,a,2i2.2)')                                         &
        runnamesv(1:lengbf),'_',loc_x,loc_y
    lengbf = lengbf + 5
    runnamesv = inifile
    WRITE(inifile,'(a,a,2i2.2)')                                        &
        runnamesv(1:lenfil),'_',loc_x,loc_y
    lenfil = lenfil + 5
  END IF

  tim = tpresent

  runnamesv = runname
  tstopsv = tstop
  thisdmpsv = thisdmp
  mapprojsv = mapproj
  trulat1sv = trulat1
  trulat2sv = trulat2
  trulonsv  = trulon
  latitudsv = latitud
  ctrlatsv = ctrlat
  ctrlonsv = ctrlon
  monthsv = month
  daysv = day
  yearsv = year
  hoursv = hour
  minutesv = minute
  secondsv = second
  dxsv = dx
  dysv = dy
  strhoptsv = strhopt
  dzsv = dz
  dzminsv = dzmin
  zrefsfcsv = zrefsfc
  dlayer1sv = dlayer1
  dlayer2sv = dlayer2
  strhtunesv = strhtune
  zflatsv = zflat

  nocmnt_old = nocmnt
  DO i=1,nocmnt_old
    cmnt_old(i)=cmnt(i)
  END DO

!  blocking inserted for ordering i/o for message passing
  DO i=0,nprocs-1,max_fopen
    IF(myproc >= i.AND.myproc <= i+max_fopen-1)THEN

      CALL dtaread(nx,ny,nz,nstyps, inifmt, nchin ,inigbf,lengbf,       &
                 inifile,lenfil,time_tmp,                               &
                 x,y,z,zp,uprt,vprt,w(1,1,1,tim),ptprt(1,1,1,tim),      &
                 pprt(1,1,1,tim),qvprt,qc(1,1,1,tim),qr(1,1,1,tim),     &
                 qi(1,1,1,tim),qs(1,1,1,tim),qh(1,1,1,tim),             &
                 tke(1,1,1,tim),kmh,kmv,                                &
                 ubar,vbar,wbar,ptbar,pbar,rhostr,qvbar,                &
                 soiltyp,stypfrct,vegtyp,lai,roufns,veg,                &
                 tsfc, tsoil, wetsfc, wetdp,wetcanp,snowdpth,           &
                 raing,rainc,prcrate,                                   &
                 radfrc,radsw,rnflx,                                    &
                 usflx,vsflx,ptsflx,qvsflx,                             &
                 ireturn, tem6,tem7,tem8)

    END IF
    IF (mp_opt > 0) CALL mpbarrier
  END DO
!

!-----------------------------------------------------------------------
!
!  To restore the original runname and comments, etc.
!  from ARPS input file.
!
!-----------------------------------------------------------------------

  eps = 0.0001

  IF(mapproj /= mapprojsv .OR.ABS(trulat1sv-trulat1) > eps              &
        .OR.ABS(trulat2sv-trulat2) > eps                                &
        .OR.ABS(trulonsv -trulon ) > eps) THEN
    WRITE(6,'(/2(/5x,a),2(/5x,3(a,f13.3))/)')                           &
        'WARNING: Map projection parameters in the input data do not ', &
        'match those specified in input file. Values in data used.',    &
        'In data,  trulat1=',trulat1  ,', trulat2=',trulat2,            &
        ', trulon=',trulon,                                             &
        'In input, trulat1=',trulat1sv,', trulat2=',trulat2sv,          &
        ', trulon=',trulonsv
  END IF

  IF(ABS(ctrlatsv-ctrlat) > eps .OR. ABS(ctrlonsv-ctrlon) > eps ) THEN
    WRITE(6,'(/3(/5x,a),2(/5x,2(a,f13.3))/)')                           &
        'WARNING: Central latitude and/or longitude of the grid ',      &
        'in the input data do not match those specified in input file.', &
        'Values in data used.',                                         &
        'In data , ctrlat=',ctrlat  ,', ctrlon=',ctrlon,                &
        'In input, ctrlat=',ctrlatsv,', ctrlon=',ctrlonsv
  END IF

  CALL ctim2abss(year,month,day,hour,minute,second, abstsec0)
  CALL ctim2abss(yearsv,monthsv,daysv,hoursv,minutesv,secondsv,         &
                 abstsec1)

  abstsec0 = abstsec0 + nint(time_tmp)
  abstsec1 = abstsec1 + nint(tstart)

  IF ( abstsec0 /= abstsec1 ) THEN
    WRITE(6,'(a,2(/a,1x,i4.4,5(a,i2.2),a,f10.3))')                      &
        'WARNING: Data time is inconsistent with user input time.',     &
        '         Data initime:', year,'-',month,'-',day,'.',          &
                                 hour,':',minute,':',second,            &
                                 '     model data time: ', tstart,      &
        '         User initime:', yearsv,'-',monthsv,'-',daysv,'.',  &
                                 hoursv,':',minutesv,':',secondsv,      &
                                 ' model starting time: ', time_tmp

    IF ( timeopt == 1 ) THEN
      WRITE(6,'(a)') 'Program continues using user specified time'
      year = yearsv
      month = monthsv
      day = daysv
      hour = hoursv
      minute = minutesv
      second = secondsv
    ELSE IF ( timeopt == 2 ) THEN
      WRITE(6,'(a)') 'Program continues using data time'
      tstart = time_tmp
    ELSE
      WRITE(6,'(a)') 'Program stopped in subroutine EXTINIT'
      CALL arpsstop ("arpsstop called from extinit due to timeopt",1)
    END IF
  ELSE
    year = yearsv
    month = monthsv
    day = daysv
    hour = hoursv
    minute = minutesv
    second = secondsv
    IF (myproc == 0) &
    WRITE(6,'(a,i4.4,5(a,i2.2),3x,a,f10.3,a)')                          &
         'Use specified initial time: ',                                &
          year,'-',month,'-',day,'.',hour,':',minute,':',second,        &
         'and model starting time: ', tstart, 'seconds'
  END IF

  runname = runnamesv
  tstop = tstopsv
  thisdmp = thisdmpsv
  dx = dxsv
  dy = dysv
  strhopt = strhoptsv
  dz = dzsv
  dzmin = dzminsv
  zrefsfc = zrefsfcsv
  dlayer1 = dlayer1sv
  dlayer2 = dlayer2sv
  strhtune = strhtunesv
  zflat = zflatsv

  nocmnt = nocmnt_old
  DO i=1,nocmnt
    cmnt(i)=cmnt_old(i)
  END DO

  IF( ireturn /= 0) THEN

    WRITE(6,'(3(/1x,a))')                                               &
        'Error occured when reading history data ',                     &
        inifile(1:lenfil)//' or '//inigbf(1:lengbf),                    &
        'Program stopped in EXTINIT.'
      CALL arpsstop ("arpsstop called from extinit in reading file",1)

  END IF

!
!-----------------------------------------------------------------------
!
!  Calculate rhostr & jacobians, set hterain.
!
!-----------------------------------------------------------------------
!

  CALL jacob(nx,ny,nz,x,y,z,zp,j1,j2,j3)

  DO k= 1,nz-1
    DO j= 1,ny-1
      DO i= 1,nx-1
        tvbar=(ptbar(i,j,k)*((pbar(i,j,k)*p0inv)**rddcp))*              &
            ((1.0+rvdrd*qvbar(i,j,k))/(1.0+qvbar(i,j,k)))
        rhostr(i,j,k)= ABS(j3(i,j,k))*pbar(i,j,k)/(rd*tvbar)
!
!       rhostr(i,j,k)= abs(j3(i,j,k))*pbar(i,j,k)/
!    :         ( Rd * ptbar(i,j,k)*(pbar(i,j,k)*p0inv)**rddcp )
!
      END DO
    END DO
  END DO

  DO j=1,ny-1
    DO i=1,nx-1
      rhostr(i,j,   1)=rhostr(i,j,2)
      rhostr(i,j,nz-1)=rhostr(i,j,nz-2)
    END DO
  END DO

  DO i=1,nx
    DO j=1,ny
      hterain(i,j) = zp(i,j,2)
    END DO
  END DO


!
!-----------------------------------------------------------------------
!
!  Put the 3-d variables into their 4-d arrays.
!
!-----------------------------------------------------------------------
!
  tim = tpresent

  DO k=1,nz-1
    DO j=1,ny-1
      DO i=1,nx
        u(i,j,k,tim)=uprt(i,j,k)+ubar(i,j,k)
      END DO
    END DO
  END DO


  DO k=1,nz-1
    DO j=1,ny
      DO i=1,nx-1
        v(i,j,k,tim)=vprt(i,j,k)+vbar(i,j,k)
      END DO
    END DO
  END DO

  DO k=1,nz-1
    DO j=1,ny-1
      DO i=1,nx-1
        qv(i,j,k,tim)=qvprt(i,j,k)+qvbar(i,j,k)
      END DO
    END DO
  END DO

!
!-----------------------------------------------------------------------
!
!  Set the future values of variables to their current values.
!  This is done primarily for safety reasons.
!
!-----------------------------------------------------------------------
!
  IF(nts > 1) THEN
    DO k=1,nz
      DO j=1,ny
        DO i=1,nx
          u    (i,j,k,tfuture) = u    (i,j,k,tpresent)
          v    (i,j,k,tfuture) = v    (i,j,k,tpresent)
          w    (i,j,k,tfuture) = w    (i,j,k,tpresent)
          ptprt(i,j,k,tfuture) = ptprt(i,j,k,tpresent)
          pprt (i,j,k,tfuture) = pprt (i,j,k,tpresent)
          qv   (i,j,k,tfuture) = qv   (i,j,k,tpresent)
          qc   (i,j,k,tfuture) = qc   (i,j,k,tpresent)
          qr   (i,j,k,tfuture) = qr   (i,j,k,tpresent)
          qi   (i,j,k,tfuture) = qi   (i,j,k,tpresent)
          qs   (i,j,k,tfuture) = qs   (i,j,k,tpresent)
          qh   (i,j,k,tfuture) = qh   (i,j,k,tpresent)
          tke  (i,j,k,tfuture) = tke  (i,j,k,tpresent)

          u    (i,j,k,tpast  ) = u    (i,j,k,tpresent)
          v    (i,j,k,tpast  ) = v    (i,j,k,tpresent)
          w    (i,j,k,tpast  ) = w    (i,j,k,tpresent)
          ptprt(i,j,k,tpast  ) = ptprt(i,j,k,tpresent)
          pprt (i,j,k,tpast  ) = pprt (i,j,k,tpresent)
          qv   (i,j,k,tpast  ) = qv   (i,j,k,tpresent)
          qc   (i,j,k,tpast  ) = qc   (i,j,k,tpresent)
          qr   (i,j,k,tpast  ) = qr   (i,j,k,tpresent)
          qi   (i,j,k,tpast  ) = qi   (i,j,k,tpresent)
          qs   (i,j,k,tpast  ) = qs   (i,j,k,tpresent)
          qh   (i,j,k,tpast  ) = qh   (i,j,k,tpresent)
          tke  (i,j,k,tpast  ) = tke  (i,j,k,tpresent)
        END DO
      END DO
    END DO
  END IF
!
  DO k=1,nz
    DO j=1,ny
      DO i=1,nx
        ptcumsrc(i,j,k)=0.0
        qcumsrc(i,j,k,1)=0.0
        qcumsrc(i,j,k,2)=0.0
        qcumsrc(i,j,k,3)=0.0
        qcumsrc(i,j,k,4)=0.0
        qcumsrc(i,j,k,5)=0.0
      END DO
    END DO
  END DO

!-----------------------------------------------------------------------
!
!  Print out the max/min of initial arrays read in.
!
!-----------------------------------------------------------------------
!
  tim = tpresent

  IF (myproc ==0) &
    WRITE(6,'(1x,a/)') 'Min. and max. of the data arrays read in:'

  CALL a3dmax0(x,1,nx,1,nx,1,1,1,1, 1,1,1,1, amax,amin)
  IF (myproc ==0) &
    WRITE(6,'(1x,2(a,e13.6))') 'xmin    = ', amin,',  xmax    =',amax

  CALL a3dmax0(y,1,ny,1,ny,1,1,1,1, 1,1,1,1, amax,amin)
  IF (myproc ==0) &
    WRITE(6,'(1x,2(a,e13.6))') 'ymin    = ', amin,',  ymax    =',amax

  CALL a3dmax0(z,1,nz,1,nz,1,1,1,1, 1,1,1,1, amax,amin)
  IF (myproc ==0) &
    WRITE(6,'(1x,2(a,e13.6))') 'zmin    = ', amin,',  zmax    =',amax

  CALL a3dmax0(zp,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz,                    &
               amax,amin)
  IF (myproc ==0) &
    WRITE(6,'(1x,2(a,e13.6))') 'zpmin   = ', amin,',  zpmax   =',amax

  CALL a3dmax0(j3,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,                  &
               amax,amin)
  IF (myproc ==0) &
    WRITE(6,'(1x,2(a,e13.6))') 'j3min   = ', amin,',  j3max   =',amax

  CALL a3dmax0(ubar,1,nx,1,nx,1,ny,1,ny-1,1,nz,1,nz-1,                  &
               amax,amin)
  IF (myproc ==0) &
    WRITE(6,'(1x,2(a,e13.6))') 'ubarmin = ', amin,',  ubarmax =',amax

  CALL a3dmax0(vbar,1,nx,1,nx-1,1,ny,1,ny,1,nz,1,nz-1,                  &
               amax,amin)
  IF (myproc ==0) &
    WRITE(6,'(1x,2(a,e13.6))') 'vbarmin = ', amin,',  vbarmax =',amax

  CALL a3dmax0(ptbar,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,               &
               amax,amin)
  IF (myproc ==0) &
    WRITE(6,'(1x,2(a,e13.6))') 'ptbarmin= ', amin,',  ptbarmax=',amax

  CALL a3dmax0(pbar,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,                &
               amax,amin)
  IF (myproc ==0) &
    WRITE(6,'(1x,2(a,e13.6))') 'pbarmin = ', amin,',  pbarmax =',amax

  CALL a3dmax0(rhostr,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,              &
               amax,amin)
  IF (myproc ==0) &
    WRITE(6,'(1x,2(a,e13.6))') 'rhostrmin=', amin,', rhostrmax=',amax

  CALL a3dmax0(qvbar,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,               &
               amax,amin)
  IF (myproc ==0) &
    WRITE(6,'(1x,2(a,e13.6))') 'qvbarmin= ', amin,',  qvbarmax=',amax

  CALL a3dmax0(               amax,amin)
  IF (myproc ==0) &
    WRITE(6,'(1x,2(a,e13.6))') 'umin    = ', amin,',  umax    =',amax

  CALL a3dmax0(               amax,amin)
  IF (myproc ==0) &
    WRITE(6,'(1x,2(a,e13.6))') 'vmin    = ', amin,',  vmax    =',amax

  CALL a3dmax0(               amax,amin)
  IF (myproc ==0) &
    WRITE(6,'(1x,2(a,e13.6))') 'wmin    = ', amin,',  wmax    =',amax

  CALL a3dmax0(ptprt(1,1,1,tim),1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,    &
               amax,amin)
  IF (myproc ==0) &
    WRITE(6,'(1x,2(a,e13.6))') 'ptprtmin= ', amin,',  ptprtmax=',amax

  CALL a3dmax0(pprt(1,1,1,tim),1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,     &
               amax,amin)
  IF (myproc ==0) &
    WRITE(6,'(1x,2(a,e13.6))') 'pprtmin = ', amin,',  pprtmax =',amax

  CALL a3dmax0(qv(1,1,1,tim),1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,       &
               amax,amin)
  IF (myproc ==0) &
    WRITE(6,'(1x,2(a,e13.6))') 'qvmin   = ', amin,',  qvmax   =',amax

  CALL a3dmax0(qc(1,1,1,tim),1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,       &
               amax,amin)
  IF (myproc ==0) &
    WRITE(6,'(1x,2(a,e13.6))') 'qcmin   = ', amin,',  qcmax   =',amax

  CALL a3dmax0(qr(1,1,1,tim),1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,       &
               amax,amin)
  IF (myproc ==0) &
    WRITE(6,'(1x,2(a,e13.6))') 'qrmin   = ', amin,',  qrmax   =',amax

  CALL a3dmax0(qi(1,1,1,tim),1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,       &
               amax,amin)
  IF (myproc ==0) &
    WRITE(6,'(1x,2(a,e13.6))') 'qimin   = ', amin,',  qimax   =',amax

  CALL a3dmax0(qs(1,1,1,tim),1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,       &
               amax,amin)
  IF (myproc ==0) &
    WRITE(6,'(1x,2(a,e13.6))') 'qsmin   = ', amin,',  qsmax   =',amax

  CALL a3dmax0(qh(1,1,1,tim),1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,       &
               amax,amin)
  IF (myproc ==0) &
    WRITE(6,'(1x,2(a,e13.6))') 'qhmin   = ', amin,',  qhmax   =',amax

  CALL a3dmax0(tke(1,1,1,tim),1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,      &
               amax,amin)
  IF (myproc ==0) &
    WRITE(6,'(1x,2(a,e13.6))') 'tkemin  = ', amin,',  tkemax  =',amax

!  IF ( sfcphy.gt.0 ) THEN

!    CALL a3dmax0(tsfc,1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1, amax,amin)
!    write(6,'(1x,2(a,e13.6))') 'tsfcmin = ',amin,', tsfcmax  =',amax

!    CALL a3dmax0(tsoil,1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1, amax,amin)
!    write(6,'(1x,2(a,e13.6))') 'tsoilmin= ',amin,', tsoilmax =',amax

!    CALL a3dmax0(wetsfc,1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1, amax,amin)
!    write(6,'(1x,2(a,e13.6))') 'wetsmin = ',amin,', wetsmax  =',amax

!    CALL a3dmax0(wetdp,1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1, amax,amin)
!    write(6,'(1x,2(a,e13.6))') 'wetdmin = ',amin,', wetdmax  =',amax

!    CALL a3dmax0(wetcanp,1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1, amax,amin)
!    write(6,'(1x,2(a,e13.6))') 'wetcmin = ',amin,', wetcmax  =',amax
!  ENDIF

  CALL a3dmax0(ptcumsrc,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,            &
               amax,amin)
  IF (myproc ==0) &
    WRITE(6,'(1x,2(a,e13.6))') 'ptcummin= ', amin,',  ptcummax=',amax

  CALL a3dmax0(qcumsrc(1,1,1,1),1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,    &
               amax,amin)
  IF (myproc ==0) &
    WRITE(6,'(1x,2(a,e13.6))') 'qvcummin= ', amin,',  qvcummax=',amax
  CALL a3dmax0(qcumsrc(1,1,1,2),1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,    &
               amax,amin)
  IF (myproc ==0) &
    WRITE(6,'(1x,2(a,e13.6))') 'qccummin= ', amin,',  qccummax=',amax
  CALL a3dmax0(qcumsrc(1,1,1,3),1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,    &
               amax,amin)
  IF (myproc ==0) &
    WRITE(6,'(1x,2(a,e13.6))') 'qrcummin= ', amin,',  qrcummax=',amax
  CALL a3dmax0(qcumsrc(1,1,1,4),1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,    &
               amax,amin)
  IF (myproc ==0) &
    WRITE(6,'(1x,2(a,e13.6))') 'qicummin= ', amin,',  qicummax=',amax
  CALL a3dmax0(qcumsrc(1,1,1,5),1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,    &
               amax,amin)
  IF (myproc ==0) &
    WRITE(6,'(1x,2(a,e13.6))') 'qscummin= ', amin,',  qscummax=',amax

  RETURN

END SUBROUTINE extinit
!
!
!##################################################################
!##################################################################
!######                                                      ######
!######                 SUBROUTINE INITSFC                   ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!
!


SUBROUTINE initsfc(nx,ny,nz,nstyps,                                     & 1,79
           pbar,pprt,ptbar,ptprt,qvbar,qv,                              &
           soiltyp,stypfrct,vegtyp,lai,roufns,veg,                      &
           ndvi,                                                        &
           tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth,qvsfc)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Initialize the surface characteristics data and soil model
!  variables according to option parameters sfcdat and soilinit.
!
!  The surface and soil data files are sequential binary files.
!
!  Their structures should be:
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Yuhe Liu
!  02/17/94
!
!  MODIFICATION:
!
!  02/07/1995 (Yuhe Liu)
!  Added a new 2-D permanent array, veg(nx,ny), to the soil model and
!  at the same time delete the table data array veg(14).
!
!  01/29/1995 (V. Wong and X. Song)
!  Add a flag wtrexist in "initsfc".
!
!  12/04/1997 (Yuhe Liu)
!  Set the soil variables to their default value read in from input
!  file if they are not included in soilinit file.
!
!  12/09/1998 (Donghai Wang)
!  Added the snow cover.
!
!  2000/01/07 (Gene Bassett)
!  Modified the behavior for the case where soil data read in from file
!  and also present in history file.  If soilinit=2, initopt=3 and
!  sfcin=1, for soil variables that are not in soildata file the values
!  in the history file are used.
!
!-----------------------------------------------------------------------
!
!  INPUT and OUTPUT:
!
!    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 z-direction (sfc/top)
!
!    pbar     Base state pressure
!    pprt     Preturbation pressure
!
!    ptbar    Base state potential temperature
!    ptprt    Preturbation potential temperature
!
!    qvbar    Base state water vapor mixing ratio
!    qv       Water vapor mixing ratio
!
!    soiltyp  Soil type at the horizontal grid points
!    stypfrct  Soil type fraction
!    vegtyp   Vegetation type at the horizontal grid points
!    lai      Leaf Area Index
!    roufns   Surface roughness
!    veg      Vegetation fraction
!
!    tsfc     Temperature at ground (K) (in top 1 cm layer)
!    tsoil    Deep soil temperature (K) (in deep 1 m layer)
!    wetsfc   Surface soil moisture in the top 1 cm layer
!    wetdp    Deep soil moisture in the deep 1 m layer
!    wetcanp  Canopy water amount
!    qvsfc    Effective S.H. at sfc.
!
!  TEMPORARY WORKING ARRAY
!
!
!-----------------------------------------------------------------------
!
!  Variable Declarations. (Local Variables)
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!
  REAL :: pi
  PARAMETER ( pi = 3.141592654)

  INTEGER :: nx, ny, nz

  REAL :: pbar   (nx,ny,nz) ! Base state pressure (Pascal)
  REAL :: pprt   (nx,ny,nz) ! Perturbation pressure (Pascal)
  REAL :: ptbar  (nx,ny,nx) ! Base state potential temperature (K)
  REAL :: ptprt  (nx,ny,nz) ! Perturbation potential temperature (K)
  REAL :: qvbar  (nx,ny,nz) ! Base state specific humidity (kg/kg)
  REAL :: qv     (nx,ny,nz) ! Total specific humidity (kg/kg)

  INTEGER :: nstyps                  ! Number of soil types
  INTEGER :: soiltyp(nx,ny,nstyps)   ! Soil type at each point
  REAL :: stypfrct(nx,ny,nstyps)  ! Fraction of soil types
  INTEGER :: vegtyp (nx,ny)          ! Vegetation type at each point
  REAL :: lai    (nx,ny)          ! Leaf Area Index
  REAL :: roufns (nx,ny)          ! Surface roughness
  REAL :: veg    (nx,ny)          ! Vegetation fraction
  REAL :: ndvi   (nx,ny)          ! NDVI

  REAL :: tsfc   (nx,ny,0:nstyps) ! Ground surface temperature (K)
  REAL :: tsoil  (nx,ny,0:nstyps) ! Deep soil temperature (K)
  REAL :: wetsfc (nx,ny,0:nstyps) ! Ground surface soil moisture
  REAL :: wetdp  (nx,ny,0:nstyps) ! Deep soil moisture
  REAL :: wetcanp(nx,ny,0:nstyps) ! Canopy moisture
  REAL :: snowdpth(nx,ny)         ! Snow depth (m)
  REAL :: qvsfc  (nx,ny,0:nstyps) ! Effective surface specific
                                  ! humidity (kg/kg)

!
!-----------------------------------------------------------------------
!
!  Define local variables
!
!-----------------------------------------------------------------------
!
  REAL :: psfc        ! Surface pressure
  REAL :: rhgs        ! The relative humidity at ground surface
  REAL :: qvsatts     ! Saturated specific humidity at surface

  REAL :: wrmax       ! Maximum value for canopy moisture, wetcanp

  REAL :: pterm       ! Temporary variable, real positive term flag

  REAL :: p0inv       ! Inverse of p0 (1000 mb)

  REAL :: tema        ! Temporary variable
  REAL :: temb        ! Temporary variable
!
!-----------------------------------------------------------------------
!
!  Global constants and parameters, most of them specify the
!  model run options.
!
!-----------------------------------------------------------------------
!
  INCLUDE 'globcst.inc'
  INCLUDE 'grid.inc'          ! Grid & map parameters.
  INCLUDE 'indtflg.inc'
  INCLUDE 'phycst.inc'
  INCLUDE 'soilcst.inc'
  INCLUDE 'bndry.inc'
  INCLUDE 'mp.inc'            ! Message passing parameters.
!
!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------
!
  INTEGER :: i,j,is,imid,jmid
  INTEGER :: tsfcin
  INTEGER :: tsoilin
  INTEGER :: wsfcin
  INTEGER :: wdpin
  INTEGER :: wcanpin
  INTEGER :: snowdin

  INTEGER :: mptag, astat
  INTEGER, allocatable :: mpitmp(:)
  REAL, allocatable :: mprtmp(:)
!
!-----------------------------------------------------------------------
!
!  Function f_qvsat and inline directive for Cray PVP
!
!-----------------------------------------------------------------------
!
  REAL :: f_qvsat

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

!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  IF (mp_opt > 0 .AND. (sfcdat == 1 .OR. soilinit == 1)) THEN
    ALLOCATE(mpitmp(MAX(nx,ny)),stat=astat)
    IF (astat /= 0) THEN
      WRITE (6,*)                                                       &
          "INITSFC: ERROR allocating temporary array mpitmp."
      CALL arpsstop ("arpsstop called from initsfc allocating mpitmp",1)
    END IF
    ALLOCATE(mprtmp(MAX(nx,ny)),stat=astat)
    IF (astat /= 0) THEN
      WRITE (6,*)                                                       &
          "INITSFC: ERROR allocating temporary array mprtmp."
      CALL arpsstop ("arpsstop called from initsfc allocating mprtmp",1)
    END IF
  END IF

  p0inv = 1.0/p0

  IF ( initopt == 2 ) THEN
    WRITE (6, '(a,a/a/a/a/a)')                                          &
        'Model initialization option was set to restart. ',             &
        'All the surface and soil arrays should have been read',        &
        'in from the restart file. No more initialization is ',         &
        'done in subroutine INITSFC.'

    wtrexist=0
    DO j = 1, ny
      DO i = 1, nx
        DO is = 1,nstyp
          IF (soiltyp(i,j,is) == 13) THEN
            wtrexist=1
            RETURN
          END IF
        END DO
      END DO
    END DO

    RETURN

  END IF

  IF ( sfcdat == 1 ) THEN

    nstyp = 1

    DO j=1, ny-1
      DO i=1, nx-1
        soiltyp(i,j,1) = styp
        vegtyp (i,j) = vtyp
        lai    (i,j) = lai0
        roufns (i,j) = roufns0
        veg    (i,j) = veg0
      END DO
    END DO

  ELSE IF (sfcdat == 2 .OR. (sfcdat == 3.AND.landin /= 1) ) THEN
!
!-----------------------------------------------------------------------
!
!  Read the surface property data from file sfcdtfl (passed
!  in globcst.inc).
!
!-----------------------------------------------------------------------
!
!  blocking inserted for ordering i/o for message passing
    DO i=0,nprocs-1,max_fopen
      IF(myproc >= i.AND.myproc <= i+max_fopen-1)THEN

        CALL readsfcdt( nx,ny,nstyps,sfcdtfl,dx,dy,                     &
               mapproj,trulat1,trulat2,trulon,sclfct,ctrlat,ctrlon,     &
                      soiltyp,stypfrct,vegtyp,lai,roufns,veg,ndvi )

      END IF
      IF (mp_opt > 0) CALL mpbarrier
    END DO

  ELSE IF(sfcdat == 3 .AND. landin == 1) THEN

    WRITE(6,'(1x,a/,a/)')                                               &
        'Surface property data in the initialization file was used.',   &
        'No more initialization performed.'

  ELSE

    WRITE(6,'(1x,a,i3,a/)')                                             &
        'Invalid surface data input option. sfcdat =',sfcdat,           &
        '. Program stopped in INITSFC.'
      CALL arpsstop ("arpsstop called from initsfc with sfcdat option",1)

  END IF

  print *, ' nstyps = ',nstyps
  imid=(nx/2)+1
  jmid=(ny/2)+1
  DO is=1,nstyps
    print *, '  Sample soil ( ',imid,jmid,') = ',                    &
             soiltyp(imid,jmid,is),stypfrct(imid,jmid,is)
  END DO

  IF ( nstyp == 1 ) THEN
    DO j=1,ny
      DO i=1,nx
        stypfrct(i,j,1) = 1.0
      END DO
    END DO
  END IF

  IF ( soilinit == 1 ) THEN
!
!-----------------------------------------------------------------------
!
!  All surface variables are specified uniformly in horizontal from
!  the input namelist.
!
!-----------------------------------------------------------------------
    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)) )

        IF( soiltyp(i,j,1) == 13) THEN  ! For water
          tsfc(i,j,0) = ptswtr0*(psfc*p0inv)**rddcp
        ELSE                          ! For land
          tsfc(i,j,0) = ptslnd0*(psfc*p0inv)**rddcp
        END IF

        tsoil  (i,j,0) = tsoil0
        wetsfc (i,j,0) = wetsfc0
        wetdp  (i,j,0) = wetdp0
        wetcanp(i,j,0) = wetcanp0
        snowdpth(i,j)  = snowdpth0
      END DO
    END DO

    DO is=1,nstyp
      DO j=1,ny-1
        DO i=1,nx-1
          tsfc   (i,j,is) = tsfc   (i,j,0)
          tsoil  (i,j,is) = tsoil  (i,j,0)
          wetsfc (i,j,is) = wetsfc (i,j,0)
          wetdp  (i,j,is) = wetdp  (i,j,0)
          wetcanp(i,j,is) = wetcanp(i,j,0)
        END DO
      END DO
    END DO

  ELSE IF ((soilinit == 2).OR.                                          &
            (soilinit == 3.AND.sfcin /= 1) ) THEN
!
!-----------------------------------------------------------------------
!
!  Read the soil variables from file soilinfl. soilinfl is
!  passed in globcst.inc.
!
!-----------------------------------------------------------------------
!
!  blocking inserted for ordering i/o for message passing
    DO i=0,nprocs-1,max_fopen
      IF(myproc >= i.AND.myproc <= i+max_fopen-1)THEN

        CALL readsoil(nx,ny,nstyps,soilinfl,dx,dy,                      &
               mapproj,trulat1,trulat2,trulon,sclfct,ctrlat,ctrlon,     &
                    tsfcin,tsoilin,wsfcin,wdpin,wcanpin,snowdin,        &
                    tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth,soiltyp)

      END IF
      IF (mp_opt > 0) CALL mpbarrier
    END DO

    IF (sfcin /= 1) THEN
      IF ( tsfcin /= 1 ) THEN
        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)) )

              IF( soiltyp(i,j,1) == 13) THEN  ! For water
                tsfc(i,j,is) = ptswtr0*(psfc*p0inv)**rddcp
              ELSE                          ! For land
                tsfc(i,j,is) = ptslnd0*(psfc*p0inv)**rddcp
              END IF
            END DO
          END DO
        END DO
      END IF

      IF ( tsoilin /= 1 ) THEN
        DO is=1,nstyp
          DO j=1, ny-1
            DO i=1, nx-1
              tsoil(i,j,is) = tsoil0
            END DO
          END DO
        END DO
      END IF

      IF ( wsfcin /= 1 ) THEN
        DO is=1,nstyp
          DO j=1, ny-1
            DO i=1, nx-1
              wetsfc(i,j,is) = wetsfc0
            END DO
          END DO
        END DO
      END IF
      IF ( wdpin /= 1 ) THEN
        DO is=1,nstyp
          DO j=1, ny-1
            DO i=1, nx-1
              wetdp(i,j,is) = wetdp0
            END DO
          END DO
        END DO
      END IF
      IF ( wcanpin /= 1 ) THEN
        DO is=1,nstyp
          DO j=1, ny-1
            DO i=1, nx-1
              wetcanp(i,j,is) = wetcanp0
            END DO
          END DO
        END DO
      END IF
      IF ( snowdin /= 1 ) THEN
        DO j=1, ny-1
          DO i=1, nx-1
            snowdpth(i,j) = snowdpth0
          END DO
        END DO
      END IF
    END IF

  ELSE IF(soilinit == 3.AND.sfcin == 1) THEN

    WRITE(6,'(1x,a/,a/)')                                               &
        'Surface variable in the initialization file was used.',        &
        'No more initialization performed.'

    IF ( nstyp == 1 ) THEN
      DO j=1,ny-1
        DO i=1,nx-1
          tsfc   (i,j,1) = tsfc   (i,j,0)
          tsoil  (i,j,1) = tsoil  (i,j,0)
          wetsfc (i,j,1) = wetsfc (i,j,0)
          wetdp  (i,j,1) = wetdp  (i,j,0)
          wetcanp(i,j,1) = wetcanp(i,j,0)
        END DO
      END DO
    END IF
  ELSE IF(soilinit == 4) THEN
    p0inv=1./p0
    DO j=1,ny-1
      DO i=1,nx-1
        tema=0.5*((ptbar(i,j,2)+ptprt(i,j,2))                           &
                 +(ptbar(i,j,1)+ptprt(i,j,1)))
        psfc=0.5*((pbar(i,j,2)+pprt(i,j,2))                             &
                 +(pbar(i,j,1)+pprt(i,j,1)))

        temb = tema*(psfc*p0inv)**rddcp

        tsfc   (i,j,0) = temb + tsprt
        tsoil  (i,j,0) = temb + t2prt
        wetsfc (i,j,0) = 0.0
        wetdp  (i,j,0) = 0.0
        wetcanp(i,j,0) = wetcanp0
        snowdpth(i,j)  = snowdpth0
      END DO
    END DO

    DO is=1,nstyp
      DO j=1,ny-1
        DO i=1,nx-1
          tsfc   (i,j,is) = tsfc (i,j,0)
          tsoil  (i,j,is) = tsoil(i,j,0)
          wetsfc (i,j,is) = wgrat*wsat(soiltyp(i,j,is))
          wetdp  (i,j,is) = w2rat*wsat(soiltyp(i,j,is))
          wetcanp(i,j,is) = wetcanp(i,j,0)

          wetsfc(i,j,0 ) = wetsfc(i,j,0)+stypfrct(i,j,is)*wetsfc(i,j,is)
          wetdp (i,j,0 ) = wetdp(i,j,0)+stypfrct(i,j,is)*wetdp(i,j,is)
        END DO
      END DO
    END DO

  ELSE

    WRITE(6,'(1x,a,i3,a/)')                                             &
        'Invalid surface variable input option. soilinit =',soilinit,   &
        '. Program stopped in INITSFC.'
      CALL arpsstop ("arpsstop called from initsfc with soilinit option",1)

  END IF

  IF ( moist /= 0 ) THEN
    DO is=1,nstyp
      DO j = 1, ny-1
        DO i = 1, nx-1
          wrmax = 0.2*1.0E-3*veg(i,j)*lai(i,j)

          wetsfc (i,j,is) = MAX( wetsfc(i,j,is), 0.0 )
          wetsfc (i,j,is) = MIN( wetsfc(i,j,is), wsat(soiltyp(i,j,is)) )
          wetdp  (i,j,is) = MAX( wetdp(i,j,is), 0.0 )
          wetdp  (i,j,is) = MIN( wetdp(i,j,is), wsat(soiltyp(i,j,is)) )
          wetcanp(i,j,is) = MAX( wetcanp(i,j,is), 0.0 )
          wetcanp(i,j,is) = MIN( wetcanp(i,j,is), wrmax )
        END DO
      END DO
    END DO
    DO j = 1, ny-1
      DO i = 1, nx-1
        wetsfc (i,j,0) = MAX( wetsfc(i,j,0), 0.0 )
        wetdp  (i,j,0) = MAX( wetdp(i,j,0), 0.0 )
        wetcanp(i,j,0) = MAX( wetcanp(i,j,0), 0.0 )
      END DO
    END DO
  ELSE
    DO is=0,nstyp
      DO j = 1, ny-1
        DO i = 1, nx-1
          wetsfc(i,j,is)  = 0.0
          wetdp(i,j,is)   = 0.0
          wetcanp(i,j,is) = 0.0
        END DO
      END DO
    END DO
  END IF

  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 IF ( soiltyp(i,j,is) > 0 ) THEN

          pterm = .5+SIGN(.5,wetsfc(i,j,is)-1.1*wfc(soiltyp(i,j,is)))

          rhgs = pterm + (1.-pterm)                                     &
                       * 0.25 * ( 1.-COS(wetsfc(i,j,is)*pi              &
                                    /(1.1*wfc(soiltyp(i,j,is)))) )**2
        ELSE

          rhgs = 0.0

        END IF

        qvsfc(i,j,is) = rhgs * qvsatts

      END DO
    END DO
  END DO

!
!-----------------------------------------------------------------------
!
!  Ensure soil values are consistent with each other.
!
!-----------------------------------------------------------------------
!
  DO j=1,ny-1
    DO i=1,nx-1
      tsfc   (i,j,0) =  0.0
      tsoil  (i,j,0) =  0.0
      wetsfc (i,j,0) =  0.0
      wetdp  (i,j,0) =  0.0
      wetcanp(i,j,0) =  0.0
    END DO
  END DO

  DO is = 1,nstyp
    DO j=1,ny-1
      DO i=1,nx-1
        tsfc   (i,j,0) = tsfc   (i,j,0)                                 &
                       + tsfc   (i,j,is) * stypfrct(i,j,is)
        tsoil  (i,j,0) = tsoil  (i,j,0)                                 &
                       + tsoil  (i,j,is) * stypfrct(i,j,is)
        wetsfc (i,j,0) = wetsfc (i,j,0)                                 &
                       + wetsfc (i,j,is) * stypfrct(i,j,is)
        wetdp  (i,j,0) = wetdp  (i,j,0)                                 &
                       + wetdp  (i,j,is) * stypfrct(i,j,is)
        wetcanp(i,j,0) = wetcanp(i,j,0)                                 &
                       + wetcanp(i,j,is) * stypfrct(i,j,is)
      END DO
    END DO
  END DO

  DO j=1,ny-1
    DO i=1,nx-1
      IF ((vegtyp(i,j) == 14) .OR. (tsfc(i,j,0) > 273.16)) THEN
        snowdpth(i,j) = 0.  ! Make sure that snow cover is consistent
      END IF                 ! with vegetation type and soil temperature.
    END DO
  END DO

  DO j = 1, ny-1
    DO i = 1, nx-1
      qvsfc(i,j,0) = 0.0
    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

  IF ( sfcdat == 1 ) THEN
    DO is=1,nstyp
      IF (mp_opt > 0) THEN
        CALL acct_interrupt(mp_acct)
        CALL mpsend1diew(soiltyp(1,1,is),nx,ny,ebc,wbc,0,mptag,mpitmp)
        CALL mprecv1diew(soiltyp(1,1,is),nx,ny,ebc,wbc,0,mptag,mpitmp)
        CALL mpsend1dins(soiltyp(1,1,is),nx,ny,nbc,sbc,0,mptag,mpitmp)
        CALL mprecv1dins(soiltyp(1,1,is),nx,ny,nbc,sbc,0,mptag,mpitmp)
        CALL mpsend1dew(stypfrct(1,1,is),nx,ny,ebc,wbc,0,mptag,mprtmp)
        CALL mprecv1dew(stypfrct(1,1,is),nx,ny,ebc,wbc,0,mptag,mprtmp)
        CALL mpsend1dns(stypfrct(1,1,is),nx,ny,nbc,sbc,0,mptag,mprtmp)
        CALL mprecv1dns(stypfrct(1,1,is),nx,ny,nbc,sbc,0,mptag,mprtmp)
      END IF
      CALL acct_interrupt(bc_acct)
      CALL bcis2d(nx,ny, soiltyp (1,1,is), ebc,wbc,nbc,sbc)
      CALL bcs2d (nx,ny, stypfrct(1,1,is), ebc,wbc,nbc,sbc)
      CALL acct_stop_inter
    END DO

    IF (mp_opt > 0) THEN
      CALL acct_interrupt(mp_acct)
      CALL mpsend1diew(vegtyp,nx,ny,ebc,wbc,0,mptag,mpitmp)
      CALL mprecv1diew(vegtyp,nx,ny,ebc,wbc,0,mptag,mpitmp)
      CALL mpsend1dins(vegtyp,nx,ny,nbc,sbc,0,mptag,mpitmp)
      CALL mprecv1dins(vegtyp,nx,ny,nbc,sbc,0,mptag,mpitmp)
      CALL mpsend1dew(lai,nx,ny,ebc,wbc,0,mptag,mprtmp)
      CALL mprecv1dew(lai,nx,ny,ebc,wbc,0,mptag,mprtmp)
      CALL mpsend1dns(lai,nx,ny,nbc,sbc,0,mptag,mprtmp)
      CALL mprecv1dns(lai,nx,ny,nbc,sbc,0,mptag,mprtmp)
      CALL mpsend1dew(roufns,nx,ny,ebc,wbc,0,mptag,mprtmp)
      CALL mprecv1dew(roufns,nx,ny,ebc,wbc,0,mptag,mprtmp)
      CALL mpsend1dns(roufns,nx,ny,nbc,sbc,0,mptag,mprtmp)
      CALL mprecv1dns(roufns,nx,ny,nbc,sbc,0,mptag,mprtmp)
      CALL mpsend1dew(veg,nx,ny,ebc,wbc,0,mptag,mprtmp)
      CALL mprecv1dew(veg,nx,ny,ebc,wbc,0,mptag,mprtmp)
      CALL mpsend1dns(veg,nx,ny,nbc,sbc,0,mptag,mprtmp)
      CALL mprecv1dns(veg,nx,ny,nbc,sbc,0,mptag,mprtmp)
      CALL mpsend1dew(snowdpth,nx,ny,ebc,wbc,0,mptag,mprtmp)
      CALL mprecv1dew(snowdpth,nx,ny,ebc,wbc,0,mptag,mprtmp)
      CALL mpsend1dns(snowdpth,nx,ny,nbc,sbc,0,mptag,mprtmp)
      CALL mprecv1dns(snowdpth,nx,ny,nbc,sbc,0,mptag,mprtmp)
    END IF
    CALL acct_interrupt(bc_acct)
    CALL bcis2d(nx,ny, vegtyp, ebc,wbc,nbc,sbc)
    CALL bcs2d (nx,ny, lai,    ebc,wbc,nbc,sbc)
    CALL bcs2d (nx,ny, roufns, ebc,wbc,nbc,sbc)
    CALL bcs2d (nx,ny, veg,    ebc,wbc,nbc,sbc)
    CALL bcs2d (nx,ny, snowdpth,ebc,wbc,nbc,sbc)
    CALL acct_stop_inter
  END IF
!call test_dump (roufns,"XXXAinit_roufns",nx,ny,1,0,1)

  IF ( soilinit == 1 ) THEN
    DO is=0,nstyp
      IF (mp_opt > 0) THEN
        CALL acct_interrupt(mp_acct)
        CALL mpsend1dew(tsfc(1,1,is),nx,ny,ebc,wbc,0,mptag,mprtmp)
        CALL mprecv1dew(tsfc(1,1,is),nx,ny,ebc,wbc,0,mptag,mprtmp)
        CALL mpsend1dns(tsfc(1,1,is),nx,ny,nbc,sbc,0,mptag,mprtmp)
        CALL mprecv1dns(tsfc(1,1,is),nx,ny,nbc,sbc,0,mptag,mprtmp)
        CALL mpsend1dew(tsoil(1,1,is),nx,ny,ebc,wbc,0,mptag,mprtmp)
        CALL mprecv1dew(tsoil(1,1,is),nx,ny,ebc,wbc,0,mptag,mprtmp)
        CALL mpsend1dns(tsoil(1,1,is),nx,ny,nbc,sbc,0,mptag,mprtmp)
        CALL mprecv1dns(tsoil(1,1,is),nx,ny,nbc,sbc,0,mptag,mprtmp)
        CALL mpsend1dew(wetsfc(1,1,is),nx,ny,ebc,wbc,0,mptag,mprtmp)
        CALL mprecv1dew(wetsfc(1,1,is),nx,ny,ebc,wbc,0,mptag,mprtmp)
        CALL mpsend1dns(wetsfc(1,1,is),nx,ny,nbc,sbc,0,mptag,mprtmp)
        CALL mprecv1dns(wetsfc(1,1,is),nx,ny,nbc,sbc,0,mptag,mprtmp)
        CALL mpsend1dew(wetdp(1,1,is),nx,ny,ebc,wbc,0,mptag,mprtmp)
        CALL mprecv1dew(wetdp(1,1,is),nx,ny,ebc,wbc,0,mptag,mprtmp)
        CALL mpsend1dns(wetdp(1,1,is),nx,ny,nbc,sbc,0,mptag,mprtmp)
        CALL mprecv1dns(wetdp(1,1,is),nx,ny,nbc,sbc,0,mptag,mprtmp)
        CALL mpsend1dew(wetcanp(1,1,is),nx,ny,ebc,wbc,0,mptag,mprtmp)
        CALL mprecv1dew(wetcanp(1,1,is),nx,ny,ebc,wbc,0,mptag,mprtmp)
        CALL mpsend1dns(wetcanp(1,1,is),nx,ny,nbc,sbc,0,mptag,mprtmp)
        CALL mprecv1dns(wetcanp(1,1,is),nx,ny,nbc,sbc,0,mptag,mprtmp)
        CALL mpsend1dew(qvsfc(1,1,is),nx,ny,ebc,wbc,0,mptag,mprtmp)
        CALL mprecv1dew(qvsfc(1,1,is),nx,ny,ebc,wbc,0,mptag,mprtmp)
        CALL mpsend1dns(qvsfc(1,1,is),nx,ny,nbc,sbc,0,mptag,mprtmp)
        CALL mprecv1dns(qvsfc(1,1,is),nx,ny,nbc,sbc,0,mptag,mprtmp)
      END IF
      CALL acct_interrupt(bc_acct)
      CALL bcs2d (nx,ny, tsfc   (1,1,is), ebc,wbc,nbc,sbc)
      CALL bcs2d (nx,ny, tsoil  (1,1,is), ebc,wbc,nbc,sbc)
      CALL bcs2d (nx,ny, wetsfc (1,1,is), ebc,wbc,nbc,sbc)
      CALL bcs2d (nx,ny, wetdp  (1,1,is), ebc,wbc,nbc,sbc)
      CALL bcs2d (nx,ny, wetcanp(1,1,is), ebc,wbc,nbc,sbc)
      CALL bcs2d (nx,ny, qvsfc  (1,1,is), ebc,wbc,nbc,sbc)
      CALL acct_stop_inter
    END DO
  END IF

  IF (mp_opt > 0 .AND. (sfcdat == 1 .OR. soilinit == 1)) THEN
    DEALLOCATE (mpitmp,stat=astat)
    DEALLOCATE (mprtmp,stat=astat)
  END IF

  wtrexist=0
  DO j = 1, ny
    DO i = 1, nx
      DO is = 1,nstyp
        IF (soiltyp(i,j,is) == 13) THEN
          wtrexist=1
          RETURN
        END IF
      END DO
    END DO
  END DO

  RETURN

  WRITE (6,'(/a,i2/a)')                                           &
         '     Read error in surface data file '                        &
         //soilinfl//' with the I/O unit ',sfcunit,                     &
         'The model will STOP in subroutine INITSFC.'

  CALL arpsstop ("arpsstop called from initsfc reading sfc data file",1)

  WRITE (6,'(/a,i2/a)')                                           &
         '     Read error in surface initial data file '                &
         //soilinfl//' with the I/O unit ',sfcunit,                     &
         'The model will STOP in subroutine INITSFC.'

  CALL arpsstop ("arpsstop called from initsfc reading sfc data file-2",1)

END SUBROUTINE initsfc
!wdt update - init0 no longer necessary (arrays set to 0 after allocation)
!
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE FLZERO                     ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE flzero(a,n) 103

!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Fill in an entire array with zeros.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Ming Xue
!  10/10/1991.
!
!  MODIFICATION HISTORY:
!
!  5/02/92 (M. Xue)
!  Added full documentation.
!
!  5/03/92 (M. Xue)
!  Further documentation.
!
!  9/10/94 (Weygandt & Y. Lu)
!  Cleaned up documentation.
!
!-----------------------------------------------------------------------
!
!  INPUT :
!
!    n        The dimension of 1-D array a.
!
!  OUTPUT:
!
!    a        1-D array filled with zeros.
!
!-----------------------------------------------------------------------
!

!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!
!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------
!
  INTEGER :: i, n

  REAL :: a(n)                 ! 1-D array filled with zeros
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!

  DO i=1,n
    a(i)=0.0
  END DO

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


SUBROUTINE initlktb 2
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Initializes arrays used for lookup table functions.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Adwait Sathye
!  06/02/94
!
!  MODIFICATION HISTORY:
!
!  2/17/97 (J. Zong)
!  Corrected definitions of the power functions in loop 100.
!
!  2/24/97 (Jinxing Zong, Ming Xue and Yuhe Liu)
!  Defined five pwr arrays for lookup tables to replace fractional
!  power calculations in Tao ice microphysics.
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!  OUTPUT:
!     pwr arrays filled with evenly spaced data points
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!
!-----------------------------------------------------------------------
!
!  INCLUDE FILES
!
!-----------------------------------------------------------------------
!
  INCLUDE 'phycst.inc'
!
!-----------------------------------------------------------------------
!
!  LOCAL VARIABLES
!
!-----------------------------------------------------------------------
!
  INTEGER :: i
  REAL :: temper
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!

  DO i = 0, 1001

    pwr2046(i) =  ( 0.05 * i / 1000.0 ) ** 0.2046
    pwr525(i)  =  ( 0.05 * i / 1000.0 ) ** 0.525
    pwr875(i)  =  ( 0.05 * i / 1000.0 ) ** 0.875
    pwr1364(i) =  ( 0.05 * i / 1000000.0 ) ** 0.1364

  END DO

!-----------------------------------------------------------------------
!
!  Build a look up table for the Latent heat of vaporation, calculated
!  in the temperature range of (-100,50) degree Celcius
!  according to formula
!
!  Lv = lathv * (273.15/tbar)**(0.167+3.67E-4*tbar)
!
!  Lv for t(K) is given by
!
!  index = max( min(151, NINT(t-172.15)), 1)
!  Lv = latheatv(index)
!
!  where lathv = 2.500780e6 is set in phycst.inc, and t is the
!
!-----------------------------------------------------------------------

  DO i=-100,50
    temper =  i+273.15
    latheatv(i+101)=lathv * (273.15/temper)**(0.167+3.67E-4*temper)
  END DO

!-----------------------------------------------------------------------
!
!  Build lookup tables for T**0.81 and (rhobar/mu/psi**2)**one6th,
!  where mu and psi stand for dynamic viscosity of air and diffusivity
!  of water vapor in air, respectively. The temerature ranges from
!  -100 to 50 degree Celcius, pressure from 0 to 1500 hPa, and air
!  density from 0 to 1.5 kg/m**3 in the calculations of the tables.
!  With the above ranges of T, p and air density, (rhobar/mu/psi**2)
!  is between 0.0 and 3000.0.
!
!-----------------------------------------------------------------------

  DO i = 0, 151
    pwr81(i) = ( 173.15 + i ) ** 0.81
  END DO

  DO i = 0, 10001
    pwr1666(i) = ( 3000.0 * i / 10000.0 ) ** 0.1666667
  END DO


!-----------------------------------------------------------------------
!
!  Lookup tables for (rhobar*q)**a, where q can be qr, qs or qh.
!  rhobar is in g/cm***3. rhobar*q ranges from 0 to 50.e-6. The
!  increment of the tables is 50.e-9.
!
!-----------------------------------------------------------------------

  DO i = 0, 10001

    pwr2    (i) =  ( 0.05 * i / 10000000.0 ) ** 0.2
    pwr0625 (i) =  ( 0.05 * i / 10000000.0 ) ** 0.0625
    pwr15625(i) =  ( 0.05 * i / 10000000.0 ) ** 0.15625

  END DO

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


SUBROUTINE gtsinlat(nx,ny, x,y, sinlat, xs, ys, temxy) 1,1
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Calculate the sin of the lattitude of each grid point, to be used
!  in the calculation of latitude-dependent Coriolis parameters.
!
!-----------------------------------------------------------------------
!
!
!  AUTHOR: Ming Xue
!  3/21/95.
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    nx       Number of grid points in the x-direction (east/west)
!    ny       Number of grid points in the y-direction (north/south)
!    x        x-coordinate of grid points in computational space (m)
!    y        y-coordinate of grid points in computational space (m)
!
!  OUTPUT:
!
!    sinlat   Sin of latitude at each grid point
!
!  WORK ARRAYS:
!
!    xs       x-coordinate at scalar points.
!    ys       y-coordinate at scalar points.
!
!-----------------------------------------------------------------------
!

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

  INTEGER :: nx,ny          ! Dimensions of model grid.
  REAL :: x     (nx)        ! The x-coord. of the u points.
  REAL :: y     (ny)        ! The y-coord. of the v points.

  REAL :: sinlat(nx,ny)     ! Sin of latitude at each grid point

  REAL :: xs    (nx)        ! x-coord. of scalar points.
  REAL :: ys    (ny)        ! y-coord. of scalar points.

  REAL :: temxy (nx,ny)     ! Work array.

  REAL :: r2d
  INTEGER :: i,j
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  DO i=1,nx-1
    xs(i) = (x(i)+x(i+1))*0.5
  END DO
  xs(nx) = -xs(nx-2)+2.0*xs(nx-1)

  DO j=1,ny-1
    ys(j) = (y(j)+y(j+1))*0.5
  END DO
  ys(ny) = -ys(ny-2)+2.0*ys(ny-1)
!
  CALL xytoll(nx,ny,xs,ys,sinlat,temxy)

  r2d = ATAN(1.0)*4.0/180.0

  DO j=1,ny
    DO i=1,nx
      sinlat(i,j) = SIN( r2d * sinlat(i,j) )
    END DO
  END DO

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


SUBROUTINE setppi(nx,ny,nz,nt,tlvl,pprt,pbar,ppi) 2
!
  IMPLICIT NONE
  INTEGER :: nx,ny,nz,nt,tlvl
  REAL :: pprt(nx,ny,nz,nt)
  REAL :: pbar(nx,ny,nz)
  REAL :: ppi (nx,ny,nz)
!
  INCLUDE 'phycst.inc'
!
  INTEGER :: i,j,k
  REAL :: p0inv
!
  p0inv=1./p0
  DO k=1,nz-1
    DO j=1,ny-1
      DO i=1,nx-1
        ppi(i,j,k)=((pprt(i,j,k,tlvl)+pbar(i,j,k))*p0inv)**rddcp
      END DO
    END DO
  END DO
!
  RETURN
END SUBROUTINE setppi