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

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

dynspg_ts.F90 in NEMO/trunk/src/OCE/DYN – NEMO

Context Navigation

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

Last change on this file since 15033 was 15018, checked in by smasson, 3 years ago
trunk: supress jpim1/jpjm1 in dynspg_ts, #2670
Property svn:keywords set to `Id`
File size: 77.1 KB

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

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

Download in other formats: