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