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

SUBROUTINE rstout(nx,ny,nz, nstyps,exbcbufsz,                           & 2,28
           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,                                    &
           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)
!EMK END
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Dump out a model restart file at a specified model time. Only permanent
!  arrays in the model (which are saved between time steps) need to be
!  dumped for a model restart. For time dependent variables, two time
!  levels (time tpast and tfuture) are needed for a model restart so
!  fields at both time levels are dumped out.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Ming Xue
!  4/01/1992.
!
!  MODIFICATION HISTORY:
!
!  5/06/92 (M. Xue)
!  Added full documentation.
!
!  5/06/92 (M. Xue)
!  Included grid and terrain data in the restart dump.
!
!  6/2/92 (M. Xue and H. Jin)
!  Further facelift.
!
!  2/10/93 (K. Droegemeier)
!  Cleaned up documentation.
!
!  02/07/1995 (Yuhe Liu)
!  Added a new 2-D permanent array, veg(nx,ny), to the argument list
!
!  05/05/1995 (M. Xue)
!  Added rainc and raing into the restart data dump.
!
!  08/22/1995 (M. Xue)
!  Added ptcumsrc and qvcumsrc into the restart data dump.
!
!  08/30/1995 (Yuhe Liu)
!  Added the external boundary data into the restart dump
!
!  2/2/96 (Donghai Wang & Yuhe Liu)
!  Added a 3-D array, mapfct, for map projection factor.
!
!  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.
!
!  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.
!
!  12/09/1998 (Donghai Wang)
!  Added the snow cover.
!
!  13 March 2002 (Eric Kemp)
!  Added arrays for WRF BMJ cumulus scheme.
!
!  April 2002 (Fanyou Kong)
!  Added cnvctopt=5 option for new WRF K-F (KF_ETA) scheme
!
!-----------------------------------------------------------------------
!
!  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
!
!    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
!    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      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)
!    rhostr   Base state density (kg/m**3)times j3.
!    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  Terrain height (m)
!
!    mapfct   Map factors at scalar, u and v points
!
!    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)
!    qvsfc    Effective S.H. at sfc.
!    tsoil     Deep soil temperature (K) (in deep 1 m layer)
!    wetsfc   Surface soil moisture
!    wetdp    Deep soil moisture
!    wetcanp  Canopy water amount
!    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 scale rainfall
!    rainc    Convective rainfall
!
!  OUTPUT:
!
!    None
!
!  WORK ARRAY:
!
!    tem1    Temporary work array.
!
!-----------------------------------------------------------------------
!

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

  INCLUDE 'timelvls.inc'

  INTEGER :: nx,ny,nz          ! 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 :: ptprt (nx,ny,nz,nt)  ! Perturbation potential temperature (K)
  REAL :: pprt  (nx,ny,nz,nt)  ! Perturbation pressure (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 :: 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 :: rhostr(nx,ny,nz)     ! Base state air density (kg/m**3) times j3.
  REAL :: qvbar (nx,ny,nz)     ! Base state water vapor specific humidity
                               ! (kg/kg)

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

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

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

  INTEGER :: nstyps                   ! Number of soil types
  INTEGER :: soiltyp (nx,ny,nstyps)   ! Soil type
  REAL :: stypfrct(nx,ny,nstyps)
  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 1cm layer)
  REAL :: tsoil  (nx,ny,0:nstyps)     ! Deep soil temperature(K)
                                      ! (in deep 1m layer)
  REAL :: wetsfc (nx,ny,0:nstyps)     ! Surface soil moisture
                                      ! in the top 1cm 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 specific humidity
                                      ! at the surface (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 :: 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(IN) :: cldefi(nx,ny)    ! BMJ cloud efficiency
  REAL,INTENT(IN) :: xland(nx,ny)     ! BMJ land mask
                                      ! (1.0 = land, 2.0 = sea)
  REAL,INTENT(IN) :: 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 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.
!
!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------
!
  INTEGER :: ijk, nxyz
  INTEGER :: basrstout ! Control parameter for the base state
                       ! array output
  INTEGER :: grdrstout ! Control parameter for the grid array output
  INTEGER :: icerstout ! Control parameter for the ice variable output
  INTEGER :: sfcrstout ! Control parameter for the surface variable
                       ! output
  INTEGER :: prcrstout ! Control parameter for the precip. rate and rain output
  INTEGER :: rcumout  ! Control parameter for ptcumsrc and qcumsrc output
  INTEGER :: exbcout  ! Control parameter for external boundary output
  INTEGER :: mapfout  ! Control parameter for map factor output
  INTEGER :: radrstout ! Control parameter for radiation forcing output
  INTEGER :: kfrsout  ! Control parameter for Kain-Fritsch output
  INTEGER :: bmjsout  ! Control parameter for WRF BMJ output

  INTEGER :: idummy
  INTEGER :: istat
  INTEGER :: lrstof
  REAL :: rdummy
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'globcst.inc'
  INCLUDE 'grid.inc'          ! Grid & map parameters.
  INCLUDE 'bndry.inc'
  INCLUDE 'exbc.inc'
!
!-----------------------------------------------------------------------
!
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!

!
!-----------------------------------------------------------------------
!
!  Get a name for the restart data file.
!
!-----------------------------------------------------------------------
!

  CALL gtrstfn(runname(1:lfnkey),dirname,ldirnam,curtim,                &
               mgrid,nestgrd, rstoutf, lrstof )

  CALL getunit( rstount )

  OPEN(UNIT=rstount,FILE=trim(rstoutf(1:lrstof)),FORM='unformatted',    &
       STATUS='new',IOSTAT=istat)

  IF( istat /= 0) THEN

    WRITE(6,'(/a,i2,/a/)')                                              &
        ' Error occured when opening restart output file '              &
        //rstoutf(1:lrstof)//                                           &
        ' using FORTRAN unit ',rstount,' Program stopped in RSTOUT.'
    CALL arpsstop('arpsstop called from RSTOUT problem opening file',1)

  END IF

  WRITE(6,'('' DUMPING OUT RESTART FILE AT TIME '',F10.2,               &
  &       ''(s) in FILE '',a,'' using fortran channel no '', i2)')      &
          curtim, rstoutf(1:lrstof),rstount
!
!-----------------------------------------------------------------------
!
!  Write out the restart data:
!
!-----------------------------------------------------------------------
!
  WRITE(rstount) curtim
  WRITE(rstount) nx,ny,nz

  basrstout = 1
  grdrstout = 1
  icerstout = ice
  mapfout  = 1

  prcrstout = 0
  IF ( moist /= 0 ) prcrstout = 1

  sfcrstout = 0
  IF( sfcphy /= 0 ) sfcrstout = 1

  rcumout=0
  IF ( cnvctopt /= 0 ) rcumout=1

  exbcout = 0
  IF ( lbcopt == 2 ) exbcout = 1

  radrstout = 0
  IF ( radopt > 0 ) radrstout = 1

  kfrsout=0
  IF ( cnvctopt == 3 .OR. cnvctopt == 5) kfrsout=1

!EMK BMJ
  bmjsout=0
  IF ( cnvctopt == 4 ) bmjsout=1
!EMK BMJ

  idummy = 0
!EMK BMJ
!  WRITE(rstount) basrstout,grdrstout,icerstout,sfcrstout,prcrstout,     &
!                 rcumout,exbcout,mapfout,radrstout,nstyp,               &
!                 kfrsout,bmjsout,rayklow,idummy,idummy,idummy,          &
!                 idummy,idummy,idummy,idummy,idummy
  WRITE(rstount) basrstout,grdrstout,icerstout,sfcrstout,prcrstout,     &
                 rcumout,exbcout,mapfout,radrstout,nstyp,               &
                 kfrsout,rayklow,bmjsout,idummy,idummy,idummy,          &
                 idummy,idummy,idummy,idummy,idummy
!EMK END

  rdummy = 0.0
  WRITE(rstount) dx,dy,dz,umove,vmove,                                  &
                 xgrdorg,ygrdorg,trulat1,trulat2,trulon,                &
                 sclfct,latitud,ctrlat,ctrlon,rdummy,                   &
                 rdummy,rdummy,rdummy,rdummy,rdummy

  IF( grdrstout == 1) THEN

    WRITE(rstount) x
    WRITE(rstount) y
    WRITE(rstount) z
    WRITE(rstount) zp

  END IF

  IF( basrstout == 1) THEN

    WRITE(rstount) ubar
    WRITE(rstount) vbar
    WRITE(rstount) ptbar
    WRITE(rstount) pbar
    WRITE(rstount) rhostr
    WRITE(rstount) qvbar

  END IF

  CALL cpyary3d(nx,ny,nz,u    (1,1,1,tpast), tem1)
  WRITE(rstount) tem1

  CALL cpyary3d(nx,ny,nz,v    (1,1,1,tpast), tem1)
  WRITE(rstount) tem1

  CALL cpyary3d(nx,ny,nz,w    (1,1,1,tpast), tem1)
  WRITE(rstount) tem1

  CALL cpyary3d(nx,ny,nz,ptprt(1,1,1,tpast), tem1)
  WRITE(rstount) tem1

  CALL cpyary3d(nx,ny,nz,pprt (1,1,1,tpast), tem1)
  WRITE(rstount) tem1

  CALL cpyary3d(nx,ny,nz,qv   (1,1,1,tpast), tem1)
  WRITE(rstount) tem1

  CALL cpyary3d(nx,ny,nz,qc   (1,1,1,tpast), tem1)
  WRITE(rstount) tem1

  CALL cpyary3d(nx,ny,nz,qr   (1,1,1,tpast), tem1)
  WRITE(rstount) tem1


  IF( icerstout /= 0) THEN

    CALL cpyary3d(nx,ny,nz,qi   (1,1,1,tpast), tem1)
    WRITE(rstount) tem1

    CALL cpyary3d(nx,ny,nz,qs   (1,1,1,tpast), tem1)
    WRITE(rstount) tem1

    CALL cpyary3d(nx,ny,nz,qh   (1,1,1,tpast), tem1)
    WRITE(rstount) tem1

  END IF

  CALL cpyary3d(nx,ny,nz,tke  (1,1,1,tpast), tem1)
  WRITE(rstount) tem1

  CALL cpyary3d(nx,ny,nz,u    (1,1,1,tpresent), tem1)
  WRITE(rstount) tem1

  CALL cpyary3d(nx,ny,nz,v    (1,1,1,tpresent), tem1)
  WRITE(rstount) tem1

  CALL cpyary3d(nx,ny,nz,w    (1,1,1,tpresent), tem1)
  WRITE(rstount) tem1

  CALL cpyary3d(nx,ny,nz,ptprt(1,1,1,tpresent), tem1)
  WRITE(rstount) tem1

  CALL cpyary3d(nx,ny,nz,pprt (1,1,1,tpresent), tem1)
  WRITE(rstount) tem1

  CALL cpyary3d(nx,ny,nz,qv   (1,1,1,tpresent), tem1)
  WRITE(rstount) tem1

  CALL cpyary3d(nx,ny,nz,qc   (1,1,1,tpresent), tem1)
  WRITE(rstount) tem1

  CALL cpyary3d(nx,ny,nz,qr   (1,1,1,tpresent), tem1)
  WRITE(rstount) tem1


  IF( icerstout /= 0) THEN

    CALL cpyary3d(nx,ny,nz,qi   (1,1,1,tpresent), tem1)
    WRITE(rstount) tem1

    CALL cpyary3d(nx,ny,nz,qs   (1,1,1,tpresent), tem1)
    WRITE(rstount) tem1

    CALL cpyary3d(nx,ny,nz,qh   (1,1,1,tpresent), tem1)
    WRITE(rstount) tem1

  END IF

  CALL cpyary3d(nx,ny,nz,tke  (1,1,1,tpresent), tem1)
  WRITE(rstount) tem1

  WRITE(rstount) udteb
  WRITE(rstount) udtwb

  WRITE(rstount) vdtnb
  WRITE(rstount) vdtsb

  WRITE(rstount) pdteb
  WRITE(rstount) pdtwb
  WRITE(rstount) pdtnb
  WRITE(rstount) pdtsb

  IF ( sfcrstout /= 0 ) THEN

    PRINT *,'write out sfc/soil variables:'

    WRITE(rstount) soiltyp
    WRITE(rstount) stypfrct
    WRITE(rstount) vegtyp
    WRITE(rstount) lai
    WRITE(rstount) roufns
    WRITE(rstount) veg

    WRITE(rstount) tsfc
    WRITE(rstount) qvsfc
    WRITE(rstount) tsoil
    WRITE(rstount) wetsfc
    WRITE(rstount) wetdp
    WRITE(rstount) wetcanp
    WRITE(rstount) snowdpth

  END IF


  IF ( prcrstout /= 0 ) THEN
    WRITE(rstount) raing
    WRITE(rstount) rainc
    WRITE(rstount) prcrate
  END IF

  IF ( rcumout /=  0 ) THEN

    WRITE(rstount) ptcumsrc
    WRITE(rstount) qcumsrc

  END IF

  IF ( exbcout /= 0 ) THEN
    WRITE(rstount) abstfcst0, abstfcst,                                 &
                   ubcrd,vbcrd,wbcrd,ptbcrd,prbcrd,                     &
                   qvbcrd,qcbcrd,qrbcrd,qibcrd,qsbcrd,qhbcrd

    WRITE(rstount) exbcbuf
  END IF

  IF ( mapfout == 1 ) THEN
    WRITE(rstount) mapfct
  END IF

  IF ( radrstout == 1 ) THEN
    WRITE(rstount) radfrc
    WRITE(rstount) radsw
    WRITE(rstount) rnflx
  END IF

  IF ( kfrsout /= 0 ) THEN

    WRITE(rstount) w0avg
    WRITE(rstount) nca
    WRITE(rstount) kfraincv

  END IF

  IF ( bmjsout /= 0 ) THEN

    WRITE(rstount) cldefi
    WRITE(rstount) xland
    WRITE(rstount) bmjraincv

  END IF

  CLOSE (UNIT=rstount)
  CALL retunit( rstount )
!
!-----------------------------------------------------------------------
!
!  Compress the restart file using system command.
!
!-----------------------------------------------------------------------
!
  IF( filcmprs == 1 ) CALL cmprs( rstoutf(1:lrstof) )

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

!EMK BMJ

SUBROUTINE rstin(nx,ny,nz,nts,nstyps,exbcbufsz,                         & 1,84
           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
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Read in data from a restart file to initialize u,v,w,prprt,pprt,
!  qv,qc,qr,qi,qs,qh and tke at time tpast and tpresent, the base state
!  variables ubar,vbar,ptbar,pbar,rhostr,qvbar, and the time tendencies
!    of variables at the lateral boundaries.
!
!  Fields at tfuture are set to the values at tpresent.
!
!  This subroutine also sets the value of tstart.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Ming Xue
!  4/01/1992.
!
!  MODIFICATION HISTORY:
!
!  5/06/92 (M. Xue)
!  Added full documentation.
!
!  10/15/1992 (M. Xue)
!  Reading of grid and base state arrays added.
!
!  2/10/93 (K. Droegemeier)
!  Cleaned up documentation.
!
!  9/7/93 (Ming Xue)
!  Changed cpyary to cpyary3d.
!
!  9/7/93 (A. Shapiro & Ming Xue)
!  Adjustment to tpast values after umove and vmove are changed.
!
!  02/07/1995 (Yuhe Liu)
!  Added a new 2-D permanent array, veg(nx,ny), to the argument list
!
!  05/05/1995 (M. Xue)
!  Added rainc and raing into the restart data dump.
!
!  08/22/1995 (M. Xue)
!  Added ptcumsrc and qvcumsrc into the restart data dump.
!
!  08/30/1995 (Yuhe Liu)
!  Added the external boundary data into the restart dump
!
!  9/10/1995 (M. Xue)
!  When umove or vmove in arps40.input is 999.0, (umove,vmove)
!  in the restart data is used. No adjustment will be
!  made to the wind fields in this case.
!
!  2/2/96 (Donghai Wang & Yuhe Liu)
!  Added a 3-D array, mapfct, for map projection factor.
!
!  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.
!
!  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.
!
!  12/09/1998 (Donghai Wang)
!  Added the snow cover.
!
!  13 March 2002 (Eric Kemp)
!  Added arrays for WRF BMJ cumulus scheme.
!
!-----------------------------------------------------------------------
!
!  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      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)
!    rhostr   Base state density (kg/m**3) times j3.
!    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  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
!
!    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)
!    qvsfc    Effective S.H. at sfc.
!    tsoil     Deep soil temperature (K) (in deep 1 m layer)
!    wetsfc   Surface soil moisture
!    wetdp    Deep soil moisture
!    wetcanp  Canopy water amount
!    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 scale rainfall
!    rainc    Convective rainfall
!
!    tstart   The time when the time integration starts, which is set to
!             the time of the restart data
!
!    tem1     Temporary work array
!    tem2     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          ! Number of grid points in 3 directions

  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) ! 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 :: rhostr(nx,ny,nz)     ! Base state air density (kg/m**3) time j3.
  REAL :: qvbar (nx,ny,nz)     ! Base state water vapor specific humidity
                               ! (kg/kg)

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

  REAL :: hterain(nx,ny)       ! 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.

  INTEGER :: nstyps                   ! Number of soil types
  INTEGER :: soiltyp (nx,ny,nstyps)   ! Soil type
  REAL :: stypfrct(nx,ny,nstyps)
  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 1cm layer)
  REAL :: qvsfc  (nx,ny,0:nstyps) ! Effective S. H. at the surface (kg/kg)
  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(OUT) :: cldefi(nx,ny)    ! BMJ cloud efficiency
  REAL,INTENT(OUT) :: xland(nx,ny)     ! BMJ land mask
  REAL,INTENT(OUT) :: bmjraincv(nx,ny) ! BMJ convective rainfall (cm)
                                       ! (1.0 = land, 2.0 = sea)
!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 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.

  INTEGER :: grdrstin ! Parameter indicating if the restart data contains
                      ! the grid variables.
  INTEGER :: basrstin ! Parameter indicating if the restart data contains
                      ! the base state variables.
  INTEGER :: icerstin ! Parameter indicating if the restart data contains
                      ! the ice variables.
  INTEGER :: sfcrstin ! Parameter indicating if the restart data contains
                      ! the surface variables.
  INTEGER :: prcrsin  ! Parameter indicating if the restart data contains
                      ! precipitation rate and rainfall
  INTEGER :: rcumin   ! Parameter indicating if the cumulus source terms
                      ! data are present.
  INTEGER :: exbcin   ! Parameter indicating if the external boundary
                      ! data are present.
  INTEGER :: mapfin   ! Parameter indicating if the map factor
                      ! data are present.
  INTEGER :: radrstin ! Parameter indicating if the radiation forcing
                      ! arrays are present.
  INTEGER :: kfrsin   ! Parameter indicating if k-f variable exists
  INTEGER :: bmjsin   ! Parameter indicating if BMJ variable exists

  REAL :: umoveold    ! The domain translation speed of in restart data
  REAL :: vmoveold    ! The domain translation speed of in restart data
  REAL :: uchange     ! Change in domain translation speed
  REAL :: vchange     ! Change in domain translation speed
!
!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------
!
  INTEGER :: tim
  INTEGER :: i, j, k, n, ijk
  INTEGER :: nxin, nyin, nzin, nxyz
  INTEGER :: istat, idummy
  REAL :: datatim,rdummy,dxin,dyin,dzin
  REAL :: amin, amax
  LOGICAL :: fexist,cmprsed
  INTEGER :: lrstfn

  CHARACTER (LEN=80) :: savename
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'globcst.inc'
  INCLUDE 'grid.inc'          ! Grid & map parameters.
  INCLUDE 'bndry.inc'
  INCLUDE 'exbc.inc'
  INCLUDE 'mp.inc'            ! Message passing parameters.

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

  IF(nts == 3 ) THEN
    tpast    = 1
    tpresent = 2
    tfuture  = 3
  ELSE IF(nts == 2 ) THEN
    tpast    = 1
    tpresent = 2
    tfuture  = 2
  ELSE
    tpast    = 1
    tpresent = 1
    tfuture  = 1
  END IF

!   Added a wrapper for scheduling the number of open files...
!   The wrapper will conclude prior to the maxmin checking....
!   note call to jacob is moved outside the file open control loop.
!   due to message passing code in the call to jacob.

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

      CALL getunit( rstiunt )

      lrstfn = 80
      CALL strlnth( rstinf, lrstfn)

      IF (mp_opt > 0) THEN
        savename(1:80) = rstinf(1:80)
        WRITE(rstinf, '(a,a,2i2.2)') trim(savename),'_',loc_x,loc_y
        lrstfn = lrstfn + 5
      END IF

      cmprsed = .false.

      INQUIRE(FILE=rstinf(1:lrstfn), EXIST = fexist )
      IF( fexist ) GO TO 100

      INQUIRE(FILE=rstinf(1:lrstfn)//'.Z', EXIST = fexist )
      IF( fexist ) THEN
        cmprsed = .true.
        CALL uncmprs( rstinf(1:lrstfn)//'.Z' )
        GO TO 100
      END IF

      INQUIRE(FILE=rstinf(1:lrstfn)//'.gz', EXIST = fexist )
      IF( fexist ) THEN
        cmprsed = .true.
        CALL uncmprs( rstinf(1:lrstfn)//'.gz' )
        GO TO 100
      END IF

      CALL wrtcomment('File '//rstinf(1:lrstfn)//                       &
          ' or its compressed version not found.',1)
      CALL arpsstop('arpsstop called from RSTIN compressed file not '// &
           'found',1)

      100   CONTINUE

      OPEN(UNIT=rstiunt,FILE=trim(rstinf(1:lrstfn)),                    &
          FORM='unformatted',STATUS='old',IOSTAT=istat)

      IF (mp_opt > 0) THEN
        rstinf(1:80) = savename(1:80)
        lrstfn = lrstfn - 5
      END IF

      IF( istat /= 0) THEN

        WRITE(6,'(/1x,a,i2/)')                                          &
            'Error occured when opening restart input file '//          &
            rstinf(1:lrstfn)//                                          &
            ' using FORTRAN unit ',rstiunt
      CALL arpsstop('arpsstop called from RSTIN restart file not'//     &
           'found',1)

      END IF

      WRITE(6,'(/1x,a,/1x,a,i2/)')                                      &
          'This is a restart run. Input was read from restart file ',   &
          rstinf(1:lrstfn)//' using fortran unit ',rstiunt
!
!
!-----------------------------------------------------------------------
!
!  Read in the restart data:
!
!-----------------------------------------------------------------------
!
      READ(rstiunt,ERR=999) datatim

      tstart = datatim
      WRITE(6,'(a,f8.1)') ' Restart data is at time ',  datatim

      READ(rstiunt,ERR=999) nxin,nyin,nzin

      IF((nx /= nxin).OR.(ny /= nyin).OR.(nz /= nzin)) THEN

        WRITE(6,'(a,/a,i5,a,i5,a,i5,/a,i5,a,i5,a,i5)')                    &
            ' Array dimension(s) in the restart data inconsistent with ', &
            ' model definitions, dimensions in input data were nx=',nxin, &
            ', ny=',nyin,', nz=',nzin,' the model definitions were nx=',  &
            nx,' ny= ', ny, ' nz= ',nz
        WRITE(6,'(a)') ' Job stopped in subroutine rstin.'
        CALL arpsstop('arpsstop called from RSTIN dimensions '//          &
             'inconsistent',1)

      END IF

!EMK BMJ
!      READ(rstiunt) basrstin,grdrstin,icerstin,sfcrstin,prcrsin,        &
!                rcumin,exbcin,mapfin,radrstin,nstyp,                    &
!                kfrsin,bmjsin,rayklow,idummy,idummy,idummy,             &
!                idummy,idummy,idummy,idummy,idummy
      READ(rstiunt) basrstin,grdrstin,icerstin,sfcrstin,prcrsin,        &
                rcumin,exbcin,mapfin,radrstin,nstyp,                    &
                kfrsin,rayklow,bmjsin,idummy,idummy,                    &
                idummy,idummy,idummy,idummy,idummy
!EMK END

      READ(rstiunt) dx,dy,dz,umoveold,vmoveold,                         &
                xgrdorg,ygrdorg,trulat1,trulat2,trulon,                 &
                sclfct,latitud,ctrlat,ctrlon,rdummy,                    &
                rdummy,rdummy,rdummy,rdummy,rdummy


      IF( grdrstin == 1) THEN

        READ(rstiunt) x
        READ(rstiunt) y
        READ(rstiunt) z
        READ(rstiunt) zp

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

      END IF

      IF( basrstin == 1) THEN

        READ(rstiunt) ubar
        READ(rstiunt) vbar
        READ(rstiunt) ptbar
        READ(rstiunt) pbar
        READ(rstiunt) rhostr
        READ(rstiunt) qvbar

        WRITE(6,'(/1x,a/,1x,a/)')                                       &
            'Base state arrays are read in from restart data set',      &
            'the base state set in INIBASE is superceded.'

      END IF


      tim = tpast

      READ(rstiunt,ERR=999) tem1
      CALL cpyary3d(nx,ny,nz,tem1,u    (1,1,1,tim))

      READ(rstiunt,ERR=999) tem1
      CALL cpyary3d(nx,ny,nz,tem1,v    (1,1,1,tim))

      READ(rstiunt,ERR=999) tem1
      CALL cpyary3d(nx,ny,nz,tem1,w    (1,1,1,tim))

      READ(rstiunt,ERR=999) tem1
      CALL cpyary3d(nx,ny,nz,tem1,ptprt(1,1,1,tim))

      READ(rstiunt,ERR=999) tem1
      CALL cpyary3d(nx,ny,nz,tem1,pprt (1,1,1,tim))

      READ(rstiunt,ERR=999) tem1
      CALL cpyary3d(nx,ny,nz,tem1,qv   (1,1,1,tim))

      READ(rstiunt,ERR=999) tem1
      CALL cpyary3d(nx,ny,nz,tem1,qc   (1,1,1,tim))

      READ(rstiunt,ERR=999) tem1
      CALL cpyary3d(nx,ny,nz,tem1,qr   (1,1,1,tim))

      IF( icerstin /= 0) THEN

        READ(rstiunt,ERR=999) tem1
        CALL cpyary3d(nx,ny,nz,tem1,qi   (1,1,1,tim))

        READ(rstiunt,ERR=999) tem1
        CALL cpyary3d(nx,ny,nz,tem1,qs   (1,1,1,tim))

        READ(rstiunt,ERR=999) tem1
        CALL cpyary3d(nx,ny,nz,tem1,qh   (1,1,1,tim))

      ELSE

        DO k=1,nz
          DO j=1,ny
            DO i=1,nx
              qi(i,j,k,tim) = 0.0
              qs(i,j,k,tim) = 0.0
              qh(i,j,k,tim) = 0.0
            END DO
          END DO
        END DO

      END IF

      READ(rstiunt,ERR=999) tem1
      CALL cpyary3d(nx,ny,nz,tem1,tke   (1,1,1,tim))

      tim = tpresent

      READ(rstiunt,ERR=999) tem1
      CALL cpyary3d(nx,ny,nz,tem1,u    (1,1,1,tim))

      READ(rstiunt,ERR=999) tem1
      CALL cpyary3d(nx,ny,nz,tem1,v    (1,1,1,tim))

      READ(rstiunt,ERR=999) tem1
      CALL cpyary3d(nx,ny,nz,tem1,w    (1,1,1,tim))

      READ(rstiunt,ERR=999) tem1
      CALL cpyary3d(nx,ny,nz,tem1,ptprt(1,1,1,tim))

      READ(rstiunt,ERR=999) tem1
      CALL cpyary3d(nx,ny,nz,tem1,pprt (1,1,1,tim))

      READ(rstiunt,ERR=999) tem1
      CALL cpyary3d(nx,ny,nz,tem1,qv   (1,1,1,tim))

      READ(rstiunt,ERR=999) tem1
      CALL cpyary3d(nx,ny,nz,tem1,qc   (1,1,1,tim))

      READ(rstiunt,ERR=999) tem1
      CALL cpyary3d(nx,ny,nz,tem1,qr   (1,1,1,tim))

      IF( icerstin /= 0) THEN

        READ(rstiunt,ERR=999) tem1
        CALL cpyary3d(nx,ny,nz,tem1,qi   (1,1,1,tim))

        READ(rstiunt,ERR=999) tem1
        CALL cpyary3d(nx,ny,nz,tem1,qs   (1,1,1,tim))

        READ(rstiunt,ERR=999) tem1
        CALL cpyary3d(nx,ny,nz,tem1,qh   (1,1,1,tim))

      ELSE

        DO k=1,nz
          DO j=1,ny
            DO i=1,nx
              qi(i,j,k,tim) = 0.0
              qs(i,j,k,tim) = 0.0
              qh(i,j,k,tim) = 0.0
            END DO
          END DO
        END DO

      END IF

      READ(rstiunt,ERR=999) tem1
      CALL cpyary3d(nx,ny,nz,tem1,tke  (1,1,1,tim))

      READ(rstiunt,ERR=999) udteb
      READ(rstiunt,ERR=999) udtwb

      READ(rstiunt,ERR=999) vdtnb
      READ(rstiunt,ERR=999) vdtsb

      READ(rstiunt,ERR=999) pdteb
      READ(rstiunt,ERR=999) pdtwb
      READ(rstiunt,ERR=999) pdtnb
      READ(rstiunt,ERR=999) pdtsb

!-----------------------------------------------------------------------
!
!  Set the future values of variables to their current values.
!  This is done primarily for safety reasons since the arrays at
!  tfuture will be overwritten by the new values during the
!  time integration.
!
!-----------------------------------------------------------------------
!
      CALL cpyary3d(nx,ny,nz,u (1,1,1,tpresent) , u (1,1,1,tfuture))
      CALL cpyary3d(nx,ny,nz,v (1,1,1,tpresent) , v (1,1,1,tfuture))
      CALL cpyary3d(nx,ny,nz,w (1,1,1,tpresent) , w (1,1,1,tfuture))
      CALL cpyary3d(nx,ny,nz,ptprt(1,1,1,tpresent),                     &
          ptprt(1,1,1,tfuture))
      CALL cpyary3d(nx,ny,nz,pprt (1,1,1,tpresent),                     &
          pprt (1,1,1,tfuture))
      CALL cpyary3d(nx,ny,nz,qv(1,1,1,tpresent) , qv(1,1,1,tfuture))
      CALL cpyary3d(nx,ny,nz,qc(1,1,1,tpresent) , qc(1,1,1,tfuture))
      CALL cpyary3d(nx,ny,nz,qr(1,1,1,tpresent) , qr(1,1,1,tfuture))
      CALL cpyary3d(nx,ny,nz,qi(1,1,1,tpresent) , qi(1,1,1,tfuture))
      CALL cpyary3d(nx,ny,nz,qs(1,1,1,tpresent) , qs(1,1,1,tfuture))
      CALL cpyary3d(nx,ny,nz,qh(1,1,1,tpresent) , qh(1,1,1,tfuture))
      CALL cpyary3d(nx,ny,nz,tke(1,1,1,tpresent) ,                      &
          tke(1,1,1,tfuture))

      IF ( sfcrstin /= 0 ) THEN

        PRINT *,'read in sfc/soil variables:'

        READ(rstiunt,ERR=999) soiltyp
        READ(rstiunt,ERR=999) stypfrct

        READ(rstiunt,ERR=999) vegtyp
        READ(rstiunt,ERR=999) lai
        READ(rstiunt,ERR=999) roufns
        READ(rstiunt,ERR=999) veg

        READ(rstiunt,ERR=999) tsfc
        READ(rstiunt,ERR=999) qvsfc
        READ(rstiunt,ERR=999) tsoil
        READ(rstiunt,ERR=999) wetsfc
        READ(rstiunt,ERR=999) wetdp
        READ(rstiunt,ERR=999) wetcanp
        READ(rstiunt,ERR=999) snowdpth

      END IF

      IF ( prcrsin /= 0 ) THEN
        READ(rstiunt,ERR=999) raing
        READ(rstiunt,ERR=999) rainc
        READ(rstiunt,ERR=999) prcrate
      END IF

      IF ( rcumin /= 0 ) THEN
        READ(rstiunt,ERR=999) ptcumsrc
        READ(rstiunt,ERR=999) qcumsrc
      END IF

      IF ( exbcin /= 0 ) THEN
        IF ( lbcopt == 2 ) THEN
          READ(rstiunt,ERR=999) abstfcst0, abstfcst,                    &
                   ubcrd,vbcrd,wbcrd,ptbcrd,prbcrd,                     &
                   qvbcrd,qcbcrd,qrbcrd,qibcrd,qsbcrd,qhbcrd

          READ(rstiunt,ERR=999) exbcbuf
        ELSE
          WRITE(6,'(a/a/a/a)')                                          &
              'WARNING: The restart file contains EXBC arrays, while',  &
              '         the this run does not have EXBC option.',       &
              '         Therefore, the results from restart run may be', &
              '         alterred. The program will continue.'

          READ(rstiunt,ERR=999)

          READ(rstiunt,ERR=999)
          READ(rstiunt,ERR=999)
          READ(rstiunt,ERR=999)
          READ(rstiunt,ERR=999)
          READ(rstiunt,ERR=999)
          READ(rstiunt,ERR=999)
          READ(rstiunt,ERR=999)
          READ(rstiunt,ERR=999)
          READ(rstiunt,ERR=999)
          READ(rstiunt,ERR=999)
          READ(rstiunt,ERR=999)
        END IF
      END IF

      IF ( mapfin == 1 ) THEN
        READ(rstiunt,ERR=999) mapfct
      END IF

      IF ( radrstin == 1 ) THEN
        READ(rstiunt,ERR=999) radfrc
        READ(rstiunt,ERR=999) radsw
        READ(rstiunt,ERR=999) rnflx
      END IF

      IF ( kfrsin /= 0) THEN
        READ(rstiunt,ERR=999) w0avg
        READ(rstiunt,ERR=999) nca
        READ(rstiunt,ERR=999) kfraincv
      END IF

      IF ( bmjsin /= 0) THEN
        READ(rstiunt,ERR=999) cldefi
        READ(rstiunt,ERR=999) xland
        READ(rstiunt,ERR=999) bmjraincv
      END IF

      CLOSE (UNIT=rstiunt)
      CALL retunit( rstiunt )
!
!-----------------------------------------------------------------------
!
!  Reset the model u and v velocity values using the new
!  domain translation speed.
!
!-----------------------------------------------------------------------
!

      IF( nint(umove) == 999 .OR. nint(vmove) == 999 ) THEN

        umove = umoveold
        vmove = vmoveold

      ELSE IF (umoveold /= umove .OR. vmoveold /= vmove ) THEN

        WRITE(6,'(3(/1x,a)/)')                                          &
            'ATTENTION: UMOVE or VMOVE in the input file were different ', &
            'from those in the restart file. Subroutine ADJUVMV is called', &
            'to adjust the time-dependent variables for option grdtrns!=0.'

        IF ( grdtrns /= 0 ) THEN
          uchange = umove - umoveold
          vchange = vmove - vmoveold

          CALL adjuvmv(nx,ny,nz,                                        &
              ubar,vbar,u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh,qvbar,       &
              uchange, vchange, tem1, tem2)

        END IF
      END IF

    END IF  ! end of FOPEN wrapper for file read/write...
    IF (mp_opt > 0) CALL mpbarrier
  END DO

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

  WRITE(6,'(/1x,a/,1x,a/)')                                             &
      'Grid definition arrays are read in from initialization data',    &
      'those set in INIGRD are superceded.'

!
!-----------------------------------------------------------------------
!
!  Print out the domain-wide max/min of output variables.
!
!-----------------------------------------------------------------------

  WRITE(6,'(/1x,a/)')                                                   &
      'Min. and max. of the data arrays read in from restart data:'

  CALL a3dmax0(x,1,nx,1,nx,1,1,1,1, 1,1,1,1, amax,amin)
  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)
  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)
  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)
  WRITE(6,'(1x,2(a,e13.6))') 'zpmin   = ', amin,',  zpmax   =',amax

  CALL a3dmax0(hterain,1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1, amax,amin)
  WRITE(6,'(1x,2(a,e13.6))') 'hmin    = ', amin,',  hmax    =',amax

  CALL a3dmax0(ubar,1,nx,1,nx,1,ny,1,ny-1,1,nz,1,nz-1, amax,amin)
  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)
  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)
  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)
  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)
  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)
  WRITE(6,'(1x,2(a,e13.6))') 'qvbarmin= ', amin,',  qvbarmax=',amax


  DO i = 1,2

    IF( i == 1) THEN
      WRITE(6,'(/1x,a/)') 'Min/max of fields at tpresent:'
      tim = tpresent
    ELSE
      WRITE(6,'(/1x,a/)') 'Min/max of fields at tpast:'
      tim = tpast
    END IF

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

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

    CALL a3dmax0(                 amax,amin)
    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)
    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)
    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)
    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)
    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)
    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)
    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)
    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)
    WRITE(6,'(1x,2(a,e13.6))') 'qhmin   = ', amin,',  qhmax   =',       &
          amax

  END DO

  WRITE(6,'(/1x,a/)')                                                   &
        'Min/max of fields for other one time level arrays:'

  CALL a3dmax0(tsfc(1,1,0),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,1,0),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,1,0),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,1,0),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,1,0),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

  CALL a3dmax0(raing,1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1, amax,amin)
  WRITE(6,'(1x,2(a,e13.6))') 'raingmin= ', amin,', raingmax =',amax

  CALL a3dmax0(rainc,1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1, amax,amin)
  WRITE(6,'(1x,2(a,e13.6))') 'raincmin= ', amin,', raincmax =',amax

  CALL a3dmax0(ptcumsrc,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,            &
               amax,amin)
  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)
  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)
  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)
  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)
  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)
  WRITE(6,'(1x,2(a,e13.6))') 'qscummin= ', amin,',  qscummax=',amax
!
!-----------------------------------------------------------------------
!
!  Compress the restart file if it was originally compressed.
!
!-----------------------------------------------------------------------
!
  IF( cmprsed .AND. filcmprs == 1 ) THEN
    CALL cmprs( rstinf(1:lrstfn) )
  END IF

  RETURN

  999   CONTINUE
  WRITE(6,'(a)') ' Error reading restart data '//rstinf
  WRITE(6,'(a,i3,a)') ' Fortran channel ',rstiunt,' was used.'
  WRITE(6,'(a)') ' Job stopped in subroutine rstin!'
  CLOSE (UNIT=rstiunt)

  CALL arpsstop('arpsstop called from RSTIN error reading restart'//   &
                 'file ',1)

END SUBROUTINE rstin

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


SUBROUTINE cpyary3d(nx,ny,nz, ain, aout) 60
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Copy the contents of array 'ain' into 'aout'.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Ming Xue
!  11/10/1992.
!
!  MODIFICATION HISTORY:
!
!  2/11/93 (K. Droegemeier)
!  Added full documentation.
!
!-----------------------------------------------------------------------
!
!  INPUT :
!
!    ain      Input array
!    nx       1st Dimension of input and output arrays.
!    ny       2nd Dimension of input and output arrays.
!    nz       3rd Dimension of input and output arrays.
!
!  OUTPUT:
!
!    aout     Output array.
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
  INTEGER :: nx,ny,nz
  REAL :: ain (nx,ny,nz) ! Input array to be copied in aout.
  REAL :: aout(nx,ny,nz) ! Array whose value will be copied fron ain.
  INTEGER :: i,j,k

!-----------------------------------------------------------------------
!
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  DO k=1,nz
    DO j=1,ny
      DO i=1,nx
        aout(i,j,k) = ain(i,j,k)
      END DO
    END DO
  END DO

  RETURN
END SUBROUTINE cpyary3d