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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

dynspg_ts.F90 in NEMO/branches/2020/temporary_r4_trunk/src/OCE/DYN – NEMO

Context Navigation

source: NEMO/branches/2020/temporary_r4_trunk/src/OCE/DYN/dynspg_ts.F90 @ 13470

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

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

Download in other formats: