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

Context Navigation

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

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

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

Download in other formats: