pkg/gmredi/gmredi_calc_eigs.F

C $Header: /u/gcmpack/MITgcm/pkg/gmredi/gmredi_calc_eigs.F,v 1.1 2013/06/21 17:23:31 m_bates 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     writediag,
     O     Kdef, 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 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 N2(  1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
      _RL Kdef(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL vec(GM_K3D_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,info
      _RL small
      _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 vecint(GM_K3D_NModes,1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL val(   1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL fCori2(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:GM_K3D_NModes)
      INTEGER idx
#ifdef ALLOW_DIAGNOSTICS
      _RL vec_diag(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
#endif
      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,GM_K3D_NModes
          vec(m,i,j,k) = zeroRL
         ENDDO
        ENDDO
       ENDDO
      ENDDO

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
             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
             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
             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
             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

#ifdef use_lapack

      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)
C         Find the second largest eigenvalue (the Rossby radius)
C         and the first M baroclinic modes (eigenvectors)
          DO m=0,GM_K3D_NModes
           eigval(m) = -HUGE(zeroRL)
          ENDDO

          DO k=1,kLow(i,j)
           eigval(0) = MAX(eigval(0),eigR(k))
          ENDDO
          DO m=1,MIN(GM_K3D_NModes,klow(i,j))
C          DO m=1,GM_K3D_NModes
           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,GM_K3D_NModes
            vec(m,i,j,k)=zeroRL
           ENDDO
          ENDDO

        ENDIF
       ENDDO
      ENDDO

C     Normalise the eigenvectors
      DO j=1-OLy,sNy+OLy
       DO i=1-OLx,sNx+OLx
        DO m=1,GM_K3D_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,GM_K3D_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,GM_K3D_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,GM_K3D_NModes
          vec(m,i,j,k) = vec(m,i,j,k)/(vecint(m,i,j)+small)
         ENDDO
        ENDDO
       ENDDO
      ENDDO

#else
C      CALL EIGENVAL(a3d,b3d,c3d,bi,bj,val)
C      CALL EIGENVEC(a3d,b3d,c3d,val,1,bi,bj,vec)
#endif

      DO k=1,Nr
       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
           Kdef(i,j) = SQRT(ABS(val(i,j)))
         ELSE
           Kdef(i,j) = zeroRL
         ENDIF
        ENDDO
       ENDDO
      ENDDO

#ifdef ALLOW_DIAGNOSTICS
C     Diagnostics
      IF ( useDiagnostics.AND.writediag ) THEN
        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)
        CALL DIAGNOSTICS_FILL(val, 'GM_VAL  ',0, 1,0,1,1,myThid)
        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)
      ENDIF
#endif

#endif /* GM_K3D */

      RETURN
      END
1	jmc	1.2	C $Header: /u/gcmpack/MITgcm/pkg/gmredi/gmredi_calc_eigs.F,v 1.1 2013/06/21 17:23:31 m_bates Exp $
2	m_bates	1.1	C $Name: $
3
4			#include "GMREDI_OPTIONS.h"
5	jmc	1.2
6	m_bates	1.1	C !ROUTINE: EIGENVAL
7			C !INTERFACE:
8			SUBROUTINE GMREDI_CALC_EIGS(
9			I iMin, iMax, jMin, jMax,
10			I bi, bj, N2, myThid, kLow,
11			I mask, hfac, recip_hfac,
12			I writediag,
13			O Kdef, vec)
14
15			C !DESCRIPTION: \bv
16			C ==========================================================
17			C \| SUBROUTINE GMREDI_CALC_URMS
18	jmc	1.2	C \| o Calculate the vertical structure of the rms eddy
19	m_bates	1.1	C \| velocity based on baroclinic modal decomposition
20			C ==========================================================
21			C \ev
22
23			IMPLICIT NONE
24
25			C == Global variables ==
26			#include "SIZE.h"
27			#include "GRID.h"
28			#include "EEPARAMS.h"
29			#include "PARAMS.h"
30			#include "GMREDI.h"
31
32			C !INPUT/OUTPUT PARAMETERS:
33			C == Routine arguments ==
34			C bi, bj :: tile indices
35			LOGICAL writediag
36			INTEGER iMin,iMax,jMin,jMax
37			INTEGER bi, bj
38			INTEGER myThid
39	jmc	1.2	INTEGER kLow(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
40			_RL mask(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
41			_RL hfac(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
42			_RL recip_hfac(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
43			_RL N2( 1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
44			_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)
46
47			#ifdef GM_K3D
48	m_bates	1.1	C !LOCAL VARIABLES:
49			C == Local variables ==
50			INTEGER i,j,k,kk,m,info
51			_RL small
52	jmc	1.2	_RL a3d( 1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
53			_RL b3d( 1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
54			_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)
56			_RL val( 1-OLx:sNx+OLx,1-OLy:sNy+OLy)
57			_RL fCori2(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
58	m_bates	1.1	_RL eigR(Nr),eigI(Nr),vecs(Nr,Nr),dummy(1,Nr),work(Nr*Nr)
59			_RL array(Nr,Nr)
60			_RL eigval(0:GM_K3D_NModes)
61			INTEGER idx
62			#ifdef ALLOW_DIAGNOSTICS
63	jmc	1.2	_RL vec_diag(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
64	m_bates	1.1	#endif
65			small = TINY(zeroRL)
66
67			C Square of the Coriolis parameter
68	jmc	1.2	DO i=1-OLx,sNx+OLx
69			DO j=1-OLy,sNy+OLy
70	m_bates	1.1	fCori2(i,j) = fCori(i,j,bi,bj)*fCori(i,j,bi,bj)
71			ENDDO
72			ENDDO
73
74			DO k=1,Nr
75	jmc	1.2	DO j=1-OLy,sNy+OLy
76			DO i=1-OLx,sNx+OLx
77	m_bates	1.1	DO m=1,GM_K3D_NModes
78			vec(m,i,j,k) = zeroRL
79			ENDDO
80			ENDDO
81			ENDDO
82			ENDDO
83
84			C Calculate the tridiagonal operator matrix for
85			C f^2 d/dz 1/N^2 d/dz
86			C a3d is the lower off-diagonal, b3d is the diagonal
87			C and c3d is the upper off-diagonal
88			DO k=1,Nr
89	jmc	1.2	DO j=1-OLy,sNy+OLy
90			DO i=1-OLx,sNx+OLx
91	m_bates	1.1	IF (kLow(i,j) .GT. 0) THEN
92			IF (k.EQ.1) THEN
93			a3d(i,j,k) = zeroRL
94			c3d(i,j,k) = fCori2(i,j)*recip_hFac(i,j,k)
95			& recip_drC(k+1)recip_drF(k)/N2(i,j,k+1)
96			b3d(i,j,k) = -c3d(i,j,k)
97
98			ELSEIF (k.LT.kLow(i,j)) THEN
99			IF (k+1.GT.Nr) THEN
100	jmc	1.2	PRINT '(a4,x,3(a1,i2),a1,3(x,i2))', 'kp1:',
101	m_bates	1.1	& '(',i,',',j,',',k,')',
102			& kLow(i,j), k+1, Nr
103			ENDIF
104			a3d(i,j,k) = fCori2(i,j)*recip_hFac(i,j,k)
105			& recip_drF(k)recip_drC(k)/N2(i,j,k)
106			c3d(i,j,k) = fCori2(i,j)*recip_hFac(i,j,k)
107			& recip_drF(k)recip_drC(k+1)/N2(i,j,k+1)
108			b3d(i,j,k) = -a3d(i,j,k)-c3d(i,j,k)
109
110			ELSEIF (k.EQ.kLow(i,j)) THEN
111			IF (k.GT.Nr) THEN
112	jmc	1.2	PRINT '(a4,x,3(a1,i2),a1,3(x,i2))', 'k :',
113	m_bates	1.1	& '(',i,',',j,',',k,')',
114			& kLow(i,j), k, Nr
115			ENDIF
116			a3d(i,j,k) = fCori2(i,j)*recip_hFac(i,j,k)
117			& recip_drF(k)recip_drC(k)/N2(i,j,k)
118			c3d(i,j,k) = zeroRL
119			b3d(i,j,k) = -a3d(i,j,k)
120
121			ELSE
122			a3d(i,j,k) = zeroRL
123			b3d(i,j,k) = zeroRL
124			c3d(i,j,k) = zeroRL
125			ENDIF
126
127			ELSE
128			a3d(i,j,k) = zeroRL
129			b3d(i,j,k) = zeroRL
130			c3d(i,j,k) = zeroRL
131			ENDIF
132			ENDDO
133			ENDDO
134			ENDDO
135
136			#ifdef use_lapack
137	jmc	1.2
138			DO j=1-OLy,sNy+OLy
139			DO i=1-OLx,sNx+OLx
140	m_bates	1.1	IF (kLow(i,j).GT.0) THEN
141			DO kk=1,Nr
142			DO k=1,Nr
143			array(k,kk) = zeroRL
144			ENDDO
145			ENDDO
146	jmc	1.2
147	m_bates	1.1	k=1
148			array(k,k) = b3d(i,j,k)
149			array(k,k+1) = c3d(i,j,k)
150			DO k=2,Nr-1
151			array(k,k-1) = a3d(i,j,k)
152			array(k,k) = b3d(i,j,k)
153			array(k,k+1) = c3d(i,j,k)
154			ENDDO
155			k=Nr
156			array(k,k-1) = a3d(i,j,k)
157			array(k,k) = b3d(i,j,k)
158	jmc	1.2
159	m_bates	1.1	CALL DGEEV('N','V',Nr,array,Nr,eigR,eigI,dummy,1,
160			& vecs,Nr,work,Nr*Nr,info)
161	jmc	1.2	C Find the second largest eigenvalue (the Rossby radius)
162	m_bates	1.1	C and the first M baroclinic modes (eigenvectors)
163			DO m=0,GM_K3D_NModes
164			eigval(m) = -HUGE(zeroRL)
165			ENDDO
166	jmc	1.2
167	m_bates	1.1	DO k=1,kLow(i,j)
168			eigval(0) = MAX(eigval(0),eigR(k))
169			ENDDO
170			DO m=1,MIN(GM_K3D_NModes,klow(i,j))
171			C DO m=1,GM_K3D_NModes
172			DO k=1,kLow(i,j)
173			IF (eigR(k).LT.eigval(m-1)) THEN
174			eigval(m) = MAX(eigval(m),eigR(k))
175			IF (eigval(m).EQ.eigR(k)) idx=k
176			ENDIF
177			ENDDO
178			IF(vecs(1,idx).LT.zeroRL) THEN
179			DO k=1,Nr
180			vec(m,i,j,k) = -vecs(k,idx)
181			ENDDO
182			ELSE
183			DO k=1,Nr
184			vec(m,i,j,k) = vecs(k,idx)
185			ENDDO
186			ENDIF
187			ENDDO
188			val(i,j) = eigval(1)
189			ELSE
190			val(i,j)=zeroRL
191			DO k=1,Nr
192			DO m=1,GM_K3D_NModes
193			vec(m,i,j,k)=zeroRL
194			ENDDO
195			ENDDO
196	jmc	1.2
197	m_bates	1.1	ENDIF
198			ENDDO
199			ENDDO
200
201			C Normalise the eigenvectors
202	jmc	1.2	DO j=1-OLy,sNy+OLy
203			DO i=1-OLx,sNx+OLx
204	m_bates	1.1	DO m=1,GM_K3D_NModes
205			vecint(m,i,j) = zeroRL
206			ENDDO
207			ENDDO
208			ENDDO
209
210			DO k=1,Nr
211	jmc	1.2	DO j=1-OLy,sNy+OLy
212			DO i=1-OLx,sNx+OLx
213	m_bates	1.1	DO m=1,GM_K3D_NModes
214	jmc	1.2	vecint(m,i,j) = vecint(m,i,j) +
215	m_bates	1.1	& mask(i,j,k)drF(k)hfac(i,j,k)
216			& vec(m,i,j,k)vec(m,i,j,k)
217			ENDDO
218			ENDDO
219			ENDDO
220			ENDDO
221
222	jmc	1.2	DO j=1-OLy,sNy+OLy
223			DO i=1-OLx,sNx+OLx
224	m_bates	1.1	DO m=1,GM_K3D_NModes
225			vecint(m,i,j) = SQRT(vecint(m,i,j))
226			ENDDO
227			ENDDO
228			ENDDO
229
230			DO k=1,Nr
231	jmc	1.2	DO j=1-OLy,sNy+OLy
232			DO i=1-OLx,sNx+OLx
233	m_bates	1.1	DO m=1,GM_K3D_NModes
234			vec(m,i,j,k) = vec(m,i,j,k)/(vecint(m,i,j)+small)
235			ENDDO
236			ENDDO
237			ENDDO
238			ENDDO
239
240	jmc	1.2	#else
241	m_bates	1.1	C CALL EIGENVAL(a3d,b3d,c3d,bi,bj,val)
242			C CALL EIGENVEC(a3d,b3d,c3d,val,1,bi,bj,vec)
243			#endif
244
245			DO k=1,Nr
246	jmc	1.2	DO j=1-OLy,sNy+OLy
247			DO i=1-OLx,sNx+OLx
248	m_bates	1.1	IF (kLow(i,j).GT.2 .AND. val(i,j).NE.zeroRL) THEN
249			Kdef(i,j) = SQRT(ABS(val(i,j)))
250			ELSE
251			Kdef(i,j) = zeroRL
252			ENDIF
253			ENDDO
254			ENDDO
255			ENDDO
256
257			#ifdef ALLOW_DIAGNOSTICS
258			C Diagnostics
259			IF ( useDiagnostics.AND.writediag ) THEN
260			CALL DIAGNOSTICS_FILL(a3d, 'GM_A3D ',0,Nr,0,1,1,myThid)
261			CALL DIAGNOSTICS_FILL(b3d, 'GM_B3D ',0,Nr,0,1,1,myThid)
262			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)
264			DO k=1,Nr
265	jmc	1.2	DO j=1-OLy,sNy+OLy
266			DO i=1-OLx,sNx+OLx
267	m_bates	1.1	vec_diag(i,j,k) = vec(1,i,j,k)
268			ENDDO
269			ENDDO
270			ENDDO
271
272			CALL DIAGNOSTICS_FILL(vec_diag, 'GM_VEC ',0,Nr,0,1,1,myThid)
273			ENDIF
274			#endif
275	jmc	1.2
276			#endif /* GM_K3D */
277
278			RETURN
279	m_bates	1.1	END