diff --git a/src/ZSOLVER/tridiag_thermo.f90 b/src/ZSOLVER/tridiag_thermo.f90
new file mode 100644
index 0000000000000000000000000000000000000000..95e86abe02d31200869309a613f7b86efe447094
--- /dev/null
+++ b/src/ZSOLVER/tridiag_thermo.f90
@@ -0,0 +1,464 @@
+!MNH_LIC Copyright 2003-2020 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
+!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
+!MNH_LIC for details. version 1.
+!-----------------------------------------------------------------
+! ###################
+ MODULE MODI_TRIDIAG_THERMO
+! ###################
+!
+INTERFACE
+!
+ SUBROUTINE TRIDIAG_THERMO(KKA,KKU,KKL,PVARM,PF,PDFDDTDZ,PTSTEP,PIMPL, &
+ PDZZ,PRHODJ,PVARP )
+!
+INTEGER, INTENT(IN) :: KKA !near ground array index
+INTEGER, INTENT(IN) :: KKU !uppest atmosphere array index
+INTEGER, INTENT(IN) :: KKL !vert. levels type 1=MNH -1=AR
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PVARM ! variable at t-1 at mass point
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PF ! flux in dT/dt=-dF/dz at flux point
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDFDDTDZ! dF/d(dT/dz) at flux point
+REAL, INTENT(IN) :: PTSTEP ! Double time step
+REAL, INTENT(IN) :: PIMPL ! implicit weight
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZZ ! Dz at flux point
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODJ ! (dry rho)*J at mass point
+!
+REAL, DIMENSION(:,:,:), INTENT(OUT):: PVARP ! variable at t+1 at mass point
+!
+END SUBROUTINE TRIDIAG_THERMO
+!
+END INTERFACE
+!
+END MODULE MODI_TRIDIAG_THERMO
+!
+! #################################################
+ SUBROUTINE TRIDIAG_THERMO(KKA,KKU,KKL,PVARM,PF,PDFDDTDZ,PTSTEP,PIMPL, &
+ PDZZ,PRHODJ,PVARP )
+! #################################################
+!
+!
+!!**** *TRIDIAG_THERMO* - routine to solve a time implicit scheme
+!!
+!!
+!! PURPOSE
+!! -------
+! The purpose of this routine is to give a field PVARP at t+1, by
+! solving an implicit TRIDIAGonal system obtained by the
+! discretization of the vertical turbulent diffusion. It should be noted
+! that the degree of implicitness can be varied (PIMPL parameter) and that
+! the function of F(dT/dz) must have been linearized.
+! PVARP is localized at a mass point.
+!
+!!** METHOD
+!! ------
+!!
+!! [T(+) - T(-)]/2Dt = -d{ F + dF/d(dT/dz) * [impl*dT/dz(+) + expl* dT/dz(-)] }/dz
+!!
+!! It is discretized as follows:
+!!
+!! PRHODJ(k)*PVARP(k)/PTSTEP
+!! =
+!! PRHODJ(k)*PVARM(k)/PTSTEP
+!! - (PRHODJ(k+1)+PRHODJ(k) )/2. * PF(k+1)/PDZZ(k+1)
+!! + (PRHODJ(k) +PRHODJ(k-1))/2. * PF(k) /PDZZ(k)
+!! + (PRHODJ(k+1)+PRHODJ(k) )/2. * ZEXPL* PDFDDTDZ(k+1) * PVARM(k+1)/PDZZ(k+1)**2
+!! - (PRHODJ(k+1)+PRHODJ(k) )/2. * PIMPL* PDFDDTDZ(k+1) * PVARP(k+1)/PDZZ(k+1)**2
+!! - (PRHODJ(k+1)+PRHODJ(k) )/2. * ZEXPL* PDFDDTDZ(k+1) * PVARM(k) /PDZZ(k+1)**2
+!! + (PRHODJ(k+1)+PRHODJ(k) )/2. * PIMPL* PDFDDTDZ(k+1) * PVARP(k) /PDZZ(k+1)**2
+!! - (PRHODJ(k) +PRHODJ(k-1))/2. * ZEXPL* PDFDDTDZ(k) * PVARM(k) /PDZZ(k)**2
+!! + (PRHODJ(k) +PRHODJ(k-1))/2. * PIMPL* PDFDDTDZ(k) * PVARP(k) /PDZZ(k)**2
+!! + (PRHODJ(k) +PRHODJ(k-1))/2. * ZEXPL* PDFDDTDZ(k) * PVARM(k-1)/PDZZ(k)**2
+!! - (PRHODJ(k) +PRHODJ(k-1))/2. * PIMPL* PDFDDTDZ(k) * PVARP(k-1)/PDZZ(k)**2
+!!
+!!
+!! The system to solve is:
+!!
+!! A*PVARP(k-1) + B*PVARP(k) + C*PVARP(k+1) = Y(k)
+!!
+!!
+!! The RHS of the linear system in PVARP writes:
+!!
+!! y(k) = PRHODJ(k)*PVARM(k)/PTSTEP
+!! - (PRHODJ(k+1)+PRHODJ(k) )/2. * PF(k+1)/PDZZ(k+1)
+!! + (PRHODJ(k) +PRHODJ(k-1))/2. * PF(k) /PDZZ(k)
+!! + (PRHODJ(k+1)+PRHODJ(k) )/2. * ZEXPL* PDFDDTDZ(k+1) * PVARM(k+1)/PDZZ(k+1)**2
+!! - (PRHODJ(k+1)+PRHODJ(k) )/2. * ZEXPL* PDFDDTDZ(k+1) * PVARM(k) /PDZZ(k+1)**2
+!! - (PRHODJ(k) +PRHODJ(k-1))/2. * ZEXPL* PDFDDTDZ(k) * PVARM(k) /PDZZ(k)**2
+!! + (PRHODJ(k) +PRHODJ(k-1))/2. * ZEXPL* PDFDDTDZ(k) * PVARM(k-1)/PDZZ(k)**2
+!!
+!!
+!! Then, the classical TRIDIAGonal algorithm is used to invert the
+!! implicit operator. Its matrix is given by:
+!!
+!! ( b(ikb) c(ikb) 0 0 0 0 0 0 )
+!! ( 0 a(ikb+1) b(ikb+1) c(ikb+1) 0 ... 0 0 0 )
+!! ( 0 0 a(ikb+2) b(ikb+2) c(ikb+2). 0 0 0 )
+!! .......................................................................
+!! ( 0 ... 0 a(k) b(k) c(k) 0 ... 0 0 )
+!! .......................................................................
+!! ( 0 0 0 0 0 ...a(ike-1) b(ike-1) c(ike-1))
+!! ( 0 0 0 0 0 ... 0 a(ike) b(ike) )
+!!
+!! ikb and ike represent the first and the last inner mass levels of the
+!! model. The coefficients are:
+!!
+!! a(k) = + (PRHODJ(k) +PRHODJ(k-1))/2. * PIMPL* PDFDDTDZ(k) /PDZZ(k)**2
+!! b(k) = PRHODJ(k) / PTSTEP
+!! - (PRHODJ(k+1)+PRHODJ(k) )/2. * PIMPL* PDFDDTDZ(k+1)/PDZZ(k+1)**2
+!! - (PRHODJ(k) +PRHODJ(k-1))/2. * PIMPL* PDFDDTDZ(k) /PDZZ(k)**2
+!! c(k) = + (PRHODJ(k+1)+PRHODJ(k) )/2. * PIMPL* PDFDDTDZ(k+1)/PDZZ(k+1)**2
+!!
+!! for all k /= ikb or ike
+!!
+!!
+!! b(ikb) = PRHODJ(ikb) / PTSTEP
+!! -(PRHODJ(ikb+1)+PRHODJ(ikb))/2.*PIMPL*PDFDDTDZ(ikb+1)/PDZZ(ikb+1)**2
+!! c(ikb) = +(PRHODJ(ikb+1)+PRHODJ(ikb))/2.*PIMPL*PDFDDTDZ(ikb+1)/PDZZ(ikb+1)**2
+!!
+!! b(ike) = PRHODJ(ike) / PTSTEP
+!! -(PRHODJ(ike)+PRHODJ(ike-1))/2.*PIMPL*PDFDDTDZ(ike)/PDZZ(ike)**2
+!! a(ike) = +(PRHODJ(ike)+PRHODJ(ike-1))/2.*PIMPL*PDFDDTDZ(ike)/PDZZ(ike)**2
+!!
+!!
+!! EXTERNAL
+!! --------
+!!
+!! NONE
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!!
+!! REFERENCE
+!! ---------
+!! Press et al: Numerical recipes (1986) Cambridge Univ. Press
+!!
+!! AUTHOR
+!! ------
+!! V. Masson * Meteo-France *
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 04/2003 (from tridiag.f90)
+!! M.Moge 04/2016 Use openACC directives to port the TURB part of Meso-NH on GPU
+!! ---------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+!
+USE MODD_PARAMETERS, ONLY : JPVEXT_TURB
+
+use mode_mppdb
+
+#ifndef MNH_OPENACC
+USE MODI_SHUMAN
+#else
+USE MODI_SHUMAN_DEVICE
+#endif
+#ifdef MNH_BITREP
+USE MODI_BITREP
+#endif
+!
+#ifdef MNH_OPENACC
+USE MODE_MNH_ZWORK , ONLY : MNH_ALLOCATE_ZT3D , MNH_REL_ZT3D , MNH_ALLOCATE_ZT2D, &
+ MNH_CHECK_IN_ZT3D,MNH_CHECK_OUT_ZT3D
+#endif
+!
+IMPLICIT NONE
+!
+!
+!* 0.1 declarations of arguments
+!
+INTEGER, INTENT(IN) :: KKA !near ground array index
+INTEGER, INTENT(IN) :: KKU !uppest atmosphere array index
+INTEGER, INTENT(IN) :: KKL !vert. levels type 1=MNH -1=ARO
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PVARM ! variable at t-1 at mass point
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PF ! flux in dT/dt=-dF/dz at flux point
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDFDDTDZ! dF/d(dT/dz) at flux point
+REAL, INTENT(IN) :: PTSTEP ! Double time step
+REAL, INTENT(IN) :: PIMPL ! implicit weight
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZZ ! Dz at flux point
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODJ ! (dry rho)*J at mass point
+!
+REAL, DIMENSION(:,:,:), INTENT(OUT):: PVARP ! variable at t+1 at mass point
+!
+!* 0.2 declarations of local variables
+!
+REAL, DIMENSION(:,:,:), pointer , contiguous :: ZRHODJ_DFDDTDZ_O_DZ2
+REAL, DIMENSION(:,:,:), pointer , contiguous :: ZMZM_RHODJ
+REAL, DIMENSION(:,:,:), pointer , contiguous :: ZA, ZB, ZC
+REAL, DIMENSION(:,:,:), pointer , contiguous :: ZY ,ZGAM
+ ! RHS of the equation, 3D work array
+INTEGER :: IZRHODJ_DFDDTDZ_O_DZ2,IZMZM_RHODJ,IZA,IZB,IZC,IZY,IZGAM
+REAL, DIMENSION(:,:), pointer , contiguous :: ZBET
+INTEGER :: IZBET
+ ! 2D work array
+INTEGER :: JI,JJ,JK ! loop counter
+INTEGER :: IKB,IKE ! inner vertical limits
+INTEGER :: IKT ! array size in k direction
+INTEGER :: IKTB,IKTE ! start, end of k loops in physical domain
+!
+!
+#ifdef MNH_OPENACC
+REAL, DIMENSION(:,:,:), pointer , contiguous :: ZTMP1_DEVICE
+INTEGER :: IZTMP1_DEVICE
+#endif
+
+INTEGER :: JIU,JJU,JKU
+
+! ---------------------------------------------------------------------------
+
+!$acc data present( PVARM, PF, PDFDDTDZ, PDZZ, PRHODJ, PVARP )
+
+if ( mppdb_initialized ) then
+ !Check all in arrays
+ call Mppdb_check( pvarm, "Tridiag_thermo beg:pvarm" )
+ call Mppdb_check( pf, "Tridiag_thermo beg:pf" )
+ call Mppdb_check( pdfddtdz, "Tridiag_thermo beg:pdfddtdz" )
+ call Mppdb_check( pdzz, "Tridiag_thermo beg:pdzz" )
+ call Mppdb_check( prhodj, "Tridiag_thermo beg:prhodj" )
+end if
+
+JIU = size( pvarm, 1 )
+JJU = size( pvarm, 2 )
+JKU = size( pvarm, 3 )
+
+#ifndef MNH_OPENACC
+allocate( zrhodj_dfddtdz_o_dz2(JIU,JJU,JKU ) )
+allocate( zmzm_rhodj (JIU,JJU,JKU ) )
+allocate( za (JIU,JJU,JKU ) )
+allocate( zb (JIU,JJU,JKU ) )
+allocate( zc (JIU,JJU,JKU ) )
+allocate( zy (JIU,JJU,JKU ) )
+allocate( zgam (JIU,JJU,JKU ) )
+allocate( zbet (JIU,JJU ) )
+#else
+CALL MNH_CHECK_IN_ZT3D("TRIDIAG_THERMO")
+izrhodj_dfddtdz_o_dz2 = MNH_ALLOCATE_ZT3D( zrhodj_dfddtdz_o_dz2,JIU,JJU,JKU )
+izmzm_rhodj = MNH_ALLOCATE_ZT3D( zmzm_rhodj ,JIU,JJU,JKU )
+iza = MNH_ALLOCATE_ZT3D( za ,JIU,JJU,JKU )
+izb = MNH_ALLOCATE_ZT3D( zb ,JIU,JJU,JKU )
+izc = MNH_ALLOCATE_ZT3D( zc ,JIU,JJU,JKU )
+izy = MNH_ALLOCATE_ZT3D( zy ,JIU,JJU,JKU )
+izgam = MNH_ALLOCATE_ZT3D( zgam ,JIU,JJU,JKU )
+izbet = MNH_ALLOCATE_ZT2D( zbet ,JIU,JJU )
+#endif
+
+#ifdef MNH_OPENACC
+iztmp1_device = MNH_ALLOCATE_ZT3D( ztmp1_device,JIU,JJU,JKU )
+#endif
+
+!$acc data present( zrhodj_dfddtdz_o_dz2, zmzm_rhodj, za, zb, zc, zy, zgam, zbet, ztmp1_device )
+!
+!* 1. Preliminaries
+! -------------
+!
+IKT=SIZE(PVARM,3)
+IKTB=1+JPVEXT_TURB
+IKTE=IKT-JPVEXT_TURB
+IKB=KKA+JPVEXT_TURB*KKL
+IKE=KKU-JPVEXT_TURB*KKL
+!
+#ifndef MNH_OPENACC
+ZMZM_RHODJ = MZM(PRHODJ)
+#else
+CALL MZM_DEVICE(PRHODJ,ZMZM_RHODJ)
+#endif
+!$acc kernels ! async
+#ifndef MNH_BITREP
+ZRHODJ_DFDDTDZ_O_DZ2 = ZMZM_RHODJ*PDFDDTDZ/PDZZ**2
+#else
+!$acc loop independent collapse(3)
+DO CONCURRENT ( JI=1:JIU,JJ=1:JJU,JK=1:JKU)
+ ZRHODJ_DFDDTDZ_O_DZ2(JI,JJ,JK) = ZMZM_RHODJ(JI,JJ,JK)*PDFDDTDZ(JI,JJ,JK)/BR_P2(PDZZ(JI,JJ,JK))
+END DO !CONCURRENT
+#endif
+!$acc end kernels
+!
+!$acc kernels ! async
+ZA=0.
+ZB=0.
+ZC=0.
+ZY=0.
+!$acc end kernels
+! acc wait
+!
+!
+!* 2. COMPUTE THE RIGHT HAND SIDE
+! ---------------------------
+!
+!$acc kernels ! async
+!$acc loop independent collapse(2)
+DO CONCURRENT ( JI=1:JIU,JJ=1:JJU)
+ZY(JI,JJ,IKB) = PRHODJ(JI,JJ,IKB)*PVARM(JI,JJ,IKB)/PTSTEP &
+ - ZMZM_RHODJ(JI,JJ,IKB+KKL) * PF(JI,JJ,IKB+KKL)/PDZZ(JI,JJ,IKB+KKL) &
+ + ZMZM_RHODJ(JI,JJ,IKB ) * PF(JI,JJ,IKB )/PDZZ(JI,JJ,IKB ) &
+ + ZRHODJ_DFDDTDZ_O_DZ2(JI,JJ,IKB+KKL) * PIMPL * PVARM(JI,JJ,IKB+KKL) &
+ - ZRHODJ_DFDDTDZ_O_DZ2(JI,JJ,IKB+KKL) * PIMPL * PVARM(JI,JJ,IKB )
+END DO !CONCURRENT
+!$acc end kernels
+!
+!$acc kernels ! async
+!$acc loop independent collapse(3)
+DO CONCURRENT ( JI=1:JIU,JJ=1:JJU,JK=IKTB+1:IKTE-1)
+ ZY(JI,JJ,JK) = PRHODJ(JI,JJ,JK)*PVARM(JI,JJ,JK)/PTSTEP &
+ - ZMZM_RHODJ(JI,JJ,JK+KKL) * PF(JI,JJ,JK+KKL)/PDZZ(JI,JJ,JK+KKL) &
+ + ZMZM_RHODJ(JI,JJ,JK ) * PF(JI,JJ,JK )/PDZZ(JI,JJ,JK ) &
+ + ZRHODJ_DFDDTDZ_O_DZ2(JI,JJ,JK+KKL) * PIMPL * PVARM(JI,JJ,JK+KKL) &
+ - ZRHODJ_DFDDTDZ_O_DZ2(JI,JJ,JK+KKL) * PIMPL * PVARM(JI,JJ,JK ) &
+ - ZRHODJ_DFDDTDZ_O_DZ2(JI,JJ,JK ) * PIMPL * PVARM(JI,JJ,JK ) &
+ + ZRHODJ_DFDDTDZ_O_DZ2(JI,JJ,JK ) * PIMPL * PVARM(JI,JJ,JK-KKL)
+END DO !CONCURRENT
+!$acc end kernels
+!
+!$acc kernels ! async
+!$acc loop independent collapse(2)
+DO CONCURRENT ( JI=1:JIU,JJ=1:JJU)
+ZY(JI,JJ,IKE) = PRHODJ(JI,JJ,IKE)*PVARM(JI,JJ,IKE)/PTSTEP &
+ - ZMZM_RHODJ(JI,JJ,IKE+KKL) * PF(JI,JJ,IKE+KKL)/PDZZ(JI,JJ,IKE+KKL) &
+ + ZMZM_RHODJ(JI,JJ,IKE ) * PF(JI,JJ,IKE )/PDZZ(JI,JJ,IKE ) &
+ - ZRHODJ_DFDDTDZ_O_DZ2(JI,JJ,IKE ) * PIMPL * PVARM(JI,JJ,IKE ) &
+ + ZRHODJ_DFDDTDZ_O_DZ2(JI,JJ,IKE ) * PIMPL * PVARM(JI,JJ,IKE-KKL)
+END DO !CONCURRENT
+!$acc end kernels
+!
+! acc wait
+!
+!* 3. INVERSION OF THE TRIDIAGONAL SYSTEM
+! -----------------------------------
+!
+IF ( PIMPL > 1.E-10 ) THEN
+!
+!* 3.1 arrays A, B, C
+! --------------
+!
+!$acc kernels ! async
+!$acc loop independent collapse(2)
+DO CONCURRENT ( JI=1:JIU,JJ=1:JJU)
+ ZB(JI,JJ,IKB) = PRHODJ(JI,JJ,IKB)/PTSTEP &
+ - ZRHODJ_DFDDTDZ_O_DZ2(JI,JJ,IKB+KKL) * PIMPL
+END DO !CONCURRENT
+!$acc end kernels
+!
+!$acc kernels ! async
+!$acc loop independent collapse(2)
+DO CONCURRENT ( JI=1:JIU,JJ=1:JJU)
+ ZC(JI,JJ,IKB) = ZRHODJ_DFDDTDZ_O_DZ2(JI,JJ,IKB+KKL) * PIMPL
+END DO !CONCURRENT
+!$acc end kernels
+!
+!$acc kernels ! async
+!$acc loop independent collapse(3)
+DO CONCURRENT ( JI=1:JIU,JJ=1:JJU,JK=IKTB+1:IKTE-1)
+ ZA(JI,JJ,JK) = ZRHODJ_DFDDTDZ_O_DZ2(JI,JJ,JK) * PIMPL
+ ZB(JI,JJ,JK) = PRHODJ(JI,JJ,JK)/PTSTEP &
+ - ZRHODJ_DFDDTDZ_O_DZ2(JI,JJ,JK+KKL) * PIMPL &
+ - ZRHODJ_DFDDTDZ_O_DZ2(JI,JJ,JK) * PIMPL
+ ZC(JI,JJ,JK) = ZRHODJ_DFDDTDZ_O_DZ2(JI,JJ,JK+KKL) * PIMPL
+END DO !CONCURRENT
+!$acc end kernels
+!
+!$acc kernels ! async
+!$acc loop independent collapse(2)
+DO CONCURRENT ( JI=1:JIU,JJ=1:JJU)
+ ZA(JI,JJ,IKE) = ZRHODJ_DFDDTDZ_O_DZ2(JI,JJ,IKE ) * PIMPL
+ ZB(JI,JJ,IKE) = PRHODJ(JI,JJ,IKE)/PTSTEP &
+ - ZRHODJ_DFDDTDZ_O_DZ2(JI,JJ,IKE ) * PIMPL
+END DO !CONCURRENT
+!$acc end kernels
+!
+! acc wait
+!
+!
+!* 3.2 going up
+! --------
+!
+!$acc kernels
+!$acc loop independent collapse(2)
+DO CONCURRENT ( JI=1:JIU,JJ=1:JJU)
+ ZBET(JI,JJ) = ZB(JI,JJ,IKB) ! bet = b(ikb)
+ PVARP(JI,JJ,IKB) = ZY(JI,JJ,IKB) / ZBET(JI,JJ)
+END DO !CONCURRENT
+!
+!$acc loop seq
+DO JK = IKB+KKL,IKE-KKL,KKL
+ !$acc loop independent collapse(2)
+ ! acc loop gang, vector collapse(2)
+ DO CONCURRENT ( JI=1:JIU,JJ=1:JJU)
+ !DO JJ=1,JJU
+ ! DO JI=1,JIU
+ ZGAM(JI,JJ,JK) = ZC(JI,JJ,JK-KKL) / ZBET(JI,JJ)
+ ! gam(k) = c(k-1) / bet
+ ZBET(JI,JJ) = ZB(JI,JJ,JK) - ZA(JI,JJ,JK) * ZGAM(JI,JJ,JK)
+ ! bet = b(k) - a(k)* gam(k)
+ PVARP(JI,JJ,JK)= ( ZY(JI,JJ,JK) - ZA(JI,JJ,JK) * PVARP(JI,JJ,JK-KKL) ) / ZBET(JI,JJ)
+ ! res(k) = (y(k) -a(k)*res(k-1))/ bet
+ ! END DO
+ !END DO
+ END DO !CONCURRENT
+END DO
+! special treatment for the last level
+!$acc loop independent collapse(2)
+DO CONCURRENT ( JI=1:JIU,JJ=1:JJU)
+ ZGAM(JI,JJ,IKE) = ZC(JI,JJ,IKE-KKL) / ZBET(JI,JJ)
+ ! gam(k) = c(k-1) / bet
+ ZBET(JI,JJ) = ZB(JI,JJ,IKE) - ZA(JI,JJ,IKE) * ZGAM(JI,JJ,IKE)
+ ! bet = b(k) - a(k)* gam(k)
+ PVARP(JI,JJ,IKE)= ( ZY(JI,JJ,IKE) - ZA(JI,JJ,IKE) * PVARP(JI,JJ,IKE-KKL) ) / ZBET(JI,JJ)
+ ! res(k) = (y(k) -a(k)*res(k-1))/ bet
+END DO !CONCURRENT
+!
+!* 3.3 going down
+! ----------
+!
+!$acc loop seq
+DO JK = IKE-KKL,IKB,-1*KKL
+ !$acc loop independent collapse(2)
+ ! acc loop gang, vector collapse(2)
+ DO CONCURRENT ( JI=1:JIU,JJ=1:JJU)
+ PVARP(JI,JJ,JK) = PVARP(JI,JJ,JK) - ZGAM(JI,JJ,JK+KKL) * PVARP(JI,JJ,JK+KKL)
+ END DO !CONCURRENT
+END DO
+!$acc end kernels
+!
+ELSE
+!
+!$acc kernels
+!$acc loop independent collapse(3)
+DO CONCURRENT ( JI=1:JIU,JJ=1:JJU,JK=IKTB:IKTE)
+ PVARP(JI,JJ,JK) = ZY(JI,JJ,JK) * PTSTEP / PRHODJ(JI,JJ,JK)
+END DO !CONCURRENT
+!$acc end kernels
+!
+END IF
+!
+!
+!* 4. FILL THE UPPER AND LOWER EXTERNAL VALUES
+! ----------------------------------------
+!
+!$acc kernels
+!$acc loop independent collapse(2)
+DO CONCURRENT ( JI=1:JIU,JJ=1:JJU)
+ PVARP(JI,JJ,KKA)=PVARP(JI,JJ,IKB)
+ PVARP(JI,JJ,KKU)=PVARP(JI,JJ,IKE)
+END DO !CONCURRENT
+!$acc end kernels
+
+if ( mppdb_initialized ) then
+ !Check all out arrays
+ call Mppdb_check( pvarp, "Tridiag_thermo end:pvarp" )
+end if
+
+!$acc end data
+
+#ifndef MNH_OPENACC
+deallocate (zrhodj_dfddtdz_o_dz2,zmzm_rhodj,za,zb,zc,zy,zgam,zbet)
+#else
+CALL MNH_REL_ZT3D(IZRHODJ_DFDDTDZ_O_DZ2,IZMZM_RHODJ,IZA,IZB,IZC,IZY,IZGAM,&
+ IZBET,iztmp1_device)
+CALL MNH_CHECK_OUT_ZT3D("TRIDIAG_THERMO")
+#endif
+
+!$acc end data
+
+!-------------------------------------------------------------------------------
+!
+END SUBROUTINE TRIDIAG_THERMO