/[MITgcm]/MITgcm/pkg/gmredi/gmredi_calc_eigs.F
ViewVC logotype

Diff of /MITgcm/pkg/gmredi/gmredi_calc_eigs.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.2 by jmc, Sat Jun 22 14:44:54 2013 UTC revision 1.6 by m_bates, Tue May 20 11:42:28 2014 UTC
# Line 9  C     !INTERFACE: Line 9  C     !INTERFACE:
9       I     iMin, iMax, jMin, jMax,       I     iMin, iMax, jMin, jMax,
10       I     bi, bj, N2, myThid, kLow,       I     bi, bj, N2, myThid, kLow,
11       I     mask, hfac, recip_hfac,       I     mask, hfac, recip_hfac,
12       I     writediag,       I     rlow, nmodes, writediag,
13       O     Kdef, vec)       O     Rmid, vec)
14    
15  C     !DESCRIPTION: \bv  C     !DESCRIPTION: \bv
16  C     *==========================================================*  C     *==========================================================*
# Line 36  C     bi, bj    :: tile indices Line 36  C     bi, bj    :: tile indices
36        INTEGER iMin,iMax,jMin,jMax        INTEGER iMin,iMax,jMin,jMax
37        INTEGER bi, bj        INTEGER bi, bj
38        INTEGER myThid        INTEGER myThid
39        INTEGER kLow(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        INTEGER nmodes
40        _RL mask(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        INTEGER kLow(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
41        _RL hfac(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL mask(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
42        _RL recip_hfac(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL hfac(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
43        _RL N2(  1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL recip_hfac(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
44        _RL Kdef(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rlow(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
45        _RL vec(GM_K3D_NModes,1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL N2(  1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
46          _RL Rmid(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
47          _RL vec(nmodes,1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
48    
49  #ifdef GM_K3D  #ifdef GM_K3D
50  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
51  C     == Local variables ==  C     == Local variables ==
52        INTEGER i,j,k,kk,m,info        INTEGER i,j,k,kk,m
53        _RL small  # ifdef HAVE_LAPACK
54        _RL a3d(   1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  C     info     :: error code from LAPACK
55        _RL b3d(   1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  C     idx      :: index used for sorting the eigenvalues
56        _RL c3d(   1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  C     a3d      :: lower diagonal of eigenvalue problem
57        _RL vecint(GM_K3D_NModes,1-OLx:sNx+OLx,1-OLy:sNy+OLy)  C     b3d      :: diagonal of eigenvalue problem
58        _RL val(   1-OLx:sNx+OLx,1-OLy:sNy+OLy)  C     c3d      :: upper diagonal of eigenvalue problem
59        _RL fCori2(1-OLx:sNx+OLx,1-OLy:sNy+OLy)  C     val      :: Eigenvalue (wavenumber) of the first baroclinic mode
60    C     eigR     :: Real component of all the eigenvalues in a water column
61    C     eigI     :: Imaginary component of all the eigenvalues in a water column
62    C     vecs     :: All the eigenvectors of a water column
63    C     dummy    :: Redundant variable for calling lapack
64    C     work     :: Work array for lapack
65    C     array    :: Array containing the matrix with a,b,c
66    C     eigval   :: Vector containing the eigenvalues
67          INTEGER info
68          INTEGER idx
69          _RL a3d(   1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
70          _RL b3d(   1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
71          _RL c3d(   1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
72          _RL val(   1-Olx:sNx+Olx,1-Oly:sNy+Oly)
73        _RL eigR(Nr),eigI(Nr),vecs(Nr,Nr),dummy(1,Nr),work(Nr*Nr)        _RL eigR(Nr),eigI(Nr),vecs(Nr,Nr),dummy(1,Nr),work(Nr*Nr)
74        _RL array(Nr,Nr)        _RL array(Nr,Nr)
75        _RL eigval(0:GM_K3D_NModes)        _RL eigval(0:nmodes)
76        INTEGER idx  # else
77  #ifdef ALLOW_DIAGNOSTICS  C     drNr     :: distance from bottom of cell to cell centre at Nr
78        _RL vec_diag(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  C     BuoyFreq :: buoyancy frequency, SQRT(N2)
79  #endif  C     intN     :: Vertical integral of BuoyFreq to each grid cell centre
80    C     intN0    :: Vertical integral of BuoyFreq to z=0
81    C     c1       :: intN0/pi
82    C     nEigs    :: number of eigenvalues/vectors to calculate
83          _RL drNr
84          _RL BuoyFreq(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr+1)
85          _RL intN(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
86          _RL intN0(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
87          _RL c1(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
88          INTEGER nEigs(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
89    # endif
90    C     small    :: a small number (used to avoid floating point exceptions)
91    C     vecint   :: vertical integral of eigenvector and/or square of eigenvector
92    C     fCori2   :: square of the Coriolis parameter
93          _RL small
94          _RL vecint(nmodes,1-Olx:sNx+Olx,1-Oly:sNy+Oly)
95          _RL fCori2(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
96          CHARACTER*(MAX_LEN_MBUF) msgBuf
97    
98        small = TINY(zeroRL)        small = TINY(zeroRL)
99    
100  C     Square of the Coriolis parameter  C     Square of the Coriolis parameter
# Line 72  C     Square of the Coriolis parameter Line 105  C     Square of the Coriolis parameter
105        ENDDO        ENDDO
106    
107        DO k=1,Nr        DO k=1,Nr
108         DO j=1-OLy,sNy+OLy         DO j=1-Oly,sNy+Oly
109          DO i=1-OLx,sNx+OLx          DO i=1-Olx,sNx+Olx
110           DO m=1,GM_K3D_NModes           DO m=1,nmodes
111            vec(m,i,j,k) = zeroRL            vec(m,i,j,k) = zeroRL
112           ENDDO           ENDDO
113          ENDDO          ENDDO
114         ENDDO         ENDDO
115        ENDDO        ENDDO
116    
117    # ifdef HAVE_LAPACK
118  C     Calculate the tridiagonal operator matrix for  C     Calculate the tridiagonal operator matrix for
119  C     f^2 d/dz 1/N^2 d/dz  C     f^2 d/dz 1/N^2 d/dz
120  C     a3d is the lower off-diagonal, b3d is the diagonal  C     a3d is the lower off-diagonal, b3d is the diagonal
# Line 96  C     and c3d is the upper off-diagonal Line 130  C     and c3d is the upper off-diagonal
130               b3d(i,j,k) = -c3d(i,j,k)               b3d(i,j,k) = -c3d(i,j,k)
131    
132             ELSEIF (k.LT.kLow(i,j)) THEN             ELSEIF (k.LT.kLow(i,j)) THEN
              IF (k+1.GT.Nr) THEN  
                PRINT '(a4,x,3(a1,i2),a1,3(x,i2))', 'kp1:',  
      &              '(',i,',',j,',',k,')',  
      &              kLow(i,j), k+1, Nr  
              ENDIF  
133               a3d(i,j,k) = fCori2(i,j)*recip_hFac(i,j,k)               a3d(i,j,k) = fCori2(i,j)*recip_hFac(i,j,k)
134       &                    *recip_drF(k)*recip_drC(k)/N2(i,j,k)       &                    *recip_drF(k)*recip_drC(k)/N2(i,j,k)
135               c3d(i,j,k) = fCori2(i,j)*recip_hFac(i,j,k)               c3d(i,j,k) = fCori2(i,j)*recip_hFac(i,j,k)
# Line 108  C     and c3d is the upper off-diagonal Line 137  C     and c3d is the upper off-diagonal
137               b3d(i,j,k) = -a3d(i,j,k)-c3d(i,j,k)               b3d(i,j,k) = -a3d(i,j,k)-c3d(i,j,k)
138    
139             ELSEIF (k.EQ.kLow(i,j)) THEN             ELSEIF (k.EQ.kLow(i,j)) THEN
              IF (k.GT.Nr) THEN  
                PRINT '(a4,x,3(a1,i2),a1,3(x,i2))', 'k  :',  
      &              '(',i,',',j,',',k,')',  
      &              kLow(i,j), k, Nr  
              ENDIF  
140               a3d(i,j,k) = fCori2(i,j)*recip_hFac(i,j,k)               a3d(i,j,k) = fCori2(i,j)*recip_hFac(i,j,k)
141       &                    *recip_drF(k)*recip_drC(k)/N2(i,j,k)       &                    *recip_drF(k)*recip_drC(k)/N2(i,j,k)
142               c3d(i,j,k) = zeroRL               c3d(i,j,k) = zeroRL
# Line 133  C     and c3d is the upper off-diagonal Line 157  C     and c3d is the upper off-diagonal
157         ENDDO         ENDDO
158        ENDDO        ENDDO
159    
160  #ifdef use_lapack        DO j=1-Oly,sNy+Oly
161           DO i=1-Olx,sNx+Olx
       DO j=1-OLy,sNy+OLy  
        DO i=1-OLx,sNx+OLx  
162          IF (kLow(i,j).GT.0) THEN          IF (kLow(i,j).GT.0) THEN
163            DO kk=1,Nr            DO kk=1,Nr
164             DO k=1,Nr             DO k=1,Nr
# Line 158  C     and c3d is the upper off-diagonal Line 180  C     and c3d is the upper off-diagonal
180    
181            CALL DGEEV('N','V',Nr,array,Nr,eigR,eigI,dummy,1,            CALL DGEEV('N','V',Nr,array,Nr,eigR,eigI,dummy,1,
182       &         vecs,Nr,work,Nr*Nr,info)       &         vecs,Nr,work,Nr*Nr,info)
183  C         Find the second largest eigenvalue (the Rossby radius)            IF( info.LT.0 ) THEN
184                WRITE(msgBuf,'(A,x,2(A1,I2),A1,x,A,I4)')
185         &           'GMREDI_CALC_EIGS problem with arguments for DGEEV at',
186         &           '(',i,',',j,')', 'error code =',info
187                CALL PRINT_ERROR( msgBuf , myThid )
188                
189              ELSEIF(info.GT.0 ) THEN
190                WRITE(msgBuf,'(A,x,2(A1,I2),A1,x,A,I4)')
191         &           'GMREDI_CALC_EIGS problems with eigensolver DGEEV at',
192         &           '(',i,',',j,')', 'error code =',info
193                CALL PRINT_ERROR( msgBuf , myThid )
194                
195              ENDIF
196    
197    C         Find the second largest eigenvalue (the Rossby radius)
198  C         and the first M baroclinic modes (eigenvectors)  C         and the first M baroclinic modes (eigenvectors)
199            DO m=0,GM_K3D_NModes            DO m=0,nmodes
200             eigval(m) = -HUGE(zeroRL)             eigval(m) = -HUGE(zeroRL)
201            ENDDO            ENDDO
202    
203            DO k=1,kLow(i,j)            DO k=1,kLow(i,j)
204             eigval(0) = MAX(eigval(0),eigR(k))             eigval(0) = MAX(eigval(0),eigR(k))
205            ENDDO            ENDDO
206            DO m=1,MIN(GM_K3D_NModes,klow(i,j))            DO m=1,MIN(nmodes,klow(i,j)-1)
 C          DO m=1,GM_K3D_NModes  
207             DO k=1,kLow(i,j)             DO k=1,kLow(i,j)
208              IF (eigR(k).LT.eigval(m-1)) THEN              IF (eigR(k).LT.eigval(m-1)) THEN
209                eigval(m) = MAX(eigval(m),eigR(k))                eigval(m) = MAX(eigval(m),eigR(k))
# Line 189  C          DO m=1,GM_K3D_NModes Line 224  C          DO m=1,GM_K3D_NModes
224          ELSE          ELSE
225            val(i,j)=zeroRL            val(i,j)=zeroRL
226            DO k=1,Nr            DO k=1,Nr
227             DO m=1,GM_K3D_NModes             DO m=1,nmodes
228              vec(m,i,j,k)=zeroRL              vec(m,i,j,k)=zeroRL
229             ENDDO             ENDDO
230            ENDDO            ENDDO
# Line 198  C          DO m=1,GM_K3D_NModes Line 233  C          DO m=1,GM_K3D_NModes
233         ENDDO         ENDDO
234        ENDDO        ENDDO
235    
236  C     Normalise the eigenvectors        DO k=1,Nr
237        DO j=1-OLy,sNy+OLy         DO j=1-Oly,sNy+Oly
238         DO i=1-OLx,sNx+OLx          DO i=1-Olx,sNx+Olx
239          DO m=1,GM_K3D_NModes           IF (kLow(i,j).GT.2 .AND. val(i,j).NE.zeroRL) THEN
240           vecint(m,i,j) = zeroRL             Rmid(i,j) = 1.0/(SQRT(ABS(val(i,j)))+small)
241             ELSE
242               Rmid(i,j) = zeroRL
243             ENDIF
244          ENDDO          ENDDO
245         ENDDO         ENDDO
246        ENDDO        ENDDO
247    
248    # else      
249        DO k=1,Nr        DO k=1,Nr
250         DO j=1-OLy,sNy+OLy         DO j=1-Oly,sNy+Oly
251          DO i=1-OLx,sNx+OLx          DO i=1-Olx,sNx+Olx
252           DO m=1,GM_K3D_NModes           BuoyFreq(i,j,k) = mask(i,j,k)*SQRT(N2(i,j,k))
253            vecint(m,i,j) = vecint(m,i,j) +          ENDDO
254       &         mask(i,j,k)*drF(k)*hfac(i,j,k)         ENDDO
255       &         *vec(m,i,j,k)*vec(m,i,j,k)        ENDDO
256          k=Nr+1
257          DO j=1-Oly,sNy+Oly
258           DO i=1-Olx,sNx+Olx
259            BuoyFreq(i,j,k) = zeroRL
260           ENDDO
261          ENDDO
262    
263    C     integrate N using something like Simpson s rule (but not quite)
264    C     drC*( (N(k+1)+N(k+2))/2 + (N(k)+N(k+1))/2 )/2
265    C     when k=Nr, say that N(k+2)=0 and N(k+1)=0
266          k=Nr
267    C     drNr is like drC(Nr+1)/2
268          drNr = rC(Nr)-rF(Nr+1)
269          DO j=1-Oly,sNy+Oly
270           DO i=1-Olx,sNx+Olx
271            intN(i,j,k) = op5*BuoyFreq(i,j,k)*drNr
272           ENDDO
273          ENDDO
274          DO k=Nr-1,1,-1
275           DO j=1-Oly,sNy+Oly
276            DO i=1-Olx,sNx+Olx
277             intN(i,j,k) = intN(i,j,k+1)
278         &        + drC(k)*( op25*BuoyFreq(i,j,k+2) + op5*BuoyFreq(i,j,k)
279         &                 + op25*BuoyFreq(i,j,k+1) )
280            ENDDO
281           ENDDO
282          ENDDO
283          
284    C     intN integrates to z=rC(1).  We want to integrate to z=0.
285    C     Assume that N(0)=0 and N(1)=0.
286    C     drC(1) is like half a grid cell.
287          DO j=1-Oly,sNy+Oly
288           DO i=1-Olx,sNx+Olx
289            intN0(i,j) = intN(i,j,1)
290         &       + drC(1)*op5*BuoyFreq(i,j,2)
291           ENDDO
292          ENDDO
293          
294          DO j=1-Oly,sNy+Oly
295           DO i=1-Olx,sNx+Olx
296            c1(i,j) = intN0(i,j)/pi
297            Rmid(i,j) = c1(i,j)/ABS(fCori(i,j,bi,bj))
298           ENDDO
299          ENDDO      
300    
301          DO j=1-Oly,sNy+Oly
302           DO i=1-Olx,sNx+Olx
303            nEigs(i,j) = MIN(klow(i,j),nmodes)
304           ENDDO
305          ENDDO
306          DO k=1,Nr
307           DO j=1-Oly,sNy+Oly
308            DO i=1-Olx,sNx+Olx
309             IF (mask(i,j,k).NE.0.0) THEN
310               DO m=1,nEigs(i,j)
311                vec(m,i,j,k) = -COS(intN(i,j,k)/(m*c1(i,j)))
312               ENDDO
313             ENDIF
314            ENDDO
315           ENDDO
316          ENDDO
317    
318    C     The WKB approximation for the baroclinic mode does not always
319    C     integrate to zero so we adjust it.
320          DO j=1-Oly,sNy+Oly
321           DO i=1-Olx,sNx+Olx
322            DO m=1,nEigs(i,j)
323             vecint(m,i,j) = zeroRL
324            ENDDO
325           ENDDO
326          ENDDO
327          DO k=1,Nr
328           DO j=1-Oly,sNy+Oly
329            DO i=1-Olx,sNx+Olx
330             DO m=1,nEigs(i,j)
331             vecint(m,i,j) = vecint(m,i,j) + hfac(i,j,k)*vec(m,i,j,k)*drF(k)
332             ENDDO
333            ENDDO
334           ENDDO
335          ENDDO
336          DO j=1-Oly,sNy+Oly
337           DO i=1-Olx,sNx+Olx
338            DO m=1,nEigs(i,j)
339             vecint(m,i,j) = vecint(m,i,j)/(-rlow(i,j)+small)
340            ENDDO
341           ENDDO
342          ENDDO
343          DO k=1,Nr
344           DO j=1-Oly,sNy+Oly
345            DO i=1-Olx,sNx+Olx
346             DO m=1,nEigs(i,j)
347              vec(m,i,j,k) = vec(m,i,j,k) - vecint(m,i,j)
348           ENDDO           ENDDO
349          ENDDO          ENDDO
350         ENDDO         ENDDO
351        ENDDO        ENDDO
352    # endif
353    
354        DO j=1-OLy,sNy+OLy  C     Normalise the eigenvectors
355         DO i=1-OLx,sNx+OLx        DO j=1-Oly,sNy+Oly
356          DO m=1,GM_K3D_NModes         DO i=1-Olx,sNx+Olx
357           vecint(m,i,j) = SQRT(vecint(m,i,j))          DO m=1,nmodes
358             vecint(m,i,j) = zeroRL
359          ENDDO          ENDDO
360         ENDDO         ENDDO
361        ENDDO        ENDDO
362    
363        DO k=1,Nr        DO k=1,Nr
364         DO j=1-OLy,sNy+OLy         DO j=1-Oly,sNy+Oly
365          DO i=1-OLx,sNx+OLx          DO i=1-Olx,sNx+Olx
366           DO m=1,GM_K3D_NModes           DO m=1,nmodes
367            vec(m,i,j,k) = vec(m,i,j,k)/(vecint(m,i,j)+small)            vecint(m,i,j) = vecint(m,i,j) +
368         &         mask(i,j,k)*drF(k)*hfac(i,j,k)
369         &         *vec(m,i,j,k)*vec(m,i,j,k)
370           ENDDO           ENDDO
371          ENDDO          ENDDO
372         ENDDO         ENDDO
373        ENDDO        ENDDO
374    
375  #else        DO j=1-Oly,sNy+Oly
376  C      CALL EIGENVAL(a3d,b3d,c3d,bi,bj,val)         DO i=1-Olx,sNx+Olx
377  C      CALL EIGENVEC(a3d,b3d,c3d,val,1,bi,bj,vec)          DO m=1,nmodes
378  #endif           vecint(m,i,j) = SQRT(vecint(m,i,j))
379            ENDDO
380           ENDDO
381          ENDDO
382    
383        DO k=1,Nr        DO k=1,Nr
384         DO j=1-OLy,sNy+OLy         DO j=1-Oly,sNy+Oly
385          DO i=1-OLx,sNx+OLx          DO i=1-Olx,sNx+Olx
386           IF (kLow(i,j).GT.2 .AND. val(i,j).NE.zeroRL) THEN           DO m=1,nmodes
387             Kdef(i,j) = SQRT(ABS(val(i,j)))            vec(m,i,j,k) = vec(m,i,j,k)/(vecint(m,i,j)+small)
388           ELSE           ENDDO
            Kdef(i,j) = zeroRL  
          ENDIF  
389          ENDDO          ENDDO
390         ENDDO         ENDDO
391        ENDDO        ENDDO
392    
393  #ifdef ALLOW_DIAGNOSTICS  # ifdef ALLOW_DIAGNOSTICS
394  C     Diagnostics  C     Diagnostics
395        IF ( useDiagnostics.AND.writediag ) THEN        IF ( useDiagnostics.AND.writediag ) THEN
396    #  ifdef HAVE_LAPACK        
397          CALL DIAGNOSTICS_FILL(a3d, 'GM_A3D  ',0,Nr,0,1,1,myThid)          CALL DIAGNOSTICS_FILL(a3d, 'GM_A3D  ',0,Nr,0,1,1,myThid)
398          CALL DIAGNOSTICS_FILL(b3d, 'GM_B3D  ',0,Nr,0,1,1,myThid)          CALL DIAGNOSTICS_FILL(b3d, 'GM_B3D  ',0,Nr,0,1,1,myThid)
399          CALL DIAGNOSTICS_FILL(c3d, 'GM_C3D  ',0,Nr,0,1,1,myThid)          CALL DIAGNOSTICS_FILL(c3d, 'GM_C3D  ',0,Nr,0,1,1,myThid)
400          CALL DIAGNOSTICS_FILL(val, 'GM_VAL  ',0, 1,0,1,1,myThid)  #  endif
         DO k=1,Nr  
          DO j=1-OLy,sNy+OLy  
           DO i=1-OLx,sNx+OLx  
            vec_diag(i,j,k) = vec(1,i,j,k)  
           ENDDO  
          ENDDO  
         ENDDO  
   
         CALL DIAGNOSTICS_FILL(vec_diag, 'GM_VEC  ',0,Nr,0,1,1,myThid)  
401        ENDIF        ENDIF
402  #endif  # endif
403    
404  #endif /* GM_K3D */  #endif /* GM_K3D */
405    
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.6
  ViewVC Help
Powered by ViewVC 1.1.22