Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

do_oceanic_phys.F@ 22758

Last change on this file since 22758 was 22758, checked in by Mathieu Morlighem, 7 years ago
merged trunk-jpl and trunk for revision 22757
File size: 36.6 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format