/[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.16 - (show annotations) (download)
Tue Dec 14 04:36:33 2004 UTC (19 years, 5 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint57d_post, checkpoint57b_post, checkpoint57c_pre, checkpoint57c_post
Changes since 1.15: +4 -4 lines
use new S/R diagnostics_fill to fill-in diagnostics.

1 C $Header: /u/gcmpack/MITgcm/model/src/do_oceanic_phys.F,v 1.15 2004/12/10 20:15:10 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 #endif /* ALLOW_AUTODIFF_TAMC */
15
16 CBOP
17 C !ROUTINE: DO_OCEANIC_PHYS
18 C !INTERFACE:
19 SUBROUTINE DO_OCEANIC_PHYS(myTime, myIter, myThid)
20 C !DESCRIPTION: \bv
21 C *==========================================================*
22 C | SUBROUTINE DO_OCEANIC_PHYS
23 C | o Controlling routine for oceanic physics and
24 C | parameterization
25 C *==========================================================*
26 C | o originally, part of S/R thermodynamics
27 C *==========================================================*
28 C \ev
29
30 C !USES:
31 IMPLICIT NONE
32 C == Global variables ===
33 #include "SIZE.h"
34 #include "EEPARAMS.h"
35 #include "PARAMS.h"
36 #include "DYNVARS.h"
37 #include "GRID.h"
38 c #include "GAD.h"
39 c #ifdef ALLOW_PTRACERS
40 c #include "PTRACERS_SIZE.h"
41 c #include "PTRACERS.h"
42 c #endif
43 c #ifdef ALLOW_TIMEAVE
44 c #include "TIMEAVE_STATV.h"
45 c #endif
46
47 #ifdef ALLOW_AUTODIFF_TAMC
48 # include "tamc.h"
49 # include "tamc_keys.h"
50 # include "FFIELDS.h"
51 # include "EOS.h"
52 # ifdef ALLOW_KPP
53 # include "KPP.h"
54 # endif
55 # ifdef ALLOW_GMREDI
56 # include "GMREDI.h"
57 # endif
58 # ifdef ALLOW_EBM
59 # include "EBM.h"
60 # endif
61 # ifdef EXACT_CONSERV
62 # include "SURFACE.h"
63 # endif
64 #endif /* ALLOW_AUTODIFF_TAMC */
65
66 C !INPUT/OUTPUT PARAMETERS:
67 C == Routine arguments ==
68 C myTime - Current time in simulation
69 C myIter - Current iteration number in simulation
70 C myThid - Thread number for this instance of the routine.
71 _RL myTime
72 INTEGER myIter
73 INTEGER myThid
74
75 C !LOCAL VARIABLES:
76 C == Local variables
77 C rhoK, rhoKM1 - Density at current level, and level above
78 C useVariableK = T when vertical diffusion is not constant
79 C iMin, iMax - Ranges and sub-block indices on which calculations
80 C jMin, jMax are applied.
81 C bi, bj
82 C k, kup, - Index for layer above and below. kup and kDown
83 C kDown, km1 are switched with layer to be the appropriate
84 C index into fVerTerm.
85 _RL rhokm1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
86 _RL rhok (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
87 _RL sigmaX (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
88 _RL sigmaY (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
89 _RL sigmaR (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
90 _RL kp1Msk
91 LOGICAL useVariableK
92 INTEGER iMin, iMax
93 INTEGER jMin, jMax
94 INTEGER bi, bj
95 INTEGER i, j
96 INTEGER k, km1, kup, kDown
97 INTEGER iTracer, ip
98
99 CEOP
100
101 #ifdef ALLOW_AUTODIFF_TAMC
102 C-- dummy statement to end declaration part
103 itdkey = 1
104 #endif /* ALLOW_AUTODIFF_TAMC */
105
106 #ifdef ALLOW_DEBUG
107 IF ( debugLevel .GE. debLevB )
108 & CALL DEBUG_ENTER('DO_OCEANIC_PHYS',myThid)
109 #endif
110
111 #ifdef ALLOW_THSICE
112 IF ( useThSIce .AND. fluidIsWater ) THEN
113 #ifdef ALLOW_DEBUG
114 IF ( debugLevel .GE. debLevB )
115 & CALL DEBUG_CALL('THSICE_MAIN',myThid)
116 #endif
117 C-- Step forward Therm.Sea-Ice variables
118 C and modify forcing terms including effects from ice
119 CALL TIMER_START('THSICE_MAIN [DO_OCEANIC_PHYS]', myThid)
120 CALL THSICE_MAIN( myTime, myIter, myThid )
121 CALL TIMER_STOP( 'THSICE_MAIN [DO_OCEANIC_PHYS]', myThid)
122 ENDIF
123 #endif /* ALLOW_THSICE */
124
125 C-- Freeze water at the surface
126 #ifdef ALLOW_AUTODIFF_TAMC
127 CADJ STORE theta = comlev1, key = ikey_dynamics
128 #endif
129 IF ( allowFreezing
130 & .AND. .NOT. useSEAICE
131 & .AND. .NOT. useThSIce ) THEN
132 CALL FREEZE_SURFACE( myTime, myIter, myThid )
133 ENDIF
134
135 #ifdef COMPONENT_MODULE
136 # ifndef ALLOW_AIM
137 C-- Apply imported data (from coupled interface) to forcing fields
138 C jmc: do not know precisely where to put this call (bf or af thSIce ?)
139 IF ( useCoupler ) THEN
140 CALL OCN_APPLY_IMPORT( .TRUE., myTime, myIter, myThid )
141 ENDIF
142 # endif
143 #endif /* COMPONENT_MODULE */
144
145 #ifdef ALLOW_AUTODIFF_TAMC
146 C-- HPF directive to help TAMC
147 CHPF$ INDEPENDENT
148 #endif /* ALLOW_AUTODIFF_TAMC */
149 DO bj=myByLo(myThid),myByHi(myThid)
150 #ifdef ALLOW_AUTODIFF_TAMC
151 C-- HPF directive to help TAMC
152 CHPF$ INDEPENDENT
153 #endif /* ALLOW_AUTODIFF_TAMC */
154 DO bi=myBxLo(myThid),myBxHi(myThid)
155
156 #ifdef ALLOW_AUTODIFF_TAMC
157 act1 = bi - myBxLo(myThid)
158 max1 = myBxHi(myThid) - myBxLo(myThid) + 1
159 act2 = bj - myByLo(myThid)
160 max2 = myByHi(myThid) - myByLo(myThid) + 1
161 act3 = myThid - 1
162 max3 = nTx*nTy
163 act4 = ikey_dynamics - 1
164 itdkey = (act1 + 1) + act2*max1
165 & + act3*max1*max2
166 & + act4*max1*max2*max3
167 #endif /* ALLOW_AUTODIFF_TAMC */
168
169 C-- Set up work arrays with valid (i.e. not NaN) values
170 C These inital values do not alter the numerical results. They
171 C just ensure that all memory references are to valid floating
172 C point numbers. This prevents spurious hardware signals due to
173 C uninitialised but inert locations.
174
175 DO j=1-OLy,sNy+OLy
176 DO i=1-OLx,sNx+OLx
177 rhok (i,j) = 0. _d 0
178 rhoKM1 (i,j) = 0. _d 0
179 ENDDO
180 ENDDO
181
182 DO k=1,Nr
183 DO j=1-OLy,sNy+OLy
184 DO i=1-OLx,sNx+OLx
185 C This is currently also used by IVDC and Diagnostics
186 sigmaX(i,j,k) = 0. _d 0
187 sigmaY(i,j,k) = 0. _d 0
188 sigmaR(i,j,k) = 0. _d 0
189 #ifdef ALLOW_AUTODIFF_TAMC
190 cph all the following init. are necessary for TAF
191 cph although some of these are re-initialised later.
192 IVDConvCount(i,j,k,bi,bj) = 0.
193 # ifdef ALLOW_GMREDI
194 Kwx(i,j,k,bi,bj) = 0. _d 0
195 Kwy(i,j,k,bi,bj) = 0. _d 0
196 Kwz(i,j,k,bi,bj) = 0. _d 0
197 # ifdef GM_NON_UNITY_DIAGONAL
198 Kux(i,j,k,bi,bj) = 0. _d 0
199 Kvy(i,j,k,bi,bj) = 0. _d 0
200 # endif
201 # ifdef GM_EXTRA_DIAGONAL
202 Kuz(i,j,k,bi,bj) = 0. _d 0
203 Kvz(i,j,k,bi,bj) = 0. _d 0
204 # endif
205 # ifdef GM_BOLUS_ADVEC
206 GM_PsiX(i,j,k,bi,bj) = 0. _d 0
207 GM_PsiY(i,j,k,bi,bj) = 0. _d 0
208 # endif
209 # ifdef GM_VISBECK_VARIABLE_K
210 VisbeckK(i,j,bi,bj) = 0. _d 0
211 # endif
212 # endif /* ALLOW_GMREDI */
213 #endif /* ALLOW_AUTODIFF_TAMC */
214 ENDDO
215 ENDDO
216 ENDDO
217
218 iMin = 1-OLx
219 iMax = sNx+OLx
220 jMin = 1-OLy
221 jMax = sNy+OLy
222
223 #ifdef ALLOW_AUTODIFF_TAMC
224 CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
225 CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
226 CADJ STORE totphihyd(:,:,:,bi,bj)
227 CADJ & = comlev1_bibj, key=itdkey, byte=isbyte
228 # ifdef ALLOW_KPP
229 CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
230 CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
231 # endif
232 # ifdef EXACT_CONSERV
233 CADJ STORE pmepr(:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
234 # endif
235 #endif /* ALLOW_AUTODIFF_TAMC */
236
237 #ifdef ALLOW_DEBUG
238 IF ( debugLevel .GE. debLevB )
239 & CALL DEBUG_MSG('ENTERING UPWARD K LOOP',myThid)
240 #endif
241
242 C-- Start of diagnostic loop
243 DO k=Nr,1,-1
244
245 #ifdef ALLOW_AUTODIFF_TAMC
246 C? Patrick, is this formula correct now that we change the loop range?
247 C? Do we still need this?
248 cph kkey formula corrected.
249 cph Needed for rhok, rhokm1, in the case useGMREDI.
250 kkey = (itdkey-1)*Nr + k
251 #endif /* ALLOW_AUTODIFF_TAMC */
252
253 C-- Calculate gradients of potential density for isoneutral
254 C slope terms (e.g. GM/Redi tensor or IVDC diffusivity)
255 c IF ( k.GT.1 .AND. (useGMRedi.OR.ivdc_kappa.NE.0.) ) THEN
256 IF ( useGMRedi .OR. (k.GT.1 .AND. ivdc_kappa.NE.0.) ) THEN
257 #ifdef ALLOW_DEBUG
258 IF ( debugLevel .GE. debLevB )
259 & CALL DEBUG_CALL('FIND_RHO',myThid)
260 #endif
261 #ifdef ALLOW_AUTODIFF_TAMC
262 CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
263 CADJ STORE salt (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
264 #endif /* ALLOW_AUTODIFF_TAMC */
265 CALL FIND_RHO(
266 I bi, bj, iMin, iMax, jMin, jMax, k, k,
267 I theta, salt,
268 O rhoK,
269 I myThid )
270
271 IF (k.GT.1) THEN
272 #ifdef ALLOW_AUTODIFF_TAMC
273 CADJ STORE theta(:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
274 CADJ STORE salt (:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
275 #endif /* ALLOW_AUTODIFF_TAMC */
276 CALL FIND_RHO(
277 I bi, bj, iMin, iMax, jMin, jMax, k-1, k,
278 I theta, salt,
279 O rhoKm1,
280 I myThid )
281 ENDIF
282 #ifdef ALLOW_DEBUG
283 IF ( debugLevel .GE. debLevB )
284 & CALL DEBUG_CALL('GRAD_SIGMA',myThid)
285 #endif
286 CALL GRAD_SIGMA(
287 I bi, bj, iMin, iMax, jMin, jMax, k,
288 I rhoK, rhoKm1, rhoK,
289 O sigmaX, sigmaY, sigmaR,
290 I myThid )
291 ENDIF
292
293 #ifdef ALLOW_AUTODIFF_TAMC
294 ctest# ifndef GM_EXCLUDE_CLIPPING
295 CADJ STORE rhok (:,:) = comlev1_bibj_k , key=kkey, byte=isbyte
296 ctest# endif
297 CADJ STORE rhokm1 (:,:) = comlev1_bibj_k , key=kkey, byte=isbyte
298 #endif /* ALLOW_AUTODIFF_TAMC */
299 C-- Implicit Vertical Diffusion for Convection
300 c ==> should use sigmaR !!!
301 IF (k.GT.1 .AND. ivdc_kappa.NE.0.) THEN
302 #ifdef ALLOW_DEBUG
303 IF ( debugLevel .GE. debLevB )
304 & CALL DEBUG_CALL('CALC_IVDC',myThid)
305 #endif
306 CALL CALC_IVDC(
307 I bi, bj, iMin, iMax, jMin, jMax, k,
308 I rhoKm1, rhoK,
309 I myTime, myIter, myThid)
310 ENDIF
311
312 C-- end of diagnostic k loop (Nr:1)
313 ENDDO
314
315 #ifdef ALLOW_DIAGNOSTICS
316 IF ( useDiagnostics .AND.
317 & (useGMRedi .OR. ivdc_kappa.NE.0.) ) THEN
318 CALL DIAGNOSTICS_FILL (sigmaR, 'DRHODR ', 0, Nr,
319 & 2, bi, bj, myThid)
320 ENDIF
321 #endif
322
323 #ifdef ALLOW_OBCS
324 C-- Calculate future values on open boundaries
325 IF (useOBCS) THEN
326 #ifdef ALLOW_DEBUG
327 IF ( debugLevel .GE. debLevB )
328 & CALL DEBUG_CALL('OBCS_CALC',myThid)
329 #endif
330 CALL OBCS_CALC( bi, bj, myTime+deltaTclock, myIter+1,
331 I uVel, vVel, wVel, theta, salt,
332 I myThid )
333 ENDIF
334 #endif /* ALLOW_OBCS */
335
336 #ifndef ALLOW_AUTODIFF_TAMC
337 IF ( fluidIsWater ) THEN
338 #endif
339 C-- Determines forcing terms based on external fields
340 C relaxation terms, etc.
341 #ifdef ALLOW_DEBUG
342 IF ( debugLevel .GE. debLevB )
343 & CALL DEBUG_CALL('EXTERNAL_FORCING_SURF',myThid)
344 #endif
345 CALL EXTERNAL_FORCING_SURF(
346 I bi, bj, iMin, iMax, jMin, jMax,
347 I myTime, myIter, myThid )
348 #ifndef ALLOW_AUTODIFF_TAMC
349 ENDIF
350 #endif
351
352 #ifdef ALLOW_AUTODIFF_TAMC
353 cph needed for KPP
354 CADJ STORE surfaceForcingU(:,:,bi,bj)
355 CADJ & = comlev1_bibj, key=itdkey, byte=isbyte
356 CADJ STORE surfaceForcingV(:,:,bi,bj)
357 CADJ & = comlev1_bibj, key=itdkey, byte=isbyte
358 CADJ STORE surfaceForcingS(:,:,bi,bj)
359 CADJ & = comlev1_bibj, key=itdkey, byte=isbyte
360 CADJ STORE surfaceForcingT(:,:,bi,bj)
361 CADJ & = comlev1_bibj, key=itdkey, byte=isbyte
362 CADJ STORE surfaceForcingTice(:,:,bi,bj)
363 CADJ & = comlev1_bibj, key=itdkey, byte=isbyte
364
365 #endif /* ALLOW_AUTODIFF_TAMC */
366
367 #ifdef ALLOW_GMREDI
368
369 #ifdef ALLOW_AUTODIFF_TAMC
370 # ifndef GM_EXCLUDE_CLIPPING
371 cph storing here is needed only for one GMREDI_OPTIONS:
372 cph define GM_BOLUS_ADVEC
373 cph keep it although TAF says you dont need to.
374 cph but I've avoided the #ifdef for now, in case more things change
375 CADJ STORE sigmaX(:,:,:) = comlev1_bibj, key=itdkey, byte=isbyte
376 CADJ STORE sigmaY(:,:,:) = comlev1_bibj, key=itdkey, byte=isbyte
377 CADJ STORE sigmaR(:,:,:) = comlev1_bibj, key=itdkey, byte=isbyte
378 # endif
379 #endif /* ALLOW_AUTODIFF_TAMC */
380
381 C-- Calculate iso-neutral slopes for the GM/Redi parameterisation
382 IF (useGMRedi) THEN
383 #ifdef ALLOW_DEBUG
384 IF ( debugLevel .GE. debLevB )
385 & CALL DEBUG_CALL('GMREDI_CALC_TENSOR',myThid)
386 #endif
387 CALL GMREDI_CALC_TENSOR(
388 I bi, bj, iMin, iMax, jMin, jMax,
389 I sigmaX, sigmaY, sigmaR,
390 I myThid )
391 #ifdef ALLOW_AUTODIFF_TAMC
392 ELSE
393 CALL GMREDI_CALC_TENSOR_DUMMY(
394 I bi, bj, iMin, iMax, jMin, jMax,
395 I sigmaX, sigmaY, sigmaR,
396 I myThid )
397 #endif /* ALLOW_AUTODIFF_TAMC */
398 ENDIF
399
400 #endif /* ALLOW_GMREDI */
401
402 #ifdef ALLOW_KPP
403 C-- Compute KPP mixing coefficients
404 IF (useKPP) THEN
405 #ifdef ALLOW_DEBUG
406 IF ( debugLevel .GE. debLevB )
407 & CALL DEBUG_CALL('KPP_CALC',myThid)
408 #endif
409 CALL KPP_CALC(
410 I bi, bj, myTime, myThid )
411 #ifdef ALLOW_AUTODIFF_TAMC
412 ELSE
413 CALL KPP_CALC_DUMMY(
414 I bi, bj, myTime, myThid )
415 #endif /* ALLOW_AUTODIFF_TAMC */
416 ENDIF
417
418 #endif /* ALLOW_KPP */
419
420 #ifdef ALLOW_PP81
421 C-- Compute PP81 mixing coefficients
422 IF (usePP81) THEN
423 #ifdef ALLOW_DEBUG
424 IF ( debugLevel .GE. debLevB )
425 & CALL DEBUG_CALL('PP81_CALC',myThid)
426 #endif
427 CALL PP81_CALC(
428 I bi, bj, myTime, myThid )
429 ENDIF
430 #endif /* ALLOW_PP81 */
431
432 #ifdef ALLOW_MY82
433 C-- Compute MY82 mixing coefficients
434 IF (useMY82) THEN
435 #ifdef ALLOW_DEBUG
436 IF ( debugLevel .GE. debLevB )
437 & CALL DEBUG_CALL('MY82_CALC',myThid)
438 #endif
439 CALL MY82_CALC(
440 I bi, bj, myTime, myThid )
441 ENDIF
442 #endif /* ALLOW_MY82 */
443
444 #ifdef ALLOW_GGL90
445 C-- Compute GGL90 mixing coefficients
446 IF (useGGL90) THEN
447 #ifdef ALLOW_DEBUG
448 IF ( debugLevel .GE. debLevB )
449 & CALL DEBUG_CALL('GGL90_CALC',myThid)
450 #endif
451 CALL GGL90_CALC(
452 I bi, bj, myTime, myThid )
453 ENDIF
454 #endif /* ALLOW_GGL90 */
455
456 C-- end bi,bj loops.
457 ENDDO
458 ENDDO
459
460 #ifdef ALLOW_DEBUG
461 IF ( debugLevel .GE. debLevB )
462 & CALL DEBUG_LEAVE('DO_OCEANIC_PHYS',myThid)
463 #endif
464
465 RETURN
466 END

  ViewVC Help
Powered by ViewVC 1.1.22