/[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.7 - (hide annotations) (download)
Tue Sep 7 17:29:14 2004 UTC (19 years, 8 months ago) by edhill
Branch: MAIN
Changes since 1.6: +2 -2 lines
 o remove single quotes (eg. "don't"-->"do not") so that the on-line code
   browser does not get confused

1 edhill 1.7 C $Header: /u/gcmpack/MITgcm/model/src/do_oceanic_phys.F,v 1.6 2004/09/02 09:13:49 mlosch 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_PASSIVE_TRACER
40     c #include "TR1.h"
41     c #endif
42     c #ifdef ALLOW_PTRACERS
43 jmc 1.3 c #include "PTRACERS_SIZE.h"
44 jmc 1.1 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 jmc 1.5 IF ( debugLevel .GE. debLevB )
103 jmc 1.1 & CALL DEBUG_ENTER('DO_OCEANIC_PHYS',myThid)
104     #endif
105    
106 jmc 1.5 #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 edhill 1.7 C jmc: do not know precisely where to put this call (bf or af thSIce ?)
133 jmc 1.5 IF ( useCoupler ) THEN
134     CALL OCN_APPLY_IMPORT( .TRUE., myTime, myIter, myThid )
135     ENDIF
136     # endif
137     #endif /* COMPONENT_MODULE */
138    
139 jmc 1.1 #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_AUTODIFF_TAMC
316     cph avoids recomputation of integrate_for_w
317     CADJ STORE wvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
318     #endif /* ALLOW_AUTODIFF_TAMC */
319    
320     #ifdef ALLOW_OBCS
321     C-- Calculate future values on open boundaries
322     IF (useOBCS) THEN
323     #ifdef ALLOW_DEBUG
324     IF ( debugLevel .GE. debLevB )
325     & CALL DEBUG_CALL('OBCS_CALC',myThid)
326     #endif
327     CALL OBCS_CALC( bi, bj, myTime+deltaT, myIter+1,
328     I uVel, vVel, wVel, theta, salt,
329     I myThid )
330     ENDIF
331     #endif /* ALLOW_OBCS */
332    
333     #ifndef ALLOW_AUTODIFF_TAMC
334     IF ( buoyancyRelation(1:7) .EQ. 'OCEANIC' ) THEN
335     #endif
336     C-- Determines forcing terms based on external fields
337     C relaxation terms, etc.
338     #ifdef ALLOW_DEBUG
339     IF ( debugLevel .GE. debLevB )
340     & CALL DEBUG_CALL('EXTERNAL_FORCING_SURF',myThid)
341     #endif
342     CALL EXTERNAL_FORCING_SURF(
343     I bi, bj, iMin, iMax, jMin, jMax,
344     I myTime, myIter, myThid )
345     #ifndef ALLOW_AUTODIFF_TAMC
346     ENDIF
347     #endif
348    
349     #ifdef ALLOW_AUTODIFF_TAMC
350     cph needed for KPP
351 jmc 1.4 CADJ STORE surfaceForcingU(:,:,bi,bj)
352 jmc 1.1 CADJ & = comlev1_bibj, key=itdkey, byte=isbyte
353 jmc 1.4 CADJ STORE surfaceForcingV(:,:,bi,bj)
354 jmc 1.1 CADJ & = comlev1_bibj, key=itdkey, byte=isbyte
355 jmc 1.4 CADJ STORE surfaceForcingS(:,:,bi,bj)
356 jmc 1.1 CADJ & = comlev1_bibj, key=itdkey, byte=isbyte
357 jmc 1.4 CADJ STORE surfaceForcingT(:,:,bi,bj)
358 jmc 1.1 CADJ & = comlev1_bibj, key=itdkey, byte=isbyte
359     # ifdef ALLOW_SEAICE
360 jmc 1.4 CADJ STORE surfaceForcingTice(:,:,bi,bj)
361 jmc 1.1 CADJ & = comlev1_bibj, key=itdkey, byte=isbyte
362     # endif
363     #endif /* ALLOW_AUTODIFF_TAMC */
364    
365     #ifdef ALLOW_GMREDI
366    
367     #ifdef ALLOW_AUTODIFF_TAMC
368     cph storing here is needed only for one GMREDI_OPTIONS:
369     cph define GM_BOLUS_ADVEC
370     cph but I've avoided the #ifdef for now, in case more things change
371     CADJ STORE sigmaX(:,:,:) = comlev1_bibj, key=itdkey, byte=isbyte
372     CADJ STORE sigmaY(:,:,:) = comlev1_bibj, key=itdkey, byte=isbyte
373     CADJ STORE sigmaR(:,:,:) = comlev1_bibj, key=itdkey, byte=isbyte
374     #endif /* ALLOW_AUTODIFF_TAMC */
375    
376     C-- Calculate iso-neutral slopes for the GM/Redi parameterisation
377     IF (useGMRedi) THEN
378     #ifdef ALLOW_DEBUG
379     IF ( debugLevel .GE. debLevB )
380     & CALL DEBUG_CALL('GMREDI_CALC_TENSOR',myThid)
381     #endif
382     CALL GMREDI_CALC_TENSOR(
383     I bi, bj, iMin, iMax, jMin, jMax,
384     I sigmaX, sigmaY, sigmaR,
385     I myThid )
386     #ifdef ALLOW_AUTODIFF_TAMC
387     ELSE
388     CALL GMREDI_CALC_TENSOR_DUMMY(
389     I bi, bj, iMin, iMax, jMin, jMax,
390     I sigmaX, sigmaY, sigmaR,
391     I myThid )
392     #endif /* ALLOW_AUTODIFF_TAMC */
393     ENDIF
394    
395     #ifdef ALLOW_AUTODIFF_TAMC
396     CADJ STORE Kwx(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
397     CADJ STORE Kwy(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
398     CADJ STORE Kwz(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
399     #endif /* ALLOW_AUTODIFF_TAMC */
400    
401     #endif /* ALLOW_GMREDI */
402    
403     #ifdef ALLOW_KPP
404     C-- Compute KPP mixing coefficients
405     IF (useKPP) THEN
406     #ifdef ALLOW_DEBUG
407     IF ( debugLevel .GE. debLevB )
408     & CALL DEBUG_CALL('KPP_CALC',myThid)
409     #endif
410     CALL KPP_CALC(
411     I bi, bj, myTime, myThid )
412     #ifdef ALLOW_AUTODIFF_TAMC
413     ELSE
414     CALL KPP_CALC_DUMMY(
415     I bi, bj, myTime, myThid )
416     #endif /* ALLOW_AUTODIFF_TAMC */
417     ENDIF
418    
419     #ifdef ALLOW_AUTODIFF_TAMC
420     CADJ STORE KPPghat (:,:,:,bi,bj)
421     CADJ & , KPPfrac (:,: ,bi,bj)
422     CADJ & = comlev1_bibj, key=itdkey, byte=isbyte
423     #endif /* ALLOW_AUTODIFF_TAMC */
424    
425     #endif /* ALLOW_KPP */
426    
427 mlosch 1.6 #ifdef ALLOW_PP81
428     C-- Compute PP81 mixing coefficients
429     IF (usePP81) THEN
430     #ifdef ALLOW_DEBUG
431     IF ( debugLevel .GE. debLevB )
432     & CALL DEBUG_CALL('PP81_CALC',myThid)
433     #endif
434     CALL PP81_CALC(
435     I bi, bj, myTime, myThid )
436     ENDIF
437     #endif /* ALLOW_PP81 */
438    
439     #ifdef ALLOW_MY82
440     C-- Compute MY82 mixing coefficients
441     IF (useMY82) THEN
442     #ifdef ALLOW_DEBUG
443     IF ( debugLevel .GE. debLevB )
444     & CALL DEBUG_CALL('MY82_CALC',myThid)
445     #endif
446     CALL MY82_CALC(
447     I bi, bj, myTime, myThid )
448     ENDIF
449     #endif /* ALLOW_MY82 */
450    
451 jmc 1.1 #ifdef ALLOW_AUTODIFF_TAMC
452     CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
453     CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
454     CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
455     CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
456     #ifdef ALLOW_PASSIVE_TRACER
457     CADJ STORE tr1 (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
458     #endif
459     #ifdef ALLOW_PTRACERS
460     cph-- moved to forward_step to avoid key computation
461     cphCADJ STORE ptracer(:,:,:,bi,bj,itracer) = comlev1_bibj,
462     cphCADJ & key=itdkey, byte=isbyte
463     #endif
464     #endif /* ALLOW_AUTODIFF_TAMC */
465    
466     C-- end bi,bj loops.
467     ENDDO
468     ENDDO
469    
470     #ifdef ALLOW_DEBUG
471     IF ( debugLevel .GE. debLevB )
472     & CALL DEBUG_LEAVE('DO_OCEANIC_PHYS',myThid)
473     #endif
474    
475     RETURN
476     END

  ViewVC Help
Powered by ViewVC 1.1.22