Context Navigation

traadv_tvd.F90 @ 7771

Last change on this file since 7771 was 7771, checked in by frrh, 7 years ago

Apply optimisations to various areas of code replacing the use of
allocated pointers with straightforward direct ALLOCATE and DEALLOCATE
operations.

These optimisations largely have an impact in models featuring MEDUSA,
i.e. those with significant numbers of tracers, although they are
expected to have a small impact in all configurations.

Code developed and tested in NEMO branch branches/UKMO/dev_r5518_optim_GO6_alloc
Tested in stand-alone GO6-GSI8, GO6-GSI8-MEDUSA and UKESM coupled models.
NEMO ticket #1821 documents this change further.

File size: 35.0 KB

Line
1	MODULE traadv_tvd
2	!!==============================================================================
3	!! * MODULE traadv_tvd *
4	!! Ocean tracers: horizontal & vertical advective trend
5	!!==============================================================================
6	!! History : OPA ! 1995-12 (L. Mortier) Original code
7	!! ! 2000-01 (H. Loukos) adapted to ORCA
8	!! ! 2000-10 (MA Foujols E.Kestenare) include file not routine
9	!! ! 2000-12 (E. Kestenare M. Levy) fix bug in trtrd indexes
10	!! ! 2001-07 (E. Durand G. Madec) adaptation to ORCA config
11	!! 8.5 ! 2002-06 (G. Madec) F90: Free form and module
12	!! NEMO 1.0 ! 2004-01 (A. de Miranda, G. Madec, J.M. Molines ): advective bbl
13	!! 2.0 ! 2008-04 (S. Cravatte) add the i-, j- & k- trends computation
14	!! - ! 2009-11 (V. Garnier) Surface pressure gradient organization
15	!! 3.3 ! 2010-05 (C. Ethe, G. Madec) merge TRC-TRA + switch from velocity to transport
16	!!----------------------------------------------------------------------
17
18	!!----------------------------------------------------------------------
19	!! tra_adv_tvd : update the tracer trend with the 3D advection trends using a TVD scheme
20	!! nonosc : compute monotonic tracer fluxes by a non-oscillatory algorithm
21	!!----------------------------------------------------------------------
22	USE oce ! ocean dynamics and active tracers
23	USE dom_oce ! ocean space and time domain
24	USE trc_oce ! share passive tracers/Ocean variables
25	USE trd_oce ! trends: ocean variables
26	USE trdtra ! tracers trends
27	USE dynspg_oce ! choice/control of key cpp for surface pressure gradient
28	USE diaptr ! poleward transport diagnostics
29	!
30	USE lib_mpp ! MPP library
31	USE lbclnk ! ocean lateral boundary condition (or mpp link)
32	USE in_out_manager ! I/O manager
33	USE wrk_nemo ! Memory Allocation
34	USE timing ! Timing
35	USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
36	USE iom
37
38	IMPLICIT NONE
39	PRIVATE
40
41	PUBLIC tra_adv_tvd ! routine called by traadv.F90
42	PUBLIC tra_adv_tvd_zts ! routine called by traadv.F90
43
44	LOGICAL :: l_trd ! flag to compute trends
45	LOGICAL :: l_trans ! flag to output vertically integrated transports
46
47	!! * Substitutions
48	# include "domzgr_substitute.h90"
49	# include "vectopt_loop_substitute.h90"
50	!!----------------------------------------------------------------------
51	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
52	!! $Id$
53	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
54	!!----------------------------------------------------------------------
55	CONTAINS
56
57	SUBROUTINE tra_adv_tvd ( kt, kit000, cdtype, p2dt, pun, pvn, pwn, &
58	& ptb, ptn, pta, kjpt )
59	!!----------------------------------------------------------------------
60	!! * ROUTINE tra_adv_tvd *
61	!!
62	!! ** Purpose : Compute the now trend due to total advection of
63	!! tracers and add it to the general trend of tracer equations
64	!!
65	!! ** Method : TVD scheme, i.e. 2nd order centered scheme with
66	!! corrected flux (monotonic correction)
67	!! note: - this advection scheme needs a leap-frog time scheme
68	!!
69	!! ** Action : - update (pta) with the now advective tracer trends
70	!! - save the trends
71	!!----------------------------------------------------------------------
72	USE oce , ONLY: zwx => ua , zwy => va ! (ua,va) used as workspace
73	!
74	INTEGER , INTENT(in ) :: kt ! ocean time-step index
75	INTEGER , INTENT(in ) :: kit000 ! first time step index
76	CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator)
77	INTEGER , INTENT(in ) :: kjpt ! number of tracers
78	REAL(wp), DIMENSION( jpk ), INTENT(in ) :: p2dt ! vertical profile of tracer time-step
79	REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pun, pvn, pwn ! 3 ocean velocity components
80	REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb, ptn ! before and now tracer fields
81	REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend
82	!
83	INTEGER :: ji, jj, jk, jn ! dummy loop indices
84	INTEGER :: ik
85	REAL(wp) :: z2dtt, zbtr, ztra ! local scalar
86	REAL(wp) :: zfp_ui, zfp_vj, zfp_wk ! - -
87	REAL(wp) :: zfm_ui, zfm_vj, zfm_wk ! - -
88	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zwi, zwz
89	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrdx, ztrdy, ztrdz, zptry
90	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z2d
91	!!----------------------------------------------------------------------
92	!
93	IF( nn_timing == 1 ) CALL timing_start('tra_adv_tvd')
94	!
95	ALLOCATE(zwi(1:jpi, 1:jpj, 1:jpk))
96	ALLOCATE(zwz(1:jpi, 1:jpj, 1:jpk))
97
98	!
99	IF( kt == kit000 ) THEN
100	IF(lwp) WRITE(numout,*)
101	IF(lwp) WRITE(numout,*) 'tra_adv_tvd : TVD advection scheme on ', cdtype
102	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
103	!
104	ENDIF
105	l_trd = .FALSE.
106	l_trans = .FALSE.
107	IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE.
108	IF( cdtype == 'TRA' .AND. (iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") ) ) l_trans = .TRUE.
109	!
110	IF( l_trd .OR. l_trans ) THEN
111	ALLOCATE(ztrdx(1:jpi, 1:jpj, 1:jpk))
112	ALLOCATE(ztrdy(1:jpi, 1:jpj, 1:jpk))
113	ALLOCATE(ztrdz(1:jpi, 1:jpj, 1:jpk))
114	ztrdx(:,:,:) = 0.e0 ; ztrdy(:,:,:) = 0.e0 ; ztrdz(:,:,:) = 0.e0
115	ALLOCATE(z2d(1:jpi, 1:jpj))
116	ENDIF
117	!
118	IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN
119	ALLOCATE(zptry(1:jpi, 1:jpj, 1:jpk))
120	zptry(:,:,:) = 0._wp
121	ENDIF
122	!
123	zwi(:,:,:) = 0.e0 ;
124	!
125	! ! ===========
126	DO jn = 1, kjpt ! tracer loop
127	! ! ===========
128	! 1. Bottom and k=1 value : flux set to zero
129	! ----------------------------------
130	zwx(:,:,jpk) = 0.e0 ; zwz(:,:,jpk) = 0.e0
131	zwy(:,:,jpk) = 0.e0 ; zwi(:,:,jpk) = 0.e0
132
133	zwz(:,:,1 ) = 0._wp
134	! 2. upstream advection with initial mass fluxes & intermediate update
135	! --------------------------------------------------------------------
136	! upstream tracer flux in the i and j direction
137	DO jk = 1, jpkm1
138	DO jj = 1, jpjm1
139	DO ji = 1, fs_jpim1 ! vector opt.
140	! upstream scheme
141	zfp_ui = pun(ji,jj,jk) + ABS( pun(ji,jj,jk) )
142	zfm_ui = pun(ji,jj,jk) - ABS( pun(ji,jj,jk) )
143	zfp_vj = pvn(ji,jj,jk) + ABS( pvn(ji,jj,jk) )
144	zfm_vj = pvn(ji,jj,jk) - ABS( pvn(ji,jj,jk) )
145	zwx(ji,jj,jk) = 0.5 * ( zfp_ui * ptb(ji,jj,jk,jn) + zfm_ui * ptb(ji+1,jj ,jk,jn) )
146	zwy(ji,jj,jk) = 0.5 * ( zfp_vj * ptb(ji,jj,jk,jn) + zfm_vj * ptb(ji ,jj+1,jk,jn) )
147	END DO
148	END DO
149	END DO
150
151	! upstream tracer flux in the k direction
152	! Interior value
153	DO jk = 2, jpkm1
154	DO jj = 1, jpj
155	DO ji = 1, jpi
156	zfp_wk = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) )
157	zfm_wk = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) )
158	zwz(ji,jj,jk) = 0.5 * ( zfp_wk * ptb(ji,jj,jk,jn) + zfm_wk * ptb(ji,jj,jk-1,jn) ) * wmask(ji,jj,jk)
159	END DO
160	END DO
161	END DO
162	! Surface value
163	IF( lk_vvl ) THEN
164	IF ( ln_isfcav ) THEN
165	DO jj = 1, jpj
166	DO ji = 1, jpi
167	zwz(ji,jj, mikt(ji,jj) ) = 0.e0 ! volume variable
168	END DO
169	END DO
170	ELSE
171	zwz(:,:,1) = 0.e0 ! volume variable
172	END IF
173	ELSE
174	IF ( ln_isfcav ) THEN
175	DO jj = 1, jpj
176	DO ji = 1, jpi
177	zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptb(ji,jj,mikt(ji,jj),jn) ! linear free surface
178	END DO
179	END DO
180	ELSE
181	zwz(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn) ! linear free surface
182	END IF
183	ENDIF
184
185	! total advective trend
186	DO jk = 1, jpkm1
187	z2dtt = p2dt(jk)
188	DO jj = 2, jpjm1
189	DO ji = fs_2, fs_jpim1 ! vector opt.
190	! total intermediate advective trends
191	ztra = - ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) &
192	& + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) &
193	& + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) / e1e2t(ji,jj)
194	! update and guess with monotonic sheme
195	pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra / fse3t_n(ji,jj,jk) * tmask(ji,jj,jk)
196	zwi(ji,jj,jk) = ( fse3t_b(ji,jj,jk) * ptb(ji,jj,jk,jn) + z2dtt * ztra ) / fse3t_a(ji,jj,jk) * tmask(ji,jj,jk)
197	END DO
198	END DO
199	END DO
200	! ! Lateral boundary conditions on zwi (unchanged sign)
201	CALL lbc_lnk( zwi, 'T', 1. )
202
203	! ! trend diagnostics (contribution of upstream fluxes)
204	IF( l_trd .OR. l_trans ) THEN
205	! store intermediate advective trends
206	ztrdx(:,:,:) = zwx(:,:,:) ; ztrdy(:,:,:) = zwy(:,:,:) ; ztrdz(:,:,:) = zwz(:,:,:)
207	END IF
208	! ! "Poleward" heat and salt transports (contribution of upstream fluxes)
209	IF( cdtype == 'TRA' .AND. ln_diaptr ) zptry(:,:,:) = zwy(:,:,:)
210
211	! 3. antidiffusive flux : high order minus low order
212	! --------------------------------------------------
213	! antidiffusive flux on i and j
214	DO jk = 1, jpkm1
215	DO jj = 1, jpjm1
216	DO ji = 1, fs_jpim1 ! vector opt.
217	zwx(ji,jj,jk) = 0.5 * pun(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji+1,jj,jk,jn) ) - zwx(ji,jj,jk)
218	zwy(ji,jj,jk) = 0.5 * pvn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj+1,jk,jn) ) - zwy(ji,jj,jk)
219	END DO
220	END DO
221	END DO
222
223	! antidiffusive flux on k
224	! Interior value
225	DO jk = 2, jpkm1
226	DO jj = 1, jpj
227	DO ji = 1, jpi
228	zwz(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj,jk-1,jn) ) - zwz(ji,jj,jk)
229	END DO
230	END DO
231	END DO
232	! surface value
233	IF ( ln_isfcav ) THEN
234	DO jj = 1, jpj
235	DO ji = 1, jpi
236	zwz(ji,jj,mikt(ji,jj)) = 0.e0
237	END DO
238	END DO
239	ELSE
240	zwz(:,:,1) = 0.e0
241	END IF
242	CALL lbc_lnk( zwx, 'U', -1. ) ; CALL lbc_lnk( zwy, 'V', -1. ) ! Lateral bondary conditions
243	CALL lbc_lnk( zwz, 'W', 1. )
244
245	! 4. monotonicity algorithm
246	! -------------------------
247	CALL nonosc( ptb(:,:,:,jn), zwx, zwy, zwz, zwi, p2dt )
248
249
250	! 5. final trend with corrected fluxes
251	! ------------------------------------
252	DO jk = 1, jpkm1
253	DO jj = 2, jpjm1
254	DO ji = fs_2, fs_jpim1 ! vector opt.
255	zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
256	! total advective trends
257	ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) &
258	& + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) &
259	& + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) )
260	! add them to the general tracer trends
261	pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra * tmask(ji,jj,jk)
262	END DO
263	END DO
264	END DO
265
266	! ! trend diagnostics (contribution of upstream fluxes)
267	IF( l_trd .OR. l_trans ) THEN
268	ztrdx(:,:,:) = ztrdx(:,:,:) + zwx(:,:,:) ! <<< Add to previously computed
269	ztrdy(:,:,:) = ztrdy(:,:,:) + zwy(:,:,:) ! <<< Add to previously computed
270	ztrdz(:,:,:) = ztrdz(:,:,:) + zwz(:,:,:) ! <<< Add to previously computed
271	ENDIF
272
273	IF( l_trd ) THEN
274	CALL trd_tra( kt, cdtype, jn, jptra_xad, ztrdx, pun, ptn(:,:,:,jn) )
275	CALL trd_tra( kt, cdtype, jn, jptra_yad, ztrdy, pvn, ptn(:,:,:,jn) )
276	CALL trd_tra( kt, cdtype, jn, jptra_zad, ztrdz, pwn, ptn(:,:,:,jn) )
277	END IF
278
279	IF( l_trans .AND. jn==jp_tem ) THEN
280	z2d(:,:) = 0._wp
281	DO jk = 1, jpkm1
282	DO jj = 2, jpjm1
283	DO ji = fs_2, fs_jpim1 ! vector opt.
284	z2d(ji,jj) = z2d(ji,jj) + ztrdx(ji,jj,jk)
285	END DO
286	END DO
287	END DO
288	CALL lbc_lnk( z2d, 'U', -1. )
289	CALL iom_put( "uadv_heattr", rau0_rcp * z2d ) ! heat transport in i-direction
290	!
291	z2d(:,:) = 0._wp
292	DO jk = 1, jpkm1
293	DO jj = 2, jpjm1
294	DO ji = fs_2, fs_jpim1 ! vector opt.
295	z2d(ji,jj) = z2d(ji,jj) + ztrdy(ji,jj,jk)
296	END DO
297	END DO
298	END DO
299	CALL lbc_lnk( z2d, 'V', -1. )
300	CALL iom_put( "vadv_heattr", rau0_rcp * z2d ) ! heat transport in j-direction
301	ENDIF
302	! "Poleward" heat and salt transports (contribution of upstream fluxes)
303	IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN
304	zptry(:,:,:) = zptry(:,:,:) + zwy(:,:,:) ! <<< Add to previously computed
305	CALL dia_ptr_ohst_components( jn, 'adv', zptry(:,:,:) )
306	ENDIF
307	!
308	END DO
309	!
310	DEALLOCATE( zwi )
311	DEALLOCATE( zwz )
312	IF( l_trd .OR. l_trans ) THEN
313	DEALLOCATE( ztrdx )
314	DEALLOCATE( ztrdy )
315	DEALLOCATE( ztrdz )
316	DEALLOCATE( z2d )
317	ENDIF
318	IF( cdtype == 'TRA' .AND. ln_diaptr ) DEALLOCATE( zptry )
319	!
320	IF( nn_timing == 1 ) CALL timing_stop('tra_adv_tvd')
321	!
322	END SUBROUTINE tra_adv_tvd
323
324	SUBROUTINE tra_adv_tvd_zts ( kt, kit000, cdtype, p2dt, pun, pvn, pwn, &
325	& ptb, ptn, pta, kjpt )
326	!!----------------------------------------------------------------------
327	!! * ROUTINE tra_adv_tvd_zts *
328	!!
329	!! ** Purpose : Compute the now trend due to total advection of
330	!! tracers and add it to the general trend of tracer equations
331	!!
332	!! ** Method : TVD ZTS scheme, i.e. 2nd order centered scheme with
333	!! corrected flux (monotonic correction). This version use sub-
334	!! timestepping for the vertical advection which increases stability
335	!! when vertical metrics are small.
336	!! note: - this advection scheme needs a leap-frog time scheme
337	!!
338	!! ** Action : - update (pta) with the now advective tracer trends
339	!! - save the trends
340	!!----------------------------------------------------------------------
341	USE oce , ONLY: zwx => ua , zwy => va ! (ua,va) used as workspace
342	!
343	INTEGER , INTENT(in ) :: kt ! ocean time-step index
344	INTEGER , INTENT(in ) :: kit000 ! first time step index
345	CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator)
346	INTEGER , INTENT(in ) :: kjpt ! number of tracers
347	REAL(wp), DIMENSION( jpk ), INTENT(in ) :: p2dt ! vertical profile of tracer time-step
348	REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pun, pvn, pwn ! 3 ocean velocity components
349	REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb, ptn ! before and now tracer fields
350	REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend
351	!
352	REAL(wp), DIMENSION( jpk ) :: zts ! length of sub-timestep for vertical advection
353	REAL(wp), DIMENSION( jpk ) :: zr_p2dt ! reciprocal of tracer timestep
354	INTEGER :: ji, jj, jk, jl, jn ! dummy loop indices
355	INTEGER :: jnzts = 5 ! number of sub-timesteps for vertical advection
356	INTEGER :: jtb, jtn, jta ! sub timestep pointers for leap-frog/euler forward steps
357	INTEGER :: jtaken ! toggle for collecting appropriate fluxes from sub timesteps
358	REAL(wp) :: z_rzts ! Fractional length of Euler forward sub-timestep for vertical advection
359	REAL(wp) :: z2dtt, zbtr, ztra ! local scalar
360	REAL(wp) :: zfp_ui, zfp_vj, zfp_wk ! - -
361	REAL(wp) :: zfm_ui, zfm_vj, zfm_wk ! - -
362	REAL(wp), ALLOCATABLE, DIMENSION(:,: ) :: zwx_sav , zwy_sav
363	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zwi, zwz, zhdiv, zwz_sav, zwzts
364	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrdx, ztrdy, ztrdz
365	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zptry
366	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) :: ztrs
367	!!----------------------------------------------------------------------
368	!
369	IF( nn_timing == 1 ) CALL timing_start('tra_adv_tvd_zts')
370	!
371	ALLOCATE(zwx_sav(1:jpi, 1:jpj))
372	ALLOCATE(zwy_sav(1:jpi, 1:jpj))
373	ALLOCATE(zwi(1:jpi, 1:jpj, 1:jpk))
374	ALLOCATE(zwz(1:jpi, 1:jpj, 1:jpk))
375	ALLOCATE(zhdiv(1:jpi, 1:jpj, 1:jpk))
376	ALLOCATE(zwz_sav(1:jpi, 1:jpj, 1:jpk))
377	ALLOCATE(zwzts(1:jpi, 1:jpj, 1:jpk))
378	ALLOCATE(ztrs(1:jpi, 1:jpj, 1:jpk, 1:kjpt+1))
379	!
380	IF( kt == kit000 ) THEN
381	IF(lwp) WRITE(numout,*)
382	IF(lwp) WRITE(numout,*) 'tra_adv_tvd_zts : TVD ZTS advection scheme on ', cdtype
383	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
384	ENDIF
385	!
386	l_trd = .FALSE.
387	IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE.
388	!
389	IF( l_trd ) THEN
390	ALLOCATE(ztrdx(1:jpi, 1:jpj, 1:jpk))
391	ALLOCATE(ztrdy(1:jpi, 1:jpj, 1:jpk))
392	ALLOCATE(ztrdz(1:jpi, 1:jpj, 1:jpk))
393	ztrdx(:,:,:) = 0._wp ; ztrdy(:,:,:) = 0._wp ; ztrdz(:,:,:) = 0._wp
394	ENDIF
395	!
396	IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN
397	ALLOCATE(zptry(1:jpi, 1:jpj, 1:jpk))
398	zptry(:,:,:) = 0._wp
399	ENDIF
400	!
401	zwi(:,:,:) = 0._wp
402	z_rzts = 1._wp / REAL( jnzts, wp )
403	zr_p2dt(:) = 1._wp / p2dt(:)
404	!
405	! ! ===========
406	DO jn = 1, kjpt ! tracer loop
407	! ! ===========
408	! 1. Bottom value : flux set to zero
409	! ----------------------------------
410	zwx(:,:,jpk) = 0._wp ; zwz(:,:,jpk) = 0._wp
411	zwy(:,:,jpk) = 0._wp ; zwi(:,:,jpk) = 0._wp
412
413	! 2. upstream advection with initial mass fluxes & intermediate update
414	! --------------------------------------------------------------------
415	! upstream tracer flux in the i and j direction
416	DO jk = 1, jpkm1
417	DO jj = 1, jpjm1
418	DO ji = 1, fs_jpim1 ! vector opt.
419	! upstream scheme
420	zfp_ui = pun(ji,jj,jk) + ABS( pun(ji,jj,jk) )
421	zfm_ui = pun(ji,jj,jk) - ABS( pun(ji,jj,jk) )
422	zfp_vj = pvn(ji,jj,jk) + ABS( pvn(ji,jj,jk) )
423	zfm_vj = pvn(ji,jj,jk) - ABS( pvn(ji,jj,jk) )
424	zwx(ji,jj,jk) = 0.5_wp * ( zfp_ui * ptb(ji,jj,jk,jn) + zfm_ui * ptb(ji+1,jj ,jk,jn) )
425	zwy(ji,jj,jk) = 0.5_wp * ( zfp_vj * ptb(ji,jj,jk,jn) + zfm_vj * ptb(ji ,jj+1,jk,jn) )
426	END DO
427	END DO
428	END DO
429
430	! upstream tracer flux in the k direction
431	! Interior value
432	DO jk = 2, jpkm1
433	DO jj = 1, jpj
434	DO ji = 1, jpi
435	zfp_wk = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) )
436	zfm_wk = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) )
437	zwz(ji,jj,jk) = 0.5_wp * ( zfp_wk * ptb(ji,jj,jk,jn) + zfm_wk * ptb(ji,jj,jk-1,jn) )
438	END DO
439	END DO
440	END DO
441	! Surface value
442	IF( lk_vvl ) THEN
443	IF ( ln_isfcav ) THEN
444	DO jj = 1, jpj
445	DO ji = 1, jpi
446	zwz(ji,jj, mikt(ji,jj) ) = 0.e0 ! volume variable + isf
447	END DO
448	END DO
449	ELSE
450	zwz(:,:,1) = 0.e0 ! volume variable + no isf
451	END IF
452	ELSE
453	IF ( ln_isfcav ) THEN
454	DO jj = 1, jpj
455	DO ji = 1, jpi
456	zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptb(ji,jj,mikt(ji,jj),jn) ! linear free surface + isf
457	END DO
458	END DO
459	ELSE
460	zwz(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn) ! linear free surface + no isf
461	END IF
462	ENDIF
463
464	! total advective trend
465	DO jk = 1, jpkm1
466	z2dtt = p2dt(jk)
467	DO jj = 2, jpjm1
468	DO ji = fs_2, fs_jpim1 ! vector opt.
469	! total intermediate advective trends
470	ztra = - ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) &
471	& + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) &
472	& + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) / e1e2t(ji,jj)
473	! update and guess with monotonic sheme
474	pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra / fse3t_n(ji,jj,jk) * tmask(ji,jj,jk)
475	zwi(ji,jj,jk) = ( fse3t_b(ji,jj,jk) * ptb(ji,jj,jk,jn) + z2dtt * ztra ) / fse3t_a(ji,jj,jk) * tmask(ji,jj,jk)
476	END DO
477	END DO
478	END DO
479	! ! Lateral boundary conditions on zwi (unchanged sign)
480	CALL lbc_lnk( zwi, 'T', 1. )
481
482	! ! trend diagnostics (contribution of upstream fluxes)
483	IF( l_trd ) THEN
484	! store intermediate advective trends
485	ztrdx(:,:,:) = zwx(:,:,:) ; ztrdy(:,:,:) = zwy(:,:,:) ; ztrdz(:,:,:) = zwz(:,:,:)
486	END IF
487	! ! "Poleward" heat and salt transports (contribution of upstream fluxes)
488	IF( cdtype == 'TRA' .AND. ln_diaptr ) zptry(:,:,:) = zwy(:,:,:)
489
490	! 3. antidiffusive flux : high order minus low order
491	! --------------------------------------------------
492	! antidiffusive flux on i and j
493	!
494	DO jk = 1, jpkm1
495	!
496	DO jj = 1, jpjm1
497	DO ji = 1, fs_jpim1 ! vector opt.
498	zwx_sav(ji,jj) = zwx(ji,jj,jk)
499	zwy_sav(ji,jj) = zwy(ji,jj,jk)
500
501	zwx(ji,jj,jk) = 0.5_wp * pun(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji+1,jj,jk,jn) )
502	zwy(ji,jj,jk) = 0.5_wp * pvn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj+1,jk,jn) )
503	END DO
504	END DO
505
506	DO jj = 2, jpjm1 ! partial horizontal divergence
507	DO ji = fs_2, fs_jpim1
508	zhdiv(ji,jj,jk) = ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk) &
509	& + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk) )
510	END DO
511	END DO
512
513	DO jj = 1, jpjm1
514	DO ji = 1, fs_jpim1 ! vector opt.
515	zwx(ji,jj,jk) = zwx(ji,jj,jk) - zwx_sav(ji,jj)
516	zwy(ji,jj,jk) = zwy(ji,jj,jk) - zwy_sav(ji,jj)
517	END DO
518	END DO
519	END DO
520
521	! antidiffusive flux on k
522	zwz(:,:,1) = 0._wp ! Surface value
523	zwz_sav(:,:,:) = zwz(:,:,:)
524	!
525	ztrs(:,:,:,1) = ptb(:,:,:,jn)
526	ztrs(:,:,1,2) = ptb(:,:,1,jn)
527	ztrs(:,:,1,3) = ptb(:,:,1,jn)
528	zwzts(:,:,:) = 0._wp
529
530	DO jl = 1, jnzts ! Start of sub timestepping loop
531
532	IF( jl == 1 ) THEN ! Euler forward to kick things off
533	jtb = 1 ; jtn = 1 ; jta = 2
534	zts(:) = p2dt(:) * z_rzts
535	jtaken = MOD( jnzts + 1 , 2) ! Toggle to collect every second flux
536	! starting at jl =1 if jnzts is odd;
537	! starting at jl =2 otherwise
538	ELSEIF( jl == 2 ) THEN ! First leapfrog step
539	jtb = 1 ; jtn = 2 ; jta = 3
540	zts(:) = 2._wp * p2dt(:) * z_rzts
541	ELSE ! Shuffle pointers for subsequent leapfrog steps
542	jtb = MOD(jtb,3) + 1
543	jtn = MOD(jtn,3) + 1
544	jta = MOD(jta,3) + 1
545	ENDIF
546	DO jk = 2, jpkm1 ! Interior value
547	DO jj = 2, jpjm1
548	DO ji = fs_2, fs_jpim1
549	zwz(ji,jj,jk) = 0.5_wp * pwn(ji,jj,jk) * ( ztrs(ji,jj,jk,jtn) + ztrs(ji,jj,jk-1,jtn) )
550	IF( jtaken == 0 ) zwzts(ji,jj,jk) = zwzts(ji,jj,jk) + zwz(ji,jj,jk)*zts(jk) ! Accumulate time-weighted vertcal flux
551	END DO
552	END DO
553	END DO
554
555	jtaken = MOD( jtaken + 1 , 2 )
556
557	DO jk = 2, jpkm1 ! Interior value
558	DO jj = 2, jpjm1
559	DO ji = fs_2, fs_jpim1
560	zbtr = 1._wp / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
561	! total advective trends
562	ztra = - zbtr * ( zhdiv(ji,jj,jk) + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) )
563	ztrs(ji,jj,jk,jta) = ztrs(ji,jj,jk,jtb) + zts(jk) * ztra
564	END DO
565	END DO
566	END DO
567
568	END DO
569
570	DO jk = 2, jpkm1 ! Anti-diffusive vertical flux using average flux from the sub-timestepping
571	DO jj = 2, jpjm1
572	DO ji = fs_2, fs_jpim1
573	zwz(ji,jj,jk) = zwzts(ji,jj,jk) * zr_p2dt(jk) - zwz_sav(ji,jj,jk)
574	END DO
575	END DO
576	END DO
577	CALL lbc_lnk( zwx, 'U', -1. ) ; CALL lbc_lnk( zwy, 'V', -1. ) ! Lateral bondary conditions
578	CALL lbc_lnk( zwz, 'W', 1. )
579
580	! 4. monotonicity algorithm
581	! -------------------------
582	CALL nonosc( ptb(:,:,:,jn), zwx, zwy, zwz, zwi, p2dt )
583
584
585	! 5. final trend with corrected fluxes
586	! ------------------------------------
587	DO jk = 1, jpkm1
588	DO jj = 2, jpjm1
589	DO ji = fs_2, fs_jpim1 ! vector opt.
590	zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
591	! total advective trends
592	ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) &
593	& + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) &
594	& + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) )
595	! add them to the general tracer trends
596	pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
597	END DO
598	END DO
599	END DO
600
601	! ! trend diagnostics (contribution of upstream fluxes)
602	IF( l_trd ) THEN
603	ztrdx(:,:,:) = ztrdx(:,:,:) + zwx(:,:,:) ! <<< Add to previously computed
604	ztrdy(:,:,:) = ztrdy(:,:,:) + zwy(:,:,:) ! <<< Add to previously computed
605	ztrdz(:,:,:) = ztrdz(:,:,:) + zwz(:,:,:) ! <<< Add to previously computed
606
607	CALL trd_tra( kt, cdtype, jn, jptra_xad, ztrdx, pun, ptn(:,:,:,jn) )
608	CALL trd_tra( kt, cdtype, jn, jptra_yad, ztrdy, pvn, ptn(:,:,:,jn) )
609	CALL trd_tra( kt, cdtype, jn, jptra_zad, ztrdz, pwn, ptn(:,:,:,jn) )
610	END IF
611	! ! "Poleward" heat and salt transports (contribution of upstream fluxes)
612	IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN
613	zptry(:,:,:) = zptry(:,:,:) + zwy(:,:,:)
614	CALL dia_ptr_ohst_components( jn, 'adv', zptry(:,:,:) )
615	ENDIF
616	!
617	END DO
618	!
619	DEALLOCATE(zwi)
620	DEALLOCATE(zwz)
621	DEALLOCATE(zhdiv)
622	DEALLOCATE(zwz_sav)
623	DEALLOCATE(zwzts)
624	DEALLOCATE(ztrs )
625	DEALLOCATE(zwx_sav)
626	DEALLOCATE(zwy_sav )
627
628	IF( l_trd ) THEN
629	DEALLOCATE(ztrdx)
630	DEALLOCATE(ztrdy)
631	DEALLOCATE(ztrdz)
632	END IF
633
634	IF( cdtype == 'TRA' .AND. ln_diaptr ) DEALLOCATE(zptry )
635	!
636	IF( nn_timing == 1 ) CALL timing_stop('tra_adv_tvd_zts')
637	!
638	END SUBROUTINE tra_adv_tvd_zts
639
640
641	SUBROUTINE nonosc( pbef, paa, pbb, pcc, paft, p2dt )
642	!!---------------------------------------------------------------------
643	!! * ROUTINE nonosc *
644	!!
645	!! ** Purpose : compute monotonic tracer fluxes from the upstream
646	!! scheme and the before field by a nonoscillatory algorithm
647	!!
648	!! ** Method : ... ???
649	!! warning : pbef and paft must be masked, but the boundaries
650	!! conditions on the fluxes are not necessary zalezak (1979)
651	!! drange (1995) multi-dimensional forward-in-time and upstream-
652	!! in-space based differencing for fluid
653	!!----------------------------------------------------------------------
654	REAL(wp), DIMENSION(jpk) , INTENT(in ) :: p2dt ! vertical profile of tracer time-step
655	REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(in ) :: pbef, paft ! before & after field
656	REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(inout) :: paa, pbb, pcc ! monotonic fluxes in the 3 directions
657	!
658	INTEGER :: ji, jj, jk ! dummy loop indices
659	INTEGER :: ikm1 ! local integer
660	REAL(wp) :: zpos, zneg, zbt, za, zb, zc, zbig, zrtrn, z2dtt ! local scalars
661	REAL(wp) :: zau, zbu, zcu, zav, zbv, zcv, zup, zdo ! - -
662	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zbetup, zbetdo, zbup, zbdo
663	!!----------------------------------------------------------------------
664	!
665	IF( nn_timing == 1 ) CALL timing_start('nonosc')
666	!
667	ALLOCATE(zbetup(1:jpi, 1:jpj, 1:jpk))
668	ALLOCATE(zbetdo(1:jpi, 1:jpj, 1:jpk))
669	ALLOCATE(zbup(1:jpi, 1:jpj, 1:jpk))
670	ALLOCATE(zbdo(1:jpi, 1:jpj, 1:jpk))
671	!
672	zbig = 1.e+40_wp
673	zrtrn = 1.e-15_wp
674	zbetup(:,:,:) = 0._wp ; zbetdo(:,:,:) = 0._wp
675
676	! Search local extrema
677	! --------------------
678	! max/min of pbef & paft with large negative/positive value (-/+zbig) inside land
679	zbup = MAX( pbef * tmask - zbig * ( 1._wp - tmask ), &
680	& paft * tmask - zbig * ( 1._wp - tmask ) )
681	zbdo = MIN( pbef * tmask + zbig * ( 1._wp - tmask ), &
682	& paft * tmask + zbig * ( 1._wp - tmask ) )
683
684	DO jk = 1, jpkm1
685	ikm1 = MAX(jk-1,1)
686	z2dtt = p2dt(jk)
687	DO jj = 2, jpjm1
688	DO ji = fs_2, fs_jpim1 ! vector opt.
689
690	! search maximum in neighbourhood
691	zup = MAX( zbup(ji ,jj ,jk ), &
692	& zbup(ji-1,jj ,jk ), zbup(ji+1,jj ,jk ), &
693	& zbup(ji ,jj-1,jk ), zbup(ji ,jj+1,jk ), &
694	& zbup(ji ,jj ,ikm1), zbup(ji ,jj ,jk+1) )
695
696	! search minimum in neighbourhood
697	zdo = MIN( zbdo(ji ,jj ,jk ), &
698	& zbdo(ji-1,jj ,jk ), zbdo(ji+1,jj ,jk ), &
699	& zbdo(ji ,jj-1,jk ), zbdo(ji ,jj+1,jk ), &
700	& zbdo(ji ,jj ,ikm1), zbdo(ji ,jj ,jk+1) )
701
702	! positive part of the flux
703	zpos = MAX( 0., paa(ji-1,jj ,jk ) ) - MIN( 0., paa(ji ,jj ,jk ) ) &
704	& + MAX( 0., pbb(ji ,jj-1,jk ) ) - MIN( 0., pbb(ji ,jj ,jk ) ) &
705	& + MAX( 0., pcc(ji ,jj ,jk+1) ) - MIN( 0., pcc(ji ,jj ,jk ) )
706
707	! negative part of the flux
708	zneg = MAX( 0., paa(ji ,jj ,jk ) ) - MIN( 0., paa(ji-1,jj ,jk ) ) &
709	& + MAX( 0., pbb(ji ,jj ,jk ) ) - MIN( 0., pbb(ji ,jj-1,jk ) ) &
710	& + MAX( 0., pcc(ji ,jj ,jk ) ) - MIN( 0., pcc(ji ,jj ,jk+1) )
711
712	! up & down beta terms
713	zbt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) / z2dtt
714	zbetup(ji,jj,jk) = ( zup - paft(ji,jj,jk) ) / ( zpos + zrtrn ) * zbt
715	zbetdo(ji,jj,jk) = ( paft(ji,jj,jk) - zdo ) / ( zneg + zrtrn ) * zbt
716	END DO
717	END DO
718	END DO
719	CALL lbc_lnk( zbetup, 'T', 1. ) ; CALL lbc_lnk( zbetdo, 'T', 1. ) ! lateral boundary cond. (unchanged sign)
720
721	! 3. monotonic flux in the i & j direction (paa & pbb)
722	! ----------------------------------------
723	DO jk = 1, jpkm1
724	DO jj = 2, jpjm1
725	DO ji = fs_2, fs_jpim1 ! vector opt.
726	zau = MIN( 1._wp, zbetdo(ji,jj,jk), zbetup(ji+1,jj,jk) )
727	zbu = MIN( 1._wp, zbetup(ji,jj,jk), zbetdo(ji+1,jj,jk) )
728	zcu = ( 0.5 + SIGN( 0.5 , paa(ji,jj,jk) ) )
729	paa(ji,jj,jk) = paa(ji,jj,jk) * ( zcu * zau + ( 1._wp - zcu) * zbu )
730
731	zav = MIN( 1._wp, zbetdo(ji,jj,jk), zbetup(ji,jj+1,jk) )
732	zbv = MIN( 1._wp, zbetup(ji,jj,jk), zbetdo(ji,jj+1,jk) )
733	zcv = ( 0.5 + SIGN( 0.5 , pbb(ji,jj,jk) ) )
734	pbb(ji,jj,jk) = pbb(ji,jj,jk) * ( zcv * zav + ( 1._wp - zcv) * zbv )
735
736	! monotonic flux in the k direction, i.e. pcc
737	! -------------------------------------------
738	za = MIN( 1., zbetdo(ji,jj,jk+1), zbetup(ji,jj,jk) )
739	zb = MIN( 1., zbetup(ji,jj,jk+1), zbetdo(ji,jj,jk) )
740	zc = ( 0.5 + SIGN( 0.5 , pcc(ji,jj,jk+1) ) )
741	pcc(ji,jj,jk+1) = pcc(ji,jj,jk+1) * ( zc * za + ( 1._wp - zc) * zb )
742	END DO
743	END DO
744	END DO
745	CALL lbc_lnk( paa, 'U', -1. ) ; CALL lbc_lnk( pbb, 'V', -1. ) ! lateral boundary condition (changed sign)
746	!
747	DEALLOCATE(zbetup)
748	DEALLOCATE(zbetdo)
749	DEALLOCATE(zbup)
750	DEALLOCATE(zbdo)
751	!
752	IF( nn_timing == 1 ) CALL timing_stop('nonosc')
753	!
754	END SUBROUTINE nonosc
755
756	!!======================================================================
757	END MODULE traadv_tvd

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

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

Context Navigation

source: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_tvd.F90 @ 7771

Download in other formats: