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/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

Context Navigation

source: branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 5969

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

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

Download in other formats: