/[MITgcm]/MITgcm/model/src/do_oceanic_phys.F
ViewVC logotype

Contents of /MITgcm/model/src/do_oceanic_phys.F

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


Revision 1.120 - (show annotations) (download)
Fri Nov 9 22:42:00 2012 UTC (11 years, 6 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint64a, checkpoint64c, checkpoint64b
Changes since 1.119: +3 -20 lines
remove reset of addMass (now done in load_fields_driver.F)

1 C $Header: /u/gcmpack/MITgcm/model/src/do_oceanic_phys.F,v 1.119 2012/08/06 14:22:53 gforget Exp $
2 C $Name: $
3
4 #include "PACKAGES_CONFIG.h"
5 #include "CPP_OPTIONS.h"
6
7 #ifdef ALLOW_AUTODIFF_TAMC
8 # ifdef ALLOW_GMREDI
9 # include "GMREDI_OPTIONS.h"
10 # endif
11 # ifdef ALLOW_KPP
12 # include "KPP_OPTIONS.h"
13 # endif
14 # ifdef ALLOW_SEAICE
15 # include "SEAICE_OPTIONS.h"
16 # endif
17 # ifdef ALLOW_EXF
18 # include "EXF_OPTIONS.h"
19 # endif
20 #endif /* ALLOW_AUTODIFF_TAMC */
21
22 CBOP
23 C !ROUTINE: DO_OCEANIC_PHYS
24 C !INTERFACE:
25 SUBROUTINE DO_OCEANIC_PHYS(myTime, myIter, myThid)
26 C !DESCRIPTION: \bv
27 C *==========================================================*
28 C | SUBROUTINE DO_OCEANIC_PHYS
29 C | o Controlling routine for oceanic physics and
30 C | parameterization
31 C *==========================================================*
32 C | o originally, part of S/R thermodynamics
33 C *==========================================================*
34 C \ev
35
36 C !USES:
37 IMPLICIT NONE
38 C == Global variables ===
39 #include "SIZE.h"
40 #include "EEPARAMS.h"
41 #include "PARAMS.h"
42 #include "GRID.h"
43 #include "DYNVARS.h"
44 #ifdef ALLOW_TIMEAVE
45 #include "TIMEAVE_STATV.h"
46 #endif
47 #ifdef ALLOW_BALANCE_FLUXES
48 #include "FFIELDS.h"
49 #endif
50
51 #ifdef ALLOW_AUTODIFF_TAMC
52 # include "AUTODIFF_MYFIELDS.h"
53 # include "tamc.h"
54 # include "tamc_keys.h"
55 #ifndef ALLOW_BALANCE_FLUXES
56 # include "FFIELDS.h"
57 #endif
58 # include "SURFACE.h"
59 # include "EOS.h"
60 # ifdef ALLOW_KPP
61 # include "KPP.h"
62 # endif
63 # ifdef ALLOW_GGL90
64 # include "GGL90.h"
65 # endif
66 # ifdef ALLOW_GMREDI
67 # include "GMREDI.h"
68 # endif
69 # ifdef ALLOW_EBM
70 # include "EBM.h"
71 # endif
72 # ifdef ALLOW_EXF
73 # include "ctrl.h"
74 # include "EXF_FIELDS.h"
75 # ifdef ALLOW_BULKFORMULAE
76 # include "EXF_CONSTANTS.h"
77 # endif
78 # endif
79 # ifdef ALLOW_SEAICE
80 # include "SEAICE_SIZE.h"
81 # include "SEAICE.h"
82 # endif
83 # ifdef ALLOW_THSICE
84 # include "THSICE_VARS.h"
85 # endif
86 # ifdef ALLOW_SALT_PLUME
87 # include "SALT_PLUME.h"
88 # endif
89 #endif /* ALLOW_AUTODIFF_TAMC */
90
91 C !INPUT/OUTPUT PARAMETERS:
92 C == Routine arguments ==
93 C myTime :: Current time in simulation
94 C myIter :: Current iteration number in simulation
95 C myThid :: Thread number for this instance of the routine.
96 _RL myTime
97 INTEGER myIter
98 INTEGER myThid
99
100 C !LOCAL VARIABLES:
101 C == Local variables
102 C rhoK, rhoKm1 :: Density at current level, and level above
103 C iMin, iMax :: Ranges and sub-block indices on which calculations
104 C jMin, jMax are applied.
105 C bi, bj :: tile indices
106 C i,j,k :: loop indices
107 _RL rhoKp1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
108 _RL rhoKm1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
109 _RL sigmaX (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
110 _RL sigmaY (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
111 _RL sigmaR (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
112 INTEGER iMin, iMax
113 INTEGER jMin, jMax
114 INTEGER bi, bj
115 INTEGER i, j, k
116 INTEGER doDiagsRho
117 #ifdef ALLOW_DIAGNOSTICS
118 LOGICAL DIAGNOSTICS_IS_ON
119 EXTERNAL DIAGNOSTICS_IS_ON
120 #endif /* ALLOW_DIAGNOSTICS */
121
122 CEOP
123
124 #ifdef ALLOW_AUTODIFF_TAMC
125 C-- dummy statement to end declaration part
126 itdkey = 1
127 #endif /* ALLOW_AUTODIFF_TAMC */
128
129 #ifdef ALLOW_DEBUG
130 IF (debugMode) CALL DEBUG_ENTER('DO_OCEANIC_PHYS',myThid)
131 #endif
132
133 doDiagsRho = 0
134 #ifdef ALLOW_DIAGNOSTICS
135 IF ( useDiagnostics .AND. fluidIsWater ) THEN
136 IF ( DIAGNOSTICS_IS_ON('MXLDEPTH',myThid) )
137 & doDiagsRho = doDiagsRho + 1
138 IF ( DIAGNOSTICS_IS_ON('DRHODR ',myThid) )
139 & doDiagsRho = doDiagsRho + 2
140 IF ( DIAGNOSTICS_IS_ON('WdRHO_P ',myThid) )
141 & doDiagsRho = doDiagsRho + 4
142 IF ( DIAGNOSTICS_IS_ON('WdRHOdP ',myThid) )
143 & doDiagsRho = doDiagsRho + 8
144 ENDIF
145 #endif /* ALLOW_DIAGNOSTICS */
146
147 #ifdef ALLOW_OBCS
148 IF (useOBCS) THEN
149 C-- Calculate future values on open boundaries
150 C-- moved before SEAICE_MODEL call since SEAICE_MODEL needs seaice-obcs fields
151 # ifdef ALLOW_AUTODIFF_TAMC
152 CADJ STORE theta = comlev1, key=ikey_dynamics, kind=isbyte
153 CADJ STORE salt = comlev1, key=ikey_dynamics, kind=isbyte
154 # endif
155 # ifdef ALLOW_DEBUG
156 IF (debugMode) CALL DEBUG_CALL('OBCS_CALC',myThid)
157 # endif
158 CALL OBCS_CALC( myTime+deltaTClock, myIter+1,
159 I uVel, vVel, wVel, theta, salt, myThid )
160 ENDIF
161 #endif /* ALLOW_OBCS */
162
163 #ifdef ALLOW_AUTODIFF_TAMC
164 # ifdef ALLOW_SALT_PLUME
165 DO bj=myByLo(myThid),myByHi(myThid)
166 DO bi=myBxLo(myThid),myBxHi(myThid)
167 DO j=1-OLy,sNy+OLy
168 DO i=1-OLx,sNx+OLx
169 saltPlumeDepth(i,j,bi,bj) = 0. _d 0
170 saltPlumeFlux(i,j,bi,bj) = 0. _d 0
171 ENDDO
172 ENDDO
173 ENDDO
174 ENDDO
175 # endif
176 #endif /* ALLOW_AUTODIFF_TAMC */
177
178 #ifdef ALLOW_FRAZIL
179 IF ( useFRAZIL ) THEN
180 C-- Freeze water in the ocean interior and let it rise to the surface
181 CALL FRAZIL_CALC_RHS( myTime, myIter, myThid )
182 ENDIF
183 #endif /* ALLOW_FRAZIL */
184
185 #ifdef ALLOW_SEAICE
186 IF ( useSEAICE ) THEN
187 # ifdef ALLOW_AUTODIFF_TAMC
188 cph-adj-test(
189 CADJ STORE area = comlev1, key=ikey_dynamics, kind=isbyte
190 CADJ STORE hsnow = comlev1, key=ikey_dynamics, kind=isbyte
191 CADJ STORE heff = comlev1, key=ikey_dynamics, kind=isbyte
192 CADJ STORE empmr,qsw,theta = comlev1, key = ikey_dynamics,
193 CADJ & kind = isbyte
194 cph-adj-test)
195 CADJ STORE atemp,aqh,precip = comlev1, key = ikey_dynamics,
196 CADJ & kind = isbyte
197 CADJ STORE swdown,lwdown = comlev1, key = ikey_dynamics,
198 CADJ & kind = isbyte
199 cph# ifdef EXF_READ_EVAP
200 CADJ STORE evap = comlev1, key = ikey_dynamics,
201 CADJ & kind = isbyte
202 cph# endif
203 CADJ STORE uvel,vvel = comlev1, key = ikey_dynamics,
204 CADJ & kind = isbyte
205 # ifdef SEAICE_CGRID
206 CADJ STORE stressdivergencex = comlev1, key = ikey_dynamics,
207 CADJ & kind = isbyte
208 CADJ STORE stressdivergencey = comlev1, key = ikey_dynamics,
209 CADJ & kind = isbyte
210 # endif
211 # ifdef SEAICE_ALLOW_DYNAMICS
212 CADJ STORE uice = comlev1, key = ikey_dynamics,
213 CADJ & kind = isbyte
214 CADJ STORE vice = comlev1, key = ikey_dynamics,
215 CADJ & kind = isbyte
216 # ifdef SEAICE_ALLOW_EVP
217 CADJ STORE seaice_sigma1 = comlev1, key = ikey_dynamics,
218 CADJ & kind = isbyte
219 CADJ STORE seaice_sigma2 = comlev1, key = ikey_dynamics,
220 CADJ & kind = isbyte
221 CADJ STORE seaice_sigma12 = comlev1, key = ikey_dynamics,
222 CADJ & kind = isbyte
223 # endif
224 # endif
225 cph# ifdef SEAICE_SALINITY
226 CADJ STORE salt = comlev1, key = ikey_dynamics,
227 CADJ & kind = isbyte
228 cph# endif
229 # ifdef ATMOSPHERIC_LOADING
230 CADJ STORE pload = comlev1, key = ikey_dynamics,
231 CADJ & kind = isbyte
232 CADJ STORE siceload = comlev1, key = ikey_dynamics,
233 CADJ & kind = isbyte
234 # endif
235 # ifdef NONLIN_FRSURF
236 CADJ STORE recip_hfacc = comlev1, key = ikey_dynamics,
237 CADJ & kind = isbyte
238 # endif
239 # ifdef ANNUAL_BALANCE
240 CADJ STORE balance_itcount = comlev1, key = ikey_dynamics,
241 CADJ & kind = isbyte
242 # endif /* ANNUAL_BALANCE */
243 # endif
244 # ifdef ALLOW_DEBUG
245 IF (debugMode) CALL DEBUG_CALL('SEAICE_MODEL',myThid)
246 # endif
247 CALL TIMER_START('SEAICE_MODEL [DO_OCEANIC_PHYS]', myThid)
248 CALL SEAICE_MODEL( myTime, myIter, myThid )
249 CALL TIMER_STOP ('SEAICE_MODEL [DO_OCEANIC_PHYS]', myThid)
250 # ifdef ALLOW_COST
251 CALL SEAICE_COST_SENSI ( myTime, myIter, myThid )
252 # endif
253 ENDIF
254 #endif /* ALLOW_SEAICE */
255
256 #ifdef ALLOW_AUTODIFF_TAMC
257 CADJ STORE sst, sss = comlev1, key = ikey_dynamics,
258 CADJ & kind = isbyte
259 CADJ STORE qsw = comlev1, key = ikey_dynamics,
260 CADJ & kind = isbyte
261 # ifdef ALLOW_SEAICE
262 CADJ STORE area = comlev1, key = ikey_dynamics,
263 CADJ & kind = isbyte
264 # endif
265 #endif
266
267 #if (defined ALLOW_THSICE) && !(defined ALLOW_ATM2D)
268 IF ( useThSIce .AND. fluidIsWater ) THEN
269 # ifdef ALLOW_AUTODIFF_TAMC
270 cph(
271 # ifdef NONLIN_FRSURF
272 CADJ STORE uice,vice = comlev1, key = ikey_dynamics,
273 CADJ & kind = isbyte
274 CADJ STORE salt,theta = comlev1, key = ikey_dynamics,
275 CADJ & kind = isbyte
276 CADJ STORE qnet,qsw, empmr = comlev1, key = ikey_dynamics,
277 CADJ & kind = isbyte
278 CADJ STORE hFac_surfC = comlev1, key = ikey_dynamics,
279 CADJ & kind = isbyte
280 # endif
281 # endif
282 # ifdef ALLOW_DEBUG
283 IF (debugMode) CALL DEBUG_CALL('THSICE_MAIN',myThid)
284 # endif
285 C-- Step forward Therm.Sea-Ice variables
286 C and modify forcing terms including effects from ice
287 CALL TIMER_START('THSICE_MAIN [DO_OCEANIC_PHYS]', myThid)
288 CALL THSICE_MAIN( myTime, myIter, myThid )
289 CALL TIMER_STOP( 'THSICE_MAIN [DO_OCEANIC_PHYS]', myThid)
290 ENDIF
291 #endif /* ALLOW_THSICE */
292
293 #ifdef ALLOW_SHELFICE
294 # ifdef ALLOW_AUTODIFF_TAMC
295 CADJ STORE salt, theta = comlev1, key = ikey_dynamics,
296 CADJ & kind = isbyte
297 # endif
298 IF ( useShelfIce .AND. fluidIsWater ) THEN
299 #ifdef ALLOW_DEBUG
300 IF (debugMode) CALL DEBUG_CALL('SHELFICE_THERMODYNAMICS',myThid)
301 #endif
302 C compute temperature and (virtual) salt flux at the
303 C shelf-ice ocean interface
304 CALL TIMER_START('SHELFICE_THERMODYNAMICS [DO_OCEANIC_PHYS]',
305 & myThid)
306 CALL SHELFICE_THERMODYNAMICS( myTime, myIter, myThid )
307 CALL TIMER_STOP( 'SHELFICE_THERMODYNAMICS [DO_OCEANIC_PHYS]',
308 & myThid)
309 ENDIF
310 #endif /* ALLOW_SHELFICE */
311
312 #ifdef ALLOW_ICEFRONT
313 IF ( useICEFRONT .AND. fluidIsWater ) THEN
314 #ifdef ALLOW_DEBUG
315 IF (debugMode) CALL DEBUG_CALL('ICEFRONT_THERMODYNAMICS',myThid)
316 #endif
317 C compute temperature and (virtual) salt flux at the
318 C ice-front ocean interface
319 CALL TIMER_START('ICEFRONT_THERMODYNAMICS [DO_OCEANIC_PHYS]',
320 & myThid)
321 CALL ICEFRONT_THERMODYNAMICS( myTime, myIter, myThid )
322 CALL TIMER_STOP( 'ICEFRONT_THERMODYNAMICS [DO_OCEANIC_PHYS]',
323 & myThid)
324 ENDIF
325 #endif /* ALLOW_ICEFRONT */
326
327 #ifdef ALLOW_SALT_PLUME
328 IF ( useSALT_PLUME ) THEN
329 CALL SALT_PLUME_DO_EXCH( myTime, myIter, myThid )
330 ENDIF
331 #endif /* ALLOW_SALT_PLUME */
332
333 C-- Freeze water at the surface
334 IF ( allowFreezing ) THEN
335 #ifdef ALLOW_AUTODIFF_TAMC
336 CADJ STORE theta = comlev1, key = ikey_dynamics,
337 CADJ & kind = isbyte
338 #endif
339 CALL FREEZE_SURFACE( myTime, myIter, myThid )
340 ENDIF
341
342 #ifdef ALLOW_OCN_COMPON_INTERF
343 C-- Apply imported data (from coupled interface) to forcing fields
344 C jmc: do not know precisely where to put this call (bf or af thSIce ?)
345 IF ( useCoupler ) THEN
346 CALL OCN_APPLY_IMPORT( .TRUE., myTime, myIter, myThid )
347 ENDIF
348 #endif /* ALLOW_OCN_COMPON_INTERF */
349
350 #ifdef ALLOW_BALANCE_FLUXES
351 C balance fluxes
352 IF ( balanceEmPmR .AND. (.NOT.useSeaice .OR. useThSIce) )
353 & CALL REMOVE_MEAN_RS( 1, EmPmR, maskInC, maskInC, rA, drF,
354 & 'EmPmR', myTime, myThid )
355 IF ( balanceQnet .AND. (.NOT.useSeaice .OR. useThSIce) )
356 & CALL REMOVE_MEAN_RS( 1, Qnet, maskInC, maskInC, rA, drF,
357 & 'Qnet ', myTime, myThid )
358 #endif /* ALLOW_BALANCE_FLUXES */
359
360 #ifdef ALLOW_AUTODIFF_TAMC
361 C-- HPF directive to help TAMC
362 CHPF$ INDEPENDENT
363 #else /* ALLOW_AUTODIFF_TAMC */
364 C if fluid is not water, by-pass find_rho, gmredi, surfaceForcing
365 C and all vertical mixing schemes, but keep OBCS_CALC
366 IF ( fluidIsWater ) THEN
367 #endif /* ALLOW_AUTODIFF_TAMC */
368 DO bj=myByLo(myThid),myByHi(myThid)
369 #ifdef ALLOW_AUTODIFF_TAMC
370 C-- HPF directive to help TAMC
371 CHPF$ INDEPENDENT
372 #endif /* ALLOW_AUTODIFF_TAMC */
373 DO bi=myBxLo(myThid),myBxHi(myThid)
374
375 #ifdef ALLOW_AUTODIFF_TAMC
376 act1 = bi - myBxLo(myThid)
377 max1 = myBxHi(myThid) - myBxLo(myThid) + 1
378 act2 = bj - myByLo(myThid)
379 max2 = myByHi(myThid) - myByLo(myThid) + 1
380 act3 = myThid - 1
381 max3 = nTx*nTy
382 act4 = ikey_dynamics - 1
383 itdkey = (act1 + 1) + act2*max1
384 & + act3*max1*max2
385 & + act4*max1*max2*max3
386 #endif /* ALLOW_AUTODIFF_TAMC */
387
388 C-- Set up work arrays with valid (i.e. not NaN) values
389 C These inital values do not alter the numerical results. They
390 C just ensure that all memory references are to valid floating
391 C point numbers. This prevents spurious hardware signals due to
392 C uninitialised but inert locations.
393
394 #ifdef ALLOW_AUTODIFF_TAMC
395 DO j=1-OLy,sNy+OLy
396 DO i=1-OLx,sNx+OLx
397 rhoKm1 (i,j) = 0. _d 0
398 rhoKp1 (i,j) = 0. _d 0
399 ENDDO
400 ENDDO
401 #endif /* ALLOW_AUTODIFF_TAMC */
402
403 DO k=1,Nr
404 DO j=1-OLy,sNy+OLy
405 DO i=1-OLx,sNx+OLx
406 C This is currently used by GMRedi, IVDC, MXL-depth and Diagnostics
407 sigmaX(i,j,k) = 0. _d 0
408 sigmaY(i,j,k) = 0. _d 0
409 sigmaR(i,j,k) = 0. _d 0
410 #ifdef ALLOW_AUTODIFF_TAMC
411 cph all the following init. are necessary for TAF
412 cph although some of these are re-initialised later.
413 rhoInSitu(i,j,k,bi,bj) = 0.
414 IVDConvCount(i,j,k,bi,bj) = 0.
415 # ifdef ALLOW_GMREDI
416 Kwx(i,j,k,bi,bj) = 0. _d 0
417 Kwy(i,j,k,bi,bj) = 0. _d 0
418 Kwz(i,j,k,bi,bj) = 0. _d 0
419 # ifdef GM_NON_UNITY_DIAGONAL
420 Kux(i,j,k,bi,bj) = 0. _d 0
421 Kvy(i,j,k,bi,bj) = 0. _d 0
422 # endif
423 # ifdef GM_EXTRA_DIAGONAL
424 Kuz(i,j,k,bi,bj) = 0. _d 0
425 Kvz(i,j,k,bi,bj) = 0. _d 0
426 # endif
427 # ifdef GM_BOLUS_ADVEC
428 GM_PsiX(i,j,k,bi,bj) = 0. _d 0
429 GM_PsiY(i,j,k,bi,bj) = 0. _d 0
430 # endif
431 # ifdef GM_VISBECK_VARIABLE_K
432 VisbeckK(i,j,bi,bj) = 0. _d 0
433 # endif
434 # endif /* ALLOW_GMREDI */
435 # ifdef ALLOW_KPP
436 KPPdiffKzS(i,j,k,bi,bj) = 0. _d 0
437 KPPdiffKzT(i,j,k,bi,bj) = 0. _d 0
438 # endif /* ALLOW_KPP */
439 # ifdef ALLOW_GGL90
440 GGL90viscArU(i,j,k,bi,bj) = 0. _d 0
441 GGL90viscArV(i,j,k,bi,bj) = 0. _d 0
442 GGL90diffKr(i,j,k,bi,bj) = 0. _d 0
443 # endif /* ALLOW_GGL90 */
444 #endif /* ALLOW_AUTODIFF_TAMC */
445 ENDDO
446 ENDDO
447 ENDDO
448
449 iMin = 1-OLx
450 iMax = sNx+OLx
451 jMin = 1-OLy
452 jMax = sNy+OLy
453
454 #ifdef ALLOW_AUTODIFF_TAMC
455 CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey,
456 CADJ & kind = isbyte
457 CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey,
458 CADJ & kind = isbyte
459 CADJ STORE totphihyd(:,:,:,bi,bj)
460 CADJ & = comlev1_bibj, key=itdkey,
461 CADJ & kind = isbyte
462 # ifdef ALLOW_KPP
463 CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey,
464 CADJ & kind = isbyte
465 CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey,
466 CADJ & kind = isbyte
467 # endif
468 # ifdef ALLOW_SALT_PLUME
469 CADJ STORE saltplumedepth(:,:,bi,bj) = comlev1_bibj, key=itdkey,
470 CADJ & kind = isbyte
471 # endif
472 #endif /* ALLOW_AUTODIFF_TAMC */
473
474 C-- Always compute density (stored in common block) here; even when it is not
475 C needed here, will be used anyway in calc_phi_hyd (data flow easier this way)
476 #ifdef ALLOW_DEBUG
477 IF (debugMode) CALL DEBUG_CALL('FIND_RHO_2D (xNr)',myThid)
478 #endif
479 #ifdef ALLOW_AUTODIFF_TAMC
480 IF ( fluidIsWater ) THEN
481 #endif /* ALLOW_AUTODIFF_TAMC */
482 #ifdef ALLOW_DOWN_SLOPE
483 IF ( useDOWN_SLOPE ) THEN
484 DO k=1,Nr
485 CALL DWNSLP_CALC_RHO(
486 I theta, salt,
487 O rhoInSitu(1-OLx,1-OLy,k,bi,bj),
488 I k, bi, bj, myTime, myIter, myThid )
489 ENDDO
490 ENDIF
491 #endif /* ALLOW_DOWN_SLOPE */
492 #ifdef ALLOW_BBL
493 IF ( useBBL ) THEN
494 C pkg/bbl requires in-situ bbl density for depths equal to and deeper than the bbl.
495 C To reduce computation and storage requirement, these densities are stored in the
496 C dry grid boxes of rhoInSitu. See BBL_CALC_RHO for details.
497 DO k=Nr,1,-1
498 CALL BBL_CALC_RHO(
499 I theta, salt,
500 O rhoInSitu,
501 I k, bi, bj, myTime, myIter, myThid )
502
503 ENDDO
504 ENDIF
505 #endif /* ALLOW_BBL */
506 IF ( .NOT. ( useDOWN_SLOPE .OR. useBBL ) ) THEN
507 DO k=1,Nr
508 CALL FIND_RHO_2D(
509 I iMin, iMax, jMin, jMax, k,
510 I theta(1-OLx,1-OLy,k,bi,bj),
511 I salt (1-OLx,1-OLy,k,bi,bj),
512 O rhoInSitu(1-OLx,1-OLy,k,bi,bj),
513 I k, bi, bj, myThid )
514 ENDDO
515 ENDIF
516 #ifdef ALLOW_AUTODIFF_TAMC
517 ELSE
518 C- fluid is not water:
519 DO k=1,Nr
520 DO j=1-OLy,sNy+OLy
521 DO i=1-OLx,sNx+OLx
522 rhoInSitu(i,j,k,bi,bj) = 0.
523 ENDDO
524 ENDDO
525 ENDDO
526 ENDIF
527 #endif /* ALLOW_AUTODIFF_TAMC */
528
529 #ifdef ALLOW_DEBUG
530 IF (debugMode) CALL DEBUG_MSG('ENTERING UPWARD K LOOP',myThid)
531 #endif
532
533 C-- Start of diagnostic loop
534 DO k=Nr,1,-1
535
536 #ifdef ALLOW_AUTODIFF_TAMC
537 C? Patrick, is this formula correct now that we change the loop range?
538 C? Do we still need this?
539 cph kkey formula corrected.
540 cph Needed for rhoK, rhoKm1, in the case useGMREDI.
541 kkey = (itdkey-1)*Nr + k
542 #endif /* ALLOW_AUTODIFF_TAMC */
543
544 c#ifdef ALLOW_AUTODIFF_TAMC
545 cCADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey,
546 cCADJ & kind = isbyte
547 cCADJ STORE salt(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey,
548 cCADJ & kind = isbyte
549 c#endif /* ALLOW_AUTODIFF_TAMC */
550
551 C-- Calculate gradients of potential density for isoneutral
552 C slope terms (e.g. GM/Redi tensor or IVDC diffusivity)
553 IF ( useGMRedi .OR. (k.GT.1 .AND. ivdc_kappa.NE.0.)
554 & .OR. usePP81 .OR. useMY82 .OR. useGGL90
555 & .OR. useSALT_PLUME .OR. doDiagsRho.GE.1 ) THEN
556 IF (k.GT.1) THEN
557 #ifdef ALLOW_AUTODIFF_TAMC
558 CADJ STORE theta(:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey,
559 CADJ & kind = isbyte
560 CADJ STORE salt (:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey,
561 CADJ & kind = isbyte
562 CADJ STORE rhokm1 (bi,bj) = comlev1_bibj_k, key=kkey,
563 CADJ & kind = isbyte
564 #endif /* ALLOW_AUTODIFF_TAMC */
565 CALL FIND_RHO_2D(
566 I iMin, iMax, jMin, jMax, k,
567 I theta(1-OLx,1-OLy,k-1,bi,bj),
568 I salt (1-OLx,1-OLy,k-1,bi,bj),
569 O rhoKm1,
570 I k-1, bi, bj, myThid )
571 ENDIF
572 #ifdef ALLOW_DEBUG
573 IF (debugMode) CALL DEBUG_CALL('GRAD_SIGMA',myThid)
574 #endif
575 cph Avoid variable aliasing for adjoint !!!
576 DO j=jMin,jMax
577 DO i=iMin,iMax
578 rhoKp1(i,j) = rhoInSitu(i,j,k,bi,bj)
579 ENDDO
580 ENDDO
581 CALL GRAD_SIGMA(
582 I bi, bj, iMin, iMax, jMin, jMax, k,
583 I rhoInSitu(1-OLx,1-OLy,k,bi,bj), rhoKm1, rhoKp1,
584 O sigmaX, sigmaY, sigmaR,
585 I myThid )
586 #ifdef ALLOW_AUTODIFF_TAMC
587 #ifdef GMREDI_WITH_STABLE_ADJOINT
588 cgf zero out adjoint fields to stabilize pkg/gmredi adjoint
589 cgf -> cuts adjoint dependency from slope to state
590 CALL ZERO_ADJ_LOC( Nr, sigmaX, myThid)
591 CALL ZERO_ADJ_LOC( Nr, sigmaY, myThid)
592 CALL ZERO_ADJ_LOC( Nr, sigmaR, myThid)
593 #endif
594 #endif /* ALLOW_AUTODIFF_TAMC */
595 ENDIF
596
597 C-- Implicit Vertical Diffusion for Convection
598 IF (k.GT.1 .AND. ivdc_kappa.NE.0.) THEN
599 #ifdef ALLOW_DEBUG
600 IF (debugMode) CALL DEBUG_CALL('CALC_IVDC',myThid)
601 #endif
602 CALL CALC_IVDC(
603 I bi, bj, iMin, iMax, jMin, jMax, k,
604 I sigmaR,
605 I myTime, myIter, myThid)
606 ENDIF
607
608 #ifdef ALLOW_DIAGNOSTICS
609 IF ( doDiagsRho.GE.4 ) THEN
610 CALL DIAGS_RHO_L( doDiagsRho, k, bi, bj,
611 I rhoInSitu(1-OLx,1-OLy,1,bi,bj),
612 I rhoKm1, wVel,
613 I myTime, myIter, myThid )
614 ENDIF
615 #endif
616
617 C-- end of diagnostic k loop (Nr:1)
618 ENDDO
619
620 #ifdef ALLOW_AUTODIFF_TAMC
621 CADJ STORE IVDConvCount(:,:,:,bi,bj)
622 CADJ & = comlev1_bibj, key=itdkey,
623 CADJ & kind = isbyte
624 #endif
625
626 C-- Diagnose Mixed Layer Depth:
627 IF ( useGMRedi .OR. MOD(doDiagsRho,2).EQ.1 ) THEN
628 CALL CALC_OCE_MXLAYER(
629 I rhoInSitu(1-OLx,1-OLy,1,bi,bj), sigmaR,
630 I bi, bj, myTime, myIter, myThid )
631 ENDIF
632
633 #ifdef ALLOW_SALT_PLUME
634 IF ( useSALT_PLUME ) THEN
635 CALL SALT_PLUME_CALC_DEPTH(
636 I rhoInSitu(1-OLx,1-OLy,1,bi,bj), sigmaR,
637 I bi, bj, myTime, myIter, myThid )
638 ENDIF
639 #endif /* ALLOW_SALT_PLUME */
640
641 #ifdef ALLOW_DIAGNOSTICS
642 IF ( MOD(doDiagsRho,4).GE.2 ) THEN
643 CALL DIAGNOSTICS_FILL (sigmaR, 'DRHODR ', 0, Nr,
644 & 2, bi, bj, myThid)
645 ENDIF
646 #endif /* ALLOW_DIAGNOSTICS */
647
648 C-- Determines forcing terms based on external fields
649 C relaxation terms, etc.
650 #ifdef ALLOW_DEBUG
651 IF (debugMode) CALL DEBUG_CALL('EXTERNAL_FORCING_SURF',myThid)
652 #endif
653 #ifdef ALLOW_AUTODIFF_TAMC
654 CADJ STORE EmPmR(:,:,bi,bj)
655 CADJ & = comlev1_bibj, key=itdkey,
656 CADJ & kind = isbyte
657 # ifdef EXACT_CONSERV
658 CADJ STORE PmEpR(:,:,bi,bj)
659 CADJ & = comlev1_bibj, key=itdkey,
660 CADJ & kind = isbyte
661 # endif
662 # ifdef NONLIN_FRSURF
663 CADJ STORE hFac_surfC(:,:,bi,bj)
664 CADJ & = comlev1_bibj, key=itdkey,
665 CADJ & kind = isbyte
666 CADJ STORE recip_hFacC(:,:,:,bi,bj)
667 CADJ & = comlev1_bibj, key=itdkey,
668 CADJ & kind = isbyte
669 # if (defined (ALLOW_PTRACERS))
670 CADJ STORE surfaceForcingS(:,:,bi,bj) = comlev1_bibj, key=itdkey,
671 CADJ & kind = isbyte
672 CADJ STORE surfaceForcingT(:,:,bi,bj) = comlev1_bibj, key=itdkey,
673 CADJ & kind = isbyte
674 # endif /* ALLOW_PTRACERS */
675 # endif /* NONLIN_FRSURF */
676 #endif
677 CALL EXTERNAL_FORCING_SURF(
678 I bi, bj, iMin, iMax, jMin, jMax,
679 I myTime, myIter, myThid )
680 #ifdef ALLOW_AUTODIFF_TAMC
681 # ifdef EXACT_CONSERV
682 cph-test
683 cphCADJ STORE PmEpR(:,:,bi,bj)
684 cphCADJ & = comlev1_bibj, key=itdkey,
685 cphCADJ & kind = isbyte
686 # endif
687 #endif
688
689 #ifdef ALLOW_AUTODIFF_TAMC
690 cph needed for KPP
691 CADJ STORE surfaceForcingU(:,:,bi,bj)
692 CADJ & = comlev1_bibj, key=itdkey,
693 CADJ & kind = isbyte
694 CADJ STORE surfaceForcingV(:,:,bi,bj)
695 CADJ & = comlev1_bibj, key=itdkey,
696 CADJ & kind = isbyte
697 CADJ STORE surfaceForcingS(:,:,bi,bj)
698 CADJ & = comlev1_bibj, key=itdkey,
699 CADJ & kind = isbyte
700 CADJ STORE surfaceForcingT(:,:,bi,bj)
701 CADJ & = comlev1_bibj, key=itdkey,
702 CADJ & kind = isbyte
703 CADJ STORE surfaceForcingTice(:,:,bi,bj)
704 CADJ & = comlev1_bibj, key=itdkey,
705 CADJ & kind = isbyte
706 #endif /* ALLOW_AUTODIFF_TAMC */
707
708 #ifdef ALLOW_KPP
709 C-- Compute KPP mixing coefficients
710 IF (useKPP) THEN
711 #ifdef ALLOW_DEBUG
712 IF (debugMode) CALL DEBUG_CALL('KPP_CALC',myThid)
713 #endif
714 CALL TIMER_START('KPP_CALC [DO_OCEANIC_PHYS]', myThid)
715 CALL KPP_CALC(
716 I bi, bj, myTime, myIter, myThid )
717 CALL TIMER_STOP ('KPP_CALC [DO_OCEANIC_PHYS]', myThid)
718 #ifdef ALLOW_AUTODIFF_TAMC
719 ELSE
720 CALL KPP_CALC_DUMMY(
721 I bi, bj, myTime, myIter, myThid )
722 #endif /* ALLOW_AUTODIFF_TAMC */
723 ENDIF
724
725 #endif /* ALLOW_KPP */
726
727 #ifdef ALLOW_PP81
728 C-- Compute PP81 mixing coefficients
729 IF (usePP81) THEN
730 #ifdef ALLOW_DEBUG
731 IF (debugMode) CALL DEBUG_CALL('PP81_CALC',myThid)
732 #endif
733 CALL PP81_CALC(
734 I bi, bj, sigmaR, myTime, myIter, myThid )
735 ENDIF
736 #endif /* ALLOW_PP81 */
737
738 #ifdef ALLOW_MY82
739 C-- Compute MY82 mixing coefficients
740 IF (useMY82) THEN
741 #ifdef ALLOW_DEBUG
742 IF (debugMode) CALL DEBUG_CALL('MY82_CALC',myThid)
743 #endif
744 CALL MY82_CALC(
745 I bi, bj, sigmaR, myTime, myIter, myThid )
746 ENDIF
747 #endif /* ALLOW_MY82 */
748
749 #ifdef ALLOW_GGL90
750 #ifdef ALLOW_AUTODIFF_TAMC
751 CADJ STORE GGL90TKE (:,:,:,bi,bj) = comlev1_bibj, key=itdkey,
752 CADJ & kind = isbyte
753 #endif /* ALLOW_AUTODIFF_TAMC */
754 C-- Compute GGL90 mixing coefficients
755 IF (useGGL90) THEN
756 #ifdef ALLOW_DEBUG
757 IF (debugMode) CALL DEBUG_CALL('GGL90_CALC',myThid)
758 #endif
759 CALL TIMER_START('GGL90_CALC [DO_OCEANIC_PHYS]', myThid)
760 CALL GGL90_CALC(
761 I bi, bj, sigmaR, myTime, myIter, myThid )
762 CALL TIMER_STOP ('GGL90_CALC [DO_OCEANIC_PHYS]', myThid)
763 ENDIF
764 #endif /* ALLOW_GGL90 */
765
766 #ifdef ALLOW_TIMEAVE
767 IF ( taveFreq.GT. 0. _d 0 ) THEN
768 CALL TIMEAVE_SURF_FLUX( bi, bj, myTime, myIter, myThid)
769 ENDIF
770 IF (taveFreq.GT.0. .AND. ivdc_kappa.NE.0.) THEN
771 CALL TIMEAVE_CUMULATE(ConvectCountTave, IVDConvCount,
772 I Nr, deltaTClock, bi, bj, myThid)
773 ENDIF
774 #endif /* ALLOW_TIMEAVE */
775
776 #ifdef ALLOW_GMREDI
777 #ifdef ALLOW_AUTODIFF_TAMC
778 # ifndef GM_EXCLUDE_CLIPPING
779 cph storing here is needed only for one GMREDI_OPTIONS:
780 cph define GM_BOLUS_ADVEC
781 cph keep it although TAF says you dont need to.
782 cph but I have avoided the #ifdef for now, in case more things change
783 CADJ STORE sigmaX(:,:,:) = comlev1_bibj, key=itdkey,
784 CADJ & kind = isbyte
785 CADJ STORE sigmaY(:,:,:) = comlev1_bibj, key=itdkey,
786 CADJ & kind = isbyte
787 CADJ STORE sigmaR(:,:,:) = comlev1_bibj, key=itdkey,
788 CADJ & kind = isbyte
789 # endif
790 #endif /* ALLOW_AUTODIFF_TAMC */
791
792 C-- Calculate iso-neutral slopes for the GM/Redi parameterisation
793 IF (useGMRedi) THEN
794 #ifdef ALLOW_DEBUG
795 IF (debugMode) CALL DEBUG_CALL('GMREDI_CALC_TENSOR',myThid)
796 #endif
797 CALL GMREDI_CALC_TENSOR(
798 I iMin, iMax, jMin, jMax,
799 I sigmaX, sigmaY, sigmaR,
800 I bi, bj, myTime, myIter, myThid )
801 #ifdef ALLOW_AUTODIFF_TAMC
802 ELSE
803 CALL GMREDI_CALC_TENSOR_DUMMY(
804 I iMin, iMax, jMin, jMax,
805 I sigmaX, sigmaY, sigmaR,
806 I bi, bj, myTime, myIter, myThid )
807 #endif /* ALLOW_AUTODIFF_TAMC */
808 ENDIF
809 #endif /* ALLOW_GMREDI */
810
811 #ifdef ALLOW_DOWN_SLOPE
812 IF ( useDOWN_SLOPE ) THEN
813 C-- Calculate Downsloping Flow for Down_Slope parameterization
814 IF ( usingPCoords ) THEN
815 CALL DWNSLP_CALC_FLOW(
816 I bi, bj, kSurfC, rhoInSitu,
817 I myTime, myIter, myThid )
818 ELSE
819 CALL DWNSLP_CALC_FLOW(
820 I bi, bj, kLowC, rhoInSitu,
821 I myTime, myIter, myThid )
822 ENDIF
823 ENDIF
824 #endif /* ALLOW_DOWN_SLOPE */
825
826 C-- end bi,bj loops.
827 ENDDO
828 ENDDO
829
830 #ifdef ALLOW_BALANCE_RELAX
831 # ifdef ALLOW_AUTODIFF_TAMC
832 CADJ STORE SSSrlx = comlev1, key=ikey_dynamics, kind=isbyte
833 CADJ STORE SSSrlxTile = comlev1, key=ikey_dynamics, kind=isbyte
834 CADJ STORE SSSrlxGlob = comlev1, key=ikey_dynamics, kind=isbyte
835 CADJ STORE SSTrlx = comlev1, key=ikey_dynamics, kind=isbyte
836 CADJ STORE SSTrlxTile = comlev1, key=ikey_dynamics, kind=isbyte
837 CADJ STORE SSTrlxGlob = comlev1, key=ikey_dynamics, kind=isbyte
838 # endif /* ALLOW_AUTODIFF_TAMC */
839 CALL BALANCE_RELAX( myTime, myIter, myThid )
840 #endif /* ALLOW_BALANCE_RELAX */
841
842 #ifndef ALLOW_AUTODIFF_TAMC
843 C--- if fluid Is Water: end
844 ENDIF
845 #endif
846
847 #ifdef ALLOW_BBL
848 IF ( useBBL ) THEN
849 CALL BBL_CALC_RHS(
850 I myTime, myIter, myThid )
851 ENDIF
852 #endif /* ALLOW_BBL */
853
854 #ifdef ALLOW_MYPACKAGE
855 IF ( useMYPACKAGE ) THEN
856 CALL MYPACKAGE_CALC_RHS(
857 I myTime, myIter, myThid )
858 ENDIF
859 #endif /* ALLOW_MYPACKAGE */
860
861 #ifdef ALLOW_GMREDI
862 IF ( useGMRedi ) THEN
863 CALL GMREDI_DO_EXCH( myTime, myIter, myThid )
864 ENDIF
865 #endif /* ALLOW_GMREDI */
866
867 #ifdef ALLOW_KPP
868 IF (useKPP) THEN
869 CALL KPP_DO_EXCH( myThid )
870 ENDIF
871 #endif /* ALLOW_KPP */
872
873 #ifdef ALLOW_DIAGNOSTICS
874 IF ( fluidIsWater .AND. useDiagnostics ) THEN
875 CALL DIAGS_RHO_G(
876 I rhoInSitu, uVel, vVel, wVel,
877 I myTime, myIter, myThid )
878 CALL DIAGS_OCEANIC_SURF_FLUX( myTime, myIter, myThid )
879 ENDIF
880 IF ( ivdc_kappa.NE.0 .AND. useDiagnostics ) THEN
881 CALL DIAGNOSTICS_FILL( IVDConvCount, 'CONVADJ ',
882 & 0, Nr, 0, 1, 1, myThid )
883 ENDIF
884 #endif
885
886 #ifdef ALLOW_DEBUG
887 IF (debugMode) CALL DEBUG_LEAVE('DO_OCEANIC_PHYS',myThid)
888 #endif
889
890 RETURN
891 END

  ViewVC Help
Powered by ViewVC 1.1.22