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

Last change on this file since 14207 was 14207, checked in by smasson, 4 years ago
trunk: minor typo in dynspg_ts...
Property svn:keywords set to `Id`
File size: 76.6 KB

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

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

Download in other formats: