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

Context Navigation

source: NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/DYN/dynspg_ts.F90 @ 12252

Last change on this file since 12252 was 12252, checked in by smasson, 4 years ago
rev12240_dev_r11943_MERGE_2019: same as [12251], merge trunk 12072:12248, all sette tests ok, GYRE_PISCES, AMM12, ISOMIP, VORTEX intentical to 12236
Property svn:keywords set to `Id`
File size: 82.1 KB

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

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

Download in other formats: