Context Navigation

traadv_fct.F90 @ 12396

Last change on this file since 12396 was 12377, checked in by acc, 4 years ago

The big one. Merging all 2019 developments from the option 1 branch back onto the trunk.

This changeset reproduces 2019/dev_r11943_MERGE_2019 on the trunk using a 2-URL merge
onto a working copy of the trunk. I.e.:

svn merge --ignore-ancestry \

svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/trunk \
svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/branches/2019/dev_r11943_MERGE_2019 ./

The --ignore-ancestry flag avoids problems that may otherwise arise from the fact that
the merge history been trunk and branch may have been applied in a different order but
care has been taken before this step to ensure that all applicable fixes and updates
are present in the merge branch.

The trunk state just before this step has been branched to releases/release-4.0-HEAD
and that branch has been immediately tagged as releases/release-4.0.2. Any fixes
or additions in response to tickets on 4.0, 4.0.1 or 4.0.2 should be done on
releases/release-4.0-HEAD. From now on future 'point' releases (e.g. 4.0.2) will
remain unchanged with periodic releases as needs demand. Note release-4.0-HEAD is a
transitional naming convention. Future full releases, say 4.2, will have a release-4.2
branch which fulfills this role and the first point release (e.g. 4.2.0) will be made
immediately following the release branch creation.

2020 developments can be started from any trunk revision later than this one.

Property svn:keywords set to Id

File size: 32.5 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 "do_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, pU, pV, pW, &
56	& Kbb, Kmm, pt, kjpt, Krhs, 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 pt(:,:,:,:,Krhs) with the now advective tracer trends
69	!! - send trends to trdtra module for further diagnostics (l_trdtra=T)
70	!! - poleward advective heat and salt transport (ln_diaptr=T)
71	!!----------------------------------------------------------------------
72	INTEGER , INTENT(in ) :: kt ! ocean time-step index
73	INTEGER , INTENT(in ) :: Kbb, Kmm, Krhs ! ocean time level indices
74	INTEGER , INTENT(in ) :: kit000 ! first time step index
75	CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator)
76	INTEGER , INTENT(in ) :: kjpt ! number of tracers
77	INTEGER , INTENT(in ) :: kn_fct_h ! order of the FCT scheme (=2 or 4)
78	INTEGER , INTENT(in ) :: kn_fct_v ! order of the FCT scheme (=2 or 4)
79	REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step
80	REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pU, pV, pW ! 3 ocean volume flux components
81	REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) :: pt ! tracers and RHS of tracer equation
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. ( iom_use( 'sophtadv' ) .OR. iom_use( 'sophtadv' ) ) ) 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_3D_00_00( 1, jpkm1 )
131	zwdia(ji,jj,jk) = 1._wp + p2dt * ( MAX( wi(ji,jj,jk ) , 0._wp ) - MIN( wi(ji,jj,jk+1) , 0._wp ) ) / e3t(ji,jj,jk,Krhs)
132	zwinf(ji,jj,jk) = p2dt * MIN( wi(ji,jj,jk ) , 0._wp ) / e3t(ji,jj,jk,Krhs)
133	zwsup(ji,jj,jk) = -p2dt * MAX( wi(ji,jj,jk+1) , 0._wp ) / e3t(ji,jj,jk,Krhs)
134	END_3D
135	END IF
136	!
137	DO jn = 1, kjpt !== loop over the tracers ==!
138	!
139	! !== upstream advection with initial mass fluxes & intermediate update ==!
140	! !* upstream tracer flux in the i and j direction
141	DO_3D_10_10( 1, jpkm1 )
142	! upstream scheme
143	zfp_ui = pU(ji,jj,jk) + ABS( pU(ji,jj,jk) )
144	zfm_ui = pU(ji,jj,jk) - ABS( pU(ji,jj,jk) )
145	zfp_vj = pV(ji,jj,jk) + ABS( pV(ji,jj,jk) )
146	zfm_vj = pV(ji,jj,jk) - ABS( pV(ji,jj,jk) )
147	zwx(ji,jj,jk) = 0.5 * ( zfp_ui * pt(ji,jj,jk,jn,Kbb) + zfm_ui * pt(ji+1,jj ,jk,jn,Kbb) )
148	zwy(ji,jj,jk) = 0.5 * ( zfp_vj * pt(ji,jj,jk,jn,Kbb) + zfm_vj * pt(ji ,jj+1,jk,jn,Kbb) )
149	END_3D
150	! !* upstream tracer flux in the k direction *!
151	DO_3D_11_11( 2, jpkm1 )
152	zfp_wk = pW(ji,jj,jk) + ABS( pW(ji,jj,jk) )
153	zfm_wk = pW(ji,jj,jk) - ABS( pW(ji,jj,jk) )
154	zwz(ji,jj,jk) = 0.5 * ( zfp_wk * pt(ji,jj,jk,jn,Kbb) + zfm_wk * pt(ji,jj,jk-1,jn,Kbb) ) * wmask(ji,jj,jk)
155	END_3D
156	IF( ln_linssh ) THEN ! top ocean value (only in linear free surface as zwz has been w-masked)
157	IF( ln_isfcav ) THEN ! top of the ice-shelf cavities and at the ocean surface
158	DO_2D_11_11
159	zwz(ji,jj, mikt(ji,jj) ) = pW(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kbb) ! linear free surface
160	END_2D
161	ELSE ! no cavities: only at the ocean surface
162	zwz(:,:,1) = pW(:,:,1) * pt(:,:,1,jn,Kbb)
163	ENDIF
164	ENDIF
165	!
166	DO_3D_00_00( 1, jpkm1 )
167	! ! total intermediate advective trends
168	ztra = - ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) &
169	& + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) &
170	& + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) * r1_e1e2t(ji,jj)
171	! ! update and guess with monotonic sheme
172	pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) + ztra / e3t(ji,jj,jk,Kmm) * tmask(ji,jj,jk)
173	zwi(ji,jj,jk) = ( e3t(ji,jj,jk,Kbb) * pt(ji,jj,jk,jn,Kbb) + p2dt * ztra ) / e3t(ji,jj,jk,Krhs) * tmask(ji,jj,jk)
174	END_3D
175
176	IF ( ll_zAimp ) THEN
177	CALL tridia_solver( zwdia, zwsup, zwinf, zwi, zwi , 0 )
178	!
179	ztw(:,:,1) = 0._wp ; ztw(:,:,jpk) = 0._wp ;
180	DO_3D_00_00( 2, jpkm1 )
181	zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) )
182	zfm_wk = wi(ji,jj,jk) - ABS( wi(ji,jj,jk) )
183	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)
184	zwz(ji,jj,jk) = zwz(ji,jj,jk) + ztw(ji,jj,jk) ! update vertical fluxes
185	END_3D
186	DO_3D_00_00( 1, jpkm1 )
187	pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( ztw(ji,jj,jk) - ztw(ji ,jj ,jk+1) ) &
188	& * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
189	END_3D
190	!
191	END IF
192	!
193	IF( l_trd .OR. l_hst ) THEN ! trend diagnostics (contribution of upstream fluxes)
194	ztrdx(:,:,:) = zwx(:,:,:) ; ztrdy(:,:,:) = zwy(:,:,:) ; ztrdz(:,:,:) = zwz(:,:,:)
195	END IF
196	! ! "Poleward" heat and salt transports (contribution of upstream fluxes)
197	IF( l_ptr ) zptry(:,:,:) = zwy(:,:,:)
198	!
199	! !== anti-diffusive flux : high order minus low order ==!
200	!
201	SELECT CASE( kn_fct_h ) !* horizontal anti-diffusive fluxes
202	!
203	CASE( 2 ) !- 2nd order centered
204	DO_3D_10_10( 1, jpkm1 )
205	zwx(ji,jj,jk) = 0.5_wp * pU(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji+1,jj,jk,jn,Kmm) ) - zwx(ji,jj,jk)
206	zwy(ji,jj,jk) = 0.5_wp * pV(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji,jj+1,jk,jn,Kmm) ) - zwy(ji,jj,jk)
207	END_3D
208	!
209	CASE( 4 ) !- 4th order centered
210	zltu(:,:,jpk) = 0._wp ! Bottom value : flux set to zero
211	zltv(:,:,jpk) = 0._wp
212	DO jk = 1, jpkm1 ! Laplacian
213	DO_2D_10_10
214	ztu(ji,jj,jk) = ( pt(ji+1,jj ,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * umask(ji,jj,jk)
215	ztv(ji,jj,jk) = ( pt(ji ,jj+1,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * vmask(ji,jj,jk)
216	END_2D
217	DO_2D_00_00
218	zltu(ji,jj,jk) = ( ztu(ji,jj,jk) + ztu(ji-1,jj,jk) ) * r1_6
219	zltv(ji,jj,jk) = ( ztv(ji,jj,jk) + ztv(ji,jj-1,jk) ) * r1_6
220	END_2D
221	END DO
222	CALL lbc_lnk_multi( 'traadv_fct', zltu, 'T', 1. , zltv, 'T', 1. ) ! Lateral boundary cond. (unchanged sgn)
223	!
224	DO_3D_10_10( 1, jpkm1 )
225	zC2t_u = pt(ji,jj,jk,jn,Kmm) + pt(ji+1,jj ,jk,jn,Kmm) ! 2 x C2 interpolation of T at u- & v-points
226	zC2t_v = pt(ji,jj,jk,jn,Kmm) + pt(ji ,jj+1,jk,jn,Kmm)
227	! ! C4 minus upstream advective fluxes
228	zwx(ji,jj,jk) = 0.5_wp * pU(ji,jj,jk) * ( zC2t_u + zltu(ji,jj,jk) - zltu(ji+1,jj,jk) ) - zwx(ji,jj,jk)
229	zwy(ji,jj,jk) = 0.5_wp * pV(ji,jj,jk) * ( zC2t_v + zltv(ji,jj,jk) - zltv(ji,jj+1,jk) ) - zwy(ji,jj,jk)
230	END_3D
231	!
232	CASE( 41 ) !- 4th order centered ==>> !!gm coding attempt need to be tested
233	ztu(:,:,jpk) = 0._wp ! Bottom value : flux set to zero
234	ztv(:,:,jpk) = 0._wp
235	DO_3D_10_10( 1, jpkm1 )
236	ztu(ji,jj,jk) = ( pt(ji+1,jj ,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * umask(ji,jj,jk)
237	ztv(ji,jj,jk) = ( pt(ji ,jj+1,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * vmask(ji,jj,jk)
238	END_3D
239	CALL lbc_lnk_multi( 'traadv_fct', ztu, 'U', -1. , ztv, 'V', -1. ) ! Lateral boundary cond. (unchanged sgn)
240	!
241	DO_3D_00_00( 1, jpkm1 )
242	zC2t_u = pt(ji,jj,jk,jn,Kmm) + pt(ji+1,jj ,jk,jn,Kmm) ! 2 x C2 interpolation of T at u- & v-points (x2)
243	zC2t_v = pt(ji,jj,jk,jn,Kmm) + pt(ji ,jj+1,jk,jn,Kmm)
244	! ! C4 interpolation of T at u- & v-points (x2)
245	zC4t_u = zC2t_u + r1_6 * ( ztu(ji-1,jj ,jk) - ztu(ji+1,jj ,jk) )
246	zC4t_v = zC2t_v + r1_6 * ( ztv(ji ,jj-1,jk) - ztv(ji ,jj+1,jk) )
247	! ! C4 minus upstream advective fluxes
248	zwx(ji,jj,jk) = 0.5_wp * pU(ji,jj,jk) * zC4t_u - zwx(ji,jj,jk)
249	zwy(ji,jj,jk) = 0.5_wp * pV(ji,jj,jk) * zC4t_v - zwy(ji,jj,jk)
250	END_3D
251	!
252	END SELECT
253	!
254	SELECT CASE( kn_fct_v ) !* vertical anti-diffusive fluxes (w-masked interior values)
255	!
256	CASE( 2 ) !- 2nd order centered
257	DO_3D_00_00( 2, jpkm1 )
258	zwz(ji,jj,jk) = ( pW(ji,jj,jk) * 0.5_wp * ( pt(ji,jj,jk,jn,Kmm) + pt(ji,jj,jk-1,jn,Kmm) ) &
259	& - zwz(ji,jj,jk) ) * wmask(ji,jj,jk)
260	END_3D
261	!
262	CASE( 4 ) !- 4th order COMPACT
263	CALL interp_4th_cpt( pt(:,:,:,jn,Kmm) , ztw ) ! zwt = COMPACT interpolation of T at w-point
264	DO_3D_00_00( 2, jpkm1 )
265	zwz(ji,jj,jk) = ( pW(ji,jj,jk) * ztw(ji,jj,jk) - zwz(ji,jj,jk) ) * wmask(ji,jj,jk)
266	END_3D
267	!
268	END SELECT
269	IF( ln_linssh ) THEN ! top ocean value: high order = upstream ==>> zwz=0
270	zwz(:,:,1) = 0._wp ! only ocean surface as interior zwz values have been w-masked
271	ENDIF
272	!
273	IF ( ll_zAimp ) THEN
274	DO_3D_00_00( 1, jpkm1 )
275	! ! total intermediate advective trends
276	ztra = - ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) &
277	& + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) &
278	& + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) * r1_e1e2t(ji,jj)
279	ztw(ji,jj,jk) = zwi(ji,jj,jk) + p2dt * ztra / e3t(ji,jj,jk,Krhs) * tmask(ji,jj,jk)
280	END_3D
281	!
282	CALL tridia_solver( zwdia, zwsup, zwinf, ztw, ztw , 0 )
283	!
284	DO_3D_00_00( 2, jpkm1 )
285	zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) )
286	zfm_wk = wi(ji,jj,jk) - ABS( wi(ji,jj,jk) )
287	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)
288	END_3D
289	END IF
290	!
291	CALL lbc_lnk_multi( 'traadv_fct', zwi, 'T', 1., zwx, 'U', -1. , zwy, 'V', -1., zwz, 'W', 1. )
292	!
293	! !== monotonicity algorithm ==!
294	!
295	CALL nonosc( Kmm, pt(:,:,:,jn,Kbb), zwx, zwy, zwz, zwi, p2dt )
296	!
297	! !== final trend with corrected fluxes ==!
298	!
299	DO_3D_00_00( 1, jpkm1 )
300	ztra = - ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) &
301	& + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) &
302	& + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) * r1_e1e2t(ji,jj)
303	pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) + ztra / e3t(ji,jj,jk,Kmm)
304	zwi(ji,jj,jk) = zwi(ji,jj,jk) + p2dt * ztra / e3t(ji,jj,jk,Krhs) * tmask(ji,jj,jk)
305	END_3D
306	!
307	IF ( ll_zAimp ) THEN
308	!
309	ztw(:,:,1) = 0._wp ; ztw(:,:,jpk) = 0._wp
310	DO_3D_00_00( 2, jpkm1 )
311	zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) )
312	zfm_wk = wi(ji,jj,jk) - ABS( wi(ji,jj,jk) )
313	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)
314	zwz(ji,jj,jk) = zwz(ji,jj,jk) + ztw(ji,jj,jk) ! Update vertical fluxes for trend diagnostic
315	END_3D
316	DO_3D_00_00( 1, jpkm1 )
317	pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( ztw(ji,jj,jk) - ztw(ji ,jj ,jk+1) ) &
318	& * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
319	END_3D
320	END IF
321	!
322	IF( l_trd .OR. l_hst ) THEN ! trend diagnostics // heat/salt transport
323	ztrdx(:,:,:) = ztrdx(:,:,:) + zwx(:,:,:) ! <<< add anti-diffusive fluxes
324	ztrdy(:,:,:) = ztrdy(:,:,:) + zwy(:,:,:) ! to upstream fluxes
325	ztrdz(:,:,:) = ztrdz(:,:,:) + zwz(:,:,:) !
326	!
327	IF( l_trd ) THEN ! trend diagnostics
328	CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_xad, ztrdx, pU, pt(:,:,:,jn,Kmm) )
329	CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_yad, ztrdy, pV, pt(:,:,:,jn,Kmm) )
330	CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_zad, ztrdz, pW, pt(:,:,:,jn,Kmm) )
331	ENDIF
332	! ! heat/salt transport
333	IF( l_hst ) CALL dia_ar5_hst( jn, 'adv', ztrdx(:,:,:), ztrdy(:,:,:) )
334	!
335	ENDIF
336	IF( l_ptr ) THEN ! "Poleward" transports
337	zptry(:,:,:) = zptry(:,:,:) + zwy(:,:,:) ! <<< add anti-diffusive fluxes
338	CALL dia_ptr_hst( jn, 'adv', zptry(:,:,:) )
339	ENDIF
340	!
341	END DO ! end of tracer loop
342	!
343	IF ( ll_zAimp ) THEN
344	DEALLOCATE( zwdia, zwinf, zwsup )
345	ENDIF
346	IF( l_trd .OR. l_hst ) THEN
347	DEALLOCATE( ztrdx, ztrdy, ztrdz )
348	ENDIF
349	IF( l_ptr ) THEN
350	DEALLOCATE( zptry )
351	ENDIF
352	!
353	END SUBROUTINE tra_adv_fct
354
355
356	SUBROUTINE nonosc( Kmm, pbef, paa, pbb, pcc, paft, p2dt )
357	!!---------------------------------------------------------------------
358	!! * ROUTINE nonosc *
359	!!
360	!! ** Purpose : compute monotonic tracer fluxes from the upstream
361	!! scheme and the before field by a nonoscillatory algorithm
362	!!
363	!! ** Method : ... ???
364	!! warning : pbef and paft must be masked, but the boundaries
365	!! conditions on the fluxes are not necessary zalezak (1979)
366	!! drange (1995) multi-dimensional forward-in-time and upstream-
367	!! in-space based differencing for fluid
368	!!----------------------------------------------------------------------
369	INTEGER , INTENT(in ) :: Kmm ! time level index
370	REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step
371	REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(in ) :: pbef, paft ! before & after field
372	REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(inout) :: paa, pbb, pcc ! monotonic fluxes in the 3 directions
373	!
374	INTEGER :: ji, jj, jk ! dummy loop indices
375	INTEGER :: ikm1 ! local integer
376	REAL(wp) :: zpos, zneg, zbt, za, zb, zc, zbig, zrtrn ! local scalars
377	REAL(wp) :: zau, zbu, zcu, zav, zbv, zcv, zup, zdo ! - -
378	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zbetup, zbetdo, zbup, zbdo
379	!!----------------------------------------------------------------------
380	!
381	zbig = 1.e+40_wp
382	zrtrn = 1.e-15_wp
383	zbetup(:,:,:) = 0._wp ; zbetdo(:,:,:) = 0._wp
384
385	! Search local extrema
386	! --------------------
387	! max/min of pbef & paft with large negative/positive value (-/+zbig) inside land
388	zbup = MAX( pbef * tmask - zbig * ( 1._wp - tmask ), &
389	& paft * tmask - zbig * ( 1._wp - tmask ) )
390	zbdo = MIN( pbef * tmask + zbig * ( 1._wp - tmask ), &
391	& paft * tmask + zbig * ( 1._wp - tmask ) )
392
393	DO jk = 1, jpkm1
394	ikm1 = MAX(jk-1,1)
395	DO_2D_00_00
396
397	! search maximum in neighbourhood
398	zup = MAX( zbup(ji ,jj ,jk ), &
399	& zbup(ji-1,jj ,jk ), zbup(ji+1,jj ,jk ), &
400	& zbup(ji ,jj-1,jk ), zbup(ji ,jj+1,jk ), &
401	& zbup(ji ,jj ,ikm1), zbup(ji ,jj ,jk+1) )
402
403	! search minimum in neighbourhood
404	zdo = MIN( zbdo(ji ,jj ,jk ), &
405	& zbdo(ji-1,jj ,jk ), zbdo(ji+1,jj ,jk ), &
406	& zbdo(ji ,jj-1,jk ), zbdo(ji ,jj+1,jk ), &
407	& zbdo(ji ,jj ,ikm1), zbdo(ji ,jj ,jk+1) )
408
409	! positive part of the flux
410	zpos = MAX( 0., paa(ji-1,jj ,jk ) ) - MIN( 0., paa(ji ,jj ,jk ) ) &
411	& + MAX( 0., pbb(ji ,jj-1,jk ) ) - MIN( 0., pbb(ji ,jj ,jk ) ) &
412	& + MAX( 0., pcc(ji ,jj ,jk+1) ) - MIN( 0., pcc(ji ,jj ,jk ) )
413
414	! negative part of the flux
415	zneg = MAX( 0., paa(ji ,jj ,jk ) ) - MIN( 0., paa(ji-1,jj ,jk ) ) &
416	& + MAX( 0., pbb(ji ,jj ,jk ) ) - MIN( 0., pbb(ji ,jj-1,jk ) ) &
417	& + MAX( 0., pcc(ji ,jj ,jk ) ) - MIN( 0., pcc(ji ,jj ,jk+1) )
418
419	! up & down beta terms
420	zbt = e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm) / p2dt
421	zbetup(ji,jj,jk) = ( zup - paft(ji,jj,jk) ) / ( zpos + zrtrn ) * zbt
422	zbetdo(ji,jj,jk) = ( paft(ji,jj,jk) - zdo ) / ( zneg + zrtrn ) * zbt
423	END_2D
424	END DO
425	CALL lbc_lnk_multi( 'traadv_fct', zbetup, 'T', 1. , zbetdo, 'T', 1. ) ! lateral boundary cond. (unchanged sign)
426
427	! 3. monotonic flux in the i & j direction (paa & pbb)
428	! ----------------------------------------
429	DO_3D_00_00( 1, jpkm1 )
430	zau = MIN( 1._wp, zbetdo(ji,jj,jk), zbetup(ji+1,jj,jk) )
431	zbu = MIN( 1._wp, zbetup(ji,jj,jk), zbetdo(ji+1,jj,jk) )
432	zcu = ( 0.5 + SIGN( 0.5 , paa(ji,jj,jk) ) )
433	paa(ji,jj,jk) = paa(ji,jj,jk) * ( zcu * zau + ( 1._wp - zcu) * zbu )
434
435	zav = MIN( 1._wp, zbetdo(ji,jj,jk), zbetup(ji,jj+1,jk) )
436	zbv = MIN( 1._wp, zbetup(ji,jj,jk), zbetdo(ji,jj+1,jk) )
437	zcv = ( 0.5 + SIGN( 0.5 , pbb(ji,jj,jk) ) )
438	pbb(ji,jj,jk) = pbb(ji,jj,jk) * ( zcv * zav + ( 1._wp - zcv) * zbv )
439
440	! monotonic flux in the k direction, i.e. pcc
441	! -------------------------------------------
442	za = MIN( 1., zbetdo(ji,jj,jk+1), zbetup(ji,jj,jk) )
443	zb = MIN( 1., zbetup(ji,jj,jk+1), zbetdo(ji,jj,jk) )
444	zc = ( 0.5 + SIGN( 0.5 , pcc(ji,jj,jk+1) ) )
445	pcc(ji,jj,jk+1) = pcc(ji,jj,jk+1) * ( zc * za + ( 1._wp - zc) * zb )
446	END_3D
447	CALL lbc_lnk_multi( 'traadv_fct', paa, 'U', -1. , pbb, 'V', -1. ) ! lateral boundary condition (changed sign)
448	!
449	END SUBROUTINE nonosc
450
451
452	SUBROUTINE interp_4th_cpt_org( pt_in, pt_out )
453	!!----------------------------------------------------------------------
454	!! * ROUTINE interp_4th_cpt_org *
455	!!
456	!! ** Purpose : Compute the interpolation of tracer at w-point
457	!!
458	!! ** Method : 4th order compact interpolation
459	!!----------------------------------------------------------------------
460	REAL(wp),DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pt_in ! now tracer fields
461	REAL(wp),DIMENSION(jpi,jpj,jpk), INTENT( out) :: pt_out ! now tracer field interpolated at w-pts
462	!
463	INTEGER :: ji, jj, jk ! dummy loop integers
464	REAL(wp),DIMENSION(jpi,jpj,jpk) :: zwd, zwi, zws, zwrm, zwt
465	!!----------------------------------------------------------------------
466
467	DO_3D_11_11( 3, jpkm1 )
468	zwd (ji,jj,jk) = 4._wp
469	zwi (ji,jj,jk) = 1._wp
470	zws (ji,jj,jk) = 1._wp
471	zwrm(ji,jj,jk) = 3._wp * ( pt_in(ji,jj,jk-1) + pt_in(ji,jj,jk) )
472	!
473	IF( tmask(ji,jj,jk+1) == 0._wp) THEN ! Switch to second order centered at bottom
474	zwd (ji,jj,jk) = 1._wp
475	zwi (ji,jj,jk) = 0._wp
476	zws (ji,jj,jk) = 0._wp
477	zwrm(ji,jj,jk) = 0.5 * ( pt_in(ji,jj,jk-1) + pt_in(ji,jj,jk) )
478	ENDIF
479	END_3D
480	!
481	jk = 2 ! Switch to second order centered at top
482	DO_2D_11_11
483	zwd (ji,jj,jk) = 1._wp
484	zwi (ji,jj,jk) = 0._wp
485	zws (ji,jj,jk) = 0._wp
486	zwrm(ji,jj,jk) = 0.5 * ( pt_in(ji,jj,jk-1) + pt_in(ji,jj,jk) )
487	END_2D
488	!
489	! !== tridiagonal solve ==!
490	DO_2D_11_11
491	zwt(ji,jj,2) = zwd(ji,jj,2)
492	END_2D
493	DO_3D_11_11( 3, jpkm1 )
494	zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) /zwt(ji,jj,jk-1)
495	END_3D
496	!
497	DO_2D_11_11
498	pt_out(ji,jj,2) = zwrm(ji,jj,2)
499	END_2D
500	DO_3D_11_11( 3, jpkm1 )
501	pt_out(ji,jj,jk) = zwrm(ji,jj,jk) - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1)
502	END_3D
503
504	DO_2D_11_11
505	pt_out(ji,jj,jpkm1) = pt_out(ji,jj,jpkm1) / zwt(ji,jj,jpkm1)
506	END_2D
507	DO_3DS_11_11( jpk-2, 2, -1 )
508	pt_out(ji,jj,jk) = ( pt_out(ji,jj,jk) - zws(ji,jj,jk) * pt_out(ji,jj,jk+1) ) / zwt(ji,jj,jk)
509	END_3D
510	!
511	END SUBROUTINE interp_4th_cpt_org
512
513
514	SUBROUTINE interp_4th_cpt( pt_in, pt_out )
515	!!----------------------------------------------------------------------
516	!! * ROUTINE interp_4th_cpt *
517	!!
518	!! ** Purpose : Compute the interpolation of tracer at w-point
519	!!
520	!! ** Method : 4th order compact interpolation
521	!!----------------------------------------------------------------------
522	REAL(wp),DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pt_in ! field at t-point
523	REAL(wp),DIMENSION(jpi,jpj,jpk), INTENT( out) :: pt_out ! field interpolated at w-point
524	!
525	INTEGER :: ji, jj, jk ! dummy loop integers
526	INTEGER :: ikt, ikb ! local integers
527	REAL(wp),DIMENSION(jpi,jpj,jpk) :: zwd, zwi, zws, zwrm, zwt
528	!!----------------------------------------------------------------------
529	!
530	! !== build the three diagonal matrix & the RHS ==!
531	!
532	DO_3D_00_00( 3, jpkm1 )
533	zwd (ji,jj,jk) = 3._wp * wmask(ji,jj,jk) + 1._wp ! diagonal
534	zwi (ji,jj,jk) = wmask(ji,jj,jk) ! lower diagonal
535	zws (ji,jj,jk) = wmask(ji,jj,jk) ! upper diagonal
536	zwrm(ji,jj,jk) = 3._wp * wmask(ji,jj,jk) & ! RHS
537	& * ( pt_in(ji,jj,jk) + pt_in(ji,jj,jk-1) )
538	END_3D
539	!
540	!!gm
541	! SELECT CASE( kbc ) !* boundary condition
542	! CASE( np_NH ) ! Neumann homogeneous at top & bottom
543	! CASE( np_CEN2 ) ! 2nd order centered at top & bottom
544	! END SELECT
545	!!gm
546	!
547	IF ( ln_isfcav ) THEN ! set level two values which may not be set in ISF case
548	zwd(:,:,2) = 1._wp ; zwi(:,:,2) = 0._wp ; zws(:,:,2) = 0._wp ; zwrm(:,:,2) = 0._wp
549	END IF
550	!
551	DO_2D_00_00
552	ikt = mikt(ji,jj) + 1 ! w-point below the 1st wet point
553	ikb = MAX(mbkt(ji,jj), 2) ! - above the last wet point
554	!
555	zwd (ji,jj,ikt) = 1._wp ! top
556	zwi (ji,jj,ikt) = 0._wp
557	zws (ji,jj,ikt) = 0._wp
558	zwrm(ji,jj,ikt) = 0.5_wp * ( pt_in(ji,jj,ikt-1) + pt_in(ji,jj,ikt) )
559	!
560	zwd (ji,jj,ikb) = 1._wp ! bottom
561	zwi (ji,jj,ikb) = 0._wp
562	zws (ji,jj,ikb) = 0._wp
563	zwrm(ji,jj,ikb) = 0.5_wp * ( pt_in(ji,jj,ikb-1) + pt_in(ji,jj,ikb) )
564	END_2D
565	!
566	! !== tridiagonal solver ==!
567	!
568	DO_2D_00_00
569	zwt(ji,jj,2) = zwd(ji,jj,2)
570	END_2D
571	DO_3D_00_00( 3, jpkm1 )
572	zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) /zwt(ji,jj,jk-1)
573	END_3D
574	!
575	DO_2D_00_00
576	pt_out(ji,jj,2) = zwrm(ji,jj,2)
577	END_2D
578	DO_3D_00_00( 3, jpkm1 )
579	pt_out(ji,jj,jk) = zwrm(ji,jj,jk) - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1)
580	END_3D
581
582	DO_2D_00_00
583	pt_out(ji,jj,jpkm1) = pt_out(ji,jj,jpkm1) / zwt(ji,jj,jpkm1)
584	END_2D
585	DO_3DS_00_00( jpk-2, 2, -1 )
586	pt_out(ji,jj,jk) = ( pt_out(ji,jj,jk) - zws(ji,jj,jk) * pt_out(ji,jj,jk+1) ) / zwt(ji,jj,jk)
587	END_3D
588	!
589	END SUBROUTINE interp_4th_cpt
590
591
592	SUBROUTINE tridia_solver( pD, pU, pL, pRHS, pt_out , klev )
593	!!----------------------------------------------------------------------
594	!! * ROUTINE tridia_solver *
595	!!
596	!! ** Purpose : solve a symmetric 3diagonal system
597	!!
598	!! ** Method : solve M.t_out = RHS(t) where M is a tri diagonal matrix ( jpk*jpk )
599	!!
600	!! ( D_1 U_1 0 0 0 )( t_1 ) ( RHS_1 )
601	!! ( L_2 D_2 U_2 0 0 )( t_2 ) ( RHS_2 )
602	!! ( 0 L_3 D_3 U_3 0 )( t_3 ) = ( RHS_3 )
603	!! ( ... )( ... ) ( ... )
604	!! ( 0 0 0 L_k D_k )( t_k ) ( RHS_k )
605	!!
606	!! M is decomposed in the product of an upper and lower triangular matrix.
607	!! The tri-diagonals matrix is given as input 3D arrays: pD, pU, pL
608	!! (i.e. the Diagonal, the Upper diagonal, and the Lower diagonal).
609	!! The solution is pta.
610	!! The 3d array zwt is used as a work space array.
611	!!----------------------------------------------------------------------
612	REAL(wp),DIMENSION(:,:,:), INTENT(in ) :: pD, pU, PL ! 3-diagonal matrix
613	REAL(wp),DIMENSION(:,:,:), INTENT(in ) :: pRHS ! Right-Hand-Side
614	REAL(wp),DIMENSION(:,:,:), INTENT( out) :: pt_out !!gm field at level=F(klev)
615	INTEGER , INTENT(in ) :: klev ! =1 pt_out at w-level
616	! ! =0 pt at t-level
617	INTEGER :: ji, jj, jk ! dummy loop integers
618	INTEGER :: kstart ! local indices
619	REAL(wp),DIMENSION(jpi,jpj,jpk) :: zwt ! 3D work array
620	!!----------------------------------------------------------------------
621	!
622	kstart = 1 + klev
623	!
624	DO_2D_00_00
625	zwt(ji,jj,kstart) = pD(ji,jj,kstart)
626	END_2D
627	DO_3D_00_00( kstart+1, jpkm1 )
628	zwt(ji,jj,jk) = pD(ji,jj,jk) - pL(ji,jj,jk) * pU(ji,jj,jk-1) /zwt(ji,jj,jk-1)
629	END_3D
630	!
631	DO_2D_00_00
632	pt_out(ji,jj,kstart) = pRHS(ji,jj,kstart)
633	END_2D
634	DO_3D_00_00( kstart+1, jpkm1 )
635	pt_out(ji,jj,jk) = pRHS(ji,jj,jk) - pL(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1)
636	END_3D
637
638	DO_2D_00_00
639	pt_out(ji,jj,jpkm1) = pt_out(ji,jj,jpkm1) / zwt(ji,jj,jpkm1)
640	END_2D
641	DO_3DS_00_00( jpk-2, kstart, -1 )
642	pt_out(ji,jj,jk) = ( pt_out(ji,jj,jk) - pU(ji,jj,jk) * pt_out(ji,jj,jk+1) ) / zwt(ji,jj,jk)
643	END_3D
644	!
645	END SUBROUTINE tridia_solver
646
647	!!======================================================================
648	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/trunk/src/OCE/TRA/traadv_fct.F90 @ 12396

Download in other formats: