integer,allocatable :: ni(:),nj(:),nk(:)              ! for multiblocks

- This is allocated to nblocks (1 in the use case)
- So ni(1) = nj(1) = nk(1) = 128 (in current use case)


       real(dp),dimension(NImax,NJmax,NKmax,nblocks,nvars) :: PHI, PHID
       real(dp),dimension(NImax) :: F,xout
       real(dp),dimension(NImax) :: AMlocal,APlocal,AClocal

- NImax,NJmax,NKmax=128
- nblocks,nvars < 10 (1, 6 in the use case)

---------

!$acc parallel private(F,xout,APlocal,AMlocal,AClocal)
!$acc loop independent seq
       do prim=1,nvars
!$acc loop independent seq
        do nbl=1,nblocks
!$acc loop independent vector
         do k=1,128 !NK(nbl)
!$acc loop
          do j=1,128 !NJ(nbl)
            F(:)=0.d0
            xout(:)=0.d0
           do i=3,NI(nbl)-2
             PHID(i,j,k,nbl,prim)=(b4*(PHI(i+2,j,k,nbl,prim)- &
                                  PHI(i-2,j,k,nbl,prim)))+ &
             (a2*(PHI(i+1,j,k,nbl,prim)-PHI(i-1,j,k,nbl,prim)))
           enddo
             PHID(1,j,k,nbl,prim)=a_1*PHI(1,j,k,nbl,prim) &
                                 +b_1*PHI(2,j,k,nbl,prim) &
                                 +c_1*PHI(3,j,k,nbl,prim) &
                                 +d_1*PHI(4,j,k,nbl,prim) &
                                 +e_1*PHI(5,j,k,nbl,prim) &
                                 +f_1*PHI(6,j,k,nbl,prim)

             PHID(2,j,k,nbl,prim)=a_2*PHI(1,j,k,nbl,prim) &
                                 +b_2*PHI(2,j,k,nbl,prim) &
                                 +c_2*PHI(3,j,k,nbl,prim) &
                                 +d_2*PHI(4,j,k,nbl,prim) &
                                 +e_2*PHI(5,j,k,nbl,prim) &
                                 +f_2*PHI(6,j,k,nbl,prim)

          PHID(NI(nbl)-1,j,k,nbl,prim)=a_M*PHI(NI(nbl),j,k,nbl,prim) &
                                    +b_M*PHI(NI(nbl)-1,j,k,nbl,prim) &
                                    +c_M*PHI(NI(nbl)-2,j,k,nbl,prim) &
                                    +d_M*PHI(NI(nbl)-3,j,k,nbl,prim) &
                                    +e_M*PHI(NI(nbl)-4,j,k,nbl,prim) &
                                    +f_M*PHI(NI(nbl)-5,j,k,nbl,prim)

            PHID(NI(nbl),j,k,nbl,prim)=a_N*PHI(NI(nbl),j,k,nbl,prim) &
                                    +b_N*PHI(NI(nbl)-1,j,k,nbl,prim) &
                                    +c_N*PHI(NI(nbl)-2,j,k,nbl,prim) &
                                    +d_N*PHI(NI(nbl)-3,j,k,nbl,prim) &
                                    +e_N*PHI(NI(nbl)-4,j,k,nbl,prim) &
                                    +f_N*PHI(NI(nbl)-5,j,k,nbl,prim)
!---------- build the coefficients of the diagonals and the right hand side of the tri-diagonal matrix ------------------

              F(1:NI(nbl))=PHID(1:NI(nbl),j,k,nbl,prim)
              APlocal(1)=alpha_1
              AMlocal(1)=0.d0
              APlocal(2)=alpha_22
              AMlocal(2)=alpha_21
              APlocal(3:NI(nbl)-2)=alpha
              AMlocal(3:NI(nbl)-2)=alpha
              APlocal(NI(nbl)-1)=alpha_M2
              AMlocal(NI(nbl)-1)=alpha_M1
              APlocal(NI(nbl))=0.d0
              AMlocal(NI(nbl))=alpha_N
              AClocal(1:NI(nbl))=1.d0

              do i=2,NI(nbl)
                AClocal(i) = AClocal(i) - &
                              AMlocal(i)/AClocal(i-1)*APlocal(i-1)
                F(i) = F(i) - AMlocal(i)/AClocal(i-1)*F(i-1)
              enddo

              xout(NI(nbl)) = F(NI(nbl))/AClocal(NI(nbl))
              do i=NI(nbl)-1,1,-1
                xout(i) = (F(i)-APlocal(i)*xout(i+1))/AClocal(i)
              enddo

             PHID(1:NI(nbl),j,k,nbl,prim)=xout(1:NI(nbl))
          enddo
!$acc end loop
         enddo
!$acc end loop
        enddo
!$acc end loop
       enddo
!$acc end loop


!$acc end parallel