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

Context Navigation

source: branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 9045

Last change on this file since 9045 was 9045, checked in by timgraham, 7 years ago
Add exlicit loops
Property svn:keywords set to `Id`
File size: 71.5 KB

Line
1	MODULE dynspg_ts
2
3	!! Includes ROMS wd scheme with diagnostic outputs ; un and ua updates are commented out !
4
5	!!======================================================================
6	!! * MODULE dynspg_ts *
7	!! Ocean dynamics: surface pressure gradient trend, split-explicit scheme
8	!!======================================================================
9	!! History : 1.0 ! 2004-12 (L. Bessieres, G. Madec) Original code
10	!! - ! 2005-11 (V. Garnier, G. Madec) optimization
11	!! - ! 2006-08 (S. Masson) distributed restart using iom
12	!! 2.0 ! 2007-07 (D. Storkey) calls to BDY routines
13	!! - ! 2008-01 (R. Benshila) change averaging method
14	!! 3.2 ! 2009-07 (R. Benshila, G. Madec) Complete revisit associated to vvl reactivation
15	!! 3.3 ! 2010-09 (D. Storkey, E. O'Dea) update for BDY for Shelf configurations
16	!! 3.3 ! 2011-03 (R. Benshila, R. Hordoir, P. Oddo) update calculation of ub_b
17	!! 3.5 ! 2013-07 (J. Chanut) Switch to Forward-backward time stepping
18	!! 3.6 ! 2013-11 (A. Coward) Update for z-tilde compatibility
19	!! 3.7 ! 2015-11 (J. Chanut) free surface simplification
20	!! - ! 2016-12 (G. Madec, E. Clementi) update for Stoke-Drift divergence
21	!! 4.0 ! 2017-05 (G. Madec) drag coef. defined at t-point (zdfdrg.F90)
22	!!---------------------------------------------------------------------
23
24	!!----------------------------------------------------------------------
25	!! dyn_spg_ts : compute surface pressure gradient trend using a time-splitting scheme
26	!! dyn_spg_ts_init: initialisation of the time-splitting scheme
27	!! ts_wgt : set time-splitting weights for temporal averaging (or not)
28	!! ts_rst : read/write time-splitting fields in restart file
29	!!----------------------------------------------------------------------
30	USE oce ! ocean dynamics and tracers
31	USE dom_oce ! ocean space and time domain
32	USE sbc_oce ! surface boundary condition: ocean
33	USE zdf_oce ! vertical physics: variables
34	USE zdfdrg ! vertical physics: top/bottom drag coef.
35	USE sbcisf ! ice shelf variable (fwfisf)
36	USE sbcapr ! surface boundary condition: atmospheric pressure
37	USE dynadv , ONLY: ln_dynadv_vec
38	USE phycst ! physical constants
39	USE dynvor ! vorticity term
40	USE wet_dry ! wetting/drying flux limter
41	USE bdy_oce ! open boundary
42	USE bdytides ! open boundary condition data
43	USE bdydyn2d ! open boundary conditions on barotropic variables
44	USE sbctide ! tides
45	USE updtide ! tide potential
46	USE sbcwave ! surface wave
47	USE diatmb ! Top,middle,bottom output
48	#if defined key_agrif
49	USE agrif_opa_interp ! agrif
50	#endif
51	#if defined key_asminc
52	USE asminc ! Assimilation increment
53	#endif
54	!
55	USE in_out_manager ! I/O manager
56	USE lib_mpp ! distributed memory computing library
57	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
58	USE prtctl ! Print control
59	USE iom ! IOM library
60	USE restart ! only for lrst_oce
61	USE timing ! Timing
62	USE diatmb ! Top,middle,bottom output
63	#if defined key_agrif
64	USE agrif_opa_interp ! agrif
65	USE agrif_oce
66	#endif
67	#if defined key_asminc
68	USE asminc ! Assimilation increment
69	#endif
70
71	IMPLICIT NONE
72	PRIVATE
73
74	PUBLIC dyn_spg_ts ! routine called in dynspg.F90
75	PUBLIC dyn_spg_ts_alloc ! " " " "
76	PUBLIC dyn_spg_ts_init ! " " " "
77	PUBLIC ts_rst ! " " " "
78
79	!! Time filtered arrays at baroclinic time step:
80	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: un_adv , vn_adv !: Advection vel. at "now" barocl. step
81	INTEGER, SAVE :: icycle ! Number of barotropic sub-steps for each internal step nn_baro <= 2.5 nn_baro
82	REAL(wp),SAVE :: rdtbt ! Barotropic time step
83	!
84	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: wgtbtp1, wgtbtp2 ! 1st & 2nd weights used in time filtering of barotropic fields
85	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zwz ! ff_f/h at F points
86	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ftnw, ftne ! triad of coriolis parameter
87	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ftsw, ftse ! (only used with een vorticity scheme)
88
89	REAL(wp) :: r1_12 = 1._wp / 12._wp ! local ratios
90	REAL(wp) :: r1_8 = 0.125_wp !
91	REAL(wp) :: r1_4 = 0.25_wp !
92	REAL(wp) :: r1_2 = 0.5_wp !
93
94	!! * Substitutions
95	# include "vectopt_loop_substitute.h90"
96	!!----------------------------------------------------------------------
97	!! NEMO/OPA 4.0 , NEMO Consortium (2017)
98	!! $Id$
99	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
100	!!----------------------------------------------------------------------
101	CONTAINS
102
103	INTEGER FUNCTION dyn_spg_ts_alloc()
104	!!----------------------------------------------------------------------
105	!! * routine dyn_spg_ts_alloc *
106	!!----------------------------------------------------------------------
107	INTEGER :: ierr(3)
108	!!----------------------------------------------------------------------
109	ierr(:) = 0
110	!
111	ALLOCATE( wgtbtp1(3nn_baro), wgtbtp2(3nn_baro), zwz(jpi,jpj), STAT=ierr(1) )
112	!
113	IF( ln_dynvor_een ) ALLOCATE( ftnw(jpi,jpj) , ftne(jpi,jpj) , &
114	& ftsw(jpi,jpj) , ftse(jpi,jpj) , STAT=ierr(2) )
115	!
116	ALLOCATE( un_adv(jpi,jpj), vn_adv(jpi,jpj) , STAT=ierr(3) )
117	!
118	dyn_spg_ts_alloc = MAXVAL( ierr(:) )
119	!
120	IF( lk_mpp ) CALL mpp_sum( dyn_spg_ts_alloc )
121	IF( dyn_spg_ts_alloc /= 0 ) CALL ctl_warn('dyn_spg_ts_alloc: failed to allocate arrays')
122	!
123	END FUNCTION dyn_spg_ts_alloc
124
125
126	SUBROUTINE dyn_spg_ts( kt )
127	!!----------------------------------------------------------------------
128	!!
129	!! ** Purpose : - Compute the now trend due to the explicit time stepping
130	!! of the quasi-linear barotropic system, and add it to the
131	!! general momentum trend.
132	!!
133	!! ** Method : - split-explicit schem (time splitting) :
134	!! Barotropic variables are advanced from internal time steps
135	!! "n" to "n+1" if ln_bt_fw=T
136	!! or from
137	!! "n-1" to "n+1" if ln_bt_fw=F
138	!! thanks to a generalized forward-backward time stepping (see ref. below).
139	!!
140	!! ** Action :
141	!! -Update the filtered free surface at step "n+1" : ssha
142	!! -Update filtered barotropic velocities at step "n+1" : ua_b, va_b
143	!! -Compute barotropic advective fluxes at step "n" : un_adv, vn_adv
144	!! These are used to advect tracers and are compliant with discrete
145	!! continuity equation taken at the baroclinic time steps. This
146	!! ensures tracers conservation.
147	!! - (ua, va) momentum trend updated with barotropic component.
148	!!
149	!! References : Shchepetkin and McWilliams, Ocean Modelling, 2005.
150	!!---------------------------------------------------------------------
151	INTEGER, INTENT(in) :: kt ! ocean time-step index
152	!
153	INTEGER :: ji, jj, jk, jn ! dummy loop indices
154	LOGICAL :: ll_fw_start ! =T : forward integration
155	LOGICAL :: ll_init ! =T : special startup of 2d equations
156	LOGICAL :: ll_tmp1, ll_tmp2 ! local logical variables used in W/D
157	INTEGER :: ikbu, iktu, noffset ! local integers
158	INTEGER :: ikbv, iktv ! - -
159	REAL(wp) :: r1_2dt_b, z2dt_bf ! local scalars
160	REAL(wp) :: zx1, zx2, zu_spg, zhura ! - -
161	REAL(wp) :: zy1, zy2, zv_spg, zhvra ! - -
162	REAL(wp) :: za0, za1, za2, za3 ! - -
163	REAL(wp) :: zmdi, zztmp ! - -
164	REAL(wp), DIMENSION(jpi,jpj) :: zsshp2_e, zhf
165	REAL(wp), DIMENSION(jpi,jpj) :: zwx, zu_trd, zu_frc, zssh_frc
166	REAL(wp), DIMENSION(jpi,jpj) :: zwy, zv_trd, zv_frc, zhdiv
167	REAL(wp), DIMENSION(jpi,jpj) :: zsshu_a, zhup2_e, zhust_e
168	REAL(wp), DIMENSION(jpi,jpj) :: zsshv_a, zhvp2_e, zhvst_e
169	REAL(wp), DIMENSION(jpi,jpj) :: zCdU_u, zCdU_v ! top/bottom stress at u- & v-points
170	!
171	REAL(wp) :: zwdramp ! local scalar - only used if ln_wd_dl = .True.
172
173	INTEGER :: iwdg, jwdg, kwdg ! short-hand values for the indices of the output point
174
175	REAL(wp) :: zepsilon, zgamma ! - -
176	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zcpx, zcpy ! Wetting/Dying gravity filter coef.
177	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: ztwdmask, zuwdmask, zvwdmask ! ROMS wetting and drying masks at t,u,v points
178	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zuwdav2, zvwdav2 ! averages over the sub-steps of zuwdmask and zvwdmask
179	!!----------------------------------------------------------------------
180	!
181	IF( ln_timing ) CALL timing_start('dyn_spg_ts')
182	!
183	IF( ln_wd_il ) ALLOCATE( zcpx(jpi,jpj), zcpy(jpi,jpj) )
184	! !* Allocate temporary arrays
185	IF( ln_wd_dl ) ALLOCATE( ztwdmask(jpi,jpj), zuwdmask(jpi,jpj), zvwdmask(jpi,jpj), zuwdav2(jpi,jpj), zvwdav2(jpi,jpj))
186	!
187	zmdi=1.e+20 ! missing data indicator for masking
188	!
189	zwdramp = r_rn_wdmin1 ! simplest ramp
190	! zwdramp = 1._wp / (rn_wdmin2 - rn_wdmin1) ! more general ramp
191	! ! reciprocal of baroclinic time step
192	IF( kt == nit000 .AND. neuler == 0 ) THEN ; z2dt_bf = rdt
193	ELSE ; z2dt_bf = 2.0_wp * rdt
194	ENDIF
195	r1_2dt_b = 1.0_wp / z2dt_bf
196	!
197	ll_init = ln_bt_av ! if no time averaging, then no specific restart
198	ll_fw_start = .FALSE.
199	! ! time offset in steps for bdy data update
200	IF( .NOT.ln_bt_fw ) THEN ; noffset = - nn_baro
201	ELSE ; noffset = 0
202	ENDIF
203	!
204	IF( kt == nit000 ) THEN !* initialisation
205	!
206	IF(lwp) WRITE(numout,*)
207	IF(lwp) WRITE(numout,*) 'dyn_spg_ts : surface pressure gradient trend'
208	IF(lwp) WRITE(numout,*) '~~~~~~~~~~ free surface with time splitting'
209	IF(lwp) WRITE(numout,*)
210	!
211	IF( neuler == 0 ) ll_init=.TRUE.
212	!
213	IF( ln_bt_fw .OR. neuler == 0 ) THEN
214	ll_fw_start =.TRUE.
215	noffset = 0
216	ELSE
217	ll_fw_start =.FALSE.
218	ENDIF
219	!
220	! Set averaging weights and cycle length:
221	CALL ts_wgt( ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2 )
222	!
223	ENDIF
224	!
225	IF( ln_isfcav ) THEN ! top+bottom friction (ocean cavities)
226	DO jj = 2, jpjm1
227	DO ji = fs_2, fs_jpim1 ! vector opt.
228	zCdU_u(ji,jj) = 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) + rCdU_top(ji+1,jj)+rCdU_top(ji,jj) )
229	zCdU_v(ji,jj) = 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) + rCdU_top(ji,jj+1)+rCdU_top(ji,jj) )
230	END DO
231	END DO
232	ELSE ! bottom friction only
233	DO jj = 2, jpjm1
234	DO ji = fs_2, fs_jpim1 ! vector opt.
235	zCdU_u(ji,jj) = 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) )
236	zCdU_v(ji,jj) = 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) )
237	END DO
238	END DO
239	ENDIF
240	!
241	! Set arrays to remove/compute coriolis trend.
242	! Do it once at kt=nit000 if volume is fixed, else at each long time step.
243	! Note that these arrays are also used during barotropic loop. These are however frozen
244	! although they should be updated in the variable volume case. Not a big approximation.
245	! To remove this approximation, copy lines below inside barotropic loop
246	! and update depths at T-F points (ht and zhf resp.) at each barotropic time step
247	!
248	IF( kt == nit000 .OR. .NOT.ln_linssh ) THEN
249	IF( ln_dynvor_een ) THEN !== EEN scheme ==!
250	SELECT CASE( nn_een_e3f ) !* ff_f/e3 at F-point
251	CASE ( 0 ) ! original formulation (masked averaging of e3t divided by 4)
252	DO jj = 1, jpjm1
253	DO ji = 1, jpim1
254	zwz(ji,jj) = ( ht_n(ji ,jj+1) + ht_n(ji+1,jj+1) + &
255	& ht_n(ji ,jj ) + ht_n(ji+1,jj ) ) * 0.25_wp
256	IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj)
257	END DO
258	END DO
259	CASE ( 1 ) ! new formulation (masked averaging of e3t divided by the sum of mask)
260	DO jj = 1, jpjm1
261	DO ji = 1, jpim1
262	zwz(ji,jj) = ( ht_n(ji ,jj+1) + ht_n(ji+1,jj+1) + &
263	& ht_n(ji ,jj ) + ht_n(ji+1,jj ) ) &
264	& / ( MAX( 1._wp, tmask(ji ,jj+1, 1) + tmask(ji+1,jj+1, 1) + &
265	& tmask(ji ,jj , 1) + tmask(ji+1,jj , 1) ) )
266	IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj)
267	END DO
268	END DO
269	END SELECT
270	CALL lbc_lnk( zwz, 'F', 1._wp )
271	!
272	ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp
273	DO jj = 2, jpj
274	DO ji = 2, jpi
275	ftne(ji,jj) = zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1)
276	ftnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj ) + zwz(ji ,jj )
277	ftse(ji,jj) = zwz(ji ,jj ) + zwz(ji ,jj-1) + zwz(ji-1,jj-1)
278	ftsw(ji,jj) = zwz(ji ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj )
279	END DO
280	END DO
281	!
282	ELSE !== all other schemes (ENE, ENS, MIX)
283	zwz(:,:) = 0._wp
284	zhf(:,:) = 0._wp
285
286	!!gm assume 0 in both cases (xhich is almost surely WRONG ! ) as hvatf has been removed
287	!!gm A priori a better value should be something like :
288	!!gm zhf(i,j) = masked sum of ht(i,j) , ht(i+1,j) , ht(i,j+1) , (i+1,j+1)
289	!!gm divided by the sum of the corresponding mask
290	!!gm
291	!!
292	IF( .NOT.ln_sco ) THEN
293
294	!!gm agree the JC comment : this should be done in a much clear way
295
296	! JC: It not clear yet what should be the depth at f-points over land in z-coordinate case
297	! Set it to zero for the time being
298	! IF( rn_hmin < 0._wp ) THEN ; jk = - INT( rn_hmin ) ! from a nb of level
299	! ELSE ; jk = MINLOC( gdepw_0, mask = gdepw_0 > rn_hmin, dim = 1 ) ! from a depth
300	! ENDIF
301	! zhf(:,:) = gdepw_0(:,:,jk+1)
302	ELSE
303	!zhf(:,:) = hbatf(:,:)
304	DO jj = 1, jpjm1
305	DO ji = 1, jpim1
306	zhf(ji,jj) = MAX( 0._wp, &
307	& ( ht_0(ji ,jj )*tmask(ji ,jj ,1) + &
308	& ht_0(ji+1,jj )*tmask(ji+1,jj ,1) + &
309	& ht_0(ji ,jj+1)*tmask(ji ,jj+1,1) + &
310	& ht_0(ji+1,jj+1)*tmask(ji+1,jj+1,1) ) / &
311	& ( tmask(ji ,jj ,1) + tmask(ji+1,jj ,1) +&
312	& tmask(ji ,jj+1,1) + tmask(ji+1,jj+1,1) +&
313	& rsmall ) )
314	END DO
315	END DO
316	END IF
317
318	DO jj = 1, jpjm1
319	zhf(:,jj) = zhf(:,jj) * (1._wp- umask(:,jj,1) * umask(:,jj+1,1))
320	END DO
321	!!gm end
322
323	DO jk = 1, jpkm1
324	DO jj = 1, jpjm1
325	zhf(:,jj) = zhf(:,jj) + e3f_n(:,jj,jk) * umask(:,jj,jk) * umask(:,jj+1,jk)
326	END DO
327	END DO
328	CALL lbc_lnk( zhf, 'F', 1._wp )
329	! JC: TBC. hf should be greater than 0
330	DO jj = 1, jpj
331	DO ji = 1, jpi
332	IF( zhf(ji,jj) /= 0._wp ) zwz(ji,jj) = 1._wp / zhf(ji,jj) ! zhf is actually hf here but it saves an array
333	END DO
334	END DO
335	zwz(:,:) = ff_f(:,:) * zwz(:,:)
336	ENDIF
337	ENDIF
338	!
339	! If forward start at previous time step, and centered integration,
340	! then update averaging weights:
341	IF (.NOT.ln_bt_fw .AND.( neuler==0 .AND. kt==nit000+1 ) ) THEN
342	ll_fw_start=.FALSE.
343	CALL ts_wgt( ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2 )
344	ENDIF
345
346	! -----------------------------------------------------------------------------
347	! Phase 1 : Coupling between general trend and barotropic estimates (1st step)
348	! -----------------------------------------------------------------------------
349	!
350	!
351	! !* e3*d/dt(Ua) (Vertically integrated)
352	! ! --------------------------------------------------
353	zu_frc(:,:) = 0._wp
354	zv_frc(:,:) = 0._wp
355	!
356	DO jk = 1, jpkm1
357	zu_frc(:,:) = zu_frc(:,:) + e3u_n(:,:,jk) * ua(:,:,jk) * umask(:,:,jk)
358	zv_frc(:,:) = zv_frc(:,:) + e3v_n(:,:,jk) * va(:,:,jk) * vmask(:,:,jk)
359	END DO
360	!
361	zu_frc(:,:) = zu_frc(:,:) * r1_hu_n(:,:)
362	zv_frc(:,:) = zv_frc(:,:) * r1_hv_n(:,:)
363	!
364	!
365	! !* baroclinic momentum trend (remove the vertical mean trend)
366	DO jk = 1, jpkm1 ! -----------------------------------------------------------
367	DO jj = 2, jpjm1
368	DO ji = fs_2, fs_jpim1 ! vector opt.
369	ua(ji,jj,jk) = ua(ji,jj,jk) - zu_frc(ji,jj) * umask(ji,jj,jk)
370	va(ji,jj,jk) = va(ji,jj,jk) - zv_frc(ji,jj) * vmask(ji,jj,jk)
371	END DO
372	END DO
373	END DO
374
375	!!gm Question here when removing the Vertically integrated trends, we remove the vertically integrated NL trends on momentum....
376	!!gm Is it correct to do so ? I think so...
377
378
379	! !* barotropic Coriolis trends (vorticity scheme dependent)
380	! ! --------------------------------------------------------
381	zwx(:,:) = un_b(:,:) * hu_n(:,:) * e2u(:,:) ! now fluxes
382	zwy(:,:) = vn_b(:,:) * hv_n(:,:) * e1v(:,:)
383	!
384	IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN ! energy conserving or mixed scheme
385	DO jj = 2, jpjm1
386	DO ji = fs_2, fs_jpim1 ! vector opt.
387	zy1 = ( zwy(ji,jj-1) + zwy(ji+1,jj-1) ) * r1_e1u(ji,jj)
388	zy2 = ( zwy(ji,jj ) + zwy(ji+1,jj ) ) * r1_e1u(ji,jj)
389	zx1 = ( zwx(ji-1,jj) + zwx(ji-1,jj+1) ) * r1_e2v(ji,jj)
390	zx2 = ( zwx(ji ,jj) + zwx(ji ,jj+1) ) * r1_e2v(ji,jj)
391	! energy conserving formulation for planetary vorticity term
392	zu_trd(ji,jj) = r1_4 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
393	zv_trd(ji,jj) = - r1_4 * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 )
394	END DO
395	END DO
396	!
397	ELSEIF ( ln_dynvor_ens ) THEN ! enstrophy conserving scheme
398	DO jj = 2, jpjm1
399	DO ji = fs_2, fs_jpim1 ! vector opt.
400	zy1 = r1_8 * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) &
401	& + zwy(ji ,jj ) + zwy(ji+1,jj ) ) * r1_e1u(ji,jj)
402	zx1 = - r1_8 * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) &
403	& + zwx(ji ,jj ) + zwx(ji ,jj+1) ) * r1_e2v(ji,jj)
404	zu_trd(ji,jj) = zy1 * ( zwz(ji ,jj-1) + zwz(ji,jj) )
405	zv_trd(ji,jj) = zx1 * ( zwz(ji-1,jj ) + zwz(ji,jj) )
406	END DO
407	END DO
408	!
409	ELSEIF ( ln_dynvor_een ) THEN ! enstrophy and energy conserving scheme
410	DO jj = 2, jpjm1
411	DO ji = fs_2, fs_jpim1 ! vector opt.
412	zu_trd(ji,jj) = + r1_12 * r1_e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) &
413	& + ftnw(ji+1,jj) * zwy(ji+1,jj ) &
414	& + ftse(ji,jj ) * zwy(ji ,jj-1) &
415	& + ftsw(ji+1,jj) * zwy(ji+1,jj-1) )
416	zv_trd(ji,jj) = - r1_12 * r1_e2v(ji,jj) * ( ftsw(ji,jj+1) * zwx(ji-1,jj+1) &
417	& + ftse(ji,jj+1) * zwx(ji ,jj+1) &
418	& + ftnw(ji,jj ) * zwx(ji-1,jj ) &
419	& + ftne(ji,jj ) * zwx(ji ,jj ) )
420	END DO
421	END DO
422	!
423	ENDIF
424	!
425	! !* Right-Hand-Side of the barotropic momentum equation
426	! ! ----------------------------------------------------
427	IF( .NOT.ln_linssh ) THEN ! Variable volume : remove surface pressure gradient
428	IF( ln_wd_il ) THEN ! Calculating and applying W/D gravity filters
429	DO jj = 2, jpjm1
430	DO ji = 2, jpim1
431	ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji+1,jj) ) > &
432	& MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) .AND. &
433	& MAX( sshn(ji,jj) + ht_0(ji,jj) , sshn(ji+1,jj) + ht_0(ji+1,jj) ) &
434	& > rn_wdmin1 + rn_wdmin2
435	ll_tmp2 = ( ABS( sshn(ji+1,jj) - sshn(ji ,jj)) > 1.E-12 ).AND.( &
436	& MAX( sshn(ji,jj) , sshn(ji+1,jj) ) > &
437	& MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 )
438
439	IF(ll_tmp1) THEN
440	zcpx(ji,jj) = 1.0_wp
441	ELSE IF(ll_tmp2) THEN
442	! no worries about sshn(ji+1,jj) - sshn(ji ,jj) = 0, it won't happen ! here
443	zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_0(ji+1,jj) - sshn(ji,jj) - ht_0(ji,jj)) &
444	& / (sshn(ji+1,jj) - sshn(ji ,jj)) )
445	zcpx(ji,jj) = max(min( zcpx(ji,jj) , 1.0_wp),0.0_wp)
446
447	ELSE
448	zcpx(ji,jj) = 0._wp
449	END IF
450
451	ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji,jj+1) ) > &
452	& MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) .AND. &
453	& MAX( sshn(ji,jj) + ht_0(ji,jj) , sshn(ji,jj+1) + ht_0(ji,jj+1) ) &
454	& > rn_wdmin1 + rn_wdmin2
455	ll_tmp2 = ( ABS( sshn(ji,jj) - sshn(ji,jj+1)) > 1.E-12 ).AND.( &
456	& MAX( sshn(ji,jj) , sshn(ji,jj+1) ) > &
457	& MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 )
458
459	IF(ll_tmp1) THEN
460	zcpy(ji,jj) = 1.0_wp
461	ELSE IF(ll_tmp2) THEN
462	! no worries about sshn(ji,jj+1) - sshn(ji,jj ) = 0, it won't happen ! here
463	zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_0(ji,jj+1) - sshn(ji,jj) - ht_0(ji,jj)) &
464	& / (sshn(ji,jj+1) - sshn(ji,jj )) )
465	zcpy(ji,jj) = max(min( zcpy(ji,jj) , 1.0_wp),0.0_wp)
466
467	ELSE
468	zcpy(ji,jj) = 0._wp
469	END IF
470	END DO
471	END DO
472
473	DO jj = 2, jpjm1
474	DO ji = 2, jpim1
475	zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) &
476	& * r1_e1u(ji,jj) * zcpx(ji,jj) * wdrampu(ji,jj) !jth
477	zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji ,jj+1) - sshn(ji ,jj ) ) &
478	& * r1_e2v(ji,jj) * zcpy(ji,jj) * wdrampv(ji,jj) !jth
479
480	END DO
481	END DO
482
483	ELSE
484	!
485	DO jj = 2, jpjm1
486	DO ji = fs_2, fs_jpim1 ! vector opt.
487	zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) * r1_e1u(ji,jj)
488	zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji ,jj+1) - sshn(ji ,jj ) ) * r1_e2v(ji,jj)
489	END DO
490	END DO
491	ENDIF
492	!
493	ENDIF
494	!
495	DO jj = 2, jpjm1 ! Remove coriolis term (and possibly spg) from barotropic trend
496	DO ji = fs_2, fs_jpim1
497	zu_frc(ji,jj) = zu_frc(ji,jj) - zu_trd(ji,jj) * ssumask(ji,jj)
498	zv_frc(ji,jj) = zv_frc(ji,jj) - zv_trd(ji,jj) * ssvmask(ji,jj)
499	END DO
500	END DO
501	!
502	! ! Add bottom stress contribution from baroclinic velocities:
503	IF (ln_bt_fw) THEN
504	DO jj = 2, jpjm1
505	DO ji = fs_2, fs_jpim1 ! vector opt.
506	ikbu = mbku(ji,jj)
507	ikbv = mbkv(ji,jj)
508	zwx(ji,jj) = un(ji,jj,ikbu) - un_b(ji,jj) ! NOW bottom baroclinic velocities
509	zwy(ji,jj) = vn(ji,jj,ikbv) - vn_b(ji,jj)
510	END DO
511	END DO
512	ELSE
513	DO jj = 2, jpjm1
514	DO ji = fs_2, fs_jpim1 ! vector opt.
515	ikbu = mbku(ji,jj)
516	ikbv = mbkv(ji,jj)
517	zwx(ji,jj) = ub(ji,jj,ikbu) - ub_b(ji,jj) ! BEFORE bottom baroclinic velocities
518	zwy(ji,jj) = vb(ji,jj,ikbv) - vb_b(ji,jj)
519	END DO
520	END DO
521	ENDIF
522	!
523	! Note that the "unclipped" bottom friction parameter is used even with explicit drag
524	IF( ln_wd_il ) THEN
525	DO jj = 2, jpjm1
526	DO ji = fs_2, fs_jpim1 ! vector opt.
527	zu_frc(ji,jj) = zu_frc(ji,jj) + MAX(r1_hu_n(ji,jj) * zCdU_u(ji,jj),-1._wp / rdtbt) * zwx(ji,jj) * wdrampu(ji,jj)
528	zv_frc(ji,jj) = zv_frc(ji,jj) + MAX(r1_hv_n(ji,jj) * zCdU_v(ji,jj),-1._wp / rdtbt) * zwy(ji,jj) * wdrampv(ji,jj)
529	END DO
530	END DO
531	ELSE
532	DO jj = 2, jpjm1
533	DO ji = fs_2, fs_jpim1 ! vector opt.
534	zu_frc(ji,jj) = zu_frc(ji,jj) + r1_hu_n(ji,jj) * zCdU_u(ji,jj) * zwx(ji,jj)
535	zv_frc(ji,jj) = zv_frc(ji,jj) + r1_hv_n(ji,jj) * zCdU_v(ji,jj) * zwy(ji,jj)
536	END DO
537	END DO
538	END IF
539	!
540	! ! Add top stress contribution from baroclinic velocities:
541	IF( ln_bt_fw ) THEN
542	DO jj = 2, jpjm1
543	DO ji = fs_2, fs_jpim1 ! vector opt.
544	zu_frc(ji,jj) = zu_frc(ji,jj) + MAX( r1_hu_n(ji,jj) * 0.5( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) , zztmp ) zwx(ji,jj)
545	zv_frc(ji,jj) = zv_frc(ji,jj) + MAX( r1_hv_n(ji,jj) * 0.5( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) , zztmp ) zwy(ji,jj)
546	END DO
547	END DO
548	ELSE
549	DO jj = 2, jpjm1
550	DO ji = fs_2, fs_jpim1 ! vector opt.
551	zu_frc(ji,jj) = zu_frc(ji,jj) + r1_hu_n(ji,jj) * 0.5( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) zwx(ji,jj)
552	zv_frc(ji,jj) = zv_frc(ji,jj) + r1_hv_n(ji,jj) * 0.5( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) zwy(ji,jj)
553	END DO
554	END DO
555	ENDIF
556	!
557	IF( ln_isfcav ) THEN ! Add TOP stress contribution from baroclinic velocities:
558	IF( ln_bt_fw ) THEN
559	DO jj = 2, jpjm1
560	DO ji = fs_2, fs_jpim1 ! vector opt.
561	iktu = miku(ji,jj)
562	iktv = mikv(ji,jj)
563	zwx(ji,jj) = un(ji,jj,iktu) - un_b(ji,jj) ! NOW top baroclinic velocities
564	zwy(ji,jj) = vn(ji,jj,iktv) - vn_b(ji,jj)
565	END DO
566	END DO
567	ELSE
568	DO jj = 2, jpjm1
569	DO ji = fs_2, fs_jpim1 ! vector opt.
570	iktu = miku(ji,jj)
571	iktv = mikv(ji,jj)
572	zwx(ji,jj) = ub(ji,jj,iktu) - ub_b(ji,jj) ! BEFORE top baroclinic velocities
573	zwy(ji,jj) = vb(ji,jj,iktv) - vb_b(ji,jj)
574	END DO
575	END DO
576	ENDIF
577	!
578	! Note that the "unclipped" top friction parameter is used even with explicit drag
579	DO jj = 2, jpjm1
580	DO ji = fs_2, fs_jpim1 ! vector opt.
581	zu_frc(ji,jj) = zu_frc(ji,jj) + r1_hu_n(ji,jj) * 0.5( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) zwx(ji,jj)
582	zv_frc(ji,jj) = zv_frc(ji,jj) + r1_hv_n(ji,jj) * 0.5( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) zwy(ji,jj)
583	END DO
584	END DO
585	ENDIF
586	!
587	IF( ln_bt_fw ) THEN ! Add wind forcing
588	DO jj = 2, jpjm1
589	DO ji = fs_2, fs_jpim1 ! vector opt.
590	zu_frc(ji,jj) = zu_frc(ji,jj) + r1_rau0 * utau(ji,jj) * r1_hu_n(ji,jj)
591	zv_frc(ji,jj) = zv_frc(ji,jj) + r1_rau0 * vtau(ji,jj) * r1_hv_n(ji,jj)
592	END DO
593	END DO
594	ELSE
595	zztmp = r1_rau0 * r1_2
596	DO jj = 2, jpjm1
597	DO ji = fs_2, fs_jpim1 ! vector opt.
598	zu_frc(ji,jj) = zu_frc(ji,jj) + zztmp * ( utau_b(ji,jj) + utau(ji,jj) ) * r1_hu_n(ji,jj)
599	zv_frc(ji,jj) = zv_frc(ji,jj) + zztmp * ( vtau_b(ji,jj) + vtau(ji,jj) ) * r1_hv_n(ji,jj)
600	END DO
601	END DO
602	ENDIF
603	!
604	IF( ln_apr_dyn ) THEN ! Add atm pressure forcing
605	IF( ln_bt_fw ) THEN
606	DO jj = 2, jpjm1
607	DO ji = fs_2, fs_jpim1 ! vector opt.
608	zu_spg = grav * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) ) * r1_e1u(ji,jj)
609	zv_spg = grav * ( ssh_ib (ji ,jj+1) - ssh_ib (ji,jj) ) * r1_e2v(ji,jj)
610	zu_frc(ji,jj) = zu_frc(ji,jj) + zu_spg
611	zv_frc(ji,jj) = zv_frc(ji,jj) + zv_spg
612	END DO
613	END DO
614	ELSE
615	zztmp = grav * r1_2
616	DO jj = 2, jpjm1
617	DO ji = fs_2, fs_jpim1 ! vector opt.
618	zu_spg = zztmp * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) &
619	& + ssh_ibb(ji+1,jj ) - ssh_ibb(ji,jj) ) * r1_e1u(ji,jj)
620	zv_spg = zztmp * ( ssh_ib (ji ,jj+1) - ssh_ib (ji,jj) &
621	& + ssh_ibb(ji ,jj+1) - ssh_ibb(ji,jj) ) * r1_e2v(ji,jj)
622	zu_frc(ji,jj) = zu_frc(ji,jj) + zu_spg
623	zv_frc(ji,jj) = zv_frc(ji,jj) + zv_spg
624	END DO
625	END DO
626	ENDIF
627	ENDIF
628	! !* Right-Hand-Side of the barotropic ssh equation
629	! ! -----------------------------------------------
630	! ! Surface net water flux and rivers
631	IF (ln_bt_fw) THEN
632	zssh_frc(:,:) = r1_rau0 * ( emp(:,:) - rnf(:,:) + fwfisf(:,:) )
633	ELSE
634	zztmp = r1_rau0 * r1_2
635	zssh_frc(:,:) = zztmp * ( emp(:,:) + emp_b(:,:) - rnf(:,:) - rnf_b(:,:) &
636	& + fwfisf(:,:) + fwfisf_b(:,:) )
637	ENDIF
638	!
639	IF( ln_sdw ) THEN ! Stokes drift divergence added if necessary
640	zssh_frc(:,:) = zssh_frc(:,:) + div_sd(:,:)
641	ENDIF
642	!
643	#if defined key_asminc
644	! ! Include the IAU weighted SSH increment
645	IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN
646	zssh_frc(:,:) = zssh_frc(:,:) - ssh_iau(:,:)
647	ENDIF
648	#endif
649	! !* Fill boundary data arrays for AGRIF
650	! ! ------------------------------------
651	#if defined key_agrif
652	IF( .NOT.Agrif_Root() ) CALL agrif_dta_ts( kt )
653	#endif
654	!
655	! -----------------------------------------------------------------------
656	! Phase 2 : Integration of the barotropic equations
657	! -----------------------------------------------------------------------
658	!
659	! ! ==================== !
660	! ! Initialisations !
661	! ! ==================== !
662	! Initialize barotropic variables:
663	IF( ll_init )THEN
664	sshbb_e(:,:) = 0._wp
665	ubb_e (:,:) = 0._wp
666	vbb_e (:,:) = 0._wp
667	sshb_e (:,:) = 0._wp
668	ub_e (:,:) = 0._wp
669	vb_e (:,:) = 0._wp
670	ENDIF
671
672	!
673	IF (ln_bt_fw) THEN ! FORWARD integration: start from NOW fields
674	sshn_e(:,:) = sshn(:,:)
675	un_e (:,:) = un_b(:,:)
676	vn_e (:,:) = vn_b(:,:)
677	!
678	hu_e (:,:) = hu_n(:,:)
679	hv_e (:,:) = hv_n(:,:)
680	hur_e (:,:) = r1_hu_n(:,:)
681	hvr_e (:,:) = r1_hv_n(:,:)
682	ELSE ! CENTRED integration: start from BEFORE fields
683	sshn_e(:,:) = sshb(:,:)
684	un_e (:,:) = ub_b(:,:)
685	vn_e (:,:) = vb_b(:,:)
686	!
687	hu_e (:,:) = hu_b(:,:)
688	hv_e (:,:) = hv_b(:,:)
689	hur_e (:,:) = r1_hu_b(:,:)
690	hvr_e (:,:) = r1_hv_b(:,:)
691	ENDIF
692	!
693	!
694	!
695	! Initialize sums:
696	ua_b (:,:) = 0._wp ! After barotropic velocities (or transport if flux form)
697	va_b (:,:) = 0._wp
698	ssha (:,:) = 0._wp ! Sum for after averaged sea level
699	un_adv(:,:) = 0._wp ! Sum for now transport issued from ts loop
700	vn_adv(:,:) = 0._wp
701	! ! ==================== !
702
703	IF (ln_wd_dl) THEN
704	zuwdmask(:,:) = 0._wp ! set to zero for definiteness (not sure this is necessary)
705	zvwdmask(:,:) = 0._wp !
706	zuwdav2(:,:) = 0._wp
707	zvwdav2(:,:) = 0._wp
708	END IF
709
710
711	DO jn = 1, icycle ! sub-time-step loop !
712	! ! ==================== !
713	! !* Update the forcing (BDY and tides)
714	! ! ------------------
715	! Update only tidal forcing at open boundaries
716	IF( ln_bdy .AND. ln_tide ) CALL bdy_dta_tides( kt, kit=jn, time_offset= noffset+1 )
717	IF( ln_tide_pot .AND. ln_tide ) CALL upd_tide ( kt, kit=jn, time_offset= noffset )
718	!
719	! Set extrapolation coefficients for predictor step:
720	IF ((jn<3).AND.ll_init) THEN ! Forward
721	za1 = 1._wp
722	za2 = 0._wp
723	za3 = 0._wp
724	ELSE ! AB3-AM4 Coefficients: bet=0.281105
725	za1 = 1.781105_wp ! za1 = 3/2 + bet
726	za2 = -1.06221_wp ! za2 = -(1/2 + 2*bet)
727	za3 = 0.281105_wp ! za3 = bet
728	ENDIF
729
730	! Extrapolate barotropic velocities at step jit+0.5:
731	ua_e(:,:) = za1 * un_e(:,:) + za2 * ub_e(:,:) + za3 * ubb_e(:,:)
732	va_e(:,:) = za1 * vn_e(:,:) + za2 * vb_e(:,:) + za3 * vbb_e(:,:)
733
734	IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only)
735	! ! ------------------
736	! Extrapolate Sea Level at step jit+0.5:
737	zsshp2_e(:,:) = za1 * sshn_e(:,:) + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:)
738
739	! set wetting & drying mask at tracer points for this barotropic sub-step
740	IF ( ln_wd_dl ) THEN
741
742	IF ( ln_wd_dl_rmp ) THEN
743	DO jj = 1, jpj
744	DO ji = 1, jpi ! vector opt.
745	IF ( zsshp2_e(ji,jj) + ht_0(ji,jj) > 2._wp * rn_wdmin1 ) THEN
746	! IF ( zsshp2_e(ji,jj) + ht_0(ji,jj) > rn_wdmin2 ) THEN
747	ztwdmask(ji,jj) = 1._wp
748	ELSE IF ( zsshp2_e(ji,jj) + ht_0(ji,jj) > rn_wdmin1 ) THEN
749	ztwdmask(ji,jj) = (tanh(50._wp( ( zsshp2_e(ji,jj) + ht_0(ji,jj) - rn_wdmin1 )r_rn_wdmin1)) )
750	ELSE
751	ztwdmask(ji,jj) = 0._wp
752	END IF
753	END DO
754	END DO
755	ELSE
756	DO jj = 1, jpj
757	DO ji = 1, jpi ! vector opt.
758	IF ( zsshp2_e(ji,jj) + ht_0(ji,jj) > rn_wdmin1 ) THEN
759	ztwdmask(ji,jj) = 1._wp
760	ELSE
761	ztwdmask(ji,jj) = 0._wp
762	END IF
763	END DO
764	END DO
765	END IF
766
767	END IF
768
769
770	DO jj = 2, jpjm1 ! Sea Surface Height at u- & v-points
771	DO ji = 2, fs_jpim1 ! Vector opt.
772	zwx(ji,jj) = r1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj) &
773	& * ( e1e2t(ji ,jj) * zsshp2_e(ji ,jj) &
774	& + e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj) )
775	zwy(ji,jj) = r1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj) &
776	& * ( e1e2t(ji,jj ) * zsshp2_e(ji,jj ) &
777	& + e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1) )
778	END DO
779	END DO
780	CALL lbc_lnk_multi( zwx, 'U', 1._wp, zwy, 'V', 1._wp )
781	!
782	zhup2_e (:,:) = hu_0(:,:) + zwx(:,:) ! Ocean depth at U- and V-points
783	zhvp2_e (:,:) = hv_0(:,:) + zwy(:,:)
784	ELSE
785	zhup2_e (:,:) = hu_n(:,:)
786	zhvp2_e (:,:) = hv_n(:,:)
787	ENDIF
788	! !* after ssh
789	! ! -----------
790	! One should enforce volume conservation at open boundaries here
791	! considering fluxes below:
792	!
793	zwx(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:) ! fluxes at jn+0.5
794	zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:)
795	!
796	#if defined key_agrif
797	! Set fluxes during predictor step to ensure volume conservation
798	IF( .NOT.Agrif_Root() .AND. ln_bt_fw ) THEN
799	IF((nbondi == -1).OR.(nbondi == 2)) THEN
800	DO jj = 1, jpj
801	zwx(2:nbghostcells+1,jj) = ubdy_w(jj) * e2u(2:nbghostcells+1,jj)
802	END DO
803	ENDIF
804	IF((nbondi == 1).OR.(nbondi == 2)) THEN
805	DO jj=1,jpj
806	zwx(nlci-nbghostcells-1:nlci-2,jj) = ubdy_e(jj) * e2u(nlci-nbghostcells-1:nlci-2,jj)
807	END DO
808	ENDIF
809	IF((nbondj == -1).OR.(nbondj == 2)) THEN
810	DO ji=1,jpi
811	zwy(ji,2:nbghostcells+1) = vbdy_s(ji) * e1v(ji,2:nbghostcells+1)
812	END DO
813	ENDIF
814	IF((nbondj == 1).OR.(nbondj == 2)) THEN
815	DO ji=1,jpi
816	zwy(ji,nlcj-nbghostcells-1:nlcj-2) = vbdy_n(ji) * e1v(ji,nlcj-nbghostcells-1:nlcj-2)
817	END DO
818	ENDIF
819	ENDIF
820	#endif
821	IF( ln_wd_il ) CALL wad_lmt_bt(zwx, zwy, sshn_e, zssh_frc, rdtbt)
822
823	IF ( ln_wd_dl ) THEN
824
825
826	! un_e and vn_e are set to zero at faces where the direction of the flow is from dry cells
827
828	DO jj = 1, jpjm1
829	DO ji = 1, jpim1
830	IF ( zwx(ji,jj) > 0.0 ) THEN
831	zuwdmask(ji, jj) = ztwdmask(ji ,jj)
832	ELSE
833	zuwdmask(ji, jj) = ztwdmask(ji+1,jj)
834	END IF
835	zwx(ji, jj) = zuwdmask(ji,jj)*zwx(ji, jj)
836	un_e(ji,jj) = zuwdmask(ji,jj)*un_e(ji,jj)
837
838	IF ( zwy(ji,jj) > 0.0 ) THEN
839	zvwdmask(ji, jj) = ztwdmask(ji, jj )
840	ELSE
841	zvwdmask(ji, jj) = ztwdmask(ji, jj+1)
842	END IF
843	zwy(ji, jj) = zvwdmask(ji,jj)*zwy(ji,jj)
844	vn_e(ji,jj) = zvwdmask(ji,jj)*vn_e(ji,jj)
845	END DO
846	END DO
847
848
849	END IF
850
851	! Sum over sub-time-steps to compute advective velocities
852	za2 = wgtbtp2(jn)
853	un_adv(:,:) = un_adv(:,:) + za2 * zwx(:,:) * r1_e2u(:,:)
854	vn_adv(:,:) = vn_adv(:,:) + za2 * zwy(:,:) * r1_e1v(:,:)
855
856	! sum over sub-time-steps to decide which baroclinic velocities to set to zero (zuwdav2 is only used when ln_wd_dl_bc = True)
857	IF ( ln_wd_dl_bc ) THEN
858	zuwdav2(:,:) = zuwdav2(:,:) + za2 * zuwdmask(:,:)
859	zvwdav2(:,:) = zvwdav2(:,:) + za2 * zvwdmask(:,:)
860	END IF
861
862	! Set next sea level:
863	DO jj = 2, jpjm1
864	DO ji = fs_2, fs_jpim1 ! vector opt.
865	zhdiv(ji,jj) = ( zwx(ji,jj) - zwx(ji-1,jj) &
866	& + zwy(ji,jj) - zwy(ji,jj-1) ) * r1_e1e2t(ji,jj)
867	END DO
868	END DO
869	ssha_e(:,:) = ( sshn_e(:,:) - rdtbt * ( zssh_frc(:,:) + zhdiv(:,:) ) ) * ssmask(:,:)
870
871	CALL lbc_lnk( ssha_e, 'T', 1._wp )
872
873	! Duplicate sea level across open boundaries (this is only cosmetic if linssh=T)
874	IF( ln_bdy ) CALL bdy_ssh( ssha_e )
875	#if defined key_agrif
876	IF( .NOT.Agrif_Root() ) CALL agrif_ssh_ts( jn )
877	#endif
878	!
879	! Sea Surface Height at u-,v-points (vvl case only)
880	IF( .NOT.ln_linssh ) THEN
881	DO jj = 2, jpjm1
882	DO ji = 2, jpim1 ! NO Vector Opt.
883	zsshu_a(ji,jj) = r1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj) &
884	& * ( e1e2t(ji ,jj ) * ssha_e(ji ,jj ) &
885	& + e1e2t(ji+1,jj ) * ssha_e(ji+1,jj ) )
886	zsshv_a(ji,jj) = r1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj) &
887	& * ( e1e2t(ji ,jj ) * ssha_e(ji ,jj ) &
888	& + e1e2t(ji ,jj+1) * ssha_e(ji ,jj+1) )
889	END DO
890	END DO
891	CALL lbc_lnk_multi( zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp )
892	ENDIF
893	!
894	! Half-step back interpolation of SSH for surface pressure computation:
895	!----------------------------------------------------------------------
896	IF ((jn==1).AND.ll_init) THEN
897	za0=1._wp ! Forward-backward
898	za1=0._wp
899	za2=0._wp
900	za3=0._wp
901	ELSEIF ((jn==2).AND.ll_init) THEN ! AB2-AM3 Coefficients; bet=0 ; gam=-1/6 ; eps=1/12
902	za0= 1.0833333333333_wp ! za0 = 1-gam-eps
903	za1=-0.1666666666666_wp ! za1 = gam
904	za2= 0.0833333333333_wp ! za2 = eps
905	za3= 0._wp
906	ELSE ! AB3-AM4 Coefficients; bet=0.281105 ; eps=0.013 ; gam=0.0880
907	IF (rn_bt_alpha==0._wp) THEN
908	za0=0.614_wp ! za0 = 1/2 + gam + 2*eps
909	za1=0.285_wp ! za1 = 1/2 - 2gam - 3eps
910	za2=0.088_wp ! za2 = gam
911	za3=0.013_wp ! za3 = eps
912	ELSE
913	zepsilon = 0.00976186_wp - 0.13451357_wp * rn_bt_alpha
914	zgamma = 0.08344500_wp - 0.51358400_wp * rn_bt_alpha
915	za0 = 0.5_wp + zgamma + 2._wp * rn_bt_alpha + 2._wp * zepsilon
916	za1 = 1._wp - za0 - zgamma - zepsilon
917	za2 = zgamma
918	za3 = zepsilon
919	ENDIF
920	ENDIF
921	!
922	zsshp2_e(:,:) = za0 * ssha_e(:,:) + za1 * sshn_e (:,:) &
923	& + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:)
924	IF( ln_wd_il ) THEN ! Calculating and applying W/D gravity filters
925	DO jj = 2, jpjm1
926	DO ji = 2, jpim1
927	ll_tmp1 = MIN( zsshp2_e(ji,jj) , zsshp2_e(ji+1,jj) ) > &
928	& MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) .AND. &
929	& MAX( zsshp2_e(ji,jj) + ht_0(ji,jj) , zsshp2_e(ji+1,jj) + ht_0(ji+1,jj) ) &
930	& > rn_wdmin1 + rn_wdmin2
931	ll_tmp2 = (ABS(zsshp2_e(ji,jj) - zsshp2_e(ji+1,jj)) > 1.E-12 ).AND.( &
932	& MAX( zsshp2_e(ji,jj) , zsshp2_e(ji+1,jj) ) > &
933	& MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 )
934
935	IF(ll_tmp1) THEN
936	zcpx(ji,jj) = 1.0_wp
937	ELSE IF(ll_tmp2) THEN
938	! no worries about zsshp2_e(ji+1,jj) - zsshp2_e(ji ,jj) = 0, it won't happen ! here
939	zcpx(ji,jj) = ABS( (zsshp2_e(ji+1,jj) + ht_0(ji+1,jj) - zsshp2_e(ji,jj) - ht_0(ji,jj)) &
940	& / (zsshp2_e(ji+1,jj) - zsshp2_e(ji ,jj)) )
941	ELSE
942	zcpx(ji,jj) = 0._wp
943	END IF
944
945	ll_tmp1 = MIN( zsshp2_e(ji,jj) , zsshp2_e(ji,jj+1) ) > &
946	& MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) .AND. &
947	& MAX( zsshp2_e(ji,jj) + ht_0(ji,jj) , zsshp2_e(ji,jj+1) + ht_0(ji,jj+1) ) &
948	& > rn_wdmin1 + rn_wdmin2
949	ll_tmp2 = (ABS(zsshp2_e(ji,jj) - zsshp2_e(ji,jj+1)) > 1.E-12 ).AND.( &
950	& MAX( zsshp2_e(ji,jj) , zsshp2_e(ji,jj+1) ) > &
951	& MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 )
952
953	IF(ll_tmp1) THEN
954	zcpy(ji,jj) = 1.0_wp
955	ELSE IF(ll_tmp2) THEN
956	! no worries about zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj ) = 0, it won't happen ! here
957	zcpy(ji,jj) = ABS( (zsshp2_e(ji,jj+1) + ht_0(ji,jj+1) - zsshp2_e(ji,jj) - ht_0(ji,jj)) &
958	& / (zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj )) )
959	ELSE
960	zcpy(ji,jj) = 0._wp
961	END IF
962	END DO
963	END DO
964	END IF
965	!
966	! Compute associated depths at U and V points:
967	IF( .NOT.ln_linssh .AND. .NOT.ln_dynadv_vec ) THEN !* Vector form
968	!
969	DO jj = 2, jpjm1
970	DO ji = 2, jpim1
971	zx1 = r1_2 * ssumask(ji ,jj) * r1_e1e2u(ji ,jj) &
972	& * ( e1e2t(ji ,jj ) * zsshp2_e(ji ,jj) &
973	& + e1e2t(ji+1,jj ) * zsshp2_e(ji+1,jj ) )
974	zy1 = r1_2 * ssvmask(ji ,jj) * r1_e1e2v(ji ,jj ) &
975	& * ( e1e2t(ji ,jj ) * zsshp2_e(ji ,jj ) &
976	& + e1e2t(ji ,jj+1) * zsshp2_e(ji ,jj+1) )
977	zhust_e(ji,jj) = hu_0(ji,jj) + zx1
978	zhvst_e(ji,jj) = hv_0(ji,jj) + zy1
979	END DO
980	END DO
981
982	ENDIF
983	!
984	! Add Coriolis trend:
985	! zwz array below or triads normally depend on sea level with ln_linssh=F and should be updated
986	! at each time step. We however keep them constant here for optimization.
987	! Recall that zwx and zwy arrays hold fluxes at this stage:
988	! zwx(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:) ! fluxes at jn+0.5
989	! zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:)
990	!
991	IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN !== energy conserving or mixed scheme ==!
992	DO jj = 2, jpjm1
993	DO ji = fs_2, fs_jpim1 ! vector opt.
994	zy1 = ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) ) * r1_e1u(ji,jj)
995	zy2 = ( zwy(ji ,jj ) + zwy(ji+1,jj ) ) * r1_e1u(ji,jj)
996	zx1 = ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) ) * r1_e2v(ji,jj)
997	zx2 = ( zwx(ji ,jj ) + zwx(ji ,jj+1) ) * r1_e2v(ji,jj)
998	zu_trd(ji,jj) = r1_4 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
999	zv_trd(ji,jj) =-r1_4 * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 )
1000	END DO
1001	END DO
1002	!
1003	ELSEIF ( ln_dynvor_ens ) THEN !== enstrophy conserving scheme ==!
1004	DO jj = 2, jpjm1
1005	DO ji = fs_2, fs_jpim1 ! vector opt.
1006	zy1 = r1_8 * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) &
1007	& + zwy(ji ,jj ) + zwy(ji+1,jj ) ) * r1_e1u(ji,jj)
1008	zx1 = - r1_8 * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) &
1009	& + zwx(ji ,jj ) + zwx(ji ,jj+1) ) * r1_e2v(ji,jj)
1010	zu_trd(ji,jj) = zy1 * ( zwz(ji ,jj-1) + zwz(ji,jj) )
1011	zv_trd(ji,jj) = zx1 * ( zwz(ji-1,jj ) + zwz(ji,jj) )
1012	END DO
1013	END DO
1014	!
1015	ELSEIF ( ln_dynvor_een ) THEN !== energy and enstrophy conserving scheme ==!
1016	DO jj = 2, jpjm1
1017	DO ji = fs_2, fs_jpim1 ! vector opt.
1018	zu_trd(ji,jj) = + r1_12 * r1_e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) &
1019	& + ftnw(ji+1,jj) * zwy(ji+1,jj ) &
1020	& + ftse(ji,jj ) * zwy(ji ,jj-1) &
1021	& + ftsw(ji+1,jj) * zwy(ji+1,jj-1) )
1022	zv_trd(ji,jj) = - r1_12 * r1_e2v(ji,jj) * ( ftsw(ji,jj+1) * zwx(ji-1,jj+1) &
1023	& + ftse(ji,jj+1) * zwx(ji ,jj+1) &
1024	& + ftnw(ji,jj ) * zwx(ji-1,jj ) &
1025	& + ftne(ji,jj ) * zwx(ji ,jj ) )
1026	END DO
1027	END DO
1028	!
1029	ENDIF
1030	!
1031	! Add tidal astronomical forcing if defined
1032	IF ( ln_tide .AND. ln_tide_pot ) THEN
1033	DO jj = 2, jpjm1
1034	DO ji = fs_2, fs_jpim1 ! vector opt.
1035	zu_spg = grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) * r1_e1u(ji,jj)
1036	zv_spg = grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) * r1_e2v(ji,jj)
1037	zu_trd(ji,jj) = zu_trd(ji,jj) + zu_spg
1038	zv_trd(ji,jj) = zv_trd(ji,jj) + zv_spg
1039	END DO
1040	END DO
1041	ENDIF
1042	!
1043	! Add bottom stresses:
1044	!jth do implicitly instead
1045	IF ( .NOT. ll_wd ) THEN ! Revert to explicit for bit comparison tests in non wad runs
1046	DO jj = 2, jpjm1
1047	DO ji = fs_2, fs_jpim1 ! vector opt.
1048	zu_trd(ji,jj) = zu_trd(ji,jj) + zCdU_u(ji,jj) * un_e(ji,jj) * hur_e(ji,jj)
1049	zv_trd(ji,jj) = zv_trd(ji,jj) + zCdU_v(ji,jj) * vn_e(ji,jj) * hvr_e(ji,jj)
1050	END DO
1051	END DO
1052	ENDIF
1053	!
1054	! Surface pressure trend:
1055	IF( ln_wd_il ) THEN
1056	DO jj = 2, jpjm1
1057	DO ji = 2, jpim1
1058	! Add surface pressure gradient
1059	zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj)
1060	zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj)
1061	zwx(ji,jj) = (1._wp - rn_scal_load) * zu_spg * zcpx(ji,jj)
1062	zwy(ji,jj) = (1._wp - rn_scal_load) * zv_spg * zcpy(ji,jj)
1063	END DO
1064	END DO
1065	ELSE
1066	DO jj = 2, jpjm1
1067	DO ji = fs_2, fs_jpim1 ! vector opt.
1068	! Add surface pressure gradient
1069	zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj)
1070	zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj)
1071	zwx(ji,jj) = (1._wp - rn_scal_load) * zu_spg
1072	zwy(ji,jj) = (1._wp - rn_scal_load) * zv_spg
1073	END DO
1074	END DO
1075	END IF
1076
1077	!
1078	! Set next velocities:
1079	IF( ln_dynadv_vec .OR. ln_linssh ) THEN !* Vector form
1080	DO jj = 2, jpjm1
1081	DO ji = fs_2, fs_jpim1 ! vector opt.
1082	ua_e(ji,jj) = ( un_e(ji,jj) &
1083	& + rdtbt * ( zwx(ji,jj) &
1084	& + zu_trd(ji,jj) &
1085	& + zu_frc(ji,jj) ) &
1086	& ) * ssumask(ji,jj)
1087
1088	va_e(ji,jj) = ( vn_e(ji,jj) &
1089	& + rdtbt * ( zwy(ji,jj) &
1090	& + zv_trd(ji,jj) &
1091	& + zv_frc(ji,jj) ) &
1092	& ) * ssvmask(ji,jj)
1093
1094	!jth implicit bottom friction:
1095	IF ( ll_wd ) THEN ! revert to explicit for bit comparison tests in non wad runs
1096	ua_e(ji,jj) = ua_e(ji,jj) /(1.0 - rdtbt * zCdU_u(ji,jj) * hur_e(ji,jj))
1097	va_e(ji,jj) = va_e(ji,jj) /(1.0 - rdtbt * zCdU_v(ji,jj) * hvr_e(ji,jj))
1098	ENDIF
1099
1100	END DO
1101	END DO
1102	!
1103	ELSE !* Flux form
1104	DO jj = 2, jpjm1
1105	DO ji = fs_2, fs_jpim1 ! vector opt.
1106
1107	zhura = hu_0(ji,jj) + zsshu_a(ji,jj)
1108	zhvra = hv_0(ji,jj) + zsshv_a(ji,jj)
1109
1110	zhura = ssumask(ji,jj)/(zhura + 1._wp - ssumask(ji,jj))
1111	zhvra = ssvmask(ji,jj)/(zhvra + 1._wp - ssvmask(ji,jj))
1112
1113	ua_e(ji,jj) = ( hu_e(ji,jj) * un_e(ji,jj) &
1114	& + rdtbt * ( zhust_e(ji,jj) * zwx(ji,jj) &
1115	& + zhup2_e(ji,jj) * zu_trd(ji,jj) &
1116	& + hu_n(ji,jj) * zu_frc(ji,jj) ) &
1117	& ) * zhura
1118
1119	va_e(ji,jj) = ( hv_e(ji,jj) * vn_e(ji,jj) &
1120	& + rdtbt * ( zhvst_e(ji,jj) * zwy(ji,jj) &
1121	& + zhvp2_e(ji,jj) * zv_trd(ji,jj) &
1122	& + hv_n(ji,jj) * zv_frc(ji,jj) ) &
1123	& ) * zhvra
1124	END DO
1125	END DO
1126	ENDIF
1127
1128
1129	IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only)
1130	hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:)
1131	hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:)
1132	hur_e(:,:) = ssumask(:,:) / ( hu_e(:,:) + 1._wp - ssumask(:,:) )
1133	hvr_e(:,:) = ssvmask(:,:) / ( hv_e(:,:) + 1._wp - ssvmask(:,:) )
1134	!
1135	ENDIF
1136	! !* domain lateral boundary
1137	CALL lbc_lnk_multi( ua_e, 'U', -1._wp, va_e , 'V', -1._wp )
1138	!
1139	! ! open boundaries
1140	IF( ln_bdy ) CALL bdy_dyn2d( jn, ua_e, va_e, un_e, vn_e, hur_e, hvr_e, ssha_e )
1141	#if defined key_agrif
1142	IF( .NOT.Agrif_Root() ) CALL agrif_dyn_ts( jn ) ! Agrif
1143	#endif
1144	! !* Swap
1145	! ! ----
1146	ubb_e (:,:) = ub_e (:,:)
1147	ub_e (:,:) = un_e (:,:)
1148	un_e (:,:) = ua_e (:,:)
1149	!
1150	vbb_e (:,:) = vb_e (:,:)
1151	vb_e (:,:) = vn_e (:,:)
1152	vn_e (:,:) = va_e (:,:)
1153	!
1154	sshbb_e(:,:) = sshb_e(:,:)
1155	sshb_e (:,:) = sshn_e(:,:)
1156	sshn_e (:,:) = ssha_e(:,:)
1157
1158	! !* Sum over whole bt loop
1159	! ! ----------------------
1160	za1 = wgtbtp1(jn)
1161	IF( ln_dynadv_vec .OR. ln_linssh ) THEN ! Sum velocities
1162	ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:)
1163	va_b (:,:) = va_b (:,:) + za1 * va_e (:,:)
1164	ELSE ! Sum transports
1165	IF ( .NOT.ln_wd_dl ) THEN
1166	ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:) * hu_e (:,:)
1167	va_b (:,:) = va_b (:,:) + za1 * va_e (:,:) * hv_e (:,:)
1168	ELSE
1169	ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:) * hu_e (:,:) * zuwdmask(:,:)
1170	va_b (:,:) = va_b (:,:) + za1 * va_e (:,:) * hv_e (:,:) * zvwdmask(:,:)
1171	END IF
1172	ENDIF
1173	! ! Sum sea level
1174	ssha(:,:) = ssha(:,:) + za1 * ssha_e(:,:)
1175
1176	! ! ==================== !
1177	END DO ! end loop !
1178	! ! ==================== !
1179	! -----------------------------------------------------------------------------
1180	! Phase 3. update the general trend with the barotropic trend
1181	! -----------------------------------------------------------------------------
1182	!
1183	! Set advection velocity correction:
1184	IF (ln_bt_fw) THEN
1185	zwx(:,:) = un_adv(:,:)
1186	zwy(:,:) = vn_adv(:,:)
1187	IF( .NOT.( kt == nit000 .AND. neuler==0 ) ) THEN
1188	un_adv(:,:) = r1_2 * ( ub2_b(:,:) + zwx(:,:) - atfp * un_bf(:,:) )
1189	vn_adv(:,:) = r1_2 * ( vb2_b(:,:) + zwy(:,:) - atfp * vn_bf(:,:) )
1190	!
1191	! Update corrective fluxes for next time step:
1192	un_bf(:,:) = atfp * un_bf(:,:) + (zwx(:,:) - ub2_b(:,:))
1193	vn_bf(:,:) = atfp * vn_bf(:,:) + (zwy(:,:) - vb2_b(:,:))
1194	ELSE
1195	un_bf(:,:) = 0._wp
1196	vn_bf(:,:) = 0._wp
1197	END IF
1198	! Save integrated transport for next computation
1199	ub2_b(:,:) = zwx(:,:)
1200	vb2_b(:,:) = zwy(:,:)
1201	ENDIF
1202
1203
1204	!
1205	! Update barotropic trend:
1206	IF( ln_dynadv_vec .OR. ln_linssh ) THEN
1207	DO jk=1,jpkm1
1208	ua(:,:,jk) = ua(:,:,jk) + ( ua_b(:,:) - ub_b(:,:) ) * r1_2dt_b
1209	va(:,:,jk) = va(:,:,jk) + ( va_b(:,:) - vb_b(:,:) ) * r1_2dt_b
1210	END DO
1211	ELSE
1212	! At this stage, ssha has been corrected: compute new depths at velocity points
1213	DO jj = 1, jpjm1
1214	DO ji = 1, jpim1 ! NO Vector Opt.
1215	zsshu_a(ji,jj) = r1_2 * umask(ji,jj,1) * r1_e1e2u(ji,jj) &
1216	& * ( e1e2t(ji ,jj) * ssha(ji ,jj) &
1217	& + e1e2t(ji+1,jj) * ssha(ji+1,jj) )
1218	zsshv_a(ji,jj) = r1_2 * vmask(ji,jj,1) * r1_e1e2v(ji,jj) &
1219	& * ( e1e2t(ji,jj ) * ssha(ji,jj ) &
1220	& + e1e2t(ji,jj+1) * ssha(ji,jj+1) )
1221	END DO
1222	END DO
1223	CALL lbc_lnk_multi( zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) ! Boundary conditions
1224	!
1225	DO jk=1,jpkm1
1226	ua(:,:,jk) = ua(:,:,jk) + r1_hu_n(:,:) * ( ua_b(:,:) - ub_b(:,:) * hu_b(:,:) ) * r1_2dt_b
1227	va(:,:,jk) = va(:,:,jk) + r1_hv_n(:,:) * ( va_b(:,:) - vb_b(:,:) * hv_b(:,:) ) * r1_2dt_b
1228	END DO
1229	! Save barotropic velocities not transport:
1230	ua_b(:,:) = ua_b(:,:) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - ssumask(:,:) )
1231	va_b(:,:) = va_b(:,:) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - ssvmask(:,:) )
1232	ENDIF
1233
1234
1235	! Correct velocities so that the barotropic velocity equals (un_adv, vn_adv) (in all cases)
1236	DO jk = 1, jpkm1
1237	un(:,:,jk) = ( un(:,:,jk) + un_adv(:,:)r1_hu_n(:,:) - un_b(:,:) ) umask(:,:,jk)
1238	vn(:,:,jk) = ( vn(:,:,jk) + vn_adv(:,:)r1_hv_n(:,:) - vn_b(:,:) ) vmask(:,:,jk)
1239	END DO
1240
1241	IF ( ln_wd_dl .and. ln_wd_dl_bc) THEN
1242	DO jk = 1, jpkm1
1243	un(:,:,jk) = ( un_adv(:,:) + zuwdav2(:,:)(un(:,:,jk) - un_adv(:,:)) ) umask(:,:,jk)
1244	vn(:,:,jk) = ( vn_adv(:,:) + zvwdav2(:,:)(vn(:,:,jk) - vn_adv(:,:)) ) vmask(:,:,jk)
1245	END DO
1246	END IF
1247
1248
1249	CALL iom_put( "ubar", un_adv(:,:)*r1_hu_n(:,:) ) ! barotropic i-current
1250	CALL iom_put( "vbar", vn_adv(:,:)*r1_hv_n(:,:) ) ! barotropic i-current
1251	!
1252	#if defined key_agrif
1253	! Save time integrated fluxes during child grid integration
1254	! (used to update coarse grid transports at next time step)
1255	!
1256	IF( .NOT.Agrif_Root() .AND. ln_bt_fw ) THEN
1257	IF( Agrif_NbStepint() == 0 ) THEN
1258	ub2_i_b(:,:) = 0._wp
1259	vb2_i_b(:,:) = 0._wp
1260	END IF
1261	!
1262	za1 = 1._wp / REAL(Agrif_rhot(), wp)
1263	ub2_i_b(:,:) = ub2_i_b(:,:) + za1 * ub2_b(:,:)
1264	vb2_i_b(:,:) = vb2_i_b(:,:) + za1 * vb2_b(:,:)
1265	ENDIF
1266	#endif
1267	! !* write time-spliting arrays in the restart
1268	IF( lrst_oce .AND.ln_bt_fw ) CALL ts_rst( kt, 'WRITE' )
1269	!
1270	IF( ln_wd_il ) DEALLOCATE( zcpx, zcpy )
1271	IF( ln_wd_dl ) DEALLOCATE( ztwdmask, zuwdmask, zvwdmask, zuwdav2, zvwdav2 )
1272	!
1273	IF( ln_diatmb ) THEN
1274	CALL iom_put( "baro_u" , un_bumask(:,:,1)+zmdi(1-umask(:,:,1 ) ) ) ! Barotropic U Velocity
1275	CALL iom_put( "baro_v" , vn_bvmask(:,:,1)+zmdi(1-vmask(:,:,1 ) ) ) ! Barotropic V Velocity
1276	ENDIF
1277	IF( ln_timing ) CALL timing_stop('dyn_spg_ts')
1278	!
1279	END SUBROUTINE dyn_spg_ts
1280
1281
1282	SUBROUTINE ts_wgt( ll_av, ll_fw, jpit, zwgt1, zwgt2)
1283	!!---------------------------------------------------------------------
1284	!! * ROUTINE ts_wgt *
1285	!!
1286	!! ** Purpose : Set time-splitting weights for temporal averaging (or not)
1287	!!----------------------------------------------------------------------
1288	LOGICAL, INTENT(in) :: ll_av ! temporal averaging=.true.
1289	LOGICAL, INTENT(in) :: ll_fw ! forward time splitting =.true.
1290	INTEGER, INTENT(inout) :: jpit ! cycle length
1291	REAL(wp), DIMENSION(3*nn_baro), INTENT(inout) :: zwgt1, & ! Primary weights
1292	zwgt2 ! Secondary weights
1293
1294	INTEGER :: jic, jn, ji ! temporary integers
1295	REAL(wp) :: za1, za2
1296	!!----------------------------------------------------------------------
1297
1298	zwgt1(:) = 0._wp
1299	zwgt2(:) = 0._wp
1300
1301	! Set time index when averaged value is requested
1302	IF (ll_fw) THEN
1303	jic = nn_baro
1304	ELSE
1305	jic = 2 * nn_baro
1306	ENDIF
1307
1308	! Set primary weights:
1309	IF (ll_av) THEN
1310	! Define simple boxcar window for primary weights
1311	! (width = nn_baro, centered around jic)
1312	SELECT CASE ( nn_bt_flt )
1313	CASE( 0 ) ! No averaging
1314	zwgt1(jic) = 1._wp
1315	jpit = jic
1316
1317	CASE( 1 ) ! Boxcar, width = nn_baro
1318	DO jn = 1, 3*nn_baro
1319	za1 = ABS(float(jn-jic))/float(nn_baro)
1320	IF (za1 < 0.5_wp) THEN
1321	zwgt1(jn) = 1._wp
1322	jpit = jn
1323	ENDIF
1324	ENDDO
1325
1326	CASE( 2 ) ! Boxcar, width = 2 * nn_baro
1327	DO jn = 1, 3*nn_baro
1328	za1 = ABS(float(jn-jic))/float(nn_baro)
1329	IF (za1 < 1._wp) THEN
1330	zwgt1(jn) = 1._wp
1331	jpit = jn
1332	ENDIF
1333	ENDDO
1334	CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_bt_flt' )
1335	END SELECT
1336
1337	ELSE ! No time averaging
1338	zwgt1(jic) = 1._wp
1339	jpit = jic
1340	ENDIF
1341
1342	! Set secondary weights
1343	DO jn = 1, jpit
1344	DO ji = jn, jpit
1345	zwgt2(jn) = zwgt2(jn) + zwgt1(ji)
1346	END DO
1347	END DO
1348
1349	! Normalize weigths:
1350	za1 = 1._wp / SUM(zwgt1(1:jpit))
1351	za2 = 1._wp / SUM(zwgt2(1:jpit))
1352	DO jn = 1, jpit
1353	zwgt1(jn) = zwgt1(jn) * za1
1354	zwgt2(jn) = zwgt2(jn) * za2
1355	END DO
1356	!
1357	END SUBROUTINE ts_wgt
1358
1359
1360	SUBROUTINE ts_rst( kt, cdrw )
1361	!!---------------------------------------------------------------------
1362	!! * ROUTINE ts_rst *
1363	!!
1364	!! ** Purpose : Read or write time-splitting arrays in restart file
1365	!!----------------------------------------------------------------------
1366	INTEGER , INTENT(in) :: kt ! ocean time-step
1367	CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag
1368	!
1369	!!----------------------------------------------------------------------
1370	!
1371	IF( TRIM(cdrw) == 'READ' ) THEN
1372	CALL iom_get( numror, jpdom_autoglo, 'ub2_b' , ub2_b (:,:) )
1373	CALL iom_get( numror, jpdom_autoglo, 'vb2_b' , vb2_b (:,:) )
1374	CALL iom_get( numror, jpdom_autoglo, 'un_bf' , un_bf (:,:) )
1375	CALL iom_get( numror, jpdom_autoglo, 'vn_bf' , vn_bf (:,:) )
1376	IF( .NOT.ln_bt_av ) THEN
1377	CALL iom_get( numror, jpdom_autoglo, 'sshbb_e' , sshbb_e(:,:) )
1378	CALL iom_get( numror, jpdom_autoglo, 'ubb_e' , ubb_e(:,:) )
1379	CALL iom_get( numror, jpdom_autoglo, 'vbb_e' , vbb_e(:,:) )
1380	CALL iom_get( numror, jpdom_autoglo, 'sshb_e' , sshb_e(:,:) )
1381	CALL iom_get( numror, jpdom_autoglo, 'ub_e' , ub_e(:,:) )
1382	CALL iom_get( numror, jpdom_autoglo, 'vb_e' , vb_e(:,:) )
1383	ENDIF
1384	#if defined key_agrif
1385	! Read time integrated fluxes
1386	IF ( .NOT.Agrif_Root() ) THEN
1387	CALL iom_get( numror, jpdom_autoglo, 'ub2_i_b' , ub2_i_b(:,:) )
1388	CALL iom_get( numror, jpdom_autoglo, 'vb2_i_b' , vb2_i_b(:,:) )
1389	ENDIF
1390	#endif
1391	!
1392	ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN
1393	CALL iom_rstput( kt, nitrst, numrow, 'ub2_b' , ub2_b (:,:) )
1394	CALL iom_rstput( kt, nitrst, numrow, 'vb2_b' , vb2_b (:,:) )
1395	CALL iom_rstput( kt, nitrst, numrow, 'un_bf' , un_bf (:,:) )
1396	CALL iom_rstput( kt, nitrst, numrow, 'vn_bf' , vn_bf (:,:) )
1397	!
1398	IF (.NOT.ln_bt_av) THEN
1399	CALL iom_rstput( kt, nitrst, numrow, 'sshbb_e' , sshbb_e(:,:) )
1400	CALL iom_rstput( kt, nitrst, numrow, 'ubb_e' , ubb_e(:,:) )
1401	CALL iom_rstput( kt, nitrst, numrow, 'vbb_e' , vbb_e(:,:) )
1402	CALL iom_rstput( kt, nitrst, numrow, 'sshb_e' , sshb_e(:,:) )
1403	CALL iom_rstput( kt, nitrst, numrow, 'ub_e' , ub_e(:,:) )
1404	CALL iom_rstput( kt, nitrst, numrow, 'vb_e' , vb_e(:,:) )
1405	ENDIF
1406	#if defined key_agrif
1407	! Save time integrated fluxes
1408	IF ( .NOT.Agrif_Root() ) THEN
1409	CALL iom_rstput( kt, nitrst, numrow, 'ub2_i_b' , ub2_i_b(:,:) )
1410	CALL iom_rstput( kt, nitrst, numrow, 'vb2_i_b' , vb2_i_b(:,:) )
1411	ENDIF
1412	#endif
1413	ENDIF
1414	!
1415	END SUBROUTINE ts_rst
1416
1417
1418	SUBROUTINE dyn_spg_ts_init
1419	!!---------------------------------------------------------------------
1420	!! * ROUTINE dyn_spg_ts_init *
1421	!!
1422	!! ** Purpose : Set time splitting options
1423	!!----------------------------------------------------------------------
1424	INTEGER :: ji ,jj ! dummy loop indices
1425	REAL(wp) :: zxr2, zyr2, zcmax ! local scalar
1426	REAL(wp), DIMENSION(jpi,jpj) :: zcu
1427	!!----------------------------------------------------------------------
1428	!
1429	! Max courant number for ext. grav. waves
1430	!
1431	DO jj = 1, jpj
1432	DO ji =1, jpi
1433	zxr2 = r1_e1t(ji,jj) * r1_e1t(ji,jj)
1434	zyr2 = r1_e2t(ji,jj) * r1_e2t(ji,jj)
1435	zcu(ji,jj) = SQRT( grav * MAX(ht_0(ji,jj),0._wp) * (zxr2 + zyr2) )
1436	END DO
1437	END DO
1438	!
1439	zcmax = MAXVAL( zcu(:,:) )
1440	IF( lk_mpp ) CALL mpp_max( zcmax )
1441
1442	! Estimate number of iterations to satisfy a max courant number= rn_bt_cmax
1443	IF( ln_bt_auto ) nn_baro = CEILING( rdt / rn_bt_cmax * zcmax)
1444
1445	rdtbt = rdt / REAL( nn_baro , wp )
1446	zcmax = zcmax * rdtbt
1447	! Print results
1448	IF(lwp) WRITE(numout,*)
1449	IF(lwp) WRITE(numout,*) 'dyn_spg_ts : split-explicit free surface'
1450	IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
1451	IF( ln_bt_auto ) THEN
1452	IF(lwp) WRITE(numout,*) ' ln_ts_auto=.true. Automatically set nn_baro '
1453	IF(lwp) WRITE(numout,*) ' Max. courant number allowed: ', rn_bt_cmax
1454	ELSE
1455	IF(lwp) WRITE(numout,*) ' ln_ts_auto=.false.: Use nn_baro in namelist '
1456	ENDIF
1457
1458	IF(ln_bt_av) THEN
1459	IF(lwp) WRITE(numout,*) ' ln_bt_av=.true. => Time averaging over nn_baro time steps is on '
1460	ELSE
1461	IF(lwp) WRITE(numout,*) ' ln_bt_av=.false. => No time averaging of barotropic variables '
1462	ENDIF
1463	!
1464	!
1465	IF(ln_bt_fw) THEN
1466	IF(lwp) WRITE(numout,*) ' ln_bt_fw=.true. => Forward integration of barotropic variables '
1467	ELSE
1468	IF(lwp) WRITE(numout,*) ' ln_bt_fw =.false.=> Centred integration of barotropic variables '
1469	ENDIF
1470	!
1471	#if defined key_agrif
1472	! Restrict the use of Agrif to the forward case only
1473	!!! IF( .NOT.ln_bt_fw .AND. .NOT.Agrif_Root() ) CALL ctl_stop( 'AGRIF not implemented if ln_bt_fw=.FALSE.' )
1474	#endif
1475	!
1476	IF(lwp) WRITE(numout,*) ' Time filter choice, nn_bt_flt: ', nn_bt_flt
1477	SELECT CASE ( nn_bt_flt )
1478	CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' Dirac'
1479	CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' Boxcar: width = nn_baro'
1480	CASE( 2 ) ; IF(lwp) WRITE(numout,) ' Boxcar: width = 2nn_baro'
1481	CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_bt_flt: should 0,1,2' )
1482	END SELECT
1483	!
1484	IF(lwp) WRITE(numout,*) ' '
1485	IF(lwp) WRITE(numout,*) ' nn_baro = ', nn_baro
1486	IF(lwp) WRITE(numout,*) ' Barotropic time step [s] is :', rdtbt
1487	IF(lwp) WRITE(numout,*) ' Maximum Courant number is :', zcmax
1488	!
1489	IF(lwp) WRITE(numout,*) ' Time diffusion parameter rn_bt_alpha: ', rn_bt_alpha
1490	IF ((ln_bt_av.AND.nn_bt_flt/=0).AND.(rn_bt_alpha>0._wp)) THEN
1491	CALL ctl_stop( 'dynspg_ts ERROR: if rn_bt_alpha > 0, remove temporal averaging' )
1492	ENDIF
1493	!
1494	IF( .NOT.ln_bt_av .AND. .NOT.ln_bt_fw ) THEN
1495	CALL ctl_stop( 'dynspg_ts ERROR: No time averaging => only forward integration is possible' )
1496	ENDIF
1497	IF( zcmax>0.9_wp ) THEN
1498	CALL ctl_stop( 'dynspg_ts ERROR: Maximum Courant number is greater than 0.9: Inc. nn_baro !' )
1499	ENDIF
1500	!
1501	END SUBROUTINE dyn_spg_ts_init
1502
1503	!!======================================================================
1504	END MODULE dynspg_ts

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

Download in other formats: