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_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/DYN – NEMO

Context Navigation

source: NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/DYN/dynspg_ts.F90 @ 11240

Last change on this file since 11240 was 11240, checked in by girrmann, 5 years ago
dev_r10984_HPC-13 : remove a communication in dyn_spg_ts, see #2285
Property svn:keywords set to `Id`
File size: 80.5 KB

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

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

Download in other formats: