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

Last change on this file since 11529 was 11529, checked in by mathiot, 5 years ago
ENHANCE-02_ISF: fixes for ice sheet coupling (ticket #2142)
Property svn:keywords set to `Id`
File size: 77.6 KB

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

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

Download in other formats: