text:
      subroutine clgen(j)

!!    ~ ~ ~ PURPOSE ~ ~ ~
!!    this subroutine calculates the daylength, distribution of
!!    radiation throughout the day and maximum radiation for day

!!    ~ ~ ~ INCOMING VARIABLES ~ ~ ~
!!    name        |units         |definition
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!!    iida        |julian date   |day being simulated (current julian date)
!!    ipet        |none          |code for potential ET method
!!                               |0 Priestley-Taylor method
!!                               |1 Penman/Monteith method
!!                               |2 Hargreaves method
!!    j           |none          |HRU number
!!    latcos(:)   |none          |Cos(Latitude) for HRU
!!    latsin(:)   |none          |Sin(Latitude) for HRU
!!    npcp(:)     |none          |prior day category
!!                               |1 dry day
!!                               |2 wet day
!!    subp(:)     |mm H2O        |precipitation for the day in HRU
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~

!!    ~ ~ ~ OUTGOING VARIABLES ~ ~ ~
!!    name        |units         |definition
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!!    dayl(:)     |hours         |day length
!!    frad(:,:)   |none          |fraction of solar radiation occuring during
!!                               |hour in day in HRU
!!    hru_rmx(:)  |MJ/m^2        |maximum possible radiation for the day in HRU
!!    npcp(:)     |none          |prior day category
!!                               |1 dry day
!!                               |2 wet day
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~

!!    ~ ~ ~ LOCAL DEFINITIONS ~ ~ ~
!!    name        |units         |definition
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!!    ch          |none          |(-Tan(sd)*Tan(lat))
!!    cosrho(:)   |none          |Cos(zenith angle for hour)
!!    dd          |none          |relative distance of the earth from the sun
!!    h           |none          |Acos(-Tan(sd)*Tan(lat))
!!    ii          |none          |counter
!!    sd          |radians       |solar declination: latitude at which the sun
!!                               |is directly overhead at noon
!!    totrho      |none          |sum of cosrho values for all hours of day
!!    yc          |none          |Cos(sd)*Cos(lat)
!!    ys          |none          |Sin(sd)*Sin(lat)
!!    w           |none          |hour angle
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~

!!    ~ ~ ~ SUBROUTINES/FUNCTIONS CALLED ~ ~ ~
!!    Intrinsic: Sin, Tan, Acos, Cos, Real

!!    ~ ~ ~ ~ ~ ~ END SPECIFICATIONS ~ ~ ~ ~ ~ ~

      use parm

      integer, intent (in) :: j
      integer :: ii
      real*8 :: sd, ch, h, ys, yc, dd, w, cosrho(nstep), totrho

      !! Reset prior day category for precipitation
      if (subp(j) >= 0.1) then
        npcp(j) = 2
      else
        npcp(j) = 1
      end if

      !! Calculate Daylength !!
      !! calculate solar declination: equation 2.1.2 in SWAT manual
      sd = 0.
      sd = Asin(.4 * Sin((dfloat(iida) - 82.) / 58.09))  !!365/2pi = 58.09

      !! calculate the relative distance of the earth from the sun
      !! the eccentricity of the orbit
      !! equation 2.1.1 in SWAT manual
      dd = 0.
      dd = 1.0 + 0.033 * Cos(dfloat(iida) / 58.09)

      !!daylength = 2 * Acos(-Tan(sd) * Tan(lat)) / omega
      !!where the angular velocity of the earth's rotation, omega, is equal
      !! to 15 deg/hr or 0.2618 rad/hr and 2/0.2618 = 7.6374
      !! equation 2.1.6 in SWAT manual

      ch = 0.
      h = 0.
      ch = -latsin(hru_sub(j)) * Tan(sd) / latcos(hru_sub(j))
      if (ch > 1.) then    !! ch will be >= 1. if latitude exceeds +/- 66.5 deg in winter
        h = 0.
      elseif (ch >= -1.) then
        h = Acos(ch)
      else
        h = 3.1416         !! latitude exceeds +/- 66.5 deg in summer
      endif
      dayl(j) = 7.6394 * h

      !! Calculate Potential (maximum) Radiation !!
      !! equation 2.2.7 in SWAT manual
      ys = 0.
      yc = 0.
      ys = latsin(hru_sub(j)) * Sin(sd)
      yc = latcos(hru_sub(j)) * Cos(sd)
      hru_rmx(j) = 30. * dd * (h * ys + yc * Sin(h))

      !! Calculate fraction of radiation recieved during each hour in day
      !! this calculation assumes solar noon (when the angle between the
      !! observer on the earth to the sun and a line normal to the earth's
      !! at that position is at a minimum) falls at 12 noon in day.
      !! equation 2.2.10 in SWAT manual

      cosrho = 0.
      totrho = 0.
      do ii = 1, nstep
        !!angular velocity times hour away from solar noon. To calculate
        !!radiation for an hour, the hour angle for the midpoint of the
        !!time period is used. time = 0. at solar noon with positive values
        !! in the morning and negative in the evening
        w = 0.
        w = (12.5 - dfloat(ii)) * 0.2618 * idt / 60.   !!0.2618 rad/hr
        cosrho(ii) = ys + yc * Cos(w)
        if (cosrho(ii) <= 0.) cosrho(ii) = 0.
        totrho = totrho + cosrho(ii)
      end do
      if (totrho > 0.001) then
        do ii = 1, nstep
          frad(j,ii) = cosrho(ii) / totrho
        end do
      end if

      return
      end

ext: f

max_stars_repo_name: GISWAT/erosion-sediment

max_issues_repo_name: GISWAT/erosion-sediment

max_forks_repo_name: GISWAT/erosion-sediment


text:
  subroutine skosgdr (photonEG, e, a0, sigkosgdr)

! ======================================================================
!
!   Giant dipole resonance cross section.
!
!    "Last" change: 13-AUG-2003 by NVMokhov
!    Edited by A. J. Sierk, LANL T-16, October, 2003.
!    Edited by AJS, LANL T-2, December, 2011.
!    Edited by LMK, XCP-3, July 2013 (included error protection).
!
! ======================================================================

    use, intrinsic:: iso_fortran_env, only: int32, real64

    implicit none
    class(PhotonEventGenerator), intent(inout) :: photonEG
    real(real64),   intent(in   ) :: e
    real(real64),   intent(in   ) :: a0
    real(real64),   intent(  out) :: sigkosgdr

    real(real64) :: aln, g1, g2, g4,g8, r1, r2, r4, r8, s2, t1, t2, &
         & t4, t8, temp, z

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

    z = log(e)
    aln = log(a0)
    temp = a0
    if (temp < divZerLim .and. temp > -divZerLim) then
       temp = divZerLim
       write(photonEG%io%message,1000) "112"
       call photonEG%io%print(4, 3, photonEG%io%message)
    end if
    s2 = one + (two/temp)**4
    r1 = (3.2d0 + 0.75d0*aln)/s2
    r2 = (four   + 0.125d0*aln)/s2
    r4 = 3.8d0  + 0.05d0*aln
    r8 = 3.65d0 - 0.05d0*aln
    t1 = (6.6d0 - 0.5d0*aln)/s2
    t2 = 3.4d0/s2
    t4 = 3.8d0  - 0.25d0*aln
    t8 = 3.5d0  - 0.16d0*aln

    temp = one + exp(thr*one*(t1 - z))
    g1 = exp(one*(r1 - z))/(temp)

    temp = one + exp(thr*two*(t2 - z))
    g2 = exp(two*(r2 - z))/(temp)

    temp = one + exp(thr*four*(t4 - z))
    g4 = exp(four*(r4 - z))/(temp)

    temp = one + exp(thr*two*four*(t8 - z))
    g8 = exp(two*four*(r8 - z))/(temp)

    sigkosgdr = g1 + g2 + g4 + g8
    return

! ======================================================================
1000 format("Divide by zero error prevented in ", &
          & "'photoEventGenerator.f90' line(s) ", A)
! ======================================================================
  end subroutine skosgdr

max_stars_repo_name: lanl/generalized-spallation-model

max_stars_repo_licenses: Intel Unlicense

max_issues_repo_name: lanl/generalized-spallation-model

max_issues_repo_licenses: Intel Unlicense

max_forks_repo_name: lanl/generalized-spallation-model

max_forks_repo_licenses: Intel Unlicense


text:
module integrals
  use log
  private

  ! Holder for integral data
  type ElectronicSystem
     integer,allocatable :: atm(:,:)
     integer,allocatable :: bas(:,:)
     integer :: natm, nelec, charge
     integer :: nbas, nprim, norb
     double precision, allocatable :: C_AO(:,:), E_AO(:)
     double precision, allocatable :: pop_AO(:)
     integer :: naos
  end type ElectronicSystem

  type UnrolledMOs
     character(len=2), allocatable :: atoms(:)
     double precision, allocatable :: Z(:)
     double precision, allocatable :: xyz(:,:)
     double precision, allocatable :: coeffs(:,:)
     double precision, allocatable :: exps(:)
     integer, allocatable :: icnt(:), ityp(:)

     integer :: natm, nprim, nmos
  end type UnrolledMOs

  integer :: current_n_atoms = 0

  ! libcint parameters
  integer,parameter :: ATM_SLOTS  = 6
  integer,parameter :: BAS_SLOTS  = 8
  integer,parameter :: AOS_SLOTS  = 5

  integer,parameter :: CHARGE_OF  = 1
  integer,parameter :: PTR_COORD  = 2
  integer,parameter :: NUC_MOD_OF = 3
  integer,parameter :: PTR_ZETA   = 4

  integer,parameter :: ATOM_OF    = 1
  integer,parameter :: ANG_OF     = 2
  integer,parameter :: NPRIM_OF   = 3
  integer,parameter :: NCTR_OF    = 4
  integer,parameter :: KAPPA_OF   = 5
  integer,parameter :: PTR_EXP    = 6
  integer,parameter :: PTR_COEFF  = 7

  ! AOS parameters
  integer,parameter :: PTR_ATOM = 1
  integer,parameter :: NAOS_OF = 2
  integer,parameter :: PTR_AOS_E = 3
  integer,parameter :: PTR_AOS_C = 4
  integer,parameter :: PTR_AOS_OCC = 5

  ! some libcint functions
  integer,external :: CINTcgto_spheric
  integer,external :: CINTcgto_cart
  double precision,external :: CINTgto_norm

  ! data from the python script
  include 'data-header.f90'

  ! spherical - cartesian conversions
  double precision   :: s2cd(6,5),s2cf(10,7),s2cg(15,9),s2ch(21,11)

  public :: integrals_initialize, integrals_init_from_file, integrals_build_basis
  public :: integrals_build_atomic_orbitals, integrals_overlaps
  public :: integrals_symm_ao2bas, integrals_symm_bas2ao
  public :: integrals_MO_AO_transform
  public :: integrals_write_to_wfn, integrals_unroll
  public :: ElectronicSystem, UnrolledMOs

contains

  subroutine integrals_initialize
    implicit none
    include 'data-content.f90'
    call set_s2c_coef
  end subroutine integrals_initialize

  subroutine integrals_init_from_file(unit, system, charge)
    use constants
    implicit none
    ! IO
    integer,intent(in)  :: unit
    integer,intent(in)  :: charge
    type(ElectronicSystem), intent(out) :: system

    ! Workspace
    character(len=2)  :: atomchar
    integer, allocatable :: atm(:,:)
    integer :: off, i, j
    integer :: info

    read(2,*, iostat=info) system%natm
    if (info .ne. 0) then
       call log_err('integrals_init_from_file', 'could not read atom number from file')
       error stop -1
    end if

    if (system%natm + current_n_atoms > MAX_NATOMS) then
       ! TODO : logger
       write(*,*) 'integrals_init_from_file:'
       print '(a,i8,/)', 'program has been compiled with MAX_NATOMS = ', MAX_NATOMS
       print '(a,i8,a,/)', 'your integral system now contains ', system%natm, ' atoms'
       print '(a/)', 'to run such a large system, recompile with a modified MAX_NATOMS parameter.'
       error stop -1
    end if

    read(2,*)                    !skip comment

    system%nelec = -charge
    system%charge = charge

    ! Set the pointer for the geometry
    off = ENV_COORD_PTR + 3*current_n_atoms

    allocate(atm(ATM_SLOTS, system%natm))
    do i = 1, system%natm
       atm(PTR_COORD, i) = off
       read(2,*, iostat=info) atomchar, ENV(off + 1 : off + 3)
       if (info .ne. 0) then
          call log_err('integrals_init_from_file', 'could not read atom or coordinate')
          error stop -1
       end if

       ENV(off+1) = ENV(off+1) / conv_bohr
       ENV(off+2) = ENV(off+2) / conv_bohr
       ENV(off+3) = ENV(off+3) / conv_bohr

       off = off + 3

       do j = 1, 110
          if (atomchar .eq. atomname(j)) then
             atm(CHARGE_OF, i) = j
             system%nelec = system%nelec + j
             exit
          end if
       end do

       if (j .eq. 110) then
          call log_err('integrals_init_from_file', 'atomic name is invalid =' // atomchar)
          error stop -1
       end if
    end do

    system%atm = atm

    ! Increment the overall atom counter
    current_n_atoms = current_n_atoms + system%natm

    ! Write out
    if (verbosity .ge. 0) then
       write(*,'(a,i4)') '          Charge:         ', system%charge
       write(*,'(a,i4)') '          N. of atoms:    ', system%natm
       write(*,'(a,i4)') '          N. of electrons:', system%nelec
    end if

    if (verbosity > 0) then
       write(*,'(a,i4)') '          N. atoms so far in evaluator:', current_n_atoms
    end if
  end subroutine integrals_init_from_file

  subroutine integrals_build_basis(system)
    use constants, only:atomname
    implicit none
    type(ElectronicSystem), intent(inout) :: system

    ! workspace
    integer :: atom_Z, iatom, bas_ptr
    integer :: nbas=0, nprim=0, norb=0,  di
    integer :: ishell, iprims, iaos
    integer, allocatable :: bas(:,:)

    ! Count all the basis functions
    if (verbosity .ge. 0) then
       write(*,*) '                 nshells   norbs    nprims'
    end if
    do iatom=1,system%natm
       atom_Z = system%atm(CHARGE_OF,iatom)
       bas_ptr = paos(PTR_ATOM, atom_Z)
       if (bas_ptr < 1) then
          print '(2a)', 'integrals_build_basis: no basis function for element ',atomname(atom_Z)
          error stop -1
       end if
       ishell = 0                       ! nshells
       iprims = 0                       ! primitives
       iaos = 0                       ! aos
       basis_loop:do while (.true.)
          if (pbas(ATOM_OF, bas_ptr) .eq. atom_Z) then
             ishell = ishell + 1
             di = CINTcgto_spheric(bas_ptr-1, pbas)
             iaos = iaos + di
             iprims = iprims + pbas(NPRIM_OF, bas_ptr) * di

             bas_ptr = bas_ptr + 1
          else
             exit basis_loop
          end if
       end do basis_loop

       nbas = nbas + ishell
       norb = norb + iaos
       nprim = nprim + iprims

       if (verbosity > 0) then
          write(*,'(a,i3,2a,i3,a,i3,a,i3,a,i3)') &
               '      ', iatom,atomname(atom_Z),&
               '       ', ishell,&
               '       ', iaos, &
               '       ', iprims
       end if
    end do

    ! allocate shells
    allocate(bas(BAS_SLOTS, nbas))

    ! build basis and atoms
    ishell = 1
    do iatom=1, system%natm
       atom_Z = system%atm(CHARGE_OF, iatom)
       bas_ptr = paos(PTR_ATOM, atom_Z)
       basis_loop2:do while (.true.)
          if (pbas(ATOM_OF, bas_ptr) .eq. atom_Z) then
             bas(:,ishell) = pbas(:,bas_ptr)
             bas(ATOM_OF,ishell) = iatom-1 ! zero-based indexing

             if (bas(NCTR_OF,ishell) .ne. 1) then
                call log_err('integrals_build_basis', 'uncontracted shells not supported yet!')
                error stop -1
             end if

             bas_ptr = bas_ptr + 1
             ishell = ishell + 1
          else
             exit basis_loop2
          end if
       end do basis_loop2
    end do

    if (verbosity .ge. 0) then
       write(*,'(a,i4,a,i4,a,i4)') '       Tot.      ', &
            nbas, '      ', norb,  '      ', nprim
    end if

    ! Write results to system
    system%bas = bas
    system%nbas = nbas
    system%nprim = nprim
    system%norb = norb
  end subroutine integrals_build_basis

  subroutine integrals_build_atomic_orbitals(system)
    ! Load atomic orbitals for system.
    use constants, only: atomname
    implicit none
    type(ElectronicSystem), intent(inout) :: system

    ! workspace
    integer :: i, atom_Z, n, ii, jj, kk, k, off, j

    ! ouptuts
    double precision, allocatable :: C_AO(:,:), E_AO(:), pop_AO(:)
    integer :: naos=0

    do i=1,system%natm
       naos = naos + paos(NAOS_OF, system%atm(CHARGE_OF, i))
    end do

    if (verbosity > 0) then
       write(*,*) '          atom    ao energies '
    end if
   allocate(C_AO(naos, system%norb))
   allocate(E_AO(naos))
   allocate(pop_AO(naos))

   C_AO = 0.0
   k = 1
   do i=1,system%natm
      atom_Z = system%atm(CHARGE_OF, i)
      n = paos(NAOS_OF,atom_Z)
      ii = paos(PTR_AOS_C, atom_Z) + 1
      jj = paos(PTR_AOS_E, atom_Z) + 1
      kk = paos(PTR_AOS_OCC, atom_Z) + 1

      if (verbosity > 0) write(*,'(a,i3,a)')  '          ',i,atomname(atom_Z)
      ! where the atom group is
      off = k                   ! row offset

      do j=0,n-1
         C_AO(off:off+n-1,k) = ENV(ii:ii+n -1)
         E_AO(k) = ENV(jj+j)
         pop_AO(k) = ENV(kk+j)
         ii = ii + n
         k = k + 1
         if (verbosity >0) write(*, '(a, f8.2)')      '                     ', ENV(jj+j)
      end do
   end do

   system%C_AO = C_AO
   system%E_AO = E_AO
   system%pop_AO = pop_AO
   system%naos = naos
 end subroutine integrals_build_atomic_orbitals

  subroutine integrals_overlaps(system, S)
    implicit none
    type(ElectronicSystem), intent(in) :: system
    double precision, allocatable, intent(out) :: S(:,:)

    ! workspace
    double precision, allocatable :: buf1e(:,:), val
    integer :: ii, jj, k, i,j, di, dj, i2,j2
    integer :: shls(4)

    if(allocated(S)) deallocate(S)
    allocate(S(1:system%norb, 1:system%norb))

    ii = 1
    k = 0                           !todo only L or U part of matrix?
    do i=0, system%nbas-1           ! 0 based index
       shls(1) = i; di = CINTcgto_spheric(i, system%bas)

       jj = 1
       do j=0, system%nbas-1    ! 0 based index
          shls(2) = j; dj = CINTcgto_spheric(j, system%bas)

          call cint1e_ovlp_sph(S(ii:ii+di-1,jj:jj+dj-1), shls, &
               system%atm, system%natm, system%bas, system%nbas, ENV)

!!$          allocate(buf1e(di,dj))
!!$          call cint1e_ovlp_sph(buf1e(1:di,1:dj), shls, &
!!$               system%atm, system%natm, system%bas, system%nbas, ENV)
!!$
!!$          do i2=1,di
!!$             do j2=1,dj
!!$                S(ii+i2-1,jj+j2-1) = buf1e(i2,j2)
!!$             end do
!!$          end do
!!$
!!$          deallocate(buf1e)
          k = k + di*dj
          jj = jj + dj
       end do

       ii = ii + di
    end do

    ! Compute norm of basis
    val = 0.0
    do i=1,system%norb
       val = val + S(i,i)
    end do

    if (verbosity.ge.0) then
      write(*,*) '     computed overlaps:',k
      write(*,'(a, f6.4)') '      basis trace = ', val/dble(system%norb)
    end if
    ! todo check if basis trace == norbs
  end subroutine integrals_overlaps

  subroutine integrals_symm_bas2ao(system, O, O_proj)
    ! Project a (real and symmetric) operator in the Gaussian basis to the AO basis.
    implicit none
    type(ElectronicSystem), intent(in) :: system
    double precision, intent(in) :: O(:,:)
    double precision,allocatable, intent(out) :: O_proj(:,:)

    ! workspace
    integer :: naos, norb
    double precision, allocatable :: tmp(:,:)

    naos = system%naos
    norb = system%norb
    allocate(O_proj(naos,naos))
    allocate(tmp(naos,norb))

    ! Compute: tmp = C_AO(naos,norb) * O(norb,norb)
    call dsymm('R','U',naos, norb, 1.0d0, O, norb, system%C_AO, naos, 0.0d0, tmp, naos)

    ! Compute: O_proj = tmp * C_AO^T
    call dgemm('n', 't', naos, naos, norb, 1.0d0, tmp, naos, system%C_AO, naos, 0.0d0, O_proj, naos)
  end subroutine integrals_symm_bas2ao

  subroutine integrals_symm_ao2bas(system, O, O_proj)
    ! Project a (real and symmetric) operator in the AO basis to the Gaussian basis.
    implicit none
    type(ElectronicSystem), intent(in) :: system
    double precision, intent(in) :: O(:,:)
    double precision,allocatable, intent(out) :: O_proj(:,:)

    ! workspace
    integer :: naos, norb
    double precision, allocatable :: tmp(:,:)

    naos = system%naos
    norb = system%norb
    allocate(O_proj(naos,naos))
    allocate(tmp(naos,norb))

    ! Compute: tmp(naos, norb) = O(naos,naos) * C_AO(naos,norb)
    call dsymm('L','U', naos, norb, 1.0d0, O, norb, system%C_AO, naos, 0.0d0, tmp, naos)

    ! Compute: O_proj(norb, norb) = C(naos,norb)^T tmp(norb, naos)
    call dgemm('t', 'n', norb, norb, naos, 1.0d0, system%C_AO, naos, tmp, norb, 0.0d0, O_proj, norb)
  end subroutine integrals_symm_ao2bas

  subroutine integrals_MO_AO_transform(system, MO_AOs, MO_bas)
    ! Project MOs in the AO basis to the original basis
    implicit none
    type(ElectronicSystem), intent(in) :: system
    double precision, intent(in) :: MO_AOs(:,:)
    double precision,allocatable, intent(out) :: MO_bas(:,:)

    ! workspace
    integer :: naos, norb

    naos = system%naos
    norb = system%norb

   ! Compute: MO = C_AO(naos, norb).T * eigenvectors(naos, naos)
   ! Note that the eigenvectors are returned in the old Hamiltonian matrix.
    if (allocated(MO_bas)) deallocate(MO_bas)
    allocate(MO_bas(norb,naos))
    call dgemm('t', 'n', norb, naos, naos, 1.0d0, system%C_AO, naos, MO_AOs, naos, 0.0d0, MO_bas, norb)
  end subroutine integrals_MO_AO_transform

  subroutine integrals_unroll(system, C_MO, unrolled_MOs)
    use constants, only: atomname
    implicit none
    type(ElectronicSystem), intent(in) :: system
    double precision, intent(in) :: C_MO(1:,1:)
    type(UnrolledMOs), intent(out) :: unrolled_MOs

    ! workspace
    double precision, allocatable ::  coeffs(:), coeffs_cart(:)
    integer :: nbas, natm, naos, i, j, k, ii, jj, imo, n
    integer :: ncart, nsph

    naos = system%naos
    natm = system%natm
    nbas = system%nbas

    ! compute how many cartesian basis we will have
    ncart = 0
    do i=1, nbas
       ncart = ncart + CINTcgto_cart(i-1, system%bas) * system%bas(NPRIM_OF, i)
    end do

    unrolled_MOs%nmos = naos
    unrolled_MOs%natm = natm
    unrolled_MOs%nprim = ncart

    allocate(unrolled_MOs%atoms(natm))
    allocate(unrolled_MOs%Z(natm))
    allocate(unrolled_MOs%xyz(natm, 3))

    do i=1,natm
       unrolled_MOs%atoms(i) = atomname(system%atm(CHARGE_OF,i))
       unrolled_MOs%Z(i) = dble(system%atm(CHARGE_OF,i))
       unrolled_MOs%xyz(i,:) = env(system%atm(PTR_COORD,i)+1 : system%atm(PTR_COORD,i)+3)
    end do

    allocate(unrolled_MOs%icnt(ncart))
    allocate(unrolled_MOs%ityp(ncart))
    allocate(unrolled_MOs%exps(ncart))

    ii = 1                          ! Prim index
    jj = 1                          ! AO index
    do i=1,nbas                                  !shells
       do j=1, system%bas(NPRIM_OF, i)             !primitives
          do k=1,CINTcgto_cart(i-1, system%bas)       !number of cartesian functions
             unrolled_MOs%icnt(ii) = system%bas(ATOM_OF,i) + 1
             unrolled_MOs%exps(ii) = env(system%bas(PTR_EXP,i) + j)
             ! this follows the .wfn / AIM  standard (or should anyway!)
             unrolled_MOs%ityp(ii) = angular_start(system%bas(ANG_OF, i)) + k

             ! update primitive counter
             ii = ii + 1
          end do

          ! update AO index
          jj = jj + 1
       end do
    end do

    allocate(coeffs(1:11))
    allocate(coeffs_cart(1:21))
    allocate(unrolled_MOs%coeffs(ncart, naos))

    do imo=1,naos
       ii = 1
       jj = 1
       do i=1,nbas                                  !shells
          nsph = CINTcgto_spheric(i-1, system%bas)
          do j=1, system%bas(NPRIM_OF, i)             !primitives
             ! ii is the spheric primitive index, jj is the cartesian primitive index

             ! coefficient of primitives * norm * MO coefficients
             coeffs(1:nsph) = cintgto_norm(system%bas(ANG_OF,i), &
                  env(system%bas(PTR_EXP,i) + j) ) &
                  * env(system%bas(PTR_COEFF,i) + j) * C_MO(ii : ii + nsph - 1, imo)

             ! convert to cartesian representation
             ncart = CINTcgto_cart(i-1, system%bas)
             call sph2car(system%bas(ANG_OF,i), coeffs(1:nsph), coeffs_cart(1:ncart))

             unrolled_MOs%coeffs(jj : jj + ncart - 1, imo) = coeffs_cart(1:ncart)
             jj = jj + ncart
          end do

          ii = ii + nsph
       end do
    end do

  contains
    pure integer function angular_start(i)
      ! Return the starting index for the basis type in a wfn file of a basis
      ! function with angular momentum i
      integer, intent(in) :: i

      select case(i)
      case(0)
         angular_start=0
      case(1)
         angular_start=1
      case(2)
         angular_start=4
      case(3)
         angular_start=10
      case(4)
         angular_start=20
      case(5)
         angular_start=35
      case default
         angular_start = -1
      end select
    end function angular_start
  end subroutine integrals_unroll

  subroutine integrals_write_to_wfn(iwfn, mos, E_MO, occ)
    use constants, only: atomname
    implicit none
    integer, intent(in) :: iwfn
    double precision, intent(in) :: E_MO(1:)
    double precision, intent(in) :: occ(1:)
    type(UnrolledMOs), intent(in) :: mos

    ! workspace
    integer :: i, j, k

    ! This next bit is adapted from molden2aim
    ! thanks open source !
    write(iwfn,*) ' Extended Huckel Theory orbitals'
    write(iwfn,"('GAUSSIAN',8x,i7,' MOL ORBITALS',i7,' PRIMITIVES',i9,' NUCLEI')") &
         mos%nmos, mos%nprim, mos%natm

    do i=1,mos%natm
       write(iwfn, "(2x,A3,i3,4x,'(CENTRE',i3,') ',3f12.8,'  CHARGE =',f5.1)") &
            mos%atoms(i), i, i, &
            mos%xyz(i,1:3), mos%Z(i)
    end do

    write(iwfn,"('CENTRE ASSIGNMENTS  ',20i3)")(mos%icnt(i), i=1, mos%nprim)
    write(iwfn,"('TYPE ASSIGNMENTS    ',20i3)")(mos%ityp(i), i=1, mos%nprim)
    write(iwfn,"('EXPONENTS ',5d14.7)")(mos%exps(i), i=1, mos%nprim)

    do i=1,mos%nmos
       write(iwfn,"('MO',i5,5x,'MO 0.0',8x,'OCC NO =',f13.7,'  ORB. ENERGY =',f12.6)") &
            i, occ(i), E_MO(i)

       write(iwfn,"(5d16.8)")(mos%coeffs(j,i), j=1,mos%nprim)
    end do

    write(iwfn,"('END DATA')")
    write(iwfn,"(' THE  HF ENERGY =',f20.12,' THE VIRIAL(-V/T)=',f13.8)")0.d0,2.d0
    close(iwfn)
  end subroutine integrals_write_to_wfn

  ! from molden2aim -> convert spherical basis to cartesian
  Subroutine sph2car(lq,fi,fo)
    implicit real(kind=8) (a-h,o-z)
    integer :: lq, i, j
    dimension         :: fi(*),fo(*)

    select case(lq)
    case(0)    ! S
       fo(1)=fi(1)
    case(1)    ! P
       fo(1:3)=fi(1:3)
    case(2)    ! D
       fo(1:6)=0.0d0
       do j=1,5
          do i=1,6
             fo(i)=fo(i)+s2cd(i,j)*fi(j)
          end do
       end do
    case(3)    ! F
       fo(1:10)=0.0d0
       do j=1,7
          do i=1,10
             fo(i)=fo(i)+s2cf(i,j)*fi(j)
          end do
       end do
    case(4)    ! G
       fo(1:15)=0.0d0
       do j=1,9
          do i=1,15
             fo(i)=fo(i)+s2cg(i,j)*fi(j)
          end do
       end do
    case(5)    ! H
       fo(1:21)=0.0d0
       do j=1,11
          do i=1,21
             fo(i)=fo(i)+s2ch(i,j)*fi(j)
          end do
       end do
    case default
       write(*,*) 'sph2car: angular momentum invalid', lq
       error stop -1
    end select

    return
  end Subroutine sph2car

Subroutine set_s2c_coef()
 implicit none

   s2cd = 0.0D0
   s2cd( 1, 1) = -0.5000000000000000D+00
   s2cd( 2, 1) = -0.5000000000000000D+00
   s2cd( 3, 1) =  1.0D0
   s2cd( 5, 2) =  1.0D0
   s2cd( 6, 3) =  1.0D0
   s2cd( 1, 4) =  0.8660254037844387D+00
   s2cd( 2, 4) = -0.8660254037844387D+00
   s2cd( 4, 5) =  1.0D0

   s2cf = 0.0D0
   s2cf( 3, 1) =  1.0D0
   s2cf( 6, 1) = -0.6708203932499368D+00
   s2cf( 9, 1) = -0.6708203932499368D+00
   s2cf( 1, 2) = -0.6123724356957946D+00
   s2cf( 4, 2) = -0.2738612787525830D+00
   s2cf( 7, 2) =  0.1095445115010332D+01
   s2cf( 2, 3) = -0.6123724356957946D+00
   s2cf( 5, 3) = -0.2738612787525830D+00
   s2cf( 8, 3) =  0.1095445115010332D+01
   s2cf( 6, 4) =  0.8660254037844386D+00
   s2cf( 9, 4) = -0.8660254037844386D+00
   s2cf(10, 5) =  1.0D0
   s2cf( 1, 6) =  0.7905694150420949D+00
   s2cf( 4, 6) = -0.1060660171779821D+01
   s2cf( 2, 7) = -0.7905694150420949D+00
   s2cf( 5, 7) =  0.1060660171779821D+01

   s2cg = 0.0D0
   s2cg( 1, 1) =  0.3750000000000000D+00
   s2cg( 2, 1) =  0.3750000000000000D+00
   s2cg( 3, 1) =  1.0D0
   s2cg(10, 1) =  0.2195775164134200D+00
   s2cg(11, 1) = -0.8783100656536799D+00
   s2cg(12, 1) = -0.8783100656536799D+00
   s2cg( 5, 2) = -0.8964214570007953D+00
   s2cg( 8, 2) =  0.1195228609334394D+01
   s2cg(14, 2) = -0.4008918628686366D+00
   s2cg( 7, 3) = -0.8964214570007953D+00
   s2cg( 9, 3) =  0.1195228609334394D+01
   s2cg(13, 3) = -0.4008918628686366D+00
   s2cg( 1, 4) = -0.5590169943749475D+00
   s2cg( 2, 4) =  0.5590169943749475D+00
   s2cg(11, 4) =  0.9819805060619659D+00
   s2cg(12, 4) = -0.9819805060619659D+00
   s2cg( 4, 5) = -0.4225771273642583D+00
   s2cg( 6, 5) = -0.4225771273642583D+00
   s2cg(15, 5) =  0.1133893419027682D+01
   s2cg( 5, 6) =  0.7905694150420950D+00
   s2cg(14, 6) = -0.1060660171779821D+01
   s2cg( 7, 7) = -0.7905694150420950D+00
   s2cg(13, 7) =  0.1060660171779821D+01
   s2cg( 1, 8) =  0.7395099728874520D+00
   s2cg( 2, 8) =  0.7395099728874520D+00
   s2cg(10, 8) = -0.1299038105676658D+01
   s2cg( 4, 9) =  0.1118033988749895D+01
   s2cg( 6, 9) = -0.1118033988749895D+01

   s2ch = 0.0D0
   s2ch( 1, 1) =  1.0D0
   s2ch( 3, 1) = -0.1091089451179962D+01
   s2ch( 5, 1) =  0.625D+00
   s2ch(12, 1) = -0.1091089451179962D+01
   s2ch(14, 1) =  0.3659625273557000D+00
   s2ch(19, 1) =  0.625D+00
   s2ch( 7, 2) =  0.1290994448735806D+01
   s2ch( 9, 2) = -0.5669467095138407D+00
   s2ch(11, 2) =  0.1613743060919757D+00
   s2ch(16, 2) = -0.1267731382092775D+01
   s2ch(18, 2) =  0.2112885636821291D+00
   s2ch(21, 2) =  0.4841229182759271D+00
   s2ch( 2, 3) =  0.1290994448735806D+01
   s2ch( 4, 3) = -0.1267731382092775D+01
   s2ch( 6, 3) =  0.4841229182759271D+00
   s2ch(13, 3) = -0.5669467095138407D+00
   s2ch(15, 3) =  0.2112885636821291D+00
   s2ch(20, 3) =  0.1613743060919757D+00
   s2ch( 3, 4) = -0.1118033988749895D+01
   s2ch( 5, 4) =  0.8539125638299665D+00
   s2ch(12, 4) =  0.1118033988749895D+01
   s2ch(19, 4) = -0.8539125638299665D+00
   s2ch( 8, 5) =  0.1290994448735806D+01
   s2ch(10, 5) = -0.6454972243679028D+00
   s2ch(17, 5) = -0.6454972243679028D+00
   s2ch( 9, 6) = -0.1224744871391589D+01
   s2ch(11, 6) =  0.5229125165837972D+00
   s2ch(16, 6) =  0.9128709291752769D+00
   s2ch(18, 6) =  0.2282177322938192D+00
   s2ch(21, 6) = -0.5229125165837972D+00
   s2ch( 4, 7) = -0.9128709291752769D+00
   s2ch( 6, 7) =  0.5229125165837972D+00
   s2ch(13, 7) =  0.1224744871391589D+01
   s2ch(15, 7) = -0.2282177322938192D+00
   s2ch(20, 7) = -0.5229125165837972D+00
   s2ch( 5, 8) =  0.7395099728874520D+00
   s2ch(14, 8) = -0.1299038105676658D+01
   s2ch(19, 8) =  0.7395099728874520D+00
   s2ch(10, 9) = -0.1118033988749895D+01
   s2ch(17, 9) =  0.1118033988749895D+01
   s2ch(11,10) =  0.1169267933366857D+01
   s2ch(18,10) = -0.1530931089239487D+01
   s2ch(21,10) =  0.7015607600201140D+00
   s2ch( 6,11) =  0.7015607600201140D+00
   s2ch(15,11) = -0.1530931089239487D+01
   s2ch(20,11) =  0.1169267933366857D+01
end Subroutine set_s2c_coef

end module integrals

max_stars_repo_name: clavigne/huckster

max_issues_repo_name: clavigne/huckster

max_forks_repo_name: clavigne/huckster


text:
c
c***************************************************************************
c
      subroutine invert(rho_in, theta_in, capr_out, capz_out,
     . nnodex, nnodey, nxmx, nymx, rho, theta, capr, capz, myid,
     . drhodr, drhodz, dthedr, dthedz, xj)

      implicit none
      integer myid, nnodex, nnodey, nxmx, nymx, n, m, nmin, mmin
      real rho_in, theta_in, capr_out, capz_out, diff_rho, diff_theta,
     .   diff, diffmin
      real  capr(nxmx), capz(nymx), rk, zk,
     .     rho(nxmx, nymx), theta(nxmx, nymx),
     .  drhodr(nxmx, nymx), drhodz(nxmx, nymx),
     .  dthedr(nxmx, nymx), dthedz(nxmx, nymx)
      real  xj, thetak, rhok, dthedrk, drhodzk, drhodrk, dthedzk

      diffmin = 1.0

      do n = 1, nnodex
         do m = 1, nnodey
	    diff_rho = rho_in - rho(n,m)
	    diff_theta = theta_in - theta(n,m)
	    diff = sqrt(diff_rho**2 + diff_theta**2)

	    if(diff .le. diffmin)then
	       diffmin = diff
	       nmin = n
	       mmin = m
	    end if

	 end do
      end do

      n = nmin
      m = mmin

      rk = capr(n)
      zk = capz(m)
      rhok = rho(n,m)
      thetak = theta(n,m)
      drhodrk = drhodr(n,m)
      drhodzk = drhodz(n,m)
      dthedrk = dthedr(n,m)
      dthedzk = dthedz(n,m)

      xj = dthedrk * drhodzk - drhodrk * dthedzk

      if(xj .eq. 0)then
         capr_out = rk
	 capz_out = zk
	 return
      end if

      capr_out = rk + 1. / xj * (drhodzk * (theta_in - thetak)
     .                         - dthedzk * (rho_in   - rhok)  )
      capz_out = zk + 1. / xj * (dthedrk * (rho_in   - rhok)
     .                         - drhodrk * (theta_in - thetak))

      return
      end

c
c***************************************************************************
c
      subroutine checkx(myid, i, j, nkx1, nkx2, nky1, nky2,
     1   fdk, fek, ffk, u, v, w,
     1   xb, xc, xd, nkdim1, nkdim2, mkdim1, mkdim2,
     1   nnodex, nnodey, nxmx, nymx, xprime, yprime,
     1   idiag, jdiag, isweep)

      implicit none

      integer nkx1, nkx2, nky1, nky2, nkdim1, nkdim2, mkdim1, mkdim2,
     1        nnodex, nnodey, nxmx, nymx, idiag, jdiag, isweep,
     1        i, j, n, m, myid

      complex fdk(nkdim1 : nkdim2, mkdim1 : mkdim2),
     1        fek(nkdim1 : nkdim2, mkdim1 : mkdim2),
     1        ffk(nkdim1 : nkdim2, mkdim1 : mkdim2),
     1        xb(nxmx, nymx), term1,
     1        xc(nxmx, nymx), term2,
     1        xd(nxmx, nymx), term3,
     1        xlhs, xrhs
      complex u(nkdim1 : nkdim2, mkdim1 : mkdim2),
     1        v(nkdim1 : nkdim2, mkdim1 : mkdim2),
     1        w(nkdim1 : nkdim2, mkdim1 : mkdim2)

      real    xprime(nxmx), yprime(nymx)

      if (myid .eq. 0) then
         write(15, 910)
         write(15, 930)
         write(15, 920)
         write(15, 940)
         write(15, 920)

         write(6, 910)
         write(6, 930)
         write(6, 920)
         write(6, 940)
         write(6, 920)

         term1 = 0.0
         term2 = 0.0
         term3 = 0.0

         do n = nkx1, nkx2
            do m = nky1, nky2

            term1 = term1 + fdk(n, m) * u(n, m)
            term2 = term2 + fek(n, m) * v(n, m)
            term3 = term3 + ffk(n, m) * w(n, m)

            end do
         end do

         xlhs = term1 + term2 + term3
         xrhs = xb(i, j)
         write(15, 900) i, j, xprime(i), yprime(j), xlhs, xrhs
         write(6,  900) i, j, xprime(i), yprime(j), xlhs, xrhs

	 write(15, 901) term1, term2, term3
         write(6,  901) term1, term2, term3

	 write(15, 901) u(64,64), v(64,64), w(64,64)
         write(6,  901) u(64,64), v(64,64), w(64,64)

	 write(15, 901) fdk(64,64), fek(64,64), ffk(64,64)
         write(6,  901) fdk(64,64), fek(64,64), ffk(64,64)

      endif

      return
  900 format(i5, i5, 1p9e12.3)
  901 format(10x, 1p9e12.3)
  940 format(1h , 4h   i, 5h    j, 4x, 6h x(i) , 6x, 6h y(j) ,
     1                    6x,  6hre lhs, 6x, 6him lhs, 6x,
     1                         6hre rhs, 6x, 6him rhs)
  910 format (1h1)
  920 format (1h0)
  930 format(3x, 14h u diagnostics)
      end

c
c***************************************************************************
c
      subroutine checky(myid, i, j, nkx1, nkx2, nky1, nky2,
     1   fgk, fak, fpk, u, v, w,
     1   xb, xc, xd, nkdim1, nkdim2, mkdim1, mkdim2,
     1   nnodex, nnodey, nxmx, nymx, xprime, yprime,
     1   idiag, jdiag, isweep)

      implicit none

      integer nkx1, nkx2, nky1, nky2, nkdim1, nkdim2, mkdim1, mkdim2,
     1        nnodex, nnodey, nxmx, nymx, idiag, jdiag, isweep,
     1        i, j, n, m, myid

      complex fgk(nkdim1 : nkdim2, mkdim1 : mkdim2),
     1        fak(nkdim1 : nkdim2, mkdim1 : mkdim2),
     1        fpk(nkdim1 : nkdim2, mkdim1 : mkdim2),
     1        xb(nxmx, nymx), term1,
     1        xc(nxmx, nymx), term2,
     1        xd(nxmx, nymx), term3,
     1        xlhs, xrhs
      complex u(nkdim1 : nkdim2, mkdim1 : mkdim2),
     1        v(nkdim1 : nkdim2, mkdim1 : mkdim2),
     1        w(nkdim1 : nkdim2, mkdim1 : mkdim2)

      real    xprime(nxmx), yprime(nymx)

      if (myid .eq. 0) then
         write(15, 910)
         write(15, 930)
         write(15, 920)
         write(15, 940)
         write(15, 920)

         write(6, 910)
         write(6, 930)
         write(6, 920)
         write(6, 940)
         write(6, 920)

         term1 = 0.0
         term2 = 0.0
         term3 = 0.0

         do n = nkx1, nkx2
            do m = nky1, nky2

            term1 = term1 + fgk(n, m) * u(n, m)
            term2 = term2 + fak(n, m) * v(n, m)
            term3 = term3 + fpk(n, m) * w(n, m)

            end do
         end do

         xlhs = term1 + term2 + term3
         xrhs = xc(i, j)
         write(15, 900) i, j, xprime(i), yprime(j), xlhs, xrhs
         write(6,  900) i, j, xprime(i), yprime(j), xlhs, xrhs

         write(15, 901) term1, term2, term3
         write(6,  901) term1, term2, term3

	 write(15, 901) u(64,64), v(64,64), w(64,64)
         write(6,  901) u(64,64), v(64,64), w(64,64)

	 write(15, 901) fgk(64,64), fak(64,64), fpk(64,64)
         write(6,  901) fgk(64,64), fak(64,64), fpk(64,64)
      endif

      return
  900 format(i5, i5, 1p9e12.3)
  901 format(10x, 1p9e12.3)
  940 format(1h , 4h   i, 5h    j, 4x, 6h x(i) , 6x, 6h y(j) ,
     1                    6x,  6hre lhs, 6x, 6him lhs, 6x,
     1                         6hre rhs, 6x, 6him rhs)
  910 format (1h1)
  920 format (1h0)
  930 format(3x, 14h v diagnostics)
      end

c
c***************************************************************************
c
      subroutine checkz(myid, i, j, nkx1, nkx2, nky1, nky2,
     1   frk, fqk, fsk, u, v, w,
     1   xb, xc, xd, nkdim1, nkdim2, mkdim1, mkdim2,
     1   nnodex, nnodey, nxmx, nymx, xprime, yprime,
     1   idiag, jdiag, isweep)

      implicit none

      integer nkx1, nkx2, nky1, nky2, nkdim1, nkdim2, mkdim1, mkdim2,
     1        nnodex, nnodey, nxmx, nymx, idiag, jdiag, isweep,
     1        i, j, n, m, myid

      complex frk(nkdim1 : nkdim2, mkdim1 : mkdim2),
     1        fqk(nkdim1 : nkdim2, mkdim1 : mkdim2),
     1        fsk(nkdim1 : nkdim2, mkdim1 : mkdim2),
     1        xb(nxmx, nymx), term1,
     1        xc(nxmx, nymx), term2,
     1        xd(nxmx, nymx), term3,
     1        xlhs, xrhs
      complex u(nkdim1 : nkdim2, mkdim1 : mkdim2),
     1        v(nkdim1 : nkdim2, mkdim1 : mkdim2),
     1        w(nkdim1 : nkdim2, mkdim1 : mkdim2)

      real    xprime(nxmx), yprime(nymx)

      if (myid .eq. 0) then
         write(15, 910)
         write(15, 930)
         write(15, 920)
         write(15, 940)
         write(15, 920)

         write(6, 910)
         write(6, 930)
         write(6, 920)
         write(6, 940)
         write(6, 920)

         term1 = 0.0
         term2 = 0.0
         term3 = 0.0

         do n = nkx1, nkx2
            do m = nky1, nky2

            term1 = term1 + frk(n, m) * u(n, m)
            term2 = term2 + fqk(n, m) * v(n, m)
            term3 = term3 + fsk(n, m) * w(n, m)

            end do
         end do

         xlhs = term1 + term2 + term3
         xrhs = xd(i, j)
         write(15, 900) i, j, xprime(i), yprime(j), xlhs, xrhs
         write(6,  900) i, j, xprime(i), yprime(j), xlhs, xrhs

	 write(15, 901) term1, term2, term3
         write(6,  901) term1, term2, term3

	 write(15, 901) u(64,64), v(64,64), w(64,64)
         write(6,  901) u(64,64), v(64,64), w(64,64)

	 write(15, 901) frk(64,64), fqk(64,64), fsk(64,64)
         write(6,  901) frk(64,64), fqk(64,64), fsk(64,64)

      endif

      return
  900 format(i5, i5, 1p9e12.3)
  901 format(10x, 1p9e12.3)
  940 format(1h , 4h   i, 5h    j, 4x, 6h x(i) , 6x, 6h y(j) ,
     1                    6x,  6hre lhs, 6x, 6him lhs, 6x,
     1                         6hre rhs, 6x, 6him rhs)
  910 format (1h1)
  920 format (1h0)
  930 format(3x, 14h w diagnostics)
      end

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

max_stars_repo_name: ORNL-Fusion/aorsa

max_issues_repo_name: ORNL-Fusion/aorsa

max_forks_repo_name: ORNL-Fusion/aorsa


text:
less_toxic(n1,i1).
less_toxic(ee1,jj1).
less_toxic(m1,dd1).
less_toxic(y1,d1).
less_toxic(cc1,v1).
less_toxic(m1,i1).
less_toxic(w1,l1).
less_toxic(bb1,l1).
less_toxic(j1,d1).
less_toxic(q1,i1).
less_toxic(ee1,i1).
less_toxic(m1,z1).
less_toxic(cc1,jj1).
less_toxic(cc1,t1).
less_toxic(kk1,d1).
less_toxic(n1,t1).
less_toxic(kk1,a1).
less_toxic(y1,h1).
less_toxic(y1,f1).
less_toxic(n1,z1).
less_toxic(t1,c1).
less_toxic(ll1,cc1).
less_toxic(bb1,a1).
less_toxic(ee1,w1).
less_toxic(j1,v1).
less_toxic(m1,l1).
less_toxic(ii1,i1).
less_toxic(k1,w1).
less_toxic(j1,aa1).
less_toxic(l1,f1).
less_toxic(kk1,h1).
less_toxic(k1,l1).
less_toxic(x1,w1).
less_toxic(s1,b1).
less_toxic(bb1,t1).
less_toxic(bb1,z1).
less_toxic(ll1,h1).
less_toxic(ff1,a1).
less_toxic(r1,c1).
less_toxic(g1,t1).
less_toxic(q1,f1).
less_toxic(ll1,z1).
less_toxic(y1,e1).
less_toxic(h1,d1).
less_toxic(w1,v1).
less_toxic(v1,f1).
less_toxic(o1,i1).

max_stars_repo_name: manoelfranca/cilppp

max_issues_repo_name: manoelfranca/cilppp

max_forks_repo_name: manoelfranca/cilppp


text:
      subroutine shape_deform_b(ibeg, jbeg, iend, jend, nhhp,
     1                          idim, jdim, r, rw, drw, drwb,
     2                          nwp, xw, xwb, idx, param_type)
      implicit none
      include 'dim.h'
      integer ibeg, jbeg, iend, jend, nhhp, idim, jdim, nwp
      real    r   (idim,jdim,NDIM),
     2        rw  (NDIM,*),
     3        drw (NDIM,*),
     4        xw  (*)
      real    drwb(NDIM,*),
     4        xwb (*)
      integer idx(*)
      integer param_type

      integer i, j, nc, iw, iinc, jinc
      real    x1, y1, x2, y2, ln, xp, yp, xf, yf, l, t,
     1        nx, ny, dh, dhb, A, B, C, D, h
      real    HicksHenne
      real    kulfan

c     Set increment to +1 or -1
      iinc = 1
      jinc = 1
      if(iend.lt.ibeg) iinc = -1
      if(jend.lt.jbeg) jinc = -1

c     (x1,y1) = first point, (x2,y2) = second point on curve
      if(param_type .eq. 1)then
         x1 = r(ibeg,jbeg,1)
         y1 = r(ibeg,jbeg,2)
         x2 = r(iend,jend,1)
         y2 = r(iend,jend,2)
      elseif(param_type .eq. 2)then
         x1 = 1.0
         y1 = 0.0
         x2 = 0.0
         y2 = 0.0
      elseif(param_type .eq. 3)then
         print*,'shape_deform_b.f:  param_type=3 is not finished'
         stop
      endif
      ln = sqrt( (x2-x1)**2 + (y2-y1)**2 )
      if(ln.eq.0.0 .and. nhhp.ne.0)then
         print*,'shape_deform: ln is <= 0'
         stop
      else if(ln.eq.0.0 .and. nhhp.eq.0)then
         ln = 1.0
      else if(ln.eq.0.0)then
         print*,'shape_deform: OOPS'
         print*,'ln =',ln
         print*,'nhhp =',nhhp
         stop
      endif
      nx = (y2-y1)/ln
      ny =-(x2-x1)/ln

c     Equation of line joining (x1,y1) to (x2,y2) -> Ax + By + C = 0
      A  = (y2 - y1)
      B  =-(x2 - x1)
      C  = -(y2 - y1)*x1 + (x2 - x1)*y1

      nc          = 0
      xwb(1:nhhp) = 0.0

      do i=ibeg,iend,iinc
         do j=jbeg,jend,jinc
            xp        = r(i,j,1)
            yp        = r(i,j,2)

            if(nhhp.eq.0) goto 100

c           Equation of perpendicular line through (xp,yp)
c           -B x + A y + D = 0
            D = B * xp - A * yp
c           Solve for intersection point (xf,yf) = foot of perpendicular
            xf = -(A*C - B*D)/(A**2 + B**2)
            yf = -(B*C + A*D)/(A**2 + B**2)
            l  = sqrt( (xf-x1)**2 + (yf-y1)**2 )
            t         = l/ln
            nc        = nc + 1
            iw        = idx(nc)
            print*,nc,iw
            if(param_type .eq. 1)then
               dh        = HicksHenne(nhhp, xw, t)
c              Reverse sweep
               dhb       = drwb(1,iw)*nx + drwb(2,iw)*ny
               call HicksHenne_b(nhhp, xw, xwb, dhb, t)
            elseif(param_type .eq. 2)then
               t = xp
               h = kulfan(nhhp, xw, t)
c              Reverse sweep
               dhb = drwb(2,iw)
               call kulfan_b(nhhp, xw, xwb, t, dhb)
            endif

100         continue
         enddo
      enddo

      return
      end


! { dg-do compile }
! Tests the fix for PR31630, in which the association of the argument
! of 'cmp' did not work.
!
! Contributed by Francois-Xavier Coudert <fxcoudert@gcc.gnu.org>
!
module box_module
  type box
    integer :: m = 0
  end type box
end module box_module

module sort_box_module
contains

  subroutine heapsort_box(cmp)
    interface
       subroutine cmp(a)
         use box_module
         type(box) :: a
       end subroutine cmp
    end interface
    optional :: cmp
  end subroutine heapsort_box

end module sort_box_module

module boxarray_module
  use box_module
  implicit none

  type boxarray
    type(box), allocatable :: bxs(:)
  end type boxarray
contains

  subroutine boxarray_build_l(ba)
    type(boxarray) :: ba
    allocate(ba%bxs(1))
  end subroutine boxarray_build_l

  subroutine boxarray_sort()
    use sort_box_module
    call heapsort_box
  end subroutine boxarray_sort

end module boxarray_module


text:
Student Political Action Committee, or SPAC, is an organization dedicated to strengthening relationships between the Associated Students of UC Davis and the undergraduate student body which it serves.

SPAC is committed to empowering students not otherwise involved in ASUCD to better their campus, lobbying the association to maximize student governments effect on the undergraduate population, discovering and developing new political talent on campus, and writing meaningful legislation.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


subroutine par
  integer ix, jx

  !$acc parallel
  !$acc loop tile (*,*) ! { dg-error "not enough DO loops for tiled" }
  do ix = 1, 30
  end do

  !$acc loop tile (*,*)
  do ix = 1, 30
     do jx = 1, ix ! { dg-error "tiled loops don.t form rectangular" }
     end do
  end do

  !$acc loop tile (*)
  do ix = 1, 30
     do jx = 1, ix
     end do
  end do
  !$acc end parallel
end subroutine par


text:
program boj1259
  implicit none
  integer :: n
  do
    read*, n
    if (n .eq. 0) then
      exit
    end if
    if (is_palindrome(n)) then
      print '(A)', 'yes'
    else
      print '(A)', 'no'
    end if
  end do

  contains
  logical function is_palindrome(i)
    integer :: arr(100), last, i
    logical :: ok
    last = -1
    ok = .true.
    do while (i > 0)
      arr(last + 1) = mod(i, 10)
      i = i / 10
      last = last + 1
    end do
    do i = 0, last / 2
      ok = ok .and. (arr(i) == arr(last - i))
    end do
    is_palindrome = ok
  end function is_palindrome
end program boj1259

max_stars_repo_name: doublik/HYU

max_issues_repo_name: doublik/HYU

max_forks_repo_name: doublik/HYU


text:
      subroutine mdrgex(adhfs,adhfsd,adhfsx,advfs,advfsd,advfsx,
     $ axhfs,axhfsd,axhfsx,axlks,axlksd,axlksx,axvfs,axvfsd,axvfsx,
     $ cdrs,cdrsd,cdrsx,ipch,ipchmx,ipcv,ipcvmx,narxmx,nbasp,nbaspd,
     $ nbtmax,ndrs,ndrsd,ndrsx,ndrsmx,ndrsr,ndrsrd,ndrsrx,nern1,
     $ nern2,noutpt,nst,nstmax,ntprmx,nttyo)
c
c     This subroutine folds the reactions and reaction properties for
c     generic ion exchangers into the 'd' set. This is generally done
c     after those reactions and properties have been manipulated.
c
c     Note that the modified arrays must first be written into a set
c     of scratch arrays. This is because the ion exchanger sections
c     must generally go in the middle of the new arrays, and their
c     lengths are generally such that the positioning of the start of
c     the terminal "unmodified" sections is modified.
c
c     This subroutine is called by:
c
c       EQLIB/chsgex.f
c
c-----------------------------------------------------------------------
c
c     Principal input:
c
c       adhfs  = standard array of coefficients for computing
c                  pressure derivatives of enthalpy functions as a
c                  function of temperature
c       advfs  = standard array of coefficients for computing
c                  pressure derivatives of volume functions as a
c                  function of temperature
c       axhfs  = standard array of coefficients for computing
c                  enthalpy functions as a function of temperature
c       axlks  = standard array of coefficients for computing
c                  equilibrium constants as a function of temperature
c       axvfs  = standard array of coefficients for computing
c                  volume functions as a function of temperature
c       cdrs   = standard array of reaction coefficients
c       narxmx = number of coefficient elements of axlks per species
c                  per temperature range
c       nbtmax = maximum number of basis species
c       ndrs   = standard  array of indices of species corresponding to
c                  reaction coefficients
c       ndrsr  = standard pointer array for indices of species appearing
c                  in reactions
c       nstmax = maximum number of species
c       ntprmx = maximum number of temperature ranges
c
c     Principal output:
c
c       adhfsd = 'data file' array of coefficients for computing
c                  pressure derivatives of enthalpy functions as a
c                  function of temperature
c       advfsd = 'data file' array of coefficients for computing
c                  pressure derivatives of volume functions as a
c                  function of temperature
c       axhfsd = 'data file' array of coefficients for computing
c                  enthalpy functions as a function of temperature
c       axlksd = 'data file' array of coefficients for computing
c                  equilibrium constants as a function of temperature
c       axvfsd = 'data file' array of coefficients for computing
c                  volume functions as a function of temperature
c       cdrsd  = 'data file' array of reaction coefficients
c       ndrsd  = 'data file'  array of indices of species corresponding
c                  to reaction coefficients
c       ndrsrd = 'data file' pointer array for indices of species
c                  appearing in reactions
c
c-----------------------------------------------------------------------
c
      implicit none
c
c-----------------------------------------------------------------------
c
c     Calling sequence variable declarations.
c
      integer ipchmx,ipcvmx,narxmx,nbtmax,ndrsmx,nstmax,ntprmx
c
      integer noutpt,nttyo
c
      integer nbasp(nbtmax),nbaspd(nbtmax),ndrs(ndrsmx),ndrsd(ndrsmx),
     $ ndrsx(ndrsmx),ndrsr(2,nstmax),ndrsrd(2,nstmax),ndrsrx(2,nstmax)
c
      integer ipch,ipcv,nern1,nern2,nst
c
      real*8 adhfs(narxmx,ntprmx,ipchmx,nstmax),
     $ adhfsd(narxmx,ntprmx,ipchmx,nstmax),
     $ adhfsx(narxmx,ntprmx,ipchmx,nstmax),
     $ advfs(narxmx,ntprmx,ipcvmx,nstmax),
     $ advfsd(narxmx,ntprmx,ipcvmx,nstmax),
     $ advfsx(narxmx,ntprmx,ipcvmx,nstmax),
     $ axhfs(narxmx,ntprmx,nstmax),axhfsd(narxmx,ntprmx,nstmax),
     $ axhfsx(narxmx,ntprmx,nstmax),axlks(narxmx,ntprmx,nstmax),
     $ axlksd(narxmx,ntprmx,nstmax),axlksx(narxmx,ntprmx,nstmax),
     $ axvfs(narxmx,ntprmx,nstmax),axvfsd(narxmx,ntprmx,nstmax),
     $ axvfsx(narxmx,ntprmx,nstmax),cdrs(ndrsmx),cdrsd(ndrsmx),
     $ cdrsx(ndrsmx)
c
c-----------------------------------------------------------------------
c
c     Local variable declarations.
c
      integer i,ipc,j,jhmax,jvmax,j2,j3,n,nmax,nb,nn,nr1,nr2,ns,nt
c
      integer ilnobl
c
      character*8 ux8a,ux8b
c
c-----------------------------------------------------------------------
c
      nmax = narxmx*ntprmx*(nern1 -1)
      jhmax = narxmx*ntprmx*ipchmx*(nern1 -1)
      jvmax = narxmx*ntprmx*ipcvmx*(nern1 -1)
c
c* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
c     Coefficients for equilibrium constants.
c
      if (nern1 .gt. 1) call copyaa(axlksd,axlksx,nmax)
c
      do ns = nern1,nern2
        do j = 1,ntprmx
          do i = 1,narxmx
            axlksx(i,j,ns) = axlks(i,j,ns)
          enddo
        enddo
      enddo
c
      if (nern2 .lt. nst) then
        do ns = nern2 + 1,nst
          do j = 1,ntprmx
            do i = 1,narxmx
              axlksx(i,j,ns) = axlksd(i,j,ns)
            enddo
          enddo
        enddo
      endif
c
c* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
c     Coefficients related to enthalpy.
c
      if (ipch .ge. 0) then
c
c       Coefficients for enthalpy functions.
c
        if (nern1 .gt. 1) call copyaa(axhfsd,axhfsx,nmax)
c
        do ns = nern1,nern2
          do j = 1,ntprmx
            do i = 1,narxmx
              axhfsx(i,j,ns) = axhfs(i,j,ns)
            enddo
          enddo
        enddo
c
        if (nern2 .lt. nst) then
          do ns = nern2 + 1,nst
            do j = 1,ntprmx
              do i = 1,narxmx
                axhfsx(i,j,ns) = axhfsd(i,j,ns)
              enddo
            enddo
          enddo
        endif
c
        if (ipch .ge. 1) then
c
c         Coefficients for derivatives of enthalpy functions.
c
          if (nern1 .gt. 1) call copyaa(adhfsd,adhfsx,jhmax)
c
          do ns = nern1,nern2
            do ipc = 1,ipchmx
              do j = 1,ntprmx
                do i = 1,narxmx
                  adhfsx(i,j,ipc,ns) = adhfs(i,j,ipc,ns)
                enddo
              enddo
            enddo
          enddo
c
          if (nern2 .lt. nst) then
            do ns = nern2 + 1,nst
              do ipc = 1,ipchmx
                do j = 1,ntprmx
                  do i = 1,narxmx
                    adhfsx(i,j,ipc,ns) = adhfsd(i,j,ipc,ns)
                  enddo
                enddo
              enddo
            enddo
          endif
        endif
      endif
c
c* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
c     Coefficients related to volume.
c
      if (ipcv .ge. 0) then
c
c         Coefficients for volume functions.
c
        if (nern1 .gt. 1) call copyaa(axvfsd,axvfsx,nmax)
c
        do ns = nern1,nern2
          do j = 1,ntprmx
            do i = 1,narxmx
              axvfsx(i,j,ns) = axvfs(i,j,ns)
            enddo
          enddo
        enddo
c
        if (nern2 .lt. nst) then
          do ns = nern2 + 1,nst
            do j = 1,ntprmx
              do i = 1,narxmx
                axvfsx(i,j,ns) = axvfsd(i,j,ns)
              enddo
            enddo
          enddo
        endif
c
        if (ipcv .ge. 1) then
c
c         Coefficients for derivatives of volume functions.
c
          if (nern1 .gt. 1) call copyaa(advfsd,advfsx,jvmax)
c
          do ns = nern1,nern2
            do ipc = 1,ipcvmx
              do j = 1,ntprmx
                do i = 1,narxmx
                  advfsx(i,j,ipc,ns) = advfs(i,j,ipc,ns)
                enddo
              enddo
            enddo
          enddo
c
          if (nern2 .lt. nst) then
            do ns = nern2 + 1,nst
              do ipc = 1,ipcvmx
                do j = 1,ntprmx
                  do i = 1,narxmx
                    advfsx(i,j,ipc,ns) = advfsd(i,j,ipc,ns)
                  enddo
                enddo
              enddo
            enddo
          endif
        endif
      endif
c
c* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
c     Index range pointers for reactions.
c
      if (nern1 .gt. 1) then
        nmax = 2*(nern1 - 1)
        call copyia(ndrsrd,ndrsrx,nmax)
      endif
c
      do ns = nern1,nern2
        ndrsrx(1,ns) = ndrsr(1,ns)
        ndrsrx(1,ns) = ndrsr(1,ns)
      enddo
c
      if (nern2 .lt. nst) then
        do ns = nern2 + 1,nst
          ndrsrx(1,ns) = ndrsrd(1,ns)
          ndrsrx(1,ns) = ndrsrd(1,ns)
        enddo
      endif
c
c* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
c     Reaction coefficients and indices of associated species.
c
      nt = 0
      if (nern1 .gt. 1) nt = ndrsrd(2,nern1 - 1)
      if (nern2 .ge. nern1)
     $  nt = nt + ndrsr(2,nern2) - ndrsr(1,nern1) + 1
      if (nst .gt. nern2)
     $ nt = nt + ndrsrd(2,nst) - ndrsrd(1,nern2 + 1) + 1
c
      if (nt .gt. ndrsmx) then
        write (ux8a,'(i5)') ndrsmx
        write (ux8b,'(i5)') nt
        j2 = ilnobl(ux8a)
        j3 = ilnobl(ux8b)
        write (noutpt,1000) ux8a(1:j2),ux8b(1:j3)
        write (nttyo,1000)  ux8a(1:j2),ux8b(1:j3)
 1000   format(/' * Error - (EQLIB/mdrgex.f) Have insufficient array',
     $  ' space to expand the',/7x,'reaction coefficient arrays to',
     $  ' accommodate changes to the',/7x,'ion exchanger section',
     $  " of the 'd' set of reaction data. Increase",/7x,'the',
     $  ' dimensioning parameter ndrspa from ',a,' to at least ',
     $  a,'.')
        stop
      endif
c
      if (nern1 .gt. 1) then
        nmax = ndrsrd(2,nern1 - 1)
        call copyaa(cdrsd,cdrsx,nmax)
        call copyia(ndrsd,ndrsx,nmax)
      endif
c
      if (nern1 .gt. 1) then
        nn = ndrsrd(2,nern1 - 1)
      else
        nn = 0
      endif
      do ns = nern1,nern2
        nr1 = ndrsr(1,ns)
        nr2 = ndrsr(2,ns)
        do n = nr1,nr2
          nn = nn + 1
          ndrsx(nn) = ndrs(n)
          cdrsx(nn) = cdrs(n)
        enddo
      enddo
c
      nn = ndrsrx(2,nern2)
      do ns = nern2 + 1,nst
        nr1 = ndrsrd(1,ns)
        nr2 = ndrsrd(2,ns)
        do n = nr1,nr2
          nn = nn + 1
          ndrsx(nn) = ndrsd(n)
          cdrsx(nn) = cdrsd(n)
        enddo
      enddo
c
c* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
      do nb = 1,nbtmax
        ns = nbasp(nb)
        if (ns.ge.nern1 .and. ns.le.nern2) nbaspd(nb) = nbasp(nb)
      enddo
c
c* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
c     Copy the scratch arrays into the corresponding 'd' set arrays.
c
c     Copy the new reactions in the 'x' set into the standard arrays
c     in the 'd' set.
c
c     Calling sequence substitutions:
c       adhfsd for adhfs
c       advfsd for advfs
c       axhfsd for axhfs
c       axlksd for axlks
c       axvfsd for avhfs
c       cdrsd for cdrs
c       ndrsd for ndrs
c       ndrsrd for ndrsr
c
      call cdrscx(adhfsd,adhfsx,advfsd,advfsx,axhfsd,axhfsx,axlksd,
     $ axlksx,axvfsd,axvfsx,cdrsd,cdrsx,ipch,ipchmx,ipcv,ipcvmx,
     $ narxmx,ndrsd,ndrsmx,ndrsx,ndrsrd,ndrsrx,nstmax,ntprmx)
c
c* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
      end


text:
           program demo_log10
           use, intrinsic :: iso_fortran_env, only : real_kinds, &
           & real32, real64, real128
           implicit none
           real(kind=real64) :: x = 10.0_real64
             x = log10(x)
           end program demo_log10

max_stars_repo_name: urbanjost/fortran-intrinsic-manpages

max_issues_repo_name: urbanjost/fortran-intrinsic-manpages

max_forks_repo_name: urbanjost/fortran-intrinsic-manpages


program basic
!! QUICK PROTOTYPE: JUST THE BARE ESSENTIALS
use M_CLI,  only : commandline, check_commandline
implicit none

integer            ::  x,y,z      ;  namelist  /args/  x,y,z
integer            ::  casen      ;  namelist  /args/  casen
integer            ::  ints(4)=0  ;  namelist  /args/  ints
character(len=20)  ::  string     ;  namelist  /args/  string

character(len=:),allocatable :: command
character(len=4096)          :: cmd
logical                      :: nq=.false.
integer                      :: e

   command=repeat(' ',4096)
   call get_command(command)
   command=trim(command)//' '
   write(*,'(*(g0,1x))')'COMMAND=',command

   call get_command_argument(0,cmd)
   e=len_trim(cmd)+1
   command='-x 10 -y 20 -z 30 -ints 11,22,33 -casen 0 -string "My string,""again"""'
   call readcli()

   select case(casen)
   case(0)
      call printit(all([x,y,z,ints].eq.[10,20,30,11,22,33,0]))
      call execute_command_line(cmd(:e)//'-x 4 -y 5 -z 6 -casen 1')
   case(1)
      call printit(all([x,y,z].eq.[4,5,6]))
      call execute_command_line(cmd(:e)//'-x 40 -y 50 -z 60 -casen 2')
   case(2)
      call printit(all([x,y,z].eq.[40,50,60]))
      call execute_command_line(cmd(:e)//'-x 400 -y 500 -z 600 -ints -1,-2,-3 -casen 3')
   case(3)
      call printit(all([x,y,z,ints].eq.[400,500,600,-1,-2,-3,0]))
      call execute_command_line(cmd(:e)//'-x 400 -y 500 -z 600 -ints -1,-2,-3 -casen 900')
   case(900)
      write(*,*)'USAGE'
      call execute_command_line(cmd(:e)//'--casen 901 --usage')
      call execute_command_line(cmd(:e)//'--casen 901')
   case(901)
      write(*,*)'HELP'
      call execute_command_line(cmd(:e)//'--casen 902 --help')
      call execute_command_line(cmd(:e)//'--casen 902')
   case(902)
      write(*,*)'VERSION'
      call execute_command_line(cmd(:e)//'--casen 999 --version')
      call execute_command_line(cmd(:e)//'--casen 999')
   case(7)
   case(8)
   case(9)
   case(999)
   case default
      write(*,'(a)')'default - should not get here'
      stop
   end select

contains

subroutine printit(testit)
logical testit
   write(*,'(*(g0,1x))',advance='no')merge('PASSED:','FAILED:',testit)
   write(*,nml=args)
   write(*,'(/,a)')repeat('=',132)
end subroutine printit

subroutine readcli()
character(len=:),allocatable :: readme ! stores updated namelist
character(len=256)           :: message
integer                      :: ios
   !!write(*,'(*(g0,1x))')'SET DEFAULTS',trim(command)
   readme=commandline(command,noquote=nq)
   !!write(*,'(*(g0,1x))')'NAMELIST STRING=',trim(readme)
   read(readme,nml=args,iostat=ios,iomsg=message)
   call check_commandline(ios,message)
end subroutine readcli

end program basic


text:
module hifuku
	implicit none
	integer, parameter :: imax=1000,jmax=1000
	real, parameter :: dx=10.0,dy=10.0,pie=3.1415927
	integer,dimension(imax,jmax) :: VAL,VISIT
	real,dimension(imax,jmax) :: DEM

contains
    subroutine over_view(is,js,ie0,je0,ie1,je1)
	integer is,js,ie0,je0,ie1,je1,iend,jend
	integer x,y,xs,xe,ys,ye
	integer inc,check
	real xx,yy
	real zs,ze0,ze1,za,zb
	real a,b,c

	zs=DEM(is,js) ; ze0=DEM(ie0,je0) ; ze1=DEM(ie1,je1)
	a=((zs-ze0)*(js-je1)-(zs-ze1)*(js-je0))/((is-ie0)*(js-je1)-(is-ie1)*(js-je0))
	b=((zs-ze0)*(is-ie1)-(zs-ze1)*(is-ie0))/((js-je0)*(is-ie1)-(js-je1)*(is-ie0))
	c=zs-a*is-b*js
	if (je0 .eq. je1) then
		if (js .lt. je0) then
			inc=1 ; jend=jmax
		else
			inc=-1 ; jend=1
		endif
		do y=je0+inc,jend,inc
			yy=real(y)
			if (ie0.lt.ie1) then
				xs=is+(ie0-is)*(yy-js)/(je0-js)+1.0
				xe=is+(ie1-is)*(yy-js)/(je1-js)
			else
				xe=is+(ie0-is)*(yy-js)/(je0-js)
				xs=is+(ie1-is)*(yy-js)/(je1-js)+1.0
			endif
			if(xs.gt.imax) exit
			if(xe.lt.1) exit
			if(xe.gt.imax) xe=imax
			if(xs.lt.1) xs=1
			check=1
			do x=xs,xe
				za=DEM(x,y)
				zb=a*x+b*y+c
				if (zb .ge. za) then
					VAL(x,y)=0 ; check=0
				endif
			enddo
			if(check.eq.1) exit
		enddo
	endif
	if (ie0 .eq. ie1) then
		if (is .lt. ie0) then
			inc=1 ; iend=imax
		else
			inc=-1 ; iend=1
		endif
		do x=ie0+inc,iend,inc
			xx=real(x)
			if(je0.lt.je1) then
				ys=js+(je0-js)*(xx-is)/(ie0-is)+1.0
				ye=js+(je1-js)*(xx-is)/(ie1-is)
			else
				ye=js+(je0-js)*(xx-is)/(ie0-is)
				ys=js+(je1-js)*(xx-is)/(ie1-is)+1.0
			endif
			if(ys.gt.jmax) exit
			if(ye.lt.1) exit
			if(ye.gt.jmax) ye=jmax
			if(ys.lt.1) ys=1
			check=1
			do y=ys,ye
				za=DEM(x,y)
				zb=a*x+b*y+c
				if (zb .ge. za) then
					VAL(x,y)=0 ; check=0
				endif
			enddo
			if(check.eq.1) exit
		enddo
	endif
	return
    end subroutine over_view

    subroutine view(is,js)
	integer is,js
	integer xe,ye,xe0,ye0
	integer dir,count

	VISIT=0 ; VAL=1
	VISIT(is,js)=1
	xe=is+1 ; ye=is
	xe0=is+1 ; ye0=js ; VISIT(xe0,ye0)=1
	dir=1 ; count=2
	do while(count.le.((imax-4)*(jmax-4)-1))
	    xe=xe0 ; ye=ye0
	    select case(dir)
	      case(0)
	        xe=xe0+1
	      case(1)
	        ye=ye0+1
	      case(2)
	        xe=xe0-1
	      case(3)
	        ye=ye0-1
            end select
	    if (VISIT(xe,ye).eq.1) then
		xe=xe0 ; ye=ye0
		dir=mod(dir+3,4)
	    	select case(dir)
	   	   case(0)
	   	     xe=xe0+1
	 	   case(1)
	 	     ye=ye0+1
	 	   case(2)
	   	     xe=xe0-1
	   	   case(3)
	   	     ye=ye0-1
            	end select
	    endif
	    if ((xe.gt.2).and.(xe.lt.(imax-1)).and.(ye.gt.2).and.(ye.lt.(jmax-1))) then
		if ((VAL(xe0,ye0).eq.1).or.(VAL(xe,ye).eq.1)) call over_view(is,js,xe0,ye0,xe,ye)
	    	VISIT(xe,ye)=1
	    	count=count+1
	    endif
	    xe0=xe ; ye0=ye
	    dir=mod(dir+1,4)
	end do
    end subroutine view
end module hifuku

max_stars_repo_name: y-iikura/ViewShed

max_issues_repo_name: y-iikura/ViewShed

max_forks_repo_name: y-iikura/ViewShed


! Empty submodule for checking disambiguation of
! nodes at the same vertical level in the tree
submodule ( parent ) sibling
contains
  module function sibling_function() result(sibling_stuff)
    logical :: sibling_stuff
    sibling_stuff=.true.
  end function
end submodule sibling


text:
submodule (BandSpectrum_mod) EnergyFluence_smod

    implicit none

    real(RK)                        :: mv_alphaPlusOne
    real(RK)                        :: mv_alphaPlusTwoOverEpk

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

contains

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    !> \brief
    !> Integrate the Band differential spectrum over the input energy range in units of the input energy.
    !>
    !> \param[in]   lowerLim        :   The lower limit energy (in units of [keV]) of the integration.
    !> \param[in]   upperLim        :   The upper limit energy (in units of [keV]) of the integration.
    !> \param[in]   epk             :   The spectral peak energy in units of [keV].
    !> \param[in]   alpha           :   The lower spectral exponent of the Band model.
    !> \param[in]   beta            :   The upper spectral exponent of the Band model.
    !> \param[in]   tolerance       :   The relative accuracy tolerance of the integration.
    !> \param[out]  energyFluence   :   The fluence in units of the input energy (keV) within the input energy range.
    !> \param[out]  Err             :   An object of class [Err_type](@ref err_mod::err_type) containing error-handling information.
    module subroutine getEnergyFluence(lowerLim,upperLim,epk,alpha,beta,tolerance,energyFluence,Err)
#if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
        !DEC$ ATTRIBUTES DLLEXPORT :: getEnergyFluence
#endif

        use Constants_mod, only: IK, RK, HUGE_RK
        use QuadPackSPR_mod, only: qag
        use Err_mod, only: Err_type
        implicit none
        real(RK), intent(in)            :: lowerLim, upperLim, epk, alpha, beta, tolerance
        real(RK), intent(out)           :: energyFluence
        type(Err_type), intent(out)     :: Err
        character(*), parameter         :: PROCEDURE_NAME = "@getEnergyFluence()"
        real(RK)                        :: ebrk, alphaPlusTwo
        real(RK)                        :: thisUpperLim, betaPlusTwo
        real(RK)                        :: alphaMinusBeta, coef
        real(RK)                        :: abserr
        integer(IK)                     :: neval
        integer(IK)                     :: ierr

        Err%occurred = .false.

        if (lowerLim>=upperLim) then
            energyFluence = 0._RK
            return
        end if

        ! check if the photon indices are consistent with the mathematical rules
        if (alpha<beta .or. alpha<-2._RK) then
            energyFluence = -HUGE_RK
            Err%occurred = .true.
            Err%msg = MODULE_NAME // PROCEDURE_NAME // ": Error occurred: alpha<beta .or. alpha<-2._RK"
            return
        end if

        ! integrate the spectrum
        alphaPlusTwo = alpha + 2._RK
        alphaMinusBeta = alpha - beta
        ebrk = epk*alphaMinusBeta/alphaPlusTwo

        if (lowerLim>ebrk) then

            ! there is only the high energy component in the energyFluence
            betaPlusTwo = beta + 2._RK
            coef = ebrk**alphaMinusBeta * exp(-alphaMinusBeta);
            energyFluence = coef * ( upperLim**betaPlusTwo - lowerLim**betaPlusTwo ) / betaPlusTwo
            return

!#if defined OS_IS_WSL && defined CODECOV_ENABLED
!! LCOV_EXCL_START
!#endif

        elseif (lowerLim<ebrk) then

            mv_alphaPlusTwoOverEpk = alphaPlusTwo / epk
            thisUpperLim = min(upperLim,ebrk)
            mv_alphaPlusOne = alpha + 1._RK
            ! use brute-force integration
            call qag( f             = getBandCompLowEnergy  &
                    , a             = lowerLim              &
                    , b             = thisUpperLim          &
                    , epsabs        = 0._RK                 &
                    , epsrel        = tolerance             &
                    , key           = 1_IK                  &
                    , result        = energyFluence         &
                    , abserr        = abserr                &
                    , neval         = neval                 &
                    , ier           = ierr                  &
                    )

            if (ierr/=0_IK) then
                energyFluence = -HUGE_RK
                Err%occurred = .true.
                Err%stat = ierr
                Err%msg = MODULE_NAME // PROCEDURE_NAME // ": Error occurred at QuadPack routine. Check the error code to identify the root cause."
                return
            end if

            if (upperLim>ebrk) then ! add the remaining part of the energyFluence from the high-energy component
                betaPlusTwo = beta + 2._RK
                alphaMinusBeta = alpha - beta
                coef = ebrk**alphaMinusBeta * exp(-alphaMinusBeta);
                energyFluence = energyFluence + coef * ( upperLim**betaPlusTwo - ebrk**betaPlusTwo ) / betaPlusTwo
                return
            end if

        end if

    end subroutine getEnergyFluence

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    pure function getBandCompLowEnergy(energy) result(bandCompLow)
#if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
        !DEC$ ATTRIBUTES DLLEXPORT :: getBandCompLowEnergy
#endif
        implicit none
        real(RK), intent(in)    :: energy
        real(RK)                :: bandCompLow
        bandCompLow = energy**mv_alphaPlusOne * exp(-mv_alphaPlusTwoOverEpk*energy)
    end function getBandCompLowEnergy

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

end submodule EnergyFluence_smod ! LCOV_EXCL_LINE

max_stars_repo_name: cdslaborg/paramonte

max_issues_repo_name: cdslaborg/paramonte

max_forks_repo_name: cdslaborg/paramonte


! ******************************************************************************************************************************** !
! ocnlite_box.f90
!
! MISCELLANEOUS ROUTINES
! ******************************************************************************************************************************** !

MODULE ocnlite_box

  USE ocnlite_lib
  IMPLICIT NONE
  SAVE

CONTAINS

  ! ****************************************************************************************************************************** !
  ! DUMMY
  SUBROUTINE sub_DUMMY()
    ! #### INSERT CODE ... ####################################################################################################### !
    ! ############################################################################################################################ !
  END SUBROUTINE sub_DUMMY
  ! ****************************************************************************************************************************** !

END MODULE ocnlite_box


text:
      subroutine diag_allocate(pointer,size_block,size_column,
     :            logical_success,size)

      pointer(pointer,array(1))
      logical logical_success;
      integer size_block, iblock, size_column,ierr_mem,size;

      logical_success = .true.

cgd   pointer=malloc(size_block*size);
      call alloc(pointer,size_block,size);
      if (pointer == 0) then
        call free(pointer)
cgd     pointer = malloc(size_column*size);
        call alloc(pointer,size_column,size);
        if (pointer.ne.0) logical_success = .false.
        if (pointer == 0) then
            write(0,*) 'Insufficient Memory: Stop in diag_alloc_memory'
        end if
      end if

      end subroutine diag_allocate


text:
subroutine VorIVel(VFlag,CalcDer,G,X1,Y1,Z1,X2,Y2,Z2,XP,YP,ZP,UP,VP,WP,DUDX)

    use vortex

    real :: DVFDX, RRAT
    integer :: VFlag       ! 0 for bound vortex element, 1 for trailing wake element, 2 for spanwise wake element (sets core radius)
    integer :: CalcDer     ! 1 to calc dudx

    AX = X2-X1
    AY = Y2-Y1
    AZ = Z2-Z1
    A2 = AX**2 + AY**2 + AZ**2

    BX = X2-XP
    BY = Y2-YP
    BZ = Z2-ZP
    B = SQRT(BX**2 + BY**2 + BZ**2)

    CX = X1-XP
    CY = Y1-YP
    CZ = Z1-ZP
    C = SQRT(CX**2 + CY**2 + CZ**2)

    if (B == 0.0) then
        B = 1.0
    end if
    if (C == 0.0) then
        C = 1.0
    end if

    ADB = (AX*BX + AY*BY + AZ*BZ)
    ADBDB = ADB/B
    ADC = (AX*CX + AY*CY + AZ*CZ)
    ADCDC = ADC/C

    CCAX = CY*AZ - AY*CZ
    CCAY = AX*CZ - CX*AZ
    CCAZ = CX*AY - AX*CY

    CCAV = CCAX**2 + CCAY**2 + CCAZ**2

    ! Select the vortex core radius for bound vortex, trailing
    ! or spanwise wake elements
    if (VFlag == 1) then
        VRAD2=VRAD2_T
    else if (VFlag == 2) then
        VRAD2=VRAD2_S
    else
        VRAD2=VRAD2_B
    end if

    if (CCAV >= vCutOffRad) then
        if (ivtxcor > 0 .and. CCAV < A2*VRAD2) then
            if (ivtxcor == 1) then
                ! Constant velocity (approx) core
                CCAV = A2*VRAD2 ! limit denominator to value at core radius
                VF   = (ADBDB-ADCDC)*G/(12.56637*CCAV)
                if (CalcDer == 1) then
                    DVFDX = G/12.56637*((AX/C-CX*ADC/C**3-AX/B+BX*ADB/B**3)/CCAV)
                end if
            else if (ivtxcor == 2) then
                ! Linear velocity (approx) core
                RRAT=sqrt(CCAV/(A2*VRAD2)) ! radius to core radius ratio
                CCAV = A2*VRAD2 ! limit denominator to value at core radius
                VF   = RRAT*(ADBDB-ADCDC)*G/(12.56637*CCAV)
                if (CalcDer == 1) then
                    DVFDX = G/12.56637*((AX/C-CX*ADC/C**3-AX/B+BX*ADB/B**3)/CCAV)
                end if
            end if
        else
            VF   = (ADBDB-ADCDC)*G/(12.56637*CCAV)
            if (CalcDer == 1) then
                DVFDX = G/12.56637*((AX/C-CX*ADC/C**3-AX/B+BX*ADB/B**3)/CCAV-2*(ADBDB-ADCDC)/CCAV**2*(AZ*CCAY-AY*CCAZ))
            end if
        end if
    else
        VF   = 0.0
        DVFDX = 0.0
    end if

    UP = UP+CCAX*VF
    VP = VP+CCAY*VF
    WP = WP+CCAZ*VF
    if (CalcDer == 1) then
        DUDX = DUDX+CCAX*DVFDX
    end if

    return
end subroutine VorIVel

max_stars_repo_name: ebranlard/CACTUS

max_issues_repo_name: ebranlard/CACTUS

max_forks_repo_name: ebranlard/CACTUS


text:
	SUBROUTINE GQCTYP  ( szctyp, ictwid, iret )
C************************************************************************
C* GQCTYP								*
C* 									*
C* This subroutine returns the cloud type symbol attributes, including  *
C* the cloud symbol size and line width. 				*
C* 									*
C* GQCTYP  ( SZCTYP, ICTWID, IRET )					*
C*									*
C* Output parameters:							*
C* 	SZCTYP		REAL		Cloud symbol size 		*
C*	ICTWID		INTEGER		Line width			*
C*	IRET		INTEGER		Return code			*
C**									*
C* Log:									*
C* S. Schotz/GSC	 3/90						*
C* M. desJardins/NMC	10/91	Fixed return code			*
C* L. Williams/EAi       3/94   Removed blank comments from header      *
C* A. Hardy/GSC          6/98   Cleaned up prolog                       *
C************************************************************************
	INCLUDE		'FUNCCODE.PRM'
	INCLUDE		'ERROR.PRM'
C*
	INTEGER		isend (2)
C------------------------------------------------------------------------
C*	Load input parameters into buffer and write them to the mailbox.
C
	isend (1) = 2
	isend (2) = FQCTYP
C
	CALL GPUT  ( isend, 2, iret )
	IF  ( iret .ne. NORMAL )  RETURN
C
C*	Get output parameters
C
	CALL GGET  ( iret, 1, ierr )
	IF  ( ierr .ne. NORMAL )  THEN
	    iret = ierr
	    RETURN
	END IF
C
	CALL GGETR  ( szctyp, 1, ierr )
	IF  ( iret .eq. NORMAL )  THEN
	    CALL GGET  ( ictwid, 1, ierr )
	END IF
	IF  ( ierr .ne. NORMAL )  iret = ierr
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
! DART software - Copyright UCAR. This open source software is provided
! by UCAR, "as is", without charge, subject to all terms of use at
! http://www.image.ucar.edu/DAReS/DART/DART_download
!
! $Id$

module dart_time_io_mod

!> \defgroup dart_time_io_mod dart_time_io_mod
!> Default routines for netCDF reading and writing dart model time.
!> If your model uses a different name for the time dimension
!> or has a different way of handing/storing time, it must provide
!> a custom read_model_time() and write_model_time() routine.
!> @{

use types_mod,            only : r8, digits12
use time_manager_mod,     only : time_type, set_time, get_time, print_time, &
                                 set_calendar_type, set_date, get_calendar_string, &
                                 operator(+)

use utilities_mod,        only : E_MSG, E_ERR, error_handler, to_upper
use netcdf_utilities_mod, only : nc_check, nc_open_file_readonly, nc_close_file
use typeSizes
use netcdf

implicit none
private

public :: read_model_time, write_model_time

! version controlled file description for error handling, do not edit
character(len=*), parameter :: source   = "$URL$"
character(len=*), parameter :: revision = "$Revision$"
character(len=*), parameter :: revdate  = "$Date$"

character(len=512) :: string1, string2, string3

contains

!--------------------------------------------------------------------
!> Make a stab at reading the time from the input file
!> and setting the calendar.  This routine assumes the file hasn't
!> been opened yet, and closes it when it's done.

function read_model_time(filename)

character(len=*), intent(in) :: filename
type(time_type)              :: read_model_time

integer :: ncid, ios, numdims, xtype, VarID
integer :: ntimes, seconds, days
integer :: year, month, day, hour, minute, second
type(time_type) :: base_time, delta_time

character(len=*), parameter :: routine = 'read_model_time'
real(digits12) :: model_time, time_array(1)

integer, dimension(NF90_MAX_VAR_DIMS) :: dimIDs
character(len=256) :: file_calendar, dart_calendar
character(len=256) :: unitstring

! this is used in many error messages below.  set it here, and
! don't reuse string3 here, please.
write(string3,*)'You may need to supply a model-specific "read_model_time()" to read the time.'

ncid = nc_open_file_readonly(filename, routine)

ios = nf90_inq_varid(ncid, "time", VarID)
if (ios /= NF90_NOERR) then
   write(string1,*)'Expecting to be able to read the model time from a "time" variable'
   write(string2,*)'in file "'//trim(filename)//'"'
   call error_handler(E_ERR,'read_model_time', string1, &
              source, revision, revdate, text2=string2, text3=string3)
endif

ios = nf90_inquire_variable(ncid, VarID, xtype=xtype, dimids=dimIDs, ndims=numdims)
call nc_check(ios, 'read_model_time', 'inquire_variable "time" from "'//trim(filename)//'"')

if (numdims > 1) then
   write(string1,*)'Expecting the "time" variable to be a single dimension.'
   write(string2,*)'in file "'//trim(filename)//'"'
   call error_handler(E_ERR,'read_model_time', string1, &
              source, revision, revdate, text2=string2, text3=string3)
endif

if (numdims == 0) then
   ios = nf90_get_var(ncid, VarID, model_time)
   call nc_check(ios, 'read_model_time','get_var scalar time' )
else

   ! Since the time variable is known to have only 1 dimension, we know it is the first one.
   ios = nf90_inquire_dimension(ncid, dimIDs(1), len=ntimes)
   call nc_check(ios, 'read_model_time', 'inquire_dimension for time dimension from "'//trim(filename) )

   ! read the last one
   ios = nf90_get_var(ncid, VarID, time_array, start=(/ntimes/), count=(/1/))
   call nc_check(ios, 'read_model_time','get_var time' )
   model_time = time_array(1)
endif

! try to handle the calendar in a generic way

call get_calendar_string(dart_calendar)

ios = nf90_get_att(ncid, VarID, 'calendar', file_calendar)
if (ios /= NF90_NOERR) file_calendar = 'NO_CALENDAR'
call to_upper(file_calendar)

if (dart_calendar /= file_calendar ) then
   ! allow NO_CALENDAR, NO CALENDAR, and NONE to be synonyms.
   ! replace file_calendar with what dart uses to simplify the tests below.
   if (dart_calendar == 'NO_CALENDAR' .and. (file_calendar == 'NONE'        .or. &
                                             file_calendar == 'NO CALENDAR')) then
      file_calendar = 'NO_CALENDAR'
   else
      write(string1,*)'inconsistent calendar types between DART program and input file.'
      write(string2,*)'DART initialized with: ', trim(dart_calendar), ' File uses: ', trim(file_calendar)
      call error_handler(E_ERR, 'read_model_time:', string1, source,revision, revdate, &
                      text2=string2, text3=string3)
   endif
endif

if ( dart_calendar == 'NO_CALENDAR' ) then

   !> assumes time variable is real and fractional days.  if this isn't true,
   !> user has to supply their own read time routine.

   days    = floor(model_time)
   seconds = (model_time-real(days,digits12))*86400
   read_model_time = set_time(seconds, days)

else if ( dart_calendar == 'GREGORIAN' ) then

   ios = nf90_get_att(ncid, VarID, 'units', unitstring)
   call nc_check(ios, 'read_model_time', 'get_att time units "'//trim(filename)//'"')

   if (unitstring(1:10) == 'days since') then

      read(unitstring,'(11x,i4,5(1x,i2))',iostat=ios)year,month,day,hour,minute,second
      if (ios /= 0) then
         write(string1,*)'Unable to interpret time unit attribute. Error status was ',ios
         write(string2,*)'expected "days since YYYY-MM-DD HH:MM:SS", got "'//trim(unitstring)//'"'
         call error_handler(E_ERR, 'read_model_time:', string1, &
                source, revision, revdate, text2=string2, text3=string3)
      endif

      ! This is the start of their calendar
      base_time = set_date(year, month, day, hour, minute, second)

      days    = floor(model_time)
      seconds = (model_time-real(days,digits12))*86400

      delta_time = set_time(seconds, days)
      read_model_time = base_time + delta_time

   else if (unitstring(1:13) == 'seconds since') then

      read(unitstring,'(14x,i4,5(1x,i2))',iostat=ios)year,month,day,hour,minute,second
      if (ios /= 0) then
         write(string1,*)'Unable to interpret time unit attribute. Error status was ',ios
         write(string2,*)'expected "seconds since YYYY-MM-DD HH:MM:SS", got "'//trim(unitstring)//'"'
         call error_handler(E_ERR, 'read_model_time:', string1, &
                source, revision, revdate, text2=string2, text3=string3)
      endif

      ! This is the start of their calendar
      base_time  = set_date(year, month, day, hour, minute, second)
      delta_time = set_time(seconds)

      read_model_time = base_time + delta_time

   else

      write(string1, *) 'looking for "days since" or "seconds since" in the "units" attribute'
      call error_handler(E_ERR, 'read_model_time:', 'unable to set base time for gregorian calendar', &
                         source, revision, revdate, text2=string1, text3=string3)

   endif
else
   call error_handler(E_ERR, 'read_model_time:', &
    'calendar type "'//trim(dart_calendar)//'" unsupported by default read_model_time() routine', &
                      source, revision, revdate, text2=string3)
endif

!>@todo FIXME: do we really want this to print from any
!> task without being asked?  i vote no.
call print_time(read_model_time,'read_model_time')

!>@todo FIXME:
! make print_date() return without error if calendar is no_calendar,
! and then add a call to print_date() here.  (also vote no.)

call nc_close_file(ncid, routine)

end function read_model_time

!--------------------------------------------------------------------
!> Write time to a netcdf file.
!> This routine assumes the file comes in already opened, and doesn't close it.

subroutine write_model_time(ncid, dart_time)

integer,             intent(in) :: ncid
type(time_type),     intent(in) :: dart_time

character(len=*), parameter :: routine = 'write_model_time'
integer  :: ios
integer  :: numdims, ntimes
integer  :: VarID
integer  :: dart_days, dart_seconds
integer  :: unlimitedDimId

logical :: has_unlimited, time_is_unlimited

real(digits12) :: model_time

integer, dimension(NF90_MAX_VAR_DIMS) :: dimIDs
character     (len=NF90_MAX_NAME)     :: dart_calendar, var_calendar

! If there is no unlimited dimension, unlimitedDimID = -1
ios = nf90_inquire(ncid, unlimitedDimId=unlimitedDimId )
call nc_check(ios,routine,'checking unlimited dimension')

has_unlimited = (unlimitedDimID /= -1)

! this is used in many error messages below.  set it here, and
! don't reuse string3 here, please.
write(string3,*)'You may need to supply a model-specific "write_model_time()" to write the time.'

! see what kind of calendar dart is currently running with.
call get_calendar_string(dart_calendar)

ios = nf90_inq_varid(ncid, "time", VarID)

! if the file doesn't already have a "time" variable, we make one
if (ios /= NF90_NOERR) then

   call error_handler(E_MSG, routine, 'no variable "time" found in file', &
              source, revision, revdate, text2='creating one')

   ! begin define mode
   ios = nf90_Redef(ncid)
   call nc_check(ios, routine, "redef")

   ! check to see if there is a time dimension
   ! if it does not exist create it
   ios = nf90_inq_dimid(ncid, "time", dimIDs(1))
   if (ios /= NF90_NOERR) then

      ! If there is already an unlimited dimension, just make a
      ! time dimension of 'normal' size. If there is no unlimited dim already
      ! make the time variable 'unlimited'.

      if (has_unlimited) then
         ios = nf90_def_dim(ncid, "time", 1, dimIDs(1))
         call nc_check(ios, routine, 'def_dim singleton dimension time')
      else
         ios = nf90_def_dim(ncid, "time", nf90_unlimited, dimIDs(1))
         call nc_check(ios, routine, 'def_dim unlimited dimension time')
         has_unlimited = .true.
         unlimitedDimID = dimIDs(1)
      endif

   endif

   ! make the time variable be dimensioned time(time) which is the
   ! netCDF convention for coordinate variables (variables with the
   ! same name as a dimension).
   ios = nf90_def_var(ncid, name="time", xtype=nf90_double, dimids=dimIDs(1), varid=VarID)
   call nc_check(ios, routine, "time def_var")

   ! define time attributes consistent with CF convention
   ios = nf90_put_att(ncid, VarID, "long_name", "valid time of the model state")
   call nc_check(ios, routine, "time long_name")

   if (dart_calendar == 'NO_CALENDAR') then
      ios = nf90_put_att(ncid, VarID, "calendar", "none")
      call nc_check(ios, routine, "calendar long_name")

      ! ncview (actually, probably udunits2) crashes or errors out or
      ! displays misleading plot axes if you use 'days since ...' as the units.
      ! if you simply use 'days' it works much better.

      ios = nf90_put_att(ncid, VarID, "units", "days")
      call nc_check(ios, routine, "units long_name")

   else if (dart_calendar == 'GREGORIAN') then
      ios = nf90_put_att(ncid, VarID, "calendar", "gregorian")
      call nc_check(ios, routine, "calendar long_name")

      ios = nf90_put_att(ncid, VarID, "units", "days since 1601-01-01 00:00:00")
      call nc_check(ios, routine, "units long_name")
   else
      write(string1,*) 'calendar type "'//trim(dart_calendar)// &
                       &'" unsupported by default write_model_time() routine'
      call error_handler(E_ERR, routine, string1, &
                 source, revision, revdate, text2=string3)
   endif

   ! end define mode
   call nc_check( nf90_Enddef(ncid),routine, "Enddef" )
endif

! See if the existing time variable has a calendar and start date to consider

ios = nf90_get_att(ncid, VarID, 'calendar', var_calendar)
if (ios /= NF90_NOERR) var_calendar = 'NO_CALENDAR'
call to_upper(var_calendar)

if (dart_calendar /= var_calendar ) then
   ! allow NO_CALENDAR, NO CALENDAR, and NONE to be synonyms.
   ! replace var_calendar with what dart uses to simplify the tests below.
   if (dart_calendar == 'NO_CALENDAR' .and. (var_calendar == 'NONE'        .or. &
                                             var_calendar == 'NO CALENDAR')) then
      var_calendar = 'NO_CALENDAR'
   else
      write(string1,*)'inconsistent calendar types between DART program and input file.'
      write(string2,*)'DART initialized with: ', trim(dart_calendar), ' File uses: ', trim(var_calendar)
      call error_handler(E_ERR, routine, string1, source, revision, revdate, &
                      text2=string2, text3=string3)
   endif
endif

! convert time to something that netCDF can store, fractional days
call get_time(dart_time, dart_seconds, dart_days)
model_time = real(dart_days,digits12) + real(dart_seconds,digits12)/86400.0_digits12

! need to know how long the time variable is and hammer the last time
ios = nf90_inquire_variable(ncid, VarID, dimids=dimIDs, ndims=numdims)
call nc_check(ios, routine, 'inquire number of dimensions of variable "time"')

if (numdims == 0) then ! variable is a scalar
   ios = nf90_put_var(ncid, VarID, model_time)
   call nc_check(ios, routine, 'put_var scalar "model_time"')
   return
endif

if (numdims > 1) then
   write(string1,*)'Expecting the "time" variable to be a single dimension.'
   call error_handler(E_ERR, routine, string1, &
              source, revision, revdate, text2=string3)
endif

ios = nf90_inquire_dimension(ncid, dimIDs(1), len=ntimes)
call nc_check(ios, routine, 'inquire_dimension for time dimension')

if (dimIDs(1) == unlimitedDimID) time_is_unlimited = .true.

if (ntimes == 0 .and. time_is_unlimited) then
   ntimes = ntimes + 1
elseif (ntimes == 0) then
   write(string1,*)'"time" variable has length 0 but is not the unlimited dimension.'
   call error_handler(E_ERR, routine, string1, &
              source, revision, revdate, text2=string3)
endif

! write dart days and seconds files to netcdf file
ios = nf90_put_var(ncid, VarID, model_time, start=(/ ntimes /))
call nc_check(ios, routine, "put_var model_time")

end subroutine write_model_time

!--------------------------------------------------------------------
!> @}
end module dart_time_io_mod

! <next few lines under version control, do not edit>
! $URL$
! $Id$
! $Revision$
! $Date$

max_stars_repo_name: hkershaw-brown/feature-preprocess

max_issues_repo_name: hkershaw-brown/feature-preprocess

max_forks_repo_name: hkershaw-brown/feature-preprocess


text:
      SUBROUTINE EVOLVE(IPAIR,LEVEL)
*
*
*       Binary diagnostics.
*       -------------------
*
      INCLUDE 'common6.h'
*     COMMON/GAMDOT/  DGAM
      REAL*8  XCM(3),VCM(3),WORK(10)
      INTEGER  IFLAG(KMAX),IOUT(LMAX)
*
*
      ZKT = -BE(1)/FLOAT(NZERO)
      M = 0
      IF (LEVEL.EQ.0) GO TO 500
*     NBSTAT = NBSTAT + 1
*       Prepare diagnostic output for regularized binary.
*
      DO 10 IP = 1,NPAIRS
          IFLAG(IP) = 0
   10 CONTINUE
      JPAIR = IPAIR
*
   20 CALL RESOLV(JPAIR,1)
      ICM = N + JPAIR
      J2 = 2*JPAIR
      J1 = J2 - 1
      SEMI = -0.5*BODY(ICM)/H(JPAIR)
      ZN = SQRT(BODY(ICM)/ABS(SEMI)**3)
      RDOT = (X(1,J1) - X(1,J2))*(XDOT(1,J1) - XDOT(1,J2))
     &     + (X(2,J1) - X(2,J2))*(XDOT(2,J1) - XDOT(2,J2))
     &     + (X(3,J1) - X(3,J2))*(XDOT(3,J1) - XDOT(3,J2))
      ECC = (1.0 - R(JPAIR)/SEMI)**2 + RDOT**2/(SEMI*BODY(ICM))
      ECC = SQRT(ECC)
      APO = ABS(SEMI)*(1.0 + ECC)
      RI = SQRT((X(1,ICM) - RDENS(1))**2 + (X(2,ICM) - RDENS(2))**2 +
     &                                     (X(3,ICM) - RDENS(3))**2)
      ZK = 0.5*BODY(ICM)*(XDOT(1,ICM)**2+XDOT(2,ICM)**2+XDOT(3,ICM)**2)
      ZK = ZK/ZKT
      EB = BODY(J1)*BODY(J2)*H(JPAIR)/BODY(ICM)
      EB = -EB/ZKT
      TK = TWOPI/ZN
      NK = TCR0/TK
*
*       Print time interval since last binary output and set new time.
      DT = 1.0E+06
      K = KZ(4)
      DO 40 L = 1,K
          IF (TIME - TLASTB(L).LT.DT) DT = TIME - TLASTB(L)
   40 CONTINUE
      TLASTB(LEVEL) = TIME
*
      IF (LEVEL.LT.0) DGAM = 1.
      GS = SIGN(GAMMA(JPAIR),DGAM)
      IF (NAME(ICM).GT.0) THEN
          WRITE (4,50)  LEVEL, TIME, NAME(J1), NAME(J2), EB, SEMI, ECC,
     &                  NK, R(JPAIR), GS, LIST(1,J1), RI, ZK, DT, TK, M
      ELSE
          WRITE (4,50)  LEVEL, TIME, NZERO + NAME(J1), NAME(J2), EB,
     &                  SEMI, ECC, NK, R(JPAIR), GS, LIST(1,J1), RI, ZK,
     &                  DT, TK, M, NAME(J1), -(NAME(ICM) + NZERO)
      END IF
*
   50 FORMAT (I3,F8.2,I6,I5,F7.1,F9.5,F6.2,I7,F9.5,F6.2,
     &                                         I5,F7.2,F6.1,2F10.6,5I5)
      M = M + 1
      IFLAG(JPAIR) = -1
*
      NNB = LIST(1,J1)
      DO 60 L = 1,NNB
          IOUT(L) = 0
   60 CONTINUE
      RJMIN2 = 1.E+06
*
      NNB1 = NNB + 1
      DO 150 L = 2,NNB1
          J = LIST(L,J1)
          RIJ2 = (X(1,ICM) - X(1,J))**2 + (X(2,ICM) - X(2,J))**2
     &                                  + (X(3,ICM) - X(3,J))**2
          VIJ2 = (XDOT(1,ICM) - XDOT(1,J))**2 +
     &           (XDOT(2,ICM) - XDOT(2,J))**2 +
     &           (XDOT(3,ICM) - XDOT(3,J))**2
          BB = BODY(ICM) + BODY(J)
          BBI = 1.0/BB
          RIJ = SQRT(RIJ2)
          EB = 0.5*VIJ2 - BB/RIJ
          SEMIJ = -0.5*BB/EB
          RDOT = (X(1,ICM) - X(1,J))*(XDOT(1,ICM) - XDOT(1,J)) +
     &           (X(2,ICM) - X(2,J))*(XDOT(2,ICM) - XDOT(2,J)) +
     &           (X(3,ICM) - X(3,J))*(XDOT(3,ICM) - XDOT(3,J))
          ECCJ = (1.0 - RIJ/SEMIJ)**2 + RDOT**2/(SEMIJ*BB)
          ECCJ = SQRT(ECCJ)
          P = SEMIJ*(1.0 - ECCJ)
          RDOT = RDOT/RIJ
          EB = -BODY(J)*BODY(ICM)*EB*BBI/ZKT
          GEFF = 2.0*BODY(J)*(ABS(SEMI)/RIJ)**3/BODY(ICM)
          IF (GEFF.GT.99.9) GEFF = 99.9
*       Exclude apocentres > A with A = 2/N for hard binary.
          APO = MIN(APO,2.0D0/FLOAT(N))
*
*       Print small impact parameter & significant effective perturbation.
          IF (P.LT.3.*APO.AND.GEFF.GT.GPRINT(1)) THEN
              WRITE (4,100)  NAME(J), EB, SEMIJ, ECCJ, RIJ, GEFF, P,RDOT
  100         FORMAT (16X,I6,F7.1,F9.5,F6.2,F16.5,F6.2,F9.5,F7.2)
              IOUT(L-1) = 1
              IF (J.GT.N.AND.IFLAG(J-N).EQ.0) IFLAG(J-N) = 1
          END IF
*
*       Save parameters for the nearest perturber.
          IF (RIJ2.LT.RJMIN2) THEN
              RJMIN2 = RIJ2
              LMIN = L
              WORK(1) = EB
              WORK(2) = SEMIJ
              WORK(3) = ECCJ
              WORK(4) = RIJ
              WORK(5) = GEFF
              WORK(6) = P
              WORK(7) = RDOT
          END IF
  150 CONTINUE
*
      IF (IOUT(LMIN-1).EQ.0) THEN
          JMIN = LIST(LMIN,J1)
          WRITE (4,100)  NAME(JMIN), (WORK(K),K=1,7)
          IF (JMIN.GT.N.AND.IFLAG(JMIN-N).EQ.0) IFLAG(JMIN-N) = 1
      END IF
*
      DO 190 IP = 1,NPAIRS
          IF(IFLAG(IP).GT.0) THEN
              JPAIR = IP
*       Only include full output for active binaries.
              IF (GAMMA(JPAIR).GT.GMAX) GO TO 20
          END IF
  190 CONTINUE
*
*       Set output times to facilitate DTCRIT estimate in routine KSINT.
      KK = KZ(4)
      DO 200 K = 1,KK
          TLASTB(K) = TIME
  200 CONTINUE
      GO TO 900
*
*       Search for binaries of single particles & c.m. bodies.
  500 SIMAX = 0.01*TCR0
      IF (IPAIR.EQ.0) THEN
          ILOW = IFIRST
          IHIGH = NTOT
      ELSE
          ILOW = IPAIR
          IHIGH = IPAIR
      END IF
      NAM = IPAIR
      IF (NAM.GT.0) NAM = NAME(NAM)
*
      DO 600 I = ILOW,IHIGH
          IF (STEP(I).GT.SIMAX) GO TO 600
          J2 = 0
          RIJ = 1.0E+06
          NNB1 = LIST(1,I) + 1
*
          DO 250 L = 2,NNB1
              J = LIST(L,I)
              IF (STEP(J).GT.SIMAX) GO TO 250
              RIJ2 = (X(1,I) - X(1,J))**2 + (X(2,I) - X(2,J))**2 +
     &                                      (X(3,I) - X(3,J))**2
              IF (RIJ2.GT.RIJ) GO TO 250
              RIJ = RIJ2
              J2 = J
  250     CONTINUE
*
          IF (J2.LT.I.AND.IPAIR.EQ.0) GO TO 600
          RIJ = SQRT(RIJ)
          VIJ2 = (XDOT(1,I) - XDOT(1,J2))**2 +
     &           (XDOT(2,I) - XDOT(2,J2))**2 +
     &           (XDOT(3,I) - XDOT(3,J2))**2
          BODYIJ = BODY(I) + BODY(J2)
          BBI = 1.0/BODYIJ
          EREL = 0.5*VIJ2 - BODYIJ/RIJ
          IF (EREL.GT.-0.1*ECLOSE) GO TO 600
          SEMI = -0.5*BODYIJ/EREL
          ZN = SQRT(BODYIJ/SEMI**3)
          RDOT = (X(1,I) - X(1,J2))*(XDOT(1,I) - XDOT(1,J2)) +
     &           (X(2,I) - X(2,J2))*(XDOT(2,I) - XDOT(2,J2)) +
     &           (X(3,I) - X(3,J2))*(XDOT(3,I) - XDOT(3,J2))
          ECC = (1.0 - RIJ/SEMI)**2 + RDOT**2/(SEMI*BODYIJ)
          ECC = SQRT(ECC)
*
          DO 520 K = 1,3
              XCM(K) = (BODY(I)*X(K,I) + BODY(J2)*X(K,J2))*BBI
              VCM(K) = (BODY(I)*XDOT(K,I) + BODY(J2)*XDOT(K,J2))*BBI
  520     CONTINUE
*
          RI = SQRT((XCM(1) - RDENS(1))**2 + (XCM(2) - RDENS(2))**2 +
     &                                       (XCM(3) - RDENS(3))**2)
          ZK = 0.5*BODYIJ*(VCM(1)**2 + VCM(2)**2 + VCM(3)**2)
          ZK = ZK/ZKT
          EB = BODY(I)*BODY(J2)*EREL*BBI
          EB = -EB/ZKT
          TK = TWOPI/ZN
          NK = TCR0/TK
          GB = 0.0
*
*       Estimate the relative perturbation at apocentre.
          DT = T0(I) - T0(J2)
          DO 540 K = 1,3
              GB = GB + (F(K,I) - (F(K,J2) + 3.0*FDOT(K,J2)*DT))**2
  540     CONTINUE
          GB = 2.0*SQRT(GB)*RIJ**2*BBI
          GB = GB*(SEMI*(1.0 + ECC)/RIJ)**3
          IF (GB.GT.99.9) GB = 99.9
*
*       Set time interval since last binary output.
          DT = 1.0E+06
          K = KZ(4)
          DO 550 L = 1,K
              IF (TIME - TLASTB(L).LT.DT) DT = TIME - TLASTB(L)
  550     CONTINUE
*
*       Only output soft binaries for dominant motion at apocentre.
          IF (EB.LT.1.0.AND.GB.GT.0.5) GO TO 600
*
          IF (I.LE.N.AND.J2.LE.N) THEN
              WRITE (4,50)  LEVEL, TIME, NAME(I), NAME(J2), EB, SEMI,
     &                      ECC, NK, RIJ, GB, LIST(1,I), RI, ZK, DT, TK,
     &                      NAM
          ELSE IF (I.LE.N.AND.J2.GT.N) THEN
              JJ2 = 2*(J2-N)
              JJ1 = JJ2-1
              WRITE (4,50)  LEVEL, TIME, NAME(I), NAME(J2), EB, SEMI,
     &                      ECC, NK, RIJ, GB, LIST(1,I), RI, ZK, DT, TK,
     &                      NAM, NAME(JJ1), NAME(JJ2)
          ELSE IF (I.GT.N.AND.J2.LE.N) THEN
              II2 = 2*(I-N)
              II1 = II2-1
              WRITE (4,50)  LEVEL, TIME, NAME(I), NAME(J2), EB, SEMI,
     &                      ECC, NK, RIJ, GB, LIST(1,I), RI, ZK, DT, TK,
     &                      NAM, NAME(II1), NAME(II2)
          ELSE
              II2 = 2*(I-N)
              II1 = II2-1
              JJ2 = 2*(J2-N)
              JJ1 = JJ2-1
              WRITE (4,50)  LEVEL, TIME, NAME(I), NAME(J2), EB, SEMI,
     &                      ECC, NK, RIJ, GB, LIST(1,I), RI, ZK, DT, TK,
     &                      NAM, NAME(II1), NAME(II2), NAME(JJ1),
     &                      NAME(JJ2)
          END IF
  600 CONTINUE
*
      IF (IPAIR.EQ.0) THEN
          TLASTS = TIME
      ELSE
          TLASTT = TIME
      END IF
*
  900 RETURN
*
      END