text:
! Initialize a divergence-free velocity field, add the gradient of a
! scalar, project, and check if we get the original field back.

subroutine varden()

  use variables
  use fabio_module
  use multifab_module, only: multifab, build, multifab_build_edge, destroy
  use ml_layout_module, only: ml_layout
  use ml_boxarray_module, only: ml_boxarray
  use runtime_init_module, only: runtime_init, runtime_close
  use probin_module, only: test_set, pmask, max_levs, &
       prob_lo, prob_hi, drdxfac, nodal, run_prefix, project_type
  use geometry, only: spherical, init_spherical, &
       init_center, destroy_geometry, &
       dr_fine, nr_fine, &
       init_cutoff, init_radial, init_multilevel, initialize_dx
  use box_util_module, only: read_a_hgproj_grid
  use bl_IO_module
  use bl_constants_module
  use define_bc_module, only: bc_tower, bc_tower_level_build, bc_tower_destroy, initialize_bc
  use test_projection_module
  use proj_parameters
  use hgproject_module
  use macproject_module

  implicit none

  integer :: i, n, nlevs, dm, comp
  integer :: ng_s

  type(ml_boxarray) :: mba
  type(ml_layout)   :: mla
  type(bc_tower)    :: the_bc_tower

  type(multifab), allocatable :: uold(:)
  type(multifab), allocatable :: umid(:)
  type(multifab), allocatable :: unew(:)
  type(multifab), allocatable :: gphi(:)

  type(multifab), allocatable :: umac_old(:,:)
  type(multifab), allocatable :: umac_mid(:,:)
  type(multifab), allocatable :: umac_new(:,:)
  type(multifab), allocatable :: gphi_mac(:,:)
  type(multifab), allocatable :: utemp(:)

  type(multifab), allocatable :: rhohalf(:)
  type(multifab), allocatable :: beta0(:)
  type(multifab), allocatable :: pi(:)
  type(multifab), allocatable :: macpi(:)
  type(multifab), allocatable :: gpi(:)

  real(kind=dp_t), pointer :: dx(:,:)

  real(kind=dp_t) :: lenx,leny,lenz,max_dist

  real(kind=dp_t) :: time, dt

  real(dp_t), parameter :: SMALL = 1.d-13

  character (len=20), allocatable :: plot_names(:)

  !---------------------------------------------------------------------------
  ! initialize some things
  !---------------------------------------------------------------------------

  ! this reads in the inputs file and sets the values of the parameters
  ! in the probin_module
  call runtime_init()

  ! set spherical in the geometry module
  call init_spherical()

  ! set the center() array in geometry with the coordinates
  ! of the center of the domain
  call init_center()

  ! initialize the integer keys used for variable names
  call init_variables()

  time = ZERO

  !---------------------------------------------------------------------------
  ! setup the grid
  !---------------------------------------------------------------------------

  ! we only use fixed grids
  if (test_set /= '') then
     call read_a_hgproj_grid(mba, test_set)
     call ml_layout_build(mla,mba,pmask)
  else
     call bl_error('test_projection only implemented with fixed grids!')
  endif

  ! check for proper nesting
  if (.not. ml_boxarray_properly_nested(mla%mba, 3, pmask)) then
     call bl_error('fixed_grids not properly nested')
  end if

  dm = mla%dim
  nlevs = mla%nlevel

  ! initialize dx
  call initialize_dx(dx,mba,nlevs)

  ! initialize boundary conditions
  call initialize_bc(the_bc_tower,nlevs,pmask)
  do n = 1,nlevs
     call bc_tower_level_build(the_bc_tower,n,mla%la(n))
  end do

  !---------------------------------------------------------------------------
  ! sanity checks
  !---------------------------------------------------------------------------
  if (nlevs .ne. max_levs) then
     call bl_error('varden.f90: nlevs .ne. max_levs not supported yet')
  end if

  if (dm .ne. get_dim(mla%mba)) then
     call bl_error('dm_in not properly set in inputs file')
  endif

  if (spherical .eq. 1 ) then
     call bl_error('spherical = 1 not supported')
  endif

  if (abs(dx(1,1) - dx(1,2)) > SMALL) then
     call bl_error('zones must be square')
  endif

  !---------------------------------------------------------------------------
  ! allocate arrays
  !---------------------------------------------------------------------------
  ng_s = 4

  if (project_type == 1) then

     ! HG projection.  Velocities are cell-centered

     allocate(uold(nlevs), umid(nlevs), unew(nlevs), gphi(nlevs))

     do n = 1, nlevs
        call build(uold(n),  mla%la(n), dm, ng_s)
        call setval(uold(n), ZERO, all=.true.)

        call build(umid(n),  mla%la(n), dm, ng_s)
        call setval(umid(n), ZERO, all=.true.)

        call build(unew(n),  mla%la(n), dm, ng_s)
        call setval(unew(n), ZERO, all=.true.)

        call build(gphi(n),  mla%la(n), dm, ng_s)
        call setval(gphi(n), ZERO, all=.true.)
     enddo

  else if (project_type == 2) then

     ! MAC projection.  Velocities are nodal in respective dimension

     allocate(umac_old(nlevs,dm), umac_mid(nlevs,dm), umac_new(nlevs,dm), &
              gphi_mac(nlevs,dm))

     do n = 1, nlevs
        do comp = 1, dm
           call multifab_build_edge(umac_old(n,comp),  mla%la(n), 1, ng_s, comp)
           call setval(umac_old(n,comp), ZERO, all=.true.)

           call multifab_build_edge(umac_mid(n,comp),  mla%la(n), 1, ng_s, comp)
           call setval(umac_mid(n,comp), ZERO, all=.true.)

           call multifab_build_edge(umac_new(n,comp),  mla%la(n), 1, ng_s, comp)
           call setval(umac_new(n,comp), ZERO, all=.true.)

           call multifab_build_edge(gphi_mac(n,comp),  mla%la(n), 1, ng_s, comp)
           call setval(gphi_mac(n,comp), ZERO, all=.true.)
        enddo
     enddo

     ! some initialization stuff wants a cell-centered multifab, also,
     ! we will need a container to map the MAC velocity to a
     ! cell-centered velocity for output.
     allocate(utemp(nlevs))

     do n = 1, nlevs
        call build(utemp(n),  mla%la(n), dm, ng_s)
        call setval(utemp(n), ZERO, all=.true.)
     enddo

  else

     call bl_error("ERROR: invalid project_type")

  endif

  ! some stats
  if (parallel_IOProcessor()) then
     print *, 'number of processors = ', parallel_nprocs()
     print *, 'number of dimensions = ', dm
     do n = 1, nlevs
        print *, 'level: ', n
        print *, '   number of boxes = ', mla%la(n)%lap%nboxes
        print *, '   maximum zones   = ', (extent(mla%mba%pd(n),i),i=1,dm)
     end do
     print *, ''
  end if

  !------------------------------------------------------------------------
  ! some multilevel base state initialization
  !------------------------------------------------------------------------

  ! setup some base state coordinate information: nr_fine and dr_fine in
  ! geometry
  if (spherical .eq. 1) then

     ! for spherical, we will now require that dr_fine = dx
     dr_fine = dx(nlevs,1) / dble(drdxfac)

     lenx = HALF * (prob_hi(1) - prob_lo(1))
     leny = HALF * (prob_hi(2) - prob_lo(2))
     lenz = HALF * (prob_hi(3) - prob_lo(3))

     max_dist = sqrt(lenx**2 + leny**2 + lenz**2)
     nr_fine = int(max_dist / dr_fine) + 1

  else

     nr_fine = extent(mla%mba%pd(nlevs),dm)
     dr_fine = (prob_hi(dm)-prob_lo(dm)) / dble(nr_fine)

  end if

  ! create numdisjointchunks, r_start_coord, r_end_coord
  if (project_type == 1) then
     call init_multilevel(uold)
  else
     call init_multilevel(utemp)
  endif

  ! now that we have nr_fine and dr_fine we can create nr, dr,
  ! r_cc_loc, r_edge_loc
  call init_radial(nlevs,mba)

  ! allocate the cutoff coordinate arrays
  call init_cutoff(nlevs)

  !---------------------------------------------------------------------------
  ! initialize velocity field
  !---------------------------------------------------------------------------
  if (project_type == 1) then
     call init_velocity(uold, dx, mla, the_bc_tower%bc_tower_array)
  else
     call init_mac_velocity(umac_old, dx, mla, the_bc_tower%bc_tower_array)
  endif

  ! output the initial velocity field
  allocate(plot_names(dm))
  plot_names(1) = "x-velocity"
  if (dm >= 2) plot_names(2) = "y-velocity"
  if (dm == 3) plot_names(3) = "z-velocity"

  if (project_type == 1) then

     call fabio_ml_write(uold, mla%mba%rr(:,1), trim(run_prefix) // "u_init", &
                         names=plot_names)

     ! copy the velocity field over to the intermediate state, umid
     do n = 1, nlevs
        call multifab_copy(umid(n), uold(n), nghost(uold(n)))
     enddo

  else

     ! cannot write out a MAC field -- convert to cell-centered
     call convert_MAC_to_cc(umac_old, utemp)

     call fabio_ml_write(utemp, mla%mba%rr(:,1), trim(run_prefix) // "u_init", &
                         names=plot_names)

     ! copy the velocity field over to the intermediate state, umid
     do n = 1, nlevs
        do comp = 1, dm
           call multifab_copy_c(umac_mid(n,comp), 1, umac_old(n,comp), 1, 1, &
                                nghost(umac_old(n,comp)))
        enddo
     enddo

  end if

  !---------------------------------------------------------------------------
  ! 'pollute' the velocity field by adding the gradient of a scalar
  !---------------------------------------------------------------------------

  if (project_type == 1) then
     call add_grad_scalar(umid, gphi, dx, mla, the_bc_tower%bc_tower_array)

     call fabio_ml_write(umid, mla%mba%rr(:,1), trim(run_prefix) // "u_plus_grad_phi", &
                         names=plot_names)

     call fabio_ml_write(gphi, mla%mba%rr(:,1), trim(run_prefix) // "grad_phi", &
                         names=plot_names)

     ! copy the velocity field over to the final state, unew
     do n = 1, nlevs
        call multifab_copy(unew(n), umid(n), nghost(uold(n)))
     enddo

  else
     call add_grad_scalar_mac(umac_mid, gphi_mac, dx, mla, the_bc_tower%bc_tower_array)

     ! cannot write out a MAC field -- convert to cell-centered
     call convert_MAC_to_cc(umac_mid, utemp)

     call fabio_ml_write(utemp, mla%mba%rr(:,1), trim(run_prefix) // "u_plus_grad_phi", &
                         names=plot_names)

     call convert_MAC_to_cc(gphi_mac, utemp)

     call fabio_ml_write(utemp, mla%mba%rr(:,1), trim(run_prefix) // "grad_phi", &
                         names=plot_names)

     ! copy the velocity field over to the final state, unew
     do n = 1, nlevs
        do comp = 1, dm
           call multifab_copy_c(umac_new(n,comp), 1, umac_mid(n,comp), 1, 1, &
                                nghost(umac_mid(n,comp)))
        enddo
     enddo
  endif

  !---------------------------------------------------------------------------
  ! project out the divergent portion of the velocity field
  !---------------------------------------------------------------------------
  if (project_type == 1) then

     ! hgprojection -- here pi is nodal and u is cell-centered

     allocate(rhohalf(nlevs), pi(nlevs), gpi(nlevs), beta0(nlevs))

     do n = 1, nlevs

        ! build the density used in the projection -- we are just doing
        ! constant density, so set it to 1
        call build(rhohalf(n),  mla%la(n), 1, ng_s)
        call setval(rhohalf(n), ONE, all=.true.)

        ! build pi
        call build(pi(n),  mla%la(n), 1, 1, nodal)
        call setval(pi(n), ZERO, all=.true.)

        ! build gpi
        call build(gpi(n),  mla%la(n), dm, 1)
        call setval(gpi(n), ZERO, all=.true.)

        ! build the coefficient in the divergence.  We are doing
        ! divergence-free (incompressible), so set beta0 = 1
        call build(beta0(n),  mla%la(n), 1, 1)
        call setval(beta0(n), ONE, all=.true.)

     enddo

     ! hgproject takes dt -- it has no meaning for our projection type
     dt = ONE

     call hgproject(initial_projection_comp, &
                    mla, &
                    unew, unew, rhohalf, &
                    pi, gpi, &
                    dx, dt, the_bc_tower, &
                    beta0)

     call fabio_ml_write(unew, mla%mba%rr(:,1), trim(run_prefix) // "u_new", &
                         names=plot_names)

  else

     ! mac projection -- here pi is cell-centered and u is MAC
     allocate(rhohalf(nlevs), macpi(nlevs))

     do n = 1, nlevs

        ! build the density used in the projection -- we are just doing
        ! constant density, so set it to 1
        call build(rhohalf(n),  mla%la(n), 1, ng_s)
        call setval(rhohalf(n), ONE, all=.true.)

        ! build macpi -- macproject assumes it has 1 ghostcell
        call build(macpi(n),  mla%la(n), 1, 1)
        call setval(macpi(n), ZERO, all=.true.)

     enddo

     call macproject(mla, &
                     umac_new, macpi, rhohalf, &
                     dx, the_bc_tower)

     ! convert to cell-centered for output
     call convert_MAC_to_cc(umac_new, utemp)

     call fabio_ml_write(utemp, mla%mba%rr(:,1), trim(run_prefix) // "u_new", &
                         names=plot_names)

  endif

  !---------------------------------------------------------------------------
  ! clean-up
  !---------------------------------------------------------------------------

  if (project_type == 1) then
     do n = 1, nlevs
        call destroy(uold(n))
        call destroy(umid(n))
        call destroy(unew(n))
        call destroy(gphi(n))

        call destroy(rhohalf(n))
        call destroy(pi(n))
        call destroy(gpi(n))
        call destroy(beta0(n))
     enddo

     deallocate(uold, umid, unew, gphi)
     deallocate(rhohalf, pi, gpi, beta0)

  else
     do n = 1, nlevs
        do comp = 1, dm
           call destroy(umac_old(n,comp))
           call destroy(umac_mid(n,comp))
           call destroy(umac_new(n,comp))
           call destroy(gphi_mac(n,comp))
        enddo

        call destroy(utemp(n))

        call destroy(rhohalf(n))
        call destroy(macpi(n))

     enddo

     deallocate(umac_old, umac_mid, umac_new, gphi_mac, utemp)
     deallocate(rhohalf, macpi)
  endif

  call destroy(mla)
  call destroy(mba)

  deallocate(plot_names)

  call bc_tower_destroy(the_bc_tower)

  call runtime_close()

  call destroy_geometry()

end subroutine varden

max_stars_repo_name: sailoridy/MAESTRO

max_issues_repo_name: sailoridy/MAESTRO

max_forks_repo_name: sailoridy/MAESTRO


text:
-makelib xcelium_lib/xilinx_vip -sv \
  "C:/Xilinx/Vivado/2020.2/data/xilinx_vip/hdl/axi4stream_vip_axi4streampc.sv" \
  "C:/Xilinx/Vivado/2020.2/data/xilinx_vip/hdl/axi_vip_axi4pc.sv" \
  "C:/Xilinx/Vivado/2020.2/data/xilinx_vip/hdl/xil_common_vip_pkg.sv" \
  "C:/Xilinx/Vivado/2020.2/data/xilinx_vip/hdl/axi4stream_vip_pkg.sv" \
  "C:/Xilinx/Vivado/2020.2/data/xilinx_vip/hdl/axi_vip_pkg.sv" \
  "C:/Xilinx/Vivado/2020.2/data/xilinx_vip/hdl/axi4stream_vip_if.sv" \
  "C:/Xilinx/Vivado/2020.2/data/xilinx_vip/hdl/axi_vip_if.sv" \
  "C:/Xilinx/Vivado/2020.2/data/xilinx_vip/hdl/clk_vip_if.sv" \
  "C:/Xilinx/Vivado/2020.2/data/xilinx_vip/hdl/rst_vip_if.sv" \
-endlib
-makelib xcelium_lib/xpm -sv \
  "C:/Xilinx/Vivado/2020.2/data/ip/xpm/xpm_cdc/hdl/xpm_cdc.sv" \
  "C:/Xilinx/Vivado/2020.2/data/ip/xpm/xpm_memory/hdl/xpm_memory.sv" \
-endlib
-makelib xcelium_lib/xpm \
  "C:/Xilinx/Vivado/2020.2/data/ip/xpm/xpm_VCOMP.vhd" \
-endlib
-makelib xcelium_lib/axi_infrastructure_v1_1_0 \
  "../../../../lab3.gen/sources_1/bd/lab3/ipshared/ec67/hdl/axi_infrastructure_v1_1_vl_rfs.v" \
-endlib
-makelib xcelium_lib/axi_vip_v1_1_8 -sv \
  "../../../../lab3.gen/sources_1/bd/lab3/ipshared/94c3/hdl/axi_vip_v1_1_vl_rfs.sv" \
-endlib
-makelib xcelium_lib/processing_system7_vip_v1_0_10 -sv \
  "../../../../lab3.gen/sources_1/bd/lab3/ipshared/34f8/hdl/processing_system7_vip_v1_0_vl_rfs.sv" \
-endlib
-makelib xcelium_lib/xil_defaultlib \
  "../../../../lab3.gen/sources_1/bd/lab3/ip/lab3_processing_system7_0_0/sim/lab3_processing_system7_0_0.v" \
-endlib
-makelib xcelium_lib/axi_lite_ipif_v3_0_4 \
  "../../../../lab3.gen/sources_1/bd/lab3/ipshared/66ea/hdl/axi_lite_ipif_v3_0_vh_rfs.vhd" \
-endlib
-makelib xcelium_lib/lib_cdc_v1_0_2 \
  "../../../../lab3.gen/sources_1/bd/lab3/ipshared/ef1e/hdl/lib_cdc_v1_0_rfs.vhd" \
-endlib
-makelib xcelium_lib/interrupt_control_v3_1_4 \
  "../../../../lab3.gen/sources_1/bd/lab3/ipshared/a040/hdl/interrupt_control_v3_1_vh_rfs.vhd" \
-endlib
-makelib xcelium_lib/axi_gpio_v2_0_24 \
  "../../../../lab3.gen/sources_1/bd/lab3/ipshared/4318/hdl/axi_gpio_v2_0_vh_rfs.vhd" \
-endlib
-makelib xcelium_lib/xil_defaultlib \
  "../../../../lab3.gen/sources_1/bd/lab3/ip/lab3_axi_gpio_0_0/sim/lab3_axi_gpio_0_0.vhd" \
-endlib
-makelib xcelium_lib/generic_baseblocks_v2_1_0 \
  "../../../../lab3.gen/sources_1/bd/lab3/ipshared/b752/hdl/generic_baseblocks_v2_1_vl_rfs.v" \
-endlib
-makelib xcelium_lib/axi_register_slice_v2_1_22 \
  "../../../../lab3.gen/sources_1/bd/lab3/ipshared/af2c/hdl/axi_register_slice_v2_1_vl_rfs.v" \
-endlib
-makelib xcelium_lib/fifo_generator_v13_2_5 \
  "../../../../lab3.gen/sources_1/bd/lab3/ipshared/276e/simulation/fifo_generator_vlog_beh.v" \
-endlib
-makelib xcelium_lib/fifo_generator_v13_2_5 \
  "../../../../lab3.gen/sources_1/bd/lab3/ipshared/276e/hdl/fifo_generator_v13_2_rfs.vhd" \
-endlib
-makelib xcelium_lib/fifo_generator_v13_2_5 \
  "../../../../lab3.gen/sources_1/bd/lab3/ipshared/276e/hdl/fifo_generator_v13_2_rfs.v" \
-endlib
-makelib xcelium_lib/axi_data_fifo_v2_1_21 \
  "../../../../lab3.gen/sources_1/bd/lab3/ipshared/54c0/hdl/axi_data_fifo_v2_1_vl_rfs.v" \
-endlib
-makelib xcelium_lib/axi_crossbar_v2_1_23 \
  "../../../../lab3.gen/sources_1/bd/lab3/ipshared/bc0a/hdl/axi_crossbar_v2_1_vl_rfs.v" \
-endlib
-makelib xcelium_lib/xil_defaultlib \
  "../../../../lab3.gen/sources_1/bd/lab3/ip/lab3_xbar_0/sim/lab3_xbar_0.v" \
-endlib
-makelib xcelium_lib/proc_sys_reset_v5_0_13 \
  "../../../../lab3.gen/sources_1/bd/lab3/ipshared/8842/hdl/proc_sys_reset_v5_0_vh_rfs.vhd" \
-endlib
-makelib xcelium_lib/xil_defaultlib \
  "../../../../lab3.gen/sources_1/bd/lab3/ip/lab3_rst_ps7_0_111M_0/sim/lab3_rst_ps7_0_111M_0.vhd" \
  "../../../../lab3.gen/sources_1/bd/lab3/VGA_controller_v2.vhd" \
  "../../../../lab3.gen/sources_1/bd/lab3/ipshared/c5a7/hdl/VGA_controller_v2_v1_0_S00_AXI.vhd" \
  "../../../../lab3.gen/sources_1/bd/lab3/ipshared/c5a7/hdl/VGA_controller_v2_v1_0.vhd" \
  "../../../../lab3.gen/sources_1/bd/lab3/ip/lab3_VGA_controller_v2_0_0/sim/lab3_VGA_controller_v2_0_0.vhd" \
-endlib
-makelib xcelium_lib/axi_protocol_converter_v2_1_22 \
  "../../../../lab3.gen/sources_1/bd/lab3/ipshared/5cee/hdl/axi_protocol_converter_v2_1_vl_rfs.v" \
-endlib
-makelib xcelium_lib/xil_defaultlib \
  "../../../../lab3.gen/sources_1/bd/lab3/ip/lab3_auto_pc_0/sim/lab3_auto_pc_0.v" \
-endlib
-makelib xcelium_lib/xil_defaultlib \
  "../../../../lab3.gen/sources_1/bd/lab3/sim/lab3.vhd" \
-endlib
-makelib xcelium_lib/xil_defaultlib \
  glbl.v
-endlib

ext: f


!------------------------------------------------------------------------------
! HLRS - Numerical Methods and Libraries, Nobelstraße 19, 70569 Stumgart
!------------------------------------------------------------------------------
!
! MODULE: strings
!
!> @author George Benthien
!> @author Johannes Gebert, gebert@hlrs.de, HLRS/NUM
!
!> @Descrition:
!> Several functions to get some string work done. Mainly developed by
!> George Benthien
!> http://numerical.recipes/forum/showthread.php?t=1659
!> http://www.gbenthien.net/strings/index.html
!> slight modifications by Johannes Gebert on 06.03.2021
!------------------------------------------------------------------------------
MODULE strings

USE global_std

CONTAINS

subroutine parse(str,delims,args,nargs)

! Parses the string 'str' into arguments args(1), ..., args(nargs) based on
! the delimiters contained in the string 'delims'. Preceding a delimiter in
! 'str' by a backslash (\) makes this particular instance not a delimiter.
! The integer output variable nargs contains the number of arguments found.

character(len=*) :: str,delims
character(len=len_trim(str)) :: strsav
character(len=mcl),dimension(:) :: args
integer(kind=ik) :: nargs

strsav=str
call compact(str)
na=size(args)
do i=1,na
  args(i)=' '
end do
nargs=0
lenstr=len_trim(str)
if(lenstr==0) return
k=0

do
   if(len_trim(str) == 0) exit
   nargs=nargs+1
   call split(str,delims,args(nargs))
   call removebksl(args(nargs))
end do
str=strsav

end subroutine parse

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

subroutine compact(str)

! Converts multiple spaces and tabs to single spaces; deletes control characters;
! removes initial spaces.

character(len=*):: str
character(len=1):: ch
character(len=len_trim(str)):: outstr

str=adjustl(str)
lenstr=len_trim(str)
outstr=' '
isp=0
k=0

do i=1,lenstr
  ch=str(i:i)
  ich=iachar(ch)

  select case(ich)

    case(9,32)     ! space or tab character
      if(isp==0) then
        k=k+1
        outstr(k:k)=' '
      end if
      isp=1

    case(33:)      ! not a space, quote, or control character
      k=k+1
      outstr(k:k)=ch
      isp=0

  end select

end do

str=adjustl(outstr)

end subroutine compact

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

subroutine removesp(str)

! Removes spaces, tabs, and control characters in string str

character(len=*):: str
character(len=1):: ch
character(len=len_trim(str))::outstr

str=adjustl(str)
lenstr=len_trim(str)
outstr=' '
k=0

do i=1,lenstr
  ch=str(i:i)
  ich=iachar(ch)
  select case(ich)
    case(0:32)  ! space, tab, or control character
         cycle
    case(33:)
      k=k+1
      outstr(k:k)=ch
  end select
end do

str=adjustl(outstr)

end subroutine removesp

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

subroutine value_dr(str,rnum,ios)

! Converts number string to a double precision real number

character(len=*)::str
real(rk)::rnum
integer :: ios

ilen=len_trim(str)
ipos=scan(str,'Ee')
if(.not.is_digit(str(ilen:ilen)) .and. ipos/=0) then
   ios=3
   return
end if
read(str,*,iostat=ios) rnum

end subroutine value_dr

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

subroutine value_sr(str,rnum,ios)

! Converts number string to a single precision real number

character(len=*)::str
real(rk) :: rnum
real(rk) :: rnumd

call value_dr(str,rnumd,ios)
if( abs(rnumd) > huge(rnum) ) then
  ios=15
  return
end if
if( abs(rnumd) < tiny(rnum) ) rnum=0.0_rk
rnum=rnumd

end subroutine value_sr

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

subroutine value_di(str,inum,ios)

! Converts number string to a double precision integer value

character(len=*)::str
integer(ik) :: inum
real(rk) :: rnum

call value_dr(str,rnum,ios)
if(abs(rnum)>huge(inum)) then
  ios=15
  return
end if
inum=nint(rnum,ik)

end subroutine value_di

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

subroutine value_si(str,inum,ios)

! Converts number string to a single precision integer value

character(len=*)::str
integer(ik) :: inum
real(rk) :: rnum

call value_dr(str,rnum,ios)
if(abs(rnum)>huge(inum)) then
  ios=15
  return
end if
inum=nint(rnum,ik)

end subroutine value_si

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

subroutine shiftstr(str,n)

! Shifts characters in in the string 'str' n positions (positive values
! denote a right shift and negative values denote a left shift). Characters
! that are shifted off the end are lost. Positions opened up by the shift
! are replaced by spaces.

character(len=*):: str

lenstr=len(str)
nabs=iabs(n)
if(nabs>=lenstr) then
  str=repeat(' ',lenstr)
  return
end if
if(n<0) str=str(nabs+1:)//repeat(' ',nabs)  ! shift left
if(n>0) str=repeat(' ',nabs)//str(:lenstr-nabs)  ! shift right
return

end subroutine shiftstr

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

subroutine insertstr(str,strins,loc)

! Inserts the string 'strins' into the string 'str' at position 'loc'.
! Characters in 'str' starting at position 'loc' are shifted right to
! make room for the inserted string. Trailing spaces of 'strins' are
! removed prior to insertion

character(len=*):: str,strins
character(len=len(str))::tempstr

lenstrins=len_trim(strins)
tempstr=str(loc:)
call shiftstr(tempstr,lenstrins)
tempstr(1:lenstrins)=strins(1:lenstrins)
str(loc:)=tempstr
return

end subroutine insertstr

!----------------------------------------------------------------------------*

subroutine write_dr(rnum,str,fmt)

! Writes double precision real number rnum to string str using format fmt

real(rk) :: rnum
character(len=*) :: str,fmt
character(len=80) :: formt

formt='('//trim(fmt)//')'
write(str,formt) rnum
str=adjustl(str)

end subroutine write_dr

!----------------------------------------------------------------------------*

subroutine write_sr(rnum,str,fmt)

! Writes single precision real number rnum to string str using format fmt

real(rk) :: rnum
character(len=*) :: str,fmt
character(len=80) :: formt

formt='('//trim(fmt)//')'
write(str,formt) rnum
str=adjustl(str)

end subroutine write_sr

!----------------------------------------------------------------------------*

subroutine write_di(inum,str,fmt)

! Writes double precision integer inum to string str using format fmt

integer(ik) :: inum
character(len=*) :: str,fmt
character(len=80) :: formt

formt='('//trim(fmt)//')'
write(str,formt) inum
str=adjustl(str)

end subroutine write_di

!----------------------------------------------------------------------------*

subroutine write_si(inum,str,fmt)

! Writes single precision integer inum to string str using format fmt

integer(ik) :: inum
character(len=*) :: str,fmt
character(len=80) :: formt

formt='('//trim(fmt)//')'
write(str,formt) inum
str=adjustl(str)

end subroutine write_si

!----------------------------------------------------------------------------*

subroutine trimzero(str)

! Deletes nonsignificant trailing zeroes from number string str. If number
! string ends in a decimal point, one trailing zero is added.

character(len=*) :: str
character :: ch
character(len=10) :: exp

ipos=scan(str,'eE')
if(ipos>0) then
   exp=str(ipos:)
   str=str(1:ipos-1)
endif
lstr=len_trim(str)
do i=lstr,1,-1
   ch=str(i:i)
   if(ch=='0') cycle
   if(ch=='.') then
      str=str(1:i)//'0'
      if(ipos>0) str=trim(str)//trim(exp)
      exit
   endif
   str=str(1:i)
   exit
end do
if(ipos>0) str=trim(str)//trim(exp)

end subroutine trimzero

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

function is_digit(ch) result(res)

! Returns .true. if ch is a digit (0,1,...,9) and .false. otherwise

character :: ch
logical :: res

select case(ch)
case('0':'9')
  res=.true.
case default
  res=.false.
end select
return

end function is_digit

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

subroutine split(str,delims,before,sep)

! Routine finds the first instance of a character from 'delims' in the
! the string 'str'. The characters before the found delimiter are
! output in 'before'. The characters after the found delimiter are
! output in 'str'. The optional output character 'sep' contains the
! found delimiter. A delimiter in 'str' is treated like an ordinary
! character if it is preceded by a backslash (\). If the backslash
! character is desired in 'str', then precede it with another backslash.

character(len=*) :: str,delims
character,optional :: sep
logical :: pres
character :: ch,cha
character(len=mcl) :: before

pres=present(sep)
str=adjustl(str)
call compact(str)
lenstr=len_trim(str)
if(lenstr == 0) return        ! string str is empty
k=0
ibsl=0                        ! backslash initially inactive
before=' '
do i=1,lenstr
   ch=str(i:i)
   if(ibsl == 1) then          ! backslash active
      k=k+1
      before(k:k)=ch
      ibsl=0
      cycle
   end if
   if(ch == '\\' ) then          ! backslash with backslash inactive
      k=k+1
      before(k:k)=ch
      ibsl=1
      cycle
   end if
   ipos=index(delims,ch)
   if(ipos == 0) then          ! character is not a delimiter
      k=k+1
      before(k:k)=ch
      cycle
   end if
   if(ch /= ' ') then          ! character is a delimiter that is not a space
      str=str(i+1:)
      if(pres) sep=ch
      exit
   end if
   cha=str(i+1:i+1)            ! character is a space delimiter
   iposa=index(delims,cha)
   if(iposa > 0) then          ! next character is a delimiter
      str=str(i+2:)
      if(pres) sep=cha
      exit
   else
      str=str(i+1:)
      if(pres) sep=ch
      exit
   end if
end do
if(i >= lenstr) str=''
str=adjustl(str)              ! remove initial spaces
return

end subroutine split

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

subroutine removebksl(str)

! Removes backslash (\) characters. Double backslashes (\\) are replaced
! by a single backslash.

character(len=*):: str
character(len=1):: ch
character(len=len_trim(str))::outstr

str=adjustl(str)
lenstr=len_trim(str)
outstr=' '
k=0
ibsl=0                        ! backslash initially inactive

do i=1,lenstr
  ch=str(i:i)
  if(ibsl == 1) then          ! backslash active
   k=k+1
   outstr(k:k)=ch
   ibsl=0
   cycle
  end if
  if(ch == '\\') then          ! backslash with backslash inactive
   ibsl=1
   cycle
  end if
  k=k+1
  outstr(k:k)=ch              ! non-backslash with backslash inactive
end do

str=adjustl(outstr)

end subroutine removebksl

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

END MODULE strings


text:
C FILE: FIB1.F
        SUBROUTINE FIB(A,N)
C
C CALCULATE FIRST N FIBONACCI NUMBERS
C
        INTEGER N
        REAL*8 A(N)
        DO I=1,N
            IF (I.EQ.1) THEN
                A(I) = 0.0D0
            ELSEIF (I.EQ.2) THEN
                A(I) = 1.0D0
            ELSE
                A(I) = A(I-1) + A(I-2)
            ENDIF
        ENDDO
        END
C END FILE FIB1.

max_stars_repo_name: kbroman/UW-Madison-swc-boot-camps

max_issues_repo_name: kbroman/UW-Madison-swc-boot-camps

max_forks_repo_name: kbroman/UW-Madison-swc-boot-camps


!-------------------------------------------------------------------------------
!> module Atmosphere / Physics Cloud Microphysics
!!
!! @par Description
!!          Container for mod_atmos_phy_mp
!!
!! @author Team SCALE
!!
!! @par History
!! @li      2014-05-04 (H.Yashiro)    [new]
!! @li      2015-09-08 (Y.Sato)       [add] Add ATMOS_PHY_MP_EVAPORATE
!! @li      2019-01-11 (S.Shima)      [add] namelist variable use_tendency_mp
!<
!-------------------------------------------------------------------------------
#include "inc_openmp.h"
module mod_atmos_phy_mp_vars
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision
  use scale_stdio
  use scale_prof
  use scale_grid_index
  use scale_tracer
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public procedure
  !
  public :: ATMOS_PHY_MP_vars_setup
  public :: ATMOS_PHY_MP_vars_fillhalo
  public :: ATMOS_PHY_MP_vars_restart_read
  public :: ATMOS_PHY_MP_vars_restart_write

  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  logical,               public :: ATMOS_PHY_MP_RESTART_OUTPUT       = .false.                !< output restart file?

  character(len=H_LONG), public :: ATMOS_PHY_MP_RESTART_IN_BASENAME  = ''                     !< basename of the restart file
  character(len=H_LONG), public :: ATMOS_PHY_MP_RESTART_OUT_BASENAME = ''                     !< basename of the output file
  character(len=H_MID),  public :: ATMOS_PHY_MP_RESTART_OUT_TITLE    = 'ATMOS_PHY_MP restart' !< title    of the output file
  character(len=H_MID),  public :: ATMOS_PHY_MP_RESTART_OUT_DTYPE    = 'DEFAULT'              !< REAL4 or REAL8

  real(RP), public, allocatable :: ATMOS_PHY_MP_DENS_t(:,:,:)    ! tendency DENS [kg/m3/s]
  real(RP), public, allocatable :: ATMOS_PHY_MP_MOMZ_t(:,:,:)    ! tendency MOMZ [kg/m2/s2]
  real(RP), public, allocatable :: ATMOS_PHY_MP_MOMX_t(:,:,:)    ! tendency MOMX [kg/m2/s2]
  real(RP), public, allocatable :: ATMOS_PHY_MP_MOMY_t(:,:,:)    ! tendency MOMY [kg/m2/s2]
  real(RP), public, allocatable :: ATMOS_PHY_MP_RHOT_t(:,:,:)    ! tendency RHOT [K*kg/m3/s]
  real(RP), public, allocatable :: ATMOS_PHY_MP_RHOQ_t(:,:,:,:)  ! tendency rho*QTRC [kg/kg/s]

  real(RP), public, allocatable :: ATMOS_PHY_MP_EVAPORATE(:,:,:) ! number concentration of evaporated cloud [/m3]
  real(RP), public, allocatable :: ATMOS_PHY_MP_SFLX_rain(:,:)   ! precipitation flux (liquid) [kg/m2/s]
  real(RP), public, allocatable :: ATMOS_PHY_MP_SFLX_snow(:,:)   ! precipitation flux (solid)  [kg/m2/s]

  logical,               public :: use_tendency_mp = .true. !< use tendency for cloud microphysics?

  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
  integer,                private, parameter :: VMAX = 2       !< number of the variables
  integer,                private, parameter :: I_SFLX_rain = 1
  integer,                private, parameter :: I_SFLX_snow = 2

  character(len=H_SHORT), private            :: VAR_NAME(VMAX) !< name  of the variables
  character(len=H_MID),   private            :: VAR_DESC(VMAX) !< desc. of the variables
  character(len=H_SHORT), private            :: VAR_UNIT(VMAX) !< unit  of the variables

  data VAR_NAME / 'SFLX_rain', &
                  'SFLX_snow'  /
  data VAR_DESC / 'precipitation flux (liquid)', &
                  'precipitation flux (solid)'   /
  data VAR_UNIT / 'kg/m2/s', &
                  'kg/m2/s' /

  !-----------------------------------------------------------------------------
contains
  !-----------------------------------------------------------------------------
  !> Setup
  subroutine ATMOS_PHY_MP_vars_setup
    use scale_process, only: &
       PRC_MPIstop
    use scale_const, only: &
       UNDEF => CONST_UNDEF
    implicit none

    NAMELIST / PARAM_ATMOS_PHY_MP_VARS / &
       ATMOS_PHY_MP_RESTART_IN_BASENAME,  &
       ATMOS_PHY_MP_RESTART_OUTPUT,       &
       ATMOS_PHY_MP_RESTART_OUT_BASENAME, &
       ATMOS_PHY_MP_RESTART_OUT_TITLE,    &
       ATMOS_PHY_MP_RESTART_OUT_DTYPE,    &
       use_tendency_mp

    integer :: ierr
    integer :: iv
    !---------------------------------------------------------------------------

    if( IO_L ) write(IO_FID_LOG,*)
    if( IO_L ) write(IO_FID_LOG,*) '++++++ Module[VARS] / Categ[ATMOS_PHY_MP] / Origin[SCALE-LES]'

    allocate( ATMOS_PHY_MP_DENS_t(KA,IA,JA)    )
    allocate( ATMOS_PHY_MP_MOMZ_t(KA,IA,JA)    )
    allocate( ATMOS_PHY_MP_MOMX_t(KA,IA,JA)    )
    allocate( ATMOS_PHY_MP_MOMY_t(KA,IA,JA)    )
    allocate( ATMOS_PHY_MP_RHOT_t(KA,IA,JA)    )
    allocate( ATMOS_PHY_MP_RHOQ_t(KA,IA,JA,QA) )
    allocate( ATMOS_PHY_MP_EVAPORATE(KA,IA,JA)    )
    ! tentative approach
    ATMOS_PHY_MP_DENS_t(:,:,:)   = 0.0_RP
    ATMOS_PHY_MP_MOMZ_t(:,:,:)   = 0.0_RP
    ATMOS_PHY_MP_MOMX_t(:,:,:)   = 0.0_RP
    ATMOS_PHY_MP_MOMY_t(:,:,:)   = 0.0_RP
    ATMOS_PHY_MP_RHOT_t(:,:,:)   = 0.0_RP
    ATMOS_PHY_MP_RHOQ_t(:,:,:,:) = 0.0_RP
    ATMOS_PHY_MP_EVAPORATE(:,:,:) = 0.0_RP

    allocate( ATMOS_PHY_MP_SFLX_rain(IA,JA) )
    allocate( ATMOS_PHY_MP_SFLX_snow(IA,JA) )
    ATMOS_PHY_MP_SFLX_rain(:,:) = UNDEF
    ATMOS_PHY_MP_SFLX_snow(:,:) = UNDEF

    !--- read namelist
    rewind(IO_FID_CONF)
    read(IO_FID_CONF,nml=PARAM_ATMOS_PHY_MP_VARS,iostat=ierr)
    if( ierr < 0 ) then !--- missing
       if( IO_L ) write(IO_FID_LOG,*) '*** Not found namelist. Default used.'
    elseif( ierr > 0 ) then !--- fatal error
       write(*,*) 'xxx Not appropriate names in namelist PARAM_ATMOS_PHY_MP_VARS. Check!'
       call PRC_MPIstop
    endif
    if( IO_LNML ) write(IO_FID_LOG,nml=PARAM_ATMOS_PHY_MP_VARS)

    if( IO_L ) write(IO_FID_LOG,*)
    if( IO_L ) write(IO_FID_LOG,*) '*** [ATMOS_PHY_MP] prognostic/diagnostic variables'
    if( IO_L ) write(IO_FID_LOG,'(1x,A,A15,A,A32,3(A))') &
               '***       |','VARNAME        ','|', 'DESCRIPTION                     ','[', 'UNIT            ',']'
    do iv = 1, VMAX
       if( IO_L ) write(IO_FID_LOG,'(1x,A,i3,A,A15,A,A32,3(A))') &
                  '*** NO.',iv,'|',VAR_NAME(iv),'|',VAR_DESC(iv),'[',VAR_UNIT(iv),']'
    enddo

    if( IO_L ) write(IO_FID_LOG,*)
    if ( ATMOS_PHY_MP_RESTART_IN_BASENAME /= '' ) then
       if( IO_L ) write(IO_FID_LOG,*) '*** Restart input?  : ', trim(ATMOS_PHY_MP_RESTART_IN_BASENAME)
    else
       if( IO_L ) write(IO_FID_LOG,*) '*** Restart input?  : NO'
    endif
    if (       ATMOS_PHY_MP_RESTART_OUTPUT             &
         .AND. ATMOS_PHY_MP_RESTART_OUT_BASENAME /= '' ) then
       if( IO_L ) write(IO_FID_LOG,*) '*** Restart output? : ', trim(ATMOS_PHY_MP_RESTART_OUT_BASENAME)
    else
       if( IO_L ) write(IO_FID_LOG,*) '*** Restart output? : NO'
       ATMOS_PHY_MP_RESTART_OUTPUT = .false.
    endif

    return
  end subroutine ATMOS_PHY_MP_vars_setup

  !-----------------------------------------------------------------------------
  !> HALO Communication
  subroutine ATMOS_PHY_MP_vars_fillhalo
    use scale_comm, only: &
       COMM_vars8, &
       COMM_wait
    implicit none
    !---------------------------------------------------------------------------

    call COMM_vars8( ATMOS_PHY_MP_SFLX_rain(:,:), 1 )
    call COMM_vars8( ATMOS_PHY_MP_SFLX_snow(:,:), 2 )
    call COMM_wait ( ATMOS_PHY_MP_SFLX_rain(:,:), 1 )
    call COMM_wait ( ATMOS_PHY_MP_SFLX_snow(:,:), 2 )

    return
  end subroutine ATMOS_PHY_MP_vars_fillhalo

  !-----------------------------------------------------------------------------
  !> Read restart
  subroutine ATMOS_PHY_MP_vars_restart_read
    use scale_fileio, only: &
       FILEIO_read
    use scale_les_statistics, only: &
       STAT_total
    implicit none

    real(RP) :: total
    !---------------------------------------------------------------------------

    if( IO_L ) write(IO_FID_LOG,*)
    if( IO_L ) write(IO_FID_LOG,*) '*** Input restart file (ATMOS_PHY_MP) ***'

    if ( ATMOS_PHY_MP_RESTART_IN_BASENAME /= '' ) then
       if( IO_L ) write(IO_FID_LOG,*) '*** basename: ', trim(ATMOS_PHY_MP_RESTART_IN_BASENAME)

       call FILEIO_read( ATMOS_PHY_MP_SFLX_rain(:,:),                                & ! [OUT]
                         ATMOS_PHY_MP_RESTART_IN_BASENAME, VAR_NAME(1), 'XY', step=1 ) ! [IN]
       call FILEIO_read( ATMOS_PHY_MP_SFLX_snow(:,:),                                & ! [OUT]
                         ATMOS_PHY_MP_RESTART_IN_BASENAME, VAR_NAME(2), 'XY', step=1 ) ! [IN]

       call ATMOS_PHY_MP_vars_fillhalo

       call STAT_total( total, ATMOS_PHY_MP_SFLX_rain(:,:), VAR_NAME(1) )
       call STAT_total( total, ATMOS_PHY_MP_SFLX_snow(:,:), VAR_NAME(2) )
    else
       if( IO_L ) write(IO_FID_LOG,*) '*** restart file for ATMOS_PHY_MP is not specified.'
    endif

    return
  end subroutine ATMOS_PHY_MP_vars_restart_read

  !-----------------------------------------------------------------------------
  !> Write restart
  subroutine ATMOS_PHY_MP_vars_restart_write
    use scale_time, only: &
       TIME_gettimelabel
    use scale_fileio, only: &
       FILEIO_write
    use scale_les_statistics, only: &
       STAT_total
    implicit none

    character(len=15)     :: timelabel
    character(len=H_LONG) :: basename

    real(RP) :: total
    !---------------------------------------------------------------------------

    if ( ATMOS_PHY_MP_RESTART_OUT_BASENAME /= '' ) then

       call TIME_gettimelabel( timelabel )
       write(basename,'(A,A,A)') trim(ATMOS_PHY_MP_RESTART_OUT_BASENAME), '_', trim(timelabel)

       if( IO_L ) write(IO_FID_LOG,*)
       if( IO_L ) write(IO_FID_LOG,*) '*** Output restart file (ATMOS_PHY_MP) ***'
       if( IO_L ) write(IO_FID_LOG,*) '*** basename: ', trim(basename)

       call ATMOS_PHY_MP_vars_fillhalo

       call STAT_total( total, ATMOS_PHY_MP_SFLX_rain(:,:), VAR_NAME(1) )
       call STAT_total( total, ATMOS_PHY_MP_SFLX_snow(:,:), VAR_NAME(1) )

       call FILEIO_write( ATMOS_PHY_MP_SFLX_rain(:,:), basename,       ATMOS_PHY_MP_RESTART_OUT_TITLE, & ! [IN]
                          VAR_NAME(1), VAR_DESC(1), VAR_UNIT(1), 'XY', ATMOS_PHY_MP_RESTART_OUT_DTYPE  ) ! [IN]
       call FILEIO_write( ATMOS_PHY_MP_SFLX_snow(:,:), basename,       ATMOS_PHY_MP_RESTART_OUT_TITLE, & ! [IN]
                          VAR_NAME(2), VAR_DESC(2), VAR_UNIT(2), 'XY', ATMOS_PHY_MP_RESTART_OUT_DTYPE  ) ! [IN]

    endif

    return
  end subroutine ATMOS_PHY_MP_vars_restart_write

end module mod_atmos_phy_mp_vars


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

!> Electronegativity equilibration charge model published in
!>
!> E. Caldeweyher, S. Ehlert, A. Hansen, H. Neugebauer, S. Spicher, C. Bannwarth
!> and S. Grimme, *J. Chem. Phys.*, **2019**, 150, 154122.
!> DOI: [10.1063/1.5090222](https://dx.doi.org/10.1063/1.5090222)
module multicharge_param_eeq2019
   use mctc_env, only : wp
   use mctc_io_symbols, only : to_number
   implicit none
   private

   public :: get_eeq_chi, get_eeq_eta, get_eeq_rad, get_eeq_kcn

   !> Element-specific electronegativity for the electronegativity equilibration charges.
   interface get_eeq_chi
      module procedure get_eeq_chi_sym
      module procedure get_eeq_chi_num
   end interface get_eeq_chi

   !> Element-specific chemical hardnesses for the electronegativity equilibration charges
   interface get_eeq_eta
      module procedure :: get_eeq_eta_sym
      module procedure :: get_eeq_eta_num
   end interface get_eeq_eta

   !> Element-specific CN scaling constant for the electronegativity equilibration charges.
   interface get_eeq_kcn
      module procedure :: get_eeq_kcn_sym
      module procedure :: get_eeq_kcn_num
   end interface get_eeq_kcn

   !> Element-specific charge widths for the electronegativity equilibration charges.
   interface get_eeq_rad
      module procedure :: get_eeq_rad_sym
      module procedure :: get_eeq_rad_num
   end interface get_eeq_rad

   !> Maximum atomic number allowed in EEQ calculations
   integer, parameter :: max_elem = 86

   !> Element-specific electronegativity for the electronegativity equilibration charges.
   real(wp), parameter :: eeq_chi(max_elem) = [&
      & 1.23695041_wp, 1.26590957_wp, 0.54341808_wp, 0.99666991_wp, 1.26691604_wp, &
      & 1.40028282_wp, 1.55819364_wp, 1.56866440_wp, 1.57540015_wp, 1.15056627_wp, &
      & 0.55936220_wp, 0.72373742_wp, 1.12910844_wp, 1.12306840_wp, 1.52672442_wp, &
      & 1.40768172_wp, 1.48154584_wp, 1.31062963_wp, 0.40374140_wp, 0.75442607_wp, &
      & 0.76482096_wp, 0.98457281_wp, 0.96702598_wp, 1.05266584_wp, 0.93274875_wp, &
      & 1.04025281_wp, 0.92738624_wp, 1.07419210_wp, 1.07900668_wp, 1.04712861_wp, &
      & 1.15018618_wp, 1.15388455_wp, 1.36313743_wp, 1.36485106_wp, 1.39801837_wp, &
      & 1.18695346_wp, 0.36273870_wp, 0.58797255_wp, 0.71961946_wp, 0.96158233_wp, &
      & 0.89585296_wp, 0.81360499_wp, 1.00794665_wp, 0.92613682_wp, 1.09152285_wp, &
      & 1.14907070_wp, 1.13508911_wp, 1.08853785_wp, 1.11005982_wp, 1.12452195_wp, &
      & 1.21642129_wp, 1.36507125_wp, 1.40340000_wp, 1.16653482_wp, 0.34125098_wp, &
      & 0.58884173_wp, 0.68441115_wp, 0.56999999_wp, 0.56999999_wp, 0.56999999_wp, &
      & 0.56999999_wp, 0.56999999_wp, 0.56999999_wp, 0.56999999_wp, 0.56999999_wp, &
      & 0.56999999_wp, 0.56999999_wp, 0.56999999_wp, 0.56999999_wp, 0.56999999_wp, &
      & 0.56999999_wp, 0.87936784_wp, 1.02761808_wp, 0.93297476_wp, 1.10172128_wp, &
      & 0.97350071_wp, 1.16695666_wp, 1.23997927_wp, 1.18464453_wp, 1.14191734_wp, &
      & 1.12334192_wp, 1.01485321_wp, 1.12950808_wp, 1.30804834_wp, 1.33689961_wp, &
      & 1.27465977_wp]

   !> Element-specific chemical hardnesses for the electronegativity equilibration charges.
   real(wp), parameter :: eeq_eta(max_elem) = [&
      &-0.35015861_wp, 1.04121227_wp, 0.09281243_wp, 0.09412380_wp, 0.26629137_wp, &
      & 0.19408787_wp, 0.05317918_wp, 0.03151644_wp, 0.32275132_wp, 1.30996037_wp, &
      & 0.24206510_wp, 0.04147733_wp, 0.11634126_wp, 0.13155266_wp, 0.15350650_wp, &
      & 0.15250997_wp, 0.17523529_wp, 0.28774450_wp, 0.42937314_wp, 0.01896455_wp, &
      & 0.07179178_wp,-0.01121381_wp,-0.03093370_wp, 0.02716319_wp,-0.01843812_wp, &
      &-0.15270393_wp,-0.09192645_wp,-0.13418723_wp,-0.09861139_wp, 0.18338109_wp, &
      & 0.08299615_wp, 0.11370033_wp, 0.19005278_wp, 0.10980677_wp, 0.12327841_wp, &
      & 0.25345554_wp, 0.58615231_wp, 0.16093861_wp, 0.04548530_wp,-0.02478645_wp, &
      & 0.01909943_wp, 0.01402541_wp,-0.03595279_wp, 0.01137752_wp,-0.03697213_wp, &
      & 0.08009416_wp, 0.02274892_wp, 0.12801822_wp,-0.02078702_wp, 0.05284319_wp, &
      & 0.07581190_wp, 0.09663758_wp, 0.09547417_wp, 0.07803344_wp, 0.64913257_wp, &
      & 0.15348654_wp, 0.05054344_wp, 0.11000000_wp, 0.11000000_wp, 0.11000000_wp, &
      & 0.11000000_wp, 0.11000000_wp, 0.11000000_wp, 0.11000000_wp, 0.11000000_wp, &
      & 0.11000000_wp, 0.11000000_wp, 0.11000000_wp, 0.11000000_wp, 0.11000000_wp, &
      & 0.11000000_wp,-0.02786741_wp, 0.01057858_wp,-0.03892226_wp,-0.04574364_wp, &
      &-0.03874080_wp,-0.03782372_wp,-0.07046855_wp, 0.09546597_wp, 0.21953269_wp, &
      & 0.02522348_wp, 0.15263050_wp, 0.08042611_wp, 0.01878626_wp, 0.08715453_wp, &
      & 0.10500484_wp]

   !> Element-specific CN scaling constant for the electronegativity equilibration charges.
   real(wp), parameter :: eeq_kcn(max_elem) = [&
      & 0.04916110_wp, 0.10937243_wp,-0.12349591_wp,-0.02665108_wp,-0.02631658_wp, &
      & 0.06005196_wp, 0.09279548_wp, 0.11689703_wp, 0.15704746_wp, 0.07987901_wp, &
      &-0.10002962_wp,-0.07712863_wp,-0.02170561_wp,-0.04964052_wp, 0.14250599_wp, &
      & 0.07126660_wp, 0.13682750_wp, 0.14877121_wp,-0.10219289_wp,-0.08979338_wp, &
      &-0.08273597_wp,-0.01754829_wp,-0.02765460_wp,-0.02558926_wp,-0.08010286_wp, &
      &-0.04163215_wp,-0.09369631_wp,-0.03774117_wp,-0.05759708_wp, 0.02431998_wp, &
      &-0.01056270_wp,-0.02692862_wp, 0.07657769_wp, 0.06561608_wp, 0.08006749_wp, &
      & 0.14139200_wp,-0.05351029_wp,-0.06701705_wp,-0.07377246_wp,-0.02927768_wp, &
      &-0.03867291_wp,-0.06929825_wp,-0.04485293_wp,-0.04800824_wp,-0.01484022_wp, &
      & 0.07917502_wp, 0.06619243_wp, 0.02434095_wp,-0.01505548_wp,-0.03030768_wp, &
      & 0.01418235_wp, 0.08953411_wp, 0.08967527_wp, 0.07277771_wp,-0.02129476_wp, &
      &-0.06188828_wp,-0.06568203_wp,-0.11000000_wp,-0.11000000_wp,-0.11000000_wp, &
      &-0.11000000_wp,-0.11000000_wp,-0.11000000_wp,-0.11000000_wp,-0.11000000_wp, &
      &-0.11000000_wp,-0.11000000_wp,-0.11000000_wp,-0.11000000_wp,-0.11000000_wp, &
      &-0.11000000_wp,-0.03585873_wp,-0.03132400_wp,-0.05902379_wp,-0.02827592_wp, &
      &-0.07606260_wp,-0.02123839_wp, 0.03814822_wp, 0.02146834_wp, 0.01580538_wp, &
      &-0.00894298_wp,-0.05864876_wp,-0.01817842_wp, 0.07721851_wp, 0.07936083_wp, &
      & 0.05849285_wp]

   !> Element-specific charge widths for the electronegativity equilibration charges.
   real(wp), parameter :: eeq_rad(max_elem) = [&
      & 0.55159092_wp, 0.66205886_wp, 0.90529132_wp, 1.51710827_wp, 2.86070364_wp, &
      & 1.88862966_wp, 1.32250290_wp, 1.23166285_wp, 1.77503721_wp, 1.11955204_wp, &
      & 1.28263182_wp, 1.22344336_wp, 1.70936266_wp, 1.54075036_wp, 1.38200579_wp, &
      & 2.18849322_wp, 1.36779065_wp, 1.27039703_wp, 1.64466502_wp, 1.58859404_wp, &
      & 1.65357953_wp, 1.50021521_wp, 1.30104175_wp, 1.46301827_wp, 1.32928147_wp, &
      & 1.02766713_wp, 1.02291377_wp, 0.94343886_wp, 1.14881311_wp, 1.47080755_wp, &
      & 1.76901636_wp, 1.98724061_wp, 2.41244711_wp, 2.26739524_wp, 2.95378999_wp, &
      & 1.20807752_wp, 1.65941046_wp, 1.62733880_wp, 1.61344972_wp, 1.63220728_wp, &
      & 1.60899928_wp, 1.43501286_wp, 1.54559205_wp, 1.32663678_wp, 1.37644152_wp, &
      & 1.36051851_wp, 1.23395526_wp, 1.65734544_wp, 1.53895240_wp, 1.97542736_wp, &
      & 1.97636542_wp, 2.05432381_wp, 3.80138135_wp, 1.43893803_wp, 1.75505957_wp, &
      & 1.59815118_wp, 1.76401732_wp, 1.63999999_wp, 1.63999999_wp, 1.63999999_wp, &
      & 1.63999999_wp, 1.63999999_wp, 1.63999999_wp, 1.63999999_wp, 1.63999999_wp, &
      & 1.63999999_wp, 1.63999999_wp, 1.63999999_wp, 1.63999999_wp, 1.63999999_wp, &
      & 1.63999999_wp, 1.47055223_wp, 1.81127084_wp, 1.40189963_wp, 1.54015481_wp, &
      & 1.33721475_wp, 1.57165422_wp, 1.04815857_wp, 1.78342098_wp, 2.79106396_wp, &
      & 1.78160840_wp, 2.47588882_wp, 2.37670734_wp, 1.76613217_wp, 2.66172302_wp, &
      & 2.82773085_wp]

contains

!> Get electronegativity for species with a given symbol
elemental function get_eeq_chi_sym(symbol) result(chi)

   !> Element symbol
   character(len=*), intent(in) :: symbol

   !> electronegativity
   real(wp) :: chi

   chi = get_eeq_chi(to_number(symbol))

end function get_eeq_chi_sym

!> Get electronegativity for species with a given atomic number
elemental function get_eeq_chi_num(number) result(chi)

   !> Atomic number
   integer, intent(in) :: number

   !> electronegativity
   real(wp) :: chi

   if (number > 0 .and. number <= size(eeq_chi, dim=1)) then
      chi = eeq_chi(number)
   else
      chi = -1.0_wp
   end if

end function get_eeq_chi_num

!> Get hardness for species with a given symbol
elemental function get_eeq_eta_sym(symbol) result(eta)

   !> Element symbol
   character(len=*), intent(in) :: symbol

   !> hardness
   real(wp) :: eta

   eta = get_eeq_eta(to_number(symbol))

end function get_eeq_eta_sym

!> Get hardness for species with a given atomic number
elemental function get_eeq_eta_num(number) result(eta)

   !> Atomic number
   integer, intent(in) :: number

   !> hardness
   real(wp) :: eta

   if (number > 0 .and. number <= size(eeq_eta, dim=1)) then
      eta = eeq_eta(number)
   else
      eta = -1.0_wp
   end if

end function get_eeq_eta_num

!> Get CN scaling for species with a given symbol
elemental function get_eeq_kcn_sym(symbol) result(kcn)

   !> Element symbol
   character(len=*), intent(in) :: symbol

   !> CN scaling
   real(wp) :: kcn

   kcn = get_eeq_kcn(to_number(symbol))

end function get_eeq_kcn_sym

!> Get CN scaling for species with a given atomic number
elemental function get_eeq_kcn_num(number) result(kcn)

   !> Atomic number
   integer, intent(in) :: number

   !> CN scaling
   real(wp) :: kcn

   if (number > 0 .and. number <= size(eeq_kcn, dim=1)) then
      kcn = eeq_kcn(number)
   else
      kcn = -1.0_wp
   end if

end function get_eeq_kcn_num

!> Get charge width for species with a given symbol
elemental function get_eeq_rad_sym(symbol) result(rad)

   !> Element symbol
   character(len=*), intent(in) :: symbol

   !> charge width
   real(wp) :: rad

   rad = get_eeq_rad(to_number(symbol))

end function get_eeq_rad_sym

!> Get charge width for species with a given atomic number
elemental function get_eeq_rad_num(number) result(rad)

   !> Atomic number
   integer, intent(in) :: number

   !> Charge width
   real(wp) :: rad

   if (number > 0 .and. number <= size(eeq_rad, dim=1)) then
      rad = eeq_rad(number)
   else
      rad = -1.0_wp
   end if

end function get_eeq_rad_num

end module multicharge_param_eeq2019

max_stars_repo_name: loriab/multicharge

max_issues_repo_name: loriab/multicharge

max_forks_repo_name: loriab/multicharge


text:
Paul David Terry (B.A. 2001) is the first Director of the UCSD Pep Band (a contending ensemble during Picnic Days Battle of the Bands).

Paul studied music and biology at the University of California, Davis.  While at Davis he performed with the UCDSO, UCD Chorus and Chamber Singers, as well as, moon lighted on occasion with the Concert Band, wind ensemble, early music ensemble, various theater productions  and of course the CalAggie Marching Banduh!

He proceeded Trevor Benson as the academic peer advisor for the Department of Music from 1998  2001, sharing this position with Marnie Efishoff (1998  1999) and Mia Anderson (1999  2001)  He served as the Emcee of Blue Polo Shirt Day and Copreseident of Impresario Society from 1998  2001.

Paul frequently participates with the Sethsingers! at the annual UCD Alumni Chorus concerts and serves as the President of the CAAA SAC (Sacramento Aggies Chapter) and the CAAA LGBTA (Lesbian, Gay, Bisexual, Transgender, & Ally) Chapter.

In his spare time he enjoys brewing ginger beer, hiking to California waterfalls, calm nights on the beach with a setting wet moon and glass of wine, sorrel, tea lattes with lapsang souchong, and solo car rides performing his one man shows of Porgy and Bess, Salome, Fiddler on the Roof, Verdis Macbeth, and Les Mis.

Related Links:
  UCSD Pep Band http://www.ucsdpepband.com
  CalAggie Alumni Marching Banduh! & Kiddie Band http://camb.ucdavis.edu/camb/cambaa/index.php

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
module mathability

  implicit none
  public square

contains

  function square(array) result (squared)
    integer, intent (in) :: array(:)
    integer :: squared(size(array))
    integer :: i

    do i = 1, size(array)
       squared(i) = array(i)*array(i)
    enddo
  end function square

end module mathability

max_stars_repo_name: mdpiper/fortran-explorations

max_issues_repo_name: mdpiper/fortran-explorations

max_forks_repo_name: mdpiper/fortran-explorations


! MIT License
!
! Copyright (c) 2020 SHEMAT-Suite
!
! Permission is hereby granted, free of charge, to any person obtaining a copy
! of this software and associated documentation files (the "Software"), to deal
! in the Software without restriction, including without limitation the rights
! to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
! copies of the Software, and to permit persons to whom the Software is
! furnished to do so, subject to the following conditions:
!
! The above copyright notice and this permission notice shall be included in all
! copies or substantial portions of the Software.
!
! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
! IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
! AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
! SOFTWARE.

!>    @brief writes covariance matrix
!>    @param[in] ident inverse iteration counter
!>    @param[in] ismpl local sample index
      SUBROUTINE write_covariances(ident,ismpl)
        use arrays
        use mod_genrl
        use mod_genrlc
        use mod_inverse
        use mod_data
        use mod_linfos
        IMPLICIT NONE
        integer :: ismpl
        integer :: i, j

        INTEGER ident, lout
        INTEGER i1_, i2_, i3_, i4_

        CHARACTER filename*80, snumber*8
        INTEGER lblank
        EXTERNAL lblank

        IF ((covar<=0) .OR. write_disable) RETURN
#ifdef NOOUT
        RETURN
#endif

        IF (mpara<1) THEN
          WRITE(*,*) 'error : in "write_covariance" need "mpara">0 !'
          STOP
        END IF

        CALL omp_new_file_handler(lout,16)

        CALL chln(project,i1_,i2_)
        IF (ident>=0) THEN
          WRITE(snumber,'(1I7)') ident
        ELSE IF (ident==-1) THEN
          WRITE(snumber,'(A8)') 'final'
        ELSE IF (ident==-2) THEN
          WRITE(snumber,'(A8)') 'debug'
        ELSE IF (ident==-3) THEN
          WRITE(snumber,'(A8)') 'ens_mean'
        ELSE IF (ident==-4) THEN
          WRITE(snumber,'(A8)') 'mean'
        ELSE IF (ident==-5) THEN
          WRITE(snumber,'(A8)') 'ens_mean'
        END IF
        CALL chln(snumber,i3_,i4_)
        filename = project(i1_:i2_) // '_cov-p_' // &
          snumber(i3_:i4_) // '.dat'

        OPEN(lout,file=filename,status='unknown',blank='null')

        IF (linfos(2)>=0) THEN
          WRITE(*,'(3A)') '  [W] : Parameter-Covariances to "', &
            filename(1:lblank(filename)), '"'
        END IF

        IF (mpara>200) WRITE(*,*) &
          'warning: line to short in "write_covarinaces" !'

        DO i = 1, mpara
          WRITE(lout,'(200(e15.7,1x))') (covar_p(i,j),j=1,mpara)
        END DO

        CLOSE(lout)
        CALL omp_del_file_handler(lout)

        RETURN
      END


text:
      subroutine gensmethh2o
!***********************************************************************
!  Copyright, 2004,  The  Regents  of the  University of California.
!  This program was prepared by the Regents of the University of
!  California at Los Alamos National Laboratory (the University) under
!  contract No. W-7405-ENG-36 with the U.S. Department of Energy (DOE).
!  All rights in the program are reserved by the DOE and the University.
!  Permission is granted to the public to copy and use this software
!  without charge, provided that this Notice and any statement of
!  authorship are reproduced on all copies. Neither the U.S. Government
!  nor the University makes any warranty, express or implied, or
!  assumes any liability or responsibility for the use of this software.
!**********************************************************************
!D1
!D1  PURPOSE
!D1
!D1  This subroutine solves the non-isothermal water-methane equations
!D1  full jacobian (unsymmetric, 3n by 3n).
!D1
!***********************************************************************
!D2
!D2 REVISION HISTORY
!D2
!D2 FEHM Version 2.30
!D2
!D2 $Log:   /pvcs.config/fehm90/src/gensmethh2o.f_a  $
!D2
!***********************************************************************
!D3
!D3  REQUIREMENTS TRACEABILITY
!D3
!D3  2.3.2 Heat- and mass-transfer equations
!D3  2.5.2 Solve nonlinear equation set at each time step
!D3
!***********************************************************************
!D4
!D4  SPECIAL COMMENTS AND REFERENCES
!D4
!D4 Requirements from SDN: 10086-RD-2.20-00
!D4   SOFTWARE REQUIREMENTS DOCUMENT (RD) for the
!D4   FEHM Application Version 2.20
!D4
!***********************************************************************

      use comai
      use comgi
      use comfi
      use comei
      use comdi
      use comci
      use combi
      use comhi
      use davidi
      use comdti
      use comcouple
      use commeth
      implicit none

      integer nmatd, index(6)
      integer i, j, icoupl_iter, id, idl, idofm
      integer jj,ii,idiag
      integer iprint,inorm
      integer neqp1, nrhs1, nrhs2, nrhs3, nsizea, nsizea1
      real*8, allocatable :: dumz(:)
      real*8, allocatable :: dumn(:)
      real*8, allocatable :: sto5(:)
      real*8  facr, fdum2, tollr, tolls
      parameter (iprint=0,inorm = 1)

      neqp1=neq+1
      nmatd=nelm(neqp1)-neqp1
      nsizea=idof*idof*nmatd
      do i=1,nsizea
         a(i)=0.0
      enddo

      do i=1,neq
       do j=1,idof
         bp(i+nrhs(j))=0.0
       enddo
      enddo

c
c    call to icectr to get permeability reduction factor
c
        call icectr(10,0)
c
c    call to ther_meth_h2o to get rel perms for water and methane
c
        call ther_meth_h2o(0,0)
c
c    hydrate componet
c
c  no rearrangement needed for hydrate since full 6 by 6
c  ie geneq_hyd assumes structure
c
c     add contribution heat-mass of hydrate
c     C of water equation
c     C of energy equation
c     care not to use derivative terms already used by water
c
         call ther_meth_h2o(3,0)
c
      do id = 1,neq
          call geneq_hyd(id)
      enddo
c
c
c     water component (gaz)
c
c
c     add contribution heat-mass of water
c     C of water equation
c     C of energy equation
c
         call ther_meth_h2o(1,0)
c
      do id = 1,neq
          call geneq_h2o(id)
      enddo
c
c     add contribution methane
c     C of methane equation
c     C of energy equation (contribution)
c     new 10-28-2002 gaz
         call ther_meth_h2o(2,0)
c
c
c     decide on equation type)
c
      do id = 1,neq
          call geneq_meth(id)
      enddo
c
c    call to ther_meth_h2o to release memory for rel perms
c
        call ther_meth_h2o(-1,0)
c
c   enforce equilibrium conditions as necessary
c
      call methh2o_combine(1,idofm)

c     print out before manipulating or normalizing arrays
      if(iprint.ne.0) then
       id =1
       write(*,*) 'node = ',id
       write(*,*) 'residuals'
       do jj=1,idof
        write(*,*) jj
        write(*,*) bp(id+nrhs(jj))
        write(*,*) 'diagonal of a'
	idiag=nelmdg(id)-neqp1
        write(*,200) (a(idiag+nmat(6*(jj-1)+ii)),ii=1,6)
       enddo
       write(*,*)
       id =neq
       write(*,*) 'node = ',id
       write(*,*) 'residuals'
       do jj=1,idof
        write(*,*) jj
        write(*,*) bp(id+nrhs(jj))
        write(*,*) 'diagonal of a'
	idiag=nelmdg(id)-neqp1
        write(*,200) (a(idiag+nmat(6*(jj-1)+ii)),ii=1,6)
       enddo
       if(iprint.eq.2) stop
      endif
200    format(1x,1p,6(g14.4))

      if(islord.ne.0) call switch(nmat,nmatb,islord,idofm,1,nsizea1)
      if(islord.ne.0) call switchb(bp,nrhs,nrhsb,islord,idofm,1,neq)

c
c     normalize the equations
c
      allocate(dumn(108))

c  gaz 2-27-03
       if(inorm.ne.0) then
        call normal_dof(neq,a,bp,nelm,nmat,nrhs,nelmdg
     &      ,index,idofm,dumn(1),dumn(17),dumn(33),0,fdum2)
       else
	 fdum2 = 0.0
         do i=1,neq
          do j=1,idofm
           fdum2 = fdum2+bp(i+nrhs(j))*bp(i+nrhs(j))
          enddo
         enddo
       endif
c     call mmrelblk ("dumn", "dpdp", ipdumn, icode)
      deallocate(dumn)
c
c     return and/or stop if singular
c
      if(index(1).lt.0 ) then
         if (iout .ne. 0) then
            write(iout,*) '   '
            write(iout, 100)
            write(iout,*) '   '
            write(iout,'(i8,3g12.3)')
     &           abs(nmatb(1)),(cord(abs(nmatb(1)),i),i=1,3)
         end if
         if(iptty.gt.0)then
            write(iptty,*) '   '
            write(iptty, 100)
            write(iptty,*) '   '
            write(iptty,'(i8,3g12.3)')
     &           abs(nmatb(1)),(cord(abs(nmatb(1)),i),i=1,3)
         endif
         iad=maxit
         return
      endif
 100   format('* singular matrix found during normalization *')

c     find residual
c
      fdum=sqrt(fdum2)
      if(fdum.eq.0.0) go to 999
      mink=n
      if(iad.eq.0) then
         f0=max(fdum*epe,tmch)
      endif
      if(fdum1.lt.0.0.and.iad.ne.0) then
         do i=1,neq
          do j=1,idofm
            if(abs(bp(i+nrhs(j))).gt.tmch) go to 99
          enddo
         enddo
         fdum=-1.0
         go to 999
 99      continue
         f0=-1.0
      endif
      if(f0.gt.0.0) then
         if(fdum.le.f0.and.iad.ne.0) goto 999
         facr=1.0
         tolls=max(facr*f0,min(g1*fdum,g2*fdum2))
         tollr=tolls*g3
         if(g3*tmch.gt.tollr) tollr=g3*tmch
       else if(f0.le.0.0) then
         tollr=tmch *g3
       endif

c
c     set maximum iterations in solve_new
c
      iter=maxsolve
      call storage_derivatives(1,1)
      allocate(dumz(neq*5*4))
c
c     set maximum recoupling iterations if necesary
c
      if(icoupl.gt.0) then
        icoupl_iter = icoupl
      else if(icoupl.lt.0) then
        icoupl_iter =
     &   max(float(abs(icoupl)),
     &   (float(iad)/float(abs(maxit)))*float(abs(2*icoupl)))
      endif
      if(irdof.eq.0) then
c
c     full dof solution
c

             if(gdpm_flag.eq.0) then
              call solve_new(neq,a,b,bp,nmat,nb,nrhs,nelm,nop,north
     &             ,tollr,irb,iirb,npvt,gmres,dumz,piv
     &             ,h,c,ss,g,y,iter,iback,idofm,iptty,maxor,accm)
             else
             call solve_dual(neq_primary,neq,a,b,bp,nmat,nb,nrhs
     &            ,nelm,nelm_primary,nop,north,tollr,irb,iirb
     &            ,npvt,gmres,dumz,piv
     &            ,h,c,ss,g,y,iter,iback,idofm,iptty,maxor
     &            ,igdpm,maxgdpmlayers,ngdpm_layers,nelmdg,accm
     &            ,mdof_sol)
             endif

      else if(irdof.eq.-idofm) then
c     reduced degree of freedom with full grmres
         call solve_rdof(neq,a,b,bp,nmat,nb,nrhs,nelm,nop,north
     &        ,tollr,irb,iirb,npvt,gmres,dumz,piv
     &        ,h,c,ss,g,y,iter,iback,-idofm,iptty
     &        ,maxor,icoupl_iter,tollr,overf,accm)
      else
c
c     approximate solution
c
         allocate(sto5(neq*6))
         call rdof_new(neq,a,b,bp,nmat,nb,nrhs,nelm,nop,north
     &        ,tollr,irb,iirb,npvt,gmres,dumz,piv
     &        ,h,c,ss,g,y,iter,iback,idofm,iptty,maxor
     &        ,overf,irdof,icoupl_iter,0,sto5,accm)
         deallocate(sto5)
      endif
      itert=itert+iter
      itotals=itotals+iter
      minkt=minkt+mink
      deallocate(dumz)
      call storage_derivatives(0,1)

 999  continue

      if(islord.ne.0) call switch(nmat,nmatb,islord,idofm,2,nsizea1)
      if(islord.ne.0) call switchb(bp,nrhs,nrhsb,islord,idofm,2,neq)

c   enforce thermal equilibrium( set t2 = t1)

      call methh2o_combine(2,idofm)

      return
      end

max_stars_repo_name: satkarra/FEHM

max_issues_repo_name: satkarra/FEHM

max_forks_repo_name: satkarra/FEHM


text:
subroutine gmres(x,d,dold,w,bg,dg,hg,ien,fext,id,&
                 ne_local,ien_local,node_local,nn_local,&
                 global_com,nn_global_com,local_com,nn_local_com,send_address,ad_length,&
                 fden,fvis,I_fluid,rngface)
    use fluid_variables, only: nsd,nn,ne,nen,ndf,inner,outer,neface
    use solid_variables, only: nn_solid
    use mpi_variables
    use pml_variables
    implicit none
    include 'mpif.h'
	real* 8 x(nsd,nn)
    integer id(ndf,nn)
	real* 8 d(ndf,nn), dold(ndf,nn),hg(ne),fext(ndf,nn)
    integer ien(nen,ne)
	real* 8 bg(ndf*(nn+nn_PML)), dg(ndf*(nn+nn_PML)), w(ndf*(nn+nn_PML))
	real* 8 Hm(inner+1,inner) !Henssenberg matrix
!=========================================================
! reduce dimension to save memory
	real* 8 Vm(ndf*(nn_local+nn_PML_local), inner+1) ! Krylov space matrix
!        real* 8 Vm(ndf*nn, inner+1) ! Krylov space matrix
!=========================================================
    integer i,j,iouter,icount,INFO
	real* 8 e1(inner+1)
	real* 8 x0(ndf*(nn+nn_PML))
	real* 8 beta(inner+1)
	real* 8 eps
	real* 8 r0(ndf*(nn+nn_PML))
	real* 8 rnorm, rnorm0, err
    real* 8 dv(ndf*(nn+nn_PML))
	real* 8 Vy(ndf*(nn+nn_PML))
    real* 8 vloc(ndf,nn+nn_PML), avloc(ndf*(nn+nn_PML))
!==============================
! reduce dimension save memory
	real* 8 temp
!==============================
    character(1) TRAN
	real* 8 workls(2*inner)
	real* 8 av_tmp(ndf,nn+nn_PML)
    integer rngface(neface,ne)
!---------------------------------------
    real(8) fden(nn)
    real(8) fvis(nn)
    real(8) I_fluid(nn)
!============================
! MPI varibalbes
    integer ne_local ! # of element on each processor
    integer ien_local(ne_local) ! subregion-->wholeregion element index
    integer ie_local ! loop parameter
    integer node_local(nn_local)
    integer nn_local
    integer jcount
    integer kcount
    integer node
    integer nn_global_com
    integer global_com(nn_global_com)  ! global node index for communication
    integer nn_local_com
    integer local_com(nn_local_com)  ! local index in the communication region on each processor
    integer ad_length
    integer send_address(ad_length,2)
    integer flag
    real(8) dg_sent(ndf*(nn+nn_PML))
    real(8) space1(inner)
    real(8) space2(inner)
    integer time_arrary_0(8)
    integer time_arrary_1(8)
    real(8) start_time
    real(8) end_time
    real(8) linerr
!---------------------------------------------
    eps = 1.0e-6
    linerr = 1.0e-6
    e1(:) = 0.0
    e1(1) = 1.0
    x0(:) = 0.0
    iouter = 1
    r0(:) = bg(:)
    TRAN = 'N'
    av_tmp(:,:) = 0
    avloc(:) = 0
!	w(:) = 1
    call getnormpml_pa(r0,ndf,nn,node_local,nn_local,rnorm0)
    rnorm = sqrt(rnorm0)

!!!!!!!!!!!!!!!start outer loop!!!!!!!!!
do 111, while((iouter .le. outer) .and. (rnorm .ge. linerr))

!==============================
    Vm(:,:) = 0.0d0
! Already is a local variable defined on each proc
!===============================
    do icount = 1, nn_local
        node=node_local(icount)
        do jcount=1,ndf
            Vm((icount-1)*ndf+jcount,1) = r0((node-1)*ndf+jcount)/rnorm
        enddo
        if (seqcPML(node) > 0) then
            do jcount=1,ndf
                Vm((nn_local+seqcPMLlocal(icount)-1)*ndf+jcount,1) = r0((nn+seqcPML(node)-1)*ndf+jcount)/rnorm
            enddo
        endif
    enddo ! get V1

    beta(:) = rnorm*e1(:) ! get beta*e1
    Hm(:,:) = 0.0
!!!!!!!!!!!!!!!!start inner loop!!!!!!!!!!!!!
    do j=1,inner
        do icount=1, nn_local
            node=node_local(icount)
            do jcount=1,ndf
                dv((node-1)*ndf+jcount) = eps/w((node-1)*ndf+jcount)*Vm((icount-1)*ndf+jcount,j)
            enddo
            if (seqcPML(node) > 0) then
                do jcount=1,ndf
                    dv((nn+seqcPML(node)-1)*ndf+jcount) = eps*Vm((nn_local+seqcPMLlocal(icount)-1)*ndf+jcount,j) &
                                                          /w((nn+seqcPML(node)-1)*ndf+jcount)
                enddo
            endif
        enddo!!!!!!!!!!calcule eps*inv(P)*V1
!===============================================
!		vloc(:,:) = 0.0d0
! Clear the matrix locally to avoid too many loops
        do icount=1,nn_local
            node=node_local(icount)
            vloc(1:ndf,node)=0.0d0
            if (seqcPML(node) > 0) then
                vloc(1:ndf,nn+seqcPML(node))=0.0d0
            endif
        enddo
! clear all the processor internal nodes
        do icount=1,nn_local_com
            node=global_com(local_com(icount))
            vloc(1:ndf,node)=0.0d0
            if (seqcPML(node) > 0) then
                vloc(1:ndf,nn+seqcPML(node))=0.0d0
            endif
        enddo
! clear the processor boundary nodes
!===============================================
        call equalpml_pa(dv,vloc,ndf,nn,node_local,nn_local)
!============================
        do icount=1,nn_local
            node=node_local(icount)
            vloc(1:ndf,node)=vloc(1:ndf,node)+d(1:ndf,node)
            if (seqcPML(node) > 0) then
                vloc(1:ndf,nn+seqcPML(node))=vloc(1:ndf,nn+seqcPML(node))+qv(1:ndf,node)
            endif
        enddo
! Let vloc=vloc+d first then communicate, and then it should same # of loop (avoiding loop at the whole domain)
!=============================
        call communicate_respml_ad_sub(vloc,ndf,nn,send_address,ad_length)
!----------------------------------------------------------------------------------------------
!vloc(:,:)=vloc(:,:)+d(:,:)
!----------------------------------------------------------------------------------------------
! Clear matrix local first interal then boundary processor nodes
        do icount=1,nn_local
            node=node_local(icount)
            av_tmp(1:ndf,node)=0.0d0
            if (seqcPML(node) > 0) then
                av_tmp(1:ndf,nn+seqcPML(node))=0.0d0
            endif
        enddo
        do icount=1,nn_local_com
            node=global_com(local_com(icount))
            av_tmp(1:ndf,node)=0.0d0
            if (seqcPML(node) > 0) then
                av_tmp(1:ndf,nn+seqcPML(node))=0.0d0
            endif
        enddo
!		av_tmp(:,:) = 0.0d0
!======================================
!!!!!! av_tmp is a temporary residual
        call blockgmresnew(x,vloc,dold,av_tmp,hg,ien,fext,ne_local,ien_local,node_local,nn_local,&
                           fden,fvis,I_fluid,rngface)
        call communicate_respml_ad_sub(av_tmp,ndf,nn,send_address,ad_length)
!===================
! avloc(:)=0.0d0
!==================
        call equalpml_pa(av_tmp,avloc,ndf,nn,node_local,nn_local)
        do icount=1, nn_local
            node=node_local(icount)
            do jcount=1,ndf
                kcount=(node-1)*ndf+jcount
                avloc(kcount) = (-avloc(kcount)+bg(kcount))/eps ! get Av,bg=-r(u)
            enddo
            if (seqcPML(node) > 0) then
                do jcount=1,ndf
                    kcount=(nn+seqcPML(node)-1)*ndf+jcount
                    avloc(kcount) = (-avloc(kcount)+bg(kcount))/eps ! get Av,bg=-r(u)
                enddo
            endif
        enddo
        call setidpml_pa(avloc,ndf,nn,nn_PML,id,node_local,nn_local)
        space1(:)=0.0d0
        space2(:)=0.0d0
        end_time=mpi_wtime()
        do i=1,j
            do icount=1,nn_local
                node=node_local(icount)
                do jcount=1,ndf
                    kcount=(node-1)*ndf+jcount
                    space1(i)=space1(i)+avloc(kcount)*Vm((icount-1)*ndf+jcount,i)
                enddo
                if (seqcPML(node) > 0) then
                    do jcount=1,ndf
                        kcount=(nn+seqcPML(node)-1)*ndf+jcount
                        space1(i)=space1(i)+avloc(kcount)*Vm((nn_local+seqcPMLlocal(icount)-1)*ndf+jcount,i)
                    enddo
                endif
            enddo
! Do the vector product of v_i * v_i set it to be h_i,j
        enddo  ! construct AVj and hi,j
        call mpi_barrier(mpi_comm_world,ierror)
        call mpi_allreduce(space1(1),space2(1),j,mpi_double_precision,mpi_sum,mpi_comm_world,ierror)
        end_time=mpi_wtime()-end_time
        Hm(1:j,j)=space2(1:j)
        do icount = 1, nn_local
            node=node_local(icount)
            do jcount=1,ndf
                kcount=(icount-1)*ndf+jcount
                do i=1,j
                    Vm(kcount,j+1) = Vm(kcount,j+1)-Hm(i,j)*Vm(kcount,i)
                enddo
                Vm(kcount,j+1)=Vm(kcount,j+1)+avloc((node-1)*ndf+jcount)
            enddo
            if (seqcPML(node) > 0) then
                do jcount=1,ndf
                    kcount=(nn_local+seqcPMLlocal(icount)-1)*ndf+jcount
                    do i=1,j
                        Vm(kcount,j+1) = Vm(kcount,j+1)-Hm(i,j)*Vm(kcount,i)
                    enddo
                    Vm(kcount,j+1)=Vm(kcount,j+1)+avloc((nn+seqcPML(node)-1)*ndf+jcount)
                enddo
            endif
        enddo  ! construct v(j+1)
        temp=0.0d0
        do icount = 1, nn_local
            node=node_local(icount)
            do jcount=1,ndf
                kcount=(icount-1)*ndf+jcount
                temp=temp+Vm(kcount,j+1)*Vm(kcount,j+1)
            enddo
            if (seqcPML(node) > 0) then
                do jcount=1,ndf
                    kcount=(nn_local+seqcPMLlocal(icount)-1)*ndf+jcount
                    temp=temp+Vm(kcount,j+1)*Vm(kcount,j+1)
                enddo
            endif
        enddo
        call mpi_barrier(mpi_comm_world,ierror)
        call mpi_allreduce(temp,rnorm0,1,mpi_double_precision,mpi_sum,mpi_comm_world,ierror)

        Hm(j+1,j) = sqrt(rnorm0)
        do icount = 1, nn_local
            node=node_local(icount)
            do jcount=1,ndf
                kcount=(icount-1)*ndf+jcount
                Vm(kcount,j+1)=Vm(kcount,j+1)/Hm(j+1,j)
            enddo
            if (seqcPML(node) > 0) then
                do jcount=1,ndf
                    kcount=(nn_local+seqcPMLlocal(icount)-1)*ndf+jcount
                    Vm(kcount,j+1)=Vm(kcount,j+1)/Hm(j+1,j)
                enddo
            endif
        enddo
    enddo  ! end inner loop
!====================================================================================
! Use lapack to solve the LS problem
!	call DGELS(TRAN,inner+1,inner,1,Hm,inner+1,beta,inner+1,workls,2*inner,INFO)
!===================================================================================
! Use Givens rotation to solve the LS problem
    call mpi_barrier(mpi_comm_world,ierror)
    call givens(Hm,inner,beta)
!!!!!!!!!!!!beta(1:inner) is ym, the solution!!!!!!!!!!!!!!!!!!!!!!!!!
    Vy(:) = 0
    do icount=1, nn_local
        node=node_local(icount)
        do jcount=1,ndf
            kcount=(node-1)*ndf+jcount
            do i=1,inner
                Vy(kcount)=Vy(kcount)+Vm((icount-1)*ndf+jcount,i)*beta(i)
            enddo
            x0(kcount)=x0(kcount)+Vy(kcount)
        enddo
        if (seqcPML(node) > 0) then
            do jcount=1,ndf
                kcount=(nn+seqcPML(node)-1)*ndf+jcount
                do i=1,inner
                    Vy(kcount)=Vy(kcount)+Vm((nn_local+seqcPMLlocal(icount)-1)*ndf+jcount,i)*beta(i)
                enddo
                x0(kcount)=x0(kcount)+Vy(kcount)
            enddo
        endif
    enddo ! calculate Xm
    do icount = 1, nn_local
        node=node_local(icount)
        do jcount=1,ndf
            kcount=(node-1)*ndf+jcount
            dv(kcount) = eps/w(kcount)*x0(kcount)
        enddo
        if (seqcPML(node) > 0) then
            do jcount=1,ndf
                kcount=(nn+seqcPML(node)-1)*ndf+jcount
                dv(kcount) = eps/w(kcount)*x0(kcount)
            enddo
        endif
    enddo
!============================
! Clear matirx locally
    do icount=1,nn_local
        node=node_local(icount)
        vloc(1:ndf,node)=0.0d0
        if (seqcPML(node) > 0) then
            vloc(1:ndf,nn+seqcPML(node))=0.0d0
        endif
    enddo
! clear all the processor internal nodes
    do icount=1,nn_local_com
        node=global_com(local_com(icount))
        vloc(1:ndf,node)=0.0d0
        if (seqcPML(node) > 0) then
            vloc(1:ndf,nn+seqcPML(node))=0.0d0
        endif
    enddo
! clear the processor boundary nodes
!	vloc(:,:) = 0
!==============================
    call equalpml_pa(dv,vloc,ndf,nn,node_local,nn_local)
    do icount=1,nn_local
        node=node_local(icount)
        vloc(1:ndf,node)=vloc(1:ndf,node)+d(1:ndf,node)
        if (seqcPML(node) > 0) then
            vloc(1:ndf,nn+seqcPML(node))=vloc(1:ndf,nn+seqcPML(node))+qv(1:ndf,node)
        endif
    enddo
    call communicate_respml_ad_sub(vloc,ndf,nn,send_address,ad_length)
!==============================
!	vloc(:,:) = vloc(:,:)+d(:,:)
!===============================
! Clear matrix local first interal then boundary processor nodes
    do icount=1,nn_local
        node=node_local(icount)
        av_tmp(1:ndf,node)=0.0d0
        if (seqcPML(node) > 0) then
            av_tmp(1:ndf,nn+seqcPML(node))=0.0d0
        endif
    enddo
    do icount=1,nn_local_com
        node=global_com(local_com(icount))
        av_tmp(1:ndf,node)=0.0d0
        if (seqcPML(node) > 0) then
            av_tmp(1:ndf,nn+seqcPML(node))=0.0d0
        endif
    enddo
!	av_tmp(:,:) = 0.0d0
!================================
    call blockgmresnew(x,vloc,dold,av_tmp,hg,ien,fext,ne_local,ien_local,node_local,nn_local,&
                       fden,fvis,I_fluid,rngface)
    call communicate_respml_ad_sub(av_tmp,ndf,nn,send_address,ad_length)
!==================
!avloc(:)=0.0d0
!==================
    call equalpml_pa(av_tmp,avloc,ndf,nn,node_local,nn_local)
    call setidpml_pa(avloc,ndf,nn,nn_PML,id,node_local,nn_local)
    do icount=1, nn_local
        node=node_local(icount)
        do jcount=1,ndf
            kcount=(node-1)*ndf+jcount
            avloc(kcount) = (-avloc(kcount)+bg(kcount))/eps
        enddo
        if (seqcPML(node) > 0) then
            do jcount=1,ndf
                kcount=(nn+seqcPML(node)-1)*ndf+jcount
                avloc(kcount) = (-avloc(kcount)+bg(kcount))/eps
            enddo
        endif
    enddo
!!!!!!!!!!calculate AXm
    do icount=1,nn_local
        node=node_local(icount)
        do jcount=1,ndf
            kcount=(node-1)*ndf+jcount
            r0(kcount) = bg(kcount)-avloc(kcount)
        enddo
        if (seqcPML(node) > 0) then
            do jcount=1,ndf
                kcount=(nn+seqcPML(node)-1)*ndf+jcount
                r0(kcount) = bg(kcount)-avloc(kcount)
            enddo
        endif
    enddo !update r0=f-AX0

    call getnormpml_pa(r0,ndf,nn,node_local,nn_local,rnorm0)
    err = sqrt(rnorm0)
    rnorm = sqrt(rnorm0)
    iouter = iouter + 1
    if (myid == 0) write(*,*) '||r0||=||f-AX0||=',err

111 continue  ! end outer loop

    dg_sent(1:ndf*(nn+nn_PML))=0.0d0
    do icount=1, nn_local
        node=node_local(icount)
        do jcount=1,ndf
            kcount=(node-1)*ndf+jcount
            dg_sent(kcount) = 1/w(kcount)*x0(kcount) ! unscaled x0
        enddo
        if (seqcPML(node) > 0) then
            do jcount=1,ndf
                kcount=(nn+seqcPML(node)-1)*ndf+jcount
                dg_sent(kcount) = 1/w(kcount)*x0(kcount) ! unscaled x0
            enddo
        endif
    enddo
    dg(1:ndf*(nn+nn_PML))=0.0d0
    call mpi_barrier(mpi_comm_world,ierror)
    call mpi_allreduce(dg_sent(1),dg(1),ndf*(nn+nn_PML),mpi_double_precision,mpi_sum,mpi_comm_world,ierror)
!	call mpi_bcast(dg(1),ndf*nn,mpi_double_precision,0,mpi_comm_world,ierror)
    call mpi_barrier(mpi_comm_world,ierror)

return
end subroutine gmres

max_stars_repo_name: biofluids/IFEM-archive

max_issues_repo_name: biofluids/IFEM-archive

max_forks_repo_name: biofluids/IFEM-archive


text:
c
c  ...  This is a slightly modified version of the gauss
c  ...  quadrature routine from netlib.  All references to
c  ...  d1mach() have been changed to d1maci(), to avoid
c  ...  potential conflicts with dvode
c
*DECK D9LGMC
      DOUBLE PRECISION FUNCTION D9LGMC (X)
C***BEGIN PROLOGUE  D9LGMC
C***SUBSIDIARY
C***PURPOSE  Compute the log Gamma correction factor so that
C            LOG(DGAMMA(X)) = LOG(SQRT(2*PI)) + (X-5.)*LOG(X) - X
C            + D9LGMC(X).
C***LIBRARY   SLATEC (FNLIB)
C***CATEGORY  C7E
C***TYPE      DOUBLE PRECISION (R9LGMC-S, D9LGMC-D, C9LGMC-C)
C***KEYWORDS  COMPLETE GAMMA FUNCTION, CORRECTION TERM, FNLIB,
C             LOG GAMMA, LOGARITHM, SPECIAL FUNCTIONS
C***AUTHOR  Fullerton, W., (LANL)
C***DESCRIPTION
C
C Compute the log gamma correction factor for X .GE. 10. so that
C LOG (DGAMMA(X)) = LOG(SQRT(2*PI)) + (X-.5)*LOG(X) - X + D9lGMC(X)
C
C Series for ALGM       on the interval  0.          to  1.00000E-02
C                                        with weighted error   1.28E-31
C                                         log weighted error  30.89
C                               significant figures required  29.81
C                                    decimal places required  31.48
C
C***REFERENCES  (NONE)
C***ROUTINES CALLED  D1MACI, DCSEVL, INITDS, XERMSG
C***REVISION HISTORY  (YYMMDD)
C   770601  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   890531  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
C   900720  Routine changed from user-callable to subsidiary.  (WRB)
C***END PROLOGUE  D9LGMC
      DOUBLE PRECISION X, ALGMCS(15), XBIG, XMAX, DCSEVL, D1MACI
      LOGICAL FIRST
      SAVE ALGMCS, NALGM, XBIG, XMAX, FIRST
      DATA ALGMCS(  1) / +.1666389480 4518632472 0572965082 2 D+0      /
      DATA ALGMCS(  2) / -.1384948176 0675638407 3298605913 5 D-4      /
      DATA ALGMCS(  3) / +.9810825646 9247294261 5717154748 7 D-8      /
      DATA ALGMCS(  4) / -.1809129475 5724941942 6330626671 9 D-10     /
      DATA ALGMCS(  5) / +.6221098041 8926052271 2601554341 6 D-13     /
      DATA ALGMCS(  6) / -.3399615005 4177219443 0333059966 6 D-15     /
      DATA ALGMCS(  7) / +.2683181998 4826987489 5753884666 6 D-17     /
      DATA ALGMCS(  8) / -.2868042435 3346432841 4462239999 9 D-19     /
      DATA ALGMCS(  9) / +.3962837061 0464348036 7930666666 6 D-21     /
      DATA ALGMCS( 10) / -.6831888753 9857668701 1199999999 9 D-23     /
      DATA ALGMCS( 11) / +.1429227355 9424981475 7333333333 3 D-24     /
      DATA ALGMCS( 12) / -.3547598158 1010705471 9999999999 9 D-26     /
      DATA ALGMCS( 13) / +.1025680058 0104709120 0000000000 0 D-27     /
      DATA ALGMCS( 14) / -.3401102254 3167487999 9999999999 9 D-29     /
      DATA ALGMCS( 15) / +.1276642195 6300629333 3333333333 3 D-30     /
      DATA FIRST /.TRUE./
C***FIRST EXECUTABLE STATEMENT  D9LGMC
      IF (FIRST) THEN
         NALGM = INITDS (ALGMCS, 15, REAL(D1MACI(3)) )
         XBIG = 1.0D0/SQRT(D1MACI(3))
         XMAX = EXP (MIN(LOG(D1MACI(2)/12.D0), -LOG(12.D0*D1MACI(1))))
      ENDIF
      FIRST = .FALSE.
C
      IF (X .LT. 10.D0) CALL XERMSG ('SLATEC', 'D9LGMC',
     +   'X MUST BE GE 10', 1, 2)
      IF (X.GE.XMAX) GO TO 20
C
      D9LGMC = 1.D0/(12.D0*X)
      IF (X.LT.XBIG) D9LGMC = DCSEVL (2.0D0*(10.D0/X)**2-1.D0, ALGMCS,
     1  NALGM) / X
      RETURN
C
 20   D9LGMC = 0.D0
      CALL XERMSG ('SLATEC', 'D9LGMC', 'X SO BIG D9LGMC UNDERFLOWS', 2,
     +   1)
      RETURN
C
      END
*DECK DCSEVL
      DOUBLE PRECISION FUNCTION DCSEVL (X, CS, N)
C***BEGIN PROLOGUE  DCSEVL
C***PURPOSE  Evaluate a Chebyshev series.
C***LIBRARY   SLATEC (FNLIB)
C***CATEGORY  C3A2
C***TYPE      DOUBLE PRECISION (CSEVL-S, DCSEVL-D)
C***KEYWORDS  CHEBYSHEV SERIES, FNLIB, SPECIAL FUNCTIONS
C***AUTHOR  Fullerton, W., (LANL)
C***DESCRIPTION
C
C  Evaluate the N-term Chebyshev series CS at X.  Adapted from
C  a method presented in the paper by Broucke referenced below.
C
C       Input Arguments --
C  X    value at which the series is to be evaluated.
C  CS   array of N terms of a Chebyshev series.  In evaluating
C       CS, only half the first coefficient is summed.
C  N    number of terms in array CS.
C
C***REFERENCES  R. Broucke, Ten subroutines for the manipulation of
C                 Chebyshev series, Algorithm 446, Communications of
C                 the A.C.M. 16, (1973) pp. 254-256.
C               L. Fox and I. B. Parker, Chebyshev Polynomials in
C                 Numerical Analysis, Oxford University Press, 1968,
C                 page 56.
C***ROUTINES CALLED  D1MACI, XERMSG
C***REVISION HISTORY  (YYMMDD)
C   770401  DATE WRITTEN
C   890831  Modified array declarations.  (WRB)
C   890831  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
C   900329  Prologued revised extensively and code rewritten to allow
C           X to be slightly outside interval (-1,+1).  (WRB)
C   920501  Reformatted the REFERENCES section.  (WRB)
C***END PROLOGUE  DCSEVL
      DOUBLE PRECISION B0, B1, B2, CS(*), ONEPL, TWOX, X, D1MACI
      LOGICAL FIRST
      SAVE FIRST, ONEPL
      DATA FIRST /.TRUE./
C***FIRST EXECUTABLE STATEMENT  DCSEVL
      IF (FIRST) ONEPL = 1.0D0 + D1MACI(4)
      FIRST = .FALSE.
      IF (N .LT. 1) CALL XERMSG ('SLATEC', 'DCSEVL',
     +   'NUMBER OF TERMS .LE. 0', 2, 2)
      IF (N .GT. 1000) CALL XERMSG ('SLATEC', 'DCSEVL',
     +   'NUMBER OF TERMS .GT. 1000', 3, 2)
      IF (ABS(X) .GT. ONEPL) CALL XERMSG ('SLATEC', 'DCSEVL',
     +   'X OUTSIDE THE INTERVAL (-1,+1)', 1, 1)
C
      B1 = 0.0D0
      B0 = 0.0D0
      TWOX = 2.0D0*X
      DO 10 I = 1,N
         B2 = B1
         B1 = B0
         NI = N + 1 - I
         B0 = TWOX*B1 - B2 + CS(NI)
   10 CONTINUE
C
      DCSEVL = 0.5D0*(B0-B2)
C
      RETURN
      END
*DECK DGAMLM
      SUBROUTINE DGAMLM (XMIN, XMAX)
C***BEGIN PROLOGUE  DGAMLM
C***PURPOSE  Compute the minimum and maximum bounds for the argument in
C            the Gamma function.
C***LIBRARY   SLATEC (FNLIB)
C***CATEGORY  C7A, R2
C***TYPE      DOUBLE PRECISION (GAMLIM-S, DGAMLM-D)
C***KEYWORDS  COMPLETE GAMMA FUNCTION, FNLIB, LIMITS, SPECIAL FUNCTIONS
C***AUTHOR  Fullerton, W., (LANL)
C***DESCRIPTION
C
C Calculate the minimum and maximum legal bounds for X in gamma(X).
C XMIN and XMAX are not the only bounds, but they are the only non-
C trivial ones to calculate.
C
C             Output Arguments --
C XMIN   double precision minimum legal value of X in gamma(X).  Any
C        smaller value of X might result in underflow.
C XMAX   double precision maximum legal value of X in gamma(X).  Any
C        larger value of X might cause overflow.
C
C***REFERENCES  (NONE)
C***ROUTINES CALLED  D1MACI, XERMSG
C***REVISION HISTORY  (YYMMDD)
C   770601  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   890531  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
C***END PROLOGUE  DGAMLM
      DOUBLE PRECISION XMIN, XMAX, ALNBIG, ALNSML, XLN, XOLD, D1MACI
C***FIRST EXECUTABLE STATEMENT  DGAMLM
      ALNSML = LOG(D1MACI(1))
      XMIN = -ALNSML
      DO 10 I=1,10
        XOLD = XMIN
        XLN = LOG(XMIN)
        XMIN = XMIN - XMIN*((XMIN+0.5D0)*XLN - XMIN - 0.2258D0 + ALNSML)
     1    / (XMIN*XLN+0.5D0)
        IF (ABS(XMIN-XOLD).LT.0.005D0) GO TO 20
 10   CONTINUE
      CALL XERMSG ('SLATEC', 'DGAMLM', 'UNABLE TO FIND XMIN', 1, 2)
C
 20   XMIN = -XMIN + 0.01D0
C
      ALNBIG = LOG (D1MACI(2))
      XMAX = ALNBIG
      DO 30 I=1,10
        XOLD = XMAX
        XLN = LOG(XMAX)
        XMAX = XMAX - XMAX*((XMAX-0.5D0)*XLN - XMAX + 0.9189D0 - ALNBIG)
     1    / (XMAX*XLN-0.5D0)
        IF (ABS(XMAX-XOLD).LT.0.005D0) GO TO 40
 30   CONTINUE
      CALL XERMSG ('SLATEC', 'DGAMLM', 'UNABLE TO FIND XMAX', 2, 2)
C
 40   XMAX = XMAX - 0.01D0
      XMIN = MAX (XMIN, -XMAX+1.D0)
C
      RETURN
      END
*DECK DGAMMA
      DOUBLE PRECISION FUNCTION DGAMMA (X)
C***BEGIN PROLOGUE  DGAMMA
C***PURPOSE  Compute the complete Gamma function.
C***LIBRARY   SLATEC (FNLIB)
C***CATEGORY  C7A
C***TYPE      DOUBLE PRECISION (GAMMA-S, DGAMMA-D, CGAMMA-C)
C***KEYWORDS  COMPLETE GAMMA FUNCTION, FNLIB, SPECIAL FUNCTIONS
C***AUTHOR  Fullerton, W., (LANL)
C***DESCRIPTION
C
C DGAMMA(X) calculates the double precision complete Gamma function
C for double precision argument X.
C
C Series for GAM        on the interval  0.          to  1.00000E+00
C                                        with weighted error   5.79E-32
C                                         log weighted error  31.24
C                               significant figures required  30.00
C                                    decimal places required  32.05
C
C***REFERENCES  (NONE)
C***ROUTINES CALLED  D1MACI, D9LGMC, DCSEVL, DGAMLM, INITDS, XERMSG
C***REVISION HISTORY  (YYMMDD)
C   770601  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   890911  Removed unnecessary intrinsics.  (WRB)
C   890911  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
C   920618  Removed space from variable name.  (RWC, WRB)
C***END PROLOGUE  DGAMMA
      DOUBLE PRECISION X, GAMCS(42), DXREL, PI, SINPIY, SQ2PIL, XMAX,
     1  XMIN, Y, D9LGMC, DCSEVL, D1MACI
      LOGICAL FIRST
C
      SAVE GAMCS, PI, SQ2PIL, NGAM, XMIN, XMAX, DXREL, FIRST
      DATA GAMCS(  1) / +.8571195590 9893314219 2006239994 2 D-2      /
      DATA GAMCS(  2) / +.4415381324 8410067571 9131577165 2 D-2      /
      DATA GAMCS(  3) / +.5685043681 5993633786 3266458878 9 D-1      /
      DATA GAMCS(  4) / -.4219835396 4185605010 1250018662 4 D-2      /
      DATA GAMCS(  5) / +.1326808181 2124602205 8400679635 2 D-2      /
      DATA GAMCS(  6) / -.1893024529 7988804325 2394702388 6 D-3      /
      DATA GAMCS(  7) / +.3606925327 4412452565 7808221722 5 D-4      /
      DATA GAMCS(  8) / -.6056761904 4608642184 8554829036 5 D-5      /
      DATA GAMCS(  9) / +.1055829546 3022833447 3182350909 3 D-5      /
      DATA GAMCS( 10) / -.1811967365 5423840482 9185589116 6 D-6      /
      DATA GAMCS( 11) / +.3117724964 7153222777 9025459316 9 D-7      /
      DATA GAMCS( 12) / -.5354219639 0196871408 7408102434 7 D-8      /
      DATA GAMCS( 13) / +.9193275519 8595889468 8778682594 0 D-9      /
      DATA GAMCS( 14) / -.1577941280 2883397617 6742327395 3 D-9      /
      DATA GAMCS( 15) / +.2707980622 9349545432 6654043308 9 D-10     /
      DATA GAMCS( 16) / -.4646818653 8257301440 8166105893 3 D-11     /
      DATA GAMCS( 17) / +.7973350192 0074196564 6076717535 9 D-12     /
      DATA GAMCS( 18) / -.1368078209 8309160257 9949917230 9 D-12     /
      DATA GAMCS( 19) / +.2347319486 5638006572 3347177168 8 D-13     /
      DATA GAMCS( 20) / -.4027432614 9490669327 6657053469 9 D-14     /
      DATA GAMCS( 21) / +.6910051747 3721009121 3833697525 7 D-15     /
      DATA GAMCS( 22) / -.1185584500 2219929070 5238712619 2 D-15     /
      DATA GAMCS( 23) / +.2034148542 4963739552 0102605193 2 D-16     /
      DATA GAMCS( 24) / -.3490054341 7174058492 7401294910 8 D-17     /
      DATA GAMCS( 25) / +.5987993856 4853055671 3505106602 6 D-18     /
      DATA GAMCS( 26) / -.1027378057 8722280744 9006977843 1 D-18     /
      DATA GAMCS( 27) / +.1762702816 0605298249 4275966074 8 D-19     /
      DATA GAMCS( 28) / -.3024320653 7353062609 5877211204 2 D-20     /
      DATA GAMCS( 29) / +.5188914660 2183978397 1783355050 6 D-21     /
      DATA GAMCS( 30) / -.8902770842 4565766924 4925160106 6 D-22     /
      DATA GAMCS( 31) / +.1527474068 4933426022 7459689130 6 D-22     /
      DATA GAMCS( 32) / -.2620731256 1873629002 5732833279 9 D-23     /
      DATA GAMCS( 33) / +.4496464047 8305386703 3104657066 6 D-24     /
      DATA GAMCS( 34) / -.7714712731 3368779117 0390152533 3 D-25     /
      DATA GAMCS( 35) / +.1323635453 1260440364 8657271466 6 D-25     /
      DATA GAMCS( 36) / -.2270999412 9429288167 0231381333 3 D-26     /
      DATA GAMCS( 37) / +.3896418998 0039914493 2081663999 9 D-27     /
      DATA GAMCS( 38) / -.6685198115 1259533277 9212799999 9 D-28     /
      DATA GAMCS( 39) / +.1146998663 1400243843 4761386666 6 D-28     /
      DATA GAMCS( 40) / -.1967938586 3451346772 9510399999 9 D-29     /
      DATA GAMCS( 41) / +.3376448816 5853380903 3489066666 6 D-30     /
      DATA GAMCS( 42) / -.5793070335 7821357846 2549333333 3 D-31     /
      DATA PI / 3.1415926535 8979323846 2643383279 50 D0 /
      DATA SQ2PIL / 0.9189385332 0467274178 0329736405 62 D0 /
      DATA FIRST /.TRUE./
C***FIRST EXECUTABLE STATEMENT  DGAMMA
      IF (FIRST) THEN
         NGAM = INITDS (GAMCS, 42, 0.1*REAL(D1MACI(3)) )
C
         CALL DGAMLM (XMIN, XMAX)
         DXREL = SQRT(D1MACI(4))
      ENDIF
      FIRST = .FALSE.
C
      Y = ABS(X)
      IF (Y.GT.10.D0) GO TO 50
C
C COMPUTE GAMMA(X) FOR -XBND .LE. X .LE. XBND.  REDUCE INTERVAL AND FIND
C GAMMA(1+Y) FOR 0.0 .LE. Y .LT. 1.0 FIRST OF ALL.
C
      N = X
      IF (X.LT.0.D0) N = N - 1
      Y = X - N
      N = N - 1
      DGAMMA = 0.9375D0 + DCSEVL (2.D0*Y-1.D0, GAMCS, NGAM)
      IF (N.EQ.0) RETURN
C
      IF (N.GT.0) GO TO 30
C
C COMPUTE GAMMA(X) FOR X .LT. 1.0
C
      N = -N
      IF (X .EQ. 0.D0) CALL XERMSG ('SLATEC', 'DGAMMA', 'X IS 0', 4, 2)
      IF (X .LT. 0.0 .AND. X+N-2 .EQ. 0.D0) CALL XERMSG ('SLATEC',
     +   'DGAMMA', 'X IS A NEGATIVE INTEGER', 4, 2)
      IF (X .LT. (-0.5D0) .AND. ABS((X-AINT(X-0.5D0))/X) .LT. DXREL)
     +   CALL XERMSG ('SLATEC', 'DGAMMA',
     +   'ANSWER LT HALF PRECISION BECAUSE X TOO NEAR NEGATIVE INTEGER',
     +   1, 1)
C
      DO 20 I=1,N
        DGAMMA = DGAMMA/(X+I-1 )
 20   CONTINUE
      RETURN
C
C GAMMA(X) FOR X .GE. 2.0 AND X .LE. 10.0
C
 30   DO 40 I=1,N
        DGAMMA = (Y+I) * DGAMMA
 40   CONTINUE
      RETURN
C
C GAMMA(X) FOR ABS(X) .GT. 10.0.  RECALL Y = ABS(X).
C
 50   IF (X .GT. XMAX) CALL XERMSG ('SLATEC', 'DGAMMA',
     +   'X SO BIG GAMMA OVERFLOWS', 3, 2)
C
      DGAMMA = 0.D0
      IF (X .LT. XMIN) CALL XERMSG ('SLATEC', 'DGAMMA',
     +   'X SO SMALL GAMMA UNDERFLOWS', 2, 1)
      IF (X.LT.XMIN) RETURN
C
      DGAMMA = EXP ((Y-0.5D0)*LOG(Y) - Y + SQ2PIL + D9LGMC(Y) )
      IF (X.GT.0.D0) RETURN
C
      IF (ABS((X-AINT(X-0.5D0))/X) .LT. DXREL) CALL XERMSG ('SLATEC',
     +   'DGAMMA',
     +   'ANSWER LT HALF PRECISION, X TOO NEAR NEGATIVE INTEGER', 1, 1)
C
      SINPIY = SIN (PI*Y)
      IF (SINPIY .EQ. 0.D0) CALL XERMSG ('SLATEC', 'DGAMMA',
     +   'X IS A NEGATIVE INTEGER', 4, 2)
C
      DGAMMA = -PI/(Y*SINPIY*DGAMMA)
C
      RETURN
      END
*DECK INITDS
      FUNCTION INITDS (OS, NOS, ETA)
C***BEGIN PROLOGUE  INITDS
C***PURPOSE  Determine the number of terms needed in an orthogonal
C            polynomial series so that it meets a specified accuracy.
C***LIBRARY   SLATEC (FNLIB)
C***CATEGORY  C3A2
C***TYPE      DOUBLE PRECISION (INITS-S, INITDS-D)
C***KEYWORDS  CHEBYSHEV, FNLIB, INITIALIZE, ORTHOGONAL POLYNOMIAL,
C             ORTHOGONAL SERIES, SPECIAL FUNCTIONS
C***AUTHOR  Fullerton, W., (LANL)
C***DESCRIPTION
C
C  Initialize the orthogonal series, represented by the array OS, so
C  that INITDS is the number of terms needed to insure the error is no
C  larger than ETA.  Ordinarily, ETA will be chosen to be one-tenth
C  machine precision.
C
C             Input Arguments --
C   OS     double precision array of NOS coefficients in an orthogonal
C          series.
C   NOS    number of coefficients in OS.
C   ETA    single precision scalar containing requested accuracy of
C          series.
C
C***REFERENCES  (NONE)
C***ROUTINES CALLED  XERMSG
C***REVISION HISTORY  (YYMMDD)
C   770601  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   890831  Modified array declarations.  (WRB)
C   891115  Modified error message.  (WRB)
C   891115  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
C***END PROLOGUE  INITDS
      DOUBLE PRECISION OS(*)
C***FIRST EXECUTABLE STATEMENT  INITDS
      IF (NOS .LT. 1) CALL XERMSG ('SLATEC', 'INITDS',
     +   'Number of coefficients is less than 1', 2, 1)
C
      ERR = 0.
      DO 10 II = 1,NOS
        I = NOS + 1 - II
        ERR = ERR + ABS(REAL(OS(I)))
        IF (ERR.GT.ETA) GO TO 20
   10 CONTINUE
C
   20 IF (I .EQ. NOS) CALL XERMSG ('SLATEC', 'INITDS',
     +   'Chebyshev series too short for specified accuracy', 1, 1)
      INITDS = I
C
      RETURN
      END