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

Context Navigation

source: branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 8093

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

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

Download in other formats: