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