Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

dynspg_ts.F90 @ 12784

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

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

New URL for NEMO forge! http://forge.nemo-ocean.eu

Context Navigation

source: branches/UKMO/AMM15_v3_6_STABLE_package_reanalysis4/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 12784

Download in other formats: