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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

dynspg_ts.F90 in branches/UKMO/r5936_CO6_CO5_shelfdiagnostic/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

Context Navigation

source: branches/UKMO/r5936_CO6_CO5_shelfdiagnostic/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 7113

Last change on this file since 7113 was 7113, checked in by jcastill, 7 years ago
Remove again the svn keywords, as it did not work before
File size: 53.6 KB

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

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

Download in other formats: