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/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

Context Navigation

source: branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 7432

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

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

Download in other formats: