so_semi_roms.h

Go to the documentation of this file.

      MODULE roms_kernel_mod
!
!git $Id$
!================================================== Hernan G. Arango ===
!  Copyright (c) 2002-2025 The ROMS Group            Andrew M. Moore   !
!    Licensed under a MIT/X style license                              !
!    See License_ROMS.md                                               !
!=======================================================================
!                                                                      !
!  ROMS Stochastic Optimals, Seminorm Estimation Driver:               !
!                                                                      !
!  This driver computes the eigenvectors of the stochastic optimals    !
!  operator with respect the seminorm of the chosen functional. The    !
!  stochastic  optimals  provide information about the influence of    !
!  stochastic  variations  (biases) in  ocean forcing.  They can be    !
!  used to build forecast ensembles.                                   !
!                                                                      !
!  These routines control the initialization,  time-stepping,  and     !
!  finalization of  ROMS  model following ESMF conventions:            !
!                                                                      !
!     ROMS_initialize                                                  !
!     ROMS_run                                                         !
!     ROMS_finalize                                                    !
!                                                                      !
!  Reference:                                                          !
!                                                                      !
!    Moore, A.M. et al., 2004: A comprehensive ocean prediction and    !
!      analysis system based on the tangent linear and adjoint of a    !
!      regional ocean model, Ocean Modelling, 7, 227-258.              !
!                                                                      !
!=======================================================================
!
      USE mod_param
      USE mod_parallel
      USE mod_arrays
      USE mod_forces
      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 mod_storage
!
      USE propagator_mod
!
      USE ad_def_his_mod,    ONLY : ad_def_his
      USE ad_wrt_his_mod,    ONLY : ad_wrt_his
      USE close_io_mod,      ONLY : close_file, close_inp, close_out
#ifdef CHECKPOINTING
      USE def_gst_mod,       ONLY : def_gst
      USE get_gst_mod,       ONLY : get_gst
#endif
      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 packing_mod,       ONLY : ad_unpack
      USE packing_mod,       ONLY : r_norm2
      USE stdout_mod,        ONLY : set_stdoutunit, stdout_unit
      USE strings_mod,       ONLY : founderror
#ifdef CHECKPOINTING
      USE wrt_gst_mod,       ONLY : wrt_gst
#endif
      USE wrt_rst_mod,       ONLY : wrt_rst
!
      implicit none
!
      PRIVATE :: iram_error
      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 :: 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 all model state arrays.
!
        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 for all grids first in order to compute
!  the size of the state vector, Nstate.  This size is computed in
!  routine "wpoints".
!-----------------------------------------------------------------------
 
#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 perturbation tangent linear model.
!
      DO ng=1,ngrids
        lreadfwd(ng)=.true.
        CALL ad_initial (ng)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END DO
!
!  Allocate arrays associated with Generalized Stability Theory (GST)
!  analysis.
!
      CALL allocate_storage
!
!  Initialize various IO flags.
!
      nrun=0
      DO ng=1,ngrids
        ldefadj(ng)=.true.
        lwrtadj(ng)=.true.
        ldeftlm(ng)=.true.
        lwrttlm(ng)=.true.
        lwrtper(ng)=.false.
        lcycletlm(ng)=.false.
        lcycleadj(ng)=.false.
      END DO
!
!  Initialize ARPACK parameters.
!
      lrvec=.true.                ! Compute Ritz vectors
      bmat='I'                    ! standard eigenvalue problem
      which='LM'                  ! compute NEV largest eigenvalues
      howmany='A'                 ! compute NEV Ritz vectors
      DO ng=1,ngrids
        ido(ng)=0                 ! reverse communication flag
        info(ng)=0                ! random initial residual vector
        iparam(1,ng)=1            ! exact shifts
        iparam(3,ng)=maxitergst   ! maximum number of Arnoldi iterations
        iparam(4,ng)=1            ! block size in the recurrence
        iparam(7,ng)=1            ! type of eigenproblem being solved
      END DO
!
!  ARPACK debugging parameters.
!
      logfil=stdout               ! output logical unit
      ndigit=-3                   ! number of decimal digits
      msaupd=1                    ! iterations, timings, Ritz
      msaup2=1                    ! norms, Ritz values
      msaitr=0
      mseigt=0
      msapps=0
      msgets=0
      mseupd=0
!
!  Determine size of the eigenproblem (Nsize) and size of work space
!  array SworkL (LworkL).
!
      DO ng=1,ngrids
        nconv(ng)=0
        nsize(ng)=nend(ng)-nstr(ng)+1
      END DO
 
#ifdef CHECKPOINTING
!
!  If restart, read in checkpointing data GST restart NetCDF file.
!  Otherwise, create checkpointing restart NetCDF file.
!
      DO ng=1,ngrids
        IF (lrstgst) THEN
          CALL get_gst (ng, iadm)
          ido(ng)=-2
          laup2(1)=.false.        ! cnorm
          laup2(2)=.false.        ! getv0
          laup2(3)=.false.        ! initv
          laup2(4)=.false.        ! update
          laup2(5)=.true.         ! ushift
        ELSE
          CALL def_gst (ng, iadm)
        END IF
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END DO
#endif
!
      RETURN
      END SUBROUTINE roms_initialize
!
      SUBROUTINE roms_run (RunInterval)
!
!=======================================================================
!                                                                      !
!  This routine computes the eigenvectors of the stochastic optimal    !
!  matrix, defined as:                                                 !
!                                                                      !
!       S = INTEGRAL [ R'(t,To) X R(t,To) dt ]   from To to T          !
!                                                                      !
!  where T is the  forecast interval from the initial state at To,     !
!  R(t,To) is  the  forward  tangent  propagator of the linearized     !
!  dynamical model,  R'(t,To) is the adjoint of  R(t,To),  and the     !
!  matrix X defines the norm of interest. Here, X is a seminorm of     !
!  the chosen functional.                                              !
!                                                                      !
!  The above operator is symmetric and the  ARPACK library is used     !
!  to select eigenvectors and eigenvalues:                             !
!                                                                      !
!  Lehoucq, R.B., D.C. Sorensen, and C. Yang, 1997:  ARPACK user's     !
!    guide:  solution  of  large  scale  eigenvalue  problems with     !
!    implicit restarted Arnoldi Methods, Rice University, 140p.        !
!                                                                      !
!=======================================================================
!
!  Imported variable declarations
!
      real(dp), intent(in) :: RunInterval            ! seconds
!
!  Local variable declarations.
!
      logical :: ITERATE
#ifdef CHECKPOINTING
      logical :: LwrtGST
#endif
!
      integer :: Fcount, Is, Ie
      integer :: i, iter, ng, tile
      integer :: NconvRitz(Ngrids)
!
      real(r8) :: Enorm
!
      TYPE (T_GST), allocatable :: ad_state(:)
      TYPE (T_GST), allocatable :: state(:)
!
      character (len=55) :: string
 
      character (len=*), parameter :: MyFile =                          &
     &  __FILE__//", ROMS_run"
!
!-----------------------------------------------------------------------
!  Implicit Restarted Arnoldi Method (IRAM) for the computation of
!  optimal perturbation Ritz eigenfunctions.
!-----------------------------------------------------------------------
!
!  Allocate nested grid pointers for state vectors.
!
      IF (.not.allocated(ad_state)) THEN
        allocate ( ad_state(ngrids) )
      END IF
      IF (.not.allocated(state)) THEN
        allocate ( state(ngrids) )
      END IF
!
!  Iterate until either convergence or maximum iterations has been
!  exceeded.
!
      iter=0
      iterate=.true.
#ifdef CHECKPOINTING
      lwrtgst=.true.
#endif
!
      iter_loop : DO WHILE (iterate)
        iter=iter+1
!
!  Reverse communication interface.
!
        DO ng=1,ngrids
#ifdef PROFILE
          CALL wclock_on (ng, iadm, 38, __line__, myfile)
#endif
#ifdef DISTRIBUTE
          CALL pdsaupd (ocn_comm_world,                                 &
     &                  ido(ng), bmat, nsize(ng), which, nev,           &
     &                  ritz_tol,                                       &
     &                  storage(ng)%resid(nstr(ng)), ncv,               &
     &                  storage(ng)%Bvec(nstr(ng),1), nsize(ng),        &
     &                  iparam(1,ng), ipntr(1,ng),                      &
     &                  storage(ng)%SworkD,                             &
     &                  sworkl(1,ng), lworkl, info(ng))
#else
          CALL dsaupd (ido(ng), bmat, nsize(ng), which, nev,            &
     &                 ritz_tol,                                        &
     &                 storage(ng)%resid, ncv,                          &
     &                 storage(ng)%Bvec, nsize(ng),                     &
     &                 iparam(1,ng), ipntr(1,ng),                       &
     &                 storage(ng)%SworkD,                              &
     &                 sworkl(1,ng), lworkl, info(ng))
#endif
          nconv(ng)=iaup2(4)
#ifdef PROFILE
          CALL wclock_off (ng, iadm, 38, __line__, myfile)
#endif
#ifdef CHECKPOINTING
!
!  If appropriate, write out check point data into GST restart NetCDF
!  file. Notice that the restart data is always saved if MaxIterGST
!  is reached without convergence. It is also saved when convergence
!  is achieved (ido=99).
!
          IF ((mod(iter,ngst).eq.0).or.(iter.ge.maxitergst).or.         &
     &        (any(ido.eq.99))) THEN
            CALL wrt_gst (ng, iadm)
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
          END IF
#endif
        END DO
!
!  Terminate computations if maximum number of iterations is reached.
!  This will faciliate splitting the analysis in several computational
!  cycles using the restart option.
!
        IF ((iter.ge.maxitergst).and.any(ido.ne.99)) THEN
          iterate=.false.
          EXIT iter_loop
        END IF
!
!  Perform matrix-vector operation:  R'(t,0)XR(0,t)u
!
        IF (any(abs(ido).eq.1)) THEN
          DO ng=1,ngrids
            fcount=adm(ng)%load
            adm(ng)%Nrec(fcount)=0
            fcount=tlm(ng)%load
            tlm(ng)%Nrec(fcount)=0
            adm(ng)%Rindex=0
            tlm(ng)%Rindex=0
          END DO
!
!  Set state vectors to process by the propagator via pointer
!  equivalence.
!
          DO ng=1,ngrids
            IF (ASSOCIATED(state(ng)%vector)) THEN
              nullify (state(ng)%vector)
            END IF
            is=ipntr(1,ng)
            ie=is+nsize(ng)-1
            state(ng)%vector => storage(ng)%SworkD(is:ie)
 
            IF (ASSOCIATED(ad_state(ng)%vector)) THEN
              nullify (ad_state(ng)%vector)
            END IF
            is=ipntr(2,ng)
            ie=is+nsize(ng)-1
            ad_state(ng)%vector => storage(ng)%SworkD(is:ie)
          END DO
 
          CALL propagator_so_semi (runinterval, state, ad_state)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        ELSE
          IF (any(info.ne.0)) THEN
            DO ng=1,ngrids
              IF (info(ng).ne.0) THEN
                IF (master) THEN
                  CALL iram_error (info(ng), string)
                  WRITE (stdout,10) 'DSAUPD', trim(string),             &
     &                              ', info = ', info(ng)
                END IF
                RETURN
              END IF
            END DO
          ELSE
!
!  Compute Ritz vectors (the only choice left is IDO=99).  They are
!  generated in ARPACK in decreasing magnitude of its eigenvalue.
!  The most significant is first.
!
            DO ng=1,ngrids
              nconvritz(ng)=iparam(5,ng)
              IF (master) THEN
                WRITE (stdout,20) 'Number of converged Ritz values:',   &
     &                            iparam(5,ng)
                WRITE (stdout,20) 'Number of Arnoldi iterations:',      &
     &                            iparam(3,ng)
              END IF
#ifdef PROFILE
              CALL wclock_on (ng, iadm, 38, __line__, myfile)
#endif
#ifdef DISTRIBUTE
              CALL pdseupd (ocn_comm_world,                             &
     &                      lrvec, howmany, pick(1,ng),                 &
     &                      rvaluer(1,ng),                              &
     &                      storage(ng)%Rvector(nstr(ng),1),            &
     &                      nsize(ng), sigmar,                          &
     &                      bmat, nsize(ng), which, nev, ritz_tol,      &
     &                      storage(ng)%resid(nstr(ng)), ncv,           &
     &                      storage(ng)%Bvec(nstr(ng),1), nsize(ng),    &
     &                      iparam(1,ng), ipntr(1,ng),                  &
     &                      storage(ng)%SworkD,                         &
     &                      sworkl(1,ng), lworkl, info(ng))
#else
              CALL dseupd (lrvec, howmany, pick(1,ng),                  &
     &                     rvaluer(1,ng),                               &
     &                     storage(ng)%Rvector, nsize(ng),              &
     &                     sigmar, bmat, nsize, which, nev, ritz_tol,   &
     &                     storage(ng)%resid, ncv,                      &
     &                     storage(ng)%Bvec, nsize(ng),                 &
     &                     iparam(1,ng), ipntr(1,ng),                   &
     &                     storage(ng)%SworkD,                          &
     &                     sworkl(1,ng), lworkl, info(ng))
#endif
#ifdef PROFILE
              CALL wclock_off (ng, iadm, 38, __line__, myfile)
#endif
            END DO
 
            IF (any(info.ne.0)) THEN
              DO ng=1,ngrids
                IF (info(ng).ne.0) THEN
                  IF (master) THEN
                    CALL iram_error (info(ng), string)
                    WRITE (stdout,10) 'DSEUPD', trim(string),           &
     &                                ', info = ', info(ng)
                  END IF
                  RETURN
                END IF
              END DO
            ELSE
!
!  Close existing adjoint NetCDF file.
!
             sourcefile=myfile
             DO ng=1,ngrids
               CALL close_file (ng, iadm, adm(ng), adm(ng)%name)
               IF (founderror(exit_flag, noerror,                       &
     &                        __line__, myfile)) RETURN
             END DO
!
!  Clear forcing arrays.
!
             DO ng=1,ngrids
               DO tile=first_tile(ng),last_tile(ng),+1
                 CALL initialize_forces (ng, tile, 0)
               END DO
             END DO
!
!  Activate writing of each eigenvector into tangent history NetCDF
!  file. The "ocean_time" is the eigenvector number.
!
              nrun=0
 
              DO i=1,maxval(nconvritz)
                DO ng=1,ngrids
                  IF ((i.eq.1).or.lmultigst) THEN
                    fcount=adm(ng)%load
                    adm(ng)%Nrec(fcount)=0
                    adm(ng)%Rindex=0
                    lwrtadj(ng)=.true.
                    lwrtstate2d(ng)=.true.
                  END IF
                  IF (lmultigst) THEN
                    ldefadj(ng)=.true.
                    WRITE (adm(ng)%name,30) trim(adm(ng)%base), i
                    CALL ad_def_his (ng, ldefadj(ng))
                    IF (founderror(exit_flag, noerror,                  &
     &                             __line__, myfile)) RETURN
                  ELSE IF (.not.lmultigst.and.(i.eq.1)) THEN
                    ldefadj(ng)=.true.
                    adm(ng)%name=trim(adm(ng)%base)//'_ritz.nc'
                    CALL ad_def_his (ng, ldefadj(ng))
                    IF (founderror(exit_flag, noerror,                  &
     &                             __line__, myfile)) RETURN
                  END IF
                END DO
!
!  Compute and write the Ritz eigenvectors.
!
                DO ng=1,ngrids
                  is=nstr(ng)
                  ie=nend(ng)
                  IF (ASSOCIATED(state(ng)%vector)) THEN
                    nullify (state(ng)%vector)
                  END IF
 
                  IF (ASSOCIATED(ad_state(ng)%vector)) THEN
                    nullify (ad_state(ng)%vector)
                  END IF
                  state(ng)%vector => storage(ng)%Rvector(is:ie,i)
                  ad_state(ng)%vector => sworkr(is:ie)
                END DO
 
                CALL propagator_so_semi (runinterval, state, ad_state)
                IF (founderror(exit_flag, noerror,                      &
     &                         __line__, myfile)) RETURN
!
!  Unpack surface forcing eigenvectors from Rvector and write into
!  nonlinear history file.
!
                DO ng=1,ngrids
                  DO tile=first_tile(ng),last_tile(ng),+1
                    CALL ad_unpack (ng, tile, nstr(ng), nend(ng),       &
     &                              state(ng)%vector)
                  END DO
#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
!
!  Compute Euclidean norm.
!
                DO ng=1,ngrids
                  CALL r_norm2 (ng, iadm, nstr(ng), nend(ng),           &
     &                          -rvaluer(i,ng),                         &
     &                          state(ng)%vector,                       &
     &                          ad_state(ng)%vector, enorm)
                  norm(i,ng)=enorm
                  IF (master) THEN
                    WRITE (stdout,40) i, norm(i,ng), rvaluer(i,ng), i
                  END IF
                END DO
!
!  Write out Ritz eigenvalues and Ritz eigenvector Euclidean norm
!  (residual) to nonlinear history NetCDF file.
!
                sourcefile=myfile
                DO ng=1,ngrids
                  SELECT CASE (adm(ng)%IOtype)
                    CASE (io_nf90)
                      CALL netcdf_put_fvar (ng, iadm,                   &
     &                                      adm(ng)%name,               &
     &                                      'Ritz_rvalue',              &
     &                                      rvaluer(i:,ng),             &
     &                                      start = (/adm(ng)%Rindex/), &
     &                                      total = (/1/),              &
     &                                      ncid = adm(ng)%ncid)
 
                      IF (founderror(exit_flag, noerror,                &
     &                               __line__, myfile)) RETURN
!
                      CALL netcdf_put_fvar (ng, iadm, adm(ng)%name,     &
     &                                      'Ritz_norm',                &
     &                                      norm(i:,ng),                &
     &                                      start = (/adm(ng)%Rindex/), &
     &                                      total = (/1/),              &
     &                                      ncid = adm(ng)%ncid)
 
                      IF (founderror(exit_flag, noerror,                &
     &                               __line__, myfile)) RETURN
 
                      IF ((i.eq.1).or.lmultigst) THEN
                        CALL netcdf_put_fvar (ng, iadm,                 &
     &                                        adm(ng)%name,             &
     &                                        'SO_trace',               &
     &                                        trnorm(ng),               &
     &                                        start = (/0/),            &
     &                                        total = (/0/),            &
     &                                        ncid = adm(ng)%ncid)
 
                        IF (founderror(exit_flag, noerror,              &
     &                                 __line__, myfile)) RETURN
                      END IF
 
#if defined PIO_LIB && defined DISTRIBUTE
                    CASE (io_pio)
                      CALL pio_netcdf_put_fvar (ng, iadm,               &
     &                                          adm(ng)%name,           &
     &                                          'Ritz_rvalue',          &
     &                                          rvaluer(i:,ng),         &
     &                                      start = (/adm(ng)%Rindex/), &
     &                                      total = (/1/),              &
     &                                      piofile = adm(ng)%pioFile)
 
                      IF (founderror(exit_flag, noerror,                &
     &                               __line__, myfile)) RETURN
!
                      CALL pio_netcdf_put_fvar (ng, iadm,               &
     &                                          adm(ng)%name,           &
     &                                          'Ritz_norm',            &
     &                                          norm(i:,ng),            &
     &                                      start = (/adm(ng)%Rindex/), &
     &                                      total = (/1/),              &
     &                                      piofile = adm(ng)%pioFile)
 
                      IF (founderror(exit_flag, noerror,                &
     &                               __line__, myfile)) RETURN
 
                      IF ((i.eq.1).or.lmultigst) THEN
                        CALL pio_netcdf_put_fvar (ng, iadm,             &
     &                                            adm(ng)%name,         &
     &                                            'SO_trace',           &
     &                                            trnorm(ng),           &
     &                                        start = (/0/),            &
     &                                        total = (/0/),            &
     &                                        piofile = adm(ng)%pioFile)
 
                        IF (founderror(exit_flag, noerror,              &
     &                                 __line__, myfile)) RETURN
                      END IF
#endif
                  END SELECT
!
                  IF (lmultigst) THEN
                    CALL close_file (ng, iadm, adm(ng), adm(ng)%name)
                    IF (founderror(exit_flag, noerror,                  &
     &                             __line__, myfile)) RETURN
                  END IF
                END DO
              END DO
            END IF
          END IF
          iterate=.false.
        END IF
 
      END DO iter_loop
!
 10   FORMAT (/,1x,'Error in ',a,1x,a,a,1x,i5,/)
 20   FORMAT (/,a,1x,i2,/)
 30   FORMAT (a,'_',i3.3,'.nc')
 40   FORMAT (1x,i4.4,'-th residual',1p,e14.6,0p,                       &
     &        '  Ritz value',1pe14.6,0p,2x,i4.4)
!
      RETURN
      END SUBROUTINE roms_run
!
      SUBROUTINE roms_finalize
!
!=======================================================================
!                                                                      !
!  This routine terminates ROMS nonlinear and adjoint models           !
!  execution.                                                          !
!                                                                      !
!=======================================================================
!
!  Local variable declarations.
!
      integer :: Fcount, ng, thread
!
      character (len=*), parameter :: MyFile =                          &
     &  __FILE__//", ROMS_finalize"
!
!-----------------------------------------------------------------------
!  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
      END SUBROUTINE roms_finalize
!
      SUBROUTINE iram_error (info, string)
!
!=======================================================================
!                                                                      !
!  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
      END SUBROUTINE iram_error
 
      END MODULE roms_kernel_mod