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/dev_r11879_ENHANCE-05_SimonM-Harmonic_Analysis/src/OCE/DYN – NEMO

Context Navigation

source: NEMO/branches/2019/dev_r11879_ENHANCE-05_SimonM-Harmonic_Analysis/src/OCE/DYN/dynspg_ts.F90 @ 12117

Last change on this file since 12117 was 12117, checked in by smueller, 4 years ago
Merging of the developments in ^{/NEMO/branches/2019/dev_r10742_ENHANCE-12_SimonM-Tides@12096 with respect to}/NEMO/trunk@12072 into ^{/NEMO/branches/2019/dev_r11879_ENHANCE-05_SimonM-Harmonic_Analysis (tickets #2175 and #2194)}
Property svn:keywords set to `Id`
File size: 81.6 KB

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

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

Download in other formats: