/[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.1 - (hide annotations) (download)
Tue Jul 6 00:58:40 2004 UTC (19 years, 11 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint54a_post
take out 1rst part of thermodynamics (KPP,GM,slope ..) into do_oceanic_phys.F

1 jmc 1.1 C $Header: $
2     C $Name: $
3    
4     #include "PACKAGES_CONFIG.h"
5     #include "CPP_OPTIONS.h"
6     c #ifdef ALLOW_PTRACERS
7     c # include "PTRACERS_OPTIONS.h"
8     c #endif
9    
10     #ifdef ALLOW_AUTODIFF_TAMC
11     # ifdef ALLOW_GMREDI
12     # include "GMREDI_OPTIONS.h"
13     # endif
14     # ifdef ALLOW_KPP
15     # include "KPP_OPTIONS.h"
16     # endif
17     #endif /* ALLOW_AUTODIFF_TAMC */
18    
19     CBOP
20     C !ROUTINE: DO_OCEANIC_PHYS
21     C !INTERFACE:
22     SUBROUTINE DO_OCEANIC_PHYS(myTime, myIter, myThid)
23     C !DESCRIPTION: \bv
24     C *==========================================================*
25     C | SUBROUTINE DO_OCEANIC_PHYS
26     C | o Controlling routine for oceanic physics and
27     C | parameterization
28     C *==========================================================*
29     C | o originally, part of S/R thermodynamics
30     C *==========================================================*
31     C \ev
32    
33     C !USES:
34     IMPLICIT NONE
35     C == Global variables ===
36     #include "SIZE.h"
37     #include "EEPARAMS.h"
38     #include "PARAMS.h"
39     #include "DYNVARS.h"
40     #include "GRID.h"
41     c #include "GAD.h"
42     c #ifdef ALLOW_PASSIVE_TRACER
43     c #include "TR1.h"
44     c #endif
45     c #ifdef ALLOW_PTRACERS
46     c #include "PTRACERS.h"
47     c #endif
48     c #ifdef ALLOW_TIMEAVE
49     c #include "TIMEAVE_STATV.h"
50     c #endif
51    
52     #ifdef ALLOW_AUTODIFF_TAMC
53     # include "tamc.h"
54     # include "tamc_keys.h"
55     # include "FFIELDS.h"
56     # include "EOS.h"
57     # ifdef ALLOW_KPP
58     # include "KPP.h"
59     # endif
60     # ifdef ALLOW_GMREDI
61     # include "GMREDI.h"
62     # endif
63     # ifdef ALLOW_EBM
64     # include "EBM.h"
65     # endif
66     #endif /* ALLOW_AUTODIFF_TAMC */
67    
68     C !INPUT/OUTPUT PARAMETERS:
69     C == Routine arguments ==
70     C myTime - Current time in simulation
71     C myIter - Current iteration number in simulation
72     C myThid - Thread number for this instance of the routine.
73     _RL myTime
74     INTEGER myIter
75     INTEGER myThid
76    
77     C !LOCAL VARIABLES:
78     C == Local variables
79     C rhoK, rhoKM1 - Density at current level, and level above
80     C useVariableK = T when vertical diffusion is not constant
81     C iMin, iMax - Ranges and sub-block indices on which calculations
82     C jMin, jMax are applied.
83     C bi, bj
84     C k, kup, - Index for layer above and below. kup and kDown
85     C kDown, km1 are switched with layer to be the appropriate
86     C index into fVerTerm.
87     _RL rhokm1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
88     _RL rhok (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
89     _RL sigmaX (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
90     _RL sigmaY (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
91     _RL sigmaR (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
92     _RL kp1Msk
93     LOGICAL useVariableK
94     INTEGER iMin, iMax
95     INTEGER jMin, jMax
96     INTEGER bi, bj
97     INTEGER i, j
98     INTEGER k, km1, kup, kDown
99     INTEGER iTracer, ip
100    
101     CEOP
102    
103     #ifdef ALLOW_DEBUG
104     IF ( debugLevel .GE. debLevB )
105     & CALL DEBUG_ENTER('DO_OCEANIC_PHYS',myThid)
106     #endif
107    
108     #ifdef ALLOW_AUTODIFF_TAMC
109     C-- dummy statement to end declaration part
110     ikey = 1
111     itdkey = 1
112     #endif /* ALLOW_AUTODIFF_TAMC */
113    
114     #ifdef ALLOW_AUTODIFF_TAMC
115     C-- HPF directive to help TAMC
116     CHPF$ INDEPENDENT
117     #endif /* ALLOW_AUTODIFF_TAMC */
118    
119     DO bj=myByLo(myThid),myByHi(myThid)
120    
121     #ifdef ALLOW_AUTODIFF_TAMC
122     C-- HPF directive to help TAMC
123     CHPF$ INDEPENDENT, NEW (rTrans,fVerT,fVerS
124     CHPF$& ,utrans,vtrans,xA,yA
125     CHPF$& )
126     #endif /* ALLOW_AUTODIFF_TAMC */
127    
128     DO bi=myBxLo(myThid),myBxHi(myThid)
129    
130     #ifdef ALLOW_AUTODIFF_TAMC
131     act1 = bi - myBxLo(myThid)
132     max1 = myBxHi(myThid) - myBxLo(myThid) + 1
133     act2 = bj - myByLo(myThid)
134     max2 = myByHi(myThid) - myByLo(myThid) + 1
135     act3 = myThid - 1
136     max3 = nTx*nTy
137     act4 = ikey_dynamics - 1
138     itdkey = (act1 + 1) + act2*max1
139     & + act3*max1*max2
140     & + act4*max1*max2*max3
141     #endif /* ALLOW_AUTODIFF_TAMC */
142    
143     C-- Set up work arrays with valid (i.e. not NaN) values
144     C These inital values do not alter the numerical results. They
145     C just ensure that all memory references are to valid floating
146     C point numbers. This prevents spurious hardware signals due to
147     C uninitialised but inert locations.
148    
149     DO j=1-OLy,sNy+OLy
150     DO i=1-OLx,sNx+OLx
151     rhok (i,j) = 0. _d 0
152     rhoKM1 (i,j) = 0. _d 0
153     ENDDO
154     ENDDO
155    
156     DO k=1,Nr
157     DO j=1-OLy,sNy+OLy
158     DO i=1-OLx,sNx+OLx
159     C This is currently also used by IVDC and Diagnostics
160     sigmaX(i,j,k) = 0. _d 0
161     sigmaY(i,j,k) = 0. _d 0
162     sigmaR(i,j,k) = 0. _d 0
163     #ifdef ALLOW_AUTODIFF_TAMC
164     cph all the following init. are necessary for TAF
165     cph although some of these are re-initialised later.
166     IVDConvCount(i,j,k,bi,bj) = 0.
167     # ifdef ALLOW_GMREDI
168     Kwx(i,j,k,bi,bj) = 0. _d 0
169     Kwy(i,j,k,bi,bj) = 0. _d 0
170     Kwz(i,j,k,bi,bj) = 0. _d 0
171     # ifdef GM_NON_UNITY_DIAGONAL
172     Kux(i,j,k,bi,bj) = 0. _d 0
173     Kvy(i,j,k,bi,bj) = 0. _d 0
174     # endif
175     # ifdef GM_EXTRA_DIAGONAL
176     Kuz(i,j,k,bi,bj) = 0. _d 0
177     Kvz(i,j,k,bi,bj) = 0. _d 0
178     # endif
179     # ifdef GM_BOLUS_ADVEC
180     GM_PsiX(i,j,k,bi,bj) = 0. _d 0
181     GM_PsiY(i,j,k,bi,bj) = 0. _d 0
182     # endif
183     # ifdef GM_VISBECK_VARIABLE_K
184     VisbeckK(i,j,bi,bj) = 0. _d 0
185     # endif
186     # endif /* ALLOW_GMREDI */
187     #endif /* ALLOW_AUTODIFF_TAMC */
188     ENDDO
189     ENDDO
190     ENDDO
191    
192     iMin = 1-OLx
193     iMax = sNx+OLx
194     jMin = 1-OLy
195     jMax = sNy+OLy
196    
197     #ifdef ALLOW_AUTODIFF_TAMC
198     CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
199     CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
200     CADJ STORE totphihyd
201     CADJ & = comlev1_bibj, key=itdkey, byte=isbyte
202     #ifdef ALLOW_KPP
203     CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
204     CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
205     #endif
206     #endif /* ALLOW_AUTODIFF_TAMC */
207    
208     #ifdef ALLOW_DEBUG
209     IF ( debugLevel .GE. debLevB )
210     & CALL DEBUG_MSG('ENTERING UPWARD K LOOP',myThid)
211     #endif
212    
213     C-- Start of diagnostic loop
214     DO k=Nr,1,-1
215    
216     #ifdef ALLOW_AUTODIFF_TAMC
217     C? Patrick, is this formula correct now that we change the loop range?
218     C? Do we still need this?
219     cph kkey formula corrected.
220     cph Needed for rhok, rhokm1, in the case useGMREDI.
221     kkey = (itdkey-1)*Nr + k
222     #endif /* ALLOW_AUTODIFF_TAMC */
223    
224     C-- Calculate gradients of potential density for isoneutral
225     C slope terms (e.g. GM/Redi tensor or IVDC diffusivity)
226     c IF ( k.GT.1 .AND. (useGMRedi.OR.ivdc_kappa.NE.0.) ) THEN
227     IF ( useGMRedi .OR. (k.GT.1 .AND. ivdc_kappa.NE.0.) ) THEN
228     #ifdef ALLOW_DEBUG
229     IF ( debugLevel .GE. debLevB )
230     & CALL DEBUG_CALL('FIND_RHO',myThid)
231     #endif
232     #ifdef ALLOW_AUTODIFF_TAMC
233     CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
234     CADJ STORE salt (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
235     #endif /* ALLOW_AUTODIFF_TAMC */
236     CALL FIND_RHO(
237     I bi, bj, iMin, iMax, jMin, jMax, k, k,
238     I theta, salt,
239     O rhoK,
240     I myThid )
241    
242     IF (k.GT.1) THEN
243     #ifdef ALLOW_AUTODIFF_TAMC
244     CADJ STORE theta(:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
245     CADJ STORE salt (:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
246     #endif /* ALLOW_AUTODIFF_TAMC */
247     CALL FIND_RHO(
248     I bi, bj, iMin, iMax, jMin, jMax, k-1, k,
249     I theta, salt,
250     O rhoKm1,
251     I myThid )
252     ENDIF
253     #ifdef ALLOW_DEBUG
254     IF ( debugLevel .GE. debLevB )
255     & CALL DEBUG_CALL('GRAD_SIGMA',myThid)
256     #endif
257     CALL GRAD_SIGMA(
258     I bi, bj, iMin, iMax, jMin, jMax, k,
259     I rhoK, rhoKm1, rhoK,
260     O sigmaX, sigmaY, sigmaR,
261     I myThid )
262     ENDIF
263    
264     #ifdef ALLOW_AUTODIFF_TAMC
265     CADJ STORE rhok (:,:) = comlev1_bibj_k , key=kkey, byte=isbyte
266     CADJ STORE rhokm1 (:,:) = comlev1_bibj_k , key=kkey, byte=isbyte
267     #endif /* ALLOW_AUTODIFF_TAMC */
268     C-- Implicit Vertical Diffusion for Convection
269     c ==> should use sigmaR !!!
270     IF (k.GT.1 .AND. ivdc_kappa.NE.0.) THEN
271     #ifdef ALLOW_DEBUG
272     IF ( debugLevel .GE. debLevB )
273     & CALL DEBUG_CALL('CALC_IVDC',myThid)
274     #endif
275     CALL CALC_IVDC(
276     I bi, bj, iMin, iMax, jMin, jMax, k,
277     I rhoKm1, rhoK,
278     I myTime, myIter, myThid)
279     ENDIF
280    
281     C-- end of diagnostic k loop (Nr:1)
282     ENDDO
283    
284     #ifdef ALLOW_AUTODIFF_TAMC
285     cph avoids recomputation of integrate_for_w
286     CADJ STORE wvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
287     #endif /* ALLOW_AUTODIFF_TAMC */
288    
289     #ifdef ALLOW_OBCS
290     C-- Calculate future values on open boundaries
291     IF (useOBCS) THEN
292     #ifdef ALLOW_DEBUG
293     IF ( debugLevel .GE. debLevB )
294     & CALL DEBUG_CALL('OBCS_CALC',myThid)
295     #endif
296     CALL OBCS_CALC( bi, bj, myTime+deltaT, myIter+1,
297     I uVel, vVel, wVel, theta, salt,
298     I myThid )
299     ENDIF
300     #endif /* ALLOW_OBCS */
301    
302     #ifndef ALLOW_AUTODIFF_TAMC
303     IF ( buoyancyRelation(1:7) .EQ. 'OCEANIC' ) THEN
304     #endif
305     C-- Determines forcing terms based on external fields
306     C relaxation terms, etc.
307     #ifdef ALLOW_DEBUG
308     IF ( debugLevel .GE. debLevB )
309     & CALL DEBUG_CALL('EXTERNAL_FORCING_SURF',myThid)
310     #endif
311     CALL EXTERNAL_FORCING_SURF(
312     I bi, bj, iMin, iMax, jMin, jMax,
313     I myTime, myIter, myThid )
314     #ifndef ALLOW_AUTODIFF_TAMC
315     ENDIF
316     #endif
317    
318     #ifdef ALLOW_AUTODIFF_TAMC
319     cph needed for KPP
320     CADJ STORE surfacetendencyU(:,:,bi,bj)
321     CADJ & = comlev1_bibj, key=itdkey, byte=isbyte
322     CADJ STORE surfacetendencyV(:,:,bi,bj)
323     CADJ & = comlev1_bibj, key=itdkey, byte=isbyte
324     CADJ STORE surfacetendencyS(:,:,bi,bj)
325     CADJ & = comlev1_bibj, key=itdkey, byte=isbyte
326     CADJ STORE surfacetendencyT(:,:,bi,bj)
327     CADJ & = comlev1_bibj, key=itdkey, byte=isbyte
328     # ifdef ALLOW_SEAICE
329     CADJ STORE surfacetendencyTice(:,:,bi,bj)
330     CADJ & = comlev1_bibj, key=itdkey, byte=isbyte
331     # endif
332     #endif /* ALLOW_AUTODIFF_TAMC */
333    
334     #ifdef ALLOW_GMREDI
335    
336     #ifdef ALLOW_AUTODIFF_TAMC
337     cph storing here is needed only for one GMREDI_OPTIONS:
338     cph define GM_BOLUS_ADVEC
339     cph but I've avoided the #ifdef for now, in case more things change
340     CADJ STORE sigmaX(:,:,:) = comlev1_bibj, key=itdkey, byte=isbyte
341     CADJ STORE sigmaY(:,:,:) = comlev1_bibj, key=itdkey, byte=isbyte
342     CADJ STORE sigmaR(:,:,:) = comlev1_bibj, key=itdkey, byte=isbyte
343     #endif /* ALLOW_AUTODIFF_TAMC */
344    
345     C-- Calculate iso-neutral slopes for the GM/Redi parameterisation
346     IF (useGMRedi) THEN
347     #ifdef ALLOW_DEBUG
348     IF ( debugLevel .GE. debLevB )
349     & CALL DEBUG_CALL('GMREDI_CALC_TENSOR',myThid)
350     #endif
351     CALL GMREDI_CALC_TENSOR(
352     I bi, bj, iMin, iMax, jMin, jMax,
353     I sigmaX, sigmaY, sigmaR,
354     I myThid )
355     #ifdef ALLOW_AUTODIFF_TAMC
356     ELSE
357     CALL GMREDI_CALC_TENSOR_DUMMY(
358     I bi, bj, iMin, iMax, jMin, jMax,
359     I sigmaX, sigmaY, sigmaR,
360     I myThid )
361     #endif /* ALLOW_AUTODIFF_TAMC */
362     ENDIF
363    
364     #ifdef ALLOW_AUTODIFF_TAMC
365     CADJ STORE Kwx(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
366     CADJ STORE Kwy(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
367     CADJ STORE Kwz(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
368     #endif /* ALLOW_AUTODIFF_TAMC */
369    
370     #endif /* ALLOW_GMREDI */
371    
372     #ifdef ALLOW_KPP
373     C-- Compute KPP mixing coefficients
374     IF (useKPP) THEN
375     #ifdef ALLOW_DEBUG
376     IF ( debugLevel .GE. debLevB )
377     & CALL DEBUG_CALL('KPP_CALC',myThid)
378     #endif
379     CALL KPP_CALC(
380     I bi, bj, myTime, myThid )
381     #ifdef ALLOW_AUTODIFF_TAMC
382     ELSE
383     CALL KPP_CALC_DUMMY(
384     I bi, bj, myTime, myThid )
385     #endif /* ALLOW_AUTODIFF_TAMC */
386     ENDIF
387    
388     #ifdef ALLOW_AUTODIFF_TAMC
389     CADJ STORE KPPghat (:,:,:,bi,bj)
390     CADJ & , KPPfrac (:,: ,bi,bj)
391     CADJ & = comlev1_bibj, key=itdkey, byte=isbyte
392     #endif /* ALLOW_AUTODIFF_TAMC */
393    
394     #endif /* ALLOW_KPP */
395    
396     #ifdef ALLOW_AUTODIFF_TAMC
397     CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
398     CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
399     CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
400     CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
401     #ifdef ALLOW_PASSIVE_TRACER
402     CADJ STORE tr1 (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
403     #endif
404     #ifdef ALLOW_PTRACERS
405     cph-- moved to forward_step to avoid key computation
406     cphCADJ STORE ptracer(:,:,:,bi,bj,itracer) = comlev1_bibj,
407     cphCADJ & key=itdkey, byte=isbyte
408     #endif
409     #endif /* ALLOW_AUTODIFF_TAMC */
410    
411     #ifdef ALLOW_AIM
412     C AIM - atmospheric intermediate model, physics package code.
413     IF ( useAIM ) THEN
414     #ifdef ALLOW_DEBUG
415     IF ( debugLevel .GE. debLevB )
416     & CALL DEBUG_CALL('AIM_DO_PHYSICS',myThid)
417     #endif
418     CALL TIMER_START('AIM_DO_PHYSICS [THERMODYNAMICS]', myThid)
419     CALL AIM_DO_PHYSICS( bi, bj, myTime, myIter, myThid )
420     CALL TIMER_STOP( 'AIM_DO_PHYSICS [THERMODYNAMICS]', myThid)
421     ENDIF
422     #endif /* ALLOW_AIM */
423    
424     C-- end bi,bj loops.
425     ENDDO
426     ENDDO
427    
428     #ifdef ALLOW_AIM
429     IF ( useAIM ) THEN
430     CALL AIM_AIM2DYN_EXCHANGES( myTime, myIter, myThid )
431     ENDIF
432     #endif /* ALLOW_AIM */
433    
434     #ifdef ALLOW_DEBUG
435     IF ( debugLevel .GE. debLevB )
436     & CALL DEBUG_LEAVE('DO_OCEANIC_PHYS',myThid)
437     #endif
438    
439     RETURN
440     END

  ViewVC Help
Powered by ViewVC 1.1.22