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_r11943_MERGE_2019/src/OCE/DYN – NEMO

Context Navigation

source: NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/DYN/dynspg_ts.F90 @ 12250

Last change on this file since 12250 was 12250, checked in by smasson, 4 years ago
rev12240_dev_r11943_MERGE_2019: same as [12246], add modifications from dev_r12114_ticket_2263, results unchanged except SPITZ12 as explained in #2263
Property svn:keywords set to `Id`
File size: 82.4 KB

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

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

Download in other formats: