hessian_op_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 Optimal Pertubations (Singular Vectors) Driver:                !
!                                                                      !
!  This driver computes the singular vectors of the propagator R(0,t)  !
!  which measure the  fastest  growing of all possible  perturbations  !
!  over a given time interval.  They  are  usually  used to study the  !
!  dynamics,  sensitivity,  and stability of the ocean circulation to  !
!  naturally occurring pertubations in the prediction system.          !
!                                                                      !
!  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_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 mod_storage
      USE mod_netcdf
!
      USE propagator_mod
!
      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
#ifdef DISTRIBUTE
      USE distribute_mod,    ONLY : mp_bcastf, mp_bcasti
#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 : 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 (itlm)
        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, itlm, 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 tangent linear 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 tl_initial (ng)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END DO
!
!-----------------------------------------------------------------------
!  Read in Lanczos algorithm coefficients (cg_beta, cg_delta) from
!  file LCZ(ng)%name NetCDF (I4D-Var adjoint file), as computed in the
!  I4D-Var Lanczos data assimilation algorithm for the first outer
!  loop.  They are needed here, in routine "tl_inner2state", to compute
!  the tangent linear model initial conditions as the weighted sum
!  of the Lanczos vectors. The weighting coefficient are computed
!  by solving a tri-diagonal system that uses "cg_beta" and "cg_gamma".
!-----------------------------------------------------------------------
!
      sourcefile=myfile
      DO ng=1,ngrids
        SELECT CASE (lcz(ng)%IOtype)
          CASE (io_nf90)
            CALL netcdf_get_fvar (ng, iadm, lcz(ng)%name,               &
     &                            'cg_beta', cg_beta)
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
            CALL netcdf_get_fvar (ng, iadm, lcz(ng)%name,               &
     &                            'cg_delta', cg_delta)
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
#if defined PIO_LIB && defined DISTRIBUTE
          CASE (io_pio)
            CALL pio_netcdf_get_fvar (ng, iadm, lcz(ng)%name,           &
     &                                'cg_beta', cg_beta)
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
            CALL pio_netcdf_get_fvar (ng, iadm, lcz(ng)%name,           &
     &                                'cg_delta', cg_delta)
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
#endif
        END SELECT
      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.
        nadj(ng)=ntimes(ng)
        ntlm(ng)=ntimes(ng)
      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)=ninner
      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, itlm)
          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, itlm)
        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 singular vectors of R(0,t) by a single    !
!  integration  of  a  perturbation "u"  forward  in time  with the    !
!  tangent linear model over [0,t],  multiplication  of the  result    !
!  by "X",  followed by an integration of the  result  backwards in    !
!  time with the  adjoint model over [t,0].  This is  equivalmet to    !
!  the matrix-vector operation:                                        !
!                                                                      !
!       transpose[R(t,0)] X R(0,t) u                                   !
!                                                                      !
!  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, i, iter, ng
      integer :: NconvRitz(Ngrids)
!
      real(r8) :: Enorm
 
      real(r8), dimension(2) :: my_norm, my_value
!
      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, itlm, 38, __line__, myfile)
#endif
          IF (master) THEN
            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))
            nconv(ng)=iaup2(4)
          END IF
#ifdef PROFILE
          CALL wclock_off (ng, itlm, 38, __line__, myfile)
#endif
#ifdef DISTRIBUTE
!
!  Broadcast various Arnoldi iteration variables to all member in the
!  group.  The Arnoldi problem solved here is very small so the ARPACK
!  library is run by the master node and all relevant variables are
!  broadcasted to all nodes.
!
          CALL mp_bcasti (ng, itlm, nconv(ng))
          CALL mp_bcasti (ng, itlm, ido(ng))
          CALL mp_bcasti (ng, itlm, info(ng))
          CALL mp_bcasti (ng, itlm, iparam(:,ng))
          CALL mp_bcasti (ng, itlm, ipntr(:,ng))
          CALL mp_bcastf (ng, itlm, storage(ng)%SworkD)
#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, itlm)
            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_hop (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, itlm, 38, __line__, myfile)
#endif
              IF (master) THEN
                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))
              END IF
#ifdef DISTRIBUTE
!
!  Broadcast various Arnoldi iteration variables to all member in the
!  group.
!
              CALL mp_bcasti (ng, itlm, info(ng))
              CALL mp_bcasti (ng, itlm, iparam(:,ng))
              CALL mp_bcasti (ng, itlm, ipntr(:,ng))
              CALL mp_bcastf (ng, itlm, rvaluer(:,ng))
              CALL mp_bcastf (ng, itlm, storage(ng)%Rvector)
#endif
#ifdef PROFILE
              CALL wclock_off (ng, itlm, 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
!
!  Activate writing of each eigenvector into the  adjoint and tangent
!  linear history NetCDF files. 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
                    fcount=tlm(ng)%load
                    tlm(ng)%Nrec(fcount)=0
                    adm(ng)%Rindex=0
                    tlm(ng)%Rindex=0
                  END IF
                  IF (lmultigst) THEN
                    ldefadj(ng)=.true.
                    ldeftlm(ng)=.true.
                    WRITE (adm(ng)%name,30) trim(adm(ng)%base), i
                    WRITE (tlm(ng)%name,30) trim(tlm(ng)%base), i
                  END IF
                END DO
!
!  Compute and write Ritz eigenvectors.
!
                DO ng=1,ngrids
                  is=1
                  ie=ninner
                  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_hop (runinterval, state, ad_state)
                IF (founderror(exit_flag, noerror,                      &
     &                         __line__, myfile)) RETURN
!
                DO ng=1,ngrids
                  CALL r_norm2 (ng, itlm, 1, ninner,                    &
     &                          -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 NetCDF file(s).  Notice that we write the same value
!  twice in the TLM file for the initial and final perturbation of
!  the eigenvector.
!
                sourcefile=myfile
                DO ng=1,ngrids
                  my_norm(1)=norm(i,ng)
                  my_norm(2)=my_norm(1)
                  my_value(1)=rvaluer(i,ng)
                  my_value(2)=my_value(1)
!
                  IF (lwrttlm(ng)) THEN
                    SELECT CASE (tlm(ng)%IOtype)
                      CASE (io_nf90)
                        CALL netcdf_put_fvar (ng, itlm,                 &
     &                                        tlm(ng)%name,             &
     &                                        'Ritz_rvalue',            &
     &                                        my_value,                 &
     &                                    start = (/tlm(ng)%Rindex-1/), &
     &                                    total = (/2/),                &
     &                                    ncid = tlm(ng)%ncid)
 
                        IF (founderror(exit_flag, noerror,              &
     &                                 __line__, myfile)) RETURN
!
                        CALL netcdf_put_fvar (ng, itlm,                 &
     &                                        tlm(ng)%name,             &
     &                                        'Ritz_norm',              &
     &                                        my_norm,                  &
     &                                    start = (/tlm(ng)%Rindex-1/), &
     &                                    total = (/2/),                &
     &                                    ncid = tlm(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,             &
     &                                            tlm(ng)%name,         &
     &                                            'Ritz_rvalue',        &
     &                                            my_value,             &
     &                                    start = (/tlm(ng)%Rindex-1/), &
     &                                    total = (/2/),                &
     &                                    piofile = tlm(ng)%pioFile)
 
                        IF (founderror(exit_flag, noerror,              &
     &                                 __line__, myfile)) RETURN
!
                        CALL pio_netcdf_put_fvar (ng, itlm,             &
     &                                            tlm(ng)%name,         &
     &                                            'Ritz_norm',          &
     &                                            my_norm,              &
     &                                    start = (/tlm(ng)%Rindex-1/), &
     &                                    total = (/2/),                &
     &                                    piofile = tlm(ng)%pioFile)
 
                        IF (founderror(exit_flag, noerror,              &
     &                             __line__, myfile)) RETURN
#endif
                    END SELECT
!
                    IF (lmultigst) THEN
                      CALL close_file (ng, itlm, tlm(ng), tlm(ng)%name)
                      IF (founderror(exit_flag, noerror,                &
     &                               __line__, myfile)) RETURN
                    END IF
                  END IF
!
                  IF (lwrtadj(ng)) THEN
                    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 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
#endif
                    END SELECT
!
                    IF (lmultigst) THEN
                      CALL netcdf_close (ng, iadm, adm(ng)%ncid,        &
     &                                   adm(ng)%name)
                      IF (founderror(exit_flag, noerror,                &
     &                               __line__, myfile)) RETURN
                    END IF
                  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, itlm, 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, itlm)
      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