Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

traadv_fct.F90 @ 12088

Last change on this file since 12088 was 12088, checked in by deazer, 4 years ago
Updating GO8 Package branch to bring in required BDY bug fixes frouse with CO8 The mirror branch is already updated to have this change, where we merge in the mirror to the package branch
File size: 38.0 KB

Line
1	MODULE traadv_fct
2	!!==============================================================================
3	!! * MODULE traadv_fct *
4	!! Ocean tracers: horizontal & vertical advective trend (2nd/4th order Flux Corrected Transport method)
5	!!==============================================================================
6	!! History : 3.7 ! 2015-09 (L. Debreu, G. Madec) original code (inspired from traadv_tvd.F90)
7	!!----------------------------------------------------------------------
8
9	!!----------------------------------------------------------------------
10	!! tra_adv_fct : update the tracer trend with a 3D advective trends using a 2nd or 4th order FCT scheme
11	!! with sub-time-stepping in the vertical direction
12	!! nonosc : compute monotonic tracer fluxes by a non-oscillatory algorithm
13	!! interp_4th_cpt : 4th order compact scheme for the vertical component of the advection
14	!!----------------------------------------------------------------------
15	USE oce ! ocean dynamics and active tracers
16	USE dom_oce ! ocean space and time domain
17	USE trc_oce ! share passive tracers/Ocean variables
18	USE trd_oce ! trends: ocean variables
19	USE trdtra ! tracers trends
20	USE diaptr ! poleward transport diagnostics
21	USE diaar5 ! AR5 diagnostics
22	USE phycst , ONLY : rau0_rcp
23	USE zdf_oce , ONLY : ln_zad_Aimp
24	!
25	USE in_out_manager ! I/O manager
26	USE iom !
27	USE lib_mpp ! MPP library
28	USE lbclnk ! ocean lateral boundary condition (or mpp link)
29	USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
30
31	IMPLICIT NONE
32	PRIVATE
33
34	PUBLIC tra_adv_fct ! called by traadv.F90
35	PUBLIC interp_4th_cpt ! called by traadv_cen.F90
36
37	LOGICAL :: l_trd ! flag to compute trends
38	LOGICAL :: l_ptr ! flag to compute poleward transport
39	LOGICAL :: l_hst ! flag to compute heat/salt transport
40	REAL(wp) :: r1_6 = 1._wp / 6._wp ! =1/6
41
42	! ! tridiag solver associated indices:
43	INTEGER, PARAMETER :: np_NH = 0 ! Neumann homogeneous boundary condition
44	INTEGER, PARAMETER :: np_CEN2 = 1 ! 2nd order centered boundary condition
45
46	!! * Substitutions
47	# include "vectopt_loop_substitute.h90"
48	!!----------------------------------------------------------------------
49	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
50	!! $Id$
51	!! Software governed by the CeCILL license (see ./LICENSE)
52	!!----------------------------------------------------------------------
53	CONTAINS
54
55	SUBROUTINE tra_adv_fct( kt, kit000, cdtype, p2dt, pun, pvn, pwn, &
56	& ptb, ptn, pta, kjpt, kn_fct_h, kn_fct_v )
57	!!----------------------------------------------------------------------
58	!! * ROUTINE tra_adv_fct *
59	!!
60	!! ** Purpose : Compute the now trend due to total advection of tracers
61	!! and add it to the general trend of tracer equations
62	!!
63	!! ** Method : - 2nd or 4th FCT scheme on the horizontal direction
64	!! (choice through the value of kn_fct)
65	!! - on the vertical the 4th order is a compact scheme
66	!! - corrected flux (monotonic correction)
67	!!
68	!! ** Action : - update pta with the now advective tracer trends
69	!! - send trends to trdtra module for further diagnostics (l_trdtra=T)
70	!! - htr_adv, str_adv : poleward advective heat and salt transport (ln_diaptr=T)
71	!!----------------------------------------------------------------------
72	INTEGER , INTENT(in ) :: kt ! ocean time-step index
73	INTEGER , INTENT(in ) :: kit000 ! first time step index
74	CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator)
75	INTEGER , INTENT(in ) :: kjpt ! number of tracers
76	INTEGER , INTENT(in ) :: kn_fct_h ! order of the FCT scheme (=2 or 4)
77	INTEGER , INTENT(in ) :: kn_fct_v ! order of the FCT scheme (=2 or 4)
78	REAL(wp) , INTENT(in ) :: p2dt ! 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	REAL(wp) :: ztra ! local scalar
85	REAL(wp) :: zfp_ui, zfp_vj, zfp_wk, zC2t_u, zC4t_u ! - -
86	REAL(wp) :: zfm_ui, zfm_vj, zfm_wk, zC2t_v, zC4t_v ! - -
87	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwi, zwx, zwy, zwz, ztu, ztv, zltu, zltv, ztw
88	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ztrdx, ztrdy, ztrdz, zptry
89	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zwinf, zwdia, zwsup
90	LOGICAL :: ll_zAimp ! flag to apply adaptive implicit vertical advection
91	!!----------------------------------------------------------------------
92	!
93	IF( kt == kit000 ) THEN
94	IF(lwp) WRITE(numout,*)
95	IF(lwp) WRITE(numout,*) 'tra_adv_fct : FCT advection scheme on ', cdtype
96	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
97	ENDIF
98	!
99	l_trd = .FALSE. ! set local switches
100	l_hst = .FALSE.
101	l_ptr = .FALSE.
102	ll_zAimp = .FALSE.
103	IF( ( cdtype =='TRA' .AND. l_trdtra ) .OR. ( cdtype =='TRC' .AND. l_trdtrc ) ) l_trd = .TRUE.
104	IF( cdtype =='TRA' .AND. ln_diaptr ) l_ptr = .TRUE.
105	IF( cdtype =='TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. &
106	& iom_use("uadv_salttr") .OR. iom_use("vadv_salttr") ) ) l_hst = .TRUE.
107	!
108	IF( l_trd .OR. l_hst ) THEN
109	ALLOCATE( ztrdx(jpi,jpj,jpk), ztrdy(jpi,jpj,jpk), ztrdz(jpi,jpj,jpk) )
110	ztrdx(:,:,:) = 0._wp ; ztrdy(:,:,:) = 0._wp ; ztrdz(:,:,:) = 0._wp
111	ENDIF
112	!
113	IF( l_ptr ) THEN
114	ALLOCATE( zptry(jpi,jpj,jpk) )
115	zptry(:,:,:) = 0._wp
116	ENDIF
117	! ! surface & bottom value : flux set to zero one for all
118	zwz(:,:, 1 ) = 0._wp
119	zwx(:,:,jpk) = 0._wp ; zwy(:,:,jpk) = 0._wp ; zwz(:,:,jpk) = 0._wp
120	!
121	zwi(:,:,:) = 0._wp
122	!
123	! If adaptive vertical advection, check if it is needed on this PE at this time
124	IF( ln_zad_Aimp ) THEN
125	IF( MAXVAL( ABS( wi(:,:,:) ) ) > 0._wp ) ll_zAimp = .TRUE.
126	END IF
127	! If active adaptive vertical advection, build tridiagonal matrix
128	IF( ll_zAimp ) THEN
129	ALLOCATE(zwdia(jpi,jpj,jpk), zwinf(jpi,jpj,jpk),zwsup(jpi,jpj,jpk))
130	DO jk = 1, jpkm1
131	DO jj = 2, jpjm1
132	DO ji = fs_2, fs_jpim1 ! vector opt. (ensure same order of calculation as below if wi=0.)
133	zwdia(ji,jj,jk) = 1._wp + p2dt * ( MAX( wi(ji,jj,jk ) , 0._wp ) - MIN( wi(ji,jj,jk+1) , 0._wp ) ) / e3t_a(ji,jj,jk)
134	zwinf(ji,jj,jk) = p2dt * MIN( wi(ji,jj,jk ) , 0._wp ) / e3t_a(ji,jj,jk)
135	zwsup(ji,jj,jk) = -p2dt * MAX( wi(ji,jj,jk+1) , 0._wp ) / e3t_a(ji,jj,jk)
136	END DO
137	END DO
138	END DO
139	END IF
140	!
141	DO jn = 1, kjpt !== loop over the tracers ==!
142	!
143	! !== upstream advection with initial mass fluxes & intermediate update ==!
144	! !* upstream tracer flux in the i and j direction
145	DO jk = 1, jpkm1
146	DO jj = 1, jpjm1
147	DO ji = 1, fs_jpim1 ! vector opt.
148	! upstream scheme
149	zfp_ui = pun(ji,jj,jk) + ABS( pun(ji,jj,jk) )
150	zfm_ui = pun(ji,jj,jk) - ABS( pun(ji,jj,jk) )
151	zfp_vj = pvn(ji,jj,jk) + ABS( pvn(ji,jj,jk) )
152	zfm_vj = pvn(ji,jj,jk) - ABS( pvn(ji,jj,jk) )
153	zwx(ji,jj,jk) = 0.5 * ( zfp_ui * ptb(ji,jj,jk,jn) + zfm_ui * ptb(ji+1,jj ,jk,jn) )
154	zwy(ji,jj,jk) = 0.5 * ( zfp_vj * ptb(ji,jj,jk,jn) + zfm_vj * ptb(ji ,jj+1,jk,jn) )
155	END DO
156	END DO
157	END DO
158	! !* upstream tracer flux in the k direction *!
159	DO jk = 2, jpkm1 ! Interior value ( multiplied by wmask)
160	DO jj = 1, jpj
161	DO ji = 1, jpi
162	zfp_wk = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) )
163	zfm_wk = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) )
164	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)
165	END DO
166	END DO
167	END DO
168	IF( ln_linssh ) THEN ! top ocean value (only in linear free surface as zwz has been w-masked)
169	IF( ln_isfcav ) THEN ! top of the ice-shelf cavities and at the ocean surface
170	DO jj = 1, jpj
171	DO ji = 1, jpi
172	zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptb(ji,jj,mikt(ji,jj),jn) ! linear free surface
173	END DO
174	END DO
175	ELSE ! no cavities: only at the ocean surface
176	zwz(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn)
177	ENDIF
178	ENDIF
179	!
180	DO jk = 1, jpkm1 !* trend and after field with monotonic scheme
181	DO jj = 2, jpjm1
182	DO ji = fs_2, fs_jpim1 ! vector opt.
183	! ! total intermediate advective trends
184	ztra = - ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) &
185	& + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) &
186	& + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) * r1_e1e2t(ji,jj)
187	! ! update and guess with monotonic sheme
188	pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra / e3t_n(ji,jj,jk) * tmask(ji,jj,jk)
189	zwi(ji,jj,jk) = ( e3t_b(ji,jj,jk) * ptb(ji,jj,jk,jn) + p2dt * ztra ) / e3t_a(ji,jj,jk) * tmask(ji,jj,jk)
190	END DO
191	END DO
192	END DO
193
194	IF ( ll_zAimp ) THEN
195	CALL tridia_solver( zwdia, zwsup, zwinf, zwi, zwi , 0 )
196	!
197	ztw(:,:,1) = 0._wp ; ztw(:,:,jpk) = 0._wp ;
198	DO jk = 2, jpkm1 ! Interior value ( multiplied by wmask)
199	DO jj = 2, jpjm1
200	DO ji = fs_2, fs_jpim1 ! vector opt.
201	zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) )
202	zfm_wk = wi(ji,jj,jk) - ABS( wi(ji,jj,jk) )
203	ztw(ji,jj,jk) = 0.5 * e1e2t(ji,jj) * ( zfp_wk * zwi(ji,jj,jk) + zfm_wk * zwi(ji,jj,jk-1) ) * wmask(ji,jj,jk)
204	zwz(ji,jj,jk) = zwz(ji,jj,jk) + ztw(ji,jj,jk) ! update vertical fluxes
205	END DO
206	END DO
207	END DO
208	DO jk = 1, jpkm1
209	DO jj = 2, jpjm1
210	DO ji = fs_2, fs_jpim1 ! vector opt.
211	pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) - ( ztw(ji,jj,jk) - ztw(ji ,jj ,jk+1) ) &
212	& * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
213	END DO
214	END DO
215	END DO
216	!
217	END IF
218	!
219	IF( l_trd .OR. l_hst ) THEN ! trend diagnostics (contribution of upstream fluxes)
220	ztrdx(:,:,:) = zwx(:,:,:) ; ztrdy(:,:,:) = zwy(:,:,:) ; ztrdz(:,:,:) = zwz(:,:,:)
221	END IF
222	! ! "Poleward" heat and salt transports (contribution of upstream fluxes)
223	IF( l_ptr ) zptry(:,:,:) = zwy(:,:,:)
224	!
225	! !== anti-diffusive flux : high order minus low order ==!
226	!
227	SELECT CASE( kn_fct_h ) !* horizontal anti-diffusive fluxes
228	!
229	CASE( 2 ) !- 2nd order centered
230	DO jk = 1, jpkm1
231	DO jj = 1, jpjm1
232	DO ji = 1, fs_jpim1 ! vector opt.
233	zwx(ji,jj,jk) = 0.5_wp * pun(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji+1,jj,jk,jn) ) - zwx(ji,jj,jk)
234	zwy(ji,jj,jk) = 0.5_wp * pvn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj+1,jk,jn) ) - zwy(ji,jj,jk)
235	END DO
236	END DO
237	END DO
238	!
239	CASE( 4 ) !- 4th order centered
240	zltu(:,:,jpk) = 0._wp ! Bottom value : flux set to zero
241	zltv(:,:,jpk) = 0._wp
242	DO jk = 1, jpkm1 ! Laplacian
243	DO jj = 1, jpjm1 ! 1st derivative (gradient)
244	DO ji = 1, fs_jpim1 ! vector opt.
245	ztu(ji,jj,jk) = ( ptn(ji+1,jj ,jk,jn) - ptn(ji,jj,jk,jn) ) * umask(ji,jj,jk)
246	ztv(ji,jj,jk) = ( ptn(ji ,jj+1,jk,jn) - ptn(ji,jj,jk,jn) ) * vmask(ji,jj,jk)
247	END DO
248	END DO
249	DO jj = 2, jpjm1 ! 2nd derivative * 1/ 6
250	DO ji = fs_2, fs_jpim1 ! vector opt.
251	zltu(ji,jj,jk) = ( ztu(ji,jj,jk) + ztu(ji-1,jj,jk) ) * r1_6
252	zltv(ji,jj,jk) = ( ztv(ji,jj,jk) + ztv(ji,jj-1,jk) ) * r1_6
253	END DO
254	END DO
255	END DO
256	CALL lbc_lnk_multi( 'traadv_fct', zltu, 'T', 1. , zltv, 'T', 1. ) ! Lateral boundary cond. (unchanged sgn)
257	!
258	DO jk = 1, jpkm1 ! Horizontal advective fluxes
259	DO jj = 1, jpjm1
260	DO ji = 1, fs_jpim1 ! vector opt.
261	zC2t_u = ptn(ji,jj,jk,jn) + ptn(ji+1,jj ,jk,jn) ! 2 x C2 interpolation of T at u- & v-points
262	zC2t_v = ptn(ji,jj,jk,jn) + ptn(ji ,jj+1,jk,jn)
263	! ! C4 minus upstream advective fluxes
264	zwx(ji,jj,jk) = 0.5_wp * pun(ji,jj,jk) * ( zC2t_u + zltu(ji,jj,jk) - zltu(ji+1,jj,jk) ) - zwx(ji,jj,jk)
265	zwy(ji,jj,jk) = 0.5_wp * pvn(ji,jj,jk) * ( zC2t_v + zltv(ji,jj,jk) - zltv(ji,jj+1,jk) ) - zwy(ji,jj,jk)
266	END DO
267	END DO
268	END DO
269	!
270	CASE( 41 ) !- 4th order centered ==>> !!gm coding attempt need to be tested
271	ztu(:,:,jpk) = 0._wp ! Bottom value : flux set to zero
272	ztv(:,:,jpk) = 0._wp
273	DO jk = 1, jpkm1 ! 1st derivative (gradient)
274	DO jj = 1, jpjm1
275	DO ji = 1, fs_jpim1 ! vector opt.
276	ztu(ji,jj,jk) = ( ptn(ji+1,jj ,jk,jn) - ptn(ji,jj,jk,jn) ) * umask(ji,jj,jk)
277	ztv(ji,jj,jk) = ( ptn(ji ,jj+1,jk,jn) - ptn(ji,jj,jk,jn) ) * vmask(ji,jj,jk)
278	END DO
279	END DO
280	END DO
281	CALL lbc_lnk_multi( 'traadv_fct', ztu, 'U', -1. , ztv, 'V', -1. ) ! Lateral boundary cond. (unchanged sgn)
282	!
283	DO jk = 1, jpkm1 ! Horizontal advective fluxes
284	DO jj = 2, jpjm1
285	DO ji = 2, fs_jpim1 ! vector opt.
286	zC2t_u = ptn(ji,jj,jk,jn) + ptn(ji+1,jj ,jk,jn) ! 2 x C2 interpolation of T at u- & v-points (x2)
287	zC2t_v = ptn(ji,jj,jk,jn) + ptn(ji ,jj+1,jk,jn)
288	! ! C4 interpolation of T at u- & v-points (x2)
289	zC4t_u = zC2t_u + r1_6 * ( ztu(ji-1,jj ,jk) - ztu(ji+1,jj ,jk) )
290	zC4t_v = zC2t_v + r1_6 * ( ztv(ji ,jj-1,jk) - ztv(ji ,jj+1,jk) )
291	! ! C4 minus upstream advective fluxes
292	zwx(ji,jj,jk) = 0.5_wp * pun(ji,jj,jk) * zC4t_u - zwx(ji,jj,jk)
293	zwy(ji,jj,jk) = 0.5_wp * pvn(ji,jj,jk) * zC4t_v - zwy(ji,jj,jk)
294	END DO
295	END DO
296	END DO
297	!
298	END SELECT
299	!
300	SELECT CASE( kn_fct_v ) !* vertical anti-diffusive fluxes (w-masked interior values)
301	!
302	CASE( 2 ) !- 2nd order centered
303	DO jk = 2, jpkm1
304	DO jj = 2, jpjm1
305	DO ji = fs_2, fs_jpim1
306	zwz(ji,jj,jk) = ( pwn(ji,jj,jk) * 0.5_wp * ( ptn(ji,jj,jk,jn) + ptn(ji,jj,jk-1,jn) ) &
307	& - zwz(ji,jj,jk) ) * wmask(ji,jj,jk)
308	END DO
309	END DO
310	END DO
311	!
312	CASE( 4 ) !- 4th order COMPACT
313	CALL interp_4th_cpt( ptn(:,:,:,jn) , ztw ) ! zwt = COMPACT interpolation of T at w-point
314	DO jk = 2, jpkm1
315	DO jj = 2, jpjm1
316	DO ji = fs_2, fs_jpim1
317	zwz(ji,jj,jk) = ( pwn(ji,jj,jk) * ztw(ji,jj,jk) - zwz(ji,jj,jk) ) * wmask(ji,jj,jk)
318	END DO
319	END DO
320	END DO
321	!
322	END SELECT
323	IF( ln_linssh ) THEN ! top ocean value: high order = upstream ==>> zwz=0
324	zwz(:,:,1) = 0._wp ! only ocean surface as interior zwz values have been w-masked
325	ENDIF
326	!
327	IF ( ll_zAimp ) THEN
328	DO jk = 1, jpkm1 !* trend and after field with monotonic scheme
329	DO jj = 2, jpjm1
330	DO ji = fs_2, fs_jpim1 ! vector opt.
331	! ! total intermediate advective trends
332	ztra = - ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) &
333	& + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) &
334	& + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) * r1_e1e2t(ji,jj)
335	ztw(ji,jj,jk) = zwi(ji,jj,jk) + p2dt * ztra / e3t_a(ji,jj,jk) * tmask(ji,jj,jk)
336	END DO
337	END DO
338	END DO
339	!
340	CALL tridia_solver( zwdia, zwsup, zwinf, ztw, ztw , 0 )
341	!
342	DO jk = 2, jpkm1 ! Interior value ( multiplied by wmask)
343	DO jj = 2, jpjm1
344	DO ji = fs_2, fs_jpim1 ! vector opt.
345	zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) )
346	zfm_wk = wi(ji,jj,jk) - ABS( wi(ji,jj,jk) )
347	zwz(ji,jj,jk) = zwz(ji,jj,jk) + 0.5 * e1e2t(ji,jj) * ( zfp_wk * ztw(ji,jj,jk) + zfm_wk * ztw(ji,jj,jk-1) ) * wmask(ji,jj,jk)
348	END DO
349	END DO
350	END DO
351	END IF
352	!
353	CALL lbc_lnk_multi( 'traadv_fct', zwi, 'T', 1., zwx, 'U', -1. , zwy, 'V', -1., zwz, 'W', 1. )
354	!
355	! !== monotonicity algorithm ==!
356	!
357	CALL nonosc( ptb(:,:,:,jn), zwx, zwy, zwz, zwi, p2dt )
358	!
359	! !== final trend with corrected fluxes ==!
360	!
361	DO jk = 1, jpkm1
362	DO jj = 2, jpjm1
363	DO ji = fs_2, fs_jpim1 ! vector opt.
364	ztra = - ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) &
365	& + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) &
366	& + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) * r1_e1e2t(ji,jj)
367	pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra / e3t_n(ji,jj,jk)
368	zwi(ji,jj,jk) = zwi(ji,jj,jk) + p2dt * ztra / e3t_a(ji,jj,jk) * tmask(ji,jj,jk)
369	END DO
370	END DO
371	END DO
372	!
373	IF ( ll_zAimp ) THEN
374	!
375	ztw(:,:,1) = 0._wp ; ztw(:,:,jpk) = 0._wp
376	DO jk = 2, jpkm1 ! Interior value ( multiplied by wmask)
377	DO jj = 2, jpjm1
378	DO ji = fs_2, fs_jpim1 ! vector opt.
379	zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) )
380	zfm_wk = wi(ji,jj,jk) - ABS( wi(ji,jj,jk) )
381	ztw(ji,jj,jk) = - 0.5 * e1e2t(ji,jj) * ( zfp_wk * zwi(ji,jj,jk) + zfm_wk * zwi(ji,jj,jk-1) ) * wmask(ji,jj,jk)
382	zwz(ji,jj,jk) = zwz(ji,jj,jk) + ztw(ji,jj,jk) ! Update vertical fluxes for trend diagnostic
383	END DO
384	END DO
385	END DO
386	DO jk = 1, jpkm1
387	DO jj = 2, jpjm1
388	DO ji = fs_2, fs_jpim1 ! vector opt.
389	pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) - ( ztw(ji,jj,jk) - ztw(ji ,jj ,jk+1) ) &
390	& * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
391	END DO
392	END DO
393	END DO
394	END IF
395	!
396	IF( l_trd .OR. l_hst ) THEN ! trend diagnostics // heat/salt transport
397	ztrdx(:,:,:) = ztrdx(:,:,:) + zwx(:,:,:) ! <<< add anti-diffusive fluxes
398	ztrdy(:,:,:) = ztrdy(:,:,:) + zwy(:,:,:) ! to upstream fluxes
399	ztrdz(:,:,:) = ztrdz(:,:,:) + zwz(:,:,:) !
400	!
401	IF( l_trd ) THEN ! trend diagnostics
402	CALL trd_tra( kt, cdtype, jn, jptra_xad, ztrdx, pun, ptn(:,:,:,jn) )
403	CALL trd_tra( kt, cdtype, jn, jptra_yad, ztrdy, pvn, ptn(:,:,:,jn) )
404	CALL trd_tra( kt, cdtype, jn, jptra_zad, ztrdz, pwn, ptn(:,:,:,jn) )
405	ENDIF
406	! ! heat/salt transport
407	IF( l_hst ) CALL dia_ar5_hst( jn, 'adv', ztrdx(:,:,:), ztrdy(:,:,:) )
408	!
409	ENDIF
410	IF( l_ptr ) THEN ! "Poleward" transports
411	zptry(:,:,:) = zptry(:,:,:) + zwy(:,:,:) ! <<< add anti-diffusive fluxes
412	CALL dia_ptr_hst( jn, 'adv', zptry(:,:,:) )
413	ENDIF
414	!
415	END DO ! end of tracer loop
416	!
417	IF ( ll_zAimp ) THEN
418	DEALLOCATE( zwdia, zwinf, zwsup )
419	ENDIF
420	IF( l_trd .OR. l_hst ) THEN
421	DEALLOCATE( ztrdx, ztrdy, ztrdz )
422	ENDIF
423	IF( l_ptr ) THEN
424	DEALLOCATE( zptry )
425	ENDIF
426	!
427	END SUBROUTINE tra_adv_fct
428
429
430	SUBROUTINE nonosc( pbef, paa, pbb, pcc, paft, p2dt )
431	!!---------------------------------------------------------------------
432	!! * ROUTINE nonosc *
433	!!
434	!! ** Purpose : compute monotonic tracer fluxes from the upstream
435	!! scheme and the before field by a nonoscillatory algorithm
436	!!
437	!! ** Method : ... ???
438	!! warning : pbef and paft must be masked, but the boundaries
439	!! conditions on the fluxes are not necessary zalezak (1979)
440	!! drange (1995) multi-dimensional forward-in-time and upstream-
441	!! in-space based differencing for fluid
442	!!----------------------------------------------------------------------
443	REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step
444	REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(in ) :: pbef, paft ! before & after field
445	REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(inout) :: paa, pbb, pcc ! monotonic fluxes in the 3 directions
446	!
447	INTEGER :: ji, jj, jk ! dummy loop indices
448	INTEGER :: ikm1 ! local integer
449	REAL(wp) :: zpos, zneg, zbt, za, zb, zc, zbig, zrtrn ! local scalars
450	REAL(wp) :: zau, zbu, zcu, zav, zbv, zcv, zup, zdo ! - -
451	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zbetup, zbetdo, zbup, zbdo
452	!!----------------------------------------------------------------------
453	!
454	zbig = 1.e+40_wp
455	zrtrn = 1.e-15_wp
456	zbetup(:,:,:) = 0._wp ; zbetdo(:,:,:) = 0._wp
457
458	! Search local extrema
459	! --------------------
460	! max/min of pbef & paft with large negative/positive value (-/+zbig) inside land
461	zbup = MAX( pbef * tmask - zbig * ( 1._wp - tmask ), &
462	& paft * tmask - zbig * ( 1._wp - tmask ) )
463	zbdo = MIN( pbef * tmask + zbig * ( 1._wp - tmask ), &
464	& paft * tmask + zbig * ( 1._wp - tmask ) )
465
466	DO jk = 1, jpkm1
467	ikm1 = MAX(jk-1,1)
468	DO jj = 2, jpjm1
469	DO ji = fs_2, fs_jpim1 ! vector opt.
470
471	! search maximum in neighbourhood
472	zup = MAX( zbup(ji ,jj ,jk ), &
473	& zbup(ji-1,jj ,jk ), zbup(ji+1,jj ,jk ), &
474	& zbup(ji ,jj-1,jk ), zbup(ji ,jj+1,jk ), &
475	& zbup(ji ,jj ,ikm1), zbup(ji ,jj ,jk+1) )
476
477	! search minimum in neighbourhood
478	zdo = MIN( zbdo(ji ,jj ,jk ), &
479	& zbdo(ji-1,jj ,jk ), zbdo(ji+1,jj ,jk ), &
480	& zbdo(ji ,jj-1,jk ), zbdo(ji ,jj+1,jk ), &
481	& zbdo(ji ,jj ,ikm1), zbdo(ji ,jj ,jk+1) )
482
483	! positive part of the flux
484	zpos = MAX( 0., paa(ji-1,jj ,jk ) ) - MIN( 0., paa(ji ,jj ,jk ) ) &
485	& + MAX( 0., pbb(ji ,jj-1,jk ) ) - MIN( 0., pbb(ji ,jj ,jk ) ) &
486	& + MAX( 0., pcc(ji ,jj ,jk+1) ) - MIN( 0., pcc(ji ,jj ,jk ) )
487
488	! negative part of the flux
489	zneg = MAX( 0., paa(ji ,jj ,jk ) ) - MIN( 0., paa(ji-1,jj ,jk ) ) &
490	& + MAX( 0., pbb(ji ,jj ,jk ) ) - MIN( 0., pbb(ji ,jj-1,jk ) ) &
491	& + MAX( 0., pcc(ji ,jj ,jk ) ) - MIN( 0., pcc(ji ,jj ,jk+1) )
492
493	! up & down beta terms
494	zbt = e1e2t(ji,jj) * e3t_n(ji,jj,jk) / p2dt
495	zbetup(ji,jj,jk) = ( zup - paft(ji,jj,jk) ) / ( zpos + zrtrn ) * zbt
496	zbetdo(ji,jj,jk) = ( paft(ji,jj,jk) - zdo ) / ( zneg + zrtrn ) * zbt
497	END DO
498	END DO
499	END DO
500	CALL lbc_lnk_multi( 'traadv_fct', zbetup, 'T', 1. , zbetdo, 'T', 1. ) ! lateral boundary cond. (unchanged sign)
501
502	! 3. monotonic flux in the i & j direction (paa & pbb)
503	! ----------------------------------------
504	DO jk = 1, jpkm1
505	DO jj = 2, jpjm1
506	DO ji = fs_2, fs_jpim1 ! vector opt.
507	zau = MIN( 1._wp, zbetdo(ji,jj,jk), zbetup(ji+1,jj,jk) )
508	zbu = MIN( 1._wp, zbetup(ji,jj,jk), zbetdo(ji+1,jj,jk) )
509	zcu = ( 0.5 + SIGN( 0.5 , paa(ji,jj,jk) ) )
510	paa(ji,jj,jk) = paa(ji,jj,jk) * ( zcu * zau + ( 1._wp - zcu) * zbu )
511
512	zav = MIN( 1._wp, zbetdo(ji,jj,jk), zbetup(ji,jj+1,jk) )
513	zbv = MIN( 1._wp, zbetup(ji,jj,jk), zbetdo(ji,jj+1,jk) )
514	zcv = ( 0.5 + SIGN( 0.5 , pbb(ji,jj,jk) ) )
515	pbb(ji,jj,jk) = pbb(ji,jj,jk) * ( zcv * zav + ( 1._wp - zcv) * zbv )
516
517	! monotonic flux in the k direction, i.e. pcc
518	! -------------------------------------------
519	za = MIN( 1., zbetdo(ji,jj,jk+1), zbetup(ji,jj,jk) )
520	zb = MIN( 1., zbetup(ji,jj,jk+1), zbetdo(ji,jj,jk) )
521	zc = ( 0.5 + SIGN( 0.5 , pcc(ji,jj,jk+1) ) )
522	pcc(ji,jj,jk+1) = pcc(ji,jj,jk+1) * ( zc * za + ( 1._wp - zc) * zb )
523	END DO
524	END DO
525	END DO
526	CALL lbc_lnk_multi( 'traadv_fct', paa, 'U', -1. , pbb, 'V', -1. ) ! lateral boundary condition (changed sign)
527	!
528	END SUBROUTINE nonosc
529
530
531	SUBROUTINE interp_4th_cpt_org( pt_in, pt_out )
532	!!----------------------------------------------------------------------
533	!! * ROUTINE interp_4th_cpt_org *
534	!!
535	!! ** Purpose : Compute the interpolation of tracer at w-point
536	!!
537	!! ** Method : 4th order compact interpolation
538	!!----------------------------------------------------------------------
539	REAL(wp),DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pt_in ! now tracer fields
540	REAL(wp),DIMENSION(jpi,jpj,jpk), INTENT( out) :: pt_out ! now tracer field interpolated at w-pts
541	!
542	INTEGER :: ji, jj, jk ! dummy loop integers
543	REAL(wp),DIMENSION(jpi,jpj,jpk) :: zwd, zwi, zws, zwrm, zwt
544	!!----------------------------------------------------------------------
545
546	DO jk = 3, jpkm1 !== build the three diagonal matrix ==!
547	DO jj = 1, jpj
548	DO ji = 1, jpi
549	zwd (ji,jj,jk) = 4._wp
550	zwi (ji,jj,jk) = 1._wp
551	zws (ji,jj,jk) = 1._wp
552	zwrm(ji,jj,jk) = 3._wp * ( pt_in(ji,jj,jk-1) + pt_in(ji,jj,jk) )
553	!
554	IF( tmask(ji,jj,jk+1) == 0._wp) THEN ! Switch to second order centered at bottom
555	zwd (ji,jj,jk) = 1._wp
556	zwi (ji,jj,jk) = 0._wp
557	zws (ji,jj,jk) = 0._wp
558	zwrm(ji,jj,jk) = 0.5 * ( pt_in(ji,jj,jk-1) + pt_in(ji,jj,jk) )
559	ENDIF
560	END DO
561	END DO
562	END DO
563	!
564	jk = 2 ! Switch to second order centered at top
565	DO jj = 1, jpj
566	DO ji = 1, jpi
567	zwd (ji,jj,jk) = 1._wp
568	zwi (ji,jj,jk) = 0._wp
569	zws (ji,jj,jk) = 0._wp
570	zwrm(ji,jj,jk) = 0.5 * ( pt_in(ji,jj,jk-1) + pt_in(ji,jj,jk) )
571	END DO
572	END DO
573	!
574	! !== tridiagonal solve ==!
575	DO jj = 1, jpj ! first recurrence
576	DO ji = 1, jpi
577	zwt(ji,jj,2) = zwd(ji,jj,2)
578	END DO
579	END DO
580	DO jk = 3, jpkm1
581	DO jj = 1, jpj
582	DO ji = 1, jpi
583	zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) /zwt(ji,jj,jk-1)
584	END DO
585	END DO
586	END DO
587	!
588	DO jj = 1, jpj ! second recurrence: Zk = Yk - Ik / Tk-1 Zk-1
589	DO ji = 1, jpi
590	pt_out(ji,jj,2) = zwrm(ji,jj,2)
591	END DO
592	END DO
593	DO jk = 3, jpkm1
594	DO jj = 1, jpj
595	DO ji = 1, jpi
596	pt_out(ji,jj,jk) = zwrm(ji,jj,jk) - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1)
597	END DO
598	END DO
599	END DO
600
601	DO jj = 1, jpj ! third recurrence: Xk = (Zk - Sk Xk+1 ) / Tk
602	DO ji = 1, jpi
603	pt_out(ji,jj,jpkm1) = pt_out(ji,jj,jpkm1) / zwt(ji,jj,jpkm1)
604	END DO
605	END DO
606	DO jk = jpk-2, 2, -1
607	DO jj = 1, jpj
608	DO ji = 1, jpi
609	pt_out(ji,jj,jk) = ( pt_out(ji,jj,jk) - zws(ji,jj,jk) * pt_out(ji,jj,jk+1) ) / zwt(ji,jj,jk)
610	END DO
611	END DO
612	END DO
613	!
614	END SUBROUTINE interp_4th_cpt_org
615
616
617	SUBROUTINE interp_4th_cpt( pt_in, pt_out )
618	!!----------------------------------------------------------------------
619	!! * ROUTINE interp_4th_cpt *
620	!!
621	!! ** Purpose : Compute the interpolation of tracer at w-point
622	!!
623	!! ** Method : 4th order compact interpolation
624	!!----------------------------------------------------------------------
625	REAL(wp),DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pt_in ! field at t-point
626	REAL(wp),DIMENSION(jpi,jpj,jpk), INTENT( out) :: pt_out ! field interpolated at w-point
627	!
628	INTEGER :: ji, jj, jk ! dummy loop integers
629	INTEGER :: ikt, ikb ! local integers
630	REAL(wp),DIMENSION(jpi,jpj,jpk) :: zwd, zwi, zws, zwrm, zwt
631	!!----------------------------------------------------------------------
632	!
633	! !== build the three diagonal matrix & the RHS ==!
634	!
635	DO jk = 3, jpkm1 ! interior (from jk=3 to jpk-1)
636	DO jj = 2, jpjm1
637	DO ji = fs_2, fs_jpim1
638	zwd (ji,jj,jk) = 3._wp * wmask(ji,jj,jk) + 1._wp ! diagonal
639	zwi (ji,jj,jk) = wmask(ji,jj,jk) ! lower diagonal
640	zws (ji,jj,jk) = wmask(ji,jj,jk) ! upper diagonal
641	zwrm(ji,jj,jk) = 3._wp * wmask(ji,jj,jk) & ! RHS
642	& * ( pt_in(ji,jj,jk) + pt_in(ji,jj,jk-1) )
643	END DO
644	END DO
645	END DO
646	!
647	!!gm
648	! SELECT CASE( kbc ) !* boundary condition
649	! CASE( np_NH ) ! Neumann homogeneous at top & bottom
650	! CASE( np_CEN2 ) ! 2nd order centered at top & bottom
651	! END SELECT
652	!!gm
653	!
654	IF ( ln_isfcav ) THEN ! set level two values which may not be set in ISF case
655	zwd(:,:,2) = 1._wp ; zwi(:,:,2) = 0._wp ; zws(:,:,2) = 0._wp ; zwrm(:,:,2) = 0._wp
656	END IF
657	!
658	DO jj = 2, jpjm1 ! 2nd order centered at top & bottom
659	DO ji = fs_2, fs_jpim1
660	ikt = mikt(ji,jj) + 1 ! w-point below the 1st wet point
661	ikb = MAX(mbkt(ji,jj), 2) ! - above the last wet point
662	!
663	zwd (ji,jj,ikt) = 1._wp ! top
664	zwi (ji,jj,ikt) = 0._wp
665	zws (ji,jj,ikt) = 0._wp
666	zwrm(ji,jj,ikt) = 0.5_wp * ( pt_in(ji,jj,ikt-1) + pt_in(ji,jj,ikt) )
667	!
668	zwd (ji,jj,ikb) = 1._wp ! bottom
669	zwi (ji,jj,ikb) = 0._wp
670	zws (ji,jj,ikb) = 0._wp
671	zwrm(ji,jj,ikb) = 0.5_wp * ( pt_in(ji,jj,ikb-1) + pt_in(ji,jj,ikb) )
672	END DO
673	END DO
674	!
675	! !== tridiagonal solver ==!
676	!
677	DO jj = 2, jpjm1 !* 1st recurrence: Tk = Dk - Ik Sk-1 / Tk-1
678	DO ji = fs_2, fs_jpim1
679	zwt(ji,jj,2) = zwd(ji,jj,2)
680	END DO
681	END DO
682	DO jk = 3, jpkm1
683	DO jj = 2, jpjm1
684	DO ji = fs_2, fs_jpim1
685	zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) /zwt(ji,jj,jk-1)
686	END DO
687	END DO
688	END DO
689	!
690	DO jj = 2, jpjm1 !* 2nd recurrence: Zk = Yk - Ik / Tk-1 Zk-1
691	DO ji = fs_2, fs_jpim1
692	pt_out(ji,jj,2) = zwrm(ji,jj,2)
693	END DO
694	END DO
695	DO jk = 3, jpkm1
696	DO jj = 2, jpjm1
697	DO ji = fs_2, fs_jpim1
698	pt_out(ji,jj,jk) = zwrm(ji,jj,jk) - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1)
699	END DO
700	END DO
701	END DO
702
703	DO jj = 2, jpjm1 !* 3d recurrence: Xk = (Zk - Sk Xk+1 ) / Tk
704	DO ji = fs_2, fs_jpim1
705	pt_out(ji,jj,jpkm1) = pt_out(ji,jj,jpkm1) / zwt(ji,jj,jpkm1)
706	END DO
707	END DO
708	DO jk = jpk-2, 2, -1
709	DO jj = 2, jpjm1
710	DO ji = fs_2, fs_jpim1
711	pt_out(ji,jj,jk) = ( pt_out(ji,jj,jk) - zws(ji,jj,jk) * pt_out(ji,jj,jk+1) ) / zwt(ji,jj,jk)
712	END DO
713	END DO
714	END DO
715	!
716	END SUBROUTINE interp_4th_cpt
717
718
719	SUBROUTINE tridia_solver( pD, pU, pL, pRHS, pt_out , klev )
720	!!----------------------------------------------------------------------
721	!! * ROUTINE tridia_solver *
722	!!
723	!! ** Purpose : solve a symmetric 3diagonal system
724	!!
725	!! ** Method : solve M.t_out = RHS(t) where M is a tri diagonal matrix ( jpk*jpk )
726	!!
727	!! ( D_1 U_1 0 0 0 )( t_1 ) ( RHS_1 )
728	!! ( L_2 D_2 U_2 0 0 )( t_2 ) ( RHS_2 )
729	!! ( 0 L_3 D_3 U_3 0 )( t_3 ) = ( RHS_3 )
730	!! ( ... )( ... ) ( ... )
731	!! ( 0 0 0 L_k D_k )( t_k ) ( RHS_k )
732	!!
733	!! M is decomposed in the product of an upper and lower triangular matrix.
734	!! The tri-diagonals matrix is given as input 3D arrays: pD, pU, pL
735	!! (i.e. the Diagonal, the Upper diagonal, and the Lower diagonal).
736	!! The solution is pta.
737	!! The 3d array zwt is used as a work space array.
738	!!----------------------------------------------------------------------
739	REAL(wp),DIMENSION(:,:,:), INTENT(in ) :: pD, pU, PL ! 3-diagonal matrix
740	REAL(wp),DIMENSION(:,:,:), INTENT(in ) :: pRHS ! Right-Hand-Side
741	REAL(wp),DIMENSION(:,:,:), INTENT( out) :: pt_out !!gm field at level=F(klev)
742	INTEGER , INTENT(in ) :: klev ! =1 pt_out at w-level
743	! ! =0 pt at t-level
744	INTEGER :: ji, jj, jk ! dummy loop integers
745	INTEGER :: kstart ! local indices
746	REAL(wp),DIMENSION(jpi,jpj,jpk) :: zwt ! 3D work array
747	!!----------------------------------------------------------------------
748	!
749	kstart = 1 + klev
750	!
751	DO jj = 2, jpjm1 !* 1st recurrence: Tk = Dk - Ik Sk-1 / Tk-1
752	DO ji = fs_2, fs_jpim1
753	zwt(ji,jj,kstart) = pD(ji,jj,kstart)
754	END DO
755	END DO
756	DO jk = kstart+1, jpkm1
757	DO jj = 2, jpjm1
758	DO ji = fs_2, fs_jpim1
759	zwt(ji,jj,jk) = pD(ji,jj,jk) - pL(ji,jj,jk) * pU(ji,jj,jk-1) /zwt(ji,jj,jk-1)
760	END DO
761	END DO
762	END DO
763	!
764	DO jj = 2, jpjm1 !* 2nd recurrence: Zk = Yk - Ik / Tk-1 Zk-1
765	DO ji = fs_2, fs_jpim1
766	pt_out(ji,jj,kstart) = pRHS(ji,jj,kstart)
767	END DO
768	END DO
769	DO jk = kstart+1, jpkm1
770	DO jj = 2, jpjm1
771	DO ji = fs_2, fs_jpim1
772	pt_out(ji,jj,jk) = pRHS(ji,jj,jk) - pL(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1)
773	END DO
774	END DO
775	END DO
776
777	DO jj = 2, jpjm1 !* 3d recurrence: Xk = (Zk - Sk Xk+1 ) / Tk
778	DO ji = fs_2, fs_jpim1
779	pt_out(ji,jj,jpkm1) = pt_out(ji,jj,jpkm1) / zwt(ji,jj,jpkm1)
780	END DO
781	END DO
782	DO jk = jpk-2, kstart, -1
783	DO jj = 2, jpjm1
784	DO ji = fs_2, fs_jpim1
785	pt_out(ji,jj,jk) = ( pt_out(ji,jj,jk) - pU(ji,jj,jk) * pt_out(ji,jj,jk+1) ) / zwt(ji,jj,jk)
786	END DO
787	END DO
788	END DO
789	!
790	END SUBROUTINE tridia_solver
791
792	!!======================================================================
793	END MODULE traadv_fct

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.1_GO8_package/src/OCE/TRA/traadv_fct.F90 @ 12088

Download in other formats: