Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

dynspg_ts.F90 @ 9499

Last change on this file since 9499 was 9499, checked in by davestorkey, 6 years ago
branches/UKMO/dev_merge_2017_CICE_interface : clear SVN keywords.
File size: 72.3 KB

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

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

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

Context Navigation

source: branches/UKMO/dev_merge_2017_CICE_interface/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 9499

Download in other formats: