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/2019/ENHANCE-02_ISF_nemo/src/OCE/DYN – NEMO

Context Navigation

source: NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/DYN/dynspg_ts.F90 @ 11521

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

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

Download in other formats: