/[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.106 - (show annotations) (download)
Wed Aug 3 00:31:48 2011 UTC (12 years, 10 months ago) by jmc
Branch: MAIN
Changes since 1.105: +14 -14 lines
move call to MYPACKAGE_CALC_RHS outside bi,bj loop (+ remove bi,bj arguments)

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

  ViewVC Help
Powered by ViewVC 1.1.22