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

Context Navigation

source: NEMO/trunk/src/OCE/DYN/dynspg_ts.F90 @ 14086

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

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

Download in other formats: