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

Context Navigation

source: branches/UKMO/dev_r5518_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 9188

Last change on this file since 9188 was 9188, checked in by kingr, 7 years ago
Small change to zdftke/zdfgls to avoid merge conflicts with AMM15 branch.
File size: 55.8 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 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	! Arrays below are saved to allow testing of the "no time averaging" option
71	! If this option is not retained, these could be replaced by temporary arrays
72	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: sshbb_e, sshb_e, & ! Instantaneous barotropic arrays
73	ubb_e, ub_e, &
74	vbb_e, vb_e
75
76	!! * Substitutions
77	# include "domzgr_substitute.h90"
78	# include "vectopt_loop_substitute.h90"
79	!!----------------------------------------------------------------------
80	!! NEMO/OPA 3.5 , NEMO Consortium (2013)
81	!! $Id$
82	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
83	!!----------------------------------------------------------------------
84	CONTAINS
85
86	INTEGER FUNCTION dyn_spg_ts_alloc()
87	!!----------------------------------------------------------------------
88	!! * routine dyn_spg_ts_alloc *
89	!!----------------------------------------------------------------------
90	INTEGER :: ierr(3)
91	!!----------------------------------------------------------------------
92	ierr(:) = 0
93
94	ALLOCATE( sshb_e(jpi,jpj), sshbb_e(jpi,jpj), &
95	& ub_e(jpi,jpj) , vb_e(jpi,jpj) , &
96	& ubb_e(jpi,jpj) , vbb_e(jpi,jpj) , STAT= ierr(1) )
97
98	ALLOCATE( wgtbtp1(3nn_baro), wgtbtp2(3nn_baro), zwz(jpi,jpj), STAT= ierr(2) )
99
100	IF( ln_dynvor_een .or. ln_dynvor_een_old ) ALLOCATE( ftnw(jpi,jpj) , ftne(jpi,jpj) , &
101	& ftsw(jpi,jpj) , ftse(jpi,jpj) , STAT=ierr(3) )
102
103	dyn_spg_ts_alloc = MAXVAL(ierr(:))
104
105	IF( lk_mpp ) CALL mpp_sum( dyn_spg_ts_alloc )
106	IF( dyn_spg_ts_alloc /= 0 ) CALL ctl_warn('dynspg_oce_alloc: failed to allocate arrays')
107	!
108	END FUNCTION dyn_spg_ts_alloc
109
110	SUBROUTINE dyn_spg_ts( kt )
111	!!----------------------------------------------------------------------
112	!!
113	!! ** Purpose :
114	!! -Compute the now trend due to the explicit time stepping
115	!! of the quasi-linear barotropic system.
116	!!
117	!! ** Method :
118	!! Barotropic variables are advanced from internal time steps
119	!! "n" to "n+1" if ln_bt_fw=T
120	!! or from
121	!! "n-1" to "n+1" if ln_bt_fw=F
122	!! thanks to a generalized forward-backward time stepping (see ref. below).
123	!!
124	!! ** Action :
125	!! -Update the filtered free surface at step "n+1" : ssha
126	!! -Update filtered barotropic velocities at step "n+1" : ua_b, va_b
127	!! -Compute barotropic advective velocities at step "n" : un_adv, vn_adv
128	!! These are used to advect tracers and are compliant with discrete
129	!! continuity equation taken at the baroclinic time steps. This
130	!! ensures tracers conservation.
131	!! -Update 3d trend (ua, va) with barotropic component.
132	!!
133	!! References : Shchepetkin, A.F. and J.C. McWilliams, 2005:
134	!! The regional oceanic modeling system (ROMS):
135	!! a split-explicit, free-surface,
136	!! topography-following-coordinate oceanic model.
137	!! Ocean Modelling, 9, 347-404.
138	!!---------------------------------------------------------------------
139	!
140	INTEGER, INTENT(in) :: kt ! ocean time-step index
141	!
142	LOGICAL :: ll_fw_start ! if true, forward integration
143	LOGICAL :: ll_init ! if true, special startup of 2d equations
144	INTEGER :: ji, jj, jk, jn ! dummy loop indices
145	INTEGER :: ikbu, ikbv, noffset ! local integers
146	REAL(wp) :: zraur, z1_2dt_b, z2dt_bf ! local scalars
147	REAL(wp) :: zx1, zy1, zx2, zy2 ! - -
148	REAL(wp) :: z1_12, z1_8, z1_4, z1_2 ! - -
149	REAL(wp) :: zu_spg, zv_spg ! - -
150	REAL(wp) :: zhura, zhvra ! - -
151	REAL(wp) :: za0, za1, za2, za3 ! - -
152	!
153	REAL(wp), POINTER, DIMENSION(:,:) :: zun_e, zvn_e, zsshp2_e
154	REAL(wp), POINTER, DIMENSION(:,:) :: zu_trd, zv_trd, zu_frc, zv_frc, zssh_frc
155	REAL(wp), POINTER, DIMENSION(:,:) :: zu_sum, zv_sum, zwx, zwy, zhdiv
156	REAL(wp), POINTER, DIMENSION(:,:) :: zhup2_e, zhvp2_e, zhust_e, zhvst_e
157	REAL(wp), POINTER, DIMENSION(:,:) :: zsshu_a, zsshv_a
158	REAL(wp), POINTER, DIMENSION(:,:) :: zhf
159	!!----------------------------------------------------------------------
160	!
161	IF( nn_timing == 1 ) CALL timing_start('dyn_spg_ts')
162	!
163	! !* Allocate temporary arrays
164	CALL wrk_alloc( jpi, jpj, zsshp2_e, zhdiv )
165	CALL wrk_alloc( jpi, jpj, zu_trd, zv_trd, zun_e, zvn_e )
166	CALL wrk_alloc( jpi, jpj, zwx, zwy, zu_sum, zv_sum, zssh_frc, zu_frc, zv_frc)
167	CALL wrk_alloc( jpi, jpj, zhup2_e, zhvp2_e, zhust_e, zhvst_e)
168	CALL wrk_alloc( jpi, jpj, zsshu_a, zsshv_a )
169	CALL wrk_alloc( jpi, jpj, zhf )
170	!
171	! !* Local constant initialization
172	z1_12 = 1._wp / 12._wp
173	z1_8 = 0.125_wp
174	z1_4 = 0.25_wp
175	z1_2 = 0.5_wp
176	zraur = 1._wp / rau0
177	!
178	IF( kt == nit000 .AND. neuler == 0 ) THEN ! reciprocal of baroclinic time step
179	z2dt_bf = rdt
180	ELSE
181	z2dt_bf = 2.0_wp * rdt
182	ENDIF
183	z1_2dt_b = 1.0_wp / z2dt_bf
184	!
185	ll_init = ln_bt_av ! if no time averaging, then no specific restart
186	ll_fw_start = .FALSE.
187	!
188	! time offset in steps for bdy data update
189	IF (.NOT.ln_bt_fw) THEN ; noffset=-2*nn_baro ; ELSE ; noffset = 0 ; ENDIF
190	!
191	IF( kt == nit000 ) THEN !* initialisation
192	!
193	IF(lwp) WRITE(numout,*)
194	IF(lwp) WRITE(numout,*) 'dyn_spg_ts : surface pressure gradient trend'
195	IF(lwp) WRITE(numout,*) '~~~~~~~~~~ free surface with time splitting'
196	IF(lwp) WRITE(numout,*)
197	!
198	IF (neuler==0) ll_init=.TRUE.
199	!
200	IF (ln_bt_fw.OR.(neuler==0)) THEN
201	ll_fw_start=.TRUE.
202	noffset = 0
203	ELSE
204	ll_fw_start=.FALSE.
205	ENDIF
206	!
207	! Set averaging weights and cycle length:
208	CALL ts_wgt(ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2)
209	!
210	!
211	ENDIF
212	!
213	! Set arrays to remove/compute coriolis trend.
214	! Do it once at kt=nit000 if volume is fixed, else at each long time step.
215	! Note that these arrays are also used during barotropic loop. These are however frozen
216	! although they should be updated in the variable volume case. Not a big approximation.
217	! To remove this approximation, copy lines below inside barotropic loop
218	! and update depths at T-F points (ht and zhf resp.) at each barotropic time step
219	!
220	IF ( kt == nit000 .OR. lk_vvl ) THEN
221	IF ( ln_dynvor_een_old ) THEN
222	DO jj = 1, jpjm1
223	DO ji = 1, jpim1
224	zwz(ji,jj) = ( ht(ji ,jj+1) + ht(ji+1,jj+1) + &
225	& ht(ji ,jj ) + ht(ji+1,jj ) ) / 4._wp
226	IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = 1._wp / zwz(ji,jj)
227	END DO
228	END DO
229	CALL lbc_lnk( zwz, 'F', 1._wp )
230	zwz(:,:) = ff(:,:) * zwz(:,:)
231
232	ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp
233	DO jj = 2, jpj
234	DO ji = fs_2, jpi ! vector opt.
235	ftne(ji,jj) = zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1)
236	ftnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj ) + zwz(ji ,jj )
237	ftse(ji,jj) = zwz(ji ,jj ) + zwz(ji ,jj-1) + zwz(ji-1,jj-1)
238	ftsw(ji,jj) = zwz(ji ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj )
239	END DO
240	END DO
241	ELSE IF ( ln_dynvor_een ) THEN
242	DO jj = 1, jpjm1
243	DO ji = 1, jpim1
244	zwz(ji,jj) = ( ht(ji ,jj+1) + ht(ji+1,jj+1) + &
245	& ht(ji ,jj ) + ht(ji+1,jj ) ) &
246	& / ( MAX( 1.0_wp, tmask(ji ,jj+1, 1) + tmask(ji+1,jj+1, 1) + &
247	& tmask(ji ,jj , 1) + tmask(ji+1,jj , 1) ) )
248	IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = 1._wp / zwz(ji,jj)
249	END DO
250	END DO
251	CALL lbc_lnk( zwz, 'F', 1._wp )
252	zwz(:,:) = ff(:,:) * zwz(:,:)
253
254	ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp
255	DO jj = 2, jpj
256	DO ji = fs_2, jpi ! vector opt.
257	ftne(ji,jj) = zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1)
258	ftnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj ) + zwz(ji ,jj )
259	ftse(ji,jj) = zwz(ji ,jj ) + zwz(ji ,jj-1) + zwz(ji-1,jj-1)
260	ftsw(ji,jj) = zwz(ji ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj )
261	END DO
262	END DO
263	ELSE
264	zwz(:,:) = 0._wp
265	zhf(:,:) = 0.
266	IF ( .not. ln_sco ) THEN
267	! JC: It not clear yet what should be the depth at f-points over land in z-coordinate case
268	! Set it to zero for the time being
269	! IF( rn_hmin < 0._wp ) THEN ; jk = - INT( rn_hmin ) ! from a nb of level
270	! ELSE ; jk = MINLOC( gdepw_0, mask = gdepw_0 > rn_hmin, dim = 1 ) ! from a depth
271	! ENDIF
272	! zhf(:,:) = gdepw_0(:,:,jk+1)
273	ELSE
274	zhf(:,:) = hbatf(:,:)
275	END IF
276
277	DO jj = 1, jpjm1
278	zhf(:,jj) = zhf(:,jj)(1._wp- umask(:,jj,1) umask(:,jj+1,1))
279	END DO
280
281	DO jk = 1, jpkm1
282	DO jj = 1, jpjm1
283	zhf(:,jj) = zhf(:,jj) + fse3f_n(:,jj,jk) * umask(:,jj,jk) * umask(:,jj+1,jk)
284	END DO
285	END DO
286	CALL lbc_lnk( zhf, 'F', 1._wp )
287	! JC: TBC. hf should be greater than 0
288	DO jj = 1, jpj
289	DO ji = 1, jpi
290	IF( zhf(ji,jj) /= 0._wp ) zwz(ji,jj) = 1._wp / zhf(ji,jj) ! zhf is actually hf here but it saves an array
291	END DO
292	END DO
293	zwz(:,:) = ff(:,:) * zwz(:,:)
294	ENDIF
295	ENDIF
296	!
297	! If forward start at previous time step, and centered integration,
298	! then update averaging weights:
299	IF ((.NOT.ln_bt_fw).AND.((neuler==0).AND.(kt==nit000+1))) THEN
300	ll_fw_start=.FALSE.
301	CALL ts_wgt(ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2)
302	ENDIF
303
304	! -----------------------------------------------------------------------------
305	! Phase 1 : Coupling between general trend and barotropic estimates (1st step)
306	! -----------------------------------------------------------------------------
307	!
308	!
309	! !* e3*d/dt(Ua) (Vertically integrated)
310	! ! --------------------------------------------------
311	zu_frc(:,:) = 0._wp
312	zv_frc(:,:) = 0._wp
313	!
314	DO jk = 1, jpkm1
315	zu_frc(:,:) = zu_frc(:,:) + fse3u_n(:,:,jk) * ua(:,:,jk) * umask(:,:,jk)
316	zv_frc(:,:) = zv_frc(:,:) + fse3v_n(:,:,jk) * va(:,:,jk) * vmask(:,:,jk)
317	END DO
318	!
319	zu_frc(:,:) = zu_frc(:,:) * hur(:,:)
320	zv_frc(:,:) = zv_frc(:,:) * hvr(:,:)
321	!
322	!
323	! !* baroclinic momentum trend (remove the vertical mean trend)
324	DO jk = 1, jpkm1 ! -----------------------------------------------------------
325	DO jj = 2, jpjm1
326	DO ji = fs_2, fs_jpim1 ! vector opt.
327	ua(ji,jj,jk) = ua(ji,jj,jk) - zu_frc(ji,jj) * umask(ji,jj,jk)
328	va(ji,jj,jk) = va(ji,jj,jk) - zv_frc(ji,jj) * vmask(ji,jj,jk)
329	END DO
330	END DO
331	END DO
332	! !* barotropic Coriolis trends (vorticity scheme dependent)
333	! ! --------------------------------------------------------
334	zwx(:,:) = un_b(:,:) * hu(:,:) * e2u(:,:) ! now fluxes
335	zwy(:,:) = vn_b(:,:) * hv(:,:) * e1v(:,:)
336	!
337	IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN ! energy conserving or mixed scheme
338	DO jj = 2, jpjm1
339	DO ji = fs_2, fs_jpim1 ! vector opt.
340	zy1 = ( zwy(ji,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj)
341	zy2 = ( zwy(ji,jj ) + zwy(ji+1,jj ) ) / e1u(ji,jj)
342	zx1 = ( zwx(ji-1,jj) + zwx(ji-1,jj+1) ) / e2v(ji,jj)
343	zx2 = ( zwx(ji ,jj) + zwx(ji ,jj+1) ) / e2v(ji,jj)
344	! energy conserving formulation for planetary vorticity term
345	zu_trd(ji,jj) = z1_4 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
346	zv_trd(ji,jj) =-z1_4 * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 )
347	END DO
348	END DO
349	!
350	ELSEIF ( ln_dynvor_ens ) THEN ! enstrophy conserving scheme
351	DO jj = 2, jpjm1
352	DO ji = fs_2, fs_jpim1 ! vector opt.
353	zy1 = z1_8 * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) &
354	& + zwy(ji ,jj ) + zwy(ji+1,jj ) ) / e1u(ji,jj)
355	zx1 = - z1_8 * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) &
356	& + zwx(ji ,jj ) + zwx(ji ,jj+1) ) / e2v(ji,jj)
357	zu_trd(ji,jj) = zy1 * ( zwz(ji ,jj-1) + zwz(ji,jj) )
358	zv_trd(ji,jj) = zx1 * ( zwz(ji-1,jj ) + zwz(ji,jj) )
359	END DO
360	END DO
361	!
362	ELSEIF ( ln_dynvor_een .or. ln_dynvor_een_old ) THEN ! enstrophy and energy conserving scheme
363	DO jj = 2, jpjm1
364	DO ji = fs_2, fs_jpim1 ! vector opt.
365	zu_trd(ji,jj) = + z1_12 / e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) &
366	& + ftnw(ji+1,jj) * zwy(ji+1,jj ) &
367	& + ftse(ji,jj ) * zwy(ji ,jj-1) &
368	& + ftsw(ji+1,jj) * zwy(ji+1,jj-1) )
369	zv_trd(ji,jj) = - z1_12 / e2v(ji,jj) * ( ftsw(ji,jj+1) * zwx(ji-1,jj+1) &
370	& + ftse(ji,jj+1) * zwx(ji ,jj+1) &
371	& + ftnw(ji,jj ) * zwx(ji-1,jj ) &
372	& + ftne(ji,jj ) * zwx(ji ,jj ) )
373	END DO
374	END DO
375	!
376	ENDIF
377	!
378	! !* Right-Hand-Side of the barotropic momentum equation
379	! ! ----------------------------------------------------
380	IF( lk_vvl ) THEN ! Variable volume : remove surface pressure gradient
381	DO jj = 2, jpjm1
382	DO ji = fs_2, fs_jpim1 ! vector opt.
383	zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) / e1u(ji,jj)
384	zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji ,jj+1) - sshn(ji ,jj ) ) / e2v(ji,jj)
385	END DO
386	END DO
387	ENDIF
388
389	DO jj = 2, jpjm1 ! Remove coriolis term (and possibly spg) from barotropic trend
390	DO ji = fs_2, fs_jpim1
391	zu_frc(ji,jj) = zu_frc(ji,jj) - zu_trd(ji,jj) * umask(ji,jj,1)
392	zv_frc(ji,jj) = zv_frc(ji,jj) - zv_trd(ji,jj) * vmask(ji,jj,1)
393	END DO
394	END DO
395	!
396	! ! Add bottom stress contribution from baroclinic velocities:
397	IF (ln_bt_fw) THEN
398	DO jj = 2, jpjm1
399	DO ji = fs_2, fs_jpim1 ! vector opt.
400	ikbu = mbku(ji,jj)
401	ikbv = mbkv(ji,jj)
402	zwx(ji,jj) = un(ji,jj,ikbu) - un_b(ji,jj) ! NOW bottom baroclinic velocities
403	zwy(ji,jj) = vn(ji,jj,ikbv) - vn_b(ji,jj)
404	END DO
405	END DO
406	ELSE
407	DO jj = 2, jpjm1
408	DO ji = fs_2, fs_jpim1 ! vector opt.
409	ikbu = mbku(ji,jj)
410	ikbv = mbkv(ji,jj)
411	zwx(ji,jj) = ub(ji,jj,ikbu) - ub_b(ji,jj) ! BEFORE bottom baroclinic velocities
412	zwy(ji,jj) = vb(ji,jj,ikbv) - vb_b(ji,jj)
413	END DO
414	END DO
415	ENDIF
416	!
417	! Note that the "unclipped" bottom friction parameter is used even with explicit drag
418	zu_frc(:,:) = zu_frc(:,:) + hur(:,:) * bfrua(:,:) * zwx(:,:)
419	zv_frc(:,:) = zv_frc(:,:) + hvr(:,:) * bfrva(:,:) * zwy(:,:)
420	!
421	IF (ln_bt_fw) THEN ! Add wind forcing
422	zu_frc(:,:) = zu_frc(:,:) + zraur * utau(:,:) * hur(:,:)
423	zv_frc(:,:) = zv_frc(:,:) + zraur * vtau(:,:) * hvr(:,:)
424	ELSE
425	zu_frc(:,:) = zu_frc(:,:) + zraur * z1_2 * ( utau_b(:,:) + utau(:,:) ) * hur(:,:)
426	zv_frc(:,:) = zv_frc(:,:) + zraur * z1_2 * ( vtau_b(:,:) + vtau(:,:) ) * hvr(:,:)
427	ENDIF
428	!
429	IF ( ln_apr_dyn ) THEN ! Add atm pressure forcing
430	IF (ln_bt_fw) THEN
431	DO jj = 2, jpjm1
432	DO ji = fs_2, fs_jpim1 ! vector opt.
433	zu_spg = grav * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) ) /e1u(ji,jj)
434	zv_spg = grav * ( ssh_ib (ji ,jj+1) - ssh_ib (ji,jj) ) /e2v(ji,jj)
435	zu_frc(ji,jj) = zu_frc(ji,jj) + zu_spg
436	zv_frc(ji,jj) = zv_frc(ji,jj) + zv_spg
437	END DO
438	END DO
439	ELSE
440	DO jj = 2, jpjm1
441	DO ji = fs_2, fs_jpim1 ! vector opt.
442	zu_spg = grav * z1_2 * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) &
443	& + ssh_ibb(ji+1,jj ) - ssh_ibb(ji,jj) ) /e1u(ji,jj)
444	zv_spg = grav * z1_2 * ( ssh_ib (ji ,jj+1) - ssh_ib (ji,jj) &
445	& + ssh_ibb(ji ,jj+1) - ssh_ibb(ji,jj) ) /e2v(ji,jj)
446	zu_frc(ji,jj) = zu_frc(ji,jj) + zu_spg
447	zv_frc(ji,jj) = zv_frc(ji,jj) + zv_spg
448	END DO
449	END DO
450	ENDIF
451	ENDIF
452	! !* Right-Hand-Side of the barotropic ssh equation
453	! ! -----------------------------------------------
454	! ! Surface net water flux and rivers
455	IF (ln_bt_fw) THEN
456	zssh_frc(:,:) = zraur * ( emp(:,:) - rnf(:,:) + rdivisf * fwfisf(:,:) )
457	ELSE
458	zssh_frc(:,:) = zraur * z1_2 * ( emp(:,:) + emp_b(:,:) - rnf(:,:) - rnf_b(:,:) &
459	& + rdivisf * ( fwfisf(:,:) + fwfisf_b(:,:) ) )
460	ENDIF
461	#if defined key_asminc
462	! ! Include the IAU weighted SSH increment
463	IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN
464	zssh_frc(:,:) = zssh_frc(:,:) - ssh_iau(:,:)
465	ENDIF
466	#endif
467	! !* Fill boundary data arrays for AGRIF
468	! ! ------------------------------------
469	#if defined key_agrif
470	IF( .NOT.Agrif_Root() ) CALL agrif_dta_ts( kt )
471	#endif
472	!
473	! -----------------------------------------------------------------------
474	! Phase 2 : Integration of the barotropic equations
475	! -----------------------------------------------------------------------
476	!
477	! ! ==================== !
478	! ! Initialisations !
479	! ! ==================== !
480	! Initialize barotropic variables:
481	IF( ll_init )THEN
482	sshbb_e(:,:) = 0._wp
483	ubb_e (:,:) = 0._wp
484	vbb_e (:,:) = 0._wp
485	sshb_e (:,:) = 0._wp
486	ub_e (:,:) = 0._wp
487	vb_e (:,:) = 0._wp
488	ENDIF
489	!
490	IF (ln_bt_fw) THEN ! FORWARD integration: start from NOW fields
491	sshn_e(:,:) = sshn (:,:)
492	zun_e (:,:) = un_b (:,:)
493	zvn_e (:,:) = vn_b (:,:)
494	!
495	hu_e (:,:) = hu (:,:)
496	hv_e (:,:) = hv (:,:)
497	hur_e (:,:) = hur (:,:)
498	hvr_e (:,:) = hvr (:,:)
499	ELSE ! CENTRED integration: start from BEFORE fields
500	sshn_e(:,:) = sshb (:,:)
501	zun_e (:,:) = ub_b (:,:)
502	zvn_e (:,:) = vb_b (:,:)
503	!
504	hu_e (:,:) = hu_b (:,:)
505	hv_e (:,:) = hv_b (:,:)
506	hur_e (:,:) = hur_b(:,:)
507	hvr_e (:,:) = hvr_b(:,:)
508	ENDIF
509	!
510	!
511	!
512	! Initialize sums:
513	ua_b (:,:) = 0._wp ! After barotropic velocities (or transport if flux form)
514	va_b (:,:) = 0._wp
515	ssha (:,:) = 0._wp ! Sum for after averaged sea level
516	zu_sum(:,:) = 0._wp ! Sum for now transport issued from ts loop
517	zv_sum(:,:) = 0._wp
518	! ! ==================== !
519	DO jn = 1, icycle ! sub-time-step loop !
520	! ! ==================== !
521	! !* Update the forcing (BDY and tides)
522	! ! ------------------
523	! Update only tidal forcing at open boundaries
524	#if defined key_tide
525	IF ( lk_bdy .AND. lk_tide ) CALL bdy_dta_tides( kt, kit=jn, time_offset=(noffset+1) )
526	IF ( ln_tide_pot .AND. lk_tide ) CALL upd_tide( kt, kit=jn, koffset=noffset )
527	#endif
528	!
529	! Set extrapolation coefficients for predictor step:
530	IF ((jn<3).AND.ll_init) THEN ! Forward
531	za1 = 1._wp
532	za2 = 0._wp
533	za3 = 0._wp
534	ELSE ! AB3-AM4 Coefficients: bet=0.281105
535	za1 = 1.781105_wp ! za1 = 3/2 + bet
536	za2 = -1.06221_wp ! za2 = -(1/2 + 2*bet)
537	za3 = 0.281105_wp ! za3 = bet
538	ENDIF
539
540	! Extrapolate barotropic velocities at step jit+0.5:
541	ua_e(:,:) = za1 * zun_e(:,:) + za2 * ub_e(:,:) + za3 * ubb_e(:,:)
542	va_e(:,:) = za1 * zvn_e(:,:) + za2 * vb_e(:,:) + za3 * vbb_e(:,:)
543
544	IF( lk_vvl ) THEN !* Update ocean depth (variable volume case only)
545	! ! ------------------
546	! Extrapolate Sea Level at step jit+0.5:
547	zsshp2_e(:,:) = za1 * sshn_e(:,:) + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:)
548	!
549	DO jj = 2, jpjm1 ! Sea Surface Height at u- & v-points
550	DO ji = 2, fs_jpim1 ! Vector opt.
551	zwx(ji,jj) = z1_2 * umask(ji,jj,1) * r1_e12u(ji,jj) &
552	& * ( e12t(ji ,jj) * zsshp2_e(ji ,jj) &
553	& + e12t(ji+1,jj) * zsshp2_e(ji+1,jj) )
554	zwy(ji,jj) = z1_2 * vmask(ji,jj,1) * r1_e12v(ji,jj) &
555	& * ( e12t(ji,jj ) * zsshp2_e(ji,jj ) &
556	& + e12t(ji,jj+1) * zsshp2_e(ji,jj+1) )
557	END DO
558	END DO
559	CALL lbc_lnk_multi( zwx, 'U', 1._wp, zwy, 'V', 1._wp )
560	!
561	zhup2_e (:,:) = hu_0(:,:) + zwx(:,:) ! Ocean depth at U- and V-points
562	zhvp2_e (:,:) = hv_0(:,:) + zwy(:,:)
563	ELSE
564	zhup2_e (:,:) = hu(:,:)
565	zhvp2_e (:,:) = hv(:,:)
566	ENDIF
567	! !* after ssh
568	! ! -----------
569	! One should enforce volume conservation at open boundaries here
570	! considering fluxes below:
571	!
572	zwx(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:) ! fluxes at jn+0.5
573	zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:)
574	!
575	#if defined key_agrif
576	! Set fluxes during predictor step to ensure
577	! volume conservation
578	IF( (.NOT.Agrif_Root()).AND.ln_bt_fw ) THEN
579	IF((nbondi == -1).OR.(nbondi == 2)) THEN
580	DO jj=1,jpj
581	zwx(2,jj) = ubdy_w(jj) * e2u(2,jj)
582	END DO
583	ENDIF
584	IF((nbondi == 1).OR.(nbondi == 2)) THEN
585	DO jj=1,jpj
586	zwx(nlci-2,jj) = ubdy_e(jj) * e2u(nlci-2,jj)
587	END DO
588	ENDIF
589	IF((nbondj == -1).OR.(nbondj == 2)) THEN
590	DO ji=1,jpi
591	zwy(ji,2) = vbdy_s(ji) * e1v(ji,2)
592	END DO
593	ENDIF
594	IF((nbondj == 1).OR.(nbondj == 2)) THEN
595	DO ji=1,jpi
596	zwy(ji,nlcj-2) = vbdy_n(ji) * e1v(ji,nlcj-2)
597	END DO
598	ENDIF
599	ENDIF
600	#endif
601	!
602	! Sum over sub-time-steps to compute advective velocities
603	za2 = wgtbtp2(jn)
604	zu_sum (:,:) = zu_sum (:,:) + za2 * zwx (:,:) / e2u (:,:)
605	zv_sum (:,:) = zv_sum (:,:) + za2 * zwy (:,:) / e1v (:,:)
606	!
607	! Set next sea level:
608	DO jj = 2, jpjm1
609	DO ji = fs_2, fs_jpim1 ! vector opt.
610	zhdiv(ji,jj) = ( zwx(ji,jj) - zwx(ji-1,jj) &
611	& + zwy(ji,jj) - zwy(ji,jj-1) ) * r1_e12t(ji,jj)
612	END DO
613	END DO
614	ssha_e(:,:) = ( sshn_e(:,:) - rdtbt * ( zssh_frc(:,:) + zhdiv(:,:) ) ) * tmask(:,:,1)
615	CALL lbc_lnk( ssha_e, 'T', 1._wp )
616
617	#if defined key_bdy
618	! Duplicate sea level across open boundaries (this is only cosmetic if lk_vvl=.false.)
619	IF (lk_bdy) CALL bdy_ssh( ssha_e )
620	#endif
621	#if defined key_agrif
622	IF( .NOT.Agrif_Root() ) CALL agrif_ssh_ts( jn )
623	#endif
624	!
625	! Sea Surface Height at u-,v-points (vvl case only)
626	IF ( lk_vvl ) THEN
627	DO jj = 2, jpjm1
628	DO ji = 2, jpim1 ! NO Vector Opt.
629	zsshu_a(ji,jj) = z1_2 * umask(ji,jj,1) * r1_e12u(ji,jj) &
630	& * ( e12t(ji ,jj ) * ssha_e(ji ,jj ) &
631	& + e12t(ji+1,jj ) * ssha_e(ji+1,jj ) )
632	zsshv_a(ji,jj) = z1_2 * vmask(ji,jj,1) * r1_e12v(ji,jj) &
633	& * ( e12t(ji ,jj ) * ssha_e(ji ,jj ) &
634	& + e12t(ji ,jj+1) * ssha_e(ji ,jj+1) )
635	END DO
636	END DO
637	CALL lbc_lnk_multi( zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp )
638	ENDIF
639	!
640	! Half-step back interpolation of SSH for surface pressure computation:
641	!----------------------------------------------------------------------
642	IF ((jn==1).AND.ll_init) THEN
643	za0=1._wp ! Forward-backward
644	za1=0._wp
645	za2=0._wp
646	za3=0._wp
647	ELSEIF ((jn==2).AND.ll_init) THEN ! AB2-AM3 Coefficients; bet=0 ; gam=-1/6 ; eps=1/12
648	za0= 1.0833333333333_wp ! za0 = 1-gam-eps
649	za1=-0.1666666666666_wp ! za1 = gam
650	za2= 0.0833333333333_wp ! za2 = eps
651	za3= 0._wp
652	ELSE ! AB3-AM4 Coefficients; bet=0.281105 ; eps=0.013 ; gam=0.0880
653	za0=0.614_wp ! za0 = 1/2 + gam + 2*eps
654	za1=0.285_wp ! za1 = 1/2 - 2gam - 3eps
655	za2=0.088_wp ! za2 = gam
656	za3=0.013_wp ! za3 = eps
657	ENDIF
658
659	zsshp2_e(:,:) = za0 * ssha_e(:,:) + za1 * sshn_e (:,:) &
660	& + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:)
661
662	!
663	! Compute associated depths at U and V points:
664	IF ( lk_vvl.AND.(.NOT.ln_dynadv_vec) ) THEN
665	!
666	DO jj = 2, jpjm1
667	DO ji = 2, jpim1
668	zx1 = z1_2 * umask(ji ,jj,1) * r1_e12u(ji ,jj) &
669	& * ( e12t(ji ,jj ) * zsshp2_e(ji ,jj) &
670	& + e12t(ji+1,jj ) * zsshp2_e(ji+1,jj ) )
671	zy1 = z1_2 * vmask(ji ,jj,1) * r1_e12v(ji ,jj ) &
672	& * ( e12t(ji ,jj ) * zsshp2_e(ji ,jj ) &
673	& + e12t(ji ,jj+1) * zsshp2_e(ji ,jj+1) )
674	zhust_e(ji,jj) = hu_0(ji,jj) + zx1
675	zhvst_e(ji,jj) = hv_0(ji,jj) + zy1
676	END DO
677	END DO
678	ENDIF
679	!
680	! Add Coriolis trend:
681	! zwz array below or triads normally depend on sea level with key_vvl and should be updated
682	! at each time step. We however keep them constant here for optimization.
683	! Recall that zwx and zwy arrays hold fluxes at this stage:
684	! zwx(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:) ! fluxes at jn+0.5
685	! zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:)
686	!
687	IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN !== energy conserving or mixed scheme ==!
688	DO jj = 2, jpjm1
689	DO ji = fs_2, fs_jpim1 ! vector opt.
690	zy1 = ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj)
691	zy2 = ( zwy(ji ,jj ) + zwy(ji+1,jj ) ) / e1u(ji,jj)
692	zx1 = ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) ) / e2v(ji,jj)
693	zx2 = ( zwx(ji ,jj ) + zwx(ji ,jj+1) ) / e2v(ji,jj)
694	zu_trd(ji,jj) = z1_4 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
695	zv_trd(ji,jj) =-z1_4 * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 )
696	END DO
697	END DO
698	!
699	ELSEIF ( ln_dynvor_ens ) THEN !== enstrophy conserving scheme ==!
700	DO jj = 2, jpjm1
701	DO ji = fs_2, fs_jpim1 ! vector opt.
702	zy1 = z1_8 * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) &
703	& + zwy(ji ,jj ) + zwy(ji+1,jj ) ) / e1u(ji,jj)
704	zx1 = - z1_8 * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) &
705	& + zwx(ji ,jj ) + zwx(ji ,jj+1) ) / e2v(ji,jj)
706	zu_trd(ji,jj) = zy1 * ( zwz(ji ,jj-1) + zwz(ji,jj) )
707	zv_trd(ji,jj) = zx1 * ( zwz(ji-1,jj ) + zwz(ji,jj) )
708	END DO
709	END DO
710	!
711	ELSEIF ( ln_dynvor_een .or. ln_dynvor_een_old ) THEN !== energy and enstrophy conserving scheme ==!
712	DO jj = 2, jpjm1
713	DO ji = fs_2, fs_jpim1 ! vector opt.
714	zu_trd(ji,jj) = + z1_12 / e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) &
715	& + ftnw(ji+1,jj) * zwy(ji+1,jj ) &
716	& + ftse(ji,jj ) * zwy(ji ,jj-1) &
717	& + ftsw(ji+1,jj) * zwy(ji+1,jj-1) )
718	zv_trd(ji,jj) = - z1_12 / e2v(ji,jj) * ( ftsw(ji,jj+1) * zwx(ji-1,jj+1) &
719	& + ftse(ji,jj+1) * zwx(ji ,jj+1) &
720	& + ftnw(ji,jj ) * zwx(ji-1,jj ) &
721	& + ftne(ji,jj ) * zwx(ji ,jj ) )
722	END DO
723	END DO
724	!
725	ENDIF
726	!
727	! Add tidal astronomical forcing if defined
728	IF ( lk_tide.AND.ln_tide_pot ) THEN
729	DO jj = 2, jpjm1
730	DO ji = fs_2, fs_jpim1 ! vector opt.
731	zu_spg = grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) / e1u(ji,jj)
732	zv_spg = grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) / e2v(ji,jj)
733	zu_trd(ji,jj) = zu_trd(ji,jj) + zu_spg
734	zv_trd(ji,jj) = zv_trd(ji,jj) + zv_spg
735	END DO
736	END DO
737	ENDIF
738	!
739	! Add bottom stresses:
740	zu_trd(:,:) = zu_trd(:,:) + bfrua(:,:) * zun_e(:,:) * hur_e(:,:)
741	zv_trd(:,:) = zv_trd(:,:) + bfrva(:,:) * zvn_e(:,:) * hvr_e(:,:)
742	!
743	! Surface pressure trend:
744	DO jj = 2, jpjm1
745	DO ji = fs_2, fs_jpim1 ! vector opt.
746	! Add surface pressure gradient
747	zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) / e1u(ji,jj)
748	zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) / e2v(ji,jj)
749	zwx(ji,jj) = zu_spg
750	zwy(ji,jj) = zv_spg
751	END DO
752	END DO
753	!
754	! Set next velocities:
755	IF( ln_dynadv_vec .OR. (.NOT. lk_vvl) ) THEN ! Vector form
756	DO jj = 2, jpjm1
757	DO ji = fs_2, fs_jpim1 ! vector opt.
758	ua_e(ji,jj) = ( zun_e(ji,jj) &
759	& + rdtbt * ( zwx(ji,jj) &
760	& + zu_trd(ji,jj) &
761	& + zu_frc(ji,jj) ) &
762	& ) * umask(ji,jj,1)
763
764	va_e(ji,jj) = ( zvn_e(ji,jj) &
765	& + rdtbt * ( zwy(ji,jj) &
766	& + zv_trd(ji,jj) &
767	& + zv_frc(ji,jj) ) &
768	& ) * vmask(ji,jj,1)
769	END DO
770	END DO
771
772	ELSE ! Flux form
773	DO jj = 2, jpjm1
774	DO ji = fs_2, fs_jpim1 ! vector opt.
775
776	zhura = umask(ji,jj,1)/(hu_0(ji,jj) + zsshu_a(ji,jj) + 1._wp - umask(ji,jj,1))
777	zhvra = vmask(ji,jj,1)/(hv_0(ji,jj) + zsshv_a(ji,jj) + 1._wp - vmask(ji,jj,1))
778
779	ua_e(ji,jj) = ( hu_e(ji,jj) * zun_e(ji,jj) &
780	& + rdtbt * ( zhust_e(ji,jj) * zwx(ji,jj) &
781	& + zhup2_e(ji,jj) * zu_trd(ji,jj) &
782	& + hu(ji,jj) * zu_frc(ji,jj) ) &
783	& ) * zhura
784
785	va_e(ji,jj) = ( hv_e(ji,jj) * zvn_e(ji,jj) &
786	& + rdtbt * ( zhvst_e(ji,jj) * zwy(ji,jj) &
787	& + zhvp2_e(ji,jj) * zv_trd(ji,jj) &
788	& + hv(ji,jj) * zv_frc(ji,jj) ) &
789	& ) * zhvra
790	END DO
791	END DO
792	ENDIF
793	!
794	IF( lk_vvl ) THEN !* Update ocean depth (variable volume case only)
795	! ! ----------------------------------------------
796	hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:)
797	hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:)
798	hur_e(:,:) = umask(:,:,1) / ( hu_e(:,:) + 1._wp - umask(:,:,1) )
799	hvr_e(:,:) = vmask(:,:,1) / ( hv_e(:,:) + 1._wp - vmask(:,:,1) )
800	!
801	ENDIF
802	! !* domain lateral boundary
803	! ! -----------------------
804	!
805	CALL lbc_lnk_multi( ua_e, 'U', -1._wp, va_e , 'V', -1._wp )
806
807	#if defined key_bdy
808	! open boundaries
809	IF( lk_bdy ) CALL bdy_dyn2d( jn, ua_e, va_e, zun_e, zvn_e, hur_e, hvr_e, ssha_e )
810	#endif
811	#if defined key_agrif
812	IF( .NOT.Agrif_Root() ) CALL agrif_dyn_ts( jn ) ! Agrif
813	#endif
814	! !* Swap
815	! ! ----
816	ubb_e (:,:) = ub_e (:,:)
817	ub_e (:,:) = zun_e (:,:)
818	zun_e (:,:) = ua_e (:,:)
819	!
820	vbb_e (:,:) = vb_e (:,:)
821	vb_e (:,:) = zvn_e (:,:)
822	zvn_e (:,:) = va_e (:,:)
823	!
824	sshbb_e(:,:) = sshb_e(:,:)
825	sshb_e (:,:) = sshn_e(:,:)
826	sshn_e (:,:) = ssha_e(:,:)
827
828	! !* Sum over whole bt loop
829	! ! ----------------------
830	za1 = wgtbtp1(jn)
831	IF (( ln_dynadv_vec ).OR. (.NOT. lk_vvl)) THEN ! Sum velocities
832	ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:)
833	va_b (:,:) = va_b (:,:) + za1 * va_e (:,:)
834	ELSE ! Sum transports
835	ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:) * hu_e (:,:)
836	va_b (:,:) = va_b (:,:) + za1 * va_e (:,:) * hv_e (:,:)
837	ENDIF
838	! ! Sum sea level
839	ssha(:,:) = ssha(:,:) + za1 * ssha_e(:,:)
840	! ! ==================== !
841	END DO ! end loop !
842	! ! ==================== !
843	! -----------------------------------------------------------------------------
844	! Phase 3. update the general trend with the barotropic trend
845	! -----------------------------------------------------------------------------
846	!
847	! At this stage ssha holds a time averaged value
848	! ! Sea Surface Height at u-,v- and f-points
849	IF( lk_vvl ) THEN ! (required only in key_vvl case)
850	DO jj = 1, jpjm1
851	DO ji = 1, jpim1 ! NO Vector Opt.
852	zsshu_a(ji,jj) = z1_2 * umask(ji,jj,1) * r1_e12u(ji,jj) &
853	& * ( e12t(ji ,jj) * ssha(ji ,jj) &
854	& + e12t(ji+1,jj) * ssha(ji+1,jj) )
855	zsshv_a(ji,jj) = z1_2 * vmask(ji,jj,1) * r1_e12v(ji,jj) &
856	& * ( e12t(ji,jj ) * ssha(ji,jj ) &
857	& + e12t(ji,jj+1) * ssha(ji,jj+1) )
858	END DO
859	END DO
860	CALL lbc_lnk_multi( zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) ! Boundary conditions
861	ENDIF
862	!
863	! Set advection velocity correction:
864	IF (((kt==nit000).AND.(neuler==0)).OR.(.NOT.ln_bt_fw)) THEN
865	un_adv(:,:) = zu_sum(:,:)*hur(:,:)
866	vn_adv(:,:) = zv_sum(:,:)*hvr(:,:)
867	ELSE
868	un_adv(:,:) = z1_2 * ( ub2_b(:,:) + zu_sum(:,:)) * hur(:,:)
869	vn_adv(:,:) = z1_2 * ( vb2_b(:,:) + zv_sum(:,:)) * hvr(:,:)
870	END IF
871
872	IF (ln_bt_fw) THEN ! Save integrated transport for next computation
873	ub2_b(:,:) = zu_sum(:,:)
874	vb2_b(:,:) = zv_sum(:,:)
875	ENDIF
876	!
877	! Update barotropic trend:
878	IF (( ln_dynadv_vec ).OR. (.NOT. lk_vvl)) THEN
879	DO jk=1,jpkm1
880	ua(:,:,jk) = ua(:,:,jk) + ( ua_b(:,:) - ub_b(:,:) ) * z1_2dt_b
881	va(:,:,jk) = va(:,:,jk) + ( va_b(:,:) - vb_b(:,:) ) * z1_2dt_b
882	END DO
883	ELSE
884	DO jk=1,jpkm1
885	ua(:,:,jk) = ua(:,:,jk) + hur(:,:) * ( ua_b(:,:) - ub_b(:,:) * hu_b(:,:) ) * z1_2dt_b
886	va(:,:,jk) = va(:,:,jk) + hvr(:,:) * ( va_b(:,:) - vb_b(:,:) * hv_b(:,:) ) * z1_2dt_b
887	END DO
888	! Save barotropic velocities not transport:
889	ua_b (:,:) = ua_b(:,:) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - umask(:,:,1) )
890	va_b (:,:) = va_b(:,:) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - vmask(:,:,1) )
891	ENDIF
892	!
893	DO jk = 1, jpkm1
894	! Correct velocities:
895	un(:,:,jk) = ( un(:,:,jk) + un_adv(:,:) - un_b(:,:) )*umask(:,:,jk)
896	vn(:,:,jk) = ( vn(:,:,jk) + vn_adv(:,:) - vn_b(:,:) )*vmask(:,:,jk)
897	!
898	END DO
899	!
900	#if defined key_agrif
901	! Save time integrated fluxes during child grid integration
902	! (used to update coarse grid transports)
903	! Useless with 2nd order momentum schemes
904	!
905	IF ( (.NOT.Agrif_Root()).AND.(ln_bt_fw) ) THEN
906	IF ( Agrif_NbStepint().EQ.0 ) THEN
907	ub2_i_b(:,:) = 0.e0
908	vb2_i_b(:,:) = 0.e0
909	END IF
910	!
911	za1 = 1._wp / REAL(Agrif_rhot(), wp)
912	ub2_i_b(:,:) = ub2_i_b(:,:) + za1 * ub2_b(:,:)
913	vb2_i_b(:,:) = vb2_i_b(:,:) + za1 * vb2_b(:,:)
914	ENDIF
915	!
916	!
917	#endif
918	!
919	! !* write time-spliting arrays in the restart
920	IF(lrst_oce .AND.ln_bt_fw) CALL ts_rst( kt, 'WRITE' )
921	!
922	CALL wrk_dealloc( jpi, jpj, zsshp2_e, zhdiv )
923	CALL wrk_dealloc( jpi, jpj, zu_trd, zv_trd, zun_e, zvn_e )
924	CALL wrk_dealloc( jpi, jpj, zwx, zwy, zu_sum, zv_sum, zssh_frc, zu_frc, zv_frc )
925	CALL wrk_dealloc( jpi, jpj, zhup2_e, zhvp2_e, zhust_e, zhvst_e )
926	CALL wrk_dealloc( jpi, jpj, zsshu_a, zsshv_a )
927	CALL wrk_dealloc( jpi, jpj, zhf )
928	!
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	INTEGER :: ios ! Local integer output status for namelist read
1073	REAL(wp) :: zxr2, zyr2, zcmax
1074	REAL(wp), POINTER, DIMENSION(:,:) :: zcu
1075	!!
1076	NAMELIST/namsplit/ ln_bt_fw, ln_bt_av, ln_bt_nn_auto, &
1077	& nn_baro, rn_bt_cmax, nn_bt_flt
1078	!!----------------------------------------------------------------------
1079	!
1080	REWIND( numnam_ref ) ! Namelist namsplit in reference namelist : time splitting parameters
1081	READ ( numnam_ref, namsplit, IOSTAT = ios, ERR = 901)
1082	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsplit in reference namelist', lwp )
1083
1084	REWIND( numnam_cfg ) ! Namelist namsplit in configuration namelist : time splitting parameters
1085	READ ( numnam_cfg, namsplit, IOSTAT = ios, ERR = 902 )
1086	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsplit in configuration namelist', lwp )
1087	IF(lwm) WRITE ( numond, namsplit )
1088	!
1089	! ! Max courant number for ext. grav. waves
1090	!
1091	CALL wrk_alloc( jpi, jpj, zcu )
1092	!
1093	IF (lk_vvl) THEN
1094	DO jj = 1, jpj
1095	DO ji =1, jpi
1096	zxr2 = 1./(e1t(ji,jj)*e1t(ji,jj))
1097	zyr2 = 1./(e2t(ji,jj)*e2t(ji,jj))
1098	zcu(ji,jj) = sqrt(gravht_0(ji,jj)(zxr2 + zyr2) )
1099	END DO
1100	END DO
1101	ELSE
1102	DO jj = 1, jpj
1103	DO ji =1, jpi
1104	zxr2 = 1./(e1t(ji,jj)*e1t(ji,jj))
1105	zyr2 = 1./(e2t(ji,jj)*e2t(ji,jj))
1106	zcu(ji,jj) = sqrt(gravht(ji,jj)(zxr2 + zyr2) )
1107	END DO
1108	END DO
1109	ENDIF
1110
1111	zcmax = MAXVAL(zcu(:,:))
1112	IF( lk_mpp ) CALL mpp_max( zcmax )
1113
1114	! Estimate number of iterations to satisfy a max courant number= rn_bt_cmax
1115	IF (ln_bt_nn_auto) nn_baro = CEILING( rdt / rn_bt_cmax * zcmax)
1116
1117	rdtbt = rdt / FLOAT(nn_baro)
1118	zcmax = zcmax * rdtbt
1119	! Print results
1120	IF(lwp) WRITE(numout,*)
1121	IF(lwp) WRITE(numout,*) 'dyn_spg_ts : split-explicit free surface'
1122	IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
1123	IF( ln_bt_nn_auto ) THEN
1124	IF(lwp) WRITE(numout,*) ' ln_ts_nn_auto=.true. Automatically set nn_baro '
1125	IF(lwp) WRITE(numout,*) ' Max. courant number allowed: ', rn_bt_cmax
1126	ELSE
1127	IF(lwp) WRITE(numout,*) ' ln_ts_nn_auto=.false.: Use nn_baro in namelist '
1128	ENDIF
1129
1130	IF(ln_bt_av) THEN
1131	IF(lwp) WRITE(numout,*) ' ln_bt_av=.true. => Time averaging over nn_baro time steps is on '
1132	ELSE
1133	IF(lwp) WRITE(numout,*) ' ln_bt_av=.false. => No time averaging of barotropic variables '
1134	ENDIF
1135	!
1136	!
1137	IF(ln_bt_fw) THEN
1138	IF(lwp) WRITE(numout,*) ' ln_bt_fw=.true. => Forward integration of barotropic variables '
1139	ELSE
1140	IF(lwp) WRITE(numout,*) ' ln_bt_fw =.false.=> Centred integration of barotropic variables '
1141	ENDIF
1142	!
1143	#if defined key_agrif
1144	! Restrict the use of Agrif to the forward case only
1145	IF ((.NOT.ln_bt_fw ).AND.(.NOT.Agrif_Root())) CALL ctl_stop( 'AGRIF not implemented if ln_bt_fw=.FALSE.' )
1146	#endif
1147	!
1148	IF(lwp) WRITE(numout,*) ' Time filter choice, nn_bt_flt: ', nn_bt_flt
1149	SELECT CASE ( nn_bt_flt )
1150	CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' Dirac'
1151	CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' Boxcar: width = nn_baro'
1152	CASE( 2 ) ; IF(lwp) WRITE(numout,) ' Boxcar: width = 2nn_baro'
1153	CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_bt_flt: should 0,1,2' )
1154	END SELECT
1155	!
1156	IF(lwp) WRITE(numout,*) ' '
1157	IF(lwp) WRITE(numout,*) ' nn_baro = ', nn_baro
1158	IF(lwp) WRITE(numout,*) ' Barotropic time step [s] is :', rdtbt
1159	IF(lwp) WRITE(numout,*) ' Maximum Courant number is :', zcmax
1160	!
1161	IF ((.NOT.ln_bt_av).AND.(.NOT.ln_bt_fw)) THEN
1162	CALL ctl_stop( 'dynspg_ts ERROR: No time averaging => only forward integration is possible' )
1163	ENDIF
1164	IF ( zcmax>0.9_wp ) THEN
1165	CALL ctl_stop( 'dynspg_ts ERROR: Maximum Courant number is greater than 0.9: Inc. nn_baro !' )
1166	ENDIF
1167	!
1168	CALL wrk_dealloc( jpi, jpj, zcu )
1169	!
1170	END SUBROUTINE dyn_spg_ts_init
1171
1172	#else
1173	!!---------------------------------------------------------------------------
1174	!! Default case : Empty module No split explicit free surface
1175	!!---------------------------------------------------------------------------
1176	CONTAINS
1177	INTEGER FUNCTION dyn_spg_ts_alloc() ! Dummy function
1178	dyn_spg_ts_alloc = 0
1179	END FUNCTION dyn_spg_ts_alloc
1180	SUBROUTINE dyn_spg_ts( kt ) ! Empty routine
1181	INTEGER, INTENT(in) :: kt
1182	WRITE(,) 'dyn_spg_ts: You should not have seen this print! error?', kt
1183	END SUBROUTINE dyn_spg_ts
1184	SUBROUTINE ts_rst( kt, cdrw ) ! Empty routine
1185	INTEGER , INTENT(in) :: kt ! ocean time-step
1186	CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag
1187	WRITE(,) 'ts_rst : You should not have seen this print! error?', kt, cdrw
1188	END SUBROUTINE ts_rst
1189	SUBROUTINE dyn_spg_ts_init( kt ) ! Empty routine
1190	INTEGER , INTENT(in) :: kt ! ocean time-step
1191	WRITE(,) 'dyn_spg_ts_init : You should not have seen this print! error?', kt
1192	END SUBROUTINE dyn_spg_ts_init
1193	#endif
1194
1195	!!======================================================================
1196	END MODULE dynspg_ts
1197
1198
1199

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

Download in other formats: