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

Context Navigation

source: NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/src/OCE/DYN/dynspg_ts.F90 @ 12166

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

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

Download in other formats: