/[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.133 - (show annotations) (download)
Tue Sep 24 01:42:29 2013 UTC (10 years, 8 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint64q, checkpoint64p, checkpoint64s, checkpoint64r, checkpoint64u, checkpoint64t, checkpoint64o
Changes since 1.132: +3 -1 lines
also fill-up some (selected) surface diags if fluidIsAir

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

  ViewVC Help
Powered by ViewVC 1.1.22