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 branches/UKMO/ROMS_WAD_7832/NEMOGCM/CONFIG/CS15mini/MY_SRC – NEMO

Context Navigation

source: branches/UKMO/ROMS_WAD_7832/NEMOGCM/CONFIG/CS15mini/MY_SRC/dynspg_ts.F90 @ 8544

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

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

Download in other formats: