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/2019/dev_r11756_SI3restart_XIOS/src/OCE/DYN – NEMO

Context Navigation

source: NEMO/branches/2019/dev_r11756_SI3restart_XIOS/src/OCE/DYN/dynspg_ts.F90 @ 11837

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

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

Download in other formats: