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 @ 11949

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

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

Download in other formats: