/[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.8 - (hide annotations) (download)
Wed Sep 8 01:50:24 2004 UTC (19 years, 8 months ago) by jmc
Branch: MAIN
Changes since 1.7: +9 -1 lines
add diagnostics od d.Rho/dr (=SigmaR).

1 jmc 1.8 C $Header: /u/gcmpack/MITgcm/model/src/do_oceanic_phys.F,v 1.7 2004/09/07 17:29:14 edhill 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 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_AUTODIFF_TAMC
324     cph avoids recomputation of integrate_for_w
325     CADJ STORE wvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
326     #endif /* ALLOW_AUTODIFF_TAMC */
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 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     # ifdef ALLOW_SEAICE
368 jmc 1.4 CADJ STORE surfaceForcingTice(:,:,bi,bj)
369 jmc 1.1 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     cph storing here is needed only for one GMREDI_OPTIONS:
377     cph define GM_BOLUS_ADVEC
378     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     CADJ STORE sigmaR(:,:,:) = comlev1_bibj, key=itdkey, byte=isbyte
382     #endif /* ALLOW_AUTODIFF_TAMC */
383    
384     C-- Calculate iso-neutral slopes for the GM/Redi parameterisation
385     IF (useGMRedi) THEN
386     #ifdef ALLOW_DEBUG
387     IF ( debugLevel .GE. debLevB )
388     & CALL DEBUG_CALL('GMREDI_CALC_TENSOR',myThid)
389     #endif
390     CALL GMREDI_CALC_TENSOR(
391     I bi, bj, iMin, iMax, jMin, jMax,
392     I sigmaX, sigmaY, sigmaR,
393     I myThid )
394     #ifdef ALLOW_AUTODIFF_TAMC
395     ELSE
396     CALL GMREDI_CALC_TENSOR_DUMMY(
397     I bi, bj, iMin, iMax, jMin, jMax,
398     I sigmaX, sigmaY, sigmaR,
399     I myThid )
400     #endif /* ALLOW_AUTODIFF_TAMC */
401     ENDIF
402    
403     #ifdef ALLOW_AUTODIFF_TAMC
404     CADJ STORE Kwx(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
405     CADJ STORE Kwy(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
406     CADJ STORE Kwz(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
407     #endif /* ALLOW_AUTODIFF_TAMC */
408    
409     #endif /* ALLOW_GMREDI */
410    
411     #ifdef ALLOW_KPP
412     C-- Compute KPP mixing coefficients
413     IF (useKPP) THEN
414     #ifdef ALLOW_DEBUG
415     IF ( debugLevel .GE. debLevB )
416     & CALL DEBUG_CALL('KPP_CALC',myThid)
417     #endif
418     CALL KPP_CALC(
419     I bi, bj, myTime, myThid )
420     #ifdef ALLOW_AUTODIFF_TAMC
421     ELSE
422     CALL KPP_CALC_DUMMY(
423     I bi, bj, myTime, myThid )
424     #endif /* ALLOW_AUTODIFF_TAMC */
425     ENDIF
426    
427     #ifdef ALLOW_AUTODIFF_TAMC
428     CADJ STORE KPPghat (:,:,:,bi,bj)
429     CADJ & , KPPfrac (:,: ,bi,bj)
430     CADJ & = comlev1_bibj, key=itdkey, byte=isbyte
431     #endif /* ALLOW_AUTODIFF_TAMC */
432    
433     #endif /* ALLOW_KPP */
434    
435 mlosch 1.6 #ifdef ALLOW_PP81
436     C-- Compute PP81 mixing coefficients
437     IF (usePP81) THEN
438     #ifdef ALLOW_DEBUG
439     IF ( debugLevel .GE. debLevB )
440     & CALL DEBUG_CALL('PP81_CALC',myThid)
441     #endif
442     CALL PP81_CALC(
443     I bi, bj, myTime, myThid )
444     ENDIF
445     #endif /* ALLOW_PP81 */
446    
447     #ifdef ALLOW_MY82
448     C-- Compute MY82 mixing coefficients
449     IF (useMY82) THEN
450     #ifdef ALLOW_DEBUG
451     IF ( debugLevel .GE. debLevB )
452     & CALL DEBUG_CALL('MY82_CALC',myThid)
453     #endif
454     CALL MY82_CALC(
455     I bi, bj, myTime, myThid )
456     ENDIF
457     #endif /* ALLOW_MY82 */
458    
459 jmc 1.1 #ifdef ALLOW_AUTODIFF_TAMC
460     CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
461     CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
462     CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
463     CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
464     #ifdef ALLOW_PASSIVE_TRACER
465     CADJ STORE tr1 (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
466     #endif
467     #ifdef ALLOW_PTRACERS
468     cph-- moved to forward_step to avoid key computation
469     cphCADJ STORE ptracer(:,:,:,bi,bj,itracer) = comlev1_bibj,
470     cphCADJ & key=itdkey, byte=isbyte
471     #endif
472     #endif /* ALLOW_AUTODIFF_TAMC */
473    
474     C-- end bi,bj loops.
475     ENDDO
476     ENDDO
477    
478     #ifdef ALLOW_DEBUG
479     IF ( debugLevel .GE. debLevB )
480     & CALL DEBUG_LEAVE('DO_OCEANIC_PHYS',myThid)
481     #endif
482    
483     RETURN
484     END

  ViewVC Help
Powered by ViewVC 1.1.22