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_r12558_HPC-08_epico_Extra_Halo/src/OCE/DYN – NEMO

Context Navigation

source: NEMO/branches/2020/dev_r12558_HPC-08_epico_Extra_Halo/src/OCE/DYN/dynspg_ts.F90 @ 12738

Last change on this file since 12738 was 12738, checked in by smasson, 4 years ago
Extra_Halo: iom cleaning/update to work only with unknown, global or local (without halos) domains, see #2366
Property svn:keywords set to `Id`
File size: 78.6 KB

Line
1	MODULE dynspg_ts
2
3	!! Includes ROMS wd scheme with diagnostic outputs ; puu(:,:,:,Kmm) and puu(:,:,:,Krhs) updates are commented out !
4
5	!!======================================================================
6	!! * MODULE dynspg_ts *
7	!! Ocean dynamics: surface pressure gradient trend, split-explicit scheme
8	!!======================================================================
9	!! History : 1.0 ! 2004-12 (L. Bessieres, G. Madec) Original code
10	!! - ! 2005-11 (V. Garnier, G. Madec) optimization
11	!! - ! 2006-08 (S. Masson) distributed restart using iom
12	!! 2.0 ! 2007-07 (D. Storkey) calls to BDY routines
13	!! - ! 2008-01 (R. Benshila) change averaging method
14	!! 3.2 ! 2009-07 (R. Benshila, G. Madec) Complete revisit associated to vvl reactivation
15	!! 3.3 ! 2010-09 (D. Storkey, E. O'Dea) update for BDY for Shelf configurations
16	!! 3.3 ! 2011-03 (R. Benshila, R. Hordoir, P. Oddo) update calculation of ub_b
17	!! 3.5 ! 2013-07 (J. Chanut) Switch to Forward-backward time stepping
18	!! 3.6 ! 2013-11 (A. Coward) Update for z-tilde compatibility
19	!! 3.7 ! 2015-11 (J. Chanut) free surface simplification
20	!! - ! 2016-12 (G. Madec, E. Clementi) update for Stoke-Drift divergence
21	!! 4.0 ! 2017-05 (G. Madec) drag coef. defined at t-point (zdfdrg.F90)
22	!!---------------------------------------------------------------------
23
24	!!----------------------------------------------------------------------
25	!! dyn_spg_ts : compute surface pressure gradient trend using a time-splitting scheme
26	!! dyn_spg_ts_init: initialisation of the time-splitting scheme
27	!! ts_wgt : set time-splitting weights for temporal averaging (or not)
28	!! ts_rst : read/write time-splitting fields in restart file
29	!!----------------------------------------------------------------------
30	USE oce ! ocean dynamics and tracers
31	USE dom_oce ! ocean space and time domain
32	USE sbc_oce ! surface boundary condition: ocean
33	USE isf_oce ! ice shelf variable (fwfisf)
34	USE zdf_oce ! vertical physics: variables
35	USE zdfdrg ! vertical physics: top/bottom drag coef.
36	USE sbcapr ! surface boundary condition: atmospheric pressure
37	USE dynadv , ONLY: ln_dynadv_vec
38	USE dynvor ! vortivity scheme indicators
39	USE phycst ! physical constants
40	USE dynvor ! vorticity term
41	USE wet_dry ! wetting/drying flux limter
42	USE bdy_oce ! open boundary
43	USE bdyvol ! open boundary volume conservation
44	USE bdytides ! open boundary condition data
45	USE bdydyn2d ! open boundary conditions on barotropic variables
46	USE tide_mod !
47	USE sbcwave ! surface wave
48	#if defined key_agrif
49	USE agrif_oce_interp ! agrif
50	USE agrif_oce
51	#endif
52	#if defined key_asminc
53	USE asminc ! Assimilation increment
54	#endif
55	!
56	USE in_out_manager ! I/O manager
57	USE lib_mpp ! distributed memory computing library
58	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
59	USE prtctl ! Print control
60	USE iom ! IOM library
61	USE restart ! only for lrst_oce
62
63	USE iom ! to remove
64
65	IMPLICIT NONE
66	PRIVATE
67
68	PUBLIC dyn_spg_ts ! called by dyn_spg
69	PUBLIC dyn_spg_ts_init ! - - dyn_spg_init
70
71	!! Time filtered arrays at baroclinic time step:
72	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: un_adv , vn_adv !: Advection vel. at "now" barocl. step
73	!
74	INTEGER, SAVE :: icycle ! Number of barotropic sub-steps for each internal step nn_e <= 2.5 nn_e
75	REAL(wp),SAVE :: rDt_e ! Barotropic time step
76	!
77	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: wgtbtp1, wgtbtp2 ! 1st & 2nd weights used in time filtering of barotropic fields
78	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zwz ! ff_f/h at F points
79	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ftnw, ftne ! triad of coriolis parameter
80	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ftsw, ftse ! (only used with een vorticity scheme)
81
82	REAL(wp) :: r1_12 = 1._wp / 12._wp ! local ratios
83	REAL(wp) :: r1_8 = 0.125_wp !
84	REAL(wp) :: r1_4 = 0.25_wp !
85	REAL(wp) :: r1_2 = 0.5_wp !
86
87	!! * Substitutions
88	# include "do_loop_substitute.h90"
89	!!----------------------------------------------------------------------
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( ln_rstart .AND. ln_bt_fw .AND. (.NOT.l_1st_euler) ) THEN !* Read the restart file
893	CALL iom_get( numror, jpdom_auto, 'ub2_b' , ub2_b (:,:), cd_type = 'U', psgn = -1._wp, ldxios = lrxios )
894	CALL iom_get( numror, jpdom_auto, 'vb2_b' , vb2_b (:,:), cd_type = 'V', psgn = -1._wp, ldxios = lrxios )
895	CALL iom_get( numror, jpdom_auto, 'un_bf' , un_bf (:,:), cd_type = 'U', psgn = -1._wp, ldxios = lrxios )
896	CALL iom_get( numror, jpdom_auto, 'vn_bf' , vn_bf (:,:), cd_type = 'V', psgn = -1._wp, ldxios = lrxios )
897	IF( .NOT.ln_bt_av ) THEN
898	CALL iom_get( numror, jpdom_auto, 'sshbb_e' , sshbb_e(:,:), cd_type = 'T', psgn = 1._wp, ldxios = lrxios )
899	CALL iom_get( numror, jpdom_auto, 'ubb_e' , ubb_e(:,:), cd_type = 'U', psgn = -1._wp, ldxios = lrxios )
900	CALL iom_get( numror, jpdom_auto, 'vbb_e' , vbb_e(:,:), cd_type = 'V', psgn = -1._wp, ldxios = lrxios )
901	CALL iom_get( numror, jpdom_auto, 'sshb_e' , sshb_e(:,:), cd_type = 'T', psgn = 1._wp, ldxios = lrxios )
902	CALL iom_get( numror, jpdom_auto, 'ub_e' , ub_e(:,:), cd_type = 'U', psgn = -1._wp, ldxios = lrxios )
903	CALL iom_get( numror, jpdom_auto, 'vb_e' , vb_e(:,:), cd_type = 'V', psgn = -1._wp, ldxios = lrxios )
904	ENDIF
905	#if defined key_agrif
906	! Read time integrated fluxes
907	IF ( .NOT.Agrif_Root() ) THEN
908	CALL iom_get( numror, jpdom_auto, 'ub2_i_b' , ub2_i_b(:,:), cd_type = 'U', psgn = -1._wp, ldxios = lrxios )
909	CALL iom_get( numror, jpdom_auto, 'vb2_i_b' , vb2_i_b(:,:), cd_type = 'V', psgn = -1._wp, ldxios = lrxios )
910	ENDIF
911	#endif
912	ELSE !* Start from rest
913	IF(lwp) WRITE(numout,*)
914	IF(lwp) WRITE(numout,*) ' ==>>> start from rest: set barotropic values to 0'
915	ub2_b (:,:) = 0._wp ; vb2_b (:,:) = 0._wp ! used in the 1st interpol of agrif
916	un_adv(:,:) = 0._wp ; vn_adv(:,:) = 0._wp ! used in the 1st interpol of agrif
917	un_bf (:,:) = 0._wp ; vn_bf (:,:) = 0._wp ! used in the 1st update of agrif
918	#if defined key_agrif
919	IF ( .NOT.Agrif_Root() ) THEN
920	ub2_i_b(:,:) = 0._wp ; vb2_i_b(:,:) = 0._wp ! used in the 1st update of agrif
921	ENDIF
922	#endif
923	ENDIF
924	!
925	ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN ! Create restart file
926	! ! -------------------
927	IF(lwp) WRITE(numout,*) '---- ts_rst ----'
928	IF( lwxios ) CALL iom_swap( cwxios_context )
929	CALL iom_rstput( kt, nitrst, numrow, 'ub2_b' , ub2_b (:,:), ldxios = lwxios )
930	CALL iom_rstput( kt, nitrst, numrow, 'vb2_b' , vb2_b (:,:), ldxios = lwxios )
931	CALL iom_rstput( kt, nitrst, numrow, 'un_bf' , un_bf (:,:), ldxios = lwxios )
932	CALL iom_rstput( kt, nitrst, numrow, 'vn_bf' , vn_bf (:,:), ldxios = lwxios )
933	!
934	IF (.NOT.ln_bt_av) THEN
935	CALL iom_rstput( kt, nitrst, numrow, 'sshbb_e' , sshbb_e(:,:), ldxios = lwxios )
936	CALL iom_rstput( kt, nitrst, numrow, 'ubb_e' , ubb_e(:,:), ldxios = lwxios )
937	CALL iom_rstput( kt, nitrst, numrow, 'vbb_e' , vbb_e(:,:), ldxios = lwxios )
938	CALL iom_rstput( kt, nitrst, numrow, 'sshb_e' , sshb_e(:,:), ldxios = lwxios )
939	CALL iom_rstput( kt, nitrst, numrow, 'ub_e' , ub_e(:,:), ldxios = lwxios )
940	CALL iom_rstput( kt, nitrst, numrow, 'vb_e' , vb_e(:,:), ldxios = lwxios )
941	ENDIF
942	#if defined key_agrif
943	! Save time integrated fluxes
944	IF ( .NOT.Agrif_Root() ) THEN
945	CALL iom_rstput( kt, nitrst, numrow, 'ub2_i_b' , ub2_i_b(:,:), ldxios = lwxios )
946	CALL iom_rstput( kt, nitrst, numrow, 'vb2_i_b' , vb2_i_b(:,:), ldxios = lwxios )
947	ENDIF
948	#endif
949	IF( lwxios ) CALL iom_swap( cxios_context )
950	ENDIF
951	!
952	END SUBROUTINE ts_rst
953
954
955	SUBROUTINE dyn_spg_ts_init
956	!!---------------------------------------------------------------------
957	!! * ROUTINE dyn_spg_ts_init *
958	!!
959	!! ** Purpose : Set time splitting options
960	!!----------------------------------------------------------------------
961	INTEGER :: ji ,jj ! dummy loop indices
962	REAL(wp) :: zxr2, zyr2, zcmax ! local scalar
963	REAL(wp), DIMENSION(jpi,jpj) :: zcu
964	!!----------------------------------------------------------------------
965	!
966	! Max courant number for ext. grav. waves
967	!
968	DO_2D_11_11
969	zxr2 = r1_e1t(ji,jj) * r1_e1t(ji,jj)
970	zyr2 = r1_e2t(ji,jj) * r1_e2t(ji,jj)
971	zcu(ji,jj) = SQRT( grav * MAX(ht_0(ji,jj),0._wp) * (zxr2 + zyr2) )
972	END_2D
973	!
974	zcmax = MAXVAL( zcu(:,:) )
975	CALL mpp_max( 'dynspg_ts', zcmax )
976
977	! Estimate number of iterations to satisfy a max courant number= rn_bt_cmax
978	IF( ln_bt_auto ) nn_e = CEILING( rn_Dt / rn_bt_cmax * zcmax)
979
980	rDt_e = rn_Dt / REAL( nn_e , wp )
981	zcmax = zcmax * rDt_e
982	! Print results
983	IF(lwp) WRITE(numout,*)
984	IF(lwp) WRITE(numout,*) 'dyn_spg_ts_init : split-explicit free surface'
985	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~'
986	IF( ln_bt_auto ) THEN
987	IF(lwp) WRITE(numout,*) ' ln_ts_auto =.true. Automatically set nn_e '
988	IF(lwp) WRITE(numout,*) ' Max. courant number allowed: ', rn_bt_cmax
989	ELSE
990	IF(lwp) WRITE(numout,*) ' ln_ts_auto=.false.: Use nn_e in namelist nn_e = ', nn_e
991	ENDIF
992
993	IF(ln_bt_av) THEN
994	IF(lwp) WRITE(numout,*) ' ln_bt_av =.true. ==> Time averaging over nn_e time steps is on '
995	ELSE
996	IF(lwp) WRITE(numout,*) ' ln_bt_av =.false. => No time averaging of barotropic variables '
997	ENDIF
998	!
999	!
1000	IF(ln_bt_fw) THEN
1001	IF(lwp) WRITE(numout,*) ' ln_bt_fw=.true. => Forward integration of barotropic variables '
1002	ELSE
1003	IF(lwp) WRITE(numout,*) ' ln_bt_fw =.false.=> Centred integration of barotropic variables '
1004	ENDIF
1005	!
1006	#if defined key_agrif
1007	! Restrict the use of Agrif to the forward case only
1008	!!! IF( .NOT.ln_bt_fw .AND. .NOT.Agrif_Root() ) CALL ctl_stop( 'AGRIF not implemented if ln_bt_fw=.FALSE.' )
1009	#endif
1010	!
1011	IF(lwp) WRITE(numout,*) ' Time filter choice, nn_bt_flt: ', nn_bt_flt
1012	SELECT CASE ( nn_bt_flt )
1013	CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' Dirac'
1014	CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' Boxcar: width = nn_e'
1015	CASE( 2 ) ; IF(lwp) WRITE(numout,) ' Boxcar: width = 2nn_e'
1016	CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_bt_flt: should 0,1, or 2' )
1017	END SELECT
1018	!
1019	IF(lwp) WRITE(numout,*) ' '
1020	IF(lwp) WRITE(numout,*) ' nn_e = ', nn_e
1021	IF(lwp) WRITE(numout,*) ' Barotropic time step [s] is :', rDt_e
1022	IF(lwp) WRITE(numout,*) ' Maximum Courant number is :', zcmax
1023	!
1024	IF(lwp) WRITE(numout,*) ' Time diffusion parameter rn_bt_alpha: ', rn_bt_alpha
1025	IF ((ln_bt_av.AND.nn_bt_flt/=0).AND.(rn_bt_alpha>0._wp)) THEN
1026	CALL ctl_stop( 'dynspg_ts ERROR: if rn_bt_alpha > 0, remove temporal averaging' )
1027	ENDIF
1028	!
1029	IF( .NOT.ln_bt_av .AND. .NOT.ln_bt_fw ) THEN
1030	CALL ctl_stop( 'dynspg_ts ERROR: No time averaging => only forward integration is possible' )
1031	ENDIF
1032	IF( zcmax>0.9_wp ) THEN
1033	CALL ctl_stop( 'dynspg_ts ERROR: Maximum Courant number is greater than 0.9: Inc. nn_e !' )
1034	ENDIF
1035	!
1036	! ! Allocate time-splitting arrays
1037	IF( dyn_spg_ts_alloc() /= 0 ) CALL ctl_stop('STOP', 'dyn_spg_init: failed to allocate dynspg_ts arrays' )
1038	!
1039	! ! read restart when needed
1040	CALL ts_rst( nit000, 'READ' )
1041	!
1042	IF( lwxios ) THEN
1043	! define variables in restart file when writing with XIOS
1044	CALL iom_set_rstw_var_active('ub2_b')
1045	CALL iom_set_rstw_var_active('vb2_b')
1046	CALL iom_set_rstw_var_active('un_bf')
1047	CALL iom_set_rstw_var_active('vn_bf')
1048	!
1049	IF (.NOT.ln_bt_av) THEN
1050	CALL iom_set_rstw_var_active('sshbb_e')
1051	CALL iom_set_rstw_var_active('ubb_e')
1052	CALL iom_set_rstw_var_active('vbb_e')
1053	CALL iom_set_rstw_var_active('sshb_e')
1054	CALL iom_set_rstw_var_active('ub_e')
1055	CALL iom_set_rstw_var_active('vb_e')
1056	ENDIF
1057	#if defined key_agrif
1058	! Save time integrated fluxes
1059	IF ( .NOT.Agrif_Root() ) THEN
1060	CALL iom_set_rstw_var_active('ub2_i_b')
1061	CALL iom_set_rstw_var_active('vb2_i_b')
1062	ENDIF
1063	#endif
1064	ENDIF
1065	!
1066	END SUBROUTINE dyn_spg_ts_init
1067
1068
1069	SUBROUTINE dyn_cor_2D_init( Kmm )
1070	!!---------------------------------------------------------------------
1071	!! * ROUTINE dyn_cor_2D_init *
1072	!!
1073	!! ** Purpose : Set time splitting options
1074	!! Set arrays to remove/compute coriolis trend.
1075	!! Do it once during initialization if volume is fixed, else at each long time step.
1076	!! Note that these arrays are also used during barotropic loop. These are however frozen
1077	!! although they should be updated in the variable volume case. Not a big approximation.
1078	!! To remove this approximation, copy lines below inside barotropic loop
1079	!! and update depths at T-F points (ht and zhf resp.) at each barotropic time step
1080	!!
1081	!! Compute zwz = f / ( height of the water colomn )
1082	!!----------------------------------------------------------------------
1083	INTEGER, INTENT(in) :: Kmm ! Time index
1084	INTEGER :: ji ,jj, jk ! dummy loop indices
1085	REAL(wp) :: z1_ht
1086	REAL(wp), DIMENSION(jpi,jpj) :: zhf
1087	!!----------------------------------------------------------------------
1088	!
1089	SELECT CASE( nvor_scheme )
1090	CASE( np_EEN ) != EEN scheme using e3f (energy & enstrophy scheme)
1091	SELECT CASE( nn_een_e3f ) !* ff_f/e3 at F-point
1092	CASE ( 0 ) ! original formulation (masked averaging of e3t divided by 4)
1093	DO_2D_10_10
1094	zwz(ji,jj) = ( ht(ji ,jj+1) + ht(ji+1,jj+1) + &
1095	& ht(ji ,jj ) + ht(ji+1,jj ) ) * 0.25_wp
1096	IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj)
1097	END_2D
1098	CASE ( 1 ) ! new formulation (masked averaging of e3t divided by the sum of mask)
1099	DO_2D_10_10
1100	zwz(ji,jj) = ( ht (ji ,jj+1) + ht (ji+1,jj+1) &
1101	& + ht (ji ,jj ) + ht (ji+1,jj ) ) &
1102	& / ( MAX( 1._wp, ssmask(ji ,jj+1) + ssmask(ji+1,jj+1) &
1103	& + ssmask(ji ,jj ) + ssmask(ji+1,jj ) ) )
1104	IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj)
1105	END_2D
1106	END SELECT
1107	CALL lbc_lnk( 'dynspg_ts', zwz, 'F', 1._wp )
1108	!
1109	ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp
1110	DO_2D_01_01
1111	ftne(ji,jj) = zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1)
1112	ftnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj ) + zwz(ji ,jj )
1113	ftse(ji,jj) = zwz(ji ,jj ) + zwz(ji ,jj-1) + zwz(ji-1,jj-1)
1114	ftsw(ji,jj) = zwz(ji ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj )
1115	END_2D
1116	!
1117	CASE( np_EET ) != EEN scheme using e3t (energy conserving scheme)
1118	ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp
1119	DO_2D_01_01
1120	z1_ht = ssmask(ji,jj) / ( ht(ji,jj) + 1._wp - ssmask(ji,jj) )
1121	ftne(ji,jj) = ( ff_f(ji-1,jj ) + ff_f(ji ,jj ) + ff_f(ji ,jj-1) ) * z1_ht
1122	ftnw(ji,jj) = ( ff_f(ji-1,jj-1) + ff_f(ji-1,jj ) + ff_f(ji ,jj ) ) * z1_ht
1123	ftse(ji,jj) = ( ff_f(ji ,jj ) + ff_f(ji ,jj-1) + ff_f(ji-1,jj-1) ) * z1_ht
1124	ftsw(ji,jj) = ( ff_f(ji ,jj-1) + ff_f(ji-1,jj-1) + ff_f(ji-1,jj ) ) * z1_ht
1125	END_2D
1126	!
1127	CASE( np_ENE, np_ENS , np_MIX ) != all other schemes (ENE, ENS, MIX) except ENT !
1128	!
1129	zwz(:,:) = 0._wp
1130	zhf(:,:) = 0._wp
1131
1132	!!gm assume 0 in both cases (which is almost surely WRONG ! ) as hvatf has been removed
1133	!!gm A priori a better value should be something like :
1134	!!gm zhf(i,j) = masked sum of ht(i,j) , ht(i+1,j) , ht(i,j+1) , (i+1,j+1)
1135	!!gm divided by the sum of the corresponding mask
1136	!!gm
1137	!!
1138	IF( .NOT.ln_sco ) THEN
1139
1140	!!gm agree the JC comment : this should be done in a much clear way
1141
1142	! JC: It not clear yet what should be the depth at f-points over land in z-coordinate case
1143	! Set it to zero for the time being
1144	! IF( rn_hmin < 0._wp ) THEN ; jk = - INT( rn_hmin ) ! from a nb of level
1145	! ELSE ; jk = MINLOC( gdepw_0, mask = gdepw_0 > rn_hmin, dim = 1 ) ! from a depth
1146	! ENDIF
1147	! zhf(:,:) = gdepw_0(:,:,jk+1)
1148	!
1149	ELSE
1150	!
1151	!zhf(:,:) = hbatf(:,:)
1152	DO_2D_10_10
1153	zhf(ji,jj) = ( ht_0 (ji,jj ) + ht_0 (ji+1,jj ) &
1154	& + ht_0 (ji,jj+1) + ht_0 (ji+1,jj+1) ) &
1155	& / MAX( ssmask(ji,jj ) + ssmask(ji+1,jj ) &
1156	& + ssmask(ji,jj+1) + ssmask(ji+1,jj+1) , 1._wp )
1157	END_2D
1158	ENDIF
1159	!
1160	DO jj = 1, jpjm1
1161	zhf(:,jj) = zhf(:,jj) * (1._wp- umask(:,jj,1) * umask(:,jj+1,1))
1162	END DO
1163	!
1164	DO jk = 1, jpkm1
1165	DO jj = 1, jpjm1
1166	zhf(:,jj) = zhf(:,jj) + e3f(:,jj,jk) * umask(:,jj,jk) * umask(:,jj+1,jk)
1167	END DO
1168	END DO
1169	CALL lbc_lnk( 'dynspg_ts', zhf, 'F', 1._wp )
1170	! JC: TBC. hf should be greater than 0
1171	DO_2D_11_11
1172	IF( zhf(ji,jj) /= 0._wp ) zwz(ji,jj) = 1._wp / zhf(ji,jj)
1173	END_2D
1174	zwz(:,:) = ff_f(:,:) * zwz(:,:)
1175	END SELECT
1176
1177	END SUBROUTINE dyn_cor_2d_init
1178
1179
1180
1181	SUBROUTINE dyn_cor_2d( phu, phv, punb, pvnb, zhU, zhV, zu_trd, zv_trd )
1182	!!---------------------------------------------------------------------
1183	!! * ROUTINE dyn_cor_2d *
1184	!!
1185	!! ** Purpose : Compute u and v coriolis trends
1186	!!----------------------------------------------------------------------
1187	INTEGER :: ji ,jj ! dummy loop indices
1188	REAL(wp) :: zx1, zx2, zy1, zy2, z1_hu, z1_hv ! - -
1189	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: phu, phv, punb, pvnb, zhU, zhV
1190	REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: zu_trd, zv_trd
1191	!!----------------------------------------------------------------------
1192	SELECT CASE( nvor_scheme )
1193	CASE( np_ENT ) ! enstrophy conserving scheme (f-point)
1194	DO_2D_00_00
1195	z1_hu = ssumask(ji,jj) / ( phu(ji,jj) + 1._wp - ssumask(ji,jj) )
1196	z1_hv = ssvmask(ji,jj) / ( phv(ji,jj) + 1._wp - ssvmask(ji,jj) )
1197	zu_trd(ji,jj) = + r1_4 * r1_e1e2u(ji,jj) * z1_hu &
1198	& * ( e1e2t(ji+1,jj)ht(ji+1,jj)ff_t(ji+1,jj) * ( pvnb(ji+1,jj) + pvnb(ji+1,jj-1) ) &
1199	& + e1e2t(ji ,jj)ht(ji ,jj)ff_t(ji ,jj) * ( pvnb(ji ,jj) + pvnb(ji ,jj-1) ) )
1200	!
1201	zv_trd(ji,jj) = - r1_4 * r1_e1e2v(ji,jj) * z1_hv &
1202	& * ( e1e2t(ji,jj+1)ht(ji,jj+1)ff_t(ji,jj+1) * ( punb(ji,jj+1) + punb(ji-1,jj+1) ) &
1203	& + e1e2t(ji,jj )ht(ji,jj )ff_t(ji,jj ) * ( punb(ji,jj ) + punb(ji-1,jj ) ) )
1204	END_2D
1205	!
1206	CASE( np_ENE , np_MIX ) ! energy conserving scheme (t-point) ENE or MIX
1207	DO_2D_00_00
1208	zy1 = ( zhV(ji,jj-1) + zhV(ji+1,jj-1) ) * r1_e1u(ji,jj)
1209	zy2 = ( zhV(ji,jj ) + zhV(ji+1,jj ) ) * r1_e1u(ji,jj)
1210	zx1 = ( zhU(ji-1,jj) + zhU(ji-1,jj+1) ) * r1_e2v(ji,jj)
1211	zx2 = ( zhU(ji ,jj) + zhU(ji ,jj+1) ) * r1_e2v(ji,jj)
1212	! energy conserving formulation for planetary vorticity term
1213	zu_trd(ji,jj) = r1_4 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
1214	zv_trd(ji,jj) = - r1_4 * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 )
1215	END_2D
1216	!
1217	CASE( np_ENS ) ! enstrophy conserving scheme (f-point)
1218	DO_2D_00_00
1219	zy1 = r1_8 * ( zhV(ji ,jj-1) + zhV(ji+1,jj-1) &
1220	& + zhV(ji ,jj ) + zhV(ji+1,jj ) ) * r1_e1u(ji,jj)
1221	zx1 = - r1_8 * ( zhU(ji-1,jj ) + zhU(ji-1,jj+1) &
1222	& + zhU(ji ,jj ) + zhU(ji ,jj+1) ) * r1_e2v(ji,jj)
1223	zu_trd(ji,jj) = zy1 * ( zwz(ji ,jj-1) + zwz(ji,jj) )
1224	zv_trd(ji,jj) = zx1 * ( zwz(ji-1,jj ) + zwz(ji,jj) )
1225	END_2D
1226	!
1227	CASE( np_EET , np_EEN ) ! energy & enstrophy scheme (using e3t or e3f)
1228	DO_2D_00_00
1229	zu_trd(ji,jj) = + r1_12 * r1_e1u(ji,jj) * ( ftne(ji,jj ) * zhV(ji ,jj ) &
1230	& + ftnw(ji+1,jj) * zhV(ji+1,jj ) &
1231	& + ftse(ji,jj ) * zhV(ji ,jj-1) &
1232	& + ftsw(ji+1,jj) * zhV(ji+1,jj-1) )
1233	zv_trd(ji,jj) = - r1_12 * r1_e2v(ji,jj) * ( ftsw(ji,jj+1) * zhU(ji-1,jj+1) &
1234	& + ftse(ji,jj+1) * zhU(ji ,jj+1) &
1235	& + ftnw(ji,jj ) * zhU(ji-1,jj ) &
1236	& + ftne(ji,jj ) * zhU(ji ,jj ) )
1237	END_2D
1238	!
1239	END SELECT
1240	!
1241	END SUBROUTINE dyn_cor_2D
1242
1243
1244	SUBROUTINE wad_tmsk( pssh, ptmsk )
1245	!!----------------------------------------------------------------------
1246	!! * ROUTINE wad_lmt *
1247	!!
1248	!! ** Purpose : set wetting & drying mask at tracer points
1249	!! for the current barotropic sub-step
1250	!!
1251	!! ** Method : ???
1252	!!
1253	!! ** Action : ptmsk : wetting & drying t-mask
1254	!!----------------------------------------------------------------------
1255	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pssh !
1256	REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: ptmsk !
1257	!
1258	INTEGER :: ji, jj ! dummy loop indices
1259	!!----------------------------------------------------------------------
1260	!
1261	IF( ln_wd_dl_rmp ) THEN
1262	DO_2D_11_11
1263	IF ( pssh(ji,jj) + ht_0(ji,jj) > 2._wp * rn_wdmin1 ) THEN
1264	! IF ( pssh(ji,jj) + ht_0(ji,jj) > rn_wdmin2 ) THEN
1265	ptmsk(ji,jj) = 1._wp
1266	ELSEIF( pssh(ji,jj) + ht_0(ji,jj) > rn_wdmin1 ) THEN
1267	ptmsk(ji,jj) = TANH( 50._wp( ( pssh(ji,jj) + ht_0(ji,jj) - rn_wdmin1 )r_rn_wdmin1) )
1268	ELSE
1269	ptmsk(ji,jj) = 0._wp
1270	ENDIF
1271	END_2D
1272	ELSE
1273	DO_2D_11_11
1274	IF ( pssh(ji,jj) + ht_0(ji,jj) > rn_wdmin1 ) THEN ; ptmsk(ji,jj) = 1._wp
1275	ELSE ; ptmsk(ji,jj) = 0._wp
1276	ENDIF
1277	END_2D
1278	ENDIF
1279	!
1280	END SUBROUTINE wad_tmsk
1281
1282
1283	SUBROUTINE wad_Umsk( pTmsk, phU, phV, pu, pv, pUmsk, pVmsk )
1284	!!----------------------------------------------------------------------
1285	!! * ROUTINE wad_lmt *
1286	!!
1287	!! ** Purpose : set wetting & drying mask at tracer points
1288	!! for the current barotropic sub-step
1289	!!
1290	!! ** Method : ???
1291	!!
1292	!! ** Action : ptmsk : wetting & drying t-mask
1293	!!----------------------------------------------------------------------
1294	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pTmsk ! W & D t-mask
1295	REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: phU, phV, pu, pv ! ocean velocities and transports
1296	REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pUmsk, pVmsk ! W & D u- and v-mask
1297	!
1298	INTEGER :: ji, jj ! dummy loop indices
1299	!!----------------------------------------------------------------------
1300	!
1301	DO_2D_11_10
1302	IF ( phU(ji,jj) > 0._wp ) THEN ; pUmsk(ji,jj) = pTmsk(ji ,jj)
1303	ELSE ; pUmsk(ji,jj) = pTmsk(ji+1,jj)
1304	ENDIF
1305	phU(ji,jj) = pUmsk(ji,jj)*phU(ji,jj)
1306	pu (ji,jj) = pUmsk(ji,jj)*pu (ji,jj)
1307	END_2D
1308	!
1309	DO_2D_10_11
1310	IF ( phV(ji,jj) > 0._wp ) THEN ; pVmsk(ji,jj) = pTmsk(ji,jj )
1311	ELSE ; pVmsk(ji,jj) = pTmsk(ji,jj+1)
1312	ENDIF
1313	phV(ji,jj) = pVmsk(ji,jj)*phV(ji,jj)
1314	pv (ji,jj) = pVmsk(ji,jj)*pv (ji,jj)
1315	END_2D
1316	!
1317	END SUBROUTINE wad_Umsk
1318
1319
1320	SUBROUTINE wad_spg( pshn, zcpx, zcpy )
1321	!!---------------------------------------------------------------------
1322	!! * ROUTINE wad_sp *
1323	!!
1324	!! ** Purpose :
1325	!!----------------------------------------------------------------------
1326	INTEGER :: ji ,jj ! dummy loop indices
1327	LOGICAL :: ll_tmp1, ll_tmp2
1328	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pshn
1329	REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: zcpx, zcpy
1330	!!----------------------------------------------------------------------
1331	DO_2D_00_00
1332	ll_tmp1 = MIN( pshn(ji,jj) , pshn(ji+1,jj) ) > &
1333	& MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) .AND. &
1334	& MAX( pshn(ji,jj) + ht_0(ji,jj) , pshn(ji+1,jj) + ht_0(ji+1,jj) ) &
1335	& > rn_wdmin1 + rn_wdmin2
1336	ll_tmp2 = ( ABS( pshn(ji+1,jj) - pshn(ji ,jj)) > 1.E-12 ).AND.( &
1337	& MAX( pshn(ji,jj) , pshn(ji+1,jj) ) > &
1338	& MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 )
1339	IF(ll_tmp1) THEN
1340	zcpx(ji,jj) = 1.0_wp
1341	ELSEIF(ll_tmp2) THEN
1342	! no worries about pshn(ji+1,jj) - pshn(ji ,jj) = 0, it won't happen ! here
1343	zcpx(ji,jj) = ABS( (pshn(ji+1,jj) + ht_0(ji+1,jj) - pshn(ji,jj) - ht_0(ji,jj)) &
1344	& / (pshn(ji+1,jj) - pshn(ji ,jj)) )
1345	zcpx(ji,jj) = max(min( zcpx(ji,jj) , 1.0_wp),0.0_wp)
1346	ELSE
1347	zcpx(ji,jj) = 0._wp
1348	ENDIF
1349	!
1350	ll_tmp1 = MIN( pshn(ji,jj) , pshn(ji,jj+1) ) > &
1351	& MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) .AND. &
1352	& MAX( pshn(ji,jj) + ht_0(ji,jj) , pshn(ji,jj+1) + ht_0(ji,jj+1) ) &
1353	& > rn_wdmin1 + rn_wdmin2
1354	ll_tmp2 = ( ABS( pshn(ji,jj) - pshn(ji,jj+1)) > 1.E-12 ).AND.( &
1355	& MAX( pshn(ji,jj) , pshn(ji,jj+1) ) > &
1356	& MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 )
1357
1358	IF(ll_tmp1) THEN
1359	zcpy(ji,jj) = 1.0_wp
1360	ELSE IF(ll_tmp2) THEN
1361	! no worries about pshn(ji,jj+1) - pshn(ji,jj ) = 0, it won't happen ! here
1362	zcpy(ji,jj) = ABS( (pshn(ji,jj+1) + ht_0(ji,jj+1) - pshn(ji,jj) - ht_0(ji,jj)) &
1363	& / (pshn(ji,jj+1) - pshn(ji,jj )) )
1364	zcpy(ji,jj) = MAX( 0._wp , MIN( zcpy(ji,jj) , 1.0_wp ) )
1365	ELSE
1366	zcpy(ji,jj) = 0._wp
1367	ENDIF
1368	END_2D
1369
1370	END SUBROUTINE wad_spg
1371
1372
1373
1374	SUBROUTINE dyn_drg_init( Kbb, Kmm, puu, pvv, puu_b ,pvv_b, pu_RHSi, pv_RHSi, pCdU_u, pCdU_v )
1375	!!----------------------------------------------------------------------
1376	!! * ROUTINE dyn_drg_init *
1377	!!
1378	!! ** Purpose : - add the baroclinic top/bottom drag contribution to
1379	!! the baroclinic part of the barotropic RHS
1380	!! - compute the barotropic drag coefficients
1381	!!
1382	!! ** Method : computation done over the INNER domain only
1383	!!----------------------------------------------------------------------
1384	INTEGER , INTENT(in ) :: Kbb, Kmm ! ocean time level indices
1385	REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(in ) :: puu, pvv ! ocean velocities and RHS of momentum equation
1386	REAL(wp), DIMENSION(jpi,jpj,jpt) , INTENT(in ) :: puu_b, pvv_b ! barotropic velocities at main time levels
1387	REAL(wp), DIMENSION(jpi,jpj) , INTENT(inout) :: pu_RHSi, pv_RHSi ! baroclinic part of the barotropic RHS
1388	REAL(wp), DIMENSION(jpi,jpj) , INTENT( out) :: pCdU_u , pCdU_v ! barotropic drag coefficients
1389	!
1390	INTEGER :: ji, jj ! dummy loop indices
1391	INTEGER :: ikbu, ikbv, iktu, iktv
1392	REAL(wp) :: zztmp
1393	REAL(wp), DIMENSION(jpi,jpj) :: zu_i, zv_i
1394	!!----------------------------------------------------------------------
1395	!
1396	! !== Set the barotropic drag coef. ==!
1397	!
1398	IF( ln_isfcav ) THEN ! top+bottom friction (ocean cavities)
1399
1400	DO_2D_00_00
1401	pCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) + rCdU_top(ji+1,jj)+rCdU_top(ji,jj) )
1402	pCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) + rCdU_top(ji,jj+1)+rCdU_top(ji,jj) )
1403	END_2D
1404	ELSE ! bottom friction only
1405	DO_2D_00_00
1406	pCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) )
1407	pCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) )
1408	END_2D
1409	ENDIF
1410	!
1411	! !== BOTTOM stress contribution from baroclinic velocities ==!
1412	!
1413	IF( ln_bt_fw ) THEN ! FORWARD integration: use NOW bottom baroclinic velocities
1414
1415	DO_2D_00_00
1416	ikbu = mbku(ji,jj)
1417	ikbv = mbkv(ji,jj)
1418	zu_i(ji,jj) = puu(ji,jj,ikbu,Kmm) - puu_b(ji,jj,Kmm)
1419	zv_i(ji,jj) = pvv(ji,jj,ikbv,Kmm) - pvv_b(ji,jj,Kmm)
1420	END_2D
1421	ELSE ! CENTRED integration: use BEFORE bottom baroclinic velocities
1422
1423	DO_2D_00_00
1424	ikbu = mbku(ji,jj)
1425	ikbv = mbkv(ji,jj)
1426	zu_i(ji,jj) = puu(ji,jj,ikbu,Kbb) - puu_b(ji,jj,Kbb)
1427	zv_i(ji,jj) = pvv(ji,jj,ikbv,Kbb) - pvv_b(ji,jj,Kbb)
1428	END_2D
1429	ENDIF
1430	!
1431	IF( ln_wd_il ) THEN ! W/D : use the "clipped" bottom friction !!gm explain WHY, please !
1432	zztmp = -1._wp / rDt_e
1433	DO_2D_00_00
1434	pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + zu_i(ji,jj) * wdrampu(ji,jj) * MAX( &
1435	& r1_hu(ji,jj,Kmm) * r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) , zztmp )
1436	pv_RHSi(ji,jj) = pv_RHSi(ji,jj) + zv_i(ji,jj) * wdrampv(ji,jj) * MAX( &
1437	& r1_hv(ji,jj,Kmm) * r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) , zztmp )
1438	END_2D
1439	ELSE ! use "unclipped" drag (even if explicit friction is used in 3D calculation)
1440
1441	DO_2D_00_00
1442	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)
1443	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)
1444	END_2D
1445	END IF
1446	!
1447	! !== TOP stress contribution from baroclinic velocities ==! (no W/D case)
1448	!
1449	IF( ln_isfcav ) THEN
1450	!
1451	IF( ln_bt_fw ) THEN ! FORWARD integration: use NOW top baroclinic velocity
1452
1453	DO_2D_00_00
1454	iktu = miku(ji,jj)
1455	iktv = mikv(ji,jj)
1456	zu_i(ji,jj) = puu(ji,jj,iktu,Kmm) - puu_b(ji,jj,Kmm)
1457	zv_i(ji,jj) = pvv(ji,jj,iktv,Kmm) - pvv_b(ji,jj,Kmm)
1458	END_2D
1459	ELSE ! CENTRED integration: use BEFORE top baroclinic velocity
1460
1461	DO_2D_00_00
1462	iktu = miku(ji,jj)
1463	iktv = mikv(ji,jj)
1464	zu_i(ji,jj) = puu(ji,jj,iktu,Kbb) - puu_b(ji,jj,Kbb)
1465	zv_i(ji,jj) = pvv(ji,jj,iktv,Kbb) - pvv_b(ji,jj,Kbb)
1466	END_2D
1467	ENDIF
1468	!
1469	! ! use "unclipped" top drag (even if explicit friction is used in 3D calculation)
1470
1471	DO_2D_00_00
1472	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)
1473	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)
1474	END_2D
1475	!
1476	ENDIF
1477	!
1478	END SUBROUTINE dyn_drg_init
1479
1480	SUBROUTINE ts_bck_interp( jn, ll_init, & ! <== in
1481	& za0, za1, za2, za3 ) ! ==> out
1482	!!----------------------------------------------------------------------
1483	INTEGER ,INTENT(in ) :: jn ! index of sub time step
1484	LOGICAL ,INTENT(in ) :: ll_init !
1485	REAL(wp),INTENT( out) :: za0, za1, za2, za3 ! Half-step back interpolation coefficient
1486	!
1487	REAL(wp) :: zepsilon, zgamma ! - -
1488	!!----------------------------------------------------------------------
1489	! ! set Half-step back interpolation coefficient
1490	IF ( jn==1 .AND. ll_init ) THEN !* Forward-backward
1491	za0 = 1._wp
1492	za1 = 0._wp
1493	za2 = 0._wp
1494	za3 = 0._wp
1495	ELSEIF( jn==2 .AND. ll_init ) THEN !* AB2-AM3 Coefficients; bet=0 ; gam=-1/6 ; eps=1/12
1496	za0 = 1.0833333333333_wp ! za0 = 1-gam-eps
1497	za1 =-0.1666666666666_wp ! za1 = gam
1498	za2 = 0.0833333333333_wp ! za2 = eps
1499	za3 = 0._wp
1500	ELSE !* AB3-AM4 Coefficients; bet=0.281105 ; eps=0.013 ; gam=0.0880
1501	IF( rn_bt_alpha == 0._wp ) THEN ! Time diffusion
1502	za0 = 0.614_wp ! za0 = 1/2 + gam + 2*eps
1503	za1 = 0.285_wp ! za1 = 1/2 - 2gam - 3eps
1504	za2 = 0.088_wp ! za2 = gam
1505	za3 = 0.013_wp ! za3 = eps
1506	ELSE ! no time diffusion
1507	zepsilon = 0.00976186_wp - 0.13451357_wp * rn_bt_alpha
1508	zgamma = 0.08344500_wp - 0.51358400_wp * rn_bt_alpha
1509	za0 = 0.5_wp + zgamma + 2._wp * rn_bt_alpha + 2._wp * zepsilon
1510	za1 = 1._wp - za0 - zgamma - zepsilon
1511	za2 = zgamma
1512	za3 = zepsilon
1513	ENDIF
1514	ENDIF
1515	END SUBROUTINE ts_bck_interp
1516
1517
1518	!!======================================================================
1519	END MODULE dynspg_ts

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

Download in other formats: