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

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

dynspg_ts.F90 in NEMO/branches/2020/dev_r14116_HPC-04_mcastril_Mixed_Precision_implementation_final/src/OCE/DYN – NEMO

Context Navigation

source: NEMO/branches/2020/dev_r14116_HPC-04_mcastril_Mixed_Precision_implementation_final/src/OCE/DYN/dynspg_ts.F90 @ 14286

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

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

Download in other formats: