/[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.12 - (hide annotations) (download)
Wed Oct 13 07:05:50 2004 UTC (19 years, 7 months ago) by heimbach
Branch: MAIN
CVS Tags: checkpoint55f_post, checkpoint55e_post
Changes since 1.11: +12 -14 lines
o some delicate re-shuffle of store directives to avoid one
  extra call of do_oceanic_physics
o NB: this may break global_ocean adjoint temporarily,
  but it is clear how to fix it. Will do later, need this now.

1 heimbach 1.12 C $Header: /u/gcmpack/MITgcm/model/src/do_oceanic_phys.F,v 1.11 2004/09/27 14:56:42 heimbach 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     IF ( useThSIce .AND. buoyancyRelation(1:7) .EQ. 'OCEANIC' ) THEN
113     #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    
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 heimbach 1.12 CADJ STORE totphihyd(:,:,:,bi,bj)
232 jmc 1.1 CADJ & = comlev1_bibj, key=itdkey, byte=isbyte
233 heimbach 1.10 # ifdef ALLOW_KPP
234 jmc 1.1 CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
235     CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
236 heimbach 1.10 # endif
237     # ifdef EXACT_CONSERV
238     CADJ STORE pmepr(:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
239     # endif
240 jmc 1.1 #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 heimbach 1.12 ctest# ifndef GM_EXCLUDE_CLIPPING
300 jmc 1.1 CADJ STORE rhok (:,:) = comlev1_bibj_k , key=kkey, byte=isbyte
301 heimbach 1.12 ctest# endif
302 jmc 1.1 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 jmc 1.8 #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 jmc 1.1 #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 heimbach 1.11 CALL OBCS_CALC( bi, bj, myTime+deltaTclock, myIter+1,
336 jmc 1.1 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 jmc 1.4 CADJ STORE surfaceForcingU(:,:,bi,bj)
360 jmc 1.1 CADJ & = comlev1_bibj, key=itdkey, byte=isbyte
361 jmc 1.4 CADJ STORE surfaceForcingV(:,:,bi,bj)
362 jmc 1.1 CADJ & = comlev1_bibj, key=itdkey, byte=isbyte
363 jmc 1.4 CADJ STORE surfaceForcingS(:,:,bi,bj)
364 jmc 1.1 CADJ & = comlev1_bibj, key=itdkey, byte=isbyte
365 jmc 1.4 CADJ STORE surfaceForcingT(:,:,bi,bj)
366 jmc 1.1 CADJ & = comlev1_bibj, key=itdkey, byte=isbyte
367 jmc 1.4 CADJ STORE surfaceForcingTice(:,:,bi,bj)
368 jmc 1.1 CADJ & = comlev1_bibj, key=itdkey, byte=isbyte
369     #endif /* ALLOW_AUTODIFF_TAMC */
370    
371     #ifdef ALLOW_GMREDI
372    
373     #ifdef ALLOW_AUTODIFF_TAMC
374 heimbach 1.10 # ifndef GM_EXCLUDE_CLIPPING
375 jmc 1.1 cph storing here is needed only for one GMREDI_OPTIONS:
376     cph define GM_BOLUS_ADVEC
377 heimbach 1.10 cph keep it although TAF says you dont need to.
378 jmc 1.1 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 heimbach 1.12 CADJ STORE sigmaR(:,:,:) = comlev1_bibj, key=itdkey, byte=isbyte
382 heimbach 1.10 # endif
383 jmc 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
384    
385     C-- Calculate iso-neutral slopes for the GM/Redi parameterisation
386     IF (useGMRedi) THEN
387     #ifdef ALLOW_DEBUG
388     IF ( debugLevel .GE. debLevB )
389     & CALL DEBUG_CALL('GMREDI_CALC_TENSOR',myThid)
390     #endif
391     CALL GMREDI_CALC_TENSOR(
392     I bi, bj, iMin, iMax, jMin, jMax,
393     I sigmaX, sigmaY, sigmaR,
394     I myThid )
395     #ifdef ALLOW_AUTODIFF_TAMC
396     ELSE
397     CALL GMREDI_CALC_TENSOR_DUMMY(
398     I bi, bj, iMin, iMax, jMin, jMax,
399     I sigmaX, sigmaY, sigmaR,
400     I myThid )
401     #endif /* ALLOW_AUTODIFF_TAMC */
402     ENDIF
403    
404     #endif /* ALLOW_GMREDI */
405    
406     #ifdef ALLOW_KPP
407     C-- Compute KPP mixing coefficients
408     IF (useKPP) THEN
409     #ifdef ALLOW_DEBUG
410     IF ( debugLevel .GE. debLevB )
411     & CALL DEBUG_CALL('KPP_CALC',myThid)
412     #endif
413     CALL KPP_CALC(
414     I bi, bj, myTime, myThid )
415     #ifdef ALLOW_AUTODIFF_TAMC
416     ELSE
417     CALL KPP_CALC_DUMMY(
418     I bi, bj, myTime, myThid )
419     #endif /* ALLOW_AUTODIFF_TAMC */
420     ENDIF
421    
422     #endif /* ALLOW_KPP */
423    
424 mlosch 1.6 #ifdef ALLOW_PP81
425     C-- Compute PP81 mixing coefficients
426     IF (usePP81) THEN
427     #ifdef ALLOW_DEBUG
428     IF ( debugLevel .GE. debLevB )
429     & CALL DEBUG_CALL('PP81_CALC',myThid)
430     #endif
431     CALL PP81_CALC(
432     I bi, bj, myTime, myThid )
433     ENDIF
434     #endif /* ALLOW_PP81 */
435    
436     #ifdef ALLOW_MY82
437     C-- Compute MY82 mixing coefficients
438     IF (useMY82) THEN
439     #ifdef ALLOW_DEBUG
440     IF ( debugLevel .GE. debLevB )
441     & CALL DEBUG_CALL('MY82_CALC',myThid)
442     #endif
443     CALL MY82_CALC(
444     I bi, bj, myTime, myThid )
445     ENDIF
446     #endif /* ALLOW_MY82 */
447    
448 mlosch 1.9 #ifdef ALLOW_GGL90
449     C-- Compute GGL90 mixing coefficients
450     IF (useGGL90) THEN
451     #ifdef ALLOW_DEBUG
452     IF ( debugLevel .GE. debLevB )
453     & CALL DEBUG_CALL('GGL90_CALC',myThid)
454     #endif
455     CALL GGL90_CALC(
456     I bi, bj, myTime, myThid )
457     ENDIF
458     #endif /* ALLOW_GGL90 */
459    
460 jmc 1.1 C-- end bi,bj loops.
461     ENDDO
462     ENDDO
463    
464     #ifdef ALLOW_DEBUG
465     IF ( debugLevel .GE. debLevB )
466     & CALL DEBUG_LEAVE('DO_OCEANIC_PHYS',myThid)
467     #endif
468    
469     RETURN
470     END

  ViewVC Help
Powered by ViewVC 1.1.22