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

Last change on this file since 9043 was 9043, checked in by acc, 6 years ago
Branch 2017/dev_merge_2017. revert renamed constants in dynspg_ts.F90
Property svn:keywords set to `Id`
File size: 71.0 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) = 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) + rCdU_top(ji+1,jj)+rCdU_top(ji,jj) )
229	zCdU_v(ji,jj) = 0.5*( 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) = 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) )
236	zCdU_v(ji,jj) = 0.5*( 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	zu_frc(:,:) = zu_frc(:,:) + MAX(r1_hu_n(:,:) * zCdU_u(:,:),-1._wp / rdtbt) * zwx(:,:) * wdrampu(:,:)
526	zv_frc(:,:) = zv_frc(:,:) + MAX(r1_hv_n(:,:) * zCdU_v(:,:),-1._wp / rdtbt) * zwy(:,:) * wdrampv(:,:)
527	ELSE
528	zu_frc(:,:) = zu_frc(:,:) + r1_hu_n(:,:) * zCdU_u(:,:) * zwx(:,:)
529	zv_frc(:,:) = zv_frc(:,:) + r1_hv_n(:,:) * zCdU_v(:,:) * zwy(:,:)
530	END IF
531	!
532	! ! Add top stress contribution from baroclinic velocities:
533	IF( ln_bt_fw ) THEN
534	DO jj = 2, jpjm1
535	DO ji = fs_2, fs_jpim1 ! vector opt.
536	zu_frc(ji,jj) = zu_frc(ji,jj) + MAX( r1_hu_n(ji,jj) * 0.5( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) , zztmp ) zwx(ji,jj)
537	zv_frc(ji,jj) = zv_frc(ji,jj) + MAX( r1_hv_n(ji,jj) * 0.5( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) , zztmp ) zwy(ji,jj)
538	END DO
539	END DO
540	ELSE
541	DO jj = 2, jpjm1
542	DO ji = fs_2, fs_jpim1 ! vector opt.
543	zu_frc(ji,jj) = zu_frc(ji,jj) + r1_hu_n(ji,jj) * 0.5( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) zwx(ji,jj)
544	zv_frc(ji,jj) = zv_frc(ji,jj) + r1_hv_n(ji,jj) * 0.5( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) zwy(ji,jj)
545	END DO
546	END DO
547	ENDIF
548	!
549	IF( ln_isfcav ) THEN ! Add TOP stress contribution from baroclinic velocities:
550	IF( ln_bt_fw ) THEN
551	DO jj = 2, jpjm1
552	DO ji = fs_2, fs_jpim1 ! vector opt.
553	iktu = miku(ji,jj)
554	iktv = mikv(ji,jj)
555	zwx(ji,jj) = un(ji,jj,iktu) - un_b(ji,jj) ! NOW top baroclinic velocities
556	zwy(ji,jj) = vn(ji,jj,iktv) - vn_b(ji,jj)
557	END DO
558	END DO
559	ELSE
560	DO jj = 2, jpjm1
561	DO ji = fs_2, fs_jpim1 ! vector opt.
562	iktu = miku(ji,jj)
563	iktv = mikv(ji,jj)
564	zwx(ji,jj) = ub(ji,jj,iktu) - ub_b(ji,jj) ! BEFORE top baroclinic velocities
565	zwy(ji,jj) = vb(ji,jj,iktv) - vb_b(ji,jj)
566	END DO
567	END DO
568	ENDIF
569	!
570	! Note that the "unclipped" top friction parameter is used even with explicit drag
571	DO jj = 2, jpjm1
572	DO ji = fs_2, fs_jpim1 ! vector opt.
573	zu_frc(ji,jj) = zu_frc(ji,jj) + r1_hu_n(ji,jj) * 0.5( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) zwx(ji,jj)
574	zv_frc(ji,jj) = zv_frc(ji,jj) + r1_hv_n(ji,jj) * 0.5( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) zwy(ji,jj)
575	END DO
576	END DO
577	ENDIF
578	!
579	IF( ln_bt_fw ) THEN ! Add wind forcing
580	DO jj = 2, jpjm1
581	DO ji = fs_2, fs_jpim1 ! vector opt.
582	zu_frc(ji,jj) = zu_frc(ji,jj) + r1_rau0 * utau(ji,jj) * r1_hu_n(ji,jj)
583	zv_frc(ji,jj) = zv_frc(ji,jj) + r1_rau0 * vtau(ji,jj) * r1_hv_n(ji,jj)
584	END DO
585	END DO
586	ELSE
587	zztmp = r1_rau0 * r1_2
588	DO jj = 2, jpjm1
589	DO ji = fs_2, fs_jpim1 ! vector opt.
590	zu_frc(ji,jj) = zu_frc(ji,jj) + zztmp * ( utau_b(ji,jj) + utau(ji,jj) ) * r1_hu_n(ji,jj)
591	zv_frc(ji,jj) = zv_frc(ji,jj) + zztmp * ( vtau_b(ji,jj) + vtau(ji,jj) ) * r1_hv_n(ji,jj)
592	END DO
593	END DO
594	ENDIF
595	!
596	IF( ln_apr_dyn ) THEN ! Add atm pressure forcing
597	IF( ln_bt_fw ) THEN
598	DO jj = 2, jpjm1
599	DO ji = fs_2, fs_jpim1 ! vector opt.
600	zu_spg = grav * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) ) * r1_e1u(ji,jj)
601	zv_spg = grav * ( ssh_ib (ji ,jj+1) - ssh_ib (ji,jj) ) * r1_e2v(ji,jj)
602	zu_frc(ji,jj) = zu_frc(ji,jj) + zu_spg
603	zv_frc(ji,jj) = zv_frc(ji,jj) + zv_spg
604	END DO
605	END DO
606	ELSE
607	zztmp = grav * r1_2
608	DO jj = 2, jpjm1
609	DO ji = fs_2, fs_jpim1 ! vector opt.
610	zu_spg = zztmp * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) &
611	& + ssh_ibb(ji+1,jj ) - ssh_ibb(ji,jj) ) * r1_e1u(ji,jj)
612	zv_spg = zztmp * ( ssh_ib (ji ,jj+1) - ssh_ib (ji,jj) &
613	& + ssh_ibb(ji ,jj+1) - ssh_ibb(ji,jj) ) * r1_e2v(ji,jj)
614	zu_frc(ji,jj) = zu_frc(ji,jj) + zu_spg
615	zv_frc(ji,jj) = zv_frc(ji,jj) + zv_spg
616	END DO
617	END DO
618	ENDIF
619	ENDIF
620	! !* Right-Hand-Side of the barotropic ssh equation
621	! ! -----------------------------------------------
622	! ! Surface net water flux and rivers
623	IF (ln_bt_fw) THEN
624	zssh_frc(:,:) = r1_rau0 * ( emp(:,:) - rnf(:,:) + fwfisf(:,:) )
625	ELSE
626	zztmp = r1_rau0 * r1_2
627	zssh_frc(:,:) = zztmp * ( emp(:,:) + emp_b(:,:) - rnf(:,:) - rnf_b(:,:) &
628	& + fwfisf(:,:) + fwfisf_b(:,:) )
629	ENDIF
630	!
631	IF( ln_sdw ) THEN ! Stokes drift divergence added if necessary
632	zssh_frc(:,:) = zssh_frc(:,:) + div_sd(:,:)
633	ENDIF
634	!
635	#if defined key_asminc
636	! ! Include the IAU weighted SSH increment
637	IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN
638	zssh_frc(:,:) = zssh_frc(:,:) - ssh_iau(:,:)
639	ENDIF
640	#endif
641	! !* Fill boundary data arrays for AGRIF
642	! ! ------------------------------------
643	#if defined key_agrif
644	IF( .NOT.Agrif_Root() ) CALL agrif_dta_ts( kt )
645	#endif
646	!
647	! -----------------------------------------------------------------------
648	! Phase 2 : Integration of the barotropic equations
649	! -----------------------------------------------------------------------
650	!
651	! ! ==================== !
652	! ! Initialisations !
653	! ! ==================== !
654	! Initialize barotropic variables:
655	IF( ll_init )THEN
656	sshbb_e(:,:) = 0._wp
657	ubb_e (:,:) = 0._wp
658	vbb_e (:,:) = 0._wp
659	sshb_e (:,:) = 0._wp
660	ub_e (:,:) = 0._wp
661	vb_e (:,:) = 0._wp
662	ENDIF
663
664	!
665	IF (ln_bt_fw) THEN ! FORWARD integration: start from NOW fields
666	sshn_e(:,:) = sshn(:,:)
667	un_e (:,:) = un_b(:,:)
668	vn_e (:,:) = vn_b(:,:)
669	!
670	hu_e (:,:) = hu_n(:,:)
671	hv_e (:,:) = hv_n(:,:)
672	hur_e (:,:) = r1_hu_n(:,:)
673	hvr_e (:,:) = r1_hv_n(:,:)
674	ELSE ! CENTRED integration: start from BEFORE fields
675	sshn_e(:,:) = sshb(:,:)
676	un_e (:,:) = ub_b(:,:)
677	vn_e (:,:) = vb_b(:,:)
678	!
679	hu_e (:,:) = hu_b(:,:)
680	hv_e (:,:) = hv_b(:,:)
681	hur_e (:,:) = r1_hu_b(:,:)
682	hvr_e (:,:) = r1_hv_b(:,:)
683	ENDIF
684	!
685	!
686	!
687	! Initialize sums:
688	ua_b (:,:) = 0._wp ! After barotropic velocities (or transport if flux form)
689	va_b (:,:) = 0._wp
690	ssha (:,:) = 0._wp ! Sum for after averaged sea level
691	un_adv(:,:) = 0._wp ! Sum for now transport issued from ts loop
692	vn_adv(:,:) = 0._wp
693	! ! ==================== !
694
695	IF (ln_wd_dl) THEN
696	zuwdmask(:,:) = 0._wp ! set to zero for definiteness (not sure this is necessary)
697	zvwdmask(:,:) = 0._wp !
698	zuwdav2(:,:) = 0._wp
699	zvwdav2(:,:) = 0._wp
700	END IF
701
702
703	DO jn = 1, icycle ! sub-time-step loop !
704	! ! ==================== !
705	! !* Update the forcing (BDY and tides)
706	! ! ------------------
707	! Update only tidal forcing at open boundaries
708	IF( ln_bdy .AND. ln_tide ) CALL bdy_dta_tides( kt, kit=jn, time_offset= noffset+1 )
709	IF( ln_tide_pot .AND. ln_tide ) CALL upd_tide ( kt, kit=jn, time_offset= noffset )
710	!
711	! Set extrapolation coefficients for predictor step:
712	IF ((jn<3).AND.ll_init) THEN ! Forward
713	za1 = 1._wp
714	za2 = 0._wp
715	za3 = 0._wp
716	ELSE ! AB3-AM4 Coefficients: bet=0.281105
717	za1 = 1.781105_wp ! za1 = 3/2 + bet
718	za2 = -1.06221_wp ! za2 = -(1/2 + 2*bet)
719	za3 = 0.281105_wp ! za3 = bet
720	ENDIF
721
722	! Extrapolate barotropic velocities at step jit+0.5:
723	ua_e(:,:) = za1 * un_e(:,:) + za2 * ub_e(:,:) + za3 * ubb_e(:,:)
724	va_e(:,:) = za1 * vn_e(:,:) + za2 * vb_e(:,:) + za3 * vbb_e(:,:)
725
726	IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only)
727	! ! ------------------
728	! Extrapolate Sea Level at step jit+0.5:
729	zsshp2_e(:,:) = za1 * sshn_e(:,:) + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:)
730
731	! set wetting & drying mask at tracer points for this barotropic sub-step
732	IF ( ln_wd_dl ) THEN
733
734	IF ( ln_wd_dl_rmp ) THEN
735	DO jj = 1, jpj
736	DO ji = 1, jpi ! vector opt.
737	IF ( zsshp2_e(ji,jj) + ht_0(ji,jj) > 2._wp * rn_wdmin1 ) THEN
738	! IF ( zsshp2_e(ji,jj) + ht_0(ji,jj) > rn_wdmin2 ) THEN
739	ztwdmask(ji,jj) = 1._wp
740	ELSE IF ( zsshp2_e(ji,jj) + ht_0(ji,jj) > rn_wdmin1 ) THEN
741	ztwdmask(ji,jj) = (tanh(50._wp( ( zsshp2_e(ji,jj) + ht_0(ji,jj) - rn_wdmin1 )r_rn_wdmin1)) )
742	ELSE
743	ztwdmask(ji,jj) = 0._wp
744	END IF
745	END DO
746	END DO
747	ELSE
748	DO jj = 1, jpj
749	DO ji = 1, jpi ! vector opt.
750	IF ( zsshp2_e(ji,jj) + ht_0(ji,jj) > rn_wdmin1 ) THEN
751	ztwdmask(ji,jj) = 1._wp
752	ELSE
753	ztwdmask(ji,jj) = 0._wp
754	END IF
755	END DO
756	END DO
757	END IF
758
759	END IF
760
761
762	DO jj = 2, jpjm1 ! Sea Surface Height at u- & v-points
763	DO ji = 2, fs_jpim1 ! Vector opt.
764	zwx(ji,jj) = r1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj) &
765	& * ( e1e2t(ji ,jj) * zsshp2_e(ji ,jj) &
766	& + e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj) )
767	zwy(ji,jj) = r1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj) &
768	& * ( e1e2t(ji,jj ) * zsshp2_e(ji,jj ) &
769	& + e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1) )
770	END DO
771	END DO
772	CALL lbc_lnk_multi( zwx, 'U', 1._wp, zwy, 'V', 1._wp )
773	!
774	zhup2_e (:,:) = hu_0(:,:) + zwx(:,:) ! Ocean depth at U- and V-points
775	zhvp2_e (:,:) = hv_0(:,:) + zwy(:,:)
776	ELSE
777	zhup2_e (:,:) = hu_n(:,:)
778	zhvp2_e (:,:) = hv_n(:,:)
779	ENDIF
780	! !* after ssh
781	! ! -----------
782	! One should enforce volume conservation at open boundaries here
783	! considering fluxes below:
784	!
785	zwx(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:) ! fluxes at jn+0.5
786	zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:)
787	!
788	#if defined key_agrif
789	! Set fluxes during predictor step to ensure volume conservation
790	IF( .NOT.Agrif_Root() .AND. ln_bt_fw ) THEN
791	IF((nbondi == -1).OR.(nbondi == 2)) THEN
792	DO jj = 1, jpj
793	zwx(2:nbghostcells+1,jj) = ubdy_w(jj) * e2u(2:nbghostcells+1,jj)
794	END DO
795	ENDIF
796	IF((nbondi == 1).OR.(nbondi == 2)) THEN
797	DO jj=1,jpj
798	zwx(nlci-nbghostcells-1:nlci-2,jj) = ubdy_e(jj) * e2u(nlci-nbghostcells-1:nlci-2,jj)
799	END DO
800	ENDIF
801	IF((nbondj == -1).OR.(nbondj == 2)) THEN
802	DO ji=1,jpi
803	zwy(ji,2:nbghostcells+1) = vbdy_s(ji) * e1v(ji,2:nbghostcells+1)
804	END DO
805	ENDIF
806	IF((nbondj == 1).OR.(nbondj == 2)) THEN
807	DO ji=1,jpi
808	zwy(ji,nlcj-nbghostcells-1:nlcj-2) = vbdy_n(ji) * e1v(ji,nlcj-nbghostcells-1:nlcj-2)
809	END DO
810	ENDIF
811	ENDIF
812	#endif
813	IF( ln_wd_il ) CALL wad_lmt_bt(zwx, zwy, sshn_e, zssh_frc, rdtbt)
814
815	IF ( ln_wd_dl ) THEN
816
817
818	! un_e and vn_e are set to zero at faces where the direction of the flow is from dry cells
819
820	DO jj = 1, jpjm1
821	DO ji = 1, jpim1
822	IF ( zwx(ji,jj) > 0.0 ) THEN
823	zuwdmask(ji, jj) = ztwdmask(ji ,jj)
824	ELSE
825	zuwdmask(ji, jj) = ztwdmask(ji+1,jj)
826	END IF
827	zwx(ji, jj) = zuwdmask(ji,jj)*zwx(ji, jj)
828	un_e(ji,jj) = zuwdmask(ji,jj)*un_e(ji,jj)
829
830	IF ( zwy(ji,jj) > 0.0 ) THEN
831	zvwdmask(ji, jj) = ztwdmask(ji, jj )
832	ELSE
833	zvwdmask(ji, jj) = ztwdmask(ji, jj+1)
834	END IF
835	zwy(ji, jj) = zvwdmask(ji,jj)*zwy(ji,jj)
836	vn_e(ji,jj) = zvwdmask(ji,jj)*vn_e(ji,jj)
837	END DO
838	END DO
839
840
841	END IF
842
843	! Sum over sub-time-steps to compute advective velocities
844	za2 = wgtbtp2(jn)
845	un_adv(:,:) = un_adv(:,:) + za2 * zwx(:,:) * r1_e2u(:,:)
846	vn_adv(:,:) = vn_adv(:,:) + za2 * zwy(:,:) * r1_e1v(:,:)
847
848	! sum over sub-time-steps to decide which baroclinic velocities to set to zero (zuwdav2 is only used when ln_wd_dl_bc = True)
849	IF ( ln_wd_dl_bc ) THEN
850	zuwdav2(:,:) = zuwdav2(:,:) + za2 * zuwdmask(:,:)
851	zvwdav2(:,:) = zvwdav2(:,:) + za2 * zvwdmask(:,:)
852	END IF
853
854	! Set next sea level:
855	DO jj = 2, jpjm1
856	DO ji = fs_2, fs_jpim1 ! vector opt.
857	zhdiv(ji,jj) = ( zwx(ji,jj) - zwx(ji-1,jj) &
858	& + zwy(ji,jj) - zwy(ji,jj-1) ) * r1_e1e2t(ji,jj)
859	END DO
860	END DO
861	ssha_e(:,:) = ( sshn_e(:,:) - rdtbt * ( zssh_frc(:,:) + zhdiv(:,:) ) ) * ssmask(:,:)
862
863	CALL lbc_lnk( ssha_e, 'T', 1._wp )
864
865	! Duplicate sea level across open boundaries (this is only cosmetic if linssh=T)
866	IF( ln_bdy ) CALL bdy_ssh( ssha_e )
867	#if defined key_agrif
868	IF( .NOT.Agrif_Root() ) CALL agrif_ssh_ts( jn )
869	#endif
870	!
871	! Sea Surface Height at u-,v-points (vvl case only)
872	IF( .NOT.ln_linssh ) THEN
873	DO jj = 2, jpjm1
874	DO ji = 2, jpim1 ! NO Vector Opt.
875	zsshu_a(ji,jj) = r1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj) &
876	& * ( e1e2t(ji ,jj ) * ssha_e(ji ,jj ) &
877	& + e1e2t(ji+1,jj ) * ssha_e(ji+1,jj ) )
878	zsshv_a(ji,jj) = r1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj) &
879	& * ( e1e2t(ji ,jj ) * ssha_e(ji ,jj ) &
880	& + e1e2t(ji ,jj+1) * ssha_e(ji ,jj+1) )
881	END DO
882	END DO
883	CALL lbc_lnk_multi( zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp )
884	ENDIF
885	!
886	! Half-step back interpolation of SSH for surface pressure computation:
887	!----------------------------------------------------------------------
888	IF ((jn==1).AND.ll_init) THEN
889	za0=1._wp ! Forward-backward
890	za1=0._wp
891	za2=0._wp
892	za3=0._wp
893	ELSEIF ((jn==2).AND.ll_init) THEN ! AB2-AM3 Coefficients; bet=0 ; gam=-1/6 ; eps=1/12
894	za0= 1.0833333333333_wp ! za0 = 1-gam-eps
895	za1=-0.1666666666666_wp ! za1 = gam
896	za2= 0.0833333333333_wp ! za2 = eps
897	za3= 0._wp
898	ELSE ! AB3-AM4 Coefficients; bet=0.281105 ; eps=0.013 ; gam=0.0880
899	IF (rn_bt_alpha==0._wp) THEN
900	za0=0.614_wp ! za0 = 1/2 + gam + 2*eps
901	za1=0.285_wp ! za1 = 1/2 - 2gam - 3eps
902	za2=0.088_wp ! za2 = gam
903	za3=0.013_wp ! za3 = eps
904	ELSE
905	zepsilon = 0.00976186_wp - 0.13451357_wp * rn_bt_alpha
906	zgamma = 0.08344500_wp - 0.51358400_wp * rn_bt_alpha
907	za0 = 0.5_wp + zgamma + 2._wp * rn_bt_alpha + 2._wp * zepsilon
908	za1 = 1._wp - za0 - zgamma - zepsilon
909	za2 = zgamma
910	za3 = zepsilon
911	ENDIF
912	ENDIF
913	!
914	zsshp2_e(:,:) = za0 * ssha_e(:,:) + za1 * sshn_e (:,:) &
915	& + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:)
916	IF( ln_wd_il ) THEN ! Calculating and applying W/D gravity filters
917	DO jj = 2, jpjm1
918	DO ji = 2, jpim1
919	ll_tmp1 = MIN( zsshp2_e(ji,jj) , zsshp2_e(ji+1,jj) ) > &
920	& MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) .AND. &
921	& MAX( zsshp2_e(ji,jj) + ht_0(ji,jj) , zsshp2_e(ji+1,jj) + ht_0(ji+1,jj) ) &
922	& > rn_wdmin1 + rn_wdmin2
923	ll_tmp2 = (ABS(zsshp2_e(ji,jj) - zsshp2_e(ji+1,jj)) > 1.E-12 ).AND.( &
924	& MAX( zsshp2_e(ji,jj) , zsshp2_e(ji+1,jj) ) > &
925	& MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 )
926
927	IF(ll_tmp1) THEN
928	zcpx(ji,jj) = 1.0_wp
929	ELSE IF(ll_tmp2) THEN
930	! no worries about zsshp2_e(ji+1,jj) - zsshp2_e(ji ,jj) = 0, it won't happen ! here
931	zcpx(ji,jj) = ABS( (zsshp2_e(ji+1,jj) + ht_0(ji+1,jj) - zsshp2_e(ji,jj) - ht_0(ji,jj)) &
932	& / (zsshp2_e(ji+1,jj) - zsshp2_e(ji ,jj)) )
933	ELSE
934	zcpx(ji,jj) = 0._wp
935	END IF
936
937	ll_tmp1 = MIN( zsshp2_e(ji,jj) , zsshp2_e(ji,jj+1) ) > &
938	& MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) .AND. &
939	& MAX( zsshp2_e(ji,jj) + ht_0(ji,jj) , zsshp2_e(ji,jj+1) + ht_0(ji,jj+1) ) &
940	& > rn_wdmin1 + rn_wdmin2
941	ll_tmp2 = (ABS(zsshp2_e(ji,jj) - zsshp2_e(ji,jj+1)) > 1.E-12 ).AND.( &
942	& MAX( zsshp2_e(ji,jj) , zsshp2_e(ji,jj+1) ) > &
943	& MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 )
944
945	IF(ll_tmp1) THEN
946	zcpy(ji,jj) = 1.0_wp
947	ELSE IF(ll_tmp2) THEN
948	! no worries about zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj ) = 0, it won't happen ! here
949	zcpy(ji,jj) = ABS( (zsshp2_e(ji,jj+1) + ht_0(ji,jj+1) - zsshp2_e(ji,jj) - ht_0(ji,jj)) &
950	& / (zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj )) )
951	ELSE
952	zcpy(ji,jj) = 0._wp
953	END IF
954	END DO
955	END DO
956	END IF
957	!
958	! Compute associated depths at U and V points:
959	IF( .NOT.ln_linssh .AND. .NOT.ln_dynadv_vec ) THEN !* Vector form
960	!
961	DO jj = 2, jpjm1
962	DO ji = 2, jpim1
963	zx1 = r1_2 * ssumask(ji ,jj) * r1_e1e2u(ji ,jj) &
964	& * ( e1e2t(ji ,jj ) * zsshp2_e(ji ,jj) &
965	& + e1e2t(ji+1,jj ) * zsshp2_e(ji+1,jj ) )
966	zy1 = r1_2 * ssvmask(ji ,jj) * r1_e1e2v(ji ,jj ) &
967	& * ( e1e2t(ji ,jj ) * zsshp2_e(ji ,jj ) &
968	& + e1e2t(ji ,jj+1) * zsshp2_e(ji ,jj+1) )
969	zhust_e(ji,jj) = hu_0(ji,jj) + zx1
970	zhvst_e(ji,jj) = hv_0(ji,jj) + zy1
971	END DO
972	END DO
973
974	ENDIF
975	!
976	! Add Coriolis trend:
977	! zwz array below or triads normally depend on sea level with ln_linssh=F and should be updated
978	! at each time step. We however keep them constant here for optimization.
979	! Recall that zwx and zwy arrays hold fluxes at this stage:
980	! zwx(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:) ! fluxes at jn+0.5
981	! zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:)
982	!
983	IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN !== energy conserving or mixed scheme ==!
984	DO jj = 2, jpjm1
985	DO ji = fs_2, fs_jpim1 ! vector opt.
986	zy1 = ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) ) * r1_e1u(ji,jj)
987	zy2 = ( zwy(ji ,jj ) + zwy(ji+1,jj ) ) * r1_e1u(ji,jj)
988	zx1 = ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) ) * r1_e2v(ji,jj)
989	zx2 = ( zwx(ji ,jj ) + zwx(ji ,jj+1) ) * r1_e2v(ji,jj)
990	zu_trd(ji,jj) = r1_4 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
991	zv_trd(ji,jj) =-r1_4 * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 )
992	END DO
993	END DO
994	!
995	ELSEIF ( ln_dynvor_ens ) THEN !== enstrophy conserving scheme ==!
996	DO jj = 2, jpjm1
997	DO ji = fs_2, fs_jpim1 ! vector opt.
998	zy1 = r1_8 * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) &
999	& + zwy(ji ,jj ) + zwy(ji+1,jj ) ) * r1_e1u(ji,jj)
1000	zx1 = - r1_8 * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) &
1001	& + zwx(ji ,jj ) + zwx(ji ,jj+1) ) * r1_e2v(ji,jj)
1002	zu_trd(ji,jj) = zy1 * ( zwz(ji ,jj-1) + zwz(ji,jj) )
1003	zv_trd(ji,jj) = zx1 * ( zwz(ji-1,jj ) + zwz(ji,jj) )
1004	END DO
1005	END DO
1006	!
1007	ELSEIF ( ln_dynvor_een ) THEN !== energy and enstrophy conserving scheme ==!
1008	DO jj = 2, jpjm1
1009	DO ji = fs_2, fs_jpim1 ! vector opt.
1010	zu_trd(ji,jj) = + r1_12 * r1_e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) &
1011	& + ftnw(ji+1,jj) * zwy(ji+1,jj ) &
1012	& + ftse(ji,jj ) * zwy(ji ,jj-1) &
1013	& + ftsw(ji+1,jj) * zwy(ji+1,jj-1) )
1014	zv_trd(ji,jj) = - r1_12 * r1_e2v(ji,jj) * ( ftsw(ji,jj+1) * zwx(ji-1,jj+1) &
1015	& + ftse(ji,jj+1) * zwx(ji ,jj+1) &
1016	& + ftnw(ji,jj ) * zwx(ji-1,jj ) &
1017	& + ftne(ji,jj ) * zwx(ji ,jj ) )
1018	END DO
1019	END DO
1020	!
1021	ENDIF
1022	!
1023	! Add tidal astronomical forcing if defined
1024	IF ( ln_tide .AND. ln_tide_pot ) THEN
1025	DO jj = 2, jpjm1
1026	DO ji = fs_2, fs_jpim1 ! vector opt.
1027	zu_spg = grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) * r1_e1u(ji,jj)
1028	zv_spg = grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) * r1_e2v(ji,jj)
1029	zu_trd(ji,jj) = zu_trd(ji,jj) + zu_spg
1030	zv_trd(ji,jj) = zv_trd(ji,jj) + zv_spg
1031	END DO
1032	END DO
1033	ENDIF
1034	!
1035	! Add bottom stresses:
1036	!jth do implicitly instead
1037	IF ( .NOT. ll_wd ) THEN ! Revert to explicit for bit comparison tests in non wad runs
1038	zu_trd(:,:) = zu_trd(:,:) + zCdU_u(:,:) * un_e(:,:) * hur_e(:,:)
1039	zv_trd(:,:) = zv_trd(:,:) + zCdU_v(:,:) * vn_e(:,:) * hvr_e(:,:)
1040	ENDIF
1041	!
1042	! Surface pressure trend:
1043	IF( ln_wd_il ) THEN
1044	DO jj = 2, jpjm1
1045	DO ji = 2, jpim1
1046	! Add surface pressure gradient
1047	zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj)
1048	zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj)
1049	zwx(ji,jj) = (1._wp - rn_scal_load) * zu_spg * zcpx(ji,jj)
1050	zwy(ji,jj) = (1._wp - rn_scal_load) * zv_spg * zcpy(ji,jj)
1051	END DO
1052	END DO
1053	ELSE
1054	DO jj = 2, jpjm1
1055	DO ji = fs_2, fs_jpim1 ! vector opt.
1056	! Add surface pressure gradient
1057	zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj)
1058	zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj)
1059	zwx(ji,jj) = (1._wp - rn_scal_load) * zu_spg
1060	zwy(ji,jj) = (1._wp - rn_scal_load) * zv_spg
1061	END DO
1062	END DO
1063	END IF
1064
1065	!
1066	! Set next velocities:
1067	IF( ln_dynadv_vec .OR. ln_linssh ) THEN !* Vector form
1068	DO jj = 2, jpjm1
1069	DO ji = fs_2, fs_jpim1 ! vector opt.
1070	ua_e(ji,jj) = ( un_e(ji,jj) &
1071	& + rdtbt * ( zwx(ji,jj) &
1072	& + zu_trd(ji,jj) &
1073	& + zu_frc(ji,jj) ) &
1074	& ) * ssumask(ji,jj)
1075
1076	va_e(ji,jj) = ( vn_e(ji,jj) &
1077	& + rdtbt * ( zwy(ji,jj) &
1078	& + zv_trd(ji,jj) &
1079	& + zv_frc(ji,jj) ) &
1080	& ) * ssvmask(ji,jj)
1081
1082	!jth implicit bottom friction:
1083	IF ( ll_wd ) THEN ! revert to explicit for bit comparison tests in non wad runs
1084	ua_e(ji,jj) = ua_e(ji,jj) /(1.0 - rdtbt * zCdU_u(ji,jj) * hur_e(ji,jj))
1085	va_e(ji,jj) = va_e(ji,jj) /(1.0 - rdtbt * zCdU_v(ji,jj) * hvr_e(ji,jj))
1086	ENDIF
1087
1088	END DO
1089	END DO
1090	!
1091	ELSE !* Flux form
1092	DO jj = 2, jpjm1
1093	DO ji = fs_2, fs_jpim1 ! vector opt.
1094
1095	zhura = hu_0(ji,jj) + zsshu_a(ji,jj)
1096	zhvra = hv_0(ji,jj) + zsshv_a(ji,jj)
1097
1098	zhura = ssumask(ji,jj)/(zhura + 1._wp - ssumask(ji,jj))
1099	zhvra = ssvmask(ji,jj)/(zhvra + 1._wp - ssvmask(ji,jj))
1100
1101	ua_e(ji,jj) = ( hu_e(ji,jj) * un_e(ji,jj) &
1102	& + rdtbt * ( zhust_e(ji,jj) * zwx(ji,jj) &
1103	& + zhup2_e(ji,jj) * zu_trd(ji,jj) &
1104	& + hu_n(ji,jj) * zu_frc(ji,jj) ) &
1105	& ) * zhura
1106
1107	va_e(ji,jj) = ( hv_e(ji,jj) * vn_e(ji,jj) &
1108	& + rdtbt * ( zhvst_e(ji,jj) * zwy(ji,jj) &
1109	& + zhvp2_e(ji,jj) * zv_trd(ji,jj) &
1110	& + hv_n(ji,jj) * zv_frc(ji,jj) ) &
1111	& ) * zhvra
1112	END DO
1113	END DO
1114	ENDIF
1115
1116
1117	IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only)
1118	hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:)
1119	hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:)
1120	hur_e(:,:) = ssumask(:,:) / ( hu_e(:,:) + 1._wp - ssumask(:,:) )
1121	hvr_e(:,:) = ssvmask(:,:) / ( hv_e(:,:) + 1._wp - ssvmask(:,:) )
1122	!
1123	ENDIF
1124	! !* domain lateral boundary
1125	CALL lbc_lnk_multi( ua_e, 'U', -1._wp, va_e , 'V', -1._wp )
1126	!
1127	! ! open boundaries
1128	IF( ln_bdy ) CALL bdy_dyn2d( jn, ua_e, va_e, un_e, vn_e, hur_e, hvr_e, ssha_e )
1129	#if defined key_agrif
1130	IF( .NOT.Agrif_Root() ) CALL agrif_dyn_ts( jn ) ! Agrif
1131	#endif
1132	! !* Swap
1133	! ! ----
1134	ubb_e (:,:) = ub_e (:,:)
1135	ub_e (:,:) = un_e (:,:)
1136	un_e (:,:) = ua_e (:,:)
1137	!
1138	vbb_e (:,:) = vb_e (:,:)
1139	vb_e (:,:) = vn_e (:,:)
1140	vn_e (:,:) = va_e (:,:)
1141	!
1142	sshbb_e(:,:) = sshb_e(:,:)
1143	sshb_e (:,:) = sshn_e(:,:)
1144	sshn_e (:,:) = ssha_e(:,:)
1145
1146	! !* Sum over whole bt loop
1147	! ! ----------------------
1148	za1 = wgtbtp1(jn)
1149	IF( ln_dynadv_vec .OR. ln_linssh ) THEN ! Sum velocities
1150	ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:)
1151	va_b (:,:) = va_b (:,:) + za1 * va_e (:,:)
1152	ELSE ! Sum transports
1153	IF ( .NOT.ln_wd_dl ) THEN
1154	ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:) * hu_e (:,:)
1155	va_b (:,:) = va_b (:,:) + za1 * va_e (:,:) * hv_e (:,:)
1156	ELSE
1157	ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:) * hu_e (:,:) * zuwdmask(:,:)
1158	va_b (:,:) = va_b (:,:) + za1 * va_e (:,:) * hv_e (:,:) * zvwdmask(:,:)
1159	END IF
1160	ENDIF
1161	! ! Sum sea level
1162	ssha(:,:) = ssha(:,:) + za1 * ssha_e(:,:)
1163
1164	! ! ==================== !
1165	END DO ! end loop !
1166	! ! ==================== !
1167	! -----------------------------------------------------------------------------
1168	! Phase 3. update the general trend with the barotropic trend
1169	! -----------------------------------------------------------------------------
1170	!
1171	! Set advection velocity correction:
1172	IF (ln_bt_fw) THEN
1173	zwx(:,:) = un_adv(:,:)
1174	zwy(:,:) = vn_adv(:,:)
1175	IF( .NOT.( kt == nit000 .AND. neuler==0 ) ) THEN
1176	un_adv(:,:) = r1_2 * ( ub2_b(:,:) + zwx(:,:) - atfp * un_bf(:,:) )
1177	vn_adv(:,:) = r1_2 * ( vb2_b(:,:) + zwy(:,:) - atfp * vn_bf(:,:) )
1178	!
1179	! Update corrective fluxes for next time step:
1180	un_bf(:,:) = atfp * un_bf(:,:) + (zwx(:,:) - ub2_b(:,:))
1181	vn_bf(:,:) = atfp * vn_bf(:,:) + (zwy(:,:) - vb2_b(:,:))
1182	ELSE
1183	un_bf(:,:) = 0._wp
1184	vn_bf(:,:) = 0._wp
1185	END IF
1186	! Save integrated transport for next computation
1187	ub2_b(:,:) = zwx(:,:)
1188	vb2_b(:,:) = zwy(:,:)
1189	ENDIF
1190
1191
1192	!
1193	! Update barotropic trend:
1194	IF( ln_dynadv_vec .OR. ln_linssh ) THEN
1195	DO jk=1,jpkm1
1196	ua(:,:,jk) = ua(:,:,jk) + ( ua_b(:,:) - ub_b(:,:) ) * r1_2dt_b
1197	va(:,:,jk) = va(:,:,jk) + ( va_b(:,:) - vb_b(:,:) ) * r1_2dt_b
1198	END DO
1199	ELSE
1200	! At this stage, ssha has been corrected: compute new depths at velocity points
1201	DO jj = 1, jpjm1
1202	DO ji = 1, jpim1 ! NO Vector Opt.
1203	zsshu_a(ji,jj) = r1_2 * umask(ji,jj,1) * r1_e1e2u(ji,jj) &
1204	& * ( e1e2t(ji ,jj) * ssha(ji ,jj) &
1205	& + e1e2t(ji+1,jj) * ssha(ji+1,jj) )
1206	zsshv_a(ji,jj) = r1_2 * vmask(ji,jj,1) * r1_e1e2v(ji,jj) &
1207	& * ( e1e2t(ji,jj ) * ssha(ji,jj ) &
1208	& + e1e2t(ji,jj+1) * ssha(ji,jj+1) )
1209	END DO
1210	END DO
1211	CALL lbc_lnk_multi( zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) ! Boundary conditions
1212	!
1213	DO jk=1,jpkm1
1214	ua(:,:,jk) = ua(:,:,jk) + r1_hu_n(:,:) * ( ua_b(:,:) - ub_b(:,:) * hu_b(:,:) ) * r1_2dt_b
1215	va(:,:,jk) = va(:,:,jk) + r1_hv_n(:,:) * ( va_b(:,:) - vb_b(:,:) * hv_b(:,:) ) * r1_2dt_b
1216	END DO
1217	! Save barotropic velocities not transport:
1218	ua_b(:,:) = ua_b(:,:) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - ssumask(:,:) )
1219	va_b(:,:) = va_b(:,:) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - ssvmask(:,:) )
1220	ENDIF
1221
1222
1223	! Correct velocities so that the barotropic velocity equals (un_adv, vn_adv) (in all cases)
1224	DO jk = 1, jpkm1
1225	un(:,:,jk) = ( un(:,:,jk) + un_adv(:,:)r1_hu_n(:,:) - un_b(:,:) ) umask(:,:,jk)
1226	vn(:,:,jk) = ( vn(:,:,jk) + vn_adv(:,:)r1_hv_n(:,:) - vn_b(:,:) ) vmask(:,:,jk)
1227	END DO
1228
1229	IF ( ln_wd_dl .and. ln_wd_dl_bc) THEN
1230	DO jk = 1, jpkm1
1231	un(:,:,jk) = ( un_adv(:,:) + zuwdav2(:,:)(un(:,:,jk) - un_adv(:,:)) ) umask(:,:,jk)
1232	vn(:,:,jk) = ( vn_adv(:,:) + zvwdav2(:,:)(vn(:,:,jk) - vn_adv(:,:)) ) 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: