/[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.111 - (show annotations) (download)
Mon Dec 19 11:47:04 2011 UTC (12 years, 5 months ago) by mlosch
Branch: MAIN
CVS Tags: checkpoint63h, checkpoint63i
Changes since 1.110: +2 -3 lines
use sigmaR in calc_ivdc rather than rhokp1 and rhokm1

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

  ViewVC Help
Powered by ViewVC 1.1.22