tl_initial.F

Go to the documentation of this file.

#include "cppdefs.h"
#ifdef TANGENT
      SUBROUTINE tl_initial (ng)
!
!git $Id$
!================================================== Hernan G. Arango ===
!  Copyright (c) 2002-2025 The ROMS Group                              !
!    Licensed under a MIT/X style license                              !
!    See License_ROMS.md                                               !
!=======================================================================
!                                                                      !
!  This routine initializes all tangent linear model variables.        !
!                                                                      !
!=======================================================================
!
      USE mod_param
      USE mod_parallel
# ifdef BBL_MODEL_NOT_YET
      USE mod_bbl
# endif
# ifdef FOUR_DVAR
#  ifdef SOLVE3D
      USE mod_coupling
#  endif
#  ifdef SP4DVAR
      USE mod_forces
#  endif
      USE mod_fourdvar
# endif
      USE mod_grid
      USE mod_iounits
      USE mod_ncparam
      USE mod_ocean
      USE mod_scalars
      USE mod_stepping
!
# ifdef ADJUST_BOUNDARY
      USE mod_boundary,        ONLY : initialize_boundary
# endif
# if defined ADJUST_STFLUX || defined ADJUST_WSTRESS
      USE mod_forces,          ONLY : initialize_forces
# endif
!
      USE analytical_mod
      USE close_io_mod,        ONLY : close_inp
      USE dateclock_mod,       ONLY : time_string
# ifdef DISTRIBUTE
      USE distribute_mod,      ONLY : mp_bcasti
# endif
      USE get_state_mod,       ONLY : get_state
# ifdef WET_DRY
      USE get_wetdry_mod,      ONLY : get_wetdry
# endif
# ifdef TLM_CHECK
      USE ini_adjust_mod,      ONLY : tl_ini_perturb
# endif
# if defined SOLVE3D && defined NONLINEAR
      USE ini_fields_mod,      ONLY : set_zeta_timeavg
# endif
      USE ini_hmixcoef_mod,    ONLY : ini_hmixcoef
# ifdef I4DVAR_ANA_SENSITIVITY
      USE ini_lanczos_mod,     ONLY : ini_lanczos
# endif
# if defined WAV_COUPLING_NOT_YET && defined MCT_LIB
      USE mct_coupler_mod,     ONLY : ocn2wav_coupling
# endif
# ifdef OBSERVATIONS
      USE obs_initial_mod,     ONLY : obs_initial
# endif
# ifdef SOLVE3D
      USE set_depth_mod,       ONLY : set_depth
# endif
# ifdef MASKING
      USE set_masks_mod,       ONLY : set_masks
# endif
      USE stiffness_mod,       ONLY : stiffness
      USE strings_mod,         ONLY : founderror
# if defined RBL4DVAR           || defined R4DVAR      || \
     defined sensitivity_4dvar  || defined tl_rbl4dvar || \
     defined tl_r4dvar
      USE tl_set_depth_mod,    ONLY : tl_bath
# endif
# ifdef FOUR_DVAR
      USE tl_def_ini_mod,      ONLY : tl_def_ini
# endif
# ifdef SOLVE3D
      USE tl_omega_mod,        ONLY : tl_omega
      USE tl_rho_eos_mod,      ONLY : tl_rho_eos
      USE tl_set_depth_mod,    ONLY : tl_set_depth
      USE tl_set_massflux_mod, ONLY : tl_set_massflux
# endif
# ifdef FOUR_DVAR
      USE tl_wrt_ini_mod,      ONLY : tl_wrt_ini
# endif
# ifdef WET_DRY
      USE wetdry_mod,          ONLY : wetdry
# endif
# if defined PROPAGATOR || \
    (defined masking    && (defined read_water || defined write_water))
      USE wpoints_mod,         ONLY : wpoints
# endif
!
      implicit none
!
!  Imported variable declarations.
!
      integer, intent(in) :: ng
!
!  Local variable declarations.
!
      logical :: update = .false.
!
      integer :: Fcount, IniRec, Tindex, wrtRec
      integer :: thread, tile
!
      real(dp) :: my_dstart
!
      character (len=*), parameter :: MyFile =                          &
     &  __FILE__
!
!=======================================================================
!   Initialize model variables.
!=======================================================================
!
      IF (master) THEN
# if defined PERTURBATION
        WRITE (stdout,10) nrun
 10     FORMAT (/,' <<<< Ensemble/Perturbation Run: ',i5.5,' >>>>',/)
# elif defined I4DVAR      || defined RBL4DVAR          || \
       defined r4dvar      || defined sensitivity_4dvar || \
       defined tl_rbl4dvar || defined tl_r4dvar
        WRITE (stdout,10) outer, inner
 10     FORMAT (/,' <<<< 4D Variational Data Assimilation, ',           &
     &          'Outer = ',i3.3, ', Inner = ',i3.3,' >>>>',/)
# endif
        WRITE (stdout,20) 'TL_INITIAL: Configuring and ',               &
     &                    'initializing tangent linear model ...'
 20     FORMAT (/,1x,a,a,/)
      END IF
!
!-----------------------------------------------------------------------
!  Initialize time stepping indices and counters.
!-----------------------------------------------------------------------
!
      iif(ng)=1
      indx1(ng)=1
      next_kstp(ng)=1
      kstp(ng)=1
      krhs(ng)=1
      knew(ng)=1
      predictor_2d_step(ng)=.false.
!
      iic(ng)=0
      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
      ntstart(ng)=int((time(ng)-dstart*day2sec)/dt(ng))+1
      ntend(ng)=ntstart(ng)+ntimes(ng)-1
      ntfirst(ng)=ntstart(ng)
      step_counter(ng)=0
!
      inirec=nrrec(ng)
      tindex=1
!
!  Initialize global diagnostics variables.
!
      avgke=0.0_dp
      avgpe=0.0_dp
      avgkp=0.0_dp
      volume=0.0_dp
 
# ifdef PROFILE
!
!-----------------------------------------------------------------------
!  Start time wall clocks.
!-----------------------------------------------------------------------
!
      DO thread=thread_range
        CALL wclock_on (ng, itlm, 2, __line__, myfile)
      END DO
# endif
 
# if defined OPT_OBSERVATIONS
!
!-----------------------------------------------------------------------
!  Initialize.
!-----------------------------------------------------------------------
!
!  Set initial conditions time record to read.
!
       inirec=1
 
# elif defined FOUR_DVAR    && \
     !(defined HESSIAN_SV   || defined HESSIAN_SO || \
       defined hessian_fsv) || defined tlm_check
!
!-----------------------------------------------------------------------
!  If variational data assimilation, reset several IO switches and
!  variables.
!-----------------------------------------------------------------------
 
#  ifndef I4DVAR_ANA_SENSITIVITY
#   ifdef I4DVAR
!
!  Set switch to create (TRUE) tangent linear initial conditions and
!  history NetCDF files or append (FALSE) to existing files. Then,
!  create tangent linear model initialization file and write zero
!  initial conditions for records 1 and 2.
!
      IF ((nrun.eq.erstr).and.(inner.eq.0)) THEN
        ldefitl(ng)=.true.
        CALL tl_def_ini (ng)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
#    ifdef DISTRIBUTE
        CALL tl_wrt_ini (ng, myrank, tindex, 1)
#    else
        CALL tl_wrt_ini (ng, -1, tindex, 1)
#    endif
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
#    ifdef DISTRIBUTE
        CALL tl_wrt_ini (ng, myrank, tindex, 2)
#    else
        CALL tl_wrt_ini (ng, -1, tindex, 2)
#    endif
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END IF
#   endif
#   ifndef WEAK_CONSTRAINT
!
!  Set switch to create tangent linear history file.
!
      IF (nrun.eq.erstr) THEN
        ldeftlm(ng)=.true.
      END IF
#   endif
!
!  Set record to read from initial tangent linear NetCDF file.
!
      inirec=itl(ng)%Rindex
 
#   ifdef ADJUST_BOUNDARY
!
!  Initialize open boundary counter for storage arrays.
!
      obccount(ng)=0
#   endif
#   if defined ADJUST_STFLUX || defined ADJUST_WSTRESS
!
!  Initialize surface forcing counter for storage arrays.
!
      sfcount(ng)=0
#   endif
#   if !defined WEAK_CONSTRAINT
!
!  Reset tangent linear model history time record counters. These
!  counters are reset in every iteration pass. This file is created
!  on the first iteration pass.
!
      tlm(ng)%Rindex=0
      fcount=tlm(ng)%Fcount
      tlm(ng)%Nrec(fcount)=0
      lwrttlm(ng)=.true.
#   endif
#  endif
# endif
!
!-----------------------------------------------------------------------
!  Set application grid, metrics, and associated variables and
!  parameters.
!-----------------------------------------------------------------------
!
      IF (setgridconfig(ng)) THEN
        CALL set_grid (ng, itlm)
        setgridconfig(ng)=.false.
      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 tile=first_tile(ng),last_tile(ng),+1
        CALL ini_hmixcoef (ng, tile, itlm)
      END DO
 
# ifdef ANA_SPONGE
!
!-----------------------------------------------------------------------
!  Increase horizontal mixing coefficients in sponge areas using
!  analytical functions.
!-----------------------------------------------------------------------
!
      IF (lsponge(ng)) THEN
        DO tile=first_tile(ng),last_tile(ng),+1
          CALL ana_sponge (ng, tile, itlm)
        END DO
      END IF
# endif
!
!=======================================================================
!  Initialize model state variables and forcing.  This part is
!  executed for each ensemble/perturbation/iteration pass.
!=======================================================================
 
# if defined FOUR_DVAR              && \
    !defined I4DVAR_ANA_SENSITIVITY && \
    !(defined HESSIAN_SV            || \
      defined hessian_so            || \
      defined hessian_fsv)
#  if defined OPT_OBSERVATIONS || defined TLM_CHECK || \
      defined weak_constraint
!
!  Clear tangent linear state variables.
!
      DO tile=first_tile(ng),last_tile(ng),+1
        CALL initialize_ocean (ng, tile, itlm)
#   ifdef SOLVE3D
        CALL initialize_coupling (ng, tile, 0)
#   endif
#   ifdef SP4DVAR
        CALL initialize_ocean (ng, tile, inlm)
        CALL initialize_forces (ng, tile, inlm)
#   endif
      END DO
 
#  else
#   ifndef WEAK_CONSTRAINT
!
!-----------------------------------------------------------------------
!  If first interation of the inner loop, clear all tangent linear
!  variables. In incrementatal 4DVAR, the tangent linear model is
!  started from rest on the first pass of the inner loop for each
!  outer loop iteration.
!-----------------------------------------------------------------------
!
      IF (inner.eq.0) THEN
        DO tile=first_tile(ng),last_tile(ng),+1
#    ifdef ADJUST_BOUNDARY
          CALL initialize_boundary (ng, tile, itlm)
#    endif
#    if defined ADJUST_STFLUX || defined ADJUST_WSTRESS
          CALL initialize_forces (ng, tile, itlm)
#    endif
          CALL initialize_ocean (ng, tile, itlm)
        END DO
!
!  Rewrite tangent linear initial NetCDF (record 1) with zero initial
!  conditions since the model needs to be started from at the first
!  pass of the inner loop.
!
        IF (nrun.gt.1) THEN
          wrtrec=1
#    ifdef DISTRIBUTE
          CALL tl_wrt_ini (ng, myrank, tindex, wrtrec)
#    else
          CALL tl_wrt_ini (ng, -1, tindex, wrtrec)
#    endif
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END IF
      END IF
#   endif
#  endif
# endif
 
# if defined RBL4DVAR          || defined R4DVAR      || \
     defined sensitivity_4dvar || defined tl_rbl4dvar || \
     defined tl_r4dvar
!
!-----------------------------------------------------------------------
!  Initialize tangent linear bathymetry to zero.
!-----------------------------------------------------------------------
!
      DO tile=first_tile(ng),last_tile(ng),+1
        CALL tl_bath (ng, tile)
      END DO
# endif
!
!-----------------------------------------------------------------------
!  Set tangent linear model state variables initial conditions.
!-----------------------------------------------------------------------
 
# ifdef ANA_INITIAL
!
!  Analytical initial conditions for momentum and active tracers.
!
      IF (nrrec(ng).eq.0) THEN
        DO tile=first_tile(ng),last_tile(ng),+1
          CALL ana_initial (ng, tile, itlm)
        END DO
      END IF
# endif
 
# if defined ANA_PASSIVE && defined SOLVE3D
!
!  Analytical initial conditions for inert passive tracers.
!
      IF (nrrec(ng).eq.0) THEN
        DO tile=first_tile(ng),last_tile(ng),+1
          CALL ana_passive (ng, tile, itlm)
        END DO
      END IF
# endif
 
# if defined ANA_BIOLOGY && defined SOLVE3D
!
!  Analytical initial conditions for biology.
!
      IF (nrrec(ng).eq.0) THEN
        DO tile=first_tile(ng),last_tile(ng),+1
          CALL ana_biology (ng, tile, itlm)
        END DO
      END IF
# endif
 
# if defined ANA_SEDIMENT_NOT_YET && defined SOLVE3D
!
!  Analytical initial conditions for sediment.
!
      IF (nrrec(ng).eq.0) THEN
        DO tile=first_tile(ng),last_tile(ng),+1
          CALL ana_sediment (ng, tile, itlm)
        END DO
      END IF
# endif
 
# ifdef I4DVAR_ANA_SENSITIVITY
!
!  Initialize with the weighted sum of all Lanczos vectors computed
!  from the first outer loop of the I4D-VAR Lanczos algorithm.
!
      DO tile=first_tile(ng),last_tile(ng),+1
        CALL ini_lanczos (ng, tile, lnew(ng), tindex)
      END DO
 
# else
!
!  Read in initial conditions for initial or restart NetCDF file.
!
#  ifdef INI_FILE
      CALL get_state (ng, itlm, 1, itl(ng), inirec, tindex)
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      time(ng)=io_time                       ! needed for shared-memory
#  else
      IF (nrrec(ng).ne.0) THEN
        CALL get_state (ng, itlm, 1, itl(ng), inirec, tindex)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        time(ng)=io_time                     ! needed for shared-memory
      END IF
#  endif
# endif
 
# ifdef WET_DRY
!
!-----------------------------------------------------------------------
!  Process initial wet/dry masks.
!-----------------------------------------------------------------------
!
!  If restart, read in wet/dry masks.
!
      IF (nrrec(ng).ne.0) THEN
#  ifdef DISTRIBUTE
        CALL get_wetdry (ng, myrank, itlm, inirec(ng))
        CALL mp_bcasti (ng, itlm, exit_flag)
#  else
        CALL get_wetdry (ng, -1, itlm, inirec(ng))
#  endif
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      ELSE
        DO tile=first_tile(ng),last_tile(ng),+1
          CALL wetdry (ng, tile, tindex(ng), .true.)
        END DO
      END IF
# endif
 
# ifdef OBSERVATIONS
!
!-----------------------------------------------------------------------
!  Open observations NetCDF file and initialize various variables
!  needed for processing the tangent linear state solution at
!  observation locations. Need to be done after processing initial
!  conditions since the correct initial time is needed to determine
!  the first "ObsTime" to process.
!-----------------------------------------------------------------------
!
      CALL obs_initial (ng, itlm, .false.)
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
# endif
 
# if defined ANA_PERTURB && defined SANITY_CHECK
!
!-----------------------------------------------------------------------
!  Perturb tangent linear initial conditions with analitical
!  expressions.
!-----------------------------------------------------------------------
!
      DO tile=first_tile(ng),last_tile(ng),+1
        CALL ana_perturb (ng, tile, itlm)
      END DO
# endif
 
# ifdef TLM_CHECK
!
!-----------------------------------------------------------------------
!  Perturb tangent linear state variable according to the outer loop
!  iteration with the steepest descent direction of the gradient
!  (adjoint state).
!-----------------------------------------------------------------------
!
      DO tile=first_tile(ng),last_tile(ng),+1
        CALL tl_ini_perturb (ng, tile, lnew(ng), tindex)
      END DO
#  endif
 
# ifdef SOLVE3D
!!
!!----------------------------------------------------------------------
!!  Compute initial time-evolving depths.
!!----------------------------------------------------------------------
!!
!!    DO tile=first_tile(ng),last_tile(ng),+1
!!      CALL tl_set_depth (ng, tile, iTLM)
!!    END DO
!!
!!----------------------------------------------------------------------
!!  Compute initial horizontal mass fluxes, Hz*u/n and Hz*v/m.
!!----------------------------------------------------------------------
!!
!!    DO tile=first_tile(ng),last_tile(ng),+1
!!      CALL tl_set_massflux (ng, tile, iTLM)
!!    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 tile=first_tile(ng),last_tile(ng),+1
!!      CALL tl_omega (ng, tile, iTLM)
!!      CALL tl_rho_eos (ng, tile, iTLM)
!!    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.
!-----------------------------------------------------------------------
!
      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, itlm)
        END DO
      END IF
#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.
!-----------------------------------------------------------------------
!
      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
!
!-----------------------------------------------------------------------
!  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.
!-----------------------------------------------------------------------
!
      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
 
# if defined SOLVE3D && defined TLM_DRIVER
!
!-----------------------------------------------------------------------
!  Compute nonlinear initial time-evolving depths.
!-----------------------------------------------------------------------
!
      DO tile=first_tile(ng),last_tile(ng),+1
#  ifdef NONLINEAR
        CALL set_zeta_timeavg (ng, tile, itlm)
#  endif
        CALL set_depth (ng, tile, itlm)
      END DO
# endif
 
# ifdef MASKING
!
!-----------------------------------------------------------------------
!  Set internal I/O mask arrays.
!-----------------------------------------------------------------------
!
      DO tile=first_tile(ng),last_tile(ng),+1
        CALL set_masks (ng, tile, itlm)
      END DO
# endif
 
# if defined PROPAGATOR || \
    (defined masking    && (defined read_water || defined write_water ))
!
!-----------------------------------------------------------------------
!  Set variables associated with the processing water points and/or
!  size of packed state arrays.
!-----------------------------------------------------------------------
!
      DO tile=first_tile(ng),last_tile(ng),+1
        CALL wpoints (ng, tile, itlm)
      END DO
# endif
 
# ifdef WEAK_CONSTRAINT
!
!-----------------------------------------------------------------------
!  If available, read in first TLM impulse forcing and its application
!  time. In true weak constraint applications, the impulse records
!  after the initial are associated with the model error and are
!  processed with different statistics. If there is only one (initial)
!  impulse forcing available, the assimilation tis similar to strong
!  constraint but in observation space.
!-----------------------------------------------------------------------
!
      IF (nadj(ng).lt.ntimes(ng)) THEN
        inirec=1
        CALL get_state (ng, 7, 7, tlf(ng), inirec, 1)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END IF
# endif
 
# if defined ANA_DRAG && defined UV_DRAG_GRID
!
!-----------------------------------------------------------------------
!  Set analytical spatially varying bottom friction parameter.
!-----------------------------------------------------------------------
!
      IF (nrun.eq.erstr) THEN
        DO tile=first_tile(ng),last_tile(ng),+1
          CALL ana_drag (ng, tile, itlm)
        END DO
      END IF
# endif
!
!-----------------------------------------------------------------------
!  Compute grid stiffness.
!-----------------------------------------------------------------------
!
      IF (lstiffness) THEN
        lstiffness=.false.
        DO tile=first_tile(ng),last_tile(ng),+1
          CALL stiffness (ng, tile, itlm)
        END DO
      END IF
 
# if defined FLOATS_NOT_YET || defined STATIONS
!
!-----------------------------------------------------------------------
!  If applicable, convert initial locations to fractional grid
!  coordinates.
!-----------------------------------------------------------------------
!
      CALL grid_coords (ng, itlm)
# endif
 
# if defined WAV_COUPLING_NOT_YET && defined MCT_LIB
!
!-----------------------------------------------------------------------
!  Read in initial forcing from coupled wave model.
!-----------------------------------------------------------------------
!
      DO tile=first_tile(ng),last_tile(ng),+1
        CALL ocn2wav_coupling (ng, tile)
      END DO
# endif
!
!-----------------------------------------------------------------------
!  Initialize time-stepping counter.
!-----------------------------------------------------------------------
!
      iic(ng)=ntstart(ng)
      CALL time_string (time(ng), time_code(ng))
 
# ifdef PROFILE
!
!-----------------------------------------------------------------------
!  Turn off initiialization time wall clock.
!-----------------------------------------------------------------------
!
      DO thread=thread_range
        CALL wclock_off (ng, itlm, 2, __line__, myfile)
      END DO
# endif
!
      RETURN
      END SUBROUTINE tl_initial
#else
      SUBROUTINE tl_initial
      RETURN
      END SUBROUTINE tl_initial
#endif