Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

dynspg_ts.F90 @ 13325

Last change on this file since 13325 was 13325, checked in by andmirek, 4 years ago
Ticket ##2482: Use explicit index in dynspg_ts.F90
File size: 77.3 KB

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

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

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

Context Navigation

source: NEMO/branches/UKMO/NEMO_4.0_mirror_SI3_GPU/src/OCE/DYN/dynspg_ts.F90 @ 13325

Download in other formats: