Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

dynspg_ts.F90

Last change on this file was 15592, checked in by techene, 2 years ago
#2605 : remove obsolete lbc_lnk_multi
Property svn:keywords set to `Id`
File size: 82.2 KB

Line
1	MODULE dynspg_ts
2
3	!! Includes ROMS wd scheme with diagnostic outputs ; puu(:,:,:,Kmm) and puu(:,:,:,Krhs) updates are commented out !
4
5	!!======================================================================
6	!! * MODULE dynspg_ts *
7	!! Ocean dynamics: surface pressure gradient trend, split-explicit scheme
8	!!======================================================================
9	!! History : 1.0 ! 2004-12 (L. Bessieres, G. Madec) Original code
10	!! - ! 2005-11 (V. Garnier, G. Madec) optimization
11	!! - ! 2006-08 (S. Masson) distributed restart using iom
12	!! 2.0 ! 2007-07 (D. Storkey) calls to BDY routines
13	!! - ! 2008-01 (R. Benshila) change averaging method
14	!! 3.2 ! 2009-07 (R. Benshila, G. Madec) Complete revisit associated to vvl reactivation
15	!! 3.3 ! 2010-09 (D. Storkey, E. O'Dea) update for BDY for Shelf configurations
16	!! 3.3 ! 2011-03 (R. Benshila, R. Hordoir, P. Oddo) update calculation of ub_b
17	!! 3.5 ! 2013-07 (J. Chanut) Switch to Forward-backward time stepping
18	!! 3.6 ! 2013-11 (A. Coward) Update for z-tilde compatibility
19	!! 3.7 ! 2015-11 (J. Chanut) free surface simplification
20	!! - ! 2016-12 (G. Madec, E. Clementi) update for Stoke-Drift divergence
21	!! 4.0 ! 2017-05 (G. Madec) drag coef. defined at t-point (zdfdrg.F90)
22	!!---------------------------------------------------------------------
23
24	!!----------------------------------------------------------------------
25	!! dyn_spg_ts : compute surface pressure gradient trend using a time-splitting scheme
26	!! dyn_spg_ts_init: initialisation of the time-splitting scheme
27	!! ts_wgt : set time-splitting weights for temporal averaging (or not)
28	!! ts_rst : read/write time-splitting fields in restart file
29	!!----------------------------------------------------------------------
30	USE oce ! ocean dynamics and tracers
31	USE dom_oce ! ocean space and time domain
32	USE sbc_oce ! surface boundary condition: ocean
33	USE isf_oce ! ice shelf variable (fwfisf)
34	USE zdf_oce ! vertical physics: variables
35	USE zdfdrg ! vertical physics: top/bottom drag coef.
36	USE sbcapr ! surface boundary condition: atmospheric pressure
37	USE dynadv , ONLY: ln_dynadv_vec
38	USE dynvor ! vortivity scheme indicators
39	USE phycst ! physical constants
40	USE dynvor ! vorticity term
41	USE wet_dry ! wetting/drying flux limter
42	USE bdy_oce ! open boundary
43	USE bdyvol ! open boundary volume conservation
44	USE bdytides ! open boundary condition data
45	USE bdydyn2d ! open boundary conditions on barotropic variables
46	USE tide_mod !
47	USE sbcwave ! surface wave
48	#if defined key_agrif
49	USE agrif_oce_interp ! agrif
50	USE agrif_oce
51	#endif
52	#if defined key_asminc
53	USE asminc ! Assimilation increment
54	#endif
55	!
56	USE in_out_manager ! I/O manager
57	USE lib_mpp ! distributed memory computing library
58	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
59	USE prtctl ! Print control
60	USE iom ! IOM library
61	USE restart ! only for lrst_oce
62
63	USE iom ! to remove
64
65	IMPLICIT NONE
66	PRIVATE
67
68	PUBLIC dyn_spg_ts ! called by dyn_spg
69	PUBLIC dyn_spg_ts_init ! - - dyn_spg_init
70	PUBLIC dyn_drg_init ! called by rk3ssh
71
72	!! Time filtered arrays at baroclinic time step:
73	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: un_adv , vn_adv !: Advection vel. at "now" barocl. step
74	!
75	INTEGER, SAVE :: icycle ! Number of barotropic sub-steps for each internal step nn_e <= 2.5 nn_e
76	REAL(wp),SAVE :: rDt_e ! Barotropic time step
77	!
78	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: wgtbtp1, wgtbtp2 ! 1st & 2nd weights used in time filtering of barotropic fields
79	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zwz ! ff_f/h at F points
80	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ftnw, ftne ! triad of coriolis parameter
81	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ftsw, ftse ! (only used with een vorticity scheme)
82
83	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: sshe_rhs ! RHS of ssh Eq.
84	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: Ue_rhs, Ve_rhs ! RHS of barotropic velocity Eq.
85	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: CdU_u, CdU_v ! top/bottom stress at u- & v-points
86
87	REAL(wp) :: r1_12 = 1._wp / 12._wp ! local ratios
88	REAL(wp) :: r1_8 = 0.125_wp !
89	REAL(wp) :: r1_4 = 0.25_wp !
90	REAL(wp) :: r1_2 = 0.5_wp !
91
92	!! * Substitutions
93	# include "do_loop_substitute.h90"
94	# include "domzgr_substitute.h90"
95	!!----------------------------------------------------------------------
96	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
97	!! $Id$
98	!! Software governed by the CeCILL license (see ./LICENSE)
99	!!----------------------------------------------------------------------
100	CONTAINS
101
102	INTEGER FUNCTION dyn_spg_ts_alloc()
103	!!----------------------------------------------------------------------
104	!! * routine dyn_spg_ts_alloc *
105	!!----------------------------------------------------------------------
106	INTEGER :: ierr(3)
107	!!----------------------------------------------------------------------
108	ierr(:) = 0
109	!
110	ALLOCATE( wgtbtp1(3nn_e), wgtbtp2(3nn_e), zwz(jpi,jpj), STAT=ierr(1) )
111	IF( ln_dynvor_een .OR. ln_dynvor_eeT ) &
112	& ALLOCATE( ftnw(jpi,jpj) , ftne(jpi,jpj) , ftsw(jpi,jpj) , ftse(jpi,jpj), STAT=ierr(2) )
113	!
114	ALLOCATE( un_adv(jpi,jpj), vn_adv(jpi,jpj) , STAT=ierr(3) )
115	!
116	dyn_spg_ts_alloc = MAXVAL( ierr(:) )
117	!
118	CALL mpp_sum( 'dynspg_ts', dyn_spg_ts_alloc )
119	IF( dyn_spg_ts_alloc /= 0 ) CALL ctl_stop( 'STOP', 'dyn_spg_ts_alloc: failed to allocate arrays' )
120	!
121	END FUNCTION dyn_spg_ts_alloc
122
123
124	SUBROUTINE dyn_spg_ts( kt, Kbb, Kmm, Krhs, puu, pvv, pssh, puu_b, pvv_b, Kaa )
125	!!----------------------------------------------------------------------
126	!!
127	!! ** Purpose : - Compute the now trend due to the explicit time stepping
128	!! of the quasi-linear barotropic system, and add it to the
129	!! general momentum trend.
130	!!
131	!! ** Method : - split-explicit schem (time splitting) :
132	!! Barotropic variables are advanced from internal time steps
133	!! "n" to "n+1" if ln_bt_fw=T
134	!! or from
135	!! "n-1" to "n+1" if ln_bt_fw=F
136	!! thanks to a generalized forward-backward time stepping (see ref. below).
137	!!
138	!! ** Action :
139	!! -Update the filtered free surface at step "n+1" : pssh(:,:,Kaa)
140	!! -Update filtered barotropic velocities at step "n+1" : puu_b(:,:,:,Kaa), vv_b(:,:,:,Kaa)
141	!! -Compute barotropic advective fluxes at step "n" : un_adv, vn_adv
142	!! These are used to advect tracers and are compliant with discrete
143	!! continuity equation taken at the baroclinic time steps. This
144	!! ensures tracers conservation.
145	!! - (puu(:,:,:,Krhs), pvv(:,:,:,Krhs)) momentum trend updated with barotropic component.
146	!!
147	!! References : Shchepetkin and McWilliams, Ocean Modelling, 2005.
148	!!---------------------------------------------------------------------
149	INTEGER , INTENT( in ) :: kt ! ocean time-step index
150	INTEGER , INTENT( in ) :: Kbb, Kmm, Krhs, Kaa ! ocean time level indices
151	REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation
152	REAL(wp), DIMENSION(jpi,jpj ,jpt), INTENT(inout) :: pssh, puu_b, pvv_b ! SSH and barotropic velocities at main time levels
153	!
154	INTEGER :: ji, jj, jk, jn ! dummy loop indices
155	LOGICAL :: ll_fw_start ! =T : forward integration
156	LOGICAL :: ll_init ! =T : special startup of 2d equations
157	INTEGER :: noffset ! local integers : time offset for bdy update
158	REAL(wp) :: z1_hu , z1_hv ! local scalars
159	REAL(wp) :: zzsshu, zzsshv ! - -
160	REAL(wp) :: za0, za1, za2, za3 ! - -
161	REAL(wp) :: zztmp, zldg ! - -
162	REAL(wp) :: zhu_bck, zhv_bck, zhdiv ! - -
163	REAL(wp) :: zun_save, zvn_save ! - -
164	REAL(wp), DIMENSION(jpi,jpj) :: zu_trd, zu_frc, zu_spg, zssh_frc
165	REAL(wp), DIMENSION(jpi,jpj) :: zv_trd, zv_frc, zv_spg
166	REAL(wp), DIMENSION(jpi,jpj) :: zsshu_a, zhup2_e, zhtp2_e
167	REAL(wp), DIMENSION(jpi,jpj) :: zsshv_a, zhvp2_e, zsshp2_e
168	REAL(wp), DIMENSION(jpi,jpj) :: zCdU_u, zCdU_v ! top/bottom stress at u- & v-points
169	REAL(wp), DIMENSION(jpi,jpj) :: zhU, zhV ! fluxes
170	!
171	REAL(wp) :: zwdramp ! local scalar - only used if ln_wd_dl = .True.
172
173	INTEGER :: iwdg, jwdg, kwdg ! short-hand values for the indices of the output point
174
175	REAL(wp) :: zepsilon, zgamma ! - -
176	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zcpx, zcpy ! Wetting/Dying gravity filter coef.
177	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: ztwdmask, zuwdmask, zvwdmask ! ROMS wetting and drying masks at t,u,v points
178	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zuwdav2, zvwdav2 ! averages over the sub-steps of zuwdmask and zvwdmask
179	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z2d ! 2D workspace
180	REAL(wp) :: zt0substep ! Time of day at the beginning of the time substep
181	!!----------------------------------------------------------------------
182	!
183	IF( ln_wd_il ) ALLOCATE( zcpx(jpi,jpj), zcpy(jpi,jpj) )
184	! !* Allocate temporary arrays
185	IF( ln_wd_dl ) ALLOCATE( ztwdmask(jpi,jpj), zuwdmask(jpi,jpj), zvwdmask(jpi,jpj), zuwdav2(jpi,jpj), zvwdav2(jpi,jpj))
186	!
187	zwdramp = r_rn_wdmin1 ! simplest ramp
188	! zwdramp = 1._wp / (rn_wdmin2 - rn_wdmin1) ! more general ramp
189	! ! inverse of baroclinic time step
190	!
191	ll_init = ln_bt_av ! if no time averaging, then no specific restart
192	ll_fw_start = .FALSE.
193	! ! time offset in steps for bdy data update
194	IF( .NOT.ln_bt_fw ) THEN ; noffset = - nn_e
195	ELSE ; noffset = 0
196	ENDIF
197	!
198	IF( kt == nit000 ) THEN !* initialisation
199	!
200	IF(lwp) WRITE(numout,*)
201	IF(lwp) WRITE(numout,*) 'dyn_spg_ts : surface pressure gradient trend'
202	IF(lwp) WRITE(numout,*) '~~~~~~~~~~ free surface with time splitting'
203	IF(lwp) WRITE(numout,*)
204	!
205	IF( l_1st_euler ) ll_init=.TRUE.
206	!
207	IF( ln_bt_fw .OR. l_1st_euler ) THEN
208	ll_fw_start =.TRUE.
209	noffset = 0
210	ELSE
211	ll_fw_start =.FALSE.
212	ENDIF
213	! ! Set averaging weights and cycle length:
214	CALL ts_wgt( ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2 )
215	!
216	ELSEIF( kt == nit000 + 1 ) THEN !* initialisation 2nd time-step
217	!
218	IF( .NOT.ln_bt_fw ) THEN
219	! If we did an Euler timestep on the first timestep we need to reset ll_fw_start
220	! and the averaging weights. We don't have an easy way of telling whether we did
221	! an Euler timestep on the first timestep (because l_1st_euler is reset to .false.
222	! at the end of the first timestep) so just do this in all cases.
223	ll_fw_start = .FALSE.
224	CALL ts_wgt( ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2 )
225	ENDIF
226	!
227	ENDIF
228	!
229	! -----------------------------------------------------------------------------
230	! Phase 1 : Coupling between general trend and barotropic estimates (1st step)
231	! -----------------------------------------------------------------------------
232	#if defined key_RK3
233	! !========================================!
234	! !== Phase 1 for RK3 time integration ==!
235	! !========================================!
236	!
237	! ! Currently, RK3 requires the forward mode
238	IF( kt == nit000 ) THEN
239	IF( .NOT.ln_bt_fw ) CALL ctl_stop( 'dyn_spg_ts: RK3 requires forward mode (ln_bt_fw=T)' )
240	ENDIF
241	!
242	! ! set values computed in RK3_ssh
243	zssh_frc(:,:) = sshe_rhs(:,:)
244	zu_frc(:,:) = Ue_rhs(:,:)
245	zv_frc(:,:) = Ve_rhs(:,:)
246	zCdU_u (:,:) = CdU_u (:,:)
247	zCdU_v (:,:) = CdU_v (:,:)
248
249	!!gm ==>>> !!ts ISSUe her on en discute changer les cas ENS ENE et triad ?
250	IF( kt == nit000 .OR. .NOT. ln_linssh ) CALL dyn_cor_2D_init( Kmm ) ! Set zwz, the barotropic Coriolis force coefficient
251	! ! recompute zwz = f/depth at every time step for (.NOT.ln_linssh) as the water colomn height changes
252	!
253
254	#else
255	! !========================================!
256	! !== Phase 1 for MLF time integration ==!
257	! !========================================!
258	!
259	!
260	!
261	! != zu_frc = 1/H e3*d/dt(Ua) =! (Vertical mean of Ua, the 3D trends)
262	! ! --------------------------- !
263	# if defined key_qco
264	zu_frc(:,:) = SUM( e3u_0(:,:,: ) * puu(:,:,:,Krhs) * umask(:,:,:), DIM=3 ) * r1_hu_0(:,:)
265	zv_frc(:,:) = SUM( e3v_0(:,:,: ) * pvv(:,:,:,Krhs) * vmask(:,:,:), DIM=3 ) * r1_hv_0(:,:)
266	# else
267	zu_frc(:,:) = SUM( e3u(:,:,:,Kmm) * puu(:,:,:,Krhs) * umask(:,:,:), DIM=3 ) * r1_hu(:,:,Kmm)
268	zv_frc(:,:) = SUM( e3v(:,:,:,Kmm) * pvv(:,:,:,Krhs) * vmask(:,:,:), DIM=3 ) * r1_hv(:,:,Kmm)
269	# endif
270	!
271	!
272	! != U(Krhs) => baroclinic trend =! (remove its vertical mean)
273	DO jk = 1, jpkm1 ! ----------------------------- !
274	puu(:,:,jk,Krhs) = ( puu(:,:,jk,Krhs) - zu_frc(:,:) ) * umask(:,:,jk)
275	pvv(:,:,jk,Krhs) = ( pvv(:,:,jk,Krhs) - zv_frc(:,:) ) * vmask(:,:,jk)
276	END DO
277
278	!!gm Question here when removing the Vertically integrated trends, we remove the vertically integrated NL trends on momentum....
279	!!gm Is it correct to do so ? I think so...
280
281	! != remove 2D Coriolis trend =!
282	! ! -------------------------- !
283	!
284	IF( kt == nit000 .OR. .NOT. ln_linssh ) CALL dyn_cor_2D_init( Kmm ) ! Set zwz, the barotropic Coriolis force coefficient
285	! ! recompute zwz = f/depth at every time step for (.NOT.ln_linssh) as the water colomn height changes
286	!
287	zhU(:,:) = puu_b(:,:,Kmm) * hu(:,:,Kmm) * e2u(:,:) ! now fluxes
288	zhV(:,:) = pvv_b(:,:,Kmm) * hv(:,:,Kmm) * e1v(:,:) ! NB: FULL domain : put a value in last row and column
289	!
290	CALL dyn_cor_2D( ht(:,:), hu(:,:,Kmm), hv(:,:,Kmm), puu_b(:,:,Kmm), pvv_b(:,:,Kmm), zhU, zhV, & ! <<== in
291	& zu_trd, zv_trd ) ! ==>> out
292	!
293	DO_2D( 0, 0, 0, 0 ) ! Remove coriolis term (and possibly spg) from barotropic trend
294	zu_frc(ji,jj) = zu_frc(ji,jj) - zu_trd(ji,jj) * ssumask(ji,jj)
295	zv_frc(ji,jj) = zv_frc(ji,jj) - zv_trd(ji,jj) * ssvmask(ji,jj)
296	END_2D
297	!
298	! != Add bottom stress contribution from baroclinic velocities =!
299	! ! ----------------------------------------------------------- !
300	CALL dyn_drg_init( Kbb, Kmm, puu , pvv , puu_b , pvv_b , & ! <<= IN
301	& zu_frc, zv_frc, zCdU_u, zCdU_v ) ! =>> OUT
302	!
303	! != Add atmospheric pressure forcing =!
304	! ! ---------------------------------- !
305	IF( ln_apr_dyn ) THEN
306	IF( ln_bt_fw ) THEN ! FORWARD integration: use kt+1/2 pressure (NOW+1/2)
307	DO_2D( 0, 0, 0, 0 )
308	zu_frc(ji,jj) = zu_frc(ji,jj) + grav * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) ) * r1_e1u(ji,jj)
309	zv_frc(ji,jj) = zv_frc(ji,jj) + grav * ( ssh_ib (ji ,jj+1) - ssh_ib (ji,jj) ) * r1_e2v(ji,jj)
310	END_2D
311	ELSE ! CENTRED integration: use kt-1/2 + kt+1/2 pressure (NOW)
312	zztmp = grav * r1_2
313	DO_2D( 0, 0, 0, 0 )
314	zu_frc(ji,jj) = zu_frc(ji,jj) + zztmp * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) &
315	& + ssh_ibb(ji+1,jj ) - ssh_ibb(ji,jj) ) * r1_e1u(ji,jj)
316	zv_frc(ji,jj) = zv_frc(ji,jj) + zztmp * ( ssh_ib (ji ,jj+1) - ssh_ib (ji,jj) &
317	& + ssh_ibb(ji ,jj+1) - ssh_ibb(ji,jj) ) * r1_e2v(ji,jj)
318	END_2D
319	ENDIF
320	ENDIF
321	!
322	! != Add wind forcing =!
323	! ! ------------------ !
324	IF( ln_bt_fw ) THEN
325	DO_2D( 0, 0, 0, 0 )
326	zu_frc(ji,jj) = zu_frc(ji,jj) + r1_rho0 * utau(ji,jj) * r1_hu(ji,jj,Kmm)
327	zv_frc(ji,jj) = zv_frc(ji,jj) + r1_rho0 * vtau(ji,jj) * r1_hv(ji,jj,Kmm)
328	END_2D
329	ELSE
330	zztmp = r1_rho0 * r1_2
331	DO_2D( 0, 0, 0, 0 )
332	zu_frc(ji,jj) = zu_frc(ji,jj) + zztmp * ( utau_b(ji,jj) + utau(ji,jj) ) * r1_hu(ji,jj,Kmm)
333	zv_frc(ji,jj) = zv_frc(ji,jj) + zztmp * ( vtau_b(ji,jj) + vtau(ji,jj) ) * r1_hv(ji,jj,Kmm)
334	END_2D
335	ENDIF
336	!
337	! !----------------!
338	! !== sssh_frc ==! Right-Hand-Side of the barotropic ssh equation (over the FULL domain)
339	! !----------------!
340	! != Net water flux forcing applied to a water column =!
341	! ! --------------------------------------------------- !
342	IF (ln_bt_fw) THEN ! FORWARD integration: use kt+1/2 fluxes (NOW+1/2)
343	zssh_frc(:,:) = r1_rho0 * ( emp(:,:) - rnf(:,:) - fwfisf_cav(:,:) - fwfisf_par(:,:) )
344	ELSE ! CENTRED integration: use kt-1/2 + kt+1/2 fluxes (NOW)
345	zztmp = r1_rho0 * r1_2
346	zssh_frc(:,:) = zztmp * ( emp(:,:) + emp_b(:,:) &
347	& - rnf(:,:) - rnf_b(:,:) &
348	& - fwfisf_cav(:,:) - fwfisf_cav_b(:,:) &
349	& - fwfisf_par(:,:) - fwfisf_par_b(:,:) )
350	ENDIF
351	! != Add Stokes drift divergence =! (if exist)
352	IF( ln_sdw ) THEN ! ----------------------------- !
353	zssh_frc(:,:) = zssh_frc(:,:) + div_sd(:,:)
354	ENDIF
355	!
356	! ! ice sheet coupling
357	IF ( ln_isf .AND. ln_isfcpl ) THEN
358	!
359	! ice sheet coupling
360	IF( ln_rstart .AND. kt == nit000 ) THEN
361	zssh_frc(:,:) = zssh_frc(:,:) + risfcpl_ssh(:,:)
362	END IF
363	!
364	! conservation option
365	IF( ln_isfcpl_cons ) THEN
366	zssh_frc(:,:) = zssh_frc(:,:) + risfcpl_cons_ssh(:,:)
367	END IF
368	!
369	END IF
370	!
371	# if defined key_asminc
372	! != Add the IAU weighted SSH increment =!
373	! ! ------------------------------------ !
374	IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN
375	zssh_frc(:,:) = zssh_frc(:,:) - ssh_iau(:,:)
376	ENDIF
377	# endif
378
379	! !== END of Phase 1 for MLF time integration ==!
380	#endif
381
382
383	! != Fill boundary data arrays for AGRIF
384	! ! ------------------------------------
385	#if defined key_agrif
386	IF( .NOT.Agrif_Root() ) CALL agrif_dta_ts( kt )
387	#endif
388	!
389	! -----------------------------------------------------------------------
390	! Phase 2 : Integration of the barotropic equations
391	! -----------------------------------------------------------------------
392	!
393	! ! ==================== !
394	! ! Initialisations !
395	! ! ==================== !
396	! Initialize barotropic variables:
397	IF( ll_init )THEN
398	sshbb_e(:,:) = 0._wp
399	ubb_e (:,:) = 0._wp
400	vbb_e (:,:) = 0._wp
401	sshb_e (:,:) = 0._wp
402	ub_e (:,:) = 0._wp
403	vb_e (:,:) = 0._wp
404	ENDIF
405	!
406	IF( ln_linssh ) THEN ! mid-step ocean depth is fixed (hup2_e=hu_n=hu_0)
407	zhup2_e(:,:) = hu_0(:,:)
408	zhvp2_e(:,:) = hv_0(:,:)
409	zhtp2_e(:,:) = ht_0(:,:)
410	ENDIF
411	!
412	IF( ln_bt_fw ) THEN ! FORWARD integration: start from NOW fields
413	sshn_e(:,:) = pssh (:,:,Kmm)
414	un_e (:,:) = puu_b(:,:,Kmm)
415	vn_e (:,:) = pvv_b(:,:,Kmm)
416	!
417	hu_e (:,:) = hu(:,:,Kmm)
418	hv_e (:,:) = hv(:,:,Kmm)
419	hur_e (:,:) = r1_hu(:,:,Kmm)
420	hvr_e (:,:) = r1_hv(:,:,Kmm)
421	ELSE ! CENTRED integration: start from BEFORE fields
422	sshn_e(:,:) = pssh (:,:,Kbb)
423	un_e (:,:) = puu_b(:,:,Kbb)
424	vn_e (:,:) = pvv_b(:,:,Kbb)
425	!
426	hu_e (:,:) = hu(:,:,Kbb)
427	hv_e (:,:) = hv(:,:,Kbb)
428	hur_e (:,:) = r1_hu(:,:,Kbb)
429	hvr_e (:,:) = r1_hv(:,:,Kbb)
430	ENDIF
431	!
432	! Initialize sums:
433	puu_b (:,:,Kaa) = 0._wp ! After barotropic velocities (or transport if flux form)
434	pvv_b (:,:,Kaa) = 0._wp
435	pssh (:,:,Kaa) = 0._wp ! Sum for after averaged sea level
436	un_adv(:,:) = 0._wp ! Sum for now transport issued from ts loop
437	vn_adv(:,:) = 0._wp
438	!
439	IF( ln_wd_dl ) THEN
440	zuwdmask(:,:) = 0._wp ! set to zero for definiteness (not sure this is necessary)
441	zvwdmask(:,:) = 0._wp !
442	zuwdav2 (:,:) = 0._wp
443	zvwdav2 (:,:) = 0._wp
444	END IF
445
446	! ! ==================== !
447	DO jn = 1, icycle ! sub-time-step loop !
448	! ! ==================== !
449	!
450	l_full_nf_update = jn == icycle ! false: disable full North fold update (performances) for jn = 1 to icycle-1
451	!
452	! !== Update the forcing ==! (BDY and tides)
453	!
454	IF( ln_bdy .AND. ln_tide ) CALL bdy_dta_tides( kt, kit=jn, pt_offset= REAL(noffset+1,wp) )
455	! Update tide potential at the beginning of current time substep
456	IF( ln_tide_pot .AND. ln_tide ) THEN
457	zt0substep = REAL(nsec_day, wp) - 0.5_wprn_Dt + (jn + noffset - 1) rn_Dt / REAL(nn_e, wp)
458	CALL upd_tide(zt0substep, Kmm)
459	END IF
460	!
461	! !== extrapolation at mid-step ==! (jn+1/2)
462	!
463	! !* Set extrapolation coefficients for predictor step:
464	IF ((jn<3).AND.ll_init) THEN ! Forward
465	za1 = 1._wp
466	za2 = 0._wp
467	za3 = 0._wp
468	ELSE ! AB3-AM4 Coefficients: bet=0.281105
469	za1 = 1.781105_wp ! za1 = 3/2 + bet
470	za2 = -1.06221_wp ! za2 = -(1/2 + 2*bet)
471	za3 = 0.281105_wp ! za3 = bet
472	ENDIF
473	!
474	! !* Extrapolate barotropic velocities at mid-step (jn+1/2)
475	!-- m+1/2 m m-1 m-2 --!
476	!-- u = (3/2+beta) u -(1/2+2beta) u + beta u --!
477	!-------------------------------------------------------------------------!
478	ua_e(:,:) = za1 * un_e(:,:) + za2 * ub_e(:,:) + za3 * ubb_e(:,:)
479	va_e(:,:) = za1 * vn_e(:,:) + za2 * vb_e(:,:) + za3 * vbb_e(:,:)
480
481	IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only)
482	! ! ------------------
483	! Extrapolate Sea Level at step jit+0.5:
484	!-- m+1/2 m m-1 m-2 --!
485	!-- ssh = (3/2+beta) ssh -(1/2+2beta) ssh + beta ssh --!
486	!--------------------------------------------------------------------------------!
487	zsshp2_e(:,:) = za1 * sshn_e(:,:) + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:)
488
489	! set wetting & drying mask at tracer points for this barotropic mid-step
490	IF( ln_wd_dl ) CALL wad_tmsk( zsshp2_e, ztwdmask )
491	!
492	! ! ocean t-depth at mid-step
493	zhtp2_e(:,:) = ht_0(:,:) + zsshp2_e(:,:)
494	!
495	! ! ocean u- and v-depth at mid-step (separate DO-loops remove the need of a lbc_lnk)
496	#if defined key_qcoTest_FluxForm
497	! ! 'key_qcoTest_FluxForm' : simple ssh average
498	DO_2D( 1, 0, 1, 1 ) ! not jpi-column
499	zhup2_e(ji,jj) = hu_0(ji,jj) + r1_2 * ( zsshp2_e(ji,jj) + zsshp2_e(ji+1,jj ) ) * ssumask(ji,jj)
500	END_2D
501	DO_2D( 1, 1, 1, 0 )
502	zhvp2_e(ji,jj) = hv_0(ji,jj) + r1_2 * ( zsshp2_e(ji,jj) + zsshp2_e(ji ,jj+1) ) * ssvmask(ji,jj)
503	END_2D
504	#else
505	! ! no 'key_qcoTest_FluxForm' : surface weighted ssh average
506	DO_2D( 1, 0, 1, 1 ) ! not jpi-column
507	zhup2_e(ji,jj) = hu_0(ji,jj) + r1_2 * r1_e1e2u(ji,jj) &
508	& * ( e1e2t(ji ,jj) * zsshp2_e(ji ,jj) &
509	& + e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj) ) * ssumask(ji,jj)
510	END_2D
511	DO_2D( 1, 1, 1, 0 ) ! not jpj-row
512	zhvp2_e(ji,jj) = hv_0(ji,jj) + r1_2 * r1_e1e2v(ji,jj) &
513	& * ( e1e2t(ji,jj ) * zsshp2_e(ji,jj ) &
514	& + e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1) ) * ssvmask(ji,jj)
515	END_2D
516	#endif
517	!
518	ENDIF
519	!
520	! !== after SSH ==! (jn+1)
521	!
522	! ! update (ua_e,va_e) to enforce volume conservation at open boundaries
523	! ! values of zhup2_e and zhvp2_e on the halo are not needed in bdy_vol2d
524	IF( ln_bdy .AND. ln_vol ) CALL bdy_vol2d( kt, jn, ua_e, va_e, zhup2_e, zhvp2_e )
525	!
526	! ! resulting flux at mid-step (not over the full domain)
527	DO_2D( 1, 0, 1, 1 ) ! not jpi-column
528	zhU(ji,jj) = e2u(ji,jj) * ua_e(ji,jj) * zhup2_e(ji,jj)
529	END_2D
530	DO_2D( 1, 1, 1, 0 ) ! not jpj-row
531	zhV(ji,jj) = e1v(ji,jj) * va_e(ji,jj) * zhvp2_e(ji,jj)
532	END_2D
533	!
534	#if defined key_agrif
535	! Set fluxes during predictor step to ensure volume conservation
536	IF( .NOT.Agrif_Root() .AND. ln_bt_fw ) CALL agrif_dyn_ts_flux( jn, zhU, zhV )
537	#endif
538	IF( ln_wd_il ) CALL wad_lmt_bt(zhU, zhV, sshn_e, zssh_frc, rDt_e) !!gm wad_lmt_bt use of lbc_lnk on zhU, zhV
539
540	IF( ln_wd_dl ) THEN ! un_e and vn_e are set to zero at faces where
541	! ! the direction of the flow is from dry cells
542	CALL wad_Umsk( ztwdmask, zhU, zhV, un_e, vn_e, zuwdmask, zvwdmask ) ! not jpi colomn for U, not jpj row for V
543	!
544	ENDIF
545	!
546	!
547	! Compute Sea Level at step jit+1
548	!-- m+1 m m+1/2 --!
549	!-- ssh = ssh - delta_t' * [ frc + div( flux ) ] --!
550	!-------------------------------------------------------------------------!
551	DO_2D( 0, 0, 0, 0 )
552	zhdiv = ( zhU(ji,jj) - zhU(ji-1,jj) + zhV(ji,jj) - zhV(ji,jj-1) ) * r1_e1e2t(ji,jj)
553	ssha_e(ji,jj) = ( sshn_e(ji,jj) - rDt_e * ( zssh_frc(ji,jj) + zhdiv ) ) * ssmask(ji,jj)
554	END_2D
555	!
556	CALL lbc_lnk( 'dynspg_ts', ssha_e, 'T', 1._wp, zhU, 'U', -1._wp, zhV, 'V', -1._wp )
557	!
558	! Duplicate sea level across open boundaries (this is only cosmetic if linssh=T)
559	IF( ln_bdy ) CALL bdy_ssh( ssha_e )
560	#if defined key_agrif
561	IF( .NOT.Agrif_Root() ) CALL agrif_ssh_ts( jn )
562	#endif
563	!
564	! ! Sum over sub-time-steps to compute advective velocities
565	za2 = wgtbtp2(jn) ! zhU, zhV hold fluxes extrapolated at jn+0.5
566	un_adv(:,:) = un_adv(:,:) + za2 * zhU(:,:) * r1_e2u(:,:)
567	vn_adv(:,:) = vn_adv(:,:) + za2 * zhV(:,:) * r1_e1v(:,:)
568	! sum over sub-time-steps to decide which baroclinic velocities to set to zero (zuwdav2 is only used when ln_wd_dl_bc=True)
569	IF ( ln_wd_dl_bc ) THEN
570	DO_2D( 1, 0, 1, 1 ) ! not jpi-column
571	zuwdav2(ji,jj) = zuwdav2(ji,jj) + za2 * zuwdmask(ji,jj)
572	END_2D
573	DO_2D( 1, 1, 1, 0 ) ! not jpj-row
574	zvwdav2(ji,jj) = zvwdav2(ji,jj) + za2 * zvwdmask(ji,jj)
575	END_2D
576	END IF
577	!
578	!
579	! Sea Surface Height at u-,v-points (vvl case only)
580	IF( .NOT.ln_linssh ) THEN
581	#if defined key_qcoTest_FluxForm
582	! ! 'key_qcoTest_FluxForm' : simple ssh average
583	DO_2D( 1, 0, 1, 1 )
584	zsshu_a(ji,jj) = r1_2 * ( ssha_e(ji,jj) + ssha_e(ji+1,jj ) ) * ssumask(ji,jj)
585	END_2D
586	DO_2D( 1, 1, 1, 0 )
587	zsshv_a(ji,jj) = r1_2 * ( ssha_e(ji,jj) + ssha_e(ji ,jj+1) ) * ssvmask(ji,jj)
588	END_2D
589	#else
590	DO_2D( 0, 0, 0, 0 )
591	zsshu_a(ji,jj) = r1_2 * r1_e1e2u(ji,jj) * ( e1e2t(ji ,jj ) * ssha_e(ji ,jj ) &
592	& + e1e2t(ji+1,jj ) * ssha_e(ji+1,jj ) ) * ssumask(ji,jj)
593	zsshv_a(ji,jj) = r1_2 * r1_e1e2v(ji,jj) * ( e1e2t(ji ,jj ) * ssha_e(ji ,jj ) &
594	& + e1e2t(ji ,jj+1) * ssha_e(ji ,jj+1) ) * ssvmask(ji,jj)
595	END_2D
596	#endif
597	ENDIF
598	!
599	! Half-step back interpolation of SSH for surface pressure computation at step jit+1/2
600	!-- m+1/2 m+1 m m-1 m-2 --!
601	!-- ssh' = za0 * ssh + za1 * ssh + za2 * ssh + za3 * ssh --!
602	!------------------------------------------------------------------------------------------!
603	CALL ts_bck_interp( jn, ll_init, za0, za1, za2, za3 ) ! coeficients of the interpolation
604	zsshp2_e(:,:) = za0 * ssha_e(:,:) + za1 * sshn_e (:,:) &
605	& + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:)
606	!
607	! ! Surface pressure gradient
608	zldg = ( 1._wp - rn_scal_load ) * grav ! local factor
609	DO_2D( 0, 0, 0, 0 )
610	zu_spg(ji,jj) = - zldg * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj)
611	zv_spg(ji,jj) = - zldg * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj)
612	END_2D
613	IF( ln_wd_il ) THEN ! W/D : gravity filters applied on pressure gradient
614	CALL wad_spg( zsshp2_e, zcpx, zcpy ) ! Calculating W/D gravity filters
615	DO_2D( 0, 0, 0, 0 )
616	zu_spg(ji,jj) = zu_spg(ji,jj) * zcpx(ji,jj)
617	zv_spg(ji,jj) = zv_spg(ji,jj) * zcpy(ji,jj)
618	END_2D
619	ENDIF
620	!
621	! Add Coriolis trend:
622	! zwz array below or triads normally depend on sea level with ln_linssh=F and should be updated
623	! at each time step. We however keep them constant here for optimization.
624	! Recall that zhU and zhV hold fluxes at jn+0.5 (extrapolated not backward interpolated)
625	CALL dyn_cor_2D( zhtp2_e, zhup2_e, zhvp2_e, ua_e, va_e, zhU, zhV, zu_trd, zv_trd )
626	!
627	! Add tidal astronomical forcing if defined
628	IF ( ln_tide .AND. ln_tide_pot ) THEN
629	DO_2D( 0, 0, 0, 0 )
630	zu_trd(ji,jj) = zu_trd(ji,jj) + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) * r1_e1u(ji,jj)
631	zv_trd(ji,jj) = zv_trd(ji,jj) + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) * r1_e2v(ji,jj)
632	END_2D
633	ENDIF
634	!
635	! Add bottom stresses:
636	!jth do implicitly instead
637	IF ( .NOT. ll_wd ) THEN ! Revert to explicit for bit comparison tests in non wad runs
638	DO_2D( 0, 0, 0, 0 )
639	zu_trd(ji,jj) = zu_trd(ji,jj) + zCdU_u(ji,jj) * un_e(ji,jj) * hur_e(ji,jj)
640	zv_trd(ji,jj) = zv_trd(ji,jj) + zCdU_v(ji,jj) * vn_e(ji,jj) * hvr_e(ji,jj)
641	END_2D
642	ENDIF
643	!
644	! Set next velocities:
645	! Compute barotropic speeds at step jit+1 (h : total height of the water colomn)
646	!-- VECTOR FORM
647	!-- m+1 m / m+1/2 \ --!
648	!-- u = u + delta_t' * \ (1-r)g grad_x( ssh') - f * k vect u + frc / --!
649	!-- --!
650	!-- FLUX FORM --!
651	!-- m+1 __1__ / m m / m+1/2 m+1/2 m+1/2 n \ \ --!
652	!-- u = m+1 \| h * u + delta_t' * \ h * (1-r)g grad_x( ssh') - h * f * k vect u + h * frc / \| --!
653	!-- h \ / --!
654	!------------------------------------------------------------------------------------------------------------------------!
655	IF( ln_dynadv_vec .OR. ln_linssh ) THEN !* Vector form
656	DO_2D( 0, 0, 0, 0 )
657	ua_e(ji,jj) = ( un_e(ji,jj) &
658	& + rDt_e * ( zu_spg(ji,jj) &
659	& + zu_trd(ji,jj) &
660	& + zu_frc(ji,jj) ) &
661	& ) * ssumask(ji,jj)
662
663	va_e(ji,jj) = ( vn_e(ji,jj) &
664	& + rDt_e * ( zv_spg(ji,jj) &
665	& + zv_trd(ji,jj) &
666	& + zv_frc(ji,jj) ) &
667	& ) * ssvmask(ji,jj)
668	END_2D
669	!
670	ELSE !* Flux form
671	DO_2D( 0, 0, 0, 0 )
672	! ! hu_e, hv_e hold depth at jn, zhup2_e, zhvp2_e hold extrapolated depth at jn+1/2
673	! ! backward interpolated depth used in spg terms at jn+1/2
674	#if defined key_qcoTest_FluxForm
675	! ! 'key_qcoTest_FluxForm' : simple ssh average
676	zhu_bck = hu_0(ji,jj) + r1_2 * ( zsshp2_e(ji,jj) + zsshp2_e(ji+1,jj ) ) * ssumask(ji,jj)
677	zhv_bck = hv_0(ji,jj) + r1_2 * ( zsshp2_e(ji,jj) + zsshp2_e(ji ,jj+1) ) * ssvmask(ji,jj)
678	#else
679	zhu_bck = hu_0(ji,jj) + r1_2r1_e1e2u(ji,jj) ( e1e2t(ji ,jj) * zsshp2_e(ji ,jj) &
680	& + e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj) ) * ssumask(ji,jj)
681	zhv_bck = hv_0(ji,jj) + r1_2r1_e1e2v(ji,jj) ( e1e2t(ji,jj ) * zsshp2_e(ji,jj ) &
682	& + e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1) ) * ssvmask(ji,jj)
683	#endif
684	! ! inverse depth at jn+1
685	z1_hu = ssumask(ji,jj) / ( hu_0(ji,jj) + zsshu_a(ji,jj) + 1._wp - ssumask(ji,jj) )
686	z1_hv = ssvmask(ji,jj) / ( hv_0(ji,jj) + zsshv_a(ji,jj) + 1._wp - ssvmask(ji,jj) )
687	!
688	ua_e(ji,jj) = ( hu_e (ji,jj) * un_e (ji,jj) &
689	& + rDt_e * ( zhu_bck * zu_spg (ji,jj) & !
690	& + zhup2_e(ji,jj) * zu_trd (ji,jj) & !
691	& + hu(ji,jj,Kmm) * zu_frc (ji,jj) ) ) * z1_hu
692	!
693	va_e(ji,jj) = ( hv_e (ji,jj) * vn_e (ji,jj) &
694	& + rDt_e * ( zhv_bck * zv_spg (ji,jj) & !
695	& + zhvp2_e(ji,jj) * zv_trd (ji,jj) & !
696	& + hv(ji,jj,Kmm) * zv_frc (ji,jj) ) ) * z1_hv
697	END_2D
698	ENDIF
699	!jth implicit bottom friction:
700	IF ( ll_wd ) THEN ! revert to explicit for bit comparison tests in non wad runs
701	DO_2D( 0, 0, 0, 0 )
702	ua_e(ji,jj) = ua_e(ji,jj) / ( 1._wp - rDt_e * zCdU_u(ji,jj) * hur_e(ji,jj) )
703	va_e(ji,jj) = va_e(ji,jj) / ( 1._wp - rDt_e * zCdU_v(ji,jj) * hvr_e(ji,jj) )
704	END_2D
705	ENDIF
706
707	IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only)
708	DO_2D( 0, 0, 0, 0 )
709	hu_e (ji,jj) = hu_0(ji,jj) + zsshu_a(ji,jj)
710	hur_e(ji,jj) = ssumask(ji,jj) / ( hu_e(ji,jj) + 1._wp - ssumask(ji,jj) )
711	hv_e (ji,jj) = hv_0(ji,jj) + zsshv_a(ji,jj)
712	hvr_e(ji,jj) = ssvmask(ji,jj) / ( hv_e(ji,jj) + 1._wp - ssvmask(ji,jj) )
713	END_2D
714	ENDIF
715	!
716	IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only)
717	CALL lbc_lnk( 'dynspg_ts', ua_e , 'U', -1._wp, va_e , 'V', -1._wp &
718	& , hu_e , 'U', 1._wp, hv_e , 'V', 1._wp &
719	& , hur_e, 'U', 1._wp, hvr_e, 'V', 1._wp )
720	ELSE
721	CALL lbc_lnk( 'dynspg_ts', ua_e , 'U', -1._wp, va_e , 'V', -1._wp )
722	ENDIF
723	! ! open boundaries
724	IF( ln_bdy ) CALL bdy_dyn2d( jn, ua_e, va_e, un_e, vn_e, hur_e, hvr_e, ssha_e )
725	#if defined key_agrif
726	IF( .NOT.Agrif_Root() ) CALL agrif_dyn_ts( jn ) ! Agrif
727	#endif
728	! !* Swap
729	! ! ----
730	ubb_e (:,:) = ub_e (:,:)
731	ub_e (:,:) = un_e (:,:)
732	un_e (:,:) = ua_e (:,:)
733	!
734	vbb_e (:,:) = vb_e (:,:)
735	vb_e (:,:) = vn_e (:,:)
736	vn_e (:,:) = va_e (:,:)
737	!
738	sshbb_e(:,:) = sshb_e(:,:)
739	sshb_e (:,:) = sshn_e(:,:)
740	sshn_e (:,:) = ssha_e(:,:)
741	!
742	! !* Sum over whole bt loop (except in weight average)
743	! ! ----------------------
744	IF( ln_bt_av ) THEN
745	za1 = wgtbtp1(jn)
746	IF( ln_dynadv_vec .OR. ln_linssh ) THEN ! Sum velocities
747	puu_b (:,:,Kaa) = puu_b (:,:,Kaa) + za1 * ua_e (:,:)
748	pvv_b (:,:,Kaa) = pvv_b (:,:,Kaa) + za1 * va_e (:,:)
749	ELSE ! Sum transports
750	IF ( .NOT.ln_wd_dl ) THEN
751	puu_b (:,:,Kaa) = puu_b (:,:,Kaa) + za1 * ua_e (:,:) * hu_e (:,:)
752	pvv_b (:,:,Kaa) = pvv_b (:,:,Kaa) + za1 * va_e (:,:) * hv_e (:,:)
753	ELSE
754	puu_b (:,:,Kaa) = puu_b (:,:,Kaa) + za1 * ua_e (:,:) * hu_e (:,:) * zuwdmask(:,:)
755	pvv_b (:,:,Kaa) = pvv_b (:,:,Kaa) + za1 * va_e (:,:) * hv_e (:,:) * zvwdmask(:,:)
756	ENDIF
757	ENDIF
758	! ! Sum sea level
759	pssh(:,:,Kaa) = pssh(:,:,Kaa) + za1 * ssha_e(:,:)
760	ENDIF
761	!
762	! ! ==================== !
763	END DO ! end loop !
764	! ! ==================== !
765
766
767	! -----------------------------------------------------------------------------
768	! Phase 3. update the general trend with the barotropic trend
769	! -----------------------------------------------------------------------------
770	!
771	IF(.NOT.ln_bt_av ) THEN !* Update Kaa barotropic external mode
772	puu_b(:,:,Kaa) = ua_e (:,:)
773	pvv_b(:,:,Kaa) = va_e (:,:)
774	pssh (:,:,Kaa) = ssha_e(:,:)
775	ENDIF
776
777	#if defined key_RK3
778	! !* RK3 case
779	!
780	IF(.NOT.ln_dynadv_vec .AND. ln_bt_av ) THEN ! at this stage, pssh(:,:,:,Krhs) has been corrected: compute new depths at velocity points
781	!
782	# if defined key_qcoTest_FluxForm
783	! ! 'key_qcoTest_FluxForm' : simple ssh average
784	DO_2D( 0, 0, 0, 0 )
785	zzsshu = r1_2 * ( pssh(ji,jj,Kaa) + pssh(ji+1,jj ,Kaa) ) * ssumask(ji,jj)
786	zzsshv = r1_2 * ( pssh(ji,jj,Kaa) + pssh(ji ,jj+1,Kaa) ) * ssvmask(ji,jj)
787	!
788	! ! Save barotropic velocities (not transport)
789	puu_b(ji,jj,Kaa) = puu_b(ji,jj,Kaa) / ( hu_0(ji,jj) + zzsshu + 1._wp - ssumask(ji,jj) )
790	pvv_b(ji,jj,Kaa) = pvv_b(ji,jj,Kaa) / ( hv_0(ji,jj) + zzsshv + 1._wp - ssvmask(ji,jj) )
791	END_2D
792	# else
793	DO_2D( 0, 0, 0, 0 )
794	zzsshu = r1_2 * r1_e1e2u(ji,jj) * ( e1e2t(ji ,jj) * pssh(ji ,jj,Kaa) &
795	& + e1e2t(ji+1,jj) * pssh(ji+1,jj,Kaa) ) * ssumask(ji,jj)
796	zzsshv = r1_2 * r1_e1e2v(ji,jj) * ( e1e2t(ji,jj ) * pssh(ji,jj ,Kaa) &
797	& + e1e2t(ji,jj+1) * pssh(ji,jj+1,Kaa) ) * ssvmask(ji,jj)
798	!
799	! ! Save barotropic velocities (not transport)
800	puu_b(ji,jj,Kaa) = puu_b(ji,jj,Kaa) / ( hu_0(ji,jj) + zzsshu + 1._wp - ssumask(ji,jj) )
801	pvv_b(ji,jj,Kaa) = pvv_b(ji,jj,Kaa) / ( hv_0(ji,jj) + zzsshv + 1._wp - ssvmask(ji,jj) )
802	END_2D
803	# endif
804	!
805	CALL lbc_lnk( 'dynspg_ts', puu_b, 'U', -1._wp, pvv_b, 'V', -1._wp ) ! Boundary conditions
806	!
807	ENDIF
808	!
809	! advective transport from N to N+1 (i.e. Kbb to Kaa)
810	ub2_b(:,:) = un_adv(:,:) ! Save integrated transport for next computation (NOT USED)
811	vb2_b(:,:) = vn_adv(:,:)
812	!
813	IF( iom_use("ubar") ) THEN ! RK3 single first: hu[N+1/2] = 1/2 ( hu[N] + hu[N+1] )
814	ALLOCATE( z2d(jpi,jpj) )
815	z2d(:,:) = 2._wp / ( hu_e(:,:) + hu(:,:,Kbb) + 1._wp - ssumask(:,:) )
816	CALL iom_put( "ubar", un_adv(:,:)*z2d(:,:) ) ! barotropic i-current
817	z2d(:,:) = 2._wp / ( hv_e(:,:) + hv(:,:,Kbb) + 1._wp - ssvmask(:,:) )
818	CALL iom_put( "vbar", vn_adv(:,:)*z2d(:,:) ) ! barotropic i-current
819	DEALLOCATE( z2d )
820	ENDIF
821	!
822	! !== END Phase 3 for RK3 (forward mode) ==!
823
824	#else
825	! !* MLF case
826	!
827	! Set advective velocity correction:
828	IF( ln_bt_fw ) THEN
829	IF( .NOT.( kt == nit000 .AND. l_1st_euler ) ) THEN
830	DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
831	zun_save = un_adv(ji,jj)
832	zvn_save = vn_adv(ji,jj)
833	! ! apply the previously computed correction
834	un_adv(ji,jj) = r1_2 * ( ub2_b(ji,jj) + zun_save - rn_atfp * un_bf(ji,jj) )
835	vn_adv(ji,jj) = r1_2 * ( vb2_b(ji,jj) + zvn_save - rn_atfp * vn_bf(ji,jj) )
836	! ! Update corrective fluxes for next time step
837	un_bf(ji,jj) = rn_atfp * un_bf(ji,jj) + ( zun_save - ub2_b(ji,jj) )
838	vn_bf(ji,jj) = rn_atfp * vn_bf(ji,jj) + ( zvn_save - vb2_b(ji,jj) )
839	! ! Save integrated transport for next computation
840	ub2_b(ji,jj) = zun_save
841	vb2_b(ji,jj) = zvn_save
842	END_2D
843	ELSE
844	un_bf(:,:) = 0._wp ! corrective fluxes for next time step set to zero
845	vn_bf(:,:) = 0._wp
846	ub2_b(:,:) = un_adv(:,:) ! Save integrated transport for next computation
847	vb2_b(:,:) = vn_adv(:,:)
848	ENDIF
849	ENDIF
850	!
851	! Update barotropic trend:
852	IF( ln_dynadv_vec .OR. ln_linssh ) THEN
853	DO jk=1,jpkm1
854	puu(:,:,jk,Krhs) = puu(:,:,jk,Krhs) + ( puu_b(:,:,Kaa) - puu_b(:,:,Kbb) ) * r1_Dt
855	pvv(:,:,jk,Krhs) = pvv(:,:,jk,Krhs) + ( pvv_b(:,:,Kaa) - pvv_b(:,:,Kbb) ) * r1_Dt
856	END DO
857	ELSE
858	IF(.NOT.ln_bt_av ) THEN ! (puu_b,pvv_b)_Kaa is a velocity (hu,hv)_Kaa = (hu_e,hv_e)
859	!
860	DO jk=1,jpkm1
861	puu(:,:,jk,Krhs) = puu(:,:,jk,Krhs) + r1_hu(:,:,Kmm) &
862	& * ( puu_b(:,:,Kaa)hu_e(:,:) - puu_b(:,:,Kbb) hu(:,:,Kbb) ) * r1_Dt
863	pvv(:,:,jk,Krhs) = pvv(:,:,jk,Krhs) + r1_hv(:,:,Kmm) &
864	& * ( pvv_b(:,:,Kaa)hv_e(:,:) - pvv_b(:,:,Kbb) hv(:,:,Kbb) ) * r1_Dt
865	END DO
866	!
867	ELSE ! at this stage, pssh(:,:,:,Krhs) has been corrected: compute new depths at velocity points
868	!
869	# if defined key_qcoTest_FluxForm
870	! ! 'key_qcoTest_FluxForm' : simple ssh average
871	DO_2D( 1, 0, 1, 0 )
872	zsshu_a(ji,jj) = r1_2 * ( pssh(ji,jj,Kaa) + pssh(ji+1,jj ,Kaa) ) * ssumask(ji,jj)
873	zsshv_a(ji,jj) = r1_2 * ( pssh(ji,jj,Kaa) + pssh(ji ,jj+1,Kaa) ) * ssvmask(ji,jj)
874	END_2D
875	# else
876	DO_2D( 1, 0, 1, 0 )
877	zsshu_a(ji,jj) = r1_2 * r1_e1e2u(ji,jj) * ( e1e2t(ji ,jj) * pssh(ji ,jj,Kaa) &
878	& + e1e2t(ji+1,jj) * pssh(ji+1,jj,Kaa) ) * ssumask(ji,jj)
879	zsshv_a(ji,jj) = r1_2 * r1_e1e2v(ji,jj) * ( e1e2t(ji,jj ) * pssh(ji,jj ,Kaa) &
880	& + e1e2t(ji,jj+1) * pssh(ji,jj+1,Kaa) ) * ssvmask(ji,jj)
881	END_2D
882	# endif
883	CALL lbc_lnk( 'dynspg_ts', zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) ! Boundary conditions
884	!
885	DO jk=1,jpkm1
886	puu(:,:,jk,Krhs) = puu(:,:,jk,Krhs) + r1_hu(:,:,Kmm) &
887	& * ( puu_b(:,:,Kaa) - puu_b(:,:,Kbb) * hu(:,:,Kbb) ) * r1_Dt
888	pvv(:,:,jk,Krhs) = pvv(:,:,jk,Krhs) + r1_hv(:,:,Kmm) &
889	& * ( pvv_b(:,:,Kaa) - pvv_b(:,:,Kbb) * hv(:,:,Kbb) ) * r1_Dt
890	END DO
891	! Save barotropic velocities not transport:
892	puu_b(:,:,Kaa) = puu_b(:,:,Kaa) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - ssumask(:,:) )
893	pvv_b(:,:,Kaa) = pvv_b(:,:,Kaa) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - ssvmask(:,:) )
894	ENDIF
895	ENDIF
896
897
898	! Correct velocities so that the barotropic velocity equals (un_adv, vn_adv) (in all cases)
899	DO jk = 1, jpkm1
900	puu(:,:,jk,Kmm) = ( puu(:,:,jk,Kmm) + un_adv(:,:)r1_hu(:,:,Kmm) - puu_b(:,:,Kmm) ) umask(:,:,jk)
901	pvv(:,:,jk,Kmm) = ( pvv(:,:,jk,Kmm) + vn_adv(:,:)r1_hv(:,:,Kmm) - pvv_b(:,:,Kmm) ) vmask(:,:,jk)
902	END DO
903
904	IF( ln_wd_dl .AND. ln_wd_dl_bc) THEN
905	DO jk = 1, jpkm1
906	puu(:,:,jk,Kmm) = ( un_adv(:,:)*r1_hu(:,:,Kmm) &
907	& + zuwdav2(:,:)(puu(:,:,jk,Kmm) - un_adv(:,:)r1_hu(:,:,Kmm)) ) * umask(:,:,jk)
908	pvv(:,:,jk,Kmm) = ( vn_adv(:,:)*r1_hv(:,:,Kmm) &
909	& + zvwdav2(:,:)(pvv(:,:,jk,Kmm) - vn_adv(:,:)r1_hv(:,:,Kmm)) ) * vmask(:,:,jk)
910	END DO
911	ENDIF
912
913	CALL iom_put( "ubar", un_adv(:,:)*r1_hu(:,:,Kmm) ) ! barotropic i-current
914	CALL iom_put( "vbar", vn_adv(:,:)*r1_hv(:,:,Kmm) ) ! barotropic i-current
915
916	! !== END Phase 3 for MLF time integration ==!
917	#endif
918
919	!
920	#if defined key_agrif
921	! Save time integrated fluxes during child grid integration
922	! (used to update coarse grid transports at next time step)
923	!
924	IF( .NOT.Agrif_Root() .AND. ln_bt_fw .AND. ln_agrif_2way ) THEN
925	IF( Agrif_NbStepint() == 0 ) THEN
926	ub2_i_b(:,:) = 0._wp
927	vb2_i_b(:,:) = 0._wp
928	END IF
929	!
930	za1 = 1._wp / REAL(Agrif_rhot(), wp)
931	ub2_i_b(:,:) = ub2_i_b(:,:) + za1 * ub2_b(:,:)
932	vb2_i_b(:,:) = vb2_i_b(:,:) + za1 * vb2_b(:,:)
933	ENDIF
934	#endif
935	! !* write time-spliting arrays in the restart
936	IF( lrst_oce .AND.ln_bt_fw ) CALL ts_rst( kt, 'WRITE' )
937	!
938	IF( ln_wd_il ) DEALLOCATE( zcpx, zcpy )
939	IF( ln_wd_dl ) DEALLOCATE( ztwdmask, zuwdmask, zvwdmask, zuwdav2, zvwdav2 )
940	!
941	CALL iom_put( "baro_u" , puu_b(:,:,Kmm) ) ! Barotropic U Velocity
942	CALL iom_put( "baro_v" , pvv_b(:,:,Kmm) ) ! Barotropic V Velocity
943	!
944	END SUBROUTINE dyn_spg_ts
945
946
947	SUBROUTINE ts_wgt( ll_av, ll_fw, Kpit, zwgt1, zwgt2)
948	!!---------------------------------------------------------------------
949	!! * ROUTINE ts_wgt *
950	!!
951	!! ** Purpose : Set time-splitting weights for temporal averaging (or not)
952	!!----------------------------------------------------------------------
953	LOGICAL, INTENT(in ) :: ll_av ! temporal averaging=.true.
954	LOGICAL, INTENT(in ) :: ll_fw ! forward time splitting =.true.
955	INTEGER, INTENT(inout) :: Kpit ! cycle length
956	!!
957	INTEGER :: jic, jn, ji ! local integers
958	REAL(wp) :: za1, za2 ! loca scalars
959	REAL(wp), DIMENSION(3*nn_e), INTENT(inout) :: zwgt1, zwgt2 ! Primary & Secondary weights
960	!!----------------------------------------------------------------------
961	!
962	zwgt1(:) = 0._wp
963	zwgt2(:) = 0._wp
964	!
965	! !== Set time index when averaged value is requested ==!
966	IF (ll_fw) THEN ; jic = nn_e
967	ELSE ; jic = 2 * nn_e
968	ENDIF
969	!
970	! !== Set primary weights ==!
971	!
972	IF (ll_av) THEN != Define simple boxcar window for primary weights
973	! ! (width = nn_e, centered around jic)
974	SELECT CASE( nn_bt_flt )
975	!
976	CASE( 0 ) ! No averaging
977	zwgt1(jic) = 1._wp
978	Kpit = jic
979	!
980	CASE( 1 ) ! Boxcar, width = nn_e
981	DO jn = 1, 3*nn_e
982	za1 = ABS( REAL( jn-jic, wp) ) / REAL( nn_e, wp )
983	IF( za1 < 0.5_wp ) THEN
984	zwgt1(jn) = 1._wp
985	Kpit = jn
986	ENDIF
987	END DO
988	!
989	CASE( 2 ) ! Boxcar, width = 2 * nn_e
990	DO jn = 1, 3*nn_e
991	za1 = ABS(REAL( jn-jic, wp) ) / REAL( nn_e, wp )
992	IF( za1 < 1._wp ) THEN
993	zwgt1(jn) = 1._wp
994	Kpit = jn
995	ENDIF
996	END DO
997	!
998	CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_bt_flt' )
999	!
1000	END SELECT
1001
1002	ELSE != No time averaging
1003	zwgt1(jic) = 1._wp
1004	Kpit = jic
1005	ENDIF
1006	!
1007	! !== Set secondary weights ==!
1008	DO jn = 1, Kpit
1009	DO ji = jn, Kpit
1010	zwgt2(jn) = zwgt2(jn) + zwgt1(ji)
1011	END DO
1012	END DO
1013	!
1014	! !== Normalize weigths ==!
1015	za1 = 1._wp / SUM( zwgt1(1:Kpit) )
1016	za2 = 1._wp / SUM( zwgt2(1:Kpit) )
1017	DO jn = 1, Kpit
1018	zwgt1(jn) = zwgt1(jn) * za1
1019	zwgt2(jn) = zwgt2(jn) * za2
1020	END DO
1021	!
1022	END SUBROUTINE ts_wgt
1023
1024
1025	SUBROUTINE ts_rst( kt, cdrw )
1026	!!---------------------------------------------------------------------
1027	!! * ROUTINE ts_rst *
1028	!!
1029	!! ** Purpose : Read or write time-splitting arrays in restart file
1030	!!----------------------------------------------------------------------
1031	INTEGER , INTENT(in) :: kt ! ocean time-step
1032	CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag
1033	!!----------------------------------------------------------------------
1034	!
1035	IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialise
1036	! ! ---------------
1037	IF( ln_rstart .AND. ln_bt_fw .AND. .NOT.l_1st_euler ) THEN !* Read the restart file
1038	CALL iom_get( numror, jpdom_auto, 'ub2_b' , ub2_b (:,:), cd_type = 'U', psgn = -1._wp )
1039	CALL iom_get( numror, jpdom_auto, 'vb2_b' , vb2_b (:,:), cd_type = 'V', psgn = -1._wp )
1040	CALL iom_get( numror, jpdom_auto, 'un_bf' , un_bf (:,:), cd_type = 'U', psgn = -1._wp )
1041	CALL iom_get( numror, jpdom_auto, 'vn_bf' , vn_bf (:,:), cd_type = 'V', psgn = -1._wp )
1042	IF( .NOT.ln_bt_av ) THEN
1043	CALL iom_get( numror, jpdom_auto, 'sshbb_e' , sshbb_e(:,:), cd_type = 'T', psgn = 1._wp )
1044	CALL iom_get( numror, jpdom_auto, 'ubb_e' , ubb_e(:,:), cd_type = 'U', psgn = -1._wp )
1045	CALL iom_get( numror, jpdom_auto, 'vbb_e' , vbb_e(:,:), cd_type = 'V', psgn = -1._wp )
1046	CALL iom_get( numror, jpdom_auto, 'sshb_e' , sshb_e(:,:), cd_type = 'T', psgn = 1._wp )
1047	CALL iom_get( numror, jpdom_auto, 'ub_e' , ub_e(:,:), cd_type = 'U', psgn = -1._wp )
1048	CALL iom_get( numror, jpdom_auto, 'vb_e' , vb_e(:,:), cd_type = 'V', psgn = -1._wp )
1049	ENDIF
1050	#if defined key_agrif
1051	! Read time integrated fluxes
1052	IF ( .NOT.Agrif_Root() ) THEN
1053	CALL iom_get( numror, jpdom_auto, 'ub2_i_b' , ub2_i_b(:,:), cd_type = 'U', psgn = -1._wp )
1054	CALL iom_get( numror, jpdom_auto, 'vb2_i_b' , vb2_i_b(:,:), cd_type = 'V', psgn = -1._wp )
1055	ELSE
1056	ub2_i_b(:,:) = 0._wp ; vb2_i_b(:,:) = 0._wp ! used in the 1st update of agrif
1057	ENDIF
1058	#endif
1059	ELSE !* Start from rest
1060	IF(lwp) WRITE(numout,*)
1061	IF(lwp) WRITE(numout,*) ' ==>>> start from rest: set barotropic values to 0'
1062	ub2_b (:,:) = 0._wp ; vb2_b (:,:) = 0._wp ! used in the 1st interpol of agrif
1063	un_adv (:,:) = 0._wp ; vn_adv (:,:) = 0._wp ! used in the 1st interpol of agrif
1064	un_bf (:,:) = 0._wp ; vn_bf (:,:) = 0._wp ! used in the 1st update of agrif
1065	#if defined key_agrif
1066	ub2_i_b(:,:) = 0._wp ; vb2_i_b(:,:) = 0._wp ! used in the 1st update of agrif
1067	#endif
1068	ENDIF
1069	!
1070	ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN ! Create restart file
1071	! ! -------------------
1072	IF(lwp) WRITE(numout,*) '---- ts_rst ----'
1073	CALL iom_rstput( kt, nitrst, numrow, 'ub2_b' , ub2_b (:,:) )
1074	CALL iom_rstput( kt, nitrst, numrow, 'vb2_b' , vb2_b (:,:) )
1075	CALL iom_rstput( kt, nitrst, numrow, 'un_bf' , un_bf (:,:) )
1076	CALL iom_rstput( kt, nitrst, numrow, 'vn_bf' , vn_bf (:,:) )
1077	!
1078	IF (.NOT.ln_bt_av) THEN
1079	CALL iom_rstput( kt, nitrst, numrow, 'sshbb_e' , sshbb_e(:,:) )
1080	CALL iom_rstput( kt, nitrst, numrow, 'ubb_e' , ubb_e(:,:) )
1081	CALL iom_rstput( kt, nitrst, numrow, 'vbb_e' , vbb_e(:,:) )
1082	CALL iom_rstput( kt, nitrst, numrow, 'sshb_e' , sshb_e(:,:) )
1083	CALL iom_rstput( kt, nitrst, numrow, 'ub_e' , ub_e(:,:) )
1084	CALL iom_rstput( kt, nitrst, numrow, 'vb_e' , vb_e(:,:) )
1085	ENDIF
1086	#if defined key_agrif
1087	! Save time integrated fluxes
1088	IF ( .NOT.Agrif_Root() ) THEN
1089	CALL iom_rstput( kt, nitrst, numrow, 'ub2_i_b' , ub2_i_b(:,:) )
1090	CALL iom_rstput( kt, nitrst, numrow, 'vb2_i_b' , vb2_i_b(:,:) )
1091	ENDIF
1092	#endif
1093	ENDIF
1094	!
1095	END SUBROUTINE ts_rst
1096
1097
1098	SUBROUTINE dyn_spg_ts_init
1099	!!---------------------------------------------------------------------
1100	!! * ROUTINE dyn_spg_ts_init *
1101	!!
1102	!! ** Purpose : Set time splitting options
1103	!!----------------------------------------------------------------------
1104	INTEGER :: ji ,jj ! dummy loop indices
1105	REAL(wp) :: zxr2, zyr2, zcmax ! local scalar
1106	REAL(wp), DIMENSION(jpi,jpj) :: zcu
1107	!!----------------------------------------------------------------------
1108	!
1109	! Max courant number for ext. grav. waves
1110	!
1111	DO_2D( 0, 0, 0, 0 )
1112	zxr2 = r1_e1t(ji,jj) * r1_e1t(ji,jj)
1113	zyr2 = r1_e2t(ji,jj) * r1_e2t(ji,jj)
1114	zcu(ji,jj) = SQRT( grav * MAX(ht_0(ji,jj),0._wp) * (zxr2 + zyr2) )
1115	END_2D
1116	!
1117	zcmax = MAXVAL( zcu(Nis0:Nie0,Njs0:Nje0) )
1118	CALL mpp_max( 'dynspg_ts', zcmax )
1119
1120	! Estimate number of iterations to satisfy a max courant number= rn_bt_cmax
1121	IF( ln_bt_auto ) nn_e = CEILING( rn_Dt / rn_bt_cmax * zcmax)
1122
1123	rDt_e = rn_Dt / REAL( nn_e , wp )
1124	zcmax = zcmax * rDt_e
1125	! Print results
1126	IF(lwp) WRITE(numout,*)
1127	IF(lwp) WRITE(numout,*) 'dyn_spg_ts_init : split-explicit free surface'
1128	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~'
1129	IF( ln_bt_auto ) THEN
1130	IF(lwp) WRITE(numout,*) ' ln_ts_auto =.true. Automatically set nn_e '
1131	IF(lwp) WRITE(numout,*) ' Max. courant number allowed: ', rn_bt_cmax
1132	ELSE
1133	IF(lwp) WRITE(numout,*) ' ln_ts_auto=.false.: Use nn_e in namelist nn_e = ', nn_e
1134	ENDIF
1135	!
1136	IF(ln_bt_av) THEN
1137	IF(lwp) WRITE(numout,*) ' ln_bt_av =.true. ==> Time averaging over nn_e time steps is on '
1138	ELSE
1139	IF(lwp) WRITE(numout,*) ' ln_bt_av =.false. => No time averaging of barotropic variables '
1140	ENDIF
1141	!
1142	IF(ln_bt_fw) THEN
1143	IF(lwp) WRITE(numout,*) ' ln_bt_fw=.true. => Forward integration of barotropic variables '
1144	ELSE
1145	IF(lwp) WRITE(numout,*) ' ln_bt_fw =.false.=> Centred integration of barotropic variables '
1146	ENDIF
1147	!
1148	#if defined key_agrif
1149	! Restrict the use of Agrif to the forward case only
1150	IF( .NOT.ln_bt_fw .AND. .NOT.Agrif_Root() ) CALL ctl_stop( 'AGRIF not implemented if ln_bt_fw=.FALSE.' )
1151	#endif
1152	!
1153	IF(lwp) WRITE(numout,*) ' Time filter choice, nn_bt_flt: ', nn_bt_flt
1154	SELECT CASE ( nn_bt_flt )
1155	CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' Dirac'
1156	CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' Boxcar: width = nn_e'
1157	CASE( 2 ) ; IF(lwp) WRITE(numout,) ' Boxcar: width = 2nn_e'
1158	CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_bt_flt: should 0,1, or 2' )
1159	END SELECT
1160	!
1161	IF(lwp) WRITE(numout,*) ' '
1162	IF(lwp) WRITE(numout,*) ' nn_e = ', nn_e
1163	IF(lwp) WRITE(numout,*) ' Barotropic time step [s] is :', rDt_e
1164	IF(lwp) WRITE(numout,*) ' Maximum Courant number is :', zcmax
1165	!
1166	IF(lwp) WRITE(numout,*) ' Time diffusion parameter rn_bt_alpha: ', rn_bt_alpha
1167	IF ((ln_bt_av.AND.nn_bt_flt/=0).AND.(rn_bt_alpha>0._wp)) THEN
1168	CALL ctl_stop( 'dynspg_ts ERROR: if rn_bt_alpha > 0, remove temporal averaging' )
1169	ENDIF
1170	!
1171	IF( .NOT.ln_bt_av .AND. .NOT.ln_bt_fw ) THEN
1172	CALL ctl_stop( 'dynspg_ts ERROR: No time averaging => only forward integration is possible' )
1173	ENDIF
1174	IF( zcmax>0.9_wp ) THEN
1175	CALL ctl_stop( 'dynspg_ts ERROR: Maximum Courant number is greater than 0.9: Inc. nn_e !' )
1176	ENDIF
1177	!
1178	! ! Allocate time-splitting arrays
1179	IF( dyn_spg_ts_alloc() /= 0 ) CALL ctl_stop('STOP', 'dyn_spg_init: failed to allocate dynspg_ts arrays' )
1180	!
1181	! ! read restart when needed
1182	CALL ts_rst( nit000, 'READ' )
1183	!
1184	END SUBROUTINE dyn_spg_ts_init
1185
1186
1187	SUBROUTINE dyn_cor_2D_init( Kmm )
1188	!!---------------------------------------------------------------------
1189	!! * ROUTINE dyn_cor_2D_init *
1190	!!
1191	!! ** Purpose : Set time splitting options
1192	!! Set arrays to remove/compute coriolis trend.
1193	!! Do it once during initialization if volume is fixed, else at each long time step.
1194	!! Note that these arrays are also used during barotropic loop. These are however frozen
1195	!! although they should be updated in the variable volume case. Not a big approximation.
1196	!! To remove this approximation, copy lines below inside barotropic loop
1197	!! and update depths at T- points (ht) at each barotropic time step
1198	!!
1199	!! Compute zwz = f/h (potential planetary voricity)
1200	!! Compute ftne, ftnw, ftse, ftsw (triad of potential planetary voricity)
1201	!!----------------------------------------------------------------------
1202	INTEGER, INTENT(in) :: Kmm ! Time index
1203	INTEGER :: ji ,jj, jk ! dummy loop indices
1204	REAL(wp) :: z1_ht
1205	!!----------------------------------------------------------------------
1206	!
1207	SELECT CASE( nvor_scheme )
1208	CASE( np_EEN, np_ENE, np_ENS , np_MIX ) != schemes using the same e3f definition
1209	SELECT CASE( nn_e3f_typ ) !* ff_f/e3 at F-point
1210	CASE ( 0 ) ! original formulation (masked averaging of e3t divided by 4)
1211	DO_2D( 0, 0, 0, 0 )
1212	zwz(ji,jj) = ( ht(ji,jj+1) + ht(ji+1,jj+1) &
1213	& + ht(ji,jj ) + ht(ji+1,jj ) ) * 0.25_wp
1214	IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj)
1215	END_2D
1216	CASE ( 1 ) ! new formulation (masked averaging of e3t divided by the sum of mask)
1217	DO_2D( 0, 0, 0, 0 )
1218	zwz(ji,jj) = ( ht(ji,jj+1) + ht(ji+1,jj+1) &
1219	& + ht(ji,jj ) + ht(ji+1,jj ) ) &
1220	& / ( MAX(ssmask(ji,jj+1) + ssmask(ji+1,jj+1) &
1221	& + ssmask(ji,jj ) + ssmask(ji+1,jj ) , 1._wp ) )
1222	IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj)
1223	END_2D
1224	END SELECT
1225	CALL lbc_lnk( 'dynspg_ts', zwz, 'F', 1._wp )
1226	END SELECT
1227	!
1228	SELECT CASE( nvor_scheme )
1229	CASE( np_EEN )
1230	!
1231	ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp
1232	DO_2D( 0, 1, 0, 1 )
1233	ftne(ji,jj) = zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1)
1234	ftnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj ) + zwz(ji ,jj )
1235	ftse(ji,jj) = zwz(ji ,jj ) + zwz(ji ,jj-1) + zwz(ji-1,jj-1)
1236	ftsw(ji,jj) = zwz(ji ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj )
1237	END_2D
1238	!
1239	CASE( np_EET ) != EEN scheme using e3t energy conserving scheme
1240	ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp
1241	DO_2D( 0, 1, 0, 1 )
1242	z1_ht = ssmask(ji,jj) / ( ht(ji,jj) + 1._wp - ssmask(ji,jj) )
1243	ftne(ji,jj) = ( ff_f(ji-1,jj ) + ff_f(ji ,jj ) + ff_f(ji ,jj-1) ) * z1_ht
1244	ftnw(ji,jj) = ( ff_f(ji-1,jj-1) + ff_f(ji-1,jj ) + ff_f(ji ,jj ) ) * z1_ht
1245	ftse(ji,jj) = ( ff_f(ji ,jj ) + ff_f(ji ,jj-1) + ff_f(ji-1,jj-1) ) * z1_ht
1246	ftsw(ji,jj) = ( ff_f(ji ,jj-1) + ff_f(ji-1,jj-1) + ff_f(ji-1,jj ) ) * z1_ht
1247	END_2D
1248	!
1249	END SELECT
1250	!
1251	END SUBROUTINE dyn_cor_2D_init
1252
1253
1254	SUBROUTINE dyn_cor_2D( pht, phu, phv, punb, pvnb, zhU, zhV, zu_trd, zv_trd )
1255	!!---------------------------------------------------------------------
1256	!! * ROUTINE dyn_cor_2D *
1257	!!
1258	!! ** Purpose : Compute u and v coriolis trends
1259	!!----------------------------------------------------------------------
1260	INTEGER :: ji ,jj ! dummy loop indices
1261	REAL(wp) :: zx1, zx2, zy1, zy2, z1_hu, z1_hv ! - -
1262	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pht, phu, phv, punb, pvnb, zhU, zhV
1263	REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: zu_trd, zv_trd
1264	!!----------------------------------------------------------------------
1265	SELECT CASE( nvor_scheme )
1266	CASE( np_ENT ) ! enstrophy conserving scheme (f-point)
1267	DO_2D( 0, 0, 0, 0 )
1268	z1_hu = ssumask(ji,jj) / ( phu(ji,jj) + 1._wp - ssumask(ji,jj) )
1269	z1_hv = ssvmask(ji,jj) / ( phv(ji,jj) + 1._wp - ssvmask(ji,jj) )
1270	zu_trd(ji,jj) = + r1_4 * r1_e1e2u(ji,jj) * z1_hu &
1271	& * ( e1e2t(ji+1,jj)pht(ji+1,jj)ff_t(ji+1,jj) * ( pvnb(ji+1,jj) + pvnb(ji+1,jj-1) ) &
1272	& + e1e2t(ji ,jj)pht(ji ,jj)ff_t(ji ,jj) * ( pvnb(ji ,jj) + pvnb(ji ,jj-1) ) )
1273	!
1274	zv_trd(ji,jj) = - r1_4 * r1_e1e2v(ji,jj) * z1_hv &
1275	& * ( e1e2t(ji,jj+1)pht(ji,jj+1)ff_t(ji,jj+1) * ( punb(ji,jj+1) + punb(ji-1,jj+1) ) &
1276	& + e1e2t(ji,jj )pht(ji,jj )ff_t(ji,jj ) * ( punb(ji,jj ) + punb(ji-1,jj ) ) )
1277	END_2D
1278	!
1279	CASE( np_ENE , np_MIX ) ! energy conserving scheme (t-point) ENE or MIX
1280	DO_2D( 0, 0, 0, 0 )
1281	zy1 = ( zhV(ji,jj-1) + zhV(ji+1,jj-1) ) * r1_e1u(ji,jj)
1282	zy2 = ( zhV(ji,jj ) + zhV(ji+1,jj ) ) * r1_e1u(ji,jj)
1283	zx1 = ( zhU(ji-1,jj) + zhU(ji-1,jj+1) ) * r1_e2v(ji,jj)
1284	zx2 = ( zhU(ji ,jj) + zhU(ji ,jj+1) ) * r1_e2v(ji,jj)
1285	! energy conserving formulation for planetary vorticity term
1286	zu_trd(ji,jj) = r1_4 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
1287	zv_trd(ji,jj) = - r1_4 * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 )
1288	END_2D
1289	!
1290	CASE( np_ENS ) ! enstrophy conserving scheme (f-point)
1291	DO_2D( 0, 0, 0, 0 )
1292	zy1 = r1_8 * ( zhV(ji ,jj-1) + zhV(ji+1,jj-1) &
1293	& + zhV(ji ,jj ) + zhV(ji+1,jj ) ) * r1_e1u(ji,jj)
1294	zx1 = - r1_8 * ( zhU(ji-1,jj ) + zhU(ji-1,jj+1) &
1295	& + zhU(ji ,jj ) + zhU(ji ,jj+1) ) * r1_e2v(ji,jj)
1296	zu_trd(ji,jj) = zy1 * ( zwz(ji ,jj-1) + zwz(ji,jj) )
1297	zv_trd(ji,jj) = zx1 * ( zwz(ji-1,jj ) + zwz(ji,jj) )
1298	END_2D
1299	!
1300	CASE( np_EET , np_EEN ) ! energy & enstrophy scheme (using e3t or e3f)
1301	DO_2D( 0, 0, 0, 0 )
1302	zu_trd(ji,jj) = + r1_12 * r1_e1u(ji,jj) * ( ftne(ji,jj ) * zhV(ji ,jj ) &
1303	& + ftnw(ji+1,jj) * zhV(ji+1,jj ) &
1304	& + ftse(ji,jj ) * zhV(ji ,jj-1) &
1305	& + ftsw(ji+1,jj) * zhV(ji+1,jj-1) )
1306	zv_trd(ji,jj) = - r1_12 * r1_e2v(ji,jj) * ( ftsw(ji,jj+1) * zhU(ji-1,jj+1) &
1307	& + ftse(ji,jj+1) * zhU(ji ,jj+1) &
1308	& + ftnw(ji,jj ) * zhU(ji-1,jj ) &
1309	& + ftne(ji,jj ) * zhU(ji ,jj ) )
1310	END_2D
1311	!
1312	END SELECT
1313	!
1314	END SUBROUTINE dyn_cor_2D
1315
1316
1317	SUBROUTINE wad_tmsk( pssh, ptmsk )
1318	!!----------------------------------------------------------------------
1319	!! * ROUTINE wad_lmt *
1320	!!
1321	!! ** Purpose : set wetting & drying mask at tracer points
1322	!! for the current barotropic sub-step
1323	!!
1324	!! ** Method : ???
1325	!!
1326	!! ** Action : ptmsk : wetting & drying t-mask
1327	!!----------------------------------------------------------------------
1328	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pssh !
1329	REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: ptmsk !
1330	!
1331	INTEGER :: ji, jj ! dummy loop indices
1332	!!----------------------------------------------------------------------
1333	!
1334	IF( ln_wd_dl_rmp ) THEN
1335	DO_2D( 1, 1, 1, 1 )
1336	IF ( pssh(ji,jj) + ht_0(ji,jj) > 2._wp * rn_wdmin1 ) THEN
1337	! IF ( pssh(ji,jj) + ht_0(ji,jj) > rn_wdmin2 ) THEN
1338	ptmsk(ji,jj) = 1._wp
1339	ELSEIF( pssh(ji,jj) + ht_0(ji,jj) > rn_wdmin1 ) THEN
1340	ptmsk(ji,jj) = TANH( 50._wp( ( pssh(ji,jj) + ht_0(ji,jj) - rn_wdmin1 )r_rn_wdmin1) )
1341	ELSE
1342	ptmsk(ji,jj) = 0._wp
1343	ENDIF
1344	END_2D
1345	ELSE
1346	DO_2D( 1, 1, 1, 1 )
1347	IF ( pssh(ji,jj) + ht_0(ji,jj) > rn_wdmin1 ) THEN ; ptmsk(ji,jj) = 1._wp
1348	ELSE ; ptmsk(ji,jj) = 0._wp
1349	ENDIF
1350	END_2D
1351	ENDIF
1352	!
1353	END SUBROUTINE wad_tmsk
1354
1355
1356	SUBROUTINE wad_Umsk( pTmsk, phU, phV, pu, pv, pUmsk, pVmsk )
1357	!!----------------------------------------------------------------------
1358	!! * ROUTINE wad_lmt *
1359	!!
1360	!! ** Purpose : set wetting & drying mask at tracer points
1361	!! for the current barotropic sub-step
1362	!!
1363	!! ** Method : ???
1364	!!
1365	!! ** Action : ptmsk : wetting & drying t-mask
1366	!!----------------------------------------------------------------------
1367	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pTmsk ! W & D t-mask
1368	REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: phU, phV, pu, pv ! ocean velocities and transports
1369	REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pUmsk, pVmsk ! W & D u- and v-mask
1370	!
1371	INTEGER :: ji, jj ! dummy loop indices
1372	!!----------------------------------------------------------------------
1373	!
1374	DO_2D( 1, 0, 1, 1 ) ! not jpi-column
1375	IF ( phU(ji,jj) > 0._wp ) THEN ; pUmsk(ji,jj) = pTmsk(ji ,jj)
1376	ELSE ; pUmsk(ji,jj) = pTmsk(ji+1,jj)
1377	ENDIF
1378	phU(ji,jj) = pUmsk(ji,jj)*phU(ji,jj)
1379	pu (ji,jj) = pUmsk(ji,jj)*pu (ji,jj)
1380	END_2D
1381	!
1382	DO_2D( 1, 1, 1, 0 ) ! not jpj-row
1383	IF ( phV(ji,jj) > 0._wp ) THEN ; pVmsk(ji,jj) = pTmsk(ji,jj )
1384	ELSE ; pVmsk(ji,jj) = pTmsk(ji,jj+1)
1385	ENDIF
1386	phV(ji,jj) = pVmsk(ji,jj)*phV(ji,jj)
1387	pv (ji,jj) = pVmsk(ji,jj)*pv (ji,jj)
1388	END_2D
1389	!
1390	END SUBROUTINE wad_Umsk
1391
1392
1393	SUBROUTINE wad_spg( pshn, zcpx, zcpy )
1394	!!---------------------------------------------------------------------
1395	!! * ROUTINE wad_sp *
1396	!!
1397	!! ** Purpose :
1398	!!----------------------------------------------------------------------
1399	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pshn
1400	REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: zcpx, zcpy
1401	!
1402	INTEGER :: ji ,jj ! dummy loop indices
1403	LOGICAL :: ll_tmp1, ll_tmp2
1404	!!----------------------------------------------------------------------
1405	!
1406	DO_2D( 0, 0, 0, 0 )
1407	ll_tmp1 = MIN( pshn(ji,jj) , pshn(ji+1,jj) ) > &
1408	& MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) .AND. &
1409	& MAX( pshn(ji,jj) + ht_0(ji,jj) , pshn(ji+1,jj) + ht_0(ji+1,jj) ) &
1410	& > rn_wdmin1 + rn_wdmin2
1411	ll_tmp2 = ( ABS( pshn(ji+1,jj) - pshn(ji ,jj)) > 1.E-12 ).AND.( &
1412	& MAX( pshn(ji,jj) , pshn(ji+1,jj) ) > &
1413	& MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 )
1414	IF(ll_tmp1) THEN
1415	zcpx(ji,jj) = 1.0_wp
1416	ELSEIF(ll_tmp2) THEN
1417	! no worries about pshn(ji+1,jj) - pshn(ji ,jj) = 0, it won't happen ! here
1418	zcpx(ji,jj) = ABS( (pshn(ji+1,jj) + ht_0(ji+1,jj) - pshn(ji,jj) - ht_0(ji,jj)) &
1419	& / (pshn(ji+1,jj) - pshn(ji ,jj)) )
1420	zcpx(ji,jj) = max(min( zcpx(ji,jj) , 1.0_wp),0.0_wp)
1421	ELSE
1422	zcpx(ji,jj) = 0._wp
1423	ENDIF
1424	!
1425	ll_tmp1 = MIN( pshn(ji,jj) , pshn(ji,jj+1) ) > &
1426	& MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) .AND. &
1427	& MAX( pshn(ji,jj) + ht_0(ji,jj) , pshn(ji,jj+1) + ht_0(ji,jj+1) ) &
1428	& > rn_wdmin1 + rn_wdmin2
1429	ll_tmp2 = ( ABS( pshn(ji,jj) - pshn(ji,jj+1)) > 1.E-12 ).AND.( &
1430	& MAX( pshn(ji,jj) , pshn(ji,jj+1) ) > &
1431	& MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 )
1432
1433	IF(ll_tmp1) THEN
1434	zcpy(ji,jj) = 1.0_wp
1435	ELSE IF(ll_tmp2) THEN
1436	! no worries about pshn(ji,jj+1) - pshn(ji,jj ) = 0, it won't happen ! here
1437	zcpy(ji,jj) = ABS( (pshn(ji,jj+1) + ht_0(ji,jj+1) - pshn(ji,jj) - ht_0(ji,jj)) &
1438	& / (pshn(ji,jj+1) - pshn(ji,jj )) )
1439	zcpy(ji,jj) = MAX( 0._wp , MIN( zcpy(ji,jj) , 1.0_wp ) )
1440	ELSE
1441	zcpy(ji,jj) = 0._wp
1442	ENDIF
1443	END_2D
1444	!
1445	END SUBROUTINE wad_spg
1446
1447
1448	SUBROUTINE dyn_drg_init( Kbb, Kmm, puu, pvv, puu_b ,pvv_b, pu_RHSi, pv_RHSi, pCdU_u, pCdU_v )
1449	!!----------------------------------------------------------------------
1450	!! * ROUTINE dyn_drg_init *
1451	!!
1452	!! ** Purpose : - add the baroclinic top/bottom drag contribution to
1453	!! the baroclinic part of the barotropic RHS
1454	!! - compute the barotropic drag coefficients
1455	!!
1456	!! ** Method : computation done over the INNER domain only
1457	!!----------------------------------------------------------------------
1458	INTEGER , INTENT(in ) :: Kbb, Kmm ! ocean time level indices
1459	REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(in ) :: puu, pvv ! ocean velocities and RHS of momentum equation
1460	REAL(wp), DIMENSION(jpi,jpj,jpt) , INTENT(in ) :: puu_b, pvv_b ! barotropic velocities at main time levels
1461	REAL(wp), DIMENSION(jpi,jpj) , INTENT(inout) :: pu_RHSi, pv_RHSi ! baroclinic part of the barotropic RHS
1462	REAL(wp), DIMENSION(jpi,jpj) , INTENT( out) :: pCdU_u , pCdU_v ! barotropic drag coefficients
1463	!
1464	INTEGER :: ji, jj ! dummy loop indices
1465	INTEGER :: ikbu, ikbv, iktu, iktv
1466	REAL(wp) :: zztmp
1467	REAL(wp), DIMENSION(jpi,jpj) :: zu_i, zv_i
1468	!!----------------------------------------------------------------------
1469	!
1470	! !== Set the barotropic drag coef. ==!
1471	!
1472	IF( ln_isfcav.OR.ln_drgice_imp ) THEN ! top+bottom friction (ocean cavities)
1473
1474	DO_2D( 0, 0, 0, 0 )
1475	pCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) + rCdU_top(ji+1,jj)+rCdU_top(ji,jj) )
1476	pCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) + rCdU_top(ji,jj+1)+rCdU_top(ji,jj) )
1477	END_2D
1478	ELSE ! bottom friction only
1479	DO_2D( 0, 0, 0, 0 )
1480	pCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) )
1481	pCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) )
1482	END_2D
1483	ENDIF
1484	!
1485	! !== BOTTOM stress contribution from baroclinic velocities ==!
1486	!
1487	IF( ln_bt_fw ) THEN ! FORWARD integration: use NOW bottom baroclinic velocities
1488	DO_2D( 0, 0, 0, 0 )
1489	ikbu = mbku(ji,jj)
1490	ikbv = mbkv(ji,jj)
1491	zu_i(ji,jj) = puu(ji,jj,ikbu,Kmm) - puu_b(ji,jj,Kmm)
1492	zv_i(ji,jj) = pvv(ji,jj,ikbv,Kmm) - pvv_b(ji,jj,Kmm)
1493	END_2D
1494	ELSE ! CENTRED integration: use BEFORE bottom baroclinic velocities
1495	DO_2D( 0, 0, 0, 0 )
1496	ikbu = mbku(ji,jj)
1497	ikbv = mbkv(ji,jj)
1498	zu_i(ji,jj) = puu(ji,jj,ikbu,Kbb) - puu_b(ji,jj,Kbb)
1499	zv_i(ji,jj) = pvv(ji,jj,ikbv,Kbb) - pvv_b(ji,jj,Kbb)
1500	END_2D
1501	ENDIF
1502	!
1503	IF( ln_wd_il ) THEN ! W/D : use the "clipped" bottom friction !!gm explain WHY, please !
1504	zztmp = -1._wp / rDt_e
1505	DO_2D( 0, 0, 0, 0 )
1506	pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + zu_i(ji,jj) * wdrampu(ji,jj) * MAX( &
1507	& r1_hu(ji,jj,Kmm) * r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) , zztmp )
1508	pv_RHSi(ji,jj) = pv_RHSi(ji,jj) + zv_i(ji,jj) * wdrampv(ji,jj) * MAX( &
1509	& r1_hv(ji,jj,Kmm) * r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) , zztmp )
1510	END_2D
1511	ELSE ! use "unclipped" drag (even if explicit friction is used in 3D calculation)
1512	DO_2D( 0, 0, 0, 0 )
1513	pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + r1_hu(ji,jj,Kmm) * r1_2( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) zu_i(ji,jj)
1514	pv_RHSi(ji,jj) = pv_RHSi(ji,jj) + r1_hv(ji,jj,Kmm) * r1_2( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) zv_i(ji,jj)
1515	END_2D
1516	END IF
1517	!
1518	! !== TOP stress contribution from baroclinic velocities ==! (no W/D case)
1519	!
1520	IF( ln_isfcav.OR.ln_drgice_imp ) THEN
1521	!
1522	IF( ln_bt_fw ) THEN ! FORWARD integration: use NOW top baroclinic velocity
1523	DO_2D( 0, 0, 0, 0 )
1524	iktu = miku(ji,jj)
1525	iktv = mikv(ji,jj)
1526	zu_i(ji,jj) = puu(ji,jj,iktu,Kmm) - puu_b(ji,jj,Kmm)
1527	zv_i(ji,jj) = pvv(ji,jj,iktv,Kmm) - pvv_b(ji,jj,Kmm)
1528	END_2D
1529	ELSE ! CENTRED integration: use BEFORE top baroclinic velocity
1530	DO_2D( 0, 0, 0, 0 )
1531	iktu = miku(ji,jj)
1532	iktv = mikv(ji,jj)
1533	zu_i(ji,jj) = puu(ji,jj,iktu,Kbb) - puu_b(ji,jj,Kbb)
1534	zv_i(ji,jj) = pvv(ji,jj,iktv,Kbb) - pvv_b(ji,jj,Kbb)
1535	END_2D
1536	ENDIF
1537	! ! use "unclipped" top drag (even if explicit friction is used in 3D calculation)
1538	DO_2D( 0, 0, 0, 0 )
1539	pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + r1_hu(ji,jj,Kmm) * r1_2( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) zu_i(ji,jj)
1540	pv_RHSi(ji,jj) = pv_RHSi(ji,jj) + r1_hv(ji,jj,Kmm) * r1_2( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) zv_i(ji,jj)
1541	END_2D
1542	!
1543	ENDIF
1544	!
1545	END SUBROUTINE dyn_drg_init
1546
1547
1548	SUBROUTINE ts_bck_interp( jn, ll_init, & ! <== in
1549	& za0, za1, za2, za3 ) ! ==> out
1550	!!----------------------------------------------------------------------
1551	INTEGER ,INTENT(in ) :: jn ! index of sub time step
1552	LOGICAL ,INTENT(in ) :: ll_init !
1553	REAL(wp),INTENT( out) :: za0, za1, za2, za3 ! Half-step back interpolation coefficient
1554	!
1555	REAL(wp) :: zepsilon, zgamma ! - -
1556	!!----------------------------------------------------------------------
1557	! ! set Half-step back interpolation coefficient
1558	IF ( jn==1 .AND. ll_init ) THEN !* Forward-backward
1559	za0 = 1._wp
1560	za1 = 0._wp
1561	za2 = 0._wp
1562	za3 = 0._wp
1563	ELSEIF( jn==2 .AND. ll_init ) THEN !* AB2-AM3 Coefficients; bet=0 ; gam=-1/6 ; eps=1/12
1564	za0 = 1.0833333333333_wp ! za0 = 1-gam-eps
1565	za1 =-0.1666666666666_wp ! za1 = gam
1566	za2 = 0.0833333333333_wp ! za2 = eps
1567	za3 = 0._wp
1568	ELSE !* AB3-AM4 Coefficients; bet=0.281105 ; eps=0.013 ; gam=0.0880
1569	IF( rn_bt_alpha == 0._wp ) THEN ! Time diffusion
1570	za0 = 0.614_wp ! za0 = 1/2 + gam + 2*eps
1571	za1 = 0.285_wp ! za1 = 1/2 - 2gam - 3eps
1572	za2 = 0.088_wp ! za2 = gam
1573	za3 = 0.013_wp ! za3 = eps
1574	ELSE ! no time diffusion
1575	zepsilon = 0.00976186_wp - 0.13451357_wp * rn_bt_alpha
1576	zgamma = 0.08344500_wp - 0.51358400_wp * rn_bt_alpha
1577	za0 = 0.5_wp + zgamma + 2._wp * rn_bt_alpha + 2._wp * zepsilon
1578	za1 = 1._wp - za0 - zgamma - zepsilon
1579	za2 = zgamma
1580	za3 = zepsilon
1581	ENDIF
1582	ENDIF
1583	END SUBROUTINE ts_bck_interp
1584
1585	!!======================================================================
1586	END MODULE dynspg_ts

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

New URL for NEMO forge! http://forge.nemo-ocean.eu

Context Navigation

source: NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/DYN/dynspg_ts.F90

Download in other formats: