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

Context Navigation

source: branches/2013/dev_r3867_MERCATOR1_DYN/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 3970

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

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

Download in other formats: