Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

dynspg_ts.F90

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

Download in other formats: