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 branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

Context Navigation

source: branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 4288

Last change on this file since 4288 was 4288, checked in by cbricaud, 10 years ago
correction for timesplitting
Property svn:keywords set to `Id`
File size: 55.9 KB

Line
1	MODULE dynspg_ts
2	!!======================================================================
3	!! History : 1.0 ! 2004-12 (L. Bessieres, G. Madec) Original code
4	!! - ! 2005-11 (V. Garnier, G. Madec) optimization
5	!! - ! 2006-08 (S. Masson) distributed restart using iom
6	!! 2.0 ! 2007-07 (D. Storkey) calls to BDY routines
7	!! - ! 2008-01 (R. Benshila) change averaging method
8	!! 3.2 ! 2009-07 (R. Benshila, G. Madec) Complete revisit associated to vvl reactivation
9	!! 3.3 ! 2010-09 (D. Storkey, E. O'Dea) update for BDY for Shelf configurations
10	!! 3.3 ! 2011-03 (R. Benshila, R. Hordoir, P. Oddo) update calculation of ub_b
11	!! 3.5 ! 2013-07 (J. Chanut) Switch to Forward-backward time stepping
12	!! 3.6 ! 2013-11 (A. Coward) Update for z-tilde compatibility
13	!!---------------------------------------------------------------------
14	#if defined key_dynspg_ts \|\| defined key_esopa
15	!!----------------------------------------------------------------------
16	!! 'key_dynspg_ts' free surface cst volume with time splitting
17	!!----------------------------------------------------------------------
18	!! dyn_spg_ts : compute surface pressure gradient trend using a time-
19	!! splitting scheme and add to the general trend
20	!!----------------------------------------------------------------------
21	USE oce ! ocean dynamics and tracers
22	USE dom_oce ! ocean space and time domain
23	USE sbc_oce ! surface boundary condition: ocean
24	USE dynspg_oce ! surface pressure gradient variables
25	USE phycst ! physical constants
26	USE domvvl ! variable volume
27	USE dynvor ! vorticity term
28	USE bdy_par ! for lk_bdy
29	USE bdy_oce ! Lateral open boundary condition
30	USE bdytides ! open boundary condition data
31	USE bdydyn2d ! open boundary conditions on barotropic variables
32	USE sbctide ! tides
33	USE updtide ! tide potential
34	USE lib_mpp ! distributed memory computing library
35	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
36	USE prtctl ! Print control
37	USE in_out_manager ! I/O manager
38	USE iom ! IOM library
39	USE restart ! only for lrst_oce
40	USE zdf_oce ! Vertical diffusion
41	USE wrk_nemo ! Memory Allocation
42	USE timing ! Timing
43	USE sbcapr ! surface boundary condition: atmospheric pressure
44	USE dynadv, ONLY: ln_dynadv_vec
45	#if defined key_agrif
46	USE agrif_opa_interp ! agrif
47	#endif
48
49
50	IMPLICIT NONE
51	PRIVATE
52
53	PUBLIC dyn_spg_ts ! routine called in dynspg.F90
54	PUBLIC dyn_spg_ts_alloc ! " " " "
55	PUBLIC dyn_spg_ts_init ! " " " "
56
57	! Potential namelist parameters below to be read in dyn_spg_ts_init
58	LOGICAL, PUBLIC, PARAMETER :: ln_bt_fw=.TRUE. !: Forward integration of barotropic sub-stepping
59	LOGICAL, PRIVATE, PARAMETER :: ln_bt_av=.TRUE. !: Time averaging of barotropic variables
60	LOGICAL, PRIVATE, PARAMETER :: ln_bt_nn_auto=.FALSE. !: Set number of iterations automatically
61	INTEGER, PRIVATE, PARAMETER :: nn_bt_flt=1 !: Filter choice
62	REAL(wp), PRIVATE, PARAMETER :: rn_bt_cmax=0.8_wp !: Max. courant number (used if ln_bt_nn_auto=T)
63	! End namelist parameters
64
65	INTEGER, SAVE :: icycle ! Number of barotropic sub-steps for each internal step nn_baro <= 2.5 nn_baro
66	REAL(wp),SAVE :: rdtbt ! Barotropic time step
67
68	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: &
69	wgtbtp1, & ! Primary weights used for time filtering of barotropic variables
70	wgtbtp2 ! Secondary weights used for time filtering of barotropic variables
71
72	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zwz ! ff/h at F points
73	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ftnw, ftne ! triad of coriolis parameter
74	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ftsw, ftse ! (only used with een vorticity scheme)
75
76	! Arrays below are saved to allow testing of the "no time averaging" option
77	! If this option is not retained, these could be replaced by temporary arrays
78	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: sshbb_e, sshb_e, & ! Instantaneous barotropic arrays
79	ubb_e, ub_e, &
80	vbb_e, vb_e
81
82	!! * Substitutions
83	# include "domzgr_substitute.h90"
84	# include "vectopt_loop_substitute.h90"
85	!!----------------------------------------------------------------------
86	!! NEMO/OPA 3.5 , NEMO Consortium (2013)
87	!! $Id: dynspg_ts.F90
88	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
89	!!----------------------------------------------------------------------
90	CONTAINS
91
92	INTEGER FUNCTION dyn_spg_ts_alloc()
93	!!----------------------------------------------------------------------
94	!! * routine dyn_spg_ts_alloc *
95	!!----------------------------------------------------------------------
96	INTEGER :: ierr(3)
97	!!----------------------------------------------------------------------
98	ierr(:) = 0
99
100	ALLOCATE( sshb_e(jpi,jpj), sshbb_e(jpi,jpj), &
101	& ub_e(jpi,jpj) , vb_e(jpi,jpj) , &
102	& ubb_e(jpi,jpj) , vbb_e(jpi,jpj) , STAT= ierr(1) )
103
104	ALLOCATE( wgtbtp1(3nn_baro), wgtbtp2(3nn_baro), zwz(jpi,jpj), STAT= ierr(2) )
105
106	IF( ln_dynvor_een ) ALLOCATE( ftnw(jpi,jpj) , ftne(jpi,jpj) , &
107	& ftsw(jpi,jpj) , ftse(jpi,jpj) , STAT=ierr(3) )
108
109	dyn_spg_ts_alloc = MAXVAL(ierr(:))
110
111	IF( lk_mpp ) CALL mpp_sum( dyn_spg_ts_alloc )
112	IF( dyn_spg_ts_alloc /= 0 ) CALL ctl_warn('dynspg_oce_alloc: failed to allocate arrays')
113	!
114	END FUNCTION dyn_spg_ts_alloc
115
116	SUBROUTINE dyn_spg_ts( kt )
117	!!----------------------------------------------------------------------
118	!!
119	!! ** Purpose :
120	!! -Compute the now trend due to the explicit time stepping
121	!! of the quasi-linear barotropic system. Barotropic variables are
122	!! advanced from internal time steps "n" to "n+1" (if ln_bt_cen=F)
123	!! or from "n-1" to "n+1" time steps (if ln_bt_cen=T) with a
124	!! generalized forward-backward (see ref. below) time stepping.
125	!! -Update the free surface at step "n+1" (ssha, zsshu_a, zsshv_a).
126	!! -Compute barotropic advective velocities at step "n" to be used
127	!! to advect tracers latter on. These are compliant with discrete
128	!! continuity equation taken at the baroclinic time steps, thus
129	!! ensuring tracers conservation.
130	!!
131	!! ** Method :
132	!!
133	!! ** Action : - Update barotropic velocities: ua_b, va_b
134	!! - Update trend (ua,va) with barotropic component
135	!! - Update ssha, zsshu_a, zsshv_a
136	!! - Update barotropic advective velocity at kt=now
137	!!
138	!! References : Shchepetkin, A.F. and J.C. McWilliams, 2005:
139	!! The regional oceanic modeling system (ROMS):
140	!! a split-explicit, free-surface,
141	!! topography-following-coordinate oceanic model.
142	!! Ocean Modelling, 9, 347-404.
143	!!---------------------------------------------------------------------
144	!
145	INTEGER, INTENT(in) :: kt ! ocean time-step index
146	!
147	LOGICAL :: ll_fw_start ! if true, forward integration
148	LOGICAL :: ll_init ! if true, special startup of 2d equations
149	INTEGER :: ji, jj, jk, jn ! dummy loop indices
150	INTEGER :: ikbu, ikbv, noffset ! local integers
151	REAL(wp) :: zraur, z1_2dt_b, z2dt_bf ! local scalars
152	REAL(wp) :: zx1, zy1, zx2, zy2 ! - -
153	REAL(wp) :: z1_12, z1_8, z1_4, z1_2 ! - -
154	REAL(wp) :: zu_spg, zv_spg ! - -
155	REAL(wp) :: zhura, zhvra ! - -
156	REAL(wp) :: za0, za1, za2, za3 ! - -
157	!
158	REAL(wp), POINTER, DIMENSION(:,:) :: zun_e, zvn_e, zsshp2_e
159	REAL(wp), POINTER, DIMENSION(:,:) :: zu_trd, zv_trd, zu_frc, zv_frc, zssh_frc
160	REAL(wp), POINTER, DIMENSION(:,:) :: zu_sum, zv_sum, zwx, zwy, zhdiv
161	REAL(wp), POINTER, DIMENSION(:,:) :: zhup2_e, zhvp2_e, zhust_e, zhvst_e
162	REAL(wp), POINTER, DIMENSION(:,:) :: zhur_b, zhvr_b
163	REAL(wp), POINTER, DIMENSION(:,:) :: zsshu_a, zsshv_a
164	REAL(wp), POINTER, DIMENSION(:,:) :: zht, zhf
165	!!----------------------------------------------------------------------
166	!
167	IF( nn_timing == 1 ) CALL timing_start('dyn_spg_ts')
168	!
169	! !* Allocate temporay arrays
170	CALL wrk_alloc( jpi, jpj, zsshp2_e, zhdiv )
171	CALL wrk_alloc( jpi, jpj, zu_trd, zv_trd, zun_e, zvn_e )
172	CALL wrk_alloc( jpi, jpj, zwx, zwy, zu_sum, zv_sum, zssh_frc, zu_frc, zv_frc)
173	CALL wrk_alloc( jpi, jpj, zhup2_e, zhvp2_e, zhust_e, zhvst_e)
174	CALL wrk_alloc( jpi, jpj, zhur_b, zhvr_b )
175	CALL wrk_alloc( jpi, jpj, zsshu_a, zsshv_a )
176	CALL wrk_alloc( jpi, jpj, zht, zhf )
177	!
178	! !* Local constant initialization
179	z1_12 = 1._wp / 12._wp
180	z1_8 = 0.125_wp
181	z1_4 = 0.25_wp
182	z1_2 = 0.5_wp
183	zraur = 1._wp / rau0
184	!
185	IF( kt == nit000 .AND. neuler == 0 ) THEN ! reciprocal of baroclinic time step
186	z2dt_bf = rdt
187	ELSE
188	z2dt_bf = 2.0_wp * rdt
189	ENDIF
190	z1_2dt_b = 1.0_wp / z2dt_bf
191	!
192	ll_init = ln_bt_av ! if no time averaging, then no specific restart
193	ll_fw_start = .FALSE.
194	!
195	! time offset in steps for bdy data update
196	IF (.NOT.ln_bt_fw) THEN ; noffset=-2*nn_baro ; ELSE ; noffset = 0 ; ENDIF
197	!
198	IF( kt == nit000 ) THEN !* initialisation
199	!
200	IF(lwp) WRITE(numout,*)
201	IF(lwp) WRITE(numout,*) 'dyn_spg_ts : surface pressure gradient trend'
202	IF(lwp) WRITE(numout,*) '~~~~~~~~~~ free surface with time splitting'
203	IF(lwp) WRITE(numout,) ' Number of sub cycle in 1 time-step (2 rdt) : icycle = ', 2nn_baro
204	!
205	IF (neuler==0) ll_init=.TRUE.
206	!
207	IF (ln_bt_fw.OR.(neuler==0)) THEN
208	ll_fw_start=.TRUE.
209	noffset = 0
210	ELSE
211	ll_fw_start=.FALSE.
212	ENDIF
213	!
214	! Set averaging weights and cycle length:
215	CALL ts_wgt(ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2)
216	!
217	IF ((neuler/=0).AND.(ln_bt_fw)) CALL ts_rst( nit000, 'READ' )
218	!
219	ENDIF
220	!
221	! Set arrays to remove/compute coriolis trend.
222	! Do it once at kt=nit000 if volume is fixed, else at each long time step.
223	! Note that these arrays are also used during barotropic loop. These are however frozen
224	! although they should be updated in variable volume case. Not a big approximation.
225	! To remove this approximation, copy lines below inside barotropic loop
226	! and update depths at T-F points (ht and hf resp.) at each barotropic time step
227	!
228	IF ( kt == nit000 .OR. lk_vvl ) THEN
229	IF ( ln_dynvor_een ) THEN
230	! JC: Simplification needed below: define ht_0 even when volume is fixed
231	IF (lk_vvl) THEN
232	zht(:,:) = (ht_0(:,:) + sshn(:,:)) * tmask(:,:,1)
233	ELSE
234	zht(:,:) = 0.
235	DO jk = 1, jpkm1
236	zht(:,:) = zht(:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk)
237	END DO
238	ENDIF
239
240	DO jj = 1, jpjm1
241	DO ji = 1, jpim1
242	zwz(ji,jj) = ( zht(ji ,jj+1) + zht(ji+1,jj+1) + &
243	& zht(ji ,jj ) + zht(ji+1,jj ) ) &
244	& / ( MAX( 1.0_wp, tmask(ji ,jj+1, 1) + tmask(ji+1,jj+1, 1) + &
245	& tmask(ji ,jj , 1) + tmask(ji+1,jj , 1) ) )
246	IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = 1._wp / zwz(ji,jj)
247	END DO
248	END DO
249	CALL lbc_lnk( zwz, 'F', 1._wp )
250	zwz(:,:) = ff(:,:) * zwz(:,:)
251
252	ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp
253	DO jj = 2, jpj
254	DO ji = fs_2, jpi ! vector opt.
255	ftne(ji,jj) = zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1)
256	ftnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj ) + zwz(ji ,jj )
257	ftse(ji,jj) = zwz(ji ,jj ) + zwz(ji ,jj-1) + zwz(ji-1,jj-1)
258	ftsw(ji,jj) = zwz(ji ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj )
259	END DO
260	END DO
261	ELSE
262	zwz(:,:) = 0._wp
263	zht(:,:) = 0.
264	IF ( .not. ln_sco ) THEN
265	! IF( rn_hmin < 0._wp ) THEN ; jk = - INT( rn_hmin ) ! from a nb of level
266	! ELSE ; jk = MINLOC( gdepw_0, mask = gdepw_0 > rn_hmin, dim = 1 ) ! from a depth
267	! ENDIF
268	! zht(:,:) = gdepw_0(:,:,jk+1)
269	ELSE
270	zht(:,:) = hbatf(:,:)
271	END IF
272
273	DO jj = 1, jpjm1
274	zht(:,jj) = zht(:,jj)(1._wp- umask(:,jj,1) umask(:,jj+1,1))
275	END DO
276
277	DO jk = 1, jpkm1
278	DO jj = 1, jpjm1
279	zht(:,jj) = zht(:,jj) + fse3f_n(:,jj,jk) * umask(:,jj,jk) * umask(:,jj+1,jk)
280	END DO
281	END DO
282	CALL lbc_lnk( zht, 'F', 1._wp )
283	! JC: TBC. hf should be greater than 0
284	DO jj = 1, jpj
285	DO ji = 1, jpi
286	IF( zht(ji,jj) /= 0._wp ) zwz(ji,jj) = 1._wp / zht(ji,jj) ! zht is actually hf here but it saves an array
287	END DO
288	END DO
289	zwz(:,:) = ff(:,:) * zwz(:,:)
290	ENDIF
291	ENDIF
292	!
293	! If forward start at previous time step, and centered integration,
294	! then update averaging weights:
295	IF ((.NOT.ln_bt_fw).AND.((neuler==0).AND.(kt==nit000+1))) THEN
296	ll_fw_start=.FALSE.
297	CALL ts_wgt(ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2)
298	ENDIF
299
300	! before inverse water column height at u- and v- points
301	IF( lk_vvl ) THEN
302	zhur_b(:,:) = 0.
303	zhvr_b(:,:) = 0.
304	DO jk = 1, jpk
305	zhur_b(:,:) = zhur_b(:,:) + fse3u_b(:,:,jk) * umask(:,:,jk)
306	zhvr_b(:,:) = zhvr_b(:,:) + fse3v_b(:,:,jk) * vmask(:,:,jk)
307	END DO
308	zhur_b(:,:) = umask(:,:,1) / ( zhur_b(:,:) + 1. - umask(:,:,1) )
309	zhvr_b(:,:) = vmask(:,:,1) / ( zhvr_b(:,:) + 1. - vmask(:,:,1) )
310	ELSE
311	zhur_b(:,:) = hur(:,:)
312	zhvr_b(:,:) = hvr(:,:)
313	ENDIF
314
315	! -----------------------------------------------------------------------------
316	! Phase 1 : Coupling between general trend and barotropic estimates (1st step)
317	! -----------------------------------------------------------------------------
318	!
319	! Some vertical sums (at now and before time steps) below could be suppressed
320	! if one swap barotropic arrays somewhere
321	!
322	! !* e3d/dt(Ua), e3Ub, e3*Vn (Vertically integrated)
323	! ! --------------------------------------------------
324	zu_frc(:,:) = 0._wp ; ub_b(:,:) = 0._wp ; un_b(:,:) = 0._wp
325	zv_frc(:,:) = 0._wp ; vb_b(:,:) = 0._wp ; vn_b(:,:) = 0._wp
326	!
327	DO jk = 1, jpkm1
328	#if defined key_vectopt_loop
329	DO jj = 1, 1 !Vector opt. => forced unrolling
330	DO ji = 1, jpij
331	#else
332	DO jj = 1, jpj
333	DO ji = 1, jpi
334	#endif
335	! ! now trend:
336	zu_frc(ji,jj) = zu_frc(ji,jj) + fse3u(ji,jj,jk) * ua(ji,jj,jk) * umask(ji,jj,jk)
337	zv_frc(ji,jj) = zv_frc(ji,jj) + fse3v(ji,jj,jk) * va(ji,jj,jk) * vmask(ji,jj,jk)
338	! ! now bt transp:
339	un_b(ji,jj) = un_b(ji,jj) + fse3u(ji,jj,jk) * un(ji,jj,jk) * umask(ji,jj,jk)
340	vn_b(ji,jj) = vn_b(ji,jj) + fse3v(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk)
341	! ! before bt transp:
342	ub_b(ji,jj) = ub_b(ji,jj) + fse3u_b(ji,jj,jk) * ub(ji,jj,jk) * umask(ji,jj,jk)
343	vb_b(ji,jj) = vb_b(ji,jj) + fse3v_b(ji,jj,jk) * vb(ji,jj,jk) * vmask(ji,jj,jk)
344	END DO
345	END DO
346	END DO
347	!
348	zu_frc(:,:) = zu_frc(:,:) * hur(:,:)
349	zv_frc(:,:) = zv_frc(:,:) * hvr(:,:)
350	!
351	IF( lk_vvl ) THEN
352	ub_b(:,:) = ub_b(:,:) * zhur_b(:,:)
353	vb_b(:,:) = vb_b(:,:) * zhvr_b(:,:)
354	ELSE
355	ub_b(:,:) = ub_b(:,:) * hur(:,:)
356	vb_b(:,:) = vb_b(:,:) * hvr(:,:)
357	ENDIF
358	!
359	! !* baroclinic momentum trend (remove the vertical mean trend)
360	DO jk = 1, jpkm1 ! -----------------------------------------------------------
361	DO jj = 2, jpjm1
362	DO ji = fs_2, fs_jpim1 ! vector opt.
363	ua(ji,jj,jk) = ua(ji,jj,jk) - zu_frc(ji,jj) * umask(ji,jj,jk)
364	va(ji,jj,jk) = va(ji,jj,jk) - zv_frc(ji,jj) * vmask(ji,jj,jk)
365	END DO
366	END DO
367	END DO
368	! !* barotropic Coriolis trends (vorticity scheme dependent)
369	! ! --------------------------------------------------------
370	zwx(:,:) = un_b(:,:) * e2u(:,:) ! now transport
371	zwy(:,:) = vn_b(:,:) * e1v(:,:)
372	!
373	IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN ! energy conserving or mixed scheme
374	DO jj = 2, jpjm1
375	DO ji = fs_2, fs_jpim1 ! vector opt.
376	zy1 = ( zwy(ji,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj)
377	zy2 = ( zwy(ji,jj ) + zwy(ji+1,jj ) ) / e1u(ji,jj)
378	zx1 = ( zwx(ji-1,jj) + zwx(ji-1,jj+1) ) / e2v(ji,jj)
379	zx2 = ( zwx(ji ,jj) + zwx(ji ,jj+1) ) / e2v(ji,jj)
380	! energy conserving formulation for planetary vorticity term
381	zu_trd(ji,jj) = z1_4 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
382	zv_trd(ji,jj) =-z1_4 * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 )
383	END DO
384	END DO
385	!
386	ELSEIF ( ln_dynvor_ens ) THEN ! enstrophy conserving scheme
387	DO jj = 2, jpjm1
388	DO ji = fs_2, fs_jpim1 ! vector opt.
389	zy1 = z1_8 * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) &
390	& + zwy(ji ,jj ) + zwy(ji+1,jj ) ) / e1u(ji,jj)
391	zx1 = - z1_8 * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) &
392	& + zwx(ji ,jj ) + zwx(ji ,jj+1) ) / e2v(ji,jj)
393	zu_trd(ji,jj) = zy1 * ( zwz(ji ,jj-1) + zwz(ji,jj) )
394	zv_trd(ji,jj) = zx1 * ( zwz(ji-1,jj ) + zwz(ji,jj) )
395	END DO
396	END DO
397	!
398	ELSEIF ( ln_dynvor_een ) THEN ! enstrophy and energy conserving scheme
399	DO jj = 2, jpjm1
400	DO ji = fs_2, fs_jpim1 ! vector opt.
401	zu_trd(ji,jj) = + z1_12 / e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) &
402	& + ftnw(ji+1,jj) * zwy(ji+1,jj ) &
403	& + ftse(ji,jj ) * zwy(ji ,jj-1) &
404	& + ftsw(ji+1,jj) * zwy(ji+1,jj-1) )
405	zv_trd(ji,jj) = - z1_12 / e2v(ji,jj) * ( ftsw(ji,jj+1) * zwx(ji-1,jj+1) &
406	& + ftse(ji,jj+1) * zwx(ji ,jj+1) &
407	& + ftnw(ji,jj ) * zwx(ji-1,jj ) &
408	& + ftne(ji,jj ) * zwx(ji ,jj ) )
409	END DO
410	END DO
411	!
412	ENDIF
413	!
414	un_b (:,:) = un_b(:,:) * hur(:,:) ! Revert now transport to barotropic velocities
415	vn_b (:,:) = vn_b(:,:) * hvr(:,:)
416	! !* Right-Hand-Side of the barotropic momentum equation
417	! ! ----------------------------------------------------
418	IF( lk_vvl ) THEN ! Variable volume : remove surface pressure gradient
419	DO jj = 2, jpjm1
420	DO ji = fs_2, fs_jpim1 ! vector opt.
421	zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) / e1u(ji,jj)
422	zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji ,jj+1) - sshn(ji ,jj ) ) / e2v(ji,jj)
423	END DO
424	END DO
425	ENDIF
426
427	DO jj = 2, jpjm1 ! Remove coriolis term (and possibly spg) from barotropic trend
428	DO ji = fs_2, fs_jpim1
429	zu_frc(ji,jj) = zu_frc(ji,jj) - zu_trd(ji,jj) * umask(ji,jj,1)
430	zv_frc(ji,jj) = zv_frc(ji,jj) - zv_trd(ji,jj) * vmask(ji,jj,1)
431	END DO
432	END DO
433	!
434	! ! Add bottom stress contribution from baroclinic velocities:
435	IF (ln_bt_fw) THEN
436	DO jj = 2, jpjm1
437	DO ji = fs_2, fs_jpim1 ! vector opt.
438	ikbu = mbku(ji,jj)
439	ikbv = mbkv(ji,jj)
440	zwx(ji,jj) = un(ji,jj,ikbu) - un_b(ji,jj) ! NOW bottom baroclinic velocities
441	zwy(ji,jj) = vn(ji,jj,ikbv) - vn_b(ji,jj)
442	END DO
443	END DO
444	ELSE
445	DO jj = 2, jpjm1
446	DO ji = fs_2, fs_jpim1 ! vector opt.
447	ikbu = mbku(ji,jj)
448	ikbv = mbkv(ji,jj)
449	zwx(ji,jj) = ub(ji,jj,ikbu) - ub_b(ji,jj) ! BEFORE bottom baroclinic velocities
450	zwy(ji,jj) = vb(ji,jj,ikbv) - vb_b(ji,jj)
451	END DO
452	END DO
453	ENDIF
454	!
455	! Note that the "unclipped" bottom friction parameter is used even with explicit drag
456	zu_frc(:,:) = zu_frc(:,:) + hur(:,:) * bfrua(:,:) * zwx(:,:)
457	zv_frc(:,:) = zv_frc(:,:) + hvr(:,:) * bfrva(:,:) * zwy(:,:)
458	!
459	IF (ln_bt_fw) THEN ! Add wind forcing
460	zu_frc(:,:) = zu_frc(:,:) + zraur * utau(:,:) * hur(:,:)
461	zv_frc(:,:) = zv_frc(:,:) + zraur * vtau(:,:) * hvr(:,:)
462	ELSE
463	zu_frc(:,:) = zu_frc(:,:) + zraur * z1_2 * ( utau_b(:,:) + utau(:,:) ) * hur(:,:)
464	zv_frc(:,:) = zv_frc(:,:) + zraur * z1_2 * ( vtau_b(:,:) + vtau(:,:) ) * hvr(:,:)
465	ENDIF
466	!
467	IF ( ln_apr_dyn ) THEN ! Add atm pressure forcing
468	IF (ln_bt_fw) THEN
469	DO jj = 2, jpjm1
470	DO ji = fs_2, fs_jpim1 ! vector opt.
471	zu_spg = grav * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) ) /e1u(ji,jj)
472	zv_spg = grav * ( ssh_ib (ji ,jj+1) - ssh_ib (ji,jj) ) /e2v(ji,jj)
473	zu_frc(ji,jj) = zu_frc(ji,jj) + zu_spg
474	zv_frc(ji,jj) = zv_frc(ji,jj) + zv_spg
475	END DO
476	END DO
477	ELSE
478	DO jj = 2, jpjm1
479	DO ji = fs_2, fs_jpim1 ! vector opt.
480	zu_spg = grav * z1_2 * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) &
481	& + ssh_ibb(ji+1,jj ) - ssh_ibb(ji,jj) ) /e1u(ji,jj)
482	zv_spg = grav * z1_2 * ( ssh_ib (ji ,jj+1) - ssh_ib (ji,jj) &
483	& + ssh_ibb(ji ,jj+1) - ssh_ibb(ji,jj) ) /e2v(ji,jj)
484	zu_frc(ji,jj) = zu_frc(ji,jj) + zu_spg
485	zv_frc(ji,jj) = zv_frc(ji,jj) + zv_spg
486	END DO
487	END DO
488	ENDIF
489	ENDIF
490	! !* Right-Hand-Side of the barotropic ssh equation
491	! ! -----------------------------------------------
492	! ! Surface net water flux and rivers
493	IF (ln_bt_fw) THEN
494	zssh_frc(:,:) = zraur * ( emp(:,:) - rnf(:,:) )
495	ELSE
496	zssh_frc(:,:) = zraur * z1_2 * (emp(:,:) + emp_b(:,:) - rnf(:,:) - rnf_b(:,:))
497	ENDIF
498	#if defined key_asminc
499	! ! Include the IAU weighted SSH increment
500	IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN
501	zssh_frc(:,:) = zssh_frc(:,:) + ssh_iau(:,:)
502	ENDIF
503	#endif
504	!
505	! -----------------------------------------------------------------------
506	! Phase 2 : Integration of the barotropic equations
507	! -----------------------------------------------------------------------
508	!
509	! ! ==================== !
510	! ! Initialisations !
511	! ! ==================== !
512	! Initialize barotropic variables:
513	IF (ln_bt_fw) THEN ! FORWARD integration: start from NOW fields
514	sshn_e (:,:) = sshn (:,:)
515	zun_e (:,:) = un_b (:,:)
516	zvn_e (:,:) = vn_b (:,:)
517	ELSE ! CENTERED integration: start from BEFORE fields
518	sshn_e (:,:) = sshb (:,:)
519	zun_e (:,:) = ub_b (:,:)
520	zvn_e (:,:) = vb_b (:,:)
521	ENDIF
522	!
523	! Initialize depths:
524	IF ( lk_vvl.AND.(.NOT.ln_bt_fw) ) THEN
525	hu_e (:,:) = umask(:,:,1) / ( zhur_b(:,:) + 1._wp - umask(:,:,1) )
526	hv_e (:,:) = vmask(:,:,1) / ( zhvr_b(:,:) + 1._wp - vmask(:,:,1) )
527	hur_e (:,:) = zhur_b(:,:)
528	hvr_e (:,:) = zhvr_b(:,:)
529	ELSE
530	hu_e (:,:) = hu (:,:)
531	hv_e (:,:) = hv (:,:)
532	hur_e (:,:) = hur (:,:)
533	hvr_e (:,:) = hvr (:,:)
534	ENDIF
535	!
536	IF (.NOT.lk_vvl) THEN ! Depths at jn+0.5:
537	zhup2_e (:,:) = hu(:,:)
538	zhvp2_e (:,:) = hv(:,:)
539	ENDIF
540	!
541	! Initialize sums:
542	ua_b (:,:) = 0._wp ! After barotropic velocities (or transport if flux form)
543	va_b (:,:) = 0._wp
544	ssha (:,:) = 0._wp ! Sum for after averaged sea level
545	zu_sum(:,:) = 0._wp ! Sum for now transport issued from ts loop
546	zv_sum(:,:) = 0._wp
547	! ! ==================== !
548	DO jn = 1, icycle ! sub-time-step loop !
549	! ! ==================== !
550	! !* Update the forcing (BDY and tides)
551	! ! ------------------
552	! Update only tidal forcing at open boundaries
553	#if defined key_tide
554	IF ( lk_bdy .AND. lk_tide ) CALL bdy_dta_tides( kt, kit=jn, time_offset=(noffset+1) )
555	IF ( ln_tide_pot .AND. lk_tide ) CALL upd_tide( kt, kit=jn, koffset=noffset )
556	#endif
557	!
558	! Set extrapolation coefficients for predictor step:
559	IF ((jn<3).AND.ll_init) THEN ! Forward
560	za1 = 1._wp
561	za2 = 0._wp
562	za3 = 0._wp
563	ELSE ! AB3-AM4 Coefficients: bet=0.281105
564	za1 = 1.781105_wp ! za1 = 3/2 + bet
565	za2 = -1.06221_wp ! za2 = -(1/2 + 2*bet)
566	za3 = 0.281105_wp ! za3 = bet
567	ENDIF
568
569	! Extrapolate barotropic velocities at step jit+0.5:
570	ua_e(:,:) = za1 * zun_e(:,:) + za2 * ub_e(:,:) + za3 * ubb_e(:,:)
571	va_e(:,:) = za1 * zvn_e(:,:) + za2 * vb_e(:,:) + za3 * vbb_e(:,:)
572
573	IF( lk_vvl ) THEN !* Update ocean depth (variable volume case only)
574	! ! ------------------
575	! Extrapolate Sea Level at step jit+0.5:
576	zsshp2_e(:,:) = za1 * sshn_e(:,:) + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:)
577	!
578	DO jj = 2, jpjm1 ! Sea Surface Height at u- & v-points
579	DO ji = 2, fs_jpim1 ! Vector opt.
580	zwx(ji,jj) = z1_2 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) ) &
581	& * ( e1t(ji ,jj) * e2t(ji ,jj) * zsshp2_e(ji ,jj) &
582	& + e1t(ji+1,jj) * e2t(ji+1,jj) * zsshp2_e(ji+1,jj) )
583	zwy(ji,jj) = z1_2 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) ) &
584	& * ( e1t(ji,jj ) * e2t(ji,jj ) * zsshp2_e(ji,jj ) &
585	& + e1t(ji,jj+1) * e2t(ji,jj+1) * zsshp2_e(ji,jj+1) )
586	END DO
587	END DO
588	CALL lbc_lnk( zwx, 'U', 1._wp ) ; CALL lbc_lnk( zwy, 'V', 1._wp )
589	!
590	zhup2_e (:,:) = hu_0(:,:) + zwx(:,:) ! Ocean depth at U- and V-points
591	zhvp2_e (:,:) = hv_0(:,:) + zwy(:,:)
592	ENDIF
593	! !* after ssh
594	! ! -----------
595	! One should enforce volume conservation at open boundaries here
596	! considering fluxes below:
597	!
598	zwx(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:) ! fluxes at jn+0.5
599	zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:)
600	DO jj = 2, jpjm1
601	DO ji = fs_2, fs_jpim1 ! vector opt.
602	zhdiv(ji,jj) = ( zwx(ji,jj) - zwx(ji-1,jj) &
603	& + zwy(ji,jj) - zwy(ji,jj-1) &
604	& ) / ( e1t(ji,jj) * e2t(ji,jj) )
605	END DO
606	END DO
607	!
608	! Sum over sub-time-steps to compute advective velocities
609	za2 = wgtbtp2(jn)
610	zu_sum (:,:) = zu_sum (:,:) + za2 * ua_e (:,:) * zhup2_e (:,:)
611	zv_sum (:,:) = zv_sum (:,:) + za2 * va_e (:,:) * zhvp2_e (:,:)
612	!
613	! Set next sea level:
614	ssha_e(:,:) = ( sshn_e(:,:) - rdtbt * ( zssh_frc(:,:) + zhdiv(:,:) ) ) * tmask(:,:,1)
615	CALL lbc_lnk( ssha_e, 'T', 1._wp )
616
617	#if defined key_bdy
618	! Duplicate sea level across open boundaries (this is only cosmetic if lk_vvl=.false.)
619	IF (lk_bdy) CALL bdy_ssh( ssha_e )
620	#endif
621	#if defined key_agrif
622	IF( .NOT.Agrif_Root() ) CALL agrif_ssh_ts( jn )
623	#endif
624	!
625	! Sea Surface Height at u-,v-points (vvl case only)
626	IF ( lk_vvl ) THEN
627	DO jj = 2, jpjm1
628	DO ji = 2, jpim1 ! NO Vector Opt.
629	zsshu_a(ji,jj) = z1_2 * umask(ji,jj,1) / ( e1u(ji ,jj) * e2u(ji ,jj) ) &
630	& * ( e1t(ji ,jj) * e2t(ji ,jj) * ssha_e(ji ,jj) &
631	& + e1t(ji+1,jj) * e2t(ji+1,jj) * ssha_e(ji+1,jj) )
632	zsshv_a(ji,jj) = z1_2 * vmask(ji,jj,1) / ( e1v(ji,jj ) * e2v(ji,jj ) ) &
633	& * ( e1t(ji,jj ) * e2t(ji,jj ) * ssha_e(ji,jj ) &
634	& + e1t(ji,jj+1) * e2t(ji,jj+1) * ssha_e(ji,jj+1) )
635	END DO
636	END DO
637	CALL lbc_lnk( zsshu_a, 'U', 1._wp ) ; CALL lbc_lnk( zsshv_a, 'V', 1._wp )
638	ENDIF
639	!
640	! Half-step back interpolation of SSH for surface pressure computation:
641	!----------------------------------------------------------------------
642	IF ((jn==1).AND.ll_init) THEN
643	za0=1._wp ! Forward-backward
644	za1=0._wp
645	za2=0._wp
646	za3=0._wp
647	ELSEIF ((jn==2).AND.ll_init) THEN ! AB2-AM3 Coefficients; bet=0 ; gam=-1/6 ; eps=1/12
648	za0= 1.0833333333333_wp ! za0 = 1-gam-eps
649	za1=-0.1666666666666_wp ! za1 = gam
650	za2= 0.0833333333333_wp ! za2 = eps
651	za3= 0._wp
652	ELSE ! AB3-AM4 Coefficients; bet=0.281105 ; eps=0.013 ; gam=0.0880
653	za0=0.614_wp ! za0 = 1/2 + gam + 2*eps
654	za1=0.285_wp ! za1 = 1/2 - 2gam - 3eps
655	za2=0.088_wp ! za2 = gam
656	za3=0.013_wp ! za3 = eps
657	ENDIF
658
659	zsshp2_e(:,:) = za0 * ssha_e(:,:) + za1 * sshn_e (:,:) &
660	& + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:)
661
662	!
663	! Compute associated depths at U and V points:
664	IF ( lk_vvl.AND.(.NOT.ln_dynadv_vec) ) THEN
665	!
666	DO jj = 2, jpjm1
667	DO ji = 2, jpim1
668	zx1 = z1_2 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) ) &
669	& * ( e1t(ji ,jj) * e2t(ji ,jj) * zsshp2_e(ji ,jj) &
670	& + e1t(ji+1,jj) * e2t(ji+1,jj) * zsshp2_e(ji+1,jj) )
671	zy1 = z1_2 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) ) &
672	& * ( e1t(ji,jj ) * e2t(ji,jj ) * zsshp2_e(ji,jj ) &
673	& + e1t(ji,jj+1) * e2t(ji,jj+1) * zsshp2_e(ji,jj+1) )
674	zhust_e(ji,jj) = hu_0(ji,jj) + zx1
675	zhvst_e(ji,jj) = hv_0(ji,jj) + zy1
676	END DO
677	END DO
678	ENDIF
679	!
680	! Add Coriolis trend:
681	! zwz array below or triads normally depend on sea level with key_vvl and should be updated
682	! at each time step. We however keep them constant here for optimization.
683	! Recall that zwx and zwy arrays hold fluxes at this stage:
684	! zwx(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:) ! fluxes at jn+0.5
685	! zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:)
686	!
687	IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN !== energy conserving or mixed scheme ==!
688	DO jj = 2, jpjm1
689	DO ji = fs_2, fs_jpim1 ! vector opt.
690	zy1 = ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj)
691	zy2 = ( zwy(ji ,jj ) + zwy(ji+1,jj ) ) / e1u(ji,jj)
692	zx1 = ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) ) / e2v(ji,jj)
693	zx2 = ( zwx(ji ,jj ) + zwx(ji ,jj+1) ) / e2v(ji,jj)
694	zu_trd(ji,jj) = z1_4 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
695	zv_trd(ji,jj) =-z1_4 * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 )
696	END DO
697	END DO
698	!
699	ELSEIF ( ln_dynvor_ens ) THEN !== enstrophy conserving scheme ==!
700	DO jj = 2, jpjm1
701	DO ji = fs_2, fs_jpim1 ! vector opt.
702	zy1 = z1_8 * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) &
703	& + zwy(ji ,jj ) + zwy(ji+1,jj ) ) / e1u(ji,jj)
704	zx1 = - z1_8 * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) &
705	& + zwx(ji ,jj ) + zwx(ji ,jj+1) ) / e2v(ji,jj)
706	zu_trd(ji,jj) = zy1 * ( zwz(ji ,jj-1) + zwz(ji,jj) )
707	zv_trd(ji,jj) = zx1 * ( zwz(ji-1,jj ) + zwz(ji,jj) )
708	END DO
709	END DO
710	!
711	ELSEIF ( ln_dynvor_een ) THEN !== energy and enstrophy conserving scheme ==!
712	DO jj = 2, jpjm1
713	DO ji = fs_2, fs_jpim1 ! vector opt.
714	zu_trd(ji,jj) = + z1_12 / e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) &
715	& + ftnw(ji+1,jj) * zwy(ji+1,jj ) &
716	& + ftse(ji,jj ) * zwy(ji ,jj-1) &
717	& + ftsw(ji+1,jj) * zwy(ji+1,jj-1) )
718	zv_trd(ji,jj) = - z1_12 / e2v(ji,jj) * ( ftsw(ji,jj+1) * zwx(ji-1,jj+1) &
719	& + ftse(ji,jj+1) * zwx(ji ,jj+1) &
720	& + ftnw(ji,jj ) * zwx(ji-1,jj ) &
721	& + ftne(ji,jj ) * zwx(ji ,jj ) )
722	END DO
723	END DO
724	!
725	ENDIF
726	!
727	! Add tidal astronomical forcing if defined
728	IF ( lk_tide.AND.ln_tide_pot ) THEN
729	DO jj = 2, jpjm1
730	DO ji = fs_2, fs_jpim1 ! vector opt.
731	zu_spg = grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) / e1u(ji,jj)
732	zv_spg = grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) / e2v(ji,jj)
733	zu_trd(ji,jj) = zu_trd(ji,jj) + zu_spg
734	zv_trd(ji,jj) = zv_trd(ji,jj) + zv_spg
735	END DO
736	END DO
737	ENDIF
738	!
739	! Add bottom stresses:
740	zu_trd(:,:) = zu_trd(:,:) + bfrua(:,:) * zun_e(:,:) * hur_e(:,:)
741	zv_trd(:,:) = zv_trd(:,:) + bfrva(:,:) * zvn_e(:,:) * hvr_e(:,:)
742	!
743	! Surface pressure trend:
744	DO jj = 2, jpjm1
745	DO ji = fs_2, fs_jpim1 ! vector opt.
746	! Add surface pressure gradient
747	zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) / e1u(ji,jj)
748	zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) / e2v(ji,jj)
749	zwx(ji,jj) = zu_spg
750	zwy(ji,jj) = zv_spg
751	END DO
752	END DO
753	!
754	! Set next velocities:
755	IF( ln_dynadv_vec .OR. (.NOT. lk_vvl) ) THEN ! Vector form
756	DO jj = 2, jpjm1
757	DO ji = fs_2, fs_jpim1 ! vector opt.
758	ua_e(ji,jj) = ( zun_e(ji,jj) &
759	& + rdtbt * ( zwx(ji,jj) &
760	& + zu_trd(ji,jj) &
761	& + zu_frc(ji,jj) ) &
762	& ) * umask(ji,jj,1)
763
764	va_e(ji,jj) = ( zvn_e(ji,jj) &
765	& + rdtbt * ( zwy(ji,jj) &
766	& + zv_trd(ji,jj) &
767	& + zv_frc(ji,jj) ) &
768	& ) * vmask(ji,jj,1)
769	END DO
770	END DO
771
772	ELSE ! Flux form
773	DO jj = 2, jpjm1
774	DO ji = fs_2, fs_jpim1 ! vector opt.
775
776	zhura = umask(ji,jj,1)/(hu_0(ji,jj) + zsshu_a(ji,jj) + 1._wp - umask(ji,jj,1))
777	zhvra = vmask(ji,jj,1)/(hv_0(ji,jj) + zsshv_a(ji,jj) + 1._wp - vmask(ji,jj,1))
778
779	ua_e(ji,jj) = ( hu_e(ji,jj) * zun_e(ji,jj) &
780	& + rdtbt * ( zhust_e(ji,jj) * zwx(ji,jj) &
781	& + zhup2_e(ji,jj) * zu_trd(ji,jj) &
782	& + hu(ji,jj) * zu_frc(ji,jj) ) &
783	& ) * zhura
784
785	va_e(ji,jj) = ( hv_e(ji,jj) * zvn_e(ji,jj) &
786	& + rdtbt * ( zhvst_e(ji,jj) * zwy(ji,jj) &
787	& + zhvp2_e(ji,jj) * zv_trd(ji,jj) &
788	& + hv(ji,jj) * zv_frc(ji,jj) ) &
789	& ) * zhvra
790	END DO
791	END DO
792	ENDIF
793	!
794	IF( lk_vvl ) THEN !* Update ocean depth (variable volume case only)
795	! ! ----------------------------------------------
796	hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:)
797	hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:)
798	hur_e(:,:) = umask(:,:,1) / ( hu_e(:,:) + 1._wp - umask(:,:,1) )
799	hvr_e(:,:) = vmask(:,:,1) / ( hv_e(:,:) + 1._wp - vmask(:,:,1) )
800	!
801	ENDIF
802	! !* domain lateral boundary
803	! ! -----------------------
804	!
805	CALL lbc_lnk( ua_e , 'U', -1._wp ) ! local domain boundaries
806	CALL lbc_lnk( va_e , 'V', -1._wp )
807
808	#if defined key_bdy
809
810	pssh => ssha_e
811	phur => hur_e
812	phvr => hvr_e
813	pua2d => ua_e
814	pva2d => va_e
815	pub2d => zun_e
816	pvb2d => zvn_e
817
818	IF( lk_bdy ) CALL bdy_dyn2d( kt ) ! open boundaries
819	#endif
820	#if defined key_agrif
821	IF( .NOT.Agrif_Root() ) CALL agrif_dyn_ts( kt, jn ) ! Agrif
822	#endif
823	! !* Swap
824	! ! ----
825	ubb_e (:,:) = ub_e (:,:)
826	ub_e (:,:) = zun_e (:,:)
827	zun_e (:,:) = ua_e (:,:)
828	!
829	vbb_e (:,:) = vb_e (:,:)
830	vb_e (:,:) = zvn_e (:,:)
831	zvn_e (:,:) = va_e (:,:)
832	!
833	sshbb_e(:,:) = sshb_e(:,:)
834	sshb_e (:,:) = sshn_e(:,:)
835	sshn_e (:,:) = ssha_e(:,:)
836
837	! !* Sum over whole bt loop
838	! ! ----------------------
839	za1 = wgtbtp1(jn)
840	IF (( ln_dynadv_vec ).OR. (.NOT. lk_vvl)) THEN ! Sum velocities
841	ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:)
842	va_b (:,:) = va_b (:,:) + za1 * va_e (:,:)
843	ELSE ! Sum transports
844	ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:) * hu_e (:,:)
845	va_b (:,:) = va_b (:,:) + za1 * va_e (:,:) * hv_e (:,:)
846	ENDIF
847	! ! Sum sea level
848	ssha(:,:) = ssha(:,:) + za1 * ssha_e(:,:)
849	! ! ==================== !
850	END DO ! end loop !
851	! ! ==================== !
852	! -----------------------------------------------------------------------------
853	! Phase 3. update the general trend with the barotropic trend
854	! -----------------------------------------------------------------------------
855	!
856	! At this stage ssha holds a time averaged value
857	! ! Sea Surface Height at u-,v- and f-points
858	IF( lk_vvl ) THEN ! (required only in key_vvl case)
859	DO jj = 1, jpjm1
860	DO ji = 1, jpim1 ! NO Vector Opt.
861	zsshu_a(ji,jj) = z1_2 * umask(ji,jj,1) / ( e1u(ji ,jj) * e2u(ji ,jj) ) &
862	& * ( e1t(ji ,jj) * e2t(ji ,jj) * ssha(ji ,jj) &
863	& + e1t(ji+1,jj) * e2t(ji+1,jj) * ssha(ji+1,jj) )
864	zsshv_a(ji,jj) = z1_2 * vmask(ji,jj,1) / ( e1v(ji,jj ) * e2v(ji,jj ) ) &
865	& * ( e1t(ji,jj ) * e2t(ji,jj ) * ssha(ji,jj ) &
866	& + e1t(ji,jj+1) * e2t(ji,jj+1) * ssha(ji,jj+1) )
867	END DO
868	END DO
869	CALL lbc_lnk( zsshu_a, 'U', 1._wp ) ; CALL lbc_lnk( zsshv_a, 'V', 1._wp ) ! Boundary conditions
870	ENDIF
871	!
872	! Set advection velocity correction:
873	IF (((kt==nit000).AND.(neuler==0)).OR.(.NOT.ln_bt_fw)) THEN
874	un_adv(:,:) = zu_sum(:,:)*hur(:,:)
875	vn_adv(:,:) = zv_sum(:,:)*hvr(:,:)
876	ELSE
877	un_adv(:,:) = z1_2 * ( ub2_b(:,:) + zu_sum(:,:)) * hur(:,:)
878	vn_adv(:,:) = z1_2 * ( vb2_b(:,:) + zv_sum(:,:)) * hvr(:,:)
879	END IF
880
881	IF (ln_bt_fw) THEN ! Save integrated transport for next computation
882	ub2_b(:,:) = zu_sum(:,:)
883	vb2_b(:,:) = zv_sum(:,:)
884	ENDIF
885	!
886	! Update barotropic trend:
887	IF (( ln_dynadv_vec ).OR. (.NOT. lk_vvl)) THEN
888	DO jk=1,jpkm1
889	ua(:,:,jk) = ua(:,:,jk) + ( ua_b(:,:) - ub_b(:,:) ) * z1_2dt_b
890	va(:,:,jk) = va(:,:,jk) + ( va_b(:,:) - vb_b(:,:) ) * z1_2dt_b
891	END DO
892	ELSE
893	hu_e (:,:) = umask(:,:,1) / ( zhur_b(:,:) + 1._wp - umask(:,:,1) )
894	hv_e (:,:) = vmask(:,:,1) / ( zhvr_b(:,:) + 1._wp - vmask(:,:,1) )
895	DO jk=1,jpkm1
896	ua(:,:,jk) = ua(:,:,jk) + hur(:,:) * ( ua_b(:,:) - ub_b(:,:) * hu_e(:,:) ) * z1_2dt_b
897	va(:,:,jk) = va(:,:,jk) + hvr(:,:) * ( va_b(:,:) - vb_b(:,:) * hv_e(:,:) ) * z1_2dt_b
898	END DO
899	! Save barotropic velocities not transport:
900	ua_b (:,:) = ua_b(:,:) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - umask(:,:,1) )
901	va_b (:,:) = va_b(:,:) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - vmask(:,:,1) )
902	ENDIF
903	!
904	DO jk = 1, jpkm1
905	! Correct velocities:
906	un(:,:,jk) = ( un(:,:,jk) + un_adv(:,:) - un_b(:,:) )*umask(:,:,jk)
907	vn(:,:,jk) = ( vn(:,:,jk) + vn_adv(:,:) - vn_b(:,:) )*vmask(:,:,jk)
908	!
909	END DO
910	!
911	! !* write time-spliting arrays in the restart
912	IF(lrst_oce .AND.ln_bt_fw) CALL ts_rst( kt, 'WRITE' )
913	!
914	CALL wrk_dealloc( jpi, jpj, zsshp2_e, zhdiv )
915	CALL wrk_dealloc( jpi, jpj, zu_trd, zv_trd, zun_e, zvn_e )
916	CALL wrk_dealloc( jpi, jpj, zwx, zwy, zu_sum, zv_sum, zssh_frc, zu_frc, zv_frc )
917	CALL wrk_dealloc( jpi, jpj, zhup2_e, zhvp2_e, zhust_e, zhvst_e )
918	CALL wrk_dealloc( jpi, jpj, zhur_b, zhvr_b )
919	CALL wrk_dealloc( jpi, jpj, zsshu_a, zsshv_a )
920	CALL wrk_dealloc( jpi, jpj, zht, zhf )
921	!
922	IF( nn_timing == 1 ) CALL timing_stop('dyn_spg_ts')
923	!
924	END SUBROUTINE dyn_spg_ts
925
926	SUBROUTINE ts_wgt( ll_av, ll_fw, jpit, zwgt1, zwgt2)
927	!!---------------------------------------------------------------------
928	!! * ROUTINE ts_wgt *
929	!!
930	!! ** Purpose : Set time-splitting weights for temporal averaging (or not)
931	!!----------------------------------------------------------------------
932	LOGICAL, INTENT(in) :: ll_av ! temporal averaging=.true.
933	LOGICAL, INTENT(in) :: ll_fw ! forward time splitting =.true.
934	INTEGER, INTENT(inout) :: jpit ! cycle length
935	REAL(wp), DIMENSION(3*nn_baro), INTENT(inout) :: zwgt1, & ! Primary weights
936	zwgt2 ! Secondary weights
937
938	INTEGER :: jic, jn, ji ! temporary integers
939	REAL(wp) :: za1, za2
940	!!----------------------------------------------------------------------
941
942	zwgt1(:) = 0._wp
943	zwgt2(:) = 0._wp
944
945	! Set time index when averaged value is requested
946	IF (ll_fw) THEN
947	jic = nn_baro
948	ELSE
949	jic = 2 * nn_baro
950	ENDIF
951
952	! Set primary weights:
953	IF (ll_av) THEN
954	! Define simple boxcar window for primary weights
955	! (width = nn_baro, centered around jic)
956	SELECT CASE ( nn_bt_flt )
957	CASE( 0 ) ! No averaging
958	zwgt1(jic) = 1._wp
959	jpit = jic
960
961	CASE( 1 ) ! Boxcar, width = nn_baro
962	DO jn = 1, 3*nn_baro
963	za1 = ABS(float(jn-jic))/float(nn_baro)
964	IF (za1 < 0.5_wp) THEN
965	zwgt1(jn) = 1._wp
966	jpit = jn
967	ENDIF
968	ENDDO
969
970	CASE( 2 ) ! Boxcar, width = 2 * nn_baro
971	DO jn = 1, 3*nn_baro
972	za1 = ABS(float(jn-jic))/float(nn_baro)
973	IF (za1 < 1._wp) THEN
974	zwgt1(jn) = 1._wp
975	jpit = jn
976	ENDIF
977	ENDDO
978	CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_bt_flt' )
979	END SELECT
980
981	ELSE ! No time averaging
982	zwgt1(jic) = 1._wp
983	jpit = jic
984	ENDIF
985
986	! Set secondary weights
987	DO jn = 1, jpit
988	DO ji = jn, jpit
989	zwgt2(jn) = zwgt2(jn) + zwgt1(ji)
990	END DO
991	END DO
992
993	! Normalize weigths:
994	za1 = 1._wp / SUM(zwgt1(1:jpit))
995	za2 = 1._wp / SUM(zwgt2(1:jpit))
996	DO jn = 1, jpit
997	zwgt1(jn) = zwgt1(jn) * za1
998	zwgt2(jn) = zwgt2(jn) * za2
999	END DO
1000	!
1001	END SUBROUTINE ts_wgt
1002
1003	SUBROUTINE ts_rst( kt, cdrw )
1004	!!---------------------------------------------------------------------
1005	!! * ROUTINE ts_rst *
1006	!!
1007	!! ** Purpose : Read or write time-splitting arrays in restart file
1008	!!----------------------------------------------------------------------
1009	INTEGER , INTENT(in) :: kt ! ocean time-step
1010	CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag
1011	!
1012	!!----------------------------------------------------------------------
1013	!
1014	IF( TRIM(cdrw) == 'READ' ) THEN
1015	CALL iom_get( numror, jpdom_autoglo, 'ub2_b' , ub2_b (:,:) )
1016	CALL iom_get( numror, jpdom_autoglo, 'vb2_b' , vb2_b (:,:) )
1017	IF( .NOT.ln_bt_av .AND. iom_varid( numror, 'sshbb_e', ldstop = .FALSE. ) > 0) THEN
1018	CALL iom_get( numror, jpdom_autoglo, 'sshbb_e' , sshbb_e(:,:) )
1019	CALL iom_get( numror, jpdom_autoglo, 'ubb_e' , ubb_e(:,:) )
1020	CALL iom_get( numror, jpdom_autoglo, 'vbb_e' , vbb_e(:,:) )
1021	CALL iom_get( numror, jpdom_autoglo, 'sshb_e' , sshb_e(:,:) )
1022	CALL iom_get( numror, jpdom_autoglo, 'ub_e' , ub_e(:,:) )
1023	CALL iom_get( numror, jpdom_autoglo, 'vb_e' , vb_e(:,:) )
1024	ELSE
1025	sshbb_e = sshn_b ! ACC GUESS WORK
1026	ubb_e = ub_b
1027	vbb_e = vb_b
1028	sshb_e = sshn_b
1029	ub_e = ub_b
1030	vb_e = vb_b
1031	ENDIF
1032	!
1033	ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN
1034	CALL iom_rstput( kt, nitrst, numrow, 'ub2_b' , ub2_b (:,:) )
1035	CALL iom_rstput( kt, nitrst, numrow, 'vb2_b' , vb2_b (:,:) )
1036	!
1037	IF (.NOT.ln_bt_av) THEN
1038	CALL iom_rstput( kt, nitrst, numrow, 'sshbb_e' , sshbb_e(:,:) )
1039	CALL iom_rstput( kt, nitrst, numrow, 'ubb_e' , ubb_e(:,:) )
1040	CALL iom_rstput( kt, nitrst, numrow, 'vbb_e' , vbb_e(:,:) )
1041	CALL iom_rstput( kt, nitrst, numrow, 'sshb_e' , sshb_e(:,:) )
1042	CALL iom_rstput( kt, nitrst, numrow, 'ub_e' , ub_e(:,:) )
1043	CALL iom_rstput( kt, nitrst, numrow, 'vb_e' , vb_e(:,:) )
1044	ENDIF
1045	ENDIF
1046	!
1047	END SUBROUTINE ts_rst
1048
1049	SUBROUTINE dyn_spg_ts_init( kt )
1050	!!---------------------------------------------------------------------
1051	!! * ROUTINE dyn_spg_ts_init *
1052	!!
1053	!! ** Purpose : Set time splitting options
1054	!!----------------------------------------------------------------------
1055	INTEGER , INTENT(in) :: kt ! ocean time-step
1056	!
1057	INTEGER :: ji ,jj, jk
1058	REAL(wp) :: zxr2, zyr2, zcmax
1059	REAL(wp), POINTER, DIMENSION(:,:) :: zcu, zht
1060	!!
1061	! NAMELIST/namsplit/ ln_bt_fw, ln_bt_av, ln_bt_nn_auto, &
1062	! & nn_baro, rn_bt_cmax, nn_bt_flt
1063	!!----------------------------------------------------------------------
1064	! REWIND( numnam ) !* Namelist namsplit: split-explicit free surface
1065	! READ ( numnam, namsplit )
1066	! ! Max courant number for ext. grav. waves
1067	!
1068	CALL wrk_alloc( jpi, jpj, zcu, zht )
1069	!
1070	! JC: Simplification needed below: define ht_0 even when volume is fixed
1071	IF (lk_vvl) THEN
1072	zht(:,:) = ht_0(:,:) * tmask(:,:,1)
1073	ELSE
1074	zht(:,:) = 0.
1075	DO jk = 1, jpkm1
1076	zht(:,:) = zht(:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk)
1077	END DO
1078	ENDIF
1079
1080	DO jj = 1, jpj
1081	DO ji =1, jpi
1082	zxr2 = 1./(e1t(ji,jj)*e1t(ji,jj))
1083	zyr2 = 1./(e2t(ji,jj)*e2t(ji,jj))
1084	zcu(ji,jj) = sqrt(gravzht(ji,jj)(zxr2 + zyr2) )
1085	END DO
1086	END DO
1087
1088	zcmax = MAXVAL(zcu(:,:))
1089	IF( lk_mpp ) CALL mpp_max( zcmax )
1090
1091	! Estimate number of iterations to satisfy a max courant number=0.8
1092	IF (ln_bt_nn_auto) nn_baro = CEILING( rdt / rn_bt_cmax * zcmax)
1093
1094	rdtbt = rdt / FLOAT(nn_baro)
1095	zcmax = zcmax * rdtbt
1096	! Print results
1097	IF(lwp) WRITE(numout,*)
1098	IF(lwp) WRITE(numout,*) 'dyn_spg_ts : split-explicit free surface'
1099	IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
1100	IF( ln_bt_nn_auto ) THEN
1101	IF(lwp) WRITE(numout,*) ' ln_ts_nn_auto=.true. Automatically set nn_baro '
1102	IF(lwp) WRITE(numout,*) ' Max. courant number allowed: ', rn_bt_cmax
1103	ELSE
1104	IF(lwp) WRITE(numout,*) ' ln_ts_nn_auto=.false.: Use nn_baro in namelist '
1105	ENDIF
1106	IF(lwp) WRITE(numout,*) ' nn_baro = ', nn_baro
1107	IF(lwp) WRITE(numout,*) ' Barotropic time step [s] is :', rdtbt
1108	IF(lwp) WRITE(numout,*) ' Maximum Courant number is :', zcmax
1109
1110	IF(ln_bt_av) THEN
1111	IF(lwp) WRITE(numout,*) ' ln_bt_av=.true. => Time averaging over nn_baro time steps is on '
1112	ELSE
1113	IF(lwp) WRITE(numout,*) ' ln_bt_av=.false. => No time averaging of barotropic variables '
1114	ENDIF
1115	!
1116	!
1117	IF(ln_bt_fw) THEN
1118	IF(lwp) WRITE(numout,*) ' ln_bt_fw=.true. => Forward integration of barotropic variables '
1119	ELSE
1120	IF(lwp) WRITE(numout,*) ' ln_bt_fw =.false.=> Centered integration of barotropic variables '
1121	ENDIF
1122	!
1123	IF(lwp) WRITE(numout,*) ' Time filter choice, nn_bt_flt: ', nn_bt_flt
1124	SELECT CASE ( nn_bt_flt )
1125	CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' Dirac'
1126	CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' Boxcar: width = nn_baro'
1127	CASE( 2 ) ; IF(lwp) WRITE(numout,) ' Boxcar: width = 2nn_baro'
1128	CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_bt_flt: should 0,1,2' )
1129	END SELECT
1130	!
1131	IF ((.NOT.ln_bt_av).AND.(.NOT.ln_bt_fw)) THEN
1132	CALL ctl_stop( 'dynspg_ts ERROR: No time averaging => only forward integration is possible' )
1133	ENDIF
1134	IF ( zcmax>0.9_wp ) THEN
1135	CALL ctl_stop( 'dynspg_ts ERROR: Maximum Courant number is greater than 0.9: Inc. nn_baro !' )
1136	ENDIF
1137	!
1138	CALL wrk_dealloc( jpi, jpj, zcu, zht )
1139	!
1140	END SUBROUTINE dyn_spg_ts_init
1141
1142	#else
1143	!!---------------------------------------------------------------------------
1144	!! Default case : Empty module No standard free surface constant volume
1145	!!---------------------------------------------------------------------------
1146
1147	USE par_kind
1148	LOGICAL, PUBLIC, PARAMETER :: ln_bt_fw=.FALSE. ! Forward integration of barotropic sub-stepping
1149	CONTAINS
1150	INTEGER FUNCTION dyn_spg_ts_alloc() ! Dummy function
1151	dyn_spg_ts_alloc = 0
1152	END FUNCTION dyn_spg_ts_alloc
1153	SUBROUTINE dyn_spg_ts( kt ) ! Empty routine
1154	INTEGER, INTENT(in) :: kt
1155	WRITE(,) 'dyn_spg_ts: You should not have seen this print! error?', kt
1156	END SUBROUTINE dyn_spg_ts
1157	SUBROUTINE ts_rst( kt, cdrw ) ! Empty routine
1158	INTEGER , INTENT(in) :: kt ! ocean time-step
1159	CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag
1160	WRITE(,) 'ts_rst : You should not have seen this print! error?', kt, cdrw
1161	END SUBROUTINE ts_rst
1162	SUBROUTINE dyn_spg_ts_init( kt ) ! Empty routine
1163	INTEGER , INTENT(in) :: kt ! ocean time-step
1164	WRITE(,) 'dyn_spg_ts_init : You should not have seen this print! error?', kt
1165	END SUBROUTINE dyn_spg_ts_init
1166	#endif
1167
1168	!!======================================================================
1169	END MODULE dynspg_ts

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

Download in other formats: