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

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

dynspg_ts.F90 in NEMO/trunk/src/OCE/DYN – NEMO

Context Navigation

source: NEMO/trunk/src/OCE/DYN/dynspg_ts.F90 @ 13286

Last change on this file since 13286 was 13286, checked in by smasson, 4 years ago
trunk: merge extra halos branch in trunk, see #2366
Property svn:keywords set to `Id`
File size: 78.9 KB

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

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

Download in other formats: