/[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.10 - (hide annotations) (download)
Fri Sep 17 23:02:00 2004 UTC (19 years, 8 months ago) by heimbach
Branch: MAIN
CVS Tags: checkpoint55c_post, checkpoint55b_post, checkpoint55, checkpoint55a_post
Changes since 1.9: +15 -39 lines
o bringing adjoint up to date for sheduled c55

1 heimbach 1.10 C $Header: /u/gcmpack/MITgcm/model/src/do_oceanic_phys.F,v 1.9 2004/09/16 09:35:11 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_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     #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 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.10 # ifndef GM_EXCLUDE_CLIPPING
300 jmc 1.1 CADJ STORE rhok (:,:) = comlev1_bibj_k , key=kkey, byte=isbyte
301 heimbach 1.10 # 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     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 heimbach 1.10 # ifndef GM_EXCLUDE_CLIPPING
377 jmc 1.1 cph storing here is needed only for one GMREDI_OPTIONS:
378     cph define GM_BOLUS_ADVEC
379 heimbach 1.10 cph keep it although TAF says you dont need to.
380 jmc 1.1 cph but I've avoided the #ifdef for now, in case more things change
381     CADJ STORE sigmaX(:,:,:) = comlev1_bibj, key=itdkey, byte=isbyte
382     CADJ STORE sigmaY(:,:,:) = comlev1_bibj, key=itdkey, byte=isbyte
383 heimbach 1.10 cnewCADJ STORE sigmaR(:,:,:) = comlev1_bibj, key=itdkey, byte=isbyte
384     # endif
385 jmc 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
386    
387     C-- Calculate iso-neutral slopes for the GM/Redi parameterisation
388     IF (useGMRedi) THEN
389     #ifdef ALLOW_DEBUG
390     IF ( debugLevel .GE. debLevB )
391     & CALL DEBUG_CALL('GMREDI_CALC_TENSOR',myThid)
392     #endif
393     CALL GMREDI_CALC_TENSOR(
394     I bi, bj, iMin, iMax, jMin, jMax,
395     I sigmaX, sigmaY, sigmaR,
396     I myThid )
397     #ifdef ALLOW_AUTODIFF_TAMC
398     ELSE
399     CALL GMREDI_CALC_TENSOR_DUMMY(
400     I bi, bj, iMin, iMax, jMin, jMax,
401     I sigmaX, sigmaY, sigmaR,
402     I myThid )
403     #endif /* ALLOW_AUTODIFF_TAMC */
404     ENDIF
405    
406     #endif /* ALLOW_GMREDI */
407    
408     #ifdef ALLOW_KPP
409     C-- Compute KPP mixing coefficients
410     IF (useKPP) THEN
411     #ifdef ALLOW_DEBUG
412     IF ( debugLevel .GE. debLevB )
413     & CALL DEBUG_CALL('KPP_CALC',myThid)
414     #endif
415     CALL KPP_CALC(
416     I bi, bj, myTime, myThid )
417     #ifdef ALLOW_AUTODIFF_TAMC
418     ELSE
419     CALL KPP_CALC_DUMMY(
420     I bi, bj, myTime, myThid )
421     #endif /* ALLOW_AUTODIFF_TAMC */
422     ENDIF
423    
424     #endif /* ALLOW_KPP */
425    
426 mlosch 1.6 #ifdef ALLOW_PP81
427     C-- Compute PP81 mixing coefficients
428     IF (usePP81) THEN
429     #ifdef ALLOW_DEBUG
430     IF ( debugLevel .GE. debLevB )
431     & CALL DEBUG_CALL('PP81_CALC',myThid)
432     #endif
433     CALL PP81_CALC(
434     I bi, bj, myTime, myThid )
435     ENDIF
436     #endif /* ALLOW_PP81 */
437    
438     #ifdef ALLOW_MY82
439     C-- Compute MY82 mixing coefficients
440     IF (useMY82) THEN
441     #ifdef ALLOW_DEBUG
442     IF ( debugLevel .GE. debLevB )
443     & CALL DEBUG_CALL('MY82_CALC',myThid)
444     #endif
445     CALL MY82_CALC(
446     I bi, bj, myTime, myThid )
447     ENDIF
448     #endif /* ALLOW_MY82 */
449    
450 mlosch 1.9 #ifdef ALLOW_GGL90
451     C-- Compute GGL90 mixing coefficients
452     IF (useGGL90) THEN
453     #ifdef ALLOW_DEBUG
454     IF ( debugLevel .GE. debLevB )
455     & CALL DEBUG_CALL('GGL90_CALC',myThid)
456     #endif
457     CALL GGL90_CALC(
458     I bi, bj, myTime, myThid )
459     ENDIF
460     #endif /* ALLOW_GGL90 */
461    
462 jmc 1.1 C-- end bi,bj loops.
463     ENDDO
464     ENDDO
465    
466     #ifdef ALLOW_DEBUG
467     IF ( debugLevel .GE. debLevB )
468     & CALL DEBUG_LEAVE('DO_OCEANIC_PHYS',myThid)
469     #endif
470    
471     RETURN
472     END

  ViewVC Help
Powered by ViewVC 1.1.22