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/2021/ENHANCE-01_davestorkey_fix_3D_momentum_trends/src/OCE/DYN – NEMO

Context Navigation

source: NEMO/branches/2021/ENHANCE-01_davestorkey_fix_3D_momentum_trends/src/OCE/DYN/dynspg_ts.F90 @ 14701

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

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

Download in other formats: