/[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.1 by m_bates, Fri Jun 21 17:23:31 2013 UTC revision 1.2 by jmc, Sat Jun 22 14:44:54 2013 UTC
# Line 2  C $Header$ Line 2  C $Header$
2  C $Name$  C $Name$
3    
4  #include "GMREDI_OPTIONS.h"  #include "GMREDI_OPTIONS.h"
5    
6  C     !ROUTINE: EIGENVAL  C     !ROUTINE: EIGENVAL
7  C     !INTERFACE:  C     !INTERFACE:
8        SUBROUTINE GMREDI_CALC_EIGS(        SUBROUTINE GMREDI_CALC_EIGS(
# Line 14  C     !INTERFACE: Line 15  C     !INTERFACE:
15  C     !DESCRIPTION: \bv  C     !DESCRIPTION: \bv
16  C     *==========================================================*  C     *==========================================================*
17  C     | SUBROUTINE GMREDI_CALC_URMS  C     | SUBROUTINE GMREDI_CALC_URMS
18  C     | o Calculate the vertical structure of the rms eddy  C     | o Calculate the vertical structure of the rms eddy
19  C     |   velocity based on baroclinic modal decomposition  C     |   velocity based on baroclinic modal decomposition
20  C     *==========================================================*  C     *==========================================================*
21  C     \ev  C     \ev
# Line 35  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 kLow(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
40        _RL mask(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL mask(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
41        _RL hfac(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL hfac(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
42        _RL recip_hfac(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL recip_hfac(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
43        _RL N2(  1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL N2(  1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
44        _RL Kdef(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL Kdef(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
45        _RL vec(GM_K3D_NModes,1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL vec(GM_K3D_NModes,1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
46          
47    #ifdef GM_K3D
48  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
49  C     == Local variables ==  C     == Local variables ==
50        INTEGER i,j,k,kk,m,info        INTEGER i,j,k,kk,m,info
51        _RL small        _RL small
52        _RL a3d(   1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL a3d(   1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
53        _RL b3d(   1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL b3d(   1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
54        _RL c3d(   1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL c3d(   1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
55        _RL vecint(GM_K3D_NModes,1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL vecint(GM_K3D_NModes,1-OLx:sNx+OLx,1-OLy:sNy+OLy)
56        _RL val(   1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL val(   1-OLx:sNx+OLx,1-OLy:sNy+OLy)
57        _RL fCori2(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL fCori2(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
58        _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)
59        _RL array(Nr,Nr)        _RL array(Nr,Nr)
60        _RL eigval(0:GM_K3D_NModes)        _RL eigval(0:GM_K3D_NModes)
61        INTEGER idx        INTEGER idx
62  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
63        _RL vec_diag(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL vec_diag(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
64  #endif  #endif
65        small = TINY(zeroRL)        small = TINY(zeroRL)
66    
67  C     Square of the Coriolis parameter  C     Square of the Coriolis parameter
68        DO i=1-Olx,sNx+Olx        DO i=1-OLx,sNx+OLx
69         DO j=1-Oly,sNy+Oly         DO j=1-OLy,sNy+OLy
70          fCori2(i,j) = fCori(i,j,bi,bj)*fCori(i,j,bi,bj)          fCori2(i,j) = fCori(i,j,bi,bj)*fCori(i,j,bi,bj)
71         ENDDO         ENDDO
72        ENDDO        ENDDO
73    
74        DO k=1,Nr        DO k=1,Nr
75         DO j=1-Oly,sNy+Oly         DO j=1-OLy,sNy+OLy
76          DO i=1-Olx,sNx+Olx          DO i=1-OLx,sNx+OLx
77           DO m=1,GM_K3D_NModes           DO m=1,GM_K3D_NModes
78            vec(m,i,j,k) = zeroRL            vec(m,i,j,k) = zeroRL
79           ENDDO           ENDDO
# Line 84  C     f^2 d/dz 1/N^2 d/dz Line 86  C     f^2 d/dz 1/N^2 d/dz
86  C     a3d is the lower off-diagonal, b3d is the diagonal  C     a3d is the lower off-diagonal, b3d is the diagonal
87  C     and c3d is the upper off-diagonal  C     and c3d is the upper off-diagonal
88        DO k=1,Nr        DO k=1,Nr
89         DO j=1-Oly,sNy+Oly         DO j=1-OLy,sNy+OLy
90          DO i=1-Olx,sNx+Olx          DO i=1-OLx,sNx+OLx
91           IF (kLow(i,j) .GT. 0) THEN           IF (kLow(i,j) .GT. 0) THEN
92             IF (k.EQ.1) THEN             IF (k.EQ.1) THEN
93               a3d(i,j,k) = zeroRL               a3d(i,j,k) = zeroRL
# Line 95  C     and c3d is the upper off-diagonal Line 97  C     and c3d is the upper off-diagonal
97    
98             ELSEIF (k.LT.kLow(i,j)) THEN             ELSEIF (k.LT.kLow(i,j)) THEN
99               IF (k+1.GT.Nr) THEN               IF (k+1.GT.Nr) THEN
100                 PRINT '(a4,x,3(a1,i2),a1,3(x,i2))', 'kp1:',                 PRINT '(a4,x,3(a1,i2),a1,3(x,i2))', 'kp1:',
101       &              '(',i,',',j,',',k,')',       &              '(',i,',',j,',',k,')',
102       &              kLow(i,j), k+1, Nr       &              kLow(i,j), k+1, Nr
103               ENDIF               ENDIF
# Line 107  C     and c3d is the upper off-diagonal Line 109  C     and c3d is the upper off-diagonal
109    
110             ELSEIF (k.EQ.kLow(i,j)) THEN             ELSEIF (k.EQ.kLow(i,j)) THEN
111               IF (k.GT.Nr) THEN               IF (k.GT.Nr) THEN
112                 PRINT '(a4,x,3(a1,i2),a1,3(x,i2))', 'k  :',                 PRINT '(a4,x,3(a1,i2),a1,3(x,i2))', 'k  :',
113       &              '(',i,',',j,',',k,')',       &              '(',i,',',j,',',k,')',
114       &              kLow(i,j), k, Nr       &              kLow(i,j), k, Nr
115               ENDIF               ENDIF
# Line 132  C     and c3d is the upper off-diagonal Line 134  C     and c3d is the upper off-diagonal
134        ENDDO        ENDDO
135    
136  #ifdef use_lapack  #ifdef use_lapack
137          
138        DO j=1-Oly,sNy+Oly        DO j=1-OLy,sNy+OLy
139         DO i=1-Olx,sNx+Olx         DO i=1-OLx,sNx+OLx
140          IF (kLow(i,j).GT.0) THEN          IF (kLow(i,j).GT.0) THEN
141            DO kk=1,Nr            DO kk=1,Nr
142             DO k=1,Nr             DO k=1,Nr
143              array(k,kk) = zeroRL              array(k,kk) = zeroRL
144             ENDDO             ENDDO
145            ENDDO            ENDDO
146              
147            k=1            k=1
148            array(k,k)   = b3d(i,j,k)            array(k,k)   = b3d(i,j,k)
149            array(k,k+1) = c3d(i,j,k)            array(k,k+1) = c3d(i,j,k)
# Line 153  C     and c3d is the upper off-diagonal Line 155  C     and c3d is the upper off-diagonal
155            k=Nr            k=Nr
156            array(k,k-1) = a3d(i,j,k)            array(k,k-1) = a3d(i,j,k)
157            array(k,k)   = b3d(i,j,k)            array(k,k)   = b3d(i,j,k)
158              
159            CALL DGEEV('N','V',Nr,array,Nr,eigR,eigI,dummy,1,            CALL DGEEV('N','V',Nr,array,Nr,eigR,eigI,dummy,1,
160       &         vecs,Nr,work,Nr*Nr,info)       &         vecs,Nr,work,Nr*Nr,info)
161  C         Find the second largest eigenvalue (the Rossby radius)  C         Find the second largest eigenvalue (the Rossby radius)
162  C         and the first M baroclinic modes (eigenvectors)  C         and the first M baroclinic modes (eigenvectors)
163            DO m=0,GM_K3D_NModes            DO m=0,GM_K3D_NModes
164             eigval(m) = -HUGE(zeroRL)             eigval(m) = -HUGE(zeroRL)
165            ENDDO            ENDDO
166              
167            DO k=1,kLow(i,j)            DO k=1,kLow(i,j)
168             eigval(0) = MAX(eigval(0),eigR(k))             eigval(0) = MAX(eigval(0),eigR(k))
169            ENDDO            ENDDO
# Line 191  C          DO m=1,GM_K3D_NModes Line 193  C          DO m=1,GM_K3D_NModes
193              vec(m,i,j,k)=zeroRL              vec(m,i,j,k)=zeroRL
194             ENDDO             ENDDO
195            ENDDO            ENDDO
196              
197          ENDIF          ENDIF
198         ENDDO         ENDDO
199        ENDDO        ENDDO
200    
201  C     Normalise the eigenvectors  C     Normalise the eigenvectors
202        DO j=1-Oly,sNy+Oly        DO j=1-OLy,sNy+OLy
203         DO i=1-Olx,sNx+Olx         DO i=1-OLx,sNx+OLx
204          DO m=1,GM_K3D_NModes          DO m=1,GM_K3D_NModes
205           vecint(m,i,j) = zeroRL           vecint(m,i,j) = zeroRL
206          ENDDO          ENDDO
# Line 206  C     Normalise the eigenvectors Line 208  C     Normalise the eigenvectors
208        ENDDO        ENDDO
209    
210        DO k=1,Nr        DO k=1,Nr
211         DO j=1-Oly,sNy+Oly         DO j=1-OLy,sNy+OLy
212          DO i=1-Olx,sNx+Olx          DO i=1-OLx,sNx+OLx
213           DO m=1,GM_K3D_NModes           DO m=1,GM_K3D_NModes
214            vecint(m,i,j) = vecint(m,i,j) +            vecint(m,i,j) = vecint(m,i,j) +
215       &         mask(i,j,k)*drF(k)*hfac(i,j,k)       &         mask(i,j,k)*drF(k)*hfac(i,j,k)
216       &         *vec(m,i,j,k)*vec(m,i,j,k)       &         *vec(m,i,j,k)*vec(m,i,j,k)
217           ENDDO           ENDDO
# Line 217  C     Normalise the eigenvectors Line 219  C     Normalise the eigenvectors
219         ENDDO         ENDDO
220        ENDDO        ENDDO
221    
222        DO j=1-Oly,sNy+Oly        DO j=1-OLy,sNy+OLy
223         DO i=1-Olx,sNx+Olx         DO i=1-OLx,sNx+OLx
224          DO m=1,GM_K3D_NModes          DO m=1,GM_K3D_NModes
225           vecint(m,i,j) = SQRT(vecint(m,i,j))           vecint(m,i,j) = SQRT(vecint(m,i,j))
226          ENDDO          ENDDO
# Line 226  C     Normalise the eigenvectors Line 228  C     Normalise the eigenvectors
228        ENDDO        ENDDO
229    
230        DO k=1,Nr        DO k=1,Nr
231         DO j=1-Oly,sNy+Oly         DO j=1-OLy,sNy+OLy
232          DO i=1-Olx,sNx+Olx          DO i=1-OLx,sNx+OLx
233           DO m=1,GM_K3D_NModes           DO m=1,GM_K3D_NModes
234            vec(m,i,j,k) = vec(m,i,j,k)/(vecint(m,i,j)+small)            vec(m,i,j,k) = vec(m,i,j,k)/(vecint(m,i,j)+small)
235           ENDDO           ENDDO
# Line 235  C     Normalise the eigenvectors Line 237  C     Normalise the eigenvectors
237         ENDDO         ENDDO
238        ENDDO        ENDDO
239    
240  #else        #else
241  C      CALL EIGENVAL(a3d,b3d,c3d,bi,bj,val)  C      CALL EIGENVAL(a3d,b3d,c3d,bi,bj,val)
242  C      CALL EIGENVEC(a3d,b3d,c3d,val,1,bi,bj,vec)  C      CALL EIGENVEC(a3d,b3d,c3d,val,1,bi,bj,vec)
243  #endif  #endif
244    
245        DO k=1,Nr        DO k=1,Nr
246         DO j=1-Oly,sNy+Oly         DO j=1-OLy,sNy+OLy
247          DO i=1-Olx,sNx+Olx          DO i=1-OLx,sNx+OLx
248           IF (kLow(i,j).GT.2 .AND. val(i,j).NE.zeroRL) THEN           IF (kLow(i,j).GT.2 .AND. val(i,j).NE.zeroRL) THEN
249             Kdef(i,j) = SQRT(ABS(val(i,j)))             Kdef(i,j) = SQRT(ABS(val(i,j)))
250           ELSE           ELSE
# Line 260  C     Diagnostics Line 262  C     Diagnostics
262          CALL DIAGNOSTICS_FILL(c3d, 'GM_C3D  ',0,Nr,0,1,1,myThid)          CALL DIAGNOSTICS_FILL(c3d, 'GM_C3D  ',0,Nr,0,1,1,myThid)
263          CALL DIAGNOSTICS_FILL(val, 'GM_VAL  ',0, 1,0,1,1,myThid)          CALL DIAGNOSTICS_FILL(val, 'GM_VAL  ',0, 1,0,1,1,myThid)
264          DO k=1,Nr          DO k=1,Nr
265           DO j=1-Oly,sNy+Oly           DO j=1-OLy,sNy+OLy
266            DO i=1-Olx,sNx+Olx            DO i=1-OLx,sNx+OLx
267             vec_diag(i,j,k) = vec(1,i,j,k)             vec_diag(i,j,k) = vec(1,i,j,k)
268            ENDDO            ENDDO
269           ENDDO           ENDDO
# Line 270  C     Diagnostics Line 272  C     Diagnostics
272          CALL DIAGNOSTICS_FILL(vec_diag, 'GM_VEC  ',0,Nr,0,1,1,myThid)          CALL DIAGNOSTICS_FILL(vec_diag, 'GM_VEC  ',0,Nr,0,1,1,myThid)
273        ENDIF        ENDIF
274  #endif  #endif
275    
276    #endif /* GM_K3D */
277    
278          RETURN
279        END        END
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.2
  ViewVC Help
Powered by ViewVC 1.1.22