Context Navigation

← Previous Revision
Next Revision →
Normal
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

Rev	Line
[22532]	1	#include "PACKAGES_CONFIG.h"
	2	#include "CPP_OPTIONS.h"
[26744]	3	#ifdef ALLOW_MOM_COMMON
	4	# include "MOM_COMMON_OPTIONS.h"
	5	#endif
[22532]	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 \|
[26744]	78	C \|-- OBCS_ADJUST
	79	C \|
[22532]	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 \|
[26744]	127	C \|-- GGL90_EXCHANGES
	128	C \|
[22532]	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"
[26744]	152	# ifdef ALLOW_AUTODIFF_TAMC
	153	# include "tamc.h"
	154	# endif
[22532]	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
[26744]	206	C rhoKp1,rhoKm1 :: Density at current level, and @ level minus one
[22532]	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
[26744]	212	C kSrf :: surface index
[22532]	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
[26744]	221	INTEGER i, j, k, kSrf
[22532]	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 */
[26744]	232	#ifdef ALLOW_AUTODIFF_TAMC
	233	INTEGER act1, act2, act3, act4
	234	INTEGER max1, max2, max3
	235	INTEGER kkey, itdkey
	236	#endif
[22532]	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
[26744]	244	kSrf = 1
	245	IF ( usingPCoords ) kSrf = Nr
	246
[22532]	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
[26744]	305	adjustColdSST_diag(i,j,bi,bj) = 0. _d 0
	306	# ifdef ALLOW_SALT_PLUME
[22532]	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
[26744]	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
[22532]	341	# ifdef NONLIN_FRSURF
[26744]	342	CADJ STORE hFac_surfC = comlev1, key=ikey_dynamics, kind=isbyte
[22532]	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
[26744]	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
[22532]	363	#if (defined ALLOW_EXF) && (defined ALLOW_ATM_TEMP)
[26744]	364	CADJ STORE evap = comlev1, key=ikey_dynamics, kind=isbyte
[22532]	365	#endif
[26744]	366	# endif /* ALLOW_AUTODIFF_TAMC */
	367	#ifdef ALLOW_AUTODIFF_TAMC
	368	CADJ STORE phiHydLow= comlev1, key=ikey_dynamics, byte=isbyte
[22532]	369	#endif
	370	IF ( useSEAICE ) THEN
	371	# ifdef ALLOW_AUTODIFF_TAMC
[26744]	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
[22532]	378	# endif
	379	# ifdef SEAICE_VARIABLE_SALINITY
[26744]	380	CADJ STORE hsalt = comlev1, key=ikey_dynamics, kind=isbyte
[22532]	381	# endif
	382	# ifdef ATMOSPHERIC_LOADING
[26744]	383	CADJ STORE pload, siceload = comlev1, key=ikey_dynamics, kind=isbyte
[22532]	384	# endif
	385	# ifdef NONLIN_FRSURF
[26744]	386	CADJ STORE recip_hfacc = comlev1, key=ikey_dynamics, kind=isbyte
[22532]	387	# endif
	388	# ifdef ANNUAL_BALANCE
[26744]	389	CADJ STORE balance_itcount = comlev1, key=ikey_dynamics, kind=isbyte
[22532]	390	# endif /* ANNUAL_BALANCE */
	391	# ifdef ALLOW_THSICE
[26744]	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
[22532]	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)
[26744]	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
[22532]	412	# ifdef ALLOW_COST
	413	CALL SEAICE_COST_SENSI ( myTime, myIter, myThid )
	414	# endif
[26744]	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 */
[22532]	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
[26744]	437	CADJ STORE sst, sss = comlev1, key=ikey_dynamics, kind=isbyte
	438	CADJ STORE qsw = comlev1, key=ikey_dynamics, kind=isbyte
[22532]	439	# ifdef ALLOW_SEAICE
[26744]	440	CADJ STORE area = comlev1, key=ikey_dynamics, kind=isbyte
[22532]	441	# endif
	442	#endif
	443
	444	#ifdef ALLOW_CPL_ISSM
[22681]	445	IF ( useCoupler) CALL CPL_ISSM( myTime, myIter, myThid )
[22532]	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
[26744]	454	CADJ STORE salt, theta = comlev1, key=ikey_dynamics, kind=isbyte
	455	CADJ STORE uvel, vvel = comlev1, key=ikey_dynamics, kind=isbyte
[22532]	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:
[26744]	485	CALL SALT_PLUME_DO_EXCH( myTime, myIter, myThid )
[22532]	486	ENDIF
	487	#endif /* ALLOW_SALT_PLUME */
	488
	489	C-- Freeze water at the surface
	490	IF ( allowFreezing ) THEN
	491	#ifdef ALLOW_AUTODIFF_TAMC
[26744]	492	CADJ STORE theta = comlev1, key=ikey_dynamics, kind=isbyte
[22532]	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
[26744]	505	CADJ STORE salt, theta = comlev1, key=ikey_dynamics, kind=isbyte
[22532]	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
[26744]	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
[22532]	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
[26744]	566	#if (defined (ALLOW_SIGMAR_COST_CONTRIBUTION) \|\| defined (ALLOW_LEITH_QG))
	567	sigmaRfield(i,j,k,bi,bj) = 0. _d 0
	568	#endif
[22532]	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
[26744]	579	#ifdef ALLOW_AUTODIFF
[22532]	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
[26744]	658	CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, kind=isbyte
	659	CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, kind=isbyte
[22532]	660	CADJ STORE totphihyd(:,:,:,bi,bj)
[26744]	661	CADJ & = comlev1_bibj, key=itdkey, kind=isbyte
[22532]	662	# ifdef ALLOW_KPP
[26744]	663	CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, kind=isbyte
	664	CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, kind=isbyte
[22532]	665	# endif
	666	# ifdef ALLOW_SALT_PLUME
[26744]	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
[22532]	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
[26744]	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
[22532]	776	#endif /* ALLOW_AUTODIFF_TAMC */
[26744]	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(
[22532]	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 )
[26744]	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
[22532]	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 )
[26744]	811
	812	#if (defined (ALLOW_SIGMAR_COST_CONTRIBUTION) \|\| defined (ALLOW_LEITH_QG))
[22532]	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
[26744]	818	#endif /* ALLOW_SIGMAR_COST_CONTRIBUTION or ALLOW_LEITH_QG */
	819
[22532]	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)
[26744]	856	CADJ & = comlev1_bibj, key=itdkey, kind=isbyte
[22532]	857	#endif
	858
	859	C-- Diagnose Mixed Layer Depth:
	860	IF ( calcGMRedi .OR. MOD(doDiagsRho,2).EQ.1 ) THEN
	861	CALL CALC_OCE_MXLAYER(
[26744]	862	I rhoInSitu(1-OLx,1-OLy,kSrf,bi,bj), sigmaR,
[22532]	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(
[26744]	869	I rhoInSitu(1-OLx,1-OLy,kSrf,bi,bj), sigmaR,
[22532]	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(
[26744]	876	I 1, bi, bj, recip_hFacC(1-OLx,1-OLy,kSrf,bi,bj),
[22532]	877	I theta, 0,
	878	I myTime, myIter, myThid )
	879	CALL SALT_PLUME_APPLY(
[26744]	880	I 2, bi, bj, recip_hFacC(1-OLx,1-OLy,kSrf,bi,bj),
[22532]	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)
[26744]	904	CADJ & = comlev1_bibj, key=itdkey, kind=isbyte
[22532]	905	CADJ STORE surfaceForcingV(:,:,bi,bj)
[26744]	906	CADJ & = comlev1_bibj, key=itdkey, kind=isbyte
[22532]	907	CADJ STORE surfaceForcingS(:,:,bi,bj)
[26744]	908	CADJ & = comlev1_bibj, key=itdkey, kind=isbyte
[22532]	909	CADJ STORE surfaceForcingT(:,:,bi,bj)
[26744]	910	CADJ & = comlev1_bibj, key=itdkey, kind=isbyte
[22532]	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
[26744]	966	CADJ STORE GGL90TKE(:,:,:,bi,bj)
	967	CADJ & = comlev1_bibj, key=itdkey, kind=isbyte
[22532]	968	#endif /* ALLOW_AUTODIFF_TAMC */
	969	C-- Compute GGL90 mixing coefficients
[26744]	970	IF ( useGGL90 .AND. Nr.GT.1 ) THEN
[22532]	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
[26744]	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
[22532]	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
[26744]	1082	#ifdef ALLOW_GGL90
	1083	IF ( useGGL90 )
	1084	& CALL GGL90_EXCHANGES( myThid )
	1085	#endif /* ALLOW_GGL90 */
	1086
[22532]	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