text:
! MD5 of template: 9d30ae671698cd14111e8b34a690eedf
! ------------------------------------------------------------------------------
! Copyright (c) 2009-13, Thomas P. Robitaille
!
! All rights reserved.
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions are met:
!
!  * Redistributions of source code must retain the above copyright notice, this
!    list of conditions and the following disclaimer.
!
!  * 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.
!
! 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.
! ------------------------------------------------------------------------------

module type_pdf2d

  use lib_array, only : locate, interp2d
  use lib_random, only : random

  implicit none

  private

  integer,parameter :: sp = selected_real_kind(p=6,r=37)
  integer,parameter :: dp = selected_real_kind(p=15,r=307)

  ! We define a 2-d PDF by a set of nx x ny probabilities defined at nodes given by x(1...nx) and y(1...ny). We start by defining a binned PDF that has dimensions nx-1 x ny-1 and defines the probability in the rectangles defined by the nodes. The value of the probability in each bin is simply the average of the four probability points defining the rectangle. Once we have this, we compute the CDF along the x-direction (this also has shape nx-1 x ny-1). We then use the total probabilities in the x-direction to define a PDF in the y direction with ny-1 elements (pdfy). We construct the corresponding CDF (cdfy) and sample a y bin from this. We then sample the x bin from the 2-d CDF along the slice defined by the y bin. Finally, we sample from the function:
  !
  ! z = b1 + b2 * x + b3 * y + b4 * x * y
  !
  ! which is defined by the bilinear interpolation of the four points.
  !
  ! TODO: support for log sampling?

  public :: pdf2d_dp
  type pdf2d_dp
     integer :: nx, ny
     real(dp),allocatable :: x(:), y(:)
     real(dp),allocatable :: prob(:,:)
     real(dp),allocatable :: pdf(:,:)
     real(dp),allocatable :: cdf(:,:)
     real(dp),allocatable :: cdfy(:)
     real(dp),allocatable :: pdfy(:)
     logical :: normalized = .false.
  end type pdf2d_dp

  public :: pdf2d_sp
  type pdf2d_sp
     integer :: nx, ny
     real(sp),allocatable :: x(:), y(:)
     real(sp),allocatable :: prob(:,:)
     real(sp),allocatable :: pdf(:,:)
     real(sp),allocatable :: cdf(:,:)
     real(sp),allocatable :: cdfy(:)
     real(sp),allocatable :: pdfy(:)
     logical :: normalized = .false.
  end type pdf2d_sp

  public :: set_pdf2d
  interface set_pdf2d
     module procedure set_pdf2d_sp
     module procedure set_pdf2d_dp
  end interface set_pdf2d

  public :: sample_pdf2d
  interface sample_pdf2d
     module procedure sample_pdf2d_sp
     module procedure sample_pdf2d_dp
  end interface sample_pdf2d

  public :: interpolate_pdf2d
  interface interpolate_pdf2d
     module procedure interpolate_pdf2d_sp
     module procedure interpolate_pdf2d_dp
  end interface interpolate_pdf2d

contains

  type(pdf2d_dp) function set_pdf2d_dp(x, y, prob) result(p)

    ! Initialize a 2-d PDF object
    !
    ! Parameters
    ! ----------
    ! x : 1-d array (size nx)
    !     The x values at which the probabilites are defined
    ! y : 1-d array (size ny)
    !     The y values at which the probabilites are defined
    ! prob : 2-d array (size nx x ny)
    !     The probabilities defined at (x,y)
    !
    ! Returns
    ! -------
    ! p : pdf2d_dp
    !     The 2-d PDF object

    implicit none

    real(dp),intent(in) :: x(:), y(:), prob(:,:)
    real(dp),allocatable :: area(:,:)

    real(dp) :: norm

    integer :: i, j

    p%nx = size(x)
    p%ny = size(y)

    allocate(p%x(p%nx))
    allocate(p%y(p%ny))

    p%x = x
    p%y = y

    if(size(prob,1) /= size(x)) stop "incorrect dimensions for prob"
    if(size(prob,2) /= size(y)) stop "incorrect dimensions for prob"

    allocate(p%prob(p%nx, p%ny))

    p%prob = prob

    ! Compute area of each rectangle
    allocate(area(p%nx-1, p%ny-1))
    do i=1,p%nx - 1
       do j=1,p%ny - 1
          area(i, j) = (x(i+1) - x(i)) * (y(j+1) - y(j))
       end do
    end do

    ! Compute binned PDF as average of four neighboring points times area of rectangle

    allocate(p%pdf(p%nx-1,p%ny-1))

    p%pdf = (prob(1:p%nx-1,1:p%ny-1) &
         & + prob(1:p%nx-1,2:p%ny) &
         & + prob(2:p%nx,1:p%ny-1) &
         & + prob(2:p%nx,2:p%ny)) &
         & * area * 0.25_dp

    ! Find total probability
    norm = sum(p%pdf)

    ! Normalize unbinned probability
    p%prob = p%prob / norm

    ! Normalize PDF
    p%pdf  = p%pdf / norm
    p%normalized = .true.

    ! Compute 2-d CDF along x direction

    allocate(p%cdf(p%nx-1,p%ny-1))

    p%cdf(1,:) = p%pdf(1,:)
    do i=2,p%nx - 1
       p%cdf(i,:) = p%cdf(i-1,:) + p%pdf(i,:)
    end do

    ! Find PDF and CDF in y direction
    allocate(p%pdfy(p%ny-1))
    allocate(p%cdfy(p%ny-1))
    p%pdfy = p%cdf(p%nx-1,:)
    p%cdfy(1) = p%pdfy(1)
    do j=2,p%ny - 1
       p%cdfy(j) = p%cdfy(j-1) + p%pdfy(j)
    end do
    p%cdfy = p%cdfy / p%cdfy(p%ny-1)

    ! Normalize 2-d CDF (has to be done after calculating the y-axis PDF/CDF)
    do i=1,p%nx - 1
       p%cdf(i,:) = p%cdf(i,:) / p%cdf(p%nx-1,:)
    end do

  end function set_pdf2d_dp

  subroutine sample_pdf2d_dp(p, x, y, xi_alt)

    ! Sample a 2-d PDF
    !
    ! Parameters
    ! ----------
    ! p : pdf2d_dp
    !     The 2-d PDF object to sample from
    ! xi_alt : real(dp) array with 4 elements, optional
    !     Random numbers to use for the sampling
    !
    ! Returns
    ! -------
    ! x, y : real(dp)
    !     The x- and y-position sampled

    implicit none

    type(pdf2d_dp), intent(in) :: p
    real(dp),intent(out) :: x, y
    real(dp),intent(in),optional :: xi_alt(4)
    real(dp) :: xi(4)
    integer :: xbin, ybin
    real(dp) :: b1, b2, b3, b4
    real(dp) :: a, b, c, delta
    integer :: i

    ! Sample random numbers if not specified

    if(present(xi_alt)) then
       xi = xi_alt
    else
       do i=1,4
          call random(xi(i))
       end do
    end if

    ! Find y bin
    if(xi(1) < p%cdfy(1)) then
       ybin = 1
    else
       ybin = locate(p%cdfy, xi(1)) + 1
    end if

    ! Find x bin
    if(xi(2) < p%cdf(1, ybin)) then
       xbin = 1
    else
       xbin = locate(p%cdf(:,ybin), xi(2)) + 1
    end if

    ! Now sample the position within the rectangle. We first find the normalized position in the range [0:1,0:1]. We do this by sampling from a function given by the plane:
    !
    ! z = b1 + b2 * x + b3 * y + b4 * x * y
    !
    ! which is the bilinear interpolation of the points. The coefficients are given by:

    b1 = p%prob(xbin, ybin)
    b2 = p%prob(xbin + 1, ybin) - b1
    b3 = p%prob(xbin, ybin + 1) - b1
    b4 = p%prob(xbin + 1, ybin + 1) - b2 - b3 - b1
    !
    !    print *,'-----'
    !    print *,p%prob(xbin, ybin)
    !    print *,p%prob(xbin+1, ybin)
    !    print *,p%prob(xbin, ybin+1)
    !    print *,p%prob(xbin+1, ybin+1)
    !    print *,b1,b2,b3,b4

    ! We now construct the cumulative PDF in the y direction, and sample from that. The solution is a second-order polynomial with coefficients:

    a = 0.5_dp * (b3 + 0.5_dp * b4)
    b = (b1 + 0.5_dp * b2)
    c = - xi(3) * (a + b)
    delta = sqrt(b * b - 4._dp * a * c)

    ! We have to choose the solution in the range [0:1]

    if (a < 0) then
       if (b > delta) then
          y = (-b + delta) / a * 0.5
       else
          y = (-b - delta) / a * 0.5
       end if
    else if (a > 0) then
       if (-b < delta) then
          y = (-b + delta) / a * 0.5
       else
          y = (-b - delta) / a * 0.5
       end if
    else  ! a == 0 so not a quadratic
       y = - c / b
    end if

    !     print *,a,b,c,y

    ! Now that we have y, we can sample x, and the solution is also a polynomial with coefficients

    a = 0.5_dp * (b2 + b4 * y)
    b = b1 + b3 * y
    c = - xi(4) * (a + b)
    delta = sqrt(b * b - 4._dp * a * c)

    ! We have to choose the solution in the range [0:1]

    if (a < 0) then
       if (b > delta) then
          x = (-b + delta) / a * 0.5
       else
          x = (-b - delta) / a * 0.5
       end if
    else if (a > 0) then
       if (-b < delta) then
          x = (-b + delta) / a * 0.5
       else
          x = (-b - delta) / a * 0.5
       end if
    else  ! a == 0 so not a quadratic
       x = - c / b
    end if

    ! x and y are in relative units - now scale to correct positions
    x = x * (p%x(xbin+1) - p%x(xbin)) + p%x(xbin)
    y = y * (p%y(ybin+1) - p%y(ybin)) + p%y(ybin)

  end subroutine sample_pdf2d_dp

  real(dp) function interpolate_pdf2d_dp(p, x, y, bounds_error, fill_value) result(prob)

    ! Interpolate a 2-d PDF
    !
    ! Parameters
    ! ----------
    ! p : pdf2d_dp
    !     The PDF to interpolate
    ! x, y : real(dp)
    !     Position at which to interpolate the 2-d PDF
    ! bounds_error : logical, optional
    !     Whether to raise an error if the interpolation is out of bounds
    ! fill_value : real(dp)
    !     The value to use for out-of-bounds interpolation if bounds_error = .false.
    !
    ! Returns
    ! -------
    ! prob : real(dp)
    !     The probability at the position requested

    implicit none
    type(pdf2d_dp),intent(in) :: p
    real(dp),intent(in) :: x, y
    logical,intent(in),optional :: bounds_error
    real(dp),intent(in),optional :: fill_value
    if(.not.p%normalized) stop "[interpolate_pdf] PDF is not normalized"
    prob = interp2d(p%x, p%y, p%prob, x, y, bounds_error, fill_value)
  end function interpolate_pdf2d_dp

  type(pdf2d_sp) function set_pdf2d_sp(x, y, prob) result(p)

    ! Initialize a 2-d PDF object
    !
    ! Parameters
    ! ----------
    ! x : 1-d array (size nx)
    !     The x values at which the probabilites are defined
    ! y : 1-d array (size ny)
    !     The y values at which the probabilites are defined
    ! prob : 2-d array (size nx x ny)
    !     The probabilities defined at (x,y)
    !
    ! Returns
    ! -------
    ! p : pdf2d_sp
    !     The 2-d PDF object

    implicit none

    real(sp),intent(in) :: x(:), y(:), prob(:,:)
    real(sp),allocatable :: area(:,:)

    real(sp) :: norm

    integer :: i, j

    p%nx = size(x)
    p%ny = size(y)

    allocate(p%x(p%nx))
    allocate(p%y(p%ny))

    p%x = x
    p%y = y

    if(size(prob,1) /= size(x)) stop "incorrect dimensions for prob"
    if(size(prob,2) /= size(y)) stop "incorrect dimensions for prob"

    allocate(p%prob(p%nx, p%ny))

    p%prob = prob

    ! Compute area of each rectangle
    allocate(area(p%nx-1, p%ny-1))
    do i=1,p%nx - 1
       do j=1,p%ny - 1
          area(i, j) = (x(i+1) - x(i)) * (y(j+1) - y(j))
       end do
    end do

    ! Compute binned PDF as average of four neighboring points times area of rectangle

    allocate(p%pdf(p%nx-1,p%ny-1))

    p%pdf = (prob(1:p%nx-1,1:p%ny-1) &
         & + prob(1:p%nx-1,2:p%ny) &
         & + prob(2:p%nx,1:p%ny-1) &
         & + prob(2:p%nx,2:p%ny)) &
         & * area * 0.25_dp

    ! Find total probability
    norm = sum(p%pdf)

    ! Normalize unbinned probability
    p%prob = p%prob / norm

    ! Normalize PDF
    p%pdf  = p%pdf / norm
    p%normalized = .true.

    ! Compute 2-d CDF along x direction

    allocate(p%cdf(p%nx-1,p%ny-1))

    p%cdf(1,:) = p%pdf(1,:)
    do i=2,p%nx - 1
       p%cdf(i,:) = p%cdf(i-1,:) + p%pdf(i,:)
    end do

    ! Find PDF and CDF in y direction
    allocate(p%pdfy(p%ny-1))
    allocate(p%cdfy(p%ny-1))
    p%pdfy = p%cdf(p%nx-1,:)
    p%cdfy(1) = p%pdfy(1)
    do j=2,p%ny - 1
       p%cdfy(j) = p%cdfy(j-1) + p%pdfy(j)
    end do
    p%cdfy = p%cdfy / p%cdfy(p%ny-1)

    ! Normalize 2-d CDF (has to be done after calculating the y-axis PDF/CDF)
    do i=1,p%nx - 1
       p%cdf(i,:) = p%cdf(i,:) / p%cdf(p%nx-1,:)
    end do

  end function set_pdf2d_sp

  subroutine sample_pdf2d_sp(p, x, y, xi_alt)

    ! Sample a 2-d PDF
    !
    ! Parameters
    ! ----------
    ! p : pdf2d_sp
    !     The 2-d PDF object to sample from
    ! xi_alt : real(sp) array with 4 elements, optional
    !     Random numbers to use for the sampling
    !
    ! Returns
    ! -------
    ! x, y : real(sp)
    !     The x- and y-position sampled

    implicit none

    type(pdf2d_sp), intent(in) :: p
    real(sp),intent(out) :: x, y
    real(sp),intent(in),optional :: xi_alt(4)
    real(sp) :: xi(4)
    integer :: xbin, ybin
    real(sp) :: b1, b2, b3, b4
    real(sp) :: a, b, c, delta
    integer :: i

    ! Sample random numbers if not specified

    if(present(xi_alt)) then
       xi = xi_alt
    else
       do i=1,4
          call random(xi(i))
       end do
    end if

    ! Find y bin
    if(xi(1) < p%cdfy(1)) then
       ybin = 1
    else
       ybin = locate(p%cdfy, xi(1)) + 1
    end if

    ! Find x bin
    if(xi(2) < p%cdf(1, ybin)) then
       xbin = 1
    else
       xbin = locate(p%cdf(:,ybin), xi(2)) + 1
    end if

    ! Now sample the position within the rectangle. We first find the normalized position in the range [0:1,0:1]. We do this by sampling from a function given by the plane:
    !
    ! z = b1 + b2 * x + b3 * y + b4 * x * y
    !
    ! which is the bilinear interpolation of the points. The coefficients are given by:

    b1 = p%prob(xbin, ybin)
    b2 = p%prob(xbin + 1, ybin) - b1
    b3 = p%prob(xbin, ybin + 1) - b1
    b4 = p%prob(xbin + 1, ybin + 1) - b2 - b3 - b1
    !
    !    print *,'-----'
    !    print *,p%prob(xbin, ybin)
    !    print *,p%prob(xbin+1, ybin)
    !    print *,p%prob(xbin, ybin+1)
    !    print *,p%prob(xbin+1, ybin+1)
    !    print *,b1,b2,b3,b4

    ! We now construct the cumulative PDF in the y direction, and sample from that. The solution is a second-order polynomial with coefficients:

    a = 0.5_sp * (b3 + 0.5_sp * b4)
    b = (b1 + 0.5_sp * b2)
    c = - xi(3) * (a + b)
    delta = sqrt(b * b - 4._sp * a * c)

    ! We have to choose the solution in the range [0:1]

    if (a < 0) then
       if (b > delta) then
          y = (-b + delta) / a * 0.5
       else
          y = (-b - delta) / a * 0.5
       end if
    else if (a > 0) then
       if (-b < delta) then
          y = (-b + delta) / a * 0.5
       else
          y = (-b - delta) / a * 0.5
       end if
    else  ! a == 0 so not a quadratic
       y = - c / b
    end if

    !     print *,a,b,c,y

    ! Now that we have y, we can sample x, and the solution is also a polynomial with coefficients

    a = 0.5_sp * (b2 + b4 * y)
    b = b1 + b3 * y
    c = - xi(4) * (a + b)
    delta = sqrt(b * b - 4._sp * a * c)

    ! We have to choose the solution in the range [0:1]

    if (a < 0) then
       if (b > delta) then
          x = (-b + delta) / a * 0.5
       else
          x = (-b - delta) / a * 0.5
       end if
    else if (a > 0) then
       if (-b < delta) then
          x = (-b + delta) / a * 0.5
       else
          x = (-b - delta) / a * 0.5
       end if
    else  ! a == 0 so not a quadratic
       x = - c / b
    end if

    ! x and y are in relative units - now scale to correct positions
    x = x * (p%x(xbin+1) - p%x(xbin)) + p%x(xbin)
    y = y * (p%y(ybin+1) - p%y(ybin)) + p%y(ybin)

  end subroutine sample_pdf2d_sp

  real(sp) function interpolate_pdf2d_sp(p, x, y, bounds_error, fill_value) result(prob)

    ! Interpolate a 2-d PDF
    !
    ! Parameters
    ! ----------
    ! p : pdf2d_sp
    !     The PDF to interpolate
    ! x, y : real(sp)
    !     Position at which to interpolate the 2-d PDF
    ! bounds_error : logical, optional
    !     Whether to raise an error if the interpolation is out of bounds
    ! fill_value : real(sp)
    !     The value to use for out-of-bounds interpolation if bounds_error = .false.
    !
    ! Returns
    ! -------
    ! prob : real(sp)
    !     The probability at the position requested

    implicit none
    type(pdf2d_sp),intent(in) :: p
    real(sp),intent(in) :: x, y
    logical,intent(in),optional :: bounds_error
    real(sp),intent(in),optional :: fill_value
    if(.not.p%normalized) stop "[interpolate_pdf] PDF is not normalized"
    prob = interp2d(p%x, p%y, p%prob, x, y, bounds_error, fill_value)
  end function interpolate_pdf2d_sp

end module type_pdf2d

max_stars_repo_name: jfitz/code-stat

max_issues_repo_name: jfitz/code-stat

max_forks_repo_name: jfitz/code-stat


text:
      subroutine repl(x,y,n,ndim)
c====================
c
c    x=y
c
      complex*16 x(ndim,ndim),y(ndim,ndim)
c
      do i=1,n
      do j=1,n
        x(i,j)=y(i,j)
      end do
      end do
      return
      end
      subroutine replt(x,y,n,ndim)
c====================
c
c    x=y^T
c
      complex*16 x(ndim,ndim),y(ndim,ndim)
c
      do i=1,n
      do j=1,n
        x(i,j)=y(j,i)
      end do
      end do
      return
      end
      subroutine replrel(x,y,n,ndim)
c=======================
c
c    x(k,k')=(-)**(l-l')* y(k',k)*
c
      complex*16 x(ndim,ndim),y(ndim,ndim)
      dimension ldex(50)
c
      data ldex/0,0,
     *          1,1,1,1,1,1,
     *          2,2,2,2,2,2,2,2,2,2,
     *          3,3,3,3,3,3,3,3,3,3,3,3,3,3,
     *          4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4/
c
      do i=1,n
        li=ldex(i)
      do j=1,n
        lj=ldex(j)
        x(i,j)=dconjg(y(j,i))
        if(mod(iabs(li-lj),2).eq.1) x(i,j)=-x(i,j)
      end do
      end do
      return
      end
      subroutine compmat(x,y,n,ndim,tol,ic)
c=======================
c
c check whether  x=y
c
      implicit real*8 (a-h,o-z)
      complex*16 x(ndim,ndim),y(ndim,ndim),xx,dx
      data tol0/1.0d-12/
c
      if(tol.le.0.d0) tol=tol0
      ic=0
      xnorm=0.d0
      dxnorm=0.d0
      do i=1,n
      do j=1,n
        xx=x(i,j)
        dx=x(i,j)-y(i,j)
        xxre=dreal(xx)
        xxim=dimag(xx)
        dxre=dreal(dx)
        dxim=dimag(dx)
        xnorm=xnorm+xxre*xxre+xxim*xxim
        dxnorm=dxnorm+dxre*dxre+dxim*dxim
      end do
      end do
      xnorm=dsqrt(xnorm)
      dxnorm=dsqrt(dxnorm)
      if(dxnorm.lt.tol*xnorm) ic=1
      return
      end
      subroutine addmat(x,y,n,ndim)
c====================
c
c    x=x+y
c
      complex*16 x(ndim,ndim),y(ndim,ndim)
c
      do i=1,n
      do j=1,n
        x(i,j)=x(i,j)+y(i,j)
      end do
      end do
      return
      end
      subroutine addmat1(a,b,c,n,ndim)
c====================
c
c    c=a+b
c
      complex*16 a(ndim,ndim),b(ndim,ndim),c(ndim,ndim)
c
      do i=1,n
      do j=1,n
        c(i,j)=a(i,j)+b(i,j)
      end do
      end do
      return
      end
      subroutine submat(x,y,n,ndim)
c====================
c
c    x=x-y
c
      complex*16 x(ndim,ndim),y(ndim,ndim)
c
      do i=1,n
      do j=1,n
        x(i,j)=x(i,j)-y(i,j)
      end do
      end do
      return
      end
      subroutine submat1(a,b,c,n,ndim)
c====================
c
c    c=a-b
c
      complex*16 a(ndim,ndim),b(ndim,ndim),c(ndim,ndim)
c
      do i=1,n
      do j=1,n
        c(i,j)=a(i,j)-b(i,j)
      end do
      end do
      return
      end
      subroutine symmat(x,n,ndim)
c======================
c Symmetrize matrix x:
c    x -> 1/2*(x+xT)
c
      complex*16 x(ndim,ndim)
c
      do i=1,n
      do j=i+1,n
        x(i,j)=0.5d0*(x(i,j)+x(j,i))
        x(j,i)=x(i,j)
      end do
      end do
      return
      end
      subroutine doubmt(amat,bmat,ndim,ndimp)
c======================
c
c  amat=amat*bmat
c
      implicit real*8 (a-h,o-z)
c
      complex*16 amat(ndimp,ndimp),bmat(ndimp,ndimp)
      complex*16 help(100),sum
c
      do i=1,ndim
        do k=1,ndim
          sum=dcmplx(0.d0,0.d0)
          do j=1,ndim
            sum=sum+bmat(j,k)*amat(i,j)
          end do
          help(k)=sum
        end do
        do k=1,ndim
          amat(i,k)=help(k)
        end do
      end do
c
      return
      end
      subroutine doubmt1(amat,bmat,cmat,ndim,ndimp)
c======================
c
c  cmat=amat*bmat
c
      implicit real*8 (a-h,o-z)
c
      complex*16 amat(ndimp,ndimp),bmat(ndimp,ndimp),cmat(ndimp,ndimp)
      complex*16 sum
c
      do i=1,ndim
        do k=1,ndim
          sum=dcmplx(0.d0,0.d0)
          do j=1,ndim
            sum=sum+bmat(j,k)*amat(i,j)
          end do
          cmat(i,k)=sum
        end do
      end do
c
      return
      end
      function trdbmt(amat,bmat,ndim,ndimp)
c===============================
c
c  trdbmt=Tr(amat*bmat)
c
      integer i,j,ndim,ndimp
      complex*16 trdbmt,amat(ndimp,ndimp),bmat(ndimp,ndimp)
c
      trdbmt=dcmplx(0.d0,0.d0)
      do i=1,ndim
      do j=1,ndim
        trdbmt=trdbmt+bmat(j,i)*amat(i,j)
      end do
      end do
c
      return
      end
      subroutine tripmt(u,b,ust,ndi1,ndi2,ndim)
c =====================
c
c vectorized routine for triple product of rectangular matrices
c
      implicit real*8 (a-h,o-z)
      parameter(nndim=100)
      complex*16 u,ust,b,c,x
      dimension u(ndim,ndim),ust(ndim,ndim),b(ndim,ndim)
      dimension c(nndim,nndim)
c
c     left product
c
      do i=1,ndi1
      do j=1,ndi2
c
      x= (0.0d0,0.0d0)
c
      do k=1,ndi1
        x = x + b(k,j)*u(i,k)
      end do
c
      c(i,j)=x
      end do
      end do
c
c     right product
c
      do i=1,ndi1
      do j=1,ndi2
c
      x = (0.d0,0.0d0)
c
      do k=1,ndi2
        x = x + ust(k,j)*c(i,k)
      end do
c
      b(i,j)=x
      end do
      end do
c
      return
      end
      subroutine tripmt1(u,b,ust,b1,ndi1,ndi2,ndim)
c =====================
c
c vectorized routine for triple product of rectangular matrices
c
      implicit real*8 (a-h,o-z)
      parameter(nndim=100)
      complex*16 u,ust,b,b1,c,x
      dimension u(ndim,ndim),ust(ndim,ndim)
      dimension b(ndim,ndim),b1(ndim,ndim)
      dimension c(nndim,nndim)
c
c     left product
c
      do i=1,ndi1
      do j=1,ndi2
c
      x= (0.0d0,0.0d0)
c
      do k=1,ndi1
        x = x + b(k,j)*u(i,k)
      end do
c
      c(i,j)=x
      end do
      end do
c
c     right product
c
      do i=1,ndi1
      do j=1,ndi2
c
      x = (0.d0,0.0d0)
c
      do k=1,ndi2
        x = x + ust(k,j)*c(i,k)
      end do
c
      b1(i,j)=x
      end do
      end do
c
      return
      end
      subroutine outmat(mat,n,m,ndim,nper)
c=====================
      implicit real*8 (a-h,o-z)
      complex*16 mat(ndim,ndim)
      complex*16 mat1(50,50)
c
      do i=1,n
       do j=1,m
        r1=dreal(mat(i,j))
        r2=dimag(mat(i,j))
        if(dabs(r1).lt.1.d-15) r1=0.d0
        if(dabs(r2).lt.1.d-15) r2=0.d0
        mat1(i,j)=dcmplx(r1,r2)
       end do
      end do
c
      write(nper,*) ' real part'
      do i=1,n
        write(nper,10) (dreal(mat1(i,j)),j=1,m)
      end do
      write(nper,*) ' imaginary part'
      do i=1,n
        write(nper,10) (dimag(mat1(i,j)),j=1,m)
      end do
c
  10  format(9(1pd14.6))
c
      return
      end
      subroutine outmat1(mat,n,m,ndim,tol,nper)
c=======================
      implicit real*8 (a-h,o-z)
      complex*16 mat(ndim,ndim)
c
      do j=1,m
      do i=1,n
        r1=dreal(mat(i,j))
        r2=dimag(mat(i,j))
        if(dabs(r1).ge.tol.or.dabs(r2).ge.tol) then
        write(nper,'(2i4,5x,1pd20.10,1pd20.10)') i,j,r1,r2
        end if
      end do
      end do
c
      return
      end

ext: f


MODULE nrtype
	INTEGER, PARAMETER :: I4B = SELECTED_INT_KIND(9)
	INTEGER, PARAMETER :: I2B = SELECTED_INT_KIND(4)
	INTEGER, PARAMETER :: I1B = SELECTED_INT_KIND(2)
	INTEGER, PARAMETER :: SP = KIND(1.0)
	INTEGER, PARAMETER :: DP = KIND(1.0D0)
	INTEGER, PARAMETER :: SPC = KIND((1.0,1.0))
	INTEGER, PARAMETER :: DPC = KIND((1.0D0,1.0D0))
	INTEGER, PARAMETER :: LGT = KIND(.true.)
	REAL(SP), PARAMETER :: PI=3.141592653589793238462643383279502884197_sp
	REAL(SP), PARAMETER :: PIO2=1.57079632679489661923132169163975144209858_sp
	REAL(SP), PARAMETER :: TWOPI=6.283185307179586476925286766559005768394_sp
	REAL(SP), PARAMETER :: SQRT2=1.41421356237309504880168872420969807856967_sp
	REAL(SP), PARAMETER :: EULER=0.5772156649015328606065120900824024310422_sp
	REAL(DP), PARAMETER :: PI_D=3.141592653589793238462643383279502884197_dp
	REAL(DP), PARAMETER :: PIO2_D=1.57079632679489661923132169163975144209858_dp
	REAL(DP), PARAMETER :: TWOPI_D=6.283185307179586476925286766559005768394_dp
	TYPE sprs2_sp
		INTEGER(I4B) :: n,len
		REAL(SP), DIMENSION(:), POINTER :: val
		INTEGER(I4B), DIMENSION(:), POINTER :: irow
		INTEGER(I4B), DIMENSION(:), POINTER :: jcol
	END TYPE sprs2_sp
	TYPE sprs2_dp
		INTEGER(I4B) :: n,len
		REAL(DP), DIMENSION(:), POINTER :: val
		INTEGER(I4B), DIMENSION(:), POINTER :: irow
		INTEGER(I4B), DIMENSION(:), POINTER :: jcol
	END TYPE sprs2_dp
END MODULE nrtype

module dimensions
use nrtype
    implicit none
    !nkk: grid point cash on hand (coh) grid
    !f_t:family type
    !h_t:health type
    !g_t: gender type
    !PI_t: permanent income type
    integer,parameter::f_t=2,clusters=4,L_gender=2,L_PI=5,generations=21,nzz=3,variables=12,groups=4,parameters_to_est=11,obs=9, &
                        wealth_q=3,L_PI2=10,min_obs=39
    integer,parameter:: moment_conditions=L_PI*obs*groups+L_PI*f_t+f_t*obs*groups+L_PI*clusters+2*f_t*clusters
    integer,parameter::samples_per_i=300,nkk=500
end module dimensions

module grids
    use dimensions; use nrtype
    implicit none
    real(SP)::coh_min=0.0_sp,coh_max=1200.0_sp
    real(SP),dimension(nkk)::coh_grid
    real(SP),dimension(nzz,1)::ep_grid
    real(SP),dimension(nzz,1)::xi_grid
    real(SP),dimension(nzz,1)::pr_ep,pr0_xi
    real(SP),dimension(nzz,nzz)::pr_pxi
end module grids

module targets
    use dimensions; use nrtype
    implicit none
    real(SP),dimension(L_PI,obs,groups)::data_NW_PI,data_NW_PI1,data_NW_PIb,data_NW_PI1b
    real(SP),dimension(2,obs)::data_NW_h_ut
    real(SP),dimension(L_PI,f_t)::data_beq100_IC
    real(SP),dimension(L_PI,clusters)::data_lfc_PI
    real(SP),dimension(f_t,clusters)::data_lfc_IC,data_govmd_IC
    real(SP),dimension(f_t,obs,groups)::data_NW_IC,data_NW_IC1
    real(SP),dimension(moment_conditions,1)::data_moments_new,data_moments1,data_moments
end module targets

module structural_p2
    use dimensions; use nrtype
    implicit none
    real(SP):: sigma,sigma_beq,nu,beta,omega,share_p,subs_p
    real(SP),dimension(clusters):: kappa_h,delta_h,u_bar_no_f,x_bar
    real(SP),dimension(clusters,f_t)::u_bar
    real(SP),dimension(clusters-1)::alpha_mu
    real(SP),dimension(clusters)::c_bar,l_bar
    real(SP),dimension(f_t):: lambda,delta
    real(SP),dimension(clusters):: mu
    real(SP),dimension(moment_conditions,moment_conditions):: W_opt,Phi
    real(SP),dimension(parameters_to_est)::se
    end module structural_p2

module structural_p1
    use dimensions; use nrtype
    implicit none
    integer,parameter::DIM=12+11+3
    real(SP),dimension(DIM,1)::med_coef
    real(SP)::rho,sigma2_ep,sigma2_ze
    real(SP),parameter::r=1.02_sp**2_sp-1.0_sp
    real(SP)::p_fc=12.0_sp/1000.0_sp
    real(SP),dimension(clusters+1)::p_nh=(/8.0_sp/365.0_sp*443.0_sp,8.0_sp/365.0_sp*437.0_sp,8.0/365.0_sp*463.0_sp,8.0/365.0_sp*527.0_sp,8.0/365.0_sp*527.0_sp/2.0_sp/)
    real(SP),parameter::p_or=12.0_sp/1000.0_sp
    real(SP),dimension(L_PI2,L_gender)::b
    real(SP),dimension(f_t,clusters)::l_ic,l_ic_or
    real(SP),dimension(f_t,clusters+1,2)::pr_nh
    real(SP),dimension(generations,L_gender,L_PI2,clusters,nzz,nzz)::m_exp_all,m_exp_all_or
    real(SP),dimension(clusters+1,clusters+1,generations,L_PI2,L_gender)::H_av
    real(SP)::load_ltci=0.32_sp,benefit_LTCI=5.0_sp !hours of care per day provided by insurance
    real(SP),dimension(L_PI,2,clusters)::price_ltci

end module structural_p1

module files_savings
implicit none
    character(LEN=76)::fdir_med="C:\Users\jbueren\Google Drive\JMP\Code\Health Dynamics\medical_expenditures\"
    character(LEN=64)::fdir_inputs="C:\Users\jbueren\Google Drive\JMP\Code\Structural Model\moments\"
    character(LEN=63)::path_h="C:\Users\jbueren\Google Drive\JMP\Code\Health Dynamics\Results\"
    !character(LEN=74)::fdir_med="C:\Users\Jesus\Google Drive\JMP\Code\Health Dynamics\medical_expenditures\"
    !character(LEN=62)::fdir_inputs="C:\Users\Jesus\Google Drive\JMP\Code\Structural Model\moments\"
    end module

module simulation_input
use nrtype; use dimensions
implicit none
    integer,parameter::indv=8472
    integer,dimension(indv):: gender_i, generation_i, PI_q_i, PI_q_i2, group_i,iwendy_i
    real(SP),dimension(indv,f_t):: IC_pr_i
    real(SP),dimension(indv):: x_i
    real(SP),dimension(indv,obs,clusters):: s_h_i,f_h_i
    real(SP),dimension(indv,obs):: dead_i,nhmliv_i
    integer,dimension(indv,obs,2):: h_i
    integer::ind_no_f=0
end module

module optimization
use nrtype; use structural_p2; use dimensions
implicit none
    real(SP)::obj_fct_s
    real(SP),dimension(moment_conditions,1)::r1,f
    integer::real_moments
    real(SP),dimension(moment_conditions,moment_conditions)::matrix_s
    real(SP),dimension(L_PI,obs,groups)::moments_NW_PI1
    real(SP),dimension(f_t,obs,groups)::moments_NW_IC1
    real(SP),dimension(wealth_q,obs,groups)::moments_NW_all1
    real(SP):: obj_fct_min=1.0_sp/0.0_sp
end module

module pdfs
    use nrtype;use dimensions;use simulation_input
    implicit none
    real(SP),dimension(L_PI,obs,groups,nkk)::pdf_all_NW_PI
end module

module MD_reform
use dimensions;use nrtype
implicit none
real(SP),dimension(nkk,clusters+1,nzz,L_gender,L_PI2,f_t)::V_70_or,V_70_new,V_70
integer::ind_or=1
real(SP)::p_sub=0.0_sp
end module

module HRS_data
    use dimensions; use nrtype; use simulation_input
    implicit none
    integer,dimension(indv,obs):: IC_q,f1nhmliv
    real(SP),dimension(indv,obs):: fc_h,NW,ic_h,beq100,govmd
end module HRS_data


SUBROUTINE RelaxationModule(E,P,time)
!
! Apply relaxation zone
!
! By Allan P. Engsig-Karup.
USE GlobalVariables, ONLY: RelaxZones, FineGrid, SFsol, g, relaxTransientTime, relaxNo, &
	GhostGridX, GhostGridY, curvilinearONOFF
USE Precision
USE Constants
IMPLICIT NONE
!REAL(KIND=long), DIMENSION(FineGrid%Nx,FineGrid%Ny) :: E, P
! GD: change
REAL(KIND=long), DIMENSION(FineGrid%Nx+2*GhostGridX,FineGrid%Ny+2*GhostGridY) :: E, P
REAL(KIND=long) :: time, FAC
INTEGER :: i, j, k
REAL(KIND=long) :: tmpx(FineGrid%Nx+2*GhostGridX), tmpy(FineGrid%Ny+2*GhostGridY)

IF (time<relaxTransientTime) THEN
  FAC = time/relaxTransientTime
ELSE
  FAC = one
ENDIF
IF (curvilinearONOFF==1) THEN ! curvilinear coordinates... FIXME for make it work whatever transformation
  ! FIXME: use Lx and Ly? Depends how x and y are defined... save e and n directly ?
  DO j=1,FineGrid%Nx+2*GhostGridX
    !tmpx(j) = REAL(j-2,long)*(FineGrid%y(2,1)-FineGrid%y(1,1)) !/REAL((FineGrid%Nx+2*GhostGridX-1),long)*(8+PI)*Lx; !FIXME find a way to define in all situation...
    tmpx(j) = FineGrid%y(j,1)
    tmpx(j) = REAL(j,long)
  ENDDO
  DO j=1,FineGrid%Ny+2*GhostGridY
    !tmpy(j) = -REAL(j-2,long)*(FineGrid%x(1,2)-FineGrid%x(1,1)) !/REAL((FineGrid%Ny+2*GhostGridY-1),long)*(Ly-Lx); !FIXME find a way to define in all situation...
    tmpy(j) = FineGrid%x(1,j)
    tmpy(j) = REAL(j,long)
  ENDDO
ELSE
  tmpx(:) = FineGrid%x(:,1)
  tmpy(:) = FineGrid%y(1,:)
ENDIF
! FOR EACH RELAXATION ZONE
DO i = 1, relaxNo
  IF (RelaxZones(i)%XorY=='X' .AND. RelaxZones(i)%XorYgen=='X') THEN
	 IF (RelaxZones(i)%WavegenONOFF==0) THEN
	   DO j = RelaxZones(i)%idx(3),RelaxZones(i)%idx(4)
         E(RelaxZones(i)%idx(1):RelaxZones(i)%idx(2),j) = E(RelaxZones(i)%idx(1):RelaxZones(i)%idx(2),j)*(RelaxZones(i)%gam)
         P(RelaxZones(i)%idx(1):RelaxZones(i)%idx(2),j) = P(RelaxZones(i)%idx(1):RelaxZones(i)%idx(2),j)*(RelaxZones(i)%gam)
 	   END DO
	 ELSE
	  IF (RelaxZones(i)%degrees==zero) THEN
      	CALL stream_func_wave_finite(RelaxZones(i)%idx(2)-RelaxZones(i)%idx(1)+1,&
		     FineGrid%x(RelaxZones(i)%idx(1):RelaxZones(i)%idx(2),1),time &
			 ,SFsol%n_four_modes,SFsol%zz,SFsol%yy,SFsol%k,g,RelaxZones(i)%Ea,RelaxZones(i)%Pa)
        !CALL stream_func_wave_finite(RelaxZones(i)%idx(2)-RelaxZones(i)%idx(1)+1,&
		!     tmpx(RelaxZones(i)%idx(1):RelaxZones(i)%idx(2)),time &
		!	 ,SFsol%n_four_modes,SFsol%zz,SFsol%yy,SFsol%k,g,RelaxZones(i)%Ea,RelaxZones(i)%Pa)
	    DO j = RelaxZones(i)%idx(3),RelaxZones(i)%idx(4)
          E(RelaxZones(i)%idx(1):RelaxZones(i)%idx(2),j) = E(RelaxZones(i)%idx(1):RelaxZones(i)%idx(2),j)*(RelaxZones(i)%gam) &
		      + FAC*RelaxZones(i)%Ea*(one-RelaxZones(i)%gam)
          P(RelaxZones(i)%idx(1):RelaxZones(i)%idx(2),j) = P(RelaxZones(i)%idx(1):RelaxZones(i)%idx(2),j)*(RelaxZones(i)%gam) &
		      + FAC*RelaxZones(i)%Pa*(one-RelaxZones(i)%gam)
 	    END DO
	  ELSE
	    DO j = RelaxZones(i)%idx(3),RelaxZones(i)%idx(4)
          CALL stream_func_wave_finite(RelaxZones(i)%idx(2)-RelaxZones(i)%idx(1)+1,&
		       tmpx(RelaxZones(i)%idx(1):RelaxZones(i)%idx(2))&
			   *COS(RelaxZones(i)%degrees/180.0_long*pi)+tmpy(j)*SIN(RelaxZones(i)%degrees/180.0_long*pi),&
			   time,SFsol%n_four_modes,SFsol%zz,SFsol%yy,SFsol%k,g,RelaxZones(i)%Ea,RelaxZones(i)%Pa)
          E(RelaxZones(i)%idx(1):RelaxZones(i)%idx(2),j) = E(RelaxZones(i)%idx(1):RelaxZones(i)%idx(2),j)*(RelaxZones(i)%gam) &
		       + FAC*RelaxZones(i)%Ea*(one-RelaxZones(i)%gam)
          P(RelaxZones(i)%idx(1):RelaxZones(i)%idx(2),j) = P(RelaxZones(i)%idx(1):RelaxZones(i)%idx(2),j)*(RelaxZones(i)%gam) &
		       + FAC*RelaxZones(i)%Pa*(one-RelaxZones(i)%gam)
 	    END DO
	  ENDIF
	 ENDIF
  ELSE IF (RelaxZones(i)%XorY=='Y' .AND. RelaxZones(i)%XorYgen=='Y') THEN ! Y-direction
	 IF (RelaxZones(i)%WavegenONOFF==0) THEN
	   DO j = RelaxZones(i)%idx(1),RelaxZones(i)%idx(2)
         E(j,RelaxZones(i)%idx(3):RelaxZones(i)%idx(4)) = E(j,RelaxZones(i)%idx(3):RelaxZones(i)%idx(4))*(RelaxZones(i)%gam)
         P(j,RelaxZones(i)%idx(3):RelaxZones(i)%idx(4)) = P(j,RelaxZones(i)%idx(3):RelaxZones(i)%idx(4))*(RelaxZones(i)%gam)
 	   END DO
	 ELSE
	   IF (RelaxZones(i)%degrees==zero) THEN
         CALL stream_func_wave_finite(RelaxZones(i)%idx(4)-RelaxZones(i)%idx(3)+1,&
		      tmpy(RelaxZones(i)%idx(3):RelaxZones(i)%idx(4)),time,&
			  SFsol%n_four_modes,SFsol%zz,SFsol%yy,SFsol%k,g,RelaxZones(i)%Ea,RelaxZones(i)%Pa)
	     DO j = RelaxZones(i)%idx(1),RelaxZones(i)%idx(2)
           E(j,RelaxZones(i)%idx(3):RelaxZones(i)%idx(4)) = E(j,RelaxZones(i)%idx(3):RelaxZones(i)%idx(4))*(RelaxZones(i)%gam) &
		        + FAC*RelaxZones(i)%Ea*(one-RelaxZones(i)%gam)
           P(j,RelaxZones(i)%idx(3):RelaxZones(i)%idx(4)) = P(j,RelaxZones(i)%idx(3):RelaxZones(i)%idx(4))*(RelaxZones(i)%gam) &
		        + FAC*RelaxZones(i)%Pa*(one-RelaxZones(i)%gam)
 	     END DO
	   ELSE
	     DO j = RelaxZones(i)%idx(1),RelaxZones(i)%idx(2)
           CALL stream_func_wave_finite(RelaxZones(i)%idx(4)-RelaxZones(i)%idx(3)+1,&
		        -tmpy(RelaxZones(i)%idx(3):RelaxZones(i)%idx(4))&
				*SIN(RelaxZones(i)%degrees/180.0_long*pi)+tmpx(j)*COS(RelaxZones(i)%degrees/180.0_long*pi),time,&
				SFsol%n_four_modes,SFsol%zz,SFsol%yy,SFsol%k,g,RelaxZones(i)%Ea,RelaxZones(i)%Pa)
           E(j,RelaxZones(i)%idx(3):RelaxZones(i)%idx(4)) = E(j,RelaxZones(i)%idx(3):RelaxZones(i)%idx(4))*(RelaxZones(i)%gam) &
		        + FAC*RelaxZones(i)%Ea*(one-RelaxZones(i)%gam)
           P(j,RelaxZones(i)%idx(3):RelaxZones(i)%idx(4)) = P(j,RelaxZones(i)%idx(3):RelaxZones(i)%idx(4))*(RelaxZones(i)%gam) &
		        + FAC*RelaxZones(i)%Pa*(one-RelaxZones(i)%gam)
 	     END DO
	   ENDIF
	 ENDIF
  ELSE IF (RelaxZones(i)%XorY=='X' .AND. RelaxZones(i)%XorYgen=='Y') THEN
	 IF (RelaxZones(i)%WavegenONOFF==0) THEN
	   DO j = RelaxZones(i)%idx(3),RelaxZones(i)%idx(4)
         E(RelaxZones(i)%idx(1):RelaxZones(i)%idx(2),j) = E(RelaxZones(i)%idx(1):RelaxZones(i)%idx(2),j)*(RelaxZones(i)%gam)
         P(RelaxZones(i)%idx(1):RelaxZones(i)%idx(2),j) = P(RelaxZones(i)%idx(1):RelaxZones(i)%idx(2),j)*(RelaxZones(i)%gam)
 	   END DO
	 ELSE
	  IF (RelaxZones(i)%degrees==zero) THEN
         CALL stream_func_wave_finite(RelaxZones(i)%idx(4)-RelaxZones(i)%idx(3)+1,&
			tmpy(RelaxZones(i)%idx(3):RelaxZones(i)%idx(4)),time,&
			SFsol%n_four_modes,SFsol%zz,SFsol%yy,SFsol%k,g,RelaxZones(i)%Ea,RelaxZones(i)%Pa)
	     DO j = RelaxZones(i)%idx(1),RelaxZones(i)%idx(2)
		   k = j - RelaxZones(i)%idx(1) + 1
           E(j,RelaxZones(i)%idx(3):RelaxZones(i)%idx(4)) = E(j,RelaxZones(i)%idx(3):RelaxZones(i)%idx(4))*(RelaxZones(i)%gam) &
		       + FAC*RelaxZones(i)%Ea*(one-RelaxZones(i)%gam)
           P(j,RelaxZones(i)%idx(3):RelaxZones(i)%idx(4)) = P(j,RelaxZones(i)%idx(3):RelaxZones(i)%idx(4))*(RelaxZones(i)%gam) &
		       + FAC*RelaxZones(i)%Pa*(one-RelaxZones(i)%gam)
 	     END DO
	   ELSE
	     DO j = RelaxZones(i)%idx(1),RelaxZones(i)%idx(2)
		   k = j - RelaxZones(i)%idx(1) + 1
           CALL stream_func_wave_finite(RelaxZones(i)%idx(4)-RelaxZones(i)%idx(3)+1,&
			    -tmpy(RelaxZones(i)%idx(3):RelaxZones(i)%idx(4))&
				*SIN(RelaxZones(i)%degrees/180.0_long*pi)+tmpx(j)*COS(RelaxZones(i)%degrees/180.0_long*pi),time,&
				SFsol%n_four_modes,SFsol%zz,SFsol%yy,SFsol%k,g,RelaxZones(i)%Ea,RelaxZones(i)%Pa)
           E(j,RelaxZones(i)%idx(3):RelaxZones(i)%idx(4)) = E(j,RelaxZones(i)%idx(3):RelaxZones(i)%idx(4))*(RelaxZones(i)%gam(k)) &
		       + FAC*RelaxZones(i)%Ea*(one-RelaxZones(i)%gam(k))
           P(j,RelaxZones(i)%idx(3):RelaxZones(i)%idx(4)) = P(j,RelaxZones(i)%idx(3):RelaxZones(i)%idx(4))*(RelaxZones(i)%gam(k)) &
		       + FAC*RelaxZones(i)%Pa*(one-RelaxZones(i)%gam(k))
 	     END DO
	   ENDIF
	 ENDIF
  ELSE IF (RelaxZones(i)%XorY=='Y' .AND. RelaxZones(i)%XorYgen=='X') THEN ! Y-direction
	 IF (RelaxZones(i)%WavegenONOFF==0) THEN
	   DO j = RelaxZones(i)%idx(1),RelaxZones(i)%idx(2)
         E(j,RelaxZones(i)%idx(3):RelaxZones(i)%idx(4)) = E(j,RelaxZones(i)%idx(3):RelaxZones(i)%idx(4))*(RelaxZones(i)%gam)
         P(j,RelaxZones(i)%idx(3):RelaxZones(i)%idx(4)) = P(j,RelaxZones(i)%idx(3):RelaxZones(i)%idx(4))*(RelaxZones(i)%gam)
 	   END DO
	 ELSE
	   IF (RelaxZones(i)%degrees==zero) THEN
        CALL stream_func_wave_finite(RelaxZones(i)%idx(2)-RelaxZones(i)%idx(1)+1,&
		   tmpx(RelaxZones(i)%idx(1):RelaxZones(i)%idx(2)),time,&
		   SFsol%n_four_modes,SFsol%zz,SFsol%yy,SFsol%k,g,RelaxZones(i)%Ea,RelaxZones(i)%Pa)
	    DO j = RelaxZones(i)%idx(3),RelaxZones(i)%idx(4)
		  k = j - RelaxZones(i)%idx(3) + 1
          E(RelaxZones(i)%idx(1):RelaxZones(i)%idx(2),j) = E(RelaxZones(i)%idx(1):RelaxZones(i)%idx(2),j)*(RelaxZones(i)%gam(k)) &
		     + FAC*RelaxZones(i)%Ea*(one-RelaxZones(i)%gam(k))
          P(RelaxZones(i)%idx(1):RelaxZones(i)%idx(2),j) = P(RelaxZones(i)%idx(1):RelaxZones(i)%idx(2),j)*(RelaxZones(i)%gam(k)) &
		     + FAC*RelaxZones(i)%Pa*(one-RelaxZones(i)%gam(k))
 	    END DO
	  ELSE
	    DO j = RelaxZones(i)%idx(3),RelaxZones(i)%idx(4)
		  k = j - RelaxZones(i)%idx(3) + 1
          CALL stream_func_wave_finite(RelaxZones(i)%idx(2)-RelaxZones(i)%idx(1)+1,&
		     tmpx(RelaxZones(i)%idx(1):RelaxZones(i)%idx(2))&
			 *COS(RelaxZones(i)%degrees/180.0_long*pi)+tmpy(j)*SIN(RelaxZones(i)%degrees/180.0_long*pi),time,&
			 SFsol%n_four_modes,SFsol%zz,SFsol%yy,SFsol%k,g,RelaxZones(i)%Ea,RelaxZones(i)%Pa)
          E(RelaxZones(i)%idx(1):RelaxZones(i)%idx(2),j) = E(RelaxZones(i)%idx(1):RelaxZones(i)%idx(2),j)*(RelaxZones(i)%gam(k)) &
		     + FAC*RelaxZones(i)%Ea*(one-RelaxZones(i)%gam(k))
          P(RelaxZones(i)%idx(1):RelaxZones(i)%idx(2),j) = P(RelaxZones(i)%idx(1):RelaxZones(i)%idx(2),j)*(RelaxZones(i)%gam(k)) &
		     + FAC*RelaxZones(i)%Pa*(one-RelaxZones(i)%gam(k))
 	    END DO
	  ENDIF
	 ENDIF
   ENDIF
END DO

END SUBROUTINE RelaxationModule


text:
c ********************************************************************
c *                                                                  *
c *  routines to support crack growth by smcs criterion              *
c *                                                                  *
c ********************************************************************
c
c     ****************************************************************
c     *                                                              *
c     *                      subroutine dam_print_elem3              *
c     *                                                              *
c     *                       written by : rhd                       *
c     *                                                              *
c     *                   last modified : 6/18/2019 rhd              *
c     *                                                              *
c     *     This routine prints out the status of SMCS elements      *
c     *     marked as killable at the beginning of a load step       *
c     *                                                              *
c     ****************************************************************
c
c
      subroutine dam_print_elem3( step, iter )
      use global_data, only : out, noelem, nonode, ltmstp, stname
      use elem_extinct_data, only : dam_state, dam_print_list,
     &                              old_plast_strain
      use main_data, only :   output_packets, packet_file_no
      use damage_data, only : dam_ptr, load_size_control_crk_grth,
     &                        num_print_list, num_elements_killed,
     &                        num_kill_elem, print_status,
     &                        print_top_list, num_top_list,
     &                        smcs_states, smcs_stream, smcs_text,
     &                        smcs_states_intlst

c
      implicit none
c
c              parameters
c
      integer :: step, iter
c
c              local declarations
c
      integer :: elem_loop, elem_ptr, element, local_count,
     &           list_count, list_length
      integer, allocatable :: element_list(:), sort_index_vec(:)
      double precision :: dummy, eps_plas, eps_crit, sig_mean,
     &    sig_mises, d_eps_plas, max_d_eps_plas, ddummy1, ddummy2,
     &    ddummy3, triaxiality, mean_zeta, ratio, mean_omega
      double precision, allocatable :: strain_ratios(:)
      double precision, parameter :: zero = 0.0d0, one = 1.0d0,
     &                               eps_plas_tol = 1.0d-06,
     &                               ratio_tol = 0.001d0
      logical :: ldummy, debug, skip_element
      character(len=1) :: cflag
c
      write(out,*) ' '
      write(out,*) ' ********************************************* '
      write(out,*) ' ***         Killable element status       *** '
      write(out,*) ' ********************************************* '
      write(out,*)
      if( load_size_control_crk_grth ) write(out,9020)
      if( .not. load_size_control_crk_grth ) write(out,9030)
c
c              Support for printing elements with largest ratio
c              of plastic strain / critical strain. build sorted
c              list of elements.
c
      if( print_top_list ) then
        allocate( strain_ratios(num_kill_elem),
     &            element_list(num_kill_elem),
     &            sort_index_vec(num_kill_elem) )
        strain_ratios  = zero
        sort_index_vec = 0
        list_count     = 0
        do element = 1, noelem
          elem_ptr = dam_ptr(element)
          if( elem_ptr == 0 ) cycle  ! element not killable
          if( dam_state(elem_ptr) .ne. 0 ) cycle ! already killed
          list_count = list_count + 1
          element_list(list_count) = element
        end do
        call dam_print_elem3_fill
        call  hpsort_eps_epw( list_count, strain_ratios,
     &                        sort_index_vec, 1.0d-07 )
      end if
c
c              print element status for smcs. follow user-specified
c              print order or the sorted list by descending
c              critical strain ratio
c
      max_d_eps_plas = zero
      local_count = 0
      if( print_status )   list_length = num_print_list
      if( print_top_list ) list_length = min( num_top_list,
     &                                        list_count )
c
      do elem_loop = 1, list_length
c
         if( print_top_list ) element =
     &              element_list(sort_index_vec(elem_loop))
         if( print_status ) then
            element = dam_print_list(elem_loop)
            elem_ptr = dam_ptr(element)
            if( elem_ptr == 0 ) cycle   ! not killable
            if( dam_state(elem_ptr) .ne. 0 ) cycle ! already killed
         end if
c
c              get current smcs parameters for element and print.
c              skip printing elements with very small plastic strain
c
         call dam_param_3_get_values( element, debug, eps_plas,
     &                   eps_crit, sig_mean, sig_mises,
     &                   triaxiality, mean_zeta, mean_omega,
     &                   2, ldummy )
         skip_element = eps_plas < eps_plas_tol
         if( skip_element ) cycle
         cflag = " "
         if( triaxiality < one ) cflag = "*"
         if( load_size_control_crk_grth ) then
            d_eps_plas = eps_plas - old_plast_strain(dam_ptr(element))
            max_d_eps_plas = max(max_d_eps_plas, d_eps_plas)
            ratio = eps_plas / eps_crit
            if( ratio >= ratio_tol )
     &       write(out,9110) element, eps_plas, eps_crit,
     &              eps_plas / eps_crit, sig_mean,
     &              sig_mises, d_eps_plas, triaxiality, mean_zeta,
     &              mean_omega, cflag
         else
            ratio = eps_plas / eps_crit
            if( ratio >= ratio_tol )
     &         write(out,9100) element, eps_plas, eps_crit,
     &              ratio, sig_mean,
     &              sig_mises, triaxiality, mean_zeta, mean_omega,
     &              cflag
         end if
         local_count = local_count + 1
      end do  ! elem_loop
c
      write(out,9220) num_elements_killed
      if( local_count > 0 ) then
        write(out,9200); write(out,9210)
      end if
      if( load_size_control_crk_grth ) write(out,9010) max_d_eps_plas
c
c              output an smcs states file for step if requested
c
      if( smcs_states ) call dam_print_elem3_states_file
c
c
c              packet output section. output results for all elements
c              not yet killed.
c              ------------------------------------------------------
c
      if( .not. output_packets ) then
         call dam_print_elem3_rel
         return
      else
         call dam_print_elem3_packets
         write(out,9300) list_count
         call dam_print_elem3_rel
      end if
c
c              flat file of state variables for SMCS
c
c
      return
c
 9010 format(/,5x,'>> Maximum change in plastic strain for this step:',
     &     f10.7,/)
 9020 format(2x,'element   eps-pls    eps-crit  ratio  ',
     & 'sig-mean  sig-mises  d(eps-pls)    triaxiality (T)   zeta',
     & '   omega',/,2x,
     &          '-------   -------    --------  -----  ',
     & '--------  ---------  ----------    ---------------   ----',
     & '   -----')
 9030 format(2x,'element   eps-pls    eps-crit  ratio  ',
     &       'sig-mean  sig-mises  triaxiality (T)   zeta',
     &       '   omega',/,2x,
     &          '-------   -------    --------  -----  ',
     &       '--------  ---------  ---------------   ----   -----')
 9100 format(1x,i8,2(2x,f9.6),2x,f5.3,1x,f9.3,2x,f9.3,4x,f6.2,8x,
     &       f6.2,2x,f6.2,2x,a1)
 9110 format(1x,i8,2(2x,f9.6),2x,f5.3,1x,f9.3,2x,f9.3,3x,f9.7,
     &  7x,f6.2,7x,f6.2,2x,f6.2,2x,a1)
 9200 format(10x,"* -- Triaxiality (sig_mean/sig_mises) < 1.0")
 9210 format(10x,"-- Triaxiality is mean over history weighted by ",
     & "plastic strain",/,
     & 10x,"-- Zeta is mean of normalized Lode angle weighted by ",
     & "plastic strain",/,
     & 10x,"-- Omega is mean of Lode angle parameter weighted by ",
     & "plastic strain",
     & /,10x,"-- Elements with ratio < 0.001 not shown")
 9220 format(10x,"Total elements now killed: ",i5)
 9300 format(10x,"-- Packet data written for element count: ",i7)
c
      contains
c     ========
c     ****************************************************************
c     *                                                              *
c     *               subroutine dam_print_elem3_rel                 *
c     *                                                              *
c     ****************************************************************
c
c
      subroutine dam_print_elem3_rel
      implicit none
c
      if( allocated( strain_ratios ) )  deallocate( strain_ratios )
      if( allocated( element_list ) )   deallocate( element_list )
      if( allocated( sort_index_vec ) ) deallocate( sort_index_vec )
c
      write(out,*) " "
c
      return
c
      end subroutine dam_print_elem3_rel
c
c     ****************************************************************
c     *                                                              *
c     *               subroutine dam_print_elem3_fill                *
c     *                                                              *
c     ****************************************************************
c
      subroutine dam_print_elem3_fill
      implicit none
      integer :: k, element
c
      do k = 1, list_count
         element  = element_list(k)
         call dam_param_3_get_values( element, debug, eps_plas,
     &                   eps_crit, sig_mean, sig_mises,
     &                   triaxiality, mean_zeta, mean_omega,
     &                   2, ldummy )
         strain_ratios(k) =  eps_plas / eps_crit
      end do  ! elem_loop
c
      return
c
      end subroutine dam_print_elem3_fill
c     ****************************************************************
c     *                                                              *
c     *               subroutine dam_print_elem3_packets             *
c     *                                                              *
c     ****************************************************************
c
      subroutine dam_print_elem3_packets
      implicit none
c
      if( .not. allocated( element_list ) )
     &    allocate( element_list(num_kill_elem) )
      list_count     = 0
      do element = 1, noelem
        elem_ptr = dam_ptr(element)
        if( elem_ptr == 0 ) cycle  ! element not killable
        if( dam_state(elem_ptr) .ne. 0 ) cycle ! already killed
        list_count = list_count + 1
        element_list(list_count) = element
      end do
      if( list_count == 0 ) then  ! all elements killed
         call dam_print_elem3_rel
         return
      end if
c
      if( load_size_control_crk_grth )
     &      write( packet_file_no ) 23, list_count, step, iter
      if( .not. load_size_control_crk_grth )
     &        write( packet_file_no ) 22, list_count, step, iter
c
      do elem_loop = 1, list_count
         element  = element_list(elem_loop)
         call dam_param_3_get_values( element, debug, eps_plas,
     &                   eps_crit, sig_mean, sig_mises,
     &                   triaxiality, mean_zeta, mean_omega,
     &                   2, ldummy )
         if( load_size_control_crk_grth ) then
            d_eps_plas = eps_plas - old_plast_strain(dam_ptr(element))
            max_d_eps_plas = max(max_d_eps_plas, d_eps_plas)
            write(packet_file_no) element, eps_plas, eps_crit,
     &              sig_mean, sig_mises, d_eps_plas, triaxiality,
     &              mean_zeta, mean_omega
         else
            write(packet_file_no) element, eps_plas, eps_crit,
     &              sig_mean, sig_mises, triaxiality, mean_zeta,
     &              mean_omega
         end if
      end do   ! elem_loop
c
      return
      end subroutine dam_print_elem3_packets
c     ****************************************************************
c     *                                                              *
c     *          subroutine dam_print_elem3_states_file              *
c     *                                                              *
c     ****************************************************************
c
      subroutine dam_print_elem3_states_file
      implicit none
c
      integer :: header_file_number, flat_file_number, step_num,
     &           now_len, lenlst
      logical ::  get_values, output_this_step
      logical, external :: scan_entry_in_list
      double precision :: evalues(8)
      character(len=24) :: sdate_time_tmp
      character(len=20) :: form_type, access_type
      character(len=80) :: flat_name
c
c              user requested output of element states file with smcs
c              values. write file if this step in user list
c
      step_num = ltmstp  ! current step number
      lenlst = sizeof( smcs_states_intlst )
      output_this_step = scan_entry_in_list( step_num,
     &                      smcs_states_intlst, lenlst )
      if( .not. output_this_step ) return
c
c              open/write/close the states header file. just overwrites
c              an existing file
c
      open( newunit=header_file_number,file="states_header_smcs",
     &      access="sequential", form="formatted" )
      write(header_file_number,9100)
      close(unit=header_file_number,status="keep")
c
c              make name of the smcs states results file for this
c              load(time) step
c                wem#######_stream_smcs   wem#######_text_smcs
c
      flat_name(1:3) ='wem'
      write(flat_name(4:),9200) step_num
c
      if( smcs_text ) then
         flat_name(11:) = '_text'
         now_len = len_trim( flat_name )
         flat_name(now_len+1:) = "_" // "smcs"
         access_type = 'sequential'
         form_type   = 'formatted'
      end if
c
      if( smcs_stream ) then
         flat_name(11:) = '_stream'
         now_len = len_trim( flat_name )
         flat_name(now_len+1:) = "_" // "smcs"
         access_type = 'stream'
         form_type   = 'unformatted'
      end if
c
      open( newunit=flat_file_number,file=flat_name, status='unknown',
     &     access=access_type, form=form_type )
c
      if( smcs_text ) then
          write(flat_file_number,9000)
          write(flat_file_number,9014)
          write(flat_file_number,9020) stname
          write(flat_file_number,9030) nonode, noelem
          call fdate( sdate_time_tmp )
          write(flat_file_number,9040) sdate_time_tmp
          write(flat_file_number,9050) step_num
          write(flat_file_number,9000)
      end if
c
c              write a result record for every element. nonkillable
c              elements have a large critical plastic strain.
c
      do element = 1, noelem
        elem_ptr = dam_ptr(element)
        if( elem_ptr == 0 ) then ! element not killable
          get_values = .true.
        elseif( dam_state(elem_ptr) .ne. 0 ) then !already killed
           evalues = zero
           get_values = .false.
        else ! killable elem not yet killed
           get_values = .true.
        end if
        if( get_values ) then
           call dam_param_3_get_values( element, debug, eps_plas,
     &           eps_crit, sig_mean, sig_mises, triaxiality,
     &           mean_zeta, mean_omega, 3, ldummy )
           evalues(1) = eps_plas
           evalues(2) = eps_crit
           evalues(3) = eps_plas / eps_crit
           evalues(4) = sig_mean
           evalues(5) = sig_mises
           evalues(6) = triaxiality
           evalues(7) = mean_zeta
           evalues(8) = mean_omega
        end if
        where( abs( evalues ) .lt. eps_plas_tol ) evalues = zero
       if( smcs_stream ) write(flat_file_number) evalues
       if( smcs_text )   write(flat_file_number,9110) evalues
c
      end do ! over all elements
c
      close(unit=flat_file_number,status="keep")

      return
c
 9000 format('#')
 9014 format('#  WARP3D element results: states-smcs')
 9020 format('#  Structure name: ',a8 )
 9030 format('#  Model nodes, elements: ',2i8)
 9040 format('#  ',a24)
 9050 format('#  Load(time) step: ',i8 )
 9060 format(2i8)
 9100 format(
     &"#",/,
     &"#  header file for smcs state output",/,
     &"#",/,
     &"#  8 character state labels and longer descriptors",/,
     &"#  material model number, number of state variables",/,
     &"#",/,
     &"     8",/,
     &"      1  epspls   Eq. plastic strain",/,
     &"      2  epscrit  Critical strain",/,
     &"      3  ratio    epspls/epscrit",/,
     &"      4  mean     mean stress",/,
     &"      5  mises    mean mises",/,
     &"      6  triax    mean triaxiality",/,
     &"      7  xi       mean xi",/,
     &"      8  omega    mean omega" )
 9110 format(8e15.6)
 9200 format( i7.7 )
c
      end subroutine dam_print_elem3_states_file
c
      end subroutine dam_print_elem3


text:
Address(Kent Drive) is a quiet residential street in Central Davis.   It is part of the college streets in the towns street name themes, named for Kent State University (it was named prior to the Kent State massacre). It runs between Apple Lane and Sycamore Lane. Proceeding north, Kent is the second to last street to intersect Sycamore before West Covell Boulevard Covell.

First Lady Rosalynn Carter once had dinner at then mayor Tom Tomasi’s house on Kent Drive in the late 70’s.
Kent Drive is the only street in Davis to have a private tennis court.

Current locally famous residents of Kent include:
  Ken Wagstaff
  Don Saylor

Former:
  Thomas J Tomasi

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
      REAL FUNCTION O2ABS(TEMP,PRES,VAPDEN,FREQ)
C  Copyright (c) 2003 Massachusetts Institute of Technology
C
C     PURPOSE: RETURNS ABSORPTION COEFFICIENT DUE TO OXYGEN IN AIR,
C              IN NEPERS/KM
C
C      5/1/95  P. Rosenkranz
C      11/5/97  P. Rosenkranz - 1- line modification.
c      12/16/98 pwr - updated submm freq's and intensities from HITRAN96
c      8/21/02  pwr - revised width at 425
c      3/20/03  pwr - 1- line mixing and width revised
C
c     IMPLICIT NONE
C
C     ARGUMENTS:
      REAL TEMP,PRES,VAPDEN,FREQ
C
C     NAME    UNITS    DESCRIPTION        VALID RANGE
C
C     TEMP    KELVIN   TEMPERATURE        UNCERTAIN, but believed to be
c                                          valid for atmosphere
C     PRES   MILLIBARS PRESSURE           3 TO 1000
C     VAPDEN  G/M**3   WATER VAPOR DENSITY  (ENTERS LINEWIDTH CALCULATION
C                      DUE TO GREATER BROADENING EFFICIENCY OF H2O)
C     FREQ    GHZ      FREQUENCY          0 TO 900
C
C     REFERENCES FOR EQUATIONS AND COEFFICIENTS:
C     P.W. Rosenkranz, CHAP. 2 and appendix, in ATMOSPHERIC REMOTE SENSING
C      BY MICROWAVE RADIOMETRY (M.A. Janssen, ed., 1993).
C     H.J. Liebe et al, JQSRT V.48, pp.629-643 (1992).
c     M.J. Schwartz, Ph.D. thesis, M.I.T. (1998).
c     A.F. Krupnov et al, J. Mol. Spect. v.215, pp.309-311 (2002).
C     M.Yu. Tretyakov et al, J. Mol. Spect. (2003 preprint).
C     SUBMILLIMETER LINE INTENSITIES FROM HITRAN96.
c
c     This version differs from Liebe's MPM92 in these significant respects:
c     1. The 1- line has the width and mixing coefficient measured by
c      Tretyakov et al.
c     2. It modifies the 1- line width temperature dependence to (1/T)**0.9
c     3. It uses the same temperature dependence (X) for submillimeter
c      line widths as in the 60 GHz band: (1/T)**0.8
c     4. The 425 GHz line width is from Krupnov et al.
C
c     Local variables:
      REAL TH,TH1,B,PRESWV,PRESDA,DEN,DENS,DFNR,SUM,STR,Y,SF1,SF2,FCEN
      INTEGER K
      REAL X,WB300,W300(40),F(40),Y300(40),S300(40),V(40),BE(40)
      COMMON /O2COM/ X,WB300,W300,F,Y300,S300,V,BE
C      LINES ARE ARRANGED 1-,1+,3-,3+,ETC. IN SPIN-ROTATION SPECTRUM
      DATA F/118.7503, 56.2648, 62.4863, 58.4466, 60.3061, 59.5910,
     2  59.1642, 60.4348, 58.3239, 61.1506, 57.6125, 61.8002,
     3  56.9682, 62.4112, 56.3634, 62.9980, 55.7838, 63.5685,
     4  55.2214, 64.1278, 54.6712, 64.6789, 54.1300, 65.2241,
     5  53.5957, 65.7648, 53.0669, 66.3021, 52.5424, 66.8368,
     6  52.0214, 67.3696, 51.5034, 67.9009, 368.4984, 424.7632,
     7  487.2494, 715.3931, 773.8397, 834.1458/
        DATA S300/.2936E-14,.8079E-15, .2480E-14,.2228E-14,
     &  .3351E-14,.3292E-14, .3721E-14,.3891E-14,
     &  .3640E-14,.4005E-14, .3227E-14,.3715E-14,
     &  .2627E-14,.3156E-14, .1982E-14,.2477E-14,
     &  .1391E-14,.1808E-14, .9124E-15,.1230E-14,
     &  .5603E-15,.7842E-15, .3228E-15,.4689E-15,
     &  .1748E-15,.2632E-15, .8898E-16,.1389E-15,
     &  .4264E-16,.6899E-16, .1924E-16,.3229E-16,
     &  .8191E-17,.1423E-16, .6494E-15, .7083E-14, .3025E-14,
     &  .1835E-14, .1158E-13, .3993E-14/
      DATA BE/.009,.015, .083,.084, 2*.212, 2*.391, 2*.626,
     & 2*.915, 2*1.260, 1.660,1.665, 2.119,2.115, 2.624,2.625,
     & 2*3.194, 2*3.814, 2*4.484, 2*5.224, 2*6.004, 2*6.844,
     & 2*7.744, .048, .044, .049, .145, .141, .145/
C      WIDTHS IN MHZ/MB
      DATA WB300/.56/, X/.8/
      DATA W300/1.67, 1.646, 1.468, 1.449, 1.382, 1.360,
     & 1.319, 1.297, 1.266, 1.248, 1.221, 1.207, 1.181, 1.171,
     & 1.144, 1.139, 1.110, 1.108, 1.079, 1.078, 2*1.05,
     & 2*1.02,2*1.00,2*.97,2*.94,2*.92,2*.89, 3*1.64, 3*1.81/
      DATA Y300/  -0.036,  0.2408, -0.3486,  0.5227,
     & -0.5430,  0.5877, -0.3970,  0.3237, -0.1348,  0.0311,
     &  0.0725, -0.1663,  0.2832, -0.3629,  0.3970, -0.4599,
     &  0.4695, -0.5199,  0.5187, -0.5597,  0.5903, -0.6246,
     &  0.6656, -0.6942,  0.7086, -0.7325,  0.7348, -0.7546,
     &  0.7702, -0.7864,  0.8083, -0.8210,  0.8439, -0.8529, 6*0./
      DATA V/  0.0079, -0.0978,  0.0844, -0.1273,
     &  0.0699, -0.0776,  0.2309, -0.2825,  0.0436, -0.0584,
     &  0.6056, -0.6619,  0.6451, -0.6759,  0.6547, -0.6675,
     &  0.6135, -0.6139,  0.2952, -0.2895,  0.2654, -0.2590,
     &  0.3750, -0.3680,  0.5085, -0.5002,  0.6206, -0.6091,
     &  0.6526, -0.6393,  0.6640, -0.6475,  0.6729, -0.6545, 6*0./
C
      TH = 300./TEMP
      TH1 = TH-1.
      B = TH**X
      PRESWV = VAPDEN*TEMP/217.
      PRESDA = PRES -PRESWV
      DEN = .001*(PRESDA*B + 1.1*PRESWV*TH)
      DENS = .001*(PRESDA*TH**.9 + 1.1*PRESWV*TH)
      DFNR = WB300*DEN
      SUM = 1.6E-17*FREQ*FREQ*DFNR/(TH*(FREQ*FREQ + DFNR*DFNR))
      DO 32 K=1,40
      IF(K.EQ.1) THEN !exception for 1- line
        DF = W300(1)*DENS
      ELSE
        DF = W300(K)*DEN
      ENDIF
      FCEN = F(K)
      Y = .001*PRES*B*(Y300(K)+V(K)*TH1)
      STR = S300(K)*EXP(-BE(K)*TH1)
      SF1 = (DF + (FREQ-FCEN)*Y)/((FREQ-FCEN)**2 + DF*DF)
      SF2 = (DF - (FREQ+FCEN)*Y)/((FREQ+FCEN)**2 + DF*DF)
32    SUM = SUM + STR*(SF1+SF2)*(FREQ/F(K))**2
      O2ABS = .5034E12*SUM*PRESDA*TH**3/3.14159
      O2ABS = AMAX1(O2ABS,0.)
      RETURN
      END

max_stars_repo_name: hsbadr/crtm

max_issues_repo_name: hsbadr/crtm

max_forks_repo_name: hsbadr/crtm


text:
subroutine diffuse_scalar3D (field,fluxb,fluxt,tkh,rho,rhow,flux)

use grid
use params, only: docolumn,dowallx,dowally,dosgs
use sgs, only: grdf_x,grdf_y,grdf_z
implicit none
! input
real field(dimx1_s:dimx2_s, dimy1_s:dimy2_s, nzm)	! scalar
real tkh(0:nxp1,1-YES3D:nyp1,nzm)	! eddy conductivity
real fluxb(nx,ny)		! bottom flux
real fluxt(nx,ny)		! top flux
real rho(nzm)
real rhow(nz)
real flux(nz)
! local
real flx(0:nx,0:ny,0:nzm)
real dfdt(nx,ny,nz)
real rdx2,rdy2,rdz2,rdz,rdx5,rdy5,rdz5,tmp
real dxy,dxz,dyx,dyz,dzx,dzy,tkx,tky,tkz,rhoi
integer i,j,k,ib,ic,jb,jc,kc,kb

if(.not.dosgs) return

rdx2=1./(dx*dx)
rdy2=1./(dy*dy)
rdz2=1./(dz*dz)
rdz=1./dz
dxy=dx/dy
dxz=dx/dz
dyx=dy/dx
dyz=dy/dz
dzx=dz/dx
dzy=dz/dy

dfdt(:,:,:)=0.

!-----------------------------------------
if(dowallx) then

  if(mod(rank,nsubdomains_x).eq.0) then
    do k=1,nzm
     do j=1,ny
         field(0,j,k) = field(1,j,k)
     end do
    end do
  end if
  if(mod(rank,nsubdomains_x).eq.nsubdomains_x-1) then
    do k=1,nzm
     do j=1,ny
         field(nx+1,j,k) = field(nx,j,k)
     end do
    end do
  end if

end if

if(dowally) then

  if(rank.lt.nsubdomains_x) then
    do k=1,nzm
       do i=1,nx
         field(i,1-YES3D,k) = field(i,1,k)
       end do
    end do
  end if
  if(rank.gt.nsubdomains-nsubdomains_x-1) then
    do k=1,nzm
       do i=1,nx
         field(i,ny+YES3D,k) = field(i,ny,k)
       end do
    end do
  end if

end if

!-----------------------------------------

!  Horizontal diffusion:

do k=1,nzm

 rdx5=0.5*rdx2  * grdf_x(k)
 rdy5=0.5*rdy2  * grdf_y(k)

 do j=1,ny
  do i=0,nx
    ic=i+1
    tkx=rdx5*(tkh(i,j,k)+tkh(ic,j,k))
    flx(i,j,k)=-tkx*(field(ic,j,k)-field(i,j,k))
  end do
  do i=1,nx
    ib=i-1
    dfdt(i,j,k)=dfdt(i,j,k)-(flx(i,j,k)-flx(ib,j,k))
  end do
 end do

 do j=0,ny
  jc=j+1
  do i=1,nx
   tky=rdy5*(tkh(i,j,k)+tkh(i,jc,k))
   flx(i,j,k)=-tky*(field(i,jc,k)-field(i,j,k))
  end do
 end do
 do j=1,ny
  jb=j-1
  do i=1,nx
    dfdt(i,j,k)=dfdt(i,j,k)-(flx(i,j,k)-flx(i,jb,k))
  end do
 end do

end do ! k

!  Vertical diffusion:

flux(1) = 0.
tmp=1./adzw(nz)
do j=1,ny
 do i=1,nx
   flx(i,j,0)=fluxb(i,j)*rdz*rhow(1)
   flx(i,j,nzm)=fluxt(i,j)*rdz*tmp*rhow(nz)
   flux(1) = flux(1) + flx(i,j,0)
 end do
end do

do k=1,nzm-1
 kc=k+1
 flux(kc)=0.
 rhoi = rhow(kc)/adzw(kc)
 rdz5=0.5*rdz2 * grdf_z(k)
 do j=1,ny
  do i=1,nx
    tkz=rdz5*(tkh(i,j,k)+tkh(i,j,kc))
    flx(i,j,k)=-tkz*(field(i,j,kc)-field(i,j,k))*rhoi
    flux(kc) = flux(kc) + flx(i,j,k)
  end do
 end do
end do

do k=1,nzm
 kb=k-1
 rhoi = 1./(adz(k)*rho(k))
 do j=1,ny
  do i=1,nx
   dfdt(i,j,k)=dtn*(dfdt(i,j,k)-(flx(i,j,k)-flx(i,j,kb))*rhoi)
   field(i,j,k)=field(i,j,k)+dfdt(i,j,k)
  end do
 end do
end do

end subroutine diffuse_scalar3D

max_stars_repo_name: kbren/uwnet

max_issues_repo_name: kbren/uwnet

max_forks_repo_name: kbren/uwnet


text:
       SUBROUTINE W3FT12(COEF,WORK,GRID,TRIGS)                          00010000
C$$$   SUBPROGRAM  DOCUMENTATION  BLOCK                                 00020000
C                                                                       00030000
C SUBPROGRAM: W3FT12         FAST FOURIER FOR 2.5 DEGREE GRID           00040000
C   AUTHOR: SELA,JOE         ORG: W323       DATE: 80-11-21             00050000
C                                                                       00060000
C ABSTRACT: FAST FOURIER TO COMPUTE 145 GRID VALUES AT DESIRED          00070000
C   LATITUDE FROM 31 COMPLEX FOURIER COEFFICIENTS. THIS SUBROUTINE      00080000
C   IS SPECIAL PURPOSE FOR CONVERTING COEFFICIENTS TO A 2.5 DEGREE      00090000
C   LAT,LON GRID.                                                       00100000
C                                                                       00110000
C PROGRAM HISTORY LOG:                                                  00120000
C   80-11-21  JOE SELA                                                  00130000
C   84-06-21  R.E.JONES   CHANGE TO IBM VS FORTRAN                      00140000
C   93-04-12  R.E.JONES   CHANGE TO CRAY CFT77 FORTRAN                  00150000
C                                                                       00160000
C USAGE:  CALL W3FT12(COEF,WORK,GRID,TRIGS)                             00170000
C                                                                       00180000
C   INPUT VARIABLES:                                                    00190000
C     NAMES  INTERFACE DESCRIPTION OF VARIABLES AND TYPES               00200000
C     ------ --------- -----------------------------------------------  00210000
C     COEF   ARG LIST  31 COMPLEX FOURIER COEFFICIENTS.                 00220000
C     TRIGS  ARG LIST  216 TRIG FUNCTIONS ASSUMED PRECOMPUTED           00230000
C                      BY W3FA13 BEFORE FIRST CALL TO W3FT12.           00240000
C     WORK   ARG LIST  144 REAL WORK SPACE                              00250000
C                                                                       00260000
C   OUTPUT VARIABLES:                                                   00270000
C     NAMES  INTERFACE DESCRIPTION OF VARIABLES AND TYPES               00280000
C     ------ --------- -----------------------------------------------  00290000
C     GRID   ARG LIST  145 GRID VALUES, GRID(1)=GRID(145)               00300000
C                                                                       00310000
C ATTRIBUTES:                                                           00320000
C   LANGUAGE: CRAY CFT77 FORTRAN 77                                     00330000
C   MACHINE:  CRAY C916-128, Y-MP8/864, Y-MP EL92/256                   00340002
C                                                                       00350000
C$$$                                                                    00360000
C                                                                       00370000
       REAL            COEF( 62 )                                       00380000
       REAL            GRID(145)                                        00390000
       REAL            TRIGS(216)                                       00400000
       REAL            WORK(144)                                        00410000
C                                                                       00411003
       SAVE                                                             00412003
C                                                                       00420000
       DATA  SIN60/0.866025403784437/                                   00430000
C                                                                       00440000
      GRID(1) = COEF(1)                                                 00450000
      GRID(2) = COEF(1)                                                 00460000
      K = 147                                                           00470000
      J = 143                                                           00480000
      DO 100 I=3, 61 ,2                                                 00490000
      TEMP = COEF(I)*TRIGS(K+1) - COEF(I+1)*TRIGS(K)                    00500000
      GRID(I) = COEF(I) - TEMP                                          00510000
      GRID(J) = COEF(I) + TEMP                                          00520000
      TEMP = COEF(I)*TRIGS(K) + COEF(I+1)*TRIGS(K+1)                    00530000
      GRID(I+1) = TEMP - COEF(I+1)                                      00540000
      GRID(J+1) = TEMP + COEF(I+1)                                      00550000
      K = K + 2                                                         00560000
      J = J - 2                                                         00570000
100   CONTINUE                                                          00580000
      DO 110 I= 63 , 84                                                 00590000
      GRID(I) = 0.0                                                     00600000
110   CONTINUE                                                          00610000
C                                                                       00620000
      A0 = GRID(1) + GRID(73)                                           00630000
      A2 = GRID(1) - GRID(73)                                           00640000
      B0 = GRID(2) + GRID(74)                                           00650000
      B2 = GRID(2) - GRID(74)                                           00660000
      A1 = GRID(37) + GRID(109)                                         00670000
      A3 = GRID(37) - GRID(109)                                         00680000
      B1 = GRID(38) + GRID(110)                                         00690000
      B3 = GRID(38) - GRID(110)                                         00700000
      WORK(1) = A0 + A1                                                 00710000
      WORK(5) = A0 - A1                                                 00720000
      WORK(2) = B0 + B1                                                 00730000
      WORK(6) = B0 - B1                                                 00740000
      WORK(3) = A2 - B3                                                 00750000
      WORK(7) = A2 + B3                                                 00760000
      WORK(4) = B2 + A3                                                 00770000
      WORK(8) = B2 - A3                                                 00780000
      KB = 3                                                            00790000
      KC = 5                                                            00800000
      KD = 7                                                            00810000
      J = 75                                                            00820000
      K = 39                                                            00830000
      L = 111                                                           00840000
      M = 9                                                             00850000
      DO 300 I=3,35,2                                                   00860000
      A0 = GRID(I) + GRID(J)                                            00870000
      A2 = GRID(I) - GRID(J)                                            00880000
      B0 = GRID(I+1) + GRID(J+1)                                        00890000
      B2 = GRID(I+1) - GRID(J+1)                                        00900000
      A1 = GRID(K) + GRID(L)                                            00910000
      A3 = GRID(K) - GRID(L)                                            00920000
      B1 = GRID(K+1) + GRID(L+1)                                        00930000
      B3 = GRID(K+1) - GRID(L+1)                                        00940000
      WORK(M  ) = A0 + A1                                               00950000
      WORK(M+4) = A0 - A1                                               00960000
      WORK(M+1) = B0 + B1                                               00970000
      WORK(M+5) = B0 - B1                                               00980000
      WORK(M+2) = A2 - B3                                               00990000
      WORK(M+6) = A2 + B3                                               01000000
      WORK(M+3) = B2 + A3                                               01010000
      WORK(M+7) = B2 - A3                                               01020000
      TEMP = WORK(M+2)*TRIGS(KB) - WORK(M+3)*TRIGS(KB+1)                01030000
      WORK(M+3) = WORK(M+2)*TRIGS(KB+1) + WORK(M+3)*TRIGS(KB)           01040000
      WORK(M+2) = TEMP                                                  01050000
      TEMP = WORK(M+4)*TRIGS(KC) - WORK(M+5)*TRIGS(KC+1)                01060000
      WORK(M+5) = WORK(M+4)*TRIGS(KC+1) + WORK(M+5)*TRIGS(KC)           01070000
      WORK(M+4) = TEMP                                                  01080000
      TEMP = WORK(M+6)*TRIGS(KD) - WORK(M+7)*TRIGS(KD+1)                01090000
      WORK(M+7) = WORK(M+6)*TRIGS(KD+1) + WORK(M+7)*TRIGS(KD)           01100000
      WORK(M+6) = TEMP                                                  01110000
      J = J + 2                                                         01120000
      K = K + 2                                                         01130000
      L = L + 2                                                         01140000
      KB = KB + 2                                                       01150000
      KC = KC + 4                                                       01160000
      KD = KD + 6                                                       01170000
      M = M + 8                                                         01180000
300   CONTINUE                                                          01190000
C                                                                       01200000
      I = 1                                                             01210000
      J = 1                                                             01220000
      K = 73                                                            01230000
      DO 440 L=1,4                                                      01240000
      GRID(I) = WORK(J) + WORK(K)                                       01250000
      GRID(I+8) = WORK(J) - WORK(K)                                     01260000
      GRID(I+1) = WORK(J+1) + WORK(K+1)                                 01270000
      GRID(I+9) = WORK(J+1) - WORK(K+1)                                 01280000
      I = I + 2                                                         01290000
      J = J + 2                                                         01300000
      K = K + 2                                                         01310000
440   CONTINUE                                                          01320000
      DO 500 KB=9,65,8                                                  01330000
      I = I + 8                                                         01340000
      DO 460 L=1,4                                                      01350000
      GRID(I) = WORK(J) + WORK(K)                                       01360000
      GRID(I+8) = WORK(J) - WORK(K)                                     01370000
      GRID(I+1) = WORK(J+1) + WORK(K+1)                                 01380000
      GRID(I+9) = WORK(J+1) - WORK(K+1)                                 01390000
      TEMP = GRID(I+8)*TRIGS(KB) - GRID(I+9)*TRIGS(KB+1)                01400000
      GRID(I+9) = GRID(I+8)*TRIGS(KB+1) + GRID(I+9)*TRIGS(KB)           01410000
      GRID(I+8) = TEMP                                                  01420000
      I = I + 2                                                         01430000
      J = J + 2                                                         01440000
      K = K + 2                                                         01450000
460   CONTINUE                                                          01460000
500   CONTINUE                                                          01470000
C                                                                       01480000
      I = 1                                                             01490000
      L = 1                                                             01500000
      KC = 1                                                            01510000
      J = 49                                                            01520000
      K = 97                                                            01530000
      M = 17                                                            01540000
      N = 33                                                            01550000
      DO 660 LL=1,8                                                     01560000
      A1 = GRID(J) + GRID(K)                                            01570000
      A3 = SIN60*(GRID(J)-GRID(K))                                      01580000
      B1 = GRID(J+1) + GRID(K+1)                                        01590000
      B3 = SIN60*(GRID(J+1)-GRID(K+1))                                  01600000
      WORK(L) = GRID(I) + A1                                            01610000
      A2 = GRID(I) - 0.5*A1                                             01620000
      WORK(L+1) = GRID(I+1) + B1                                        01630000
      B2 = GRID(I+1) - 0.5*B1                                           01640000
      WORK(N) = A2 + B3                                                 01650000
      WORK(M) = A2 - B3                                                 01660000
      WORK(M+1) = B2 + A3                                               01670000
      WORK(N+1) = B2 - A3                                               01680000
      I = I + 2                                                         01690000
      J = J + 2                                                         01700000
      K = K + 2                                                         01710000
      L = L + 2                                                         01720000
      M = M + 2                                                         01730000
      N = N + 2                                                         01740000
660   CONTINUE                                                          01750000
      DO 700 KB=17,33,16                                                01760000
      L = L + 32                                                        01770000
      M = M + 32                                                        01780000
      N = N + 32                                                        01790000
      KC = KC + 32                                                      01800000
      DO 680 LL=1,8                                                     01810000
      A1 = GRID(J) + GRID(K)                                            01820000
      A3 = SIN60*(GRID(J)-GRID(K))                                      01830000
      B1 = GRID(J+1) + GRID(K+1)                                        01840000
      B3 = SIN60*(GRID(J+1)-GRID(K+1))                                  01850000
      WORK(L) = GRID(I) + A1                                            01860000
      A2 = GRID(I) - 0.5*A1                                             01870000
      WORK(L+1) = GRID(I+1) + B1                                        01880000
      B2 = GRID(I+1) - 0.5*B1                                           01890000
      WORK(N) = A2 + B3                                                 01900000
      WORK(M) = A2 - B3                                                 01910000
      WORK(M+1) = B2 + A3                                               01920000
      WORK(N+1) = B2 - A3                                               01930000
      TEMP = WORK(M)*TRIGS(KB) - WORK(M+1)*TRIGS(KB+1)                  01940000
      WORK(M+1) = WORK(M)*TRIGS(KB+1) + WORK(M+1)*TRIGS(KB)             01950000
      WORK(M) = TEMP                                                    01960000
      TEMP = WORK(N)*TRIGS(KC) - WORK(N+1)*TRIGS(KC+1)                  01970000
      WORK(N+1) = WORK(N)*TRIGS(KC+1) + WORK(N+1)*TRIGS(KC)             01980000
      WORK(N) = TEMP                                                    01990000
      I = I + 2                                                         02000000
      J = J + 2                                                         02010000
      K = K + 2                                                         02020000
      L = L + 2                                                         02030000
      M = M + 2                                                         02040000
      N = N + 2                                                         02050000
680   CONTINUE                                                          02060000
700   CONTINUE                                                          02070000
C                                                                       02080000
      J = 49                                                            02090000
      K = 97                                                            02100000
      L = 144                                                           02110000
      M = 96                                                            02120000
      N = 48                                                            02130000
      DO 900 I=1,47,2                                                   02140000
      A1 = WORK(J) + WORK(K)                                            02150000
      A3 = SIN60 * (WORK(J)-WORK(K))                                    02160000
      B3 = SIN60 * (WORK(J+1)-WORK(K+1))                                02170000
      B1 = WORK(J+1) + WORK(K+1)                                        02180000
      GRID(L+1) = WORK(I) + A1                                          02190000
      A2 = WORK(I) - 0.5*A1                                             02200000
      B2 = WORK(I+1) - 0.5*B1                                           02210000
      GRID(L) = WORK(I+1) + B1                                          02220000
      GRID(N+1) = A2 + B3                                               02230000
      GRID(M+1) = A2 - B3                                               02240000
      GRID(M) = B2 + A3                                                 02250000
      GRID(N) = B2 - A3                                                 02260000
      J = J + 2                                                         02270000
      K = K + 2                                                         02280000
      L = L - 2                                                         02290000
      M = M - 2                                                         02300000
      N = N - 2                                                         02310000
900   CONTINUE                                                          02320000
      GRID(1) = GRID(145)                                               02330000
C                                                                       02340000
      RETURN                                                            02350000
      END                                                               02360000

max_stars_repo_name: GMAO-SI-Team/jedi-stack

max_issues_repo_name: GMAO-SI-Team/jedi-stack

max_forks_repo_name: GMAO-SI-Team/jedi-stack


text:
module global
	character (len=9) c_starttime
	integer cycle_num
        character (len=32)in_file
end module

max_stars_repo_name: timkphd/examples

max_issues_repo_name: timkphd/examples

max_forks_repo_name: timkphd/examples


text:
!
! ::::::::::::::::::::::: INTFIL ::::::::::::::::::::::::::::::::;
!  INTFIL: interpolates values for a patch at the specified level and
!  location, using values from grids at LEVEL and coarser, if nec.
!
!  take the intersection of a grid patch with corners at ilo,ihi,jlo,jhi
!  and all grids mptr at LEVEL.  If there is a non-null intersection
!  copy the solution vaues from mptr (at TIME) into VAL array.
!  assumes patch at same levels do straight copy, not skipping
!  every intrat or doing any interpolation here,
!  assume called in correct order of levels, so that when copying
!  is ok to overwrite.
!
!  N.B.: there are no dummy points around patch, since
!        this is not an official "grid" but a "patch".
!
!  used array marks when point filled. filpatch checks if points left over
!  after intersections at specified level.
! :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::;
!
!
!> Fill values for a patch at the specified level and
!! location, using values from grids at level **level** ONLY!
!! Leave cells that are not filled unchanged.
!!
!! Search the intersection of a grid patch with corners at **ilo**,**ihi**,
!! **jlo**,**jhi** and all grids **mptr** on level **level**.
!! If there is a non-null intersection,
!! copy the solution \f$q\f$ and auxiliary values from grid **mptr** into **val** array.
!!
!! Assume called in correct order of levels, so that when copying
!! is ok to overwrite.
!!
!! It uses array, **flaguse**, to indicate whether each cell is filled.
!! All cells outside computational domain will be marked as "used" as well
!! and will be processed later by [bc2amr()](@ref bc2amr) in [filrecur()](@ref filrecur).
!!
!! **Input**:
!! * a patch that needs to be filled
!! * global indices that describe the patch to be filled
!! * the grid **msrc** that might contain with the patch
!! * relative position of the patch to grid **msrc** (nrowst, ncolst)
!!
!! **Output**:
!! * data array **val** that stores the new values
!! * flagging array, **flaguse** that indicates whether a cell is filled
!!
!! ## Algorithm
!!
!! If **msrc** \f$ = -1 \f$, this patch is not associated with any grid.
!! So the ''group of grids'', **S**, below contains all level **level** grids.
!!
!! If **msrc** \f$ \neq -1 \f$, this patch is contained in grid **msrc**.
!! So the ''group of grids'', **S**, below only contains
!! level **level** grids that intersect with grid **msrc**.
!! Size of **S** is smaller in this case.
!!
!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!!  procedure intfil()
!!      for each grid mptr in a group of grids, S, do
!!          if grid mptr intersect with the patch
!!              if grid mptr has values at the time needed
!!                  copy values inside the intersection from grid mptr to val
!!                  mark cells in this intersection as "used"
!!              else
!!                  interpolate from values on grid mptr at two different time
!!                  fill the intersection with interpolated values
!!                  mark cells in this intersection as "used"
!!              end if
!!          end if
!!      end for
!!      mark any portion of the patch that is out of whole computational domain as "used",
!!      which will be processed by bc2amr elsewhere later
!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!!
!!
!! \param val array where values are copied to. The array covers the whole grid **msrc**
!! \param mi size of **val** array in *i* direction
!! \param mj size of **val** array in *j* direction
!! \param time simulation time of the values in **val** array
!! \param flaguse marks indicating whether a position in **val** is filled
!! \param nrowst local *i* index (relative to lower-left corner of grid **msrc**) of the left-most cell of the patch to be filled
!! \param ncolst local *j* index (relative to lower-left corner of grid **msrc**) of the lower-most cell of the patch to be filled
!! \param ilo global *i* index of left-most cell  of the patch to be filled
!! \param ihi global *i* index of right-most cell of the patch to be filled
!! \param jlo global *j* index of lower-most cell of the patch to be filled
!! \param jhi global *j* index of upper-most cell of the patch to be filled
!! \param level This patch is on level **level**
!! \param nvar number of equations for the system
!! \param naux number of auxiliary variables
!! \param msrc index of the grid that contains this patch
!!
!! \callgraph
!! \callergraph

subroutine intfil(val,mi,mj,time,flaguse,nrowst,ncolst,ilo,ihi,jlo,jhi,level,nvar,naux,msrc)

    use amr_module, only: possk, mxnest, iregsz, jregsz, nghost, outunit, alloc
    use amr_module, only: node, lstart, store1, store2, levelptr, timemult,gridNbor
    use amr_module, only: rnode, ndilo, ndihi, ndjlo, ndjhi, nextfree
    use amr_module, only: bndListNum, bndListSt
    use amr_module, only: listStart, listOfGrids, bndList, numgrids

    implicit none

    ! Input
    integer, intent(in) :: mi, mj, nrowst, ncolst, ilo, ihi, jlo, jhi, level, nvar, naux,msrc
    real(kind=8), intent(in) :: time

    ! In/Out
    integer(kind=1), intent(in out) :: flaguse(ilo:ihi, jlo:jhi)
    real(kind=8), intent(in out) :: val(nvar,mi,mj)

    ! Locals
    integer :: imlo, jmlo, imhi, jmhi, nx, ny, mitot, mjtot
    integer :: ixlo, ixhi, jxlo, jxhi, locold, locnew, nextSpot
    integer :: icount, bndNum, bndLoc, levSt
    integer :: ivar, i, j, mptr, mstart, loc, numg
    real(kind=8) :: dt, alphac, alphai
    logical :: t_interpolate

    integer :: patch_rect(4)

    real(kind=8), parameter :: t_epsilon = 1.0d-4

    ! Formats for error statements
    character(len=*), parameter :: missing_error = &
            "(' time wanted ',e15.7,' not available from grid ',i4,'level',i4)"
    character(len=*), parameter :: time_error = &
            "(' trying to interpolate from previous time values ',/," // &
            "' for a patch with corners ilo,ihi,jlo,jhi:'" // &
            ",/,2x,4i10,/," // &
            "' from source grid ',i4,' at level ',i4,/," // &
            "' time wanted ',e24.16,' source time is ',e24.16,/," // &
            "' alphai, t_epsilon ',2e24.16)"

    patch_rect = [ilo,ihi,jlo,jhi]

    ! Note that we need a non-dimensionalized t epspatch_rect(1)n as it was a problem
    ! in tsunami tests ( instead of t_epsilon   = dt / 10.d0 )

    ! Time step at this level
    dt = possk(level)

    ! Initialize the flagging where we set things
    flaguse = 0

    ! Loop through all grids at this level, initialize to first
!    mptr = lstart(level)
!    do while (mptr /= 0)
     if (msrc .eq. -1) then
         numg = numgrids(level)
         levSt = listStart(level)
     else
         bndLoc = node(bndListSt,msrc)  ! index of first grid in bndList
         bndNum = node(bndListNum,msrc)
         nextSpot = node(bndListSt, msrc) ! initialize
         numg = bndNum
     endif

     do icount = 1, numg

         if (msrc .eq. -1) then
            mptr = listOfGrids(levSt+icount-1)
         else
            mptr = bndList(nextSpot,gridNbor)
         endif

        ! Check if grid mptr and patch intersect
        imlo = node(ndilo, mptr)
        jmlo = node(ndjlo, mptr)
        imhi = node(ndihi, mptr)
        jmhi = node(ndjhi, mptr)

        nx = node(ndihi,mptr) - node(ndilo,mptr) + 1
        ny = node(ndjhi,mptr) - node(ndjlo,mptr) + 1

        mitot = nx + 2 * nghost
        mjtot = ny + 2 * nghost

        ixlo = max(imlo,patch_rect(1))
        ixhi = min(imhi,patch_rect(2))
        jxlo = max(jmlo,patch_rect(3))
        jxhi = min(jmhi,patch_rect(4))

        ! Check to see if grid and patch interesect, if not continue to next
        ! grid in the list
        if (ixlo <= ixhi .and. jxlo <= jxhi) then

            ! grids intersect. figure out what time to use.
            ! alphai = 1 for new time; 0 for old time
            alphac = (rnode(timemult,mptr) - time)/dt
            alphai = 1.d0 - alphac

            if ((alphai < -t_epsilon) .or. (alphai > 1.d0 + t_epsilon)) then
                write(outunit,missing_error) time, mptr, level
                print missing_error, time, mptr, level
                write(outunit,'(A,E24.16)') 'Line 80', dt
                write(outunit,time_error) patch_rect,mptr,level,time,rnode(timemult,mptr),alphai,t_epsilon
                print time_error, patch_rect,mptr,level,time,rnode(timemult,mptr),alphai,t_epsilon
                call outtre(mstart,.false.,nvar,naux)
                stop
            endif

            ! Check if we should interpolate in time
            t_interpolate = .false.
            if (abs(alphai - 1.d0) < t_epsilon) then
                ! need no interpolation
                loc = node(store1,mptr)
            else if (dabs(alphai) .lt. t_epsilon) then
                loc = node(store2,mptr)
                if (level == mxnest) then
                    write(outunit,'(A,E24.16)') 'Line 95', dt
                    write(outunit,time_error) patch_rect,mptr,level,time, &
                                              rnode(timemult,mptr),alphai,t_epsilon
                    write(*,time_error) patch_rect,mptr,level,time, &
                                        rnode(timemult,mptr),alphai,t_epsilon
                    stop
                endif
            else
                locold  = node(store2,mptr)
                locnew  = node(store1,mptr)
                t_interpolate = .true.

                ! If we are at the maximum level nesting, abort
                if (level == mxnest) then
                    write(outunit,'(A,E24.16)') 'Line 107',dt
                    write(outunit,time_error) patch_rect,mptr,level,time,rnode(timemult,mptr),alphai,t_epsilon
                    print time_error, patch_rect,mptr,level,time,rnode(timemult,mptr),alphai,t_epsilon
                    stop
                endif
            endif

            ! Actual interpolation
            if (.not. t_interpolate) then
                ! No time interp. copy the solution values
                do ivar = 1, nvar
                    do j = jxlo, jxhi
                        do i = ixlo, ixhi
                            val(ivar,i-patch_rect(1)+nrowst,j-jlo+ncolst) = &
                                alloc(iadd(ivar,i-imlo+nghost+1,j-jmlo+nghost+1))
                            flaguse(i,j) = 1
                        end do
                    end do
                end do
            else
                ! Linear interpolation in time
                do ivar = 1, nvar
                    do j = jxlo, jxhi
                        do i = ixlo, ixhi
                            val(ivar,i-patch_rect(1)+nrowst,j-jlo+ncolst) = &
                                alloc(iadnew(ivar,i-imlo+nghost+1,j-jmlo+nghost+1))*alphai + &
                                alloc(iadold(ivar,i-imlo+nghost+1,j-jmlo+nghost+1))*alphac
                            flaguse(i,j) = 1
                        end do
                    end do
                end do
            endif

        endif

        ! Get next grid
!        mptr = node(levelptr, mptr)
         if (msrc .ne. -1) then
            nextSpot = bndList(nextSpot,nextfree)
         endif

    end do

    ! Set used array points which intersect domain boundary to be equal to 1,
    ! they will be set elsewhere, namely boundary conditions and periodic
    ! domains
    if (jhi >= jregsz(level)) then
        flaguse(patch_rect(1):ihi, max(jregsz(level),jlo):jhi) = 1
    endif

    if (jlo < 0) then
        flaguse(patch_rect(1):ihi, jlo:min(-1,ncolst + jhi - jlo )) = 1
    endif

    if (patch_rect(1) < 0) then
        flaguse(patch_rect(1):min(-1,nrowst + ihi - patch_rect(1)), jlo:jhi) = 1
    endif

    if (ihi >= iregsz(level)) then
        flaguse(max(iregsz(level),patch_rect(1)):ihi, jlo:jhi) = 1
    endif

contains

    integer pure function iadd(ivar,i,j)
        implicit none
        integer, intent(in) :: ivar, i, j
        iadd = loc + ivar-1 + nvar*((j-1)*mitot+i-1)
    end function iadd

    integer pure function iadnew(ivar,i,j)
        implicit none
        integer, intent(in) :: ivar, i, j
        iadnew = locnew + ivar-1 + nvar*((j-1)*mitot+i-1)
    end function iadnew

    integer pure function iadold(ivar,i,j)
        implicit none
        integer, intent(in) :: ivar, i, j
        iadold = locold + ivar-1 + nvar*((j-1)*mitot+i-1)
    end function iadold

end subroutine intfil

max_stars_repo_name: navravi/amrclaw

max_issues_repo_name: navravi/amrclaw

max_forks_repo_name: BrisaDavis/amrclaw