symmetry.h

Go to the documentation of this file.

      MODULE roms_kernel_mod
!
!git $Id$
!================================================== Hernan G. Arango ===
!  Copyright (c) 2002-2025 The ROMS Group                              !
!    Licensed under a MIT/X style license                              !
!    See License_ROMS.md                                               !
!=======================================================================
!                                                                      !
!  ROMS Representer Driver:                                            !
!                                                                      !
!  This driver executes ROMS weak constraint 4DVAR inner loop and      !
!  checks the symmetry of the H R R' H' operator,  where R' and H' are !
!  transpose of R and H, respectively. The  R' H'  term is computed by !
!  integrating the adjoint model backwards while the  H R  is computed !
!  integrating forward the tangent linear model.                       !
!                                                                      !
!  It controls the initialization,  time-stepping, and finalization of !
!  the model execution following ESMF conventions:                     !
!                                                                      !
!     ROMS_initialize                                                  !
!     ROMS_run                                                         !
!     ROMS_finalize                                                    !
!                                                                      !
!=======================================================================
!
      USE mod_param
      USE mod_parallel
      USE mod_arrays
      USE mod_fourdvar
      USE mod_iounits
      USE mod_ncparam
      USE mod_ocean
      USE mod_scalars
      USE mod_stepping
!
      USE close_io_mod,      ONLY : close_inp, close_out
      USE convolve_mod,      ONLY : error_covariance
      USE def_impulse_mod,   ONLY : def_impulse
      USE def_norm_mod,      ONLY : def_norm
#ifdef DISTRIBUTE
      USE distribute_mod,    ONLY : mp_bcastf
#endif
      USE get_state_mod,     ONLY : get_state
      USE inp_par_mod,       ONLY : inp_par
      USE normalization_mod, ONLY : normalization
      USE stdout_mod,        ONLY : set_stdoutunit, stdout_unit
      USE strings_mod,       ONLY : founderror
      USE tl_def_ini_mod,    ONLY : tl_def_ini
      USE wrt_impulse_mod,   ONLY : wrt_impulse
      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 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 :: STDrec, Tindex
      integer :: chunk_size, ng, thread
#ifdef _OPENMP
      integer :: my_threadnum
#endif
!
      character (len=*), parameter :: MyFile =                          &
     &  __FILE__//", ROMS_initialize"
 
#ifdef DISTRIBUTE
!
!-----------------------------------------------------------------------
!  Set distribute-memory (mpi) world communictor.
!-----------------------------------------------------------------------
!
      IF (PRESENT(mpicomm)) THEN
        ocn_comm_world=mpicomm
      ELSE
        ocn_comm_world=mpi_comm_world
      END IF
      CALL mpi_comm_rank (ocn_comm_world, myrank, myerror)
      CALL mpi_comm_size (ocn_comm_world, mysize, myerror)
#endif
!
!-----------------------------------------------------------------------
!  On first pass, initialize model parameters a variables for all
!  nested/composed grids.  Notice that the logical switch "first"
!  is used to allow multiple calls to this routine during ensemble
!  configurations.
!-----------------------------------------------------------------------
!
      IF (first) THEN
        first=.false.
!
!  Initialize parallel control switches. These scalars switches are
!  independent from standard input parameters.
!
        CALL initialize_parallel
!
!  Set the ROMS standard output unit to write verbose execution info.
!  Notice that the default standard out unit in Fortran is 6.
!
!  In some applications like coupling or disjointed mpi-communications,
!  it is advantageous to write standard output to a specific filename
!  instead of the default Fortran standard output unit 6. If that is
!  the case, it opens such formatted file for writing.
!
        IF (set_stdoutunit) THEN
          stdout=stdout_unit(master)
          set_stdoutunit=.false.
        END IF
!
!  Read in model tunable parameters from standard input. Allocate and
!  initialize variables in several modules after the number of nested
!  grids and dimension parameters are known.
!
        CALL inp_par (inlm)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Set domain decomposition tile partition range.  This range is
!  computed only once since the "first_tile" and "last_tile" values
!  are private for each parallel thread/node.
!
#if defined _OPENMP
      mythread=my_threadnum()
#elif defined DISTRIBUTE
      mythread=myrank
#else
      mythread=0
#endif
      DO ng=1,ngrids
        chunk_size=(ntilex(ng)*ntilee(ng)+numthreads-1)/numthreads
        first_tile(ng)=mythread*chunk_size
        last_tile(ng)=first_tile(ng)+chunk_size-1
      END DO
!
!  Initialize internal wall clocks. Notice that the timings does not
!  includes processing standard input because several parameters are
!  needed to allocate clock variables.
!
        IF (master) THEN
          WRITE (stdout,10)
 10       FORMAT (/,' Process Information:',/)
        END IF
!
        DO ng=1,ngrids
          DO thread=thread_range
            CALL wclock_on (ng, inlm, 0, __line__, myfile)
          END DO
        END DO
!
!  Allocate and initialize modules variables.
!
        CALL roms_allocate_arrays (allocate_vars)
        CALL roms_initialize_arrays
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
      END IF
!
!-----------------------------------------------------------------------
!  Read in standard deviation factors for error covariance.
!-----------------------------------------------------------------------
!
!  Initial conditions standard deviation. They are loaded in Tindex=1
!  of the e_var(...,Tindex) state variables.
!
      stdrec=1
      tindex=1
      DO ng=1,ngrids
        CALL get_state (ng, 10, 10, std(1,ng), stdrec, tindex)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END DO
!
!  Model error standard deviation. They are loaded in Tindex=2
!  of the e_var(...,Tindex) state variables.
!
      stdrec=1
      tindex=2
      IF (nsa.eq.2) THEN
        DO ng=1,ngrids
          CALL get_state (ng, 11, 11, std(2,ng), stdrec, tindex)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END DO
      END IF
 
#ifdef ADJUST_BOUNDARY
!
!  Open boundary conditions standard deviation.
!
      stdrec=1
      tindex=1
      DO ng=1,ngrids
        CALL get_state (ng, 12, 12, std(3,ng), stdrec, tindex)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END DO
#endif
#if defined ADJUST_WSTRESS || defined ADJUST_STFLUX
!
!  Surface forcing standard deviation.
!
      stdrec=1
      tindex=1
      DO ng=1,ngrids
        CALL get_state (ng, 13, 13, std(4,ng), stdrec, tindex)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END DO
#endif
!
      RETURN
      END SUBROUTINE roms_initialize
!
      SUBROUTINE roms_run (RunInterval)
!
!=======================================================================
!                                                                      !
!  This routine time-steps ROMS weak constraint 4DVAR inner loop.      !
!                                                                      !
!=======================================================================
!
!  Imported variable declarations
!
      real(dp), intent(in) :: RunInterval            ! seconds
!
!  Local variable declarations.
!
      logical :: BOUNDED_TL, Lposterior
!
      integer :: i, j, ng, tile
      integer :: Lbck, Lini, Lstate, NRMrec, Rec1, Rec2
      integer :: Fcount
      integer :: IperAD, JperAD, KperAD, ivarAD
      integer :: IoutTL, JoutTL, KoutTL, ivarTL
#ifdef DISTRIBUTE
      integer :: Istr, Iend, Jstr, Jend
#endif
      integer, dimension(Ngrids) :: Nrec
 
      integer, allocatable :: StateVar(:)
!
      real(r8), allocatable :: R(:,:,:), Rerr(:,:)
!
      character (len=8 ) :: driver
      character (len=24) :: frmt
 
      character (len=*), parameter :: MyFile =                          &
     &  __FILE__//", ROMS_run"
!
!-----------------------------------------------------------------------
!  Run nonlinear model and compute basic state tracjectory.
!-----------------------------------------------------------------------
!
      lini=1                ! NLM initial conditions record in INI
      lbck=1                ! background record in INI
      rec1=1
      rec2=2
      driver='symmetry'
!
!  Initialize and set nonlinear model initial conditions.
!
      DO ng=1,ngrids
        wrtnlmod(ng)=.false.
        wrtrpmod(ng)=.false.
        wrttlmod(ng)=.false.
      END DO
      CALL initial
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Run nonlinear model and compute basic state trajectory.
!
      DO ng=1,ngrids
        IF (master) THEN
          WRITE (stdout,10) 'NL', ng, ntstart(ng), ntend(ng)
        END IF
      END DO
!
#ifdef SOLVE3D
      CALL main3d (runinterval)
#else
      CALL main2d (runinterval)
#endif
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Set structure for the nonlinear forward trajectory to be processed
!  by the tangent linear and adjoint models. Also, set switches to
!  process the FWD structure in routine "check_multifile". Notice that
!  it is possible to split solution into multiple NetCDF files to reduce
!  their size.
!
      CALL edit_multifile ('HIS2FWD')
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      DO ng=1,ngrids
        lreadfwd(ng)=.true.
      END DO
 
#ifdef FORWARD_FLUXES
!
!  Set the BLK structure to contain the nonlinear model surface fluxes
!  needed by the tangent linear and adjoint models. Also, set switches
!  to process that structure in routine "check_multifile". Notice that
!  it is possible to split the solution into multiple NetCDF files to
!  reduce their size.
!
!  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
!
!-----------------------------------------------------------------------
!  Compute or read in background-error correlations normalization
!  factors.
!-----------------------------------------------------------------------
!
!  If computing, write out factors to NetCDF. This is an expensive
!  computation and needs to be computed once for a particular
!  application grid.
!
      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
!
!  Define TLM 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
!
!=======================================================================
!  Perturb the adjoint state variable, one at the time, with a delta
!  function at the specified perturbation point. Integrate the adjoint
!  model backwards and the force the tangent linear model with the
!  result. Then, store tangent linear solution at the requested point
!  for each state variable.
!=======================================================================
!
!  Determine variables to perturb.
!
#ifdef SOLVE3D
      allocate ( statevar(maxval(nstatevar)-2) )
      statevar(1)=isfsur
      statevar(2)=isuvel
      statevar(3)=isvvel
      DO i=1,mt
        statevar(3+i)=istvar(i)
      END DO
      lstate=maxval(statevar)
#else
      allocate ( statevar(maxval(nstatevar)) )
      statevar(1)=isfsur
      statevar(2)=isubar
      statevar(3)=isvbar
      lstate=3
#endif
      allocate ( r(lstate,lstate,ngrids) )
      allocate ( rerr(lstate,lstate) )
!
!  Initialize sample representer matrix to report.
!
      r(1:lstate,1:lstate,1:ngrids)=0.0_r8
!
!  Set point to perturb and point to report.
!
      ivartl=int(user(1))
      ivarad=int(user(2))
      iouttl=int(user(3))
      iperad=int(user(4))
      jouttl=int(user(5))
      jperad=int(user(6))
      kouttl=int(user(7))
      kperad=int(user(8))
!
!  Open tangent linear initial conditions NetCDF file (ITL struc) and
!  inquire about its variables IDs.
!
      DO ng=1,ngrids
        ldefitl(ng)=.false.
        CALL tl_def_ini (ng)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END DO
!
!  Determine number of perturbation iterations.
!
      erstr=1
#ifdef SOLVE3D
      erend=minval(nstatevar)-2
#else
      erend=minval(nstatevar)
#endif
!
      pert_loop : DO nrun=erstr,erend
        user(2)=statevar(nrun)
!
!-----------------------------------------------------------------------
!  Time-step the adjoint model.
!-----------------------------------------------------------------------
!
!  Perturb the adjoint model at specified point.
!
        admodel=.true.
        tlmodel=.false.
        DO ng=1,ngrids
          lold(ng)=1
          CALL ad_initial (ng)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
          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,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
!
!-----------------------------------------------------------------------
!  Convolve adjoint trajectory with model-error covariance and convert
!  to impulse forcing.
!-----------------------------------------------------------------------
!
        lposterior=.false.
        CALL error_covariance (itlm, driver, -1, -1,                    &
     &                         lbck, lini, lold, lnew,                  &
     &                         rec1, rec2, lposterior)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Convert convolved adjoint solution to impulse forcing. Write out
!  impulse forcing into TLF(ng)%name NetCDF file. To facilitate the
!  forcing by the TLM and RPM, the forcing is process and written in
!  increasing time coordinates.
!
        DO ng=1,ngrids
          tlf(ng)%Rindex=0
#ifdef DISTRIBUTE
          tile=myrank
#else
          tile=-1
#endif
          CALL wrt_impulse (ng, tile, iadm, adm(ng)%name)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END DO
!
!-----------------------------------------------------------------------
!  Integrate tangent linear model forced by the convolved adjoint
!  trajectory (impulse forcing) to compute R_n * PSI at observation
!  points.
!-----------------------------------------------------------------------
!
!  Initialize tangent linear model from rest. The initial contribution
!  from the adjoint model will be added as impulse forcing in
!  "tl_forcing".
!
        admodel=.false.
        tlmodel=.false.
        DO ng=1,ngrids
          wrttlmod(ng)=.false.
          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 R_n * PSI at observation points which
!  are used in the conjugate gradient algorithm.
!
        DO ng=1,ngrids
          IF (master) THEN
            WRITE (stdout,10) 'TL', ng, ntstart(ng), ntend(ng)
          END IF
        END DO
!
#ifdef SOLVE3D
        CALL tl_main3d (runinterval)
#else
        CALL tl_main2d (runinterval)
#endif
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Extract solution at requested points.
!
        DO ng=1,ngrids
#ifdef DISTRIBUTE
          istr=bounds(ng)%Istr(myrank)
          iend=bounds(ng)%Iend(myrank)
          jstr=bounds(ng)%Jstr(myrank)
          jend=bounds(ng)%Jend(myrank)
          bounded_tl=((istr.le.iouttl).and.(iouttl.le.iend)).and.       &
     &               ((jstr.le.jouttl).and.(jouttl.le.jend))
#else
          bounded_tl=.true.
#endif
          r(1,nrun,ng)=ocean(ng)%tl_zeta(iouttl,jouttl,knew(ng))
#ifndef SOLVE3D
          r(2,nrun,ng)=ocean(ng)%tl_ubar(iouttl,jouttl,knew(ng))
          r(3,nrun,ng)=ocean(ng)%tl_vbar(iouttl,jouttl,knew(ng))
#else
          r(2,nrun,ng)=ocean(ng)%tl_u(iouttl,jouttl,kouttl,nstp(ng))
          r(3,nrun,ng)=ocean(ng)%tl_v(iouttl,jouttl,kouttl,nstp(ng))
          DO i=1,nt(ng)
            r(i+3,nrun,ng)=ocean(ng)%tl_t(iouttl,jouttl,kouttl,         &
     &                                    nstp(ng),i)
          END DO
#endif
        END DO
 
      END DO pert_loop
!
!  Report sampled representer matrix and report symmetry.
!
      DO ng=1,ngrids
#ifdef DISTRIBUTE
        CALL mp_bcastf (ng, itlm, r(:,:,ng))
#endif
        IF (master) THEN
          WRITE (stdout,20) 'Representer Matrix Symmetry Test: ',       &
     &                      'Perturbing Point: ',                       &
     &                      iperad, jperad, kperad,                     &
     &                      'Sampling   Point: ',                       &
     &                      iouttl, jouttl, kouttl
          WRITE (stdout,30) 'Sampled Representer Matrix: '
          WRITE (frmt,'(a,i0,a)') '(', lstate, '(1x,1p,e14.7,0p))'
!!        WRITE (frmt,'(a)') '(*(1x,1p,e14.7,0p))'        ! FORTRAN 2008
          WRITE (stdout,'(a,a)') 'My Format: ', frmt
          DO i=1,lstate
            DO j=1,lstate
              rerr(i,j)=r(i,j,ng)-r(j,i,ng)
            END DO
          END DO
          DO i=1,lstate
            WRITE (stdout,frmt) (r(i,j,ng),j=1,lstate)
          END DO
          WRITE (stdout,30) 'Representer Matrix Symmetry Error: '
          DO i=1,lstate
            WRITE (stdout,frmt) (rerr(i,j),j=1,lstate)
          END DO
        END IF
      END DO
!
 10   FORMAT (/,1x,a,1x,'ROMS: started time-stepping:',                 &
     &        ' (Grid: ',i2.2,' TimeSteps: ',i8.8,' - ',i8.8,')',/)
 20   FORMAT (/,1x,a,/,/,3x,a,' i = ',i4.4,' j = ',i4.4,' k = ',i4.4,   &
     &                 /,3x,a,' i = ',i4.4,' j = ',i4.4,' k = ',i4.4)
 30   FORMAT (/,1x,a,/)
!
      RETURN
      END SUBROUTINE roms_run
!
      SUBROUTINE roms_finalize
!
!=======================================================================
!                                                                      !
!  This routine terminates ROMS nonlinear model execution.             !
!                                                                      !
!=======================================================================
!
      USE mod_param
      USE mod_parallel
      USE mod_iounits
      USE mod_ncparam
      USE mod_scalars
!
!  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, 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