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

Context Navigation

source: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 4354

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

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

Download in other formats: