/[MITgcm]/MITgcm_contrib/MPMice/beaufort/code/do_oceanic_phys.F
ViewVC logotype

Annotation of /MITgcm_contrib/MPMice/beaufort/code/do_oceanic_phys.F

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

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