/[MITgcm]/MITgcm_contrib/atnguyen/code_21Dec2012_saltplume/do_oceanic_phys.F
ViewVC logotype

Annotation of /MITgcm_contrib/atnguyen/code_21Dec2012_saltplume/do_oceanic_phys.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph

Revision 1.1 - (hide annotations) (download)
Tue Apr 22 04:32:25 2014 UTC (11 years, 3 months ago) by atn
Branch: MAIN 
in progress:
1. sort out bi,bj loop
2. skip remove saltplumeflux in external_forcing_surf, (thus) skip kpp
3. move saltplume outside k-loop in [salt,temp]_integrate.F
4. add diagnostics
1 atn 1.1 C $Header: /u/gcmpack/MITgcm/model/src/do_oceanic_phys.F,v 1.134 2014/04/04 20:54:11 jmc Exp $
2     C $Name: $
3    
4     #include "PACKAGES_CONFIG.h"
5     #include "CPP_OPTIONS.h"
6     #ifdef ALLOW_AUTODIFF
7     # include "AUTODIFF_OPTIONS.h"
8     #endif
9     #ifdef ALLOW_CTRL
10     # include "CTRL_OPTIONS.h"
11     #endif
12    
13     #ifdef ALLOW_AUTODIFF
14     # ifdef ALLOW_GMREDI
15     # include "GMREDI_OPTIONS.h"
16     # endif
17     # ifdef ALLOW_KPP
18     # include "KPP_OPTIONS.h"
19     # endif
20     # ifdef ALLOW_SEAICE
21     # include "SEAICE_OPTIONS.h"
22     # endif
23     # ifdef ALLOW_EXF
24     # include "EXF_OPTIONS.h"
25     # endif
26     #endif /* ALLOW_AUTODIFF */
27    
28     CBOP
29     C !ROUTINE: DO_OCEANIC_PHYS
30     C !INTERFACE:
31     SUBROUTINE DO_OCEANIC_PHYS(myTime, myIter, myThid)
32     C !DESCRIPTION: \bv
33     C *==========================================================*
34     C | SUBROUTINE DO_OCEANIC_PHYS
35     C | o Controlling routine for oceanic physics and
36     C | parameterization
37     C *==========================================================*
38     C | o originally, part of S/R thermodynamics
39     C *==========================================================*
40     C \ev
41    
42     C !USES:
43     IMPLICIT NONE
44     C == Global variables ===
45     #include "SIZE.h"
46     #include "EEPARAMS.h"
47     #include "PARAMS.h"
48     #include "GRID.h"
49     #include "DYNVARS.h"
50     #ifdef ALLOW_TIMEAVE
51     # include "TIMEAVE_STATV.h"
52     #endif
53     #ifdef ALLOW_OFFLINE
54     # include "OFFLINE_SWITCH.h"
55     #endif
56    
57     #ifdef ALLOW_AUTODIFF
58     # include "AUTODIFF_MYFIELDS.h"
59     # include "tamc.h"
60     # include "tamc_keys.h"
61     # include "FFIELDS.h"
62     # include "SURFACE.h"
63     # include "EOS.h"
64     # ifdef ALLOW_KPP
65     # include "KPP.h"
66     # endif
67     # ifdef ALLOW_GGL90
68     # include "GGL90.h"
69     # endif
70     # ifdef ALLOW_GMREDI
71     # include "GMREDI.h"
72     # endif
73     # ifdef ALLOW_EBM
74     # include "EBM.h"
75     # endif
76     # ifdef ALLOW_EXF
77     # include "ctrl.h"
78     # include "EXF_FIELDS.h"
79     # ifdef ALLOW_BULKFORMULAE
80     # include "EXF_CONSTANTS.h"
81     # endif
82     # endif
83     # ifdef ALLOW_SEAICE
84     # include "SEAICE_SIZE.h"
85     # include "SEAICE.h"
86     # include "SEAICE_PARAMS.h"
87     # endif
88     # ifdef ALLOW_THSICE
89     # include "THSICE_VARS.h"
90     # endif
91     # ifdef ALLOW_SALT_PLUME
92     # include "SALT_PLUME.h"
93     # endif
94     #endif /* ALLOW_AUTODIFF */
95    
96     C !INPUT/OUTPUT PARAMETERS:
97     C == Routine arguments ==
98     C myTime :: Current time in simulation
99     C myIter :: Current iteration number in simulation
100     C myThid :: Thread number for this instance of the routine.
101     _RL myTime
102     INTEGER myIter
103     INTEGER myThid
104    
105     C !LOCAL VARIABLES:
106     C == Local variables
107     C rhoK, rhoKm1 :: Density at current level, and level above
108     C iMin, iMax :: Ranges and sub-block indices on which calculations
109     C jMin, jMax are applied.
110     C bi, bj :: tile indices
111     C i,j,k :: loop indices
112     _RL rhoKp1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
113     _RL rhoKm1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
114     _RL sigmaX (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
115     _RL sigmaY (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
116     _RL sigmaR (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
117     INTEGER iMin, iMax
118     INTEGER jMin, jMax
119     INTEGER bi, bj
120     INTEGER i, j, k
121     INTEGER doDiagsRho
122     LOGICAL calcGMRedi, calcKPP, calcConvect
123     #ifdef ALLOW_DIAGNOSTICS
124     LOGICAL DIAGNOSTICS_IS_ON
125     EXTERNAL DIAGNOSTICS_IS_ON
126     #endif /* ALLOW_DIAGNOSTICS */
127    
128     CEOP
129    
130     #ifdef ALLOW_AUTODIFF_TAMC
131     C-- dummy statement to end declaration part
132     itdkey = 1
133     #endif /* ALLOW_AUTODIFF_TAMC */
134    
135     #ifdef ALLOW_DEBUG
136     IF (debugMode) CALL DEBUG_ENTER('DO_OCEANIC_PHYS',myThid)
137     #endif
138    
139     doDiagsRho = 0
140     #ifdef ALLOW_DIAGNOSTICS
141     IF ( useDiagnostics .AND. fluidIsWater ) THEN
142     IF ( DIAGNOSTICS_IS_ON('MXLDEPTH',myThid) )
143     & doDiagsRho = doDiagsRho + 1
144     IF ( DIAGNOSTICS_IS_ON('DRHODR ',myThid) )
145     & doDiagsRho = doDiagsRho + 2
146     IF ( DIAGNOSTICS_IS_ON('WdRHO_P ',myThid) )
147     & doDiagsRho = doDiagsRho + 4
148     IF ( DIAGNOSTICS_IS_ON('WdRHOdP ',myThid) )
149     & doDiagsRho = doDiagsRho + 8
150     ENDIF
151     #endif /* ALLOW_DIAGNOSTICS */
152    
153     calcGMRedi = useGMRedi
154     calcKPP = useKPP
155     calcConvect = ivdc_kappa.NE.0.
156     #ifdef ALLOW_OFFLINE
157     IF ( useOffLine ) THEN
158     calcGMRedi = useGMRedi .AND. .NOT.offlineLoadGMRedi
159     calcKPP = useKPP .AND. .NOT.offlineLoadKPP
160     calcConvect=calcConvect.AND. .NOT.offlineLoadConvec
161     ENDIF
162     #endif /* ALLOW_OFFLINE */
163    
164     #ifdef ALLOW_OBCS
165     IF (useOBCS) THEN
166     C-- Calculate future values on open boundaries
167     C-- moved before SEAICE_MODEL call since SEAICE_MODEL needs seaice-obcs fields
168     # ifdef ALLOW_AUTODIFF_TAMC
169     CADJ STORE theta = comlev1, key=ikey_dynamics, kind=isbyte
170     CADJ STORE salt = comlev1, key=ikey_dynamics, kind=isbyte
171     # endif
172     # ifdef ALLOW_DEBUG
173     IF (debugMode) CALL DEBUG_CALL('OBCS_CALC',myThid)
174     # endif
175     CALL OBCS_CALC( myTime+deltaTClock, myIter+1,
176     I uVel, vVel, wVel, theta, salt, myThid )
177     ENDIF
178     #endif /* ALLOW_OBCS */
179    
180     #ifdef ALLOW_AUTODIFF
181     # ifdef ALLOW_SALT_PLUME
182     DO bj=myByLo(myThid),myByHi(myThid)
183     DO bi=myBxLo(myThid),myBxHi(myThid)
184     DO j=1-OLy,sNy+OLy
185     DO i=1-OLx,sNx+OLx
186     saltPlumeDepth(i,j,bi,bj) = 0. _d 0
187     saltPlumeFlux(i,j,bi,bj) = 0. _d 0
188     ENDDO
189     ENDDO
190     ENDDO
191     ENDDO
192     # endif
193     #endif /* ALLOW_AUTODIFF */
194    
195     #ifdef ALLOW_FRAZIL
196     IF ( useFRAZIL ) THEN
197     C-- Freeze water in the ocean interior and let it rise to the surface
198     CALL FRAZIL_CALC_RHS( myTime, myIter, myThid )
199     ENDIF
200     #endif /* ALLOW_FRAZIL */
201    
202     #ifndef OLD_THSICE_CALL_SEQUENCE
203     #if (defined ALLOW_THSICE) && !(defined ALLOW_ATM2D)
204     IF ( useThSIce .AND. fluidIsWater ) THEN
205     # ifdef ALLOW_AUTODIFF_TAMC
206     CADJ STORE uice,vice = comlev1, key = ikey_dynamics,
207     CADJ & kind = isbyte
208     CADJ STORE iceMask,iceHeight = comlev1, key = ikey_dynamics,
209     CADJ & kind = isbyte
210     CADJ STORE snowHeight, Tsrf = comlev1, key = ikey_dynamics,
211     CADJ & kind = isbyte
212     CADJ STORE Qice1, Qice2 = comlev1, key = ikey_dynamics,
213     CADJ & kind = isbyte
214     CADJ STORE sHeating, snowAge = comlev1, key = ikey_dynamics,
215     CADJ & kind = isbyte
216     CADJ STORE hocemxl = comlev1, key = ikey_dynamics,
217     CADJ & kind = isbyte
218     CADJ STORE icflxsw = comlev1, key = ikey_dynamics,
219     CADJ & kind = isbyte
220     CADJ STORE salt,theta = comlev1, key = ikey_dynamics,
221     CADJ & kind = isbyte
222     CADJ STORE uvel,vvel = comlev1, key = ikey_dynamics,
223     CADJ & kind = isbyte
224     CADJ STORE qnet,qsw, empmr = comlev1, key = ikey_dynamics,
225     CADJ & kind = isbyte
226     CADJ STORE atemp,aqh,precip = comlev1, key = ikey_dynamics,
227     CADJ & kind = isbyte
228     CADJ STORE swdown,lwdown = comlev1, key = ikey_dynamics,
229     CADJ & kind = isbyte
230     # ifdef NONLIN_FRSURF
231     CADJ STORE hFac_surfC = comlev1, key = ikey_dynamics,
232     CADJ & kind = isbyte
233     # endif
234     # endif /* ALLOW_AUTODIFF_TAMC */
235     # ifdef ALLOW_DEBUG
236     IF (debugMode) CALL DEBUG_CALL('THSICE_MAIN',myThid)
237     # endif
238     C-- Step forward Therm.Sea-Ice variables
239     C and modify forcing terms including effects from ice
240     CALL TIMER_START('THSICE_MAIN [DO_OCEANIC_PHYS]', myThid)
241     CALL THSICE_MAIN( myTime, myIter, myThid )
242     CALL TIMER_STOP( 'THSICE_MAIN [DO_OCEANIC_PHYS]', myThid)
243     ENDIF
244     #endif /* ALLOW_THSICE */
245     #endif /* ndef OLD_THSICE_CALL_SEQUENCE */
246    
247     #ifdef ALLOW_SEAICE
248     # ifdef ALLOW_AUTODIFF
249     CADJ STORE area = comlev1, key=ikey_dynamics, kind=isbyte
250     CADJ STORE fu,fv = comlev1, key=ikey_dynamics, kind=isbyte
251     CADJ STORE qnet = comlev1, key=ikey_dynamics, kind=isbyte
252     CADJ STORE qsw = comlev1, key=ikey_dynamics, kind=isbyte
253     CADJ STORE theta = comlev1, key=ikey_dynamics, kind=isbyte
254     CADJ STORE salt = comlev1, key=ikey_dynamics, kind=isbyte
255     #if (defined ALLOW_EXF) && (defined ALLOW_ATM_TEMP)
256     CADJ STORE evap = comlev1, key=ikey_dynamics, kind=isbyte
257     #endif
258     IF ( .NOT.useSEAICE .AND. SEAICEadjMODE .EQ. -1 ) THEN
259     CALL SEAICE_FAKE( myTime, myIter, myThid )
260     ENDIF
261     CADJ STORE area = comlev1, key=ikey_dynamics, kind=isbyte
262     CADJ STORE fu,fv = comlev1, key=ikey_dynamics, kind=isbyte
263     CADJ STORE qnet = comlev1, key=ikey_dynamics, kind=isbyte
264     CADJ STORE qsw = comlev1, key=ikey_dynamics, kind=isbyte
265     CADJ STORE theta = comlev1, key=ikey_dynamics, kind=isbyte
266     CADJ STORE salt = comlev1, key=ikey_dynamics, kind=isbyte
267     #if (defined ALLOW_EXF) && (defined ALLOW_ATM_TEMP)
268     CADJ STORE evap = comlev1, key=ikey_dynamics, kind=isbyte
269     #endif
270     # endif /* ALLOW_AUTODIFF */
271     #endif /* ALLOW_SEAICE */
272    
273     #ifdef ALLOW_SEAICE
274     IF ( useSEAICE ) THEN
275     # ifdef ALLOW_AUTODIFF_TAMC
276     cph-adj-test(
277     CADJ STORE area = comlev1, key=ikey_dynamics, kind=isbyte
278     CADJ STORE hsnow = comlev1, key=ikey_dynamics, kind=isbyte
279     CADJ STORE heff = comlev1, key=ikey_dynamics, kind=isbyte
280     CADJ STORE tices = comlev1, key=ikey_dynamics, kind=isbyte
281     CADJ STORE empmr, qnet = comlev1, key=ikey_dynamics, kind=isbyte
282     CADJ STORE qsw,saltflux = comlev1, key=ikey_dynamics, kind=isbyte
283     CADJ STORE fu, fv = comlev1, key=ikey_dynamics, kind=isbyte
284     cCADJ STORE theta = comlev1, key=ikey_dynamics, kind=isbyte
285     cCADJ STORE salt = comlev1, key=ikey_dynamics, kind=isbyte
286     cph-adj-test)
287     c#ifdef ALLOW_EXF
288     CADJ STORE atemp,aqh,precip = comlev1, key = ikey_dynamics,
289     CADJ & kind = isbyte
290     CADJ STORE swdown,lwdown = comlev1, key = ikey_dynamics,
291     CADJ & kind = isbyte
292     CADJ STORE evap = comlev1, key = ikey_dynamics,
293     CADJ & kind = isbyte
294     CADJ STORE uwind,vwind = comlev1, key = ikey_dynamics,
295     CADJ & kind = isbyte
296     c#endif
297     CADJ STORE uvel,vvel = comlev1, key = ikey_dynamics,
298     CADJ & kind = isbyte
299     # ifdef SEAICE_CGRID
300     CADJ STORE stressdivergencex = comlev1, key = ikey_dynamics,
301     CADJ & kind = isbyte
302     CADJ STORE stressdivergencey = comlev1, key = ikey_dynamics,
303     CADJ & kind = isbyte
304     # endif
305     # ifdef SEAICE_ALLOW_DYNAMICS
306     CADJ STORE uice = comlev1, key = ikey_dynamics,
307     CADJ & kind = isbyte
308     CADJ STORE vice = comlev1, key = ikey_dynamics,
309     CADJ & kind = isbyte
310     CADJ STORE dwatn = comlev1, key = ikey_dynamics,
311     CADJ & kind = isbyte
312     # ifdef SEAICE_ALLOW_EVP
313     CADJ STORE seaice_sigma1 = comlev1, key = ikey_dynamics,
314     CADJ & kind = isbyte
315     CADJ STORE seaice_sigma2 = comlev1, key = ikey_dynamics,
316     CADJ & kind = isbyte
317     CADJ STORE seaice_sigma12 = comlev1, key = ikey_dynamics,
318     CADJ & kind = isbyte
319     # endif
320     # endif
321     # ifdef SEAICE_VARIABLE_SALINITY
322     CADJ STORE hsalt = comlev1, key = ikey_dynamics,
323     CADJ & kind = isbyte
324     # endif
325     # ifdef ATMOSPHERIC_LOADING
326     CADJ STORE pload = comlev1, key = ikey_dynamics,
327     CADJ & kind = isbyte
328     CADJ STORE siceload = comlev1, key = ikey_dynamics,
329     CADJ & kind = isbyte
330     # endif
331     # ifdef NONLIN_FRSURF
332     CADJ STORE recip_hfacc = comlev1, key = ikey_dynamics,
333     CADJ & kind = isbyte
334     # endif
335     # ifdef ANNUAL_BALANCE
336     CADJ STORE balance_itcount = comlev1, key = ikey_dynamics,
337     CADJ & kind = isbyte
338     # endif /* ANNUAL_BALANCE */
339     # ifdef ALLOW_THSICE
340     C-- store thSIce vars before advection (called from SEAICE_MODEL) update them:
341     CADJ STORE iceMask,iceHeight = comlev1, key = ikey_dynamics,
342     CADJ & kind = isbyte
343     CADJ STORE snowHeight, hOceMxL = comlev1, key = ikey_dynamics,
344     CADJ & kind = isbyte
345     CADJ STORE Qice1, Qice2 = comlev1, key = ikey_dynamics,
346     CADJ & kind = isbyte
347     # endif /* ALLOW_THSICE */
348     # endif /* ALLOW_AUTODIFF_TAMC */
349     # ifdef ALLOW_DEBUG
350     IF (debugMode) CALL DEBUG_CALL('SEAICE_MODEL',myThid)
351     # endif
352     CALL TIMER_START('SEAICE_MODEL [DO_OCEANIC_PHYS]', myThid)
353     CALL SEAICE_MODEL( myTime, myIter, myThid )
354     CALL TIMER_STOP ('SEAICE_MODEL [DO_OCEANIC_PHYS]', myThid)
355     # ifdef ALLOW_COST
356     CALL SEAICE_COST_SENSI ( myTime, myIter, myThid )
357     # endif
358     ENDIF
359     #endif /* ALLOW_SEAICE */
360    
361     #ifdef ALLOW_AUTODIFF_TAMC
362     CADJ STORE sst, sss = comlev1, key = ikey_dynamics,
363     CADJ & kind = isbyte
364     CADJ STORE qsw = comlev1, key = ikey_dynamics,
365     CADJ & kind = isbyte
366     # ifdef ALLOW_SEAICE
367     CADJ STORE area = comlev1, key = ikey_dynamics,
368     CADJ & kind = isbyte
369     # endif
370     #endif
371    
372     #ifdef OLD_THSICE_CALL_SEQUENCE
373     #if (defined ALLOW_THSICE) && !(defined ALLOW_ATM2D)
374     IF ( useThSIce .AND. fluidIsWater ) THEN
375     # ifdef ALLOW_AUTODIFF_TAMC
376     cph(
377     # ifdef NONLIN_FRSURF
378     CADJ STORE uice,vice = comlev1, key = ikey_dynamics,
379     CADJ & kind = isbyte
380     CADJ STORE salt,theta = comlev1, key = ikey_dynamics,
381     CADJ & kind = isbyte
382     CADJ STORE qnet,qsw, empmr = comlev1, key = ikey_dynamics,
383     CADJ & kind = isbyte
384     CADJ STORE hFac_surfC = comlev1, key = ikey_dynamics,
385     CADJ & kind = isbyte
386     # endif
387     # endif
388     # ifdef ALLOW_DEBUG
389     IF (debugMode) CALL DEBUG_CALL('THSICE_MAIN',myThid)
390     # endif
391     C-- Step forward Therm.Sea-Ice variables
392     C and modify forcing terms including effects from ice
393     CALL TIMER_START('THSICE_MAIN [DO_OCEANIC_PHYS]', myThid)
394     CALL THSICE_MAIN( myTime, myIter, myThid )
395     CALL TIMER_STOP( 'THSICE_MAIN [DO_OCEANIC_PHYS]', myThid)
396     ENDIF
397     #endif /* ALLOW_THSICE */
398     #endif /* OLD_THSICE_CALL_SEQUENCE */
399    
400     #ifdef ALLOW_SHELFICE
401     # ifdef ALLOW_AUTODIFF_TAMC
402     CADJ STORE salt, theta = comlev1, key = ikey_dynamics,
403     CADJ & kind = isbyte
404     # endif
405     IF ( useShelfIce .AND. fluidIsWater ) THEN
406     #ifdef ALLOW_DEBUG
407     IF (debugMode) CALL DEBUG_CALL('SHELFICE_THERMODYNAMICS',myThid)
408     #endif
409     C compute temperature and (virtual) salt flux at the
410     C shelf-ice ocean interface
411     CALL TIMER_START('SHELFICE_THERMODYNAMICS [DO_OCEANIC_PHYS]',
412     & myThid)
413     CALL SHELFICE_THERMODYNAMICS( myTime, myIter, myThid )
414     CALL TIMER_STOP( 'SHELFICE_THERMODYNAMICS [DO_OCEANIC_PHYS]',
415     & myThid)
416     ENDIF
417     #endif /* ALLOW_SHELFICE */
418    
419     #ifdef ALLOW_ICEFRONT
420     IF ( useICEFRONT .AND. fluidIsWater ) THEN
421     #ifdef ALLOW_DEBUG
422     IF (debugMode) CALL DEBUG_CALL('ICEFRONT_THERMODYNAMICS',myThid)
423     #endif
424     C compute temperature and (virtual) salt flux at the
425     C ice-front ocean interface
426     CALL TIMER_START('ICEFRONT_THERMODYNAMICS [DO_OCEANIC_PHYS]',
427     & myThid)
428     CALL ICEFRONT_THERMODYNAMICS( myTime, myIter, myThid )
429     CALL TIMER_STOP( 'ICEFRONT_THERMODYNAMICS [DO_OCEANIC_PHYS]',
430     & myThid)
431     ENDIF
432     #endif /* ALLOW_ICEFRONT */
433    
434     #ifdef ALLOW_SALT_PLUME
435     IF ( useSALT_PLUME ) THEN
436     Catn: exchanging saltPlumeFlux:
437     CALL SALT_PLUME_DO_EXCH( myTime, myIter, myThid )
438     ENDIF
439     #endif /* ALLOW_SALT_PLUME */
440    
441     C-- Freeze water at the surface
442     IF ( allowFreezing ) THEN
443     #ifdef ALLOW_AUTODIFF_TAMC
444     CADJ STORE theta = comlev1, key = ikey_dynamics,
445     CADJ & kind = isbyte
446     #endif
447     CALL FREEZE_SURFACE( myTime, myIter, myThid )
448     ENDIF
449    
450     #ifdef ALLOW_OCN_COMPON_INTERF
451     C-- Apply imported data (from coupled interface) to forcing fields
452     C jmc: do not know precisely where to put this call (bf or af thSIce ?)
453     IF ( useCoupler ) THEN
454     CALL OCN_APPLY_IMPORT( .TRUE., myTime, myIter, myThid )
455     ENDIF
456     #endif /* ALLOW_OCN_COMPON_INTERF */
457    
458     iMin = 1-OLx
459     iMax = sNx+OLx
460     jMin = 1-OLy
461     jMax = sNy+OLy
462    
463     C--- Determines forcing terms based on external fields
464     C relaxation terms, etc.
465     #ifdef ALLOW_DEBUG
466     IF (debugMode) CALL DEBUG_CALL('EXTERNAL_FORCING_SURF',myThid)
467     #endif
468     #ifdef ALLOW_AUTODIFF
469     CADJ STORE salt, theta = comlev1, key = ikey_dynamics,
470     CADJ & kind = isbyte
471     #else /* ALLOW_AUTODIFF */
472     C-- if fluid is not water, by-pass surfaceForcing, find_rho, gmredi
473     C and all vertical mixing schemes, but keep OBCS_CALC
474     IF ( fluidIsWater ) THEN
475     #endif /* ALLOW_AUTODIFF */
476     CALL EXTERNAL_FORCING_SURF(
477     I iMin, iMax, jMin, jMax,
478     I myTime, myIter, myThid )
479    
480     #ifdef ALLOW_AUTODIFF_TAMC
481     C-- HPF directive to help TAMC
482     CHPF$ INDEPENDENT
483     #endif /* ALLOW_AUTODIFF_TAMC */
484     DO bj=myByLo(myThid),myByHi(myThid)
485     #ifdef ALLOW_AUTODIFF_TAMC
486     C-- HPF directive to help TAMC
487     CHPF$ INDEPENDENT
488     #endif /* ALLOW_AUTODIFF_TAMC */
489     DO bi=myBxLo(myThid),myBxHi(myThid)
490    
491     #ifdef ALLOW_AUTODIFF_TAMC
492     act1 = bi - myBxLo(myThid)
493     max1 = myBxHi(myThid) - myBxLo(myThid) + 1
494     act2 = bj - myByLo(myThid)
495     max2 = myByHi(myThid) - myByLo(myThid) + 1
496     act3 = myThid - 1
497     max3 = nTx*nTy
498     act4 = ikey_dynamics - 1
499     itdkey = (act1 + 1) + act2*max1
500     & + act3*max1*max2
501     & + act4*max1*max2*max3
502     #endif /* ALLOW_AUTODIFF_TAMC */
503    
504     C-- Set up work arrays with valid (i.e. not NaN) values
505     C These inital values do not alter the numerical results. They
506     C just ensure that all memory references are to valid floating
507     C point numbers. This prevents spurious hardware signals due to
508     C uninitialised but inert locations.
509     DO k=1,Nr
510     DO j=1-OLy,sNy+OLy
511     DO i=1-OLx,sNx+OLx
512     C This is currently used by GMRedi, IVDC, MXL-depth and Diagnostics
513     sigmaX(i,j,k) = 0. _d 0
514     sigmaY(i,j,k) = 0. _d 0
515     sigmaR(i,j,k) = 0. _d 0
516     ENDDO
517     ENDDO
518     ENDDO
519    
520     #ifdef ALLOW_AUTODIFF
521     DO j=1-OLy,sNy+OLy
522     DO i=1-OLx,sNx+OLx
523     rhoKm1 (i,j) = 0. _d 0
524     rhoKp1 (i,j) = 0. _d 0
525     ENDDO
526     ENDDO
527     cph all the following init. are necessary for TAF
528     cph although some of these are re-initialised later.
529     DO k=1,Nr
530     DO j=1-OLy,sNy+OLy
531     DO i=1-OLx,sNx+OLx
532     rhoInSitu(i,j,k,bi,bj) = 0.
533     # ifdef ALLOW_GGL90
534     GGL90viscArU(i,j,k,bi,bj) = 0. _d 0
535     GGL90viscArV(i,j,k,bi,bj) = 0. _d 0
536     GGL90diffKr(i,j,k,bi,bj) = 0. _d 0
537     # endif /* ALLOW_GGL90 */
538     ENDDO
539     ENDDO
540     ENDDO
541     #ifdef ALLOW_OFFLINE
542     IF ( calcConvect ) THEN
543     #endif
544     DO k=1,Nr
545     DO j=1-OLy,sNy+OLy
546     DO i=1-OLx,sNx+OLx
547     IVDConvCount(i,j,k,bi,bj) = 0.
548     ENDDO
549     ENDDO
550     ENDDO
551     #ifdef ALLOW_OFFLINE
552     ENDIF
553     IF ( calcGMRedi ) THEN
554     #endif
555     # ifdef ALLOW_GMREDI
556     DO k=1,Nr
557     DO j=1-OLy,sNy+OLy
558     DO i=1-OLx,sNx+OLx
559     Kwx(i,j,k,bi,bj) = 0. _d 0
560     Kwy(i,j,k,bi,bj) = 0. _d 0
561     Kwz(i,j,k,bi,bj) = 0. _d 0
562     # ifdef GM_NON_UNITY_DIAGONAL
563     Kux(i,j,k,bi,bj) = 0. _d 0
564     Kvy(i,j,k,bi,bj) = 0. _d 0
565     # endif
566     # ifdef GM_EXTRA_DIAGONAL
567     Kuz(i,j,k,bi,bj) = 0. _d 0
568     Kvz(i,j,k,bi,bj) = 0. _d 0
569     # endif
570     # ifdef GM_BOLUS_ADVEC
571     GM_PsiX(i,j,k,bi,bj) = 0. _d 0
572     GM_PsiY(i,j,k,bi,bj) = 0. _d 0
573     # endif
574     # ifdef GM_VISBECK_VARIABLE_K
575     VisbeckK(i,j,bi,bj) = 0. _d 0
576     # endif
577     ENDDO
578     ENDDO
579     ENDDO
580     # endif /* ALLOW_GMREDI */
581     #ifdef ALLOW_OFFLINE
582     ENDIF
583     IF ( calcKPP ) THEN
584     #endif
585     # ifdef ALLOW_KPP
586     DO k=1,Nr
587     DO j=1-OLy,sNy+OLy
588     DO i=1-OLx,sNx+OLx
589     KPPdiffKzS(i,j,k,bi,bj) = 0. _d 0
590     KPPdiffKzT(i,j,k,bi,bj) = 0. _d 0
591     ENDDO
592     ENDDO
593     ENDDO
594     # endif /* ALLOW_KPP */
595     #ifdef ALLOW_OFFLINE
596     ENDIF
597     #endif
598     #endif /* ALLOW_AUTODIFF */
599    
600     #ifdef ALLOW_AUTODIFF_TAMC
601     CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey,
602     CADJ & kind = isbyte
603     CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey,
604     CADJ & kind = isbyte
605     CADJ STORE totphihyd(:,:,:,bi,bj)
606     CADJ & = comlev1_bibj, key=itdkey,
607     CADJ & kind = isbyte
608     # ifdef ALLOW_KPP
609     CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey,
610     CADJ & kind = isbyte
611     CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey,
612     CADJ & kind = isbyte
613     # endif
614     # ifdef ALLOW_SALT_PLUME
615     CADJ STORE saltplumedepth(:,:,bi,bj) = comlev1_bibj, key=itdkey,
616     CADJ & kind = isbyte
617     # endif
618     #endif /* ALLOW_AUTODIFF_TAMC */
619    
620     C-- Always compute density (stored in common block) here; even when it is not
621     C needed here, will be used anyway in calc_phi_hyd (data flow easier this way)
622     #ifdef ALLOW_DEBUG
623     IF (debugMode) CALL DEBUG_CALL('FIND_RHO_2D (xNr)',myThid)
624     #endif
625     #ifdef ALLOW_AUTODIFF
626     IF ( fluidIsWater ) THEN
627     #endif /* ALLOW_AUTODIFF */
628     #ifdef ALLOW_DOWN_SLOPE
629     IF ( useDOWN_SLOPE ) THEN
630     DO k=1,Nr
631     CALL DWNSLP_CALC_RHO(
632     I theta, salt,
633     O rhoInSitu(1-OLx,1-OLy,k,bi,bj),
634     I k, bi, bj, myTime, myIter, myThid )
635     ENDDO
636     ENDIF
637     #endif /* ALLOW_DOWN_SLOPE */
638     #ifdef ALLOW_BBL
639     IF ( useBBL ) THEN
640     C pkg/bbl requires in-situ bbl density for depths equal to and deeper than the bbl.
641     C To reduce computation and storage requirement, these densities are stored in the
642     C dry grid boxes of rhoInSitu. See BBL_CALC_RHO for details.
643     DO k=Nr,1,-1
644     CALL BBL_CALC_RHO(
645     I theta, salt,
646     O rhoInSitu,
647     I k, bi, bj, myTime, myIter, myThid )
648    
649     ENDDO
650     ENDIF
651     #endif /* ALLOW_BBL */
652     IF ( .NOT. ( useDOWN_SLOPE .OR. useBBL ) ) THEN
653     DO k=1,Nr
654     CALL FIND_RHO_2D(
655     I iMin, iMax, jMin, jMax, k,
656     I theta(1-OLx,1-OLy,k,bi,bj),
657     I salt (1-OLx,1-OLy,k,bi,bj),
658     O rhoInSitu(1-OLx,1-OLy,k,bi,bj),
659     I k, bi, bj, myThid )
660     ENDDO
661     ENDIF
662     #ifdef ALLOW_AUTODIFF
663     ELSE
664     C- fluid is not water:
665     DO k=1,Nr
666     DO j=1-OLy,sNy+OLy
667     DO i=1-OLx,sNx+OLx
668     rhoInSitu(i,j,k,bi,bj) = 0.
669     ENDDO
670     ENDDO
671     ENDDO
672     ENDIF
673     #endif /* ALLOW_AUTODIFF */
674    
675     #ifdef ALLOW_DEBUG
676     IF (debugMode) CALL DEBUG_MSG('ENTERING UPWARD K LOOP',myThid)
677     #endif
678    
679     C-- Start of diagnostic loop
680     DO k=Nr,1,-1
681    
682     #ifdef ALLOW_AUTODIFF_TAMC
683     C? Patrick, is this formula correct now that we change the loop range?
684     C? Do we still need this?
685     cph kkey formula corrected.
686     cph Needed for rhoK, rhoKm1, in the case useGMREDI.
687     kkey = (itdkey-1)*Nr + k
688     #endif /* ALLOW_AUTODIFF_TAMC */
689    
690     c#ifdef ALLOW_AUTODIFF_TAMC
691     cCADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey,
692     cCADJ & kind = isbyte
693     cCADJ STORE salt(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey,
694     cCADJ & kind = isbyte
695     c#endif /* ALLOW_AUTODIFF_TAMC */
696    
697     C-- Calculate gradients of potential density for isoneutral
698     C slope terms (e.g. GM/Redi tensor or IVDC diffusivity)
699     IF ( calcGMRedi .OR. (k.GT.1 .AND. calcConvect)
700     & .OR. usePP81 .OR. useMY82 .OR. useGGL90
701     & .OR. useSALT_PLUME .OR. doDiagsRho.GE.1 ) THEN
702     IF (k.GT.1) THEN
703     #ifdef ALLOW_AUTODIFF_TAMC
704     CADJ STORE theta(:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey,
705     CADJ & kind = isbyte
706     CADJ STORE salt (:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey,
707     CADJ & kind = isbyte
708     CADJ STORE rhokm1 (bi,bj) = comlev1_bibj_k, key=kkey,
709     CADJ & kind = isbyte
710     #endif /* ALLOW_AUTODIFF_TAMC */
711     CALL FIND_RHO_2D(
712     I iMin, iMax, jMin, jMax, k,
713     I theta(1-OLx,1-OLy,k-1,bi,bj),
714     I salt (1-OLx,1-OLy,k-1,bi,bj),
715     O rhoKm1,
716     I k-1, bi, bj, myThid )
717     ENDIF
718     #ifdef ALLOW_DEBUG
719     IF (debugMode) CALL DEBUG_CALL('GRAD_SIGMA',myThid)
720     #endif
721     cph Avoid variable aliasing for adjoint !!!
722     DO j=jMin,jMax
723     DO i=iMin,iMax
724     rhoKp1(i,j) = rhoInSitu(i,j,k,bi,bj)
725     ENDDO
726     ENDDO
727     CALL GRAD_SIGMA(
728     I bi, bj, iMin, iMax, jMin, jMax, k,
729     I rhoInSitu(1-OLx,1-OLy,k,bi,bj), rhoKm1, rhoKp1,
730     O sigmaX, sigmaY, sigmaR,
731     I myThid )
732     #ifdef ALLOW_AUTODIFF
733     #ifdef GMREDI_WITH_STABLE_ADJOINT
734     cgf zero out adjoint fields to stabilize pkg/gmredi adjoint
735     cgf -> cuts adjoint dependency from slope to state
736     CALL ZERO_ADJ_LOC( Nr, sigmaX, myThid)
737     CALL ZERO_ADJ_LOC( Nr, sigmaY, myThid)
738     CALL ZERO_ADJ_LOC( Nr, sigmaR, myThid)
739     #endif
740     #endif /* ALLOW_AUTODIFF */
741     ENDIF
742    
743     C-- Implicit Vertical Diffusion for Convection
744     IF (k.GT.1 .AND. calcConvect) THEN
745     #ifdef ALLOW_DEBUG
746     IF (debugMode) CALL DEBUG_CALL('CALC_IVDC',myThid)
747     #endif
748     CALL CALC_IVDC(
749     I bi, bj, iMin, iMax, jMin, jMax, k,
750     I sigmaR,
751     I myTime, myIter, myThid)
752     ENDIF
753    
754     #ifdef ALLOW_DIAGNOSTICS
755     IF ( doDiagsRho.GE.4 ) THEN
756     CALL DIAGS_RHO_L( doDiagsRho, k, bi, bj,
757     I rhoInSitu(1-OLx,1-OLy,1,bi,bj),
758     I rhoKm1, wVel,
759     I myTime, myIter, myThid )
760     ENDIF
761     #endif
762    
763     C-- end of diagnostic k loop (Nr:1)
764     ENDDO
765    
766     #ifdef ALLOW_AUTODIFF_TAMC
767     CADJ STORE IVDConvCount(:,:,:,bi,bj)
768     CADJ & = comlev1_bibj, key=itdkey,
769     CADJ & kind = isbyte
770     #endif
771    
772     C-- Diagnose Mixed Layer Depth:
773     IF ( calcGMRedi .OR. MOD(doDiagsRho,2).EQ.1 ) THEN
774     CALL CALC_OCE_MXLAYER(
775     I rhoInSitu(1-OLx,1-OLy,1,bi,bj), sigmaR,
776     I bi, bj, myTime, myIter, myThid )
777     ENDIF
778    
779     #ifdef ALLOW_SALT_PLUME
780     IF ( useSALT_PLUME ) THEN
781     CALL SALT_PLUME_CALC_DEPTH(
782     I rhoInSitu(1-OLx,1-OLy,1,bi,bj), sigmaR,
783     I bi, bj, myTime, myIter, myThid )
784     #ifdef SALT_PLUME_VOLUME
785     CALL SALT_PLUME_VOLFRAC(
786     I bi, bj, myTime, myIter, myThid )
787     #endif /* SALT_PLUME_VOLUME */
788     ENDIF
789     #endif /* ALLOW_SALT_PLUME */
790    
791     #ifdef ALLOW_DIAGNOSTICS
792     IF ( MOD(doDiagsRho,4).GE.2 ) THEN
793     CALL DIAGNOSTICS_FILL (sigmaR, 'DRHODR ', 0, Nr,
794     & 2, bi, bj, myThid)
795     ENDIF
796     #endif /* ALLOW_DIAGNOSTICS */
797    
798     C-- This is where EXTERNAL_FORCING_SURF(bi,bj) used to be called;
799     C now called earlier, before bi,bj loop.
800    
801     #ifdef ALLOW_AUTODIFF_TAMC
802     cph needed for KPP
803     CADJ STORE surfaceForcingU(:,:,bi,bj)
804     CADJ & = comlev1_bibj, key=itdkey,
805     CADJ & kind = isbyte
806     CADJ STORE surfaceForcingV(:,:,bi,bj)
807     CADJ & = comlev1_bibj, key=itdkey,
808     CADJ & kind = isbyte
809     CADJ STORE surfaceForcingS(:,:,bi,bj)
810     CADJ & = comlev1_bibj, key=itdkey,
811     CADJ & kind = isbyte
812     CADJ STORE surfaceForcingT(:,:,bi,bj)
813     CADJ & = comlev1_bibj, key=itdkey,
814     CADJ & kind = isbyte
815     CADJ STORE surfaceForcingTice(:,:,bi,bj)
816     CADJ & = comlev1_bibj, key=itdkey,
817     CADJ & kind = isbyte
818     #endif /* ALLOW_AUTODIFF_TAMC */
819    
820     #ifdef ALLOW_KPP
821     C-- Compute KPP mixing coefficients
822     IF ( calcKPP ) THEN
823     #ifdef ALLOW_DEBUG
824     IF (debugMode) CALL DEBUG_CALL('KPP_CALC',myThid)
825     #endif
826     CALL TIMER_START('KPP_CALC [DO_OCEANIC_PHYS]', myThid)
827     CALL KPP_CALC(
828     I bi, bj, myTime, myIter, myThid )
829     CALL TIMER_STOP ('KPP_CALC [DO_OCEANIC_PHYS]', myThid)
830     #if (defined ALLOW_AUTODIFF) && !(defined ALLOW_OFFLINE)
831     ELSE
832     CALL KPP_CALC_DUMMY(
833     I bi, bj, myTime, myIter, myThid )
834     #endif /* ALLOW_AUTODIFF and not ALLOW_OFFLINE */
835     ENDIF
836     #endif /* ALLOW_KPP */
837    
838     #ifdef ALLOW_PP81
839     C-- Compute PP81 mixing coefficients
840     IF (usePP81) THEN
841     #ifdef ALLOW_DEBUG
842     IF (debugMode) CALL DEBUG_CALL('PP81_CALC',myThid)
843     #endif
844     CALL PP81_CALC(
845     I bi, bj, sigmaR, myTime, myIter, myThid )
846     ENDIF
847     #endif /* ALLOW_PP81 */
848    
849     #ifdef ALLOW_MY82
850     C-- Compute MY82 mixing coefficients
851     IF (useMY82) THEN
852     #ifdef ALLOW_DEBUG
853     IF (debugMode) CALL DEBUG_CALL('MY82_CALC',myThid)
854     #endif
855     CALL MY82_CALC(
856     I bi, bj, sigmaR, myTime, myIter, myThid )
857     ENDIF
858     #endif /* ALLOW_MY82 */
859    
860     #ifdef ALLOW_GGL90
861     #ifdef ALLOW_AUTODIFF_TAMC
862     CADJ STORE GGL90TKE (:,:,:,bi,bj) = comlev1_bibj, key=itdkey,
863     CADJ & kind = isbyte
864     #endif /* ALLOW_AUTODIFF_TAMC */
865     C-- Compute GGL90 mixing coefficients
866     IF (useGGL90) THEN
867     #ifdef ALLOW_DEBUG
868     IF (debugMode) CALL DEBUG_CALL('GGL90_CALC',myThid)
869     #endif
870     CALL TIMER_START('GGL90_CALC [DO_OCEANIC_PHYS]', myThid)
871     CALL GGL90_CALC(
872     I bi, bj, sigmaR, myTime, myIter, myThid )
873     CALL TIMER_STOP ('GGL90_CALC [DO_OCEANIC_PHYS]', myThid)
874     ENDIF
875     #endif /* ALLOW_GGL90 */
876    
877     #ifdef ALLOW_TIMEAVE
878     IF ( taveFreq.GT. 0. _d 0 ) THEN
879     CALL TIMEAVE_SURF_FLUX( bi, bj, myTime, myIter, myThid)
880     ENDIF
881     IF ( taveFreq.GT.0. .AND. calcConvect ) THEN
882     CALL TIMEAVE_CUMULATE(ConvectCountTave, IVDConvCount,
883     I Nr, deltaTClock, bi, bj, myThid)
884     ENDIF
885     #endif /* ALLOW_TIMEAVE */
886    
887     #ifdef ALLOW_GMREDI
888     #ifdef ALLOW_AUTODIFF_TAMC
889     # ifndef GM_EXCLUDE_CLIPPING
890     cph storing here is needed only for one GMREDI_OPTIONS:
891     cph define GM_BOLUS_ADVEC
892     cph keep it although TAF says you dont need to.
893     cph but I have avoided the #ifdef for now, in case more things change
894     CADJ STORE sigmaX(:,:,:) = comlev1_bibj, key=itdkey,
895     CADJ & kind = isbyte
896     CADJ STORE sigmaY(:,:,:) = comlev1_bibj, key=itdkey,
897     CADJ & kind = isbyte
898     CADJ STORE sigmaR(:,:,:) = comlev1_bibj, key=itdkey,
899     CADJ & kind = isbyte
900     # endif
901     #endif /* ALLOW_AUTODIFF_TAMC */
902    
903     C-- Calculate iso-neutral slopes for the GM/Redi parameterisation
904     IF ( calcGMRedi ) THEN
905     #ifdef ALLOW_DEBUG
906     IF (debugMode) CALL DEBUG_CALL('GMREDI_CALC_TENSOR',myThid)
907     #endif
908     CALL GMREDI_CALC_TENSOR(
909     I iMin, iMax, jMin, jMax,
910     I sigmaX, sigmaY, sigmaR,
911     I bi, bj, myTime, myIter, myThid )
912     #if (defined ALLOW_AUTODIFF) && !(defined ALLOW_OFFLINE)
913     ELSE
914     CALL GMREDI_CALC_TENSOR_DUMMY(
915     I iMin, iMax, jMin, jMax,
916     I sigmaX, sigmaY, sigmaR,
917     I bi, bj, myTime, myIter, myThid )
918     #endif /* ALLOW_AUTODIFF and not ALLOW_OFFLINE */
919     ENDIF
920     #endif /* ALLOW_GMREDI */
921    
922     #ifdef ALLOW_DOWN_SLOPE
923     IF ( useDOWN_SLOPE ) THEN
924     C-- Calculate Downsloping Flow for Down_Slope parameterization
925     IF ( usingPCoords ) THEN
926     CALL DWNSLP_CALC_FLOW(
927     I bi, bj, kSurfC, rhoInSitu,
928     I myTime, myIter, myThid )
929     ELSE
930     CALL DWNSLP_CALC_FLOW(
931     I bi, bj, kLowC, rhoInSitu,
932     I myTime, myIter, myThid )
933     ENDIF
934     ENDIF
935     #endif /* ALLOW_DOWN_SLOPE */
936    
937     C-- end bi,bj loops.
938     ENDDO
939     ENDDO
940    
941     #ifndef ALLOW_AUTODIFF
942     C--- if fluid Is Water: end
943     ENDIF
944     #endif
945    
946     #ifdef ALLOW_BBL
947     IF ( useBBL ) THEN
948     CALL BBL_CALC_RHS(
949     I myTime, myIter, myThid )
950     ENDIF
951     #endif /* ALLOW_BBL */
952    
953     #ifdef ALLOW_MYPACKAGE
954     IF ( useMYPACKAGE ) THEN
955     CALL MYPACKAGE_CALC_RHS(
956     I myTime, myIter, myThid )
957     ENDIF
958     #endif /* ALLOW_MYPACKAGE */
959    
960     #ifdef ALLOW_GMREDI
961     IF ( calcGMRedi ) THEN
962     CALL GMREDI_DO_EXCH( myTime, myIter, myThid )
963     ENDIF
964     #endif /* ALLOW_GMREDI */
965    
966     #ifdef ALLOW_KPP
967     IF ( calcKPP ) THEN
968     CALL KPP_DO_EXCH( myThid )
969     ENDIF
970     #endif /* ALLOW_KPP */
971    
972     #ifdef ALLOW_DIAGNOSTICS
973     IF ( fluidIsWater .AND. useDiagnostics ) THEN
974     CALL DIAGS_RHO_G(
975     I rhoInSitu, uVel, vVel, wVel,
976     I myTime, myIter, myThid )
977     ENDIF
978     IF ( useDiagnostics ) THEN
979     CALL DIAGS_OCEANIC_SURF_FLUX( myTime, myIter, myThid )
980     ENDIF
981     IF ( calcConvect .AND. useDiagnostics ) THEN
982     CALL DIAGNOSTICS_FILL( IVDConvCount, 'CONVADJ ',
983     & 0, Nr, 0, 1, 1, myThid )
984     ENDIF
985     #endif
986    
987     #ifdef ALLOW_ECCO
988     CALL ECCO_PHYS( myThid )
989     #endif
990    
991     #ifdef ALLOW_DEBUG
992     IF (debugMode) CALL DEBUG_LEAVE('DO_OCEANIC_PHYS',myThid)
993     #endif
994    
995     RETURN
996     END
  ViewVC Help
Powered by ViewVC 1.1.22