      SUBROUTINE wrt_rst (ng)
!
!svn $Id: wrt_rst.F 331 2009-03-12 00:34:51Z arango $
!================================================== Hernan G. Arango ===
!  Copyright (c) 2002-2009 The ROMS/TOMS Group                         !
!    Licensed under a MIT/X style license                              !
!    See License_ROMS.txt                                              !
!=======================================================================
!                                                                      !
!  This routine writes fields into restart NetCDF file.                !
!                                                                      !
!=======================================================================
!
      USE mod_param
      USE mod_parallel
      USE mod_grid
      USE mod_iounits
      USE mod_mixing
      USE mod_ncparam
      USE mod_netcdf
      USE mod_ocean
      USE mod_scalars
      USE mod_stepping
!
      USE nf_fwrite2d_mod, ONLY : nf_fwrite2d
      USE nf_fwrite3d_mod, ONLY : nf_fwrite3d
      USE nf_fwrite4d_mod, ONLY : nf_fwrite4d
!
      implicit none
!
!  Imported variable declarations.
!
      integer, intent(in) :: ng
!
!  Local variable declarations.
!
      integer :: LBi, UBi, LBj, UBj
      integer :: gfactor, gtype, i, itrc, status, varid
      integer :: ntmp(1)
      real(r8) :: scale
!
      LBi=LBOUND(GRID(ng)%h,DIM=1)
      UBi=UBOUND(GRID(ng)%h,DIM=1)
      LBj=LBOUND(GRID(ng)%h,DIM=2)
      UBj=UBOUND(GRID(ng)%h,DIM=2)
!
      SourceFile='wrt_rst.F'
!
!-----------------------------------------------------------------------
!  Write out restart fields.
!-----------------------------------------------------------------------
!
      IF (exit_flag.ne.NoError) RETURN
!
!  Set grid type factor to write full (gfactor=1) fields or water
!  points (gfactor=-1) fields only.
!
      gfactor=1
!
!  Set time record index.
!
      tRSTindx(ng)=tRSTindx(ng)+1
      NrecRST(ng)=NrecRST(ng)+1
!
!  If requested, set time index to recycle time records in restart
!  file.
!
      IF (LcycleRST(ng)) THEN
        tRSTindx(ng)=MOD(tRSTindx(ng)-1,2)+1
      END IF
!
!  Write out time-stepping indices.
!
      ntmp(1)=1+MOD((iic(ng)-1)-ntstart(ng),2)
      CALL netcdf_put_ivar (ng, iNLM, RSTname(ng), 'nstp',              &
     &                      ntmp, (/tRSTindx(ng)/), (/1/),              &
     &                      ncid = ncRSTid(ng))
      IF (exit_flag.ne.NoError) RETURN
      CALL netcdf_put_ivar (ng, iNLM, RSTname(ng), 'nrhs',              &
     &                      ntmp, (/tRSTindx(ng)/), (/1/),              &
     &                      ncid = ncRSTid(ng))
      IF (exit_flag.ne.NoError) RETURN
      ntmp(1)=3-ntmp(1)
      CALL netcdf_put_ivar (ng, iNLM, RSTname(ng), 'nnew',              &
     &                      ntmp, (/tRSTindx(ng)/), (/1/),              &
     &                      ncid = ncRSTid(ng))
      IF (exit_flag.ne.NoError) RETURN
      CALL netcdf_put_ivar (ng, iNLM, RSTname(ng), 'kstp',              &
     &                      kstp(ng:), (/tRSTindx(ng)/), (/1/),         &
     &                      ncid = ncRSTid(ng))
      IF (exit_flag.ne.NoError) RETURN
      CALL netcdf_put_ivar (ng, iNLM, RSTname(ng), 'krhs',              &
     &                      krhs(ng:), (/tRSTindx(ng)/), (/1/),         &
     &                      ncid = ncRSTid(ng))
      IF (exit_flag.ne.NoError) RETURN
      CALL netcdf_put_ivar (ng, iNLM, RSTname(ng), 'knew',              &
     &                      knew(ng:), (/tRSTindx(ng)/), (/1/),         &
     &                      ncid = ncRSTid(ng))
      IF (exit_flag.ne.NoError) RETURN
!
!  Write out model time (s).
!
      CALL netcdf_put_fvar (ng, iNLM, RSTname(ng),                      &
     &                      TRIM(Vname(idtime,ng)), time(ng:),          &
     &                      (/tRSTindx(ng)/), (/1/),                    &
     &                      ncid = ncRSTid(ng),                         &
     &                      varid = rstVid(idtime,ng))
      IF (exit_flag.ne.NoError) RETURN
!
!  Write out free-surface (m).
!
      scale=1.0_r8
      gtype=gfactor*r3dvar
      status=nf_fwrite3d(ng, iNLM, ncRSTid(ng), rstVid(idFsur,ng),      &
     &                   tRSTindx(ng), gtype,                           &
     &                   LBi, UBi, LBj, UBj, 1, 3, scale,               &
     &                   GRID(ng) % rmask,                              &
     &                   OCEAN(ng) % zeta)
      IF (status.ne.nf90_noerr) THEN
        IF (Master) THEN
          WRITE (stdout,10) TRIM(Vname(1,idFsur)), tRSTindx(ng)
        END IF
        exit_flag=3
        ioerror=status
        RETURN
      END IF
!
!  Write out RHS of free-surface equation.
!
      scale=1.0_r8
      gtype=gfactor*r3dvar
      status=nf_fwrite3d(ng, iNLM, ncRSTid(ng), rstVid(idRzet,ng),      &
     &                   tRSTindx(ng), gtype,                           &
     &                   LBi, UBi, LBj, UBj, 1, 2, scale,               &
     &                   GRID(ng) % rmask,                              &
     &                   OCEAN(ng) % rzeta)
      IF (status.ne.nf90_noerr) THEN
        IF (Master) THEN
          WRITE (stdout,10) TRIM(Vname(1,idRzet)), tRSTindx(ng)
        END IF
        exit_flag=3
        ioerror=status
        RETURN
      END IF
!
!  Write out 2D momentum component (m/s) in the XI-direction.
!
      scale=1.0_r8
      gtype=gfactor*u3dvar
      status=nf_fwrite3d(ng, iNLM, ncRSTid(ng), rstVid(idUbar,ng),      &
     &                   tRSTindx(ng), gtype,                           &
     &                   LBi, UBi, LBj, UBj, 1, 3, scale,               &
     &                   GRID(ng) % umask,                              &
     &                   OCEAN(ng) % ubar)
      IF (status.ne.nf90_noerr) THEN
        IF (Master) THEN
          WRITE (stdout,10) TRIM(Vname(1,idUbar)), tRSTindx(ng)
        END IF
        exit_flag=3
        ioerror=status
        RETURN
      END IF
!
!  Write out RHS of 2D momentum equation in the XI-direction.
!
      scale=1.0_r8
      gtype=gfactor*u3dvar
      status=nf_fwrite3d(ng, iNLM, ncRSTid(ng), rstVid(idRu2d,ng),      &
     &                   tRSTindx(ng), gtype,                           &
     &                   LBi, UBi, LBj, UBj, 1, 2, scale,               &
     &                   GRID(ng) % umask,                              &
     &                   OCEAN(ng) % rubar)
      IF (status.ne.nf90_noerr) THEN
        IF (Master) THEN
          WRITE (stdout,10) TRIM(Vname(1,idRu2d)), tRSTindx(ng)
        END IF
        exit_flag=3
        ioerror=status
        RETURN
      END IF
!
!  Write out 2D momentum component (m/s) in the ETA-direction.
!
      scale=1.0_r8
      gtype=gfactor*v3dvar
      status=nf_fwrite3d(ng, iNLM, ncRSTid(ng), rstVid(idVbar,ng),      &
     &                   tRSTindx(ng), gtype,                           &
     &                   LBi, UBi, LBj, UBj, 1, 3, scale,               &
     &                   GRID(ng) % vmask,                              &
     &                   OCEAN(ng) % vbar)
      IF (status.ne.nf90_noerr) THEN
        IF (Master) THEN
          WRITE (stdout,10) TRIM(Vname(1,idVbar)), tRSTindx(ng)
        END IF
        exit_flag=3
        ioerror=status
        RETURN
      END IF
!
!  Write out RHS of 2D momentum equation in the ETA-direction.
!
      scale=1.0_r8
      gtype=gfactor*v3dvar
      status=nf_fwrite3d(ng, iNLM, ncRSTid(ng), rstVid(idRv2d,ng),      &
     &                   tRSTindx(ng), gtype,                           &
     &                   LBi, UBi, LBj, UBj, 1, 2, scale,               &
     &                   GRID(ng) % vmask,                              &
     &                   OCEAN(ng) % rvbar)
      IF (status.ne.nf90_noerr) THEN
        IF (Master) THEN
          WRITE (stdout,10) TRIM(Vname(1,idRv2d)), tRSTindx(ng)
        END IF
        exit_flag=3
        ioerror=status
        RETURN
      END IF
!
!  Write out 3D momentum component (m/s) in the XI-direction.
!
      scale=1.0_r8
      gtype=gfactor*u3dvar
      status=nf_fwrite4d(ng, iNLM, ncRSTid(ng), rstVid(idUvel,ng),      &
     &                   tRSTindx(ng), gtype,                           &
     &                   LBi, UBi, LBj, UBj, 1, N(ng), 1, 2, scale,     &
     &                   GRID(ng) % umask,                              &
     &                   OCEAN(ng) % u)
      IF (status.ne.nf90_noerr) THEN
        IF (Master) THEN
          WRITE (stdout,10) TRIM(Vname(1,idUvel)), tRSTindx(ng)
        END IF
        exit_flag=3
        ioerror=status
        RETURN
      END IF
!
!  Write out RHS of 3D momentum equation in the XI-direction.
!
      scale=1.0_r8
      gtype=gfactor*u3dvar
      status=nf_fwrite4d(ng, iNLM, ncRSTid(ng), rstVid(idRu3d,ng),      &
     &                   tRSTindx(ng), gtype,                           &
     &                   LBi, UBi, LBj, UBj, 0, N(ng), 1, 2, scale,     &
     &                   GRID(ng) % umask,                              &
     &                   OCEAN(ng) % ru)
      IF (status.ne.nf90_noerr) THEN
        IF (Master) THEN
          WRITE (stdout,10) TRIM(Vname(1,idRu3d)), tRSTindx(ng)
        END IF
        exit_flag=3
        ioerror=status
        RETURN
      END IF
!
!  Write out momentum component (m/s) in the ETA-direction.
!
      scale=1.0_r8
      gtype=gfactor*v3dvar
      status=nf_fwrite4d(ng, iNLM, ncRSTid(ng), rstVid(idVvel,ng),      &
     &                   tRSTindx(ng), gtype,                           &
     &                   LBi, UBi, LBj, UBj, 1, N(ng), 1, 2, scale,     &
     &                   GRID(ng) % vmask,                              &
     &                   OCEAN(ng) % v)
      IF (status.ne.nf90_noerr) THEN
        IF (Master) THEN
          WRITE (stdout,10) TRIM(Vname(1,idVvel)), tRSTindx(ng)
        END IF
        exit_flag=3
        ioerror=status
        RETURN
      END IF
!
!  Write out RHS of 3D momentum equation in the ETA-direction.
!
      scale=1.0_r8
      gtype=gfactor*v3dvar
      status=nf_fwrite4d(ng, iNLM, ncRSTid(ng), rstVid(idRv3d,ng),      &
     &                   tRSTindx(ng), gtype,                           &
     &                   LBi, UBi, LBj, UBj, 0, N(ng), 1, 2, scale,     &
     &                   GRID(ng) % vmask,                              &
     &                   OCEAN(ng) % rv)
      IF (status.ne.nf90_noerr) THEN
        IF (Master) THEN
          WRITE (stdout,10) TRIM(Vname(1,idRv3d)), tRSTindx(ng)
        END IF
        exit_flag=3
        ioerror=status
        RETURN
      END IF
!
!  Write out tracer type variables.
!
      DO itrc=1,NT(ng)
        scale=1.0_r8
        gtype=gfactor*r3dvar
        status=nf_fwrite4d(ng, iNLM, ncRSTid(ng), rstTid(itrc,ng),      &
     &                     tRSTindx(ng), gtype,                         &
     &                     LBi, UBi, LBj, UBj, 1, N(ng), 1, 2, scale,   &
     &                     GRID(ng) % rmask,                            &
     &                     OCEAN(ng) % t(:,:,:,:,itrc))
        IF (status.ne.nf90_noerr) THEN
          IF (Master) THEN
            WRITE (stdout,10) TRIM(Vname(1,idTvar(itrc))), tRSTindx(ng)
          END IF
          exit_flag=3
          ioerror=status
          RETURN
        END IF
      END DO
!
!  Write out density anomaly.
!
      scale=1.0_r8
      gtype=gfactor*r3dvar
      status=nf_fwrite3d(ng, iNLM, ncRSTid(ng), rstVid(idDano,ng),      &
     &                   tRSTindx(ng), gtype,                           &
     &                   LBi, UBi, LBj, UBj, 1, N(ng), scale,           &
     &                   GRID(ng) % rmask,                              &
     &                   OCEAN(ng) % rho)
      IF (status.ne.nf90_noerr) THEN
        IF (Master) THEN
          WRITE (stdout,10) TRIM(Vname(1,idDano)), tRSTindx(ng)
        END IF
        exit_flag=3
        ioerror=status
        RETURN
      END IF
!
!  Write out vertical viscosity coefficient.
!
      scale=1.0_r8
      gtype=gfactor*w3dvar
      status=nf_fwrite3d(ng, iNLM, ncRSTid(ng), rstVid(idVvis,ng),      &
     &                   tRSTindx(ng), gtype,                           &
     &                   LBi, UBi, LBj, UBj, 0, N(ng), scale,           &
     &                   GRID(ng) % rmask,                              &
     &                   MIXING(ng) % Akv)
      IF (status.ne.nf90_noerr) THEN
        IF (Master) THEN
          WRITE (stdout,10) TRIM(Vname(1,idVvis)), tRSTindx(ng)
        END IF
        exit_flag=3
        ioerror=status
        RETURN
      END IF
!
!  Write out vertical diffusion coefficient for potential temperature.
!
      scale=1.0_r8
      gtype=gfactor*w3dvar
      status=nf_fwrite3d(ng, iNLM, ncRSTid(ng), rstVid(idTdif,ng),      &
     &                   tRSTindx(ng), gtype,                           &
     &                   LBi, UBi, LBj, UBj, 0, N(ng), scale,           &
     &                   GRID(ng) % rmask,                              &
     &                   MIXING(ng) % Akt(:,:,:,itemp))
      IF (status.ne.nf90_noerr) THEN
        IF (Master) THEN
          WRITE (stdout,10) TRIM(Vname(1,idTdif)), tRSTindx(ng)
        END IF
        exit_flag=3
        ioerror=status
        RETURN
      END IF
!
!  Write out vertical diffusion coefficient for salinity.
!
      scale=1.0_r8
      gtype=gfactor*w3dvar
      status=nf_fwrite3d(ng, iNLM, ncRSTid(ng), rstVid(idSdif,ng),      &
     &                   tRSTindx(ng), gtype,                           &
     &                   LBi, UBi, LBj, UBj, 0, N(ng), scale,           &
     &                   GRID(ng) % rmask,                              &
     &                   MIXING(ng) % Akt(:,:,:,isalt))
      IF (status.ne.nf90_noerr) THEN
        IF (Master) THEN
          WRITE (stdout,10) TRIM(Vname(1,idSdif)), tRSTindx(ng)
        END IF
        exit_flag=3
        ioerror=status
        RETURN
      END IF
!
!  Write out turbulent kinetic energy.
!
      scale=1.0_r8
      gtype=gfactor*w3dvar
      status=nf_fwrite4d(ng, iNLM, ncRSTid(ng), rstVid(idMtke,ng),      &
     &                   tRSTindx(ng), gtype,                           &
     &                   LBi, UBi, LBj, UBj, 0, N(ng), 1, 2, scale,     &
     &                   GRID(ng) % rmask,                              &
     &                   MIXING(ng) % tke)
      IF (status.ne.nf90_noerr) THEN
        IF (Master) THEN
          WRITE (stdout,10) TRIM(Vname(1,idMtke)), tRSTindx(ng)
        END IF
        exit_flag=3
        ioerror=status
        RETURN
      END IF
!
!  Define turbulent kinetic energy time length scale.
!
      scale=1.0_r8
      gtype=gfactor*w3dvar
      status=nf_fwrite4d(ng, iNLM, ncRSTid(ng), rstVid(idMtls,ng),      &
     &                   tRSTindx(ng), gtype,                           &
     &                   LBi, UBi, LBj, UBj, 0, N(ng), 1, 2, scale,     &
     &                   GRID(ng) % rmask,                              &
     &                   MIXING(ng) % gls)
      IF (status.ne.nf90_noerr) THEN
        IF (Master) THEN
          WRITE (stdout,10) TRIM(Vname(1,idMtls)), tRSTindx(ng)
        END IF
        exit_flag=3
        ioerror=status
        RETURN
      END IF
!
!  Define vertical mixing turbulent length scale.
!
      scale=1.0_r8
      gtype=gfactor*w3dvar
      status=nf_fwrite3d(ng, iNLM, ncRSTid(ng), rstVid(idVmLS,ng),      &
     &                   tRSTindx(ng), gtype,                           &
     &                   LBi, UBi, LBj, UBj, 0, N(ng), scale,           &
     &                   GRID(ng) % rmask,                              &
     &                   MIXING(ng) % Lscale)
      IF (status.ne.nf90_noerr) THEN
        IF (Master) THEN
          WRITE (stdout,10) TRIM(Vname(1,idVmLS)), tRSTindx(ng)
        END IF
        exit_flag=3
        ioerror=status
        RETURN
      END IF
!
!  Define turbulent kinetic energy vertical diffusion coefficient.
!
      scale=1.0_r8
      gtype=gfactor*w3dvar
      status=nf_fwrite3d(ng, iNLM, ncRSTid(ng), rstVid(idVmKK,ng),      &
     &                   tRSTindx(ng), gtype,                           &
     &                   LBi, UBi, LBj, UBj, 0, N(ng), scale,           &
     &                   GRID(ng) % rmask,                              &
     &                   MIXING(ng) % Akk)
      IF (status.ne.nf90_noerr) THEN
        IF (Master) THEN
          WRITE (stdout,10) TRIM(Vname(1,idVmKK)), tRSTindx(ng)
        END IF
        exit_flag=3
        ioerror=status
        RETURN
      END IF
!
!  Define turbulent length scale vertical diffusion coefficient.
!
      scale=1.0_r8
      gtype=gfactor*w3dvar
      status=nf_fwrite3d(ng, iNLM, ncRSTid(ng), rstVid(idVmKP,ng),      &
     &                   tRSTindx(ng), gtype,                           &
     &                   LBi, UBi, LBj, UBj, 0, N(ng), scale,           &
     &                   GRID(ng) % rmask,                              &
     &                   MIXING(ng) % Akp)
      IF (status.ne.nf90_noerr) THEN
        IF (Master) THEN
          WRITE (stdout,10) TRIM(Vname(1,idVmKP)), tRSTindx(ng)
        END IF
        exit_flag=3
        ioerror=status
        RETURN
      END IF
!
!-----------------------------------------------------------------------
!  Synchronize restart NetCDF file to disk.
!-----------------------------------------------------------------------
!
      CALL netcdf_sync (ng, iNLM, RSTname(ng), ncRSTid(ng))
      IF (exit_flag.ne.NoError) RETURN
      IF (Master) WRITE (stdout,20) kstp(ng), nrhs(ng), tRSTindx(ng)
!
  10  FORMAT (/,' WRT_RST - error while writing variable: ',a,/,11x,    &
     &        'into restart NetCDF file for time record: ',i4)
  20  FORMAT (6x,'WRT_RST   - wrote re-start fields (Index=', i1,       &
     &        ',',i1,') into time record = ',i7.7)
      RETURN
      END SUBROUTINE wrt_rst
