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/NEMO/OPA_SRC/DYN – NEMO

Context Navigation

source: branches/UKMO/ROMS_WAD_7832/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 8865

Last change on this file since 8865 was 8865, checked in by deazer, 6 years ago

Moving Changes from CS15mini config into NEMO main src
Updating TEST configs to run wit this version of the code
all sette tests pass at this revision other than AGRIF
Includes updates to dynnxt and tranxt required for 3D rives in WAD case to be conservative.

Next commit will update naming conventions and tidy the code.

Property svn:keywords set to Id

File size: 72.4 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$
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	! un_e and vn_e are set to zero at faces where the direction of the flow is from dry cells
785
786	DO jj = 1, jpjm1
787	DO ji = 1, jpim1
788	IF ( zwx(ji,jj) > 0.0 ) THEN
789	zuwdmask(ji, jj) = ztwdmask(ji ,jj)
790	ELSE
791	zuwdmask(ji, jj) = ztwdmask(ji+1,jj)
792	END IF
793	zwx(ji, jj) = zuwdmask(ji,jj)*zwx(ji, jj)
794	un_e(ji,jj) = zuwdmask(ji,jj)*un_e(ji,jj)
795
796	IF ( zwy(ji,jj) > 0.0 ) THEN
797	zvwdmask(ji, jj) = ztwdmask(ji, jj )
798	ELSE
799	zvwdmask(ji, jj) = ztwdmask(ji, jj+1)
800	END IF
801	zwy(ji, jj) = zvwdmask(ji,jj)*zwy(ji,jj)
802	vn_e(ji,jj) = zvwdmask(ji,jj)*vn_e(ji,jj)
803	END DO
804	END DO
805
806	IF ( ln_wd_diag ) THEN
807	WRITE(numout, *) 'zuwdmask: (i,j) = ', zuwdmask(iwdg,jwdg)
808	WRITE(numout, *) 'zwx: (i,j) = ', zwx(iwdg,jwdg)
809	WRITE(numout, *) 'e2u: (i,j) = ', e2u(iwdg,jwdg)
810	WRITE(numout, *) 'ua_e: (i,j) = ', ua_e(iwdg,jwdg)
811	WRITE(numout, *) 'un_e: (i,j) = ', un_e(iwdg,jwdg)
812	WRITE(numout, *) 'zhup2_e: (i,j) = ', zhup2_e(iwdg,jwdg)
813	WRITE(numout, *) 'zvwdmask: (i,j) = ', zvwdmask(iwdg,jwdg)
814	WRITE(numout, *) 'zwy: (i,j) = ', zwy(iwdg,jwdg)
815	END IF
816
817	END IF
818
819	! Sum over sub-time-steps to compute advective velocities
820	za2 = wgtbtp2(jn)
821	un_adv(:,:) = un_adv(:,:) + za2 * zwx(:,:) * r1_e2u(:,:)
822	vn_adv(:,:) = vn_adv(:,:) + za2 * zwy(:,:) * r1_e1v(:,:)
823
824	! sum over sub-time-steps to decide which baroclinic velocities to set to zero (zuwdav2 is only used when ln_rwd_bc = True)
825	IF ( ln_rwd_bc ) THEN
826	zuwdav2(:,:) = zuwdav2(:,:) + za2 * zuwdmask(:,:)
827	zvwdav2(:,:) = zvwdav2(:,:) + za2 * zvwdmask(:,:)
828
829	IF ( ln_wd_diag ) THEN
830	WRITE(numout, *) 'za2, r1_e2u(i,j) = ', za2, r1_e2u(iwdg,jwdg)
831	WRITE(numout, *) 'un_adv: (i,j) = ', un_adv(iwdg,jwdg)
832	WRITE(numout, *) 'zuwdav2: (i,j) = ', zuwdav2(iwdg,jwdg)
833	WRITE(numout, *) 'zvwdav2: (i,j) = ', zvwdav2(iwdg,jwdg)
834	END IF
835
836	END IF
837
838	! Set next sea level:
839	DO jj = 2, jpjm1
840	DO ji = fs_2, fs_jpim1 ! vector opt.
841	zhdiv(ji,jj) = ( zwx(ji,jj) - zwx(ji-1,jj) &
842	& + zwy(ji,jj) - zwy(ji,jj-1) ) * r1_e1e2t(ji,jj)
843	END DO
844	END DO
845	ssha_e(:,:) = ( sshn_e(:,:) - rdtbt * ( zssh_frc(:,:) + zhdiv(:,:) ) ) * ssmask(:,:)
846
847	CALL lbc_lnk( ssha_e, 'T', 1._wp )
848
849	! Duplicate sea level across open boundaries (this is only cosmetic if linssh=T)
850	IF( ln_bdy ) CALL bdy_ssh( ssha_e )
851	#if defined key_agrif
852	IF( .NOT.Agrif_Root() ) CALL agrif_ssh_ts( jn )
853	#endif
854	!
855	! Sea Surface Height at u-,v-points (vvl case only)
856	IF( .NOT.ln_linssh ) THEN
857	DO jj = 2, jpjm1
858	DO ji = 2, jpim1 ! NO Vector Opt.
859	zsshu_a(ji,jj) = z1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj) &
860	& * ( e1e2t(ji ,jj ) * ssha_e(ji ,jj ) &
861	& + e1e2t(ji+1,jj ) * ssha_e(ji+1,jj ) )
862	zsshv_a(ji,jj) = z1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj) &
863	& * ( e1e2t(ji ,jj ) * ssha_e(ji ,jj ) &
864	& + e1e2t(ji ,jj+1) * ssha_e(ji ,jj+1) )
865	END DO
866	END DO
867	CALL lbc_lnk_multi( zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp )
868	ENDIF
869	!
870	! Half-step back interpolation of SSH for surface pressure computation:
871	!----------------------------------------------------------------------
872	IF ((jn==1).AND.ll_init) THEN
873	za0=1._wp ! Forward-backward
874	za1=0._wp
875	za2=0._wp
876	za3=0._wp
877	ELSEIF ((jn==2).AND.ll_init) THEN ! AB2-AM3 Coefficients; bet=0 ; gam=-1/6 ; eps=1/12
878	za0= 1.0833333333333_wp ! za0 = 1-gam-eps
879	za1=-0.1666666666666_wp ! za1 = gam
880	za2= 0.0833333333333_wp ! za2 = eps
881	za3= 0._wp
882	ELSE ! AB3-AM4 Coefficients; bet=0.281105 ; eps=0.013 ; gam=0.0880
883	za0=0.614_wp ! za0 = 1/2 + gam + 2*eps
884	za1=0.285_wp ! za1 = 1/2 - 2gam - 3eps
885	za2=0.088_wp ! za2 = gam
886	za3=0.013_wp ! za3 = eps
887	ENDIF
888	!
889	zsshp2_e(:,:) = za0 * ssha_e(:,:) + za1 * sshn_e (:,:) &
890	& + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:)
891	IF( ln_wd ) THEN ! Calculating and applying W/D gravity filters
892	DO jj = 2, jpjm1
893	DO ji = 2, jpim1
894	ll_tmp1 = MIN( zsshp2_e(ji,jj) , zsshp2_e(ji+1,jj) ) > &
895	& MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) .AND. &
896	& MAX( zsshp2_e(ji,jj) + ht_0(ji,jj), zsshp2_e(ji+1,jj) + ht_0(ji+1,jj) ) &
897	& > rn_wdmin1 + rn_wdmin2
898	ll_tmp2 = (ABS(zsshp2_e(ji,jj) - zsshp2_e(ji+1,jj)) > 1.E-12 ).AND.( &
899	& MAX( zsshp2_e(ji,jj) , zsshp2_e(ji+1,jj) ) > &
900	& MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 )
901
902	IF(ll_tmp1) THEN
903	zcpx(ji,jj) = 1.0_wp
904	ELSE IF(ll_tmp2) THEN
905	! no worries about zsshp2_e(ji+1,jj) - zsshp2_e(ji ,jj) = 0, it won't happen ! here
906	zcpx(ji,jj) = ABS( (zsshp2_e(ji+1,jj) + ht_0(ji+1,jj) - zsshp2_e(ji,jj) - ht_0(ji,jj)) &
907	& / (zsshp2_e(ji+1,jj) - zsshp2_e(ji ,jj)) )
908	ELSE
909	zcpx(ji,jj) = 0._wp
910	END IF
911
912	ll_tmp1 = MIN( zsshp2_e(ji,jj) , zsshp2_e(ji,jj+1) ) > &
913	& MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) .AND. &
914	& MAX( zsshp2_e(ji,jj) + ht_0(ji,jj), zsshp2_e(ji,jj+1) + ht_0(ji,jj+1) ) &
915	& > rn_wdmin1 + rn_wdmin2
916	ll_tmp2 = (ABS(zsshp2_e(ji,jj) - zsshp2_e(ji,jj+1)) > 1.E-12 ).AND.( &
917	& MAX( zsshp2_e(ji,jj) , zsshp2_e(ji,jj+1) ) > &
918	& MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 )
919
920	IF(ll_tmp1) THEN
921	zcpy(ji,jj) = 1.0_wp
922	ELSE IF(ll_tmp2) THEN
923	! no worries about zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj ) = 0, it won't happen ! here
924	zcpy(ji,jj) = ABS( (zsshp2_e(ji,jj+1) + ht_0(ji,jj+1) - zsshp2_e(ji,jj) - ht_0(ji,jj)) &
925	& / (zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj )) )
926	ELSE
927	zcpy(ji,jj) = 0._wp
928	END IF
929	END DO
930	END DO
931	END IF
932	!
933	! Compute associated depths at U and V points:
934	IF( .NOT.ln_linssh .AND. .NOT.ln_dynadv_vec ) THEN !* Vector form
935	!
936	DO jj = 2, jpjm1
937	DO ji = 2, jpim1
938	zx1 = z1_2 * ssumask(ji ,jj) * r1_e1e2u(ji ,jj) &
939	& * ( e1e2t(ji ,jj ) * zsshp2_e(ji ,jj) &
940	& + e1e2t(ji+1,jj ) * zsshp2_e(ji+1,jj ) )
941	zy1 = z1_2 * ssvmask(ji ,jj) * r1_e1e2v(ji ,jj ) &
942	& * ( e1e2t(ji ,jj ) * zsshp2_e(ji ,jj ) &
943	& + e1e2t(ji ,jj+1) * zsshp2_e(ji ,jj+1) )
944	zhust_e(ji,jj) = hu_0(ji,jj) + zx1
945	zhvst_e(ji,jj) = hv_0(ji,jj) + zy1
946	END DO
947	END DO
948
949	ENDIF
950	!
951	! Add Coriolis trend:
952	! zwz array below or triads normally depend on sea level with ln_linssh=F and should be updated
953	! at each time step. We however keep them constant here for optimization.
954	! Recall that zwx and zwy arrays hold fluxes at this stage:
955	! zwx(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:) ! fluxes at jn+0.5
956	! zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:)
957	!
958	IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN !== energy conserving or mixed scheme ==!
959	DO jj = 2, jpjm1
960	DO ji = fs_2, fs_jpim1 ! vector opt.
961	zy1 = ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) ) * r1_e1u(ji,jj)
962	zy2 = ( zwy(ji ,jj ) + zwy(ji+1,jj ) ) * r1_e1u(ji,jj)
963	zx1 = ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) ) * r1_e2v(ji,jj)
964	zx2 = ( zwx(ji ,jj ) + zwx(ji ,jj+1) ) * r1_e2v(ji,jj)
965	zu_trd(ji,jj) = z1_4 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
966	zv_trd(ji,jj) =-z1_4 * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 )
967	END DO
968	END DO
969	!
970	ELSEIF ( ln_dynvor_ens ) THEN !== enstrophy conserving scheme ==!
971	DO jj = 2, jpjm1
972	DO ji = fs_2, fs_jpim1 ! vector opt.
973	zy1 = z1_8 * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) &
974	& + zwy(ji ,jj ) + zwy(ji+1,jj ) ) * r1_e1u(ji,jj)
975	zx1 = - z1_8 * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) &
976	& + zwx(ji ,jj ) + zwx(ji ,jj+1) ) * r1_e2v(ji,jj)
977	zu_trd(ji,jj) = zy1 * ( zwz(ji ,jj-1) + zwz(ji,jj) )
978	zv_trd(ji,jj) = zx1 * ( zwz(ji-1,jj ) + zwz(ji,jj) )
979	END DO
980	END DO
981	!
982	ELSEIF ( ln_dynvor_een ) THEN !== energy and enstrophy conserving scheme ==!
983	DO jj = 2, jpjm1
984	DO ji = fs_2, fs_jpim1 ! vector opt.
985	zu_trd(ji,jj) = + z1_12 * r1_e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) &
986	& + ftnw(ji+1,jj) * zwy(ji+1,jj ) &
987	& + ftse(ji,jj ) * zwy(ji ,jj-1) &
988	& + ftsw(ji+1,jj) * zwy(ji+1,jj-1) )
989	zv_trd(ji,jj) = - z1_12 * r1_e2v(ji,jj) * ( ftsw(ji,jj+1) * zwx(ji-1,jj+1) &
990	& + ftse(ji,jj+1) * zwx(ji ,jj+1) &
991	& + ftnw(ji,jj ) * zwx(ji-1,jj ) &
992	& + ftne(ji,jj ) * zwx(ji ,jj ) )
993	END DO
994	END DO
995	!
996	ENDIF
997	!
998	! Add tidal astronomical forcing if defined
999	IF ( ln_tide .AND. ln_tide_pot ) THEN
1000	DO jj = 2, jpjm1
1001	DO ji = fs_2, fs_jpim1 ! vector opt.
1002	zu_spg = grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) * r1_e1u(ji,jj)
1003	zv_spg = grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) * r1_e2v(ji,jj)
1004	zu_trd(ji,jj) = zu_trd(ji,jj) + zu_spg
1005	zv_trd(ji,jj) = zv_trd(ji,jj) + zv_spg
1006	END DO
1007	END DO
1008	ENDIF
1009	!
1010	! Add bottom stresses:
1011	!jth do implicitly instead
1012	! zu_trd(:,:) = zu_trd(:,:) + bfrua(:,:) * un_e(:,:) * hur_e(:,:)
1013	! zv_trd(:,:) = zv_trd(:,:) + bfrva(:,:) * vn_e(:,:) * hvr_e(:,:)
1014	!
1015	! Add top stresses:
1016	zu_trd(:,:) = zu_trd(:,:) + tfrua(:,:) * un_e(:,:) * hur_e(:,:)
1017	zv_trd(:,:) = zv_trd(:,:) + tfrva(:,:) * vn_e(:,:) * hvr_e(:,:)
1018	!
1019	! Surface pressure trend:
1020
1021	IF( ln_wd ) THEN
1022	DO jj = 2, jpjm1
1023	DO ji = 2, jpim1
1024	! Add surface pressure gradient
1025	zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj)
1026	zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj)
1027	zwx(ji,jj) = zu_spg * zcpx(ji,jj)
1028	zwy(ji,jj) = zv_spg * zcpy(ji,jj)
1029	END DO
1030	END DO
1031	ELSE
1032	DO jj = 2, jpjm1
1033	DO ji = fs_2, fs_jpim1 ! vector opt.
1034	! Add surface pressure gradient
1035	zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj)
1036	zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj)
1037	zwx(ji,jj) = zu_spg
1038	zwy(ji,jj) = zv_spg
1039	END DO
1040	END DO
1041	END IF
1042
1043	!
1044	! Set next velocities:
1045	IF( ln_dynadv_vec .OR. ln_linssh ) THEN !* Vector form
1046	DO jj = 2, jpjm1
1047	DO ji = fs_2, fs_jpim1 ! vector opt.
1048	ua_e(ji,jj) = ( un_e(ji,jj) &
1049	& + rdtbt * ( zwx(ji,jj) &
1050	& + zu_trd(ji,jj) &
1051	& + zu_frc(ji,jj) ) &
1052	& ) * ssumask(ji,jj)
1053
1054	va_e(ji,jj) = ( vn_e(ji,jj) &
1055	& + rdtbt * ( zwy(ji,jj) &
1056	& + zv_trd(ji,jj) &
1057	& + zv_frc(ji,jj) ) &
1058	& ) * ssvmask(ji,jj)
1059
1060	!jth implicit bottom friction:
1061	ua_e(ji,jj) = ua_e(ji,jj) /(1.0 - rdtbt * bfrua(ji,jj) * hur_e(ji,jj))
1062	va_e(ji,jj) = va_e(ji,jj) /(1.0 - rdtbt * bfrva(ji,jj) * hvr_e(ji,jj))
1063
1064	END DO
1065	END DO
1066	!
1067	ELSE !* Flux form
1068	DO jj = 2, jpjm1
1069	DO ji = fs_2, fs_jpim1 ! vector opt.
1070
1071	IF( ln_wd ) THEN
1072	zhura = hu_0(ji,jj) + zsshu_a(ji,jj)
1073	zhvra = hv_0(ji,jj) + zsshv_a(ji,jj)
1074	ELSE
1075	zhura = hu_0(ji,jj) + zsshu_a(ji,jj)
1076	zhvra = hv_0(ji,jj) + zsshv_a(ji,jj)
1077	END IF
1078	zhura = ssumask(ji,jj)/(zhura + 1._wp - ssumask(ji,jj))
1079	zhvra = ssvmask(ji,jj)/(zhvra + 1._wp - ssvmask(ji,jj))
1080
1081	ua_e(ji,jj) = ( hu_e(ji,jj) * un_e(ji,jj) &
1082	& + rdtbt * ( zhust_e(ji,jj) * zwx(ji,jj) &
1083	& + zhup2_e(ji,jj) * zu_trd(ji,jj) &
1084	& + hu_n(ji,jj) * zu_frc(ji,jj) ) &
1085	& ) * zhura
1086
1087	va_e(ji,jj) = ( hv_e(ji,jj) * vn_e(ji,jj) &
1088	& + rdtbt * ( zhvst_e(ji,jj) * zwy(ji,jj) &
1089	& + zhvp2_e(ji,jj) * zv_trd(ji,jj) &
1090	& + hv_n(ji,jj) * zv_frc(ji,jj) ) &
1091	& ) * zhvra
1092	END DO
1093	END DO
1094	ENDIF
1095
1096	! if ln_rwd: ua_e and va_e should not be masked ; they are used to determine the direction of flow into all cells
1097
1098	! IF ( ln_rwd) THEN
1099	! IF ( ln_wd_diag ) THEN
1100	! WRITE(numout, *) 'ua_e: (i,j) = ', ua_e(iwdg,jwdg)
1101	! WRITE(numout, *) 'va_e: (i,j) = ', va_e(iwdg,jwdg)
1102	! END IF
1103	! ua_e(:,:) = ua_e(:,:) * zuwdmask(:,:)
1104	! va_e(:,:) = va_e(:,:) * zvwdmask(:,:)
1105	! IF ( ln_wd_diag ) THEN
1106	! WRITE(numout, *) 'ua_e: (i,j) = ', ua_e(iwdg,jwdg)
1107	! WRITE(numout, *) 'va_e: (i,j) = ', va_e(iwdg,jwdg)
1108	! END IF
1109	! END IF
1110
1111
1112	IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only)
1113	IF( ln_wd ) THEN
1114	hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:)
1115	hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:)
1116	ELSE
1117	hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:)
1118	hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:)
1119	END IF
1120	hur_e(:,:) = ssumask(:,:) / ( hu_e(:,:) + 1._wp - ssumask(:,:) )
1121	hvr_e(:,:) = ssvmask(:,:) / ( hv_e(:,:) + 1._wp - ssvmask(:,:) )
1122	!
1123	ENDIF
1124	! !* domain lateral boundary
1125	CALL lbc_lnk_multi( ua_e, 'U', -1._wp, va_e , 'V', -1._wp )
1126	!
1127	! ! open boundaries
1128	IF( ln_bdy ) CALL bdy_dyn2d( jn, ua_e, va_e, un_e, vn_e, hur_e, hvr_e, ssha_e )
1129	#if defined key_agrif
1130	IF( .NOT.Agrif_Root() ) CALL agrif_dyn_ts( jn ) ! Agrif
1131	#endif
1132	! !* Swap
1133	! ! ----
1134	ubb_e (:,:) = ub_e (:,:)
1135	ub_e (:,:) = un_e (:,:)
1136	un_e (:,:) = ua_e (:,:)
1137	!
1138	vbb_e (:,:) = vb_e (:,:)
1139	vb_e (:,:) = vn_e (:,:)
1140	vn_e (:,:) = va_e (:,:)
1141	!
1142	sshbb_e(:,:) = sshb_e(:,:)
1143	sshb_e (:,:) = sshn_e(:,:)
1144	sshn_e (:,:) = ssha_e(:,:)
1145
1146	! !* Sum over whole bt loop
1147	! ! ----------------------
1148	za1 = wgtbtp1(jn)
1149	IF( ln_dynadv_vec .OR. ln_linssh ) THEN ! Sum velocities
1150	ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:)
1151	va_b (:,:) = va_b (:,:) + za1 * va_e (:,:)
1152	ELSE ! Sum transports
1153	IF (.NOT.ln_rwd) THEN
1154	ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:) * hu_e (:,:)
1155	va_b (:,:) = va_b (:,:) + za1 * va_e (:,:) * hv_e (:,:)
1156	ELSE
1157	ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:) * hu_e (:,:) * zuwdmask(:,:)
1158	va_b (:,:) = va_b (:,:) + za1 * va_e (:,:) * hv_e (:,:) * zvwdmask(:,:)
1159	END IF
1160	ENDIF
1161	! ! Sum sea level
1162	ssha(:,:) = ssha(:,:) + za1 * ssha_e(:,:)
1163
1164	! ! ==================== !
1165	END DO ! end loop !
1166	! ! ==================== !
1167	! -----------------------------------------------------------------------------
1168	! Phase 3. update the general trend with the barotropic trend
1169	! -----------------------------------------------------------------------------
1170	!
1171	! Set advection velocity correction:
1172	zwx(:,:) = un_adv(:,:)
1173	zwy(:,:) = vn_adv(:,:)
1174	IF( ( kt == nit000 .AND. neuler==0 ) .OR. .NOT.ln_bt_fw ) THEN
1175	un_adv(:,:) = zwx(:,:) * r1_hu_n(:,:)
1176	vn_adv(:,:) = zwy(:,:) * r1_hv_n(:,:)
1177	ELSE
1178	un_adv(:,:) = z1_2 * ( ub2_b(:,:) + zwx(:,:) ) * r1_hu_n(:,:)
1179	vn_adv(:,:) = z1_2 * ( vb2_b(:,:) + zwy(:,:) ) * r1_hv_n(:,:)
1180	END IF
1181
1182	IF( ln_bt_fw ) THEN ! Save integrated transport for next computation
1183	ub2_b(:,:) = zwx(:,:)
1184	vb2_b(:,:) = zwy(:,:)
1185	ENDIF
1186
1187	IF ( ln_wd_diag ) THEN
1188	WRITE(numout, *) 'ub2_b: (i,j) = ', ub2_b(iwdg,jwdg)
1189	WRITE(numout, *) 'r1_hu_n: (i,j) = ', r1_hu_n(iwdg,jwdg)
1190	WRITE(numout, *) 'zwx: (i,j) = ', zwx(iwdg,jwdg)
1191	WRITE(numout, *) 'un_adv: (i,j) = ', un_adv(iwdg,jwdg)
1192	END IF
1193
1194	!
1195	! Update barotropic trend:
1196	IF( ln_dynadv_vec .OR. ln_linssh ) THEN
1197	DO jk=1,jpkm1
1198	ua(:,:,jk) = ua(:,:,jk) + ( ua_b(:,:) - ub_b(:,:) ) * z1_2dt_b
1199	va(:,:,jk) = va(:,:,jk) + ( va_b(:,:) - vb_b(:,:) ) * z1_2dt_b
1200	END DO
1201	ELSE
1202	! At this stage, ssha has been corrected: compute new depths at velocity points
1203	DO jj = 1, jpjm1
1204	DO ji = 1, jpim1 ! NO Vector Opt.
1205	zsshu_a(ji,jj) = z1_2 * umask(ji,jj,1) * r1_e1e2u(ji,jj) &
1206	& * ( e1e2t(ji ,jj) * ssha(ji ,jj) &
1207	& + e1e2t(ji+1,jj) * ssha(ji+1,jj) )
1208	zsshv_a(ji,jj) = z1_2 * vmask(ji,jj,1) * r1_e1e2v(ji,jj) &
1209	& * ( e1e2t(ji,jj ) * ssha(ji,jj ) &
1210	& + e1e2t(ji,jj+1) * ssha(ji,jj+1) )
1211	END DO
1212	END DO
1213	CALL lbc_lnk_multi( zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) ! Boundary conditions
1214	!
1215	DO jk=1,jpkm1
1216	ua(:,:,jk) = ua(:,:,jk) + r1_hu_n(:,:) * ( ua_b(:,:) - ub_b(:,:) * hu_b(:,:) ) * z1_2dt_b
1217	va(:,:,jk) = va(:,:,jk) + r1_hv_n(:,:) * ( va_b(:,:) - vb_b(:,:) * hv_b(:,:) ) * z1_2dt_b
1218	END DO
1219	! Save barotropic velocities not transport:
1220	ua_b(:,:) = ua_b(:,:) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - ssumask(:,:) )
1221	va_b(:,:) = va_b(:,:) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - ssvmask(:,:) )
1222	ENDIF
1223
1224	IF ( ln_wd_diag ) THEN
1225	WRITE(numout, *) 'ua_b: (i,j) A = ', ua_b(iwdg,jwdg)
1226	WRITE(numout, *) 'va_b: (i,j) B = ', va_b(iwdg,jwdg)
1227	END IF
1228
1229	! temporary debugging code
1230	IF ( ln_wd_diag ) THEN
1231	WRITE(numout, *) 'ua: (i,j,k) B = ', ua(iwdg,jwdg,kwdg)
1232	WRITE(numout, *) 'ua_b: (i,j) B = ', ua_b(iwdg,jwdg)
1233	WRITE(numout, *) 'un: (i,j,k) = ', un(iwdg,jwdg,kwdg)
1234	WRITE(numout, *) 'un_b: (i,j) = ', un_b(iwdg,jwdg)
1235	WRITE(numout, *) 'un_adv: (i,j) = ', un_adv(iwdg,jwdg)
1236	WRITE(numout, *) 'va: (i,j,k) = ', va(iwdg,jwdg,kwdg)
1237	WRITE(numout, *) 'va_b: (i,j,k) = ', va_b(iwdg,jwdg)
1238	WRITE(numout, *) 'vn: (i,j,k) = ', vn(iwdg,jwdg,kwdg)
1239	WRITE(numout, *) 'vn_b: (i,j) = ', vn_b(iwdg,jwdg)
1240	WRITE(numout, *) 'vn_adv: (i,j) = ', vn_adv(iwdg,jwdg)
1241	END IF
1242
1243	! Correct velocities so that the barotropic velocity equals (un_adv, vn_adv) (in all cases)
1244	DO jk = 1, jpkm1
1245	un(:,:,jk) = ( un(:,:,jk) + un_adv(:,:) - un_b(:,:) ) * umask(:,:,jk)
1246	vn(:,:,jk) = ( vn(:,:,jk) + vn_adv(:,:) - vn_b(:,:) ) * vmask(:,:,jk)
1247	END DO
1248
1249	IF ( ln_rwd .and. ln_rwd_bc) THEN
1250	DO jk = 1, jpkm1
1251	un(:,:,jk) = ( un_adv(:,:) + zuwdav2(:,:)(un(:,:,jk) - un_adv(:,:)) ) umask(:,:,jk)
1252	vn(:,:,jk) = ( vn_adv(:,:) + zvwdav2(:,:)(vn(:,:,jk) - vn_adv(:,:)) ) vmask(:,:,jk)
1253	END DO
1254	END IF
1255
1256	IF ( ln_wd_diag ) THEN
1257	WRITE(numout, *) 'ua: (i,j,k) = ', ua(iwdg,jwdg,kwdg)
1258	WRITE(numout, *) 'ua_b: (i,j,k) = ', ua_b(iwdg,jwdg)
1259	WRITE(numout, *) 'un: (i,j,k) = ', un(iwdg,jwdg,kwdg)
1260	WRITE(numout, *) 'va: (i,j,k) = ', va(iwdg,jwdg,kwdg)
1261	WRITE(numout, *) 'va_b: (i,j,k) = ', va_b(iwdg,jwdg)
1262	WRITE(numout, *) 'vn: (i,j,k) = ', vn(iwdg,jwdg,kwdg)
1263	END IF
1264
1265	CALL iom_put( "ubar", un_adv(:,:) ) ! barotropic i-current
1266	CALL iom_put( "vbar", vn_adv(:,:) ) ! barotropic i-current
1267	!
1268	#if defined key_agrif
1269	! Save time integrated fluxes during child grid integration
1270	! (used to update coarse grid transports at next time step)
1271	!
1272	IF( .NOT.Agrif_Root() .AND. ln_bt_fw ) THEN
1273	IF( Agrif_NbStepint() == 0 ) THEN
1274	ub2_i_b(:,:) = 0._wp
1275	vb2_i_b(:,:) = 0._wp
1276	END IF
1277	!
1278	za1 = 1._wp / REAL(Agrif_rhot(), wp)
1279	ub2_i_b(:,:) = ub2_i_b(:,:) + za1 * ub2_b(:,:)
1280	vb2_i_b(:,:) = vb2_i_b(:,:) + za1 * vb2_b(:,:)
1281	ENDIF
1282	#endif
1283	! !* write time-spliting arrays in the restart
1284	IF( lrst_oce .AND.ln_bt_fw ) CALL ts_rst( kt, 'WRITE' )
1285	!
1286	CALL wrk_dealloc( jpi,jpj, zsshp2_e, zhdiv )
1287	CALL wrk_dealloc( jpi,jpj, zu_trd, zv_trd )
1288	CALL wrk_dealloc( jpi,jpj, zwx, zwy, zssh_frc, zu_frc, zv_frc )
1289	CALL wrk_dealloc( jpi,jpj, zhup2_e, zhvp2_e, zhust_e, zhvst_e )
1290	CALL wrk_dealloc( jpi,jpj, zsshu_a, zsshv_a )
1291	CALL wrk_dealloc( jpi,jpj, zhf )
1292	IF( ln_wd ) CALL wrk_dealloc( jpi, jpj, zcpx, zcpy )
1293	IF( ln_rwd ) CALL wrk_dealloc( jpi, jpj, ztwdmask, zuwdmask, zvwdmask, zuwdav2, zvwdav2 )
1294	!
1295	IF ( ln_diatmb ) THEN
1296	CALL iom_put( "baro_u" , un_bumask(:,:,1)+zmdi(1-umask(:,:,1 ) ) ) ! Barotropic U Velocity
1297	CALL iom_put( "baro_v" , vn_bvmask(:,:,1)+zmdi(1-vmask(:,:,1 ) ) ) ! Barotropic V Velocity
1298	ENDIF
1299	IF( nn_timing == 1 ) CALL timing_stop('dyn_spg_ts')
1300	!
1301	END SUBROUTINE dyn_spg_ts
1302
1303
1304	SUBROUTINE ts_wgt( ll_av, ll_fw, jpit, zwgt1, zwgt2)
1305	!!---------------------------------------------------------------------
1306	!! * ROUTINE ts_wgt *
1307	!!
1308	!! ** Purpose : Set time-splitting weights for temporal averaging (or not)
1309	!!----------------------------------------------------------------------
1310	LOGICAL, INTENT(in) :: ll_av ! temporal averaging=.true.
1311	LOGICAL, INTENT(in) :: ll_fw ! forward time splitting =.true.
1312	INTEGER, INTENT(inout) :: jpit ! cycle length
1313	REAL(wp), DIMENSION(3*nn_baro), INTENT(inout) :: zwgt1, & ! Primary weights
1314	zwgt2 ! Secondary weights
1315
1316	INTEGER :: jic, jn, ji ! temporary integers
1317	REAL(wp) :: za1, za2
1318	!!----------------------------------------------------------------------
1319
1320	zwgt1(:) = 0._wp
1321	zwgt2(:) = 0._wp
1322
1323	! Set time index when averaged value is requested
1324	IF (ll_fw) THEN
1325	jic = nn_baro
1326	ELSE
1327	jic = 2 * nn_baro
1328	ENDIF
1329
1330	! Set primary weights:
1331	IF (ll_av) THEN
1332	! Define simple boxcar window for primary weights
1333	! (width = nn_baro, centered around jic)
1334	SELECT CASE ( nn_bt_flt )
1335	CASE( 0 ) ! No averaging
1336	zwgt1(jic) = 1._wp
1337	jpit = jic
1338
1339	CASE( 1 ) ! Boxcar, width = nn_baro
1340	DO jn = 1, 3*nn_baro
1341	za1 = ABS(float(jn-jic))/float(nn_baro)
1342	IF (za1 < 0.5_wp) THEN
1343	zwgt1(jn) = 1._wp
1344	jpit = jn
1345	ENDIF
1346	ENDDO
1347
1348	CASE( 2 ) ! Boxcar, width = 2 * nn_baro
1349	DO jn = 1, 3*nn_baro
1350	za1 = ABS(float(jn-jic))/float(nn_baro)
1351	IF (za1 < 1._wp) THEN
1352	zwgt1(jn) = 1._wp
1353	jpit = jn
1354	ENDIF
1355	ENDDO
1356	CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_bt_flt' )
1357	END SELECT
1358
1359	ELSE ! No time averaging
1360	zwgt1(jic) = 1._wp
1361	jpit = jic
1362	ENDIF
1363
1364	! Set secondary weights
1365	DO jn = 1, jpit
1366	DO ji = jn, jpit
1367	zwgt2(jn) = zwgt2(jn) + zwgt1(ji)
1368	END DO
1369	END DO
1370
1371	! Normalize weigths:
1372	za1 = 1._wp / SUM(zwgt1(1:jpit))
1373	za2 = 1._wp / SUM(zwgt2(1:jpit))
1374	DO jn = 1, jpit
1375	zwgt1(jn) = zwgt1(jn) * za1
1376	zwgt2(jn) = zwgt2(jn) * za2
1377	END DO
1378	!
1379	END SUBROUTINE ts_wgt
1380
1381
1382	SUBROUTINE ts_rst( kt, cdrw )
1383	!!---------------------------------------------------------------------
1384	!! * ROUTINE ts_rst *
1385	!!
1386	!! ** Purpose : Read or write time-splitting arrays in restart file
1387	!!----------------------------------------------------------------------
1388	INTEGER , INTENT(in) :: kt ! ocean time-step
1389	CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag
1390	!
1391	!!----------------------------------------------------------------------
1392	!
1393	IF( TRIM(cdrw) == 'READ' ) THEN
1394	CALL iom_get( numror, jpdom_autoglo, 'ub2_b' , ub2_b (:,:) )
1395	CALL iom_get( numror, jpdom_autoglo, 'vb2_b' , vb2_b (:,:) )
1396	IF( .NOT.ln_bt_av ) THEN
1397	CALL iom_get( numror, jpdom_autoglo, 'sshbb_e' , sshbb_e(:,:) )
1398	CALL iom_get( numror, jpdom_autoglo, 'ubb_e' , ubb_e(:,:) )
1399	CALL iom_get( numror, jpdom_autoglo, 'vbb_e' , vbb_e(:,:) )
1400	CALL iom_get( numror, jpdom_autoglo, 'sshb_e' , sshb_e(:,:) )
1401	CALL iom_get( numror, jpdom_autoglo, 'ub_e' , ub_e(:,:) )
1402	CALL iom_get( numror, jpdom_autoglo, 'vb_e' , vb_e(:,:) )
1403	ENDIF
1404	#if defined key_agrif
1405	! Read time integrated fluxes
1406	IF ( .NOT.Agrif_Root() ) THEN
1407	CALL iom_get( numror, jpdom_autoglo, 'ub2_i_b' , ub2_i_b(:,:) )
1408	CALL iom_get( numror, jpdom_autoglo, 'vb2_i_b' , vb2_i_b(:,:) )
1409	ENDIF
1410	#endif
1411	!
1412	ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN
1413	CALL iom_rstput( kt, nitrst, numrow, 'ub2_b' , ub2_b (:,:) )
1414	CALL iom_rstput( kt, nitrst, numrow, 'vb2_b' , vb2_b (:,:) )
1415	!
1416	IF (.NOT.ln_bt_av) THEN
1417	CALL iom_rstput( kt, nitrst, numrow, 'sshbb_e' , sshbb_e(:,:) )
1418	CALL iom_rstput( kt, nitrst, numrow, 'ubb_e' , ubb_e(:,:) )
1419	CALL iom_rstput( kt, nitrst, numrow, 'vbb_e' , vbb_e(:,:) )
1420	CALL iom_rstput( kt, nitrst, numrow, 'sshb_e' , sshb_e(:,:) )
1421	CALL iom_rstput( kt, nitrst, numrow, 'ub_e' , ub_e(:,:) )
1422	CALL iom_rstput( kt, nitrst, numrow, 'vb_e' , vb_e(:,:) )
1423	ENDIF
1424	#if defined key_agrif
1425	! Save time integrated fluxes
1426	IF ( .NOT.Agrif_Root() ) THEN
1427	CALL iom_rstput( kt, nitrst, numrow, 'ub2_i_b' , ub2_i_b(:,:) )
1428	CALL iom_rstput( kt, nitrst, numrow, 'vb2_i_b' , vb2_i_b(:,:) )
1429	ENDIF
1430	#endif
1431	ENDIF
1432	!
1433	END SUBROUTINE ts_rst
1434
1435
1436	SUBROUTINE dyn_spg_ts_init
1437	!!---------------------------------------------------------------------
1438	!! * ROUTINE dyn_spg_ts_init *
1439	!!
1440	!! ** Purpose : Set time splitting options
1441	!!----------------------------------------------------------------------
1442	INTEGER :: ji ,jj ! dummy loop indices
1443	REAL(wp) :: zxr2, zyr2, zcmax ! local scalar
1444	REAL(wp), POINTER, DIMENSION(:,:) :: zcu
1445	!!----------------------------------------------------------------------
1446	!
1447	! Max courant number for ext. grav. waves
1448	!
1449	CALL wrk_alloc( jpi,jpj, zcu )
1450	!
1451	DO jj = 1, jpj
1452	DO ji =1, jpi
1453	zxr2 = r1_e1t(ji,jj) * r1_e1t(ji,jj)
1454	zyr2 = r1_e2t(ji,jj) * r1_e2t(ji,jj)
1455	zcu(ji,jj) = SQRT( grav * MAX(ht_0(ji,jj),0._wp) * (zxr2 + zyr2) )
1456	END DO
1457	END DO
1458	!
1459	zcmax = MAXVAL( zcu(:,:) )
1460	IF( lk_mpp ) CALL mpp_max( zcmax )
1461
1462	! Estimate number of iterations to satisfy a max courant number= rn_bt_cmax
1463	IF( ln_bt_auto ) nn_baro = CEILING( rdt / rn_bt_cmax * zcmax)
1464
1465	rdtbt = rdt / REAL( nn_baro , wp )
1466	zcmax = zcmax * rdtbt
1467	! Print results
1468	IF(lwp) WRITE(numout,*)
1469	IF(lwp) WRITE(numout,*) 'dyn_spg_ts : split-explicit free surface'
1470	IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
1471	IF( ln_bt_auto ) THEN
1472	IF(lwp) WRITE(numout,*) ' ln_ts_auto=.true. Automatically set nn_baro '
1473	IF(lwp) WRITE(numout,*) ' Max. courant number allowed: ', rn_bt_cmax
1474	ELSE
1475	IF(lwp) WRITE(numout,*) ' ln_ts_auto=.false.: Use nn_baro in namelist '
1476	ENDIF
1477
1478	IF(ln_bt_av) THEN
1479	IF(lwp) WRITE(numout,*) ' ln_bt_av=.true. => Time averaging over nn_baro time steps is on '
1480	ELSE
1481	IF(lwp) WRITE(numout,*) ' ln_bt_av=.false. => No time averaging of barotropic variables '
1482	ENDIF
1483	!
1484	!
1485	IF(ln_bt_fw) THEN
1486	IF(lwp) WRITE(numout,*) ' ln_bt_fw=.true. => Forward integration of barotropic variables '
1487	ELSE
1488	IF(lwp) WRITE(numout,*) ' ln_bt_fw =.false.=> Centred integration of barotropic variables '
1489	ENDIF
1490	!
1491	#if defined key_agrif
1492	! Restrict the use of Agrif to the forward case only
1493	IF( .NOT.ln_bt_fw .AND. .NOT.Agrif_Root() ) CALL ctl_stop( 'AGRIF not implemented if ln_bt_fw=.FALSE.' )
1494	#endif
1495	!
1496	IF(lwp) WRITE(numout,*) ' Time filter choice, nn_bt_flt: ', nn_bt_flt
1497	SELECT CASE ( nn_bt_flt )
1498	CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' Dirac'
1499	CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' Boxcar: width = nn_baro'
1500	CASE( 2 ) ; IF(lwp) WRITE(numout,) ' Boxcar: width = 2nn_baro'
1501	CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_bt_flt: should 0,1,2' )
1502	END SELECT
1503	!
1504	IF(lwp) WRITE(numout,*) ' '
1505	IF(lwp) WRITE(numout,*) ' nn_baro = ', nn_baro
1506	IF(lwp) WRITE(numout,*) ' Barotropic time step [s] is :', rdtbt
1507	IF(lwp) WRITE(numout,*) ' Maximum Courant number is :', zcmax
1508	!
1509	IF( .NOT.ln_bt_av .AND. .NOT.ln_bt_fw ) THEN
1510	CALL ctl_stop( 'dynspg_ts ERROR: No time averaging => only forward integration is possible' )
1511	ENDIF
1512	IF( zcmax>0.9_wp ) THEN
1513	CALL ctl_stop( 'dynspg_ts ERROR: Maximum Courant number is greater than 0.9: Inc. nn_baro !' )
1514	ENDIF
1515	!
1516	CALL wrk_dealloc( jpi,jpj, zcu )
1517	!
1518	END SUBROUTINE dyn_spg_ts_init
1519
1520	!!======================================================================
1521	END MODULE dynspg_ts

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

Download in other formats: