/[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.113 - (show annotations) (download)
Fri Mar 2 01:45:22 2012 UTC (12 years, 3 months ago) by dimitri
Branch: MAIN
Changes since 1.112: +8 -14 lines
adding pkg/frazil - see frazil_description.tex for details

1 C $Header: /u/gcmpack/MITgcm/model/src/do_oceanic_phys.F,v 1.112 2012/02/04 14:52:08 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_FRAZIL
190 IF ( useFRAZIL ) THEN
191 C-- Freeze water in the ocean interior and let it rise to the surface
192 CALL FRAZIL_CALC_RHS( myTime, myIter, myThid )
193 ENDIF
194 #endif /* ALLOW_FRAZIL */
195
196 #ifdef ALLOW_SEAICE
197 IF ( useSEAICE ) THEN
198 # ifdef ALLOW_AUTODIFF_TAMC
199 cph-adj-test(
200 CADJ STORE area = comlev1, key=ikey_dynamics, kind=isbyte
201 CADJ STORE hsnow = comlev1, key=ikey_dynamics, kind=isbyte
202 CADJ STORE heff = comlev1, key=ikey_dynamics, kind=isbyte
203 CADJ STORE empmr,qsw,theta = comlev1, key = ikey_dynamics,
204 CADJ & kind = isbyte
205 cph-adj-test)
206 CADJ STORE atemp,aqh,precip = comlev1, key = ikey_dynamics,
207 CADJ & kind = isbyte
208 CADJ STORE swdown,lwdown = comlev1, key = ikey_dynamics,
209 CADJ & kind = isbyte
210 cph# ifdef EXF_READ_EVAP
211 CADJ STORE evap = comlev1, key = ikey_dynamics,
212 CADJ & kind = isbyte
213 cph# endif
214 CADJ STORE uvel,vvel = comlev1, key = ikey_dynamics,
215 CADJ & kind = isbyte
216 # ifdef SEAICE_CGRID
217 CADJ STORE stressdivergencex = comlev1, key = ikey_dynamics,
218 CADJ & kind = isbyte
219 CADJ STORE stressdivergencey = comlev1, key = ikey_dynamics,
220 CADJ & kind = isbyte
221 # endif
222 # ifdef SEAICE_ALLOW_DYNAMICS
223 CADJ STORE uice = comlev1, key = ikey_dynamics,
224 CADJ & kind = isbyte
225 CADJ STORE vice = comlev1, key = ikey_dynamics,
226 CADJ & kind = isbyte
227 # ifdef SEAICE_ALLOW_EVP
228 CADJ STORE seaice_sigma1 = comlev1, key = ikey_dynamics,
229 CADJ & kind = isbyte
230 CADJ STORE seaice_sigma2 = comlev1, key = ikey_dynamics,
231 CADJ & kind = isbyte
232 CADJ STORE seaice_sigma12 = comlev1, key = ikey_dynamics,
233 CADJ & kind = isbyte
234 # endif
235 # endif
236 cph# ifdef SEAICE_SALINITY
237 CADJ STORE salt = comlev1, key = ikey_dynamics,
238 CADJ & kind = isbyte
239 cph# endif
240 # ifdef ATMOSPHERIC_LOADING
241 CADJ STORE pload = comlev1, key = ikey_dynamics,
242 CADJ & kind = isbyte
243 CADJ STORE siceload = comlev1, key = ikey_dynamics,
244 CADJ & kind = isbyte
245 # endif
246 # ifdef NONLIN_FRSURF
247 CADJ STORE recip_hfacc = comlev1, key = ikey_dynamics,
248 CADJ & kind = isbyte
249 # endif
250 # ifdef ANNUAL_BALANCE
251 CADJ STORE balance_itcount = comlev1, key = ikey_dynamics,
252 CADJ & kind = isbyte
253 # endif /* ANNUAL_BALANCE */
254 # endif
255 # ifdef ALLOW_DEBUG
256 IF (debugMode) CALL DEBUG_CALL('SEAICE_MODEL',myThid)
257 # endif
258 CALL TIMER_START('SEAICE_MODEL [DO_OCEANIC_PHYS]', myThid)
259 CALL SEAICE_MODEL( myTime, myIter, myThid )
260 CALL TIMER_STOP ('SEAICE_MODEL [DO_OCEANIC_PHYS]', myThid)
261 # ifdef ALLOW_COST
262 CALL SEAICE_COST_SENSI ( myTime, myIter, myThid )
263 # endif
264 ENDIF
265 #endif /* ALLOW_SEAICE */
266
267 #ifdef ALLOW_AUTODIFF_TAMC
268 CADJ STORE sst, sss = comlev1, key = ikey_dynamics,
269 CADJ & kind = isbyte
270 CADJ STORE qsw = comlev1, key = ikey_dynamics,
271 CADJ & kind = isbyte
272 # ifdef ALLOW_SEAICE
273 CADJ STORE area = comlev1, key = ikey_dynamics,
274 CADJ & kind = isbyte
275 # endif
276 #endif
277
278 #if (defined ALLOW_THSICE) && !(defined ALLOW_ATM2D)
279 IF ( useThSIce .AND. fluidIsWater ) THEN
280 # ifdef ALLOW_AUTODIFF_TAMC
281 cph(
282 # ifdef NONLIN_FRSURF
283 CADJ STORE uice,vice = comlev1, key = ikey_dynamics,
284 CADJ & kind = isbyte
285 CADJ STORE salt,theta = comlev1, key = ikey_dynamics,
286 CADJ & kind = isbyte
287 CADJ STORE qnet,qsw, empmr = comlev1, key = ikey_dynamics,
288 CADJ & kind = isbyte
289 CADJ STORE hFac_surfC = comlev1, key = ikey_dynamics,
290 CADJ & kind = isbyte
291 # endif
292 # endif
293 # ifdef ALLOW_DEBUG
294 IF (debugMode) CALL DEBUG_CALL('THSICE_MAIN',myThid)
295 # endif
296 C-- Step forward Therm.Sea-Ice variables
297 C and modify forcing terms including effects from ice
298 CALL TIMER_START('THSICE_MAIN [DO_OCEANIC_PHYS]', myThid)
299 CALL THSICE_MAIN( myTime, myIter, myThid )
300 CALL TIMER_STOP( 'THSICE_MAIN [DO_OCEANIC_PHYS]', myThid)
301 ENDIF
302 #endif /* ALLOW_THSICE */
303
304 #ifdef ALLOW_SHELFICE
305 # ifdef ALLOW_AUTODIFF_TAMC
306 CADJ STORE salt, theta = comlev1, key = ikey_dynamics,
307 CADJ & kind = isbyte
308 # endif
309 IF ( useShelfIce .AND. fluidIsWater ) THEN
310 #ifdef ALLOW_DEBUG
311 IF (debugMode) CALL DEBUG_CALL('SHELFICE_THERMODYNAMICS',myThid)
312 #endif
313 C compute temperature and (virtual) salt flux at the
314 C shelf-ice ocean interface
315 CALL TIMER_START('SHELFICE_THERMODYNAMICS [DO_OCEANIC_PHYS]',
316 & myThid)
317 CALL SHELFICE_THERMODYNAMICS( myTime, myIter, myThid )
318 CALL TIMER_STOP( 'SHELFICE_THERMODYNAMICS [DO_OCEANIC_PHYS]',
319 & myThid)
320 ENDIF
321 #endif /* ALLOW_SHELFICE */
322
323 #ifdef ALLOW_ICEFRONT
324 IF ( useICEFRONT .AND. fluidIsWater ) THEN
325 #ifdef ALLOW_DEBUG
326 IF (debugMode) CALL DEBUG_CALL('ICEFRONT_THERMODYNAMICS',myThid)
327 #endif
328 C compute temperature and (virtual) salt flux at the
329 C ice-front ocean interface
330 CALL TIMER_START('ICEFRONT_THERMODYNAMICS [DO_OCEANIC_PHYS]',
331 & myThid)
332 CALL ICEFRONT_THERMODYNAMICS( myTime, myIter, myThid )
333 CALL TIMER_STOP( 'ICEFRONT_THERMODYNAMICS [DO_OCEANIC_PHYS]',
334 & myThid)
335 ENDIF
336 #endif /* ALLOW_ICEFRONT */
337
338 #ifdef ALLOW_SALT_PLUME
339 IF ( useSALT_PLUME ) THEN
340 CALL SALT_PLUME_DO_EXCH( myTime, myIter, myThid )
341 ENDIF
342 #endif /* ALLOW_SALT_PLUME */
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