text:
  function energy (evapObj, iz, ia)

! ======================================================================
!
!****************This routine was originally in the LAHET code********
!  ENERGY
!   Calculates mass excess using OLD Cameron formula when outside
!   Wapstra table!!!!!  (Based on O-16 = 0.00 binding energies!)
!
!    Modified to use Moller-Nix Masses for N, Z >= 8; macroscopic
!    M-N masses outside the table, and Wapstra tables only for N,Z
!   < 8 and masses not measured.  AJS
!
!   For light nuclei outside the Wapstra systematics table, relies on
!   the old Cameron calculation.....VERY IMPROBABLE!
!  Z = 1, A > 12;  Z = 2, A > 14;  Z = 3, A > 17;  Z = 4, A > 22;
!  Z = 5, A > 25;  Z = 6, A > 28;  Z = 7, A > 31.
!
! =====================================================================
! <variables>
!     ia   :   mass number of nucleus               (IN)
!     iz   :   charge of nucleus                    (IN)
!  energy  :   Mass excess [MeV]                    (OUT)
!
!   Last change: 13-AUG-2003 by NVMokhov
!   Modified by A. J. Sierk, LANL T-16, October, 2003.
!   Edited by AJS, LANL T-2, December, 2011.
!   Edited by CMJ, XCP-3, July 2018 (Evap. class creation)
!
! ======================================================================

    use, intrinsic :: iso_fortran_env, only: int32, real64
    use evaporationParams, only: zro, ato3rd
    use evaporationFissionData, only: wapsm, inn, iiz, sz, sn

    implicit none
    class(Evaporation), intent(inout) :: evapObj
    integer(int32),     intent(in   ) :: iz
    integer(int32),     intent(in   ) :: ia
    real(real64)                      :: energy

    integer(int32) :: in
    real(real64)   :: a, a3, z
    logical        :: nsmal, zsmal

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

    energy = 1.d10
    in = ia - iz

    ! Limits in program
    ! IA = 4, IZ = 1, IN = 3
    if ( (iz < zro .OR. iiz < iz) .OR. (in < zro .OR. inn < in) ) then
       ! Outside function limits, warn user
       write(evapObj%io%message,1000) ia, iz, in
       call evapObj%io%print(2, 3, evapObj%io%message)
    endif
    if (iz < 0 .or. in < 0) return

    zsmal = iz < 8 .or. in < 8
    if (iz == 0) then
       nsmal = .true.
       if (in == 0) return
    elseif (iz < 8) then
       nsmal = in < evapObj%evapMolnix%nmina(iz) .or. in > evapObj%evapMolnix%nmaxa(iz)
    elseif (iz >= 8) then
       nsmal = in < evapObj%evapMolnix%nmin(iz-7) .or. in > evapObj%evapMolnix%nmax(iz-7)
    endif

    ! Obtain mass excess values
    if (.not.zsmal .or. (zsmal .and. .not.nsmal)) then
       energy = evapObj%evapMolnix%defineEnergy (iz, in, 2)
    elseif (zsmal .and. nsmal) then
       energy = wapsm(iz,in)
       if (energy == 0.d0) then
          if (iz == 6 .and. ia == 12) return
          if (iz > 0 .and. in > 0) then
!   Cameron's mass formula for nuclides outside the Wapstra table:
!   ONLY for light nuclei outside the Wapstra table in mass.tbl.
             a = dble(ia)
             z = dble(iz)
             a3 = ato3rd(ia)
             energy = evapObj%cam (a, z, a3) + sz(iz) + sn(in)
          endif
       endif
    endif

    return
! ======================================================================
1000 format("Secondary array exceeded in 'energy.f90', ", &
          & "line 55 (A=", i3, ", Z = ", i3, " [N = ", i3, "]).")
! ======================================================================
  end function energy

  function cam (evapObj, a, z, b)

! ======================================================================
!
!   Cameron mass-excess formula; extracted from inline function from
!   ENERGY. Formula originally in LAHET. b = A**(1/3).
!
!   Written by A. J. Sierk, LANL T-16, September, 2003.
!   Edited by AJS, LANL T-2, December, 2011.
!   Edited by LMK, XCP-3, July 2013 (included error protection).
!   Edited by CMJ, XCP-3, July 2018 (Evap. class creation)
!
! ======================================================================

    use, intrinsic :: iso_fortran_env, only: int32, real64
    use evaporationParams, only: one, two, ato3rd

    implicit none
    class(Evaporation), intent(inout) :: evapObj
    real(real64),       intent(in   ) :: a
    real(real64),       intent(in   ) :: z
    real(real64),       intent(in   ) :: b
    real(real64)                      :: cam

    integer(int32) :: iz
    real(real64)   :: abi, zft, tempa, tempb

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

    iz = nint(z)
    tempa = a
    if ( tempa < div0Lim .and. tempa > -div0Lim) then
       tempa = div0Lim
       write(evapObj%io%message,1000) "134, 139-141"
       call evapObj%io%print(4, 3, evapObj%io%message)
    end if
    abi = one - two*z/tempa
    zft = ato3rd(iz)**4
    tempb = b
    if ( tempb < div0Lim .and. tempb > -div0Lim) then
       tempb = div0Lim
       write(evapObj%io%message,1000) "139-141"
       call evapObj%io%print(4, 3, evapObj%io%message)
    end if
    cam = 8.071323d0*tempa - 0.782354d0*z - 17.0354d0*tempa*(one - &
         & 1.84619d0*abi**2) + 25.8357d0*tempb**2*(one - &
         & 1.71219d0*abi**2)*(one - 0.62025d0/tempb**2)**2 + &
         & 0.779d0*z*(z - one)*(one - 1.5849d0/tempb**2 + &
         & 1.2273d0/tempa + 1.5772d0/(tempa*tempb))/tempb - 0.4323d0 * &
         & zft*(one - 0.57811d0/tempb - &
         & 0.14518d0/tempb**2 + 0.49597d0/tempa)/tempb

    return

! ======================================================================
1000 format("Divide by zero error prevented in 'energy.f90,' line ", A)
! ======================================================================
  end function cam

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:
        SUBROUTINE DSPRCL ( BLOCK )
      INCLUDE 'DSIOF.COM'
      INCLUDE 'XNSTRN.COM'
        INTEGER BLOCK( 15 )
        INTEGER IDIV( 4 )
        DATA    IDIV / 1, 2, 1, 2 /
        BLOCK(  2 ) = IAND( IBASE( INDCBP ), MASKQ4 )
        BLOCK(  3 ) = IAND( IBASE( INDCBP ), MASKQ3 )
        BLOCK(  3 ) = BLOCK( 3 ) / MULQ3
        BLOCK( 11 ) = NWRDEL( BLOCK( 2 ) )
        BLOCK( 12 ) = IBASE( INDCBP+1 )
        BLOCK( 14 ) = IDIV( BLOCK( 2 ) )
        RETURN
        END

ext: f


text:
      PROGRAM xjacobi
C     driver for routine jacobi
      INTEGER NP,NMAT
      PARAMETER(NP=10,NMAT=3)
      INTEGER i,ii,j,jj,k,kk,l,ll,nrot,num(3)
      REAL ratio,d(NP),v(NP,NP),r(NP)
      REAL a(3,3),b(5,5),c(10,10),e(NP,NP)
      DATA num/3,5,10/
      DATA a/1.0,2.0,3.0,2.0,2.0,3.0,3.0,3.0,3.0/
      DATA b/-2.0,-1.0,0.0,1.0,2.0,-1.0,-1.0,0.0,1.0,2.0,
     *     0.0,0.0,0.0,1.0,2.0,1.0,1.0,1.0,1.0,2.0,
     *     2.0,2.0,2.0,2.0,2.0/
      DATA c/5.0,4.3,3.0,2.0,1.0,0.0,-1.0,-2.0,-3.0,-4.0,
     *     4.3,5.1,4.0,3.0,2.0,1.0,0.0,-1.0,-2.0,-3.0,
     *     3.0,4.0,5.0,4.0,3.0,2.0,1.0,0.0,-1.0,-2.0,
     *     2.0,3.0,4.0,5.0,4.0,3.0,2.0,1.0,0.0,-1.0,
     *     1.0,2.0,3.0,4.0,5.0,4.0,3.0,2.0,1.0,0.0,
     *     0.0,1.0,2.0,3.0,4.0,5.0,4.0,3.0,2.0,1.0,
     *     -1.0,0.0,1.0,2.0,3.0,4.0,5.0,4.0,3.0,2.0,
     *     -2.0,-1.0,0.0,1.0,2.0,3.0,4.0,5.0,4.0,3.0,
     *     -3.0,-2.0,-1.0,0.0,1.0,2.0,3.0,4.0,5.0,4.0,
     *     -4.0,-3.0,-2.0,-1.0,0.0,1.0,2.0,3.0,4.0,5.0/
      do 24 i=1,NMAT
        if (i.eq.1) then
          do 12 ii=1,3
            do 11 jj=1,3
              e(ii,jj)=a(ii,jj)
11          continue
12        continue
          call jacobi(e,3,NP,d,v,nrot)
        else if (i.eq.2) then
          do 14 ii=1,5
            do 13 jj=1,5
              e(ii,jj)=b(ii,jj)
13          continue
14        continue
          call jacobi(e,5,NP,d,v,nrot)
        else if (i.eq.3) then
          do 16 ii=1,10
            do 15 jj=1,10
              e(ii,jj)=c(ii,jj)
15          continue
16        continue
          call jacobi(e,10,NP,d,v,nrot)
        endif
        write(*,'(/1x,a,i2)') 'Matrix Number ',i
        write(*,'(1x,a,i3)') 'Number of JACOBI rotations: ',nrot
        write(*,'(/1x,a)') 'Eigenvalues:'
        do 17 j=1,num(i)
          write(*,'(1x,5f12.6)') d(j)
17      continue
        write(*,'(/1x,a)') 'Eigenvectors:'
        do 18 j=1,num(i)
          write(*,'(1x,t5,a,i3)') 'Number',j
          write(*,'(1x,5f12.6)') (v(k,j),k=1,num(i))
18      continue
C     eigenvector test
        write(*,'(/1x,a)') 'Eigenvector Test'
        do 23 j=1,num(i)
          do 21 l=1,num(i)
            r(l)=0.0
            do 19 k=1,num(i)
              if (k.gt.l) then
                kk=l
                ll=k
              else
                kk=k
                ll=l
              endif
              if (i.eq.1) then
                r(l)=r(l)+a(ll,kk)*v(k,j)
              else if (i.eq.2) then
                r(l)=r(l)+b(ll,kk)*v(k,j)
              else if (i.eq.3) then
                r(l)=r(l)+c(ll,kk)*v(k,j)
              endif
19          continue
21        continue
          write(*,'(/1x,a,i3)') 'Vector Number',j
          write(*,'(/1x,t7,a,t18,a,t31,a)')
     *         'Vector','Mtrx*Vec.','Ratio'
          do 22 l=1,num(i)
            ratio=r(l)/v(l,j)
            write(*,'(1x,3f12.6)') v(l,j),r(l),ratio
22        continue
23      continue
        write(*,*) 'press RETURN to continue...'
        read(*,*)
24    continue
      END

ext: for

max_stars_repo_path: NR-Functions/Numerical Recipes- Example & Functions/Examples/xjacobi.for


text:
!##############################################################################
!# ****************************************************************************
!# <name> kktsystemhierarchy </name>
!# ****************************************************************************
!#
!# <purpose>
!# Realises a hierarchy of solutions of a KKT system.
!# </purpose>
!##############################################################################

module kktsystemhierarchy

  use fsystem
  use genoutput
  use spatialdiscretisation
  use timediscretisation

  use structuresdiscretisation
  use structuresoptcontrol
  use structuresgeneral

  use fespacehierarchybase
  use spacetimehierarchy

  use kktsystemspaces
  use kktsystem

  use structuresoperatorasm
  use structuresboundaryconditions
  use spacematvecassembly

  implicit none

  private

!<types>

!<typeblock>

  ! A hierarchy of solutions of a KKT system
  type t_kktsystemHierarchy

    ! A space-time hierarchy based on the primal space
    type(t_spaceTimeHierarchy), pointer :: p_rspaceTimeHierPrimal => null()

    ! A space-time hierarchy based on the dual space
    type(t_spaceTimeHierarchy), pointer :: p_rspaceTimeHierDual => null()

    ! A space-time hierarchy based on the control space
    type(t_spaceTimeHierarchy), pointer :: p_rspaceTimeHierControl => null()

    ! Hierarchy of space-time operator assembly structures for the discretisation
    type(t_spacetimeOpAsmHierarchy), pointer :: p_roperatorAsmHier => null()

    ! Number of levels in the hierarchy
    integer :: nlevels = 0

    ! Array of KKT system (solutions) for all the levels
    type(t_kktsystem), dimension(:), pointer :: p_RkktSystems => null()

    ! Alternative solution on the maximum level or NULL if there is none.
    type(t_kktsystem), pointer :: p_rkktSystemMaxLevel => null()

  end type

!</typeblock>

  public :: t_kktsystemHierarchy

!<typeblock>

  ! A hierarchy of directional derivatives of a KKT system
  type t_kktsystemDirDerivHierarchy

    ! A connected KKT system hierarchy specifying the evaluation
    ! point of the derivative.
    type(t_kktsystemHierarchy), pointer :: p_rkktsystemHierarchy => null()

    ! Number of levels in the hierarchy
    integer :: nlevels = 0

    ! Array of KKT system (solutions) for all the levels
    type(t_kktsystemDirDeriv), dimension(:), pointer :: p_RkktSysDirDeriv => null()

  end type

!</typeblock>

  public :: t_kktsystemDirDerivHierarchy

!</types>

  ! Initialises a primal solution according to level ilevel in the
  ! KKT hierarchy.
  public :: kkth_initPrimalSol

  ! Initialises a primal solution according to level ilevel in the
  ! KKT hierarchy.
  public :: kkth_initDualSol

  ! Initialises a primal solution according to level ilevel in the
  ! KKT hierarchy.
  public :: kkth_initControl

  ! Initialises a KKT system structure for level ilevel of the hierarchy
  public :: kkth_initKKTSystem

  ! Initialises a hierarchy of KKT systems.
  public :: kkth_initHierarchy

  ! Cleans up a hierarchy of KKT systems.
  public :: kkth_doneHierarchy

  ! Get the KKT system at a level
  public :: kkth_getKKTsystem

  ! Initialises a KKT system derivative structure for level ilevel of the hierarchy
  public :: kkth_initKKTSystemDirDeriv

  ! Initialises a hierarchy of KKT system directional derivatives
  public :: kkth_initHierarchyDirDeriv

  ! Cleans up a hierarchy of KKT systems directional derivatives
  public :: kkth_doneHierarchyDirDeriv

  ! Get the KKT system directional derivative at a level
  public :: kkth_getKKTsystemDirDeriv

contains

  ! ***************************************************************************

!<subroutine>

  subroutine kkth_initPrimalSol (rvector,rkktsystemHierarchy,ilevel)

!<description>
  ! Initialises a primal solution according to level ilevel in the
  ! KKT hierarchy.
!</description>

!<input>
  ! Hierarchy of KKT system solutions
  type(t_kktsystemHierarchy), intent(in) :: rkktsystemHierarchy

  ! Level, the KKT system should be created corresponding to.
  ! If this is set to <= 0, rkktsystem will be created at level NLMAX+ilevel
  integer, intent(in) :: ilevel
!</input>

!<output>
  ! Vector to be created
  type(t_primalSpace), intent(out) :: rvector
!</output>

!</subroutine>

    ! local variables
    integer :: ilev
    type(t_feSpaceLevel), pointer :: p_rspaceDiscr
    type(t_timeDiscretisation), pointer :: p_rtimeDiscr

    ilev = ilevel
    if (ilevel .le. 0) ilev = rkktsystemHierarchy%nlevels+ilevel

    ! Get the spatial and time discretisations of that level
    call sth_getLevel (rkktsystemHierarchy%p_rspaceTimeHierPrimal,ilev,&
        p_rspaceDiscr,p_rtimeDiscr)

    ! Create the vector
    call kktsp_initPrimalVector (rvector,&
        p_rspaceDiscr%p_rdiscretisation,p_rtimeDiscr)

  end subroutine

  ! ***************************************************************************

!<subroutine>

  subroutine kkth_initDualSol (rvector,rkktsystemHierarchy,ilevel)

!<description>
  ! Initialises a primal solution according to level ilevel in the
  ! KKT hierarchy.
!</description>

!<input>
  ! Hierarchy of KKT system solutions
  type(t_kktsystemHierarchy), intent(in) :: rkktsystemHierarchy

  ! Level, the KKT system should be created corresponding to.
  ! If this is set to <= 0, rkktsystem will be created at level NLMAX+ilevel
  integer, intent(in) :: ilevel
!</input>

!<output>
  ! Vector to be created
  type(t_dualSpace), intent(out) :: rvector
!</output>

!</subroutine>

    ! local variables
    integer :: ilev
    type(t_feSpaceLevel), pointer :: p_rspaceDiscr
    type(t_timeDiscretisation), pointer :: p_rtimeDiscr

    ilev = ilevel
    if (ilevel .le. 0) ilev = rkktsystemHierarchy%nlevels+ilevel

    ! Get the spatial and time discretisations of that level
    call sth_getLevel (rkktsystemHierarchy%p_rspaceTimeHierDual,ilev,&
        p_rspaceDiscr,p_rtimeDiscr)

    ! Create the vector
    call kktsp_initDualVector (rvector,&
        p_rspaceDiscr%p_rdiscretisation,p_rtimeDiscr)

  end subroutine

  ! ***************************************************************************

!<subroutine>

  subroutine kkth_initControl (rvector,rkktsystemHierarchy,ilevel)

!<description>
  ! Initialises a primal solution according to level ilevel in the
  ! KKT hierarchy.
!</description>

!<input>
  ! Hierarchy of KKT system solutions
  type(t_kktsystemHierarchy), intent(in) :: rkktsystemHierarchy

  ! Level, the KKT system should be created corresponding to.
  ! If this is set to <= 0, rkktsystem will be created at level NLMAX+ilevel
  integer, intent(in) :: ilevel
!</input>

!<output>
  ! Vector to be created
  type(t_controlSpace), intent(out) :: rvector
!</output>

!</subroutine>

    ! local variables
    integer :: ilev
    type(t_feSpaceLevel), pointer :: p_rspaceDiscr
    type(t_timeDiscretisation), pointer :: p_rtimeDiscr

    ilev = ilevel
    if (ilevel .le. 0) ilev = rkktsystemHierarchy%nlevels+ilevel

    ! Get the spatial and time discretisations of that level
    call sth_getLevel (rkktsystemHierarchy%p_rspaceTimeHierControl,ilev,&
        p_rspaceDiscr,p_rtimeDiscr)

    ! Create the vector
    call kktsp_initControlVector (rvector,&
        p_rspaceDiscr%p_rdiscretisation,p_rtimeDiscr)

  end subroutine

  ! ***************************************************************************

!<subroutine>

  subroutine kkth_initKKTSystem (rkktsystem,rkktsystemHierarchy,ilevel,roperatorAsmHier,roptcBDC)

!<description>
  ! Initialises a KKT system structure for level ilevel of the hierarchy
!</description>

!<input>
  ! Hierarchy of KKT system solutions
  type(t_kktsystemHierarchy), intent(in) :: rkktsystemHierarchy

  ! Level, the KKT system should be created corresponding to.
  ! If this is set to <= 0, rkktsystem will be created at level NLMAX+ilevel
  integer, intent(in) :: ilevel

  ! Hierarchy of all possible space-time operators
  type(t_spacetimeOpAsmHierarchy), intent(in), target :: roperatorAsmHier

  ! Boudary conditions of the KKT system
  type(t_optcBDC), intent(in), target :: roptcBDC
!</input>

!<output>
  ! KKT sytem to be created
  type(t_kktsystem), intent(out) :: rkktsystem
!</output>

!</subroutine>

    ! local variables
    integer :: ilev, ispacelevel, itimelevel

    ilev = ilevel
    if (ilevel .le. 0) ilev = rkktsystemHierarchy%nlevels+ilevel

    ! Get the spatial and time discretisations of that level
    call sth_getLevel (rkktsystemHierarchy%p_rspaceTimeHierPrimal,ilev,&
        ispacelevel=ispacelevel, itimelevel=itimelevel)

    ! Initialise the KKT system on level i
    call kkt_initKKTsystem (rkktsystem,roperatorAsmHier,&
        ispacelevel, itimelevel, roptcBDC)

  end subroutine

  ! ***************************************************************************

!<subroutine>

  subroutine kkth_initHierarchy (rkktsystemHierarchy,roperatorAsmHier,&
      rspaceTimeHierPrimal,rspaceTimeHierDual,rspaceTimeHierControl,ballocate,&
      roptcBDC,rkktSystemMaxLevel)

!<description>
  ! Initialises a hierarchy of KKT systems.
!</description>

!<input>
  ! Hierarchy of space-time operator assembly structures.
  type(t_spacetimeOpAsmHierarchy), intent(in), target :: roperatorAsmHier

  ! Space-time hierarchy in the primal space
  type(t_spacetimeHierarchy), intent(in), target :: rspaceTimeHierPrimal

  ! Space-time hierarchy in the dual space
  type(t_spacetimeHierarchy), intent(in), target :: rspaceTimeHierDual

  ! Space-time hierarchy in the control space
  type(t_spacetimeHierarchy), intent(in), target :: rspaceTimeHierControl

  ! Whether or not to allocate memory for the solutions
  logical, intent(in) :: ballocate

  ! Boudary conditions of the KKT system
  type(t_optcBDC), intent(in), target :: roptcBDC

  ! OPTIONAL: A KKT system solution on the maximum level.
  ! If specified, this solution is used on the maximum level, and
  ! no memory is allocated in the hierarchy for the maximum level.
  type(t_kktSystem), intent(in), optional, target :: rkktSystemMaxLevel
!</input>

!<output>
  ! Hierarchy of KKT system solutions
  type(t_kktsystemHierarchy), intent(out) :: rkktsystemHierarchy
!</output>

!</subroutine>

    ! local variables
    integer :: i,nlevelallocate

    ! Remember the structures
    rkktsystemHierarchy%p_rspaceTimeHierPrimal => rspaceTimeHierPrimal
    rkktsystemHierarchy%p_rspaceTimeHierDual => rspaceTimeHierDual
    rkktsystemHierarchy%p_rspaceTimeHierControl => rspaceTimeHierControl
    rkktsystemHierarchy%p_roperatorAsmHier => roperatorAsmHier

    ! Numer of levels in the hierarchy?
    rkktsystemHierarchy%nlevels = rspaceTimeHierPrimal%nlevels

    ! Alternative solution on the max. level?
    nlevelallocate = rkktsystemHierarchy%nlevels
    if (present(rkktSystemMaxLevel)) then
      rkktsystemHierarchy%p_rkktSystemMaxLevel => rkktSystemMaxLevel
      nlevelallocate = nlevelallocate - 1
    end if

    nullify(rkktsystemHierarchy%p_RkktSystems)

    if (ballocate) then

      ! Create the KKT system solutions
      allocate(rkktsystemHierarchy%p_RkktSystems(nlevelallocate))
      do i=1,nlevelallocate
        ! Initialise the KKT system on level i
        call kkth_initKKTSystem (&
            rkktsystemHierarchy%p_RkktSystems(i),rkktsystemHierarchy,i,&
            roperatorAsmHier,roptcBDC)
      end do

    end if

  end subroutine

  ! ***************************************************************************

!<subroutine>

  subroutine kkth_doneHierarchy (rkktsystemHierarchy)

!<description>
  ! Cleans up a hierarchy of KKT systems.
!</description>

!<inputoutput>
  ! Hierarchy of KKT system solutions to be cleaned up
  type(t_kktsystemHierarchy), intent(inout) :: rkktsystemHierarchy
!</inputoutput>

!</subroutine>

    ! local variables
    integer :: i

    ! Remember the structures
    nullify(rkktsystemHierarchy%p_rspaceTimeHierPrimal)
    nullify(rkktsystemHierarchy%p_rspaceTimeHierDual)
    nullify(rkktsystemHierarchy%p_rspaceTimeHierControl)
    nullify(rkktsystemHierarchy%p_roperatorAsmHier)
    nullify(rkktsystemHierarchy%p_rkktSystemMaxLevel)

    if (associated(rkktsystemHierarchy%p_RkktSystems)) then
      ! Clean up the solutions on every level
      do i=size(rkktsystemHierarchy%p_RkktSystems),1,-1
        call kkt_doneKKTsystem (rkktsystemHierarchy%p_RkktSystems(i))
      end do

      deallocate(rkktsystemHierarchy%p_RkktSystems)
    end if

    rkktsystemHierarchy%nlevels = 0

  end subroutine

  ! ***************************************************************************

!<subroutine>

  subroutine kkth_getKKTsystem (rkktsystemHierarchy,ilevel,p_rkktSystem)

!<description>
  ! Returns a pointer to the KKT system on level ilevel.
!</description>

!<input>
  ! Hierarchy of KKT system solutions.
  type(t_kktsystemHierarchy), intent(in) :: rkktsystemHierarchy

  ! Level to which a pointer should be set up.
  ! If the value is <= 0, the KKT system at level NLMAX+ilevel is returned.
  integer, intent(in) :: ilevel
!</input>

!<inputoutput>
  ! Pointer to the KKT system on level ilevel
  type(t_kktsystem), pointer :: p_rkktSystem
!</inputoutput>

!</subroutine>

    integer :: ilev

    ilev = ilevel
    if (ilev .le. 0) ilev = rkktsystemHierarchy%nlevels+ilev

    ! Get the pointer
    if ((ilev .eq. rkktsystemHierarchy%nlevels) .and. &
        associated(rkktsystemHierarchy%p_rkktSystemMaxLevel)) then
      ! Take the alternative solution on the maximum level
      p_rkktSystem => rkktsystemHierarchy%p_rkktSystemMaxLevel
    else
      p_rkktSystem => rkktsystemHierarchy%p_RkktSystems(ilev)
    end if

  end subroutine

  ! ***************************************************************************

!<subroutine>

  subroutine kkth_initKKTSystemDirDeriv (&
      rkktsystemDirDeriv,rkktsystemDirDerivHier,ilevel)

!<description>
  ! Initialises a KKT system derivative for level ilevel of the hierarchy
!</description>

!<input>
  ! Hierarchy structure for the directional derivative.
  type(t_kktsystemDirDerivHierarchy), intent(inout) :: rkktsystemDirDerivHier

  ! Level, the KKT system derivative should be created corresponding to.
  ! If this is set to <= 0, rkktsystem will be created at level NLMAX+ilevel
  integer, intent(in) :: ilevel
!</input>

!<output>
  ! KKT sytem derivative structure to be created
  type(t_kktsystemDirDeriv), intent(out) :: rkktsystemDirDeriv
!</output>

!</subroutine>

    ! local variables
    integer :: ilev
    type(t_kktSystem), pointer :: p_rkktsystem

    ilev = ilevel
    if (ilevel .le. 0) ilev = rkktsystemDirDerivHier%nlevels+ilevel

    ! Get the KKT system solution on that level
    call kkth_getKKTsystem (&
        rkktsystemDirDerivHier%p_rkktsystemHierarchy,ilevel,p_rkktsystem)

    ! Initialise the KKT system on level i
    call kkt_initKKTsystemDirDeriv (rkktsystemDirDeriv,p_rkktsystem)

  end subroutine

  ! ***************************************************************************

!<subroutine>

  subroutine kkth_initHierarchyDirDeriv (rkktsystemDirDerivHier,rkktsystemHierarchy)

!<description>
  ! Initialises a hierarchy of KKT system directional derivatives
!</description>

!<input>
  ! Evaluation point of the derivative
  type(t_kktsystemHierarchy), intent(in), target :: rkktsystemHierarchy
!</input>

!<output>
  ! Structure for the directional derivative in the point rkktsystemHierarchy.
  type(t_kktsystemDirDerivHierarchy), intent(out) :: rkktsystemDirDerivHier
!</output>

!</subroutine>

    ! local variables
    integer :: i

    rkktsystemDirDerivHier%p_rkktsystemHierarchy => rkktsystemHierarchy
    rkktsystemDirDerivHier%nlevels = rkktsystemHierarchy%nlevels

    ! Allocate memory for the solutions
    allocate(rkktsystemDirDerivHier%p_RkktSysDirDeriv(rkktsystemDirDerivHier%nlevels))
    do i=1,rkktsystemDirDerivHier%nlevels
      ! Initialise the derivative on level i
      call kkth_initKKTSystemDirDeriv (&
          rkktsystemDirDerivHier%p_RkktSysDirDeriv(i),rkktsystemDirDerivHier,i)
    end do

  end subroutine

  ! ***************************************************************************

!<subroutine>

  subroutine kkth_doneHierarchyDirDeriv (rkktsystemDirDerivHier)

!<description>
  ! Cleans up a hierarchy of KKT systems.
!</description>

!<inputoutput>
  ! Hierarchy structure for the directional derivative.
  type(t_kktsystemDirDerivHierarchy), intent(inout) :: rkktsystemDirDerivHier
!</inputoutput>

!</subroutine>

    ! local variables
    integer :: i

    ! Release memory
    if (associated(rkktsystemDirDerivHier%p_RkktSysDirDeriv)) then
      ! Clean up the solutions on every level
      do i=size(rkktsystemDirDerivHier%p_RkktSysDirDeriv),1,-1
        call kkt_doneKKTsystemDirDeriv (rkktsystemDirDerivHier%p_RkktSysDirDeriv(i))
      end do

      deallocate(rkktsystemDirDerivHier%p_RkktSysDirDeriv)
    end if

    ! Final cleanup
    nullify(rkktsystemDirDerivHier%p_rkktsystemHierarchy)
    rkktsystemDirDerivHier%nlevels = 0

  end subroutine

  ! ***************************************************************************

!<subroutine>

  subroutine kkth_getKKTsystemDirDeriv (rkktsystemDirDerivHier,ilevel,p_rkktSystemDirDeriv)

!<description>
  ! Returns a pointer to the KKT system on level ilevel.
!</description>

!<input>
  ! Hierarchy structure for the directional derivative.
  type(t_kktsystemDirDerivHierarchy), intent(inout) :: rkktsystemDirDerivHier

  ! Level to which a pointer should be set up.
  ! If the value is <= 0, the KKT system at level NLMAX+ilevel is returned.
  integer, intent(in) :: ilevel
!</input>

!<inputoutput>
  ! Pointer to the KKT system on level ilevel
  type(t_kktsystemDirDeriv), pointer :: p_rkktSystemDirDeriv
!</inputoutput>

!</subroutine>

    integer :: ilev

    ilev = ilevel
    if (ilev .le. 0) ilev = rkktsystemDirDerivHier%nlevels+ilev

    ! Get the pointer
    p_rkktSystemDirDeriv => rkktsystemDirDerivHier%p_RkktSysDirDeriv(ilev)

  end subroutine

end module


text:
      subroutine test_mpi3rma_acc_replace_lock( base, array_size, DIM1, win, comm )

      use utilities, only: check_fompi_status

#if INCLUDE_TYPE==1
      implicit none

      include 'fompif.h'
#endif
#if INCLUDE_TYPE==2
      use fompi

      implicit none
#endif

      integer, intent(in) :: array_size
      integer, dimension(array_size), intent(inout) :: base
      integer, intent(in) :: DIM1
      integer, intent(in) :: win
      integer, intent(in) :: comm

      integer :: i
      integer :: rank
      integer :: commsize
      integer :: mpierr
      integer(kind=MPI_ADDRESS_KIND) :: start

      logical :: flag
      integer(kind=MPI_ADDRESS_KIND) :: create_flavor

      integer(kind=MPI_ADDRESS_KIND), dimension(:), allocatable :: base_adr

      character(50) :: filename

      write (filename,*) "acc_replace_lock.f90"

      call MPI_Comm_rank( comm, rank, mpierr )
      call MPI_Comm_size( comm, commsize, mpierr )

      call foMPI_Win_lock( foMPI_LOCK_EXCLUSIVE, rank, 0, win, mpierr )
      call check_fompi_status(mpierr, MPI_SUCCESS, filename, 34, rank)

      do i=1,DIM1
        base(i) = rank+i
      enddo
      if (rank .eq. 0) then
        do i=1,DIM1
          base(DIM1+i) = commsize+i
        enddo
      endif

      call foMPI_Win_unlock( rank, win, mpierr )
      call check_fompi_status(mpierr, MPI_SUCCESS, filename, 46, rank)

      call foMPI_Win_get_attr( win, foMPI_WIN_CREATE_FLAVOR, create_flavor, flag, mpierr )
      call check_fompi_status(mpierr, MPI_SUCCESS, filename, 49, rank)
      if (flag .neqv. .true.) then
        write (*,*) "Didn't found foMPI_WIN_CREATE_FLAVOR keyval."
        STOP
      endif

      if ( create_flavor .eq. foMPI_WIN_FLAVOR_DYNAMIC ) then
        allocate( base_adr(commsize), stat=mpierr )
        call foMPI_Get_address( base(1), start, mpierr )
        call MPI_Allgather( start, 1, MPI_AINT, base_adr(1), 1, MPI_AINT, comm, mpierr ) ! the init phase is finished
      else
        start = 0
        call MPI_Barrier( comm, mpierr ) ! the init phase is finished
      endif

      if (rank .eq. 0) then
        do i=0,commsize-1
          call foMPI_Win_lock( foMPI_LOCK_SHARED, i, 0, win, mpierr )
          call check_fompi_status(mpierr, MPI_SUCCESS, filename, 67, rank)

          if ( create_flavor .eq. foMPI_WIN_FLAVOR_DYNAMIC ) then
            call foMPI_Accumulate( base(DIM1+1), DIM1, MPI_INTEGER, i, base_adr(i+1), DIM1, MPI_INTEGER, foMPI_REPLACE, &
              win, mpierr )
          else
            call foMPI_Accumulate( base(DIM1+1), DIM1, MPI_INT, i, start, DIM1, MPI_INTEGER, foMPI_REPLACE, win, mpierr )
          endif
          call check_fompi_status(mpierr, MPI_SUCCESS, filename, 74, rank)

          call foMPI_Win_unlock( i, win, mpierr )
          call check_fompi_status(mpierr, MPI_SUCCESS, filename, 77, rank)
        enddo
      endif

      call MPI_Barrier( comm, mpierr ) !wait for all accesses to finish

      if ( create_flavor .eq. foMPI_WIN_FLAVOR_DYNAMIC ) then
        deallocate( base_adr, stat=mpierr )
      endif

      do i=1,DIM1
        if( base(i) .ne. (commsize+i) ) then
          write (*,*) "Error on rank ", rank, " in ", trim(filename), &
            " line ", 88, ": element ", i, " is ", base(i), &
            " instead of ", commsize+i, "."
          STOP
        endif
      enddo

      call MPI_Barrier( comm, mpierr )

      if (rank .eq. 0) then
        write (*,*) "No Error"
      endif

      end subroutine test_mpi3rma_acc_replace_lock

max_stars_repo_name: spcl/foMPI-NA

max_issues_repo_name: spcl/foMPI-NA

max_forks_repo_name: spcl/foMPI-NA


text:
! { dg-do compile }
! Program to test the initialisation range of enumerators
! and kind values check

program main
  implicit none
  enum, bind (c)
    enumerator :: pp, qq = 4294967295, rr ! { dg-error "too big for its kind" }
  end enum  ! { dg-error "has no ENUMERATORS" }

  enum, bind (c)
    enumerator :: p , q = 4294967299_8, r  ! { dg-error "Arithmetic overflow" }
  end enum  ! { dg-error "has no ENUMERATORS" }

end program main

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text: #redirect Mace Ranch Park

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
      subroutine  dated   ( idatex,jtimex )
!***********************************************************************
!  Copyright, 2004,  The  Regents  of the  University of California.
!  This program was prepared by the Regents of the University of
!  California at Los Alamos National Laboratory (the University) under
!  contract No. W-7405-ENG-36 with the U.S. Department of Energy (DOE).
!  All rights in the program are reserved by the DOE and the University.
!  Permission is granted to the public to copy and use this software
!  without charge, provided that this Notice and any statement of
!  authorship are reproduced on all copies. Neither the U.S. Government
!  nor the University makes any warranty, express or implied, or
!  assumes any liability or responsibility for the use of this software.
!***********************************************************************
CD1
CD1 PURPOSE
CD1
CD1 To determine the current date and time for the sun.
CD1
C**********************************************************************
CD2
CD2 REVISION HISTORY
CD2
CD2 Revision                    ECD
CD2 Date         Programmer     Number  Comments
CD2
CD2 05-20-92     G. Zyvoloski   00022   Initial implementation.
CD2                                     However, previous non-YMP
CD2                                     versions of FEHM exist, and
CD2                                     the current version differs
CD2                                     from these in minor ways.
CD2
CD2 $Log:   /pvcs.config/fehm90/src/dated.f_a  $
!D2
!D2    Rev 2.5   06 Jan 2004 10:42:48   pvcs
!D2 FEHM Version 2.21, STN 10086-2.21-00, Qualified October 2003
!D2
!D2    Rev 2.4   29 Jan 2003 09:02:10   pvcs
!D2 FEHM Version 2.20, STN 10086-2.20-00
!D2
!D2    Rev 2.3   14 Nov 2001 13:03:48   pvcs
!D2 FEHM Version 2.12, STN 10086-2.12-00
!D2 Update the GoldSim / FEHM interface
!D2
!D2    Rev 2.2   06 Jun 2001 13:35:56   pvcs
!D2 FEHM Version 2.11, STN 10086-2.11-00
!D2
!D2    Rev 2.2   06 Jun 2001 08:26:06   pvcs
!D2 Update for extended dispersion tensor model
!D2
!D2    Rev 2.1   30 Nov 2000 11:59:14   pvcs
!D2 FEHM Version 2.10, STN 10086-2.10-00
!D2
!D2    Rev 2.0   Fri May 07 14:40:28 1999   pvcs
!D2 FEHM Version 2.0, SC-194 (Fortran 90)
CD2
CD2    Rev 1.11   Tue May 07 15:59:40 1996   llt
CD2 added p for prototype to version label
CD2
CD2    Rev 1.10   Tue May 07 12:58:44 1996   llt
CD2 new version number
CD2
CD2    Rev 1.9   Mon Jan 29 14:25:50 1996   hend
CD2 Updated Requirements Traceability
CD2
CD2    Rev 1.8   09/13/95 09:33:34   llt
CD2 put "p" in version number - FEHMN 95-09-09p - for prototype
CD2
CD2    Rev 1.7   09/11/95 09:11:14   llt
CD2 updated version number
CD2
CD2    Rev 1.6   05/01/95 09:55:16   llt
CD2 protype version
CD2
CD2    Rev 1.5   05/01/95 09:21:28   llt
CD2 new version
CD2
CD2    Rev 1.4   04/25/95 09:26:38   llt
CD2 corrected log history information
CD2
CD2    Rev 1.3   02/01/95 19:02:20   llt
CD2 added prototype to version number
CD2
CD2    Rev 1.2   02/01/95 14:57:56   llt
CD2 new verno
CD2
CD2    Rev 1.1   01/28/95 15:22:08   llt
CD2 updated to new version
CD2
CD2    Rev 1.0   01/27/95 15:40:26   llt
CD2 remade archive to correct pvcs problem
CD2
C**********************************************************************
CD3
CD3 INTERFACES
CD3
CD3 Formal Calling Parameters
CD3
CD3 Identifier     Type    Use      Description
CD3
CD3 idatex         char     O       String containing the date
CD3                                    information
CD3 jtimex         char     O       String containing the time
CD3                                    information
CD3
CD3 Interface Tables
CD3
CD3 None
CD3
CD3 Files
CD3
CD3 None
CD3
C**********************************************************************
CD4
CD4 GLOBAL OBJECTS
CD4
CD4 Global Constants
CD4
CD4 None
CD4
CD4 Global Types
CD4
CD4 NONE
CD4
CD4 Global Variables
CD4
CD4 verno    FEHMN Version number
CD4
CD4 Global Subprograms
CD4
CD4 Name      Type       Description
CD4
CD4 idate     N/A        SUN system routine for determining the date
CD4 itime     N/A        SUN system routine for determining the time
CD4
C**********************************************************************
CD5
CD5 LOCAL IDENTIFIERS
CD5
CD5 Local Constants
CD5
CD5 None
CD5
CD5 Local Types
CD5
CD5 None
CD5
CD5 Local Variables
CD5
CD5 Identifier   Type        Description
CD5
CD5 jjdate       int         Returned date information in integer form
CD5 jjtime       int         Returned time information in integer form
CD5
CD5 Local Subprograms
CD5
CD5 None
CD5
C**********************************************************************
CD6
CD6 FUNCTIONAL DESCRIPTION
CD6
CD6 N/A
CD6
C**********************************************************************
CD7
CD7 ASSUMPTIONS AND LIMITATIONS
CD7
CD7 None
CD7
C**********************************************************************
CD8
CD8 SPECIAL COMMENTS
CD8
CD8 This routine is used for timing information only.
CD8
C**********************************************************************
CD9
CD9 REQUIREMENTS TRACEABILITY
CD9
CD9 Not Applicable.  See special Comments.
CD9
C**********************************************************************
CDA
CDA REFERENCES
CDA
CDA See FEHMN SRS, MMS, and SDD
CDA
C**********************************************************************
CPS
CPS PSEUDOCODE
CPS
CPS BEGIN dated
CPS
CPS idate - determine current date
CPS Perform internal write of date information into string
CPS itime - determine current time
CPS Perform internal write of time information into string
CPS set version number
CPS
CPS END dated
CPS
C**********************************************************************
      use comai
      implicit none

      character*8 current_date
      character*10 current_time
      character*11  idatex
      character*8   jtimex

c      integer*4  jjdate (3), jjtime (3)
c      idatex='        '
c      jtimex='        '

c      jjdate(1)=0
c      jjdate(2)=0
c      jjdate(3)=0

c      jjtime(1)=0
c      jjtime(2)=0
c      jjtime(3)=0

c      call  idate( jjdate )
c      jjdate(3) =  mod( jjdate(3),100 )
c      write(idatex,'(i2.2,1h/,i2.2,1h/,i2.2)')  jjdate(3) , jjdate(2)
c     *                                                   , jjdate(1)
c      call  itime( jjtime )
c      write(jtimex,'(i2.2,1h:,i2.2,1h:,i2.2)')  jjtime
      call date_and_time(current_date,current_time)
      jtimex(1:2) = current_time(1:2)
      jtimex(3:3) = ":"
      jtimex(4:5) = current_time(3:4)
      jtimex(6:6) = ":"
      jtimex(7:8) = current_time(5:6)
      idatex(1:2) = current_date(5:6)
      idatex(3:3) = "/"
      idatex(4:5) = current_date(7:8)
      idatex(6:6) = "/"
      idatex(7:8) = current_date(1:2)
      idatex(9:11) = current_date(3:4)
! Version number passed to GoldSim
      vernum = 3.4
! Code version identifier
      verno = "FEHM V3.4.0lnUbuntu16 19-09-12 QA:QA"

      end

max_stars_repo_name: satkarra/FEHM

max_issues_repo_name: satkarra/FEHM

max_forks_repo_name: satkarra/FEHM


text:
*
*===program crint======================================================*
*
C      OPTIONS/ EXTEND_SOURCE
C      SUBROUTINE CRINT
      SUBROUTINE DT_PRODUCEEVENT(ENERGY_SL, NPARTICLES)

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      REAL ENERGY_SL
      INTEGER INIT
      REAL ne,etest,prob,slump
      SAVE

* Call the init sub routine in the first event
      DATA INIT /0/

      PARAMETER (NMXHKK=200000)

      COMMON /DTIONT/ LINP,LOUT,LDAT

      COMMON /DTEVT1/ NHKK,NEVHKK,ISTHKK(NMXHKK),IDHKK(NMXHKK),
     &                JMOHKK(2,NMXHKK),JDAHKK(2,NMXHKK),
     &                PHKK(5,NMXHKK),VHKK(4,NMXHKK),WHKK(4,NMXHKK)

*     event flag
      COMMON /DTEVNO/ NEVENT, ICASCA

      IF(INIT.EQ.0) THEN
         OPEN (UNIT = 50, file = "my.input")
	 LINP = 50
         CALL DT_DTUINI(NEVTS,EPN,NPMASS,NPCHAR,NTMASS,NTCHAR,IDP,IEMU)
*        Init called, make sure it's not called again
         INIT = 1
      ENDIF
*-----------------------------------------------------------------------
*     generation of one event
      NEVENT = 1
      KKMAT = -1

*   If an energy-range has been defined with the ENERGY input-card the
*   laboratory energy ELAB can be set to any value within that range,..
C        ELAB = DT_RNDM(EPN)*(EPN-0.5D7)+0.5D7

*   ..otherwise it has to coincide with EPN.
C        ELAB = EPN

      ELAB = ENERGY_SL

*   sampling of one event

*     TEST

      CALL DT_KKINC(NPMASS,NPCHAR,NTMASS,NTCHAR,IDP,ELAB,KKMAT,IREJ)

      IF (IREJ.NE.0) RETURN

c     Return the number of particles produced

c     Fill the particle info
      CALL DT_GETPARTICLES(NPARTICLES)

      END

      SUBROUTINE DT_GETPARTICLES(NPARTICLES)

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      INTEGER pid,qch,q_sum,Ntpc,Nfinal,NACCEPT,IPART,RES
      DOUBLE PRECISION yrap,pt,mass,mt,etot
      DOUBLE PRECISION pt_cut_tpc
      PARAMETER(pt_cut_tpc=0.050)

      SAVE
*
* COMMON /DTEVT1/ :
*                   NHKK         number of entries in common block
*                   NEVHKK       number of the event
*                   ISTHKK(i)    status code for entry i
*                   IDHKK(i)     identifier for the entry
*                                (for particles: identifier according
*                                 to the PDG numbering scheme)
*                   JMOHKK(1,i)  pointer to the entry of the first mother
*                                of entry i
*                   JMOHKK(2,i)  pointer to the entry of the second mother
*                                of entry i
*                   JDAHKK(1,i)  pointer to the entry of the first daughter
*                                of entry i
*                   JDAHKK(2,i)  pointer to the entry of the second daughter
*                                of entry i
*                   PHKK(1..3,i) 3-momentum
*                   PHKK(4,i)    energy
*                   PHKK(5,i)    mass
*
* event history

      PARAMETER (NMXHKK=200000)

      COMMON /DTEVT1/ NHKK,NEVHKK,ISTHKK(NMXHKK),IDHKK(NMXHKK),
     &                JMOHKK(2,NMXHKK),JDAHKK(2,NMXHKK),
     &                PHKK(5,NMXHKK),VHKK(4,NMXHKK),WHKK(4,NMXHKK)

* extended event history
      COMMON /DTEVT2/ IDRES(NMXHKK),IDXRES(NMXHKK),NOBAM(NMXHKK),
     &                IDBAM(NMXHKK),IDCH(NMXHKK),NPOINT(10),
     &                IHIST(2,NMXHKK)

      DOUBLE PRECISION SLPX, SLPY, SLPZ, SLE, SLM
      INTEGER SLPID, SLCHARGE
      COMMON /DPMJETPARTICLE/ SLPX(NMXHKK), SLPY(NMXHKK), SLPZ(NMXHKK),
     &       SLE(NMXHKK), SLM(NMXHKK), SLPID(NMXHKK), SLCHARGE(NMXHKK)

C     >> Set Counter to Zero

      Nfinal=0

      DO 42 I=1, NHKK
c      I = IPART

CC       >> Remove all non-final-state particles
        IF(.not.(ISTHKK(I).eq.1.or.ISTHKK(I).eq.-1.or.
     $ISTHKK(I).eq.1001)) GOTO 42

C	>> Find Particle Charge, qch
        IF((ABS(ISTHKK(I)).eq.1).and.(IDHKK(I).ne.80000))THEN
C         >> final state ptcles except nuclei

          qch=IPHO_CHR3(IDHKK(I),1)/3
        ELSEIF(IDHKK(I).eq.80000)THEN
C         >> final state nuclei
          qch=IDXRES(I)
        ELSE
C         >> not a final state particle, qch not interesting
          qch=-999
        ENDIF

	Nfinal = Nfinal + 1
	SLPX(Nfinal) = PHKK(1,I)
        SLPY(Nfinal) = PHKK(2,I)
        SLPZ(Nfinal) = PHKK(3,I)
        SLE(Nfinal) = PHKK(4,I)
        SLM(Nfinal) = PHKK(5,I)
        SLPID(Nfinal) = IDHKK(I)
        SLCHARGE(Nfinal) = qch

 42     CONTINUE
        NPARTICLES = Nfinal

      END

      SUBROUTINE DT_USRHIS(MODE)
c Dummy to make the linker happy
      END

max_stars_repo_name: AllaMaevskaya/AliRoot

max_issues_repo_name: AllaMaevskaya/AliRoot

max_forks_repo_name: AllaMaevskaya/AliRoot


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

module unit

    ! This module provides two subroutines: claim_unit( u ) and release_unit( u ) to
    ! get an unused unit number and set it free again.

    ! The unit numbers obtained by claim_unit( u ) are (currently) in the range [10000,90000]
    ! so the calling program should never use any unit number in this range explicitly.

    ! It also provides the constants stderr, stdin and stdout, wich should be set to the compiler's
    ! error, in and out preconnected units.

    ! LIMITATIONS
    ! (should be OK for all normal purposes,
    !  hardly anybody needs 10000 open files at once, and chance is, that
    !  the OS won't allow it anyway...)
    ! This will fail, after 10001 units have been claimed (but not released).
    ! Performance may become bad when lots of units have been claimed (but not released).

    use util

    implicit none

    ! access specifiers

    private
    public :: claim_unit, release_unit

    ! type defs

    type, private :: unit_t
        private
        integer :: number
        type(unit_t), pointer :: next
    end type

    ! globals constants

    integer, dimension(2), private, parameter :: unit_range = (/10000,20000/)

    ! global variables

    integer, private :: current = 10000
    type(unit_t), pointer :: used_units => null()

  contains

    subroutine claim_unit(u)

        integer, intent(out) :: u
        integer :: count

        ! get an unused unit number

        count = 0
        do while ( isin( used_units, current ) )
            current = current + 1
            if (current > unit_range(2)) current = unit_range(1)
            count = count+1
            if ( count > (unit_range(2)-unit_range(1)+1) ) then
                stop 'unit.f90: too many claimed units'
            end if
        end do
        call push( used_units, current )
        u = current

    end subroutine claim_unit

    subroutine release_unit(u)

        integer, intent(in) :: u

        ! free a unit number

        call remove_by_val( used_units, u )

    end subroutine release_unit

! stuff below here should be in separate module
! cause it is only a very basic (and slow) implementation of a list container
! should be provided by some standard library

    subroutine push( list, i )

        type(unit_t), pointer :: list
        integer, intent(in) :: i
        type(unit_t), pointer :: current, new

        ! add an item to the end of the list

        allocate( new )
        new%number = i
        new%next => null()

        if (.not. associated( list ) ) then
            list => new
        else
            current => list
            do while (associated(current%next))
                current => current%next
            end do
            current%next => new
        end if

    end subroutine push

    subroutine remove_by_val( list, i )

        type(unit_t), pointer :: list
        integer, intent(in) :: i

        ! remove all items with value i from the list

        type(unit_t), pointer :: current, del

        if ( .not. associated(list) ) return

        do while (list%number == i)
            current => list
            list => current%next
            deallocate(current)
            if ( .not. associated(list) ) exit
        end do

        if ( .not. associated(list) ) return

        current => list
        do while (associated(current%next))
            if (current%next%number == i) then
                del => current%next
                current%next => current%next%next
                deallocate( del )
            else
                current => current%next
            end if
        end do

    end subroutine remove_by_val

    function isin( list, i ) result(yn)

        type(unit_t), pointer :: list
        integer, intent(in) :: i
        logical :: yn

        ! check if i is in list

        type(unit_t), pointer :: current

        current => list
        yn = .false.

        do while ( associated(current) )
            if ( current % number == i ) then
                yn = .true.
                return
            end if
            current => current%next
        end do

    end function isin

end module unit

max_stars_repo_name: emolch/kiwi

max_issues_repo_name: alirezaniki/kiwi

max_forks_repo_name: emolch/kiwi


text:
Cedar waxwings are the species B. cedrorum and live all over the United States and sometimes in South America. They can be seen yearround in Davis. They eat berries, fruit, and insects. According to Wikipedia, they make very highpitched whistles and buzzy trills often represented as see or sree.

A tree between the Wells Fargo parking lot and the Davis Enterprise lot is full of these cute lil birdies.

These are pretty cool looking birds and I like the noises they make, but they crap all over my car and sometimes on me. I cant park in the driveway anymore. Its like a rain of a bird poo when theyre around. Users/MikeyCrews


text:
import base

let main to-List(apply-pow2n-times(n20, map(inc), list-with-100-zeros))

max_stars_repo_name: rinor/Formality

max_issues_repo_name: rinor/Formality

max_forks_repo_name: rinor/Formality


text:
! -*- mode: F90; mode: font-lock -*-
! ------------------------------------------------------------------------------
! $Id$
! ------------------------------------------------------------------------------
! Module potential_module
! ------------------------------------------------------------------------------
! Code area 3: operators
! ------------------------------------------------------------------------------

!!****h* Conquest/potential_module *
!!  NAME
!!   potential_module
!!  PURPOSE
!!   Hold routines and data associated with local potential (this
!!   means the total potential on the grid at the moment, but may grow
!!   to include neutral atom potential in the future).  Will hold
!!   local potential expansion ideas of Ozaki in the near future.
!!  AUTHOR
!!   D. R. Bowler
!!  CREATION DATE
!!   2006/08/02 17:28 dave
!!  MODIFICATION HISTORY
!!   2011/04/01 L.Tong
!!    Added potential_dn for spin polarisaed calculations
!!   2012/03/14 L.Tong
!!    Changed spin implementation. Now potential has a second spin
!!    index to denote different channel
!!  SOURCE
!!
module potential_module

  use datatypes

  implicit none
  save

  real(double), allocatable, dimension(:,:) :: potential

!!***

end module potential_module

max_stars_repo_name: negf/CONQUEST-release

max_issues_repo_name: negf/CONQUEST-release

max_forks_repo_name: negf/CONQUEST-release


text:
These people (they can be women too!) kick your butt for misbehaving/or being too drunk (stumbling, passing out, throwing up) at bars. And they take your money and make sure you dont have fake IDs.

Once your fake ID is confiscated, just leave, if you pester the bouncers to try to get it back, the Davis Police Department may be called..

Attacking a bouncer is one of the stupidest things a person can do. Not only will you anger the bouncer cause he has to stop discussing his AWESOME bouncer judo with all other bouncers, he will attempt to fight back, and then all the other bouncers will pounce upon you and fight (read reasonable defend themselves cause a working somewhere make 242 p.c. not apply) you as well. You will be arrested for battery when they lie about everything that goes down. If you get kicked out of a bar, let it slide and just go to another bars bar ...

Seriously thou, dont act a fool, maintain composure, and bouncers will usually have nothing but a scowly version of love for you

  Sophias  Sophias Thai Kitchen needs bouncers? Users/JimStewart Its a Thai Bar as well Daubert
   Justin
   Jim
   Ben
   Travis
   John

  The G Street Wunderbar
   Thomas
   James
   JoshLawlson
   Stephen
   Jay
   Mike
   Zak
   Bouncer emeritus: Owen (he used to work there when it was the G St. Pub, and hes my brother) Users/jsbmeb

  Froggys
  Vincent
  Marlon
  Steven

  Bistro 33 Bistro 33
   Derrick
  Aaron
  Colin
  Wade
  Austin

  de Veres Irish Pub de Veres
   Ed

 KetMoree
  Ben (Big Bad Brunz)

 Tres Hermanas
  D

 Little Prague

 Uncle Vitos Craig works the door, and is a champion

 Luigis

 Woodstocks Pizza  Users/WesP Wes  Formerly of Little Prague Bar, Wes is a well known face in downtown Davis and has now moved over to Woodstocks as of 5/14/2012.


text:
      subroutine disp_ono(yplot)

      use size_mod

      implicit none

      integer nmax, npoints, nline, i, nsum, n, kmax, kn, jmid
      integer pgopen, pgbeg, ier, nlevmax, iflag, numb
      integer nxmx, nymx, nnodex, nnodey, j, jplot

      parameter (nxmx = nmodesmax)
      parameter (nymx = mmodesmax)
      parameter (nmax = nmodesmax)

      parameter (nlevmax = 101)

      real xgrid(nmax), ygrid(nmax), f(nmax), dx, xmin, xmax, capf(nmax)
      real capr(nmax), ff(101), rhoplasm, dy, yplot, ybottom, yprime
      complex xkperp2_slow(nxmx, nymx), xkperp2_fast(nxmx, nymx)
      complex xkperp_slow(nxmx, nymx), xkperp_fast(nxmx, nymx)
      complex P(nxmx, nymx)
      real xkprl(nxmx, nymx)

      complex xkperp2_slow_plot(nxmx), xkperp2_fast_plot(nxmx)
      complex xkperp_slow_plot(nxmx), xkperp_fast_plot(nxmx)

      real rho(nxmx, nymx)

      real f3(nmax)
      real y(nmax), k, capk, c, xdeg, yh, L, kr, nr
      real factrl, xn, an, diff1, diff2

      real ans_simpson, ans_trapezoidal, ans_analytic, ymax, ymin
      real dummy, time, t1, t2, second1, pi, x, a, e, dydx, sum

      character*32 title
      character*32 titlx
      character*32 tityl
      character*32 tityr

      character*32 titx
      character*32 tity
      character*32 titz

      integer nblack, nred, nyellow, ngreen, nblue, ncyan, nmagenta,
     .   nwhite, norange

      common/boundcom/rhoplasm

      open(unit=130, file='Ono_disp', status='unknown',
     .      form='formatted')

      nwhite = 0
      nblack = 1
      nred = 2
      ngreen = 3
      nblue = 4
      ncyan = 5
      nmagenta = 6
      nyellow = 7
      norange = 8

      numb = 100

      pi = 3.14159
      e =  2.71828
      a = 3.0
      k = .0015
      capk = 6000.
      f3 = 0.

      read(130, 309) nnodex, nnodey, jmid
      read(130, 310) rhoplasm
      read(130, 310) (capr(i), i = 1, nnodex)
      read(130, 310) (y(j), j = 1, nnodey)
      read(130, 310) ((xkperp2_slow(i, j), i = 1, nnodex),
     .      j = 1, nnodey)
      read(130, 310) ((xkperp2_fast(i, j), i = 1, nnodex),
     .      j = 1, nnodey)
      read(130, 310) ((rho(i, j), i = 1, nnodex), j = 1, nnodey)
      read(130, 310) ((xkprl(i,j), i = 1, nnodex), j = 1, nnodey)
      read(130, 310) ((P(i, j), i = 1, nnodex), j = 1, nnodey)

      xgrid = capr
      ygrid = y
      npoints = nnodex

      xmin = xgrid(1)
      xmax = xgrid(nnodex)

      jplot = jmid

      dy = y(2) - y(1)
      ybottom = y(1)
c      yplot = 0.05
      yprime = yplot - ybottom
      jplot = yprime / dy + 1
      write(6, *) "jplot = ", jplot
      write(6, *) "yplot = ", yplot

!     -------------------------------
!     Open the pgplot graphics device
!     -------------------------------
      IER = PGBEG(0, 'disper.ps/vcps', 1, 1)
      IF (IER.NE.1) STOP

      call PGSCH (1.5)

!     -------------------------------------------
!     Plot contours of k**2 - slow root
!     -------------------------------------------
      titx = 'R (m)'
      tity = 'Z (m)'
      title = 'Re(k**2) - slow root'
      call ezconc(capr, y, real(xkperp2_slow), ff, nnodex, nnodey, numb,
     .   nxmx, nymx, nlevmax, title, titx, tity, iflag)
      if (iflag .eq. 0) call boundary(capr, y, rho, ff, nnodex,
     .   nnodey, numb,
     .   nxmx, nymx, nlevmax, title, titx, tity)

      title = 'k_parallel'
      call ezconc(capr, y, xkprl, ff, nnodex, nnodey, numb,
     .   nxmx, nymx, nlevmax, title, titx, tity, iflag)
      if (iflag .eq. 0) call boundary(capr, y, rho, ff, nnodex,
     .   nnodey, numb,
     .   nxmx, nymx, nlevmax, title, titx, tity)

      title = 'Re(P)'
      call ezconc(capr, y, real(P), ff, nnodex, nnodey, numb,
     .   nxmx, nymx, nlevmax, title, titx, tity, iflag)
      if (iflag .eq. 0) call boundary(capr, y, rho, ff, nnodex,
     .   nnodey, numb,
     .   nxmx, nymx, nlevmax, title, titx, tity)

      title = 'Im(P)'
      call ezconc(capr, y, aimag(P), ff, nnodex, nnodey, numb,
     .   nxmx, nymx, nlevmax, title, titx, tity, iflag)
      if (iflag .eq. 0) call boundary(capr, y, rho, ff, nnodex,
     .   nnodey, numb,
     .   nxmx, nymx, nlevmax, title, titx, tity)

      title = 'Im(k**2) - slow root'
      call ezconc(capr, y, aimag(xkperp2_slow), ff, nnodex, nnodey,numb,
     .   nxmx, nymx, nlevmax, title, titx, tity, iflag)
      if (iflag .eq. 0) call boundary(capr, y, rho, ff, nnodex,
     .   nnodey, numb,
     .   nxmx, nymx, nlevmax, title, titx, tity)

!     ------------------------
!     Plot cut in y at yplot
!     ------------------------
      do i = 1, nnodex
         xkperp2_slow_plot(i) = xkperp2_slow(i, jplot)
         xkperp2_fast_plot(i) = xkperp2_fast(i, jplot)
      end do

      xkperp_slow_plot = csqrt(xkperp2_slow_plot)
      xkperp_fast_plot = csqrt(xkperp2_fast_plot)

      titlx  = 'R (m)'
      tityl = 'kperp2 (m-2)'
      tityr = ' '
      title = 'fast and slow roots'

      call ezplot2_s(title, titlx, tityl, tityr, xgrid,
     .   real(xkperp2_slow_plot), real(xkperp2_fast_plot), f3,
     .   npoints, nmax, ymax, ymin, xmin, xmax)

      tityl = 'kperp (m-1)'

      call ezplot2_b(title, titlx, tityl, tityr, xgrid,
     .   real(xkperp_slow_plot), real(xkperp_fast_plot), f3,
     .   npoints, nmax, ymax, ymin, xmin, xmax)

c      title = 'slow root'
c     call ezplot1_s(title, titlx, tityl, tityr, xgrid,
c     .   real(xkperp2_slow_plot), f3,
c     .   npoints, nmax, ymax, ymin, xmin, xmax, nred)

      tityl = 'kperp2 (m-2)'
      title = 'fast root'
      call ezplot1_s(title, titlx, tityl, tityr, xgrid,
     .   real(xkperp2_fast_plot), f3,
     .   npoints, nmax, ymax, ymin, xmin, xmax, nblue)

!     -------------------------------------------
!     Plot contours of k**2 - fast root
!     -------------------------------------------

      title = 'Re(k**2) - fast root'
      call ezconc(capr, y, real(xkperp2_fast), ff, nnodex, nnodey, numb,
     .   nxmx, nymx, nlevmax, title, titx, tity, iflag)
      if (iflag .eq. 0) call boundary(capr, y, rho, ff, nnodex,
     .   nnodey, numb,
     .   nxmx, nymx, nlevmax, title, titx, tity)

      title = 'Im(k**2) - fast root'
      call ezconc(capr, y, aimag(xkperp2_fast), ff, nnodex, nnodey,numb,
     .   nxmx, nymx, nlevmax, title, titx, tity, iflag)
      if (iflag .eq. 0) call boundary(capr, y, rho, ff, nnodex,
     .   nnodey, numb,
     .   nxmx, nymx, nlevmax, title, titx, tity)

      xkperp_slow = csqrt(xkperp2_slow)
      xkperp_fast = csqrt(xkperp2_fast)

      do i = 1, nnodex
         xkperp_slow_plot(i) = xkperp_slow(i, jplot)
         xkperp_fast_plot(i) = xkperp_fast(i, jplot)
      end do

!     -------------------------------------------------
!     write the function xkperp_fast_plot(i) to screen
!     -------------------------------------------------
      do i = 1, npoints
	 write(6, '(i10, 1p5e15.5)') i, xgrid(i), xkperp_fast_plot(i)
      end do

      title = 'Re(k) - fast and slow roots'
      titlx  = 'R (m)'
      tityl = 'Re(k) (m-1)'
      tityr = ' '

      call ezplot2_sym(title, titlx, tityl, tityr, xgrid,
     .   real(xkperp_slow_plot), real(xkperp_fast_plot), f3,
     .   npoints, nmax, 0.0, 0.0, xmin, xmax)

      title = 'Re(k) - fast and slow (blowup)'
      titlx  = 'R (m)'
      tityl = 'Re(k) (m-1)'
      tityr = ' '

      call ezplot2_sym(title, titlx, tityl, tityr, xgrid,
     .   real(xkperp_slow_plot), real(xkperp_fast_plot), f3,
     .   npoints, nmax, 2000., -2000., xmin, xmax)

      title = 'Re(k) - slow root'
      call ezconc(capr, y, real(xkperp_slow), ff, nnodex, nnodey, numb,
     .   nxmx, nymx, nlevmax, title, titx, tity, iflag)
      if (iflag .eq. 0) call boundary(capr, y, rho, ff, nnodex,
     .   nnodey, numb,
     .   nxmx, nymx, nlevmax, title, titx, tity)

!     -------------------------
!     Close the graphics device
!     -------------------------
      call pgclos

      close (130)

  310 format(1p6e12.4)
 3310 format(1p6e18.10)
 8310 format(1p6e14.6)
  309 format(10i10)

      return
      end
!
!*******************************************************************************
!

      subroutine ezplot1_s(title, titx, titll, titlr, x1, y1, y3,
     .   nr, nrmax,
     .   ymax, ymin, xmin, xmax, ncolor)

      implicit none

      real xzmax,xzmin,xnmin,xnmax,rhomin,rhomax
      real x1(nrmax), y1(nrmax), y3(nrmax)
      real y1max,y2max,y3max,y1min,y2min,y3min
      real ymin, ymax, xmin, xmax

      character*32 title
      character*32 titx
      character*32 titll
      character*32 titlr

      integer nr, nrmax

      integer nplot1,ncollab, ncolion,ncolbox, ncyan, ncolelec,
     .   ncolln2, ncollin, ncolbrd

      integer nblack, nred, nyellow, ngreen, naqua, npink, nwheat,
     .   ngrey, nbrown, nblue, nblueviolet, ncyan1, ncolor
      integer nturquoise, nmagenta, nsalmon, nwhite, norange, ncolln3

      nwhite = 0
      nblack = 1
      nred = 2
      ngreen = 3
      nblue = 4
      ncyan = 5
      nmagenta = 6
      nyellow = 7
      norange = 8

      ncolln3 = ngreen

      call a1mnmx(y1, nrmax, nr, y1min, y1max)
      if(y1max .eq. 0.0 .and. y1min .eq. 0.0)return

      ymax = y1max
      ymin = y1min

      rhomax = xmax
      rhomin = xmin

! Advance plotter to a new page, define coordinate range of graph and draw axes

      CALL PGPAGE
      CALL PGSVP (0.15,0.85,0.15,0.85)
      CALL PGSWIN (rhomin, rhomax, ymin, ymax)
      CALL PGBOX   ('BCNST', 0.0, 0,  'BCNST', 0.0, 0)

      CALL PGSCI(nblack)

! Label the axes (note use of \u and \d for raising exponent).

      call pglab(titx, titll, title)

! Plot the line graph.

      CALL PGSCI(ncolor)
      call pgline(nr, x1, y1)

  300 format (1p9e11.3)

      CALL PGSCI(nblack)
      call pgline(nr, x1, y3)

      return
      end

!
!********************************************************************
!
      subroutine ezplot2_s(title, titlb, titll, titlr, x1, y1, y2, y3,
     .   nr, nrmax, ymax, ymin, xmin, xmax)

      implicit none

      real xzmax,xzmin,xnmin,xnmax,rhomin,rhomax
      real x1(nrmax),y1(nrmax), y2(nrmax), y3(nrmax)
      real y1max,y2max,y3max,y1min,y2min,y3min
      real ymin,ymax, xmin, xmax

      character*32 title
      character*32 titll
      character*32 titlr
      character*32 titlb

      integer nr, nrmax, n

      integer nplot1,ncollab, ncolion,ncolbox, ncyan, ncolelec,
     .   ncolln2, ncollin, ncolbrd
      integer nblack,nred,nyellow, ngreen,naqua,npink, nwheat,
     .   ngrey,nbrown,nblue,nblueviolet,ncyan1
      integer nturquoise,nmagenta,nsalmon,nwhite,ncolln3, norange

      nwhite = 0
      nblack = 1
      nred = 2
      ngreen = 3
      nblue = 4
      ncyan = 5
      nmagenta = 6
      nyellow = 7
      norange = 8

      call a1mnmx(y1, nrmax, nr, y1min, y1max)

      if(y1max .eq. 0.0 .and. y1min .eq. 0.0)return

      ymax = y1max
      ymin = y1min

      rhomax = xmax
      rhomin = xmin

      CALL PGPAGE
      CALL PGSVP (0.15,0.85,0.15,0.85)
      CALL PGSWIN (rhomin, rhomax, ymin, ymax)
      CALL PGBOX  ('BCNST', 0.0, 0, 'BCNST', 0.0, 0)

      CALL PGSCI(nblack)

! Label the axes (note use of \u and \d for raising exponent).

      call pglab(titlb, titll, title)

      CALL PGSCI(nred)
      call pgline(nr, x1, y1)

      CALL PGSCI(nblue)
      call pgline(nr, x1, y2)

      CALL PGSCI(nblack)
      call pgline(nr, x1, y3)

  300 format (1p9e11.3)
      return
      end

!
!********************************************************************
!
      subroutine ezplot2_b(title, titlb, titll, titlr, x1, y1, y2, y3,
     .   nr, nrmax, ymax, ymin, xmin, xmax)

      implicit none

      real xzmax,xzmin,xnmin,xnmax,rhomin,rhomax
      real x1(nrmax),y1(nrmax), y2(nrmax), y3(nrmax)
      real y1max,y2max,y3max,y1min,y2min,y3min
      real ymin,ymax, xmin, xmax

      character*32 title
      character*32 titll
      character*32 titlr
      character*32 titlb

      integer nr, nrmax, n

      integer nplot1,ncollab, ncolion,ncolbox, ncyan, ncolelec,
     .   ncolln2, ncollin, ncolbrd
      integer nblack,nred,nyellow, ngreen,naqua,npink, nwheat,
     .   ngrey,nbrown,nblue,nblueviolet,ncyan1
      integer nturquoise,nmagenta,nsalmon,nwhite,ncolln3, norange

      nwhite = 0
      nblack = 1
      nred = 2
      ngreen = 3
      nblue = 4
      ncyan = 5
      nmagenta = 6
      nyellow = 7
      norange = 8

      call a1mnmx(y2, nrmax, nr, y2min, y2max)

      if(y1max .eq. 0.0 .and. y1min .eq. 0.0)return

      ymax = y2max * 2.0
      ymin = y2min

      rhomax = xmax
      rhomin = xmin

      CALL PGPAGE
      CALL PGSVP (0.15,0.85,0.15,0.85)
      CALL PGSWIN (rhomin, rhomax, ymin, ymax)
      CALL PGBOX  ('BCNST', 0.0, 0, 'BCNST', 0.0, 0)

      CALL PGSCI(nblack)

! Label the axes (note use of \u and \d for raising exponent).

      call pglab(titlb, titll, title)

      CALL PGSCI(nred)
      call pgline(nr, x1, y1)

      CALL PGSCI(nblue)
      call pgline(nr, x1, y2)

      CALL PGSCI(nblack)
      call pgline(nr, x1, y3)

  300 format (1p9e11.3)
      return
      end

!
!*****************************************************************
!
      subroutine ezplot2_sym(title, titlb, titll, titlr, x1, y1, y2, y3,
     .   nr, nrmax, ymax_giv, ymin_giv, xmin, xmax)

      implicit none

      real xzmax,xzmin,xnmin,xnmax,rhomin,rhomax
      real x1(nrmax),y1(nrmax), y2(nrmax), y3(nrmax)
      real y1max,y2max,y3max,y1min,y2min,y3min
      real ymin,ymax, xmin, xmax, ymax_giv, ymin_giv

      character*32 title
      character*32 titll
      character*32 titlr
      character*32 titlb

      integer nr, nrmax, n

      integer nplot1,ncollab, ncolion,ncolbox, ncyan, ncolelec,
     .   ncolln2, ncollin, ncolbrd
      integer nblack,nred,nyellow, ngreen,naqua,npink, nwheat,
     .   ngrey,nbrown,nblue,nblueviolet,ncyan1
      integer nturquoise,nmagenta,nsalmon,nwhite,ncolln3, norange

      nwhite = 0
      nblack = 1
      nred = 2
      ngreen = 3
      nblue = 4
      ncyan = 5
      nmagenta = 6
      nyellow = 7
      norange = 8

      call a1mnmx(y1, nrmax, nr, y1min, y1max)

!      write(6, 300) y2min, y2max
      if(y1max .eq. 0.0 .and. y1min .eq. 0.0)return

      ymax = y1max * 1.1
      ymin = y1min

      if(ymax_giv .ne. 0.0)ymax = ymax_giv
      ymin = -ymax

      rhomax = xmax
      rhomin = xmin

      CALL PGPAGE
      CALL PGSVP (0.15,0.85,0.15,0.85)
      CALL PGSWIN (rhomin, rhomax, ymin, ymax)
      CALL PGBOX  ('BCNST', 0.0, 0, 'BCNST', 0.0, 0)

      CALL PGSCI(nblack)

! Label the axes (note use of \u and \d for raising exponent).

      call pglab(titlb, titll, title)

      CALL PGSCI(nred)
      call pgline(nr, x1, y1)
      call pgline(nr, x1, -y1)

      CALL PGSCI(nblue)
      call pgline(nr, x1, y2)
      call pgline(nr, x1, -y2)

      CALL PGSCI(nblack)
      call pgline(nr, x1, y3)

  300 format (1p9e11.3)
      return
      end

!
!********************************************************************
!

max_stars_repo_name: ORNL-Fusion/aorsa

max_issues_repo_name: ORNL-Fusion/aorsa

max_forks_repo_name: ORNL-Fusion/aorsa


text:
La Raza PreLaw Student Association (LRPLSA) serves as a conduit for a collective School of Law prelaw student voice. They work on issues of discrimination such as race, religion, and sex affecting Latin@ and other disadvantaged communities. La Raza PreLaw is a nonprofit Student Organizations student organization with a progressive agenda that seeks to protect and advance the civil rights of individuals who are underrepresented while undergraduates undergrads at the UC Davis University of California, Davis and in respective communities. La Raza PreLaw also seeks to bring about positive change by addressing legal issues dealing with, but not limited to immigration, schools education, Medical Services healthcare, and employment issues. LRPSA works to promote and sustain the Academic Help academic success of any prospective law student at UC Davis.

Meetings

Our regular meeting time is every Monday during the school year from 67pm in Olson 105 (20132014)

Board Members (20142015)

  CoChairs: Astrid Robles, Sara Ramirez
  Secretary: Monserrat Garcia
  Publicity: Joaquin Pacheco
  Academic Chair: Alejandro
  Treasurer: Adriana Reyna
  Social and Banquet Chair: Tiffany Avila
  Fundraising Chair: Sabina Garcia
  La Raza Cultural Days: vacant
  Political Chair: Antonio Alvarado

Contact Information (20142015)

Please contact one of the cochairs for more information to join!

ajrobles@ucdavis.edu or syramirez@ucdavis.edu

La Raza PreLaw email: razaprelaw@ucdavis.edu

Facebook page: http://www.facebook.com/groups/333486620019955/

Follow us on instagram @larazaprelaw

Website page:http://lrplsaucdavis.webs.com/

History

Board Members (20132014)

  CoChairs: Anthony Rascon, Veronica Serrato
  Secretary: Liliana Robles
  Publicity: Astrid Robles
  Academic Chair: Dania Herrera
  Treasurer: Sabina Garcia
  Social and Banquet Chair: Melissa Castillo
  Fundraising Chair: Sara Ramirez
  La Raza Cultural Days: Sara Ramirez
  Political Chair: Antonio Alvarado

Board Members (20122013)

  CoChairs: Alejandra Duran, Edith Santibáñez
  Secretary: Dianna Rivas
  Publicity: Maricela Hernandez
  Academic Chair: Dania Herrera
  Treasurer: Anthony Rascon
  Social and Banquet Chair: Veronica Serrato
  Fundraising Chair/ La gran tardeada chair: Sara Hernandez
  La Raza Cultural Days: Marissa Cuevas
  Political Chair:vacant

Board Members (20112012)

  CoChairs: Arasele Torrez, Kevin Finkelberg
  Secretary: Alejandra Duran
  Publicity: Edith Santibañez
  Academic Chair: Kristy Ayala
  Treasurer: Dianna Rivas
  Social and Banquet Chair: Faustino Valdivia
  Fundraising Chair: Sara Hernandez
  Political Chair: Angel Castro

Board Members (20092010)

  CoChairs: Andrea Delgadillo, Sara Kohgadai
  Club Financial Officer: Baryo Dee
  Secretary: Maggie Tapia
  Publicity: Efren Barragan
  Academic Chair: Carlos Castillo
  Social and Banquet Chair: Arasele Torrez
  DREAM Scholarship CoChairs: Vacant


text:
! mpiexec -np 4 bin/Debug/WGFM_test2.exe
program matrices_test
    use meshread_mod
    use EFIE_RWG_matrices
    use MFIEDuffy_functions
    use linalg_mod
    use spherePEC_mod
    use WMFIE_matrices
    use WGFM_matrices
    use mpi
    implicit none

    real(8), parameter :: pi = 4d0 * atan(1d0)
    complex(8), parameter :: IU = (0d0, 1d0)

    ! Variables
    type(Mesh) :: msh
    integer, parameter :: nNS = 6           ! number of quadrature points for EFIE

    complex(8), allocatable :: K1Mat(:, :)   ! MFIE K_1 operator
    complex(8), allocatable :: K2Mat(:, :)   ! MFIE K_2 operator
    complex(8), allocatable :: T1Mat(:, :)   ! EFIE T_1 operator
    complex(8), allocatable :: T2Mat(:, :)   ! EFIE T_2 operator

    complex(8), allocatable :: A1(:, :)
    complex(8), allocatable :: A2(:, :)
    complex(8), allocatable :: A3(:, :)
    complex(8), allocatable :: Z1(:, :)
    complex(8), allocatable :: Z2(:, :)
    complex(8), allocatable :: Z3(:, :)

    ! MPI variables
    integer :: ierr, N_procs, id

    ! Parameters
    real(8), parameter :: ep0 = 1d0 / (35950207149.4727056d0 * PI)  ! vacuum permittivity [F/m]
    real(8), parameter :: mu0 = PI * 4d-7                           ! vacuum permeability [H/m]

    real(8), parameter :: k = 0.1d0 * PI          ! free-space wavenumber [1/m]
    real(8), parameter :: ep1 = ep0             ! upper half-space permittivity [F/m]
    real(8), parameter :: mu1 = mu0             ! upper half-space permeability [H/m]
    real(8), parameter :: ep2 = ep0 * 2.5d0     ! lower half-space permittivity [F/m]
    real(8), parameter :: mu2 = mu0 * 5.3       ! lower half-space permeability [H/m]

    real(8), parameter :: w = k / sqrt(ep0 * mu0)   ! angular frequency [rad/s]
    real(8), parameter :: k1 = w * sqrt(ep1 * mu1)  ! upper half-space wavenumber [1/m]
    real(8), parameter :: k2 = w * sqrt(ep2 * mu2)  ! lower half-space wavenumber [1/m]

    character(len=100) :: file_msh = 'meshes/half_sphere_plane/half_sphere_plane_a5_h04_t20.msh'
    !character(len=100) :: file_msh = 'meshes/convergence_sphere/simple_mesh.msh'

    ! Set window (just to test the program)
    call set_circ_window(Lr=4d0, c_in=0.7d0)

    ! MPI init
    call mpi_init(ierr)
    call mpi_comm_size(MPI_COMM_WORLD, N_procs, ierr)
    call mpi_comm_rank(MPI_COMM_WORLD, id, ierr)
    print *, "Total processes:", N_procs, "Process id:", id

    ! Load mesh
    call load_gmsh(file_msh, msh)

    ! Compute MFIE matrices
    if (id == 0) print *, "MFIE"
    call genMFIESauterSchwabMat(K1Mat, msh, k1)
    call genMFIESauterSchwabMat(K2Mat, msh, k2)

    ! Compute EFIE matrices
    if (id == 0) print *, "EFIE"
    call genEFIESauterSchwabMat_RWG(T1Mat, msh, k1, nNS, nNS)
    call genEFIESauterSchwabMat_RWG(T2Mat, msh, k2, nNS, nNS)

    ! Compute WGFM matrices
    if (id == 0) print *, "WGFM"
    call genWGFMMat(msh, k1, k2, ep1, ep2, mu1, mu2, Z1, Z2, Z3)

    if (id == 0) then
        ! Compute WGFM matrices
        allocate(A1(msh%nbEdg, msh%nbEdg))
        allocate(A2(msh%nbEdg, msh%nbEdg))
        allocate(A3(msh%nbEdg, msh%nbEdg))
        A1 = mu2 * K2Mat - mu1 * K1Mat
        A2 = k2**2 * T2Mat - k1**2 * T1Mat
        A3 = ep2 * K2Mat - ep1 * K1Mat

        ! Compute matrices norm
        print *, "size", msh%nbEdg
        print *, "k", k
        print *, "w", w
        print *, "A1", sqrt(sum(abs(A1) ** 2))
        print *, "A2", sqrt(sum(abs(A2) ** 2))
        print *, "A3", sqrt(sum(abs(A3) ** 2))
        print *, "Z1", sqrt(sum(abs(Z1) ** 2))
        print *, "Z2", sqrt(sum(abs(Z2) ** 2))
        print *, "Z3", sqrt(sum(abs(Z3) ** 2))
        print *, "A1-Z1", sqrt(sum(abs(A1-Z1) ** 2))
        print *, "A2-Z2", sqrt(sum(abs(A2-Z2) ** 2))
        print *, "A3-Z3", sqrt(sum(abs(A3-Z3) ** 2))

    end if

    call MPI_Finalize(ierr)
end program

max_stars_repo_name: Riarrieta/WindowedGreenFunctionMethod

max_issues_repo_name: Riarrieta/WindowedGreenFunctionMethod

max_forks_repo_name: Riarrieta/WindowedGreenFunctionMethod


text:
      SUBROUTINE getlim
      USE vmec_main
      USE realspace
      USE vsvd
      USE mgrid_mod, ONLY: nlim, limitr, rlim, zlim, reslim, seplim
      IMPLICIT NONE
C-----------------------------------------------
C   L o c a l   P a r a m e t e r s
C-----------------------------------------------
      REAL(rprec), PARAMETER :: resup = 3.0_dp
      REAL(rprec), PARAMETER :: resdn = 0.5_dp
      REAL(rprec), PARAMETER :: eps = 0.005_dp
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER :: ntheta_2pi, nphi_plane, i,
     1   nthtot, iexpand, ishrink, ionlim, n,
     2   limpts, nonlim, nexpand, nshrink, ilim0, nlim0
      REAL(rprec), DIMENSION(2*ntheta1) ::
     1   rbdy, zbdy, rubdy, zubdy
      REAL(rprec) :: fshrink, distmax, fexpand
C-----------------------------------------------

c
c     DETERMINES WHEN PLASMA TOUCHES LIMITER
c     USE DOUBLE THE NO. OF THETA POINTS FOR INCREASED RESOLUTION
c
      ntheta_2pi = ntheta1
      nphi_plane = 1                 !Pick a phi plane (phi = 0 for now)

      rbdy(:ntheta3*2-1:2) = r1(ns*nphi_plane:ns*ntheta3*nphi_plane:ns*
     1   nphi_plane,0) + r1(ns*nphi_plane:ns*ntheta3*nphi_plane:ns*
     2   nphi_plane,1)
      zbdy(:ntheta3*2-1:2) = z1(ns*nphi_plane:ns*ntheta3*nphi_plane:ns*
     1   nphi_plane,0) + z1(ns*nphi_plane:ns*ntheta3*nphi_plane:ns*
     2   nphi_plane,1)
      rubdy(:ntheta3*2-1:2) = ru(ns*nphi_plane:ns*ntheta3*nphi_plane:ns*
     1   nphi_plane,0) + ru(ns*nphi_plane:ns*ntheta3*nphi_plane:ns*
     2   nphi_plane,1)
      zubdy(:ntheta3*2-1:2) = zu(ns*nphi_plane:ns*ntheta3*nphi_plane:ns*
     1   nphi_plane,0) + zu(ns*nphi_plane:ns*ntheta3*nphi_plane:ns*
     2   nphi_plane,1)

      IF (.not.lasym) THEN
!     FOR NOW, THIS ONLY WORKS FOR NZETA=1 (PHI=0 PLANE)
!     TO EXTEND TO OTHER PHI PLANES, MUST USE IREFLECT(NZETA)
      DO i = 1, ntheta_2pi - ntheta2
         rbdy(2*(ntheta2+i)-1) = rbdy(2*(ntheta1-ntheta2-i)+3)
      END DO
      DO i = 1, ntheta_2pi - ntheta2
         zbdy(2*(ntheta2+i)-1) = -zbdy(2*(ntheta1-ntheta2-i)+3)
      END DO
      DO i = 1, ntheta_2pi - ntheta2
         rubdy(2*(ntheta2+i)-1) = -rubdy(2*(ntheta1-ntheta2-i)+3)
      END DO
      DO i = 1, ntheta_2pi - ntheta2
         zubdy(2*(ntheta2+i)-1) = zubdy(2*(ntheta1-ntheta2-i)+3)
      END DO
      END IF

!
!     FIND EVEN INDEXED POINTS BY INTERPOLATION
!
      nthtot = 2*ntheta_2pi
      rbdy(nthtot) = .5_dp*(rbdy(1)+rbdy(nthtot-1))
      zbdy(nthtot) = .5_dp*(zbdy(1)+zbdy(nthtot-1))
      rubdy(nthtot) = .5_dp*(rubdy(1)+rubdy(nthtot-1))
      zubdy(nthtot) = .5_dp*(zubdy(1)+zubdy(nthtot-1))
      rbdy(2:(ntheta_2pi-1)*2:2) = .5_dp*(rbdy(3:ntheta_2pi*2-1:2) +
     1      rbdy(:ntheta_2pi*2-3:2))
      zbdy(2:(ntheta_2pi-1)*2:2) = .5_dp*(zbdy(3:ntheta_2pi*2-1:2) +
     1      zbdy(:ntheta_2pi*2-3:2))
      rubdy(2:(ntheta_2pi-1)*2:2) = .5_dp*(rubdy(3:ntheta_2pi*2-1:2)+
     1   rubdy(:ntheta_2pi*2-3:2))
      zubdy(2:(ntheta_2pi-1)*2:2) = .5_dp*(zubdy(3:ntheta_2pi*2-1:2)+
     1   zubdy(:ntheta_2pi*2-3:2))

      fshrink = 0.0_dp
      distmax = SUM(rbdy(:nthtot)**2) + SUM(zbdy(:nthtot)**2)
      fexpand = distmax
      iexpand = 0
      ishrink = 0
      ionlim = 0

      DO n = 1, nlim
         limpts = limitr(n)
         CALL cauchy (rbdy, zbdy, rubdy, zubdy, rlim(:,n), zlim(:,n),
     1      reslim(:,n), seplim(:,n), distmax, nthtot, limpts)

           DO i = 1,limpts
c       LIMITER POINT ON PLASMA
            IF( (ABS(reslim(i,n)-resdn).lt.eps) )then
!    .gt.      .or. (ABS(reslim(i,n)).gt.resup) )then
              ionlim = i
              nonlim = n
c       LIMITER POINT OUTSIDE PLASMA
            ELSE IF( reslim(i,n).lt.RESDN )then
              IF( seplim(i,n).le.fexpand )then
                fexpand = seplim(i,n)
                iexpand = i
                nexpand = n
              ENDIF
c       LIMITER POINT INSIDE PLASMA
            ELSE IF( reslim(i,n).ge.RESDN )then
              IF( seplim(i,n).gt.fshrink )then
                fshrink =  seplim(i,n)
                ishrink = i
                nshrink = n
              ENDIF
            ENDIF
          ENDDO
      END DO

c
c       LOGIC: IF THERE IS A LIMITER POINT INSIDE PLASMA, THEN MUST
c       SHRINK CONTOUR. OTHERWISE, IF THERE IS AT LEAST ONE LIMITER
c       POINT ON CONTOUR, AND ALL THE REST OUTSIDE, DO NOT CHANGE ANYTHING.
c       FINALLY, IF ALL LIMITER POINTS ARE OUTSIDE PLASMA, EXPAND PLASMA
c       TO OSCULATE WITH LIMITER
c
      IF (ishrink .gt. 0) THEN
         gphifac = -SQRT(fshrink)
         ilim0 = ishrink
         nlim0 = nshrink
      ELSE IF (ionlim .gt. 0) THEN
         gphifac = 0.
         ilim0 = ionlim
         nlim0 = nonlim
      ELSE IF (iexpand .gt. 0) THEN
         gphifac = SQRT(fexpand)
         ilim0 = iexpand
         nlim0 = nexpand
      ENDIF

      dlim_min = gphifac
      rlim_min = rlim(ilim0,nlim0)
      zlim_min = zlim(ilim0,nlim0)

c     overALL damping in time, rsfac/SQRT(rsfac) = SQRT(rsfac)
      gphifac = gphifac/r01(ns)
      IF (ABS(gphifac) .gt. 0.04*one)
     1   gphifac = 0.04*gphifac/ABS(gphifac)
      gphifac = 0.20*gphifac/SQRT(rsfac)

      END SUBROUTINE getlim


text:
#ifdef M3D_C1
!*******************************************************************************
!  @file korc_m3d_c1_interface.f90
!  @brief Interface for the m3d_c1 interpolation library.
!*******************************************************************************

MODULE korc_m3d_c1
  USE, INTRINSIC :: iso_c_binding
  USE korc_types
  USE korc_input
  USE korc_HDF5

  IMPLICIT NONE

  INTEGER (C_INT), PARAMETER :: FIO_SUCCESS        = 0
  INTEGER (C_INT), PARAMETER :: FIO_OUT_OF_BOUNDS  = 10002
  INTEGER (C_INT), PARAMETER :: FIO_NO_DATA        = 10006

  INTEGER (C_INT), PARAMETER :: FIO_M3DC1_SOURCE   = 3

  INTEGER (C_INT), PARAMETER :: FIO_TIMESLICE      = 1

  INTEGER (C_INT), PARAMETER :: FIO_SPECIES        = 3

  INTEGER (C_INT), PARAMETER :: FIO_ELECTRON       = 1
  INTEGER (C_INT), PARAMETER :: FIO_MAIN_ION       = -1

  INTEGER (C_INT), PARAMETER :: FIO_DENSITY        = 102
  INTEGER (C_INT), PARAMETER :: FIO_TEMPERATURE    = 103

  INTEGER (C_INT), PARAMETER :: FIO_ELECTRIC_FIELD = 1001
  INTEGER (C_INT), PARAMETER :: FIO_MAGNETIC_FIELD = 1003
  INTEGER (C_INT), PARAMETER :: FIO_VECTOR_POTENTIAL = 1002

  INTEGER (C_INT), PARAMETER :: FIO_TIME = 7001

  INTERFACE
     INTEGER (C_INT) FUNCTION fio_add_field(icfield, ifield, op, fac)      &
          BIND(C, NAME='fio_add_field')
       USE, INTRINSIC :: iso_c_binding

       IMPLICIT NONE

       INTEGER (C_INT), VALUE, INTENT(IN) :: icfield
       INTEGER (C_INT), VALUE, INTENT(IN) :: ifield
       INTEGER (C_INT), VALUE, INTENT(IN) :: op
       REAL (C_DOUBLE), VALUE, INTENT(IN) :: fac
     END FUNCTION fio_add_field
  END INTERFACE

  INTERFACE
     INTEGER (C_INT) FUNCTION fio_allocate_search_hint(isrc, hint)         &
          BIND(C, NAME='fio_allocate_search_hint')
       USE, INTRINSIC :: iso_c_binding

       IMPLICIT NONE

       INTEGER (C_INT), VALUE, INTENT(IN) :: isrc
       TYPE (C_PTR), INTENT(OUT)          :: hint
     END FUNCTION fio_allocate_search_hint
  END INTERFACE

  INTERFACE
     INTEGER (C_INT) FUNCTION fio_close_field(ifield)                      &
          BIND(C, NAME='fio_close_field')
       USE, INTRINSIC :: iso_c_binding

       IMPLICIT NONE

       INTEGER (C_INT), VALUE, INTENT(IN) :: ifield
     END FUNCTION fio_close_field
  END INTERFACE

  INTERFACE
     INTEGER (C_INT) FUNCTION fio_close_series(iseries)                    &
          BIND(C, NAME='fio_close_series')
       USE, INTRINSIC :: iso_c_binding

       IMPLICIT NONE

       INTEGER (C_INT), VALUE, INTENT(IN) :: iseries
     END FUNCTION fio_close_series
  END INTERFACE

  INTERFACE
     INTEGER (C_INT) FUNCTION fio_close_source(isrc)                       &
          BIND(C, NAME='fio_close_source')
       USE, INTRINSIC :: iso_c_binding

       IMPLICIT NONE

       INTEGER (C_INT), VALUE, INTENT(IN) :: isrc
     END FUNCTION fio_close_source
  END INTERFACE

  INTERFACE
     INTEGER (C_INT) FUNCTION fio_create_compound_field(ifield)            &
          BIND(C, NAME='fio_create_compound_field')
       USE, INTRINSIC :: iso_c_binding

       IMPLICIT NONE

       INTEGER (C_INT), INTENT(IN) :: ifield
     END FUNCTION fio_create_compound_field
  END INTERFACE

  INTERFACE
     INTEGER (C_INT) FUNCTION fio_deallocate_search_hint(isrc, hint)       &
          BIND(C, NAME='fio_deallocate_search_hint')
       USE, INTRINSIC :: iso_c_binding

       IMPLICIT NONE

       INTEGER (C_INT), VALUE, INTENT(IN) :: isrc
       TYPE (C_PTR), INTENT(INOUT)        :: hint
     END FUNCTION fio_deallocate_search_hint
  END INTERFACE

  INTERFACE
     INTEGER (C_INT) FUNCTION fio_eval_field(ifield, x, v, hint)         &
          BIND(C, NAME='fio_eval_field')
       USE, INTRINSIC :: iso_c_binding

       IMPLICIT NONE

       INTEGER (C_INT), VALUE, INTENT(IN) :: ifield
       REAL (C_DOUBLE), INTENT(IN)        :: x
       REAL (C_DOUBLE), INTENT(OUT)       :: v
       TYPE (C_PTR), VALUE, INTENT(IN)    :: hint
     END FUNCTION fio_eval_field
  END INTERFACE

  INTERFACE
     INTEGER (C_INT) FUNCTION fio_eval_field_deriv(ifield, x, v, hint)   &
          BIND(C, NAME='fio_eval_field_deriv')
       USE, INTRINSIC :: iso_c_binding

       IMPLICIT NONE

       INTEGER (C_INT), VALUE, INTENT(IN) :: ifield
       REAL (C_DOUBLE), INTENT(IN)        :: x
       REAL (C_DOUBLE), INTENT(OUT)       :: v
       TYPE (C_PTR), VALUE, INTENT(IN)    :: hint
     END FUNCTION fio_eval_field_deriv
  END INTERFACE

  INTERFACE
     INTEGER (C_INT) FUNCTION fio_eval_series(iseries, x, v)             &
          BIND(C, NAME='fio_eval_series')
       USE, INTRINSIC :: iso_c_binding

       IMPLICIT NONE

       INTEGER (C_INT), VALUE, INTENT(IN) :: iseries
       REAL (C_DOUBLE), INTENT(IN)        :: x
       REAL (C_DOUBLE), INTENT(OUT)       :: v
     END FUNCTION fio_eval_series
  END INTERFACE

  INTERFACE
     INTEGER (C_INT) FUNCTION fio_get_options(isrc)                        &
          BIND(C, NAME='fio_get_options')
       USE, INTRINSIC :: iso_c_binding

       IMPLICIT NONE

       INTEGER (C_INT), VALUE, INTENT(IN) :: isrc
     END FUNCTION fio_get_options
  END INTERFACE

  INTERFACE
     INTEGER (C_INT) FUNCTION fio_get_available_fields(isrc, n, f)         &
          BIND(C, NAME='fio_get_available_fields')
       USE, INTRINSIC :: iso_c_binding

       IMPLICIT NONE

       INTEGER (C_INT), VALUE, INTENT(IN)         :: isrc
       INTEGER (C_INT), INTENT(OUT)               :: n
       INTEGER (C_INT), DIMENSION(:), INTENT(OUT) :: f
     END FUNCTION fio_get_available_fields
  END INTERFACE

  INTERFACE
     INTEGER (C_INT) FUNCTION fio_open_source(itype, filename, handle)     &
          BIND(C, NAME='fio_open_source')
       USE, INTRINSIC :: iso_c_binding

       IMPLICIT NONE

       INTEGER (C_INT), VALUE, INTENT(IN)        :: itype
       CHARACTER (kind=C_CHAR,len=1), INTENT(IN) :: filename
       INTEGER (C_INT), INTENT(OUT)              :: handle
     END FUNCTION fio_open_source
  END INTERFACE

  INTERFACE
     INTEGER (C_INT) FUNCTION fio_get_field(isrc, type, handle)            &
          BIND(C, NAME='fio_get_field')
       USE, INTRINSIC :: iso_c_binding

       IMPLICIT NONE

       INTEGER (C_INT), VALUE, INTENT(in) :: isrc
       INTEGER (C_INT), VALUE, INTENT(IN) :: type
       INTEGER (C_INT), INTENT(INOUT)     :: handle
     END FUNCTION fio_get_field
  END INTERFACE

  INTERFACE
     INTEGER (C_INT) FUNCTION fio_set_int_option(iopt, v)                  &
          BIND(C, NAME='fio_set_int_option')
       USE, INTRINSIC :: iso_c_binding

       IMPLICIT NONE

       INTEGER (C_INT), VALUE, INTENT(in) :: iopt
       INTEGER (C_INT), VALUE, INTENT(in) :: v
     END FUNCTION fio_set_int_option
  END INTERFACE

  INTERFACE
     INTEGER (C_INT) FUNCTION fio_get_real_field_parameter(ifield, t, p) &
          BIND(C, NAME='fio_get_real_field_parameter')
       USE, INTRINSIC :: iso_c_binding

       IMPLICIT NONE

       integer(c_int), intent(in), value :: ifield, t
       real(c_double) :: p
     END FUNCTION fio_get_real_field_parameter
  END INTERFACE

CONTAINS

  SUBROUTINE initialize_m3d_c1(params, F, P, spp,init)

    TYPE(KORC_PARAMS), INTENT(INOUT)           :: params
    TYPE(FIELDS), INTENT(INOUT)                :: F
    TYPE(PROFILES), INTENT(INOUT)              :: P
    TYPE(SPECIES), DIMENSION(:), INTENT(INOUT) :: spp
    LOGICAL, INTENT(IN)  :: init

    INTEGER                                    :: ii
    INTEGER                                    :: pp
    INTEGER                                    :: status
    INTEGER                                    :: isrc
    real(c_double)  ::  time0,time1
    INTEGER (C_INT)                         :: M3D_C1_tmp

    if (init) then
       F%Efield = Efield
       F%PSIp_lim=PSIp_lim
       F%PSIp_0=PSIp_0
       F%ReInterp_2x1t=ReInterp_2x1t

       if (params%proceed) then
          call load_prev_iter(params)
          F%ind0_2x1t=params%prev_iter_2x1t+1
       else
          F%ind0_2x1t=ind0_2x1t
       end if

       F%ind_2x1t=F%ind0_2x1t

    end if

    status = fio_open_source(FIO_M3DC1_SOURCE,           &
         TRIM(params%magnetic_field_filename)            &
         // C_NULL_CHAR, isrc)

    status = fio_get_options(isrc)

    if (.not.F%ReInterp_2x1t) then
       status = fio_set_int_option(FIO_TIMESLICE, params%time_slice)
    else
       status = fio_set_int_option(FIO_TIMESLICE, F%ind_2x1t)
    end if

    status = fio_get_field(isrc, FIO_MAGNETIC_FIELD, F%M3D_C1_B)
    status = fio_get_field(isrc, FIO_ELECTRIC_FIELD, F%M3D_C1_E)
    status = fio_get_field(isrc, FIO_VECTOR_POTENTIAL, F%M3D_C1_A)

    status = fio_get_real_field_parameter(F%M3D_C1_B, FIO_TIME, time0)

    write(output_unit_write,*) 'FIO present time',time0

    if (F%ReInterp_2x1t) then
       status = fio_set_int_option(FIO_TIMESLICE, F%ind_2x1t+1)

       status = fio_get_field(isrc, FIO_MAGNETIC_FIELD, M3D_C1_tmp)

       status = fio_get_real_field_parameter(M3D_C1_tmp, FIO_TIME, time1)
       write(output_unit_write,*) 'FIO next time',time1

       status = fio_set_int_option(FIO_TIMESLICE, F%ind_2x1t)

       params%snapshot_frequency=time1-time0

       !write(6,*) 'snapshot_frequency',params%snapshot_frequency
       !write(6,*) 'dt',params%dt

       if (.not.init) then

          params%t_skip = FLOOR(params%snapshot_frequency/params%cpp%time/ &
               params%dt,ip)

       end if

    end if

    if (.not.F%Efield) F%M3D_C1_E=-1

    status = fio_set_int_option(FIO_SPECIES, FIO_ELECTRON);
    status = fio_get_field(isrc, FIO_DENSITY, P%M3D_C1_ne);
    status = fio_get_field(isrc, FIO_TEMPERATURE, P%M3D_C1_te);

    status = fio_set_int_option(FIO_SPECIES, FIO_MAIN_ION);
    status = fio_get_field(isrc, FIO_DENSITY, P%M3D_C1_ni);

    if (init) then
    !  Hardcode Bo to one for now until a better method of determining the a
    !  characteristic magnetic field value.
       F%Bo = 1.0
       F%Eo = 1.0
       F%Ro = 1.0
       F%Zo = 1.0

       do ii = 1, params%num_species

          do pp = 1, spp(ii)%ppp
             status = fio_allocate_search_hint(isrc, spp(ii)%vars%hint(pp))
             !spp(ii)%vars%hint(pp)=c_null_ptr
          end do

          spp(ii)%vars%cart = .false.
       end do
    end if

    if (params%mpi_params%rank .EQ. 0) then
       write(output_unit_write,*) 'Calculate B',F%M3D_C1_B
       write(output_unit_write,*) 'Calculate E',F%M3D_C1_E
       write(output_unit_write,*) 'Calculate A',F%M3D_C1_A
       write(output_unit_write,*) 'Calculate ne',P%M3D_C1_ne
       write(output_unit_write,*) 'Calculate Te',P%M3D_C1_te
       write(output_unit_write,*) 'Calculate ni',P%M3D_C1_ni
    end if

  END SUBROUTINE initialize_m3d_c1

  SUBROUTINE initialize_m3d_c1_imp(params,F,P,num_imp,init)

    TYPE(KORC_PARAMS), INTENT(IN)           :: params
    TYPE(FIELDS), INTENT(IN)                :: F
    TYPE(PROFILES), INTENT(INOUT)              :: P
    INTEGER, INTENT(IN)				:: num_imp
    LOGICAL, INTENT(IN)  :: init

    INTEGER                                    :: ii
    INTEGER                                    :: status
    INTEGER                                    :: isrc
    INTEGER,ALLOCATABLE,DIMENSION(:)          :: Zo
    INTEGER  :: A,Zo1

    status = fio_open_source(FIO_M3DC1_SOURCE,           &
         TRIM(params%magnetic_field_filename)            &
         // C_NULL_CHAR, isrc)

    status = fio_get_options(isrc)

    if (.not.F%ReInterp_2x1t) then
       status = fio_set_int_option(FIO_TIMESLICE, params%time_slice)
    else
       status = fio_set_int_option(FIO_TIMESLICE, F%ind_2x1t)
    end if

    if (init) ALLOCATE(P%M3D_C1_nimp(num_imp))

    !write(6,*) size(params%Zj)
    !write(6,*) size(params%Zj(ubound(params%Zj)))
    Zo=int(params%Zj(ubound(params%Zj)))
    Zo1=Zo(1)

    if (Zo1.eq.18) then
       A=40
    else if (Zo1.eq.6) then
       A=12
    end if

    do ii=1,num_imp
       status = fio_set_int_option(FIO_SPECIES, &
            FIO_MAKE_SPECIES(A, Zo1, Zo1+1-ii));
       status = fio_get_field(isrc, FIO_DENSITY, P%M3D_C1_nimp(ii));
    end do

    if (params%mpi_params%rank .EQ. 0) then
       do ii=1,num_imp
          write(output_unit_write,*) 'Calculate nimp_',ii,P%M3D_C1_nimp(ii)
       end do
    end if

  END SUBROUTINE initialize_m3d_c1_imp

  FUNCTION FIO_MAKE_SPECIES(m, p, e)
    INTEGER, INTENT(IN) :: m
    INTEGER, INTENT(IN) :: p
    INTEGER, INTENT(IN) :: e
    INTEGER(C_INT) :: FIO_MAKE_SPECIES

    FIO_MAKE_SPECIES = e + p*256 + (m-p)*65536
  END FUNCTION FIO_MAKE_SPECIES

END MODULE korc_m3d_c1
#endif

max_stars_repo_name: ORNL-Fusion/KORC

max_issues_repo_name: ORNL-Fusion/KORC

max_forks_repo_name: ORNL-Fusion/KORC


text:
c     Subroutine to find distinct record numbers
c     AJ_Kettle, 02Oct2019

      SUBROUTINE find_distinct_recnum_gsom(
     +   l_channel,s_avec_recordnumber,
     +   l_distinct_recnum,s_dvec_recordnumber)

      IMPLICIT NONE
c************************************************************************
c     Variables passed into program

      INTEGER             :: l_channel
      INTEGER             :: l_distinct_recnum

      CHARACTER(LEN=*)    :: s_avec_recordnumber(l_channel)
      CHARACTER(LEN=*)    :: s_dvec_recordnumber(l_channel)
c*****
c     Variables used in subroutine

      INTEGER             :: i,j,k,ii,jj,kk

c     intermendiate info on valid recnums
      INTEGER             :: l_validvec
      CHARACTER(LEN=2)    :: s_validvec_recnum(l_channel)

      INTEGER             :: i_dvec_count(l_channel)
c************************************************************************
c      print*,'just entered find_distinct_recnum_gsom'
c*****
c     Initialize
      l_distinct_recnum=0
      DO i=1,l_channel
       s_dvec_recordnumber(i)=''
       i_dvec_count(i)       =0
      ENDDO
c*****
      ii=0

      DO i=1,l_channel
       IF (LEN_TRIM(s_avec_recordnumber(i)).GT.0) THEN
        ii=ii+1
        s_validvec_recnum(ii)=s_avec_recordnumber(i)
       ENDIF
      ENDDO

      l_validvec=ii
c*****
c     Case of no valid numbers
      IF (l_validvec.EQ.0) THEN
       GOTO 10
      ENDIF

c      print*,'l_validvec=',l_validvec
c      print*,'s_validvec_recnum=',(s_validvec_recnum(i),i=1,l_validvec)
c*****
c     Case of single record_num

      IF (l_validvec.EQ.1) THEN
       l_distinct_recnum     =1
       s_dvec_recordnumber(1)=s_validvec_recnum(1)
       i_dvec_count(1)       =1
       GOTO 10
      ENDIF
c*****
c     For greater numbers find distinct record number

      l_distinct_recnum=0
      l_distinct_recnum=l_distinct_recnum+1
      s_dvec_recordnumber(l_distinct_recnum)=s_validvec_recnum(1)
      i_dvec_count(l_distinct_recnum)=i_dvec_count(l_distinct_recnum)+1

      DO i=2,l_validvec
       DO j=1,l_distinct_recnum
        IF (TRIM(s_validvec_recnum(i)).EQ.TRIM(s_dvec_recordnumber(j)))
     +     THEN
         i_dvec_count(j)=i_dvec_count(j)+1
         GOTO 20
        ENDIF
       ENDDO

c      If here then must augment
       l_distinct_recnum=l_distinct_recnum+1
       s_dvec_recordnumber(l_distinct_recnum)=s_validvec_recnum(1)
       i_dvec_count(l_distinct_recnum)=i_dvec_count(l_distinct_recnum)+1

 20    CONTINUE
      ENDDO

c      print*,'l_distinct_recnum=',l_distinct_recnum
c      print*,'s_dvec_recordnumber=',
c     +   (s_dvec_recordnumber(i),i=1,l_distinct_recnum)
c      print*,'i_dvec_count=',
c     +   (i_dvec_count(i),i=1,l_distinct_recnum)

c*****
 10   CONTINUE

c      print*,'just leaving find_distinct_recnum_gsom'

c      STOP 'find_distinct_recnum_gsom'

      RETURN
      END

max_stars_repo_name: glamod/glamod-nuim

max_issues_repo_name: glamod/glamod-nuim

max_forks_repo_name: glamod/glamod-nuim


text:
C  *********************************************************
C  *                                                       *
C  *    TEST NUMBER: 09.01.09.01/04                        *
C  *    TEST TITLE : Error 405 handled by user             *
C  *                                                       *
C  *    PHIGS Validation Tests, produced by NIST           *
C  *                                                       *
C  *********************************************************

      CALL INITGL ('09.01.09.01/04')
      CALL E09405 (1)
      CALL ENDIT
      END