text:
module mkscalforce_module
  !
  ! this module contains the 2d and 3d routines that make the forcing
  ! terms for the scalar equations.
  !
  ! mkrhohforce computes the  wtilde dp/dr + psi source for rho*h evolution
  !
  ! mktempforce computes the source terms that appear in the
  ! temperature evolution equation.  Note, this formulation is only
  ! used in the prediction of the temperature on the edges, when
  ! enthalpy_pred_type = predict_T_then_rhohprime or predict_T_then_h.
  ! The edge temperatures are then converted to enthalpy for the full update.
  !

  use bl_constants_module
  use multifab_module
  use ml_layout_module
  use define_bc_module
  use ml_restrict_fill_module
  use pred_parameters
  use probin_module, only: enthalpy_pred_type

  implicit none

  private

  public :: mkrhohforce, mkhprimeforce, mktempforce

contains

  subroutine mkrhohforce(mla, scal_force, is_prediction, thermal, umac, &
                         p0_1, p0_2, rho0_1, rho0_2, &
                         psi, dx, add_thermal, the_bc_level)

    ! note, in the prediction of the interface states, we will set
    ! both p0/rho0_1 and p0/rho0_2 to the same old value.  In the computation
    ! of the rhoh_force for the update, they will be used to time-center.

    use bl_prof_module
    use variables, only: foextrap_comp, rhoh_comp
    use geometry, only: spherical, nr_fine, nlevs_radial
    use fill_3d_module
    use make_grav_module

    type(multifab) , intent(inout) :: scal_force(:)
    logical        , intent(in   ) :: is_prediction
    type(multifab) , intent(in   ) :: thermal(:)
    type(multifab) , intent(in   ) :: umac(:,:)
    real(kind=dp_t), intent(in   ) ::   p0_1(:,0:)
    real(kind=dp_t), intent(in   ) ::   p0_2(:,0:)
    real(kind=dp_t), intent(in   ) :: rho0_1(:,0:)
    real(kind=dp_t), intent(in   ) :: rho0_2(:,0:)
    real(kind=dp_t), intent(in   ) ::      psi(:,0:)
    real(kind=dp_t), intent(in   ) :: dx(:,:)
    logical        , intent(in   ) :: add_thermal
    type(ml_layout), intent(in   ) :: mla
    type(bc_level) , intent(in   ) :: the_bc_level(:)

    ! local
    integer                  :: i,n,comp,ng_f,ng_um,ng_th,ng_p0,ng_pm
    integer                  :: lo(mla%dim),hi(mla%dim),dm,nlevs

    real(kind=dp_t), pointer :: ump(:,:,:,:)
    real(kind=dp_t), pointer :: vmp(:,:,:,:)
    real(kind=dp_t), pointer :: wmp(:,:,:,:)
    real(kind=dp_t), pointer :: fp(:,:,:,:)
    real(kind=dp_t), pointer :: tp(:,:,:,:)
    real(kind=dp_t), pointer :: pp(:,:,:,:)
    real(kind=dp_t), pointer :: pmx(:,:,:,:)
    real(kind=dp_t), pointer :: pmy(:,:,:,:)
    real(kind=dp_t), pointer :: pmz(:,:,:,:)

    type(multifab) :: p0_cart(mla%nlevel)
    type(multifab) :: p0mac(mla%nlevel,mla%dim)

    real(kind=dp_t) ::   p0(nlevs_radial,0:nr_fine-1)
    real(kind=dp_t) :: rho0(nlevs_radial,0:nr_fine-1)
    real(kind=dp_t) :: grav(nlevs_radial,0:nr_fine-1)

    type(bl_prof_timer), save :: bpt

    call build(bpt, "mkrhohforce")

    ! if we are doing the prediction, then it only makes sense to be in
    ! this routine if the quantity we are predicting is rhoh', h, or rhoh
    if (is_prediction .and. &
         .not. (enthalpy_pred_type .eq. predict_rhohprime .or. &
                enthalpy_pred_type .eq. predict_h .or. &
                enthalpy_pred_type .eq. predict_rhoh) ) then
       call bl_error("ERROR: should only call mkrhohforce when predicting rhoh', h, or rhoh")
    endif

    dm = mla%dim
    nlevs = mla%nlevel

    rho0 = HALF*(rho0_1 + rho0_2)
    p0 = HALF * (p0_1 + p0_2)

    if (spherical .eq. 1) then
       do n = 1,nlevs
          call multifab_build(p0_cart(n),mla%la(n),1,1)
          do comp=1,dm
             call multifab_build_edge(p0mac(n,comp),mla%la(n),1,1,comp)
          end do
       end do

       call put_1d_array_on_cart(p0,p0_cart,foextrap_comp,.false.,.false.,&
                                 dx,the_bc_level,mla)

       call make_s0mac(mla,p0,p0mac,dx,foextrap_comp,the_bc_level)
    end if

    ng_f  = nghost(scal_force(1))
    ng_um = nghost(umac(1,1))
    ng_th = nghost(thermal(1))
    ng_pm = nghost(p0mac(1,1))
    ng_p0 = nghost(p0_cart(1))

    call make_grav_cell(grav,rho0)

    do n=1,nlevs

       do i=1, nfabs(scal_force(n))
          fp => dataptr(scal_force(n), i)
          tp  => dataptr(thermal(n),i)
          lo = lwb(get_box(scal_force(n),i))
          hi = upb(get_box(scal_force(n),i))
          select case (dm)
          case (1)
             ump => dataptr(umac(n,1),i)
             call mkrhohforce_1d(n,fp(:,1,1,rhoh_comp), ng_f, is_prediction, &
                                 ump(:,1,1,1), ng_um, &
                                 tp(:,1,1,1), ng_th, lo, hi, p0(n,:), &
                                 rho0(n,:), grav(n,:), psi(n,:), add_thermal)
          case (2)
             vmp => dataptr(umac(n,2),i)
             call mkrhohforce_2d(n,fp(:,:,1,rhoh_comp), ng_f, is_prediction, &
                                 vmp(:,:,1,1), ng_um, &
                                 tp(:,:,1,1), ng_th, lo, hi, p0(n,:),  &
                                 rho0(n,:), grav(n,:), psi(n,:), add_thermal)
          case(3)
             wmp => dataptr(umac(n,3), i)
             if (spherical .eq. 0) then
                call mkrhohforce_3d(n,fp(:,:,:,rhoh_comp), ng_f, is_prediction, &
                                    wmp(:,:,:,1), ng_um, &
                                    tp(:,:,:,1), ng_th, lo, hi, p0(n,:), &
                                    rho0(n,:), grav(n,:), psi(n,:), add_thermal)
             else
                ump => dataptr(umac(n,1),i)
                vmp => dataptr(umac(n,2),i)
                wmp => dataptr(umac(n,3),i)
                pp  => dataptr(p0_cart(n),i)
                pmx => dataptr(p0mac(n,1),i)
                pmy => dataptr(p0mac(n,2),i)
                pmz => dataptr(p0mac(n,3),i)
                call mkrhohforce_3d_sphr(fp(:,:,:,rhoh_comp), ng_f, is_prediction, &
                                         ump(:,:,:,1), vmp(:,:,:,1), wmp(:,:,:,1), ng_um, &
                                         tp(:,:,:,1), ng_th, pp(:,:,:,1), ng_p0, &
                                         pmx(:,:,:,1), pmy(:,:,:,1), pmz(:,:,:,1), ng_pm, &
                                         lo, hi, dx(n,:), psi(1,:), add_thermal)
             end if
          end select
       end do

    end do

    if (spherical .eq. 1) then
       do n=1,nlevs
          call destroy(p0_cart(n))
          do comp=1,dm
             call destroy(p0mac(n,comp))
          end do
       end do
    end if

    ! restrict data and fill all ghost cells
    call ml_restrict_and_fill(nlevs,scal_force,mla%mba%rr,the_bc_level, &
                              icomp=rhoh_comp, &
                              bcomp=foextrap_comp, &
                              nc=1, &
                              ng=scal_force(1)%ng)

    call destroy(bpt)

  end subroutine mkrhohforce

  subroutine mkrhohforce_1d(n,rhoh_force,ng_f,is_prediction, &
                            umac,ng_um,thermal,ng_th,lo,hi, &
                            p0,rho0,grav,psi,add_thermal)

    use geometry, only: dr, nr, base_cutoff_density_coord

    ! compute the source terms for the non-reactive part of the enthalpy equation {w dp0/dr}

    integer,         intent(in   ) :: n,lo(:),hi(:),ng_f,ng_um,ng_th
    logical,         intent(in   ) :: is_prediction
    real(kind=dp_t), intent(  out) :: rhoh_force(lo(1)-ng_f :)
    real(kind=dp_t), intent(in   ) ::       umac(lo(1)-ng_um:)
    real(kind=dp_t), intent(in   ) ::    thermal(lo(1)-ng_th:)
    real(kind=dp_t), intent(in   ) ::   p0(0:)
    real(kind=dp_t), intent(in   ) :: rho0(0:)
    real(kind=dp_t), intent(in   ) :: grav(0:)
    real(kind=dp_t), intent(in   ) :: psi(0:)
    logical        , intent(in   ) :: add_thermal

    real(kind=dp_t) :: gradp0, uadv
    integer :: i

    ! Add wtilde d(p0)/dr
    do i = lo(1),hi(1)

       if (i .lt. base_cutoff_density_coord(n)) then
          gradp0 = rho0(i) * grav(i)
       else if (i.eq.nr(n)-1) then
          ! NOTE: this should be zero since p0 is constant up here
          gradp0 = ( p0(i) - p0(i-1) ) / dr(n)
       else
          ! NOTE: this should be zero since p0 is constant up here
          gradp0 = ( p0(i+1) - p0(i) ) / dr(n)
       end if

       uadv = HALF*(umac(i)+umac(i+1))
       rhoh_force(i) = uadv * gradp0
    enddo

    ! psi should always be in the force if we are doing the final update
    ! For prediction, it should not be in the force if we are predicting
    ! (rho h)', but should be there if we are predicting h or rhoh
    if ((is_prediction .AND. enthalpy_pred_type == predict_h) .OR. &
        (is_prediction .AND. enthalpy_pred_type == predict_rhoh) .OR. &
        (.NOT. is_prediction)) then
       do i = lo(1),hi(1)
          rhoh_force(i) =  rhoh_force(i)  + psi(i)
       end do
    endif

    if (add_thermal) then
       do i=lo(1),hi(1)
          rhoh_force(i) = rhoh_force(i) + thermal(i)
       end do
    end if

  end subroutine mkrhohforce_1d

  subroutine mkrhohforce_2d(n,rhoh_force,ng_f,is_prediction, &
                            vmac,ng_um,thermal,ng_th,lo,hi, &
                            p0,rho0,grav,psi,add_thermal)

    use geometry, only: dr, nr, base_cutoff_density_coord

    ! compute the source terms for the non-reactive part of the enthalpy equation {w dp0/dr}

    integer,         intent(in   ) :: n,lo(:),hi(:),ng_f,ng_um,ng_th
    logical,         intent(in   ) :: is_prediction
    real(kind=dp_t), intent(  out) :: rhoh_force(lo(1)-ng_f :,lo(2)-ng_f :)
    real(kind=dp_t), intent(in   ) ::       vmac(lo(1)-ng_um:,lo(2)-ng_um:)
    real(kind=dp_t), intent(in   ) ::    thermal(lo(1)-ng_th:,lo(2)-ng_th:)
    real(kind=dp_t), intent(in   ) ::   p0(0:)
    real(kind=dp_t), intent(in   ) :: rho0(0:)
    real(kind=dp_t), intent(in   ) :: grav(0:)
    real(kind=dp_t), intent(in   ) :: psi(0:)
    logical        , intent(in   ) :: add_thermal

    real(kind=dp_t) :: gradp0, vadv
    integer :: i,j

    ! Add wtilde d(p0)/dr
    do j = lo(2),hi(2)

       if (j .lt. base_cutoff_density_coord(n)) then
          gradp0 = rho0(j) * grav(j)
       else if (j.eq.nr(n)-1) then
          ! NOTE: this should be zero since p0 is constant up here
          gradp0 = ( p0(j) - p0(j-1) ) / dr(n)
       else
          ! NOTE: this should be zero since p0 is constant up here
          gradp0 = ( p0(j+1) - p0(j) ) / dr(n)
       end if

       do i = lo(1),hi(1)
          vadv = HALF*(vmac(i,j)+vmac(i,j+1))
          rhoh_force(i,j) =  vadv * gradp0
       end do
    enddo

    ! psi should always be in the force if we are doing the final update
    ! For prediction, it should not be in the force if we are predicting
    ! (rho h)', but should be there if we are predicting h or rhoh
    if ((is_prediction .AND. enthalpy_pred_type == predict_h) .OR. &
        (is_prediction .AND. enthalpy_pred_type == predict_rhoh) .OR. &
        (.NOT. is_prediction)) then
       do j = lo(2),hi(2)
          do i = lo(1),hi(1)
             rhoh_force(i,j) =  rhoh_force(i,j)  + psi(j)
          end do
       end do
    endif

    if (add_thermal) then
       do j=lo(2),hi(2)
          do i=lo(1),hi(1)
             rhoh_force(i,j) = rhoh_force(i,j) + thermal(i,j)
          end do
       end do
    end if

  end subroutine mkrhohforce_2d

  subroutine mkrhohforce_3d(n,rhoh_force,ng_f,is_prediction, &
                            wmac,ng_um,thermal,ng_th,lo,hi,&
                            p0,rho0,grav,psi,add_thermal)

    use geometry, only: dr, nr, base_cutoff_density_coord

    ! compute the source terms for the non-reactive part of the enthalpy equation {w dp0/dr}

    integer,         intent(in   ) :: n,lo(:),hi(:),ng_f,ng_um,ng_th
    logical,         intent(in   ) :: is_prediction
    real(kind=dp_t), intent(  out) :: rhoh_force(lo(1)-ng_f :,lo(2)-ng_f :,lo(3)-ng_f :)
    real(kind=dp_t), intent(in   ) ::       wmac(lo(1)-ng_um:,lo(2)-ng_um:,lo(3)-ng_um:)
    real(kind=dp_t), intent(in   ) ::    thermal(lo(1)-ng_th:,lo(2)-ng_th:,lo(3)-ng_th:)
    real(kind=dp_t), intent(in   ) ::   p0(0:)
    real(kind=dp_t), intent(in   ) :: rho0(0:)
    real(kind=dp_t), intent(in   ) :: grav(0:)
    real(kind=dp_t), intent(in   ) :: psi(0:)
    logical        , intent(in   ) :: add_thermal

    real(kind=dp_t) :: gradp0,wadv
    integer :: i,j,k
    !
    ! Add wtilde d(p0)/dr
    !
    !$OMP PARALLEL DO PRIVATE(i,j,k,gradp0,wadv)
    do k = lo(3),hi(3)

       if (k .lt. base_cutoff_density_coord(n)) then
          gradp0 = rho0(k) * grav(k)
       else if (k.eq.nr(n)-1) then
          ! NOTE: this should be zero since p0 is constant up here
          gradp0 = ( p0(k) - p0(k-1) ) / dr(n)
       else
          ! NOTE: this should be zero since p0 is constant up here
          gradp0 = ( p0(k+1) - p0(k) ) / dr(n)
       end if

       do j = lo(2),hi(2)
          do i = lo(1),hi(1)
             wadv = HALF*(wmac(i,j,k)+wmac(i,j,k+1))
             rhoh_force(i,j,k) = wadv * gradp0
          end do
       end do
    enddo
    !$OMP END PARALLEL DO

    ! psi should always be in the force if we are doing the final update
    ! For prediction, it should not be in the force if we are predicting
    ! (rho h)', but should be there if we are predicting h or rhoh
    if ((is_prediction .AND. enthalpy_pred_type == predict_h) .OR. &
        (is_prediction .AND. enthalpy_pred_type == predict_rhoh) .OR. &
        (.NOT. is_prediction)) then
       !$OMP PARALLEL DO PRIVATE(i,j,k)
       do k = lo(3),hi(3)
          do j = lo(2),hi(2)
             do i = lo(1),hi(1)
                rhoh_force(i,j,k) = rhoh_force(i,j,k) + psi(k)
             end do
          end do
       enddo
       !$OMP END PARALLEL DO
    endif

    if (add_thermal) then
       !$OMP PARALLEL DO PRIVATE(i,j,k)
       do k=lo(3),hi(3)
          do j=lo(2),hi(2)
             do i=lo(1),hi(1)
                rhoh_force(i,j,k) = rhoh_force(i,j,k) + thermal(i,j,k)
             end do
          end do
       end do
       !$OMP END PARALLEL DO
    end if

  end subroutine mkrhohforce_3d

  subroutine mkrhohforce_3d_sphr(rhoh_force,ng_f,is_prediction, &
                                 umac,vmac,wmac,ng_um,thermal,ng_th, &
                                 p0_cart,ng_p0, &
                                 p0macx,p0macy,p0macz, &
                                 ng_pm,lo,hi,dx,psi,add_thermal)

    use fill_3d_module

    ! compute the source terms for the non-reactive part of the enthalpy equation {w dp0/dr}

    integer,         intent(in   ) :: lo(:),hi(:),ng_f,ng_um,ng_th,ng_p0,ng_pm
    logical,         intent(in   ) :: is_prediction
    real(kind=dp_t), intent(  out) :: rhoh_force(lo(1)-ng_f :,lo(2)-ng_f :,lo(3)-ng_f :)
    real(kind=dp_t), intent(in   ) ::       umac(lo(1)-ng_um:,lo(2)-ng_um:,lo(3)-ng_um:)
    real(kind=dp_t), intent(in   ) ::       vmac(lo(1)-ng_um:,lo(2)-ng_um:,lo(3)-ng_um:)
    real(kind=dp_t), intent(in   ) ::       wmac(lo(1)-ng_um:,lo(2)-ng_um:,lo(3)-ng_um:)
    real(kind=dp_t), intent(in   ) ::    thermal(lo(1)-ng_th:,lo(2)-ng_th:,lo(3)-ng_th:)
    real(kind=dp_t), intent(in   ) ::    p0_cart(lo(1)-ng_p0:,lo(2)-ng_p0:,lo(3)-ng_p0:)
    real(kind=dp_t), intent(in   ) ::     p0macx(lo(1)-ng_pm:,lo(2)-ng_pm:,lo(3)-ng_pm:)
    real(kind=dp_t), intent(in   ) ::     p0macy(lo(1)-ng_pm:,lo(2)-ng_pm:,lo(3)-ng_pm:)
    real(kind=dp_t), intent(in   ) ::     p0macz(lo(1)-ng_pm:,lo(2)-ng_pm:,lo(3)-ng_pm:)
    real(kind=dp_t), intent(in   ) :: dx(:)
    real(kind=dp_t), intent(in   ) :: psi(0:)
    logical        , intent(in   ) :: add_thermal

    real(kind=dp_t), allocatable :: psi_cart(:,:,:,:)

    real(kind=dp_t) :: divup, p0divu
    integer         :: i,j,k
    !
    ! Here we make u grad p = div (u p) - p div (u)
    !
    !$OMP PARALLEL DO PRIVATE(i,j,k,divup,p0divu)
    do k = lo(3),hi(3)
       do j = lo(2),hi(2)
          do i = lo(1),hi(1)

             divup = (umac(i+1,j,k) * p0macx(i+1,j,k) - umac(i,j,k) * p0macx(i,j,k)) / dx(1) + &
                     (vmac(i,j+1,k) * p0macy(i,j+1,k) - vmac(i,j,k) * p0macy(i,j,k)) / dx(2) + &
                     (wmac(i,j,k+1) * p0macz(i,j,k+1) - wmac(i,j,k) * p0macz(i,j,k)) / dx(3)

             p0divu = ( (umac(i+1,j,k) - umac(i,j,k)) / dx(1) + &
                        (vmac(i,j+1,k) - vmac(i,j,k)) / dx(2) + &
                        (wmac(i,j,k+1) - wmac(i,j,k)) / dx(3) ) * p0_cart(i,j,k)

             rhoh_force(i,j,k) = divup - p0divu

          end do
       end do
    end do
    !$OMP END PARALLEL DO
    !
    ! psi should always be in the force if we are doing the final update
    ! For prediction, it should not be in the force if we are predicting
    ! (rho h)', but should be there if we are predicting h or rhoh
    !
    if ((is_prediction .AND. enthalpy_pred_type == predict_h) .OR. &
        (is_prediction .AND. enthalpy_pred_type == predict_rhoh) .OR. &
        (.NOT. is_prediction)) then

       allocate(psi_cart(lo(1):hi(1),lo(2):hi(2),lo(3):hi(3),1))

       call put_1d_array_on_cart_3d_sphr(.false.,.false.,psi,psi_cart,lo,hi,dx,0)

       !$OMP PARALLEL DO PRIVATE(i,j,k)
       do k = lo(3),hi(3)
          do j = lo(2),hi(2)
             do i = lo(1),hi(1)
                rhoh_force(i,j,k) = rhoh_force(i,j,k) + psi_cart(i,j,k,1)
             enddo
          enddo
       enddo
       !$OMP END PARALLEL DO
       deallocate(psi_cart)
    endif

    if (add_thermal) then
       !$OMP PARALLEL DO PRIVATE(i,j,k)
       do k=lo(3),hi(3)
          do j=lo(2),hi(2)
             do i=lo(1),hi(1)
                rhoh_force(i,j,k) = rhoh_force(i,j,k) + thermal(i,j,k)
             end do
          end do
       end do
       !$OMP END PARALLEL DO
    end if

  end subroutine mkrhohforce_3d_sphr

  subroutine mkhprimeforce(mla, sold, scal_force, is_prediction, thermal, umac, &
                           p0_old, h0_old, psi, dx, add_thermal, the_bc_level)

    use bl_prof_module
    use variables, only: foextrap_comp, rhoh_comp, rho_comp
    use geometry, only: spherical
    use fill_3d_module, only: put_1d_array_on_cart

    type(multifab) , intent(inout) :: scal_force(:)
    type(multifab) , intent(in   ) :: sold(:)
    logical        , intent(in   ) :: is_prediction
    type(multifab) , intent(in   ) :: thermal(:)
    type(multifab) , intent(in   ) :: umac(:,:)
    real(kind=dp_t), intent(in   ) :: p0_old(:,0:)
    real(kind=dp_t), intent(in   ) :: h0_old(:,0:)
    real(kind=dp_t), intent(in   ) :: psi(:,0:)
    real(kind=dp_t), intent(in   ) :: dx(:,:)
    logical        , intent(in   ) :: add_thermal
    type(ml_layout), intent(in   ) :: mla
    type(bc_level) , intent(in   ) :: the_bc_level(:)

    ! local
    integer                  :: i,n
    integer                  :: lo(mla%dim),hi(mla%dim),dm,nlevs
    integer                  :: ng_f,ng_um,ng_th,ng_s
    real(kind=dp_t), pointer :: ump(:,:,:,:)
    real(kind=dp_t), pointer :: vmp(:,:,:,:)
    real(kind=dp_t), pointer :: wmp(:,:,:,:)
    real(kind=dp_t), pointer :: fp(:,:,:,:)
    real(kind=dp_t), pointer :: tp(:,:,:,:)
    real(kind=dp_t), pointer :: pp(:,:,:,:)
    real(kind=dp_t), pointer :: hp(:,:,:,:)
    real(kind=dp_t), pointer :: sop(:,:,:,:)

    type(multifab)  :: p0_cart(mla%nlevel)
    type(multifab)  :: h0_cart(mla%nlevel)

    type(bl_prof_timer), save :: bpt

    call build(bpt, "mkhprimeforce")

    ! if we are doing the prediction, then it only makes sense to be in
    ! this routine if the quantity we are predicting is h'
    if (is_prediction .and. enthalpy_pred_type .ne. predict_hprime) then
       call bl_error("ERROR: should only call mkhprimeforce when predicting h'")
    endif

    dm = mla%dim
    nlevs = mla%nlevel

    ng_f  = nghost(scal_force(1))
    ng_um = nghost(umac(1,1))
    ng_th = nghost(thermal(1))
    ng_s  = nghost(sold(1))

    if (spherical .eq. 1) then
       do n = 1,nlevs
          call multifab_build(p0_cart(n),mla%la(n),1,1)
          call multifab_build(h0_cart(n),mla%la(n),1,1)
       end do
    end if

    do n=1,nlevs

       if (spherical .eq. 1) then
          call put_1d_array_on_cart(p0_old,p0_cart,foextrap_comp,.false.,.false.,&
                                    dx,the_bc_level,mla)
          call put_1d_array_on_cart(h0_old,h0_cart,foextrap_comp,.false.,.false.,&
                                    dx,the_bc_level,mla)
       end if

       do i=1, nfabs(scal_force(n))
          fp => dataptr(scal_force(n), i)
          ump => dataptr(umac(n,1),i)
          vmp => dataptr(umac(n,2),i)
          tp  => dataptr(thermal(n),i)
          sop => dataptr(sold(n), i)
          lo = lwb(get_box(scal_force(n),i))
          hi = upb(get_box(scal_force(n),i))
          select case (dm)
          case (2)
             call bl_error("mkscalforce: mkhprime force not written for 2d")
          case(3)
             wmp  => dataptr(umac(n,3), i)
             if (spherical .eq. 0) then
             call bl_error("mkscalforce: mkhprime force not written for 3d plane parallel")
             else
                pp  => dataptr(p0_cart(n),i)
                hp  => dataptr(h0_cart(n),i)
                call mkhprimeforce_3d_sphr(fp(:,:,:,rhoh_comp), ng_f, is_prediction, &
                                           ump(:,:,:,1), vmp(:,:,:,1), wmp(:,:,:,1), ng_um, &
                                           tp(:,:,:,1), ng_th, &
                                           sop(:,:,:,rho_comp), ng_s, &
                                           pp(:,:,:,1), hp(:,:,:,1), &
                                           lo, hi, dx(n,:), &
                                           psi(1,:), add_thermal)
             end if
          end select
       end do

    end do

    if (spherical .eq. 1) then
       do n = 1,nlevs
          call destroy(p0_cart(n))
          call destroy(h0_cart(n))
       end do
    end if

    ! restrict data and fill all ghost cells
    call ml_restrict_and_fill(nlevs,scal_force,mla%mba%rr,the_bc_level, &
                              icomp=rhoh_comp, &
                              bcomp=foextrap_comp, &
                              nc=1, &
                              ng=scal_force(1)%ng)

    call destroy(bpt)

  end subroutine mkhprimeforce

  subroutine mkhprimeforce_3d_sphr(rhoh_force,ng_f,is_prediction,umac,vmac,wmac,ng_um, &
                                   thermal,ng_th,rhoold,ng_s,p0_cart,h0_cart,lo,hi, &
                                   dx,psi,add_thermal)

    use fill_3d_module

    ! compute the source terms for the non-reactive part of the enthalpy equation {w dp0/dr}

    integer,         intent(in   ) :: lo(:),hi(:),ng_f,ng_um,ng_th,ng_s
    logical,         intent(in   ) :: is_prediction
    real(kind=dp_t), intent(  out) :: rhoh_force(lo(1)-ng_f :,lo(2)-ng_f :,lo(3)-ng_f :)
    real(kind=dp_t), intent(in   ) ::       umac(lo(1)-ng_um:,lo(2)-ng_um:,lo(3)-ng_um:)
    real(kind=dp_t), intent(in   ) ::       vmac(lo(1)-ng_um:,lo(2)-ng_um:,lo(3)-ng_um:)
    real(kind=dp_t), intent(in   ) ::       wmac(lo(1)-ng_um:,lo(2)-ng_um:,lo(3)-ng_um:)
    real(kind=dp_t), intent(in   ) ::    thermal(lo(1)-ng_th:,lo(2)-ng_th:,lo(3)-ng_th:)
    real(kind=dp_t), intent(in   ) ::     rhoold(lo(1)-ng_s :,lo(2)-ng_s :,lo(3)-ng_s :)
    real(kind=dp_t), intent(in   ) ::    p0_cart(lo(1)-   1 :,lo(2)-   1 :,lo(3)-   1 :)
    real(kind=dp_t), intent(in   ) ::    h0_cart(lo(1)-   1 :,lo(2)-   1 :,lo(3)-   1 :)
    real(kind=dp_t), intent(in   ) :: dx(:)
    real(kind=dp_t), intent(in   ) :: psi(0:)
    logical        , intent(in   ) :: add_thermal

    real(kind=dp_t), allocatable :: psi_cart(:,:,:,:)

    real(kind=dp_t) :: p0_lox,p0_hix,p0_loy,p0_hiy,p0_loz,p0_hiz
    real(kind=dp_t) :: h0_lox,h0_hix,h0_loy,h0_hiy,h0_loz,h0_hiz
    real(kind=dp_t) :: divup, p0divu
    real(kind=dp_t) :: divuh, h0divu
    integer         :: i,j,k
    !
    ! Here we make u grad p = div (u p) - p div (u)
    !
    !$OMP PARALLEL DO PRIVATE(i,j,k,p0_lox,p0_hix,p0_loy,p0_hiy,p0_loz,p0_hiz,divup,p0divu)
    do k = lo(3),hi(3)
       do j = lo(2),hi(2)
          do i = lo(1),hi(1)

             p0_lox = HALF * (p0_cart(i,j,k) + p0_cart(i-1,j,k))
             p0_hix = HALF * (p0_cart(i,j,k) + p0_cart(i+1,j,k))
             p0_loy = HALF * (p0_cart(i,j,k) + p0_cart(i,j-1,k))
             p0_hiy = HALF * (p0_cart(i,j,k) + p0_cart(i,j+1,k))
             p0_loz = HALF * (p0_cart(i,j,k) + p0_cart(i,j,k-1))
             p0_hiz = HALF * (p0_cart(i,j,k) + p0_cart(i,j,k+1))

             divup = (umac(i+1,j,k) * p0_hix - umac(i,j,k) * p0_lox) / dx(1) + &
                     (vmac(i,j+1,k) * p0_hiy - vmac(i,j,k) * p0_loy) / dx(2) + &
                     (wmac(i,j,k+1) * p0_hiz - wmac(i,j,k) * p0_loz) / dx(3)

             p0divu = ( (umac(i+1,j,k) - umac(i,j,k)) / dx(1) + &
                        (vmac(i,j+1,k) - vmac(i,j,k)) / dx(2) + &
                        (wmac(i,j,k+1) - wmac(i,j,k)) / dx(3) ) * p0_cart(i,j,k)

             rhoh_force(i,j,k) = (divup - p0divu) / rhoold(i,j,k)

          end do
       end do
    end do
    !$OMP END PARALLEL DO
    !
    ! Here we make u grad h_0 = div (u h_0) - h_0 div (u)
    !
    !$OMP PARALLEL DO PRIVATE(i,j,k,h0_lox,h0_hix,h0_loy,h0_hiy,h0_loz,h0_hiz,divuh,h0divu)
    do k = lo(3),hi(3)
       do j = lo(2),hi(2)
          do i = lo(1),hi(1)

             h0_lox = HALF * (h0_cart(i,j,k) + h0_cart(i-1,j,k))
             h0_hix = HALF * (h0_cart(i,j,k) + h0_cart(i+1,j,k))
             h0_loy = HALF * (h0_cart(i,j,k) + h0_cart(i,j-1,k))
             h0_hiy = HALF * (h0_cart(i,j,k) + h0_cart(i,j+1,k))
             h0_loz = HALF * (h0_cart(i,j,k) + h0_cart(i,j,k-1))
             h0_hiz = HALF * (h0_cart(i,j,k) + h0_cart(i,j,k+1))

             divuh = (umac(i+1,j,k) * h0_hix - umac(i,j,k) * h0_lox) / dx(1) + &
                     (vmac(i,j+1,k) * h0_hiy - vmac(i,j,k) * h0_loy) / dx(2) + &
                     (wmac(i,j,k+1) * h0_hiz - wmac(i,j,k) * h0_loz) / dx(3)

             h0divu = ( (umac(i+1,j,k) - umac(i,j,k)) / dx(1) + &
                        (vmac(i,j+1,k) - vmac(i,j,k)) / dx(2) + &
                        (wmac(i,j,k+1) - wmac(i,j,k)) / dx(3) ) * h0_cart(i,j,k)

             rhoh_force(i,j,k) = rhoh_force(i,j,k) - divuh + h0divu

          end do
       end do
    end do

    ! psi should always be in the force if we are doing the final update
    ! For prediction, it should not be in the force if we are predicting
    ! (rho h)', but should be there if we are predicting h or rhoh
    if (.NOT. is_prediction) then

       allocate(psi_cart(lo(1):hi(1),lo(2):hi(2),lo(3):hi(3),1))
       call put_1d_array_on_cart_3d_sphr(.false.,.false.,psi,psi_cart,lo,hi,dx,0)

       !$OMP PARALLEL DO PRIVATE(i,j,k)
       do k = lo(3),hi(3)
          do j = lo(2),hi(2)
             do i = lo(1),hi(1)
                rhoh_force(i,j,k) = rhoh_force(i,j,k) + psi_cart(i,j,k,1)
             enddo
          enddo
       enddo
       !$OMP END PARALLEL DO
       deallocate(psi_cart)
    endif

    if (add_thermal) then
       !$OMP PARALLEL DO PRIVATE(i,j,k)
       do k=lo(3),hi(3)
          do j=lo(2),hi(2)
             do i=lo(1),hi(1)
                rhoh_force(i,j,k) = rhoh_force(i,j,k) + thermal(i,j,k)
             end do
          end do
       end do
       !$OMP END PARALLEL DO
    end if

  end subroutine mkhprimeforce_3d_sphr

  subroutine mktempforce(mla,temp_force,umac,s,thermal,p0_old,psi,dx,the_bc_level)

    use bl_prof_module
    use variables, only: foextrap_comp, temp_comp
    use geometry, only: spherical
    use fill_3d_module, only: put_1d_array_on_cart

    type(ml_layout), intent(in   ) :: mla
    type(multifab) , intent(inout) :: temp_force(:)
    type(multifab) , intent(in   ) :: umac(:,:)
    type(multifab) , intent(in   ) :: s(:)
    type(multifab) , intent(in   ) :: thermal(:)
    real(kind=dp_t), intent(in   ) :: p0_old(:,0:)
    real(kind=dp_t), intent(in   ) :: psi(:,0:)
    real(kind=dp_t), intent(in   ) :: dx(:,:)
    type(bc_level) , intent(in   ) :: the_bc_level(:)

    ! local
    integer         :: i,n,ng_f,ng_um,ng_s,ng_th
    type(multifab)  :: p0_cart(mla%nlevel)
    integer         :: lo(mla%dim),hi(mla%dim),dm,nlevs

    real(kind=dp_t), pointer :: tp(:,:,:,:)
    real(kind=dp_t), pointer :: ump(:,:,:,:)
    real(kind=dp_t), pointer :: vmp(:,:,:,:)
    real(kind=dp_t), pointer :: wmp(:,:,:,:)
    real(kind=dp_t), pointer :: sp(:,:,:,:)
    real(kind=dp_t), pointer :: fp(:,:,:,:)
    real(kind=dp_t), pointer :: pp(:,:,:,:)

    type(bl_prof_timer), save :: bpt

    call build(bpt, "mktempforce")

    ! it only makes sense to be in
    ! this routine if the quantity we are predicting is T or T'
    if (.not. (enthalpy_pred_type .eq. predict_T_then_rhohprime .or. &
               enthalpy_pred_type .eq. predict_T_then_h .or. &
               enthalpy_pred_type .eq. predict_Tprime_then_h) ) then
       call bl_error("ERROR: should only call mkrhohforce when predicting T or T'")
    endif

    dm = mla%dim
    nlevs = mla%nlevel

    ng_f  = nghost(temp_force(1))
    ng_um = nghost(umac(1,1))
    ng_s  = nghost(s(1))
    ng_th = nghost(thermal(1))

    if (spherical .eq. 1) then
       do n = 1,nlevs
          call multifab_build(p0_cart(n),mla%la(n),1,1)
       end do
       call put_1d_array_on_cart(p0_old,p0_cart,foextrap_comp,.false.,.false.,&
                                 dx,the_bc_level,mla)
    end if

    do n=1,nlevs

       do i=1, nfabs(temp_force(n))
          fp  => dataptr(temp_force(n),i)
          ump => dataptr(umac(n,1),i)
          sp  => dataptr(s(n),i)
          tp  => dataptr(thermal(n),i)

          lo  =  lwb(get_box(s(n),i))
          hi  =  upb(get_box(s(n),i))

          select case (dm)
          case (1)
             call mktempforce_1d(n, fp(:,1,1,temp_comp), ng_f, sp(:,1,1,:), ng_s, &
                                 ump(:,1,1,1), ng_um, tp(:,1,1,1), ng_th, lo, hi, &
                                 p0_old(n,:), psi(n,:))

          case (2)
             vmp => dataptr(umac(n,2),i)
             call mktempforce_2d(n, fp(:,:,1,temp_comp), ng_f, sp(:,:,1,:), ng_s, &
                                 vmp(:,:,1,1), ng_um, tp(:,:,1,1), ng_th, lo, hi, &
                                 p0_old(n,:), psi(n,:))

          case(3)
             vmp => dataptr(umac(n,2),i)
             wmp => dataptr(umac(n,3),i)
             if (spherical .eq. 1) then
                pp  => dataptr(p0_cart(n),i)
                call mktempforce_3d_sphr(fp(:,:,:,temp_comp), ng_f, sp(:,:,:,:), ng_s, &
                                         ump(:,:,:,1), vmp(:,:,:,1), wmp(:,:,:,1), ng_um, &
                                         tp(:,:,:,1), ng_th, lo, hi, &
                                         pp(:,:,:,1), psi(1,:), dx(n,:))
             else
                call mktempforce_3d(n, fp(:,:,:,temp_comp), ng_f, sp(:,:,:,:), ng_s, &
                                    wmp(:,:,:,1), ng_um, tp(:,:,:,1), ng_th, lo, hi, &
                                    p0_old(n,:), psi(n,:))
             end if
          end select
       end do

    end do

    if (spherical .eq. 1) then
       do n = 1,nlevs
          call destroy(p0_cart(n))
       end do
    end if

    ! restrict data and fill all ghost cells
    call ml_restrict_and_fill(nlevs,temp_force,mla%mba%rr,the_bc_level, &
                              icomp=temp_comp, &
                              bcomp=foextrap_comp, &
                              nc=1, &
                              ng=temp_force(1)%ng)

    call destroy(bpt)

  end subroutine mktempforce

  subroutine mktempforce_1d(n, temp_force, ng_f, s, ng_s, wmac, ng_um, thermal, ng_th, &
                            lo, hi, p0_old, psi)

    use geometry, only: dr, nr
    use variables, only: temp_comp, rho_comp, spec_comp
    use eos_module, only: eos_input_rt, eos
    use eos_type_module
    use network, only: nspec

    ! compute the source terms for temperature

    integer,         intent(in   ) :: lo(:),hi(:),ng_f,ng_s,ng_um,ng_th
    integer,         intent(in   ) :: n
    real(kind=dp_t), intent(  out) :: temp_force(lo(1)-ng_f :)
    real(kind=dp_t), intent(in   ) ::          s(lo(1)-ng_s :,:)
    real(kind=dp_t), intent(in   ) ::       wmac(lo(1)-ng_um:)
    real(kind=dp_t), intent(in   ) ::    thermal(lo(1)-ng_th:)
    real(kind=dp_t), intent(in   ) :: p0_old(0:), psi(0:)

    integer :: i

    real(kind=dp_t) :: gradp0, wadv, dhdp

    type (eos_t) :: eos_state
    integer :: pt_index(MAX_SPACEDIM)

    do i = lo(1),hi(1)

       if (i.eq.0) then
          gradp0 = ( p0_old(i+1) - p0_old(i) ) / dr(n)
       else if (i.eq.nr(n)-1) then
          gradp0 = ( p0_old(i) - p0_old(i-1) ) / dr(n)
       else
          gradp0 = HALF * ( p0_old(i+1) - p0_old(i-1) ) / dr(n)
       end if

       eos_state%T     = s(i,temp_comp)
       eos_state%rho   = s(i,rho_comp)
       eos_state%xn(:) = s(i,spec_comp:spec_comp+nspec-1) / s(i,rho_comp)

       pt_index(:) = (/i, -1, -1/)

       ! dens, temp, xmass inputs
       call eos(eos_input_rt, eos_state, pt_index)

       dhdp = ONE / s(i,rho_comp) + ( s(i,rho_comp) * eos_state%dedr - &
                                      eos_state%p / s(i,rho_comp) ) &
                                     / ( s(i,rho_comp) * eos_state%dpdr )

       wadv = HALF*(wmac(i)+wmac(i+1))

       temp_force(i) =  thermal(i) + (ONE - s(i,rho_comp) * dhdp) * &
                          (wadv * gradp0 + psi(i))
       temp_force(i) = temp_force(i) / (eos_state%cp * s(i,rho_comp))

    end do

  end subroutine mktempforce_1d

  subroutine mktempforce_2d(n, temp_force, ng_f, s, ng_s, wmac, ng_um, thermal, ng_th, &
                            lo, hi, p0_old, psi)

    use geometry, only: dr, nr
    use variables, only: temp_comp, rho_comp, spec_comp
    use eos_module, only: eos_input_rt, eos
    use eos_type_module
    use network, only: nspec

    ! compute the source terms for temperature

    integer,         intent(in   ) :: lo(:),hi(:),ng_f,ng_s,ng_um,ng_th
    integer,         intent(in   ) :: n
    real(kind=dp_t), intent(  out) :: temp_force(lo(1)-ng_f :,lo(2)-ng_f :)
    real(kind=dp_t), intent(in   ) ::          s(lo(1)-ng_s :,lo(2)-ng_s :,:)
    real(kind=dp_t), intent(in   ) ::       wmac(lo(1)-ng_um:,lo(2)-ng_um:)
    real(kind=dp_t), intent(in   ) ::    thermal(lo(1)-ng_th:,lo(2)-ng_th:)
    real(kind=dp_t), intent(in   ) :: p0_old(0:), psi(0:)

    integer :: i,j

    real(kind=dp_t) :: gradp0, wadv, dhdp

    type (eos_t) :: eos_state
    integer :: pt_index(MAX_SPACEDIM)

    do j = lo(2),hi(2)

       if (j.eq.0) then
          gradp0 = ( p0_old(j+1) - p0_old(j) ) / dr(n)
       else if (j.eq.nr(n)-1) then
          gradp0 = ( p0_old(j) - p0_old(j-1) ) / dr(n)
       else
          gradp0 = HALF * ( p0_old(j+1) - p0_old(j-1) ) / dr(n)
       end if

       do i = lo(1),hi(1)

          eos_state%T     = s(i,j,temp_comp)
          eos_state%rho   = s(i,j,rho_comp)
          eos_state%xn(:) = s(i,j,spec_comp:spec_comp+nspec-1) / s(i,j,rho_comp)

          pt_index(:) = (/i, j, -1/)

          ! dens, temp, xmass inputs
         call eos(eos_input_rt, eos_state, pt_index)

         dhdp = ONE / s(i,j,rho_comp) + ( s(i,j,rho_comp) * eos_state%dedr - &
                                          eos_state%p / s(i,j,rho_comp) ) &
                                        / ( s(i,j,rho_comp) * eos_state%dpdr )

         wadv = HALF*(wmac(i,j)+wmac(i,j+1))

         temp_force(i,j) =  thermal(i,j) + (ONE - s(i,j,rho_comp) * dhdp) * &
                            (wadv * gradp0 + psi(j))
         temp_force(i,j) = temp_force(i,j) / (eos_state%cp * s(i,j,rho_comp))

       end do
    end do

  end subroutine mktempforce_2d

  subroutine mktempforce_3d(n, temp_force, ng_f, s, ng_s, wmac, ng_um, thermal, ng_th, &
                            lo, hi, p0_old, psi)

    use geometry,  only: dr, nr
    use variables, only: temp_comp, rho_comp, spec_comp
    use eos_module, only: eos_input_rt, eos
    use eos_type_module
    use network, only: nspec

    ! compute the source terms for temperature

    integer,         intent(in   ) :: lo(:),hi(:),n,ng_f,ng_s,ng_um,ng_th
    real(kind=dp_t), intent(  out) :: temp_force(lo(1)-ng_f :,lo(2)-ng_f :,lo(3)-ng_f :)
    real(kind=dp_t), intent(in   ) ::          s(lo(1)-ng_s :,lo(2)-ng_s :,lo(3)-ng_s :,:)
    real(kind=dp_t), intent(in   ) ::       wmac(lo(1)-ng_um:,lo(2)-ng_um:,lo(3)-ng_um:)
    real(kind=dp_t), intent(in   ) ::    thermal(lo(1)-ng_th:,lo(2)-ng_th:,lo(3)-ng_th:)
    real(kind=dp_t), intent(in   ) :: p0_old(0:), psi(0:)

    integer         :: i,j,k
    real(kind=dp_t) :: dhdp, gradp0, wadv

    type (eos_t) :: eos_state
    integer :: pt_index(MAX_SPACEDIM)

    !$OMP PARALLEL DO PRIVATE(i,j,k,gradp0,eos_state,pt_index,dhdp,wadv)
    do k = lo(3),hi(3)
       if (k.eq.0) then
          gradp0 = ( p0_old(k+1) - p0_old(k) ) / dr(n)
       else if (k.eq.nr(n)-1) then
          gradp0 = ( p0_old(k) - p0_old(k-1) ) / dr(n)
       else
          gradp0 = HALF * ( p0_old(k+1) - p0_old(k-1) ) / dr(n)
       end if

       do j = lo(2),hi(2)
          do i = lo(1),hi(1)

             eos_state%T     = s(i,j,k,temp_comp)
             eos_state%rho   = s(i,j,k,rho_comp)
             eos_state%xn(:) = s(i,j,k,spec_comp:spec_comp+nspec-1) / s(i,j,k,rho_comp)

             pt_index(:) = (/i, j, k/)

             ! dens, temp, xmass inputs
             call eos(eos_input_rt, eos_state, pt_index)

             dhdp = ONE / s(i,j,k,rho_comp) + ( s(i,j,k,rho_comp) * eos_state%dedr - &
                  eos_state%p / s(i,j,k,rho_comp) ) / ( s(i,j,k,rho_comp) * eos_state%dpdr )

             wadv = HALF * (wmac(i,j,k+1) + wmac(i,j,k))

             temp_force(i,j,k) =  thermal(i,j,k) + &
                  (ONE - s(i,j,k,rho_comp) * dhdp) * (wadv * gradp0 + psi(k))

             temp_force(i,j,k) = temp_force(i,j,k) / (eos_state%cp * s(i,j,k,rho_comp))

          end do
       end do
    end do
    !$OMP END PARALLEL DO

  end subroutine mktempforce_3d

  subroutine mktempforce_3d_sphr(temp_force, ng_f, s, ng_s, umac, vmac, wmac, ng_um, &
                                 thermal, ng_th, lo, hi, p0_cart, psi, dx)

    use fill_3d_module
    use variables, only: temp_comp, rho_comp, spec_comp
    use eos_module, only: eos_input_rt, eos
    use eos_type_module
    use network, only: nspec

    ! compute the source terms for temperature

    integer,         intent(in   ) :: lo(:),hi(:),ng_f,ng_s,ng_um,ng_th
    real(kind=dp_t), intent(  out) :: temp_force(lo(1)-ng_f :,lo(2)-ng_f :,lo(3)-ng_f :)
    real(kind=dp_t), intent(in   ) ::          s(lo(1)-ng_s :,lo(2)-ng_s :,lo(3)-ng_s :,:)
    real(kind=dp_t), intent(in   ) ::       umac(lo(1)-ng_um:,lo(2)-ng_um:,lo(3)-ng_um:)
    real(kind=dp_t), intent(in   ) ::       vmac(lo(1)-ng_um:,lo(2)-ng_um:,lo(3)-ng_um:)
    real(kind=dp_t), intent(in   ) ::       wmac(lo(1)-ng_um:,lo(2)-ng_um:,lo(3)-ng_um:)
    real(kind=dp_t), intent(in   ) ::    thermal(lo(1)-ng_th:,lo(2)-ng_th:,lo(3)-ng_th:)
    real(kind=dp_t), intent(in   ) ::    p0_cart(lo(1)-1    :,lo(2)-1    :,lo(3)-1    :)
    real(kind=dp_t), intent(in   ) :: psi(0:)
    real(kind=dp_t), intent(in   ) :: dx(:)

    integer :: i,j,k
    real(kind=dp_t) :: p0_lox,p0_hix,p0_loy,p0_hiy,p0_loz,p0_hiz
    real(kind=dp_t) :: divup,p0divu,ugradp,dhdp
    real(kind=dp_t), allocatable :: psi_cart(:,:,:,:)

    type (eos_t) :: eos_state
    integer :: pt_index(MAX_SPACEDIM)

    allocate(psi_cart(lo(1):hi(1),lo(2):hi(2),lo(3):hi(3),1))

    call put_1d_array_on_cart_3d_sphr(.false.,.false.,psi,psi_cart,lo,hi,dx,0)

    !$OMP PARALLEL DO PRIVATE(i,j,k,dhdp,p0_lox,p0_hix,p0_loy,p0_hiy,p0_loz,p0_hiz,divup,p0divu,ugradp,eos_state,pt_index)
    do k = lo(3),hi(3)
       do j = lo(2),hi(2)
          do i = lo(1),hi(1)

             eos_state%T     = s(i,j,k,temp_comp)
             eos_state%rho   = s(i,j,k,rho_comp)
             eos_state%xn(:) = s(i,j,k,spec_comp:spec_comp+nspec-1) / s(i,j,k,rho_comp)

             pt_index(:) = (/i, j, k/)

             ! dens, temp, xmass inputs
             call eos(eos_input_rt, eos_state, pt_index)

             dhdp = ONE / s(i,j,k,rho_comp) + ( s(i,j,k,rho_comp) * eos_state%dedr - &
                                                eos_state%p / s(i,j,k,rho_comp) ) &
                                                / ( s(i,j,k,rho_comp) * eos_state%dpdr )

             p0_lox = HALF * (p0_cart(i,j,k) + p0_cart(i-1,j,k))
             p0_hix = HALF * (p0_cart(i,j,k) + p0_cart(i+1,j,k))
             p0_loy = HALF * (p0_cart(i,j,k) + p0_cart(i,j-1,k))
             p0_hiy = HALF * (p0_cart(i,j,k) + p0_cart(i,j+1,k))
             p0_loz = HALF * (p0_cart(i,j,k) + p0_cart(i,j,k-1))
             p0_hiz = HALF * (p0_cart(i,j,k) + p0_cart(i,j,k+1))

             divup = (umac(i+1,j,k) * p0_hix - umac(i,j,k) * p0_lox) / dx(1) + &
                     (vmac(i,j+1,k) * p0_hiy - vmac(i,j,k) * p0_loy) / dx(2) + &
                     (wmac(i,j,k+1) * p0_hiz - wmac(i,j,k) * p0_loz) / dx(3)

             p0divu = ( (umac(i+1,j,k) - umac(i,j,k)) / dx(1) + &
                        (vmac(i,j+1,k) - vmac(i,j,k)) / dx(2) + &
                        (wmac(i,j,k+1) - wmac(i,j,k)) / dx(3) ) * p0_cart(i,j,k)

             ugradp = divup - p0divu

             temp_force(i,j,k) =  thermal(i,j,k) + &
                  (ONE - s(i,j,k,rho_comp) * dhdp) * (ugradp + psi_cart(i,j,k,1))

             temp_force(i,j,k) = temp_force(i,j,k) / (eos_state%cp * s(i,j,k,rho_comp))

          end do
       end do
    end do
    !$OMP END PARALLEL DO

    deallocate(psi_cart)

  end subroutine mktempforce_3d_sphr

end module mkscalforce_module

max_stars_repo_name: sailoridy/MAESTRO

max_issues_repo_name: sailoridy/MAESTRO

max_forks_repo_name: sailoridy/MAESTRO


text:
! { dg-do run }
!
! executing simple SELECT TYPE statements
!
! Contributed by Janus Weil <janus@gcc.gnu.org>

  type :: t1
    integer :: i
  end type t1

  type, extends(t1) :: t2
    integer :: j
  end type t2

  type, extends(t1) :: t3
    real :: r
  end type

  class(t1), pointer :: cp
  type(t1), target :: a
  type(t2), target :: b
  type(t3), target :: c
  integer :: i

  cp => a
  i = 0

  select type (cp)
  type is (t1)
    i = 1
  type is (t2)
    i = 2
  class is (t1)
    i = 3
  end select

  if (i /= 1) call abort()

  cp => b
  i = 0

  select type (cp)
  type is (t1)
    i = 1
  type is (t2)
    i = 2
  class is (t2)
    i = 3
  end select

  if (i /= 2) call abort()

  cp => c
  i = 0

  select type (cp)
  type is (t1)
    i = 1
  type is (t2)
    i = 2
  class default
    i = 3
  end select

  if (i /= 3) call abort()

end

max_stars_repo_name: maurizioabba/rose

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
subroutine voronoiproj(leniw,lenrw,colg,nrow,rw,dres,goxd,dvxd,gozd,dvzd,depz,&
                       nx,ny,nz,nd,ncells,hvratio,damp,iproj,dv)
      use lsmrModule, only:lsmr

      implicit none
      integer leniw,lenrw
      integer nx,ny,nz
!      integer iw(leniw)
      integer colg(lenrw),nrow(nd)
      real depz(nz)
      real rw(lenrw)
      integer ncells,acells
      real dv(*),dres(*)
      real goxd,gozd,dvxd,dvzd
      real damp
      real hvratio,cmb
      integer ndim,nd
      integer iproj

      real,parameter:: radius = 6371.0,ftol = 0.1,pi = 3.141592654
      integer ii,ix,iy,iz
      real,dimension(:),allocatable:: grow,gcol,subrow,dis,dws,xunknown
      real,dimension(:),allocatable:: lat,lon,rad,theta,phi,rrad,xpts,ypts,zpts
      real,dimension(:),allocatable :: rw_p,rwgp,norm
      integer,dimension(:),allocatable:: iw_p,row,col,iwgp,colgp
      integer idx
      integer maxnar,nzid
      integer iseed(4)
      real xs,ys,zs
      real rx

      real atol,btol
      real conlim
      integer istop
      integer itnlim
      real acond
      real anorm
      real arnorm
      real rnorm
      real xnorm
      integer localSize,nout,itn
      integer leniw_p,lenrw_p,leniwgp,lenrwgp
      integer start

      allocate(lat(nx-2),lon(ny-2),rad(nz-1))
      ndim = (nx-2)*(ny-2)*(nz-1)

      do ii = 1,nx-2
      lat(ii) = (goxd-(ii-1)*dvxd)*pi/180
      enddo

      do ii = 1,ny-2
      lon(ii) = (gozd+(ii-1)*dvzd)*pi/180
      enddo

      !cmb = radius - depz(nz-1)*hvratio
      do ii = 1,nz-1
      rad(ii) = radius-depz(ii)*hvratio
      !rad(ii) = cmb+depz(ii)*hvratio
      enddo

      allocate(theta(ncells),phi(ncells),rrad(ncells),norm(ncells))
      allocate(xpts(ncells),ypts(ncells),zpts(ncells),dis(ncells),xunknown(ncells))
      allocate(rw_p(ndim))
      allocate(iw_p(2*ndim+1),row(ndim),col(ndim),dws(ndim))

      iseed(1:3) = (/38,62,346/)
      iseed(4) = 2*iproj+1
      call slarnv(1,iseed,ncells,theta)
      theta = (gozd+theta*(ny-3)*dvzd)*pi/180
      call slarnv(1,iseed,ncells,phi)
      phi = pi/2-(goxd-phi*(nx-3)*dvxd)*pi/180
      call slarnv(1,iseed,ncells,rrad)
      rrad = radius-rrad*depz(nz-1)*hvratio

      ! adaptive cells based on dws, assume 1/2 of all ncells are used
      ! as adaptive cells
      dws = 0
      do ii = 1,lenrw
      dws(colg(ii)) = dws(colg(ii))+abs(rw(ii))
      enddo
      acells = int(ncells/2.0)
      do ii = ncells-acells,ncells
      call random_index(idx,dws)
      ix = mod(idx,nx-2)
      iy = idx/(nx-2)
      iz = idx/((nx-2)*(ny-2))
      theta(ii) = (gozd+(ix+1)*dvzd)*pi/180
      phi(ii) = pi/2-(goxd-(iy+1)*dvxd)*pi/180
      rrad(ii) = radius-depz(iz+1)*hvratio
      enddo

      xpts = rrad*sin(phi)*cos(theta)
      ypts = rrad*sin(phi)*sin(theta)
      zpts = rrad*cos(phi)

      idx = 0
      do iz = 1,nz-1
      do iy = 1,ny-2
      do ix = 1,nx-2
      xs = rrad(iz)*sin(pi/2-lat(ix))*cos(lon(iy))
      ys = rrad(iz)*sin(pi/2-lat(ix))*sin(lon(iy))
      zs = rrad(iz)*cos(pi/2-lat(ix))
      dis =  (xpts-xs)**2+(ypts-ys)**2+(zpts-zs)**2
      idx = idx+1
      col(idx) = (iz-1)*(nx-2)*(ny-2)+(iy-1)*(nx-2)+ix
      row(idx) = minloc(dis,1)
      enddo
      enddo
      enddo
      rw_p = 1.0
      leniw_p = 2*ndim+1
      lenrw_p = ndim
      iw_p(1) = ndim
      iw_p(2:ndim+1) = row
      iw_p(ndim+2:2*ndim+1) = col

      allocate(grow(ndim),gcol(nd),subrow(ncells))
      maxnar = int(0.6*nd*ncells)
      allocate(iwgp(maxnar*2+1),colgp(maxnar),rwgp(maxnar))

      nzid = 0
      do ii = 1,nd
      grow = 0
      start = sum(nrow(1:ii-1))
      do ix = 1,nrow(ii)
      grow(colg(start+ix)) = rw(start+ix)
      enddo
      !gcol = 0
      !gcol(ii) = 1.0
      !call aprod(2,nd,ndim,grow,gcol,leniw,lenrw,iw,rw)
      subrow = 0
      call aprod(1,ncells,ndim,grow,subrow,leniw_p,lenrw_p,iw_p,rw_p)
      do ix = 1,ncells
      if(abs(subrow(ix))>ftol) then
          nzid = nzid+1
          rwgp(nzid) = subrow(ix)
          iwgp(1+nzid) = ii
          colgp(nzid) = ix
      endif
      enddo
      enddo
      leniwgp = nzid*2+1
      lenrwgp = nzid
      iwgp(1) = lenrwgp
      iwgp(nzid+2:nzid*2+1) = colgp(1:nzid)

      norm = 0
      do ii=1,nzid
      norm(iwgp(1+ii+nzid)) = norm(iwgp(1+ii+nzid))+rwgp(ii)**2
      enddo

      do ii =1,ncells
      norm(ii) = sqrt(norm(ii)/nd+0.01)
      enddo

      do ii =1,nzid
      rwgp(ii) = rwgp(ii)/norm(iwgp(1+ii+nzid))
      enddo

      conlim = 50
      itnlim = 100
      atol = 1e-3/((dvxd+dvzd)*111.19/2.0*0.1) !1e-2
      btol = 1e-3/(dvxd*nx*111.19/3.0)!1e-3
      istop = 0
      anorm = 0.0
      acond = 0.0
      arnorm = 0.0
      xnorm = 0.0
      localSize = int(ncells/4)
      !damp = dampvel
      ! using lsmr to solve for the projection coefficients
      !print*, 'LSMR beginning ...'

      nout = -1
      !nout = 36
      !open(nout,file='lsmrout_sub.txt')

      call LSMR(nd, ncells, leniwgp, lenrwgp,iwgp,rwgp,dres,damp,&
      atol, btol, conlim, itnlim, localSize,nout,&
      xunknown, istop, itn, anorm, acond,rnorm, arnorm, xnorm)
      !close(nout)
      do ii = 1,ncells
        xunknown(ii) = xunknown(ii)/norm(ii)
      enddo
      norm = (norm**2-0.01)*nd
      do ii = 1,ncells
      if (norm(ii)<0.01) then
          call random_number(rx)
          xunknown(ii) = xunknown(ii)+rx-0.5
      endif
      enddo

      dv(1:ndim) = 0
      call aprod(2,ncells,ndim,dv,xunknown,leniw_p,lenrw_p,iw_p,rw_p)
      deallocate(grow,gcol,subrow)
      deallocate(theta,phi,rrad,dws,norm)
      deallocate(xpts,ypts,zpts,dis,xunknown)
      deallocate(iw_p,rw_p,row,col)
      deallocate(lat,lon,rad)
      deallocate(iwgp,colgp,rwgp)

contains
subroutine random_index( idx, weights )
    integer :: idx
    real, intent(in) :: weights(:)

    real x, wsum, prob

    wsum = sum( weights )

    call random_number( x )

    prob = 0
    do idx = 1, size( weights )
        prob = prob + weights( idx ) / wsum   !! 0 < prob < 1
        if ( x <= prob ) exit
    enddo
end subroutine random_index

end subroutine

max_stars_repo_name: ChetGoerzen/DSurfTomo

max_issues_repo_name: ChetGoerzen/DSurfTomo

max_forks_repo_name: ChetGoerzen/DSurfTomo


text:
** Copyright (c) 1989, NVIDIA CORPORATION.  All rights reserved.
**
** Licensed under the Apache License, Version 2.0 (the "License");
** you may not use this file except in compliance with the License.
** You may obtain a copy of the License at
**
**     http://www.apache.org/licenses/LICENSE-2.0
**
** Unless required by applicable law or agreed to in writing, software
** distributed under the License is distributed on an "AS IS" BASIS,
** WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
** See the License for the specific language governing permissions and
** limitations under the License.

*   Intrinsics and generics: trigonometric and hyperbolic functions.

	program p
	implicit complex (c), double precision (d), double complex(z)
	parameter(n=55)
	integer rslts(n), expect(n), ctoi, dtoi
	parameter (d_dr=0.174532925199432957692D-1)
	parameter (r_dr=0.174532925199432957692E-1)
	parameter (d_rd=0.572957795130823208769D+2)
	parameter r_rd=0.572957795130823208769E+2

	dtoi(d) = d * 1000 + .499
	ctoi(c) = 1000 * (real(c) + aimag(c))
	ztoi(z) = 1000 * (real(z) + dimag(z))
	d_dtor(d) = d_dr*d
	d_rtod(d) = d_rd*d
	r_dtor(r) = r_dr*r
	r_rtod(r) = r_rd*r

	data x3, xx3 / 2 * 3.0 /,  d1, cx / 1.0d0, (1.0,2.0) /
	data zx/(1.0d0,2.0d0)/

c --- tests 1 - 4:   SIN

	rslts(1) = sin(.5) * 1000
	rslts(2) = dsin(d1) * 1000
	rslts(3) = dtoi( sin(dble(.5)) )
	rslts(4) = ctoi( csin( cx ) )

c --- tests 5 - 8:   COS

	rslts(5) = nint(cos(0.0))
	rslts(6) = dtoi( dcos(.34907d0) )
	rslts(7) = .001 + cos(x3 - xx3) * 1000
	rslts(8) = ctoi( ccos( cx ) )

c --- tests 9 - 11:  TAN

	rslts(9) = tan(1 / (x3 + 1)) * 1000
	rslts(10) = 1000 * tan(dble(x3) - 2)
	rslts(11) = dtan(.25d0) * 1000

c --- tests 12 - 14: ASIN

	rslts(12) = nint(asin(.47942) * 100)
	rslts(13) = nint(dasin(dble(x3) - 2) * 10)
	rslts(14) = asin(1.0d0) * 10

c --- tests 15 - 17: ACOS

	rslts(15) = nint(acos(1.0) * 100)
	rslts(16) = 100 * dacos(.93969 * 1d0) + .2
	rslts(17) = acos(1.0d0) + x3 + .01

c --- tests 18 - 20: ATAN

	rslts(18) = atan(0.0) * 100
	rslts(19) = dtoi( atan(.54630d0) )
	rslts(20) = dtoi( -datan(.54630d0) )

c --- tests 21 - 23: ATAN2

	rslts(21) = atan2(.54630, 1.0) * 1000 + .1
	rslts(22) = atan2(1.09260d0, x3-1d0) * 1000 + .1
	rslts(23) = datan2(.25534d0/2, dble(x3/6)) * 1000 + .1

c --- tests 24 - 26: SINH

	rslts(24) = sinh(3.0)
	rslts(25) = sinh( dble(x3) + 1)
	rslts(26) = dsinh(dble(x3) + 1)

c --- tests 27 - 29: COSH

	rslts(27) = -cosh(2.09)
	rslts(28) = cosh(dble(2.09))
	rslts(29) = dcosh(1d1)

c --- tests 30 - 32: TANH

	rslts(30) = 1 / tanh(.02d0)
	rslts(31) = tanh(x3) * 1000
	rslts(32) = dtanh( dble(x3) ) * 1000

c --- tests 33 - 36:   SIN

	rslts(33) = sind(r_rtod(.5)) * 1000
	rslts(34) = dsind(d_rtod(d1)) * 1000
	rslts(35) = dtoi( sind(dble(r_rtod(.5))) )
	rslts(36) = ztoi( cdsin( zx ) )

c --- tests 37 - 40:   COS

	rslts(37) = nint(cosd(0.0))
	rslts(38) = dtoi( dcosd(d_rtod(.34907d0)) )
	rslts(39) = .001 + cosd(r_rtod(x3 - xx3)) * 1000
	rslts(40) = ztoi( cdcos( zx ) )

c --- tests 41 - 43:  TAN

	rslts(41) = tand(r_rtod(1 / (x3 + 1))) * 1000
	rslts(42) = 1000 * tand(d_rtod(dble(x3) - 2))
	rslts(43) = dtand(d_rtod(.25d0)) * 1000

c --- tests 44 - 46: ASIN

	rslts(44) = nint(r_dtor(asind(.47942)) * 100)
	rslts(45) = nint(d_dtor(dasind(dble(x3) - 2)) * 10)
	rslts(46) = d_dtor(asind(1.0d0)) * 10

c --- tests 47 - 49: ACOS

	rslts(47) = nint(r_dtor(acosd(1.0)) * 100)
	rslts(48) = 100 * d_dtor(dacosd(.93969 * 1d0)) + .2
	rslts(49) = d_dtor(acosd(1.0d0)) + x3 + .01

c --- tests 50 - 52: ATAN

	rslts(50) = r_dtor(atand(0.0)) * 100
	rslts(51) = dtoi( d_dtor(atand(.54630d0)) )
	rslts(52) = dtoi( -d_dtor(datand(.54630d0)) )

c --- tests 53 - 55: ATAN2

	rslts(53) = r_dtor(atan2d(.54630, 1.0)) * 1000 + .1
	rslts(54) = d_dtor(atan2d(1.09260d0, x3-1d0)) * 1000 + .1
	rslts(55) = d_dtor(datan2d(.25534d0/2, dble(x3/6))) * 1000 + .1

c --- check results:

	call check(rslts, expect, n)

	data expect / 479, 841, 479, 5125,
     +                1, 940, 1000, -1019,
c           --- tests 9 - 14:
     +                255, 1557, 255, 50, 16, 15,
c           --- tests 15 - 20:
     +                0, 35, 3, 0, 500, -499,
c           --- tests 21 - 26:
     +                500, 500, 250, 10, 27, 27,
c           --- tests 27 - 32:
     +                -4, 4, 11013, 50, 995, 995,
c           --- tests 33 - 55 copy 1-23
     +		 479, 841, 479, 5125,
     +           1, 940, 1000, -1019,
     +           255, 1557, 255, 50, 16, 15,
     +           0, 35, 3, 0, 500, -499,
     +           500, 500, 250 /
	end

ext: f

max_stars_repo_name: kammerdienerb/flang

max_issues_repo_name: kammerdienerb/flang

max_forks_repo_name: kammerdienerb/flang


!***********************************/
!	Name:内挿元となるドナー要素の検索・登録プログラム
!	Alias:Register4InterpolationMain
!	Description:ドナー要素の共有面積も計算する
!	Type:
!	Input:
!	Output:
!	Note:
!	Author:Akitaka Toyota
!	Date:2017.12.19
!	Update:2017.12.19
!	Other:
!***********************************/
subroutine Register4InterpolationMain(LG,iGrid,Geom)
    use StructVar_Mod
!    use LoopVar_Mod
    use ConstantVar_Mod , only:CoreNumberOfCPU
    use RegistVar_Mod4OMP

    implicit none
    type(LocalGrids), intent(inout) :: LG
    type(Geometry), intent(inout) :: Geom
    integer, intent(in) :: iGrid
    type(CountUpVertex) :: CUV

    integer :: iDonorGrid, iTimeCombination
    integer :: iTargetCell
    integer :: iKindOfDonor, iFindDonor
    integer :: iGlobalCellNumber
    integer :: iLoop

    type(PrivateVarOfRegist) :: PVR
    allocate(PVR%RG2L%CenterDistanceVector(3),PVR%RG2L%iCenter(2),PVR%RG2L%iCenter2(2))
    allocate(PVR%RIG%CellCenterCoords(3),PVR%RIG%GlobalGridNum(3))
    allocate(PVR%RSDILG%MWDIAG%AuxiliaryGridNum(3),PVR%RSDILG%SSDIG%CenterDistanceVector(3))
!内挿用のDonor格子を見つけて共有面積を保存する
!    write(6,*) "Registoration of Interpolation-Cells"
        call AllocCUV !CUVの確保

        LG%UG(iGrid)%GM%Interpolated = 0 !時間積分をスキップするフラグの初期化
        LG%GG%GASG%ContentsAS2U = 0 !initialize

    !格子ごとのOverlap情報の作成
    do iDonorGrid = 1, LG%TotalGridsNumber
    !1.Lp格子と，n+1段階のLx格子(x<p)またはn段階のLy格子(p=<y)格子について接触判定を行う
    !物体を扱えないため：La<Li<Lbについて，LiとLaはN+1同士，LiとLbはN+1とNとで接触判定を行うため，1つとして同じ接触判定はない
    !n+1同士の接触判定は両側から見たものとして登録することで処理を減らすAtoB is BtoA

        if(iDonorGrid < iGrid) then
            iTimeCombination = 0 !N+1同士
        else
            iTimeCombination = 1 !N+1とN
        end if

        call OContactDetermination(LG%UG(iGrid),LG%UG(iDonorGrid),LG%MG(iGrid),LG%MG(iDonorGrid),iTimeCombination)

    end do

!$omp parallel num_threads(CoreNumberOfCPU),shared(LG,iGrid,Geom),firstprivate(iTargetCell,iKindOfDonor,iFindDonor,CUV,iGlobalCellNumber,PVR,iLoop,iTimeCombination)
!$omp do
!内挿用のDonor格子を見つけて共有面積を保存する
    do iTargetCell=1, LG%UG(iGrid)%GI%RealCells !Li格子の全実要素について

        call InputCUV(1,0,iTargetCell,LG%UG(iGrid),LG%MG(iGrid)%NSR,CUV,LG%MG(iGrid))
        !グローバル格子における位置判定+登録を行う
        call RegisterIntoGlobal(iGrid,LG%GG,iGlobalCellNumber,Geom,CUV,iLoop,PVR%RIG) !CutCell入れるなら消す？

        iFindDonor = 0
        do iKindOfDonor=1, 3
            if(iKindOfDonor <= 2 .or. iFindDonor == 0) then !上位格子+物体あり格子からの補間が終わった時点でまだ補間要素が見つかっていない時下位格子の検索を行う

                do iDonorGrid=1, LG%TotalGridsNumber !計算領域中の全格子について
                    if(iDonorGrid < iGrid) then
                        iTimeCombination = 0 !N+1同士
                    else
                        iTimeCombination = 1 !N+1とN
                    end if

                    if(LG%UG(iGrid)%Overlap(iTimeCombination+1,iDonorGrid) == iKindOfDonor) then !iKind=1で上位格子，iKind=2で物体ありの下位格子，iKind=3で下位格子
                        call RoughSearchDonorInLocalGrid(iTimeCombination,LG%MG(iDonorGrid),LG%UG(iDonorGrid),LG%MG(iGrid),CUV,iFindDonor,iLoop,PVR%RSDILG)
                        call InputCUV(1,0,iTargetCell,LG%UG(iGrid),LG%MG(iGrid)%NSR,CUV,LG%MG(iGrid)) !Donor格子が変わるたびにその座標系に合わせる必要があるため，毎回グローバル座標系の値に戻す
                    end if

                end do

            else !上位格子+物体あり下位格子の検索が終わった時点で補間用要素が見つかっているとき = 積分不要
                LG%UG(iGrid)%GM%Interpolated(iTargetCell,1,1) = 1
            end if
        end do

    !この時点で要素と被る可能性のあるすべてのローカル格子との接触判定が終わっている
    !そのため，この時点でiFindDonorが0であるならばグローバル格子からの補間が必要となる
        if(LG%UG(iGrid)%GM%CellType(iTargetCell,1,1) == 99) then
            call RegisterGlobal2Local(iGlobalCellNumber,LG%MG(iGrid),iFindDonor,1,CUV,Geom,9999,iLoop,PVR%RG2L)
            iFindDonor = 1

        end if

        if(iFindDonor == 0) then

            call RegisterGlobal2Local(iGlobalCellNumber,LG%MG(iGrid),iFindDonor,0,CUV,Geom,9999,iLoop,PVR%RG2L)

        end if

    end do
!$omp end do
!$omp end parallel
    deallocate(PVR%RG2L%TempDE,PVR%RIG%Temp,PVR%RSDILG%SSDIG%ASDCN%TempDE,PVR%RSDILG%TempS2U)

    do iTargetCell = LG%UG(iGrid)%GI%RealCells+1, LG%UG(iGrid)%GI%AllCells !すべての仮想セルについて
        !仮想セルのセル中心座標から下位格子のデータを拾ってくる
        !仮想セルの補間は当初想定していなかったために大きめの修正なしでCUVが使えない．そのため，暫定的にグローバルから拾うだけに留める
        CUV%TargetCoords(0,1:3) = LG%MG(iGrid)%NSR%C%NextCoordsInG(iTargetCell,1:3)
        CUV%DonorCoords(0,1:3) = LG%MG(iGrid)%NSR%C%NextCoordsInG(LG%UG(iGrid)%VC%Cell(iTargetCell,1),1:3)
        CUV%TargetCoords(1,1:3) = LG%UG(igrid)%GM%Normal(LG%UG(iGrid)%VC%Edge(iTargetCell),1:3) !Normal Vector On Boundary Surface

        call RegisterGlobal2Local(iGlobalCellNumber,LG%MG(iGrid),9999,2,CUV,Geom,iTargetCell,iLoop,PVR%RG2L)
    end do

    !if(LG%UG(iGrid)%InternalBoundary /= 0) then
    !    call ZeroFill
    !end if

    !VC Loop:To Search G-Cell what Overlap with OuterBoundary Of L-Cell
    !and nullify Outer Boundary G-Cell's Interpolation relationship
    call ExcludeImperfectOverlapped(LG%UG(iGrid),LG%MG(iGrid),LG%GG,Geom)
    call ReRegistBoundaryCell(LG%UG(iGrid),LG%MG(iGrid),Geom)

return
contains

    subroutine AllocCUV
    implicit none

        allocate(CUV%AdditionVertex(10,3))
        allocate(CUV%AuxiliaryCellNumber(9))
        allocate(CUV%DonorCoords(0:3,4))
        allocate(CUV%iDonorNumber(0:3))
        allocate(CUV%iTargetNumber(0:3))
        allocate(CUV%TargetCoords(0:3,4))
        CUV%TargetCoords(:,4) = 1.0d0

    return
    end subroutine AllocCUV

    subroutine ZeroFill !実際使おうとするとかなり面倒なので却下
    use FrequentOperation
    implicit none
    integer :: iCenterX, iCenterY, iGlobalNum, GetGlobalNumber2D
    type(DonorElement), pointer :: TempDE
    double precision, allocatable :: CenterCoords(:), DistanceVector(:)

    allocate(CenterCoords(3),DistanceVector(3))

        do iCenterY = 1, Geom%CellNumber(2)
            do iCenterX = 1, Geom%CellNumber(1)
                call GetCellCenterCoords(iCenterX,iCenterY,Geom%Bound(1,1:2),2.0d0*Geom%Width(1,1,1:2),CenterCoords(1:2))
                CenterCoords(3) = 0.0d0

                DistanceVector = CenterCoords(1:3) - LG%MG(iGrid)%RC%GravityCenter(2,1:3)

                if(LG%UG(iGrid)%InternalRadius > sqrt(dot_product(DistanceVector,DistanceVector))) then !物体の内部とかぶっているなら(クソ雑判定)
                    iGlobalNum = GetGlobalNumber2D(CenterCoords(1:2),2.0d0*Geom%Width(1,1,1:2),Geom%Bound(1,1:2),Geom%CellNumber(1))

                    Geom%Interpolated(iCenterX,iCenterY,1) = 1
                    nullify(LG%GG%GDE(iGlobalNum)%DE)
                    allocate(TempDE) !新規ノードの確保
                        TempDE%Cell = 0 !ノードへの記録
                        TempDE%Grid = -1 !物体の内側
                        TempDE%SharedArea = 1.0d0

                        TempDE%Next => LG%GG%GDE(iGlobalNum)%DE !前のノードへのリンク作成
                        LG%GG%GDE(iGlobalNum)%DE => TempDE !データの保存

                end if
            end do
        end do

        return
    end subroutine ZeroFill

end subroutine Register4InterpolationMain


text:
      subroutine ini_optcld(Qxv,Qxi,Qsv,Qsi,gv,gi,Qextref,TAUREF_CLD)

!     Legacy Mars GCM v24
!     Mars Climate Modeling Center
!     NASA Ames Research Center

      use grid_h
      use radinc_h
      use dtcommon_h
      use cldcommon_h

      implicit none

c  Arguments
c  ---------

      real*8  Qxv(L_LEVELS+1,L_NSPECTV)
      real*8  Qsv(L_LEVELS+1,L_NSPECTV)
      real*8  gv(L_LEVELS+1,L_NSPECTV)

      real*8  Qxi(L_LEVELS+1,L_NSPECTI)
      real*8  Qsi(L_LEVELS+1,L_NSPECTI)
      real*8  gi(L_LEVELS+1,L_NSPECTI)

      real*8  Qextref(L_LEVELS+1)
      real*8  TAUREF_CLD(L_LEVELS+1)

c  Local variables
c  ---------------

      integer i,k

c Initialyze various variables
c ----------------------------

      DO K = 1, L_LEVELS+1
        Qextref(K)    = 1.
      ENDDO

      call settozero(L_LEVELS+1,TAUREF_CLD)

      call settozero(nlonv*(L_LEVELS+1),Qxv)
      call settozero(nlonv*(L_LEVELS+1),Qsv)
      call settozero(nlonv*(L_LEVELS+1),gv)

      call settozero(nloni*(L_LEVELS+1),Qxi)
      call settozero(nloni*(L_LEVELS+1),Qsi)
      call settozero(nloni*(L_LEVELS+1),gi)

      return

      end

max_stars_repo_name: nasa/legacy-mars-global-climate-model

max_issues_repo_name: nasa/legacy-mars-global-climate-model

max_forks_repo_name: nasa/legacy-mars-global-climate-model


text:
MODULE bc_variables

  implicit none

  PRIVATE
  PUBLIC :: n_neighbor, fdinfo_send, fdinfo_recv
  PUBLIC :: www, sbuf, rbuf
  PUBLIC :: TYPE_MAIN, zero, Md

  include 'mpif.h'

  integer :: n_neighbor(6)
  integer,allocatable :: fdinfo_send(:,:,:),fdinfo_recv(:,:,:)

#ifdef _DRSDFT_
  integer,parameter :: TYPE_MAIN=MPI_REAL8
  real(8),allocatable :: www(:,:,:,:)
  real(8),allocatable :: sbuf(:,:,:),rbuf(:,:,:)
  real(8),parameter :: zero=0.0d0
#else
  integer,parameter :: TYPE_MAIN=MPI_COMPLEX16
  complex(8),allocatable :: www(:,:,:,:)
  complex(8),allocatable :: sbuf(:,:,:),rbuf(:,:,:)
  complex(8),parameter :: zero=(0.d0,0.d0)
#endif

  integer :: Md

END MODULE bc_variables

max_stars_repo_name: j-iwata/RSDFT

max_issues_repo_name: j-iwata/RSDFT

max_forks_repo_name: j-iwata/RSDFT


text:
!----------------------------------------------------------------------------- best with 100 columns

function matVecMul()
  use modSparse
  integer::matVecMul
  type(linEq)::p
  integer::iA(7),jA(7)
  double precision::A(7)
  double precision::y(4),x(4)

  matVecMul=0
  call p%initLinEq(4,7)
  iA=[1,2,2,3,3,4,4]
  jA=[1,2,1,3,4,4,2]
  A=[1d0,2d0,-1d0,5d0,-4d0,2d0,-1d0]
  x=[1d0,3d0/2d0,14d0/5d0,11d0/4d0]
  call p%setCOO(iA,jA,A,7)
  call p%mulVec(x,y)
  if(norm2(y-[1d0,2d0,3d0,4d0])>1d-9)then
    matVecMul=1
  end if
end function

max_stars_repo_name: mianzhi/fosolvers

max_issues_repo_name: mianzhi/fosolvers

max_forks_repo_name: mianzhi/fosolvers


text:
       PROGRAM  POLYSCRB
C#######################################################################
C  NAME OF ROUTINE
C      POLYSCRB  ( POLYgon SCRiBe )
C
C  PURPOSE
C      THIS IS THE STANDARD MAIN PROGRAM USE FOR TAE/VICAR PROGRAMS.
C      THIS MODULE CALLS SUBROUTINE MAIN44 TO ENTER INTO THE BODY OF THE
C      PROGRAM.
C  PREPARED FOR USE ON MIPL SYSTEM BY
C      STEVE POHORSKY   INFORMATICS GENERAL CORPORATION     AUGUST 1983
C  FOR
C      RON ALLEY and JERRY SOLOMON
C
C  ENVIRONMENT
C      VAX 11/780    VMS 3.0 with TAE/VICAR1 EXECUTIVE       FORTRAN-77
C
C  REVISION HISTORY
C
C  CALLING SEQUENCE (TAE COMMAND LINE)
C      The following command line formats show the major allowable forms:
C      POLYSCRB INP=a OUT=b SIZE=(sl,ss,nl,ns) optional parameters
C      POLYSCRB INP=a OUT=b SL=sl SS=ss NL=nl NS=ns optional parameters
C      POLYSCRB a b (sl,ss,nl,ns) optional parameters
C
C       Here 'a' represents the input file name,
C       'b' represents the output image file name.
C
C  INPUT PARAMETERS (listed by keyword)
C      INP    - Input file name.
C      OUT    - Output file name followed by work file name.
C      SIZE   - Standard Vicar size field:  (SL,SS,NL,NS)
C               SL = Starting line number.
C               SS = Starting sample number.
C               NL = Number of lines.
C               NS = Number of samples.
C      HALF   - halfword output flag.
C      DN     - Data number for polygon borders
C      DIM    - Specifies 2 or 3-dimensional graphics-1 file  6/2/88 NDR
C      BDN    - Data number for background areas
C      PDN    - Data number for polygon corners (endpoints of line segments)
C      EXCLOSE- Flag to check for polygons which close exactly
C               when determining the end of sequences of points in the
C               input file.
C      APCLOSE- Distances in checking for polygons which close approximately.
C      OUTSIDE- Flag to expand the polygon border by one pixel
C               around the outside of the polygon.
C      IGNOR  - coordinates of a null point.
C
C  OUTPUT PARAMETERS
C      The output image produced is written to the output file.
C  PROGRAM LIMITATIONS
C      1. The output image can be byte or halfword data.
C      2. Maximum number of samples is 500,000 per line.  The TAE/VICAR
C         executive may impose further size limitations.
C  SUBROUTINES CALLED
C      MAIN44
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      INCLUDE 'VICMAIN_FOR'
      SUBROUTINE MAIN44
C#######################################################################
C  NAME OF ROUTINE
C     MAIN44  (name of top level subroutine by VICAR convention)
C
C  PURPOSE
C      MAIN44 takes the data from a standard Graphics-1 vector
C      input file and produces from it an output image file of polygon borders
C      (line segments) that are scribed (drawn) against a uniform
C      background.
C
C  CONVERTED FOR USE ON MIPL SYSTEM BY
C      STEVE POHORSKY   INFORMATICS GENERAL CORPORATION     AUGUST 1983
C  FOR
C      RON ALLEY and JERRY SOLOMON
C
C  ORIGINAL POLYSCRB PROGRAM BY
C      AL ZOBRIST
C  ENVIRONMENT
C      VAX 11/780    VMS 3.0 with TAE/VICAR1 EXECUTIVE       FORTRAN-77
C
C  REVISION HISTORY
C     8-83  SP   CONVERTED FROM IBM VICAR VERSION: MISCELLANEOUS CLEANUP.
C     8-83  SP   INCREASED WORK BUFFER SIZE TO 500000 BYTES.
C     8-83  SP   CORRECTED LABEL PROCESSING OF INPUT FILE BY USING LABELC.
C     8-83  SP   CHANGED LABEL PROCESSING OF OUTPUT FILE TO INCLUDE CORRECT
C                DATA FORMAT IN LABEL.
C     8-83  SP   IMPLEMENTED BORIS GOKHMAN'S CHANGES TO THE IBM POLYSCRB WHICH
C                ALLOW POLYSCRB TO PRODUCE HALFWORD OUTPUT. ADDED HALF KEYWORD.
C     8-83  SP   CORRECTED INITIALIZATION OF OX2 AND OY2 WHEN IGNOR PARAMETERS
C                ENTERED.
C     8-83  SP   CHANGED REFERENCES TO DN(4), BDN(4), PDN(4) TO  DN(1),
C                BDN(1), PDN(4) BECAUSE VAX STORES LOW-ORDER BYTE IN FIRST
C                BYTE.
C     5-84  SP   USED XVSIZE INSTEAD OF PARAM(2),PARAM(4) FOR SIZE
C                PARAMETER BECAUSE PARAM DOES NOT KNOW SIZE PARAMETER DOES NOT
C                REFER TO INPUT FILE.
C     6-88  NR   ADDED DIM OPTION FOR SCRIBING 3-DIMENSIONAL FILES
C     1-95  AS   ...CRI... MSTP S/W CONVERSION (VICAR PORTING)
C
C  CALLING SEQUENCE
C      Standard subroutine call and return.
C  INPUT AND OUTPUT PARAMETERS
C      SEE UNDER PROGRAM POLYSCRBMAIN.
C
C  CALLED BY
C      POLYSCRBMAIN
C  SUBROUTINES CALLED
C      The library routines  ABEND, KSCAN, LABELC, OPEN,
C      PARAM, PRNT, XVMESSAGE, RDCHEK, READ, WCHECK, WRITE.
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C                 ORIGINAL HEADER COMMENTS.
C
C  IBIS ROUTINE POLYSCRB
C
C  POLYSCRB IS USED TO TRANSFORM A STANDARD POLYGON FILE INTO AN IMAGE
C  FILE OF POLYGON BORDERS.  THE INPUT DATA IN COORDINATE-POINT FORMAT
C  IS USED TO FORM A RASTER BASE CONTAINING POLYGON OUTLINES.  BOTH THE
C  DN VALUES FOR THE POLYGON BORDERS AND THE BACKGROUND AREAS MAY BE
C  SPECIFIED.  PARAMETERS FOR CHAINING AND CLOSURE MAY BE INVOKED.
C
C  USER PARAMETERS:
C
C  DN,N -    SPECIFIES THE DN VALUE WHICH WILL BE USED FOR SCRIBING POLY
C            BOUNDARIES.  THE DEFAULT IS 255.
C  EXCLOSE - INDICATES THAT THE BEGINNING AND ENDING NODES OF EACH POLYG
C            MUST MATCH EXACTLY.
C  APCLOSE,X,Y - THE FLOATING POINT NUMBERS X,Y SPECIFY THE MINIMUM
C            ACCEPTABLE X AND Y DISTANCES IN THE APPROXIMATE POLYGON CLO
C            OPTION.
C  BDN,N -   SPECIFIES THE DN VALUE WHICH WILL BE USED AS BACKGROUND.
C            THE DEFAULT IS ZERO.
C  OUTSIDE - SPECIFIES THAT THE BOUNDARY OF RIGHT HAND CODED POLYGONS AR
C            EXPANDED BY ONE PIXEL.
C  IGNOR,X,Y - CAUSES ALL OCCURENCES OF THE POINT (X,Y) TO BE TREATED AS
C            NULL POINTS.  A NULL POINT INVOKES AUTOMATIC POLYGON CLOSUR
C            THE DEFAULT IS (0.,0.).
C  PDN,P -   SPECIFIES A DN VALUE TO BE PLACED AT THE ENDS OF EACH LINE
C            SEGMENT OR THE CORNERS OF EACH POLYGON.
C
C REVISIONS
C   04 NOV 81   ...ALZ...   FIX KSCAN EXIT FOR NO PARAMS AGAIN
C   10 AUG 81   ...ALZ...   FIX KSCAN EXIT FOR NO PARAMS
C   14 FEB 80   ...ALZ...    REPLACE CALL END BY RETURN
C   2 FEB 79   ...ALZ...    INITIAL RELEASE
C
C                 END OF ORIGINAL HEADER COMMENTS.
C
      IMPLICIT NONE

      INTEGER WORKSIZE_PAR, HWRKSIZE_PAR
      PARAMETER ( WORKSIZE_PAR = 500000 )   ! NUMBER OF BYTES IN WORK BUFFER
      PARAMETER ( HWRKSIZE_PAR = WORKSIZE_PAR/2 )

      REAL ABS,FLOAT,AMIN1,AMAX1,Z
      REAL FXL,FXU,FYL,FYU,LINES,X1,X2,Y1,Y2,X,Y,XD,YD,TST(4),T,T1,T2
      REAL XINC,YINC,TL,TU,XQ,YQ,XFIRST,YFIRST,LEPS,SEPS
      REAL PX,PY,OX2,OY2,SQ,XT,YT,PMJ,PMN,CYCLE,XIGNOR,YIGNOR
      REAL TEMP_REAL(2), maxX, maxY, max
      BYTE WORK,BDN
      LOGICAL INLINE,CONTIN,OUTSID,LPRNT,HALF,NOPRINT,XVPTST,USEZVALUE

      INTEGER   GDIM, IDN, IBDN, IPDN, PTR, WSIZ, COUNT, CLOSR,TEMP_INT
      INTEGER   NUM_VALUES, DEFAULT_USED, ERROR_NUM, SIZE(4), NINT, I
      INTEGER   STATUS, ISL, ISS, NL, NS, SLINE, NLINE, NLINEU, SSAMP
      INTEGER   NSAMP, SPIXEL, NSAMPU, NSAMP4, LU, LD, LL, WSUP, ISC
      INTEGER   JSC, IST, JST, ISZ, JSZ, ISTZ, JSTZ, NSQT, INDEX, LINDEX
      INTEGER   ZINDEX, XBIT, YBIT, F2H_BUF(12), OUT_UNIT
      INTEGER*2 HWORK(HWRKSIZE_PAR),  HDN,  HBDN, HPDN
      LOGICAL   ALL_ZERO_SW, EOF_SW, SAMEPDN
      INTEGER   RDGR, GETGR, SETGR, CLGR

      include 'fortport'

      COMMON /COM/WORK(WORKSIZE_PAR)
      COMMON /COM1/LINES(900)

       ! NOTE THAT ON THE VAX THE FIRST BYTE (NOT THE FOURTH) OF A WORD
       ! IS THE LOW ORDER BYTE.

      DATA XIGNOR,YIGNOR/2*0./,IPDN/255/
      DATA IDN,IBDN/255,0/,CONTIN,OUTSID/2*.FALSE./,CLOSR/-1/
      DATA OX2,OY2,XFIRST,YFIRST/4*0./,SEPS,LEPS/.1,.1/,LPRNT/.TRUE./

C
C=================START OF EXECUTABLE CODE===============================
C
      CALL IFMESSAGE('POLYSCRB version 2-JAN-95')
      CALL XVEACTION('SA',' ')
C
C  INITIALIZE, GET PARAMETERS, OPEN FILES.
C
      CALL XVTRANS_SET(F2H_BUF,'FULL','HALF',STATUS)
      IF (STATUS.NE.1) CALL MABEND('BUFFER SETUP UNSUCCESSFUL')

      WSIZ = WORKSIZE_PAR                ! SIZE OF WORK BUFFER IN PIXELS.
      HALF = .FALSE.

	if (xvptst('HALF')) then
          HALF = .TRUE.
          WSIZ = HWRKSIZE_PAR
        END IF

	CALL XVP ('DIM',GDIM,COUNT)

	call xvparm ('DN', temp_int, num_values, default_used,1)
	USEZVALUE  = (DEFAULT_USED .EQ. 1) .AND. (GDIM .EQ. 3)
	if (.NOT. USEZVALUE) idn = temp_int

	call xvparm ('BDN', temp_int, num_values, default_used,1)
	if (default_used .eq. 0) ibdn = temp_int

        IPDN = IDN
	call xvparm ('PDN', temp_int, num_values, default_used,1)
	SAMEPDN = (DEFAULT_USED .EQ.  1) .AND. (USEZVALUE)
	if (.NOT. SAMEPDN) ipdn = temp_int

	outsid = xvptst('OUTSIDE')

	if (xvptst('EXCLOSE')) closr = 0

	call xvparm ('IGNOR', temp_real, num_values, default_used,2)
	if (default_used .eq. 0) then
	    xignor = temp_real(1)
	    yignor = temp_real(2)
            OX2    = XIGNOR
            OY2    = YIGNOR
        END IF

	call xvparm ('APCLOSE', temp_real, num_values, default_used,2)
	if (default_used .eq. 0) then
	    CLOSR = 1

!	 If one value is supplied for the approximate closure (num_values
!	 = 1), then assign it to the line and sample closure. If a second
!	 value is supplied, then assign it to the sample closure.
	    if (num_values .ge. 1) then
		leps = temp_real(1)
		seps = temp_real(num_values)
	    endif
	endif

      NOPRINT=XVPTST('NOPRINT')

! Open the input graphics file. We'll open the first output file, call it
! the first graphics file, and we'll use 2 coordinates per point.
!********************************************************************
!********** UNLESS DIM = 3 ******************************************
!********************************************************************

	status= rdgr(1, 1, GDIM)
        if (status.ne.1) call signalgr(1,status,1)

! Get the SIZE parameters
	call xvparm ('SIZE', size, num_values, default_used,4)
	if (default_used .eq. 0) then
	    isl = size(1)
	    iss = size(2)
	    nl = size(3)
	    ns = size(4)
	else
	    call xvp ('SL', isl, num_values)
	    call xvp ('SS', iss, num_values)
	    call xvp ('NL', nl, num_values)
	    call xvp ('NS', ns, num_values)
	endif

! If number of lines and number of samples are defaulted, and thus have a
! value of 0, then run through the graphics file and get the values that
! will enclose the maximum graphics coordinates. Subtract the starting
! line and starting sample from this maximum to get values for NL and NS
! that will display all of the graphics when SL = SS = 1. Graphics with
! line (x) values less than SL or sample (y) values less than SS will be
! chopped. We'll leave the graphics file after the last coordinate set
! because there's a SETGR call made before the graphics file is re-used.
	if (nl .le. 0  .or.  ns .le. 0) then
	    maxX = isl
	    maxY = iss
	    eof_sw = .false.
	    do while (.not. eof_sw)		! go thru file until eof
	       	status= getgr(1, all_zero_sw, eof_sw, x, y, Z)
                if (status.ne.1) call signalgr(1,status,1)
		if (.not. all_zero_sw .and. .not. eof_sw) then
		    maxX = max(maxX, x)
		    maxY = max(maxY, y)
		endif
	    enddo

!     Convert from absolute maximums to window widths
	    maxX = maxX - isl + 1
	    maxY = maxY - iss + 1

!     Modify the defaults, rounding the real-valued window widths to integers
	    if (nl .le. 0) nl = nint(maxX)
	    if (ns .le. 0) ns = nint(maxY)
	endif

	sline = isl
	nline = nl
	nlineu = nline + sline - 1

        SSAMP = ISS
        NSAMP = NS

         SPIXEL = SSAMP
         NSAMPU = NSAMP+SSAMP-1
         NSAMP4 = NSAMP+4-MOD(NSAMP,4)  ! NSAMP4 PROVIDES MARGIN FOR SCRIBING
					! AT BORDERS,ALSO FOR FULLWORD ALIG

! Open the output image file, setting it to BYTE or HALF, depending on
! whether half is false or true. Use an internal buffer that matches the
! image type so that we can use as big an addressable buffer as possible.

! Set the image type to "IMAGE" to avoid confusing some other programs-NDR.
	call xvunit (out_unit, 'OUT', 1, error_num,' ')
	if (half) then
	    call xvopen (out_unit, error_num, 'U_FORMAT', 'HALF',
     +		 'O_FORMAT', 'HALF', 'U_NL', nline, 'U_NS', nsamp,
     +		 'OP', 'WRITE','TYPE','IMAGE',' ')
	else
	    call xvopen (out_unit, error_num, 'U_FORMAT', 'BYTE',
     +		 'O_FORMAT', 'BYTE', 'U_NL', nline, 'U_NS', nsamp,
     +		 'OP', 'WRITE','TYPE','IMAGE',' ')
	endif

C
C  SET UP ONE 'STRIP' OR 'LAYER' OF IMAGE FOR SCRIBING THE LINE DATA SET.
C  AS MUCH AS CAN BE HELD BY WORK VECTOR.  HERE X VALUES REFER
C  TO LINE POSITON, AND Y VALUES REFER TO PIXEL POSITION.
C
      LU = SLINE-1
      LD = WSIZ/NSAMP4              ! NUMBER OF LINES IN STRIP.
      WSUP = LD*NSAMP4
      FYL = FLOAT(SPIXEL)-.6        ! LEFT EDGE OF WINDOW.
      FYU = FLOAT(NSAMPU)+.6        ! RIGHT EDGE OF WINDOW.

C  LOOP TO BUILD OUTPUT FILE ONE STRIP AT A TIME.

 100  LL = LU+1                    ! FIRST LINE OF STRIP.
      LU = MIN0(LU+LD,NLINEU)      ! LAST LINE OF STRIP.
      IF (LL.GT.LU)  GOTO  8000    ! ALL DONE IF WE'RE PAST LAST LINE OF IMAGE.

      FXL = FLOAT(LL)-.6           ! TOP OF STRIP.
      FXU = FLOAT(LU)+.6           ! BOTTOM OF STRIP.

      IF ( HALF )   THEN
         CALL XVTRANS(F2H_BUF,IBDN,HBDN,1)
         DO I = 1, HWRKSIZE_PAR    ! FILL BUFFER WITH BACKGROUND PIXELS.
           HWORK(I) = HBDN
         END DO
      ELSE
         BDN=INT2BYTE(IBDN)
         DO I = 1, WORKSIZE_PAR
           WORK(I) = BDN
         END DO
      END IF

! Reposition the input graphics file to the first coordinate set in
! the file
	status= setgr(1, 1)
        if (status.ne.1) call signalgr(1,status,1)

C
C  READ THE ENTIRE LINES DATA SET.  LINES THAT CROSS THE STRIP ARE PROCES
C
      do while (.true.)

      X1 = OX2               ! TERMINATING POINT OF LAST SEGMENT IS
      Y1 = OY2               ! BEGINNING POINT OF THIS SEGMENT.
      status= getgr(1, all_zero_sw, eof_sw, x2, y2, Z) ! coords of next point
      if (status.ne.1) call signalgr(1,status,1)
      if (eof_sw) go to 91

      IF (USEZVALUE) IDN = NINT(Z)
      IF (SAMEPDN) IPDN = IDN

      OX2 = X2	! save coordinates of end of this segment for use as
      OY2 = Y2	! the beginning of the next segment

C  CHECK FOR NULL POINT.

      IF (ABS(X1-XIGNOR)+ABS(Y1-YIGNOR).LE.1.E-6) GO TO 88
      IF (ABS(X2-XIGNOR)+ABS(Y2-YIGNOR).LE.1.E-6) GO TO 88

C   CLOSR = -1   NO CLOSURE OPTION
C            0   EXCLOSE              CHECK FOR POLYGON CLOSURE.
C            1   APCLOSE              CHECK FOR POLYGON CLOSURE.

      IF (CLOSR) 15,8,7
 7    IF (ABS(X1-XFIRST).GT.LEPS) GO TO 9
      IF (ABS(Y1-YFIRST).GT.SEPS) GO TO 9
      GO TO 88
 8    IF (X1.EQ.XFIRST.AND.Y1.EQ.YFIRST) GO TO 88
 9    IF (CONTIN) GO TO 15
      CONTIN = .TRUE.
      XFIRST = X1
      YFIRST = Y1
 15   IF (.NOT.OUTSID) GO TO 16
      PX = Y1-Y2
      PY = X2-X1
      SQ = SQRT(PX*PX+PY*PY)       ! JUST LIKE IN THE ANALYTIC GEOMETRY
      IF (SQ.LE..01) GO TO 16      ! BOOK.
      PX = PX/SQ
      PY = PY/SQ
      X1 = X1+PX                   ! SHIFT LINE SEGMENT ONE PIXEL
      Y1 = Y1+PY                   ! IN PERPINDICULAR DIRECTION.
      X2 = X2+PX		   !        ^ SPELLING -- CHANGE TO "E"
      Y2 = Y2+PY
C  QUICK TEST FOR INTERSECTION
 16   IF ((X1-FXL).LT.0..AND.(X2-FXL).LT.0.) GO TO 89
      IF ((X1-FXU).GT.0..AND.(X2-FXU).GT.0.) GO TO 89
C  MORE COMPLETE TEST FOR INTERSECTION.  ALSO FIND INITIAL SQUARE.

      IF (ABS(X2).GT.1.E6.OR.ABS(Y2).GT.1.E6) GO TO 401  ! CRAZY NUMBERS?

      IF (X1-X2) 19,17,18
 17   IF (Y1-Y2) 19,19,18
 18   X = X1
      X1 = X2                ! SWAP ENDPOINTS SO X2 >= X1
      X2 = X
      Y = Y1
      Y1 = Y2
      Y2 = Y

C  THE STRIP IS A RECTANGLE. FIND WHERE THE LINE THROUGH THE ENDPOINTS
C  INTERSECTS THE LINES THROUGH EACH OF THE FOUR SIDES OF THE STRIP.
C

 19   XD = (X2-X1)+1.141593E-9   ! OBSCURE NUMBERS DESIGNED TO PREVENT
      YD = (Y2-Y1)+1.141593E-9   ! OVERFLOW ON DIVIDE.
      IF ( ABS(YD)  .LE. 1.E-10 )   YD = 1.E-10        ! EXTRA CARE.
      TST(1) = (FXL+.1-X1)/XD
      TST(2) = (FYL+.1-Y1)/YD
      TST(3) = (FXU-.1-X1)/XD
      TST(4) = (FYU-.1-Y1)/YD
      TL = 3.
      TU = -3.
      INLINE = .FALSE.
      DO 4 I=1,4
      X = XD*TST(I)+X1
      Y = YD*TST(I)+Y1
      IF (X.LT.FXL.OR.X.GT.FXU) GO TO 4
      IF (Y.LT.FYL.OR.Y.GT.FYU) GO TO 4
      INLINE = INLINE.OR.(TST(I).GE.0..AND.TST(I).LE.1.)
      TL = AMIN1(TL,TST(I))
      TU = AMAX1(TU,TST(I))
 4    CONTINUE
      IF (TL.GT.2.) GO TO 89
      T = AMAX1(TL,0.)
      TU = AMIN1(TU,1.)
C
C  CHECK FOR LINE SEGMENT BETWEEN POINTS INTERSECTING A SIDE OF THE STRIP
C  (INLINE = .TRUE. ) OR FOR AN ENDPOINT IN THE INTERIOR OF STRIP.
C
      IF (INLINE) GO TO 5
      IF (X1.GE.FXL.AND.X1.LE.FXU.AND.Y1.GE.FYL.AND.Y1.LE.FYU) GO TO 5
      IF (X2.GE.FXL.AND.X2.LE.FXU.AND.Y2.GE.FYL.AND.Y2.LE.FYU) GO TO 5
      GO TO 89
 5    X = XD*T+X1
      Y = YD*T+Y1
      ISC = MAX0(MIN0(INT(X+.5),LU),LL)-LL+1
      JSC = MAX0(MIN0(INT(Y+.5),NSAMPU),SPIXEL)-SPIXEL+1
      XT = XD*TU+X1
      YT = YD*TU+Y1
      IST = MAX0(MIN0(INT(XT+.5),LU),LL)-LL+1
      JST = MAX0(MIN0(INT(YT+.5),NSAMPU),SPIXEL)-SPIXEL+1
      NSQT = MAX0(IABS(IST-ISC),IABS(JST-JSC))
      INDEX = (ISC-1)*NSAMP4+JSC
      LINDEX = (IST-1)*NSAMP4+JST

      IF (HALF)  THEN
          CALL XVTRANS(F2H_BUF,IDN,HDN,1)
          CALL XVTRANS(F2H_BUF,IPDN,HPDN,1)
          HWORK(INDEX) = HDN
          IF (T.LT..0001) HWORK(INDEX) = HPDN
      ELSE
          WORK(INDEX) = INT2BYTE(IDN)
          IF (T.LT..0001) WORK(INDEX) = INT2BYTE(IPDN)
      END IF

      IF (NSQT.EQ.0) GO TO 89
      XQ = FLOAT(ISC+LL-1)
      YQ = FLOAT(JSC+SPIXEL-1)
C  FINE ADJUSTMENT OF TRAJECTORY OF LINE.
      XINC = ABS(1./XD)
      YINC = ABS(1./YD)
      XBIT = NSAMP4
      YBIT = 1
      IF (XD.LT.0) XBIT = -XBIT
      IF (YD.LT.0) YBIT = -YBIT
      IF (XINC-YINC) 27,27,28
 27   YBIT = YBIT+XBIT
      YINC = YINC-XINC
      CYCLE = (XINC+YINC)/XINC
      PMN = Y-YQ
      IF (YD.LT.0.) PMN = -PMN
      PMJ = X-XQ
      IF (XD.LT.0.) PMJ = -PMJ
      T1 = XINC*.5
      T2 = YINC*(.5-PMN+PMJ/CYCLE)
      GO TO 29
 28   XBIT = XBIT+YBIT
      XINC = XINC-YINC
      CYCLE = (XINC+YINC)/YINC
      PMN = X-XQ
      IF (XD.LT.0.) PMN = -PMN
      PMJ = Y-YQ
      IF (YD.LT.0.) PMJ = -PMJ
      T1 = XINC*(.5-PMN+PMJ/CYCLE)
      T2 = YINC*.5
 29   CONTINUE
C  THIS LOOP SCRIBES THE LINE
      DO 10 I=1,NSQT
      IF (T2-T1) 11,11,12
 11   INDEX = INDEX+YBIT
      IF (INDEX.GT.WSUP) INDEX = INDEX-NSAMP4

      IF (HALF)  THEN
         CALL XVTRANS(F2H_BUF,IDN,HDN,1)
         HWORK(INDEX) = HDN
      ELSE
         WORK(INDEX) = INT2BYTE(IDN)
      END IF

      T2 = T2+YINC
      GO TO 10
 12   INDEX = INDEX+XBIT
      IF (INDEX.GT.WSUP) INDEX = INDEX-NSAMP4

      IF (HALF)  THEN
         CALL XVTRANS(F2H_BUF,IDN,HDN,1)
         HWORK(INDEX) = HDN
      ELSE
         WORK(INDEX) = INT2BYTE(IDN)
      END IF

      T1 = T1+XINC
 10   CONTINUE
      ISZ = INDEX/NSAMP4+1
      JSZ = INDEX-(ISZ-1)*NSAMP4
      ISTZ = (IST+ISZ)/2
      JSTZ = (JST+JSZ)/2
      ZINDEX = (ISTZ-1)*NSAMP4+JSTZ

      IF (HALF)  THEN
         CALL XVTRANS(F2H_BUF,IDN,HDN,1)
         CALL XVTRANS(F2H_BUF,IPDN,HPDN,1)
         HWORK(ZINDEX) = HDN
         HWORK(LINDEX) = HDN
         IF (TU.GT..9999) HWORK(LINDEX) = HPDN
      ELSE
         WORK(ZINDEX) = INT2BYTE(IDN)
         WORK(LINDEX) = INT2BYTE(IDN)
         IF (TU.GT..9999) WORK(LINDEX) = INT2BYTE(IPDN)
      END IF

      GO TO 89
 88   CONTIN = .FALSE.
      XFIRST = 0.
 89   CONTINUE
 90   CONTINUE
      end do  ! end of loop thorough all line segments in input file
C
C
C  WRITE OUT THE SCRIBED STRIP.
 91   PTR = 1
      DO 20 I=LL,LU

      IF (HALF)  THEN
	  call xvwrit (out_unit, hwork(ptr), error_num,' ')
      ELSE
	  call xvwrit (out_unit, work(ptr), error_num,' ')
      END IF

 20   PTR = PTR+NSAMP4
      GO TO 100

 401  IF(.NOT.NOPRINT) CALL XVMESSAGE('BAD DATA',' ')

! Close the input graphics file and the output image file
	status= clgr(1)
        if (status.ne.1) call signalgr(1,status,1)
	call xvclose (out_unit, error_num,' ')

	CALL ABEND      ! ABNORMAL END. (NO RETURN FROM ABEND.)

C

8000  CONTINUE

! Close the input graphics file and the output image file
	status= clgr(1)
        if (status.ne.1) call signalgr(1,status,1)
	call xvclose (out_unit, error_num,' ')

      IF(.NOT.NOPRINT) CALL XVMESSAGE('POLYSCRB NORMAL TERMINATION',' ')
      RETURN          ! NORMAL END.
      END

max_stars_repo_name: NASA-AMMOS/VICAR

max_issues_repo_name: NASA-AMMOS/VICAR

max_forks_repo_name: NASA-AMMOS/VICAR


text:
 SUBROUTINE dumpin(ircon,endst,actmsh,endmsh)

! writes restart file

 USE param_db,ONLY: mlin
 USE c_input,ONLY: openfi
 USE ctrl_db
 USE curv_db
 USE damp_db
 USE esets_db,ONLY: eset, msets, esets, nelms, dump_esets
 USE ifx_db
 USE kinc_db, ONLY: dump_kinc, naris, nvelr
 USE mat_dba
 USE nsets_db, ONLY: dump_nsdb
 USE npo_db
 USE outp_db, ONLY: dump_outp
 USE meshmo_db
 USE gvar_db,ONLY: nrst,static
 USE ift_db,ONLY: dump_ift
 USE flc_db,ONLY: dump_flc
 USE cms_db,ONLY: dump_cm
 USE loa_db,ONLY: dump_lo
 USE export_db,ONLY: dump_exp_data
 USE sms_db,ONLY: dump_sms_data
 USE static_db,ONLY: dump_static
 USE nsld_db,ONLY: nnsld

 IMPLICIT NONE

   ! Dummy variables
   INTEGER(kind=4),INTENT(IN OUT):: ircon
   LOGICAL,INTENT(IN):: endst,  & !.TRUE. if the end of the strategy has been reached
                        actmsh, & !.TRUE. if last strategy modified the mesh
                        endmsh    !.TRUE. if end estrategy mesh modification is required
   ! Local variables
   CHARACTER(len=mlin) :: namei
   CHARACTER(len=1),SAVE ::  ni(0:9)=(/'0','1','2','3','4','5','6','7','8','9'/)
   INTEGER (kind=4) :: ircod,j2,j3,j4,flag
   REAL (kind=8) :: cpuac

   INTERFACE
     INCLUDE 'contac.h'
     INCLUDE 'elemnt.h'
   END INTERFACE

   !       ircon defines the file where to WRITE
   ircon = ircon+1
   !       ircod defines the old file to delete (if exist)
   ircod = ircon-nrst  !keeps latest 'nrst' files for restart

   j2 = MOD(ircon,1000)/100
   j3 = MOD(ircon,100)/10
   j4 = MOD(ircon,10)
   namei = '.'//ni(j2)//ni(j3)//ni(j4)   !file name
   IF (endst) THEN  !file corresponds to the end of the strategy
     flag = 0  ! file corresponds to end of the strategy
   ELSE
     flag = 1  ! file corresponds to period NREST
   END IF

   CALL timuse(cpuac)    !compute actual clock time
   CALL openfi(nfile=50,rest=namei,flag=flag,oflag=flag)  !open file to write in

   WRITE(50,ERR=9999) ircon, cpuac
   WRITE(50,ERR=9999) actmsh, endmsh

   CALL dump_ctrl        !dumps CTRL_DB
   CALL dump_esets       !dumps ESETS_DB
   CALL dump_nsdb  ! dumps nodal sets database
   CALL dump_kinc(nnsld)  ! dumps kinematic constrains parameters
   CALL dump_npo   ! dumps nodal data base
   CALL dump_outp  ! dumps output control parameters

   CALL dump_curv    ! dumps curve database
   CALL dump_flc()   ! dumps FLC curve database
   CALL dump_mate    ! dumps material database

   ! dump element databases
   CALL elemnt ('DUMPIN')

   WRITE(50,ERR=9999) 'ENDSET'

   CALL dump_kc (naris, nvelr, ndofn)
   CALL dump_ifx
   CALL dump_cm (ndofn)
   CALL dump_sms_data ()
   CALL dump_lo (ndime,ndofn,nload)
   CALL dump_damp (neq)
   CALL dump_meshmo

   IF(ANY(cactive)) CALL contac('DUMPIN',0)

   IF(therm) THEN
     CALL dump_ift
   END IF
   CALL dump_exp_data
   WRITE(50,ERR=9999) static
   IF( static )CALL dump_static (neq,ndime,npoin)

   CLOSE(50,STATUS='keep')
   CALL flushf(50)

   IF (flag /= -1) THEN
     CALL manage_last_file_names(1)   !Stores the name of last restart file
     IF (ircod > 0 .AND. flag /= -1) THEN
       j2 = MOD(ircod,1000)/100
       j3 = MOD(ircod,100)/10
       j4 = MOD(ircod,10)
       namei = '.'//ni(j2)//ni(j3)//ni(j4)
       CALL del_old_files(flag=3,rest=namei)
     END IF
   END IF

 RETURN
 9999 CALL runen2('')
 END SUBROUTINE dumpin

max_stars_repo_name: jerebenitez/IFE-simpact-openfoam

max_issues_repo_name: jerebenitez/IFE-simpact-openfoam

max_forks_repo_name: jerebenitez/IFE-simpact-openfoam


text:
program while_02
implicit none
integer :: i, j
i = 0
j = 0
do while (i < 10)
    i = i + 1
    j = j + i
end do
if (j /= 55) error stop
if (i /= 10) error stop

i = 0
j = 0
do while (i < 10)
    i = i + 1
    if (i == 2) exit
    j = j + i
end do
if (j /= 1) error stop
if (i /= 2) error stop

i = 0
j = 0
do while (i < 10)
    i = i + 1
    if (i == 2) cycle
    j = j + i
end do
if (j /= 53) error stop
if (i /= 10) error stop
end

max_stars_repo_name: Thirumalai-Shaktivel/lfortran

max_issues_repo_name: Thirumalai-Shaktivel/lfortran

max_forks_repo_name: Thirumalai-Shaktivel/lfortran


text:
      program prog
      dimension nums(2)
      data nums /x'fe', x'dcba'/
      write(6, *) x'A'
      write(6, *) x'AB'
      write(6, *) x'ABC'
      write(6, *) x'ABCD'
      write(6, *) x'7FFFFFFF'
      write(6, *) nums
      end


text:
! 012 LAPACK_BLAS_AUX ##############################################################################################################

      SUBROUTINE  DSCAL(N,DA,DX,INCX)

      USE PENTIUM_II_KIND, ONLY         :  BYTE, LONG, DOUBLE
      USE IOUNT1, ONLY                  :  ERR, F04, F06, SC1, WRT_LOG
      USE SCONTR, ONLY                  :  BLNK_SUB_NAM, FATAL_ERR
      USE TIMDAT, ONLY                  :  HOUR, MINUTE, SEC,
     &                                     SFRAC, TSEC
      USE SUBR_BEGEND_LEVELS, ONLY      :  LAPACK_BEGEND

      USE OUTA_HERE_Interface

c
c     scales a vector by a constant.
c     uses unrolled loops for increment equal to one.
c     jack dongarra, linpack, 3/11/78.
c     modified 3/93 to return if incx .le. 0.
c     modified 12/3/93, array(1) declarations changed to array(*)
c
      REAL(DOUBLE) da,dx(*)
      integer i,incx,m,mp1,n,nincx
c
      if( n.le.0 .or. incx.le.0 )return
      if(incx.eq.1)go to 20
c
c        code for increment not equal to 1
c
      nincx = n*incx
      do 10 i = 1,nincx,incx
        dx(i) = da*dx(i)
   10 continue
      return
c
c        code for increment equal to 1
c
c
c        clean-up loop
c
   20 m = mod(n,5)
      if( m .eq. 0 ) go to 40
      do 30 i = 1,m
        dx(i) = da*dx(i)
   30 continue
      if( n .lt. 5 ) return
   40 mp1 = m + 1
      do 50 i = mp1,n,5
        dx(i) = da*dx(i)
        dx(i + 1) = da*dx(i + 1)
        dx(i + 2) = da*dx(i + 2)
        dx(i + 3) = da*dx(i + 3)
        dx(i + 4) = da*dx(i + 4)
   50 continue
      return

      END SUBROUTINE DSCAL