Context Navigation

dynspg_ts.F90 @ 15487

Last change on this file since 15487 was 15487, checked in by dford, 3 years ago
Remove some comments and add some spaces.
File size: 80.6 KB

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

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

New URL for NEMO forge! http://forge.nemo-ocean.eu

Context Navigation

source: NEMO/branches/UKMO/NEMO_4.0.4_generic_obs/src/OCE/DYN/dynspg_ts.F90 @ 15487

Download in other formats: