Context Navigation

dynspg_ts.F90 @ 8831

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

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

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

Context Navigation

source: branches/2017/dev_r8600_xios_read_write_v2/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 8831

Download in other formats: