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_r13327_KERNEL-06_2_techene_e3/src/OCE/DYN – NEMO

Context Navigation

source: NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/DYN/dynspg_ts.F90 @ 14018

Last change on this file since 14018 was 14018, checked in by techene, 3 years ago
#2385 branch updated with trunk 13970
Property svn:keywords set to `Id`
File size: 77.7 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( 0, 0, 0, 0 ) ! SPG with the application of W/D gravity filters
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( 0, 0, 0, 0 )
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( 0, 0, 0, 0 ) ! Remove coriolis term (and possibly spg) from barotropic trend
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( 0, 0, 0, 0 )
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( 0, 0, 0, 0 )
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 wind forcing =!
309	! ! ------------------ !
310	IF( ln_bt_fw ) THEN
311	DO_2D( 0, 0, 0, 0 )
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( 0, 0, 0, 0 )
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_0(:,:)
389	zhvp2_e(:,:) = hv_0(:,:)
390	zhtp2_e(:,:) = ht_0(:,:)
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	#if defined key_qcoTest_FluxForm
478	! ! 'key_qcoTest_FluxForm' : simple ssh average
479	DO_2D( 1, 1, 1, 0 ) ! not jpi-column
480	zhup2_e(ji,jj) = hu_0(ji,jj) + r1_2 * ( zsshp2_e(ji,jj) + zsshp2_e(ji+1,jj ) ) * ssumask(ji,jj)
481	END_2D
482	DO_2D( 1, 0, 1, 1 )
483	zhvp2_e(ji,jj) = hv_0(ji,jj) + r1_2 * ( zsshp2_e(ji,jj) + zsshp2_e(ji ,jj+1) ) * ssvmask(ji,jj)
484	END_2D
485	#else
486	! ! no 'key_qcoTest_FluxForm' : surface weighted ssh average
487	DO_2D( 1, 1, 1, 0 ) ! not jpi-column
488	zhup2_e(ji,jj) = hu_0(ji,jj) + r1_2 * r1_e1e2u(ji,jj) &
489	& * ( e1e2t(ji ,jj) * zsshp2_e(ji ,jj) &
490	& + e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj) ) * ssumask(ji,jj)
491	END_2D
492	DO_2D( 1, 0, 1, 1 ) ! not jpj-row
493	zhvp2_e(ji,jj) = hv_0(ji,jj) + r1_2 * r1_e1e2v(ji,jj) &
494	& * ( e1e2t(ji,jj ) * zsshp2_e(ji,jj ) &
495	& + e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1) ) * ssvmask(ji,jj)
496	END_2D
497	#endif
498	!
499	ENDIF
500	!
501	! !== after SSH ==! (jn+1)
502	!
503	! ! update (ua_e,va_e) to enforce volume conservation at open boundaries
504	! ! values of zhup2_e and zhvp2_e on the halo are not needed in bdy_vol2d
505	IF( ln_bdy .AND. ln_vol ) CALL bdy_vol2d( kt, jn, ua_e, va_e, zhup2_e, zhvp2_e )
506	!
507	! ! resulting flux at mid-step (not over the full domain)
508	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
509	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
510	!
511	#if defined key_agrif
512	! Set fluxes during predictor step to ensure volume conservation
513	IF( .NOT.Agrif_Root() .AND. ln_bt_fw ) CALL agrif_dyn_ts_flux( jn, zhU, zhV )
514	#endif
515	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
516
517	IF( ln_wd_dl ) THEN ! un_e and vn_e are set to zero at faces where
518	! ! the direction of the flow is from dry cells
519	CALL wad_Umsk( ztwdmask, zhU, zhV, un_e, vn_e, zuwdmask, zvwdmask ) ! not jpi colomn for U, not jpj row for V
520	!
521	ENDIF
522	!
523	!
524	! Compute Sea Level at step jit+1
525	!-- m+1 m m+1/2 --!
526	!-- ssh = ssh - delta_t' * [ frc + div( flux ) ] --!
527	!-------------------------------------------------------------------------!
528	DO_2D( 0, 0, 0, 0 )
529	zhdiv = ( zhU(ji,jj) - zhU(ji-1,jj) + zhV(ji,jj) - zhV(ji,jj-1) ) * r1_e1e2t(ji,jj)
530	ssha_e(ji,jj) = ( sshn_e(ji,jj) - rDt_e * ( zssh_frc(ji,jj) + zhdiv ) ) * ssmask(ji,jj)
531	END_2D
532	!
533	CALL lbc_lnk_multi( 'dynspg_ts', ssha_e, 'T', 1._wp, zhU, 'U', -1._wp, zhV, 'V', -1._wp )
534	!
535	! Duplicate sea level across open boundaries (this is only cosmetic if linssh=T)
536	IF( ln_bdy ) CALL bdy_ssh( ssha_e )
537	#if defined key_agrif
538	IF( .NOT.Agrif_Root() ) CALL agrif_ssh_ts( jn )
539	#endif
540	!
541	! ! Sum over sub-time-steps to compute advective velocities
542	za2 = wgtbtp2(jn) ! zhU, zhV hold fluxes extrapolated at jn+0.5
543	un_adv(:,:) = un_adv(:,:) + za2 * zhU(:,:) * r1_e2u(:,:)
544	vn_adv(:,:) = vn_adv(:,:) + za2 * zhV(:,:) * r1_e1v(:,:)
545	! sum over sub-time-steps to decide which baroclinic velocities to set to zero (zuwdav2 is only used when ln_wd_dl_bc=True)
546	IF ( ln_wd_dl_bc ) THEN
547	zuwdav2(1:jpim1,1:jpj ) = zuwdav2(1:jpim1,1:jpj ) + za2 * zuwdmask(1:jpim1,1:jpj ) ! not jpi-column
548	zvwdav2(1:jpi ,1:jpjm1) = zvwdav2(1:jpi ,1:jpjm1) + za2 * zvwdmask(1:jpi ,1:jpjm1) ! not jpj-row
549	END IF
550	!
551	!
552	! Sea Surface Height at u-,v-points (vvl case only)
553	IF( .NOT.ln_linssh ) THEN
554	#if defined key_qcoTest_FluxForm
555	! ! 'key_qcoTest_FluxForm' : simple ssh average
556	DO_2D( 1, 1, 1, 0 )
557	zsshu_a(ji,jj) = r1_2 * ( ssha_e(ji,jj) + ssha_e(ji+1,jj ) ) * ssumask(ji,jj)
558	END_2D
559	DO_2D( 1, 0, 1, 1 )
560	zsshv_a(ji,jj) = r1_2 * ( ssha_e(ji,jj) + ssha_e(ji ,jj+1) ) * ssvmask(ji,jj)
561	END_2D
562	#else
563	DO_2D( 0, 0, 0, 0 )
564	zsshu_a(ji,jj) = r1_2 * r1_e1e2u(ji,jj) * ( e1e2t(ji ,jj ) * ssha_e(ji ,jj ) &
565	& + e1e2t(ji+1,jj ) * ssha_e(ji+1,jj ) ) * ssumask(ji,jj)
566	zsshv_a(ji,jj) = r1_2 * r1_e1e2v(ji,jj) * ( e1e2t(ji ,jj ) * ssha_e(ji ,jj ) &
567	& + e1e2t(ji ,jj+1) * ssha_e(ji ,jj+1) ) * ssvmask(ji,jj)
568	END_2D
569	#endif
570	ENDIF
571	!
572	! Half-step back interpolation of SSH for surface pressure computation at step jit+1/2
573	!-- m+1/2 m+1 m m-1 m-2 --!
574	!-- ssh' = za0 * ssh + za1 * ssh + za2 * ssh + za3 * ssh --!
575	!------------------------------------------------------------------------------------------!
576	CALL ts_bck_interp( jn, ll_init, za0, za1, za2, za3 ) ! coeficients of the interpolation
577	zsshp2_e(:,:) = za0 * ssha_e(:,:) + za1 * sshn_e (:,:) &
578	& + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:)
579	!
580	! ! Surface pressure gradient
581	zldg = ( 1._wp - rn_scal_load ) * grav ! local factor
582	DO_2D( 0, 0, 0, 0 )
583	zu_spg(ji,jj) = - zldg * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj)
584	zv_spg(ji,jj) = - zldg * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj)
585	END_2D
586	IF( ln_wd_il ) THEN ! W/D : gravity filters applied on pressure gradient
587	CALL wad_spg( zsshp2_e, zcpx, zcpy ) ! Calculating W/D gravity filters
588	zu_spg(2:jpim1,2:jpjm1) = zu_spg(2:jpim1,2:jpjm1) * zcpx(2:jpim1,2:jpjm1)
589	zv_spg(2:jpim1,2:jpjm1) = zv_spg(2:jpim1,2:jpjm1) * zcpy(2:jpim1,2:jpjm1)
590	ENDIF
591	!
592	! Add Coriolis trend:
593	! zwz array below or triads normally depend on sea level with ln_linssh=F and should be updated
594	! at each time step. We however keep them constant here for optimization.
595	! Recall that zhU and zhV hold fluxes at jn+0.5 (extrapolated not backward interpolated)
596	CALL dyn_cor_2d( zhtp2_e, zhup2_e, zhvp2_e, ua_e, va_e, zhU, zhV, zu_trd, zv_trd )
597	!
598	! Add tidal astronomical forcing if defined
599	IF ( ln_tide .AND. ln_tide_pot ) THEN
600	DO_2D( 0, 0, 0, 0 )
601	zu_trd(ji,jj) = zu_trd(ji,jj) + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) * r1_e1u(ji,jj)
602	zv_trd(ji,jj) = zv_trd(ji,jj) + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) * r1_e2v(ji,jj)
603	END_2D
604	ENDIF
605	!
606	! Add bottom stresses:
607	!jth do implicitly instead
608	IF ( .NOT. ll_wd ) THEN ! Revert to explicit for bit comparison tests in non wad runs
609	DO_2D( 0, 0, 0, 0 )
610	zu_trd(ji,jj) = zu_trd(ji,jj) + zCdU_u(ji,jj) * un_e(ji,jj) * hur_e(ji,jj)
611	zv_trd(ji,jj) = zv_trd(ji,jj) + zCdU_v(ji,jj) * vn_e(ji,jj) * hvr_e(ji,jj)
612	END_2D
613	ENDIF
614	!
615	! Set next velocities:
616	! Compute barotropic speeds at step jit+1 (h : total height of the water colomn)
617	!-- VECTOR FORM
618	!-- m+1 m / m+1/2 \ --!
619	!-- u = u + delta_t' * \ (1-r)g grad_x( ssh') - f * k vect u + frc / --!
620	!-- --!
621	!-- FLUX FORM --!
622	!-- m+1 __1__ / m m / m+1/2 m+1/2 m+1/2 n \ \ --!
623	!-- u = m+1 \| h * u + delta_t' * \ h * (1-r)g grad_x( ssh') - h * f * k vect u + h * frc / \| --!
624	!-- h \ / --!
625	!------------------------------------------------------------------------------------------------------------------------!
626	IF( ln_dynadv_vec .OR. ln_linssh ) THEN !* Vector form
627	DO_2D( 0, 0, 0, 0 )
628	ua_e(ji,jj) = ( un_e(ji,jj) &
629	& + rDt_e * ( zu_spg(ji,jj) &
630	& + zu_trd(ji,jj) &
631	& + zu_frc(ji,jj) ) &
632	& ) * ssumask(ji,jj)
633
634	va_e(ji,jj) = ( vn_e(ji,jj) &
635	& + rDt_e * ( zv_spg(ji,jj) &
636	& + zv_trd(ji,jj) &
637	& + zv_frc(ji,jj) ) &
638	& ) * ssvmask(ji,jj)
639	END_2D
640	!
641	ELSE !* Flux form
642	DO_2D( 0, 0, 0, 0 )
643	! ! hu_e, hv_e hold depth at jn, zhup2_e, zhvp2_e hold extrapolated depth at jn+1/2
644	! ! backward interpolated depth used in spg terms at jn+1/2
645	#if defined key_qcoTest_FluxForm
646	! ! 'key_qcoTest_FluxForm' : simple ssh average
647	zhu_bck = hu_0(ji,jj) + r1_2 * ( zsshp2_e(ji,jj) + zsshp2_e(ji+1,jj ) ) * ssumask(ji,jj)
648	zhv_bck = hv_0(ji,jj) + r1_2 * ( zsshp2_e(ji,jj) + zsshp2_e(ji ,jj+1) ) * ssvmask(ji,jj)
649	#else
650	zhu_bck = hu_0(ji,jj) + r1_2r1_e1e2u(ji,jj) ( e1e2t(ji ,jj) * zsshp2_e(ji ,jj) &
651	& + e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj) ) * ssumask(ji,jj)
652	zhv_bck = hv_0(ji,jj) + r1_2r1_e1e2v(ji,jj) ( e1e2t(ji,jj ) * zsshp2_e(ji,jj ) &
653	& + e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1) ) * ssvmask(ji,jj)
654	#endif
655	! ! inverse depth at jn+1
656	z1_hu = ssumask(ji,jj) / ( hu_0(ji,jj) + zsshu_a(ji,jj) + 1._wp - ssumask(ji,jj) )
657	z1_hv = ssvmask(ji,jj) / ( hv_0(ji,jj) + zsshv_a(ji,jj) + 1._wp - ssvmask(ji,jj) )
658	!
659	ua_e(ji,jj) = ( hu_e (ji,jj) * un_e (ji,jj) &
660	& + rDt_e * ( zhu_bck * zu_spg (ji,jj) & !
661	& + zhup2_e(ji,jj) * zu_trd (ji,jj) & !
662	& + hu(ji,jj,Kmm) * zu_frc (ji,jj) ) ) * z1_hu
663	!
664	va_e(ji,jj) = ( hv_e (ji,jj) * vn_e (ji,jj) &
665	& + rDt_e * ( zhv_bck * zv_spg (ji,jj) & !
666	& + zhvp2_e(ji,jj) * zv_trd (ji,jj) & !
667	& + hv(ji,jj,Kmm) * zv_frc (ji,jj) ) ) * z1_hv
668	END_2D
669	ENDIF
670	!jth implicit bottom friction:
671	IF ( ll_wd ) THEN ! revert to explicit for bit comparison tests in non wad runs
672	DO_2D( 0, 0, 0, 0 )
673	ua_e(ji,jj) = ua_e(ji,jj) / ( 1._wp - rDt_e * zCdU_u(ji,jj) * hur_e(ji,jj) )
674	va_e(ji,jj) = va_e(ji,jj) / ( 1._wp - rDt_e * zCdU_v(ji,jj) * hvr_e(ji,jj) )
675	END_2D
676	ENDIF
677
678	IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only)
679	hu_e (2:jpim1,2:jpjm1) = hu_0(2:jpim1,2:jpjm1) + zsshu_a(2:jpim1,2:jpjm1)
680	hur_e(2:jpim1,2:jpjm1) = ssumask(2:jpim1,2:jpjm1) / ( hu_e(2:jpim1,2:jpjm1) + 1._wp - ssumask(2:jpim1,2:jpjm1) )
681	hv_e (2:jpim1,2:jpjm1) = hv_0(2:jpim1,2:jpjm1) + zsshv_a(2:jpim1,2:jpjm1)
682	hvr_e(2:jpim1,2:jpjm1) = ssvmask(2:jpim1,2:jpjm1) / ( hv_e(2:jpim1,2:jpjm1) + 1._wp - ssvmask(2:jpim1,2:jpjm1) )
683	ENDIF
684	!
685	IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only)
686	CALL lbc_lnk_multi( 'dynspg_ts', ua_e , 'U', -1._wp, va_e , 'V', -1._wp &
687	& , hu_e , 'U', 1._wp, hv_e , 'V', 1._wp &
688	& , hur_e, 'U', 1._wp, hvr_e, 'V', 1._wp )
689	ELSE
690	CALL lbc_lnk_multi( 'dynspg_ts', ua_e , 'U', -1._wp, va_e , 'V', -1._wp )
691	ENDIF
692	! ! open boundaries
693	IF( ln_bdy ) CALL bdy_dyn2d( jn, ua_e, va_e, un_e, vn_e, hur_e, hvr_e, ssha_e )
694	#if defined key_agrif
695	IF( .NOT.Agrif_Root() ) CALL agrif_dyn_ts( jn ) ! Agrif
696	#endif
697	! !* Swap
698	! ! ----
699	ubb_e (:,:) = ub_e (:,:)
700	ub_e (:,:) = un_e (:,:)
701	un_e (:,:) = ua_e (:,:)
702	!
703	vbb_e (:,:) = vb_e (:,:)
704	vb_e (:,:) = vn_e (:,:)
705	vn_e (:,:) = va_e (:,:)
706	!
707	sshbb_e(:,:) = sshb_e(:,:)
708	sshb_e (:,:) = sshn_e(:,:)
709	sshn_e (:,:) = ssha_e(:,:)
710
711	! !* Sum over whole bt loop
712	! ! ----------------------
713	za1 = wgtbtp1(jn)
714	IF( ln_dynadv_vec .OR. ln_linssh ) THEN ! Sum velocities
715	puu_b (:,:,Kaa) = puu_b (:,:,Kaa) + za1 * ua_e (:,:)
716	pvv_b (:,:,Kaa) = pvv_b (:,:,Kaa) + za1 * va_e (:,:)
717	ELSE ! Sum transports
718	IF ( .NOT.ln_wd_dl ) THEN
719	puu_b (:,:,Kaa) = puu_b (:,:,Kaa) + za1 * ua_e (:,:) * hu_e (:,:)
720	pvv_b (:,:,Kaa) = pvv_b (:,:,Kaa) + za1 * va_e (:,:) * hv_e (:,:)
721	ELSE
722	puu_b (:,:,Kaa) = puu_b (:,:,Kaa) + za1 * ua_e (:,:) * hu_e (:,:) * zuwdmask(:,:)
723	pvv_b (:,:,Kaa) = pvv_b (:,:,Kaa) + za1 * va_e (:,:) * hv_e (:,:) * zvwdmask(:,:)
724	END IF
725	ENDIF
726	! ! Sum sea level
727	pssh(:,:,Kaa) = pssh(:,:,Kaa) + za1 * ssha_e(:,:)
728
729	! ! ==================== !
730	END DO ! end loop !
731	! ! ==================== !
732	! -----------------------------------------------------------------------------
733	! Phase 3. update the general trend with the barotropic trend
734	! -----------------------------------------------------------------------------
735	!
736	! Set advection velocity correction:
737	IF (ln_bt_fw) THEN
738	IF( .NOT.( kt == nit000 .AND. l_1st_euler ) ) THEN
739	DO_2D( 1, 1, 1, 1 )
740	zun_save = un_adv(ji,jj)
741	zvn_save = vn_adv(ji,jj)
742	! ! apply the previously computed correction
743	un_adv(ji,jj) = r1_2 * ( ub2_b(ji,jj) + zun_save - rn_atfp * un_bf(ji,jj) )
744	vn_adv(ji,jj) = r1_2 * ( vb2_b(ji,jj) + zvn_save - rn_atfp * vn_bf(ji,jj) )
745	! ! Update corrective fluxes for next time step
746	un_bf(ji,jj) = rn_atfp * un_bf(ji,jj) + ( zun_save - ub2_b(ji,jj) )
747	vn_bf(ji,jj) = rn_atfp * vn_bf(ji,jj) + ( zvn_save - vb2_b(ji,jj) )
748	! ! Save integrated transport for next computation
749	ub2_b(ji,jj) = zun_save
750	vb2_b(ji,jj) = zvn_save
751	END_2D
752	ELSE
753	un_bf(:,:) = 0._wp ! corrective fluxes for next time step set to zero
754	vn_bf(:,:) = 0._wp
755	ub2_b(:,:) = un_adv(:,:) ! Save integrated transport for next computation
756	vb2_b(:,:) = vn_adv(:,:)
757	END IF
758	ENDIF
759
760
761	!
762	! Update barotropic trend:
763	IF( ln_dynadv_vec .OR. ln_linssh ) THEN
764	DO jk=1,jpkm1
765	puu(:,:,jk,Krhs) = puu(:,:,jk,Krhs) + ( puu_b(:,:,Kaa) - puu_b(:,:,Kbb) ) * r1_Dt_b
766	pvv(:,:,jk,Krhs) = pvv(:,:,jk,Krhs) + ( pvv_b(:,:,Kaa) - pvv_b(:,:,Kbb) ) * r1_Dt_b
767	END DO
768	ELSE
769	! At this stage, pssh(:,:,:,Krhs) has been corrected: compute new depths at velocity points
770	#if defined key_qcoTest_FluxForm
771	! ! 'key_qcoTest_FluxForm' : simple ssh average
772	DO_2D( 1, 0, 1, 0 )
773	zsshu_a(ji,jj) = r1_2 * ( pssh(ji,jj,Kaa) + pssh(ji+1,jj ,Kaa) ) * ssumask(ji,jj)
774	zsshv_a(ji,jj) = r1_2 * ( pssh(ji,jj,Kaa) + pssh(ji ,jj+1,Kaa) ) * ssvmask(ji,jj)
775	END_2D
776	#else
777	DO_2D( 1, 0, 1, 0 )
778	zsshu_a(ji,jj) = r1_2 * r1_e1e2u(ji,jj) * ( e1e2t(ji ,jj) * pssh(ji ,jj,Kaa) &
779	& + e1e2t(ji+1,jj) * pssh(ji+1,jj,Kaa) ) * ssumask(ji,jj)
780	zsshv_a(ji,jj) = r1_2 * r1_e1e2v(ji,jj) * ( e1e2t(ji,jj ) * pssh(ji,jj ,Kaa) &
781	& + e1e2t(ji,jj+1) * pssh(ji,jj+1,Kaa) ) * ssvmask(ji,jj)
782	END_2D
783	#endif
784	CALL lbc_lnk_multi( 'dynspg_ts', zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) ! Boundary conditions
785	!
786	DO jk=1,jpkm1
787	puu(:,:,jk,Krhs) = puu(:,:,jk,Krhs) + r1_hu(:,:,Kmm) &
788	& * ( puu_b(:,:,Kaa) - puu_b(:,:,Kbb) * hu(:,:,Kbb) ) * r1_Dt_b
789	pvv(:,:,jk,Krhs) = pvv(:,:,jk,Krhs) + r1_hv(:,:,Kmm) &
790	& * ( pvv_b(:,:,Kaa) - pvv_b(:,:,Kbb) * hv(:,:,Kbb) ) * r1_Dt_b
791	END DO
792	! Save barotropic velocities not transport:
793	puu_b(:,:,Kaa) = puu_b(:,:,Kaa) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - ssumask(:,:) )
794	pvv_b(:,:,Kaa) = pvv_b(:,:,Kaa) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - ssvmask(:,:) )
795	ENDIF
796
797
798	! Correct velocities so that the barotropic velocity equals (un_adv, vn_adv) (in all cases)
799	DO jk = 1, jpkm1
800	puu(:,:,jk,Kmm) = ( puu(:,:,jk,Kmm) + un_adv(:,:)r1_hu(:,:,Kmm) - puu_b(:,:,Kmm) ) umask(:,:,jk)
801	pvv(:,:,jk,Kmm) = ( pvv(:,:,jk,Kmm) + vn_adv(:,:)r1_hv(:,:,Kmm) - pvv_b(:,:,Kmm) ) vmask(:,:,jk)
802	END DO
803
804	IF ( ln_wd_dl .and. ln_wd_dl_bc) THEN
805	DO jk = 1, jpkm1
806	puu(:,:,jk,Kmm) = ( un_adv(:,:)*r1_hu(:,:,Kmm) &
807	& + zuwdav2(:,:)(puu(:,:,jk,Kmm) - un_adv(:,:)r1_hu(:,:,Kmm)) ) * umask(:,:,jk)
808	pvv(:,:,jk,Kmm) = ( vn_adv(:,:)*r1_hv(:,:,Kmm) &
809	& + zvwdav2(:,:)(pvv(:,:,jk,Kmm) - vn_adv(:,:)r1_hv(:,:,Kmm)) ) * vmask(:,:,jk)
810	END DO
811	END IF
812
813
814	CALL iom_put( "ubar", un_adv(:,:)*r1_hu(:,:,Kmm) ) ! barotropic i-current
815	CALL iom_put( "vbar", vn_adv(:,:)*r1_hv(:,:,Kmm) ) ! barotropic i-current
816	!
817	#if defined key_agrif
818	! Save time integrated fluxes during child grid integration
819	! (used to update coarse grid transports at next time step)
820	!
821	IF( .NOT.Agrif_Root() .AND. ln_bt_fw ) THEN
822	IF( Agrif_NbStepint() == 0 ) THEN
823	ub2_i_b(:,:) = 0._wp
824	vb2_i_b(:,:) = 0._wp
825	END IF
826	!
827	za1 = 1._wp / REAL(Agrif_rhot(), wp)
828	ub2_i_b(:,:) = ub2_i_b(:,:) + za1 * ub2_b(:,:)
829	vb2_i_b(:,:) = vb2_i_b(:,:) + za1 * vb2_b(:,:)
830	ENDIF
831	#endif
832	! !* write time-spliting arrays in the restart
833	IF( lrst_oce .AND.ln_bt_fw ) CALL ts_rst( kt, 'WRITE' )
834	!
835	IF( ln_wd_il ) DEALLOCATE( zcpx, zcpy )
836	IF( ln_wd_dl ) DEALLOCATE( ztwdmask, zuwdmask, zvwdmask, zuwdav2, zvwdav2 )
837	!
838	CALL iom_put( "baro_u" , puu_b(:,:,Kmm) ) ! Barotropic U Velocity
839	CALL iom_put( "baro_v" , pvv_b(:,:,Kmm) ) ! Barotropic V Velocity
840	!
841	END SUBROUTINE dyn_spg_ts
842
843
844	SUBROUTINE ts_wgt( ll_av, ll_fw, jpit, zwgt1, zwgt2)
845	!!---------------------------------------------------------------------
846	!! * ROUTINE ts_wgt *
847	!!
848	!! ** Purpose : Set time-splitting weights for temporal averaging (or not)
849	!!----------------------------------------------------------------------
850	LOGICAL, INTENT(in) :: ll_av ! temporal averaging=.true.
851	LOGICAL, INTENT(in) :: ll_fw ! forward time splitting =.true.
852	INTEGER, INTENT(inout) :: jpit ! cycle length
853	REAL(wp), DIMENSION(3*nn_e), INTENT(inout) :: zwgt1, & ! Primary weights
854	zwgt2 ! Secondary weights
855
856	INTEGER :: jic, jn, ji ! temporary integers
857	REAL(wp) :: za1, za2
858	!!----------------------------------------------------------------------
859
860	zwgt1(:) = 0._wp
861	zwgt2(:) = 0._wp
862
863	! Set time index when averaged value is requested
864	IF (ll_fw) THEN
865	jic = nn_e
866	ELSE
867	jic = 2 * nn_e
868	ENDIF
869
870	! Set primary weights:
871	IF (ll_av) THEN
872	! Define simple boxcar window for primary weights
873	! (width = nn_e, centered around jic)
874	SELECT CASE ( nn_bt_flt )
875	CASE( 0 ) ! No averaging
876	zwgt1(jic) = 1._wp
877	jpit = jic
878
879	CASE( 1 ) ! Boxcar, width = nn_e
880	DO jn = 1, 3*nn_e
881	za1 = ABS(float(jn-jic))/float(nn_e)
882	IF (za1 < 0.5_wp) THEN
883	zwgt1(jn) = 1._wp
884	jpit = jn
885	ENDIF
886	ENDDO
887
888	CASE( 2 ) ! Boxcar, width = 2 * nn_e
889	DO jn = 1, 3*nn_e
890	za1 = ABS(float(jn-jic))/float(nn_e)
891	IF (za1 < 1._wp) THEN
892	zwgt1(jn) = 1._wp
893	jpit = jn
894	ENDIF
895	ENDDO
896	CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_bt_flt' )
897	END SELECT
898
899	ELSE ! No time averaging
900	zwgt1(jic) = 1._wp
901	jpit = jic
902	ENDIF
903
904	! Set secondary weights
905	DO jn = 1, jpit
906	DO ji = jn, jpit
907	zwgt2(jn) = zwgt2(jn) + zwgt1(ji)
908	END DO
909	END DO
910
911	! Normalize weigths:
912	za1 = 1._wp / SUM(zwgt1(1:jpit))
913	za2 = 1._wp / SUM(zwgt2(1:jpit))
914	DO jn = 1, jpit
915	zwgt1(jn) = zwgt1(jn) * za1
916	zwgt2(jn) = zwgt2(jn) * za2
917	END DO
918	!
919	END SUBROUTINE ts_wgt
920
921
922	SUBROUTINE ts_rst( kt, cdrw )
923	!!---------------------------------------------------------------------
924	!! * ROUTINE ts_rst *
925	!!
926	!! ** Purpose : Read or write time-splitting arrays in restart file
927	!!----------------------------------------------------------------------
928	INTEGER , INTENT(in) :: kt ! ocean time-step
929	CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag
930	!!----------------------------------------------------------------------
931	!
932	IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialise
933	! ! ---------------
934	IF( ln_rstart .AND. ln_bt_fw .AND. .NOT.l_1st_euler ) THEN !* Read the restart file
935	CALL iom_get( numror, jpdom_auto, 'ub2_b' , ub2_b (:,:), cd_type = 'U', psgn = -1._wp )
936	CALL iom_get( numror, jpdom_auto, 'vb2_b' , vb2_b (:,:), cd_type = 'V', psgn = -1._wp )
937	CALL iom_get( numror, jpdom_auto, 'un_bf' , un_bf (:,:), cd_type = 'U', psgn = -1._wp )
938	CALL iom_get( numror, jpdom_auto, 'vn_bf' , vn_bf (:,:), cd_type = 'V', psgn = -1._wp )
939	IF( .NOT.ln_bt_av ) THEN
940	CALL iom_get( numror, jpdom_auto, 'sshbb_e' , sshbb_e(:,:), cd_type = 'T', psgn = 1._wp )
941	CALL iom_get( numror, jpdom_auto, 'ubb_e' , ubb_e(:,:), cd_type = 'U', psgn = -1._wp )
942	CALL iom_get( numror, jpdom_auto, 'vbb_e' , vbb_e(:,:), cd_type = 'V', psgn = -1._wp )
943	CALL iom_get( numror, jpdom_auto, 'sshb_e' , sshb_e(:,:), cd_type = 'T', psgn = 1._wp )
944	CALL iom_get( numror, jpdom_auto, 'ub_e' , ub_e(:,:), cd_type = 'U', psgn = -1._wp )
945	CALL iom_get( numror, jpdom_auto, 'vb_e' , vb_e(:,:), cd_type = 'V', psgn = -1._wp )
946	ENDIF
947	#if defined key_agrif
948	! Read time integrated fluxes
949	IF ( .NOT.Agrif_Root() ) THEN
950	CALL iom_get( numror, jpdom_auto, 'ub2_i_b' , ub2_i_b(:,:), cd_type = 'U', psgn = -1._wp )
951	CALL iom_get( numror, jpdom_auto, 'vb2_i_b' , vb2_i_b(:,:), cd_type = 'V', psgn = -1._wp )
952	ELSE
953	ub2_i_b(:,:) = 0._wp ; vb2_i_b(:,:) = 0._wp ! used in the 1st update of agrif
954	ENDIF
955	#endif
956	ELSE !* Start from rest
957	IF(lwp) WRITE(numout,*)
958	IF(lwp) WRITE(numout,*) ' ==>>> start from rest: set barotropic values to 0'
959	ub2_b (:,:) = 0._wp ; vb2_b (:,:) = 0._wp ! used in the 1st interpol of agrif
960	un_adv (:,:) = 0._wp ; vn_adv (:,:) = 0._wp ! used in the 1st interpol of agrif
961	un_bf (:,:) = 0._wp ; vn_bf (:,:) = 0._wp ! used in the 1st update of agrif
962	#if defined key_agrif
963	ub2_i_b(:,:) = 0._wp ; vb2_i_b(:,:) = 0._wp ! used in the 1st update of agrif
964	#endif
965	ENDIF
966	!
967	ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN ! Create restart file
968	! ! -------------------
969	IF(lwp) WRITE(numout,*) '---- ts_rst ----'
970	CALL iom_rstput( kt, nitrst, numrow, 'ub2_b' , ub2_b (:,:) )
971	CALL iom_rstput( kt, nitrst, numrow, 'vb2_b' , vb2_b (:,:) )
972	CALL iom_rstput( kt, nitrst, numrow, 'un_bf' , un_bf (:,:) )
973	CALL iom_rstput( kt, nitrst, numrow, 'vn_bf' , vn_bf (:,:) )
974	!
975	IF (.NOT.ln_bt_av) THEN
976	CALL iom_rstput( kt, nitrst, numrow, 'sshbb_e' , sshbb_e(:,:) )
977	CALL iom_rstput( kt, nitrst, numrow, 'ubb_e' , ubb_e(:,:) )
978	CALL iom_rstput( kt, nitrst, numrow, 'vbb_e' , vbb_e(:,:) )
979	CALL iom_rstput( kt, nitrst, numrow, 'sshb_e' , sshb_e(:,:) )
980	CALL iom_rstput( kt, nitrst, numrow, 'ub_e' , ub_e(:,:) )
981	CALL iom_rstput( kt, nitrst, numrow, 'vb_e' , vb_e(:,:) )
982	ENDIF
983	#if defined key_agrif
984	! Save time integrated fluxes
985	IF ( .NOT.Agrif_Root() ) THEN
986	CALL iom_rstput( kt, nitrst, numrow, 'ub2_i_b' , ub2_i_b(:,:) )
987	CALL iom_rstput( kt, nitrst, numrow, 'vb2_i_b' , vb2_i_b(:,:) )
988	ENDIF
989	#endif
990	ENDIF
991	!
992	END SUBROUTINE ts_rst
993
994
995	SUBROUTINE dyn_spg_ts_init
996	!!---------------------------------------------------------------------
997	!! * ROUTINE dyn_spg_ts_init *
998	!!
999	!! ** Purpose : Set time splitting options
1000	!!----------------------------------------------------------------------
1001	INTEGER :: ji ,jj ! dummy loop indices
1002	REAL(wp) :: zxr2, zyr2, zcmax ! local scalar
1003	REAL(wp), DIMENSION(jpi,jpj) :: zcu
1004	!!----------------------------------------------------------------------
1005	!
1006	! Max courant number for ext. grav. waves
1007	!
1008	DO_2D( 0, 0, 0, 0 )
1009	zxr2 = r1_e1t(ji,jj) * r1_e1t(ji,jj)
1010	zyr2 = r1_e2t(ji,jj) * r1_e2t(ji,jj)
1011	zcu(ji,jj) = SQRT( grav * MAX(ht_0(ji,jj),0._wp) * (zxr2 + zyr2) )
1012	END_2D
1013	!
1014	zcmax = MAXVAL( zcu(Nis0:Nie0,Njs0:Nje0) )
1015	CALL mpp_max( 'dynspg_ts', zcmax )
1016
1017	! Estimate number of iterations to satisfy a max courant number= rn_bt_cmax
1018	IF( ln_bt_auto ) nn_e = CEILING( rn_Dt / rn_bt_cmax * zcmax)
1019
1020	rDt_e = rn_Dt / REAL( nn_e , wp )
1021	zcmax = zcmax * rDt_e
1022	! Print results
1023	IF(lwp) WRITE(numout,*)
1024	IF(lwp) WRITE(numout,*) 'dyn_spg_ts_init : split-explicit free surface'
1025	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~'
1026	IF( ln_bt_auto ) THEN
1027	IF(lwp) WRITE(numout,*) ' ln_ts_auto =.true. Automatically set nn_e '
1028	IF(lwp) WRITE(numout,*) ' Max. courant number allowed: ', rn_bt_cmax
1029	ELSE
1030	IF(lwp) WRITE(numout,*) ' ln_ts_auto=.false.: Use nn_e in namelist nn_e = ', nn_e
1031	ENDIF
1032
1033	IF(ln_bt_av) THEN
1034	IF(lwp) WRITE(numout,*) ' ln_bt_av =.true. ==> Time averaging over nn_e time steps is on '
1035	ELSE
1036	IF(lwp) WRITE(numout,*) ' ln_bt_av =.false. => No time averaging of barotropic variables '
1037	ENDIF
1038	!
1039	!
1040	IF(ln_bt_fw) THEN
1041	IF(lwp) WRITE(numout,*) ' ln_bt_fw=.true. => Forward integration of barotropic variables '
1042	ELSE
1043	IF(lwp) WRITE(numout,*) ' ln_bt_fw =.false.=> Centred integration of barotropic variables '
1044	ENDIF
1045	!
1046	#if defined key_agrif
1047	! Restrict the use of Agrif to the forward case only
1048	!!! IF( .NOT.ln_bt_fw .AND. .NOT.Agrif_Root() ) CALL ctl_stop( 'AGRIF not implemented if ln_bt_fw=.FALSE.' )
1049	#endif
1050	!
1051	IF(lwp) WRITE(numout,*) ' Time filter choice, nn_bt_flt: ', nn_bt_flt
1052	SELECT CASE ( nn_bt_flt )
1053	CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' Dirac'
1054	CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' Boxcar: width = nn_e'
1055	CASE( 2 ) ; IF(lwp) WRITE(numout,) ' Boxcar: width = 2nn_e'
1056	CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_bt_flt: should 0,1, or 2' )
1057	END SELECT
1058	!
1059	IF(lwp) WRITE(numout,*) ' '
1060	IF(lwp) WRITE(numout,*) ' nn_e = ', nn_e
1061	IF(lwp) WRITE(numout,*) ' Barotropic time step [s] is :', rDt_e
1062	IF(lwp) WRITE(numout,*) ' Maximum Courant number is :', zcmax
1063	!
1064	IF(lwp) WRITE(numout,*) ' Time diffusion parameter rn_bt_alpha: ', rn_bt_alpha
1065	IF ((ln_bt_av.AND.nn_bt_flt/=0).AND.(rn_bt_alpha>0._wp)) THEN
1066	CALL ctl_stop( 'dynspg_ts ERROR: if rn_bt_alpha > 0, remove temporal averaging' )
1067	ENDIF
1068	!
1069	IF( .NOT.ln_bt_av .AND. .NOT.ln_bt_fw ) THEN
1070	CALL ctl_stop( 'dynspg_ts ERROR: No time averaging => only forward integration is possible' )
1071	ENDIF
1072	IF( zcmax>0.9_wp ) THEN
1073	CALL ctl_stop( 'dynspg_ts ERROR: Maximum Courant number is greater than 0.9: Inc. nn_e !' )
1074	ENDIF
1075	!
1076	! ! Allocate time-splitting arrays
1077	IF( dyn_spg_ts_alloc() /= 0 ) CALL ctl_stop('STOP', 'dyn_spg_init: failed to allocate dynspg_ts arrays' )
1078	!
1079	! ! read restart when needed
1080	CALL ts_rst( nit000, 'READ' )
1081	!
1082	END SUBROUTINE dyn_spg_ts_init
1083
1084
1085	SUBROUTINE dyn_cor_2D_init( Kmm )
1086	!!---------------------------------------------------------------------
1087	!! * ROUTINE dyn_cor_2D_init *
1088	!!
1089	!! ** Purpose : Set time splitting options
1090	!! Set arrays to remove/compute coriolis trend.
1091	!! Do it once during initialization if volume is fixed, else at each long time step.
1092	!! Note that these arrays are also used during barotropic loop. These are however frozen
1093	!! although they should be updated in the variable volume case. Not a big approximation.
1094	!! To remove this approximation, copy lines below inside barotropic loop
1095	!! and update depths at T- points (ht) at each barotropic time step
1096	!!
1097	!! Compute zwz = f / ( height of the water colomn )
1098	!!----------------------------------------------------------------------
1099	INTEGER, INTENT(in) :: Kmm ! Time index
1100	INTEGER :: ji ,jj, jk ! dummy loop indices
1101	REAL(wp) :: z1_ht
1102	!!----------------------------------------------------------------------
1103	!
1104	SELECT CASE( nvor_scheme )
1105	CASE( np_EEN, np_ENE, np_ENS , np_MIX ) != schemes using the same e3f definition
1106	SELECT CASE( nn_e3f_typ ) !* ff_f/e3 at F-point
1107	CASE ( 0 ) ! original formulation (masked averaging of e3t divided by 4)
1108	DO_2D( 0, 0, 0, 0 )
1109	zwz(ji,jj) = ( ht(ji,jj+1) + ht(ji+1,jj+1) &
1110	& + ht(ji,jj ) + ht(ji+1,jj ) ) * 0.25_wp
1111	IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj)
1112	END_2D
1113	CASE ( 1 ) ! new formulation (masked averaging of e3t divided by the sum of mask)
1114	DO_2D( 0, 0, 0, 0 )
1115	zwz(ji,jj) = ( ht(ji,jj+1) + ht(ji+1,jj+1) &
1116	& + ht(ji,jj ) + ht(ji+1,jj ) ) &
1117	& / ( MAX(ssmask(ji,jj+1) + ssmask(ji+1,jj+1) &
1118	& + ssmask(ji,jj ) + ssmask(ji+1,jj ) , 1._wp ) )
1119	IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj)
1120	END_2D
1121	END SELECT
1122	CALL lbc_lnk( 'dynspg_ts', zwz, 'F', 1._wp )
1123	END SELECT
1124	!
1125	SELECT CASE( nvor_scheme )
1126	CASE( np_EEN )
1127	!
1128	ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp
1129	DO_2D( 0, 1, 0, 1 )
1130	ftne(ji,jj) = zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1)
1131	ftnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj ) + zwz(ji ,jj )
1132	ftse(ji,jj) = zwz(ji ,jj ) + zwz(ji ,jj-1) + zwz(ji-1,jj-1)
1133	ftsw(ji,jj) = zwz(ji ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj )
1134	END_2D
1135	!
1136	CASE( np_EET ) != EEN scheme using e3t energy conserving scheme
1137	ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp
1138	DO_2D( 0, 1, 0, 1 )
1139	z1_ht = ssmask(ji,jj) / ( ht(ji,jj) + 1._wp - ssmask(ji,jj) )
1140	ftne(ji,jj) = ( ff_f(ji-1,jj ) + ff_f(ji ,jj ) + ff_f(ji ,jj-1) ) * z1_ht
1141	ftnw(ji,jj) = ( ff_f(ji-1,jj-1) + ff_f(ji-1,jj ) + ff_f(ji ,jj ) ) * z1_ht
1142	ftse(ji,jj) = ( ff_f(ji ,jj ) + ff_f(ji ,jj-1) + ff_f(ji-1,jj-1) ) * z1_ht
1143	ftsw(ji,jj) = ( ff_f(ji ,jj-1) + ff_f(ji-1,jj-1) + ff_f(ji-1,jj ) ) * z1_ht
1144	END_2D
1145	!
1146	END SELECT
1147
1148	END SUBROUTINE dyn_cor_2d_init
1149
1150
1151	SUBROUTINE dyn_cor_2d( pht, phu, phv, punb, pvnb, zhU, zhV, zu_trd, zv_trd )
1152	!!---------------------------------------------------------------------
1153	!! * ROUTINE dyn_cor_2d *
1154	!!
1155	!! ** Purpose : Compute u and v coriolis trends
1156	!!----------------------------------------------------------------------
1157	INTEGER :: ji ,jj ! dummy loop indices
1158	REAL(wp) :: zx1, zx2, zy1, zy2, z1_hu, z1_hv ! - -
1159	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pht, phu, phv, punb, pvnb, zhU, zhV
1160	REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: zu_trd, zv_trd
1161	!!----------------------------------------------------------------------
1162	SELECT CASE( nvor_scheme )
1163	CASE( np_ENT ) ! enstrophy conserving scheme (f-point)
1164	DO_2D( 0, 0, 0, 0 )
1165	z1_hu = ssumask(ji,jj) / ( phu(ji,jj) + 1._wp - ssumask(ji,jj) )
1166	z1_hv = ssvmask(ji,jj) / ( phv(ji,jj) + 1._wp - ssvmask(ji,jj) )
1167	zu_trd(ji,jj) = + r1_4 * r1_e1e2u(ji,jj) * z1_hu &
1168	& * ( e1e2t(ji+1,jj)pht(ji+1,jj)ff_t(ji+1,jj) * ( pvnb(ji+1,jj) + pvnb(ji+1,jj-1) ) &
1169	& + e1e2t(ji ,jj)pht(ji ,jj)ff_t(ji ,jj) * ( pvnb(ji ,jj) + pvnb(ji ,jj-1) ) )
1170	!
1171	zv_trd(ji,jj) = - r1_4 * r1_e1e2v(ji,jj) * z1_hv &
1172	& * ( e1e2t(ji,jj+1)pht(ji,jj+1)ff_t(ji,jj+1) * ( punb(ji,jj+1) + punb(ji-1,jj+1) ) &
1173	& + e1e2t(ji,jj )pht(ji,jj )ff_t(ji,jj ) * ( punb(ji,jj ) + punb(ji-1,jj ) ) )
1174	END_2D
1175	!
1176	CASE( np_ENE , np_MIX ) ! energy conserving scheme (t-point) ENE or MIX
1177	DO_2D( 0, 0, 0, 0 )
1178	zy1 = ( zhV(ji,jj-1) + zhV(ji+1,jj-1) ) * r1_e1u(ji,jj)
1179	zy2 = ( zhV(ji,jj ) + zhV(ji+1,jj ) ) * r1_e1u(ji,jj)
1180	zx1 = ( zhU(ji-1,jj) + zhU(ji-1,jj+1) ) * r1_e2v(ji,jj)
1181	zx2 = ( zhU(ji ,jj) + zhU(ji ,jj+1) ) * r1_e2v(ji,jj)
1182	! energy conserving formulation for planetary vorticity term
1183	zu_trd(ji,jj) = r1_4 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
1184	zv_trd(ji,jj) = - r1_4 * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 )
1185	END_2D
1186	!
1187	CASE( np_ENS ) ! enstrophy conserving scheme (f-point)
1188	DO_2D( 0, 0, 0, 0 )
1189	zy1 = r1_8 * ( zhV(ji ,jj-1) + zhV(ji+1,jj-1) &
1190	& + zhV(ji ,jj ) + zhV(ji+1,jj ) ) * r1_e1u(ji,jj)
1191	zx1 = - r1_8 * ( zhU(ji-1,jj ) + zhU(ji-1,jj+1) &
1192	& + zhU(ji ,jj ) + zhU(ji ,jj+1) ) * r1_e2v(ji,jj)
1193	zu_trd(ji,jj) = zy1 * ( zwz(ji ,jj-1) + zwz(ji,jj) )
1194	zv_trd(ji,jj) = zx1 * ( zwz(ji-1,jj ) + zwz(ji,jj) )
1195	END_2D
1196	!
1197	CASE( np_EET , np_EEN ) ! energy & enstrophy scheme (using e3t or e3f)
1198	DO_2D( 0, 0, 0, 0 )
1199	zu_trd(ji,jj) = + r1_12 * r1_e1u(ji,jj) * ( ftne(ji,jj ) * zhV(ji ,jj ) &
1200	& + ftnw(ji+1,jj) * zhV(ji+1,jj ) &
1201	& + ftse(ji,jj ) * zhV(ji ,jj-1) &
1202	& + ftsw(ji+1,jj) * zhV(ji+1,jj-1) )
1203	zv_trd(ji,jj) = - r1_12 * r1_e2v(ji,jj) * ( ftsw(ji,jj+1) * zhU(ji-1,jj+1) &
1204	& + ftse(ji,jj+1) * zhU(ji ,jj+1) &
1205	& + ftnw(ji,jj ) * zhU(ji-1,jj ) &
1206	& + ftne(ji,jj ) * zhU(ji ,jj ) )
1207	END_2D
1208	!
1209	END SELECT
1210	!
1211	END SUBROUTINE dyn_cor_2D
1212
1213
1214	SUBROUTINE wad_tmsk( pssh, ptmsk )
1215	!!----------------------------------------------------------------------
1216	!! * ROUTINE wad_lmt *
1217	!!
1218	!! ** Purpose : set wetting & drying mask at tracer points
1219	!! for the current barotropic sub-step
1220	!!
1221	!! ** Method : ???
1222	!!
1223	!! ** Action : ptmsk : wetting & drying t-mask
1224	!!----------------------------------------------------------------------
1225	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pssh !
1226	REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: ptmsk !
1227	!
1228	INTEGER :: ji, jj ! dummy loop indices
1229	!!----------------------------------------------------------------------
1230	!
1231	IF( ln_wd_dl_rmp ) THEN
1232	DO_2D( 1, 1, 1, 1 )
1233	IF ( pssh(ji,jj) + ht_0(ji,jj) > 2._wp * rn_wdmin1 ) THEN
1234	! IF ( pssh(ji,jj) + ht_0(ji,jj) > rn_wdmin2 ) THEN
1235	ptmsk(ji,jj) = 1._wp
1236	ELSEIF( pssh(ji,jj) + ht_0(ji,jj) > rn_wdmin1 ) THEN
1237	ptmsk(ji,jj) = TANH( 50._wp( ( pssh(ji,jj) + ht_0(ji,jj) - rn_wdmin1 )r_rn_wdmin1) )
1238	ELSE
1239	ptmsk(ji,jj) = 0._wp
1240	ENDIF
1241	END_2D
1242	ELSE
1243	DO_2D( 1, 1, 1, 1 )
1244	IF ( pssh(ji,jj) + ht_0(ji,jj) > rn_wdmin1 ) THEN ; ptmsk(ji,jj) = 1._wp
1245	ELSE ; ptmsk(ji,jj) = 0._wp
1246	ENDIF
1247	END_2D
1248	ENDIF
1249	!
1250	END SUBROUTINE wad_tmsk
1251
1252
1253	SUBROUTINE wad_Umsk( pTmsk, phU, phV, pu, pv, pUmsk, pVmsk )
1254	!!----------------------------------------------------------------------
1255	!! * ROUTINE wad_lmt *
1256	!!
1257	!! ** Purpose : set wetting & drying mask at tracer points
1258	!! for the current barotropic sub-step
1259	!!
1260	!! ** Method : ???
1261	!!
1262	!! ** Action : ptmsk : wetting & drying t-mask
1263	!!----------------------------------------------------------------------
1264	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pTmsk ! W & D t-mask
1265	REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: phU, phV, pu, pv ! ocean velocities and transports
1266	REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pUmsk, pVmsk ! W & D u- and v-mask
1267	!
1268	INTEGER :: ji, jj ! dummy loop indices
1269	!!----------------------------------------------------------------------
1270	!
1271	DO_2D( 1, 1, 1, 0 ) ! not jpi-column
1272	IF ( phU(ji,jj) > 0._wp ) THEN ; pUmsk(ji,jj) = pTmsk(ji ,jj)
1273	ELSE ; pUmsk(ji,jj) = pTmsk(ji+1,jj)
1274	ENDIF
1275	phU(ji,jj) = pUmsk(ji,jj)*phU(ji,jj)
1276	pu (ji,jj) = pUmsk(ji,jj)*pu (ji,jj)
1277	END_2D
1278	!
1279	DO_2D( 1, 0, 1, 1 ) ! not jpj-row
1280	IF ( phV(ji,jj) > 0._wp ) THEN ; pVmsk(ji,jj) = pTmsk(ji,jj )
1281	ELSE ; pVmsk(ji,jj) = pTmsk(ji,jj+1)
1282	ENDIF
1283	phV(ji,jj) = pVmsk(ji,jj)*phV(ji,jj)
1284	pv (ji,jj) = pVmsk(ji,jj)*pv (ji,jj)
1285	END_2D
1286	!
1287	END SUBROUTINE wad_Umsk
1288
1289
1290	SUBROUTINE wad_spg( pshn, zcpx, zcpy )
1291	!!---------------------------------------------------------------------
1292	!! * ROUTINE wad_sp *
1293	!!
1294	!! ** Purpose :
1295	!!----------------------------------------------------------------------
1296	INTEGER :: ji ,jj ! dummy loop indices
1297	LOGICAL :: ll_tmp1, ll_tmp2
1298	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pshn
1299	REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: zcpx, zcpy
1300	!!----------------------------------------------------------------------
1301	DO_2D( 0, 0, 0, 0 )
1302	ll_tmp1 = MIN( pshn(ji,jj) , pshn(ji+1,jj) ) > &
1303	& MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) .AND. &
1304	& MAX( pshn(ji,jj) + ht_0(ji,jj) , pshn(ji+1,jj) + ht_0(ji+1,jj) ) &
1305	& > rn_wdmin1 + rn_wdmin2
1306	ll_tmp2 = ( ABS( pshn(ji+1,jj) - pshn(ji ,jj)) > 1.E-12 ).AND.( &
1307	& MAX( pshn(ji,jj) , pshn(ji+1,jj) ) > &
1308	& MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 )
1309	IF(ll_tmp1) THEN
1310	zcpx(ji,jj) = 1.0_wp
1311	ELSEIF(ll_tmp2) THEN
1312	! no worries about pshn(ji+1,jj) - pshn(ji ,jj) = 0, it won't happen ! here
1313	zcpx(ji,jj) = ABS( (pshn(ji+1,jj) + ht_0(ji+1,jj) - pshn(ji,jj) - ht_0(ji,jj)) &
1314	& / (pshn(ji+1,jj) - pshn(ji ,jj)) )
1315	zcpx(ji,jj) = max(min( zcpx(ji,jj) , 1.0_wp),0.0_wp)
1316	ELSE
1317	zcpx(ji,jj) = 0._wp
1318	ENDIF
1319	!
1320	ll_tmp1 = MIN( pshn(ji,jj) , pshn(ji,jj+1) ) > &
1321	& MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) .AND. &
1322	& MAX( pshn(ji,jj) + ht_0(ji,jj) , pshn(ji,jj+1) + ht_0(ji,jj+1) ) &
1323	& > rn_wdmin1 + rn_wdmin2
1324	ll_tmp2 = ( ABS( pshn(ji,jj) - pshn(ji,jj+1)) > 1.E-12 ).AND.( &
1325	& MAX( pshn(ji,jj) , pshn(ji,jj+1) ) > &
1326	& MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 )
1327
1328	IF(ll_tmp1) THEN
1329	zcpy(ji,jj) = 1.0_wp
1330	ELSE IF(ll_tmp2) THEN
1331	! no worries about pshn(ji,jj+1) - pshn(ji,jj ) = 0, it won't happen ! here
1332	zcpy(ji,jj) = ABS( (pshn(ji,jj+1) + ht_0(ji,jj+1) - pshn(ji,jj) - ht_0(ji,jj)) &
1333	& / (pshn(ji,jj+1) - pshn(ji,jj )) )
1334	zcpy(ji,jj) = MAX( 0._wp , MIN( zcpy(ji,jj) , 1.0_wp ) )
1335	ELSE
1336	zcpy(ji,jj) = 0._wp
1337	ENDIF
1338	END_2D
1339
1340	END SUBROUTINE wad_spg
1341
1342
1343	SUBROUTINE dyn_drg_init( Kbb, Kmm, puu, pvv, puu_b ,pvv_b, pu_RHSi, pv_RHSi, pCdU_u, pCdU_v )
1344	!!----------------------------------------------------------------------
1345	!! * ROUTINE dyn_drg_init *
1346	!!
1347	!! ** Purpose : - add the baroclinic top/bottom drag contribution to
1348	!! the baroclinic part of the barotropic RHS
1349	!! - compute the barotropic drag coefficients
1350	!!
1351	!! ** Method : computation done over the INNER domain only
1352	!!----------------------------------------------------------------------
1353	INTEGER , INTENT(in ) :: Kbb, Kmm ! ocean time level indices
1354	REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(in ) :: puu, pvv ! ocean velocities and RHS of momentum equation
1355	REAL(wp), DIMENSION(jpi,jpj,jpt) , INTENT(in ) :: puu_b, pvv_b ! barotropic velocities at main time levels
1356	REAL(wp), DIMENSION(jpi,jpj) , INTENT(inout) :: pu_RHSi, pv_RHSi ! baroclinic part of the barotropic RHS
1357	REAL(wp), DIMENSION(jpi,jpj) , INTENT( out) :: pCdU_u , pCdU_v ! barotropic drag coefficients
1358	!
1359	INTEGER :: ji, jj ! dummy loop indices
1360	INTEGER :: ikbu, ikbv, iktu, iktv
1361	REAL(wp) :: zztmp
1362	REAL(wp), DIMENSION(jpi,jpj) :: zu_i, zv_i
1363	!!----------------------------------------------------------------------
1364	!
1365	! !== Set the barotropic drag coef. ==!
1366	!
1367	IF( ln_isfcav.OR.ln_drgice_imp ) THEN ! top+bottom friction (ocean cavities)
1368
1369	DO_2D( 0, 0, 0, 0 )
1370	pCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) + rCdU_top(ji+1,jj)+rCdU_top(ji,jj) )
1371	pCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) + rCdU_top(ji,jj+1)+rCdU_top(ji,jj) )
1372	END_2D
1373	ELSE ! bottom friction only
1374	DO_2D( 0, 0, 0, 0 )
1375	pCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) )
1376	pCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) )
1377	END_2D
1378	ENDIF
1379	!
1380	! !== BOTTOM stress contribution from baroclinic velocities ==!
1381	!
1382	IF( ln_bt_fw ) THEN ! FORWARD integration: use NOW bottom baroclinic velocities
1383
1384	DO_2D( 0, 0, 0, 0 )
1385	ikbu = mbku(ji,jj)
1386	ikbv = mbkv(ji,jj)
1387	zu_i(ji,jj) = puu(ji,jj,ikbu,Kmm) - puu_b(ji,jj,Kmm)
1388	zv_i(ji,jj) = pvv(ji,jj,ikbv,Kmm) - pvv_b(ji,jj,Kmm)
1389	END_2D
1390	ELSE ! CENTRED integration: use BEFORE bottom baroclinic velocities
1391
1392	DO_2D( 0, 0, 0, 0 )
1393	ikbu = mbku(ji,jj)
1394	ikbv = mbkv(ji,jj)
1395	zu_i(ji,jj) = puu(ji,jj,ikbu,Kbb) - puu_b(ji,jj,Kbb)
1396	zv_i(ji,jj) = pvv(ji,jj,ikbv,Kbb) - pvv_b(ji,jj,Kbb)
1397	END_2D
1398	ENDIF
1399	!
1400	IF( ln_wd_il ) THEN ! W/D : use the "clipped" bottom friction !!gm explain WHY, please !
1401	zztmp = -1._wp / rDt_e
1402	DO_2D( 0, 0, 0, 0 )
1403	pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + zu_i(ji,jj) * wdrampu(ji,jj) * MAX( &
1404	& r1_hu(ji,jj,Kmm) * r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) , zztmp )
1405	pv_RHSi(ji,jj) = pv_RHSi(ji,jj) + zv_i(ji,jj) * wdrampv(ji,jj) * MAX( &
1406	& r1_hv(ji,jj,Kmm) * r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) , zztmp )
1407	END_2D
1408	ELSE ! use "unclipped" drag (even if explicit friction is used in 3D calculation)
1409
1410	DO_2D( 0, 0, 0, 0 )
1411	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)
1412	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)
1413	END_2D
1414	END IF
1415	!
1416	! !== TOP stress contribution from baroclinic velocities ==! (no W/D case)
1417	!
1418	IF( ln_isfcav.OR.ln_drgice_imp ) THEN
1419	!
1420	IF( ln_bt_fw ) THEN ! FORWARD integration: use NOW top baroclinic velocity
1421
1422	DO_2D( 0, 0, 0, 0 )
1423	iktu = miku(ji,jj)
1424	iktv = mikv(ji,jj)
1425	zu_i(ji,jj) = puu(ji,jj,iktu,Kmm) - puu_b(ji,jj,Kmm)
1426	zv_i(ji,jj) = pvv(ji,jj,iktv,Kmm) - pvv_b(ji,jj,Kmm)
1427	END_2D
1428	ELSE ! CENTRED integration: use BEFORE top baroclinic velocity
1429
1430	DO_2D( 0, 0, 0, 0 )
1431	iktu = miku(ji,jj)
1432	iktv = mikv(ji,jj)
1433	zu_i(ji,jj) = puu(ji,jj,iktu,Kbb) - puu_b(ji,jj,Kbb)
1434	zv_i(ji,jj) = pvv(ji,jj,iktv,Kbb) - pvv_b(ji,jj,Kbb)
1435	END_2D
1436	ENDIF
1437	!
1438	! ! use "unclipped" top drag (even if explicit friction is used in 3D calculation)
1439
1440	DO_2D( 0, 0, 0, 0 )
1441	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)
1442	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)
1443	END_2D
1444	!
1445	ENDIF
1446	!
1447	END SUBROUTINE dyn_drg_init
1448
1449	SUBROUTINE ts_bck_interp( jn, ll_init, & ! <== in
1450	& za0, za1, za2, za3 ) ! ==> out
1451	!!----------------------------------------------------------------------
1452	INTEGER ,INTENT(in ) :: jn ! index of sub time step
1453	LOGICAL ,INTENT(in ) :: ll_init !
1454	REAL(wp),INTENT( out) :: za0, za1, za2, za3 ! Half-step back interpolation coefficient
1455	!
1456	REAL(wp) :: zepsilon, zgamma ! - -
1457	!!----------------------------------------------------------------------
1458	! ! set Half-step back interpolation coefficient
1459	IF ( jn==1 .AND. ll_init ) THEN !* Forward-backward
1460	za0 = 1._wp
1461	za1 = 0._wp
1462	za2 = 0._wp
1463	za3 = 0._wp
1464	ELSEIF( jn==2 .AND. ll_init ) THEN !* AB2-AM3 Coefficients; bet=0 ; gam=-1/6 ; eps=1/12
1465	za0 = 1.0833333333333_wp ! za0 = 1-gam-eps
1466	za1 =-0.1666666666666_wp ! za1 = gam
1467	za2 = 0.0833333333333_wp ! za2 = eps
1468	za3 = 0._wp
1469	ELSE !* AB3-AM4 Coefficients; bet=0.281105 ; eps=0.013 ; gam=0.0880
1470	IF( rn_bt_alpha == 0._wp ) THEN ! Time diffusion
1471	za0 = 0.614_wp ! za0 = 1/2 + gam + 2*eps
1472	za1 = 0.285_wp ! za1 = 1/2 - 2gam - 3eps
1473	za2 = 0.088_wp ! za2 = gam
1474	za3 = 0.013_wp ! za3 = eps
1475	ELSE ! no time diffusion
1476	zepsilon = 0.00976186_wp - 0.13451357_wp * rn_bt_alpha
1477	zgamma = 0.08344500_wp - 0.51358400_wp * rn_bt_alpha
1478	za0 = 0.5_wp + zgamma + 2._wp * rn_bt_alpha + 2._wp * zepsilon
1479	za1 = 1._wp - za0 - zgamma - zepsilon
1480	za2 = zgamma
1481	za3 = zepsilon
1482	ENDIF
1483	ENDIF
1484	END SUBROUTINE ts_bck_interp
1485
1486
1487	!!======================================================================
1488	END MODULE dynspg_ts

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

Download in other formats: