/[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.10 - (show annotations) (download)
Fri Sep 17 23:02:00 2004 UTC (19 years, 8 months ago) by heimbach
Branch: MAIN
CVS Tags: checkpoint55c_post, checkpoint55b_post, checkpoint55, checkpoint55a_post
Changes since 1.9: +15 -39 lines
o bringing adjoint up to date for sheduled c55

1 C $Header: /u/gcmpack/MITgcm/model/src/do_oceanic_phys.F,v 1.9 2004/09/16 09:35:11 mlosch 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_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 # 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 # ifndef GM_EXCLUDE_CLIPPING
300 CADJ STORE rhok (:,:) = comlev1_bibj_k , key=kkey, byte=isbyte
301 # 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+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 # ifndef GM_EXCLUDE_CLIPPING
377 cph storing here is needed only for one GMREDI_OPTIONS:
378 cph define GM_BOLUS_ADVEC
379 cph keep it although TAF says you dont need to.
380 cph but I've avoided the #ifdef for now, in case more things change
381 CADJ STORE sigmaX(:,:,:) = comlev1_bibj, key=itdkey, byte=isbyte
382 CADJ STORE sigmaY(:,:,:) = comlev1_bibj, key=itdkey, byte=isbyte
383 cnewCADJ STORE sigmaR(:,:,:) = comlev1_bibj, key=itdkey, byte=isbyte
384 # endif
385 #endif /* ALLOW_AUTODIFF_TAMC */
386
387 C-- Calculate iso-neutral slopes for the GM/Redi parameterisation
388 IF (useGMRedi) THEN
389 #ifdef ALLOW_DEBUG
390 IF ( debugLevel .GE. debLevB )
391 & CALL DEBUG_CALL('GMREDI_CALC_TENSOR',myThid)
392 #endif
393 CALL GMREDI_CALC_TENSOR(
394 I bi, bj, iMin, iMax, jMin, jMax,
395 I sigmaX, sigmaY, sigmaR,
396 I myThid )
397 #ifdef ALLOW_AUTODIFF_TAMC
398 ELSE
399 CALL GMREDI_CALC_TENSOR_DUMMY(
400 I bi, bj, iMin, iMax, jMin, jMax,
401 I sigmaX, sigmaY, sigmaR,
402 I myThid )
403 #endif /* ALLOW_AUTODIFF_TAMC */
404 ENDIF
405
406 #endif /* ALLOW_GMREDI */
407
408 #ifdef ALLOW_KPP
409 C-- Compute KPP mixing coefficients
410 IF (useKPP) THEN
411 #ifdef ALLOW_DEBUG
412 IF ( debugLevel .GE. debLevB )
413 & CALL DEBUG_CALL('KPP_CALC',myThid)
414 #endif
415 CALL KPP_CALC(
416 I bi, bj, myTime, myThid )
417 #ifdef ALLOW_AUTODIFF_TAMC
418 ELSE
419 CALL KPP_CALC_DUMMY(
420 I bi, bj, myTime, myThid )
421 #endif /* ALLOW_AUTODIFF_TAMC */
422 ENDIF
423
424 #endif /* ALLOW_KPP */
425
426 #ifdef ALLOW_PP81
427 C-- Compute PP81 mixing coefficients
428 IF (usePP81) THEN
429 #ifdef ALLOW_DEBUG
430 IF ( debugLevel .GE. debLevB )
431 & CALL DEBUG_CALL('PP81_CALC',myThid)
432 #endif
433 CALL PP81_CALC(
434 I bi, bj, myTime, myThid )
435 ENDIF
436 #endif /* ALLOW_PP81 */
437
438 #ifdef ALLOW_MY82
439 C-- Compute MY82 mixing coefficients
440 IF (useMY82) THEN
441 #ifdef ALLOW_DEBUG
442 IF ( debugLevel .GE. debLevB )
443 & CALL DEBUG_CALL('MY82_CALC',myThid)
444 #endif
445 CALL MY82_CALC(
446 I bi, bj, myTime, myThid )
447 ENDIF
448 #endif /* ALLOW_MY82 */
449
450 #ifdef ALLOW_GGL90
451 C-- Compute GGL90 mixing coefficients
452 IF (useGGL90) THEN
453 #ifdef ALLOW_DEBUG
454 IF ( debugLevel .GE. debLevB )
455 & CALL DEBUG_CALL('GGL90_CALC',myThid)
456 #endif
457 CALL GGL90_CALC(
458 I bi, bj, myTime, myThid )
459 ENDIF
460 #endif /* ALLOW_GGL90 */
461
462 C-- end bi,bj loops.
463 ENDDO
464 ENDDO
465
466 #ifdef ALLOW_DEBUG
467 IF ( debugLevel .GE. debLevB )
468 & CALL DEBUG_LEAVE('DO_OCEANIC_PHYS',myThid)
469 #endif
470
471 RETURN
472 END

  ViewVC Help
Powered by ViewVC 1.1.22