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_r12377_KERNEL-06_techene_e3/src/OCE/DYN – NEMO

Context Navigation

source: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DYN/dynspg_ts.F90 @ 12724

Last change on this file since 12724 was 12724, checked in by techene, 4 years ago
branch KERNEL-06 : merge with trunk@12698 #2385 - in duplcated files : changes to comply to the new trunk variables and some loop bug fixes
Property svn:keywords set to `Id`
File size: 78.5 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( 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_00_00
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_00_00
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_00_00
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_00_00
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_00_00
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 atmospheric pressure forcing =!
309	! ! ---------------------------------- !
310	IF( ln_bt_fw ) THEN ! Add wind forcing
311	DO_2D_00_00
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_00_00
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(:,:,Kmm)
389	zhvp2_e(:,:) = hv(:,:,Kmm)
390	zhtp2_e(:,:) = ht(:,:)
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	DO_2D_11_10
478	zhup2_e(ji,jj) = hu_0(ji,jj) + r1_2 * r1_e1e2u(ji,jj) &
479	& * ( e1e2t(ji ,jj) * zsshp2_e(ji ,jj) &
480	& + e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj) ) * ssumask(ji,jj)
481	END_2D
482	DO_2D_10_11
483	zhvp2_e(ji,jj) = hv_0(ji,jj) + r1_2 * r1_e1e2v(ji,jj) &
484	& * ( e1e2t(ji,jj ) * zsshp2_e(ji,jj ) &
485	& + e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1) ) * ssvmask(ji,jj)
486	END_2D
487	!
488	ENDIF
489	!
490	! !== after SSH ==! (jn+1)
491	!
492	! ! update (ua_e,va_e) to enforce volume conservation at open boundaries
493	! ! values of zhup2_e and zhvp2_e on the halo are not needed in bdy_vol2d
494	IF( ln_bdy .AND. ln_vol ) CALL bdy_vol2d( kt, jn, ua_e, va_e, zhup2_e, zhvp2_e )
495	!
496	! ! resulting flux at mid-step (not over the full domain)
497	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
498	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
499	!
500	#if defined key_agrif
501	! Set fluxes during predictor step to ensure volume conservation
502	IF( .NOT.Agrif_Root() .AND. ln_bt_fw ) CALL agrif_dyn_ts_flux( jn, zhU, zhV )
503	#endif
504	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
505
506	IF( ln_wd_dl ) THEN ! un_e and vn_e are set to zero at faces where
507	! ! the direction of the flow is from dry cells
508	CALL wad_Umsk( ztwdmask, zhU, zhV, un_e, vn_e, zuwdmask, zvwdmask ) ! not jpi colomn for U, not jpj row for V
509	!
510	ENDIF
511	!
512	!
513	! Compute Sea Level at step jit+1
514	!-- m+1 m m+1/2 --!
515	!-- ssh = ssh - delta_t' * [ frc + div( flux ) ] --!
516	!-------------------------------------------------------------------------!
517	DO_2D_00_00
518	zhdiv = ( zhU(ji,jj) - zhU(ji-1,jj) + zhV(ji,jj) - zhV(ji,jj-1) ) * r1_e1e2t(ji,jj)
519	ssha_e(ji,jj) = ( sshn_e(ji,jj) - rDt_e * ( zssh_frc(ji,jj) + zhdiv ) ) * ssmask(ji,jj)
520	END_2D
521	!
522	CALL lbc_lnk_multi( 'dynspg_ts', ssha_e, 'T', 1._wp, zhU, 'U', -1._wp, zhV, 'V', -1._wp )
523	!
524	! ! Sum over sub-time-steps to compute advective velocities
525	za2 = wgtbtp2(jn) ! zhU, zhV hold fluxes extrapolated at jn+0.5
526	un_adv(:,:) = un_adv(:,:) + za2 * zhU(:,:) * r1_e2u(:,:)
527	vn_adv(:,:) = vn_adv(:,:) + za2 * zhV(:,:) * r1_e1v(:,:)
528	! sum over sub-time-steps to decide which baroclinic velocities to set to zero (zuwdav2 is only used when ln_wd_dl_bc=True)
529	IF ( ln_wd_dl_bc ) THEN
530	zuwdav2(1:jpim1,1:jpj ) = zuwdav2(1:jpim1,1:jpj ) + za2 * zuwdmask(1:jpim1,1:jpj ) ! not jpi-column
531	zvwdav2(1:jpi ,1:jpjm1) = zvwdav2(1:jpi ,1:jpjm1) + za2 * zvwdmask(1:jpi ,1:jpjm1) ! not jpj-row
532	END IF
533	!
534	! Duplicate sea level across open boundaries (this is only cosmetic if linssh=T)
535	IF( ln_bdy ) CALL bdy_ssh( ssha_e )
536	#if defined key_agrif
537	IF( .NOT.Agrif_Root() ) CALL agrif_ssh_ts( jn )
538	#endif
539	!
540	! Sea Surface Height at u-,v-points (vvl case only)
541	IF( .NOT.ln_linssh ) THEN
542	DO_2D_00_00
543	zsshu_a(ji,jj) = r1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj) &
544	& * ( e1e2t(ji ,jj ) * ssha_e(ji ,jj ) &
545	& + e1e2t(ji+1,jj ) * ssha_e(ji+1,jj ) )
546	zsshv_a(ji,jj) = r1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj) &
547	& * ( e1e2t(ji ,jj ) * ssha_e(ji ,jj ) &
548	& + e1e2t(ji ,jj+1) * ssha_e(ji ,jj+1) )
549	END_2D
550	ENDIF
551	!
552	! Half-step back interpolation of SSH for surface pressure computation at step jit+1/2
553	!-- m+1/2 m+1 m m-1 m-2 --!
554	!-- ssh' = za0 * ssh + za1 * ssh + za2 * ssh + za3 * ssh --!
555	!------------------------------------------------------------------------------------------!
556	CALL ts_bck_interp( jn, ll_init, za0, za1, za2, za3 ) ! coeficients of the interpolation
557	zsshp2_e(:,:) = za0 * ssha_e(:,:) + za1 * sshn_e (:,:) &
558	& + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:)
559	!
560	! ! Surface pressure gradient
561	zldg = ( 1._wp - rn_scal_load ) * grav ! local factor
562	DO_2D_00_00
563	zu_spg(ji,jj) = - zldg * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj)
564	zv_spg(ji,jj) = - zldg * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj)
565	END_2D
566	IF( ln_wd_il ) THEN ! W/D : gravity filters applied on pressure gradient
567	CALL wad_spg( zsshp2_e, zcpx, zcpy ) ! Calculating W/D gravity filters
568	zu_spg(2:jpim1,2:jpjm1) = zu_spg(2:jpim1,2:jpjm1) * zcpx(2:jpim1,2:jpjm1)
569	zv_spg(2:jpim1,2:jpjm1) = zv_spg(2:jpim1,2:jpjm1) * zcpy(2:jpim1,2:jpjm1)
570	ENDIF
571	!
572	! Add Coriolis trend:
573	! zwz array below or triads normally depend on sea level with ln_linssh=F and should be updated
574	! at each time step. We however keep them constant here for optimization.
575	! Recall that zhU and zhV hold fluxes at jn+0.5 (extrapolated not backward interpolated)
576	CALL dyn_cor_2d( zhup2_e, zhvp2_e, ua_e, va_e, zhU, zhV, zu_trd, zv_trd )
577	!
578	! Add tidal astronomical forcing if defined
579	IF ( ln_tide .AND. ln_tide_pot ) THEN
580	DO_2D_00_00
581	zu_trd(ji,jj) = zu_trd(ji,jj) + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) * r1_e1u(ji,jj)
582	zv_trd(ji,jj) = zv_trd(ji,jj) + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) * r1_e2v(ji,jj)
583	END_2D
584	ENDIF
585	!
586	! Add bottom stresses:
587	!jth do implicitly instead
588	IF ( .NOT. ll_wd ) THEN ! Revert to explicit for bit comparison tests in non wad runs
589	DO_2D_00_00
590	zu_trd(ji,jj) = zu_trd(ji,jj) + zCdU_u(ji,jj) * un_e(ji,jj) * hur_e(ji,jj)
591	zv_trd(ji,jj) = zv_trd(ji,jj) + zCdU_v(ji,jj) * vn_e(ji,jj) * hvr_e(ji,jj)
592	END_2D
593	ENDIF
594	!
595	! Set next velocities:
596	! Compute barotropic speeds at step jit+1 (h : total height of the water colomn)
597	!-- VECTOR FORM
598	!-- m+1 m / m+1/2 \ --!
599	!-- u = u + delta_t' * \ (1-r)g grad_x( ssh') - f * k vect u + frc / --!
600	!-- --!
601	!-- FLUX FORM --!
602	!-- m+1 __1__ / m m / m+1/2 m+1/2 m+1/2 n \ \ --!
603	!-- u = m+1 \| h * u + delta_t' * \ h * (1-r)g grad_x( ssh') - h * f * k vect u + h * frc / \| --!
604	!-- h \ / --!
605	!------------------------------------------------------------------------------------------------------------------------!
606	IF( ln_dynadv_vec .OR. ln_linssh ) THEN !* Vector form
607	DO_2D_00_00
608	ua_e(ji,jj) = ( un_e(ji,jj) &
609	& + rDt_e * ( zu_spg(ji,jj) &
610	& + zu_trd(ji,jj) &
611	& + zu_frc(ji,jj) ) &
612	& ) * ssumask(ji,jj)
613
614	va_e(ji,jj) = ( vn_e(ji,jj) &
615	& + rDt_e * ( zv_spg(ji,jj) &
616	& + zv_trd(ji,jj) &
617	& + zv_frc(ji,jj) ) &
618	& ) * ssvmask(ji,jj)
619	END_2D
620	!
621	ELSE !* Flux form
622	DO_2D_00_00
623	! ! hu_e, hv_e hold depth at jn, zhup2_e, zhvp2_e hold extrapolated depth at jn+1/2
624	! ! backward interpolated depth used in spg terms at jn+1/2
625	zhu_bck = hu_0(ji,jj) + r1_2r1_e1e2u(ji,jj) ( e1e2t(ji ,jj) * zsshp2_e(ji ,jj) &
626	& + e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj) ) * ssumask(ji,jj)
627	zhv_bck = hv_0(ji,jj) + r1_2r1_e1e2v(ji,jj) ( e1e2t(ji,jj ) * zsshp2_e(ji,jj ) &
628	& + e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1) ) * ssvmask(ji,jj)
629	! ! inverse depth at jn+1
630	z1_hu = ssumask(ji,jj) / ( hu_0(ji,jj) + zsshu_a(ji,jj) + 1._wp - ssumask(ji,jj) )
631	z1_hv = ssvmask(ji,jj) / ( hv_0(ji,jj) + zsshv_a(ji,jj) + 1._wp - ssvmask(ji,jj) )
632	!
633	ua_e(ji,jj) = ( hu_e (ji,jj) * un_e (ji,jj) &
634	& + rDt_e * ( zhu_bck * zu_spg (ji,jj) & !
635	& + zhup2_e(ji,jj) * zu_trd (ji,jj) & !
636	& + hu(ji,jj,Kmm) * zu_frc (ji,jj) ) ) * z1_hu
637	!
638	va_e(ji,jj) = ( hv_e (ji,jj) * vn_e (ji,jj) &
639	& + rDt_e * ( zhv_bck * zv_spg (ji,jj) & !
640	& + zhvp2_e(ji,jj) * zv_trd (ji,jj) & !
641	& + hv(ji,jj,Kmm) * zv_frc (ji,jj) ) ) * z1_hv
642	END_2D
643	ENDIF
644	!jth implicit bottom friction:
645	IF ( ll_wd ) THEN ! revert to explicit for bit comparison tests in non wad runs
646	DO_2D_00_00
647	ua_e(ji,jj) = ua_e(ji,jj) /(1.0 - rDt_e * zCdU_u(ji,jj) * hur_e(ji,jj))
648	va_e(ji,jj) = va_e(ji,jj) /(1.0 - rDt_e * zCdU_v(ji,jj) * hvr_e(ji,jj))
649	END_2D
650	ENDIF
651
652	IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only)
653	hu_e (2:jpim1,2:jpjm1) = hu_0(2:jpim1,2:jpjm1) + zsshu_a(2:jpim1,2:jpjm1)
654	hur_e(2:jpim1,2:jpjm1) = ssumask(2:jpim1,2:jpjm1) / ( hu_e(2:jpim1,2:jpjm1) + 1._wp - ssumask(2:jpim1,2:jpjm1) )
655	hv_e (2:jpim1,2:jpjm1) = hv_0(2:jpim1,2:jpjm1) + zsshv_a(2:jpim1,2:jpjm1)
656	hvr_e(2:jpim1,2:jpjm1) = ssvmask(2:jpim1,2:jpjm1) / ( hv_e(2:jpim1,2:jpjm1) + 1._wp - ssvmask(2:jpim1,2:jpjm1) )
657	CALL lbc_lnk_multi( 'dynspg_ts', ua_e , 'U', -1._wp, va_e , 'V', -1._wp &
658	& , hu_e , 'U', 1._wp, hv_e , 'V', 1._wp &
659	& , hur_e, 'U', 1._wp, hvr_e, 'V', 1._wp )
660	ELSE
661	CALL lbc_lnk_multi( 'dynspg_ts', ua_e , 'U', -1._wp, va_e , 'V', -1._wp )
662	ENDIF
663	!
664	!
665	! ! open boundaries
666	IF( ln_bdy ) CALL bdy_dyn2d( jn, ua_e, va_e, un_e, vn_e, hur_e, hvr_e, ssha_e )
667	#if defined key_agrif
668	IF( .NOT.Agrif_Root() ) CALL agrif_dyn_ts( jn ) ! Agrif
669	#endif
670	! !* Swap
671	! ! ----
672	ubb_e (:,:) = ub_e (:,:)
673	ub_e (:,:) = un_e (:,:)
674	un_e (:,:) = ua_e (:,:)
675	!
676	vbb_e (:,:) = vb_e (:,:)
677	vb_e (:,:) = vn_e (:,:)
678	vn_e (:,:) = va_e (:,:)
679	!
680	sshbb_e(:,:) = sshb_e(:,:)
681	sshb_e (:,:) = sshn_e(:,:)
682	sshn_e (:,:) = ssha_e(:,:)
683
684	! !* Sum over whole bt loop
685	! ! ----------------------
686	za1 = wgtbtp1(jn)
687	IF( ln_dynadv_vec .OR. ln_linssh ) THEN ! Sum velocities
688	puu_b (:,:,Kaa) = puu_b (:,:,Kaa) + za1 * ua_e (:,:)
689	pvv_b (:,:,Kaa) = pvv_b (:,:,Kaa) + za1 * va_e (:,:)
690	ELSE ! Sum transports
691	IF ( .NOT.ln_wd_dl ) THEN
692	puu_b (:,:,Kaa) = puu_b (:,:,Kaa) + za1 * ua_e (:,:) * hu_e (:,:)
693	pvv_b (:,:,Kaa) = pvv_b (:,:,Kaa) + za1 * va_e (:,:) * hv_e (:,:)
694	ELSE
695	puu_b (:,:,Kaa) = puu_b (:,:,Kaa) + za1 * ua_e (:,:) * hu_e (:,:) * zuwdmask(:,:)
696	pvv_b (:,:,Kaa) = pvv_b (:,:,Kaa) + za1 * va_e (:,:) * hv_e (:,:) * zvwdmask(:,:)
697	END IF
698	ENDIF
699	! ! Sum sea level
700	pssh(:,:,Kaa) = pssh(:,:,Kaa) + za1 * ssha_e(:,:)
701
702	! ! ==================== !
703	END DO ! end loop !
704	! ! ==================== !
705	! -----------------------------------------------------------------------------
706	! Phase 3. update the general trend with the barotropic trend
707	! -----------------------------------------------------------------------------
708	!
709	! Set advection velocity correction:
710	IF (ln_bt_fw) THEN
711	IF( .NOT.( kt == nit000 .AND. l_1st_euler ) ) THEN
712	DO_2D_11_11
713	zun_save = un_adv(ji,jj)
714	zvn_save = vn_adv(ji,jj)
715	! ! apply the previously computed correction
716	un_adv(ji,jj) = r1_2 * ( ub2_b(ji,jj) + zun_save - rn_atfp * un_bf(ji,jj) )
717	vn_adv(ji,jj) = r1_2 * ( vb2_b(ji,jj) + zvn_save - rn_atfp * vn_bf(ji,jj) )
718	! ! Update corrective fluxes for next time step
719	un_bf(ji,jj) = rn_atfp * un_bf(ji,jj) + ( zun_save - ub2_b(ji,jj) )
720	vn_bf(ji,jj) = rn_atfp * vn_bf(ji,jj) + ( zvn_save - vb2_b(ji,jj) )
721	! ! Save integrated transport for next computation
722	ub2_b(ji,jj) = zun_save
723	vb2_b(ji,jj) = zvn_save
724	END_2D
725	ELSE
726	un_bf(:,:) = 0._wp ! corrective fluxes for next time step set to zero
727	vn_bf(:,:) = 0._wp
728	ub2_b(:,:) = un_adv(:,:) ! Save integrated transport for next computation
729	vb2_b(:,:) = vn_adv(:,:)
730	END IF
731	ENDIF
732
733
734	!
735	! Update barotropic trend:
736	IF( ln_dynadv_vec .OR. ln_linssh ) THEN
737	DO jk=1,jpkm1
738	puu(:,:,jk,Krhs) = puu(:,:,jk,Krhs) + ( puu_b(:,:,Kaa) - puu_b(:,:,Kbb) ) * r1_Dt_b
739	pvv(:,:,jk,Krhs) = pvv(:,:,jk,Krhs) + ( pvv_b(:,:,Kaa) - pvv_b(:,:,Kbb) ) * r1_Dt_b
740	END DO
741	ELSE
742	! At this stage, pssh(:,:,:,Krhs) has been corrected: compute new depths at velocity points
743	DO_2D_10_10
744	zsshu_a(ji,jj) = r1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj) &
745	& * ( e1e2t(ji ,jj) * pssh(ji ,jj,Kaa) &
746	& + e1e2t(ji+1,jj) * pssh(ji+1,jj,Kaa) )
747	zsshv_a(ji,jj) = r1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj) &
748	& * ( e1e2t(ji,jj ) * pssh(ji,jj ,Kaa) &
749	& + e1e2t(ji,jj+1) * pssh(ji,jj+1,Kaa) )
750	END_2D
751	CALL lbc_lnk_multi( 'dynspg_ts', zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) ! Boundary conditions
752	!
753	DO jk=1,jpkm1
754	puu(:,:,jk,Krhs) = puu(:,:,jk,Krhs) + r1_hu(:,:,Kmm) * ( puu_b(:,:,Kaa) - puu_b(:,:,Kbb) * hu(:,:,Kbb) ) * r1_Dt_b
755	pvv(:,:,jk,Krhs) = pvv(:,:,jk,Krhs) + r1_hv(:,:,Kmm) * ( pvv_b(:,:,Kaa) - pvv_b(:,:,Kbb) * hv(:,:,Kbb) ) * r1_Dt_b
756	END DO
757	! Save barotropic velocities not transport:
758	puu_b(:,:,Kaa) = puu_b(:,:,Kaa) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - ssumask(:,:) )
759	pvv_b(:,:,Kaa) = pvv_b(:,:,Kaa) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - ssvmask(:,:) )
760	ENDIF
761
762
763	! Correct velocities so that the barotropic velocity equals (un_adv, vn_adv) (in all cases)
764	DO jk = 1, jpkm1
765	puu(:,:,jk,Kmm) = ( puu(:,:,jk,Kmm) + un_adv(:,:)r1_hu(:,:,Kmm) - puu_b(:,:,Kmm) ) umask(:,:,jk)
766	pvv(:,:,jk,Kmm) = ( pvv(:,:,jk,Kmm) + vn_adv(:,:)r1_hv(:,:,Kmm) - pvv_b(:,:,Kmm) ) vmask(:,:,jk)
767	END DO
768
769	IF ( ln_wd_dl .and. ln_wd_dl_bc) THEN
770	DO jk = 1, jpkm1
771	puu(:,:,jk,Kmm) = ( un_adv(:,:)*r1_hu(:,:,Kmm) &
772	& + zuwdav2(:,:)(puu(:,:,jk,Kmm) - un_adv(:,:)r1_hu(:,:,Kmm)) ) * umask(:,:,jk)
773	pvv(:,:,jk,Kmm) = ( vn_adv(:,:)*r1_hv(:,:,Kmm) &
774	& + zvwdav2(:,:)(pvv(:,:,jk,Kmm) - vn_adv(:,:)r1_hv(:,:,Kmm)) ) * vmask(:,:,jk)
775	END DO
776	END IF
777
778
779	CALL iom_put( "ubar", un_adv(:,:)*r1_hu(:,:,Kmm) ) ! barotropic i-current
780	CALL iom_put( "vbar", vn_adv(:,:)*r1_hv(:,:,Kmm) ) ! barotropic i-current
781	!
782	#if defined key_agrif
783	! Save time integrated fluxes during child grid integration
784	! (used to update coarse grid transports at next time step)
785	!
786	IF( .NOT.Agrif_Root() .AND. ln_bt_fw ) THEN
787	IF( Agrif_NbStepint() == 0 ) THEN
788	ub2_i_b(:,:) = 0._wp
789	vb2_i_b(:,:) = 0._wp
790	END IF
791	!
792	za1 = 1._wp / REAL(Agrif_rhot(), wp)
793	ub2_i_b(:,:) = ub2_i_b(:,:) + za1 * ub2_b(:,:)
794	vb2_i_b(:,:) = vb2_i_b(:,:) + za1 * vb2_b(:,:)
795	ENDIF
796	#endif
797	! !* write time-spliting arrays in the restart
798	IF( lrst_oce .AND.ln_bt_fw ) CALL ts_rst( kt, 'WRITE' )
799	!
800	IF( ln_wd_il ) DEALLOCATE( zcpx, zcpy )
801	IF( ln_wd_dl ) DEALLOCATE( ztwdmask, zuwdmask, zvwdmask, zuwdav2, zvwdav2 )
802	!
803	CALL iom_put( "baro_u" , puu_b(:,:,Kmm) ) ! Barotropic U Velocity
804	CALL iom_put( "baro_v" , pvv_b(:,:,Kmm) ) ! Barotropic V Velocity
805	!
806	END SUBROUTINE dyn_spg_ts
807
808
809	SUBROUTINE ts_wgt( ll_av, ll_fw, jpit, zwgt1, zwgt2)
810	!!---------------------------------------------------------------------
811	!! * ROUTINE ts_wgt *
812	!!
813	!! ** Purpose : Set time-splitting weights for temporal averaging (or not)
814	!!----------------------------------------------------------------------
815	LOGICAL, INTENT(in) :: ll_av ! temporal averaging=.true.
816	LOGICAL, INTENT(in) :: ll_fw ! forward time splitting =.true.
817	INTEGER, INTENT(inout) :: jpit ! cycle length
818	REAL(wp), DIMENSION(3*nn_e), INTENT(inout) :: zwgt1, & ! Primary weights
819	zwgt2 ! Secondary weights
820
821	INTEGER :: jic, jn, ji ! temporary integers
822	REAL(wp) :: za1, za2
823	!!----------------------------------------------------------------------
824
825	zwgt1(:) = 0._wp
826	zwgt2(:) = 0._wp
827
828	! Set time index when averaged value is requested
829	IF (ll_fw) THEN
830	jic = nn_e
831	ELSE
832	jic = 2 * nn_e
833	ENDIF
834
835	! Set primary weights:
836	IF (ll_av) THEN
837	! Define simple boxcar window for primary weights
838	! (width = nn_e, centered around jic)
839	SELECT CASE ( nn_bt_flt )
840	CASE( 0 ) ! No averaging
841	zwgt1(jic) = 1._wp
842	jpit = jic
843
844	CASE( 1 ) ! Boxcar, width = nn_e
845	DO jn = 1, 3*nn_e
846	za1 = ABS(float(jn-jic))/float(nn_e)
847	IF (za1 < 0.5_wp) THEN
848	zwgt1(jn) = 1._wp
849	jpit = jn
850	ENDIF
851	ENDDO
852
853	CASE( 2 ) ! Boxcar, width = 2 * nn_e
854	DO jn = 1, 3*nn_e
855	za1 = ABS(float(jn-jic))/float(nn_e)
856	IF (za1 < 1._wp) THEN
857	zwgt1(jn) = 1._wp
858	jpit = jn
859	ENDIF
860	ENDDO
861	CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_bt_flt' )
862	END SELECT
863
864	ELSE ! No time averaging
865	zwgt1(jic) = 1._wp
866	jpit = jic
867	ENDIF
868
869	! Set secondary weights
870	DO jn = 1, jpit
871	DO ji = jn, jpit
872	zwgt2(jn) = zwgt2(jn) + zwgt1(ji)
873	END DO
874	END DO
875
876	! Normalize weigths:
877	za1 = 1._wp / SUM(zwgt1(1:jpit))
878	za2 = 1._wp / SUM(zwgt2(1:jpit))
879	DO jn = 1, jpit
880	zwgt1(jn) = zwgt1(jn) * za1
881	zwgt2(jn) = zwgt2(jn) * za2
882	END DO
883	!
884	END SUBROUTINE ts_wgt
885
886
887	SUBROUTINE ts_rst( kt, cdrw )
888	!!---------------------------------------------------------------------
889	!! * ROUTINE ts_rst *
890	!!
891	!! ** Purpose : Read or write time-splitting arrays in restart file
892	!!----------------------------------------------------------------------
893	INTEGER , INTENT(in) :: kt ! ocean time-step
894	CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag
895	!!----------------------------------------------------------------------
896	!
897	IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialise
898	! ! ---------------
899	IF( ln_rstart .AND. ln_bt_fw .AND. (.NOT.l_1st_euler) ) THEN !* Read the restart file
900	CALL iom_get( numror, jpdom_autoglo, 'ub2_b' , ub2_b (:,:), ldxios = lrxios )
901	CALL iom_get( numror, jpdom_autoglo, 'vb2_b' , vb2_b (:,:), ldxios = lrxios )
902	CALL iom_get( numror, jpdom_autoglo, 'un_bf' , un_bf (:,:), ldxios = lrxios )
903	CALL iom_get( numror, jpdom_autoglo, 'vn_bf' , vn_bf (:,:), ldxios = lrxios )
904	IF( .NOT.ln_bt_av ) THEN
905	CALL iom_get( numror, jpdom_autoglo, 'sshbb_e' , sshbb_e(:,:), ldxios = lrxios )
906	CALL iom_get( numror, jpdom_autoglo, 'ubb_e' , ubb_e(:,:), ldxios = lrxios )
907	CALL iom_get( numror, jpdom_autoglo, 'vbb_e' , vbb_e(:,:), ldxios = lrxios )
908	CALL iom_get( numror, jpdom_autoglo, 'sshb_e' , sshb_e(:,:), ldxios = lrxios )
909	CALL iom_get( numror, jpdom_autoglo, 'ub_e' , ub_e(:,:), ldxios = lrxios )
910	CALL iom_get( numror, jpdom_autoglo, 'vb_e' , vb_e(:,:), ldxios = lrxios )
911	ENDIF
912	#if defined key_agrif
913	! Read time integrated fluxes
914	IF ( .NOT.Agrif_Root() ) THEN
915	CALL iom_get( numror, jpdom_autoglo, 'ub2_i_b' , ub2_i_b(:,:), ldxios = lrxios )
916	CALL iom_get( numror, jpdom_autoglo, 'vb2_i_b' , vb2_i_b(:,:), ldxios = lrxios )
917	ENDIF
918	#endif
919	ELSE !* Start from rest
920	IF(lwp) WRITE(numout,*)
921	IF(lwp) WRITE(numout,*) ' ==>>> start from rest: set barotropic values to 0'
922	ub2_b (:,:) = 0._wp ; vb2_b (:,:) = 0._wp ! used in the 1st interpol of agrif
923	un_adv(:,:) = 0._wp ; vn_adv(:,:) = 0._wp ! used in the 1st interpol of agrif
924	un_bf (:,:) = 0._wp ; vn_bf (:,:) = 0._wp ! used in the 1st update of agrif
925	#if defined key_agrif
926	IF ( .NOT.Agrif_Root() ) THEN
927	ub2_i_b(:,:) = 0._wp ; vb2_i_b(:,:) = 0._wp ! used in the 1st update of agrif
928	ENDIF
929	#endif
930	ENDIF
931	!
932	ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN ! Create restart file
933	! ! -------------------
934	IF(lwp) WRITE(numout,*) '---- ts_rst ----'
935	IF( lwxios ) CALL iom_swap( cwxios_context )
936	CALL iom_rstput( kt, nitrst, numrow, 'ub2_b' , ub2_b (:,:), ldxios = lwxios )
937	CALL iom_rstput( kt, nitrst, numrow, 'vb2_b' , vb2_b (:,:), ldxios = lwxios )
938	CALL iom_rstput( kt, nitrst, numrow, 'un_bf' , un_bf (:,:), ldxios = lwxios )
939	CALL iom_rstput( kt, nitrst, numrow, 'vn_bf' , vn_bf (:,:), ldxios = lwxios )
940	!
941	IF (.NOT.ln_bt_av) THEN
942	CALL iom_rstput( kt, nitrst, numrow, 'sshbb_e' , sshbb_e(:,:), ldxios = lwxios )
943	CALL iom_rstput( kt, nitrst, numrow, 'ubb_e' , ubb_e(:,:), ldxios = lwxios )
944	CALL iom_rstput( kt, nitrst, numrow, 'vbb_e' , vbb_e(:,:), ldxios = lwxios )
945	CALL iom_rstput( kt, nitrst, numrow, 'sshb_e' , sshb_e(:,:), ldxios = lwxios )
946	CALL iom_rstput( kt, nitrst, numrow, 'ub_e' , ub_e(:,:), ldxios = lwxios )
947	CALL iom_rstput( kt, nitrst, numrow, 'vb_e' , vb_e(:,:), ldxios = lwxios )
948	ENDIF
949	#if defined key_agrif
950	! Save time integrated fluxes
951	IF ( .NOT.Agrif_Root() ) THEN
952	CALL iom_rstput( kt, nitrst, numrow, 'ub2_i_b' , ub2_i_b(:,:), ldxios = lwxios )
953	CALL iom_rstput( kt, nitrst, numrow, 'vb2_i_b' , vb2_i_b(:,:), ldxios = lwxios )
954	ENDIF
955	#endif
956	IF( lwxios ) CALL iom_swap( cxios_context )
957	ENDIF
958	!
959	END SUBROUTINE ts_rst
960
961
962	SUBROUTINE dyn_spg_ts_init
963	!!---------------------------------------------------------------------
964	!! * ROUTINE dyn_spg_ts_init *
965	!!
966	!! ** Purpose : Set time splitting options
967	!!----------------------------------------------------------------------
968	INTEGER :: ji ,jj ! dummy loop indices
969	REAL(wp) :: zxr2, zyr2, zcmax ! local scalar
970	REAL(wp), DIMENSION(jpi,jpj) :: zcu
971	!!----------------------------------------------------------------------
972	!
973	! Max courant number for ext. grav. waves
974	!
975	DO_2D_11_11
976	zxr2 = r1_e1t(ji,jj) * r1_e1t(ji,jj)
977	zyr2 = r1_e2t(ji,jj) * r1_e2t(ji,jj)
978	zcu(ji,jj) = SQRT( grav * MAX(ht_0(ji,jj),0._wp) * (zxr2 + zyr2) )
979	END_2D
980	!
981	zcmax = MAXVAL( zcu(:,:) )
982	CALL mpp_max( 'dynspg_ts', zcmax )
983
984	! Estimate number of iterations to satisfy a max courant number= rn_bt_cmax
985	IF( ln_bt_auto ) nn_e = CEILING( rn_Dt / rn_bt_cmax * zcmax)
986
987	rDt_e = rn_Dt / REAL( nn_e , wp )
988	zcmax = zcmax * rDt_e
989	! Print results
990	IF(lwp) WRITE(numout,*)
991	IF(lwp) WRITE(numout,*) 'dyn_spg_ts_init : split-explicit free surface'
992	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~'
993	IF( ln_bt_auto ) THEN
994	IF(lwp) WRITE(numout,*) ' ln_ts_auto =.true. Automatically set nn_e '
995	IF(lwp) WRITE(numout,*) ' Max. courant number allowed: ', rn_bt_cmax
996	ELSE
997	IF(lwp) WRITE(numout,*) ' ln_ts_auto=.false.: Use nn_e in namelist nn_e = ', nn_e
998	ENDIF
999
1000	IF(ln_bt_av) THEN
1001	IF(lwp) WRITE(numout,*) ' ln_bt_av =.true. ==> Time averaging over nn_e time steps is on '
1002	ELSE
1003	IF(lwp) WRITE(numout,*) ' ln_bt_av =.false. => No time averaging of barotropic variables '
1004	ENDIF
1005	!
1006	!
1007	IF(ln_bt_fw) THEN
1008	IF(lwp) WRITE(numout,*) ' ln_bt_fw=.true. => Forward integration of barotropic variables '
1009	ELSE
1010	IF(lwp) WRITE(numout,*) ' ln_bt_fw =.false.=> Centred integration of barotropic variables '
1011	ENDIF
1012	!
1013	#if defined key_agrif
1014	! Restrict the use of Agrif to the forward case only
1015	!!! IF( .NOT.ln_bt_fw .AND. .NOT.Agrif_Root() ) CALL ctl_stop( 'AGRIF not implemented if ln_bt_fw=.FALSE.' )
1016	#endif
1017	!
1018	IF(lwp) WRITE(numout,*) ' Time filter choice, nn_bt_flt: ', nn_bt_flt
1019	SELECT CASE ( nn_bt_flt )
1020	CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' Dirac'
1021	CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' Boxcar: width = nn_e'
1022	CASE( 2 ) ; IF(lwp) WRITE(numout,) ' Boxcar: width = 2nn_e'
1023	CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_bt_flt: should 0,1, or 2' )
1024	END SELECT
1025	!
1026	IF(lwp) WRITE(numout,*) ' '
1027	IF(lwp) WRITE(numout,*) ' nn_e = ', nn_e
1028	IF(lwp) WRITE(numout,*) ' Barotropic time step [s] is :', rDt_e
1029	IF(lwp) WRITE(numout,*) ' Maximum Courant number is :', zcmax
1030	!
1031	IF(lwp) WRITE(numout,*) ' Time diffusion parameter rn_bt_alpha: ', rn_bt_alpha
1032	IF ((ln_bt_av.AND.nn_bt_flt/=0).AND.(rn_bt_alpha>0._wp)) THEN
1033	CALL ctl_stop( 'dynspg_ts ERROR: if rn_bt_alpha > 0, remove temporal averaging' )
1034	ENDIF
1035	!
1036	IF( .NOT.ln_bt_av .AND. .NOT.ln_bt_fw ) THEN
1037	CALL ctl_stop( 'dynspg_ts ERROR: No time averaging => only forward integration is possible' )
1038	ENDIF
1039	IF( zcmax>0.9_wp ) THEN
1040	CALL ctl_stop( 'dynspg_ts ERROR: Maximum Courant number is greater than 0.9: Inc. nn_e !' )
1041	ENDIF
1042	!
1043	! ! Allocate time-splitting arrays
1044	IF( dyn_spg_ts_alloc() /= 0 ) CALL ctl_stop('STOP', 'dyn_spg_init: failed to allocate dynspg_ts arrays' )
1045	!
1046	! ! read restart when needed
1047	CALL ts_rst( nit000, 'READ' )
1048	!
1049	IF( lwxios ) THEN
1050	! define variables in restart file when writing with XIOS
1051	CALL iom_set_rstw_var_active('ub2_b')
1052	CALL iom_set_rstw_var_active('vb2_b')
1053	CALL iom_set_rstw_var_active('un_bf')
1054	CALL iom_set_rstw_var_active('vn_bf')
1055	!
1056	IF (.NOT.ln_bt_av) THEN
1057	CALL iom_set_rstw_var_active('sshbb_e')
1058	CALL iom_set_rstw_var_active('ubb_e')
1059	CALL iom_set_rstw_var_active('vbb_e')
1060	CALL iom_set_rstw_var_active('sshb_e')
1061	CALL iom_set_rstw_var_active('ub_e')
1062	CALL iom_set_rstw_var_active('vb_e')
1063	ENDIF
1064	#if defined key_agrif
1065	! Save time integrated fluxes
1066	IF ( .NOT.Agrif_Root() ) THEN
1067	CALL iom_set_rstw_var_active('ub2_i_b')
1068	CALL iom_set_rstw_var_active('vb2_i_b')
1069	ENDIF
1070	#endif
1071	ENDIF
1072	!
1073	END SUBROUTINE dyn_spg_ts_init
1074
1075
1076	SUBROUTINE dyn_cor_2D_init( Kmm )
1077	!!---------------------------------------------------------------------
1078	!! * ROUTINE dyn_cor_2D_init *
1079	!!
1080	!! ** Purpose : Set time splitting options
1081	!! Set arrays to remove/compute coriolis trend.
1082	!! Do it once during initialization if volume is fixed, else at each long time step.
1083	!! Note that these arrays are also used during barotropic loop. These are however frozen
1084	!! although they should be updated in the variable volume case. Not a big approximation.
1085	!! To remove this approximation, copy lines below inside barotropic loop
1086	!! and update depths at T-F points (ht and zhf resp.) at each barotropic time step
1087	!!
1088	!! Compute zwz = f / ( height of the water colomn )
1089	!!----------------------------------------------------------------------
1090	INTEGER, INTENT(in) :: Kmm ! Time index
1091	INTEGER :: ji ,jj, jk ! dummy loop indices
1092	REAL(wp) :: z1_ht
1093	REAL(wp), DIMENSION(jpi,jpj) :: zhf
1094	!!----------------------------------------------------------------------
1095	!
1096	SELECT CASE( nvor_scheme )
1097	CASE( np_EEN ) != EEN scheme using e3f energy & enstrophy scheme
1098	SELECT CASE( nn_een_e3f ) !* ff_f/e3 at F-point
1099	CASE ( 0 ) ! original formulation (masked averaging of e3t divided by 4)
1100	DO_2D_10_10
1101	zwz(ji,jj) = ( ht(ji ,jj+1) + ht(ji+1,jj+1) + &
1102	& ht(ji ,jj ) + ht(ji+1,jj ) ) * 0.25_wp
1103	IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj)
1104	END_2D
1105	CASE ( 1 ) ! new formulation (masked averaging of e3t divided by the sum of mask)
1106	DO_2D_10_10
1107	zwz(ji,jj) = ( ht (ji ,jj+1) + ht (ji+1,jj+1) &
1108	& + ht (ji ,jj ) + ht (ji+1,jj ) ) &
1109	& / ( MAX( 1._wp, ssmask(ji ,jj+1) + ssmask(ji+1,jj+1) &
1110	& + ssmask(ji ,jj ) + ssmask(ji+1,jj ) ) )
1111	IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj)
1112	END_2D
1113	END SELECT
1114	CALL lbc_lnk( 'dynspg_ts', zwz, 'F', 1._wp )
1115	!
1116	ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp
1117	DO_2D_01_01
1118	ftne(ji,jj) = zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1)
1119	ftnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj ) + zwz(ji ,jj )
1120	ftse(ji,jj) = zwz(ji ,jj ) + zwz(ji ,jj-1) + zwz(ji-1,jj-1)
1121	ftsw(ji,jj) = zwz(ji ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj )
1122	END_2D
1123	!
1124	CASE( np_EET ) != EEN scheme using e3t energy conserving scheme
1125	ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp
1126	DO_2D_01_01
1127	z1_ht = ssmask(ji,jj) / ( ht(ji,jj) + 1._wp - ssmask(ji,jj) )
1128	ftne(ji,jj) = ( ff_f(ji-1,jj ) + ff_f(ji ,jj ) + ff_f(ji ,jj-1) ) * z1_ht
1129	ftnw(ji,jj) = ( ff_f(ji-1,jj-1) + ff_f(ji-1,jj ) + ff_f(ji ,jj ) ) * z1_ht
1130	ftse(ji,jj) = ( ff_f(ji ,jj ) + ff_f(ji ,jj-1) + ff_f(ji-1,jj-1) ) * z1_ht
1131	ftsw(ji,jj) = ( ff_f(ji ,jj-1) + ff_f(ji-1,jj-1) + ff_f(ji-1,jj ) ) * z1_ht
1132	END_2D
1133	!
1134	CASE( np_ENE, np_ENS , np_MIX ) != all other schemes (ENE, ENS, MIX) except ENT !
1135	!
1136	zwz(:,:) = 0._wp
1137	zhf(:,:) = 0._wp
1138
1139	!!gm assume 0 in both cases (which is almost surely WRONG ! ) as hvatf has been removed
1140	!!gm A priori a better value should be something like :
1141	!!gm zhf(i,j) = masked sum of ht(i,j) , ht(i+1,j) , ht(i,j+1) , (i+1,j+1)
1142	!!gm divided by the sum of the corresponding mask
1143	!!gm
1144	!!
1145	IF( .NOT.ln_sco ) THEN
1146
1147	!!gm agree the JC comment : this should be done in a much clear way
1148
1149	! JC: It not clear yet what should be the depth at f-points over land in z-coordinate case
1150	! Set it to zero for the time being
1151	! IF( rn_hmin < 0._wp ) THEN ; jk = - INT( rn_hmin ) ! from a nb of level
1152	! ELSE ; jk = MINLOC( gdepw_0, mask = gdepw_0 > rn_hmin, dim = 1 ) ! from a depth
1153	! ENDIF
1154	! zhf(:,:) = gdepw_0(:,:,jk+1)
1155	!
1156	ELSE
1157	!
1158	!zhf(:,:) = hbatf(:,:)
1159	DO_2D_10_10
1160	zhf(ji,jj) = ( ht_0 (ji,jj ) + ht_0 (ji+1,jj ) &
1161	& + ht_0 (ji,jj+1) + ht_0 (ji+1,jj+1) ) &
1162	& / MAX( ssmask(ji,jj ) + ssmask(ji+1,jj ) &
1163	& + ssmask(ji,jj+1) + ssmask(ji+1,jj+1) , 1._wp )
1164	END_2D
1165	ENDIF
1166	!
1167	DO jj = 1, jpjm1
1168	zhf(:,jj) = zhf(:,jj) * (1._wp- umask(:,jj,1) * umask(:,jj+1,1))
1169	END DO
1170	!
1171	DO jk = 1, jpkm1
1172	DO jj = 1, jpjm1
1173	zhf(:,jj) = zhf(:,jj) + e3f(:,jj,jk) * umask(:,jj,jk) * umask(:,jj+1,jk)
1174	END DO
1175	END DO
1176	CALL lbc_lnk( 'dynspg_ts', zhf, 'F', 1._wp )
1177	! JC: TBC. hf should be greater than 0
1178	DO_2D_11_11
1179	IF( zhf(ji,jj) /= 0._wp ) zwz(ji,jj) = 1._wp / zhf(ji,jj)
1180	END_2D
1181	zwz(:,:) = ff_f(:,:) * zwz(:,:)
1182	END SELECT
1183
1184	END SUBROUTINE dyn_cor_2d_init
1185
1186
1187
1188	SUBROUTINE dyn_cor_2d( phu, phv, punb, pvnb, zhU, zhV, zu_trd, zv_trd )
1189	!!---------------------------------------------------------------------
1190	!! * ROUTINE dyn_cor_2d *
1191	!!
1192	!! ** Purpose : Compute u and v coriolis trends
1193	!!----------------------------------------------------------------------
1194	INTEGER :: ji ,jj ! dummy loop indices
1195	REAL(wp) :: zx1, zx2, zy1, zy2, z1_hu, z1_hv ! - -
1196	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: phu, phv, punb, pvnb, zhU, zhV
1197	REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: zu_trd, zv_trd
1198	!!----------------------------------------------------------------------
1199	SELECT CASE( nvor_scheme )
1200	CASE( np_ENT ) ! enstrophy conserving scheme (f-point)
1201	DO_2D_00_00
1202	z1_hu = ssumask(ji,jj) / ( phu(ji,jj) + 1._wp - ssumask(ji,jj) )
1203	z1_hv = ssvmask(ji,jj) / ( phv(ji,jj) + 1._wp - ssvmask(ji,jj) )
1204	zu_trd(ji,jj) = + r1_4 * r1_e1e2u(ji,jj) * z1_hu &
1205	& * ( e1e2t(ji+1,jj)ht(ji+1,jj)ff_t(ji+1,jj) * ( pvnb(ji+1,jj) + pvnb(ji+1,jj-1) ) &
1206	& + e1e2t(ji ,jj)ht(ji ,jj)ff_t(ji ,jj) * ( pvnb(ji ,jj) + pvnb(ji ,jj-1) ) )
1207	!
1208	zv_trd(ji,jj) = - r1_4 * r1_e1e2v(ji,jj) * z1_hv &
1209	& * ( e1e2t(ji,jj+1)ht(ji,jj+1)ff_t(ji,jj+1) * ( punb(ji,jj+1) + punb(ji-1,jj+1) ) &
1210	& + e1e2t(ji,jj )ht(ji,jj )ff_t(ji,jj ) * ( punb(ji,jj ) + punb(ji-1,jj ) ) )
1211	END_2D
1212	!
1213	CASE( np_ENE , np_MIX ) ! energy conserving scheme (t-point) ENE or MIX
1214	DO_2D_00_00
1215	zy1 = ( zhV(ji,jj-1) + zhV(ji+1,jj-1) ) * r1_e1u(ji,jj)
1216	zy2 = ( zhV(ji,jj ) + zhV(ji+1,jj ) ) * r1_e1u(ji,jj)
1217	zx1 = ( zhU(ji-1,jj) + zhU(ji-1,jj+1) ) * r1_e2v(ji,jj)
1218	zx2 = ( zhU(ji ,jj) + zhU(ji ,jj+1) ) * r1_e2v(ji,jj)
1219	! energy conserving formulation for planetary vorticity term
1220	zu_trd(ji,jj) = r1_4 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
1221	zv_trd(ji,jj) = - r1_4 * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 )
1222	END_2D
1223	!
1224	CASE( np_ENS ) ! enstrophy conserving scheme (f-point)
1225	DO_2D_00_00
1226	zy1 = r1_8 * ( zhV(ji ,jj-1) + zhV(ji+1,jj-1) &
1227	& + zhV(ji ,jj ) + zhV(ji+1,jj ) ) * r1_e1u(ji,jj)
1228	zx1 = - r1_8 * ( zhU(ji-1,jj ) + zhU(ji-1,jj+1) &
1229	& + zhU(ji ,jj ) + zhU(ji ,jj+1) ) * r1_e2v(ji,jj)
1230	zu_trd(ji,jj) = zy1 * ( zwz(ji ,jj-1) + zwz(ji,jj) )
1231	zv_trd(ji,jj) = zx1 * ( zwz(ji-1,jj ) + zwz(ji,jj) )
1232	END_2D
1233	!
1234	CASE( np_EET , np_EEN ) ! energy & enstrophy scheme (using e3t or e3f)
1235	DO_2D_00_00
1236	zu_trd(ji,jj) = + r1_12 * r1_e1u(ji,jj) * ( ftne(ji,jj ) * zhV(ji ,jj ) &
1237	& + ftnw(ji+1,jj) * zhV(ji+1,jj ) &
1238	& + ftse(ji,jj ) * zhV(ji ,jj-1) &
1239	& + ftsw(ji+1,jj) * zhV(ji+1,jj-1) )
1240	zv_trd(ji,jj) = - r1_12 * r1_e2v(ji,jj) * ( ftsw(ji,jj+1) * zhU(ji-1,jj+1) &
1241	& + ftse(ji,jj+1) * zhU(ji ,jj+1) &
1242	& + ftnw(ji,jj ) * zhU(ji-1,jj ) &
1243	& + ftne(ji,jj ) * zhU(ji ,jj ) )
1244	END_2D
1245	!
1246	END SELECT
1247	!
1248	END SUBROUTINE dyn_cor_2D
1249
1250
1251	SUBROUTINE wad_tmsk( pssh, ptmsk )
1252	!!----------------------------------------------------------------------
1253	!! * ROUTINE wad_lmt *
1254	!!
1255	!! ** Purpose : set wetting & drying mask at tracer points
1256	!! for the current barotropic sub-step
1257	!!
1258	!! ** Method : ???
1259	!!
1260	!! ** Action : ptmsk : wetting & drying t-mask
1261	!!----------------------------------------------------------------------
1262	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pssh !
1263	REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: ptmsk !
1264	!
1265	INTEGER :: ji, jj ! dummy loop indices
1266	!!----------------------------------------------------------------------
1267	!
1268	IF( ln_wd_dl_rmp ) THEN
1269	DO_2D_11_11
1270	IF ( pssh(ji,jj) + ht_0(ji,jj) > 2._wp * rn_wdmin1 ) THEN
1271	! IF ( pssh(ji,jj) + ht_0(ji,jj) > rn_wdmin2 ) THEN
1272	ptmsk(ji,jj) = 1._wp
1273	ELSEIF( pssh(ji,jj) + ht_0(ji,jj) > rn_wdmin1 ) THEN
1274	ptmsk(ji,jj) = TANH( 50._wp( ( pssh(ji,jj) + ht_0(ji,jj) - rn_wdmin1 )r_rn_wdmin1) )
1275	ELSE
1276	ptmsk(ji,jj) = 0._wp
1277	ENDIF
1278	END_2D
1279	ELSE
1280	DO_2D_11_11
1281	IF ( pssh(ji,jj) + ht_0(ji,jj) > rn_wdmin1 ) THEN ; ptmsk(ji,jj) = 1._wp
1282	ELSE ; ptmsk(ji,jj) = 0._wp
1283	ENDIF
1284	END_2D
1285	ENDIF
1286	!
1287	END SUBROUTINE wad_tmsk
1288
1289
1290	SUBROUTINE wad_Umsk( pTmsk, phU, phV, pu, pv, pUmsk, pVmsk )
1291	!!----------------------------------------------------------------------
1292	!! * ROUTINE wad_lmt *
1293	!!
1294	!! ** Purpose : set wetting & drying mask at tracer points
1295	!! for the current barotropic sub-step
1296	!!
1297	!! ** Method : ???
1298	!!
1299	!! ** Action : ptmsk : wetting & drying t-mask
1300	!!----------------------------------------------------------------------
1301	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pTmsk ! W & D t-mask
1302	REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: phU, phV, pu, pv ! ocean velocities and transports
1303	REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pUmsk, pVmsk ! W & D u- and v-mask
1304	!
1305	INTEGER :: ji, jj ! dummy loop indices
1306	!!----------------------------------------------------------------------
1307	!
1308	DO_2D_11_10
1309	IF ( phU(ji,jj) > 0._wp ) THEN ; pUmsk(ji,jj) = pTmsk(ji ,jj)
1310	ELSE ; pUmsk(ji,jj) = pTmsk(ji+1,jj)
1311	ENDIF
1312	phU(ji,jj) = pUmsk(ji,jj)*phU(ji,jj)
1313	pu (ji,jj) = pUmsk(ji,jj)*pu (ji,jj)
1314	END_2D
1315	!
1316	DO_2D_10_11
1317	IF ( phV(ji,jj) > 0._wp ) THEN ; pVmsk(ji,jj) = pTmsk(ji,jj )
1318	ELSE ; pVmsk(ji,jj) = pTmsk(ji,jj+1)
1319	ENDIF
1320	phV(ji,jj) = pVmsk(ji,jj)*phV(ji,jj)
1321	pv (ji,jj) = pVmsk(ji,jj)*pv (ji,jj)
1322	END_2D
1323	!
1324	END SUBROUTINE wad_Umsk
1325
1326
1327	SUBROUTINE wad_spg( pshn, zcpx, zcpy )
1328	!!---------------------------------------------------------------------
1329	!! * ROUTINE wad_sp *
1330	!!
1331	!! ** Purpose :
1332	!!----------------------------------------------------------------------
1333	INTEGER :: ji ,jj ! dummy loop indices
1334	LOGICAL :: ll_tmp1, ll_tmp2
1335	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pshn
1336	REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: zcpx, zcpy
1337	!!----------------------------------------------------------------------
1338	DO_2D_00_00
1339	ll_tmp1 = MIN( pshn(ji,jj) , pshn(ji+1,jj) ) > &
1340	& MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) .AND. &
1341	& MAX( pshn(ji,jj) + ht_0(ji,jj) , pshn(ji+1,jj) + ht_0(ji+1,jj) ) &
1342	& > rn_wdmin1 + rn_wdmin2
1343	ll_tmp2 = ( ABS( pshn(ji+1,jj) - pshn(ji ,jj)) > 1.E-12 ).AND.( &
1344	& MAX( pshn(ji,jj) , pshn(ji+1,jj) ) > &
1345	& MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 )
1346	IF(ll_tmp1) THEN
1347	zcpx(ji,jj) = 1.0_wp
1348	ELSEIF(ll_tmp2) THEN
1349	! no worries about pshn(ji+1,jj) - pshn(ji ,jj) = 0, it won't happen ! here
1350	zcpx(ji,jj) = ABS( (pshn(ji+1,jj) + ht_0(ji+1,jj) - pshn(ji,jj) - ht_0(ji,jj)) &
1351	& / (pshn(ji+1,jj) - pshn(ji ,jj)) )
1352	zcpx(ji,jj) = max(min( zcpx(ji,jj) , 1.0_wp),0.0_wp)
1353	ELSE
1354	zcpx(ji,jj) = 0._wp
1355	ENDIF
1356	!
1357	ll_tmp1 = MIN( pshn(ji,jj) , pshn(ji,jj+1) ) > &
1358	& MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) .AND. &
1359	& MAX( pshn(ji,jj) + ht_0(ji,jj) , pshn(ji,jj+1) + ht_0(ji,jj+1) ) &
1360	& > rn_wdmin1 + rn_wdmin2
1361	ll_tmp2 = ( ABS( pshn(ji,jj) - pshn(ji,jj+1)) > 1.E-12 ).AND.( &
1362	& MAX( pshn(ji,jj) , pshn(ji,jj+1) ) > &
1363	& MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 )
1364
1365	IF(ll_tmp1) THEN
1366	zcpy(ji,jj) = 1.0_wp
1367	ELSE IF(ll_tmp2) THEN
1368	! no worries about pshn(ji,jj+1) - pshn(ji,jj ) = 0, it won't happen ! here
1369	zcpy(ji,jj) = ABS( (pshn(ji,jj+1) + ht_0(ji,jj+1) - pshn(ji,jj) - ht_0(ji,jj)) &
1370	& / (pshn(ji,jj+1) - pshn(ji,jj )) )
1371	zcpy(ji,jj) = MAX( 0._wp , MIN( zcpy(ji,jj) , 1.0_wp ) )
1372	ELSE
1373	zcpy(ji,jj) = 0._wp
1374	ENDIF
1375	END_2D
1376
1377	END SUBROUTINE wad_spg
1378
1379
1380
1381	SUBROUTINE dyn_drg_init( Kbb, Kmm, puu, pvv, puu_b ,pvv_b, pu_RHSi, pv_RHSi, pCdU_u, pCdU_v )
1382	!!----------------------------------------------------------------------
1383	!! * ROUTINE dyn_drg_init *
1384	!!
1385	!! ** Purpose : - add the baroclinic top/bottom drag contribution to
1386	!! the baroclinic part of the barotropic RHS
1387	!! - compute the barotropic drag coefficients
1388	!!
1389	!! ** Method : computation done over the INNER domain only
1390	!!----------------------------------------------------------------------
1391	INTEGER , INTENT(in ) :: Kbb, Kmm ! ocean time level indices
1392	REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(in ) :: puu, pvv ! ocean velocities and RHS of momentum equation
1393	REAL(wp), DIMENSION(jpi,jpj,jpt) , INTENT(in ) :: puu_b, pvv_b ! barotropic velocities at main time levels
1394	REAL(wp), DIMENSION(jpi,jpj) , INTENT(inout) :: pu_RHSi, pv_RHSi ! baroclinic part of the barotropic RHS
1395	REAL(wp), DIMENSION(jpi,jpj) , INTENT( out) :: pCdU_u , pCdU_v ! barotropic drag coefficients
1396	!
1397	INTEGER :: ji, jj ! dummy loop indices
1398	INTEGER :: ikbu, ikbv, iktu, iktv
1399	REAL(wp) :: zztmp
1400	REAL(wp), DIMENSION(jpi,jpj) :: zu_i, zv_i
1401	!!----------------------------------------------------------------------
1402	!
1403	! !== Set the barotropic drag coef. ==!
1404	!
1405	IF( ln_isfcav ) THEN ! top+bottom friction (ocean cavities)
1406
1407	DO_2D_00_00
1408	pCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) + rCdU_top(ji+1,jj)+rCdU_top(ji,jj) )
1409	pCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) + rCdU_top(ji,jj+1)+rCdU_top(ji,jj) )
1410	END_2D
1411	ELSE ! bottom friction only
1412	DO_2D_00_00
1413	pCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) )
1414	pCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) )
1415	END_2D
1416	ENDIF
1417	!
1418	! !== BOTTOM stress contribution from baroclinic velocities ==!
1419	!
1420	IF( ln_bt_fw ) THEN ! FORWARD integration: use NOW bottom baroclinic velocities
1421
1422	DO_2D_00_00
1423	ikbu = mbku(ji,jj)
1424	ikbv = mbkv(ji,jj)
1425	zu_i(ji,jj) = puu(ji,jj,ikbu,Kmm) - puu_b(ji,jj,Kmm)
1426	zv_i(ji,jj) = pvv(ji,jj,ikbv,Kmm) - pvv_b(ji,jj,Kmm)
1427	END_2D
1428	ELSE ! CENTRED integration: use BEFORE bottom baroclinic velocities
1429
1430	DO_2D_00_00
1431	ikbu = mbku(ji,jj)
1432	ikbv = mbkv(ji,jj)
1433	zu_i(ji,jj) = puu(ji,jj,ikbu,Kbb) - puu_b(ji,jj,Kbb)
1434	zv_i(ji,jj) = pvv(ji,jj,ikbv,Kbb) - pvv_b(ji,jj,Kbb)
1435	END_2D
1436	ENDIF
1437	!
1438	IF( ln_wd_il ) THEN ! W/D : use the "clipped" bottom friction !!gm explain WHY, please !
1439	zztmp = -1._wp / rDt_e
1440	DO_2D_00_00
1441	pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + zu_i(ji,jj) * wdrampu(ji,jj) * MAX( &
1442	& r1_hu(ji,jj,Kmm) * r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) , zztmp )
1443	pv_RHSi(ji,jj) = pv_RHSi(ji,jj) + zv_i(ji,jj) * wdrampv(ji,jj) * MAX( &
1444	& r1_hv(ji,jj,Kmm) * r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) , zztmp )
1445	END_2D
1446	ELSE ! use "unclipped" drag (even if explicit friction is used in 3D calculation)
1447
1448	DO_2D_00_00
1449	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)
1450	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)
1451	END_2D
1452	END IF
1453	!
1454	! !== TOP stress contribution from baroclinic velocities ==! (no W/D case)
1455	!
1456	IF( ln_isfcav ) THEN
1457	!
1458	IF( ln_bt_fw ) THEN ! FORWARD integration: use NOW top baroclinic velocity
1459
1460	DO_2D_00_00
1461	iktu = miku(ji,jj)
1462	iktv = mikv(ji,jj)
1463	zu_i(ji,jj) = puu(ji,jj,iktu,Kmm) - puu_b(ji,jj,Kmm)
1464	zv_i(ji,jj) = pvv(ji,jj,iktv,Kmm) - pvv_b(ji,jj,Kmm)
1465	END_2D
1466	ELSE ! CENTRED integration: use BEFORE top baroclinic velocity
1467
1468	DO_2D_00_00
1469	iktu = miku(ji,jj)
1470	iktv = mikv(ji,jj)
1471	zu_i(ji,jj) = puu(ji,jj,iktu,Kbb) - puu_b(ji,jj,Kbb)
1472	zv_i(ji,jj) = pvv(ji,jj,iktv,Kbb) - pvv_b(ji,jj,Kbb)
1473	END_2D
1474	ENDIF
1475	!
1476	! ! use "unclipped" top drag (even if explicit friction is used in 3D calculation)
1477
1478	DO_2D_00_00
1479	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)
1480	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)
1481	END_2D
1482	!
1483	ENDIF
1484	!
1485	END SUBROUTINE dyn_drg_init
1486
1487	SUBROUTINE ts_bck_interp( jn, ll_init, & ! <== in
1488	& za0, za1, za2, za3 ) ! ==> out
1489	!!----------------------------------------------------------------------
1490	INTEGER ,INTENT(in ) :: jn ! index of sub time step
1491	LOGICAL ,INTENT(in ) :: ll_init !
1492	REAL(wp),INTENT( out) :: za0, za1, za2, za3 ! Half-step back interpolation coefficient
1493	!
1494	REAL(wp) :: zepsilon, zgamma ! - -
1495	!!----------------------------------------------------------------------
1496	! ! set Half-step back interpolation coefficient
1497	IF ( jn==1 .AND. ll_init ) THEN !* Forward-backward
1498	za0 = 1._wp
1499	za1 = 0._wp
1500	za2 = 0._wp
1501	za3 = 0._wp
1502	ELSEIF( jn==2 .AND. ll_init ) THEN !* AB2-AM3 Coefficients; bet=0 ; gam=-1/6 ; eps=1/12
1503	za0 = 1.0833333333333_wp ! za0 = 1-gam-eps
1504	za1 =-0.1666666666666_wp ! za1 = gam
1505	za2 = 0.0833333333333_wp ! za2 = eps
1506	za3 = 0._wp
1507	ELSE !* AB3-AM4 Coefficients; bet=0.281105 ; eps=0.013 ; gam=0.0880
1508	IF( rn_bt_alpha == 0._wp ) THEN ! Time diffusion
1509	za0 = 0.614_wp ! za0 = 1/2 + gam + 2*eps
1510	za1 = 0.285_wp ! za1 = 1/2 - 2gam - 3eps
1511	za2 = 0.088_wp ! za2 = gam
1512	za3 = 0.013_wp ! za3 = eps
1513	ELSE ! no time diffusion
1514	zepsilon = 0.00976186_wp - 0.13451357_wp * rn_bt_alpha
1515	zgamma = 0.08344500_wp - 0.51358400_wp * rn_bt_alpha
1516	za0 = 0.5_wp + zgamma + 2._wp * rn_bt_alpha + 2._wp * zepsilon
1517	za1 = 1._wp - za0 - zgamma - zepsilon
1518	za2 = zgamma
1519	za3 = zepsilon
1520	ENDIF
1521	ENDIF
1522	END SUBROUTINE ts_bck_interp
1523
1524
1525	!!======================================================================
1526	END MODULE dynspg_ts

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

Download in other formats: