text:
module diffusion_module

  use amrex_fort_module, only : rt => amrex_real
  implicit none

  public

contains

  subroutine ca_tempdiffextrap(lo, hi, tdif, t_lo, t_hi) &
       bind(C, name="ca_tempdiffextrap")

    ! this routine extrapolates the temperature diffusion term into the
    ! ghostcells

    use prob_params_module, only: dg

    use amrex_fort_module, only : rt => amrex_real
    implicit none

    integer         , intent(in   ) :: lo(3), hi(3)
    integer         , intent(in   ) :: t_lo(3), t_hi(3)
    real(rt)        , intent(inout) :: tdif(t_lo(1):t_hi(1),t_lo(2):t_hi(2),t_lo(3):t_hi(3))

    ! Local variables

    integer :: i, j, k

    ! left side
    i = lo(1)-1*dg(1)
    do k = lo(3), hi(3)
       do j = lo(2), hi(2)
          tdif(i,j,k) = tdif(i+1*dg(1),j,k)
       end do
    end do

    ! right side
    i = hi(1)+1*dg(1)
    do k = lo(3), hi(3)
       do j = lo(2), hi(2)
          tdif(i,j,k) = tdif(i-1*dg(1),j,k)
       end do
    end do

    ! bottom side
    j = lo(2)-1*dg(2)
    do k = lo(3), hi(3)
       do i = lo(1), hi(1)
          tdif(i,j,k) = tdif(i,j+1*dg(2),k)
       end do
    end do

    ! top side
    j = hi(2)+1*dg(2)
    do k = lo(3), hi(3)
       do i = lo(1), hi(1)
          tdif(i,j,k) = tdif(i,j-1*dg(2),k)
       end do
    end do

    ! down side
    k = lo(3)-1*dg(3)
    do j = lo(2), hi(2)
       do i = lo(1), hi(1)
          tdif(i,j,k) = tdif(i,j,k+1*dg(3))
       end do
    end do

    ! up side
    k = hi(3)+1*dg(3)
    do j = lo(2), hi(2)
       do i = lo(1), hi(1)
          tdif(i,j,k) = tdif(i,j,k-1*dg(3))
       end do
    end do

    ! k-edges
    i = lo(1)-1*dg(1)
    j = lo(2)-1*dg(2)
    do k = lo(3), hi(3)
       tdif(i,j,k) = tdif(i+1*dg(1),j+1*dg(2),k)
    end do

    i = lo(1)-1*dg(1)
    j = hi(2)+1*dg(2)
    do k = lo(3), hi(3)
       tdif(i,j,k) = tdif(i+1*dg(1),j-1*dg(2),k)
    end do

    i = hi(1)+1*dg(1)
    j = lo(2)-1*dg(2)
    do k = lo(3), hi(3)
       tdif(i,j,k) = tdif(i-1*dg(1),j+1*dg(2),k)
    end do

    i = hi(1)+1*dg(1)
    j = hi(2)+1*dg(2)
    do k = lo(3), hi(3)
       tdif(i,j,k) = tdif(i-1*dg(1),j-1*dg(2),k)
    end do

    ! j-edges
    i = lo(1)-1*dg(1)
    k = lo(3)-1*dg(3)
    do j = lo(2), hi(2)
       tdif(i,j,k) = tdif(i+1*dg(1),j,k+1*dg(3))
    end do

    i = lo(1)-1*dg(1)
    k = hi(3)+1*dg(3)
    do j = lo(2), hi(2)
       tdif(i,j,k) = tdif(i+1*dg(1),j,k-1*dg(3))
    end do

    i = hi(1)+1*dg(1)
    k = lo(3)-1*dg(3)
    do j = lo(2), hi(2)
       tdif(i,j,k) = tdif(i-1*dg(1),j,k+1*dg(3))
    end do

    i = hi(1)+1*dg(1)
    k = hi(3)+1*dg(3)
    do j = lo(2), hi(2)
       tdif(i,j,k) = tdif(i-1*dg(1),j,k-1*dg(3))
    end do

    ! i-edges
    j = lo(2)-1*dg(2)
    k = lo(3)-1*dg(3)
    do i = lo(1), hi(1)
       tdif(i,j,k) = tdif(i,j+1*dg(2),k+1*dg(3))
    end do

    j = lo(2)-1*dg(2)
    k = hi(3)+1*dg(3)
    do i = lo(1), hi(1)
       tdif(i,j,k) = tdif(i,j+1*dg(2),k-1*dg(3))
    end do

    j = hi(2)+1*dg(2)
    k = lo(3)-1*dg(3)
    do i = lo(1), hi(1)
       tdif(i,j,k) = tdif(i,j-1*dg(2),k+1*dg(3))
    end do

    j = hi(2)+1*dg(2)
    k = hi(3)+1*dg(3)
    do i = lo(1), hi(1)
       tdif(i,j,k) = tdif(i,j-1*dg(2),k-1*dg(3))
    end do

    ! corners
    i = lo(1)-1*dg(1)
    j = lo(2)-1*dg(2)
    k = lo(3)-1*dg(3)
    tdif(i,j,k) = tdif(i+1*dg(1),j+1*dg(2),k+1*dg(3))

    i = lo(1)-1*dg(1)
    j = hi(2)+1*dg(2)
    k = lo(3)-1*dg(3)
    tdif(i,j,k) = tdif(i+1*dg(1),j-1*dg(2),k+1*dg(3))

    i = hi(1)+1*dg(1)
    j = hi(2)+1*dg(2)
    k = lo(3)-1*dg(3)
    tdif(i,j,k) = tdif(i-1*dg(1),j-1*dg(2),k+1*dg(3))

    i = hi(1)+1*dg(1)
    j = lo(2)-1*dg(2)
    k = lo(3)-1*dg(3)
    tdif(i,j,k) = tdif(i-1*dg(1),j+1*dg(2),k+1*dg(3))

    i = lo(1)-1*dg(1)
    j = lo(2)-1*dg(2)
    k = hi(3)+1*dg(3)
    tdif(i,j,k) = tdif(i+1*dg(1),j+1*dg(2),k-1*dg(3))

    i = lo(1)-1*dg(1)
    j = hi(2)+1*dg(2)
    k = hi(3)+1*dg(3)
    tdif(i,j,k) = tdif(i+1*dg(1),j-1*dg(2),k-1*dg(3))

    i = hi(1)+1*dg(1)
    j = lo(2)-1*dg(2)
    k = hi(3)+1*dg(3)
    tdif(i,j,k) = tdif(i-1*dg(1),j+1*dg(2),k-1*dg(3))

    i = hi(1)+1*dg(1)
    j = hi(2)+1*dg(2)
    k = hi(3)+1*dg(3)
    tdif(i,j,k) = tdif(i-1*dg(1),j-1*dg(2),k-1*dg(3))

  end subroutine ca_tempdiffextrap

  ! This routine fills the species diffusion coefficients on the edges of a zone
  ! by calling the cell-centered conductivity routine and averaging to
  ! the interfaces

  subroutine ca_fill_spec_coeff(lo,hi, &
       state,s_lo,s_hi, &
       coefx,cx_lo,cx_hi, &
       coefy,cy_lo,cy_hi, &
       coefz,cz_lo,cz_hi) &
       bind(C, name="ca_fill_spec_coeff")

    use bl_constants_module
    use network, only: nspec, naux
    use meth_params_module, only : NVAR, URHO, UTEMP, UEINT, UFS, UFX, diffuse_cutoff_density
    use prob_params_module, only : dg
    use conductivity_module
    use eos_type_module
    use eos_module, only : eos
    use amrex_fort_module, only : rt => amrex_real
    implicit none

    integer         , intent(in   ) :: lo(3), hi(3)
    integer         , intent(in   ) :: s_lo(3), s_hi(3)
    integer         , intent(in   ) :: cx_lo(3), cx_hi(3), cy_lo(3), cy_hi(3), cz_lo(3), cz_hi(3)
    real(rt)        , intent(in   ) :: state(s_lo(1):s_hi(1),s_lo(2):s_hi(2),s_lo(3):s_hi(3),NVAR)
    real(rt)        , intent(inout) :: coefx(cx_lo(1):cx_hi(1),cx_lo(2):cx_hi(2),cx_lo(3):cx_hi(3))
    real(rt)        , intent(inout) :: coefy(cy_lo(1):cy_hi(1),cy_lo(2):cy_hi(2),cy_lo(3):cy_hi(3))
    real(rt)        , intent(inout) :: coefz(cz_lo(1):cz_hi(1),cz_lo(2):cz_hi(2),cz_lo(3):cz_hi(3))

    ! local variables
    integer          :: i, j, k
    real(rt)         :: coef_cc(lo(1)-1:hi(1)+1,lo(2)-1:hi(2)+1,lo(3)-1:hi(3)+1)

    type (eos_t) :: eos_state
    real(rt)         :: coeff

    ! fill the cell-centered diffusion coefficient

    do k = lo(3)-1*dg(3),hi(3)+1*dg(3)
       do j = lo(2)-1*dg(2),hi(2)+1*dg(2)
          do i = lo(1)-1*dg(1),hi(1)+1*dg(1)
             eos_state%rho    = state(i,j,k,URHO)
             eos_state%T      = state(i,j,k,UTEMP)
             eos_state%e      = state(i,j,k,UEINT)/state(i,j,k,URHO)
             eos_state%xn(:)  = state(i,j,k,UFS:UFS-1+nspec)/ state(i,j,k,URHO)
             eos_state%aux(:) = state(i,j,k,UFX:UFX-1+naux)/ state(i,j,k,URHO)
             call eos(eos_input_re,eos_state)

             if (eos_state%rho > diffuse_cutoff_density) then
                call conductivity(eos_state, coeff)
                coeff = coeff / eos_state%cp
             else
                coeff = ZERO
             endif

             coef_cc(i,j,k) = coeff
          enddo
       enddo
    enddo

    ! average to the interfaces
    do k = lo(3),hi(3)
       do j = lo(2),hi(2)
          do i = lo(1),hi(1)+1*dg(1)
             coefx(i,j,k) = 0.5e0_rt * (coef_cc(i,j,k) + coef_cc(i-1*dg(1),j,k))
          end do
       end do
    enddo

    do k = lo(3),hi(3)
       do j = lo(2),hi(2)+1*dg(2)
          do i = lo(1),hi(1)
             coefy(i,j,k) = 0.5e0_rt * (coef_cc(i,j,k) + coef_cc(i,j-1*dg(2),k))
          end do
       end do
    enddo

    do k = lo(3),hi(3)+1*dg(3)
       do j = lo(2),hi(2)
          do i = lo(1),hi(1)
             coefz(i,j,k) = 0.5e0_rt * (coef_cc(i,j,k) + coef_cc(i,j,k-1*dg(3)))
          end do
       end do
    enddo

  end subroutine ca_fill_spec_coeff

  ! This routine fills the thermal conductivity on the edges of a zone
  ! by calling the cell-centered conductivity routine and averaging to
  ! the interfaces

  subroutine ca_fill_temp_cond(lo,hi, &
       state,s_lo,s_hi, &
       coefx,cx_lo,cx_hi, &
       coefy,cy_lo,cy_hi, &
       coefz,cz_lo,cz_hi) &
       bind(C, name="ca_fill_temp_cond")

    use bl_constants_module
    use network, only: nspec, naux
    use meth_params_module, only : NVAR, URHO, UTEMP, UEINt, UFS, UFX, &
         diffuse_cutoff_density, diffuse_cond_scale_fac, small_temp
    use prob_params_module, only : dg
    use conductivity_module
    use eos_type_module
    use eos_module, only : eos
    use amrex_fort_module, only : rt => amrex_real
    implicit none

    integer         , intent(in   ) :: lo(3), hi(3)
    integer         , intent(in   ) :: s_lo(3), s_hi(3)
    integer         , intent(in   ) :: cx_lo(3), cx_hi(3), cy_lo(3), cy_hi(3), cz_lo(3), cz_hi(3)
    real(rt)        , intent(in   ) :: state(s_lo(1):s_hi(1),s_lo(2):s_hi(2),s_lo(3):s_hi(3),NVAR)
    real(rt)        , intent(inout) :: coefx(cx_lo(1):cx_hi(1),cx_lo(2):cx_hi(2),cx_lo(3):cx_hi(3))
    real(rt)        , intent(inout) :: coefy(cy_lo(1):cy_hi(1),cy_lo(2):cy_hi(2),cy_lo(3):cy_hi(3))
    real(rt)        , intent(inout) :: coefz(cz_lo(1):cz_hi(1),cz_lo(2):cz_hi(2),cz_lo(3):cz_hi(3))

    ! local variables
    integer          :: i, j, k
    real(rt)         :: coef_cc(lo(1)-1:hi(1)+1,lo(2)-1:hi(2)+1,lo(3)-1:hi(3)+1)

    type (eos_t) :: eos_state
    real(rt)         :: cond

    ! fill the cell-centered conductivity

    do k = lo(3)-1*dg(3),hi(3)+1*dg(3)
       do j = lo(2)-1*dg(2),hi(2)+1*dg(2)
          do i = lo(1)-1*dg(1),hi(1)+1*dg(1)

             eos_state%rho    = state(i,j,k,URHO)
             eos_state%T      = state(i,j,k,UTEMP)   ! needed as an initial guess
             eos_state%e      = state(i,j,k,UEINT)/state(i,j,k,URHO)
             eos_state%xn(:)  = state(i,j,k,UFS:UFS-1+nspec)/ state(i,j,k,URHO)
             eos_state%aux(:) = state(i,j,k,UFX:UFX-1+naux)/ state(i,j,k,URHO)

             if (eos_state%e < ZERO) then
                eos_state%T = small_temp
                call eos(eos_input_rt,eos_state)
             else
                call eos(eos_input_re,eos_state)
             endif

             if (eos_state%rho > diffuse_cutoff_density) then

                call conductivity(eos_state, cond)
             else
                cond = ZERO
             endif
             coef_cc(i,j,k) = diffuse_cond_scale_fac * cond
          enddo
       enddo
    enddo

    ! average to the interfaces
    do k = lo(3),hi(3)
       do j = lo(2),hi(2)
          do i = lo(1),hi(1)+1*dg(1)
             coefx(i,j,k) = 0.5e0_rt * (coef_cc(i,j,k) + coef_cc(i-1*dg(1),j,k))
          end do
       end do
    enddo

    do k = lo(3),hi(3)
       do j = lo(2),hi(2)+1*dg(2)
          do i = lo(1),hi(1)
             coefy(i,j,k) = 0.5e0_rt * (coef_cc(i,j,k) + coef_cc(i,j-1*dg(2),k))
          end do
       end do
    enddo

    do k = lo(3),hi(3)+1*dg(3)
       do j = lo(2),hi(2)
          do i = lo(1),hi(1)
             coefz(i,j,k) = 0.5e0_rt * (coef_cc(i,j,k) + coef_cc(i,j,k-1*dg(3)))
          end do
       end do
    enddo

  end subroutine ca_fill_temp_cond

  ! This routine fills the coefficient of grad(enthalpy) on the edges of a zone
  ! by calling the cell-centered conductivity routine and averaging to
  ! the interfaces

  subroutine ca_fill_enth_cond(lo,hi, &
       state,s_lo,s_hi, &
       coefx,cx_lo,cx_hi, &
       coefy,cy_lo,cy_hi, &
       coefz,cz_lo,cz_hi) &
       bind(C, name="ca_fill_enth_cond")

    use bl_constants_module
    use network, only: nspec, naux
    use meth_params_module, only : NVAR, URHO, UTEMP, UEINT, UFS, UFX, diffuse_cutoff_density
    use prob_params_module, only : dg
    use conductivity_module
    use eos_type_module
    use eos_module, only : eos

    use amrex_fort_module, only : rt => amrex_real
    implicit none

    integer         , intent(in   ) :: lo(3), hi(3)
    integer         , intent(in   ) :: s_lo(3), s_hi(3)
    integer         , intent(in   ) :: cx_lo(3), cx_hi(3), cy_lo(3), cy_hi(3), cz_lo(3), cz_hi(3)
    real(rt)        , intent(in   ) :: state(s_lo(1):s_hi(1),s_lo(2):s_hi(2),s_lo(3):s_hi(3),NVAR)
    real(rt)        , intent(inout) :: coefx(cx_lo(1):cx_hi(1),cx_lo(2):cx_hi(2),cx_lo(3):cx_hi(3))
    real(rt)        , intent(inout) :: coefy(cy_lo(1):cy_hi(1),cy_lo(2):cy_hi(2),cy_lo(3):cy_hi(3))
    real(rt)        , intent(inout) :: coefz(cz_lo(1):cz_hi(1),cz_lo(2):cz_hi(2),cz_lo(3):cz_hi(3))

    ! local variables
    integer          :: i, j, k
    real(rt)         :: coef_cc(lo(1)-1:hi(1)+1,lo(2)-1:hi(2)+1,lo(3)-1:hi(3)+1)

    type (eos_t) :: eos_state
    real(rt)         :: cond

    ! fill the cell-centered conductivity

    do k = lo(3)-1*dg(3),hi(3)+1*dg(3)
       do j = lo(2)-1*dg(2),hi(2)+1*dg(2)
          do i = lo(1)-1*dg(1),hi(1)+1*dg(1)
             eos_state%rho    = state(i,j,k,URHO)
             eos_state%T      = state(i,j,k,UTEMP)
             eos_state%e      = state(i,j,k,UEINT)/state(i,j,k,URHO)
             eos_state%xn(:)  = state(i,j,k,UFS:UFS-1+nspec)/ state(i,j,k,URHO)
             eos_state%aux(:) = state(i,j,k,UFX:UFX-1+naux)/ state(i,j,k,URHO)
             call eos(eos_input_re,eos_state)

             if (eos_state%rho > diffuse_cutoff_density) then
                call conductivity(eos_state, cond)
                cond = cond / eos_state%cp
             else
                cond = ZERO
             endif

             coef_cc(i,j,k) = cond
          enddo
       enddo
    enddo

    ! average to the interfaces
    do k = lo(3),hi(3)
       do j = lo(2),hi(2)
          do i = lo(1),hi(1)+1*dg(1)
             coefx(i,j,k) = 0.5e0_rt * (coef_cc(i,j,k) + coef_cc(i-1*dg(1),j,k))
          end do
       end do
    enddo

    do k = lo(3),hi(3)
       do j = lo(2),hi(2)+1*dg(2)
          do i = lo(1),hi(1)
             coefy(i,j,k) = 0.5e0_rt * (coef_cc(i,j,k) + coef_cc(i,j-1*dg(2),k))
          end do
       end do
    enddo

    do k = lo(3),hi(3)+1*dg(3)
       do j = lo(2),hi(2)
          do i = lo(1),hi(1)
             coefz(i,j,k) = 0.5e0_rt * (coef_cc(i,j,k) + coef_cc(i,j,k-1*dg(3)))
          end do
       end do
    enddo

  end subroutine ca_fill_enth_cond

  ! This routine fills the viscous coefficients "mu" on the edges of a zone
  ! by calling the cell-centered coefficient routine and averaging to
  ! the interfaces

  subroutine ca_fill_first_visc_coeff(lo,hi, &
       state,s_lo,s_hi, &
       coefx,cx_lo,cx_hi, &
       coefy,cy_lo,cy_hi, &
       coefz,cz_lo,cz_hi) bind(C, name="ca_fill_first_visc_coeff")

    use bl_constants_module
    use network, only: nspec, naux
    use meth_params_module, only : NVAR, URHO, UTEMP, UFS, UFX, diffuse_cutoff_density
    use prob_params_module, only : dg
    use viscosity_module
    use eos_type_module

    use amrex_fort_module, only : rt => amrex_real
    implicit none

    integer         , intent(in   ) :: lo(3), hi(3)
    integer         , intent(in   ) :: s_lo(3), s_hi(3)
    integer         , intent(in   ) :: cx_lo(3), cx_hi(3), cy_lo(3), cy_hi(3), cz_lo(3), cz_hi(3)
    real(rt)        , intent(in   ) :: state(s_lo(1):s_hi(1),s_lo(2):s_hi(2),s_lo(3):s_hi(3),NVAR)
    real(rt)        , intent(inout) :: coefx(cx_lo(1):cx_hi(1),cx_lo(2):cx_hi(2),cx_lo(3):cx_hi(3))
    real(rt)        , intent(inout) :: coefy(cy_lo(1):cy_hi(1),cy_lo(2):cy_hi(2),cy_lo(3):cy_hi(3))
    real(rt)        , intent(inout) :: coefz(cz_lo(1):cz_hi(1),cz_lo(2):cz_hi(2),cz_lo(3):cz_hi(3))

    ! local variables
    integer          :: i, j, k
    real(rt)         :: coef_cc(lo(1)-1:hi(1)+1,lo(2)-1:hi(2)+1,lo(3)-1:hi(3)+1)

    type (eos_t) :: eos_state
    real(rt)         :: coeff

    ! fill the cell-centered viscous coefficient

    do k = lo(3)-1*dg(3),hi(3)+1*dg(3)
       do j = lo(2)-1*dg(2),hi(2)+1*dg(2)
          do i = lo(1)-1*dg(1),hi(1)+1*dg(1)
             eos_state%rho    = state(i,j,k,URHO)
             eos_state%T      = state(i,j,k,UTEMP)
             eos_state%xn(:)  = state(i,j,k,UFS:UFS-1+nspec)
             eos_state%aux(:) = state(i,j,k,UFX:UFX-1+naux)

             if (eos_state%rho > diffuse_cutoff_density) then
                call viscous_coeff(eos_state, coeff)
             else
                coeff = ZERO
             endif
             coef_cc(i,j,k) = 2.e0_rt * coeff
          enddo
       enddo
    enddo

    ! average to the interfaces
    do k = lo(3),hi(3)
       do j = lo(2),hi(2)
          do i = lo(1),hi(1)+1*dg(1)
             coefx(i,j,k) = 0.5e0_rt * (coef_cc(i,j,k) + coef_cc(i-1*dg(1),j,k))
          end do
       end do
    enddo

    do k = lo(3),hi(3)
       do j = lo(2),hi(2)+1*dg(2)
          do i = lo(1),hi(1)
             coefy(i,j,k) = 0.5e0_rt * (coef_cc(i,j,k) + coef_cc(i,j-1*dg(2),k))
          end do
       end do
    enddo

    do k = lo(3),hi(3)+1*dg(3)
       do j = lo(2),hi(2)
          do i = lo(1),hi(1)
             coefz(i,j,k) = 0.5e0_rt * (coef_cc(i,j,k) + coef_cc(i,j,k-1*dg(3)))
          end do
       end do
    enddo

  end subroutine ca_fill_first_visc_coeff

  subroutine ca_fill_secnd_visc_coeff(lo,hi, &
       state,s_lo,s_hi, &
       coefx,cx_lo,cx_hi, &
       coefy,cy_lo,cy_hi, &
       coefz,cz_lo,cz_hi) bind(C, name="ca_fill_secnd_visc_coeff")

    use bl_constants_module
    use network, only: nspec, naux
    use meth_params_module, only : NVAR, URHO, UTEMP, UFS, UFX, diffuse_cutoff_density
    use prob_params_module, only : dg
    use viscosity_module
    use eos_type_module

    use amrex_fort_module, only : rt => amrex_real
    implicit none

    integer         , intent(in   ) :: lo(3), hi(3)
    integer         , intent(in   ) :: s_lo(3), s_hi(3)
    integer         , intent(in   ) :: cx_lo(3), cx_hi(3), cy_lo(3), cy_hi(3), cz_lo(3), cz_hi(3)
    real(rt)        , intent(in   ) :: state(s_lo(1):s_hi(1),s_lo(2):s_hi(2),s_lo(3):s_hi(3),NVAR)
    real(rt)        , intent(inout) :: coefx(cx_lo(1):cx_hi(1),cx_lo(2):cx_hi(2),cx_lo(3):cx_hi(3))
    real(rt)        , intent(inout) :: coefy(cy_lo(1):cy_hi(1),cy_lo(2):cy_hi(2),cy_lo(3):cy_hi(3))
    real(rt)        , intent(inout) :: coefz(cz_lo(1):cz_hi(1),cz_lo(2):cz_hi(2),cz_lo(3):cz_hi(3))

    ! local variables
    integer          :: i, j, k
    real(rt)         :: coef_cc(lo(1)-1:hi(1)+1,lo(2)-1:hi(2)+1,lo(3)-1:hi(3)+1)

    type (eos_t) :: eos_state
    real(rt)         :: bulk_visc, mu, twothirds

    bulk_visc = 0.e0_rt
    twothirds = 2.e0_rt / 3.e0_rt

    ! fill the cell-centered viscous coefficient

    do k = lo(3)-1*dg(3),hi(3)+1*dg(3)
       do j = lo(2)-1*dg(2),hi(2)+1*dg(2)
          do i = lo(1)-1*dg(1),hi(1)+1*dg(1)
             eos_state%rho    = state(i,j,k,URHO)
             eos_state%T      = state(i,j,k,UTEMP)
             eos_state%xn(:)  = state(i,j,k,UFS:UFS-1+nspec)
             eos_state%aux(:) = state(i,j,k,UFX:UFX-1+naux)

             if (eos_state%rho > diffuse_cutoff_density) then
                call viscous_coeff(eos_state, mu)
             else
                mu = ZERO
             endif
             !          coef_cc(i,j,k) = coeff

             coef_cc(i,j,k) = (bulk_visc - twothirds*mu)

          enddo
       enddo
    enddo

    ! average to the interfaces
    do k = lo(3),hi(3)
       do j = lo(2),hi(2)
          do i = lo(1),hi(1)+1*dg(1)
             coefx(i,j,k) = 0.5e0_rt * (coef_cc(i,j,k) + coef_cc(i-1*dg(1),j,k))
          end do
       end do
    enddo

    do k = lo(3),hi(3)
       do j = lo(2),hi(2)+1*dg(2)
          do i = lo(1),hi(1)
             coefy(i,j,k) = 0.5e0_rt * (coef_cc(i,j,k) + coef_cc(i,j-1*dg(2),k))
          end do
       end do
    enddo

    do k = lo(3),hi(3)+1*dg(3)
       do j = lo(2),hi(2)
          do i = lo(1),hi(1)
             coefz(i,j,k) = 0.5e0_rt * (coef_cc(i,j,k) + coef_cc(i,j,k-1*dg(3)))
          end do
       end do
    enddo

  end subroutine ca_fill_secnd_visc_coeff

  subroutine ca_compute_div_tau_u(lo,hi,&
       div_tau_u,d_lo,d_hi, &
       state,s_lo,s_hi,dx,coord_type) bind(C, name="ca_compute_div_tau_u")

    use bl_constants_module
    use network, only: nspec, naux
    use meth_params_module, only : NVAR, URHO, UMX, UTEMP, UFS, UFX
    use viscosity_module
    use eos_type_module

    use amrex_fort_module, only : rt => amrex_real
    implicit none

    integer         , intent(in   ) ::   lo(3),   hi(3)
    integer         , intent(in   ) :: d_lo(3), d_hi(3)
    integer         , intent(in   ) :: s_lo(3), s_hi(3)
    real(rt)        , intent(  out) :: div_tau_u(d_lo(1):d_hi(1),d_lo(2):d_hi(2),d_lo(3):d_hi(3))
    real(rt)        , intent(in   ) :: state    (s_lo(1):s_hi(1),s_lo(2):s_hi(2),s_lo(3):s_hi(3),NVAR)
    real(rt)        , intent(in   ) :: dx(3)
    integer         , intent(in   ) :: coord_type

    ! local variables
    real(rt)         :: twothirds
    type (eos_t)     :: eos_state

    ! These are cell-centered
    real(rt)         :: rm,rc,rp,vm,vc,vp,mu_temp
    real(rt)        , allocatable :: mu(:), kappa(:)

    ! These are edge-centered
    real(rt)         :: mu_em, kap_em, tau1m, tau2m, x_em, v_em
    real(rt)         :: mu_ep, kap_ep, tau1p, tau2p, x_ep, v_ep

    integer          :: i, j, k, imin

    twothirds = 2.e0_rt / 3.e0_rt

    allocate(   mu(lo(1)-2:hi(1)+2))
    allocate(kappa(lo(1)-2:hi(1)+2))

    ! Fill cell-centered values of div(tau u) as a viscous source term for the energy equation

    ! NOTE THIS IS ONLY CORRECT IN 1-D SPHERICAL !!!!!!

    k = lo(3)
    j = lo(2)

    do i = lo(1)-2,hi(1)+2
       eos_state%rho    = state(i,j,k,URHO)
       eos_state%T      = state(i,j,k,UTEMP)
       eos_state%xn(:)  = state(i,j,k,UFS:UFS-1+nspec)
       eos_state%aux(:) = state(i,j,k,UFX:UFX-1+naux)

       call viscous_coeff(eos_state, mu_temp)
       mu(i) = mu_temp

       kappa(i) = 0.e0_rt
    end do

    imin = max(lo(1)-1,0)

    if (coord_type .eq. 2) then

       do i = imin,hi(1)+1

          ! These are cell-centered
          rm = (dble(i)-0.5e0_rt) * dx(1)
          rc = (dble(i)+0.5e0_rt) * dx(1)
          rp = (dble(i)+1.5e0_rt) * dx(1)

          ! These are cell-centered
          vp = state(i+1,j,k,UMX)/state(i+1,j,k,URHO)
          vc = state(i  ,j,k,UMX)/state(i  ,j,k,URHO)
          vm = state(i-1,j,k,UMX)/state(i-1,j,k,URHO)

          ! These are edge-centered
          mu_em = (mu(i-1)+mu(i)) * 0.5e0_rt
          mu_ep = (mu(i+1)+mu(i)) * 0.5e0_rt

          ! These are edge-centered
          kap_em = (kappa(i-1)+kappa(i)) * 0.5e0_rt
          kap_ep = (kappa(i+1)+kappa(i)) * 0.5e0_rt

          ! These are edge-centered
          tau1p = 2.e0_rt * mu_ep * (vp - vc) /dx(1)
          tau1m = 2.e0_rt * mu_em * (vc - vm) /dx(1)

          ! These are edge-centered
          tau2p = (kap_ep - twothirds*mu_ep) * (rp**2*vp - rc**2* vc) /dx(1)
          tau2m = (kap_em - twothirds*mu_em) * (rc**2*vc - rm**2* vm) /dx(1)

          ! These are edge-centered
          x_ep = dble(i+1)*dx(1)
          x_em = dble(i  )*dx(1)

          ! These are edge-centered
          v_ep = (vc + vp) * 0.5e0_rt
          v_em = (vc + vm) * 0.5e0_rt

          div_tau_u(i,j,k) = ( (x_ep**2*tau1p*v_ep) - (x_em**2*tau1m*v_em) &
               +(tau2p*v_ep)-(tau2m*v_em) ) / (rc**2 * dx(1))

       enddo

       if (lo(1) .eq. 0) &
            div_tau_u(lo(1)-1,j,k) = div_tau_u(lo(1),j,k)

    else

       call bl_abort("compute_div_tau_u currently hard-wired for 1-d spherical")

    end if

    deallocate(mu,kappa)

  end subroutine ca_compute_div_tau_u

end module diffusion_module


text:
	SUBROUTINE UPDGDG  ( iret )
C************************************************************************
C* UPDGDG								*
C* 									*
C* This subroutine updates the common areas for graphs in the G 	*
C* coordinate system.							*
C* 									*
C* UPDGDG ( IRET )							*
C* 									*
C* Output parameters:							*
C*	IRET		INTEGER		Return code			*
C**									*
C* Log:									*
C* M. desJardins/GSFC	9/85	GEMPLT Version 3.1			*
C* M. desJardins/GSFC	6/86	Corrected computations of all scalings	*
C* M. desJardins/GSFC	10/86	Added polar coordinate system		*
C* M. desJardins/GSFC	 6/88	Documentation				*
C************************************************************************
	INCLUDE		'DEVCHR.CMN'
	INCLUDE		'XYDEF.CMN'
	INCLUDE		'GEMPRM.PRM'
C-------------------------------------------------------------------------
C*	Do computations for graph mode.
C
	IF  ( ( igmode .eq. 2 ) .and. ggset ) THEN
C*
	    IF  ( ( jgytyp .eq. 1 ) .or. ( jgytyp .eq. 5 ) ) THEN
		yb = gybmg
		yt = gytmg
	      ELSE IF  ( jgytyp .eq. 2 )  THEN
		yb = ALOG ( gybmg )
		yt = ALOG ( gytmg )
	      ELSE IF  ( jgytyp .eq. 3 )  THEN
		yb = gybmg ** RKAPPA
		yt = gytmg ** RKAPPA
	    END IF
C
	    IF  ( ( jgxtyp .eq. 1 ) .or. ( jgxtyp .eq. 5 ) ) THEN
		xl = gxlmg
		xr = gxrmg
	      ELSE IF  ( jgxtyp .eq. 2 )  THEN
		xl = ALOG ( gxlmg )
		xr = ALOG ( gxrmg )
	      ELSE IF  ( jgxtyp .eq. 3 )  THEN
		xl = gxlmg ** RKAPPA
		xr = gxrmg ** RKAPPA
	    END IF
C
C*	    Save grid to linear ( 0. - 1. ) scaling terms for grid coord sys.
C
	    gxx1 = 1.0 / ( gpxr - gpxl )
	    gxx0 =   - ( gpxl * gxx1 )
	    gyy1 = 1.0 / ( gpyt - gpyb )
	    gyy0 =   - ( gpyb * gyy1 )
C
C*	    Compute scaling from (0.-1.) grid sys to values at points.
C
	    gglx1 = 1.0 / ( xr - xl )
	    gglx0 = - (gglx1 * xl)
	    ggly1 = 1.0 / ( yt - yb )
	    ggly0 = - ( ggly1 * yb )
	END IF
C*
	iret = 0
C*
	RETURN
	END

ext: f

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
program cdn_bhz_2d
   USE CDMFT_ED !
   USE SCIFOR
   USE DMFT_TOOLS
   !
   USE MPI
   !
   implicit none
   integer                                                                :: Nx,Ny,Nlso,iloop,Nb,Nkx,iw,iii,jjj,kkk
   integer,dimension(2):: recover
   logical                                                                :: converged
   real(8)                                                                :: ts,Mh,lambda,wmixing,observable
   !Bath:
   real(8),allocatable                                                    :: Bath(:),Bath_prev(:)
   !The local hybridization function:
   complex(8),allocatable                                                 :: Hloc(:,:)
   complex(8),allocatable,dimension(:,:,:,:,:,:,:)                        :: Gmats,Greal,Smats,Sreal,Weiss,Weiss_old
   character(len=16)                                                      :: finput
   real(8),allocatable                                                    :: wt(:)
   complex(8),allocatable                                                 :: wm(:),wr(:)
   complex(8),allocatable                                                 :: Hk(:,:,:),Smats_lso(:,:,:)
   complex(8),dimension(:,:,:,:,:,:),allocatable                          :: observable_matrix
   !SYMMETRIES TEST
   real(8),dimension(:),allocatable                                       :: lambdasym_vector
   complex(8),dimension(:,:,:,:,:,:,:),allocatable                        :: Hsym_basis
   !MPI VARIABLES (local use -> ED code has its own set of MPI variables)
   integer                                                                :: comm
   integer                                                                :: rank
   integer                                                                :: mpi_size
   logical                                                                :: master,hermiticize

   !Init MPI: use of MPI overloaded functions in SciFor
   call init_MPI(comm,.true.)
   rank   = get_Rank_MPI(comm)
   master = get_Master_MPI(comm)

   !

   !Parse input variables
   call parse_cmd_variable(finput,"FINPUT",default='inputED.conf')
   call parse_input_variable(wmixing,"wmixing",finput,default=1.d0,comment="Mixing bath parameter")
   call parse_input_variable(ts,"TS",finput,default=0.25d0,comment="hopping parameter")
   call parse_input_variable(Mh,"Mh",finput,default=1.d0,comment="crystal field splitting")
   call parse_input_variable(lambda,"lambda",finput,default=0.3d0,comment="spin-orbit coupling")
   call parse_input_variable(Nx,"Nx",finput,default=2,comment="Number of cluster sites in x direction")
   call parse_input_variable(Ny,"Ny",finput,default=2,comment="Number of cluster sites in y direction")
   call parse_input_variable(Nkx,"Nkx",finput,default=10,comment="Number of kx point for BZ integration")
   call parse_input_variable(hermiticize,"HERMITICIZE",finput,default=.true.,comment="are bath replicas hermitian")
   !
   call ed_read_input(trim(finput),comm)
   !
   !Add dmft control variables
   !
   call add_ctrl_var(beta,"BETA")
   call add_ctrl_var(Norb,"Norb")
   call add_ctrl_var(Nspin,"Nspin")
   call add_ctrl_var(xmu,"xmu")
   call add_ctrl_var(wini,"wini")
   call add_ctrl_var(wfin,"wfin")
   call add_ctrl_var(eps,"eps")

   !set global variables
   !if (Nspin/=1.or.Norb/=1) stop "You are using too many spin-orbitals"
   Ny=1
   !Nky=Nkx
   Nlat=Nx*Ny
   Nlso=Nlat*Nspin*Norb
   if(.not.allocated(wm))allocate(wm(Lmats))
   wm     = xi*pi/beta*real(2*arange(1,Lmats)-1,8)

   if(ED_VERBOSE > 0)call naming_convention()
   !Allocate Weiss Field:
   allocate(Weiss(Nlat,Nlat,Nspin,Nspin,Norb,Norb,Lmats),Weiss_old(Nlat,Nlat,Nspin,Nspin,Norb,Norb,Lmats))
   allocate(Gmats(Nlat,Nlat,Nspin,Nspin,Norb,Norb,Lmats),Greal(Nlat,Nlat,Nspin,Nspin,Norb,Norb,Lreal))
   allocate(Smats(Nlat,Nlat,Nspin,Nspin,Norb,Norb,Lmats),Sreal(Nlat,Nlat,Nspin,Nspin,Norb,Norb,Lreal))
   allocate(Smats_lso(Nlso,Nlso,Lmats))

   !Build Hk and Hloc
   call generate_hk_hloc()

   !CUSTOM OBSERVABLES: n1, n2, Ekin
   allocate(observable_matrix(Nlat,Nlat,Nspin,Nspin,Norb,Norb))
   call init_custom_observables(3,Hk)
   observable_matrix=zero
   do iii=1,Nlat
      observable_matrix(iii,iii,1,1,1,1)=one/Nlat
      observable_matrix(iii,iii,Nspin,Nspin,1,1)=one/Nlat
   enddo
   call add_custom_observable("n1",nnn2lso(observable_matrix))
   observable_matrix=zero
   do iii=1,Nlat
      observable_matrix(iii,iii,1,1,2,2)=one/Nlat
      observable_matrix(iii,iii,Nspin,Nspin,2,2)=one/Nlat
   enddo
   call add_custom_observable("n2",nnn2lso(observable_matrix))
   call add_custom_observable("Ekin",Hk)

   !SETUP BATH STEP 1
   allocate(lambdasym_vector(3))
   allocate(Hsym_basis(Nlat,Nlat,Nspin,Nspin,Norb,Norb,3))
   !
   lambdasym_vector(1)=Mh
   Hsym_basis(:,:,:,:,:,:,1)=lso2nnn(hloc_model(Nlso,1.d0,0.d0,0.d0))
   !
   lambdasym_vector(2)=ts
   Hsym_basis(:,:,:,:,:,:,2)=lso2nnn(hloc_model(Nlso,0.d0,1.d0,0.d0))
   !
   lambdasym_vector(3)=lambda
   Hsym_basis(:,:,:,:,:,:,3)=lso2nnn(hloc_model(Nlso,0.d0,0.d0,1.d0))
   !
   !SETUP BATH STEP 2 and SETUP SOLVER
   call ed_set_Hreplica(Hsym_basis,lambdasym_vector)
   Nb=ed_get_bath_dimension(Hsym_basis)
   allocate(bath(Nb))
   allocate(bath_prev(Nb))
   call ed_init_solver(comm,bath)

   !DMFT loop
   iloop=0;converged=.false.
   do while(.not.converged.AND.iloop<nloop)
      iloop=iloop+1
      if(master)call start_loop(iloop,nloop,"DMFT-loop")

      !Solve the EFFECTIVE IMPURITY PROBLEM (first w/ a guess for the bath)
      call ed_solve(comm,bath,lso2nnn(Hloc))
      call ed_get_sigma_matsubara(Smats)
      call ed_get_sigma_realaxis(Sreal)

      !Compute the local gfs:
      call dmft_gloc_matsubara(Hk,Gmats,Smats)
      if(master)call dmft_print_gf_matsubara(Gmats,"Gloc",iprint=4)
      !
      !Get the Weiss field/Delta function to be fitted
      call dmft_self_consistency(Gmats,Smats,Weiss,lso2nnn(Hloc),cg_scheme)
      call Bcast_MPI(comm,Weiss)
      !
      !MIXING:
      !if(iloop>1)Weiss = wmixing*Weiss + (1.d0-wmixing)*Weiss_Old
      !Weiss_old=Weiss
      !
      !Perform the SELF-CONSISTENCY by fitting the new bath
      if(master)then
         call ed_chi2_fitgf(Weiss,bath)
         !
         !MIXING:
         !call adaptive_mix(Bath(Nbath+1:),Bath_fitted(Nbath+1:)-Bath(Nbath+1:),wmixing,iloop)
         if(iloop>1)Bath = wmixing*Bath + (1.d0-wmixing)*Bath_Prev
         Bath_Prev=Bath
         !
         !Check convergence (if required change chemical potential)
         converged = check_convergence(Weiss(:,:,1,1,1,1,:),dmft_error,nsuccess,nloop)
      endif
      !
      call Bcast_MPI(comm,bath)
      call Bcast_MPI(comm,converged)
      call Bcast_MPI(comm,xmu)
      !
      if(master)call end_loop
   enddo

   !Compute the local gfs:
   call dmft_gloc_realaxis(Hk,Greal,Sreal)
   if(master)call dmft_print_gf_realaxis(Greal,"Gloc",iprint=4)

   !Compute the Kinetic Energy:
   do iw=1,Lmats
      Smats_lso(:,:,iw)=nnn2lso(Smats(:,:,:,:,:,:,iw))
   enddo
   call dmft_kinetic_energy(Hk(:,:,:),Smats_lso)

   call finalize_MPI()

contains

   !+------------------------------------------------------------------+
   !PURPOSE  : Hloc for the 2d BHZ model
   !+------------------------------------------------------------------+

   function Hloc_model(N,Mh_,ts_,lambda_) result (H0)
      integer                                               :: N,ilat,jlat,ispin,iorb,jorb,ind1,ind2
      real(8)                                               :: Mh_,ts_,lambda_
      complex(8),dimension(Nlat,Nlat,Nspin,Nspin,Norb,Norb) :: hopping_matrix
      complex(8),dimension(N,N)                             :: H0
      !
      hopping_matrix=zero
      !
      do ispin=1,Nspin
         do ilat=1,Nx
            do jlat=1,Ny
               ind1=indices2N([ilat,jlat])
               hopping_matrix(ind1,ind1,ispin,ispin,:,:)= t_m(mh_)
               if(ilat<Nx)then
                  ind2=indices2N([ilat+1,jlat])
                  hopping_matrix(ind2,ind1,ispin,ispin,:,:)= t_x(ts_,lambda_,ispin)
               endif
               if(ilat>1)then
                  ind2=indices2N([ilat-1,jlat])
                  hopping_matrix(ind2,ind1,ispin,ispin,:,:)= dconjg(transpose(t_x(ts_,lambda_,ispin)))
               endif
               if(jlat<Ny)then
                  ind2=indices2N([ilat,jlat+1])
                  hopping_matrix(ind2,ind1,ispin,ispin,:,:)= t_y(ts_,lambda_)
               endif
               if(jlat>1)then
                  ind2=indices2N([ilat,jlat-1])
                  hopping_matrix(ind2,ind1,ispin,ispin,:,:)= transpose(t_y(ts_,lambda_))
               endif
            enddo
         enddo
      enddo
      !
      H0=nnn2lso(hopping_matrix)
      !
   end function hloc_model

   function hk_model(kpoint,N) result(Hk)
      integer                                                      :: N,ilat,jlat,ispin,iorb,jorb,i,j,ind1,ind2
      real(8),dimension(:)                                         :: kpoint
      complex(8),dimension(Nlat,Nlat,Nspin,Nspin,Norb,Norb)        :: hopping_matrix
      complex(8),dimension(N,N)                                    :: hk
      !
      hopping_matrix=zero
      !
      do ispin=1,Nspin
         do ilat=1,Ny
            ind1=indices2N([1,ilat])
            ind2=indices2N([Nx,ilat])
            hopping_matrix(ind2,ind1,ispin,ispin,:,:)=hopping_matrix(ind2,ind1,ispin,ispin,:,:) + dconjg(transpose(t_x(ts,lambda,ispin)))*exp(xi*kpoint(1)*Nx)
            hopping_matrix(ind1,ind2,ispin,ispin,:,:)=hopping_matrix(ind1,ind2,ispin,ispin,:,:) + t_x(ts,lambda,ispin)*exp(-xi*kpoint(1)*Nx)
         enddo
      enddo
      !
      Hk=nnn2lso(hopping_matrix)+hloc_model(N,Mh,ts,lambda)
      !
   end function hk_model

   !AUXILLIARY HOPPING MATRIX CONSTRUCTORS

   function t_m(mass) result(tmpmat)
      complex(8),dimension(Norb,Norb) :: tmpmat
      real(8)                         :: mass
      !
      tmpmat=zero
      tmpmat=mass*pauli_sigma_z
      !
   end function t_m

   function t_x(hop1,hop2,spinsign) result(tmpmat)
      complex(8),dimension(Norb,Norb) :: tmpmat
      real(8)                         :: hop1,hop2,sz
      integer                         :: spinsign
      !
      tmpmat=zero
      sz=(-1.d0)**(spinsign+1)
      tmpmat=-hop1*pauli_sigma_z+0.5d0*sz*xi*hop2*pauli_sigma_x
      !
   end function t_x

   function t_y(hop1,hop2) result(tmpmat)
      complex(8),dimension(Norb,Norb) :: tmpmat
      real(8)                         :: hop1,hop2
      !
      tmpmat=zero
      tmpmat=-hop1*pauli_sigma_z
      tmpmat(1,2)=-hop2*0.5d0
      tmpmat(2,1)=hop2*0.5d0
      !
   end function t_y

   !-------------------------------------------------------------------------------------------
   !PURPOSE: generate Hloc and Hk
   !-------------------------------------------------------------------------------------------

   subroutine generate_hk_hloc()
      integer                                     :: ik
      real(8),dimension(Nkx,2)                    :: kgrid
      real(8),dimension(2)                        :: e1,e2,bk1,bk2
      real(8)                                     :: bklen
      !
      e1 = [1d0, 0d0]
      e2 = [0d0, 1d0]
      call TB_set_ei(eix=e1,eiy=e2)
      bklen=2d0*pi
      bk1=bklen*[1d0, 0d0]
      bk2=bklen*[0d0, 1d0]
      call TB_set_bk(bkx=bk1,bky=bk2)
      !
      call TB_build_kgrid([Nkx,1],kgrid)
      kgrid(:,1)=kgrid(:,1)/Nx
      !
      if(allocated(hk))deallocate(hk)
      if(allocated(Wt))deallocate(Wt)
      if(allocated(hloc))deallocate(Hloc)
      !
      allocate(Hk(Nlso,Nlso,Nkx),Wt(Nkx),Hloc(Nlso,Nlso))
      hk=zero
      wt=zero
      hloc=zero
      !
      ! SEVER !
      call TB_build_model(Hk,hk_model,Nlso,kgrid)
      Wt = 1d0/(Nkx)
      Hloc=hloc_model(Nlso,Mh,ts,lambda)
      !
   end subroutine generate_hk_hloc

   !+------------------------------------------------------------------+
   !PURPOSE  : Auxilliary reshape functions
   !+------------------------------------------------------------------+

   function lso2nnn(Hlso) result(Hnnn)
      complex(8),dimension(Nlat*Nspin*Norb,Nlat*Nspin*Norb) :: Hlso
      complex(8),dimension(Nlat,Nlat,Nspin,Nspin,Norb,Norb) :: Hnnn
      integer                                               :: ilat,jlat
      integer                                               :: iorb,jorb
      integer                                               :: ispin,jspin
      integer                                               :: is,js
      Hnnn=zero
      do ilat=1,Nlat
         do jlat=1,Nlat
            do ispin=1,Nspin
               do jspin=1,Nspin
                  do iorb=1,Norb
                     do jorb=1,Norb
                        is = iorb + (ilat-1)*Norb + (ispin-1)*Norb*Nlat
                        js = jorb + (jlat-1)*Norb + (jspin-1)*Norb*Nlat
                        Hnnn(ilat,jlat,ispin,jspin,iorb,jorb) = Hlso(is,js)
                     enddo
                  enddo
               enddo
            enddo
         enddo
      enddo
   end function lso2nnn

   function nnn2lso(Hnnn) result(Hlso)
      complex(8),dimension(Nlat,Nlat,Nspin,Nspin,Norb,Norb) :: Hnnn
      complex(8),dimension(Nlat*Nspin*Norb,Nlat*Nspin*Norb) :: Hlso
      integer                                               :: ilat,jlat
      integer                                               :: iorb,jorb
      integer                                               :: ispin,jspin
      integer                                               :: is,js
      Hlso=zero
      do ilat=1,Nlat
         do jlat=1,Nlat
            do ispin=1,Nspin
               do jspin=1,Nspin
                  do iorb=1,Norb
                     do jorb=1,Norb
                        is = iorb + (ilat-1)*Norb + (ispin-1)*Norb*Nlat
                        js = jorb + (jlat-1)*Norb + (jspin-1)*Norb*Nlat
                        Hlso(is,js) = Hnnn(ilat,jlat,ispin,jspin,iorb,jorb)
                     enddo
                  enddo
               enddo
            enddo
         enddo
      enddo
   end function nnn2lso

   function indices2N(indices) result(N)
      integer,dimension(2)         :: indices
      integer                      :: N,i
      !
      !
      N=Nx*(indices(2)-1)+indices(1)
   end function indices2N

   function N2indices(N) result(indices)
      integer,dimension(2)         :: indices
      integer                      :: N,i
      !
      indices(1)=mod(N,Nx)
      if(indices(1)==0)then
         indices(1)=Nx
         indices(2)=(N-Nx)/Nx+1
      else
         indices(2)=N/Nx+1
      endif
   end function N2indices

   subroutine naming_convention()
      integer                       :: i,j
      integer,dimension(Nx,Ny)      :: matrix
      !
      do j=1,Ny
         do i=1,Nx
            matrix(i,j)=indices2N([i,j])
         enddo
      enddo
      !
      write(LOGfile,"(A)")"The unique index of each site (on the cartesian plane) is as follows:"
      write(LOGfile,"(A)")" "
      do j=1,Ny
         write(LOGfile,"(20(I2,2x))")(matrix(i,Ny+1-j),i =1,Nx)
      enddo
      write(LOGfile,"(A)")" "
   end subroutine naming_convention
   !
end program cdn_bhz_2d

max_stars_repo_name: QcmPlab/CDMFT-LANC-ED

max_issues_repo_name: QcmPlab/CDMFT-LANC-ED

max_forks_repo_name: QcmPlab/CDMFT-LANC-ED


text:
! @@name:	flush_nolist.1f
! @@type:	F-fixed
! @@compilable:	yes
! @@linkable:	yes
! @@expect:	success
      SUBROUTINE F1(Q)
        COMMON /DATA/ X, P
        INTEGER, TARGET  :: X
        INTEGER, POINTER :: P
        INTEGER Q

        Q = 1
!$OMP   FLUSH
        ! X, P and Q are flushed
        ! because they are shared and accessible
      END SUBROUTINE F1

      SUBROUTINE F2(Q)
        COMMON /DATA/ X, P
        INTEGER, TARGET  :: X
        INTEGER, POINTER :: P
        INTEGER Q

!$OMP   BARRIER
          Q = 2
!$OMP   BARRIER
          ! a barrier implies a flush
          ! X, P and Q are flushed
          ! because they are shared and accessible
      END SUBROUTINE F2

      INTEGER FUNCTION G(N)
        COMMON /DATA/ X, P
        INTEGER, TARGET  :: X
        INTEGER, POINTER :: P
        INTEGER N
        INTEGER I, J, SUM

        I = 1
        SUM = 0
        P = 1
!$OMP   PARALLEL REDUCTION(+: SUM) NUM_THREADS(10)
          CALL F1(J)
            ! I, N and SUM were not flushed
            !   because they were not accessible in F1
            ! J was flushed because it was accessible
          SUM = SUM + J

          CALL F2(J)
            ! I, N, and SUM were not flushed
            !   because they were not accessible in f2
            ! J was flushed because it was accessible
          SUM = SUM + I + J + P + N
!$OMP   END PARALLEL

        G = SUM
      END FUNCTION G

      PROGRAM FLUSH_NOLIST
        COMMON /DATA/ X, P
        INTEGER, TARGET  :: X
        INTEGER, POINTER :: P
        INTEGER RESULT, G

        P => X
        RESULT = G(7)
        PRINT *, RESULT
      END PROGRAM FLUSH_NOLIST

max_stars_repo_name: kammerdienerb/flang

max_issues_repo_name: kammerdienerb/flang

max_forks_repo_name: kammerdienerb/flang


text:
! GenMask.f90
!
! Author : Joseph Schoonover
! E-mail : jschoonover@lanl.gov, schoonover.numerics@gmail.com
!
! Copyright 2017 Joseph Schoonover, Los Alamos National Laboratory
!
! Redistribution and use in source and binary forms, with or without
! modification,
! are permitted provided that the following conditions are met:
!
! 1. Redistributions of source code must retain the above copyright notice, this
! list of conditions and the following disclaimer.
!
! 2. Redistributions in binary form must reproduce the above copyright notice,
! this list of conditions and the following disclaimer in the
! documentation and/or other materials provided with the distribution.
!
! 3. Neither the name of the copyright holder nor the names of its contributors
! may be used to endorse or promote products derived from this
!  software without specific prior written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
! AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR
! CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
! EXEMPLARY,  OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
! PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR
! BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
! LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
! NEGLIGENCE  OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
! SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
!
! ////////////////////////////////////////////////////////////////////////////////////////////////

PROGRAM GenMask

! src/common/
USE CommonRoutines
! src/POP/
USE POP_Params_Class
USE POP_Mesh_Class

IMPLICIT NONE

   TYPE( POP_Params )   :: params
   TYPE( POP_Mesh )     :: mesh
   INTEGER              :: i, j
   INTEGER, ALLOCATABLE :: maskField(:,:)
   CHARACTER(400)       :: ncfile
   REAL(prec)           :: x, y, r

      CALL params % Build( )

      CALL mesh % Load( TRIM(params % meshFile)  )

      ALLOCATE( maskField(1:mesh % nX,1:mesh % nY) )

      maskField = 0

      DO j = 1, mesh % nY
         DO i = 1, mesh % nX

            x = mesh % tLon(i,j)
            y = mesh % tLat(i,j)

            IF( x >= 180.0_prec )THEN
               x = x -360.0_prec
            ENDIF
            ! Build a circular region around the Agulhas
            r = sqrt( (x-20.0_prec)**2 + (y+40.0_prec)**2 )
            IF( r <= 30.0_prec )THEN
               IF( r > 29.5_prec )THEN
                  maskfield(i,j) = -1 ! Prescribed Points
               ELSE
                  maskfield(i,j) = 1  ! Interior Points
               ENDIF
            ENDIF

         ENDDO
      ENDDO

      CALL WriteMaskField( mesh, maskField, TRIM(params % maskFile) )
      CALL mesh % Trash( )

CONTAINS
 SUBROUTINE WriteMaskField( mesh, maskfield, maskfile )

   IMPLICIT NONE
   TYPE( POP_Mesh ), INTENT(inout)      :: mesh
   INTEGER, INTENT(in)                  :: maskfield(1:mesh % nX, 1:mesh % nY)
   CHARACTER(*), INTENT(in)             :: maskfile
   ! Local
   INTEGER :: start(1:2), recCount(1:2)
   INTEGER :: ncid, varid, x_dimid, y_dimid

      start    = (/1, 1/)
      recCount = (/mesh % nX, mesh % nY/)

      CALL Check( nf90_create( PATH=TRIM(maskfile),&
                               CMODE=OR(nf90_clobber,nf90_64bit_offset),&
                               NCID=ncid ) )
      CALL Check( nf90_def_dim( ncid, "nlon", mesh % nX, x_dimid ) )
      CALL Check( nf90_def_dim( ncid, "nlat", mesh % nY, y_dimid ) )
      CALL Check( nf90_def_var( ncid, "mask",NF90_INT,&
                               (/ x_dimid, y_dimid /),&
                                varid ) )

      CALL Check( nf90_put_att( ncid, varid, "long_name", "Domain Mask" ) )
      CALL Check( nf90_put_att( ncid, varid, "units", "" ) )
!      CALL Check( nf90_put_att( ncid, varid, "_FillValue", fillValue) )
!      CALL Check( nf90_put_att( ncid, varid, "missing_value", fillValue) )

      CALL Check( nf90_enddef(ncid) )

      CALL Check( nf90_put_var( ncid, &
                                varid, &
                                maskfield, &
                                start, recCount ) )

      CALL Check( nf90_close( ncid ) )

 END SUBROUTINE WriteMaskField

END PROGRAM GenMask

max_stars_repo_name: lanl/FEOTS

max_issues_repo_name: lanl/FEOTS

max_forks_repo_name: lanl/FEOTS


text:
      subroutine FracCheckElemNearestCrack(y, z, ix, id, u, ele, nd, cflag)
!!!!!!!!            FracCheckElemNearestCrack(y, z, ix, id, u, ele, ndcircle, cflag)
!!!!!c     compare distance from current element to cracked elment
!!!!!c     with preset distance
!!!
!!!      parameter (nume=40000)
!!!	  common/bk06/nprnt,mprint,itmpop,numelt,jprint,idump,locstr
!!!	  common /stressflag/ ElemFractCode(nume),ElemDecayCount(nume)
!!!	  common/cracktip/ tipele(2,nume), ntp
!!!
!!!      dimension y(*), z(*), ix(4,*), id(2,*), u(*)
!!!      real y,z,u
!!!      integer ix,id
!!!
!!!	  integer tipele, ntp, ele, nd, cflag
!!!	  real y3, z3, y1, z1, ds0, ds
!!!	  real ye, ze, yt, zt
!!!	  integer ntel, ElemFractCode, numelt
!!!
!!!	  cflag=1
!!!
!!!	  y3=y(ix(3,1))
!!!	  z3=z(ix(3,1))
!!!	  y1=y(ix(1,1))
!!!	  z1=z(ix(1,1))
!!!	  ds0=sqrt((y3-y1)**2+(z3-z1)**2)*nd   ! nd is a coefficient
!!!
!!!	  ye=(y(ix(1,ele))+y(ix(3,ele)))/2.0
!!!	  ze=(z(ix(1,ele))+z(ix(3,ele)))/2.0
!!!
!!!	  !-- if it's on the crack tip of another existing crack
!!!	  do i=1, ntp
!!!	      if(ele==tipele(2,i)) then
!!!		      go to 30
!!!		  end if
!!!	  end do
!!!
!!!!c	  do i=1, ntp
!!!!c	      if(ele/=tipele(2,i)) then   !element not at the crack tip
!!!!c		      ntel=tipele(1,i)
!!!!c		      yt=(y(ix(1,ntel))+u(id(1,ix(1,ntel)))
!!!!c     >            +y(ix(3,ntel))+u(id(1,ix(3,ntel))))/2.0
!!!!c	          zt=(z(ix(1,ntel))+u(id(2,ix(1,ntel)))
!!!!c     >            +z(ix(3,ntel))+u(id(2,ix(3,ntel))))/2.0
!!!!c	          ds=sqrt((ye-yt)**2+(ze-zt)**2)
!!!!c			  if(ds<ds0) then
!!!!c			      cflag=0
!!!!c			  end if
!!!!c		  end if
!!!!c	  end do
!!!
!!!	  do i=1,numelt
!!!	      if (ElemFractCode(i) > 0) then ! if it's decaying or decayed .. it's cracking
!!!		      ntel=i
!!!			  yt=(y(ix(1,ntel))+y(ix(3,ntel)))/2.0
!!!	          zt=(z(ix(1,ntel))+z(ix(3,ntel)))/2.0
!!!	          ds=sqrt((ye-yt)**2+(ze-zt)**2)
!!!			  if(ds<ds0) then   !- if this elem i
!!!			      cflag=0
!!!				  go to 30
!!!			  end if
!!!		  end if
!!!	  end do
!!!
   30 end

max_stars_repo_name: imohame/LabCode

max_issues_repo_name: imohame/LabCode

max_forks_repo_name: imohame/LabCode