/[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.102 - (show annotations) (download)
Wed Apr 20 13:45:14 2011 UTC (13 years, 1 month ago) by dimitri
Branch: MAIN
Changes since 1.101: +6 -1 lines
Add allowInteriorFreezing option to check for water below freezing point
at depth and bring the negative heat anomaly to the surface level.
Modified Files: doc/tag-index, model/inc/PARAMS.h,
 model/src/do_oceanic_phys.F, ini_parms.F, set_defaults.F, thermodynamics.F
Added File: model/src/freeze_interior.F

1 C $Header: /u/gcmpack/MITgcm/model/src/do_oceanic_phys.F,v 1.101 2011/03/05 17:51:19 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_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 # ifdef SEAICE_SALINITY
224 CADJ STORE salt = comlev1, key = ikey_dynamics,
225 CADJ & kind = isbyte
226 # 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 IF ( allowInteriorFreezing ) THEN
327 CALL FREEZE_INTERIOR( myTime, myIter, myThid )
328 ENDIF
329
330 C-- Freeze water at the surface
331 #ifdef ALLOW_AUTODIFF_TAMC
332 CADJ STORE theta = comlev1, key = ikey_dynamics,
333 CADJ & kind = isbyte
334 #endif
335 IF ( allowFreezing ) THEN
336 CALL FREEZE_SURFACE( myTime, myIter, myThid )
337 ENDIF
338
339 #ifdef ALLOW_OCN_COMPON_INTERF
340 C-- Apply imported data (from coupled interface) to forcing fields
341 C jmc: do not know precisely where to put this call (bf or af thSIce ?)
342 IF ( useCoupler ) THEN
343 CALL OCN_APPLY_IMPORT( .TRUE., myTime, myIter, myThid )
344 ENDIF
345 #endif /* ALLOW_OCN_COMPON_INTERF */
346
347 #ifdef ALLOW_BALANCE_FLUXES
348 C balance fluxes
349 IF ( balanceEmPmR )
350 & CALL REMOVE_MEAN_RS( 1, EmPmR, maskInC, maskInC, rA, drF,
351 & 'EmPmR', myTime, myThid )
352 IF ( balanceQnet )
353 & CALL REMOVE_MEAN_RS( 1, Qnet, maskInC, maskInC, rA, drF,
354 & 'Qnet ', myTime, myThid )
355 #endif /* ALLOW_BALANCE_FLUXES */
356
357 #ifdef ALLOW_AUTODIFF_TAMC
358 C-- HPF directive to help TAMC
359 CHPF$ INDEPENDENT
360 #endif /* ALLOW_AUTODIFF_TAMC */
361 DO bj=myByLo(myThid),myByHi(myThid)
362 #ifdef ALLOW_AUTODIFF_TAMC
363 C-- HPF directive to help TAMC
364 CHPF$ INDEPENDENT
365 #endif /* ALLOW_AUTODIFF_TAMC */
366 DO bi=myBxLo(myThid),myBxHi(myThid)
367
368 #ifdef ALLOW_AUTODIFF_TAMC
369 act1 = bi - myBxLo(myThid)
370 max1 = myBxHi(myThid) - myBxLo(myThid) + 1
371 act2 = bj - myByLo(myThid)
372 max2 = myByHi(myThid) - myByLo(myThid) + 1
373 act3 = myThid - 1
374 max3 = nTx*nTy
375 act4 = ikey_dynamics - 1
376 itdkey = (act1 + 1) + act2*max1
377 & + act3*max1*max2
378 & + act4*max1*max2*max3
379 #else /* ALLOW_AUTODIFF_TAMC */
380 C if fluid is not water, by-pass find_rho, gmredi, surfaceForcing
381 C and all vertical mixing schemes, but keep OBCS_CALC
382 IF ( fluidIsWater ) THEN
383 #endif /* ALLOW_AUTODIFF_TAMC */
384
385 C-- Set up work arrays with valid (i.e. not NaN) values
386 C These inital values do not alter the numerical results. They
387 C just ensure that all memory references are to valid floating
388 C point numbers. This prevents spurious hardware signals due to
389 C uninitialised but inert locations.
390
391 #ifdef ALLOW_AUTODIFF_TAMC
392 DO j=1-OLy,sNy+OLy
393 DO i=1-OLx,sNx+OLx
394 rhoKm1 (i,j) = 0. _d 0
395 rhoKp1 (i,j) = 0. _d 0
396 ENDDO
397 ENDDO
398 #endif /* ALLOW_AUTODIFF_TAMC */
399
400 DO k=1,Nr
401 DO j=1-OLy,sNy+OLy
402 DO i=1-OLx,sNx+OLx
403 C This is currently used by GMRedi, IVDC, MXL-depth and Diagnostics
404 sigmaX(i,j,k) = 0. _d 0
405 sigmaY(i,j,k) = 0. _d 0
406 sigmaR(i,j,k) = 0. _d 0
407 #ifdef ALLOW_AUTODIFF_TAMC
408 cph all the following init. are necessary for TAF
409 cph although some of these are re-initialised later.
410 c rhoInSitu(i,j,k,bi,bj) = 0.
411 IVDConvCount(i,j,k,bi,bj) = 0.
412 # ifdef ALLOW_GMREDI
413 Kwx(i,j,k,bi,bj) = 0. _d 0
414 Kwy(i,j,k,bi,bj) = 0. _d 0
415 Kwz(i,j,k,bi,bj) = 0. _d 0
416 # ifdef GM_NON_UNITY_DIAGONAL
417 Kux(i,j,k,bi,bj) = 0. _d 0
418 Kvy(i,j,k,bi,bj) = 0. _d 0
419 # endif
420 # ifdef GM_EXTRA_DIAGONAL
421 Kuz(i,j,k,bi,bj) = 0. _d 0
422 Kvz(i,j,k,bi,bj) = 0. _d 0
423 # endif
424 # ifdef GM_BOLUS_ADVEC
425 GM_PsiX(i,j,k,bi,bj) = 0. _d 0
426 GM_PsiY(i,j,k,bi,bj) = 0. _d 0
427 # endif
428 # ifdef GM_VISBECK_VARIABLE_K
429 VisbeckK(i,j,bi,bj) = 0. _d 0
430 # endif
431 # endif /* ALLOW_GMREDI */
432 # ifdef ALLOW_KPP
433 KPPdiffKzS(i,j,k,bi,bj) = 0. _d 0
434 KPPdiffKzT(i,j,k,bi,bj) = 0. _d 0
435 # endif /* ALLOW_KPP */
436 # ifdef ALLOW_GGL90
437 GGL90viscArU(i,j,k,bi,bj) = 0. _d 0
438 GGL90viscArV(i,j,k,bi,bj) = 0. _d 0
439 GGL90diffKr(i,j,k,bi,bj) = 0. _d 0
440 # endif /* ALLOW_GGL90 */
441 #endif /* ALLOW_AUTODIFF_TAMC */
442 ENDDO
443 ENDDO
444 ENDDO
445
446 iMin = 1-OLx
447 iMax = sNx+OLx
448 jMin = 1-OLy
449 jMax = sNy+OLy
450
451 #ifdef ALLOW_AUTODIFF_TAMC
452 CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey,
453 CADJ & kind = isbyte
454 CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey,
455 CADJ & kind = isbyte
456 CADJ STORE totphihyd(:,:,:,bi,bj)
457 CADJ & = comlev1_bibj, key=itdkey,
458 CADJ & kind = isbyte
459 # ifdef ALLOW_KPP
460 CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey,
461 CADJ & kind = isbyte
462 CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey,
463 CADJ & kind = isbyte
464 # endif
465 #endif /* ALLOW_AUTODIFF_TAMC */
466
467 #ifdef ALLOW_DEBUG
468 IF (debugMode) CALL DEBUG_MSG('ENTERING UPWARD K LOOP',myThid)
469 #endif
470
471 C-- Start of diagnostic loop
472 DO k=Nr,1,-1
473
474 #ifdef ALLOW_AUTODIFF_TAMC
475 C? Patrick, is this formula correct now that we change the loop range?
476 C? Do we still need this?
477 cph kkey formula corrected.
478 cph Needed for rhoK, rhoKm1, in the case useGMREDI.
479 kkey = (itdkey-1)*Nr + k
480 #endif /* ALLOW_AUTODIFF_TAMC */
481
482 C-- Always compute density (stored in common block) here; even when it is not
483 C needed here, will be used anyway in calc_phi_hyd (data flow easier this way)
484 #ifdef ALLOW_DEBUG
485 IF (debugMode) CALL DEBUG_CALL('FIND_RHO_2D',myThid)
486 #endif
487 #ifdef ALLOW_AUTODIFF_TAMC
488 IF ( fluidIsWater ) THEN
489 CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey,
490 CADJ & kind = isbyte
491 CADJ STORE salt(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey,
492 CADJ & kind = isbyte
493 #endif /* ALLOW_AUTODIFF_TAMC */
494 #ifdef ALLOW_DOWN_SLOPE
495 IF ( useDOWN_SLOPE ) THEN
496 CALL DWNSLP_CALC_RHO(
497 I theta, salt,
498 O rhoInSitu(1-OLx,1-OLy,k,bi,bj),
499 I k, bi, bj, myTime, myIter, myThid )
500 ELSE
501 #endif /* ALLOW_DOWN_SLOPE */
502 CALL FIND_RHO_2D(
503 I iMin, iMax, jMin, jMax, k,
504 I theta(1-OLx,1-OLy,k,bi,bj),
505 I salt (1-OLx,1-OLy,k,bi,bj),
506 O rhoInSitu(1-OLx,1-OLy,k,bi,bj),
507 I k, bi, bj, myThid )
508 #ifdef ALLOW_DOWN_SLOPE
509 ENDIF
510 #endif /* ALLOW_DOWN_SLOPE */
511 #ifdef ALLOW_AUTODIFF_TAMC
512 ELSE
513 C- fluid is not water:
514 DO j=1-OLy,sNy+OLy
515 DO i=1-OLx,sNx+OLx
516 rhoInSitu(i,j,k,bi,bj) = 0.
517 ENDDO
518 ENDDO
519 ENDIF
520 #endif /* ALLOW_AUTODIFF_TAMC */
521
522 C-- Calculate gradients of potential density for isoneutral
523 C slope terms (e.g. GM/Redi tensor or IVDC diffusivity)
524 IF ( useGMRedi .OR. (k.GT.1 .AND. ivdc_kappa.NE.0.)
525 & .OR. useSALT_PLUME .OR. doDiagsRho.GE.1 ) THEN
526 IF (k.GT.1) THEN
527 #ifdef ALLOW_AUTODIFF_TAMC
528 CADJ STORE theta(:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey,
529 CADJ & kind = isbyte
530 CADJ STORE salt (:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey,
531 CADJ & kind = isbyte
532 CADJ STORE rhokm1 (bi,bj) = comlev1_bibj_k, key=kkey,
533 CADJ & kind = isbyte
534 #endif /* ALLOW_AUTODIFF_TAMC */
535 CALL FIND_RHO_2D(
536 I iMin, iMax, jMin, jMax, k,
537 I theta(1-OLx,1-OLy,k-1,bi,bj),
538 I salt (1-OLx,1-OLy,k-1,bi,bj),
539 O rhoKm1,
540 I k-1, bi, bj, myThid )
541 ENDIF
542 #ifdef ALLOW_DEBUG
543 IF (debugMode) CALL DEBUG_CALL('GRAD_SIGMA',myThid)
544 #endif
545 cph Avoid variable aliasing for adjoint !!!
546 DO j=jMin,jMax
547 DO i=iMin,iMax
548 rhoKp1(i,j) = rhoInSitu(i,j,k,bi,bj)
549 ENDDO
550 ENDDO
551 CALL GRAD_SIGMA(
552 I bi, bj, iMin, iMax, jMin, jMax, k,
553 I rhoInSitu(1-OLx,1-OLy,k,bi,bj), rhoKm1, rhoKp1,
554 O sigmaX, sigmaY, sigmaR,
555 I myThid )
556 #ifdef ALLOW_AUTODIFF_TAMC
557 #ifdef GMREDI_WITH_STABLE_ADJOINT
558 cgf zero out adjoint fields to stabilize pkg/gmredi adjoint
559 cgf -> cuts adjoint dependency from slope to state
560 CALL ZERO_ADJ_LOC( Nr, sigmaX, myThid)
561 CALL ZERO_ADJ_LOC( Nr, sigmaY, myThid)
562 CALL ZERO_ADJ_LOC( Nr, sigmaR, myThid)
563 #endif
564 #endif /* ALLOW_AUTODIFF_TAMC */
565 ENDIF
566
567 C-- Implicit Vertical Diffusion for Convection
568 c ==> should use sigmaR !!!
569 IF (k.GT.1 .AND. ivdc_kappa.NE.0.) THEN
570 #ifdef ALLOW_DEBUG
571 IF (debugMode) CALL DEBUG_CALL('CALC_IVDC',myThid)
572 #endif
573 CALL CALC_IVDC(
574 I bi, bj, iMin, iMax, jMin, jMax, k,
575 I rhoKm1, rhoInSitu(1-OLx,1-OLy,k,bi,bj),
576 I myTime, myIter, myThid)
577 ENDIF
578
579 #ifdef ALLOW_DIAGNOSTICS
580 IF ( MOD(doDiagsRho,2).EQ.1 ) THEN
581 CALL DIAGS_RHO_L( k, bi, bj,
582 I rhoInSitu(1-OLx,1-OLy,k,bi,bj),
583 I rhoKm1, wVel,
584 I myTime, myIter, myThid )
585 ENDIF
586 #endif
587
588 C-- end of diagnostic k loop (Nr:1)
589 ENDDO
590
591 #ifdef ALLOW_AUTODIFF_TAMC
592 CADJ STORE IVDConvCount(:,:,:,bi,bj)
593 CADJ & = comlev1_bibj, key=itdkey,
594 CADJ & kind = isbyte
595 #endif
596
597 C-- Diagnose Mixed Layer Depth:
598 IF ( useGMRedi .OR. doDiagsRho.GE.4 ) THEN
599 CALL CALC_OCE_MXLAYER(
600 I rhoInSitu(1-OLx,1-OLy,1,bi,bj), sigmaR,
601 I bi, bj, myTime, myIter, myThid )
602 ENDIF
603
604 #ifdef ALLOW_SALT_PLUME
605 IF ( useSALT_PLUME ) THEN
606 CALL SALT_PLUME_CALC_DEPTH(
607 I rhoInSitu(1-OLx,1-OLy,1,bi,bj), sigmaR,
608 I bi, bj, myTime, myIter, myThid )
609 ENDIF
610 #endif /* ALLOW_SALT_PLUME */
611
612 #ifdef ALLOW_DIAGNOSTICS
613 IF ( MOD(doDiagsRho,4).GE.2 ) THEN
614 CALL DIAGNOSTICS_FILL (sigmaR, 'DRHODR ', 0, Nr,
615 & 2, bi, bj, myThid)
616 ENDIF
617 #endif /* ALLOW_DIAGNOSTICS */
618
619 C-- Determines forcing terms based on external fields
620 C relaxation terms, etc.
621 #ifdef ALLOW_DEBUG
622 IF (debugMode) CALL DEBUG_CALL('EXTERNAL_FORCING_SURF',myThid)
623 #endif
624 #ifdef ALLOW_AUTODIFF_TAMC
625 CADJ STORE EmPmR(:,:,bi,bj)
626 CADJ & = comlev1_bibj, key=itdkey,
627 CADJ & kind = isbyte
628 # ifdef EXACT_CONSERV
629 CADJ STORE PmEpR(:,:,bi,bj)
630 CADJ & = comlev1_bibj, key=itdkey,
631 CADJ & kind = isbyte
632 # endif
633 # ifdef NONLIN_FRSURF
634 CADJ STORE hFac_surfC(:,:,bi,bj)
635 CADJ & = comlev1_bibj, key=itdkey,
636 CADJ & kind = isbyte
637 CADJ STORE recip_hFacC(:,:,:,bi,bj)
638 CADJ & = comlev1_bibj, key=itdkey,
639 CADJ & kind = isbyte
640 # if (defined (ALLOW_PTRACERS))
641 CADJ STORE surfaceForcingS(:,:,bi,bj) = comlev1_bibj, key=itdkey,
642 CADJ & kind = isbyte
643 CADJ STORE surfaceForcingT(:,:,bi,bj) = comlev1_bibj, key=itdkey,
644 CADJ & kind = isbyte
645 # endif /* ALLOW_PTRACERS */
646 # endif /* NONLIN_FRSURF */
647 #endif
648 CALL EXTERNAL_FORCING_SURF(
649 I bi, bj, iMin, iMax, jMin, jMax,
650 I myTime, myIter, myThid )
651 #ifdef ALLOW_AUTODIFF_TAMC
652 # ifdef EXACT_CONSERV
653 cph-test
654 cphCADJ STORE PmEpR(:,:,bi,bj)
655 cphCADJ & = comlev1_bibj, key=itdkey,
656 cphCADJ & kind = isbyte
657 # endif
658 #endif
659
660 #ifdef ALLOW_AUTODIFF_TAMC
661 cph needed for KPP
662 CADJ STORE surfaceForcingU(:,:,bi,bj)
663 CADJ & = comlev1_bibj, key=itdkey,
664 CADJ & kind = isbyte
665 CADJ STORE surfaceForcingV(:,:,bi,bj)
666 CADJ & = comlev1_bibj, key=itdkey,
667 CADJ & kind = isbyte
668 CADJ STORE surfaceForcingS(:,:,bi,bj)
669 CADJ & = comlev1_bibj, key=itdkey,
670 CADJ & kind = isbyte
671 CADJ STORE surfaceForcingT(:,:,bi,bj)
672 CADJ & = comlev1_bibj, key=itdkey,
673 CADJ & kind = isbyte
674 CADJ STORE surfaceForcingTice(:,:,bi,bj)
675 CADJ & = comlev1_bibj, key=itdkey,
676 CADJ & kind = isbyte
677 #endif /* ALLOW_AUTODIFF_TAMC */
678
679 #ifdef ALLOW_KPP
680 C-- Compute KPP mixing coefficients
681 IF (useKPP) THEN
682 #ifdef ALLOW_DEBUG
683 IF (debugMode) CALL DEBUG_CALL('KPP_CALC',myThid)
684 #endif
685 CALL TIMER_START('KPP_CALC [DO_OCEANIC_PHYS]', myThid)
686 CALL KPP_CALC(
687 I bi, bj, myTime, myIter, myThid )
688 CALL TIMER_STOP ('KPP_CALC [DO_OCEANIC_PHYS]', myThid)
689 #ifdef ALLOW_AUTODIFF_TAMC
690 ELSE
691 CALL KPP_CALC_DUMMY(
692 I bi, bj, myTime, myIter, myThid )
693 #endif /* ALLOW_AUTODIFF_TAMC */
694 ENDIF
695
696 #endif /* ALLOW_KPP */
697
698 #ifdef ALLOW_PP81
699 C-- Compute PP81 mixing coefficients
700 IF (usePP81) THEN
701 #ifdef ALLOW_DEBUG
702 IF (debugMode) CALL DEBUG_CALL('PP81_CALC',myThid)
703 #endif
704 CALL PP81_CALC(
705 I bi, bj, myTime, myThid )
706 ENDIF
707 #endif /* ALLOW_PP81 */
708
709 #ifdef ALLOW_MY82
710 C-- Compute MY82 mixing coefficients
711 IF (useMY82) THEN
712 #ifdef ALLOW_DEBUG
713 IF (debugMode) CALL DEBUG_CALL('MY82_CALC',myThid)
714 #endif
715 CALL MY82_CALC(
716 I bi, bj, myTime, myThid )
717 ENDIF
718 #endif /* ALLOW_MY82 */
719
720 #ifdef ALLOW_GGL90
721 #ifdef ALLOW_AUTODIFF_TAMC
722 CADJ STORE GGL90TKE (:,:,:,bi,bj) = comlev1_bibj, key=itdkey,
723 CADJ & kind = isbyte
724 #endif /* ALLOW_AUTODIFF_TAMC */
725 C-- Compute GGL90 mixing coefficients
726 IF (useGGL90) THEN
727 #ifdef ALLOW_DEBUG
728 IF (debugMode) CALL DEBUG_CALL('GGL90_CALC',myThid)
729 #endif
730 CALL TIMER_START('GGL90_CALC [DO_OCEANIC_PHYS]', myThid)
731 CALL GGL90_CALC(
732 I bi, bj, myTime, myIter, myThid )
733 CALL TIMER_STOP ('GGL90_CALC [DO_OCEANIC_PHYS]', myThid)
734 ENDIF
735 #endif /* ALLOW_GGL90 */
736
737 #ifdef ALLOW_TIMEAVE
738 IF ( taveFreq.GT. 0. _d 0 ) THEN
739 CALL TIMEAVE_SURF_FLUX( bi, bj, myTime, myIter, myThid)
740 ENDIF
741 IF (taveFreq.GT.0. .AND. ivdc_kappa.NE.0.) THEN
742 CALL TIMEAVE_CUMULATE(ConvectCountTave, IVDConvCount,
743 I Nr, deltaTclock, bi, bj, myThid)
744 ENDIF
745 #endif /* ALLOW_TIMEAVE */
746
747 #ifdef ALLOW_GMREDI
748 #ifdef ALLOW_AUTODIFF_TAMC
749 # ifndef GM_EXCLUDE_CLIPPING
750 cph storing here is needed only for one GMREDI_OPTIONS:
751 cph define GM_BOLUS_ADVEC
752 cph keep it although TAF says you dont need to.
753 cph but I have avoided the #ifdef for now, in case more things change
754 CADJ STORE sigmaX(:,:,:) = comlev1_bibj, key=itdkey,
755 CADJ & kind = isbyte
756 CADJ STORE sigmaY(:,:,:) = comlev1_bibj, key=itdkey,
757 CADJ & kind = isbyte
758 CADJ STORE sigmaR(:,:,:) = comlev1_bibj, key=itdkey,
759 CADJ & kind = isbyte
760 # endif
761 #endif /* ALLOW_AUTODIFF_TAMC */
762
763 C-- Calculate iso-neutral slopes for the GM/Redi parameterisation
764 IF (useGMRedi) THEN
765 #ifdef ALLOW_DEBUG
766 IF (debugMode) CALL DEBUG_CALL('GMREDI_CALC_TENSOR',myThid)
767 #endif
768 CALL GMREDI_CALC_TENSOR(
769 I iMin, iMax, jMin, jMax,
770 I sigmaX, sigmaY, sigmaR,
771 I bi, bj, myTime, myIter, myThid )
772 #ifdef ALLOW_AUTODIFF_TAMC
773 ELSE
774 CALL GMREDI_CALC_TENSOR_DUMMY(
775 I iMin, iMax, jMin, jMax,
776 I sigmaX, sigmaY, sigmaR,
777 I bi, bj, myTime, myIter, myThid )
778 #endif /* ALLOW_AUTODIFF_TAMC */
779 ENDIF
780 #endif /* ALLOW_GMREDI */
781
782 #ifdef ALLOW_DOWN_SLOPE
783 IF ( useDOWN_SLOPE ) THEN
784 C-- Calculate Downsloping Flow for Down_Slope parameterization
785 IF ( usingPCoords ) THEN
786 CALL DWNSLP_CALC_FLOW(
787 I bi, bj, kSurfC, rhoInSitu,
788 I myTime, myIter, myThid )
789 ELSE
790 CALL DWNSLP_CALC_FLOW(
791 I bi, bj, kLowC, rhoInSitu,
792 I myTime, myIter, myThid )
793 ENDIF
794 ENDIF
795 #endif /* ALLOW_DOWN_SLOPE */
796
797 #ifndef ALLOW_AUTODIFF_TAMC
798 C--- if fluid Is Water: end
799 ENDIF
800 #endif
801
802 #ifdef ALLOW_MYPACKAGE
803 IF ( useMYPACKAGE ) THEN
804 CALL MYPACKAGE_CALC_RHS(
805 I bi, bj, myTime, myIter, myThid )
806 ENDIF
807 #endif /* ALLOW_MYPACKAGE */
808
809 C-- end bi,bj loops.
810 ENDDO
811 ENDDO
812
813 #ifdef ALLOW_GMREDI
814 IF ( useGMRedi ) THEN
815 CALL GMREDI_DO_EXCH( myTime, myIter, myThid )
816 ENDIF
817 #endif /* ALLOW_GMREDI */
818
819 #ifdef ALLOW_KPP
820 IF (useKPP) THEN
821 CALL KPP_DO_EXCH( myThid )
822 ENDIF
823 #endif /* ALLOW_KPP */
824
825 #ifdef ALLOW_DIAGNOSTICS
826 IF ( fluidIsWater .AND. useDiagnostics ) THEN
827 CALL DIAGS_RHO_G(
828 I rhoInSitu, uVel, vVel,
829 I myTime, myIter, myThid )
830 CALL DIAGS_OCEANIC_SURF_FLUX( myTime, myIter, myThid )
831 ENDIF
832 IF ( ivdc_kappa.NE.0 .AND. useDiagnostics ) THEN
833 CALL DIAGNOSTICS_FILL( IVDConvCount, 'CONVADJ ',
834 & 0, Nr, 0, 1, 1, myThid )
835 ENDIF
836 #endif
837
838 #ifdef ALLOW_DEBUG
839 IF (debugMode) CALL DEBUG_LEAVE('DO_OCEANIC_PHYS',myThid)
840 #endif
841
842 RETURN
843 END

  ViewVC Help
Powered by ViewVC 1.1.22