/[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.107 - (show annotations) (download)
Sun Aug 7 07:08:15 2011 UTC (12 years, 9 months ago) by dimitri
Branch: MAIN
Changes since 1.106: +21 -6 lines
o adding package bbl (Bottom Boundary Layer)
  description in MITgcm/pkg/bbl/bbl_description.tex
  example/test experiment in MITgcm_contrib/bbl

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

  ViewVC Help
Powered by ViewVC 1.1.22