jedi_roms.h

Go to the documentation of this file.

      MODULE roms_kernel_mod
!
!git $Id$
!================================================== Hernan G. Arango ===
!  Copyright (c) 2002-2025 The ROMS Group                              !
!    Licensed under a MIT/X style license                              !
!    See License_ROMS.md                                               !
!=======================================================================
!                                                                      !
!  ROMSJEDI Interface: Joint Effort for Data-assimilation Integration  !
!                                                                      !
!  The routines in this driver control the initialization, time-       !
!  stepping, and finalization of ROMSJEDI interface following          !
!  ESMF/NUOPC conventions:                                             !
!                                                                      !
!     ROMS_initializeP1         Phase 1 ROMSJEDI Initialization        !
!     ROMS_initializeP2         Phase 2 ROMSJEDI Initialization        !
!     ROMS_initializeP3         Phase 3 ROMSJEDI Initialization        !
!     ROMS_run                                                         !
!     ROMS_finalize                                                    !
!                                                                      !
!=======================================================================
!
      USE mod_param
      USE mod_parallel
      USE mod_arrays
      USE mod_fourdvar
      USE mod_iounits
      USE mod_ncparam
#ifdef NESTING
      USE mod_nesting
#endif
      USE mod_ocean
      USE mod_scalars
      USE mod_stepping
!
      USE analytical_mod
!
#ifdef TANGENT
      USE ad_post_initial_mod, ONLY : ad_post_initial
#endif
      USE close_io_mod,        ONLY : close_inp, close_out
      USE dateclock_mod,       ONLY : time_string
#ifdef DISTRIBUTE
      USE distribute_mod,      ONLY : mp_bcasti
#endif
#ifdef WET_DRY
      USE get_wetdry_mod,      ONLY : get_wetdry
#endif
      USE ini_hmixcoef_mod,    ONLY : ini_hmixcoef
      USE inp_par_mod,         ONLY : inp_par
#ifdef NESTING
      USE nesting_mod,         ONLY : nesting
#endif
#ifdef SOLVE3D
      USE set_depth_mod,       ONLY : set_depth0, set_depth
      USE omega_mod,           ONLY : omega
      USE rho_eos_mod,         ONLY : rho_eos
      USE set_massflux_mod,    ONLY : set_massflux
#endif
      USE post_initial_mod,    ONLY : post_initial
#ifdef MASKING
      USE set_masks_mod,       ONLY : set_masks
#endif
      USE stiffness_mod,       ONLY : stiffness
      USE strings_mod,         ONLY : founderror
      USE stdout_mod,          ONLY : set_stdoutunit, stdout_unit
#ifdef TANGENT
      USE tl_post_initial_mod, ONLY : tl_post_initial
      USE tl_set_depth_mod,    ONLY : tl_bath
#endif
#ifdef WET_DRY
      USE wetdry_mod,          ONLY : wetdry
#endif
      USE wrt_rst_mod,         ONLY : wrt_rst
!
      implicit none
!
      PUBLIC  :: roms_initializep1
      PUBLIC  :: roms_initializep2
      PUBLIC  :: roms_initializep3
      PUBLIC  :: roms_run
      PUBLIC  :: roms_finalize
!
      PRIVATE :: nlm_initial
#ifdef TANGENT
      PRIVATE :: tlm_initial
#endif
#ifdef ADJOINT
      PRIVATE :: adm_initial
#endif
!
      CONTAINS
!
      SUBROUTINE roms_initializep1 (first, mpiCOMM, kernel)
!
!=======================================================================
!                                                                      !
!  ROMSJEDI phase 2 initialization. It allocates and initializes ROMS  !
!  state variables and internal parameters. It reads standard input    !
!  parameters.                                                         !
!                                                                      !
!=======================================================================
!
!  Imported variable declarations.
!
      logical, intent(inout) :: first
!
      integer, intent(in), optional :: mpicomm
      integer, intent(in), optional :: kernel
!
!  Local variable declarations.
!
      logical :: allocate_vars = .true.
 
#ifdef DISTRIBUTE
      integer :: myerror, mysize
#endif
      integer :: phase, chunk_size, ng, thread, tile
!
      character (len=*), parameter :: myfile =                          &
     &  __FILE__//", ROMS_initialize"
 
#ifdef DISTRIBUTE
!
!-----------------------------------------------------------------------
!  Set distribute-memory (mpi) world communictor.
!-----------------------------------------------------------------------
!
      IF (PRESENT(mpicomm)) THEN
        ocn_comm_world=mpicomm
      ELSE
        ocn_comm_world=mpi_comm_world
      END IF
      CALL mpi_comm_rank (ocn_comm_world, myrank, myerror)
      CALL mpi_comm_size (ocn_comm_world, mysize, myerror)
#endif
!
!-----------------------------------------------------------------------
!  On first pass, initialize model parameters a variables for all
!  nested/composed grids.  Notice that the logical switch "first"
!  is used to allow multiple calls to this routine during ensemble
!  configurations.
!-----------------------------------------------------------------------
!
      IF (first) THEN
        first=.false.
!
!  Initialize parallel control switches. These scalars switches are
!  independent from standard input parameters.
!
        CALL initialize_parallel
!
!  Set the ROMS standard output unit to write verbose execution info.
!  Notice that the default standard out unit in Fortran is 6.
!
!  In some applications like coupling or disjointed mpi-communications,
!  it is advantageous to write standard output to a specific filename
!  instead of the default Fortran standard output unit 6. If that is
!  the case, it opens such formatted file for writing.
!
        IF (set_stdoutunit) THEN
          stdout=stdout_unit(master)
          set_stdoutunit=.false.
        END IF
!
!  Read in model tunable parameters from standard input. Allocate and
!  initialize variables in several modules after the number of nested
!  grids and dimension parameters are known.
!
        CALL inp_par (inlm)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Initialize counters.
!
        nrun=1                ! run counter
!
!  Set domain decomposition tile partition range.  This range is
!  computed only once since the "first_tile" and "last_tile" values
!  are private for each parallel thread/node.
!
#if defined DISTRIBUTE
        mythread=myrank
#else
        mythread=0
#endif
        DO ng=1,ngrids
          chunk_size=(ntilex(ng)*ntilee(ng)+numthreads-1)/numthreads
          first_tile(ng)=mythread*chunk_size
          last_tile(ng)=first_tile(ng)+chunk_size-1
        END DO
!
!  Initialize internal wall clocks. Notice that the timings does not
!  includes processing standard input because several parameters are
!  needed to allocate clock variables.
!
        IF (master) THEN
          WRITE (stdout,10)
 10       FORMAT (/,' Process Information:',/)
        END IF
!
        DO ng=1,ngrids
          DO thread=thread_range
            CALL wclock_on (ng, inlm, 0, __line__, myfile)
          END DO
        END DO
!
!  Allocate and initialize modules variables.
!
        CALL roms_allocate_arrays (allocate_vars)
        CALL roms_initialize_arrays
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
      END IF
!
!-----------------------------------------------------------------------
!  Initialize ROMS, phase 1.
!-----------------------------------------------------------------------
!
      phase=1
!
      SELECT CASE (kernel)
        CASE (inlm)
          CALL nlm_initial (phase)
#ifdef TANGENT
        CASE (itlm)
          CALL tlm_initial (phase)
#endif
#ifdef ADJOINT
        CASE (iadm)
          CALL adm_initial (phase)
#endif
      END SELECT
!
!-----------------------------------------------------------------------
!  Set ROMS application grid configuration. It is done once.
!-----------------------------------------------------------------------
!
      DO ng=1,ngrids
        IF (setgridconfig(ng)) THEN
          CALL set_grid (ng, inlm)
          setgridconfig(ng)=.false.
        END IF
      END DO
!
      RETURN
      END SUBROUTINE roms_initializep1
!
      SUBROUTINE roms_initializep2 (kernel)
!
!=======================================================================
!                                                                      !
!  ROMSJEDI phase 2 initialization. It requires the initial state      !
!  (set elsewhere) to complete the full initialization of the kernel.  !
!  It computes depths, density, and horizontal mass fluxes.            !
!                                                                      !
!=======================================================================
!
!  Imported variable declarations
!
      integer, intent(in) :: kernel
!
!  Local variable declarations.
!
      integer :: phase
!
!-----------------------------------------------------------------------
!  Initialize ROMSJEDI, Phase 2.
!-----------------------------------------------------------------------
!
      phase=2
!
      SELECT CASE (kernel)
        CASE (inlm)
          CALL nlm_initial (phase)
#ifdef TANGENT
        CASE (itlm)
          CALL tlm_initial (phase)
#endif
#ifdef ADJOINT
        CASE (iadm)
          CALL adm_initial (phase)
#endif
      END SELECT
!
      RETURN
      END SUBROUTINE roms_initializep2
!
      SUBROUTINE roms_initializep3 (kernel)
!
!=======================================================================
!                                                                      !
!  ROMSJEDI phase 3 initialization. It computes the initial depths     !
!  and level thicknesses from the time-averaged free-surface, Zt_avg1, !
!  which are needed to calculate initial vertically integrated         !
!  momentum.                                                           !
!                                                                      !
!  Also, it initializes the state variables for all time levels and    !
!  applies lateral boundary conditions.                                !
!                                                                      !
!  Notice that the second data snapshot needs to be processed because  !
!  it is necessary for the initial lateral boundary conditions.        !
!                                                                      !
!=======================================================================
!
!  Imported variable declarations
!
      integer, intent(in) :: kernel
!
!  Local variable declarations.
!
      integer :: phase
!
!-----------------------------------------------------------------------
!  Initialize ROMSJEDI, Phase 3.
!-----------------------------------------------------------------------
!
      phase=3
!
      SELECT CASE (kernel)
        CASE (inlm)
          CALL nlm_initial (phase)
#ifdef TANGENT
        CASE (itlm)
          CALL tlm_initial (phase)
#endif
#ifdef ADJOINT
        CASE (iadm)
          CALL adm_initial (phase)
#endif
      END SELECT
!
      RETURN
      END SUBROUTINE roms_initializep3
!
      SUBROUTINE roms_run (RunInterval, kernel)
!
!=======================================================================
!                                                                      !
!  This routine advance ROMS kernel for the specified time interval.   !
!                                                                      !
!  On Input:                                                           !
!                                                                      !
!     RunInterval     Execution time stepping window (seconds)         !
!     kernel          Dynamical/numerical kernel (integer)             !
!                                                                      !
!=======================================================================
!
!  Imported variable declarations
!
      integer, intent(in), optional :: kernel
!
      real(dp), intent(in) :: runinterval
!
!  Local variable declarations.
!
      integer :: ng
      integer :: nstrstep, nendstep, extra
!
      real(dp) :: endtime, nexttime
!
      character (len=2), dimension(4) :: mid = (/'NL','TL','RP','AD'/)
 
      character (len=*), parameter :: myfile =                          &
     &  __FILE__//", ROMS_run"
!
      sourcefile=myfile
!
!=======================================================================
!  Advance ROMS dynamical/numerical kernel: NLM, TLM, or ADM
!=======================================================================
!
!  OOPS will advance ROMS kernels by smaller intervals, usually a single
!  timestep. The strategy here is different to that used for coupling
!  since ROMS delayed output. The delayed ouput last half timestep will
!  affect the OOPS trajectory logic needed to save the NLM background
!  fields needed to linearize the TLM and ADM kernels.
!
      myruninterval=runinterval
      IF (master) WRITE (stdout,'(1x)')
!
      SELECT CASE (kernel)
        CASE (inlm)
          DO ng=1,ngrids
            nexttime=time(ng)+runinterval
            endtime=initime(ng)+ntimes(ng)*dt(ng)
            extra=1
            step_counter(ng)=0
            nstrstep=iic(ng)-1
            nendstep=nstrstep+int((myruninterval)/dt(ng))-extra
            IF (iic(ng).eq.ntstart(ng)) THEN
              processinputdata(ng)=.false.         ! because Phase 3
            ELSE
              processinputdata(ng)=.true.
            END IF
            IF (master) WRITE (stdout,10) mid(kernel), ng,              &
     &                                    nstrstep, max(0,nendstep)
          END DO
          IF (master) WRITE (stdout,'(1x)')
        CASE (itlm)
          DO ng=1,ngrids
            nexttime=time(ng)+runinterval
            endtime=initime(ng)+ntimes(ng)*dt(ng)
            extra=1
            step_counter(ng)=0
            nstrstep=iic(ng)-1
            nendstep=nstrstep+int((myruninterval)/dt(ng))-extra
            IF (iic(ng).eq.ntstart(ng)) THEN
              processinputdata(ng)=.false.         ! because Phase 3
            ELSE
              processinputdata(ng)=.true.
            END IF
            IF (master) WRITE (stdout,10) mid(kernel), ng,              &
     &                                    nstrstep, max(0,nendstep)
          END DO
          IF (master) WRITE (stdout,'(1x)')
        CASE (iadm)
          DO ng=1,ngrids
            nexttime=time(ng)-runinterval
            endtime=initime(ng)+ntimes(ng)*dt(ng)
            extra=1
            step_counter(ng)=0
            nstrstep=iic(ng)-1
            nendstep=nstrstep-int((myruninterval)/dt(ng))+extra
            IF (master) WRITE (stdout,10) mid(kernel), ng,              &
     &                                    nstrstep, max(0,nendstep)
          END DO
          IF (master) WRITE (stdout,'(1x)')
      END SELECT
!
!  Time-step ROMS kernel.
!
      SELECT CASE (kernel)
        CASE (inlm)
          CALL main3d (runinterval)
#ifdef TANGENT
        CASE (itlm)
          CALL tl_main3d (runinterval)
#endif
#ifdef ADJOINT
        CASE (iadm)
          CALL ad_main3d (runinterval)
#endif
      END SELECT
!
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
 10   FORMAT (1x,a,1x,'ROMS: started time-stepping:',                   &
     &        ' (Grid: ',i2.2,' TimeSteps: ',i0,' - ',i0,')')
!
      RETURN
      END SUBROUTINE roms_run
!
      SUBROUTINE roms_finalize
!
!=======================================================================
!                                                                      !
!  This routine terminates ROMS W4D-RBLanczos execution.               !
!                                                                      !
!=======================================================================
!
!  Local variable declarations.
!
      integer :: fcount, ng, thread
!
      character (len=*), parameter :: myfile =                          &
     &  __FILE__//", ROMS_finalize"
!
!-----------------------------------------------------------------------
!  If blowing-up, save latest model state into RESTART NetCDF file.
!-----------------------------------------------------------------------
!
!  If cycling restart records, write solution into record 3.
!
      IF (exit_flag.eq.1) THEN
        DO ng=1,ngrids
          IF (lwrtrst(ng)) THEN
            IF (master) WRITE (stdout,10)
 10         FORMAT (/,' Blowing-up: Saving latest model state into ',   &
     &                ' RESTART file',/)
            fcount=rst(ng)%load
            IF (lcyclerst(ng).and.(rst(ng)%Nrec(fcount).ge.2)) THEN
              rst(ng)%Rindex=2
              lcyclerst(ng)=.false.
            END IF
            blowup=exit_flag
            exit_flag=noerror
#ifdef DISTRIBUTE
            CALL wrt_rst (ng, myrank)
#else
            CALL wrt_rst (ng, -1)
#endif
          END IF
        END DO
      END IF
!
!-----------------------------------------------------------------------
!  Stop model and time profiling clocks, report memory requirements,
!  and close output NetCDF files.
!-----------------------------------------------------------------------
!
!  Stop time clocks.
!
      IF (master) THEN
        WRITE (stdout,20)
 20     FORMAT (/,'Elapsed wall CPU time for each process (seconds):',/)
      END IF
!
      DO ng=1,ngrids
        DO thread=thread_range
          CALL wclock_off (ng, inlm, 0, __line__, myfile)
        END DO
      END DO
!
!  Report dynamic memory and automatic memory requirements.
!
      CALL memory
!
!  Close IO files.
!
      DO ng=1,ngrids
        CALL close_inp (ng, inlm)
      END DO
      CALL close_out
!
      RETURN
      END SUBROUTINE roms_finalize
!
      SUBROUTINE nlm_initial (Phase)
!
!=======================================================================
!                                                                      !
!  ROMSJEDI nonlinear model (NLM) kernel initialization Phases:        !
!                                                                      !
!     Phase 1: Set time-stepping parameters                            !
!                                                                      !
!     Phase 2: Computes initial depths, density, horizontal mass       !
!              fluxes, and other configuration arrays. It reads        !
!              forcing snapshots.                                      !
!                                                                      !
!     Phase 3: Computes the initial depths and level thicknesses from  !
!              the initial time-averaged free-surface field, Zt_avg1.  !
!              Additionally, it initializes the state variables for    !
!              all time levels and applies lateral boundary conditions.!
!                                                                      !
!=======================================================================
!
!  Imported variable declarations
!
      integer, intent(in) :: phase
!
!  Local variable declarations.
!
      integer :: fcount
      integer :: ng, thread, tile
#ifdef NESTING
      integer :: ig, nl
      integer :: cr, i, m
#endif
      integer, dimension(Ngrids) :: inirec, tindex
!
      character (len=*), parameter :: myfile =                          &
     &  __FILE__//", nlm_initial"
 
#ifdef PROFILE
!
!-----------------------------------------------------------------------
!  Start time wall clocks.
!-----------------------------------------------------------------------
!
      DO ng=1,ngrids
        DO thread=thread_range
          CALL wclock_on (ng, inlm, 2, __line__, myfile)
        END DO
      END DO
#endif
!
!=======================================================================
!  ROMSJEDI NLM kernel, Phase 1 Initialization.
!=======================================================================
!
      phase1 : IF (phase.eq.1) THEN
!
        IF (master) THEN
          WRITE (stdout,10) 'NLM_INITIAL: Phase 1 Initialization: ',    &
     &                      'Configuring ROMS nonlinear kernel ...'
        END IF
!
!  Initialize time stepping indices and counters.
!
        DO ng=1,ngrids
          iif(ng)=1
          indx1(ng)=1
          kstp(ng)=1
          krhs(ng)=1
          knew(ng)=1
          predictor_2d_step(ng)=.false.
!
          iic(ng)=0
# ifdef JEDI
          jic(ng)=0
# endif
          nstp(ng)=1
          nrhs(ng)=1
          nnew(ng)=1
#ifdef FLOATS
          nf(ng)=0
          nfp1(ng)=1
          nfm1(ng)=4
          nfm2(ng)=3
          nfm3(ng)=2
#endif
!
          inirec(ng)=nrrec(ng)
          tindex(ng)=1
!
          synchro_flag(ng)=.true.
          first_time(ng)=0
          tdays(ng)=dstart
          time(ng)=tdays(ng)*day2sec
#ifdef JEDI
          time4jedi(ng)=time(ng)
#endif
          ntstart(ng)=int((time(ng)-dstart*day2sec)/dt(ng))+1
          ntfirst(ng)=ntstart(ng)
          step_counter(ng)=0
        END DO
!
!  Initialize global diagnostics variables.
!
        avgke=0.0_dp
        avgpe=0.0_dp
        avgkp=0.0_dp
        volume=0.0_dp
!
!  Reset output history files time record counters. These counters are
!  reset on every iteration pass. This file is created on the first
!  iteration pass.
!
        DO ng=1,ngrids
          ldefhis(ng)=.true.
          lwrthis(ng)=.true.
          his(ng)%Rindex=0
          fcount=his(ng)%Fcount
          his(ng)%Nrec(fcount)=0
 
          ldefqck(ng)=.true.
          lwrtqck(ng)=.true.
          qck(ng)%Rindex=0
          fcount=qck(ng)%Fcount
          qck(ng)%Nrec(fcount)=0
 
          ldefrst(ng)=.true.
          lwrtrst(ng)=.true.
          rst(ng)%Rindex=0
          fcount=rst(ng)%Fcount
          rst(ng)%Nrec(fcount)=0
        END DO
 
      END IF phase1
!
!=======================================================================
!  ROMSJEDI NLM kernel, Phase 2 Initialization.
!=======================================================================
!
      phase2: IF (phase.eq.2) THEN
!
        IF (master) THEN
          WRITE (stdout,10) 'NLM_INITIAL: Phase 2 Initialization: ',    &
     &                      'Get/Set required applications fields ...'
        END IF
!
!-----------------------------------------------------------------------
!  Initialize horizontal mixing coefficients. If applicable, scale
!  mixing coefficients according to the grid size (smallest area).
#ifndef ANA_SPONGE
!  Also increase their values in sponge areas using the "visc_factor"
!  and/or "diff_factor" read from input Grid NetCDF file.
#endif
!-----------------------------------------------------------------------
!
        DO ng=1,ngrids
          DO tile=first_tile(ng),last_tile(ng),+1
            CALL ini_hmixcoef (ng, tile, inlm)
          END DO
        END DO
 
#ifdef ANA_SPONGE
!
!-----------------------------------------------------------------------
!  Increase horizontal mixing coefficients in sponge areas using
!  analytical functions.
!-----------------------------------------------------------------------
!
        DO ng=1,ngrids
          IF (lsponge(ng)) THEN
            DO tile=first_tile(ng),last_tile(ng),+1
              CALL ana_sponge (ng, tile, inlm)
            END DO
          END IF
        END DO
#endif
 
#ifdef WET_DRY
!
!-----------------------------------------------------------------------
!  Process initial wet/dry masks.
!-----------------------------------------------------------------------
!
        DO ng=1,ngrids
          DO tile=first_tile(ng),last_tile(ng),+1
            CALL wetdry (ng, tile, tindex(ng), .true.)
          END DO
        END DO
#endif
 
#ifdef SOLVE3D
!
!-----------------------------------------------------------------------
!  Compute time independent (Zt_avg1=0) and initial time dependent
!  depths and level thicknesses.
!-----------------------------------------------------------------------
!
        DO ng=1,ngrids
          DO tile=first_tile(ng),last_tile(ng),+1
            CALL set_depth0 (ng, tile, inlm)
            CALL set_depth  (ng, tile, inlm)
          END DO
        END DO
!
!-----------------------------------------------------------------------
!  Compute initial horizontal mass fluxes, Hz*u/n and Hz*v/m.
!-----------------------------------------------------------------------
!
        DO ng=1,ngrids
          DO tile=first_tile(ng),last_tile(ng),+1
            CALL set_massflux (ng, tile, inlm)
          END DO
        END DO
!
!-----------------------------------------------------------------------
!  Compute initial S-coordinates vertical velocity. Compute initial
!  density anomaly from potential temperature and salinity via equation
!  of state for seawater.  Also compute other equation of state related
!  quatities.
!-----------------------------------------------------------------------
!
        DO ng=1,ngrids
          DO tile=first_tile(ng),last_tile(ng),+1
            CALL omega (ng, tile, inlm)
            CALL rho_eos (ng, tile, inlm)
          END DO
        END DO
#endif
 
#ifdef ANA_PSOURCE
!
!-----------------------------------------------------------------------
!  Set point Sources/Sinks position, direction, special flag, and mass
!  transport nondimensional shape profile with analytcal expressions.
!  Point sources are at U- and V-points. We need to get their positions
!  to process internal Land/Sea masking arrays during initialization.
!-----------------------------------------------------------------------
!
        DO ng=1,ngrids
          IF (luvsrc(ng).or.lwsrc(ng).or.any(ltracersrc(:,ng))) THEN
            DO tile=first_tile(ng),last_tile(ng),+1
              CALL ana_psource (ng, tile, inlm)
            END DO
          END IF
        END DO
#endif
!
!-----------------------------------------------------------------------
!  If applicable, close all input boundary, climatology, and forcing
!  NetCDF files and set associated parameters to the closed state. This
!  step is essential in iterative algorithms that run the full TLM
!  repetitively. Then, Initialize several parameters in their file
!  structure, so the appropriate input single or multi-file is selected
!  during initialization/restart.
!-----------------------------------------------------------------------
!
        DO ng=1,ngrids
          CALL close_inp (ng, inlm)
          CALL check_multifile (ng, inlm)
#ifdef DISTRIBUTE
          CALL mp_bcasti (ng, inlm, exit_flag)
#endif
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END DO
!
!  If applicable, read in input data.
!
        DO ng=1,ngrids
          CALL get_idata (ng)
          CALL get_data (ng)
#ifdef DISTRIBUTE
          CALL mp_bcasti (ng, inlm, exit_flag)
#endif
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END DO
 
#ifdef MASKING
!
!-----------------------------------------------------------------------
!  Set internal I/O mask arrays.
!-----------------------------------------------------------------------
!
        DO ng=1,ngrids
          DO tile=first_tile(ng),last_tile(ng),+1
            CALL set_masks (ng, tile, inlm)
          END DO
        END DO
#endif
 
#ifdef NESTING
# if defined MASKING || defined WET_DRY
!
!-----------------------------------------------------------------------
!  If nesting and Land/Sea masking, scale horizontal interpolation
!  weights to account for land contact points.
!-----------------------------------------------------------------------
!
        DO ng=1,ngrids
          CALL nesting (ng, inlm, nmask)
        END DO
# endif
!
!-----------------------------------------------------------------------
!  If nesting, process state fields initial conditions in the contact
!  regions.
!-----------------------------------------------------------------------
!
!  Free-surface and 2D-momentum.
!
        DO nl=1,nestlayers
          DO ig=1,gridsinlayer(nl)
            ng=gridnumber(ig,nl)
            IF (any(compositegrid(:,ng))) THEN
              CALL nesting (ng, inlm, nfsic)      ! free-surface
# ifndef SOLVE3D
              CALL nesting (ng, inlm, n2dic)      ! 2d momentum
# endif
            END IF
          END DO
        END DO
 
# ifdef SOLVE3D
!
!  Determine vertical indices and vertical interpolation weights in
!  the contact zone using initial unperturbed depth arrays.
!
        DO ng=1,ngrids
          CALL nesting (ng, inlm, nzwgt)
        END DO
!
!  3D-momentum and tracers.
!
        DO nl=1,nestlayers
          DO ig=1,gridsinlayer(nl)
            ng=gridnumber(ig,nl)
            IF (any(compositegrid(:,ng))) THEN
              CALL nesting (ng, inlm, n3dic)      ! 3D momentum
              CALL nesting (ng, inlm, ntvic)      ! Tracer variables
            END IF
          END DO
        END DO
# endif
#endif
 
#if defined ANA_DRAG && defined UV_DRAG_GRID
!
!-----------------------------------------------------------------------
!  Set analytical spatially varying bottom friction parameter.
!-----------------------------------------------------------------------
!
        IF (nrun.eq.erstr) THEN
          DO ng=1,ngrids
            DO tile=first_tile(ng),last_tile(ng),+1
              CALL ana_drag (ng, tile, inlm)
            END DO
          END DO
        END IF
#endif
!
!-----------------------------------------------------------------------
!  Compute grid stiffness.
!-----------------------------------------------------------------------
!
        IF (lstiffness) THEN
          lstiffness=.false.
          DO ng=1,ngrids
            DO tile=first_tile(ng),last_tile(ng),+1
              CALL stiffness (ng, tile, inlm)
            END DO
          END DO
        END IF
 
#if defined FLOATS || defined STATIONS
!
!-----------------------------------------------------------------------
!  If applicable, convert initial locations to fractional grid
!  coordinates.
!-----------------------------------------------------------------------
!
        DO ng=1,ngrids
          CALL grid_coords (ng, inlm)
        END DO
#endif
      END IF phase2
!
!=======================================================================
!  ROMSJEDI NLM kernel, Phase 3 Initialization.
!=======================================================================
!
      phase3: IF (phase.eq.3) THEN
!
        IF (master) THEN
          WRITE (stdout,10) 'NLM_INITIAL: Phase 3 Initialization: ',    &
     &                      'State Lateral Boundary Conditions ...'
        END IF
!
!-----------------------------------------------------------------------
!  Initialize time-stepping counter and time/date string. Save NLM
!  initial conditions time.
!
!  Notice that it is allowed to modify the "simulation length" in the
!  roms-jedi YAML file, which will affect the values of "ntimes" and
!  "ntend".
!-----------------------------------------------------------------------
!
        DO ng=1,ngrids
          initime(ng)=time(ng)
          iic(ng)=ntstart(ng)
          ntend(ng)=ntstart(ng)+ntimes(ng)-1
#ifdef JEDI
          jic(ng)=ntstart(ng)
          time4jedi(ng)=time(ng)
#endif
          CALL time_string (time(ng), time_code(ng))
        END DO
!
!-----------------------------------------------------------------------
!  Read in required data, if any, from input NetCDF files.
!-----------------------------------------------------------------------
!
        DO ng=1,ngrids
          CALL get_data (ng)
          IF (founderror(exit_flag, noerror,                            &
     &                   __line__, myfile)) RETURN
        END DO
!
!-----------------------------------------------------------------------
!  If applicable, process input data: time interpolate between data
!  snapshots.
!-----------------------------------------------------------------------
!
        DO ng=1,ngrids
          DO tile=first_tile(ng),last_tile(ng),+1
            CALL set_data (ng, tile)
          END DO
        END DO
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!-----------------------------------------------------------------------
!  Computes the initial depths and level thicknesses from the initial
!  time=averaged free-surface field, Zt_avg1 . Additionally, it
!  initializes the nonlinear state variables for all time levels and
!  applies lateral boundary conditions.
!-----------------------------------------------------------------------
!
        DO ng=1,ngrids
          CALL post_initial (ng, inlm)
        END DO
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
      END IF phase3
 
#ifdef PROFILE
!
!-----------------------------------------------------------------------
!  Turn off initialization time wall clock.
!-----------------------------------------------------------------------
!
      DO ng=1,ngrids
        DO thread=thread_range
          CALL wclock_off (ng, inlm, 2, __line__, myfile)
        END DO
      END DO
#endif
!
  10  FORMAT (/,1x,a,a,/,1x,'***********',/)
!
      RETURN
      END SUBROUTINE nlm_initial
 
#ifdef TANGENT
!
      SUBROUTINE tlm_initial (Phase)
!
!=======================================================================
!                                                                      !
!  ROMSJEDI tangent linear model (TLM) kernel initialization Phases:   !
!                                                                      !
!     Phase 1: Set time-stepping parameters                            !
!                                                                      !
!     Phase 2: Computes masks and other configuration arrays. It reads !
!              initial forcing snapshots.                              !
!                                                                      !
!=======================================================================
!
!  Imported variable declarations
!
      integer, intent(in) :: phase
!
!  Local variable declarations.
!
      integer :: fcount
      integer :: ng, thread, tile
      integer :: lstr, ifile
 
      integer, dimension(Ngrids) :: inirec, tindex
 
# ifdef GENERIC_DSTART
!
      real(dp) :: my_dstart
# endif
!
      character (len=10) :: suffix
 
      character (len=*), parameter :: myfile =                          &
     &  __FILE__//", tlm_initial"
 
# ifdef PROFILE
!
!-----------------------------------------------------------------------
!  Start time wall clocks.
!-----------------------------------------------------------------------
!
      DO ng=1,ngrids
        DO thread=thread_range
          CALL wclock_on (ng, itlm, 2, __line__, myfile)
        END DO
      END DO
# endif
!
!=======================================================================
!  ROMSJEDI TLM kernel, Phase 1 Initialization.
!=======================================================================
!
      phase1 : IF (phase.eq.1) THEN
!
        IF (master) THEN
          WRITE (stdout,10) 'TLM_INITIAL: Phase 1 Initialization: ',    &
     &                      'Configuring ROMS tangent linear kernel ...'
        END IF
!
!  Initialize time stepping indices and counters.
!
        DO ng=1,ngrids
          iif(ng)=1
          indx1(ng)=1
          kstp(ng)=1
          krhs(ng)=1
          knew(ng)=1
          predictor_2d_step(ng)=.false.
!
          iic(ng)=0
# ifdef JEDI
          jic(ng)=0
# endif
          nstp(ng)=1
          nrhs(ng)=1
          nnew(ng)=1
# ifdef FLOATS_NOT_YET
          nf(ng)=0
          nfp1(ng)=1
          nfm1(ng)=4
          nfm2(ng)=3
          nfm3(ng)=2
# endif
!
          synchro_flag(ng)=.true.
          first_time(ng)=0
          IF (any(tl_volcons(:,ng))) THEN
            tl_ubar_xs=0.0_r8
          END IF
 
# if defined GENERIC_DSTART
!
!  Rarely, the tangent linear model is initialized from a NetCDF file,
!  so we do not know its actual initialization time. Usually, it is
!  computed from DSTART, implying that its value is correct in the ROMS
!  input script. Therefore, the user needs to check and update its value
!  to every time that ROMS is executed. Alternatively, if available, we
!  can use the initialization time from the nonlinear model, INItime.
!  This variable is assigned when computing or processing the basic
!  state trajectory needed to linearize the adjoint model.
!
          IF (initime(ng).lt.0.0_dp) THEN
            my_dstart=dstart                  ! ROMS input script
          ELSE
            my_dstart=initime(ng)/86400.0_dp  ! NLM IC time is known
          END IF
          tdays(ng)=my_dstart
# else
          tdays(ng)=dstart
# endif
          time(ng)=tdays(ng)*day2sec
# ifdef JEDI
          time4jedi(ng)=time(ng)
# endif
          ntstart(ng)=int((time(ng)-dstart*day2sec)/dt(ng))+1
          ntend(ng)=ntstart(ng)+ntimes(ng)-1
          ntfirst(ng)=ntstart(ng)
!
          inirec(ng)=nrrec(ng)
          tindex(ng)=1
        END DO
!
!  Initialize global diagnostics variables.
!
        avgke=0.0_dp
        avgpe=0.0_dp
        avgkp=0.0_dp
        volume=0.0_dp
 
      END IF phase1
!
!=======================================================================
!  ROMSJEDI TLM kernel, Phase 2 Initialization.
!=======================================================================
!
      phase2: IF (phase.eq.2) THEN
!
        IF (master) THEN
          WRITE (stdout,10) 'TLM_INITIAL: Phase 2 Initialization: ',   &
     &                      'Get/Set required applications fields ...'
        END IF
!
!-----------------------------------------------------------------------
!  Initialize horizontal mixing coefficients. If applicable, scale
!  mixing coefficients according to the grid size (smallest area).
# ifndef ANA_SPONGE
!  Also increase their values in sponge areas using the "visc_factor"
!  and/or "diff_factor" read from input Grid NetCDF file.
# endif
!-----------------------------------------------------------------------
!
        DO ng=1,ngrids
          DO tile=first_tile(ng),last_tile(ng),+1
            CALL ini_hmixcoef (ng, tile, itlm)
          END DO
        END DO
 
# ifdef ANA_SPONGE
!
!-----------------------------------------------------------------------
!  Increase horizontal mixing coefficients in sponge areas using
!  analytical functions.
!-----------------------------------------------------------------------
!
        DO ng=1,ngrids
          IF (lsponge(ng)) THEN
            DO tile=first_tile(ng),last_tile(ng),+1
              CALL ana_sponge (ng, tile, itlm)
            END DO
          END IF
        END DO
# endif
!
!-----------------------------------------------------------------------
!  Initialize tangent linear bathymetry to zero.
!-----------------------------------------------------------------------
!
        DO ng=1,ngrids
          DO tile=first_tile(ng),last_tile(ng),+1
            CALL tl_bath (ng, tile)
          END DO
        END DO
 
# ifdef WET_DRY
!
!-----------------------------------------------------------------------
!  Process initial wet/dry masks.
!-----------------------------------------------------------------------
!
        DO ng=1,ngrids
          DO tile=first_tile(ng),last_tile(ng),+1
            CALL wetdry (ng, tile, tindex(ng), .true.)
          END DO
        END DO
# endif
 
# ifdef FORWARD_FLUXES
!
!-----------------------------------------------------------------------
!  Set the BLK structure to contain the nonlinear model surface fluxes
!  needed by the tangent linear and adjoint models. Also, set switches
!  to process that structure in routine "check_multifile". Notice that
!  it is possible to split the solution into multiple NetCDF files to
!  reduce their size.
!-----------------------------------------------------------------------
!
!  Set the nonlinear model quicksave-history file as the basic state for
!  the surface fluxes computed in "bulk_flux", which may be available at
!  more frequent intervals while avoiding large files. Since the 4D-Var
!  increment phase is calculated by a different executable and needs to
!  know some of the QCK structure information.
!
      DO ng=1,ngrids
        WRITE (qck(ng)%name,"(a,'.nc')") trim(qck(ng)%head)
        lstr=len_trim(qck(ng)%name)
        qck(ng)%base=qck(ng)%name(1:lstr-3)
        IF (qck(ng)%Nfiles.gt.1) THEN
          DO ifile=1,qck(ng)%Nfiles
            WRITE (suffix,"('_',i4.4,'.nc')") ifile
            qck(ng)%files(ifile)=trim(qck(ng)%base)//trim(suffix)
          END DO
          qck(ng)%name=trim(qck(ng)%files(1))
        ELSE
          qck(ng)%files(1)=trim(qck(ng)%name)
        END IF
      END DO
!
!  The switch LreadFRC is deactivated because all the atmospheric
!  forcing, including shortwave radiation, is read from the NLM
!  surface fluxes or is assigned during ESM coupling.  Such fluxes
!  are available from the QCK structure. There is no need for reading
!  and processing from the FRC structure input forcing-files.
!
      CALL edit_multifile ('QCK2BLK')
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      DO ng=1,ngrids
        lreadblk(ng)=.true.
        lreadfrc(ng)=.false.
        lreadqck(ng)=.false.
      END DO
# endif
!
!-----------------------------------------------------------------------
!  If applicable, close all input boundary, climatology, and forcing
!  NetCDF files and set associated parameters to the closed state. This
!  step is essential in iterative algorithms that run the full TLM
!  repetitively. Then, Initialize several parameters in their file
!  structure, so the appropriate input single or multi-file is selected
!  during initialization/restart.
!-----------------------------------------------------------------------
!
        DO ng=1,ngrids
          CALL close_inp (ng, itlm)
          CALL check_multifile (ng, itlm)
# ifdef DISTRIBUTE
          CALL mp_bcasti (ng, itlm, exit_flag)
# endif
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END DO
!
!-----------------------------------------------------------------------
!  Read in initial forcing, climatology and assimilation data from
!  input NetCDF files.  It loads the first relevant data record for
!  the time-interpolation between snapshots.
!-----------------------------------------------------------------------
!
        DO ng=1,ngrids
          CALL tl_get_idata (ng)
          CALL tl_get_data (ng)
# ifdef DISTRIBUTE
          CALL mp_bcasti (ng, itlm, exit_flag)
# endif
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END DO
 
# ifdef MASKING
!
!-----------------------------------------------------------------------
!  Set internal I/O mask arrays.
!-----------------------------------------------------------------------
!
        DO ng=1,ngrids
          DO tile=first_tile(ng),last_tile(ng),+1
            CALL set_masks (ng, tile, itlm)
          END DO
        END DO
# endif
 
# if defined ANA_DRAG && defined UV_DRAG_GRID
!
!-----------------------------------------------------------------------
!  Set analytical spatially varying bottom friction parameter.
!-----------------------------------------------------------------------
!
        IF (nrun.eq.erstr) THEN
          DO ng=1,ngrids
            DO tile=first_tile(ng),last_tile(ng),+1
              CALL ana_drag (ng, tile, itlm)
            END DO
          END DO
        END IF
# endif
!
!-----------------------------------------------------------------------
!  Compute grid stiffness.
!-----------------------------------------------------------------------
!
        IF (lstiffness) THEN
          lstiffness=.false.
          DO ng=1,ngrids
            DO tile=first_tile(ng),last_tile(ng),+1
              CALL stiffness (ng, tile, itlm)
            END DO
          END DO
        END IF
 
# if defined FLOATS_NOT_YET || defined STATIONS
!
!-----------------------------------------------------------------------
!  If applicable, convert initial locations to fractional grid
!  coordinates.
!-----------------------------------------------------------------------
!
        DO ng=1,ngrids
          CALL grid_coords (ng, itlm)
        END DO
# endif
      END IF phase2
!
!=======================================================================
!  ROMSJEDI TLM kernel, Phase 3 Initialization.
!=======================================================================
!
      phase3: IF (phase.eq.3) THEN
!
        IF (master) THEN
          WRITE (stdout,10) 'TLM_INITIAL: Phase 3 Initialization: ',    &
     &                      'State Lateral Boundary Conditions ...'
        END IF
!
!-----------------------------------------------------------------------
!  Initialize time-stepping counter and time/date string.
!
!  Notice that it is allowed to modify the "simulation length" in the
!  roms-jedi YAML file, which will affect the values of "ntimes" and
!  "ntend".
!-----------------------------------------------------------------------
!
!  Subsract one time unit to avoid special case due to initialization
!  in the main time-stepping routine.
!
        DO ng=1,ngrids
          iic(ng)=ntstart(ng)
          ntend(ng)=ntstart(ng)+ntimes(ng)-1
#ifdef JEDI
          jic(ng)=ntstart(ng)
          time4jedi(ng)=time(ng)
#endif
          CALL time_string (time(ng), time_code(ng))
          ldeftlm(ng)=.true.
          lwrttlm(ng)=.true.
        END DO
!
!-----------------------------------------------------------------------
!  Read in required data, if any, from input NetCDF files.
!-----------------------------------------------------------------------
!
        DO ng=1,ngrids
          CALL tl_get_data (ng)
          IF (founderror(exit_flag, noerror,                            &
     &                   __line__, myfile)) RETURN
        END DO
!
!-----------------------------------------------------------------------
!  If applicable, process input data: time interpolate between data
!  snapshots.
!-----------------------------------------------------------------------
!
        DO ng=1,ngrids
          DO tile=first_tile(ng),last_tile(ng),+1
            CALL tl_set_data (ng, tile)
          END DO
        END DO
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!-----------------------------------------------------------------------
!  Computes the initial depths and level thicknesses from the initial
!  time-averaged free-surface field, Zt_avg1. Additionally, it
!  initializes the nonlinear state variables for all time levels and
!  applies lateral boundary conditions.
!-----------------------------------------------------------------------
!
        DO ng=1,ngrids
          CALL tl_post_initial (ng, itlm)
        END DO
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
      END IF phase3
 
# ifdef PROFILE
!
!-----------------------------------------------------------------------
!  Turn off initialization time wall clock.
!-----------------------------------------------------------------------
!
      DO ng=1,ngrids
        DO thread=thread_range
          CALL wclock_off (ng, itlm, 2, __line__, myfile)
        END DO
      END DO
# endif
!
  10  FORMAT (/,1x,a,a,/,1x,'***********',/)
!
      RETURN
      END SUBROUTINE tlm_initial
#endif
 
#ifdef ADJOINT
!
      SUBROUTINE adm_initial (Phase)
!
!=======================================================================
!                                                                      !
!  ROMSJEDI adjoint model (ADM) kernel initialization Phases:          !
!                                                                      !
!     Phase 1: Set time-stepping parameters                            !
!                                                                      !
!     Phase 2: Computes masks and other configuration arrays. It reads !
!              initial forcing snapshots.                              !
!                                                                      !
!=======================================================================
!
!  Imported variable declarations
!
      integer, intent(in) :: phase
!
!  Local variable declarations.
!
      integer :: fcount
      integer :: ng, thread, tile
      integer :: lstr, ifile
 
      integer, dimension(Ngrids) :: inirec, tindex
 
# ifdef GENERIC_DSTART
!
      real(dp) :: my_dstart
# endif
!
      character (len=10) :: suffix
 
      character (len=*), parameter :: myfile =                          &
     &  __FILE__//", adm_initial"
 
# ifdef PROFILE
!
!-----------------------------------------------------------------------
!  Start time wall clocks.
!-----------------------------------------------------------------------
!
      DO ng=1,ngrids
        DO thread=thread_range
          CALL wclock_on (ng, iadm, 2, __line__, myfile)
        END DO
      END DO
# endif
!
!=======================================================================
!  ROMSJEDI ADM kernel, Phase 1 Initialization.
!=======================================================================
!
      phase1 : IF (phase.eq.1) THEN
!
        IF (master) THEN
          WRITE (stdout,10) 'ADM_INITIAL: Phase 1 Initialization: ',    &
     &                      'Configuring ROMS nonlinear kernel ...'
        END IF
!
!  Initialize time stepping indices and counters.
!
        DO ng=1,ngrids
          iif(ng)=1
          indx1(ng)=1
          kstp(ng)=1
          krhs(ng)=3
          knew(ng)=2
          predictor_2d_step(ng)=.false.
!
          iic(ng)=0
# ifdef JEDI
          jic(ng)=0
# endif
# ifdef SOLVE3D
          nstp(ng)=1
          nnew(ng)=2
          nrhs(ng)=nstp(ng)
# endif
# ifdef FLOATS_NOT_YET
          nf(ng)=0
          nfp1(ng)=1
          nfm1(ng)=4
          nfm2(ng)=3
          nfm3(ng)=2
# endif
!
          synchro_flag(ng)=.true.
          first_time(ng)=0
          ad_ubar_xs=0.0_r8
 
# ifdef GENERIC_DSTART
!
!  Rarely, the adjoint model is initialized from a NetCDF file, so we do
!  not know its actual initialization time. Usually, it is computed from
!  DSTART, implying that its value is correct in the ROMS input script.
!  Therefore, the user needs to check and update its value to every time
!  that ROMS is executed. Alternatively, if available, we can use the
!  initialization time from the nonlinear model, INItime. This variable
!  is assigned when computing or processing the basic state trajectory
!  needed to linearize the adjoint model.
!
          IF (initime(ng).lt.0.0_dp) THEN
            my_dstart=dstart                ! ROMS input script
          ELSE
            my_dstart=initime(ng)/86400.0_dp    ! NLM IC time is known
          END IF
          tdays(ng)=my_dstart+dt(ng)*real(ntimes(ng),r8)*sec2day
          time(ng)=tdays(ng)*day2sec
          ntstart(ng)=int((time(ng)-dstart*day2sec)/dt(ng))+1
          ntend(ng)=ntstart(ng)-ntimes(ng)
          ntfirst(ng)=ntend(ng)
# else
          tdays(ng)=dstart+                                             &
     &              dt(ng)*real(ntimes(ng)-ntfirst(ng)+1,r8)*sec2day
          time(ng)=tdays(ng)*day2sec
          ntstart(ng)=ntimes(ng)+1
          ntend(ng)=ntfirst(ng)
          ntfirst(ng)=ntend(ng)
# endif
# ifdef JEDI
          time4jedi(ng)=time(ng)
# endif
          inirec(ng)=nrrec(ng)
          tindex(ng)=1
        END DO
!
!  Initialize global diagnostics variables.
!
        avgke=0.0_dp
        avgpe=0.0_dp
        avgkp=0.0_dp
        volume=0.0_dp
 
      END IF phase1
!
!=======================================================================
!  ROMSJEDI ADM kernel, Phase 2 Initialization.
!=======================================================================
!
      phase2: IF (phase.eq.2) THEN
!
        IF (master) THEN
          WRITE (stdout,10) 'ADM_INITIAL: Phase 2 Initialization: ',    &
     &                      'Get/Set required applications fields ...'
        END IF
!
!-----------------------------------------------------------------------
!  Initialize horizontal mixing coefficients. If applicable, scale
!  mixing coefficients according to the grid size (smallest area).
# ifndef ANA_SPONGE
!  Also increase their values in sponge areas using the "visc_factor"
!  and/or "diff_factor" read from input Grid NetCDF file.
# endif
!-----------------------------------------------------------------------
!
        DO ng=1,ngrids
          DO tile=first_tile(ng),last_tile(ng),+1
            CALL ini_hmixcoef (ng, tile, iadm)
          END DO
        END DO
 
# ifdef ANA_SPONGE
!
!-----------------------------------------------------------------------
!  Increase horizontal mixing coefficients in sponge areas using
!  analytical functions.
!-----------------------------------------------------------------------
!
        DO ng=1,ngrids
          IF (lsponge(ng)) THEN
            DO tile=first_tile(ng),last_tile(ng),+1
              CALL ana_sponge (ng, tile, iadm)
            END DO
          END IF
        END DO
# endif
 
# ifdef WET_DRY
!
!-----------------------------------------------------------------------
!  Process initial wet/dry masks.
!-----------------------------------------------------------------------
!
        DO ng=1,ngrids
          DO tile=first_tile(ng),last_tile(ng),+1
            CALL wetdry (ng, tile, tindex(ng), .true.)
          END DO
        END DO
# endif
 
# ifdef FORWARD_FLUXES
!
!-----------------------------------------------------------------------
!  Set the BLK structure to contain the nonlinear model surface fluxes
!  needed by the tangent linear and adjoint models. Also, set switches
!  to process that structure in routine "check_multifile". Notice that
!  it is possible to split the solution into multiple NetCDF files to
!  reduce their size.
!-----------------------------------------------------------------------
!
!  Set the nonlinear model quicksave-history file as the basic state for
!  the surface fluxes computed in "bulk_flux", which may be available at
!  more frequent intervals while avoiding large files. Since the 4D-Var
!  increment phase is calculated by a different executable and needs to
!  know some of the QCK structure information.
!
      DO ng=1,ngrids
        WRITE (qck(ng)%name,"(a,'.nc')") trim(qck(ng)%head)
        lstr=len_trim(qck(ng)%name)
        qck(ng)%base=qck(ng)%name(1:lstr-3)
        IF (qck(ng)%Nfiles.gt.1) THEN
          DO ifile=1,qck(ng)%Nfiles
            WRITE (suffix,"('_',i4.4,'.nc')") ifile
            qck(ng)%files(ifile)=trim(qck(ng)%base)//trim(suffix)
          END DO
          qck(ng)%name=trim(qck(ng)%files(1))
        ELSE
          qck(ng)%files(1)=trim(qck(ng)%name)
        END IF
      END DO
!
!  The switch LreadFRC is deactivated because all the atmospheric
!  forcing, including shortwave radiation, is read from the NLM
!  surface fluxes or is assigned during ESM coupling.  Such fluxes
!  are available from the QCK structure. There is no need for reading
!  and processing from the FRC structure input forcing-files.
!
      CALL edit_multifile ('QCK2BLK')
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      DO ng=1,ngrids
        lreadblk(ng)=.true.
        lreadfrc(ng)=.false.
        lreadqck(ng)=.false.
      END DO
# endif
!
!-----------------------------------------------------------------------
!  If applicable, close all input boundary, climatology, and forcing
!  NetCDF files and set associated parameters to the closed state. This
!  step is essential in iterative algorithms that run the full TLM
!  repetitively. Then, Initialize several parameters in their file
!  structure, so the appropriate input single or multi-file is selected
!  during initialization/restart.
!-----------------------------------------------------------------------
!
        DO ng=1,ngrids
          CALL close_inp (ng, iadm)
          CALL check_multifile (ng, iadm)
# ifdef DISTRIBUTE
          CALL mp_bcasti (ng, iadm, exit_flag)
# endif
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END DO
!
!-----------------------------------------------------------------------
!  Read in initial forcing, climatology and assimilation data from
!  input NetCDF files.  It loads the first relevant data record for
!  the time-interpolation between snapshots.
!-----------------------------------------------------------------------
!
        DO ng=1,ngrids
          CALL ad_get_idata (ng)
          CALL ad_get_data (ng)
# ifdef DISTRIBUTE
          CALL mp_bcasti (ng, iadm, exit_flag)
# endif
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END DO
 
# ifdef MASKING
!
!-----------------------------------------------------------------------
!  Set internal I/O mask arrays.
!-----------------------------------------------------------------------
!
        DO ng=1,ngrids
          DO tile=first_tile(ng),last_tile(ng),+1
            CALL set_masks (ng, tile, iadm)
          END DO
        END DO
# endif
 
# if defined ANA_DRAG && defined UV_DRAG_GRID
!
!-----------------------------------------------------------------------
!  Set analytical spatially varying bottom friction parameter.
!-----------------------------------------------------------------------
!
        IF (nrun.eq.erstr) THEN
          DO ng=1,ngrids
            DO tile=first_tile(ng),last_tile(ng),+1
              CALL ana_drag (ng, tile, iadm)
            END DO
          END DO
        END IF
# endif
!
!-----------------------------------------------------------------------
!  Compute grid stiffness.
!-----------------------------------------------------------------------
!
        IF (lstiffness) THEN
          lstiffness=.false.
          DO ng=1,ngrids
            DO tile=first_tile(ng),last_tile(ng),+1
              CALL stiffness (ng, tile, iadm)
            END DO
          END DO
        END IF
 
# if defined FLOATS_NOT_YET || defined STATIONS
!
!-----------------------------------------------------------------------
!  If applicable, convert initial locations to fractional grid
!  coordinates.
!-----------------------------------------------------------------------
!
        DO ng=1,ngrids
          CALL grid_coords (ng, iadm)
        END DO
# endif
!
!-----------------------------------------------------------------------
!  Initialize time-stepping counter.
!-----------------------------------------------------------------------
!
        DO ng=1,ngrids
!!        ntstart(ng)=ntstart(ng)-1
          iic(ng)=ntstart(ng)
# ifdef JEDI
          jic(ng)=ntstart(ng)+1
          time4jedi(ng)=time(ng)+dt(ng)
# endif
          CALL time_string (time(ng), time_code(ng))
          ldefadj(ng)=.true.
          lwrtadj(ng)=.true.
        END DO
 
      END IF phase2
 
# ifdef PROFILE
!
!-----------------------------------------------------------------------
!  Turn off initialization time wall clock.
!-----------------------------------------------------------------------
!
      DO ng=1,ngrids
        DO thread=thread_range
          CALL wclock_off (ng, iadm, 2, __line__, myfile)
        END DO
      END DO
# endif
!
  10  FORMAT (/,1x,a,a,/,1x,'***********',/)
!
      RETURN
      END SUBROUTINE adm_initial
#endif
 
      END MODULE roms_kernel_mod