!======================================================================
! CROCO is a branch of ROMS developped at IRD, INRIA, 
! Ifremer, CNRS and Univ. Toulouse III  in France
! The two other branches from UCLA (Shchepetkin et al)
! and Rutgers University (Arango et al) are under MIT/X style license.
! CROCO specific routines (nesting) are under CeCILL-C license.
!
! CROCO website : http://www.croco-ocean.org
!======================================================================
!
#ifndef AVRH
# include "cppdefs.h"
!
      subroutine def_his (ncid, total_rec, ierr)
#else
      subroutine def_avg (ncid, total_rec, ierr)
#endif
!
! Create/open averages/history netCDF file. In the case when a new
! netCDF file is created, define all variables, their dimensions and
! attributes. In the case when a previously existing netCDF file is
! to be opened for addition of new data, verify that all dimensions
! of the file are consistent with the present model configuration
! and all necessary variables exist. Save netCDF IDs for all needed
! variables. Also determine size of the unlimited dimension.
!
! The difference between def_his and def_avg is as follows: they
! have different netCDF file name (hisname/avgname); netCDF file ID
! (passed as argument); time record index (hisindx/avgindx); array
! of switches which variables to write (wrthis/wrtavg); and different
! sets of netCDF variable IDs (hisTime...hisHbl/avgTime...avgHbl);
! and the first attribute of each variable, long_name, has prefix
! 'averaged'. Because most of the code is identical for both
! routines, the second one is generated from the first entirely
! by CPP.
!
#if defined MUSTANG
      USE comsubstance, ONLY : nvpc
      USE comMUSTANG, ONLY : hisMust, avgMust, outMust, outMust_nbvar,
     &     out2DMust, out3DsedMust, vname_Must, nk_nivsed_out
#endif
#if defined OBSTRUCTION
      USE com_OBSTRUCTIONS, ONLY : hisObst, avgObst, 
     &     obst_nbvar, outObst, out2DObst, vname_obst
#endif


      implicit none
      logical create_new_file
      integer ncid, total_rec, ierr, rec, lstr,lvar,lenstr, timedim
     &      , r2dgrd(3),  u2dgrd(3), v2dgrd(3), auxil(2), checkdims
     &      , b3dgrd(4)
#ifdef NC4PAR
     &      , csize,cmode
#endif
#if defined SEDIMENT || defined BBL
     &      , indxWrk
#endif
#ifdef SOLVE3D
     &      , r3dgrd(4),  u3dgrd(4), v3dgrd(4), w3dgrd(4), itrc,itrv
#endif
#ifdef MRL_WCI
     &      , iwkb
#endif
#if defined MUSTANG || defined OBSTRUCTION
     &      , indx
#endif
#include "param.h"
#include "scalars.h"
#include "ncscrum.h"
#include "netcdf.inc"
#ifdef NC4PAR
# include "mpi_cpl.h"
      include 'mpif.h'
#endif
      character*70 text

#ifndef AVRH
# define ncname hisname
# define rec_per_file nrpfhis
# define wrt wrthis
# define vidTime hisTime
# define vidTime2 hisTime2
# define vidTstep hisTstep
# define vidZ hisZ
# define vidUb hisUb
# define vidVb hisVb
# define vidU hisU
# define vidV hisV
# define vidT hisT
# define vidR hisR
# define vidbvf hisbvf
# define vidO hisO
# define vidW hisW
# define vidVisc hisVisc
# define vidDiff hisDiff
# define vidAkv hisAkv
# define vidAkt hisAkt
# define vidAks hisAks
# define vidHbl hisHbl
# define vidHbbl hisHbbl
# define vidTke hisTke
# define vidGls hisGls
# define vidLsc hisLsc
# define vidHm hisHm
# define vidWAVE hisWAVE
# define vidSUP hisSUP
# define vidUST2D hisUST2D
# define vidVST2D hisVST2D
# define vidUST hisUST
# define vidVST hisVST
# define vidWST hisWST
# define vidAkb hisAkb
# define vidAkw hisAkw
# define vidKVF hisKVF
# define vidCALP hisCALP
# define vidKAPS hisKAPS
# define vidHel hisHel
# define vidChC hisChC
# define vidU10 hisU10
# define vidKvO2 hisKvO2
# define vidO2sat hisO2sat
# define vidAOU hisAOU
# define vidwind10 hiswind10
# define vidSed hisSed
# define vidBBL hisBBL
# define vidBostr hisBostr
# define vidBustr hisBustr
# define vidBvstr hisBvstr
# define vidWstr hisWstr
# define vidUWstr hisUWstr
# define vidVWstr hisVWstr
# define vidShflx hisShflx
# ifdef SALINITY
#  define vidSwflx hisSwflx
# endif
# define vidShflx_rsw hisShflx_rsw
# ifdef BULK_FLUX
#  define vidShflx_rlw hisShflx_rlw
#  define vidShflx_lat hisShflx_lat
#  define vidShflx_sen hisShflx_sen
# endif
# if defined TEMPERATURE && defined BHFLUX
#  define vidBhflx hisBhflx
# endif
# if defined SALINITY && defined BWFLUX
#  define vidBwflx hisBwflx
# endif
# ifdef SST_SKIN
#  define vidSST_skin hisSST_skin
# endif
# define vidMust hisMust
# define vidObst hisObst

#else

# define ncname avgname
# define rec_per_file nrpfavg
# define wrt wrtavg
# define vidTime avgTime
# define vidTime2 avgTime2
# define vidTstep avgTstep
# define vidZ avgZ
# define vidUb avgUb
# define vidVb avgVb
# define vidU avgU
# define vidV avgV
# define vidT avgT
# define vidR avgR
# define vidbvf avgbvf
# define vidO avgO
# define vidW avgW
# define vidVisc avgVisc
# define vidDiff avgDiff
# define vidAkv avgAkv
# define vidAkt avgAkt
# define vidAks avgAks
# define vidHbl avgHbl
# define vidHbbl avgHbbl
# define vidTke avgTke
# define vidGls avgGls
# define vidLsc avgLsc
# define vidHel avgHel
# define vidChC avgChC
# define vidU10 avgU10
# define vidKvO2 avgKvO2
# define vidO2sat avgO2sat
# define vidAOU avgAOU
# define vidwind10 avgwind10
# define vidSed avgSed
# define vidBBL avgBBL
# define vidBostr avgBostr
# define vidBustr avgBustr
# define vidBvstr avgBvstr
# define vidWstr avgWstr
# define vidUWstr avgUWstr
# define vidVWstr avgVWstr
# define vidHm avgHm
# define vidWAVE avgWAVE
# define vidSUP avgSUP
# define vidUST2D avgUST2D
# define vidVST2D avgVST2D
# define vidUST avgUST
# define vidWST avgWST
# define vidVST avgVST
# define vidAkb avgAkb
# define vidAkw avgAkw
# define vidKVF avgKVF
# define vidCALP avgCALP
# define vidKAPS avgKAPS
# define vidShflx avgShflx
# define vidSwflx avgSwflx
# define vidShflx_rsw avgShflx_rsw
# if defined BULK_FLUX
#  define vidShflx_rlw avgShflx_rlw
#  define vidShflx_lat avgShflx_lat
#  define vidShflx_sen avgShflx_sen
# endif
# ifdef BHFLUX
#  define vidBhflx avgBhflx
# endif
# if defined BWFLUX && defined SALINITY
#  define vidBwflx avgBwflx
# endif
# ifdef SST_SKIN
#  define vidSST_skin avgSST_skin
# endif

# define vidMust avgMust
# define vidObst avgObst
#endif

# ifdef M3FAST_HIS
#  define vidRnbq hisRnbq
#  define vidUnbq hisUnbq
#  define vidVnbq hisVnbq
#  define vidWnbq hisWnbq
#  define vidCnbq hisCnbq
!
! M3FAST-kernel related variables:
!
!****** Init parameters ******
!
        wrt(indxRnbq)=.true.
        wrt(indxUnbq)=.true.
        wrt(indxVnbq)=.true.
        wrt(indxWnbq)=.true.
        wrt(indxCnbq)=.true.
        
        vname(1,indxRnbq)="Rnbq"
        vname(2,indxRnbq)=vname(2,indxR)
        vname(3,indxRnbq)=vname(3,indxR)
        vname(4,indxRnbq)=vname(4,indxR)
        vname(1,indxUnbq)="Unbq"
        vname(2,indxUnbq)=vname(2,indxU)
        vname(3,indxUnbq)=vname(3,indxU)
        vname(4,indxUnbq)=vname(4,indxU)
        vname(1,indxVnbq)="Vnbq"
        vname(2,indxVnbq)=vname(2,indxV)
        vname(3,indxVnbq)=vname(3,indxV)
        vname(4,indxVnbq)=vname(4,indxV)
        vname(1,indxWnbq)="Wnbq"
        vname(2,indxWnbq)=vname(2,indxW)
        vname(3,indxWnbq)=vname(3,indxW)
        vname(4,indxWnbq)=vname(4,indxW)
        vname(1,indxCnbq)="Cnbq"
        vname(2,indxCnbq)=vname(1,indxCnbq)
        vname(3,indxCnbq)="ms-1"
        vname(4,indxCnbq)="Sound speed, scalar."
!
# endif
!
! Put time record index into file name. In  the case when model
! output is to be arranged into sequence of named files, the naming
! convention is as follows: 'his_root.INDEX.[MPI_node.]nc', where
! INDEX is an integer number such that (i) it is divisible by the
! specified number of records per file; and (ii)
!
!      INDEX + record_within_the_file = total_record
!
! where, 1 =< record_within_the_file =< records_per_file, so that
! total_record changes continuously throughout the sequence of files.
!
      ierr=0
      lstr=lenstr(ncname)
      if (rec_per_file.gt.0) then
        lvar=total_rec-(1+mod(total_rec-1, rec_per_file))
        call insert_time_index (ncname, lstr, lvar, ierr)
        if (ierr .ne. 0) goto 99
      endif
!
! Decide whether to create a new file, or open existing one.
! Overall the whole code below is organized into 3-way switch,
!
! 10  if (create_new_file) then
!        .... create new file, save netCDF ids for all variables;
!     elseif (ncid.eq.-1) then
!        .... try to open existing file and check its dimensions
!       if (cannot be opened or rejected) then
!         create_new_file=.true.
!         goto 10
!       endif   and prepare
!        .... prepare the file for adding new data,
!        .... find and save netCDF ids for all variables
!     else
!        .... just open, no checking, all ids are assumed to be
!        .... already known (MPI single file output only).
!     endif
!
! which is designed to implement flexible opening policy:
! if ldefhis=.true., it forces creation of a new file [if the
! file already exists, it will be overwritten]; on the other hand,
! ldefhis=.false., it is assumed that the file already exists and
! an attempt to open it is made; if the attempt is successful, the
! file is prepared for appending hew data; if it fails, a new file
! is created.
!

#ifdef DO_NOT_OVERWRITE
      if (lvar.eq.total_rec-1) then
        create_new_file=ldefhis
      else
        create_new_file=.false.
      endif
#else
      create_new_file=ldefhis
#endif
      if (ncid.ne.-1) create_new_file=.false.
#if defined MPI & !defined PARALLEL_FILES & !defined NC4PAR
      if (mynode.gt.0) create_new_file=.false.
#endif
!
!=================================================================
!
! Create new history/averages file:    Put global attributes
! ====== === ======= ======== =====    and define all variables.
!
!=================================================================
!
  10  if (create_new_file) then
#ifndef NC4PAR
        ierr=nf_create(ncname(1:lstr),NF_CLOBBER, ncid)
#else
        cmode = ior(nf_netcdf4,nf_classic_model)
        cmode = ior(cmode, nf_mpiio)
        csize=xi_rho*eta_rho/NNODES
        MPI_master_only write(stdout,*)
     &       'CREATE HIS/AVG NC4 PARALLEL FILE'
        ierr=nf_create_par(ncname(1:lstr),cmode,
     &        MPI_COMM_WORLD,MPI_INFO_NULL,ncid)
#endif
        if (ierr .ne. nf_noerr) then
          MPI_master_only write(stdout,'(/3(1x,A)/)')
     &           'ERROR in def_his/avg:',
     &           'Cannot create netCDF file:', ncname(1:lstr)
          goto 99                                         !--> ERROR
        endif
        if (rec_per_file.eq.0) total_rec=0
!
! Put global attributes.
! --- ------ -----------
!
        call put_global_atts (ncid, ierr)
!
! Define dimensions of staggered fields.
! ------ ---------- -- --------- -------
!
        ierr=nf_def_dim (ncid, 'xi_rho',   xi_rho,   r2dgrd(1))
        ierr=nf_def_dim (ncid, 'xi_u',     xi_u,     u2dgrd(1))
        ierr=nf_def_dim (ncid, 'eta_rho',  eta_rho,  r2dgrd(2))
        ierr=nf_def_dim (ncid, 'eta_v',    eta_v,    v2dgrd(2))
#ifdef SEDIMENT
        ierr=nf_def_dim (ncid, 's_b',      NLAY,     b3dgrd(3))
#endif
#ifdef  MUSTANG
        ierr=nf_def_dim (ncid, 's_b',      nk_nivsed_out,  b3dgrd(3))
#endif
#ifdef SOLVE3D
        ierr=nf_def_dim (ncid, 's_rho',    N,        r3dgrd(3))
        ierr=nf_def_dim (ncid, 's_w',      N+1,      w3dgrd(3))
#endif
        ierr=nf_def_dim (ncid, 'time', nf_unlimited, timedim)
#if defined OUTPUTS_SURFACE && ! defined XIOS
        ierr=nf_def_dim (ncid, 'auxil',    6,        auxil(1))
#else
        ierr=nf_def_dim (ncid, 'auxil',    4,        auxil(1))
#endif
        auxil(2)=timedim

        r2dgrd(3)=timedim           ! Free surface

        u2dgrd(2)=r2dgrd(2)         ! 2D UBAR-type
        u2dgrd(3)=timedim

        v2dgrd(1)=r2dgrd(1)         ! 2D VBAR-type
        v2dgrd(3)=timedim

        b3dgrd(1)=r2dgrd(1)         !
        b3dgrd(2)=r2dgrd(2)         ! 3D BED-type
        b3dgrd(4)=timedim           !

#ifdef SOLVE3D
        r3dgrd(1)=r2dgrd(1)         !
        r3dgrd(2)=r2dgrd(2)         ! 3D RHO-type
        r3dgrd(4)=timedim           !

        u3dgrd(1)=u2dgrd(1)         !
        u3dgrd(2)=r2dgrd(2)         ! 3D U-type
        u3dgrd(3)=r3dgrd(3)         !
        u3dgrd(4)=timedim

        v3dgrd(1)=r2dgrd(1)         !
        v3dgrd(2)=v2dgrd(2)         ! 3D V-type
        v3dgrd(3)=r3dgrd(3)         !
        v3dgrd(4)=timedim

        w3dgrd(1)=r2dgrd(1)         !
        w3dgrd(2)=r2dgrd(2)         ! 3D W-type
        w3dgrd(4)=timedim           !
#endif
#if (defined PUT_GRID_INTO_HISTORY && !defined AVRH)\
 || (defined PUT_GRID_INTO_AVERAGES && defined AVRH)
!
! Define grid variables.
! ------ ---- ----------
!
        if (total_rec.le.1) then
# ifdef SOLVE3D
          call def_grid_3d(ncid, r2dgrd, u2dgrd, v2dgrd
     &                                  ,r3dgrd, w3dgrd)
# else
          call def_grid_2d(ncid, r2dgrd, u2dgrd, v2dgrd)
# endif
        endif
#endif

!
! Define evolving model variables.
! ------ -------- ----- ----------
!
! Time step number and time record indices:
!
        ierr=nf_def_var (ncid, 'time_step', nf_int, 2, auxil,
     &                                                 vidTstep)
#ifdef NC4PAR
        ierr=nf_var_par_access(ncid,vidTstep,nf_collective)
#endif
        ierr=nf_put_att_text (ncid, vidTstep, 'long_name', 48,
     &       'time step and record numbers from initialization')
!
! Time.
!
        lvar=lenstr(vname(1,indxTime))
        ierr=nf_def_var (ncid, vname(1,indxTime)(1:lvar),
     &                            NF_DOUBLE, 1, timedim, vidTime)

#ifdef NC4PAR
        ierr=nf_var_par_access(ncid,vidTime,nf_collective)
#endif
#ifndef AVRH
        lvar=lenstr(vname(2,indxTime))
        ierr=nf_put_att_text (ncid, vidTime, 'long_name', lvar,
     &                                vname(2,indxTime)(1:lvar))
#else
        text='averaged '/ /vname(2,indxTime)
        lvar=lenstr(text)
        ierr=nf_put_att_text (ncid, vidTime, 'long_name', lvar,
     &                                             text(1:lvar))
#endif
        lvar=lenstr(vname(3,indxTime))
        ierr=nf_put_att_text (ncid, vidTime, 'units',  lvar,
     &                                vname(3,indxTime)(1:lvar))
        lvar=lenstr(vname(4,indxTime))
        ierr=nf_put_att_text (ncid, vidTime, 'field',  lvar,
     &                                vname(4,indxTime)(1:lvar))

        call nf_add_attribute(ncid, vidTime, vname(:,indxTime), 5,
     &       NF_FOUT, ierr)

!
! Time2.
!
        lvar=lenstr(vname(1,indxTime2))
        ierr=nf_def_var (ncid, vname(1,indxTime2)(1:lvar),
     &                            NF_DOUBLE, 1, timedim, vidTime2)
#ifdef NC4PAR
        ierr=nf_var_par_access(ncid,vidTime2,nf_collective)
#endif
#ifndef AVRH
        lvar=lenstr(vname(2,indxTime2))
        ierr=nf_put_att_text (ncid, vidTime2, 'long_name', lvar,
     &                                vname(2,indxTime2)(1:lvar))
#else
        text='averaged '/ /vname(2,indxTime2)
        lvar=lenstr(text)
        ierr=nf_put_att_text (ncid, vidTime2, 'long_name', lvar,
     &                                             text(1:lvar))
#endif
        lvar=lenstr(vname(3,indxTime2))
        ierr=nf_put_att_text (ncid, vidTime2, 'units',  lvar,
     &                                vname(3,indxTime2)(1:lvar))
        lvar=lenstr(vname(4,indxTime2))
        ierr=nf_put_att_text (ncid, vidTime2, 'field',  lvar,
     &                                vname(4,indxTime2)(1:lvar))

        call nf_add_attribute(ncid, vidTime2, vname(:,indxTime2), 5,
     &       NF_FOUT, ierr)

!=======================================================================
!
! Define primary variables
!
!=======================================================================

!
!-----------------------------------------------------------------------
! Free-surface.
!-----------------------------------------------------------------------
!
        if (wrt(indxZ)) then
          lvar=lenstr(vname(1,indxZ))
          ierr=nf_def_var (ncid, vname(1,indxZ)(1:lvar),
     &                              NF_FOUT, 3, r2dgrd, vidZ)
#ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidZ,nf_collective)
#endif
#ifndef AVRH
          lvar=lenstr(vname(2,indxZ))
          ierr=nf_put_att_text (ncid, vidZ, 'long_name', lvar,
     &                                  vname(2,indxZ)(1:lvar))
#else
          text='averaged '/ /vname(2,indxZ)
          lvar=lenstr(text)
          ierr=nf_put_att_text (ncid, vidZ, 'long_name', lvar,
     &                                            text(1:lvar))
#endif
          lvar=lenstr(vname(3,indxZ))
          ierr=nf_put_att_text (ncid, vidZ, 'units',     lvar,
     &                                  vname(3,indxZ)(1:lvar))
          lvar=lenstr(vname(4,indxZ))
          ierr=nf_put_att_text (ncid, vidZ, 'field',     lvar,
     &                                  vname(4,indxZ)(1:lvar))

          call nf_add_attribute(ncid, vidZ, vname(:,indxZ), 
     &         5, NF_FOUT, ierr)

        endif
!
!-----------------------------------------------------------------------
! 2D momenta in XI- and ETA-directions.
!-----------------------------------------------------------------------
!
        if (wrt(indxUb)) then
          lvar=lenstr(vname(1,indxUb))
          ierr=nf_def_var (ncid, vname(1,indxUb)(1:lvar),
     &                             NF_FOUT, 3, u2dgrd, vidUb)
#ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidUb,nf_collective)
#endif
#ifndef AVRH
          lvar=lenstr(vname(2,indxUb))
          ierr=nf_put_att_text (ncid, vidUb, 'long_name', lvar,
     &                                  vname(2,indxUb)(1:lvar))
#else
          text='averaged '/ /vname(2,indxUb)
          lvar=lenstr(text)
          ierr=nf_put_att_text (ncid, vidUb, 'long_name', lvar,
     &                                             text(1:lvar))
#endif
          lvar=lenstr(vname(3,indxUb))
          ierr=nf_put_att_text (ncid, vidUb, 'units',     lvar,
     &                                  vname(3,indxUb)(1:lvar))
          lvar=lenstr(vname(4,indxUb))
          ierr=nf_put_att_text (ncid, vidUb, 'field',    lvar,
     &                                  vname(4,indxUb)(1:lvar))

          call nf_add_attribute(ncid, vidUb, vname(:,indxUb), 
     &           5, NF_FOUT, ierr)
        endif

        if (wrt(indxVb)) then
          lvar=lenstr(vname(1,indxVb))
          ierr=nf_def_var (ncid, vname(1,indxVb)(1:lvar),
     &                              NF_FOUT, 3, v2dgrd, vidVb)
#ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidVb,nf_collective)
#endif
#ifndef AVRH
          lvar=lenstr(vname(2,indxVb))
          ierr=nf_put_att_text (ncid, vidVb, 'long_name', lvar,
     &                                  vname(2,indxVb)(1:lvar))
#else
          text='averaged '/ /vname(2,indxVb)
          lvar=lenstr(text)
          ierr=nf_put_att_text (ncid, vidVb, 'long_name', lvar,
     &                                             text(1:lvar))
#endif
          lvar=lenstr(vname(3,indxVb))
          ierr=nf_put_att_text (ncid, vidVb, 'units',     lvar,
     &                                  vname(3,indxVb)(1:lvar))
          lvar=lenstr(vname(4,indxVb))
          ierr=nf_put_att_text (ncid, vidVb, 'field',     lvar,
     &                                  vname(4,indxVb)(1:lvar))

          call nf_add_attribute(ncid, vidVb, vname(:,indxVb), 
     &        5, NF_FOUT, ierr)
        endif

#ifdef MORPHODYN
!
!-----------------------------------------------------------------------
!  Time-dependent bathymetry
!-----------------------------------------------------------------------
!
        if (wrt(indxHm)) then
          lvar=lenstr(vname(1,indxHm))
          ierr=nf_def_var (ncid, vname(1,indxHm)(1:lvar),
     &                              NF_FOUT, 3, r2dgrd, vidHm)
# ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidHm,nf_collective)
# endif
# ifndef AVRH
          lvar=lenstr(vname(2,indxHm))
          ierr=nf_put_att_text (ncid, vidHm, 'long_name', lvar,
     &                                  vname(2,indxHm)(1:lvar))
# else
          text='averaged '/ /vname(2,indxHm)
          lvar=lenstr(text)
          ierr=nf_put_att_text (ncid, vidHm, 'long_name', lvar,
     &                                             text(1:lvar))
# endif
          lvar=lenstr(vname(3,indxHm))
          ierr=nf_put_att_text (ncid, vidHm, 'units',     lvar,
     &                                  vname(3,indxHm)(1:lvar))
          lvar=lenstr(vname(4,indxHm))
          ierr=nf_put_att_text (ncid, vidHm, 'field',     lvar,
     &                                  vname(4,indxHm)(1:lvar))

          call nf_add_attribute(ncid, vidHm, vname(:,indxHm), 
     &         5, NF_FOUT, ierr)
        endif
#endif /* MORPHODYN */

#ifdef SOLVE3D
!
!-----------------------------------------------------------------------
! 3D momenta in XI- and ETA-directions.
!-----------------------------------------------------------------------
!
        if (wrt(indxU)) then
          lvar=lenstr(vname(1,indxU))
          ierr=nf_def_var (ncid, vname(1,indxU)(1:lvar),
     &                             NF_FOUT, 4, u3dgrd, vidU)
# ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidU,nf_collective)
# endif
# ifndef AVRH
          lvar=lenstr(vname(2,indxU))
          ierr=nf_put_att_text (ncid, vidU, 'long_name', lvar,
     &                                  vname(2,indxU)(1:lvar))
# else
          text='averaged '/ /vname(2,indxU)
          lvar=lenstr(text)
          ierr=nf_put_att_text (ncid, vidU, 'long_name', lvar,
     &                                            text(1:lvar))
# endif
          lvar=lenstr(vname(3,indxU))
          ierr=nf_put_att_text (ncid, vidU, 'units',     lvar,
     &                                  vname(3,indxU)(1:lvar))
          lvar=lenstr(vname(4,indxU))
          ierr=nf_put_att_text (ncid, vidU, 'field',     lvar,
     &                                  vname(4,indxU)(1:lvar))

          call nf_add_attribute(ncid, vidU, vname(:,indxU), 
     &         5, NF_FOUT, ierr)
        endif

        if (wrt(indxV)) then
          lvar=lenstr(vname(1,indxV))
          ierr=nf_def_var (ncid, vname(1,indxV)(1:lvar),
     &                             NF_FOUT, 4, v3dgrd, vidV)
# ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidV,nf_collective)
# endif
# ifndef AVRH
          lvar=lenstr(vname(2,indxV))
          ierr=nf_put_att_text (ncid, vidV, 'long_name', lvar,
     &                                  vname(2,indxV)(1:lvar))
# else
          text='averaged '/ /vname(2,indxV)
          lvar=lenstr(text)
          ierr=nf_put_att_text (ncid, vidV, 'long_name', lvar,
     &                                            text(1:lvar))
# endif
          lvar=lenstr(vname(3,indxV))
          ierr=nf_put_att_text (ncid, vidV, 'units',     lvar,
     &                                  vname(3,indxV)(1:lvar))
          lvar=lenstr(vname(4,indxV))
          ierr=nf_put_att_text (ncid, vidV, 'field',     lvar,
     &                                  vname(4,indxV)(1:lvar))

          call nf_add_attribute(ncid, vidV, vname(:,indxV), 
     &         5, NF_FOUT, ierr)
        endif
!
!-----------------------------------------------------------------------
! Tracer variables.
!-----------------------------------------------------------------------
!
# ifdef TRACERS
        do itrc=1,NT
          if (wrt(indxV+itrc)) then
            lvar=lenstr(vname(1,indxV+itrc))
            ierr=nf_def_var (ncid, vname(1,indxV+itrc)(1:lvar),
     &                             NF_FOUT, 4, r3dgrd, vidT(itrc))
#  ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidT(itrc),nf_collective)
#  endif
#  ifndef AVRH
            lvar=lenstr(vname(2,indxV+itrc))
            ierr=nf_put_att_text (ncid, vidT(itrc), 'long_name',
     &                         lvar, vname(2,indxV+itrc)(1:lvar))
#  else
            text='averaged '/ /vname(2,indxV+itrc)
            lvar=lenstr(text)
            ierr=nf_put_att_text (ncid, vidT(itrc), 'long_name',
     &                                          lvar, text(1:lvar))
#  endif
            lvar=lenstr(vname(3,indxV+itrc))
            ierr=nf_put_att_text (ncid, vidT(itrc), 'units', lvar,
     &                               vname(3,indxV+itrc)(1:lvar))
            lvar=lenstr(vname(4,indxV+itrc))
            ierr=nf_put_att_text (ncid, vidT(itrc), 'field', lvar,
     &                               vname(4,indxV+itrc)(1:lvar))

            call nf_add_attribute(ncid,vidT(itrc),vname(:,indxV+itrc),5,
     &           NF_FOUT, ierr)
          endif
        enddo
# endif
!
!-----------------------------------------------------------------------
! Density anomaly.
!-----------------------------------------------------------------------
!
        if (wrt(indxR)) then
          lvar=lenstr(vname(1,indxR))
          ierr=nf_def_var (ncid, vname(1,indxR)(1:lvar),
     &                             NF_FOUT, 4, r3dgrd, vidR)
# ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidR,nf_collective)
# endif
# ifndef AVRH
          lvar=lenstr(vname(2,indxR))
          ierr=nf_put_att_text (ncid, vidR, 'long_name', lvar,
     &                                  vname(2,indxR)(1:lvar))
# else
          text='averaged '/ /vname(2,indxR)
          lvar=lenstr(text)
          ierr=nf_put_att_text (ncid, vidR, 'long_name', lvar,
     &                                            text(1:lvar))
# endif
          lvar=lenstr(vname(3,indxR))
          ierr=nf_put_att_text (ncid, vidR, 'units',     lvar,
     &                                  vname(3,indxR)(1:lvar))
          lvar=lenstr(vname(4,indxR))
          ierr=nf_put_att_text (ncid, vidR, 'field',     lvar,
     &                                  vname(4,indxR)(1:lvar))

          call nf_add_attribute(ncid, vidR, vname(:,indxR), 
     &         5, NF_FOUT, ierr)
        endif
!
!-----------------------------------------------------------------------
! Brunt Vaisala Frequency.
!-----------------------------------------------------------------------
!
# if defined ANA_VMIX || defined BVF_MIXING \
  || defined LMD_MIXING || defined LMD_SKPP || defined LMD_BKPP \
  || defined GLS_MIXING
        if (wrt(indxbvf)) then
          lvar=lenstr(vname(1,indxbvf))
          ierr=nf_def_var (ncid, vname(1,indxbvf)(1:lvar),
     &                             NF_FOUT, 4, w3dgrd, vidbvf)
#  ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidbvf,nf_collective)
#  endif
#  ifndef AVRH
          lvar=lenstr(vname(2,indxbvf))
          ierr=nf_put_att_text (ncid, vidbvf, 'long_name', lvar,
     &                                  vname(2,indxbvf)(1:lvar))
#  else
          text='averaged '/ /vname(2,indxbvf)
          lvar=lenstr(text)
          ierr=nf_put_att_text (ncid, vidbvf, 'long_name', lvar,
     &                                            text(1:lvar))
#  endif
          lvar=lenstr(vname(3,indxbvf))
          ierr=nf_put_att_text (ncid, vidbvf, 'units',     lvar,
     &                                  vname(3,indxbvf)(1:lvar))
          lvar=lenstr(vname(4,indxbvf))
          ierr=nf_put_att_text (ncid, vidR, 'field',     lvar,
     &                                  vname(4,indxbvf)(1:lvar))

          call nf_add_attribute(ncid, vidbvf, vname(:,indxbvf), 
     &         5, NF_FOUT, ierr)
        endif
# endif
!
!-----------------------------------------------------------------------
! BIO_BioEBUS : AOU and pCO2
!-----------------------------------------------------------------------
!
# if defined BIOLOGY && defined BIO_BioEBUS
        if (wrt(indxAOU)) then
          lvar=lenstr(vname(1,indxAOU))
          ierr=nf_def_var (ncid, vname(1,indxAOU)(1:lvar),
     &                           NF_FOUT, 4, r3dgrd, vidAOU)
#  ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidAOU,nf_collective)
#  endif
#  ifndef AVRH
          lvar=lenstr(vname(2,indxAOU))
          ierr=nf_put_att_text (ncid, vidAOU, 'long_name',lvar,
     &                                vname(2,indxAOU)(1:lvar))
#  else
          text='averaged '/ /vname(2,indxAOU)
          lvar=lenstr(text)
          ierr=nf_put_att_text (ncid, vidAOU, 'long_name',lvar,
     &                                          text(1:lvar))
#  endif
          lvar=lenstr(vname(3,indxAOU))
          ierr=nf_put_att_text (ncid, vidAOU, 'units',lvar,
     &                                vname(3,indxAOU)(1:lvar))
          lvar=lenstr(vname(4,indxAOU))
          ierr=nf_put_att_text (ncid, vidAOU, 'field',lvar,
     &                                vname(4,indxAOU)(1:lvar))

          call nf_add_attribute(ncid, vidAOU, vname(:,indxAOU), 5,
     &    NF_FOUT, ierr)
        endif

        if (wrt(indxWIND10)) then
          lvar=lenstr(vname(1,indxWIND10))
          ierr=nf_def_var (ncid, vname(1,indxWIND10)(1:lvar),
     &                         NF_FOUT, 3, r2dgrd, vidwind10)
#  ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidwind10,nf_collective)
#  endif
#  ifndef AVRH
          lvar=lenstr(vname(2,indxWIND10))
          ierr=nf_put_att_text (ncid,vidwind10,'long_name',lvar,
     &                              vname(2,indxWIND10)(1:lvar))
#  else
          text='averaged '/ /vname(2,indxWIND10)
          lvar=lenstr(text)
          ierr=nf_put_att_text (ncid,vidwind10,'long_name',lvar,
     &                                        text(1:lvar))
#  endif
          lvar=lenstr(vname(3,indxWIND10))
          ierr=nf_put_att_text (ncid,vidwind10,'units',lvar,
     &                              vname(3,indxWIND10)(1:lvar))
          lvar=lenstr(vname(4,indxWIND10))
          ierr=nf_put_att_text (ncid,vidwind10,'field',lvar,
     &                              vname(4,indxWIND10)(1:lvar))

          call nf_add_attribute(ncid, vidwind10, vname(:,indxWIND10), 5,
     &    NF_FOUT, ierr)
        endif
# endif /* BIO_BioEBUS */
!
!-----------------------------------------------------------------------
! S-coordinate "omega" vertical velocity.
!-----------------------------------------------------------------------
!
        if (wrt(indxO)) then
          lvar=lenstr(vname(1,indxO))
          ierr=nf_def_var (ncid, vname(1,indxO)(1:lvar),
     &                             NF_FOUT, 4, w3dgrd, vidO)
# ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidO,nf_collective)
# endif
# ifndef AVRH
          lvar=lenstr(vname(2,indxO))
          ierr=nf_put_att_text (ncid, vidO, 'long_name', lvar,
     &                                  vname(2,indxO)(1:lvar))
# else
          text='averaged '/ /vname(2,indxO)
          lvar=lenstr(text)
          ierr=nf_put_att_text (ncid, vidO, 'long_name', lvar,
     &                                            text(1:lvar))
# endif
          lvar=lenstr(vname(3,indxO))
          ierr=nf_put_att_text (ncid, vidO, 'units',     lvar,
     &                                  vname(3,indxO)(1:lvar))
          lvar=lenstr(vname(4,indxO))
          ierr=nf_put_att_text (ncid, vidO, 'field',     lvar,
     &                                  vname(4,indxO)(1:lvar))

          call nf_add_attribute(ncid, vidO, vname(:,indxO), 
     &          5, NF_FOUT, ierr)
        endif
!
!-----------------------------------------------------------------------
! True vertical velocity.
!-----------------------------------------------------------------------
!
        if (wrt(indxW)) then
          lvar=lenstr(vname(1,indxW))
# ifdef NBQ
          ierr=nf_def_var (ncid, vname(1,indxW)(1:lvar),
     &                               NF_FOUT, 4, w3dgrd, vidW)
# else
          ierr=nf_def_var (ncid, vname(1,indxW)(1:lvar),
     &                               NF_FOUT, 4, r3dgrd, vidW)
# endif
# ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidW,nf_collective)
# endif
# ifndef AVRH
          lvar=lenstr(vname(2,indxW))
          ierr=nf_put_att_text (ncid, vidW, 'long_name', lvar,
     &                                  vname(2,indxW)(1:lvar))
# else
          text='averaged '/ /vname(2,indxW)
          lvar=lenstr(text)
          ierr=nf_put_att_text (ncid, vidW, 'long_name', lvar,
     &                                            text(1:lvar))
# endif
          lvar=lenstr(vname(3,indxW))
          ierr=nf_put_att_text (ncid, vidW, 'units',     lvar,
     &                                  vname(3,indxW)(1:lvar))
          lvar=lenstr(vname(4,indxW))
          ierr=nf_put_att_text (ncid, vidW, 'field',     lvar,
     &                                  vname(4,indxW)(1:lvar))

          call nf_add_attribute(ncid, vidW, vname(:,indxW), 
     &        5, NF_FOUT, ierr)
        endif
        
#ifdef M3FAST_HIS

!****** Rho *****
!
 
        if (wrt(indxRnbq)) then
          lvar=lenstr(vname(1,indxRnbq))
          ierr=nf_def_var (ncid, vname(1,indxRnbq)(1:lvar),
     &                               NF_FOUT, 4, r3dgrd, vidRnbq)
# ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidRnbq,nf_collective)
# endif
# ifndef AVRH
          lvar=lenstr(vname(2,indxRnbq))
          ierr=nf_put_att_text (ncid, vidRnbq, 'long_name', lvar,
     &                                  vname(2,indxRnbq)(1:lvar))
# else
          text='averaged '/ /vname(2,indxRnbq)
          lvar=lenstr(text)
          ierr=nf_put_att_text (ncid, vidRnbq, 'long_name', lvar,
     &                                            text(1:lvar))
# endif
          lvar=lenstr(vname(3,indxRnbq))
          ierr=nf_put_att_text (ncid, vidRnbq, 'units',     lvar,
     &                                  vname(3,indxRnbq)(1:lvar))
          lvar=lenstr(vname(4,indxRnbq))
          ierr=nf_put_att_text (ncid, vidRnbq, 'field',     lvar,
     &                                  vname(4,indxRnbq)(1:lvar))

          call nf_add_attribute(ncid, vidRnbq,
     &         indxRnbq, 5, NF_FOUT, ierr)

        endif
!
!****** U *****
!
        if (wrt(indxUnbq)) then
          
          lvar=lenstr(vname(1,indxUnbq))
          ierr=nf_def_var (ncid, vname(1,indxUnbq)(1:lvar),
     &                               NF_FOUT, 4, u3dgrd, vidUnbq)
# ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidUnbq,nf_collective)
# endif
# ifndef AVRH
          lvar=lenstr(vname(2,indxUnbq))
          ierr=nf_put_att_text (ncid, vidUnbq, 'long_name', lvar,
     &                                  vname(2,indxUnbq)(1:lvar))
# else
          text='averaged '/ /vname(2,indxUnbq)
          lvar=lenstr(text)
          ierr=nf_put_att_text (ncid, vidUnbq, 'long_name', lvar,
     &                                            text(1:lvar))
# endif
          lvar=lenstr(vname(3,indxUnbq))
          ierr=nf_put_att_text (ncid, vidUnbq, 'units',     lvar,
     &                                  vname(3,indxUnbq)(1:lvar))
          lvar=lenstr(vname(4,indxUnbq))
          ierr=nf_put_att_text (ncid, vidUnbq, 'field',     lvar,
     &                                  vname(4,indxUnbq)(1:lvar))

          call nf_add_attribute(ncid, vidUnbq, 
     &        indxUnbq, 5, NF_FOUT, ierr)

        endif
!
!****** V *****
!
        if (wrt(indxVnbq)) then
          lvar=lenstr(vname(1,indxVnbq))
          ierr=nf_def_var (ncid, vname(1,indxVnbq)(1:lvar),
     &                               NF_FOUT, 4, v3dgrd, vidVnbq)
# ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidVnbq,nf_collective)
# endif
# ifndef AVRH
          lvar=lenstr(vname(2,indxVnbq))
          ierr=nf_put_att_text (ncid, vidVnbq, 'long_name', lvar,
     &                                  vname(2,indxVnbq)(1:lvar))
# else
          text='averaged '/ /vname(2,indxVnbq)
          lvar=lenstr(text)
          ierr=nf_put_att_text (ncid, vidVnbq, 'long_name', lvar,
     &                                            text(1:lvar))
# endif
          lvar=lenstr(vname(3,indxVnbq))
          ierr=nf_put_att_text (ncid, vidVnbq, 'units',     lvar,
     &                                  vname(3,indxVnbq)(1:lvar))
          lvar=lenstr(vname(4,indxVnbq))
          ierr=nf_put_att_text (ncid, vidVnbq, 'field',     lvar,
     &                                  vname(4,indxVnbq)(1:lvar))

          call nf_add_attribute(ncid, vidVnbq, 
     &         indxVnbq, 5, NF_FOUT, ierr)

        endif
!
!****** W *****
!
        if (wrt(indxWnbq)) then
          lvar=lenstr(vname(1,indxWnbq))
          ierr=nf_def_var (ncid, vname(1,indxWnbq)(1:lvar),
     &                               NF_FOUT, 4, w3dgrd, vidWnbq)
# ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidWnbq,nf_collective)
# endif
# ifndef AVRH
          lvar=lenstr(vname(2,indxWnbq))
          ierr=nf_put_att_text (ncid, vidWnbq, 'long_name', lvar,
     &                                  vname(2,indxWnbq)(1:lvar))
# else
          text='averaged '/ /vname(2,indxWnbq)
          lvar=lenstr(text)
          ierr=nf_put_att_text (ncid, vidWnbq, 'long_name', lvar,
     &                                            text(1:lvar))
# endif
          lvar=lenstr(vname(3,indxWnbq))
          ierr=nf_put_att_text (ncid, vidWnbq, 'units',     lvar,
     &                                  vname(3,indxWnbq)(1:lvar))
          lvar=lenstr(vname(4,indxWnbq))
          ierr=nf_put_att_text (ncid, vidWnbq, 'field',     lvar,
     &                                  vname(4,indxWnbq)(1:lvar))

          call nf_add_attribute(ncid, vidWnbq, 
     &          indxWnbq, 5, NF_FOUT, ierr)

        endif
!
!****** Csound *****
!
        if (wrt(indxCnbq)) then
          lvar=lenstr(vname(1,indxCnbq))
          ierr=nf_def_var (ncid, vname(1,indxCnbq)(1:lvar),
     &                               NF_FOUT, 4, r3dgrd, vidCnbq)
# ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidCnbq,nf_collective)
# endif
# ifndef AVRH
          lvar=lenstr(vname(2,indxCnbq))
          ierr=nf_put_att_text (ncid, vidCnbq, 'long_name', lvar,
     &                                  vname(2,indxCnbq)(1:lvar))
# else
          text='averaged '/ /vname(2,indxCnbq)
          lvar=lenstr(text)
          ierr=nf_put_att_text (ncid, vidCnbq, 'long_name', lvar,
     &                                            text(1:lvar))
# endif
          lvar=lenstr(vname(3,indxCnbq))
          ierr=nf_put_att_text (ncid, vidCnbq, 'units',     lvar,
     &                                  vname(3,indxCnbq)(1:lvar))
          lvar=lenstr(vname(4,indxCnbq))
          ierr=nf_put_att_text (ncid, vidCnbq, 'field',     lvar,
     &                                  vname(4,indxCnbq)(1:lvar))

          call nf_add_attribute(ncid, vidCnbq,
     &         indxCnbq, 5, NF_FOUT, ierr)

        endif
#endif

#endif /* SOLVE3D */

!================================================================
!
! Define secondary variables
!
!================================================================
        if (wrt(indxBostr)) then
          lvar=lenstr(vname(1,indxBostr))
          ierr=nf_def_var (ncid, vname(1,indxBostr)(1:lvar),
     &                             NF_FOUT, 3, r2dgrd, vidBostr)
#ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidBostr,nf_collective)
#endif
#ifndef AVRH
          lvar=lenstr(vname(2,indxBostr))
          ierr=nf_put_att_text (ncid, vidBostr, 'long_name', lvar,
     &                                 vname(2,indxBostr)(1:lvar))
#else
          text='averaged '/ /vname(2,indxBostr)
          lvar=lenstr(text)
          ierr=nf_put_att_text (ncid, vidBostr, 'long_name', lvar,
     &                                               text(1:lvar))
#endif
          lvar=lenstr(vname(3,indxBostr))
          ierr=nf_put_att_text (ncid, vidBostr, 'units',     lvar,
     &                                 vname(3,indxBostr)(1:lvar))

          call nf_add_attribute(ncid, vidBostr,vname(:,indxBostr),5,
     &                          NF_FOUT, ierr)
        endif



        if (wrt(indxBustr)) then
          lvar=lenstr(vname(1,indxBustr))
          ierr=nf_def_var (ncid, vname(1,indxBustr)(1:lvar),
     &                             NF_FOUT, 3, u2dgrd, vidBustr)
#ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidBustr,nf_collective)
#endif
#ifndef AVRH
          lvar=lenstr(vname(2,indxBustr))
          ierr=nf_put_att_text (ncid, vidBustr, 'long_name', lvar,
     &                                 vname(2,indxBustr)(1:lvar))
#else
          text='averaged '/ /vname(2,indxBustr)
          lvar=lenstr(text)
          ierr=nf_put_att_text (ncid, vidBustr, 'long_name', lvar,
     &                                               text(1:lvar))
#endif
          lvar=lenstr(vname(3,indxBustr))
          ierr=nf_put_att_text (ncid, vidBustr, 'units',     lvar,
     &                                 vname(3,indxBustr)(1:lvar))

          call nf_add_attribute(ncid, vidBustr,vname(:,indxBustr),5,
     &                          NF_FOUT, ierr)

        endif


        if (wrt(indxBvstr)) then
          lvar=lenstr(vname(1,indxBvstr))
          ierr=nf_def_var (ncid, vname(1,indxBvstr)(1:lvar),
     &                             NF_FOUT, 3, v2dgrd, vidBvstr)
#ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidBvstr,nf_collective)
#endif
#ifndef AVRH
          lvar=lenstr(vname(2,indxBvstr))
          ierr=nf_put_att_text (ncid, vidBvstr, 'long_name', lvar,
     &                                 vname(2,indxBvstr)(1:lvar))
#else
          text='averaged '/ /vname(2,indxBvstr)
          lvar=lenstr(text)
          ierr=nf_put_att_text (ncid, vidBvstr, 'long_name', lvar,
     &                                               text(1:lvar))
#endif
          lvar=lenstr(vname(3,indxBvstr))
          ierr=nf_put_att_text (ncid, vidBvstr, 'units',     lvar,
     &                                 vname(3,indxBvstr)(1:lvar))

          call nf_add_attribute(ncid, vidBvstr,vname(:,indxBvstr),5,
     &                          NF_FOUT, ierr)

        endif
!
        if (wrt(indxWstr)) then
          lvar=lenstr(vname(1,indxWstr))
          ierr=nf_def_var (ncid, vname(1,indxWstr)(1:lvar),
     &                             NF_FOUT, 3, r2dgrd, vidWstr)
#ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidWstr,nf_collective)
#endif
#ifndef AVRH
          lvar=lenstr(vname(2,indxWstr))
          ierr=nf_put_att_text (ncid, vidWstr, 'long_name', lvar,
     &                                 vname(2,indxWstr)(1:lvar))
#else
          text='averaged '/ /vname(2,indxWstr)
          lvar=lenstr(text)
          ierr=nf_put_att_text (ncid, vidWstr, 'long_name', lvar,
     &                                               text(1:lvar))
#endif
          lvar=lenstr(vname(3,indxWstr))
          ierr=nf_put_att_text (ncid, vidWstr, 'units',     lvar,
     &                                 vname(3,indxWstr)(1:lvar))

          call nf_add_attribute(ncid, vidWstr, vname(:,indxWstr), 5,
     &                          NF_FOUT, ierr)
        endif
!
        if (wrt(indxUWstr)) then
          lvar=lenstr(vname(1,indxUWstr))
#ifdef STRESS_AT_RHO_POINTS
          ierr=nf_def_var (ncid, vname(1,indxUWstr)(1:lvar),
     &                             NF_FOUT, 3, r2dgrd, vidUWstr)
#else
          ierr=nf_def_var (ncid, vname(1,indxUWstr)(1:lvar),
     &                             NF_FOUT, 3, u2dgrd, vidUWstr)
#endif
#ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidUWstr,nf_collective)
#endif
#ifndef AVRH
          lvar=lenstr(vname(2,indxUWstr))
          ierr=nf_put_att_text (ncid, vidUWstr, 'long_name', lvar,
     &                                 vname(2,indxUWstr)(1:lvar))
#else
          text='averaged '/ /vname(2,indxUWstr)
          lvar=lenstr(text)
          ierr=nf_put_att_text (ncid, vidUWstr, 'long_name', lvar,
     &                                               text(1:lvar))
#endif
          lvar=lenstr(vname(3,indxUWstr))
          ierr=nf_put_att_text (ncid, vidUWstr, 'units',     lvar,
     &                                 vname(3,indxUWstr)(1:lvar))

          call nf_add_attribute(ncid, vidUWstr, vname(:,indxUWstr), 5,
     &                          NF_FOUT, ierr)
        endif
!
        if (wrt(indxVWstr)) then
          lvar=lenstr(vname(1,indxVWstr))
#ifdef STRESS_AT_RHO_POINTS
          ierr=nf_def_var (ncid, vname(1,indxVWstr)(1:lvar),
     &                             NF_FOUT, 3, r2dgrd, vidVWstr)
#else
          ierr=nf_def_var (ncid, vname(1,indxVWstr)(1:lvar),
     &                             NF_FOUT, 3, v2dgrd, vidVWstr)
#endif
#ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidVWstr,nf_collective)
#endif
#ifndef AVRH
          lvar=lenstr(vname(2,indxVWstr))
          ierr=nf_put_att_text (ncid, vidVWstr, 'long_name', lvar,
     &                                 vname(2,indxVWstr)(1:lvar))
#else
          text='averaged '/ /vname(2,indxVWstr)
          lvar=lenstr(text)
          ierr=nf_put_att_text (ncid, vidVWstr, 'long_name', lvar,
     &                                               text(1:lvar))
#endif
          lvar=lenstr(vname(3,indxVWstr))
          ierr=nf_put_att_text (ncid, vidVWstr, 'units',     lvar,
     &                                 vname(3,indxVWstr)(1:lvar))

          call nf_add_attribute(ncid, vidVWstr, vname(:,indxVWstr), 5,
     &                          NF_FOUT, ierr)
        endif

#ifdef SOLVE3D
# ifdef VIS_COEF_3D
!
! Horizontal viscosity coefficient.
!
        if (wrt(indxVisc)) then
          lvar=lenstr(vname(1,indxVisc))
          ierr=nf_def_var (ncid, vname(1,indxVisc)(1:lvar),
     &                             NF_FOUT, 4, r3dgrd, vidVisc)
#  ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidVisc,nf_collective)
#  endif
#  ifndef AVRH
          lvar=lenstr(vname(2,indxVisc))
          ierr=nf_put_att_text (ncid, vidVisc, 'long_name', lvar,
     &                                  vname(2,indxVisc)(1:lvar))
#  else
          text='averaged '/ /vname(2,indxVisc)
          lvar=lenstr(text)
          ierr=nf_put_att_text (ncid, vidVisc, 'long_name', lvar,
     &                                              text(1:lvar))
#  endif
          lvar=lenstr(vname(3,indxVisc))
          ierr=nf_put_att_text (ncid, vidVisc, 'units',     lvar,
     &                                  vname(3,indxVisc)(1:lvar))
          lvar=lenstr(vname(4,indxVisc))
          ierr=nf_put_att_text (ncid, vidVisc, 'field',     lvar,
     &                                  vname(4,indxVisc)(1:lvar))

          call nf_add_attribute(ncid, vidVisc, vname(:,indxVisc), 5,
     &                          NF_FOUT, ierr)
        endif
# endif
# ifdef DIF_COEF_3D
!
! Horizontal diffusivity coefficient.
!
        if (wrt(indxDiff)) then
          lvar=lenstr(vname(1,indxDiff))
          ierr=nf_def_var (ncid, vname(1,indxDiff)(1:lvar),
     &                             NF_FOUT, 4, r3dgrd, vidDiff)
#  ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidDiff,nf_collective)
#  endif
#  ifndef AVRH
          lvar=lenstr(vname(2,indxDiff))
          ierr=nf_put_att_text (ncid, vidDiff, 'long_name', lvar,
     &                                  vname(2,indxDiff)(1:lvar))
#  else
          text='averaged '/ /vname(2,indxDiff)
          lvar=lenstr(text)
          ierr=nf_put_att_text (ncid, vidDiff, 'long_name', lvar,
     &                                              text(1:lvar))
#  endif
          lvar=lenstr(vname(4,indxDiff))
          ierr=nf_put_att_text (ncid, vidDiff, 'field',     lvar,
     &                                  vname(4,indxDiff)(1:lvar))

          call nf_add_attribute(ncid, vidDiff, vname(:,indxDiff), 5,
     &                          NF_FOUT, ierr)
        endif
# endif
!
!     Vertical viscosity coefficient.
!
        if (wrt(indxAkv)) then
          lvar=lenstr(vname(1,indxAkv))
          ierr=nf_def_var (ncid, vname(1,indxAkv)(1:lvar),
     &                             NF_FOUT, 4, w3dgrd, vidAkv)
# ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidAkv,nf_collective)
# endif
# ifndef AVRH
          lvar=lenstr(vname(2,indxAkv))
          ierr=nf_put_att_text (ncid, vidAkv, 'long_name', lvar,
     &                                  vname(2,indxAkv)(1:lvar))
# else
          text='averaged '/ /vname(2,indxAkv)
          lvar=lenstr(text)
          ierr=nf_put_att_text (ncid, vidAkv, 'long_name', lvar,
     &                                              text(1:lvar))
# endif
          lvar=lenstr(vname(3,indxAkv))
          ierr=nf_put_att_text (ncid, vidAkv, 'units',     lvar,
     &                                  vname(3,indxAkv)(1:lvar))
          lvar=lenstr(vname(4,indxAkv))
          ierr=nf_put_att_text (ncid, vidAkv, 'field',     lvar,
     &                                  vname(4,indxAkv)(1:lvar))

          call nf_add_attribute(ncid, vidAkv, vname(:,indxAkv), 
     &                          5, NF_FOUT, ierr)
        endif
!
! Vertical diffusion coefficient for potential temperature.
!
# ifdef TEMPERATURE
        if (wrt(indxAkt)) then
          lvar=lenstr(vname(1,indxAkt))
          ierr=nf_def_var (ncid, vname(1,indxAkt)(1:lvar),
     &                              NF_FOUT, 4, w3dgrd, vidAkt)
#  ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidAkt,nf_collective)
#  endif
#  ifndef AVRH
          lvar=lenstr(vname(2,indxAkt))
          ierr=nf_put_att_text (ncid, vidAkt, 'long_name', lvar,
     &                                  vname(2,indxAkt)(1:lvar))
#  else
          text='averaged '/ /vname(2,indxAkt)
          lvar=lenstr(text)
          ierr=nf_put_att_text (ncid, vidAkt, 'long_name', lvar,
     &                                              text(1:lvar))
#  endif
          lvar=lenstr(vname(3,indxAkt))
          ierr=nf_put_att_text (ncid, vidAkt, 'units',     lvar,
     &                                  vname(3,indxAkt)(1:lvar))
          lvar=lenstr(vname(4,indxAkt))
          ierr=nf_put_att_text (ncid, vidAkt, 'field',     lvar,
     &                                  vname(4,indxAkt)(1:lvar))

          call nf_add_attribute(ncid, vidAkt, vname(:,indxAkt), 
     &                           5, NF_FOUT, ierr)
        endif
# endif
# ifdef SALINITY
!
! Vertical diffusion coefficient for salinity.
!
        if (wrt(indxAks)) then
          lvar=lenstr(vname(1,indxAks))
          ierr=nf_def_var (ncid, vname(1,indxAks)(1:lvar),
     &                              NF_FOUT, 4, w3dgrd, vidAks)
#  ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidAks,nf_collective)
#  endif
#  ifndef AVRH
          lvar=lenstr(vname(2,indxAks))
          ierr=nf_put_att_text (ncid, vidAks, 'long_name', lvar,
     &                                  vname(2,indxAks)(1:lvar))
#  else
          text='averaged '/ /vname(2,indxAks)
          lvar=lenstr(text)
          ierr=nf_put_att_text (ncid, vidAks, 'long_name', lvar,
     &                                              text(1:lvar))
#  endif
          lvar=lenstr(vname(3,indxAks))
          ierr=nf_put_att_text (ncid, vidAks, 'units',     lvar,
     &                                  vname(3,indxAks)(1:lvar))
          lvar=lenstr(vname(4,indxAks))
          ierr=nf_put_att_text (ncid, vidAks, 'field',     lvar,
     &                                  vname(4,indxAks)(1:lvar))

          call nf_add_attribute(ncid, vidAks, vname(:,indxAks), 
     &                          5, NF_FOUT, ierr)
        endif
# endif


# if defined LMD_SKPP || defined GLS_MIXING
!
! Depth of planetary boundary layer.
!
        if (wrt(indxHbl)) then
          lvar=lenstr(vname(1,indxHbl))
          ierr=nf_def_var (ncid, vname(1,indxHbl)(1:lvar),
     &                             NF_FOUT, 3, r2dgrd, vidHbl)
#  ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidHbl,nf_collective)
#  endif
#  ifndef AVRH
          lvar=lenstr(vname(2,indxHbl))
          ierr=nf_put_att_text (ncid, vidHbl, 'long_name', lvar,
     &                                  vname(2,indxHbl)(1:lvar))
#  else
          text='averaged '/ /vname(2,indxHbl)
          lvar=lenstr(text)
          ierr=nf_put_att_text (ncid, vidHbl, 'long_name', lvar,
     &                                              text(1:lvar))
#  endif
          lvar=lenstr(vname(3,indxHbl))
          ierr=nf_put_att_text (ncid, vidHbl, 'units',     lvar,
     &                                  vname(3,indxHbl)(1:lvar))
          lvar=lenstr(vname(4,indxHbl))
          ierr=nf_put_att_text (ncid, vidHbl, 'field',     lvar,
     &                                  vname(4,indxHbl)(1:lvar))

          call nf_add_attribute(ncid, vidHbl, vname(:,indxHbl), 
     &                          5, NF_FOUT, ierr)
        endif
# endif
# ifdef LMD_BKPP
!
! Depth of bottom planetary boundary layer.
!
        if (wrt(indxHbbl)) then
          lvar=lenstr(vname(1,indxHbbl))
          ierr=nf_def_var (ncid, vname(1,indxHbbl)(1:lvar),
     &                             NF_FOUT, 3, r2dgrd, vidHbbl)
#  ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidHbbl,nf_collective)
#  endif
#  ifndef AVRH
          lvar=lenstr(vname(2,indxHbbl))
          ierr=nf_put_att_text (ncid, vidHbbl, 'long_name', lvar,
     &                                  vname(2,indxHbbl)(1:lvar))
#  else
          text='averaged '/ /vname(2,indxHbbl)
          lvar=lenstr(text)
          ierr=nf_put_att_text (ncid, vidHbbl, 'long_name', lvar,
     &                                              text(1:lvar))
#  endif
          lvar=lenstr(vname(3,indxHbbl))
          ierr=nf_put_att_text (ncid, vidHbbl, 'units',     lvar,
     &                                  vname(3,indxHbbl)(1:lvar))
          lvar=lenstr(vname(4,indxHbbl))
          ierr=nf_put_att_text (ncid, vidHbbl, 'field',     lvar,
     &                                  vname(4,indxHbbl)(1:lvar))

          call nf_add_attribute(ncid, vidHbbl, vname(:,indxHbbl), 5,
     &                          NF_FOUT, ierr)
        endif
# endif

# ifdef GLS_MIXING
!
! Turbulent Kinetic Energy.
!
        if (wrt(indxTke)) then
          lvar=lenstr(vname(1,indxTke))
          ierr=nf_def_var (ncid, vname(1,indxTke)(1:lvar),
     &                              NF_FOUT, 4, w3dgrd, vidTke)
#  ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidTke,nf_collective)
#  endif
#  ifndef AVRH
          lvar=lenstr(vname(2,indxTke))
          ierr=nf_put_att_text (ncid, vidTke, 'long_name', lvar,
     &                                  vname(2,indxTke)(1:lvar))
#  else
          text='averaged '/ /vname(2,indxTke)
          lvar=lenstr(text)
          ierr=nf_put_att_text (ncid, vidTke, 'long_name', lvar,
     &                                              text(1:lvar))
#  endif
          lvar=lenstr(vname(3,indxTke))
          ierr=nf_put_att_text (ncid, vidTke, 'units',     lvar,
     &                                  vname(3,indxTke)(1:lvar))
          lvar=lenstr(vname(4,indxTke))
          ierr=nf_put_att_text (ncid, vidTke, 'field',     lvar,
     &                                  vname(4,indxTke)(1:lvar))

          call nf_add_attribute(ncid, vidTke, vname(:,indxTke), 
     &                          5, NF_FOUT, ierr)
        endif
!
! Generic Length Scale.
!
        if (wrt(indxGls)) then
          lvar=lenstr(vname(1,indxGls))
          ierr=nf_def_var (ncid, vname(1,indxGls)(1:lvar),
     &                              NF_FOUT, 4, w3dgrd, vidGls)
#  ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidGls,nf_collective)
#  endif
#  ifndef AVRH
          lvar=lenstr(vname(2,indxGls))
          ierr=nf_put_att_text (ncid, vidGls, 'long_name', lvar,
     &                                  vname(2,indxGls)(1:lvar))
#  else
          text='averaged '/ /vname(2,indxGls)
          lvar=lenstr(text)
          ierr=nf_put_att_text (ncid, vidGls, 'long_name', lvar,
     &                                              text(1:lvar))
#  endif
          lvar=lenstr(vname(3,indxGls))
          ierr=nf_put_att_text (ncid, vidGls, 'units',     lvar,
     &                                  vname(3,indxGls)(1:lvar))
          lvar=lenstr(vname(4,indxGls))
          ierr=nf_put_att_text (ncid, vidGls, 'field',     lvar,
     &                                  vname(4,indxGls)(1:lvar))

          call nf_add_attribute(ncid, vidGls, vname(:,indxGls), 
     &                          5, NF_FOUT, ierr)
        endif
!
! Vertical mixing turbulent length scale.
!
        if (wrt(indxLsc)) then
          lvar=lenstr(vname(1,indxLsc))
          ierr=nf_def_var (ncid, vname(1,indxLsc)(1:lvar),
     &                              NF_FOUT, 4, w3dgrd, vidLsc)
#  ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidLsc,nf_collective)
#  endif
#  ifndef AVRH
          lvar=lenstr(vname(2,indxLsc))
          ierr=nf_put_att_text (ncid, vidLsc, 'long_name', lvar,
     &                                  vname(2,indxLsc)(1:lvar))
#  else
          text='averaged '/ /vname(2,indxLsc)
          lvar=lenstr(text)
          ierr=nf_put_att_text (ncid, vidLsc, 'long_name', lvar,
     &                                              text(1:lvar))
#  endif
          lvar=lenstr(vname(3,indxLsc))
          ierr=nf_put_att_text (ncid, vidLsc, 'units',     lvar,
     &                                  vname(3,indxLsc)(1:lvar))
          lvar=lenstr(vname(4,indxLsc))
          ierr=nf_put_att_text (ncid, vidLsc, 'field',     lvar,
     &                                  vname(4,indxLsc)(1:lvar))

          call nf_add_attribute(ncid, vidLsc, vname(:,indxLsc),
     &                          5, NF_FOUT, ierr)
        endif
# endif /* GLS_MIXING */

!
! Shflx [W/m2]
!
# ifdef TEMPERATURE
        if (wrt(indxShflx)) then
          lvar=lenstr(vname(1,indxShflx))
          ierr=nf_def_var (ncid, vname(1,indxShflx)(1:lvar),
     &                             NF_FOUT, 3, r2dgrd, vidShflx)
#  ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidShflx,nf_collective)
#  endif
#  ifndef AVRH
          lvar=lenstr(vname(2,indxShflx))
          ierr=nf_put_att_text (ncid, vidShflx, 'long_name', lvar,
     &                                 vname(2,indxShflx)(1:lvar))
#  else
          text='averaged '/ /vname(2,indxShflx)
          lvar=lenstr(text)
          ierr=nf_put_att_text (ncid, vidShflx, 'long_name', lvar,
     &                                               text(1:lvar))
#  endif
          lvar=lenstr(vname(3,indxShflx))
          ierr=nf_put_att_text (ncid, vidShflx, 'units',     lvar,
     &                                 vname(3,indxShflx)(1:lvar))

          call nf_add_attribute(ncid, vidShflx, vname(:,indxShflx), 5,
     &                          NF_FOUT, ierr)
        endif
# endif
!
! Swflx [cm/day]
!
# ifdef SALINITY
        if (wrt(indxSwflx)) then
          lvar=lenstr(vname(1,indxSwflx))
          ierr=nf_def_var (ncid, vname(1,indxSwflx)(1:lvar),
     &                             NF_FOUT, 3, r2dgrd, vidSwflx)
#  ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidSwflx,nf_collective)
#  endif
#  ifndef AVRH
          lvar=lenstr(vname(2,indxSwflx))
          ierr=nf_put_att_text (ncid, vidSwflx, 'long_name', lvar,
     &                                 vname(2,indxSwflx)(1:lvar))
#  else
          text='averaged '/ /vname(2,indxSwflx)
          lvar=lenstr(text)
          ierr=nf_put_att_text (ncid, vidSwflx, 'long_name', lvar,
     &                                               text(1:lvar))
#  endif
          lvar=lenstr(vname(3,indxSwflx))
          ierr=nf_put_att_text (ncid, vidSwflx, 'units',     lvar,
     &                                 vname(3,indxSwflx)(1:lvar))

          call nf_add_attribute(ncid, vidSwflx, vname(:,indxSwflx), 5,
     &                          NF_FOUT, ierr)
        endif
# endif

# if defined BHFLUX && defined TEMPERATURE
!
! Bhflx [W/m2]
!
        if (wrt(indxBhflx)) then
          lvar=lenstr(vname(1,indxBhflx))
          ierr=nf_def_var (ncid, vname(1,indxBhflx)(1:lvar),
     &                             NF_FOUT, 3, r2dgrd, vidBhflx)
#  ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidBhflx,nf_collective)
#  endif
#  ifndef AVRH
          lvar=lenstr(vname(2,indxBhflx))
          ierr=nf_put_att_text (ncid, vidBhflx, 'long_name', lvar,
     &                                 vname(2,indxBhflx)(1:lvar))
#  else
          text='averaged '/ /vname(2,indxBhflx)
          lvar=lenstr(text)
          ierr=nf_put_att_text (ncid, vidBhflx, 'long_name', lvar,
     &                                               text(1:lvar))
#  endif
          lvar=lenstr(vname(3,indxBhflx))
          ierr=nf_put_att_text (ncid, vidBhflx, 'units',     lvar,
     &                                 vname(3,indxBhflx)(1:lvar))

          call nf_add_attribute(ncid, vidBhflx, vname(:,indxBhflx), 5,
     &                          NF_FOUT, ierr)

        endif
# endif

# if defined BWFLUX && defined SALINITY
!
! Bwflx [cm/day]
!
        if (wrt(indxBwflx)) then
          lvar=lenstr(vname(1,indxBwflx))
          ierr=nf_def_var (ncid, vname(1,indxBwflx)(1:lvar),
     &                             NF_FOUT, 3, r2dgrd, vidBwflx)
#  ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidBwflx,nf_collective)
#  endif
#  ifndef AVRH
          lvar=lenstr(vname(2,indxBwflx))
          ierr=nf_put_att_text (ncid, vidBwflx, 'long_name', lvar,
     &                                 vname(2,indxBwflx)(1:lvar))
#  else
          text='averaged '/ /vname(2,indxBwflx)
          lvar=lenstr(text)
          ierr=nf_put_att_text (ncid, vidBwflx, 'long_name', lvar,
     &                                               text(1:lvar))
#  endif
          lvar=lenstr(vname(3,indxBwflx))
          ierr=nf_put_att_text (ncid, vidBwflx, 'units',     lvar,
     &                                 vname(3,indxBwflx)(1:lvar))

          call nf_add_attribute(ncid, vidBwflx, vname(:,indxBwflx), 5,
     &                          NF_FOUT, ierr)

        endif
# endif
!
! Short wave Shflx_rsw [W/m2]
!
# ifdef TEMPERATURE
      if (wrt(indxShflx_rsw)) then
        lvar=lenstr(vname(1,indxShflx_rsw))
        ierr=nf_def_var (ncid, vname(1,indxShflx_rsw)(1:lvar),
     &                           NF_FOUT, 3, r2dgrd, vidShflx_rsw)
#  ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidShflx_rsw,nf_collective)
#  endif
#  ifndef AVRH
          lvar=lenstr(vname(2,indxShflx_rsw))
          ierr=nf_put_att_text (ncid, vidShflx_rsw, 'long_name', lvar,
     &                                 vname(2,indxShflx_rsw)(1:lvar))
#  else
          text='averaged '/ /vname(2,indxShflx_rsw)
          lvar=lenstr(text)
          ierr=nf_put_att_text (ncid, vidShflx_rsw, 'long_name', lvar,
     &                                               text(1:lvar))
#  endif
          lvar=lenstr(vname(3,indxShflx_rsw))
          ierr=nf_put_att_text (ncid, vidShflx_rsw, 'units',     lvar,
     &                                 vname(3,indxShflx_rsw)(1:lvar))

        call nf_add_attribute(ncid, vidShflx_rsw, 
     &               vname(:,indxShflx_rsw), 5, NF_FOUT, ierr)

      endif
# endif


# if defined BULK_FLUX && defined TEMPERATURE
!
! Long wave Shflx [W/m2]
!
        if (wrt(indxShflx_rlw)) then
          lvar=lenstr(vname(1,indxShflx_rlw))
          ierr=nf_def_var (ncid, vname(1,indxShflx_rlw)(1:lvar),
     &                           NF_FOUT, 3, r2dgrd, vidShflx_rlw)
#  ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidShflx_rlw,nf_collective)
#  endif
#  ifndef AVRH
          lvar=lenstr(vname(2,indxShflx_rlw))
          ierr=nf_put_att_text (ncid, vidShflx_rlw, 'long_name', lvar,
     &                               vname(2,indxShflx_rlw)(1:lvar))
#  else
          text='averaged '/ /vname(2,indxShflx_rlw)
          lvar=lenstr(text)
          ierr=nf_put_att_text (ncid, vidShflx_rlw, 'long_name', lvar,
     &                                             text(1:lvar))
#  endif
          lvar=lenstr(vname(3,indxShflx_rlw))
          ierr=nf_put_att_text (ncid, vidShflx_rlw, 'units',     lvar,
     &                               vname(3,indxShflx_rlw)(1:lvar))

          call nf_add_attribute(ncid, vidShflx_rlw, 
     &         vname(:,indxShflx_rlw),5, NF_FOUT, ierr)
        endif
!
! Latent Shflx [W/m2]
!
      if (wrt(indxShflx_lat)) then
        lvar=lenstr(vname(1,indxShflx_lat))
        ierr=nf_def_var (ncid, vname(1,indxShflx_lat)(1:lvar),
     &                           NF_FOUT, 3, r2dgrd, vidShflx_lat)
#  ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidShflx_lat,nf_collective)
#  endif
#  ifndef AVRH
          lvar=lenstr(vname(2,indxShflx_lat))
          ierr=nf_put_att_text (ncid, vidShflx_lat, 'long_name', lvar,
     &                               vname(2,indxShflx_lat)(1:lvar))
#  else
          text='averaged '/ /vname(2,indxShflx_lat)
          lvar=lenstr(text)
          ierr=nf_put_att_text (ncid, vidShflx_lat, 'long_name', lvar,
     &                                               text(1:lvar))
#  endif
          lvar=lenstr(vname(3,indxShflx_lat))
          ierr=nf_put_att_text (ncid, vidShflx_lat, 'units',     lvar,
     &                                 vname(3,indxShflx_lat)(1:lvar))

          call nf_add_attribute(ncid, vidShflx_lat, 
     &         vname(:,indxShflx_lat), 5, NF_FOUT, ierr)
        endif
!
! Sensible Shflx [W/m2]
!
      if (wrt(indxShflx_sen)) then
        lvar=lenstr(vname(1,indxShflx_sen))
        ierr=nf_def_var (ncid, vname(1,indxShflx_sen)(1:lvar),
     &                           NF_FOUT, 3, r2dgrd, vidShflx_sen)
#  ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidShflx_sen,nf_collective)
#  endif
#  ifndef AVRH
          lvar=lenstr(vname(2,indxShflx_sen))
          ierr=nf_put_att_text (ncid, vidShflx_sen, 'long_name', lvar,
     &                                 vname(2,indxShflx_sen)(1:lvar))
#  else
          text='averaged '/ /vname(2,indxShflx_sen)
          lvar=lenstr(text)
          ierr=nf_put_att_text (ncid, vidShflx_sen, 'long_name', lvar,
     &                                               text(1:lvar))
#  endif
          lvar=lenstr(vname(3,indxShflx_sen))
          ierr=nf_put_att_text (ncid, vidShflx_sen, 'units',     lvar,
     &                                 vname(3,indxShflx_sen)(1:lvar))

        call nf_add_attribute(ncid, vidShflx_sen, 
     &       vname(:,indxShflx_sen), 5, NF_FOUT, ierr)

      endif
# endif  /* BULK_FLUX && TEMPERATURE  */

# if defined BIOLOGY & !defined PISCES
!
! Depth of the euphotic layer.
!
        if (wrt(indxHel)) then
          lvar=lenstr(vname(1,indxHel))
          ierr=nf_def_var (ncid, vname(1,indxHel)(1:lvar),
     &                             NF_FOUT, 3, r2dgrd, vidHel)
#  ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidHel,nf_collective)
#  endif
#  ifndef AVRH
          lvar=lenstr(vname(2,indxHel))
          ierr=nf_put_att_text (ncid, vidHel, 'long_name', lvar,
     &                                vname(2,indxHel)(1:lvar))
#  else
          text='averaged '/ /vname(2,indxHel)
          lvar=lenstr(text)
          ierr=nf_put_att_text (ncid, vidHel, 'long_name', lvar,
     &                                              text(1:lvar))
#  endif
          lvar=lenstr(vname(3,indxHel))
          ierr=nf_put_att_text (ncid, vidHel, 'units',     lvar,
     &                                  vname(3,indxHel)(1:lvar))
          lvar=lenstr(vname(4,indxHel))
          ierr=nf_put_att_text (ncid, vidHel, 'field',     lvar,
     &                                  vname(4,indxHel)(1:lvar))

          call nf_add_attribute(ncid, vidHel, vname(:,indxHel), 
     &                          5, NF_FOUT, ierr)
        endif
!
! Chlorophyll/Carbon ratio
!
#  ifdef BIO_NChlPZD
        if (wrt(indxChC)) then
          lvar=lenstr(vname(1,indxChC))
          ierr=nf_def_var (ncid, vname(1,indxChC)(1:lvar),
     &                             NF_FOUT, 4, r3dgrd, vidChC)
#   ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidChC,nf_collective)
#   endif
#   ifndef AVRH
          lvar=lenstr(vname(2,indxChC))
          ierr=nf_put_att_text (ncid, vidChC, 'long_name', lvar,
     &                                  vname(2,indxChC)(1:lvar))
#   else
          text='averaged '/ /vname(2,indxChC)
          lvar=lenstr(text)
          ierr=nf_put_att_text (ncid, vidChC, 'long_name', lvar,
     &                                              text(1:lvar))
#   endif
          lvar=lenstr(vname(3,indxChC))
          ierr=nf_put_att_text (ncid, vidChC, 'units',     lvar,
     &                                  vname(3,indxChC)(1:lvar))
          lvar=lenstr(vname(4,indxChC))
          ierr=nf_put_att_text (ncid, vidChC, 'field',     lvar,
     &                                  vname(4,indxChC)(1:lvar))

          call nf_add_attribute(ncid, vidChC, vname(:,indxChC), 
     &                          5, NF_FOUT, ierr)
        endif
#   ifdef OXYGEN
!
! Wind speed at 10 m
!
        if (wrt(indxU10)) then
          lvar=lenstr(vname(1,indxU10))
          ierr=nf_def_var (ncid, vname(1,indxU10)(1:lvar),
     &                             NF_FOUT, 3, r2dgrd, vidU10)
#    ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidU10,nf_collective)
#    endif
#    ifndef AVRH
          lvar=lenstr(vname(2,indxU10))
          ierr=nf_put_att_text (ncid, vidU10, 'long_name', lvar,
     &                                  vname(2,indxU10)(1:lvar))
#    else
          text='averaged '/ /vname(2,indxU10)
          lvar=lenstr(text)
          ierr=nf_put_att_text (ncid, vidU10, 'long_name', lvar,
     &                                              text(1:lvar))
#    endif
          lvar=lenstr(vname(3,indxU10))
          ierr=nf_put_att_text (ncid, vidU10, 'units',     lvar,
     &                                  vname(3,indxU10)(1:lvar))
          lvar=lenstr(vname(4,indxU10))
          ierr=nf_put_att_text (ncid, vidU10, 'field',     lvar,
     &                                  vname(4,indxU10)(1:lvar))

          call nf_add_attribute(ncid, vidU10, vname(:,indxU10), 
     &                          5, NF_FOUT, ierr)
        endif
!
! Gas exchange coefficient of O2
!
        if (wrt(indxKvO2)) then
          lvar=lenstr(vname(1,indxKvO2))
          ierr=nf_def_var (ncid, vname(1,indxKvO2)(1:lvar),
     &                             NF_FOUT, 3, r2dgrd, vidKvO2)
#    ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidKvO2,nf_collective)
#    endif
#    ifndef AVRH
          lvar=lenstr(vname(2,indxKvO2))
          ierr=nf_put_att_text (ncid, vidKvO2, 'long_name', lvar,
     &                                  vname(2,indxKvO2)(1:lvar))
#    else
          text='averaged '/ /vname(2,indxKvO2)
          lvar=lenstr(text)
          ierr=nf_put_att_text (ncid, vidKvO2, 'long_name', lvar,
     &                                              text(1:lvar))
#    endif
          lvar=lenstr(vname(3,indxKvO2))
          ierr=nf_put_att_text (ncid, vidKvO2, 'units',     lvar,
     &                                  vname(3,indxKvO2)(1:lvar))
          lvar=lenstr(vname(4,indxKvO2))
          ierr=nf_put_att_text (ncid, vidKvO2, 'field',     lvar,
     &                                  vname(4,indxKvO2)(1:lvar))

          call nf_add_attribute(ncid, vidKvO2, vname(:,indxKvO2), 5,
     &         NF_FOUT, ierr)
        endif
!
! Saturation concentration of O2
!
        if (wrt(indxO2sat)) then
          lvar=lenstr(vname(1,indxO2sat))
          ierr=nf_def_var (ncid, vname(1,indxO2sat)(1:lvar),
     &                             NF_FOUT, 3, r2dgrd, vidO2sat)
#    ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidO2sat,nf_collective)
#    endif
#    ifndef AVRH
          lvar=lenstr(vname(2,indxO2sat))
          ierr=nf_put_att_text (ncid, vidO2sat, 'long_name', lvar,
     &                                  vname(2,indxO2sat)(1:lvar))
#    else
          text='averaged '/ /vname(2,indxO2sat)
          lvar=lenstr(text)
          ierr=nf_put_att_text (ncid, vidO2sat, 'long_name', lvar,
     &                                              text(1:lvar))
#    endif
          lvar=lenstr(vname(3,indxO2sat))
          ierr=nf_put_att_text (ncid, vidO2sat, 'units',     lvar,
     &                                  vname(3,indxO2sat)(1:lvar))
          lvar=lenstr(vname(4,indxO2sat))
          ierr=nf_put_att_text (ncid, vidO2sat, 'field',     lvar,
     &                                  vname(4,indxO2sat)(1:lvar))

          call nf_add_attribute(ncid, vidO2sat, vname(:,indxO2sat), 5,
     &         NF_FOUT, ierr)
        endif
#   endif /* OXYGEN */
#  endif
# endif
!
# ifdef SEDIMENT
!
! bed_thick  Sediment bed layer thickness (m)
! bed_poros  Bed porosity (-)
! bed_frac_sand/silt
!           volume fraction of size class ised in each bed layer (-)
! dflx_ ; e_flx_
!            Deposition/Erosion fluxes
! bdlu_ ; bdlv_
!
        indxWrk=indxSed-1
        if (wrt(indxWrk)) then
          lvar=lenstr(vname(1,indxWrk))
          ierr=nf_def_var (ncid, vname(1,indxWrk)(1:lvar),
     &                                NF_FOUT, 3, r2dgrd, vidSed(1))
#  ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidSed(1),nf_collective)
#  endif
#  ifndef AVRH
          lvar=lenstr(vname(2,indxWrk))
          ierr=nf_put_att_text (ncid, vidSed(1), 'long_name', lvar,
     &                                  vname(2,indxWrk)(1:lvar))
#  else
          text='averaged '/ /vname(2,indxWrk)
          lvar=lenstr(text)
          ierr=nf_put_att_text (ncid, vidSed(1), 'long_name', lvar,
     &                                            text(1:lvar))
#  endif
          lvar=lenstr(vname(3,indxWrk))
          ierr=nf_put_att_text (ncid, vidSed(1), 'units',     lvar,
     &                                  vname(3,indxWrk)(1:lvar))
          lvar=lenstr(vname(4,indxWrk))
          ierr=nf_put_att_text (ncid, vidSed(1), 'field',     lvar,
     &                                  vname(4,indxWrk)(1:lvar))

          call nf_add_attribute(ncid, vidSed(1), vname(:,indxWrk),
     &          5, NF_FOUT, ierr)
        endif
!
        do itrc=1,2+NST
          indxWrk=indxSed+itrc-1
          if (wrt(indxWrk)) then
            lvar=lenstr(vname(1,indxWrk))
            ierr=nf_def_var (ncid, vname(1,indxWrk)(1:lvar),
     &                             NF_FOUT, 4, b3dgrd, vidSed(itrc+1))
#  ifdef NC4PAR
            ierr=nf_var_par_access(ncid,vidSed(itrc+1),nf_collective)
#  endif
#  ifndef AVRH
            lvar=lenstr(vname(2,indxWrk))
            ierr=nf_put_att_text (ncid,vidSed(itrc+1),'long_name', lvar,
     &                                    vname(2,indxWrk)(1:lvar))
#  else
            text='averaged '/ /vname(2,indxWrk)
            lvar=lenstr(text)
            ierr=nf_put_att_text (ncid,vidSed(itrc+1),'long_name', lvar,
     &                                                text(1:lvar))
#  endif
            lvar=lenstr(vname(3,indxWrk))
            ierr=nf_put_att_text (ncid, vidSed(itrc+1), 'units',   lvar,
     &                                    vname(3,indxWrk)(1:lvar))

            call nf_add_attribute(ncid, vidSed(itrc+1), 
     &          vname(:,indxWrk), 5, NF_FOUT, ierr)
          endif
        enddo

        do itrc=2+NST+1,2+NST
#  ifdef SUSPLOAD
     &                 +2*NST
#  endif
#  ifdef BEDLOAD
     &                 +2*NST
#  endif
          indxWrk=indxSed+itrc-1
          if (wrt(indxWrk)) then
            lvar=lenstr(vname(1,indxWrk))
            ierr=nf_def_var (ncid, vname(1,indxWrk)(1:lvar),
     &                             NF_FOUT, 3, r2dgrd, vidSed(itrc+1))
#  ifdef NC4PAR
            ierr=nf_var_par_access(ncid,vidSed(itrc+1),nf_collective)
#  endif
#  ifndef AVRH
            lvar=lenstr(vname(2,indxWrk))
            ierr=nf_put_att_text (ncid,vidSed(itrc+1),'long_name', lvar,
     &                                    vname(2,indxWrk)(1:lvar))
#  else
            text='averaged '/ /vname(2,indxWrk)
            lvar=lenstr(text)
            ierr=nf_put_att_text (ncid,vidSed(itrc+1),'long_name', lvar,
     &                                                text(1:lvar))
#  endif
            lvar=lenstr(vname(3,indxWrk))
            ierr=nf_put_att_text (ncid,vidSed(itrc+1),'units',     lvar,
     &                                    vname(3,indxWrk)(1:lvar))

            call nf_add_attribute(ncid, vidSed(itrc+1), 
     &            vname(:,indxWrk), 5, NF_FOUT, ierr)
          endif
        enddo
!
#  if defined MIXED_BED || defined COHESIVE_BED
        indxWrk=indxBTCR
        if (wrt(indxWrk)) then
          lvar=lenstr(vname(1,indxWrk))
          ierr=nf_def_var (ncid, vname(1,indxWrk)(1:lvar),
     &                             NF_FOUT, 4, b3dgrd, vidSed(itrc+1))
#   ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidSed(itrc+1),nf_collective)
#   endif
#   ifndef AVRH
          lvar=lenstr(vname(2,indxWrk))
          ierr=nf_put_att_text (ncid, vidSed(itrc+1), 'long_name', lvar,
     &                                  vname(2,indxWrk)(1:lvar))
#   else
          text='averaged '/ /vname(2,indxWrk)
          lvar=lenstr(text)
          ierr=nf_put_att_text (ncid, vidSed(itrc+1), 'long_name', lvar,
     &                                              text(1:lvar))
#   endif
          lvar=lenstr(vname(3,indxWrk))
          ierr=nf_put_att_text (ncid, vidSed(itrc+1), 'units',     lvar,
     &                                  vname(3,indxWrk)(1:lvar))

          call nf_add_attribute(ncid, vidSed(itrc+1), 
     &         vname(:,indxWrk), 5, NF_FOUT, ierr)
        endif
#  endif
!
# endif /* SEDIMENT */

# if defined BBL && !defined AVRH
!
! Abed       Wind-induced, bed wave excursion amplitude (m)
! Hripple    Bed ripple height (m)
! Lripple    Bed ripple length (m)
! Zbnot      Physical hydraulic bottom roughness (grain roughness) (m)
! Zbapp      Effective apparent hydraulic bottom roughness
! bostrw     Wave-induced kinematic bottom stress (N/m2)
!
        do itrc=1,6
          indxWrk=indxBBL+itrc-1
          if (wrt(indxWrk)) then
            lvar=lenstr(vname(1,indxWrk))
            ierr=nf_def_var (ncid, vname(1,indxWrk)(1:lvar),
     &                               NF_FOUT, 3, r2dgrd, vidBBL(itrc))
#  ifdef NC4PAR
            ierr=nf_var_par_access(ncid,vidBBL(itrc),nf_collective)
#  endif
#  ifndef AVRH
            lvar=lenstr(vname(2,indxWrk))
            ierr=nf_put_att_text (ncid, vidBBL(itrc), 'long_name', lvar,
     &                                    vname(2,indxWrk)(1:lvar))
#  else
            text='averaged '/ /vname(2,indxWrk)
            lvar=lenstr(text)
            ierr=nf_put_att_text (ncid, vidBBL(itrc), 'long_name', lvar,
     &                                                text(1:lvar))
#  endif
            lvar=lenstr(vname(3,indxWrk))
            ierr=nf_put_att_text (ncid, vidBBL(itrc), 'units',     lvar,
     &                                    vname(3,indxWrk)(1:lvar))

            call nf_add_attribute(ncid, vidBBL(itrc), 
     &           vname(:,indxWrk), 5, NF_FOUT, ierr)
          endif
        enddo
# endif /* BBL */
# if defined SST_SKIN && defined TEMPERATURE
! Skin temperature
        if (wrt(indxT)) then
          lvar=lenstr(vname(1,indxSST_skin))
          ierr=nf_def_var (ncid, vname(1,indxSST_skin)(1:lvar),
     &                           NF_FOUT, 3, r2dgrd, vidSST_skin)
#  ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidSST_skin,nf_collective)
#  endif
#  ifndef AVRH
          lvar=lenstr(vname(2,indxSST_skin))
          ierr=nf_put_att_text (ncid, vidSST_skin, 'long_name', lvar,
     &                                 vname(2,indxSST_skin)(1:lvar))
#  else
          text='averaged '/ /vname(2,indxSST_skin)
          lvar=lenstr(text)
          ierr=nf_put_att_text (ncid, vidSST_skin, 'long_name', lvar,
     &                                                 text(1:lvar))
#  endif
          lvar=lenstr(vname(3,indxSST_skin))
          ierr=nf_put_att_text (ncid, vidSST_skin, 'units',     lvar,
     &                                 vname(3,indxSST_skin)(1:lvar))
          lvar=lenstr(vname(5,indxSST_skin))

          call nf_add_attribute(ncid, vidSST_skin, 
     &         vname(:,indxSST_skin), 5, NF_FOUT , ierr)
        endif
# endif /* SST_SKIN */
!
!-----------------------------------------------------------------------
!  MUSTANG sediment model output
!-----------------------------------------------------------------------
!
# ifdef MUSTANG
        do indx=1,outMust_nbvar
            if (outMust(indx)) then 
               if (out2DMust(indx)) then
                   ierr=nf_def_var (ncid, vname_Must(1,indx),
     &                 NF_FOUT, 3, r2dgrd, vidMust(indx))  
               elseif (out3DsedMust(indx)) then
                    ierr=nf_def_var (ncid, vname_Must(1,indx),
     &                 NF_FOUT, 4, b3dgrd, vidMust(indx))
               else
                   ierr=nf_def_var (ncid, vname_Must(1,indx),
     &                 NF_FOUT, 4, r3dgrd, vidMust(indx)) 
               endif
#   ifdef NC4PAR
                ierr=nf_var_par_access(ncid,vidMust(indx),
     &                           nf_collective)
#   endif
                lvar=lenstr(vname_Must(2,indx))
                ierr=nf_put_att_text (ncid, vidMust(indx), 'long_name',
     &                                  lvar,vname_Must(2,indx)(1:lvar))
                lvar=lenstr(vname_Must(3,indx))
                ierr=nf_put_att_text (ncid, vidMust(indx), 'units',
     &                                  lvar,vname_Must(3,indx)(1:lvar))
                lvar=lenstr(vname_Must(4,indx))
                ierr=nf_put_att_text (ncid, vidMust(indx), 'field',
     &                                  lvar,vname_Must(4,indx)(1:lvar))
                call nf_add_attribute(ncid, vidMust(indx),
     &              vname_Must(:,indx), 5, NF_FOUT, ierr)
            endif
        enddo

# endif /* MUSTANG */
!
!-----------------------------------------------------------------------
!  OBSTRUCTION model output
!-----------------------------------------------------------------------
!
# ifdef OBSTRUCTION

        do itrv=1,17*obst_nbvar+15
            indx = itrv
            if (outObst(indx)) then 
               if (out2DObst(indx)) then
                   ierr=nf_def_var (ncid, vname_obst(1,indx),
     &                 NF_FOUT, 3, r2dgrd, vidObst(indx))  
               else
                   ierr=nf_def_var (ncid, vname_obst(1,indx),
     &                 NF_FOUT, 4, r3dgrd, vidObst(indx)) 
               endif
#   ifdef NC4PAR
                ierr=nf_var_par_access(ncid,vidObst(indx),
     &                           nf_collective)
#   endif
                lvar=lenstr(vname_obst(2,indx))
                ierr=nf_put_att_text (ncid, vidObst(indx), 'long_name',
     &                                  lvar,vname_obst(2,indx)(1:lvar))
                lvar=lenstr(vname_obst(3,indx))
                ierr=nf_put_att_text (ncid, vidObst(indx), 'units',
     &                                  lvar,vname_obst(3,indx)(1:lvar))
                lvar=lenstr(vname_obst(4,indx))
                ierr=nf_put_att_text (ncid, vidObst(indx), 'field',
     &                                  lvar,vname_obst(4,indx)(1:lvar))
                call nf_add_attribute(ncid, vidObst(indx),
     &              vname_obst(:,indx), 5, NF_FOUT, ierr)
            endif
        enddo

# endif /* OBSTRUCTION */

#endif /* SOLVE3D */

#ifdef WAVE_IO
!
! Wave prognostic/diagnostic variables in XI/ETA direction.
! =========================================================
!
# ifndef WAVE_ROLLER
        do iwkb=1,7
# else
        do iwkb=1,9
# endif
          if (wrt(indxHRM+iwkb-1)) then
            lvar=lenstr(vname(1,indxHRM+iwkb-1))
            ierr=nf_def_var (ncid, vname(1,indxHRM+iwkb-1)(1:lvar),
     &                            NF_FOUT, 3, r2dgrd, vidWAVE(iwkb))
# ifdef NC4PAR
            ierr=nf_var_par_access(ncid,vidWAVE(iwkb),nf_collective)
# endif
# ifdef AVRH
            text='averaged '/ /vname(2,indxHRM+iwkb-1)
# else
            text=vname(2,indxHRM+iwkb-1)
# endif
            lvar=lenstr(text)
            ierr=nf_put_att_text (ncid, vidWAVE(iwkb), 'long_name',
     &                                         lvar, text(1:lvar))
            lvar=lenstr(vname(3,indxHRM+iwkb-1))
            ierr=nf_put_att_text (ncid, vidWAVE(iwkb), 'units', lvar,
     &                               vname(3,indxHRM+iwkb-1)(1:lvar))
          endif
        enddo
#endif /* WAVE_IO */
!
#ifdef MRL_WCI
!
! Wave-current interaction diagnostic variables.
! =============================================
!
! sup: quasi-static sea-level response (wave set-up/down)
!
        if (wrt(indxSUP)) then
          lvar=lenstr(vname(1,indxSUP))
          ierr=nf_def_var (ncid, vname(1,indxSUP)(1:lvar), NF_FOUT,
     &                                           3, r2dgrd, vidSUP)
# ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidSUP,nf_collective)
# endif

# ifndef AVRH
          lvar=lenstr(vname(2,indxSUP))
          ierr=nf_put_att_text (ncid, vidSUP, 'long_name', lvar,
     &                                  vname(2,indxSUP)(1:lvar))
# else
          text='averaged '/ /vname(2,indxSUP)
          lvar=lenstr(text)
          ierr=nf_put_att_text (ncid, vidSUP, 'long_name', lvar,
     &                                              text(1:lvar))
# endif
          lvar=lenstr(vname(3,indxSUP))
          ierr=nf_put_att_text (ncid, vidSUP, 'units',     lvar,
     &                                 vname(3,indxSUP)(1:lvar))
        endif
!
! ust2d & vst2d: 2D, depth-averaged Stokes drift velocities.
!
        if (wrt(indxUST2D)) then
          lvar=lenstr(vname(1,indxUST2D))
          ierr=nf_def_var (ncid, vname(1,indxUST2D)(1:lvar), NF_FOUT,
     &                                           3, u2dgrd, vidUST2D)
# ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidUST2D,nf_collective)
# endif
# ifndef AVRH
          lvar=lenstr(vname(2,indxUST2D))
          ierr=nf_put_att_text (ncid, vidUST2D, 'long_name', lvar,
     &                                  vname(2,indxUST2D)(1:lvar))
# else
          text='averaged '/ /vname(2,indxUST2D)
          lvar=lenstr(text)
          ierr=nf_put_att_text (ncid, vidUST2D, 'long_name', lvar,
     &                                              text(1:lvar))
# endif
          lvar=lenstr(vname(3,indxUST2D))
          ierr=nf_put_att_text (ncid, vidUST2D, 'units',     lvar,
     &                                 vname(3,indxUST2D)(1:lvar))
        endif
!
        if (wrt(indxVST2D)) then
          lvar=lenstr(vname(1,indxVST2D))
          ierr=nf_def_var (ncid, vname(1,indxVST2D)(1:lvar), NF_FOUT,
     &                                           3, v2dgrd, vidVST2D)
# ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidVST2D,nf_collective)
# endif
# ifndef AVRH
          lvar=lenstr(vname(2,indxVST2D))
          ierr=nf_put_att_text (ncid, vidVST2D, 'long_name', lvar,
     &                                  vname(2,indxVST2D)(1:lvar))
# else
          text='averaged '/ /vname(2,indxVST2D)
          lvar=lenstr(text)
          ierr=nf_put_att_text (ncid, vidVST2D, 'long_name', lvar,
     &                                              text(1:lvar))
# endif
          lvar=lenstr(vname(3,indxVST2D))
          ierr=nf_put_att_text (ncid, vidVST2D, 'units',     lvar,
     &                                  vname(3,indxVST2D)(1:lvar))
        endif
# ifdef SOLVE3D
!
! ust3d & vst3d: 3D Stokes drift velocities.
!
        if (wrt(indxUST)) then
          lvar=lenstr(vname(1,indxUST))
          ierr=nf_def_var (ncid, vname(1,indxUST)(1:lvar), NF_FOUT,
     &                                           4, u3dgrd, vidUST)
#  ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidUST,nf_collective)
#  endif
#  ifndef AVRH
          lvar=lenstr(vname(2,indxUST))
          ierr=nf_put_att_text (ncid, vidUST, 'long_name', lvar,
     &                                  vname(2,indxUST)(1:lvar))
#  else
          text='averaged '/ /vname(2,indxUST)
          lvar=lenstr(text)
          ierr=nf_put_att_text (ncid, vidUST, 'long_name', lvar,
     &                                              text(1:lvar))
#  endif
          lvar=lenstr(vname(3,indxUST))
          ierr=nf_put_att_text (ncid, vidUST, 'units',     lvar,
     &                                 vname(3,indxUST)(1:lvar))
        endif
!
        if (wrt(indxVST)) then
          lvar=lenstr(vname(1,indxVST))
          ierr=nf_def_var (ncid, vname(1,indxVST)(1:lvar), NF_FOUT,
     &                                           4, v3dgrd, vidVST)
#  ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidVST,nf_collective)
#  endif
#  ifndef AVRH
          lvar=lenstr(vname(2,indxVST))
          ierr=nf_put_att_text (ncid, vidVST, 'long_name', lvar,
     &                                  vname(2,indxVST)(1:lvar))
#  else
          text='averaged '/ /vname(2,indxUST)
          lvar=lenstr(text)
          ierr=nf_put_att_text (ncid, vidVST, 'long_name', lvar,
     &                                              text(1:lvar))
#  endif
          lvar=lenstr(vname(3,indxVST))
          ierr=nf_put_att_text (ncid, vidVST, 'units',     lvar,
     &                                 vname(3,indxVST)(1:lvar))
        endif
!
! wst: vertical Stokes drift velocity (m/s) at rho-point
!
        if (wrt(indxWST)) then
          lvar=lenstr(vname(1,indxWST))
          ierr=nf_def_var (ncid, vname(1,indxWST)(1:lvar), NF_FOUT,
     &                                           4, r3dgrd, vidWST)
#  ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidWST,nf_collective)
#  endif
#  ifndef AVRH
          lvar=lenstr(vname(2,indxWST))
          ierr=nf_put_att_text (ncid, vidWST, 'long_name', lvar,
     &                                  vname(2,indxWST)(1:lvar))
#  else
          text='averaged '/ /vname(2,indxWST)
          lvar=lenstr(text)
          ierr=nf_put_att_text (ncid, vidWST, 'long_name', lvar,
     &                                              text(1:lvar))
#  endif
          lvar=lenstr(vname(3,indxWST))
          ierr=nf_put_att_text (ncid, vidWST, 'units',     lvar,
     &                                 vname(3,indxWST)(1:lvar))
        endif
!
!
! Akb: Vertical eddy viscosity coefficient due to depth-induced wave breaking.
!
        if (wrt(indxAkb)) then
          lvar=lenstr(vname(1,indxAkb))
          ierr=nf_def_var (ncid, vname(1,indxAkb)(1:lvar), NF_FOUT,
     &                                           4, w3dgrd, vidAkb)
#  ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidAkb,nf_collective)
#  endif
#  ifndef AVRH
          lvar=lenstr(vname(2,indxAkb))
          ierr=nf_put_att_text (ncid, vidAkb, 'long_name', lvar,
     &                                  vname(2,indxAkb)(1:lvar))
#  else
          text='averaged '/ /vname(2,indxAkb)
          lvar=lenstr(text)
          ierr=nf_put_att_text (ncid, vidAkb, 'long_name', lvar,
     &                                              text(1:lvar))
#  endif
          lvar=lenstr(vname(3,indxAkb))
          ierr=nf_put_att_text (ncid, vidAkb, 'units',     lvar,
     &                                  vname(3,indxAkb)(1:lvar))
        endif
!
! Akw: Vertical eddy diffusivity coefficient due to primary waves.
!
        if (wrt(indxAkw)) then
          lvar=lenstr(vname(1,indxAkw))
          ierr=nf_def_var (ncid, vname(1,indxAkw)(1:lvar), NF_FOUT,
     &                                           4, w3dgrd, vidAkw)
#  ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidAkw,nf_collective)
#  endif
#  ifndef AVRH
          lvar=lenstr(vname(2,indxAkw))
          ierr=nf_put_att_text (ncid, vidAkw, 'long_name', lvar,
     &                                  vname(2,indxAkw)(1:lvar))
#  else
          text='averaged '/ /vname(2,indxAkw)
          lvar=lenstr(text)
          ierr=nf_put_att_text (ncid, vidAkw, 'long_name', lvar,
     &                                              text(1:lvar))
#  endif
          lvar=lenstr(vname(3,indxAkw))
          ierr=nf_put_att_text (ncid, vidAkw, 'units',     lvar,
     &                                  vname(3,indxAkw)(1:lvar))
        endif
!
! kvf: vertical vortex force term (u^St du/dz)
!
        if (wrt(indxKVF)) then
          lvar=lenstr(vname(1,indxKVF))
          ierr=nf_def_var (ncid, vname(1,indxKVF)(1:lvar), NF_FOUT,
     &                                           4, r3dgrd, vidKVF)
#  ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidKVF,nf_collective)
#  endif
#  ifndef AVRH
          lvar=lenstr(vname(2,indxKVF))
          ierr=nf_put_att_text (ncid, vidKVF, 'long_name', lvar,
     &                                  vname(2,indxKVF)(1:lvar))
#  else
          text='averaged '/ /vname(2,indxKVF)
          lvar=lenstr(text)
          ierr=nf_put_att_text (ncid, vidKVF, 'long_name', lvar,
     &                                              text(1:lvar))
#  endif
          lvar=lenstr(vname(3,indxKVF))
          ierr=nf_put_att_text (ncid, vidKVF, 'units',     lvar,
     &                                 vname(3,indxKVF)(1:lvar))
        endif
!
! calP: surface pressure correction term appeared in prsgrd term
!
        if (wrt(indxCALP)) then
          lvar=lenstr(vname(1,indxCALP))
          ierr=nf_def_var (ncid, vname(1,indxCALP)(1:lvar), NF_FOUT,
     &                                           3, r2dgrd, vidCALP)
#  ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidCALP,nf_collective)
#  endif
#  ifndef AVRH
          lvar=lenstr(vname(2,indxCALP))
          ierr=nf_put_att_text (ncid, vidCALP, 'long_name', lvar,
     &                                  vname(2,indxCALP)(1:lvar))
#  else
          text='averaged '/ /vname(2,indxCALP)
          lvar=lenstr(text)
          ierr=nf_put_att_text (ncid, vidCALP, 'long_name', lvar,
     &                                              text(1:lvar))
#  endif
          lvar=lenstr(vname(3,indxCALP))
          ierr=nf_put_att_text (ncid, vidCALP, 'units',     lvar,
     &                                 vname(3,indxCALP)(1:lvar))
        endif
!
! Kapsrf: surface Bernoulli head term appeared in prsgrd term
!
        if (wrt(indxKAPS)) then
          lvar=lenstr(vname(1,indxKAPS))
          ierr=nf_def_var (ncid, vname(1,indxKAPS)(1:lvar), NF_FOUT,
     &                                           3, r2dgrd, vidKAPS)
#  ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidKAPS,nf_collective)
#  endif
#  ifndef AVRH
          lvar=lenstr(vname(2,indxKAPS))
          ierr=nf_put_att_text (ncid, vidKAPS, 'long_name', lvar,
     &                                  vname(2,indxKAPS)(1:lvar))
#  else
          text='averaged '/ /vname(2,indxKAPS)
          lvar=lenstr(text)
          ierr=nf_put_att_text (ncid, vidKAPS, 'long_name', lvar,
     &                                              text(1:lvar))
#  endif
          lvar=lenstr(vname(3,indxKAPS))
          ierr=nf_put_att_text (ncid, vidKAPS, 'units',     lvar,
     &                                 vname(3,indxKAPS)(1:lvar))
        endif
# endif  /* SOLVE3D */
#endif   /* MRL_WCI */

!
! Leave definition mode.
!
!
        ierr=nf_enddef(ncid)
        MPI_master_only write(stdout,'(6x,4A,1x,A,i4)')
     &                'DEF_HIS/AVG - Created ',
     &                'new netCDF file ''', ncname(1:lstr), '''.'
     &                 MYID
!
!===================================================================
! Open an existing file and prepare for appending data.
! ==== == ======== ==== === ======= === ========= =====
! Inquire about the dimensions and variables. Check for
! consistency with model dimensions. In the case when file
! is rejected (whether it cannot be opened, or something
! is wrong with its dimensions) create a new file.
!
! After that verify that all necessary variables are already
! defined there and find their netCDF IDs.
!===================================================================
!
      elseif (ncid.eq.-1) then
#ifndef NC4PAR
        ierr=nf_open (ncname(1:lstr), nf_write, ncid)
#else
        ierr=nf_open_par (ncname(1:lstr), IOR(nf_write, nf_mpiio),
     &     MPI_COMM_WORLD, MPI_INFO_NULL, ncid)
#endif
        if (ierr .eq. nf_noerr) then
          ierr=checkdims (ncid, ncname, lstr, rec)
          if (ierr .eq. nf_noerr) then
            if (rec_per_file.eq.0) then
              ierr=rec+1 - total_rec
            else
              ierr=rec+1 - (1+mod(total_rec-1, rec_per_file))
            endif
            if (ierr.gt.0) then
              MPI_master_only write( stdout,
     &                 '(/1x,A,I5,1x,A/8x,3A,I5,/8x,A,I5,1x,A/)'
     &            ) 'DEF_HIS/AVG WARNING: Actual number of records',
     &               rec,  'in netCDF file',  '''',  ncname(1:lstr),
     &             ''' exceeds the record number from restart data',
     &             rec+1-ierr,'/', total_rec,', restart is assumed.'
              rec=rec-ierr
            elseif (rec_per_file.eq.0) then
              total_rec=rec+1           ! <-- set to the next record
#if defined MPI & !defined PARALLEL_FILES
              if (mynode.gt.0) total_rec=total_rec-1
#endif
            endif
            ierr=nf_noerr
          endif
        endif
!
        if (ierr .ne. nf_noerr) then
#if defined MPI & !defined PARALLEL_FILES & !defined NC4PAR
          if (mynode.eq.0) then
            create_new_file=.true.
            goto 10
          else
            MPI_master_only write(stdout,'(/1x,4A,2x,A,I4/)')
     &                  'DEF_HIS/AVG ERROR: ',
     &                  'Cannot open file ''', ncname(1:lstr), '''.'
     &                   MYID
            goto 99                                       !--> ERROR
          endif
#else
          create_new_file=.true.
          goto 10
#endif
        endif
!
! Find netCDF IDs of evolving model variables:
! ---- ------ --- -- -------- ----- ----------
!
! Time step indices:
!
        ierr=nf_inq_varid (ncid, 'time_step', vidTstep)
        if (ierr .ne. nf_noerr) then
          write(stdout,1) 'time_step', ncname(1:lstr)
          goto 99                                         !--> ERROR
        endif
#ifdef NC4PAR
        ierr=nf_var_par_access(ncid,vidTstep,nf_collective)
#endif
!
! Time.
!
        lvar=lenstr(vname(1,indxTime))
        ierr=nf_inq_varid (ncid,vname(1,indxTime)(1:lvar),vidTime)
        if (ierr .ne. nf_noerr) then
          write(stdout,1) vname(1,indxTime)(1:lvar), ncname(1:lstr)
          goto 99                                         !--> ERROR
        endif
#ifdef NC4PAR
        ierr=nf_var_par_access(ncid,vidTime,nf_collective)
#endif

!
! Time2.
!
        lvar=lenstr(vname(1,indxTime2))
        ierr=nf_inq_varid (ncid,vname(1,indxTime2)(1:lvar),vidTime2)
        if (ierr .ne. nf_noerr) then
          write(stdout,1) vname(1,indxTime2)(1:lvar), ncname(1:lstr)
          goto 99                                         !--> ERROR
        endif
#ifdef NC4PAR
        ierr=nf_var_par_access(ncid,vidTime2,nf_collective)
#endif

!
! Free-surface.
!
        if (wrt(indxZ)) then
          lvar=lenstr(vname(1,indxZ))
          ierr=nf_inq_varid (ncid, vname(1,indxZ)(1:lvar), vidZ)
          if (ierr .ne. nf_noerr) then
            write(stdout,1) vname(1,indxZ)(1:lvar), ncname(1:lstr)
            goto 99                                       !--> ERROR
          endif
#ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidZ,nf_collective)
#endif
        endif
!
! 2D momenta in XI- and ETA-directions.
!
        if (wrt(indxUb)) then
          lvar=lenstr(vname(1,indxUb))
          ierr=nf_inq_varid (ncid, vname(1,indxUb)(1:lvar), vidUb)
          if (ierr .ne. nf_noerr) then
            write(stdout,1) vname(1,indxUb)(1:lvar), ncname(1:lstr)
            goto 99                                       !--> ERROR
          endif
#ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidUb,nf_collective)
#endif
        endif

        if (wrt(indxVb)) then
          lvar=lenstr(vname(1,indxVb))
          ierr=nf_inq_varid (ncid, vname(1,indxVb)(1:lvar), vidVb)
          if (ierr .ne. nf_noerr) then
            write(stdout,1) vname(1,indxVb)(1:lvar), ncname(1:lstr)
            goto 99                                       !--> ERROR
          endif
#ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidVb,nf_collective)
#endif
        endif
!
! time evolving bathymetry
!
#ifdef MORPHODYN
        if (wrt(indxHm)) then
          lvar=lenstr(vname(1,indxHm))
          ierr=nf_inq_varid (ncid, vname(1,indxHm)(1:lvar), vidHm)
          if (ierr .ne. nf_noerr) then
            write(stdout,1) vname(1,indxHm)(1:lvar), ncname(1:lstr)
            goto 99                                       !--> ERROR
          endif
# ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidHm,nf_collective)
# endif
        endif
#endif
!
! bostr     Kinematic bottom stress (N/m2)
!
        if (wrt(indxBostr)) then
          lvar=lenstr(vname(1,indxBostr))
          ierr=nf_inq_varid (ncid,vname(1,indxBostr)(1:lvar),
     &                                                   vidBostr)
          if (ierr .ne. nf_noerr) then
            write(stdout,1) vname(1,indxBostr)(1:lvar), ncname(1:lstr)
            goto 99                                       !--> ERROR
          endif
#ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidBostr,nf_collective)
#endif
        endif


!
! bustr     Kinematic U bottom stress (N/m2)
!
        if (wrt(indxBustr)) then
          lvar=lenstr(vname(1,indxBustr))
          ierr=nf_inq_varid (ncid,vname(1,indxBustr)(1:lvar),
     &                                                   vidBustr)
          if (ierr .ne. nf_noerr) then
          write(stdout,1) vname(1,indxBustr)(1:lvar), ncname(1:lstr)
            goto 99                                       !--> ERROR
          endif
#ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidBustr,nf_collective)
#endif
        endif


!
! bvstr     Kinematic V bottom stress (N/m2)
!
        if (wrt(indxBvstr)) then
          lvar=lenstr(vname(1,indxBvstr))
          ierr=nf_inq_varid (ncid,vname(1,indxBvstr)(1:lvar),
     &                                                   vidBvstr)
          if (ierr .ne. nf_noerr) then
          write(stdout,1) vname(1,indxBvstr)(1:lvar), ncname(1:lstr)
            goto 99                                       !--> ERROR
          endif
#ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidBvstr,nf_collective)
#endif
        endif
!
! wstr     Kinematic surface wind stress (N/m2)
!
        if (wrt(indxWstr)) then
          lvar=lenstr(vname(1,indxWstr))
          ierr=nf_inq_varid (ncid,vname(1,indxWstr)(1:lvar),
     &                                                   vidWstr)
          if (ierr .ne. nf_noerr) then
            write(stdout,1) vname(1,indxWstr)(1:lvar), ncname(1:lstr)
            goto 99                                       !--> ERROR
          endif
#ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidWstr,nf_collective)
#endif
        endif
!
! uwstr     Kinematic surface u wind stress component (N/m2)
!
        if (wrt(indxUWstr)) then
          lvar=lenstr(vname(1,indxUWstr))
          ierr=nf_inq_varid (ncid,vname(1,indxUWstr)(1:lvar),
     &                                                   vidUWstr)
          if (ierr .ne. nf_noerr) then
            write(stdout,1) vname(1,indxUWstr)(1:lvar), ncname(1:lstr)
            goto 99                                       !--> ERROR
          endif
#ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidUWstr,nf_collective)
#endif
        endif
!
! vwstr     Kinematic surface v wind stress component (N/m2)
!
        if (wrt(indxVWstr)) then
          lvar=lenstr(vname(1,indxVWstr))
          ierr=nf_inq_varid (ncid,vname(1,indxVWstr)(1:lvar),
     &                                                   vidVWstr)
          if (ierr .ne. nf_noerr) then
            write(stdout,1) vname(1,indxVWstr)(1:lvar), ncname(1:lstr)
            goto 99                                       !--> ERROR
          endif
#ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidVWstr,nf_collective)
#endif
        endif
#ifdef SOLVE3D
!
! 3D momenta in XI- and ETA-directions.
!
#ifdef M3FAST_HIS
!
! M3FAST-kernel related variables:
!
        if (wrt(indxRnbq)) then
          lvar=lenstr(vname(1,indxRnbq))
          ierr=nf_inq_varid (ncid, vname(1,indxRnbq)(1:lvar), vidRnbq)
          if (ierr .ne. nf_noerr) then
            write(stdout,1) vname(1,indxRnbq)(1:lvar), ncname(1:lstr)
            goto 99                                       !--> ERROR
          endif
        endif

       if (wrt(indxUnbq)) then
          lvar=lenstr(vname(1,indxUnbq))
          ierr=nf_inq_varid (ncid, vname(1,indxUnbq)(1:lvar), vidUnbq)
          if (ierr .ne. nf_noerr) then
            write(stdout,1) vname(1,indxUnbq)(1:lvar), ncname(1:lstr)
            goto 99                                       !--> ERROR
          endif
        endif

        if (wrt(indxVnbq)) then
          lvar=lenstr(vname(1,indxVnbq))
          ierr=nf_inq_varid (ncid, vname(1,indxVnbq)(1:lvar), vidVnbq)
          if (ierr .ne. nf_noerr) then
            write(stdout,1) vname(1,indxVnbq)(1:lvar), ncname(1:lstr)
            goto 99                                       !--> ERROR
          endif
        endif
        
        if (wrt(indxWnbq)) then
          lvar=lenstr(vname(1,indxWnbq))
          ierr=nf_inq_varid (ncid, vname(1,indxWnbq)(1:lvar), vidWnbq)
          if (ierr .ne. nf_noerr) then
            write(stdout,1) vname(1,indxWnbq)(1:lvar), ncname(1:lstr)
            goto 99                                       !--> ERROR
          endif
        endif

        if (wrt(indxCnbq)) then
          lvar=lenstr(vname(1,indxCnbq))
          ierr=nf_inq_varid (ncid, vname(1,indxCnbq)(1:lvar), vidCnbq)
          if (ierr .ne. nf_noerr) then
            write(stdout,1) vname(1,indxCnbq)(1:lvar), ncname(1:lstr)
            goto 99                                       !--> ERROR
          endif
        endif
#endif        
        if (wrt(indxU)) then
          lvar=lenstr(vname(1,indxU))
          ierr=nf_inq_varid (ncid, vname(1,indxU)(1:lvar), vidU)
          if (ierr .ne. nf_noerr) then
            write(stdout,1) vname(1,indxU)(1:lvar), ncname(1:lstr)
            goto 99                                       !--> ERROR
          endif
# ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidU,nf_collective)
# endif
        endif
!
        if (wrt(indxV)) then
          lvar=lenstr(vname(1,indxV))
          ierr=nf_inq_varid (ncid, vname(1,indxV)(1:lvar), vidV)
          if (ierr .ne. nf_noerr) then
            write(stdout,1) vname(1,indxV)(1:lvar), ncname(1:lstr)
            goto 99                                       !--> ERROR
          endif
# ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidV,nf_collective)
# endif
        endif
!
! Tracer variables.
!
# ifdef TRACERS
        do itrc=1,NT
          if (wrt(indxV+itrc)) then
            lvar=lenstr(vname(1,indxV+itrc))
            ierr=nf_inq_varid (ncid, vname(1,indxV+itrc)(1:lvar),
     &                                                 vidT(itrc))
            if (ierr .ne. nf_noerr) then
              write(stdout,1) vname(1,indxV+itrc)(1:lvar),
     &                                       ncname(1:lstr)
              goto 99                                     !--> ERROR
            endif
#  ifdef NC4PAR
            ierr=nf_var_par_access(ncid,vidT(itrc),nf_collective)
#  endif
          endif
        enddo
# endif
!
! Density anomaly.
!
        if (wrt(indxR)) then
          lvar=lenstr(vname(1,indxR))
          ierr=nf_inq_varid (ncid, vname(1,indxR)(1:lvar), vidR)
          if (ierr .ne. nf_noerr) then
            write(stdout,1) vname(1,indxR)(1:lvar), ncname(1:lstr)
            goto 99                                       !--> ERROR
          endif
# ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidR,nf_collective)
# endif
        endif
# if defined ANA_VMIX || defined BVF_MIXING \
  || defined LMD_MIXING || defined LMD_SKPP || defined LMD_BKPP \
  || defined GLS_MIXING
!-----------------------------------------------------------------------------
!
! Brunt Vaisala Frequency.
!
        if (wrt(indxbvf)) then
          lvar=lenstr(vname(1,indxbvf))
          ierr=nf_inq_varid (ncid, vname(1,indxbvf)(1:lvar), vidbvf)
          if (ierr .ne. nf_noerr) then
            write(stdout,1) vname(1,indxbvf)(1:lvar), ncname(1:lstr)
            goto 99                                       !--> ERROR
          endif
#  ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidbvf,nf_collective)
#  endif
        endif
# endif
!-----------------------------------------------------------------------------
! AOU and pCO2
!
# if defined BIOLOGY && defined BIO_BioEBUS
        if (wrt(indxAOU)) then
          lvar=lenstr(vname(1,indxAOU))
          ierr=nf_inq_varid (ncid, vname(1,indxAOU)(1:lvar),
     &                                               vidAOU)
          if (ierr .ne. nf_noerr) then
            write(stdout,1) vname(1,indxAOU)(1:lvar), ncname(1:lstr)
            goto 99                                       !--> ERROR
          endif
#  ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidAOU,nf_collective)
#  endif
        endif
        if (wrt(indxWIND10)) then
          lvar=lenstr(vname(1,indxWIND10))
          ierr=nf_inq_varid (ncid, vname(1,indxWIND10)(1:lvar),
     &                                             vidwind10)
          if (ierr .ne. nf_noerr) then
            write(stdout,1) vname(1,indxWIND10)(1:lvar),ncname(1:lstr)
            goto 99                                       !--> ERROR
          endif
#  ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidwind10,nf_collective)
#  endif
        endif
# endif
!
! S-coordinate "omega" vertical velocity.
!
        if (wrt(indxO)) then
          lvar=lenstr(vname(1,indxO))
          ierr=nf_inq_varid (ncid, vname(1,indxO)(1:lvar), vidO)
          if (ierr .ne. nf_noerr) then
            write(stdout,1) vname(1,indxO)(1:lvar), ncname(1:lstr)
            goto 99                                       !--> ERROR
          endif
# ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidO,nf_collective)
# endif
        endif
!
! True vertical velocity.
!
        if (wrt(indxW)) then
          lvar=lenstr(vname(1,indxW))
          ierr=nf_inq_varid (ncid, vname(1,indxW)(1:lvar), vidW)
          if (ierr .ne. nf_noerr) then
            write(stdout,1) vname(1,indxW)(1:lvar), ncname(1:lstr)
            goto 99                                       !--> ERROR
          endif
# ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidW,nf_collective)
# endif
        endif
# ifdef VIS_COEF_3D
!
! Horizontal viscosity coefficient.
!
        if (wrt(indxVisc)) then
          lvar=lenstr(vname(1,indxVisc))
          ierr=nf_inq_varid (ncid, vname(1,indxVisc)(1:lvar), vidVisc)
          if (ierr .ne. nf_noerr) then
            write(stdout,1) vname(1,indxVisc)(1:lvar), ncname(1:lstr)
            goto 99                                       !--> ERROR
          endif
#  ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidVisc,nf_collective)
#  endif
        endif
# endif
# ifdef DIF_COEF_3D
!
! Horizontal diffusivity coefficient.
!
        if (wrt(indxDiff)) then
          lvar=lenstr(vname(1,indxDiff))
          ierr=nf_inq_varid (ncid, vname(1,indxDiff)(1:lvar), vidDiff)
          if (ierr .ne. nf_noerr) then
            write(stdout,1) vname(1,indxDiff)(1:lvar), ncname(1:lstr)
            goto 99                                       !--> ERROR
          endif
#  ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidDiff,nf_collective)
#  endif
        endif
# endif
!
! Vertical viscosity coefficient.
!
        if (wrt(indxAkv)) then
          lvar=lenstr(vname(1,indxAkv))
          ierr=nf_inq_varid (ncid, vname(1,indxAkv)(1:lvar), vidAkv)
          if (ierr .ne. nf_noerr) then
            write(stdout,1) vname(1,indxAkv)(1:lvar), ncname(1:lstr)
            goto 99                                       !--> ERROR
          endif
# ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidAkv,nf_collective)
# endif
        endif
!
! Vertical diffusion coefficient for potential temperature.
!
# ifdef TEMPERATURE
        if (wrt(indxAkt)) then
          lvar=lenstr(vname(1,indxAkt))
          ierr=nf_inq_varid (ncid,vname(1,indxAkt)(1:lvar), vidAkt)
          if (ierr .ne. nf_noerr) then
            write(stdout,1) vname(1,indxAkt)(1:lvar), ncname(1:lstr)
            goto 99                                       !--> ERROR
          endif
#  ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidAkt,nf_collective)
#  endif
        endif
# endif

# ifdef SALINITY
!
! Vertical diffusion coefficient for salinity.
!
        if (wrt(indxAks)) then
          lvar=lenstr(vname(1,indxAks))
          ierr=nf_inq_varid (ncid,vname(1,indxAks)(1:lvar), vidAks)
          if (ierr .ne. nf_noerr) then
            write(stdout,1) vname(1,indxAks)(1:lvar), ncname(1:lstr)
            goto 99                                       !--> ERROR
          endif
#  ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidAks,nf_collective)
#  endif
        endif
# endif


# if defined LMD_SKPP || defined GLS_MIXING
!
! Depth of surface planetary boundary layer.
!
        if (wrt(indxHbl)) then
          lvar=lenstr(vname(1,indxHbl))
          ierr=nf_inq_varid (ncid,vname(1,indxHbl)(1:lvar), vidHbl)
          if (ierr .ne. nf_noerr) then
            write(stdout,1) vname(1,indxHbl)(1:lvar), ncname(1:lstr)
            goto 99                                       !--> ERROR
          endif
#  ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidHbl,nf_collective)
#  endif
        endif
# endif
# ifdef LMD_BKPP
!
! Depth of bottom planetary boundary layer.
!
        if (wrt(indxHbbl)) then
          lvar=lenstr(vname(1,indxHbbl))
          ierr=nf_inq_varid (ncid,vname(1,indxHbbl)(1:lvar), vidHbbl)
          if (ierr .ne. nf_noerr) then
            write(stdout,1) vname(1,indxHbbl)(1:lvar), ncname(1:lstr)
            goto 99                                       !--> ERROR
          endif
#  ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidHbbl,nf_collective)
#  endif
        endif
# endif

# ifdef GLS_MIXING
!
! Turbulent Kinetic Energy.
!
        if (wrt(indxTke)) then
          lvar=lenstr(vname(1,indxTke))
          ierr=nf_inq_varid (ncid,vname(1,indxTke)(1:lvar), vidTke)
          if (ierr .ne. nf_noerr) then
            write(stdout,1) vname(1,indxTke)(1:lvar), ncname(1:lstr)
            goto 99                                       !--> ERROR
          endif
#  ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidTke,nf_collective)
#  endif
        endif
!
! Generic Length Scale.
!
        if (wrt(indxGls)) then
          lvar=lenstr(vname(1,indxGls))
          ierr=nf_inq_varid (ncid,vname(1,indxGls)(1:lvar), vidGls)
          if (ierr .ne. nf_noerr) then
            write(stdout,1) vname(1,indxGls)(1:lvar), ncname(1:lstr)
            goto 99                                       !--> ERROR
          endif
#  ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidGls,nf_collective)
#  endif
        endif
!
! Turbulent Length Scale.
!
        if (wrt(indxLsc)) then
          lvar=lenstr(vname(1,indxLsc))
          ierr=nf_inq_varid (ncid,vname(1,indxLsc)(1:lvar), vidLsc)
          if (ierr .ne. nf_noerr) then
            write(stdout,1) vname(1,indxLsc)(1:lvar), ncname(1:lstr)
            goto 99                                       !--> ERROR
          endif
#  ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidLsc,nf_collective)
#  endif
        endif
# endif

!
! Shflx     Surface net heat flux [W/m2]
!
# ifdef TEMPERATURE
        if (wrt(indxShflx)) then
          lvar=lenstr(vname(1,indxShflx))
          ierr=nf_inq_varid (ncid,vname(1,indxShflx)(1:lvar),
     &                                                   vidShflx)
          if (ierr .ne. nf_noerr) then
            write(stdout,1) vname(1,indxShflx)(1:lvar), ncname(1:lstr)
            goto 99                                       !--> ERROR
          endif
#  ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidShflx,nf_collective)
#  endif
        endif
#  ifdef BHFLUX
!
! Bhflx     Bottom net heat flux [W/m2]
!
        if (wrt(indxBhflx)) then
          lvar=lenstr(vname(1,indxBhflx))
          ierr=nf_inq_varid (ncid,vname(1,indxBhflx)(1:lvar),
     &                                                   vidBhflx)
          if (ierr .ne. nf_noerr) then
          write(stdout,1) vname(1,indxBhflx)(1:lvar), ncname(1:lstr)
            goto 99                                       !--> ERROR
          endif
#   ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidBhflx,nf_collective)
#   endif
        endif
#  endif
# endif

# ifdef SALINITY
!
! Swflx    E-P flux [cm/day]
!
        if (wrt(indxSwflx)) then
          lvar=lenstr(vname(1,indxSwflx))
          ierr=nf_inq_varid (ncid,vname(1,indxSwflx)(1:lvar),
     &                                                   vidSwflx)
          if (ierr .ne. nf_noerr) then
            write(stdout,1) vname(1,indxSwflx)(1:lvar), ncname(1:lstr)
            goto 99                                       !--> ERROR
          endif
#  ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidSwflx,nf_collective)
#  endif
        endif
#  ifdef BWFLUX
!
! Bwflx    Freshwater at bottom flux [cm/day]
!
        if (wrt(indxBwflx)) then
          lvar=lenstr(vname(1,indxBwflx))
          ierr=nf_inq_varid (ncid,vname(1,indxBwflx)(1:lvar),
     &                                                   vidBwflx)
          if (ierr .ne. nf_noerr) then
          write(stdout,1) vname(1,indxBwflx)(1:lvar), ncname(1:lstr)
            goto 99                                       !--> ERROR
          endif
#   ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidBwflx,nf_collective)
#   endif
        endif
#  endif
# endif
!
! Shflx     Surface shortwave radiation [W/m2]
!
# ifdef TEMPERATURE
        if (wrt(indxShflx_rsw)) then
          lvar=lenstr(vname(1,indxShflx_rsw))
          ierr=nf_inq_varid (ncid,vname(1,indxShflx_rsw)(1:lvar),
     &                                                vidShflx_rsw)
          if (ierr .ne. nf_noerr) then
            write(stdout,1) vname(1,indxShflx_rsw)(1:lvar), ncname(1:lstr)
            goto 99                                       !--> ERROR
          endif
#  ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidShflx_rsw,nf_collective)
#  endif
        endif
# endif
# if defined BULK_FLUX && defined TEMPERATURE
!
        if (wrt(indxShflx_rlw)) then
          lvar=lenstr(vname(1,indxShflx_rlw))
          ierr=nf_inq_varid (ncid,vname(1,indxShflx_rlw)(1:lvar),
     &                                               vidShflx_rlw)
          if (ierr .ne. nf_noerr) then
            write(stdout,1) vname(1,indxShflx_rlw)(1:lvar), ncname(1:lstr)
            goto 99                                       !--> ERROR
          endif
#  ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidShflx_rlw,nf_collective)
#  endif
        endif
!
        if (wrt(indxShflx_lat)) then
          lvar=lenstr(vname(1,indxShflx_lat))
          ierr=nf_inq_varid (ncid,vname(1,indxShflx_lat)(1:lvar),
     &                                               vidShflx_lat)
          if (ierr .ne. nf_noerr) then
            write(stdout,1) vname(1,indxShflx_lat)(1:lvar), ncname(1:lstr)
            goto 99                                       !--> ERROR
          endif
#  ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidShflx_lat,nf_collective)
#  endif
        endif
!
        if (wrt(indxShflx_sen)) then
          lvar=lenstr(vname(1,indxShflx_sen))
          ierr=nf_inq_varid (ncid,vname(1,indxShflx_sen)(1:lvar),
     &                                               vidShflx_sen)
          if (ierr .ne. nf_noerr) then
            write(stdout,1) vname(1,indxShflx_sen)(1:lvar), ncname(1:lstr)
            goto 99                                       !--> ERROR
          endif
#  ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidShflx_sen,nf_collective)
#  endif
        endif
!
# endif
!
# if defined BIOLOGY & !defined PISCES
!
! Depth of the euphotic layer.
!
        if (wrt(indxHel)) then
          lvar=lenstr(vname(1,indxHel))
          ierr=nf_inq_varid (ncid,vname(1,indxHel)(1:lvar), vidHel)
          if (ierr .ne. nf_noerr) then
            write(stdout,1) vname(1,indxHel)(1:lvar), ncname(1:lstr)
            goto 99                                       !--> ERROR
          endif
#  ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidHel,nf_collective)
#  endif
        endif
!
! Chlorophyll/Carbon ratio.
!
#  ifdef BIO_NChlPZD
        if (wrt(indxChC)) then
          lvar=lenstr(vname(1,indxChC))
          ierr=nf_inq_varid (ncid,vname(1,indxChC)(1:lvar), vidChC)
          if (ierr .ne. nf_noerr) then
            write(stdout,1) vname(1,indxChC)(1:lvar), ncname(1:lstr)
            goto 99                                       !--> ERROR
          endif
#   ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidChC,nf_collective)
#   endif
        endif
#   ifdef OXYGEN
!
! Wind speed at 10 m
        if (wrt(indxU10)) then
          lvar=lenstr(vname(1,indxU10))
          ierr=nf_inq_varid (ncid,vname(1,indxU10)(1:lvar), vidU10)
          if (ierr .ne. nf_noerr) then
            write(stdout,1) vname(1,indxU10)(1:lvar), ncname(1:lstr)
            goto 99                                       !--> ERROR
          endif
#    ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidU10,nf_collective)
#    endif
        endif
!
! Gas exchange coefficient of O2
        if (wrt(indxKvO2)) then
          lvar=lenstr(vname(1,indxKvO2))
          ierr=nf_inq_varid (ncid,vname(1,indxKvO2)(1:lvar), vidKvO2)
          if (ierr .ne. nf_noerr) then
            write(stdout,1) vname(1,indxKvO2)(1:lvar), ncname(1:lstr)
            goto 99                                       !--> ERROR
          endif
#    ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidKvO2,nf_collective)
#    endif
        endif
!
! Saturation concentration of O2
        if (wrt(indxO2sat)) then
          lvar=lenstr(vname(1,indxO2sat))
          ierr=nf_inq_varid (ncid,vname(1,indxO2sat)(1:lvar), vidO2sat)
          if (ierr .ne. nf_noerr) then
            write(stdout,1) vname(1,indxO2sat)(1:lvar), ncname(1:lstr)
            goto 99                                       !--> ERROR
          endif
#    ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidO2sat,nf_collective)
#    endif
        endif
#   endif /* OXYGEN */
#  endif
# endif
!
!-----------------------------------------------------------------------
!  USGS sediment model output
!-----------------------------------------------------------------------
!
# ifdef SEDIMENT
!
! Sediment bed layer thickness, porosity, volume fraction of
! size class ised in sediment layer, bed fluxes (2+5*NST r2dgrd variables)
!
        do itrc=1,3+  NST
#  ifdef SUSPLOAD
     &             +2*NST
#  endif
#  ifdef BEDLOAD
     &             +2*NST
#  endif
#  if defined MIXED_BED || defined COHESIVE_BED
     &             +1
#  endif
          indxWrk=indxSed+itrc-2
          if (wrt(indxWrk)) then
            lvar=lenstr(vname(1,indxWrk))
            ierr=nf_inq_varid (ncid,vname(1,indxWrk)(1:lvar),
     &                                                vidSed(itrc))
            if (ierr .ne. nf_noerr) then
              write(stdout,1) vname(1,indxWrk)(1:lvar), ncname(1:lstr)
              goto 99                                       !--> ERROR
            endif
          endif
#  ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidSed(itrc),nf_collective)
#  endif
        enddo
# endif /* SEDIMENT */
# ifdef BBL
!
! Abed       Wind-induced, bed wave excursion amplitude (m)
! Hripple    Bed ripple height (m)
! Lripple    Bed ripple length (m)
! Zbnot      Physical hydraulic bottom roughness (grain roughness) (m)
! Zbapp      Effective apparent hydraulic bottom roughness
! bostrw     Wave-induced kinematic bottom stress (N/m2)
!
        do itrc=1,6
          indxWrk=indxBBL+itrc-1
          if (wrt(indxWrk)) then
            lvar=lenstr(vname(1,indxWrk))
            ierr=nf_inq_varid (ncid,vname(1,indxWrk)(1:lvar),
     &                                                vidBBL(itrc))
            if (ierr .ne. nf_noerr) then
              write(stdout,1) vname(1,indxWrk)(1:lvar), ncname(1:lstr)
              goto 99                                       !--> ERROR
            endif
#  ifdef NC4PAR
            ierr=nf_var_par_access(ncid,vidBBL(itrc),nf_collective)
#  endif
          endif
        enddo
# endif /* BBL */
# if defined SST_SKIN && defined TEMPERATURE
!
! sst_skin     surface skin temperature (degC)
!
       if (wrt(indxT)) then
         lvar=lenstr(vname(1,indxSST_skin))
         ierr=nf_inq_varid (ncid,vname(1,indxSST_skin)(1:lvar),
     &                                               vidSST_skin)
         if (ierr .ne. nf_noerr) then
           write(stdout,1) vname(1,indxSST_skin)(1:lvar), ncname(1:lstr)
           goto 99                                       !--> ERROR
         endif
#  ifdef NC4PAR
         ierr=nf_var_par_access(ncid,vidSST_skin,nf_collective)
#  endif
       endif
# endif /* SST_SKIN */
!
# ifdef MUSTANG
       do indx=1,outMust_nbvar
          if (outMust(indx)) then 
             lvar=lenstr(vname_Must(1,indx))
             ierr=nf_inq_varid (ncid, vname_Must(1,indx)(1:lvar),
     &                             vidMust(indx))
#  ifdef NC4PAR
             ierr=nf_var_par_access(ncid,vidMust(indx),
     &                             nf_collective)
#  endif
          endif
      enddo ! indx
# endif /* MUSTANG */

# ifdef OBSTRUCTION
      do itrv=1,17*obst_nbvar+15
         indx=itrv
         if (outObst(indx)) then 
            lvar=lenstr(vname_obst(1,indx))
            ierr=nf_inq_varid (ncid, vname_obst(1,indx)(1:lvar),
     &                             vidObst(itrv))
#  ifdef NC4PAR
            ierr=nf_var_par_access(ncid,vidObst(itrv),
     &                             nf_collective)
#  endif
         endif
      enddo ! itrv

# endif /* OBSTRUCTION */
!
#endif /* SOLVE3D */
!
#ifdef WAVE_IO
!
! Wave action density, wavenumber in XI/ETA directions
!
# ifndef WAVE_ROLLER
        do iwkb=1,7
# else
        do iwkb=1,9
# endif
          if (wrt(indxHRM+iwkb-1)) then
            lvar=lenstr(vname(1,indxHRM+iwkb-1))
            ierr=nf_inq_varid (ncid, vname(1,indxHRM+iwkb-1)(1:lvar),
     &                                                  vidWAVE(iwkb))
            if (ierr .ne. nf_noerr) then
              write(stdout,1) vname(1,indxHRM+iwkb-1)(1:lvar),
     &                                                 ncname(1:lstr)
              goto 99                                     !--> ERROR
            endif
# ifdef NC4PAR
            ierr=nf_var_par_access(ncid,vidWAVE(iwkb),nf_collective)
# endif
          endif
        enddo
#endif /* WAVE_IO */

#ifdef MRL_WCI
!
! sup: quasi-static sea-level response (wave set-up/down)
!
        if (wrt(indxSUP)) then
          lvar=lenstr(vname(1,indxSUP))
          ierr=nf_inq_varid (ncid, vname(1,indxSUP)(1:lvar), vidSUP)
          if (ierr .ne. nf_noerr) then
            write(stdout,1) vname(1,indxSUP)(1:lvar), ncname(1:lstr)
            goto 99                                       !--> ERROR
          endif
# ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidSUP,nf_collective)
# endif
        endif
!
! ust2d & vst2d: 2D, depth-averaged Stokes drift velocities.
!
        if (wrt(indxUST2D)) then
          lvar=lenstr(vname(1,indxUST2D))
          ierr=nf_inq_varid (ncid, vname(1,indxUST2D)(1:lvar), vidUST2D)
          if (ierr .ne. nf_noerr) then
            write(stdout,1) vname(1,indxUST2D)(1:lvar), ncname(1:lstr)
            goto 99                                       !--> ERROR
          endif
# ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidUST2D,nf_collective)
# endif
        endif
        if (wrt(indxVST2D)) then
          lvar=lenstr(vname(1,indxVST2D))
          ierr=nf_inq_varid (ncid, vname(1,indxVST2D)(1:lvar), vidVST2D)
          if (ierr .ne. nf_noerr) then
            write(stdout,1) vname(1,indxVST2D)(1:lvar), ncname(1:lstr)
            goto 99                                       !--> ERROR
          endif
# ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidVST2D,nf_collective)
# endif
        endif
# ifdef SOLVE3D
!
! ust3d & vst3d: 3D Stokes drift velocities.
!
        if (wrt(indxUST)) then
          lvar=lenstr(vname(1,indxUST))
          ierr=nf_inq_varid (ncid, vname(1,indxUST)(1:lvar), vidUST)
          if (ierr .ne. nf_noerr) then
            write(stdout,1) vname(1,indxUST)(1:lvar), ncname(1:lstr)
            goto 99                                       !--> ERROR
          endif
#  ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidUST,nf_collective)
#  endif
        endif

        if (wrt(indxVST)) then
          lvar=lenstr(vname(1,indxVST))
          ierr=nf_inq_varid (ncid, vname(1,indxVST)(1:lvar), vidVST)
          if (ierr .ne. nf_noerr) then
            write(stdout,1) vname(1,indxVST)(1:lvar), ncname(1:lstr)
            goto 99                                       !--> ERROR
          endif
#  ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidVST,nf_collective)
#  endif
        endif
!
! wst: vertical Stokes drift velocity (m/s) at rho-point
!
        if (wrt(indxWST)) then
          lvar=lenstr(vname(1,indxWST))
          ierr=nf_inq_varid (ncid, vname(1,indxWST)(1:lvar), vidWST)
          if (ierr .ne. nf_noerr) then
            write(stdout,1) vname(1,indxWST)(1:lvar), ncname(1:lstr)
            goto 99                                       !--> ERROR
          endif
#  ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidWST,nf_collective)
#  endif
        endif
!
! Akb: Vertical eddy viscosity coefficient due to depth-induced wave breaking.
!
        if (wrt(indxAkb)) then
          lvar=lenstr(vname(1,indxAkb))
          ierr=nf_inq_varid (ncid, vname(1,indxAkb)(1:lvar), vidAkb)
          if (ierr .ne. nf_noerr) then
            write(stdout,1) vname(1,indxAkb)(1:lvar), ncname(1:lstr)
            goto 99                                       !--> ERROR
          endif
#  ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidAkb,nf_collective)
#  endif
        endif
!
! Akw: Vertical eddy diffusivity coefficient due to primary waves.
!
        if (wrt(indxAkw)) then
          lvar=lenstr(vname(1,indxAkw))
          ierr=nf_inq_varid (ncid, vname(1,indxAkw)(1:lvar), vidAkw)
          if (ierr .ne. nf_noerr) then
            write(stdout,1) vname(1,indxAkw)(1:lvar), ncname(1:lstr)
            goto 99                                       !--> ERROR
          endif
#  ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidAkw,nf_collective)
#  endif
        endif
!
! kvf: vertical vortex force term (u^St du/dz)
!
        if (wrt(indxKVF)) then
          lvar=lenstr(vname(1,indxKVF))
          ierr=nf_inq_varid (ncid, vname(1,indxKVF)(1:lvar), vidKVF)
          if (ierr .ne. nf_noerr) then
            write(stdout,1) vname(1,indxKVF)(1:lvar), ncname(1:lstr)
            goto 99                                       !--> ERROR
          endif
#  ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidKVF,nf_collective)
#  endif
        endif
!
! calP: surface pressure correction term appeared in prsgrd term
!
        if (wrt(indxCALP)) then
          lvar=lenstr(vname(1,indxCALP))
          ierr=nf_inq_varid (ncid, vname(1,indxCALP)(1:lvar), vidCALP)
          if (ierr .ne. nf_noerr) then
            write(stdout,1) vname(1,indxCALP)(1:lvar), ncname(1:lstr)
            goto 99                                       !--> ERROR
          endif
#  ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidCALP,nf_collective)
#  endif
        endif
!
! Kapsrf: surface Bernoulli head term appeared in prsgrd term
!
        if (wrt(indxKAPS)) then
          lvar=lenstr(vname(1,indxKAPS))
          ierr=nf_inq_varid (ncid, vname(1,indxKAPS)(1:lvar), vidKAPS)
          if (ierr .ne. nf_noerr) then
            write(stdout,1) vname(1,indxKAPS)(1:lvar), ncname(1:lstr)
            goto 99                                       !--> ERROR
          endif
#  ifdef NC4PAR
          ierr=nf_var_par_access(ncid,vidKAPS,nf_collective)
#  endif
        endif
# endif  /* SOLVE3D */
#endif   /* MRL_WCI */

      MPI_master_only write(*,'(6x,2A,i4,1x,A,i4)')
     &                     'DEF_HIS/AVG -- Opened ',
     &                     'existing file  from record =', rec
     &                      MYID

#if defined MPI & !defined PARALLEL_FILES & !defined NC4PAR
      else
        ierr=nf_open (ncname(1:lstr), nf_write, ncid)
        if (ierr .ne. nf_noerr) then
          MPI_master_only write(stdout,'(/1x,4A,2x,A,I4/)')
     &                'DEF_HIS/AVG ERROR: ',
     &                'Cannot open file ''', ncname(1:lstr), '''.'
     &                 MYID
          goto 99                                         !--> ERROR
        endif
#endif
      endif             !<--  create_new_file

      ierr=nf_set_fill (ncid, nf_nofill, lvar)
      if (ierr .ne. nf_noerr) then
        write(*,'(6x,2A,i4,1x,A,i4)') 'DEF_HIS/AVG ERROR: Cannot ',
     &    'switch to ''nf_nofill'' more; netCDF error code =', ierr
      endif
   1  format(/1x,'DEF_HIS/AVG ERROR: Cannot find variable ''',
     &                   A, ''' in netCDF file ''', A, '''.'/)
#if (defined PUT_GRID_INTO_HISTORY && !defined AVRH)\
 || (defined PUT_GRID_INTO_AVERAGES && defined AVRH)
!
! Write grid variables.
! ----- ---- ----------
!
        if (total_rec.le.1) call wrt_grid (ncid, ncname, lstr)
#endif

  99  return
      end

#undef ncname
#undef rec_per_file
#undef wrt
#undef vidTime
#undef vidTime2
#undef vidTstep
#undef vidZ
#undef vidUb
#undef vidVb
#undef vidU
#undef vidV
#undef vidT
#undef vidR
#undef vidbvf
#undef vidO
#undef vidW
#undef vidW
#undef vidVisc
#undef vidDiff
#undef vidAkv
#undef vidAkt
#undef vidAks
#undef vidHbl
#undef vidHbbl
#undef vidTke
#undef vidGls
#undef vidLsc
#undef vidHel
#undef vidChC
#undef vidU10
#undef vidKvO2
#undef vidO2sat
#undef vidAOU
#undef vidwind10
#undef vidSed
#undef vidBBL
#undef vidBostr
#undef vidBustr
#undef vidBvstr
#undef vidWstr
#undef vidUWstr
#undef vidVWstr
#undef vidShflx
#undef vidSwflx
#undef vidShflx_rsw
#if defined BULK_FLUX
# undef vidShflx_rlw
# undef vidShflx_lat
# undef vidShflx_sen
#endif
#if defined BHFLUX && defined TEMPERATURE
# undef vidBhflx
#endif
#if defined BWFLUX & defined SALINITY
# undef vidBwflx
#endif
#undef vidSST_skin
#undef vidWAVE
#undef vidSUP
#undef vidUST2D
#undef vidVST2D
#undef vidUST
#undef vidVST
#undef vidWST
#undef vidAkb
#undef vidAkw
#undef vidKVF
#undef vidCALP
#undef vidKAPS
#undef vidHm
#undef vidMust
#undef vidObst

#ifdef AVERAGES
# ifndef AVRH
#  define AVRH
#  include "def_his.F"
# endif
#endif
