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

Context Navigation

source: branches/2017/dev_METO_2017/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 8987

Last change on this file since 8987 was 8987, checked in by timgraham, 6 years ago
Merge in XIOS read/write branch again
Property svn:keywords set to `Id`
File size: 69.4 KB

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

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

Download in other formats: