Context Navigation

dynspg_ts.F90 @ 15540

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

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

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

Context Navigation

source: NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/DYN/dynspg_ts.F90 @ 15540

Download in other formats: