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

Context Navigation

source: branches/UKMO/ROMS_WAD_7832/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 8870

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

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

Download in other formats: