Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

dynnxt.F90 @ 15422

Last change on this file since 15422 was 15422, checked in by jcastill, 3 years ago
Changes tested so that they can merged with the CO9 Met Office branch - jpmax_harmo should be 34 with FES14 tides, but the last components are not used anyway
File size: 22.7 KB

Line
1	MODULE dynnxt
2	!!=========================================================================
3	!! * MODULE dynnxt *
4	!! Ocean dynamics: time stepping
5	!!=========================================================================
6	!! History : OPA ! 1987-02 (P. Andrich, D. L Hostis) Original code
7	!! ! 1990-10 (C. Levy, G. Madec)
8	!! 7.0 ! 1993-03 (M. Guyon) symetrical conditions
9	!! 8.0 ! 1997-02 (G. Madec & M. Imbard) opa, release 8.0
10	!! 8.2 ! 1997-04 (A. Weaver) Euler forward step
11	!! - ! 1997-06 (G. Madec) lateral boudary cond., lbc routine
12	!! NEMO 1.0 ! 2002-08 (G. Madec) F90: Free form and module
13	!! - ! 2002-10 (C. Talandier, A-M. Treguier) Open boundary cond.
14	!! 2.0 ! 2005-11 (V. Garnier) Surface pressure gradient organization
15	!! 2.3 ! 2007-07 (D. Storkey) Calls to BDY routines.
16	!! 3.2 ! 2009-06 (G. Madec, R.Benshila) re-introduce the vvl option
17	!! 3.3 ! 2010-09 (D. Storkey, E.O'Dea) Bug fix for BDY module
18	!! 3.3 ! 2011-03 (P. Oddo) Bug fix for time-splitting+(BDY-OBC) and not VVL
19	!! 3.5 ! 2013-07 (J. Chanut) Compliant with time splitting changes
20	!! 3.6 ! 2014-04 (G. Madec) add the diagnostic of the time filter trends
21	!! 3.7 ! 2015-11 (J. Chanut) Free surface simplification
22	!!-------------------------------------------------------------------------
23
24	!!-------------------------------------------------------------------------
25	!! dyn_nxt : obtain the next (after) horizontal velocity
26	!!-------------------------------------------------------------------------
27	USE oce ! ocean dynamics and tracers
28	USE dom_oce ! ocean space and time domain
29	USE sbc_oce ! Surface boundary condition: ocean fields
30	USE sbcrnf ! river runoffs
31	USE sbcisf ! ice shelf
32	USE phycst ! physical constants
33	USE dynadv ! dynamics: vector invariant versus flux form
34	USE dynspg_ts ! surface pressure gradient: split-explicit scheme
35	USE dynspg
36	USE domvvl ! variable volume
37	USE bdy_oce , ONLY : ln_bdy
38	USE bdydta ! ocean open boundary conditions
39	USE bdydyn ! ocean open boundary conditions
40	USE bdyvol ! ocean open boundary condition (bdy_vol routines)
41	USE trd_oce ! trends: ocean variables
42	USE trddyn ! trend manager: dynamics
43	USE trdken ! trend manager: kinetic energy
44	!
45	USE in_out_manager ! I/O manager
46	USE iom ! I/O manager library
47	USE lbclnk ! lateral boundary condition (or mpp link)
48	USE lib_mpp ! MPP library
49	USE prtctl ! Print control
50	USE timing ! Timing
51	USE zdfdrg , ONLY : ln_drgice_imp, rCdU_top
52	#if defined key_agrif
53	USE agrif_oce_interp
54	#endif
55
56	IMPLICIT NONE
57	PRIVATE
58
59	PUBLIC dyn_nxt ! routine called by step.F90
60
61	!! Substitution
62	# include "vectopt_loop_substitute.h90"
63	!!----------------------------------------------------------------------
64	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
65	!! $Id$
66	!! Software governed by the CeCILL license (see ./LICENSE)
67	!!----------------------------------------------------------------------
68	CONTAINS
69
70	SUBROUTINE dyn_nxt ( kt )
71	!!----------------------------------------------------------------------
72	!! * ROUTINE dyn_nxt *
73	!!
74	!! ** Purpose : Finalize after horizontal velocity. Apply the boundary
75	!! condition on the after velocity, achieve the time stepping
76	!! by applying the Asselin filter on now fields and swapping
77	!! the fields.
78	!!
79	!! ** Method : * Ensure after velocities transport matches time splitting
80	!! estimate (ln_dynspg_ts=T)
81	!!
82	!! * Apply lateral boundary conditions on after velocity
83	!! at the local domain boundaries through lbc_lnk call,
84	!! at the one-way open boundaries (ln_bdy=T),
85	!! at the AGRIF zoom boundaries (lk_agrif=T)
86	!!
87	!! * Apply the time filter applied and swap of the dynamics
88	!! arrays to start the next time step:
89	!! (ub,vb) = (un,vn) + atfp [ (ub,vb) + (ua,va) - 2 (un,vn) ]
90	!! (un,vn) = (ua,va).
91	!! Note that with flux form advection and non linear free surface,
92	!! the time filter is applied on thickness weighted velocity.
93	!! As a result, dyn_nxt MUST be called after tra_nxt.
94	!!
95	!! ** Action : ub,vb filtered before horizontal velocity of next time-step
96	!! un,vn now horizontal velocity of next time-step
97	!!----------------------------------------------------------------------
98	INTEGER, INTENT( in ) :: kt ! ocean time-step index
99	!
100	INTEGER :: ji, jj, jk ! dummy loop indices
101	INTEGER :: ikt ! local integers
102	REAL(wp) :: zue3a, zue3n, zue3b, zuf, zcoef ! local scalars
103	REAL(wp) :: zve3a, zve3n, zve3b, zvf, z1_2dt ! - -
104	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zue, zve
105	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zutau, zvtau
106	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ze3u_f, ze3v_f, zua, zva
107	!!----------------------------------------------------------------------
108	!
109	IF( ln_timing ) CALL timing_start('dyn_nxt')
110	IF( ln_dynspg_ts ) ALLOCATE( zue(jpi,jpj) , zve(jpi,jpj) )
111	IF( l_trddyn ) ALLOCATE( zua(jpi,jpj,jpk) , zva(jpi,jpj,jpk) )
112	!
113	IF( kt == nit000 ) THEN
114	IF(lwp) WRITE(numout,*)
115	IF(lwp) WRITE(numout,*) 'dyn_nxt : time stepping'
116	IF(lwp) WRITE(numout,*) '~~~~~~~'
117	ENDIF
118
119	IF ( ln_dynspg_ts ) THEN
120	! Ensure below that barotropic velocities match time splitting estimate
121	! Compute actual transport and replace it with ts estimate at "after" time step
122	zue(:,:) = e3u_a(:,:,1) * ua(:,:,1) * umask(:,:,1)
123	zve(:,:) = e3v_a(:,:,1) * va(:,:,1) * vmask(:,:,1)
124	DO jk = 2, jpkm1
125	zue(:,:) = zue(:,:) + e3u_a(:,:,jk) * ua(:,:,jk) * umask(:,:,jk)
126	zve(:,:) = zve(:,:) + e3v_a(:,:,jk) * va(:,:,jk) * vmask(:,:,jk)
127	END DO
128	DO jk = 1, jpkm1
129	ua(:,:,jk) = ( ua(:,:,jk) - zue(:,:) * r1_hu_a(:,:) + ua_b(:,:) ) * umask(:,:,jk)
130	va(:,:,jk) = ( va(:,:,jk) - zve(:,:) * r1_hv_a(:,:) + va_b(:,:) ) * vmask(:,:,jk)
131	END DO
132	!
133	IF( .NOT.ln_bt_fw ) THEN
134	! Remove advective velocity from "now velocities"
135	! prior to asselin filtering
136	! In the forward case, this is done below after asselin filtering
137	! so that asselin contribution is removed at the same time
138	DO jk = 1, jpkm1
139	un(:,:,jk) = ( un(:,:,jk) - un_adv(:,:)r1_hu_n(:,:) + un_b(:,:) )umask(:,:,jk)
140	vn(:,:,jk) = ( vn(:,:,jk) - vn_adv(:,:)r1_hv_n(:,:) + vn_b(:,:) )vmask(:,:,jk)
141	END DO
142	ENDIF
143	ENDIF
144
145	! Update after velocity on domain lateral boundaries
146	! --------------------------------------------------
147	# if defined key_agrif
148	CALL Agrif_dyn( kt ) !* AGRIF zoom boundaries
149	# endif
150	!
151	CALL lbc_lnk_multi( 'dynnxt', ua, 'U', -1., va, 'V', -1. ) !* local domain boundaries
152	!
153	! !* BDY open boundaries
154	IF( ln_bdy .AND. ln_dynspg_exp ) CALL bdy_dyn( kt )
155	IF( ln_bdy .AND. ln_dynspg_ts ) CALL bdy_dyn( kt, dyn3d_only=.true. )
156
157	!!$ Do we need a call to bdy_vol here??
158	!
159	IF( l_trddyn ) THEN ! prepare the atf trend computation + some diagnostics
160	z1_2dt = 1._wp / (2. * rdt) ! Euler or leap-frog time step
161	IF( neuler == 0 .AND. kt == nit000 ) z1_2dt = 1._wp / rdt
162	!
163	! ! Kinetic energy and Conversion
164	IF( ln_KE_trd ) CALL trd_dyn( ua, va, jpdyn_ken, kt )
165	!
166	IF( ln_dyn_trd ) THEN ! 3D output: total momentum trends
167	zua(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) * z1_2dt
168	zva(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) * z1_2dt
169	CALL iom_put( "utrd_tot", zua ) ! total momentum trends, except the asselin time filter
170	CALL iom_put( "vtrd_tot", zva )
171	ENDIF
172	!
173	zua(:,:,:) = un(:,:,:) ! save the now velocity before the asselin filter
174	zva(:,:,:) = vn(:,:,:) ! (caution: there will be a shift by 1 timestep in the
175	! ! computation of the asselin filter trends)
176	ENDIF
177
178	! Time filter and swap of dynamics arrays
179	! ------------------------------------------
180	IF( neuler == 0 .AND. kt == nit000 ) THEN !* Euler at first time-step: only swap
181	DO jk = 1, jpkm1
182	ub(:,:,jk) = un(:,:,jk) ! ub <-- un
183	vb(:,:,jk) = vn(:,:,jk)
184	un(:,:,jk) = ua(:,:,jk) ! un <-- ua
185	vn(:,:,jk) = va(:,:,jk)
186	END DO
187	! limit velocities
188	IF (ln_ulimit) THEN
189	call dyn_limit_velocity (kt)
190	ENDIF
191	IF( .NOT.ln_linssh ) THEN ! e3._b <-- e3._n
192	!!gm BUG ???? I don't understand why it is not : e3._n <-- e3._a
193	DO jk = 1, jpkm1
194	! e3t_b(:,:,jk) = e3t_n(:,:,jk)
195	! e3u_b(:,:,jk) = e3u_n(:,:,jk)
196	! e3v_b(:,:,jk) = e3v_n(:,:,jk)
197	!
198	e3t_n(:,:,jk) = e3t_a(:,:,jk)
199	e3u_n(:,:,jk) = e3u_a(:,:,jk)
200	e3v_n(:,:,jk) = e3v_a(:,:,jk)
201	END DO
202	!!gm BUG end
203	ENDIF
204	!
205
206	ELSE !* Leap-Frog : Asselin filter and swap
207	! ! =============!
208	IF( ln_linssh ) THEN ! Fixed volume !
209	! ! =============!
210	DO jk = 1, jpkm1
211	DO jj = 1, jpj
212	DO ji = 1, jpi
213	zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2._wp * un(ji,jj,jk) + ua(ji,jj,jk) )
214	zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2._wp * vn(ji,jj,jk) + va(ji,jj,jk) )
215	!
216	ub(ji,jj,jk) = zuf ! ub <-- filtered velocity
217	vb(ji,jj,jk) = zvf
218	un(ji,jj,jk) = ua(ji,jj,jk) ! un <-- ua
219	vn(ji,jj,jk) = va(ji,jj,jk)
220	END DO
221	END DO
222	END DO
223	! limit velocities
224	IF (ln_ulimit) THEN
225	call dyn_limit_velocity (kt)
226	ENDIF
227	! ! ================!
228	ELSE ! Variable volume !
229	! ! ================!
230	! Before scale factor at t-points
231	! (used as a now filtered scale factor until the swap)
232	! ----------------------------------------------------
233	DO jk = 1, jpkm1
234	e3t_b(:,:,jk) = e3t_n(:,:,jk) + atfp * ( e3t_b(:,:,jk) - 2._wp * e3t_n(:,:,jk) + e3t_a(:,:,jk) )
235	END DO
236	! Add volume filter correction: compatibility with tracer advection scheme
237	! => time filter + conservation correction (only at the first level)
238	zcoef = atfp * rdt * r1_rau0
239
240	DO jk = 1, jpkm1
241	e3t_b(:,:,jk) = e3t_b(:,:,jk) - zcoef * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,jk) &
242	& * e3t_n(:,:,jk) / ( ht_n(:,:) + 1._wp - ssmask(:,:) )
243	END DO
244
245	IF ( ln_rnf ) THEN
246	DO jk = 1, jpkm1
247	e3t_b(:,:,jk) = e3t_b(:,:,jk) + zcoef * ( rnf_b(:,:) - rnf(:,:) ) * tmask(:,:,jk) &
248	& * e3t_n(:,:,jk) / ( ht_n(:,:) + 1._wp - ssmask(:,:) )
249	END DO
250	ENDIF
251
252	IF ( ln_isf ) THEN
253	DO jk = 1, jpkm1
254	e3t_b(:,:,jk) = e3t_b(:,:,jk) - zcoef * ( fwfisf_b(:,:) - fwfisf(:,:) ) * tmask(:,:,jk) &
255	& * e3t_n(:,:,jk) / ( ht_n(:,:) + 1._wp - ssmask(:,:) )
256	END DO
257	ENDIF
258	!
259	IF( ln_dynadv_vec ) THEN ! Asselin filter applied on velocity
260	! Before filtered scale factor at (u/v)-points
261	CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' )
262	CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' )
263	DO jk = 1, jpkm1
264	DO jj = 1, jpj
265	DO ji = 1, jpi
266	zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2._wp * un(ji,jj,jk) + ua(ji,jj,jk) )
267	zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2._wp * vn(ji,jj,jk) + va(ji,jj,jk) )
268	!
269	ub(ji,jj,jk) = zuf ! ub <-- filtered velocity
270	vb(ji,jj,jk) = zvf
271	un(ji,jj,jk) = ua(ji,jj,jk) ! un <-- ua
272	vn(ji,jj,jk) = va(ji,jj,jk)
273	END DO
274	END DO
275	END DO
276	! limit velocities
277	IF (ln_ulimit) THEN
278	call dyn_limit_velocity (kt)
279	ENDIF
280	!
281	ELSE ! Asselin filter applied on thickness weighted velocity
282	!
283	ALLOCATE( ze3u_f(jpi,jpj,jpk) , ze3v_f(jpi,jpj,jpk) )
284	! Before filtered scale factor at (u/v)-points stored in ze3u_f, ze3v_f
285	CALL dom_vvl_interpol( e3t_b(:,:,:), ze3u_f, 'U' )
286	CALL dom_vvl_interpol( e3t_b(:,:,:), ze3v_f, 'V' )
287	DO jk = 1, jpkm1
288	DO jj = 1, jpj
289	DO ji = 1, jpi
290	zue3a = e3u_a(ji,jj,jk) * ua(ji,jj,jk)
291	zve3a = e3v_a(ji,jj,jk) * va(ji,jj,jk)
292	zue3n = e3u_n(ji,jj,jk) * un(ji,jj,jk)
293	zve3n = e3v_n(ji,jj,jk) * vn(ji,jj,jk)
294	zue3b = e3u_b(ji,jj,jk) * ub(ji,jj,jk)
295	zve3b = e3v_b(ji,jj,jk) * vb(ji,jj,jk)
296	!
297	zuf = ( zue3n + atfp * ( zue3b - 2._wp * zue3n + zue3a ) ) / ze3u_f(ji,jj,jk)
298	zvf = ( zve3n + atfp * ( zve3b - 2._wp * zve3n + zve3a ) ) / ze3v_f(ji,jj,jk)
299	!
300	ub(ji,jj,jk) = zuf ! ub <-- filtered velocity
301	vb(ji,jj,jk) = zvf
302	un(ji,jj,jk) = ua(ji,jj,jk) ! un <-- ua
303	vn(ji,jj,jk) = va(ji,jj,jk)
304	END DO
305	END DO
306	END DO
307	! limit velocities
308	IF (ln_ulimit) THEN
309	call dyn_limit_velocity (kt)
310	ENDIF
311	e3u_b(:,:,1:jpkm1) = ze3u_f(:,:,1:jpkm1) ! e3u_b <-- filtered scale factor
312	e3v_b(:,:,1:jpkm1) = ze3v_f(:,:,1:jpkm1)
313	!
314	DEALLOCATE( ze3u_f , ze3v_f )
315	ENDIF
316	!
317	ENDIF
318	!
319	IF( ln_dynspg_ts .AND. ln_bt_fw ) THEN
320	! Revert "before" velocities to time split estimate
321	! Doing it here also means that asselin filter contribution is removed
322	zue(:,:) = e3u_b(:,:,1) * ub(:,:,1) * umask(:,:,1)
323	zve(:,:) = e3v_b(:,:,1) * vb(:,:,1) * vmask(:,:,1)
324	DO jk = 2, jpkm1
325	zue(:,:) = zue(:,:) + e3u_b(:,:,jk) * ub(:,:,jk) * umask(:,:,jk)
326	zve(:,:) = zve(:,:) + e3v_b(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk)
327	END DO
328	DO jk = 1, jpkm1
329	ub(:,:,jk) = ub(:,:,jk) - (zue(:,:) * r1_hu_n(:,:) - un_b(:,:)) * umask(:,:,jk)
330	vb(:,:,jk) = vb(:,:,jk) - (zve(:,:) * r1_hv_n(:,:) - vn_b(:,:)) * vmask(:,:,jk)
331	END DO
332	ENDIF
333	!
334	ENDIF ! neuler =/0
335	!
336	! Set "now" and "before" barotropic velocities for next time step:
337	! JC: Would be more clever to swap variables than to make a full vertical
338	! integration
339	!
340	!
341	IF(.NOT.ln_linssh ) THEN
342	hu_b(:,:) = e3u_b(:,:,1) * umask(:,:,1)
343	hv_b(:,:) = e3v_b(:,:,1) * vmask(:,:,1)
344	DO jk = 2, jpkm1
345	hu_b(:,:) = hu_b(:,:) + e3u_b(:,:,jk) * umask(:,:,jk)
346	hv_b(:,:) = hv_b(:,:) + e3v_b(:,:,jk) * vmask(:,:,jk)
347	END DO
348	r1_hu_b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1._wp - ssumask(:,:) )
349	r1_hv_b(:,:) = ssvmask(:,:) / ( hv_b(:,:) + 1._wp - ssvmask(:,:) )
350	ENDIF
351	!
352	un_b(:,:) = e3u_a(:,:,1) * un(:,:,1) * umask(:,:,1)
353	ub_b(:,:) = e3u_b(:,:,1) * ub(:,:,1) * umask(:,:,1)
354	vn_b(:,:) = e3v_a(:,:,1) * vn(:,:,1) * vmask(:,:,1)
355	vb_b(:,:) = e3v_b(:,:,1) * vb(:,:,1) * vmask(:,:,1)
356	DO jk = 2, jpkm1
357	un_b(:,:) = un_b(:,:) + e3u_a(:,:,jk) * un(:,:,jk) * umask(:,:,jk)
358	ub_b(:,:) = ub_b(:,:) + e3u_b(:,:,jk) * ub(:,:,jk) * umask(:,:,jk)
359	vn_b(:,:) = vn_b(:,:) + e3v_a(:,:,jk) * vn(:,:,jk) * vmask(:,:,jk)
360	vb_b(:,:) = vb_b(:,:) + e3v_b(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk)
361	END DO
362	un_b(:,:) = un_b(:,:) * r1_hu_a(:,:)
363	vn_b(:,:) = vn_b(:,:) * r1_hv_a(:,:)
364	ub_b(:,:) = ub_b(:,:) * r1_hu_b(:,:)
365	vb_b(:,:) = vb_b(:,:) * r1_hv_b(:,:)
366	!
367	IF( .NOT.ln_dynspg_ts ) THEN ! output the barotropic currents
368	CALL iom_put( "ubar", un_b(:,:) )
369	CALL iom_put( "vbar", vn_b(:,:) )
370	ENDIF
371	IF( l_trddyn ) THEN ! 3D output: asselin filter trends on momentum
372	zua(:,:,:) = ( ub(:,:,:) - zua(:,:,:) ) * z1_2dt
373	zva(:,:,:) = ( vb(:,:,:) - zva(:,:,:) ) * z1_2dt
374	CALL trd_dyn( zua, zva, jpdyn_atf, kt )
375	ENDIF
376	!
377	IF ( iom_use("utau") ) THEN
378	IF ( ln_drgice_imp.OR.ln_isfcav ) THEN
379	ALLOCATE(zutau(jpi,jpj))
380	DO jj = 2, jpjm1
381	DO ji = 2, jpim1
382	jk = miku(ji,jj)
383	zutau(ji,jj) = utau(ji,jj) &
384	& + 0.5_wp * rau0 * (rCdU_top(ji+1,jj)+rCdU_top(ji,jj)) * ua(ji,jj,jk)
385	END DO
386	END DO
387	CALL lbc_lnk( 'dynnxt' , zutau, 'U', -1.)
388	CALL iom_put( "utau", zutau(:,:) )
389	DEALLOCATE(zutau)
390	ELSE
391	CALL iom_put( "utau", utau(:,:) )
392	ENDIF
393	ENDIF
394	!
395	IF ( iom_use("vtau") ) THEN
396	IF ( ln_drgice_imp.OR.ln_isfcav ) THEN
397	ALLOCATE(zvtau(jpi,jpj))
398	DO jj = 2, jpjm1
399	DO ji = 2, jpim1
400	jk = mikv(ji,jj)
401	zvtau(ji,jj) = vtau(ji,jj) &
402	& + 0.5_wp * rau0 * (rCdU_top(ji,jj+1)+rCdU_top(ji,jj)) * va(ji,jj,jk)
403	END DO
404	END DO
405	CALL lbc_lnk( 'dynnxt' , zvtau, 'V', -1.)
406	CALL iom_put( "vtau", zvtau(:,:) )
407	DEALLOCATE(zvtau)
408	ELSE
409	CALL iom_put( "vtau", vtau(:,:) )
410	ENDIF
411	ENDIF
412	!
413	IF(ln_ctl) CALL prt_ctl( tab3d_1=un, clinfo1=' nxt - Un: ', mask1=umask, &
414	& tab3d_2=vn, clinfo2=' Vn: ' , mask2=vmask )
415	!
416	IF( ln_dynspg_ts ) DEALLOCATE( zue, zve )
417	IF( l_trddyn ) DEALLOCATE( zua, zva )
418	IF( ln_timing ) CALL timing_stop('dyn_nxt')
419	!
420	END SUBROUTINE dyn_nxt
421
422	SUBROUTINE dyn_limit_velocity (kt)
423	! limits maximum values of un and vn by volume courant number
424	INTEGER, INTENT( in ) :: kt ! ocean time-step index
425	!
426	INTEGER :: ji, jj, jk ! dummy loop indices
427	REAL(wp) :: zzu,zplim,zmlim,isp,ism,zcn,ze3e1,zzcn,zcnn,idivp,idivm
428
429	! limit fluxes
430	zcn =cn_ulimit !0.9 ! maximum velocity inverse courant number
431	zcnn = cnn_ulimit !0.54 ! how much to reduce cn by in divergen flow
432
433	DO jk = 1, jpkm1
434	DO jj = 2, jpjm1
435	DO ji = fs_2, fs_jpim1 ! vect. opt.
436	! U direction
437	zzu = un(ji,jj,jk)
438	ze3e1 = e3u_n(ji ,jj,jk) * e2u(ji ,jj)
439	! ips is 1 if flow is positive othersize zero
440	isp = 0.5 * (sign(1.0_wp,zzu) + 1.0_wp )
441	ism = -0.5 * (sign(1.0_wp,zzu) - 1.0_wp )
442	! idev is 1 if divergent flow otherwise zero
443	idivp = -isp * 0.5 * (sign(1.0_wp, un(ji-1,jj,jk)) - 1.0_wp )
444	idivm = ism * 0.5 * (sign(1.0_wp, un(ji+1,jj,jk)) + 1.0_wp )
445	zzcn = (idivp+idivm)*(zcnn-1.0_wp)+1.0_wp
446	zzcn = zzcn * zcn
447	zplim = zzcn * (e3t_n(ji ,jj,jk) * e1t(ji ,jj) * e2t(ji ,jj)) / &
448	(2.0rdt ze3e1)*umask(ji,jj,jk)
449	zmlim = -zzcn * (e3t_n(ji+1,jj,jk) * e1t(ji+1,jj) * e2t(ji+1,jj)) / &
450	(2.0rdt ze3e1)*umask(ji,jj,jk)
451	! limit currents
452	un(ji,jj,jk) = min ( zzu,zplim) * isp + max(zzu,zmlim) *ism
453
454	! V direction
455	zzu = vn(ji,jj,jk)
456	ze3e1 = e3v_n(ji ,jj,jk) * e1v(ji ,jj)
457	isp = 0.5 * (sign(1.0_wp,zzu) + 1.0_wp )
458	ism = -0.5 * (sign(1.0_wp,zzu) - 1.0_wp )
459	! idev is 1 if divergent flow otherwise zero
460	idivp = -isp * 0.5 * (sign(1.0_wp, vn(ji,jj-1,jk)) - 1.0_wp )
461	idivm = ism * 0.5 * (sign(1.0_wp, vn(ji,jj+1,jk)) + 1.0_wp )
462	zzcn = (idivp+idivm)*(zcnn-1.0_wp)+1.0_wp
463	zzcn = zzcn * zcn
464	zplim = zzcn * (e3t_n(ji,jj ,jk) * e1t(ji,jj ) * e2t(ji,jj )) / &
465	(2.0rdt ze3e1)*vmask(ji,jj,jk)
466	zmlim = -zzcn * (e3t_n(ji,jj+1,jk) * e1t(ji,jj+1) * e2t(ji,jj+1)) / &
467	(2.0rdt ze3e1)*vmask(ji,jj,jk)
468	! limit currents
469	vn(ji,jj,jk) = min ( zzu,zplim) * isp + max(zzu,zmlim) *ism
470	ENDDO
471	ENDDO
472	ENDDO
473	CALL lbc_lnk_multi( 'dynnxt', un(:,:,:), 'U', -1., vn(:,:,:), 'V', -1. )
474
475	END SUBROUTINE dyn_limit_velocity
476
477	!!=========================================================================
478	END MODULE dynnxt

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

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

Context Navigation

source: NEMO/branches/UKMO/r14075_India_uncoupled/src/OCE/DYN/dynnxt.F90 @ 15422

Download in other formats: