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

Context Navigation

source: branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 6069

Last change on this file since 6069 was 6069, checked in by timgraham, 8 years ago
Merge of dev_MetOffice_merge_2015 into branch (only NEMO directory for now).
Property svn:keywords set to `Id`
File size: 54.8 KB

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

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

Download in other formats: