text:
module mass_flux_utilities_module

  use multifab_module
  use ml_layout_module
  use probin_common_module, only: k_B, molmass, rhobar
  use probin_multispecies_module, only: nspecies, use_lapack, fraction_tolerance, &
       is_ideal_mixture, inverse_type, is_nonisothermal, chi_iterations
  use matrix_utilities
  use F95_LAPACK

  implicit none

  private

  public :: correct_rho_with_drho, &
            convert_cons_to_prim, &
            compute_rhotot, &
            compute_Gama, &
            compute_chi, &
            compute_rhoWchi, &
            compute_Lonsager, &
            compute_baro_coef

contains

  subroutine correct_rho_with_drho(mla,rho,drho)

   type(ml_layout), intent(in   )  :: mla
   type(multifab) , intent(inout)  :: rho(:)
   type(multifab) , intent(inout)  :: drho(:)

   ! local variables
   integer :: lo(rho(1)%dim), hi(rho(1)%dim)
   integer :: n,i,ng_1,ng_2,dm,nlevs

   ! pointer for rho(nspecies), rhotot(1), molarconc(nspecies)
   real(kind=dp_t), pointer        :: dp1(:,:,:,:)  ! for rho
   real(kind=dp_t), pointer        :: dp2(:,:,:,:)  ! for drho

   type(mfiter) :: mfi
   type(box) :: tilebox
   integer :: tlo(mla%dim), thi(mla%dim)

   type(bl_prof_timer), save :: bpt

   call build(bpt,"correct_rho_with_drho")

   dm    = mla%dim     ! dimensionality
   ng_1  = rho(1)%ng   ! number of ghost cells
   ng_2  = drho(1)%ng
   nlevs = mla%nlevel  ! number of levels

   !$omp parallel private(n,i,mfi,tilebox,tlo,thi,dp1,dp2,lo,hi)

    ! loop over all boxes
    do n=1,nlevs
       call mfiter_build(mfi, rho(n), tiling=.true.)

       do while (more_tile(mfi))
          i = get_fab_index(mfi)

          tilebox = get_growntilebox(mfi,rho(n)%ng)
          tlo = lwb(tilebox)
          thi = upb(tilebox)

!       do i=1,nfabs(rho(n))
          dp1 => dataptr(rho(n),i)
          dp2 => dataptr(drho(n),i)
          lo = lwb(get_box(rho(n),i))
          hi = upb(get_box(rho(n),i))
          select case(dm)
          case (2)
             call correct_rho_with_drho_2d(dp1(:,:,1,:),dp2(:,:,1,:),ng_1,ng_2,lo,hi,tlo,thi)
          case (3)
             call correct_rho_with_drho_3d(dp1(:,:,:,:),dp2(:,:,:,:),ng_1,ng_2,lo,hi,tlo,thi)
          end select
       end do
    end do
    !$omp end parallel

    call destroy(bpt)

  end subroutine correct_rho_with_drho

  subroutine correct_rho_with_drho_2d(rho,drho,ng_1,ng_2,glo,ghi,tlo,thi)

    integer          :: glo(2), ghi(2), ng_1, ng_2, tlo(2), thi(2)
    real(kind=dp_t)  ::  rho(glo(1)-ng_1:,glo(2)-ng_1:,:) ! density- last dim for #species
    real(kind=dp_t)  :: drho(glo(1)-ng_2:,glo(2)-ng_2:,:) ! total density in each cell

    ! local variables
    integer          :: i,j

    ! for specific box, now start loops over alloted cells
    do j=tlo(2), thi(2)
       do i=tlo(1), thi(1)

         call correct_rho_with_drho_local(rho(i,j,:),drho(i,j,:))

       end do
    end do

  end subroutine correct_rho_with_drho_2d

  subroutine correct_rho_with_drho_3d(rho,drho,ng_1,ng_2,glo,ghi,tlo,thi)

    integer          :: glo(3), ghi(3), ng_1, ng_2,tlo(:),thi(:)
    real(kind=dp_t)  ::  rho(glo(1)-ng_1:,glo(2)-ng_1:,glo(3)-ng_1:,:) ! density- last dim for #species
    real(kind=dp_t)  :: drho(glo(1)-ng_2:,glo(2)-ng_2:,glo(3)-ng_2:,:) ! total density in each cell

    ! local variables
    integer          :: i,j,k

    ! for specific box, now start loops over alloted cells
    do k=tlo(3), thi(3)
       do j=tlo(2), thi(2)
          do i=tlo(1), thi(1)

             call correct_rho_with_drho_local(rho(i,j,k,:),drho(i,j,k,:))

          end do
       end do
    end do

  end subroutine correct_rho_with_drho_3d

  subroutine correct_rho_with_drho_local(rho,drho)

    real(kind=dp_t), intent(inout) :: rho(nspecies)    ! density- last dim for #species
    real(kind=dp_t), intent(out)   :: drho(nspecies)   ! total density in each cell

    ! local variables
    real(kind=dp_t)  :: rhotot_local
    real(kind=dp_t)  :: rho_old(nspecies)

    ! make a copy of the input rho's
    rho_old = rho

    ! total rho in the cell
    rhotot_local = sum(rho)

    ! make sure each rho is greater than +epsilon
    ! to prevent division by zero in case species
    ! density, molar concentration or total density = 0.
    rho(1:nspecies) = max(rho(1:nspecies), fraction_tolerance*rhotot_local)

    ! keep track of how much rho changed so we can reset it back later
    drho = rho - rho_old

  end subroutine correct_rho_with_drho_local

  subroutine convert_cons_to_prim(mla,rho,rhotot,molarconc,molmtot)

   type(ml_layout), intent(in   )  :: mla
   type(multifab) , intent(in   )  :: rho(:)
   type(multifab) , intent(inout)  :: rhotot(:)
   type(multifab) , intent(inout)  :: molarconc(:)
   type(multifab) , intent(inout)  :: molmtot(:)

   ! local variables
   integer :: lo(rho(1)%dim), hi(rho(1)%dim)
   integer :: n,i,dm,nlevs,ng_1,ng_2,ng_3,ng_4

   ! pointer for rho(nspecies), rhotot(1), molarconc(nspecies)
   real(kind=dp_t), pointer        :: dp1(:,:,:,:)  ! for rho
   real(kind=dp_t), pointer        :: dp2(:,:,:,:)  ! for rhotot
   real(kind=dp_t), pointer        :: dp3(:,:,:,:)  ! for molar concentrations
   real(kind=dp_t), pointer        :: dp4(:,:,:,:)  ! for total molar concentrations

   type(mfiter) :: mfi
   type(box) :: tilebox
   integer :: tlo(mla%dim), thi(mla%dim)

   type(bl_prof_timer), save :: bpt

   call build(bpt,"convert_cons_to_prim")

   dm    = mla%dim     ! dimensionality
   ng_1  = rho(1)%ng   ! number of ghost cells
   ng_2  = rhotot(1)%ng
   ng_3  = molarconc(1)%ng
   ng_4  = molmtot(1)%ng
   nlevs = mla%nlevel  ! number of levels

   if (ng_3 .ne. ng_4) then
      call bl_error('convert_cons_to_prim: ng for molarconc and molmass differ')
   end if

   !$omp parallel private(n,i,mfi,tilebox,tlo,thi,dp1,dp2,dp3,dp4,lo,hi)

    ! loop over all boxes
    do n=1,nlevs
       call mfiter_build(mfi, rho(n), tiling=.true.)

       do while (more_tile(mfi))
          i = get_fab_index(mfi)

          tilebox = get_growntilebox(mfi,molarconc(n)%ng)
          tlo = lwb(tilebox)
          thi = upb(tilebox)

 !      do i=1,nfabs(rho(n))
          dp1 => dataptr(rho(n),i)
          dp2 => dataptr(rhotot(n),i)
          dp3 => dataptr(molarconc(n),i)
          dp4 => dataptr(molmtot(n),i)
          lo = lwb(get_box(rho(n),i))
          hi = upb(get_box(rho(n),i))
          select case(dm)
          case (2)
             call compute_molconc_rhotot_2d(dp1(:,:,1,:),dp2(:,:,1,1), &
                                            dp3(:,:,1,:),dp4(:,:,1,1), &
                                            ng_1,ng_2,ng_3,ng_4,lo,hi,tlo,thi)
          case (3)
             call compute_molconc_rhotot_3d(dp1(:,:,:,:),dp2(:,:,:,1), &
                                            dp3(:,:,:,:),dp4(:,:,:,1), &
                                            ng_1,ng_2,ng_3,ng_4,lo,hi,tlo,thi)
          end select
       end do
    end do
    !$omp end parallel

    call destroy(bpt)

  end subroutine convert_cons_to_prim

  subroutine compute_molconc_rhotot_2d(rho,rhotot,molarconc,molmtot, &
                                       ng_1,ng_2,ng_3,ng_4,glo,ghi,tlo,thi)

    integer          :: glo(2), ghi(2), ng_1, ng_2, ng_3, ng_4
    integer          :: tlo(2), thi(2)
    real(kind=dp_t)  ::       rho(glo(1)-ng_1:,glo(2)-ng_1:,:)       ! density- last dim for #species
    real(kind=dp_t)  ::    rhotot(glo(1)-ng_2:,glo(2)-ng_2:)     ! total density in each cell
    real(kind=dp_t)  :: molarconc(glo(1)-ng_3:,glo(2)-ng_3:,:) ! molar concentration
    real(kind=dp_t)  ::   molmtot(glo(1)-ng_4:,glo(2)-ng_4:)     ! total molar mass

    ! local variables
    integer          :: i,j

    ! for specific box, now start loops over alloted cells
    do j=tlo(2), thi(2)
       do i=tlo(1), thi(1)

         call compute_molconc_rhotot_local(rho(i,j,:),rhotot(i,j), &
                                           molarconc(i,j,:),molmtot(i,j))

       end do
    end do

  end subroutine compute_molconc_rhotot_2d

  subroutine compute_molconc_rhotot_3d(rho,rhotot,molarconc,molmtot, &
                                       ng_1,ng_2,ng_3,ng_4,glo,ghi,tlo,thi)

    integer          :: glo(3), ghi(3), ng_1, ng_2, ng_3, ng_4
    integer          :: tlo(3), thi(3)
    real(kind=dp_t)  ::       rho(glo(1)-ng_1:,glo(2)-ng_1:,glo(3)-ng_1:,:) ! density- last dim for #species
    real(kind=dp_t)  ::    rhotot(glo(1)-ng_2:,glo(2)-ng_2:,glo(3)-ng_2:)   ! total density in each cell
    real(kind=dp_t)  :: molarconc(glo(1)-ng_3:,glo(2)-ng_3:,glo(3)-ng_3:,:) ! molar concentration
    real(kind=dp_t)  ::   molmtot(glo(1)-ng_4:,glo(2)-ng_4:,glo(3)-ng_4:)   ! total molar mass

    ! local variables
    integer          :: i,j,k

    ! for specific box, now start loops over alloted cells
    do k=tlo(3), thi(3)
       do j=tlo(2), thi(2)
          do i=tlo(1), thi(1)

             call compute_molconc_rhotot_local(rho(i,j,k,:),rhotot(i,j,k),&
                                               molarconc(i,j,k,:),molmtot(i,j,k))

          end do
       end do
    end do

  end subroutine compute_molconc_rhotot_3d

  subroutine compute_molconc_rhotot_local(rho,rhotot,molarconc,molmtot)

    real(kind=dp_t), intent(in)   :: rho(nspecies)       ! density- last dim for #species
    real(kind=dp_t), intent(out)  :: rhotot             ! total density in each cell
    real(kind=dp_t), intent(out)  :: molarconc(nspecies) ! molar concentration
    real(kind=dp_t), intent(out)  :: molmtot             ! total molar mass

    ! local variables
    integer          :: n
    real(kind=dp_t), dimension(nspecies) :: W            ! mass fraction w_i = rho_i/rho
    real(kind=dp_t)  :: Sum_woverm, rhotot_local

    ! calculate total density inside each cell
    rhotot_local=0.d0
    do n=1, nspecies
       rhotot_local = rhotot_local + rho(n)
    end do
    rhotot = rhotot_local

    ! calculate mass fraction and total molar mass (1/m=Sum(w_i/m_i))
    Sum_woverm=0.d0
    do n=1, nspecies
       W(n) = rho(n)/rhotot
       Sum_woverm = Sum_woverm + W(n)/molmass(n)
    end do
    molmtot = 1.0d0/Sum_woverm

    ! calculate molar concentrations in each cell (x_i=m*w_i/m_i)
    do n=1, nspecies
       molarconc(n) = molmtot*W(n)/molmass(n)
    end do

  end subroutine compute_molconc_rhotot_local

  subroutine compute_rhotot(mla,rho,rhotot)

   type(ml_layout), intent(in   )  :: mla
   type(multifab) , intent(in   )  :: rho(:)
   type(multifab) , intent(inout)  :: rhotot(:)

   ! local variables
   integer :: lo(mla%dim), hi(mla%dim)
   integer :: n,i,dm,nlevs,ng_1,ng_2

   ! pointer for rho(nspecies), rhotot(1), molarconc(nspecies)
   real(kind=dp_t), pointer        :: dp1(:,:,:,:)   ! for rho
   real(kind=dp_t), pointer        :: dp2(:,:,:,:)  ! for rhotot

   type(mfiter) :: mfi
   type(box) :: tilebox
   integer :: tlo(mla%dim), thi(mla%dim)

   type(bl_prof_timer), save :: bpt

   call build(bpt, "compute_rhotot")

   dm    = mla%dim     ! dimensionality
   ng_1  = rho(1)%ng   ! number of ghost cells
   ng_2  = rhotot(1)%ng
   nlevs = mla%nlevel  ! number of levels

   !$omp parallel private(n,i,mfi,tilebox,tlo,thi,dp1,dp2,lo,hi)

    ! loop over all boxes
    do n=1,nlevs
       call mfiter_build(mfi, rho(n), tiling=.true.)

       do while (more_tile(mfi))
          i = get_fab_index(mfi)

          tilebox = get_tilebox(mfi)
          tlo = lwb(tilebox)
          thi = upb(tilebox)

!       do i=1,nfabs(rho(n))
          dp1 => dataptr(rho(n),i)
          dp2 => dataptr(rhotot(n),i)
          lo = lwb(get_box(rho(n),i))
          hi = upb(get_box(rho(n),i))

          select case(dm)
          case (2)
             call compute_rhotot_2d(dp1(:,:,1,:),dp2(:,:,1,1),&
                                    ng_1,ng_2,lo,hi,tlo,thi)
          case (3)
             call compute_rhotot_3d(dp1(:,:,:,:),dp2(:,:,:,1),&
                                    ng_1,ng_2,lo,hi,tlo,thi)
          end select
       end do
    end do
    !$omp end parallel

    call destroy(bpt)

  end subroutine compute_rhotot

  subroutine compute_rhotot_2d(rho,rhotot,ng_1,ng_2,glo,ghi,tlo,thi)

    integer          :: glo(2), ghi(2), ng_1, ng_2, tlo(2), thi(2)
    real(kind=dp_t)  ::    rho(glo(1)-ng_1:,glo(2)-ng_1:,:) ! density- last dim for #species
    real(kind=dp_t)  :: rhotot(glo(1)-ng_2:,glo(2)-ng_2:)   ! total density in each cell

    ! local variables
    integer          :: i,j,n

    ! for specific box, now start loops over alloted cells
    do j=tlo(2),thi(2)
       do i=tlo(1),thi(1)

          ! calculate total density inside each cell
          rhotot(i,j)=0.d0
          do n=1, nspecies
             rhotot(i,j) = rhotot(i,j) + rho(i,j,n)
          end do

       end do
    end do

  end subroutine compute_rhotot_2d

  subroutine compute_rhotot_3d(rho,rhotot,ng_1,ng_2,glo,ghi,tlo,thi)

    integer          :: glo(3), ghi(3), ng_1, ng_2, tlo(3), thi(3)
    real(kind=dp_t)  ::    rho(glo(1)-ng_1:,glo(2)-ng_1:,glo(3)-ng_1:,:) ! density- last dim for #species
    real(kind=dp_t)  :: rhotot(glo(1)-ng_2:,glo(2)-ng_2:,glo(3)-ng_2:)   ! total density in each cell

    ! local variables
    integer          :: i,j,k,n

    ! for specific box, now start loops over alloted cells
    do k=tlo(3),thi(3)
       do j=tlo(2),thi(2)
          do i=tlo(1),thi(1)

             ! calculate total density inside each cell
             rhotot(i,j,k)=0.d0
             do n=1,nspecies
                rhotot(i,j,k) = rhotot(i,j,k) + rho(i,j,k,n)
             end do

          end do
       end do
    end do

  end subroutine compute_rhotot_3d

  subroutine compute_Gama(mla,rho,rhotot,molarconc,molmtot,Hessian,Gama)

    type(ml_layout), intent(in   )  :: mla
    type(multifab),  intent(in   )  :: rho(:)
    type(multifab),  intent(in   )  :: rhotot(:)
    type(multifab),  intent(in   )  :: molarconc(:)
    type(multifab),  intent(in   )  :: molmtot(:)
    type(multifab),  intent(in   )  :: Hessian(:)    ! Hessian matrix
    type(multifab),  intent(inout)  :: Gama(:)       ! Non-ideality coeficient

    ! local variables
    integer :: lo(rho(1)%dim), hi(rho(1)%dim)
    integer :: n,i,dm,nlevs,ng_0,ng_1,ng_2,ng_3,ng_4,ng_5

    ! assign pointers for multifabs to be passed
    real(kind=dp_t), pointer        :: dp0(:,:,:,:)  ! for rho
    real(kind=dp_t), pointer        :: dp1(:,:,:,:)  ! for rhotot
    real(kind=dp_t), pointer        :: dp2(:,:,:,:)  ! for molarconc
    real(kind=dp_t), pointer        :: dp3(:,:,:,:)  ! for molmtot
    real(kind=dp_t), pointer        :: dp4(:,:,:,:)  ! for Hessian
    real(kind=dp_t), pointer        :: dp5(:,:,:,:)  ! for Gama

    type(mfiter) :: mfi
    type(box) :: tilebox
    integer :: tlo(mla%dim), thi(mla%dim)

    type(bl_prof_timer), save :: bpt

    call build(bpt, "compute_Gama")

    dm    = mla%dim     ! dimensionality
    ng_0  = rho(1)%ng   ! number of ghost cells
    ng_1  = rhotot(1)%ng
    ng_2  = molarconc(1)%ng
    ng_3  = molmtot(1)%ng
    ng_4  = Hessian(1)%ng
    ng_5  = Gama(1)%ng
    nlevs = mla%nlevel  ! number of levels

   !$omp parallel private(n,i,mfi,tilebox,tlo,thi) &
   !$omp private(dp0,dp1,dp2,dp3,dp4,dp5,lo,hi)

    ! loop over all boxes
    do n=1,nlevs
       call mfiter_build(mfi, Gama(n), tiling=.true.)

       do while (more_tile(mfi))
          i = get_fab_index(mfi)

          tilebox = get_growntilebox(mfi,rho(n)%ng)
          tlo = lwb(tilebox)
          thi = upb(tilebox)
!       do i=1,nfabs(rho(n))
          dp0 => dataptr(rho(n),i)
          dp1 => dataptr(rhotot(n),i)
          dp2 => dataptr(molarconc(n),i)
          dp3 => dataptr(molmtot(n),i)
          dp4 => dataptr(Hessian(n),i)
          dp5 => dataptr(Gama(n),i)
          lo = lwb(get_box(rho(n),i))
          hi = upb(get_box(rho(n),i))
          select case(dm)
          case (2)
             call compute_Gama_2d(dp0(:,:,1,:),dp1(:,:,1,1),dp2(:,:,1,:),&
                                  dp3(:,:,1,1),dp4(:,:,1,:),dp5(:,:,1,:), &
                                  ng_0,ng_1,ng_2,ng_3,ng_4,ng_5,lo,hi,tlo,thi)
          case (3)
             call compute_Gama_3d(dp0(:,:,:,:),dp1(:,:,:,1),dp2(:,:,:,:),&
                                  dp3(:,:,:,1),dp4(:,:,:,:),dp5(:,:,:,:), &
                                  ng_0,ng_1,ng_2,ng_3,ng_4,ng_5,lo,hi,tlo,thi)
          end select
       end do
    end do
    !$omp end parallel

    call destroy(bpt)

  end subroutine compute_Gama

  subroutine compute_Gama_2d(rho,rhotot,molarconc,molmtot,Hessian,Gama, &
                             ng_0,ng_1,ng_2,ng_3,ng_4,ng_5,glo,ghi,tlo,thi)

    integer          :: glo(2), ghi(2), ng_0,ng_1,ng_2,ng_3,ng_4,ng_5
    integer          :: tlo(2), thi(2)
    real(kind=dp_t)  ::       rho(glo(1)-ng_0:,glo(2)-ng_0:,:) ! density; last dimension for species
    real(kind=dp_t)  ::    rhotot(glo(1)-ng_1:,glo(2)-ng_1:)   ! total density in each cell
    real(kind=dp_t)  :: molarconc(glo(1)-ng_2:,glo(2)-ng_2:,:) ! molar concentration
    real(kind=dp_t)  ::   molmtot(glo(1)-ng_3:,glo(2)-ng_3:)   ! total molar mass
    real(kind=dp_t)  ::   Hessian(glo(1)-ng_4:,glo(2)-ng_4:,:) ! last dimension for nspecies^2
    real(kind=dp_t)  ::      Gama(glo(1)-ng_5:,glo(2)-ng_5:,:) ! last dimension for nspecies^2

    ! local varialbes
    integer          :: i,j

    ! for specific box, now start loops over alloted cells
    do j=tlo(2),thi(2)
       do i=tlo(1),thi(1)

          call compute_Gama_local(rho(i,j,:),rhotot(i,j),molarconc(i,j,:),&
                                  molmtot(i,j),Hessian(i,j,:),Gama(i,j,:))

       end do
    end do

  end subroutine compute_Gama_2d

  subroutine compute_Gama_3d(rho,rhotot,molarconc,molmtot,Hessian,Gama, &
                             ng_0,ng_1,ng_2,ng_3,ng_4,ng_5,glo,ghi,tlo,thi)

    integer          :: glo(3), ghi(3), ng_0,ng_1,ng_2,ng_3,ng_4,ng_5
    integer          :: tlo(3), thi(3)
    real(kind=dp_t)  ::       rho(glo(1)-ng_0:,glo(2)-ng_0:,glo(3)-ng_0:,:) ! density; last dimension for species
    real(kind=dp_t)  ::    rhotot(glo(1)-ng_1:,glo(2)-ng_1:,glo(3)-ng_1:)   ! total density in each cell
    real(kind=dp_t)  :: molarconc(glo(1)-ng_2:,glo(2)-ng_2:,glo(3)-ng_2:,:) ! molar concentration;
    real(kind=dp_t)  ::   molmtot(glo(1)-ng_3:,glo(2)-ng_3:,glo(3)-ng_3:)   ! total molar mass
    real(kind=dp_t)  ::   Hessian(glo(1)-ng_4:,glo(2)-ng_4:,glo(3)-ng_4:,:) ! last dimension for nspecies^2
    real(kind=dp_t)  ::      Gama(glo(1)-ng_5:,glo(2)-ng_5:,glo(3)-ng_5:,:) ! last dimension for nspecies^2

    ! local varialbes
    integer          :: i,j,k

    ! for specific box, now start loops over alloted cells
    do k=tlo(3),thi(3)
       do j=tlo(2),thi(2)
          do i=tlo(1),thi(1)

             call compute_Gama_local(rho(i,j,k,:),rhotot(i,j,k),molarconc(i,j,k,:),&
                                     molmtot(i,j,k),Hessian(i,j,k,:),Gama(i,j,k,:))
          end do
       end do
    end do

  end subroutine compute_Gama_3d

  subroutine compute_Gama_local(rho,rhotot,molarconc,molmtot,Hessian,Gama)

    real(kind=dp_t), intent(in)   :: rho(nspecies)
    real(kind=dp_t), intent(in)   :: rhotot
    real(kind=dp_t), intent(in)   :: molarconc(nspecies)
    real(kind=dp_t), intent(in)   :: molmtot
    real(kind=dp_t), intent(in)   :: Hessian(nspecies,nspecies)
    real(kind=dp_t), intent(out)  :: Gama(nspecies,nspecies)

    ! local variables
    integer                                       :: row,column
    real(kind=dp_t), dimension(nspecies,nspecies) :: I, X_xxT

    ! free the memory
    I=0; X_xxT=0;

    ! populate I and X_xxT, where X = molmtot*W*M^(-1)
    do row=1, nspecies
       do column=1, row-1
          I(row, column)    = 0.d0
          I(column, row)    = I(row, column) ! symmetric

          if(.not. is_ideal_mixture) then
             X_xxT(row,column)   = -molarconc(row)*molarconc(column)  ! form x*transpose(x) off diagonals
             X_xxT(column, row)  = X_xxT(row, column)                 ! symmetric
          end if
       end do

       I(row, row) = 1.d0        ! unit matrix for ideal mixture
       if(.not. is_ideal_mixture) then
          X_xxT(row,row) = molarconc(row) - molarconc(row)**2
       end if
    end do

    ! compute Gama
    Gama = I + matmul(X_xxT, Hessian)

  end subroutine compute_Gama_local

  subroutine compute_chi(mla,rho,rhotot,molarconc,chi,D_bar,D_therm,Temp,zeta_by_Temp)

    type(ml_layout), intent(in   )  :: mla
    type(multifab) , intent(in   )  :: rho(:)
    type(multifab) , intent(in   )  :: rhotot(:)
    type(multifab) , intent(in   )  :: molarconc(:)
    type(multifab) , intent(inout)  :: chi(:)
    type(multifab) , intent(in   )  :: D_bar(:)
    type(multifab) , intent(in   )  :: D_therm(:)
    type(multifab) , intent(in   )  :: Temp(:)
    type(multifab) , intent(inout)  :: zeta_by_Temp(:)

    ! local variables
    integer :: lo(rho(1)%dim), hi(rho(1)%dim)
    integer :: n,i,dm,nlevs,ng_0,ng_1,ng_2,ng_3,ng_4,ng_5,ng_6,ng_7

    ! pointer for rho(nspecies), molarconc(nspecies)
    real(kind=dp_t), pointer        :: dp0(:,:,:,:)  ! for rho
    real(kind=dp_t), pointer        :: dp1(:,:,:,:)  ! for rhotot
    real(kind=dp_t), pointer        :: dp2(:,:,:,:)  ! for molarconc
    real(kind=dp_t), pointer        :: dp3(:,:,:,:)  ! for chi
    real(kind=dp_t), pointer        :: dp4(:,:,:,:)  ! for D_bar
    real(kind=dp_t), pointer        :: dp5(:,:,:,:)  ! for Temp
    real(kind=dp_t), pointer        :: dp6(:,:,:,:)  ! for zeta_by_Temp
    real(kind=dp_t), pointer        :: dp7(:,:,:,:)  ! for D_therm

    type(mfiter) :: mfi
    type(box) :: tilebox
    integer :: tlo(mla%dim), thi(mla%dim)

    type(bl_prof_timer), save :: bpt

    call build(bpt,"compute_chi")

    dm = mla%dim        ! dimensionality
    ng_0 = rho(1)%ng    ! number of ghost cells
    ng_1 = rhotot(1)%ng
    ng_2 = molarconc(1)%ng
    ng_3 = chi(1)%ng
    ng_4 = D_bar(1)%ng
    ng_5 = Temp(1)%ng
    ng_6 = zeta_by_Temp(1)%ng
    ng_7 = D_therm(1)%ng
    nlevs = mla%nlevel  ! number of levels

    !$omp parallel private(n,i,mfi,tilebox,tlo,thi) &
    !$omp private(dp0,dp1,dp2,dp3,dp4,dp5,dp6,dp7,lo,hi)

    ! loop over all boxes
    do n=1,nlevs
       call mfiter_build(mfi, rho(n), tiling=.true.)

       do while (more_tile(mfi))
          i = get_fab_index(mfi)

          tilebox = get_growntilebox(mfi,chi(n)%ng)
          tlo = lwb(tilebox)
          thi = upb(tilebox)

!       do i=1,nfabs(rho(n))
          dp0  => dataptr(rho(n), i)
          dp1 => dataptr(rhotot(n), i)
          dp2 => dataptr(molarconc(n), i)
          dp3 => dataptr(chi(n), i)
          dp4 => dataptr(D_bar(n), i)
          dp5 => dataptr(Temp(n), i)
          dp6 => dataptr(zeta_by_Temp(n), i)
          dp7 => dataptr(D_therm(n), i)
          lo  =  lwb(get_box(rho(n), i))
          hi  =  upb(get_box(rho(n), i))
          select case(dm)
          case (2)
             call compute_chi_2d(dp0(:,:,1,:),dp1(:,:,1,1),dp2(:,:,1,:),dp3(:,:,1,:),&
                                 dp4(:,:,1,:),dp5(:,:,1,1),dp6(:,:,1,:),dp7(:,:,1,:),&
                                 ng_0,ng_1,ng_2,ng_3,ng_4,ng_5,ng_6,ng_7,lo,hi,tlo,thi)
          case (3)
             call compute_chi_3d(dp0(:,:,:,:),dp1(:,:,:,1),dp2(:,:,:,:),dp3(:,:,:,:),&
                                 dp4(:,:,:,:),dp5(:,:,:,1),dp6(:,:,:,:),dp7(:,:,:,:),&
                                 ng_0,ng_1,ng_2,ng_3,ng_4,ng_5,ng_6,ng_7,lo,hi,tlo,thi)
          end select
       end do
    end do
    !$omp end parallel

    call destroy(bpt)

  end subroutine compute_chi

  subroutine compute_chi_2d(rho,rhotot,molarconc,chi,D_bar,Temp,zeta_by_Temp,D_therm, &
                            ng_0,ng_1,ng_2,ng_3,ng_4,ng_5,ng_6,ng_7,glo,ghi,tlo,thi)

    integer          :: glo(2), ghi(2), ng_0,ng_1,ng_2,ng_3,ng_4,ng_5,ng_6,ng_7
    integer          :: tlo(2), thi(2)
    real(kind=dp_t)  ::          rho(glo(1)-ng_0:,glo(2)-ng_0:,:) ! density; last dimension for species
    real(kind=dp_t)  ::       rhotot(glo(1)-ng_1:,glo(2)-ng_1:)   ! total density in each cell
    real(kind=dp_t)  ::    molarconc(glo(1)-ng_2:,glo(2)-ng_2:,:) ! molar concentration
    real(kind=dp_t)  ::          chi(glo(1)-ng_3:,glo(2)-ng_3:,:) ! last dimension for nspecies^2
    real(kind=dp_t)  ::        D_bar(glo(1)-ng_4:,glo(2)-ng_4:,:) ! MS diff-coefs
    real(kind=dp_t)  ::         Temp(glo(1)-ng_5:,glo(2)-ng_5:)   ! Temperature
    real(kind=dp_t)  :: zeta_by_Temp(glo(1)-ng_6:,glo(2)-ng_6:,:) ! zeta/T
    real(kind=dp_t)  ::      D_therm(glo(1)-ng_7:,glo(2)-ng_7:,:) ! thermo diff-coefs

    ! local variables
    integer          :: i,j,row
    real(kind=dp_t), dimension(nspecies,nspecies) :: chilocal

    ! for specific box, now start loops over alloted cells
    do j=tlo(2),thi(2)
       do i=tlo(1),thi(1)

          call compute_chi_local(rho(i,j,:),rhotot(i,j),molarconc(i,j,:),&
                                 chi(i,j,:),D_bar(i,j,:),Temp(i,j),zeta_by_Temp(i,j,:),&
                                 D_therm(i,j,:))

          ! print chi for one cell
          if(.false.) then
          if(i.eq.32 .and. j.eq.16) then
            if(use_lapack) then
              print*, 'print chi via inverse/p-inverse'
            else
              print*, 'print chi via iterative methods'
            end if
            call set_Xij(chilocal, chi(i,j,:))
            do row=1, nspecies
               print*, chilocal(row, :)
            end do
          end if
          end if

       end do
    end do

  end subroutine compute_chi_2d

  subroutine compute_chi_3d(rho,rhotot,molarconc,chi,D_bar,Temp,zeta_by_Temp,D_therm, &
                            ng_0,ng_1,ng_2,ng_3,ng_4,ng_5,ng_6,ng_7,glo,ghi,tlo,thi)

    integer          :: glo(3), ghi(3), ng_0,ng_1,ng_2,ng_3,ng_4,ng_5,ng_6,ng_7
    integer          :: tlo(3), thi(3)
    real(kind=dp_t)  ::          rho(glo(1)-ng_0:,glo(2)-ng_0:,glo(3)-ng_0:,:) ! density; last dimension for species
    real(kind=dp_t)  ::       rhotot(glo(1)-ng_1:,glo(2)-ng_1:,glo(3)-ng_1:)   ! total density in each cell
    real(kind=dp_t)  ::    molarconc(glo(1)-ng_2:,glo(2)-ng_2:,glo(3)-ng_2:,:) ! molar concentration;
    real(kind=dp_t)  ::          chi(glo(1)-ng_3:,glo(2)-ng_3:,glo(3)-ng_3:,:) ! last dimension for nspecies^2
    real(kind=dp_t)  ::        D_bar(glo(1)-ng_4:,glo(2)-ng_4:,glo(3)-ng_4:,:) ! SM diffusion constants
    real(kind=dp_t)  ::         Temp(glo(1)-ng_5:,glo(2)-ng_5:,glo(3)-ng_5:)   ! Temperature
    real(kind=dp_t)  :: zeta_by_Temp(glo(1)-ng_6:,glo(2)-ng_6:,glo(3)-ng_6:,:) ! zeta/T
    real(kind=dp_t)  ::      D_therm(glo(1)-ng_7:,glo(2)-ng_7:,glo(3)-ng_7:,:) ! thermo diffusion constants

    ! local variables
    integer          :: i,j,k

    ! for specific box, now start loops over alloted cells
    do k=tlo(3),thi(3)
       do j=tlo(2),thi(2)
          do i=tlo(1),thi(1)

             call compute_chi_local(rho(i,j,k,:),rhotot(i,j,k),molarconc(i,j,k,:),&
                                    chi(i,j,k,:),D_bar(i,j,k,:),Temp(i,j,k),zeta_by_Temp(i,j,k,:),&
                                    D_therm(i,j,k,:))

          end do
       end do
    end do

  end subroutine compute_chi_3d

  subroutine compute_chi_local(rho,rhotot,molarconc,chi,D_bar,Temp,zeta_by_Temp,D_therm)

    real(kind=dp_t), intent(inout)  :: rho(nspecies)
    real(kind=dp_t), intent(in)     :: rhotot
    real(kind=dp_t), intent(inout)  :: molarconc(nspecies)
    real(kind=dp_t), intent(out)    :: chi(nspecies,nspecies)
    real(kind=dp_t), intent(in)     :: D_bar(nspecies,nspecies)
    real(kind=dp_t), intent(in)     :: Temp
    real(kind=dp_t), intent(inout)  :: zeta_by_Temp(nspecies)
    real(kind=dp_t), intent(in)     :: D_therm(nspecies)

    ! local variables
    integer                         :: row,column
    real(kind=dp_t)                 :: Sum_knoti

    ! vectors and matrices to be used by LAPACK
    real(kind=dp_t), dimension(nspecies,nspecies) :: Lambda
    real(kind=dp_t), dimension(nspecies)          :: W

    ! free up memory
    Lambda        = 0.d0
    W             = 0.d0

    ! compute Lambda_ij matrix and massfraction W_i = rho_i/rho; molarconc is
    ! expressed in terms of molmtot,mi,rhotot etc.
    do row=1, nspecies
       do column=1, row-1
          Lambda(row, column) = -molarconc(row)*molarconc(column)/D_bar(row,column)
          Lambda(column, row) = Lambda(row, column)
       end do
       W(row) = rho(row)/rhotot
    end do

    ! compute Lambda_ii
    do row=1, nspecies
       Sum_knoti = 0.d0
       do column=1, nspecies
          if(column.ne.row) then
             Sum_knoti = Sum_knoti - Lambda(row,column)
          end if
          Lambda(row,row) = Sum_knoti
       end do
    end do

    ! compute zeta_by_Temp for thermodiffusion
    if(is_nonisothermal) then
       do row=1, nspecies
          Sum_knoti = 0.d0
          do column=1, nspecies
             if(column.ne.row) then
                Sum_knoti = Sum_knoti + Lambda(row,column)*(D_therm(row)-D_therm(column))
             end if
             zeta_by_Temp(row) = Sum_knoti/Temp
          end do
       end do
    end if

    ! compute chi either selecting inverse/pseudoinverse or iterative methods
    if(use_lapack) then
       call compute_chi_lapack(Lambda(:,:),chi(:,:),W(:))
    else
       call Dbar2chi_iterative(chi_iterations,D_bar(:,:),molarconc(:),chi(:,:))
    end if

  end subroutine compute_chi_local

  subroutine compute_chi_lapack(Lambda,chi,W)

    real(kind=dp_t)  :: Lambda(nspecies,nspecies)
    real(kind=dp_t)  :: chi(nspecies,nspecies)
    real(kind=dp_t)  :: W(nspecies)

    ! local variables
    integer          :: row,column,info
    real(kind=dp_t)  :: alpha

    ! vectors and matrices to be used by LAPACK
    real(kind=dp_t), dimension(nspecies,nspecies) :: Sdag,chilocal
    real(kind=dp_t), dimension(nspecies,nspecies) :: U, UT, V, VT
    real(kind=dp_t), dimension(nspecies)          :: S, work
    integer,         dimension(nspecies)          :: ipiv

    ! free up the memory
    Sdag     = 0.d0
    U        = 0.d0
    UT       = 0.d0
    V        = 0.d0
    VT       = 0.d0
    S        = 0.d0
    work     = 0.d0
    alpha    = 0.d0
    chilocal = 0.d0

    ! calculate trace(Lambda)
    alpha = 0.d0
    do row=1, nspecies
       alpha = alpha + Lambda(row,row)
    end do

    ! calculate Lambda + alpha*W*WT (Equation 6)
    do row=1, nspecies
       do column=1, nspecies
          chilocal(row,column) = alpha*W(row)*W(column) + Lambda(row,column)
       end do
    end do

    !===============================================================
    ! select LAPACK inversion type, 1=inverse, 2=pseudo inverse
    !===============================================================
    select case(inverse_type)

    case(1)
    !==========================================================
    ! Using Inverse
    !==========================================================

    ! compute chilocal inverse
    call dgetrf(nspecies, nspecies, chilocal, nspecies, ipiv, info)
    call dgetri(nspecies, chilocal, nspecies, ipiv, work, nspecies, info)

    ! populate chi with B^(-1)
    chi = chilocal

    case(2)
    !==========================================================
    ! Using pseudoinverse
    !==========================================================

    ! SVD decomposition of chilocal = U * S * VTranspose; note that chilocal
    ! is changed. also V=(VT)T, UT = (U)T are needed for pseudoinverse of chilocal.
    call la_gesvd(chilocal, S, U, VT)
    V = transpose(VT)
    UT = transpose(U)

    ! populate diagonal matrix Sdag = 1/S with diagonal=0 below a chosen tolerance
    do row=1, nspecies
       do column=1,nspecies
          Sdag(row,column) = 0.0d0
       end do

       if(S(row).gt.fraction_tolerance*sum(S)) then
          Sdag(row,row) = 1.0d0/S(row)
       else
          Sdag(row,row) = 0.0d0
       end if
    end do

    ! compute chi = V*Sdag*UT, the pseudoinverse of chilocal
    chi = matmul(V, matmul(Sdag, UT))

    end select
    !===============================================================

    ! compute chi as equation (6)
    do row=1, nspecies
       do column=1, nspecies
          chi(row,column) = chi(row,column) - 1.0d0/alpha
       end do
    end do

  end subroutine compute_chi_lapack

  subroutine compute_Lonsager(mla,rho,rhotot,molarconc,molmtot,chi,Gama,Lonsager)

    type(ml_layout), intent(in   )  :: mla
    type(multifab) , intent(in   )  :: rho(:)
    type(multifab) , intent(in   )  :: rhotot(:)
    type(multifab) , intent(in   )  :: molarconc(:)
    type(multifab) , intent(in   )  :: molmtot(:)
    type(multifab) , intent(in   )  :: chi(:)
    type(multifab) , intent(in   )  :: Gama(:)
    type(multifab) , intent(inout)  :: Lonsager(:)

    ! local variables
    integer :: lo(rho(1)%dim), hi(rho(1)%dim)
    integer :: n,i,dm,nlevs,ng_0,ng_1,ng_2,ng_3,ng_4,ng_5,ng_6

    ! pointer for rho(nspecies), molarconc(nspecies)
    real(kind=dp_t), pointer        :: dp0(:,:,:,:)  ! for rho
    real(kind=dp_t), pointer        :: dp1(:,:,:,:)  ! for rhotot
    real(kind=dp_t), pointer        :: dp2(:,:,:,:)  ! for molarconc
    real(kind=dp_t), pointer        :: dp3(:,:,:,:)  ! for molmtot
    real(kind=dp_t), pointer        :: dp4(:,:,:,:)  ! for chi
    real(kind=dp_t), pointer        :: dp5(:,:,:,:)  ! for Gama
    real(kind=dp_t), pointer        :: dp6(:,:,:,:)  ! for Lonsager

    type(bl_prof_timer), save :: bpt

    call build(bpt, "compute_Lonsager")

    dm = mla%dim        ! dimensionality
    ng_0 = rho(1)%ng    ! number of ghost cells
    ng_1 = rhotot(1)%ng
    ng_2 = molarconc(1)%ng
    ng_3 = molmtot(1)%ng
    ng_4 = chi(1)%ng
    ng_5 = Gama(1)%ng
    ng_6 = Lonsager(1)%ng
    nlevs = mla%nlevel  ! number of levels

    ! loop over all boxes
    do n=1,nlevs
       do i=1,nfabs(rho(n))
          dp0 => dataptr(rho(n), i)
          dp1 => dataptr(rhotot(n), i)
          dp2 => dataptr(molarconc(n),i)
          dp3 => dataptr(molmtot(n),i)
          dp4 => dataptr(chi(n), i)
          dp5 => dataptr(Gama(n), i)
          dp6 => dataptr(Lonsager(n), i)
          lo  =  lwb(get_box(rho(n), i))
          hi  =  upb(get_box(rho(n), i))

          select case(dm)
          case (2)
             call compute_Lonsager_2d(dp0(:,:,1,:),dp1(:,:,1,1),dp2(:,:,1,:),&
                                      dp3(:,:,1,1),dp4(:,:,1,:),dp5(:,:,1,:),dp6(:,:,1,:),&
                                      ng_0,ng_1,ng_2,ng_3,ng_4,ng_5,ng_6,lo,hi)
          case (3)
             call compute_Lonsager_3d(dp0(:,:,:,:),dp1(:,:,:,1),dp2(:,:,:,:),&
                                      dp3(:,:,:,1),dp4(:,:,:,:),dp5(:,:,:,:),dp6(:,:,:,:),&
                                      ng_0,ng_1,ng_2,ng_3,ng_4,ng_5,ng_6,lo,hi)
          end select
       end do
    end do

    call destroy(bpt)

  end subroutine compute_Lonsager

  subroutine compute_Lonsager_2d(rho,rhotot,molarconc,molmtot,chi,Gama,Lonsager, &
                                 ng_0,ng_1,ng_2,ng_3,ng_4,ng_5,ng_6,lo,hi)

    integer          :: lo(2), hi(2), ng_0,ng_1,ng_2,ng_3,ng_4,ng_5,ng_6
    real(kind=dp_t)  ::       rho(lo(1)-ng_0:,lo(2)-ng_0:,:) ! density; last dimension for species
    real(kind=dp_t)  ::    rhotot(lo(1)-ng_1:,lo(2)-ng_1:)   ! total density in each cell
    real(kind=dp_t)  :: molarconc(lo(1)-ng_2:,lo(2)-ng_2:,:) ! molar concentration
    real(kind=dp_t)  ::   molmtot(lo(1)-ng_3:,lo(2)-ng_3:)   ! total molar mass
    real(kind=dp_t)  ::       chi(lo(1)-ng_4:,lo(2)-ng_4:,:) ! last dimension for nspecies^2
    real(kind=dp_t)  ::      Gama(lo(1)-ng_5:,lo(2)-ng_5:,:) ! non-ideality coeficient
    real(kind=dp_t)  ::  Lonsager(lo(1)-ng_6:,lo(2)-ng_6:,:) ! last dimension for nspecies^2

    ! local variables
    integer          :: i,j,row
    real(kind=dp_t), dimension(nspecies,nspecies) :: Lonsager_local

    ! for specific box, now start loops over alloted cells
    do j=lo(2)-ng_6,hi(2)+ng_6
       do i=lo(1)-ng_6,hi(1)+ng_6

          call compute_Lonsager_local(rho(i,j,:),rhotot(i,j),molarconc(i,j,:),&
                                      molmtot(i,j),chi(i,j,:),Gama(i,j,:),Lonsager(i,j,:))

          if(.false.) then
          if(i.eq.7 .and. j.eq.14) then
          call set_Xij(Lonsager_local, Lonsager(i,j,:))
            write(*,*), "Lonsager"
            do row=1, nspecies
               write(*,*), Lonsager_local(row,:)
               !write(*,*), sum(Lonsager_local(row,:))
            end do
            do row=1, nspecies
               !write(*,*), sum(Lonsager_local(:,row))
            end do
          end if
          end if

       end do
    end do

  end subroutine compute_Lonsager_2d

  subroutine compute_Lonsager_3d(rho,rhotot,molarconc,molmtot,chi,Gama,Lonsager, &
                                 ng_0,ng_1,ng_2,ng_3,ng_4,ng_5,ng_6,lo,hi)

    integer          :: lo(3), hi(3), ng_0,ng_1,ng_2,ng_3,ng_4,ng_5,ng_6
    real(kind=dp_t)  ::       rho(lo(1)-ng_0:,lo(2)-ng_0:,lo(3)-ng_0:,:) ! density; last dimension for species
    real(kind=dp_t)  ::    rhotot(lo(1)-ng_1:,lo(2)-ng_1:,lo(3)-ng_1:)   ! total density in each cell
    real(kind=dp_t)  :: molarconc(lo(1)-ng_2:,lo(2)-ng_2:,lo(3)-ng_2:,:) ! molar concentration
    real(kind=dp_t)  ::   molmtot(lo(1)-ng_3:,lo(2)-ng_3:,lo(3)-ng_3:)   ! total molar mass
    real(kind=dp_t)  ::       chi(lo(1)-ng_4:,lo(2)-ng_4:,lo(3)-ng_4:,:) ! last dimension for nspecies^2
    real(kind=dp_t)  ::      Gama(lo(1)-ng_5:,lo(2)-ng_5:,lo(3)-ng_5:,:) ! non-ideality coeficient
    real(kind=dp_t)  ::  Lonsager(lo(1)-ng_6:,lo(2)-ng_6:,lo(3)-ng_6:,:) ! last dimension for nspecies^2

    ! local variables
    integer          :: i,j,k

    ! for specific box, now start loops over alloted cells
    do k=lo(3)-ng_6,hi(3)+ng_6
       do j=lo(2)-ng_6,hi(2)+ng_6
          do i=lo(1)-ng_6,hi(1)+ng_6

             call compute_Lonsager_local(rho(i,j,k,:),rhotot(i,j,k),molarconc(i,j,k,:),&
                                         molmtot(i,j,k),chi(i,j,k,:),Gama(i,j,k,:),Lonsager(i,j,k,:))

         end do
      end do
    end do

  end subroutine compute_Lonsager_3d

subroutine compute_Lonsager_local(rho,rhotot,molarconc,molmtot,chi,Gama,Lonsager)

    real(kind=dp_t), intent(in)   :: rho(nspecies)
    real(kind=dp_t), intent(in)   :: rhotot
    real(kind=dp_t), intent(in)   :: molarconc(nspecies)      ! molar concentration
    real(kind=dp_t), intent(in)   :: molmtot                  ! total molar mass
    real(kind=dp_t), intent(in)   :: chi(nspecies,nspecies)   ! rank conversion done
    real(kind=dp_t), intent(in)   :: Gama(nspecies,nspecies)
    real(kind=dp_t), intent(out)  :: Lonsager(nspecies,nspecies)

    ! local variables
    integer                              :: row,column,info
    real(kind=dp_t), dimension(nspecies) :: W
    real(kind=dp_t)                      :: rcond

    type(bl_prof_timer), save :: bpt

    call build(bpt,"compute_Lonsager_local")

    ! compute massfraction W_i = rho_i/rho;
    do row=1, nspecies
       W(row) = rho(row)/rhotot
    end do

    ! compute Onsager matrix L
    do column=1, nspecies
       do row=1, nspecies
          Lonsager(row, column) = molmtot*rhotot*W(row)*chi(row,column)*W(column)/k_B
       end do
    end do

    ! compute cell-centered Cholesky factor of Lonsager
    if(use_lapack) then

       call dpotrf_f95(Lonsager,'L', rcond, 'I', info)

       ! remove all upper-triangular entries and NXN entry that lapack doesn't set to zero
       do row=1, nspecies
          do column=row+1, nspecies
             Lonsager(row, column) = 0.0d0
          end do
       end do
       Lonsager(nspecies, nspecies) = 0.0d0

    else
       call choldc(Lonsager,nspecies)
    end if

    call destroy(bpt)

  end subroutine compute_Lonsager_local

  subroutine compute_rhoWchi(mla,rho,rhotot,chi,rhoWchi)

    type(ml_layout), intent(in   )  :: mla
    type(multifab) , intent(in   )  :: rho(:)
    type(multifab) , intent(in   )  :: rhotot(:)
    type(multifab) , intent(in   )  :: chi(:)
    type(multifab) , intent(inout)  :: rhoWchi(:)

    ! local variables
    integer :: lo(mla%dim), hi(mla%dim)
    integer :: n,i,dm,nlevs,ng_1,ng_2,ng_3,ng_4

    ! pointer for rho(nspecies), molarconc(nspecies)
    real(kind=dp_t), pointer        :: dp1(:,:,:,:)  ! for rho
    real(kind=dp_t), pointer        :: dp2(:,:,:,:)  ! for rhotot
    real(kind=dp_t), pointer        :: dp3(:,:,:,:)  ! for chi
    real(kind=dp_t), pointer        :: dp4(:,:,:,:)  ! for rhoWchi

    type(mfiter) :: mfi
    type(box) :: tilebox
    integer :: tlo(mla%dim), thi(mla%dim)

    type(bl_prof_timer), save :: bpt

    call build(bpt,"compute_rhoWchi")

    dm = mla%dim        ! dimensionality
    ng_1 = rho(1)%ng    ! number of ghost cells
    ng_2 = rhotot(1)%ng
    ng_3 = chi(1)%ng
    ng_4 = rhoWchi(1)%ng
    nlevs = mla%nlevel  ! number of levels

    !$omp parallel private(n,i,mfi,tilebox,tlo,thi) &
    !$omp private(dp1,dp2,dp3,dp4,lo,hi)

    ! loop over all boxes
    do n=1,nlevs
       call mfiter_build(mfi, rho(n), tiling=.true.)

       do while (more_tile(mfi))
          i = get_fab_index(mfi)

          tilebox = get_growntilebox(mfi,rhoWchi(n)%ng)
          tlo = lwb(tilebox)
          thi = upb(tilebox)

!       do i=1,nfabs(rho(n))
          dp1 => dataptr(rho(n), i)
          dp2 => dataptr(rhotot(n), i)
          dp3 => dataptr(chi(n), i)
          dp4 => dataptr(rhoWchi(n), i)
          lo  =  lwb(get_box(rho(n), i))
          hi  =  upb(get_box(rho(n), i))

          select case(dm)
          case (2)
             call compute_rhoWchi_2d(dp1(:,:,1,:),dp2(:,:,1,1),dp3(:,:,1,:),dp4(:,:,1,:), &
                                     ng_1,ng_2,ng_3,ng_4,lo,hi,tlo,thi)
          case (3)
             call compute_rhoWchi_3d(dp1(:,:,:,:),dp2(:,:,:,1),dp3(:,:,:,:),dp4(:,:,:,:), &
                                     ng_1,ng_2,ng_3,ng_4,lo,hi,tlo,thi)
          end select
       end do
    end do
    !$omp end parallel

    call destroy(bpt)

  end subroutine compute_rhoWchi

  subroutine compute_rhoWchi_2d(rho,rhotot,chi,rhoWchi,ng_1,ng_2,ng_3,ng_4,glo,ghi,tlo,thi)

    integer          :: glo(2), ghi(2), ng_1,ng_2,ng_3,ng_4,tlo(2),thi(2)
    real(kind=dp_t)  ::     rho(glo(1)-ng_1:,glo(2)-ng_1:,:) ! density; last dimension for species
    real(kind=dp_t)  ::  rhotot(glo(1)-ng_2:,glo(2)-ng_2:)   ! total density in each cell
    real(kind=dp_t)  ::     chi(glo(1)-ng_3:,glo(2)-ng_3:,:) ! last dimension for nspecies^2
    real(kind=dp_t)  :: rhoWchi(glo(1)-ng_4:,glo(2)-ng_4:,:) ! last dimension for nspecies^2

    ! local variables
    integer          :: i,j,row,column
    real(kind=dp_t), dimension(nspecies,nspecies) :: rhoWchiloc

    ! for specific box, now start loops over alloted cells
    do j=tlo(2),thi(2)
       do i=tlo(1),thi(1)

          call compute_rhoWchi_local(rho(i,j,:),rhotot(i,j),chi(i,j,:),rhoWchi(i,j,:))

          if(.false.) then
          if(i.eq.7 .and. j.eq.14) then
          call set_Xij(rhoWchiloc, rhoWchi(i,j,:))
             do row=1, nspecies
                do column=1, nspecies
                   print*, rhoWchiloc(row,column)
                end do
                print*, ''
             end do
          end if
          end if

       end do
    end do

  end subroutine compute_rhoWchi_2d

  subroutine compute_rhoWchi_3d(rho,rhotot,chi,rhoWchi,ng_1,ng_2,ng_3,ng_4,glo,ghi,tlo,thi)

    integer          :: glo(3), ghi(3), ng_1,ng_2,ng_3,ng_4,tlo(3),thi(3)
    real(kind=dp_t)  ::     rho(glo(1)-ng_1:,glo(2)-ng_1:,glo(3)-ng_1:,:) ! density; last dimension for species
    real(kind=dp_t)  ::  rhotot(glo(1)-ng_2:,glo(2)-ng_2:,glo(3)-ng_2:)   ! total density in each cell
    real(kind=dp_t)  ::     chi(glo(1)-ng_3:,glo(2)-ng_3:,glo(3)-ng_3:,:) ! last dimension for nspecies^2
    real(kind=dp_t)  :: rhoWchi(glo(1)-ng_4:,glo(2)-ng_4:,glo(3)-ng_4:,:) ! last dimension for nspecies^2

    ! local variables
    integer          :: i,j,k

    ! for specific box, now start loops over alloted cells
    do k=tlo(3),thi(3)
       do j=tlo(2),thi(2)
          do i=tlo(1),thi(1)

             call compute_rhoWchi_local(rho(i,j,k,:),rhotot(i,j,k),chi(i,j,k,:),rhoWchi(i,j,k,:))

         end do
      end do
    end do

  end subroutine compute_rhoWchi_3d

  subroutine compute_rhoWchi_local(rho,rhotot,chi,rhoWchi)

    real(kind=dp_t), intent(in)   :: rho(nspecies)
    real(kind=dp_t), intent(in)   :: rhotot
    real(kind=dp_t), intent(in)   :: chi(nspecies,nspecies)   ! rank conversion done
    real(kind=dp_t), intent(out)  :: rhoWchi(nspecies,nspecies)

    ! local variables
    integer                              :: row,column
    real(kind=dp_t), dimension(nspecies) :: W !,chiw

    ! compute massfraction W_i = rho_i/rho;
    do row=1, nspecies
       W(row) = rho(row)/rhotot
    end do

    ! populate -rho*W*chi = -rho_i*chi
    do row=1, nspecies
       do column=1, nspecies
          rhoWchi(row,column) = -rho(row)*chi(row,column)
       end do
    end do

  end subroutine compute_rhoWchi_local

  subroutine set_Xij(Xout_ij, Xin_ij)

    real(kind=dp_t), dimension(nspecies,nspecies), intent(in)  :: Xin_ij
    real(kind=dp_t), dimension(nspecies,nspecies), intent(out) :: Xout_ij

    type(bl_prof_timer), save :: bpt

    call build(bpt, "set_Xij")

    ! reshape array into matrix without doing index algebra
    Xout_ij = Xin_ij

    call destroy(bpt)

  end subroutine set_Xij

  subroutine compute_baro_coef(mla,baro_coef,rho,rhotot,Temp)

    type(ml_layout), intent(in   )  :: mla
    type(multifab) , intent(inout)  :: baro_coef(:)
    type(multifab) , intent(in   )  :: rho(:)
    type(multifab) , intent(in   )  :: rhotot(:)
    type(multifab) , intent(in   )  :: Temp(:)

    ! local
    integer :: lo(mla%dim), hi(mla%dim)
    integer :: n,i,dm,nlevs,ng_1,ng_2,ng_3,ng_4

    real(kind=dp_t), pointer        :: dp1(:,:,:,:)  ! for baro_coef
    real(kind=dp_t), pointer        :: dp2(:,:,:,:)  ! for rho
    real(kind=dp_t), pointer        :: dp3(:,:,:,:)  ! for rhotot
    real(kind=dp_t), pointer        :: dp4(:,:,:,:)  ! for Temp

    type(bl_prof_timer), save :: bpt

    call build(bpt,"compute_baro_coef")

    dm = mla%dim
    nlevs = mla%nlevel

    ng_1 = baro_coef(1)%ng
    ng_2 = rho(1)%ng
    ng_3 = rhotot(1)%ng
    ng_4 = Temp(1)%ng

    do n=1,nlevs
       do i=1,nfabs(rho(n))
          dp1 => dataptr(baro_coef(n),i)
          dp2 => dataptr(rho(n),i)
          dp3 => dataptr(rhotot(n),i)
          dp4 => dataptr(Temp(n),i)
          lo = lwb(get_box(rho(n),i))
          hi = upb(get_box(rho(n),i))
          select case(dm)
          case (2)
             call compute_baro_coef_2d(dp1(:,:,1,:),dp2(:,:,1,:),dp3(:,:,1,1),dp4(:,:,1,1), &
                                        ng_1,ng_2,ng_3,ng_4,lo,hi)
          case (3)
             call bl_error("compute_baro_coef_3d not written yet")
          end select
       end do
    end do

    call destroy(bpt)

  end subroutine compute_baro_coef

  subroutine compute_baro_coef_2d(baro_coef,rho,rhotot,Temp,ng_1,ng_2,ng_3,ng_4,lo,hi)

    integer        , intent(in   ) :: lo(2), hi(2), ng_1, ng_2, ng_3, ng_4
    real(kind=dp_t), intent(inout) :: baro_coef(lo(1)-ng_1:,lo(2)-ng_1:,:)
    real(kind=dp_t), intent(in   ) ::       rho(lo(1)-ng_2:,lo(2)-ng_2:,:)
    real(kind=dp_t), intent(in   ) ::    rhotot(lo(1)-ng_3:,lo(2)-ng_3:)
    real(kind=dp_t), intent(in   ) ::      Temp(lo(1)-ng_4:,lo(2)-ng_4:)

    ! local variables
    integer :: i,j,comp
    real(kind=dp_t) :: n

    do j=lo(2)-ng_1,hi(2)+ng_1
       do i=lo(1)-ng_1,hi(1)+ng_1

          n = 0.d0
          do comp=1,nspecies
             n = n + rho(i,j,comp)/molmass(comp)
          end do

          do comp=1,nspecies
             baro_coef(i,j,comp) = rho(i,j,comp)/rhobar(comp) / (n*k_B*Temp(i,j))
          end do

       end do
    end do

  end subroutine compute_baro_coef_2d

end module mass_flux_utilities_module

max_stars_repo_name: BoxLib-Codes/LowMachMixtures

max_issues_repo_name: BoxLib-Codes/LowMachMixtures

max_forks_repo_name: BoxLib-Codes/LowMachMixtures


text:
MODULE definitions_module

  IMPLICIT NONE
  TYPE field_type
    REAL(KIND=8),    DIMENSION(:,:), ALLOCATABLE :: density0,density1
  END TYPE field_type

  TYPE tile_type
    TYPE(field_type):: field
  END TYPE tile_type

  TYPE chunk_type
    TYPE(tile_type), DIMENSION(:), ALLOCATABLE :: tiles
  END TYPE chunk_type
  TYPE(chunk_type)       :: chunk
END MODULE definitions_module
SUBROUTINE start
  use definitions_module
  ALLOCATE(chunk%tiles (1))
  ALLOCATE(chunk%tiles(1)%field%density0  (-1:962,-1:962))
  ALLOCATE(chunk%tiles(1)%field%density1  (-1:962,-1:962))

 print '(Z)', loc(chunk%tiles(1))
 print '(Z)', loc(chunk%tiles(1)%field)
 print '(Z)', loc(chunk%tiles(1)%field%density0)
 print '(Z)', loc(chunk%tiles(1)%field%density1)

 call flush(6)
!$omp target data map(                   &
!$omp   chunk%tiles(1)%field%density0,   &
!$omp   chunk%tiles(1)%field%density1)

  DO tile=1,2
    chunk%tiles(1)%field%density0(tile,-1) =0
    chunk%tiles(1)%field%density1(tile,-1) =0
  ENDDO
!$omp end target data

END SUBROUTINE start
program foobar
   call start
end program

max_stars_repo_name: raramakr/aomp

max_issues_repo_name: raramakr/aomp

max_forks_repo_name: raramakr/aomp


text:
! Program to test the transpose intrinsic
program intrinsic_transpose
   integer, dimension (3, 3) :: a, b
   complex(kind=8), dimension (2, 2) :: c, d
   complex(kind=4), dimension (2, 2) :: e

   a = 0
   b = reshape ((/1, 2, 3, 4, 5, 6, 7, 8, 9/), (/3, 3/))
   a = transpose (b)
   if (any (a .ne. reshape ((/1, 4, 7, 2, 5, 8, 3, 6, 9/), (/3, 3/)))) &
      call abort
   c = (0.0, 0.0)
   d = reshape ((/(1d0,2d0), (3d0, 4d0), (5d0, 6d0), (7d0, 8d0)/), (/2, 2/))
   c = transpose (d);
   if (any (c .ne. reshape ((/(1d0, 2d0), (5d0, 6d0), &
                              (3d0, 4d0), (7d0, 8d0)/), (/2, 2/)))) &
    call abort ();

   e = reshape ((/(1.0,2.0), (3.0, 4.0), (5.0, 6.0), (7.0, 8.0)/), (/2, 2/))
   e = transpose (e);
   if (any (e .ne. reshape ((/(1.0, 2.0), (5.0, 6.0), &
                              (3.0, 4.0), (7.0, 8.0)/), (/2, 2/)))) &
    call abort ();
end program

max_issues_repo_name: vidkidz/crossbridge

max_forks_repo_name: vidkidz/crossbridge


text:
module helloworld
    implicit none

    contains

    subroutine hello
        write(*,*) "Hello world!"
    end subroutine hello
end module helloworld

max_stars_repo_name: jaimedelacruz/ISPy

max_issues_repo_name: jaimedelacruz/ISPy

max_forks_repo_name: jaimedelacruz/ISPy


module analysis_mod

  use real_type_mod
  use variables_mod
  use output_mod
  use parameters_mod
  use compute_fft_mod
  use master_function_mod

  implicit none

contains

  subroutine analysis_create_all(out_dir)

    character(len=200), intent(in) :: out_dir

    real(kind=wp), dimension(:), allocatable :: freq
    integer :: istat, i
    real(kind=wp) :: pi

    pi = 4.0_wp * atan(1.0_wp)

    allocate(freq(np/2 + 1), &
             STAT = istat)
    if (istat .ne. 0) then
      write(*,*) 'Allocation error in analysis_create_all'
      stop
    end if
    do i = 1, np/2 + 1
      freq(i) = 2.0_wp * pi * dble(i-1) / (xmax - xmin)
    end do

    call create_output_file(out_dir, 'n_up',    N_UP_FILE)
    call create_output_file(out_dir, 'n_down',  N_DOWN_FILE)
    call create_output_file(out_dir, 'mu_up',   MU_UP_FILE)
    call create_output_file(out_dir, 'mu_down', MU_DOWN_FILE)
    call create_output_file(out_dir, 'n_sqd',   N_SQD_FILE)
    call create_output_file(out_dir, 'v_up',    V_UP_FILE)
!!$    call create_output_file(out_dir, 'c_cc',    C_CC_FILE)
!!$    call create_output_file(out_dir, 'c_ent',   C_ENT_FILE)
!!$    call create_output_file(out_dir, 'cs2',     CS2_FILE)
    call create_output_file(out_dir, 'constraints', CONSTRAINTS_FILE)

    if (do_fft) then
      if (vtk_output) then
        call create_output_timefreq_vtk_file(out_dir, 'fft', FFT_FILE, &
             np, freq, fft_out_total, fft_out_dt)
      else
        call create_output_file(out_dir, 'fft',    FFT_FILE)
      end if
    end if

    deallocate(freq)

  end subroutine analysis_create_all

  subroutine analysis_close_all()

    close(N_UP_FILE   )
    close(N_DOWN_FILE )
    close(MU_UP_FILE  )
    close(MU_DOWN_FILE)
    close(N_SQD_FILE  )
    close(V_UP_FILE   )
!!$    close(C_CC_FILE   )
!!$    close(C_ENT_FILE  )
!!$    close(CS2_FILE    )
    close(CONSTRAINTS_FILE)

    if (do_fft) then
      close(FFT_FILE)
    end if

  end subroutine analysis_close_all

  subroutine analysis_output_all(x, n_up, n_down, mu_up, mu_down, mu_down_p, &
       n_sqd, t, dt, dx)

    real(kind=wp), dimension(1-gz:np+gz), intent(in) :: x
    real(kind=wp), dimension(NF, NV, 1-gz:np+gz), intent(in) :: &
         n_up, n_down, mu_up, mu_down, mu_down_p
    real(kind=wp), dimension(NFF, 1-gz:np+gz), intent(in) :: &
         n_sqd
    real(kind=wp), intent(in) :: t, dt, dx

    real(kind=wp), dimension(:, :), allocatable :: v_up, cs2, constraints
    real(kind=wp), dimension(:), allocatable :: c_cc, c_ent
    real(kind=wp), dimension(NFF) :: dLdn
    real(kind=wp) :: d2Ldn2
    real(kind=wp), dimension(NF) :: d2Ldn2_2

    integer :: istat, f, f1, f2, i

    allocate(v_up(NFF, 1-gz:np+gz), &
             cs2 (2, 1-gz:np+gz), &
             c_cc(     1-gz:np+gz), &
             c_ent(    1-gz:np+gz), &
             constraints(NF, 1-gz:np+gz), &
             STAT=istat)

    call output2(N_UP_FILE   , x, n_up   , NF, NV, np, gz, t)
    call output2(N_DOWN_FILE , x, n_down , NF, NV, np, gz, t)
    call output2(MU_UP_FILE  , x, mu_up  , NF, NV, np, gz, t)
    call output2(MU_DOWN_FILE, x, mu_down, NF, NV, np, gz, t)
    call output (N_SQD_FILE  , x, n_sqd  , NFF   , np, gz, t)

    do f = 1, NF
      v_up(f, :) = n_up(f, 2, :) / n_up(f, 1, :)
    end do
    f = NF
    do f1 = 1, NF
      do f2 = f1 + 1, NF
        f = f + 1
        v_up(f, :) = sqrt(1.0_wp - &
             n_sqd(f1, :) * n_sqd(f2, :) / n_sqd(f, :)**2)
      end do
    end do

    call output (V_UP_FILE   , x, v_up   , NFF   , np, gz, t)

!!$    These quantities are from Lars et al's two stream instability paper
!!$    but in the end did not turn out to be particularly useful
!!$    do i = 1-gz, np+gz
!!$      call MasterFunction_Derivs(n_sqd, dLdn)
!!$      call MasterFunction_SecondDerivs(n_sqd(:, i), d2Ldn2, d2Ldn2_2)
!!$      c_cc(i) = 1.0_wp / ( dLdn(1) * dLdn(2) ) * &
!!$           (2.0_wp * sqrt(n_sqd(1, i) * n_sqd(2, i)) * d2Ldn2)**2
!!$      c_ent(i) = dLdn(3)**2 / ( dLdn(1) * dLdn(2) )
!!$      do f1 = 1, NF
!!$        cs2(f1, i) = sqrt(n_sqd(f1, i)) / dLdn(f1) * d2Ldn2_2(f1) + 1
!!$      end do
!!$    end do
!!$
!!$    call output (C_CC_FILE   , x, c_cc   , 1     , np, gz, t)
!!$    call output (C_ENT_FILE  , x, c_ent  , 1     , np, gz, t)
!!$    call output (CS2_FILE    , x, cs2    , 2     , np, gz, t)

!!$    Compute the constraint
!!$    d_t mu^X_x = d_x mu^X_t
!!$    using first order differencing in time and
!!$    central second order differencing in space

!!$    This can only be done after the first timestep
    if (t .ge. 0.5_wp * dt) then
      constraints = 0.0_wp
      do f = 1, NF
        do i = 1, np
          constraints(f, i) = (mu_down(f, 2, i) - mu_down_p(f, 2, i)) / dt
          constraints(f, i) = constraints(f, i) - 0.5_wp * &
               (mu_down(f, 1, i+1) - mu_down(f, 1, i-1)) / dx
        end do
      end do

      call output (CONSTRAINTS_FILE    , x, constraints, NF, np, gz, t)
    end if

    deallocate(v_up, c_cc, c_ent, cs2, constraints)

  end subroutine analysis_output_all

  subroutine analysis_fft(n_sqd, n_up, t)

    real(kind=wp), dimension(NF, NV, 1-gz:np+gz), intent(in) :: n_up
    real(kind=wp), dimension(NFF, 1-gz:np+gz), intent(in) :: n_sqd
    real(kind=wp), intent(in) :: t

    real(kind=wp), dimension(:), allocatable :: freq
    complex(kind=wp), dimension(:, :), allocatable :: fft
    integer :: istat, i, ff, vv
    real(kind=wp) :: pi

    pi = 4.0_wp * atan(1.0_wp)

    allocate(freq(np/2 + 1), &
             fft(NF*NV, np/2 + 1), &
             STAT = istat)
    if (istat .ne. 0) then
      write(*,*) 'Allocation error in analysis_fft'
      stop
    end if

    do ff = 1, NF
      do vv = 1, NV
        call compute_fft(np, n_up(ff, vv, 1:np), fft(NF*(ff-1)+vv, :))
      end do
    end do

    if (vtk_output) then
      call output_timefreq_vtk(FFT_FILE, fft, np)
    else
      do i = 1, np/2 + 1
        freq(i) = 2.0_wp * pi * dble(i-1) / (xmax - xmin)
      end do
      call output3(FFT_FILE, freq, fft, np, NV, NF, t)
    end if

    deallocate(freq, fft)

  end subroutine analysis_fft

end module analysis_mod


PROGRAM TrapezoidalRuleForIntegration
    IMPLICIT NONE
    REAL:: a, b, Area
    INTEGER::  N      ! N is number of segment
    WRITE(*, fmt = '(/A)', ADVANCE = 'NO') "Enter the number of segment: "
    READ(*, *) N
    WRITE(*, fmt = '(/A)') "Enter the value of a & b"
    READ(*, *) a, b

    CALL TrapezoidalRule(a, b, Area, N)
    WRITE(*, *) "Area under the curve is: ", Area

    CONTAINS
    SUBROUTINE TrapezoidalRule(a, b, Area, N)
        IMPLICIT NONE
        REAL:: a, b, Area, Xi
        INTEGER:: N, i
        ! N is number of segment
        IF(N .LT. 1) THEN
            Area = 0.0
            RETURN
        ENDIF
        Area = (f(a) + f(b)) / 2.0
        DO i = 1, N-1
            Xi = a + i * (b-a)/N
            Area = Area + f(Xi)
        ENDDO
        Area = ((b-a)/N) * Area
    END SUBROUTINE TrapezoidalRule

    REAL FUNCTION f(x)
        IMPLICIT NONE
        REAL:: x
        f = 1.0/x
        return
    END FUNCTION f

END PROGRAM TrapezoidalRuleForIntegration


text:
C    @(#)ldptisec.f	20.14 5/27/99
C****************************************************************
C
C     File: ldptisec.f
C
C     Purpose: Routine to load PTI Multi-section Line Grouping data
C              from raw data file
C
c     Return code:  n = 0 : Success
c                   n = 1 " Error
c
C     Author: Walt Powell  Date: 21 May 1996
C     Called by: load_pti.f
C
C****************************************************************
      integer function ldptisec (xbuf, options, numver, error)
      integer error, options(*), numver
      character *(*) xbuf

      include 'ipfinc/parametr.inc'

      include 'ipfinc/prt.inc'
      include 'ipfinc/pti_data.inc'
      include 'ipfinc/blank.inc'
      include 'ipfinc/bus.inc'
      include 'ipfinc/branch.inc'

      character word(26)*10, id*1, busname*8, temp(2)*120, newtemp*120,
     &          bus1*8, bus2*8, bus3*8
      logical first, found, finished
      integer gtptinam, find_bus, status, nsec, section(4,13),
     &        ptr, pold, rename_bus, count, bptr1, bptr2, ptr9, ptr10,
     &        gtpieqiv

      ldptisec = 0
      error = 0
      last = index (xbuf, '/') - 1
      if (last .le. 0) last = lastch(xbuf)
      call ptiscan (xbuf(1:last), word, nwrd, ',',  ' ')

      do i = nwrd+1, 24
        word(i) = ' '
      enddo
c
c     Set defaults
c
      if (word(3) .eq. ' ') word(3) = '&1'
      if (word(4) .eq. ' ') word(4) = '1'

      status = gtptinam (word(1), npti1, kp1, isgn1, bus1, basekv1)
      if (npti1 .eq. 0 .and. bus1 .eq. ' ') then
        ldptisec = 1
        go to 900
      endif
      if (bus1 .eq. ' ') then
        write (errbuf(1), 100) npti1
  100   format (' Terminal 1 Bus (', i6,
     &    ') for Multi-section line grouping record (', i6,
     &    ') is not in raw data file')
        errbuf(2) = ' (' // xbuf(1:60) // ')'
        call prterx ('W',2)
        error = 1
        go to 900
      endif
      k1 = find_bus (bus1, basekv1)
      if (k1 .le. 0) then
        write (errbuf(1), 110) npti1, bus1, basekv1
  110   format (' Terminal 1 bus (', i6,
     &    1x, a8, f7.1,
     &   ') for Multi-section line grouping record is not in system')
        errbuf(2) = ' (' // xbuf(1:60) // ')'
        call prterx ('W',2)
        error = 1
        go to 900
      endif
      status = gtptinam (word(2), npti2, kp2, isgn2, bus2, basekv2)
      if (bus2 .eq. ' ') then
        write (errbuf(1), 120) npti2
  120   format (' Terminal 2 Bus (', i6,
     &    ') for Multi-section line grouping record (', i6,
     &    ') is not in raw data file')
        errbuf(2) = ' (' // xbuf(1:60) // ')'
        call prterx ('W',2)
        error = 1
        go to 900
      endif
      k2 = find_bus (bus2, basekv2)
      if (k2 .le. 0) then
        write (errbuf(1), 130) npti2, bus2, basekv2
  130   format (' Terminal 2 bus (', i6, 1x, a8, f7.1,
     &   ') for Multi-section line grouping record is not in system')
        errbuf(2) = ' (' // xbuf(1:60) // ')'
        call prterx ('W',2)
        error = 1
        go to 900
      endif
      id = word(4)
      nsec = 1
c
c     First section is terminal bus 1
c
      section(1,nsec) = kp1
      section(2,nsec) = ichar(id)
      section(3,nsec) = k1
      section(4,nsec) = 0
      i = 1
      finished = .false.
      do while (i .le. 11 .and. .not. finished)
        j = 2*i + 3
        if (word(j) .eq. ' ') then
          finished = .true.
        else
          status = gtptinam (word(j), npti3, kp3, isgn3, bus3,
     &                       basekv3)
          section(1,i+1) = kp3
          section(2,i+1) = ichar(word(j+1)(1:1))
          section(3,i+1) = 0
          section(4,i+1) = 0
          if (kp3 .le. 0) then
            write (errbuf(1), 140) section(1,nsec)
  140       format (' Passive Bus (', i6,
     &   ') for Multi-section line grouping record is not in raw data fi
     &le')
            errbuf(2) = ' (' // xbuf(1:60) // ')'
            call prterx ('W',2)
            error = 1
            go to 900
          endif
          nsec = i + 1
          k3 = find_bus (bus3, basekv3)
          if (k3 .le. 0) then
            write (errbuf(1), 150) section(1,j), bus3, basekv3
  150       format (' Passive bus (', i6, 1x, a8, f7.1,
     &   ') for Multi-section line grouping record is not in system')
            errbuf(2) = ' (' // xbuf(1:60) // ')'
            call prterx ('W',2)
            error = 1
            go to 900
          endif
          section(3,nsec) = k3
          count = 0
          ptr = kbsdta(16,k3)
          do while (ptr .gt. 0)
            count = count + 1
            ptr = brnch_nxt(ptr)
          enddo
          if (count .ne. 2) then
            write (errbuf(1), 152) k3, bus(k3), base(k3)
  152       format (' Passive bus (', i6, 1x, a8, f7.1,
     &   ') for Multi-section line grouping record does not have exacty
     &two branches')
            errbuf(2) = ' (' // xbuf(1:60) // ')'
            call prterx ('W',2)
            error = 1
            go to 900
          endif
          i = i + 1
        endif
      enddo
c
c     Last section is terminal bus 2
c
      nsec = nsec + 1
      section(1,nsec) = kp2
      section(2,nsec) = ichar(id)
      section(3,nsec) = k2
      section(4,nsec) = 0
c
c     Add an entry for the pi-equivalent section
C
      if (ltot+1 .ge. MAXBRN) then
        write (errbuf(1), 170) MAXBRN
  170   format ('More than ',i5,
     &          ' branch records. Overflow occurred at branch:')
        errbuf(2) = ' (' // xbuf(1:60) // ')'
        call prterx ('E',2)
        ltot = 1
        error = 1
        go to 900
      endif

      ltot = ltot + 1
      call lkbrdata (ltot, k1, k2, bptr1, bptr2, error)
      kx(bptr1) = k1
      ky(bptr1) = k2
      brid(bptr1) = id
      brsect(bptr1) = 0
      brtype(bptr1) = 1

      kx(bptr2) = k2
      ky(bptr2) = k1
      brid(bptr2) = id
      brsect(bptr2) = 0
      brtype(bptr2) = 1

      do i = 1, 18
        brnch(i,ltot) = 0.0
      enddo
      call putchr (1, id, kbrnch(13,ltot))
      do i = 1, 3
        rateln(i,ltot) = 0.0
      enddo
      if (isgn1 .lt. 0) then
        kbrnch(15,ltot) = 2
      else
        kbrnch(15,ltot) = 1
      endif
c
c     Now delete all passive section buses and regroup their
c     associated branch data
c
      do i = 1, nsec-1
        k1x = section(3,i)
        k2x = section(3,i+1)
c
c       Delete branch pointers k1x - k2x
c
        pold = 0
        ptr = kbsdta(16,k1x)
        found = .false.
        do while (ptr .gt. 0 .and. .not. found)
          if (ky(ptr) .eq. k2x) then
            found = .true.
            section(4,i) = brnch_ptr(ptr)
            if (pold .eq. 0) then
              kbsdta(16,k1x) = brnch_nxt(ptr)
            else
              brnch_nxt(pold) = brnch_nxt(ptr)
            endif
          else
            pold = ptr
            ptr = brnch_nxt(ptr)
          endif
        enddo
        if (.not. found) then
          write (errbuf(1), 180) bus(k1x), base(k1x), bus(k2x),
     &      base(k2x)
  180     format (' Error regrouping multi-section branch (', a8,
     &      f7.1, 1x, a8, f7.1, ')')
          errbuf(2) = ' (' // xbuf(1:60) // ')'
          call prterx ('W',2)
          error = 1
          go to 900
        endif
c
c       Delete branch pointers k2x - k1x
c
        pold = 0
        ptr = kbsdta(16,k2x)
        found = .false.
        do while (ptr .gt. 0 .and. .not. found)
          if (ky(ptr) .eq. k1x) then
            found = .true.
            if (pold .eq. 0) then
              kbsdta(16,k2x) = brnch_nxt(ptr)
            else
              brnch_nxt(pold) = brnch_nxt(ptr)
            endif
          else
            pold = ptr
            ptr = brnch_nxt(ptr)
          endif
        enddo
        if (.not. found) then
          write (errbuf(1), 180) bus(k2x), base(k2x), bus(k1x),
     &      base(k1x)
          errbuf(2) = ' (' // xbuf(1:60) // ')'
          call prterx ('W',2)
          error = 1
          go to 900
        endif
c
c       Now relink branch (k1 - k2 - id - sect)
c
        nbr = section(4,i)
        inbr = iabs(nbr)
        call lkbrdata (nbr, k1, k2, bptr1, bptr2, error)
        kx(bptr1) = k1
        ky(bptr1) = k2
        brid(bptr1) = id
        brsect(bptr1) = i
        brtype(bptr1) = kbrnch(1,inbr)

        kx(bptr2) = k2
        ky(bptr2) = k1
        brid(bptr2) = id
        brsect(bptr2) = i
        brtype(bptr2) = kbrnch(1,inbr)
c
c       Transpose data if not oriented k1-k2
c
        if (brnch_ptr(bptr1) .lt. 0) call txposebr (bptr1, bptr2)

      enddo

      do i = 2, nsec-1
        k1x = section(3,i)

c       Remove bus from bus hash table by renaming
c
        busname = srtlst
        basekv = 9999.0
        status = rename_bus (k1x, busname, basekv)
c
c       Flag k1x as deleted
c
        do j = 1, 11
          kbsdta(j,k1x) = 0
        enddo
        capcor(1,k1x) = 0.0
        capcor(2,k1x) = -9.0e10
      enddo
c
c     If more that 9 sections, consolidate the 9-th and 10-th sections
c
      if (nsec .gt. 10) then
        write (errbuf(1), 190) nsec-1
  190   format ('More than ',i5, ' branch sections. 9-th and 10-th secti
     &ons consolidated')
        errbuf(2) = ' (' // xbuf(1:60) // ')'
        call prterx ('W', 2)

        ptr9 = kbsdta(16,k1)
        found = .false.
        do while (ptr9 .gt. 0 .and. .not. found)
          if (ky(ptr9) .eq. k2 .and. brid(ptr9) .eq. id .and.
     &        brsect(ptr9) .eq. 9) then
            found = .true.
            call bcdbrn (ptr9, temp(1), ierror)
          else
            ptr9 = brnch_nxt(ptr9)
          endif
        enddo

        ptr10 = kbsdta(16,k1)
        found = .false.
        do while (ptr10 .gt. 0 .and. .not. found)
          if (ky(ptr10) .eq. k2 .and. brid(ptr10) .eq. id .and.
     &        brsect(ptr10) .eq. 10) then
            found = .true.
            brsect(ptr10) = 9
            call bcdbrn (ptr10, temp(2), ierror)
            brsect(ptr10) = 10
          else
            ptr10 = brnch_nxt(ptr10)
          endif
        enddo

        if (ptr9 .eq. 0 .or. ptr10 .eq. 0) then
          write (errbuf(1), 200) bus(k1), base(k1), bus(k2), base(k2),
     &      id
  200     format ('Error consolidating sections 9 and 10 in branch ',
     &      a8, f6.1, 1x, a8, f7.1, 1x, a)
          errbuf(2) = ' (' // xbuf(1:60) // ')'
          call prterx ('E',2)
          error = 1
          go to 900
        endif
        numtemp = 2
        status = gtpieqiv (numtemp, temp, newtemp)
        if (status .eq. 0) then
          write (outbuf, 210)
  210     format ('Old sections 9 and 10:')
          call prtout (1)
          write (outbuf, 220) temp(1)(1:80)
  220     format (' (', a, ')')
          call prtout (1)
          write (outbuf, 220) temp(2)(1:80)
          call prtout (1)
          write (outbuf, 230)
  230     format ('New equivalent section 9:')
          call prtout (1)
          write (outbuf, 220) newtemp(1:80)
          call prtout (1)
          nbr = iabs(brnch_ptr(ptr9))
          read (newtemp, 240, err=900) (brnch(k,nbr),k=5,10)
  240     format (bz, 38x, 6f6.5)
c
c         Now delete section 10, proceding from k1 to k2
c
          pold = 0
          ptr = kbsdta(16,k1)
          found = .false.
          do while (ptr .gt. 0 .and. .not. found)
            if (ky(ptr) .eq. k2 .and. brid(ptr) .eq. id .and.
     &          brsect(ptr) .eq. 10) then
              found = .true.
              if (pold .eq. 0) then
                kbsdta(16,k1) = brnch_nxt(ptr)
              else
                brnch_ptr(pold) = brnch_nxt(ptr)
              endif
            else
              pold = ptr
              brnch_nxt(ptr) = ptr
            endif
          enddo
c
c         Now delete section 10, proceding from k2 to k1
c
          pold = 0
          ptr = kbsdta(16,k2)
          found = .false.
          do while (ptr .gt. 0 .and. .not. found)
            if (ky(ptr) .eq. k1 .and. brid(ptr) .eq. id .and.
     &          brsect(ptr) .eq. 10) then
              found = .true.
              if (pold .eq. 0) then
                kbsdta(16,k2) = brnch_nxt(ptr)
              else
                brnch_nxt(pold) = brnch_nxt(ptr)
              endif
            else
              pold = ptr
              brnch_nxt(ptr) = ptr
            endif
          enddo
c
c         Now delete section 10 in its entirety
c
          nbr = iabs(brnch_ptr(ptr10))
          do i = 1, 18
            kbrnch(i,nbr) = 0
          enddo
c
c         Flag ptr as deleted
c
          brnch_ptr(ptr10) = 0
          brnch_nxt(ptr10) = 0

        endif

      endif

  900 continue
      return
      end

ext: f

max_stars_repo_name: mbheinen/bpa-ipf-tsp

max_issues_repo_name: cuihantao/bpa-ipf-tsp

max_forks_repo_name: mbheinen/bpa-ipf-tsp


module heap_map_int_int_mod
    ! This module include
    ! double_heap_sort(key_arr, value_arr) O(NlogN)
    ! heap map
    ! 最小値を素早く返します。
    ! 最大値を素早く返したい場合
    !  1. 値の符号を変えてpush
    !  2. popやtopで取得した値は符号を元に戻すのを忘れずに。
    use,intrinsic :: iso_fortran_env
    implicit none
    private
    public:: double_heap_sort
    type, public:: heap_map_int_int
        integer(int32),private:: len
        integer(int64), allocatable,private:: key(:), val(:)
    contains
        procedure:: size => hm_size
        procedure:: remain => hm_remain
        procedure:: top => hm_top
        procedure:: push => hm_push
        procedure:: pop => hm_pop
        procedure:: to_array => hm_to_array
    end type
    interface heap_map_int_int
        module procedure init_hm, init_hm_with_len, init_hm_with_ar
    end interface
contains
    pure elemental function to_int64(num) result(ret)
        class(*), intent(in) :: num
        integer(int64):: ret

        select type (num)
            type is (integer(int64))
                ret = num
            type is (integer(int32))
                ret = int(num, int64)
            class default
                ret = 0_int64
                ret = 1_int64/ret
        end select
    end function

    function init_hm() result(hm)
        type(heap_map_int_int):: hm
        hm%len=0
        allocate(hm%key(1))
        allocate(hm%val(1))
    end function

    function init_hm_with_len(n) result(hm)
        integer(int32),intent(in):: n
        type(heap_map_int_int):: hm
        hm%len=0
        allocate(hm%key(n))
        allocate(hm%val(n))
    end function

    function init_hm_with_ar(k_ar, v_ar) result(hm)
        class(*),intent(in):: k_ar(:), v_ar(:)
        integer(int32):: i, l
        type(heap_map_int_int):: hm

        l = min(size(k_ar), size(v_ar))
        hm = heap_map_int_int(l)
        do i=1,l
            call hm%push(to_int64(k_ar(i)), to_int64(v_ar(i)))
        end do
    end function

    function hm_size(hm) result(ret)
        class(heap_map_int_int):: hm
        integer(int32):: ret

        ret = hm%len
    end function

    function hm_remain(hm) result(ret)
        class(heap_map_int_int):: hm
        logical:: ret

        ret = hm%size() > 0
    end function

    subroutine hm_push(hm,k,v)
        class(heap_map_int_int):: hm
        class(*):: k,v

        if (.not. allocated(hm%key)) then
            hm%len=0
            allocate(hm%key(1))
            allocate(hm%val(1))
        end if

        if (hm%len+1 >= size(hm%key)) call add(hm)
        hm%len=hm%len+1
        hm%key(hm%len) = to_int64(k)
        hm%val(hm%len) = to_int64(v)
        call heap_up(hm,hm%len)
    end subroutine

    subroutine hm_top(hm,k,v)
        class(heap_map_int_int):: hm
        integer(int64),intent(out):: k, v

        k=hm%key(1)
        v=hm%val(1)
    end subroutine

    subroutine hm_pop(hm,k,v)
        class(heap_map_int_int):: hm
        integer(int64),intent(out):: k, v

        k=hm%key(1)
        v=hm%val(1)
        hm%key(1) = hm%key(hm%len)
        hm%val(1) = hm%val(hm%len)
        hm%len=hm%len-1
        call heap_down(hm,1_int32)
    end subroutine

    subroutine add(hm)
        class(heap_map_int_int):: hm

        call add_array(hm%key)
        call add_array(hm%val)
    end subroutine

    subroutine add_array(ar)
        integer(int64),allocatable,intent(inout):: ar(:)
        integer(int64),allocatable:: tmp(:)
        integer(int32):: l

        l = size(ar)
        allocate(tmp(1:2*l))
        tmp(1:l) = ar(1:l)
        call move_alloc(tmp, ar)
    end subroutine

    subroutine heap_up(hm,ind)
        class(heap_map_int_int):: hm
        integer(int32),value:: ind
        integer(int32):: c

        do while(ind > 1)
            c = ind/2
            if (hm%key(ind) > hm%key(c)) return
            call kv_swap(hm, ind, c)
            ind=c
        end do
    end subroutine

    subroutine heap_down(hm,ind)
        class(heap_map_int_int):: hm
        integer(int32),value:: ind
        integer(int32):: c1,c2,c

        do while(ind*2 <= hm%len)
            c1 = ind*2; c2 = c1+1; c = c1
            if (c2 <= hm%len) c = merge(c1, c2, hm%key(c1) < hm%key(c2))
            if (hm%key(c) > hm%key(ind)) return
            call kv_swap(hm, c, ind)
            ind = c
        end do
    end subroutine

    subroutine hm_to_array(hm, k_ar, v_ar)
        class(heap_map_int_int):: hm
        integer(int64):: k_ar(hm%len), v_ar(hm%len)
        integer(int64):: i

        i=1
        do while(hm_remain(hm))
            call hm_pop(hm,k_ar(i), v_ar(i))
            i=i+1
        end do
    end subroutine

    subroutine double_heap_sort(k_ar, v_ar)
        type(heap_map_int_int):: hm
        integer(int64):: k_ar(:), v_ar(:)

        hm = heap_map_int_int(k_ar, v_ar)
        call hm_to_array(hm,k_ar,v_ar)
    end subroutine

    subroutine kv_swap(hm, i1, i2)
        type(heap_map_int_int),intent(inout):: hm
        integer(int32),intent(in):: i1,i2

        call swap(hm%key(i1),hm%key(i2))
        call swap(hm%val(i1),hm%val(i2))
    contains
        subroutine swap(x, y)
            integer(int64),intent(inout):: x,y
            x=xor(x,y)
            y=xor(x,y)
            x=xor(x,y)
        end subroutine
    end subroutine
end module