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_r10973_AGRIF-01_jchanut_small_jpi_jpj/src/OCE/DYN – NEMO

Context Navigation

source: NEMO/branches/2019/dev_r10973_AGRIF-01_jchanut_small_jpi_jpj/src/OCE/DYN/dynspg_ts.F90 @ 11564

Last change on this file since 11564 was 11564, checked in by jchanut, 5 years ago
#2199, merged with trunk
Property svn:keywords set to `Id`
File size: 79.9 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( 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 dyn_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, jpjm1 ! not jpj-row
471	DO ji = 1, jpi
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 ) CALL agrif_dyn_ts_flux( jn, zhU, zhV )
493	#endif
494	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
495
496	IF( ln_wd_dl ) THEN ! un_e and vn_e are set to zero at faces where
497	! ! the direction of the flow is from dry cells
498	CALL wad_Umsk( ztwdmask, zhU, zhV, un_e, vn_e, zuwdmask, zvwdmask ) ! not jpi colomn for U, not jpj row for V
499	!
500	ENDIF
501	!
502	!
503	! Compute Sea Level at step jit+1
504	!-- m+1 m m+1/2 --!
505	!-- ssh = ssh - delta_t' * [ frc + div( flux ) ] --!
506	!-------------------------------------------------------------------------!
507	DO jj = 2, jpjm1 ! INNER domain
508	DO ji = 2, jpim1
509	zhdiv = ( zhU(ji,jj) - zhU(ji-1,jj) + zhV(ji,jj) - zhV(ji,jj-1) ) * r1_e1e2t(ji,jj)
510	ssha_e(ji,jj) = ( sshn_e(ji,jj) - rdtbt * ( zssh_frc(ji,jj) + zhdiv ) ) * ssmask(ji,jj)
511	END DO
512	END DO
513	!
514	CALL lbc_lnk_multi( 'dynspg_ts', ssha_e, 'T', 1._wp, zhU, 'U', -1._wp, zhV, 'V', -1._wp )
515	!
516	! ! Sum over sub-time-steps to compute advective velocities
517	za2 = wgtbtp2(jn) ! zhU, zhV hold fluxes extrapolated at jn+0.5
518	un_adv(:,:) = un_adv(:,:) + za2 * zhU(:,:) * r1_e2u(:,:)
519	vn_adv(:,:) = vn_adv(:,:) + za2 * zhV(:,:) * r1_e1v(:,:)
520	! sum over sub-time-steps to decide which baroclinic velocities to set to zero (zuwdav2 is only used when ln_wd_dl_bc=True)
521	IF ( ln_wd_dl_bc ) THEN
522	zuwdav2(1:jpim1,1:jpj ) = zuwdav2(1:jpim1,1:jpj ) + za2 * zuwdmask(1:jpim1,1:jpj ) ! not jpi-column
523	zvwdav2(1:jpi ,1:jpjm1) = zvwdav2(1:jpi ,1:jpjm1) + za2 * zvwdmask(1:jpi ,1:jpjm1) ! not jpj-row
524	END IF
525	!
526	! Duplicate sea level across open boundaries (this is only cosmetic if linssh=T)
527	IF( ln_bdy ) CALL bdy_ssh( ssha_e )
528	#if defined key_agrif
529	IF( .NOT.Agrif_Root() ) CALL agrif_ssh_ts( jn )
530	#endif
531	!
532	! Sea Surface Height at u-,v-points (vvl case only)
533	IF( .NOT.ln_linssh ) THEN
534	DO jj = 2, jpjm1 ! INNER domain, will be extended to whole domain later
535	DO ji = 2, jpim1 ! NO Vector Opt.
536	zsshu_a(ji,jj) = r1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj) &
537	& * ( e1e2t(ji ,jj ) * ssha_e(ji ,jj ) &
538	& + e1e2t(ji+1,jj ) * ssha_e(ji+1,jj ) )
539	zsshv_a(ji,jj) = r1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj) &
540	& * ( e1e2t(ji ,jj ) * ssha_e(ji ,jj ) &
541	& + e1e2t(ji ,jj+1) * ssha_e(ji ,jj+1) )
542	END DO
543	END DO
544	ENDIF
545	!
546	! Half-step back interpolation of SSH for surface pressure computation at step jit+1/2
547	!-- m+1/2 m+1 m m-1 m-2 --!
548	!-- ssh' = za0 * ssh + za1 * ssh + za2 * ssh + za3 * ssh --!
549	!------------------------------------------------------------------------------------------!
550	CALL ts_bck_interp( jn, ll_init, za0, za1, za2, za3 ) ! coeficients of the interpolation
551	zsshp2_e(:,:) = za0 * ssha_e(:,:) + za1 * sshn_e (:,:) &
552	& + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:)
553	!
554	! ! Surface pressure gradient
555	zldg = ( 1._wp - rn_scal_load ) * grav ! local factor
556	DO jj = 2, jpjm1
557	DO ji = 2, jpim1
558	zu_spg(ji,jj) = - zldg * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj)
559	zv_spg(ji,jj) = - zldg * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj)
560	END DO
561	END DO
562	IF( ln_wd_il ) THEN ! W/D : gravity filters applied on pressure gradient
563	CALL wad_spg( zsshp2_e, zcpx, zcpy ) ! Calculating W/D gravity filters
564	zu_spg(2:jpim1,2:jpjm1) = zu_spg(2:jpim1,2:jpjm1) * zcpx(2:jpim1,2:jpjm1)
565	zv_spg(2:jpim1,2:jpjm1) = zv_spg(2:jpim1,2:jpjm1) * zcpy(2:jpim1,2:jpjm1)
566	ENDIF
567	!
568	! Add Coriolis trend:
569	! zwz array below or triads normally depend on sea level with ln_linssh=F and should be updated
570	! at each time step. We however keep them constant here for optimization.
571	! Recall that zhU and zhV hold fluxes at jn+0.5 (extrapolated not backward interpolated)
572	CALL dyn_cor_2d( zhup2_e, zhvp2_e, ua_e, va_e, zhU, zhV, zu_trd, zv_trd )
573	!
574	! Add tidal astronomical forcing if defined
575	IF ( ln_tide .AND. ln_tide_pot ) THEN
576	DO jj = 2, jpjm1
577	DO ji = fs_2, fs_jpim1 ! vector opt.
578	zu_trd(ji,jj) = zu_trd(ji,jj) + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) * r1_e1u(ji,jj)
579	zv_trd(ji,jj) = zv_trd(ji,jj) + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) * r1_e2v(ji,jj)
580	END DO
581	END DO
582	ENDIF
583	!
584	! Add bottom stresses:
585	!jth do implicitly instead
586	IF ( .NOT. ll_wd ) THEN ! Revert to explicit for bit comparison tests in non wad runs
587	DO jj = 2, jpjm1
588	DO ji = fs_2, fs_jpim1 ! vector opt.
589	zu_trd(ji,jj) = zu_trd(ji,jj) + zCdU_u(ji,jj) * un_e(ji,jj) * hur_e(ji,jj)
590	zv_trd(ji,jj) = zv_trd(ji,jj) + zCdU_v(ji,jj) * vn_e(ji,jj) * hvr_e(ji,jj)
591	END DO
592	END DO
593	ENDIF
594	!
595	! Set next velocities:
596	! Compute barotropic speeds at step jit+1 (h : total height of the water colomn)
597	!-- VECTOR FORM
598	!-- m+1 m / m+1/2 \ --!
599	!-- u = u + delta_t' * \ (1-r)g grad_x( ssh') - f * k vect u + frc / --!
600	!-- --!
601	!-- FLUX FORM --!
602	!-- m+1 __1__ / m m / m+1/2 m+1/2 m+1/2 n \ \ --!
603	!-- u = m+1 \| h * u + delta_t' * \ h * (1-r)g grad_x( ssh') - h * f * k vect u + h * frc / \| --!
604	!-- h \ / --!
605	!------------------------------------------------------------------------------------------------------------------------!
606	IF( ln_dynadv_vec .OR. ln_linssh ) THEN !* Vector form
607	DO jj = 2, jpjm1
608	DO ji = fs_2, fs_jpim1 ! vector opt.
609	ua_e(ji,jj) = ( un_e(ji,jj) &
610	& + rdtbt * ( zu_spg(ji,jj) &
611	& + zu_trd(ji,jj) &
612	& + zu_frc(ji,jj) ) &
613	& ) * ssumask(ji,jj)
614
615	va_e(ji,jj) = ( vn_e(ji,jj) &
616	& + rdtbt * ( zv_spg(ji,jj) &
617	& + zv_trd(ji,jj) &
618	& + zv_frc(ji,jj) ) &
619	& ) * ssvmask(ji,jj)
620	END DO
621	END DO
622	!
623	ELSE !* Flux form
624	DO jj = 2, jpjm1
625	DO ji = 2, jpim1
626	! ! hu_e, hv_e hold depth at jn, zhup2_e, zhvp2_e hold extrapolated depth at jn+1/2
627	! ! backward interpolated depth used in spg terms at jn+1/2
628	zhu_bck = hu_0(ji,jj) + r1_2r1_e1e2u(ji,jj) ( e1e2t(ji ,jj) * zsshp2_e(ji ,jj) &
629	& + e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj) ) * ssumask(ji,jj)
630	zhv_bck = hv_0(ji,jj) + r1_2r1_e1e2v(ji,jj) ( e1e2t(ji,jj ) * zsshp2_e(ji,jj ) &
631	& + e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1) ) * ssvmask(ji,jj)
632	! ! inverse depth at jn+1
633	z1_hu = ssumask(ji,jj) / ( hu_0(ji,jj) + zsshu_a(ji,jj) + 1._wp - ssumask(ji,jj) )
634	z1_hv = ssvmask(ji,jj) / ( hv_0(ji,jj) + zsshv_a(ji,jj) + 1._wp - ssvmask(ji,jj) )
635	!
636	ua_e(ji,jj) = ( hu_e (ji,jj) * un_e (ji,jj) &
637	& + rdtbt * ( zhu_bck * zu_spg (ji,jj) & !
638	& + zhup2_e(ji,jj) * zu_trd (ji,jj) & !
639	& + hu_n (ji,jj) * zu_frc (ji,jj) ) ) * z1_hu
640	!
641	va_e(ji,jj) = ( hv_e (ji,jj) * vn_e (ji,jj) &
642	& + rdtbt * ( zhv_bck * zv_spg (ji,jj) & !
643	& + zhvp2_e(ji,jj) * zv_trd (ji,jj) & !
644	& + hv_n (ji,jj) * zv_frc (ji,jj) ) ) * z1_hv
645	END DO
646	END DO
647	ENDIF
648	!jth implicit bottom friction:
649	IF ( ll_wd ) THEN ! revert to explicit for bit comparison tests in non wad runs
650	DO jj = 2, jpjm1
651	DO ji = fs_2, fs_jpim1 ! vector opt.
652	ua_e(ji,jj) = ua_e(ji,jj) /(1.0 - rdtbt * zCdU_u(ji,jj) * hur_e(ji,jj))
653	va_e(ji,jj) = va_e(ji,jj) /(1.0 - rdtbt * zCdU_v(ji,jj) * hvr_e(ji,jj))
654	END DO
655	END DO
656	ENDIF
657
658	IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only)
659	hu_e (2:jpim1,2:jpjm1) = hu_0(2:jpim1,2:jpjm1) + zsshu_a(2:jpim1,2:jpjm1)
660	hur_e(2:jpim1,2:jpjm1) = ssumask(2:jpim1,2:jpjm1) / ( hu_e(2:jpim1,2:jpjm1) + 1._wp - ssumask(2:jpim1,2:jpjm1) )
661	hv_e (2:jpim1,2:jpjm1) = hv_0(2:jpim1,2:jpjm1) + zsshv_a(2:jpim1,2:jpjm1)
662	hvr_e(2:jpim1,2:jpjm1) = ssvmask(2:jpim1,2:jpjm1) / ( hv_e(2:jpim1,2:jpjm1) + 1._wp - ssvmask(2:jpim1,2:jpjm1) )
663	CALL lbc_lnk_multi( 'dynspg_ts', ua_e , 'U', -1._wp, va_e , 'V', -1._wp &
664	& , hu_e , 'U', -1._wp, hv_e , 'V', -1._wp &
665	& , hur_e, 'U', -1._wp, hvr_e, 'V', -1._wp )
666	ELSE
667	CALL lbc_lnk_multi( 'dynspg_ts', ua_e , 'U', -1._wp, va_e , 'V', -1._wp )
668	ENDIF
669	!
670	!
671	! ! open boundaries
672	IF( ln_bdy ) CALL bdy_dyn2d( jn, ua_e, va_e, un_e, vn_e, hur_e, hvr_e, ssha_e )
673	#if defined key_agrif
674	IF( .NOT.Agrif_Root() ) CALL agrif_dyn_ts( jn ) ! Agrif
675	#endif
676	! !* Swap
677	! ! ----
678	ubb_e (:,:) = ub_e (:,:)
679	ub_e (:,:) = un_e (:,:)
680	un_e (:,:) = ua_e (:,:)
681	!
682	vbb_e (:,:) = vb_e (:,:)
683	vb_e (:,:) = vn_e (:,:)
684	vn_e (:,:) = va_e (:,:)
685	!
686	sshbb_e(:,:) = sshb_e(:,:)
687	sshb_e (:,:) = sshn_e(:,:)
688	sshn_e (:,:) = ssha_e(:,:)
689
690	! !* Sum over whole bt loop
691	! ! ----------------------
692	za1 = wgtbtp1(jn)
693	IF( ln_dynadv_vec .OR. ln_linssh ) THEN ! Sum velocities
694	ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:)
695	va_b (:,:) = va_b (:,:) + za1 * va_e (:,:)
696	ELSE ! Sum transports
697	IF ( .NOT.ln_wd_dl ) THEN
698	ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:) * hu_e (:,:)
699	va_b (:,:) = va_b (:,:) + za1 * va_e (:,:) * hv_e (:,:)
700	ELSE
701	ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:) * hu_e (:,:) * zuwdmask(:,:)
702	va_b (:,:) = va_b (:,:) + za1 * va_e (:,:) * hv_e (:,:) * zvwdmask(:,:)
703	END IF
704	ENDIF
705	! ! Sum sea level
706	ssha(:,:) = ssha(:,:) + za1 * ssha_e(:,:)
707
708	! ! ==================== !
709	END DO ! end loop !
710	! ! ==================== !
711	! -----------------------------------------------------------------------------
712	! Phase 3. update the general trend with the barotropic trend
713	! -----------------------------------------------------------------------------
714	!
715	! Set advection velocity correction:
716	IF (ln_bt_fw) THEN
717	IF( .NOT.( kt == nit000 .AND. neuler==0 ) ) THEN
718	DO jj = 1, jpj
719	DO ji = 1, jpi
720	zun_save = un_adv(ji,jj)
721	zvn_save = vn_adv(ji,jj)
722	! ! apply the previously computed correction
723	un_adv(ji,jj) = r1_2 * ( ub2_b(ji,jj) + zun_save - atfp * un_bf(ji,jj) )
724	vn_adv(ji,jj) = r1_2 * ( vb2_b(ji,jj) + zvn_save - atfp * vn_bf(ji,jj) )
725	! ! Update corrective fluxes for next time step
726	un_bf(ji,jj) = atfp * un_bf(ji,jj) + ( zun_save - ub2_b(ji,jj) )
727	vn_bf(ji,jj) = atfp * vn_bf(ji,jj) + ( zvn_save - vb2_b(ji,jj) )
728	! ! Save integrated transport for next computation
729	ub2_b(ji,jj) = zun_save
730	vb2_b(ji,jj) = zvn_save
731	END DO
732	END DO
733	ELSE
734	un_bf(:,:) = 0._wp ! corrective fluxes for next time step set to zero
735	vn_bf(:,:) = 0._wp
736	ub2_b(:,:) = un_adv(:,:) ! Save integrated transport for next computation
737	vb2_b(:,:) = vn_adv(:,:)
738	END IF
739	ENDIF
740
741
742	!
743	! Update barotropic trend:
744	IF( ln_dynadv_vec .OR. ln_linssh ) THEN
745	DO jk=1,jpkm1
746	ua(:,:,jk) = ua(:,:,jk) + ( ua_b(:,:) - ub_b(:,:) ) * r1_2dt_b
747	va(:,:,jk) = va(:,:,jk) + ( va_b(:,:) - vb_b(:,:) ) * r1_2dt_b
748	END DO
749	ELSE
750	! At this stage, ssha has been corrected: compute new depths at velocity points
751	DO jj = 1, jpjm1
752	DO ji = 1, jpim1 ! NO Vector Opt.
753	zsshu_a(ji,jj) = r1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj) &
754	& * ( e1e2t(ji ,jj) * ssha(ji ,jj) &
755	& + e1e2t(ji+1,jj) * ssha(ji+1,jj) )
756	zsshv_a(ji,jj) = r1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj) &
757	& * ( e1e2t(ji,jj ) * ssha(ji,jj ) &
758	& + e1e2t(ji,jj+1) * ssha(ji,jj+1) )
759	END DO
760	END DO
761	CALL lbc_lnk_multi( 'dynspg_ts', zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) ! Boundary conditions
762	!
763	DO jk=1,jpkm1
764	ua(:,:,jk) = ua(:,:,jk) + r1_hu_n(:,:) * ( ua_b(:,:) - ub_b(:,:) * hu_b(:,:) ) * r1_2dt_b
765	va(:,:,jk) = va(:,:,jk) + r1_hv_n(:,:) * ( va_b(:,:) - vb_b(:,:) * hv_b(:,:) ) * r1_2dt_b
766	END DO
767	! Save barotropic velocities not transport:
768	ua_b(:,:) = ua_b(:,:) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - ssumask(:,:) )
769	va_b(:,:) = va_b(:,:) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - ssvmask(:,:) )
770	ENDIF
771
772
773	! Correct velocities so that the barotropic velocity equals (un_adv, vn_adv) (in all cases)
774	DO jk = 1, jpkm1
775	un(:,:,jk) = ( un(:,:,jk) + un_adv(:,:)r1_hu_n(:,:) - un_b(:,:) ) umask(:,:,jk)
776	vn(:,:,jk) = ( vn(:,:,jk) + vn_adv(:,:)r1_hv_n(:,:) - vn_b(:,:) ) vmask(:,:,jk)
777	END DO
778
779	IF ( ln_wd_dl .and. ln_wd_dl_bc) THEN
780	DO jk = 1, jpkm1
781	un(:,:,jk) = ( un_adv(:,:)*r1_hu_n(:,:) &
782	& + zuwdav2(:,:)(un(:,:,jk) - un_adv(:,:)r1_hu_n(:,:)) ) * umask(:,:,jk)
783	vn(:,:,jk) = ( vn_adv(:,:)*r1_hv_n(:,:) &
784	& + zvwdav2(:,:)(vn(:,:,jk) - vn_adv(:,:)r1_hv_n(:,:)) ) * vmask(:,:,jk)
785	END DO
786	END IF
787
788
789	CALL iom_put( "ubar", un_adv(:,:)*r1_hu_n(:,:) ) ! barotropic i-current
790	CALL iom_put( "vbar", vn_adv(:,:)*r1_hv_n(:,:) ) ! barotropic i-current
791	!
792	#if defined key_agrif
793	! Save time integrated fluxes during child grid integration
794	! (used to update coarse grid transports at next time step)
795	!
796	IF( .NOT.Agrif_Root() .AND. ln_bt_fw ) THEN
797	IF( Agrif_NbStepint() == 0 ) THEN
798	ub2_i_b(:,:) = 0._wp
799	vb2_i_b(:,:) = 0._wp
800	END IF
801	!
802	za1 = 1._wp / REAL(Agrif_rhot(), wp)
803	ub2_i_b(:,:) = ub2_i_b(:,:) + za1 * ub2_b(:,:)
804	vb2_i_b(:,:) = vb2_i_b(:,:) + za1 * vb2_b(:,:)
805	ENDIF
806	#endif
807	! !* write time-spliting arrays in the restart
808	IF( lrst_oce .AND.ln_bt_fw ) CALL ts_rst( kt, 'WRITE' )
809	!
810	IF( ln_wd_il ) DEALLOCATE( zcpx, zcpy )
811	IF( ln_wd_dl ) DEALLOCATE( ztwdmask, zuwdmask, zvwdmask, zuwdav2, zvwdav2 )
812	!
813	IF( ln_diatmb ) THEN
814	CALL iom_put( "baro_u" , un_bssumask(:,:)+zmdi(1.-ssumask(:,:) ) ) ! Barotropic U Velocity
815	CALL iom_put( "baro_v" , vn_bssvmask(:,:)+zmdi(1.-ssvmask(:,:) ) ) ! Barotropic V Velocity
816	ENDIF
817	!
818	END SUBROUTINE dyn_spg_ts
819
820
821	SUBROUTINE ts_wgt( ll_av, ll_fw, jpit, zwgt1, zwgt2)
822	!!---------------------------------------------------------------------
823	!! * ROUTINE ts_wgt *
824	!!
825	!! ** Purpose : Set time-splitting weights for temporal averaging (or not)
826	!!----------------------------------------------------------------------
827	LOGICAL, INTENT(in) :: ll_av ! temporal averaging=.true.
828	LOGICAL, INTENT(in) :: ll_fw ! forward time splitting =.true.
829	INTEGER, INTENT(inout) :: jpit ! cycle length
830	REAL(wp), DIMENSION(3*nn_baro), INTENT(inout) :: zwgt1, & ! Primary weights
831	zwgt2 ! Secondary weights
832
833	INTEGER :: jic, jn, ji ! temporary integers
834	REAL(wp) :: za1, za2
835	!!----------------------------------------------------------------------
836
837	zwgt1(:) = 0._wp
838	zwgt2(:) = 0._wp
839
840	! Set time index when averaged value is requested
841	IF (ll_fw) THEN
842	jic = nn_baro
843	ELSE
844	jic = 2 * nn_baro
845	ENDIF
846
847	! Set primary weights:
848	IF (ll_av) THEN
849	! Define simple boxcar window for primary weights
850	! (width = nn_baro, centered around jic)
851	SELECT CASE ( nn_bt_flt )
852	CASE( 0 ) ! No averaging
853	zwgt1(jic) = 1._wp
854	jpit = jic
855
856	CASE( 1 ) ! Boxcar, width = nn_baro
857	DO jn = 1, 3*nn_baro
858	za1 = ABS(float(jn-jic))/float(nn_baro)
859	IF (za1 < 0.5_wp) THEN
860	zwgt1(jn) = 1._wp
861	jpit = jn
862	ENDIF
863	ENDDO
864
865	CASE( 2 ) ! Boxcar, width = 2 * nn_baro
866	DO jn = 1, 3*nn_baro
867	za1 = ABS(float(jn-jic))/float(nn_baro)
868	IF (za1 < 1._wp) THEN
869	zwgt1(jn) = 1._wp
870	jpit = jn
871	ENDIF
872	ENDDO
873	CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_bt_flt' )
874	END SELECT
875
876	ELSE ! No time averaging
877	zwgt1(jic) = 1._wp
878	jpit = jic
879	ENDIF
880
881	! Set secondary weights
882	DO jn = 1, jpit
883	DO ji = jn, jpit
884	zwgt2(jn) = zwgt2(jn) + zwgt1(ji)
885	END DO
886	END DO
887
888	! Normalize weigths:
889	za1 = 1._wp / SUM(zwgt1(1:jpit))
890	za2 = 1._wp / SUM(zwgt2(1:jpit))
891	DO jn = 1, jpit
892	zwgt1(jn) = zwgt1(jn) * za1
893	zwgt2(jn) = zwgt2(jn) * za2
894	END DO
895	!
896	END SUBROUTINE ts_wgt
897
898
899	SUBROUTINE ts_rst( kt, cdrw )
900	!!---------------------------------------------------------------------
901	!! * ROUTINE ts_rst *
902	!!
903	!! ** Purpose : Read or write time-splitting arrays in restart file
904	!!----------------------------------------------------------------------
905	INTEGER , INTENT(in) :: kt ! ocean time-step
906	CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag
907	!!----------------------------------------------------------------------
908	!
909	IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialise
910	! ! ---------------
911	IF( ln_rstart .AND. ln_bt_fw .AND. (neuler/=0) ) THEN !* Read the restart file
912	CALL iom_get( numror, jpdom_autoglo, 'ub2_b' , ub2_b (:,:), ldxios = lrxios )
913	CALL iom_get( numror, jpdom_autoglo, 'vb2_b' , vb2_b (:,:), ldxios = lrxios )
914	CALL iom_get( numror, jpdom_autoglo, 'un_bf' , un_bf (:,:), ldxios = lrxios )
915	CALL iom_get( numror, jpdom_autoglo, 'vn_bf' , vn_bf (:,:), ldxios = lrxios )
916	IF( .NOT.ln_bt_av ) THEN
917	CALL iom_get( numror, jpdom_autoglo, 'sshbb_e' , sshbb_e(:,:), ldxios = lrxios )
918	CALL iom_get( numror, jpdom_autoglo, 'ubb_e' , ubb_e(:,:), ldxios = lrxios )
919	CALL iom_get( numror, jpdom_autoglo, 'vbb_e' , vbb_e(:,:), ldxios = lrxios )
920	CALL iom_get( numror, jpdom_autoglo, 'sshb_e' , sshb_e(:,:), ldxios = lrxios )
921	CALL iom_get( numror, jpdom_autoglo, 'ub_e' , ub_e(:,:), ldxios = lrxios )
922	CALL iom_get( numror, jpdom_autoglo, 'vb_e' , vb_e(:,:), ldxios = lrxios )
923	ENDIF
924	#if defined key_agrif
925	! Read time integrated fluxes
926	IF ( .NOT.Agrif_Root() ) THEN
927	CALL iom_get( numror, jpdom_autoglo, 'ub2_i_b' , ub2_i_b(:,:), ldxios = lrxios )
928	CALL iom_get( numror, jpdom_autoglo, 'vb2_i_b' , vb2_i_b(:,:), ldxios = lrxios )
929	ENDIF
930	#endif
931	ELSE !* Start from rest
932	IF(lwp) WRITE(numout,*)
933	IF(lwp) WRITE(numout,*) ' ==>>> start from rest: set barotropic values to 0'
934	ub2_b (:,:) = 0._wp ; vb2_b (:,:) = 0._wp ! used in the 1st interpol of agrif
935	un_adv(:,:) = 0._wp ; vn_adv(:,:) = 0._wp ! used in the 1st interpol of agrif
936	un_bf (:,:) = 0._wp ; vn_bf (:,:) = 0._wp ! used in the 1st update of agrif
937	#if defined key_agrif
938	IF ( .NOT.Agrif_Root() ) THEN
939	ub2_i_b(:,:) = 0._wp ; vb2_i_b(:,:) = 0._wp ! used in the 1st update of agrif
940	ENDIF
941	#endif
942	ENDIF
943	!
944	ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN ! Create restart file
945	! ! -------------------
946	IF(lwp) WRITE(numout,*) '---- ts_rst ----'
947	IF( lwxios ) CALL iom_swap( cwxios_context )
948	CALL iom_rstput( kt, nitrst, numrow, 'ub2_b' , ub2_b (:,:), ldxios = lwxios )
949	CALL iom_rstput( kt, nitrst, numrow, 'vb2_b' , vb2_b (:,:), ldxios = lwxios )
950	CALL iom_rstput( kt, nitrst, numrow, 'un_bf' , un_bf (:,:), ldxios = lwxios )
951	CALL iom_rstput( kt, nitrst, numrow, 'vn_bf' , vn_bf (:,:), ldxios = lwxios )
952	!
953	IF (.NOT.ln_bt_av) THEN
954	CALL iom_rstput( kt, nitrst, numrow, 'sshbb_e' , sshbb_e(:,:), ldxios = lwxios )
955	CALL iom_rstput( kt, nitrst, numrow, 'ubb_e' , ubb_e(:,:), ldxios = lwxios )
956	CALL iom_rstput( kt, nitrst, numrow, 'vbb_e' , vbb_e(:,:), ldxios = lwxios )
957	CALL iom_rstput( kt, nitrst, numrow, 'sshb_e' , sshb_e(:,:), ldxios = lwxios )
958	CALL iom_rstput( kt, nitrst, numrow, 'ub_e' , ub_e(:,:), ldxios = lwxios )
959	CALL iom_rstput( kt, nitrst, numrow, 'vb_e' , vb_e(:,:), ldxios = lwxios )
960	ENDIF
961	#if defined key_agrif
962	! Save time integrated fluxes
963	IF ( .NOT.Agrif_Root() ) THEN
964	CALL iom_rstput( kt, nitrst, numrow, 'ub2_i_b' , ub2_i_b(:,:), ldxios = lwxios )
965	CALL iom_rstput( kt, nitrst, numrow, 'vb2_i_b' , vb2_i_b(:,:), ldxios = lwxios )
966	ENDIF
967	#endif
968	IF( lwxios ) CALL iom_swap( cxios_context )
969	ENDIF
970	!
971	END SUBROUTINE ts_rst
972
973
974	SUBROUTINE dyn_spg_ts_init
975	!!---------------------------------------------------------------------
976	!! * ROUTINE dyn_spg_ts_init *
977	!!
978	!! ** Purpose : Set time splitting options
979	!!----------------------------------------------------------------------
980	INTEGER :: ji ,jj ! dummy loop indices
981	REAL(wp) :: zxr2, zyr2, zcmax ! local scalar
982	REAL(wp), DIMENSION(jpi,jpj) :: zcu
983	INTEGER :: inum
984	!!----------------------------------------------------------------------
985	!
986	! Max courant number for ext. grav. waves
987	!
988	DO jj = 1, jpj
989	DO ji =1, jpi
990	zxr2 = r1_e1t(ji,jj) * r1_e1t(ji,jj)
991	zyr2 = r1_e2t(ji,jj) * r1_e2t(ji,jj)
992	zcu(ji,jj) = SQRT( grav * MAX(ht_0(ji,jj),0._wp) * (zxr2 + zyr2) )
993	END DO
994	END DO
995	!
996	zcmax = MAXVAL( zcu(:,:) )
997	CALL mpp_max( 'dynspg_ts', zcmax )
998
999	! Estimate number of iterations to satisfy a max courant number= rn_bt_cmax
1000	IF( ln_bt_auto ) nn_baro = CEILING( rdt / rn_bt_cmax * zcmax)
1001
1002	rdtbt = rdt / REAL( nn_baro , wp )
1003	zcmax = zcmax * rdtbt
1004	! Print results
1005	IF(lwp) WRITE(numout,*)
1006	IF(lwp) WRITE(numout,*) 'dyn_spg_ts_init : split-explicit free surface'
1007	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~'
1008	IF( ln_bt_auto ) THEN
1009	IF(lwp) WRITE(numout,*) ' ln_ts_auto =.true. Automatically set nn_baro '
1010	IF(lwp) WRITE(numout,*) ' Max. courant number allowed: ', rn_bt_cmax
1011	ELSE
1012	IF(lwp) WRITE(numout,*) ' ln_ts_auto=.false.: Use nn_baro in namelist nn_baro = ', nn_baro
1013	ENDIF
1014
1015	IF(ln_bt_av) THEN
1016	IF(lwp) WRITE(numout,*) ' ln_bt_av =.true. ==> Time averaging over nn_baro time steps is on '
1017	ELSE
1018	IF(lwp) WRITE(numout,*) ' ln_bt_av =.false. => No time averaging of barotropic variables '
1019	ENDIF
1020	!
1021	!
1022	IF(ln_bt_fw) THEN
1023	IF(lwp) WRITE(numout,*) ' ln_bt_fw=.true. => Forward integration of barotropic variables '
1024	ELSE
1025	IF(lwp) WRITE(numout,*) ' ln_bt_fw =.false.=> Centred integration of barotropic variables '
1026	ENDIF
1027	!
1028	#if defined key_agrif
1029	! Restrict the use of Agrif to the forward case only
1030	!!! IF( .NOT.ln_bt_fw .AND. .NOT.Agrif_Root() ) CALL ctl_stop( 'AGRIF not implemented if ln_bt_fw=.FALSE.' )
1031	#endif
1032	!
1033	IF(lwp) WRITE(numout,*) ' Time filter choice, nn_bt_flt: ', nn_bt_flt
1034	SELECT CASE ( nn_bt_flt )
1035	CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' Dirac'
1036	CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' Boxcar: width = nn_baro'
1037	CASE( 2 ) ; IF(lwp) WRITE(numout,) ' Boxcar: width = 2nn_baro'
1038	CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_bt_flt: should 0,1, or 2' )
1039	END SELECT
1040	!
1041	IF(lwp) WRITE(numout,*) ' '
1042	IF(lwp) WRITE(numout,*) ' nn_baro = ', nn_baro
1043	IF(lwp) WRITE(numout,*) ' Barotropic time step [s] is :', rdtbt
1044	IF(lwp) WRITE(numout,*) ' Maximum Courant number is :', zcmax
1045	!
1046	IF(lwp) WRITE(numout,*) ' Time diffusion parameter rn_bt_alpha: ', rn_bt_alpha
1047	IF ((ln_bt_av.AND.nn_bt_flt/=0).AND.(rn_bt_alpha>0._wp)) THEN
1048	CALL ctl_stop( 'dynspg_ts ERROR: if rn_bt_alpha > 0, remove temporal averaging' )
1049	ENDIF
1050	!
1051	IF( .NOT.ln_bt_av .AND. .NOT.ln_bt_fw ) THEN
1052	CALL ctl_stop( 'dynspg_ts ERROR: No time averaging => only forward integration is possible' )
1053	ENDIF
1054	IF( zcmax>0.9_wp ) THEN
1055	CALL ctl_stop( 'dynspg_ts ERROR: Maximum Courant number is greater than 0.9: Inc. nn_baro !' )
1056	ENDIF
1057	!
1058	! ! Allocate time-splitting arrays
1059	IF( dyn_spg_ts_alloc() /= 0 ) CALL ctl_stop('STOP', 'dyn_spg_init: failed to allocate dynspg_ts arrays' )
1060	!
1061	! ! read restart when needed
1062	CALL ts_rst( nit000, 'READ' )
1063	!
1064	IF( lwxios ) THEN
1065	! define variables in restart file when writing with XIOS
1066	CALL iom_set_rstw_var_active('ub2_b')
1067	CALL iom_set_rstw_var_active('vb2_b')
1068	CALL iom_set_rstw_var_active('un_bf')
1069	CALL iom_set_rstw_var_active('vn_bf')
1070	!
1071	IF (.NOT.ln_bt_av) THEN
1072	CALL iom_set_rstw_var_active('sshbb_e')
1073	CALL iom_set_rstw_var_active('ubb_e')
1074	CALL iom_set_rstw_var_active('vbb_e')
1075	CALL iom_set_rstw_var_active('sshb_e')
1076	CALL iom_set_rstw_var_active('ub_e')
1077	CALL iom_set_rstw_var_active('vb_e')
1078	ENDIF
1079	#if defined key_agrif
1080	! Save time integrated fluxes
1081	IF ( .NOT.Agrif_Root() ) THEN
1082	CALL iom_set_rstw_var_active('ub2_i_b')
1083	CALL iom_set_rstw_var_active('vb2_i_b')
1084	ENDIF
1085	#endif
1086	ENDIF
1087	!
1088	END SUBROUTINE dyn_spg_ts_init
1089
1090
1091	SUBROUTINE dyn_cor_2d_init
1092	!!---------------------------------------------------------------------
1093	!! * ROUTINE dyn_cor_2d_init *
1094	!!
1095	!! ** Purpose : Set time splitting options
1096	!! Set arrays to remove/compute coriolis trend.
1097	!! Do it once during initialization if volume is fixed, else at each long time step.
1098	!! Note that these arrays are also used during barotropic loop. These are however frozen
1099	!! although they should be updated in the variable volume case. Not a big approximation.
1100	!! To remove this approximation, copy lines below inside barotropic loop
1101	!! and update depths at T-F points (ht and zhf resp.) at each barotropic time step
1102	!!
1103	!! Compute zwz = f / ( height of the water colomn )
1104	!!----------------------------------------------------------------------
1105	INTEGER :: ji ,jj, jk ! dummy loop indices
1106	REAL(wp) :: z1_ht
1107	REAL(wp), DIMENSION(jpi,jpj) :: zhf
1108	!!----------------------------------------------------------------------
1109	!
1110	SELECT CASE( nvor_scheme )
1111	CASE( np_EEN ) != EEN scheme using e3f (energy & enstrophy scheme)
1112	SELECT CASE( nn_een_e3f ) !* ff_f/e3 at F-point
1113	CASE ( 0 ) ! original formulation (masked averaging of e3t divided by 4)
1114	DO jj = 1, jpjm1
1115	DO ji = 1, jpim1
1116	zwz(ji,jj) = ( ht_n(ji ,jj+1) + ht_n(ji+1,jj+1) + &
1117	& ht_n(ji ,jj ) + ht_n(ji+1,jj ) ) * 0.25_wp
1118	IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj)
1119	END DO
1120	END DO
1121	CASE ( 1 ) ! new formulation (masked averaging of e3t divided by the sum of mask)
1122	DO jj = 1, jpjm1
1123	DO ji = 1, jpim1
1124	zwz(ji,jj) = ( ht_n (ji ,jj+1) + ht_n (ji+1,jj+1) &
1125	& + ht_n (ji ,jj ) + ht_n (ji+1,jj ) ) &
1126	& / ( MAX( 1._wp, ssmask(ji ,jj+1) + ssmask(ji+1,jj+1) &
1127	& + ssmask(ji ,jj ) + ssmask(ji+1,jj ) ) )
1128	IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj)
1129	END DO
1130	END DO
1131	END SELECT
1132	CALL lbc_lnk( 'dynspg_ts', zwz, 'F', 1._wp )
1133	!
1134	ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp
1135	DO jj = 2, jpj
1136	DO ji = 2, jpi
1137	ftne(ji,jj) = zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1)
1138	ftnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj ) + zwz(ji ,jj )
1139	ftse(ji,jj) = zwz(ji ,jj ) + zwz(ji ,jj-1) + zwz(ji-1,jj-1)
1140	ftsw(ji,jj) = zwz(ji ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj )
1141	END DO
1142	END DO
1143	!
1144	CASE( np_EET ) != EEN scheme using e3t (energy conserving scheme)
1145	ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp
1146	DO jj = 2, jpj
1147	DO ji = 2, jpi
1148	z1_ht = ssmask(ji,jj) / ( ht_n(ji,jj) + 1._wp - ssmask(ji,jj) )
1149	ftne(ji,jj) = ( ff_f(ji-1,jj ) + ff_f(ji ,jj ) + ff_f(ji ,jj-1) ) * z1_ht
1150	ftnw(ji,jj) = ( ff_f(ji-1,jj-1) + ff_f(ji-1,jj ) + ff_f(ji ,jj ) ) * z1_ht
1151	ftse(ji,jj) = ( ff_f(ji ,jj ) + ff_f(ji ,jj-1) + ff_f(ji-1,jj-1) ) * z1_ht
1152	ftsw(ji,jj) = ( ff_f(ji ,jj-1) + ff_f(ji-1,jj-1) + ff_f(ji-1,jj ) ) * z1_ht
1153	END DO
1154	END DO
1155	!
1156	CASE( np_ENE, np_ENS , np_MIX ) != all other schemes (ENE, ENS, MIX) except ENT !
1157	!
1158	zwz(:,:) = 0._wp
1159	zhf(:,:) = 0._wp
1160
1161	!!gm assume 0 in both cases (which is almost surely WRONG ! ) as hvatf has been removed
1162	!!gm A priori a better value should be something like :
1163	!!gm zhf(i,j) = masked sum of ht(i,j) , ht(i+1,j) , ht(i,j+1) , (i+1,j+1)
1164	!!gm divided by the sum of the corresponding mask
1165	!!gm
1166	!!
1167	IF( .NOT.ln_sco ) THEN
1168
1169	!!gm agree the JC comment : this should be done in a much clear way
1170
1171	! JC: It not clear yet what should be the depth at f-points over land in z-coordinate case
1172	! Set it to zero for the time being
1173	! IF( rn_hmin < 0._wp ) THEN ; jk = - INT( rn_hmin ) ! from a nb of level
1174	! ELSE ; jk = MINLOC( gdepw_0, mask = gdepw_0 > rn_hmin, dim = 1 ) ! from a depth
1175	! ENDIF
1176	! zhf(:,:) = gdepw_0(:,:,jk+1)
1177	!
1178	ELSE
1179	!
1180	!zhf(:,:) = hbatf(:,:)
1181	DO jj = 1, jpjm1
1182	DO ji = 1, jpim1
1183	zhf(ji,jj) = ( ht_0 (ji,jj ) + ht_0 (ji+1,jj ) &
1184	& + ht_0 (ji,jj+1) + ht_0 (ji+1,jj+1) ) &
1185	& / MAX( ssmask(ji,jj ) + ssmask(ji+1,jj ) &
1186	& + ssmask(ji,jj+1) + ssmask(ji+1,jj+1) , 1._wp )
1187	END DO
1188	END DO
1189	ENDIF
1190	!
1191	DO jj = 1, jpjm1
1192	zhf(:,jj) = zhf(:,jj) * (1._wp- umask(:,jj,1) * umask(:,jj+1,1))
1193	END DO
1194	!
1195	DO jk = 1, jpkm1
1196	DO jj = 1, jpjm1
1197	zhf(:,jj) = zhf(:,jj) + e3f_n(:,jj,jk) * umask(:,jj,jk) * umask(:,jj+1,jk)
1198	END DO
1199	END DO
1200	CALL lbc_lnk( 'dynspg_ts', zhf, 'F', 1._wp )
1201	! JC: TBC. hf should be greater than 0
1202	DO jj = 1, jpj
1203	DO ji = 1, jpi
1204	IF( zhf(ji,jj) /= 0._wp ) zwz(ji,jj) = 1._wp / zhf(ji,jj)
1205	END DO
1206	END DO
1207	zwz(:,:) = ff_f(:,:) * zwz(:,:)
1208	END SELECT
1209
1210	END SUBROUTINE dyn_cor_2d_init
1211
1212
1213
1214	SUBROUTINE dyn_cor_2d( hu_n, hv_n, un_b, vn_b, zhU, zhV, zu_trd, zv_trd )
1215	!!---------------------------------------------------------------------
1216	!! * ROUTINE dyn_cor_2d *
1217	!!
1218	!! ** Purpose : Compute u and v coriolis trends
1219	!!----------------------------------------------------------------------
1220	INTEGER :: ji ,jj ! dummy loop indices
1221	REAL(wp) :: zx1, zx2, zy1, zy2, z1_hu, z1_hv ! - -
1222	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: hu_n, hv_n, un_b, vn_b, zhU, zhV
1223	REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: zu_trd, zv_trd
1224	!!----------------------------------------------------------------------
1225	SELECT CASE( nvor_scheme )
1226	CASE( np_ENT ) ! enstrophy conserving scheme (f-point)
1227	DO jj = 2, jpjm1
1228	DO ji = 2, jpim1
1229	z1_hu = ssumask(ji,jj) / ( hu_n(ji,jj) + 1._wp - ssumask(ji,jj) )
1230	z1_hv = ssvmask(ji,jj) / ( hv_n(ji,jj) + 1._wp - ssvmask(ji,jj) )
1231	zu_trd(ji,jj) = + r1_4 * r1_e1e2u(ji,jj) * z1_hu &
1232	& * ( 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) ) &
1233	& + e1e2t(ji ,jj)ht_n(ji ,jj)ff_t(ji ,jj) * ( vn_b(ji ,jj) + vn_b(ji ,jj-1) ) )
1234	!
1235	zv_trd(ji,jj) = - r1_4 * r1_e1e2v(ji,jj) * z1_hv &
1236	& * ( 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) ) &
1237	& + e1e2t(ji,jj )ht_n(ji,jj )ff_t(ji,jj ) * ( un_b(ji,jj ) + un_b(ji-1,jj ) ) )
1238	END DO
1239	END DO
1240	!
1241	CASE( np_ENE , np_MIX ) ! energy conserving scheme (t-point) ENE or MIX
1242	DO jj = 2, jpjm1
1243	DO ji = fs_2, fs_jpim1 ! vector opt.
1244	zy1 = ( zhV(ji,jj-1) + zhV(ji+1,jj-1) ) * r1_e1u(ji,jj)
1245	zy2 = ( zhV(ji,jj ) + zhV(ji+1,jj ) ) * r1_e1u(ji,jj)
1246	zx1 = ( zhU(ji-1,jj) + zhU(ji-1,jj+1) ) * r1_e2v(ji,jj)
1247	zx2 = ( zhU(ji ,jj) + zhU(ji ,jj+1) ) * r1_e2v(ji,jj)
1248	! energy conserving formulation for planetary vorticity term
1249	zu_trd(ji,jj) = r1_4 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
1250	zv_trd(ji,jj) = - r1_4 * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 )
1251	END DO
1252	END DO
1253	!
1254	CASE( np_ENS ) ! enstrophy conserving scheme (f-point)
1255	DO jj = 2, jpjm1
1256	DO ji = fs_2, fs_jpim1 ! vector opt.
1257	zy1 = r1_8 * ( zhV(ji ,jj-1) + zhV(ji+1,jj-1) &
1258	& + zhV(ji ,jj ) + zhV(ji+1,jj ) ) * r1_e1u(ji,jj)
1259	zx1 = - r1_8 * ( zhU(ji-1,jj ) + zhU(ji-1,jj+1) &
1260	& + zhU(ji ,jj ) + zhU(ji ,jj+1) ) * r1_e2v(ji,jj)
1261	zu_trd(ji,jj) = zy1 * ( zwz(ji ,jj-1) + zwz(ji,jj) )
1262	zv_trd(ji,jj) = zx1 * ( zwz(ji-1,jj ) + zwz(ji,jj) )
1263	END DO
1264	END DO
1265	!
1266	CASE( np_EET , np_EEN ) ! energy & enstrophy scheme (using e3t or e3f)
1267	DO jj = 2, jpjm1
1268	DO ji = fs_2, fs_jpim1 ! vector opt.
1269	zu_trd(ji,jj) = + r1_12 * r1_e1u(ji,jj) * ( ftne(ji,jj ) * zhV(ji ,jj ) &
1270	& + ftnw(ji+1,jj) * zhV(ji+1,jj ) &
1271	& + ftse(ji,jj ) * zhV(ji ,jj-1) &
1272	& + ftsw(ji+1,jj) * zhV(ji+1,jj-1) )
1273	zv_trd(ji,jj) = - r1_12 * r1_e2v(ji,jj) * ( ftsw(ji,jj+1) * zhU(ji-1,jj+1) &
1274	& + ftse(ji,jj+1) * zhU(ji ,jj+1) &
1275	& + ftnw(ji,jj ) * zhU(ji-1,jj ) &
1276	& + ftne(ji,jj ) * zhU(ji ,jj ) )
1277	END DO
1278	END DO
1279	!
1280	END SELECT
1281	!
1282	END SUBROUTINE dyn_cor_2D
1283
1284
1285	SUBROUTINE wad_tmsk( pssh, ptmsk )
1286	!!----------------------------------------------------------------------
1287	!! * ROUTINE wad_lmt *
1288	!!
1289	!! ** Purpose : set wetting & drying mask at tracer points
1290	!! for the current barotropic sub-step
1291	!!
1292	!! ** Method : ???
1293	!!
1294	!! ** Action : ptmsk : wetting & drying t-mask
1295	!!----------------------------------------------------------------------
1296	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pssh !
1297	REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: ptmsk !
1298	!
1299	INTEGER :: ji, jj ! dummy loop indices
1300	!!----------------------------------------------------------------------
1301	!
1302	IF( ln_wd_dl_rmp ) THEN
1303	DO jj = 1, jpj
1304	DO ji = 1, jpi
1305	IF ( pssh(ji,jj) + ht_0(ji,jj) > 2._wp * rn_wdmin1 ) THEN
1306	! IF ( pssh(ji,jj) + ht_0(ji,jj) > rn_wdmin2 ) THEN
1307	ptmsk(ji,jj) = 1._wp
1308	ELSEIF( pssh(ji,jj) + ht_0(ji,jj) > rn_wdmin1 ) THEN
1309	ptmsk(ji,jj) = TANH( 50._wp( ( pssh(ji,jj) + ht_0(ji,jj) - rn_wdmin1 )r_rn_wdmin1) )
1310	ELSE
1311	ptmsk(ji,jj) = 0._wp
1312	ENDIF
1313	END DO
1314	END DO
1315	ELSE
1316	DO jj = 1, jpj
1317	DO ji = 1, jpi
1318	IF ( pssh(ji,jj) + ht_0(ji,jj) > rn_wdmin1 ) THEN ; ptmsk(ji,jj) = 1._wp
1319	ELSE ; ptmsk(ji,jj) = 0._wp
1320	ENDIF
1321	END DO
1322	END DO
1323	ENDIF
1324	!
1325	END SUBROUTINE wad_tmsk
1326
1327
1328	SUBROUTINE wad_Umsk( pTmsk, phU, phV, pu, pv, pUmsk, pVmsk )
1329	!!----------------------------------------------------------------------
1330	!! * ROUTINE wad_lmt *
1331	!!
1332	!! ** Purpose : set wetting & drying mask at tracer points
1333	!! for the current barotropic sub-step
1334	!!
1335	!! ** Method : ???
1336	!!
1337	!! ** Action : ptmsk : wetting & drying t-mask
1338	!!----------------------------------------------------------------------
1339	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pTmsk ! W & D t-mask
1340	REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: phU, phV, pu, pv ! ocean velocities and transports
1341	REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pUmsk, pVmsk ! W & D u- and v-mask
1342	!
1343	INTEGER :: ji, jj ! dummy loop indices
1344	!!----------------------------------------------------------------------
1345	!
1346	DO jj = 1, jpj
1347	DO ji = 1, jpim1 ! not jpi-column
1348	IF ( phU(ji,jj) > 0._wp ) THEN ; pUmsk(ji,jj) = pTmsk(ji ,jj)
1349	ELSE ; pUmsk(ji,jj) = pTmsk(ji+1,jj)
1350	ENDIF
1351	phU(ji,jj) = pUmsk(ji,jj)*phU(ji,jj)
1352	pu (ji,jj) = pUmsk(ji,jj)*pu (ji,jj)
1353	END DO
1354	END DO
1355	!
1356	DO jj = 1, jpjm1 ! not jpj-row
1357	DO ji = 1, jpi
1358	IF ( phV(ji,jj) > 0._wp ) THEN ; pVmsk(ji,jj) = pTmsk(ji,jj )
1359	ELSE ; pVmsk(ji,jj) = pTmsk(ji,jj+1)
1360	ENDIF
1361	phV(ji,jj) = pVmsk(ji,jj)*phV(ji,jj)
1362	pv (ji,jj) = pVmsk(ji,jj)*pv (ji,jj)
1363	END DO
1364	END DO
1365	!
1366	END SUBROUTINE wad_Umsk
1367
1368
1369	SUBROUTINE wad_spg( sshn, zcpx, zcpy )
1370	!!---------------------------------------------------------------------
1371	!! * ROUTINE wad_sp *
1372	!!
1373	!! ** Purpose :
1374	!!----------------------------------------------------------------------
1375	INTEGER :: ji ,jj ! dummy loop indices
1376	LOGICAL :: ll_tmp1, ll_tmp2
1377	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: sshn
1378	REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: zcpx, zcpy
1379	!!----------------------------------------------------------------------
1380	DO jj = 2, jpjm1
1381	DO ji = 2, jpim1
1382	ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji+1,jj) ) > &
1383	& MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) .AND. &
1384	& MAX( sshn(ji,jj) + ht_0(ji,jj) , sshn(ji+1,jj) + ht_0(ji+1,jj) ) &
1385	& > rn_wdmin1 + rn_wdmin2
1386	ll_tmp2 = ( ABS( sshn(ji+1,jj) - sshn(ji ,jj)) > 1.E-12 ).AND.( &
1387	& MAX( sshn(ji,jj) , sshn(ji+1,jj) ) > &
1388	& MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 )
1389	IF(ll_tmp1) THEN
1390	zcpx(ji,jj) = 1.0_wp
1391	ELSEIF(ll_tmp2) THEN
1392	! no worries about sshn(ji+1,jj) - sshn(ji ,jj) = 0, it won't happen ! here
1393	zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_0(ji+1,jj) - sshn(ji,jj) - ht_0(ji,jj)) &
1394	& / (sshn(ji+1,jj) - sshn(ji ,jj)) )
1395	zcpx(ji,jj) = max(min( zcpx(ji,jj) , 1.0_wp),0.0_wp)
1396	ELSE
1397	zcpx(ji,jj) = 0._wp
1398	ENDIF
1399	!
1400	ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji,jj+1) ) > &
1401	& MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) .AND. &
1402	& MAX( sshn(ji,jj) + ht_0(ji,jj) , sshn(ji,jj+1) + ht_0(ji,jj+1) ) &
1403	& > rn_wdmin1 + rn_wdmin2
1404	ll_tmp2 = ( ABS( sshn(ji,jj) - sshn(ji,jj+1)) > 1.E-12 ).AND.( &
1405	& MAX( sshn(ji,jj) , sshn(ji,jj+1) ) > &
1406	& MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 )
1407
1408	IF(ll_tmp1) THEN
1409	zcpy(ji,jj) = 1.0_wp
1410	ELSE IF(ll_tmp2) THEN
1411	! no worries about sshn(ji,jj+1) - sshn(ji,jj ) = 0, it won't happen ! here
1412	zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_0(ji,jj+1) - sshn(ji,jj) - ht_0(ji,jj)) &
1413	& / (sshn(ji,jj+1) - sshn(ji,jj )) )
1414	zcpy(ji,jj) = MAX( 0._wp , MIN( zcpy(ji,jj) , 1.0_wp ) )
1415	ELSE
1416	zcpy(ji,jj) = 0._wp
1417	ENDIF
1418	END DO
1419	END DO
1420
1421	END SUBROUTINE wad_spg
1422
1423
1424
1425	SUBROUTINE dyn_drg_init( pu_RHSi, pv_RHSi, pCdU_u, pCdU_v )
1426	!!----------------------------------------------------------------------
1427	!! * ROUTINE dyn_drg_init *
1428	!!
1429	!! ** Purpose : - add the baroclinic top/bottom drag contribution to
1430	!! the baroclinic part of the barotropic RHS
1431	!! - compute the barotropic drag coefficients
1432	!!
1433	!! ** Method : computation done over the INNER domain only
1434	!!----------------------------------------------------------------------
1435	REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pu_RHSi, pv_RHSi ! baroclinic part of the barotropic RHS
1436	REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: pCdU_u , pCdU_v ! barotropic drag coefficients
1437	!
1438	INTEGER :: ji, jj ! dummy loop indices
1439	INTEGER :: ikbu, ikbv, iktu, iktv
1440	REAL(wp) :: zztmp
1441	REAL(wp), DIMENSION(jpi,jpj) :: zu_i, zv_i
1442	!!----------------------------------------------------------------------
1443	!
1444	! !== Set the barotropic drag coef. ==!
1445	!
1446	IF( ln_isfcav ) THEN ! top+bottom friction (ocean cavities)
1447
1448	DO jj = 2, jpjm1
1449	DO ji = 2, jpim1 ! INNER domain
1450	pCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) + rCdU_top(ji+1,jj)+rCdU_top(ji,jj) )
1451	pCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) + rCdU_top(ji,jj+1)+rCdU_top(ji,jj) )
1452	END DO
1453	END DO
1454	ELSE ! bottom friction only
1455	DO jj = 2, jpjm1
1456	DO ji = 2, jpim1 ! INNER domain
1457	pCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) )
1458	pCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) )
1459	END DO
1460	END DO
1461	ENDIF
1462	!
1463	! !== BOTTOM stress contribution from baroclinic velocities ==!
1464	!
1465	IF( ln_bt_fw ) THEN ! FORWARD integration: use NOW bottom baroclinic velocities
1466
1467	DO jj = 2, jpjm1
1468	DO ji = 2, jpim1 ! INNER domain
1469	ikbu = mbku(ji,jj)
1470	ikbv = mbkv(ji,jj)
1471	zu_i(ji,jj) = un(ji,jj,ikbu) - un_b(ji,jj)
1472	zv_i(ji,jj) = vn(ji,jj,ikbv) - vn_b(ji,jj)
1473	END DO
1474	END DO
1475	ELSE ! CENTRED integration: use BEFORE bottom baroclinic velocities
1476
1477	DO jj = 2, jpjm1
1478	DO ji = 2, jpim1 ! INNER domain
1479	ikbu = mbku(ji,jj)
1480	ikbv = mbkv(ji,jj)
1481	zu_i(ji,jj) = ub(ji,jj,ikbu) - ub_b(ji,jj)
1482	zv_i(ji,jj) = vb(ji,jj,ikbv) - vb_b(ji,jj)
1483	END DO
1484	END DO
1485	ENDIF
1486	!
1487	IF( ln_wd_il ) THEN ! W/D : use the "clipped" bottom friction !!gm explain WHY, please !
1488	zztmp = -1._wp / rdtbt
1489	DO jj = 2, jpjm1
1490	DO ji = 2, jpim1 ! INNER domain
1491	pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + zu_i(ji,jj) * wdrampu(ji,jj) * MAX( &
1492	& r1_hu_n(ji,jj) * r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) , zztmp )
1493	pv_RHSi(ji,jj) = pv_RHSi(ji,jj) + zv_i(ji,jj) * wdrampv(ji,jj) * MAX( &
1494	& r1_hv_n(ji,jj) * r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) , zztmp )
1495	END DO
1496	END DO
1497	ELSE ! use "unclipped" drag (even if explicit friction is used in 3D calculation)
1498
1499	DO jj = 2, jpjm1
1500	DO ji = 2, jpim1 ! INNER domain
1501	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)
1502	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)
1503	END DO
1504	END DO
1505	END IF
1506	!
1507	! !== TOP stress contribution from baroclinic velocities ==! (no W/D case)
1508	!
1509	IF( ln_isfcav ) THEN
1510	!
1511	IF( ln_bt_fw ) THEN ! FORWARD integration: use NOW top baroclinic velocity
1512
1513	DO jj = 2, jpjm1
1514	DO ji = 2, jpim1 ! INNER domain
1515	iktu = miku(ji,jj)
1516	iktv = mikv(ji,jj)
1517	zu_i(ji,jj) = un(ji,jj,iktu) - un_b(ji,jj)
1518	zv_i(ji,jj) = vn(ji,jj,iktv) - vn_b(ji,jj)
1519	END DO
1520	END DO
1521	ELSE ! CENTRED integration: use BEFORE top baroclinic velocity
1522
1523	DO jj = 2, jpjm1
1524	DO ji = 2, jpim1 ! INNER domain
1525	iktu = miku(ji,jj)
1526	iktv = mikv(ji,jj)
1527	zu_i(ji,jj) = ub(ji,jj,iktu) - ub_b(ji,jj)
1528	zv_i(ji,jj) = vb(ji,jj,iktv) - vb_b(ji,jj)
1529	END DO
1530	END DO
1531	ENDIF
1532	!
1533	! ! use "unclipped" top drag (even if explicit friction is used in 3D calculation)
1534
1535	DO jj = 2, jpjm1
1536	DO ji = 2, jpim1 ! INNER domain
1537	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)
1538	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)
1539	END DO
1540	END DO
1541	!
1542	ENDIF
1543	!
1544	END SUBROUTINE dyn_drg_init
1545
1546	SUBROUTINE ts_bck_interp( jn, ll_init, & ! <== in
1547	& za0, za1, za2, za3 ) ! ==> out
1548	!!----------------------------------------------------------------------
1549	INTEGER ,INTENT(in ) :: jn ! index of sub time step
1550	LOGICAL ,INTENT(in ) :: ll_init !
1551	REAL(wp),INTENT( out) :: za0, za1, za2, za3 ! Half-step back interpolation coefficient
1552	!
1553	REAL(wp) :: zepsilon, zgamma ! - -
1554	!!----------------------------------------------------------------------
1555	! ! set Half-step back interpolation coefficient
1556	IF ( jn==1 .AND. ll_init ) THEN !* Forward-backward
1557	za0 = 1._wp
1558	za1 = 0._wp
1559	za2 = 0._wp
1560	za3 = 0._wp
1561	ELSEIF( jn==2 .AND. ll_init ) THEN !* AB2-AM3 Coefficients; bet=0 ; gam=-1/6 ; eps=1/12
1562	za0 = 1.0833333333333_wp ! za0 = 1-gam-eps
1563	za1 =-0.1666666666666_wp ! za1 = gam
1564	za2 = 0.0833333333333_wp ! za2 = eps
1565	za3 = 0._wp
1566	ELSE !* AB3-AM4 Coefficients; bet=0.281105 ; eps=0.013 ; gam=0.0880
1567	IF( rn_bt_alpha == 0._wp ) THEN ! Time diffusion
1568	za0 = 0.614_wp ! za0 = 1/2 + gam + 2*eps
1569	za1 = 0.285_wp ! za1 = 1/2 - 2gam - 3eps
1570	za2 = 0.088_wp ! za2 = gam
1571	za3 = 0.013_wp ! za3 = eps
1572	ELSE ! no time diffusion
1573	zepsilon = 0.00976186_wp - 0.13451357_wp * rn_bt_alpha
1574	zgamma = 0.08344500_wp - 0.51358400_wp * rn_bt_alpha
1575	za0 = 0.5_wp + zgamma + 2._wp * rn_bt_alpha + 2._wp * zepsilon
1576	za1 = 1._wp - za0 - zgamma - zepsilon
1577	za2 = zgamma
1578	za3 = zepsilon
1579	ENDIF
1580	ENDIF
1581	END SUBROUTINE ts_bck_interp
1582
1583
1584	!!======================================================================
1585	END MODULE dynspg_ts

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

Download in other formats: