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

Context Navigation

source: branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 9112

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

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

Download in other formats: