Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

dynspg_ts.F90 @ 10094

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

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

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

Context Navigation

source: NEMO/branches/UKMO/dev_r10037_dynvor_EEUV/src/OCE/DYN/dynspg_ts.F90 @ 10094

Download in other formats: