Context Navigation

dynspg_ts.F90 @ 15194

Last change on this file since 15194 was 15194, checked in by davestorkey, 3 years ago
UKMO/NEMO_4.0.4_momentum_trends: Bug fixes so that it works with explicit friction (ln_drgimp=F).
File size: 88.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	!! 4.0 ! 2017-05 (G. Madec) drag coef. defined at t-point (zdfdrg.F90)
22	!!---------------------------------------------------------------------
23
24	!!----------------------------------------------------------------------
25	!! dyn_spg_ts : compute surface pressure gradient trend using a time-splitting scheme
26	!! dyn_spg_ts_init: initialisation of the time-splitting scheme
27	!! ts_wgt : set time-splitting weights for temporal averaging (or not)
28	!! ts_rst : read/write time-splitting fields in restart file
29	!!----------------------------------------------------------------------
30	USE oce ! ocean dynamics and tracers
31	USE dom_oce ! ocean space and time domain
32	USE sbc_oce ! surface boundary condition: ocean
33	USE zdf_oce ! vertical physics: variables
34	USE zdfdrg ! vertical physics: top/bottom drag coef.
35	USE sbcisf ! ice shelf variable (fwfisf)
36	USE sbcapr ! surface boundary condition: atmospheric pressure
37	USE dynadv , ONLY: ln_dynadv_vec
38	USE dynvor ! vortivity scheme indicators
39	USE phycst ! physical constants
40	USE dynvor ! vorticity term
41	USE wet_dry ! wetting/drying flux limter
42	USE bdy_oce ! open boundary
43	USE bdyvol ! open boundary volume conservation
44	USE bdytides ! open boundary condition data
45	USE bdydyn2d ! open boundary conditions on barotropic variables
46	USE sbctide ! tides
47	USE updtide ! tide potential
48	USE sbcwave ! surface wave
49	#if defined key_agrif
50	USE agrif_oce_interp ! agrif
51	USE agrif_oce
52	#endif
53	#if defined key_asminc
54	USE asminc ! Assimilation increment
55	#endif
56	!
57	USE in_out_manager ! I/O manager
58	USE lib_mpp ! distributed memory computing library
59	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
60	USE prtctl ! Print control
61	USE iom ! IOM library
62	USE restart ! only for lrst_oce
63	USE trd_oce ! trends: ocean variables
64	USE trddyn ! trend manager: dynamics
65
66	USE iom ! to remove
67
68	IMPLICIT NONE
69	PRIVATE
70
71	PUBLIC dyn_spg_ts ! called by dyn_spg
72	PUBLIC dyn_spg_ts_init ! - - dyn_spg_init
73
74	!! Time filtered arrays at baroclinic time step:
75	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: un_adv , vn_adv !: Advection vel. at "now" barocl. step
76	!
77	INTEGER, SAVE :: icycle ! Number of barotropic sub-steps for each internal step nn_baro <= 2.5 nn_baro
78	REAL(wp),SAVE :: rdtbt ! Barotropic time step
79	!
80	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: wgtbtp1, wgtbtp2 ! 1st & 2nd weights used in time filtering of barotropic fields
81	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zwz ! ff_f/h at F points
82	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ftnw, ftne ! triad of coriolis parameter
83	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ftsw, ftse ! (only used with een vorticity scheme)
84
85	REAL(wp) :: r1_12 = 1._wp / 12._wp ! local ratios
86	REAL(wp) :: r1_8 = 0.125_wp !
87	REAL(wp) :: r1_4 = 0.25_wp !
88	REAL(wp) :: r1_2 = 0.5_wp !
89
90	!! * Substitutions
91	# include "vectopt_loop_substitute.h90"
92	!!----------------------------------------------------------------------
93	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
94	!! $Id$
95	!! Software governed by the CeCILL license (see ./LICENSE)
96	!!----------------------------------------------------------------------
97	CONTAINS
98
99	INTEGER FUNCTION dyn_spg_ts_alloc()
100	!!----------------------------------------------------------------------
101	!! * routine dyn_spg_ts_alloc *
102	!!----------------------------------------------------------------------
103	INTEGER :: ierr(3)
104	!!----------------------------------------------------------------------
105	ierr(:) = 0
106	!
107	ALLOCATE( wgtbtp1(3nn_baro), wgtbtp2(3nn_baro), zwz(jpi,jpj), STAT=ierr(1) )
108	IF( ln_dynvor_een .OR. ln_dynvor_eeT ) &
109	& ALLOCATE( ftnw(jpi,jpj) , ftne(jpi,jpj) , ftsw(jpi,jpj) , ftse(jpi,jpj), STAT=ierr(2) )
110	!
111	ALLOCATE( un_adv(jpi,jpj), vn_adv(jpi,jpj) , STAT=ierr(3) )
112	!
113	dyn_spg_ts_alloc = MAXVAL( ierr(:) )
114	!
115	CALL mpp_sum( 'dynspg_ts', dyn_spg_ts_alloc )
116	IF( dyn_spg_ts_alloc /= 0 ) CALL ctl_stop( 'STOP', 'dyn_spg_ts_alloc: failed to allocate arrays' )
117	!
118	END FUNCTION dyn_spg_ts_alloc
119
120
121	SUBROUTINE dyn_spg_ts( kt )
122	!!----------------------------------------------------------------------
123	!!
124	!! ** Purpose : - Compute the now trend due to the explicit time stepping
125	!! of the quasi-linear barotropic system, and add it to the
126	!! general momentum trend.
127	!!
128	!! ** Method : - split-explicit schem (time splitting) :
129	!! Barotropic variables are advanced from internal time steps
130	!! "n" to "n+1" if ln_bt_fw=T
131	!! or from
132	!! "n-1" to "n+1" if ln_bt_fw=F
133	!! thanks to a generalized forward-backward time stepping (see ref. below).
134	!!
135	!! ** Action :
136	!! -Update the filtered free surface at step "n+1" : ssha
137	!! -Update filtered barotropic velocities at step "n+1" : ua_b, va_b
138	!! -Compute barotropic advective fluxes at step "n" : un_adv, vn_adv
139	!! These are used to advect tracers and are compliant with discrete
140	!! continuity equation taken at the baroclinic time steps. This
141	!! ensures tracers conservation.
142	!! - (ua, va) momentum trend updated with barotropic component.
143	!!
144	!! References : Shchepetkin and McWilliams, Ocean Modelling, 2005.
145	!!---------------------------------------------------------------------
146	INTEGER, INTENT(in) :: kt ! ocean time-step index
147	!
148	INTEGER :: ji, jj, jk, jn ! dummy loop indices
149	LOGICAL :: ll_fw_start ! =T : forward integration
150	LOGICAL :: ll_init ! =T : special startup of 2d equations
151	INTEGER :: noffset ! local integers : time offset for bdy update
152	REAL(wp) :: r1_2dt_b, z1_hu, z1_hv ! local scalars
153	REAL(wp) :: za0, za1, za2, za3 ! - -
154	REAL(wp) :: zmdi, zztmp, zldg ! - -
155	REAL(wp) :: zhu_bck, zhv_bck, zhdiv ! - -
156	REAL(wp) :: zun_save, zvn_save ! - -
157	REAL(wp), DIMENSION(jpi,jpj) :: zu_trd, zu_frc, zu_spg, zssh_frc
158	REAL(wp), DIMENSION(jpi,jpj) :: zv_trd, zv_frc, zv_spg
159	REAL(wp), DIMENSION(jpi,jpj) :: zsshu_a, zhup2_e, zhtp2_e
160	REAL(wp), DIMENSION(jpi,jpj) :: zsshv_a, zhvp2_e, zsshp2_e
161	REAL(wp), DIMENSION(jpi,jpj) :: zCdU_u, zCdU_v ! top/bottom stress at u- & v-points
162	REAL(wp), DIMENSION(jpi,jpj) :: zhU, zhV ! fluxes
163	!
164	REAL(wp) :: zwdramp ! local scalar - only used if ln_wd_dl = .True.
165
166	INTEGER :: iwdg, jwdg, kwdg ! short-hand values for the indices of the output point
167
168	REAL(wp) :: zepsilon, zgamma ! - -
169	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zcpx, zcpy ! Wetting/Dying gravity filter coef.
170	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: ztwdmask, zuwdmask, zvwdmask ! ROMS wetting and drying masks at t,u,v points
171	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zuwdav2, zvwdav2 ! averages over the sub-steps of zuwdmask and zvwdmask
172	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zspgtrdu, zspgtrdv, zpvotrdu, zpvotrdv ! SPG and PVO trends (if l_trddyn)
173	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: ztautrdu, ztautrdv, zbfrtrdu, zbfrtrdv ! TAU and BFR trends (if l_trddyn)
174	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: ztfrtrdu, ztfrtrdv, ztottrdu, ztottrdv ! TFR and TOT trends (if l_trddyn)
175	!!----------------------------------------------------------------------
176	!
177	IF( ln_wd_il ) ALLOCATE( zcpx(jpi,jpj), zcpy(jpi,jpj) )
178	! !* Allocate temporary arrays
179	IF( ln_wd_dl ) ALLOCATE( ztwdmask(jpi,jpj), zuwdmask(jpi,jpj), zvwdmask(jpi,jpj), zuwdav2(jpi,jpj), zvwdav2(jpi,jpj))
180	!
181	IF( l_trddyn ) THEN
182	ALLOCATE( zspgtrdu(jpi,jpj), zspgtrdv(jpi,jpj), zpvotrdu(jpi,jpj), zpvotrdv(jpi,jpj), &
183	& ztautrdu(jpi,jpj), ztautrdv(jpi,jpj), zbfrtrdu(jpi,jpj), zbfrtrdv(jpi,jpj), &
184	& ztottrdu(jpi,jpj), ztottrdv(jpi,jpj) )
185	zspgtrdu(:,:) = 0._wp
186	zspgtrdv(:,:) = 0._wp
187	zpvotrdu(:,:) = 0._wp
188	zpvotrdv(:,:) = 0._wp
189	ztautrdu(:,:) = 0._wp
190	ztautrdv(:,:) = 0._wp
191	zbfrtrdu(:,:) = 0._wp
192	zbfrtrdv(:,:) = 0._wp
193	ztottrdu(:,:) = 0._wp
194	ztottrdv(:,:) = 0._wp
195	IF( ln_isfcav.OR.ln_drgice_imp ) THEN ! top+bottom friction (ocean cavities)
196	ALLOCATE( ztfrtrdu(jpi,jpj), ztfrtrdv(jpi,jpj) )
197	ztfrtrdu(:,:) = 0._wp
198	ztfrtrdv(:,:) = 0._wp
199	ENDIF
200	ENDIF
201	!
202	zu_trd(:,:) = 0._wp
203	zv_trd(:,:) = 0._wp
204	zu_spg(:,:) = 0._wp
205	zv_spg(:,:) = 0._wp
206	!
207	zmdi=1.e+20 ! missing data indicator for masking
208	!
209	zwdramp = r_rn_wdmin1 ! simplest ramp
210	! zwdramp = 1._wp / (rn_wdmin2 - rn_wdmin1) ! more general ramp
211	! ! inverse of baroclinic time step
212	IF( kt == nit000 .AND. neuler == 0 ) THEN ; r1_2dt_b = 1._wp / ( rdt )
213	ELSE ; r1_2dt_b = 1._wp / ( 2._wp * rdt )
214	ENDIF
215	!
216	ll_init = ln_bt_av ! if no time averaging, then no specific restart
217	ll_fw_start = .FALSE.
218	! ! time offset in steps for bdy data update
219	IF( .NOT.ln_bt_fw ) THEN ; noffset = - nn_baro
220	ELSE ; noffset = 0
221	ENDIF
222	!
223	IF( kt == nit000 ) THEN !* initialisation
224	!
225	IF(lwp) WRITE(numout,*)
226	IF(lwp) WRITE(numout,*) 'dyn_spg_ts : surface pressure gradient trend'
227	IF(lwp) WRITE(numout,*) '~~~~~~~~~~ free surface with time splitting'
228	IF(lwp) WRITE(numout,*)
229	!
230	IF( neuler == 0 ) ll_init=.TRUE.
231	!
232	IF( ln_bt_fw .OR. neuler == 0 ) THEN
233	ll_fw_start =.TRUE.
234	noffset = 0
235	ELSE
236	ll_fw_start =.FALSE.
237	ENDIF
238	! ! Set averaging weights and cycle length:
239	CALL ts_wgt( ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2 )
240	!
241	ENDIF
242	!
243	! If forward start at previous time step, and centered integration,
244	! then update averaging weights:
245	IF (.NOT.ln_bt_fw .AND.( neuler==0 .AND. kt==nit000+1 ) ) THEN
246	ll_fw_start=.FALSE.
247	CALL ts_wgt( ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2 )
248	ENDIF
249	!
250
251	! -----------------------------------------------------------------------------
252	! Phase 1 : Coupling between general trend and barotropic estimates (1st step)
253	! -----------------------------------------------------------------------------
254	!
255	!
256	! != zu_frc = 1/H e3*d/dt(Ua) =! (Vertical mean of Ua, the 3D trends)
257	! ! --------------------------- !
258	zu_frc(:,:) = SUM( e3u_n(:,:,:) * ua(:,:,:) * umask(:,:,:) , DIM=3 ) * r1_hu_n(:,:)
259	zv_frc(:,:) = SUM( e3v_n(:,:,:) * va(:,:,:) * vmask(:,:,:) , DIM=3 ) * r1_hv_n(:,:)
260	!
261	!
262	! != Ua => baroclinic trend =! (remove its vertical mean)
263	DO jk = 1, jpkm1 ! ------------------------ !
264	ua(:,:,jk) = ( ua(:,:,jk) - zu_frc(:,:) ) * umask(:,:,jk)
265	va(:,:,jk) = ( va(:,:,jk) - zv_frc(:,:) ) * vmask(:,:,jk)
266	END DO
267
268	!!gm Question here when removing the Vertically integrated trends, we remove the vertically integrated NL trends on momentum....
269	!!gm Is it correct to do so ? I think so...
270
271	! != remove 2D Coriolis and pressure gradient trends =!
272	! ! ------------------------------------------------- !
273	!
274	IF( kt == nit000 .OR. .NOT. ln_linssh ) CALL dyn_cor_2D_init ! Set zwz, the barotropic Coriolis force coefficient
275	! ! recompute zwz = f/depth at every time step for (.NOT.ln_linssh) as the water colomn height changes
276	!
277	! !* 2D Coriolis trends
278	zhU(:,:) = un_b(:,:) * hu_n(:,:) * e2u(:,:) ! now fluxes
279	zhV(:,:) = vn_b(:,:) * hv_n(:,:) * e1v(:,:) ! NB: FULL domain : put a value in last row and column
280	!
281	CALL dyn_cor_2d( hu_n, hv_n, un_b, vn_b, zhU, zhV, & ! <<== in
282	& zu_trd, zv_trd ) ! ==>> out
283	!
284	IF( l_trddyn ) THEN
285	! send correction to baroclinic planetary vorticity trend to trd_dyn
286	CALL trd_dyn( zu_trd, zv_trd, jpdyn_pvo_corr, kt )
287	ENDIF
288	!
289	IF( .NOT.ln_linssh ) THEN !* surface pressure gradient (variable volume only)
290	!
291	IF( l_trddyn ) THEN
292	zspgtrdu(:,:) = zu_trd(:,:)
293	zspgtrdv(:,:) = zv_trd(:,:)
294	ENDIF
295	!
296	IF( ln_wd_il ) THEN ! W/D : limiter applied to spgspg
297	CALL wad_spg( sshn, zcpx, zcpy ) ! Calculating W/D gravity filters, zcpx and zcpy
298	DO jj = 2, jpjm1
299	DO ji = 2, jpim1 ! SPG with the application of W/D gravity filters
300	zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) &
301	& * r1_e1u(ji,jj) * zcpx(ji,jj) * wdrampu(ji,jj) !jth
302	zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji ,jj+1) - sshn(ji ,jj ) ) &
303	& * r1_e2v(ji,jj) * zcpy(ji,jj) * wdrampv(ji,jj) !jth
304	END DO
305	END DO
306	ELSE ! now suface pressure gradient
307	DO jj = 2, jpjm1
308	DO ji = fs_2, fs_jpim1 ! vector opt.
309	zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) * r1_e1u(ji,jj)
310	zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji ,jj+1) - sshn(ji ,jj ) ) * r1_e2v(ji,jj)
311	END DO
312	END DO
313	ENDIF
314	!
315	IF( l_trddyn ) THEN
316	zspgtrdu(:,:) = zu_trd(:,:) - zspgtrdu(:,:)
317	zspgtrdv(:,:) = zv_trd(:,:) - zspgtrdv(:,:)
318	! send correction to HPG trend to trd_dyn
319	CALL trd_dyn( zspgtrdu, zspgtrdv, jpdyn_hpg_corr, kt )
320	! reset temporary arrays for use later
321	zspgtrdu(:,:) = 0._wp
322	zspgtrdv(:,:) = 0._wp
323	ENDIF
324	!
325	ENDIF
326	!
327	DO jj = 2, jpjm1 ! Remove coriolis term (and possibly spg) from barotropic trend
328	DO ji = fs_2, fs_jpim1
329	zu_frc(ji,jj) = zu_frc(ji,jj) - zu_trd(ji,jj) * ssumask(ji,jj)
330	zv_frc(ji,jj) = zv_frc(ji,jj) - zv_trd(ji,jj) * ssvmask(ji,jj)
331	END DO
332	END DO
333	!
334	! != Add bottom stress contribution from baroclinic velocities =!
335	! ! ----------------------------------------------------------- !
336	IF( l_trddyn ) THEN
337	!
338	! Output constant forcing terms (excluding top and bottom stresses) as diagnostics.
339	CALL trd_dyn( zu_frc, zv_frc, jpdyn_frc2d, kt )
340	!
341	CALL dyn_drg_init( zu_frc, zv_frc, zCdU_u, zCdU_v, & ! also provide the barotropic drag coefficients
342	& ztfrtrdu, ztfrtrdv, zbfrtrdu, zbfrtrdv )
343	!
344	ELSE
345	!
346	CALL dyn_drg_init( zu_frc, zv_frc, zCdU_u, zCdU_v ) ! also provide the barotropic drag coefficients
347	!
348	ENDIF
349	!
350	! != Add atmospheric pressure forcing =!
351	! ! ---------------------------------- !
352	IF( ln_apr_dyn ) THEN
353	IF( ln_bt_fw ) THEN ! FORWARD integration: use kt+1/2 pressure (NOW+1/2)
354	DO jj = 2, jpjm1
355	DO ji = fs_2, fs_jpim1 ! vector opt.
356	zu_frc(ji,jj) = zu_frc(ji,jj) + grav * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) ) * r1_e1u(ji,jj)
357	zv_frc(ji,jj) = zv_frc(ji,jj) + grav * ( ssh_ib (ji ,jj+1) - ssh_ib (ji,jj) ) * r1_e2v(ji,jj)
358	END DO
359	END DO
360	ELSE ! CENTRED integration: use kt-1/2 + kt+1/2 pressure (NOW)
361	zztmp = grav * r1_2
362	DO jj = 2, jpjm1
363	DO ji = fs_2, fs_jpim1 ! vector opt.
364	zu_frc(ji,jj) = zu_frc(ji,jj) + zztmp * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) &
365	& + ssh_ibb(ji+1,jj ) - ssh_ibb(ji,jj) ) * r1_e1u(ji,jj)
366	zv_frc(ji,jj) = zv_frc(ji,jj) + zztmp * ( ssh_ib (ji ,jj+1) - ssh_ib (ji,jj) &
367	& + ssh_ibb(ji ,jj+1) - ssh_ibb(ji,jj) ) * r1_e2v(ji,jj)
368	END DO
369	END DO
370	ENDIF
371	ENDIF
372	!
373	! != Add atmospheric pressure forcing =!
374	! ! ---------------------------------- !
375	IF( l_trddyn ) THEN
376	IF( nn_ice == 2 ) THEN ! Calculate and output the (partial) ice-ocean stress if using SI3.
377	ztautrdu(:,:) = 0._wp ; ztautrdv(:,:) = 0._wp
378	IF( ln_bt_fw ) THEN
379	DO jj = 2, jpjm1
380	DO ji = fs_2, fs_jpim1 ! vector opt.
381	ztautrdu(ji,jj) = r1_rau0 * uiceoc(ji,jj) * r1_hu_n(ji,jj)
382	ztautrdv(ji,jj) = r1_rau0 * viceoc(ji,jj) * r1_hv_n(ji,jj)
383	END DO
384	END DO
385	ELSE
386	zztmp = r1_rau0 * r1_2
387	DO jj = 2, jpjm1
388	DO ji = fs_2, fs_jpim1 ! vector opt.
389	ztautrdu(ji,jj) = zztmp * ( uiceoc_b(ji,jj) + uiceoc(ji,jj) ) * r1_hu_n(ji,jj)
390	ztautrdv(ji,jj) = zztmp * ( viceoc_b(ji,jj) + viceoc(ji,jj) ) * r1_hv_n(ji,jj)
391	END DO
392	END DO
393	ENDIF
394	IF( ln_isfcav.OR.ln_drgice_imp ) THEN
395	! Save this part of the ice-ocean drag as the first installment of top friction
396	CALL trd_dyn( ztautrdu, ztautrdv, jpdyn_iceoc2d, kt )
397	ELSE
398	! In this case this is the whole top friction trend
399	CALL trd_dyn( ztautrdu, ztautrdv, jpdyn_tfr, kt )
400	ENDIF
401	ENDIF
402	! initialise fields for wind stress trends
403	ztautrdu(:,:) = zu_frc(:,:)
404	ztautrdv(:,:) = zv_frc(:,:)
405	ENDIF
406	!
407	IF( ln_bt_fw ) THEN ! Add wind forcing
408	DO jj = 2, jpjm1
409	DO ji = fs_2, fs_jpim1 ! vector opt.
410	zu_frc(ji,jj) = zu_frc(ji,jj) + r1_rau0 * utau(ji,jj) * r1_hu_n(ji,jj)
411	zv_frc(ji,jj) = zv_frc(ji,jj) + r1_rau0 * vtau(ji,jj) * r1_hv_n(ji,jj)
412	END DO
413	END DO
414	ELSE
415	zztmp = r1_rau0 * r1_2
416	DO jj = 2, jpjm1
417	DO ji = fs_2, fs_jpim1 ! vector opt.
418	zu_frc(ji,jj) = zu_frc(ji,jj) + zztmp * ( utau_b(ji,jj) + utau(ji,jj) ) * r1_hu_n(ji,jj)
419	zv_frc(ji,jj) = zv_frc(ji,jj) + zztmp * ( vtau_b(ji,jj) + vtau(ji,jj) ) * r1_hv_n(ji,jj)
420	END DO
421	END DO
422	ENDIF
423	!
424	IF( l_trddyn ) THEN
425	! wind stress trend diagnostic
426	ztautrdu(:,:) = zu_frc(:,:) - ztautrdu(:,:)
427	ztautrdv(:,:) = zv_frc(:,:) - ztautrdv(:,:)
428	CALL trd_dyn( ztautrdu, ztautrdv, jpdyn_tau2d, kt )
429	ENDIF
430	! !----------------!
431	! !== sssh_frc ==! Right-Hand-Side of the barotropic ssh equation (over the FULL domain)
432	! !----------------!
433	! != Net water flux forcing applied to a water column =!
434	! ! --------------------------------------------------- !
435	IF (ln_bt_fw) THEN ! FORWARD integration: use kt+1/2 fluxes (NOW+1/2)
436	zssh_frc(:,:) = r1_rau0 * ( emp(:,:) - rnf(:,:) + fwfisf(:,:) )
437	ELSE ! CENTRED integration: use kt-1/2 + kt+1/2 fluxes (NOW)
438	zztmp = r1_rau0 * r1_2
439	zssh_frc(:,:) = zztmp * ( emp(:,:) + emp_b(:,:) - rnf(:,:) - rnf_b(:,:) + fwfisf(:,:) + fwfisf_b(:,:) )
440	ENDIF
441	! != Add Stokes drift divergence =! (if exist)
442	IF( ln_sdw ) THEN ! ----------------------------- !
443	zssh_frc(:,:) = zssh_frc(:,:) + div_sd(:,:)
444	ENDIF
445	!
446	#if defined key_asminc
447	! != Add the IAU weighted SSH increment =!
448	! ! ------------------------------------ !
449	IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN
450	zssh_frc(:,:) = zssh_frc(:,:) - ssh_iau(:,:)
451	ENDIF
452	#endif
453	! != Fill boundary data arrays for AGRIF
454	! ! ------------------------------------
455	#if defined key_agrif
456	IF( .NOT.Agrif_Root() ) CALL agrif_dta_ts( kt )
457	#endif
458	!
459	! -----------------------------------------------------------------------
460	! Phase 2 : Integration of the barotropic equations
461	! -----------------------------------------------------------------------
462	!
463	! ! ==================== !
464	! ! Initialisations !
465	! ! ==================== !
466	! Initialize barotropic variables:
467	IF( ll_init )THEN
468	sshbb_e(:,:) = 0._wp
469	ubb_e (:,:) = 0._wp
470	vbb_e (:,:) = 0._wp
471	sshb_e (:,:) = 0._wp
472	ub_e (:,:) = 0._wp
473	vb_e (:,:) = 0._wp
474	ENDIF
475	!
476	IF( ln_linssh ) THEN ! mid-step ocean depth is fixed (hup2_e=hu_n=hu_0)
477	zhup2_e(:,:) = hu_n(:,:)
478	zhvp2_e(:,:) = hv_n(:,:)
479	zhtp2_e(:,:) = ht_n(:,:)
480	ENDIF
481	!
482	IF (ln_bt_fw) THEN ! FORWARD integration: start from NOW fields
483	sshn_e(:,:) = sshn(:,:)
484	un_e (:,:) = un_b(:,:)
485	vn_e (:,:) = vn_b(:,:)
486	!
487	hu_e (:,:) = hu_n(:,:)
488	hv_e (:,:) = hv_n(:,:)
489	hur_e (:,:) = r1_hu_n(:,:)
490	hvr_e (:,:) = r1_hv_n(:,:)
491	ELSE ! CENTRED integration: start from BEFORE fields
492	sshn_e(:,:) = sshb(:,:)
493	un_e (:,:) = ub_b(:,:)
494	vn_e (:,:) = vb_b(:,:)
495	!
496	hu_e (:,:) = hu_b(:,:)
497	hv_e (:,:) = hv_b(:,:)
498	hur_e (:,:) = r1_hu_b(:,:)
499	hvr_e (:,:) = r1_hv_b(:,:)
500	ENDIF
501	!
502	! Initialize sums:
503	ua_b (:,:) = 0._wp ! After barotropic velocities (or transport if flux form)
504	va_b (:,:) = 0._wp
505	ssha (:,:) = 0._wp ! Sum for after averaged sea level
506	un_adv(:,:) = 0._wp ! Sum for now transport issued from ts loop
507	vn_adv(:,:) = 0._wp
508	!
509	IF( ln_wd_dl ) THEN
510	zuwdmask(:,:) = 0._wp ! set to zero for definiteness (not sure this is necessary)
511	zvwdmask(:,:) = 0._wp !
512	zuwdav2 (:,:) = 0._wp
513	zvwdav2 (:,:) = 0._wp
514	END IF
515
516	! ! ==================== !
517	DO jn = 1, icycle ! sub-time-step loop !
518	! ! ==================== !
519	!
520	l_full_nf_update = jn == icycle ! false: disable full North fold update (performances) for jn = 1 to icycle-1
521	!
522	! !== Update the forcing ==! (BDY and tides)
523	!
524	IF( ln_bdy .AND. ln_tide ) CALL bdy_dta_tides( kt, kit=jn, kt_offset= noffset+1 )
525	IF( ln_tide_pot .AND. ln_tide ) CALL upd_tide ( kt, kit=jn, kt_offset= noffset )
526	!
527	! !== extrapolation at mid-step ==! (jn+1/2)
528	!
529	! !* Set extrapolation coefficients for predictor step:
530	IF ((jn<3).AND.ll_init) THEN ! Forward
531	za1 = 1._wp
532	za2 = 0._wp
533	za3 = 0._wp
534	ELSE ! AB3-AM4 Coefficients: bet=0.281105
535	za1 = 1.781105_wp ! za1 = 3/2 + bet
536	za2 = -1.06221_wp ! za2 = -(1/2 + 2*bet)
537	za3 = 0.281105_wp ! za3 = bet
538	ENDIF
539	!
540	! !* Extrapolate barotropic velocities at mid-step (jn+1/2)
541	!-- m+1/2 m m-1 m-2 --!
542	!-- u = (3/2+beta) u -(1/2+2beta) u + beta u --!
543	!-------------------------------------------------------------------------!
544	ua_e(:,:) = za1 * un_e(:,:) + za2 * ub_e(:,:) + za3 * ubb_e(:,:)
545	va_e(:,:) = za1 * vn_e(:,:) + za2 * vb_e(:,:) + za3 * vbb_e(:,:)
546
547	IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only)
548	! ! ------------------
549	! Extrapolate Sea Level at step jit+0.5:
550	!-- m+1/2 m m-1 m-2 --!
551	!-- ssh = (3/2+beta) ssh -(1/2+2beta) ssh + beta ssh --!
552	!--------------------------------------------------------------------------------!
553	zsshp2_e(:,:) = za1 * sshn_e(:,:) + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:)
554
555	! set wetting & drying mask at tracer points for this barotropic mid-step
556	IF( ln_wd_dl ) CALL wad_tmsk( zsshp2_e, ztwdmask )
557	!
558	! ! ocean t-depth at mid-step
559	zhtp2_e(:,:) = ht_0(:,:) + zsshp2_e(:,:)
560	!
561	! ! ocean u- and v-depth at mid-step (separate DO-loops remove the need of a lbc_lnk)
562	DO jj = 1, jpj
563	DO ji = 1, jpim1 ! not jpi-column
564	zhup2_e(ji,jj) = hu_0(ji,jj) + r1_2 * r1_e1e2u(ji,jj) &
565	& * ( e1e2t(ji ,jj) * zsshp2_e(ji ,jj) &
566	& + e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj) ) * ssumask(ji,jj)
567	END DO
568	END DO
569	DO jj = 1, jpjm1 ! not jpj-row
570	DO ji = 1, jpi
571	zhvp2_e(ji,jj) = hv_0(ji,jj) + r1_2 * r1_e1e2v(ji,jj) &
572	& * ( e1e2t(ji,jj ) * zsshp2_e(ji,jj ) &
573	& + e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1) ) * ssvmask(ji,jj)
574	END DO
575	END DO
576	!
577	ENDIF
578	!
579	! !== after SSH ==! (jn+1)
580	!
581	! ! update (ua_e,va_e) to enforce volume conservation at open boundaries
582	! ! values of zhup2_e and zhvp2_e on the halo are not needed in bdy_vol2d
583	IF( ln_bdy .AND. ln_vol ) CALL bdy_vol2d( kt, jn, ua_e, va_e, zhup2_e, zhvp2_e )
584	!
585	! ! resulting flux at mid-step (not over the full domain)
586	zhU(1:jpim1,1:jpj ) = e2u(1:jpim1,1:jpj ) * ua_e(1:jpim1,1:jpj ) * zhup2_e(1:jpim1,1:jpj ) ! not jpi-column
587	zhV(1:jpi ,1:jpjm1) = e1v(1:jpi ,1:jpjm1) * va_e(1:jpi ,1:jpjm1) * zhvp2_e(1:jpi ,1:jpjm1) ! not jpj-row
588	!
589	#if defined key_agrif
590	! Set fluxes during predictor step to ensure volume conservation
591	IF( .NOT.Agrif_Root() .AND. ln_bt_fw ) THEN
592	IF( l_Westedge ) THEN
593	DO jj = 1, jpj
594	zhU(2:nbghostcells+1,jj) = ubdy_w(1:nbghostcells,jj) * e2u(2:nbghostcells+1,jj)
595	zhV(2:nbghostcells+1,jj) = vbdy_w(1:nbghostcells,jj) * e1v(2:nbghostcells+1,jj)
596	END DO
597	ENDIF
598	IF( l_Eastedge ) THEN
599	DO jj=1,jpj
600	zhU(nlci-nbghostcells-1:nlci-2,jj) = ubdy_e(1:nbghostcells,jj) * e2u(nlci-nbghostcells-1:nlci-2,jj)
601	zhV(nlci-nbghostcells :nlci-1,jj) = vbdy_e(1:nbghostcells,jj) * e1v(nlci-nbghostcells :nlci-1,jj)
602	END DO
603	ENDIF
604	IF( l_Southedge ) THEN
605	DO ji=1,jpi
606	zhV(ji,2:nbghostcells+1) = vbdy_s(ji,1:nbghostcells) * e1v(ji,2:nbghostcells+1)
607	zhU(ji,2:nbghostcells+1) = ubdy_s(ji,1:nbghostcells) * e2u(ji,2:nbghostcells+1)
608	END DO
609	ENDIF
610	IF( l_Northedge ) THEN
611	DO ji=1,jpi
612	zhV(ji,nlcj-nbghostcells-1:nlcj-2) = vbdy_n(ji,1:nbghostcells) * e1v(ji,nlcj-nbghostcells-1:nlcj-2)
613	zhU(ji,nlcj-nbghostcells :nlcj-1) = ubdy_n(ji,1:nbghostcells) * e2u(ji,nlcj-nbghostcells :nlcj-1)
614	END DO
615	ENDIF
616	ENDIF
617	#endif
618	IF( ln_wd_il ) CALL wad_lmt_bt(zhU, zhV, sshn_e, zssh_frc, rdtbt) !!gm wad_lmt_bt use of lbc_lnk on zhU, zhV
619
620	IF( ln_wd_dl ) THEN ! un_e and vn_e are set to zero at faces where
621	! ! the direction of the flow is from dry cells
622	CALL wad_Umsk( ztwdmask, zhU, zhV, un_e, vn_e, zuwdmask, zvwdmask ) ! not jpi colomn for U, not jpj row for V
623	!
624	ENDIF
625	!
626	!
627	! Compute Sea Level at step jit+1
628	!-- m+1 m m+1/2 --!
629	!-- ssh = ssh - delta_t' * [ frc + div( flux ) ] --!
630	!-------------------------------------------------------------------------!
631	DO jj = 2, jpjm1 ! INNER domain
632	DO ji = 2, jpim1
633	zhdiv = ( zhU(ji,jj) - zhU(ji-1,jj) + zhV(ji,jj) - zhV(ji,jj-1) ) * r1_e1e2t(ji,jj)
634	ssha_e(ji,jj) = ( sshn_e(ji,jj) - rdtbt * ( zssh_frc(ji,jj) + zhdiv ) ) * ssmask(ji,jj)
635	END DO
636	END DO
637	!
638	CALL lbc_lnk_multi( 'dynspg_ts', ssha_e, 'T', 1._wp, zhU, 'U', -1._wp, zhV, 'V', -1._wp )
639	!
640	! ! Sum over sub-time-steps to compute advective velocities
641	za2 = wgtbtp2(jn) ! zhU, zhV hold fluxes extrapolated at jn+0.5
642	un_adv(:,:) = un_adv(:,:) + za2 * zhU(:,:) * r1_e2u(:,:)
643	vn_adv(:,:) = vn_adv(:,:) + za2 * zhV(:,:) * r1_e1v(:,:)
644	! sum over sub-time-steps to decide which baroclinic velocities to set to zero (zuwdav2 is only used when ln_wd_dl_bc=True)
645	IF ( ln_wd_dl_bc ) THEN
646	zuwdav2(1:jpim1,1:jpj ) = zuwdav2(1:jpim1,1:jpj ) + za2 * zuwdmask(1:jpim1,1:jpj ) ! not jpi-column
647	zvwdav2(1:jpi ,1:jpjm1) = zvwdav2(1:jpi ,1:jpjm1) + za2 * zvwdmask(1:jpi ,1:jpjm1) ! not jpj-row
648	END IF
649	!
650	! Duplicate sea level across open boundaries (this is only cosmetic if linssh=T)
651	IF( ln_bdy ) CALL bdy_ssh( ssha_e )
652	#if defined key_agrif
653	IF( .NOT.Agrif_Root() ) CALL agrif_ssh_ts( jn )
654	#endif
655	!
656	! Sea Surface Height at u-,v-points (vvl case only)
657	IF( .NOT.ln_linssh ) THEN
658	DO jj = 2, jpjm1 ! INNER domain, will be extended to whole domain later
659	DO ji = 2, jpim1 ! NO Vector Opt.
660	zsshu_a(ji,jj) = r1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj) &
661	& * ( e1e2t(ji ,jj ) * ssha_e(ji ,jj ) &
662	& + e1e2t(ji+1,jj ) * ssha_e(ji+1,jj ) )
663	zsshv_a(ji,jj) = r1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj) &
664	& * ( e1e2t(ji ,jj ) * ssha_e(ji ,jj ) &
665	& + e1e2t(ji ,jj+1) * ssha_e(ji ,jj+1) )
666	END DO
667	END DO
668	ENDIF
669	!
670	! Half-step back interpolation of SSH for surface pressure computation at step jit+1/2
671	!-- m+1/2 m+1 m m-1 m-2 --!
672	!-- ssh' = za0 * ssh + za1 * ssh + za2 * ssh + za3 * ssh --!
673	!------------------------------------------------------------------------------------------!
674	CALL ts_bck_interp( jn, ll_init, za0, za1, za2, za3 ) ! coeficients of the interpolation
675	zsshp2_e(:,:) = za0 * ssha_e(:,:) + za1 * sshn_e (:,:) &
676	& + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:)
677	!
678	! ! Surface pressure gradient
679	zldg = ( 1._wp - rn_scal_load ) * grav ! local factor
680	DO jj = 2, jpjm1
681	DO ji = 2, jpim1
682	zu_spg(ji,jj) = - zldg * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj)
683	zv_spg(ji,jj) = - zldg * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj)
684	END DO
685	END DO
686	IF( ln_wd_il ) THEN ! W/D : gravity filters applied on pressure gradient
687	CALL wad_spg( zsshp2_e, zcpx, zcpy ) ! Calculating W/D gravity filters
688	zu_spg(2:jpim1,2:jpjm1) = zu_spg(2:jpim1,2:jpjm1) * zcpx(2:jpim1,2:jpjm1)
689	zv_spg(2:jpim1,2:jpjm1) = zv_spg(2:jpim1,2:jpjm1) * zcpy(2:jpim1,2:jpjm1)
690	ENDIF
691	!
692	IF( l_trddyn ) THEN
693	za2 = wgtbtp2(jn)
694	zspgtrdu(:,:) = zspgtrdu(:,:) + za2 * zu_spg(:,:) * ssumask(:,:)
695	zspgtrdv(:,:) = zspgtrdv(:,:) + za2 * zv_spg(:,:) * ssvmask(:,:)
696	ENDIF
697	!
698	! Add Coriolis trend:
699	! zwz array below or triads normally depend on sea level with ln_linssh=F and should be updated
700	! at each time step. We however keep them constant here for optimization.
701	! Recall that zhU and zhV hold fluxes at jn+0.5 (extrapolated not backward interpolated)
702	CALL dyn_cor_2d( zhup2_e, zhvp2_e, ua_e, va_e, zhU, zhV, zu_trd, zv_trd )
703	!
704	IF( l_trddyn ) THEN
705	za2 = wgtbtp2(jn)
706	zpvotrdu(:,:) = zpvotrdu(:,:) + za2 * zu_trd(:,:) * ssumask(:,:)
707	zpvotrdv(:,:) = zpvotrdv(:,:) + za2 * zv_trd(:,:) * ssvmask(:,:)
708	ENDIF
709	!
710	! Add tidal astronomical forcing if defined
711	IF ( ln_tide .AND. ln_tide_pot ) THEN
712	DO jj = 2, jpjm1
713	DO ji = fs_2, fs_jpim1 ! vector opt.
714	zu_trd(ji,jj) = zu_trd(ji,jj) + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) * r1_e1u(ji,jj)
715	zv_trd(ji,jj) = zv_trd(ji,jj) + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) * r1_e2v(ji,jj)
716	END DO
717	END DO
718	ENDIF
719	!
720	! Add bottom stresses:
721	!jth do implicitly instead
722	IF ( .NOT. ll_wd ) THEN ! Revert to explicit for bit comparison tests in non wad runs
723	DO jj = 2, jpjm1
724	DO ji = fs_2, fs_jpim1 ! vector opt.
725	zu_trd(ji,jj) = zu_trd(ji,jj) + zCdU_u(ji,jj) * un_e(ji,jj) * hur_e(ji,jj)
726	zv_trd(ji,jj) = zv_trd(ji,jj) + zCdU_v(ji,jj) * vn_e(ji,jj) * hvr_e(ji,jj)
727	END DO
728	END DO
729	IF( l_trddyn ) THEN
730	za2 = wgtbtp2(jn)
731	IF( ln_isfcav.OR.ln_drgice_imp ) THEN ! top+bottom friction (ocean cavities)
732	DO jj = 2, jpjm1
733	DO ji = fs_2, fs_jpim1 ! vector opt.
734	ztfrtrdu(ji,jj) = ztfrtrdu(ji,jj) + za2 * 0.5_wp( rCdU_top(ji+1,jj)+rCdU_top(ji,jj)) un_e(ji,jj) * hur_e(ji,jj)
735	ztfrtrdv(ji,jj) = ztfrtrdv(ji,jj) + za2 * 0.5_wp( rCdU_top(ji+1,jj)+rCdU_top(ji,jj)) vn_e(ji,jj) * hvr_e(ji,jj)
736	END DO
737	END DO
738	ENDIF
739	DO jj = 2, jpjm1
740	DO ji = fs_2, fs_jpim1 ! vector opt.
741	zbfrtrdu(ji,jj) = zbfrtrdu(ji,jj) + za2 * 0.5_wp( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj)) un_e(ji,jj) * hur_e(ji,jj)
742	zbfrtrdv(ji,jj) = zbfrtrdv(ji,jj) + za2 * 0.5_wp( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj)) vn_e(ji,jj) * hvr_e(ji,jj)
743	END DO
744	END DO
745	ENDIF
746	ENDIF
747	!
748	! Set next velocities:
749	! Compute barotropic speeds at step jit+1 (h : total height of the water colomn)
750	!-- VECTOR FORM
751	!-- m+1 m / m+1/2 \ --!
752	!-- u = u + delta_t' * \ (1-r)g grad_x( ssh') - f * k vect u + frc / --!
753	!-- --!
754	!-- FLUX FORM --!
755	!-- m+1 __1__ / m m / m+1/2 m+1/2 m+1/2 n \ \ --!
756	!-- u = m+1 \| h * u + delta_t' * \ h * (1-r)g grad_x( ssh') - h * f * k vect u + h * frc / \| --!
757	!-- h \ / --!
758	!------------------------------------------------------------------------------------------------------------------------!
759	IF( ln_dynadv_vec .OR. ln_linssh ) THEN !* Vector form
760	DO jj = 2, jpjm1
761	DO ji = fs_2, fs_jpim1 ! vector opt.
762	ua_e(ji,jj) = ( un_e(ji,jj) &
763	& + rdtbt * ( zu_spg(ji,jj) &
764	& + zu_trd(ji,jj) &
765	& + zu_frc(ji,jj) ) &
766	& ) * ssumask(ji,jj)
767
768	va_e(ji,jj) = ( vn_e(ji,jj) &
769	& + rdtbt * ( zv_spg(ji,jj) &
770	& + zv_trd(ji,jj) &
771	& + zv_frc(ji,jj) ) &
772	& ) * ssvmask(ji,jj)
773	END DO
774	END DO
775	!
776	ELSE !* Flux form
777	DO jj = 2, jpjm1
778	DO ji = 2, jpim1
779	! ! hu_e, hv_e hold depth at jn, zhup2_e, zhvp2_e hold extrapolated depth at jn+1/2
780	! ! backward interpolated depth used in spg terms at jn+1/2
781	zhu_bck = hu_0(ji,jj) + r1_2r1_e1e2u(ji,jj) ( e1e2t(ji ,jj) * zsshp2_e(ji ,jj) &
782	& + e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj) ) * ssumask(ji,jj)
783	zhv_bck = hv_0(ji,jj) + r1_2r1_e1e2v(ji,jj) ( e1e2t(ji,jj ) * zsshp2_e(ji,jj ) &
784	& + e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1) ) * ssvmask(ji,jj)
785	! ! inverse depth at jn+1
786	z1_hu = ssumask(ji,jj) / ( hu_0(ji,jj) + zsshu_a(ji,jj) + 1._wp - ssumask(ji,jj) )
787	z1_hv = ssvmask(ji,jj) / ( hv_0(ji,jj) + zsshv_a(ji,jj) + 1._wp - ssvmask(ji,jj) )
788	!
789	ua_e(ji,jj) = ( hu_e (ji,jj) * un_e (ji,jj) &
790	& + rdtbt * ( zhu_bck * zu_spg (ji,jj) & !
791	& + zhup2_e(ji,jj) * zu_trd (ji,jj) & !
792	& + hu_n (ji,jj) * zu_frc (ji,jj) ) ) * z1_hu
793	!
794	va_e(ji,jj) = ( hv_e (ji,jj) * vn_e (ji,jj) &
795	& + rdtbt * ( zhv_bck * zv_spg (ji,jj) & !
796	& + zhvp2_e(ji,jj) * zv_trd (ji,jj) & !
797	& + hv_n (ji,jj) * zv_frc (ji,jj) ) ) * z1_hv
798	END DO
799	END DO
800	ENDIF
801	!jth implicit bottom friction:
802	IF ( ll_wd ) THEN ! revert to explicit for bit comparison tests in non wad runs
803	DO jj = 2, jpjm1
804	DO ji = fs_2, fs_jpim1 ! vector opt.
805	ua_e(ji,jj) = ua_e(ji,jj) /(1.0 - rdtbt * zCdU_u(ji,jj) * hur_e(ji,jj))
806	va_e(ji,jj) = va_e(ji,jj) /(1.0 - rdtbt * zCdU_v(ji,jj) * hvr_e(ji,jj))
807	END DO
808	END DO
809	ENDIF
810
811	IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only)
812	hu_e (2:jpim1,2:jpjm1) = hu_0(2:jpim1,2:jpjm1) + zsshu_a(2:jpim1,2:jpjm1)
813	hur_e(2:jpim1,2:jpjm1) = ssumask(2:jpim1,2:jpjm1) / ( hu_e(2:jpim1,2:jpjm1) + 1._wp - ssumask(2:jpim1,2:jpjm1) )
814	hv_e (2:jpim1,2:jpjm1) = hv_0(2:jpim1,2:jpjm1) + zsshv_a(2:jpim1,2:jpjm1)
815	hvr_e(2:jpim1,2:jpjm1) = ssvmask(2:jpim1,2:jpjm1) / ( hv_e(2:jpim1,2:jpjm1) + 1._wp - ssvmask(2:jpim1,2:jpjm1) )
816	CALL lbc_lnk_multi( 'dynspg_ts', ua_e , 'U', -1._wp, va_e , 'V', -1._wp &
817	& , hu_e , 'U', 1._wp, hv_e , 'V', 1._wp &
818	& , hur_e, 'U', 1._wp, hvr_e, 'V', 1._wp )
819	ELSE
820	CALL lbc_lnk_multi( 'dynspg_ts', ua_e , 'U', -1._wp, va_e , 'V', -1._wp )
821	ENDIF
822	!
823	!
824	! ! open boundaries
825	IF( ln_bdy ) CALL bdy_dyn2d( jn, ua_e, va_e, un_e, vn_e, hur_e, hvr_e, ssha_e )
826	#if defined key_agrif
827	IF( .NOT.Agrif_Root() ) CALL agrif_dyn_ts( jn ) ! Agrif
828	#endif
829	! !* Swap
830	! ! ----
831	ubb_e (:,:) = ub_e (:,:)
832	ub_e (:,:) = un_e (:,:)
833	un_e (:,:) = ua_e (:,:)
834	!
835	vbb_e (:,:) = vb_e (:,:)
836	vb_e (:,:) = vn_e (:,:)
837	vn_e (:,:) = va_e (:,:)
838	!
839	sshbb_e(:,:) = sshb_e(:,:)
840	sshb_e (:,:) = sshn_e(:,:)
841	sshn_e (:,:) = ssha_e(:,:)
842
843	! !* Sum over whole bt loop
844	! ! ----------------------
845	za1 = wgtbtp1(jn)
846	IF( ln_dynadv_vec .OR. ln_linssh ) THEN ! Sum velocities
847	ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:)
848	va_b (:,:) = va_b (:,:) + za1 * va_e (:,:)
849	ELSE ! Sum transports
850	IF ( .NOT.ln_wd_dl ) THEN
851	ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:) * hu_e (:,:)
852	va_b (:,:) = va_b (:,:) + za1 * va_e (:,:) * hv_e (:,:)
853	ELSE
854	ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:) * hu_e (:,:) * zuwdmask(:,:)
855	va_b (:,:) = va_b (:,:) + za1 * va_e (:,:) * hv_e (:,:) * zvwdmask(:,:)
856	END IF
857	ENDIF
858	! ! Sum sea level
859	ssha(:,:) = ssha(:,:) + za1 * ssha_e(:,:)
860
861	! ! ==================== !
862	END DO ! end loop !
863	! ! ==================== !
864	! -----------------------------------------------------------------------------
865	! Phase 3. update the general trend with the barotropic trend
866	! -----------------------------------------------------------------------------
867	!
868	! Set advection velocity correction:
869	IF (ln_bt_fw) THEN
870	IF( .NOT.( kt == nit000 .AND. neuler==0 ) ) THEN
871	DO jj = 1, jpj
872	DO ji = 1, jpi
873	zun_save = un_adv(ji,jj)
874	zvn_save = vn_adv(ji,jj)
875	! ! apply the previously computed correction
876	un_adv(ji,jj) = r1_2 * ( ub2_b(ji,jj) + zun_save - atfp * un_bf(ji,jj) )
877	vn_adv(ji,jj) = r1_2 * ( vb2_b(ji,jj) + zvn_save - atfp * vn_bf(ji,jj) )
878	! ! Update corrective fluxes for next time step
879	un_bf(ji,jj) = atfp * un_bf(ji,jj) + ( zun_save - ub2_b(ji,jj) )
880	vn_bf(ji,jj) = atfp * vn_bf(ji,jj) + ( zvn_save - vb2_b(ji,jj) )
881	! ! Save integrated transport for next computation
882	ub2_b(ji,jj) = zun_save
883	vb2_b(ji,jj) = zvn_save
884	END DO
885	END DO
886	ELSE
887	un_bf(:,:) = 0._wp ! corrective fluxes for next time step set to zero
888	vn_bf(:,:) = 0._wp
889	ub2_b(:,:) = un_adv(:,:) ! Save integrated transport for next computation
890	vb2_b(:,:) = vn_adv(:,:)
891	END IF
892	ENDIF
893
894
895	!
896	! Update barotropic trend:
897	IF( ln_dynadv_vec .OR. ln_linssh ) THEN
898	DO jk=1,jpkm1
899	ua(:,:,jk) = ua(:,:,jk) + ( ua_b(:,:) - ub_b(:,:) ) * r1_2dt_b
900	va(:,:,jk) = va(:,:,jk) + ( va_b(:,:) - vb_b(:,:) ) * r1_2dt_b
901	END DO
902	IF( l_trddyn ) THEN
903	ztottrdu(:,:) = ( ua_b(:,:) - ub_b(:,:) ) * r1_2dt_b
904	ztottrdv(:,:) = ( va_b(:,:) - vb_b(:,:) ) * r1_2dt_b
905	CALL trd_dyn( ztottrdu, ztottrdv, jpdyn_tot, kt )
906	ENDIF
907	ELSE
908	! At this stage, ssha has been corrected: compute new depths at velocity points
909	DO jj = 1, jpjm1
910	DO ji = 1, jpim1 ! NO Vector Opt.
911	zsshu_a(ji,jj) = r1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj) &
912	& * ( e1e2t(ji ,jj) * ssha(ji ,jj) &
913	& + e1e2t(ji+1,jj) * ssha(ji+1,jj) )
914	zsshv_a(ji,jj) = r1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj) &
915	& * ( e1e2t(ji,jj ) * ssha(ji,jj ) &
916	& + e1e2t(ji,jj+1) * ssha(ji,jj+1) )
917	END DO
918	END DO
919	CALL lbc_lnk_multi( 'dynspg_ts', zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) ! Boundary conditions
920	!
921	DO jk=1,jpkm1
922	ua(:,:,jk) = ua(:,:,jk) + r1_hu_n(:,:) * ( ua_b(:,:) - ub_b(:,:) * hu_b(:,:) ) * r1_2dt_b
923	va(:,:,jk) = va(:,:,jk) + r1_hv_n(:,:) * ( va_b(:,:) - vb_b(:,:) * hv_b(:,:) ) * r1_2dt_b
924	END DO
925	! Save barotropic velocities not transport:
926	ua_b(:,:) = ua_b(:,:) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - ssumask(:,:) )
927	va_b(:,:) = va_b(:,:) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - ssvmask(:,:) )
928	IF( l_trddyn ) THEN
929	ztottrdu(:,:) = r1_hu_n(:,:) * ( ua_b(:,:) - ub_b(:,:) * hu_b(:,:) ) * r1_2dt_b
930	ztottrdv(:,:) = r1_hv_n(:,:) * ( va_b(:,:) - vb_b(:,:) * hv_b(:,:) ) * r1_2dt_b
931	CALL trd_dyn( ztottrdu, ztottrdv, jpdyn_tot, kt )
932	ENDIF
933	ENDIF
934
935
936	! Correct velocities so that the barotropic velocity equals (un_adv, vn_adv) (in all cases)
937	DO jk = 1, jpkm1
938	un(:,:,jk) = ( un(:,:,jk) + un_adv(:,:)r1_hu_n(:,:) - un_b(:,:) ) umask(:,:,jk)
939	vn(:,:,jk) = ( vn(:,:,jk) + vn_adv(:,:)r1_hv_n(:,:) - vn_b(:,:) ) vmask(:,:,jk)
940	END DO
941
942	IF ( ln_wd_dl .and. ln_wd_dl_bc) THEN
943	! need to set lbc here because not done prior time averaging
944	CALL lbc_lnk_multi( 'dynspg_ts', zuwdav2, 'U', 1._wp, zvwdav2, 'V', 1._wp)
945	DO jk = 1, jpkm1
946	un(:,:,jk) = ( un_adv(:,:)*r1_hu_n(:,:) &
947	& + zuwdav2(:,:)(un(:,:,jk) - un_adv(:,:)r1_hu_n(:,:)) ) * umask(:,:,jk)
948	vn(:,:,jk) = ( vn_adv(:,:)*r1_hv_n(:,:) &
949	& + zvwdav2(:,:)(vn(:,:,jk) - vn_adv(:,:)r1_hv_n(:,:)) ) * vmask(:,:,jk)
950	END DO
951	END IF
952
953
954	CALL iom_put( "ubar", un_adv(:,:)*r1_hu_n(:,:) ) ! barotropic i-current
955	CALL iom_put( "vbar", vn_adv(:,:)*r1_hv_n(:,:) ) ! barotropic i-current
956	!
957	#if defined key_agrif
958	! Save time integrated fluxes during child grid integration
959	! (used to update coarse grid transports at next time step)
960	!
961	IF( .NOT.Agrif_Root() .AND. ln_bt_fw ) THEN
962	IF( Agrif_NbStepint() == 0 ) THEN
963	ub2_i_b(:,:) = 0._wp
964	vb2_i_b(:,:) = 0._wp
965	END IF
966	!
967	za1 = 1._wp / REAL(Agrif_rhot(), wp)
968	ub2_i_b(:,:) = ub2_i_b(:,:) + za1 * ub2_b(:,:)
969	vb2_i_b(:,:) = vb2_i_b(:,:) + za1 * vb2_b(:,:)
970	ENDIF
971	#endif
972	! !* write time-spliting arrays in the restart
973	IF( lrst_oce .AND.ln_bt_fw ) CALL ts_rst( kt, 'WRITE' )
974	!
975	IF( ln_wd_il ) DEALLOCATE( zcpx, zcpy )
976	IF( ln_wd_dl ) DEALLOCATE( ztwdmask, zuwdmask, zvwdmask, zuwdav2, zvwdav2 )
977	!
978	IF( l_trddyn ) THEN
979	CALL trd_dyn( zspgtrdu, zspgtrdv, jpdyn_spg, kt )
980	CALL trd_dyn( zpvotrdu, zpvotrdv, jpdyn_pvo, kt )
981	CALL trd_dyn( zbfrtrdu, zbfrtrdv, jpdyn_bfr, kt )
982	DEALLOCATE( zspgtrdu, zspgtrdv, zpvotrdu, zpvotrdv, ztautrdu, ztautrdv, zbfrtrdu, zbfrtrdv )
983	IF( ln_isfcav.OR.ln_drgice_imp ) THEN ! top+bottom friction (ocean cavities)
984	CALL trd_dyn( ztfrtrdu, ztfrtrdv, jpdyn_tfr, kt )
985	DEALLOCATE( ztfrtrdu, ztfrtrdv )
986	ENDIF
987	ENDIF
988	!
989	CALL iom_put( "baro_u" , un_bssumask(:,:)+zmdi(1.-ssumask(:,:) ) ) ! Barotropic U Velocity
990	CALL iom_put( "baro_v" , vn_bssvmask(:,:)+zmdi(1.-ssvmask(:,:) ) ) ! Barotropic V Velocity
991	!
992	END SUBROUTINE dyn_spg_ts
993
994
995	SUBROUTINE ts_wgt( ll_av, ll_fw, jpit, zwgt1, zwgt2)
996	!!---------------------------------------------------------------------
997	!! * ROUTINE ts_wgt *
998	!!
999	!! ** Purpose : Set time-splitting weights for temporal averaging (or not)
1000	!!----------------------------------------------------------------------
1001	LOGICAL, INTENT(in) :: ll_av ! temporal averaging=.true.
1002	LOGICAL, INTENT(in) :: ll_fw ! forward time splitting =.true.
1003	INTEGER, INTENT(inout) :: jpit ! cycle length
1004	REAL(wp), DIMENSION(3*nn_baro), INTENT(inout) :: zwgt1, & ! Primary weights
1005	zwgt2 ! Secondary weights
1006
1007	INTEGER :: jic, jn, ji ! temporary integers
1008	REAL(wp) :: za1, za2
1009	!!----------------------------------------------------------------------
1010
1011	zwgt1(:) = 0._wp
1012	zwgt2(:) = 0._wp
1013
1014	! Set time index when averaged value is requested
1015	IF (ll_fw) THEN
1016	jic = nn_baro
1017	ELSE
1018	jic = 2 * nn_baro
1019	ENDIF
1020
1021	! Set primary weights:
1022	IF (ll_av) THEN
1023	! Define simple boxcar window for primary weights
1024	! (width = nn_baro, centered around jic)
1025	SELECT CASE ( nn_bt_flt )
1026	CASE( 0 ) ! No averaging
1027	zwgt1(jic) = 1._wp
1028	jpit = jic
1029
1030	CASE( 1 ) ! Boxcar, width = nn_baro
1031	DO jn = 1, 3*nn_baro
1032	za1 = ABS(float(jn-jic))/float(nn_baro)
1033	IF (za1 < 0.5_wp) THEN
1034	zwgt1(jn) = 1._wp
1035	jpit = jn
1036	ENDIF
1037	ENDDO
1038
1039	CASE( 2 ) ! Boxcar, width = 2 * nn_baro
1040	DO jn = 1, 3*nn_baro
1041	za1 = ABS(float(jn-jic))/float(nn_baro)
1042	IF (za1 < 1._wp) THEN
1043	zwgt1(jn) = 1._wp
1044	jpit = jn
1045	ENDIF
1046	ENDDO
1047	CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_bt_flt' )
1048	END SELECT
1049
1050	ELSE ! No time averaging
1051	zwgt1(jic) = 1._wp
1052	jpit = jic
1053	ENDIF
1054
1055	! Set secondary weights
1056	DO jn = 1, jpit
1057	DO ji = jn, jpit
1058	zwgt2(jn) = zwgt2(jn) + zwgt1(ji)
1059	END DO
1060	END DO
1061
1062	! Normalize weigths:
1063	za1 = 1._wp / SUM(zwgt1(1:jpit))
1064	za2 = 1._wp / SUM(zwgt2(1:jpit))
1065	DO jn = 1, jpit
1066	zwgt1(jn) = zwgt1(jn) * za1
1067	zwgt2(jn) = zwgt2(jn) * za2
1068	END DO
1069	!
1070	END SUBROUTINE ts_wgt
1071
1072
1073	SUBROUTINE ts_rst( kt, cdrw )
1074	!!---------------------------------------------------------------------
1075	!! * ROUTINE ts_rst *
1076	!!
1077	!! ** Purpose : Read or write time-splitting arrays in restart file
1078	!!----------------------------------------------------------------------
1079	INTEGER , INTENT(in) :: kt ! ocean time-step
1080	CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag
1081	!!----------------------------------------------------------------------
1082	!
1083	IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialise
1084	! ! ---------------
1085	IF( ln_rstart .AND. ln_bt_fw .AND. (neuler/=0) ) THEN !* Read the restart file
1086	CALL iom_get( numror, jpdom_autoglo, 'ub2_b' , ub2_b (:,:), ldxios = lrxios )
1087	CALL iom_get( numror, jpdom_autoglo, 'vb2_b' , vb2_b (:,:), ldxios = lrxios )
1088	CALL iom_get( numror, jpdom_autoglo, 'un_bf' , un_bf (:,:), ldxios = lrxios )
1089	CALL iom_get( numror, jpdom_autoglo, 'vn_bf' , vn_bf (:,:), ldxios = lrxios )
1090	IF( .NOT.ln_bt_av ) THEN
1091	CALL iom_get( numror, jpdom_autoglo, 'sshbb_e' , sshbb_e(:,:), ldxios = lrxios )
1092	CALL iom_get( numror, jpdom_autoglo, 'ubb_e' , ubb_e(:,:), ldxios = lrxios )
1093	CALL iom_get( numror, jpdom_autoglo, 'vbb_e' , vbb_e(:,:), ldxios = lrxios )
1094	CALL iom_get( numror, jpdom_autoglo, 'sshb_e' , sshb_e(:,:), ldxios = lrxios )
1095	CALL iom_get( numror, jpdom_autoglo, 'ub_e' , ub_e(:,:), ldxios = lrxios )
1096	CALL iom_get( numror, jpdom_autoglo, 'vb_e' , vb_e(:,:), ldxios = lrxios )
1097	ENDIF
1098	#if defined key_agrif
1099	! Read time integrated fluxes
1100	IF ( .NOT.Agrif_Root() ) THEN
1101	CALL iom_get( numror, jpdom_autoglo, 'ub2_i_b' , ub2_i_b(:,:), ldxios = lrxios )
1102	CALL iom_get( numror, jpdom_autoglo, 'vb2_i_b' , vb2_i_b(:,:), ldxios = lrxios )
1103	ENDIF
1104	#endif
1105	ELSE !* Start from rest
1106	IF(lwp) WRITE(numout,*)
1107	IF(lwp) WRITE(numout,*) ' ==>>> start from rest: set barotropic values to 0'
1108	ub2_b (:,:) = 0._wp ; vb2_b (:,:) = 0._wp ! used in the 1st interpol of agrif
1109	un_adv(:,:) = 0._wp ; vn_adv(:,:) = 0._wp ! used in the 1st interpol of agrif
1110	un_bf (:,:) = 0._wp ; vn_bf (:,:) = 0._wp ! used in the 1st update of agrif
1111	#if defined key_agrif
1112	IF ( .NOT.Agrif_Root() ) THEN
1113	ub2_i_b(:,:) = 0._wp ; vb2_i_b(:,:) = 0._wp ! used in the 1st update of agrif
1114	ENDIF
1115	#endif
1116	ENDIF
1117	!
1118	ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN ! Create restart file
1119	! ! -------------------
1120	IF(lwp) WRITE(numout,*) '---- ts_rst ----'
1121	IF( lwxios ) CALL iom_swap( cwxios_context )
1122	CALL iom_rstput( kt, nitrst, numrow, 'ub2_b' , ub2_b (:,:), ldxios = lwxios )
1123	CALL iom_rstput( kt, nitrst, numrow, 'vb2_b' , vb2_b (:,:), ldxios = lwxios )
1124	CALL iom_rstput( kt, nitrst, numrow, 'un_bf' , un_bf (:,:), ldxios = lwxios )
1125	CALL iom_rstput( kt, nitrst, numrow, 'vn_bf' , vn_bf (:,:), ldxios = lwxios )
1126	!
1127	IF (.NOT.ln_bt_av) THEN
1128	CALL iom_rstput( kt, nitrst, numrow, 'sshbb_e' , sshbb_e(:,:), ldxios = lwxios )
1129	CALL iom_rstput( kt, nitrst, numrow, 'ubb_e' , ubb_e(:,:), ldxios = lwxios )
1130	CALL iom_rstput( kt, nitrst, numrow, 'vbb_e' , vbb_e(:,:), ldxios = lwxios )
1131	CALL iom_rstput( kt, nitrst, numrow, 'sshb_e' , sshb_e(:,:), ldxios = lwxios )
1132	CALL iom_rstput( kt, nitrst, numrow, 'ub_e' , ub_e(:,:), ldxios = lwxios )
1133	CALL iom_rstput( kt, nitrst, numrow, 'vb_e' , vb_e(:,:), ldxios = lwxios )
1134	ENDIF
1135	#if defined key_agrif
1136	! Save time integrated fluxes
1137	IF ( .NOT.Agrif_Root() ) THEN
1138	CALL iom_rstput( kt, nitrst, numrow, 'ub2_i_b' , ub2_i_b(:,:), ldxios = lwxios )
1139	CALL iom_rstput( kt, nitrst, numrow, 'vb2_i_b' , vb2_i_b(:,:), ldxios = lwxios )
1140	ENDIF
1141	#endif
1142	IF( lwxios ) CALL iom_swap( cxios_context )
1143	ENDIF
1144	!
1145	END SUBROUTINE ts_rst
1146
1147
1148	SUBROUTINE dyn_spg_ts_init
1149	!!---------------------------------------------------------------------
1150	!! * ROUTINE dyn_spg_ts_init *
1151	!!
1152	!! ** Purpose : Set time splitting options
1153	!!----------------------------------------------------------------------
1154	INTEGER :: ji ,jj ! dummy loop indices
1155	REAL(wp) :: zxr2, zyr2, zcmax ! local scalar
1156	REAL(wp), DIMENSION(jpi,jpj) :: zcu
1157	INTEGER :: inum
1158	!!----------------------------------------------------------------------
1159	!
1160	! Max courant number for ext. grav. waves
1161	!
1162	DO jj = 1, jpj
1163	DO ji =1, jpi
1164	zxr2 = r1_e1t(ji,jj) * r1_e1t(ji,jj)
1165	zyr2 = r1_e2t(ji,jj) * r1_e2t(ji,jj)
1166	zcu(ji,jj) = SQRT( grav * MAX(ht_0(ji,jj),0._wp) * (zxr2 + zyr2) )
1167	END DO
1168	END DO
1169	!
1170	zcmax = MAXVAL( zcu(:,:) )
1171	CALL mpp_max( 'dynspg_ts', zcmax )
1172
1173	! Estimate number of iterations to satisfy a max courant number= rn_bt_cmax
1174	IF( ln_bt_auto ) nn_baro = CEILING( rdt / rn_bt_cmax * zcmax)
1175
1176	rdtbt = rdt / REAL( nn_baro , wp )
1177	zcmax = zcmax * rdtbt
1178	! Print results
1179	IF(lwp) WRITE(numout,*)
1180	IF(lwp) WRITE(numout,*) 'dyn_spg_ts_init : split-explicit free surface'
1181	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~'
1182	IF( ln_bt_auto ) THEN
1183	IF(lwp) WRITE(numout,*) ' ln_ts_auto =.true. Automatically set nn_baro '
1184	IF(lwp) WRITE(numout,*) ' Max. courant number allowed: ', rn_bt_cmax
1185	ELSE
1186	IF(lwp) WRITE(numout,*) ' ln_ts_auto=.false.: Use nn_baro in namelist nn_baro = ', nn_baro
1187	ENDIF
1188
1189	IF(ln_bt_av) THEN
1190	IF(lwp) WRITE(numout,*) ' ln_bt_av =.true. ==> Time averaging over nn_baro time steps is on '
1191	ELSE
1192	IF(lwp) WRITE(numout,*) ' ln_bt_av =.false. => No time averaging of barotropic variables '
1193	ENDIF
1194	!
1195	!
1196	IF(ln_bt_fw) THEN
1197	IF(lwp) WRITE(numout,*) ' ln_bt_fw=.true. => Forward integration of barotropic variables '
1198	ELSE
1199	IF(lwp) WRITE(numout,*) ' ln_bt_fw =.false.=> Centred integration of barotropic variables '
1200	ENDIF
1201	!
1202	#if defined key_agrif
1203	! Restrict the use of Agrif to the forward case only
1204	!!! IF( .NOT.ln_bt_fw .AND. .NOT.Agrif_Root() ) CALL ctl_stop( 'AGRIF not implemented if ln_bt_fw=.FALSE.' )
1205	#endif
1206	!
1207	IF(lwp) WRITE(numout,*) ' Time filter choice, nn_bt_flt: ', nn_bt_flt
1208	SELECT CASE ( nn_bt_flt )
1209	CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' Dirac'
1210	CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' Boxcar: width = nn_baro'
1211	CASE( 2 ) ; IF(lwp) WRITE(numout,) ' Boxcar: width = 2nn_baro'
1212	CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_bt_flt: should 0,1, or 2' )
1213	END SELECT
1214	!
1215	IF(lwp) WRITE(numout,*) ' '
1216	IF(lwp) WRITE(numout,*) ' nn_baro = ', nn_baro
1217	IF(lwp) WRITE(numout,*) ' Barotropic time step [s] is :', rdtbt
1218	IF(lwp) WRITE(numout,*) ' Maximum Courant number is :', zcmax
1219	!
1220	IF(lwp) WRITE(numout,*) ' Time diffusion parameter rn_bt_alpha: ', rn_bt_alpha
1221	IF ((ln_bt_av.AND.nn_bt_flt/=0).AND.(rn_bt_alpha>0._wp)) THEN
1222	CALL ctl_stop( 'dynspg_ts ERROR: if rn_bt_alpha > 0, remove temporal averaging' )
1223	ENDIF
1224	!
1225	IF( .NOT.ln_bt_av .AND. .NOT.ln_bt_fw ) THEN
1226	CALL ctl_stop( 'dynspg_ts ERROR: No time averaging => only forward integration is possible' )
1227	ENDIF
1228	IF( zcmax>0.9_wp ) THEN
1229	CALL ctl_stop( 'dynspg_ts ERROR: Maximum Courant number is greater than 0.9: Inc. nn_baro !' )
1230	ENDIF
1231	!
1232	! ! Allocate time-splitting arrays
1233	IF( dyn_spg_ts_alloc() /= 0 ) CALL ctl_stop('STOP', 'dyn_spg_init: failed to allocate dynspg_ts arrays' )
1234	!
1235	! ! read restart when needed
1236	CALL ts_rst( nit000, 'READ' )
1237	!
1238	IF( lwxios ) THEN
1239	! define variables in restart file when writing with XIOS
1240	CALL iom_set_rstw_var_active('ub2_b')
1241	CALL iom_set_rstw_var_active('vb2_b')
1242	CALL iom_set_rstw_var_active('un_bf')
1243	CALL iom_set_rstw_var_active('vn_bf')
1244	!
1245	IF (.NOT.ln_bt_av) THEN
1246	CALL iom_set_rstw_var_active('sshbb_e')
1247	CALL iom_set_rstw_var_active('ubb_e')
1248	CALL iom_set_rstw_var_active('vbb_e')
1249	CALL iom_set_rstw_var_active('sshb_e')
1250	CALL iom_set_rstw_var_active('ub_e')
1251	CALL iom_set_rstw_var_active('vb_e')
1252	ENDIF
1253	#if defined key_agrif
1254	! Save time integrated fluxes
1255	IF ( .NOT.Agrif_Root() ) THEN
1256	CALL iom_set_rstw_var_active('ub2_i_b')
1257	CALL iom_set_rstw_var_active('vb2_i_b')
1258	ENDIF
1259	#endif
1260	ENDIF
1261	!
1262	END SUBROUTINE dyn_spg_ts_init
1263
1264
1265	SUBROUTINE dyn_cor_2d_init
1266	!!---------------------------------------------------------------------
1267	!! * ROUTINE dyn_cor_2d_init *
1268	!!
1269	!! ** Purpose : Set time splitting options
1270	!! Set arrays to remove/compute coriolis trend.
1271	!! Do it once during initialization if volume is fixed, else at each long time step.
1272	!! Note that these arrays are also used during barotropic loop. These are however frozen
1273	!! although they should be updated in the variable volume case. Not a big approximation.
1274	!! To remove this approximation, copy lines below inside barotropic loop
1275	!! and update depths at T-F points (ht and zhf resp.) at each barotropic time step
1276	!!
1277	!! Compute zwz = f / ( height of the water colomn )
1278	!!----------------------------------------------------------------------
1279	INTEGER :: ji ,jj, jk ! dummy loop indices
1280	REAL(wp) :: z1_ht
1281	REAL(wp), DIMENSION(jpi,jpj) :: zhf
1282	!!----------------------------------------------------------------------
1283	!
1284	SELECT CASE( nvor_scheme )
1285	CASE( np_EEN ) != EEN scheme using e3f (energy & enstrophy scheme)
1286	SELECT CASE( nn_een_e3f ) !* ff_f/e3 at F-point
1287	CASE ( 0 ) ! original formulation (masked averaging of e3t divided by 4)
1288	DO jj = 1, jpjm1
1289	DO ji = 1, jpim1
1290	zwz(ji,jj) = ( ht_n(ji ,jj+1) + ht_n(ji+1,jj+1) + &
1291	& ht_n(ji ,jj ) + ht_n(ji+1,jj ) ) * 0.25_wp
1292	IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj)
1293	END DO
1294	END DO
1295	CASE ( 1 ) ! new formulation (masked averaging of e3t divided by the sum of mask)
1296	DO jj = 1, jpjm1
1297	DO ji = 1, jpim1
1298	zwz(ji,jj) = ( ht_n (ji ,jj+1) + ht_n (ji+1,jj+1) &
1299	& + ht_n (ji ,jj ) + ht_n (ji+1,jj ) ) &
1300	& / ( MAX( 1._wp, ssmask(ji ,jj+1) + ssmask(ji+1,jj+1) &
1301	& + ssmask(ji ,jj ) + ssmask(ji+1,jj ) ) )
1302	IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj)
1303	END DO
1304	END DO
1305	END SELECT
1306	CALL lbc_lnk( 'dynspg_ts', zwz, 'F', 1._wp )
1307	!
1308	ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp
1309	DO jj = 2, jpj
1310	DO ji = 2, jpi
1311	ftne(ji,jj) = zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1)
1312	ftnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj ) + zwz(ji ,jj )
1313	ftse(ji,jj) = zwz(ji ,jj ) + zwz(ji ,jj-1) + zwz(ji-1,jj-1)
1314	ftsw(ji,jj) = zwz(ji ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj )
1315	END DO
1316	END DO
1317	!
1318	CASE( np_EET ) != EEN scheme using e3t (energy conserving scheme)
1319	ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp
1320	DO jj = 2, jpj
1321	DO ji = 2, jpi
1322	z1_ht = ssmask(ji,jj) / ( ht_n(ji,jj) + 1._wp - ssmask(ji,jj) )
1323	ftne(ji,jj) = ( ff_f(ji-1,jj ) + ff_f(ji ,jj ) + ff_f(ji ,jj-1) ) * z1_ht
1324	ftnw(ji,jj) = ( ff_f(ji-1,jj-1) + ff_f(ji-1,jj ) + ff_f(ji ,jj ) ) * z1_ht
1325	ftse(ji,jj) = ( ff_f(ji ,jj ) + ff_f(ji ,jj-1) + ff_f(ji-1,jj-1) ) * z1_ht
1326	ftsw(ji,jj) = ( ff_f(ji ,jj-1) + ff_f(ji-1,jj-1) + ff_f(ji-1,jj ) ) * z1_ht
1327	END DO
1328	END DO
1329	!
1330	CASE( np_ENE, np_ENS , np_MIX ) != all other schemes (ENE, ENS, MIX) except ENT !
1331	!
1332	zwz(:,:) = 0._wp
1333	zhf(:,:) = 0._wp
1334
1335	!!gm assume 0 in both cases (which is almost surely WRONG ! ) as hvatf has been removed
1336	!!gm A priori a better value should be something like :
1337	!!gm zhf(i,j) = masked sum of ht(i,j) , ht(i+1,j) , ht(i,j+1) , (i+1,j+1)
1338	!!gm divided by the sum of the corresponding mask
1339	!!gm
1340	!!
1341	IF( .NOT.ln_sco ) THEN
1342
1343	!!gm agree the JC comment : this should be done in a much clear way
1344
1345	! JC: It not clear yet what should be the depth at f-points over land in z-coordinate case
1346	! Set it to zero for the time being
1347	! IF( rn_hmin < 0._wp ) THEN ; jk = - INT( rn_hmin ) ! from a nb of level
1348	! ELSE ; jk = MINLOC( gdepw_0, mask = gdepw_0 > rn_hmin, dim = 1 ) ! from a depth
1349	! ENDIF
1350	! zhf(:,:) = gdepw_0(:,:,jk+1)
1351	!
1352	ELSE
1353	!
1354	!zhf(:,:) = hbatf(:,:)
1355	DO jj = 1, jpjm1
1356	DO ji = 1, jpim1
1357	zhf(ji,jj) = ( ht_0 (ji,jj ) + ht_0 (ji+1,jj ) &
1358	& + ht_0 (ji,jj+1) + ht_0 (ji+1,jj+1) ) &
1359	& / MAX( ssmask(ji,jj ) + ssmask(ji+1,jj ) &
1360	& + ssmask(ji,jj+1) + ssmask(ji+1,jj+1) , 1._wp )
1361	END DO
1362	END DO
1363	ENDIF
1364	!
1365	DO jj = 1, jpjm1
1366	zhf(:,jj) = zhf(:,jj) * (1._wp- umask(:,jj,1) * umask(:,jj+1,1))
1367	END DO
1368	!
1369	DO jk = 1, jpkm1
1370	DO jj = 1, jpjm1
1371	zhf(:,jj) = zhf(:,jj) + e3f_n(:,jj,jk) * umask(:,jj,jk) * umask(:,jj+1,jk)
1372	END DO
1373	END DO
1374	CALL lbc_lnk( 'dynspg_ts', zhf, 'F', 1._wp )
1375	! JC: TBC. hf should be greater than 0
1376	DO jj = 1, jpj
1377	DO ji = 1, jpi
1378	IF( zhf(ji,jj) /= 0._wp ) zwz(ji,jj) = 1._wp / zhf(ji,jj)
1379	END DO
1380	END DO
1381	zwz(:,:) = ff_f(:,:) * zwz(:,:)
1382	END SELECT
1383
1384	END SUBROUTINE dyn_cor_2d_init
1385
1386
1387
1388	SUBROUTINE dyn_cor_2d( hu_n, hv_n, un_b, vn_b, zhU, zhV, zu_trd, zv_trd )
1389	!!---------------------------------------------------------------------
1390	!! * ROUTINE dyn_cor_2d *
1391	!!
1392	!! ** Purpose : Compute u and v coriolis trends
1393	!!----------------------------------------------------------------------
1394	INTEGER :: ji ,jj ! dummy loop indices
1395	REAL(wp) :: zx1, zx2, zy1, zy2, z1_hu, z1_hv ! - -
1396	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: hu_n, hv_n, un_b, vn_b, zhU, zhV
1397	REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: zu_trd, zv_trd
1398	!!----------------------------------------------------------------------
1399	SELECT CASE( nvor_scheme )
1400	CASE( np_ENT ) ! enstrophy conserving scheme (f-point)
1401	DO jj = 2, jpjm1
1402	DO ji = 2, jpim1
1403	z1_hu = ssumask(ji,jj) / ( hu_n(ji,jj) + 1._wp - ssumask(ji,jj) )
1404	z1_hv = ssvmask(ji,jj) / ( hv_n(ji,jj) + 1._wp - ssvmask(ji,jj) )
1405	zu_trd(ji,jj) = + r1_4 * r1_e1e2u(ji,jj) * z1_hu &
1406	& * ( e1e2t(ji+1,jj)ht_n(ji+1,jj)ff_t(ji+1,jj) * ( vn_b(ji+1,jj) + vn_b(ji+1,jj-1) ) &
1407	& + e1e2t(ji ,jj)ht_n(ji ,jj)ff_t(ji ,jj) * ( vn_b(ji ,jj) + vn_b(ji ,jj-1) ) )
1408	!
1409	zv_trd(ji,jj) = - r1_4 * r1_e1e2v(ji,jj) * z1_hv &
1410	& * ( e1e2t(ji,jj+1)ht_n(ji,jj+1)ff_t(ji,jj+1) * ( un_b(ji,jj+1) + un_b(ji-1,jj+1) ) &
1411	& + e1e2t(ji,jj )ht_n(ji,jj )ff_t(ji,jj ) * ( un_b(ji,jj ) + un_b(ji-1,jj ) ) )
1412	END DO
1413	END DO
1414	!
1415	CASE( np_ENE , np_MIX ) ! energy conserving scheme (t-point) ENE or MIX
1416	DO jj = 2, jpjm1
1417	DO ji = fs_2, fs_jpim1 ! vector opt.
1418	zy1 = ( zhV(ji,jj-1) + zhV(ji+1,jj-1) ) * r1_e1u(ji,jj)
1419	zy2 = ( zhV(ji,jj ) + zhV(ji+1,jj ) ) * r1_e1u(ji,jj)
1420	zx1 = ( zhU(ji-1,jj) + zhU(ji-1,jj+1) ) * r1_e2v(ji,jj)
1421	zx2 = ( zhU(ji ,jj) + zhU(ji ,jj+1) ) * r1_e2v(ji,jj)
1422	! energy conserving formulation for planetary vorticity term
1423	zu_trd(ji,jj) = r1_4 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
1424	zv_trd(ji,jj) = - r1_4 * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 )
1425	END DO
1426	END DO
1427	!
1428	CASE( np_ENS ) ! enstrophy conserving scheme (f-point)
1429	DO jj = 2, jpjm1
1430	DO ji = fs_2, fs_jpim1 ! vector opt.
1431	zy1 = r1_8 * ( zhV(ji ,jj-1) + zhV(ji+1,jj-1) &
1432	& + zhV(ji ,jj ) + zhV(ji+1,jj ) ) * r1_e1u(ji,jj)
1433	zx1 = - r1_8 * ( zhU(ji-1,jj ) + zhU(ji-1,jj+1) &
1434	& + zhU(ji ,jj ) + zhU(ji ,jj+1) ) * r1_e2v(ji,jj)
1435	zu_trd(ji,jj) = zy1 * ( zwz(ji ,jj-1) + zwz(ji,jj) )
1436	zv_trd(ji,jj) = zx1 * ( zwz(ji-1,jj ) + zwz(ji,jj) )
1437	END DO
1438	END DO
1439	!
1440	CASE( np_EET , np_EEN ) ! energy & enstrophy scheme (using e3t or e3f)
1441	DO jj = 2, jpjm1
1442	DO ji = fs_2, fs_jpim1 ! vector opt.
1443	zu_trd(ji,jj) = + r1_12 * r1_e1u(ji,jj) * ( ftne(ji,jj ) * zhV(ji ,jj ) &
1444	& + ftnw(ji+1,jj) * zhV(ji+1,jj ) &
1445	& + ftse(ji,jj ) * zhV(ji ,jj-1) &
1446	& + ftsw(ji+1,jj) * zhV(ji+1,jj-1) )
1447	zv_trd(ji,jj) = - r1_12 * r1_e2v(ji,jj) * ( ftsw(ji,jj+1) * zhU(ji-1,jj+1) &
1448	& + ftse(ji,jj+1) * zhU(ji ,jj+1) &
1449	& + ftnw(ji,jj ) * zhU(ji-1,jj ) &
1450	& + ftne(ji,jj ) * zhU(ji ,jj ) )
1451	END DO
1452	END DO
1453	!
1454	END SELECT
1455	!
1456	END SUBROUTINE dyn_cor_2D
1457
1458
1459	SUBROUTINE wad_tmsk( pssh, ptmsk )
1460	!!----------------------------------------------------------------------
1461	!! * ROUTINE wad_lmt *
1462	!!
1463	!! ** Purpose : set wetting & drying mask at tracer points
1464	!! for the current barotropic sub-step
1465	!!
1466	!! ** Method : ???
1467	!!
1468	!! ** Action : ptmsk : wetting & drying t-mask
1469	!!----------------------------------------------------------------------
1470	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pssh !
1471	REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: ptmsk !
1472	!
1473	INTEGER :: ji, jj ! dummy loop indices
1474	!!----------------------------------------------------------------------
1475	!
1476	IF( ln_wd_dl_rmp ) THEN
1477	DO jj = 1, jpj
1478	DO ji = 1, jpi
1479	IF ( pssh(ji,jj) + ht_0(ji,jj) > 2._wp * rn_wdmin1 ) THEN
1480	! IF ( pssh(ji,jj) + ht_0(ji,jj) > rn_wdmin2 ) THEN
1481	ptmsk(ji,jj) = 1._wp
1482	ELSEIF( pssh(ji,jj) + ht_0(ji,jj) > rn_wdmin1 ) THEN
1483	ptmsk(ji,jj) = TANH( 50._wp( ( pssh(ji,jj) + ht_0(ji,jj) - rn_wdmin1 )r_rn_wdmin1) )
1484	ELSE
1485	ptmsk(ji,jj) = 0._wp
1486	ENDIF
1487	END DO
1488	END DO
1489	ELSE
1490	DO jj = 1, jpj
1491	DO ji = 1, jpi
1492	IF ( pssh(ji,jj) + ht_0(ji,jj) > rn_wdmin1 ) THEN ; ptmsk(ji,jj) = 1._wp
1493	ELSE ; ptmsk(ji,jj) = 0._wp
1494	ENDIF
1495	END DO
1496	END DO
1497	ENDIF
1498	!
1499	END SUBROUTINE wad_tmsk
1500
1501
1502	SUBROUTINE wad_Umsk( pTmsk, phU, phV, pu, pv, pUmsk, pVmsk )
1503	!!----------------------------------------------------------------------
1504	!! * ROUTINE wad_lmt *
1505	!!
1506	!! ** Purpose : set wetting & drying mask at tracer points
1507	!! for the current barotropic sub-step
1508	!!
1509	!! ** Method : ???
1510	!!
1511	!! ** Action : ptmsk : wetting & drying t-mask
1512	!!----------------------------------------------------------------------
1513	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pTmsk ! W & D t-mask
1514	REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: phU, phV, pu, pv ! ocean velocities and transports
1515	REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pUmsk, pVmsk ! W & D u- and v-mask
1516	!
1517	INTEGER :: ji, jj ! dummy loop indices
1518	!!----------------------------------------------------------------------
1519	!
1520	DO jj = 1, jpj
1521	DO ji = 1, jpim1 ! not jpi-column
1522	IF ( phU(ji,jj) > 0._wp ) THEN ; pUmsk(ji,jj) = pTmsk(ji ,jj)
1523	ELSE ; pUmsk(ji,jj) = pTmsk(ji+1,jj)
1524	ENDIF
1525	phU(ji,jj) = pUmsk(ji,jj)*phU(ji,jj)
1526	pu (ji,jj) = pUmsk(ji,jj)*pu (ji,jj)
1527	END DO
1528	END DO
1529	!
1530	DO jj = 1, jpjm1 ! not jpj-row
1531	DO ji = 1, jpi
1532	IF ( phV(ji,jj) > 0._wp ) THEN ; pVmsk(ji,jj) = pTmsk(ji,jj )
1533	ELSE ; pVmsk(ji,jj) = pTmsk(ji,jj+1)
1534	ENDIF
1535	phV(ji,jj) = pVmsk(ji,jj)*phV(ji,jj)
1536	pv (ji,jj) = pVmsk(ji,jj)*pv (ji,jj)
1537	END DO
1538	END DO
1539	!
1540	END SUBROUTINE wad_Umsk
1541
1542
1543	SUBROUTINE wad_spg( sshn, zcpx, zcpy )
1544	!!---------------------------------------------------------------------
1545	!! * ROUTINE wad_sp *
1546	!!
1547	!! ** Purpose :
1548	!!----------------------------------------------------------------------
1549	INTEGER :: ji ,jj ! dummy loop indices
1550	LOGICAL :: ll_tmp1, ll_tmp2
1551	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: sshn
1552	REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: zcpx, zcpy
1553	!!----------------------------------------------------------------------
1554	DO jj = 2, jpjm1
1555	DO ji = 2, jpim1
1556	ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji+1,jj) ) > &
1557	& MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) .AND. &
1558	& MAX( sshn(ji,jj) + ht_0(ji,jj) , sshn(ji+1,jj) + ht_0(ji+1,jj) ) &
1559	& > rn_wdmin1 + rn_wdmin2
1560	ll_tmp2 = ( ABS( sshn(ji+1,jj) - sshn(ji ,jj)) > 1.E-12 ).AND.( &
1561	& MAX( sshn(ji,jj) , sshn(ji+1,jj) ) > &
1562	& MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 )
1563	IF(ll_tmp1) THEN
1564	zcpx(ji,jj) = 1.0_wp
1565	ELSEIF(ll_tmp2) THEN
1566	! no worries about sshn(ji+1,jj) - sshn(ji ,jj) = 0, it won't happen ! here
1567	zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_0(ji+1,jj) - sshn(ji,jj) - ht_0(ji,jj)) &
1568	& / (sshn(ji+1,jj) - sshn(ji ,jj)) )
1569	zcpx(ji,jj) = max(min( zcpx(ji,jj) , 1.0_wp),0.0_wp)
1570	ELSE
1571	zcpx(ji,jj) = 0._wp
1572	ENDIF
1573	!
1574	ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji,jj+1) ) > &
1575	& MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) .AND. &
1576	& MAX( sshn(ji,jj) + ht_0(ji,jj) , sshn(ji,jj+1) + ht_0(ji,jj+1) ) &
1577	& > rn_wdmin1 + rn_wdmin2
1578	ll_tmp2 = ( ABS( sshn(ji,jj) - sshn(ji,jj+1)) > 1.E-12 ).AND.( &
1579	& MAX( sshn(ji,jj) , sshn(ji,jj+1) ) > &
1580	& MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 )
1581
1582	IF(ll_tmp1) THEN
1583	zcpy(ji,jj) = 1.0_wp
1584	ELSE IF(ll_tmp2) THEN
1585	! no worries about sshn(ji,jj+1) - sshn(ji,jj ) = 0, it won't happen ! here
1586	zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_0(ji,jj+1) - sshn(ji,jj) - ht_0(ji,jj)) &
1587	& / (sshn(ji,jj+1) - sshn(ji,jj )) )
1588	zcpy(ji,jj) = MAX( 0._wp , MIN( zcpy(ji,jj) , 1.0_wp ) )
1589	ELSE
1590	zcpy(ji,jj) = 0._wp
1591	ENDIF
1592	END DO
1593	END DO
1594
1595	END SUBROUTINE wad_spg
1596
1597
1598
1599	SUBROUTINE dyn_drg_init( pu_RHSi, pv_RHSi, pCdU_u, pCdU_v, &
1600	& ptfrtrdu, ptfrtrdv, pbfrtrdu, pbfrtrdv )
1601	!!----------------------------------------------------------------------
1602	!! * ROUTINE dyn_drg_init *
1603	!!
1604	!! ** Purpose : - add the baroclinic top/bottom drag contribution to
1605	!! the baroclinic part of the barotropic RHS
1606	!! - compute the barotropic drag coefficients
1607	!!
1608	!! ** Method : computation done over the INNER domain only
1609	!!----------------------------------------------------------------------
1610	REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pu_RHSi, pv_RHSi ! baroclinic part of the barotropic RHS
1611	REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: pCdU_u , pCdU_v ! barotropic drag coefficients
1612	REAL(wp), DIMENSION(jpi,jpj), INTENT(inout), OPTIONAL :: ptfrtrdu, ptfrtrdv ! top friction trends
1613	REAL(wp), DIMENSION(jpi,jpj), INTENT(inout), OPTIONAL :: pbfrtrdu, pbfrtrdv ! bottom friction trends
1614	!
1615	INTEGER :: ji, jj ! dummy loop indices
1616	INTEGER :: ikbu, ikbv, iktu, iktv
1617	REAL(wp) :: zztmp
1618	REAL(wp), DIMENSION(jpi,jpj) :: zu_i, zv_i
1619	!!----------------------------------------------------------------------
1620	!
1621	! !== Set the barotropic drag coef. ==!
1622	!
1623	IF( ln_isfcav.OR.ln_drgice_imp ) THEN ! top+bottom friction (ocean cavities)
1624
1625	DO jj = 2, jpjm1
1626	DO ji = 2, jpim1 ! INNER domain
1627	pCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) + rCdU_top(ji+1,jj)+rCdU_top(ji,jj) )
1628	pCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) + rCdU_top(ji,jj+1)+rCdU_top(ji,jj) )
1629	END DO
1630	END DO
1631	ELSE ! bottom friction only
1632	DO jj = 2, jpjm1
1633	DO ji = 2, jpim1 ! INNER domain
1634	pCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) )
1635	pCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) )
1636	END DO
1637	END DO
1638	ENDIF
1639	!
1640	! !== BOTTOM stress contribution from baroclinic velocities ==!
1641	!
1642	IF( l_trddyn ) THEN
1643	pbfrtrdu(:,:) = pu_RHSi(:,:)
1644	pbfrtrdv(:,:) = pv_RHSi(:,:)
1645	ENDIF
1646	!
1647	IF( ln_bt_fw ) THEN ! FORWARD integration: use NOW bottom baroclinic velocities
1648
1649	DO jj = 2, jpjm1
1650	DO ji = 2, jpim1 ! INNER domain
1651	ikbu = mbku(ji,jj)
1652	ikbv = mbkv(ji,jj)
1653	zu_i(ji,jj) = un(ji,jj,ikbu) - un_b(ji,jj)
1654	zv_i(ji,jj) = vn(ji,jj,ikbv) - vn_b(ji,jj)
1655	END DO
1656	END DO
1657	ELSE ! CENTRED integration: use BEFORE bottom baroclinic velocities
1658
1659	DO jj = 2, jpjm1
1660	DO ji = 2, jpim1 ! INNER domain
1661	ikbu = mbku(ji,jj)
1662	ikbv = mbkv(ji,jj)
1663	zu_i(ji,jj) = ub(ji,jj,ikbu) - ub_b(ji,jj)
1664	zv_i(ji,jj) = vb(ji,jj,ikbv) - vb_b(ji,jj)
1665	END DO
1666	END DO
1667	ENDIF
1668	!
1669	IF( ln_wd_il ) THEN ! W/D : use the "clipped" bottom friction !!gm explain WHY, please !
1670	zztmp = -1._wp / rdtbt
1671	DO jj = 2, jpjm1
1672	DO ji = 2, jpim1 ! INNER domain
1673	pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + zu_i(ji,jj) * wdrampu(ji,jj) * MAX( &
1674	& r1_hu_n(ji,jj) * r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) , zztmp )
1675	pv_RHSi(ji,jj) = pv_RHSi(ji,jj) + zv_i(ji,jj) * wdrampv(ji,jj) * MAX( &
1676	& r1_hv_n(ji,jj) * r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) , zztmp )
1677	END DO
1678	END DO
1679	ELSE ! use "unclipped" drag (even if explicit friction is used in 3D calculation)
1680
1681	DO jj = 2, jpjm1
1682	DO ji = 2, jpim1 ! INNER domain
1683	pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + r1_hu_n(ji,jj) * r1_2( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) zu_i(ji,jj)
1684	pv_RHSi(ji,jj) = pv_RHSi(ji,jj) + r1_hv_n(ji,jj) * r1_2( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) zv_i(ji,jj)
1685	END DO
1686	END DO
1687	END IF
1688	!
1689	IF( l_trddyn ) THEN
1690	pbfrtrdu(:,:) = pu_RHSi(:,:) - pbfrtrdu(:,:)
1691	pbfrtrdv(:,:) = pv_RHSi(:,:) - pbfrtrdv(:,:)
1692	ENDIF
1693	!
1694	! !== TOP stress contribution from baroclinic velocities ==! (no W/D case)
1695	!
1696	IF( ln_isfcav.OR.ln_drgice_imp ) THEN
1697	!
1698	IF( l_trddyn ) THEN
1699	ptfrtrdu(:,:) = pu_RHSi(:,:)
1700	ptfrtrdv(:,:) = pv_RHSi(:,:)
1701	ENDIF
1702	!
1703	IF( ln_bt_fw ) THEN ! FORWARD integration: use NOW top baroclinic velocity
1704
1705	DO jj = 2, jpjm1
1706	DO ji = 2, jpim1 ! INNER domain
1707	iktu = miku(ji,jj)
1708	iktv = mikv(ji,jj)
1709	zu_i(ji,jj) = un(ji,jj,iktu) - un_b(ji,jj)
1710	zv_i(ji,jj) = vn(ji,jj,iktv) - vn_b(ji,jj)
1711	END DO
1712	END DO
1713	ELSE ! CENTRED integration: use BEFORE top baroclinic velocity
1714
1715	DO jj = 2, jpjm1
1716	DO ji = 2, jpim1 ! INNER domain
1717	iktu = miku(ji,jj)
1718	iktv = mikv(ji,jj)
1719	zu_i(ji,jj) = ub(ji,jj,iktu) - ub_b(ji,jj)
1720	zv_i(ji,jj) = vb(ji,jj,iktv) - vb_b(ji,jj)
1721	END DO
1722	END DO
1723	ENDIF
1724	!
1725	! ! use "unclipped" top drag (even if explicit friction is used in 3D calculation)
1726
1727	DO jj = 2, jpjm1
1728	DO ji = 2, jpim1 ! INNER domain
1729	pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + r1_hu_n(ji,jj) * r1_2( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) zu_i(ji,jj)
1730	pv_RHSi(ji,jj) = pv_RHSi(ji,jj) + r1_hv_n(ji,jj) * r1_2( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) zv_i(ji,jj)
1731	END DO
1732	END DO
1733	!
1734	IF( l_trddyn ) THEN
1735	ptfrtrdu(:,:) = pu_RHSi(:,:) - ptfrtrdu(:,:)
1736	ptfrtrdv(:,:) = pv_RHSi(:,:) - ptfrtrdv(:,:)
1737	ENDIF
1738	!
1739	ENDIF
1740	!
1741	END SUBROUTINE dyn_drg_init
1742
1743	SUBROUTINE ts_bck_interp( jn, ll_init, & ! <== in
1744	& za0, za1, za2, za3 ) ! ==> out
1745	!!----------------------------------------------------------------------
1746	INTEGER ,INTENT(in ) :: jn ! index of sub time step
1747	LOGICAL ,INTENT(in ) :: ll_init !
1748	REAL(wp),INTENT( out) :: za0, za1, za2, za3 ! Half-step back interpolation coefficient
1749	!
1750	REAL(wp) :: zepsilon, zgamma ! - -
1751	!!----------------------------------------------------------------------
1752	! ! set Half-step back interpolation coefficient
1753	IF ( jn==1 .AND. ll_init ) THEN !* Forward-backward
1754	za0 = 1._wp
1755	za1 = 0._wp
1756	za2 = 0._wp
1757	za3 = 0._wp
1758	ELSEIF( jn==2 .AND. ll_init ) THEN !* AB2-AM3 Coefficients; bet=0 ; gam=-1/6 ; eps=1/12
1759	za0 = 1.0833333333333_wp ! za0 = 1-gam-eps
1760	za1 =-0.1666666666666_wp ! za1 = gam
1761	za2 = 0.0833333333333_wp ! za2 = eps
1762	za3 = 0._wp
1763	ELSE !* AB3-AM4 Coefficients; bet=0.281105 ; eps=0.013 ; gam=0.0880
1764	IF( rn_bt_alpha == 0._wp ) THEN ! Time diffusion
1765	za0 = 0.614_wp ! za0 = 1/2 + gam + 2*eps
1766	za1 = 0.285_wp ! za1 = 1/2 - 2gam - 3eps
1767	za2 = 0.088_wp ! za2 = gam
1768	za3 = 0.013_wp ! za3 = eps
1769	ELSE ! no time diffusion
1770	zepsilon = 0.00976186_wp - 0.13451357_wp * rn_bt_alpha
1771	zgamma = 0.08344500_wp - 0.51358400_wp * rn_bt_alpha
1772	za0 = 0.5_wp + zgamma + 2._wp * rn_bt_alpha + 2._wp * zepsilon
1773	za1 = 1._wp - za0 - zgamma - zepsilon
1774	za2 = zgamma
1775	za3 = zepsilon
1776	ENDIF
1777	ENDIF
1778	END SUBROUTINE ts_bck_interp
1779
1780
1781	!!======================================================================
1782	END MODULE dynspg_ts

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

New URL for NEMO forge! http://forge.nemo-ocean.eu

Context Navigation

source: NEMO/branches/UKMO/NEMO_4.0.4_momentum_trends/src/OCE/DYN/dynspg_ts.F90 @ 15194

Download in other formats: