/[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.14 - (show annotations) (download)
Tue Oct 19 02:39:58 2004 UTC (19 years, 7 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint55j_post, checkpoint56b_post, checkpoint55h_post, checkpoint56c_post, checkpoint57a_pre, checkpoint55i_post, checkpoint57, checkpoint56, checkpoint56a_post
Changes since 1.13: +3 -3 lines
use flags: fluidIsAir/Water, usingP/ZCoords instead of buoyancyRelation

1 C $Header: /u/gcmpack/MITgcm/model/src/do_oceanic_phys.F,v 1.13 2004/10/14 05:22:21 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
150 DO bj=myByLo(myThid),myByHi(myThid)
151
152 #ifdef ALLOW_AUTODIFF_TAMC
153 C-- HPF directive to help TAMC
154 CHPF$ INDEPENDENT, NEW (rTrans,fVerT,fVerS
155 CHPF$& ,utrans,vtrans,xA,yA
156 CHPF$& )
157 #endif /* ALLOW_AUTODIFF_TAMC */
158
159 DO bi=myBxLo(myThid),myBxHi(myThid)
160
161 #ifdef ALLOW_AUTODIFF_TAMC
162 act1 = bi - myBxLo(myThid)
163 max1 = myBxHi(myThid) - myBxLo(myThid) + 1
164 act2 = bj - myByLo(myThid)
165 max2 = myByHi(myThid) - myByLo(myThid) + 1
166 act3 = myThid - 1
167 max3 = nTx*nTy
168 act4 = ikey_dynamics - 1
169 itdkey = (act1 + 1) + act2*max1
170 & + act3*max1*max2
171 & + act4*max1*max2*max3
172 #endif /* ALLOW_AUTODIFF_TAMC */
173
174 C-- Set up work arrays with valid (i.e. not NaN) values
175 C These inital values do not alter the numerical results. They
176 C just ensure that all memory references are to valid floating
177 C point numbers. This prevents spurious hardware signals due to
178 C uninitialised but inert locations.
179
180 DO j=1-OLy,sNy+OLy
181 DO i=1-OLx,sNx+OLx
182 rhok (i,j) = 0. _d 0
183 rhoKM1 (i,j) = 0. _d 0
184 ENDDO
185 ENDDO
186
187 DO k=1,Nr
188 DO j=1-OLy,sNy+OLy
189 DO i=1-OLx,sNx+OLx
190 C This is currently also used by IVDC and Diagnostics
191 sigmaX(i,j,k) = 0. _d 0
192 sigmaY(i,j,k) = 0. _d 0
193 sigmaR(i,j,k) = 0. _d 0
194 #ifdef ALLOW_AUTODIFF_TAMC
195 cph all the following init. are necessary for TAF
196 cph although some of these are re-initialised later.
197 IVDConvCount(i,j,k,bi,bj) = 0.
198 # ifdef ALLOW_GMREDI
199 Kwx(i,j,k,bi,bj) = 0. _d 0
200 Kwy(i,j,k,bi,bj) = 0. _d 0
201 Kwz(i,j,k,bi,bj) = 0. _d 0
202 # ifdef GM_NON_UNITY_DIAGONAL
203 Kux(i,j,k,bi,bj) = 0. _d 0
204 Kvy(i,j,k,bi,bj) = 0. _d 0
205 # endif
206 # ifdef GM_EXTRA_DIAGONAL
207 Kuz(i,j,k,bi,bj) = 0. _d 0
208 Kvz(i,j,k,bi,bj) = 0. _d 0
209 # endif
210 # ifdef GM_BOLUS_ADVEC
211 GM_PsiX(i,j,k,bi,bj) = 0. _d 0
212 GM_PsiY(i,j,k,bi,bj) = 0. _d 0
213 # endif
214 # ifdef GM_VISBECK_VARIABLE_K
215 VisbeckK(i,j,bi,bj) = 0. _d 0
216 # endif
217 # endif /* ALLOW_GMREDI */
218 #endif /* ALLOW_AUTODIFF_TAMC */
219 ENDDO
220 ENDDO
221 ENDDO
222
223 iMin = 1-OLx
224 iMax = sNx+OLx
225 jMin = 1-OLy
226 jMax = sNy+OLy
227
228 #ifdef ALLOW_AUTODIFF_TAMC
229 CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
230 CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
231 CADJ STORE totphihyd(:,:,:,bi,bj)
232 CADJ & = comlev1_bibj, key=itdkey, byte=isbyte
233 # ifdef ALLOW_KPP
234 CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
235 CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
236 # endif
237 # ifdef EXACT_CONSERV
238 CADJ STORE pmepr(:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
239 # endif
240 #endif /* ALLOW_AUTODIFF_TAMC */
241
242 #ifdef ALLOW_DEBUG
243 IF ( debugLevel .GE. debLevB )
244 & CALL DEBUG_MSG('ENTERING UPWARD K LOOP',myThid)
245 #endif
246
247 C-- Start of diagnostic loop
248 DO k=Nr,1,-1
249
250 #ifdef ALLOW_AUTODIFF_TAMC
251 C? Patrick, is this formula correct now that we change the loop range?
252 C? Do we still need this?
253 cph kkey formula corrected.
254 cph Needed for rhok, rhokm1, in the case useGMREDI.
255 kkey = (itdkey-1)*Nr + k
256 #endif /* ALLOW_AUTODIFF_TAMC */
257
258 C-- Calculate gradients of potential density for isoneutral
259 C slope terms (e.g. GM/Redi tensor or IVDC diffusivity)
260 c IF ( k.GT.1 .AND. (useGMRedi.OR.ivdc_kappa.NE.0.) ) THEN
261 IF ( useGMRedi .OR. (k.GT.1 .AND. ivdc_kappa.NE.0.) ) THEN
262 #ifdef ALLOW_DEBUG
263 IF ( debugLevel .GE. debLevB )
264 & CALL DEBUG_CALL('FIND_RHO',myThid)
265 #endif
266 #ifdef ALLOW_AUTODIFF_TAMC
267 CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
268 CADJ STORE salt (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
269 #endif /* ALLOW_AUTODIFF_TAMC */
270 CALL FIND_RHO(
271 I bi, bj, iMin, iMax, jMin, jMax, k, k,
272 I theta, salt,
273 O rhoK,
274 I myThid )
275
276 IF (k.GT.1) THEN
277 #ifdef ALLOW_AUTODIFF_TAMC
278 CADJ STORE theta(:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
279 CADJ STORE salt (:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
280 #endif /* ALLOW_AUTODIFF_TAMC */
281 CALL FIND_RHO(
282 I bi, bj, iMin, iMax, jMin, jMax, k-1, k,
283 I theta, salt,
284 O rhoKm1,
285 I myThid )
286 ENDIF
287 #ifdef ALLOW_DEBUG
288 IF ( debugLevel .GE. debLevB )
289 & CALL DEBUG_CALL('GRAD_SIGMA',myThid)
290 #endif
291 CALL GRAD_SIGMA(
292 I bi, bj, iMin, iMax, jMin, jMax, k,
293 I rhoK, rhoKm1, rhoK,
294 O sigmaX, sigmaY, sigmaR,
295 I myThid )
296 ENDIF
297
298 #ifdef ALLOW_AUTODIFF_TAMC
299 ctest# ifndef GM_EXCLUDE_CLIPPING
300 CADJ STORE rhok (:,:) = comlev1_bibj_k , key=kkey, byte=isbyte
301 ctest# endif
302 CADJ STORE rhokm1 (:,:) = comlev1_bibj_k , key=kkey, byte=isbyte
303 #endif /* ALLOW_AUTODIFF_TAMC */
304 C-- Implicit Vertical Diffusion for Convection
305 c ==> should use sigmaR !!!
306 IF (k.GT.1 .AND. ivdc_kappa.NE.0.) THEN
307 #ifdef ALLOW_DEBUG
308 IF ( debugLevel .GE. debLevB )
309 & CALL DEBUG_CALL('CALC_IVDC',myThid)
310 #endif
311 CALL CALC_IVDC(
312 I bi, bj, iMin, iMax, jMin, jMax, k,
313 I rhoKm1, rhoK,
314 I myTime, myIter, myThid)
315 ENDIF
316
317 C-- end of diagnostic k loop (Nr:1)
318 ENDDO
319
320 #ifdef ALLOW_DIAGNOSTICS
321 IF ( usediagnostics .AND.
322 & (useGMRedi .OR. ivdc_kappa.NE.0.) ) THEN
323 CALL fill_diagnostics (myThid, 'DRHODR ', 0, Nr,
324 & 3, bi, bj, sigmaR)
325 ENDIF
326 #endif
327
328 #ifdef ALLOW_OBCS
329 C-- Calculate future values on open boundaries
330 IF (useOBCS) THEN
331 #ifdef ALLOW_DEBUG
332 IF ( debugLevel .GE. debLevB )
333 & CALL DEBUG_CALL('OBCS_CALC',myThid)
334 #endif
335 CALL OBCS_CALC( bi, bj, myTime+deltaTclock, myIter+1,
336 I uVel, vVel, wVel, theta, salt,
337 I myThid )
338 ENDIF
339 #endif /* ALLOW_OBCS */
340
341 #ifndef ALLOW_AUTODIFF_TAMC
342 IF ( fluidIsWater ) THEN
343 #endif
344 C-- Determines forcing terms based on external fields
345 C relaxation terms, etc.
346 #ifdef ALLOW_DEBUG
347 IF ( debugLevel .GE. debLevB )
348 & CALL DEBUG_CALL('EXTERNAL_FORCING_SURF',myThid)
349 #endif
350 CALL EXTERNAL_FORCING_SURF(
351 I bi, bj, iMin, iMax, jMin, jMax,
352 I myTime, myIter, myThid )
353 #ifndef ALLOW_AUTODIFF_TAMC
354 ENDIF
355 #endif
356
357 #ifdef ALLOW_AUTODIFF_TAMC
358 cph needed for KPP
359 CADJ STORE surfaceForcingU(:,:,bi,bj)
360 CADJ & = comlev1_bibj, key=itdkey, byte=isbyte
361 CADJ STORE surfaceForcingV(:,:,bi,bj)
362 CADJ & = comlev1_bibj, key=itdkey, byte=isbyte
363 CADJ STORE surfaceForcingS(:,:,bi,bj)
364 CADJ & = comlev1_bibj, key=itdkey, byte=isbyte
365 CADJ STORE surfaceForcingT(:,:,bi,bj)
366 CADJ & = comlev1_bibj, key=itdkey, byte=isbyte
367 CADJ STORE surfaceForcingTice(:,:,bi,bj)
368 CADJ & = comlev1_bibj, key=itdkey, byte=isbyte
369
370 #endif /* ALLOW_AUTODIFF_TAMC */
371
372 #ifdef ALLOW_GMREDI
373
374 #ifdef ALLOW_AUTODIFF_TAMC
375 # ifndef GM_EXCLUDE_CLIPPING
376 cph storing here is needed only for one GMREDI_OPTIONS:
377 cph define GM_BOLUS_ADVEC
378 cph keep it although TAF says you dont need to.
379 cph but I've avoided the #ifdef for now, in case more things change
380 CADJ STORE sigmaX(:,:,:) = comlev1_bibj, key=itdkey, byte=isbyte
381 CADJ STORE sigmaY(:,:,:) = comlev1_bibj, key=itdkey, byte=isbyte
382 CADJ STORE sigmaR(:,:,:) = comlev1_bibj, key=itdkey, byte=isbyte
383 # endif
384 #endif /* ALLOW_AUTODIFF_TAMC */
385
386 C-- Calculate iso-neutral slopes for the GM/Redi parameterisation
387 IF (useGMRedi) THEN
388 #ifdef ALLOW_DEBUG
389 IF ( debugLevel .GE. debLevB )
390 & CALL DEBUG_CALL('GMREDI_CALC_TENSOR',myThid)
391 #endif
392 CALL GMREDI_CALC_TENSOR(
393 I bi, bj, iMin, iMax, jMin, jMax,
394 I sigmaX, sigmaY, sigmaR,
395 I myThid )
396 #ifdef ALLOW_AUTODIFF_TAMC
397 ELSE
398 CALL GMREDI_CALC_TENSOR_DUMMY(
399 I bi, bj, iMin, iMax, jMin, jMax,
400 I sigmaX, sigmaY, sigmaR,
401 I myThid )
402 #endif /* ALLOW_AUTODIFF_TAMC */
403 ENDIF
404
405 #endif /* ALLOW_GMREDI */
406
407 #ifdef ALLOW_KPP
408 C-- Compute KPP mixing coefficients
409 IF (useKPP) THEN
410 #ifdef ALLOW_DEBUG
411 IF ( debugLevel .GE. debLevB )
412 & CALL DEBUG_CALL('KPP_CALC',myThid)
413 #endif
414 CALL KPP_CALC(
415 I bi, bj, myTime, myThid )
416 #ifdef ALLOW_AUTODIFF_TAMC
417 ELSE
418 CALL KPP_CALC_DUMMY(
419 I bi, bj, myTime, myThid )
420 #endif /* ALLOW_AUTODIFF_TAMC */
421 ENDIF
422
423 #endif /* ALLOW_KPP */
424
425 #ifdef ALLOW_PP81
426 C-- Compute PP81 mixing coefficients
427 IF (usePP81) THEN
428 #ifdef ALLOW_DEBUG
429 IF ( debugLevel .GE. debLevB )
430 & CALL DEBUG_CALL('PP81_CALC',myThid)
431 #endif
432 CALL PP81_CALC(
433 I bi, bj, myTime, myThid )
434 ENDIF
435 #endif /* ALLOW_PP81 */
436
437 #ifdef ALLOW_MY82
438 C-- Compute MY82 mixing coefficients
439 IF (useMY82) THEN
440 #ifdef ALLOW_DEBUG
441 IF ( debugLevel .GE. debLevB )
442 & CALL DEBUG_CALL('MY82_CALC',myThid)
443 #endif
444 CALL MY82_CALC(
445 I bi, bj, myTime, myThid )
446 ENDIF
447 #endif /* ALLOW_MY82 */
448
449 #ifdef ALLOW_GGL90
450 C-- Compute GGL90 mixing coefficients
451 IF (useGGL90) THEN
452 #ifdef ALLOW_DEBUG
453 IF ( debugLevel .GE. debLevB )
454 & CALL DEBUG_CALL('GGL90_CALC',myThid)
455 #endif
456 CALL GGL90_CALC(
457 I bi, bj, myTime, myThid )
458 ENDIF
459 #endif /* ALLOW_GGL90 */
460
461 C-- end bi,bj loops.
462 ENDDO
463 ENDDO
464
465 #ifdef ALLOW_DEBUG
466 IF ( debugLevel .GE. debLevB )
467 & CALL DEBUG_LEAVE('DO_OCEANIC_PHYS',myThid)
468 #endif
469
470 RETURN
471 END

  ViewVC Help
Powered by ViewVC 1.1.22