C $Header: /u/gcmpack/MITgcm/pkg/gmredi/gmredi_calc_eigs.F,v 1.7 2014/11/24 10:23:50 dfer Exp $
C $Name: $
#include "GMREDI_OPTIONS.h"
C !ROUTINE: EIGENVAL
C !INTERFACE:
SUBROUTINE GMREDI_CALC_EIGS(
I iMin, iMax, jMin, jMax,
I bi, bj, N2, myThid, kLow,
I mask, hfac, recip_hfac,
I rlow, nmodes, writediag,
O Rmid, vec)
C !DESCRIPTION: \bv
C *==========================================================*
C | SUBROUTINE GMREDI_CALC_URMS
C | o Calculate the vertical structure of the rms eddy
C | velocity based on baroclinic modal decomposition
C *==========================================================*
C \ev
IMPLICIT NONE
C == Global variables ==
#include "SIZE.h"
#include "GRID.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GMREDI.h"
C !INPUT/OUTPUT PARAMETERS:
C == Routine arguments ==
C bi, bj :: tile indices
LOGICAL writediag
INTEGER iMin,iMax,jMin,jMax
INTEGER bi, bj
INTEGER myThid
INTEGER nmodes
INTEGER kLow(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
_RL mask(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
_RL hfac(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
_RL recip_hfac(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
_RL rlow(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
_RL N2( 1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
_RL Rmid(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
_RL vec(nmodes,1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
#ifdef GM_K3D
C !LOCAL VARIABLES:
C == Local variables ==
INTEGER i,j,k,kk,m
# ifdef HAVE_LAPACK
C info :: error code from LAPACK
C idx :: index used for sorting the eigenvalues
C a3d :: lower diagonal of eigenvalue problem
C b3d :: diagonal of eigenvalue problem
C c3d :: upper diagonal of eigenvalue problem
C val :: Eigenvalue (wavenumber) of the first baroclinic mode
C eigR :: Real component of all the eigenvalues in a water column
C eigI :: Imaginary component of all the eigenvalues in a water column
C vecs :: All the eigenvectors of a water column
C dummy :: Redundant variable for calling lapack
C work :: Work array for lapack
C array :: Array containing the matrix with a,b,c
C eigval :: Vector containing the eigenvalues
INTEGER info
INTEGER idx
_RL a3d( 1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
_RL b3d( 1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
_RL c3d( 1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
_RL val( 1-Olx:sNx+Olx,1-Oly:sNy+Oly)
_RL eigR(Nr),eigI(Nr),vecs(Nr,Nr),dummy(1,Nr),work(Nr*Nr)
_RL array(Nr,Nr)
_RL eigval(0:nmodes)
else
C drNr :: distance from bottom of cell to cell centre at Nr
C BuoyFreq :: buoyancy frequency, SQRT(N2)
C intN :: Vertical integral of BuoyFreq to each grid cell centre
C intN0 :: Vertical integral of BuoyFreq to z=0
C c1 :: intN0/pi
C nEigs :: number of eigenvalues/vectors to calculate
_RL drNr
_RL BuoyFreq(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr+1)
_RL intN(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
_RL intN0(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
_RL c1(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
INTEGER nEigs(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
# endif
C small :: a small number (used to avoid floating point exceptions)
C vecint :: vertical integral of eigenvector and/or square of eigenvector
C fCori2 :: square of the Coriolis parameter
_RL small
_RL vecint(nmodes,1-Olx:sNx+Olx,1-Oly:sNy+Oly)
_RL fCori2(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
CHARACTER*(MAX_LEN_MBUF) msgBuf
small = TINY(zeroRL)
C Square of the Coriolis parameter
DO i=1-OLx,sNx+OLx
DO j=1-OLy,sNy+OLy
fCori2(i,j) = fCori(i,j,bi,bj)*fCori(i,j,bi,bj)
ENDDO
ENDDO
DO k=1,Nr
DO j=1-Oly,sNy+Oly
DO i=1-Olx,sNx+Olx
DO m=1,nmodes
vec(m,i,j,k) = zeroRL
ENDDO
ENDDO
ENDDO
ENDDO
# ifdef HAVE_LAPACK
C Calculate the tridiagonal operator matrix for
C f^2 d/dz 1/N^2 d/dz
C a3d is the lower off-diagonal, b3d is the diagonal
C and c3d is the upper off-diagonal
DO k=1,Nr
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
IF (kLow(i,j) .GT. 0) THEN
IF (k.EQ.1) THEN
a3d(i,j,k) = zeroRL
c3d(i,j,k) = fCori2(i,j)*recip_hFac(i,j,k)
& *recip_drC(k+1)*recip_drF(k)/N2(i,j,k+1)
b3d(i,j,k) = -c3d(i,j,k)
ELSEIF (k.LT.kLow(i,j)) THEN
a3d(i,j,k) = fCori2(i,j)*recip_hFac(i,j,k)
& *recip_drF(k)*recip_drC(k)/N2(i,j,k)
c3d(i,j,k) = fCori2(i,j)*recip_hFac(i,j,k)
& *recip_drF(k)*recip_drC(k+1)/N2(i,j,k+1)
b3d(i,j,k) = -a3d(i,j,k)-c3d(i,j,k)
ELSEIF (k.EQ.kLow(i,j)) THEN
a3d(i,j,k) = fCori2(i,j)*recip_hFac(i,j,k)
& *recip_drF(k)*recip_drC(k)/N2(i,j,k)
c3d(i,j,k) = zeroRL
b3d(i,j,k) = -a3d(i,j,k)
ELSE
a3d(i,j,k) = zeroRL
b3d(i,j,k) = zeroRL
c3d(i,j,k) = zeroRL
ENDIF
ELSE
a3d(i,j,k) = zeroRL
b3d(i,j,k) = zeroRL
c3d(i,j,k) = zeroRL
ENDIF
ENDDO
ENDDO
ENDDO
DO j=1-Oly,sNy+Oly
DO i=1-Olx,sNx+Olx
IF (kLow(i,j).GT.0) THEN
DO kk=1,Nr
DO k=1,Nr
array(k,kk) = zeroRL
ENDDO
ENDDO
k=1
array(k,k) = b3d(i,j,k)
array(k,k+1) = c3d(i,j,k)
DO k=2,Nr-1
array(k,k-1) = a3d(i,j,k)
array(k,k) = b3d(i,j,k)
array(k,k+1) = c3d(i,j,k)
ENDDO
k=Nr
array(k,k-1) = a3d(i,j,k)
array(k,k) = b3d(i,j,k)
CALL DGEEV('N','V',Nr,array,Nr,eigR,eigI,dummy,1,
& vecs,Nr,work,Nr*Nr,info)
IF( info.LT.0 ) THEN
WRITE(msgBuf,'(A,x,2(A1,I2),A1,x,A,I4)')
& 'GMREDI_CALC_EIGS problem with arguments for DGEEV at',
& '(',i,',',j,')', 'error code =',info
CALL PRINT_ERROR( msgBuf , myThid )
ELSEIF(info.GT.0 ) THEN
WRITE(msgBuf,'(A,x,2(A1,I2),A1,x,A,I4)')
& 'GMREDI_CALC_EIGS problems with eigensolver DGEEV at',
& '(',i,',',j,')', 'error code =',info
CALL PRINT_ERROR( msgBuf , myThid )
ENDIF
C Find the second largest eigenvalue (the Rossby radius)
C and the first M baroclinic modes (eigenvectors)
DO m=0,nmodes
eigval(m) = -HUGE(zeroRL)
ENDDO
DO k=1,kLow(i,j)
eigval(0) = MAX(eigval(0),eigR(k))
ENDDO
DO m=1,MIN(nmodes,klow(i,j)-1)
DO k=1,kLow(i,j)
IF (eigR(k).LT.eigval(m-1)) THEN
eigval(m) = MAX(eigval(m),eigR(k))
IF (eigval(m).EQ.eigR(k)) idx=k
ENDIF
ENDDO
IF(vecs(1,idx).LT.zeroRL) THEN
DO k=1,Nr
vec(m,i,j,k) = -vecs(k,idx)
ENDDO
ELSE
DO k=1,Nr
vec(m,i,j,k) = vecs(k,idx)
ENDDO
ENDIF
ENDDO
val(i,j) = eigval(1)
ELSE
val(i,j)=zeroRL
DO k=1,Nr
DO m=1,nmodes
vec(m,i,j,k)=zeroRL
ENDDO
ENDDO
ENDIF
ENDDO
ENDDO
DO j=1-Oly,sNy+Oly
DO i=1-Olx,sNx+Olx
IF (kLow(i,j).GT.2 .AND. val(i,j).NE.zeroRL) THEN
Rmid(i,j) = 1.0/(SQRT(ABS(val(i,j)))+small)
ELSE
Rmid(i,j) = zeroRL
ENDIF
ENDDO
ENDDO
else
DO k=1,Nr
DO j=1-Oly,sNy+Oly
DO i=1-Olx,sNx+Olx
BuoyFreq(i,j,k) = mask(i,j,k)*SQRT(N2(i,j,k))
ENDDO
ENDDO
ENDDO
k=Nr+1
DO j=1-Oly,sNy+Oly
DO i=1-Olx,sNx+Olx
BuoyFreq(i,j,k) = zeroRL
ENDDO
ENDDO
C integrate N using something like Simpson s rule (but not quite)
C drC*( (N(k+1)+N(k+2))/2 + (N(k)+N(k+1))/2 )/2
C when k=Nr, say that N(k+2)=0 and N(k+1)=0
k=Nr
C drNr is like drC(Nr+1)/2
drNr = rC(Nr)-rF(Nr+1)
DO j=1-Oly,sNy+Oly
DO i=1-Olx,sNx+Olx
intN(i,j,k) = op5*BuoyFreq(i,j,k)*drNr
ENDDO
ENDDO
DO k=Nr-1,1,-1
DO j=1-Oly,sNy+Oly
DO i=1-Olx,sNx+Olx
intN(i,j,k) = intN(i,j,k+1)
& + drC(k)*( op25*BuoyFreq(i,j,k+2) + op5*BuoyFreq(i,j,k)
& + op25*BuoyFreq(i,j,k+1) )
ENDDO
ENDDO
ENDDO
C intN integrates to z=rC(1). We want to integrate to z=0.
C Assume that N(0)=0 and N(1)=0.
C drC(1) is like half a grid cell.
DO j=1-Oly,sNy+Oly
DO i=1-Olx,sNx+Olx
intN0(i,j) = intN(i,j,1)
& + drC(1)*op5*BuoyFreq(i,j,2)
ENDDO
ENDDO
DO j=1-Oly,sNy+Oly
DO i=1-Olx,sNx+Olx
c1(i,j) = intN0(i,j)/pi
Rmid(i,j) = c1(i,j)/ABS(fCori(i,j,bi,bj))
ENDDO
ENDDO
DO j=1-Oly,sNy+Oly
DO i=1-Olx,sNx+Olx
nEigs(i,j) = MIN(klow(i,j),nmodes)
ENDDO
ENDDO
DO k=1,Nr
DO j=1-Oly,sNy+Oly
DO i=1-Olx,sNx+Olx
IF (mask(i,j,k).NE.0.0) THEN
DO m=1,nEigs(i,j)
vec(m,i,j,k) = -COS(intN(i,j,k)/(m*c1(i,j)))
ENDDO
ENDIF
ENDDO
ENDDO
ENDDO
C The WKB approximation for the baroclinic mode does not always
C integrate to zero so we adjust it.
DO j=1-Oly,sNy+Oly
DO i=1-Olx,sNx+Olx
DO m=1,nEigs(i,j)
vecint(m,i,j) = zeroRL
ENDDO
ENDDO
ENDDO
DO k=1,Nr
DO j=1-Oly,sNy+Oly
DO i=1-Olx,sNx+Olx
DO m=1,nEigs(i,j)
vecint(m,i,j) = vecint(m,i,j) + hfac(i,j,k)*vec(m,i,j,k)*drF(k)
ENDDO
ENDDO
ENDDO
ENDDO
DO j=1-Oly,sNy+Oly
DO i=1-Olx,sNx+Olx
DO m=1,nEigs(i,j)
vecint(m,i,j) = vecint(m,i,j)/(-rlow(i,j)+small)
ENDDO
ENDDO
ENDDO
DO k=1,Nr
DO j=1-Oly,sNy+Oly
DO i=1-Olx,sNx+Olx
DO m=1,nEigs(i,j)
vec(m,i,j,k) = vec(m,i,j,k) - vecint(m,i,j)
ENDDO
ENDDO
ENDDO
ENDDO
# endif
C Normalise the eigenvectors
DO j=1-Oly,sNy+Oly
DO i=1-Olx,sNx+Olx
DO m=1,nmodes
vecint(m,i,j) = zeroRL
ENDDO
ENDDO
ENDDO
DO k=1,Nr
DO j=1-Oly,sNy+Oly
DO i=1-Olx,sNx+Olx
DO m=1,nmodes
vecint(m,i,j) = vecint(m,i,j) +
& mask(i,j,k)*drF(k)*hfac(i,j,k)
& *vec(m,i,j,k)*vec(m,i,j,k)
ENDDO
ENDDO
ENDDO
ENDDO
DO j=1-Oly,sNy+Oly
DO i=1-Olx,sNx+Olx
DO m=1,nmodes
vecint(m,i,j) = SQRT(vecint(m,i,j))
ENDDO
ENDDO
ENDDO
DO k=1,Nr
DO j=1-Oly,sNy+Oly
DO i=1-Olx,sNx+Olx
DO m=1,nmodes
vec(m,i,j,k) = vec(m,i,j,k)/(vecint(m,i,j)+small)
ENDDO
ENDDO
ENDDO
ENDDO
# ifdef ALLOW_DIAGNOSTICS
C Diagnostics
IF ( useDiagnostics.AND.writediag ) THEN
# ifdef HAVE_LAPACK
CALL DIAGNOSTICS_FILL(a3d, 'GM_A3D ',0,Nr,0,1,1,myThid)
CALL DIAGNOSTICS_FILL(b3d, 'GM_B3D ',0,Nr,0,1,1,myThid)
CALL DIAGNOSTICS_FILL(c3d, 'GM_C3D ',0,Nr,0,1,1,myThid)
# endif
ENDIF
# endif
#endif /* GM_K3D */
RETURN
END