Context Navigation

dynnxt.F90 @ 15176

Last change on this file since 15176 was 15176, checked in by davestorkey, 3 years ago

UKMO/NEMO_4.0.4_momentum_trends:

Output complete budget of 2D trends.
Separate out ice-ocean stress from wind stress (for 2D budget).
Allow thickness weighting of 2D trends.

File size: 19.5 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 domvvl ! variable volume
36	USE bdy_oce , ONLY : ln_bdy
37	USE bdydta ! ocean open boundary conditions
38	USE bdydyn ! ocean open boundary conditions
39	USE bdyvol ! ocean open boundary condition (bdy_vol routines)
40	USE trd_oce ! trends: ocean variables
41	USE trddyn ! trend manager: dynamics
42	USE trdken ! trend manager: kinetic energy
43	!
44	USE in_out_manager ! I/O manager
45	USE iom ! I/O manager library
46	USE lbclnk ! lateral boundary condition (or mpp link)
47	USE lib_mpp ! MPP library
48	USE prtctl ! Print control
49	USE timing ! Timing
50	USE zdfdrg , ONLY : ln_drgice_imp, rCdU_top
51	#if defined key_agrif
52	USE agrif_oce_interp
53	#endif
54
55	IMPLICIT NONE
56	PRIVATE
57
58	PUBLIC dyn_nxt ! routine called by step.F90
59
60	!!----------------------------------------------------------------------
61	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
62	!! $Id$
63	!! Software governed by the CeCILL license (see ./LICENSE)
64	!!----------------------------------------------------------------------
65	CONTAINS
66
67	SUBROUTINE dyn_nxt ( kt )
68	!!----------------------------------------------------------------------
69	!! * ROUTINE dyn_nxt *
70	!!
71	!! ** Purpose : Finalize after horizontal velocity. Apply the boundary
72	!! condition on the after velocity, achieve the time stepping
73	!! by applying the Asselin filter on now fields and swapping
74	!! the fields.
75	!!
76	!! ** Method : * Ensure after velocities transport matches time splitting
77	!! estimate (ln_dynspg_ts=T)
78	!!
79	!! * Apply lateral boundary conditions on after velocity
80	!! at the local domain boundaries through lbc_lnk call,
81	!! at the one-way open boundaries (ln_bdy=T),
82	!! at the AGRIF zoom boundaries (lk_agrif=T)
83	!!
84	!! * Apply the time filter applied and swap of the dynamics
85	!! arrays to start the next time step:
86	!! (ub,vb) = (un,vn) + atfp [ (ub,vb) + (ua,va) - 2 (un,vn) ]
87	!! (un,vn) = (ua,va).
88	!! Note that with flux form advection and non linear free surface,
89	!! the time filter is applied on thickness weighted velocity.
90	!! As a result, dyn_nxt MUST be called after tra_nxt.
91	!!
92	!! ** Action : ub,vb filtered before horizontal velocity of next time-step
93	!! un,vn now horizontal velocity of next time-step
94	!!----------------------------------------------------------------------
95	INTEGER, INTENT( in ) :: kt ! ocean time-step index
96	!
97	INTEGER :: ji, jj, jk ! dummy loop indices
98	INTEGER :: ikt ! local integers
99	REAL(wp) :: zue3a, zue3n, zue3b, zuf, zcoef ! local scalars
100	REAL(wp) :: zve3a, zve3n, zve3b, zvf, z1_2dt ! - -
101	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zue, zve
102	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zutau, zvtau
103	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ze3u_f, ze3v_f, zua, zva
104	!!----------------------------------------------------------------------
105	!
106	IF( ln_timing ) CALL timing_start('dyn_nxt')
107	IF( ln_dynspg_ts ) ALLOCATE( zue(jpi,jpj) , zve(jpi,jpj) )
108	IF( l_trddyn ) ALLOCATE( zua(jpi,jpj,jpk) , zva(jpi,jpj,jpk) )
109	!
110	IF( kt == nit000 ) THEN
111	IF(lwp) WRITE(numout,*)
112	IF(lwp) WRITE(numout,*) 'dyn_nxt : time stepping'
113	IF(lwp) WRITE(numout,*) '~~~~~~~'
114	ENDIF
115
116	IF ( ln_dynspg_ts ) THEN
117	! Ensure below that barotropic velocities match time splitting estimate
118	! Compute actual transport and replace it with ts estimate at "after" time step
119	zue(:,:) = e3u_a(:,:,1) * ua(:,:,1) * umask(:,:,1)
120	zve(:,:) = e3v_a(:,:,1) * va(:,:,1) * vmask(:,:,1)
121	DO jk = 2, jpkm1
122	zue(:,:) = zue(:,:) + e3u_a(:,:,jk) * ua(:,:,jk) * umask(:,:,jk)
123	zve(:,:) = zve(:,:) + e3v_a(:,:,jk) * va(:,:,jk) * vmask(:,:,jk)
124	END DO
125	DO jk = 1, jpkm1
126	ua(:,:,jk) = ( ua(:,:,jk) - zue(:,:) * r1_hu_a(:,:) + ua_b(:,:) ) * umask(:,:,jk)
127	va(:,:,jk) = ( va(:,:,jk) - zve(:,:) * r1_hv_a(:,:) + va_b(:,:) ) * vmask(:,:,jk)
128	END DO
129	!
130	IF( .NOT.ln_bt_fw ) THEN
131	! Remove advective velocity from "now velocities"
132	! prior to asselin filtering
133	! In the forward case, this is done below after asselin filtering
134	! so that asselin contribution is removed at the same time
135	DO jk = 1, jpkm1
136	un(:,:,jk) = ( un(:,:,jk) - un_adv(:,:)r1_hu_n(:,:) + un_b(:,:) )umask(:,:,jk)
137	vn(:,:,jk) = ( vn(:,:,jk) - vn_adv(:,:)r1_hv_n(:,:) + vn_b(:,:) )vmask(:,:,jk)
138	END DO
139	ENDIF
140	ENDIF
141
142	! Update after velocity on domain lateral boundaries
143	! --------------------------------------------------
144	# if defined key_agrif
145	CALL Agrif_dyn( kt ) !* AGRIF zoom boundaries
146	# endif
147	!
148	CALL lbc_lnk_multi( 'dynnxt', ua, 'U', -1., va, 'V', -1. ) !* local domain boundaries
149	!
150	! !* BDY open boundaries
151	IF( ln_bdy .AND. ln_dynspg_exp ) CALL bdy_dyn( kt )
152	IF( ln_bdy .AND. ln_dynspg_ts ) CALL bdy_dyn( kt, dyn3d_only=.true. )
153
154	!!$ Do we need a call to bdy_vol here??
155	!
156	IF( l_trddyn ) THEN ! prepare the atf trend computation + some diagnostics
157	z1_2dt = 1._wp / (2. * rdt) ! Euler or leap-frog time step
158	IF( neuler == 0 .AND. kt == nit000 ) z1_2dt = 1._wp / rdt
159	!
160	! ! Kinetic energy and Conversion
161	IF( ln_KE_trd ) CALL trd_dyn( ua, va, jpdyn_ken, kt )
162	!
163	IF( ln_dyn_trd ) THEN ! 3D output: total momentum trends
164	IF( ln_linssh ) THEN
165	zua(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) * z1_2dt
166	zva(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) * z1_2dt
167	ELSE
168	zua(:,:,:) = ( e3u_a(:,:,:)ua(:,:,:) - e3u_b(:,:,:)ub(:,:,:) ) * z1_2dt / e3u_n(:,:,:)
169	zva(:,:,:) = ( e3v_a(:,:,:)va(:,:,:) - e3v_b(:,:,:)vb(:,:,:) ) * z1_2dt / e3v_n(:,:,:)
170	ENDIF
171	CALL trd_dyn( zua, zva, jpdyn_tot, kt ) ! total momentum trends, except the asselin time filter
172	ENDIF
173	!
174	zua(:,:,:) = un(:,:,:) ! save the now velocity before the asselin filter
175	zva(:,:,:) = vn(:,:,:) ! (caution: there will be a shift by 1 timestep in the
176	! ! computation of the asselin filter trends)
177	ENDIF
178
179	! Time filter and swap of dynamics arrays
180	! ------------------------------------------
181	IF( neuler == 0 .AND. kt == nit000 ) THEN !* Euler at first time-step: only swap
182	DO jk = 1, jpkm1
183	ub(:,:,jk) = un(:,:,jk) ! ub <-- un
184	vb(:,:,jk) = vn(:,:,jk)
185	un(:,:,jk) = ua(:,:,jk) ! un <-- ua
186	vn(:,:,jk) = va(:,:,jk)
187	END DO
188	IF( .NOT.ln_linssh ) THEN ! e3._b <-- e3._n
189	!!gm BUG ???? I don't understand why it is not : e3._n <-- e3._a
190	DO jk = 1, jpkm1
191	! e3t_b(:,:,jk) = e3t_n(:,:,jk)
192	! e3u_b(:,:,jk) = e3u_n(:,:,jk)
193	! e3v_b(:,:,jk) = e3v_n(:,:,jk)
194	!
195	e3t_n(:,:,jk) = e3t_a(:,:,jk)
196	e3u_n(:,:,jk) = e3u_a(:,:,jk)
197	e3v_n(:,:,jk) = e3v_a(:,:,jk)
198	END DO
199	!!gm BUG end
200	ENDIF
201	!
202
203	ELSE !* Leap-Frog : Asselin filter and swap
204	! ! =============!
205	IF( ln_linssh ) THEN ! Fixed volume !
206	! ! =============!
207	DO jk = 1, jpkm1
208	DO jj = 1, jpj
209	DO ji = 1, jpi
210	zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2._wp * un(ji,jj,jk) + ua(ji,jj,jk) )
211	zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2._wp * vn(ji,jj,jk) + va(ji,jj,jk) )
212	!
213	ub(ji,jj,jk) = zuf ! ub <-- filtered velocity
214	vb(ji,jj,jk) = zvf
215	un(ji,jj,jk) = ua(ji,jj,jk) ! un <-- ua
216	vn(ji,jj,jk) = va(ji,jj,jk)
217	END DO
218	END DO
219	END DO
220	! ! ================!
221	ELSE ! Variable volume !
222	! ! ================!
223	! Before scale factor at t-points
224	! (used as a now filtered scale factor until the swap)
225	! ----------------------------------------------------
226	DO jk = 1, jpkm1
227	e3t_b(:,:,jk) = e3t_n(:,:,jk) + atfp * ( e3t_b(:,:,jk) - 2._wp * e3t_n(:,:,jk) + e3t_a(:,:,jk) )
228	END DO
229	! Add volume filter correction: compatibility with tracer advection scheme
230	! => time filter + conservation correction (only at the first level)
231	zcoef = atfp * rdt * r1_rau0
232
233	DO jk = 1, jpkm1
234	e3t_b(:,:,jk) = e3t_b(:,:,jk) - zcoef * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,jk) &
235	& * e3t_n(:,:,jk) / ( ht_n(:,:) + 1._wp - ssmask(:,:) )
236	END DO
237
238	IF ( ln_rnf ) THEN
239	DO jk = 1, jpkm1
240	e3t_b(:,:,jk) = e3t_b(:,:,jk) + zcoef * ( rnf_b(:,:) - rnf(:,:) ) * tmask(:,:,jk) &
241	& * e3t_n(:,:,jk) / ( ht_n(:,:) + 1._wp - ssmask(:,:) )
242	END DO
243	ENDIF
244
245	IF ( ln_isf ) THEN
246	DO jk = 1, jpkm1
247	e3t_b(:,:,jk) = e3t_b(:,:,jk) - zcoef * ( fwfisf_b(:,:) - fwfisf(:,:) ) * tmask(:,:,jk) &
248	& * e3t_n(:,:,jk) / ( ht_n(:,:) + 1._wp - ssmask(:,:) )
249	END DO
250	ENDIF
251	!
252	IF( ln_dynadv_vec ) THEN ! Asselin filter applied on velocity
253	! Before filtered scale factor at (u/v)-points
254	CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' )
255	CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' )
256	DO jk = 1, jpkm1
257	DO jj = 1, jpj
258	DO ji = 1, jpi
259	zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2._wp * un(ji,jj,jk) + ua(ji,jj,jk) )
260	zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2._wp * vn(ji,jj,jk) + va(ji,jj,jk) )
261	!
262	ub(ji,jj,jk) = zuf ! ub <-- filtered velocity
263	vb(ji,jj,jk) = zvf
264	un(ji,jj,jk) = ua(ji,jj,jk) ! un <-- ua
265	vn(ji,jj,jk) = va(ji,jj,jk)
266	END DO
267	END DO
268	END DO
269	!
270	ELSE ! Asselin filter applied on thickness weighted velocity
271	!
272	ALLOCATE( ze3u_f(jpi,jpj,jpk) , ze3v_f(jpi,jpj,jpk) )
273	! Before filtered scale factor at (u/v)-points stored in ze3u_f, ze3v_f
274	CALL dom_vvl_interpol( e3t_b(:,:,:), ze3u_f, 'U' )
275	CALL dom_vvl_interpol( e3t_b(:,:,:), ze3v_f, 'V' )
276	DO jk = 1, jpkm1
277	DO jj = 1, jpj
278	DO ji = 1, jpi
279	zue3a = e3u_a(ji,jj,jk) * ua(ji,jj,jk)
280	zve3a = e3v_a(ji,jj,jk) * va(ji,jj,jk)
281	zue3n = e3u_n(ji,jj,jk) * un(ji,jj,jk)
282	zve3n = e3v_n(ji,jj,jk) * vn(ji,jj,jk)
283	zue3b = e3u_b(ji,jj,jk) * ub(ji,jj,jk)
284	zve3b = e3v_b(ji,jj,jk) * vb(ji,jj,jk)
285	!
286	zuf = ( zue3n + atfp * ( zue3b - 2._wp * zue3n + zue3a ) ) / ze3u_f(ji,jj,jk)
287	zvf = ( zve3n + atfp * ( zve3b - 2._wp * zve3n + zve3a ) ) / ze3v_f(ji,jj,jk)
288	!
289	ub(ji,jj,jk) = zuf ! ub <-- filtered velocity
290	vb(ji,jj,jk) = zvf
291	un(ji,jj,jk) = ua(ji,jj,jk) ! un <-- ua
292	vn(ji,jj,jk) = va(ji,jj,jk)
293	END DO
294	END DO
295	END DO
296	e3u_b(:,:,1:jpkm1) = ze3u_f(:,:,1:jpkm1) ! e3u_b <-- filtered scale factor
297	e3v_b(:,:,1:jpkm1) = ze3v_f(:,:,1:jpkm1)
298	!
299	DEALLOCATE( ze3u_f , ze3v_f )
300	ENDIF
301	!
302	ENDIF
303	!
304	IF( ln_dynspg_ts .AND. ln_bt_fw ) THEN
305	! Revert "before" velocities to time split estimate
306	! Doing it here also means that asselin filter contribution is removed
307	zue(:,:) = e3u_b(:,:,1) * ub(:,:,1) * umask(:,:,1)
308	zve(:,:) = e3v_b(:,:,1) * vb(:,:,1) * vmask(:,:,1)
309	DO jk = 2, jpkm1
310	zue(:,:) = zue(:,:) + e3u_b(:,:,jk) * ub(:,:,jk) * umask(:,:,jk)
311	zve(:,:) = zve(:,:) + e3v_b(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk)
312	END DO
313	DO jk = 1, jpkm1
314	ub(:,:,jk) = ub(:,:,jk) - (zue(:,:) * r1_hu_n(:,:) - un_b(:,:)) * umask(:,:,jk)
315	vb(:,:,jk) = vb(:,:,jk) - (zve(:,:) * r1_hv_n(:,:) - vn_b(:,:)) * vmask(:,:,jk)
316	END DO
317	ENDIF
318	!
319	ENDIF ! neuler =/0
320	!
321	! Set "now" and "before" barotropic velocities for next time step:
322	! JC: Would be more clever to swap variables than to make a full vertical
323	! integration
324	!
325	!
326	IF(.NOT.ln_linssh ) THEN
327	hu_b(:,:) = e3u_b(:,:,1) * umask(:,:,1)
328	hv_b(:,:) = e3v_b(:,:,1) * vmask(:,:,1)
329	DO jk = 2, jpkm1
330	hu_b(:,:) = hu_b(:,:) + e3u_b(:,:,jk) * umask(:,:,jk)
331	hv_b(:,:) = hv_b(:,:) + e3v_b(:,:,jk) * vmask(:,:,jk)
332	END DO
333	r1_hu_b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1._wp - ssumask(:,:) )
334	r1_hv_b(:,:) = ssvmask(:,:) / ( hv_b(:,:) + 1._wp - ssvmask(:,:) )
335	ENDIF
336	!
337	un_b(:,:) = e3u_a(:,:,1) * un(:,:,1) * umask(:,:,1)
338	ub_b(:,:) = e3u_b(:,:,1) * ub(:,:,1) * umask(:,:,1)
339	vn_b(:,:) = e3v_a(:,:,1) * vn(:,:,1) * vmask(:,:,1)
340	vb_b(:,:) = e3v_b(:,:,1) * vb(:,:,1) * vmask(:,:,1)
341	DO jk = 2, jpkm1
342	un_b(:,:) = un_b(:,:) + e3u_a(:,:,jk) * un(:,:,jk) * umask(:,:,jk)
343	ub_b(:,:) = ub_b(:,:) + e3u_b(:,:,jk) * ub(:,:,jk) * umask(:,:,jk)
344	vn_b(:,:) = vn_b(:,:) + e3v_a(:,:,jk) * vn(:,:,jk) * vmask(:,:,jk)
345	vb_b(:,:) = vb_b(:,:) + e3v_b(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk)
346	END DO
347	un_b(:,:) = un_b(:,:) * r1_hu_a(:,:)
348	vn_b(:,:) = vn_b(:,:) * r1_hv_a(:,:)
349	ub_b(:,:) = ub_b(:,:) * r1_hu_b(:,:)
350	vb_b(:,:) = vb_b(:,:) * r1_hv_b(:,:)
351	!
352	IF( .NOT.ln_dynspg_ts ) THEN ! output the barotropic currents
353	CALL iom_put( "ubar", un_b(:,:) )
354	CALL iom_put( "vbar", vn_b(:,:) )
355	ENDIF
356	IF( l_trddyn ) THEN ! 3D output: asselin filter trends on momentum
357	zua(:,:,:) = ( ub(:,:,:) - zua(:,:,:) ) * z1_2dt
358	zva(:,:,:) = ( vb(:,:,:) - zva(:,:,:) ) * z1_2dt
359	CALL trd_dyn( zua, zva, jpdyn_atf, kt )
360	ENDIF
361	!
362	IF ( iom_use("utau") ) THEN
363	IF ( ln_drgice_imp.OR.ln_isfcav ) THEN
364	ALLOCATE(zutau(jpi,jpj))
365	DO jj = 2, jpjm1
366	DO ji = 2, jpim1
367	jk = miku(ji,jj)
368	zutau(ji,jj) = utau(ji,jj) &
369	& + 0.5_wp * rau0 * (rCdU_top(ji+1,jj)+rCdU_top(ji,jj)) * ua(ji,jj,jk)
370	END DO
371	END DO
372	CALL lbc_lnk( 'dynnxt' , zutau, 'U', -1.)
373	CALL iom_put( "utau", zutau(:,:) )
374	DEALLOCATE(zutau)
375	ELSE
376	CALL iom_put( "utau", utau(:,:) )
377	ENDIF
378	ENDIF
379	!
380	IF ( iom_use("vtau") ) THEN
381	IF ( ln_drgice_imp.OR.ln_isfcav ) THEN
382	ALLOCATE(zvtau(jpi,jpj))
383	DO jj = 2, jpjm1
384	DO ji = 2, jpim1
385	jk = mikv(ji,jj)
386	zvtau(ji,jj) = vtau(ji,jj) &
387	& + 0.5_wp * rau0 * (rCdU_top(ji,jj+1)+rCdU_top(ji,jj)) * va(ji,jj,jk)
388	END DO
389	END DO
390	CALL lbc_lnk( 'dynnxt' , zvtau, 'V', -1.)
391	CALL iom_put( "vtau", zvtau(:,:) )
392	DEALLOCATE(zvtau)
393	ELSE
394	CALL iom_put( "vtau", vtau(:,:) )
395	ENDIF
396	ENDIF
397	!
398	IF(ln_ctl) CALL prt_ctl( tab3d_1=un, clinfo1=' nxt - Un: ', mask1=umask, &
399	& tab3d_2=vn, clinfo2=' Vn: ' , mask2=vmask )
400	!
401	IF( ln_dynspg_ts ) DEALLOCATE( zue, zve )
402	IF( l_trddyn ) DEALLOCATE( zua, zva )
403	IF( ln_timing ) CALL timing_stop('dyn_nxt')
404	!
405	END SUBROUTINE dyn_nxt
406
407	!!=========================================================================
408	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/NEMO_4.0.4_momentum_trends/src/OCE/DYN/dynnxt.F90 @ 15176

Download in other formats: