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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

dynspg_ts.F90 in branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

Context Navigation

source: branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 9019

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

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

Download in other formats: