split_r4dvar_roms.h

Go to the documentation of this file.

      MODULE roms_kernel_mod
!
!git $Id$
!=================================================== Andrew M. Moore ===
!  Copyright (c) 2002-2025 The ROMS Group           Hernan G. Arango   !
!    Licensed under a MIT/X style license                              !
!    See License_ROMS.md                                               !
!=======================================================================
!                                                                      !
!  ROMS Strong/Weak Constraint Split 4-Dimensional Variational Data    !
!            Assimilation Driver: Indirect Representer Approach        !
!                                 (R4D-Var).                           !
!                                                                      !
!  This driver is used for the dual formulation (observation space),   !
!  strong or weak constraint 4D-Var where errors may be considered     !
!  in both model and observations.                                     !
!                                                                      !
!  The R4D-Var algorithm is split into multiple executables to         !
!  facilitate various configurations:                                  !
!                                                                      !
!    (1) Executable A computes ROMS nonlinear trajectory used to       !
!        linearize the tangent linear and adjoint models used in       !
!        the iterations of the inner loop for the minimization of      !
!        the cost function. It allows the nonlinear trajectory to      !
!        be part of a coupling system and or include nested grids.     !
!        It calls either the RBL4D-Var "background" or "analysis"      !
!        routines.                                                     !
!                                                                      !
!    (2) Executable B calls either R4D-Var "increment" or              !
!        "posterior_error". The R4D-Var increment is obtained by       !
!        minimizing the cost function over Ninner loops. It is         !
!        possible to use a coarser grid resolution in the inner        !
!        loop.  If so, the finer background trajectory needs to        !
!        be interpolated into the coarser grid. Then, at the end       !
!        of inner loops, the coarse grid increment needs to be         !
!        interpolated to the finer grid.  The increment phase          !
!        may be run at a lower precision.                              !
!                                                                      !
!  The routines in this driver control the initialization,  time-      !
!  stepping, and finalization of ROMS model following ESMF/NUOPC       !
!  conventions:                                                        !
!                                                                      !
!     ROMS_initialize                                                  !
!     ROMS_run                                                         !
!     ROMS_finalize                                                    !
!                                                                      !
!  References:                                                         !
!                                                                      !
!    Moore, A.M., H.G. Arango, G. Broquet, B.S. Powell, A.T. Weaver,   !
!      and J. Zavala-Garay, 2011: The Regional Ocean Modeling System   !
!      (ROMS)  4-dimensional variational data assimilations systems,   !
!      Part I - System overview and formulation, Prog. Oceanogr., 91,  !
!      34-49, doi:10.1016/j.pocean.2011.05.004.                        !
!                                                                      !
!    Moore, A.M., H.G. Arango, G. Broquet, C. Edward, M. Veneziani,    !
!      B. Powell, D. Foley, J.D. Doyle, D. Costa, and P. Robinson,     !
!      2011: The Regional Ocean Modeling System (ROMS) 4-dimensional   !
!      variational data assimilations systems, Part II - Performance   !
!      and application to the California Current System, Prog.         !
!      Oceanogr., 91, 50-73, doi:10.1016/j.pocean.2011.05.003.         !
!                                                                      !
!=======================================================================
!
      USE mod_param
      USE mod_parallel
      USE mod_arrays
      USE mod_fourdvar
      USE mod_iounits
      USE mod_ncparam
      USE mod_netcdf
#if defined PIO_LIB && defined DISTRIBUTE
      USE mod_pio_netcdf
#endif
      USE mod_scalars
      USE mod_stepping
!
      USE r4dvar_mod
!
      USE close_io_mod,      ONLY : close_inp, close_out
      USE def_dai_mod,       ONLY : def_dai
      USE get_state_mod,     ONLY : get_state
      USE inp_par_mod,       ONLY : inp_par
#ifdef MCT_LIB
# ifdef ATM_COUPLING
      USE mct_coupler_mod,   ONLY : initialize_ocn2atm_coupling
# endif
# ifdef WAV_COUPLING
      USE mct_coupler_mod,   ONLY : initialize_ocn2wav_coupling
# endif
#endif
      USE stats_modobs_mod,  ONLY : stats_modobs
      USE stdinp_mod,        ONLY : getpar_s
      USE stdout_mod,        ONLY : set_stdoutunit, stdout_unit
      USE strings_mod,       ONLY : founderror, uppercase
      USE wrt_dai_mod,       ONLY : wrt_dai
      USE wrt_rst_mod,       ONLY : wrt_rst
!
      implicit none
!
      PUBLIC :: roms_initialize
      PUBLIC :: roms_run
      PUBLIC :: roms_finalize
!
      CONTAINS
!
      SUBROUTINE roms_initialize (first, mpiCOMM)
!
!=======================================================================
!                                                                      !
!  This routine allocates and initializes ROMS state variables and     !
!  internal parameters. It reads standard input 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
!
!  Get 4D-Var phase from APARNAM input script file.
!
        CALL getpar_s (master, aparnam, 'APARNAM')
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
        CALL getpar_s (master, phase4dvar, 'Phase4DVAR', aparnam)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Determine ROMS standard output append switch. It is only relevant if
!  "ROMS_STDINP" is activated. The standard output is created in the
!  "background" phase and open to append in the other phases. Set
!  switch so the "stiffness" routine is only called in the "background"
!  phase.
!
        IF (index(trim(uppercase(phase4dvar)),'BACKG').ne.0) THEN
          lappend=.false.
          lstiffness=.true.
        ELSE
          lappend=.true.
          lstiffness=.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. The 'Nrun' counter will be recomputed in the
!  RBL4D-Var phases to process the obervation operator correctly.
!
        nrun=1                ! run counter
        erstr=1               ! ensemble start counter
        erend=nouter          ! ensemble end 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 _OPENMP
        mythread=my_threadnum()
#elif defined DISTRIBUTE
        mythread=myrank
#else
        mythread=0
#endif
        DO ng=1,ngrids
          chunk_size=(ntilex(ng)*ntilee(ng)+numthreads-1)/numthreads
          first_tile(ng)=mythread*chunk_size
          last_tile(ng)=first_tile(ng)+chunk_size-1
        END DO
!
!  Initialize internal wall clocks. Notice that the timings does not
!  includes processing standard input because several parameters are
!  needed to allocate clock variables.
!
        IF (master) THEN
          WRITE (stdout,10)
 10       FORMAT (/,' Process Information:',/)
        END IF
!
        DO ng=1,ngrids
          DO thread=thread_range
            CALL wclock_on (ng, inlm, 0, __line__, myfile)
          END DO
        END DO
!
!  Allocate and initialize modules variables.
!
        CALL roms_allocate_arrays (allocate_vars)
        CALL roms_initialize_arrays
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
      END IF
 
#if defined MCT_LIB && (defined ATM_COUPLING || defined WAV_COUPLING)
!
!-----------------------------------------------------------------------
!  Initialize coupling streams between model(s).
!-----------------------------------------------------------------------
!
      DO ng=1,ngrids
# ifdef ATM_COUPLING
        CALL initialize_ocn2atm_coupling (ng, myrank)
# endif
# ifdef WAV_COUPLING
        CALL initialize_ocn2wav_coupling (ng, myrank)
# endif
      END DO
#endif
!
!-----------------------------------------------------------------------
!  Set application grid, metrics, and associated variables. Then,
!  proccess background and mode prior error covariance standard
!  deviations and normalization coefficients.
!-----------------------------------------------------------------------
!
      DO ng=1,ngrids
#ifdef STD_MODEL
        lwrtstd(ng)=.true.
        IF (index(trim(uppercase(phase4dvar)),'BACKG').ne.0) THEN
          ldefstd(ng)=.true.
          lreadstd(ng)=.false.
        ELSE
          ldefstd(ng)=.false.
          lreadstd(ng)=.true.
        END IF
#else
        lreadstd(ng)=.true.
#endif
        CALL prior_error (ng)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        setgridconfig(ng)=.false.
      END DO
!
      RETURN
      END SUBROUTINE roms_initialize
!
      SUBROUTINE roms_run (RunInterval)
!
!=======================================================================
!                                                                      !
!  This subroutine runs the Strong or Weak constraint, Indirect        !
!  Representers 4D-Var data assimilation (R4D-Var) algorithm. It       !
!  time-steps ROMS nonlinear, representer, tangent linear, and         !
!  adjoint kernels.                                                    !
!                                                                      !
!  On Input:                                                           !
!                                                                      !
!     RunInterval     Execution time stepping window (seconds)         !
!                                                                      !
!=======================================================================
!
!  Imported variable declarations
!
      real(dp), intent(in) :: RunInterval
!
!  Local variable declarations.
!
      integer :: my_outer, ng
!
      character (len=*), parameter :: MyFile =                          &
     &  __FILE__//", ROMS_run"
!
!=======================================================================
!  Run Split R4D-Var Data Assimilation algorithm.
!=======================================================================
!
!  Initialize several global 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
!
      ldone=.false.         ! 4D-Var cycle finish switch
!
!  Select R4D-Var phase to execute.
!
      SELECT CASE (uppercase(phase4dvar(1:6)))
 
!
!  Compute nonlinear background state trajectory, Xb(t)|n-1. Interpolate
!  the background at the observation locations, and compute the quality
!  control accept/reject flag, ObsScale. The background state is used
!  to linearize the tangent linear and adjoint models during the
!  minimization.
!
        CASE ('BACKGR')
 
          my_outer=0
          outer=0
          inner=0
 
          CALL background (outer, runinterval)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Compute 4D-Var data assimilation increment, dXa, by iterating over
!  the inner loops, and minimizing the cost function.
!
        CASE ('INCREM')
 
          my_outer=outerloop
          outer=outerloop
          inner=0
 
          CALL increment (my_outer, runinterval)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Compute 4D-Var data assimilation analysis, Xa = Xb + dXa.  Set
!  nonlinear model initial conditions for next outer loop.
!
        CASE ('ANALYS')
 
          my_outer=outerloop
          outer=outerloop
          inner=ninner
 
          CALL analysis (my_outer, runinterval)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
 
#if defined POSTERIOR_ERROR_I || \
    defined posterior_error_f || \
    defined posterior_eofs
!
!  Compute full (diagonal) posterior analysis error covariance matrix.
!  (NOTE: Currently, this code only works for a single outer-loop).
!
        CASE ('POST_E')
 
          CALL posterior_error (runinterval)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
#endif
!
!  Issue an error if incorrect 4D-Var phase.
!
        CASE DEFAULT
 
          IF (master) THEN
            WRITE (stdout,10) trim(phase4dvar)
 10         FORMAT (' ROMS_run - illegal 4D-Var phase: ''',a,'''')
          END IF
          exit_flag=5
          RETURN
 
      END SELECT
!
!  Set finish R4D-Var cycle switch.
!
      IF ((outerloop.eq.nouter).and.                                    &
#if defined POSTERIOR_ERROR_I || \
    defined posterior_error_f || \
    defined posterior_eofs
     &    (index(trim(uppercase(phase4dvar)),'POST_E').ne.0)) THEN
        ldone=.true.
#else
     &    (index(trim(uppercase(phase4dvar)),'ANALYS').ne.0)) THEN
        ldone=.true.
#endif
      END IF
!
      RETURN
      END SUBROUTINE roms_run
!
      SUBROUTINE roms_finalize
!
!=======================================================================
!                                                                      !
!  This routine terminates ROMS R4D-Var execution.                     !
!                                                                      !
!=======================================================================
!
!  Local variable declarations.
!
      integer :: Fcount, InpRec, Nfiles, Tindex
      integer :: ifile, lstr, ng, thread
!
      character (len=10) :: suffix
 
      character (len=*), parameter :: MyFile =                          &
     &  __FILE__//", ROMS_finalize"
!
!-----------------------------------------------------------------------
!  Write out 4D-Var analysis fields that can be used as the initial
!  conditions for the next data assimilation cycle. Here, use the
!  last record of the RPM for the final outer loop.
!-----------------------------------------------------------------------
!
      IF (ldone.and.(exit_flag.eq.noerror)) THEN
        DO ng=1,ngrids
          ldefdai(ng)=.true.
          CALL def_dai (ng)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
          WRITE (tlm(ng)%name,10) trim(fwd(ng)%head), nouter
 10       FORMAT (a,'_outer',i0,'.nc')
          lstr=len_trim(tlm(ng)%name)
          tlm(ng)%base=tlm(ng)%name(1:lstr-3)
          IF (tlm(ng)%Nfiles.gt.1) THEN
            nfiles=tlm(ng)%Nfiles
            DO ifile=1,nfiles
              WRITE (suffix,"('_',i4.4,'.nc')") ifile
              tlm(ng)%files(ifile)=trim(tlm(ng)%base)//trim(suffix)
            END DO
            tlm(ng)%name=trim(tlm(ng)%files(nfiles))
          ELSE
            tlm(ng)%files(1)=trim(tlm(ng)%name)
          END IF
!
          SELECT CASE (tlm(ng)%IOtype)
            CASE (io_nf90)
              CALL netcdf_get_dim (ng, irpm, tlm(ng)%name,              &
     &                             dimname = 'ocean_time',              &
     &                             dimsize = inprec)
 
#if defined PIO_LIB && defined DISTRIBUTE
            CASE (io_pio)
              CALL pio_netcdf_get_dim (ng, irpm, tlm(ng)%name,          &
     &                                 dimname = 'ocean_time',          &
     &                                 dimsize = inprec)
#endif
          END SELECT
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
          tindex=1
          CALL get_state (ng, irpm, 1, tlm(ng), inprec, tindex)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
          kout=tindex
          nout=tindex
#ifdef DISTRIBUTE
          CALL wrt_dai (ng, myrank)
#else
          CALL wrt_dai (ng, -1)
#endif
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END DO
      END IF
!
!-----------------------------------------------------------------------
!  Compute and report model-observation comparison statistics.
!-----------------------------------------------------------------------
!
      IF (ldone.or.(exit_flag.eq.1)) THEN
        DO ng=1,ngrids
#ifdef DISTRIBUTE
          CALL stats_modobs (ng, myrank)
#else
          CALL stats_modobs (ng, -1)
#endif
        END DO
      END IF
!
!-----------------------------------------------------------------------
!  If blowing-up, save latest model state into RESTART NetCDF file.
!-----------------------------------------------------------------------
!
!  If cycling restart records, write solution into record 3.
!
      IF (exit_flag.eq.1) THEN
        DO ng=1,ngrids
          IF (lwrtrst(ng)) THEN
            IF (master) WRITE (stdout,20)
 20         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,30)
 30     FORMAT (/,'Elapsed wall CPU time for each process (seconds):',/)
      END IF
!
      DO ng=1,ngrids
        DO thread=thread_range
          CALL wclock_off (ng, inlm, 0, __line__, myfile)
        END DO
      END DO
!
!  Report dynamic memory and automatic memory requirements.
!
      CALL memory
!
!  Close IO files.
!
      DO ng=1,ngrids
        CALL close_inp (ng, inlm)
      END DO
      CALL close_out
!
      RETURN
      END SUBROUTINE roms_finalize
 
      END MODULE roms_kernel_mod