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/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/DYN – NEMO

Context Navigation

source: NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/DYN/dynspg_ts.F90 @ 11245

Last change on this file since 11245 was 11245, checked in by girrmann, 5 years ago
dev_r10984_HPC-13 : bugfix following [11241], see #2285
Property svn:keywords set to `Id`
File size: 81.0 KB

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

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

Download in other formats: