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/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/DYN – NEMO

Context Navigation

source: NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/DYN/dynspg_ts.F90 @ 13427

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

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

Download in other formats: