rp_def_ini.F

Go to the documentation of this file.

#include "cppdefs.h"
      MODULE rp_def_ini_mod
#if defined TL_IOMS && defined FOUR_DVAR
!
!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 opens existing representer model initial conditions     !
!  file and either define new new variables or inquire its contents.   !
!  It either use the standard NetCDF library or the Parallel-IO (PIO)  !
!  library.                                                            !
!                                                                      !
!=======================================================================
!
      USE mod_param
      USE mod_parallel
      USE mod_iounits
      USE mod_ncparam
      USE mod_scalars
!
      USE def_dim_mod, ONLY : def_dim
      USE def_var_mod, ONLY : def_var
      USE strings_mod, ONLY : founderror
!
      implicit none
!
      PUBLIC  :: rp_def_ini
      PRIVATE :: rp_def_ini_nf90
# if defined PIO_LIB && defined DISTRIBUTE
      PRIVATE :: rp_def_ini_pio
# endif
!
      CONTAINS
!
!***********************************************************************
      SUBROUTINE rp_def_ini (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 (irp(ng)%IOtype)
        CASE (io_nf90)
          CALL rp_def_ini_nf90 (ng)
 
# if defined PIO_LIB && defined DISTRIBUTE
        CASE (io_pio)
          CALL rp_def_ini_pio (ng)
# endif
        CASE DEFAULT
          IF (master) WRITE (stdout,10) irp(ng)%IOtype
          exit_flag=3
      END SELECT
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
  10  FORMAT (' RP_DEF_INI - Illegal output type, io_type = ',i0,       &
     &        /,14x,'Check KeyWord ''OUT_LIB'' in ''roms.in''.')
!
      RETURN
      END SUBROUTINE rp_def_ini
!
!***********************************************************************
      SUBROUTINE rp_def_ini_nf90 (ng)
!***********************************************************************
!
      USE mod_netcdf
!
!  Imported variable declarations.
!
      integer, intent(in) :: ng
!
!  Local variable declarations.
!
      logical :: ldefine = .false.
# ifdef ADJUST_BOUNDARY
      logical :: got_iorj, got_boundary, got_obc_adjust
# endif
# if defined ADJUST_STFLUX || defined ADJUST_WSTRESS
      logical :: got_frc_adjust
# endif
      logical :: got_var(nv)
!
      integer, parameter :: natt = 25
 
      integer :: i, j, ifield, itrc, status
      integer :: fcount
# ifdef ADJUST_BOUNDARY
      integer :: iorjdim, brecdim
# endif
# if defined ADJUST_STFLUX || defined ADJUST_WSTRESS
      integer :: frecdim
# endif
      integer :: dimids(ndimid)
# ifdef ADJUST_BOUNDARY
      integer :: t2dobc(4)
#  ifdef SOLVE3D
      integer :: t3dobc(5)
#  endif
# endif
# if defined ADJUST_STFLUX || defined ADJUST_WSTRESS
      integer :: t4dfrc(4), u4dfrc(4), v4dfrc(4)
# endif
!
      real(r8) :: aval(6)
!
      character (len=256)    :: ncname
      character (len=MaxLen) :: vinfo(natt)
 
      character (len=*), parameter :: myfile =                          &
     &  __FILE__//", rp_def_ini_nf90"
!
      sourcefile=myfile
 
# if defined ADJUST_BOUNDARY || \
     defined adjust_stflux   || defined adjust_wstress
!
!=======================================================================
!  Open existing representer model initial conditions and define new
!  variables.
!=======================================================================
!
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      ncname=irp(ng)%name
!
!  Open representer model initialization file for read/write.
!
      CALL netcdf_open (ng, irpm, ncname, 1, irp(ng)%ncid)
      IF (founderror(exit_flag, noerror, __line__, myfile)) THEN
        IF (master) WRITE (stdout,10) trim(ncname)
        RETURN
      END IF
!
!  Inquire about the dimensions and check for consistency.
!
      CALL netcdf_check_dim (ng, irpm, ncname,                          &
     &                       ncid = irp(ng)%ncid)
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Inquire about the variables.
!
      CALL netcdf_inq_var (ng, irpm, ncname,                            &
     &                     ncid = irp(ng)%ncid)
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Check if surface forcing variables have been already defined.
!
      DO i=1,nv
        got_var(i)=.false.
      END DO
!
      DO i=1,n_var
#  ifdef ADJUST_BOUNDARY
        IF (trim(var_name(i)).eq.                                       &
     &      trim(vname(1,idsbry(isfsur)))) THEN
          got_var(idsbry(isfsur))=.true.
          irp(ng)%Vid(idsbry(isfsur))=var_id(i)
        ELSE IF (trim(var_name(i)).eq.                                  &
     &           trim(vname(1,idsbry(isubar)))) THEN
          got_var(idsbry(isubar))=.true.
          irp(ng)%Vid(idsbry(isubar))=var_id(i)
        ELSE IF (trim(var_name(i)).eq.                                  &
     &           trim(vname(1,idsbry(isvbar)))) THEN
          got_var(idsbry(isvbar))=.true.
          irp(ng)%Vid(idsbry(isvbar))=var_id(i)
#   ifdef SOLVE3D
        ELSE IF (trim(var_name(i)).eq.                                  &
     &           trim(vname(1,idsbry(isuvel)))) THEN
          got_var(idsbry(isuvel))=.true.
          irp(ng)%Vid(idsbry(isuvel))=var_id(i)
        ELSE IF (trim(var_name(i)).eq.                                  &
     &           trim(vname(1,idsbry(isvvel)))) THEN
          got_var(idsbry(isvvel))=.true.
          irp(ng)%Vid(idsbry(isvvel))=var_id(i)
#   endif
        END IF
#   ifdef SOLVE3D
        DO itrc=1,nt(ng)
          IF (trim(var_name(i)).eq.                                     &
     &        trim(vname(1,idsbry(istvar(itrc))))) THEN
            got_var(idsbry(istvar(itrc)))=.true.
            irp(ng)%Vid(idsbry(istvar(itrc)))=var_id(i)
          END IF
        END DO
#   endif
#  endif
#  ifdef ADJUST_WSTRESS
        IF (trim(var_name(i)).eq.trim(vname(1,idusms))) THEN
          got_var(idusms)=.true.
          irp(ng)%Vid(idusms)=var_id(i)
        ELSE IF (trim(var_name(i)).eq.trim(vname(1,idvsms))) THEN
          got_var(idvsms)=.true.
          irp(ng)%Vid(idvsms)=var_id(i)
        END IF
#  endif
#  if defined ADJUST_STFLUX && defined SOLVE3D
        DO itrc=1,nt(ng)
          IF (lstflux(itrc,ng)) THEN
            IF (trim(var_name(i)).eq.trim(vname(1,idtsur(itrc)))) THEN
              got_var(idtsur(itrc))=.true.
              irp(ng)%Vid(idtsur(itrc))=var_id(i)
            END IF
          END IF
        END DO
#  endif
      END DO
 
#  ifdef ADJUST_BOUNDARY
      IF (.not.got_var(idsbry(isfsur))) ldefine=.true.
      IF (.not.got_var(idsbry(isubar))) ldefine=.true.
      IF (.not.got_var(idsbry(isvbar))) ldefine=.true.
#   ifdef SOLVE3D
      IF (.not.got_var(idsbry(isuvel))) ldefine=.true.
      IF (.not.got_var(idsbry(isvvel))) ldefine=.true.
      DO itrc=1,nt(ng)
        IF (.not.got_var(idsbry(istvar(itrc)))) ldefine=.true.
      END DO
#   endif
#  endif
#  ifdef ADJUST_WSTRESS
      IF (.not.got_var(idusms)) ldefine=.true.
      IF (.not.got_var(idvsms)) ldefine=.true.
#  endif
#  if defined ADJUST_STFLUX && defined SOLVE3D
      DO itrc=1,nt(ng)
        IF (lstflux(itrc,ng)) THEN
          IF (.not.got_var(idtsur(itrc))) ldefine=.true.
        END IF
      END DO
#  endif
!
!  Put existing file into define mode so new variables can be added.
!
      IF (ldefine) THEN
        CALL netcdf_redef (ng, irpm, ncname, irp(ng)%ncid)
        IF (founderror(exit_flag, noerror, __line__, myfile)) THEN
          IF (master) WRITE (stdout,20) trim(ncname)
          RETURN
        END IF
      END IF
!
!-----------------------------------------------------------------------
!  Define the dimensions of staggered fields.
!-----------------------------------------------------------------------
!
      define: IF (ldefine) THEN
 
#  ifdef ADJUST_BOUNDARY
        got_iorj=.false.
        got_boundary=.false.
        got_obc_adjust=.false.
#  endif
#  if defined ADJUST_STFLUX || defined ADJUST_WSTRESS
        got_frc_adjust=.false.
#  endif
        DO i=1,n_dim
          SELECT CASE (trim(adjustl(dim_name(i))))
            CASE ('xi_rho')
              dimids( 1)=dim_id(i)
            CASE ('xi_u')
              dimids( 2)=dim_id(i)
            CASE ('xi_v')
              dimids( 3)=dim_id(i)
            CASE ('eta_rho')
              dimids( 5)=dim_id(i)
            CASE ('eta_u')
              dimids( 6)=dim_id(i)
            CASE ('eta_v')
              dimids( 7)=dim_id(i)
#  ifdef SOLVE3D
            CASE ('s_rho')
              dimids( 9)=dim_id(i)
            CASE ('s_w')
              dimids(10)=dim_id(i)
#  endif
#  ifdef ADJUST_BOUNDARY
            CASE ('boundary')
              dimids(14)=dim_id(i)
              got_boundary=.true.
            CASE ('IorJ')
              iorjdim=dim_id(i)
              got_iorj=.true.
#  endif
#  if defined ADJUST_STFLUX || defined ADJUST_WSTRESS
            CASE ('frc_adjust')
              frecdim=dim_id(i)
              got_frc_adjust=.true.
#  endif
#  ifdef ADJUST_BOUNDARY
            CASE ('obc_adjust')
              brecdim=dim_id(i)
              got_obc_adjust=.true.
#  endif
          END SELECT
        END DO
 
        dimids(12)=rec_id
#  ifdef ADJUST_BOUNDARY
        IF (.not.got_boundary) THEN
          status=def_dim(ng, irpm, irp(ng)%ncid, ncname, 'boundary',    &
     &                   4, dimids(14))
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END IF
        IF (.not.got_iorj) THEN
          status=def_dim(ng, irpm, irp(ng)%ncid, ncname, 'IorJ',        &
     &                   iobounds(ng)%IorJ, iorjdim)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END IF
#  endif
#  if defined ADJUST_STFLUX || defined ADJUST_WSTRESS
        IF (.not.got_frc_adjust) THEN
          status=def_dim(ng, irpm, irp(ng)%ncid, ncname, 'frc_adjust',  &
     &                   nfrec(ng), frecdim)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END IF
#  endif
#  ifdef ADJUST_BOUNDARY
        IF (.not.got_obc_adjust) THEN
          status=def_dim(ng, irpm, irp(ng)%ncid, ncname, 'obc_adjust',  &
     &                   nbrec(ng), brecdim)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END IF
#  endif
!
!  Define dimension vectors for staggered tracer type variables.
!
#  ifdef ADJUST_BOUNDARY
        t2dobc(1)=iorjdim
        t2dobc(2)=dimids(14)
        t2dobc(3)=brecdim
        t2dobc(4)=dimids(12)
#   ifdef SOLVE3D
        t3dobc(1)=iorjdim
        t3dobc(2)=dimids( 9)
        t3dobc(3)=dimids(14)
        t3dobc(4)=brecdim
        t3dobc(5)=dimids(12)
#   endif
#  endif
#  ifdef ADJUST_STFLUX
        t4dfrc(1)=dimids( 1)
        t4dfrc(2)=dimids( 5)
        t4dfrc(3)=frecdim
        t4dfrc(4)=dimids(12)
#  endif
#  ifdef ADJUST_WSTRESS
!
!  Define dimension vectors for staggered u-momentum type variables.
!
        u4dfrc(1)=dimids( 2)
        u4dfrc(2)=dimids( 6)
        u4dfrc(3)=frecdim
        u4dfrc(4)=dimids(12)
#  endif
#  ifdef ADJUST_WSTRESS
!
!  Define dimension vectors for staggered v-momentum type variables.
!
        v4dfrc(1)=dimids( 3)
        v4dfrc(2)=dimids( 7)
        v4dfrc(3)=frecdim
        v4dfrc(4)=dimids(12)
#  endif
!
!  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 additional variables. Notice that these variables have their
!  own fixed time-dimension to allow 4DVAR adjustments at other times
!  in addition to initialization time.
!-----------------------------------------------------------------------
 
#  ifdef ADJUST_BOUNDARY
!
!  Define free-surface open boundaries.
!
        IF (any(lobc(:,isfsur,ng))) THEN
          ifield=idsbry(isfsur)
          vinfo( 1)=vname(1,ifield)
          vinfo( 2)=vname(2,ifield)
          vinfo( 3)=vname(3,ifield)
          vinfo(14)=vname(4,ifield)
          vinfo(16)=vname(1,idtime)
          vinfo(21)=vname(6,ifield)
          aval(5)=real(iinfo(1,ifield,ng),r8)
          status=def_var(ng, irpm, irp(ng)%ncid, irp(ng)%Vid(ifield),   &
     &                   nf_fout, 4, t2dobc, aval, vinfo, ncname,       &
     &                   setfillval = .false.)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END IF
!
!  Define 2D U-momentum component open boundaries.
!
        IF (any(lobc(:,isubar,ng))) THEN
          ifield=idsbry(isubar)
          vinfo( 1)=vname(1,ifield)
          vinfo( 2)=vname(2,ifield)
          vinfo( 3)=vname(3,ifield)
          vinfo(14)=vname(4,ifield)
          vinfo(16)=vname(1,idtime)
          vinfo(21)=vname(6,ifield)
          aval(5)=real(iinfo(1,ifield,ng),r8)
          status=def_var(ng, irpm, irp(ng)%ncid, irp(ng)%Vid(ifield),   &
     &                   nf_fout, 4, t2dobc, aval, vinfo, ncname,       &
     &                   setfillval = .false.)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END IF
!
!  Define 2D V-momentum component open boundaries.
!
        IF (any(lobc(:,isvbar,ng))) THEN
          ifield=idsbry(isvbar)
          vinfo( 1)=vname(1,ifield)
          vinfo( 2)=vname(2,ifield)
          vinfo( 3)=vname(3,ifield)
          vinfo(14)=vname(4,ifield)
          vinfo(16)=vname(1,idtime)
          vinfo(21)=vname(6,ifield)
          aval(5)=real(iinfo(1,ifield,ng),r8)
          status=def_var(ng, irpm, irp(ng)%ncid, irp(ng)%Vid(ifield),   &
     &                   nf_fout, 4, t2dobc, aval, vinfo, ncname,       &
     &                   setfillval = .false.)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END IF
 
#   ifdef SOLVE3D
!
!  Define 3D U-momentum component open boundaries.
!
        IF (any(lobc(:,isuvel,ng))) THEN
          ifield=idsbry(isuvel)
          vinfo( 1)=vname(1,ifield)
          vinfo( 2)=vname(2,ifield)
          vinfo( 3)=vname(3,ifield)
          vinfo(14)=vname(4,ifield)
          vinfo(16)=vname(1,idtime)
          vinfo(21)=vname(6,ifield)
          aval(5)=real(iinfo(1,ifield,ng),r8)
          status=def_var(ng, irpm, irp(ng)%ncid, irp(ng)%Vid(ifield),   &
     &                   nf_fout, 5, t3dobc, aval, vinfo, ncname,       &
     &                   setfillval = .false.)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END IF
!
!  Define 3D V-momentum component open boundaries.
!
        IF (any(lobc(:,isvvel,ng))) THEN
          ifield=idsbry(isvvel)
          vinfo( 1)=vname(1,ifield)
          vinfo( 2)=vname(2,ifield)
          vinfo( 3)=vname(3,ifield)
          vinfo(14)=vname(4,ifield)
          vinfo(16)=vname(1,idtime)
          vinfo(21)=vname(6,ifield)
          aval(5)=real(iinfo(1,ifield,ng),r8)
          status=def_var(ng, irpm, irp(ng)%ncid, irp(ng)%Vid(ifield),   &
     &                   nf_fout, 5, t3dobc, aval, vinfo, ncname,       &
     &                   setfillval = .false.)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END IF
!
!  Define tracer type variables open boundaries.
!
        DO itrc=1,nt(ng)
          IF (any(lobc(:,istvar(itrc),ng))) THEN
            ifield=idsbry(istvar(itrc))
            vinfo( 1)=vname(1,ifield)
            vinfo( 2)=vname(2,ifield)
            vinfo( 3)=vname(3,ifield)
            vinfo(14)=vname(4,ifield)
            vinfo(16)=vname(1,idtime)
            vinfo(21)=vname(6,ifield)
            aval(5)=real(iinfo(1,ifield,ng),r8)
            status=def_var(ng, irpm, irp(ng)%ncid, irp(ng)%Vid(ifield), &
     &                     nf_fout, 5, t3dobc, aval, vinfo, ncname,     &
     &                     setfillval = .false.)
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
          END IF
        END DO
#   endif
#  endif
#  ifdef ADJUST_WSTRESS
!
!  Define surface U-momentum stress.
!
        IF (.not.got_var(idusms)) THEN
          vinfo( 1)=vname(1,idusms)
          vinfo( 2)=vname(2,idusms)
          vinfo( 3)=vname(3,idusms)
#   if defined WRITE_WATER && defined MASKING
          vinfo(20)='mask_u'
#   endif
          vinfo(21)=vname(6,idusms)
          vinfo(22)='coordinates'
          aval(5)=real(u2dvar,r8)
          status=def_var(ng, irpm, irp(ng)%ncid, irp(ng)%Vid(idusms),   &
     &                   nf_fout, 4, u4dfrc, aval, vinfo, ncname)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END IF
!
!  Define surface V-momentum stress.
!
        IF (.not.got_var(idvsms)) THEN
          vinfo( 1)=vname(1,idvsms)
          vinfo( 2)=vname(2,idvsms)
          vinfo( 3)=vname(3,idvsms)
#   if defined WRITE_WATER && defined MASKING
          vinfo(20)='mask_v'
#   endif
          vinfo(21)=vname(6,idvsms)
          vinfo(22)='coordinates'
          aval(5)=real(v2dvar,r8)
          status=def_var(ng, irpm, irp(ng)%ncid, irp(ng)%Vid(idvsms),   &
     &                   nf_fout, 4, v4dfrc, aval, vinfo, ncname)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END IF
#  endif
#  if defined ADJUST_STFLUX && defined SOLVE3D
!
!  Define surface tracer flux.
!
        DO itrc=1,nt(ng)
          IF (.not.got_var(idtsur(itrc)).and.lstflux(itrc,ng)) THEN
            vinfo( 1)=vname(1,idtsur(itrc))
            vinfo( 2)=trim(vname(2,idtsur(itrc)))
            vinfo( 3)=vname(3,idtsur(itrc))
            IF (itrc.eq.itemp) THEN
              vinfo(11)='upward flux, cooling'
              vinfo(12)='downward flux, heating'
            ELSE IF (itrc.eq.isalt) THEN
              vinfo(11)='upward flux, freshening (net precipitation)'
              vinfo(12)='downward flux, salting (net evaporation)'
            END IF
#   if defined WRITE_WATER && defined MASKING
            vinfo(20)='mask_rho'
#   endif
            vinfo(21)=vname(6,idtsur(itrc))
            vinfo(22)='coordinates'
            aval(5)=real(r2dvar,r8)
            status=def_var(ng, irpm, irp(ng)%ncid,                      &
     &                     irp(ng)%Vid(idtsur(itrc)), nf_fout,          &
     &                     4, t4dfrc, aval, vinfo, ncname)
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
          END IF
        END DO
#  endif
!
!-----------------------------------------------------------------------
!  Leave definition mode.
!-----------------------------------------------------------------------
!
        CALL netcdf_enddef (ng, irpm, ncname, irp(ng)%ncid)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
      END IF define
# endif
!
!=======================================================================
!  Open an existing tangent linear initial file, check its contents, and
!  prepare for appending data.
!=======================================================================
!
      IF (.not.ldefirp(ng)) THEN
        ncname=irp(ng)%name
 
# if !(defined ADJUST_STFLUX || defined ADJUST_WSTRESS)
!
!  Open representer model initialization file for read/write.
!
        CALL netcdf_open (ng, irpm, trim(ncname), 1, irp(ng)%ncid)
        IF (founderror(exit_flag, noerror, __line__, myfile)) THEN
          WRITE (stdout,10) trim(ncname)
          RETURN
        END IF
!
!  Inquire about the dimensions and check for consistency.
!
        CALL netcdf_check_dim (ng, irpm, ncname,                        &
     &                         ncid = irp(ng)%ncid)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Inquire about the variables.
!
        CALL netcdf_inq_var (ng, irpm, ncname,                          &
     &                         ncid = irp(ng)%ncid)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Initialize logical switches.
!
        DO i=1,nv
          got_var(i)=.false.
        END DO
# endif
!
!  Scan variable list from input NetCDF and activate switches for
!  initialization variables. Get variable IDs.
!
        DO i=1,n_var
          IF (trim(var_name(i)).eq.trim(vname(1,idtime))) THEN
            got_var(idtime)=.true.
            irp(ng)%Vid(idtime)=var_id(i)
          ELSE IF (trim(var_name(i)).eq.trim(vname(1,idfsur))) THEN
            got_var(idfsur)=.true.
            irp(ng)%Vid(idfsur)=var_id(i)
          ELSE IF (trim(var_name(i)).eq.trim(vname(1,idubar))) THEN
            got_var(idubar)=.true.
            irp(ng)%Vid(idubar)=var_id(i)
          ELSE IF (trim(var_name(i)).eq.trim(vname(1,idvbar))) THEN
            got_var(idvbar)=.true.
            irp(ng)%Vid(idvbar)=var_id(i)
# ifdef ADJUST_BOUNDARY
          ELSE IF (trim(var_name(i)).eq.                                &
     &             trim(vname(1,idsbry(isfsur)))) THEN
            got_var(idsbry(isfsur))=.true.
            irp(ng)%Vid(idsbry(isfsur))=var_id(i)
          ELSE IF (trim(var_name(i)).eq.                                &
     &             trim(vname(1,idsbry(isubar)))) THEN
            got_var(idsbry(isubar))=.true.
            irp(ng)%Vid(idsbry(isubar))=var_id(i)
          ELSE IF (trim(var_name(i)).eq.                                &
     &             trim(vname(1,idsbry(isvbar)))) THEN
            got_var(idsbry(isvbar))=.true.
            irp(ng)%Vid(idsbry(isvbar))=var_id(i)
# endif
# ifdef ADJUST_WSTRESS
          ELSE IF (trim(var_name(i)).eq.trim(vname(1,idusms))) THEN
            got_var(idusms)=.true.
            irp(ng)%Vid(idusms)=var_id(i)
          ELSE IF (trim(var_name(i)).eq.trim(vname(1,idvsms))) THEN
            got_var(idvsms)=.true.
            irp(ng)%Vid(idvsms)=var_id(i)
# endif
# ifdef SOLVE3D
          ELSE IF (trim(var_name(i)).eq.trim(vname(1,iduvel))) THEN
            got_var(iduvel)=.true.
            irp(ng)%Vid(iduvel)=var_id(i)
          ELSE IF (trim(var_name(i)).eq.trim(vname(1,idvvel))) THEN
            got_var(idvvel)=.true.
            irp(ng)%Vid(idvvel)=var_id(i)
#  ifdef ADJUST_BOUNDARY
          ELSE IF (trim(var_name(i)).eq.                                &
     &             trim(vname(1,idsbry(isuvel)))) THEN
            got_var(idsbry(isuvel))=.true.
            irp(ng)%Vid(idsbry(isuvel))=var_id(i)
          ELSE IF (trim(var_name(i)).eq.                                &
     &             trim(vname(1,idsbry(isvvel)))) THEN
            got_var(idsbry(isvvel))=.true.
            irp(ng)%Vid(idsbry(isvvel))=var_id(i)
#  endif
#  if defined BVF_MIXING || defined LMD_MIXING || \
      defined gls_mixing || defined my25_mixing
          ELSE IF (trim(var_name(i)).eq.trim(vname(1,idvvis))) THEN
            got_var(idvvis)=.true.
            irp(ng)%Vid(idvvis)=var_id(i)
          ELSE IF (trim(var_name(i)).eq.trim(vname(1,idtdif))) THEN
            got_var(idtdif)=.true.
            irp(ng)%Vid(idtdif)=var_id(i)
          ELSE IF (trim(var_name(i)).eq.trim(vname(1,idsdif))) THEN
            got_var(idsdif)=.true.
            irp(ng)%Vid(idsdif)=var_id(i)
#  endif
# 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.
              irp(ng)%Tid(itrc)=var_id(i)
#  ifdef ADJUST_BOUNDARY
            ELSE IF (trim(var_name(i)).eq.                              &
     &               trim(vname(1,idsbry(istvar(itrc))))) THEN
              got_var(idsbry(istvar(itrc)))=.true.
              irp(ng)%Vid(idsbry(istvar(itrc)))=var_id(i)
#  endif
#  ifdef ADJUST_STFLUX
            ELSE IF (trim(var_name(i)).eq.                              &
     &               trim(vname(1,idtsur(itrc)))) THEN
              got_var(idtsur(itrc))=.true.
              irp(ng)%Vid(idtsur(itrc))=var_id(i)
#  endif
            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,30) trim(vname(1,idtime)),          &
     &                                  trim(ncname)
          exit_flag=3
          RETURN
        END IF
        IF (.not.got_var(idfsur)) THEN
          IF (master) WRITE (stdout,30) trim(vname(1,idfsur)),          &
     &                                  trim(ncname)
          exit_flag=3
          RETURN
        END IF
        IF (.not.got_var(idubar)) THEN
          IF (master) WRITE (stdout,30) trim(vname(1,idubar)),          &
     &                                  trim(ncname)
          exit_flag=3
          RETURN
        END IF
        IF (.not.got_var(idvbar)) THEN
          IF (master) WRITE (stdout,30) trim(vname(1,idvbar)),          &
     &                                  trim(ncname)
          exit_flag=3
          RETURN
        END IF
# ifdef ADJUST_BOUNDARY
        IF (.not.got_var(idsbry(isfsur)).and.                           &
     &      any(lobc(:,isfsur,ng))) THEN
          IF (master) WRITE (stdout,30) trim(vname(1,idsbry(isfsur))),  &
     &                                  trim(ncname)
          exit_flag=3
          RETURN
        END IF
        IF (.not.got_var(idsbry(isubar)).and.                           &
     &      any(lobc(:,isubar,ng))) THEN
          IF (master) WRITE (stdout,30) trim(vname(1,idsbry(isubar))),  &
     &                                  trim(ncname)
          exit_flag=3
          RETURN
        END IF
        IF (.not.got_var(idsbry(isvbar)).and.                           &
     &      any(lobc(:,isvbar,ng))) THEN
          IF (master) WRITE (stdout,30) trim(vname(1,idsbry(isvbar))),  &
     &                                  trim(ncname)
          exit_flag=3
          RETURN
        END IF
# endif
# ifdef ADJUST_WSTRESS
        IF (.not.got_var(idusms)) THEN
          IF (master) WRITE (stdout,30) trim(vname(1,idusms)),          &
     &                                  trim(ncname)
          exit_flag=3
          RETURN
        END IF
        IF (.not.got_var(idvsms)) THEN
          IF (master) WRITE (stdout,30) trim(vname(1,idvsms)),          &
     &                                  trim(ncname)
          exit_flag=3
          RETURN
        END IF
# endif
# ifdef SOLVE3D
        IF (.not.got_var(iduvel)) THEN
          IF (master) WRITE (stdout,30) trim(vname(1,iduvel)),          &
     &                                  trim(ncname)
          exit_flag=3
          RETURN
        END IF
        IF (.not.got_var(idvvel)) THEN
          IF (master) WRITE (stdout,30) trim(vname(1,idvvel)),          &
     &                                  trim(ncname)
          exit_flag=3
          RETURN
        END IF
#  ifdef ADJUST_BOUNDARY
        IF (.not.got_var(idsbry(isuvel)).and.                           &
     &      any(lobc(:,isuvel,ng))) THEN
          IF (master) WRITE (stdout,30) trim(vname(1,idsbry(isuvel))),  &
     &                                  trim(ncname)
          exit_flag=3
          RETURN
        END IF
        IF (.not.got_var(idsbry(isvvel)).and.                           &
     &      any(lobc(:,isvvel,ng))) THEN
          IF (master) WRITE (stdout,30) trim(vname(1,idsbry(isvvel))),  &
     &                                  trim(ncname)
          exit_flag=3
          RETURN
        END IF
#  endif
#  ifdef ADJUST_BOUNDARY
        IF (.not.got_var(idsbry(isuvel)).and.                           &
     &      any(lobc(:,isuvel,ng))) THEN
          IF (master) WRITE (stdout,30) trim(vname(1,idsbry(isuvel))),  &
     &                                  trim(ncname)
          exit_flag=3
          RETURN
        END IF
        IF (.not.got_var(idsbry(isvvel)).and.                           &
     &      any(lobc(:,isvvel,ng))) THEN
          IF (master) WRITE (stdout,30) trim(vname(1,idsbry(isvvel))),  &
     &                                  trim(ncname)
          exit_flag=3
          RETURN
        END IF
#  endif
# endif
# ifdef SOLVE3D
        DO itrc=1,nt(ng)
          IF (.not.got_var(idtvar(itrc))) THEN
            IF (master) WRITE (stdout,30) trim(vname(1,idtvar(itrc))),  &
     &                                    trim(ncname)
            exit_flag=3
            RETURN
          END IF
#  ifdef ADJUST_BOUNDARY
          IF (.not.got_var(idsbry(istvar(itrc))).and.                   &
     &        any(lobc(:,istvar(itrc),ng))) THEN
            IF (master) WRITE (stdout,30)                               &
     &                        trim(vname(1,idsbry(istvar(itrc)))),      &
     &                        trim(ncname)
            exit_flag=3
            RETURN
          END IF
#  endif
#  ifdef ADJUST_STFLUX
          IF (.not.got_var(idtsur(itrc)).and.lstflux(itrc,ng)) THEN
            IF (master) WRITE (stdout,30) trim(vname(1,idtsur(itrc))),  &
     &                                    trim(ncname)
            exit_flag=3
            RETURN
          END IF
#  endif
        END DO
# endif
!
!  Set unlimited time record dimension to the appropriate value.
!
        irp(ng)%Rindex=rec_size
        fcount=irp(ng)%Fcount
        irp(ng)%Nrec(fcount)=rec_size
      END IF
!
  10  FORMAT (/,' RP_DEF_INI_NF90 - unable to open initial NetCDF',     &
     &        ' file: ',a)
  20  FORMAT (/,' RP_DEF_INI_NF90 - unable to put in define mode',      &
     &        ' initial NetCDF file: ',a)
  30  FORMAT (/,' RP_DEF_INI_NF90 - unable to find variable: ',a,2x,    &
     &        ' in file: ',a)
!
      RETURN
      END SUBROUTINE rp_def_ini_nf90
 
# if defined PIO_LIB && defined DISTRIBUTE
!
!***********************************************************************
      SUBROUTINE rp_def_ini_pio (ng)
!***********************************************************************
!
      USE mod_pio_netcdf
!
!  Imported variable declarations.
!
      integer, intent(in) :: ng
!
!  Local variable declarations.
!
      logical :: ldefine = .false.
#  ifdef ADJUST_BOUNDARY
      logical :: got_iorj, got_boundary, got_obc_adjust
#  endif
#  if defined ADJUST_STFLUX || defined ADJUST_WSTRESS
      logical :: got_frc_adjust
#  endif
      logical :: got_var(nv)
!
      integer, parameter :: natt = 25
 
      integer :: i, j, ifield, itrc, status
      integer :: fcount
#  ifdef ADJUST_BOUNDARY
      integer :: iorjdim, brecdim
#  endif
#  if defined ADJUST_STFLUX || defined ADJUST_WSTRESS
      integer :: frecdim
#  endif
      integer :: dimids(ndimid)
#  ifdef ADJUST_BOUNDARY
      integer :: t2dobc(4)
#   ifdef SOLVE3D
      integer :: t3dobc(5)
#   endif
#  endif
#  if defined ADJUST_STFLUX || defined ADJUST_WSTRESS
      integer :: t4dfrc(4), u4dfrc(4), v4dfrc(4)
#  endif
!
      real(r8) :: aval(6)
!
      character (len=256)    :: ncname
      character (len=MaxLen) :: vinfo(natt)
 
      character (len=*), parameter :: myfile =                          &
     &  __FILE__//", rp_def_ini_pio"
!
      sourcefile=myfile
 
#  if defined ADJUST_BOUNDARY || \
     defined adjust_stflux   || defined adjust_wstress
!
!=======================================================================
!  Open existing representer model initial conditions and define new
!  variables.
!=======================================================================
!
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      ncname=irp(ng)%name
!
!  Open representer model initialization file for read/write.
!
      CALL pio_netcdf_open (ng, irpm, ncname, 1, irp(ng)%pioFile)
      IF (founderror(exit_flag, noerror, __line__, myfile)) THEN
        IF (master) WRITE (stdout,10) trim(ncname)
        RETURN
      END IF
!
!  Inquire about the dimensions and check for consistency.
!
      CALL pio_netcdf_check_dim (ng, irpm, ncname,                      &
     &                           piofile = irp(ng)%pioFile)
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Inquire about the variables.
!
      CALL pio_netcdf_inq_var (ng, irpm, ncname,                        &
     &                         piofile = irp(ng)%pioFile)
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Check if surface forcing variables have been already defined.
!
      DO i=1,nv
        got_var(i)=.false.
      END DO
!
      DO i=1,n_var
#   ifdef ADJUST_BOUNDARY
        IF (trim(var_name(i)).eq.                                       &
     &      trim(vname(1,idsbry(isfsur)))) THEN
          got_var(idsbry(isfsur))=.true.
          irp(ng)%Vid(idsbry(isfsur))=var_id(i)
        ELSE IF (trim(var_name(i)).eq.                                  &
     &           trim(vname(1,idsbry(isubar)))) THEN
          got_var(idsbry(isubar))=.true.
          irp(ng)%Vid(idsbry(isubar))=var_id(i)
        ELSE IF (trim(var_name(i)).eq.                                  &
     &           trim(vname(1,idsbry(isvbar)))) THEN
          got_var(idsbry(isvbar))=.true.
          irp(ng)%Vid(idsbry(isvbar))=var_id(i)
#    ifdef SOLVE3D
        ELSE IF (trim(var_name(i)).eq.                                  &
     &           trim(vname(1,idsbry(isuvel)))) THEN
          got_var(idsbry(isuvel))=.true.
          irp(ng)%Vid(idsbry(isuvel))=var_id(i)
        ELSE IF (trim(var_name(i)).eq.                                  &
     &           trim(vname(1,idsbry(isvvel)))) THEN
          got_var(idsbry(isvvel))=.true.
          irp(ng)%Vid(idsbry(isvvel))=var_id(i)
#    endif
        END IF
#    ifdef SOLVE3D
        DO itrc=1,nt(ng)
          IF (trim(var_name(i)).eq.                                     &
     &        trim(vname(1,idsbry(istvar(itrc))))) THEN
            got_var(idsbry(istvar(itrc)))=.true.
            irp(ng)%Vid(idsbry(istvar(itrc)))=var_id(i)
          END IF
        END DO
#    endif
#   endif
#   ifdef ADJUST_WSTRESS
        IF (trim(var_name(i)).eq.trim(vname(1,idusms))) THEN
          got_var(idusms)=.true.
          irp(ng)%Vid(idusms)=var_id(i)
        ELSE IF (trim(var_name(i)).eq.trim(vname(1,idvsms))) THEN
          got_var(idvsms)=.true.
          irp(ng)%Vid(idvsms)=var_id(i)
        END IF
#   endif
#   if defined ADJUST_STFLUX && defined SOLVE3D
        DO itrc=1,nt(ng)
          IF (lstflux(itrc,ng)) THEN
            IF (trim(var_name(i)).eq.trim(vname(1,idtsur(itrc)))) THEN
              got_var(idtsur(itrc))=.true.
              irp(ng)%Vid(idtsur(itrc))=var_id(i)
            END IF
          END IF
        END DO
#   endif
      END DO
 
#   ifdef ADJUST_BOUNDARY
      IF (.not.got_var(idsbry(isfsur))) ldefine=.true.
      IF (.not.got_var(idsbry(isubar))) ldefine=.true.
      IF (.not.got_var(idsbry(isvbar))) ldefine=.true.
#    ifdef SOLVE3D
      IF (.not.got_var(idsbry(isuvel))) ldefine=.true.
      IF (.not.got_var(idsbry(isvvel))) ldefine=.true.
      DO itrc=1,nt(ng)
        IF (.not.got_var(idsbry(istvar(itrc)))) ldefine=.true.
      END DO
#    endif
#   endif
#   ifdef ADJUST_WSTRESS
      IF (.not.got_var(idusms)) ldefine=.true.
      IF (.not.got_var(idvsms)) ldefine=.true.
#   endif
#   if defined ADJUST_STFLUX && defined SOLVE3D
      DO itrc=1,nt(ng)
        IF (lstflux(itrc,ng)) THEN
          IF (.not.got_var(idtsur(itrc))) ldefine=.true.
        END IF
      END DO
#   endif
!
!  Put existing file into define mode so new variables can be added.
!
      IF (ldefine) THEN
        CALL pio_netcdf_redef (ng, irpm, ncname, irp(ng)%pioFile)
        IF (founderror(exit_flag, noerror, __line__, myfile)) THEN
          IF (master) WRITE (stdout,20) trim(ncname)
          RETURN
        END IF
      END IF
!
!-----------------------------------------------------------------------
!  Define the dimensions of staggered fields.
!-----------------------------------------------------------------------
!
      define: IF (ldefine) THEN
 
#   ifdef ADJUST_BOUNDARY
        got_iorj=.false.
        got_boundary=.false.
        got_obc_adjust=.false.
#   endif
#   if defined ADJUST_STFLUX || defined ADJUST_WSTRESS
        got_frc_adjust=.false.
#   endif
        DO i=1,n_dim
          SELECT CASE (trim(adjustl(dim_name(i))))
            CASE ('xi_rho')
              dimids( 1)=dim_id(i)
            CASE ('xi_u')
              dimids( 2)=dim_id(i)
            CASE ('xi_v')
              dimids( 3)=dim_id(i)
            CASE ('eta_rho')
              dimids( 5)=dim_id(i)
            CASE ('eta_u')
              dimids( 6)=dim_id(i)
            CASE ('eta_v')
              dimids( 7)=dim_id(i)
#   ifdef SOLVE3D
            CASE ('s_rho')
              dimids( 9)=dim_id(i)
            CASE ('s_w')
              dimids(10)=dim_id(i)
#   endif
#   ifdef ADJUST_BOUNDARY
            CASE ('boundary')
              dimids(14)=dim_id(i)
              got_boundary=.true.
            CASE ('IorJ')
              iorjdim=dim_id(i)
              got_iorj=.true.
#   endif
#   if defined ADJUST_STFLUX || defined ADJUST_WSTRESS
            CASE ('frc_adjust')
              frecdim=dim_id(i)
              got_frc_adjust=.true.
#   endif
#   ifdef ADJUST_BOUNDARY
            CASE ('obc_adjust')
              brecdim=dim_id(i)
              got_obc_adjust=.true.
#   endif
          END SELECT
        END DO
!
        dimids(12)=rec_id
#   ifdef ADJUST_BOUNDARY
        IF (.not.got_boundary) THEN
          status=def_dim(ng, irpm, irp(ng)%pioFile, ncname,             &
     &                   'boundary', 4, dimids(14))
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END IF
        IF (.not.got_iorj) THEN
          status=def_dim(ng, irpm, irp(ng)%pioFile, ncname,             &
     &                   'IorJ', iobounds(ng)%IorJ, iorjdim)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END IF
#   endif
#   if defined ADJUST_STFLUX || defined ADJUST_WSTRESS
        IF (.not.got_frc_adjust) THEN
          status=def_dim(ng, irpm, irp(ng)%pioFile, ncname,             &
     &                   'frc_adjust', nfrec(ng), frecdim)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END IF
#   endif
#   ifdef ADJUST_BOUNDARY
        IF (.not.got_obc_adjust) THEN
          status=def_dim(ng, irpm, irp(ng)%pioFile, ncname,             &
     &                   'obc_adjust', nbrec(ng), brecdim)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END IF
#   endif
!
!  Define dimension vectors for staggered tracer type variables.
!
#   ifdef ADJUST_BOUNDARY
        t2dobc(1)=iorjdim
        t2dobc(2)=dimids(14)
        t2dobc(3)=brecdim
        t2dobc(4)=dimids(12)
#    ifdef SOLVE3D
        t3dobc(1)=iorjdim
        t3dobc(2)=dimids( 9)
        t3dobc(3)=dimids(14)
        t3dobc(4)=brecdim
        t3dobc(5)=dimids(12)
#    endif
#   endif
#   ifdef ADJUST_STFLUX
        t4dfrc(1)=dimids( 1)
        t4dfrc(2)=dimids( 5)
        t4dfrc(3)=frecdim
        t4dfrc(4)=dimids(12)
#   endif
#   ifdef ADJUST_WSTRESS
!
!  Define dimension vectors for staggered u-momentum type variables.
!
        u4dfrc(1)=dimids( 2)
        u4dfrc(2)=dimids( 6)
        u4dfrc(3)=frecdim
        u4dfrc(4)=dimids(12)
#   endif
#   ifdef ADJUST_WSTRESS
!
!  Define dimension vectors for staggered v-momentum type variables.
!
        v4dfrc(1)=dimids( 3)
        v4dfrc(2)=dimids( 7)
        v4dfrc(3)=frecdim
        v4dfrc(4)=dimids(12)
#   endif
!
!  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 additional variables. Notice that these variables have their
!  own fixed time-dimension to allow 4DVAR adjustments at other times
!  in addition to initialization time.
!-----------------------------------------------------------------------
 
#   ifdef ADJUST_BOUNDARY
!
!  Define free-surface open boundaries.
!
        IF (any(lobc(:,isfsur,ng))) THEN
          ifield=idsbry(isfsur)
          vinfo( 1)=vname(1,ifield)
          vinfo( 2)=vname(2,ifield)
          vinfo( 3)=vname(3,ifield)
          vinfo(14)=vname(4,ifield)
          vinfo(16)=vname(1,idtime)
          vinfo(21)=vname(6,ifield)
          aval(5)=real(iinfo(1,ifield,ng),r8)
          irp(ng)%pioVar(ifield)%dkind=pio_fout
          irp(ng)%pioVar(ifield)%gtype=r2dobc
!
          status=def_var(ng, irpm, irp(ng)%pioFile,                     &
     &                   irp(ng)%pioVar(ifield)%vd, pio_fout,           &
     &                   4, t2dobc, aval, vinfo, ncname,                &
     &                   setfillval = .false.)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END IF
!
!  Define 2D U-momentum component open boundaries.
!
        IF (any(lobc(:,isubar,ng))) THEN
          ifield=idsbry(isubar)
          vinfo( 1)=vname(1,ifield)
          vinfo( 2)=vname(2,ifield)
          vinfo( 3)=vname(3,ifield)
          vinfo(14)=vname(4,ifield)
          vinfo(16)=vname(1,idtime)
          vinfo(21)=vname(6,ifield)
          aval(5)=real(iinfo(1,ifield,ng),r8)
          irp(ng)%pioVar(ifield)%dkind=pio_fout
          irp(ng)%pioVar(ifield)%gtype=u2dobc
!
          status=def_var(ng, irpm, irp(ng)%pioFile,                     &
     &                   irp(ng)%pioVar(ifield)%vd, pio_fout,           &
     &                   4, t2dobc, aval, vinfo, ncname,                &
     &                   setfillval = .false.)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END IF
!
!  Define 2D V-momentum component open boundaries.
!
        IF (any(lobc(:,isvbar,ng))) THEN
          ifield=idsbry(isvbar)
          vinfo( 1)=vname(1,ifield)
          vinfo( 2)=vname(2,ifield)
          vinfo( 3)=vname(3,ifield)
          vinfo(14)=vname(4,ifield)
          vinfo(16)=vname(1,idtime)
          vinfo(21)=vname(6,ifield)
          aval(5)=real(iinfo(1,ifield,ng),r8)
          irp(ng)%pioVar(ifield)%dkind=pio_fout
          irp(ng)%pioVar(ifield)%gtype=v2dobc
!
          status=def_var(ng, irpm, irp(ng)%pioFile,                     &
     &                   irp(ng)%pioVar(ifield)%vd, pio_fout,           &
     &                   4, t2dobc, aval, vinfo, ncname,                &
     &                   setfillval = .false.)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END IF
 
#    ifdef SOLVE3D
!
!  Define 3D U-momentum component open boundaries.
!
        IF (any(lobc(:,isuvel,ng))) THEN
          ifield=idsbry(isuvel)
          vinfo( 1)=vname(1,ifield)
          vinfo( 2)=vname(2,ifield)
          vinfo( 3)=vname(3,ifield)
          vinfo(14)=vname(4,ifield)
          vinfo(16)=vname(1,idtime)
          vinfo(21)=vname(6,ifield)
          aval(5)=real(iinfo(1,ifield,ng),r8)
          irp(ng)%pioVar(ifield)%dkind=pio_fout
          irp(ng)%pioVar(ifield)%gtype=u3dobc
!
          status=def_var(ng, irpm, irp(ng)%pioFile,                     &
     &                   irp(ng)%pioVar(ifield)%vd, pio_fout,           &
     &                   5, t3dobc, aval, vinfo, ncname,                &
     &                   setfillval = .false.)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END IF
!
!  Define 3D V-momentum component open boundaries.
!
        IF (any(lobc(:,isvvel,ng))) THEN
          ifield=idsbry(isvvel)
          vinfo( 1)=vname(1,ifield)
          vinfo( 2)=vname(2,ifield)
          vinfo( 3)=vname(3,ifield)
          vinfo(14)=vname(4,ifield)
          vinfo(16)=vname(1,idtime)
          vinfo(21)=vname(6,ifield)
          aval(5)=real(iinfo(1,ifield,ng),r8)
          irp(ng)%pioVar(ifield)%dkind=pio_fout
          irp(ng)%pioVar(ifield)%gtype=v3dobc
!
          status=def_var(ng, irpm, irp(ng)%pioFile,                     &
     &                   irp(ng)%pioVar(ifield)%vd, pio_fout,           &
     &                   5, t3dobc, aval, vinfo, ncname,                &
     &                   setfillval = .false.)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END IF
!
!  Define tracer type variables open boundaries.
!
        DO itrc=1,nt(ng)
          IF (any(lobc(:,istvar(itrc),ng))) THEN
            ifield=idsbry(istvar(itrc))
            vinfo( 1)=vname(1,ifield)
            vinfo( 2)=vname(2,ifield)
            vinfo( 3)=vname(3,ifield)
            vinfo(14)=vname(4,ifield)
            vinfo(16)=vname(1,idtime)
            vinfo(21)=vname(6,ifield)
            aval(5)=real(iinfo(1,ifield,ng),r8)
            irp(ng)%pioVar(ifield)%dkind=pio_fout
            irp(ng)%pioVar(ifield)%gtype=r3dobc
!
            status=def_var(ng, irpm, irp(ng)%pioFile,                   &
     &                     irp(ng)%pioVar(ifield)%vd, pio_fout,         &
     &                     5, t3dobc, aval, vinfo, ncname,              &
     &                     setfillval = .false.)
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
          END IF
        END DO
#    endif
#   endif
#   ifdef ADJUST_WSTRESS
!
!  Define surface U-momentum stress.
!
        IF (.not.got_var(idusms)) THEN
          vinfo( 1)=vname(1,idusms)
          vinfo( 2)=vname(2,idusms)
          vinfo( 3)=vname(3,idusms)
#    if defined WRITE_WATER && defined MASKING
          vinfo(20)='mask_u'
#    endif
          vinfo(21)=vname(6,idusms)
          vinfo(22)='coordinates'
          aval(5)=real(u2dvar,r8)
          irp(ng)%pioVar(idusms)%dkind=pio_fout
          irp(ng)%pioVar(idusms)%gtype=u2dvar
!
          status=def_var(ng, irpm, irp(ng)%pioFile,                     &
     &                   irp(ng)%pioVar(idusms)%vd, pio_fout,           &
     &                   4, u4dfrc, aval, vinfo, ncname)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END IF
!
!  Define surface V-momentum stress.
!
        IF (.not.got_var(idvsms)) THEN
          vinfo( 1)=vname(1,idvsms)
          vinfo( 2)=vname(2,idvsms)
          vinfo( 3)=vname(3,idvsms)
#    if defined WRITE_WATER && defined MASKING
          vinfo(20)='mask_v'
#    endif
          vinfo(21)=vname(6,idvsms)
          vinfo(22)='coordinates'
          aval(5)=real(v2dvar,r8)
          irp(ng)%pioVar(idvsms)%dkind=pio_fout
          irp(ng)%pioVar(idvsms)%gtype=v2dvar
!
          status=def_var(ng, irpm, irp(ng)%pioFile,                     &
     &                   irp(ng)%pioVar(idvsms)%vd, pio_fout,           &
     &                   4, v4dfrc, aval, vinfo, ncname)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END IF
#   endif
#   if defined ADJUST_STFLUX && defined SOLVE3D
!
!  Define surface tracer flux.
!
        DO itrc=1,nt(ng)
          IF (.not.got_var(idtsur(itrc)).and.lstflux(itrc,ng)) THEN
            vinfo( 1)=vname(1,idtsur(itrc))
            vinfo( 2)=trim(vname(2,idtsur(itrc)))
            vinfo( 3)=vname(3,idtsur(itrc))
            IF (itrc.eq.itemp) THEN
              vinfo(11)='upward flux, cooling'
              vinfo(12)='downward flux, heating'
            ELSE IF (itrc.eq.isalt) THEN
              vinfo(11)='upward flux, freshening (net precipitation)'
              vinfo(12)='downward flux, salting (net evaporation)'
            END IF
#    if defined WRITE_WATER && defined MASKING
            vinfo(20)='mask_rho'
#    endif
            vinfo(21)=vname(6,idtsur(itrc))
            vinfo(22)='coordinates'
            aval(5)=real(r2dvar,r8)
            irp(ng)%pioVar(idtsur(itrc))%dkind=pio_fout
            irp(ng)%pioVar(idtsur(itrc))%gtype=r2dvar
!
            status=def_var(ng, irpm, irp(ng)%pioFile,                   &
     &                     irp(ng)%pioVar(idtsur(itrc))%vd, pio_fout,   &
     &                     4, t4dfrc, aval, vinfo, ncname)
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
          END IF
        END DO
#   endif
!
!-----------------------------------------------------------------------
!  Leave definition mode.
!-----------------------------------------------------------------------
!
        CALL pio_netcdf_enddef (ng, irpm, ncname, irp(ng)%pioFile)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
      END IF define
#  endif
!
!=======================================================================
!  Open an existing tangent linear initial file, check its contents, and
!  prepare for appending data.
!=======================================================================
!
      IF (.not.ldefirp(ng)) THEN
        ncname=irp(ng)%name
 
#  if !(defined ADJUST_STFLUX || defined ADJUST_WSTRESS)
!
!  Open representer model initialization file for read/write.
!
        CALL pio_netcdf_open (ng, irpm, ncname, 1, irp(ng)%pioFile)
        IF (founderror(exit_flag, noerror, __line__, myfile)) THEN
          WRITE (stdout,10) trim(ncname)
          RETURN
        END IF
!
!  Inquire about the dimensions and check for consistency.
!
        CALL pio_netcdf_check_dim (ng, irpm, ncname,                    &
     &                             piofile = irp(ng)%pioFile)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Inquire about the variables.
!
        CALL pio_netcdf_inq_var (ng, irpm, ncname,                      &
     &                           piofile = irp(ng)%pioFile)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Initialize logical switches.
!
        DO i=1,nv
          got_var(i)=.false.
        END DO
#  endif
!
!  Scan variable list from input NetCDF and activate switches for
!  initialization variables. Get variable IDs.
!
        DO i=1,n_var
          IF (trim(var_name(i)).eq.trim(vname(1,idtime))) THEN
            got_var(idtime)=.true.
            irp(ng)%pioVar(idtime)%vd=var_desc(i)
            irp(ng)%pioVar(idtime)%dkind=pio_tout
            irp(ng)%pioVar(idtime)%gtype=0
          ELSE IF (trim(var_name(i)).eq.trim(vname(1,idfsur))) THEN
            got_var(idfsur)=.true.
            irp(ng)%pioVar(idfsur)%vd=var_desc(i)
            irp(ng)%pioVar(idfsur)%dkind=pio_fout
            irp(ng)%pioVar(idfsur)%gtype=r2dvar
          ELSE IF (trim(var_name(i)).eq.trim(vname(1,idubar))) THEN
            got_var(idubar)=.true.
            irp(ng)%pioVar(idubar)%vd=var_desc(i)
            irp(ng)%pioVar(idubar)%dkind=pio_fout
            irp(ng)%pioVar(idubar)%gtype=u2dvar
          ELSE IF (trim(var_name(i)).eq.trim(vname(1,idvbar))) THEN
            got_var(idvbar)=.true.
            irp(ng)%pioVar(idvbar)%vd=var_desc(i)
            irp(ng)%pioVar(idvbar)%dkind=pio_fout
            irp(ng)%pioVar(idvbar)%gtype=v2dvar
#  ifdef ADJUST_BOUNDARY
          ELSE IF (trim(var_name(i)).eq.                                &
     &             trim(vname(1,idsbry(isfsur)))) THEN
            got_var(idsbry(isfsur))=.true.
            irp(ng)%pioVar(idsbry(isfsur))%vd=var_desc(i)
            irp(ng)%pioVar(idsbry(isfsur))%dkind=pio_fout
            irp(ng)%pioVar(idsbry(isfsur))%gtype=r2dobc
          ELSE IF (trim(var_name(i)).eq.                                &
     &             trim(vname(1,idsbry(isubar)))) THEN
            got_var(idsbry(isubar))=.true.
            irp(ng)%pioVar(idsbry(isubar))%vd=var_desc(i)
            irp(ng)%pioVar(idsbry(isubar))%dkind=pio_fout
            irp(ng)%pioVar(idsbry(isubar))%gtype=u2dobc
          ELSE IF (trim(var_name(i)).eq.                                &
     &             trim(vname(1,idsbry(isvbar)))) THEN
            got_var(idsbry(isvbar))=.true.
            irp(ng)%pioVar(idsbry(isvbar))%vd=var_desc(i)
            irp(ng)%pioVar(idsbry(isvbar))%dkind=pio_fout
            irp(ng)%pioVar(idsbry(isvbar))%gtype=v2dobc
#  endif
#  ifdef ADJUST_WSTRESS
          ELSE IF (trim(var_name(i)).eq.trim(vname(1,idusms))) THEN
            got_var(idusms)=.true.
            irp(ng)%pioVar(idusms)%vd=var_desc(i)
            irp(ng)%pioVar(idusms)%dkind=pio_fout
            irp(ng)%pioVar(idusms)%gtype=u2dvar
          ELSE IF (trim(var_name(i)).eq.trim(vname(1,idvsms))) THEN
            got_var(idvsms)=.true.
            irp(ng)%pioVar(idvsms)%vd=var_desc(i)
            irp(ng)%pioVar(idvsms)%dkind=pio_fout
            irp(ng)%pioVar(idvsms)%gtype=v2dvar
#  endif
#  ifdef SOLVE3D
          ELSE IF (trim(var_name(i)).eq.trim(vname(1,iduvel))) THEN
            got_var(iduvel)=.true.
            irp(ng)%pioVar(iduvel)%vd=var_desc(i)
            irp(ng)%pioVar(iduvel)%dkind=pio_fout
            irp(ng)%pioVar(iduvel)%gtype=u3dvar
          ELSE IF (trim(var_name(i)).eq.trim(vname(1,idvvel))) THEN
            got_var(idvvel)=.true.
            irp(ng)%pioVar(idvvel)%vd=var_desc(i)
            irp(ng)%pioVar(idvvel)%dkind=pio_fout
            irp(ng)%pioVar(idvvel)%gtype=v3dvar
#   ifdef ADJUST_BOUNDARY
          ELSE IF (trim(var_name(i)).eq.                                &
     &             trim(vname(1,idsbry(isuvel)))) THEN
            got_var(idsbry(isuvel))=.true.
            irp(ng)%pioVar(idsbry(isuvel))%vd=var_desc(i)
            irp(ng)%pioVar(idsbry(isuvel))%dkind=pio_fout
            irp(ng)%pioVar(idsbry(isuvel))%gtype=u3dobc
          ELSE IF (trim(var_name(i)).eq.                                &
     &             trim(vname(1,idsbry(isvvel)))) THEN
            got_var(idsbry(isvvel))=.true.
            irp(ng)%pioVar(idsbry(isvvel))%vd=var_desc(i)
            irp(ng)%pioVar(idsbry(isvvel))%dkind=pio_fout
            irp(ng)%pioVar(idsbry(isvvel))%gtype=v3dobc
#   endif
#   if defined BVF_MIXING || defined LMD_MIXING || \
      defined gls_mixing || defined my25_mixing
          ELSE IF (trim(var_name(i)).eq.trim(vname(1,idvvis))) THEN
            got_var(idvvis)=.true.
            irp(ng)%pioVar(idvvis)%vd=var_desc(i)
            irp(ng)%pioVar(idvvis)%dkind=pio_fout
            irp(ng)%pioVar(idvvis)%gtype=w3dvar
          ELSE IF (trim(var_name(i)).eq.trim(vname(1,idtdif))) THEN
            got_var(idtdif)=.true.
            irp(ng)%pioVar(idtdif)%vd=var_desc(i)
            irp(ng)%pioVar(idtdif)%dkind=pio_fout
            irp(ng)%pioVar(idtdif)%gtype=w3dvar
          ELSE IF (trim(var_name(i)).eq.trim(vname(1,idsdif))) THEN
            got_var(idsdif)=.true.
            irp(ng)%pioVar(idsdif)%vd=var_desc(i)
            irp(ng)%pioVar(idsdif)%dkind=pio_fout
            irp(ng)%pioVar(idsdif)%gtype=w3dvar
#   endif
#  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.
              irp(ng)%pioTrc(itrc)%vd=var_desc(i)
              irp(ng)%pioTrc(itrc)%dkind=pio_fout
              irp(ng)%pioTrc(itrc)%gtype=r3dvar
#   ifdef ADJUST_BOUNDARY
            ELSE IF (trim(var_name(i)).eq.                              &
     &               trim(vname(1,idsbry(istvar(itrc))))) THEN
              got_var(idsbry(istvar(itrc)))=.true.
              irp(ng)%pioVar(idsbry(istvar(itrc)))%vd=var_desc(i)
              irp(ng)%pioVar(idsbry(istvar(itrc)))%dkind=pio_fout
              irp(ng)%pioVar(idsbry(istvar(itrc)))%gtype=r3dobc
#   endif
#   ifdef ADJUST_STFLUX
            ELSE IF (trim(var_name(i)).eq.                              &
     &               trim(vname(1,idtsur(itrc)))) THEN
              got_var(idtsur(itrc))=.true.
              irp(ng)%pioVar(idtsur(itrc))%vd=var_desc(i)
              irp(ng)%pioVar(idtsur(itrc))%dkind=pio_fout
              irp(ng)%pioVar(idtsur(itrc))%gtype=r2dvar
#   endif
            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,30) trim(vname(1,idtime)),          &
     &                                  trim(ncname)
          exit_flag=3
          RETURN
        END IF
        IF (.not.got_var(idfsur)) THEN
          IF (master) WRITE (stdout,30) trim(vname(1,idfsur)),          &
     &                                  trim(ncname)
          exit_flag=3
          RETURN
        END IF
        IF (.not.got_var(idubar)) THEN
          IF (master) WRITE (stdout,30) trim(vname(1,idubar)),          &
     &                                  trim(ncname)
          exit_flag=3
          RETURN
        END IF
        IF (.not.got_var(idvbar)) THEN
          IF (master) WRITE (stdout,30) trim(vname(1,idvbar)),          &
     &                                  trim(ncname)
          exit_flag=3
          RETURN
        END IF
#  ifdef ADJUST_BOUNDARY
        IF (.not.got_var(idsbry(isfsur)).and.                           &
     &      any(lobc(:,isfsur,ng))) THEN
          IF (master) WRITE (stdout,30) trim(vname(1,idsbry(isfsur))),  &
     &                                  trim(ncname)
          exit_flag=3
          RETURN
        END IF
        IF (.not.got_var(idsbry(isubar)).and.                           &
     &      any(lobc(:,isubar,ng))) THEN
          IF (master) WRITE (stdout,30) trim(vname(1,idsbry(isubar))),  &
     &                                  trim(ncname)
          exit_flag=3
          RETURN
        END IF
        IF (.not.got_var(idsbry(isvbar)).and.                           &
     &      any(lobc(:,isvbar,ng))) THEN
          IF (master) WRITE (stdout,30) trim(vname(1,idsbry(isvbar))),  &
     &                                  trim(ncname)
          exit_flag=3
          RETURN
        END IF
#  endif
#  ifdef ADJUST_WSTRESS
        IF (.not.got_var(idusms)) THEN
          IF (master) WRITE (stdout,30) trim(vname(1,idusms)),          &
     &                                  trim(ncname)
          exit_flag=3
          RETURN
        END IF
        IF (.not.got_var(idvsms)) THEN
          IF (master) WRITE (stdout,30) trim(vname(1,idvsms)),          &
     &                                  trim(ncname)
          exit_flag=3
          RETURN
        END IF
#  endif
#  ifdef SOLVE3D
        IF (.not.got_var(iduvel)) THEN
          IF (master) WRITE (stdout,30) trim(vname(1,iduvel)),          &
     &                                  trim(ncname)
          exit_flag=3
          RETURN
        END IF
        IF (.not.got_var(idvvel)) THEN
          IF (master) WRITE (stdout,30) trim(vname(1,idvvel)),          &
     &                                  trim(ncname)
          exit_flag=3
          RETURN
        END IF
#   ifdef ADJUST_BOUNDARY
        IF (.not.got_var(idsbry(isuvel)).and.                           &
     &      any(lobc(:,isuvel,ng))) THEN
          IF (master) WRITE (stdout,30) trim(vname(1,idsbry(isuvel))),  &
     &                                  trim(ncname)
          exit_flag=3
          RETURN
        END IF
        IF (.not.got_var(idsbry(isvvel)).and.                           &
     &      any(lobc(:,isvvel,ng))) THEN
          IF (master) WRITE (stdout,30) trim(vname(1,idsbry(isvvel))),  &
     &                                  trim(ncname)
          exit_flag=3
          RETURN
        END IF
#   endif
#   ifdef ADJUST_BOUNDARY
        IF (.not.got_var(idsbry(isuvel)).and.                           &
     &      any(lobc(:,isuvel,ng))) THEN
          IF (master) WRITE (stdout,30) trim(vname(1,idsbry(isuvel))),  &
     &                                  trim(ncname)
          exit_flag=3
          RETURN
        END IF
        IF (.not.got_var(idsbry(isvvel)).and.                           &
     &      any(lobc(:,isvvel,ng))) THEN
          IF (master) WRITE (stdout,30) trim(vname(1,idsbry(isvvel))),  &
     &                                  trim(ncname)
          exit_flag=3
          RETURN
        END IF
#   endif
#  endif
#  ifdef SOLVE3D
        DO itrc=1,nt(ng)
          IF (.not.got_var(idtvar(itrc))) THEN
            IF (master) WRITE (stdout,30) trim(vname(1,idtvar(itrc))),  &
     &                                    trim(ncname)
            exit_flag=3
            RETURN
          END IF
#   ifdef ADJUST_BOUNDARY
          IF (.not.got_var(idsbry(istvar(itrc))).and.                   &
     &        any(lobc(:,istvar(itrc),ng))) THEN
            IF (master) WRITE (stdout,30)                               &
     &                        trim(vname(1,idsbry(istvar(itrc)))),      &
     &                        trim(ncname)
            exit_flag=3
            RETURN
          END IF
#   endif
#   ifdef ADJUST_STFLUX
          IF (.not.got_var(idtsur(itrc)).and.lstflux(itrc,ng)) THEN
            IF (master) WRITE (stdout,30) trim(vname(1,idtsur(itrc))),  &
     &                                    trim(ncname)
            exit_flag=3
            RETURN
          END IF
#   endif
        END DO
#  endif
!
!  Set unlimited time record dimension to the appropriate value.
!
        irp(ng)%Rindex=rec_size
        fcount=irp(ng)%Fcount
        irp(ng)%Nrec(fcount)=rec_size
      END IF
!
  10  FORMAT (/,' RP_DEF_INI_PIO - unable to open initial NetCDF',      &
     &        ' file: ',a)
  20  FORMAT (/,' RP_DEF_INI_PIO - unable to put in define mode',       &
     &        ' initial NetCDF file: ',a)
  30  FORMAT (/,' RP_DEF_INI_PIO - unable to find variable: ',a,2x,     &
     &        ' in file: ',a)
!
      RETURN
      END SUBROUTINE rp_def_ini_pio
# endif
#endif
      END MODULE rp_def_ini_mod