1 |
C $Header: /u/gcmpack/MITgcm/pkg/gmredi/gmredi_calc_eigs.F,v 1.2 2013/06/22 14:44:54 jmc Exp $ |
2 |
C $Name: $ |
3 |
|
4 |
#include "GMREDI_OPTIONS.h" |
5 |
|
6 |
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 rlow, nmodes, writediag, |
13 |
O Rmid, vec) |
14 |
|
15 |
C !DESCRIPTION: \bv |
16 |
C *==========================================================* |
17 |
C | SUBROUTINE GMREDI_CALC_URMS |
18 |
C | o Calculate the vertical structure of the rms eddy |
19 |
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 |
INTEGER nmodes |
40 |
INTEGER kLow(1-Olx:sNx+Olx,1-Oly:sNy+Oly) |
41 |
_RL mask(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr) |
42 |
_RL hfac(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr) |
43 |
_RL recip_hfac(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr) |
44 |
_RL rlow(1-Olx:sNx+Olx,1-Oly:sNy+Oly) |
45 |
_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 |
50 |
C !LOCAL VARIABLES: |
51 |
C == Local variables == |
52 |
INTEGER i,j,k,kk,m |
53 |
# ifdef use_lapack |
54 |
C info :: error code from LAPACK |
55 |
C idx :: index used for sorting the eigenvalues |
56 |
C a3d :: lower diagonal of eigenvalue problem |
57 |
C b3d :: diagonal of eigenvalue problem |
58 |
C c3d :: upper diagonal of eigenvalue problem |
59 |
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) |
74 |
_RL array(Nr,Nr) |
75 |
_RL eigval(0:nmodes) |
76 |
# else |
77 |
C drNr :: distance from bottom of cell to cell centre at Nr |
78 |
C BuoyFreq :: buoyancy frequency, SQRT(N2) |
79 |
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) |
99 |
|
100 |
C Square of the Coriolis parameter |
101 |
DO i=1-OLx,sNx+OLx |
102 |
DO j=1-OLy,sNy+OLy |
103 |
fCori2(i,j) = fCori(i,j,bi,bj)*fCori(i,j,bi,bj) |
104 |
ENDDO |
105 |
ENDDO |
106 |
|
107 |
DO k=1,Nr |
108 |
DO j=1-Oly,sNy+Oly |
109 |
DO i=1-Olx,sNx+Olx |
110 |
DO m=1,nmodes |
111 |
vec(m,i,j,k) = zeroRL |
112 |
ENDDO |
113 |
ENDDO |
114 |
ENDDO |
115 |
ENDDO |
116 |
|
117 |
# ifdef use_lapack |
118 |
C Calculate the tridiagonal operator matrix for |
119 |
C f^2 d/dz 1/N^2 d/dz |
120 |
C a3d is the lower off-diagonal, b3d is the diagonal |
121 |
C and c3d is the upper off-diagonal |
122 |
DO k=1,Nr |
123 |
DO j=1-OLy,sNy+OLy |
124 |
DO i=1-OLx,sNx+OLx |
125 |
IF (kLow(i,j) .GT. 0) THEN |
126 |
IF (k.EQ.1) THEN |
127 |
a3d(i,j,k) = zeroRL |
128 |
c3d(i,j,k) = fCori2(i,j)*recip_hFac(i,j,k) |
129 |
& *recip_drC(k+1)*recip_drF(k)/N2(i,j,k+1) |
130 |
b3d(i,j,k) = -c3d(i,j,k) |
131 |
|
132 |
ELSEIF (k.LT.kLow(i,j)) THEN |
133 |
a3d(i,j,k) = fCori2(i,j)*recip_hFac(i,j,k) |
134 |
& *recip_drF(k)*recip_drC(k)/N2(i,j,k) |
135 |
c3d(i,j,k) = fCori2(i,j)*recip_hFac(i,j,k) |
136 |
& *recip_drF(k)*recip_drC(k+1)/N2(i,j,k+1) |
137 |
b3d(i,j,k) = -a3d(i,j,k)-c3d(i,j,k) |
138 |
|
139 |
ELSEIF (k.EQ.kLow(i,j)) THEN |
140 |
a3d(i,j,k) = fCori2(i,j)*recip_hFac(i,j,k) |
141 |
& *recip_drF(k)*recip_drC(k)/N2(i,j,k) |
142 |
c3d(i,j,k) = zeroRL |
143 |
b3d(i,j,k) = -a3d(i,j,k) |
144 |
|
145 |
ELSE |
146 |
a3d(i,j,k) = zeroRL |
147 |
b3d(i,j,k) = zeroRL |
148 |
c3d(i,j,k) = zeroRL |
149 |
ENDIF |
150 |
|
151 |
ELSE |
152 |
a3d(i,j,k) = zeroRL |
153 |
b3d(i,j,k) = zeroRL |
154 |
c3d(i,j,k) = zeroRL |
155 |
ENDIF |
156 |
ENDDO |
157 |
ENDDO |
158 |
ENDDO |
159 |
|
160 |
DO j=1-Oly,sNy+Oly |
161 |
DO i=1-Olx,sNx+Olx |
162 |
IF (kLow(i,j).GT.0) THEN |
163 |
DO kk=1,Nr |
164 |
DO k=1,Nr |
165 |
array(k,kk) = zeroRL |
166 |
ENDDO |
167 |
ENDDO |
168 |
|
169 |
k=1 |
170 |
array(k,k) = b3d(i,j,k) |
171 |
array(k,k+1) = c3d(i,j,k) |
172 |
DO k=2,Nr-1 |
173 |
array(k,k-1) = a3d(i,j,k) |
174 |
array(k,k) = b3d(i,j,k) |
175 |
array(k,k+1) = c3d(i,j,k) |
176 |
ENDDO |
177 |
k=Nr |
178 |
array(k,k-1) = a3d(i,j,k) |
179 |
array(k,k) = b3d(i,j,k) |
180 |
|
181 |
CALL DGEEV('N','V',Nr,array,Nr,eigR,eigI,dummy,1, |
182 |
& vecs,Nr,work,Nr*Nr,info) |
183 |
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) |
199 |
DO m=0,nmodes |
200 |
eigval(m) = -HUGE(zeroRL) |
201 |
ENDDO |
202 |
|
203 |
DO k=1,kLow(i,j) |
204 |
eigval(0) = MAX(eigval(0),eigR(k)) |
205 |
ENDDO |
206 |
DO m=1,MIN(nmodes,klow(i,j)) |
207 |
DO k=1,kLow(i,j) |
208 |
IF (eigR(k).LT.eigval(m-1)) THEN |
209 |
eigval(m) = MAX(eigval(m),eigR(k)) |
210 |
IF (eigval(m).EQ.eigR(k)) idx=k |
211 |
ENDIF |
212 |
ENDDO |
213 |
IF(vecs(1,idx).LT.zeroRL) THEN |
214 |
DO k=1,Nr |
215 |
vec(m,i,j,k) = -vecs(k,idx) |
216 |
ENDDO |
217 |
ELSE |
218 |
DO k=1,Nr |
219 |
vec(m,i,j,k) = vecs(k,idx) |
220 |
ENDDO |
221 |
ENDIF |
222 |
ENDDO |
223 |
val(i,j) = eigval(1) |
224 |
ELSE |
225 |
val(i,j)=zeroRL |
226 |
DO k=1,Nr |
227 |
DO m=1,nmodes |
228 |
vec(m,i,j,k)=zeroRL |
229 |
ENDDO |
230 |
ENDDO |
231 |
|
232 |
ENDIF |
233 |
ENDDO |
234 |
ENDDO |
235 |
|
236 |
DO k=1,Nr |
237 |
DO j=1-Oly,sNy+Oly |
238 |
DO i=1-Olx,sNx+Olx |
239 |
IF (kLow(i,j).GT.2 .AND. val(i,j).NE.zeroRL) THEN |
240 |
Rmid(i,j) = 1.0/(SQRT(ABS(val(i,j)))+small) |
241 |
ELSE |
242 |
Rmid(i,j) = zeroRL |
243 |
ENDIF |
244 |
ENDDO |
245 |
ENDDO |
246 |
ENDDO |
247 |
|
248 |
# else |
249 |
DO k=1,Nr |
250 |
DO j=1-Oly,sNy+Oly |
251 |
DO i=1-Olx,sNx+Olx |
252 |
BuoyFreq(i,j,k) = mask(i,j,k)*SQRT(N2(i,j,k)) |
253 |
ENDDO |
254 |
ENDDO |
255 |
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 |
349 |
ENDDO |
350 |
ENDDO |
351 |
ENDDO |
352 |
# endif |
353 |
|
354 |
C Normalise the eigenvectors |
355 |
DO j=1-Oly,sNy+Oly |
356 |
DO i=1-Olx,sNx+Olx |
357 |
DO m=1,nmodes |
358 |
vecint(m,i,j) = zeroRL |
359 |
ENDDO |
360 |
ENDDO |
361 |
ENDDO |
362 |
|
363 |
DO k=1,Nr |
364 |
DO j=1-Oly,sNy+Oly |
365 |
DO i=1-Olx,sNx+Olx |
366 |
DO m=1,nmodes |
367 |
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 |
371 |
ENDDO |
372 |
ENDDO |
373 |
ENDDO |
374 |
|
375 |
DO j=1-Oly,sNy+Oly |
376 |
DO i=1-Olx,sNx+Olx |
377 |
DO m=1,nmodes |
378 |
vecint(m,i,j) = SQRT(vecint(m,i,j)) |
379 |
ENDDO |
380 |
ENDDO |
381 |
ENDDO |
382 |
|
383 |
DO k=1,Nr |
384 |
DO j=1-Oly,sNy+Oly |
385 |
DO i=1-Olx,sNx+Olx |
386 |
DO m=1,nmodes |
387 |
vec(m,i,j,k) = vec(m,i,j,k)/(vecint(m,i,j)+small) |
388 |
ENDDO |
389 |
ENDDO |
390 |
ENDDO |
391 |
ENDDO |
392 |
|
393 |
# ifdef ALLOW_DIAGNOSTICS |
394 |
C Diagnostics |
395 |
IF ( useDiagnostics.AND.writediag ) THEN |
396 |
# ifdef use_lapack |
397 |
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) |
399 |
CALL DIAGNOSTICS_FILL(c3d, 'GM_C3D ',0,Nr,0,1,1,myThid) |
400 |
# endif |
401 |
ENDIF |
402 |
# endif |
403 |
|
404 |
#endif /* GM_K3D */ |
405 |
|
406 |
RETURN |
407 |
END |