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 NEMO/branches/UKMO/NEMO_4.0.1_MIRROR_WAD_ZENV/src/OCE/DYN – NEMO

Context Navigation

source: NEMO/branches/UKMO/NEMO_4.0.1_MIRROR_WAD_ZENV/src/OCE/DYN/dynspg_ts.F90 @ 12095

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

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

Download in other formats: