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

Context Navigation

source: branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 4254

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

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

Download in other formats: