/[MITgcm]/MITgcm/pkg/longstep/longstep_thermodynamics.F
ViewVC logotype

Contents of /MITgcm/pkg/longstep/longstep_thermodynamics.F

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


Revision 1.1 - (show annotations) (download)
Fri Jun 26 23:10:10 2009 UTC (14 years, 11 months ago) by jahn
Branch: MAIN
CVS Tags: checkpoint62, checkpoint61v, checkpoint61w, checkpoint61t, checkpoint61u, checkpoint61s, checkpoint61z, checkpoint61x, checkpoint61y
add package longstep

1 C $Header$
2 C $Name$
3
4 #include "LONGSTEP_OPTIONS.h"
5
6 #ifdef ALLOW_AUTODIFF_TAMC
7 # ifdef ALLOW_GMREDI
8 # include "GMREDI_OPTIONS.h"
9 # endif
10 #endif /* ALLOW_AUTODIFF_TAMC */
11
12 CBOP
13 C !ROUTINE: LONGSTEP_THERMODYNAMICS
14 C !INTERFACE:
15 SUBROUTINE LONGSTEP_THERMODYNAMICS(myTime, myIter, myThid)
16 C !DESCRIPTION: \bv
17 C *==========================================================*
18 C | SUBROUTINE LONGSTEP_THERMODYNAMICS
19 C | o Controlling routine for the prognostics of passive tracers
20 C | with longer time step.
21 C *===========================================================
22 C | This is a copy of THERMODYNAMICS, but only with the
23 C | parts relevant to ptracers, and dynamics fields replaced
24 C | by their longstep averages.
25 C | When THERMODYNAMICS is changed, this routine probably has
26 C | to be changed too :(
27 C *==========================================================*
28 C \ev
29
30 C !USES:
31 IMPLICIT NONE
32 C == Global variables ===
33 #include "SIZE.h"
34 #include "EEPARAMS.h"
35 #include "PARAMS.h"
36 #include "RESTART.h"
37 #include "DYNVARS.h"
38 #include "GRID.h"
39 #ifdef ALLOW_GENERIC_ADVDIFF
40 # include "GAD.h"
41 #endif
42 #include "LONGSTEP_PARAMS.h"
43 #include "LONGSTEP.h"
44 #ifdef ALLOW_PTRACERS
45 # include "PTRACERS_SIZE.h"
46 # include "PTRACERS_PARAMS.h"
47 # include "PTRACERS_FIELDS.h"
48 #endif
49 #ifdef ALLOW_TIMEAVE
50 # include "TIMEAVE_STATV.h"
51 #endif
52
53 #ifdef ALLOW_AUTODIFF_TAMC
54 # include "tamc.h"
55 # include "tamc_keys.h"
56 # include "FFIELDS.h"
57 # include "SURFACE.h"
58 # include "EOS.h"
59 # ifdef ALLOW_GMREDI
60 # include "GMREDI.h"
61 # endif
62 # ifdef ALLOW_EBM
63 # include "EBM.h"
64 # endif
65 #endif /* ALLOW_AUTODIFF_TAMC */
66
67 C !INPUT/OUTPUT PARAMETERS:
68 C == Routine arguments ==
69 C myTime - Current time in simulation
70 C myIter - Current iteration number in simulation
71 C myThid - Thread number for this instance of the routine.
72 _RL myTime
73 INTEGER myIter
74 INTEGER myThid
75
76 #ifdef ALLOW_LONGSTEP
77 C !LOCAL VARIABLES:
78 C == Local variables
79 C xA, yA - Per block temporaries holding face areas
80 C uFld, vFld, wFld - Local copy of velocity field (3 components)
81 C uTrans, vTrans, rTrans - Per block temporaries holding flow transport
82 C o uTrans: Zonal transport
83 C o vTrans: Meridional transport
84 C o rTrans: Vertical transport
85 C rTransKp1 o vertical volume transp. at interface k+1
86 C maskUp o maskUp: land/water mask for W points
87 C fVer[STUV] o fVer: Vertical flux term - note fVer
88 C is "pipelined" in the vertical
89 C so we need an fVer for each
90 C variable.
91 C kappaRT, - Total diffusion in vertical at level k, for T and S
92 C kappaRS (background + spatially varying, isopycnal term).
93 C kappaRTr - Total diffusion in vertical at level k,
94 C for each passive Tracer
95 C kappaRk - Total diffusion in vertical, all levels, 1 tracer
96 C useVariableK = T when vertical diffusion is not constant
97 C iMin, iMax - Ranges and sub-block indices on which calculations
98 C jMin, jMax are applied.
99 C bi, bj
100 C k, kup, - Index for layer above and below. kup and kDown
101 C kDown, km1 are switched with layer to be the appropriate
102 C index into fVerTerm.
103 _RS xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
104 _RS yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
105 _RL uFld (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
106 _RL vFld (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
107 _RL wFld (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
108 _RL uTrans (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
109 _RL vTrans (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
110 _RL rTrans (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
111 _RL rTransKp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
112 _RS maskUp (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
113 _RL fVerT (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
114 _RL fVerS (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
115 _RL kappaRT (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
116 _RL kappaRS (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
117 #ifdef ALLOW_PTRACERS
118 _RL fVerP (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2,PTRACERS_num)
119 _RL kappaRTr(1-Olx:sNx+Olx,1-Oly:sNy+Oly,PTRACERS_num)
120 #endif
121 _RL kappaRk (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
122 INTEGER iMin, iMax
123 INTEGER jMin, jMax
124 INTEGER bi, bj
125 INTEGER i, j
126 INTEGER k, km1, kup, kDown
127 #ifdef ALLOW_ADAMSBASHFORTH_3
128 INTEGER iterNb, m1, m2
129 #endif
130 #ifdef ALLOW_TIMEAVE
131 LOGICAL useVariableK
132 #endif
133 #ifdef ALLOW_PTRACERS
134 INTEGER iTracer, ip
135 #endif
136 INTEGER trIter
137 CEOP
138
139 #ifdef ALLOW_DEBUG
140 IF ( debugLevel .GE. debLevB )
141 & CALL DEBUG_ENTER('LONGSTEP_THERMODYNAMICS',myThid)
142 #endif
143
144 C time for a ptracer time step?
145 IF ( LS_doTimeStep ) THEN
146
147 #ifdef ALLOW_AUTODIFF_TAMC
148 C-- dummy statement to end declaration part
149 ikey = 1
150 itdkey = 1
151 #endif /* ALLOW_AUTODIFF_TAMC */
152
153 #ifdef ALLOW_AUTODIFF_TAMC
154 C-- HPF directive to help TAMC
155 CHPF$ INDEPENDENT
156 #endif /* ALLOW_AUTODIFF_TAMC */
157
158 DO bj=myByLo(myThid),myByHi(myThid)
159
160 #ifdef ALLOW_AUTODIFF_TAMC
161 C-- HPF directive to help TAMC
162 CHPF$ INDEPENDENT, NEW (rTrans,fVerT,fVerS
163 CHPF$& ,utrans,vtrans,xA,yA
164 CHPF$& ,kappaRT,kappaRS
165 CHPF$& )
166 # ifdef ALLOW_PTRACERS
167 CHPF$ INDEPENDENT, NEW (fVerP,kappaRTr)
168 # endif
169 #endif /* ALLOW_AUTODIFF_TAMC */
170
171 DO bi=myBxLo(myThid),myBxHi(myThid)
172
173 #ifdef ALLOW_AUTODIFF_TAMC
174 act1 = bi - myBxLo(myThid)
175 max1 = myBxHi(myThid) - myBxLo(myThid) + 1
176 act2 = bj - myByLo(myThid)
177 max2 = myByHi(myThid) - myByLo(myThid) + 1
178 act3 = myThid - 1
179 max3 = nTx*nTy
180 act4 = ikey_dynamics - 1
181 itdkey = (act1 + 1) + act2*max1
182 & + act3*max1*max2
183 & + act4*max1*max2*max3
184 #endif /* ALLOW_AUTODIFF_TAMC */
185
186 C-- Set up work arrays with valid (i.e. not NaN) values
187 C These inital values do not alter the numerical results. They
188 C just ensure that all memory references are to valid floating
189 C point numbers. This prevents spurious hardware signals due to
190 C uninitialised but inert locations.
191
192 DO j=1-OLy,sNy+OLy
193 DO i=1-OLx,sNx+OLx
194 xA(i,j) = 0. _d 0
195 yA(i,j) = 0. _d 0
196 uTrans(i,j) = 0. _d 0
197 vTrans(i,j) = 0. _d 0
198 rTrans (i,j) = 0. _d 0
199 rTransKp1(i,j) = 0. _d 0
200 ENDDO
201 ENDDO
202
203 DO k=1,Nr
204 DO j=1-OLy,sNy+OLy
205 DO i=1-OLx,sNx+OLx
206 C This is currently also used by IVDC and Diagnostics
207 kappaRk(i,j,k) = 0. _d 0
208 ENDDO
209 ENDDO
210 ENDDO
211
212 #ifdef ALLOW_PTRACERS
213 IF ( usePTRACERS ) THEN
214 DO ip=1,PTRACERS_num
215 DO j=1-OLy,sNy+OLy
216 DO i=1-OLx,sNx+OLx
217 fVerP (i,j,1,ip) = 0. _d 0
218 fVerP (i,j,2,ip) = 0. _d 0
219 kappaRTr(i,j,ip) = 0. _d 0
220 ENDDO
221 ENDDO
222 ENDDO
223 C- set tracer tendency to zero:
224 DO iTracer=1,PTRACERS_num
225 DO k=1,Nr
226 DO j=1-OLy,sNy+OLy
227 DO i=1-OLx,sNx+OLx
228 gPTr(i,j,k,bi,bj,itracer) = 0. _d 0
229 ENDDO
230 ENDDO
231 ENDDO
232 ENDDO
233 ENDIF
234 #endif
235
236 c iMin = 1-OLx
237 c iMax = sNx+OLx
238 c jMin = 1-OLy
239 c jMax = sNy+OLy
240
241 #ifdef ALLOW_AUTODIFF_TAMC
242 cph avoids recomputation of integrate_for_w
243 CADJ STORE wvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
244 #endif /* ALLOW_AUTODIFF_TAMC */
245
246 C-- Attention: by defining "SINGLE_LAYER_MODE" in CPP_OPTIONS.h
247 C-- MOST of THERMODYNAMICS will be disabled
248 #ifndef SINGLE_LAYER_MODE
249
250 #ifdef ALLOW_AUTODIFF_TAMC
251 CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
252 CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
253 CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
254 CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
255 # if ((defined ALLOW_DEPTH_CONTROL) || (defined NONLIN_FRSURF))
256 CADJ STORE gtnm1(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
257 CADJ STORE gsnm1(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
258 # endif
259 #endif /* ALLOW_AUTODIFF_TAMC */
260
261 #ifndef DISABLE_MULTIDIM_ADVECTION
262 C-- Some advection schemes are better calculated using a multi-dimensional
263 C method in the absence of any other terms and, if used, is done here.
264 C
265 C The CPP flag DISABLE_MULTIDIM_ADVECTION is currently unset in GAD_OPTIONS.h
266 C The default is to use multi-dimensinal advection for non-linear advection
267 C schemes. However, for the sake of efficiency of the adjoint it is necessary
268 C to be able to exclude this scheme to avoid excessive storage and
269 C recomputation. It *is* differentiable, if you need it.
270 C Edit GAD_OPTIONS.h and #define DISABLE_MULTIDIM_ADVECTION to
271 C disable this section of code.
272
273 C Since passive tracers are configurable separately from T,S we
274 C call the multi-dimensional method for PTRACERS regardless
275 C of whether multiDimAdvection is set or not.
276 #ifdef ALLOW_PTRACERS
277 IF ( usePTRACERS ) THEN
278 #ifdef ALLOW_DEBUG
279 IF ( debugLevel .GE. debLevB )
280 & CALL DEBUG_CALL('PTRACERS_ADVECTION',myThid)
281 #endif
282 CALL PTRACERS_ADVECTION(
283 I bi,bj,myIter,myTime,myThid )
284 ENDIF
285 #endif /* ALLOW_PTRACERS */
286 #endif /* DISABLE_MULTIDIM_ADVECTION */
287
288 #ifdef ALLOW_DEBUG
289 IF ( debugLevel .GE. debLevB )
290 & CALL DEBUG_MSG('ENTERING DOWNWARD K LOOP IN LONGSTEP',
291 & myThid)
292 #endif
293
294 C-- Start of thermodynamics loop
295 DO k=Nr,1,-1
296 #ifdef ALLOW_AUTODIFF_TAMC
297 C? Patrick Is this formula correct?
298 cph Yes, but I rewrote it.
299 cph Also, the kappaR? need the index and subscript k!
300 kkey = (itdkey-1)*Nr + k
301 #endif /* ALLOW_AUTODIFF_TAMC */
302
303 C-- km1 Points to level above k (=k-1)
304 C-- kup Cycles through 1,2 to point to layer above
305 C-- kDown Cycles through 2,1 to point to current layer
306
307 km1 = MAX(1,k-1)
308 kup = 1+MOD(k+1,2)
309 kDown= 1+MOD(k,2)
310
311 iMin = 1-OLx
312 iMax = sNx+OLx
313 jMin = 1-OLy
314 jMax = sNy+OLy
315
316 IF (k.EQ.Nr) THEN
317 DO j=1-Oly,sNy+Oly
318 DO i=1-Olx,sNx+Olx
319 rTransKp1(i,j) = 0. _d 0
320 ENDDO
321 ENDDO
322 ELSE
323 DO j=1-Oly,sNy+Oly
324 DO i=1-Olx,sNx+Olx
325 rTransKp1(i,j) = rTrans(i,j)
326 ENDDO
327 ENDDO
328 ENDIF
329 #ifdef ALLOW_AUTODIFF_TAMC
330 CADJ STORE rTransKp1(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
331 #endif
332
333 C-- Get temporary terms used by tendency routines :
334 C- Calculate horizontal "volume transport" through tracer cell face
335 C anelastic: uTrans,vTrans are scaled by rhoFacC (~ mass transport)
336 CALL CALC_COMMON_FACTORS (
337 I LS_uVel, LS_vVel,
338 O uFld, vFld, uTrans, vTrans, xA, yA,
339 I k,bi,bj, myThid )
340
341 C- Calculate vertical "volume transport" through tracer cell face
342 IF (k.EQ.1) THEN
343 C- Surface interface :
344 DO j=1-Oly,sNy+Oly
345 DO i=1-Olx,sNx+Olx
346 wFld(i,j) = 0. _d 0
347 maskUp(i,j) = 0. _d 0
348 rTrans(i,j) = 0. _d 0
349 ENDDO
350 ENDDO
351 ELSE
352 C- Interior interface :
353 C anelastic: rTrans is scaled by rhoFacF (~ mass transport)
354 DO j=1-Oly,sNy+Oly
355 DO i=1-Olx,sNx+Olx
356 wFld(i,j) = LS_wVel(i,j,k,bi,bj)
357 maskUp(i,j) = maskC(i,j,k-1,bi,bj)*maskC(i,j,k,bi,bj)
358 rTrans(i,j) = wFld(i,j)*rA(i,j,bi,bj)*maskUp(i,j)
359 & *deepFac2F(k)*rhoFacF(k)
360 ENDDO
361 ENDDO
362 ENDIF
363
364 #ifdef ALLOW_GMREDI
365 C-- Residual transp = Bolus transp + Eulerian transp
366 IF (useGMRedi) THEN
367 CALL GMREDI_CALC_UVFLOW(
368 U uFld, vFld, uTrans, vTrans,
369 I k, bi, bj, myThid )
370 IF (K.GE.2) THEN
371 CALL GMREDI_CALC_WFLOW(
372 U wFld, rTrans,
373 I k, bi, bj, myThid )
374 ENDIF
375 ENDIF
376 # ifdef ALLOW_AUTODIFF_TAMC
377 CADJ STORE rTrans(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
378 CADJ STORE wfld(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
379 # ifdef GM_BOLUS_ADVEC
380 CADJ STORE ufld(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
381 CADJ STORE vfld(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
382 CADJ STORE uTrans(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
383 CADJ STORE vTrans(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
384 # endif
385 # endif /* ALLOW_AUTODIFF_TAMC */
386 #endif /* ALLOW_GMREDI */
387
388 iMin = 1-OLx+2
389 iMax = sNx+OLx-1
390 jMin = 1-OLy+2
391 jMax = sNy+OLy-1
392
393 C-- Calculate active tracer tendencies (gT,gS,...)
394 C and step forward storing result in gT, gS, etc.
395 C--
396 # ifdef ALLOW_AUTODIFF_TAMC
397 # if ((defined NONLIN_FRSURF) || (defined ALLOW_DEPTH_CONTROL)) && (defined ALLOW_GMREDI)
398 # ifdef GM_NON_UNITY_DIAGONAL
399 CADJ STORE kux(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
400 CADJ STORE kvy(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
401 # endif
402 # ifdef GM_EXTRA_DIAGONAL
403 CADJ STORE kuz(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
404 CADJ STORE kvz(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
405 # endif
406 # endif
407 # endif /* ALLOW_AUTODIFF_TAMC */
408 C
409 #ifdef ALLOW_AUTODIFF_TAMC
410 # if (defined NONLIN_FRSURF) || (defined ALLOW_DEPTH_CONTROL)
411 cph-test
412 CADJ STORE uFld(:,:), vFld(:,:), wFld(:,:)
413 CADJ & = comlev1_bibj_k, key=kkey, byte=isbyte
414 CADJ STORE uTrans(:,:), vTrans(:,:)
415 CADJ & = comlev1_bibj_k, key=kkey, byte=isbyte
416 CADJ STORE xA(:,:), yA(:,:)
417 CADJ & = comlev1_bibj_k, key=kkey, byte=isbyte
418 # endif
419 #endif /* ALLOW_AUTODIFF_TAMC */
420 C
421 #ifdef ALLOW_PTRACERS
422 IF ( usePTRACERS ) THEN
423 IF ( .NOT.implicitDiffusion ) THEN
424 CALL PTRACERS_CALC_DIFF(
425 I bi,bj,iMin,iMax,jMin,jMax,k,
426 I maskUp,
427 O kappaRTr,
428 I myThid)
429 ENDIF
430 # ifdef ALLOW_AUTODIFF_TAMC
431 CADJ STORE kappaRTr(:,:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
432 # endif /* ALLOW_AUTODIFF_TAMC */
433 CALL PTRACERS_INTEGRATE(
434 I bi,bj,k,
435 I xA, yA, maskUp, uFld, vFld, wFld,
436 I uTrans, vTrans, rTrans, rTransKp1,
437 I kappaRTr,
438 U fVerP,
439 I myTime,myIter,myThid)
440 ENDIF
441 #endif /* ALLOW_PTRACERS */
442
443 C-- end of thermodynamic k loop (Nr:1)
444 ENDDO
445
446 #ifdef ALLOW_DOWN_SLOPE
447 #ifdef ALLOW_PTRACERS
448 IF ( usePTRACERS .AND. useDOWN_SLOPE ) THEN
449 CALL PTRACERS_DWNSLP_APPLY(
450 I bi, bj, myTime, myIter, myThid )
451 ENDIF
452 #endif /* ALLOW_PTRACERS */
453 #endif /* ALLOW_DOWN_SLOPE */
454
455 C All explicit advection/diffusion/sources should now be
456 C done. The updated tracer field is in gPtr. Accumalate
457 C explicit tendency and also reset gPtr to initial tracer
458 C field for implicit matrix calculation
459
460 #ifdef ALLOW_MATRIX
461 IF (useMATRIX)
462 & CALL MATRIX_STORE_TENDENCY_EXP(bi,bj, myTime,myIter,myThid)
463 #endif
464
465 iMin = 1
466 iMax = sNx
467 jMin = 1
468 jMax = sNy
469
470 #ifdef ALLOW_PTRACERS
471 IF ( usePTRACERS ) THEN
472 C-- Vertical advection/diffusion (implicit) for passive tracers
473 CALL PTRACERS_IMPLICIT(
474 U kappaRk,
475 I bi, bj, myTime, myIter, myThid )
476 ENDIF
477 #endif /* ALLOW_PTRACERS */
478
479 #endif /* SINGLE_LAYER_MODE */
480
481 C-- end bi,bj loops.
482 ENDDO
483 ENDDO
484
485 #ifdef ALLOW_DEBUG
486 If (debugMode) THEN
487 CALL DEBUG_STATS_RL(Nr,LS_uVel,'LS_Uvel (THERMODYNAMICS)',myThid)
488 CALL DEBUG_STATS_RL(Nr,LS_vVel,'LS_Vvel (THERMODYNAMICS)',myThid)
489 CALL DEBUG_STATS_RL(Nr,LS_wVel,'LS_Wvel (THERMODYNAMICS)',myThid)
490 CALL DEBUG_STATS_RL(Nr,LS_theta,'LS_Theta (THERMODYNAMICS)',
491 & myThid)
492 CALL DEBUG_STATS_RL(Nr,LS_salt,'LS_Salt (THERMODYNAMICS)',myThid)
493 #ifdef ALLOW_PTRACERS
494 IF ( usePTRACERS ) THEN
495 CALL PTRACERS_DEBUG(myThid)
496 ENDIF
497 #endif /* ALLOW_PTRACERS */
498 ENDIF
499 #endif /* ALLOW_DEBUG */
500
501 C LS_doTimeStep
502 ENDIF
503
504 #ifdef ALLOW_DEBUG
505 IF ( debugLevel .GE. debLevB )
506 & CALL DEBUG_LEAVE('LONGSTEP_THERMODYNAMICS',myThid)
507 #endif
508
509 #endif /* ALLOW_LONGSTEP */
510
511 RETURN
512 END

  ViewVC Help
Powered by ViewVC 1.1.22