text:
!> @file
!! A simple Fortran example for the ParallelIO Library.
module pioExampleADIOS

    use pio, only : PIO_init, PIO_rearr_subset, iosystem_desc_t, file_desc_t
    use pio, only : PIO_finalize, PIO_noerr, PIO_iotype_adios, PIO_createfile
    use pio, only : PIO_int,var_desc_t, PIO_redef, PIO_def_dim, PIO_def_var, PIO_enddef
    use pio, only : PIO_closefile, io_desc_t, PIO_initdecomp, PIO_write_darray
    use pio, only : PIO_freedecomp, PIO_clobber, PIO_read_darray, PIO_syncfile, PIO_OFFSET_KIND
    use pio, only : PIO_nowrite, PIO_openfile, PIO_seterrorhandling, PIO_RETURN_ERROR
    use pio, only : PIO_GLOBAL, PIO_put_att, PIO_put_var

    implicit none

    private

    include 'mpif.h'

    !> @brief Length of the data array we are using.  This is then
    !! divided among MPI processes.
    integer, parameter :: LEN = 16

    !> @brief Value used for array that will be written to netcdf file.
    integer, parameter :: VAL = 42

    !> @brief Error code if anything goes wrong.
    integer, parameter :: ERR_CODE = 99

    integer, parameter :: d0 = 3;

    !> @brief A class to hold example code and data.
    !! This class contains the data and functions to execute the
    !! example.
    type, public :: pioExampleClass

        !> @brief Rank of processor running the code.
        integer :: myRank

        !> @brief Number of processors participating in MPI communicator.
        integer :: ntasks

        !> @brief Number of processors performing I/O.
        integer :: niotasks

        !> @brief Stride in the mpi rank between io tasks.
        integer :: stride

        !> @brief Number of aggregator.
        integer :: numAggregator

        !> @brief Start index of I/O processors.
        integer :: optBase

        !> @brief The ParallelIO system set up by @ref PIO_init.
        type(iosystem_desc_t) :: pioIoSystem

        !> @brief Contains data identifying the file.
        type(file_desc_t)     :: pioFileDesc

        !> @brief The netCDF variable ID.
        type(var_desc_t)      :: pioVar

        !> @brief The netCDF variable ID.
        type(var_desc_t)      :: pioVard0

        !> @brief An io descriptor handle that is generated in @ref PIO_initdecomp.
        type(io_desc_t)       :: iodescNCells

        !> @brief Specifies the flavor of netCDF output.
        integer               :: iotype

        !> @brief The netCDF dimension ID.
        integer               :: pioDimId

        !> @brief 1-based index of start of this processors data in full data array.
        integer               :: ista

        !> @brief Size of data array for this processor.
        integer               :: isto

        !> @brief Number of elements handled by each processor.
        integer               :: arrIdxPerPe

        !> @brief The length of the dimension of the netCDF variable.
        integer, dimension(1) :: dimLen

        !> @brief Buffer to hold sample data that is written to netCDF file.
        integer, allocatable  :: dataBuffer(:)

        !> @brief Buffer to read data into.
        integer, allocatable  :: readBuffer(:)

        !> @brief Array describing the decomposition of the data.
        integer, allocatable  :: compdof(:)

        !> @brief Name of the sample netCDF file written by this example.
        character(len=255) :: fileName

    contains

        !> @brief Initialize MPI, ParallelIO, and example data.
        !! Initialize the MPI and ParallelIO libraries. Also allocate
        !! memory to write and read the sample data to the netCDF file.
        procedure,  public  :: init

        !> @brief Create the decomposition for the example.
        !! This subroutine creates the decomposition for the example.
        procedure,  public  :: createDecomp

        !> @brief Create netCDF output file.
        !! This subroutine creates the netCDF output file for the example.
        procedure,  public  :: createFile

        !> @brief Define the netCDF metadata.
        !! This subroutine defines the netCDF dimension and variable used
        !! in the output file.
        procedure,  public  :: defineVar

        !> @brief Write the sample data to the output file.
        !! This subroutine writes the sample data array to the netCDF
        !! output file.
        procedure,  public  :: writeVar

        !> @brief Read the sample data from the output file.
        !! This subroutine reads the sample data array from the netCDF
        !! output file.
        procedure,  public  :: readVar

        !> @brief Close the netCDF output file.
        !! This subroutine closes the output file used by this example.
        procedure,  public  :: closeFile

        !> @brief Clean up resources.
        !! This subroutine cleans up resources used in the example. The
        !! ParallelIO and MPI libraries are finalized, and memory
        !! allocated in this example program is freed.
        procedure,  public  :: cleanUp

        !> @brief Handle errors.
        !! This subroutine is called if there is an error.
        procedure,  private :: errorHandle

    end type pioExampleClass

contains

    !> @brief Initialize MPI, ParallelIO, and example data.
    !! Initialize the MPI and ParallelIO libraries. Also allocate
    !! memory to write and read the sample data to the netCDF file.
    subroutine init(this)

        implicit none

        class(pioExampleClass), intent(inout) :: this

        integer :: ierr,i

        !
        ! initialize MPI
        !

        call MPI_Init(ierr)
        call MPI_Comm_rank(MPI_COMM_WORLD, this%myRank, ierr)
        call MPI_Comm_size(MPI_COMM_WORLD, this%ntasks , ierr)

        !
        ! set up PIO for rest of example
        !

        this%stride        = 1
        this%numAggregator = 0
        this%optBase       = 1
        this%iotype        = PIO_iotype_adios
        this%fileName      = "examplePio_f90.bp"
        this%dimLen(1)     = LEN

        this%niotasks = this%ntasks ! keep things simple - 1 iotask per MPI process

        !write(*,*) 'this%niotasks ',this%niotasks

        call PIO_init(this%myRank,      & ! MPI rank
            MPI_COMM_WORLD,             & ! MPI communicator
            this%niotasks,              & ! Number of iotasks (ntasks/stride)
            this%numAggregator,         & ! number of aggregators to use
            this%stride,                & ! stride
            PIO_rearr_subset,           & ! do not use any form of rearrangement
            this%pioIoSystem,           & ! iosystem
            base=this%optBase)            ! base (optional argument)

        !
        ! set up some data that we will write to a netcdf file
        !

        this%arrIdxPerPe = LEN / this%ntasks

        if (this%arrIdxPerPe < 1) then
            call this%errorHandle("Not enough work to distribute among pes", ERR_CODE)
        endif

        this%ista = this%myRank * this%arrIdxPerPe + 1
        this%isto = this%ista + (this%arrIdxPerPe - 1)

        allocate(this%compdof(this%ista:this%isto))
        allocate(this%dataBuffer(this%ista:this%isto))
        allocate(this%readBuffer(this%ista:this%isto))

        this%compdof(this%ista:this%isto) = (/(i, i=this%ista,this%isto, 1)/)
        this%dataBuffer(this%ista:this%isto) = this%myRank + VAL
        this%readBuffer(this%ista:this%isto) = 0

    end subroutine init

    subroutine createDecomp(this)

        implicit none

        class(pioExampleClass), intent(inout) :: this

        integer(PIO_OFFSET_KIND) :: start(1)
        integer(PIO_OFFSET_KIND) :: count(1)

        start(1) = this%ista
        count(1) = this%arrIdxPerPe

        call PIO_initdecomp(this%pioIoSystem, PIO_int, this%dimLen, this%compdof(this%ista:this%isto), &
            this%iodescNCells)

    end subroutine createDecomp

    subroutine createFile(this)

        implicit none

        class(pioExampleClass), intent(inout) :: this

        integer :: retVal, olderrorhandler

        retVal = PIO_createfile(this%pioIoSystem, this%pioFileDesc, this%iotype, trim(this%fileName), PIO_clobber)
        call PIO_seterrorhandling(this%pioFileDesc,PIO_RETURN_ERROR,olderrorhandler)
        call this%errorHandle("Could not create "//trim(this%fileName), retVal)

    end subroutine createFile

    subroutine defineVar(this)

        implicit none

        class(pioExampleClass), intent(inout) :: this

        integer :: retVal
        character(len=100) :: str       ! character temporary

        retVal = PIO_def_dim(this%pioFileDesc, 'x', this%dimLen(1) , this%pioDimId)
        call this%errorHandle("Could not define dimension x", retVal)

        retVal = PIO_def_var(this%pioFileDesc, 'foo', PIO_int, (/this%pioDimId/), this%pioVar)
        call this%errorHandle("Could not define variable foo", retVal)

        retVal = PIO_put_att(this%pioFileDesc, this%pioVar, 'startvalue', VAL)
        str = 'Just a string attribute'
        retVal = PIO_put_att(this%pioFileDesc, PIO_GLOBAL, 'astringattr', str)

        retVal = PIO_def_var(this%pioFileDesc, 'd0', PIO_int, this%pioVard0)

        retVal = PIO_enddef(this%pioFileDesc)
        call this%errorHandle("Could not end define mode", retVal)

    end subroutine defineVar

    subroutine writeVar(this)

    implicit none

        class(pioExampleClass), intent(inout) :: this

        integer :: retVal

        call PIO_write_darray(this%pioFileDesc, this%pioVar, this%iodescNCells, this%dataBuffer(this%ista:this%isto), retVal)
        call this%errorHandle("Could not write foo", retVal)
        retVal = PIO_put_var(this%pioFileDesc, this%pioVard0, d0);
        call PIO_syncfile(this%pioFileDesc)

    end subroutine writeVar

    subroutine readVar(this)

        implicit none

        class(pioExampleClass), intent(inout) :: this

        integer :: retVal

        call PIO_read_darray(this%pioFileDesc, this%pioVar, this%iodescNCells,  this%readBuffer, retVal)
        call this%errorHandle("Could not read foo", retVal)

    end subroutine readVar

    subroutine closeFile(this)

        implicit none

        class(pioExampleClass), intent(inout) :: this

        call PIO_closefile(this%pioFileDesc)

    end subroutine closeFile

    subroutine cleanUp(this)

        implicit none

        class(pioExampleClass), intent(inout) :: this

        integer :: ierr

        deallocate(this%compdof)
        deallocate(this%dataBuffer)
        deallocate(this%readBuffer)

        call PIO_freedecomp(this%pioIoSystem, this%iodescNCells)
        call PIO_finalize(this%pioIoSystem, ierr)
        call MPI_Finalize(ierr)

    end subroutine cleanUp

    subroutine errorHandle(this, errMsg, retVal)

        implicit none

        class(pioExampleClass), intent(inout) :: this
        character(len=*),       intent(in)    :: errMsg
        integer,                intent(in)    :: retVal

        if (retVal .ne. PIO_NOERR) then
            write(*,*) retVal,errMsg
            call PIO_closefile(this%pioFileDesc)
            call mpi_abort(MPI_COMM_WORLD,0,retVal)
        end if

    end subroutine errorHandle

end module pioExampleADIOS

!> @brief Main execution of example code.
!! This is an example program for the ParallelIO library.
!!
!! This program creates a netCDF output file with the ParallelIO
!! library, then writes and reads some data to and from the file.
!!
!! This example does the following:
!!
!! - initialization initializes the MPI library, initializes the
!!   ParallelIO library with @ref PIO_init. Then allocate memory for a
!!   data array of sample data to write, and an array to read the data
!!   back into. Also allocate an array to hold decomposition
!!   information.
!!
!! - creation of decomposition by calling @ref PIO_initdecomp.
!!
!! - creation of netCDF file with @ref PIO_createfile.
!!
!! - define netCDF metadata with @ref PIO_def_dim and @ref
!!   PIO_def_var. Then end define mode with @ref PIO_enddef.
!!
!! - write the sample data with @ref PIO_write_darray. Then sync the
!!   file with @ref PIO_syncfile.
!!
!! - read the sample data with @ref PIO_read_darray.
!!
!! - close the netCDF file with @ref PIO_closefile.
!!
!! - clean up local memory, ParallelIO library resources with @ref
!!   PIO_freedecomp and @ref PIO_finalize, and MPI library resources.
!!
program main

    use pioExampleADIOS, only : pioExampleClass

    implicit none

    type(pioExampleClass) :: pioExInst

    call pioExInst%init()
    call pioExInst%createDecomp()
    call pioExInst%createFile()
    call pioExInst%defineVar()
    call pioExInst%writeVar()
    !call pioExInst%readVar()
    call pioExInst%closeFile()
    call pioExInst%cleanUp()

end program main

max_stars_repo_licenses: FTL zlib-acknowledgement RSA-MD

max_issues_repo_licenses: FTL zlib-acknowledgement RSA-MD

max_forks_repo_licenses: FTL zlib-acknowledgement RSA-MD


text:
module mod_determinant
#include <messenger.h>
    use mod_kinds,      only: rk,ik,rdouble,rsingle
    use mod_constants,  only: ONE
    implicit none

contains

    !>  Compute and return the determinant of a 3x3 matrix.
    !!
    !!  @author Nathan A. Wukie (AFRL)
    !!  @date   8/14/2017
    !!
    !!
    !----------------------------------------------------------
    function det_3x3(mat) result(det)
        real(rk),   intent(in)  :: mat(3,3)

        real(rk)    :: det

        det = mat(1,1)*mat(2,2)*mat(3,3) - mat(1,2)*mat(2,1)*mat(3,3) - &
              mat(1,1)*mat(2,3)*mat(3,2) + mat(1,3)*mat(2,1)*mat(3,2) + &
              mat(1,2)*mat(2,3)*mat(3,1) - mat(1,3)*mat(2,2)*mat(3,1)

    end function det_3x3
    !**********************************************************

    !>  Compute and return the derivative of determinant of a 3x3 matrix
    !!  wrt a given differentiated matrix
    !!
    !!  \partial det(mat) / \partial x
    !!
    !!  dmat = dmat/dx
    !!
    !!
    !!  @author Matteo Ugolotti
    !!  @date   7/17/2018
    !!
    !!
    !----------------------------------------------------------
    function ddet_3x3(mat,dmat) result(det)
        real(rk),   intent(in)  :: mat(3,3)
        real(rk),   intent(in)  :: dmat(3,3)

        real(rk)    :: det1, det2, det
        real(rk)    :: adjugate(3,3), cofactor(3,3), temp(3,3)

        !
        ! Method 1 and 2 give the same result, 1 is less expensive
        !

        !
        ! Method 1: Direct differentiation.
        !
        !det1 =  dmat(1,1)*mat(2,2)*mat(3,3) + mat(1,1)*dmat(2,2)*mat(3,3) + mat(1,1)*mat(2,2)*dmat(3,3) &
        !      - dmat(1,2)*mat(2,1)*mat(3,3) - mat(1,2)*dmat(2,1)*mat(3,3) - mat(1,2)*mat(2,1)*dmat(3,3) &
        !      - dmat(1,1)*mat(2,3)*mat(3,2) - mat(1,1)*dmat(2,3)*mat(3,2) - mat(1,1)*mat(2,3)*dmat(3,2) &
        !      + dmat(1,3)*mat(2,1)*mat(3,2) + mat(1,3)*dmat(2,1)*mat(3,2) + mat(1,3)*mat(2,1)*dmat(3,2) &
        !      + dmat(1,2)*mat(2,3)*mat(3,1) + mat(1,2)*dmat(2,3)*mat(3,1) + mat(1,2)*mat(2,3)*dmat(3,1) &
        !      - dmat(1,3)*mat(2,2)*mat(3,1) - mat(1,3)*dmat(2,2)*mat(3,1) - mat(1,3)*mat(2,2)*dmat(3,1)

        !
        ! Method 2: Jacobi's formula
        !   ref: https://en.wikipedia.org/wiki/Jacobi%27s_formula
        !

        ! Find Cofactor matrix of mat
        cofactor(1,1) =  mat(2,2)*mat(3,3)-mat(2,3)*mat(3,2)
        cofactor(1,2) = -mat(2,1)*mat(3,3)+mat(2,3)*mat(3,1)
        cofactor(1,3) =  mat(2,1)*mat(3,2)-mat(2,2)*mat(3,1)
        cofactor(2,1) = -mat(1,2)*mat(3,3)+mat(1,3)*mat(3,2)
        cofactor(2,2) =  mat(1,1)*mat(3,3)-mat(1,3)*mat(3,1)
        cofactor(2,3) = -mat(1,1)*mat(3,2)+mat(1,2)*mat(3,1)
        cofactor(3,1) =  mat(1,2)*mat(2,3)-mat(1,3)*mat(2,2)
        cofactor(3,2) = -mat(1,1)*mat(2,3)+mat(1,3)*mat(2,1)
        cofactor(3,3) =  mat(1,1)*mat(2,2)-mat(1,2)*mat(2,1)

        ! Adjugate of mat is the transpose of cofactor of mat
        adjugate = transpose(cofactor)

        ! Use Jacobi's formula to find the derivative of the determinant
        temp = matmul(adjugate,dmat)

        ! Find the trace of temp
        det2 = temp(1,1)+temp(2,2)+temp(3,3)

        ! Result
        det = det2

    end function ddet_3x3
    !**********************************************************

end module mod_determinant

max_stars_repo_name: nwukie/ChiDG

max_issues_repo_name: nwukie/ChiDG

max_forks_repo_name: nwukie/ChiDG


text:
      DOUBLE PRECISION FUNCTION D1MACH(I)
C
C  Double-precision machine constants.
C  This implementation for use in MATLAB Fortran Mex-files.
C
C  D1MACH(1) = realmin = B**(EMIN-1), the smallest positive magnitude.
C  D1MACH(2) = realmax = B**EMAX*(1 - B**(-T)), the largest magnitude.
C  D1MACH(3) = eps/2 = B**(-T), the smallest relative spacing.
C  D1MACH(4) = eps = B**(1-T), the largest relative spacing.
C  D1MACH(5) = LOG10(B)
C
c  DLAMCH
*          = 'E' or 'e',   DLAMCH := eps
*          = 'S' or 's ,   DLAMCH := sfmin
*          = 'B' or 'b',   DLAMCH := base
*          = 'P' or 'p',   DLAMCH := eps*base
*          = 'N' or 'n',   DLAMCH := t
*          = 'R' or 'r',   DLAMCH := rnd
*          = 'M' or 'm',   DLAMCH := emin
*          = 'U' or 'u',   DLAMCH := rmin
*          = 'L' or 'l',   DLAMCH := emax
*          = 'O' or 'o',   DLAMCH := rmax
*
*          where
*
*          eps   = relative machine precision
*          sfmin = safe minimum, such that 1/sfmin does not overflow
*          base  = base of the machine
*          prec  = eps*base
*          t     = number of (base) digits in the mantissa
*          rnd   = 1.0 when rounding occurs in addition, 0.0 otherwise
*          emin  = minimum exponent before (gradual) underflow
*          rmin  = underflow threshold - base**(emin-1)
*          emax  = largest exponent before overflow
*          rmax  = overflow threshold  - (base**emax)*(1-eps)
*

      double precision DLAMCH
      IF (I .EQ. 1) D1MACH = DLAMCH('U')
      IF (I .EQ. 2) D1MACH = DLAMCH('O')
      IF (I .EQ. 3) D1MACH = DLAMCH('E')
      IF (I .EQ. 4) D1MACH = DLAMCH('P')
      IF (I .EQ. 5) D1MACH = log10(DLAMCH('B'))
      RETURN
      END

ext: f

max_stars_repo_name: alpgodev/numx

max_issues_repo_name: alpgodev/numx

max_forks_repo_name: alpgodev/numx


text:
module parameters

 !Number of x, y and z grid divisions:
integer,parameter:: nx=160, ny=80, nz=40

 !Maximum number of parcels:
integer,parameter:: nm=64*nx*ny*nz

 !Box width, breadth and height:
double precision,parameter:: ellx=4.d0, elly=2.d0, ellz=1.d0

end module parameters

max_stars_repo_name: matt-frey/epic

max_issues_repo_name: matt-frey/epic

max_forks_repo_name: matt-frey/epic


text:
program preguntar_nombre
  implicit none

  integer, parameter :: MAXLEN = 250
  character(len=MAXLEN), dimension(:), allocatable, target :: lista_nom, aux_arr
  !character(len=MAXLEN), pointer :: lista_ptr(:)

  character :: c_input = ''
  character(len=MAXLEN) :: s_input = ""
  integer :: i_input = 0, i = 0, LIST_SIZE = 0

  allocate(lista_nom(1:4))

  lista_nom(1) = 'Nicky'
  lista_nom(2) = 'Jorge'
  lista_nom(3) = 'Luis'
  lista_nom(4) = 'Alejandra'

  LIST_SIZE = size(lista_nom)

  print *, "Crear nueva lista? (s/n)"
  read *, c_input

  if (c_input.eq.'s') then
      print*, "Okay, utilizando nueva lista :)"

      do i = 1, size(lista_nom)
        lista_nom(i) = ""
      end do

      deallocate(lista_nom)

      print*, "Cuántos elementos?"
      read*, i_input

      allocate(lista_nom(1:i_input))

      do i=1, i_input
        print *, "Ingrese el elemento ", i
        read*, lista_nom(i)
      end do

      print*, "Desea ingresar más datos? (s/n)"
      c_input = ''
      read*, c_input

      if (c_input.eq.'s') then
        print*, "Okay!, escribe porfavor 'fin' cuando hayas terminado"
        do while(.true.)
          LIST_SIZE = size(lista_nom) + 1
          allocate(aux_arr(LIST_SIZE))

          print*, "Escribe el elemento ", LIST_SIZE
          s_input = ""
          aux_arr = lista_nom
          deallocate(lista_nom)
          allocate(lista_nom(LIST_SIZE))
          lista_nom = aux_arr

          do i=1, size(aux_arr)
            lista_nom(i) = aux_arr(i)
          end do

          print*, "[DEBUG] lista_nom: ", lista_nom

          read*, s_input

          if (s_input.eq."fin") then
            exit
          else
            lista_nom(LIST_SIZE) = s_input
          end if

          deallocate(aux_arr)

        end do
      end if
    else
      print*, "Okay, utilizando lista predeterminada :)"
  end if

  print *, "Imprimiendo tu lista:"
  print*

  do i = 1, size(lista_nom)
    print*,lista_nom(i)

  end do

  deallocate(lista_nom)

end program preguntar_nombre

max_stars_repo_name: kotoromo/Fortran-Basico

max_issues_repo_name: kotoromo/Fortran-Basico

max_forks_repo_name: kotoromo/Fortran-Basico


text:
subroutine ca_filt_prim(lo, hi, &
     Stmp, Stmp_l1, Stmp_h1, &
     Snew, Snew_l1, Snew_h1, &
     mask, mask_l1, mask_h1, &
     filt_T, S, domlo,domhi, &
     delta,xlo,problo,time,level) bind(C, name="ca_filt_prim")

  use probdata_module, only : filter_rhomax, filter_timemax
  use network, only : naux, nspec
  use eos_module, only : eos
  use eos_type_module, only : eos_t, eos_input_rt
  use meth_params_module, only : NVAR, URHO, UMX, UEDEN, UEINT, UTEMP, UFA, UFS, UFX, &
       small_temp, small_dens, nadv
  use filter_module

  use amrex_fort_module, only : rt => amrex_real
  implicit none

  integer, intent(in) :: lo(1), hi(1), domlo(1), domhi(1), level
  integer, intent(in) :: filt_T, S
  real(rt)        , intent(in) :: delta(1), xlo(1), problo(1), time
  integer, intent(in) ::   Stmp_l1,Stmp_h1
  integer, intent(in) ::   Snew_l1,Snew_h1
  integer, intent(in) ::   mask_l1,mask_h1
  real(rt)         :: Stmp(Stmp_l1:Stmp_h1,NVAR)
  real(rt)         :: Snew(Snew_l1:Snew_h1,NVAR)
  real(rt)         :: mask(mask_l1:mask_h1)
  ! mask has three possible values: -1.e0_rt, 0.e0_rt, and 1.e0_rt.
  ! -1.e0_rt appears only in cells that are covered by neither this level nor the finer level.
  !       It can only appear in ghost cells.
  !  0.e0_rt appears only in cells that are covered by only this level, but not the finer level.
  !  1.e0_rt appears only in cells that are covered by the finer level.
  !       It can appear in either valid cells or ghost cells.

  integer :: i
  real(rt)         :: rhotmpinv
  type(eos_t) :: eos_state
  logical, allocatable :: filtered(:)

  allocate(filtered(lo(1):hi(1)))
  filtered = .false.

  if (filt_T .eq. 1) then

     do i=lo(1), hi(1)
        if (time < filter_timemax .or. Snew(i,URHO) < filter_rhomax) then
           filtered(i) = .true.
           Snew(i,URHO ) = ff1(0) *  Stmp(i  ,URHO ) &
                &        + ff1(1) * (Stmp(i-1,URHO ) + Stmp(i+1,URHO ))
           Snew(i,UMX  ) = ff1(0) *  Stmp(i  ,UMX  ) &
                &        + ff1(1) * (Stmp(i-1,UMX  ) + Stmp(i+1,UMX  ))
           Snew(i,UTEMP) = ff1(0) *  Stmp(i  ,UTEMP) &
                &        + ff1(1) * (Stmp(i-1,UTEMP) + Stmp(i+1,UTEMP))
        end if
     end do

  else if (filt_T .eq. 2) then

     do i=lo(1), hi(1)
        if (time < filter_timemax .or. Snew(i,URHO) < filter_rhomax) then
           filtered(i) = .true.
           Snew(i,URHO ) = ff2(0,S) *  Stmp(i  ,URHO ) &
                &        + ff2(1,S) * (Stmp(i-1,URHO ) + Stmp(i+1,URHO )) &
                &        + ff2(2,S) * (Stmp(i-2,URHO ) + Stmp(i+2,URHO ))
           Snew(i,UMX  ) = ff2(0,S) *  Stmp(i  ,UMX  ) &
                &        + ff2(1,S) * (Stmp(i-1,UMX  ) + Stmp(i+1,UMX  )) &
                &        + ff2(2,S) * (Stmp(i-2,UMX  ) + Stmp(i+2,UMX  ))
           Snew(i,UTEMP) = ff2(0,S) *  Stmp(i  ,UTEMP) &
                &        + ff2(1,S) * (Stmp(i-1,UTEMP) + Stmp(i+1,UTEMP)) &
                &        + ff2(2,S) * (Stmp(i-2,UTEMP) + Stmp(i+2,UTEMP))
        end if
     end do

  else if (filt_T .eq. 3) then

     do i=lo(1), hi(1)
        if (time < filter_timemax .or. Snew(i,URHO) < filter_rhomax) then
           filtered(i) = .true.
           Snew(i,URHO ) = ff3(0,S) *  Stmp(i  ,URHO ) &
                &        + ff3(1,S) * (Stmp(i-1,URHO ) + Stmp(i+1,URHO )) &
                &        + ff3(2,S) * (Stmp(i-2,URHO ) + Stmp(i+2,URHO )) &
                &        + ff3(3,S) * (Stmp(i-3,URHO ) + Stmp(i+3,URHO ))
           Snew(i,UMX  ) = ff3(0,S) *  Stmp(i  ,UMX  ) &
                &        + ff3(1,S) * (Stmp(i-1,UMX  ) + Stmp(i+1,UMX  )) &
                &        + ff3(2,S) * (Stmp(i-2,UMX  ) + Stmp(i+2,UMX  )) &
                &        + ff3(3,S) * (Stmp(i-3,UMX  ) + Stmp(i+3,UMX  ))
           Snew(i,UTEMP) = ff3(0,S) *  Stmp(i  ,UTEMP) &
                &        + ff3(1,S) * (Stmp(i-1,UTEMP) + Stmp(i+1,UTEMP)) &
                &        + ff3(2,S) * (Stmp(i-2,UTEMP) + Stmp(i+2,UTEMP)) &
                &        + ff3(3,S) * (Stmp(i-3,UTEMP) + Stmp(i+3,UTEMP))
        end if
     end do

  else if (filt_T .eq. 4) then

     do i=lo(1), hi(1)
        if (time < filter_timemax .or. Snew(i,URHO) < filter_rhomax) then
           filtered(i) = .true.
           Snew(i,URHO ) = ff4(0,S) *  Stmp(i  ,URHO ) &
                &        + ff4(1,S) * (Stmp(i-1,URHO ) + Stmp(i+1,URHO )) &
                &        + ff4(2,S) * (Stmp(i-2,URHO ) + Stmp(i+2,URHO )) &
                &        + ff4(3,S) * (Stmp(i-3,URHO ) + Stmp(i+3,URHO )) &
                &        + ff4(4,S) * (Stmp(i-4,URHO ) + Stmp(i+4,URHO ))
           Snew(i,UMX  ) = ff4(0,S) *  Stmp(i  ,UMX  ) &
                &        + ff4(1,S) * (Stmp(i-1,UMX  ) + Stmp(i+1,UMX  )) &
                &        + ff4(2,S) * (Stmp(i-2,UMX  ) + Stmp(i+2,UMX  )) &
                &        + ff4(3,S) * (Stmp(i-3,UMX  ) + Stmp(i+3,UMX  )) &
                &        + ff4(4,S) * (Stmp(i-4,UMX  ) + Stmp(i+4,UMX  ))
           Snew(i,UTEMP) = ff4(0,S) *  Stmp(i  ,UTEMP) &
                &        + ff4(1,S) * (Stmp(i-1,UTEMP) + Stmp(i+1,UTEMP)) &
                &        + ff4(2,S) * (Stmp(i-2,UTEMP) + Stmp(i+2,UTEMP)) &
                &        + ff4(3,S) * (Stmp(i-3,UTEMP) + Stmp(i+3,UTEMP)) &
                &        + ff4(4,S) * (Stmp(i-4,UTEMP) + Stmp(i+4,UTEMP))
        end if
     end do

  end if

  do i=lo(1), hi(1)
     if (filtered(i)) then

        Snew(i,URHO ) = max(Snew(i,URHO ), small_dens)
        Snew(i,UTEMP) = max(Snew(i,UTEMP), small_temp)

        eos_state % rho = Snew(i,URHO)
        eos_state % T = Snew(i,UTEMP)

        rhotmpInv = 1.e0_rt / Stmp(i,URHO)
        eos_state % xn  = Stmp(i,UFS:UFS+nspec-1) * rhotmpInv
        eos_state % aux = Stmp(i,UFX:UFX+naux-1) * rhotmpInv

        call eos(eos_input_rt, eos_state)

        Snew(i,UEINT) = Snew(i,URHO) * eos_state % e
        Snew(i,UEDEN) = Snew(i,UEINT) + &
             0.5*(Snew(i,UMX)**2)/Snew(i,URHO)
        if (nadv > 0) then
           Snew(i,UFA:UFA-1+nadv) = Stmp(i,UFA:UFA-1+nadv) * rhotmpinv * Snew(i,URHO)
        end if
        Snew(i,UFS:UFS-1+nspec) = eos_state % xn * Snew(i,URHO)
        if (naux > 0) then
           Snew(i,UFX:UFX-1+naux) = eos_state % aux * Snew(i,URHO)
        end if

     end if
  end do

  deallocate(filtered)

end subroutine ca_filt_prim

max_stars_repo_name: abigailbishop/Castro

max_issues_repo_name: abigailbishop/Castro

max_forks_repo_name: abigailbishop/Castro


text:
C------------------------------------------------------------------
C                        VICAR SUBROUTINE pu75
C
C Edit History:
C      4-94        CRI    MSTP S/W CONVERSION (VICAR PORTING)
C------------------------------------------------------------------

      SUBROUTINE PU75(U,N)

      REAL*4 U(N,2,103)
      CHARACTER*132 RBUF
      CHARACTER*132 QBUF
      CHARACTER*132 PBUF

C  U IS REAL*4 IF N=1, REAL*8 IF N=2
C
      J = 1
C
      DO II = 1, 5
         PBUF(1:132) = ' '
         QBUF(1:132) = ' '
         RBUF(1:132) = ' '
         K = 7
C
         DO I = 1, 11
            WRITE (RBUF(K-2-2:K-2),'(I3)') J
            WRITE (PBUF(K-5:K),'(F6.1)') U(1,1,J)
            WRITE (QBUF(K-5:K),'(F6.1)') U(1,2,J)
            J = J + 1
            K = K + 12
         end do
C
         CALL XVMESSAGE(RBUF(2:132),' ')
         CALL XVMESSAGE(PBUF(2:132),' ')
         CALL XVMESSAGE(QBUF(2:132),' ')
         CALL XVMESSAGE(' ',' ')
         IF (II .EQ. 5) RETURN
         PBUF(1:132) = ' '
         QBUF(1:132) = ' '
         RBUF(1:132) = ' '
         K = 7
C
         DO I = 1, 12
            WRITE (RBUF(K-2-2:K-2),'(I3)') J
            WRITE (PBUF(K-5:K),'(F6.1)') U(1,1,J)
            WRITE (QBUF(K-5:K),'(F6.1)') U(1,2,J)
            J = J + 1
            IF (I .EQ. 1) K = K - 6
            IF (I .EQ. 11) K = K - 6
            K = K + 12
         END DO
C
         CALL XVMESSAGE(RBUF(2:132),' ')
         CALL XVMESSAGE(PBUF(2:132),' ')
         CALL XVMESSAGE(QBUF(2:132),' ')
         CALL XVMESSAGE(' ',' ')
      END DO
C
      return
      END

max_stars_repo_name: NASA-AMMOS/VICAR

max_issues_repo_name: NASA-AMMOS/VICAR

max_forks_repo_name: NASA-AMMOS/VICAR


text: The fountain in the Memorial Union MU Patio patio is the http://wwwdateline.ucdavis.edu/030896/30896no2.html Students First recognition Water Fountains fountain. Donors who contributed more than a thousand dollars toward the Students First campaign were recognized with a tile.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
!
! Spline_Interpolation
!
! Module containing spline interpolation routines
!
!
! CREATION HISTORY:
!       Written by:     Paul van Delst, CIMSS/SSEC 06-Oct-2006
!                       paul.vandelst@ssec.wisc.edu
!

MODULE Spline_Interpolation

  ! ------------------
  ! Environment set up
  ! ------------------
  ! Modules used
  USE Type_Kinds,      ONLY: fp=>fp_kind
  USE Message_Handler, ONLY: SUCCESS, FAILURE, Display_Message
  USE Search_Utility,  ONLY: Value_Locate, Bisection_Search
  ! Disable all implicit typing
  IMPLICIT NONE

  ! ------------
  ! Visibilities
  ! ------------
  PRIVATE
  PUBLIC :: Spline_Interpolate
  PUBLIC :: Spline_Initialize

  ! ---------------------
  ! Procedure overloading
  ! ---------------------
  INTERFACE Spline_Interpolate
    MODULE PROCEDURE splint_scalar
    MODULE PROCEDURE splint_rank1
  END INTERFACE Spline_Interpolate

  ! -----------------
  ! Module parameters
  ! -----------------
  ! RCS Id field
  CHARACTER(*), PRIVATE, PARAMETER :: MODULE_RCS_ID = &
  ! Literal constants
  REAL(fp), PARAMETER :: ZERO = 0.0_fp
  REAL(fp), PARAMETER :: ONE  = 1.0_fp
  REAL(fp), PARAMETER :: TWO  = 2.0_fp
  REAL(fp), PARAMETER :: SIX  = 6.0_fp
  REAL(fp), PARAMETER :: TOLERANCE = EPSILON( ZERO )

CONTAINS

!------------------------------------------------------------------------------
!
! NAME:
!       Spline_Initialize
!
! PURPOSE:
!       Function is called to establish the type of interpolating
!       spline for a tabulated set of functional values,
!           Xi, Yi = F(Xi)
!
! CALLING SEQUENCE:
!       Error_Status = Spline_Initialize( x,                      &  ! Input
!                                         y,                      &  ! Input
!                                         y2,                     &  ! Output
!                                         RCS_Id     =RCS_Id,     &  ! Optional output
!                                         Message_Log=Message_Log )  ! Error messaging
! INPUT ARGUMENTS:
!       x:            The abscissa values for the tabulated function, y=F(x).
!                     Must be monotonically ascending values.
!                     UNITS:      Data dependent.
!                     TYPE:       REAL(fp)
!                     DIMENSION:  Rank-1 array
!                     ATTRIBUTES: INTENT(IN)
!
!       y:            An input vector of tabulated values, y=F(x), for which
!                     interpolates are required. Must have the same number of
!                     elements as the X input argument.
!                     UNITS:      Data dependent.
!                     TYPE:       REAL(fp)
!                     DIMENSION:  Rank-1 array
!                     ATTRIBUTES: INTENT(IN)
!
! OPTIONAL INPUT ARGUMENTS:
!       Message_Log:  Character string specifying a filename in which any
!                     messages will be logged. If not specified, or if an
!                     error occurs opening the log file, the default action
!                     is to output messages to the screen.
!                     UNITS:      N/A
!                     TYPE:       CHARACTER(*)
!                     DIMENSION:  Scalar
!                     ATTRIBUTES: OPTIONAL, INTENT(IN)
!
! OUTPUT ARGUMENTS:
!       y2:           Values of the second derivative of the interpolating
!                     function at the points defined by X. Must have the same
!                     number of elements as the X/Y input arguments.
!                     UNITS:      Data dependent.
!                     TYPE:       REAL(fp)
!                     DIMENSION:  Rank-1 array
!                     ATTRIBUTES: INTENT(IN)
!
! OPTIONAL OUTPUT ARGUMENTS:
!       RCS_Id:       Character string containing the Revision Control System
!                     Id field for the module.
!                     UNITS:      N/A
!                     TYPE:       CHARACTER(*)
!                     DIMENSION:  Scalar
!                     ATTRIBUTES: OPTIONAL, INTENT(OUT)
!
! FUNCTION RESULT:
!       Error_Status: The return value is an integer defining the error status.
!                     The error codes are defined in the Message_Handler module.
!                     If == SUCCESS the function executed normally.
!                        == FAILURE an unrecoverable error occurred.
!                     UNITS:      N/A
!                     TYPE:       INTEGER
!                     DIMENSION:  Scalar
!
! RESTRICTIONS:
!       Input X vectors *must* be monotonically ascending.
!
! PROCEDURE:
!       Adapted from SPL_INIT in "Numerical Recipes in FORTRAN; The Art of
!       Scientific Computing" 2nd ed. by W.H. Press et al., Cambridge University
!       Press, 1992, New York, New York, pp109.
!
!       SPLINE_INITIALIZE should be called only *once* to process an entire
!       tabulated function in the input arrays X and Y. Once this has been
!       done, values of the interpolated function for any value of X can be
!       obtained by calls to the separate routine SPLINE_INTERPOLATE.
!
! CREATION HISTORY:
!       Written by:   Paul van Delst, CIMSS/SSEC, 22-Mar-1996
!                     paul.vandelst@ssec.wisc.edu
!
!------------------------------------------------------------------------------

  FUNCTION Spline_Initialize( x,            &  ! Input
                              y,            &  ! Input
                              y2,           &  ! Output
                              RCS_Id,       &  ! Optional output
                              Message_Log ) &  ! Error messaging
                            RESULT ( Error_Status )
    ! Arguments
    REAL(fp), DIMENSION(:), INTENT(IN)  :: x
    REAL(fp), DIMENSION(:), INTENT(IN)  :: y
    REAL(fp), DIMENSION(:), INTENT(OUT) :: y2
    CHARACTER(*), OPTIONAL, INTENT(OUT) :: RCS_Id
    CHARACTER(*), OPTIONAL, INTENT(IN)  :: Message_Log
    ! Function result
    INTEGER :: Error_Status
    ! Local parameters
    CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'Spline_Initialize'
    ! Local variables
    INTEGER :: i
    INTEGER :: n
    REAL(fp) :: p
    REAL(fp) :: qn
    REAL(fp) :: sig
    REAL(fp) :: un
    REAL(fp) :: numer1
    REAL(fp) :: numer2
    REAL(fp) :: numerator
    REAL(fp) :: denominator
    REAL(fp) :: offset
    REAL(fp), DIMENSION(SIZE(x)) :: u

    ! ------
    ! Set up
    ! ------
    Error_Status = SUCCESS
    IF ( PRESENT( RCS_Id ) ) RCS_Id = MODULE_RCS_Id

    ! Are vector sizes consistent?
    n = SIZE(x)
    IF ( SIZE(y) /= n ) THEN
      Error_Status = FAILURE
      CALL Display_Message( ROUTINE_NAME, &
                            'Input X, Y arguments are different sizes.', &
                            Error_Status, &
                            Message_Log=Message_Log )
      RETURN
    END IF

    IF ( SIZE(y2) /= n ) THEN
      Error_Status = FAILURE
      CALL Display_Message( ROUTINE_NAME, &
                            'Output Y2 argument has inconsistent size.', &
                            Error_Status, &
                            Message_Log=Message_Log )
      RETURN
    END IF

    ! Is X-data in ascending order?
    u(1:n-1) = x(2:n) - x(1:n-1)
    IF ( ANY( u(1:n-1) < TOLERANCE ) ) THEN
      Error_Status = FAILURE
      CALL Display_Message( ROUTINE_NAME, &
                            'Input X data must be in ascending order', &
                            Error_Status, &
                            Message_Log=Message_Log )
      RETURN
    END IF

    ! ---------------------------------------
    ! Set second derivative at lower boundary
    ! ---------------------------------------
    y2(1) = ZERO
    u(1)  = ZERO

    ! -------------------------------------------------
    ! Begin decomposition loop of tridiagonal algorithm
    ! -------------------------------------------------
    DO i = 2, n - 1
      sig     = ( x(i) - x(i-1) ) / ( x(i+1) - x(i-1) )
      p       = ( sig * y2( i-1 ) ) + TWO
      y2( i ) = ( sig - ONE ) / p

      numer1    = ( y(i+1) - y(i) ) / ( x(i+1) - x(i) )
      numer2    = ( y(i) - y(i-1) ) / ( x(i) - x(i-1) )
      numerator   = SIX * ( numer1 - numer2 )
      denominator = x(i+1) - x(i-1)
      offset      = sig * u(i-1)

      u( i ) = ( ( numerator / denominator ) - offset ) / p
    END DO

    ! ---------------------------------------
    ! Set second derivative at upper boundary
    ! ---------------------------------------
    qn = ZERO
    un = ZERO
    y2( n ) = ( un - ( qn * u(n-1) ) ) / ( ( qn * y2(n-1) ) + ONE )

    ! ----------------------------------------------
    ! Backsubstitution loop of tridiagonal algorithm
    ! ----------------------------------------------
    DO i = n - 1, 1, -1
      y2( i ) = ( y2(i) * y2(i+1) ) + u(i)
    END DO

  END FUNCTION Spline_Initialize

!--------------------------------------------------------------------------------
!
! NAME:
!       Spline_Interpolate
!
! PURPOSE:
!       Function to calculate the cubic-spline interpolated value of a tabulated
!       function, X, Y=F(X), at a given value of Xi.
!
! CALLING SEQUENCE:
!       Error_Status = Spline_Interpolate( x,                      &  ! Input
!                                          y,                      &  ! Input
!                                          x_int,                  &  ! Input
!                                          y_int,                  &  ! Output
!                                          y2         =y2,         &  ! Optional input
!                                          RCS_Id     =RCS_Id,     &  ! Optional output
!                                          Message_Log=Message_Log )  ! Error messaging
!
! INPUTS:
!       x:            The abscissa values for the tabulated function, y=F(x).
!                     Must be monotonically ascending values.
!                     UNITS:      Data dependent.
!                     TYPE:       REAL(fp)
!                     DIMENSION:  Rank-1 array
!                     ATTRIBUTES: INTENT(IN)
!
!       y:            An input vector of tabulated values, y=F(x), for which
!                     interpolates are required. Must have the same number of
!                     elements as the X input argument.
!                     UNITS:      Data dependent.
!                     TYPE:       REAL(fp)
!                     DIMENSION:  Rank-1 array
!                     ATTRIBUTES: INTENT(IN)
!
!       x_int:        Abscissa values for which interpolates are required.
!                     Must be monotonically ascending or descending in the
!                     same direction as the X input argument.
!                     UNITS:      Data dependent, same as x.
!                     TYPE:       REAL(fp)
!                     DIMENSION:  Scalar or Rank-1 array
!                     ATTRIBUTES: INTENT(IN)
!
! OPTIONAL INPUTS:
!       y2:           Values of the second derivative of the interpolating
!                     function at the points defined by X. Must have the same
!                     number of elements as the X/Y input arguments.
!                     If this argument is not supplied, it is obtained via a
!                     call to SPLINE_INITIALIZE.
!                     UNITS:      Data dependent.
!                     TYPE:       REAL(fp)
!                     DIMENSION:  Rank-1 array
!                     ATTRIBUTES: INTENT(IN)
!
!       Message_Log:  Character string specifying a filename in which any
!                     messages will be logged. If not specified, or if an
!                     error occurs opening the log file, the default action
!                     is to output messages to the screen.
!                     UNITS:      N/A
!                     TYPE:       CHARACTER(*)
!                     DIMENSION:  Scalar
!                     ATTRIBUTES: OPTIONAL, INTENT(IN)
!
! OUTPUTS:
!       y_int:        Interpolate values corresponding to the X_INT abscissa
!                     values. Must have the same number of elements as the
!                     X_INT input argument.
!                     UNITS:      Data dependent, same as y.
!                     TYPE:       REAL(fp)
!                     DIMENSION:  Scalar or Rank-1 array
!                     ATTRIBUTES: INTENT(OUT)
!
! OPTIONAL OUTPUT ARGUMENTS:
!       RCS_Id:       Character string containing the Revision Control System
!                     Id field for the module.
!                     UNITS:      N/A
!                     TYPE:       CHARACTER(*)
!                     DIMENSION:  Scalar
!                     ATTRIBUTES: OPTIONAL, INTENT(OUT)
!
! FUNCTION RESULT:
!       Error_Status: The return value is an integer defining the error status.
!                     The error codes are defined in the Message_Handler module.
!                     If == SUCCESS the function executed normally.
!                        == FAILURE an unrecoverable error occurred.
!                     UNITS:      N/A
!                     TYPE:       INTEGER
!                     DIMENSION:  Scalar
!
! RESTRICTIONS:
!       Input X vectors *must* be monotonically ascending.
!
! PROCEDURE:
!       Taken from "Numerical Recipes in FORTRAN; The Art of Scientific Computing"
!       2nd ed. by W.H. Press et al., Cambridge University Press, 1992, New York,
!       New York, pp109.
!
!       Given arrays of a tabulated function, X and Y = F(X), and optionally
!       the array of the second derivatives of the interpolating function at
!       the points X, SPLINE_INTERPOLATE calculates a cubic-spline interpolated
!       value for the given value(s) of Xi.
!
! CREATION HISTORY:
!       Written by:   Paul van Delst, CIMSS/SSEC, 22-Mar-1996
!                     paul.vandelst@ssec.wisc.edu
!
!--------------------------------------------------------------------------------

  FUNCTION splint_scalar( x,            &  ! Input
                          y,            &  ! Input
                          x_int,        &  ! Input
                          y_int,        &  ! Output
                          y2,           &  ! Optional input
                          RCS_Id,       &  ! Optional output
                          Message_Log ) &  ! Error messaging
                        RESULT ( Error_Status )
    ! Arguments
    REAL(fp),               DIMENSION(:), INTENT(IN)  :: x
    REAL(fp),               DIMENSION(:), INTENT(IN)  :: y
    REAL(fp),                             INTENT(IN)  :: x_int
    REAL(fp),                             INTENT(OUT) :: y_int
    REAL(fp),     OPTIONAL, DIMENSION(:), INTENT(IN)  :: y2
    CHARACTER(*), OPTIONAL,               INTENT(OUT) :: RCS_Id
    CHARACTER(*), OPTIONAL,               INTENT(IN)  :: Message_Log
    ! Function result
    INTEGER :: Error_Status
    ! Local parameters
    CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'Spline_Interpolate(scalar)'
    ! Local variables
    CHARACTER(256) :: Message
    INTEGER :: n
    INTEGER :: i_lo, i_hi
    REAL(fp) :: a, b, dx
    REAL(fp), DIMENSION(SIZE(x)) :: y2_derivative

    ! ------
    ! Set up
    ! ------
    Error_Status = SUCCESS
    IF ( PRESENT( RCS_Id ) ) RCS_Id = MODULE_RCS_Id

    ! Are vector sizes consistent?
    n = SIZE(x)
    IF ( SIZE(y) /= n ) THEN
      Error_Status = FAILURE
      CALL Display_Message( ROUTINE_NAME, &
                            'Input X, Y arguments are different sizes.', &
                            Error_Status, &
                            Message_Log=Message_Log )
      RETURN
    END IF

    ! Check for second derivatives
    IF ( PRESENT( y2 ) ) THEN
      IF ( SIZE( y2 ) /= n ) THEN
        Error_Status = FAILURE
        CALL Display_Message( ROUTINE_NAME, &
                              'Input Y2 argument has inconsistent size.', &
                              Error_Status, &
                              Message_Log=Message_Log )
        RETURN
      END IF
      y2_derivative = y2
    ELSE
      ! Calculate them. This is inefficient so should only be done
      ! if X and Y are different each call.
      Error_Status = Spline_Initialize( x, y, y2_derivative,    &
                                        Message_Log=Message_Log )
      IF ( Error_Status /= SUCCESS ) THEN
        CALL Display_Message( ROUTINE_NAME, &
                              'Error occurred in spline initialisation function.', &
                              Error_Status, &
                              Message_Log=Message_Log )
        RETURN
      END IF
    END IF

    ! ---------------------------------------------------------
    ! Find the x-values in x that bracket the input x_int value
    ! by means of bisection. This is optimal if sequential
    ! calls to this function are at random values of x_int.
    ! ---------------------------------------------------------
    i_lo = MAX( 1, Bisection_Search( x, x_int ) )
    i_hi = MIN( i_lo + 1, n )

    ! ---------------------
    ! Calculate dx and test
    ! ---------------------
    dx = x(i_hi) - x(i_lo)
    IF ( dx < TOLERANCE ) THEN
      Error_Status = FAILURE
      WRITE( Message, '( "Error : dx < or = 0. i_lo, i_hi : ", i6, 2x, i6, &
                        &" x(i_lo), x(i_hi) : ", f12.6, 2x, f12.6 )' ) &
                      i_lo, i_hi, x(i_lo), x(i_hi)
      CALL Display_Message( ROUTINE_NAME, &
                            TRIM( Message ), &
                            Error_Status, &
                            Message_Log=Message_Log )
      RETURN
    END IF

    ! ---------------------
    ! Evaluate cubic dpline
    ! ---------------------
    a = ( x(i_hi) - x_int ) / dx
    b = ( x_int - x(i_lo) ) / dx
    y_int = ( a * y(i_lo) ) + ( b * y(i_hi) ) + &
            ( ( ( (a**3) - a ) * y2_derivative(i_lo) ) + &
              ( ( (b**3) - b ) * y2_derivative(i_hi) )   ) * &
            (dx*dx) / SIX

  END FUNCTION splint_scalar

  FUNCTION splint_rank1( x,            &  ! Input
                         y,            &  ! Input
                         x_int,        &  ! Input
                         y_int,        &  ! Output
                         y2,           &  ! Optional input
                         RCS_Id,       &  ! Optional output
                         Message_Log ) &  ! Error messaging
                       RESULT ( Error_Status )
    ! Arguments
    REAL(fp),               DIMENSION(:), INTENT(IN)  :: x      ! N
    REAL(fp),               DIMENSION(:), INTENT(IN)  :: y      ! N
    REAL(fp),               DIMENSION(:), INTENT(IN)  :: x_int  ! N_int
    REAL(fp),               DIMENSION(:), INTENT(OUT) :: y_int  ! N_int
    REAL(fp),     OPTIONAL, DIMENSION(:), INTENT(IN)  :: y2     ! N
    CHARACTER(*), OPTIONAL,               INTENT(OUT) :: RCS_Id
    CHARACTER(*), OPTIONAL,               INTENT(IN)  :: Message_Log
    ! Function result
    INTEGER :: Error_Status
    ! Local parameters
    CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'Spline_Interpolate(rank-1)'
    ! Local variables
    CHARACTER(256) :: Message
    INTEGER :: n, n_int
    INTEGER :: i
    REAL(fp)  :: a, b, dx
    INTEGER,  DIMENSION(SIZE(x_int)) :: i_lo, i_hi
    REAL(fp), DIMENSION(SIZE(x))     :: y2_derivative

    ! ------
    ! Set up
    ! ------
    Error_Status = SUCCESS
    IF ( PRESENT( RCS_Id ) ) RCS_Id = MODULE_RCS_Id

    ! Are input vector sizes consistent?
    n = SIZE(x)
    IF ( SIZE(y) /= n ) THEN
      Error_Status = FAILURE
      CALL Display_Message( ROUTINE_NAME, &
                            'Input X, Y arguments are different sizes.', &
                            Error_Status, &
                            Message_Log=Message_Log )
      RETURN
    END IF

    ! Are interpolate vector sizes consistent?
    n_int = SIZE(x_int)
    IF ( SIZE(y_int) /= n_int ) THEN
      Error_Status = FAILURE
      CALL Display_Message( ROUTINE_NAME, &
                            'Interpolate X_INT, Y_INT arguments are different sizes.', &
                            Error_Status, &
                            Message_Log=Message_Log )
      RETURN
    END IF

    ! Check for second derivatives
    IF ( PRESENT(y2) ) THEN
      IF ( SIZE(y2) /= n ) THEN
        Error_Status = FAILURE
        CALL Display_Message( ROUTINE_NAME, &
                              'Input Y2 argument has inconsistent size.', &
                              Error_Status, &
                              Message_Log=Message_Log )
        RETURN
      END IF
      y2_derivative = y2
    ELSE
      ! Calculate them.
      Error_Status = Spline_Initialize( x, y, y2_derivative,    &
                                        Message_Log=Message_Log )
      IF ( Error_Status /= SUCCESS ) THEN
        CALL Display_Message( ROUTINE_NAME, &
                              'Error occurred in spline initialisation function.', &
                              Error_Status, &
                              Message_Log=Message_Log )
        RETURN
      END IF
    END IF

    ! ---------------------------------------------------------
    ! Find the x-values in x that bracket the input x_int value
    ! by means of bisection. This is optimal if sequential
    ! calls to this function are at random values of x_int.
    ! ---------------------------------------------------------
    i_lo = Value_Locate( x, x_int )
    i_hi = i_lo + 1

    ! ----------------------------
    ! Calculate interpolated array
    ! ----------------------------
    DO i = 1, n_int

      ! Calculate dx and test
      dx = x(i_hi(i)) - x(i_lo(i))
      IF ( dx < TOLERANCE ) THEN
        Error_Status = FAILURE
        WRITE( Message, '( "Error : dx < or = 0. i_lo, i_hi : ", i6, 2x, i6, &
                          &" x(i_lo), x(i_hi) : ", f12.6, 2x, f12.6 )' ) &
                        i_lo(i), i_hi(i), x(i_lo(i)), x(i_hi(i))
        CALL Display_Message( ROUTINE_NAME, &
                              TRIM( Message ), &
                              Error_Status, &
                              Message_Log=Message_Log )
        RETURN
      END IF

      ! Evaluate cubic spline
      a = ( x(i_hi(i)) - x_int(i) ) / dx
      b = ( x_int(i) - x(i_lo(i)) ) / dx
      y_int(i) = ( a * y(i_lo(i)) ) + ( b * y(i_hi(i)) ) + &
                 ( ( ( (a**3) - a ) * y2_derivative(i_lo(i)) ) + &
                   ( ( (b**3) - b ) * y2_derivative(i_hi(i)) )   ) * &
                 (dx*dx) / SIX

    END DO

  END FUNCTION splint_rank1

END MODULE Spline_Interpolation

max_stars_repo_name: hsbadr/crtm

max_issues_repo_name: hsbadr/crtm

max_forks_repo_name: hsbadr/crtm


text:
c ******COMPDISS******************
      subroutine compDISS(lecho,ldiag,dE)

      implicit none
      include 'cdparams.fh'
      include 'cdwell0.fh'
      include 'cdwell1.fh'
      include 'cdwell2.fh'
      include 'cdisprop.fh'
      include 'cdcolls.fh'
      include 'cdrange.fh'
      include 'cdrates.fh'
      include 'cdcontrl.fh'
      include 'cdASATempPres.fh'

c     dmm 20000811
c     for printing out pops and well names we include labels for now
      include 'cdlabels.fh'
c     dmm 20000811

      integer mxSteps
      parameter (mxSteps = 1000)

c   local variables
      logical divert(mxWells)
      logical rkovar,inAcces(mxWells),noExit
      integer it,ip,iWell,iProd
      integer nEint,lecho,ldiag,nTdone,nElast(mxTpts),imostd,nvibs,
     2   ierr,ipfill
      real*8 E,fofEnum,Emin,rNoExNm,Edenom,Emax,dE,Estart,ACsave,
     2   Eoff,expEoff,rhoNum
      real*8 Qpart(mxTpts),stab(mxWells,mxTpts),ACpop(mxWells),
     2   ACarray(mxWells,mxWells),rateP(mxWells,0:mxProds),
     3   rateI(mxWells,mxWells),B(mxWells)
      real*8 rKofEI(mxWells,mxWells),rKofEP(mxWells,mxProds),
     2   ACDepI(mxWells),ACDepP(mxWells)

      real*8 rks,calRho,xNOS, calQ, tkinf, tksum, oldcalQ
      integer ipop
      real*8 tpop(mxWells),M
c tal 20000621
      equivalence (rkovar, rkold)
c tal 20000621

c******end declarations ***********

c    NOTE - this is just a modified version of compchem - however, our zero's
c    are different here E = 0 from bottom of Well

c       now energy must be defined with reference to the bottom of each
c       Well - however variable E in loop is referenced with respect to
c       offset energy Estart - hopefully this will enable us to avoid
c       underflows of the exponential terms  - renormalization will
c       occur at end <ayc 3/93>

c    setup resultant matrix for solver; init noDExit flag
      do 5 iWell = 1,nWells
         B(iWell) = 0.0d0
         noDExit(iWell) = .false.
 5    continue
      B(inpWell) = -1.d0

c     initialize density of states data block
      tksum = 0.0D0
c     dmm 20000510 commented out
c     call initRho
c     dmm 20000510

c     tal 20000621 adding switch
      if (oldrho) call initRho
c     tal 20000621

c     dmm 20000908 attempt to fix ACpop init bug ... zero out acpop
      call initCDcalc(ACpop, RateP)
c     dmm 20000908

c    find max barrier height in system, deepest Well depth
c    find minimum barrier to first exit (Product) channel
      imostd = 1
      Emax = 0.d0
      Emin = Exmax(inpWell)
      do 10 iWell = 1,nWells
         if (idpWell(iWell).gt.idpWell(imostd)) imostd = iWell
         Emax = dmax1(Exmax(iWell) + EWell(iWell),Emax)
         do 10 iProd = 1, NProds(iWell)
            Emin = dmin1(EbarLP(iWell) + EProd(iWell,iProd),Emin)
 10   continue

c     we start at minimum exit channel

c     dmm 20000510 commented out for new xNOS number-of-states
c     calculation.  No longer necessary.
c      Estart = Exmin(inpWell)
c dmm 20000510

c tal 20000621 adding switch
      if (oldrho) then
         Estart = Exmin(inpWell) + .5d0*dE
      else
         Estart = Exmin(inpWell)
      end if
c tal 20000621

c    IMPORTANT: for the low p limit, this subroutine must switch if
c    calculating a stabilization channel and there are no Product
c    channels at that E - this results in a different P-dependence
c    than the high E terms

c    we try adding an offset energy to avoid overflows in our exponentials
      Eoff = .5d0*(Emax - Estart)

      write(lecho,15) Estart,Emax,dE,dE/2.859d-3
 15   format(/,1x,'DISSOC calculations from E of ',
     2     f6.2,' to at least ',f6.2,' kcals,',/,
     3     ' in ',f5.2,' kcal (',f6.1,' cm-1) steps.')

c    initialize variables before E sum
      do 30 it = 1,nTemps
         nElast(it) = mxsteps
         do 30 iWell = 1,nWells
            do 30 iProd = 0,nProds(iWell)
               do 30 ip = -1,nPres+1
                  RK(iWell,iProd,it,ip) = 0.0d0
30    continue

c***********************************************************************
c dmm 20000510
c
c     Generate xrhoofn, or sum-of-states vector, for every well.  See
c     cdchem.f comments.  Skip if chemact calc has already been done
c     (assumes chemact calc will be done first if both are selected)

c      if(.not.chemact) then
c         call genRhos((eMax/2.859D-3),(dE/2.859D-3), temp(nTemps))
c      endif
c dmm 20000510

c tal 20000621
c      adding switch
      if ((.not.chemact).and.(.not.oldrho)) then
         call genRhos((eMax/2.859D-3),(dE/2.859D-3), temp(nTemps))
      endif
c tal 20000621

c***********************************************************************

c    calculate stabilization rate which is a function of T only
c    parameter rkovar will set up constant FE, original omega calculation
c    we use lowest E barrier for Well as Enot
c    we redirect output for stabilization parameters to ldiag

      do 35 iWell = 1,nWells
         write(ldiag,*) 'Calculating stabilization rate parameters for',
     2        ' Well #',iWell
         if (.not.rkovar) call rksset(ldiag,iWell,Exmin(iWell),
     2        nvibs,Edenom)
         do 35 it = 1,nTemps
            stab(iWell,it) = rks(ldiag,rkovar,temp(it),nvibs,
     2           Exmin(iWell),Edenom,iWell)
35    continue

c    compute partition function vs T
      do 37 it = 1,nTemps
c tal 20000621 for old vs. new chemdis
         if (oldrho) then
            Qpart(it) = oldcalQ(temp(it),inpWell)
         else
            Qpart(it) = calQ(temp(it),inpWell,dE)
         end if
c tal 20000621
 37   continue

      nEint = - 1
      nTdone = 0

c****start E loop************
c    loop is terminated when nTdone = nTemps
c    E is kcal units and is referenced from bottom of inpWell
c    note each step represents energy between E-dE/2 and E+dE/2
c    nEint is now the total # of iterations

c***********************************************************************
c     20000809 tal for printing out population vectors
c      open(23, file='ac.dat',status='unknown')
c     20000809 tal
c
c     20000809 dmm for printing out detailed balance terms
c      open(24, file='detail.txt', status = 'unknown')
c      write(24,27) temp(1), (ISname(ipop), ipop=1,nwells)
c 27   format (F8.2, 1x, 8(A10, 3x))
c***********************************************************************

c     tal 20000815 - for new rmin for dissoc
      rdiss = 0.d0
c     tal 20000815

 40   continue
      nEint = nEint + 1
      E = Estart + nEint*dE

c     get density of states for numerator
c     rhoNum = calRho(E,dE,inpWell,0)  dmm 20000510
c     rhoNum =  xNOS(E,dE,inpWell)
c     tal 20000621
      if (oldrho) then
         rhonum = calRho(E,dE,inpWell,0)
      else
         rhoNum =  xNOS(E,dE,inpWell)
      endif
c tal 20000621

c       compute all k(E)'s at this E
      call getKofEMatrix(E,dE,rKofEI,rKofEP)

c       we use EbarLP to flag inaccessible channels at this E
      do 330 iWell = 1,NWells
         if (E.le.EbarLP(iWell)) then
            inAcces(iWell) = .true.
         else
            inAcces(iWell) = .false.
         endif
 330  continue

c       start T loop
      do 500 it = 1,nTemps

c     variable nElast starts at mxsteps but is downsized
c     as determined directly below

         if (nEint.le.nElast(it)) then

c     calculate new f(E,T)
c     if using old rho method, do not offset by half-width for
c     exponential; if using new method, need offset in exponential
c     term
            if (oldrho) then
               fofEnum = rhoNum / Qpart(it)*
     $              dexp(-(E-Eoff)/(1.987d-3*temp(it)))
            else
c     calculate new f(E,T) with half-width ...
               fofEnum = rhoNum / Qpart(it)*
     2              dexp(-(E + 0.5*dE -Eoff)/(1.987d-3*temp(it)))
            endif
c***********************************************************************
c     dmm 20000623
c     TEMPORARY CALQ debugging lines ... these lines evaluate the HP
c     limit for dissociation and dump results to fort.91.  tksum will
c     be the integrated estimate of the HP limit for dissociation fed to
c     the code in fort.10.  It should not be higher than the limit!
c
c
c
c      if(.not.oldrho) then
c         tkinf = Aprod(1,1)*xNOS(E-Ewell(1)-Eprod(1,1),dE
c     $        ,inpWell)*dE*exp(-(E+0.5)/(1.987d-3*temp(it)))*
c     $        (temp(it)**rNprod(1,1))/Qpart(it)
c
c
c         tksum = tksum + tkinf
c         write(91, 111) temp(it), tkinf, tksum
c 111     format(F8.2, 1X, E16.9, 1X, E16.9)
c      endif
c     dmm 20000623
c***********************************************************************

c     setup T-dependent portion of RateI and RateP
            call fillAr1(it,rKofEI,RKofEP,RateI,RateP,
     2           ACDepI,ACDepP)

c     start P loop
            do 400 ip = 0,nPres+1

c     complete activated complex matrix AC and RateP
c     note M is factored in limits
c     also we do an incredible kludge in fillAr so the
c     low P dissoc limit works out for noExit conditions
               call fillAr2(it,ip,ACDepP,ACDepI,stab,ACarray,
     2              RateI,RateP,inAcces,noExit,rNoExNm)

c     this is for debugging matrix solver
c     if ((nEint.eq.0).or.(nEint.eq.1)) then
c     write(lecho,43) ip,it,ACarray(1,1),ACarray(1,2)
c     write(lecho,43) ip,it,ACarray(2,1),ACarray(2,2)
c     43                  format(1x,i2,' ',i2,' ',1pe14.7,' ',1pe14.7)
c     endif

c***********************************************************************
c     20000810 dmm code to print detailed balance terms
c               if((ip.gt.0).and.(ip.lt.nPres+1)) then
c                  do 101 ipop = 1, nwells
c                     tpop(ipop) = RateI(inpWell, ipop)*xNOS(E, dE
c     $                    ,inpWell)
c 101              continue
c                  write(24,102) E, (tpop(ipop), ipop=1, nWells)
c 102              format(F6.2, 1x, 8(E12.5, 1x))
c               endif
c     20000811 dmm
c***********************************************************************

c     now we are almost ready to solve matrix
c     hopefully, we've eliminated all singularities
               call solve(NWells,ACarray,B,idpWell,
     2               idpWell(imostd),ACpop,inAcces,ierr,lecho,E)

c***********************************************************************
c     20000809 tal code to print out population vectors (commented out by pey, 25/3/05)
c
c               if((ip.gt.0).and.(ip.lt.nPres+1)) then
c                  expEoff = dexp(-Eoff/(1.987d-3*temp(it)))
c                  do 111 ipop = 1, nwells
c
c                     tpop(ipop) = ACpop(ipop)*fofenum*expEoff
c     $                    *stab(inpWell,it)*pres(it,ip)/(82.1d0*temp(it)
c     $                    )
c 111              continue
c                  ipop =1
c                  write(23,112) E, (tpop(ipop),ipop=1,nwells)
c 112              format(F6.2, 1x, 8(E12.5, 1x))
c               endif
c***********************************************************************
c     20000809 tal code end

c     here's our kludge!!! recognizing the solution
c     vector goes as one over our very small M factor, we
c     multiply it out by same factor; also trip noDExit
c     flag and divert flag for latter

               if ((ip.eq.0).and.(noExit)) then
                  do 340 iWell = 1,NWells
                     ACpop(iWell) = ACpop(iWell) * rNoExNm
                     if (.not.inAcces(iWell)) then
                        noDExit(iWell) = .true.
                        divert(iWell) = .true.
                     else
                        divert(iWell) = .false.
                     endif
 340              continue

               else
                  do 345 iWell = 1,NWells
                     divert(iWell) = .false.
 345              continue
               endif

c                we save parameter for convergence criterion
               if (ip.eq.nPres+1) ACsave = ACpop(iwcmax)

c     more matrix debugging stuff
c     if ((nEint.eq.0).or.(nEint.eq.1))
c     2               write(lecho,*) ACpop(1),ACpop(2)

               if (ierr.eq.1) then
                  write(lecho,*) ' '
                  write(lecho,*) 'Failure at nEint, E, it, ip: ',
     2                 nEint,E,it,ip
                  stop
               endif

c     finally compute contribution to RK and keep
c     running sum over all levels
               do 350 iWell = 1,nWells
                  do 350 iProd = 0,nProds(iWell)
                     ipfill = ip
                     if ((ip.eq.0).and.divert(iWell).and.
     2                    (iProd.eq.0)) ipfill = -1
                     RK(iWell,iProd,it,ipfill) = RK(iWell,iProd,
     2                     it,ipfill) + fofEnum*RateP(iWell,iProd)
     3                    *ACpop(iWell)
 350           continue
 400        continue
c  end P loop

c     tal 20000815
c     setting up new criteria for dissoc ASA cases
            if (ASA) then
               M = pres(icurtemp(1),icurpres(1)) / (82.1 *
     2              temp(icurtemp(1)))
               rdiss = rdiss + fofenum*stab(inpwell,it)*M*
     2			dexp(-Eoff/(1.987d-3*temp(it)))
            end if
c     tal 20000815

c             check to see if this is our last sum vs T
c             don't need to sum any more terms if E is greater
c             and E is greater then Emax and there is convergence on
c             highest Product channel - we believe high P limit is
c             slowest to converge

            if (((E-1.d0*dE.gt.Emax).and.(RK(iwcmax,ipcmax,it,
     2           nPres+1).gt.1.d-175)
     3           .and.(fofEnum*RateP(iwcmax,ipcmax)*ACsave/
     4           RK(iwcmax,ipcmax,it,nPres+1).lt.2.d-3)).or.
     5           (nEint.ge.mxsteps)) then
               nElast(it) = nEint
               nTdone = nTdone + 1
            endif

c     when N = nElast(T) there will be no more iterations
c     on this T; when nTdone = nTemps we are done

         endif
 500  continue
c end T loop
      if (nTdone.lt.nTemps) goto 40
c*****E loop done
         close(23)
c    finally we normalize RK by stabilization rate and account for offset
c    we also factor M from the high and low p limits
      do 700 it = 1,nTemps
         expEoff = dexp(-Eoff/(1.987d-3*temp(it)))

         do 700 iWell = 1,nWells
            do 700 iProd = 0,nProds(iWell)
              do 700 ip = -1,nPres+1
                  if ((ip.le.0).or.(ip.eq.nPres+1)) then
                     RK(iWell,iProd,it,ip) =  RK(iWell,iProd,it,ip)
     2                  *expEoff*stab(inpWell,it)
                  else
                     RK(iWell,iProd,it,ip) =
     2                  RK(iWell,iProd,it,ip)*expEoff*stab(inpWell,it)
     3                  * pres(it,ip) / (82.1d0*temp(it))

c***********************************************************************
c     dmm 20000823 print out ks*M (temporary)

                     if ((iWell.eq.1).and.(iProd.eq.1)) then
                        write(*,701) temp(it), pres(it,ip), stab(inpWell
     $                       ,it)*pres(it,ip)/(82.1d0*temp(it))
 701                    format('T = ', F8.2, 1x,'P = ', F8.2, 1x,
     $                       'ks[M] = ', E13.3)

                     endif
c***********************************************************************

                  endif
                  if (RK(iWell,iProd,it,ip).lt.1.d-200)
     2               RK(iWell,iProd,it,ip) = 1.d-200
700   continue

c    echo out some integration parameters
      write(lecho,710)
710   format(/,1x,'INTEGRATION STEP RANGE',
     2   /,1x,' T (K) ',2x,'# steps',2x,'Estart (kcal) Eend (kcal)')

      do 760 it = 1,nTemps
         write(lecho,750) temp(it),nElast(it)+1,Estart,
     2      Estart+dE*nElast(it)
750      format(1x,f7.1,4x,i3,4x,2x,f7.2,4x,f7.2)
760   continue

      return
      end

c $Author$
c $Date$
c $Source$
c $Revision$
c $Log$
c Revision 1.1  2007-02-20 23:10:23  sandeeps
c Initial revision
c
c Revision 1.11  2000/08/18 13:11:06  dmmatheu
c after merge w/tlada code ... changes include pop vector labelling and
c printing, detailed balance labeling and printing
c
c Revision 1.10  2000/08/09 18:20:48  dmmatheu
c after tlada changes to print out population vectors
c
c Revision 1.9  2000/06/26 14:00:50  dmmatheu
c successful merge with tlada code
c
c Revision 1.8  2000/06/23 17:08:14  dmmatheu
c after merge with Tom Lada code
c
c Revision 1.7  2000/06/23 17:00:08  dmmatheu
c before merge with Tom Lada code
c
c Revision 1.6  2000/06/22 21:59:43  dmmatheu
c dissoc side fixed -- Q consistent with fofEnum and fofEnum fixed ...
c OK!
c
c Revision 1.5  2000/06/22 17:10:46  dmmatheu
c Note:  Before fixing "greater rate constant at high pressure than HP
c limit rate constant" in dissoc problem.  Basically this comes from the
c fact that fofEnum is evaluated with a dE that's usually too big (1
c kcal) so that either the exponential term must either be evaluated halfway
c through (at E + 0.5*dE) or the partition function must be evaluated
c "erroneously" (e.g. with grainsize dE = 1 kcal), or even both.  It's
c not clear what's best.
c
c Revision 1.4  2000/05/11 15:36:47  dmmatheu
c Doesn't work well ... final results ~2 OOM higher than they should
c be.  Likely due to use of calQ and Qpart which may be inconsistent
c with our xNOS.  There is some serious unit disagreement ugliness here
c -- I can't seem to get to the bottom of it!  Someday may need to iron
c out calRho/xNOS differences once and for all.  For now try hack
c solution (running denom approach as in cdchem.f)
c
c Revision 1.3  2000/05/10 19:20:43  dmmatheu
c before altering calRhos to put in xNOS calls ...
c

max_issues_repo_name: bslakman/RMG-Java

max_forks_repo_name: bslakman/RMG-Java


text:
      SUBROUTINE BMG2_SymStd_restrict( &
     &                       Q, QC, CI, &
     &                       Nx, Ny, Nxc, Nyc, JPN&
     &                       ) BIND(C, NAME='BMG2_SymStd_restrict')

! ======================================================================
!  --------------------
!   DESCRIPTION:
!  --------------------
!
!     BMG2_SymStd_restrict computes the restriction of a vector on the
!     fine grid, Q, to a vector on the coarse grid, QC.  The weights
!     involve the transpose of the interpolation operator from the
!     coarse grid to the fine grid.
!
! ======================================================================
! $license_flag$
! ======================================================================
!  --------------------
!   INPUT:
!  --------------------
!
!
! ======================================================================
!  --------------------
!   INPUT/OUTPUT:
!  --------------------
!
!
! ======================================================================
!  --------------------
!   OUTPUT:
!  --------------------
!
!
! ======================================================================
!  --------------------
!   LOCAL:
!  --------------------
!
!
! ======================================================================
      USE ModInterface
      IMPLICIT NONE

! -----------------------------
!     Includes
!
      INCLUDE 'BMG_constants_f90.h'
      INCLUDE 'BMG_stencils_f90.h'
      INCLUDE 'BMG_parameters_f90.h'

! ----------------------------
!     Argument Declarations
!
      INTEGER(len_t), VALUE ::  Nx, Nxc, Ny, Nyc

      INTEGER(C_INT), VALUE :: JPN
      REAL(real_t) :: CI(Nxc,Nyc,8), Q(Nx,Ny), QC(Nxc,Nyc)

! ----------------------------
!     Local Declarations
!
      INTEGER  i, ic, IPN, j, jc, PER_x, PER_y, PER_xy

! ======================================================================

! ---------------------------------------------
!     Restrict the vector Q -> QC:
! ---------------------------------------------
      IPN = IABS(JPN)
      IF(IPN.EQ.BMG_BCs_definite .OR. IPN.EQ.BMG_BCs_indef_nonper)&
     &     THEN

         j=0
         DO jc=2, Nyc-1
            j=j+2
            i=0
            DO ic=2, Nxc-1
               i=i+2
               QC(ic,jc) = CI(ic,jc,LNE)*Q(i-1,j-1)&
     &              + CI(ic,jc,LA)*Q(i,j-1)&
     &              + CI(ic+1,jc,LNW)*Q(i+1,j-1)&
     &              + CI(ic,jc,LR)*Q(i-1,j)&
     &              + Q(i,j)&
     &              + CI(ic+1,jc,LL)*Q(i+1,j)&
     &              + CI(ic,jc+1,LSE)*Q(i-1,j+1)&
     &              + CI(ic,jc+1,LB)*Q(i,j+1)&
     &              + CI(ic+1,jc+1,LSW)*Q(i+1,j+1)
            ENDDO
         ENDDO

      ELSE

         PER_x = IABS(BMG_BCs_def_per_x)
         PER_y = IABS(BMG_BCs_def_per_y)
         PER_xy = IABS(BMG_BCs_def_per_xy)

         IF((IPN.EQ.PER_y.OR.IPN.EQ.PER_xy).AND.Ny/2+1.EQ.Nyc) THEN
            DO I=1,Nx
               Q(I,1) = Q(I,Ny-1)
               Q(I,Ny) = Q(I,2)
            ENDDO
         ENDIF
         IF((IPN.EQ.PER_x.OR.IPN.EQ.PER_xy).AND.Nx/2+1.EQ.Nxc) THEN
            DO J=1,Ny
               Q(1,J) = Q(Nx-1,J)
               Q(Nx,J) = Q(2,J)
            ENDDO
         ENDIF
         j=0
         DO jc=2, Nyc-1
            j=j+2
            i=0
            DO ic=2, Nxc-1
               i=i+2
               QC(ic,jc) = CI(ic,jc,LNE)*Q(i-1,j-1)&
     &              + CI(ic,jc,LA)*Q(i,j-1)&
     &              + CI(ic+1,jc,LNW)*Q(i+1,j-1)&
     &              + CI(ic,jc,LR)*Q(i-1,j)&
     &              + Q(i,j)&
     &              + CI(ic+1,jc,LL)*Q(i+1,j)&
     &              + CI(ic,jc+1,LSE)*Q(i-1,j+1)&
     &              + CI(ic,jc+1,LB)*Q(i,j+1)&
     &              + CI(ic+1,jc+1,LSW)*Q(i+1,j+1)
            ENDDO
         ENDDO

      ENDIF

! ======================================================================

       RETURN
       END

max_stars_repo_name: cedar-framework/cedar

max_issues_repo_name: cedar-framework/cedar

max_forks_repo_name: cedar-framework/cedar