Context Navigation

dynspg_ts.F90 @ 11134

Last change on this file since 11134 was 11134, checked in by jcastill, 5 years ago
Full set of changes as in the original branch
File size: 56.1 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 dynspg_oce ! surface pressure gradient variables
26	USE phycst ! physical constants
27	USE dynvor ! vorticity term
28	USE bdy_par ! for lk_bdy
29	USE bdytides ! open boundary condition data
30	USE bdydyn2d ! open boundary conditions on barotropic variables
31	USE sbctide ! tides
32	USE updtide ! tide potential
33	USE lib_mpp ! distributed memory computing library
34	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
35	USE prtctl ! Print control
36	USE in_out_manager ! I/O manager
37	USE iom ! IOM library
38	USE restart ! only for lrst_oce
39	USE zdf_oce ! Bottom friction coefts
40	USE wrk_nemo ! Memory Allocation
41	USE timing ! Timing
42	USE sbcapr ! surface boundary condition: atmospheric pressure
43	USE diatmb ! Top,middle,bottom output
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) :: 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	! !* barotropic Coriolis trends (vorticity scheme dependent)
336	! ! --------------------------------------------------------
337	zwx(:,:) = un_b(:,:) * hu(:,:) * e2u(:,:) ! now fluxes
338	zwy(:,:) = vn_b(:,:) * hv(:,:) * e1v(:,:)
339	!
340	IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN ! energy conserving or mixed scheme
341	DO jj = 2, jpjm1
342	DO ji = fs_2, fs_jpim1 ! vector opt.
343	zy1 = ( zwy(ji,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj)
344	zy2 = ( zwy(ji,jj ) + zwy(ji+1,jj ) ) / e1u(ji,jj)
345	zx1 = ( zwx(ji-1,jj) + zwx(ji-1,jj+1) ) / e2v(ji,jj)
346	zx2 = ( zwx(ji ,jj) + zwx(ji ,jj+1) ) / e2v(ji,jj)
347	! energy conserving formulation for planetary vorticity term
348	zu_trd(ji,jj) = z1_4 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
349	zv_trd(ji,jj) =-z1_4 * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 )
350	END DO
351	END DO
352	!
353	ELSEIF ( ln_dynvor_ens ) THEN ! enstrophy conserving scheme
354	DO jj = 2, jpjm1
355	DO ji = fs_2, fs_jpim1 ! vector opt.
356	zy1 = z1_8 * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) &
357	& + zwy(ji ,jj ) + zwy(ji+1,jj ) ) / e1u(ji,jj)
358	zx1 = - z1_8 * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) &
359	& + zwx(ji ,jj ) + zwx(ji ,jj+1) ) / e2v(ji,jj)
360	zu_trd(ji,jj) = zy1 * ( zwz(ji ,jj-1) + zwz(ji,jj) )
361	zv_trd(ji,jj) = zx1 * ( zwz(ji-1,jj ) + zwz(ji,jj) )
362	END DO
363	END DO
364	!
365	ELSEIF ( ln_dynvor_een .or. ln_dynvor_een_old ) THEN ! enstrophy and energy conserving scheme
366	DO jj = 2, jpjm1
367	DO ji = fs_2, fs_jpim1 ! vector opt.
368	zu_trd(ji,jj) = + z1_12 / e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) &
369	& + ftnw(ji+1,jj) * zwy(ji+1,jj ) &
370	& + ftse(ji,jj ) * zwy(ji ,jj-1) &
371	& + ftsw(ji+1,jj) * zwy(ji+1,jj-1) )
372	zv_trd(ji,jj) = - z1_12 / e2v(ji,jj) * ( ftsw(ji,jj+1) * zwx(ji-1,jj+1) &
373	& + ftse(ji,jj+1) * zwx(ji ,jj+1) &
374	& + ftnw(ji,jj ) * zwx(ji-1,jj ) &
375	& + ftne(ji,jj ) * zwx(ji ,jj ) )
376	END DO
377	END DO
378	!
379	ENDIF
380	!
381	! !* Right-Hand-Side of the barotropic momentum equation
382	! ! ----------------------------------------------------
383	IF( lk_vvl ) THEN ! Variable volume : remove surface pressure gradient
384	DO jj = 2, jpjm1
385	DO ji = fs_2, fs_jpim1 ! vector opt.
386	zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) / e1u(ji,jj)
387	zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji ,jj+1) - sshn(ji ,jj ) ) / e2v(ji,jj)
388	END DO
389	END DO
390	ENDIF
391
392	DO jj = 2, jpjm1 ! Remove coriolis term (and possibly spg) from barotropic trend
393	DO ji = fs_2, fs_jpim1
394	zu_frc(ji,jj) = zu_frc(ji,jj) - zu_trd(ji,jj) * umask(ji,jj,1)
395	zv_frc(ji,jj) = zv_frc(ji,jj) - zv_trd(ji,jj) * vmask(ji,jj,1)
396	END DO
397	END DO
398	!
399	! ! Add bottom stress contribution from baroclinic velocities:
400	IF (ln_bt_fw) THEN
401	DO jj = 2, jpjm1
402	DO ji = fs_2, fs_jpim1 ! vector opt.
403	ikbu = mbku(ji,jj)
404	ikbv = mbkv(ji,jj)
405	zwx(ji,jj) = un(ji,jj,ikbu) - un_b(ji,jj) ! NOW bottom baroclinic velocities
406	zwy(ji,jj) = vn(ji,jj,ikbv) - vn_b(ji,jj)
407	END DO
408	END DO
409	ELSE
410	DO jj = 2, jpjm1
411	DO ji = fs_2, fs_jpim1 ! vector opt.
412	ikbu = mbku(ji,jj)
413	ikbv = mbkv(ji,jj)
414	zwx(ji,jj) = ub(ji,jj,ikbu) - ub_b(ji,jj) ! BEFORE bottom baroclinic velocities
415	zwy(ji,jj) = vb(ji,jj,ikbv) - vb_b(ji,jj)
416	END DO
417	END DO
418	ENDIF
419	!
420	! Note that the "unclipped" bottom friction parameter is used even with explicit drag
421	zu_frc(:,:) = zu_frc(:,:) + hur(:,:) * bfrua(:,:) * zwx(:,:)
422	zv_frc(:,:) = zv_frc(:,:) + hvr(:,:) * bfrva(:,:) * zwy(:,:)
423	!
424	IF (ln_bt_fw) THEN ! Add wind forcing
425	zu_frc(:,:) = zu_frc(:,:) + zraur * utau(:,:) * hur(:,:)
426	zv_frc(:,:) = zv_frc(:,:) + zraur * vtau(:,:) * hvr(:,:)
427	ELSE
428	zu_frc(:,:) = zu_frc(:,:) + zraur * z1_2 * ( utau_b(:,:) + utau(:,:) ) * hur(:,:)
429	zv_frc(:,:) = zv_frc(:,:) + zraur * z1_2 * ( vtau_b(:,:) + vtau(:,:) ) * hvr(:,:)
430	ENDIF
431	!
432	IF ( ln_apr_dyn ) THEN ! Add atm pressure forcing
433	IF (ln_bt_fw) THEN
434	DO jj = 2, jpjm1
435	DO ji = fs_2, fs_jpim1 ! vector opt.
436	zu_spg = grav * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) ) /e1u(ji,jj)
437	zv_spg = grav * ( ssh_ib (ji ,jj+1) - ssh_ib (ji,jj) ) /e2v(ji,jj)
438	zu_frc(ji,jj) = zu_frc(ji,jj) + zu_spg
439	zv_frc(ji,jj) = zv_frc(ji,jj) + zv_spg
440	END DO
441	END DO
442	ELSE
443	DO jj = 2, jpjm1
444	DO ji = fs_2, fs_jpim1 ! vector opt.
445	zu_spg = grav * z1_2 * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) &
446	& + ssh_ibb(ji+1,jj ) - ssh_ibb(ji,jj) ) /e1u(ji,jj)
447	zv_spg = grav * z1_2 * ( ssh_ib (ji ,jj+1) - ssh_ib (ji,jj) &
448	& + ssh_ibb(ji ,jj+1) - ssh_ibb(ji,jj) ) /e2v(ji,jj)
449	zu_frc(ji,jj) = zu_frc(ji,jj) + zu_spg
450	zv_frc(ji,jj) = zv_frc(ji,jj) + zv_spg
451	END DO
452	END DO
453	ENDIF
454	ENDIF
455	! !* Right-Hand-Side of the barotropic ssh equation
456	! ! -----------------------------------------------
457	! ! Surface net water flux and rivers
458	IF (ln_bt_fw) THEN
459	zssh_frc(:,:) = zraur * ( emp(:,:) - rnf(:,:) + fwfisf(:,:) )
460	ELSE
461	zssh_frc(:,:) = zraur * z1_2 * ( emp(:,:) + emp_b(:,:) - rnf(:,:) - rnf_b(:,:) &
462	& + fwfisf(:,:) + fwfisf_b(:,:) )
463	ENDIF
464	#if defined key_asminc
465	! ! Include the IAU weighted SSH increment
466	IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN
467	zssh_frc(:,:) = zssh_frc(:,:) - ssh_iau(:,:)
468	ENDIF
469	#endif
470	! !* Fill boundary data arrays with AGRIF
471	! ! -------------------------------------
472	#if defined key_agrif
473	IF( .NOT.Agrif_Root() ) CALL agrif_dta_ts( kt )
474	#endif
475	!
476	! -----------------------------------------------------------------------
477	! Phase 2 : Integration of the barotropic equations
478	! -----------------------------------------------------------------------
479	!
480	! ! ==================== !
481	! ! Initialisations !
482	! ! ==================== !
483	! Initialize barotropic variables:
484	IF( ll_init )THEN
485	sshbb_e(:,:) = 0._wp
486	ubb_e (:,:) = 0._wp
487	vbb_e (:,:) = 0._wp
488	sshb_e (:,:) = 0._wp
489	ub_e (:,:) = 0._wp
490	vb_e (:,:) = 0._wp
491	ENDIF
492	!
493	IF (ln_bt_fw) THEN ! FORWARD integration: start from NOW fields
494	sshn_e(:,:) = sshn (:,:)
495	zun_e (:,:) = un_b (:,:)
496	zvn_e (:,:) = vn_b (:,:)
497	!
498	hu_e (:,:) = hu (:,:)
499	hv_e (:,:) = hv (:,:)
500	hur_e (:,:) = hur (:,:)
501	hvr_e (:,:) = hvr (:,:)
502	ELSE ! CENTRED integration: start from BEFORE fields
503	sshn_e(:,:) = sshb (:,:)
504	zun_e (:,:) = ub_b (:,:)
505	zvn_e (:,:) = vb_b (:,:)
506	!
507	hu_e (:,:) = hu_b (:,:)
508	hv_e (:,:) = hv_b (:,:)
509	hur_e (:,:) = hur_b(:,:)
510	hvr_e (:,:) = hvr_b(:,:)
511	ENDIF
512	!
513	!
514	!
515	! Initialize sums:
516	ua_b (:,:) = 0._wp ! After barotropic velocities (or transport if flux form)
517	va_b (:,:) = 0._wp
518	ssha (:,:) = 0._wp ! Sum for after averaged sea level
519	zu_sum(:,:) = 0._wp ! Sum for now transport issued from ts loop
520	zv_sum(:,:) = 0._wp
521	! ! ==================== !
522	DO jn = 1, icycle ! sub-time-step loop !
523	! ! ==================== !
524	! !* Update the forcing (BDY and tides)
525	! ! ------------------
526	! Update only tidal forcing at open boundaries
527	#if defined key_tide
528	IF ( lk_bdy .AND. lk_tide ) CALL bdy_dta_tides( kt, kit=jn, time_offset=(noffset+1) )
529	IF ( ln_tide_pot .AND. lk_tide ) CALL upd_tide( kt, kit=jn, time_offset=noffset )
530	#endif
531	!
532	! Set extrapolation coefficients for predictor step:
533	IF ((jn<3).AND.ll_init) THEN ! Forward
534	za1 = 1._wp
535	za2 = 0._wp
536	za3 = 0._wp
537	ELSE ! AB3-AM4 Coefficients: bet=0.281105
538	za1 = 1.781105_wp ! za1 = 3/2 + bet
539	za2 = -1.06221_wp ! za2 = -(1/2 + 2*bet)
540	za3 = 0.281105_wp ! za3 = bet
541	ENDIF
542
543	! Extrapolate barotropic velocities at step jit+0.5:
544	ua_e(:,:) = za1 * zun_e(:,:) + za2 * ub_e(:,:) + za3 * ubb_e(:,:)
545	va_e(:,:) = za1 * zvn_e(:,:) + za2 * vb_e(:,:) + za3 * vbb_e(:,:)
546
547	IF( lk_vvl ) THEN !* Update ocean depth (variable volume case only)
548	! ! ------------------
549	! Extrapolate Sea Level at step jit+0.5:
550	zsshp2_e(:,:) = za1 * sshn_e(:,:) + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:)
551	!
552	DO jj = 2, jpjm1 ! Sea Surface Height at u- & v-points
553	DO ji = 2, fs_jpim1 ! Vector opt.
554	zwx(ji,jj) = z1_2 * umask(ji,jj,1) * r1_e12u(ji,jj) &
555	& * ( e12t(ji ,jj) * zsshp2_e(ji ,jj) &
556	& + e12t(ji+1,jj) * zsshp2_e(ji+1,jj) )
557	zwy(ji,jj) = z1_2 * vmask(ji,jj,1) * r1_e12v(ji,jj) &
558	& * ( e12t(ji,jj ) * zsshp2_e(ji,jj ) &
559	& + e12t(ji,jj+1) * zsshp2_e(ji,jj+1) )
560	END DO
561	END DO
562	CALL lbc_lnk_multi( zwx, 'U', 1._wp, zwy, 'V', 1._wp )
563	!
564	zhup2_e (:,:) = hu_0(:,:) + zwx(:,:) ! Ocean depth at U- and V-points
565	zhvp2_e (:,:) = hv_0(:,:) + zwy(:,:)
566	ELSE
567	zhup2_e (:,:) = hu(:,:)
568	zhvp2_e (:,:) = hv(:,:)
569	ENDIF
570	! !* after ssh
571	! ! -----------
572	! One should enforce volume conservation at open boundaries here
573	! considering fluxes below:
574	!
575	zwx(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:) ! fluxes at jn+0.5
576	zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:)
577	!
578	#if defined key_agrif
579	! Set fluxes during predictor step to ensure
580	! volume conservation
581	IF( (.NOT.Agrif_Root()).AND.ln_bt_fw ) THEN
582	IF((nbondi == -1).OR.(nbondi == 2)) THEN
583	DO jj=1,jpj
584	zwx(2,jj) = ubdy_w(jj) * e2u(2,jj)
585	END DO
586	ENDIF
587	IF((nbondi == 1).OR.(nbondi == 2)) THEN
588	DO jj=1,jpj
589	zwx(nlci-2,jj) = ubdy_e(jj) * e2u(nlci-2,jj)
590	END DO
591	ENDIF
592	IF((nbondj == -1).OR.(nbondj == 2)) THEN
593	DO ji=1,jpi
594	zwy(ji,2) = vbdy_s(ji) * e1v(ji,2)
595	END DO
596	ENDIF
597	IF((nbondj == 1).OR.(nbondj == 2)) THEN
598	DO ji=1,jpi
599	zwy(ji,nlcj-2) = vbdy_n(ji) * e1v(ji,nlcj-2)
600	END DO
601	ENDIF
602	ENDIF
603	#endif
604	!
605	! Sum over sub-time-steps to compute advective velocities
606	za2 = wgtbtp2(jn)
607	zu_sum (:,:) = zu_sum (:,:) + za2 * zwx (:,:) / e2u (:,:)
608	zv_sum (:,:) = zv_sum (:,:) + za2 * zwy (:,:) / e1v (:,:)
609	!
610	! Set next sea level:
611	DO jj = 2, jpjm1
612	DO ji = fs_2, fs_jpim1 ! vector opt.
613	zhdiv(ji,jj) = ( zwx(ji,jj) - zwx(ji-1,jj) &
614	& + zwy(ji,jj) - zwy(ji,jj-1) ) * r1_e12t(ji,jj)
615	END DO
616	END DO
617	ssha_e(:,:) = ( sshn_e(:,:) - rdtbt * ( zssh_frc(:,:) + zhdiv(:,:) ) ) * tmask(:,:,1)
618	CALL lbc_lnk( ssha_e, 'T', 1._wp )
619
620	#if defined key_bdy
621	! Duplicate sea level across open boundaries (this is only cosmetic if lk_vvl=.false.)
622	IF (lk_bdy) CALL bdy_ssh( ssha_e )
623	#endif
624	#if defined key_agrif
625	IF( .NOT.Agrif_Root() ) CALL agrif_ssh_ts( jn )
626	#endif
627	!
628	! Sea Surface Height at u-,v-points (vvl case only)
629	IF ( lk_vvl ) THEN
630	DO jj = 2, jpjm1
631	DO ji = 2, jpim1 ! NO Vector Opt.
632	zsshu_a(ji,jj) = z1_2 * umask(ji,jj,1) * r1_e12u(ji,jj) &
633	& * ( e12t(ji ,jj ) * ssha_e(ji ,jj ) &
634	& + e12t(ji+1,jj ) * ssha_e(ji+1,jj ) )
635	zsshv_a(ji,jj) = z1_2 * vmask(ji,jj,1) * r1_e12v(ji,jj) &
636	& * ( e12t(ji ,jj ) * ssha_e(ji ,jj ) &
637	& + e12t(ji ,jj+1) * ssha_e(ji ,jj+1) )
638	END DO
639	END DO
640	CALL lbc_lnk_multi( zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp )
641	ENDIF
642	!
643	! Half-step back interpolation of SSH for surface pressure computation:
644	!----------------------------------------------------------------------
645	IF ((jn==1).AND.ll_init) THEN
646	za0=1._wp ! Forward-backward
647	za1=0._wp
648	za2=0._wp
649	za3=0._wp
650	ELSEIF ((jn==2).AND.ll_init) THEN ! AB2-AM3 Coefficients; bet=0 ; gam=-1/6 ; eps=1/12
651	za0= 1.0833333333333_wp ! za0 = 1-gam-eps
652	za1=-0.1666666666666_wp ! za1 = gam
653	za2= 0.0833333333333_wp ! za2 = eps
654	za3= 0._wp
655	ELSE ! AB3-AM4 Coefficients; bet=0.281105 ; eps=0.013 ; gam=0.0880
656	za0=0.614_wp ! za0 = 1/2 + gam + 2*eps
657	za1=0.285_wp ! za1 = 1/2 - 2gam - 3eps
658	za2=0.088_wp ! za2 = gam
659	za3=0.013_wp ! za3 = eps
660	ENDIF
661
662	zsshp2_e(:,:) = za0 * ssha_e(:,:) + za1 * sshn_e (:,:) &
663	& + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:)
664
665	!
666	! Compute associated depths at U and V points:
667	IF ( lk_vvl.AND.(.NOT.ln_dynadv_vec) ) THEN
668	!
669	DO jj = 2, jpjm1
670	DO ji = 2, jpim1
671	zx1 = z1_2 * umask(ji ,jj,1) * r1_e12u(ji ,jj) &
672	& * ( e12t(ji ,jj ) * zsshp2_e(ji ,jj) &
673	& + e12t(ji+1,jj ) * zsshp2_e(ji+1,jj ) )
674	zy1 = z1_2 * vmask(ji ,jj,1) * r1_e12v(ji ,jj ) &
675	& * ( e12t(ji ,jj ) * zsshp2_e(ji ,jj ) &
676	& + e12t(ji ,jj+1) * zsshp2_e(ji ,jj+1) )
677	zhust_e(ji,jj) = hu_0(ji,jj) + zx1
678	zhvst_e(ji,jj) = hv_0(ji,jj) + zy1
679	END DO
680	END DO
681	ENDIF
682	!
683	! Add Coriolis trend:
684	! zwz array below or triads normally depend on sea level with key_vvl and should be updated
685	! at each time step. We however keep them constant here for optimization.
686	! Recall that zwx and zwy arrays hold fluxes at this stage:
687	! zwx(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:) ! fluxes at jn+0.5
688	! zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:)
689	!
690	IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN !== energy conserving or mixed scheme ==!
691	DO jj = 2, jpjm1
692	DO ji = fs_2, fs_jpim1 ! vector opt.
693	zy1 = ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj)
694	zy2 = ( zwy(ji ,jj ) + zwy(ji+1,jj ) ) / e1u(ji,jj)
695	zx1 = ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) ) / e2v(ji,jj)
696	zx2 = ( zwx(ji ,jj ) + zwx(ji ,jj+1) ) / e2v(ji,jj)
697	zu_trd(ji,jj) = z1_4 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
698	zv_trd(ji,jj) =-z1_4 * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 )
699	END DO
700	END DO
701	!
702	ELSEIF ( ln_dynvor_ens ) THEN !== enstrophy conserving scheme ==!
703	DO jj = 2, jpjm1
704	DO ji = fs_2, fs_jpim1 ! vector opt.
705	zy1 = z1_8 * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) &
706	& + zwy(ji ,jj ) + zwy(ji+1,jj ) ) / e1u(ji,jj)
707	zx1 = - z1_8 * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) &
708	& + zwx(ji ,jj ) + zwx(ji ,jj+1) ) / e2v(ji,jj)
709	zu_trd(ji,jj) = zy1 * ( zwz(ji ,jj-1) + zwz(ji,jj) )
710	zv_trd(ji,jj) = zx1 * ( zwz(ji-1,jj ) + zwz(ji,jj) )
711	END DO
712	END DO
713	!
714	ELSEIF ( ln_dynvor_een .or. ln_dynvor_een_old ) THEN !== energy and enstrophy conserving scheme ==!
715	DO jj = 2, jpjm1
716	DO ji = fs_2, fs_jpim1 ! vector opt.
717	zu_trd(ji,jj) = + z1_12 / e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) &
718	& + ftnw(ji+1,jj) * zwy(ji+1,jj ) &
719	& + ftse(ji,jj ) * zwy(ji ,jj-1) &
720	& + ftsw(ji+1,jj) * zwy(ji+1,jj-1) )
721	zv_trd(ji,jj) = - z1_12 / e2v(ji,jj) * ( ftsw(ji,jj+1) * zwx(ji-1,jj+1) &
722	& + ftse(ji,jj+1) * zwx(ji ,jj+1) &
723	& + ftnw(ji,jj ) * zwx(ji-1,jj ) &
724	& + ftne(ji,jj ) * zwx(ji ,jj ) )
725	END DO
726	END DO
727	!
728	ENDIF
729	!
730	! Add tidal astronomical forcing if defined
731	IF ( lk_tide.AND.ln_tide_pot ) THEN
732	DO jj = 2, jpjm1
733	DO ji = fs_2, fs_jpim1 ! vector opt.
734	zu_spg = grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) / e1u(ji,jj)
735	zv_spg = grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) / e2v(ji,jj)
736	zu_trd(ji,jj) = zu_trd(ji,jj) + zu_spg
737	zv_trd(ji,jj) = zv_trd(ji,jj) + zv_spg
738	END DO
739	END DO
740	ENDIF
741	!
742	! Add bottom stresses:
743	zu_trd(:,:) = zu_trd(:,:) + bfrua(:,:) * zun_e(:,:) * hur_e(:,:)
744	zv_trd(:,:) = zv_trd(:,:) + bfrva(:,:) * zvn_e(:,:) * hvr_e(:,:)
745	!
746	! Surface pressure trend:
747	DO jj = 2, jpjm1
748	DO ji = fs_2, fs_jpim1 ! vector opt.
749	! Add surface pressure gradient
750	zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) / e1u(ji,jj)
751	zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) / e2v(ji,jj)
752	zwx(ji,jj) = zu_spg
753	zwy(ji,jj) = zv_spg
754	END DO
755	END DO
756	!
757	! Set next velocities:
758	IF( ln_dynadv_vec .OR. (.NOT. lk_vvl) ) THEN ! Vector form
759	DO jj = 2, jpjm1
760	DO ji = fs_2, fs_jpim1 ! vector opt.
761	ua_e(ji,jj) = ( zun_e(ji,jj) &
762	& + rdtbt * ( zwx(ji,jj) &
763	& + zu_trd(ji,jj) &
764	& + zu_frc(ji,jj) ) &
765	& ) * umask(ji,jj,1)
766
767	va_e(ji,jj) = ( zvn_e(ji,jj) &
768	& + rdtbt * ( zwy(ji,jj) &
769	& + zv_trd(ji,jj) &
770	& + zv_frc(ji,jj) ) &
771	& ) * vmask(ji,jj,1)
772	END DO
773	END DO
774
775	ELSE ! Flux form
776	DO jj = 2, jpjm1
777	DO ji = fs_2, fs_jpim1 ! vector opt.
778
779	zhura = umask(ji,jj,1)/(hu_0(ji,jj) + zsshu_a(ji,jj) + 1._wp - umask(ji,jj,1))
780	zhvra = vmask(ji,jj,1)/(hv_0(ji,jj) + zsshv_a(ji,jj) + 1._wp - vmask(ji,jj,1))
781
782	ua_e(ji,jj) = ( hu_e(ji,jj) * zun_e(ji,jj) &
783	& + rdtbt * ( zhust_e(ji,jj) * zwx(ji,jj) &
784	& + zhup2_e(ji,jj) * zu_trd(ji,jj) &
785	& + hu(ji,jj) * zu_frc(ji,jj) ) &
786	& ) * zhura
787
788	va_e(ji,jj) = ( hv_e(ji,jj) * zvn_e(ji,jj) &
789	& + rdtbt * ( zhvst_e(ji,jj) * zwy(ji,jj) &
790	& + zhvp2_e(ji,jj) * zv_trd(ji,jj) &
791	& + hv(ji,jj) * zv_frc(ji,jj) ) &
792	& ) * zhvra
793	END DO
794	END DO
795	ENDIF
796	!
797	IF( lk_vvl ) THEN !* Update ocean depth (variable volume case only)
798	! ! ----------------------------------------------
799	hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:)
800	hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:)
801	hur_e(:,:) = umask(:,:,1) / ( hu_e(:,:) + 1._wp - umask(:,:,1) )
802	hvr_e(:,:) = vmask(:,:,1) / ( hv_e(:,:) + 1._wp - vmask(:,:,1) )
803	!
804	ENDIF
805	! !* domain lateral boundary
806	! ! -----------------------
807	!
808	CALL lbc_lnk_multi( ua_e, 'U', -1._wp, va_e , 'V', -1._wp )
809
810	#if defined key_bdy
811	! open boundaries
812	IF( lk_bdy ) CALL bdy_dyn2d( jn, ua_e, va_e, zun_e, zvn_e, hur_e, hvr_e, ssha_e )
813	#endif
814	#if defined key_agrif
815	IF( .NOT.Agrif_Root() ) CALL agrif_dyn_ts( jn ) ! Agrif
816	#endif
817	! !* Swap
818	! ! ----
819	ubb_e (:,:) = ub_e (:,:)
820	ub_e (:,:) = zun_e (:,:)
821	zun_e (:,:) = ua_e (:,:)
822	!
823	vbb_e (:,:) = vb_e (:,:)
824	vb_e (:,:) = zvn_e (:,:)
825	zvn_e (:,:) = va_e (:,:)
826	!
827	sshbb_e(:,:) = sshb_e(:,:)
828	sshb_e (:,:) = sshn_e(:,:)
829	sshn_e (:,:) = ssha_e(:,:)
830
831	! !* Sum over whole bt loop
832	! ! ----------------------
833	za1 = wgtbtp1(jn)
834	IF (( ln_dynadv_vec ).OR. (.NOT. lk_vvl)) THEN ! Sum velocities
835	ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:)
836	va_b (:,:) = va_b (:,:) + za1 * va_e (:,:)
837	ELSE ! Sum transports
838	ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:) * hu_e (:,:)
839	va_b (:,:) = va_b (:,:) + za1 * va_e (:,:) * hv_e (:,:)
840	ENDIF
841	! ! Sum sea level
842	ssha(:,:) = ssha(:,:) + za1 * ssha_e(:,:)
843	! ! ==================== !
844	END DO ! end loop !
845	! ! ==================== !
846	! -----------------------------------------------------------------------------
847	! Phase 3. update the general trend with the barotropic trend
848	! -----------------------------------------------------------------------------
849	!
850	! At this stage ssha holds a time averaged value
851	! ! Sea Surface Height at u-,v- and f-points
852	IF( lk_vvl ) THEN ! (required only in key_vvl case)
853	DO jj = 1, jpjm1
854	DO ji = 1, jpim1 ! NO Vector Opt.
855	zsshu_a(ji,jj) = z1_2 * umask(ji,jj,1) * r1_e12u(ji,jj) &
856	& * ( e12t(ji ,jj) * ssha(ji ,jj) &
857	& + e12t(ji+1,jj) * ssha(ji+1,jj) )
858	zsshv_a(ji,jj) = z1_2 * vmask(ji,jj,1) * r1_e12v(ji,jj) &
859	& * ( e12t(ji,jj ) * ssha(ji,jj ) &
860	& + e12t(ji,jj+1) * ssha(ji,jj+1) )
861	END DO
862	END DO
863	CALL lbc_lnk_multi( zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) ! Boundary conditions
864	ENDIF
865	!
866	! Set advection velocity correction:
867	IF (((kt==nit000).AND.(neuler==0)).OR.(.NOT.ln_bt_fw)) THEN
868	un_adv(:,:) = zu_sum(:,:)*hur(:,:)
869	vn_adv(:,:) = zv_sum(:,:)*hvr(:,:)
870	ELSE
871	un_adv(:,:) = z1_2 * ( ub2_b(:,:) + zu_sum(:,:)) * hur(:,:)
872	vn_adv(:,:) = z1_2 * ( vb2_b(:,:) + zv_sum(:,:)) * hvr(:,:)
873	END IF
874
875	IF (ln_bt_fw) THEN ! Save integrated transport for next computation
876	ub2_b(:,:) = zu_sum(:,:)
877	vb2_b(:,:) = zv_sum(:,:)
878	ENDIF
879	!
880	! Update barotropic trend:
881	IF (( ln_dynadv_vec ).OR. (.NOT. lk_vvl)) THEN
882	DO jk=1,jpkm1
883	ua(:,:,jk) = ua(:,:,jk) + ( ua_b(:,:) - ub_b(:,:) ) * z1_2dt_b
884	va(:,:,jk) = va(:,:,jk) + ( va_b(:,:) - vb_b(:,:) ) * z1_2dt_b
885	END DO
886	ELSE
887	DO jk=1,jpkm1
888	ua(:,:,jk) = ua(:,:,jk) + hur(:,:) * ( ua_b(:,:) - ub_b(:,:) * hu_b(:,:) ) * z1_2dt_b
889	va(:,:,jk) = va(:,:,jk) + hvr(:,:) * ( va_b(:,:) - vb_b(:,:) * hv_b(:,:) ) * z1_2dt_b
890	END DO
891	! Save barotropic velocities not transport:
892	ua_b (:,:) = ua_b(:,:) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - umask(:,:,1) )
893	va_b (:,:) = va_b(:,:) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - vmask(:,:,1) )
894	ENDIF
895	!
896	DO jk = 1, jpkm1
897	! Correct velocities:
898	un(:,:,jk) = ( un(:,:,jk) + un_adv(:,:) - un_b(:,:) )*umask(:,:,jk)
899	vn(:,:,jk) = ( vn(:,:,jk) + vn_adv(:,:) - vn_b(:,:) )*vmask(:,:,jk)
900	!
901	END DO
902	!
903	#if defined key_agrif
904	! Save time integrated fluxes during child grid integration
905	! (used to update coarse grid transports at next time step)
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 ( ln_diatmb ) THEN
932	CALL iom_put( "baro_u" , un_bumask(:,:,1)+zmdi(1-umask(:,:,1 ) ) ) ! Barotropic U Velocity
933	CALL iom_put( "baro_v" , vn_bvmask(:,:,1)+zmdi(1-vmask(:,:,1 ) ) ) ! Barotropic V Velocity
934	ENDIF
935	IF( nn_timing == 1 ) CALL timing_stop('dyn_spg_ts')
936	!
937	END SUBROUTINE dyn_spg_ts
938
939	SUBROUTINE ts_wgt( ll_av, ll_fw, jpit, zwgt1, zwgt2)
940	!!---------------------------------------------------------------------
941	!! * ROUTINE ts_wgt *
942	!!
943	!! ** Purpose : Set time-splitting weights for temporal averaging (or not)
944	!!----------------------------------------------------------------------
945	LOGICAL, INTENT(in) :: ll_av ! temporal averaging=.true.
946	LOGICAL, INTENT(in) :: ll_fw ! forward time splitting =.true.
947	INTEGER, INTENT(inout) :: jpit ! cycle length
948	REAL(wp), DIMENSION(3*nn_baro), INTENT(inout) :: zwgt1, & ! Primary weights
949	zwgt2 ! Secondary weights
950
951	INTEGER :: jic, jn, ji ! temporary integers
952	REAL(wp) :: za1, za2
953	!!----------------------------------------------------------------------
954
955	zwgt1(:) = 0._wp
956	zwgt2(:) = 0._wp
957
958	! Set time index when averaged value is requested
959	IF (ll_fw) THEN
960	jic = nn_baro
961	ELSE
962	jic = 2 * nn_baro
963	ENDIF
964
965	! Set primary weights:
966	IF (ll_av) THEN
967	! Define simple boxcar window for primary weights
968	! (width = nn_baro, centered around jic)
969	SELECT CASE ( nn_bt_flt )
970	CASE( 0 ) ! No averaging
971	zwgt1(jic) = 1._wp
972	jpit = jic
973
974	CASE( 1 ) ! Boxcar, width = nn_baro
975	DO jn = 1, 3*nn_baro
976	za1 = ABS(float(jn-jic))/float(nn_baro)
977	IF (za1 < 0.5_wp) THEN
978	zwgt1(jn) = 1._wp
979	jpit = jn
980	ENDIF
981	ENDDO
982
983	CASE( 2 ) ! Boxcar, width = 2 * nn_baro
984	DO jn = 1, 3*nn_baro
985	za1 = ABS(float(jn-jic))/float(nn_baro)
986	IF (za1 < 1._wp) THEN
987	zwgt1(jn) = 1._wp
988	jpit = jn
989	ENDIF
990	ENDDO
991	CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_bt_flt' )
992	END SELECT
993
994	ELSE ! No time averaging
995	zwgt1(jic) = 1._wp
996	jpit = jic
997	ENDIF
998
999	! Set secondary weights
1000	DO jn = 1, jpit
1001	DO ji = jn, jpit
1002	zwgt2(jn) = zwgt2(jn) + zwgt1(ji)
1003	END DO
1004	END DO
1005
1006	! Normalize weigths:
1007	za1 = 1._wp / SUM(zwgt1(1:jpit))
1008	za2 = 1._wp / SUM(zwgt2(1:jpit))
1009	DO jn = 1, jpit
1010	zwgt1(jn) = zwgt1(jn) * za1
1011	zwgt2(jn) = zwgt2(jn) * za2
1012	END DO
1013	!
1014	END SUBROUTINE ts_wgt
1015
1016	SUBROUTINE ts_rst( kt, cdrw )
1017	!!---------------------------------------------------------------------
1018	!! * ROUTINE ts_rst *
1019	!!
1020	!! ** Purpose : Read or write time-splitting arrays in restart file
1021	!!----------------------------------------------------------------------
1022	INTEGER , INTENT(in) :: kt ! ocean time-step
1023	CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag
1024	!
1025	!!----------------------------------------------------------------------
1026	!
1027	IF( TRIM(cdrw) == 'READ' ) THEN
1028	CALL iom_get( numror, jpdom_autoglo, 'ub2_b' , ub2_b (:,:) )
1029	CALL iom_get( numror, jpdom_autoglo, 'vb2_b' , vb2_b (:,:) )
1030	IF( .NOT.ln_bt_av ) THEN
1031	CALL iom_get( numror, jpdom_autoglo, 'sshbb_e' , sshbb_e(:,:) )
1032	CALL iom_get( numror, jpdom_autoglo, 'ubb_e' , ubb_e(:,:) )
1033	CALL iom_get( numror, jpdom_autoglo, 'vbb_e' , vbb_e(:,:) )
1034	CALL iom_get( numror, jpdom_autoglo, 'sshb_e' , sshb_e(:,:) )
1035	CALL iom_get( numror, jpdom_autoglo, 'ub_e' , ub_e(:,:) )
1036	CALL iom_get( numror, jpdom_autoglo, 'vb_e' , vb_e(:,:) )
1037	ENDIF
1038	#if defined key_agrif
1039	! Read time integrated fluxes
1040	IF ( .NOT.Agrif_Root() ) THEN
1041	CALL iom_get( numror, jpdom_autoglo, 'ub2_i_b' , ub2_i_b(:,:) )
1042	CALL iom_get( numror, jpdom_autoglo, 'vb2_i_b' , vb2_i_b(:,:) )
1043	ENDIF
1044	#endif
1045	!
1046	ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN
1047	CALL iom_rstput( kt, nitrst, numrow, 'ub2_b' , ub2_b (:,:) )
1048	CALL iom_rstput( kt, nitrst, numrow, 'vb2_b' , vb2_b (:,:) )
1049	!
1050	IF (.NOT.ln_bt_av) THEN
1051	CALL iom_rstput( kt, nitrst, numrow, 'sshbb_e' , sshbb_e(:,:) )
1052	CALL iom_rstput( kt, nitrst, numrow, 'ubb_e' , ubb_e(:,:) )
1053	CALL iom_rstput( kt, nitrst, numrow, 'vbb_e' , vbb_e(:,:) )
1054	CALL iom_rstput( kt, nitrst, numrow, 'sshb_e' , sshb_e(:,:) )
1055	CALL iom_rstput( kt, nitrst, numrow, 'ub_e' , ub_e(:,:) )
1056	CALL iom_rstput( kt, nitrst, numrow, 'vb_e' , vb_e(:,:) )
1057	ENDIF
1058	#if defined key_agrif
1059	! Save time integrated fluxes
1060	IF ( .NOT.Agrif_Root() ) THEN
1061	CALL iom_rstput( kt, nitrst, numrow, 'ub2_i_b' , ub2_i_b(:,:) )
1062	CALL iom_rstput( kt, nitrst, numrow, 'vb2_i_b' , vb2_i_b(:,:) )
1063	ENDIF
1064	#endif
1065	ENDIF
1066	!
1067	END SUBROUTINE ts_rst
1068
1069	SUBROUTINE dyn_spg_ts_init( kt )
1070	!!---------------------------------------------------------------------
1071	!! * ROUTINE dyn_spg_ts_init *
1072	!!
1073	!! ** Purpose : Set time splitting options
1074	!!----------------------------------------------------------------------
1075	INTEGER , INTENT(in) :: kt ! ocean time-step
1076	!
1077	INTEGER :: ji ,jj
1078	INTEGER :: ios ! Local integer output status for namelist read
1079	REAL(wp) :: zxr2, zyr2, zcmax
1080	REAL(wp), POINTER, DIMENSION(:,:) :: zcu
1081	!!
1082	NAMELIST/namsplit/ ln_bt_fw, ln_bt_av, ln_bt_nn_auto, &
1083	& nn_baro, rn_bt_cmax, nn_bt_flt
1084	!!----------------------------------------------------------------------
1085	!
1086	REWIND( numnam_ref ) ! Namelist namsplit in reference namelist : time splitting parameters
1087	READ ( numnam_ref, namsplit, IOSTAT = ios, ERR = 901)
1088	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsplit in reference namelist', lwp )
1089
1090	REWIND( numnam_cfg ) ! Namelist namsplit in configuration namelist : time splitting parameters
1091	READ ( numnam_cfg, namsplit, IOSTAT = ios, ERR = 902 )
1092	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsplit in configuration namelist', lwp )
1093	IF(lwm) WRITE ( numond, namsplit )
1094	!
1095	! ! Max courant number for ext. grav. waves
1096	!
1097	CALL wrk_alloc( jpi, jpj, zcu )
1098	!
1099	IF (lk_vvl) THEN
1100	DO jj = 1, jpj
1101	DO ji =1, jpi
1102	zxr2 = 1./(e1t(ji,jj)*e1t(ji,jj))
1103	zyr2 = 1./(e2t(ji,jj)*e2t(ji,jj))
1104	zcu(ji,jj) = sqrt(gravht_0(ji,jj)(zxr2 + zyr2) )
1105	END DO
1106	END DO
1107	ELSE
1108	DO jj = 1, jpj
1109	DO ji =1, jpi
1110	zxr2 = 1./(e1t(ji,jj)*e1t(ji,jj))
1111	zyr2 = 1./(e2t(ji,jj)*e2t(ji,jj))
1112	zcu(ji,jj) = sqrt(gravht(ji,jj)(zxr2 + zyr2) )
1113	END DO
1114	END DO
1115	ENDIF
1116
1117	zcmax = MAXVAL(zcu(:,:))
1118	IF( lk_mpp ) CALL mpp_max( zcmax )
1119
1120	! Estimate number of iterations to satisfy a max courant number= rn_bt_cmax
1121	IF (ln_bt_nn_auto) nn_baro = CEILING( rdt / rn_bt_cmax * zcmax)
1122
1123	rdtbt = rdt / FLOAT(nn_baro)
1124	zcmax = zcmax * rdtbt
1125	! Print results
1126	IF(lwp) WRITE(numout,*)
1127	IF(lwp) WRITE(numout,*) 'dyn_spg_ts : split-explicit free surface'
1128	IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
1129	IF( ln_bt_nn_auto ) THEN
1130	IF(lwp) WRITE(numout,*) ' ln_ts_nn_auto=.true. Automatically set nn_baro '
1131	IF(lwp) WRITE(numout,*) ' Max. courant number allowed: ', rn_bt_cmax
1132	ELSE
1133	IF(lwp) WRITE(numout,*) ' ln_ts_nn_auto=.false.: Use nn_baro in namelist '
1134	ENDIF
1135
1136	IF(ln_bt_av) THEN
1137	IF(lwp) WRITE(numout,*) ' ln_bt_av=.true. => Time averaging over nn_baro time steps is on '
1138	ELSE
1139	IF(lwp) WRITE(numout,*) ' ln_bt_av=.false. => No time averaging of barotropic variables '
1140	ENDIF
1141	!
1142	!
1143	IF(ln_bt_fw) THEN
1144	IF(lwp) WRITE(numout,*) ' ln_bt_fw=.true. => Forward integration of barotropic variables '
1145	ELSE
1146	IF(lwp) WRITE(numout,*) ' ln_bt_fw =.false.=> Centred integration of barotropic variables '
1147	ENDIF
1148	!
1149	#if defined key_agrif
1150	! Restrict the use of Agrif to the forward case only
1151	IF ((.NOT.ln_bt_fw ).AND.(.NOT.Agrif_Root())) CALL ctl_stop( 'AGRIF not implemented if ln_bt_fw=.FALSE.' )
1152	#endif
1153	!
1154	IF(lwp) WRITE(numout,*) ' Time filter choice, nn_bt_flt: ', nn_bt_flt
1155	SELECT CASE ( nn_bt_flt )
1156	CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' Dirac'
1157	CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' Boxcar: width = nn_baro'
1158	CASE( 2 ) ; IF(lwp) WRITE(numout,) ' Boxcar: width = 2nn_baro'
1159	CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_bt_flt: should 0,1,2' )
1160	END SELECT
1161	!
1162	IF(lwp) WRITE(numout,*) ' '
1163	IF(lwp) WRITE(numout,*) ' nn_baro = ', nn_baro
1164	IF(lwp) WRITE(numout,*) ' Barotropic time step [s] is :', rdtbt
1165	IF(lwp) WRITE(numout,*) ' Maximum Courant number is :', zcmax
1166	!
1167	IF ((.NOT.ln_bt_av).AND.(.NOT.ln_bt_fw)) THEN
1168	CALL ctl_stop( 'dynspg_ts ERROR: No time averaging => only forward integration is possible' )
1169	ENDIF
1170	IF ( zcmax>0.9_wp ) THEN
1171	CALL ctl_stop( 'dynspg_ts ERROR: Maximum Courant number is greater than 0.9: Inc. nn_baro !' )
1172	ENDIF
1173	!
1174	CALL wrk_dealloc( jpi, jpj, zcu )
1175	!
1176	END SUBROUTINE dyn_spg_ts_init
1177
1178	#else
1179	!!---------------------------------------------------------------------------
1180	!! Default case : Empty module No split explicit free surface
1181	!!---------------------------------------------------------------------------
1182	CONTAINS
1183	INTEGER FUNCTION dyn_spg_ts_alloc() ! Dummy function
1184	dyn_spg_ts_alloc = 0
1185	END FUNCTION dyn_spg_ts_alloc
1186	SUBROUTINE dyn_spg_ts( kt ) ! Empty routine
1187	INTEGER, INTENT(in) :: kt
1188	WRITE(,) 'dyn_spg_ts: You should not have seen this print! error?', kt
1189	END SUBROUTINE dyn_spg_ts
1190	SUBROUTINE ts_rst( kt, cdrw ) ! Empty routine
1191	INTEGER , INTENT(in) :: kt ! ocean time-step
1192	CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag
1193	WRITE(,) 'ts_rst : You should not have seen this print! error?', kt, cdrw
1194	END SUBROUTINE ts_rst
1195	SUBROUTINE dyn_spg_ts_init( kt ) ! Empty routine
1196	INTEGER , INTENT(in) :: kt ! ocean time-step
1197	WRITE(,) 'dyn_spg_ts_init : You should not have seen this print! error?', kt
1198	END SUBROUTINE dyn_spg_ts_init
1199	#endif
1200
1201	!!======================================================================
1202	END MODULE dynspg_ts
1203
1204
1205

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/r6232_collate_bgc_diagnostics/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 11134

Download in other formats: