Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

do_oceanic_phys.F

Last change on this file was 26744, checked in by Mathieu Morlighem, 3 years ago
merged trunk-jpl and trunk for revision 26742
File size: 34.4 KB

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

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

Download in other formats:

Original Format