/[MITgcm]/MITgcm_contrib/pvcode/fortran/calc_qgpv.F
ViewVC logotype

Annotation of /MITgcm_contrib/pvcode/fortran/calc_qgpv.F

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


Revision 1.1 - (hide annotations) (download)
Sat May 27 02:52:26 2006 UTC (19 years, 1 month ago) by cnh
Branch point for: MAIN, cnh
Initial revision

1 cnh 1.1 #include "CPP_EEOPTIONS.h"
2    
3     CStartOfInterFace
4     SUBROUTINE CALC_QGPV(
5     I bi, bj, iMin,iMax,jMin,jMax, pH, nsquare,
6     I K13, K23,
7     I myTime, myIter,
8     O q, qr, qs, KpvU, KpvV, gUpv, gVpv, dqdyU,
9     I myThid )
10    
11     C /==========================================================\
12     C | SUBROUTINE CALC_QGPV |
13     C | o Calculate the quasigeostrophic potential vorticity |
14     C | o and eddy flux terms of qgpv etc |
15     C |==========================================================|
16     C | Richard Wardle 7/98 |
17     C | Daniel Jamous 1/00 N^2 entered as an argument |
18     C | instead of being calculated from |
19     C | temperature |
20     C | Updated to take into account |
21     C | arbitrary topography |
22     C | 3/00 compute PV fluxes using a local |
23     C | N^2 (as opposed to an horizontal |
24     C | average)--see key N2local |
25     C | IMPORTANT NOTE: For now, the routine has been tested |
26     C | only when there are no interpolation of gsf and dTdz |
27     C | at horizontal boundaries, when qgpv is only made of |
28     C | qs the stretching term, and when the Kpv are constants. |
29     C | For a more general case, some more coding might be |
30     C | necessary. |
31     C \==========================================================/
32     C
33     C == Global variables ==
34     #include "SIZE.h"
35     #include "DYNVARS.h"
36     #include "GRID.h"
37     #include "EEPARAMS.h"
38     #include "PARAMS.h"
39     #include "CG2D.h"
40     #include "FFIELDS.h"
41     C
42     C ========= Local variables that we may want to export =============
43     C q - quasigeostrophic potential vorticity
44     C ==================================================================
45     C
46     _RL pH (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
47     _RL q (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
48     _RL qf (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
49     _RL qr (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
50     _RL qs (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
51     _RL KpvV (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
52     _RL KpvU (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
53     _RL gUpv (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
54     _RL gVpv (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
55     _RL nsquare (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
56     _RL K13 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
57     _RL K23 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
58     C
59     INTEGER bi,bj,iMin,iMax,jMin,jMax
60     INTEGER myThid
61     INTEGER myIter
62     _RL myTime
63     C
64     CEndOfInterface
65     C
66     C define local variables here:
67     _RL ptotal (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
68     _RL sfn (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
69     _RL gsfn (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
70     _RL tmpphi (lShare8,MAX_NO_THREADS)
71     _RL tmpphi_area (lShare8,MAX_NO_THREADS)
72     c _RL tbarxy (Nr)
73     _RL sfnbarxy (Nr)
74     _RL dTdz (Nr)
75     _RL dTdz3d (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
76     _RL dVdx (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
77     _RL dUdy (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
78     _RL dVdxbarxy (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
79     _RL dUdybarxy (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
80     _RL dgsfndz (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
81     _RL arr (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr+1)
82     _RL arrprime (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr+1)
83     _RL arr2 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
84     C
85     _RL dqdyV (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
86     _RL dqdyT (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
87     _RL dqdyU (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
88     C
89     _RL dqdxU (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
90     _RL dqdxT (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
91     _RL dqdxV (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
92     C
93     _RL pvFacT (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
94     _RL pvFacU (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
95     _RL pvFacV (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
96     C
97     _RL d2gsfndz2 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
98     _RL interp1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
99     _RL d2Tdz2 (Nr)
100     _RL interp2 (Nr)
101     C
102     _RL pmask (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
103     _RL umask (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
104     _RL vmask (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
105     C
106     _RL num_i (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
107     _RL num_b (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
108     _RL denom (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
109     _RL kbotindex (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
110    
111     INTEGER i,j,k
112     INTEGER iG, jG
113     INTEGER ku, kv, kup1, kbot
114     ctest
115     CHARACTER*(MAX_LEN_MBUF) suff
116     ctest
117     _RL tmpphiS
118     _RL Kpvref
119     _RL tmp, fac
120     _RL mldu, mldv
121     _RL int1, int2, int3, int4, int5, int6
122     _RL int2km1, signtest
123    
124     C ==================================================================
125     C-- iG and jG are the global indices
126     C ==================================================================
127     Kpvref = 1000.
128     C ==================================================================
129     C-- Compute mask at pressure and velocity points
130     C ==================================================================
131    
132     DO k=1,Nr
133     DO j=iMin,jMax
134     DO i=jMin,iMax
135     pmask(i,j,k) = 1.
136     IF (_hFacC(i,j,k,bi,bj).eq.0.) pmask(i,j,k)=0.
137     ENDDO
138     ENDDO
139     DO j=iMin,jMax
140     DO i=jMin,iMax
141     umask(i,j,k) = pmask(i-1,j,k)*pmask(i,j,k)
142     vmask(i,j,k) = pmask(i,j-1,k)*pmask(i,j,k)
143     ENDDO
144     ENDDO
145     ENDDO
146    
147     DO j=iMin,jMax
148     DO i=jMin,iMax
149     kbotindex(i,j) = 0.
150     IF ( pmask(i,j,Nr) .eq. 1. ) THEN
151     kbotindex(i,j) = float(Nr)
152     ELSE
153     DO K = Nr-1,1,-1
154     IF ( pmask(i,j,k+1) .EQ. 0. .AND.
155     & pmask(i,j,k) .EQ. 1. ) THEN
156     kbotindex(i,j) = float(k)
157     ENDIF
158     ENDDO
159     ENDIF
160     ENDDO
161     ENDDO
162    
163     C
164     C ==================================================================
165     C Compute the geostrophic streamfunction: gsf ======================
166     C ==================================================================
167     C
168     C -------v--------
169     C | |
170     C | |
171     C u x u Streamfunction is located at p points
172     C | gsfn |
173     C | |
174     C -------v--------
175     C
176     C ==================================================================
177     C Compute the total pressure field:
178     DO k=1,Nr
179     DO j=iMin,jMax
180     DO i=jMin,iMax
181     ptotal(i,j,k) = pH(i,j,k) +
182     & cg2d_x(i,j,bi,bj) * (gBaro * rhonil)
183     ENDDO
184     ENDDO
185     ENDDO
186     C Compute the streamfunction:
187     DO k=1,Nr
188     DO j=iMin,jMax
189     DO i=iMin,iMax
190     sfn(i,j,k) = ptotal(i,j,k) / ( rhonil * f0 )
191     ENDDO
192     ENDDO
193     ENDDO
194     C
195     C Compute the streamfunction:
196     DO k=1,Nr
197     DO j=iMin,jMax
198     DO i=iMin,iMax
199     sfn(i,j,k) = ptotal(i,j,k) / ( rhonil * f0 )
200     ENDDO
201     ENDDO
202     ENDDO
203     C ==================================================================
204     C 1. Evaluate the global horizontal mean of sf:
205     c write(0,*)'Evaluate the horizontal mean'
206     do k=1,Nr
207     tmpphi(1,myThid)=0.
208     tmpphi_area(1,myThid)=0.
209     do j=1,sNy
210     do i=1,sNx
211     tmpphi(1,myThid) = tmpphi(1,myThid) + sfn(i,j,k)*
212     & _rA(i,j,bi,bj)*pmask(i,j,k)
213     tmpphi_area(1,myThid) = tmpphi_area(1,myThid) +
214     & _rA(i,j,bi,bj)*pmask(i,j,k)
215     enddo
216     enddo
217     _GLOBAL_SUM_R8( tmpphi, myThid )
218     _GLOBAL_SUM_R8( tmpphi_area, myThid )
219     sfnbarxy(k)=tmpphi(1,myThid)/tmpphi_area(1,myThid)
220     c write(0,*)'level :',k,' horizontal mean :',sfnbarxy(k)
221     enddo
222     C
223     C 2. Subtract sgbarxy from sf to give the geostrophic streamfunction:
224     DO k=1,Nr
225     DO j=jMin,jMax
226     DO i=iMin,iMax
227     gsfn(i,j,k) = ( sfn(i,j,k) - sfnbarxy(k) ) * pmask(i,j,k)
228     ENDDO
229     ENDDO
230     ENDDO
231     C ==================================================================
232     C-- Interpolate gsf at horizontal boundaries:
233     C ==================================================================
234     C-- Vertical gradient of gsfn:
235     DO k=2,Nr
236     DO j=jMin,jMax
237     DO i=iMin,iMax
238     dgsfndz(i,j,k) = recip_drC(k)*( gsfn(i,j,k-1)-gsfn(i,j,k) )
239     ENDDO
240     ENDDO
241     ENDDO
242     C-- Vertical gradient of dgsfndz:
243     DO k=2,Nr
244     DO j=jMin,jMax
245     DO i=iMin,iMax
246     d2gsfndz2(i,j,k) = recip_drC(k)*(dgsfndz(i,j,k)-dgsfndz(i,j,k+1))
247     ENDDO
248     ENDDO
249     ENDDO
250     C--
251     DO k=2,3
252     DO j=jMin,jMax
253     DO i=iMin,iMax
254     d2gsfndz2(i,j,k) = d2gsfndz2(i,j,4)
255     ENDDO
256     ENDDO
257     ENDDO
258     C
259     C-- Evaluate the gradients of gsfn
260     DO j=jMin,jMax
261     DO i=iMin,iMax
262     interp1(i,j,3) = dgsfndz(i,j,4) + drC(3)*d2gsfndz2(i,j,3)
263     ENDDO
264     ENDDO
265     C
266     DO j=jMin,jMax
267     DO i=iMin,iMax
268     interp1(i,j,2) = interp1(i,j,3) + drC(2)*d2gsfndz2(i,j,2)
269     ENDDO
270     ENDDO
271     DO k=4,Nr
272     DO j=jMin,jMax
273     DO i=iMin,iMax
274     interp1(i,j,k) = dgsfndz(i,j,k)
275     ENDDO
276     ENDDO
277     ENDDO
278     C ==================================================================
279     C ==================================================================
280     C === Global mean temperature profile: dTdz3d =====================
281     C ==================================================================
282     C Evaluate the global horizontal mean profile of nsquare
283     c write(0,*)'mean N2'
284     do k=1,Nr
285     tmpphi(1,myThid)=0.
286     tmpphi_area(1,myThid)=0.
287     do j=1,sNy
288     do i=1,sNx
289     tmpphi(1,myThid)=tmpphi(1,myThid)+nsquare(i,j,k,bi,bj)*
290     & _rA(i,j,bi,bj)*pmask(i,j,k)
291     tmpphi_area(1,myThid)=tmpphi_area(1,myThid) +
292     & _rA(i,j,bi,bj)*pmask(i,j,k)
293     enddo
294     enddo
295     _GLOBAL_SUM_R8( tmpphi, myThid )
296     _GLOBAL_SUM_R8( tmpphi_area, myThid )
297     dTdz(K)=tmpphi(1,myThid)/tmpphi_area(1,myThid)
298     c write(0,*)'level :',k,' horizontal mean :',dTdz(k)
299     enddo
300     C
301     C Vertical gradient of tbarxy
302     c DO k=2,Nr
303     c dTdz(k) = recip_drC(k)*( tbarxy(k-1)-tbarxy(k) )
304     c ENDDO
305     C Avoid dividing by zero:
306     dTdz(1) = dTdz(2)
307     C ==================================================================
308     C interpolate dTdz at horizontal boundaries:
309     C ==================================================================
310     C-- Vertical gradient of dTdz:
311     DO k=2,Nr
312     DO j=jMin,jMax
313     DO i=iMin,iMax
314     d2Tdz2(k) = recip_drC(k)*(dTdz(k)-dTdz(k+1))
315     ENDDO
316     ENDDO
317     ENDDO
318     C--
319     DO k=2,3
320     DO j=jMin,jMax
321     DO i=iMin,iMax
322     d2Tdz2(k) = d2Tdz2(4)
323     ENDDO
324     ENDDO
325     ENDDO
326     C
327     C-- Evaluate the gradients of gsfn
328     DO j=jMin,jMax
329     DO i=iMin,iMax
330     interp2(3) = dTdz(4) + drC(3)*d2Tdz2(3)
331     ENDDO
332     ENDDO
333     C
334     DO j=jMin,jMax
335     DO i=iMin,iMax
336     interp2(2) = interp2(3) + drC(2)*d2Tdz2(2)
337     ENDDO
338     ENDDO
339     C
340     DO k=4,Nr
341     DO j=jMin,jMax
342     DO i=iMin,iMax
343     interp2(k) = dTdz(k)
344     ENDDO
345     ENDDO
346     ENDDO
347     C
348     C ==================================================================
349     C Create a threedimensional array of dTdz:
350     DO k=1,Nr
351     DO j=jMin,jMax
352     DO i=iMin,iMax
353     dTdz3d(i,j,k) = dTdz(k)
354     c dTdz3d(i,j,k) = interp2(k)
355     ENDDO
356     ENDDO
357     ENDDO
358     C ==================================================================
359     C ==================================================================
360     C ======== Calculate the quasigeostrophic potential vorticity ======
361     C ==================================================================
362     C
363     C q = quasigeostrophic potential vorticity
364     C
365     C -------v--------
366     C | |
367     C | |
368     C u x u qgpv is located at p points
369     C | qgpv |
370     C | |
371     C -------v--------
372     C
373     C q = f + (gsf)_xx + (gsf)_yy + fo^2 ( (gsf)_z / N^2 )_z
374     C
375     C = qf + qr + qs
376     C
377     C ==============================
378     C 1: Coriolis
379     C ==============================
380     C although fCori is defined at u-points they are at the same
381     C latitude as p-points and so do not need to be inerpolated:
382     DO j=jMin,jMax
383     DO i=iMin,iMax
384     qf(i,j) = fCori(i,j,bi,bj)
385     ENDDO
386     ENDDO
387     C
388     C ==============================
389     C 2: relative
390     C ==============================
391     C !! for now use gradients of velocity field !!
392     C ==============================
393     C-- Zonal gradient of vVel
394     C ==============================
395     DO k=1,Nr
396     DO j=jMin,jMax
397     DO i=iMin,iMax
398     dVdx(i,j,k) = _recip_dxV(i,j,bi,bj) *
399     & ( vVel(i,j,k,bi,bj) - vVel(i-1,j,k,bi,bj) )
400     ENDDO
401     ENDDO
402     ENDDO
403     C
404     C-- Interpolate term to center point
405     DO k=1,Nr
406     DO j=jMin,jMax
407     DO i=iMin,iMax
408     iG = myXGlobalLo-1+(bi-1)*sNx+I
409     C
410     dVdxbarxy(i,j,k) = 0.25 *
411     & ( dVdx(i,j ,k) + dVdx(i+1,j ,k)
412     & + dVdx(i,j+1,k) + dVdx(i+1,j+1,k) )
413     C
414     C-- free-slip b.c.'s:
415     IF ( iG .EQ. 1 ) THEN
416     dVdxbarxy(i,j,k) = 0.5 *
417     & ( dVdx(i+1,j ,k) + dVdx(i+1,j+1,k) )
418     ENDIF
419     IF ( iG .EQ. Nx-1 ) THEN
420     dVdxbarxy(i,j,k) = 0.5 *
421     & ( dVdx(i,j,k) + dVdx(i,j+1,k) )
422     ENDIF
423     C
424     ENDDO
425     ENDDO
426     ENDDO
427     C--------------------------------------------------------
428     C ==============================
429     C-- Meridional gradient of uVel
430     C ==============================
431     DO k=1,Nr
432     DO j=jMin,jMax
433     DO i=iMin,iMax
434     dUdy(i,j,k) = _recip_dyU(i,j,bi,bj) *
435     & ( uVel(i,j,k,bi,bj) - uVel(i,j-1,k,bi,bj) )
436     ENDDO
437     ENDDO
438     ENDDO
439     C
440     C-- Interpolate term to center point
441     DO k=1,Nr
442     DO j=jMin,jMax
443     DO i=iMin,iMax
444     jG = myYGlobalLo-1+(bj-1)*sNy+J
445     C
446     dUdybarxy(i,j,k) = 0.25 *
447     & ( dUdy(i,j ,k) + dUdy(i+1,j ,k)
448     & + dUdy(i,j+1,k) + dUdy(i+1,j+1,k) )
449     C
450     C-- free-slip b.c.'s
451     IF ( jG .EQ. 1 ) THEN
452     dUdybarxy(i,j,k) = 0.5 *
453     & ( dUdy(i, j+1,k) + dUdy(i+1,j+1,k) )
454     ENDIF
455     IF ( jG .EQ. Ny-1 ) THEN
456     dUdybarxy(i,j,k) = 0.5 *
457     & ( dUdy(i, j,k) + dUdy(i+1,j,k) )
458     ENDIF
459     C
460     ENDDO
461     ENDDO
462     ENDDO
463     C-- =================================
464     DO k=1,Nr
465     DO j=jMin,jMax
466     DO i=iMin,iMax
467     qr(i,j,k) = ( dVdxbarxy(i,j,k) - dUdybarxy(i,j,k) )
468     ENDDO
469     ENDDO
470     ENDDO
471     C-- =================================
472     C 3. stretching
473     C-- =================================
474     C-- ====a. Vertical gradient of gsfn:
475     DO k=2,Nr
476     DO j=jMin,jMax
477     DO i=iMin,iMax
478     dgsfndz(i,j,k) = recip_drC(k)*( gsfn(i,j,k-1)-gsfn(i,j,k) )*pmask(i,j,k)
479     ENDDO
480     ENDDO
481     ENDDO
482     C
483     C-- ====b. Evaluate (gsf)_z / N^2
484     DO k=2,Nr
485     DO j=jMin,jMax
486     DO i=iMin,iMax
487     c arr(i,j,k) = interp1(i,j,k)
488     arr(i,j,k) = dgsfndz(i,j,k)
489     & / dTdz3d(i,j,k)
490     ENDDO
491     ENDDO
492     ENDDO
493     C set to zero in the upper and lower layers:
494     C to include the pv sheets
495     DO j=jMin,jMax
496     DO i=iMin,iMax
497     arr(i,j,1) = 0.
498     arr(i,j,Nr+1) = 0.
499     ENDDO
500     ENDDO
501     C
502     C-- ====c. Evaluate z-derivative of (gsf)_z / N^2
503     DO k=1,Nr
504     DO j=jMin,jMax
505     DO i=iMin,iMax
506     arr2(i,j,k) = recip_drF(k)*( arr(i,j,k)-arr(i,j,k+1) )
507     ENDDO
508     ENDDO
509     ENDDO
510     C
511     C-- ====d. mulitply by f0^2
512     DO k=1,Nr
513     DO j=jMin,jMax
514     DO i=iMin,iMax
515     qs(i,j,k) = ( f0**2. ) * arr2(i,j,k) * pmask(i,j,k)
516     ENDDO
517     ENDDO
518     ENDDO
519     C ==================================================================
520     C ==================================================================
521     C Sum for q: the quasigeostrophic potential vorticity
522     C ==================================================================
523     DO k=1,Nr
524     DO j=jMin,jMax
525     DO i=iMin,iMax
526     c q(i,j,k)= qf(i,j) + qr(i,j,k) + qs(i,j,k)
527     q(i,j,k) = qs(i,j,k)
528     ENDDO
529     ENDDO
530     ENDDO
531     c write(0,*) (q(10,j,2),j=1,sNy)
532    
533     C ==================================================================
534     C ==================================================================
535     C--
536     IF ( pvForcing ) THEN
537     C--
538     C-- =================================================================
539     C-- ================Compute the PV gradients=========================
540     C-- =================================================================
541     C--
542     IF ( N2local ) THEN
543     C
544     C Interpolate K13 and K23 at horizontal u and v positions but vertical
545     C w level
546     DO k=1,Nr
547     DO j=jMin,jMax
548     DO i=iMin,iMax
549     arr(i,j,k) = - 0.5 *
550     & ( K23(i-1,j,k)+K23(i,j,k) ) * umask(i,j,k)
551     arrprime(i,j,k) = - 0.5 *
552     & ( K13(i,j-1,k)+K13(i,j,k) ) * vmask(i,j,k)
553     ENDDO
554     ENDDO
555     ENDDO
556     DO j=jMin,jMax
557     DO i=iMin,iMax
558     arr(i,j,Nr+1) = 0.
559     arrprime(i,j,Nr+1) = 0.
560     ENDDO
561     ENDDO
562     C Take z-derivative of arr and arrprime
563     DO k=1,Nr
564     DO j=jMin,jMax
565     DO i=iMin,iMax
566     dqdyU(i,j,k) = fCori(i,j,bi,bj) * recip_drF(k) *
567     & (arr(i,j,k)-arr(i,j,k+1))
568     dqdxV(i,j,k) = 0.5 * ( fCori(i,j-1,bi,bj)+fCori(i,j,bi,bj) ) * recip_drF(k) *
569     & (arrprime(i,j,k)-arrprime(i,j,k+1))
570     IF ( PVsheetmld ) THEN
571     c ku = int( min(mldindex(i-1,j),mldindex(i,j)) )
572     c mldu = min( mld(i-1,j),mld(i,j) )
573     ku = int( max(mldindex(i-1,j),mldindex(i,j)) )
574     mldu = max( mld(i-1,j),mld(i,j) )
575     IF ( k .le. ku .AND. mldu .ne. 0. )
576     & dqdyU(i,j,k) = fCori(i,j,bi,bj)*(-arr(i,j,ku+1))/mldu
577     c kv = int( min(mldindex(i,j-1),mldindex(i,j)) )
578     c mldv = min( mld(i,j-1),mld(i,j) )
579     kv = int( max(mldindex(i,j-1),mldindex(i,j)) )
580     mldv = max( mld(i,j-1),mld(i,j) )
581     IF ( k .le. kv .AND. mldv .ne. 0. )
582     & dqdxV(i,j,k) = 0.5*(fCori(i,j-1,bi,bj)+fCori(i,j,bi,bj))
583     & *(-arrprime(i,j,kv+1))/mldv
584     ENDIF
585     IF ( betaterm ) THEN
586     dqdyU(i,j,k) = dqdyU(i,j,k) + dfdy(i,j,bi,bj)*umask(i,j,k)
587     ENDIF
588     ENDDO
589     ENDDO
590     ENDDO
591     C
592     ELSE
593     C
594     C ========================
595     C-- Meridional gradient of q
596     C ========================
597     c write(0,*)'Meridional gradient of q'
598     c-- At V-points:
599     DO k=1,Nr
600     DO j=jMin,jMax
601     DO i=iMin,iMax
602     dqdyV(i,j,k) = _recip_dyU(i,j,bi,bj)
603     & * (q(i,j,k)-q(i,j-1,k))
604     ENDDO
605     ENDDO
606     ENDDO
607     C
608     C-- ----------------------------
609     C-- Interpolate term to U points
610     C-- ----------------------------
611     C-- ::T-points first
612     DO k=1,Nr
613     DO j=jMin,jMax
614     DO i=iMin,iMax
615     jG = myYGlobalLo-1+(bj-1)*sNy+J
616     dqdyT(i,j,k) = 0.
617     tmp = vmask(i,j,k) + vmask(i,j+1,k)
618     IF ( tmp .ne. 0. ) THEN
619     dqdyT(i,j,k) = ( dqdyV(i,j,k)*vmask(i,j,k) + dqdyV(i,j+1,k)*vmask(i,j+1,k) )/( vmask(i,j,k) + vmask(i,j+1,k) )
620     ENDIF
621     C
622     C-- bondaries:
623     c IF ( jG .EQ. 1 ) THEN
624     c dqdyT(i,j,k) = dqdyV(i, j+1,k)
625     c ENDIF
626     c IF ( jG .EQ. sNy-1 ) THEN
627     c dqdyT(i,j,k) = dqdyV(i,j,k)
628     c ENDIF
629     C--
630     ENDDO
631     ENDDO
632     ENDDO
633     C-- ::U-points
634     DO k=1,Nr
635     DO j=jMin,jMax
636     DO i=iMin,iMax
637     iG = myXGlobalLo-1+(bi-1)*sNx+I
638     dqdyU(i,j,k) = 0.
639     IF ( umask(i,j,k) .ne. 0. ) THEN
640     dqdyU(i,j,k) = ( dqdyT(i,j,k)*pmask(i,j,k) + dqdyT(i-1,j,k)*pmask(i-1,j,k) )/( pmask(i,j,k) + pmask(i-1,j,k) )
641     ENDIF
642     C
643     C-- bondaries:
644     c IF ( iG .EQ. 1 ) THEN
645     c dqdyU(i,j,k) = dqdyT(i+1,j,k)
646     c ENDIF
647     c IF ( iG .EQ. sNx-1 ) THEN
648     c dqdyU(i,j,k) = dqdyT(i,j,k)
649     c ENDIF
650     C--
651     ENDDO
652     ENDDO
653     ENDDO
654    
655     C ===================
656     C-- Zonal gradient of q
657     C ===================
658     c-- At U-points:
659     c write(0,*)'Zonal gradient of q'
660     DO k=1,Nr
661     DO j=jMin,jMax
662     DO i=iMin,iMax
663     dqdxU(i,j,k) = _recip_dxC(i,j,bi,bj)
664     & * ( q(i,j,k)-q(i-1,j,k) )
665     ENDDO
666     ENDDO
667     ENDDO
668     C
669     C-- ----------------------------
670     C-- Interpolate term to V points
671     C-- ----------------------------
672     C-- ::T-points first
673     DO k=1,Nr
674     DO j=jMin,jMax
675     DO i=iMin,iMax
676     iG = myXGlobalLo-1+(bi-1)*sNx+I
677     dqdxT(i,j,k) = 0.
678     tmp = umask(i,j,k) + umask(i+1,j,k)
679     IF ( tmp .ne. 0. ) THEN
680     dqdxT(i,j,k) = ( dqdxU(i,j,k)*umask(i,j,k) + dqdxU(i+1,j,k)*umask(i+1,j,k) )/( umask(i,j,k) + umask(i+1,j,k) )
681     ENDIF
682     C
683     C-- bondaries:
684     c IF ( iG .EQ. 1 ) THEN
685     c dqdxT(i,j,k) = dqdxU(i+1,j,k)
686     c ENDIF
687     c IF ( iG .EQ. sNx-1 ) THEN
688     c dqdxT(i,j,k) = dqdxU(i ,j,k)
689     c ENDIF
690     C--
691     ENDDO
692     ENDDO
693     ENDDO
694     C-- ::V-points
695     DO k=1,Nr
696     DO j=jMin,jMax
697     DO i=iMin,iMax
698     jG = myYGlobalLo-1+(bj-1)*sNy+J
699     dqdxV(i,j,k) = 0.
700     IF ( vmask(i,j,k) .ne. 0. ) THEN
701     dqdxV(i,j,k) = ( dqdxT(i,j,k)*pmask(i,j,k) + dqdxT(i,j-1,k)*pmask(i,j-1,k) )/( pmask(i,j,k) + pmask(i,j-1,k) )
702     ENDIF
703     C
704     C-- bondaries:
705     c IF ( jG .EQ. 1 ) THEN
706     c dqdxV(i,j,k) = dqdxT(i, j+1,k)
707     c ENDIF
708     c IF ( jG .EQ. sNy-1 ) THEN
709     c dqdxV(i,j,k) = dqdxT(i,j,k)
710     c ENDIF
711     C--
712     ENDDO
713     ENDDO
714     ENDDO
715     C
716     ENDIF
717     C--
718     C ===============
719     C-- Evaluate Kpv :
720     C ===============
721     C--
722     IF ( N2local .AND. betaterm ) THEN
723     c IF ( N2local ) THEN
724     c
725     DO i=iMin,iMax
726     DO j=jMin,jMax
727     kbot = int(kbotindex(i,j))
728     c zonal component
729     int1 = 0.
730     int2 = 0.
731     int3 = 0.
732     ku = int( max(mldindex(i-1,j),mldindex(i,j)) )
733     DO k=1,ku
734     int1 = int1 + delz(k)*dqdyU(i,j,k)
735     ENDDO
736     IF ( ku .LT. kbot ) int3 = delz(kbot)*dqdyU(i,j,kbot)
737     IF ( ku .LT. kbot-1 ) THEN
738     DO k=ku+1,kbot-1
739     int2 = int2 + delz(k)*dqdyU(i,j,k)
740     ENDDO
741     ENDIF
742     qymld(i,j) = int1
743     qyint(i,j) = int2
744     qybot(i,j) = int3
745     c meridional component
746     int4 = 0.
747     int5 = 0.
748     int6 = 0.
749     kv = int( max(mldindex(i,j-1),mldindex(i,j)) )
750     DO k=1,kv
751     int4 = int4 + delz(k)*dqdxV(i,j,k)
752     ENDDO
753     IF ( kv .LT. kbot ) int6 = delz(kbot)*dqdxV(i,j,kbot)
754     IF ( kv .LT. kbot-1 ) THEN
755     DO k=kv+1,kbot-1
756     int5 = int5 + delz(k)*dqdxV(i,j,k)
757     ENDDO
758     ENDIF
759     qxmld(i,j) = int4
760     qxint(i,j) = int5
761     qxbot(i,j) = int6
762     c compute the K
763     num_i(i,j) = int3*int4 - int1*int6
764     num_b(i,j) = int1*int5 - int2*int4
765     denom(i,j) = int2*int6 - int3*int5
766     IF ( denom(i,j) .ne. 0. ) THEN
767     ratioqy(i,j) = num_i(i,j)/denom(i,j)
768     ratio_b(i,j) = num_b(i,j)/denom(i,j)
769     ENDIF
770     IF ( abs(ratioqy(i,j)-25.) .ge. 25. .OR.
771     & abs(ratio_b(i,j)-25.) .ge. 25. ) THEN
772     fac = 0.
773     ELSE
774     fac = 1.
775     ENDIF
776     DO k=1,kbot-1
777     IF ( k .le. ku ) THEN
778     KpvU(i,j,k) = fac * Kpvref * umask(i,j,k)
779     ELSE
780     KpvU(i,j,k) = fac * ratioqy(i,j) * Kpvref * umask(i,j,k)
781     ENDIF
782     IF ( k .le. kv ) THEN
783     KpvV(i,j,k) = fac * Kpvref * vmask(i,j,k)
784     ELSE
785     KpvV(i,j,k) = fac * ratioqy(i,j) * Kpvref * vmask(i,j,k)
786     ENDIF
787     ENDDO
788     KpvU(i,j,kbot) = fac * ratio_b(i,j) * Kpvref * umask(i,j,kbot)
789     KpvV(i,j,kbot) = fac * ratio_b(i,j) * Kpvref * vmask(i,j,kbot)
790     ENDDO
791     ENDDO
792     c
793     ELSE
794     c
795     DO i=iMin,iMax
796     DO j=jMin,jMax
797     pvFacT(i,j) = 0. _d 0
798     ENDDO
799     ENDDO
800     C =========================
801     C-- Lateral variation of Kpv:
802     DO i=iMin,iMax
803     DO j=jMin,jMax
804     crmw pvFacT(i,j) = qp2(i,j,bi,bj) * 1.0 _d 0
805     pvFacT(i,j) = 1.0 _d 0
806     ENDDO
807     ENDDO
808     C--
809     C ::At U-point
810     DO i=iMin,iMax
811     DO j=jMin,jMax
812     iG = myXGlobalLo-1+(bi-1)*sNx+I
813     C--
814     pvFacU(i,j) = 0.5 * ( pvFacT(i,j) + pvFacT(i+1,j) )
815     C--
816     C-- bondaries:
817     IF ( iG .EQ. 1 ) THEN
818     pvFacU(i,j) = pvFacT(i+1,j)
819     ENDIF
820     C--
821     ENDDO
822     ENDDO
823     C--
824     C ::At V-point
825     DO i=iMin,iMax
826     DO j=jMin,jMax
827     jG = myYGlobalLo-1+(bj-1)*sNy+J
828     C--
829     pvFacV(i,j) = 0.5 * ( pvFacT(i,j) + pvFacT(i,j+1) )
830     C--
831     C-- bondaries:
832     IF ( jG .EQ. 1 ) THEN
833     pvFacV(i,j) = pvFacT(i,j+1)
834     ENDIF
835     C--
836     ENDDO
837     ENDDO
838     C--
839     C =========================
840     DO K=1,Nr
841     DO i=iMin,iMax
842     DO j=jMin,jMax
843     KpvU(i,j,K) = Kpvref * pvFacU(i,j)
844     KpvV(i,j,K) = Kpvref * pvFacV(i,j)
845     c IF ( I .EQ. 10 .AND. J .EQ. 10) THEN
846     c write(0,*) 't0', Kpvref, pvFacU(i,j), KpvU(i,j,K)
847     c ENDIF
848     ENDDO
849     ENDDO
850     ENDDO
851     c
852     ENDIF
853     C--
854     C =========================================
855     C--
856     C =======================
857     C-- Evaluate gUpv and gVpv:
858     C =======================
859     C--
860     c write(0,*)' Evaluate gUpv and gVpv'
861     C
862     DO k=1,Nr
863     DO j=jMin,jMax
864     DO i=iMin,iMax
865     gUpv(i,j,k,bi,bj) = - KpvU(i,j,K) * dqdyU(i,j,k)
866     c gUpv(i,j,k,bi,bj) = - Kpvref * ( dqdyU(i,j,k) - dfdy(i,j,bi,bj) )
867     gVpv(i,j,k,bi,bj) = KpvV(i,j,k) * dqdxV(i,j,k)
868     c gVpv(i,j,k,bi,bj) = Kpvref * dqdxV(i,j,k)
869     ENDDO
870     ENDDO
871     ENDDO
872     C
873     C--
874     ELSE !!!!!!!(pvForcing = .FALSE.)
875     C--
876     IF ( GMadv ) THEN
877     C
878     C Compute horizontal eddy-induced velocities and save them in
879     C arrays gUpv, gVpv
880     DO k=1,Nr
881     DO j=jMin,jMax
882     DO i=iMin,iMax
883     arr(i,j,k) = - Kpvref * 0.5 *
884     & ( K13(i-1,j,k)+K13(i,j,k) ) * umask(i,j,k)
885     arrprime(i,j,k) = - Kpvref * 0.5 *
886     & ( K23(i,j-1,k)+K23(i,j,k) ) * vmask(i,j,k)
887     ENDDO
888     ENDDO
889     ENDDO
890     DO j=jMin,jMax
891     DO i=iMin,iMax
892     arr(i,j,Nr+1) = 0.
893     arrprime(i,j,Nr+1) = 0.
894     ENDDO
895     ENDDO
896     C Take z-derivative of arr and arrprime
897     DO k=1,Nr
898     DO j=jMin,jMax
899     DO i=iMin,iMax
900     KpvU(i,j,k) = 0. _d 0
901     KpvV(i,j,k) = 0. _d 0
902     gUpv(i,j,k,bi,bj) = recip_drF(k) * (arr(i,j,k)-arr(i,j,k+1))
903     gVpv(i,j,k,bi,bj) = recip_drF(k) * (arrprime(i,j,k)-arrprime(i,j,k+1))
904     ENDDO
905     ENDDO
906     ENDDO
907     C
908     ELSE
909     C
910     DO k=1,Nr
911     DO j=jMin,jMax
912     DO i=iMin,iMax
913     KpvU(i,j,k) = 0. _d 0
914     KpvV(i,j,k) = 0. _d 0
915     gUpv(i,j,k,bi,bj) = 0. _d 0
916     gVpv(i,j,k,bi,bj) = 0. _d 0
917     ENDDO
918     ENDDO
919     ENDDO
920     C
921     ENDIF
922     C--
923     ENDIF
924    
925     ctest
926     _BARRIER
927     _BEGIN_MASTER( myThid )
928    
929     WRITE(suff,'(I10.10)') myIter
930    
931     C-- Write model fields
932     CALL WRITE_FLD_XY_RL('kbot.',suff,kbotindex,myIter,myThid)
933     c CALL WRITE_FLD_XY_RL('mldindex.',suff,mldindex,myIter,myThid)
934     c CALL WRITE_FLD_XY_RL('mld.',suff,mld,myIter,myThid)
935     c CALL WRITE_FLD_XY_RL('qymld.',suff,qymld,myIter,myThid)
936     c CALL WRITE_FLD_XY_RL('qyint.',suff,qyint,myIter,myThid)
937     c CALL WRITE_FLD_XY_RL('qybot.',suff,qybot,myIter,myThid)
938     c CALL WRITE_FLD_XY_RL('qxmld.',suff,qxmld,myIter,myThid)
939     c CALL WRITE_FLD_XY_RL('qxint.',suff,qxint,myIter,myThid)
940     c CALL WRITE_FLD_XY_RL('qxbot.',suff,qxbot,myIter,myThid)
941     c CALL WRITE_FLD_XY_RL('ratioqy.',suff,ratioqy,myIter,myThid)
942     c CALL WRITE_FLD_XY_RL('ratio_b.',suff,ratio_b,myIter,myThid)
943     c CALL WRITE_FLD_XYZ_RL( 'dqdyU.',suff,dqdyU,myIter,myThid)
944     c CALL WRITE_FLD_XYZ_RL( 'dqdxV.',suff,dqdxV,myIter,myThid)
945     c CALL WRITE_FLD_XYZ_RL( 'KpvU.',suff,KpvU,myIter,myThid)
946     c CALL WRITE_FLD_XYZ_RL( 'KpvV.',suff,KpvV,myIter,myThid)
947     c CALL WRITE_FLD_XYZ_RL( 'gUpv.',suff,gUpv,myIter,myThid)
948     c CALL WRITE_FLD_XYZ_RL( 'gVpv.',suff,gVpv,myIter,myThid)
949     c CALL WRITE_FLD_XY_RL('num_i.',suff,num_i,myIter,myThid)
950     c CALL WRITE_FLD_XY_RL('num_b.',suff,num_b,myIter,myThid)
951     c CALL WRITE_FLD_XY_RL('denom.',suff,denom,myIter,myThid)
952    
953     _END_MASTER( myThid )
954     _BARRIER
955     ctest
956     C--
957     C =========================================================
958     c DO k=1,Nr
959     c DO j=jMin,jMax
960     c DO i=iMin,iMax
961     c write(0,*) 'TTTT ', Kpvref, pvFacU(i,j)
962     c ENDDO
963     c ENDDO
964     c ENDDO
965     C--
966     C =========================================================
967    
968     RETURN
969     END

  ViewVC Help
Powered by ViewVC 1.1.22