Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

dynspg_ts.F90 @ 5605

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

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

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

Context Navigation

source: branches/NERC/dev_r5589_marine_glacier_termini/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 5605

Download in other formats: