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 @ 11541

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

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

Download in other formats: