Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

dynspg_ts.F90 @ 13652

Last change on this file since 13652 was 13652, checked in by cguiavarch, 4 years ago
UKMO/NEMO_4.0.2_momentum_trends : code changes (equivalent to rev 13500 of NEMO 4.0.1 branch).
File size: 84.6 KB

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

Download in other formats: