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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

traadv_fct.F90 in NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/OCE/TRA – NEMO

Context Navigation

source: NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/OCE/TRA/traadv_fct.F90 @ 13553

Last change on this file since 13553 was 13553, checked in by hadcv, 4 years ago
Merge in trunk up to [13550]
Property svn:keywords set to `Id`
File size: 35.2 KB

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

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

Download in other formats: