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/dev_r10742_ENHANCE-12_SimonM-Tides/src/OCE/DYN – NEMO

Context Navigation

source: NEMO/branches/2019/dev_r10742_ENHANCE-12_SimonM-Tides/src/OCE/DYN/dynspg_ts.F90 @ 10773

Last change on this file since 10773 was 10773, checked in by smueller, 5 years ago
Transfer of five public variables, their allocation, and two subroutines from module sbctide to module tide_mod (ticket #2194)
Property svn:keywords set to `Id`
File size: 77.2 KB

Line
1	MODULE dynspg_ts
2
3	!! Includes ROMS wd scheme with diagnostic outputs ; un and ua updates are commented out !
4
5	!!======================================================================
6	!! * MODULE dynspg_ts *
7	!! Ocean dynamics: surface pressure gradient trend, split-explicit scheme
8	!!======================================================================
9	!! History : 1.0 ! 2004-12 (L. Bessieres, G. Madec) Original code
10	!! - ! 2005-11 (V. Garnier, G. Madec) optimization
11	!! - ! 2006-08 (S. Masson) distributed restart using iom
12	!! 2.0 ! 2007-07 (D. Storkey) calls to BDY routines
13	!! - ! 2008-01 (R. Benshila) change averaging method
14	!! 3.2 ! 2009-07 (R. Benshila, G. Madec) Complete revisit associated to vvl reactivation
15	!! 3.3 ! 2010-09 (D. Storkey, E. O'Dea) update for BDY for Shelf configurations
16	!! 3.3 ! 2011-03 (R. Benshila, R. Hordoir, P. Oddo) update calculation of ub_b
17	!! 3.5 ! 2013-07 (J. Chanut) Switch to Forward-backward time stepping
18	!! 3.6 ! 2013-11 (A. Coward) Update for z-tilde compatibility
19	!! 3.7 ! 2015-11 (J. Chanut) free surface simplification
20	!! - ! 2016-12 (G. Madec, E. Clementi) update for Stoke-Drift divergence
21	!! 4.0 ! 2017-05 (G. Madec) drag coef. defined at t-point (zdfdrg.F90)
22	!!---------------------------------------------------------------------
23
24	!!----------------------------------------------------------------------
25	!! dyn_spg_ts : compute surface pressure gradient trend using a time-splitting scheme
26	!! dyn_spg_ts_init: initialisation of the time-splitting scheme
27	!! ts_wgt : set time-splitting weights for temporal averaging (or not)
28	!! ts_rst : read/write time-splitting fields in restart file
29	!!----------------------------------------------------------------------
30	USE oce ! ocean dynamics and tracers
31	USE dom_oce ! ocean space and time domain
32	USE sbc_oce ! surface boundary condition: ocean
33	USE zdf_oce ! vertical physics: variables
34	USE zdfdrg ! vertical physics: top/bottom drag coef.
35	USE sbcisf ! ice shelf variable (fwfisf)
36	USE sbcapr ! surface boundary condition: atmospheric pressure
37	USE dynadv , ONLY: ln_dynadv_vec
38	USE dynvor ! vortivity scheme indicators
39	USE phycst ! physical constants
40	USE dynvor ! vorticity term
41	USE wet_dry ! wetting/drying flux limter
42	USE bdy_oce ! open boundary
43	USE bdyvol ! open boundary volume conservation
44	USE bdytides ! open boundary condition data
45	USE bdydyn2d ! open boundary conditions on barotropic variables
46	USE sbctide ! tides
47	USE updtide ! tide potential
48	USE tide_mod !
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(:,:) )
636	ELSE
637	zztmp = r1_rau0 * r1_2
638	zssh_frc(:,:) = zztmp * ( emp(:,:) + emp_b(:,:) - rnf(:,:) - rnf_b(:,:) &
639	& + fwfisf(:,:) + fwfisf_b(:,:) )
640	ENDIF
641	!
642	IF( ln_sdw ) THEN ! Stokes drift divergence added if necessary
643	zssh_frc(:,:) = zssh_frc(:,:) + div_sd(:,:)
644	ENDIF
645	!
646	#if defined key_asminc
647	! ! Include the IAU weighted SSH increment
648	IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN
649	zssh_frc(:,:) = zssh_frc(:,:) - ssh_iau(:,:)
650	ENDIF
651	#endif
652	! !* Fill boundary data arrays for AGRIF
653	! ! ------------------------------------
654	#if defined key_agrif
655	IF( .NOT.Agrif_Root() ) CALL agrif_dta_ts( kt )
656	#endif
657	!
658	! -----------------------------------------------------------------------
659	! Phase 2 : Integration of the barotropic equations
660	! -----------------------------------------------------------------------
661	!
662	! ! ==================== !
663	! ! Initialisations !
664	! ! ==================== !
665	! Initialize barotropic variables:
666	IF( ll_init )THEN
667	sshbb_e(:,:) = 0._wp
668	ubb_e (:,:) = 0._wp
669	vbb_e (:,:) = 0._wp
670	sshb_e (:,:) = 0._wp
671	ub_e (:,:) = 0._wp
672	vb_e (:,:) = 0._wp
673	ENDIF
674
675	!
676	IF (ln_bt_fw) THEN ! FORWARD integration: start from NOW fields
677	sshn_e(:,:) = sshn(:,:)
678	un_e (:,:) = un_b(:,:)
679	vn_e (:,:) = vn_b(:,:)
680	!
681	hu_e (:,:) = hu_n(:,:)
682	hv_e (:,:) = hv_n(:,:)
683	hur_e (:,:) = r1_hu_n(:,:)
684	hvr_e (:,:) = r1_hv_n(:,:)
685	ELSE ! CENTRED integration: start from BEFORE fields
686	sshn_e(:,:) = sshb(:,:)
687	un_e (:,:) = ub_b(:,:)
688	vn_e (:,:) = vb_b(:,:)
689	!
690	hu_e (:,:) = hu_b(:,:)
691	hv_e (:,:) = hv_b(:,:)
692	hur_e (:,:) = r1_hu_b(:,:)
693	hvr_e (:,:) = r1_hv_b(:,:)
694	ENDIF
695	!
696	!
697	!
698	! Initialize sums:
699	ua_b (:,:) = 0._wp ! After barotropic velocities (or transport if flux form)
700	va_b (:,:) = 0._wp
701	ssha (:,:) = 0._wp ! Sum for after averaged sea level
702	un_adv(:,:) = 0._wp ! Sum for now transport issued from ts loop
703	vn_adv(:,:) = 0._wp
704	!
705	IF( ln_wd_dl ) THEN
706	zuwdmask(:,:) = 0._wp ! set to zero for definiteness (not sure this is necessary)
707	zvwdmask(:,:) = 0._wp !
708	zuwdav2 (:,:) = 0._wp
709	zvwdav2 (:,:) = 0._wp
710	END IF
711
712	! ! ==================== !
713	DO jn = 1, icycle ! sub-time-step loop !
714	! ! ==================== !
715	!
716	l_full_nf_update = jn == icycle ! false: disable full North fold update (performances) for jn = 1 to icycle-1
717	! ! ------------------
718	! !* Update the forcing (BDY and tides)
719	! ! ------------------
720	! Update only tidal forcing at open boundaries
721	IF( ln_bdy .AND. ln_tide ) CALL bdy_dta_tides( kt, kit=jn, time_offset= noffset+1 )
722	IF( ln_tide_pot .AND. ln_tide ) CALL upd_tide ( kt, kit=jn, time_offset= noffset )
723	!
724	! Set extrapolation coefficients for predictor step:
725	IF ((jn<3).AND.ll_init) THEN ! Forward
726	za1 = 1._wp
727	za2 = 0._wp
728	za3 = 0._wp
729	ELSE ! AB3-AM4 Coefficients: bet=0.281105
730	za1 = 1.781105_wp ! za1 = 3/2 + bet
731	za2 = -1.06221_wp ! za2 = -(1/2 + 2*bet)
732	za3 = 0.281105_wp ! za3 = bet
733	ENDIF
734
735	! Extrapolate barotropic velocities at step jit+0.5:
736	ua_e(:,:) = za1 * un_e(:,:) + za2 * ub_e(:,:) + za3 * ubb_e(:,:)
737	va_e(:,:) = za1 * vn_e(:,:) + za2 * vb_e(:,:) + za3 * vbb_e(:,:)
738
739	IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only)
740	! ! ------------------
741	! Extrapolate Sea Level at step jit+0.5:
742	zsshp2_e(:,:) = za1 * sshn_e(:,:) + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:)
743
744	! set wetting & drying mask at tracer points for this barotropic sub-step
745	IF ( ln_wd_dl ) THEN
746	!
747	IF ( ln_wd_dl_rmp ) THEN
748	DO jj = 1, jpj
749	DO ji = 1, jpi ! vector opt.
750	IF ( zsshp2_e(ji,jj) + ht_0(ji,jj) > 2._wp * rn_wdmin1 ) THEN
751	! IF ( zsshp2_e(ji,jj) + ht_0(ji,jj) > rn_wdmin2 ) THEN
752	ztwdmask(ji,jj) = 1._wp
753	ELSE IF ( zsshp2_e(ji,jj) + ht_0(ji,jj) > rn_wdmin1 ) THEN
754	ztwdmask(ji,jj) = (tanh(50._wp( ( zsshp2_e(ji,jj) + ht_0(ji,jj) - rn_wdmin1 )r_rn_wdmin1)) )
755	ELSE
756	ztwdmask(ji,jj) = 0._wp
757	END IF
758	END DO
759	END DO
760	ELSE
761	DO jj = 1, jpj
762	DO ji = 1, jpi ! vector opt.
763	IF ( zsshp2_e(ji,jj) + ht_0(ji,jj) > rn_wdmin1 ) THEN
764	ztwdmask(ji,jj) = 1._wp
765	ELSE
766	ztwdmask(ji,jj) = 0._wp
767	ENDIF
768	END DO
769	END DO
770	ENDIF
771	!
772	ENDIF
773	!
774	DO jj = 2, jpjm1 ! Sea Surface Height at u- & v-points
775	DO ji = 2, fs_jpim1 ! Vector opt.
776	zwx(ji,jj) = r1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj) &
777	& * ( e1e2t(ji ,jj) * zsshp2_e(ji ,jj) &
778	& + e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj) )
779	zwy(ji,jj) = r1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj) &
780	& * ( e1e2t(ji,jj ) * zsshp2_e(ji,jj ) &
781	& + e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1) )
782	END DO
783	END DO
784	CALL lbc_lnk_multi( 'dynspg_ts', zwx, 'U', 1._wp, zwy, 'V', 1._wp )
785	!
786	zhup2_e(:,:) = hu_0(:,:) + zwx(:,:) ! Ocean depth at U- and V-points
787	zhvp2_e(:,:) = hv_0(:,:) + zwy(:,:)
788	zhtp2_e(:,:) = ht_0(:,:) + zsshp2_e(:,:)
789	ELSE
790	zhup2_e(:,:) = hu_n(:,:)
791	zhvp2_e(:,:) = hv_n(:,:)
792	zhtp2_e(:,:) = ht_n(:,:)
793	ENDIF
794	! !* after ssh
795	! ! -----------
796	!
797	! Enforce volume conservation at open boundaries:
798	IF( ln_bdy .AND. ln_vol ) CALL bdy_vol2d( kt, jn, ua_e, va_e, zhup2_e, zhvp2_e )
799	!
800	zwx(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:) ! fluxes at jn+0.5
801	zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:)
802	!
803	#if defined key_agrif
804	! Set fluxes during predictor step to ensure volume conservation
805	IF( .NOT.Agrif_Root() .AND. ln_bt_fw ) THEN
806	IF((nbondi == -1).OR.(nbondi == 2)) THEN
807	DO jj = 1, jpj
808	zwx(2:nbghostcells+1,jj) = ubdy_w(1:nbghostcells,jj) * e2u(2:nbghostcells+1,jj)
809	zwy(2:nbghostcells+1,jj) = vbdy_w(1:nbghostcells,jj) * e1v(2:nbghostcells+1,jj)
810	END DO
811	ENDIF
812	IF((nbondi == 1).OR.(nbondi == 2)) THEN
813	DO jj=1,jpj
814	zwx(nlci-nbghostcells-1:nlci-2,jj) = ubdy_e(1:nbghostcells,jj) * e2u(nlci-nbghostcells-1:nlci-2,jj)
815	zwy(nlci-nbghostcells :nlci-1,jj) = vbdy_e(1:nbghostcells,jj) * e1v(nlci-nbghostcells :nlci-1,jj)
816	END DO
817	ENDIF
818	IF((nbondj == -1).OR.(nbondj == 2)) THEN
819	DO ji=1,jpi
820	zwy(ji,2:nbghostcells+1) = vbdy_s(ji,1:nbghostcells) * e1v(ji,2:nbghostcells+1)
821	zwx(ji,2:nbghostcells+1) = ubdy_s(ji,1:nbghostcells) * e2u(ji,2:nbghostcells+1)
822	END DO
823	ENDIF
824	IF((nbondj == 1).OR.(nbondj == 2)) THEN
825	DO ji=1,jpi
826	zwy(ji,nlcj-nbghostcells-1:nlcj-2) = vbdy_n(ji,1:nbghostcells) * e1v(ji,nlcj-nbghostcells-1:nlcj-2)
827	zwx(ji,nlcj-nbghostcells :nlcj-1) = ubdy_n(ji,1:nbghostcells) * e2u(ji,nlcj-nbghostcells :nlcj-1)
828	END DO
829	ENDIF
830	ENDIF
831	#endif
832	IF( ln_wd_il ) CALL wad_lmt_bt(zwx, zwy, sshn_e, zssh_frc, rdtbt)
833
834	IF ( ln_wd_dl ) THEN
835	!
836	! un_e and vn_e are set to zero at faces where the direction of the flow is from dry cells
837	!
838	DO jj = 1, jpjm1
839	DO ji = 1, jpim1
840	IF ( zwx(ji,jj) > 0.0 ) THEN
841	zuwdmask(ji, jj) = ztwdmask(ji ,jj)
842	ELSE
843	zuwdmask(ji, jj) = ztwdmask(ji+1,jj)
844	END IF
845	zwx(ji, jj) = zuwdmask(ji,jj)*zwx(ji, jj)
846	un_e(ji,jj) = zuwdmask(ji,jj)*un_e(ji,jj)
847
848	IF ( zwy(ji,jj) > 0.0 ) THEN
849	zvwdmask(ji, jj) = ztwdmask(ji, jj )
850	ELSE
851	zvwdmask(ji, jj) = ztwdmask(ji, jj+1)
852	END IF
853	zwy(ji, jj) = zvwdmask(ji,jj)*zwy(ji,jj)
854	vn_e(ji,jj) = zvwdmask(ji,jj)*vn_e(ji,jj)
855	END DO
856	END DO
857	!
858	ENDIF
859
860	! Sum over sub-time-steps to compute advective velocities
861	za2 = wgtbtp2(jn)
862	un_adv(:,:) = un_adv(:,:) + za2 * zwx(:,:) * r1_e2u(:,:)
863	vn_adv(:,:) = vn_adv(:,:) + za2 * zwy(:,:) * r1_e1v(:,:)
864
865	! sum over sub-time-steps to decide which baroclinic velocities to set to zero (zuwdav2 is only used when ln_wd_dl_bc = True)
866	IF ( ln_wd_dl_bc ) THEN
867	zuwdav2(:,:) = zuwdav2(:,:) + za2 * zuwdmask(:,:)
868	zvwdav2(:,:) = zvwdav2(:,:) + za2 * zvwdmask(:,:)
869	END IF
870
871	! Set next sea level:
872	DO jj = 2, jpjm1
873	DO ji = fs_2, fs_jpim1 ! vector opt.
874	zhdiv(ji,jj) = ( zwx(ji,jj) - zwx(ji-1,jj) &
875	& + zwy(ji,jj) - zwy(ji,jj-1) ) * r1_e1e2t(ji,jj)
876	END DO
877	END DO
878	ssha_e(:,:) = ( sshn_e(:,:) - rdtbt * ( zssh_frc(:,:) + zhdiv(:,:) ) ) * ssmask(:,:)
879
880	CALL lbc_lnk( 'dynspg_ts', ssha_e, 'T', 1._wp )
881
882	! Duplicate sea level across open boundaries (this is only cosmetic if linssh=T)
883	IF( ln_bdy ) CALL bdy_ssh( ssha_e )
884	#if defined key_agrif
885	IF( .NOT.Agrif_Root() ) CALL agrif_ssh_ts( jn )
886	#endif
887	!
888	! Sea Surface Height at u-,v-points (vvl case only)
889	IF( .NOT.ln_linssh ) THEN
890	DO jj = 2, jpjm1
891	DO ji = 2, jpim1 ! NO Vector Opt.
892	zsshu_a(ji,jj) = r1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj) &
893	& * ( e1e2t(ji ,jj ) * ssha_e(ji ,jj ) &
894	& + e1e2t(ji+1,jj ) * ssha_e(ji+1,jj ) )
895	zsshv_a(ji,jj) = r1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj) &
896	& * ( e1e2t(ji ,jj ) * ssha_e(ji ,jj ) &
897	& + e1e2t(ji ,jj+1) * ssha_e(ji ,jj+1) )
898	END DO
899	END DO
900	CALL lbc_lnk_multi( 'dynspg_ts', zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp )
901	ENDIF
902	!
903	! Half-step back interpolation of SSH for surface pressure computation:
904	!----------------------------------------------------------------------
905	IF ((jn==1).AND.ll_init) THEN
906	za0=1._wp ! Forward-backward
907	za1=0._wp
908	za2=0._wp
909	za3=0._wp
910	ELSEIF ((jn==2).AND.ll_init) THEN ! AB2-AM3 Coefficients; bet=0 ; gam=-1/6 ; eps=1/12
911	za0= 1.0833333333333_wp ! za0 = 1-gam-eps
912	za1=-0.1666666666666_wp ! za1 = gam
913	za2= 0.0833333333333_wp ! za2 = eps
914	za3= 0._wp
915	ELSE ! AB3-AM4 Coefficients; bet=0.281105 ; eps=0.013 ; gam=0.0880
916	IF (rn_bt_alpha==0._wp) THEN
917	za0=0.614_wp ! za0 = 1/2 + gam + 2*eps
918	za1=0.285_wp ! za1 = 1/2 - 2gam - 3eps
919	za2=0.088_wp ! za2 = gam
920	za3=0.013_wp ! za3 = eps
921	ELSE
922	zepsilon = 0.00976186_wp - 0.13451357_wp * rn_bt_alpha
923	zgamma = 0.08344500_wp - 0.51358400_wp * rn_bt_alpha
924	za0 = 0.5_wp + zgamma + 2._wp * rn_bt_alpha + 2._wp * zepsilon
925	za1 = 1._wp - za0 - zgamma - zepsilon
926	za2 = zgamma
927	za3 = zepsilon
928	ENDIF
929	ENDIF
930	!
931	zsshp2_e(:,:) = za0 * ssha_e(:,:) + za1 * sshn_e (:,:) &
932	& + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:)
933
934	IF( ln_wd_il ) THEN ! Calculating and applying W/D gravity filters
935	DO jj = 2, jpjm1
936	DO ji = 2, jpim1
937	ll_tmp1 = MIN( zsshp2_e(ji,jj) , zsshp2_e(ji+1,jj) ) > &
938	& MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) .AND. &
939	& MAX( zsshp2_e(ji,jj) + ht_0(ji,jj) , zsshp2_e(ji+1,jj) + ht_0(ji+1,jj) ) &
940	& > rn_wdmin1 + rn_wdmin2
941	ll_tmp2 = (ABS(zsshp2_e(ji,jj) - zsshp2_e(ji+1,jj)) > 1.E-12 ).AND.( &
942	& MAX( zsshp2_e(ji,jj) , zsshp2_e(ji+1,jj) ) > &
943	& MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 )
944
945	IF(ll_tmp1) THEN
946	zcpx(ji,jj) = 1.0_wp
947	ELSE IF(ll_tmp2) THEN
948	! no worries about zsshp2_e(ji+1,jj) - zsshp2_e(ji ,jj) = 0, it won't happen ! here
949	zcpx(ji,jj) = ABS( (zsshp2_e(ji+1,jj) + ht_0(ji+1,jj) - zsshp2_e(ji,jj) - ht_0(ji,jj)) &
950	& / (zsshp2_e(ji+1,jj) - zsshp2_e(ji ,jj)) )
951	ELSE
952	zcpx(ji,jj) = 0._wp
953	ENDIF
954	!
955	ll_tmp1 = MIN( zsshp2_e(ji,jj) , zsshp2_e(ji,jj+1) ) > &
956	& MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) .AND. &
957	& MAX( zsshp2_e(ji,jj) + ht_0(ji,jj) , zsshp2_e(ji,jj+1) + ht_0(ji,jj+1) ) &
958	& > rn_wdmin1 + rn_wdmin2
959	ll_tmp2 = (ABS(zsshp2_e(ji,jj) - zsshp2_e(ji,jj+1)) > 1.E-12 ).AND.( &
960	& MAX( zsshp2_e(ji,jj) , zsshp2_e(ji,jj+1) ) > &
961	& MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 )
962
963	IF(ll_tmp1) THEN
964	zcpy(ji,jj) = 1.0_wp
965	ELSEIF(ll_tmp2) THEN
966	! no worries about zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj ) = 0, it won't happen ! here
967	zcpy(ji,jj) = ABS( (zsshp2_e(ji,jj+1) + ht_0(ji,jj+1) - zsshp2_e(ji,jj) - ht_0(ji,jj)) &
968	& / (zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj )) )
969	ELSE
970	zcpy(ji,jj) = 0._wp
971	ENDIF
972	END DO
973	END DO
974	ENDIF
975	!
976	! Compute associated depths at U and V points:
977	IF( .NOT.ln_linssh .AND. .NOT.ln_dynadv_vec ) THEN !* Vector form
978	!
979	DO jj = 2, jpjm1
980	DO ji = 2, jpim1
981	zx1 = r1_2 * ssumask(ji ,jj) * r1_e1e2u(ji ,jj) &
982	& * ( e1e2t(ji ,jj ) * zsshp2_e(ji ,jj) &
983	& + e1e2t(ji+1,jj ) * zsshp2_e(ji+1,jj ) )
984	zy1 = r1_2 * ssvmask(ji ,jj) * r1_e1e2v(ji ,jj ) &
985	& * ( e1e2t(ji ,jj ) * zsshp2_e(ji ,jj ) &
986	& + e1e2t(ji ,jj+1) * zsshp2_e(ji ,jj+1) )
987	zhust_e(ji,jj) = hu_0(ji,jj) + zx1
988	zhvst_e(ji,jj) = hv_0(ji,jj) + zy1
989	END DO
990	END DO
991	!
992	ENDIF
993	!
994	! Add Coriolis trend:
995	! zwz array below or triads normally depend on sea level with ln_linssh=F and should be updated
996	! at each time step. We however keep them constant here for optimization.
997	! Recall that zwx and zwy arrays hold fluxes at this stage:
998	! zwx(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:) ! fluxes at jn+0.5
999	! zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:)
1000	!
1001	SELECT CASE( nvor_scheme )
1002	CASE( np_ENT ) ! energy conserving scheme (t-point)
1003	DO jj = 2, jpjm1
1004	DO ji = 2, jpim1 ! vector opt.
1005
1006	z1_hu = ssumask(ji,jj) / ( zhup2_e(ji,jj) + 1._wp - ssumask(ji,jj) )
1007	z1_hv = ssvmask(ji,jj) / ( zhvp2_e(ji,jj) + 1._wp - ssvmask(ji,jj) )
1008
1009	zu_trd(ji,jj) = + r1_4 * r1_e1e2u(ji,jj) * z1_hu &
1010	& * ( 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) ) &
1011	& + e1e2t(ji ,jj)zhtp2_e(ji ,jj)ff_t(ji ,jj) * ( va_e(ji ,jj) + va_e(ji ,jj-1) ) )
1012	!
1013	zv_trd(ji,jj) = - r1_4 * r1_e1e2v(ji,jj) * z1_hv &
1014	& * ( 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) ) &
1015	& + e1e2t(ji,jj )zhtp2_e(ji,jj )ff_t(ji,jj ) * ( ua_e(ji,jj ) + ua_e(ji-1,jj ) ) )
1016	END DO
1017	END DO
1018	!
1019	CASE( np_ENE, np_MIX ) ! energy conserving scheme (f-point)
1020	DO jj = 2, jpjm1
1021	DO ji = fs_2, fs_jpim1 ! vector opt.
1022	zy1 = ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) ) * r1_e1u(ji,jj)
1023	zy2 = ( zwy(ji ,jj ) + zwy(ji+1,jj ) ) * r1_e1u(ji,jj)
1024	zx1 = ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) ) * r1_e2v(ji,jj)
1025	zx2 = ( zwx(ji ,jj ) + zwx(ji ,jj+1) ) * r1_e2v(ji,jj)
1026	zu_trd(ji,jj) = r1_4 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
1027	zv_trd(ji,jj) =-r1_4 * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 )
1028	END DO
1029	END DO
1030	!
1031	CASE( np_ENS ) ! enstrophy conserving scheme (f-point)
1032	DO jj = 2, jpjm1
1033	DO ji = fs_2, fs_jpim1 ! vector opt.
1034	zy1 = r1_8 * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) &
1035	& + zwy(ji ,jj ) + zwy(ji+1,jj ) ) * r1_e1u(ji,jj)
1036	zx1 = - r1_8 * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) &
1037	& + zwx(ji ,jj ) + zwx(ji ,jj+1) ) * r1_e2v(ji,jj)
1038	zu_trd(ji,jj) = zy1 * ( zwz(ji ,jj-1) + zwz(ji,jj) )
1039	zv_trd(ji,jj) = zx1 * ( zwz(ji-1,jj ) + zwz(ji,jj) )
1040	END DO
1041	END DO
1042	!
1043	CASE( np_EET , np_EEN ) ! energy & enstrophy scheme (using e3t or e3f)
1044	DO jj = 2, jpjm1
1045	DO ji = fs_2, fs_jpim1 ! vector opt.
1046	zu_trd(ji,jj) = + r1_12 * r1_e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) &
1047	& + ftnw(ji+1,jj) * zwy(ji+1,jj ) &
1048	& + ftse(ji,jj ) * zwy(ji ,jj-1) &
1049	& + ftsw(ji+1,jj) * zwy(ji+1,jj-1) )
1050	zv_trd(ji,jj) = - r1_12 * r1_e2v(ji,jj) * ( ftsw(ji,jj+1) * zwx(ji-1,jj+1) &
1051	& + ftse(ji,jj+1) * zwx(ji ,jj+1) &
1052	& + ftnw(ji,jj ) * zwx(ji-1,jj ) &
1053	& + ftne(ji,jj ) * zwx(ji ,jj ) )
1054	END DO
1055	END DO
1056	!
1057	END SELECT
1058	!
1059	! Add tidal astronomical forcing if defined
1060	IF ( ln_tide .AND. ln_tide_pot ) THEN
1061	DO jj = 2, jpjm1
1062	DO ji = fs_2, fs_jpim1 ! vector opt.
1063	zu_spg = grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) * r1_e1u(ji,jj)
1064	zv_spg = grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) * r1_e2v(ji,jj)
1065	zu_trd(ji,jj) = zu_trd(ji,jj) + zu_spg
1066	zv_trd(ji,jj) = zv_trd(ji,jj) + zv_spg
1067	END DO
1068	END DO
1069	ENDIF
1070	!
1071	! Add bottom stresses:
1072	!jth do implicitly instead
1073	IF ( .NOT. ll_wd ) THEN ! Revert to explicit for bit comparison tests in non wad runs
1074	DO jj = 2, jpjm1
1075	DO ji = fs_2, fs_jpim1 ! vector opt.
1076	zu_trd(ji,jj) = zu_trd(ji,jj) + zCdU_u(ji,jj) * un_e(ji,jj) * hur_e(ji,jj)
1077	zv_trd(ji,jj) = zv_trd(ji,jj) + zCdU_v(ji,jj) * vn_e(ji,jj) * hvr_e(ji,jj)
1078	END DO
1079	END DO
1080	ENDIF
1081	!
1082	! Surface pressure trend:
1083	IF( ln_wd_il ) THEN
1084	DO jj = 2, jpjm1
1085	DO ji = 2, jpim1
1086	! Add surface pressure gradient
1087	zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj)
1088	zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj)
1089	zwx(ji,jj) = (1._wp - rn_scal_load) * zu_spg * zcpx(ji,jj)
1090	zwy(ji,jj) = (1._wp - rn_scal_load) * zv_spg * zcpy(ji,jj)
1091	END DO
1092	END DO
1093	ELSE
1094	DO jj = 2, jpjm1
1095	DO ji = fs_2, fs_jpim1 ! vector opt.
1096	! Add surface pressure gradient
1097	zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj)
1098	zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj)
1099	zwx(ji,jj) = (1._wp - rn_scal_load) * zu_spg
1100	zwy(ji,jj) = (1._wp - rn_scal_load) * zv_spg
1101	END DO
1102	END DO
1103	END IF
1104
1105	!
1106	! Set next velocities:
1107	IF( ln_dynadv_vec .OR. ln_linssh ) THEN !* Vector form
1108	DO jj = 2, jpjm1
1109	DO ji = fs_2, fs_jpim1 ! vector opt.
1110	ua_e(ji,jj) = ( un_e(ji,jj) &
1111	& + rdtbt * ( zwx(ji,jj) &
1112	& + zu_trd(ji,jj) &
1113	& + zu_frc(ji,jj) ) &
1114	& ) * ssumask(ji,jj)
1115
1116	va_e(ji,jj) = ( vn_e(ji,jj) &
1117	& + rdtbt * ( zwy(ji,jj) &
1118	& + zv_trd(ji,jj) &
1119	& + zv_frc(ji,jj) ) &
1120	& ) * ssvmask(ji,jj)
1121
1122	END DO
1123	END DO
1124	!
1125	ELSE !* Flux form
1126	DO jj = 2, jpjm1
1127	DO ji = fs_2, fs_jpim1 ! vector opt.
1128
1129	zhura = hu_0(ji,jj) + zsshu_a(ji,jj)
1130	zhvra = hv_0(ji,jj) + zsshv_a(ji,jj)
1131
1132	zhura = ssumask(ji,jj)/(zhura + 1._wp - ssumask(ji,jj))
1133	zhvra = ssvmask(ji,jj)/(zhvra + 1._wp - ssvmask(ji,jj))
1134
1135	ua_e(ji,jj) = ( hu_e(ji,jj) * un_e(ji,jj) &
1136	& + rdtbt * ( zhust_e(ji,jj) * zwx(ji,jj) &
1137	& + zhup2_e(ji,jj) * zu_trd(ji,jj) &
1138	& + hu_n(ji,jj) * zu_frc(ji,jj) ) &
1139	& ) * zhura
1140
1141	va_e(ji,jj) = ( hv_e(ji,jj) * vn_e(ji,jj) &
1142	& + rdtbt * ( zhvst_e(ji,jj) * zwy(ji,jj) &
1143	& + zhvp2_e(ji,jj) * zv_trd(ji,jj) &
1144	& + hv_n(ji,jj) * zv_frc(ji,jj) ) &
1145	& ) * zhvra
1146	END DO
1147	END DO
1148	ENDIF
1149	!jth implicit bottom friction:
1150	IF ( ll_wd ) THEN ! revert to explicit for bit comparison tests in non wad runs
1151	DO jj = 2, jpjm1
1152	DO ji = fs_2, fs_jpim1 ! vector opt.
1153	ua_e(ji,jj) = ua_e(ji,jj) /(1.0 - rdtbt * zCdU_u(ji,jj) * hur_e(ji,jj))
1154	va_e(ji,jj) = va_e(ji,jj) /(1.0 - rdtbt * zCdU_v(ji,jj) * hvr_e(ji,jj))
1155	END DO
1156	END DO
1157	ENDIF
1158
1159
1160	IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only)
1161	hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:)
1162	hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:)
1163	hur_e(:,:) = ssumask(:,:) / ( hu_e(:,:) + 1._wp - ssumask(:,:) )
1164	hvr_e(:,:) = ssvmask(:,:) / ( hv_e(:,:) + 1._wp - ssvmask(:,:) )
1165	!
1166	ENDIF
1167	! !* domain lateral boundary
1168	CALL lbc_lnk_multi( 'dynspg_ts', ua_e, 'U', -1._wp, va_e , 'V', -1._wp )
1169	!
1170	! ! open boundaries
1171	IF( ln_bdy ) CALL bdy_dyn2d( jn, ua_e, va_e, un_e, vn_e, hur_e, hvr_e, ssha_e )
1172	#if defined key_agrif
1173	IF( .NOT.Agrif_Root() ) CALL agrif_dyn_ts( jn ) ! Agrif
1174	#endif
1175	! !* Swap
1176	! ! ----
1177	ubb_e (:,:) = ub_e (:,:)
1178	ub_e (:,:) = un_e (:,:)
1179	un_e (:,:) = ua_e (:,:)
1180	!
1181	vbb_e (:,:) = vb_e (:,:)
1182	vb_e (:,:) = vn_e (:,:)
1183	vn_e (:,:) = va_e (:,:)
1184	!
1185	sshbb_e(:,:) = sshb_e(:,:)
1186	sshb_e (:,:) = sshn_e(:,:)
1187	sshn_e (:,:) = ssha_e(:,:)
1188
1189	! !* Sum over whole bt loop
1190	! ! ----------------------
1191	za1 = wgtbtp1(jn)
1192	IF( ln_dynadv_vec .OR. ln_linssh ) THEN ! Sum velocities
1193	ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:)
1194	va_b (:,:) = va_b (:,:) + za1 * va_e (:,:)
1195	ELSE ! Sum transports
1196	IF ( .NOT.ln_wd_dl ) THEN
1197	ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:) * hu_e (:,:)
1198	va_b (:,:) = va_b (:,:) + za1 * va_e (:,:) * hv_e (:,:)
1199	ELSE
1200	ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:) * hu_e (:,:) * zuwdmask(:,:)
1201	va_b (:,:) = va_b (:,:) + za1 * va_e (:,:) * hv_e (:,:) * zvwdmask(:,:)
1202	END IF
1203	ENDIF
1204	! ! Sum sea level
1205	ssha(:,:) = ssha(:,:) + za1 * ssha_e(:,:)
1206
1207	! ! ==================== !
1208	END DO ! end loop !
1209	! ! ==================== !
1210	! -----------------------------------------------------------------------------
1211	! Phase 3. update the general trend with the barotropic trend
1212	! -----------------------------------------------------------------------------
1213	!
1214	! Set advection velocity correction:
1215	IF (ln_bt_fw) THEN
1216	zwx(:,:) = un_adv(:,:)
1217	zwy(:,:) = vn_adv(:,:)
1218	IF( .NOT.( kt == nit000 .AND. neuler==0 ) ) THEN
1219	un_adv(:,:) = r1_2 * ( ub2_b(:,:) + zwx(:,:) - atfp * un_bf(:,:) )
1220	vn_adv(:,:) = r1_2 * ( vb2_b(:,:) + zwy(:,:) - atfp * vn_bf(:,:) )
1221	!
1222	! Update corrective fluxes for next time step:
1223	un_bf(:,:) = atfp * un_bf(:,:) + (zwx(:,:) - ub2_b(:,:))
1224	vn_bf(:,:) = atfp * vn_bf(:,:) + (zwy(:,:) - vb2_b(:,:))
1225	ELSE
1226	un_bf(:,:) = 0._wp
1227	vn_bf(:,:) = 0._wp
1228	END IF
1229	! Save integrated transport for next computation
1230	ub2_b(:,:) = zwx(:,:)
1231	vb2_b(:,:) = zwy(:,:)
1232	ENDIF
1233
1234
1235	!
1236	! Update barotropic trend:
1237	IF( ln_dynadv_vec .OR. ln_linssh ) THEN
1238	DO jk=1,jpkm1
1239	ua(:,:,jk) = ua(:,:,jk) + ( ua_b(:,:) - ub_b(:,:) ) * r1_2dt_b
1240	va(:,:,jk) = va(:,:,jk) + ( va_b(:,:) - vb_b(:,:) ) * r1_2dt_b
1241	END DO
1242	ELSE
1243	! At this stage, ssha has been corrected: compute new depths at velocity points
1244	DO jj = 1, jpjm1
1245	DO ji = 1, jpim1 ! NO Vector Opt.
1246	zsshu_a(ji,jj) = r1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj) &
1247	& * ( e1e2t(ji ,jj) * ssha(ji ,jj) &
1248	& + e1e2t(ji+1,jj) * ssha(ji+1,jj) )
1249	zsshv_a(ji,jj) = r1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj) &
1250	& * ( e1e2t(ji,jj ) * ssha(ji,jj ) &
1251	& + e1e2t(ji,jj+1) * ssha(ji,jj+1) )
1252	END DO
1253	END DO
1254	CALL lbc_lnk_multi( 'dynspg_ts', zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) ! Boundary conditions
1255	!
1256	DO jk=1,jpkm1
1257	ua(:,:,jk) = ua(:,:,jk) + r1_hu_n(:,:) * ( ua_b(:,:) - ub_b(:,:) * hu_b(:,:) ) * r1_2dt_b
1258	va(:,:,jk) = va(:,:,jk) + r1_hv_n(:,:) * ( va_b(:,:) - vb_b(:,:) * hv_b(:,:) ) * r1_2dt_b
1259	END DO
1260	! Save barotropic velocities not transport:
1261	ua_b(:,:) = ua_b(:,:) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - ssumask(:,:) )
1262	va_b(:,:) = va_b(:,:) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - ssvmask(:,:) )
1263	ENDIF
1264
1265
1266	! Correct velocities so that the barotropic velocity equals (un_adv, vn_adv) (in all cases)
1267	DO jk = 1, jpkm1
1268	un(:,:,jk) = ( un(:,:,jk) + un_adv(:,:)r1_hu_n(:,:) - un_b(:,:) ) umask(:,:,jk)
1269	vn(:,:,jk) = ( vn(:,:,jk) + vn_adv(:,:)r1_hv_n(:,:) - vn_b(:,:) ) vmask(:,:,jk)
1270	END DO
1271
1272	IF ( ln_wd_dl .and. ln_wd_dl_bc) THEN
1273	DO jk = 1, jpkm1
1274	un(:,:,jk) = ( un_adv(:,:)*r1_hu_n(:,:) &
1275	& + zuwdav2(:,:)(un(:,:,jk) - un_adv(:,:)r1_hu_n(:,:)) ) * umask(:,:,jk)
1276	vn(:,:,jk) = ( vn_adv(:,:)*r1_hv_n(:,:) &
1277	& + zvwdav2(:,:)(vn(:,:,jk) - vn_adv(:,:)r1_hv_n(:,:)) ) * vmask(:,:,jk)
1278	END DO
1279	END IF
1280
1281
1282	CALL iom_put( "ubar", un_adv(:,:)*r1_hu_n(:,:) ) ! barotropic i-current
1283	CALL iom_put( "vbar", vn_adv(:,:)*r1_hv_n(:,:) ) ! barotropic i-current
1284	!
1285	#if defined key_agrif
1286	! Save time integrated fluxes during child grid integration
1287	! (used to update coarse grid transports at next time step)
1288	!
1289	IF( .NOT.Agrif_Root() .AND. ln_bt_fw ) THEN
1290	IF( Agrif_NbStepint() == 0 ) THEN
1291	ub2_i_b(:,:) = 0._wp
1292	vb2_i_b(:,:) = 0._wp
1293	END IF
1294	!
1295	za1 = 1._wp / REAL(Agrif_rhot(), wp)
1296	ub2_i_b(:,:) = ub2_i_b(:,:) + za1 * ub2_b(:,:)
1297	vb2_i_b(:,:) = vb2_i_b(:,:) + za1 * vb2_b(:,:)
1298	ENDIF
1299	#endif
1300	! !* write time-spliting arrays in the restart
1301	IF( lrst_oce .AND.ln_bt_fw ) CALL ts_rst( kt, 'WRITE' )
1302	!
1303	IF( ln_wd_il ) DEALLOCATE( zcpx, zcpy )
1304	IF( ln_wd_dl ) DEALLOCATE( ztwdmask, zuwdmask, zvwdmask, zuwdav2, zvwdav2 )
1305	!
1306	IF( ln_diatmb ) THEN
1307	CALL iom_put( "baro_u" , un_bssumask(:,:)+zmdi(1.-ssumask(:,:) ) ) ! Barotropic U Velocity
1308	CALL iom_put( "baro_v" , vn_bssvmask(:,:)+zmdi(1.-ssvmask(:,:) ) ) ! Barotropic V Velocity
1309	ENDIF
1310	!
1311	END SUBROUTINE dyn_spg_ts
1312
1313
1314	SUBROUTINE ts_wgt( ll_av, ll_fw, jpit, zwgt1, zwgt2)
1315	!!---------------------------------------------------------------------
1316	!! * ROUTINE ts_wgt *
1317	!!
1318	!! ** Purpose : Set time-splitting weights for temporal averaging (or not)
1319	!!----------------------------------------------------------------------
1320	LOGICAL, INTENT(in) :: ll_av ! temporal averaging=.true.
1321	LOGICAL, INTENT(in) :: ll_fw ! forward time splitting =.true.
1322	INTEGER, INTENT(inout) :: jpit ! cycle length
1323	REAL(wp), DIMENSION(3*nn_baro), INTENT(inout) :: zwgt1, & ! Primary weights
1324	zwgt2 ! Secondary weights
1325
1326	INTEGER :: jic, jn, ji ! temporary integers
1327	REAL(wp) :: za1, za2
1328	!!----------------------------------------------------------------------
1329
1330	zwgt1(:) = 0._wp
1331	zwgt2(:) = 0._wp
1332
1333	! Set time index when averaged value is requested
1334	IF (ll_fw) THEN
1335	jic = nn_baro
1336	ELSE
1337	jic = 2 * nn_baro
1338	ENDIF
1339
1340	! Set primary weights:
1341	IF (ll_av) THEN
1342	! Define simple boxcar window for primary weights
1343	! (width = nn_baro, centered around jic)
1344	SELECT CASE ( nn_bt_flt )
1345	CASE( 0 ) ! No averaging
1346	zwgt1(jic) = 1._wp
1347	jpit = jic
1348
1349	CASE( 1 ) ! Boxcar, width = nn_baro
1350	DO jn = 1, 3*nn_baro
1351	za1 = ABS(float(jn-jic))/float(nn_baro)
1352	IF (za1 < 0.5_wp) THEN
1353	zwgt1(jn) = 1._wp
1354	jpit = jn
1355	ENDIF
1356	ENDDO
1357
1358	CASE( 2 ) ! Boxcar, width = 2 * nn_baro
1359	DO jn = 1, 3*nn_baro
1360	za1 = ABS(float(jn-jic))/float(nn_baro)
1361	IF (za1 < 1._wp) THEN
1362	zwgt1(jn) = 1._wp
1363	jpit = jn
1364	ENDIF
1365	ENDDO
1366	CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_bt_flt' )
1367	END SELECT
1368
1369	ELSE ! No time averaging
1370	zwgt1(jic) = 1._wp
1371	jpit = jic
1372	ENDIF
1373
1374	! Set secondary weights
1375	DO jn = 1, jpit
1376	DO ji = jn, jpit
1377	zwgt2(jn) = zwgt2(jn) + zwgt1(ji)
1378	END DO
1379	END DO
1380
1381	! Normalize weigths:
1382	za1 = 1._wp / SUM(zwgt1(1:jpit))
1383	za2 = 1._wp / SUM(zwgt2(1:jpit))
1384	DO jn = 1, jpit
1385	zwgt1(jn) = zwgt1(jn) * za1
1386	zwgt2(jn) = zwgt2(jn) * za2
1387	END DO
1388	!
1389	END SUBROUTINE ts_wgt
1390
1391
1392	SUBROUTINE ts_rst( kt, cdrw )
1393	!!---------------------------------------------------------------------
1394	!! * ROUTINE ts_rst *
1395	!!
1396	!! ** Purpose : Read or write time-splitting arrays in restart file
1397	!!----------------------------------------------------------------------
1398	INTEGER , INTENT(in) :: kt ! ocean time-step
1399	CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag
1400	!!----------------------------------------------------------------------
1401	!
1402	IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialise
1403	! ! ---------------
1404	IF( ln_rstart .AND. ln_bt_fw .AND. (neuler/=0) ) THEN !* Read the restart file
1405	CALL iom_get( numror, jpdom_autoglo, 'ub2_b' , ub2_b (:,:), ldxios = lrxios )
1406	CALL iom_get( numror, jpdom_autoglo, 'vb2_b' , vb2_b (:,:), ldxios = lrxios )
1407	CALL iom_get( numror, jpdom_autoglo, 'un_bf' , un_bf (:,:), ldxios = lrxios )
1408	CALL iom_get( numror, jpdom_autoglo, 'vn_bf' , vn_bf (:,:), ldxios = lrxios )
1409	IF( .NOT.ln_bt_av ) THEN
1410	CALL iom_get( numror, jpdom_autoglo, 'sshbb_e' , sshbb_e(:,:), ldxios = lrxios )
1411	CALL iom_get( numror, jpdom_autoglo, 'ubb_e' , ubb_e(:,:), ldxios = lrxios )
1412	CALL iom_get( numror, jpdom_autoglo, 'vbb_e' , vbb_e(:,:), ldxios = lrxios )
1413	CALL iom_get( numror, jpdom_autoglo, 'sshb_e' , sshb_e(:,:), ldxios = lrxios )
1414	CALL iom_get( numror, jpdom_autoglo, 'ub_e' , ub_e(:,:), ldxios = lrxios )
1415	CALL iom_get( numror, jpdom_autoglo, 'vb_e' , vb_e(:,:), ldxios = lrxios )
1416	ENDIF
1417	#if defined key_agrif
1418	! Read time integrated fluxes
1419	IF ( .NOT.Agrif_Root() ) THEN
1420	CALL iom_get( numror, jpdom_autoglo, 'ub2_i_b' , ub2_i_b(:,:), ldxios = lrxios )
1421	CALL iom_get( numror, jpdom_autoglo, 'vb2_i_b' , vb2_i_b(:,:), ldxios = lrxios )
1422	ENDIF
1423	#endif
1424	ELSE !* Start from rest
1425	IF(lwp) WRITE(numout,*)
1426	IF(lwp) WRITE(numout,*) ' ==>>> start from rest: set barotropic values to 0'
1427	ub2_b (:,:) = 0._wp ; vb2_b (:,:) = 0._wp ! used in the 1st interpol of agrif
1428	un_adv(:,:) = 0._wp ; vn_adv(:,:) = 0._wp ! used in the 1st interpol of agrif
1429	un_bf (:,:) = 0._wp ; vn_bf (:,:) = 0._wp ! used in the 1st update of agrif
1430	#if defined key_agrif
1431	IF ( .NOT.Agrif_Root() ) THEN
1432	ub2_i_b(:,:) = 0._wp ; vb2_i_b(:,:) = 0._wp ! used in the 1st update of agrif
1433	ENDIF
1434	#endif
1435	ENDIF
1436	!
1437	ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN ! Create restart file
1438	! ! -------------------
1439	IF(lwp) WRITE(numout,*) '---- ts_rst ----'
1440	IF( lwxios ) CALL iom_swap( cwxios_context )
1441	CALL iom_rstput( kt, nitrst, numrow, 'ub2_b' , ub2_b (:,:), ldxios = lwxios )
1442	CALL iom_rstput( kt, nitrst, numrow, 'vb2_b' , vb2_b (:,:), ldxios = lwxios )
1443	CALL iom_rstput( kt, nitrst, numrow, 'un_bf' , un_bf (:,:), ldxios = lwxios )
1444	CALL iom_rstput( kt, nitrst, numrow, 'vn_bf' , vn_bf (:,:), ldxios = lwxios )
1445	!
1446	IF (.NOT.ln_bt_av) THEN
1447	CALL iom_rstput( kt, nitrst, numrow, 'sshbb_e' , sshbb_e(:,:), ldxios = lwxios )
1448	CALL iom_rstput( kt, nitrst, numrow, 'ubb_e' , ubb_e(:,:), ldxios = lwxios )
1449	CALL iom_rstput( kt, nitrst, numrow, 'vbb_e' , vbb_e(:,:), ldxios = lwxios )
1450	CALL iom_rstput( kt, nitrst, numrow, 'sshb_e' , sshb_e(:,:), ldxios = lwxios )
1451	CALL iom_rstput( kt, nitrst, numrow, 'ub_e' , ub_e(:,:), ldxios = lwxios )
1452	CALL iom_rstput( kt, nitrst, numrow, 'vb_e' , vb_e(:,:), ldxios = lwxios )
1453	ENDIF
1454	#if defined key_agrif
1455	! Save time integrated fluxes
1456	IF ( .NOT.Agrif_Root() ) THEN
1457	CALL iom_rstput( kt, nitrst, numrow, 'ub2_i_b' , ub2_i_b(:,:), ldxios = lwxios )
1458	CALL iom_rstput( kt, nitrst, numrow, 'vb2_i_b' , vb2_i_b(:,:), ldxios = lwxios )
1459	ENDIF
1460	#endif
1461	IF( lwxios ) CALL iom_swap( cxios_context )
1462	ENDIF
1463	!
1464	END SUBROUTINE ts_rst
1465
1466
1467	SUBROUTINE dyn_spg_ts_init
1468	!!---------------------------------------------------------------------
1469	!! * ROUTINE dyn_spg_ts_init *
1470	!!
1471	!! ** Purpose : Set time splitting options
1472	!!----------------------------------------------------------------------
1473	INTEGER :: ji ,jj ! dummy loop indices
1474	REAL(wp) :: zxr2, zyr2, zcmax ! local scalar
1475	REAL(wp), DIMENSION(jpi,jpj) :: zcu
1476	!!----------------------------------------------------------------------
1477	!
1478	! Max courant number for ext. grav. waves
1479	!
1480	DO jj = 1, jpj
1481	DO ji =1, jpi
1482	zxr2 = r1_e1t(ji,jj) * r1_e1t(ji,jj)
1483	zyr2 = r1_e2t(ji,jj) * r1_e2t(ji,jj)
1484	zcu(ji,jj) = SQRT( grav * MAX(ht_0(ji,jj),0._wp) * (zxr2 + zyr2) )
1485	END DO
1486	END DO
1487	!
1488	zcmax = MAXVAL( zcu(:,:) )
1489	CALL mpp_max( 'dynspg_ts', zcmax )
1490
1491	! Estimate number of iterations to satisfy a max courant number= rn_bt_cmax
1492	IF( ln_bt_auto ) nn_baro = CEILING( rdt / rn_bt_cmax * zcmax)
1493
1494	rdtbt = rdt / REAL( nn_baro , wp )
1495	zcmax = zcmax * rdtbt
1496	! Print results
1497	IF(lwp) WRITE(numout,*)
1498	IF(lwp) WRITE(numout,*) 'dyn_spg_ts_init : split-explicit free surface'
1499	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~'
1500	IF( ln_bt_auto ) THEN
1501	IF(lwp) WRITE(numout,*) ' ln_ts_auto =.true. Automatically set nn_baro '
1502	IF(lwp) WRITE(numout,*) ' Max. courant number allowed: ', rn_bt_cmax
1503	ELSE
1504	IF(lwp) WRITE(numout,*) ' ln_ts_auto=.false.: Use nn_baro in namelist nn_baro = ', nn_baro
1505	ENDIF
1506
1507	IF(ln_bt_av) THEN
1508	IF(lwp) WRITE(numout,*) ' ln_bt_av =.true. ==> Time averaging over nn_baro time steps is on '
1509	ELSE
1510	IF(lwp) WRITE(numout,*) ' ln_bt_av =.false. => No time averaging of barotropic variables '
1511	ENDIF
1512	!
1513	!
1514	IF(ln_bt_fw) THEN
1515	IF(lwp) WRITE(numout,*) ' ln_bt_fw=.true. => Forward integration of barotropic variables '
1516	ELSE
1517	IF(lwp) WRITE(numout,*) ' ln_bt_fw =.false.=> Centred integration of barotropic variables '
1518	ENDIF
1519	!
1520	#if defined key_agrif
1521	! Restrict the use of Agrif to the forward case only
1522	!!! IF( .NOT.ln_bt_fw .AND. .NOT.Agrif_Root() ) CALL ctl_stop( 'AGRIF not implemented if ln_bt_fw=.FALSE.' )
1523	#endif
1524	!
1525	IF(lwp) WRITE(numout,*) ' Time filter choice, nn_bt_flt: ', nn_bt_flt
1526	SELECT CASE ( nn_bt_flt )
1527	CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' Dirac'
1528	CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' Boxcar: width = nn_baro'
1529	CASE( 2 ) ; IF(lwp) WRITE(numout,) ' Boxcar: width = 2nn_baro'
1530	CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_bt_flt: should 0,1, or 2' )
1531	END SELECT
1532	!
1533	IF(lwp) WRITE(numout,*) ' '
1534	IF(lwp) WRITE(numout,*) ' nn_baro = ', nn_baro
1535	IF(lwp) WRITE(numout,*) ' Barotropic time step [s] is :', rdtbt
1536	IF(lwp) WRITE(numout,*) ' Maximum Courant number is :', zcmax
1537	!
1538	IF(lwp) WRITE(numout,*) ' Time diffusion parameter rn_bt_alpha: ', rn_bt_alpha
1539	IF ((ln_bt_av.AND.nn_bt_flt/=0).AND.(rn_bt_alpha>0._wp)) THEN
1540	CALL ctl_stop( 'dynspg_ts ERROR: if rn_bt_alpha > 0, remove temporal averaging' )
1541	ENDIF
1542	!
1543	IF( .NOT.ln_bt_av .AND. .NOT.ln_bt_fw ) THEN
1544	CALL ctl_stop( 'dynspg_ts ERROR: No time averaging => only forward integration is possible' )
1545	ENDIF
1546	IF( zcmax>0.9_wp ) THEN
1547	CALL ctl_stop( 'dynspg_ts ERROR: Maximum Courant number is greater than 0.9: Inc. nn_baro !' )
1548	ENDIF
1549	!
1550	! ! Allocate time-splitting arrays
1551	IF( dyn_spg_ts_alloc() /= 0 ) CALL ctl_stop('STOP', 'dyn_spg_init: failed to allocate dynspg_ts arrays' )
1552	!
1553	! ! read restart when needed
1554	CALL ts_rst( nit000, 'READ' )
1555	!
1556	IF( lwxios ) THEN
1557	! define variables in restart file when writing with XIOS
1558	CALL iom_set_rstw_var_active('ub2_b')
1559	CALL iom_set_rstw_var_active('vb2_b')
1560	CALL iom_set_rstw_var_active('un_bf')
1561	CALL iom_set_rstw_var_active('vn_bf')
1562	!
1563	IF (.NOT.ln_bt_av) THEN
1564	CALL iom_set_rstw_var_active('sshbb_e')
1565	CALL iom_set_rstw_var_active('ubb_e')
1566	CALL iom_set_rstw_var_active('vbb_e')
1567	CALL iom_set_rstw_var_active('sshb_e')
1568	CALL iom_set_rstw_var_active('ub_e')
1569	CALL iom_set_rstw_var_active('vb_e')
1570	ENDIF
1571	#if defined key_agrif
1572	! Save time integrated fluxes
1573	IF ( .NOT.Agrif_Root() ) THEN
1574	CALL iom_set_rstw_var_active('ub2_i_b')
1575	CALL iom_set_rstw_var_active('vb2_i_b')
1576	ENDIF
1577	#endif
1578	ENDIF
1579	!
1580	END SUBROUTINE dyn_spg_ts_init
1581
1582	!!======================================================================
1583	END MODULE dynspg_ts

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

Download in other formats: