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_12905_xios_restart/src/OCE/DYN – NEMO

Context Navigation

source: NEMO/branches/2020/dev_12905_xios_restart/src/OCE/DYN/dynspg_ts.F90 @ 12950

Last change on this file since 12950 was 12950, checked in by andmirek, 4 years ago
Ticket #2462: new XIOS restart read/write interfaces
Property svn:keywords set to `Id`
File size: 77.5 KB

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

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

Download in other formats: