New URL for NEMO forge! http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

dynspg_ts.F90 in NEMO/branches/2020/dev_r13333_KERNEL-08_techene_gm_HPG_SPG/src/OCE/DYN – NEMO

Context Navigation

source: NEMO/branches/2020/dev_r13333_KERNEL-08_techene_gm_HPG_SPG/src/OCE/DYN/dynspg_ts.F90 @ 14037

Last change on this file since 14037 was 14037, checked in by ayoung, 3 years ago
Updated to trunk at 14020. Sette tests passed with change of results for configurations with non-linear ssh. Ticket #2506.
Property svn:keywords set to `Id`
File size: 77.2 KB

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

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

Download in other formats: