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 NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/OCE/DYN – NEMO

Context Navigation

source: NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/OCE/DYN/dynspg_ts.F90 @ 10420

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

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

Download in other formats: