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

Last change on this file since 9109 was 9109, checked in by jchanut, 6 years ago
Correct advective velocities for wetting and drying
Property svn:keywords set to `Id`
File size: 71.6 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	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) + MAX(r1_hu_n(ji,jj) * zCdU_u(ji,jj), zztmp ) * zwx(ji,jj) * wdrampu(ji,jj)
529	zv_frc(ji,jj) = zv_frc(ji,jj) + MAX(r1_hv_n(ji,jj) * zCdU_v(ji,jj), zztmp ) * zwy(ji,jj) * wdrampv(ji,jj)
530	END DO
531	END DO
532	ELSE
533	DO jj = 2, jpjm1
534	DO ji = fs_2, fs_jpim1 ! vector opt.
535	zu_frc(ji,jj) = zu_frc(ji,jj) + r1_hu_n(ji,jj) * zCdU_u(ji,jj) * zwx(ji,jj)
536	zv_frc(ji,jj) = zv_frc(ji,jj) + r1_hv_n(ji,jj) * zCdU_v(ji,jj) * zwy(ji,jj)
537	END DO
538	END DO
539	END IF
540	!
541	! ! Add top stress contribution from baroclinic velocities:
542	IF( ln_bt_fw ) THEN
543	zztmp = -1._wp / rdtbt
544	DO jj = 2, jpjm1
545	DO ji = fs_2, fs_jpim1 ! vector opt.
546	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)
547	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)
548	END DO
549	END DO
550	ELSE
551	DO jj = 2, jpjm1
552	DO ji = fs_2, fs_jpim1 ! vector opt.
553	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)
554	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)
555	END DO
556	END DO
557	ENDIF
558	!
559	IF( ln_isfcav ) THEN ! Add TOP stress contribution from baroclinic velocities:
560	IF( ln_bt_fw ) THEN
561	DO jj = 2, jpjm1
562	DO ji = fs_2, fs_jpim1 ! vector opt.
563	iktu = miku(ji,jj)
564	iktv = mikv(ji,jj)
565	zwx(ji,jj) = un(ji,jj,iktu) - un_b(ji,jj) ! NOW top baroclinic velocities
566	zwy(ji,jj) = vn(ji,jj,iktv) - vn_b(ji,jj)
567	END DO
568	END DO
569	ELSE
570	DO jj = 2, jpjm1
571	DO ji = fs_2, fs_jpim1 ! vector opt.
572	iktu = miku(ji,jj)
573	iktv = mikv(ji,jj)
574	zwx(ji,jj) = ub(ji,jj,iktu) - ub_b(ji,jj) ! BEFORE top baroclinic velocities
575	zwy(ji,jj) = vb(ji,jj,iktv) - vb_b(ji,jj)
576	END DO
577	END DO
578	ENDIF
579	!
580	! Note that the "unclipped" top friction parameter is used even with explicit drag
581	DO jj = 2, jpjm1
582	DO ji = fs_2, fs_jpim1 ! vector opt.
583	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)
584	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)
585	END DO
586	END DO
587	ENDIF
588	!
589	IF( ln_bt_fw ) THEN ! Add wind forcing
590	DO jj = 2, jpjm1
591	DO ji = fs_2, fs_jpim1 ! vector opt.
592	zu_frc(ji,jj) = zu_frc(ji,jj) + r1_rau0 * utau(ji,jj) * r1_hu_n(ji,jj)
593	zv_frc(ji,jj) = zv_frc(ji,jj) + r1_rau0 * vtau(ji,jj) * r1_hv_n(ji,jj)
594	END DO
595	END DO
596	ELSE
597	zztmp = r1_rau0 * r1_2
598	DO jj = 2, jpjm1
599	DO ji = fs_2, fs_jpim1 ! vector opt.
600	zu_frc(ji,jj) = zu_frc(ji,jj) + zztmp * ( utau_b(ji,jj) + utau(ji,jj) ) * r1_hu_n(ji,jj)
601	zv_frc(ji,jj) = zv_frc(ji,jj) + zztmp * ( vtau_b(ji,jj) + vtau(ji,jj) ) * r1_hv_n(ji,jj)
602	END DO
603	END DO
604	ENDIF
605	!
606	IF( ln_apr_dyn ) THEN ! Add atm pressure forcing
607	IF( ln_bt_fw ) THEN
608	DO jj = 2, jpjm1
609	DO ji = fs_2, fs_jpim1 ! vector opt.
610	zu_spg = grav * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) ) * r1_e1u(ji,jj)
611	zv_spg = grav * ( ssh_ib (ji ,jj+1) - ssh_ib (ji,jj) ) * r1_e2v(ji,jj)
612	zu_frc(ji,jj) = zu_frc(ji,jj) + zu_spg
613	zv_frc(ji,jj) = zv_frc(ji,jj) + zv_spg
614	END DO
615	END DO
616	ELSE
617	zztmp = grav * r1_2
618	DO jj = 2, jpjm1
619	DO ji = fs_2, fs_jpim1 ! vector opt.
620	zu_spg = zztmp * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) &
621	& + ssh_ibb(ji+1,jj ) - ssh_ibb(ji,jj) ) * r1_e1u(ji,jj)
622	zv_spg = zztmp * ( ssh_ib (ji ,jj+1) - ssh_ib (ji,jj) &
623	& + ssh_ibb(ji ,jj+1) - ssh_ibb(ji,jj) ) * r1_e2v(ji,jj)
624	zu_frc(ji,jj) = zu_frc(ji,jj) + zu_spg
625	zv_frc(ji,jj) = zv_frc(ji,jj) + zv_spg
626	END DO
627	END DO
628	ENDIF
629	ENDIF
630	! !* Right-Hand-Side of the barotropic ssh equation
631	! ! -----------------------------------------------
632	! ! Surface net water flux and rivers
633	IF (ln_bt_fw) THEN
634	zssh_frc(:,:) = r1_rau0 * ( emp(:,:) - rnf(:,:) + fwfisf(:,:) )
635	ELSE
636	zztmp = r1_rau0 * r1_2
637	zssh_frc(:,:) = zztmp * ( emp(:,:) + emp_b(:,:) - rnf(:,:) - rnf_b(:,:) &
638	& + fwfisf(:,:) + fwfisf_b(:,:) )
639	ENDIF
640	!
641	IF( ln_sdw ) THEN ! Stokes drift divergence added if necessary
642	zssh_frc(:,:) = zssh_frc(:,:) + div_sd(:,:)
643	ENDIF
644	!
645	#if defined key_asminc
646	! ! Include the IAU weighted SSH increment
647	IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN
648	zssh_frc(:,:) = zssh_frc(:,:) - ssh_iau(:,:)
649	ENDIF
650	#endif
651	! !* Fill boundary data arrays for AGRIF
652	! ! ------------------------------------
653	#if defined key_agrif
654	IF( .NOT.Agrif_Root() ) CALL agrif_dta_ts( kt )
655	#endif
656	!
657	! -----------------------------------------------------------------------
658	! Phase 2 : Integration of the barotropic equations
659	! -----------------------------------------------------------------------
660	!
661	! ! ==================== !
662	! ! Initialisations !
663	! ! ==================== !
664	! Initialize barotropic variables:
665	IF( ll_init )THEN
666	sshbb_e(:,:) = 0._wp
667	ubb_e (:,:) = 0._wp
668	vbb_e (:,:) = 0._wp
669	sshb_e (:,:) = 0._wp
670	ub_e (:,:) = 0._wp
671	vb_e (:,:) = 0._wp
672	ENDIF
673
674	!
675	IF (ln_bt_fw) THEN ! FORWARD integration: start from NOW fields
676	sshn_e(:,:) = sshn(:,:)
677	un_e (:,:) = un_b(:,:)
678	vn_e (:,:) = vn_b(:,:)
679	!
680	hu_e (:,:) = hu_n(:,:)
681	hv_e (:,:) = hv_n(:,:)
682	hur_e (:,:) = r1_hu_n(:,:)
683	hvr_e (:,:) = r1_hv_n(:,:)
684	ELSE ! CENTRED integration: start from BEFORE fields
685	sshn_e(:,:) = sshb(:,:)
686	un_e (:,:) = ub_b(:,:)
687	vn_e (:,:) = vb_b(:,:)
688	!
689	hu_e (:,:) = hu_b(:,:)
690	hv_e (:,:) = hv_b(:,:)
691	hur_e (:,:) = r1_hu_b(:,:)
692	hvr_e (:,:) = r1_hv_b(:,:)
693	ENDIF
694	!
695	!
696	!
697	! Initialize sums:
698	ua_b (:,:) = 0._wp ! After barotropic velocities (or transport if flux form)
699	va_b (:,:) = 0._wp
700	ssha (:,:) = 0._wp ! Sum for after averaged sea level
701	un_adv(:,:) = 0._wp ! Sum for now transport issued from ts loop
702	vn_adv(:,:) = 0._wp
703	! ! ==================== !
704
705	IF (ln_wd_dl) THEN
706	zuwdmask(:,:) = 0._wp ! set to zero for definiteness (not sure this is necessary)
707	zvwdmask(:,:) = 0._wp !
708	zuwdav2(:,:) = 0._wp
709	zvwdav2(:,:) = 0._wp
710	END IF
711
712
713	DO jn = 1, icycle ! sub-time-step loop !
714	! ! ==================== !
715	! !* Update the forcing (BDY and tides)
716	! ! ------------------
717	! Update only tidal forcing at open boundaries
718	IF( ln_bdy .AND. ln_tide ) CALL bdy_dta_tides( kt, kit=jn, time_offset= noffset+1 )
719	IF( ln_tide_pot .AND. ln_tide ) CALL upd_tide ( kt, kit=jn, time_offset= noffset )
720	!
721	! Set extrapolation coefficients for predictor step:
722	IF ((jn<3).AND.ll_init) THEN ! Forward
723	za1 = 1._wp
724	za2 = 0._wp
725	za3 = 0._wp
726	ELSE ! AB3-AM4 Coefficients: bet=0.281105
727	za1 = 1.781105_wp ! za1 = 3/2 + bet
728	za2 = -1.06221_wp ! za2 = -(1/2 + 2*bet)
729	za3 = 0.281105_wp ! za3 = bet
730	ENDIF
731
732	! Extrapolate barotropic velocities at step jit+0.5:
733	ua_e(:,:) = za1 * un_e(:,:) + za2 * ub_e(:,:) + za3 * ubb_e(:,:)
734	va_e(:,:) = za1 * vn_e(:,:) + za2 * vb_e(:,:) + za3 * vbb_e(:,:)
735
736	IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only)
737	! ! ------------------
738	! Extrapolate Sea Level at step jit+0.5:
739	zsshp2_e(:,:) = za1 * sshn_e(:,:) + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:)
740
741	! set wetting & drying mask at tracer points for this barotropic sub-step
742	IF ( ln_wd_dl ) THEN
743
744	IF ( ln_wd_dl_rmp ) THEN
745	DO jj = 1, jpj
746	DO ji = 1, jpi ! vector opt.
747	IF ( zsshp2_e(ji,jj) + ht_0(ji,jj) > 2._wp * rn_wdmin1 ) THEN
748	! IF ( zsshp2_e(ji,jj) + ht_0(ji,jj) > rn_wdmin2 ) THEN
749	ztwdmask(ji,jj) = 1._wp
750	ELSE IF ( zsshp2_e(ji,jj) + ht_0(ji,jj) > rn_wdmin1 ) THEN
751	ztwdmask(ji,jj) = (tanh(50._wp( ( zsshp2_e(ji,jj) + ht_0(ji,jj) - rn_wdmin1 )r_rn_wdmin1)) )
752	ELSE
753	ztwdmask(ji,jj) = 0._wp
754	END IF
755	END DO
756	END DO
757	ELSE
758	DO jj = 1, jpj
759	DO ji = 1, jpi ! vector opt.
760	IF ( zsshp2_e(ji,jj) + ht_0(ji,jj) > rn_wdmin1 ) THEN
761	ztwdmask(ji,jj) = 1._wp
762	ELSE
763	ztwdmask(ji,jj) = 0._wp
764	END IF
765	END DO
766	END DO
767	END IF
768
769	END IF
770
771
772	DO jj = 2, jpjm1 ! Sea Surface Height at u- & v-points
773	DO ji = 2, fs_jpim1 ! Vector opt.
774	zwx(ji,jj) = r1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj) &
775	& * ( e1e2t(ji ,jj) * zsshp2_e(ji ,jj) &
776	& + e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj) )
777	zwy(ji,jj) = r1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj) &
778	& * ( e1e2t(ji,jj ) * zsshp2_e(ji,jj ) &
779	& + e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1) )
780	END DO
781	END DO
782	CALL lbc_lnk_multi( zwx, 'U', 1._wp, zwy, 'V', 1._wp )
783	!
784	zhup2_e (:,:) = hu_0(:,:) + zwx(:,:) ! Ocean depth at U- and V-points
785	zhvp2_e (:,:) = hv_0(:,:) + zwy(:,:)
786	ELSE
787	zhup2_e (:,:) = hu_n(:,:)
788	zhvp2_e (:,:) = hv_n(:,:)
789	ENDIF
790	! !* after ssh
791	! ! -----------
792	! One should enforce volume conservation at open boundaries here
793	! considering fluxes below:
794	!
795	zwx(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:) ! fluxes at jn+0.5
796	zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:)
797	!
798	#if defined key_agrif
799	! Set fluxes during predictor step to ensure volume conservation
800	IF( .NOT.Agrif_Root() .AND. ln_bt_fw ) THEN
801	IF((nbondi == -1).OR.(nbondi == 2)) THEN
802	DO jj = 1, jpj
803	zwx(2:nbghostcells+1,jj) = ubdy_w(jj) * e2u(2:nbghostcells+1,jj)
804	END DO
805	ENDIF
806	IF((nbondi == 1).OR.(nbondi == 2)) THEN
807	DO jj=1,jpj
808	zwx(nlci-nbghostcells-1:nlci-2,jj) = ubdy_e(jj) * e2u(nlci-nbghostcells-1:nlci-2,jj)
809	END DO
810	ENDIF
811	IF((nbondj == -1).OR.(nbondj == 2)) THEN
812	DO ji=1,jpi
813	zwy(ji,2:nbghostcells+1) = vbdy_s(ji) * e1v(ji,2:nbghostcells+1)
814	END DO
815	ENDIF
816	IF((nbondj == 1).OR.(nbondj == 2)) THEN
817	DO ji=1,jpi
818	zwy(ji,nlcj-nbghostcells-1:nlcj-2) = vbdy_n(ji) * e1v(ji,nlcj-nbghostcells-1:nlcj-2)
819	END DO
820	ENDIF
821	ENDIF
822	#endif
823	IF( ln_wd_il ) CALL wad_lmt_bt(zwx, zwy, sshn_e, zssh_frc, rdtbt)
824
825	IF ( ln_wd_dl ) THEN
826
827
828	! un_e and vn_e are set to zero at faces where the direction of the flow is from dry cells
829
830	DO jj = 1, jpjm1
831	DO ji = 1, jpim1
832	IF ( zwx(ji,jj) > 0.0 ) THEN
833	zuwdmask(ji, jj) = ztwdmask(ji ,jj)
834	ELSE
835	zuwdmask(ji, jj) = ztwdmask(ji+1,jj)
836	END IF
837	zwx(ji, jj) = zuwdmask(ji,jj)*zwx(ji, jj)
838	un_e(ji,jj) = zuwdmask(ji,jj)*un_e(ji,jj)
839
840	IF ( zwy(ji,jj) > 0.0 ) THEN
841	zvwdmask(ji, jj) = ztwdmask(ji, jj )
842	ELSE
843	zvwdmask(ji, jj) = ztwdmask(ji, jj+1)
844	END IF
845	zwy(ji, jj) = zvwdmask(ji,jj)*zwy(ji,jj)
846	vn_e(ji,jj) = zvwdmask(ji,jj)*vn_e(ji,jj)
847	END DO
848	END DO
849
850
851	END IF
852
853	! Sum over sub-time-steps to compute advective velocities
854	za2 = wgtbtp2(jn)
855	un_adv(:,:) = un_adv(:,:) + za2 * zwx(:,:) * r1_e2u(:,:)
856	vn_adv(:,:) = vn_adv(:,:) + za2 * zwy(:,:) * r1_e1v(:,:)
857
858	! sum over sub-time-steps to decide which baroclinic velocities to set to zero (zuwdav2 is only used when ln_wd_dl_bc = True)
859	IF ( ln_wd_dl_bc ) THEN
860	zuwdav2(:,:) = zuwdav2(:,:) + za2 * zuwdmask(:,:)
861	zvwdav2(:,:) = zvwdav2(:,:) + za2 * zvwdmask(:,:)
862	END IF
863
864	! Set next sea level:
865	DO jj = 2, jpjm1
866	DO ji = fs_2, fs_jpim1 ! vector opt.
867	zhdiv(ji,jj) = ( zwx(ji,jj) - zwx(ji-1,jj) &
868	& + zwy(ji,jj) - zwy(ji,jj-1) ) * r1_e1e2t(ji,jj)
869	END DO
870	END DO
871	ssha_e(:,:) = ( sshn_e(:,:) - rdtbt * ( zssh_frc(:,:) + zhdiv(:,:) ) ) * ssmask(:,:)
872
873	CALL lbc_lnk( ssha_e, 'T', 1._wp )
874
875	! Duplicate sea level across open boundaries (this is only cosmetic if linssh=T)
876	IF( ln_bdy ) CALL bdy_ssh( ssha_e )
877	#if defined key_agrif
878	IF( .NOT.Agrif_Root() ) CALL agrif_ssh_ts( jn )
879	#endif
880	!
881	! Sea Surface Height at u-,v-points (vvl case only)
882	IF( .NOT.ln_linssh ) THEN
883	DO jj = 2, jpjm1
884	DO ji = 2, jpim1 ! NO Vector Opt.
885	zsshu_a(ji,jj) = r1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj) &
886	& * ( e1e2t(ji ,jj ) * ssha_e(ji ,jj ) &
887	& + e1e2t(ji+1,jj ) * ssha_e(ji+1,jj ) )
888	zsshv_a(ji,jj) = r1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj) &
889	& * ( e1e2t(ji ,jj ) * ssha_e(ji ,jj ) &
890	& + e1e2t(ji ,jj+1) * ssha_e(ji ,jj+1) )
891	END DO
892	END DO
893	CALL lbc_lnk_multi( zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp )
894	ENDIF
895	!
896	! Half-step back interpolation of SSH for surface pressure computation:
897	!----------------------------------------------------------------------
898	IF ((jn==1).AND.ll_init) THEN
899	za0=1._wp ! Forward-backward
900	za1=0._wp
901	za2=0._wp
902	za3=0._wp
903	ELSEIF ((jn==2).AND.ll_init) THEN ! AB2-AM3 Coefficients; bet=0 ; gam=-1/6 ; eps=1/12
904	za0= 1.0833333333333_wp ! za0 = 1-gam-eps
905	za1=-0.1666666666666_wp ! za1 = gam
906	za2= 0.0833333333333_wp ! za2 = eps
907	za3= 0._wp
908	ELSE ! AB3-AM4 Coefficients; bet=0.281105 ; eps=0.013 ; gam=0.0880
909	IF (rn_bt_alpha==0._wp) THEN
910	za0=0.614_wp ! za0 = 1/2 + gam + 2*eps
911	za1=0.285_wp ! za1 = 1/2 - 2gam - 3eps
912	za2=0.088_wp ! za2 = gam
913	za3=0.013_wp ! za3 = eps
914	ELSE
915	zepsilon = 0.00976186_wp - 0.13451357_wp * rn_bt_alpha
916	zgamma = 0.08344500_wp - 0.51358400_wp * rn_bt_alpha
917	za0 = 0.5_wp + zgamma + 2._wp * rn_bt_alpha + 2._wp * zepsilon
918	za1 = 1._wp - za0 - zgamma - zepsilon
919	za2 = zgamma
920	za3 = zepsilon
921	ENDIF
922	ENDIF
923	!
924	zsshp2_e(:,:) = za0 * ssha_e(:,:) + za1 * sshn_e (:,:) &
925	& + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:)
926	IF( ln_wd_il ) THEN ! Calculating and applying W/D gravity filters
927	DO jj = 2, jpjm1
928	DO ji = 2, jpim1
929	ll_tmp1 = MIN( zsshp2_e(ji,jj) , zsshp2_e(ji+1,jj) ) > &
930	& MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) .AND. &
931	& MAX( zsshp2_e(ji,jj) + ht_0(ji,jj) , zsshp2_e(ji+1,jj) + ht_0(ji+1,jj) ) &
932	& > rn_wdmin1 + rn_wdmin2
933	ll_tmp2 = (ABS(zsshp2_e(ji,jj) - zsshp2_e(ji+1,jj)) > 1.E-12 ).AND.( &
934	& MAX( zsshp2_e(ji,jj) , zsshp2_e(ji+1,jj) ) > &
935	& MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 )
936
937	IF(ll_tmp1) THEN
938	zcpx(ji,jj) = 1.0_wp
939	ELSE IF(ll_tmp2) THEN
940	! no worries about zsshp2_e(ji+1,jj) - zsshp2_e(ji ,jj) = 0, it won't happen ! here
941	zcpx(ji,jj) = ABS( (zsshp2_e(ji+1,jj) + ht_0(ji+1,jj) - zsshp2_e(ji,jj) - ht_0(ji,jj)) &
942	& / (zsshp2_e(ji+1,jj) - zsshp2_e(ji ,jj)) )
943	ELSE
944	zcpx(ji,jj) = 0._wp
945	END IF
946
947	ll_tmp1 = MIN( zsshp2_e(ji,jj) , zsshp2_e(ji,jj+1) ) > &
948	& MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) .AND. &
949	& MAX( zsshp2_e(ji,jj) + ht_0(ji,jj) , zsshp2_e(ji,jj+1) + ht_0(ji,jj+1) ) &
950	& > rn_wdmin1 + rn_wdmin2
951	ll_tmp2 = (ABS(zsshp2_e(ji,jj) - zsshp2_e(ji,jj+1)) > 1.E-12 ).AND.( &
952	& MAX( zsshp2_e(ji,jj) , zsshp2_e(ji,jj+1) ) > &
953	& MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 )
954
955	IF(ll_tmp1) THEN
956	zcpy(ji,jj) = 1.0_wp
957	ELSE IF(ll_tmp2) THEN
958	! no worries about zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj ) = 0, it won't happen ! here
959	zcpy(ji,jj) = ABS( (zsshp2_e(ji,jj+1) + ht_0(ji,jj+1) - zsshp2_e(ji,jj) - ht_0(ji,jj)) &
960	& / (zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj )) )
961	ELSE
962	zcpy(ji,jj) = 0._wp
963	END IF
964	END DO
965	END DO
966	END IF
967	!
968	! Compute associated depths at U and V points:
969	IF( .NOT.ln_linssh .AND. .NOT.ln_dynadv_vec ) THEN !* Vector form
970	!
971	DO jj = 2, jpjm1
972	DO ji = 2, jpim1
973	zx1 = r1_2 * ssumask(ji ,jj) * r1_e1e2u(ji ,jj) &
974	& * ( e1e2t(ji ,jj ) * zsshp2_e(ji ,jj) &
975	& + e1e2t(ji+1,jj ) * zsshp2_e(ji+1,jj ) )
976	zy1 = r1_2 * ssvmask(ji ,jj) * r1_e1e2v(ji ,jj ) &
977	& * ( e1e2t(ji ,jj ) * zsshp2_e(ji ,jj ) &
978	& + e1e2t(ji ,jj+1) * zsshp2_e(ji ,jj+1) )
979	zhust_e(ji,jj) = hu_0(ji,jj) + zx1
980	zhvst_e(ji,jj) = hv_0(ji,jj) + zy1
981	END DO
982	END DO
983
984	ENDIF
985	!
986	! Add Coriolis trend:
987	! zwz array below or triads normally depend on sea level with ln_linssh=F and should be updated
988	! at each time step. We however keep them constant here for optimization.
989	! Recall that zwx and zwy arrays hold fluxes at this stage:
990	! zwx(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:) ! fluxes at jn+0.5
991	! zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:)
992	!
993	IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN !== energy conserving or mixed scheme ==!
994	DO jj = 2, jpjm1
995	DO ji = fs_2, fs_jpim1 ! vector opt.
996	zy1 = ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) ) * r1_e1u(ji,jj)
997	zy2 = ( zwy(ji ,jj ) + zwy(ji+1,jj ) ) * r1_e1u(ji,jj)
998	zx1 = ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) ) * r1_e2v(ji,jj)
999	zx2 = ( zwx(ji ,jj ) + zwx(ji ,jj+1) ) * r1_e2v(ji,jj)
1000	zu_trd(ji,jj) = r1_4 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
1001	zv_trd(ji,jj) =-r1_4 * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 )
1002	END DO
1003	END DO
1004	!
1005	ELSEIF ( ln_dynvor_ens ) THEN !== enstrophy conserving scheme ==!
1006	DO jj = 2, jpjm1
1007	DO ji = fs_2, fs_jpim1 ! vector opt.
1008	zy1 = r1_8 * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) &
1009	& + zwy(ji ,jj ) + zwy(ji+1,jj ) ) * r1_e1u(ji,jj)
1010	zx1 = - r1_8 * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) &
1011	& + zwx(ji ,jj ) + zwx(ji ,jj+1) ) * r1_e2v(ji,jj)
1012	zu_trd(ji,jj) = zy1 * ( zwz(ji ,jj-1) + zwz(ji,jj) )
1013	zv_trd(ji,jj) = zx1 * ( zwz(ji-1,jj ) + zwz(ji,jj) )
1014	END DO
1015	END DO
1016	!
1017	ELSEIF ( ln_dynvor_een ) THEN !== energy and enstrophy conserving scheme ==!
1018	DO jj = 2, jpjm1
1019	DO ji = fs_2, fs_jpim1 ! vector opt.
1020	zu_trd(ji,jj) = + r1_12 * r1_e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) &
1021	& + ftnw(ji+1,jj) * zwy(ji+1,jj ) &
1022	& + ftse(ji,jj ) * zwy(ji ,jj-1) &
1023	& + ftsw(ji+1,jj) * zwy(ji+1,jj-1) )
1024	zv_trd(ji,jj) = - r1_12 * r1_e2v(ji,jj) * ( ftsw(ji,jj+1) * zwx(ji-1,jj+1) &
1025	& + ftse(ji,jj+1) * zwx(ji ,jj+1) &
1026	& + ftnw(ji,jj ) * zwx(ji-1,jj ) &
1027	& + ftne(ji,jj ) * zwx(ji ,jj ) )
1028	END DO
1029	END DO
1030	!
1031	ENDIF
1032	!
1033	! Add tidal astronomical forcing if defined
1034	IF ( ln_tide .AND. ln_tide_pot ) THEN
1035	DO jj = 2, jpjm1
1036	DO ji = fs_2, fs_jpim1 ! vector opt.
1037	zu_spg = grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) * r1_e1u(ji,jj)
1038	zv_spg = grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) * r1_e2v(ji,jj)
1039	zu_trd(ji,jj) = zu_trd(ji,jj) + zu_spg
1040	zv_trd(ji,jj) = zv_trd(ji,jj) + zv_spg
1041	END DO
1042	END DO
1043	ENDIF
1044	!
1045	! Add bottom stresses:
1046	!jth do implicitly instead
1047	IF ( .NOT. ll_wd ) THEN ! Revert to explicit for bit comparison tests in non wad runs
1048	DO jj = 2, jpjm1
1049	DO ji = fs_2, fs_jpim1 ! vector opt.
1050	zu_trd(ji,jj) = zu_trd(ji,jj) + zCdU_u(ji,jj) * un_e(ji,jj) * hur_e(ji,jj)
1051	zv_trd(ji,jj) = zv_trd(ji,jj) + zCdU_v(ji,jj) * vn_e(ji,jj) * hvr_e(ji,jj)
1052	END DO
1053	END DO
1054	ENDIF
1055	!
1056	! Surface pressure trend:
1057	IF( ln_wd_il ) THEN
1058	DO jj = 2, jpjm1
1059	DO ji = 2, jpim1
1060	! Add surface pressure gradient
1061	zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj)
1062	zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj)
1063	zwx(ji,jj) = (1._wp - rn_scal_load) * zu_spg * zcpx(ji,jj)
1064	zwy(ji,jj) = (1._wp - rn_scal_load) * zv_spg * zcpy(ji,jj)
1065	END DO
1066	END DO
1067	ELSE
1068	DO jj = 2, jpjm1
1069	DO ji = fs_2, fs_jpim1 ! vector opt.
1070	! Add surface pressure gradient
1071	zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj)
1072	zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj)
1073	zwx(ji,jj) = (1._wp - rn_scal_load) * zu_spg
1074	zwy(ji,jj) = (1._wp - rn_scal_load) * zv_spg
1075	END DO
1076	END DO
1077	END IF
1078
1079	!
1080	! Set next velocities:
1081	IF( ln_dynadv_vec .OR. ln_linssh ) THEN !* Vector form
1082	DO jj = 2, jpjm1
1083	DO ji = fs_2, fs_jpim1 ! vector opt.
1084	ua_e(ji,jj) = ( un_e(ji,jj) &
1085	& + rdtbt * ( zwx(ji,jj) &
1086	& + zu_trd(ji,jj) &
1087	& + zu_frc(ji,jj) ) &
1088	& ) * ssumask(ji,jj)
1089
1090	va_e(ji,jj) = ( vn_e(ji,jj) &
1091	& + rdtbt * ( zwy(ji,jj) &
1092	& + zv_trd(ji,jj) &
1093	& + zv_frc(ji,jj) ) &
1094	& ) * ssvmask(ji,jj)
1095
1096	!jth implicit bottom friction:
1097	IF ( ll_wd ) THEN ! revert to explicit for bit comparison tests in non wad runs
1098	ua_e(ji,jj) = ua_e(ji,jj) /(1.0 - rdtbt * zCdU_u(ji,jj) * hur_e(ji,jj))
1099	va_e(ji,jj) = va_e(ji,jj) /(1.0 - rdtbt * zCdU_v(ji,jj) * hvr_e(ji,jj))
1100	ENDIF
1101
1102	END DO
1103	END DO
1104	!
1105	ELSE !* Flux form
1106	DO jj = 2, jpjm1
1107	DO ji = fs_2, fs_jpim1 ! vector opt.
1108
1109	zhura = hu_0(ji,jj) + zsshu_a(ji,jj)
1110	zhvra = hv_0(ji,jj) + zsshv_a(ji,jj)
1111
1112	zhura = ssumask(ji,jj)/(zhura + 1._wp - ssumask(ji,jj))
1113	zhvra = ssvmask(ji,jj)/(zhvra + 1._wp - ssvmask(ji,jj))
1114
1115	ua_e(ji,jj) = ( hu_e(ji,jj) * un_e(ji,jj) &
1116	& + rdtbt * ( zhust_e(ji,jj) * zwx(ji,jj) &
1117	& + zhup2_e(ji,jj) * zu_trd(ji,jj) &
1118	& + hu_n(ji,jj) * zu_frc(ji,jj) ) &
1119	& ) * zhura
1120
1121	va_e(ji,jj) = ( hv_e(ji,jj) * vn_e(ji,jj) &
1122	& + rdtbt * ( zhvst_e(ji,jj) * zwy(ji,jj) &
1123	& + zhvp2_e(ji,jj) * zv_trd(ji,jj) &
1124	& + hv_n(ji,jj) * zv_frc(ji,jj) ) &
1125	& ) * zhvra
1126	END DO
1127	END DO
1128	ENDIF
1129
1130
1131	IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only)
1132	hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:)
1133	hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:)
1134	hur_e(:,:) = ssumask(:,:) / ( hu_e(:,:) + 1._wp - ssumask(:,:) )
1135	hvr_e(:,:) = ssvmask(:,:) / ( hv_e(:,:) + 1._wp - ssvmask(:,:) )
1136	!
1137	ENDIF
1138	! !* domain lateral boundary
1139	CALL lbc_lnk_multi( ua_e, 'U', -1._wp, va_e , 'V', -1._wp )
1140	!
1141	! ! open boundaries
1142	IF( ln_bdy ) CALL bdy_dyn2d( jn, ua_e, va_e, un_e, vn_e, hur_e, hvr_e, ssha_e )
1143	#if defined key_agrif
1144	IF( .NOT.Agrif_Root() ) CALL agrif_dyn_ts( jn ) ! Agrif
1145	#endif
1146	! !* Swap
1147	! ! ----
1148	ubb_e (:,:) = ub_e (:,:)
1149	ub_e (:,:) = un_e (:,:)
1150	un_e (:,:) = ua_e (:,:)
1151	!
1152	vbb_e (:,:) = vb_e (:,:)
1153	vb_e (:,:) = vn_e (:,:)
1154	vn_e (:,:) = va_e (:,:)
1155	!
1156	sshbb_e(:,:) = sshb_e(:,:)
1157	sshb_e (:,:) = sshn_e(:,:)
1158	sshn_e (:,:) = ssha_e(:,:)
1159
1160	! !* Sum over whole bt loop
1161	! ! ----------------------
1162	za1 = wgtbtp1(jn)
1163	IF( ln_dynadv_vec .OR. ln_linssh ) THEN ! Sum velocities
1164	ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:)
1165	va_b (:,:) = va_b (:,:) + za1 * va_e (:,:)
1166	ELSE ! Sum transports
1167	IF ( .NOT.ln_wd_dl ) THEN
1168	ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:) * hu_e (:,:)
1169	va_b (:,:) = va_b (:,:) + za1 * va_e (:,:) * hv_e (:,:)
1170	ELSE
1171	ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:) * hu_e (:,:) * zuwdmask(:,:)
1172	va_b (:,:) = va_b (:,:) + za1 * va_e (:,:) * hv_e (:,:) * zvwdmask(:,:)
1173	END IF
1174	ENDIF
1175	! ! Sum sea level
1176	ssha(:,:) = ssha(:,:) + za1 * ssha_e(:,:)
1177
1178	! ! ==================== !
1179	END DO ! end loop !
1180	! ! ==================== !
1181	! -----------------------------------------------------------------------------
1182	! Phase 3. update the general trend with the barotropic trend
1183	! -----------------------------------------------------------------------------
1184	!
1185	! Set advection velocity correction:
1186	IF (ln_bt_fw) THEN
1187	zwx(:,:) = un_adv(:,:)
1188	zwy(:,:) = vn_adv(:,:)
1189	IF( .NOT.( kt == nit000 .AND. neuler==0 ) ) THEN
1190	un_adv(:,:) = r1_2 * ( ub2_b(:,:) + zwx(:,:) - atfp * un_bf(:,:) )
1191	vn_adv(:,:) = r1_2 * ( vb2_b(:,:) + zwy(:,:) - atfp * vn_bf(:,:) )
1192	!
1193	! Update corrective fluxes for next time step:
1194	un_bf(:,:) = atfp * un_bf(:,:) + (zwx(:,:) - ub2_b(:,:))
1195	vn_bf(:,:) = atfp * vn_bf(:,:) + (zwy(:,:) - vb2_b(:,:))
1196	ELSE
1197	un_bf(:,:) = 0._wp
1198	vn_bf(:,:) = 0._wp
1199	END IF
1200	! Save integrated transport for next computation
1201	ub2_b(:,:) = zwx(:,:)
1202	vb2_b(:,:) = zwy(:,:)
1203	ENDIF
1204
1205
1206	!
1207	! Update barotropic trend:
1208	IF( ln_dynadv_vec .OR. ln_linssh ) THEN
1209	DO jk=1,jpkm1
1210	ua(:,:,jk) = ua(:,:,jk) + ( ua_b(:,:) - ub_b(:,:) ) * r1_2dt_b
1211	va(:,:,jk) = va(:,:,jk) + ( va_b(:,:) - vb_b(:,:) ) * r1_2dt_b
1212	END DO
1213	ELSE
1214	! At this stage, ssha has been corrected: compute new depths at velocity points
1215	DO jj = 1, jpjm1
1216	DO ji = 1, jpim1 ! NO Vector Opt.
1217	zsshu_a(ji,jj) = r1_2 * umask(ji,jj,1) * r1_e1e2u(ji,jj) &
1218	& * ( e1e2t(ji ,jj) * ssha(ji ,jj) &
1219	& + e1e2t(ji+1,jj) * ssha(ji+1,jj) )
1220	zsshv_a(ji,jj) = r1_2 * vmask(ji,jj,1) * r1_e1e2v(ji,jj) &
1221	& * ( e1e2t(ji,jj ) * ssha(ji,jj ) &
1222	& + e1e2t(ji,jj+1) * ssha(ji,jj+1) )
1223	END DO
1224	END DO
1225	CALL lbc_lnk_multi( zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) ! Boundary conditions
1226	!
1227	DO jk=1,jpkm1
1228	ua(:,:,jk) = ua(:,:,jk) + r1_hu_n(:,:) * ( ua_b(:,:) - ub_b(:,:) * hu_b(:,:) ) * r1_2dt_b
1229	va(:,:,jk) = va(:,:,jk) + r1_hv_n(:,:) * ( va_b(:,:) - vb_b(:,:) * hv_b(:,:) ) * r1_2dt_b
1230	END DO
1231	! Save barotropic velocities not transport:
1232	ua_b(:,:) = ua_b(:,:) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - ssumask(:,:) )
1233	va_b(:,:) = va_b(:,:) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - ssvmask(:,:) )
1234	ENDIF
1235
1236
1237	! Correct velocities so that the barotropic velocity equals (un_adv, vn_adv) (in all cases)
1238	DO jk = 1, jpkm1
1239	un(:,:,jk) = ( un(:,:,jk) + un_adv(:,:)r1_hu_n(:,:) - un_b(:,:) ) umask(:,:,jk)
1240	vn(:,:,jk) = ( vn(:,:,jk) + vn_adv(:,:)r1_hv_n(:,:) - vn_b(:,:) ) vmask(:,:,jk)
1241	END DO
1242
1243	IF ( ln_wd_dl .and. ln_wd_dl_bc) THEN
1244	DO jk = 1, jpkm1
1245	un(:,:,jk) = ( un_adv(:,:)*r1_hu_n(:,:) &
1246	& + zuwdav2(:,:)(un(:,:,jk) - un_adv(:,:)r1_hu_n(:,:)) ) * umask(:,:,jk)
1247	vn(:,:,jk) = ( vn_adv(:,:)*r1_hv_n(:,:) &
1248	& + zvwdav2(:,:)(vn(:,:,jk) - vn_adv(:,:)r1_hv_n(:,:)) ) * vmask(:,:,jk)
1249	END DO
1250	END IF
1251
1252
1253	CALL iom_put( "ubar", un_adv(:,:)*r1_hu_n(:,:) ) ! barotropic i-current
1254	CALL iom_put( "vbar", vn_adv(:,:)*r1_hv_n(:,:) ) ! barotropic i-current
1255	!
1256	#if defined key_agrif
1257	! Save time integrated fluxes during child grid integration
1258	! (used to update coarse grid transports at next time step)
1259	!
1260	IF( .NOT.Agrif_Root() .AND. ln_bt_fw ) THEN
1261	IF( Agrif_NbStepint() == 0 ) THEN
1262	ub2_i_b(:,:) = 0._wp
1263	vb2_i_b(:,:) = 0._wp
1264	END IF
1265	!
1266	za1 = 1._wp / REAL(Agrif_rhot(), wp)
1267	ub2_i_b(:,:) = ub2_i_b(:,:) + za1 * ub2_b(:,:)
1268	vb2_i_b(:,:) = vb2_i_b(:,:) + za1 * vb2_b(:,:)
1269	ENDIF
1270	#endif
1271	! !* write time-spliting arrays in the restart
1272	IF( lrst_oce .AND.ln_bt_fw ) CALL ts_rst( kt, 'WRITE' )
1273	!
1274	IF( ln_wd_il ) DEALLOCATE( zcpx, zcpy )
1275	IF( ln_wd_dl ) DEALLOCATE( ztwdmask, zuwdmask, zvwdmask, zuwdav2, zvwdav2 )
1276	!
1277	IF( ln_diatmb ) THEN
1278	CALL iom_put( "baro_u" , un_bumask(:,:,1)+zmdi(1-umask(:,:,1 ) ) ) ! Barotropic U Velocity
1279	CALL iom_put( "baro_v" , vn_bvmask(:,:,1)+zmdi(1-vmask(:,:,1 ) ) ) ! Barotropic V Velocity
1280	ENDIF
1281	IF( ln_timing ) CALL timing_stop('dyn_spg_ts')
1282	!
1283	END SUBROUTINE dyn_spg_ts
1284
1285
1286	SUBROUTINE ts_wgt( ll_av, ll_fw, jpit, zwgt1, zwgt2)
1287	!!---------------------------------------------------------------------
1288	!! * ROUTINE ts_wgt *
1289	!!
1290	!! ** Purpose : Set time-splitting weights for temporal averaging (or not)
1291	!!----------------------------------------------------------------------
1292	LOGICAL, INTENT(in) :: ll_av ! temporal averaging=.true.
1293	LOGICAL, INTENT(in) :: ll_fw ! forward time splitting =.true.
1294	INTEGER, INTENT(inout) :: jpit ! cycle length
1295	REAL(wp), DIMENSION(3*nn_baro), INTENT(inout) :: zwgt1, & ! Primary weights
1296	zwgt2 ! Secondary weights
1297
1298	INTEGER :: jic, jn, ji ! temporary integers
1299	REAL(wp) :: za1, za2
1300	!!----------------------------------------------------------------------
1301
1302	zwgt1(:) = 0._wp
1303	zwgt2(:) = 0._wp
1304
1305	! Set time index when averaged value is requested
1306	IF (ll_fw) THEN
1307	jic = nn_baro
1308	ELSE
1309	jic = 2 * nn_baro
1310	ENDIF
1311
1312	! Set primary weights:
1313	IF (ll_av) THEN
1314	! Define simple boxcar window for primary weights
1315	! (width = nn_baro, centered around jic)
1316	SELECT CASE ( nn_bt_flt )
1317	CASE( 0 ) ! No averaging
1318	zwgt1(jic) = 1._wp
1319	jpit = jic
1320
1321	CASE( 1 ) ! Boxcar, width = nn_baro
1322	DO jn = 1, 3*nn_baro
1323	za1 = ABS(float(jn-jic))/float(nn_baro)
1324	IF (za1 < 0.5_wp) THEN
1325	zwgt1(jn) = 1._wp
1326	jpit = jn
1327	ENDIF
1328	ENDDO
1329
1330	CASE( 2 ) ! Boxcar, width = 2 * nn_baro
1331	DO jn = 1, 3*nn_baro
1332	za1 = ABS(float(jn-jic))/float(nn_baro)
1333	IF (za1 < 1._wp) THEN
1334	zwgt1(jn) = 1._wp
1335	jpit = jn
1336	ENDIF
1337	ENDDO
1338	CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_bt_flt' )
1339	END SELECT
1340
1341	ELSE ! No time averaging
1342	zwgt1(jic) = 1._wp
1343	jpit = jic
1344	ENDIF
1345
1346	! Set secondary weights
1347	DO jn = 1, jpit
1348	DO ji = jn, jpit
1349	zwgt2(jn) = zwgt2(jn) + zwgt1(ji)
1350	END DO
1351	END DO
1352
1353	! Normalize weigths:
1354	za1 = 1._wp / SUM(zwgt1(1:jpit))
1355	za2 = 1._wp / SUM(zwgt2(1:jpit))
1356	DO jn = 1, jpit
1357	zwgt1(jn) = zwgt1(jn) * za1
1358	zwgt2(jn) = zwgt2(jn) * za2
1359	END DO
1360	!
1361	END SUBROUTINE ts_wgt
1362
1363
1364	SUBROUTINE ts_rst( kt, cdrw )
1365	!!---------------------------------------------------------------------
1366	!! * ROUTINE ts_rst *
1367	!!
1368	!! ** Purpose : Read or write time-splitting arrays in restart file
1369	!!----------------------------------------------------------------------
1370	INTEGER , INTENT(in) :: kt ! ocean time-step
1371	CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag
1372	!
1373	!!----------------------------------------------------------------------
1374	!
1375	IF( TRIM(cdrw) == 'READ' ) THEN
1376	CALL iom_get( numror, jpdom_autoglo, 'ub2_b' , ub2_b (:,:) )
1377	CALL iom_get( numror, jpdom_autoglo, 'vb2_b' , vb2_b (:,:) )
1378	CALL iom_get( numror, jpdom_autoglo, 'un_bf' , un_bf (:,:) )
1379	CALL iom_get( numror, jpdom_autoglo, 'vn_bf' , vn_bf (:,:) )
1380	IF( .NOT.ln_bt_av ) THEN
1381	CALL iom_get( numror, jpdom_autoglo, 'sshbb_e' , sshbb_e(:,:) )
1382	CALL iom_get( numror, jpdom_autoglo, 'ubb_e' , ubb_e(:,:) )
1383	CALL iom_get( numror, jpdom_autoglo, 'vbb_e' , vbb_e(:,:) )
1384	CALL iom_get( numror, jpdom_autoglo, 'sshb_e' , sshb_e(:,:) )
1385	CALL iom_get( numror, jpdom_autoglo, 'ub_e' , ub_e(:,:) )
1386	CALL iom_get( numror, jpdom_autoglo, 'vb_e' , vb_e(:,:) )
1387	ENDIF
1388	#if defined key_agrif
1389	! Read time integrated fluxes
1390	IF ( .NOT.Agrif_Root() ) THEN
1391	CALL iom_get( numror, jpdom_autoglo, 'ub2_i_b' , ub2_i_b(:,:) )
1392	CALL iom_get( numror, jpdom_autoglo, 'vb2_i_b' , vb2_i_b(:,:) )
1393	ENDIF
1394	#endif
1395	!
1396	ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN
1397	CALL iom_rstput( kt, nitrst, numrow, 'ub2_b' , ub2_b (:,:) )
1398	CALL iom_rstput( kt, nitrst, numrow, 'vb2_b' , vb2_b (:,:) )
1399	CALL iom_rstput( kt, nitrst, numrow, 'un_bf' , un_bf (:,:) )
1400	CALL iom_rstput( kt, nitrst, numrow, 'vn_bf' , vn_bf (:,:) )
1401	!
1402	IF (.NOT.ln_bt_av) THEN
1403	CALL iom_rstput( kt, nitrst, numrow, 'sshbb_e' , sshbb_e(:,:) )
1404	CALL iom_rstput( kt, nitrst, numrow, 'ubb_e' , ubb_e(:,:) )
1405	CALL iom_rstput( kt, nitrst, numrow, 'vbb_e' , vbb_e(:,:) )
1406	CALL iom_rstput( kt, nitrst, numrow, 'sshb_e' , sshb_e(:,:) )
1407	CALL iom_rstput( kt, nitrst, numrow, 'ub_e' , ub_e(:,:) )
1408	CALL iom_rstput( kt, nitrst, numrow, 'vb_e' , vb_e(:,:) )
1409	ENDIF
1410	#if defined key_agrif
1411	! Save time integrated fluxes
1412	IF ( .NOT.Agrif_Root() ) THEN
1413	CALL iom_rstput( kt, nitrst, numrow, 'ub2_i_b' , ub2_i_b(:,:) )
1414	CALL iom_rstput( kt, nitrst, numrow, 'vb2_i_b' , vb2_i_b(:,:) )
1415	ENDIF
1416	#endif
1417	ENDIF
1418	!
1419	END SUBROUTINE ts_rst
1420
1421
1422	SUBROUTINE dyn_spg_ts_init
1423	!!---------------------------------------------------------------------
1424	!! * ROUTINE dyn_spg_ts_init *
1425	!!
1426	!! ** Purpose : Set time splitting options
1427	!!----------------------------------------------------------------------
1428	INTEGER :: ji ,jj ! dummy loop indices
1429	REAL(wp) :: zxr2, zyr2, zcmax ! local scalar
1430	REAL(wp), DIMENSION(jpi,jpj) :: zcu
1431	!!----------------------------------------------------------------------
1432	!
1433	! Max courant number for ext. grav. waves
1434	!
1435	DO jj = 1, jpj
1436	DO ji =1, jpi
1437	zxr2 = r1_e1t(ji,jj) * r1_e1t(ji,jj)
1438	zyr2 = r1_e2t(ji,jj) * r1_e2t(ji,jj)
1439	zcu(ji,jj) = SQRT( grav * MAX(ht_0(ji,jj),0._wp) * (zxr2 + zyr2) )
1440	END DO
1441	END DO
1442	!
1443	zcmax = MAXVAL( zcu(:,:) )
1444	IF( lk_mpp ) CALL mpp_max( zcmax )
1445
1446	! Estimate number of iterations to satisfy a max courant number= rn_bt_cmax
1447	IF( ln_bt_auto ) nn_baro = CEILING( rdt / rn_bt_cmax * zcmax)
1448
1449	rdtbt = rdt / REAL( nn_baro , wp )
1450	zcmax = zcmax * rdtbt
1451	! Print results
1452	IF(lwp) WRITE(numout,*)
1453	IF(lwp) WRITE(numout,*) 'dyn_spg_ts : split-explicit free surface'
1454	IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
1455	IF( ln_bt_auto ) THEN
1456	IF(lwp) WRITE(numout,*) ' ln_ts_auto=.true. Automatically set nn_baro '
1457	IF(lwp) WRITE(numout,*) ' Max. courant number allowed: ', rn_bt_cmax
1458	ELSE
1459	IF(lwp) WRITE(numout,*) ' ln_ts_auto=.false.: Use nn_baro in namelist '
1460	ENDIF
1461
1462	IF(ln_bt_av) THEN
1463	IF(lwp) WRITE(numout,*) ' ln_bt_av=.true. => Time averaging over nn_baro time steps is on '
1464	ELSE
1465	IF(lwp) WRITE(numout,*) ' ln_bt_av=.false. => No time averaging of barotropic variables '
1466	ENDIF
1467	!
1468	!
1469	IF(ln_bt_fw) THEN
1470	IF(lwp) WRITE(numout,*) ' ln_bt_fw=.true. => Forward integration of barotropic variables '
1471	ELSE
1472	IF(lwp) WRITE(numout,*) ' ln_bt_fw =.false.=> Centred integration of barotropic variables '
1473	ENDIF
1474	!
1475	#if defined key_agrif
1476	! Restrict the use of Agrif to the forward case only
1477	!!! IF( .NOT.ln_bt_fw .AND. .NOT.Agrif_Root() ) CALL ctl_stop( 'AGRIF not implemented if ln_bt_fw=.FALSE.' )
1478	#endif
1479	!
1480	IF(lwp) WRITE(numout,*) ' Time filter choice, nn_bt_flt: ', nn_bt_flt
1481	SELECT CASE ( nn_bt_flt )
1482	CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' Dirac'
1483	CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' Boxcar: width = nn_baro'
1484	CASE( 2 ) ; IF(lwp) WRITE(numout,) ' Boxcar: width = 2nn_baro'
1485	CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_bt_flt: should 0,1,2' )
1486	END SELECT
1487	!
1488	IF(lwp) WRITE(numout,*) ' '
1489	IF(lwp) WRITE(numout,*) ' nn_baro = ', nn_baro
1490	IF(lwp) WRITE(numout,*) ' Barotropic time step [s] is :', rdtbt
1491	IF(lwp) WRITE(numout,*) ' Maximum Courant number is :', zcmax
1492	!
1493	IF(lwp) WRITE(numout,*) ' Time diffusion parameter rn_bt_alpha: ', rn_bt_alpha
1494	IF ((ln_bt_av.AND.nn_bt_flt/=0).AND.(rn_bt_alpha>0._wp)) THEN
1495	CALL ctl_stop( 'dynspg_ts ERROR: if rn_bt_alpha > 0, remove temporal averaging' )
1496	ENDIF
1497	!
1498	IF( .NOT.ln_bt_av .AND. .NOT.ln_bt_fw ) THEN
1499	CALL ctl_stop( 'dynspg_ts ERROR: No time averaging => only forward integration is possible' )
1500	ENDIF
1501	IF( zcmax>0.9_wp ) THEN
1502	CALL ctl_stop( 'dynspg_ts ERROR: Maximum Courant number is greater than 0.9: Inc. nn_baro !' )
1503	ENDIF
1504	!
1505	END SUBROUTINE dyn_spg_ts_init
1506
1507	!!======================================================================
1508	END MODULE dynspg_ts

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

Download in other formats: