mod_storage.F

Go to the documentation of this file.

#include "cppdefs.h"
      MODULE mod_storage
 
#if defined PROPAGATOR
!
!git $Id$
!================================================== Hernan G. Arango ===
!  Copyright (c) 2002-2025 The ROMS Group                              !
!    Licensed under a MIT/X style license                              !
!    See License_ROMS.md                                               !
!=======================================================================
!                                                                      !
!  ROMS Generalized Stability Theory (GST) Analysis: ARPACK            !
!                                                                      !
!  LworkL       Size of Arnoldi iterations work array SworkL.          !
!  Lrvec        ARPACK logical to compute converged Ritz values.       !
!  NCV          Number of Lanczos vectors to compute.                  !
!  NEV          Number of eigenvalues to compute.                      !
!  Bvec         Lanczos/Arnoldi basis vectors.                         !
!  RvalueR      Real Ritz eigenvalues.                                 !
!  RvalueI      Imaginary Ritz eigenvalues.                            !
!  Rvector      Real Ritz eigenvectors.                                !
!  Swork        FULL state work array used in distributed memory       !
!                 communications.                                      !
!  SworkR       FULL state work array used to compute the ROMS state   !
!                 eigenvectors.                                        !
!  SworkD       State work array for ARPACK reverse communications.    !
!  SworkEV      ARPACK work array.                                     !
!  SworkL       ARPACK work array.                                     !
!  bmat         ARPACK eigenvalue value problem identifier:            !
!                 bmat='I'     standard eigenvalue problem             !
!                 bmat='G'     generalized eigenvalue problem          !
!  howmany      ARPACK form of basis functions identifier:             !
!                 howmany='A'  compute NEV Ritz vectors                !
!                 howmany='P'  compute NEV Schur vectors               !
!                 howmany='S'  compute some Ritz vectors using select  !
!  ido          ARPACK reverse communications flag (input/output).     !
!  iparam       ARPACK input/output integer parameters.                !
!  ipntr        ARPACK pointer to mark the starting location in SworkD !
!                 and SworkK arrays used in the Arnoldi iterations.    !
!  info         ARPACK Information (input) and error flag (output).    !
!  norm         Euclidean norm.                                        !
!  pick         ARPACK logical switch of Ritz vectors to compute.      !
!  resid        Initial/final residual vector.                         !
!  sigmaR       ARPACK real part of the shifts (not referenced).       !
!  sigmaI       ARPACK imaginary part of the shifts (not referenced).  !
# ifdef SO_SEMI
!  so_state     Stochastic optimals adjoint state surface forcing      !
!                 vector sample in time.                               !
# endif
!  which        ARPACK Ritz eigenvalues to compute identifier:         !
!                 which='LA'   compute NEV largest (algebraic)         !
!                 which='SA'   compute NEV smallest (algebraic)        !
!                 which='LM'   compute NEV largest in magnitude        !
!                 which='SM'   compute NEV smallest in magnitude       !
!                 which='BE'   compute NEV from each end of spectrum   !
!                                                                      !
!=======================================================================
!
      USE mod_param
!
      implicit none
!
      PUBLIC :: allocate_storage
      PUBLIC :: deallocate_storage
!
!-----------------------------------------------------------------------
!  Define T_STORAGE structure and other module variables.
!-----------------------------------------------------------------------
!
!  State work arrays structure for nested grids.
!
      TYPE t_storage
 
        real(r8), pointer :: bvec(:,:)              ! [Nstr:Nend,NCV]
        real(r8), pointer :: rvector(:,:)           ! [Nstr:Nend,NEV+1]
        real(r8), pointer :: sworkd(:)              ! [3*Nstate]
        real(r8), pointer :: resid(:)               ! [Nstr:Nend]
# ifdef STOCHASTIC_OPT
        real(r8), pointer :: ad_work(:)             ! [Mstate]
# endif
# if defined STOCHASTIC_OPT && defined HESSIAN_SO
        real(r8), pointer :: my_state(:)            ! [Nstr:Nend]
# endif
# ifdef SO_SEMI
        real(r8), pointer :: so_state(:,:)          ! [Nstr:Nend,Nsemi]
# endif
 
      END TYPE t_storage
!
      TYPE (t_storage), allocatable :: storage(:)
!
!  ARPACK logical switches.
!
      logical :: lrvec                    ! Ritz values switch
      logical, allocatable :: pick(:,:)   ! Ritz vectors, [NCV,Ngrids]
!
!  Eigenproblem parameters.
!
      integer :: ncv                      ! number of Lanczos vectors
      integer :: nev                      ! number of eigenvalues
      integer :: lworkl                   ! size of array SworkL
 
      integer, allocatable :: lworkd(:)   ! size of array SworkD
!
!  ARPACK integer parameters and flags.
!
      integer, allocatable :: iparam(:,:) ! eigenproblem parameters
      integer, allocatable :: ipntr(:,:)  ! index for work arrays
      integer, allocatable :: ido(:)      ! reverse communication flag
      integer, allocatable :: info(:)     ! information and error flag
!
!  ARPACK floating-point parameters.
!
      real(r8) :: sigmai                  ! real part shifts
      real(r8) :: sigmar                  ! imaginary part shifts
!
!  ARPACK character parameters.
!
      character (len=1) :: bmat           ! eigenvalue problem ID
      character (len=1) :: howmany        ! basis functions ID
      character (len=2) :: which          ! Ritz values ID
!
!  ARPACK work arrays.
!
      real(r8), allocatable :: rvaluer(:,:)       ! [NEV+1,Ngrids]
      real(r8), allocatable :: rvaluei(:,:)       ! [NEV+1,Ngrids]
      real(r8), allocatable :: norm(:,:)          ! [NEV+1,Ngrids]
      real(r8), allocatable :: sworkl(:,:)        ! [LworkL,Ngrids]
      real(r8), allocatable :: sworkev(:,:)       ! [3*NCV,Ngrids]
# ifdef DISTRIBUTE
      real(r8), allocatable :: swork(:)           ! [Mstate]
# endif
      real(r8), pointer :: sworkr(:)              ! [Mstate]
!
!  ARPACK private common blocks containing parameters needed for
!  checkpoiniting. The original include files "i_aupd.h" and
!  "idaup2.h" have several parameters in their commom blocks. All
!  these values are compacted here in vector arrays to allow IO
!  manipulations during checkpointing.
!
      integer  :: iaitr(8), iaup2(8), iaupd(20)
      logical  :: laitr(5), laup2(5)
      real(r8) :: raitr(8), raup2(2)
!
      common /i_aupd/ iaupd
# ifdef DOUBLE_PRECISION
      common /idaitr/ iaitr
      common /ldaitr/ laitr
      common /rdaitr/ raitr
      common /idaup2/ iaup2
      common /ldaup2/ laup2
      common /rdaup2/ raup2
# else
      common /isaitr/ iaitr
      common /lsaitr/ laitr
      common /rsaitr/ raitr
      common /isaup2/ iaup2
      common /lsaup2/ laup2
      common /rsaup2/ raup2
# endif
!
!  ARPACK debugging common block.
!
      integer :: logfil, ndigit, mgetv0
      integer :: msaupd, msaup2, msaitr, mseigt, msapps, msgets, mseupd
      integer :: mnaupd, mnaup2, mnaitr, mneigh, mnapps, mngets, mneupd
      integer :: mcaupd, mcaup2, mcaitr, mceigh, mcapps, mcgets, mceupd
!
      common /debug/                                                    &
     &       logfil, ndigit, mgetv0,                                    &
     &       msaupd, msaup2, msaitr, mseigt, msapps, msgets, mseupd,    &
     &       mnaupd, mnaup2, mnaitr, mneigh, mnapps, mngets, mneupd,    &
     &       mcaupd, mcaup2, mcaitr, mceigh, mcapps, mcgets, mceupd
!
      CONTAINS
!
      SUBROUTINE allocate_storage
!
!=======================================================================
!                                                                      !
!  This routine allocates and initialize module variables. For now,    !
!  only non-nested applications are considered.                        !
!                                                                      !
!=======================================================================
!
      USE mod_scalars
!
!  Local variable declarations
!
      integer :: i, j, ng
 
      integer, allocatable :: mstr(:), mend(:)
!
      real(r8), parameter :: inival = 0.0_r8
!
!-----------------------------------------------------------------------
!  Allocate module variables.
!-----------------------------------------------------------------------
!
!  Allocate local dimension parameters.
!
      allocate ( mstr(ngrids) )
      allocate ( mend(ngrids) )
!
!  Allocate parameters associated with nested grids.
!
      allocate ( ido(ngrids) )
      allocate ( info(ngrids) )
      allocate ( lworkd(ngrids) )
 
      allocate ( iparam(11,ngrids) )
# if defined AFT_EIGENMODES || defined FT_EIGENMODES
      allocate ( ipntr(14,ngrids) )
# else
      allocate ( ipntr(11,ngrids) )
# endif
!
!  Determine size of work array SworkL:
!
# if defined OPT_PERTURBATION
      lworkl=ncv*(ncv+8)
# elif defined HESSIAN_SV
      lworkl=ncv*(ncv+8)
# elif defined HESSIAN_FSV
      lworkl=ncv*(ncv+8)
# elif defined FORCING_SV
      lworkl=ncv*(ncv+8)
# elif defined STOCHASTIC_OPT
      lworkl=ncv*(ncv+8)
# elif defined SO_SEMI
      lworkl=ncv*(ncv+8)
# elif defined FT_EIGENMODES || defined AFT_EIGENMODES
      lworkl=3*ncv*ncv+6*ncv
# endif
# ifdef SO_SEMI
      DO ng=1,ngrids
        nsemi(ng)=1+ntimes(ng)/nadj(ng)
      END DO
# endif
!
!  Set local dimensions of storage arrays according to the propagator
!  driver.
!
      DO ng=1,ngrids
# if defined HESSIAN_FSV || defined HESSIAN_SO || defined HESSIAN_SV
        mstr(ng)=1
        mend(ng)=ninner
        lworkd(ng)=3*ninner
# else
        mstr(ng)=nstr(ng)
        mend(ng)=nend(ng)
        lworkd(ng)=3*nstate(ng)
# endif
      END DO
!
!  Allocate structure.
!
      allocate ( storage(ngrids) )
 
      DO ng=1,ngrids
        allocate ( storage(ng) % Bvec(mstr(ng):mend(ng),ncv) )
        dmem(ng)=dmem(ng)+real((mend(ng)-mstr(ng))*ncv,r8)
 
        allocate ( storage(ng) % Rvector(mstr(ng):mend(ng),nev+1) )
        dmem(ng)=dmem(ng)+real((mend(ng)-mstr(ng))*(nev+1),r8)
 
        allocate ( storage(ng) % SworkD(lworkd(ng)) )
        dmem(ng)=dmem(ng)+real(lworkd(ng),r8)
 
        allocate ( storage(ng) % resid(mstr(ng):mend(ng)) )
        dmem(ng)=dmem(ng)+real(mend(ng)-mstr(ng),r8)
 
# ifdef STOCHASTIC_OPT
        allocate ( storage(ng) % ad_Work(maxval(mstate)) )
        dmem(ng)=dmem(ng)+real(maxval(mstate),r8)
# endif
 
# if defined STOCHASTIC_OPT && defined HESSIAN_SO
        allocate ( storage(ng) % my_state(nstr(ng):nend(ng)) )
        dmem(ng)=dmem(ng)+real(nend(ng)-nstr(ng),r8)
# endif
 
# ifdef SO_SEMI
        allocate ( storage(ng) % so_state(mstr(ng):mend(ng),nsemi(ng)) )
        dmem(ng)=dmem(ng)+real((mend(ng)-mstr(ng))*nsemi(ng),r8)
# endif
      END DO
!
!  Allocate other arrays.
!
      ng=1
!
      allocate ( pick(ncv,ngrids) )
      dmem(ng)=dmem(ng)+real(ncv*ngrids,r8)
 
      allocate ( norm(nev+1,ngrids) )
      dmem(ng)=dmem(ng)+real((nev+1)*ngrids,r8)
 
      allocate ( rvaluer(nev+1,ngrids) )
      dmem(ng)=dmem(ng)+real((nev+1)*ngrids,r8)
 
      allocate ( rvaluei(nev+1,ngrids) )
      dmem(ng)=dmem(ng)+real((nev+1)*ngrids,r8)
 
      allocate ( sworkev(3*ncv,ngrids) )
      dmem(ng)=dmem(ng)+3.0_r8*real(ncv,r8)
 
      allocate ( sworkl(lworkl,ngrids) )
      dmem(ng)=dmem(ng)+real(lworkl*ngrids,r8)
 
      allocate ( sworkr(maxval(mstate)) )
      dmem(ng)=dmem(ng)+real(maxval(mstate),r8)
 
# ifdef DISTRIBUTE
      allocate ( swork(maxval(mstate)) )
      dmem(ng)=dmem(ng)+real(maxval(mstate),r8)
# endif
!
!-----------------------------------------------------------------------
!  Initialize module variables.
!-----------------------------------------------------------------------
!
      iparam=0
      ipntr=0
 
      DO ng=1,ngrids
        DO j=1,ncv
          pick(j,ng) = .true.
          DO i=mstr(ng),mend(ng)
            storage(ng) % Bvec(i,j) = inival
          END DO
        END DO
        DO j=1,nev+1
          norm(j,ng) = inival
          rvaluer(j,ng) = inival
          rvaluei(j,ng) = inival
          DO i=mstr(ng),mend(ng)
            storage(ng) % Rvector(i,j) = inival
          END DO
        END DO
        DO i=mstr(ng),mend(ng)
          storage(ng) % resid(i) = inival
        END DO
# ifdef STOCHASTIC_OPT
        DO i=1,maxval(mstate)
          storage(ng) % ad_Work(i) = inival
        END DO
# endif
# if defined STOCHASTIC_OPT && defined HESSIAN_SO
        DO i=nstr(ng),nend(ng)
          storage(ng) % my_state(i) = inival
        END DO
# endif
# ifdef SO_SEMI
        DO j=1,nsemi(ng)
          DO i=mstr(ng),mend(ng)
            storage(ng) % so_state(i,j) = inival
          END DO
        END DO
# endif
      END DO
      DO i=1,maxval(mstate)
        sworkr(i) = inival
      END DO
# ifdef DISTRIBUTE
      DO i=1,maxval(mstate)
        swork(i) = inival
      END DO
# endif
!
      DO ng=1,ngrids
        DO i=1,lworkd(ng)
          storage(ng) % SworkD(i) = inival
        END DO
        DO i=1,3*ncv
          sworkev(i,ng) = inival
        END DO
        DO i=1,lworkl
          sworkl(i,ng) = inival
        END DO
      END DO
!
      RETURN
      END SUBROUTINE allocate_storage
!
      SUBROUTINE deallocate_storage (ng)
!
!=======================================================================
!                                                                      !
!  This routine deallocates and initialize module variables.           !
!                                                                      !
!=======================================================================
!
      USE mod_param, ONLY : ngrids
!
!  Imported variable declarations
!
      integer, intent(in) :: ng
 
# ifdef SUBOBJECT_DEALLOCATION
!
!-----------------------------------------------------------------------
!  Deallocate each variable in the derived-type T_STORAGE structure
!  separately.
!-----------------------------------------------------------------------
!
      IF (associated(storage(ng)%Bvec)) THEN
        deallocate ( storage(ng)%Bvec )
      END IF
 
      IF (associated(storage(ng)%Rvector)) THEN
        deallocate ( storage(ng)%Rvector )
      END IF
 
      IF (associated(storage(ng)%SworkD)) THEN
        deallocate ( storage(ng)%SworkD )
      END IF
 
      IF (associated(storage(ng)%resid)) THEN
        deallocate ( storage(ng)%resid )
      END IF
 
#  ifdef STOCHASTIC_OPT
      IF (associated(storage(ng)%ad_Work)) THEN
        deallocate ( storage(ng)%ad_Work )
      END IF
#  endif
 
#  if defined STOCHASTIC_OPT && defined HESSIAN_SO
      IF (associated(storage(ng)%my_state)) THEN
        deallocate ( storage(ng)%my_state )
      END IF
#  endif
 
#  ifdef SO_SEMI
      IF (associated(storage(ng)%so_state)) THEN
        deallocate ( storage(ng)%so_state )
      END IF
#  endif
# endif
!
!-----------------------------------------------------------------------
!  Deallocate derived-type STORAGE structure.
!-----------------------------------------------------------------------
!
      IF (ng.eq.ngrids) THEN
        IF (allocated(storage)) deallocate ( storage )
      END IF
!
!-----------------------------------------------------------------------
!  Deallocate other variables in module.
!-----------------------------------------------------------------------
!
      IF (allocated(mstr))    deallocate ( mstr )
      IF (allocated(mend))    deallocate ( mend )
      IF (allocated(ido))     deallocate ( ido )
      IF (allocated(info))    deallocate ( info )
      IF (allocated(lworkd))  deallocate ( lworkd )
      IF (allocated(iparam))  deallocate ( iparam )
      IF (allocated(ipntr))   deallocate ( ipntr )
      IF (allocated(pick))    deallocate ( pick )
      IF (allocated(norm))    deallocate ( norm )
      IF (allocated(rvaluer)) deallocate ( rvaluer )
      IF (allocated(rvaluei)) deallocate ( rvaluei )
      IF (allocated(sworkev)) deallocate ( sworkev )
      IF (allocated(sworkl))  deallocate ( sworkl )
      IF (allocated(sworkr))  deallocate ( sworkr )
# ifdef DISTRIBUTE
      IF (allocated(swork))   deallocate ( swork )
# endif
!
      RETURN
      END SUBROUTINE deallocate_storage
#endif
      END MODULE mod_storage