Context Navigation

traadv_fct.F90 @ 10136

Last change on this file since 10136 was 10136, checked in by dguibert, 6 years ago

bull: async/datatype

Experimental changes to enable/study/bench various mpi "optimisations":

BULL_ASYNC
BULL_DATATYPE_VECTOR/SUBARRAY

this has been applied to the nonosc subroutine (only for now).

Property svn:keywords set to Id

File size: 43.4 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	!
24	USE in_out_manager ! I/O manager
25	USE iom !
26	USE lib_mpp ! MPP library
27	USE lbclnk ! ocean lateral boundary condition (or mpp link)
28	USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
29	USE timing ! Timing
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 licence (./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	!!----------------------------------------------------------------------
90	!
91	IF( kt == kit000 ) THEN
92	IF(lwp) WRITE(numout,*)
93	IF(lwp) WRITE(numout,*) 'tra_adv_fct : FCT advection scheme on ', cdtype
94	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
95	ENDIF
96	!
97	l_trd = .FALSE. ! set local switches
98	l_hst = .FALSE.
99	l_ptr = .FALSE.
100	IF( ( cdtype =='TRA' .AND. l_trdtra ) .OR. ( cdtype =='TRC' .AND. l_trdtrc ) ) l_trd = .TRUE.
101	IF( cdtype =='TRA' .AND. ln_diaptr ) l_ptr = .TRUE.
102	IF( cdtype =='TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. &
103	& iom_use("uadv_salttr") .OR. iom_use("vadv_salttr") ) ) l_hst = .TRUE.
104	!
105	IF( l_trd .OR. l_hst ) THEN
106	ALLOCATE( ztrdx(jpi,jpj,jpk), ztrdy(jpi,jpj,jpk), ztrdz(jpi,jpj,jpk) )
107	ztrdx(:,:,:) = 0._wp ; ztrdy(:,:,:) = 0._wp ; ztrdz(:,:,:) = 0._wp
108	ENDIF
109	!
110	IF( l_ptr ) THEN
111	ALLOCATE( zptry(jpi,jpj,jpk) )
112	zptry(:,:,:) = 0._wp
113	ENDIF
114	! ! surface & bottom value : flux set to zero one for all
115	zwz(:,:, 1 ) = 0._wp
116	zwx(:,:,jpk) = 0._wp ; zwy(:,:,jpk) = 0._wp ; zwz(:,:,jpk) = 0._wp
117	!
118	zwi(:,:,:) = 0._wp
119	!
120	DO jn = 1, kjpt !== loop over the tracers ==!
121	!
122	! !== upstream advection with initial mass fluxes & intermediate update ==!
123	! !* upstream tracer flux in the i and j direction
124	DO jk = 1, jpkm1
125	DO jj = 1, jpjm1
126	DO ji = 1, fs_jpim1 ! vector opt.
127	! upstream scheme
128	zfp_ui = pun(ji,jj,jk) + ABS( pun(ji,jj,jk) )
129	zfm_ui = pun(ji,jj,jk) - ABS( pun(ji,jj,jk) )
130	zfp_vj = pvn(ji,jj,jk) + ABS( pvn(ji,jj,jk) )
131	zfm_vj = pvn(ji,jj,jk) - ABS( pvn(ji,jj,jk) )
132	zwx(ji,jj,jk) = 0.5 * ( zfp_ui * ptb(ji,jj,jk,jn) + zfm_ui * ptb(ji+1,jj ,jk,jn) )
133	zwy(ji,jj,jk) = 0.5 * ( zfp_vj * ptb(ji,jj,jk,jn) + zfm_vj * ptb(ji ,jj+1,jk,jn) )
134	END DO
135	END DO
136	END DO
137	! !* upstream tracer flux in the k direction *!
138	DO jk = 2, jpkm1 ! Interior value ( multiplied by wmask)
139	DO jj = 1, jpj
140	DO ji = 1, jpi
141	zfp_wk = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) )
142	zfm_wk = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) )
143	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)
144	END DO
145	END DO
146	END DO
147	IF( ln_linssh ) THEN ! top ocean value (only in linear free surface as zwz has been w-masked)
148	IF( ln_isfcav ) THEN ! top of the ice-shelf cavities and at the ocean surface
149	DO jj = 1, jpj
150	DO ji = 1, jpi
151	zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptb(ji,jj,mikt(ji,jj),jn) ! linear free surface
152	END DO
153	END DO
154	ELSE ! no cavities: only at the ocean surface
155	zwz(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn)
156	ENDIF
157	ENDIF
158	!
159	DO jk = 1, jpkm1 !* trend and after field with monotonic scheme
160	DO jj = 2, jpjm1
161	DO ji = fs_2, fs_jpim1 ! vector opt.
162	! ! total intermediate advective trends
163	ztra = - ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) &
164	& + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) &
165	& + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) * r1_e1e2t(ji,jj)
166	! ! update and guess with monotonic sheme
167	pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra / e3t_n(ji,jj,jk) * tmask(ji,jj,jk)
168	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)
169	END DO
170	END DO
171	END DO
172	CALL lbc_lnk("traadv_fct",zwi, 'T', 1. ) ! Lateral boundary conditions on zwi (unchanged sign)
173	!
174	IF( l_trd .OR. l_hst ) THEN ! trend diagnostics (contribution of upstream fluxes)
175	ztrdx(:,:,:) = zwx(:,:,:) ; ztrdy(:,:,:) = zwy(:,:,:) ; ztrdz(:,:,:) = zwz(:,:,:)
176	END IF
177	! ! "Poleward" heat and salt transports (contribution of upstream fluxes)
178	IF( l_ptr ) zptry(:,:,:) = zwy(:,:,:)
179	!
180	! !== anti-diffusive flux : high order minus low order ==!
181	!
182	SELECT CASE( kn_fct_h ) !* horizontal anti-diffusive fluxes
183	!
184	CASE( 2 ) !- 2nd order centered
185	DO jk = 1, jpkm1
186	DO jj = 1, jpjm1
187	DO ji = 1, fs_jpim1 ! vector opt.
188	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)
189	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)
190	END DO
191	END DO
192	END DO
193	!
194	CASE( 4 ) !- 4th order centered
195	zltu(:,:,jpk) = 0._wp ! Bottom value : flux set to zero
196	zltv(:,:,jpk) = 0._wp
197	DO jk = 1, jpkm1 ! Laplacian
198	DO jj = 1, jpjm1 ! 1st derivative (gradient)
199	DO ji = 1, fs_jpim1 ! vector opt.
200	ztu(ji,jj,jk) = ( ptn(ji+1,jj ,jk,jn) - ptn(ji,jj,jk,jn) ) * umask(ji,jj,jk)
201	ztv(ji,jj,jk) = ( ptn(ji ,jj+1,jk,jn) - ptn(ji,jj,jk,jn) ) * vmask(ji,jj,jk)
202	END DO
203	END DO
204	DO jj = 2, jpjm1 ! 2nd derivative * 1/ 6
205	DO ji = fs_2, fs_jpim1 ! vector opt.
206	zltu(ji,jj,jk) = ( ztu(ji,jj,jk) + ztu(ji-1,jj,jk) ) * r1_6
207	zltv(ji,jj,jk) = ( ztv(ji,jj,jk) + ztv(ji,jj-1,jk) ) * r1_6
208	END DO
209	END DO
210	END DO
211	CALL lbc_lnk_multi("traadv_fct",zltu, 'T', 1. , zltv, 'T', 1. ) ! Lateral boundary cond. (unchanged sgn)
212	!
213	DO jk = 1, jpkm1 ! Horizontal advective fluxes
214	DO jj = 1, jpjm1
215	DO ji = 1, fs_jpim1 ! vector opt.
216	zC2t_u = ptn(ji,jj,jk,jn) + ptn(ji+1,jj ,jk,jn) ! 2 x C2 interpolation of T at u- & v-points
217	zC2t_v = ptn(ji,jj,jk,jn) + ptn(ji ,jj+1,jk,jn)
218	! ! C4 minus upstream advective fluxes
219	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)
220	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)
221	END DO
222	END DO
223	END DO
224	!
225	CASE( 41 ) !- 4th order centered ==>> !!gm coding attempt need to be tested
226	ztu(:,:,jpk) = 0._wp ! Bottom value : flux set to zero
227	ztv(:,:,jpk) = 0._wp
228	DO jk = 1, jpkm1 ! 1st derivative (gradient)
229	DO jj = 1, jpjm1
230	DO ji = 1, fs_jpim1 ! vector opt.
231	ztu(ji,jj,jk) = ( ptn(ji+1,jj ,jk,jn) - ptn(ji,jj,jk,jn) ) * umask(ji,jj,jk)
232	ztv(ji,jj,jk) = ( ptn(ji ,jj+1,jk,jn) - ptn(ji,jj,jk,jn) ) * vmask(ji,jj,jk)
233	END DO
234	END DO
235	END DO
236	CALL lbc_lnk_multi("traadv_fct",ztu, 'U', -1. , ztv, 'V', -1. ) ! Lateral boundary cond. (unchanged sgn)
237	!
238	DO jk = 1, jpkm1 ! Horizontal advective fluxes
239	DO jj = 2, jpjm1
240	DO ji = 2, fs_jpim1 ! vector opt.
241	zC2t_u = ptn(ji,jj,jk,jn) + ptn(ji+1,jj ,jk,jn) ! 2 x C2 interpolation of T at u- & v-points (x2)
242	zC2t_v = ptn(ji,jj,jk,jn) + ptn(ji ,jj+1,jk,jn)
243	! ! C4 interpolation of T at u- & v-points (x2)
244	zC4t_u = zC2t_u + r1_6 * ( ztu(ji-1,jj ,jk) - ztu(ji+1,jj ,jk) )
245	zC4t_v = zC2t_v + r1_6 * ( ztv(ji ,jj-1,jk) - ztv(ji ,jj+1,jk) )
246	! ! C4 minus upstream advective fluxes
247	zwx(ji,jj,jk) = 0.5_wp * pun(ji,jj,jk) * zC4t_u - zwx(ji,jj,jk)
248	zwy(ji,jj,jk) = 0.5_wp * pvn(ji,jj,jk) * zC4t_v - zwy(ji,jj,jk)
249	END DO
250	END DO
251	END DO
252	!
253	END SELECT
254	!
255	SELECT CASE( kn_fct_v ) !* vertical anti-diffusive fluxes (w-masked interior values)
256	!
257	CASE( 2 ) !- 2nd order centered
258	DO jk = 2, jpkm1
259	DO jj = 2, jpjm1
260	DO ji = fs_2, fs_jpim1
261	zwz(ji,jj,jk) = ( pwn(ji,jj,jk) * 0.5_wp * ( ptn(ji,jj,jk,jn) + ptn(ji,jj,jk-1,jn) ) &
262	& - zwz(ji,jj,jk) ) * wmask(ji,jj,jk)
263	END DO
264	END DO
265	END DO
266	!
267	CASE( 4 ) !- 4th order COMPACT
268	CALL interp_4th_cpt( ptn(:,:,:,jn) , ztw ) ! zwt = COMPACT interpolation of T at w-point
269	DO jk = 2, jpkm1
270	DO jj = 2, jpjm1
271	DO ji = fs_2, fs_jpim1
272	zwz(ji,jj,jk) = ( pwn(ji,jj,jk) * ztw(ji,jj,jk) - zwz(ji,jj,jk) ) * wmask(ji,jj,jk)
273	END DO
274	END DO
275	END DO
276	!
277	END SELECT
278	IF( ln_linssh ) THEN ! top ocean value: high order = upstream ==>> zwz=0
279	zwz(:,:,1) = 0._wp ! only ocean surface as interior zwz values have been w-masked
280	ENDIF
281	!
282	CALL lbc_lnk_multi("traadv_fct",zwx, 'U', -1. , zwy, 'V', -1., zwz, 'W', 1. )
283	!
284	! !== monotonicity algorithm ==!
285	!
286	CALL nonosc( ptb(:,:,:,jn), zwx, zwy, zwz, zwi, p2dt )
287	!
288	! !== final trend with corrected fluxes ==!
289	!
290	DO jk = 1, jpkm1
291	DO jj = 2, jpjm1
292	DO ji = fs_2, fs_jpim1 ! vector opt.
293	pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) - ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) &
294	& + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) &
295	& + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) &
296	& * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
297	END DO
298	END DO
299	END DO
300	!
301	IF( l_trd .OR. l_hst ) THEN ! trend diagnostics // heat/salt transport
302	ztrdx(:,:,:) = ztrdx(:,:,:) + zwx(:,:,:) ! <<< add anti-diffusive fluxes
303	ztrdy(:,:,:) = ztrdy(:,:,:) + zwy(:,:,:) ! to upstream fluxes
304	ztrdz(:,:,:) = ztrdz(:,:,:) + zwz(:,:,:) !
305	!
306	IF( l_trd ) THEN ! trend diagnostics
307	CALL trd_tra( kt, cdtype, jn, jptra_xad, ztrdx, pun, ptn(:,:,:,jn) )
308	CALL trd_tra( kt, cdtype, jn, jptra_yad, ztrdy, pvn, ptn(:,:,:,jn) )
309	CALL trd_tra( kt, cdtype, jn, jptra_zad, ztrdz, pwn, ptn(:,:,:,jn) )
310	ENDIF
311	! ! heat/salt transport
312	IF( l_hst ) CALL dia_ar5_hst( jn, 'adv', ztrdx(:,:,:), ztrdy(:,:,:) )
313	!
314	DEALLOCATE( ztrdx, ztrdy, ztrdz )
315	ENDIF
316	IF( l_ptr ) THEN ! "Poleward" transports
317	zptry(:,:,:) = zptry(:,:,:) + zwy(:,:,:) ! <<< add anti-diffusive fluxes
318	CALL dia_ptr_hst( jn, 'adv', zptry(:,:,:) )
319	DEALLOCATE( zptry )
320	ENDIF
321	!
322	END DO ! end of tracer loop
323	!
324	END SUBROUTINE tra_adv_fct
325
326
327	SUBROUTINE nonosc( pbef, paa, pbb, pcc, paft, p2dt )
328	#ifdef SCOREP_USER_ENABLE
329	use mpi
330	#include "scorep/SCOREP_User.inc"
331	#endif
332	!!---------------------------------------------------------------------
333	!! * ROUTINE nonosc *
334	!!
335	!! ** Purpose : compute monotonic tracer fluxes from the upstream
336	!! scheme and the before field by a nonoscillatory algorithm
337	!!
338	!! ** Method : ... ???
339	!! warning : pbef and paft must be masked, but the boundaries
340	!! conditions on the fluxes are not necessary zalezak (1979)
341	!! drange (1995) multi-dimensional forward-in-time and upstream-
342	!! in-space based differencing for fluid
343	!!----------------------------------------------------------------------
344	REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step
345	REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(in ) :: pbef, paft ! before & after field
346	REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(inout) :: paa, pbb, pcc ! monotonic fluxes in the 3 directions
347	!
348	INTEGER :: ji, jj, jk ! dummy loop indices
349	INTEGER :: ikm1 ! local integer
350	REAL(wp) :: zpos, zneg, zbt, za, zb, zc, zbig, zrtrn ! local scalars
351	REAL(wp) :: zau, zbu, zcu, zav, zbv, zcv, zup, zdo ! - -
352	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zbetup, zbetdo, zbup, zbdo
353	!dir$ attributes align:64 :: zbetup, zbetdo, zbup, zbdo
354	#ifdef SCOREP_USER_ENABLE
355	integer :: ier
356	SCOREP_USER_REGION_DEFINE( reg_nonosc )
357	SCOREP_USER_REGION_DEFINE( reg_nonosc_setup )
358	SCOREP_USER_REGION_DEFINE( reg_nonosc_cb1 )
359	SCOREP_USER_REGION_DEFINE( reg_nonosc_cb2 )
360	SCOREP_USER_REGION_DEFINE( reg_nonosc_barrier )
361	SCOREP_USER_REGION_DEFINE( reg_nonosc_imbalance )
362
363	SCOREP_USER_REGION_BEGIN( reg_nonosc_barrier, "nonosc barrier", SCOREP_USER_REGION_TYPE_COMMON )
364	call MPI_Barrier(MPI_COMM_WORLD, ier)
365	SCOREP_USER_REGION_END( reg_nonosc_barrier )
366	SCOREP_USER_REGION_BEGIN( reg_nonosc, "nonosc", SCOREP_USER_REGION_TYPE_FUNCTION )
367	SCOREP_USER_REGION_BEGIN( reg_nonosc_setup, "nonosc setup", SCOREP_USER_REGION_TYPE_COMMON )
368	#endif
369	IF( ln_timing ) CALL timing_start( 'nonosc' )
370	!!----------------------------------------------------------------------
371	!
372	zbig = 1.e+40_wp
373	zrtrn = 1.e-15_wp
374	#ifndef BULL_NONOSC_INIT
375	zbetup(:,:,:) = 0._wp ; zbetdo(:,:,:) = 0._wp
376	#else
377	zbetup(:,:,jpk) = 0._wp ; zbetdo(:,:,jpk) = 0._wp
378	#endif
379
380	! Search local extrema
381	! --------------------
382	! max/min of pbef & paft with large negative/positive value (-/+zbig) inside land
383	zbup = MAX( pbef * tmask - zbig * ( 1._wp - tmask ), &
384	& paft * tmask - zbig * ( 1._wp - tmask ) )
385	zbdo = MIN( pbef * tmask + zbig * ( 1._wp - tmask ), &
386	& paft * tmask + zbig * ( 1._wp - tmask ) )
387
388	#ifdef SCOREP_USER_ENABLE
389	SCOREP_USER_REGION_END( reg_nonosc_setup )
390	#endif
391
392	#ifndef BULL_ASYNC
393	#ifdef SCOREP_USER_ENABLE
394	SCOREP_USER_REGION_BEGIN( reg_nonosc_cb1, "cb1", SCOREP_USER_REGION_TYPE_LOOP )
395	#endif
396	! loads:
397	! - zbup: ji-1/ji/ji+1, jj-1/jj/jj+1, ji/jk+1/jk-1
398	! - zbdo: "
399	! - paa: ji-1/ji
400	! - pbb: jj-1/jj
401	! - pcc: ji, jj, jk/jk+1
402	! - e1e2t, e3t_n, paft (*2): ji,jj,jk
403	!
404	! stores:
405	! - zbetup
406	! - zbetdo
407	DO jk = 1, jpkm1
408	ikm1 = MAX(jk-1,1)
409	DO jj = 2, jpjm1
410	DO ji = fs_2, fs_jpim1 ! vector opt.
411
412	! search maximum in neighbourhood
413	zup = MAX( zbup(ji ,jj ,jk ), &
414	& zbup(ji-1,jj ,jk ), zbup(ji+1,jj ,jk ), &
415	& zbup(ji ,jj-1,jk ), zbup(ji ,jj+1,jk ), &
416	& zbup(ji ,jj ,ikm1), zbup(ji ,jj ,jk+1) )
417
418	! search minimum in neighbourhood
419	zdo = MIN( zbdo(ji ,jj ,jk ), &
420	& zbdo(ji-1,jj ,jk ), zbdo(ji+1,jj ,jk ), &
421	& zbdo(ji ,jj-1,jk ), zbdo(ji ,jj+1,jk ), &
422	& zbdo(ji ,jj ,ikm1), zbdo(ji ,jj ,jk+1) )
423
424	! positive part of the flux
425	zpos = MAX( 0., paa(ji-1,jj ,jk ) ) - MIN( 0., paa(ji ,jj ,jk ) ) &
426	& + MAX( 0., pbb(ji ,jj-1,jk ) ) - MIN( 0., pbb(ji ,jj ,jk ) ) &
427	& + MAX( 0., pcc(ji ,jj ,jk+1) ) - MIN( 0., pcc(ji ,jj ,jk ) )
428
429	! negative part of the flux
430	zneg = MAX( 0., paa(ji ,jj ,jk ) ) - MIN( 0., paa(ji-1,jj ,jk ) ) &
431	& + MAX( 0., pbb(ji ,jj ,jk ) ) - MIN( 0., pbb(ji ,jj-1,jk ) ) &
432	& + MAX( 0., pcc(ji ,jj ,jk ) ) - MIN( 0., pcc(ji ,jj ,jk+1) )
433
434	! up & down beta terms
435	zbt = e1e2t(ji,jj) * e3t_n(ji,jj,jk) / p2dt
436	zbetup(ji,jj,jk) = ( zup - paft(ji,jj,jk) ) / ( zpos + zrtrn ) * zbt
437	zbetdo(ji,jj,jk) = ( paft(ji,jj,jk) - zdo ) / ( zneg + zrtrn ) * zbt
438	END DO
439	END DO
440	END DO
441	#ifdef SCOREP_USER_ENABLE
442	SCOREP_USER_REGION_END( reg_nonosc_cb1 )
443	#endif
444	CALL lbc_lnk_multi("traadv_fct",zbetup, 'T', 1. , zbetdo, 'T', 1. ) ! lateral boundary cond. (unchanged sign)
445	#else
446	call lbc_lnk_multi_async( "traadv_fct", cb1, zbetup, 'T', 1. , zbetdo, 'T', 1. ) ! lateral boundary cond. (unchanged sign)
447	#endif
448
449	#ifndef BULL_ASYNC
450	! 3. monotonic flux in the i & j direction (paa & pbb)
451	! ----------------------------------------
452	#ifdef SCOREP_USER_ENABLE
453	SCOREP_USER_REGION_BEGIN( reg_nonosc_cb2, "cb2", SCOREP_USER_REGION_TYPE_LOOP )
454	#endif
455	DO jk = 1, jpkm1
456	DO jj = 2, jpjm1
457	DO ji = fs_2, fs_jpim1 ! vector opt.
458	zau = MIN( 1._wp, zbetdo(ji,jj,jk), zbetup(ji+1,jj,jk) )
459	zbu = MIN( 1._wp, zbetup(ji,jj,jk), zbetdo(ji+1,jj,jk) )
460	zcu = ( 0.5 + SIGN( 0.5 , paa(ji,jj,jk) ) )
461	paa(ji,jj,jk) = paa(ji,jj,jk) * ( zcu * zau + ( 1._wp - zcu) * zbu )
462
463	zav = MIN( 1._wp, zbetdo(ji,jj,jk), zbetup(ji,jj+1,jk) )
464	zbv = MIN( 1._wp, zbetup(ji,jj,jk), zbetdo(ji,jj+1,jk) )
465	zcv = ( 0.5 + SIGN( 0.5 , pbb(ji,jj,jk) ) )
466	pbb(ji,jj,jk) = pbb(ji,jj,jk) * ( zcv * zav + ( 1._wp - zcv) * zbv )
467
468	! monotonic flux in the k direction, i.e. pcc
469	! -------------------------------------------
470	za = MIN( 1., zbetdo(ji,jj,jk+1), zbetup(ji,jj,jk) )
471	zb = MIN( 1., zbetup(ji,jj,jk+1), zbetdo(ji,jj,jk) )
472	zc = ( 0.5 + SIGN( 0.5 , pcc(ji,jj,jk+1) ) )
473	pcc(ji,jj,jk+1) = pcc(ji,jj,jk+1) * ( zc * za + ( 1._wp - zc) * zb )
474	END DO
475	END DO
476	END DO
477	#ifdef SCOREP_USER_ENABLE
478	SCOREP_USER_REGION_END( reg_nonosc_cb2 )
479	#endif
480	CALL lbc_lnk_multi("traadv_fct",paa, 'U', -1. , pbb, 'V', -1. ) ! lateral boundary condition (changed sign)
481	#else
482	call lbc_lnk_multi_async( "traadv_fct", cb2, paa, 'U', -1. , pbb, 'V', -1. ) ! lateral boundary condition (changed sign)
483	#endif
484	!
485	IF( ln_timing ) CALL timing_stop( 'nonosc' )
486	#ifdef SCOREP_USER_ENABLE
487	SCOREP_USER_REGION_END( reg_nonosc )
488	SCOREP_USER_REGION_BEGIN( reg_nonosc_imbalance, "nonosc imbalance", SCOREP_USER_REGION_TYPE_COMMON )
489	call MPI_Barrier(MPI_COMM_WORLD, ier)
490	SCOREP_USER_REGION_END( reg_nonosc_imbalance )
491	#endif
492	#ifdef BULL_ASYNC
493	contains
494	subroutine cb1(i0, i1, j0, j1, k0, k1, buf)
495	integer, intent(in) :: i0, i1, j0, j1, k0, k1
496	real(wp), dimension(:,:,:,:,:,:), optional, intent(out) :: buf
497	integer jji, jjj, jjk
498	real(wp) :: p2dt_inv
499	!REAL(wp), DIMENSION (40,jpj,jpk) :: paa, pbb, pcc ! monotonic fluxes in the 3 directions
500	!REAL(wp), DIMENSION (40,jpj,jpk) :: e3t_n, paft
501	!REAL(wp), DIMENSION (40,jpj) :: e1e2t
502	!REAL(wp), DIMENSION(40,jpj,jpk) :: zbetup, zbetdo, zbup, zbdo
503	!DIR$ ASSUME_ALIGNED zbup:64
504	!DIR$ ASSUME (jpi .EQ.40)
505	!DIR$ ASSUME (jpj .EQ.42)
506	!DIR$ ASSUME (jpk .EQ.75)
507
508	p2dt_inv = 1._wp * p2dt
509	if(i0 == i1) then
510	ji=i0
511	! DO jjj = j0, j1, 8
512	DO jk = k0, k1
513	ikm1 = MAX(jk-1,1)
514	!DIR$ vector always
515	DO jj = j0, j1
516	!DO jj = jjj, min(jjj+7, j1)
517	! search maximum in neighbourhood
518	zup = MAX( zbup(ji ,jj ,jk ), &
519	& zbup(ji-1,jj ,jk ), zbup(ji+1,jj ,jk ), &
520	& zbup(ji ,jj-1,jk ), zbup(ji ,jj+1,jk ), &
521	& zbup(ji ,jj ,ikm1), zbup(ji ,jj ,jk+1) )
522
523	! search minimum in neighbourhood
524	zdo = MIN( zbdo(ji ,jj ,jk ), &
525	& zbdo(ji-1,jj ,jk ), zbdo(ji+1,jj ,jk ), &
526	& zbdo(ji ,jj-1,jk ), zbdo(ji ,jj+1,jk ), &
527	& zbdo(ji ,jj ,ikm1), zbdo(ji ,jj ,jk+1) )
528
529	! positive part of the flux
530	zpos = MAX( 0., paa(ji-1,jj ,jk ) ) - MIN( 0., paa(ji ,jj ,jk ) ) &
531	& + MAX( 0., pbb(ji ,jj-1,jk ) ) - MIN( 0., pbb(ji ,jj ,jk ) ) &
532	& + MAX( 0., pcc(ji ,jj ,jk+1) ) - MIN( 0., pcc(ji ,jj ,jk ) )
533
534	! negative part of the flux
535	zneg = MAX( 0., paa(ji ,jj ,jk ) ) - MIN( 0., paa(ji-1,jj ,jk ) ) &
536	& + MAX( 0., pbb(ji ,jj ,jk ) ) - MIN( 0., pbb(ji ,jj-1,jk ) ) &
537	& + MAX( 0., pcc(ji ,jj ,jk ) ) - MIN( 0., pcc(ji ,jj ,jk+1) )
538
539	! up & down beta terms
540	zbt = e1e2t(ji,jj) * e3t_n(ji,jj,jk) * p2dt_inv
541	zbetup(ji,jj,jk) = ( zup - paft(ji,jj,jk) ) / ( zpos + zrtrn ) * zbt
542	zbetdo(ji,jj,jk) = ( paft(ji,jj,jk) - zdo ) / ( zneg + zrtrn ) * zbt
543
544	#ifdef BULL_CB_WITH_BUF
545	! zt3ew(:,jh,jk,jl,jf,1) = ARRAY_IN(nn_hls+jh,:,jk,jl,jf
546	buf(jj,1,jk,1,1,1) = zbetup(ji,jj,jk)
547	buf(jj,1,jk,1,2,1) = zbetdo(ji,jj,jk)
548	#endif
549	END DO
550	END DO
551	!end do
552	else
553	DO jk = k0, k1
554	ikm1 = MAX(jk-1,1)
555	DO jj = j0, j1
556	!DIR$ vector always
557	DO ji = i0, i1
558
559	! search maximum in neighbourhood
560	zup = MAX( zbup(ji ,jj ,jk ), &
561	& zbup(ji-1,jj ,jk ), zbup(ji+1,jj ,jk ), &
562	& zbup(ji ,jj-1,jk ), zbup(ji ,jj+1,jk ), &
563	& zbup(ji ,jj ,ikm1), zbup(ji ,jj ,jk+1) )
564
565	! search minimum in neighbourhood
566	zdo = MIN( zbdo(ji ,jj ,jk ), &
567	& zbdo(ji-1,jj ,jk ), zbdo(ji+1,jj ,jk ), &
568	& zbdo(ji ,jj-1,jk ), zbdo(ji ,jj+1,jk ), &
569	& zbdo(ji ,jj ,ikm1), zbdo(ji ,jj ,jk+1) )
570
571	! positive part of the flux
572	zpos = MAX( 0., paa(ji-1,jj ,jk ) ) - MIN( 0., paa(ji ,jj ,jk ) ) &
573	& + MAX( 0., pbb(ji ,jj-1,jk ) ) - MIN( 0., pbb(ji ,jj ,jk ) ) &
574	& + MAX( 0., pcc(ji ,jj ,jk+1) ) - MIN( 0., pcc(ji ,jj ,jk ) )
575
576	! negative part of the flux
577	zneg = MAX( 0., paa(ji ,jj ,jk ) ) - MIN( 0., paa(ji-1,jj ,jk ) ) &
578	& + MAX( 0., pbb(ji ,jj ,jk ) ) - MIN( 0., pbb(ji ,jj-1,jk ) ) &
579	& + MAX( 0., pcc(ji ,jj ,jk ) ) - MIN( 0., pcc(ji ,jj ,jk+1) )
580
581	! up & down beta terms
582	zbt = e1e2t(ji,jj) * e3t_n(ji,jj,jk) * p2dt_inv
583	zbetup(ji,jj,jk) = ( zup - paft(ji,jj,jk) ) / ( zpos + zrtrn ) * zbt
584	zbetdo(ji,jj,jk) = ( paft(ji,jj,jk) - zdo ) / ( zneg + zrtrn ) * zbt
585
586	END DO
587	END DO
588	END DO
589	endif
590
591	end subroutine
592	subroutine cb2(i0, i1, j0, j1, k0, k1, buf)
593	integer, intent(in) :: i0, i1, j0, j1, k0, k1
594	real(wp), dimension(:,:,:,:,:,:), optional, intent(out) :: buf
595	integer jji, jjj, jjk
596
597	! 3. monotonic flux in the i & j direction (paa & pbb)
598	if(i0 == i1) then
599	ji=i0
600	do jjj=j0, j1, 8
601	DO jk = k0, k1
602	!DIR$ vector always
603	!DO jj = j0, j1
604	DO jj = jjj, min(jjj+7, j1)
605	zau = MIN( 1._wp, zbetdo(ji,jj,jk), zbetup(ji+1,jj,jk) )
606	zbu = MIN( 1._wp, zbetup(ji,jj,jk), zbetdo(ji+1,jj,jk) )
607	zcu = ( 0.5 + SIGN( 0.5 , paa(ji,jj,jk) ) )
608	paa(ji,jj,jk) = paa(ji,jj,jk) * ( zcu * zau + ( 1._wp - zcu) * zbu )
609
610	zav = MIN( 1._wp, zbetdo(ji,jj,jk), zbetup(ji,jj+1,jk) )
611	zbv = MIN( 1._wp, zbetup(ji,jj,jk), zbetdo(ji,jj+1,jk) )
612	zcv = ( 0.5 + SIGN( 0.5 , pbb(ji,jj,jk) ) )
613	pbb(ji,jj,jk) = pbb(ji,jj,jk) * ( zcv * zav + ( 1._wp - zcv) * zbv )
614
615	! monotonic flux in the k direction, i.e. pcc
616	! -------------------------------------------
617	za = MIN( 1., zbetdo(ji,jj,jk+1), zbetup(ji,jj,jk) )
618	zb = MIN( 1., zbetup(ji,jj,jk+1), zbetdo(ji,jj,jk) )
619	zc = ( 0.5 + SIGN( 0.5 , pcc(ji,jj,jk+1) ) )
620	pcc(ji,jj,jk+1) = pcc(ji,jj,jk+1) * ( zc * za + ( 1._wp - zc) * zb )
621	#ifdef BULL_CB_WITH_BUF
622	! zt3ew(:,jh,jk,jl,jf,1) = ARRAY_IN(nn_hls+jh,:,jk,jl,jf
623	buf(jj,1,jk,1,1,1) = paa(ji,jj,jk)
624	buf(jj,1,jk,1,2,1) = pbb(ji,jj,jk)
625	#endif
626	END DO
627	END DO
628	end do
629	else
630	DO jk = k0, k1
631	DO jj = j0, j1
632	!DIR$ vector always
633	DO ji = i0, i1
634	zau = MIN( 1._wp, zbetdo(ji,jj,jk), zbetup(ji+1,jj,jk) )
635	zbu = MIN( 1._wp, zbetup(ji,jj,jk), zbetdo(ji+1,jj,jk) )
636	zcu = ( 0.5 + SIGN( 0.5 , paa(ji,jj,jk) ) )
637	paa(ji,jj,jk) = paa(ji,jj,jk) * ( zcu * zau + ( 1._wp - zcu) * zbu )
638
639	zav = MIN( 1._wp, zbetdo(ji,jj,jk), zbetup(ji,jj+1,jk) )
640	zbv = MIN( 1._wp, zbetup(ji,jj,jk), zbetdo(ji,jj+1,jk) )
641	zcv = ( 0.5 + SIGN( 0.5 , pbb(ji,jj,jk) ) )
642	pbb(ji,jj,jk) = pbb(ji,jj,jk) * ( zcv * zav + ( 1._wp - zcv) * zbv )
643
644	! monotonic flux in the k direction, i.e. pcc
645	! -------------------------------------------
646	za = MIN( 1., zbetdo(ji,jj,jk+1), zbetup(ji,jj,jk) )
647	zb = MIN( 1., zbetup(ji,jj,jk+1), zbetdo(ji,jj,jk) )
648	zc = ( 0.5 + SIGN( 0.5 , pcc(ji,jj,jk+1) ) )
649	pcc(ji,jj,jk+1) = pcc(ji,jj,jk+1) * ( zc * za + ( 1._wp - zc) * zb )
650	END DO
651	END DO
652	END DO
653	endif
654	end subroutine
655	#endif
656	END SUBROUTINE nonosc
657
658
659	SUBROUTINE interp_4th_cpt_org( pt_in, pt_out )
660	!!----------------------------------------------------------------------
661	!! * ROUTINE interp_4th_cpt_org *
662	!!
663	!! ** Purpose : Compute the interpolation of tracer at w-point
664	!!
665	!! ** Method : 4th order compact interpolation
666	!!----------------------------------------------------------------------
667	REAL(wp),DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pt_in ! now tracer fields
668	REAL(wp),DIMENSION(jpi,jpj,jpk), INTENT( out) :: pt_out ! now tracer field interpolated at w-pts
669	!
670	INTEGER :: ji, jj, jk ! dummy loop integers
671	REAL(wp),DIMENSION(jpi,jpj,jpk) :: zwd, zwi, zws, zwrm, zwt
672	!!----------------------------------------------------------------------
673
674	DO jk = 3, jpkm1 !== build the three diagonal matrix ==!
675	DO jj = 1, jpj
676	DO ji = 1, jpi
677	zwd (ji,jj,jk) = 4._wp
678	zwi (ji,jj,jk) = 1._wp
679	zws (ji,jj,jk) = 1._wp
680	zwrm(ji,jj,jk) = 3._wp * ( pt_in(ji,jj,jk-1) + pt_in(ji,jj,jk) )
681	!
682	IF( tmask(ji,jj,jk+1) == 0._wp) THEN ! Switch to second order centered at bottom
683	zwd (ji,jj,jk) = 1._wp
684	zwi (ji,jj,jk) = 0._wp
685	zws (ji,jj,jk) = 0._wp
686	zwrm(ji,jj,jk) = 0.5 * ( pt_in(ji,jj,jk-1) + pt_in(ji,jj,jk) )
687	ENDIF
688	END DO
689	END DO
690	END DO
691	!
692	jk = 2 ! Switch to second order centered at top
693	DO jj = 1, jpj
694	DO ji = 1, jpi
695	zwd (ji,jj,jk) = 1._wp
696	zwi (ji,jj,jk) = 0._wp
697	zws (ji,jj,jk) = 0._wp
698	zwrm(ji,jj,jk) = 0.5 * ( pt_in(ji,jj,jk-1) + pt_in(ji,jj,jk) )
699	END DO
700	END DO
701	!
702	! !== tridiagonal solve ==!
703	DO jj = 1, jpj ! first recurrence
704	DO ji = 1, jpi
705	zwt(ji,jj,2) = zwd(ji,jj,2)
706	END DO
707	END DO
708	DO jk = 3, jpkm1
709	DO jj = 1, jpj
710	DO ji = 1, jpi
711	zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) /zwt(ji,jj,jk-1)
712	END DO
713	END DO
714	END DO
715	!
716	DO jj = 1, jpj ! second recurrence: Zk = Yk - Ik / Tk-1 Zk-1
717	DO ji = 1, jpi
718	pt_out(ji,jj,2) = zwrm(ji,jj,2)
719	END DO
720	END DO
721	DO jk = 3, jpkm1
722	DO jj = 1, jpj
723	DO ji = 1, jpi
724	pt_out(ji,jj,jk) = zwrm(ji,jj,jk) - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1)
725	END DO
726	END DO
727	END DO
728
729	DO jj = 1, jpj ! third recurrence: Xk = (Zk - Sk Xk+1 ) / Tk
730	DO ji = 1, jpi
731	pt_out(ji,jj,jpkm1) = pt_out(ji,jj,jpkm1) / zwt(ji,jj,jpkm1)
732	END DO
733	END DO
734	DO jk = jpk-2, 2, -1
735	DO jj = 1, jpj
736	DO ji = 1, jpi
737	pt_out(ji,jj,jk) = ( pt_out(ji,jj,jk) - zws(ji,jj,jk) * pt_out(ji,jj,jk+1) ) / zwt(ji,jj,jk)
738	END DO
739	END DO
740	END DO
741	!
742	END SUBROUTINE interp_4th_cpt_org
743
744
745	SUBROUTINE interp_4th_cpt( pt_in, pt_out )
746	!!----------------------------------------------------------------------
747	!! * ROUTINE interp_4th_cpt *
748	!!
749	!! ** Purpose : Compute the interpolation of tracer at w-point
750	!!
751	!! ** Method : 4th order compact interpolation
752	!!----------------------------------------------------------------------
753	REAL(wp),DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pt_in ! field at t-point
754	REAL(wp),DIMENSION(jpi,jpj,jpk), INTENT( out) :: pt_out ! field interpolated at w-point
755	!
756	INTEGER :: ji, jj, jk ! dummy loop integers
757	INTEGER :: ikt, ikb ! local integers
758	REAL(wp),DIMENSION(jpi,jpj,jpk) :: zwd, zwi, zws, zwrm, zwt
759	!!----------------------------------------------------------------------
760	!
761	! !== build the three diagonal matrix & the RHS ==!
762	!
763	DO jk = 3, jpkm1 ! interior (from jk=3 to jpk-1)
764	DO jj = 2, jpjm1
765	DO ji = fs_2, fs_jpim1
766	zwd (ji,jj,jk) = 3._wp * wmask(ji,jj,jk) + 1._wp ! diagonal
767	zwi (ji,jj,jk) = wmask(ji,jj,jk) ! lower diagonal
768	zws (ji,jj,jk) = wmask(ji,jj,jk) ! upper diagonal
769	zwrm(ji,jj,jk) = 3._wp * wmask(ji,jj,jk) & ! RHS
770	& * ( pt_in(ji,jj,jk) + pt_in(ji,jj,jk-1) )
771	END DO
772	END DO
773	END DO
774	!
775	!!gm
776	! SELECT CASE( kbc ) !* boundary condition
777	! CASE( np_NH ) ! Neumann homogeneous at top & bottom
778	! CASE( np_CEN2 ) ! 2nd order centered at top & bottom
779	! END SELECT
780	!!gm
781	!
782	DO jj = 2, jpjm1 ! 2nd order centered at top & bottom
783	DO ji = fs_2, fs_jpim1
784	ikt = mikt(ji,jj) + 1 ! w-point below the 1st wet point
785	ikb = mbkt(ji,jj) ! - above the last wet point
786	!
787	zwd (ji,jj,ikt) = 1._wp ! top
788	zwi (ji,jj,ikt) = 0._wp
789	zws (ji,jj,ikt) = 0._wp
790	zwrm(ji,jj,ikt) = 0.5_wp * ( pt_in(ji,jj,ikt-1) + pt_in(ji,jj,ikt) )
791	!
792	zwd (ji,jj,ikb) = 1._wp ! bottom
793	zwi (ji,jj,ikb) = 0._wp
794	zws (ji,jj,ikb) = 0._wp
795	zwrm(ji,jj,ikb) = 0.5_wp * ( pt_in(ji,jj,ikb-1) + pt_in(ji,jj,ikb) )
796	END DO
797	END DO
798	!
799	! !== tridiagonal solver ==!
800	!
801	DO jj = 2, jpjm1 !* 1st recurrence: Tk = Dk - Ik Sk-1 / Tk-1
802	DO ji = fs_2, fs_jpim1
803	zwt(ji,jj,2) = zwd(ji,jj,2)
804	END DO
805	END DO
806	DO jk = 3, jpkm1
807	DO jj = 2, jpjm1
808	DO ji = fs_2, fs_jpim1
809	zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) /zwt(ji,jj,jk-1)
810	END DO
811	END DO
812	END DO
813	!
814	DO jj = 2, jpjm1 !* 2nd recurrence: Zk = Yk - Ik / Tk-1 Zk-1
815	DO ji = fs_2, fs_jpim1
816	pt_out(ji,jj,2) = zwrm(ji,jj,2)
817	END DO
818	END DO
819	DO jk = 3, jpkm1
820	DO jj = 2, jpjm1
821	DO ji = fs_2, fs_jpim1
822	pt_out(ji,jj,jk) = zwrm(ji,jj,jk) - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1)
823	END DO
824	END DO
825	END DO
826
827	DO jj = 2, jpjm1 !* 3d recurrence: Xk = (Zk - Sk Xk+1 ) / Tk
828	DO ji = fs_2, fs_jpim1
829	pt_out(ji,jj,jpkm1) = pt_out(ji,jj,jpkm1) / zwt(ji,jj,jpkm1)
830	END DO
831	END DO
832	DO jk = jpk-2, 2, -1
833	DO jj = 2, jpjm1
834	DO ji = fs_2, fs_jpim1
835	pt_out(ji,jj,jk) = ( pt_out(ji,jj,jk) - zws(ji,jj,jk) * pt_out(ji,jj,jk+1) ) / zwt(ji,jj,jk)
836	END DO
837	END DO
838	END DO
839	!
840	END SUBROUTINE interp_4th_cpt
841
842
843	SUBROUTINE tridia_solver( pD, pU, pL, pRHS, pt_out , klev )
844	!!----------------------------------------------------------------------
845	!! * ROUTINE tridia_solver *
846	!!
847	!! ** Purpose : solve a symmetric 3diagonal system
848	!!
849	!! ** Method : solve M.t_out = RHS(t) where M is a tri diagonal matrix ( jpk*jpk )
850	!!
851	!! ( D_1 U_1 0 0 0 )( t_1 ) ( RHS_1 )
852	!! ( L_2 D_2 U_2 0 0 )( t_2 ) ( RHS_2 )
853	!! ( 0 L_3 D_3 U_3 0 )( t_3 ) = ( RHS_3 )
854	!! ( ... )( ... ) ( ... )
855	!! ( 0 0 0 L_k D_k )( t_k ) ( RHS_k )
856	!!
857	!! M is decomposed in the product of an upper and lower triangular matrix.
858	!! The tri-diagonals matrix is given as input 3D arrays: pD, pU, pL
859	!! (i.e. the Diagonal, the Upper diagonal, and the Lower diagonal).
860	!! The solution is pta.
861	!! The 3d array zwt is used as a work space array.
862	!!----------------------------------------------------------------------
863	REAL(wp),DIMENSION(:,:,:), INTENT(in ) :: pD, pU, PL ! 3-diagonal matrix
864	REAL(wp),DIMENSION(:,:,:), INTENT(in ) :: pRHS ! Right-Hand-Side
865	REAL(wp),DIMENSION(:,:,:), INTENT( out) :: pt_out !!gm field at level=F(klev)
866	INTEGER , INTENT(in ) :: klev ! =1 pt_out at w-level
867	! ! =0 pt at t-level
868	INTEGER :: ji, jj, jk ! dummy loop integers
869	INTEGER :: kstart ! local indices
870	REAL(wp),DIMENSION(jpi,jpj,jpk) :: zwt ! 3D work array
871	!!----------------------------------------------------------------------
872	!
873	kstart = 1 + klev
874	!
875	DO jj = 2, jpjm1 !* 1st recurrence: Tk = Dk - Ik Sk-1 / Tk-1
876	DO ji = fs_2, fs_jpim1
877	zwt(ji,jj,kstart) = pD(ji,jj,kstart)
878	END DO
879	END DO
880	DO jk = kstart+1, jpkm1
881	DO jj = 2, jpjm1
882	DO ji = fs_2, fs_jpim1
883	zwt(ji,jj,jk) = pD(ji,jj,jk) - pL(ji,jj,jk) * pU(ji,jj,jk-1) /zwt(ji,jj,jk-1)
884	END DO
885	END DO
886	END DO
887	!
888	DO jj = 2, jpjm1 !* 2nd recurrence: Zk = Yk - Ik / Tk-1 Zk-1
889	DO ji = fs_2, fs_jpim1
890	pt_out(ji,jj,kstart) = pRHS(ji,jj,kstart)
891	END DO
892	END DO
893	DO jk = kstart+1, jpkm1
894	DO jj = 2, jpjm1
895	DO ji = fs_2, fs_jpim1
896	pt_out(ji,jj,jk) = pRHS(ji,jj,jk) - pL(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1)
897	END DO
898	END DO
899	END DO
900
901	DO jj = 2, jpjm1 !* 3d recurrence: Xk = (Zk - Sk Xk+1 ) / Tk
902	DO ji = fs_2, fs_jpim1
903	pt_out(ji,jj,jpkm1) = pt_out(ji,jj,jpkm1) / zwt(ji,jj,jpkm1)
904	END DO
905	END DO
906	DO jk = jpk-2, kstart, -1
907	DO jj = 2, jpjm1
908	DO ji = fs_2, fs_jpim1
909	pt_out(ji,jj,jk) = ( pt_out(ji,jj,jk) - pU(ji,jj,jk) * pt_out(ji,jj,jk+1) ) / zwt(ji,jj,jk)
910	END DO
911	END DO
912	END DO
913	!
914	END SUBROUTINE tridia_solver
915
916	!!======================================================================
917	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/2018/dev_r9759_HPC09_ESIWACE/src/OCE/TRA/traadv_fct.F90 @ 10136

Download in other formats: