/[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.8 - (show annotations) (download)
Wed Sep 8 01:50:24 2004 UTC (19 years, 8 months ago) by jmc
Branch: MAIN
Changes since 1.7: +9 -1 lines
add diagnostics od d.Rho/dr (=SigmaR).

1 C $Header: /u/gcmpack/MITgcm/model/src/do_oceanic_phys.F,v 1.7 2004/09/07 17:29:14 edhill 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: do not 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_DIAGNOSTICS
316 IF ( usediagnostics .AND.
317 & (useGMRedi .OR. ivdc_kappa.NE.0.) ) THEN
318 CALL fill_diagnostics (myThid, 'DRHODR ', 0, Nr,
319 & 3, bi, bj, sigmaR)
320 ENDIF
321 #endif
322
323 #ifdef ALLOW_AUTODIFF_TAMC
324 cph avoids recomputation of integrate_for_w
325 CADJ STORE wvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
326 #endif /* ALLOW_AUTODIFF_TAMC */
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+deltaT, 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 ( buoyancyRelation(1:7) .EQ. 'OCEANIC' ) 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 # ifdef ALLOW_SEAICE
368 CADJ STORE surfaceForcingTice(:,:,bi,bj)
369 CADJ & = comlev1_bibj, key=itdkey, byte=isbyte
370 # endif
371 #endif /* ALLOW_AUTODIFF_TAMC */
372
373 #ifdef ALLOW_GMREDI
374
375 #ifdef ALLOW_AUTODIFF_TAMC
376 cph storing here is needed only for one GMREDI_OPTIONS:
377 cph define GM_BOLUS_ADVEC
378 cph but I've avoided the #ifdef for now, in case more things change
379 CADJ STORE sigmaX(:,:,:) = comlev1_bibj, key=itdkey, byte=isbyte
380 CADJ STORE sigmaY(:,:,:) = comlev1_bibj, key=itdkey, byte=isbyte
381 CADJ STORE sigmaR(:,:,:) = comlev1_bibj, key=itdkey, byte=isbyte
382 #endif /* ALLOW_AUTODIFF_TAMC */
383
384 C-- Calculate iso-neutral slopes for the GM/Redi parameterisation
385 IF (useGMRedi) THEN
386 #ifdef ALLOW_DEBUG
387 IF ( debugLevel .GE. debLevB )
388 & CALL DEBUG_CALL('GMREDI_CALC_TENSOR',myThid)
389 #endif
390 CALL GMREDI_CALC_TENSOR(
391 I bi, bj, iMin, iMax, jMin, jMax,
392 I sigmaX, sigmaY, sigmaR,
393 I myThid )
394 #ifdef ALLOW_AUTODIFF_TAMC
395 ELSE
396 CALL GMREDI_CALC_TENSOR_DUMMY(
397 I bi, bj, iMin, iMax, jMin, jMax,
398 I sigmaX, sigmaY, sigmaR,
399 I myThid )
400 #endif /* ALLOW_AUTODIFF_TAMC */
401 ENDIF
402
403 #ifdef ALLOW_AUTODIFF_TAMC
404 CADJ STORE Kwx(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
405 CADJ STORE Kwy(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
406 CADJ STORE Kwz(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
407 #endif /* ALLOW_AUTODIFF_TAMC */
408
409 #endif /* ALLOW_GMREDI */
410
411 #ifdef ALLOW_KPP
412 C-- Compute KPP mixing coefficients
413 IF (useKPP) THEN
414 #ifdef ALLOW_DEBUG
415 IF ( debugLevel .GE. debLevB )
416 & CALL DEBUG_CALL('KPP_CALC',myThid)
417 #endif
418 CALL KPP_CALC(
419 I bi, bj, myTime, myThid )
420 #ifdef ALLOW_AUTODIFF_TAMC
421 ELSE
422 CALL KPP_CALC_DUMMY(
423 I bi, bj, myTime, myThid )
424 #endif /* ALLOW_AUTODIFF_TAMC */
425 ENDIF
426
427 #ifdef ALLOW_AUTODIFF_TAMC
428 CADJ STORE KPPghat (:,:,:,bi,bj)
429 CADJ & , KPPfrac (:,: ,bi,bj)
430 CADJ & = comlev1_bibj, key=itdkey, byte=isbyte
431 #endif /* ALLOW_AUTODIFF_TAMC */
432
433 #endif /* ALLOW_KPP */
434
435 #ifdef ALLOW_PP81
436 C-- Compute PP81 mixing coefficients
437 IF (usePP81) THEN
438 #ifdef ALLOW_DEBUG
439 IF ( debugLevel .GE. debLevB )
440 & CALL DEBUG_CALL('PP81_CALC',myThid)
441 #endif
442 CALL PP81_CALC(
443 I bi, bj, myTime, myThid )
444 ENDIF
445 #endif /* ALLOW_PP81 */
446
447 #ifdef ALLOW_MY82
448 C-- Compute MY82 mixing coefficients
449 IF (useMY82) THEN
450 #ifdef ALLOW_DEBUG
451 IF ( debugLevel .GE. debLevB )
452 & CALL DEBUG_CALL('MY82_CALC',myThid)
453 #endif
454 CALL MY82_CALC(
455 I bi, bj, myTime, myThid )
456 ENDIF
457 #endif /* ALLOW_MY82 */
458
459 #ifdef ALLOW_AUTODIFF_TAMC
460 CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
461 CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
462 CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
463 CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
464 #ifdef ALLOW_PASSIVE_TRACER
465 CADJ STORE tr1 (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
466 #endif
467 #ifdef ALLOW_PTRACERS
468 cph-- moved to forward_step to avoid key computation
469 cphCADJ STORE ptracer(:,:,:,bi,bj,itracer) = comlev1_bibj,
470 cphCADJ & key=itdkey, byte=isbyte
471 #endif
472 #endif /* ALLOW_AUTODIFF_TAMC */
473
474 C-- end bi,bj loops.
475 ENDDO
476 ENDDO
477
478 #ifdef ALLOW_DEBUG
479 IF ( debugLevel .GE. debLevB )
480 & CALL DEBUG_LEAVE('DO_OCEANIC_PHYS',myThid)
481 #endif
482
483 RETURN
484 END

  ViewVC Help
Powered by ViewVC 1.1.22