/[MITgcm]/MITgcm/model/src/do_oceanic_phys.F
ViewVC logotype

Annotation of /MITgcm/model/src/do_oceanic_phys.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.15 - (hide annotations) (download)
Fri Dec 10 20:15:10 2004 UTC (19 years, 6 months ago) by heimbach
Branch: MAIN
CVS Tags: checkpoint57a_post
Changes since 1.14: +3 -8 lines
This set of changes restores TAMC(!) compatibility.

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

  ViewVC Help
Powered by ViewVC 1.1.22