Context Navigation

do_oceanic_phys.F@ 23189

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

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

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

Download in other formats:

Original Format