!======================================================================
! 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
!======================================================================
!
#include "cppdefs.h"
#ifdef STATIONS
!
      subroutine def_sta(ncid, total_rec, ierr,res)
!
!====================================================================
!                                                                    !
!  This routine creates STATIONS data NetCDF file:                   !
!  Define its variables, dimensions,and attributes                   !
!                                                                    !
!=====================================================================
!
          implicit none
# include "param.h"
# include "mixing.h"
# include "ncscrum.h"
# include "nc_sta.h"
# include "sta.h"
# include "scalars.h"
# include "strings.h"
# include "netcdf.inc"
!
      logical create_new_file, res
      integer ncid, total_rec, rec, ierr
      integer lstr, lvar, stadim, ftimedim, twodim
      integer pgrd(2), temp(2), checkdims
      integer lenstr
# ifdef MUSTANG
      integer itrc
# endif
# ifdef SOLVE3D
      integer trcdim
#  ifdef ALL_SIGMA
      integer srho, allSigma(3)
#  endif
# endif /* SOLVE3D */
      character*65 vinfo(4)
      if (may_day_flag.ne.0) return      !-->  EXIT

# define rec_per_file nrpfsta

      ierr=0
      lstr=lenstr(staname)
      if (rec_per_file.gt.0) then
        lvar=total_rec-(1+mod(total_rec-1, nrpfsta))
        call insert_time_index (staname, lstr, lvar, ierr)
        if (ierr .ne. 0) goto 99
      endif

!  Create a new station data file.
!---------------------------------
!
      create_new_file=ldefsta
      if (ncid.ne.-1) create_new_file=.false.
# if defined MPI & !defined PARALLEL_FILES
      if (mynode.gt.0) create_new_file=.false.
# endif

!
!=================================================================
!  Create new station file
!=================================================================
!
 10   if (create_new_file) then

        lstr=lenstr(staname)
        ierr=nf_create(staname(1:lstr),NF_CLOBBER,ncid)
        if (ierr.ne.nf_noerr) then
          write(stdout,11) staname(1:lstr)
          may_day_flag=3
          return                         !-->  EXIT
        endif

!
!  Put global attributes.
!  ----------------------
!
        call put_global_atts (ncid, ierr)
        if (ierr.ne.nf_noerr) then
          write(stdout,11) staname(1:lstr)
          may_day_flag=3
          return                         !-->  EXIT
        endif

!
!  Define the dimensions of staggered fields.
!  ------------------------------------------
!
        ierr=nf_def_dim(ncid,'stanum',nstas0,stadim)
# ifdef SOLVE3D
        ierr=nf_def_dim(ncid,'tracer',NT,trcdim)
#  ifdef ALL_SIGMA
        ierr=nf_def_dim(ncid,'s_rho',N,srho)
#  endif
# endif /* SOLVE3D */
        ierr=nf_def_dim(ncid,'ftime',nf_unlimited,ftimedim)
        ierr=nf_def_dim(ncid,'two',2,twodim)
!
!  Define dimension vectors for point variables.
!  ---------------------------------------------
!
        pgrd(1)=stadim
        pgrd(2)=ftimedim
!
        temp(1)=twodim
        temp(2)=ftimedim
!
# if defined SOLVE3D && defined ALL_SIGMA
        allSigma(1)=srho
        allSigma(2)=stadim
        allSigma(3)=ftimedim
# endif

!  Define variables and their attributes.
!  --------------------------------------
!
!  Define time step and model time.

        vinfo(1)='time_step'
        vinfo(2)='time step and record numbers from initialization'
        vinfo(3)='nondimensionnal'
        vinfo(4)='time step/record number, vector, series'
        lvar=lenstr(vinfo(1))
        ierr=nf_def_var(ncid,vinfo(1)(1:lvar),nf_int,
     &                    2,temp,staTstep)
        lvar=lenstr(vinfo(2))
        ierr=nf_put_att_text(ncid,staTstep,'long_name',lvar,
     &                         vinfo(2)(1:lvar))
        lvar=lenstr(vinfo(3))
        ierr=nf_put_att_text(ncid,staTstep,'units',lvar,
     &                         vinfo(3)(1:lvar))
        lvar=lenstr(vinfo(4))
        ierr=nf_put_att_text(ncid,staTstep,'field',lvar,
     &                         vinfo(4)(1:lvar))

!
!  Define time
!
        lvar=lenstr(vname(1,indxTime))
        ierr=nf_def_var(ncid,vname(1,indxTime)(1:lvar),NF_FTYPE,
     &                    1,ftimedim,staTime)
        lvar=lenstr(vname(2,indxTime))
        ierr=nf_put_att_text(ncid,staTime,'long_name',lvar,
     &                         vname(2,indxTime)(1:lvar))
        lvar=lenstr(vname(3,indxTime))
        ierr=nf_put_att_text(ncid,staTime,'units',lvar,
     &                         vname(3,indxTime)(1:lvar))
        lvar=lenstr(vname(4,indxTime))
        ierr=nf_put_att_text(ncid,staTime,'field',lvar,
     &                         vname(4,indxTime)(1:lvar))
!
!  Define stations (lon,lat) or (x,y) and z locations.
!
       if (wrtsta(indxstaGrd)) then
# ifdef AGRIF
        vinfo(1)='grid_level'
        vinfo(2)='grid level in nested grid hierarchy'
        vinfo(3)='nondimensionnal'
        lvar=lenstr(vinfo(1))
        ierr=nf_def_var(ncid,vinfo(1)(1:lvar),nf_int,
     &                    1,stadim,staGlevel)
        lvar=lenstr(vinfo(2))
        ierr=nf_put_att_text(ncid,staGlevel,'long_name',lvar,
     &                         vinfo(2)(1:lvar))
        lvar=lenstr(vinfo(3))
        ierr=nf_put_att_text(ncid,staGlevel,'units',lvar,
     &                         vinfo(3)(1:lvar))
# endif /* AGRIF */
# ifdef SPHERICAL
        vinfo(1)='lon'
        vinfo(2)='longitude of model control points'
        vinfo(3)='degree_east'
        vinfo(4)='lon, scalar, series'
        lvar=lenstr(vinfo(1))
        ierr=nf_def_var(ncid,vinfo(1)(1:lvar),NF_FTYPE,
     &                    1,stadim,staLon)
        lvar=lenstr(vinfo(2))
        ierr=nf_put_att_text(ncid,staLon,'long_name',lvar,
     &                         vinfo(2)(1:lvar))
        lvar=lenstr(vinfo(3))
        ierr=nf_put_att_text(ncid,staLon,'units',lvar,
     &                         vinfo(3)(1:lvar))
        lvar=lenstr(vinfo(4))
        ierr=nf_put_att_text(ncid,staLon,'field',lvar,
     &                         vinfo(4)(1:lvar))
        vinfo(1)='lat'
        vinfo(2)='latitude of model control points'
        vinfo(3)='degree_north'
        vinfo(4)='lat, scalar, series'
        lvar=lenstr(vinfo(1))
        ierr=nf_def_var(ncid,vinfo(1)(1:lvar),NF_FTYPE,
     &                    1,stadim,staLat)
        lvar=lenstr(vinfo(2))
        ierr=nf_put_att_text(ncid,staLat,'long_name',lvar,
     &                         vinfo(2)(1:lvar))
        lvar=lenstr(vinfo(3))
        ierr=nf_put_att_text(ncid,staLat,'units',lvar,
     &                         vinfo(3)(1:lvar))
        lvar=lenstr(vinfo(4))
        ierr=nf_put_att_text(ncid,staLat,'field',lvar,
     &                         vinfo(4)(1:lvar))
# endif  /* SPHERICAL */
!
! Define Station X-position in the grid
!
        vinfo(1)='Xgrid'
        vinfo(2)='x-grid sta locations'
        vinfo(3)='nondimensional'
        vinfo(4)='Xgrid, scalar, series'
        lvar=lenstr(vinfo(1))
        ierr=nf_def_var(ncid,vinfo(1)(1:lvar),NF_FTYPE,
     &                    1,stadim,staXgrd)
        lvar=lenstr(vinfo(2))
        ierr=nf_put_att_text(ncid,staXgrd,'long_name',lvar,
     &                         vinfo(2)(1:lvar))
        lvar=lenstr(vinfo(3))
        ierr=nf_put_att_text(ncid,staXgrd,'units',lvar,
     &                         vinfo(3)(1:lvar))
        lvar=lenstr(vinfo(4))
        ierr=nf_put_att_text(ncid,staXgrd,'field',lvar,
     &                         vinfo(4)(1:lvar))
!
! Define float Y-position in the grid
!
        vinfo(1)='Ygrid'
        vinfo(2)='y-grid sta locations'
        vinfo(3)='nondimensional'
        vinfo(4)='Ygrid, scalar, series'
        lvar=lenstr(vinfo(1))
        ierr=nf_def_var(ncid,vinfo(1)(1:lvar),NF_FTYPE,
     &                    1,stadim,staYgrd)
        lvar=lenstr(vinfo(2))
        ierr=nf_put_att_text(ncid,staYgrd,'long_name',lvar,
     &                         vinfo(2)(1:lvar))
        lvar=lenstr(vinfo(3))
        ierr=nf_put_att_text(ncid,staYgrd,'units',lvar,
     &                         vinfo(3)(1:lvar))
        lvar=lenstr(vinfo(4))
        ierr=nf_put_att_text(ncid,staYgrd,'field',lvar,
     &                         vinfo(4)(1:lvar))
# ifdef SOLVE3D
!
!  Define depth.
!
#  ifdef ALL_SIGMA
        vinfo(1)='depth'
        vinfo(2)='depth of stations levels'
        vinfo(3)='meter'
        lvar=lenstr(vinfo(1))
        ierr=nf_def_var(ncid,vinfo(1)(1:lvar),NF_FTYPE,
     &                    3,allSigma,staDepth)
        lvar=lenstr(vinfo(2))
        ierr=nf_put_att_text(ncid,staDepth,'long_name',lvar,
     &                         vinfo(2)(1:lvar))
        lvar=lenstr(vinfo(3))
        ierr=nf_put_att_text(ncid,staDepth,'units',lvar,
     &                         vinfo(3)(1:lvar))
#  else
        vinfo(1)='depth'
        vinfo(2)='depth of point stations'
        vinfo(3)='meter'
        lvar=lenstr(vinfo(1))
        ierr=nf_def_var(ncid,vinfo(1)(1:lvar),NF_FTYPE,
     &                    1,pgrd,staDepth)
        lvar=lenstr(vinfo(2))
        ierr=nf_put_att_text(ncid,staDepth,'long_name',lvar,
     &                         vinfo(2)(1:lvar))
        lvar=lenstr(vinfo(3))
        ierr=nf_put_att_text(ncid,staDepth,'units',lvar,
     &                         vinfo(3)(1:lvar))
#  endif
!
! Define  vertical positions S or depth(m)
!
#  ifdef ALL_SIGMA
!     Do nothing in this case Zgrid is allways 1:N
#  else
        vinfo(1)='Zgrid'
        vinfo(2)='time varing sigma level of point stations'
        vinfo(3)='nondimensional'
        vinfo(4)='Zgrid, scalar, series'
        lvar=lenstr(vinfo(1))
        ierr=nf_def_var(ncid,vinfo(1)(1:lvar),NF_FTYPE,
     &                    2,pgrd,staZgrd)
        lvar=lenstr(vinfo(2))
        ierr=nf_put_att_text(ncid,staZgrd,'long_name',lvar,
     &                         vinfo(2)(1:lvar))
        lvar=lenstr(vinfo(3))
        ierr=nf_put_att_text(ncid,staZgrd,'units',lvar,
     &                         vinfo(3)(1:lvar))
        lvar=lenstr(vinfo(4))
        ierr=nf_put_att_text(ncid,staZgrd,'field',lvar,
     &                         vinfo(4)(1:lvar))
#  endif
# endif /* SOLVE3D */
       endif ! grid parameters
!
!  --- Define variables for station data ---
!
!
!  Define sea level
!
       vinfo(1)='zeta'
       vinfo(2)='sea level at station'
       vinfo(3)='meter'
       lvar=lenstr(vinfo(1))
       ierr=nf_def_var(ncid,vinfo(1)(1:lvar),NF_FTYPE,
     &                 2,pgrd,staZeta)
       lvar=lenstr(vinfo(2))
       ierr=nf_put_att_text(ncid,staZeta,'long_name',lvar,
     &                      vinfo(2)(1:lvar))
       lvar=lenstr(vinfo(3))
       ierr=nf_put_att_text(ncid,staZeta,'units',lvar,
     &                      vinfo(3)(1:lvar))
# ifdef SOLVE3D
!
!  Define temperature
!
       if (wrtsta(indxstaTemp)) then
#  ifdef ALL_SIGMA
        vinfo(1)='temp'
        vinfo(2)='temperature'
        vinfo(3)='degrees Celsius'
        lvar=lenstr(vinfo(1))
        ierr=nf_def_var(ncid,vinfo(1)(1:lvar),NF_FTYPE,
     &                    3,AllSigma,staTemp)
        lvar=lenstr(vinfo(2))
        ierr=nf_put_att_text(ncid,staTemp,'long_name',lvar,
     &                         vinfo(2)(1:lvar))
        lvar=lenstr(vinfo(3))
        ierr=nf_put_att_text(ncid,staTemp,'units',lvar,
     &                         vinfo(3)(1:lvar))
#  else
        vinfo(1)='temp'
        vinfo(2)='temperature'
        vinfo(3)='degrees Celsius'
        lvar=lenstr(vinfo(1))
        ierr=nf_def_var(ncid,vinfo(1)(1:lvar),NF_FTYPE,
     &                    2,pgrd,staTemp)
        lvar=lenstr(vinfo(2))
        ierr=nf_put_att_text(ncid,staTemp,'long_name',lvar,
     &                         vinfo(2)(1:lvar))
        lvar=lenstr(vinfo(3))
        ierr=nf_put_att_text(ncid,staTemp,'units',lvar,
     &                         vinfo(3)(1:lvar))
#  endif
       endif ! temperature
!
!  Define salinity
!
#  ifdef SALINITY
       if (wrtsta(indxstaSalt)) then
#   ifdef ALL_SIGMA
        vinfo(1)='salt'
        vinfo(2)='salinity'
        vinfo(3)='PSU'
        lvar=lenstr(vinfo(1))
        ierr=nf_def_var(ncid,vinfo(1)(1:lvar),NF_FTYPE,
     &                    3,allSigma,staSal)
        lvar=lenstr(vinfo(2))
        ierr=nf_put_att_text(ncid,staSal,'long_name',lvar,
     &                         vinfo(2)(1:lvar))
        lvar=lenstr(vinfo(3))
        ierr=nf_put_att_text(ncid,staSal,'units',lvar,
     &                         vinfo(3)(1:lvar))
#   else
        vinfo(1)='salt'
        vinfo(2)='salinity'
        vinfo(3)='PSU'
        lvar=lenstr(vinfo(1))
        ierr=nf_def_var(ncid,vinfo(1)(1:lvar),NF_FTYPE,
     &                    2,pgrd,staSal)
        lvar=lenstr(vinfo(2))
        ierr=nf_put_att_text(ncid,staSal,'long_name',lvar,
     &                         vinfo(2)(1:lvar))
        lvar=lenstr(vinfo(3))
        ierr=nf_put_att_text(ncid,staSal,'units',lvar,
     &                         vinfo(3)(1:lvar))
#   endif
       endif ! salinity
#  endif  /* SALINITY */

#  ifdef MUSTANG
       do itrc=1,NT-2
#   ifdef ALL_SIGMA
        vinfo(1)=vname(1,indxT+2+itrc-1)
        vinfo(2)='concentration'
        vinfo(3)='kg/m3'
        lvar=lenstr(vinfo(1))
        ierr=nf_def_var(ncid,vinfo(1)(1:lvar),NF_FTYPE,
     &                    3,allSigma,staMUS(itrc))
        lvar=lenstr(vinfo(2))
        ierr=nf_put_att_text(ncid,staMUS(itrc),'long_name',lvar,
     &                         vinfo(2)(1:lvar))
        lvar=lenstr(vinfo(3))
        ierr=nf_put_att_text(ncid,staMUS(itrc),'units',lvar,
     &                         vinfo(3)(1:lvar))
#   else
        vinfo(1)=vname(1,indxT+2+itrc-1)
        vinfo(2)='concentration'
        vinfo(3)='kg/m3'
        lvar=lenstr(vinfo(1))
        ierr=nf_def_var(ncid,vinfo(1)(1:lvar),NF_FTYPE,
     &                    2,pgrd,staMUS(itrc))
        lvar=lenstr(vinfo(2))
        ierr=nf_put_att_text(ncid,staMUS(itrc),'long_name',lvar,
     &                         vinfo(2)(1:lvar))
        lvar=lenstr(vinfo(3))
        ierr=nf_put_att_text(ncid,staMUS(itrc),'units',lvar,
     &                         vinfo(3)(1:lvar))
#   endif
       enddo
#  endif  /* MUSTANG*/


!
!  Define density anomaly.
!
       if (wrtsta(indxstaRho)) then
#  ifdef ALL_SIGMA
        vinfo(1)='rho'
        vinfo(2)='density anomaly'
        vinfo(3)='kilogram meter-3'
        lvar=lenstr(vinfo(1))
        ierr=nf_def_var(ncid,vinfo(1)(1:lvar),NF_FTYPE,
     &                    3,allSigma,staDen)
        lvar=lenstr(vinfo(2))
        ierr=nf_put_att_text(ncid,staDen,'long_name',lvar,
     &                         vinfo(2)(1:lvar))
        lvar=lenstr(vinfo(3))
        ierr=nf_put_att_text(ncid,staDen,'units',lvar,
     &                         vinfo(3)(1:lvar))
#  else
        vinfo(1)='rho'
        vinfo(2)='density anomaly'
        vinfo(3)='kilogram meter-3'
        lvar=lenstr(vinfo(1))
        ierr=nf_def_var(ncid,vinfo(1)(1:lvar),NF_FTYPE,
     &                    2,pgrd,staDen)
        lvar=lenstr(vinfo(2))
        ierr=nf_put_att_text(ncid,staDen,'long_name',lvar,
     &                         vinfo(2)(1:lvar))
        lvar=lenstr(vinfo(3))
        ierr=nf_put_att_text(ncid,staDen,'units',lvar,
     &                         vinfo(3)(1:lvar))
#  endif
       endif  ! density
!
!  Velocity components
!
       if (wrtsta(indxstaVel)) then
#  ifdef ALL_SIGMA
        vinfo(1)='u'
        vinfo(2)='U component velocity'
        vinfo(3)='meter/s'
        lvar=lenstr(vinfo(1))
        ierr=nf_def_var(ncid,vinfo(1)(1:lvar),NF_FTYPE,
     &                   3,allSigma,staU)
        lvar=lenstr(vinfo(2))
        ierr=nf_put_att_text(ncid,staU,'long_name',lvar,
     &                         vinfo(2)(1:lvar))
        lvar=lenstr(vinfo(3))
        ierr=nf_put_att_text(ncid,staU,'units',lvar,
     &                         vinfo(3)(1:lvar))
        vinfo(1)='v'
        vinfo(2)='V component velocity'
        vinfo(3)='meter/s'
        lvar=lenstr(vinfo(1))
        ierr=nf_def_var(ncid,vinfo(1)(1:lvar),NF_FTYPE,
     &                    3,allSigma,staV)
        lvar=lenstr(vinfo(2))
        ierr=nf_put_att_text(ncid,staV,'long_name',lvar,
     &                         vinfo(2)(1:lvar))
        lvar=lenstr(vinfo(3))
        ierr=nf_put_att_text(ncid,staV,'units',lvar,
     &                         vinfo(3)(1:lvar))
#  else
        vinfo(1)='u'
        vinfo(2)='U component velocity'
        vinfo(3)='meter/s'
        lvar=lenstr(vinfo(1))
        ierr=nf_def_var(ncid,vinfo(1)(1:lvar),NF_FTYPE,
     &                   2,pgrd,staU)
        lvar=lenstr(vinfo(2))
        ierr=nf_put_att_text(ncid,staU,'long_name',lvar,
     &                         vinfo(2)(1:lvar))
        lvar=lenstr(vinfo(3))
        ierr=nf_put_att_text(ncid,staU,'units',lvar,
     &                         vinfo(3)(1:lvar))
        vinfo(1)='v'
        vinfo(2)='V component velocity'
        vinfo(3)='meter/s'
        lvar=lenstr(vinfo(1))
        ierr=nf_def_var(ncid,vinfo(1)(1:lvar),NF_FTYPE,
     &                    2,pgrd,staV)
        lvar=lenstr(vinfo(2))
        ierr=nf_put_att_text(ncid,staV,'long_name',lvar,
     &                         vinfo(2)(1:lvar))
        lvar=lenstr(vinfo(3))
        ierr=nf_put_att_text(ncid,staV,'units',lvar,
     &                         vinfo(3)(1:lvar))
#  endif
       endif ! velocity
# else
       if (wrtsta(indxstaVel)) then
        vinfo(1)='u'
        vinfo(2)='U barotrope component velocity'
        vinfo(3)='meter/s'
        lvar=lenstr(vinfo(1))
        ierr=nf_def_var(ncid,vinfo(1)(1:lvar),NF_FTYPE,
     &                   2,pgrd,staU)
        lvar=lenstr(vinfo(2))
        ierr=nf_put_att_text(ncid,staU,'long_name',lvar,
     &                         vinfo(2)(1:lvar))
        lvar=lenstr(vinfo(3))
        ierr=nf_put_att_text(ncid,staU,'units',lvar,
     &                         vinfo(3)(1:lvar))
        vinfo(1)='v'
        vinfo(2)='V barotrope component velocity'
        vinfo(3)='meter/s'
        lvar=lenstr(vinfo(1))
        ierr=nf_def_var(ncid,vinfo(1)(1:lvar),NF_FTYPE,
     &                    2,pgrd,staV)
        lvar=lenstr(vinfo(2))
        ierr=nf_put_att_text(ncid,staV,'long_name',lvar,
     &                         vinfo(2)(1:lvar))
        lvar=lenstr(vinfo(3))
        ierr=nf_put_att_text(ncid,staV,'units',lvar,
     &                         vinfo(3)(1:lvar))
       endif
# endif  /* SOLVE3D */


!  Leave definition mode.
!  ----------------------
!
        ierr=nf_enddef(ncid)
        res=.true. ! marker that a file has been created
        write(stdout,'(6x,4A,1x,A,i3)') 'DEF_STA -- Created ',
     &                'new netCDF file ''', staname(1:lstr), '''.'
     &                 MYID


      elseif (ncid.eq.-1) then
!
!===================================================================
! Open an existing file and prepare for appending data.
!===================================================================
!
        ierr=nf_open (staname(1:lstr), nf_write, ncid)

        if (ierr. eq. nf_noerr) then
          ierr=checkdims (ncid, staname, 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_STA WARNING: Actual number of records',
     &               rec,  'in netCDF file',  '''',  staname(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
          if (mynode.eq.0) then
            create_new_file=.true.
            goto 10
          else
            write(stdout,'(/1x,4A,2x,A,I4/)') 'DEF_HIS/AVG ERROR: ',
     &                  'Cannot open file ''', staname(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', staTstep)
        if (ierr .ne. nf_noerr) then
          write(stdout,1) 'time_step', staname(1:lstr)
          goto 99                                         !--> ERROR
        endif
!
! Time.
!
        lvar=lenstr(vname(1,indxTime))
        ierr=nf_inq_varid (ncid,vname(1,indxTime)(1:lvar),staTime)
        if (ierr .ne. nf_noerr) then
          write(stdout,1) vname(1,indxTime)(1:lvar), staname(1:lstr)
          goto 99                                         !--> ERROR
        endif
!
!  Define stations (lon,lat) or (x,y) and z locations.
!
       if (wrtsta(indxstaGrd)) then
# ifdef AGRIF
        vinfo(1)='grid_level'
        lvar=lenstr(vinfo(1))
        ierr=nf_inq_varid(ncid,vinfo(1)(1:lvar),staGlevel)
        if (ierr .ne. nf_noerr) then
          write(stdout,1) vinfo(1)(1:lvar), staname(1:lstr)
          goto 99                                         !--> ERROR
        endif
# endif /* AGRIF */
# ifdef SPHERICAL
        vinfo(1)='lon'
        lvar=lenstr(vinfo(1))
        ierr=nf_inq_varid(ncid,vinfo(1)(1:lvar),staLon)
        if (ierr .ne. nf_noerr) then
          write(stdout,1) vinfo(1)(1:lvar), staname(1:lstr)
          goto 99                                         !--> ERROR
        endif
        vinfo(1)='lat'
        lvar=lenstr(vinfo(1))
        ierr=nf_inq_varid(ncid,vinfo(1)(1:lvar),staLat)
        if (ierr .ne. nf_noerr) then
          write(stdout,1) vinfo(1)(1:lvar), staname(1:lstr)
          goto 99                                         !--> ERROR
        endif
# endif  /* SPHERICAL */
!
! Define Station X-position in the grid
!
        vinfo(1)='Xgrid'
        lvar=lenstr(vinfo(1))
        ierr=nf_inq_varid(ncid,vinfo(1)(1:lvar),staXgrd)
        if (ierr .ne. nf_noerr) then
          write(stdout,1) vinfo(1)(1:lvar), staname(1:lstr)
          goto 99                                         !--> ERROR
        endif
!
! Define float Y-position in the grid
!
        vinfo(1)='Ygrid'
        lvar=lenstr(vinfo(1))
        ierr=nf_inq_varid(ncid,vinfo(1)(1:lvar),staYgrd)
        if (ierr .ne. nf_noerr) then
          write(stdout,1) vinfo(1)(1:lvar), staname(1:lstr)
          goto 99                                         !--> ERROR
        endif
# ifdef SOLVE3D
!
!  Define depth.
!
        vinfo(1)='depth'
        lvar=lenstr(vinfo(1))
        ierr=nf_inq_varid(ncid,vinfo(1)(1:lvar),staDepth)
!
! Define  vertical positions S or depth(m)
!
#  ifndef ALL_SIGMA
        vinfo(1)='Zgrid'
        lvar=lenstr(vinfo(1))
        ierr=nf_inq_varid(ncid,vinfo(1)(1:lvar),staZgrd)
        if (ierr .ne. nf_noerr) then
          write(stdout,1) vinfo(1)(1:lvar), staname(1:lstr)
          goto 99                                         !--> ERROR
        endif
#  endif
# endif /* SOLVE3D */
       endif ! grid parameters
!
! --- Define veriables for station data ---
!
!
!  Define sea level
!
        vinfo(1)='zeta'
        lvar=lenstr(vinfo(1))
        ierr=nf_inq_varid(ncid,vinfo(1)(1:lvar),staZeta)
        if (ierr .ne. nf_noerr) then
          write(stdout,1) vinfo(1)(1:lvar), staname(1:lstr)
          goto 99                                         !--> ERROR
        endif
# ifdef SOLVE3D
!
!  Define temperature
!
       if (wrtsta(indxstaTemp)) then
        vinfo(1)='temp'
        lvar=lenstr(vinfo(1))
        ierr=nf_inq_varid(ncid,vinfo(1)(1:lvar),staTemp)
        if (ierr .ne. nf_noerr) then
          write(stdout,1) vinfo(1)(1:lvar), staname(1:lstr)
          goto 99                                         !--> ERROR
        endif
       endif
!
!  Define salinity
!
#  ifdef SALINITY
       if (wrtsta(indxstaSalt)) then
        vinfo(1)='salt'
        lvar=lenstr(vinfo(1))
        ierr=nf_inq_varid(ncid,vinfo(1)(1:lvar),staSal)
        if (ierr .ne. nf_noerr) then
          write(stdout,1) vinfo(1)(1:lvar), staname(1:lstr)
          goto 99                                         !--> ERROR
        endif
       endif
#  endif  /* SALINITY */
#  ifdef MUSTANG
       do itrc=1,NT-2
        vinfo(1)=vname(1,indxT+2+itrc-1)
        lvar=lenstr(vinfo(1))
        ierr=nf_inq_varid(ncid,vinfo(1)(1:lvar),staMUS(itrc))
        if (ierr .ne. nf_noerr) then
          write(stdout,1) vinfo(1)(1:lvar), staname(1:lstr)
          goto 99                                         !--> ERROR
        endif
       enddo
#  endif
!
!  Define density anomaly.
!
       if (wrtsta(indxstaRho)) then
        vinfo(1)='rho'
        lvar=lenstr(vinfo(1))
        ierr=nf_inq_varid(ncid,vinfo(1)(1:lvar),staDen)
        if (ierr .ne. nf_noerr) then
          write(stdout,1) vinfo(1)(1:lvar), staname(1:lstr)
          goto 99                                         !--> ERROR
        endif
       endif
!
!  Velocity components
!
       if (wrtsta(indxstaVel)) then
        vinfo(1)='v'
        lvar=lenstr(vinfo(1))
        ierr=nf_inq_varid(ncid,vinfo(1)(1:lvar),staV)
        if (ierr .ne. nf_noerr) then
          write(stdout,1) vinfo(1)(1:lvar), staname(1:lstr)
          goto 99                                         !--> ERROR
        endif
!
        vinfo(1)='u'
        lvar=lenstr(vinfo(1))
        ierr=nf_inq_varid(ncid,vinfo(1)(1:lvar),staU)
        if (ierr .ne. nf_noerr) then
          write(stdout,1) vinfo(1)(1:lvar), staname(1:lstr)
          goto 99                                         !--> ERROR
        endif
       endif
# else
!
!  Barotropic velocity components
!
       if (wrtsta(indxstaVel)) then
        vinfo(1)='v'
        lvar=lenstr(vinfo(1))
        ierr=nf_inq_varid(ncid,vinfo(1)(1:lvar),staV)
        if (ierr .ne. nf_noerr) then
          write(stdout,1) vinfo(1)(1:lvar), staname(1:lstr)
          goto 99                                         !--> ERROR
        endif
!
        vinfo(1)='u'
        lvar=lenstr(vinfo(1))
        ierr=nf_inq_varid(ncid,vinfo(1)(1:lvar),staU)
        if (ierr .ne. nf_noerr) then
          write(stdout,1) vinfo(1)(1:lvar), staname(1:lstr)
          goto 99                                         !--> ERROR
        endif
       endif
# endif  /* SOLVE3D */


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

# if defined MPI & !defined PARALLEL_FILES
      else

      if (nstas.gt.0) then
        ierr=nf_open (staname(1:lstr), nf_write, ncid)
        if (ierr .ne. nf_noerr) then
          MPI_master_only write(stdout,'(/1x,4A,2x,A,I4/)')
     &                'DEF_STA ERROR: ',
     &                'Cannot open file ''', staname(1:lstr), '''.'
     &                 MYID
       endif
          goto 99                                         !--> ERROR
        endif
# endif
      endif  ! create or open

      ierr=nf_set_fill (ncid, nf_nofill, lvar)
      if (ierr .ne. nf_noerr) then
        write(*,'(6x,2A,i4,1x,A,i4)') 'DEF_STA ERROR: Cannot ',
     &    'switch to ''nf_nofill'' more; netCDF error code =', ierr
      endif

   1  format(/1x,'DEF_STA ERROR: Cannot find variable ''',
     &                   A, ''' in netCDF file ''', A, '''.'/)

  11  format(/' DEF_STA - unable to create station file: ',a)
  20  format(/' DEF_STA - error while writing variable: ',a,
     &        /,15x,'into station file: ',a)

 99   return
      end
#else
      subroutine def_sta_empty
      return
      end
#endif /* STATIONS */
