array_modes.F

Go to the documentation of this file.

#include "cppdefs.h"
      MODULE array_modes_mod
#if defined WEAK_CONSTRAINT && (defined ARRAY_MODES || defined CLIPPING)
!
!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                                               !
!=======================================================================
!                                                                      !
!  These routines are used to process the requested eigenvector(s)     !
!  of the  stabilized representer matrix  when computing the array     !
!  modes or clipped spectrum analysis.                                 !
!                                                                      !
!  Reference:                                                          !
!                                                                      !
!  Bennett, A.F, 2002: Inverse Modeling of the Ocean and Atmosphere,   !
!     Cambridge University Press, p 49-51.                             !
!                                                                      !
!======================================================================!
!
      implicit none
 
      PRIVATE
# ifdef ARRAY_MODES
      PUBLIC :: rep_check
      PUBLIC :: rep_eigen
# endif
# ifdef CLIPPING
      PUBLIC :: rep_clip
# endif
 
      CONTAINS
 
# ifdef ARRAY_MODES
      SUBROUTINE rep_check (ng, model, outLoop, NinnLoop)
!
!=======================================================================
!                                                                      !
!  This routine checks the dot-product of the array mode with the      !
!  corresponding eigenvector of the stabilized representer matrix.     !
!                                                                      !
!=======================================================================
!
      USE mod_param
      USE mod_parallel
      USE mod_fourdvar
      USE mod_iounits
      USE mod_scalars
 
#  ifdef DISTRIBUTE
!
      USE distribute_mod, ONLY : mp_bcastf
#  endif
!
!  Imported variable declarations
!
      integer, intent(in) :: ng, model, outloop, ninnloop
!
!  Local variable declarations.
!
      integer :: iobs, innloop
 
      real(r8) :: zsum
#  ifdef RPCG
      real(r8) :: zfact
      real(r8), dimension(NinnLoop)  :: zdot
#  endif
!
!-----------------------------------------------------------------------
!  Check the dot-product of the array mode with the corresponding
!  eigenvector of the stabilized representer matrix.
!-----------------------------------------------------------------------
!
      master_thread : IF (master) THEN
 
#  ifdef RPCG
        DO innloop=1,ninnloop
          zdot(innloop)=0.0_r8
          IF (innloop.eq.1) THEN
            zfact=1.0_r8/cg_gnorm_v(outloop)
          ELSE
            zfact=1.0_r8/cg_beta(innloop,outloop)
          END IF
          DO iobs=1,ndatum(ng)
            IF (obserr(iobs).ne.0.0_r8) THEN
              zdot(innloop)=zdot(innloop)+                              &
     &                      tlmodval(iobs)*zfact*                       &
     &                      tlmodval_s(iobs,innloop,outloop)/           &
     &                      obserr(iobs)
            END IF
          END DO
        END DO
!
        zsum=0.0_r8
        DO innloop=1,ninnloop
          zsum=zsum+zdot(innloop)*cg_zv(innloop,nvct,outloop)
        END DO
#  else
        DO iobs=1,ndatum(ng)
          admodval(iobs)=0.0_r8
        END DO
!
!  Multiply desired eigenvector of Lanczos tridiagonal matrix
!  by the Lanczos vectors to obtain the corresponding eigenvector
!  of the preconditioned stabilized representer matrix.
!
        DO iobs=1,ndatum(ng)
          DO innloop=1,ninnloop
            admodval(iobs)=admodval(iobs)+                              &
     &                     cg_zv(innloop,nvct,outloop)*                 &
     &                     zcglwk(iobs,innloop,outloop)
          END DO
        END DO
!
!  Now convert ADmodVal back to physical units.
!
        DO iobs=1,ndatum(ng)
          IF (obserr(iobs).ne.0.0_r8) THEN
            admodval(iobs)=admodval(iobs)/sqrt(obserr(iobs))
          END IF
        END DO
!
!  Now compute the dot-product of the final solution with the initial
!  eigenvector.
!
        zsum=0.0_r8
        DO iobs=1,ndatum(ng)
          zsum=zsum+tlmodval(iobs)*admodval(iobs)
        END DO
#  endif
!
!  Compare the dot-product with (cg_Ritz-1).
!
        WRITE (stdout,10) zsum, cg_ritz(nvct,outloop)-1
 10     FORMAT (/,' REP CHECK: zsum = ', 1p, e14.7,2x,                  &
     &          'cg_Ritz-1 = ', 1p, e14.7)
 
      END IF master_thread
 
#  ifdef DISTRIBUTE
!
!-----------------------------------------------------------------------
!  Broadcast new solution to other nodes.
!-----------------------------------------------------------------------
!
      CALL mp_bcastf (ng, model, admodval)
#  endif
 
      RETURN
      END SUBROUTINE rep_check
 
      SUBROUTINE rep_eigen (ng, model, outLoop, NinnLoop)
!
!=======================================================================
!                                                                      !
!  This routine computes the specified eigenvector, Nvct, of the       !
!  stabilized representer matrix.                                      !
!                                                                      !
!=======================================================================
!
      USE mod_param
      USE mod_parallel
      USE mod_fourdvar
      USE mod_iounits
      USE mod_scalars
 
#  ifdef DISTRIBUTE
!
      USE distribute_mod, ONLY : mp_bcastf
#  endif
!
!  Imported variable declarations
!
      integer, intent(in) :: ng, model, outloop, ninnloop
!
!  Local variable declarations.
!
      integer :: iobs, innloop
!
!-----------------------------------------------------------------------
!  Compute specified eignevector of the stabilized representer matrix.
!-----------------------------------------------------------------------
!
      master_thread : IF (master) THEN
 
        DO iobs=1,ndatum(ng)
          admodval(iobs)=0.0_r8
        END DO
!
!  Multiply desired eigenvector of Lanczos tridiagonal matrix
!  by the Lanczos vectors to obtain the corresponding eigenvector
!  of the preconditioned stabilized representer matrix.
!
        DO iobs=1,ndatum(ng)
          DO innloop=1,ninnloop
            admodval(iobs)=admodval(iobs)+                              &
     &                     cg_zv(innloop,nvct,outloop)*                 &
     &                     zcglwk(iobs,innloop,outloop)
          END DO
        END DO
 
#  ifndef RPCG
!
!  Now convert ADmodVal back to physical units.
!
        DO iobs=1,ndatum(ng)
          IF (obserr(iobs).ne.0.0_r8) THEN
            admodval(iobs)=admodval(iobs)/sqrt(obserr(iobs))
          END IF
        END DO
#  endif
 
      END IF master_thread
 
#  ifdef DISTRIBUTE
!
!-----------------------------------------------------------------------
!  Broadcast new solution to other nodes.
!-----------------------------------------------------------------------
!
      CALL mp_bcastf (ng, model, admodval)
#  endif
 
      RETURN
      END SUBROUTINE rep_eigen
# endif
 
# ifdef CLIPPING
      SUBROUTINE rep_clip (ng, model, outLoop, NinnLoop)
!
!=======================================================================
!                                                                      !
!  This routine performs clipping of the analysis by disgarding        !
!  potentially unphysical array modes.                                 !
!                                                                      !
!=======================================================================
!
      USE mod_param
      USE mod_parallel
      USE mod_fourdvar
      USE mod_iounits
      USE mod_scalars
 
#  ifdef DISTRIBUTE
!
      USE distribute_mod, ONLY : mp_bcastf
#  endif
!
!  Imported variable declarations
!
      integer, intent(in) :: ng, model, outloop, ninnloop
!
!  Local variable declarations.
!
      integer :: iobs, ivec, innloop
 
      real(r8), dimension(NinnLoop) :: zu
 
      real(r8), dimension(Ndatum(ng)) :: innov, rsvec
!
!-----------------------------------------------------------------------
!  Clipp the analysis by disgarding potentially unphysical array modes.
!-----------------------------------------------------------------------
!
      master_thread : IF (master) THEN
!
!  First compute the dot-product of innovation vector with each
!  selected eigenvector of the stabilized representer matrix.
!  All eigenvectors < Nvct are disgarded.
!
        DO iobs=1,ndatum(ng)
          innov(iobs)=obsval(iobs)-nlmodval(iobs)
        END DO
        DO ivec=nvct,ninner
          DO iobs=1,ndatum(ng)
            rsvec(iobs)=0.0_r8
          END DO
          DO iobs=1,ndatum(ng)
            DO innloop=1,ninnloop
              rsvec(iobs)=rsvec(iobs)+                                  &
     &                    cg_zv(innloop,ivec,outloop)*                  &
     &                    zcglwk(iobs,innloop,outloop)
            END DO
          END DO
!
!  Convert RSVEC back to physical units.
!
          DO iobs=1,ndatum(ng)
            IF (obserr(iobs).ne.0.0_r8) THEN
              rsvec(iobs)=rsvec(iobs)/sqrt(obserr(iobs))
            END IF
          END DO
          zu(ivec)=0.0_r8
          DO iobs=1,ndatum(ng)
            zu(ivec)=zu(ivec)+innov(iobs)*rsvec(iobs)
          END DO
        END DO
!
!  Second divide by the eigenvalues of the stabilized representer
!  matrix.
!
        DO ivec=nvct,ninner
          zu(ivec)=zu(ivec)/cg_ritz(ivec,outloop)
        END DO
!
!  Finally form the weighted sum of the selected eigenvectors of the
!  stabilized representer matrix.
!
        DO iobs=1,ndatum(ng)
          admodval(iobs)=0.0_r8
        END DO
!
!  Multiply desired eigenvector of Lanczos tridiagonal matrix
!  by the Lanczos vectors to obtain the corresponding eigenvector
!  of the preconditioned stabilized representer matrix.
!
        DO ivec=nvct,ninner
          DO iobs=1,ndatum(ng)
            DO innloop=1,ninnloop
              admodval(iobs)=admodval(iobs)+                            &
     &                       cg_zv(innloop,ivec,outloop)*               &
     &                       zcglwk(iobs,innloop,outloop)*              &
     &                       zu(ivec)
            END DO
          END DO
        END DO
!
!  Now convert ADmodVal back to physical units.
!
        DO iobs=1,ndatum(ng)
          IF (obserr(iobs).ne.0.0_r8) THEN
            admodval(iobs)=admodval(iobs)/sqrt(obserr(iobs))
          END IF
        END DO
 
      END IF master_thread
 
#  ifdef DISTRIBUTE
!
!-----------------------------------------------------------------------
!  Broadcast new solution to other nodes.
!-----------------------------------------------------------------------
!
      CALL mp_bcastf (ng, model, admodval)
#  endif
 
      RETURN
      END SUBROUTINE rep_clip
# endif
#endif
      END MODULE array_modes_mod