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 @ 14433

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

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

Download in other formats: