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 branches/2017/dev_r8624_ENHANCE4_FREESURFACE/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

Context Navigation

source: branches/2017/dev_r8624_ENHANCE4_FREESURFACE/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 8805

Last change on this file since 8805 was 8805, checked in by jchanut, 7 years ago
Enhance4-freesurface. step 1: recover tracer conservation with split explicit free surface - #1959
Property svn:keywords set to `Id`
File size: 64.1 KB

Line
1	MODULE dynspg_ts
2	!!======================================================================
3	!! * MODULE dynspg_ts *
4	!! Ocean dynamics: surface pressure gradient trend, split-explicit scheme
5	!!======================================================================
6	!! History : 1.0 ! 2004-12 (L. Bessieres, G. Madec) Original code
7	!! - ! 2005-11 (V. Garnier, G. Madec) optimization
8	!! - ! 2006-08 (S. Masson) distributed restart using iom
9	!! 2.0 ! 2007-07 (D. Storkey) calls to BDY routines
10	!! - ! 2008-01 (R. Benshila) change averaging method
11	!! 3.2 ! 2009-07 (R. Benshila, G. Madec) Complete revisit associated to vvl reactivation
12	!! 3.3 ! 2010-09 (D. Storkey, E. O'Dea) update for BDY for Shelf configurations
13	!! 3.3 ! 2011-03 (R. Benshila, R. Hordoir, P. Oddo) update calculation of ub_b
14	!! 3.5 ! 2013-07 (J. Chanut) Switch to Forward-backward time stepping
15	!! 3.6 ! 2013-11 (A. Coward) Update for z-tilde compatibility
16	!! 3.7 ! 2015-11 (J. Chanut) free surface simplification
17	!! - ! 2016-12 (G. Madec, E. Clementi) update for Stoke-Drift divergence
18	!!---------------------------------------------------------------------
19
20	!!----------------------------------------------------------------------
21	!! dyn_spg_ts : compute surface pressure gradient trend using a time-splitting scheme
22	!! dyn_spg_ts_init: initialisation of the time-splitting scheme
23	!! ts_wgt : set time-splitting weights for temporal averaging (or not)
24	!! ts_rst : read/write time-splitting fields in restart file
25	!!----------------------------------------------------------------------
26	USE oce ! ocean dynamics and tracers
27	USE dom_oce ! ocean space and time domain
28	USE sbc_oce ! surface boundary condition: ocean
29	USE zdf_oce ! Bottom friction coefts
30	USE sbcisf ! ice shelf variable (fwfisf)
31	USE sbcapr ! surface boundary condition: atmospheric pressure
32	USE dynadv , ONLY: ln_dynadv_vec
33	USE phycst ! physical constants
34	USE dynvor ! vorticity term
35	USE wet_dry ! wetting/drying flux limter
36	USE bdy_oce ! open boundary
37	USE bdytides ! open boundary condition data
38	USE bdydyn2d ! open boundary conditions on barotropic variables
39	USE sbctide ! tides
40	USE updtide ! tide potential
41	USE sbcwave ! surface wave
42	!
43	USE in_out_manager ! I/O manager
44	USE lib_mpp ! distributed memory computing library
45	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
46	USE prtctl ! Print control
47	USE iom ! IOM library
48	USE restart ! only for lrst_oce
49	USE wrk_nemo ! Memory Allocation
50	USE timing ! Timing
51	USE diatmb ! Top,middle,bottom output
52	#if defined key_agrif
53	USE agrif_opa_interp ! agrif
54	#endif
55	#if defined key_asminc
56	USE asminc ! Assimilation increment
57	#endif
58
59
60	IMPLICIT NONE
61	PRIVATE
62
63	PUBLIC dyn_spg_ts ! routine called in dynspg.F90
64	PUBLIC dyn_spg_ts_alloc ! " " " "
65	PUBLIC dyn_spg_ts_init ! " " " "
66	PUBLIC ts_rst ! " " " "
67
68	INTEGER, SAVE :: icycle ! Number of barotropic sub-steps for each internal step nn_baro <= 2.5 nn_baro
69	REAL(wp),SAVE :: rdtbt ! Barotropic time step
70
71	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: wgtbtp1, wgtbtp2 !: 1st & 2nd weights used in time filtering of barotropic fields
72
73	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zwz !: ff_f/h at F points
74	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ftnw, ftne !: triad of coriolis parameter
75	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ftsw, ftse !: (only used with een vorticity scheme)
76
77	!! Time filtered arrays at baroclinic time step:
78	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: un_adv , vn_adv !: Advection vel. at "now" barocl. step
79
80	!! * Substitutions
81	# include "vectopt_loop_substitute.h90"
82	!!----------------------------------------------------------------------
83	!! NEMO/OPA 3.5 , NEMO Consortium (2013)
84	!! $Id$
85	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
86	!!----------------------------------------------------------------------
87	CONTAINS
88
89	INTEGER FUNCTION dyn_spg_ts_alloc()
90	!!----------------------------------------------------------------------
91	!! * routine dyn_spg_ts_alloc *
92	!!----------------------------------------------------------------------
93	INTEGER :: ierr(3)
94	!!----------------------------------------------------------------------
95	ierr(:) = 0
96	!
97	ALLOCATE( wgtbtp1(3nn_baro), wgtbtp2(3nn_baro), zwz(jpi,jpj), STAT=ierr(1) )
98	!
99	IF( ln_dynvor_een ) ALLOCATE( ftnw(jpi,jpj) , ftne(jpi,jpj) , &
100	& ftsw(jpi,jpj) , ftse(jpi,jpj) , STAT=ierr(2) )
101	!
102	ALLOCATE( un_adv(jpi,jpj), vn_adv(jpi,jpj) , STAT=ierr(3) )
103	!
104	dyn_spg_ts_alloc = MAXVAL( ierr(:) )
105	!
106	IF( lk_mpp ) CALL mpp_sum( dyn_spg_ts_alloc )
107	IF( dyn_spg_ts_alloc /= 0 ) CALL ctl_warn('dyn_spg_ts_alloc: failed to allocate arrays')
108	!
109	END FUNCTION dyn_spg_ts_alloc
110
111
112	SUBROUTINE dyn_spg_ts( kt )
113	!!----------------------------------------------------------------------
114	!!
115	!! ** Purpose : - Compute the now trend due to the explicit time stepping
116	!! of the quasi-linear barotropic system, and add it to the
117	!! general momentum trend.
118	!!
119	!! ** Method : - split-explicit schem (time splitting) :
120	!! Barotropic variables are advanced from internal time steps
121	!! "n" to "n+1" if ln_bt_fw=T
122	!! or from
123	!! "n-1" to "n+1" if ln_bt_fw=F
124	!! thanks to a generalized forward-backward time stepping (see ref. below).
125	!!
126	!! ** Action :
127	!! -Update the filtered free surface at step "n+1" : ssha
128	!! -Update filtered barotropic velocities at step "n+1" : ua_b, va_b
129	!! -Compute barotropic advective velocities at step "n" : un_adv, vn_adv
130	!! These are used to advect tracers and are compliant with discrete
131	!! continuity equation taken at the baroclinic time steps. This
132	!! ensures tracers conservation.
133	!! - (ua, va) momentum trend updated with barotropic component.
134	!!
135	!! References : Shchepetkin and McWilliams, Ocean Modelling, 2005.
136	!!---------------------------------------------------------------------
137	INTEGER, INTENT(in) :: kt ! ocean time-step index
138	!
139	LOGICAL :: ll_fw_start ! if true, forward integration
140	LOGICAL :: ll_init ! if true, special startup of 2d equations
141	LOGICAL :: ll_tmp1, ll_tmp2 ! local logical variables used in W/D
142	INTEGER :: ji, jj, jk, jn ! dummy loop indices
143	INTEGER :: ikbu, ikbv, noffset ! local integers
144	INTEGER :: iktu, iktv ! local integers
145	REAL(wp) :: zmdi
146	REAL(wp) :: zraur, z1_2dt_b, z2dt_bf ! local scalars
147	REAL(wp) :: zx1, zy1, zx2, zy2 ! - -
148	REAL(wp) :: z1_12, z1_8, z1_4, z1_2 ! - -
149	REAL(wp) :: zu_spg, zv_spg ! - -
150	REAL(wp) :: zhura, zhvra ! - -
151	REAL(wp) :: za0, za1, za2, za3 ! - -
152	!
153	REAL(wp), POINTER, DIMENSION(:,:) :: zsshp2_e
154	REAL(wp), POINTER, DIMENSION(:,:) :: zu_trd, zv_trd, zu_frc, zv_frc, zssh_frc
155	REAL(wp), POINTER, DIMENSION(:,:) :: zwx, zwy, zhdiv
156	REAL(wp), POINTER, DIMENSION(:,:) :: zhup2_e, zhvp2_e, zhust_e, zhvst_e
157	REAL(wp), POINTER, DIMENSION(:,:) :: zsshu_a, zsshv_a
158	REAL(wp), POINTER, DIMENSION(:,:) :: zhf
159	REAL(wp), POINTER, DIMENSION(:,:) :: zcpx, zcpy ! Wetting/Dying gravity filter coef.
160	!!----------------------------------------------------------------------
161	!
162	IF( nn_timing == 1 ) CALL timing_start('dyn_spg_ts')
163	!
164	! !* Allocate temporary arrays
165	CALL wrk_alloc( jpi,jpj, zsshp2_e, zhdiv )
166	CALL wrk_alloc( jpi,jpj, zu_trd, zv_trd)
167	CALL wrk_alloc( jpi,jpj, zwx, zwy, zssh_frc, zu_frc, zv_frc)
168	CALL wrk_alloc( jpi,jpj, zhup2_e, zhvp2_e, zhust_e, zhvst_e)
169	CALL wrk_alloc( jpi,jpj, zsshu_a, zsshv_a )
170	CALL wrk_alloc( jpi,jpj, zhf )
171	IF( ln_wd ) CALL wrk_alloc( jpi, jpj, zcpx, zcpy )
172	!
173	zmdi=1.e+20 ! missing data indicator for masking
174	! !* Local constant initialization
175	z1_12 = 1._wp / 12._wp
176	z1_8 = 0.125_wp
177	z1_4 = 0.25_wp
178	z1_2 = 0.5_wp
179	zraur = 1._wp / rau0
180	! ! reciprocal of baroclinic time step
181	IF( kt == nit000 .AND. neuler == 0 ) THEN ; z2dt_bf = rdt
182	ELSE ; z2dt_bf = 2.0_wp * rdt
183	ENDIF
184	z1_2dt_b = 1.0_wp / z2dt_bf
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_baro
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( neuler == 0 ) ll_init=.TRUE.
201	!
202	IF( ln_bt_fw .OR. neuler == 0 ) THEN
203	ll_fw_start =.TRUE.
204	noffset = 0
205	ELSE
206	ll_fw_start =.FALSE.
207	ENDIF
208	!
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	! Set arrays to remove/compute coriolis trend.
215	! Do it once at kt=nit000 if volume is fixed, else at each long time step.
216	! Note that these arrays are also used during barotropic loop. These are however frozen
217	! although they should be updated in the variable volume case. Not a big approximation.
218	! To remove this approximation, copy lines below inside barotropic loop
219	! and update depths at T-F points (ht and zhf resp.) at each barotropic time step
220	!
221	IF( kt == nit000 .OR. .NOT.ln_linssh ) THEN
222	IF( ln_dynvor_een ) THEN !== EEN scheme ==!
223	SELECT CASE( nn_een_e3f ) !* ff_f/e3 at F-point
224	CASE ( 0 ) ! original formulation (masked averaging of e3t divided by 4)
225	DO jj = 1, jpjm1
226	DO ji = 1, jpim1
227	zwz(ji,jj) = ( ht_n(ji ,jj+1) + ht_n(ji+1,jj+1) + &
228	& ht_n(ji ,jj ) + ht_n(ji+1,jj ) ) * 0.25_wp
229	IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj)
230	END DO
231	END DO
232	CASE ( 1 ) ! new formulation (masked averaging of e3t divided by the sum of mask)
233	DO jj = 1, jpjm1
234	DO ji = 1, jpim1
235	zwz(ji,jj) = ( ht_n(ji ,jj+1) + ht_n(ji+1,jj+1) + &
236	& ht_n(ji ,jj ) + ht_n(ji+1,jj ) ) &
237	& / ( MAX( 1._wp, tmask(ji ,jj+1, 1) + tmask(ji+1,jj+1, 1) + &
238	& tmask(ji ,jj , 1) + tmask(ji+1,jj , 1) ) )
239	IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj)
240	END DO
241	END DO
242	END SELECT
243	CALL lbc_lnk( zwz, 'F', 1._wp )
244	!
245	ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp
246	DO jj = 2, jpj
247	DO ji = 2, jpi
248	ftne(ji,jj) = zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1)
249	ftnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj ) + zwz(ji ,jj )
250	ftse(ji,jj) = zwz(ji ,jj ) + zwz(ji ,jj-1) + zwz(ji-1,jj-1)
251	ftsw(ji,jj) = zwz(ji ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj )
252	END DO
253	END DO
254	!
255	ELSE !== all other schemes (ENE, ENS, MIX)
256	zwz(:,:) = 0._wp
257	zhf(:,:) = 0._wp
258
259	!!gm assume 0 in both cases (xhich is almost surely WRONG ! ) as hvatf has been removed
260	!!gm A priori a better value should be something like :
261	!!gm zhf(i,j) = masked sum of ht(i,j) , ht(i+1,j) , ht(i,j+1) , (i+1,j+1)
262	!!gm divided by the sum of the corresponding mask
263	!!gm
264	!!
265	IF ( .not. ln_sco ) THEN
266
267	!!gm agree the JC comment : this should be done in a much clear way
268
269	! JC: It not clear yet what should be the depth at f-points over land in z-coordinate case
270	! Set it to zero for the time being
271	! IF( rn_hmin < 0._wp ) THEN ; jk = - INT( rn_hmin ) ! from a nb of level
272	! ELSE ; jk = MINLOC( gdepw_0, mask = gdepw_0 > rn_hmin, dim = 1 ) ! from a depth
273	! ENDIF
274	! zhf(:,:) = gdepw_0(:,:,jk+1)
275	ELSE
276	!zhf(:,:) = hbatf(:,:)
277	DO jj = 1, jpjm1
278	DO ji = 1, jpim1
279	zhf(ji,jj) = MAX( 0._wp, &
280	& ( ht_0(ji ,jj )*tmask(ji ,jj ,1) + &
281	& ht_0(ji+1,jj )*tmask(ji+1,jj ,1) + &
282	& ht_0(ji ,jj+1)*tmask(ji ,jj+1,1) + &
283	& ht_0(ji+1,jj+1)*tmask(ji+1,jj+1,1) ) / &
284	& ( tmask(ji ,jj ,1) + tmask(ji+1,jj ,1) +&
285	& tmask(ji ,jj+1,1) + tmask(ji+1,jj+1,1) +&
286	& rsmall ) )
287	END DO
288	END DO
289	END IF
290
291	DO jj = 1, jpjm1
292	zhf(:,jj) = zhf(:,jj) * (1._wp- umask(:,jj,1) * umask(:,jj+1,1))
293	END DO
294	!!gm end
295
296	DO jk = 1, jpkm1
297	DO jj = 1, jpjm1
298	zhf(:,jj) = zhf(:,jj) + e3f_n(:,jj,jk) * umask(:,jj,jk) * umask(:,jj+1,jk)
299	END DO
300	END DO
301	CALL lbc_lnk( zhf, 'F', 1._wp )
302	! JC: TBC. hf should be greater than 0
303	DO jj = 1, jpj
304	DO ji = 1, jpi
305	IF( zhf(ji,jj) /= 0._wp ) zwz(ji,jj) = 1._wp / zhf(ji,jj) ! zhf is actually hf here but it saves an array
306	END DO
307	END DO
308	zwz(:,:) = ff_f(:,:) * zwz(:,:)
309	ENDIF
310	ENDIF
311	!
312	! If forward start at previous time step, and centered integration,
313	! then update averaging weights:
314	IF (.NOT.ln_bt_fw .AND.( neuler==0 .AND. kt==nit000+1 ) ) THEN
315	ll_fw_start=.FALSE.
316	CALL ts_wgt(ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2)
317	ENDIF
318
319	! -----------------------------------------------------------------------------
320	! Phase 1 : Coupling between general trend and barotropic estimates (1st step)
321	! -----------------------------------------------------------------------------
322	!
323	!
324	! !* e3*d/dt(Ua) (Vertically integrated)
325	! ! --------------------------------------------------
326	zu_frc(:,:) = 0._wp
327	zv_frc(:,:) = 0._wp
328	!
329	DO jk = 1, jpkm1
330	zu_frc(:,:) = zu_frc(:,:) + e3u_n(:,:,jk) * ua(:,:,jk) * umask(:,:,jk)
331	zv_frc(:,:) = zv_frc(:,:) + e3v_n(:,:,jk) * va(:,:,jk) * vmask(:,:,jk)
332	END DO
333	!
334	zu_frc(:,:) = zu_frc(:,:) * r1_hu_n(:,:)
335	zv_frc(:,:) = zv_frc(:,:) * r1_hv_n(:,:)
336	!
337	!
338	! !* baroclinic momentum trend (remove the vertical mean trend)
339	DO jk = 1, jpkm1 ! -----------------------------------------------------------
340	DO jj = 2, jpjm1
341	DO ji = fs_2, fs_jpim1 ! vector opt.
342	ua(ji,jj,jk) = ua(ji,jj,jk) - zu_frc(ji,jj) * umask(ji,jj,jk)
343	va(ji,jj,jk) = va(ji,jj,jk) - zv_frc(ji,jj) * vmask(ji,jj,jk)
344	END DO
345	END DO
346	END DO
347
348	!!gm Question here when removing the Vertically integrated trends, we remove the vertically integrated NL trends on momentum....
349	!!gm Is it correct to do so ? I think so...
350
351
352	! !* barotropic Coriolis trends (vorticity scheme dependent)
353	! ! --------------------------------------------------------
354	zwx(:,:) = un_b(:,:) * hu_n(:,:) * e2u(:,:) ! now fluxes
355	zwy(:,:) = vn_b(:,:) * hv_n(:,:) * e1v(:,:)
356	!
357	IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN ! energy conserving or mixed scheme
358	DO jj = 2, jpjm1
359	DO ji = fs_2, fs_jpim1 ! vector opt.
360	zy1 = ( zwy(ji,jj-1) + zwy(ji+1,jj-1) ) * r1_e1u(ji,jj)
361	zy2 = ( zwy(ji,jj ) + zwy(ji+1,jj ) ) * r1_e1u(ji,jj)
362	zx1 = ( zwx(ji-1,jj) + zwx(ji-1,jj+1) ) * r1_e2v(ji,jj)
363	zx2 = ( zwx(ji ,jj) + zwx(ji ,jj+1) ) * r1_e2v(ji,jj)
364	! energy conserving formulation for planetary vorticity term
365	zu_trd(ji,jj) = z1_4 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
366	zv_trd(ji,jj) =-z1_4 * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 )
367	END DO
368	END DO
369	!
370	ELSEIF ( ln_dynvor_ens ) THEN ! enstrophy conserving scheme
371	DO jj = 2, jpjm1
372	DO ji = fs_2, fs_jpim1 ! vector opt.
373	zy1 = z1_8 * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) &
374	& + zwy(ji ,jj ) + zwy(ji+1,jj ) ) * r1_e1u(ji,jj)
375	zx1 = - z1_8 * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) &
376	& + zwx(ji ,jj ) + zwx(ji ,jj+1) ) * r1_e2v(ji,jj)
377	zu_trd(ji,jj) = zy1 * ( zwz(ji ,jj-1) + zwz(ji,jj) )
378	zv_trd(ji,jj) = zx1 * ( zwz(ji-1,jj ) + zwz(ji,jj) )
379	END DO
380	END DO
381	!
382	ELSEIF ( ln_dynvor_een ) THEN ! enstrophy and energy conserving scheme
383	DO jj = 2, jpjm1
384	DO ji = fs_2, fs_jpim1 ! vector opt.
385	zu_trd(ji,jj) = + z1_12 * r1_e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) &
386	& + ftnw(ji+1,jj) * zwy(ji+1,jj ) &
387	& + ftse(ji,jj ) * zwy(ji ,jj-1) &
388	& + ftsw(ji+1,jj) * zwy(ji+1,jj-1) )
389	zv_trd(ji,jj) = - z1_12 * r1_e2v(ji,jj) * ( ftsw(ji,jj+1) * zwx(ji-1,jj+1) &
390	& + ftse(ji,jj+1) * zwx(ji ,jj+1) &
391	& + ftnw(ji,jj ) * zwx(ji-1,jj ) &
392	& + ftne(ji,jj ) * zwx(ji ,jj ) )
393	END DO
394	END DO
395	!
396	ENDIF
397	!
398	! !* Right-Hand-Side of the barotropic momentum equation
399	! ! ----------------------------------------------------
400	IF( .NOT.ln_linssh ) THEN ! Variable volume : remove surface pressure gradient
401	IF( ln_wd ) THEN ! Calculating and applying W/D gravity filters
402	DO jj = 2, jpjm1
403	DO ji = 2, jpim1
404	ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji+1,jj) ) > &
405	& MAX( -ht_wd(ji,jj) , -ht_wd(ji+1,jj) ) .AND. &
406	& MAX( sshn(ji,jj) + ht_wd(ji,jj), sshn(ji+1,jj) + ht_wd(ji+1,jj) ) &
407	& > rn_wdmin1 + rn_wdmin2
408	ll_tmp2 = ( ABS( sshn(ji+1,jj) - sshn(ji ,jj)) > 1.E-12 ).AND.( &
409	& MAX( sshn(ji,jj) , sshn(ji+1,jj) ) > &
410	& MAX( -ht_wd(ji,jj) , -ht_wd(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 )
411
412	IF(ll_tmp1) THEN
413	zcpx(ji,jj) = 1.0_wp
414	ELSE IF(ll_tmp2) THEN
415	! no worries about sshn(ji+1,jj) - sshn(ji ,jj) = 0, it won't happen ! here
416	zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_wd(ji+1,jj) - sshn(ji,jj) - ht_wd(ji,jj)) &
417	& / (sshn(ji+1,jj) - sshn(ji ,jj)) )
418	ELSE
419	zcpx(ji,jj) = 0._wp
420	END IF
421
422	ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji,jj+1) ) > &
423	& MAX( -ht_wd(ji,jj) , -ht_wd(ji,jj+1) ) .AND. &
424	& MAX( sshn(ji,jj) + ht_wd(ji,jj), sshn(ji,jj+1) + ht_wd(ji,jj+1) ) &
425	& > rn_wdmin1 + rn_wdmin2
426	ll_tmp2 = ( ABS( sshn(ji,jj) - sshn(ji,jj+1)) > 1.E-12 ).AND.( &
427	& MAX( sshn(ji,jj) , sshn(ji,jj+1) ) > &
428	& MAX( -ht_wd(ji,jj) , -ht_wd(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 )
429
430	IF(ll_tmp1) THEN
431	zcpy(ji,jj) = 1.0_wp
432	ELSE IF(ll_tmp2) THEN
433	! no worries about sshn(ji,jj+1) - sshn(ji,jj ) = 0, it won't happen ! here
434	zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_wd(ji,jj+1) - sshn(ji,jj) - ht_wd(ji,jj)) &
435	& / (sshn(ji,jj+1) - sshn(ji,jj )) )
436	ELSE
437	zcpy(ji,jj) = 0._wp
438	END IF
439	END DO
440	END DO
441
442	DO jj = 2, jpjm1
443	DO ji = 2, jpim1
444	zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) &
445	& * r1_e1u(ji,jj) * zcpx(ji,jj)
446	zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji ,jj+1) - sshn(ji ,jj ) ) &
447	& * r1_e2v(ji,jj) * zcpy(ji,jj)
448	END DO
449	END DO
450
451	ELSE
452
453	DO jj = 2, jpjm1
454	DO ji = fs_2, fs_jpim1 ! vector opt.
455	zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) * r1_e1u(ji,jj)
456	zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji ,jj+1) - sshn(ji ,jj ) ) * r1_e2v(ji,jj)
457	END DO
458	END DO
459	ENDIF
460
461	ENDIF
462
463	DO jj = 2, jpjm1 ! Remove coriolis term (and possibly spg) from barotropic trend
464	DO ji = fs_2, fs_jpim1
465	zu_frc(ji,jj) = zu_frc(ji,jj) - zu_trd(ji,jj) * ssumask(ji,jj)
466	zv_frc(ji,jj) = zv_frc(ji,jj) - zv_trd(ji,jj) * ssvmask(ji,jj)
467	END DO
468	END DO
469	!
470	! ! Add bottom stress contribution from baroclinic velocities:
471	IF (ln_bt_fw) THEN
472	DO jj = 2, jpjm1
473	DO ji = fs_2, fs_jpim1 ! vector opt.
474	ikbu = mbku(ji,jj)
475	ikbv = mbkv(ji,jj)
476	zwx(ji,jj) = un(ji,jj,ikbu) - un_b(ji,jj) ! NOW bottom baroclinic velocities
477	zwy(ji,jj) = vn(ji,jj,ikbv) - vn_b(ji,jj)
478	END DO
479	END DO
480	ELSE
481	DO jj = 2, jpjm1
482	DO ji = fs_2, fs_jpim1 ! vector opt.
483	ikbu = mbku(ji,jj)
484	ikbv = mbkv(ji,jj)
485	zwx(ji,jj) = ub(ji,jj,ikbu) - ub_b(ji,jj) ! BEFORE bottom baroclinic velocities
486	zwy(ji,jj) = vb(ji,jj,ikbv) - vb_b(ji,jj)
487	END DO
488	END DO
489	ENDIF
490	!
491	! Note that the "unclipped" bottom friction parameter is used even with explicit drag
492	IF( ln_wd ) THEN
493	zu_frc(:,:) = zu_frc(:,:) + MAX(r1_hu_n(:,:) * bfrua(:,:),-1._wp / rdtbt) * zwx(:,:)
494	zv_frc(:,:) = zv_frc(:,:) + MAX(r1_hv_n(:,:) * bfrva(:,:),-1._wp / rdtbt) * zwy(:,:)
495	ELSE
496	zu_frc(:,:) = zu_frc(:,:) + r1_hu_n(:,:) * bfrua(:,:) * zwx(:,:)
497	zv_frc(:,:) = zv_frc(:,:) + r1_hv_n(:,:) * bfrva(:,:) * zwy(:,:)
498	END IF
499	!
500	! ! Add top stress contribution from baroclinic velocities:
501	IF( ln_bt_fw ) THEN
502	DO jj = 2, jpjm1
503	DO ji = fs_2, fs_jpim1 ! vector opt.
504	iktu = miku(ji,jj)
505	iktv = mikv(ji,jj)
506	zwx(ji,jj) = un(ji,jj,iktu) - un_b(ji,jj) ! NOW top baroclinic velocities
507	zwy(ji,jj) = vn(ji,jj,iktv) - vn_b(ji,jj)
508	END DO
509	END DO
510	ELSE
511	DO jj = 2, jpjm1
512	DO ji = fs_2, fs_jpim1 ! vector opt.
513	iktu = miku(ji,jj)
514	iktv = mikv(ji,jj)
515	zwx(ji,jj) = ub(ji,jj,iktu) - ub_b(ji,jj) ! BEFORE top baroclinic velocities
516	zwy(ji,jj) = vb(ji,jj,iktv) - vb_b(ji,jj)
517	END DO
518	END DO
519	ENDIF
520	!
521	! Note that the "unclipped" top friction parameter is used even with explicit drag
522	zu_frc(:,:) = zu_frc(:,:) + r1_hu_n(:,:) * tfrua(:,:) * zwx(:,:)
523	zv_frc(:,:) = zv_frc(:,:) + r1_hv_n(:,:) * tfrva(:,:) * zwy(:,:)
524	!
525	IF (ln_bt_fw) THEN ! Add wind forcing
526	zu_frc(:,:) = zu_frc(:,:) + zraur * utau(:,:) * r1_hu_n(:,:)
527	zv_frc(:,:) = zv_frc(:,:) + zraur * vtau(:,:) * r1_hv_n(:,:)
528	ELSE
529	zu_frc(:,:) = zu_frc(:,:) + zraur * z1_2 * ( utau_b(:,:) + utau(:,:) ) * r1_hu_n(:,:)
530	zv_frc(:,:) = zv_frc(:,:) + zraur * z1_2 * ( vtau_b(:,:) + vtau(:,:) ) * r1_hv_n(:,:)
531	ENDIF
532	!
533	IF ( ln_apr_dyn ) THEN ! Add atm pressure forcing
534	IF (ln_bt_fw) THEN
535	DO jj = 2, jpjm1
536	DO ji = fs_2, fs_jpim1 ! vector opt.
537	zu_spg = grav * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) ) * r1_e1u(ji,jj)
538	zv_spg = grav * ( ssh_ib (ji ,jj+1) - ssh_ib (ji,jj) ) * r1_e2v(ji,jj)
539	zu_frc(ji,jj) = zu_frc(ji,jj) + zu_spg
540	zv_frc(ji,jj) = zv_frc(ji,jj) + zv_spg
541	END DO
542	END DO
543	ELSE
544	DO jj = 2, jpjm1
545	DO ji = fs_2, fs_jpim1 ! vector opt.
546	zu_spg = grav * z1_2 * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) &
547	& + ssh_ibb(ji+1,jj ) - ssh_ibb(ji,jj) ) * r1_e1u(ji,jj)
548	zv_spg = grav * z1_2 * ( ssh_ib (ji ,jj+1) - ssh_ib (ji,jj) &
549	& + ssh_ibb(ji ,jj+1) - ssh_ibb(ji,jj) ) * r1_e2v(ji,jj)
550	zu_frc(ji,jj) = zu_frc(ji,jj) + zu_spg
551	zv_frc(ji,jj) = zv_frc(ji,jj) + zv_spg
552	END DO
553	END DO
554	ENDIF
555	ENDIF
556	! !* Right-Hand-Side of the barotropic ssh equation
557	! ! -----------------------------------------------
558	! ! Surface net water flux and rivers
559	IF (ln_bt_fw) THEN
560	zssh_frc(:,:) = zraur * ( emp(:,:) - rnf(:,:) + fwfisf(:,:) )
561	ELSE
562	zssh_frc(:,:) = zraur * z1_2 * ( emp(:,:) + emp_b(:,:) - rnf(:,:) - rnf_b(:,:) &
563	& + fwfisf(:,:) + fwfisf_b(:,:) )
564	ENDIF
565	!
566	IF( ln_sdw ) THEN ! Stokes drift divergence added if necessary
567	zssh_frc(:,:) = zssh_frc(:,:) + div_sd(:,:)
568	ENDIF
569	!
570	#if defined key_asminc
571	! ! Include the IAU weighted SSH increment
572	IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN
573	zssh_frc(:,:) = zssh_frc(:,:) - ssh_iau(:,:)
574	ENDIF
575	#endif
576	! !* Fill boundary data arrays for AGRIF
577	! ! ------------------------------------
578	#if defined key_agrif
579	IF( .NOT.Agrif_Root() ) CALL agrif_dta_ts( kt )
580	#endif
581	!
582	! -----------------------------------------------------------------------
583	! Phase 2 : Integration of the barotropic equations
584	! -----------------------------------------------------------------------
585	!
586	! ! ==================== !
587	! ! Initialisations !
588	! ! ==================== !
589	! Initialize barotropic variables:
590	IF( ll_init )THEN
591	sshbb_e(:,:) = 0._wp
592	ubb_e (:,:) = 0._wp
593	vbb_e (:,:) = 0._wp
594	sshb_e (:,:) = 0._wp
595	ub_e (:,:) = 0._wp
596	vb_e (:,:) = 0._wp
597	ENDIF
598
599	!
600	IF (ln_bt_fw) THEN ! FORWARD integration: start from NOW fields
601	sshn_e(:,:) = sshn(:,:)
602	un_e (:,:) = un_b(:,:)
603	vn_e (:,:) = vn_b(:,:)
604	!
605	hu_e (:,:) = hu_n(:,:)
606	hv_e (:,:) = hv_n(:,:)
607	hur_e (:,:) = r1_hu_n(:,:)
608	hvr_e (:,:) = r1_hv_n(:,:)
609	ELSE ! CENTRED integration: start from BEFORE fields
610	sshn_e(:,:) = sshb(:,:)
611	un_e (:,:) = ub_b(:,:)
612	vn_e (:,:) = vb_b(:,:)
613	!
614	hu_e (:,:) = hu_b(:,:)
615	hv_e (:,:) = hv_b(:,:)
616	hur_e (:,:) = r1_hu_b(:,:)
617	hvr_e (:,:) = r1_hv_b(:,:)
618	ENDIF
619	!
620	!
621	!
622	! Initialize sums:
623	ua_b (:,:) = 0._wp ! After barotropic velocities (or transport if flux form)
624	va_b (:,:) = 0._wp
625	ssha (:,:) = 0._wp ! Sum for after averaged sea level
626	un_adv(:,:) = 0._wp ! Sum for now transport issued from ts loop
627	vn_adv(:,:) = 0._wp
628	! ! ==================== !
629	DO jn = 1, icycle ! sub-time-step loop !
630	! ! ==================== !
631	! !* Update the forcing (BDY and tides)
632	! ! ------------------
633	! Update only tidal forcing at open boundaries
634	IF( ln_bdy .AND. ln_tide ) CALL bdy_dta_tides( kt, kit=jn, time_offset= noffset+1 )
635	IF( ln_tide_pot .AND. ln_tide ) CALL upd_tide ( kt, kit=jn, time_offset= noffset )
636	!
637	! Set extrapolation coefficients for predictor step:
638	IF ((jn<3).AND.ll_init) THEN ! Forward
639	za1 = 1._wp
640	za2 = 0._wp
641	za3 = 0._wp
642	ELSE ! AB3-AM4 Coefficients: bet=0.281105
643	za1 = 1.781105_wp ! za1 = 3/2 + bet
644	za2 = -1.06221_wp ! za2 = -(1/2 + 2*bet)
645	za3 = 0.281105_wp ! za3 = bet
646	ENDIF
647
648	! Extrapolate barotropic velocities at step jit+0.5:
649	ua_e(:,:) = za1 * un_e(:,:) + za2 * ub_e(:,:) + za3 * ubb_e(:,:)
650	va_e(:,:) = za1 * vn_e(:,:) + za2 * vb_e(:,:) + za3 * vbb_e(:,:)
651
652	IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only)
653	! ! ------------------
654	! Extrapolate Sea Level at step jit+0.5:
655	zsshp2_e(:,:) = za1 * sshn_e(:,:) + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:)
656	!
657	DO jj = 2, jpjm1 ! Sea Surface Height at u- & v-points
658	DO ji = 2, fs_jpim1 ! Vector opt.
659	zwx(ji,jj) = z1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj) &
660	& * ( e1e2t(ji ,jj) * zsshp2_e(ji ,jj) &
661	& + e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj) )
662	zwy(ji,jj) = z1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj) &
663	& * ( e1e2t(ji,jj ) * zsshp2_e(ji,jj ) &
664	& + e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1) )
665	END DO
666	END DO
667	CALL lbc_lnk_multi( zwx, 'U', 1._wp, zwy, 'V', 1._wp )
668	!
669	zhup2_e (:,:) = hu_0(:,:) + zwx(:,:) ! Ocean depth at U- and V-points
670	zhvp2_e (:,:) = hv_0(:,:) + zwy(:,:)
671	ELSE
672	zhup2_e (:,:) = hu_n(:,:)
673	zhvp2_e (:,:) = hv_n(:,:)
674	ENDIF
675	! !* after ssh
676	! ! -----------
677	! One should enforce volume conservation at open boundaries here
678	! considering fluxes below:
679	!
680	zwx(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:) ! fluxes at jn+0.5
681	zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:)
682	!
683	#if defined key_agrif
684	! Set fluxes during predictor step to ensure volume conservation
685	IF( .NOT.Agrif_Root() .AND. ln_bt_fw ) THEN
686	IF((nbondi == -1).OR.(nbondi == 2)) THEN
687	DO jj=1,jpj
688	zwx(2,jj) = ubdy_w(jj) * e2u(2,jj)
689	END DO
690	ENDIF
691	IF((nbondi == 1).OR.(nbondi == 2)) THEN
692	DO jj=1,jpj
693	zwx(nlci-2,jj) = ubdy_e(jj) * e2u(nlci-2,jj)
694	END DO
695	ENDIF
696	IF((nbondj == -1).OR.(nbondj == 2)) THEN
697	DO ji=1,jpi
698	zwy(ji,2) = vbdy_s(ji) * e1v(ji,2)
699	END DO
700	ENDIF
701	IF((nbondj == 1).OR.(nbondj == 2)) THEN
702	DO ji=1,jpi
703	zwy(ji,nlcj-2) = vbdy_n(ji) * e1v(ji,nlcj-2)
704	END DO
705	ENDIF
706	ENDIF
707	#endif
708	IF( ln_wd ) CALL wad_lmt_bt(zwx, zwy, sshn_e, zssh_frc, rdtbt)
709	!
710	! Sum over sub-time-steps to compute advective velocities
711	za2 = wgtbtp2(jn)
712	un_adv(:,:) = un_adv(:,:) + za2 * zwx(:,:) * r1_e2u(:,:)
713	vn_adv(:,:) = vn_adv(:,:) + za2 * zwy(:,:) * r1_e1v(:,:)
714	!
715	! Set next sea level:
716	DO jj = 2, jpjm1
717	DO ji = fs_2, fs_jpim1 ! vector opt.
718	zhdiv(ji,jj) = ( zwx(ji,jj) - zwx(ji-1,jj) &
719	& + zwy(ji,jj) - zwy(ji,jj-1) ) * r1_e1e2t(ji,jj)
720	END DO
721	END DO
722	ssha_e(:,:) = ( sshn_e(:,:) - rdtbt * ( zssh_frc(:,:) + zhdiv(:,:) ) ) * ssmask(:,:)
723
724	CALL lbc_lnk( ssha_e, 'T', 1._wp )
725
726	! Duplicate sea level across open boundaries (this is only cosmetic if linssh=T)
727	IF( ln_bdy ) CALL bdy_ssh( ssha_e )
728	#if defined key_agrif
729	IF( .NOT.Agrif_Root() ) CALL agrif_ssh_ts( jn )
730	#endif
731	!
732	! Sea Surface Height at u-,v-points (vvl case only)
733	IF( .NOT.ln_linssh ) THEN
734	DO jj = 2, jpjm1
735	DO ji = 2, jpim1 ! NO Vector Opt.
736	zsshu_a(ji,jj) = z1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj) &
737	& * ( e1e2t(ji ,jj ) * ssha_e(ji ,jj ) &
738	& + e1e2t(ji+1,jj ) * ssha_e(ji+1,jj ) )
739	zsshv_a(ji,jj) = z1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj) &
740	& * ( e1e2t(ji ,jj ) * ssha_e(ji ,jj ) &
741	& + e1e2t(ji ,jj+1) * ssha_e(ji ,jj+1) )
742	END DO
743	END DO
744	CALL lbc_lnk_multi( zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp )
745	ENDIF
746	!
747	! Half-step back interpolation of SSH for surface pressure computation:
748	!----------------------------------------------------------------------
749	IF ((jn==1).AND.ll_init) THEN
750	za0=1._wp ! Forward-backward
751	za1=0._wp
752	za2=0._wp
753	za3=0._wp
754	ELSEIF ((jn==2).AND.ll_init) THEN ! AB2-AM3 Coefficients; bet=0 ; gam=-1/6 ; eps=1/12
755	za0= 1.0833333333333_wp ! za0 = 1-gam-eps
756	za1=-0.1666666666666_wp ! za1 = gam
757	za2= 0.0833333333333_wp ! za2 = eps
758	za3= 0._wp
759	ELSE ! AB3-AM4 Coefficients; bet=0.281105 ; eps=0.013 ; gam=0.0880
760	za0=0.614_wp ! za0 = 1/2 + gam + 2*eps
761	za1=0.285_wp ! za1 = 1/2 - 2gam - 3eps
762	za2=0.088_wp ! za2 = gam
763	za3=0.013_wp ! za3 = eps
764	ENDIF
765	!
766	zsshp2_e(:,:) = za0 * ssha_e(:,:) + za1 * sshn_e (:,:) &
767	& + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:)
768	IF( ln_wd ) THEN ! Calculating and applying W/D gravity filters
769	DO jj = 2, jpjm1
770	DO ji = 2, jpim1
771	ll_tmp1 = MIN( zsshp2_e(ji,jj) , zsshp2_e(ji+1,jj) ) > &
772	& MAX( -ht_wd(ji,jj) , -ht_wd(ji+1,jj) ) .AND. &
773	& MAX( zsshp2_e(ji,jj) + ht_wd(ji,jj), zsshp2_e(ji+1,jj) + ht_wd(ji+1,jj) ) &
774	& > rn_wdmin1 + rn_wdmin2
775	ll_tmp2 = (ABS(zsshp2_e(ji,jj) - zsshp2_e(ji+1,jj)) > 1.E-12 ).AND.( &
776	& MAX( zsshp2_e(ji,jj) , zsshp2_e(ji+1,jj) ) > &
777	& MAX( -ht_wd(ji,jj) , -ht_wd(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 )
778
779	IF(ll_tmp1) THEN
780	zcpx(ji,jj) = 1.0_wp
781	ELSE IF(ll_tmp2) THEN
782	! no worries about zsshp2_e(ji+1,jj) - zsshp2_e(ji ,jj) = 0, it won't happen ! here
783	zcpx(ji,jj) = ABS( (zsshp2_e(ji+1,jj) + ht_wd(ji+1,jj) - zsshp2_e(ji,jj) - ht_wd(ji,jj)) &
784	& / (zsshp2_e(ji+1,jj) - zsshp2_e(ji ,jj)) )
785	ELSE
786	zcpx(ji,jj) = 0._wp
787	END IF
788
789	ll_tmp1 = MIN( zsshp2_e(ji,jj) , zsshp2_e(ji,jj+1) ) > &
790	& MAX( -ht_wd(ji,jj) , -ht_wd(ji,jj+1) ) .AND. &
791	& MAX( zsshp2_e(ji,jj) + ht_wd(ji,jj), zsshp2_e(ji,jj+1) + ht_wd(ji,jj+1) ) &
792	& > rn_wdmin1 + rn_wdmin2
793	ll_tmp2 = (ABS(zsshp2_e(ji,jj) - zsshp2_e(ji,jj+1)) > 1.E-12 ).AND.( &
794	& MAX( zsshp2_e(ji,jj) , zsshp2_e(ji,jj+1) ) > &
795	& MAX( -ht_wd(ji,jj) , -ht_wd(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 )
796
797	IF(ll_tmp1) THEN
798	zcpy(ji,jj) = 1.0_wp
799	ELSE IF(ll_tmp2) THEN
800	! no worries about zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj ) = 0, it won't happen ! here
801	zcpy(ji,jj) = ABS( (zsshp2_e(ji,jj+1) + ht_wd(ji,jj+1) - zsshp2_e(ji,jj) - ht_wd(ji,jj)) &
802	& / (zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj )) )
803	ELSE
804	zcpy(ji,jj) = 0._wp
805	END IF
806	END DO
807	END DO
808	END IF
809	!
810	! Compute associated depths at U and V points:
811	IF( .NOT.ln_linssh .AND. .NOT.ln_dynadv_vec ) THEN !* Vector form
812	!
813	DO jj = 2, jpjm1
814	DO ji = 2, jpim1
815	zx1 = z1_2 * ssumask(ji ,jj) * r1_e1e2u(ji ,jj) &
816	& * ( e1e2t(ji ,jj ) * zsshp2_e(ji ,jj) &
817	& + e1e2t(ji+1,jj ) * zsshp2_e(ji+1,jj ) )
818	zy1 = z1_2 * ssvmask(ji ,jj) * r1_e1e2v(ji ,jj ) &
819	& * ( e1e2t(ji ,jj ) * zsshp2_e(ji ,jj ) &
820	& + e1e2t(ji ,jj+1) * zsshp2_e(ji ,jj+1) )
821	zhust_e(ji,jj) = hu_0(ji,jj) + zx1
822	zhvst_e(ji,jj) = hv_0(ji,jj) + zy1
823	END DO
824	END DO
825
826	ENDIF
827	!
828	! Add Coriolis trend:
829	! zwz array below or triads normally depend on sea level with ln_linssh=F and should be updated
830	! at each time step. We however keep them constant here for optimization.
831	! Recall that zwx and zwy arrays hold fluxes at this stage:
832	! zwx(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:) ! fluxes at jn+0.5
833	! zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:)
834	!
835	IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN !== energy conserving or mixed scheme ==!
836	DO jj = 2, jpjm1
837	DO ji = fs_2, fs_jpim1 ! vector opt.
838	zy1 = ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) ) * r1_e1u(ji,jj)
839	zy2 = ( zwy(ji ,jj ) + zwy(ji+1,jj ) ) * r1_e1u(ji,jj)
840	zx1 = ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) ) * r1_e2v(ji,jj)
841	zx2 = ( zwx(ji ,jj ) + zwx(ji ,jj+1) ) * r1_e2v(ji,jj)
842	zu_trd(ji,jj) = z1_4 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
843	zv_trd(ji,jj) =-z1_4 * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 )
844	END DO
845	END DO
846	!
847	ELSEIF ( ln_dynvor_ens ) THEN !== enstrophy conserving scheme ==!
848	DO jj = 2, jpjm1
849	DO ji = fs_2, fs_jpim1 ! vector opt.
850	zy1 = z1_8 * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) &
851	& + zwy(ji ,jj ) + zwy(ji+1,jj ) ) * r1_e1u(ji,jj)
852	zx1 = - z1_8 * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) &
853	& + zwx(ji ,jj ) + zwx(ji ,jj+1) ) * r1_e2v(ji,jj)
854	zu_trd(ji,jj) = zy1 * ( zwz(ji ,jj-1) + zwz(ji,jj) )
855	zv_trd(ji,jj) = zx1 * ( zwz(ji-1,jj ) + zwz(ji,jj) )
856	END DO
857	END DO
858	!
859	ELSEIF ( ln_dynvor_een ) THEN !== energy and enstrophy conserving scheme ==!
860	DO jj = 2, jpjm1
861	DO ji = fs_2, fs_jpim1 ! vector opt.
862	zu_trd(ji,jj) = + z1_12 * r1_e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) &
863	& + ftnw(ji+1,jj) * zwy(ji+1,jj ) &
864	& + ftse(ji,jj ) * zwy(ji ,jj-1) &
865	& + ftsw(ji+1,jj) * zwy(ji+1,jj-1) )
866	zv_trd(ji,jj) = - z1_12 * r1_e2v(ji,jj) * ( ftsw(ji,jj+1) * zwx(ji-1,jj+1) &
867	& + ftse(ji,jj+1) * zwx(ji ,jj+1) &
868	& + ftnw(ji,jj ) * zwx(ji-1,jj ) &
869	& + ftne(ji,jj ) * zwx(ji ,jj ) )
870	END DO
871	END DO
872	!
873	ENDIF
874	!
875	! Add tidal astronomical forcing if defined
876	IF ( ln_tide .AND. ln_tide_pot ) THEN
877	DO jj = 2, jpjm1
878	DO ji = fs_2, fs_jpim1 ! vector opt.
879	zu_spg = grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) * r1_e1u(ji,jj)
880	zv_spg = grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) * r1_e2v(ji,jj)
881	zu_trd(ji,jj) = zu_trd(ji,jj) + zu_spg
882	zv_trd(ji,jj) = zv_trd(ji,jj) + zv_spg
883	END DO
884	END DO
885	ENDIF
886	!
887	! Add bottom stresses:
888	zu_trd(:,:) = zu_trd(:,:) + bfrua(:,:) * un_e(:,:) * hur_e(:,:)
889	zv_trd(:,:) = zv_trd(:,:) + bfrva(:,:) * vn_e(:,:) * hvr_e(:,:)
890	!
891	! Add top stresses:
892	zu_trd(:,:) = zu_trd(:,:) + tfrua(:,:) * un_e(:,:) * hur_e(:,:)
893	zv_trd(:,:) = zv_trd(:,:) + tfrva(:,:) * vn_e(:,:) * hvr_e(:,:)
894	!
895	! Surface pressure trend:
896
897	IF( ln_wd ) THEN
898	DO jj = 2, jpjm1
899	DO ji = 2, jpim1
900	! Add surface pressure gradient
901	zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj)
902	zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj)
903	zwx(ji,jj) = zu_spg * zcpx(ji,jj)
904	zwy(ji,jj) = zv_spg * zcpy(ji,jj)
905	END DO
906	END DO
907	ELSE
908	DO jj = 2, jpjm1
909	DO ji = fs_2, fs_jpim1 ! vector opt.
910	! Add surface pressure gradient
911	zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj)
912	zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj)
913	zwx(ji,jj) = zu_spg
914	zwy(ji,jj) = zv_spg
915	END DO
916	END DO
917	END IF
918
919	!
920	! Set next velocities:
921	IF( ln_dynadv_vec .OR. ln_linssh ) THEN !* Vector form
922	DO jj = 2, jpjm1
923	DO ji = fs_2, fs_jpim1 ! vector opt.
924	ua_e(ji,jj) = ( un_e(ji,jj) &
925	& + rdtbt * ( zwx(ji,jj) &
926	& + zu_trd(ji,jj) &
927	& + zu_frc(ji,jj) ) &
928	& ) * ssumask(ji,jj)
929
930	va_e(ji,jj) = ( vn_e(ji,jj) &
931	& + rdtbt * ( zwy(ji,jj) &
932	& + zv_trd(ji,jj) &
933	& + zv_frc(ji,jj) ) &
934	& ) * ssvmask(ji,jj)
935	END DO
936	END DO
937	!
938	ELSE !* Flux form
939	DO jj = 2, jpjm1
940	DO ji = fs_2, fs_jpim1 ! vector opt.
941
942	IF( ln_wd ) THEN
943	zhura = MAX(hu_0(ji,jj) + zsshu_a(ji,jj), rn_wdmin1)
944	zhvra = MAX(hv_0(ji,jj) + zsshv_a(ji,jj), rn_wdmin1)
945	ELSE
946	zhura = hu_0(ji,jj) + zsshu_a(ji,jj)
947	zhvra = hv_0(ji,jj) + zsshv_a(ji,jj)
948	END IF
949	zhura = ssumask(ji,jj)/(zhura + 1._wp - ssumask(ji,jj))
950	zhvra = ssvmask(ji,jj)/(zhvra + 1._wp - ssvmask(ji,jj))
951
952	ua_e(ji,jj) = ( hu_e(ji,jj) * un_e(ji,jj) &
953	& + rdtbt * ( zhust_e(ji,jj) * zwx(ji,jj) &
954	& + zhup2_e(ji,jj) * zu_trd(ji,jj) &
955	& + hu_n(ji,jj) * zu_frc(ji,jj) ) &
956	& ) * zhura
957
958	va_e(ji,jj) = ( hv_e(ji,jj) * vn_e(ji,jj) &
959	& + rdtbt * ( zhvst_e(ji,jj) * zwy(ji,jj) &
960	& + zhvp2_e(ji,jj) * zv_trd(ji,jj) &
961	& + hv_n(ji,jj) * zv_frc(ji,jj) ) &
962	& ) * zhvra
963	END DO
964	END DO
965	ENDIF
966	!
967	IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only)
968	IF( ln_wd ) THEN
969	hu_e (:,:) = MAX(hu_0(:,:) + zsshu_a(:,:), rn_wdmin1)
970	hv_e (:,:) = MAX(hv_0(:,:) + zsshv_a(:,:), rn_wdmin1)
971	ELSE
972	hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:)
973	hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:)
974	END IF
975	hur_e(:,:) = ssumask(:,:) / ( hu_e(:,:) + 1._wp - ssumask(:,:) )
976	hvr_e(:,:) = ssvmask(:,:) / ( hv_e(:,:) + 1._wp - ssvmask(:,:) )
977	!
978	ENDIF
979	! !* domain lateral boundary
980	CALL lbc_lnk_multi( ua_e, 'U', -1._wp, va_e , 'V', -1._wp )
981	!
982	! ! open boundaries
983	IF( ln_bdy ) CALL bdy_dyn2d( jn, ua_e, va_e, un_e, vn_e, hur_e, hvr_e, ssha_e )
984	#if defined key_agrif
985	IF( .NOT.Agrif_Root() ) CALL agrif_dyn_ts( jn ) ! Agrif
986	#endif
987	! !* Swap
988	! ! ----
989	ubb_e (:,:) = ub_e (:,:)
990	ub_e (:,:) = un_e (:,:)
991	un_e (:,:) = ua_e (:,:)
992	!
993	vbb_e (:,:) = vb_e (:,:)
994	vb_e (:,:) = vn_e (:,:)
995	vn_e (:,:) = va_e (:,:)
996	!
997	sshbb_e(:,:) = sshb_e(:,:)
998	sshb_e (:,:) = sshn_e(:,:)
999	sshn_e (:,:) = ssha_e(:,:)
1000
1001	! !* Sum over whole bt loop
1002	! ! ----------------------
1003	za1 = wgtbtp1(jn)
1004	IF( ln_dynadv_vec .OR. ln_linssh ) THEN ! Sum velocities
1005	ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:)
1006	va_b (:,:) = va_b (:,:) + za1 * va_e (:,:)
1007	ELSE ! Sum transports
1008	ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:) * hu_e (:,:)
1009	va_b (:,:) = va_b (:,:) + za1 * va_e (:,:) * hv_e (:,:)
1010	ENDIF
1011	! ! Sum sea level
1012	ssha(:,:) = ssha(:,:) + za1 * ssha_e(:,:)
1013	! ! ==================== !
1014	END DO ! end loop !
1015	! ! ==================== !
1016	! -----------------------------------------------------------------------------
1017	! Phase 3. update the general trend with the barotropic trend
1018	! -----------------------------------------------------------------------------
1019	!
1020	! Set advection velocity correction:
1021	zwx(:,:) = un_adv(:,:)
1022	zwy(:,:) = vn_adv(:,:)
1023	IF( ( kt == nit000 .AND. neuler==0 ) .OR. .NOT.ln_bt_fw ) THEN
1024	un_adv(:,:) = zwx(:,:) * r1_hu_n(:,:)
1025	vn_adv(:,:) = zwy(:,:) * r1_hv_n(:,:)
1026	!
1027	un_bf(:,:) = 0._wp
1028	vn_bf(:,:) = 0._wp
1029	ELSE
1030	un_adv(:,:) = z1_2 * ( ub2_b(:,:) + zwx(:,:) - atfp * un_bf(:,:) ) * r1_hu_n(:,:)
1031	vn_adv(:,:) = z1_2 * ( vb2_b(:,:) + zwy(:,:) - atfp * vn_bf(:,:) ) * r1_hv_n(:,:)
1032	! Update corrective fluxes for next time step:
1033	un_bf(:,:) = atfp * un_bf(:,:) + (zwx(:,:) - ub2_b(:,:))
1034	vn_bf(:,:) = atfp * vn_bf(:,:) + (zwy(:,:) - vb2_b(:,:))
1035	END IF
1036
1037	IF( ln_bt_fw ) THEN ! Save integrated transport for next computation
1038	ub2_b(:,:) = zwx(:,:)
1039	vb2_b(:,:) = zwy(:,:)
1040	ENDIF
1041	!
1042	! Update barotropic trend:
1043	IF( ln_dynadv_vec .OR. ln_linssh ) THEN
1044	DO jk=1,jpkm1
1045	ua(:,:,jk) = ua(:,:,jk) + ( ua_b(:,:) - ub_b(:,:) ) * z1_2dt_b
1046	va(:,:,jk) = va(:,:,jk) + ( va_b(:,:) - vb_b(:,:) ) * z1_2dt_b
1047	END DO
1048	ELSE
1049	! At this stage, ssha has been corrected: compute new depths at velocity points
1050	DO jj = 1, jpjm1
1051	DO ji = 1, jpim1 ! NO Vector Opt.
1052	zsshu_a(ji,jj) = z1_2 * umask(ji,jj,1) * r1_e1e2u(ji,jj) &
1053	& * ( e1e2t(ji ,jj) * ssha(ji ,jj) &
1054	& + e1e2t(ji+1,jj) * ssha(ji+1,jj) )
1055	zsshv_a(ji,jj) = z1_2 * vmask(ji,jj,1) * r1_e1e2v(ji,jj) &
1056	& * ( e1e2t(ji,jj ) * ssha(ji,jj ) &
1057	& + e1e2t(ji,jj+1) * ssha(ji,jj+1) )
1058	END DO
1059	END DO
1060	CALL lbc_lnk_multi( zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) ! Boundary conditions
1061	!
1062	DO jk=1,jpkm1
1063	ua(:,:,jk) = ua(:,:,jk) + r1_hu_n(:,:) * ( ua_b(:,:) - ub_b(:,:) * hu_b(:,:) ) * z1_2dt_b
1064	va(:,:,jk) = va(:,:,jk) + r1_hv_n(:,:) * ( va_b(:,:) - vb_b(:,:) * hv_b(:,:) ) * z1_2dt_b
1065	END DO
1066	! Save barotropic velocities not transport:
1067	ua_b(:,:) = ua_b(:,:) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - ssumask(:,:) )
1068	va_b(:,:) = va_b(:,:) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - ssvmask(:,:) )
1069	ENDIF
1070	!
1071	DO jk = 1, jpkm1
1072	! Correct velocities:
1073	un(:,:,jk) = ( un(:,:,jk) + un_adv(:,:) - un_b(:,:) ) * umask(:,:,jk)
1074	vn(:,:,jk) = ( vn(:,:,jk) + vn_adv(:,:) - vn_b(:,:) ) * vmask(:,:,jk)
1075	!
1076	END DO
1077	!
1078	CALL iom_put( "ubar", un_adv(:,:) ) ! barotropic i-current
1079	CALL iom_put( "vbar", vn_adv(:,:) ) ! barotropic i-current
1080	!
1081	#if defined key_agrif
1082	! Save time integrated fluxes during child grid integration
1083	! (used to update coarse grid transports at next time step)
1084	!
1085	IF( .NOT.Agrif_Root() .AND. ln_bt_fw ) THEN
1086	IF( Agrif_NbStepint() == 0 ) THEN
1087	ub2_i_b(:,:) = 0._wp
1088	vb2_i_b(:,:) = 0._wp
1089	END IF
1090	!
1091	za1 = 1._wp / REAL(Agrif_rhot(), wp)
1092	ub2_i_b(:,:) = ub2_i_b(:,:) + za1 * ub2_b(:,:)
1093	vb2_i_b(:,:) = vb2_i_b(:,:) + za1 * vb2_b(:,:)
1094	ENDIF
1095	#endif
1096	! !* write time-spliting arrays in the restart
1097	IF( lrst_oce .AND.ln_bt_fw ) CALL ts_rst( kt, 'WRITE' )
1098	!
1099	CALL wrk_dealloc( jpi,jpj, zsshp2_e, zhdiv )
1100	CALL wrk_dealloc( jpi,jpj, zu_trd, zv_trd )
1101	CALL wrk_dealloc( jpi,jpj, zwx, zwy, zssh_frc, zu_frc, zv_frc )
1102	CALL wrk_dealloc( jpi,jpj, zhup2_e, zhvp2_e, zhust_e, zhvst_e )
1103	CALL wrk_dealloc( jpi,jpj, zsshu_a, zsshv_a )
1104	CALL wrk_dealloc( jpi,jpj, zhf )
1105	IF( ln_wd ) CALL wrk_dealloc( jpi, jpj, zcpx, zcpy )
1106	!
1107	IF ( ln_diatmb ) THEN
1108	CALL iom_put( "baro_u" , un_bumask(:,:,1)+zmdi(1-umask(:,:,1 ) ) ) ! Barotropic U Velocity
1109	CALL iom_put( "baro_v" , vn_bvmask(:,:,1)+zmdi(1-vmask(:,:,1 ) ) ) ! Barotropic V Velocity
1110	ENDIF
1111	IF( nn_timing == 1 ) CALL timing_stop('dyn_spg_ts')
1112	!
1113	END SUBROUTINE dyn_spg_ts
1114
1115
1116	SUBROUTINE ts_wgt( ll_av, ll_fw, jpit, zwgt1, zwgt2)
1117	!!---------------------------------------------------------------------
1118	!! * ROUTINE ts_wgt *
1119	!!
1120	!! ** Purpose : Set time-splitting weights for temporal averaging (or not)
1121	!!----------------------------------------------------------------------
1122	LOGICAL, INTENT(in) :: ll_av ! temporal averaging=.true.
1123	LOGICAL, INTENT(in) :: ll_fw ! forward time splitting =.true.
1124	INTEGER, INTENT(inout) :: jpit ! cycle length
1125	REAL(wp), DIMENSION(3*nn_baro), INTENT(inout) :: zwgt1, & ! Primary weights
1126	zwgt2 ! Secondary weights
1127
1128	INTEGER :: jic, jn, ji ! temporary integers
1129	REAL(wp) :: za1, za2
1130	!!----------------------------------------------------------------------
1131
1132	zwgt1(:) = 0._wp
1133	zwgt2(:) = 0._wp
1134
1135	! Set time index when averaged value is requested
1136	IF (ll_fw) THEN
1137	jic = nn_baro
1138	ELSE
1139	jic = 2 * nn_baro
1140	ENDIF
1141
1142	! Set primary weights:
1143	IF (ll_av) THEN
1144	! Define simple boxcar window for primary weights
1145	! (width = nn_baro, centered around jic)
1146	SELECT CASE ( nn_bt_flt )
1147	CASE( 0 ) ! No averaging
1148	zwgt1(jic) = 1._wp
1149	jpit = jic
1150
1151	CASE( 1 ) ! Boxcar, width = nn_baro
1152	DO jn = 1, 3*nn_baro
1153	za1 = ABS(float(jn-jic))/float(nn_baro)
1154	IF (za1 < 0.5_wp) THEN
1155	zwgt1(jn) = 1._wp
1156	jpit = jn
1157	ENDIF
1158	ENDDO
1159
1160	CASE( 2 ) ! Boxcar, width = 2 * nn_baro
1161	DO jn = 1, 3*nn_baro
1162	za1 = ABS(float(jn-jic))/float(nn_baro)
1163	IF (za1 < 1._wp) THEN
1164	zwgt1(jn) = 1._wp
1165	jpit = jn
1166	ENDIF
1167	ENDDO
1168	CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_bt_flt' )
1169	END SELECT
1170
1171	ELSE ! No time averaging
1172	zwgt1(jic) = 1._wp
1173	jpit = jic
1174	ENDIF
1175
1176	! Set secondary weights
1177	DO jn = 1, jpit
1178	DO ji = jn, jpit
1179	zwgt2(jn) = zwgt2(jn) + zwgt1(ji)
1180	END DO
1181	END DO
1182
1183	! Normalize weigths:
1184	za1 = 1._wp / SUM(zwgt1(1:jpit))
1185	za2 = 1._wp / SUM(zwgt2(1:jpit))
1186	DO jn = 1, jpit
1187	zwgt1(jn) = zwgt1(jn) * za1
1188	zwgt2(jn) = zwgt2(jn) * za2
1189	END DO
1190	!
1191	END SUBROUTINE ts_wgt
1192
1193
1194	SUBROUTINE ts_rst( kt, cdrw )
1195	!!---------------------------------------------------------------------
1196	!! * ROUTINE ts_rst *
1197	!!
1198	!! ** Purpose : Read or write time-splitting arrays in restart file
1199	!!----------------------------------------------------------------------
1200	INTEGER , INTENT(in) :: kt ! ocean time-step
1201	CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag
1202	!
1203	!!----------------------------------------------------------------------
1204	!
1205	IF( TRIM(cdrw) == 'READ' ) THEN
1206	CALL iom_get( numror, jpdom_autoglo, 'ub2_b' , ub2_b (:,:) )
1207	CALL iom_get( numror, jpdom_autoglo, 'vb2_b' , vb2_b (:,:) )
1208	CALL iom_get( numror, jpdom_autoglo, 'un_bf' , un_bf (:,:) )
1209	CALL iom_get( numror, jpdom_autoglo, 'vn_bf' , vn_bf (:,:) )
1210	IF( .NOT.ln_bt_av ) THEN
1211	CALL iom_get( numror, jpdom_autoglo, 'sshbb_e' , sshbb_e(:,:) )
1212	CALL iom_get( numror, jpdom_autoglo, 'ubb_e' , ubb_e(:,:) )
1213	CALL iom_get( numror, jpdom_autoglo, 'vbb_e' , vbb_e(:,:) )
1214	CALL iom_get( numror, jpdom_autoglo, 'sshb_e' , sshb_e(:,:) )
1215	CALL iom_get( numror, jpdom_autoglo, 'ub_e' , ub_e(:,:) )
1216	CALL iom_get( numror, jpdom_autoglo, 'vb_e' , vb_e(:,:) )
1217	ENDIF
1218	#if defined key_agrif
1219	! Read time integrated fluxes
1220	IF ( .NOT.Agrif_Root() ) THEN
1221	CALL iom_get( numror, jpdom_autoglo, 'ub2_i_b' , ub2_i_b(:,:) )
1222	CALL iom_get( numror, jpdom_autoglo, 'vb2_i_b' , vb2_i_b(:,:) )
1223	ENDIF
1224	#endif
1225	!
1226	ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN
1227	CALL iom_rstput( kt, nitrst, numrow, 'ub2_b' , ub2_b (:,:) )
1228	CALL iom_rstput( kt, nitrst, numrow, 'vb2_b' , vb2_b (:,:) )
1229	CALL iom_rstput( kt, nitrst, numrow, 'un_bf' , un_bf (:,:) )
1230	CALL iom_rstput( kt, nitrst, numrow, 'vn_bf' , vn_bf (:,:) )
1231	!
1232	IF (.NOT.ln_bt_av) THEN
1233	CALL iom_rstput( kt, nitrst, numrow, 'sshbb_e' , sshbb_e(:,:) )
1234	CALL iom_rstput( kt, nitrst, numrow, 'ubb_e' , ubb_e(:,:) )
1235	CALL iom_rstput( kt, nitrst, numrow, 'vbb_e' , vbb_e(:,:) )
1236	CALL iom_rstput( kt, nitrst, numrow, 'sshb_e' , sshb_e(:,:) )
1237	CALL iom_rstput( kt, nitrst, numrow, 'ub_e' , ub_e(:,:) )
1238	CALL iom_rstput( kt, nitrst, numrow, 'vb_e' , vb_e(:,:) )
1239	ENDIF
1240	#if defined key_agrif
1241	! Save time integrated fluxes
1242	IF ( .NOT.Agrif_Root() ) THEN
1243	CALL iom_rstput( kt, nitrst, numrow, 'ub2_i_b' , ub2_i_b(:,:) )
1244	CALL iom_rstput( kt, nitrst, numrow, 'vb2_i_b' , vb2_i_b(:,:) )
1245	ENDIF
1246	#endif
1247	ENDIF
1248	!
1249	END SUBROUTINE ts_rst
1250
1251
1252	SUBROUTINE dyn_spg_ts_init
1253	!!---------------------------------------------------------------------
1254	!! * ROUTINE dyn_spg_ts_init *
1255	!!
1256	!! ** Purpose : Set time splitting options
1257	!!----------------------------------------------------------------------
1258	INTEGER :: ji ,jj ! dummy loop indices
1259	REAL(wp) :: zxr2, zyr2, zcmax ! local scalar
1260	REAL(wp), POINTER, DIMENSION(:,:) :: zcu
1261	!!----------------------------------------------------------------------
1262	!
1263	! Max courant number for ext. grav. waves
1264	!
1265	CALL wrk_alloc( jpi,jpj, zcu )
1266	!
1267	DO jj = 1, jpj
1268	DO ji =1, jpi
1269	zxr2 = r1_e1t(ji,jj) * r1_e1t(ji,jj)
1270	zyr2 = r1_e2t(ji,jj) * r1_e2t(ji,jj)
1271	zcu(ji,jj) = SQRT( grav * MAX(ht_0(ji,jj),0._wp) * (zxr2 + zyr2) )
1272	END DO
1273	END DO
1274	!
1275	zcmax = MAXVAL( zcu(:,:) )
1276	IF( lk_mpp ) CALL mpp_max( zcmax )
1277
1278	! Estimate number of iterations to satisfy a max courant number= rn_bt_cmax
1279	IF( ln_bt_auto ) nn_baro = CEILING( rdt / rn_bt_cmax * zcmax)
1280
1281	rdtbt = rdt / REAL( nn_baro , wp )
1282	zcmax = zcmax * rdtbt
1283	! Print results
1284	IF(lwp) WRITE(numout,*)
1285	IF(lwp) WRITE(numout,*) 'dyn_spg_ts : split-explicit free surface'
1286	IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
1287	IF( ln_bt_auto ) THEN
1288	IF(lwp) WRITE(numout,*) ' ln_ts_auto=.true. Automatically set nn_baro '
1289	IF(lwp) WRITE(numout,*) ' Max. courant number allowed: ', rn_bt_cmax
1290	ELSE
1291	IF(lwp) WRITE(numout,*) ' ln_ts_auto=.false.: Use nn_baro in namelist '
1292	ENDIF
1293
1294	IF(ln_bt_av) THEN
1295	IF(lwp) WRITE(numout,*) ' ln_bt_av=.true. => Time averaging over nn_baro time steps is on '
1296	ELSE
1297	IF(lwp) WRITE(numout,*) ' ln_bt_av=.false. => No time averaging of barotropic variables '
1298	ENDIF
1299	!
1300	!
1301	IF(ln_bt_fw) THEN
1302	IF(lwp) WRITE(numout,*) ' ln_bt_fw=.true. => Forward integration of barotropic variables '
1303	ELSE
1304	IF(lwp) WRITE(numout,*) ' ln_bt_fw =.false.=> Centred integration of barotropic variables '
1305	ENDIF
1306	!
1307	#if defined key_agrif
1308	! Restrict the use of Agrif to the forward case only
1309	IF( .NOT.ln_bt_fw .AND. .NOT.Agrif_Root() ) CALL ctl_stop( 'AGRIF not implemented if ln_bt_fw=.FALSE.' )
1310	#endif
1311	!
1312	IF(lwp) WRITE(numout,*) ' Time filter choice, nn_bt_flt: ', nn_bt_flt
1313	SELECT CASE ( nn_bt_flt )
1314	CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' Dirac'
1315	CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' Boxcar: width = nn_baro'
1316	CASE( 2 ) ; IF(lwp) WRITE(numout,) ' Boxcar: width = 2nn_baro'
1317	CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_bt_flt: should 0,1,2' )
1318	END SELECT
1319	!
1320	IF(lwp) WRITE(numout,*) ' '
1321	IF(lwp) WRITE(numout,*) ' nn_baro = ', nn_baro
1322	IF(lwp) WRITE(numout,*) ' Barotropic time step [s] is :', rdtbt
1323	IF(lwp) WRITE(numout,*) ' Maximum Courant number is :', zcmax
1324	!
1325	IF( .NOT.ln_bt_av .AND. .NOT.ln_bt_fw ) THEN
1326	CALL ctl_stop( 'dynspg_ts ERROR: No time averaging => only forward integration is possible' )
1327	ENDIF
1328	IF( zcmax>0.9_wp ) THEN
1329	CALL ctl_stop( 'dynspg_ts ERROR: Maximum Courant number is greater than 0.9: Inc. nn_baro !' )
1330	ENDIF
1331	!
1332	CALL wrk_dealloc( jpi,jpj, zcu )
1333	!
1334	END SUBROUTINE dyn_spg_ts_init
1335
1336	!!======================================================================
1337	END MODULE dynspg_ts

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

Download in other formats: