!***********************************************************************************************************************************
! UNIT:
!
!  (MODULE) MR_DEF_COORS
!
! PURPOSE:
!
!   TO DEFINE THE COORDINATES.
!
! DEFINITION OF VARIABLES:
!
!                            XYUV    :    COORDINATES AT CENTROIDS
!                      XYUU, XYVV    :    COORDINATES AT U- AND V-, RESPECTIVELY, FACIAL NODES
!                            XYOO    :    COORDINATES AT CORNER NODES
!
! RECORD OF REVISIONS:
!
!      DATE       |    PROGRAMMER    |    DESCRIPTION OF CHANGE
!      ====       |    ==========    |    =====================
!   20XX-XX-XX    |     DR. HYDE     |    ORIGINAL CODE.
!
!***********************************************************************************************************************************
  MODULE MR_DEF_COORS

    USE MR_KINDS

    IMPLICIT NONE

    REAL   (XYRD_KIND) , ALLOCATABLE , DIMENSION(:,:,:  ) :: XYUV
    REAL   (XYRD_KIND) , ALLOCATABLE , DIMENSION(:,:,:  ) :: XYUU
    REAL   (XYRD_KIND) , ALLOCATABLE , DIMENSION(:,:,:  ) :: XYVV
    REAL   (XYRD_KIND) , ALLOCATABLE , DIMENSION(:,:,:  ) :: XYOO

  END MODULE MR_DEF_COORS


text:
C
C**************************  NORMS  ************************************
C
C  calculate the nonlinear convergence and residual errors, using
C  the L2 and infinity norms
C
C***********************************************************************
C
      SUBROUTINE NORMS(N,IKMAX,PIKMAX,PINF,PL2,FINF,FL2,PNEW,POLD,TNOTI)
C
      IMPLICIT  NONE
      INTEGER   K
      INTEGER   N,IKMAX
      REAL*8    DIFF,DIFABS,FABS
      REAL*8    PIKMAX,PINF,PL2,FINF,FL2
      REAL*8    PNEW(*),POLD(*),TNOTI(*)
C
      PL2=0.0D0
      FL2=0.0D0
      PINF=0.0D0
      FINF=0.0D0
      DO K=1,N
         DIFF=PNEW(K)-POLD(K)
         DIFABS=DABS(DIFF)
         PL2=PL2 + DIFF*DIFF
         FL2=FL2 + TNOTI(K)*TNOTI(K)
         IF (DIFABS .GE. PINF) THEN
            PINF=DIFABS
            PIKMAX=DIFF
            IKMAX=K
         END IF
         FABS=DABS(TNOTI(K))
         IF (FABS .GE. FINF) THEN
            FINF=FABS
         END IF
      END DO
      PL2=DSQRT(PL2)
      FL2=DSQRT(FL2)
C
      RETURN
      END

ext: f


text:
      subroutine mkspl2zb(fun,x,nx,th,nth,fspl,nf2,
     >   ibcxmin,bcxmin,ibcxmax,bcxmax,
     >   wk,inwk,ilinx,ilinth,ier)
C
C  create a bicubic periodic spline with knots at grid points and
C  function values from the callable function `fun' passed.
C
C  use bpsplinb to setup the spline coefficients
C
C  periodic boundary condition for theta grid;
C  boundary condition data may be specified at x(1) and x(nx)
C     ibcxmin=0 ==> "not a knot" boundary condition, cf cubspl.for
C     ibcxmin=1 ==> match slope, bcxmin(ith) gives df/dx at x(1),th(ith).
C     ibcxmin=2 ==> match 2nd derivative, given at x(1),th(ith) by bcxmin(ith)
C
C     ibcxmax,bcxmax have analogous interpretation -- at x(nx)
C
      external fun                      ! passed real function(x,th)
      real x(nx)                        ! x coordinate array
      real th(nth)                      ! th coordinate array
C
      real fspl(4,nf2,nth)              ! function data / spline coeff array
      real wk(inwk)                     ! workspace for bpsplinb
C
      integer ibcxmin                   ! boundary condition indicator
      real bcxmin(nth)                  ! boundary condition data
      integer ibcxmax                   ! boundary condition indicator
      real bcxmax(nth)                  ! boundary condition data
C
      integer ilinx                     ! output =1 if x(...) evenly spaced
      integer ilinth                    ! output =1 if th(...) evenly spaced
C
      integer ier                       ! completion code from bpspline
C
C----------------------------
C
      ier=0
      if(nf2.lt.nx) then
         write(6,'('' ?mkspl2pb -- array dim error, nf2.lt.nx'')')
         ier=1
      endif
C
      do ix=1,nx
         do ith=1,nth
            fspl(1,ix,ith)=fun(x(ix),th(ith))
         enddo
      enddo
C
      call mkbicubw(x,nx,th,nth,fspl,nf2,
     >   ibcxmin,bcxmin,ibcxmax,bcxmax,
     >   -1,zdummy,-1,zdummy,
     >   wk,inwk,ilinx,ilinth,ier)
C
      return
      end

max_stars_repo_name: ORNL-Fusion/RFSciDAC-testing

max_issues_repo_name: ORNL-Fusion/RFSciDAC-testing

max_forks_repo_name: ORNL-Fusion/RFSciDAC-testing


! { dg-do run }
!
! Contributed by mrestelli@gmail.com
! Adapated by Andre Vehreschild  <vehre@gcc.gnu.org>
! Test that fix for PR69296 is working.

program p
 implicit none

 integer :: j, a(2,6), i(3,2)

  a(1,:) = (/ (     j , j=1,6) /)
  a(2,:) = (/ ( -10*j , j=1,6) /)

  i(:,1) = (/ 1 , 3 , 5 /)
  i(:,2) = (/ 4 , 5 , 6 /)

  associate( ai => a(:,i(:,1)) )
    if (any(shape(ai) /= [2, 3])) call abort()
    if (any(reshape(ai, [6]) /= [1 , -10, 3, -30, 5, -50])) call abort()
  end associate

end program p


text:
!***************************************************************************************************
! xnet_flux.f90 10/18/17
! This file containes the data structures and routines for calculating net fluxes.
!***************************************************************************************************

Module xnet_flux
  !-------------------------------------------------------------------------------------------------
  ! This module contains the data for calculating the net flux of each matched reaction pair.
  !-------------------------------------------------------------------------------------------------
  Use xnet_types, Only: dp
  Implicit None
  Real(dp), Allocatable :: flx(:,:), flx_int(:,:)
  Real(dp), Allocatable :: dcflx(:,:)
  Integer, Allocatable  :: ifl_orig(:), ifl_term(:)
  !$omp threadprivate(flx,flx_int)

Contains

  Subroutine flux_init
    !-----------------------------------------------------------------------------------------------
    ! This routine allocates the flux arrays and determines the double counting factors necessary for
    ! reactions with identical reactants.
    !-----------------------------------------------------------------------------------------------
    Use xnet_controls, Only: idiag, lun_diag, nzbatchmx
    Use xnet_match, Only: mflx, nflx
    Use xnet_types, Only: dp
    Implicit None

    ! Local variables
    Real(dp), Parameter :: fact(5) = 1.0 / (/ 1.0, 2.0, 6.0, 24.0, 4.0 /)
    Integer :: i, j, k, countf(mflx), countr(mflx)

    ! Allocate flux arrays
    Allocate (dcflx(2,mflx),ifl_orig(mflx),ifl_term(mflx))

    ! Check double counting of reactants for each flux
    Do i = 1, mflx
      countf(i) = 1 + count(nflx(2:4,i) /= 0 .and. nflx(2:4,i) == nflx(1:3,i))
      countr(i) = 1 + count(nflx(6:8,i) /= 0 .and. nflx(6:8,i) == nflx(5:7,i))

      ! countr=3 can result from a+a+a+b or a+a+b+b, which has different factor
      If ( countf(i) == 3 .and. nflx(2,i) /= nflx(3,i) ) countf(i) = 5
      If ( countr(i) == 3 .and. nflx(6,i) /= nflx(7,i) ) countr(i) = 5
      dcflx(1,i) = fact(countf(i))
      dcflx(2,i) = fact(countr(i))

      ! Determine flux origin, ifl_orig, and terminus, ifl_term
      ifl_orig(i) = nflx(count(nflx(1:4,i) /= 0 ),    i)
      ifl_term(i) = nflx(count(nflx(5:8,i) /= 0 ) + 4,i)
    EndDo

    !$omp parallel default(shared)
    Allocate (flx(mflx,nzbatchmx),flx_int(mflx,nzbatchmx))
    flx_int = 0.0
    If ( idiag >= 5 ) Then
      Write(lun_diag,"(a)") 'Flux Init'
      Write(lun_diag,"(12i5,2f6.3)") &
        & ((nflx(j,i),j=1,8),ifl_orig(i),ifl_term(i),countf(i),countr(i),(dcflx(k,i),k=1,2),i=1,mflx)
    EndIf
    !$omp end parallel

    Return
  End Subroutine flux_init

  Subroutine flux(mask_in)
    !-----------------------------------------------------------------------------------------------
    ! This routine calculates the net flux of each matched reaction pair. Positive fluxes flow in the
    ! direction of increasing binding.
    !-----------------------------------------------------------------------------------------------
    Use nuclear_data, Only: nname
    Use reaction_data, Only: csect1, csect2, csect3, csect4, n1i, n2i, n3i, n4i, nreac
    Use xnet_abundances, Only: yt
    Use xnet_conditions, Only: tdel
    Use xnet_controls, Only: idiag, lun_diag, ymin, nzbatchmx, szbatch, lzactive
    Use xnet_match, Only: ifl1, ifl2, ifl3, ifl4, iwflx, mflx, nflx
    Use xnet_types, Only: dp
    Implicit None

    ! Optional variables
    Logical, Optional, Target, Intent(in) :: mask_in(:)

    ! Local variables
    Integer :: i, k, ifl, idcfl, izb, izone
    Real(dp) :: flt, flxmin
    Logical, Pointer :: mask(:)

    If ( present(mask_in) ) Then
      mask => mask_in(:)
    Else
      mask => lzactive(:)
    EndIf
    If ( .not. any(mask(:)) ) Return

    Do izb = 1, nzbatchmx
      If ( mask(izb) ) Then

        ! Get fluxes
        flx(:,izb) = 0.0
        Do i = 1, nreac(1)
          ifl = abs(ifl1(i))
          idcfl = nint(1.5_dp - sign(0.5_dp,real(ifl1(i),dp)))
          flt = csect1(i,izb) * yt(n1i(1,i),izb)
          flx(ifl,izb) = flx(ifl,izb) + dcflx(idcfl,ifl)*sign(flt,real(ifl1(i),dp))
        EndDo
        Do i = 1, nreac(2)
          ifl = abs(ifl2(i))
          idcfl = nint(1.5_dp - sign(0.5_dp,real(ifl2(i),dp)))
          flt = csect2(i,izb) * yt(n2i(1,i),izb) * yt(n2i(2,i),izb)
          flx(ifl,izb) = flx(ifl,izb) + dcflx(idcfl,ifl)*sign(flt,real(ifl2(i),dp))
        EndDo
        Do i = 1, nreac(3)
          ifl = abs(ifl3(i))
          idcfl = nint(1.5_dp - sign(0.5_dp,real(ifl3(i),dp)))
          flt = csect3(i,izb) * product(yt(n3i(1:3,i),izb))
          flx(ifl,izb) = flx(ifl,izb) + dcflx(idcfl,ifl)*sign(flt,real(ifl3(i),dp))
        EndDo
        Do i = 1, nreac(4)
          ifl = abs(ifl4(i))
          idcfl = nint(1.5_dp - sign(0.5_dp,real(ifl4(i),dp)))
          flt = csect4(i,izb) * product(yt(n4i(1:4,i),izb))
          flx(ifl,izb) = flx(ifl,izb) + dcflx(idcfl,ifl)*sign(flt,real(ifl4(i),dp))
        EndDo

        ! Since we cut off abundances less than ymin, we should similarly cut off small fluxes
        flxmin = ymin / tdel(izb)
        Where ( abs(flx(:,izb)) > flxmin )
          flx_int(:,izb) = flx_int(:,izb) + flx(:,izb)*tdel(izb)
        EndWhere
        If ( idiag >= 5 ) Then
          izone = izb + szbatch - 1
          Write(lun_diag,"(a5,i5,es10.3)") 'Flux',izone,flxmin
          Write(lun_diag,"(i5,9a5,i5,2es23.15)") &
            & (k,(nname(nflx(i,k)),i=1,4),' <-> ',(nname(nflx(i,k)),i=5,8),iwflx(k),flx(k,izb),flx_int(k,izb),k=1,mflx)
        EndIf
      EndIf
    EndDo

    Return
  End Subroutine flux

  Subroutine flux_check(mask_in)
    !-----------------------------------------------------------------------------------------------
    ! This routine compares the fluxes to the changes in abundance.
    !-----------------------------------------------------------------------------------------------
    Use nuclear_data, Only: ny, nname
    Use xnet_abundances, Only: y, ydot, yo
    Use xnet_conditions, Only: tdel
    Use xnet_controls, Only: idiag, lun_diag, nzbatchmx, szbatch, lzactive
    Use xnet_match, Only: mflx, nflx
    Use xnet_types, Only: dp
    Implicit None

    ! Optional variables
    Logical, Optional, Target, Intent(in) :: mask_in(:)

    ! Local variables
    Real(dp) :: flx_test(0:ny)
    Integer :: i, izb, izone
    Logical, Pointer :: mask(:)

    If ( present(mask_in) ) Then
      mask => mask_in(:)
    Else
      mask => lzactive(:)
    EndIf
    If ( .not. any(mask(:)) ) Return

    Do izb = 1, nzbatchmx
      If ( mask(izb) ) Then

        ! Calculate flx_test = dy + dt * (flux_out - flux_in)
        flx_test(0) = 0.0
        !flx_test(1:ny) = -ydot(:,izb) * tdel(izb) ! Use only with convc>0
        flx_test(1:ny) = yo(:,izb) - y(:,izb)
        Do i = 1, mflx
          flx_test(nflx(1:4,i)) = flx_test(nflx(1:4,i)) - flx(i,izb)*tdel(izb)
          flx_test(nflx(5:8,i)) = flx_test(nflx(5:8,i)) + flx(i,izb)*tdel(izb)
        EndDo
        If ( idiag >= 5 ) Then
          izone = izb + szbatch - 1
          Write(lun_diag,"(a,2i5,es11.3)") "Flux Check",izone,mflx,tdel(izb)
          Write(lun_diag,"(a5,4es11.3)") &
            & (nname(i),y(i,izb),yo(i,izb)-y(i,izb),ydot(i,izb)*tdel(izb),flx_test(i),i=1,ny)
        EndIf
      EndIf
    EndDo

    Return
  End Subroutine flux_check

End Module xnet_flux

max_stars_repo_name: vishaltiwari/XNet

max_issues_repo_name: vishaltiwari/XNet

max_forks_repo_name: vishaltiwari/XNet


text:
South Hall is address(home, 38.541203, 121.748132) of Advising Services, the Internship and Career Center, Unitrans Business Office, Specialized Transportation Services (STS), and the Student Recruitment and Retention Center.

Confusingly, many key parts of Advising Services are located across campus in 160 South Silo. Most notably, all prehealth, prebusiness, and prelaw advising and services have moved to South Silo. See the Advising Services page for more details.

 First Floor Offices:
  The First Resort  Room 116
  The Academic Peer Advising Program  Room 116
  Graduate Letter Service, a.k.a. GLS  Room 111
  Career and Graduate Study Resource Center Room 114

http://content.cdlib.org/ark:/13030/tf0489n4t1/?querydavis&brandcalisphere South Hall was one of the three original dormitories (with North Hall and West Hall) on the Campus UC Davis campus, and was built in early 1912, shortly after North Hall.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
-v alu.v
+ncaccess+rwc
alutb.v

max_stars_repo_name: eda-ricercatore/PicenoDecoders

max_issues_repo_name: eda-globetrotter/MarcheProcessor

max_forks_repo_name: eda-globetrotter/MarcheProcessor


text:
!Subprograms for doing task-specific input (read namelists, data)
MODULE Input

  ! #DES: Subprograms for doing task-specific input (read namelist, input files)
  ! #DES: Also stores the input data and provides access to it for other modules

  IMPLICIT NONE

!  PRIVATE
!  PUBLIC :: ReadInput, DeallocateInputArrays

  ! #DES: The namelist settings.nml needs to contain the variables below, and must match the data in the energy files.
  ! Nstates is the number of states in the EVB potential.
  ! The particular states used to generate the mapping potential are stateA and stateB.
  ! NFepSteps is the number of steps taken from stateA to stateB in the FEP simulation.
  ! Ntimesteps is the number of timesteps in each file.
  ! Nskip tells the program to ignore the first Nskip timesteps in each file.
  ! The final state for LRA study is stepTS.
  ! nBins sets the number of bins to use in the FEP/US histogram.
  ! minPop is the minimum number of points that must be in a FEP/US histogram bin to be considered usable.
  ! fileBase is the part of the filename preceding the FEP index (in I3 format)
  ! temperature is the thermodynamic temperature at which the simulation was run in Kelvin
  ! rcCoeffA and rcCoeffB are used to construct a reaction coordinate from the energy difference
  ! RC = (rcCoeffA * E_A) + (rcCoeffB * E_B)

  NAMELIST /settings/ nStates, stateA, stateB, nFepSteps, nSkip, stepTS, nBins, minPop, fileBase, trajecBase, &
                      temperature, rcCoeffA, rcCoeffB, useEnergyGapCoupling, dGTS, dGPS, outDir
  INTEGER :: nFepSteps, nStates, nSkip, stateA, stateB, stepTS, nBins, minPop
  REAL(8) :: temperature, rcCoeffA, rcCoeffB, dGTS, dGPS
  LOGICAL :: useEnergyGapCoupling
  CHARACTER(100) :: fileBase, trajecBase, outDir

  ! #DES: The stateEnergy and coeffs arrays hold all of the simulation data (irrespective of nSkip, nTimesteps)
  REAL(8), ALLOCATABLE :: stateEnergy(:,:,:,:), coeffs(:,:,:)
  ! The trajectory array holds all the Cartesian coordinates of the set of FEP simulations - Single Precision cos of the dcd format.
  REAL(4), ALLOCATABLE :: trajectory(:,:,:,:)
  ! #DES: These arrays hold the EVB parameters
  REAL(8), ALLOCATABLE :: alpha(:), couplingConstant(:,:), couplingGaussExpFactor(:,:), couplingExpExpFactor(:,:)
  REAL(8) :: beta
  LOGICAL :: targetsPresent = .FALSE.

  ! This data is used to determine what data in stateEnergy and coeffs is to be used in the FEP/US procedure
  ! The block mask(fep,skip(fep)+1:nTimesteps(fep)) must be set to .TRUE. initially, the rest to .FALSE.

  INTEGER, ALLOCATABLE :: nTimesteps(:) ! nFepSteps - number of actual timesteps in each FEP simulation
  INTEGER, ALLOCATABLE :: skip(:)       ! nFepSteps - How many steps to ignore at the start of each step; allows a per-simulation equilibration period
  LOGICAL, ALLOCATABLE :: mask(:,:)     ! nFepSteps x nTimesteps - pointwise mask for simulation data; turn points on/off
  INTEGER :: maxTimesteps               ! Length of the longest FEP simulation (largest of nTimesteps(:)); used to allocate arrays

  ! These will eventually have to be inputs unless a logical switch for Q data is TRUE
  CHARACTER(20), ALLOCATABLE :: energyNames(:)
  INTEGER :: nEnergyTypes

  CONTAINS

    SUBROUTINE AllocateInputArrays

      ! #DES: Allocate and zero the input arrays based on the supplied namelist inputs

      IMPLICIT NONE

      ALLOCATE(stateEnergy(maxTimesteps,nFepSteps,nStates,nEnergyTypes)); stateEnergy = 0.0d0
      ALLOCATE(skip(nFepSteps));                                          skip(:) = 0
      ALLOCATE(mask(nFepSteps,maxTimesteps));                             mask(:,:) = .TRUE.
      ALLOCATE(coeffs(maxTimesteps,nFepSteps,nStates));                   coeffs(:,:,:) = 0.0d0

      ALLOCATE(alpha(nStates));                                           alpha(:)                    = 0.0d0
      ALLOCATE(couplingConstant(nStates,nStates));                        couplingConstant(:,:)       = 0.0d0
      ALLOCATE(couplingGaussExpFactor(nStates,nStates));                  couplingGaussExpFactor(:,:) = 0.0d0
      ALLOCATE(couplingExpExpFactor(nStates,nStates));                    couplingExpExpFactor(:,:)   = 0.0d0

    END SUBROUTINE AllocateInputArrays

!*

    SUBROUTINE DeallocateInputArrays

      ! #DES: Free the memory associated with the input arrays

      IMPLICIT NONE

      IF (ALLOCATED(stateEnergy))             DEALLOCATE(stateEnergy)
      IF (ALLOCATED(coeffs))                  DEALLOCATE(coeffs)
      IF (ALLOCATED(alpha))                   DEALLOCATE(alpha)
      IF (ALLOCATED(couplingConstant))        DEALLOCATE(couplingConstant)
      IF (ALLOCATED(couplingGaussExpFactor))  DEALLOCATE(couplingGaussExpFactor)
      IF (ALLOCATED(couplingExpExpFactor))    DEALLOCATE(couplingExpExpFactor)
      IF (ALLOCATED(skip))                    DEALLOCATE(skip)
      IF (ALLOCATED(nTimesteps))              DEALLOCATE(nTimesteps)
      IF (ALLOCATED(mask))                    DEALLOCATE(mask)
      IF (ALLOCATED(trajectory))              DEALLOCATE(trajectory)

    END SUBROUTINE DeallocateInputArrays

!*

    SUBROUTINE ReadInput(readCoords)

      ! #DES: Master routine for reading all needed input from files

      USE Log, ONLY : logUnit
      IMPLICIT NONE
      LOGICAL, INTENT(IN) :: readCoords

      CALL ProcessNamelist(logUnit)
      CALL ReadEnergyTypes(logUnit)
      CALL ReadSimulationLengths()
      CALL WriteSimulationLengths(logUnit)
      CALL AllocateInputArrays()
      CALL ReadEVBParameters(logUnit)
      CALL PrintEVBParameters(logUnit)
      CALL ReadConvergenceParameters(logUnit)
      CALL WriteConvergenceParameters(logUnit)
      CALL ReadFormattedEnergyFiles(logUnit)
      IF (readCoords .EQV. .TRUE.) CALL ReadCoordinateFiles(logUnit)

    END SUBROUTINE ReadInput

!*

    SUBROUTINE ReadSimulationLengths

      ! #DES: Determine the number of timesteps in each simulation

      USE FileIO, ONLY : fileLengths
      IMPLICIT NONE
      CHARACTER(500) :: fileNames(nFepSteps)
      INTEGER :: step

      ALLOCATE(nTimeSteps(nFepSteps))
      CALL CreateFileNames('en',fileBase,fileNames)
      nTimesteps(:) = fileLengths(fileNames)
      nTimeSteps(:) = nTimesteps(:) / nStates
      DO step = 1, nFepSteps
        IF (ALLOCATED(mask)) mask(step,nTimesteps(step)+1:maxTimesteps) = .FALSE.
      ENDDO
      maxTimesteps = MAXVAL(nTimesteps)

      IF (nSkip < 0 .OR. nSkip >= maxTimesteps-1) THEN
        WRITE(*,*) "ReadSimulationLength - Illegal value of nSkip: ", nSkip
        IF (nSkip < 0) STOP "nSkip < 0"
        IF (nSkip >= maxTimesteps-1) STOP "nSkip >= maxTimesteps - 1"
      ENDIF
      IF (minPop < 1 .OR. minPop > maxTimesteps)  STOP ": ReadSimulationLength: Illegal value of minPop"

    END SUBROUTINE ReadSimulationLengths

!*

    SUBROUTINE WriteSimulationLengths(logUnit)

      IMPLICIT NONE
      INTEGER, INTENT(IN) :: logUnit
      CHARACTER(500) :: fileNames(nFepSteps)
      INTEGER :: step

      CALL CreateFileNames('en',fileBase,fileNames)
      WRITE(logUnit,'(A)') "* Simulation Lengths *"; WRITE(logUnit,*)
      DO step = 1, nFepSteps
        WRITE(logUnit,'(A,I8)') TRIM(ADJUSTL(fileNames(step))), nTimesteps(step)
      ENDDO
      WRITE(logUnit,*)
      WRITE(logUnit,'(A,I8)') "Longest Simulation Length: ", maxTimesteps
      WRITE(logUnit,*)

    END SUBROUTINE WriteSimulationLengths

!*

    SUBROUTINE WriteConvergenceParameters(logUnit)

      IMPLICIT NONE
      INTEGER, INTENT(IN) :: logUnit
      CHARACTER(500) :: fileNames(nFepSteps)
      INTEGER :: step

      WRITE(logUnit,'(A)') "Convergence parameters:"
      WRITE(logUnit,*)
      WRITE(logUnit,'(A)') "Energy File  Fep Step   Skip"
      CALL createFileNames('en',fileBase,fileNames)
      DO step = 1, nFepSteps
        WRITE(logUnit,'(A,2I8)')  TRIM(ADJUSTL(fileNames(step))), step, skip(step)
      ENDDO
      WRITE(logUnit,*)

    END SUBROUTINE WriteConvergenceParameters

!*

    SUBROUTINE ReadConvergenceParameters(logUnit)

      USE FileIO, ONLY : OpenFile, CloseFile

      IMPLICIT NONE
      INTEGER, INTENT(IN) :: logUnit
      INTEGER, PARAMETER :: prmUnit = 16
      INTEGER :: step
      LOGICAL :: fileExists
      CHARACTER(100) :: discard

      WRITE(logUnit,'(A)', ADVANCE='NO') "Determining FEP Simulation Convergence Parameters... "
      CALL OpenFile(prmUnit,"skip.prm","read",success=fileExists)
      IF (fileExists) THEN
        DO step = 1, nFepSteps
          READ(prmUnit,*) discard, skip(step)
          IF (skip(step) > nTimesteps(step)) STOP "ERROR: nSkip > nTimesteps (ReadConvergenceParameters)"
          mask(step,1:skip(step)) = .FALSE.
        ENDDO
        CALL CloseFile(prmUnit)
        WRITE(logUnit,'(A)') "data read from skip.prm"
        WRITE(logUnit,*)
      ELSE
        skip(:) = nSkip
        mask(:,1:nSkip) = .FALSE.
        CALL WriteConvergenceParamFile()
        WRITE(logUnit,'(A)') "no convergence parameter file provided; template written to skip.prm"
        WRITE(logUnit,*)
      ENDIF

    END SUBROUTINE ReadConvergenceParameters

!*

    SUBROUTINE WriteConvergenceParamFile()

      ! #DES: Write an empty convergence parameter file to skip.prm for the number of FEP simulations

      USE FileIO, ONLY : OpenFile, CloseFile

      IMPLICIT NONE
      INTEGER, PARAMETER :: prmUnit = 16
      CHARACTER(500) :: fileNames(nFepSteps)
      INTEGER :: step

      CALL OpenFile(prmUnit,"skip.prm","write")
      CALL CreateFileNames('en',fileBase,fileNames)

      DO step = 1, nFepSteps
        WRITE(prmUnit,'(A,I8)') TRIM(ADJUSTL(fileNames(step))), skip(step)
      ENDDO

      CALL CloseFile(prmUnit)

    END SUBROUTINE WriteConvergenceParamFile

!*

    SUBROUTINE ReadEVBParameters(logUnit)

      ! #DES: Read a text file containing diagonal constants and of-diagonal function parameters
      ! #DES: If this cannot be found, run the calculation with all zeroes, dump a .prm file and note in the log file

      USE FileIO, ONLY : OpenFile, CloseFile

      IMPLICIT NONE
      INTEGER, INTENT(IN) :: logUnit
      INTEGER, PARAMETER :: prmUnit = 20
      INTEGER :: state, i, j, dummy
      LOGICAL :: fileExists

      WRITE(logUnit,'(A34)',ADVANCE='NO') "Determining EVB Parameter Data... "

      CALL OpenFile(prmUnit,"evb.prm","read",success=fileExists)

      IF (fileExists) THEN

        DO state = 1, nStates
          READ(prmUnit,*) alpha(state)
        ENDDO

        DO i = 1, nStates
          DO j = i+1, nstates
            READ(prmUnit,*) dummy, dummy, couplingConstant(i,j), couplingExpExpFactor(i,j), couplingGaussExpFactor(i,j)
            couplingConstant(j,i)     = couplingConstant(i,j)
            couplingExpExpFactor(j,i) = couplingExpExpFactor(i,j)
            couplingGaussExpFactor(j,i) = couplingGaussExpFactor(i,j)
          ENDDO
        ENDDO

        CALL CloseFile(prmUnit)
        WRITE(logUnit,'(A22)') "data read from evb.prm"
        WRITE(logUnit,*)

      ELSE

        IF (targetsPresent .EQV. .TRUE.) THEN
          WRITE(logUnit,'(A)') "no EVB parameter file provided: using targets to guess EVB parameters"
          WRITE(logUnit,*)
        ELSE
          CALL WriteEmptyEVBParameterFile()
          WRITE(logUnit,'(A)') "no EVB parameter file or free energy targets provided; template written to evb.prm"
          WRITE(logUnit,*)
        ENDIF

      ENDIF

    END SUBROUTINE ReadEVBParameters

!*

    SUBROUTINE PrintEVBParameters(logUnit)

      IMPLICIT NONE
      INTEGER, INTENT(IN) :: logUnit
      INTEGER :: state, i, j

      WRITE(logUnit,'(A15)') "EVB parameters:"
      WRITE(logUnit,*)

      WRITE(logUnit,'(A22)') "* Diagonal Constants *"
      WRITE(logUnit,'(A13)') "State   Alpha"
      DO state = 1, nStates
        WRITE(logUnit,'(I3,3X,F7.3)') state, alpha(state)
      ENDDO
      WRITE(logUnit,*)

      WRITE(logUnit,'(A27)') "* Off-Diagonal Parameters *"
      WRITE(logUnit,'(A)') "i   j   A   mu   eta"
      DO i = 1, nStates
        DO j = i+1, nStates
          WRITE(logUnit,'(2I4,3F7.3)') i, j, couplingConstant(i,j), couplingExpExpFactor(i,j), couplingGaussExpFactor(i,j)
        ENDDO
      ENDDO
      WRITE(logUnit,*)

      WRITE(logUnit,*) "* Coupling Coordinate *"
      IF (useEnergyGapCoupling) WRITE(logUnit,'(A)') "Using energy gap as coupling coordinate"

    END SUBROUTINE PrintEVBParameters

!*

    SUBROUTINE WriteEmptyEVBParameterFile

      ! #DES: Write an empty EVB parameters file to evb.prm for the number of states in settings.nml

      USE FileIO, ONLY : OpenFile, CloseFile

      IMPLICIT NONE
      INTEGER, PARAMETER :: prmUnit = 20
      INTEGER :: i, j

      CALL OpenFile(prmUnit,"evb.prm","write")

      DO i = 1, nStates
        WRITE(prmUnit,'(F7.3)') 0.0d0
      ENDDO

      DO i = 1, nStates
        DO j = i+1 , nStates
          WRITE(prmUnit,'(2I4,3F7.3)') i, j, 0.0d0, 0.0d0, 0.0d0
        ENDDO
      ENDDO

      CALL CloseFile(prmUnit)

    END SUBROUTINE WriteEmptyEVBParameterFile

!*

    SUBROUTINE ProcessNamelist(logUnit)

      ! #DES: Read the namelist settings.nml, check the inputs, determine dependent variables and print to log file

      IMPLICIT NONE
      INTEGER, INTENT(IN) :: logUnit
      CHARACTER(100) :: command

      CALL ReadNameList()
      CALL CheckNameList()

      !Compute dependent information
      beta = 1.0d0 / (0.0019872041d0 * temperature)
      IF (dGTS /= HUGE(0.0d0) .AND. dGPS /= HUGE(0.0d0)) targetsPresent = .TRUE.
      command = "mkdir "//TRIM(ADJUSTL(outdir))
      !CALL SYSTEM(command) ! Write functional directory creation routine

      CALL PrintNameList(logUnit)

    END SUBROUTINE ProcessNamelist

!*

    SUBROUTINE ReadNamelist

      ! #DES: Open, read and close the namelist file settings.nml

      USE FileIO, ONLY : OpenFile, CloseFile
      IMPLICIT NONE
      INTEGER, PARAMETER :: nmlUnit = 20

      CALL OpenFile(nmlUnit,"settings.nml","read")

      nStates     = 0
      stateA      = 0
      stateB      = 0
      nFepSteps   = 0
      stepTS      = 0
      temperature = -1.0d0

      ! useable defaults
      nSkip      =  0
      rcCoeffA   =  1.0d0
      rcCoeffB   = -1.0d0
      nBins      =  0
      minPop     =  1
      fileBase   =  "EnergyFile"
      outDir     =  "output"
      useEnergyGapCoupling = .TRUE.
      dGTS = HUGE(0.0d0)
      dGPS = HUGE(0.0d0)

      READ(nmlUnit,NML=settings)
      CALL CloseFile(nmlUnit)

  END SUBROUTINE ReadNameList

!*

  SUBROUTINE CheckNameList

    ! #DES: Confirm that all input values are in the correct bounds

    IMPLICIT NONE

    IF (nStates < 2)                            STOP ": ReadNameList: ILLEGAL NUMBER OF EVB STATES REQUESTED"
    IF (stateA < 1 .OR. stateA > nStates)       STOP ": ReadNameList: Illegal value of stateA"
    IF (stateB < 1 .OR. stateB > nStates)       STOP ": ReadNameList: Illegal value of stateB"
    IF (nFepSteps  < 2)                         STOP ": ReadNameList: ILLEGAL NUMBER OF FEP STEPS REQUESTED"
    IF (stepTS < 1 .OR. stepTS > nFepSteps)     STOP ": ReadNameList: Illegal value of StepTS"
    IF (nBins < 1)                              STOP ": ReadNameList: Illegal value of nBins"
    IF (fileBase == "")                         STOP ": ReadNameList: Illegal value of fileBase"
    IF (temperature < 0.0d0)                    STOP ": ReadNameList: Illegal value of temperature"

  END SUBROUTINE CheckNameList

!*

  SUBROUTINE PrintNameList(logUnit)

    ! #DES: Reprint all namelist input to the log file

    IMPLICIT NONE
    INTEGER, INTENT(IN) :: logUnit

    WRITE(logUnit,'(A)') "FEP/US parameters read from settings.nml:"
    WRITE(logUnit,*)
    WRITE(logUnit,'(A)') "* EVB Settings *"
    WRITE(logUnit,'(A,I6)')   "Number of EVB States         : ", nStates
    WRITE(logUnit,'(A,I6)')   "EVB State a                  : ", stateA
    WRITE(logUnit,'(A,I6)')   "EVB State b                  : ", stateB
    WRITE(logUnit,*)
    WRITE(logUnit,'(A)') "* MD/FEP Settings *"
    WRITE(logUnit,'(A,I6)')   "Number of Mapping Simulations: ", nFepSteps
    WRITE(logUnit,'(A,I6)')   "Simulation Steps to Skip     : ", nSkip
    WRITE(logUnit,'(A,F6.2)') "Simulation Temperature       : ", temperature
    WRITE(logUnit,'(A,A)')    "Energy File Prefix           : ", TRIM(ADJUSTL(fileBase))
    WRITE(logUnit,*)
    WRITE(logUnit,'(A)') "* FEP/US Settings *"
    WRITE(logUnit,'(A,I6)')   "Step for LRA calculations    : ", stepTS
    WRITE(logUnit,'(A,I6)')   "Number of FEP/US Bins        : ", nBins
    WRITE(logUnit,'(A,I6)')   "Population Threshold of Bins : ", minPop
    WRITE(logUnit,'(A,L6)')   "Use Energy Gap for Coupling  : ", useEnergyGapCoupling
    WRITE(logUnit,'(A,A)')    "Output Directory             : ", outDir

    IF (targetsPresent .EQV. .TRUE.) THEN
      WRITE(logUnit,'(A,F6.2)') "Target Free Energy RS -> TS  : ", dGTS
      WRITE(logUnit,'(A,F6.2)') "Target Free Energy RS -> PS  : ", dGPS
      WRITE(logUnit,*)
    ENDIF

    WRITE(logUnit,*)

  END SUBROUTINE PrintNameList

!*

  SUBROUTINE CreateFileNames(ext,prefix,fileNames)

    IMPLICIT NONE
    CHARACTER(*), INTENT(IN) :: ext, prefix
    CHARACTER(500), INTENT(OUT) :: fileNames(nFepSteps)
    CHARACTER(3) :: stepString
    INTEGER :: step

    DO step = 1, nFepSteps
      WRITE(stepString,'(I0.3)') step - 1
      fileNames(step) = TRIM(ADJUSTL(prefix))//stepString//'.'//TRIM(ADJUSTL(ext))
    ENDDO

  END SUBROUTINE CreateFileNames

!*

  SUBROUTINE ReadCoordinateFiles(logUnit)

    ! #DES: Read the trajectory from the .dcd files corresponding to each .en file

    USE DCDFiles, ONLY : ReadDCD, CountAtoms

    IMPLICIT NONE
    INTEGER, INTENT(IN) :: logUnit
    CHARACTER(500) :: fileNames(nFepSteps)
    INTEGER :: step, nAtoms, maxAtoms

    WRITE(logUnit,'(A)') "Reading Simulation Trajectory Data"
    WRITE(logUnit,*)
    CALL CreateFileNames('dcd',trajecBase,fileNames)

    maxAtoms = 0
    DO step = 1, NFepSteps
      CALL CountAtoms(fileNames(step),nAtoms)
      IF (nAtoms > maxAtoms) maxAtoms = nAtoms
    ENDDO

    ALLOCATE(trajectory(nFepSteps,3,maxAtoms,maxTimesteps))

    DO step = 1, nFepSteps
      WRITE(logUnit,*) "Reading trajectory from ", TRIM(ADJUSTL(fileNames(step)))
      CALL ReadDCD(fileNames(step),trajectory(step,:,:,:))
    ENDDO

    WRITE(logUnit,*)
    WRITE(logUnit,'(A)') "Finished Reading Simulation Trajectory Data"
    WRITE(logUnit,*)

  END SUBROUTINE ReadCoordinateFiles

!*

  SUBROUTINE ReadFormattedEnergyFiles(logUnit)

    ! #DES: Read the nFepsSteps energy files into the arrays stateEnergy and coeffs
    ! Need to edit this such that it does not need to know beforehand how many timesteps are in each fep file
    ! and if a fep file is missing/incomplete, instead of dying just have however many points were read OK in stateEnergy(fepstep,:)
    ! This saves knowing the length of simulations, allows analysis on partially completed simulations.

    USE FileIO, ONLY : OpenFile, CloseFile, fileLengths

    IMPLICIT NONE
    INTEGER, INTENT(IN) :: logUnit
    INTEGER, PARAMETER :: eneUnit = 10
    CHARACTER(500) :: fileNames(nFepSteps)
    INTEGER :: step, state, timestep, ios, type
    LOGICAL :: openSuccess

    WRITE(logUnit,*) ""
    WRITE(logUnit,'(A30)') "Reading Simulation Energy Data"
    WRITE(logUnit,*) ""

    CALL CreateFileNames('en',fileBase,fileNames)

    ReadSteps: DO step = 1, nFepSteps

      CALL OpenFile(eneUnit,fileNames(step),"read",success=openSuccess)
      IF (nTimesteps(step) == 0) CYCLE ReadSteps

      ! This would be the place to allocate space for this fep simulation's energy data
      DO timestep = 1, nTimesteps(step)
        DO state = 1, nStates

          READ(eneUnit,'(F15.8)',ADVANCE='NO',IOSTAT=ios) coeffs(timestep,step,state)
          IF (ios /= 0) THEN
            WRITE(logUnit,'(A,I0.2,A,I0.6,A,I0.4)') "Failed to read expected coefficient for state ", state, " at timestep ", timestep, " of FEP simulation ", step
            WRITE(logUNit,'(A)') "Moving to next FEP simulation"
            CYCLE ReadSteps
          ENDIF

          DO type = 1, nEnergyTypes
            READ(eneUnit,'(F15.8)',ADVANCE='NO',IOSTAT=ios) stateEnergy(timestep,step,state,type)
          ENDDO; READ(eneUnit,*)

        ENDDO
      ENDDO
      CALL CloseFile(eneUnit)

      WRITE(logUnit,'(A4,1X,I0.6,1X,A14,1X,A)') "Read", timestep-1, "timesteps from", TRIM(ADJUSTL(fileNames(step)))

    ENDDO ReadSteps

    WRITE(logUnit,*)
    WRITE(logUnit,'(A39)') "Finished Reading Simulation Energy Data"
    WRITE(logUnit,*)

  END SUBROUTINE ReadFormattedEnergyFiles

!*

  SUBROUTINE ReadEnergyTypes(logUnit)

    USE FileIO, ONLY : OpenFile, CloseFile, FileLength
    IMPLICIT NONE
    INTEGER, INTENT(IN) :: logUnit
    INTEGER, PARAMETER :: eneUnit = 22
    LOGICAL :: fileExists = .FALSE.
    INTEGER :: line

    WRITE(logUnit,'(A)') "* Determine Energy Types *"; WRITE(logUnit,*)
    CALL OpenFile(eneUnit,"EnergyNames.dat","read",success=fileExists)

    IF (fileExists .EQV. .TRUE.) THEN

      WRITE(logUnit,'(A)') "Reading energy names from EnergyNames.dat"
      ! This is a fudge to get around the way FileLength is written - add a 2nd version that just takes a unit #, counts, then rewinds
      CALL CloseFile(eneUnit)
      nEnergyTypes = FileLength("EnergyNames.dat")
      CALL OpenFile(eneUnit,"EnergyNames.dat","read",success=fileExists)

      ALLOCATE(energyNames(nEnergyTypes))
      DO line = 1, nEnergyTypes
        READ(eneUnit,'(A)') EnergyNames(line)
        IF (TRIM(ADJUSTL(energyNames(line))) == "") STOP "Error: Remove blank spaces in EnergyNames.dat"
      ENDDO

    ELSE

      WRITE(logUnit,'(A)') "No file EnergyNames.dat - using Q Defaults"
      ! No file so assume Q types
      CALL DefineQTypes()

    ENDIF

    WRITE(logUnit,'(A,1X,I4)') "Number of Types: ", nEnergyTypes
    DO line = 1, nEnergyTypes
      WRITE(logUnit,'(A)') trim(adjustl(energyNames(line)))
    ENDDO
    WRITE(logUnit,*)

  END SUBROUTINE ReadEnergyTypes
!*

  SUBROUTINE DefineQTypes

    ! #DES: Defines the energy types as found in Q output files for convenience

    IMPLICIT NONE

    nEnergyTypes = 14
    ALLOCATE(energyNames(nEnergyTypes))

    energyNames(1)  = "    Q_TOTAL"; energyNames(2)  = "     Q_BOND";
    energyNames(3)  = "    Q_ANGLE"; energyNames(4)  = "  Q_TORSION";
    energyNames(5)  = " Q_IMPROPER"; energyNames(6)  = "     Q_ELEC";
    energyNames(7)  = "      Q_VDW"; energyNames(8)  = "   Q-Q_ELEC";
    energyNames(9)  = "    Q-Q_VDW"; energyNames(10) = "   Q-P_ELEC";
    energyNames(11) = "    Q-P_VDW"; energyNames(12) = "   Q-W_ELEC";
    energyNames(13) = "    Q-W_VDW"; energyNames(14) = "    Q_RESTR";

  END SUBROUTINE DefineQTypes

!*

  SUBROUTINE WriteEnergyTypes(fileUnit)

    ! #DES: Write the names of the energy types to fileUnit in .csv format

    IMPLICIT NONE
    INTEGER, INTENT(IN) :: fileUnit
    INTEGER :: type

    DO type = 1, SIZE(energyNames)
      WRITE(fileUnit,'(A11,A1)',ADVANCE='no') adjustr(energyNames(type)), ","
    ENDDO

  END SUBROUTINE WriteEnergyTypes

END MODULE Input

max_stars_repo_name: MJLMills/fepcat

max_issues_repo_name: MJLMills/fepcat

max_forks_repo_name: MJLMills/fepcat


text:
C++:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
C       PROGRAM MASS_HUBBLE
C
C  This program, from HUBBLE computes the mass of pseudo-Hubble models
C
C JLP
C Version of 19-11-86
C--:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
        PROGRAM MASS_HUBBLE
	PARAMETER (IDIM=10000)
	IMPLICIT DOUBLE PRECISION (A-H,O-Z)
	REAL*8 RHOS(IDIM),RAD(IDIM),XMAG(IDIM)
	REAL*4 ALPHA
	CHARACTER ANS*1,NAME*40
	COMMON/MODEL/ALPHA,RR2

	PI=3.14159265358979323846
10	FORMAT(A)
	OPEN(1,FILE='mass_hubble.log',STATUS='unknown')
	WRITE(1,20)
	PRINT 20
20	FORMAT(' PROGRAMME MASS_HUBBLE ',
     1	'---VERSION OF NOV 19TH 1986---')

C***********************************************************************
40	PRINT 50
50	FORMAT(' MENU :',/,
     1	' 1. COMPUTING THE MASS OF A PSEUDO-HUBBLE MODEL',/,
     1	' WITH A SPACE DENSITY IN 1/(1+S)**ALPHA',/,
     1	' (I.E. PROJECTED DENSITY IN 1/(1+S)**(ALPHA-1) )',/,
     1	' 2. INVERTING A HUBBLE PROFILE',/,
     1	' 3. COMPUTING THE MASS FROM THE SPACE DENSITY',
     1	' OBTAINED IN 2',/,
     1	' 4. COMPUTING THE TOTAL LUMINOSITY OF A HUBBLE'
     1	' MODEL',/,
     1	' 5. GENERATION OF A CURVE WITH THE LUMINOSITY',/,
     1	' 6. COMPUTING THE MASS OF A PSEUDO-HUBBLE MODEL',/,
     1	' WITH A POTENTIAL IN 1/(1+S**2)**ALPHA',/,
     1	' 10. EXIT',/,
     1	' ENTER OPTION : ',$)
	READ(5,*) IOPT

	IF(IOPT.NE.10.AND.IOPT.NE.3)THEN
	PRINT 55
55	FORMAT(' ALPHA : ',$)
	READ(5,*) ALPHA
	ENDIF

C---------------------------------------------------------------------
C Options 1 and 3 : Computing the mass
C--------------------------------------------------------------------
	IF(IOPT.EQ.1.OR.IOPT.EQ.3.OR.IOPT.EQ.4.OR.IOPT.EQ.6)THEN
100	PRINT 101
101	FORMAT(' LIMITS FOR THE INTEGRATION',
     1	'(to exit type 0.,0.) : ',$)
	READ(5,*) XLOW,XHIGH

	IF(XHIGH.NE.0.)THEN
	XMASS=XMASS_HUBBLE(XLOW,XHIGH,IOPT)
	WRITE(1,102)ALPHA,XLOW,XHIGH,XMASS
	PRINT 102,ALPHA,XLOW,XHIGH,XMASS
102	FORMAT(' ALPHA :',F11.4,/,' LOWER LIMIT :',E11.4,
     1	' UPPER LIMIT :',E11.4,' INTEGRAL :',E12.4,/)
	GO TO 100
	ENDIF

	GO TO 40
	ENDIF

C---------------------------------------------------------------------
C Option 2 : Inversion of a Hubble profile
C--------------------------------------------------------------------
	IF(IOPT.EQ.2)THEN
	PRINT *,' NUMBER OF POINTS FOR THE MODEL (100 ?) :'
	READ(5,*) NPOINT
	WRITE(1,201)NPOINT
201	FORMAT(' INVERSION OF A HUBBLE PROFILE WITH :',I5,
     1	' POINTS')

C Setting up the radius array :
	PRINT *,' LOWER AND UPPER LIMITS IN LOG10 FOR THE RADIUS :'
	READ(5,*) XXLOW,XXHIGH
	XSTEP=(XXHIGH-XXLOW)/FLOAT(NPOINT)
	DO I=1,NPOINT
	RAD(I)=10**(XXLOW+(FLOAT(I)*XSTEP))
	END DO

	PRINT 202,RAD(1),RAD(NPOINT)
	WRITE(1,202)RAD(1),RAD(NPOINT)
202	FORMAT(' RADMIN, RADMAX :',E12.4,2X,E12.4)
	CALL INVERT_HUBBLE2(RAD,RHOS,NPOINT)

	GO TO 40
	ENDIF

C---------------------------------------------------------------------
C Option 5 : luminosity
C--------------------------------------------------------------------
	IF(IOPT.EQ.5)THEN
500	PRINT 501
501	FORMAT(' OUTPUT FILE :')
	READ(5,10) NAME
	OPEN(1,FILE=NAME,STATUS='NEW',ERR=500)
	PRINT *,' CORE RADIUS (ARCSEC.)',
     1	' AND CENTRAL SURFACE BRIGHT. (MAG/ARC.)?'
	READ(5,*) RCORE,B0

	PRINT *,' MAX RADIUS AND STEP (ARCSECONDS)'
	READ(5,*) RADMAX,STEP
	NPTS=JIDINT(RADMAX/STEP)
	WRITE(1,*)NPTS

C To call XMASS_HUBBLE I set IOPT to 4 (same option : luminosity)
	IOPT=4
	XLOW=0.
	CTE=2.D0*PI*RCORE*RCORE

	DO I=1,NPTS
	RAD(I)=FLOAT(I)*STEP
	XHIGH=RAD(I)/RCORE
	XWORK=XMASS_HUBBLE(XLOW,XHIGH,IOPT)
	XLUMINO=CTE*XWORK
	XMAG(I)=-2.5*DLOG10(XLUMINO)+B0
	WRITE(1,*)RAD(I),XMAG(I)
	END DO

	CLOSE(1)
	GO TO 40
	ENDIF

C---------------
	PRINT *,' Output in "mass_hubble.log".'
	CLOSE(1)
	STOP
	END
C---------------------------------------------------------------------------
C Function XMASS_HUBBLE
C--------------------------------------------------------------------------
	DOUBLE PRECISION FUNCTION XMASS_HUBBLE(XLOW,XHIGH,IOPT)
	IMPLICIT REAL*8(A-H,O-Z)
	REAL*8 W(6000)
	INTEGER*4 IW(1000)
	EXTERNAL FUN_HUBBLE1
	EXTERNAL FUN_HUBBLE3
	EXTERNAL FUN_HUBBLE4
	EXTERNAL FUN_HUBBLE6

C Absolute and relative accuracy :
        EPSABS=0.0001
        EPSREL=0.001

C Size of the work areas W(real) and IW(integer)
        LW=6000
        LIW=1000

C Calling a rather sophisticated integration routine D01AKF :
        IFAIL=1
	IF(IOPT.EQ.1)THEN
	CALL D01AKF(FUN_HUBBLE1,XLOW,XHIGH,EPSABS,EPSREL,AL,
     1	ABSERR,W,LW,IW,LIW,IFAIL)
	ENDIF
	IF(IOPT.EQ.3)THEN
	CALL D01AKF(FUN_HUBBLE3,XLOW,XHIGH,EPSABS,EPSREL,AL,
     1	ABSERR,W,LW,IW,LIW,IFAIL)
	ENDIF
	IF(IOPT.EQ.4)THEN
	CALL D01AKF(FUN_HUBBLE4,XLOW,XHIGH,EPSABS,EPSREL,AL,
     1	ABSERR,W,LW,IW,LIW,IFAIL)
	ENDIF
	IF(IOPT.EQ.6)THEN
	CALL D01AKF(FUN_HUBBLE6,XLOW,XHIGH,EPSABS,EPSREL,AL,
     1	ABSERR,W,LW,IW,LIW,IFAIL)
	ENDIF

	IF(IFAIL.NE.0) THEN
	PRINT *,' WARNING : FAILURE IN D01AKF,  IFAIL =',IFAIL
	END IF

	XMASS_HUBBLE=AL

	RETURN
	END
C---------------------------------------------------------------------------
C Function FUN_HUBBLE1
C Called by D01AKF for option 1 (mass)
C Generates the value to be integrated at the point X
C---------------------------------------------------------------------------
	DOUBLE PRECISION FUNCTION FUN_HUBBLE1(X)
	IMPLICIT REAL*8(A-H,O-Z)
	REAL*4 ALPHA
	COMMON/MODEL/ALPHA,RR2
	FUN_HUBBLE1=(X*X)/((1.D0+X)**ALPHA)
	RETURN
	END
C---------------------------------------------------------------------------
C Function FUN_HUBBLE3
C Called by D01AKF for option 3 (mass from spatial density)
C Generates the value to be integrated at the point X
C---------------------------------------------------------------------------
	DOUBLE PRECISION FUNCTION FUN_HUBBLE3(X)
	PARAMETER (IDIM=10000)
	IMPLICIT DOUBLE PRECISION (A-H,O-Z)
	REAL*8 SPLINE2(IDIM),K2(IDIM)
	COMMON/SPLINE2/NCAP72,SPLINE2,K2
	IFAIL=1
	CALL E02BBF(NCAP72,K2,SPLINE2,X,RESULT,IFAIL)
	  IF(IFAIL.NE.0) THEN
	  PRINT *,' WARNING : FAILURE IN E02BBF,  IFAIL =',IFAIL
	  RESULT=0.D0
	  END IF

	FUN_HUBBLE3=X*X*RESULT
	RETURN
	END
C---------------------------------------------------------------------------
C Function FUN_HUBBLE4
C Called by D01AKF for option 4 (Luminosity)
C Generates the value to be integrated at the point X
C---------------------------------------------------------------------------
	DOUBLE PRECISION FUNCTION FUN_HUBBLE4(X)
	IMPLICIT REAL*8(A-H,O-Z)
	REAL*4 ALPHA
	COMMON/MODEL/ALPHA,RR2
	FUN_HUBBLE4=X/((1.D0+X)**ALPHA)
	RETURN
	END
C---------------------------------------------------------------------------
C Function FUN_HUBBLE6
C Called by D01AKF for option 6 (mass with a potential in (1+R**2)**(-ALPHA) )
C Generates the value to be integrated at the point X
C---------------------------------------------------------------------------
	DOUBLE PRECISION FUNCTION FUN_HUBBLE6(X)
	IMPLICIT REAL*8(A-H,O-Z)
	REAL*4 ALPHA
	COMMON/MODEL/ALPHA,RR2
	RHO=((1+X*X)**(-2.-ALPHA))*(X*X*(3.D0-2.D0*(1.D0+ALPHA))+3.D0)
	FUN_HUBBLE6=X*X*RHO
	RETURN
	END
C---------------------------------------------------------------------------
C Subroutine INVERT_HUBBLE2 to compute spatial density RHOS
C Calls NAG routine D01AMF
C designed to integrate over a semi-infinite or infinite interval
C--------------------------------------------------------------------------
	SUBROUTINE INVERT_HUBBLE2(RAD,RHOS,NPOINT)
	PARAMETER (IDIM=10000)
	IMPLICIT REAL*8(A-H,O-Z)
	REAL*8 W(4000),SPLINE2(IDIM),K2(IDIM)
	REAL*8 RHOS(IDIM),RAD(IDIM)
	REAL*4 ALPHA
	INTEGER*4 IW(502)
	COMMON/MODEL/ALPHA,RR2
	COMMON/SPLINE2/NCAP72,SPLINE2,K2
	EXTERNAL FUN_HUBBLE2

C Loop on all the values contained in the array RAD(NPOINT)
	DO I=1,NPOINT

C Absolute and relative accuracy :
        EPSABS=1.E-06
        EPSREL=1.E-04
C Size of the work areas W(real) and IW(integer)
        LW=4000
        LIW=502
C Lower limit for the integration
        XLOW=0.
        IFAIL=1
C INF=1 (BOUND--> +INF) INF=-1 (-INF --> BOUND)  INF=2 (-INF --> +INF)
	INF=1
	BOUND=RAD(I)
	RR2=RAD(I)*RAD(I)

C Calling a rather sophisticated integration routine D01AMF :
        CALL D01AMF(FUN_HUBBLE2,BOUND,INF,EPSABS,EPSREL,
     1	RESULT,ABSERR,W,LW,IW,LIW,IFAIL)

	IF(IFAIL.NE.0) THEN
	PRINT *,' WARNING : FAILURE IN D01AMF,  IFAIL =',IFAIL
	ENDIF

	RHOS(I)=RESULT
	WRITE(1,200)I,RHOS(I)
200	FORMAT(' I=',I4,'RHO=',E12.4)
	END DO

C Interpolates with cubic-spline :
	CALL SPLINE_INTER(RAD,RHOS,NPOINT,SPLINE2,NCAP72,K2)

	RETURN
	END
C---------------------------------------------------------------------------
C Function FUN_HUBBLE2
C Called by D01AMF
C Generates the value to be integrated at the point X
C---------------------------------------------------------------------------
	DOUBLE PRECISION FUNCTION FUN_HUBBLE2(X)
	IMPLICIT REAL*8(A-H,O-Z)
	REAL*4 ALPHA
	COMMON/MODEL/ALPHA,RR2
	FUN_HUBBLE2=(1.D0+X*(2.D0+ALPHA))*DSQRT(X*X-RR2)
     1	/(X*X*((1.D0+X)**(2.+ALPHA)))
	RETURN
	END
C---------------------------------------------------------------------------
	include 'jlpsub:project.for'

ext: for


module Performance

use constants

implicit none

type, public :: RunPerfData

        integer(kind=long) :: nProcs
        integer(kind=long) :: nSpatialPoints
        integer(kind=long) :: nRowLocal
        integer(kind=long) :: nColLocal
        integer(kind=long) :: nRowGlobal
        integer(kind=long) :: nColGlobal
        real :: MatSizeLocal_GB
        real :: MatSizeGlobal_GB
        real :: TimeWorkList
        real :: TimeFill
        real :: TimeSolve
        real :: TimeCurrent
        real :: TimeTotal
        real :: GflopsFillLocal
        real :: GflopsFillGlobal
        real :: GflopsSolveLocal
        real :: GflopsSolveGlobal
        real :: GflopsCurrentLocal
        real :: GflopsCurrentGlobal

end type RunPerfData

type, public :: MemoryUsage

        integer :: nAllocations
        real :: MBytes

end type MemoryUsage

type(MemoryUsage) :: Mem
type(RunPerfData) :: Perf

end module Performance


text:
      subroutine fd_insolation(lyear,idayjl,isecdy,xlat, rsol)
      implicit none

c Now, during the time period of 1978 through 1998,
c the mean value of daily averages for the solar constant from six different satellites
c yielded a solar constant of 1366.1 Wm2. This same source offers a minimum - maximum range of
c the readings for 1363 - 1368 Wm2. Adjustments yielded "a solar constant" calculated to 1366.22 Wm2.
c [Royal Meteorological Institute of Belgium: Department of Aerology;
c http://remotesensing.oma.be/RadiometryPapers/article2.html.].
c	c.gentemann gentemann@remss.com
c	6/2009

	real(4), parameter :: solar_constant=1366.2   !watts/m**2
      real(4), parameter :: pi=3.141592654

c     inputs and outputs
	integer(4) lyear,idayjl,isecdy  ! year, julian day, seconds in day
	real(4)    xlat                 ! latitude of observation
	real(4)    rsol                 ! average insolation watts/m**2
	real(4) sunlat,sunlon,sundis,rr,H,sinh,Q

	call sunloc1(lyear,idayjl,isecdy, sunlat,sunlon,sundis)	  !sundis is radius/(mean radius)

      rr=-tand(xlat)*tand(sunlat);
      if(rr> 1) rr= 1.
      if(rr<-1) rr=-1.

      H = acos(rr);
      sinH=sqrt(1 - rr**2) !sqrt(1-cos(h)**2) =sin(h)
	Q=sind(xlat)*sind(sunlat)*H + cosd(xlat)*cosd(sunlat)*sinH
      rsol = (solar_constant/pi)*Q/sundis**2

      return
      end


text:
C$Procedure      CMMORE ( Command Loop---More Commands)

      LOGICAL  FUNCTION CMMORE ( COMMND )

C$ Abstract
C
C    Determine whether or not more command loop processing
C    should be performed.
C
C$ Disclaimer
C
C     THIS SOFTWARE AND ANY RELATED MATERIALS WERE CREATED BY THE
C     CALIFORNIA INSTITUTE OF TECHNOLOGY (CALTECH) UNDER A U.S.
C     GOVERNMENT CONTRACT WITH THE NATIONAL AERONAUTICS AND SPACE
C     ADMINISTRATION (NASA). THE SOFTWARE IS TECHNOLOGY AND SOFTWARE
C     PUBLICLY AVAILABLE UNDER U.S. EXPORT LAWS AND IS PROVIDED "AS-IS"
C     TO THE RECIPIENT WITHOUT WARRANTY OF ANY KIND, INCLUDING ANY
C     WARRANTIES OF PERFORMANCE OR MERCHANTABILITY OR FITNESS FOR A
C     PARTICULAR USE OR PURPOSE (AS SET FORTH IN UNITED STATES UCC
C     SECTIONS 2312-2313) OR FOR ANY PURPOSE WHATSOEVER, FOR THE
C     SOFTWARE AND RELATED MATERIALS, HOWEVER USED.
C
C     IN NO EVENT SHALL CALTECH, ITS JET PROPULSION LABORATORY, OR NASA
C     BE LIABLE FOR ANY DAMAGES AND/OR COSTS, INCLUDING, BUT NOT
C     LIMITED TO, INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND,
C     INCLUDING ECONOMIC DAMAGE OR INJURY TO PROPERTY AND LOST PROFITS,
C     REGARDLESS OF WHETHER CALTECH, JPL, OR NASA BE ADVISED, HAVE
C     REASON TO KNOW, OR, IN FACT, SHALL KNOW OF THE POSSIBILITY.
C
C     RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF
C     THE SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY
C     CALTECH AND NASA FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE
C     ACTIONS OF RECIPIENT IN THE USE OF THE SOFTWARE.
C
C$ Required_Reading
C
C     None.
C
C$ Keywords
C
C     Command Loop
C
C$ Declarations

      IMPLICIT NONE
      CHARACTER*(*)         COMMND

C$ Brief_I/O
C
C     VARIABLE  I/O  DESCRIPTION
C     --------  ---  --------------------------------------------------
C     COMMND     I   A command to be processed by CMLOOP
C
C     The function returns .TRUE. if the command is not the "exit"
C     command. If it is the exit command it returns .FALSE.
C
C$ Detailed_Input
C
C     COMMND     A commmand that should be acted on by CMLOOP
C
C
C$ Detailed_Output
C
C     The function returns .TRUE. if this is not the exit command.
C     The meaning being "there is still more to do in CMLOOP."
C
C     If the input command is equivalent to the exit command
C     (Same words when converted to uppercase) The function
C     returns .FALSE.  The intended meaning is "there is nothing
C     left for CMLOOP to do but cleanup and return."
C
C$ Parameters
C
C     None.
C
C$ Files
C
C     None.
C
C$ Exceptions
C
C     Error free.
C
C$ Particulars
C
C     This is utility function for use by CMLOOP.  It is the
C     function tested each pass through the loop to see if the
C     loop has finished its work
C
C$ Examples
C
C     See CMLOOP.  There is no other use for this function.
C
C$ Restrictions
C
C     None.
C
C$ Author_and_Institution
C
C       W.L. Taber      (JPL)
C
C$ Literature_References
C
C       None.
C
C$ Version
C
C-    Command Loop Version 1.0.0, 4-AUG-1995 (WLT)
C
C
C-&

C$ Index_Entries
C
C     More command processing required
C
C-&
C     SPICELIB Functions
C
      INTEGER               RTRIM
      LOGICAL               NECHR

C
C     Local Variables.
C
      INTEGER               WDSIZE
      PARAMETER           ( WDSIZE = 32 )

      CHARACTER*(WDSIZE)    EXIT

      INTEGER               LC
      INTEGER               I
      INTEGER               R

      LOGICAL               FIRST
      SAVE

      DATA                  FIRST / .TRUE. /

C
C     On the first pass we fetch the "exit" command and
C     spruce it up a bit for use when comparing with
C     the input command.
C
      IF ( FIRST ) THEN

         FIRST = .FALSE.

         CALL TRNLAT ( 'EXIT',       EXIT )
         CALL CMPRSS ( ' ', 1, EXIT, EXIT )
         CALL LJUST  (  EXIT,        EXIT )
         R  = RTRIM  (  EXIT )

      END IF
C
C     If the input command is shorter than the non-blank
C     length of EXIT, then this cannot be the exit command.
C     There is more to do.
C
C     Note we assign a value to CMMORE so that the compiler
C     won't have a fit about having a function unassigned.
C     The if conditions below ensure that we assign a value
C     but most compilers aren't smart enough to figure that
C     out.
C
      CMMORE = .TRUE.
      LC     = LEN(COMMND)

      IF ( LC .LT. R ) THEN
         CMMORE = .TRUE.
         RETURN
      END IF

C
C     Check to see if the input command matches the exit command.
C     We do this a character at a time.  We search from the
C     left to right, because most commands are not EXIT and this
C     allows us to quit early in the process.
C
      DO I = 1, R

         IF ( NECHR( COMMND(I:I), EXIT(I:I) ) ) THEN
            CMMORE = .TRUE.
            RETURN
         END IF

      END DO

C
C     It's looking like this might be it.  See if the rest of
C     the input command is blank.
C
      IF ( LC.EQ. R ) THEN
C
C        We've got an exact match.  There are no more commands
C        to look at.
C
         CMMORE = .FALSE.

      ELSE IF ( LC .GT. R ) THEN
C
C        There will be more commands only if the rest of the input
C        command is non-blank.
C
         CMMORE =  COMMND(R+1:) .NE. ' '

      END IF

      RETURN
      END

max_stars_repo_name: agforero/FTFramework

max_issues_repo_name: agforero/fortran-testing-framework

max_forks_repo_name: agforero/fortran-testing-framework


PROGRAM bad_nested_loops_2
INTEGER :: i, j, product
outer: DO i = 1, 3
   inner: DO j = 1, 3
      product = i * j
      WRITE (*,*) i, ' * ', j, ' = ', product
END DO outer
END PROGRAM bad_nested_loops_2


text:
! Routine for linear approximation
! Igor Lopes, February 2015
SUBROUTINE POLYLIN(NPOINT,X,F)
    IMPLICIT NONE
    ! Arguments
    INTEGER NPOINT
    REAL(8) :: XMIN, FMIN
    REAL(8),DIMENSION(4) :: X , F
    ! Locals
    REAL(8) :: A0,A1
    ! Begin algorithm
    IF(NPOINT.EQ.1)THEN
        A1=F(2)
        A0=F(1)-F(2)*X(1)
    ELSEIF(NPOINT.EQ.2)THEN
        A1=(F(2)-F(1))/(X(2)-X(1))
        A0=F(1)-A1*X(1)
    ELSE
        WRITE(*,*)'Wrong number of input points in POLYLIN'
        GOTO 99
    ENDIF
    WRITE(*,*)'Approximation coefficients:'
    WRITE(*,*)'a0=',A0
    WRITE(*,*)'a1=',A1
    WRITE(11,*)'Approximation coefficients:'
    WRITE(11,*)'a0=',A0
    WRITE(11,*)'a1=',A1
99 CONTINUE
END SUBROUTINE

max_stars_repo_name: iarlopes/GPROPT

max_issues_repo_name: iarlopes/GPROPT

max_forks_repo_name: iarlopes/GPROPT


text:
C
C
      SUBROUTINE EOS(TEMP,TEMP0,P,RHO,E,M)
C
C	PROGRAM TO USE ANEOS TO MATCH INPUT ENERGY DENSITY
C
C       This is the simplest possible algorithm, that sometimes
C       fails due to classic Newton-Raphson nonconverging loop,
C       especially near phase transitions.  If trouble is
C       encountered, use EOS_safe, a routine that backs the N-R
C       method up with bisection.
C
C
C        INPUTS:  TEMP0 = Initial estimate of temperature
C                 RHO   = Density
C                 E     = Internal energy to be matched by EOS
C                 M     = Material number
C
C       OUTPUTS:  TEMP  = Temperature (in eV) at given RHO, E
C                 P     = Pressure, P(RHO,E)
C
C       All units are cgs-eV
C
      IMPLICIT REAL*8 (A-H,O-Z)
      TEMP=TEMP0
      IT=0
 10   CONTINUE
      IT=IT+1
      IF(IT.GT.50) STOP 'TOO MANY ITERATIONS--EOS'
      CALL ANEOSD (TEMP,RHO,P,ETRIAL,S,CV,DPDT,DPDR,FKROS,CS,KPA,M)
      IF (ABS(E-ETRIAL) .LT. 1.0D-9*E) RETURN
      TEMP=TEMP+(E-ETRIAL)/CV
      GOTO 10
      END

max_stars_repo_name: chreinhardt/M-ANEOS

max_issues_repo_name: chreinhardt/M-ANEOS

max_forks_repo_name: chreinhardt/M-ANEOS


text:
module reader
!! simple file reading procedures

use, intrinsic :: iso_fortran_env, only: stderr=>error_unit
use phys_consts, only: wp, debug
use pathlib, only : get_suffix

implicit none (type, external)
private
public :: get_simsize3, get_simsize2, get_grid2, get_Efield, get_precip, get_neutral2, get_neutral3

interface ! reader_{raw,hdf5,nc4}.f90
module subroutine get_simsize2_raw(path, llon, llat)
character(*), intent(in) :: path
integer, intent(out) :: llon, llat
end subroutine get_simsize2_raw

module subroutine get_simsize2_hdf5(path, llon, llat)
character(*), intent(in) :: path
integer, intent(out) :: llon, llat
end subroutine get_simsize2_hdf5

module subroutine get_simsize2_nc4(path, llon, llat)
character(*), intent(in) :: path
integer, intent(out) :: llon, llat
end subroutine get_simsize2_nc4

module subroutine get_simsize3_raw(path, lx1, lx2all, lx3all)
character(*), intent(in) :: path
integer, intent(out) :: lx1, lx2all
integer, intent(out), optional :: lx3all
end subroutine get_simsize3_raw

module subroutine get_simsize3_hdf5(path, lx1, lx2all, lx3all)
character(*), intent(in) :: path
integer, intent(out) :: lx1, lx2all
integer, intent(out), optional :: lx3all
end subroutine get_simsize3_hdf5

module subroutine get_simsize3_nc4(path, lx1, lx2all, lx3all)
character(*), intent(in) :: path
integer, intent(out) :: lx1, lx2all
integer, intent(out), optional :: lx3all
end subroutine get_simsize3_nc4

module subroutine get_grid2_raw(path, mlonp, mlatp)
character(*), intent(in) :: path
real(wp), dimension(:), intent(out) :: mlonp, mlatp
end subroutine get_grid2_raw

module subroutine get_grid2_hdf5(path, mlonp, mlatp)
character(*), intent(in) :: path
real(wp), dimension(:), intent(out) :: mlonp, mlatp
end subroutine get_grid2_hdf5

module subroutine get_grid2_nc4(path, mlonp, mlatp)
character(*), intent(in) :: path
real(wp), dimension(:), intent(out) :: mlonp, mlatp
end subroutine get_grid2_nc4

module subroutine get_Efield_raw(path, flagdirich,E0xp,E0yp,Vminx1p,Vmaxx1p,Vminx2pslice,Vmaxx2pslice,Vminx3pslice,Vmaxx3pslice)
character(*), intent(in) :: path
integer, intent(out) :: flagdirich
real(wp), dimension(:,:), intent(out) :: E0xp,E0yp,Vminx1p,Vmaxx1p
real(wp), dimension(:), intent(out) :: Vminx2pslice,Vmaxx2pslice,Vminx3pslice,Vmaxx3pslice
end subroutine get_Efield_raw

module subroutine get_Efield_hdf5(path, flagdirich,E0xp,E0yp,Vminx1p,Vmaxx1p,Vminx2pslice,Vmaxx2pslice,Vminx3pslice,Vmaxx3pslice)
character(*), intent(in) :: path
integer, intent(out) :: flagdirich
real(wp), dimension(:,:), intent(out) :: E0xp,E0yp,Vminx1p,Vmaxx1p
real(wp), dimension(:), intent(out) :: Vminx2pslice,Vmaxx2pslice,Vminx3pslice,Vmaxx3pslice
end subroutine get_Efield_hdf5

module subroutine get_Efield_nc4(path, flagdirich,E0xp,E0yp,Vminx1p,Vmaxx1p,Vminx2pslice,Vmaxx2pslice,Vminx3pslice,Vmaxx3pslice)
character(*), intent(in) :: path
integer, intent(out) :: flagdirich
real(wp), dimension(:,:), intent(out) :: E0xp,E0yp,Vminx1p,Vmaxx1p
real(wp), dimension(:), intent(out) :: Vminx2pslice,Vmaxx2pslice,Vminx3pslice,Vmaxx3pslice
end subroutine get_Efield_nc4

module subroutine get_precip_raw(path, Qp, E0p)
character(*), intent(in) :: path
real(wp), dimension(:,:), intent(out) :: Qp, E0p
end subroutine get_precip_raw

module subroutine get_precip_hdf5(path, Qp, E0p)
character(*), intent(in) :: path
real(wp), dimension(:,:), intent(out) :: Qp, E0p
end subroutine get_precip_hdf5

module subroutine get_precip_nc4(path, Qp, E0p)
character(*), intent(in) :: path
real(wp), dimension(:,:), intent(out) :: Qp, E0p
end subroutine get_precip_nc4

module subroutine get_neutral2_raw(path, dnO,dnN2,dnO2,dvnrho,dvnz,dTn)
character(*), intent(in) :: path
real(wp), dimension(:,:,:), intent(out) :: dnO,dnN2,dnO2,dvnrho,dvnz,dTn
end subroutine get_neutral2_raw

module subroutine get_neutral2_hdf5(path, dnO,dnN2,dnO2,dvnrho,dvnz,dTn)
character(*), intent(in) :: path
real(wp), dimension(:,:,:), intent(out) :: dnO,dnN2,dnO2,dvnrho,dvnz,dTn
end subroutine get_neutral2_hdf5

module subroutine get_neutral2_nc4(path, dnO,dnN2,dnO2,dvnrho,dvnz,dTn)
character(*), intent(in) :: path
real(wp), dimension(:,:,:), intent(out) :: dnO,dnN2,dnO2,dvnrho,dvnz,dTn
end subroutine get_neutral2_nc4

module subroutine get_neutral3_raw(path, dnOall,dnN2all,dnO2all,dvnxall,dvnrhoall,dvnzall,dTnall)
character(*), intent(in) :: path
real(wp), dimension(:,:,:), intent(out) :: dnOall,dnN2all,dnO2all,dvnxall,dvnrhoall,dvnzall,dTnall
end subroutine get_neutral3_raw

module subroutine get_neutral3_hdf5(path, dnOall,dnN2all,dnO2all,dvnxall,dvnrhoall,dvnzall,dTnall)
character(*), intent(in) :: path
real(wp), dimension(:,:,:), intent(out) :: dnOall,dnN2all,dnO2all,dvnxall,dvnrhoall,dvnzall,dTnall
end subroutine get_neutral3_hdf5

module subroutine get_neutral3_nc4(path, dnOall,dnN2all,dnO2all,dvnxall,dvnrhoall,dvnzall,dTnall)
character(*), intent(in) :: path
real(wp), dimension(:,:,:), intent(out) :: dnOall,dnN2all,dnO2all,dvnxall,dvnrhoall,dvnzall,dTnall
end subroutine get_neutral3_nc4

end interface

contains

function get_filename(path, stem) result(fn)
!! given a path and stem, find the full filename
!! assuming suffix is
character(*), intent(in) :: path
character(*), intent(in), optional :: stem

character(:), allocatable :: fn, stem1
integer :: i
logical :: exists
character(*), parameter :: suffix(3) = [character(4) :: '.h5', '.nc', '.dat' ]

if(present(stem)) then
  stem1 = '/' // stem
  !> find suffix location, if it exists
  i = index(path, '.', back=.true.)

  !> find filename if path specified
  if (i > 8) then
    !! has to be nested since no short-circuit logic
    if (path(i-7:i-1) == stem) then
      fn = path  !< assuming it's the full path
      return
    endif
  endif
else
  stem1 = ''
endif

do i = 1, size(suffix)
  fn = path // stem1 // trim(suffix(i))
  inquire(file=fn, exist=exists)
  ! print *, 'reader:get_filename: ', fn, exists
  if (exists) return
enddo

write(stderr,*) 'ERROR: filename not found in ', path
fn = ''

end function get_filename

subroutine get_simsize2(path, llon, llat)
character(*), intent(in) :: path
integer, intent(out) :: llon, llat

character(:), allocatable :: fn

fn = get_filename(path, 'simsize')

select case (get_suffix(fn))
case ('.h5')
  call get_simsize2_hdf5(fn, llon, llat)
case ('.nc')
  call get_simsize2_nc4(fn, llon, llat)
case ('.dat')
  call get_simsize2_raw(fn, llon, llat)
case default
  write(stderr,*) 'reader:get_simsize2: unknown file suffix on ' // fn
  error stop 6
end select
end subroutine get_simsize2

subroutine get_simsize3(path, lx1, lx2all, lx3all)
character(*), intent(in) :: path
integer, intent(out) :: lx1, lx2all
integer, intent(out), optional :: lx3all

character(:), allocatable :: fn

fn = get_filename(path, 'simsize')

select case (get_suffix(fn))
case ('.h5')
  call get_simsize3_hdf5(fn, lx1, lx2all, lx3all)
case ('.nc')
  call get_simsize3_nc4(fn, lx1, lx2all, lx3all)
case ('.dat')
  call get_simsize3_raw(fn, lx1, lx2all, lx3all)
case default
  write(stderr,*) 'reader:get_simsize3: unknown file suffix' // fn
  error stop 6
end select
end subroutine get_simsize3

subroutine get_grid2(path, mlonp, mlatp)
character(*), intent(in) :: path
real(wp), dimension(:), intent(out) :: mlonp, mlatp

character(:), allocatable :: fn

fn = get_filename(path, 'simgrid')

select case (get_suffix(fn))
case ('.h5')
  call get_grid2_hdf5(fn, mlonp, mlatp)
case ('.nc')
  call get_grid2_nc4(fn, mlonp, mlatp)
case ('.dat')
  call get_grid2_raw(fn, mlonp, mlatp)
case default
  write(stderr,*) 'reader:get_grid2: unknown file suffix on ' // fn
  error stop 6
end select
end subroutine get_grid2

subroutine get_Efield(path, flagdirich,E0xp,E0yp,Vminx1p,Vmaxx1p,Vminx2pslice,Vmaxx2pslice,Vminx3pslice,Vmaxx3pslice)
character(*), intent(in) :: path
integer, intent(out) :: flagdirich
real(wp), dimension(:,:), intent(out) :: E0xp,E0yp,Vminx1p,Vmaxx1p
real(wp), dimension(:), intent(out) :: Vminx2pslice,Vmaxx2pslice,Vminx3pslice,Vmaxx3pslice

character(:), allocatable :: fn

fn = get_filename(path)

select case (get_suffix(fn))
case ('.h5')
  call get_Efield_hdf5(fn, flagdirich,E0xp,E0yp,Vminx1p,Vmaxx1p,Vminx2pslice,Vmaxx2pslice,Vminx3pslice,Vmaxx3pslice)
case ('.nc')
  call get_Efield_nc4(fn, flagdirich,E0xp,E0yp,Vminx1p,Vmaxx1p,Vminx2pslice,Vmaxx2pslice,Vminx3pslice,Vmaxx3pslice)
case ('.dat')
  call get_Efield_raw(fn, flagdirich,E0xp,E0yp,Vminx1p,Vmaxx1p,Vminx2pslice,Vmaxx2pslice,Vminx3pslice,Vmaxx3pslice)
case default
  write(stderr,*) 'reader:Efield: unknown file suffix on ' // fn
  error stop 6
end select

end subroutine get_Efield

subroutine get_precip(path, Qp, E0p)
!! Qp, E0p are (llon, llat)
character(*), intent(in) :: path
real(wp), dimension(:,:), intent(out) :: Qp, E0p

character(:), allocatable :: fn

fn = get_filename(path)

select case (get_suffix(fn))
case ('.h5')
  call get_precip_hdf5(fn, Qp, E0p)
case ('.nc')
  call get_precip_nc4(fn, Qp, E0p)
case ('.dat')
  call get_precip_raw(fn, Qp, E0p)
case default
  write(stderr,*) 'reader:get_precip: unknown file suffix on ' // fn
  error stop 6
end select
end subroutine get_precip

subroutine get_neutral2(path, dnO,dnN2,dnO2,dvnrho,dvnz,dTn)
character(*), intent(in) :: path
real(wp), dimension(:,:,:), intent(out) :: dnO,dnN2,dnO2,dvnrho,dvnz,dTn

character(:), allocatable :: fn

fn = get_filename(path)

select case (get_suffix(fn))
case ('.h5')
  call get_neutral2_hdf5(fn, dnO,dnN2,dnO2,dvnrho,dvnz,dTn)
case ('.nc')
  call get_neutral2_nc4(fn, dnO,dnN2,dnO2,dvnrho,dvnz,dTn)
case ('.dat')
  call get_neutral2_raw(fn, dnO,dnN2,dnO2,dvnrho,dvnz,dTn)
case default
  write(stderr,*) 'reader:get_neutral2: unknown file suffix' // get_suffix(fn)
  error stop 6
end select
end subroutine get_neutral2

subroutine get_neutral3(path, dnOall,dnN2all,dnO2all,dvnxall,dvnrhoall,dvnzall,dTnall)
character(*), intent(in) :: path
real(wp), dimension(:,:,:), intent(out) :: dnOall,dnN2all,dnO2all,dvnxall,dvnrhoall,dvnzall,dTnall

character(:), allocatable :: fn

fn = get_filename(path)

select case (get_suffix(fn))
case ('.h5')
  call get_neutral3_hdf5(fn, dnOall,dnN2all,dnO2all,dvnxall,dvnrhoall,dvnzall,dTnall)
case ('.nc')
  call get_neutral3_nc4(fn, dnOall,dnN2all,dnO2all,dvnxall,dvnrhoall,dvnzall,dTnall)
case ('.dat')
  call get_neutral3_raw(fn, dnOall,dnN2all,dnO2all,dvnxall,dvnrhoall,dvnzall,dTnall)
case default
  write(stderr,*) 'reader:get_neutral3: unknown file suffix' // get_suffix(fn)
  error stop 6
end select
end subroutine get_neutral3

end module reader

max_stars_repo_name: jklenzing/gemini

max_issues_repo_name: jklenzing/gemini

max_forks_repo_name: jklenzing/gemini


subroutine threshold_image(image, n, threshold, output)
    ! Inputs: image, n, threshold.
    ! Output: output
    !   output(i,j) is 1 if image(i,j) > threshold and 0 otherwise.

    integer n
    real(8) threshold
    real(8), dimension(n,n) :: image, output

    !f2py intent(in) :: image, threshold
    !f2py intent(hide), depend(image) :: n = shape(image, 0)
    !f2py intent(out) output

    write(*,*) "Hello from the Fortran subroutine!"

    ! Loop through columns and rows and threshold the image.
    do j=1,n
        do i=1,n
            if (image(i,j) > threshold) then
                output(i,j) = 1.0
            else
                output(i,j) = 0.0
            end if
        end do
    end do
end subroutine


module VectorIntModule
  use KindModule, only: I4B
  use SimModule, only: ustop
  use ArrayHandlersModule, only: ExpandArray
  implicit none
  private
  public :: VectorInt

  integer(I4B), parameter :: defaultInitialCapacity = 4

  ! This is a dynamic vector type for integers
  type :: VectorInt
    integer(I4B), private, allocatable :: values(:) ! the internal array for storage
    integer(I4B) :: size                            ! the number of elements (technically this stuff should be unsigned)
    integer(I4B) :: capacity                        ! the reserved storage
  contains
    procedure, pass(this) :: init                   ! allocate memory, init size and capacity
    procedure, pass(this) :: push_back              ! adds an element at the end of the vector
    procedure, pass(this) :: at                     ! random access, unsafe, no bounds checking
    procedure, pass(this) :: at_safe                ! random access with bounds checking
    procedure, pass(this) :: clear                  ! empties the vector, leaves memory unchanged
    procedure, pass(this) :: shrink_to_fit          ! reduces the allocated memory to fit the actual vector size
    procedure, pass(this) :: destroy                ! deletes the memory
    ! private
    procedure, private, pass(this) :: expand
  end type VectorInt

contains ! module routines

  subroutine init(this, capacity)
    class(VectorInt), intent(inout) :: this
    integer(I4B), intent(in), optional :: capacity ! the initial capacity, when given

    if (present(capacity)) then
      this%capacity = capacity
    else
      this%capacity = defaultInitialCapacity
    end if

    allocate(this%values(this%capacity))
    this%size = 0

  end subroutine init

  subroutine push_back(this, newValue)
    class(VectorInt), intent(inout) :: this
    integer(I4B) :: newValue
    ! check capacity
    if (this%size + 1 > this%capacity) then
      call this%expand()
    end if

    this%size = this%size + 1
    this%values(this%size) = newValue

  end subroutine push_back

  function at(this, idx) result(value)
    class(VectorInt), intent(inout) :: this
    integer(I4B), intent(in)        :: idx
    integer(I4B)                    :: value

    value = this%values(idx)

  end function at

  function at_safe(this, idx) result(value)
    class(VectorInt), intent(inout) :: this
    integer(I4B), intent(in)        :: idx
    integer(I4B)                    :: value

    if (idx > this%size) then
      write(*,*) 'VectorInt exception: access out of bounds, index ', idx, ' exceeds actual size (', this%size, ')'
      call ustop()
    end if
    value = this%at(idx)

  end function at_safe

  subroutine clear(this)
    class(VectorInt), intent(inout) :: this

    ! really, this is all there is to it...
    this%size = 0

  end subroutine clear

  subroutine shrink_to_fit(this)
    class(VectorInt), intent(inout) :: this
    ! local
    integer(I4B), allocatable :: tempValues(:)
    integer(I4B) :: i, newSize

    if (this%size == this%capacity) then
      return
    end if

    ! store temp
    newSize = this%size
    allocate(tempValues(newSize))
    do i = 1, newSize
        tempValues(i) = this%values(i)
    end do

    ! reinit
    call this%destroy()
    call this%init(newSize)

    ! copy back
    do i = 1, newSize
        call this%push_back(tempValues(i))
    end do

  end subroutine shrink_to_fit

  subroutine destroy(this)
    class(VectorInt), intent(inout) :: this

    if (allocated(this%values)) then
      deallocate(this%values)
      this%size = 0
      this%capacity = 0
    else
      write(*,*) 'VectorInt exception: cannot delete an unallocated array'
      call ustop()
    end if

  end subroutine destroy

  ! expand the array with the given strategy, at
  ! least by 1
  subroutine expand(this)
    class(VectorInt), intent(inout) :: this
    integer(I4B) :: increment

    ! expansion strategy
    increment = this%capacity/2 + 1
    call ExpandArray(this%values, increment)
    this%capacity = this%capacity + increment

  end subroutine expand

end module VectorIntModule


text:
Compassion in Action is sort of an alias band for Brilliant Red Lights because they share the same members, but play different instruments.

Web site: http://www.brilliantredlights.com/cia.htm

 Members
  Aaron Davidson  guitar and vocals
  Jon Roberts  Moog synthesizer
  Evan Michalski  drums

 Releases
Territories (2004)  13 tracks, selfrecorded, selfreleased


text:
William King (1838  ) arrived in Yolo County in 1853. In 1856 he rented 60 acres of land and raise wheat. Later in 1869 he bought 85 acres to raise barley and wheat.

1875 he moved to the village of Davis.

He served as roadmaster.

In 1879 he began a 9 year run as justice of the peace.

In 1890 he was elected to the board of county supervisors and selected as chairman.

He married Miss Rebecca M. Montgomery on March 30, 1864. They had eight children, J. K., Daisy, Catherine Scott (married to W. H. Scott), Amanda Harrington ( married to J. B. Harrington), Thomas J., Edna Cole, Nellie Lillard (married to Samuel Lillard), and Belle.

In 1901, his home was the first in town to be wired for electric service.


! This file was generated by gpufort

module zhegvdx_gpu_kernels
  use hip
  implicit none

  interface

    subroutine launch_krnl_959801_0(grid,&
        block,&
        sharedMem,&
        stream,&
        a,&
        n,&
        z) bind(c, name="launch_krnl_959801_0")
      use iso_c_binding
      use hip
      implicit none
      type(dim3),intent(IN) :: grid
      type(dim3),intent(IN) :: block
      integer(c_int),intent(IN) :: sharedMem
      type(c_ptr),value,intent(IN) :: stream
      TODO declaration not found :: a
      integer,value :: n
      TODO declaration not found :: z
    end subroutine

    subroutine launch_krnl_959801_0_auto(sharedMem,&
        stream,&
        a,&
        n,&
        z) bind(c, name="launch_krnl_959801_0_auto")
      use iso_c_binding
      use hip
      implicit none
      integer(c_int),intent(IN) :: sharedMem
      type(c_ptr),value,intent(IN) :: stream
      TODO declaration not found :: a
      integer,value :: n
      TODO declaration not found :: z
    end subroutine

  end interface

  contains

    subroutine launch_krnl_959801_0_cpu(sharedMem,&
        stream,&
        a,&
        n,&
        z)
      use iso_c_binding
      use hip
      implicit none
      integer(c_int),intent(IN) :: sharedMem
      type(c_ptr),value,intent(IN) :: stream
      TODO declaration not found :: a
      integer,value :: n
      TODO declaration not found :: z
      integer :: j
      integer :: i
      do j = 1, N
      do i = 1, N
      if (i > j) then
         Z(i, j) = A(i, j)
      endif
      end do
      end do

    end subroutine

end module zhegvdx_gpu_kernels


text:
      MODULE iniest2_I
      INTERFACE
!...Generated by Pacific-Sierra Research 77to90  4.3E  12:36:59   2/14/04
!...Modified by Charlotte Froese Fischer
!                     Gediminas Gaigalas  10/05/17
      SUBROUTINE iniest2 (NMAX, NCF, NIV, BASIS, HMX, JCOL, IROW)
      USE vast_kind_param,ONLY: DOUBLE
      INTEGER, INTENT(IN) :: NMAX
      INTEGER, INTENT(IN) :: NCF
      INTEGER, INTENT(IN) :: NIV
      REAL(DOUBLE), DIMENSION(*) :: BASIS
      REAL(DOUBLE), DIMENSION(*), INTENT(IN) :: HMX
      INTEGER, DIMENSION(0:*), INTENT(IN) :: JCOL
      INTEGER, DIMENSION(*), INTENT(IN) :: IROW
!VAST...Calls: ALLOC, DINIT, DSPEVX, DCOPY, DALLOC
      END SUBROUTINE
      END INTERFACE
      END MODULE

max_stars_repo_name: sylas/grasp-continuum

max_issues_repo_name: sylas/grasp-continuum

max_forks_repo_name: sylas/grasp-continuum


text:
!  *******************************************************************
!  *                                                                 *
!  *                       function idamax                           *
!  *                                                                 *
!  *******************************************************************
!  Double Precision Version 06/17/77 1.01
!  Written by Jack Dongarra, Linpack
!  Modified by Gordon A. Fenton, Aug. 24, 1993
!  Latest Update: Jun 9, 1999
!
!  PURPOSE   find the index of the maximum absolute element in a double
!            precision vector of length n
!
!
!  REVISION HISTORY:
!  1.01	replaced dummy dimensions with a (*) for GNU's compiler (Jun 9/99)
!------------------------------------------------------------------------
      integer function idamax( n, dx )
      real*8 dx(*), dmax
      integer i, n

      if( n .le. 0 ) then
         idamax = 0
         return
      endif
      idamax = 1
      if( n .eq. 1 )return

      dmax = dabs(dx(1))
      do 10 i = 2, n
         if( dabs(dx(i)) .gt. dmax ) then
            idamax = i
            dmax   = dabs(dx(i))
         endif
  10  continue

      return
      end

max_stars_repo_name: leemargetts/PoreFEM

max_issues_repo_name: leemargetts/PoreFEM

max_forks_repo_name: leemargetts/PoreFEM


text:
c
c ------------------------------------------------------------
c
       subroutine upbnd(listbc,val,nvar,mitot,mjtot,
     1                  maxsp,mptr)
c     1                  maxsp,iused,mptr)

      implicit double precision (a-h,o-z)

      include  "call.i"

       dimension val(mitot,mjtot,nvar),listbc(5,maxsp),
     1           iused(mitot,mjtot)

       iaddaux(i,j) = locaux + i-1 +  mitot*(j-1)
     1                 + mitot*mjtot*(mcapa-1)

c
c :::::::::::::::::::::::::::: UPBND :::::::::::::::::::::::::::::
c We now correct the coarse grid with the flux differences stored
c with each of the fine grids. We use an array   iused
c to indicate whether the flux has been updated or not for that zone.
c iused(i,j) = sum from (l=1,4) i(l)*2**(l-1), where i(l) = 1 if the
c flux for the  l-th side of the (i,j)-th cell has already been
c updated, and i(l) = 0 if not.

c if there is a capacity fn. it needs to be included in update formula
c indicated by mcapa not zero (is index of capacity fn.)
c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
c

      do 10 j=1,mjtot
      do 10 i=1,mitot
         iused(i,j) = 0.
 10   continue

      locaux = node(storeaux,mptr)
      levc   = node(nestlevel,mptr)
      area   = hxposs(levc)*hyposs(levc)

      if (uprint) write(outunit,*)" upbnding grid ",mptr

      do 40 ispot = 1,maxsp
         icrse = listbc(1,ispot)
         if (icrse.eq.0) go to 99

         jcrse = listbc(2,ispot)
         iside = listbc(3,ispot)
c        continue to use iside1/norm for debugging, but should soon delete
c        this if/then/else block needed due to new categories corresponding
c        to mapped bcs. should still only have one update per side of coarse cell though
         if (iside .lt. 5) then
           iside1 = iside
         elseif (iside .eq. 5) then
           iside1 = 2
         else  ! iside is 6
           iside1 = 4
         endif
         norm = 2**(iside1-1)
         iflag =iused(icrse,jcrse)/norm
         if (mod(iflag,2).eq.1) then
            write(6,*)" ***  double flux update CAN happen in upbnd ***"
            go to 40
         endif
         mkid = listbc(4,ispot)
         kidlst = node(ffluxptr,mkid)
         lkid = listbc(5,ispot)
c         if (mod(iside,4).gt.1) then
c         modified to include other side options
          if (iside .eq. 2 .or. iside .eq. 3 .or. iside .eq. 6) then
c           (iside .eq. 2 .or. iside .eq. 3)
            sgnm = -1.
         else
c           (iside .eq. 4 .or. iside .eq. 1)
            sgnm = 1.
         endif

c        ## debugging output
         if (uprint) then
           write(outunit,101) icrse,jcrse,
     .         (val(icrse,jcrse,ivar),ivar=1,nvar)
 101       format(" old ",1x,2i4,4e15.7)
         endif

         if (mcapa .gt. 0) then
c            # capacity array:  need to divide by capa in each cell.
c            # modify sgnm which is reset for each grid cell.
c            # Note capa is stored in aux(icrse,jcrse,mcapa)
             sgnm = sgnm / alloc(iaddaux(icrse,jcrse))
         endif

         do 20 ivar = 1,nvar
            val(icrse,jcrse,ivar) = val(icrse,jcrse,ivar) +
     1      sgnm*alloc(kidlst+nvar*(lkid-1)+ivar-1)/area
 20      continue
         iused(icrse,jcrse) = iused(icrse,jcrse) + norm

c        ## debugging output
         if (uprint) then
           write(outunit,102) mkid,
     .         (val(icrse,jcrse,ivar),ivar=1,nvar)
 102       format(" new ","(grid",i3,")",4e15.7)
         endif

 40   continue
c
 99   return
      end

max_issues_repo_name: che-wenchao/D-Claw

max_forks_repo_name: che-wenchao/D-Claw


text:
I am a fourth year majoring in Human Development.  I have four wonderful pet RATS: NoyNoy & ToyToy (Each one year old) and Skittles (2.5 years old, RIP 2/13/06) and Skuttles (2 years old, RIP 7/27/05)

20060127 23:15:12 nbsp You should add pictures of the ratties. ) Users/MichaelGiardina

20060213 20:00:41 nbsp Hi Krystal, my name is Patty McCoy.  I went to school at Birch Lane School Birch Lane Elementary until 1970 and Im interested in locating a guy I was in school with there named Gilbert Villanueva...any realtion? Users/PattyMcCoy

20060224 21:26:17 nbsp Cuuuuute Users/MichaelGiardina


text:
      subroutine rot3i (lxp1, mxp1, ileg)
      implicit double precision (a-h,o-z)

c     input:  lxp1, mxp1, ileg (lmax+1, mmax+1)
c             also beta(ileg) used from common /pdata/
c     output: dri(...ileg) in common /rotmat/

c     subroutine rot3 calculates rotation matrices for l = 0,lxp1-1

c     subroutine rot3 calculates the beta dependence of rotation
c     matrix elements using recursion of an iterated version of
c     formula (4.4.1) in edmonds.
c
c     first written:(september 17,1986) by j. mustre
c     version 2  (17 sep 86)
c     version 3  (22 feb 87) modified by j. rehr
c     version for genfmt, modified by s. zabinsky, Sept 1991
c     Initialized dri0.  Some elements may be used before being
c        initialized elsewhere -- rot3i needs to be carefully
c        checked.  S. Zabinsky, April 1993
c
c******************** warning******************************************
c     ltot must be at least lxp1 or overwriting will occur
c     nmax must be at least nm or overwriting will occur
c----------------------------------------------------------------------
c     notation dri0(l,m,n) =  drot_i(l'm'n')
c     l = l'+1, n' = n-l, m' = m-l, primes denoting subscripts
c     thus dri0(1,1,1) corresponds to the rotation matrix with
c     l' = 0, and n' and m' = 0; dri0(3,5,5) : l' = 2,n' = 2,m' = 2.
c--------------------------------------------------------------------

      include 'dim.h'
      include 'rotmat.h'
      include 'pdata.h'
c     dri0 is larger than needed for genfmt, but necessary for
c     this calculation algorithm.  Copy result into smaller
c     dri arrays (in common) at end of this routine.
      dimension  dri0 (ltot+1, 2*ltot+1, 2*ltot+1)

c     initialize dri0
      do 200 il = 1, ltot+1
         do 200 im = 1, 2*ltot+1
            do 200 in = 1, 2*ltot+1
               dri0(il,im,in) = 0
  200 continue

      nm = mxp1
      ndm = lxp1+nm-1
      xc = cos(beta(ileg)/2)
      xs = sin(beta(ileg)/2)
      s = sin(beta(ileg))
      dri0(1,1,1) = 1
      dri0(2,1,1) = xc**2
      dri0(2,1,2) = s/sqrt(2.0d0)
      dri0(2,1,3) = xs**2
      dri0(2,2,1) = -dri0(2,1,2)
      dri0(2,2,2) = cos(beta(ileg))
      dri0(2,2,3) = dri0(2,1,2)
      dri0(2,3,1) = dri0(2,1,3)
      dri0(2,3,2) = -dri0(2,2,3)
      dri0(2,3,3) = dri0(2,1,1)
      do 30  l = 3, lxp1
         ln = 2*l - 1
         lm = 2*l - 3
         if (ln .gt. ndm)  ln = ndm
         if (lm .gt. ndm)  lm = ndm
         do 20  n = 1, ln
            do 10  m = 1, lm
               t1 = (2*l-1-n) * (2*l-2-n)
               t = (2*l-1-m) * (2*l-2-m)
               f1 = sqrt (t1/t)
               f2 = sqrt ((2*l-1-n) * (n-1) / t)
               t3 = (n-2) * (n-1)
               f3 = sqrt(t3/t)
               dlnm = f1 * xc**2 * dri0(l-1,n,m)
               if (n-1 .gt. 0) dlnm = dlnm - f2*s*dri0(l-1,n-1,m)
               if (n-2 .gt. 0) dlnm = dlnm + f3*xs**2*dri0(l-1,n-2,m)
               dri0(l,n,m) = dlnm
               if (n .gt. (2*l-3))
     1            dri0(l,m,n) = (-1)**(n-m) * dri0(l,n,m)
   10       continue
            if (n .gt. (2*l-3)) then
               dri0(l,2*l-2,2*l-2) = dri0(l,2,2)
               dri0(l,2*l-1,2*l-2) = -dri0(l,1,2)
               dri0(l,2*l-2,2*l-1) = -dri0(l,2,1)
               dri0(l,2*l-1,2*l-1) = dri0(l,1,1)
            endif
   20    continue
   30 continue
   40 continue

c-----test sum rule on d
c     open (19,file='rotmat.dat',status='new',carriagecontrol='list')
c     write(19,*)  ' l, m, sum'
c     write(19,*) ' (dri0(il,im,in),in = 1,ln)'
c     do 70 il = 1,lxp1
c        l = il-1
c        ln = 2*l+1
c        if(ln.gt.ndm) ln = ndm
c        do 37 im = 1,ln
c           sum = 0
c           do 50 in = 1,ln
c              m = im-il
c              term = dri0(il,im,in)
c  50       sum = sum+term**2
c           write(19,60) l,m,sum
c           write(19,62) (dri0(il,im,in),in = 1,ln)
c  60       format(2i3,e30.20)
c  62       format(5e14.6)
c  70 continue
c     close(19)
c-----end test------------------------

c     Copy result into dri(...ileg) in /rotmat/ (zero it first...)
      do 90  il = 1, ltot+1
         do 90  m1 = 1, 2*mtot+1
            do 90  m2 = 1, 2*mtot+1
               dri(il,m1,m2,ileg) = 0
   90 continue

      do 120  il = 1, lxp1
         mx = min (il-1, mxp1-1)
         do 110  m1 = -mx, mx
            do 100  m2 = -mx, mx
               dri(il,m1+mtot+1,m2+mtot+1,ileg)=dri0(il,m1+il,m2+il)
  100       continue
  110    continue
  120 continue

      return
      end

max_stars_repo_name: keechul/ifeffit

max_stars_repo_licenses: Naumen Condor-1.1 MS-PL

max_issues_repo_name: bruceravel/ifeffit

max_issues_repo_licenses: Naumen Condor-1.1 MS-PL

max_forks_repo_name: bruceravel/ifeffit

max_forks_repo_licenses: Naumen Condor-1.1 MS-PL


! sdl2_joystick.f90
!
! Fortran 2008 interfaces to SDL_gamecontroller.h.
!
! Author:  Philipp Engel
! GitHub:  https://github.com/interkosmos/fortran-sdl2/
! Licence: ISC
module sdl2_gamecontroller
    use, intrinsic :: iso_c_binding, only: c_int, c_ptr
    use :: sdl2_stdinc
    implicit none
    private

    public :: sdl_game_controller_close
    public :: sdl_game_controller_open
    public :: sdl_is_game_controller

    interface
        ! SDL_*SDL_GameControllerOpen(int joystick_index)
        function sdl_game_controller_open(joystick_index) bind(c, name='SDL_GameControllerOpen')
            import :: c_int, c_ptr
            integer(kind=c_int), intent(in), value :: joystick_index
            type(c_ptr)                            :: sdl_game_controller_open
        end function sdl_game_controller_open

        ! SDL_bool SDL_IsGameController(int joystick_index)
        function sdl_is_game_controller(joystick_index) bind(c, name='SDL_IsGameController')
            import :: c_int, sdl_bool
            integer(kind=c_int), intent(in), value :: joystick_index
            integer(kind=sdl_bool)                 :: sdl_is_game_controller
        end function sdl_is_game_controller

        ! void SDL_GameControllerClose(SDL_GameController *gamecontroller)
        subroutine sdl_game_controller_close(game_controller) bind(c, name='SDL_GameControllerClose')
            import :: c_ptr
            type(c_ptr), intent(in), value :: game_controller
        end subroutine sdl_game_controller_close
    end interface
end module sdl2_gamecontroller


text:
      SUBROUTINE XSTORE(COUT,CON,MWORDS)
C$$$  SUBPROGRAM DOCUMENTATION BLOCK
C                .      .    .                                       .
C SUBPROGRAM:    XSTORE      STORES A CONSTANT VALUE INTO AN ARRAY
C   PRGMMR: KEYSER           ORG: W/NMC22    DATE: 07-02-92
C
C ABSTRACT: STORES AN 4-BYTE (FULLWORD) VALUE THROUGH CONSECUTIVE
C   STORAGE LOCATIONS.  (MOVING IS ACCOMPLISHED WITH A DO LOOP.)
C
C PROGRAM HISTORY LOG:
C   92-07-02  D. A. KEYSER (W/NMC22)
C   93-03-29  R.E.JONES   ADD SAVE STATEMENT
C   93-04-15  R.E.JONES   CHANGES FOR MICROSOFT FORTRAN 5.0
C
C USAGE:    CALL XSTORE(COUT,CON,MWORDS)
C   INPUT ARGUMENT LIST:
C     CON      - CONSTANT TO BE STORED INTO "MWORDS" CONSECUTIVE
C                FULLWORDS BEGINNING WITH "COUT" ARRAY
C     MWORDS   - NUMBER OF FULLWORDS IN "COUT" ARRAY TO STORE "CON";
C                MUST BE .GT. ZERO (NOT CHECKED FOR THIS)
C
C   OUTPUT ARGUMENT LIST:      (INCLUDING WORK ARRAYS)
C     COUT     - STARTING ADDRESS FOR ARRAY OF "MWORDS" FULLWORDS
C                SET TO THE CONTENTS OF THE VALUE "CON"
C
C REMARKS: THE VERSION OF THIS SUBROUTINE ON THE HDS COMMON LIBRARY
C   IS NAS-SPECIFIC SUBR. WRITTEN IN ASSEMBLY LANG. TO ALLOW FAST
C   COMPUTATION TIME.  SUBR. PLACED IN CRAY W3LIB TO ALLOW CODES TO
C   COMPILE ON BOTH THE HDS AND CRAY MACHINES.
C
C ATTRIBUTES:
C   LANGUAGE: MICROSOFT FORTRAN 5.0 OPTIMIZING COMPILER
C   MACHINE:  IBM PC, AT, PS/2, 386, 486, 586, CLONES.
C
C$$$
C
      DIMENSION  COUT(*)
C
      SAVE
C
      DO 1000  I = 1,MWORDS
        COUT(I) = CON
1000  CONTINUE
C
      RETURN
      END

max_stars_repo_name: maxinye/laps-mirror

max_stars_repo_licenses: Intel Unlicense OLDAP-2.2.1 NetCDF

max_issues_repo_name: longwosion/laps-mirror

max_issues_repo_licenses: Intel NetCDF OLDAP-2.2.1 Unlicense

max_forks_repo_name: longwosion/laps-mirror

max_forks_repo_licenses: Intel NetCDF OLDAP-2.2.1 Unlicense


text:
!
!     CalculiX - A 3-dimensional finite element program
!              Copyright (C) 1998-2018 Guido Dhondt
!
!     This program is free software; you can redistribute it and/or
!     modify it under the terms of the GNU General Public License as
!     published by the Free Software Foundation(version 2);
!
!
!     This program is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!     GNU General Public License for more details.
!
!     You should have received a copy of the GNU General Public License
!     along with this program; if not, write to the Free Software
!     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
!
      subroutine changeplastics(inpc,textpart,imat,ntmat_,npmat_,
     &        plicon,nplicon,plkcon,nplkcon,istep,istat,n,iline,ipol,
     &        inl,ipoinp,inp,ipoinpc,nelcon,ier)
!
!     reading the input deck: *CHANGE PLASTIC
!
      implicit none
!
      logical iso
!
      character*1 inpc(*)
      character*132 textpart(16)
!
      integer imat,ntmat_,ntmat,npmat_,npmat,istep,nelcon(2,*),
     &  n,key,i,nplicon(0:ntmat_,*),nplkcon(0:ntmat_,*),istat,
     &  iline,ipol,inl,ipoinp(2,*),inp(3,*),ipoinpc(0:*),ier
!
      real*8 plicon(0:2*npmat_,ntmat_,*),plkcon(0:2*npmat_,ntmat_,*),
     & temperature
!
      iso=.true.
!
      ntmat=0
      npmat=0
!
      if(istep.lt.1) then
         write(*,*) '*ERROR reading *CHANGE PLASTIC: *CHANGE PLASTIC'
         write(*,*) '       should only be used within a STEP'
         ier=1
         return
      endif
!
      if((nelcon(1,imat).ne.-51).and.
     &   (nelcon(1,imat).ne.-52)) then
         write(*,*) '*ERROR reading *CHANGE PLASTIC: *CHANGE PLASTIC'
         write(*,*) '       can only be used to change the plastic'
         write(*,*) '       definition of an elastically isotropic'
         write(*,*) '       material with *PLASTIC data'
         ier=1
         return
      endif
!
      do i=2,n
         if(textpart(i)(1:10).eq.'HARDENING=') then
            if(textpart(i)(11:19).eq.'KINEMATIC') then
               iso=.false.
            elseif(textpart(i)(11:18).eq.'COMBINED') then
               write(*,*) '*ERROR reading *CHANGE PLASTIC'
               write(*,*) '       combined hardening is not allowed'
               ier=1
               return
            elseif(textpart(i)(11:14).eq.'USER') then
               write(*,*) '*ERROR reading *CHANGE PLASTIC'
               write(*,*) '       parameter USER is not allowed'
               ier=1
               return
            endif
            exit
         else
            write(*,*)
     &     '*WARNING reading *CHANGE PLASTIC: parameter not recognized:'
            write(*,*) '         ',
     &                 textpart(i)(1:index(textpart(i),' ')-1)
            call inputwarning(inpc,ipoinpc,iline,
     &"*CHANGE PLASTIC%")
         endif
      enddo
!
      if(iso) then
!
!        isotropic hardening coefficients
!
         do
            call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl,
     &           ipoinp,inp,ipoinpc)
            if((istat.lt.0).or.(key.eq.1)) exit
            read(textpart(3)(1:20),'(f20.0)',iostat=istat) temperature
            if(istat.gt.0) then
               call inputerror(inpc,ipoinpc,iline,
     &              "*CHANGE PLASTIC%",ier)
               return
            endif
!
!           first temperature
!
            if(ntmat.eq.0) then
               npmat=0
               ntmat=ntmat+1
               if(ntmat.gt.ntmat_) then
                  write(*,*) '*ERROR reading *CHANGE PLASTIC:'
                  write(*,*) '       more temperature data points'
                  write(*,*) '       than underneath the *PLASTIC card'
                  ier=1
                  return
               endif
               nplicon(0,imat)=ntmat
               plicon(0,ntmat,imat)=temperature
!
!           new temperature
!
            elseif(plicon(0,ntmat,imat).ne.temperature) then
               npmat=0
               ntmat=ntmat+1
               if(ntmat.gt.ntmat_) then
                  write(*,*) '*ERROR reading *CHANGE PLASTIC:'
                  write(*,*) '       more temperature data points'
                  write(*,*) '       than underneath the *PLASTIC card'
                  ier=1
                  return
               endif
               nplicon(0,imat)=ntmat
               plicon(0,ntmat,imat)=temperature
            endif
            do i=1,2
               read(textpart(i)(1:20),'(f20.0)',iostat=istat)
     &              plicon(2*npmat+i,ntmat,imat)
               if(istat.gt.0) then
                  call inputerror(inpc,ipoinpc,iline,
     &                 "*CHANGE PLASTIC%",ier)
                  return
               endif
            enddo
            npmat=npmat+1
            if(npmat.gt.npmat_) then
                  write(*,*) '*ERROR reading *CHANGE PLASTIC:'
                  write(*,*) '       more stress versus equivalent'
                  write(*,*) '       plastic strain data points'
                  write(*,*) '       than underneath the *PLASTIC card'
               ier=1
               return
            endif
            nplicon(ntmat,imat)=npmat
         enddo
      else
!
!        kinematic hardening coefficients
!
         do
            call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl,
     &           ipoinp,inp,ipoinpc)
            if((istat.lt.0).or.(key.eq.1)) exit
            read(textpart(3)(1:20),'(f20.0)',iostat=istat) temperature
            if(istat.gt.0) then
               call inputerror(inpc,ipoinpc,iline,
     &              "*CHANGE PLASTIC%",ier)
               return
            endif
!
!           first temperature
!
            if(ntmat.eq.0) then
               npmat=0
               ntmat=ntmat+1
               if(ntmat.gt.ntmat_) then
                  write(*,*) '*ERROR reading *CHANGE PLASTIC:'
                  write(*,*) '       more temperature data points'
                  write(*,*) '       than underneath the *PLASTIC card'
                  ier=1
                  return
               endif
               nplkcon(0,imat)=ntmat
               plkcon(0,ntmat,imat)=temperature
!
!           new temperature
!
            elseif(plkcon(0,ntmat,imat).ne.temperature) then
               npmat=0
               ntmat=ntmat+1
               if(ntmat.gt.ntmat_) then
                  write(*,*) '*ERROR reading *CHANGE PLASTIC:'
                  write(*,*) '       more temperature data points'
                  write(*,*) '       than underneath the *PLASTIC card'
                  ier=1
                  return
               endif
               nplkcon(0,imat)=ntmat
               plkcon(0,ntmat,imat)=temperature
            endif
            do i=1,2
               read(textpart(i)(1:20),'(f20.0)',iostat=istat)
     &              plkcon(2*npmat+i,ntmat,imat)
               if(istat.gt.0) then
                  call inputerror(inpc,ipoinpc,iline,
     &                 "*CHANGE PLASTIC%",ier)
                  return
               endif
            enddo
            npmat=npmat+1
            if(npmat.gt.npmat_) then
               write(*,*) '*ERROR reading *CHANGE PLASTIC:'
               write(*,*) '       more stress versus equivalent'
               write(*,*) '       plastic strain data points'
               write(*,*) '       than underneath the *PLASTIC card'
               ier=1
               return
            endif
            nplkcon(ntmat,imat)=npmat
         enddo
      endif
!
      if(ntmat.eq.0) then
         write(*,*) '*ERROR reading *CHANGE PLASTIC:'
         write(*,*) '       *CHANGE PLASTIC card without data'
         ier=1
         return
      endif
!
      return
      end

max_stars_repo_name: alleindrach/calculix-desktop

max_issues_repo_name: alleindrach/calculix-desktop

max_forks_repo_name: alleindrach/calculix-desktop


text:
! This is the main program
program main

    use global_variables
    use user_interface
    use results

    implicit none
    double precision :: Nc_final, dpmg_final, sigmag_final, vapor_final, t_final_output
    double precision :: Nc_initial, dpmg_initial, sigmag_initial, vapor_initial
    !double precision :: num_zero_0 = 1e11, dpg_zero_0 = 10.0e-9, sigma_g_zero_0 = 1.5, vapor_zero_0 = 1.85717257286e+16
    !double precision :: num_zero_0 = 0e0, dpg_zero_0 = 0.0e-9, sigma_g_zero_0 = 1.0, vapor_zero_0 = 1.85717257286e+18
    integer :: i

        call  read_inputfile()   ! To read the inputs
        call  read_temperature_profile()

        Nc_initial = num_zero_0
        dpmg_initial = dpg_zero_0
        sigmag_initial = sigma_g_zero_0

        OPEN (UNIT=10,FILE='output.txt')
90010   FORMAT (ES11.2E3,5ES11.2E3)
        WRITE (10,90010) t0, Temperature(1), Nc_initial*1e-6, dpmg_initial*1e9, sigmag_initial, vapor_initial
        do i = 1,nlines-1
            t_final_output = time(i)
            T = Temperature(i)+273.0
            vapor_initial = vapor_conc(i)
            call  output(Nc_initial , dpmg_initial , sigmag_initial, vapor_initial, Nc_final, dpmg_final, sigmag_final, vapor_final,t_final_output)
            !call output(Nc_final, dpmg_final, sigmag_final, vapor_final, t_final_output)          ! To perform all the calculations, and write the results
            !print*, t_final_output, Nc_final, dpmg_final, sigmag_final, vapor_final

            if (Nc_final < 0.0) then
                Nc_final =0.0
            end if

            if (dpmg_final < 0.0) then
                dpmg_final =0.0
            end if

            if (sigmag_final < 1.0) then
                sigmag_final =1.0
            end if

            if (vapor_final < 0.0) then
                vapor_final =0.0
            end if

            Nc_initial = Nc_final*1e6
            dpmg_initial = dpmg_final*1e-9
            sigmag_initial = sigmag_final

            !print*,time(i), T, Nc_final, dpmg_final, sigmag_final, vapor_final

            WRITE (10,90010) z_pos(i+1), T, Nc_final, dpmg_final, sigmag_final, vapor_final

        end do
        close(unit = 10)
end program main

max_stars_repo_path: AAQRL-ADM/Moment model/Fortran/main.f90

max_issues_repo_path: AAQRL-ADM/Moment model/Fortran/main.f90

max_forks_repo_path: AAQRL-ADM/Moment model/Fortran/main.f90