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 @ 4263

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

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

Download in other formats: