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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

dynspg_ts.F90 in NEMO/branches/2019/ENHANCE-02_ISF_nemo_TEST_MERGE/src/OCE/DYN – NEMO

Context Navigation

source: NEMO/branches/2019/ENHANCE-02_ISF_nemo_TEST_MERGE/src/OCE/DYN/dynspg_ts.F90 @ 11970

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

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

Download in other formats: