def_dai.F

Go to the documentation of this file.

#include "cppdefs.h"
      MODULE def_dai_mod
#if defined FOUR_DVAR || defined ENKF_RESTART
!
!git $Id$
!================================================== Hernan G. Arango ===
!  Copyright (c) 2002-2025 The ROMS Group                              !
!    Licensed under a MIT/X style license                              !
!    See License_ROMS.md                                               !
!=======================================================================
!                                                                      !
!  This module creates Data Assimilation initial conditions (4D-Var    !
!  analysis) or Ensemble Kalman Filter (EnKF) restart file using       !
!  either the standard NetCDF library or the Parallel-IO (PIO)         !
!  library. It defines its dimensions, attributes, and variables.      !
!                                                                      !
!=======================================================================
!
      USE mod_param
      USE mod_parallel
# ifdef BIOLOGY
      USE mod_biology
# endif
# ifdef FOUR_DVAR
      USE mod_fourdvar
# endif
      USE mod_iounits
      USE mod_ncparam
      USE mod_scalars
# ifdef SEDIMENT
      USE mod_sediment
# endif
!
      USE def_dim_mod,  ONLY : def_dim
      USE def_info_mod, ONLY : def_info
      USE def_var_mod,  ONLY : def_var
      USE strings_mod,  ONLY : founderror
      USE wrt_info_mod, ONLY : wrt_info
!
      implicit none
!
      PUBLIC  :: def_dai
      PRIVATE :: def_dai_nf90
# if defined PIO_LIB && defined DISTRIBUTE
      PRIVATE :: def_dai_pio
# endif
!
      CONTAINS
!
!***********************************************************************
      SUBROUTINE def_dai (ng)
!***********************************************************************
!
!  Imported variable declarations.
!
      integer, intent(in) :: ng
!
!  Local variable declarations.
!
      character (len=*), parameter :: myfile =                          &
     &  __FILE__
!
!-----------------------------------------------------------------------
!  Create a new history file according to IO type.
!-----------------------------------------------------------------------
!
      SELECT CASE (dai(ng)%IOtype)
        CASE (io_nf90)
          CALL def_dai_nf90 (ng)
 
# if defined PIO_LIB && defined DISTRIBUTE
        CASE (io_pio)
          CALL def_dai_pio (ng)
# endif
        CASE DEFAULT
          IF (master) WRITE (stdout,10) dai(ng)%IOtype
          exit_flag=3
      END SELECT
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
  10  FORMAT (' DEF_DAI - Illegal output file type, io_type = ',i0,     &
     &        /,11x,'Check KeyWord ''OUT_LIB'' in ''roms.in''.')
!
      RETURN
      END SUBROUTINE def_dai
!
!***********************************************************************
      SUBROUTINE def_dai_nf90 (ng)
!***********************************************************************
!
      USE mod_netcdf
!
!  Imported variable declarations.
!
      integer, intent(in) :: ng
!
!  Local variable declarations.
!
      logical :: got_var(nv)
!
      integer, parameter :: natt = 25
 
      integer :: i, j, itrc, nvd3, nvd4
      integer :: recdim, status, varid
 
      integer :: dimids(ndimid)
      integer :: t2dgrd(3), u2dgrd(3), v2dgrd(3)
# ifdef SOLVE3D
      integer :: t3dgrd(4), u3dgrd(4), v3dgrd(4), w3dgrd(4)
# endif
!
      real(r8) :: aval(6)
!
      character (len=256)    :: ncname
      character (len=MaxLen) :: vinfo(natt)
 
      character (len=*), parameter :: myfile =                          &
     &  __FILE__//", def_dai_nf90"
!
      sourcefile=myfile
!
!-----------------------------------------------------------------------
!  Set and report file name.
!-----------------------------------------------------------------------
!
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      ncname=dai(ng)%name
!
      IF (master) THEN
        IF (ldefdai(ng)) THEN
          WRITE (stdout,10) ng, trim(ncname)
        ELSE
          WRITE (stdout,20) ng, trim(ncname)
        END IF
      END IF
!
!=======================================================================
!  Create a new Data Assimilation initial/restart NetCDF file.
!=======================================================================
!
      define : IF (ldefdai(ng)) THEN
        CALL netcdf_create (ng, inlm, trim(ncname), dai(ng)%ncid)
        IF (founderror(exit_flag, noerror, __line__, myfile)) THEN
          IF (master) WRITE (stdout,30) trim(ncname)
          RETURN
        END IF
!
!-----------------------------------------------------------------------
!  Define file dimensions.
!-----------------------------------------------------------------------
!
        dimids=0
!
        status=def_dim(ng, inlm, dai(ng)%ncid, ncname, 'xi_rho',        &
     &                 iobounds(ng)%xi_rho, dimids( 1))
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
        status=def_dim(ng, inlm, dai(ng)%ncid, ncname, 'xi_u',          &
     &                 iobounds(ng)%xi_u, dimids( 2))
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
        status=def_dim(ng, inlm, dai(ng)%ncid, ncname, 'xi_v',          &
     &                 iobounds(ng)%xi_v, dimids( 3))
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
        status=def_dim(ng, inlm, dai(ng)%ncid, ncname, 'xi_psi',        &
     &                 iobounds(ng)%xi_psi, dimids( 4))
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
        status=def_dim(ng, inlm, dai(ng)%ncid, ncname, 'eta_rho',       &
     &                 iobounds(ng)%eta_rho, dimids( 5))
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
        status=def_dim(ng, inlm, dai(ng)%ncid, ncname, 'eta_u',         &
     &                 iobounds(ng)%eta_u, dimids( 6))
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
        status=def_dim(ng, inlm, dai(ng)%ncid, ncname, 'eta_v',         &
     &                 iobounds(ng)%eta_v, dimids( 7))
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
        status=def_dim(ng, inlm, dai(ng)%ncid, ncname, 'eta_psi',       &
     &                 iobounds(ng)%eta_psi, dimids( 8))
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
# ifdef SOLVE3D
        status=def_dim(ng, inlm, dai(ng)%ncid, ncname, 'N',             &
     &                 n(ng), dimids( 9))
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
        status=def_dim(ng, inlm, dai(ng)%ncid, ncname, 's_rho',         &
     &                 n(ng), dimids( 9))
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
        status=def_dim(ng, inlm, dai(ng)%ncid, ncname, 's_w',           &
     &                 n(ng)+1, dimids(10))
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
        status=def_dim(ng, inlm, dai(ng)%ncid, ncname, 'tracer',        &
     &                 nt(ng), dimids(11))
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
#  ifdef SEDIMENT
        status=def_dim(ng, inlm, dai(ng)%ncid, ncname, 'NST',           &
     &                 nst, dimids(32))
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
        status=def_dim(ng, inlm, dai(ng)%ncid, ncname, 'Nbed',          &
     &                 nbed, dimids(16))
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
#  endif
#  ifdef ECOSIM
        status=def_dim(ng, inlm, dai(ng)%ncid, ncname, 'Nbands',        &
     &                 nbands, dimids(33))
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
        status=def_dim(ng, inlm, dai(ng)%ncid, ncname, 'Nphy',          &
     &                 nphy, dimids(25))
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
        status=def_dim(ng, inlm, dai(ng)%ncid, ncname, 'Nbac',          &
     &                 nbac, dimids(26))
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
        status=def_dim(ng, inlm, dai(ng)%ncid, ncname, 'Ndom',          &
     &                 ndom, dimids(27))
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
        status=def_dim(ng, inlm, dai(ng)%ncid, ncname, 'Nfec',          &
     &                 nfec, dimids(28))
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
#  endif
# endif
 
        status=def_dim(ng, inlm, dai(ng)%ncid, ncname, 'boundary',      &
     &                 4, dimids(14))
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
#ifdef FOUR_DVAR
        status=def_dim(ng, inlm, dai(ng)%ncid, ncname, 'Nstate',        &
     &                 nstatevar(ng), dimids(29))
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
#endif
 
        status=def_dim(ng, inlm, dai(ng)%ncid, ncname,                  &
     &                 trim(adjustl(vname(5,idtime))),                  &
     &                 nf90_unlimited, dimids(12))
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
        recdim=dimids(12)
!
!  Set number of dimensions for output variables.
!
        nvd3=3
        nvd4=4
!
!  Define dimension vectors for staggered tracer type variables.
!
        t2dgrd(1)=dimids( 1)
        t2dgrd(2)=dimids( 5)
        t2dgrd(3)=dimids(12)
# ifdef SOLVE3D
        t3dgrd(1)=dimids( 1)
        t3dgrd(2)=dimids( 5)
        t3dgrd(3)=dimids( 9)
        t3dgrd(4)=dimids(12)
# endif
!
!  Define dimension vectors for staggered u-momentum type variables.
!
        u2dgrd(1)=dimids( 2)
        u2dgrd(2)=dimids( 6)
        u2dgrd(3)=dimids(12)
# ifdef SOLVE3D
        u3dgrd(1)=dimids( 2)
        u3dgrd(2)=dimids( 6)
        u3dgrd(3)=dimids( 9)
        u3dgrd(4)=dimids(12)
# endif
!
!  Define dimension vectors for staggered v-momentum type variables.
!
        v2dgrd(1)=dimids( 3)
        v2dgrd(2)=dimids( 7)
        v2dgrd(3)=dimids(12)
# ifdef SOLVE3D
        v3dgrd(1)=dimids( 3)
        v3dgrd(2)=dimids( 7)
        v3dgrd(3)=dimids( 9)
        v3dgrd(4)=dimids(12)
# endif
# ifdef SOLVE3D
!
!  Define dimension vector for staggered w-momentum type variables.
!
        w3dgrd(1)=dimids( 1)
        w3dgrd(2)=dimids( 5)
        w3dgrd(3)=dimids(10)
        w3dgrd(4)=dimids(12)
# endif
!
!  Initialize unlimited time record dimension.
!
        dai(ng)%Rindex=0
!
!  Initialize local information variable arrays.
!
        DO i=1,natt
          DO j=1,len(vinfo(1))
            vinfo(i)(j:j)=' '
          END DO
        END DO
        DO i=1,6
          aval(i)=0.0_r8
        END DO
!
!-----------------------------------------------------------------------
!  Define time-recordless information variables.
!-----------------------------------------------------------------------
!
        CALL def_info (ng, inlm, dai(ng)%ncid, ncname, dimids)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
# ifdef SOLVE3D
!
!  Define time independent (unperturb) depths of RHO-points.
!
        vinfo( 1)=vname(1,idpthr)
        WRITE (vinfo( 2),40) trim(vname(2,idpthr))
        vinfo( 3)=vname(3,idpthr)
        vinfo(11)='downwards'
        vinfo(22)='coordinates'
        aval(5)=real(iinfo(1,idpthr,ng),r8)
        status=def_var(ng, inlm, dai(ng)%ncid, dai(ng)%Vid(idpthr),     &
     &                 nf_frst, nvd3, t3dgrd, aval, vinfo, ncname,      &
     &                 setfillval = .false.)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Define time independent (unperturb) depths of U-points.
!
        vinfo( 1)=vname(1,idpthu)
        WRITE (vinfo( 2),40) trim(vname(2,idpthu))
        vinfo( 3)=vname(3,idpthu)
        vinfo(11)='downwards'
        vinfo(22)='coordinates'
        aval(5)=real(iinfo(1,idpthu,ng),r8)
        status=def_var(ng, inlm, dai(ng)%ncid, dai(ng)%Vid(idpthu),     &
     &                 nf_frst, nvd3, u3dgrd, aval, vinfo, ncname,      &
     &                 setfillval = .false.)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Define time independent (unperturb) depths of V-points.
!
        vinfo( 1)=vname(1,idpthv)
        WRITE (vinfo( 2),40) trim(vname(2,idpthv))
        vinfo( 3)=vname(3,idpthv)
        vinfo(11)='downwards'
        vinfo(22)='coordinates'
        aval(5)=real(iinfo(1,idpthv,ng),r8)
        status=def_var(ng, inlm, dai(ng)%ncid, dai(ng)%Vid(idpthv),     &
     &                 nf_frst, nvd3, v3dgrd, aval, vinfo, ncname,      &
     &                 setfillval = .false.)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Define time independent (unperturb) depths of W-points.
!
        vinfo( 1)=vname(1,idpthw)
        WRITE (vinfo( 2),40) trim(vname(2,idpthw))
        vinfo( 3)=vname(3,idpthw)
        vinfo(11)='downwards'
        vinfo(22)='coordinates'
        aval(5)=real(iinfo(1,idpthw,ng),r8)
        status=def_var(ng, inlm, dai(ng)%ncid, dai(ng)%Vid(idpthw),     &
     &                 nf_frst, nvd3, w3dgrd, aval, vinfo, ncname,      &
     &                 setfillval = .false.)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
# endif
!
!-----------------------------------------------------------------------
!  Define time-varying variables.
!-----------------------------------------------------------------------
!
!  Define model time.
!
        vinfo( 1)=vname(1,idtime)
        vinfo( 2)=vname(2,idtime)
        WRITE (vinfo( 3),'(a,a)') 'seconds since ', trim(rclock%string)
        vinfo( 4)=trim(rclock%calendar)
        vinfo(14)=vname(4,idtime)
        vinfo(21)=vname(6,idtime)
        status=def_var(ng, inlm, dai(ng)%ncid, dai(ng)%Vid(idtime),     &
     &                 nf_tout, 1, (/recdim/), aval, vinfo, ncname,     &
     &                 setparaccess = .true.)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Define free-surface.
!
        vinfo( 1)=vname(1,idfsur)
        vinfo( 2)=vname(2,idfsur)
        vinfo( 3)=vname(3,idfsur)
        vinfo(14)=vname(4,idfsur)
        vinfo(16)=vname(1,idtime)
        vinfo(21)=vname(6,idfsur)
        vinfo(22)='coordinates'
        aval(5)=real(iinfo(1,idfsur,ng),r8)
        status=def_var(ng, inlm, dai(ng)%ncid, dai(ng)%Vid(idfsur),     &
# ifdef WET_DRY
     &                 nf_frst, nvd3, t2dgrd, aval, vinfo, ncname,      &
     &                 setfillval = .false.)
# else
     &                 nf_frst, nvd3, t2dgrd, aval, vinfo, ncname)
 
# endif
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Define 2D momentum in the XI-direction.
!
        vinfo( 1)=vname(1,idubar)
        vinfo( 2)=vname(2,idubar)
        vinfo( 3)=vname(3,idubar)
        vinfo(14)=vname(4,idubar)
        vinfo(16)=vname(1,idtime)
        vinfo(21)=vname(6,idubar)
        vinfo(22)='coordinates'
        aval(5)=real(iinfo(1,idubar,ng),r8)
        status=def_var(ng, inlm, dai(ng)%ncid, dai(ng)%Vid(idubar),     &
     &                 nf_frst, nvd3, u2dgrd, aval, vinfo, ncname)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Define 2D momentum in the ETA-direction.
!
        vinfo( 1)=vname(1,idvbar)
        vinfo( 2)=vname(2,idvbar)
        vinfo( 3)=vname(3,idvbar)
        vinfo(14)=vname(4,idvbar)
        vinfo(16)=vname(1,idtime)
        vinfo(21)=vname(6,idvbar)
        vinfo(22)='coordinates'
        aval(5)=real(iinfo(1,idvbar,ng),r8)
        status=def_var(ng, inlm, dai(ng)%ncid, dai(ng)%Vid(idvbar),     &
     &                 nf_frst, nvd3, v2dgrd, aval, vinfo, ncname)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
# ifdef SOLVE3D
!
!  Define 3D momentum component in the XI-direction.
!
        vinfo( 1)=vname(1,iduvel)
        vinfo( 2)=vname(2,iduvel)
        vinfo( 3)=vname(3,iduvel)
        vinfo(14)=vname(4,iduvel)
        vinfo(16)=vname(1,idtime)
        vinfo(21)=vname(6,iduvel)
        vinfo(22)='coordinates'
        aval(5)=real(iinfo(1,iduvel,ng),r8)
        status=def_var(ng, inlm, dai(ng)%ncid, dai(ng)%Vid(iduvel),     &
     &                 nf_frst, nvd4, u3dgrd, aval, vinfo, ncname)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Define 3D momentum component in the ETA-direction.
!
        vinfo( 1)=vname(1,idvvel)
        vinfo( 2)=vname(2,idvvel)
        vinfo( 3)=vname(3,idvvel)
        vinfo(14)=vname(4,idvvel)
        vinfo(16)=vname(1,idtime)
        vinfo(21)=vname(6,idvvel)
        vinfo(22)='coordinates'
        aval(5)=real(iinfo(1,idvvel,ng),r8)
        status=def_var(ng, inlm, dai(ng)%ncid, dai(ng)%Vid(idvvel),     &
     &                 nf_frst, nvd4, v3dgrd, aval, vinfo, ncname)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Define tracer type variables.
!
        DO itrc=1,nt(ng)
          vinfo( 1)=vname(1,idtvar(itrc))
          vinfo( 2)=vname(2,idtvar(itrc))
          vinfo( 3)=vname(3,idtvar(itrc))
          vinfo(14)=vname(4,idtvar(itrc))
          vinfo(16)=vname(1,idtime)
#  ifdef SEDIMENT
          DO i=1,nst
            IF (itrc.eq.idsed(i)) THEN
              WRITE (vinfo(19),50) 1000.0_r8*sd50(i,ng)
            END IF
          END DO
#  endif
          vinfo(21)=vname(6,idtvar(itrc))
          vinfo(22)='coordinates'
          aval(5)=real(r3dvar,r8)
          status=def_var(ng, inlm, dai(ng)%ncid, dai(ng)%Tid(itrc),     &
     &                   nf_frst, nvd4, t3dgrd, aval, vinfo, ncname)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END DO
!
!  Define vertical viscosity coefficient.
!
        vinfo( 1)=vname(1,idvvis)
        vinfo( 2)=vname(2,idvvis)
        vinfo( 3)=vname(3,idvvis)
        vinfo(14)=vname(4,idvvis)
        vinfo(16)=vname(1,idtime)
        vinfo(21)=vname(6,idvvis)
        vinfo(22)='coordinates'
        aval(5)=real(iinfo(1,idvvis,ng),r8)
        status=def_var(ng, inlm, dai(ng)%ncid, dai(ng)%Vid(idvvis),     &
     &                   nf_frst, nvd4, w3dgrd, aval, vinfo, ncname,    &
     &                   setfillval = .false.)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Define vertical diffusion coefficient for potential temperature.
!
        vinfo( 1)=vname(1,idtdif)
        vinfo( 2)=vname(2,idtdif)
        vinfo( 3)=vname(3,idtdif)
        vinfo(14)=vname(4,idtdif)
        vinfo(16)=vname(1,idtime)
        vinfo(21)=vname(6,idtdif)
        vinfo(22)='coordinates'
        aval(5)=real(iinfo(1,idtdif,ng),r8)
        status=def_var(ng, inlm, dai(ng)%ncid, dai(ng)%Vid(idtdif),     &
     &                 nf_frst, nvd4, w3dgrd, aval, vinfo, ncname,      &
     &                 setfillval = .false.)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
#  ifdef SALINITY
!
!  Define vertical diffusion coefficient for salinity.
!
        vinfo( 1)=vname(1,idsdif)
        vinfo( 2)=vname(2,idsdif)
        vinfo( 3)=vname(3,idsdif)
        vinfo(14)=vname(4,idsdif)
        vinfo(16)=vname(1,idtime)
        vinfo(21)=vname(6,idsdif)
        vinfo(22)='coordinates'
        aval(5)=real(iinfo(1,idsdif,ng),r8)
        status=def_var(ng, inlm, dai(ng)%ncid, dai(ng)%Vid(idsdif),     &
     &                 nf_frst, nvd4, w3dgrd, aval, vinfo, ncname,      &
     &                 setfillval = .false.)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
#  endif
# endif
!
!-----------------------------------------------------------------------
!  Leave definition mode.
!-----------------------------------------------------------------------
!
        CALL netcdf_enddef (ng, inlm, ncname, dai(ng)%ncid)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!-----------------------------------------------------------------------
!  Write out time-recordless, information variables.
!-----------------------------------------------------------------------
!
        CALL wrt_info (ng, inlm, dai(ng)%ncid, ncname)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
      END IF define
!
!=======================================================================
!  Open an existing Data Assimilation initial/restart file, check its
!  contents, and prepare for appending data.
!=======================================================================
!
      query : IF (.not.ldefdai(ng)) THEN
        ncname=dai(ng)%name
!
!  Open restart file for read/write.
!
        CALL netcdf_open (ng, inlm, ncname, 1, dai(ng)%ncid)
        IF (founderror(exit_flag, noerror, __line__, myfile)) THEN
          WRITE (stdout,60) trim(ncname)
          RETURN
        END IF
!
!  Inquire about the dimensions and check for consistency.
!
        CALL netcdf_check_dim (ng, inlm, ncname,                        &
     &                         ncid = dai(ng)%ncid)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Inquire about the variables.
!
        CALL netcdf_inq_var (ng, inlm, ncname,                          &
     &                       ncid = dai(ng)%ncid)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Initialize logical switches.
!
        DO i=1,nv
          got_var(i)=.false.
        END DO
!
!  Scan variable list from input NetCDF and activate switches for
!  Data Assimilation initial/restart variables. Get variable IDs.
!
        DO i=1,n_var
          IF (trim(var_name(i)).eq.trim(vname(1,idtime))) THEN
            got_var(idtime)=.true.
            dai(ng)%Vid(idtime)=var_id(i)
          ELSE IF (trim(var_name(i)).eq.trim(vname(1,idfsur))) THEN
            got_var(idfsur)=.true.
            dai(ng)%Vid(idfsur)=var_id(i)
          ELSE IF (trim(var_name(i)).eq.trim(vname(1,idubar))) THEN
            got_var(idubar)=.true.
            dai(ng)%Vid(idubar)=var_id(i)
          ELSE IF (trim(var_name(i)).eq.trim(vname(1,idvbar))) THEN
            got_var(idvbar)=.true.
            dai(ng)%Vid(idvbar)=var_id(i)
# ifdef SOLVE3D
          ELSE IF (trim(var_name(i)).eq.trim(vname(1,iduvel))) THEN
            got_var(iduvel)=.true.
            dai(ng)%Vid(iduvel)=var_id(i)
          ELSE IF (trim(var_name(i)).eq.trim(vname(1,idvvel))) THEN
            got_var(idvvel)=.true.
            dai(ng)%Vid(idvvel)=var_id(i)
          ELSE IF (trim(var_name(i)).eq.trim(vname(1,idvvis))) THEN
            got_var(idvvis)=.true.
            dai(ng)%Vid(idvvis)=var_id(i)
          ELSE IF (trim(var_name(i)).eq.trim(vname(1,idtdif))) THEN
            got_var(idtdif)=.true.
            dai(ng)%Vid(idtdif)=var_id(i)
          ELSE IF (trim(var_name(i)).eq.trim(vname(1,idsdif))) THEN
            got_var(idsdif)=.true.
            dai(ng)%Vid(idsdif)=var_id(i)
# endif
          END IF
# ifdef SOLVE3D
          DO itrc=1,nt(ng)
            IF (trim(var_name(i)).eq.trim(vname(1,idtvar(itrc)))) THEN
              got_var(idtvar(itrc))=.true.
              dai(ng)%Tid(itrc)=var_id(i)
            END IF
          END DO
# endif
        END DO
!
!  Check if initialization variables are available in input NetCDF
!  file.
!
        IF (.not.got_var(idtime)) THEN
          IF (master) WRITE (stdout,70) trim(vname(1,idtime)),          &
     &                                  trim(ncname)
          exit_flag=3
          RETURN
        END IF
        IF (.not.got_var(idfsur)) THEN
          IF (master) WRITE (stdout,70) trim(vname(1,idfsur)),          &
     &                                  trim(ncname)
          exit_flag=3
          RETURN
        END IF
        IF (.not.got_var(idubar)) THEN
          IF (master) WRITE (stdout,70) trim(vname(1,idubar)),          &
     &                                  trim(ncname)
          exit_flag=3
          RETURN
        END IF
        IF (.not.got_var(idvbar)) THEN
          IF (master) WRITE (stdout,70) trim(vname(1,idvbar)),          &
     &                                  trim(ncname)
          exit_flag=3
          RETURN
        END IF
# ifdef SOLVE3D
        IF (.not.got_var(iduvel)) THEN
          IF (master) WRITE (stdout,70) trim(vname(1,iduvel)),          &
     &                                  trim(ncname)
          exit_flag=3
          RETURN
        END IF
        IF (.not.got_var(idvvel)) THEN
          IF (master) WRITE (stdout,70) trim(vname(1,idvvel)),          &
     &                                  trim(ncname)
          exit_flag=3
          RETURN
        END IF
        DO itrc=1,nt(ng)
          IF (.not.got_var(idtvar(itrc))) THEN
            IF (master) WRITE (stdout,70) trim(vname(1,idtvar(itrc))),  &
     &                                    trim(ncname)
            exit_flag=3
            RETURN
          END IF
        END DO
        IF (.not.got_var(idvvis)) THEN
          IF (master) WRITE (stdout,70) trim(vname(1,idvvis)),          &
     &                                  trim(ncname)
          exit_flag=3
          RETURN
        END IF
        IF (.not.got_var(idtdif)) THEN
          IF (master) WRITE (stdout,70) trim(vname(1,idtdif)),          &
     &                                  trim(ncname)
          exit_flag=3
          RETURN
        END IF
#  ifdef SALINITY
        IF (.not.got_var(idsdif)) THEN
          IF (master) WRITE (stdout,70) trim(vname(1,idsdif)),          &
     &                                  trim(ncname)
          exit_flag=3
          RETURN
        END IF
#  endif
# endif
!
!  Set unlimited time record dimension to current value.
!
        dai(ng)%Rindex=rec_size
      END IF query
!
  10  FORMAT (/,2x,'DEF_DAI_NF90     - creating DA INI/RST file,',t56,  &
     &        'Grid ',i2.2,': ',a)
  20  FORMAT (/,2x,'DEF_DAI_NF90     - inquiring DA INI/RST file,',t56, &
     &        'Grid ',i2.2,': ',a)
  30  FORMAT (/,' DEF_DAI_NF90 - unable to create DA initial/restart',  &
     &        ' NetCDF file: ',a)
  40  FORMAT ('time independent',1x,a)
  50  FORMAT (1pe11.4,1x,'millimeter')
  60  FORMAT (/,' DEF_DAI_NF90 - unable to open DA initial/restart'     &
     &        ' NetCDF file: ',a)
  70  FORMAT (/,' DEF_DAI_NF90 - unable to find variable: ',a,2x,       &
     &        ' in DA initial/restart NetCDF file: ',a)
!
      RETURN
      END SUBROUTINE def_dai_nf90
 
# if defined PIO_LIB && defined DISTRIBUTE
!
!***********************************************************************
      SUBROUTINE def_dai_pio (ng)
!***********************************************************************
!
      USE mod_pio_netcdf
!
!  Imported variable declarations.
!
      integer, intent(in) :: ng
!
!  Local variable declarations.
!
      logical :: got_var(nv)
!
      integer, parameter :: natt = 25
 
      integer :: i, j, itrc, nvd3, nvd4
      integer :: recdim, status
 
      integer :: dimids(ndimid)
      integer :: t2dgrd(3), u2dgrd(3), v2dgrd(3)
#  ifdef SOLVE3D
      integer :: t3dgrd(4), u3dgrd(4), v3dgrd(4), w3dgrd(4)
#  endif
!
      real(r8) :: aval(6)
!
      character (len=256)    :: ncname
      character (len=MaxLen) :: vinfo(natt)
 
      character (len=*), parameter :: myfile =                          &
     &  __FILE__//", def_dai_pio"
!
      sourcefile=myfile
!
!-----------------------------------------------------------------------
!  Set and report file name.
!-----------------------------------------------------------------------
!
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      ncname=dai(ng)%name
!
      IF (master) THEN
        IF (ldefdai(ng)) THEN
          WRITE (stdout,10) ng, trim(ncname)
        ELSE
          WRITE (stdout,20) ng, trim(ncname)
        END IF
      END IF
!
!=======================================================================
!  Create a new Data Assimilation initial/restart NetCDF file.
!=======================================================================
!
      define : IF (ldefdai(ng)) THEN
        CALL pio_netcdf_create (ng, inlm, trim(ncname), dai(ng)%pioFile)
        IF (founderror(exit_flag, noerror, __line__, myfile)) THEN
          IF (master) WRITE (stdout,30) trim(ncname)
          RETURN
        END IF
!
!-----------------------------------------------------------------------
!  Define file dimensions.
!-----------------------------------------------------------------------
!
        dimids=0
!
        status=def_dim(ng, inlm, dai(ng)%pioFile, ncname, 'xi_rho',     &
     &                 iobounds(ng)%xi_rho, dimids( 1))
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
        status=def_dim(ng, inlm, dai(ng)%pioFile, ncname, 'xi_u',       &
     &                 iobounds(ng)%xi_u, dimids( 2))
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
        status=def_dim(ng, inlm, dai(ng)%pioFile, ncname, 'xi_v',       &
     &                 iobounds(ng)%xi_v, dimids( 3))
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
        status=def_dim(ng, inlm, dai(ng)%pioFile, ncname, 'xi_psi',     &
     &                 iobounds(ng)%xi_psi, dimids( 4))
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
        status=def_dim(ng, inlm, dai(ng)%pioFile, ncname, 'eta_rho',    &
     &                 iobounds(ng)%eta_rho, dimids( 5))
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
        status=def_dim(ng, inlm, dai(ng)%pioFile, ncname, 'eta_u',      &
     &                 iobounds(ng)%eta_u, dimids( 6))
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
        status=def_dim(ng, inlm, dai(ng)%pioFile, ncname, 'eta_v',      &
     &                 iobounds(ng)%eta_v, dimids( 7))
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
        status=def_dim(ng, inlm, dai(ng)%pioFile, ncname, 'eta_psi',    &
     &                 iobounds(ng)%eta_psi, dimids( 8))
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
#  ifdef SOLVE3D
        status=def_dim(ng, inlm, dai(ng)%pioFile, ncname, 'N',          &
     &                 n(ng), dimids( 9))
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
        status=def_dim(ng, inlm, dai(ng)%pioFile, ncname, 's_rho',      &
     &                 n(ng), dimids( 9))
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
        status=def_dim(ng, inlm, dai(ng)%pioFile, ncname, 's_w',        &
     &                 n(ng)+1, dimids(10))
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
        status=def_dim(ng, inlm, dai(ng)%pioFile, ncname, 'tracer',     &
     &                 nt(ng), dimids(11))
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
#   ifdef SEDIMENT
        status=def_dim(ng, inlm, dai(ng)%pioFile, ncname, 'NST',        &
     &                 nst, dimids(32))
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
        status=def_dim(ng, inlm, dai(ng)%pioFile, ncname, 'Nbed',       &
     &                 nbed, dimids(16))
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
#   endif
#   ifdef ECOSIM
        status=def_dim(ng, inlm, dai(ng)%pioFile, ncname, 'Nbands',     &
     &                 nbands, dimids(33))
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
        status=def_dim(ng, inlm, dai(ng)%pioFile, ncname, 'Nphy',       &
     &                 nphy, dimids(25))
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
        status=def_dim(ng, inlm, dai(ng)%pioFile, ncname, 'Nbac',       &
     &                 nbac, dimids(26))
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
        status=def_dim(ng, inlm, dai(ng)%pioFile, ncname, 'Ndom',       &
     &                 ndom, dimids(27))
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
        status=def_dim(ng, inlm, dai(ng)%pioFile, ncname, 'Nfec',       &
     &                 nfec, dimids(28))
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
#   endif
#  endif
 
        status=def_dim(ng, inlm, dai(ng)%pioFile, ncname, 'boundary',   &
     &                 4, dimids(14))
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
#  ifdef FOUR_DVAR
        status=def_dim(ng, inlm, dai(ng)%pioFile, ncname, 'Nstate',     &
     &                 nstatevar(ng), dimids(29))
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
#  endif
 
        status=def_dim(ng, inlm, dai(ng)%pioFile, ncname,               &
     &                 trim(adjustl(vname(5,idtime))),                  &
     &                 pio_unlimited, dimids(12))
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
        recdim=dimids(12)
!
!  Set number of dimensions for output variables.
!
        nvd3=3
        nvd4=4
!
!  Define dimension vectors for staggered tracer type variables.
!
        t2dgrd(1)=dimids( 1)
        t2dgrd(2)=dimids( 5)
        t2dgrd(3)=dimids(12)
#  ifdef SOLVE3D
        t3dgrd(1)=dimids( 1)
        t3dgrd(2)=dimids( 5)
        t3dgrd(3)=dimids( 9)
        t3dgrd(4)=dimids(12)
#  endif
!
!  Define dimension vectors for staggered u-momentum type variables.
!
        u2dgrd(1)=dimids( 2)
        u2dgrd(2)=dimids( 6)
        u2dgrd(3)=dimids(12)
#  ifdef SOLVE3D
        u3dgrd(1)=dimids( 2)
        u3dgrd(2)=dimids( 6)
        u3dgrd(3)=dimids( 9)
        u3dgrd(4)=dimids(12)
#  endif
!
!  Define dimension vectors for staggered v-momentum type variables.
!
        v2dgrd(1)=dimids( 3)
        v2dgrd(2)=dimids( 7)
        v2dgrd(3)=dimids(12)
#  ifdef SOLVE3D
        v3dgrd(1)=dimids( 3)
        v3dgrd(2)=dimids( 7)
        v3dgrd(3)=dimids( 9)
        v3dgrd(4)=dimids(12)
#  endif
#  ifdef SOLVE3D
!
!  Define dimension vector for staggered w-momentum type variables.
!
        w3dgrd(1)=dimids( 1)
        w3dgrd(2)=dimids( 5)
        w3dgrd(3)=dimids(10)
        w3dgrd(4)=dimids(12)
#  endif
!
!  Initialize unlimited time record dimension.
!
        dai(ng)%Rindex=0
!
!  Initialize local information variable arrays.
!
        DO i=1,natt
          DO j=1,len(vinfo(1))
            vinfo(i)(j:j)=' '
          END DO
        END DO
        DO i=1,6
          aval(i)=0.0_r8
        END DO
!
!-----------------------------------------------------------------------
!  Define time-recordless information variables.
!-----------------------------------------------------------------------
!
        CALL def_info (ng, inlm, dai(ng)%pioFile, ncname, dimids)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
#  ifdef SOLVE3D
!
!  Define time independent (unperturb) depths of RHO-points.
!
        vinfo( 1)=vname(1,idpthr)
        WRITE (vinfo( 2),40) trim(vname(2,idpthr))
        vinfo( 3)=vname(3,idpthr)
        vinfo(11)='downwards'
        vinfo(22)='coordinates'
        aval(5)=real(iinfo(1,idpthr,ng),r8)
        dai(ng)%pioVar(idpthr)%dkind=pio_frst
        dai(ng)%pioVar(idpthr)%gtype=r3dvar
!
        status=def_var(ng, inlm, dai(ng)%pioFile,                       &
     &                 dai(ng)%pioVar(idpthr)%vd,                       &
     &                 pio_frst, nvd3, t3dgrd, aval, vinfo, ncname,     &
     &                 setfillval = .false.)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Define time independent (unperturb) depths of U-points.
!
        vinfo( 1)=vname(1,idpthu)
        WRITE (vinfo( 2),40) trim(vname(2,idpthu))
        vinfo( 3)=vname(3,idpthu)
        vinfo(11)='downwards'
        vinfo(22)='coordinates'
        aval(5)=real(iinfo(1,idpthu,ng),r8)
        dai(ng)%pioVar(idpthu)%dkind=pio_frst
        dai(ng)%pioVar(idpthu)%gtype=u3dvar
!
        status=def_var(ng, inlm, dai(ng)%pioFile,                       &
     &                 dai(ng)%pioVar(idpthu)%vd,                       &
     &                 pio_frst, nvd3, u3dgrd, aval, vinfo, ncname,     &
     &                 setfillval = .false.)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Define time independent (unperturb) depths of V-points.
!
        vinfo( 1)=vname(1,idpthv)
        WRITE (vinfo( 2),40) trim(vname(2,idpthv))
        vinfo( 3)=vname(3,idpthv)
        vinfo(11)='downwards'
        vinfo(22)='coordinates'
        aval(5)=real(iinfo(1,idpthv,ng),r8)
        dai(ng)%pioVar(idpthv)%dkind=pio_frst
        dai(ng)%pioVar(idpthv)%gtype=v3dvar
!
        status=def_var(ng, inlm, dai(ng)%pioFile,                       &
     &                 dai(ng)%pioVar(idpthv)%vd,                       &
     &                 pio_frst, nvd3, v3dgrd, aval, vinfo, ncname,     &
     &                 setfillval = .false.)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Define time independent (unperturb) depths of W-points.
!
        vinfo( 1)=vname(1,idpthw)
        WRITE (vinfo( 2),40) trim(vname(2,idpthw))
        vinfo( 3)=vname(3,idpthw)
        vinfo(11)='downwards'
        vinfo(22)='coordinates'
        aval(5)=real(iinfo(1,idpthw,ng),r8)
        dai(ng)%pioVar(idpthw)%dkind=pio_frst
        dai(ng)%pioVar(idpthw)%gtype=w3dvar
!
        status=def_var(ng, inlm, dai(ng)%pioFile,                       &
     &                 dai(ng)%pioVar(idpthw)%vd,                       &
     &                 pio_frst, nvd3, w3dgrd, aval, vinfo, ncname,     &
     &                 setfillval = .false.)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
#  endif
!
!-----------------------------------------------------------------------
!  Define time-varying variables.
!-----------------------------------------------------------------------
!
!  Define model time.
!
        vinfo( 1)=vname(1,idtime)
        vinfo( 2)=vname(2,idtime)
        WRITE (vinfo( 3),'(a,a)') 'seconds since ', trim(rclock%string)
        vinfo( 4)=trim(rclock%calendar)
        vinfo(14)=vname(4,idtime)
        vinfo(21)=vname(6,idtime)
        dai(ng)%pioVar(idtime)%dkind=pio_tout
        dai(ng)%pioVar(idtime)%gtype=0
!
        status=def_var(ng, inlm, dai(ng)%pioFile,                       &
     &                 dai(ng)%pioVar(idtime)%vd,                       &
     &                 pio_tout, 1, (/recdim/), aval, vinfo, ncname,    &
     &                 setparaccess = .true.)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Define free-surface.
!
        vinfo( 1)=vname(1,idfsur)
        vinfo( 2)=vname(2,idfsur)
        vinfo( 3)=vname(3,idfsur)
        vinfo(14)=vname(4,idfsur)
        vinfo(16)=vname(1,idtime)
        vinfo(21)=vname(6,idfsur)
        vinfo(22)='coordinates'
        aval(5)=real(iinfo(1,idfsur,ng),r8)
        dai(ng)%pioVar(idfsur)%dkind=pio_frst
        dai(ng)%pioVar(idfsur)%gtype=r2dvar
!
        status=def_var(ng, inlm, dai(ng)%pioFile,                       &
     &                 dai(ng)%pioVar(idfsur)%vd,                       &
#  ifdef WET_DRY
     &                 pio_frst, nvd3, t2dgrd, aval, vinfo, ncname,     &
     &                 setfillval = .false.)
#  else
     &                 pio_frst, nvd3, t2dgrd, aval, vinfo, ncname)
 
#  endif
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Define 2D momentum in the XI-direction.
!
        vinfo( 1)=vname(1,idubar)
        vinfo( 2)=vname(2,idubar)
        vinfo( 3)=vname(3,idubar)
        vinfo(14)=vname(4,idubar)
        vinfo(16)=vname(1,idtime)
        vinfo(21)=vname(6,idubar)
        vinfo(22)='coordinates'
        aval(5)=real(iinfo(1,idubar,ng),r8)
        dai(ng)%pioVar(idubar)%dkind=pio_frst
        dai(ng)%pioVar(idubar)%gtype=u2dvar
!
        status=def_var(ng, inlm, dai(ng)%pioFile,                       &
     &                 dai(ng)%pioVar(idubar)%vd,                       &
     &                 pio_frst, nvd3, u2dgrd, aval, vinfo, ncname)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Define 2D momentum in the ETA-direction.
!
        vinfo( 1)=vname(1,idvbar)
        vinfo( 2)=vname(2,idvbar)
        vinfo( 3)=vname(3,idvbar)
        vinfo(14)=vname(4,idvbar)
        vinfo(16)=vname(1,idtime)
        vinfo(21)=vname(6,idvbar)
        vinfo(22)='coordinates'
        aval(5)=real(iinfo(1,idvbar,ng),r8)
        dai(ng)%pioVar(idvbar)%dkind=pio_frst
        dai(ng)%pioVar(idvbar)%gtype=v2dvar
!
        status=def_var(ng, inlm, dai(ng)%pioFile,                       &
     &                 dai(ng)%pioVar(idvbar)%vd,                       &
     &                 pio_frst, nvd3, v2dgrd, aval, vinfo, ncname)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
#  ifdef SOLVE3D
!
!  Define 3D momentum component in the XI-direction.
!
        vinfo( 1)=vname(1,iduvel)
        vinfo( 2)=vname(2,iduvel)
        vinfo( 3)=vname(3,iduvel)
        vinfo(14)=vname(4,iduvel)
        vinfo(16)=vname(1,idtime)
        vinfo(21)=vname(6,iduvel)
        vinfo(22)='coordinates'
        aval(5)=real(iinfo(1,iduvel,ng),r8)
        dai(ng)%pioVar(iduvel)%dkind=pio_frst
        dai(ng)%pioVar(iduvel)%gtype=u3dvar
!
        status=def_var(ng, inlm, dai(ng)%pioFile,                       &
     &                 dai(ng)%pioVar(iduvel)%vd,                       &
     &                 pio_frst, nvd4, u3dgrd, aval, vinfo, ncname)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Define 3D momentum component in the ETA-direction.
!
        vinfo( 1)=vname(1,idvvel)
        vinfo( 2)=vname(2,idvvel)
        vinfo( 3)=vname(3,idvvel)
        vinfo(14)=vname(4,idvvel)
        vinfo(16)=vname(1,idtime)
        vinfo(21)=vname(6,idvvel)
        vinfo(22)='coordinates'
        aval(5)=real(iinfo(1,idvvel,ng),r8)
        dai(ng)%pioVar(idvvel)%dkind=pio_frst
        dai(ng)%pioVar(idvvel)%gtype=v3dvar
!
        status=def_var(ng, inlm, dai(ng)%pioFile,                       &
     &                 dai(ng)%pioVar(idvvel)%vd,                       &
     &                 pio_frst, nvd4, v3dgrd, aval, vinfo, ncname)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Define tracer type variables.
!
        DO itrc=1,nt(ng)
          vinfo( 1)=vname(1,idtvar(itrc))
          vinfo( 2)=vname(2,idtvar(itrc))
          vinfo( 3)=vname(3,idtvar(itrc))
          vinfo(14)=vname(4,idtvar(itrc))
          vinfo(16)=vname(1,idtime)
#   ifdef SEDIMENT
          DO i=1,nst
            IF (itrc.eq.idsed(i)) THEN
              WRITE (vinfo(19),50) 1000.0_r8*sd50(i,ng)
            END IF
          END DO
#   endif
          vinfo(21)=vname(6,idtvar(itrc))
          vinfo(22)='coordinates'
          aval(5)=real(r3dvar,r8)
          dai(ng)%pioTrc(itrc)%dkind=pio_frst
          dai(ng)%pioTrc(itrc)%gtype=r3dvar
!
          status=def_var(ng, inlm, dai(ng)%pioFile,                     &
     &                   dai(ng)%pioTrc(itrc)%vd,                       &
     &                   pio_frst, nvd4, t3dgrd, aval, vinfo, ncname)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END DO
!
!  Define vertical viscosity coefficient.
!
        vinfo( 1)=vname(1,idvvis)
        vinfo( 2)=vname(2,idvvis)
        vinfo( 3)=vname(3,idvvis)
        vinfo(14)=vname(4,idvvis)
        vinfo(16)=vname(1,idtime)
        vinfo(21)=vname(6,idvvis)
        vinfo(22)='coordinates'
        aval(5)=real(iinfo(1,idvvis,ng),r8)
        dai(ng)%pioVar(idvvis)%dkind=pio_frst
        dai(ng)%pioVar(idvvis)%gtype=w3dvar
!
        status=def_var(ng, inlm, dai(ng)%pioFile,                       &
     &                 dai(ng)%pioVar(idvvis)%vd,                       &
     &                   pio_frst, nvd4, w3dgrd, aval, vinfo, ncname,   &
     &                   setfillval = .false.)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Define vertical diffusion coefficient for potential temperature.
!
        vinfo( 1)=vname(1,idtdif)
        vinfo( 2)=vname(2,idtdif)
        vinfo( 3)=vname(3,idtdif)
        vinfo(14)=vname(4,idtdif)
        vinfo(16)=vname(1,idtime)
        vinfo(21)=vname(6,idtdif)
        vinfo(22)='coordinates'
        aval(5)=real(iinfo(1,idtdif,ng),r8)
        dai(ng)%pioVar(idtdif)%dkind=pio_frst
        dai(ng)%pioVar(idtdif)%gtype=w3dvar
!
        status=def_var(ng, inlm, dai(ng)%pioFile,                       &
     &                 dai(ng)%pioVar(idtdif)%vd,                       &
     &                 pio_frst, nvd4, w3dgrd, aval, vinfo, ncname,     &
     &                 setfillval = .false.)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
#   ifdef SALINITY
!
!  Define vertical diffusion coefficient for salinity.
!
        vinfo( 1)=vname(1,idsdif)
        vinfo( 2)=vname(2,idsdif)
        vinfo( 3)=vname(3,idsdif)
        vinfo(14)=vname(4,idsdif)
        vinfo(16)=vname(1,idtime)
        vinfo(21)=vname(6,idsdif)
        vinfo(22)='coordinates'
        aval(5)=real(iinfo(1,idsdif,ng),r8)
        dai(ng)%pioVar(idsdif)%dkind=pio_frst
        dai(ng)%pioVar(idsdif)%gtype=w3dvar
!
        status=def_var(ng, inlm, dai(ng)%pioFile,                       &
     &                 dai(ng)%pioVar(idsdif)%vd,                       &
     &                 pio_frst, nvd4, w3dgrd, aval, vinfo, ncname,     &
     &                 setfillval = .false.)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
#   endif
#  endif
!
!-----------------------------------------------------------------------
!  Leave definition mode.
!-----------------------------------------------------------------------
!
        CALL pio_netcdf_enddef (ng, inlm, ncname, dai(ng)%pioFile)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!-----------------------------------------------------------------------
!  Write out time-recordless, information variables.
!-----------------------------------------------------------------------
!
        CALL wrt_info (ng, inlm, dai(ng)%pioFile, ncname)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
      END IF define
!
!=======================================================================
!  Open an existing Data Assimilation initial/restart file, check its
!  contents, and prepare for appending data.
!=======================================================================
!
      query : IF (.not.ldefdai(ng)) THEN
        ncname=dai(ng)%name
!
!  Open restart file for read/write.
!
        CALL pio_netcdf_open (ng, inlm, ncname, 1, dai(ng)%pioFile)
        IF (founderror(exit_flag, noerror, __line__, myfile)) THEN
          WRITE (stdout,60) trim(ncname)
          RETURN
        END IF
!
!  Inquire about the dimensions and check for consistency.
!
        CALL pio_netcdf_check_dim (ng, inlm, ncname,                    &
     &                             piofile = dai(ng)%pioFile)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Inquire about the variables.
!
        CALL pio_netcdf_inq_var (ng, inlm, ncname,                      &
     &                           piofile= dai(ng)%pioFile)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Initialize logical switches.
!
        DO i=1,nv
          got_var(i)=.false.
        END DO
!
!  Scan variable list from input NetCDF and activate switches for
!  Data Assimilation initial/restart variables. Get variable IDs.
!
        DO i=1,n_var
          IF (trim(var_name(i)).eq.trim(vname(1,idtime))) THEN
            got_var(idtime)=.true.
            dai(ng)%pioVar(idtime)%vd=var_desc(i)
            dai(ng)%pioVar(idtime)%dkind=pio_tout
            dai(ng)%pioVar(idtime)%gtype=0
          ELSE IF (trim(var_name(i)).eq.trim(vname(1,idfsur))) THEN
            got_var(idfsur)=.true.
            dai(ng)%pioVar(idfsur)%vd=var_desc(i)
            dai(ng)%pioVar(idfsur)%dkind=pio_frst
            dai(ng)%pioVar(idfsur)%gtype=r2dvar
          ELSE IF (trim(var_name(i)).eq.trim(vname(1,idubar))) THEN
            got_var(idubar)=.true.
            dai(ng)%pioVar(idubar)%vd=var_desc(i)
            dai(ng)%pioVar(idubar)%dkind=pio_frst
            dai(ng)%pioVar(idubar)%gtype=u2dvar
          ELSE IF (trim(var_name(i)).eq.trim(vname(1,idvbar))) THEN
            got_var(idvbar)=.true.
            dai(ng)%pioVar(idvbar)%vd=var_desc(i)
            dai(ng)%pioVar(idvbar)%dkind=pio_frst
            dai(ng)%pioVar(idvbar)%gtype=v2dvar
#  ifdef SOLVE3D
          ELSE IF (trim(var_name(i)).eq.trim(vname(1,iduvel))) THEN
            got_var(iduvel)=.true.
            dai(ng)%pioVar(iduvel)%vd=var_desc(i)
            dai(ng)%pioVar(iduvel)%dkind=pio_frst
            dai(ng)%pioVar(iduvel)%gtype=u3dvar
          ELSE IF (trim(var_name(i)).eq.trim(vname(1,idvvel))) THEN
            got_var(idvvel)=.true.
            dai(ng)%pioVar(idvvel)%vd=var_desc(i)
            dai(ng)%pioVar(idvvel)%dkind=pio_frst
            dai(ng)%pioVar(idvvel)%gtype=v3dvar
          ELSE IF (trim(var_name(i)).eq.trim(vname(1,idvvis))) THEN
            got_var(idvvis)=.true.
            dai(ng)%pioVar(idvvis)%vd=var_desc(i)
            dai(ng)%pioVar(idvvis)%dkind=pio_frst
            dai(ng)%pioVar(idvvis)%gtype=w3dvar
          ELSE IF (trim(var_name(i)).eq.trim(vname(1,idtdif))) THEN
            got_var(idtdif)=.true.
            dai(ng)%pioVar(idtdif)%vd=var_desc(i)
            dai(ng)%pioVar(idtdif)%dkind=pio_frst
            dai(ng)%pioVar(idtdif)%gtype=w3dvar
          ELSE IF (trim(var_name(i)).eq.trim(vname(1,idsdif))) THEN
            got_var(idsdif)=.true.
            dai(ng)%pioVar(idsdif)%vd=var_desc(i)
            dai(ng)%pioVar(idsdif)%dkind=pio_frst
            dai(ng)%pioVar(idsdif)%gtype=w3dvar
#  endif
          END IF
#  ifdef SOLVE3D
          DO itrc=1,nt(ng)
            IF (trim(var_name(i)).eq.trim(vname(1,idtvar(itrc)))) THEN
              got_var(idtvar(itrc))=.true.
              dai(ng)%pioTrc(itrc)%vd=var_desc(i)
              dai(ng)%pioTrc(itrc)%dkind=pio_frst
              dai(ng)%pioTrc(itrc)%gtype=r3dvar
            END IF
          END DO
#  endif
        END DO
!
!  Check if initialization variables are available in input NetCDF
!  file.
!
        IF (.not.got_var(idtime)) THEN
          IF (master) WRITE (stdout,70) trim(vname(1,idtime)),          &
     &                                  trim(ncname)
          exit_flag=3
          RETURN
        END IF
        IF (.not.got_var(idfsur)) THEN
          IF (master) WRITE (stdout,70) trim(vname(1,idfsur)),          &
     &                                  trim(ncname)
          exit_flag=3
          RETURN
        END IF
        IF (.not.got_var(idubar)) THEN
          IF (master) WRITE (stdout,70) trim(vname(1,idubar)),          &
     &                                  trim(ncname)
          exit_flag=3
          RETURN
        END IF
        IF (.not.got_var(idvbar)) THEN
          IF (master) WRITE (stdout,70) trim(vname(1,idvbar)),          &
     &                                  trim(ncname)
          exit_flag=3
          RETURN
        END IF
#  ifdef SOLVE3D
        IF (.not.got_var(iduvel)) THEN
          IF (master) WRITE (stdout,70) trim(vname(1,iduvel)),          &
     &                                  trim(ncname)
          exit_flag=3
          RETURN
        END IF
        IF (.not.got_var(idvvel)) THEN
          IF (master) WRITE (stdout,70) trim(vname(1,idvvel)),          &
     &                                  trim(ncname)
          exit_flag=3
          RETURN
        END IF
        DO itrc=1,nt(ng)
          IF (.not.got_var(idtvar(itrc))) THEN
            IF (master) WRITE (stdout,70) trim(vname(1,idtvar(itrc))),  &
     &                                    trim(ncname)
            exit_flag=3
            RETURN
          END IF
        END DO
        IF (.not.got_var(idvvis)) THEN
          IF (master) WRITE (stdout,70) trim(vname(1,idvvis)),          &
     &                                  trim(ncname)
          exit_flag=3
          RETURN
        END IF
        IF (.not.got_var(idtdif)) THEN
          IF (master) WRITE (stdout,70) trim(vname(1,idtdif)),          &
     &                                  trim(ncname)
          exit_flag=3
          RETURN
        END IF
#   ifdef SALINITY
        IF (.not.got_var(idsdif)) THEN
          IF (master) WRITE (stdout,70) trim(vname(1,idsdif)),          &
     &                                  trim(ncname)
          exit_flag=3
          RETURN
        END IF
#   endif
#  endif
!
!  Set unlimited time record dimension to current value.
!
        dai(ng)%Rindex=rec_size
      END IF query
!
  10  FORMAT (/,2x,'DEF_DAI_PIO      - creating DA INI/RST file,',t56,  &
     &        'Grid ',i2.2,': ',a)
  20  FORMAT (/,2x,'DEF_DAI_PIO      - inquiring DA INI/RST file,',t56, &
     &        'Grid ',i2.2,': ',a)
  30  FORMAT (/,' DEF_DAI_PIO - unable to create DA initial/restart',   &
     &        ' NetCDF file: ',a)
  40  FORMAT ('time independent',1x,a)
  50  FORMAT (1pe11.4,1x,'millimeter')
  60  FORMAT (/,' DEF_DAI_PIO - unable to open DA initial/restart'      &
     &        ' NetCDF file: ',a)
  70  FORMAT (/,' DEF_DAI_PIO - unable to find variable: ',a,2x,        &
     &        ' in DA initial/restart NetCDF file: ',a)
!
      RETURN
      END SUBROUTINE def_dai_pio
# endif
#endif
      END MODULE def_dai_mod