Back to home page

MITgcm
 
 

    

File indexing completed on 2022-01-29 06:11:22 UTC

view on githubraw file Latest commit 4578baf3 on 2021-12-13 15:21:55 UTC
0d1e4b948d Mich*0001 #include "GMREDI_OPTIONS.h"
f3b9070ca3 Jean*0002 
0d1e4b948d Mich*0003 C     !ROUTINE: EIGENVAL
                0004 C     !INTERFACE:
                0005       SUBROUTINE GMREDI_CALC_EIGS(
                0006      I     iMin, iMax, jMin, jMax,
                0007      I     bi, bj, N2, myThid, kLow,
                0008      I     mask, hfac, recip_hfac,
5a6ef5c2b4 Mich*0009      I     rlow, nmodes, writediag,
                0010      O     Rmid, vec)
0d1e4b948d Mich*0011 
                0012 C     !DESCRIPTION: \bv
                0013 C     *==========================================================*
                0014 C     | SUBROUTINE GMREDI_CALC_URMS
f3b9070ca3 Jean*0015 C     | o Calculate the vertical structure of the rms eddy
0d1e4b948d Mich*0016 C     |   velocity based on baroclinic modal decomposition
                0017 C     *==========================================================*
                0018 C     \ev
                0019 
                0020       IMPLICIT NONE
                0021 
                0022 C     == Global variables ==
                0023 #include "SIZE.h"
                0024 #include "GRID.h"
                0025 #include "EEPARAMS.h"
                0026 #include "PARAMS.h"
                0027 #include "GMREDI.h"
                0028 
                0029 C     !INPUT/OUTPUT PARAMETERS:
                0030 C     == Routine arguments ==
                0031 C     bi, bj    :: tile indices
                0032       LOGICAL writediag
                0033       INTEGER iMin,iMax,jMin,jMax
                0034       INTEGER bi, bj
                0035       INTEGER myThid
5a6ef5c2b4 Mich*0036       INTEGER nmodes
4578baf364 Jean*0037       INTEGER kLow(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0038       _RL mask(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0039       _RL hfac(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0040       _RL recip_hfac(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0041       _RL rlow(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0042       _RL N2(  1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0043       _RL Rmid(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0044       _RL vec(nmodes,1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
f3b9070ca3 Jean*0045 
05118ac017 Jean*0046 #ifdef GM_BATES_K3D
0d1e4b948d Mich*0047 C     !LOCAL VARIABLES:
                0048 C     == Local variables ==
5a6ef5c2b4 Mich*0049       INTEGER i,j,k,kk,m
dc60433548 Mich*0050 # ifdef HAVE_LAPACK
5a6ef5c2b4 Mich*0051 C     info     :: error code from LAPACK
                0052 C     idx      :: index used for sorting the eigenvalues
                0053 C     a3d      :: lower diagonal of eigenvalue problem
                0054 C     b3d      :: diagonal of eigenvalue problem
                0055 C     c3d      :: upper diagonal of eigenvalue problem
                0056 C     val      :: Eigenvalue (wavenumber) of the first baroclinic mode
                0057 C     eigR     :: Real component of all the eigenvalues in a water column
                0058 C     eigI     :: Imaginary component of all the eigenvalues in a water column
                0059 C     vecs     :: All the eigenvectors of a water column
                0060 C     dummy    :: Redundant variable for calling lapack
                0061 C     work     :: Work array for lapack
                0062 C     array    :: Array containing the matrix with a,b,c
                0063 C     eigval   :: Vector containing the eigenvalues
                0064       INTEGER info
                0065       INTEGER idx
4578baf364 Jean*0066       _RL a3d(   1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0067       _RL b3d(   1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0068       _RL c3d(   1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0069       _RL val(   1-OLx:sNx+OLx,1-OLy:sNy+OLy)
0d1e4b948d Mich*0070       _RL eigR(Nr),eigI(Nr),vecs(Nr,Nr),dummy(1,Nr),work(Nr*Nr)
                0071       _RL array(Nr,Nr)
5a6ef5c2b4 Mich*0072       _RL eigval(0:nmodes)
                0073 # else
                0074 C     drNr     :: distance from bottom of cell to cell centre at Nr
                0075 C     BuoyFreq :: buoyancy frequency, SQRT(N2)
                0076 C     intN     :: Vertical integral of BuoyFreq to each grid cell centre
                0077 C     intN0    :: Vertical integral of BuoyFreq to z=0
                0078 C     c1       :: intN0/pi
                0079 C     nEigs    :: number of eigenvalues/vectors to calculate
                0080       _RL drNr
4578baf364 Jean*0081       _RL BuoyFreq(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr+1)
                0082       _RL intN(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0083       _RL intN0(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0084       _RL c1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0085       INTEGER nEigs(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
5a6ef5c2b4 Mich*0086 # endif
                0087 C     small    :: a small number (used to avoid floating point exceptions)
                0088 C     vecint   :: vertical integral of eigenvector and/or square of eigenvector
                0089 C     fCori2   :: square of the Coriolis parameter
                0090       _RL small
4578baf364 Jean*0091       _RL vecint(nmodes,1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0092       _RL fCori2(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
5a6ef5c2b4 Mich*0093       CHARACTER*(MAX_LEN_MBUF) msgBuf
                0094 
0d1e4b948d Mich*0095       small = TINY(zeroRL)
                0096 
                0097 C     Square of the Coriolis parameter
f3b9070ca3 Jean*0098       DO i=1-OLx,sNx+OLx
                0099        DO j=1-OLy,sNy+OLy
0d1e4b948d Mich*0100         fCori2(i,j) = fCori(i,j,bi,bj)*fCori(i,j,bi,bj)
                0101        ENDDO
                0102       ENDDO
                0103 
                0104       DO k=1,Nr
4578baf364 Jean*0105        DO j=1-OLy,sNy+OLy
                0106         DO i=1-OLx,sNx+OLx
5a6ef5c2b4 Mich*0107          DO m=1,nmodes
0d1e4b948d Mich*0108           vec(m,i,j,k) = zeroRL
                0109          ENDDO
                0110         ENDDO
                0111        ENDDO
                0112       ENDDO
                0113 
dc60433548 Mich*0114 # ifdef HAVE_LAPACK
0d1e4b948d Mich*0115 C     Calculate the tridiagonal operator matrix for
                0116 C     f^2 d/dz 1/N^2 d/dz
                0117 C     a3d is the lower off-diagonal, b3d is the diagonal
                0118 C     and c3d is the upper off-diagonal
                0119       DO k=1,Nr
f3b9070ca3 Jean*0120        DO j=1-OLy,sNy+OLy
                0121         DO i=1-OLx,sNx+OLx
0d1e4b948d Mich*0122          IF (kLow(i,j) .GT. 0) THEN
                0123            IF (k.EQ.1) THEN
                0124              a3d(i,j,k) = zeroRL
                0125              c3d(i,j,k) = fCori2(i,j)*recip_hFac(i,j,k)
                0126      &                    *recip_drC(k+1)*recip_drF(k)/N2(i,j,k+1)
                0127              b3d(i,j,k) = -c3d(i,j,k)
                0128 
                0129            ELSEIF (k.LT.kLow(i,j)) THEN
                0130              a3d(i,j,k) = fCori2(i,j)*recip_hFac(i,j,k)
                0131      &                    *recip_drF(k)*recip_drC(k)/N2(i,j,k)
                0132              c3d(i,j,k) = fCori2(i,j)*recip_hFac(i,j,k)
                0133      &                    *recip_drF(k)*recip_drC(k+1)/N2(i,j,k+1)
                0134              b3d(i,j,k) = -a3d(i,j,k)-c3d(i,j,k)
                0135 
                0136            ELSEIF (k.EQ.kLow(i,j)) THEN
                0137              a3d(i,j,k) = fCori2(i,j)*recip_hFac(i,j,k)
                0138      &                    *recip_drF(k)*recip_drC(k)/N2(i,j,k)
                0139              c3d(i,j,k) = zeroRL
                0140              b3d(i,j,k) = -a3d(i,j,k)
                0141 
                0142            ELSE
                0143              a3d(i,j,k) = zeroRL
                0144              b3d(i,j,k) = zeroRL
                0145              c3d(i,j,k) = zeroRL
                0146            ENDIF
                0147 
                0148          ELSE
                0149            a3d(i,j,k) = zeroRL
                0150            b3d(i,j,k) = zeroRL
                0151            c3d(i,j,k) = zeroRL
                0152          ENDIF
                0153         ENDDO
                0154        ENDDO
                0155       ENDDO
                0156 
4578baf364 Jean*0157       DO j=1-OLy,sNy+OLy
                0158        DO i=1-OLx,sNx+OLx
0d1e4b948d Mich*0159         IF (kLow(i,j).GT.0) THEN
                0160           DO kk=1,Nr
                0161            DO k=1,Nr
                0162             array(k,kk) = zeroRL
                0163            ENDDO
                0164           ENDDO
f3b9070ca3 Jean*0165 
0d1e4b948d Mich*0166           k=1
                0167           array(k,k)   = b3d(i,j,k)
                0168           array(k,k+1) = c3d(i,j,k)
                0169           DO k=2,Nr-1
                0170            array(k,k-1) = a3d(i,j,k)
                0171            array(k,k)   = b3d(i,j,k)
                0172            array(k,k+1) = c3d(i,j,k)
                0173           ENDDO
                0174           k=Nr
                0175           array(k,k-1) = a3d(i,j,k)
                0176           array(k,k)   = b3d(i,j,k)
f3b9070ca3 Jean*0177 
0d1e4b948d Mich*0178           CALL DGEEV('N','V',Nr,array,Nr,eigR,eigI,dummy,1,
                0179      &         vecs,Nr,work,Nr*Nr,info)
5a6ef5c2b4 Mich*0180           IF( info.LT.0 ) THEN
4578baf364 Jean*0181             WRITE(msgBuf,'(A,x,2(A1,I2),A1,x,A,I4)')
5a6ef5c2b4 Mich*0182      &           'GMREDI_CALC_EIGS problem with arguments for DGEEV at',
                0183      &           '(',i,',',j,')', 'error code =',info
                0184             CALL PRINT_ERROR( msgBuf , myThid )
4578baf364 Jean*0185 
5a6ef5c2b4 Mich*0186           ELSEIF(info.GT.0 ) THEN
4578baf364 Jean*0187             WRITE(msgBuf,'(A,x,2(A1,I2),A1,x,A,I4)')
5a6ef5c2b4 Mich*0188      &           'GMREDI_CALC_EIGS problems with eigensolver DGEEV at',
                0189      &           '(',i,',',j,')', 'error code =',info
                0190             CALL PRINT_ERROR( msgBuf , myThid )
4578baf364 Jean*0191 
5a6ef5c2b4 Mich*0192           ENDIF
                0193 
4578baf364 Jean*0194 C         Find the second largest eigenvalue (the Rossby radius)
0d1e4b948d Mich*0195 C         and the first M baroclinic modes (eigenvectors)
5a6ef5c2b4 Mich*0196           DO m=0,nmodes
0d1e4b948d Mich*0197            eigval(m) = -HUGE(zeroRL)
                0198           ENDDO
f3b9070ca3 Jean*0199 
0d1e4b948d Mich*0200           DO k=1,kLow(i,j)
                0201            eigval(0) = MAX(eigval(0),eigR(k))
                0202           ENDDO
219a8971b6 Mich*0203           DO m=1,MIN(nmodes,klow(i,j)-1)
0d1e4b948d Mich*0204            DO k=1,kLow(i,j)
                0205             IF (eigR(k).LT.eigval(m-1)) THEN
                0206               eigval(m) = MAX(eigval(m),eigR(k))
                0207               IF (eigval(m).EQ.eigR(k)) idx=k
                0208             ENDIF
                0209            ENDDO
                0210            IF(vecs(1,idx).LT.zeroRL) THEN
                0211              DO k=1,Nr
                0212               vec(m,i,j,k) = -vecs(k,idx)
                0213              ENDDO
                0214            ELSE
                0215              DO k=1,Nr
                0216               vec(m,i,j,k) = vecs(k,idx)
                0217              ENDDO
                0218            ENDIF
                0219           ENDDO
                0220           val(i,j) = eigval(1)
                0221         ELSE
                0222           val(i,j)=zeroRL
                0223           DO k=1,Nr
5a6ef5c2b4 Mich*0224            DO m=1,nmodes
0d1e4b948d Mich*0225             vec(m,i,j,k)=zeroRL
                0226            ENDDO
                0227           ENDDO
f3b9070ca3 Jean*0228 
0d1e4b948d Mich*0229         ENDIF
                0230        ENDDO
                0231       ENDDO
                0232 
4578baf364 Jean*0233        DO j=1-OLy,sNy+OLy
                0234         DO i=1-OLx,sNx+OLx
5a6ef5c2b4 Mich*0235          IF (kLow(i,j).GT.2 .AND. val(i,j).NE.zeroRL) THEN
                0236            Rmid(i,j) = 1.0/(SQRT(ABS(val(i,j)))+small)
                0237          ELSE
                0238            Rmid(i,j) = zeroRL
                0239          ENDIF
0d1e4b948d Mich*0240         ENDDO
                0241        ENDDO
                0242 
4578baf364 Jean*0243 # else
0d1e4b948d Mich*0244       DO k=1,Nr
4578baf364 Jean*0245        DO j=1-OLy,sNy+OLy
                0246         DO i=1-OLx,sNx+OLx
5a6ef5c2b4 Mich*0247          BuoyFreq(i,j,k) = mask(i,j,k)*SQRT(N2(i,j,k))
                0248         ENDDO
                0249        ENDDO
                0250       ENDDO
                0251       k=Nr+1
4578baf364 Jean*0252       DO j=1-OLy,sNy+OLy
                0253        DO i=1-OLx,sNx+OLx
5a6ef5c2b4 Mich*0254         BuoyFreq(i,j,k) = zeroRL
                0255        ENDDO
                0256       ENDDO
                0257 
c0612052db Jean*0258 C     integrate N using something like Simpson s rule (but not quite)
5a6ef5c2b4 Mich*0259 C     drC*( (N(k+1)+N(k+2))/2 + (N(k)+N(k+1))/2 )/2
                0260 C     when k=Nr, say that N(k+2)=0 and N(k+1)=0
                0261       k=Nr
                0262 C     drNr is like drC(Nr+1)/2
                0263       drNr = rC(Nr)-rF(Nr+1)
4578baf364 Jean*0264       DO j=1-OLy,sNy+OLy
                0265        DO i=1-OLx,sNx+OLx
5a6ef5c2b4 Mich*0266         intN(i,j,k) = op5*BuoyFreq(i,j,k)*drNr
                0267        ENDDO
                0268       ENDDO
                0269       DO k=Nr-1,1,-1
4578baf364 Jean*0270        DO j=1-OLy,sNy+OLy
                0271         DO i=1-OLx,sNx+OLx
                0272          intN(i,j,k) = intN(i,j,k+1)
                0273      &        + drC(k)*( op25*BuoyFreq(i,j,k+2) + op5*BuoyFreq(i,j,k)
5a6ef5c2b4 Mich*0274      &                 + op25*BuoyFreq(i,j,k+1) )
                0275         ENDDO
                0276        ENDDO
                0277       ENDDO
4578baf364 Jean*0278 
5a6ef5c2b4 Mich*0279 C     intN integrates to z=rC(1).  We want to integrate to z=0.
                0280 C     Assume that N(0)=0 and N(1)=0.
                0281 C     drC(1) is like half a grid cell.
4578baf364 Jean*0282       DO j=1-OLy,sNy+OLy
                0283        DO i=1-OLx,sNx+OLx
                0284         intN0(i,j) = intN(i,j,1)
5a6ef5c2b4 Mich*0285      &       + drC(1)*op5*BuoyFreq(i,j,2)
                0286        ENDDO
                0287       ENDDO
4578baf364 Jean*0288 
                0289       DO j=1-OLy,sNy+OLy
                0290        DO i=1-OLx,sNx+OLx
5a6ef5c2b4 Mich*0291         c1(i,j) = intN0(i,j)/pi
                0292         Rmid(i,j) = c1(i,j)/ABS(fCori(i,j,bi,bj))
                0293        ENDDO
4578baf364 Jean*0294       ENDDO
5a6ef5c2b4 Mich*0295 
4578baf364 Jean*0296       DO j=1-OLy,sNy+OLy
                0297        DO i=1-OLx,sNx+OLx
5a6ef5c2b4 Mich*0298         nEigs(i,j) = MIN(klow(i,j),nmodes)
                0299        ENDDO
                0300       ENDDO
                0301       DO k=1,Nr
4578baf364 Jean*0302        DO j=1-OLy,sNy+OLy
                0303         DO i=1-OLx,sNx+OLx
5a6ef5c2b4 Mich*0304          IF (mask(i,j,k).NE.0.0) THEN
                0305            DO m=1,nEigs(i,j)
                0306             vec(m,i,j,k) = -COS(intN(i,j,k)/(m*c1(i,j)))
                0307            ENDDO
                0308          ENDIF
                0309         ENDDO
                0310        ENDDO
                0311       ENDDO
                0312 
4578baf364 Jean*0313 C     The WKB approximation for the baroclinic mode does not always
5a6ef5c2b4 Mich*0314 C     integrate to zero so we adjust it.
4578baf364 Jean*0315       DO j=1-OLy,sNy+OLy
                0316        DO i=1-OLx,sNx+OLx
5a6ef5c2b4 Mich*0317         DO m=1,nEigs(i,j)
                0318          vecint(m,i,j) = zeroRL
                0319         ENDDO
                0320        ENDDO
                0321       ENDDO
                0322       DO k=1,Nr
4578baf364 Jean*0323        DO j=1-OLy,sNy+OLy
                0324         DO i=1-OLx,sNx+OLx
5a6ef5c2b4 Mich*0325          DO m=1,nEigs(i,j)
                0326          vecint(m,i,j) = vecint(m,i,j) + hfac(i,j,k)*vec(m,i,j,k)*drF(k)
                0327          ENDDO
                0328         ENDDO
                0329        ENDDO
                0330       ENDDO
4578baf364 Jean*0331       DO j=1-OLy,sNy+OLy
                0332        DO i=1-OLx,sNx+OLx
5a6ef5c2b4 Mich*0333         DO m=1,nEigs(i,j)
                0334          vecint(m,i,j) = vecint(m,i,j)/(-rlow(i,j)+small)
                0335         ENDDO
                0336        ENDDO
                0337       ENDDO
                0338       DO k=1,Nr
4578baf364 Jean*0339        DO j=1-OLy,sNy+OLy
                0340         DO i=1-OLx,sNx+OLx
5a6ef5c2b4 Mich*0341          DO m=1,nEigs(i,j)
                0342           vec(m,i,j,k) = vec(m,i,j,k) - vecint(m,i,j)
0d1e4b948d Mich*0343          ENDDO
                0344         ENDDO
                0345        ENDDO
                0346       ENDDO
5a6ef5c2b4 Mich*0347 # endif
0d1e4b948d Mich*0348 
5a6ef5c2b4 Mich*0349 C     Normalise the eigenvectors
4578baf364 Jean*0350       DO j=1-OLy,sNy+OLy
                0351        DO i=1-OLx,sNx+OLx
5a6ef5c2b4 Mich*0352         DO m=1,nmodes
                0353          vecint(m,i,j) = zeroRL
0d1e4b948d Mich*0354         ENDDO
                0355        ENDDO
                0356       ENDDO
                0357 
                0358       DO k=1,Nr
4578baf364 Jean*0359        DO j=1-OLy,sNy+OLy
                0360         DO i=1-OLx,sNx+OLx
5a6ef5c2b4 Mich*0361          DO m=1,nmodes
4578baf364 Jean*0362           vecint(m,i,j) = vecint(m,i,j) +
5a6ef5c2b4 Mich*0363      &         mask(i,j,k)*drF(k)*hfac(i,j,k)
                0364      &         *vec(m,i,j,k)*vec(m,i,j,k)
0d1e4b948d Mich*0365          ENDDO
                0366         ENDDO
                0367        ENDDO
                0368       ENDDO
                0369 
4578baf364 Jean*0370       DO j=1-OLy,sNy+OLy
                0371        DO i=1-OLx,sNx+OLx
5a6ef5c2b4 Mich*0372         DO m=1,nmodes
                0373          vecint(m,i,j) = SQRT(vecint(m,i,j))
                0374         ENDDO
                0375        ENDDO
                0376       ENDDO
0d1e4b948d Mich*0377 
                0378       DO k=1,Nr
4578baf364 Jean*0379        DO j=1-OLy,sNy+OLy
                0380         DO i=1-OLx,sNx+OLx
5a6ef5c2b4 Mich*0381          DO m=1,nmodes
                0382           vec(m,i,j,k) = vec(m,i,j,k)/(vecint(m,i,j)+small)
                0383          ENDDO
0d1e4b948d Mich*0384         ENDDO
                0385        ENDDO
                0386       ENDDO
                0387 
5a6ef5c2b4 Mich*0388 # ifdef ALLOW_DIAGNOSTICS
0d1e4b948d Mich*0389 C     Diagnostics
                0390       IF ( useDiagnostics.AND.writediag ) THEN
4578baf364 Jean*0391 #  ifdef HAVE_LAPACK
0d1e4b948d Mich*0392         CALL DIAGNOSTICS_FILL(a3d, 'GM_A3D  ',0,Nr,0,1,1,myThid)
                0393         CALL DIAGNOSTICS_FILL(b3d, 'GM_B3D  ',0,Nr,0,1,1,myThid)
                0394         CALL DIAGNOSTICS_FILL(c3d, 'GM_C3D  ',0,Nr,0,1,1,myThid)
5a6ef5c2b4 Mich*0395 #  endif
0d1e4b948d Mich*0396       ENDIF
5a6ef5c2b4 Mich*0397 # endif
f3b9070ca3 Jean*0398 
05118ac017 Jean*0399 #endif /* GM_BATES_K3D */
f3b9070ca3 Jean*0400 
                0401       RETURN
0d1e4b948d Mich*0402       END
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated from https://github.com/MITgcm/MITgcm by the 2.2.1-MITgcm-0.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!