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

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

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

Context Navigation

source: branches/2017/dev_r8624_AGRIF3_VVL/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 8898

Last change on this file since 8898 was 8898, checked in by jchanut, 6 years ago
AGRIF+VVL: Merge with free surface changes - #1965
Property svn:keywords set to `Id`
File size: 64.8 KB

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

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

Download in other formats: