/[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.19 - (show annotations) (download)
Thu Sep 15 14:55:15 2005 UTC (18 years, 8 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint57t_post, checkpoint57s_post, checkpoint57v_post, checkpoint57w_post, checkpint57u_post
Changes since 1.18: +5 -1 lines
add diagnostics for Convective Adjusment.

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

  ViewVC Help
Powered by ViewVC 1.1.22