rp_wrt_ini.F

Go to the documentation of this file.

#include "cppdefs.h"
      MODULE rp_wrt_ini_mod
#if (defined TANGENT || 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 writes initial conditions into representer model file   !
!  using either the standard NetCDF library or the Parallel-IO (PIO)   !
!  library.                                                            !
!                                                                      !
!  On Input:                                                           !
!                                                                      !
!     ng         Nested grid number.                                   !
!     Tindex     State variables time index to write.                  !
!     OutRec     NetCDF file unlimited dimension record to write.      !
!                                                                      !
!  Notice that only momentum is affected by the full time-averaged     !
!  masks.  If applicable, these mask contains information about        !
!  river runoff and time-dependent wetting and drying variations.      !
!                                                                      !
!=======================================================================
!
      USE mod_param
      USE mod_parallel
# ifdef ADJUST_BOUNDARY
      USE mod_boundary
# endif
# if defined ADJUST_WSTRESS || defined ADJUST_STFLUX
      USE mod_forces
# endif
      USE mod_fourdvar
      USE mod_grid
      USE mod_iounits
      USE mod_mixing
      USE mod_ncparam
      USE mod_ocean
      USE mod_scalars
# if defined SEDIMENT_NOT_YET || defined BBL_MODEL_NOT_YET
      USE mod_sediment
# endif
      USE mod_stepping
!
      USE nf_fwrite2d_mod,     ONLY : nf_fwrite2d
# ifdef ADJUST_BOUNDARY
      USE nf_fwrite2d_bry_mod, ONLY : nf_fwrite2d_bry
# endif
# ifdef SOLVE3D
      USE nf_fwrite3d_mod,     ONLY : nf_fwrite3d
#  ifdef ADJUST_BOUNDARY
      USE nf_fwrite3d_bry_mod, ONLY : nf_fwrite3d_bry
#  endif
# endif
      USE strings_mod,         ONLY : founderror
!
      implicit none
!
      PUBLIC  :: rp_wrt_ini
      PRIVATE :: rp_wrt_ini_nf90
# if defined PIO_LIB && defined DISTRIBUTE
      PRIVATE :: rp_wrt_ini_pio
# endif
!
      CONTAINS
!
!***********************************************************************
      SUBROUTINE rp_wrt_ini (ng, tile, Tindex, OutRec)
!***********************************************************************
!
!  Imported variable declarations.
!
      integer, intent(in) :: ng, tile, tindex, outrec
!
!  Local variable declarations.
!
# ifdef ADJUST_BOUNDARY
      integer :: lbij, ubij
# endif
      integer :: lbi, ubi, lbj, ubj
!
      character (len=*), parameter :: myfile =                          &
     &  __FILE__
!
!-----------------------------------------------------------------------
!  Write out history fields according to IO type.
!-----------------------------------------------------------------------
!
# ifdef ADJUST_BOUNDARY
      lbij=bounds(ng)%LBij
      ubij=bounds(ng)%UBij
# endif
      lbi=bounds(ng)%LBi(tile)
      ubi=bounds(ng)%UBi(tile)
      lbj=bounds(ng)%LBj(tile)
      ubj=bounds(ng)%UBj(tile)
!
      SELECT CASE (irp(ng)%IOtype)
        CASE (io_nf90)
          CALL rp_wrt_ini_nf90 (ng, tile, tindex, outrec,               &
# ifdef ADJUST_BOUNDARY
     &                          lbij, ubij,                             &
# endif
     &                          lbi, ubi, lbj, ubj)
 
# if defined PIO_LIB && defined DISTRIBUTE
        CASE (io_pio)
          CALL rp_wrt_ini_pio (ng, tile, tindex, outrec,                &
#  ifdef ADJUST_BOUNDARY
     &                         lbij, ubij,                              &
#  endif
     &                         lbi, ubi, lbj, ubj)
# 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_WRT_INI - Illegal output file type, io_type = ',i0,  &
     &        /,14x,'Check KeyWord ''OUT_LIB'' in ''roms.in''.')
!
      RETURN
      END SUBROUTINE rp_wrt_ini
!
!***********************************************************************
      SUBROUTINE rp_wrt_ini_nf90 (ng, tile, Tindex, OutRec,             &
# ifdef ADJUST_BOUNDARY
     &                            LBij, UBij,                           &
# endif
     &                            LBi, UBi, LBj, UBj)
!***********************************************************************
!
      USE mod_netcdf
!
!  Imported variable declarations.
!
      integer, intent(in) :: ng, tile, tindex, outrec
# ifdef ADJUST_BOUNDARY
      integer, intent(in) :: lbij, ubij
# endif
      integer, intent(in) :: lbi, ubi, lbj, ubj
!
!  Local variable declarations.
!
      integer :: gfactor, gtype, i, itrc, status, varid
!
      real(dp) :: my_time, scale
!
      character (len=15) :: string
 
      character (len=*), parameter :: myfile =                          &
     &  __FILE__//", rp_wrt_ini_nf90"
!
      sourcefile=myfile
!
!-----------------------------------------------------------------------
!  Write out tangent linear initial conditions.
!-----------------------------------------------------------------------
!
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Report.
!
      IF (master) THEN
        IF (outrec.eq.2) THEN
          string='initial  fields'
        END IF
# ifdef SOLVE3D
        WRITE (stdout,10) string, outer, inner, tindex, tindex, outrec
# else
        WRITE (stdout,10) string, outer, inner, tindex, outrec
# endif
      END IF
!
!  Set grid type factor to write full (gfactor=1) fields or water
!  points (gfactor=-1) fields only.
!
# if defined WRITE_WATER && defined MASKING
      gfactor=-1
# else
      gfactor=1
# endif
!
!  Write out model time (s).  Use the "tdays" variable here because of
!  the management of the "time" variable due to nesting.
!
      my_time=tdays(ng)*day2sec
 
      CALL netcdf_put_fvar (ng, irpm, irp(ng)%name,                     &
     &                      trim(vname(1,idtime)), my_time,             &
     &                      (/outrec/), (/1/),                          &
     &                      ncid = irp(ng)%ncid,                        &
     &                      varid = irp(ng)%Vid(idtime))
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Write out free-surface (m)
!
      scale=1.0_dp
      gtype=gfactor*r2dvar
      status=nf_fwrite2d(ng, irpm, irp(ng)%ncid, idfsur,                &
     &                   irp(ng)%Vid(idfsur),                           &
     &                   outrec, gtype,                                 &
     &                   lbi, ubi, lbj, ubj, scale,                     &
# ifdef MASKING
     &                   grid(ng) % rmask,                              &
# endif
# ifdef WET_DRY
     &                   ocean(ng) % tl_zeta(:,:,tindex),               &
     &                   setfillval = .false.)
# else
     &                   ocean(ng) % tl_zeta(:,:,tindex))
# endif
      IF (founderror(status, nf90_noerr, __line__, myfile)) THEN
        IF (master) THEN
          WRITE (stdout,20) trim(vname(1,idfsur)), outrec
        END IF
        exit_flag=3
        ioerror=status
        RETURN
      END IF
 
# ifdef ADJUST_BOUNDARY
!
!  Write out free-surface open boundaries.
!
      IF (any(lobc(:,isfsur,ng))) THEN
        scale=1.0_dp
        status=nf_fwrite2d_bry(ng, irpm, irp(ng)%name, irp(ng)%ncid,    &
     &                         vname(1,idsbry(isfsur)),                 &
     &                         irp(ng)%Vid(idsbry(isfsur)),             &
     &                         outrec, r2dvar,                          &
     &                         lbij, ubij, nbrec(ng), scale,            &
     &                         boundary(ng) % tl_zeta_obc(lbij:,:,:,    &
     &                                                    tindex))
        IF (founderror(status, nf90_noerr, __line__, myfile)) THEN
          IF (master) THEN
            WRITE (stdout,20) trim(vname(1,idsbry(isfsur))), outrec
          END IF
          exit_flag=3
          ioerror=status
          RETURN
        END IF
      END IF
# endif
!
!  Write out 2D U-momentum component (m/s).
!
      scale=1.0_dp
      gtype=gfactor*u2dvar
      status=nf_fwrite2d(ng, irpm, irp(ng)%ncid, idubar,                &
     &                   irp(ng)%Vid(idubar),                           &
     &                   outrec, gtype,                                 &
     &                   lbi, ubi, lbj, ubj, scale,                     &
# ifdef MASKING
     &                   grid(ng) % umask_full,                         &
# endif
     &                   ocean(ng) % tl_ubar(:,:,tindex))
      IF (founderror(status, nf90_noerr, __line__, myfile)) THEN
        IF (master) THEN
          WRITE (stdout,20) trim(vname(1,idubar)), outrec
        END IF
        exit_flag=3
        ioerror=status
        RETURN
      END IF
 
# ifdef ADJUST_BOUNDARY
!
!  Write out 2D U-momentum component open boundaries.
!
      IF (any(lobc(:,isubar,ng))) THEN
        scale=1.0_dp
        status=nf_fwrite2d_bry(ng, irpm, irp(ng)%name, irp(ng)%ncid,    &
     &                         vname(1,idsbry(isubar)),                 &
     &                         irp(ng)%Vid(idsbry(isubar)),             &
     &                         outrec, u2dvar,                          &
     &                         lbij, ubij, nbrec(ng), scale,            &
     &                         boundary(ng) % tl_ubar_obc(lbij:,:,:,    &
     &                                                    tindex))
        IF (founderror(status, nf90_noerr, __line__, myfile)) THEN
          IF (master) THEN
            WRITE (stdout,20) trim(vname(1,idsbry(isubar))), outrec
          END IF
          exit_flag=3
          ioerror=status
          RETURN
        END IF
      END IF
# endif
!
!  Write out 2D V-momentum component (m/s).
!
      scale=1.0_dp
      gtype=gfactor*v2dvar
      status=nf_fwrite2d(ng, irpm, irp(ng)%ncid, idvbar,                &
     &                   irp(ng)%Vid(idvbar),                           &
     &                   outrec, gtype,                                 &
     &                   lbi, ubi, lbj, ubj, scale,                     &
# ifdef MASKING
     &                   grid(ng) % vmask_full,                         &
# endif
     &                   ocean(ng) % tl_vbar(:,:,tindex))
      IF (founderror(status, nf90_noerr, __line__, myfile)) THEN
        IF (master) THEN
          WRITE (stdout,20) trim(vname(1,idvbar)), outrec
        END IF
        exit_flag=3
        ioerror=status
        RETURN
      END IF
 
# ifdef ADJUST_BOUNDARY
!
!  Write out 2D V-momentum component open boundaries.
!
      IF (any(lobc(:,isvbar,ng))) THEN
        scale=1.0_dp
        status=nf_fwrite2d_bry(ng, irpm, irp(ng)%name, irp(ng)%ncid,    &
     &                         vname(1,idsbry(isvbar)),                 &
     &                         irp(ng)%Vid(idsbry(isvbar)),             &
     &                         outrec, v2dvar,                          &
     &                         lbij, ubij, nbrec(ng), scale,            &
     &                         boundary(ng) % tl_vbar_obc(lbij:,:,:,    &
     &                                                    tindex))
        IF (founderror(status, nf90_noerr, __line__, myfile)) THEN
          IF (master) THEN
            WRITE (stdout,20) trim(vname(1,idsbry(isvbar))), outrec
          END IF
          exit_flag=3
          ioerror=status
          RETURN
        END IF
      END IF
# endif
 
# ifdef ADJUST_WSTRESS
!
!  Write out surface U-momentum stress.  Notice that the stress has its
!  own fixed time-dimension (of size Nfrec) to allow 4DVAR adjustments
!  at other times in addition to initialization time.
!
!!    scale=rho0
      scale=1.0_dp
      gtype=gfactor*u3dvar
      status=nf_fwrite3d(ng, irpm, irp(ng)%ncid, idusms,                &
     &                   irp(ng)%Vid(idusms),                           &
     &                   outrec, gtype,                                 &
     &                   lbi, ubi, lbj, ubj, 1, nfrec(ng), scale,       &
#  ifdef MASKING
     &                   grid(ng) % umask,                              &
#  endif
     &                   forces(ng) % tl_ustr(:,:,:,tindex))
      IF (founderror(status, nf90_noerr, __line__, myfile)) THEN
        IF (master) THEN
          WRITE (stdout,20) trim(vname(1,idusms)), outrec
        END IF
        exit_flag=3
        ioerror=status
        RETURN
      END IF
!
!  Write out surface V-momentum stress.
!
!!    scale=rho0
      scale=1.0_dp
      gtype=gfactor*v3dvar
      status=nf_fwrite3d(ng, irpm, irp(ng)%ncid, idvsms,                &
     &                   irp(ng)%Vid(idvsms),                           &
     &                   outrec, gtype,                                 &
     &                   lbi, ubi, lbj, ubj, 1, nfrec(ng), scale,       &
#  ifdef MASKING
     &                   grid(ng) % vmask,                              &
#  endif
     &                   forces(ng) % tl_vstr(:,:,:,tindex))
      IF (founderror(status, nf90_noerr, __line__, myfile)) THEN
        IF (master) THEN
          WRITE (stdout,20) trim(vname(1,idvsms)), outrec
        END IF
        exit_flag=3
        ioerror=status
        RETURN
      END IF
# endif
 
# ifdef SOLVE3D
!
!  Write out 3D U-momentum component (m/s).
!
      scale=1.0_dp
      gtype=gfactor*u3dvar
      status=nf_fwrite3d(ng, irpm, irp(ng)%ncid, iduvel,                &
     &                   irp(ng)%Vid(iduvel),                           &
     &                   outrec, gtype,                                 &
     &                   lbi, ubi, lbj, ubj, 1, n(ng), scale,           &
#  ifdef MASKING
     &                   grid(ng) % umask_full,                         &
#  endif
     &                   ocean(ng) % tl_u(:,:,:,tindex))
      IF (founderror(status, nf90_noerr, __line__, myfile)) THEN
        IF (master) THEN
          WRITE (stdout,20) trim(vname(1,iduvel)), outrec
        END IF
        exit_flag=3
        ioerror=status
        RETURN
      END IF
 
#  ifdef ADJUST_BOUNDARY
!
!  Write out 3D U-momentum component open boundaries.
!
      IF (any(lobc(:,isuvel,ng))) THEN
        scale=1.0_dp
        status=nf_fwrite3d_bry(ng, irpm, irp(ng)%name, irp(ng)%ncid,    &
     &                         vname(1,idsbry(isuvel)),                 &
     &                         irp(ng)%Vid(idsbry(isuvel)),             &
     &                         outrec, u3dvar,                          &
     &                         lbij, ubij, 1, n(ng), nbrec(ng), scale,  &
     &                         boundary(ng) % tl_u_obc(lbij:,:,:,:,     &
     &                                                 tindex))
        IF (founderror(status, nf90_noerr, __line__, myfile)) THEN
          IF (master) THEN
            WRITE (stdout,20) trim(vname(1,idsbry(isuvel))), outrec
          END IF
          exit_flag=3
          ioerror=status
          RETURN
        END IF
      END IF
#  endif
!
!  Write out 3D V-momentum component (m/s).
!
      scale=1.0_dp
      gtype=gfactor*v3dvar
      status=nf_fwrite3d(ng, irpm, irp(ng)%ncid, idvvel,                &
     &                   irp(ng)%Vid(idvvel),                           &
     &                   outrec, gtype,                                 &
     &                   lbi, ubi, lbj, ubj, 1, n(ng), scale,           &
#  ifdef MASKING
     &                   grid(ng) % vmask_full,                         &
#  endif
     &                   ocean(ng) % tl_v(:,:,:,tindex))
      IF (founderror(status, nf90_noerr, __line__, myfile)) THEN
        IF (master) THEN
          WRITE (stdout,20) trim(vname(1,idvvel)), outrec
        END IF
        exit_flag=3
        ioerror=status
        RETURN
      END IF
 
#  ifdef ADJUST_BOUNDARY
!
!  Write out 3D V-momentum component open boundaries.
!
      IF (any(lobc(:,isvvel,ng))) THEN
        scale=1.0_dp
        status=nf_fwrite3d_bry(ng, irpm, irp(ng)%name, irp(ng)%ncid,    &
     &                         vname(1,idsbry(isvvel)),                 &
     &                         irp(ng)%Vid(idsbry(isvvel)),             &
     &                         outrec, v3dvar,                          &
     &                         lbij, ubij, 1, n(ng), nbrec(ng), scale,  &
     &                         boundary(ng) % tl_v_obc(lbij:,:,:,:,     &
     &                                                 tindex))
        IF (founderror(status, nf90_noerr, __line__, myfile)) THEN
          IF (master) THEN
            WRITE (stdout,20) trim(vname(1,idsbry(isvvel))), outrec
          END IF
          exit_flag=3
          ioerror=status
          RETURN
        END IF
      END IF
#  endif
!
!  Write out tracer type variables.
!
      DO itrc=1,nt(ng)
        scale=1.0_dp
        gtype=gfactor*r3dvar
        status=nf_fwrite3d(ng, irpm, irp(ng)%ncid, idtvar(itrc),        &
     &                     irp(ng)%Tid(itrc),                           &
     &                     outrec, gtype,                               &
     &                     lbi, ubi, lbj, ubj, 1, n(ng), scale,         &
#  ifdef MASKING
     &                     grid(ng) % rmask,                            &
#  endif
     &                     ocean(ng) % tl_t(:,:,:,tindex,itrc))
        IF (founderror(status, nf90_noerr, __line__, myfile)) THEN
          IF (master) THEN
            WRITE (stdout,20) trim(vname(1,idtvar(itrc))), outrec
          END IF
          exit_flag=3
          ioerror=status
          RETURN
        END IF
      END DO
 
#  ifdef ADJUST_BOUNDARY
!
!  Write out tracers open boundaries.
!
      DO itrc=1,nt(ng)
        IF (any(lobc(:,istvar(itrc),ng))) THEN
          scale=1.0_dp
          status=nf_fwrite3d_bry(ng, irpm, irp(ng)%name, irp(ng)%ncid,  &
     &                           vname(1,idsbry(istvar(itrc))),         &
     &                           irp(ng)%Vid(idsbry(istvar(itrc))),     &
     &                           outrec, r3dvar,                        &
     &                           lbij, ubij, 1, n(ng), nbrec(ng),       &
     &                           scale,                                 &
     &                           boundary(ng) % tl_t_obc(lbij:,:,:,:,   &
     &                                                   tindex,itrc))
          IF (founderror(status, nf90_noerr, __line__, myfile)) THEN
            IF (master) THEN
              WRITE (stdout,20) trim(vname(1,idsbry(istvar(itrc)))),    &
     &                          outrec
            END IF
            exit_flag=3
            ioerror=status
            RETURN
          END IF
        END IF
      END DO
#  endif
 
#  ifdef ADJUST_STFLUX
!
!  Write out surface net tracers fluxes. Notice that fluxes have their
!  own fixed time-dimension (of size Nfrec) to allow 4DVAR adjustments
!  at other times in addition to initialization time.
!
      DO itrc=1,nt(ng)
        IF (lstflux(itrc,ng)) THEN
          scale=1.0_dp
          gtype=gfactor*r3dvar
          status=nf_fwrite3d(ng, irpm, irp(ng)%ncid, idtsur(itrc),      &
     &                       irp(ng)%Vid(idtsur(itrc)),                 &
     &                       outrec, gtype,                             &
     &                       lbi, ubi, lbj, ubj, 1, nfrec(ng), scale,   &
#   ifdef MASKING
     &                       grid(ng) % rmask,                          &
#   endif
     &                       forces(ng) % tl_tflux(:,:,:,tindex,itrc))
          IF (founderror(status, nf90_noerr, __line__, myfile)) THEN
            IF (master) THEN
              WRITE (stdout,20) trim(vname(1,idtsur(itrc))),            &
     &                          outrec
            END IF
            exit_flag=3
            ioerror=status
            RETURN
          END IF
        END IF
      END DO
#  endif
# endif
!
!-----------------------------------------------------------------------
!  Synchronize tangent linear initial NetCDF file to disk to allow other
!  processes to access data immediately after it is written.
!-----------------------------------------------------------------------
!
      CALL netcdf_sync (ng, irpm, irp(ng)%name, irp(ng)%ncid)
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
  10  FORMAT (2x,'RP_WRT_INI_NF90  - writing ',a,                       &
     &        ' (Outer=',i2.2,', Inner=',i3.3,', Index=',i0,            &
# ifdef SOLVE3D
     &        ',',i0,', Rec=',i0,')')
# else
     &        ', Rec=',i0,')')
# endif
  20  FORMAT (/,' RP_WRT_INI_NF90 - error while writing variable: ',a,  &
     &        /,19x,'into representers initial file for time record:',  &
     &        1x,i4)
!
      RETURN
      END SUBROUTINE rp_wrt_ini_nf90
 
# if defined PIO_LIB && defined DISTRIBUTE
!
!***********************************************************************
      SUBROUTINE rp_wrt_ini_pio (ng, tile, Tindex, OutRec,              &
#  ifdef ADJUST_BOUNDARY
     &                           LBij, UBij,                            &
#  endif
     &                           LBi, UBi, LBj, UBj)
!***********************************************************************
!
      USE mod_pio_netcdf
!
!  Imported variable declarations.
!
      integer, intent(in) :: ng, tile, tindex, outrec
#  ifdef ADJUST_BOUNDARY
      integer, intent(in) :: lbij, ubij
#  endif
      integer, intent(in) :: lbi, ubi, lbj, ubj
!
!  Local variable declarations.
!
      integer :: i, ifield, itrc, status, varid
!
      real(dp) :: my_time, scale
!
      character (len=15) :: string
 
      character (len=*), parameter :: myfile =                          &
     &  __FILE__//", rp_wrt_ini_pio"
!
      TYPE (io_desc_t), pointer :: iodesc
!
      sourcefile=myfile
!
!-----------------------------------------------------------------------
!  Write out tangent linear initial conditions.
!-----------------------------------------------------------------------
!
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Report.
!
      IF (master) THEN
        IF (outrec.eq.2) THEN
          string='initial  fields'
        END IF
#  ifdef SOLVE3D
        WRITE (stdout,10) string, outer, inner, tindex, tindex, outrec
#  else
        WRITE (stdout,10) string, outer, inner, tindex, outrec
#  endif
      END IF
!
!  Write out model time (s).  Use the "tdays" variable here because of
!  the management of the "time" variable due to nesting.
!
      my_time=tdays(ng)*day2sec
 
      CALL pio_netcdf_put_fvar (ng, irpm, irp(ng)%name,                 &
     &                          trim(vname(1,idtime)), my_time,         &
     &                          (/outrec/), (/1/),                      &
     &                          piofile = irp(ng)%pioFile,              &
     &                          piovar = irp(ng)%pioVar(idtime)%vd)
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Write out free-surface (m)
!
      scale=1.0_dp
      IF (irp(ng)%pioVar(idfsur)%dkind.eq.pio_double) THEN
        iodesc => iodesc_dp_r2dvar(ng)
      ELSE
        iodesc => iodesc_sp_r2dvar(ng)
      END IF
!
      status=nf_fwrite2d(ng, irpm, irp(ng)%pioFile, idfsur,             &
     &                   irp(ng)%pioVar(idfsur),                        &
     &                   outrec, iodesc,                                &
     &                   lbi, ubi, lbj, ubj, scale,                     &
#  ifdef MASKING
     &                   grid(ng) % rmask,                              &
#  endif
#  ifdef WET_DRY
     &                   ocean(ng) % tl_zeta(:,:,tindex),               &
     &                   setfillval = .false.)
#  else
     &                   ocean(ng) % tl_zeta(:,:,tindex))
#  endif
      IF (founderror(status, pio_noerr, __line__, myfile)) THEN
        IF (master) THEN
          WRITE (stdout,20) trim(vname(1,idfsur)), outrec
        END IF
        exit_flag=3
        ioerror=status
        RETURN
      END IF
 
#  ifdef ADJUST_BOUNDARY
!
!  Write out free-surface open boundaries.
!
      IF (any(lobc(:,isfsur,ng))) THEN
        scale=1.0_dp
        IF (irp(ng)%pioVar(idsbry(isfsur))%dkind.eq.pio_double) THEN
          iodesc => iodesc_dp_r2dobc(ng)
        ELSE
          iodesc => iodesc_sp_r2dobc(ng)
        END IF
!
        status=nf_fwrite2d_bry(ng, irpm, irp(ng)%name,                  &
     &                         irp(ng)%pioFile,                         &
     &                         vname(1,idsbry(isfsur)),                 &
     &                         irp(ng)%pioVar(idsbry(isfsur)),          &
     &                         outrec, iodesc,                          &
     &                         lbij, ubij, nbrec(ng), scale,            &
     &                         boundary(ng) % tl_zeta_obc(lbij:,:,:,    &
     &                                                    tindex))
        IF (founderror(status, pio_noerr, __line__, myfile)) THEN
          IF (master) THEN
            WRITE (stdout,20) trim(vname(1,idsbry(isfsur))), outrec
          END IF
          exit_flag=3
          ioerror=status
          RETURN
        END IF
      END IF
#  endif
!
!  Write out 2D U-momentum component (m/s).
!
      scale=1.0_dp
      IF (irp(ng)%pioVar(idubar)%dkind.eq.pio_double) THEN
        iodesc => iodesc_dp_u2dvar(ng)
      ELSE
        iodesc => iodesc_sp_u2dvar(ng)
      END IF
!
      status=nf_fwrite2d(ng, irpm, irp(ng)%pioFile, idubar,             &
     &                   irp(ng)%pioVar(idubar),                        &
     &                   outrec, iodesc,                                &
     &                   lbi, ubi, lbj, ubj, scale,                     &
#  ifdef MASKING
     &                   grid(ng) % umask_full,                         &
#  endif
     &                   ocean(ng) % tl_ubar(:,:,tindex))
      IF (founderror(status, pio_noerr, __line__, myfile)) THEN
        IF (master) THEN
          WRITE (stdout,20) trim(vname(1,idubar)), outrec
        END IF
        exit_flag=3
        ioerror=status
        RETURN
      END IF
 
#  ifdef ADJUST_BOUNDARY
!
!  Write out 2D U-momentum component open boundaries.
!
      IF (any(lobc(:,isubar,ng))) THEN
        scale=1.0_dp
        IF (irp(ng)%pioVar(idsbry(isubar))%dkind.eq.pio_double) THEN
          iodesc => iodesc_dp_u2dobc(ng)
        ELSE
          iodesc => iodesc_sp_u2dobc(ng)
        END IF
!
        status=nf_fwrite2d_bry(ng, irpm, irp(ng)%name,                  &
     &                         irp(ng)%pioFile,                         &
     &                         vname(1,idsbry(isubar)),                 &
     &                         irp(ng)%pioVar(idsbry(isubar)),          &
     &                         outrec, iodesc,                          &
     &                         lbij, ubij, nbrec(ng), scale,            &
     &                         boundary(ng) % tl_ubar_obc(lbij:,:,:,    &
     &                                                    tindex))
        IF (founderror(status, pio_noerr, __line__, myfile)) THEN
          IF (master) THEN
            WRITE (stdout,20) trim(vname(1,idsbry(isubar))), outrec
          END IF
          exit_flag=3
          ioerror=status
          RETURN
        END IF
      END IF
#  endif
!
!  Write out 2D V-momentum component (m/s).
!
      scale=1.0_dp
      IF (irp(ng)%pioVar(idvbar)%dkind.eq.pio_double) THEN
        iodesc => iodesc_dp_v2dvar(ng)
      ELSE
        iodesc => iodesc_sp_v2dvar(ng)
      END IF
!
      status=nf_fwrite2d(ng, irpm, irp(ng)%pioFile, idvbar,             &
     &                   irp(ng)%pioVar(idvbar),                        &
     &                   outrec, iodesc,                                &
     &                   lbi, ubi, lbj, ubj, scale,                     &
#  ifdef MASKING
     &                   grid(ng) % vmask_full,                         &
#  endif
     &                   ocean(ng) % tl_vbar(:,:,tindex))
      IF (founderror(status, pio_noerr, __line__, myfile)) THEN
        IF (master) THEN
          WRITE (stdout,20) trim(vname(1,idvbar)), outrec
        END IF
        exit_flag=3
        ioerror=status
        RETURN
      END IF
 
#  ifdef ADJUST_BOUNDARY
!
!  Write out 2D V-momentum component open boundaries.
!
      IF (any(lobc(:,isvbar,ng))) THEN
        scale=1.0_dp
        IF (irp(ng)%pioVar(idsbry(isvbar))%dkind.eq.pio_double) THEN
          iodesc => iodesc_dp_v2dobc(ng)
        ELSE
          iodesc => iodesc_sp_v2dobc(ng)
        END IF
!
        status=nf_fwrite2d_bry(ng, irpm, irp(ng)%name,                  &
     &                         irp(ng)%pioFile,                         &
     &                         vname(1,idsbry(isvbar)),                 &
     &                         irp(ng)%pioVar(idsbry(isvbar)),          &
     &                         outrec, iodesc,                          &
     &                         lbij, ubij, nbrec(ng), scale,            &
     &                         boundary(ng) % tl_vbar_obc(lbij:,:,:,    &
     &                                                    tindex))
        IF (founderror(status, pio_noerr, __line__, myfile)) THEN
          IF (master) THEN
            WRITE (stdout,20) trim(vname(1,idsbry(isvbar))), outrec
          END IF
          exit_flag=3
          ioerror=status
          RETURN
        END IF
      END IF
#  endif
 
#  ifdef ADJUST_WSTRESS
!
!  Write out surface U-momentum stress.  Notice that the stress has its
!  own fixed time-dimension (of size Nfrec) to allow 4DVAR adjustments
!  at other times in addition to initialization time.
!
!!    scale=rho0
      scale=1.0_dp
      IF (irp(ng)%pioVar(idusms)%dkind.eq.pio_double) THEN
        iodesc => iodesc_dp_u2dfrc(ng)
      ELSE
        iodesc => iodesc_sp_u2dfrc(ng)
      END IF
!
      status=nf_fwrite3d(ng, irpm, irp(ng)%pioFile, idusms,             &
     &                   irp(ng)%pioVar(idusms),                        &
     &                   outrec, iodesc,                                &
     &                   lbi, ubi, lbj, ubj, 1, nfrec(ng), scale,       &
#   ifdef MASKING
     &                   grid(ng) % umask,                              &
#   endif
     &                   forces(ng) % tl_ustr(:,:,:,tindex))
      IF (founderror(status, pio_noerr, __line__, myfile)) THEN
        IF (master) THEN
          WRITE (stdout,20) trim(vname(1,idusms)), outrec
        END IF
        exit_flag=3
        ioerror=status
        RETURN
      END IF
!
!  Write out surface V-momentum stress.
!
!!    scale=rho0
      scale=1.0_dp
      IF (irp(ng)%pioVar(idvsms)%dkind.eq.pio_double) THEN
        iodesc => iodesc_dp_v2dfrc(ng)
      ELSE
        iodesc => iodesc_sp_v2dfrc(ng)
      END IF
!
      status=nf_fwrite3d(ng, irpm, irp(ng)%pioFile, idvsms,             &
     &                   irp(ng)%pioVar(idvsms),                        &
     &                   outrec, iodesc,                                &
     &                   lbi, ubi, lbj, ubj, 1, nfrec(ng), scale,       &
#   ifdef MASKING
     &                   grid(ng) % vmask,                              &
#   endif
     &                   forces(ng) % tl_vstr(:,:,:,tindex))
      IF (founderror(status, pio_noerr, __line__, myfile)) THEN
        IF (master) THEN
          WRITE (stdout,20) trim(vname(1,idvsms)), outrec
        END IF
        exit_flag=3
        ioerror=status
        RETURN
      END IF
#  endif
 
#  ifdef SOLVE3D
!
!  Write out 3D U-momentum component (m/s).
!
      scale=1.0_dp
      IF (irp(ng)%pioVar(iduvel)%dkind.eq.pio_double) THEN
        iodesc => iodesc_dp_u3dvar(ng)
      ELSE
        iodesc => iodesc_sp_u3dvar(ng)
      END IF
!
      status=nf_fwrite3d(ng, irpm, irp(ng)%pioFile, iduvel,             &
     &                   irp(ng)%pioVar(iduvel),                        &
     &                   outrec, iodesc,                                &
     &                   lbi, ubi, lbj, ubj, 1, n(ng), scale,           &
#   ifdef MASKING
     &                   grid(ng) % umask_full,                         &
#   endif
     &                   ocean(ng) % tl_u(:,:,:,tindex))
      IF (founderror(status, pio_noerr, __line__, myfile)) THEN
        IF (master) THEN
          WRITE (stdout,20) trim(vname(1,iduvel)), outrec
        END IF
        exit_flag=3
        ioerror=status
        RETURN
      END IF
 
#   ifdef ADJUST_BOUNDARY
!
!  Write out 3D U-momentum component open boundaries.
!
      IF (any(lobc(:,isuvel,ng))) THEN
        scale=1.0_dp
        IF (irp(ng)%pioVar(idsbry(isuvel))%dkind.eq.pio_double) THEN
          iodesc => iodesc_dp_u3dobc(ng)
        ELSE
          iodesc => iodesc_sp_u3dobc(ng)
        END IF
!
        status=nf_fwrite3d_bry(ng, irpm, irp(ng)%name,                  &
     &                         irp(ng)%pioFile,                         &
     &                         vname(1,idsbry(isuvel)),                 &
     &                         irp(ng)%pioVar(idsbry(isuvel)),          &
     &                         outrec, iodesc,                          &
     &                         lbij, ubij, 1, n(ng), nbrec(ng), scale,  &
     &                         boundary(ng) % tl_u_obc(lbij:,:,:,:,     &
     &                                                 tindex))
        IF (founderror(status, pio_noerr, __line__, myfile)) THEN
          IF (master) THEN
            WRITE (stdout,20) trim(vname(1,idsbry(isuvel))), outrec
          END IF
          exit_flag=3
          ioerror=status
          RETURN
        END IF
      END IF
#   endif
!
!  Write out 3D V-momentum component (m/s).
!
      scale=1.0_dp
      IF (irp(ng)%pioVar(idvvel)%dkind.eq.pio_double) THEN
        iodesc => iodesc_dp_v3dvar(ng)
      ELSE
        iodesc => iodesc_sp_v3dvar(ng)
      END IF
!
      status=nf_fwrite3d(ng, irpm, irp(ng)%pioFile, idvvel,             &
     &                   irp(ng)%pioVar(idvvel),                        &
     &                   outrec, iodesc,                                &
     &                   lbi, ubi, lbj, ubj, 1, n(ng), scale,           &
#   ifdef MASKING
     &                   grid(ng) % vmask_full,                         &
#   endif
     &                   ocean(ng) % tl_v(:,:,:,tindex))
      IF (founderror(status, pio_noerr, __line__, myfile)) THEN
        IF (master) THEN
          WRITE (stdout,20) trim(vname(1,idvvel)), outrec
        END IF
        exit_flag=3
        ioerror=status
        RETURN
      END IF
 
#   ifdef ADJUST_BOUNDARY
!
!  Write out 3D V-momentum component open boundaries.
!
      IF (any(lobc(:,isvvel,ng))) THEN
        scale=1.0_dp
        IF (irp(ng)%pioVar(idsbry(isvvel))%dkind.eq.pio_double) THEN
          iodesc => iodesc_dp_v3dobc(ng)
        ELSE
          iodesc => iodesc_sp_v3dobc(ng)
        END IF
!
        status=nf_fwrite3d_bry(ng, irpm, irp(ng)%name,                  &
     &                         irp(ng)%pioFile,                         &
     &                         vname(1,idsbry(isvvel)),                 &
     &                         irp(ng)%pioVar(idsbry(isvvel)),          &
     &                         outrec, iodesc,                          &
     &                         lbij, ubij, 1, n(ng), nbrec(ng), scale,  &
     &                         boundary(ng) % tl_v_obc(lbij:,:,:,:,     &
     &                                                 tindex))
        IF (founderror(status, pio_noerr, __line__, myfile)) THEN
          IF (master) THEN
            WRITE (stdout,20) trim(vname(1,idsbry(isvvel))), outrec
          END IF
          exit_flag=3
          ioerror=status
          RETURN
        END IF
      END IF
#   endif
!
!  Write out tracer type variables.
!
      DO itrc=1,nt(ng)
        scale=1.0_dp
        IF (irp(ng)%pioTrc(itrc)%dkind.eq.pio_double) THEN
          iodesc => iodesc_dp_r3dvar(ng)
        ELSE
          iodesc => iodesc_sp_r3dvar(ng)
        END IF
!
        status=nf_fwrite3d(ng, irpm, irp(ng)%pioFile, idtvar(itrc),     &
     &                     irp(ng)%pioTrc(itrc),                        &
     &                     outrec, iodesc,                              &
     &                     lbi, ubi, lbj, ubj, 1, n(ng), scale,         &
#   ifdef MASKING
     &                     grid(ng) % rmask,                            &
#   endif
     &                     ocean(ng) % tl_t(:,:,:,tindex,itrc))
        IF (founderror(status, pio_noerr, __line__, myfile)) THEN
          IF (master) THEN
            WRITE (stdout,20) trim(vname(1,idtvar(itrc))), outrec
          END IF
          exit_flag=3
          ioerror=status
          RETURN
        END IF
      END DO
 
#   ifdef ADJUST_BOUNDARY
!
!  Write out tracers open boundaries.
!
      DO itrc=1,nt(ng)
        IF (any(lobc(:,istvar(itrc),ng))) THEN
          scale=1.0_dp
          ifield=idsbry(istvar(itrc))
          IF (irp(ng)%pioVar(ifield)%dkind.eq.pio_double) THEN
            iodesc => iodesc_dp_r3dobc(ng)
          ELSE
            iodesc => iodesc_sp_r3dobc(ng)
          END IF
!
          status=nf_fwrite3d_bry(ng, irpm, irp(ng)%name,                &
     &                           irp(ng)%pioFile,                       &
     &                           vname(1,ifield),                       &
     &                           irp(ng)%pioVar(ifield),                &
     &                           outrec, iodesc,                        &
     &                           lbij, ubij, 1, n(ng), nbrec(ng),       &
     &                           scale,                                 &
     &                           boundary(ng) % tl_t_obc(lbij:,:,:,:,   &
     &                                                   tindex,itrc))
          IF (founderror(status, pio_noerr, __line__, myfile)) THEN
            IF (master) THEN
              WRITE (stdout,20) trim(vname(1,idsbry(istvar(itrc)))),    &
     &                          outrec
            END IF
            exit_flag=3
            ioerror=status
            RETURN
          END IF
        END IF
      END DO
#   endif
 
#   ifdef ADJUST_STFLUX
!
!  Write out surface net tracers fluxes. Notice that fluxes have their
!  own fixed time-dimension (of size Nfrec) to allow 4DVAR adjustments
!  at other times in addition to initialization time.
!
      DO itrc=1,nt(ng)
        IF (lstflux(itrc,ng)) THEN
          scale=1.0_dp
          IF (irp(ng)%pioVar(idtsur(itrc))%dkind.eq.pio_double) THEN
            iodesc => iodesc_dp_r2dfrc(ng)
          ELSE
            iodesc => iodesc_sp_r2dfrc(ng)
          END IF
!
          status=nf_fwrite3d(ng, irpm, irp(ng)%pioFile, idtsur(itrc),   &
     &                       irp(ng)%pioVar(idtsur(itrc)),              &
     &                       outrec, iodesc,                            &
     &                       lbi, ubi, lbj, ubj, 1, nfrec(ng), scale,   &
#    ifdef MASKING
     &                       grid(ng) % rmask,                          &
#    endif
     &                       forces(ng) % tl_tflux(:,:,:,tindex,itrc))
          IF (founderror(status, pio_noerr, __line__, myfile)) THEN
            IF (master) THEN
              WRITE (stdout,20) trim(vname(1,idtsur(itrc))),            &
     &                          outrec
            END IF
            exit_flag=3
            ioerror=status
            RETURN
          END IF
        END IF
      END DO
#   endif
#  endif
!
!-----------------------------------------------------------------------
!  Synchronize tangent linear initial NetCDF file to disk to allow other
!  processes to access data immediately after it is written.
!-----------------------------------------------------------------------
!
      CALL pio_netcdf_sync (ng, irpm, irp(ng)%name, irp(ng)%pioFile)
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
  10  FORMAT (2x,'RP_WRT_INI_PIO   - writing ',a,                       &
     &        ' (Outer=',i2.2,', Inner=',i3.3,', Index=',i0,            &
#  ifdef SOLVE3D
     &        ',',i0,', Rec=',i0,')')
#  else
     &        ', Rec=',i0,')')
#  endif
  20  FORMAT (/,' RP_WRT_INI_PIO - error while writing variable: ',a,   &
     &        /,18x,'into representers initial file for time record:',  &
     &        1x,i4)
!
      RETURN
      END SUBROUTINE rp_wrt_ini_pio
# endif
#endif
      END MODULE rp_wrt_ini_mod