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

Context Navigation

source: branches/2016/dev_r6393_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 7157

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

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

Download in other formats: