roms_kernel_mod Module Reference

Functions/Subroutines
subroutine, public	roms_initialize (first, mpicomm)
subroutine, public	roms_run (runinterval)
subroutine, public	roms_finalize
subroutine, private	iram_error (info, icall, string)
subroutine, private	iram_error (info, string)
subroutine, public	roms_initializep1 (first, mpicomm, kernel)
subroutine, public	roms_initializep2 (kernel)
subroutine, public	roms_initializep3 (kernel)
subroutine, public	roms_run (runinterval, kernel)
subroutine, private	nlm_initial (phase)
subroutine, private	tlm_initial (phase)
subroutine, private	adm_initial (phase)
subroutine, public	roms_run (runinterval)

Function/Subroutine Documentation

adm_initial()

subroutine, private roms_kernel_mod::adm_initial ( integer, intent(in) phase )

private

Definition at line 1419 of file jedi_roms.h.

!
!=======================================================================
!                                                                      !
!  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

References ad_get_data(), ad_get_idata(), mod_scalars::ad_ubar_xs, analytical_mod::ana_drag(), analytical_mod::ana_sponge(), mod_scalars::avgke, mod_scalars::avgkp, mod_scalars::avgpe, check_multifile(), close_io_mod::close_inp(), mod_scalars::day2sec, mod_scalars::dstart, mod_scalars::dt, edit_multifile(), mod_scalars::erstr, mod_scalars::exit_flag, mod_parallel::first_tile, mod_scalars::first_time, strings_mod::founderror(), grid_coords(), mod_param::iadm, mod_scalars::iic, mod_scalars::iif, mod_scalars::indx1, ini_hmixcoef_mod::ini_hmixcoef(), mod_scalars::initime, mod_scalars::jic, mod_stepping::knew, mod_stepping::krhs, mod_stepping::kstp, mod_parallel::last_tile, mod_scalars::ldefadj, mod_scalars::lreadblk, mod_scalars::lreadfrc, mod_scalars::lreadqck, mod_scalars::lsponge, mod_scalars::lstiffness, mod_scalars::lwrtadj, mod_parallel::master, mod_stepping::nf, mod_stepping::nfm1, mod_stepping::nfm2, mod_stepping::nfm3, mod_stepping::nfp1, mod_param::ngrids, mod_stepping::nnew, mod_scalars::noerror, mod_stepping::nrhs, mod_scalars::nrrec, mod_scalars::nrun, mod_stepping::nstp, mod_scalars::ntend, mod_scalars::ntfirst, mod_scalars::ntimes, mod_scalars::ntstart, mod_scalars::predictor_2d_step, mod_iounits::qck, mod_scalars::sec2day, set_masks_mod::set_masks(), mod_iounits::stdout, stiffness_mod::stiffness(), mod_scalars::synchro_flag, mod_scalars::tdays, mod_scalars::time, mod_scalars::time4jedi, mod_scalars::time_code, dateclock_mod::time_string(), mod_scalars::volume, wclock_off(), wclock_on(), and wetdry_mod::wetdry().

Referenced by roms_initializep1(), roms_initializep2(), and roms_initializep3().

Here is the call graph for this function:

Here is the caller graph for this function:

iram_error() [1/2]

subroutine private roms_kernel_mod::iram_error ( integer, intent(in) info, integer, intent(in) icall, character (len=*), intent(out) string )

private

Definition at line 861 of file afte_roms.h.

!
!=======================================================================
!                                                                      !
!  This routine decodes internal error messages from the Implicit      !
!  Restarted Arnoldi Method (IRAM) for the computation of optimal      !
!  perturbation Ritz eigenfunctions.                                   !
!                                                                      !
!=======================================================================
!
!  imported variable declarations.
!
      integer, intent(in) :: info, icall
!
      character (len=*), intent(out) :: string
!
!-----------------------------------------------------------------------
!  Decode error message from IRAM.
!-----------------------------------------------------------------------
!
      IF (info.eq.0)  THEN
        string='Normal exit                                            '
      ELSE IF (info.eq.1) THEN
        IF (icall.eq.1) THEN
          string='Maximum number of iterations taken                   '
        ELSE
          string='Could not reorder Schur vectors                      '
        END IF
      ELSE IF (info.eq.3) THEN
        string='No shifts could be applied during an IRAM cycle        '
      ELSE IF (info.eq.-1) THEN
        string='Nstate must be positive                                '
      ELSE IF (info.eq.-2) THEN
        string='NEV must be positive                                   '
      ELSE IF (info.eq.-3) THEN
        string='NCV must be greater NEV and less than or equal Nstate  '
      ELSE IF (info.eq.-4) THEN
        string='Maximum number of iterations must be greater than zero '
      ELSE IF (info.eq.-5) THEN
        string='WHICH must be one of LM, SM, LA, SA or BE              '
      ELSE IF (info.eq.-6) THEN
        string='BMAT must be one of I or G                             '
      ELSE IF (info.eq.-7) THEN
        string='Length of private work array SworkL is not sufficient  '
      ELSE IF (info.eq.-8) THEN
        IF (icall.eq.1) THEN
          string='Error return from LAPACK eigenvalue calculation      '
        ELSE
          string='Error in DLAHQR in the Shurn vectors calculation     '
        END IF
      ELSE IF (info.eq.-9) THEN
        IF (icall.eq.1) THEN
          string='Starting vector is zero'
        ELSE
          string='Error in DTREVC in the eigenvectors calculation      '
        END IF
      ELSE IF (info.eq.-10) THEN
        string='IPARAM(7) must be 1, 2, 3, 4, 5                        '
      ELSE IF (info.eq.-11) THEN
        string='IPARAM(7) = 1 and BMAT = G are incompatable            '
      ELSE IF (info.eq.-12) THEN
        IF (icall.eq.1) THEN
          string='IPARAM(1) must be equal to 0 or 1                    '
        ELSE
          string='HOWMANY = S not yet implemented                      '
        END IF
      ELSE IF (info.eq.-13) THEN
        string='HOWMANY must be one of A or P if Lrvec = .TRUE.        '
      ELSE IF (info.eq.-14) THEN
        string='Did not find any eigenvalues to sufficient accuaracy   '
      ELSE IF (info.eq.-15) THEN
        string='Different count of converge Ritz values in DNEUPD      '
      ELSE IF (info.eq.-9999) THEN
        string='Could not build and Arnoldi factorization              '
      END IF
!
      RETURN

References mod_storage::info.

iram_error() [2/2]

subroutine private roms_kernel_mod::iram_error ( integer, intent(in) info, character (len=*), intent(out) string )

private

Definition at line 837 of file fsv_roms.h.

!
!=======================================================================
!                                                                      !
!  This routine decodes internal error messages from the Implicit      !
!  Restarted Arnoldi Method (IRAM) for the computation of optimal      !
!  perturbation Ritz eigenfunctions.                                   !
!                                                                      !
!=======================================================================
!
!  imported variable declarations.
!
      integer, intent(in) :: info
!
      character (len=*), intent(out) :: string
!
!-----------------------------------------------------------------------
!  Decode error message from IRAM.
!-----------------------------------------------------------------------
!
      IF (info.eq.0)  THEN
        string='Normal exit                                            '
      ELSE IF (info.eq.1) THEN
        string='Maximum number of iterations taken                     '
      ELSE IF (info.eq.3) THEN
        string='No shifts could be applied during an IRAM cycle        '
      ELSE IF (info.eq.-1) THEN
        string='Nstate must be positive                                '
      ELSE IF (info.eq.-2) THEN
        string='NEV must be positive                                   '
      ELSE IF (info.eq.-3) THEN
        string='NCV must be greater NEV and less than or equal Nstate  '
      ELSE IF (info.eq.-4) THEN
        string='Maximum number of iterations must be greater than zero '
      ELSE IF (info.eq.-5) THEN
        string='WHICH must be one of LM, SM, LA, SA or BE              '
      ELSE IF (info.eq.-6) THEN
        string='BMAT must be one of I or G                             '
      ELSE IF (info.eq.-7) THEN
        string='Length of private work array SworkL is not sufficient  '
      ELSE IF (info.eq.-8) THEN
        string='Error in DSTEQR in the eigenvalue calculation          '
      ELSE IF (info.eq.-9) THEN
        string='Starting vector is zero                                '
      ELSE IF (info.eq.-10) THEN
        string='IPARAM(7) must be 1, 2, 3, 4, 5                        '
      ELSE IF (info.eq.-11) THEN
        string='IPARAM(7) = 1 and BMAT = G are incompatable            '
      ELSE IF (info.eq.-12) THEN
        string='IPARAM(1) must be equal to 0 or 1                      '
      ELSE IF (info.eq.-13) THEN
        string='NEV and WHICH = BE are incompatable                    '
      ELSE IF (info.eq.-14) THEN
        string='Did not find any eigenvalues to sufficient accuaracy   '
      ELSE IF (info.eq.-15) THEN
        string='HOWMANY must be one of A or S if RVEC = .TRUE.         '
      ELSE IF (info.eq.-16) THEN
        string='HOWMANY = S not yet implemented                        '
      ELSE IF (info.eq.-17) THEN
        string='Different count of converge Ritz values in DSEUPD      '
      ELSE IF (info.eq.-9999) THEN
        string='Could not build and Arnoldi factorization              '
      END IF
!
      RETURN

References mod_storage::info.

nlm_initial()

subroutine, private roms_kernel_mod::nlm_initial ( integer, intent(in) phase )

private

Definition at line 529 of file jedi_roms.h.

!
!=======================================================================
!                                                                      !
!  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

References analytical_mod::ana_drag(), analytical_mod::ana_psource(), analytical_mod::ana_sponge(), mod_scalars::avgke, mod_scalars::avgkp, mod_scalars::avgpe, check_multifile(), close_io_mod::close_inp(), mod_scalars::compositegrid, mod_scalars::day2sec, mod_scalars::dstart, mod_scalars::dt, mod_scalars::erstr, mod_scalars::exit_flag, mod_parallel::first_tile, mod_scalars::first_time, strings_mod::founderror(), get_data(), get_idata(), grid_coords(), mod_param::gridnumber, mod_param::gridsinlayer, mod_iounits::his, mod_scalars::iic, mod_scalars::iif, mod_scalars::indx1, ini_hmixcoef_mod::ini_hmixcoef(), mod_scalars::initime, mod_param::inlm, mod_scalars::jic, mod_stepping::knew, mod_stepping::krhs, mod_stepping::kstp, mod_parallel::last_tile, mod_scalars::ldefhis, mod_scalars::ldefqck, mod_scalars::ldefrst, mod_scalars::lsponge, mod_scalars::lstiffness, mod_scalars::ltracersrc, mod_scalars::luvsrc, mod_scalars::lwrthis, mod_scalars::lwrtqck, mod_scalars::lwrtrst, mod_scalars::lwsrc, mod_parallel::master, mod_nesting::n2dic, mod_nesting::n3dic, nesting_mod::nesting(), mod_param::nestlayers, mod_stepping::nf, mod_stepping::nfm1, mod_stepping::nfm2, mod_stepping::nfm3, mod_stepping::nfp1, mod_nesting::nfsic, mod_param::ngrids, mod_nesting::nmask, mod_stepping::nnew, mod_scalars::noerror, mod_stepping::nrhs, mod_scalars::nrrec, mod_scalars::nrun, mod_stepping::nstp, mod_scalars::ntend, mod_scalars::ntfirst, mod_scalars::ntimes, mod_scalars::ntstart, mod_nesting::ntvic, mod_nesting::nzwgt, omega_mod::omega(), post_initial_mod::post_initial(), mod_scalars::predictor_2d_step, mod_iounits::qck, rho_eos_mod::rho_eos(), mod_iounits::rst, set_data(), set_depth_mod::set_depth(), set_depth_mod::set_depth0(), set_masks_mod::set_masks(), set_massflux_mod::set_massflux(), mod_iounits::stdout, mod_scalars::step_counter, stiffness_mod::stiffness(), mod_scalars::synchro_flag, mod_scalars::tdays, mod_scalars::time, mod_scalars::time4jedi, mod_scalars::time_code, dateclock_mod::time_string(), mod_scalars::volume, wclock_off(), wclock_on(), and wetdry_mod::wetdry().

Referenced by roms_initializep1(), roms_initializep2(), and roms_initializep3().

Here is the call graph for this function:

Here is the caller graph for this function:

roms_finalize()

subroutine public roms_kernel_mod::roms_finalize

Definition at line 282 of file ad_roms.h.

!
!=======================================================================
!                                                                      !
!  This routine terminates ROMS adjoint model 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 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, iadm, 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, iadm)
      END DO
      CALL close_out
!
      RETURN

References mod_scalars::blowup, close_io_mod::close_inp(), close_io_mod::close_out(), mod_scalars::exit_flag, mod_param::iadm, mod_scalars::lcyclerst, mod_scalars::lwrtrst, mod_parallel::master, memory(), mod_parallel::myrank, mod_param::ngrids, mod_scalars::noerror, mod_iounits::rst, mod_iounits::stdout, wclock_off(), and wrt_rst_mod::wrt_rst().

Referenced by mct_driver(), myroms(), cmeps_roms_mod::roms_modeladvance(), esmf_roms_mod::roms_modeladvance(), cmeps_roms_mod::roms_setfinalize(), and esmf_roms_mod::roms_setfinalize().

Here is the call graph for this function:

Here is the caller graph for this function:

roms_initialize()

subroutine public roms_kernel_mod::roms_initialize	(	logical, intent(inout)	first,
		integer, intent(in), optional	mpicomm )

Definition at line 51 of file ad_roms.h.

!
!=======================================================================
!                                                                      !
!  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 :: 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 (iadm)
        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, iadm, 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
!
!-----------------------------------------------------------------------
!  Initialize adjoint model state variables over all nested grids, if
!  applicable.
!-----------------------------------------------------------------------
 
#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.
      END DO
#endif
!
!  Initialize adjoint model.
!
      lstiffness=.false.
      DO ng=1,ngrids
        lreadfwd(ng)=.true.
        CALL ad_initial (ng)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END DO
!
!  Initialize run or ensemble counter.
!
      nrun=1
!
!  Activate adjoint output.
!
      DO ng=1,ngrids
        ldefadj(ng)=.true.
        lwrtadj(ng)=.true.
        lcycleadj(ng)=.false.
      END DO
!
      RETURN

References ad_initial(), edit_multifile(), mod_scalars::exit_flag, mod_parallel::first_tile, strings_mod::founderror(), mod_param::iadm, mod_parallel::initialize_parallel(), inp_par_mod::inp_par(), mod_parallel::last_tile, mod_scalars::lcycleadj, mod_scalars::ldefadj, mod_scalars::lreadblk, mod_scalars::lreadfrc, mod_scalars::lreadfwd, mod_scalars::lstiffness, mod_scalars::lwrtadj, mod_parallel::master, my_threadnum(), mod_parallel::myrank, mod_parallel::mythread, mod_param::ngrids, mod_scalars::noerror, mod_scalars::nrun, mod_param::ntilee, mod_param::ntilex, mod_parallel::numthreads, mod_parallel::ocn_comm_world, mod_arrays::roms_allocate_arrays(), mod_arrays::roms_initialize_arrays(), stdout_mod::set_stdoutunit, mod_iounits::stdout, stdout_mod::stdout_unit(), and wclock_on().

Referenced by mct_driver(), myroms(), cmeps_roms_mod::roms_setinitializep2(), and esmf_roms_mod::roms_setinitializep2().

Here is the call graph for this function:

Here is the caller graph for this function:

roms_initializep1()

subroutine, public roms_kernel_mod::roms_initializep1	(	logical, intent(inout)	first,
		integer, intent(in), optional	mpicomm,
		integer, intent(in), optional	kernel )

Definition at line 95 of file jedi_roms.h.

!
!=======================================================================
!                                                                      !
!  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

References adm_initial(), mod_scalars::exit_flag, mod_parallel::first_tile, strings_mod::founderror(), mod_param::iadm, mod_parallel::initialize_parallel(), mod_param::inlm, inp_par_mod::inp_par(), mod_param::itlm, mod_parallel::last_tile, mod_parallel::master, mod_parallel::myrank, mod_parallel::mythread, mod_param::ngrids, nlm_initial(), mod_scalars::noerror, mod_scalars::nrun, mod_param::ntilee, mod_param::ntilex, mod_parallel::numthreads, mod_parallel::ocn_comm_world, mod_arrays::roms_allocate_arrays(), mod_arrays::roms_initialize_arrays(), set_grid(), stdout_mod::set_stdoutunit, mod_scalars::setgridconfig, mod_iounits::stdout, stdout_mod::stdout_unit(), tlm_initial(), and wclock_on().

Here is the call graph for this function:

roms_initializep2()

subroutine, public roms_kernel_mod::roms_initializep2

(

integer, intent(in)

kernel

)

Definition at line 249 of file jedi_roms.h.

!
!=======================================================================
!                                                                      !
!  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

References adm_initial(), mod_param::iadm, mod_param::inlm, mod_param::itlm, nlm_initial(), and tlm_initial().

Here is the call graph for this function:

roms_initializep3()

subroutine, public roms_kernel_mod::roms_initializep3

(

integer, intent(in)

kernel

)

Definition at line 289 of file jedi_roms.h.

!
!=======================================================================
!                                                                      !
!  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

References adm_initial(), mod_param::iadm, mod_param::inlm, mod_param::itlm, nlm_initial(), and tlm_initial().

Here is the call graph for this function:

roms_run() [1/3]

subroutine public roms_kernel_mod::roms_run

(

real(dp), intent(in)

runinterval

)

Definition at line 238 of file ad_roms.h.

!
!=======================================================================
!                                                                      !
!  This routine runs ROMS adjoint model backwards for the              !
!  specified time interval (seconds), RunInterval.                     !
!                                                                      !
!=======================================================================
!
!  Imported variable declarations.
!
      real(dp), intent(in) :: RunInterval            ! seconds
!
!  Local variable declarations.
!
      integer :: ng
!
      character (len=*), parameter :: MyFile =                          &
     &  __FILE__//", ROMS_run"
!
!-----------------------------------------------------------------------
!  Time-step adjoint model over all nested grids, if applicable.
!-----------------------------------------------------------------------
!
      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
!
 10   FORMAT (/,1x,a,1x,'ROMS: started time-stepping:',                 &
     &        ' (Grid: ',i2.2,' TimeSteps: ',i8.8,' - ',i8.8,')',/)
!
      RETURN

References ad_main2d(), ad_main3d(), mod_scalars::exit_flag, strings_mod::founderror(), mod_parallel::master, mod_param::ngrids, mod_scalars::noerror, mod_scalars::ntend, mod_scalars::ntstart, and mod_iounits::stdout.

Referenced by mct_driver(), myroms(), cmeps_roms_mod::roms_modeladvance(), and esmf_roms_mod::roms_modeladvance().

Here is the call graph for this function:

Here is the caller graph for this function:

roms_run() [2/3]

subroutine, public roms_kernel_mod::roms_run

(

real(dp), intent(inout)

runinterval

)

Definition at line 493 of file obs_sen_rbl4dvar_forecast.h.

!
!=======================================================================
!                                                                      !
!  This routine time-steps ROMS nonlinear, tangent linear and          !
!  adjoint models.                                                     !
!                                                                      !
!=======================================================================
!
!  Imported variable declarations
!
      real(dp), intent(inout) :: RunInterval            ! seconds
!
!  Local variable declarations.
!
      logical :: Lcgini, Linner, Lposterior, add
!
      integer :: my_inner, my_outer
      integer :: Lbck, Lini, Rec1, Rec2, ImpOrd
      integer :: i, lstr, ng, status, tile
      integer :: Fcount, NRMrec
 
      integer, dimension(Ngrids) :: indxSave
      integer, dimension(Ngrids) :: Nrec
!
      real(r8) :: str_day, end_day, dstartS, rtime
!
      character (len=1) :: charA, charB
#ifdef OBS_SPACE
      character (len=1) :: charC
#endif
      character (len=25) :: driver
      character (len=20) :: string
 
      character (len=*), parameter :: MyFile =                          &
     &  __FILE__//", ROMS_run"
!
!=======================================================================
!  Run model for all nested grids, if any.
!=======================================================================
!
!  Initialize relevant parameters.
!
      DO ng=1,ngrids
#if 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
        lold(ng)=1          ! old minimization time index
        lnew(ng)=2          ! new minimization time index
      END DO
      lini=1                ! NLM initial conditions record in INI
      lbck=2                ! background record in INI
      rec1=1
      rec2=2
      nrun=1
      outer=0
      inner=0
      erstr=1
      erend=nouter
      driver='obs_sen_rbl4dvar_forecast'
      chara='A'
      charb='B'
#ifdef OBS_SPACE
      charc='C'
#endif
!
!  Choose the desired order approximation of the forecast impact
!  calculations:
!
!          ImpOrd: 1 (1st order)
!                  2 (2nd order) or
!                  3 (3rd order)
!
!  ImpOrd=1 is the least accurate and can be in error by a factor of 2.
!
      impord=3
!
!  For this driver DSTART is assumed to be the start time of the
!  analysis cycle that yields the analysis for the subsequent forecast
!  cycle. The number of time steps in the analysis cycle is NTIMES_ANA,
!  and the number of time steps in the forecast cycle is NTIMES_FCT.
!
!  Initially we will be running the adjoint model over the forecast
!  time interval so we temporarily reset dstart and RunInterval to
!  reflect the forecast interval start and run time.
!
      dstarts=dstart
      rtime=0.0_r8
      DO ng=1,ngrids
        rtime=max(rtime, dt(ng)*ntimes_ana(ng))
      END DO
      dstart=dstart+rtime*sec2day
!
      rtime=0.0_r8
      DO ng=1,ngrids
        rtime=max(rtime, dt(ng)*ntimes_fct(ng))
        ntimes(ng)=ntimes_fct(ng)
      END DO
      runinterval=rtime
!
!-----------------------------------------------------------------------
!  Configure weak constraint RBL4D-Var algorithm.
!-----------------------------------------------------------------------
!
!  Initialize the switch to gather weak constraint forcing.
!
      DO ng=1,ngrids
        wrtforce(ng)=.false.
      END DO
!
!  Initialize and set nonlinear model initial conditions.
!
      DO ng=1,ngrids
        lreadfwd(ng)=.false.
        lreadblk(ng)=.false.
        wrtnlmod(ng)=.true.
        wrtrpmod(ng)=.false.
        wrttlmod(ng)=.false.
        lsen4dvar(ng)=.false.
#ifdef OBS_SPACE
        lobspace(ng)=.false.
# ifndef OBS_IMPACT
        ladjvar(ng)=.false.
# endif
#endif
      END DO
 
      CALL initial
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Save nonlinear initial conditions (currently in time index 1,
!  background) into record "Lbck" of INI(ng)%name NetCDF file. The
!  record "Lbck" becomes the background state record and the record
!  "Lini" becomes current nonlinear initial conditions.
!
      DO ng=1,ngrids
        ini(ng)%Rindex=1
        fcount=ini(ng)%load
        ini(ng)%Nrec(fcount)=1
#ifdef DISTRIBUTE
        CALL wrt_ini (ng, myrank, 1)
#else
        CALL wrt_ini (ng, -1, 1)
#endif
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END DO
!
!  The FWD structure is used several times and contains the nonnlinear
!  background trajectories (RBL4DVAR, FCTA, and FCTB) used to linearize
!  the tangent linear and adjoint models.  These trajectories can be
!  split into multi-files. If so, the user needs to specified such
!  multiple files in the standard input script (roms.in). Since the HIS
!  structure is not used here, copy the original FWD containing the
!  trajectory information loaded during configuration in "inp_par",
!  which has the values for the regular RBL4D-Var at specified outer
!  loop into the HIS structure.
!
      CALL edit_multifile ('FWD2HIS')
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
!  Model-error covariance normalization and stardard deviation factors.
!:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
!
!  Compute or read in the error covariance normalization factors.
!  If computing, write out factors to NetCDF. This is an expensive
!  computation that needs to be computed only once for a particular
!  application grid and decorrelation scales.
!
      DO ng=1,ngrids
        IF (any(lwrtnrm(:,ng))) THEN
          CALL def_norm (ng, inlm, 1)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
          IF (nsa.eq.2) THEN
            CALL def_norm (ng, inlm, 2)
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
          END IF
 
#ifdef ADJUST_BOUNDARY
          CALL def_norm (ng, inlm, 3)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
#endif
 
#if defined ADJUST_WSTRESS || defined ADJUST_STFLUX
          CALL def_norm (ng, inlm, 4)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
#endif
 
          DO tile=first_tile(ng),last_tile(ng),+1
            CALL normalization (ng, tile, 2)
          END DO
 
          ldefnrm(1:4,ng)=.false.
          lwrtnrm(1:4,ng)=.false.
        ELSE
          nrmrec=1
          CALL get_state (ng, 14, 14, nrm(1,ng), nrmrec, 1)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
          IF (nsa.eq.2) THEN
            CALL get_state (ng, 15, 15, nrm(2,ng), nrmrec, 2)
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
          END IF
 
#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 IF
      END DO
 
#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
          CALL get_state (ng, inlm, 2, ini(ng), lini, lini)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
          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
!
!  Define tangent linear initial conditions file.
!
      DO ng=1,ngrids
        ldefitl(ng)=.true.
        CALL tl_def_ini (ng)
        ldefitl(ng)=.false.
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END DO
!
!  Define impulse forcing NetCDF file.
!
      DO ng=1,ngrids
        ldeftlf(ng)=.true.
        CALL def_impulse (ng)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END DO
 
#ifndef OBS_SPACE
!
!  Define output 4DVAR NetCDF file containing all processed data
!  at observation locations.
!
      DO ng=1,ngrids
        ldefmod(ng)=.true.
        CALL def_mod (ng)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      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
#endif
!
!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
!  Run the adjoint model forced by forecast minus analysis difference
!  using appropriate solution trajectories.
!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
!
!  Reset the start and end times for the adjoint forcing.
!
      DO ng=1,ngrids
        str_day=tdays(ng)+ntimes(ng)*dt(ng)*sec2day
        end_day=tdays(ng)
        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,70) 'AD', dends(ng), dstrs(ng)
        END IF
      END DO
!
!  Set basic state trajectory and adjoint forcing file.
!
!     ads_xxxx_a.nc    if obs_space is off
!     ads_xxxx_b.nc
!
!     obs_xxxx_a.nc    if obs_space is on
!     obs_xxxx_b.nc
!
      sourcefile=myfile
      DO ng=1,ngrids
#ifdef OBS_SPACE
        CALL close_file (ng, inlm, obs(ng), obs(ng)%name)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
        IF (impord.ne.2) THEN
          WRITE (obs(ng)%name,90) trim(oifb(ng)%head)
        ELSE
          WRITE (obs(ng)%name,90) trim(oifa(ng)%head)
        END IF
#else
        CALL close_file (ng, inlm, ads(ng), ads(ng)%name)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
        IF (impord.ne.2) THEN
          WRITE (ads(ng)%name,90) trim(foib(ng)%head)
        ELSE
          WRITE (ads(ng)%name,90) trim(foia(ng)%head)
        END IF
#endif
      END DO
!
!  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 ('FCTB2FWD')
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      DO ng=1,ngrids
        lreadfwd(ng)=.true.
      END DO
!
!  Set structure for the nonlinear surface fluxes to be processed by
!  by the tangent linear and adjoint models. Also, set switches to
!  process the BLK structure in routine "check_multifile".  Notice that
!  it is possible to split solution into multiple NetCDF files to reduce
!  their size.
!
      CALL edit_multifile ('FCTB2BLK')
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      DO ng=1,ngrids
        lreadblk(ng)=.true.
      END DO
!
      IF (master) THEN
       WRITE (stdout,50)
      END IF
!
!  Initialize the adjoint model: initialize using dI/dxf is appropriate.
!
!  TODO: There is problem with the "haveObsMeta" flag. It is always
!        FALSE on the first call to OBS_READ, so "ObsMeta" is zero
!  for the radials. The reason this is happening is because that
!  NetCDF file had already been opened somewhere prior to calling
!  OBS_INITIAL in TL_INITIAL (and AD_INITIAL when appropriate),
!  so the obs arrays were not getting initialized properly by
!  these routines. For now, the solution is to close the NetCDF
!  file before each call to TL_INITIAL and AD_INITIAL.
!
      DO ng=1,ngrids
        lstiffness=.false.
#ifdef OBS_SPACE
        lsen4dvar(ng)=.false.
        lsenfct(ng)=.true.
        lobspace(ng)=.true.
# ifndef OBS_IMPACT
        ladjvar(ng)=.false.
# endif
#else
        lsen4dvar(ng)=.true.
        lsenfct(ng)=.false.
#endif
        CALL close_file (ng, iadm, obs(ng), obs(ng)%name)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
        CALL ad_initial (ng)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END DO
!
!  Set adjoint history NetCDF parameters.  Define adjoint history
!  file once to avoid opening to many files.
!
      DO ng=1,ngrids
        wrtforce=.true.
        IF (adm(ng)%ncid.ne.-1) ldefadj(ng)=.false.
        fcount=adm(ng)%load
        adm(ng)%Nrec(fcount)=0
        adm(ng)%Rindex=0
      END DO
 
!  Time-step adjoint model backwards.
!
      DO ng=1,ngrids
        IF (master) THEN
          WRITE (stdout,20) '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
!
!  Write out last weak-constraint forcing (WRTforce is still .TRUE.)
!  record into the adjoint history file.  Note that the weak-constraint
!  forcing is delayed by nADJ time-steps.
!
      DO ng=1,ngrids
#ifdef DISTRIBUTE
        CALL ad_wrt_his (ng, myrank)
#else
        CALL ad_wrt_his (ng, -1)
#endif
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END DO
!
!  Write out adjoint initial condition record into the adjoint
!  history file.
!
      DO ng=1,ngrids
        wrtforce(ng)=.false.
        lwrtstate2d(ng)=.false.
#ifdef DISTRIBUTE
        CALL ad_wrt_his (ng, myrank)
#else
        CALL ad_wrt_his (ng, -1)
#endif
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END DO
!
!   Retrieve adjoint solution into tl_var(1).
!
      DO ng=1,ngrids
         CALL get_state (ng, itlm, 4, adm(ng), adm(ng)%Rindex, rec1)
         IF (founderror(exit_flag, noerror, __line__,  myfile)) RETURN
      END DO
!
!  Set basic state trajectory and adjoint forcing file.
!
      IF (impord.gt.1) THEN
        sourcefile=myfile
        DO ng=1,ngrids
#ifdef OBS_SPACE
          CALL close_file (ng, inlm, obs(ng), obs(ng)%name)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
          IF (impord.eq.2) THEN
             WRITE (obs(ng)%name,90) trim(oifb(ng)%head)
          ELSE
             WRITE (obs(ng)%name,90) trim(oifa(ng)%head)
          END IF
#else
          CALL close_file (ng, inlm, ads(ng), ads(ng)%name)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
          IF (impord.eq.2) THEN
             WRITE (ads(ng)%name,90) trim(foib(ng)%head)
          ELSE
             WRITE (ads(ng)%name,90) trim(foia(ng)%head)
          END IF
#endif
        END DO
!
!  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 ('FCTA2FWD')
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        DO ng=1,ngrids
          lreadfwd(ng)=.true.
        END DO
!
!  Set structure for the nonlinear surface fluxes to be processed by
!  by the tangent linear and adjoint models. Also, set switches to
!  process the BLK structure in routine "check_multifile".  Notice that
!  it is possible to split solution into multiple NetCDF files to reduce
!  their size.
!
        CALL edit_multifile ('FCTA2BLK')
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        DO ng=1,ngrids
          lreadblk(ng)=.true.
        END DO
!
        IF (master) THEN
          WRITE (stdout,50)
        END IF
!
!  Initialize the adjoint model: initialize using dI/dxf is appropriate.
!
        DO ng=1,ngrids
          lstiffness=.false.
#ifdef OBS_SPACE
          lsen4dvar(ng)=.false.
          lsenfct(ng)=.true.
          lobspace(ng)=.true.
# ifndef OBS_IMPACT
          ladjvar(ng)=.false.
# endif
#else
          lsen4dvar(ng)=.true.
          lsenfct(ng)=.false.
#endif
          CALL close_file (ng, iadm, obs(ng), obs(ng)%name)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
          CALL ad_initial (ng)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END DO
!
!  Set adjoint history NetCDF parameters.  Define adjoint history
!  file once to avoid opening to many files.
!
        DO ng=1,ngrids
          wrtforce=.true.
          IF (adm(ng)%ncid.ne.-1) ldefadj(ng)=.false.
          fcount=adm(ng)%load
          adm(ng)%Nrec(fcount)=0
          adm(ng)%Rindex=0
        END DO
!
!  Time-step adjoint model backwards.
!
        DO ng=1,ngrids
          IF (master) THEN
            WRITE (stdout,20) '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
!
!  Write out last weak-constraint forcing (WRTforce is still .TRUE.)
!  record into the adjoint history file.  Note that the weak-constraint
!  forcing is delayed by nADJ time-steps.
!
        DO ng=1,ngrids
#ifdef DISTRIBUTE
          CALL ad_wrt_his (ng, myrank)
#else
          CALL ad_wrt_his (ng, -1)
#endif
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END DO
!
!  Write out adjoint initial condition record into the adjoint
!  history file.
!
        DO ng=1,ngrids
          wrtforce(ng)=.false.
          lwrtstate2d(ng)=.false.
#ifdef DISTRIBUTE
          CALL ad_wrt_his (ng, myrank)
#else
          CALL ad_wrt_his (ng, -1)
#endif
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END DO
!
!   Retrieve adjoint solution.
!
        DO ng=1,ngrids
          CALL get_state (ng, iadm, 4, adm(ng), adm(ng)%Rindex, rec1)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END DO
!
!   Add the retrieved adjoint solution to the previous solution saved
!   in tl_var(1).
!
        add=.true.
        DO ng=1,ngrids
          DO tile=first_tile(ng),last_tile(ng),+1
            CALL load_adtotl (ng, tile, rec1, rec1, add)
          END DO
        END DO
!
      END IF
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!  Adjoint of RBL4D-Var to compute the observation sensitivity.
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
!  Reset the start and end times for the adjoint forcing.
!
      DO ng=1,ngrids
        str_day=tdays(ng)+ntimes(ng)*dt(ng)*sec2day
        end_day=tdays(ng)
        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,70) 'AD', dends(ng), dstrs(ng)
        END IF
      END DO
!
!  WARNING: ONLY one outer loop can be used for this application.
!  =======  For more than 1 outer-loop, we require the second
!  derivative of each model operator (i.e. the tangent linear
!  of the tangent linear operator).
!
      ad_outer_loop : DO my_outer=nimpact,nimpact
        outer=my_outer
        inner=0
!
!-----------------------------------------------------------------------
!  Run the adjoint model initialized and forced by dI/dx where I is the
!  chosen function of the analysis/forecast state x.
!-----------------------------------------------------------------------
!
!  Reset DSTART and RunIterval for the analysis cycle interval using
!  NTIMES_ANA.
!
        dstart=dstarts
!
        rtime=0.0_r8
        DO ng=1,ngrids
          rtime=max(rtime, dt(ng)*ntimes_ana(ng))
          ntimes(ng)=ntimes_ana(ng)
        END DO
        runinterval=rtime
!
#if defined ADJUST_WSTRESS || defined ADJUST_STFLUX || \
    defined adjust_boundary
!
!  Call initial again to reset the forcing and obc adjustment time
!  counters.
!
        DO ng=1,ngrids
          lreadfwd(ng)=.false.
          lreadblk(ng)=.false.
        END DO
!
        CALL initial
#endif
!
!  Set basic state trajectory.
!
        DO ng=1,ngrids
          WRITE (fwd(ng)%name,10) trim(his(ng)%head), outer-1
        END DO
!
!  Set structure for the nonlinear forward trajectory (from regular
!  RBL4DVAR) 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
!
        IF (master) THEN
          WRITE (stdout,50)
        END IF
!
!  Initialize the adjoint model: initialize using dI/dxf is
!  appropriate.
!
        DO ng=1,ngrids
          lstiffness=.false.
!!AMM
#ifdef OBS_SPACE
          lsen4dvar(ng)=.true.
          lsenfct(ng)=.false.
          lobspace(ng)=.false.
# ifndef OBS_IMPACT
          ladjvar(ng)=.false.
# endif
#else
          lsen4dvar(ng)=.false.
          lsenfct(ng)=.false.
#endif
!!AMM
          CALL close_file (ng, iadm, obs(ng), obs(ng)%name)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
          CALL ad_initial (ng)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END DO
!
!    Copy tl_var(1) solution to ad_var(1).
!    It is important to clear the tl surface forcing and boundary
!    arrays first since only the initial fields must be used to
!    initialize the next run of the adjoint.
!
        add=.false.
        DO ng=1,ngrids
          DO tile=first_tile(ng),last_tile(ng),+1
#if defined ADJUST_WSTRESS || defined ADJUST_STFLUX
            CALL initialize_forces (ng, tile, itlm)
#endif
#ifdef ADJUST_BOUNDARY
            CALL initialize_boundary (ng, tile, itlm)
#endif
            CALL load_tltoad (ng, tile, rec1, rec1, add)
          END DO
        END DO
 
#if defined ADJUST_WSTRESS || defined ADJUST_STFLUX || \
    defined adjust_boundary
!
!  Clear the adjoint forcing and open boundary 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
#endif
!
!  Set adjoint history NetCDF parameters.  Define adjoint history
!  file one to avoid opening to many files.
!
        DO ng=1,ngrids
          wrtforce=.true.
          IF (adm(ng)%ncid.ne.-1) ldefadj(ng)=.false.
          fcount=adm(ng)%load
          adm(ng)%Nrec(fcount)=0
          adm(ng)%Rindex=0
        END DO
!
!  Time-step adjoint model backwards.
!  ??? What do we do in the case of model error? Save forcing for TLM?
!
        DO ng=1,ngrids
          IF (master) THEN
            WRITE (stdout,20) '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
!
!  Write out last weak-constraint forcing (WRTforce is still .TRUE.)
!  record into the adjoint history file.  Note that the weak-constraint
!  forcing is delayed by nADJ time-steps.
!
        DO ng=1,ngrids
#ifdef DISTRIBUTE
          CALL ad_wrt_his (ng, myrank)
#else
          CALL ad_wrt_his (ng, -1)
#endif
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END DO
!
!  Write out adjoint initial condition record into the adjoint
!  history file.
!
        DO ng=1,ngrids
          wrtforce(ng)=.false.
          lwrtstate2d(ng)=.false.
#ifdef DISTRIBUTE
          CALL ad_wrt_his (ng, myrank)
#else
          CALL ad_wrt_his (ng, -1)
#endif
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END DO
!
!  Convolve adjoint trajectory with error covariances.
!
        lposterior=.false.
        CALL error_covariance (itlm, driver, outer, inner,              &
     &                         lbck, lini, lold, lnew,                  &
     &                         rec1, rec2, lposterior)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Convert the current adjoint solution in ADM(ng)%name to impulse
!  forcing. Write out impulse forcing into TLF(ng)%name NetCDF file.
!  To facilitate the forcing to the TLM and RPM, the forcing is
!  processed and written in increasing time coordinates (recall that
!  the adjoint solution in ADM(ng)%name is backwards in time).
!
!  AMM: Do not know what to do in the weak constraint case yet.
!
!!      IF (Master) THEN
!!        WRITE (stdout,40) outer, inner
!!      END IF
!!      DO ng=1,Ngrids
!!        TLF(ng)%Rindex=0
#ifdef DISTRIBUTE
!!        CALL wrt_impulse (ng, MyRank, iADM, ADM(ng)%name)
#else
!!        CALL wrt_impulse (ng, -1, iADM, ADM(ng)%name)
#endif
!!        IF (FoundError(exit_flag, NoError, __LINE__, MyFile)) RETURN
!!      END DO
!
!:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
!  Integrate tangent linear model forced by the convolved adjoint
!  trajectory.
!:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
!
#ifdef OBS_SPACE
!
!  Reinitialize fourdvar arrays using the analysis interval observation
!  file.
!
        DO ng=1,ngrids
          CALL close_file (ng, inlm, obs(ng), obs(ng)%name)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
          WRITE (obs(ng)%name,80) trim(obs(ng)%head), charc
        END DO
!
        CALL deallocate_fourdvar
        CALL initialize_fourdvar
!
!-----------------------------------------------------------------------
!  If the required vectors and arrays from congrad from a previous
!  run of the assimilation cycle are available, read them here from
!  LCZ(ng)%name NetCDF file.
!-----------------------------------------------------------------------
!
        sourcefile=myfile
        DO ng=1,ngrids
          SELECT CASE (lcz(ng)%IOtype)
            CASE (io_nf90)
              CALL netcdf_get_fvar (ng, itlm, lcz(ng)%name,             &
     &                              'cg_beta', cg_beta)
              IF (founderror(exit_flag, noerror,                        &
     &                       __line__, myfile)) RETURN
 
              CALL netcdf_get_fvar (ng, itlm, lcz(ng)%name,             &
     &                              'cg_delta', cg_delta)
              IF (founderror(exit_flag, noerror,                        &
     &                       __line__, myfile)) RETURN
 
              CALL netcdf_get_fvar (ng, itlm, lcz(ng)%name,             &
     &                              'cg_Gnorm_v', cg_gnorm_v)
              IF (founderror(exit_flag, noerror,                        &
     &                       __line__, myfile)) RETURN
 
              CALL netcdf_get_fvar (ng, itlm, lcz(ng)%name,             &
     &                              'cg_dla', cg_dla)
              IF (founderror(exit_flag, noerror,                        &
     &                       __line__, myfile)) RETURN
 
              CALL netcdf_get_fvar (ng, itlm, lcz(ng)%name,             &
     &                               'cg_QG', cg_qg)
              IF (founderror(exit_flag, noerror,                        &
     &                       __line__, myfile)) RETURN
 
              CALL netcdf_get_fvar (ng, itlm, lcz(ng)%name,             &
     &                              'zgrad0', zgrad0)
              IF (founderror(exit_flag, noerror,                        &
     &                       __line__, myfile)) RETURN
 
              CALL netcdf_get_fvar (ng, itlm, lcz(ng)%name,             &
     &                              'zcglwk', zcglwk)
              IF (founderror(exit_flag, noerror,                        &
     &                       __line__, myfile)) RETURN
 
              CALL netcdf_get_fvar (ng, itlm, lcz(ng)%name,             &
     &                              'TLmodVal_S', tlmodval_s,           &
     &                              broadcast = .false.)   ! Master only
              IF (founderror(exit_flag, noerror,                        &
     &                       __line__, myfile)) RETURN
 
# ifdef RPCG
              CALL netcdf_get_fvar (ng, itlm, lcz(ng)%name,             &
     &                              'Hbk', hbk)
              IF (founderror(exit_flag, noerror,                        &
     &                       __line__, myfile)) RETURN
 
              CALL netcdf_get_fvar (ng, itlm, lcz(ng)%name,             &
     &                              'Jb0', jb0)
              IF (founderror(exit_flag, noerror,                        &
     &                       __line__, myfile)) RETURN
 
              CALL netcdf_get_fvar (ng, itlm, lcz(ng)%name,             &
     &                              'vcglwk', vcglwk)
              IF (founderror(exit_flag, noerror,                        &
     &                       __line__, myfile)) RETURN
# endif
 
# if defined PIO_LIB && defined DISTRIBUTE
            CASE (io_pio)
              CALL pio_netcdf_get_fvar (ng, itlm, lcz(ng)%name,         &
     &                                  'cg_beta', cg_beta)
              IF (founderror(exit_flag, noerror,                        &
     &                       __line__, myfile)) RETURN
 
              CALL pio_netcdf_get_fvar (ng, itlm, lcz(ng)%name,         &
     &                                  'cg_delta', cg_delta)
              IF (founderror(exit_flag, noerror,                        &
     &                       __line__, myfile)) RETURN
 
              CALL pio_netcdf_get_fvar (ng, itlm, lcz(ng)%name,         &
     &                                  'cg_Gnorm_v', cg_gnorm_v)
              IF (founderror(exit_flag, noerror,                        &
     &                       __line__, myfile)) RETURN
 
              CALL pio_netcdf_get_fvar (ng, itlm, lcz(ng)%name,         &
     &                                  'cg_dla', cg_dla)
              IF (founderror(exit_flag, noerror,                        &
     &                       __line__, myfile)) RETURN
 
              CALL pio_netcdf_get_fvar (ng, itlm, lcz(ng)%name,         &
     &                                  'cg_QG', cg_qg)
              IF (founderror(exit_flag, noerror,                        &
     &                       __line__, myfile)) RETURN
 
              CALL pio_netcdf_get_fvar (ng, itlm, lcz(ng)%name,         &
     &                                  'zgrad0', zgrad0)
              IF (founderror(exit_flag, noerror,                        &
     &                       __line__, myfile)) RETURN
 
              CALL pio_netcdf_get_fvar (ng, itlm, lcz(ng)%name,         &
     &                                  'zcglwk', zcglwk)
              IF (founderror(exit_flag, noerror,                        &
     &                       __line__, myfile)) RETURN
 
              CALL pio_netcdf_get_fvar (ng, itlm, lcz(ng)%name,         &
     &                                  'TLmodVal_S', tlmodval_s)
              IF (founderror(exit_flag, noerror,                        &
     &                       __line__, myfile)) RETURN
 
#  ifdef RPCG
              CALL pio_netcdf_get_fvar (ng, itlm, lcz(ng)%name,         &
     &                                  'Hbk', hbk)
              IF (founderror(exit_flag, noerror,                        &
     &                       __line__, myfile)) RETURN
 
              CALL pio_netcdf_get_fvar (ng, itlm, lcz(ng)%name,         &
     &                                  'Jb0', jb0)
              IF (founderror(exit_flag, noerror,                        &
     &                       __line__, myfile)) RETURN
 
              CALL pio_netcdf_get_fvar (ng, itlm, lcz(ng)%name,         &
     &                                  'vcglwk', vcglwk)
              IF (founderror(exit_flag, noerror,                        &
     &                       __line__, myfile)) RETURN
#  endif
# endif
          END SELECT
        END DO
!
!-----------------------------------------------------------------------
!  If skiping runing nonlinear model, read in observation screening and
!  quality control flag.
!-----------------------------------------------------------------------
!
        sourcefile=myfile
        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 pio_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
!
!  Define output 4DVAR NetCDF file containing all processed data
!  at observation locations.
!
        DO ng=1,ngrids
          ldefmod(ng)=.true.
          CALL def_mod (ng)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END DO
#endif /* OBS_SPACE */
!
        DO ng=1,ngrids
          wrtnlmod(ng)=.false.
          wrttlmod(ng)=.true.
          lwrttlm(ng)=.false.
        END DO
!
!  Clear tangent linear forcing arrays before entering inner-loop.
!  This is very important.
!
        DO ng=1,ngrids
          DO tile=first_tile(ng),last_tile(ng),+1
            CALL initialize_forces (ng, tile, itlm)
#ifdef ADJUST_BOUNDARY
            CALL initialize_boundary (ng, tile, itlm)
#endif
          END DO
        END DO
!
!  Set basic state trajectory.
!
        DO ng=1,ngrids
          WRITE (fwd(ng)%name,10) trim(his(ng)%head), outer-1
          lstr=len_trim(fwd(ng)%name)
          fwd(ng)%base=fwd(ng)%name(1:lstr-3)
        END DO
!
!  If weak constraint, the impulses are time-interpolated at each
!  time-steps.
!
        DO ng=1,ngrids
          IF (frcrec(ng).gt.3) THEN
            frequentimpulse(ng)=.true.
          END IF
        END DO
!
!  Initialize tangent linear model from ITL(ng)%name, record Rec1.
!
        DO ng=1,ngrids
          itl(ng)%Rindex=rec1
          CALL close_file (ng, itlm, obs(ng), obs(ng)%name)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
          CALL tl_initial (ng)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END DO
!
!  Run tangent linear model forward and force with convolved adjoint
!  trajectory impulses. Compute (HMBM'H')_n * PSI at observation points
!  which are used in the conjugate gradient algorithm.
!
        DO ng=1,ngrids
          IF (master) THEN
            WRITE (stdout,20) '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
!
        DO ng=1,ngrids
          wrtnlmod(ng)=.false.
          wrttlmod(ng)=.false.
        END DO
 
#ifdef OBS_IMPACT
!
!  Compute observation impact to the data assimilation system.
!
        DO ng=1,ngrids
# ifdef RPCG
          CALL rep_matrix (ng, itlm, outer, ninner-1)
# else
          CALL rep_matrix (ng, itlm, outer, ninner)
# endif
        END DO
#else
!
!  Set basic state trajectory for adjoint inner-loops.
!
        DO ng=1,ngrids
          WRITE (fwd(ng)%name,10) trim(his(ng)%head), outer-1
          lstr=len_trim(fwd(ng)%name)
          fwd(ng)%base=fwd(ng)%name(1:lstr-3)
        END DO
!
!  Clear tangent linear forcing arrays before entering inner-loop.
!
        DO ng=1,ngrids
          DO tile=first_tile(ng),last_tile(ng),+1
            CALL initialize_forces (ng, tile, itlm)
# ifdef ADJUST_BOUNDARY
            CALL initialize_boundary (ng, tile, itlm)
# endif
          END DO
        END DO
!
# ifdef RPCG
        ad_inner_loop : DO my_inner=ninner,0,-1
# else
        ad_inner_loop : DO my_inner=ninner,1,-1
# endif
          inner=my_inner
 
# ifdef RPCG
!
!  Retrieve NLmodVal when inner=0 for use as BCKmodVal.
!
          IF (inner.eq.0) THEN
            DO ng=1,ngrids
              SELECT CASE (dav(ng)%IOtype)
                CASE (io_nf90)
                  CALL netcdf_get_fvar (ng, itlm, dav(ng)%name,         &
     &                                  'NLmodel_value', nlmodval)
                  IF (founderror(exit_flag, noerror,                    &
     &                           __line__, myfile)) RETURN
 
#  if defined PIO_LIB && defined DISTRIBUTE
                CASE (io_pio)
                  CALL pio_netcdf_get_fvar (ng, itlm, dav(ng)%name,     &
     &                                      'NLmodel_value', nlmodval)
                  IF (founderror(exit_flag, noerror,                    &
     &                           __line__, myfile)) RETURN
#  endif
              END SELECT
            END DO
          END IF
          IF (inner.ne.ninner) THEN
            linner=.true.
          ELSE
            linner=.false.
          END IF
# endif
          IF (master) THEN
            WRITE (stdout,60) 'Adjoint of', uppercase('rbl4dvar'),      &
     &                        outer, inner
          END IF
# ifdef RPCG
!
          inner_compute : IF (linner) THEN
!
# else
!
!  Call adjoint conjugate gradient algorithm.
!
          lcgini=.false.
          DO ng=1,ngrids
            CALL ad_congrad (ng, itlm, outer, inner, ninner, lcgini)
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
          END DO
# endif
!
!  Initialize the adjoint model from rest.
!
          DO ng=1,ngrids
            lsen4dvar(ng)=.false.
            lsenfct(ng)=.true.
# ifdef OBS_SPACE
            lobspace(ng)=.true.
#  ifndef OBS_IMPACT
            ladjvar(ng)=.true.
#  endif
# endif
            CALL close_file (ng, iadm, obs(ng), obs(ng)%name)
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
            CALL ad_initial (ng)
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
            wrtmisfit(ng)=.false.
          END DO
!
!  Set adjoint history NetCDF parameters.  Define adjoint history
!  file only once to avoid opening too many files.
!
          DO ng=1,ngrids
            wrtforce(ng)=.true.
            IF (adm(ng)%ncid.ne.-1) ldefadj(ng)=.false.
            fcount=adm(ng)%load
            adm(ng)%Nrec(fcount)=0
            adm(ng)%Rindex=0
          END DO
!
!  Time-step adjoint model backwards forced with current PSI vector.
!
          DO ng=1,ngrids
            IF (master) THEN
              WRITE (stdout,20) '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
!
!  Write out last weak-constraint forcing (WRTforce is still .TRUE.)
!  record into the adjoint history file.  Note that the weak-constraint
!  forcing is delayed by nADJ time-steps.
!
          DO ng=1,ngrids
# ifdef DISTRIBUTE
            CALL ad_wrt_his (ng, myrank)
# else
            CALL ad_wrt_his (ng, -1)
# endif
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
          END DO
!
!  Write out adjoint initial condition record into the adjoint
!  history file.
!
          DO ng=1,ngrids
            wrtforce(ng)=.false.
            lwrtstate2d(ng)=.false.
# ifdef DISTRIBUTE
            CALL ad_wrt_his (ng, myrank)
# else
            CALL ad_wrt_his (ng, -1)
# endif
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
          END DO
!
!  Convolve adjoint trajectory with error covariances.
!
          lposterior=.false.
          CALL error_covariance (itlm, driver, outer, inner,            &
     &                           lbck, lini, lold, lnew,                &
     &                           rec1, rec2, lposterior)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Convert the current adjoint solution in ADM(ng)%name to impulse
!  forcing. Write out impulse forcing into TLF(ng)%name NetCDF file.
!  To facilitate the forcing to the TLM and RPM, the forcing is
!  processed and written in increasing time coordinates (recall that
!  the adjoint solution in ADM(ng)%name is backwards in time).
!
!!
!! AMM: Do not know what to do in the weak constraint case yet.
!!
!!        IF (Master) THEN
!!          WRITE (stdout,40) outer, inner
!!        END IF
!!        DO ng=1,Ngrids
!!          TLF(ng)%Rindex=0
# ifdef DISTRIBUTE
!!          CALL wrt_impulse (ng, MyRank, iADM, ADM(ng)%name)
# else
!!          CALL wrt_impulse (ng, -1, iADM, ADM(ng)%name)
# endif
!!          IF (FoundError(exit_flag, NoError, __LINE__, MyFile)) RETURN
!!        END DO
!
!:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
!  Integrate tangent linear model.
!:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
!
!  Initialize tangent linear model from initial impulse which is now
!  stored in file ITL(ng)%name.
!
          DO ng=1,ngrids
            wrtnlmod(ng)=.false.
            wrttlmod(ng)=.true.
          END DO
!
!  If weak constraint, the impulses are time-interpolated at each
!  time-steps.
!
          DO ng=1,ngrids
            IF (frcrec(ng).gt.3) THEN
              frequentimpulse(ng)=.true.
            END IF
          END DO
!
!  Initialize tangent linear model from ITL(ng)%name, record 1.
!
          DO ng=1,ngrids
            itl(ng)%Rindex=rec1
            CALL close_file (ng, itlm, obs(ng), obs(ng)%name)
!
            CALL tl_initial (ng)
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
          END DO
!
!  Set tangent linear history NetCDF parameters.  Define tangent linear
!  history file at the beggining of each inner loop  to avoid opening
!  too many NetCDF files.
!
          DO ng=1,ngrids
            IF (inner.gt.ninner) ldeftlm(ng)=.false.
            fcount=tlm(ng)%load
            tlm(ng)%Nrec(fcount)=0
            tlm(ng)%Rindex=0
          END DO
!
!  Run tangent linear model forward and force with convolved adjoint
!  trajectory impulses.
!
          DO ng=1,ngrids
            IF (master) THEN
              WRITE (stdout,20) '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
 
          DO ng=1,ngrids
            wrtnlmod(ng)=.false.
            wrttlmod(ng)=.false.
          END DO
# ifdef RPCG
          END IF inner_compute
!
          DO ng=1,ngrids
            CALL ad_rpcg_lanczos (ng, irpm, outer, inner, ninner,       &
     &                            lcgini)
          END DO
# endif
 
        END DO ad_inner_loop
# ifndef RPCG
!
!  Call adjoint conjugate gradient algorithm.
!
        inner=0
        lcgini=.true.
        DO ng=1,ngrids
          CALL ad_congrad (ng, itlm, outer, inner, ninner, lcgini)
        END DO
# endif
 
#endif /* !OBS_IMPACT */
 
#ifdef OBS_IMPACT
!
!  Write out total observation impact.
!
        sourcefile=myfile
        DO ng=1,ngrids
          SELECT CASE (dav(ng)%IOtype)
            CASE (io_nf90)
              CALL netcdf_put_fvar (ng, inlm, dav(ng)%name,             &
     &                              'ObsImpact_total', ad_obsval,       &
# ifdef IMPACT_INNER
     &                              (/1,1/), (/mobs,ninner/),           &
# else
     &                              (/1/), (/mobs/),                    &
# endif
     &                              ncid = dav(ng)%ncid)
              IF (founderror(exit_flag, noerror,                        &
     &                       __line__, myfile)) RETURN
!
              CALL netcdf_sync (ng, inlm, dav(ng)%name, dav(ng)%ncid)
              IF (founderror(exit_flag, noerror,                        &
     &                       __line__, myfile)) RETURN
 
# if defined PIO_LIB && defined DISTRIBUTE
            CASE (io_pio)
              CALL pio_netcdf_put_fvar (ng, inlm, dav(ng)%name,         &
     &                                  'ObsImpact_total', ad_obsval,   &
#  ifdef IMPACT_INNER
     &                                  (/1,1/), (/mobs,ninner/),       &
#  else
     &                                  (/1/), (/mobs/),                &
#  endif
     &                                  piofile = dav(ng)%pioFile)
              IF (founderror(exit_flag, noerror,                        &
     &                       __line__, myfile)) RETURN
!
              CALL pio_netcdf_sync (ng, inlm, dav(ng)%name,             &
     &                              dav(ng)%pioFile)
              IF (founderror(exit_flag, noerror,                        &
     &                       __line__, myfile)) RETURN
# endif
          END SELECT
        END DO
#else
!
!  Write out observation sensitivity.
!
        sourcefile=myfile
        DO ng=1,ngrids
          SELECT CASE (dav(ng)%IOtype)
            CASE (io_nf90)
              CALL netcdf_put_fvar (ng, itlm, dav(ng)%name,             &
     &                              'ObsSens_total', ad_obsval,         &
# ifdef IMPACT_INNER
     &                              (/1,1/), (/mobs,ninner/),           &
# else
     &                              (/1/), (/mobs/),                    &
# endif
     &                              ncid = dav(ng)%ncid)
              IF (founderror(exit_flag, noerror,                        &
     &                       __line__, myfile)) RETURN
!
              CALL netcdf_sync (ng, inlm, dav(ng)%name, dav(ng)%ncid)
              IF (founderror(exit_flag, noerror,                        &
     &                       __line__, myfile)) RETURN
 
# if defined PIO_LIB && defined DISTRIBUTE
            CASE (io_pio)
              CALL pio_netcdf_put_fvar (ng, itlm, dav(ng)%name,         &
     &                                  'ObsSens_total', ad_obsval,     &
#  ifdef IMPACT_INNER
     &                                  (/1,1/), (/mobs,ninner/),       &
#  else
     &                                  (/1/), (/mobs/),                &
#  endif
     &                                  piofile = dav(ng)%pioFile)
              IF (founderror(exit_flag, noerror,                        &
     &                       __line__, myfile)) RETURN
!
              CALL pio_netcdf_sync (ng, inlm, dav(ng)%name,             &
     &                              dav(ng)%pioFile)
              IF (founderror(exit_flag, noerror,                        &
     &                       __line__, myfile)) RETURN
# endif
          END SELECT
        END DO
#endif
!
!  Close tangent linear NetCDF file.
!
        sourcefile=myfile
        DO ng=1,ngrids
          CALL close_file (ng, itlm, tlm(ng), tlm(ng)%name)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END DO
 
#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.
!:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
!
        DO ng=1,ngrids
          wrtnlmod(ng)=.false.
          wrttlmod(ng)=.true.
          lwrttlm(ng)=.false.
        END DO
!
!  Clear tangent linear forcing arrays before entering inner-loop.
!  This is very important since these arrays are non-zero after
!  running the representer model and must be zero when running the
!  tangent linear model.
!
        DO ng=1,ngrids
          DO tile=first_tile(ng),last_tile(ng),+1
            CALL initialize_forces (ng, tile, itlm)
# ifdef ADJUST_BOUNDARY
            CALL initialize_boundary (ng, tile, itlm)
# endif
          END DO
        END DO
!
!  Set basic state trajectory.
!
        DO ng=1,ngrids
          WRITE (fwd(ng)%name,10) trim(his(ng)%head), outer-1
          lstr=len_trim(fwd(ng)%name)
          fwd(ng)%base=fwd(ng)%name(1:lstr-3)
        END DO
!
!  If weak constraint, the impulses are time-interpolated at each
!  time-steps.
!
        DO ng=1,ngrids
          IF (frcrec(ng).gt.3) THEN
            frequentimpulse(ng)=.true.
          END IF
        END DO
!
!  Initialize tangent linear model from ITL(ng)%name, record 1.
!
        DO ng=1,ngrids
          itl(ng)%Rindex=rec1
          CALL close_file (ng, itlm, obs(ng), obs(ng)%name)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
          CALL tl_initial (ng)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END DO
!
!  Clear tangent linear forcing arrays and boundary arrays
!  before the obs impact initial condition calculation.
!
        DO ng=1,ngrids
          DO tile=first_tile(ng),last_tile(ng),+1
            CALL initialize_forces (ng, tile, itlm)
# ifdef ADJUST_BOUNDARY
            CALL initialize_boundary (ng, tile, itlm)
# endif
          END DO
        END DO
!
!  Run tangent linear model forward and force with convolved adjoint
!  trajectory impulses. Compute (HMBM'H')_n * PSI at observation points
!  which are used in the conjugate gradient algorithm.
!
        DO ng=1,ngrids
          IF (master) THEN
            WRITE (stdout,20) '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
!
        DO ng=1,ngrids
          wrtnlmod(ng)=.false.
          wrttlmod(ng)=.false.
        END DO
!
!  Compute observation impact to the data assimilation system.
!
        DO ng=1,ngrids
# ifdef RPCG
          CALL rep_matrix (ng, itlm, outer, ninner-1)
# else
          CALL rep_matrix (ng, itlm, outer, ninner)
# endif
        END DO
!
!  Write out observation sentivity.
!
        sourcefile=myfile
        DO ng=1,ngrids
          SELECT CASE (dav(ng)%IOtype)
            CASE (io_nf90)
              CALL netcdf_put_fvar (ng, itlm, dav(ng)%name,             &
     &                              'ObsImpact_IC', ad_obsval,          &
# ifdef IMPACT_INNER
     &                              (/1,1/), (/mobs,ninner/),           &
# else
     &                              (/1/), (/mobs/),                    &
# endif
     &                              ncid = dav(ng)%ncid)
              IF (founderror(exit_flag, noerror,                        &
     &                       __line__, myfile)) RETURN
!
              CALL netcdf_sync (ng, inlm, dav(ng)%name, dav(ng)%ncid)
              IF (founderror(exit_flag, noerror,                        &
     &                       __line__, myfile)) RETURN
 
# if defined PIO_LIB && defined DISTRIBUTE
            CASE (io_pio)
              CALL pio_netcdf_put_fvar (ng, itlm, dav(ng)%name,         &
     &                                  'ObsImpact_IC', ad_obsval,      &
# ifdef IMPACT_INNER
     &                                  (/1,1/), (/mobs,ninner/),       &
# else
     &                                  (/1/), (/mobs/),                &
# endif
     &                                  piofile = dav(ng)%pioFile)
              IF (founderror(exit_flag, noerror,                        &
     &                       __line__, myfile)) RETURN
!
              CALL pio_netcdf_sync (ng, inlm, dav(ng)%name,             &
     &                              dav(ng)%pioFile)
              IF (founderror(exit_flag, noerror,                        &
     &                       __line__, myfile)) RETURN
# endif
          END SELECT
        END DO
 
# 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.
!:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
!
        DO ng=1,ngrids
          wrtnlmod(ng)=.false.
          wrttlmod(ng)=.true.
          lwrttlm(ng)=.false.
        END DO
!
!  Clear tangent linear forcing arrays before entering inner-loop.
!  This is very important since these arrays are non-zero after
!  running the representer model and must be zero when running the
!  tangent linear model.
!
        DO ng=1,ngrids
          DO tile=first_tile(ng),last_tile(ng),+1
            CALL initialize_forces (ng, tile, itlm)
#  ifdef ADJUST_BOUNDARY
            CALL initialize_boundary (ng, tile, itlm)
#  endif
          END DO
        END DO
!
!  Set basic state trajectory.
!
        DO ng=1,ngrids
          WRITE (fwd(ng)%name,10) trim(his(ng)%head), outer-1
          lstr=len_trim(fwd(ng)%name)
          fwd(ng)%base=fwd(ng)%name(1:lstr-3)
        END DO
!
!  If weak constraint, the impulses are time-interpolated at each
!  time-steps.
!
        DO ng=1,ngrids
          IF (frcrec(ng).gt.3) THEN
            frequentimpulse(ng)=.true.
          END IF
        END DO
!
!  Initialize tangent linear model from ITL(ng)%name, record 1.
!
        DO ng=1,ngrids
          itl(ng)%Rindex=rec1
          CALL close_file (ng, itlm, obs(ng), obs(ng)%name)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
          CALL tl_initial (ng)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END DO
!
!  Clear tangent initial condition arrays and boundary arrays
!  before the obs impact initial condition calculation.
!
        DO ng=1,ngrids
          DO tile=first_tile(ng),last_tile(ng),+1
            CALL initialize_ocean (ng, tile, itlm)
#  ifdef ADJUST_BOUNDARY
            CALL initialize_boundary (ng, tile, itlm)
#  endif
          END DO
        END DO
!
!  Run tangent linear model forward and force with convolved adjoint
!  trajectory impulses. Compute (HMBM'H')_n * PSI at observation points
!  which are used in the conjugate gradient algorithm.
!
        DO ng=1,ngrids
          IF (master) THEN
            WRITE (stdout,20) '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
!
        DO ng=1,ngrids
          wrtnlmod(ng)=.false.
          wrttlmod(ng)=.false.
        END DO
!
!  Compute observation impact to the data assimilation system.
!
        DO ng=1,ngrids
#  ifdef RPCG
          CALL rep_matrix (ng, itlm, outer, ninner-1)
#  else
          CALL rep_matrix (ng, itlm, outer, ninner)
#  endif
        END DO
!
!  Write out observation sentivity.
!
        sourcefile=myfile
        DO ng=1,ngrids
          SELECT CASE (dav(ng)%IOtype)
            CASE (io_nf90)
              CALL netcdf_put_fvar (ng, itlm, dav(ng)%name,             &
     &                              'ObsImpact_FC', ad_obsval,          &
#  ifdef IMPACT_INNER
     &                              (/1,1/), (/mobs,ninner/),           &
#  else
     &                              (/1/), (/mobs/),                    &
#  endif
     &                              ncid = dav(ng)%ncid)
              IF (founderror(exit_flag, noerror,                        &
     &                       __line__, myfile)) RETURN
!
              CALL netcdf_sync (ng, inlm, dav(ng)%name, dav(ng)%ncid)
              IF (founderror(exit_flag, noerror,                        &
     &                       __line__, myfile)) RETURN
 
#  if defined PIO_LIB && defined DISTRIBUTE
            CASE (io_pio)
              CALL pio_netcdf_put_fvar (ng, itlm, dav(ng)%name,         &
     &                                  'ObsImpact_FC', ad_obsval,      &
#   ifdef IMPACT_INNER
     &                                  (/1,1/), (/mobs,ninner/),       &
#   else
     &                                  (/1/), (/mobs/),                &
#   endif
     &                                  piofile = dav(ng)%pioFile)
              IF (founderror(exit_flag, noerror,                        &
     &                       __line__, myfile)) RETURN
!
              CALL pio_netcdf_sync (ng, inlm, dav(ng)%name,             &
     &                              dav(ng)%pioFile)
              IF (founderror(exit_flag, noerror,                        &
     &                       __line__, myfile)) RETURN
#  endif
          END SELECT
        END DO
# endif
 
# if defined ADJUST_BOUNDARY
!
!:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
!  Integrate tangent linear model with boundary condition increments
!  only to compute the observation impact associated with the boundary
!  conditions.
!:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
!
        DO ng=1,ngrids
          wrtnlmod(ng)=.false.
          wrttlmod(ng)=.true.
          lwrttlm(ng)=.false.
        END DO
!
!  Clear tangent linear forcing arrays before entering inner-loop.
!  This is very important since these arrays are non-zero after
!  running the representer model and must be zero when running the
!  tangent linear model.
!
        DO ng=1,ngrids
          DO tile=first_tile(ng),last_tile(ng),+1
            CALL initialize_forces (ng, tile, itlm)
#  ifdef ADJUST_BOUNDARY
            CALL initialize_boundary (ng, tile, itlm)
#  endif
          END DO
        END DO
!
!  Set basic state trajectory.
!
        DO ng=1,ngrids
          WRITE (fwd(ng)%name,10) trim(his(ng)%head), outer-1
          lstr=len_trim(fwd(ng)%name)
          fwd(ng)%base=fwd(ng)%name(1:lstr-3)
        END DO
!
!  If weak constraint, the impulses are time-interpolated at each
!  time-steps.
!
        DO ng=1,ngrids
          IF (frcrec(ng).gt.3) THEN
            frequentimpulse(ng)=.true.
          END IF
        END DO
!
!  Initialize tangent linear model from ITL(ng)%name, record Rec1.
!
        DO ng=1,ngrids
          itl(ng)%Rindex=rec1
          CALL close_file (ng, itlm, obs(ng), obs(ng)%name)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
          CALL tl_initial (ng)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END DO
!
!  Clear tangent linear initial condition and forcing arrays
!  before the obs impact initial condition calculation.
!
        DO ng=1,ngrids
          DO tile=first_tile(ng),last_tile(ng),+1
            CALL initialize_ocean (ng, tile, itlm)
            CALL initialize_forces (ng, tile, itlm)
          END DO
        END DO
!
!  Run tangent linear model forward and force with convolved adjoint
!  trajectory impulses. Compute (HMBM'H')_n * PSI at observation points
!  which are used in the conjugate gradient algorithm.
!
        DO ng=1,ngrids
          IF (master) THEN
            WRITE (stdout,20) '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
!
        DO ng=1,ngrids
          wrtnlmod(ng)=.false.
          wrttlmod(ng)=.false.
        END DO
!
!  Compute observation impact to the data assimilation system.
!
        DO ng=1,ngrids
#  ifdef RPCG
          CALL rep_matrix (ng, itlm, outer, ninner-1)
#  else
          CALL rep_matrix (ng, itlm, outer, ninner)
#  endif
        END DO
!
!  Write out observation sentivity.
!
        sourcefile=myfile
        DO ng=1,ngrids
          SELECT CASE (dav(ng)%IOtype)
            CASE (io_nf90)
              CALL netcdf_put_fvar (ng, itlm, dav(ng)%name,             &
     &                              'ObsImpact_BC', ad_obsval,          &
#  ifdef IMPACT_INNER
     &                              (/1,1/), (/mobs,ninner/),           &
#  else
     &                              (/1/), (/mobs/),                    &
#  endif
     &                              ncid = dav(ng)%ncid)
              IF (founderror(exit_flag, noerror,                        &
     &                       __line__, myfile)) RETURN
!
              CALL netcdf_sync (ng, inlm, dav(ng)%name, dav(ng)%ncid)
              IF (founderror(exit_flag, noerror,
     &                       __line__, myfile)) RETURN
 
#  if defined PIO_LIB && defined DISTRIBUTE
            CASE (io_pio)
              CALL pio_netcdf_put_fvar (ng, itlm, dav(ng)%name,         &
     &                                  'ObsImpact_BC', ad_obsval,      &
#   ifdef IMPACT_INNER
     &                                  (/1,1/), (/mobs,ninner/),       &
#   else
     &                                  (/1/), (/mobs/),                &
#   endif
     &                                  piofile = dav(ng)%pioFile)
              IF (founderror(exit_flag, noerror,                        &
     &                       __line__, myfile)) RETURN
!
              CALL pio_netcdf_sync (ng, inlm, dav(ng)%name,             &
     &                              dav(ng)%pioFile)
              IF (founderror(exit_flag, noerror,
     &                       __line__, myfile)) RETURN
#  endif
          END SELECT
        END DO
# endif
#endif /* OBS_IMPACT_SPLIT */
!
!  Close current forward NetCDF file.
!
        sourcefile=myfile
        DO ng=1,ngrids
          CALL close_file (ng, inlm, fwd(ng), fwd(ng)%name)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
          IF (his(ng)%IOtype.eq.io_nf90) THEN
            his(ng)%ncid=-1
#if defined PIO_LIB && defined DISTRIBUTE
          ELSE IF (his(ng)%IOtype.eq.io_pio) THEN
            his(ng)%pioFile%fh=-1
#endif
          END IF
        END DO
 
      END DO ad_outer_loop
!
 10   FORMAT (a,'_outer',i0,'.nc')
 20   FORMAT (/,1x,a,1x,'ROMS: started time-stepping:',                 &
     &        ' (Grid: ',i2.2,' TimeSteps: ',i8.8,' - ',i8.8,')',/)
 30   FORMAT (' (',i3.3,',',i3.3,'): ',a,' data penalty, Jdata = ',     &
     &        1p,e17.10,0p,t68,a)
 40   FORMAT (/,' Converting Convolved Adjoint Trajectory to',          &
     &          ' Impulses: Outer = ',i3.3,' Inner = ',i3.3,/)
 50   FORMAT (/,'ROMS: Started adjoint Sensitivity calculation',        &
     &          ' ...',/)
 60   FORMAT (/,'ROMS: ',a,1x,a,', Outer = ',i3.3,                      &
     &          ' Inner = ',i3.3,/)
 70   FORMAT (/,1x,a,1x,'ROMS: adjoint forcing time range: ',           &
     &        f12.4,' - ',f12.4 ,/)
 80   FORMAT (a,'_',a,'.nc')
 90   FORMAT (a,'.nc')
!
      RETURN

References ad_congrad(), ad_initial(), ad_main2d(), ad_main3d(), mod_fourdvar::ad_obsval, ad_wrt_his_mod::ad_wrt_his(), mod_iounits::adm, mod_iounits::ads, mod_scalars::balance, zeta_balance_mod::balance_ref(), zeta_balance_mod::biconj(), mod_scalars::blowup, mod_fourdvar::cg_beta, mod_fourdvar::cg_delta, mod_fourdvar::cg_dla, mod_fourdvar::cg_gnorm_v, mod_fourdvar::cg_qg, close_io_mod::close_file(), close_io_mod::close_inp(), close_io_mod::close_out(), mod_iounits::dav, mod_fourdvar::deallocate_fourdvar(), def_impulse_mod::def_impulse(), def_mod_mod::def_mod(), def_norm_mod::def_norm(), mod_scalars::dends, mod_scalars::dstart, mod_scalars::dstarts, mod_scalars::dstrs, mod_scalars::dt, edit_multifile(), mod_scalars::erend, convolve_mod::error_covariance(), mod_scalars::erstr, mod_scalars::exit_flag, mod_parallel::first_tile, strings_mod::founderror(), mod_scalars::frcrec, mod_scalars::frequentimpulse, mod_iounits::fwd, get_state_mod::get_state(), mod_iounits::his, mod_param::iadm, mod_ncparam::idobss, mod_iounits::ini, initial(), mod_boundary::initialize_boundary(), mod_forces::initialize_forces(), mod_fourdvar::initialize_fourdvar(), mod_ocean::initialize_ocean(), mod_param::inlm, mod_scalars::inner, mod_ncparam::io_nf90, mod_ncparam::io_pio, mod_param::irpm, mod_ncparam::isfsur, mod_iounits::itl, mod_param::itlm, mod_fourdvar::jb0, mod_parallel::last_tile, mod_stepping::lbinp, mod_stepping::lbout, mod_scalars::lcyclerst, mod_iounits::lcz, mod_scalars::ldefadj, mod_scalars::ldefitl, mod_scalars::ldefmod, mod_scalars::ldefnrm, mod_scalars::ldeftlf, mod_scalars::ldeftlm, mod_stepping::lfinp, mod_stepping::lfout, mod_stepping::lnew, ini_adjust_mod::load_adtotl(), ini_adjust_mod::load_tltoad(), mod_fourdvar::lobspace, mod_stepping::lold, mod_scalars::lreadblk, mod_scalars::lreadfrc, mod_scalars::lreadfwd, mod_scalars::lreadqck, mod_fourdvar::lsen4dvar, mod_fourdvar::lsenfct, mod_scalars::lstiffness, mod_scalars::lwrtnrm, mod_scalars::lwrtrst, mod_scalars::lwrtstate2d, mod_scalars::lwrttlm, mod_parallel::master, memory(), mod_fourdvar::mobs, mod_parallel::myrank, mod_netcdf::netcdf_sync(), mod_param::ngrids, mod_fourdvar::nimpact, mod_scalars::ninner, mod_fourdvar::nlmodval, mod_scalars::noerror, normalization_mod::normalization(), mod_scalars::nouter, mod_iounits::nrm, mod_scalars::nrun, mod_param::nsa, mod_scalars::ntend, mod_scalars::ntimes, mod_fourdvar::ntimes_ana, mod_fourdvar::ntimes_fct, mod_scalars::ntstart, mod_iounits::obs, mod_fourdvar::obsscale, mod_iounits::oifa, mod_iounits::oifb, mod_scalars::outer, mod_pio_netcdf::pio_netcdf_sync(), rep_matrix(), mod_iounits::rst, mod_scalars::sec2day, mod_iounits::sourcefile, stats_modobs_mod::stats_modobs(), mod_iounits::stdout, mod_scalars::tdays, tl_def_ini_mod::tl_def_ini(), tl_initial(), tl_main2d(), tl_main3d(), mod_iounits::tlm, strings_mod::uppercase(), mod_fourdvar::vcglwk, mod_ncparam::vname, wclock_off(), wrt_ini_mod::wrt_ini(), wrt_rst_mod::wrt_rst(), mod_fourdvar::wrtforce, mod_fourdvar::wrtmisfit, mod_fourdvar::wrtnlmod, mod_fourdvar::wrtobsscale, mod_fourdvar::wrtrpmod, mod_fourdvar::wrttlmod, mod_fourdvar::wrtzetaref, mod_fourdvar::zcglwk, and mod_fourdvar::zgrad0.

Here is the call graph for this function:

roms_run() [3/3]

subroutine, public roms_kernel_mod::roms_run	(	real(dp), intent(in)	runinterval,
		integer, intent(in), optional	kernel )

Definition at line 336 of file jedi_roms.h.

!
!=======================================================================
!                                                                      !
!  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

References ad_main3d(), mod_scalars::blowup, close_io_mod::close_inp(), close_io_mod::close_out(), mod_scalars::dt, mod_scalars::exit_flag, strings_mod::founderror(), mod_param::iadm, mod_scalars::iic, mod_scalars::initime, mod_param::inlm, mod_param::itlm, mod_scalars::lcyclerst, mod_scalars::lwrtrst, main3d(), mod_parallel::master, memory(), mod_parallel::myrank, mod_scalars::myruninterval, mod_param::ngrids, mod_scalars::noerror, mod_scalars::ntimes, mod_scalars::ntstart, mod_scalars::processinputdata, mod_iounits::rst, mod_iounits::sourcefile, mod_iounits::stdout, mod_scalars::step_counter, mod_scalars::time, tl_main3d(), wclock_off(), and wrt_rst_mod::wrt_rst().

Here is the call graph for this function:

tlm_initial()

subroutine, private roms_kernel_mod::tlm_initial ( integer, intent(in) phase )

private

Definition at line 999 of file jedi_roms.h.

!
!=======================================================================
!                                                                      !
!  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

References analytical_mod::ana_drag(), analytical_mod::ana_sponge(), mod_scalars::avgke, mod_scalars::avgkp, mod_scalars::avgpe, check_multifile(), close_io_mod::close_inp(), mod_scalars::day2sec, mod_scalars::dstart, mod_scalars::dt, edit_multifile(), mod_scalars::erstr, mod_scalars::exit_flag, mod_parallel::first_tile, mod_scalars::first_time, strings_mod::founderror(), grid_coords(), mod_scalars::iic, mod_scalars::iif, mod_scalars::indx1, ini_hmixcoef_mod::ini_hmixcoef(), mod_scalars::initime, mod_param::itlm, mod_scalars::jic, mod_stepping::knew, mod_stepping::krhs, mod_stepping::kstp, mod_parallel::last_tile, mod_scalars::ldeftlm, mod_scalars::lreadblk, mod_scalars::lreadfrc, mod_scalars::lreadqck, mod_scalars::lsponge, mod_scalars::lstiffness, mod_scalars::lwrttlm, mod_parallel::master, mod_stepping::nf, mod_stepping::nfm1, mod_stepping::nfm2, mod_stepping::nfm3, mod_stepping::nfp1, mod_param::ngrids, mod_stepping::nnew, mod_scalars::noerror, mod_stepping::nrhs, mod_scalars::nrrec, mod_scalars::nrun, mod_stepping::nstp, mod_scalars::ntend, mod_scalars::ntfirst, mod_scalars::ntimes, mod_scalars::ntstart, mod_scalars::predictor_2d_step, mod_iounits::qck, set_masks_mod::set_masks(), mod_iounits::stdout, stiffness_mod::stiffness(), mod_scalars::synchro_flag, mod_scalars::tdays, mod_scalars::time, mod_scalars::time4jedi, mod_scalars::time_code, dateclock_mod::time_string(), tl_set_depth_mod::tl_bath(), tl_get_data(), tl_get_idata(), tl_post_initial_mod::tl_post_initial(), tl_set_data(), mod_scalars::tl_ubar_xs, mod_scalars::tl_volcons, mod_scalars::volume, wclock_off(), wclock_on(), and wetdry_mod::wetdry().

Referenced by roms_initializep1(), roms_initializep2(), and roms_initializep3().

Here is the call graph for this function:

Here is the caller graph for this function: