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/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DYN – NEMO

Context Navigation

source: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DYN/dynspg_ts.F90 @ 12624

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

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

Download in other formats: