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/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/DYN – NEMO

Context Navigation

source: NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/DYN/dynspg_ts.F90 @ 12708

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

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

Download in other formats: