obs_sen_i4dvar_analysis.h

Go to the documentation of this file.

      MODULE roms_kernel_mod
!
!git $Id$
!=================================================== Andrew M. Moore ===
!  Copyright (c) 2002-2025 The ROMS Group           Hernan G. Arango   !
!    Licensed under a MIT/X style license                              !
!    See License_ROMS.md                                               !
!=======================================================================
!                                                                      !
!  ROMS I4D-VAR Observation Sensitivity Analysis Driver:               !
!                                                                      !
!  This driver evaluates the impact of each observation in the         !
!  4D-Var analysis increment by measuring their sensitivity over       !
!  a specified circulation functional index, J, similar to the         !
!  adjoint sensitivity driver. This is equivalent to taking the        !
!  adjoint of the 4D-Var algorithm.                                    !
!                                                                      !
!  Algorithm Outline:                                                  !
!                                                                      !
!                                                                      !
!    |----------------------------|-----------------------------|      !
!   t=t0                         t=t1                          t=t2    !
!         Assimilation window             Forecast period              !
!                                                                      !
!    |----------------------------------------------------------> NLM  !
!                                                                      !
!    <----------------------------------------------------------| ADM  !
!                                                                      !
!    |----------------------------> TLM                                !
!                                                                      !
!  (1) We begin by running an I4D-Var Lanczos calculation using k      !
!      inner-loops and 1 outer-loop for the period t=t0 to t1. We      !
!      will denote by xb(0) the background initial condition, and      !
!      the observations vector by y.  The resulting Lanczos vectors    !
!      that we save in the adjoint NetCDF file will be denoted by      !
!      q_i, where i=1,2,...,k.                                         !
!                                                                      !
!  (2) Next we run the NLM for the combined assimilation+forecast      !
!      period t=t0 to t2, where t2>t1. This represents the final       !
!      sweep of the NLM for the period t=t0 to t1 after exiting the    !
!      inner-loop in the I4D-Var plus the forecast period t=t1 to t2.  !
!      The initial condition for the NLM at t=t0 is xb(0) and not      !
!      the new estimated initial conditions. We save the basic state   !
!      trajectory, xb(t), of this NLM run for use in the adjoint       !
!      sensitivity calculation next, and for use in the TLM run        !
!      later in step (7).                                              !
!                                                                      !
!      Depending on time for which the sensitivity functional J(t)     !
!      is defined, this will dictate t2. For example, if J(t) is a     !
!      functional defined during the forecast interval t1<t<t2, then   !
!      t2>t1 for this run of the NLM. However, if J(t) is defined      !
!      during the assimilation interval t0<t<t1, then t2=t1. That is,  !
!      the definition of t2 should be flexible depending on the        !
!      choice of J.                                                    !
!                                                                      !
!  (3) The next step involves an adjoint sensitivity calculation       !
!      for the combined assimilation+forecast period t=t0 to t2. The   !
!      basic state trajectory for this calculation will be that from   !
!      the NLM run in step (2).                                        !
!                                                                      !
!  (4) After running the regular adjoint sensitivity calculation in    !
!      (3), we will have a full 3D-adjoint state vector at time t=t0.  !
!      Let's call this vector x(0). The next thing we want to do is    !
!      to compute the dot-product of x(0) with each of the Lanczos     !
!      vectors from the previous I4D-Var run. So if we ran I4D-Var     !
!      with k inner-loops we will have k Lanczos vectors which we      !
!      denote as q_i where i=1,2,...,k. So we will compute a_i=x'q_i   !
!      where x' is the transpose of the vector x(0), and a_i for       !
!      i=1,2,...,k are scalars, so there will be k of them.            !
!                                                                      !
!  (5) The next step is to invert the tridiagonal matrix associated    !
!      with the Lanczos vectors. Let's denote this matrix as T. So     !
!      what we want to solve T*b=a, where a is the k by 1 vector of    !
!      scalars a_i from step (4), and b is the k by 1 vector that we   !
!      want to find. So we solve for b by using a tridiagonal solver.  !
!                                                                      !
!  (6) The next step is to compute a weighted sum of the Lanczos       !
!      vectors. Let's call this z, where:                              !
!                                                                      !
!      z = SUM_i (b_i * q_i)    and i=1,2,...,k                        !
!                                                                      !
!      The b_i are obtained from solving the tridiagonal equation in   !
!      (5), and the q_i are the Lanczos vectors. The vector z is a     !
!      full-state vector and be used as an initial condition for the   !
!      TLM in step (8).                                                !
!                                                                      !
!  (7) Finally, we run the TLM from t=t0 to t=t1 using z from (6) as   !
!      the TLM initial condition. During this run of the TLM, we       !
!      need to read and process the observations that we used in       !
!      the I4D-Var of step (1) and write the TLM solution at the       !
!      observation points and times to the DAV(ng)%name NetCDF file.   !
!      The values that we write into this DAV(ng)%name are actually    !
!      the TLM values multiplied by error variance assigned to each    !
!      observation during the I4D-Var in step (1).                     !
!                                                                      !
!  These routines control the initialization,  time-stepping,  and     !
!  finalization of  ROMS  model following ESMF conventions:            !
!                                                                      !
!     ROMS_initialize                                                  !
!     ROMS_run                                                         !
!     ROMS_finalize                                                    !
!                                                                      !
!  References:                                                         !
!                                                                      !
!  Moore, A.M., H.G. Arango, G. Broquet, B.S. Powell, A.T. Weaver,     !
!    and J. Zavala-Garay, 2011: The Regional Ocean Modeling System     !
!    (ROMS)  4-dimensional variational data assimilations systems,     !
!    Part I - System overview and formulation, Prog. Oceanogr., 91,    !
!    34-49, doi:10.1016/j.pocean.2011.05.004.                          !
!                                                                      !
!  Moore, A.M., H.G. Arango, G. Broquet, C. Edward, M. Veneziani,      !
!    B. Powell, D. Foley, J.D. Doyle, D. Costa, and P. Robinson,       !
!    2011: The Regional Ocean Modeling System (ROMS) 4-dimensional     !
!    variational data assimilations systems, Part II - Performance     !
!    and application to the California Current System, Prog.           !
!    Oceanogr., 91, 50-73, doi:10.1016/j.pocean.2011.05.003.           !
!                                                                      !
!  Moore, A.M., H.G. Arango, G. Broquet, C. Edward, M. Veneziani,      !
!    B. Powell, D. Foley, J.D. Doyle, D. Costa, and P. Robinson,       !
!    2011: The Regional Ocean Modeling System (ROMS) 4-dimensional     !
!    variational data assimilations systems, Part III - Observation    !
!    impact and observation sensitivity in the California Current      !
!    System, Prog. Oceanogr., 91, 74-94,                               !
!    doi:10.1016/j.pocean.2011.05.005.                                 !
!                                                                      !
!=======================================================================
!
      USE mod_param
      USE mod_parallel
      USE mod_arrays
      USE mod_fourdvar
      USE mod_iounits
      USE mod_ncparam
      USE mod_netcdf
#if defined PIO_LIB && defined DISTRIBUTE
      USE mod_pio_netcdf
#endif
      USE mod_scalars
      USE mod_stepping
!
#ifdef ADJUST_BOUNDARY
      USE mod_boundary,       ONLY : initialize_boundary
#endif
#if defined OBS_IMPACT && defined OBS_IMPACT_SPLIT
      USE mod_forces,         ONLY : initialize_forces
      USE mod_ocean,          ONLY : initialize_ocean
#endif
!
#ifdef BALANCE_OPERATOR
      USE ad_balance_mod,     ONLY : ad_balance
#endif
      USE ad_convolution_mod, ONLY : ad_convolution
      USE ad_variability_mod, ONLY : ad_variability
      USE close_io_mod,       ONLY : close_inp, close_out
      USE def_mod_mod,        ONLY : def_mod
      USE get_state_mod,      ONLY : get_state
      USE inp_par_mod,        ONLY : inp_par
#ifdef MCT_LIB
# ifdef ATM_COUPLING
      USE mct_coupler_mod,    ONLY : initialize_ocn2atm_coupling
# endif
# ifdef WAV_COUPLING
      USE mct_coupler_mod,    ONLY : initialize_ocn2wav_coupling
# endif
#endif
      USE stats_modobs_mod,   ONLY : stats_modobs
      USE stdout_mod,         ONLY : set_stdoutunit, stdout_unit
      USE strings_mod,        ONLY : founderror
#ifdef BALANCE_OPERATOR
      USE tl_balance_mod,     ONLY : tl_balance
#endif
      USE tl_convolution_mod, ONLY : tl_convolution
      USE tl_variability_mod, ONLY : tl_variability
      USE wrt_rst_mod,        ONLY : wrt_rst
#if defined BALANCE_OPERATOR && defined ZETA_ELLIPTIC
      USE zeta_balance_mod,   ONLY : balance_ref, biconj
#endif
!
      implicit none
!
      PUBLIC :: roms_initialize
      PUBLIC :: roms_run
      PUBLIC :: roms_finalize
!
      CONTAINS
!
      SUBROUTINE roms_initialize (first, mpiCOMM)
!
!=======================================================================
!                                                                      !
!  This routine allocates and initializes ROMS state variables         !
!  and internal and external parameters.                               !
!                                                                      !
!=======================================================================
!
!  Imported variable declarations.
!
      logical, intent(inout) :: first
!
      integer, intent(in), optional :: mpiCOMM
!
!  Local variable declarations.
!
      logical :: allocate_vars = .true.
!
#ifdef DISTRIBUTE
      integer :: MyError, MySize
#endif
      integer :: NRMrec, STDrec, Tindex
      integer :: chunk_size, ng, thread
#ifdef _OPENMP
      integer :: my_threadnum
#endif
!
      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
!
!  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 _OPENMP
      mythread=my_threadnum()
#elif 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
 
#if defined MCT_LIB && (defined ATM_COUPLING || defined WAV_COUPLING)
!
!-----------------------------------------------------------------------
!  Initialize coupling streams between model(s).
!-----------------------------------------------------------------------
!
      DO ng=1,ngrids
# ifdef ATM_COUPLING
        CALL initialize_ocn2atm_coupling (ng, myrank)
# endif
# ifdef WAV_COUPLING
        CALL initialize_ocn2wav_coupling (ng, myrank)
# endif
      END DO
#endif
!
!-----------------------------------------------------------------------
!  Read in Lanczos algorithm coefficients (cg_beta, cg_delta) from
!  file LCZ(ng)%name NetCDF (I4D-Var adjoint file), as computed in the
!  I4D-Var Lanczos data assimilation algorithm for the first outer
!  loop.  They are needed here, in routine "ini_lanczos", to compute
!  the tangent linear model initial conditions as the weighted sum
!  of the Lanczos vectors. The weighting coefficient are computed
!  by solving a tri-diagonal system that uses cg_beta and cg_gamma.
!-----------------------------------------------------------------------
!
      sourcefile=myfile
      DO ng=1,ngrids
        SELECT CASE (lcz(ng)%IOtype)
          CASE (io_nf90)
            CALL netcdf_get_fvar (ng, iadm, lcz(ng)%name,               &
     &                            'cg_beta', cg_beta)
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
            CALL netcdf_get_fvar (ng, iadm, lcz(ng)%name,               &
     &                            'cg_delta', cg_delta)
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
#if defined PIO_LIB && defined DISTRIBUTE
          CASE (io_pio)
            CALL pio_netcdf_get_fvar (ng, iadm, lcz(ng)%name,           &
     &                                'cg_beta', cg_beta)
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
            CALL pio_netcdf_get_fvar (ng, iadm, lcz(ng)%name,           &
     &                                'cg_delta', cg_delta)
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
#endif
        END SELECT
      END DO
 
#ifdef SKIP_NLM
!
!-----------------------------------------------------------------------
!  If skiping runing nonlinear model, read in observation screening and
!  quality control flag.
!-----------------------------------------------------------------------
!
      wrtobsscale(1:ngrids)=.false.
      DO ng=1,ngrids
        SELECT CASE (lcz(ng)%IOtype)
          CASE (io_nf90)
            CALL netcdf_get_fvar (ng, itlm, lcz(ng)%name,               &
     &                            vname(1,idobss), obsscale)
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
# if defined PIO_LIB && defined DISTRIBUTE
          CASE (io_pio)
            CALL netcdf_get_fvar (ng, itlm, lcz(ng)%name,               &
     &                            vname(1,idobss), obsscale)
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
# endif
        END SELECT
      END DO
#endif
!
!-----------------------------------------------------------------------
!  Read in standard deviation factors for error covariance.
!-----------------------------------------------------------------------
!
!  Initial conditions standard deviation. They are loaded in Tindex=1
!  of the e_var(...,Tindex) state variables.
!
      stdrec=1
      tindex=1
      DO ng=1,ngrids
        CALL get_state (ng, 10, 10, std(1,ng), stdrec, tindex)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END DO
 
#ifdef ADJUST_BOUNDARY
!
!  Open boundary conditions standard deviation.
!
      stdrec=1
      tindex=1
      DO ng=1,ngrids
        CALL get_state (ng, 12, 12, std(3,ng), stdrec, tindex)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END DO
#endif
#if defined ADJUST_WSTRESS || defined ADJUST_STFLUX
!
!  Surface forcing standard deviation.
!
      stdrec=1
      tindex=1
      DO ng=1,ngrids
        CALL get_state (ng, 13, 13, std(4,ng), stdrec, tindex)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END DO
#endif
!
!-----------------------------------------------------------------------
!  Read in initial conditions, boundary conditions, and surface
!  forcing error covariance covariance normalization factors.
!-----------------------------------------------------------------------
!
      nrmrec=1
      DO ng=1,ngrids
        CALL get_state (ng, 14, 14, nrm(1,ng), nrmrec, 1)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
#ifdef ADJUST_BOUNDARY
        CALL get_state (ng, 16, 16, nrm(3,ng), nrmrec, 1)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
#endif
#if defined ADJUST_WSTRESS || defined ADJUST_STFLUX
        CALL get_state (ng, 17, 17, nrm(4,ng), nrmrec, 1)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
#endif
      END DO
!
      RETURN
      END SUBROUTINE roms_initialize
!
      SUBROUTINE roms_run (RunInterval)
!
!=======================================================================
!                                                                      !
!  This routine computes the adjoint sensitivity analysis, dJ/dS,      !
!  to the specified functional J.  The sensitivity masking arrays      !
!  Rscope, Uscope, and Vscope are used to evaluate the functional      !
!  in the desired spatial area.                                        !
!                                                                      !
!=======================================================================
!
!  Imported variable declarations
!
      real(dp), intent(in) :: RunInterval            ! seconds
!
!  Local variable declarations.
!
      logical :: Ladjoint, Lweak
!
      integer :: i, lstr, ng, tile
      integer :: Fcount, Lbck, Lini, Litl, Rec
!
      real (r8) :: str_day, end_day
!
      character (len=*), parameter :: MyFile =                          &
     &  __FILE__//", ROMS_run"
!
!=======================================================================
!  Run model for all nested grids, if any.
!=======================================================================
!
!  Initialize relevant parameters.
!
      lini=1          ! 4DVAR initial conditions record in INI
      lbck=2          ! First guess initial conditions record in INI
      litl=1          ! TLM initial conditions record
      lweak=.false.
      DO ng=1,ngrids
#if defined ADJUST_BOUNDARY || defined ADJUST_STFLUX || \
    defined adjust_wstress
        lfinp(ng)=1         ! forcing index for input
        lfout(ng)=1         ! forcing index for output history files
#endif
#ifdef ADJUST_BOUNDARY
        lbinp(ng)=1         ! boundary index for input
        lbout(ng)=1         ! boundary index for output history files
#endif
        lnew(ng)=1
      END DO
!
!  Initialize nonlinear model with the estimated initial conditions
!  from the I4D-Var Lanczos algorithm.  Notice that the LreadBLK and
!  LreadFWD switches are turned off to suppress processing of the
!  structures when "check_multifile" during
!  nonlinear model execution.
!
      DO ng=1,ngrids
#ifdef FORWARD_FLUXES
        lreadblk(ng)=.false.
#endif
        lreadfwd(ng)=.false.
        wrtnlmod(ng)=.false.
        wrttlmod(ng)=.false.
        rst(ng)%Rindex=0
        fcount=rst(ng)%load
        rst(ng)%Nrec(fcount)=0
      END DO
!
      CALL initial
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
#if defined BALANCE_OPERATOR && defined ZETA_ELLIPTIC
!
!  Compute the reference zeta and biconjugate gradient arrays
!  required for the balance of free surface.
!
      IF (balance(isfsur)) THEN
        DO ng=1,ngrids
          DO tile=first_tile(ng),last_tile(ng),+1
            CALL balance_ref (ng, tile, lini)
            CALL biconj (ng, tile, inlm, lini)
          END DO
          wrtzetaref(ng)=.true.
        END DO
      END IF
#endif
#ifndef SKIP_NLM
!
!  Run nonlinear model for the combined assimilation plus forecast
!  period, t=t0 to t2. Save nonlinear (basic state) tracjectory, xb(t),
!  needed by the adjoint model.
!
      DO ng=1,ngrids
        IF (master) THEN
          WRITE (stdout,10) 'NL', ng, ntstart(ng), ntend(ng)
        END IF
      END DO
!
# ifdef SOLVE3D
      CALL main3d (runinterval)
# else
      CALL main2d (runinterval)
# endif
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
#endif
!
!  Set structure for the nonlinear forward trajectory to be processed
!  by the tangent linear and adjoint models. Also, set switches to
!  process the FWD structure in routine "check_multifile". Notice that
!  it is possible to split solution into multiple NetCDF files to reduce
!  their size.
!
        CALL edit_multifile ('HIS2FWD')
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        DO ng=1,ngrids
          lreadfwd(ng)=.true.
        END DO
 
#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.
!
!  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
!
!  Initialize adjoint model and define sensitivity functional.
!
      lstiffness=.false.
      DO ng=1,ngrids
        CALL ad_initial (ng)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END DO
!
!  Activate adjoint output.
!
      DO ng=1,ngrids
        ldefadj(ng)=.true.
        lwrtadj(ng)=.true.
        lcycleadj(ng)=.false.
      END DO
!
!  Time-step adjoint model for the combined plus forecast period,
!  t=t2 to t0. Compute the gradient or index, dJ/dS, of the
!  sensitivity functional.
!
      DO ng=1,ngrids
        str_day=tdays(ng)
        end_day=str_day-ntimes(ng)*dt(ng)*sec2day
        IF ((dstrs(ng).eq.0.0_r8).and.(dends(ng).eq.0.0_r8)) THEN
          dstrs(ng)=end_day
          dends(ng)=str_day
        END IF
        IF (master) THEN
          WRITE (stdout,20) 'AD', ntstart(ng), ntend(ng),               &
     &                      dends(ng), dstrs(ng)
        END IF
#ifndef OBS_IMPACT
        IF ((dstrs(ng).gt.str_day).or.(dstrs(ng).lt.end_day)) THEN
          IF (master)  WRITE (stdout,30) 'DstrS = ', dstrs(ng),         &
     &                                   end_day, str_day
          exit_flag=7
          RETURN
        END IF
        IF ((dends(ng).gt.str_day).or.(dends(ng).lt.end_day)) THEN
          IF (master)  WRITE (stdout,30) 'DendS = ', dends(ng),         &
     &                                   end_day, str_day
          exit_flag=7
          RETURN
        END IF
#endif
      END DO
!
!   If DstrS=DendS=dstart, skip the adjoint model and read adjoint solution
!   from ADS netcdf file, record 1.
!
      ladjoint=.true.
      DO ng=1,ngrids
        IF ((dstrs(ng).eq.dends(ng)).and.(dstrs(ng).eq.dstart)) THEN
          rec=1
          CALL get_state (ng, iadm, 4, ads(ng), rec, lnew(ng))
          ladjoint=.false.
        END IF
      END DO
 
      IF (ladjoint) THEN
        DO ng=1,ngrids
          IF (master) THEN
            WRITE (stdout,10) 'AD', ng, ntstart(ng), ntend(ng)
          END IF
        END DO
!
#ifdef SOLVE3D
        CALL ad_main3d (runinterval)
#else
        CALL ad_main2d (runinterval)
#endif
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END IF
!
!  Load full adjoint sensitivity vector, x(0), for t=t0 into adjoint
!  state arrays at index Lnew.
!
      IF (ladjoint) THEN
        DO ng=1,ngrids
          CALL get_state (ng, iadm, 4, adm(ng), adm(ng)%Rindex,         &
     &                    lnew(ng))
        END DO
      END IF
 
#ifdef BALANCE_OPERATOR
!
!  Read in NLM reference state in readiness for the balance operator.
!
      DO ng=1,ngrids
        CALL get_state (ng, inlm, 2, ini(ng), lini, lini)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        nrhs(ng)=lini
      END DO
#endif
!
!  Convert adjoint solution to v-space since the Lanczos vectors
!  are in v-space.
!
      DO ng=1,ngrids
        DO tile=first_tile(ng),last_tile(ng),+1
#ifdef BALANCE_OPERATOR
          CALL ad_balance (ng, tile, lini, lnew(ng))
#endif
          CALL ad_variability (ng, tile, lnew(ng), lweak)
          CALL ad_convolution (ng, tile, lnew(ng), lweak, 2)
        END DO
      END DO
!
!  Check Lanczos vector input file and determine t=t1. That is, the
!  time to run the tangent linear model.  This time must be the same
!  as the I4D-Var Lanczos algorithm.
!
      sourcefile=myfile
      DO ng=1,ngrids
        SELECT CASE (lcz(ng)%IOtype)
          CASE (io_nf90)
            CALL netcdf_get_ivar (ng, iadm, lcz(ng)%name,               &
     &                            'ntimes', ntimes(ng))
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
#if defined PIO_LIB && defined DISTRIBUTE
          CASE (io_pio)
            CALL pio_netcdf_get_ivar (ng, iadm, lcz(ng)%name,           &
     &                                'ntimes', ntimes(ng))
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
#endif
        END SELECT
      END DO
 
#ifndef OBS_IMPACT
!
!  Initialize nonlinear model with the same initial conditions, xb(0),
!  Lbck record in INI(ng)%name. This is the first guess NLM initial
!  conditions used to start the I4D-Var Lanczos algorithm. Notice that
!  the LreadBLK and LreadFWD switches are turned off to suppress
!  processing of the structures when "check_multifile" during
!  nonlinear model execution.
!
      DO ng=1,ngrids
        ldefini(ng)=.false.
# ifdef FORWARD_FLUXES
        lreadblk(ng)=.false.
# endif
        lreadfwd(ng)=.false.
        wrtnlmod(ng)=.false.
        wrttlmod(ng)=.false.
        rst(ng)%Rindex=0
        ini(ng)%Rindex=lbck
        fcount=rst(ng)%load
        rst(ng)%Nrec(fcount)=0
      END DO
!
      CALL initial
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Run nonlinear model for the combined assimilation plus forecast
!  period, t=t0 to t2. Save nonlinear (basic state) tracjectory, xb(t),
!  needed by the tangent linear model.
!
      DO ng=1,ngrids
        IF (master) THEN
          WRITE (stdout,10) 'NL', ng, ntstart(ng), ntend(ng)
        END IF
      END DO
!
# ifdef SOLVE3D
      CALL main3d (runinterval)
# else
      CALL main2d (runinterval)
# endif
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Set structure for the nonlinear forward trajectory to be processed
!  by the tangent linear and adjoint models. Also, set switches to
!  process the FWD structure in routine "check_multifile".  Notice that
!  it is possible to split solution into multiple NetCDF files to reduce
!  their size.
!
      CALL edit_multifile ('HIS2FWD')
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      DO ng=1,ngrids
        lreadfwd(ng)=.true.
      END DO
 
# 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.
!
      CALL edit_multifile ('QCK2BLK')
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      DO ng=1,ngrids
        lreadblk(ng)=.true.
      END DO
# endif
#endif
!
!  Initialize tangent linear model with the weighted sum of the
!  Lanczos vectors, steps (4) to (6) from the algorithm summary
!  above.
!
      DO ng=1,ngrids
        CALL tl_initial (ng)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        ldeftlm(ng)=.true.
        lwrttlm(ng)=.true.
        wrttlmod(ng)=.true.
      END DO
!
!  Convert TL initial condition from v-space to x-space.
!
      DO ng=1,ngrids
        DO tile=first_tile(ng),last_tile(ng),+1
          CALL tl_convolution (ng, tile, litl, lweak, 2)
          CALL tl_variability (ng, tile, litl, lweak)
#ifdef BALANCE_OPERATOR
          CALL tl_balance (ng, tile, lini, litl)
#endif
        END DO
      END DO
!
!  Define output 4DVAR NetCDF file containing the sensitivity at the
!  observation locations.
!
      DO ng=1,ngrids
        ldefmod(ng)=.true.
        CALL def_mod (ng)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        wrtimpact_tot(ng)=.true.
#ifdef OBS_IMPACT_SPLIT
        wrtimpact_ic(ng)=.false.
#endif
      END DO
!
!  Write out outer loop beeing processed.
!
      sourcefile=myfile
      DO ng=1,ngrids
        SELECT CASE (dav(ng)%IOtype)
          CASE (io_nf90)
            CALL netcdf_put_ivar (ng, inlm, dav(ng)%name,               &
     &                            'Nimpact', nimpact,                   &
     &                            (/0/), (/0/),                         &
     &                            ncid = dav(ng)%ncid)
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
#if defined PIO_LIB && defined DISTRIBUTE
          CASE (io_pio)
            CALL pio_netcdf_put_ivar (ng, inlm, dav(ng)%name,           &
     &                                'Nimpact', nimpact,               &
     &                                (/0/), (/0/),                     &
     &                                piofile = dav(ng)%pioFile)
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
#endif
        END SELECT
      END DO
!
!  Run tangent linear model for the assimilation period, t=t0 to t1.
!  Read and process the 4DVAR observations.
!
      DO ng=1,ngrids
        IF (master) THEN
          WRITE (stdout,10) 'TL', ng, ntstart(ng), ntend(ng)
        END IF
      END DO
!
#ifdef SOLVE3D
      CALL tl_main3d (runinterval)
#else
      CALL tl_main2d (runinterval)
#endif
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
#if defined OBS_IMPACT && defined OBS_IMPACT_SPLIT
!
!:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
!  Integrate tangent linear model with initial condition increments
!  only to compute the observation impact associated with the initial
!  conditions.
!:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
!
!  Load full adjoint sensitivity vector, x(0), for t=t0 into adjoint
!  state arrays at index Lnew.
!
      DO ng=1,ngrids
        CALL get_state (ng, iadm, 4, adm(ng), adm(ng)%Rindex, lnew(ng))
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END DO
 
# ifdef BALANCE_OPERATOR
!
!  Read in NLM reference state in readiness for the balance operator.
!
      DO ng=1,ngrids
        CALL get_state (ng, inlm, 2, ini(ng), lini, lini)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        nrhs(ng)=lini
      END DO
# endif
!
!  Clear the adjoint forcing and boundary condition increment arrays.
!
      DO ng=1,ngrids
        DO tile=first_tile(ng),last_tile(ng),+1
          CALL initialize_forces (ng, tile, iadm)
# ifdef ADJUST_BOUNDARY
          CALL initialize_boundary (ng, tile, iadm)
# endif
        END DO
      END DO
!
!  Convert adjoint solution to v-space since the Lanczos vectors
!  are in v-space.
!
      DO ng=1,ngrids
        DO tile=first_tile(ng),last_tile(ng),+1
# ifdef BALANCE_OPERATOR
          CALL ad_balance (ng, tile, lini, lnew(ng))
# endif
          CALL ad_variability (ng, tile, lnew(ng), lweak)
          CALL ad_convolution (ng, tile, lnew(ng), lweak, 2)
        END DO
      END DO
!
!  Check Lanczos vector input file and determine t=t1. That is, the
!  time to run the tangent linear model.  This time must be the same
!  as the I4D-Var Lanczos algorithm.
!
      sourcefile=myfile
      DO ng=1,ngrids
        SELECT CASE (lcz(ng)%IOtype)
          CASE (io_nf90)
            CALL netcdf_get_ivar (ng, iadm, lcz(ng)%name,               &
     &                            'ntimes', ntimes(ng))
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
# if defined PIO_LIB && defined DISTRIBUTE
          CASE (io_pio)
            CALL pio_netcdf_get_ivar (ng, iadm, lcz(ng)%name,           &
     &                                'ntimes', ntimes(ng))
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
# endif
        END SELECT
      END DO
!
!  Initialize tangent linear model with the weighted sum of the
!  Lanczos vectors, steps (4) to (6) from the algorithm summary
!  above.
!
      DO ng=1,ngrids
        CALL tl_initial (ng)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        ldeftlm(ng)=.true.
        lwrttlm(ng)=.true.
        wrttlmod(ng)=.true.
        wrtimpact_tot(ng)=.false.
        wrtimpact_ic(ng)=.true.
# if defined ADJUST_WSTRESS || defined ADJUST_STFLUX
        wrtimpact_fc(ng)=.false.
# endif
# if defined ADJUST_BOUNDARY
        wrtimpact_bc(ng)=.false.
# endif
      END DO
!
!  Convert TL initial condition from v-space to x-space.
!
      DO ng=1,ngrids
        DO tile=first_tile(ng),last_tile(ng),+1
          CALL tl_convolution (ng, tile, litl, lweak, 2)
          CALL tl_variability (ng, tile, litl, lweak)
# ifdef BALANCE_OPERATOR
          CALL tl_balance (ng, tile, lini, litl)
# endif
        END DO
      END DO
!
!  Run tangent linear model for the assimilation period, t=t0 to t1.
!  Read and process the 4DVAR observations.
!
      DO ng=1,ngrids
        IF (master) THEN
          WRITE (stdout,10) 'TL', ng, ntstart(ng), ntend(ng)
        END IF
      END DO
!
# ifdef SOLVE3D
      CALL tl_main3d (runinterval)
# else
      CALL tl_main2d (runinterval)
# endif
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
# if defined ADJUST_WSTRESS || defined ADJUST_STFLUX
!
!:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
!  Integrate tangent linear model with surface forcing increments only
!  to compute the observation impact associated with the surface forcing.
!:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
!
!  Load full adjoint sensitivity vector, x(0), for t=t0 into adjoint
!  state arrays at index Lnew.
!
      DO ng=1,ngrids
        CALL get_state (ng, iadm, 4, adm(ng), adm(ng)%Rindex, lnew(ng))
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END DO
 
#  ifdef BALANCE_OPERATOR
!
!  Read in NLM reference state in readiness for the balance operator.
!
      DO ng=1,ngrids
        CALL get_state (ng, inlm, 2, ini(ng), lini, lini)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        nrhs(ng)=lini
      END DO
#  endif
!
!  Clear the adjoint initial condition and boundary condition increment
!  arrays.
!
      DO ng=1,ngrids
        DO tile=first_tile(ng),last_tile(ng),+1
          CALL initialize_ocean (ng, tile, iadm)
#  ifdef ADJUST_BOUNDARY
          CALL initialize_boundary (ng, tile, iadm)
#  endif
        END DO
      END DO
!
!  Convert adjoint solution to v-space since the Lanczos vectors
!  are in v-space.
!
      DO ng=1,ngrids
        DO tile=first_tile(ng),last_tile(ng),+1
#  ifdef BALANCE_OPERATOR
          CALL ad_balance (ng, tile, lini, lnew(ng))
#  endif
          CALL ad_variability (ng, tile, lnew(ng), lweak)
          CALL ad_convolution (ng, tile, lnew(ng), lweak, 2)
        END DO
      END DO
!
!  Check Lanczos vector input file and determine t=t1. That is, the
!  time to run the tangent linear model.  This time must be the same
!  as the I4D-Var Lanczos algorithm.
!
      sourcefile=myfile
      DO ng=1,ngrids
        SELECT CASE (lcz(ng)%IOtype)
          CASE (io_nf90)
            CALL netcdf_get_ivar (ng, iadm, lcz(ng)%name,               &
     &                            'ntimes', ntimes(ng))
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
#  if defined PIO_LIB && defined DISTRIBUTE
          CASE (io_pio)
            CALL pio_netcdf_get_ivar (ng, iadm, lcz(ng)%name,           &
     &                                'ntimes', ntimes(ng))
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
#  endif
        END SELECT
      END DO
!
!  Initialize tangent linear model with the weighted sum of the
!  Lanczos vectors, steps (4) to (6) from the algorithm summary
!  above.
!
      DO ng=1,ngrids
        CALL tl_initial (ng)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        ldeftlm(ng)=.true.
        lwrttlm(ng)=.true.
        wrttlmod(ng)=.true.
        wrtimpact_tot(ng)=.false.
        wrtimpact_ic(ng)=.false.
        wrtimpact_fc(ng)=.true.
#  if defined ADJUST_BOUNDARY
        wrtimpact_bc(ng)=.false.
#  endif
      END DO
!
!  Convert TL initial condition from v-space to x-space.
!
      DO ng=1,ngrids
        DO tile=first_tile(ng),last_tile(ng),+1
          CALL tl_convolution (ng, tile, litl, lweak, 2)
          CALL tl_variability (ng, tile, litl, lweak)
#  ifdef BALANCE_OPERATOR
          CALL tl_balance (ng, tile, lini, litl)
#  endif
        END DO
      END DO
!
!  Run tangent linear model for the assimilation period, t=t0 to t1.
!  Read and process the 4DVAR observations.
!
      DO ng=1,ngrids
        IF (master) THEN
          WRITE (stdout,10) 'TL', ng, ntstart(ng), ntend(ng)
        END IF
      END DO
!
#  ifdef SOLVE3D
      CALL tl_main3d (runinterval)
#  else
      CALL tl_main2d (runinterval)
#  endif
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
# endif
 
# if defined ADJUST_BOUNDARY
!
!:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
!  Integrate tangent linear model with boundary increments only
!  to compute the observation impact associated with the boundaries.
!:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
!
!  Load full adjoint sensitivity vector, x(0), for t=t0 into adjoint
!  state arrays at index Lnew.
!
      DO ng=1,ngrids
        CALL get_state (ng, iadm, 4, adm(ng), adm(ng)%Rindex, lnew(ng))
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END DO
 
#  ifdef BALANCE_OPERATOR
!
!  Read in NLM reference state in readiness for the balance operator.
!
      DO ng=1,ngrids
        CALL get_state (ng, inlm, 2, ini(ng), lini, lini)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        nrhs(ng)=lini
      END DO
#  endif
!
!  Clear the adjoint increment initial condition and forcing arrays.
!
      DO ng=1,ngrids
        DO tile=first_tile(ng),last_tile(ng),+1
          CALL initialize_ocean (ng, tile, iadm)
          CALL initialize_forces (ng, tile, iadm)
        END DO
      END DO
!
!  Convert adjoint solution to v-space since the Lanczos vectors
!  are in v-space.
!
      DO ng=1,ngrids
        DO tile=first_tile(ng),last_tile(ng),+1
#  ifdef BALANCE_OPERATOR
          CALL ad_balance (ng, tile, lini, lnew(ng))
#  endif
          CALL ad_variability (ng, tile, lnew(ng), lweak)
          CALL ad_convolution (ng, tile, lnew(ng), lweak, 2)
        END DO
      END DO
!
!  Check Lanczos vector input file and determine t=t1. That is, the
!  time to run the tangent linear model.  This time must be the same
!  as the I4D-Var Lanczos algorithm.
!
      sourcefile=myfile
      DO ng=1,ngrids
        SELECT CASE (lcz(ng)%IOtype)
          CASE (io_nf90)
            CALL netcdf_get_ivar (ng, iadm, lcz(ng)%name,               &
     &                            'ntimes', ntimes(ng))
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
#  if defined PIO_LIB && defined DISTRIBUTE
          CASE (io_pio)
            CALL pio_netcdf_get_ivar (ng, iadm, lcz(ng)%name,           &
     &                                'ntimes', ntimes(ng))
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
#  endif
        END SELECT
      END DO
!
!  Initialize tangent linear model with the weighted sum of the
!  Lanczos vectors, steps (4) to (6) from the algorithm summary
!  above.
!
      DO ng=1,ngrids
        CALL tl_initial (ng)
        IF (founderror(exit_flag, noerror, __line__,  myfile)) RETURN
 
        ldeftlm(ng)=.true.
        lwrttlm(ng)=.true.
        wrttlmod(ng)=.true.
        wrtimpact_tot(ng)=.false.
        wrtimpact_ic(ng)=.false.
#  if defined ADJUST_WSTRESS || defined ADJUST_STFLUX
        wrtimpact_fc(ng)=.false.
#  endif
        wrtimpact_bc(ng)=.true.
      END DO
!
!  Convert TL initial condition from v-space to x-space.
!
      DO ng=1,ngrids
        DO tile=first_tile(ng),last_tile(ng),+1
          CALL tl_convolution (ng, tile, litl, lweak, 2)
          CALL tl_variability (ng, tile, litl, lweak)
#  ifdef BALANCE_OPERATOR
          CALL tl_balance (ng, tile, lini, litl)
#  endif
        END DO
      END DO
!
!  Run tangent linear model for the assimilation period, t=t0 to t1.
!  Read and process the 4DVAR observations.
!
      DO ng=1,ngrids
        IF (master) THEN
          WRITE (stdout,10) 'TL', ng, ntstart(ng), ntend(ng)
        END IF
      END DO
!
#  ifdef SOLVE3D
      CALL tl_main3d (runinterval)
#  else
      CALL tl_main2d (runinterval)
#  endif
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
# endif
 
#endif
!
 10   FORMAT (/,1x,a,1x,'ROMS: started time-stepping:',                 &
     &        ' (Grid: ',i2.2,' TimeSteps: ',i8.8,' - ',i8.8,')',/)
 20   FORMAT (/,1x,a,1x,'ROMS: started time-stepping:',                 &
     &        '( TimeSteps: ',i8.8,' - ',i8.8,')',/,15x,                &
     &        'adjoint forcing time range: ',f12.4,' - ',f12.4 ,/)
 30   FORMAT (/,' Out of range adjoint forcing time, ',a,f12.4,/,       &
     &        ' It must be between ',f12.4,' and ',f12.4)
!
      RETURN
      END SUBROUTINE roms_run
!
      SUBROUTINE roms_finalize
!
!=======================================================================
!                                                                      !
!  This routine terminates ROMS nonlinear and adjoint models           !
!  execution.                                                          !
!                                                                      !
!=======================================================================
!
!  Local variable declarations.
!
      integer :: Fcount, ng, thread
!
      character (len=*), parameter :: MyFile =                          &
     &  __FILE__//", ROMS_finalize"
!
!-----------------------------------------------------------------------
!  Read and write observation variables for completeness.
!-----------------------------------------------------------------------
!
      DO ng=1,ngrids
#ifdef DISTRIBUTE
        CALL stats_modobs (ng, myrank)
#else
        CALL stats_modobs (ng, -1)
#endif
      END DO
!
!-----------------------------------------------------------------------
!  If blowing-up, save latest model state into RESTART NetCDF file.
!-----------------------------------------------------------------------
!
!  If cycling restart records, write solution into the next record.
!
      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
 
      END MODULE roms_kernel_mod