Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

dynspg_ts.F90 @ 12358

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

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

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

Context Navigation

source: branches/UKMO/AMM15_v3_6_STABLE_package_collate_PS44/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 12358

Download in other formats: