      SUBROUTINE wrt_quick (ng, tile)
!
!svn $Id$
!================================================== Hernan G. Arango ===
!  Copyright (c) 2002-2020 The ROMS/TOMS Group                         !
!    Licensed under a MIT/X style license                              !
!    See License_ROMS.txt                                              !
!=======================================================================
!                                                                      !
!  This routine writes requested model fields at requested levels      !
!  into ROMS QUICKSAVE NetCDF file.                                    !
!                                                                      !
!  Notice that only momentum is affected by the full time-averaged     !
!  masks.  If applicable, these mask contains information about        !
!  river runoff and time-dependent wetting and drying variations.      !
!                                                                      !
!=======================================================================
!
      USE mod_param
      USE mod_parallel
      USE mod_coupling
      USE mod_forces
      USE mod_grid
      USE mod_ice
      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 omega_mod,       ONLY : scale_omega
      USE uv_rotate_mod,   ONLY : uv_rotate2d
      USE uv_rotate_mod,   ONLY : uv_rotate3d
      USE strings_mod,     ONLY : FoundError
!
      implicit none
!
!  Imported variable declarations.
!
      integer, intent(in) :: ng, tile
!
!  Local variable declarations.
!
      integer :: LBi, UBi, LBj, UBj
      integer :: Fcount, gfactor, gtype, status
      integer :: i, itrc, j, k
      real(r8) :: scale
      real(r8), allocatable :: Ur2d(:,:)
      real(r8), allocatable :: Vr2d(:,:)
      real(r8), allocatable :: EBF(:,:)
      real(r8), allocatable :: dRdx(:,:), dRde(:,:)
      real(r8), allocatable :: dx(:,:), dy(:,:), frho(:,:)
      real(r8), allocatable :: tauXpsi(:,:), tauYpsi(:,:)
      real(r8) :: cff
      real(r8), allocatable :: strain(:,:)
      real(r8), allocatable :: dUdx(:,:), dUde(:,:)
      real(r8), allocatable :: dVdx(:,:), dVde(:,:)
      real(r8), allocatable :: rvor(:,:)
      real(r8), allocatable :: rdVdx(:,:), rdUde(:,:)
      real(r8), allocatable :: dz_r(:,:), dz_u(:,:), dz_v(:,:)
      integer :: k1,k2
      real(r8) :: orho0, fomn_p
      real(r8) :: dRde_pr, dRdx_pr, dRdz_pr, dUdz_pr, dVdz_pr
      real(r8), allocatable :: pvor(:,:),pvork(:,:)
      real(r8), allocatable :: pvori(:,:), pvorj(:,:)
      real(r8), allocatable :: dUdz(:,:), dVdz(:,:), dRdz(:,:)
      real(r8), allocatable :: Ur3d(:,:,:)
      real(r8), allocatable :: Vr3d(:,:,:)
      real(r8), allocatable :: Wr3d(:,:,:)
!
!-----------------------------------------------------------------------
!  Set lower and upper tile bounds and staggered variables bounds for
!  this horizontal domain partition.  Notice that if tile=-1, it will
!  set the values for the global grid.
!-----------------------------------------------------------------------
!
      integer :: Istr, IstrB, IstrP, IstrR, IstrT, IstrM, IstrU
      integer :: Iend, IendB, IendP, IendR, IendT
      integer :: Jstr, JstrB, JstrP, JstrR, JstrT, JstrM, JstrV
      integer :: Jend, JendB, JendP, JendR, JendT
      integer :: Istrm3, Istrm2, Istrm1, IstrUm2, IstrUm1
      integer :: Iendp1, Iendp2, Iendp2i, Iendp3
      integer :: Jstrm3, Jstrm2, Jstrm1, JstrVm2, JstrVm1
      integer :: Jendp1, Jendp2, Jendp2i, Jendp3
!
      Istr   =BOUNDS(ng) % Istr   (tile)
      IstrB  =BOUNDS(ng) % IstrB  (tile)
      IstrM  =BOUNDS(ng) % IstrM  (tile)
      IstrP  =BOUNDS(ng) % IstrP  (tile)
      IstrR  =BOUNDS(ng) % IstrR  (tile)
      IstrT  =BOUNDS(ng) % IstrT  (tile)
      IstrU  =BOUNDS(ng) % IstrU  (tile)
      Iend   =BOUNDS(ng) % Iend   (tile)
      IendB  =BOUNDS(ng) % IendB  (tile)
      IendP  =BOUNDS(ng) % IendP  (tile)
      IendR  =BOUNDS(ng) % IendR  (tile)
      IendT  =BOUNDS(ng) % IendT  (tile)
      Jstr   =BOUNDS(ng) % Jstr   (tile)
      JstrB  =BOUNDS(ng) % JstrB  (tile)
      JstrM  =BOUNDS(ng) % JstrM  (tile)
      JstrP  =BOUNDS(ng) % JstrP  (tile)
      JstrR  =BOUNDS(ng) % JstrR  (tile)
      JstrT  =BOUNDS(ng) % JstrT  (tile)
      JstrV  =BOUNDS(ng) % JstrV  (tile)
      Jend   =BOUNDS(ng) % Jend   (tile)
      JendB  =BOUNDS(ng) % JendB  (tile)
      JendP  =BOUNDS(ng) % JendP  (tile)
      JendR  =BOUNDS(ng) % JendR  (tile)
      JendT  =BOUNDS(ng) % JendT  (tile)
!
      Istrm3 =BOUNDS(ng) % Istrm3 (tile)            ! Istr-3
      Istrm2 =BOUNDS(ng) % Istrm2 (tile)            ! Istr-2
      Istrm1 =BOUNDS(ng) % Istrm1 (tile)            ! Istr-1
      IstrUm2=BOUNDS(ng) % IstrUm2(tile)            ! IstrU-2
      IstrUm1=BOUNDS(ng) % IstrUm1(tile)            ! IstrU-1
      Iendp1 =BOUNDS(ng) % Iendp1 (tile)            ! Iend+1
      Iendp2 =BOUNDS(ng) % Iendp2 (tile)            ! Iend+2
      Iendp2i=BOUNDS(ng) % Iendp2i(tile)            ! Iend+2 interior
      Iendp3 =BOUNDS(ng) % Iendp3 (tile)            ! Iend+3
      Jstrm3 =BOUNDS(ng) % Jstrm3 (tile)            ! Jstr-3
      Jstrm2 =BOUNDS(ng) % Jstrm2 (tile)            ! Jstr-2
      Jstrm1 =BOUNDS(ng) % Jstrm1 (tile)            ! Jstr-1
      JstrVm2=BOUNDS(ng) % JstrVm2(tile)            ! JstrV-2
      JstrVm1=BOUNDS(ng) % JstrVm1(tile)            ! JstrV-1
      Jendp1 =BOUNDS(ng) % Jendp1 (tile)            ! Jend+1
      Jendp2 =BOUNDS(ng) % Jendp2 (tile)            ! Jend+2
      Jendp2i=BOUNDS(ng) % Jendp2i(tile)            ! Jend+2 interior
      Jendp3 =BOUNDS(ng) % Jendp3 (tile)            ! Jend+3
!
      SourceFile="ROMS/Utility/wrt_quick.F"
!
      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)
!
!-----------------------------------------------------------------------
!  Write out quicksave fields.
!-----------------------------------------------------------------------
!
      IF (FoundError(exit_flag, NoError, 117,                      &
     &               "ROMS/Utility/wrt_quick.F")) RETURN
!
!  Set grid type factor to write full (gfactor=1) fields or water
!  points (gfactor=-1) fields only.
!
      gfactor=1
!
!  Set time record index.
!
      QCK(ng)%Rindex=QCK(ng)%Rindex+1
      Fcount=QCK(ng)%load
      QCK(ng)%Nrec(Fcount)=QCK(ng)%Nrec(Fcount)+1
!
!  Write out model time (s).
!
      CALL netcdf_put_fvar (ng, iNLM, QCK(ng)%name,                     &
     &                      TRIM(Vname(idtime,ng)), time(ng:),          &
     &                      (/QCK(ng)%Rindex/), (/1/),                  &
     &                      ncid = QCK(ng)%ncid,                        &
     &                      varid = QCK(ng)%Vid(idtime))
      IF (FoundError(exit_flag, NoError, 142,                      &
     &               "ROMS/Utility/wrt_quick.F")) RETURN
!
!  Write time-varying depths of RHO-points.
!
      IF (Qout(idpthR,ng)) THEN
        scale=1.0_r8
        gtype=gfactor*r3dvar
        status=nf_fwrite3d(ng, iNLM, QCK(ng)%ncid, QCK(ng)%Vid(idpthR), &
     &                     QCK(ng)%Rindex, gtype,                       &
     &                     LBi, UBi, LBj, UBj, 1, N(ng), scale,         &
     &                     GRID(ng) % rmask,                            &
     &                     GRID(ng) % z_r)
        IF (FoundError(status, nf90_noerr, 275,                    &
     &                 "ROMS/Utility/wrt_quick.F")) THEN
          IF (Master) THEN
            WRITE (stdout,10) TRIM(Vname(1,idpthR)), QCK(ng)%Rindex
          END IF
          exit_flag=3
          ioerror=status
          RETURN
        END IF
      END IF
!
!  Write time-varying depths of U-points.
!
      IF (Qout(idpthU,ng)) THEN
        scale=1.0_r8
        gtype=gfactor*u3dvar
        DO k=1,N(ng)
          DO j=Jstr-1,Jend+1
            DO i=IstrU-1,Iend+1
              GRID(ng)%z_v(i,j,k)=0.5_r8*(GRID(ng)%z_r(i-1,j,k)+        &
     &                                    GRID(ng)%z_r(i  ,j,k))
            END DO
          END DO
        END DO
        status=nf_fwrite3d(ng, iNLM, QCK(ng)%ncid, QCK(ng)%Vid(idpthU), &
     &                     QCK(ng)%Rindex, gtype,                       &
     &                     LBi, UBi, LBj, UBj, 1, N(ng), scale,         &
     &                     GRID(ng) % umask,                            &
     &                     GRID(ng) % z_v)
        IF (FoundError(status, nf90_noerr, 306,                    &
     &                 "ROMS/Utility/wrt_quick.F")) THEN
          IF (Master) THEN
            WRITE (stdout,10) TRIM(Vname(1,idpthU)), QCK(ng)%Rindex
          END IF
          exit_flag=3
          ioerror=status
          RETURN
        END IF
      END IF
!
!  Write time-varying depths of V-points.
!
      IF (Qout(idpthV,ng)) THEN
        scale=1.0_r8
        gtype=gfactor*v3dvar
        DO k=1,N(ng)
          DO j=JstrV-1,Jend+1
            DO i=Istr-1,Iend+1
              GRID(ng)%z_v(i,j,k)=0.5_r8*(GRID(ng)%z_r(i,j-1,k)+       &
     &                                    GRID(ng)%z_r(i,j  ,k))
            END DO
          END DO
        END DO
        status=nf_fwrite3d(ng, iNLM, QCK(ng)%ncid, QCK(ng)%Vid(idpthV), &
     &                     QCK(ng)%Rindex, gtype,                       &
     &                     LBi, UBi, LBj, UBj, 1, N(ng), scale,         &
     &                     GRID(ng) % vmask,                            &
     &                     GRID(ng) % z_v)
        IF (FoundError(status, nf90_noerr, 337,                    &
     &                 "ROMS/Utility/wrt_quick.F")) THEN
          IF (Master) THEN
            WRITE (stdout,10) TRIM(Vname(1,idpthV)), QCK(ng)%Rindex
          END IF
          exit_flag=3
          ioerror=status
          RETURN
        END IF
      END IF
!
!  Write time-varying depths of W-points.
!
      IF (Qout(idpthW,ng)) THEN
        scale=1.0_r8
        gtype=gfactor*w3dvar
        status=nf_fwrite3d(ng, iNLM, QCK(ng)%ncid, QCK(ng)%Vid(idpthW), &
     &                     QCK(ng)%Rindex, gtype,                       &
     &                     LBi, UBi, LBj, UBj, 0, N(ng), scale,         &
     &                     GRID(ng) % rmask,                            &
     &                     GRID(ng) % z_w)
        IF (FoundError(status, nf90_noerr, 360,                    &
     &                 "ROMS/Utility/wrt_quick.F")) THEN
          IF (Master) THEN
            WRITE (stdout,10) TRIM(Vname(1,idpthW)), QCK(ng)%Rindex
          END IF
          exit_flag=3
          ioerror=status
          RETURN
        END IF
      END IF
!
!  Write out free-surface (m)
!
      IF (Qout(idFsur,ng)) THEN
        scale=1.0_r8
        gtype=gfactor*r2dvar
        status=nf_fwrite2d(ng, iNLM, QCK(ng)%ncid, QCK(ng)%Vid(idFsur), &
     &                     QCK(ng)%Rindex, gtype,                       &
     &                     LBi, UBi, LBj, UBj, scale,                   &
     &                     GRID(ng) % rmask,                            &
     &                     OCEAN(ng) % zeta(:,:,kstp(ng)))
        IF (FoundError(status, nf90_noerr, 389,                    &
     &                 "ROMS/Utility/wrt_quick.F")) THEN
          IF (Master) THEN
            WRITE (stdout,10) TRIM(Vname(1,idFsur)), QCK(ng)%Rindex
          END IF
          exit_flag=3
          ioerror=status
          RETURN
        END IF
      END IF
!
!  Write out 2D U-momentum component (m/s).
!
      IF (Qout(idUbar,ng)) THEN
        scale=1.0_r8
        gtype=gfactor*u2dvar
        status=nf_fwrite2d(ng, iNLM, QCK(ng)%ncid, QCK(ng)%Vid(idUbar), &
     &                     QCK(ng)%Rindex, gtype,                       &
     &                     LBi, UBi, LBj, UBj, scale,                   &
     &                     GRID(ng) % umask_full,                       &
     &                     OCEAN(ng) % ubar(:,:,kstp(ng)))
        IF (FoundError(status, nf90_noerr, 412,                    &
     &                 "ROMS/Utility/wrt_quick.F")) THEN
          IF (Master) THEN
            WRITE (stdout,10) TRIM(Vname(1,idUbar)), QCK(ng)%Rindex
          END IF
          exit_flag=3
          ioerror=status
          RETURN
        END IF
      END IF
!
!  Write out 2D V-momentum component (m/s).
!
      IF (Qout(idVbar,ng)) THEN
        scale=1.0_r8
        gtype=gfactor*v2dvar
        status=nf_fwrite2d(ng, iNLM, QCK(ng)%ncid, QCK(ng)%Vid(idVbar), &
     &                     QCK(ng)%Rindex, gtype,                       &
     &                     LBi, UBi, LBj, UBj, scale,                   &
     &                     GRID(ng) % vmask_full,                       &
     &                     OCEAN(ng) % vbar(:,:,kstp(ng)))
        IF (FoundError(status, nf90_noerr, 435,                    &
     &                 "ROMS/Utility/wrt_quick.F")) THEN
          IF (Master) THEN
            WRITE (stdout,10) TRIM(Vname(1,idVbar)), QCK(ng)%Rindex
          END IF
          exit_flag=3
          ioerror=status
          RETURN
        END IF
      END IF
!
!  Write out Ekman buoyancy flux at PSI-points.
!
!psl  IF (Qout(idEBF,ng)) THEN !psl20210702: idEBF absent from varinfo.
!       IF (.not.allocated(EBF)) THEN
!         allocate (EBF(LBi:UBi,LBj:UBj))
!           EBF(LBi:UBi,LBj:UBj)=0.0_r8
!         allocate (dRdx(LBi:UBi,LBj:UBj))
!           dRdx(LBi:UBi,LBj:UBj)=0.0_r8
!         allocate (dRde(LBi:UBi,LBj:UBj))
!           dRde(LBi:UBi,LBj:UBj)=0.0_r8
!         allocate (tauXpsi(LBi:UBi,LBj:UBj))
!           tauXpsi(LBi:UBi,LBj:UBj)=0.0_r8
!         allocate (tauYpsi(LBi:UBi,LBj:UBj))
!           tauYpsi(LBi:UBi,LBj:UBj)=0.0_r8
!         allocate (dx(LBi:UBi,LBj:UBj))
!           dx(LBi:UBi,LBj:UBj)=0.0_r8
!         allocate (dy(LBi:UBi,LBj:UBj))
!           dy(LBi:UBi,LBj:UBj)=0.0_r8
!         allocate (frho(LBi:UBi,LBj:UBj))
!           frho(LBi:UBi,LBj:UBj)=0.0_r8
!       END IF
!       DO i=LBi+1,UBi-1
!         DO j=LBj+1,UBj-1
! Now do Ekman buoyancy flux
! Subtract 1 in the y direction to get sustr
!            EBF(i,j) = FORCES(ng) % sustr(i,j-1)/rho0
! Subtract 1 in the x direction to get svstr
!            EBF(i,j) = FORCES(ng) % svstr(i-1,j)/rho0
! Subtract 1 in both the x direction and the y direction to get rho, same with pm and pn
!            EBF(i,j) = OCEAN(ng) % rho(i-1,j-1,N(ng))
! Calculate tau_x on psi grid
!            tauXpsi(i,j) =                                           &
!                  (  FORCES(ng) % sustr(i,j-1)                       &
!                   + FORCES(ng) % sustr(i,j)   ) * rho0 /2.
!             EBF(i,j) = tauXpsi(i,j)
! Calculate tau_y on psi grid
!            tauYpsi(i,j) =                                           &
!                  (  FORCES(ng) % svstr(i-1,j)                       &
!                   + FORCES(ng) % svstr(i,j)   ) * rho0 /2.
!             EBF(i,j) = tauYpsi(i,j)
! Calculate dx on psi grid
!            dx(i,j) =                                                &
!                  (  1./GRID(ng) % pm(i-1,j-1)                       &
!                   + 1./GRID(ng) % pm(i  ,j-1)                       &
!                   + 1./GRID(ng) % pm(i-1,j  )                       &
!                   + 1./GRID(ng) % pm(i  ,j  ) ) /4.
!             EBF(i,j) = dx(i,j)
! Calculate dy on psi grid
!            dy(i,j) =                                                &
!                  (  1./GRID(ng) % pn(i-1,j-1)                       &
!                   + 1./GRID(ng) % pn(i  ,j-1)                       &
!                   + 1./GRID(ng) % pn(i-1,j  )                       &
!                   + 1./GRID(ng) % pn(i  ,j  ) ) /4.
!             EBF(i,j) = dy(i,j)
! Calculate f * rho on psi grid
!            frho(i,j) =                                              &
!                  (  (OCEAN(ng) % rho(i-1,j-1,N(ng)) + rho0)         &
!                      * GRID(ng) % f(i-1,j-1)                        &
!                   + (OCEAN(ng) % rho(i  ,j-1,N(ng)) + rho0)         &
!                      * GRID(ng) % f(i  ,j-1)                        &
!                   + (OCEAN(ng) % rho(i-1,j  ,N(ng)) + rho0)         &
!                      * GRID(ng) % f(i-1,j  )                        &
!                   + (OCEAN(ng) % rho(i  ,j  ,N(ng)) + rho0)         &
!                      * GRID(ng) % f(i  ,j  )                        &
!                   ) /4.
!             EBF(i,j) = frho(i,j)
! Calculate (d rho)/dx on psi grid
!            dRdx(i,j) =                                              &
!                  (   OCEAN(ng) % rho(i,j,N(ng)) -                   &
!                           OCEAN(ng) % rho(i-1,j,N(ng))              &
!                    + OCEAN(ng) % rho(i,j-1,N(ng)) -                 &
!                            OCEAN(ng) % rho(i-1,j-1,N(ng))  ) /2.    &
!                  / dx(i,j)
!             EBF(i,j) = dRdx(i,j)
! Calculate (d rho)/dy on psi grid
!            dRde(i,j) =                                              &
!                  (   OCEAN(ng) % rho(i,j,N(ng)) -                   &
!                           OCEAN(ng) % rho(i,j-1,N(ng))              &
!                    + OCEAN(ng) % rho(i-1,j,N(ng)) -                 &
!                            OCEAN(ng) % rho(i-1,j-1,N(ng))  ) /2.    &
!                  / dy(i,j)
!             EBF(i,j) = dRde(i,j)
! Now do Ekman buoyancy flux
!            EBF(i,j) = - 9.81 / rho0 *                               &
!                  (   tauYpsi(i,j) * dRdx(i,j)                     &
!                    - tauXpsi(i,j) * dRde(i,j) )                   &
!                  / frho(i,j)
!         END DO
!       END DO
!       scale=1.0_r8
!       gtype=gfactor*p2dvar
!       status=nf_fwrite2d(ng, iNLM, QCK(ng)%ncid,                      &
!    &                     QCK(ng)%Vid(idEBF),                          &
!    &                     QCK(ng)%Rindex, gtype,                       &
!    &                     LBi, UBi, LBj, UBj, scale,                   &
!    &                     GRID(ng) % pmask_full,                       &
!    &                     EBF)
!       IF (FoundError(status, nf90_noerr, 560,                    &
!    &                 "ROMS/Utility/wrt_quick.F")) THEN
!         IF (Master) THEN
!           WRITE (stdout,10) TRIM(Vname(1,idEBF)), QCK(ng)%Rindex
!         END IF
!         exit_flag=3
!         ioerror=status
!         RETURN
!       END IF
!     END IF
!-----------------------------------------------------------------------
!  Compute surface strain.
!-----------------------------------------------------------------------
!
!     IF (Qout(idStr2,ng)) THEN !psl20210702: absent from varinfo.
!       IF (.not.allocated(strain)) THEN
!         allocate (dUde(LBi:UBi,LBj:UBj))
!           dUde(LBi:UBi,LBj:UBj)=0.0_r8
!         allocate (dVdx(LBi:UBi,LBj:UBj))
!           dVdx(LBi:UBi,LBj:UBj)=0.0_r8
!         allocate (dUdx(LBi:UBi,LBj:UBj))
!           dUdx(LBi:UBi,LBj:UBj)=0.0_r8
!         allocate (dVde(LBi:UBi,LBj:UBj))
!           dVde(LBi:UBi,LBj:UBj)=0.0_r8
!         allocate (strain(LBi:UBi,LBj:UBj))
!           strain(LBi:UBi,LBj:UBj)=0.0_r8
!       END IF
! Compute surface strain
!     DO i=LBi+1,UBi-1
!       DO j=LBj+1,UBj-1
!         cff=GRID(ng)%pm(i,j) * GRID(ng)%pn(i,j)
!         dUde(i,j)=                                                    &
!    &      GRID(ng)%om_u(i,j  ) * OCEAN(ng)%u(i,j  ,N(ng),nrhs(ng)) -      &
!    &      GRID(ng)%om_u(i,j-1) * OCEAN(ng)%u(i,j-1,N(ng),nrhs(ng))
!         dVdx(i,j)=                                                    &
!    &      GRID(ng)%on_v(i  ,j) * OCEAN(ng)%v(i  ,j,N(ng),nrhs(ng)) -      &
!    &      GRID(ng)%on_v(i-1,j) * OCEAN(ng)%v(i-1,j,N(ng),nrhs(ng))
!         dVde(i,j)=                                                    &
!    &      GRID(ng)%om_v(i,j  ) * OCEAN(ng)%v(i,j  ,N(ng),nrhs(ng)) -      &
!    &      GRID(ng)%om_v(i,j-1) * OCEAN(ng)%v(i,j-1,N(ng),nrhs(ng))
!         dUdx(i,j)=                                                    &
!    &      GRID(ng)%on_u(i  ,j) * OCEAN(ng)%u(i  ,j,N(ng),nrhs(ng)) -      &
!    &      GRID(ng)%on_u(i-1,j) * OCEAN(ng)%u(i-1,j,N(ng),nrhs(ng))
!         strain(i,j)=sqrt(  (  dUdx(i,j) - dVde(i,j)  )**2 +           &
!    &                       (  dVdx(i,j) + dUde(i,j)  )**2  ) * cff
!       END DO
!     END DO
!     scale=1.0_r8
!     gtype=gfactor*p2dvar
!     status=nf_fwrite2d(ng, iNLM, QCK(ng)%ncid,                        &
!    &                   QCK(ng)%Vid(idStr2),                           &
!    &                   QCK(ng)%Rindex, gtype,                         &
!    &                   LBi, UBi, LBj, UBj, scale,                     &
!    &                     GRID(ng) % pmask_full,                       &
!    &                     strain)
!     IF (FoundError(status, nf90_noerr, 628,                      &
!    &               "ROMS/Utility/wrt_quick.F")) THEN
!         IF (Master) THEN
!           WRITE (stdout,10) TRIM(Vname(1,idStr2)), QCK(ng)%Rindex
!         END IF
!         exit_flag=3
!         ioerror=status
!         RETURN
!       END IF
!     END IF
!-----------------------------------------------------------------------
!  Compute surface relative vorticity.
!-----------------------------------------------------------------------
!
!     IF (Qout(idRVsu,ng)) THEN !psl20210702: absent from varinfo.
!       IF (.not.allocated(rvor)) THEN
!         allocate (rvor(LBi:UBi,LBj:UBj))
!           rvor(LBi:UBi,LBj:UBj)=0.0_r8
!         allocate (rdUde(LBi:UBi,LBj:UBj))
!           rdUde(LBi:UBi,LBj:UBj)=0.0_r8
!         allocate (rdVdx(LBi:UBi,LBj:UBj))
!           rdVdx(LBi:UBi,LBj:UBj)=0.0_r8
!       END IF
! Compute surface strain
!       DO i=LBi+1,UBi-1
!         DO j=LBj+1,UBj-1
!           rdUde(i,j)=                                                 &
!    &        GRID(ng)%om_u(i,j  ) * OCEAN(ng)%u(i,j  ,N(ng),nrhs(ng)) -    &
!    &        GRID(ng)%om_u(i,j-1) * OCEAN(ng)%u(i,j-1,N(ng),nrhs(ng))
!           rdVdx(i,j)=                                                 &
!    &        GRID(ng)%on_v(i  ,j) * OCEAN(ng)%v(i  ,j,N(ng),nrhs(ng)) -    &
!    &        GRID(ng)%on_v(i-1,j) * OCEAN(ng)%v(i-1,j,N(ng),nrhs(ng))
!           cff=GRID(ng)%pm(i,j) * GRID(ng)%pn(i,j)
!           rvor(i,j)=cff*(rdVdx(i,j)-rdUde(i,j))
!         END DO
!       END DO
!       scale=1.0_r8
!       gtype=gfactor*p2dvar
!       status=nf_fwrite2d(ng, iNLM, QCK(ng)%ncid,                      &
!    &                     QCK(ng)%Vid(idRVsu),                         &
!    &                     QCK(ng)%Rindex, gtype,                       &
!    &                     LBi, UBi, LBj, UBj, scale,                   &
!    &                     GRID(ng) % pmask_full,                       &
!    &                     rvor)
!       IF (FoundError(status, nf90_noerr, 682,                    &
!    &                 "ROMS/Utility/wrt_quick.F")) THEN
!         IF (Master) THEN
!           WRITE (stdout,10) TRIM(Vname(1,idRVsu)), QCK(ng)%Rindex
!         END IF
!         exit_flag=3
!         ioerror=status
!         RETURN
!       END IF
!     END IF
!-----------------------------------------------------------------------
!  Compute surface potential vorticity.
!-----------------------------------------------------------------------
!
!     IF (Qout(idPVsu,ng).or.Qout(idPVIsu,ng).or.Qout(idPVJsu,          &!psl20210702: absent from varinfo
!    &          ng).or.Qout(idPVKsu,ng) ) THEN
!       IF (.not.allocated(pvor)) THEN
!         allocate (pvor(LBi:UBi,LBj:UBj))
!                   pvor(LBi:UBi,LBj:UBj)=0.0_r8
!       END IF
!       IF (.not.allocated(pvori)) THEN
!         allocate (pvori(LBi:UBi,LBj:UBj))
!                   pvori(LBi:UBi,LBj:UBj)=0.0_r8
!       END IF
!       IF (.not.allocated(pvorj)) THEN
!         allocate (pvorj(LBi:UBi,LBj:UBj))
!                   pvorj(LBi:UBi,LBj:UBj)=0.0_r8
!       END IF
!       IF (.not.allocated(pvork)) THEN
!         allocate (pvork(LBi:UBi,LBj:UBj))
!                   pvork(LBi:UBi,LBj:UBj)=0.0_r8
!       END IF
!       IF (.not.allocated(dz_r)) THEN
!         allocate (dz_r(LBi:UBi,LBj:UBj))
!                   dz_r(LBi:UBi,LBj:UBj)=0.0_r8
!       END IF
!       IF (.not.allocated(dz_u)) THEN
!         allocate (dz_u(LBi:UBi,LBj:UBj))
!                   dz_u(LBi:UBi,LBj:UBj)=0.0_r8
!       END IF
!       IF (.not.allocated(dz_v)) THEN
!         allocate (dz_v(LBi:UBi,LBj:UBj))
!                   dz_v(LBi:UBi,LBj:UBj)=0.0_r8
!       END IF
!       IF (.not.allocated(dRdx)) THEN
!         allocate (dRdx(LBi:UBi,LBj:UBj))
!           dRdx(LBi:UBi,LBj:UBj)=0.0_r8
!       END IF
!       IF (.not.allocated(dRde)) THEN
!         allocate (dRde(LBi:UBi,LBj:UBj))
!           dRde(LBi:UBi,LBj:UBj)=0.0_r8
!       END IF
!       IF (.not.allocated(dRdz)) THEN
!         allocate (dRdz(LBi:UBi,LBj:UBj))
!           dRdz(LBi:UBi,LBj:UBj)=0.0_r8
!       END IF
!       IF (.not.allocated(dUdz)) THEN
!         allocate (dUdz(LBi:UBi,LBj:UBj))
!           dUdz(LBi:UBi,LBj:UBj)=0.0_r8
!       END IF
!       IF (.not.allocated(dVdz)) THEN
!         allocate (dVdz(LBi:UBi,LBj:UBj))
!           dVdz(LBi:UBi,LBj:UBj)=0.0_r8
!       END IF
!       IF (.not.allocated(dUde)) THEN
!         allocate (dUde(LBi:UBi,LBj:UBj))
!                   dUde(LBi:UBi,LBj:UBj)=0.0_r8
!       END IF
!       IF (.not.allocated(dVdx)) THEN
!         allocate (dVdx(LBi:UBi,LBj:UBj))
!                   dVdx(LBi:UBi,LBj:UBj)=0.0_r8
!       END IF
!
! Find dz on the various grids
!
!       DO j=Jstr,JendR
!         DO i=Istr,IendR
!           dz_r(i,j)= GRID(ng)%z_w(i,j,N(ng))-GRID(ng)%z_w(i,j,N(ng)-1)
!         END DO
!       END DO
!       DO j=Jstr,JendR
!         DO i=Istr,IendR
!           dz_u(i,j) = ( dz_r(i,j) + dz_r(i+1,j  ) ) / 2
!           dz_v(i,j) = ( dz_r(i,j) + dz_r(i  ,j+1) ) / 2
!         END DO
!       END DO
!       DO j=Jstr,JendR
!         dz_r(Istr-1,j) = dz_r(Istr,j)
!         dz_u(Istr-1,j) = dz_u(Istr,j)
!         dz_v(Istr-1,j) = dz_v(Istr,j)
!       END DO
!       DO i=Istr,IendR
!         dz_r(i,Jstr-1) = dz_r(i,Jstr)
!         dz_u(i,Jstr-1) = dz_u(i,Jstr)
!         dz_v(i,Jstr-1) = dz_v(i,Jstr)
!       END DO
!       dz_r(Istr-1,Jstr-1) = dz_r(Istr,Jstr)
!       dz_u(Istr-1,Jstr-1) = dz_u(Istr,Jstr)
!       dz_v(Istr-1,Jstr-1) = dz_v(Istr,Jstr)
!
! find the horizontal derivatives on the psi grid
!
!       DO j=Jstr,JendR
!         DO i=Istr,IendR
!           cff = ( OCEAN(ng) % rho( i  ,j  , N(ng) ) -                 &
!    &              OCEAN(ng) % rho( i-1,j  , N(ng) )    )              &
!    &             / GRID(ng) % om_r(i  ,j  )  +                        &
!    &            ( OCEAN(ng) % rho( i  ,j-1, N(ng) ) -                 &
!    &              OCEAN(ng) % rho( i-1,j-1, N(ng) )    )              &
!    &             / GRID(ng) % om_r(i  ,j-1)
!           dRdx(i,j) = cff / 2
!           cff = ( OCEAN(ng) % rho( i  ,j  , N(ng) ) -                 &
!    &              OCEAN(ng) % rho( i,  j-1, N(ng) )    )              &
!    &             / GRID(ng) % on_r(i  ,j  )  +                        &
!    &            ( OCEAN(ng) % rho( i-1,j  , N(ng) ) -                 &
!    &              OCEAN(ng) % rho( i-1,j-1, N(ng) )    )              &
!    &             / GRID(ng) % om_r(i-1,j  )
!           dRde(i,j) = cff / 2
!           dVdx(i,j)= ( OCEAN(ng) % v(i  ,j, N(ng) , nrhs(ng)) -           &
!    &                   OCEAN(ng) % v(i-1,j, N(ng) , nrhs(ng))    )        &
!    &             / GRID(ng) % om_v(i  ,j)
!           dUde(i,j)= ( OCEAN(ng) % u(i  ,j, N(ng) , nrhs(ng)) -           &
!    &                   OCEAN(ng) % u(i,j-1, N(ng) , nrhs(ng))    )        &
!    &             / GRID(ng) % on_u(i  ,j)
!           dUde(i,j)=dUde(i,j)*GRID(ng)%pmask(i,j)
!           dVdx(i,j)=dVdx(i,j)*GRID(ng)%pmask(i,j)
!           dRde(i,j)=dRde(i,j)*GRID(ng)%pmask(i,j)
!           drdx(i,j)=drdx(i,j)*GRID(ng)%pmask(i,j)
!         END DO
!       END DO
!
! vertical derivatives
!
!       DO j=Jstr,JendR
!         DO i=Istr,IendR
!           cff = ( OCEAN(ng) % rho(i  ,j  ,N(ng)  ) -                  &
!    &              OCEAN(ng) % rho(i  ,j  ,N(ng)-1)   )  /             &
!    &            (  GRID(ng) % z_r(i  ,j  ,N(ng)  ) -                  &
!    &               GRID(ng) % z_r(i  ,j  ,N(ng)-1)   )
!           cff = cff +                                                 &
!    &            ( OCEAN(ng) % rho(i-1,j  ,N(ng)  ) -                  &
!    &              OCEAN(ng) % rho(i-1,j  ,N(ng)-1)   )   /            &
!    &            (  GRID(ng) % z_r(i-1,j  ,N(ng)  ) -                  &
!    &               GRID(ng) % z_r(i-1,j  ,N(ng)-1)   )
!           cff = cff +                                                 &
!    &            ( OCEAN(ng) % rho(i  ,j-1,N(ng)  ) -                  &
!    &              OCEAN(ng) % rho(i  ,j-1,N(ng)-1)   )   /            &
!    &            (  GRID(ng) % z_r(i  ,j-1,N(ng)  ) -                  &
!    &               GRID(ng) % z_r(i  ,j-1,N(ng)-1)   )
!           cff = cff +                                                 &
!    &            ( OCEAN(ng) % rho(i-1,j-1,N(ng)  ) -                  &
!    &              OCEAN(ng) % rho(i-1,j-1,N(ng)-1)   )   /            &
!    &            (  GRID(ng) % z_r(i-1,j-1,N(ng)  ) -                  &
!    &               GRID(ng) % z_r(i-1,j-1,N(ng)-1)   )
!           dRdz(i,j) = cff / 4
!           cff = ( OCEAN(ng) % v(i  ,j  ,N(ng)  ,nrhs(ng)) -               &
!    &              OCEAN(ng) % v(i  ,j  ,N(ng)-1,nrhs(ng))   ) / dz_v(i,j)
!           cff = cff +                                                 &
!    &         ( OCEAN(ng) % v(i-1,j  ,N(ng)  ,nrhs(ng)) -                  &
!    &           OCEAN(ng) % v(i-1,j  ,N(ng)-1,nrhs(ng))   ) / dz_v(i-1,j)
!           dVdz(i,j) = cff / 2
!           cff = ( OCEAN(ng) % u(i  ,j  ,N(ng)  ,nrhs(ng)) -               &
!    &              OCEAN(ng) % u(i  ,j  ,N(ng)-1,nrhs(ng))   ) / dz_u(i,j)
!           cff = cff +                                                 &
!    &         ( OCEAN(ng) % u(i  ,j-1,N(ng)  ,nrhs(ng)) -                  &
!    &           OCEAN(ng) % u(i  ,j-1,N(ng)-1,nrhs(ng))   ) / dz_u(i,j-1)
!           dUdz(i,j) = cff / 2
!         END DO
!       END DO
!
!  Compute potential vorticity (meter-1 second-1)
!  at horizontal PSI-points and vertical RHO-points.
!
!       DO j=Jstr,JendR
!         DO i=Istr,IendR
!           fomn_p=0.25_r8*(GRID(ng)%f(i-1,j-1)+GRID(ng)%f(i-1,j)+      &
!    &                GRID(ng)%f(i,j-1)+GRID(ng)%f(i,j))
!           pvori(i,j) = -dVdz(i,j) * dRdx(i,j) / rho0
!           pvorj(i,j) =  dUdz(i,j) * dRde(i,j) / rho0
!           pvork(i,j) =  dRdz(i,j) *                                   &
!    &                     (fomn_p + dVdx(i,j) - dUde(i,j) ) / rho0
!           pvor(i,j) = pvori(i,j) + pvorj(i,j) + pvork(i,j)
!         END DO
!       END DO
!
!       IF (Qout(idPVsu,ng)) THEN
!         scale=1.0_r8
!         gtype=gfactor*p2dvar
!         status=nf_fwrite2d(ng, iNLM, QCK(ng)%ncid,                    &
!    &                     QCK(ng)%Vid(idPVsu),                         &
!    &                     QCK(ng)%Rindex, gtype,                       &
!    &                     LBi, UBi, LBj, UBj, scale,                   &
!    &                     GRID(ng) % pmask_full,                       &
!    &                     pvor)
!         IF (FoundError(status, nf90_noerr, 921,                  &
!    &                 "ROMS/Utility/wrt_quick.F")) THEN
!           IF (Master) THEN
!             WRITE (stdout,10) TRIM(Vname(1,idPVsu)), QCK(ng)%Rindex
!           END IF
!           exit_flag=3
!           ioerror=status
!           RETURN
!         END IF
!       END IF
!
!       IF (Qout(idPVIsu,ng)) THEN
!         scale=1.0_r8
!         gtype=gfactor*p2dvar
!         status=nf_fwrite2d(ng, iNLM, QCK(ng)%ncid,                    &
!    &                     QCK(ng)%Vid(idPVIsu),                        &
!    &                     QCK(ng)%Rindex, gtype,                       &
!    &                     LBi, UBi, LBj, UBj, scale,                   &
!    &                     GRID(ng) % pmask_full,                       &
!    &                     pvori)
!         IF (FoundError(status, nf90_noerr, 943,                  &
!    &                 "ROMS/Utility/wrt_quick.F")) THEN
!           IF (Master) THEN
!             WRITE (stdout,10) TRIM(Vname(1,idPVIsu)), QCK(ng)%Rindex
!           END IF
!           exit_flag=3
!           ioerror=status
!           RETURN
!         END IF
!       END IF
!
!       IF (Qout(idPVJsu,ng)) THEN
!         scale=1.0_r8
!         gtype=gfactor*p2dvar
!         status=nf_fwrite2d(ng, iNLM, QCK(ng)%ncid,                    &
!    &                     QCK(ng)%Vid(idPVJsu),                        &
!    &                     QCK(ng)%Rindex, gtype,                       &
!    &                     LBi, UBi, LBj, UBj, scale,                   &
!    &                     GRID(ng) % pmask_full,                       &
!    &                     pvorj)
!         IF (FoundError(status, nf90_noerr, 965,                  &
!    &                 "ROMS/Utility/wrt_quick.F")) THEN
!           IF (Master) THEN
!             WRITE (stdout,10) TRIM(Vname(1,idPVJsu)), QCK(ng)%Rindex
!           END IF
!           exit_flag=3
!           ioerror=status
!           RETURN
!         END IF
!       END IF
!
!       IF (Qout(idPVKsu,ng)) THEN
!         scale=1.0_r8
!         gtype=gfactor*p2dvar
!         status=nf_fwrite2d(ng, iNLM, QCK(ng)%ncid,                    &
!    &                     QCK(ng)%Vid(idPVKsu),                        &
!    &                     QCK(ng)%Rindex, gtype,                       &
!    &                     LBi, UBi, LBj, UBj, scale,                   &
!    &                     GRID(ng) % pmask_full,                       &
!    &                     pvork)
!         IF (FoundError(status, nf90_noerr, 987,                  &
!    &                 "ROMS/Utility/wrt_quick.F")) THEN
!           IF (Master) THEN
!             WRITE (stdout,10) TRIM(Vname(1,idPVKsu)), QCK(ng)%Rindex
!           END IF
!           exit_flag=3
!           ioerror=status
!           RETURN
!         END IF
!       END IF
!     END IF
!
!  Write out 2D Eastward and Northward momentum components (m/s) at
!  RHO-points.
!
      IF (Qout(idu2dE,ng).and.Qout(idv2dN,ng)) THEN
        IF (.not.allocated(Ur2d)) THEN
          allocate (Ur2d(LBi:UBi,LBj:UBj))
            Ur2d(LBi:UBi,LBj:UBj)=0.0_r8
        END IF
        IF (.not.allocated(Vr2d)) THEN
          allocate (Vr2d(LBi:UBi,LBj:UBj))
            Vr2d(LBi:UBi,LBj:UBj)=0.0_r8
        END IF
        CALL uv_rotate2d (ng, tile, .FALSE., .TRUE.,                    &
     &                    LBi, UBi, LBj, UBj,                           &
     &                    GRID(ng) % CosAngler,                         &
     &                    GRID(ng) % SinAngler,                         &
     &                    GRID(ng) % rmask_full,                        &
     &                    OCEAN(ng) % ubar(:,:,kstp(ng)),                   &
     &                    OCEAN(ng) % vbar(:,:,kstp(ng)),                   &
     &                    Ur2d, Vr2d)
        scale=1.0_r8
        gtype=gfactor*r2dvar
        status=nf_fwrite2d(ng, iNLM, QCK(ng)%ncid, QCK(ng)%Vid(idu2dE), &
     &                     QCK(ng)%Rindex, gtype,                       &
     &                     LBi, UBi, LBj, UBj, scale,                   &
     &                     GRID(ng) % rmask_full,                       &
     &                     Ur2d)
        IF (FoundError(status, nf90_noerr, 1031,                    &
     &                 "ROMS/Utility/wrt_quick.F")) THEN
          IF (Master) THEN
            WRITE (stdout,10) TRIM(Vname(1,idu2dE)), QCK(ng)%Rindex
          END IF
          exit_flag=3
          ioerror=status
          RETURN
        END IF
        status=nf_fwrite2d(ng, iNLM, QCK(ng)%ncid, QCK(ng)%Vid(idv2dN), &
     &                     QCK(ng)%Rindex, gtype,                       &
     &                     LBi, UBi, LBj, UBj, scale,                   &
     &                     GRID(ng) % rmask_full,                       &
     &                     Vr2d)
        IF (FoundError(status, nf90_noerr, 1048,                    &
     &                 "ROMS/Utility/wrt_quick.F")) THEN
          IF (Master) THEN
            WRITE (stdout,10) TRIM(Vname(1,idv2dN)), QCK(ng)%Rindex
          END IF
          exit_flag=3
          ioerror=status
          RETURN
        END IF
        deallocate (Ur2d)
        deallocate (Vr2d)
      END IF
!
!  Write out 3D U-momentum component (m/s).
!
      IF (Qout(idUvel,ng)) THEN
        scale=1.0_r8
        gtype=gfactor*u3dvar
        status=nf_fwrite3d(ng, iNLM, QCK(ng)%ncid, QCK(ng)%Vid(idUvel), &
     &                     QCK(ng)%Rindex, gtype,                       &
     &                     LBi, UBi, LBj, UBj, 1, N(ng), scale,         &
     &                     GRID(ng) % umask_full,                       &
     &                     OCEAN(ng) % u(:,:,:,nrhs(ng)))
        IF (FoundError(status, nf90_noerr, 1075,                    &
     &                 "ROMS/Utility/wrt_quick.F")) THEN
          IF (Master) THEN
            WRITE (stdout,10) TRIM(Vname(1,idUvel)), QCK(ng)%Rindex
          END IF
          exit_flag=3
          ioerror=status
          RETURN
        END IF
      END IF
!
!  Write out 3D V-momentum component (m/s).
!
      IF (Qout(idVvel,ng)) THEN
        scale=1.0_r8
        gtype=gfactor*v3dvar
        status=nf_fwrite3d(ng, iNLM, QCK(ng)%ncid, QCK(ng)%Vid(idVvel), &
     &                     QCK(ng)%Rindex, gtype,                       &
     &                     LBi, UBi, LBj, UBj, 1, N(ng), scale,         &
     &                     GRID(ng) % vmask_full,                       &
     &                     OCEAN(ng) % v(:,:,:,nrhs(ng)))
        IF (FoundError(status, nf90_noerr, 1098,                    &
     &                 "ROMS/Utility/wrt_quick.F")) THEN
          IF (Master) THEN
            WRITE (stdout,10) TRIM(Vname(1,idVvel)), QCK(ng)%Rindex
          END IF
          exit_flag=3
          ioerror=status
          RETURN
        END IF
      END IF
!
!  Write out surface U-momentum component (m/s).
!
      IF (Qout(idUsur,ng)) THEN
        scale=1.0_r8
        gtype=gfactor*u2dvar
        status=nf_fwrite2d(ng, iNLM, QCK(ng)%ncid, QCK(ng)%Vid(idUsur), &
     &                     QCK(ng)%Rindex, gtype,                       &
     &                     LBi, UBi, LBj, UBj, scale,                   &
     &                     GRID(ng) % umask_full,                       &
     &                     OCEAN(ng) % u(:,:,N(ng),nrhs(ng)))
        IF (FoundError(status, nf90_noerr, 1121,                    &
     &                 "ROMS/Utility/wrt_quick.F")) THEN
          IF (Master) THEN
            WRITE (stdout,10) TRIM(Vname(1,idUsur)), QCK(ng)%Rindex
          END IF
          exit_flag=3
          ioerror=status
          RETURN
        END IF
      END IF
!
!  Write out surface V-momentum component (m/s).
!
      IF (Qout(idVsur,ng)) THEN
        scale=1.0_r8
        gtype=gfactor*v2dvar
        status=nf_fwrite2d(ng, iNLM, QCK(ng)%ncid, QCK(ng)%Vid(idVsur), &
     &                     QCK(ng)%Rindex, gtype,                       &
     &                     LBi, UBi, LBj, UBj, scale,                   &
     &                     GRID(ng) % vmask_full,                       &
     &                     OCEAN(ng) % v(:,:,N(ng),nrhs(ng)))
        IF (FoundError(status, nf90_noerr, 1144,                    &
     &                 "ROMS/Utility/wrt_quick.F")) THEN
          IF (Master) THEN
            WRITE (stdout,10) TRIM(Vname(1,idVsur)), QCK(ng)%Rindex
          END IF
          exit_flag=3
          ioerror=status
          RETURN
        END IF
      END IF
!
!  Write out 3D Eastward and Northward momentum components (m/s) at
!  RHO-points.
!
      IF ((Qout(idu3dE,ng).and.Qout(idv3dN,ng)).or.                     &
     &    (Qout(idUsuE,ng).and.Qout(idVsuN,ng))) THEN
        IF (.not.allocated(Ur3d)) THEN
          allocate (Ur3d(LBi:UBi,LBj:UBj,N(ng)))
          Ur3d(LBi:UBi,LBj:UBj,1:N(ng))=0.0_r8
        END IF
        IF (.not.allocated(Vr3d)) THEN
          allocate (Vr3d(LBi:UBi,LBj:UBj,N(ng)))
          Vr3d(LBi:UBi,LBj:UBj,1:N(ng))=0.0_r8
        END IF
        CALL uv_rotate3d (ng, tile, .FALSE., .TRUE.,                    &
     &                    LBi, UBi, LBj, UBj, 1, N(ng),                 &
     &                    GRID(ng) % CosAngler,                         &
     &                    GRID(ng) % SinAngler,                         &
     &                    GRID(ng) % rmask_full,                        &
     &                    OCEAN(ng) % u(:,:,:,nrhs(ng)),                    &
     &                    OCEAN(ng) % v(:,:,:,nrhs(ng)),                    &
     &                    Ur3d, Vr3d)
        IF ((Qout(idu3dE,ng).and.Qout(idv3dN,ng))) THEN
          scale=1.0_r8
          gtype=gfactor*r3dvar
          status=nf_fwrite3d(ng, iNLM, QCK(ng)%ncid,                    &
     &                       QCK(ng)%Vid(idu3dE),                       &
     &                       QCK(ng)%Rindex, gtype,                     &
     &                       LBi, UBi, LBj, UBj, 1, N(ng), scale,       &
     &                       GRID(ng) % rmask_full,                     &
     &                       Ur3d)
          IF (FoundError(status, nf90_noerr, 1190,                  &
     &                   "ROMS/Utility/wrt_quick.F")) THEN
            IF (Master) THEN
              WRITE (stdout,10) TRIM(Vname(1,idu3dE)), QCK(ng)%Rindex
            END IF
            exit_flag=3
            ioerror=status
            RETURN
          END IF
          status=nf_fwrite3d(ng, iNLM, QCK(ng)%ncid,                    &
     &                       QCK(ng)%Vid(idv3dN),                       &
     &                       QCK(ng)%Rindex, gtype,                     &
     &                       LBi, UBi, LBj, UBj, 1, N(ng), scale,       &
     &                       GRID(ng) % rmask_full,                     &
     &                       Vr3d)
          IF (FoundError(status, nf90_noerr, 1208,                  &
     &                   "ROMS/Utility/wrt_quick.F")) THEN
            IF (Master) THEN
              WRITE (stdout,10) TRIM(Vname(1,idv3dN)), QCK(ng)%Rindex
            END IF
            exit_flag=3
            ioerror=status
            RETURN
          END IF
        END IF
!
!  Write out surface Eastward and Northward momentum components (m/s) at
!  RHO-points.
!
        IF ((Qout(idUsuE,ng).and.Qout(idVsuN,ng))) THEN
          scale=1.0_r8
          gtype=gfactor*r2dvar
          status=nf_fwrite2d(ng, iNLM, QCK(ng)%ncid,                    &
     &                       QCK(ng)%Vid(idUsuE),                       &
     &                       QCK(ng)%Rindex, gtype,                     &
     &                       LBi, UBi, LBj, UBj, scale,                 &
     &                       GRID(ng) % rmask_full,                     &
     &                       Ur3d(:,:,N(ng)))
          IF (FoundError(status, nf90_noerr, 1233,                  &
     &                   "ROMS/Utility/wrt_quick.F")) THEN
            IF (Master) THEN
              WRITE (stdout,10) TRIM(Vname(1,idUsuE)), QCK(ng)%Rindex
            END IF
            exit_flag=3
            ioerror=status
            RETURN
          END IF
          status=nf_fwrite2d(ng, iNLM, QCK(ng)%ncid,                    &
     &                       QCK(ng)%Vid(idVsuN),                       &
     &                       QCK(ng)%Rindex, gtype,                     &
     &                       LBi, UBi, LBj, UBj, scale,                 &
     &                       GRID(ng) % rmask_full,                     &
     &                       Vr3d(:,:,N(ng)))
          IF (FoundError(status, nf90_noerr, 1251,                  &
     &                   "ROMS/Utility/wrt_quick.F")) THEN
            IF (Master) THEN
              WRITE (stdout,10) TRIM(Vname(1,idVsuN)), QCK(ng)%Rindex
            END IF
            exit_flag=3
            ioerror=status
            RETURN
          END IF
        END IF
        deallocate (Ur3d)
        deallocate (Vr3d)
      END IF
!
!  Write out S-coordinate omega vertical velocity (m/s).
!
      IF (Qout(idOvel,ng)) THEN
        IF (.not.allocated(Wr3d)) THEN
          allocate (Wr3d(LBi:UBi,LBj:UBj,0:N(ng)))
          Wr3d(LBi:UBi,LBj:UBj,0:N(ng))=0.0_r8
        END IF
        scale=1.0_r8
        gtype=gfactor*w3dvar
        CALL scale_omega (ng, tile, LBi, UBi, LBj, UBj, 0, N(ng),       &
     &                    GRID(ng) % pm,                                &
     &                    GRID(ng) % pn,                                &
     &                    OCEAN(ng) % W,                                &
     &                    Wr3d)
        status=nf_fwrite3d(ng, iNLM, QCK(ng)%ncid, QCK(ng)%Vid(idOvel), &
     &                     QCK(ng)%Rindex, gtype,                       &
     &                     LBi, UBi, LBj, UBj, 0, N(ng), scale,         &
     &                     GRID(ng) % rmask,                            &
     &                     Wr3d)
        IF (FoundError(status, nf90_noerr, 1286,                    &
     &                 "ROMS/Utility/wrt_quick.F")) THEN
          IF (Master) THEN
            WRITE (stdout,10) TRIM(Vname(1,idOvel)), QCK(ng)%Rindex
          END IF
          exit_flag=3
          ioerror=status
          RETURN
        END IF
        deallocate (Wr3d)
      END IF
!
!  Write out vertical velocity (m/s).
!
      IF (Qout(idWvel,ng)) THEN
        scale=1.0_r8
        gtype=gfactor*w3dvar
        status=nf_fwrite3d(ng, iNLM, QCK(ng)%ncid, QCK(ng)%Vid(idWvel), &
     &                     QCK(ng)%Rindex, gtype,                       &
     &                     LBi, UBi, LBj, UBj, 0, N(ng), scale,         &
     &                     GRID(ng) % rmask,                            &
     &                     OCEAN(ng) % wvel)
        IF (FoundError(status, nf90_noerr, 1310,                    &
     &                 "ROMS/Utility/wrt_quick.F")) THEN
          IF (Master) THEN
            WRITE (stdout,10) TRIM(Vname(1,idWvel)), QCK(ng)%Rindex
          END IF
          exit_flag=3
          ioerror=status
          RETURN
        END IF
      END IF
!
!  Write out tracer type variables.
!
      DO itrc=1,NT(ng)
        IF (Qout(idTvar(itrc),ng)) THEN
          scale=1.0_r8
          gtype=gfactor*r3dvar
          status=nf_fwrite3d(ng, iNLM, QCK(ng)%ncid, QCK(ng)%Tid(itrc), &
     &                       QCK(ng)%Rindex, gtype,                     &
     &                       LBi, UBi, LBj, UBj, 1, N(ng), scale,       &
     &                       GRID(ng) % rmask,                          &
     &                       OCEAN(ng) % t(:,:,:,nrhs(ng),itrc))
          IF (FoundError(status, nf90_noerr, 1334,                  &
     &                   "ROMS/Utility/wrt_quick.F")) THEN
            IF (Master) THEN
              WRITE (stdout,10) TRIM(Vname(1,idTvar(itrc))),            &
     &                          QCK(ng)%Rindex
            END IF
            exit_flag=3
            ioerror=status
            RETURN
          END IF
        END IF
      END DO
!
!  Write out surface tracer type variables.
!
      DO itrc=1,NT(ng)
        IF (Qout(idsurT(itrc),ng)) THEN
          scale=1.0_r8
          gtype=gfactor*r2dvar
          status=nf_fwrite2d(ng, iNLM, QCK(ng)%ncid,                    &
     &                       QCK(ng)%Vid(idsurT(itrc)),                 &
     &                       QCK(ng)%Rindex, gtype,                     &
     &                       LBi, UBi, LBj, UBj, scale,                 &
     &                       GRID(ng) % rmask,                          &
     &                       OCEAN(ng) % t(:,:,N(ng),nrhs(ng),itrc))
          IF (FoundError(status, nf90_noerr, 1361,                  &
     &                   "ROMS/Utility/wrt_quick.F")) THEN
            IF (Master) THEN
              WRITE (stdout,10) TRIM(Vname(1,idsurT(itrc))),            &
     &                          QCK(ng)%Rindex
            END IF
            exit_flag=3
            ioerror=status
            RETURN
          END IF
        END IF
      END DO
!
!  Write out density anomaly.
!
      IF (Qout(idDano,ng)) THEN
        scale=1.0_r8
        gtype=gfactor*r3dvar
        status=nf_fwrite3d(ng, iNLM, QCK(ng)%ncid, QCK(ng)%Vid(idDano), &
     &                     QCK(ng)%Rindex, gtype,                       &
     &                     LBi, UBi, LBj, UBj, 1, N(ng), scale,         &
     &                     GRID(ng) % rmask,                            &
     &                     OCEAN(ng) % rho)
        IF (FoundError(status, nf90_noerr, 1386,                    &
     &                 "ROMS/Utility/wrt_quick.F")) THEN
          IF (Master) THEN
            WRITE (stdout,10) TRIM(Vname(1,idDano)), QCK(ng)%Rindex
          END IF
          exit_flag=3
          ioerror=status
          RETURN
        END IF
      END IF
!
!  Write out depth surface boundary layer.
!
      IF (Qout(idHsbl,ng)) THEN
        scale=1.0_r8
        gtype=gfactor*r2dvar
        status=nf_fwrite2d(ng, iNLM, QCK(ng)%ncid, QCK(ng)%Vid(idHsbl), &
     &                     QCK(ng)%Rindex, gtype,                       &
     &                     LBi, UBi, LBj, UBj, scale,                   &
     &                     GRID(ng) % rmask,                            &
     &                     MIXING(ng) % hsbl)
        IF (FoundError(status, nf90_noerr, 1410,                    &
     &                 "ROMS/Utility/wrt_quick.F")) THEN
          IF (Master) THEN
            WRITE (stdout,10) TRIM(Vname(1,idHsbl)), QCK(ng)%Rindex
          END IF
          exit_flag=3
          ioerror=status
          RETURN
        END IF
      END IF
!
!  Write out depth surface boundary layer.
!
      IF (Qout(idHbbl,ng)) THEN
        scale=1.0_r8
        gtype=gfactor*r2dvar
        status=nf_fwrite2d(ng, iNLM, QCK(ng)%ncid, QCK(ng)%Vid(idHbbl), &
     &                     QCK(ng)%Rindex, gtype,                       &
     &                     LBi, UBi, LBj, UBj, scale,                   &
     &                     GRID(ng) % rmask,                            &
     &                     MIXING(ng) % hbbl)
        IF (FoundError(status, nf90_noerr, 1435,                    &
     &                 "ROMS/Utility/wrt_quick.F")) THEN
          IF (Master) THEN
            WRITE (stdout,10) TRIM(Vname(1,idHbbl)), QCK(ng)%Rindex
          END IF
          exit_flag=3
          ioerror=status
          RETURN
        END IF
      END IF
!
!  Write out ice 2D momentum component (m/s) in the XI-direction.
!
      IF (Qout(idUice,ng)) THEN
        scale=1.0_r8
        gtype=gfactor*u2dvar
        status=nf_fwrite2d(ng, iNLM, QCK(ng)%ncid,                     &
     &                     QCK(ng)%Vid(idUice),                        &
     &                     QCK(ng)%Rindex, gtype,                      &
     &                     LBi, UBi, LBj, UBj, scale,                   &
     &                     GRID(ng) % umask_full,                       &
     &                     ICE(ng) % ui(:,:,liunw(ng)))
        IF (FoundError(status, nf90_noerr, 1461,                    &
     &                 "ROMS/Utility/wrt_quick.F")) THEN
          IF (Master) THEN
            WRITE (stdout,10) TRIM(Vname(1,idUice)), QCK(ng)%Rindex
          END IF
          exit_flag=3
          RETURN
        END IF
      END IF
!
!  Write out ice 2D momentum component (m/s) in the ETA-direction.
!
      IF (Qout(idVice,ng)) THEN
        scale=1.0_r8
        gtype=gfactor*v2dvar
        status=nf_fwrite2d(ng, iNLM, QCK(ng)%ncid,                      &
     &                     QCK(ng)%Vid(idVice),                         &
     &                     QCK(ng)%Rindex, gtype,                       &
     &                     LBi, UBi, LBj, UBj, scale,                   &
     &                     GRID(ng) % vmask_full,                       &
     &                     ICE(ng) % vi(:,:,liunw(ng)))
        IF (FoundError(status, nf90_noerr, 1484,                    &
     &                 "ROMS/Utility/wrt_quick.F")) THEN
          IF (Master) THEN
            WRITE (stdout,10) TRIM(Vname(1,idVice)), QCK(ng)%Rindex
          END IF
          exit_flag=3
          RETURN
        END IF
      END IF
!
!  Write out 2D Eastward ice momentum component (m/s) at RHO-points.
!
      IF (Qout(idUiceE,ng).and.Qout(idViceN,ng)) THEN
        IF (.not.allocated(Ur2d)) THEN
          allocate (Ur2d(LBi:UBi,LBj:UBj))
            Ur2d(LBi:UBi,LBj:UBj)=0.0_r8
        END IF
        IF (.not.allocated(Vr2d)) THEN
          allocate (Vr2d(LBi:UBi,LBj:UBj))
            Vr2d(LBi:UBi,LBj:UBj)=0.0_r8
        END IF
        CALL uv_rotate2d (ng, TILE, .FALSE., .TRUE.,                    &
     &                    LBi, UBi, LBj, UBj,                           &
     &                    GRID(ng) % CosAngler,                         &
     &                    GRID(ng) % SinAngler,                         &
     &                    GRID(ng) % rmask_full,                        &
     &                    ICE(ng) % ui(:,:,liunw(ng)),                      &
     &                    ICE(ng) % vi(:,:,liunw(ng)),                      &
     &                    Ur2d, Vr2d)
        scale=1.0_r8
        gtype=gfactor*r2dvar
        status=nf_fwrite2d(ng, iNLM, QCK(ng)%ncid,                      &
     &                     QCK(ng)%Vid(idUiceE),                        &
     &                     QCK(ng)%Rindex, gtype,                       &
     &                     LBi, UBi, LBj, UBj, scale,                   &
     &                     GRID(ng) % rmask_full,                       &
     &                     Ur2d)
        IF (FoundError(status, nf90_noerr, 1525,                    &
     &                 "ROMS/Utility/wrt_quick.F")) THEN
          IF (Master) THEN
            WRITE (stdout,10) TRIM(Vname(1,idUiceE)), QCK(ng)%Rindex
          END IF
          exit_flag=3
          ioerror=status
          RETURN
        END IF
!
!  Write out 2D Northward ice momentum component (m/s) at RHO-points.
!
        scale=1.0_r8
        gtype=gfactor*r2dvar
        status=nf_fwrite2d(ng, iNLM, QCK(ng)%ncid,                      &
     &                     QCK(ng)%Vid(idViceN),                        &
     &                     QCK(ng)%Rindex, gtype,                       &
     &                     LBi, UBi, LBj, UBj, scale,                   &
     &                     GRID(ng) % rmask_full,                       &
     &                     Vr2d)
        IF (FoundError(status, nf90_noerr, 1547,                    &
     &                 "ROMS/Utility/wrt_quick.F")) THEN
          IF (Master) THEN
            WRITE (stdout,10) TRIM(Vname(1,idViceN)), QCK(ng)%Rindex
          END IF
          exit_flag=3
          ioerror=status
          RETURN
        END IF
      END IF
!
! Write out ice concentration.
!
      IF (Qout(idAice,ng)) THEN
        scale=1.0_r8
        gtype=gfactor*r2dvar
        status=nf_fwrite2d(ng, iNLM, QCK(ng)%ncid,                      &
     &                     QCK(ng)%Vid(idAice),                         &
     &                     QCK(ng)%Rindex, gtype,                       &
     &                     LBi, UBi, LBj, UBj, scale,                   &
     &                     GRID(ng) % rmask_full,                       &
     &                     ICE(ng) % ai(:,:,liunw(ng)))
        IF (FoundError(status, nf90_noerr, 1572,                    &
     &                 "ROMS/Utility/wrt_quick.F")) THEN
          IF (Master) THEN
            WRITE (stdout,10) TRIM(Vname(1,idAice)), QCK(ng)%Rindex
          END IF
          exit_flag=3
          RETURN
        END IF
      END IF
!
! Write out ice thickness.
!
      IF (Qout(idHice,ng)) THEN
        scale=1.0_r8
        gtype=gfactor*r2dvar
        status=nf_fwrite2d(ng, iNLM, QCK(ng)%ncid,                      &
     &                     QCK(ng)%Vid(idHice),                         &
     &                     QCK(ng)%Rindex, gtype,                       &
     &                     LBi, UBi, LBj, UBj, scale,                   &
     &                     GRID(ng) % rmask_full,                       &
     &                     ICE(ng) % hi(:,:,liunw(ng)))
        IF (FoundError(status, nf90_noerr, 1595,                    &
     &                 "ROMS/Utility/wrt_quick.F")) THEN
          IF (Master) THEN
            WRITE (stdout,10) TRIM(Vname(1,idHice)), QCK(ng)%Rindex
          END IF
          exit_flag=3
          RETURN
        END IF
      END IF
!
! Write out ice/snow surface temperature.
!
      IF (Qout(idTice,ng)) THEN
        scale=1.0_r8
        gtype=gfactor*r2dvar
        status=nf_fwrite2d(ng, iNLM, QCK(ng)%ncid,                      &
     &                     QCK(ng)%Vid(idTice),                         &
     &                     QCK(ng)%Rindex, gtype,                       &
     &                     LBi, UBi, LBj, UBj, scale,                   &
     &                     GRID(ng) % rmask_full,                       &
     &                     ICE(ng) % tis)
        IF (FoundError(status, nf90_noerr, 1618,                    &
     &                 "ROMS/Utility/wrt_quick.F")) THEN
          IF (Master) THEN
            WRITE (stdout,10) TRIM(Vname(1,idTice)), QCK(ng)%Rindex
          END IF
          exit_flag=3
          RETURN
        END IF
      END IF
!
! Write out ice internal temperature.
!
      IF (Qout(idTimid,ng)) THEN
        scale=1.0_r8
        gtype=gfactor*r2dvar
        status=nf_fwrite2d(ng, iNLM, QCK(ng)%ncid,                      &
     &                     QCK(ng)%Vid(idTimid),                        &
     &                     QCK(ng)%Rindex, gtype,                       &
     &                     LBi, UBi, LBj, UBj, scale,                   &
     &                     GRID(ng) % rmask_full,                       &
     &                     ICE(ng) % ti(:,:,liunw(ng)))
        IF (FoundError(status, nf90_noerr, 1641,                    &
     &                 "ROMS/Utility/wrt_quick.F")) THEN
          IF (Master) THEN
            WRITE (stdout,10) TRIM(Vname(1,idTimid)), QCK(ng)%Rindex
          END IF
          exit_flag=3
          RETURN
        END IF
      END IF
!
! Write out snow thickness.
!
      IF (Qout(idHsno,ng)) THEN
        scale=1.0_r8
        gtype=gfactor*r2dvar
        status=nf_fwrite2d(ng, iNLM, QCK(ng)%ncid,                      &
     &                     QCK(ng)%Vid(idHsno),                         &
     &                     QCK(ng)%Rindex, gtype,                       &
     &                     LBi, UBi, LBj, UBj, scale,                   &
     &                     GRID(ng) % rmask_full,                       &
     &                     ICE(ng) % hsn(:,:,liunw(ng)))
        IF (FoundError(status, nf90_noerr, 1664,                    &
     &                 "ROMS/Utility/wrt_quick.F")) THEN
          IF (Master) THEN
            WRITE (stdout,10) TRIM(Vname(1,idHsno)), QCK(ng)%Rindex
          END IF
          exit_flag=3
          RETURN
        END IF
      END IF
!
! Write out ice age.
!
      IF (Qout(idAgeice,ng)) THEN
        scale=1.0_r8
        gtype=gfactor*r2dvar
        status=nf_fwrite2d(ng, iNLM, QCK(ng)%ncid,                      &
     &                     QCK(ng)%Vid(idAgeice),                       &
     &                     QCK(ng)%Rindex, gtype,                       &
     &                     LBi, UBi, LBj, UBj, scale,                   &
     &                     GRID(ng) % rmask_full,                       &
     &                     ICE(ng) % ageice(:,:,liunw(ng)))
        IF (FoundError(status, nf90_noerr, 1687,                    &
     &                 "ROMS/Utility/wrt_quick.F")) THEN
          IF (Master) THEN
            WRITE (stdout,10) TRIM(Vname(1,idAgeice)), QCK(ng)%Rindex
          END IF
          exit_flag=3
          RETURN
        END IF
      END IF
!
! Write out ice-ocean mass flux
!
      IF (Qout(idIomflx,ng)) THEN
        scale=1.0_r8
        gtype=gfactor*r2dvar
        status=nf_fwrite2d(ng, iNLM, QCK(ng)%ncid,                      &
     &                     QCK(ng)%Vid(idIomflx),                       &
     &                     QCK(ng)%Rindex, gtype,                       &
     &                     LBi, UBi, LBj, UBj, scale,                   &
     &                     GRID(ng) % rmask_full,                       &
     &                     ICE(ng) % io_mflux)
        IF (FoundError(status, nf90_noerr, 1710,                    &
     &                 "ROMS/Utility/wrt_quick.F")) THEN
          IF (Master) THEN
            WRITE (stdout,10) TRIM(Vname(1,idIomflx)), QCK(ng)%Rindex
          END IF
          exit_flag=3
          RETURN
        END IF
      END IF
!
! Write out internal ice stress component 11
!
      IF (Qout(idSig11,ng)) THEN
        scale=1.0_r8
        gtype=gfactor*r2dvar
        status=nf_fwrite2d(ng, iNLM, QCK(ng)%ncid,                      &
     &                     QCK(ng)%Vid(idSig11),                        &
     &                     QCK(ng)%Rindex, gtype,                       &
     &                     LBi, UBi, LBj, UBj, scale,                   &
     &                     GRID(ng) % rmask_full,                       &
     &                     ICE(ng) % sig11(:,:,liunw(ng)))
        IF (FoundError(status, nf90_noerr, 1733,                    &
     &                 "ROMS/Utility/wrt_quick.F")) THEN
          IF (Master) THEN
            WRITE (stdout,10) TRIM(Vname(1,idSig11)), QCK(ng)%Rindex
          END IF
          exit_flag=3
          RETURN
        END IF
      END IF
!
! Write out internal ice stress component 12
!
      IF (Qout(idSig12,ng)) THEN
        scale=1.0_r8
        gtype=gfactor*r2dvar
        status=nf_fwrite2d(ng, iNLM, QCK(ng)%ncid,                      &
     &                     QCK(ng)%Vid(idSig12),                        &
     &                     QCK(ng)%Rindex, gtype,                       &
     &                     LBi, UBi, LBj, UBj, scale,                   &
     &                     GRID(ng) % rmask_full,                       &
     &                     ICE(ng) % sig12(:,:,liunw(ng)))
        IF (FoundError(status, nf90_noerr, 1756,                    &
     &                 "ROMS/Utility/wrt_quick.F")) THEN
          IF (Master) THEN
            WRITE (stdout,10) TRIM(Vname(1,idSig12)), QCK(ng)%Rindex
          END IF
          exit_flag=3
          RETURN
        END IF
      END IF
!
! Write out internal ice stress component 22
!
      IF (Qout(idSig22,ng)) THEN
        scale=1.0_r8
        gtype=gfactor*r2dvar
        status=nf_fwrite2d(ng, iNLM, QCK(ng)%ncid,                      &
     &                     QCK(ng)%Vid(idSig22),                        &
     &                     QCK(ng)%Rindex, gtype,                       &
     &                     LBi, UBi, LBj, UBj, scale,                   &
     &                     GRID(ng) % rmask_full,                       &
     &                     ICE(ng) % sig22(:,:,liunw(ng)))
        IF (FoundError(status, nf90_noerr, 1779,                    &
     &                 "ROMS/Utility/wrt_quick.F")) THEN
          IF (Master) THEN
            WRITE (stdout,10) TRIM(Vname(1,idSig22)), QCK(ng)%Rindex
          END IF
          exit_flag=3
          RETURN
        END IF
      END IF
!
! Write out temperature of molecular sublayer under ice.
!
      IF (Qout(idT0mk,ng)) THEN
        scale=1.0_r8
        gtype=gfactor*r2dvar
        status=nf_fwrite2d(ng, iNLM, QCK(ng)%ncid,                      &
     &                     QCK(ng)%Vid(idT0mk),                         &
     &                     QCK(ng)%Rindex, gtype,                       &
     &                     LBi, UBi, LBj, UBj, scale,                   &
     &                     GRID(ng) % rmask_full,                       &
     &                     ICE(ng) % T0mk)
        IF (FoundError(status, nf90_noerr, 1802,                    &
     &                 "ROMS/Utility/wrt_quick.F")) THEN
          IF (Master) THEN
            WRITE (stdout,10) TRIM(Vname(1,idT0mk)), QCK(ng)%Rindex
          END IF
          exit_flag=3
          RETURN
        END IF
      END IF
!
! Write out salinity of molecular sublayer under ice.
!
      IF (Qout(idS0mk,ng)) THEN
        scale=1.0_r8
        gtype=gfactor*r2dvar
        status=nf_fwrite2d(ng, iNLM, QCK(ng)%ncid,                      &
     &                     QCK(ng)%Vid(idS0mk),                         &
     &                     QCK(ng)%Rindex, gtype,                       &
     &                     LBi, UBi, LBj, UBj, scale,                   &
     &                     GRID(ng) % rmask_full,                       &
     &                     ICE(ng) % S0mk)
        IF (FoundError(status, nf90_noerr, 1825,                    &
     &                 "ROMS/Utility/wrt_quick.F")) THEN
          IF (Master) THEN
            WRITE (stdout,10) TRIM(Vname(1,idS0mk)), QCK(ng)%Rindex
          END IF
          exit_flag=3
          RETURN
        END IF
      END IF
!
! Write out ice-water friction velocity
!
      IF (Qout(idTauiw,ng)) THEN
        scale=1.0_r8
        gtype=gfactor*r2dvar
        status=nf_fwrite2d(ng, iNLM, QCK(ng)%ncid,                      &
     &                     QCK(ng)%Vid(idTauiw),                        &
     &                     QCK(ng)%Rindex, gtype,                       &
     &                     LBi, UBi, LBj, UBj, scale,                   &
     &                     GRID(ng) % rmask_full,                       &
     &                     ICE(ng) % utau_iw)
        IF (FoundError(status, nf90_noerr, 1848,                    &
     &                 "ROMS/Utility/wrt_quick.F")) THEN
          IF (Master) THEN
            WRITE (stdout,10) TRIM(Vname(1,idTauiw)), QCK(ng)%Rindex
          END IF
          exit_flag=3
          RETURN
        END IF
      END IF
!
! Write out ice-water momentum transfer coefficient
!
      IF (Qout(idChuiw,ng)) THEN
        scale=1.0_r8
        gtype=gfactor*r2dvar
        status=nf_fwrite2d(ng, iNLM, QCK(ng)%ncid,                      &
     &                     QCK(ng)%Vid(idChuiw),                        &
     &                     QCK(ng)%Rindex, gtype,                       &
     &                     LBi, UBi, LBj, UBj, scale,                   &
     &                     GRID(ng) % rmask_full,                       &
     &                     ICE(ng) % chu_iw)
        IF (FoundError(status, nf90_noerr, 1871,                    &
     &                 "ROMS/Utility/wrt_quick.F")) THEN
          IF (Master) THEN
            WRITE (stdout,10) TRIM(Vname(1,idChuiw)), QCK(ng)%Rindex
          END IF
          exit_flag=3
          RETURN
        END IF
      END IF
!
!  Write out vertical viscosity coefficient.
!
      IF (Qout(idVvis,ng)) THEN
        scale=1.0_r8
        gtype=gfactor*w3dvar
        status=nf_fwrite3d(ng, iNLM, QCK(ng)%ncid, QCK(ng)%Vid(idVvis), &
     &                     QCK(ng)%Rindex, gtype,                       &
     &                     LBi, UBi, LBj, UBj, 0, N(ng), scale,         &
     &                     GRID(ng) % rmask,                            &
     &                     MIXING(ng) % Akv,                            &
     &                     SetFillVal = .FALSE.)
        IF (FoundError(status, nf90_noerr, 1895,                    &
     &                 "ROMS/Utility/wrt_quick.F")) THEN
          IF (Master) THEN
            WRITE (stdout,10) TRIM(Vname(1,idVvis)), QCK(ng)%Rindex
          END IF
          exit_flag=3
          ioerror=status
          RETURN
        END IF
      END IF
!
!  Write out vertical diffusion coefficient for potential temperature.
!
      IF (Qout(idTdif,ng)) THEN
        scale=1.0_r8
        gtype=gfactor*w3dvar
        status=nf_fwrite3d(ng, iNLM, QCK(ng)%ncid, QCK(ng)%Vid(idTdif), &
     &                     QCK(ng)%Rindex, gtype,                       &
     &                     LBi, UBi, LBj, UBj, 0, N(ng), scale,         &
     &                     GRID(ng) % rmask,                            &
     &                     MIXING(ng) % Akt(:,:,:,itemp),               &
     &                     SetFillVal = .FALSE.)
        IF (FoundError(status, nf90_noerr, 1919,                    &
     &                 "ROMS/Utility/wrt_quick.F")) THEN
          IF (Master) THEN
            WRITE (stdout,10) TRIM(Vname(1,idTdif)), QCK(ng)%Rindex
          END IF
          exit_flag=3
          ioerror=status
          RETURN
        END IF
      END IF
!
!  Write out vertical diffusion coefficient for salinity.
!
      IF (Qout(idSdif,ng)) THEN
        scale=1.0_r8
        gtype=gfactor*w3dvar
        status=nf_fwrite3d(ng, iNLM, QCK(ng)%ncid, QCK(ng)%Vid(idSdif), &
     &                     QCK(ng)%Rindex, gtype,                       &
     &                     LBi, UBi, LBj, UBj, 0, N(ng), scale,         &
     &                     GRID(ng) % rmask,                            &
     &                     MIXING(ng) % Akt(:,:,:,isalt),               &
     &                     SetFillVal = .FALSE.)
        IF (FoundError(status, nf90_noerr, 1944,                    &
     &                 "ROMS/Utility/wrt_quick.F")) THEN
          IF (Master) THEN
            WRITE (stdout,10) TRIM(Vname(1,idSdif)), QCK(ng)%Rindex
          END IF
          exit_flag=3
          ioerror=status
          RETURN
        END IF
      END IF
!
!  Write out surface air pressure.
!
      IF (Qout(idPair,ng)) THEN
        scale=1.0_r8
        gtype=gfactor*r2dvar
        status=nf_fwrite2d(ng, iNLM, QCK(ng)%ncid, QCK(ng)%Vid(idPair), &
     &                     QCK(ng)%Rindex, gtype,                       &
     &                     LBi, UBi, LBj, UBj, scale,                   &
     &                     GRID(ng) % rmask,                            &
     &                     FORCES(ng) % Pair)
        IF (FoundError(status, nf90_noerr, 2020,                    &
     &                 "ROMS/Utility/wrt_quick.F")) THEN
          IF (Master) THEN
            WRITE (stdout,10) TRIM(Vname(1,idPair)), QCK(ng)%Rindex
          END IF
          exit_flag=3
          ioerror=status
          RETURN
        END IF
      END IF
!
!  Write out surface air temperature.
!
      IF (Qout(idTair,ng)) THEN
        scale=1.0_r8
        gtype=gfactor*r2dvar
        status=nf_fwrite2d(ng, iNLM, QCK(ng)%ncid, QCK(ng)%Vid(idTair), &
     &                     QCK(ng)%Rindex, gtype,                       &
     &                     LBi, UBi, LBj, UBj, scale,                   &
     &                     GRID(ng) % rmask,                            &
     &                     FORCES(ng) % Tair)
        IF (FoundError(status, nf90_noerr, 2045,                    &
     &                 "ROMS/Utility/wrt_quick.F")) THEN
          IF (Master) THEN
            WRITE (stdout,10) TRIM(Vname(1,idTair)), QCK(ng)%Rindex
          END IF
          exit_flag=3
          ioerror=status
          RETURN
        END IF
      END IF
!
!  Write out surface winds.
!
      IF (Qout(idUair,ng)) THEN
        scale=1.0_r8
        gtype=gfactor*r2dvar
        status=nf_fwrite2d(ng, iNLM, QCK(ng)%ncid, QCK(ng)%Vid(idUair), &
     &                     QCK(ng)%Rindex, gtype,                       &
     &                     LBi, UBi, LBj, UBj, scale,                   &
     &                     GRID(ng) % rmask,                            &
     &                     FORCES(ng) % Uwind)
        IF (FoundError(status, nf90_noerr, 2070,                    &
     &                 "ROMS/Utility/wrt_quick.F")) THEN
          IF (Master) THEN
            WRITE (stdout,10) TRIM(Vname(1,idUair)), QCK(ng)%Rindex
          END IF
          exit_flag=3
          ioerror=status
          RETURN
        END IF
      END IF
      IF (Qout(idVair,ng)) THEN
        scale=1.0_r8
        gtype=gfactor*r2dvar
        status=nf_fwrite2d(ng, iNLM, QCK(ng)%ncid, QCK(ng)%Vid(idVair), &
     &                     QCK(ng)%Rindex, gtype,                       &
     &                     LBi, UBi, LBj, UBj, scale,                   &
     &                     GRID(ng) % rmask,                            &
     &                     FORCES(ng) % Vwind)
        IF (FoundError(status, nf90_noerr, 2091,                    &
     &                 "ROMS/Utility/wrt_quick.F")) THEN
          IF (Master) THEN
            WRITE (stdout,10) TRIM(Vname(1,idVair)), QCK(ng)%Rindex
          END IF
          exit_flag=3
          ioerror=status
          RETURN
        END IF
      END IF
!
!  Write out 2D Eastward surface wind (m/s) at RHO-points.
!
      IF (Qout(idUairE,ng).and.Qout(idVairN,ng)) THEN
        IF (.not.allocated(Ur2d)) THEN
          allocate (Ur2d(LBi:UBi,LBj:UBj))
            Ur2d(LBi:UBi,LBj:UBj)=0.0_r8
        END IF
        IF (.not.allocated(Vr2d)) THEN
          allocate (Vr2d(LBi:UBi,LBj:UBj))
            Vr2d(LBi:UBi,LBj:UBj)=0.0_r8
        END IF
        CALL uv_rotate2d (ng, TILE, .FALSE., .TRUE.,                    &
     &                    LBi, UBi, LBj, UBj,                           &
     &                    GRID(ng) % CosAngler,                         &
     &                    GRID(ng) % SinAngler,                         &
     &                    GRID(ng) % rmask_full,                        &
     &                    FORCES(ng) % Uwind,                           &
     &                    FORCES(ng) % Vwind,                           &
     &                    Ur2d, Vr2d)
        scale=1.0_r8
        gtype=gfactor*r2dvar
        status=nf_fwrite2d(ng, iNLM, QCK(ng)%ncid,                      &
     &                     QCK(ng)%Vid(idUairE),                        &
     &                     QCK(ng)%Rindex, gtype,                       &
     &                     LBi, UBi, LBj, UBj, scale,                   &
     &                     GRID(ng) % rmask_full,                       &
     &                     Ur2d)
        IF (FoundError(status, nf90_noerr, 2133,                    &
     &                 "ROMS/Utility/wrt_quick.F")) THEN
          IF (Master) THEN
            WRITE (stdout,10) TRIM(Vname(1,idUairE)), QCK(ng)%Rindex
          END IF
          exit_flag=3
          ioerror=status
          RETURN
        END IF
!
!  Write out 2D Northward surface wind (m/s) at RHO-points.
!
        scale=1.0_r8
        gtype=gfactor*r2dvar
        status=nf_fwrite2d(ng, iNLM, QCK(ng)%ncid,                      &
     &                     QCK(ng)%Vid(idVairN),                        &
     &                     QCK(ng)%Rindex, gtype,                       &
     &                     LBi, UBi, LBj, UBj, scale,                   &
     &                     GRID(ng) % rmask_full,                       &
     &                     Vr2d)
        IF (FoundError(status, nf90_noerr, 2155,                    &
     &                 "ROMS/Utility/wrt_quick.F")) THEN
          IF (Master) THEN
            WRITE (stdout,10) TRIM(Vname(1,idVairN)), QCK(ng)%Rindex
          END IF
          exit_flag=3
          ioerror=status
          RETURN
        END IF
      END IF
!
!  Write out surface active traces fluxes.
!
      DO itrc=1,NAT
        IF (Qout(idTsur(itrc),ng)) THEN
          IF (itrc.eq.itemp) THEN
            scale=rho0*Cp                   ! Celsius m/s to W/m2
          ELSE IF (itrc.eq.isalt) THEN
            scale=1.0_r8
          END IF
          gtype=gfactor*r2dvar
          status=nf_fwrite2d(ng, iNLM, QCK(ng)%ncid,                    &
     &                       QCK(ng)%Vid(idTsur(itrc)),                 &
     &                       QCK(ng)%Rindex, gtype,                     &
     &                       LBi, UBi, LBj, UBj, scale,                 &
     &                       GRID(ng) % rmask,                          &
     &                       FORCES(ng) % stflx(:,:,itrc))
          IF (FoundError(status, nf90_noerr, 2190,                  &
     &                   "ROMS/Utility/wrt_quick.F")) THEN
            IF (Master) THEN
              WRITE (stdout,10) TRIM(Vname(1,idTsur(itrc))),            &
     &                          QCK(ng)%Rindex
            END IF
            exit_flag=3
            ioerror=status
            RETURN
          END IF
        END IF
      END DO
!
!  Write out latent heat flux.
!
      IF (Qout(idLhea,ng)) THEN
        scale=rho0*Cp
        gtype=gfactor*r2dvar
        status=nf_fwrite2d(ng, iNLM, QCK(ng)%ncid, QCK(ng)%Vid(idLhea), &
     &                     QCK(ng)%Rindex, gtype,                       &
     &                     LBi, UBi, LBj, UBj, scale,                   &
     &                     GRID(ng) % rmask,                            &
     &                     FORCES(ng) % lhflx)
        IF (FoundError(status, nf90_noerr, 2217,                    &
     &                 "ROMS/Utility/wrt_quick.F")) THEN
          IF (Master) THEN
            WRITE (stdout,10) TRIM(Vname(1,idLhea)), QCK(ng)%Rindex
          END IF
          exit_flag=3
          ioerror=status
          RETURN
        END IF
      END IF
!
!  Write out sensible heat flux.
!
      IF (Qout(idShea,ng)) THEN
        scale=rho0*Cp
        gtype=gfactor*r2dvar
        status=nf_fwrite2d(ng, iNLM, QCK(ng)%ncid, QCK(ng)%Vid(idShea), &
     &                     QCK(ng)%Rindex, gtype,                       &
     &                     LBi, UBi, LBj, UBj, scale,                   &
     &                     GRID(ng) % rmask,                            &
     &                     FORCES(ng) % shflx)
        IF (FoundError(status, nf90_noerr, 2240,                    &
     &                 "ROMS/Utility/wrt_quick.F")) THEN
          IF (Master) THEN
            WRITE (stdout,10) TRIM(Vname(1,idShea)), QCK(ng)%Rindex
          END IF
          exit_flag=3
          ioerror=status
          RETURN
        END IF
      END IF
!
!  Write out longwave radiation flux.
!
      IF (Qout(idLrad,ng)) THEN
        scale=rho0*Cp
        gtype=gfactor*r2dvar
        status=nf_fwrite2d(ng, iNLM, QCK(ng)%ncid, QCK(ng)%Vid(idLrad), &
     &                     QCK(ng)%Rindex, gtype,                       &
     &                     LBi, UBi, LBj, UBj, scale,                   &
     &                     GRID(ng) % rmask,                            &
     &                     FORCES(ng) % lrflx)
        IF (FoundError(status, nf90_noerr, 2263,                    &
     &                 "ROMS/Utility/wrt_quick.F")) THEN
          IF (Master) THEN
            WRITE (stdout,10) TRIM(Vname(1,idLrad)), QCK(ng)%Rindex
          END IF
          exit_flag=3
          ioerror=status
          RETURN
        END IF
      END IF
!
!  Write out evaporation rate (kg/m2/s).
!
      IF (Qout(idevap,ng)) THEN
        scale=1.0_r8
        gtype=gfactor*r2dvar
        status=nf_fwrite2d(ng, iNLM, QCK(ng)%ncid, QCK(ng)%Vid(idevap), &
     &                     QCK(ng)%Rindex, gtype,                       &
     &                     LBi, UBi, LBj, UBj, scale,                   &
     &                     GRID(ng) % rmask,                            &
     &                     FORCES(ng) % evap)
        IF (FoundError(status, nf90_noerr, 2287,                    &
     &                 "ROMS/Utility/wrt_quick.F")) THEN
          IF (Master) THEN
            WRITE (stdout,10) TRIM(Vname(1,idevap)), QCK(ng)%Rindex
          END IF
          exit_flag=3
          ioerror=status
          RETURN
        END IF
      END IF
!
!  Write out precipitation rate (kg/m2/s).
!
      IF (Qout(idrain,ng)) THEN
        scale=1.0_r8
        gtype=gfactor*r2dvar
        status=nf_fwrite2d(ng, iNLM, QCK(ng)%ncid, QCK(ng)%Vid(idrain), &
     &                     QCK(ng)%Rindex, gtype,                       &
     &                     LBi, UBi, LBj, UBj, scale,                   &
     &                     GRID(ng) % rmask,                            &
     &                     FORCES(ng) % rain)
        IF (FoundError(status, nf90_noerr, 2310,                    &
     &                 "ROMS/Utility/wrt_quick.F")) THEN
          IF (Master) THEN
            WRITE (stdout,10) TRIM(Vname(1,idrain)), QCK(ng)%Rindex
          END IF
          exit_flag=3
          ioerror=status
          RETURN
        END IF
      END IF
!
!  Write out E-P (m/s).
!
      IF (Qout(idEmPf,ng)) THEN
        scale=1.0_r8
        gtype=gfactor*r2dvar
        status=nf_fwrite2d(ng, iNLM, QCK(ng)%ncid, QCK(ng)%Vid(idEmPf), &
     &                     QCK(ng)%Rindex, gtype,                       &
     &                     LBi, UBi, LBj, UBj, scale,                   &
     &                     GRID(ng) % rmask,                            &
     &                     FORCES(ng) % EminusP)
        IF (FoundError(status, nf90_noerr, 2336,                    &
     &                 "ROMS/Utility/wrt_quick.F")) THEN
          IF (Master) THEN
            WRITE (stdout,10) TRIM(Vname(1,idEmPf)), QCK(ng)%Rindex
          END IF
          exit_flag=3
          ioerror=status
          RETURN
        END IF
      END IF
!
!  Write out shortwave radiation flux.
!
      IF (Qout(idSrad,ng)) THEN
        scale=rho0*Cp
        gtype=gfactor*r2dvar
        status=nf_fwrite2d(ng, iNLM, QCK(ng)%ncid, QCK(ng)%Vid(idSrad), &
     &                     QCK(ng)%Rindex, gtype,                       &
     &                     LBi, UBi, LBj, UBj, scale,                   &
     &                     GRID(ng) % rmask,                            &
     &                     FORCES(ng) % srflx)
        IF (FoundError(status, nf90_noerr, 2361,                    &
     &                 "ROMS/Utility/wrt_quick.F")) THEN
          IF (Master) THEN
            WRITE (stdout,10) TRIM(Vname(1,idSrad)), QCK(ng)%Rindex
          END IF
          exit_flag=3
          ioerror=status
          RETURN
        END IF
      END IF
!
!  Write out surface U-momentum stress.
!
      IF (Qout(idUsms,ng)) THEN
        scale=rho0                          ! m2/s2 to Pa
        gtype=gfactor*u2dvar
        status=nf_fwrite2d(ng, iNLM, QCK(ng)%ncid, QCK(ng)%Vid(idUsms), &
     &                     QCK(ng)%Rindex, gtype,                       &
     &                     LBi, UBi, LBj, UBj, scale,                   &
     &                     GRID(ng) % umask,                            &
     &                     FORCES(ng) % sustr)
        IF (FoundError(status, nf90_noerr, 2386,                    &
     &                 "ROMS/Utility/wrt_quick.F")) THEN
          IF (Master) THEN
            WRITE (stdout,10) TRIM(Vname(1,idUsms)), QCK(ng)%Rindex
          END IF
          exit_flag=3
          ioerror=status
          RETURN
        END IF
      END IF
!
!  Write out surface V-momentum stress.
!
      IF (Qout(idVsms,ng)) THEN
        scale=rho0
        gtype=gfactor*v2dvar
        status=nf_fwrite2d(ng, iNLM, QCK(ng)%ncid, QCK(ng)%Vid(idVsms), &
     &                     QCK(ng)%Rindex, gtype,                       &
     &                     LBi, UBi, LBj, UBj, scale,                   &
     &                     GRID(ng) % vmask,                            &
     &                     FORCES(ng) % svstr)
        IF (FoundError(status, nf90_noerr, 2409,                    &
     &                 "ROMS/Utility/wrt_quick.F")) THEN
          IF (Master) THEN
            WRITE (stdout,10) TRIM(Vname(1,idVsms)), QCK(ng)%Rindex
          END IF
          exit_flag=3
          ioerror=status
          RETURN
        END IF
      END IF
!
!  Write out bottom U-momentum stress.
!
      IF (Qout(idUbms,ng)) THEN
        scale=-rho0
        gtype=gfactor*u2dvar
        status=nf_fwrite2d(ng, iNLM, QCK(ng)%ncid, QCK(ng)%Vid(idUbms), &
     &                     QCK(ng)%Rindex, gtype,                       &
     &                     LBi, UBi, LBj, UBj, scale,                   &
     &                     GRID(ng) % umask,                            &
     &                     FORCES(ng) % bustr)
        IF (FoundError(status, nf90_noerr, 2432,                    &
     &                 "ROMS/Utility/wrt_quick.F")) THEN
          IF (Master) THEN
            WRITE (stdout,10) TRIM(Vname(1,idUbms)), QCK(ng)%Rindex
          END IF
          exit_flag=3
          ioerror=status
          RETURN
        END IF
      END IF
!
!  Write out bottom V-momentum stress.
!
      IF (Qout(idVbms,ng)) THEN
        scale=-rho0
        gtype=gfactor*v2dvar
        status=nf_fwrite2d(ng, iNLM, QCK(ng)%ncid, QCK(ng)%Vid(idVbms), &
     &                     QCK(ng)%Rindex, gtype,                       &
     &                     LBi, UBi, LBj, UBj, scale,                   &
     &                     GRID(ng) % vmask,                            &
     &                     FORCES(ng) % bvstr)
        IF (FoundError(status, nf90_noerr, 2455,                    &
     &                 "ROMS/Utility/wrt_quick.F")) THEN
          IF (Master) THEN
            WRITE (stdout,10) TRIM(Vname(1,idVbms)), QCK(ng)%Rindex
          END IF
          exit_flag=3
          ioerror=status
          RETURN
        END IF
      END IF
!
!-----------------------------------------------------------------------
!  Synchronize quicksave NetCDF file to disk to allow other processes
!  to access data immediately after it is written.
!-----------------------------------------------------------------------
!
      CALL netcdf_sync (ng, iNLM, QCK(ng)%name, QCK(ng)%ncid)
      IF (FoundError(exit_flag, NoError, 3420,                      &
     &               "ROMS/Utility/wrt_quick.F")) RETURN
      IF (Master) WRITE (stdout,20) kstp(ng), nrhs(ng), QCK(ng)%Rindex
!
  10  FORMAT (/,' WRT_QUICK - error while writing variable: ',a,/,13x,  &
     &        'into quicksave NetCDF file for time record: ',i4)
  20  FORMAT (6x,'WRT_QUICK   - wrote quicksave', t39,                  &
     &        'fields (Index=',i1,',',i1,') in record = ',i7.7)
      RETURN
      END SUBROUTINE wrt_quick
