/[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.6 - (show annotations) (download)
Thu Sep 2 09:13:49 2004 UTC (19 years, 8 months ago) by mlosch
Branch: MAIN
CVS Tags: checkpoint54e_post
Changes since 1.5: +25 -1 lines
o add calls for two new packages
  - pp81 (Packanowski and Philander, 1981), Richardson number and
    stratification dependent mixing
  - my82 (Mellor and Yamada, 1982) level 2 turbulence closure scheme

1 C $Header: /u/gcmpack/MITgcm/model/src/do_oceanic_phys.F,v 1.5 2004/07/26 20:18:14 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 #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_PASSIVE_TRACER
40 c #include "TR1.h"
41 c #endif
42 c #ifdef ALLOW_PTRACERS
43 c #include "PTRACERS_SIZE.h"
44 c #include "PTRACERS.h"
45 c #endif
46 c #ifdef ALLOW_TIMEAVE
47 c #include "TIMEAVE_STATV.h"
48 c #endif
49
50 #ifdef ALLOW_AUTODIFF_TAMC
51 # include "tamc.h"
52 # include "tamc_keys.h"
53 # include "FFIELDS.h"
54 # include "EOS.h"
55 # ifdef ALLOW_KPP
56 # include "KPP.h"
57 # endif
58 # ifdef ALLOW_GMREDI
59 # include "GMREDI.h"
60 # endif
61 # ifdef ALLOW_EBM
62 # include "EBM.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_DEBUG
102 IF ( debugLevel .GE. debLevB )
103 & CALL DEBUG_ENTER('DO_OCEANIC_PHYS',myThid)
104 #endif
105
106 #ifdef ALLOW_THSICE
107 IF ( useThSIce .AND. buoyancyRelation(1:7) .EQ. 'OCEANIC' ) THEN
108 #ifdef ALLOW_DEBUG
109 IF ( debugLevel .GE. debLevB )
110 & CALL DEBUG_CALL('THSICE_MAIN',myThid)
111 #endif
112 C-- Step forward Therm.Sea-Ice variables
113 C and modify forcing terms including effects from ice
114 CALL TIMER_START('THSICE_MAIN [DO_OCEANIC_PHYS]', myThid)
115 CALL THSICE_MAIN( myTime, myIter, myThid )
116 CALL TIMER_STOP( 'THSICE_MAIN [DO_OCEANIC_PHYS]', myThid)
117 ENDIF
118 #endif /* ALLOW_THSICE */
119
120 C-- Freeze water at the surface
121 #ifdef ALLOW_AUTODIFF_TAMC
122 CADJ STORE theta = comlev1, key = ikey_dynamics
123 #endif
124 IF ( allowFreezing .AND. .NOT. useSEAICE
125 & .AND. .NOT. useThSIce ) THEN
126 CALL FREEZE_SURFACE( myTime, myIter, myThid )
127 ENDIF
128
129 #ifdef COMPONENT_MODULE
130 # ifndef ALLOW_AIM
131 C-- Apply imported data (from coupled interface) to forcing fields
132 C jmc: don't know precisely where to put this call (bf or af thSIce ?)
133 IF ( useCoupler ) THEN
134 CALL OCN_APPLY_IMPORT( .TRUE., myTime, myIter, myThid )
135 ENDIF
136 # endif
137 #endif /* COMPONENT_MODULE */
138
139 #ifdef ALLOW_AUTODIFF_TAMC
140 C-- dummy statement to end declaration part
141 ikey = 1
142 itdkey = 1
143 #endif /* ALLOW_AUTODIFF_TAMC */
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
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 #endif /* ALLOW_AUTODIFF_TAMC */
238
239 #ifdef ALLOW_DEBUG
240 IF ( debugLevel .GE. debLevB )
241 & CALL DEBUG_MSG('ENTERING UPWARD K LOOP',myThid)
242 #endif
243
244 C-- Start of diagnostic loop
245 DO k=Nr,1,-1
246
247 #ifdef ALLOW_AUTODIFF_TAMC
248 C? Patrick, is this formula correct now that we change the loop range?
249 C? Do we still need this?
250 cph kkey formula corrected.
251 cph Needed for rhok, rhokm1, in the case useGMREDI.
252 kkey = (itdkey-1)*Nr + k
253 #endif /* ALLOW_AUTODIFF_TAMC */
254
255 C-- Calculate gradients of potential density for isoneutral
256 C slope terms (e.g. GM/Redi tensor or IVDC diffusivity)
257 c IF ( k.GT.1 .AND. (useGMRedi.OR.ivdc_kappa.NE.0.) ) THEN
258 IF ( useGMRedi .OR. (k.GT.1 .AND. ivdc_kappa.NE.0.) ) THEN
259 #ifdef ALLOW_DEBUG
260 IF ( debugLevel .GE. debLevB )
261 & CALL DEBUG_CALL('FIND_RHO',myThid)
262 #endif
263 #ifdef ALLOW_AUTODIFF_TAMC
264 CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
265 CADJ STORE salt (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
266 #endif /* ALLOW_AUTODIFF_TAMC */
267 CALL FIND_RHO(
268 I bi, bj, iMin, iMax, jMin, jMax, k, k,
269 I theta, salt,
270 O rhoK,
271 I myThid )
272
273 IF (k.GT.1) THEN
274 #ifdef ALLOW_AUTODIFF_TAMC
275 CADJ STORE theta(:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
276 CADJ STORE salt (:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
277 #endif /* ALLOW_AUTODIFF_TAMC */
278 CALL FIND_RHO(
279 I bi, bj, iMin, iMax, jMin, jMax, k-1, k,
280 I theta, salt,
281 O rhoKm1,
282 I myThid )
283 ENDIF
284 #ifdef ALLOW_DEBUG
285 IF ( debugLevel .GE. debLevB )
286 & CALL DEBUG_CALL('GRAD_SIGMA',myThid)
287 #endif
288 CALL GRAD_SIGMA(
289 I bi, bj, iMin, iMax, jMin, jMax, k,
290 I rhoK, rhoKm1, rhoK,
291 O sigmaX, sigmaY, sigmaR,
292 I myThid )
293 ENDIF
294
295 #ifdef ALLOW_AUTODIFF_TAMC
296 CADJ STORE rhok (:,:) = comlev1_bibj_k , key=kkey, byte=isbyte
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_AUTODIFF_TAMC
316 cph avoids recomputation of integrate_for_w
317 CADJ STORE wvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
318 #endif /* ALLOW_AUTODIFF_TAMC */
319
320 #ifdef ALLOW_OBCS
321 C-- Calculate future values on open boundaries
322 IF (useOBCS) THEN
323 #ifdef ALLOW_DEBUG
324 IF ( debugLevel .GE. debLevB )
325 & CALL DEBUG_CALL('OBCS_CALC',myThid)
326 #endif
327 CALL OBCS_CALC( bi, bj, myTime+deltaT, myIter+1,
328 I uVel, vVel, wVel, theta, salt,
329 I myThid )
330 ENDIF
331 #endif /* ALLOW_OBCS */
332
333 #ifndef ALLOW_AUTODIFF_TAMC
334 IF ( buoyancyRelation(1:7) .EQ. 'OCEANIC' ) THEN
335 #endif
336 C-- Determines forcing terms based on external fields
337 C relaxation terms, etc.
338 #ifdef ALLOW_DEBUG
339 IF ( debugLevel .GE. debLevB )
340 & CALL DEBUG_CALL('EXTERNAL_FORCING_SURF',myThid)
341 #endif
342 CALL EXTERNAL_FORCING_SURF(
343 I bi, bj, iMin, iMax, jMin, jMax,
344 I myTime, myIter, myThid )
345 #ifndef ALLOW_AUTODIFF_TAMC
346 ENDIF
347 #endif
348
349 #ifdef ALLOW_AUTODIFF_TAMC
350 cph needed for KPP
351 CADJ STORE surfaceForcingU(:,:,bi,bj)
352 CADJ & = comlev1_bibj, key=itdkey, byte=isbyte
353 CADJ STORE surfaceForcingV(:,:,bi,bj)
354 CADJ & = comlev1_bibj, key=itdkey, byte=isbyte
355 CADJ STORE surfaceForcingS(:,:,bi,bj)
356 CADJ & = comlev1_bibj, key=itdkey, byte=isbyte
357 CADJ STORE surfaceForcingT(:,:,bi,bj)
358 CADJ & = comlev1_bibj, key=itdkey, byte=isbyte
359 # ifdef ALLOW_SEAICE
360 CADJ STORE surfaceForcingTice(:,:,bi,bj)
361 CADJ & = comlev1_bibj, key=itdkey, byte=isbyte
362 # endif
363 #endif /* ALLOW_AUTODIFF_TAMC */
364
365 #ifdef ALLOW_GMREDI
366
367 #ifdef ALLOW_AUTODIFF_TAMC
368 cph storing here is needed only for one GMREDI_OPTIONS:
369 cph define GM_BOLUS_ADVEC
370 cph but I've avoided the #ifdef for now, in case more things change
371 CADJ STORE sigmaX(:,:,:) = comlev1_bibj, key=itdkey, byte=isbyte
372 CADJ STORE sigmaY(:,:,:) = comlev1_bibj, key=itdkey, byte=isbyte
373 CADJ STORE sigmaR(:,:,:) = comlev1_bibj, key=itdkey, byte=isbyte
374 #endif /* ALLOW_AUTODIFF_TAMC */
375
376 C-- Calculate iso-neutral slopes for the GM/Redi parameterisation
377 IF (useGMRedi) THEN
378 #ifdef ALLOW_DEBUG
379 IF ( debugLevel .GE. debLevB )
380 & CALL DEBUG_CALL('GMREDI_CALC_TENSOR',myThid)
381 #endif
382 CALL GMREDI_CALC_TENSOR(
383 I bi, bj, iMin, iMax, jMin, jMax,
384 I sigmaX, sigmaY, sigmaR,
385 I myThid )
386 #ifdef ALLOW_AUTODIFF_TAMC
387 ELSE
388 CALL GMREDI_CALC_TENSOR_DUMMY(
389 I bi, bj, iMin, iMax, jMin, jMax,
390 I sigmaX, sigmaY, sigmaR,
391 I myThid )
392 #endif /* ALLOW_AUTODIFF_TAMC */
393 ENDIF
394
395 #ifdef ALLOW_AUTODIFF_TAMC
396 CADJ STORE Kwx(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
397 CADJ STORE Kwy(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
398 CADJ STORE Kwz(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
399 #endif /* ALLOW_AUTODIFF_TAMC */
400
401 #endif /* ALLOW_GMREDI */
402
403 #ifdef ALLOW_KPP
404 C-- Compute KPP mixing coefficients
405 IF (useKPP) THEN
406 #ifdef ALLOW_DEBUG
407 IF ( debugLevel .GE. debLevB )
408 & CALL DEBUG_CALL('KPP_CALC',myThid)
409 #endif
410 CALL KPP_CALC(
411 I bi, bj, myTime, myThid )
412 #ifdef ALLOW_AUTODIFF_TAMC
413 ELSE
414 CALL KPP_CALC_DUMMY(
415 I bi, bj, myTime, myThid )
416 #endif /* ALLOW_AUTODIFF_TAMC */
417 ENDIF
418
419 #ifdef ALLOW_AUTODIFF_TAMC
420 CADJ STORE KPPghat (:,:,:,bi,bj)
421 CADJ & , KPPfrac (:,: ,bi,bj)
422 CADJ & = comlev1_bibj, key=itdkey, byte=isbyte
423 #endif /* ALLOW_AUTODIFF_TAMC */
424
425 #endif /* ALLOW_KPP */
426
427 #ifdef ALLOW_PP81
428 C-- Compute PP81 mixing coefficients
429 IF (usePP81) THEN
430 #ifdef ALLOW_DEBUG
431 IF ( debugLevel .GE. debLevB )
432 & CALL DEBUG_CALL('PP81_CALC',myThid)
433 #endif
434 CALL PP81_CALC(
435 I bi, bj, myTime, myThid )
436 ENDIF
437 #endif /* ALLOW_PP81 */
438
439 #ifdef ALLOW_MY82
440 C-- Compute MY82 mixing coefficients
441 IF (useMY82) THEN
442 #ifdef ALLOW_DEBUG
443 IF ( debugLevel .GE. debLevB )
444 & CALL DEBUG_CALL('MY82_CALC',myThid)
445 #endif
446 CALL MY82_CALC(
447 I bi, bj, myTime, myThid )
448 ENDIF
449 #endif /* ALLOW_MY82 */
450
451 #ifdef ALLOW_AUTODIFF_TAMC
452 CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
453 CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
454 CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
455 CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
456 #ifdef ALLOW_PASSIVE_TRACER
457 CADJ STORE tr1 (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
458 #endif
459 #ifdef ALLOW_PTRACERS
460 cph-- moved to forward_step to avoid key computation
461 cphCADJ STORE ptracer(:,:,:,bi,bj,itracer) = comlev1_bibj,
462 cphCADJ & key=itdkey, byte=isbyte
463 #endif
464 #endif /* ALLOW_AUTODIFF_TAMC */
465
466 C-- end bi,bj loops.
467 ENDDO
468 ENDDO
469
470 #ifdef ALLOW_DEBUG
471 IF ( debugLevel .GE. debLevB )
472 & CALL DEBUG_LEAVE('DO_OCEANIC_PHYS',myThid)
473 #endif
474
475 RETURN
476 END

  ViewVC Help
Powered by ViewVC 1.1.22