/[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.104 - (show annotations) (download)
Tue May 3 18:05:03 2011 UTC (13 years ago) by heimbach
Branch: MAIN
CVS Tags: checkpoint63, checkpoint62z, checkpoint62y, checkpoint62x
Changes since 1.103: +8 -6 lines
Remove unpleasant #ifndef ALLOW_AUTODIFF_TAMC
(hard to track down later)

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

  ViewVC Help
Powered by ViewVC 1.1.22