c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE RSTOUT                     ######
c     ######                                                      ######
c     ######                Copyright (c) 1992                    ######
c     ######  Center for the Analysis and Prediction of Storms    ######
c     ######    University of Oklahoma.  All rights reserved.     ######
c     ######                                                      ######
c     ##################################################################
c     ##################################################################
c
      SUBROUTINE RSTOUT(nx,ny,nz, 
     :           u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh,
     :           udteb, udtwb, udtnb, udtsb,
     :           vdteb, vdtwb, vdtnb, vdtsb,
     :           wdteb, wdtwb, wdtnb, wdtsb,
     :           ptdteb,ptdtwb,ptdtnb,ptdtsb,
     :           pdteb ,pdtwb ,pdtnb ,pdtsb,
     :           qvdteb,qvdtwb,qvdtnb,qvdtsb,
     :           qcdteb,qcdtwb,qcdtnb,qcdtsb,
     :           qrdteb,qrdtwb,qrdtnb,qrdtsb,
     :           qidteb,qidtwb,qidtnb,qidtsb,
     :           qsdteb,qsdtwb,qsdtnb,qsdtsb,
     :           qhdteb,qhdtwb,qhdtnb,qhdtsb,
     :           ubar,vbar,ptbar,pbar,rhobar,qvbar,
     :           x,y,z,zp,hterain )     
c
c#######################################################################
c
c     PURPOSE:
c
c     Dump out a model restart file at a given model time. Only permanent 
c     arrays in the model (which are saved between time steps) need to be
c     dumped for a model restart. For time dependent variables, two time 
c     levels (time tpast and tfuture) are needed for a model restart so 
c     fields at both time levels are dumped out.
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c     4/01/1992.  
c
c     MODIFICATION HISTORY:
c
c     5/06/92 (M. Xue)
c     Added full documentation.
c
c     5/06/92 (M. Xue)
c     Included grid and terrain data in the restart dump.
c
c     6/2/92 (M. Xue and H. Jin)
c     Further facelift.
C
c#######################################################################
c
c     INPUT:
c
c
c       nx       Number of grid points in the x-direction (east/west)
c       ny       Number of grid points in the y-direction (north/south)
c       nz       Number of grid points in the vertical
c
c       u        x component of velocity at times tpast and tpresent (m/s)
c       v        y component of velocity at times tpast and tpresent (m/s)
c       w        Vertical component of Cartesian velocity at times
c       ptprt    Perturbation potential temperature at times tpast and
c                tpresent (K)
c       pprt     Perturbation pressure at times tpast and tpresent (Pascal)
c
c       qv       Water vapor specific humidity at times tpast and tpresent (kg/kg)
c       qc       Cloud water mixing ratio at times tpast and tpresent (kg/kg)
c       qr       Rainwater mixing ratio at times tpast and tpresent (kg/kg)
c       qi       Cloud ice mixing ratio at times tpast and tpresent (kg/kg)
c       qs       Snow mixing ratio at times tpast and tpresent (kg/kg)
c       qh       Hail mixing ratio at times tpast and tpresent (kg/kg)
c
c       udteb    Time tendency of u field at east boundary (m/s**2)
c       udtwb    Time tendency of u field at west boundary (m/s**2)
c       udtnb    Time tendency of u field at north boundary (m/s**2)
c       udtsb    Time tendency of u field at south boundary (m/s**2)
c
c       vdteb    Time tendency of v field at east boundary (m/s**2)
c       vdtwb    Time tendency of v field at west boundary (m/s**2)
c       vdtnb    Time tendency of v field at north boundary (m/s**2)
c       vdtsb    Time tendency of v field at south boundary (m/s**2)
c
c       wdteb    Time tendency of w field at east boundary (m/s**2)
c       wdtwb    Time tendency of w field at west boundary (m/s**2)
c       wdtnb    Time tendency of w field at north boundary (m/s**2)
c       wdtsb    Time tendency of w field at south boundary (m/s**2)
c
c       ptdteb   Time tendency of ptprt field at east boundary (K/s)
c       ptdtwb   Time tendency of ptprt field at west boundary (K/s)
c       ptdtnb   Time tendency of ptprt field at north boundary (K/s)
c       ptdtsb   Time tendency of ptprt field at south boundary (K/s)
c
c       pdteb    Time tendency of pprt field at east boundary (PASCAL/s)
c       pdtwb    Time tendency of pprt field at west boundary (PASCAL/s)
c       pdtnb    Time tendency of pprt field at north boundary (PASCAL/s)
c       pdtsb    Time tendency of pprt field at south boundary (PASCAL/s)
c
c       qvdteb   Time tendency of qv field at east boundary (kg/(kg*s))
c       qvdtwb   Time tendency of qv field at west boundary (kg/(kg*s))
c       qvdtnb   Time tendency of qv field at north boundary (kg/(kg*s))
c       qvdtsb   Time tendency of qv field at south boundary (kg/(kg*s))
c
c       qcdteb   Time tendency of qc field at east boundary (kg/(kg*s))
c       qcdtwb   Time tendency of qc field at west boundary (kg/(kg*s))
c       qcdtnb   Time tendency of qc field at north boundary (kg/(kg*s))
c       qcdtsb   Time tendency of qc field at south boundary (kg/(kg*s))
c
c       qrdteb   Time tendency of qr field at east boundary (kg/(kg*s))
c       qrdtwb   Time tendency of qr field at west boundary (kg/(kg*s))
c       qrdtnb   Time tendency of qr field at north boundary (kg/(kg*s))
c       qrdtsb   Time tendency of qr field at south boundary (kg/(kg*s))
c
c       qidteb   Time tendency of qi field at east boundary (kg/(kg*s))
c       qidtwb   Time tendency of qi field at west boundary (kg/(kg*s))
c       qidtnb   Time tendency of qi field at north boundary (kg/(kg*s))
c       qidtsb   Time tendency of qi field at south boundary (kg/(kg*s))
c
c       qsdteb   Time tendency of qs field at east boundary (kg/(kg*s))
c       qsdtwb   Time tendency of qs field at west boundary (kg/(kg*s))
c       qsdtnb   Time tendency of qs field at north boundary (kg/(kg*s))
c       qsdtsb   Time tendency of qs field at south boundary (kg/(kg*s))
c
c       qhdteb   Time tendency of qh field at east boundary (kg/(kg*s))
c       qhdtwb   Time tendency of qh field at west boundary (kg/(kg*s))
c       qhdtnb   Time tendency of qh field at north boundary (kg/(kg*s))
c       qhdtsb   Time tendency of qh field at south boundary (kg/(kg*s))
c
c       ubar     Base state zonal velocity component (m/s)
c       vbar     Base state meridional velocity component (m/s)
c       ptbar    Base state potential temperature (K)
c       pbar     Base state pressure (Pascal)
c       rhobar   Base state density (kg/m**3)
c       qvbar    Base state water vapor specific humidity (kg/kg) 
c
c       x        x coordinate of grid points in physical/comp. space (m)
c       y        y coordinate of grid points in physical/comp. space (m)
c       z        z coordinate of grid points in computational space (m)
c       zp       Vertical coordinate of grid points in physical space (m)
c       hterain  The height of terrain (m)
c
c     OUTPUT:
c
c       None
c
c####################################################################### 
c

c
c#######################################################################
c
c     Variable Declarations.
c
c#######################################################################
c
      implicit none 

      integer nt           ! The no. of t-levels of t-dependent arrays.
      integer tpast        ! Index of time level for the past time.
      integer tpresent     ! Index of time level for the present time.
      integer tfuture      ! Index of time level for the future time.

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

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

      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 ptdteb(ny,nz)        ! T-tendency of ptprt at e-boundary (K/s)
      real ptdtwb(ny,nz)        ! T-tendency of ptprt at w-boundary (K/s)
      real ptdtnb(nx,nz)        ! T-tendency of ptprt at n-boundary (K/s)
      real ptdtsb(nx,nz)        ! T-tendency of ptprt at s-boundary (K/s)

      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 qvdteb(ny,nz)        ! T-tendency of qv at e-boundary (kg/(kg*s))
      real qvdtwb(ny,nz)        ! T-tendency of qv at w-boundary (kg/(kg*s))
      real qvdtnb(nx,nz)        ! T-tendency of qv at n-boundary (kg/(kg*s))
      real qvdtsb(nx,nz)        ! T-tendency of qv at s-boundary (kg/(kg*s))

      real qcdteb(ny,nz)        ! T-tendency of qc at e-boundary (kg/(kg*s))
      real qcdtwb(ny,nz)        ! T-tendency of qc at w-boundary (kg/(kg*s))
      real qcdtnb(nx,nz)        ! T-tendency of qc at n-boundary (kg/(kg*s))
      real qcdtsb(nx,nz)        ! T-tendency of qc at s-boundary (kg/(kg*s))

      real qrdteb(ny,nz)        ! T-tendency of qr at e-boundary (kg/(kg*s))
      real qrdtwb(ny,nz)        ! T-tendency of qr at w-boundary (kg/(kg*s))
      real qrdtnb(nx,nz)        ! T-tendency of qr at n-boundary (kg/(kg*s))
      real qrdtsb(nx,nz)        ! T-tendency of qr at s-boundary (kg/(kg*s))

      real qidteb(ny,nz)        ! T-tendency of qi at e-boundary (kg/(kg*s))
      real qidtwb(ny,nz)        ! T-tendency of qi at w-boundary (kg/(kg*s))
      real qidtnb(nx,nz)        ! T-tendency of qi at n-boundary (kg/(kg*s))
      real qidtsb(nx,nz)        ! T-tendency of qi at s-boundary (kg/(kg*s))

      real qsdteb(ny,nz)        ! T-tendency of qs at e-boundary (kg/(kg*s))
      real qsdtwb(ny,nz)        ! T-tendency of qs at w-boundary (kg/(kg*s))
      real qsdtnb(nx,nz)        ! T-tendency of qs at n-boundary (kg/(kg*s))
      real qsdtsb(nx,nz)        ! T-tendency of qs at s-boundary (kg/(kg*s))

      real qhdteb(ny,nz)        ! T-tendency of qh at e-boundary (kg/(kg*s))
      real qhdtwb(ny,nz)        ! T-tendency of qh at w-boundary (kg/(kg*s))
      real qhdtnb(nx,nz)        ! T-tendency of qh at n-boundary (kg/(kg*s))
      real qhdtsb(nx,nz)        ! T-tendency of qh at s-boundary (kg/(kg*s))
c
      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 rhobar(nx,ny,nz)     ! Base state air density (kg/m**3)
      real qvbar (nx,ny,nz)     ! Base state water vapor specific humidity
                                ! (kg/kg)
c
      real x     (nx)           ! The x-coord. of the physical and
                                ! computational grid. Defined at u-point.
      real y     (ny)           ! The y-coord. of the physical and
                                ! computational grid. Defined at v-point.
      real z     (nz)           ! The z-coord. of the computational grid.
                                ! Defined at w-point on the staggered grid.
      real zp    (nx,ny,nz)     ! The physical height coordinate defined at
                                ! w-point of the staggered grid.

      real hterain(nx,ny)       ! The height of the terrain (m).
c
c#######################################################################
c
c     Misc. local variables:
c
c#######################################################################
c
      integer i, j, k
      integer istat
      integer lrstof
      character temchar*80
c
c#######################################################################
c
c     Include files: 
c
c#######################################################################
c
      include 'globcst.inc'
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c

c
c#######################################################################
c
c     Construct the file name:
c
c     The format of the restart file names is
c
c                    XXXXXX.rstnnnnnn
c     
c     where XXXXXX are 6 characters specifying this model run, which is
c     is given by character string runname. nnnnnn is a 6 digit integer
c     number indicating the time of the data in seconds.
c
c#######################################################################
c

      return
      write(rstoutf,'(a,i6.6)') runname(1:6)//'.rst',nint(curtim)
      lrstof = 16


      IF( dirname.ne.' ' ) THEN

        temchar = rstoutf(1:lrstof)
        rstoutf = dirname(1:ldirnam)//'/'//temchar
        lrstof  = lrstof + ldirnam + 1

      ENDIF

      open(unit=rstount,file=rstoutf(1:lrstof),form='unformatted',
     :     status='unknown',iostat=istat)

      IF( istat.ne.0) THEN

c        write(6,'(/a,i2,/a/)')
c     :  ' Error occured when opening restart output file '
c     :  //rstoutf(1:lrstof)//
c     :  ' using FORTRAN unit ',rstount,' Program stopped in RSTOUT.'
        STOP

      ENDIF

      write(6,'('' Dumping out restart file at time '',f8.2,
     :        ''(s) in file '',a,'' using fortran channel no '', i2)') 
     :        curtim, rstoutf(1:16),rstount
c
c#######################################################################
c
c     Write out the restart data:
c
c#######################################################################
c
      write(rstount) curtim
      write(rstount) nx,ny,nz

      write(rstount) (((u    (i,j,k,tpast),i=1,nx),j=1,ny),k=1,nz)
      write(rstount) (((v    (i,j,k,tpast),i=1,nx),j=1,ny),k=1,nz)
      write(rstount) (((w    (i,j,k,tpast),i=1,nx),j=1,ny),k=1,nz)
      write(rstount) (((ptprt(i,j,k,tpast),i=1,nx),j=1,ny),k=1,nz)
      write(rstount) (((pprt (i,j,k,tpast),i=1,nx),j=1,ny),k=1,nz)
      write(rstount) (((qv   (i,j,k,tpast),i=1,nx),j=1,ny),k=1,nz)
      write(rstount) (((qc   (i,j,k,tpast),i=1,nx),j=1,ny),k=1,nz)
      write(rstount) (((qr   (i,j,k,tpast),i=1,nx),j=1,ny),k=1,nz)

      IF( iceout. ne. 0) THEN 

        write(rstount) (((qi   (i,j,k,tpast),i=1,nx),j=1,ny),k=1,nz)
        write(rstount) (((qs   (i,j,k,tpast),i=1,nx),j=1,ny),k=1,nz)
        write(rstount) (((qh   (i,j,k,tpast),i=1,nx),j=1,ny),k=1,nz)

      ENDIF
      
      write(rstount) (((u    (i,j,k,tpresent),i=1,nx),j=1,ny),k=1,nz)
      write(rstount) (((v    (i,j,k,tpresent),i=1,nx),j=1,ny),k=1,nz)
      write(rstount) (((w    (i,j,k,tpresent),i=1,nx),j=1,ny),k=1,nz)
      write(rstount) (((ptprt(i,j,k,tpresent),i=1,nx),j=1,ny),k=1,nz)
      write(rstount) (((pprt (i,j,k,tpresent),i=1,nx),j=1,ny),k=1,nz)
      write(rstount) (((qv   (i,j,k,tpresent),i=1,nx),j=1,ny),k=1,nz)
      write(rstount) (((qc   (i,j,k,tpresent),i=1,nx),j=1,ny),k=1,nz)
      write(rstount) (((qr   (i,j,k,tpresent),i=1,nx),j=1,ny),k=1,nz)

      IF( iceout. ne. 0) THEN 

        write(rstount) (((qi   (i,j,k,tpresent),i=1,nx),j=1,ny),k=1,nz)
        write(rstount) (((qs   (i,j,k,tpresent),i=1,nx),j=1,ny),k=1,nz)
        write(rstount) (((qh   (i,j,k,tpresent),i=1,nx),j=1,ny),k=1,nz)

      ENDIF

      write(rstount) ((udteb (j,k),j=1,ny),k=1,nz)
      write(rstount) ((udtwb (j,k),j=1,ny),k=1,nz)
      write(rstount) ((udtnb (i,k),i=1,nx),k=1,nz)
      write(rstount) ((udtsb (i,k),i=1,nx),k=1,nz)

      write(rstount) ((vdteb (j,k),j=1,ny),k=1,nz)
      write(rstount) ((vdtwb (j,k),j=1,ny),k=1,nz)
      write(rstount) ((vdtnb (i,k),i=1,nx),k=1,nz)
      write(rstount) ((vdtsb (i,k),i=1,nx),k=1,nz)

      write(rstount) ((wdteb (j,k),j=1,ny),k=1,nz)
      write(rstount) ((wdtwb (j,k),j=1,ny),k=1,nz)
      write(rstount) ((wdtnb (i,k),i=1,nx),k=1,nz)
      write(rstount) ((wdtsb (i,k),i=1,nx),k=1,nz)

      write(rstount) ((ptdteb(j,k),j=1,ny),k=1,nz)
      write(rstount) ((ptdtwb(j,k),j=1,ny),k=1,nz)
      write(rstount) ((ptdtnb(i,k),i=1,nx),k=1,nz)
      write(rstount) ((ptdtsb(i,k),i=1,nx),k=1,nz)

      write(rstount) ((pdteb (j,k),j=1,ny),k=1,nz)
      write(rstount) ((pdtwb (j,k),j=1,ny),k=1,nz)
      write(rstount) ((pdtnb (i,k),i=1,nx),k=1,nz)
      write(rstount) ((pdtsb (i,k),i=1,nx),k=1,nz)

      write(rstount) ((qvdteb(j,k),j=1,ny),k=1,nz)
      write(rstount) ((qvdtwb(j,k),j=1,ny),k=1,nz)
      write(rstount) ((qvdtnb(i,k),i=1,nx),k=1,nz)
      write(rstount) ((qvdtsb(i,k),i=1,nx),k=1,nz)

      write(rstount) ((qcdteb(j,k),j=1,ny),k=1,nz)
      write(rstount) ((qcdtwb(j,k),j=1,ny),k=1,nz)
      write(rstount) ((qcdtnb(i,k),i=1,nx),k=1,nz)
      write(rstount) ((qcdtsb(i,k),i=1,nx),k=1,nz)

      write(rstount) ((qrdteb(j,k),j=1,ny),k=1,nz)
      write(rstount) ((qrdtwb(j,k),j=1,ny),k=1,nz)
      write(rstount) ((qrdtnb(i,k),i=1,nx),k=1,nz)
      write(rstount) ((qrdtsb(i,k),i=1,nx),k=1,nz)

      write(rstount) ((qidteb(j,k),j=1,ny),k=1,nz)
      write(rstount) ((qidtwb(j,k),j=1,ny),k=1,nz)
      write(rstount) ((qidtnb(i,k),i=1,nx),k=1,nz)
      write(rstount) ((qidtsb(i,k),i=1,nx),k=1,nz)

      write(rstount) ((qsdteb(j,k),j=1,ny),k=1,nz)
      write(rstount) ((qsdtwb(j,k),j=1,ny),k=1,nz)
      write(rstount) ((qsdtnb(i,k),i=1,nx),k=1,nz)
      write(rstount) ((qsdtsb(i,k),i=1,nx),k=1,nz)

      write(rstount) ((qhdteb(j,k),j=1,ny),k=1,nz)
      write(rstount) ((qhdtwb(j,k),j=1,ny),k=1,nz)
      write(rstount) ((qhdtnb(i,k),i=1,nx),k=1,nz)
      write(rstount) ((qhdtsb(i,k),i=1,nx),k=1,nz)

      write(rstount) (((ubar (i,j,k),i=1,nx),j=1,ny),k=1,nz)
      write(rstount) (((vbar (i,j,k),i=1,nx),j=1,ny),k=1,nz)
      write(rstount) (((ptbar(i,j,k),i=1,nx),j=1,ny),k=1,nz)
      write(rstount) (((pbar (i,j,k),i=1,nx),j=1,ny),k=1,nz)
      write(rstount) (((rhobar(i,j,k),i=1,nx),j=1,ny),k=1,nz)
      write(rstount) (((qvbar(i,j,k),i=1,nx),j=1,ny),k=1,nz)


      close (unit=rstount)

      RETURN
      END
c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE RSTIN                      ######
c     ######                                                      ######
c     ######                Copyright (c) 1992                    ######
c     ######  Center for the Analysis and Prediction of Storms    ######
c     ######    University of Oklahoma.  All rights reserved.     ######
c     ######                                                      ######
c     ##################################################################
c     ##################################################################
c

      SUBROUTINE RSTIN(nx,ny,nz, 
     :           u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh,
     :           udteb, udtwb, udtnb, udtsb,
     :           vdteb, vdtwb, vdtnb, vdtsb,
     :           wdteb, wdtwb, wdtnb, wdtsb,
     :           ptdteb,ptdtwb,ptdtnb,ptdtsb,
     :           pdteb ,pdtwb ,pdtnb ,pdtsb,
     :           qvdteb,qvdtwb,qvdtnb,qvdtsb,
     :           qcdteb,qcdtwb,qcdtnb,qcdtsb,
     :           qrdteb,qrdtwb,qrdtnb,qrdtsb,
     :           qidteb,qidtwb,qidtnb,qidtsb,
     :           qsdteb,qsdtwb,qsdtnb,qsdtsb,
     :           qhdteb,qhdtwb,qhdtnb,qhdtsb,
     :           ubar,vbar,ptbar,pbar,rhobar,qvbar,
     :           x,y,z,zp,hterain )     
c
c#######################################################################
c
c     PURPOSE:
c
c     Read in restart data from a restart file to initialize fields
c     u,v,w,prprt,pprt,qv,qc,qr,qi,qs,qh at time tpast and tpresent,
c     the base state variables ubar,vbar,ptbar,pbar,rhobar,qvbar,
c     and the time tendencies of variables at the lateral boundaries.
c
c     Fields at tfuture are set to the values at tpresent.
c
c     This subroutine also sets the value of tstart.
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c     4/01/1992.  
c
c     MODIFICATION HISTORY:
c
c     5/06/92 (M. Xue)
c     Added full documentation.
c
c     6/2/92 (M. Xue and H. Jin)
c     Further facelift.
C
c#######################################################################
c
c     INPUT : 
c
c       nx       Number of grid points in the x-direction (east/west)
c       ny       Number of grid points in the y-direction (north/south)
c       nz       Number of grid points in the vertical
c
c     OUTPUT:
c
c       u        x component of velocity at times tpast and tpresent (m/s)
c       v        y component of velocity at times tpast and tpresent (m/s)
c       w        Vertical component of Cartesian velocity at times
c                tpast and tpresent (m/s)
c       ptprt    Perturbation potential temperature at times tpast and
c                tpresent (K)
c       pprt     Perturbation pressure at times tpast and tpresent (Pascal)
c
c       qv       Water vapor specific humidity at times tpast and tpresent (kg/kg)
c       qc       Cloud water mixing ratio at times tpast and tpresent (kg/kg)
c       qr       Rainwater mixing ratio at times tpast and tpresent (kg/kg)
c       qi       Cloud ice mixing ratio at times tpast and tpresent (kg/kg)
c       qs       Snow mixing ratio at times tpast and tpresent (kg/kg)
c       qh       Hail mixing ratio at times tpast and tpresent (kg/kg)
c
c       udteb    Time tendency of u field at east boundary (m/s**2)
c       udtwb    Time tendency of u field at west boundary (m/s**2)
c       udtnb    Time tendency of u field at north boundary (m/s**2)
c       udtsb    Time tendency of u field at south boundary (m/s**2)
c
c       vdteb    Time tendency of v field at east boundary (m/s**2)
c       vdtwb    Time tendency of v field at west boundary (m/s**2)
c       vdtnb    Time tendency of v field at north boundary (m/s**2)
c       vdtsb    Time tendency of v field at south boundary (m/s**2)
c
c       wdteb    Time tendency of w field at east boundary (m/s**2)
c       wdtwb    Time tendency of w field at west boundary (m/s**2)
c       wdtnb    Time tendency of w field at north boundary (m/s**2)
c       wdtsb    Time tendency of w field at south boundary (m/s**2)
c
c       ptdteb   Time tendency of ptprt field at east boundary (K/s)
c       ptdtwb   Time tendency of ptprt field at west boundary (K/s)
c       ptdtnb   Time tendency of ptprt field at north boundary (K/s)
c       ptdtsb   Time tendency of ptprt field at south boundary (K/s)
c
c       pdteb    Time tendency of pprt field at east boundary (PASCAL/s)
c       pdtwb    Time tendency of pprt field at west boundary (PASCAL/s)
c       pdtnb    Time tendency of pprt field at north boundary (PASCAL/s)
c       pdtsb    Time tendency of pprt field at south boundary (PASCAL/s)
c
c       qvdteb   Time tendency of qv field at east boundary (kg/(kg*s))
c       qvdtwb   Time tendency of qv field at west boundary (kg/(kg*s))
c       qvdtnb   Time tendency of qv field at north boundary (kg/(kg*s))
c       qvdtsb   Time tendency of qv field at south boundary (kg/(kg*s))
c
c       qcdteb   Time tendency of qc field at east boundary (kg/(kg*s))
c       qcdtwb   Time tendency of qc field at west boundary (kg/(kg*s))
c       qcdtnb   Time tendency of qc field at north boundary (kg/(kg*s))
c       qcdtsb   Time tendency of qc field at south boundary (kg/(kg*s))
c
c       qrdteb   Time tendency of qr field at east boundary (kg/(kg*s))
c       qrdtwb   Time tendency of qr field at west boundary (kg/(kg*s))
c       qrdtnb   Time tendency of qr field at north boundary (kg/(kg*s))
c       qrdtsb   Time tendency of qr field at south boundary (kg/(kg*s))
c
c       qidteb   Time tendency of qi field at east boundary (kg/(kg*s))
c       qidtwb   Time tendency of qi field at west boundary (kg/(kg*s))
c       qidtnb   Time tendency of qi field at north boundary (kg/(kg*s))
c       qidtsb   Time tendency of qi field at south boundary (kg/(kg*s))
c
c       qsdteb   Time tendency of qs field at east boundary (kg/(kg*s))
c       qsdtwb   Time tendency of qs field at west boundary (kg/(kg*s))
c       qsdtnb   Time tendency of qs field at north boundary (kg/(kg*s))
c       qsdtsb   Time tendency of qs field at south boundary (kg/(kg*s))
c
c       qhdteb   Time tendency of qh field at east boundary (kg/(kg*s))
c       qhdtwb   Time tendency of qh field at west boundary (kg/(kg*s))
c       qhdtnb   Time tendency of qh field at north boundary (kg/(kg*s))
c       qhdtsb   Time tendency of qh field at south boundary (kg/(kg*s))
c
c       ubar     Base state zonal velocity component (m/s)
c       vbar     Base state meridional velocity component (m/s)
c       ptbar    Base state potential temperature (K)
c       pbar     Base state pressure (Pascal)
c       rhobar   Base state density (kg/m**3)
c       qvbar    Base state water vapor specific humidity (kg/kg)
c
c       x        x coordinate of grid points in physical/comp. space (m)
c       y        y coordinate of grid points in physical/comp. space (m)
c       z        z coordinate of grid points in computational space (m)
c       zp       Vertical coordinate of grid points in physical space (m)
c       hterain  The height of terrain (m)
c
c       tstart   The time when the time integration starts, which is set to
c                the time of the restart data
c
c#######################################################################
c

c
c#######################################################################
c
c     Variable Declarations.
c
c#######################################################################
c
      implicit none 

      integer nt           ! The no. of t-levels of t-dependent arrays.
      integer tpast        ! Index of time level for the past time.
      integer tpresent     ! Index of time level for the present time.
      integer tfuture      ! Index of time level for the future time.

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

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

      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 ptdteb(ny,nz)        ! T-tendency of ptprt at e-boundary (K/s)
      real ptdtwb(ny,nz)        ! T-tendency of ptprt at w-boundary (K/s)
      real ptdtnb(nx,nz)        ! T-tendency of ptprt at n-boundary (K/s)
      real ptdtsb(nx,nz)        ! T-tendency of ptprt at s-boundary (K/s)

      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 qvdteb(ny,nz)        ! T-tendency of qv at e-boundary (kg/(kg*s))
      real qvdtwb(ny,nz)        ! T-tendency of qv at w-boundary (kg/(kg*s))
      real qvdtnb(nx,nz)        ! T-tendency of qv at n-boundary (kg/(kg*s))
      real qvdtsb(nx,nz)        ! T-tendency of qv at s-boundary (kg/(kg*s))

      real qcdteb(ny,nz)        ! T-tendency of qc at e-boundary (kg/(kg*s))
      real qcdtwb(ny,nz)        ! T-tendency of qc at w-boundary (kg/(kg*s))
      real qcdtnb(nx,nz)        ! T-tendency of qc at n-boundary (kg/(kg*s))
      real qcdtsb(nx,nz)        ! T-tendency of qc at s-boundary (kg/(kg*s))

      real qrdteb(ny,nz)        ! T-tendency of qr at e-boundary (kg/(kg*s))
      real qrdtwb(ny,nz)        ! T-tendency of qr at w-boundary (kg/(kg*s))
      real qrdtnb(nx,nz)        ! T-tendency of qr at n-boundary (kg/(kg*s))
      real qrdtsb(nx,nz)        ! T-tendency of qr at s-boundary (kg/(kg*s))

      real qidteb(ny,nz)        ! T-tendency of qi at e-boundary (kg/(kg*s))
      real qidtwb(ny,nz)        ! T-tendency of qi at w-boundary (kg/(kg*s))
      real qidtnb(nx,nz)        ! T-tendency of qi at n-boundary (kg/(kg*s))
      real qidtsb(nx,nz)        ! T-tendency of qi at s-boundary (kg/(kg*s))

      real qsdteb(ny,nz)        ! T-tendency of qs at e-boundary (kg/(kg*s))
      real qsdtwb(ny,nz)        ! T-tendency of qs at w-boundary (kg/(kg*s))
      real qsdtnb(nx,nz)        ! T-tendency of qs at n-boundary (kg/(kg*s))
      real qsdtsb(nx,nz)        ! T-tendency of qs at s-boundary (kg/(kg*s))

      real qhdteb(ny,nz)        ! T-tendency of qh at e-boundary (kg/(kg*s))
      real qhdtwb(ny,nz)        ! T-tendency of qh at w-boundary (kg/(kg*s))
      real qhdtnb(nx,nz)        ! T-tendency of qh at n-boundary (kg/(kg*s))
      real qhdtsb(nx,nz)        ! T-tendency of qh at s-boundary (kg/(kg*s))
c
      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 rhobar(nx,ny,nz)     ! Base state air density (kg/m**3)
      real qvbar (nx,ny,nz)     ! Base state water vapor specific humidity
                                ! (kg/kg)
c
      real x     (nx)           ! The x-coord. of the physical and
                                ! computational grid. Defined at u-point.
      real y     (ny)           ! The y-coord. of the physical and
                                ! computational grid. Defined at v-point.
      real z     (nz)           ! The z-coord. of the computational grid.
                                ! Defined at w-point on the staggered grid.
      real zp    (nx,ny,nz)     ! The physical height coordinate defined at
                                ! w-point of the staggered grid.

      real hterain(nx,ny)       ! The height of the terrain (m).
c
c#######################################################################
c
c     Misc. local variables:
c
c#######################################################################
c
      integer tim
      integer i, j, k
      integer nx1, ny1, nz1
      integer istat
      real    datatim
c
c#######################################################################
c
c     Include files: 
c
c#######################################################################
c
      include 'globcst.inc'

c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c
      open(unit=rstiunt,file=rstinf,form='unformatted',
     :     status='old',iostat=istat)

      IF( istat.ne.0) THEN

        write(6,'(/1x,a,i2,/1x,a/)')
     :       'Error occured when opening restart input file '//rstinf//
     :       ' using FORTRAN unit ',rstiunt,'Program stopped in RSTIN.'
        STOP

      ENDIF

      write(6,'(/1x,a,/1x,a,i2/)')
     :     'This is a restart run. Input was read from restart file ', 
     :      rstinf//' using fortran unit ',rstiunt
c
c
c#######################################################################
c
c     Write out the restart data:
c
c#######################################################################
c
      read(rstiunt,err=999) datatim 
 
      tstart = datatim

	write(6,'(a,f8.1)') ' Restart data is at time ',  datatim
      read(rstiunt,err=999) nx1,ny1,nz1


      IF((nx.ne.nx1).or.(ny.ne.ny1).or.(nz.ne.nz1)) 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=',
     :  nx1,', ny=',ny,', nz=',nz,' the model definitions were nx=',
     :  nx,' ny= ', ny, ' nz= ',nz
        write(6,'(a)') ' Job stopped in subroutine rstin.'
        STOP

      ENDIF

      tim = tpast

      read(rstiunt,err=999) (((u    (i,j,k,tim),i=1,nx),j=1,ny),k=1,nz)
      read(rstiunt,err=999) (((v    (i,j,k,tim),i=1,nx),j=1,ny),k=1,nz)
      read(rstiunt,err=999) (((w    (i,j,k,tim),i=1,nx),j=1,ny),k=1,nz)
      read(rstiunt,err=999) (((ptprt(i,j,k,tim),i=1,nx),j=1,ny),k=1,nz)
      read(rstiunt,err=999) (((pprt (i,j,k,tim),i=1,nx),j=1,ny),k=1,nz)
      read(rstiunt,err=999) (((qv   (i,j,k,tim),i=1,nx),j=1,ny),k=1,nz)
      read(rstiunt,err=999) (((qc   (i,j,k,tim),i=1,nx),j=1,ny),k=1,nz)
      read(rstiunt,err=999) (((qr   (i,j,k,tim),i=1,nx),j=1,ny),k=1,nz)

      IF( iceout. ne. 0) THEN 

        read(rstount,err=999) (((qi (i,j,k,tim),i=1,nx),j=1,ny),k=1,nz)
        read(rstount,err=999) (((qs (i,j,k,tim),i=1,nx),j=1,ny),k=1,nz)
        read(rstount,err=999) (((qh (i,j,k,tim),i=1,nx),j=1,ny),k=1,nz)

      ENDIF

      tim = tpresent
      
      read(rstiunt,err=999) (((u    (i,j,k,tim),i=1,nx),j=1,ny),k=1,nz)
      read(rstiunt,err=999) (((v    (i,j,k,tim),i=1,nx),j=1,ny),k=1,nz)
      read(rstiunt,err=999) (((w    (i,j,k,tim),i=1,nx),j=1,ny),k=1,nz)
      read(rstiunt,err=999) (((ptprt(i,j,k,tim),i=1,nx),j=1,ny),k=1,nz)
      read(rstiunt,err=999) (((pprt (i,j,k,tim),i=1,nx),j=1,ny),k=1,nz)
      read(rstiunt,err=999) (((qv   (i,j,k,tim),i=1,nx),j=1,ny),k=1,nz)
      read(rstiunt,err=999) (((qc   (i,j,k,tim),i=1,nx),j=1,ny),k=1,nz)
      read(rstiunt,err=999) (((qr   (i,j,k,tim),i=1,nx),j=1,ny),k=1,nz)

      IF( iceout. ne. 0) THEN 

        read(rstount,err=999) (((qi (i,j,k,tim),i=1,nx),j=1,ny),k=1,nz)
        read(rstount,err=999) (((qs (i,j,k,tim),i=1,nx),j=1,ny),k=1,nz)
        read(rstount,err=999) (((qh (i,j,k,tim),i=1,nx),j=1,ny),k=1,nz)

      ENDIF

      read(rstiunt,err=999) ((udteb (j,k),j=1,ny),k=1,nz)
      read(rstiunt,err=999) ((udtwb (j,k),j=1,ny),k=1,nz)
      read(rstiunt,err=999) ((udtnb (i,k),i=1,nx),k=1,nz)
      read(rstiunt,err=999) ((udtsb (i,k),i=1,nx),k=1,nz)

      read(rstiunt,err=999) ((vdteb (j,k),j=1,ny),k=1,nz)
      read(rstiunt,err=999) ((vdtwb (j,k),j=1,ny),k=1,nz)
      read(rstiunt,err=999) ((vdtnb (i,k),i=1,nx),k=1,nz)
      read(rstiunt,err=999) ((vdtsb (i,k),i=1,nx),k=1,nz)

      read(rstiunt,err=999) ((wdteb (j,k),j=1,ny),k=1,nz)
      read(rstiunt,err=999) ((wdtwb (j,k),j=1,ny),k=1,nz)
      read(rstiunt,err=999) ((wdtnb (i,k),i=1,nx),k=1,nz)
      read(rstiunt,err=999) ((wdtsb (i,k),i=1,nx),k=1,nz)

      read(rstiunt,err=999) ((ptdteb(j,k),j=1,ny),k=1,nz)
      read(rstiunt,err=999) ((ptdtwb(j,k),j=1,ny),k=1,nz)
      read(rstiunt,err=999) ((ptdtnb(i,k),i=1,nx),k=1,nz)
      read(rstiunt,err=999) ((ptdtsb(i,k),i=1,nx),k=1,nz)

      read(rstiunt,err=999) ((pdteb (j,k),j=1,ny),k=1,nz)
      read(rstiunt,err=999) ((pdtwb (j,k),j=1,ny),k=1,nz)
      read(rstiunt,err=999) ((pdtnb (i,k),i=1,nx),k=1,nz)
      read(rstiunt,err=999) ((pdtsb (i,k),i=1,nx),k=1,nz)

      read(rstiunt,err=999) ((qvdteb(j,k),j=1,ny),k=1,nz)
      read(rstiunt,err=999) ((qvdtwb(j,k),j=1,ny),k=1,nz)
      read(rstiunt,err=999) ((qvdtnb(i,k),i=1,nx),k=1,nz)
      read(rstiunt,err=999) ((qvdtsb(i,k),i=1,nx),k=1,nz)

      read(rstiunt,err=999) ((qcdteb(j,k),j=1,ny),k=1,nz)
      read(rstiunt,err=999) ((qcdtwb(j,k),j=1,ny),k=1,nz)
      read(rstiunt,err=999) ((qcdtnb(i,k),i=1,nx),k=1,nz)
      read(rstiunt,err=999) ((qcdtsb(i,k),i=1,nx),k=1,nz)

      read(rstiunt,err=999) ((qrdteb(j,k),j=1,ny),k=1,nz)
      read(rstiunt,err=999) ((qrdtwb(j,k),j=1,ny),k=1,nz)
      read(rstiunt,err=999) ((qrdtnb(i,k),i=1,nx),k=1,nz)
      read(rstiunt,err=999) ((qrdtsb(i,k),i=1,nx),k=1,nz)

      read(rstiunt,err=999) ((qidteb(j,k),j=1,ny),k=1,nz)
      read(rstiunt,err=999) ((qidtwb(j,k),j=1,ny),k=1,nz)
      read(rstiunt,err=999) ((qidtnb(i,k),i=1,nx),k=1,nz)
      read(rstiunt,err=999) ((qidtsb(i,k),i=1,nx),k=1,nz)

      read(rstiunt,err=999) ((qsdteb(j,k),j=1,ny),k=1,nz)
      read(rstiunt,err=999) ((qsdtwb(j,k),j=1,ny),k=1,nz)
      read(rstiunt,err=999) ((qsdtnb(i,k),i=1,nx),k=1,nz)
      read(rstiunt,err=999) ((qsdtsb(i,k),i=1,nx),k=1,nz)

      read(rstiunt,err=999) ((qhdteb(j,k),j=1,ny),k=1,nz)
      read(rstiunt,err=999) ((qhdtwb(j,k),j=1,ny),k=1,nz)
      read(rstiunt,err=999) ((qhdtnb(i,k),i=1,nx),k=1,nz)
      read(rstiunt,err=999) ((qhdtsb(i,k),i=1,nx),k=1,nz)

      read(rstiunt,err=999) (((ubar (i,j,k),i=1,nx),j=1,ny),k=1,nz)
      read(rstiunt,err=999) (((vbar (i,j,k),i=1,nx),j=1,ny),k=1,nz)
      read(rstiunt,err=999) (((ptbar(i,j,k),i=1,nx),j=1,ny),k=1,nz)
      read(rstiunt,err=999) (((pbar (i,j,k),i=1,nx),j=1,ny),k=1,nz)
      read(rstiunt,err=999) (((rhobar(i,j,k),i=1,nx),j=1,ny),k=1,nz)
      read(rstiunt,err=999) (((qvbar(i,j,k),i=1,nx),j=1,ny),k=1,nz)
c
c#######################################################################
c
c     Set the future values of variables to their current values.
c     This is done primarily for safety reasons.
c
c#######################################################################
c
      DO 100 k=1,nz
      DO 100 j=1,ny
      DO 100 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)
100   CONTINUE

      IF( iceout. ne. 0) THEN

        DO 200 k=1,nz
        DO 200 j=1,ny
        DO 200 i=1,nx
          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)
200     CONTINUE

      ENDIF

      close (unit=rstiunt)

      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)

      STOP

      END
