text:
! Modified 14/01/22
    ! FOBOS-ME4.f90 takes four epochs of observation of a two body system and fits
    ! the orbital parameters (semi-major axis, eccentricity, inclination, orientation
    ! and mean anomaly). The code is currently set up to read in three separations (in milliarcseconds,
    ! column 2) and three position angles (in degrees, column 4), plus their associated respective errors,
    ! (columns 3 and 5) from a text file called "input". The input file should have the format as shown below.

    ! 0.0000     sep1 (mas)   sep1_err (mas)   pa1 (deg)   pa1_err (deg)
    ! t1 (yrs)   sep2 (mas)   sep2_err (mas)   pa2 (deg)   pa2_err (deg)
    ! t2 (yrs)   sep3 (mas)   sep3_err (mas)   pa3 (deg)   pa3_err (deg)

    ! t1 is the time between the first and second epochs of observation (in years) and
    ! t2 is the time between the second and third epochs of observation (in years)

    ! The target name, which affects the name of the output files, is given in a text file
    ! called "sys-name". The masses and distance to the objects are contained in the file
    ! "mass-dist", which has the format

    ! M1 M1_err  M2  M2_err  d   d_err

    ! All masses will be given in solar masses and distances in parsecs
! End of information

PROGRAM binariesme
    IMPLICIT NONE
    INCLUDE "omp_lib.h"

    INTEGER:: n,i,idum,kdum,match,nmatch
    INTEGER:: nt,my_id,count,loop
    REAL:: M(1:2),M_err(1:2),dtobs1,dtobs2,dtobs3,d,d_err
    REAL:: obs(1:7),prec(1:7),params(1:5),coords(1:6),sep(1:7)
    REAL:: dts(1:4),seps(1:4),sep_err(1:4),pa(1:4),pa_err(1:4)
    REAL:: dmoved1,dmoved2,dmoved3
    REAL:: pi,au,year,msun,G,d2r
    REAL:: amin,amax,vobs,vobs1,vobs2,vobs3
    REAL:: t_run,texit,er
    REAL:: ran3
    INTEGER,ALLOCATABLE:: rinit(:)
    !Params
    REAL:: obsdata(1:5,1:4)
    REAL:: properties(1:6)
    CHARACTER(len=1024):: name

    !Date and time and filenames
    CHARACTER(len=1024):: fname1
    CHARACTER(*), PARAMETER:: fpath="./"
    CHARACTER(10):: dat,tim
    INTEGER(8) :: dt(1:8)

    !Required to make ran3 parallelisable
    INTEGER :: iff,inext,inextp,ma(55)
    COMMON /ranvar/ iff,inext,inextp,ma
    !$OMP THREADPRIVATE(/ranvar/)

    OPEN(unit=1,file=fpath//"input",status="old")
    OPEN(unit=2,file=fpath//"mass-dist",status="old")
    OPEN(unit=3,file=fpath//"sys-name",status="old")

    !Read the values into their arrays
    READ(1,*) obsdata
    READ(2,*) properties
    READ(3,*) name
    !Close files after reading
    CLOSE(1)
    CLOSE(2)
    CLOSE(3)

    !Simulation parameters
    n=1e9           !Number of iterations of the main do loop before checking for matches
    kdum=625        !Random seed
    nmatch=5000     !Number of required matches before terminating the simulation
    texit=3600.     !Time after which to exit the simulation if too few matches have been found (seconds)

    nt=omp_get_max_threads()    !Set number of threads to run on as maximum

    !Constants
    pi  = 4.*ATAN(1.)
    au  = 1.5E11
    year= 365.25*24.*3600.
    msun= 1.989E30
    G   = 6.67E-11
    d2r = (0.5*pi)/90

    !Counters
    count=0
    loop =0

    !Unpacking the input data
    dts     = obsdata(1,:)
    seps    = obsdata(2,:)
    sep_err = obsdata(3,:)
    pa      = obsdata(4,:)*d2r
    pa_err  = obsdata(5,:)*d2r

    !Masses and distance
    M(1)    = properties(1)
    M_err(1)= properties(2)
    M(2)    = properties(3)
    M_err(2)= properties(4)
    d       = properties(5)
    d_err   = properties(6)

    !Assign separations and angles to points in the obs array
    obs(1) = seps(1)*(1.E-3)*d
    obs(2) = seps(2)*(1.E-3)*d
    obs(3) = pa(2)-pa(1)
    obs(4) = seps(3)*(1.E-3)*d
    obs(5) = pa(3)-pa(2)
    obs(6) = seps(4)*(1.E-3)*d
    obs(7) = pa(4)-pa(3)

    !Factor to multiply the errors by (er=1. means use 1 sigma errors,
    !er=2. would use 2 sigma errors, etc.)
    er=1.

    !Observational errors
    prec(1) = obs(1) * SQRT((sep_err(1)/seps(1))**2. + (d_err/d)**2.)*er
    prec(2) = obs(2) * SQRT((sep_err(2)/seps(2))**2. + (d_err/d)**2.)*er
    prec(3) = (pa_err(1)+pa_err(2))*er
    prec(4) = obs(4) * SQRT((sep_err(3)/seps(3))**2. + (d_err/d)**2.)*er
    prec(5) = (pa_err(2)+pa_err(3))*er
    prec(6) = obs(6) * SQRT((sep_err(4)/seps(4))**2. + (d_err/d)**2.)*er
    prec(7) = (pa_err(3)+pa_err(4))*er

    !Time between observations
    dtobs1 = dts(2)-dts(1)
    dtobs2 = dts(3)-dts(2)
    dtobs3 = dts(4)-dts(3)

    !Calculating the on sky velocities
    dmoved1 = SQRT(obs(1)**2+obs(2)**2-(2.*obs(1)*obs(2)*COS(obs(3))))
    dmoved2 = SQRT(obs(2)**2+obs(4)**2-(2.*obs(2)*obs(4)*COS(obs(5))))
    dmoved3 = SQRT(obs(4)**2+obs(6)**2-(2.*obs(4)*obs(6)*COS(obs(7))))

    !Minimum semi-major axis
    amin = MIN(obs(1),obs(2),obs(4),obs(6))*0.5

    !Observed velocities
    vobs1 = (dmoved1*au)/(dtobs1*year)
    vobs2 = (dmoved2*au)/(dtobs2*year)
    vobs3 = (dmoved3*au)/(dtobs3*year)

    !Maximum semi-major axis
    vobs = MIN(vobs1,vobs2,vobs3)
    amax = 20.*G*(M(1)+M(2))*msun/(vobs*vobs*au)

    !Getting the output file set up
    WRITE(fname1,"(A7,A9)") "params-",name
    OPEN(unit=1,file=fpath//trim(fname1),status="new")

    !Params file preamble
        CALL DATE_AND_TIME(DATE=dat,TIME=tim,VALUES=dt)
        WRITE(1,'(A)') "The results contained in this file are for testing FOBOS-ME: Few Observation Binary Orbit Solver"
        WRITE(1,'(A,A9)') "(Multiple Epoch) on the binary/exoplanetary system ",trim(name)
        WRITE(1,'(A,I2,A,I2,A,I4,A,I2,A,I2,A,I2)') "This simulation ran on  ",dt(3),"/",dt(2),"/",dt(1), &
        "  at  ",dt(5),":",dt(6),":",dt(7)
        WRITE(1,'(A,E8.1,A)') "Tested ", REAL(n),"fake systems per loop"
        WRITE(1,'(A,I3,A,I1,A)') "The random seed value was ",kdum, " and it ran on ",nt," threads."
        WRITE(1,'(A,F4.2,A,F4.2,A,F4.2)') "Masses: m1=",M(1),", m2=",M(2)
    !End preamble

    !Random seeds for different threads
    CALL omp_set_num_threads(nt)
    ALLOCATE(rinit(nt))
    do i=1,nt
        rinit(i)=int(ran3(kdum)*1000000)
    end do

    !Main do loop
1   continue
    !$OMP PARALLEL PRIVATE(my_id,idum,params,match,sep)
    if (loop == 0) then
        my_id= omp_get_thread_num()+1
        idum = rinit(my_id)
    end if
    !$OMP DO
    do i=1,n
        CALL fakesystem(amin,amax,dtobs1,dtobs2,dtobs3,vobs,M,obs,prec,idum,params,match,coords,sep)
        if (match == 1) then
            !$OMP CRITICAL
            count=count+1
            WRITE(1,*) params
            WRITE(2,*) sep
            WRITE(3,*) coords
            !$OMP END CRITICAL
        end if
    end do
    !$OMP END DO
    !$OMP END PARALLEL

    !Write to screen with how many loops of n interations have been found
    !and the number of matches so far
    loop = loop + 1
    WRITE(6,*) "Looped:",loop,"Number of matches:",count
    call cpu_time(t_run)

    !Exit the simulation with no matches if it takes too long
    if (t_run/nt < texit) then
       if (count<nmatch) goto 1
    end if

    !If simulation stops because of the time write to screen
    if (count == 0) then
        WRITE(6,*) "Timeout"
    end if

    !End of program admin
    DEALLOCATE(rinit)
    WRITE(1,*) "Number of loops:",loop," and number of matches:",count
    WRITE(1,"(A,F8.2,A)") "Run time =", t_run/nt, "s"
    WRITE(6,"(A,F8.2,A)") "Run time =", t_run/nt, "s"

    CLOSE(1)

END PROGRAM binariesme

!----------------------------------------------------------------
!Projects the binary randomly on the sky
SUBROUTINE makeobsbinary(a,ecc,inc,phi,mean,sky)
    IMPLICIT NONE

    REAL, INTENT(in) :: a,ecc,inc,phi,mean
    REAL, INTENT(out):: sky(1:2)
    REAL:: pi,Ee,old,theta,rt,r(1:3),x,y

    pi=4.*ATAN(1.)

    !Calculate eccentric anomoly
    Ee=mean
    old=mean
    do
        old=Ee
        Ee=Ee-(Ee-ecc*SIN(Ee)-mean)/(1.-ecc*COS(Ee))
        if (old-Ee<1.E-6) exit
    end do

    theta=2.*ATAN(SQRT((1.+ecc)/(1.-ecc))*TAN(0.5*Ee))
    rt = a*(1.-ecc*ecc)/(1.+ecc*COS(theta))

    !Get x and y positions
    r(1:3)=0.
    r(1)=rt*COS(theta)
    r(2)=rt*SIN(theta)
    !Orientation around z axis
    x=r(1);y=r(2)
    r(1)=x*COS(phi)-y*SIN(phi)
    r(2)=x*SIN(phi)+y*COS(phi)
    !Modify inclination
    x=r(1);y=r(2)
    r(1)=x
    r(2)=y*COS(inc)
    r(3)=y*SIN(inc)
    !Positions on the sky
    sky(1)=r(1)
    sky(2)=r(3)
END SUBROUTINE makeobsbinary

!-----------------------------------------------------------------
SUBROUTINE fakesystem(amin,amax,dtobs1,dtobs2,dtobs3,vobs,M,obs,prec,idum,params,match,coords,sep)
    IMPLICIT NONE

    REAL, INTENT(in) :: amin,amax,dtobs1,dtobs2,dtobs3,vobs,M(1:2),obs(1:7),prec(1:7)
    REAL, INTENT(out):: params(1:5),coords(1:6),sep(1:7)
    INTEGER, INTENT(in) :: idum
    INTEGER, INTENT(out):: match
    INTEGER:: i,itemp,match2,match3,match4
    REAL:: a,ecc,inc,phi
    REAL:: mean1,mean2,mean3,mean4
    REAL:: dM1,dM2,dM3
    REAL:: sky1(1:2),sky2(1:2),sky3(1:2),sky4(1:4)
    REAL:: pi,G,msun,au,T,f(1:2),v,rt

    REAL:: ran3
    EXTERNAL:: ran3

    !Constants and zeroing array
    pi=4.*ATAN(1.)
    G=6.67E-11
    msun=1.989E30
    au=1.496E11
    match=0
    match2=0
    match3=0
    match4=0
    itemp=0

    !Generate system parameters
    a = amin+(ran3(idum)*amax)
    ecc = ran3(idum)*0.99
    inc = ASIN(ran3(idum))
    phi = ran3(idum)*2.*pi
    mean1 = ran3(idum)*2.*pi

    !Period
    T  = SQRT((a*a*a)/(M(1)+M(2)))
    dM1 = (dtobs1/T)*2.*pi
    dM2 = (dtobs2/T)*2.*pi
    dM3 = (dtobs3/T)*2.*pi

    params=(/a,ecc,inc,phi,mean1/)

    !Projecting and testing first point
    CALL skyproject(a,ecc,inc,phi,mean1,rt,sky1)
    sep(1)=SQRT(sky1(1)*sky1(1)+sky1(2)*sky1(2))
    if (rt < obs(1)) goto 10
    if (ABS(obs(1)-sep(1))>prec(1)) goto 10

    rt=0.
    f=(/1.,-1./)

    !Checking second point in prograte and retrograde
    do i=1,2
        mean2 = mean1 + f(i)*dM1
        CALL skyproject(a,ecc,inc,phi,mean2,rt,sky2)
        if (rt < obs(2)) cycle
        !Velocity check
        v=G*(M(1)+M(2))*msun*((2./(rt*au))-(1./(a*au)))
        if (v < vobs) cycle
        !Separation
        sep(2)=SQRT(sky2(1)*sky2(1)+sky2(2)*sky2(2))
        if (ABS(obs(2)-sep(2)) > prec(2)) cycle
        !Angle
        sep(3)=ATAN2(sky2(2),sky2(1))-ATAN2(sky1(2),sky1(1))
        !if (ABS(ABS(obs(3))-ABS(sep(3))) > prec(3)) cycle
        if (obs(3)>(sep(3)-ABS(prec(3))).AND.obs(3)<(sep(3)+ABS(prec(3)))) then
            continue
        else
            cycle
        end if
        match2=1
        itemp=i
        exit
    end do

    if (match2==0) goto 10

    !Adding the third data point (same direction)----------------------
    mean3 = mean2 + f(itemp)*dM2
    CALL skyproject(a,ecc,inc,phi,mean3,rt,sky3)
    if (rt < obs(4)) goto 10
    !Velocity check
    v = G*(M(1)+M(2))*msun*((2./(rt*au))-(1./(a*au)))
    if (v < vobs) goto 10
    !Separation
    sep(4)=SQRT(sky3(1)*sky3(1)+sky3(2)*sky3(2))
    if (ABS(obs(4)-sep(4)) > prec(4)) goto 10
    !Angle
    sep(5)=ATAN2(sky3(2),sky3(1))-ATAN2(sky2(2),sky2(1))
    !if (ABS(ABS(obs(5))-ABS(sep(5))) > prec(5)) goto 10
    if (obs(5)>(sep(5)-ABS(prec(5))).AND.obs(5)<(sep(5)+ABS(prec(5)))) then
        continue
        match3=1
    else
        goto 10
    end if
    if (match3==0) goto 10

    !Adding the fourth data point (same direction)----------------------
    mean4 = mean3 + f(itemp)*dM3
    CALL skyproject(a,ecc,inc,phi,mean4,rt,sky4)
    if (rt < obs(6)) goto 10
    !Velocity check
    v = G*(M(1)+M(2))*msun*((2./(rt*au))-(1./(a*au)))
    if (v < vobs) goto 10
    !Separation
    sep(6)=SQRT(sky4(1)*sky4(1)+sky4(2)*sky4(2))
    if (ABS(obs(6)-sep(6)) > prec(6)) goto 10
    !Angle
    sep(7)=ATAN2(sky4(2),sky4(1))-ATAN2(sky3(2),sky3(1))
    !if (ABS(ABS(obs(5))-ABS(sep(5))) > prec(5)) goto 10
    if (obs(7)>(sep(7)-ABS(prec(7))).AND.obs(7)<(sep(7)+ABS(prec(7)))) then
        continue
        match4=1
    else
        goto 10
    end if
    if (match4==0) goto 10

    coords=(/sky1,sky2,sky3/)
    match=1

10  CONTINUE

CONTAINS
SUBROUTINE skyproject(a,ecc,inc,phi,mean,rt,sky)
    IMPLICIT NONE

    REAL, INTENT(in) :: a,ecc,inc,phi,mean
    REAL, INTENT(out):: rt,sky(1:2)
    REAL:: Ee,old,theta,r(1:3),x,y

    Ee=mean
    old=mean
    do
        old=Ee
        Ee=Ee-(Ee-ecc*SIN(Ee)-mean)/(1.-ecc*COS(Ee))
        if (old-Ee<1.e-6) exit
    end do

    theta=2.*ATAN(SQRT((1.+ecc)/(1.-ecc))*TAN(0.5*Ee))
    rt = a*(1.-ecc*ecc)/(1.+ecc*COS(theta))

    !Get x and y positions
    r(1:3)=0.
    r(1)=rt*COS(theta)
    r(2)=rt*SIN(theta)
    !Orientation around z axis
    x=r(1);y=r(2)
    r(1)=x*COS(phi)-y*SIN(phi)
    r(2)=x*SIN(phi)+y*COS(phi)
    !Modify inclination
    x=r(1);y=r(2)
    r(1)=x
    r(2)=y*COS(inc)
    r(3)=y*SIN(inc)
    !Positions on the sky
    sky(1)=r(1)
    sky(2)=r(3)
END SUBROUTINE skyproject

END SUBROUTINE fakesystem

!-----------------
SUBROUTINE cosine(a,b,thetai,thetaf,c)
    IMPLICIT NONE

    REAL, INTENT(in) :: a,b,thetai,thetaf
    REAL, INTENT(out):: c
    REAL:: pi,d2r,dtheta

    !Convert to radians
    pi=4.*ATAN(1.)
    dtheta=ABS(thetai-thetaf)
    d2r=(2.*pi)/360.
    c=SQRT(a*a + b*b - 2.*a*b*COS(dtheta*d2r))

END SUBROUTINE cosine

!======================================================
!--------------------RAN3------------------------------
REAL FUNCTION ran3(idum)
    INTEGER :: idum
    INTEGER :: MBIG,MSEED,MZ
    REAL :: FAC
    PARAMETER (MBIG=1000000000,MSEED=161803398,MZ=0,FAC=1./MBIG)
    INTEGER :: i,iff,ii,inext,inextp,k
    INTEGER :: mj,mk,ma(55)
    !SAVE iff,inext,inextp,ma
    COMMON /ranvar/ iff,inext,inextp,ma
    !$OMP THREADPRIVATE(/ranvar/)
    DATA iff /0/
    IF(idum.LT.0.OR.iff.EQ.0)THEN
      iff=1
      mj=MSEED-iabs(idum)
      mj=MOD(mj,MBIG)
      ma(55)=mj
      mk=1
      DO 11 i=1,54
        ii=MOD(21*i,55)
        ma(ii)=mk
        mk=mj-mk
        IF(mk.LT.MZ)mk=mk+MBIG
        mj=ma(ii)
    11      CONTINUE
      DO 13 k=1,4
        DO 12 i=1,55
          ma(i)=ma(i)-ma(1+MOD(i+30,55))
          IF(ma(i).LT.MZ)ma(i)=ma(i)+MBIG
    12        CONTINUE
    13      CONTINUE
      inext=0
      inextp=31
      idum=1
    ENDIF
    inext=inext+1
    IF(inext.EQ.56)inext=1
    inextp=inextp+1
    IF(inextp.EQ.56)inextp=1
    mj=ma(inext)-ma(inextp)
    IF(mj.LT.MZ)mj=mj+MBIG
    ma(inext)=mj
    ran3=mj*FAC
    RETURN
END FUNCTION ran3

max_stars_repo_name: rebeccahoughton/FOBOS

max_issues_repo_name: rebeccahoughton/FOBOS

max_forks_repo_name: rebeccahoughton/FOBOS


text:
!---------------------------------------------------------------------!
!                                PHAML                                !
!                                                                     !
! The Parallel Hierarchical Adaptive MultiLevel code for solving      !
! linear elliptic partial differential equations of the form          !
! (PUx)x + (QUy)y + RU = F on 2D polygonal domains with mixed         !
! boundary conditions, and eigenvalue problems where F is lambda*U.   !
!                                                                     !
! PHAML is public domain software.  It was produced as part of work   !
! done by the U.S. Government, and is not subject to copyright in     !
! the United States.                                                  !
!                                                                     !
!     William F. Mitchell                                             !
!     Applied and Computational Mathematics Division                  !
!     National Institute of Standards and Technology                  !
!     william.mitchell@nist.gov                                       !
!     http://math.nist.gov/phaml                                      !
!                                                                     !
!---------------------------------------------------------------------!

module sysdep

!----------------------------------------------------
! This module contains routines that are not standard Fortran 90 and may
! need to be changed for different compilers.
!----------------------------------------------------

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

implicit none
private
public my_system, my_flush

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

contains

!          ---------
subroutine my_system(cmd)
!          ---------

!----------------------------------------------------
! This routine executes a system command as if from the system command line.
! It is not supported by the standard, but I haven't had it fail on any
! compilers I have tried.
!----------------------------------------------------

!----------------------------------------------------
! Dummy arguments

character(len=*), intent(in) :: cmd
!----------------------------------------------------
! Begin executable code

call system(cmd)

end subroutine my_system

!          --------
subroutine my_flush(unit)
!          --------

!----------------------------------------------------
! This routine flushes the io buffer of the given unit.
!
! The flush statement is standard Fortran 2003.  According to
! ACM Fortran Forum August 2007, it is supported by Cray, gfortran,
! g95, IBM, and Intel, but not by NAG.  The article did not include
! Lahey, Absoft, PGI, Pathscale, or other compilers.
!
! Some compilers supply flush as a nonstandard intrinsic subroutine, so
! if the flush statement fails, you might try "call flush(unit)".
!
! If the compiler does not support flush, then just return without doing
! anything.
!
! Documentation for Lahey says it has flush as a subroutine.
!
! Documentation for Nag says it has a flush subroutine in a nonstandard
! intrinsic module, so might need
! use f90_unix_io, only:flush
!
! Documentation for Absoft makes no reference to flush.
!
! Documentation for PGI makes no reference to flush (except as an OpenMP
! directive), but the Pathscale documentation implies that PGI has it
! as a subroutine.
!
! Documentation for Pathscale says it has flush as a subroutine, but you might
! need -g77 or -pgi on the compiler command line.
!----------------------------------------------------

!----------------------------------------------------
! Dummy arguments

integer, intent(in) :: unit
!----------------------------------------------------
! Begin executable code

! Since more compilers support the nonstandard flush subroutine than the
! standard (2003) flush statement, I'll use the subroutine for now.

call flush(unit)

end subroutine my_flush

end module sysdep

max_stars_repo_name: qsnake/phaml

max_issues_repo_name: qsnake/phaml

max_forks_repo_name: qsnake/phaml


text:
c     program DRDBINOM
c>> 1998-05-12 DRDBINOM Krogh Ensure Pascal print output.
c>> 1995-12-15 DRDBINUM Krogh Initial Code.
c--D replaces "?": DR?BINOM, ?BINOM
c
c Test DBINOM for computing binomial coefficients.
c Checks a large value computed with factorials, one a little larger
c that is computed using log gamma, and then finds first that differs
c from that obtained using the Pascal triangle.
c
      integer K, N, NMAX
      parameter (NMAX = 120)
      external DBINOM
      double precision BC1, BC2, P(0:NMAX), TP, DBINOM
c            ( 150 )            ( 151 )
c      BC1 = (     )      BC2 = (     )
c            (  30 )            (  29 )
      parameter (BC1 = 32198785340494567031466236484400.D0)
      parameter (BC2 = 9880808670399701168713050486000.D0)
c
   10 format (/'  N    K',12X,'DBINOM',20X, A, 8X, ' Col. 2 - Col. 1')
   20 format (I4, I4,1P,E26.17,E26.17,1P,E16.7)
c
      print 10, ' True '
      TP =  DBINOM(150, 30)
      print 20, 150, 30, TP, BC1, BC1 - TP
      TP =  DBINOM(151, 29)
      print 20, 151, 29, TP, BC2, BC2 - TP
c
      print 10, 'Pascal'
      P(0) = 1.D0
      do 80 N = 0, NMAX
         do 40 K = 0, N
            TP = DBINOM(N, K)
            if (TP .ne. P(K)) go to 100
   40    continue
c                  Update the Pascal Triangle.
         P(N+1) = 1.D0
         do 60 K = N, 1, -1
            P(K) = P(K) + P(K-1)
   60    continue
   80 continue
  100 print 20, N, K, TP, P(K), P(K) - TP
      stop
      end

ext: f


text:
!> IMPACT
!! \author Rolf Henniger, Institute of Fluid Dynamics, ETH Zurich (henniger@ifd.mavt.ethz.ch)
!! \date Mai 2005 - Dec 2011

!> \brief module
module cmod_DivGrad2ndOOp

  use iso_c_binding

  implicit none

contains

  subroutine Op_getCDG_dir( &
      M,                    &
      N,                    &
      BL,                   &
      BU,                   &
      BCL,                  &
      BCU,                  &
      yu,                   &
      xp,                   &
      xu,                   &
      cDG ) bind( c, name='Op_getCDG_dir' )

    implicit none

    integer(c_int), intent(in)    :: M

    integer(c_int), intent(in)    :: N

    integer(c_int), intent(in)    :: BL
    integer(c_int), intent(in)    :: BU

    integer(c_int), intent(in)    :: BCL
    integer(c_int), intent(in)    :: BCU

    real(c_double), intent(in)    :: yu(0:M)

    real(c_double), intent(in)    :: xp(bl:N+bu)

    real(c_double), intent(in)    :: xu(bl:N+bu)

    real(c_double), intent(inout) :: cDG(-1:1,1:N)

    real(c_double)                :: cGp( 0:1,0:N)

    real(c_double)                :: cDu(-1:0,1:N)

    integer(c_int)                :: i

    cDG = 0.

    cGp = 0.
    cDu = 0.

    !---------------------------------------------------------------------------!
    ! Achtung: Falls nicht periodisch, wird hier am Rand ein Fehler gemacht,    !
    !          der nur über die Sonderbehandlung weiter unten korrigiert wird!  !
    !---------------------------------------------------------------------------!

    do i = 0, N
      cGp(0,i) = -1./( xp(i+1) - xp(i) )
      cGp(1,i) =  1./( xp(i+1) - xp(i) )
    end do

    do i = 1, N
      cDu(-1,i) = -1./( xu(i) - xu(i-1) )
      cDu( 0,i) =  1./( xu(i) - xu(i-1) )
    end do

    if( BCL > 0 ) then
      cGp(:,0) = 0.
      cDu(:,1) = 0.
      cDu(0,1) = 1./( yu(1) - yu(0) ) ! TEST!!! Der Schoenheit halber, s.u.
      !cDu(0,1) = 1./( xp(2) - xp(1) ) ! TEST!!! Der Schoenheit halber, s.u.
      !cDu(0,1) =  1./(xu(1) - xu (0)) !! shoudl be the same
    end if

    if( BCU > 0 ) then
      cGp(:,N) =  0.
      cDu(:,N) =  0.
      cDu(-1,N) = -1./( yu(M) - yu(M-1) )
      !cDu(-1,N) = -1./( xp(N) - xp(N-1) ) ! consistency bad, good convergence,
      ! consistency prefered
      !cDu(-1,N) = -1./( xu(N) - xu(N-1) ) ! TEST!!! Das geht in die Hose ...
      !cDu(-1,N) = -1./( yu(M+1) - yu(M) )
    end if
    !-------------------------------------------------------------------------------------------
    do i = 1, N
      cDG(-1,i) = cDu(-1,i)*cGp(0,i-1)
      cDG( 0,i) = cDu(-1,i)*cGp(1,i-1) + cDu(0,i)*cGp(0,i)
      cDG( 1,i) =                        cDu(0,i)*cGp(1,i)
    end do

    ! Faktor 2 kommt von der Approximation DH-¹D ~= DI-¹G, wobei I hier das
    ! Interpolationspolynom am Rand darstellt
    if( BCL > 0 ) then
      cDG( :,1) = 2.*cDG(:,1)
      cDG(-1,1) = 0.
    end if
    if( BCU > 0 ) then
      cDG( :,N) = 2.*cDG(:,N)
      cDG( 1,N) = 0.
    end if

    if( BCL == -2) then
      cDG( 1,1) = cDG( 1, 1 ) + cDG( -1, 1 )
      cDG(-1,1) = 0.
    end if

    if( BCU == -2 ) then
      cDG(-1,N) = cDG(-1, N ) + cDG( 1, N )
      cDG( 1,N) = 0.
    end if

  end subroutine Op_getCDG_dir

end module cmod_DivGrad2ndOOp

max_stars_repo_name: huppd/PINTimpact

max_issues_repo_name: huppd/PINTimpact

max_forks_repo_name: huppd/PINTimpact


text:
      double complex function ttbgggppp(i1,i2,i3,i4,i5,i6,i7)
      implicit none
      include 'constants.f'
      include 'sprods_com.f'
      include 'zprods_com.f'
      include 'masses.f'
      integer i1,i2,i3,i4,i5,i6,i7
      double precision s129,s1345,s6789,mtsq
      s129=s(i1,i2)+s(i1,i3)+s(i2,i3)
      s1345=s(i1,i4)+s(i1,i5)
      s6789=s(i3,i6)+s(i3,i7)
      mtsq=mt**2
      ttbgggppp =  + mtsq*s129**(-1) * ( 1/(za(i1,i2))/(za(i1,i5))*za(
     &    i5,i7)*zb(i2,i3)*zb(i4,i3) + 1/(za(i1,i2))/(za(i1,i5))/(za(i5
     &    ,i3))*za(i2,i5)*za(i5,i7)*zb(i2,i3)*zb(i2,i4) + 1/(za(i1,i2))
     &    /(za(i2,i5))*za(i5,i7)*zb(i1,i3)*zb(i4,i3) + 1/(za(i1,i2))/(
     &    za(i2,i5))/(za(i5,i3))*za(i1,i5)*za(i5,i7)*zb(i1,i3)*zb(i1,i4
     &    ) + 1/(za(i1,i2))/(za(i5,i3))*za(i5,i7)*zb(i1,i3)*zb(i2,i4)
     &     + 1/(za(i1,i2))/(za(i5,i3))*za(i5,i7)*zb(i1,i4)*zb(i2,i3) +
     &    1/(za(i1,i5))/(za(i2,i3))*za(i5,i7)*zb(i1,i2)*zb(i4,i3) - 1/(
     &    za(i1,i5))/(za(i2,i3))*za(i5,i7)*zb(i1,i3)*zb(i2,i4) - 1/(za(
     &    i1,i5))/(za(i2,i3))/(za(i2,i5))*za(i5,i3)*za(i5,i7)*zb(i1,i3)
     &    *zb(i4,i3) + 1/(za(i1,i5))/(za(i2,i3))/(za(i5,i3))*za(i2,i5)*
     &    za(i5,i7)*zb(i1,i2)*zb(i2,i4) - 1/(za(i2,i3))/(za(i2,i5))*za(
     &    i5,i7)*zb(i1,i3)*zb(i1,i4) + 1/(za(i2,i3))/(za(i5,i3))*za(i5,
     &    i7)*zb(i1,i2)*zb(i1,i4) )
      ttbgggppp = ttbgggppp + mtsq*s1345**(-1)*s6789**(-1) * ( 1/(za(i1
     &    ,i5))/(za(i2,i5))/(za(i5,i3))*za(i4,i5)*za(i5,i6)*za(i5,i7)*
     &    zb(i1,i4)*zb(i2,i4)*zb(i6,i3) + 1/(za(i1,i5))/(za(i2,i5))/(
     &    za(i5,i3))*za(i4,i5)*za(i5,i7)**2*zb(i1,i4)*zb(i2,i4)*zb(i7,
     &    i3) - 1/(za(i2,i5))/(za(i5,i3))*za(i5,i6)*za(i5,i7)*zb(i1,i2)
     &    *zb(i1,i4)*zb(i6,i3) - 1/(za(i2,i5))/(za(i5,i3))*za(i5,i7)**2
     &    *zb(i1,i2)*zb(i1,i4)*zb(i7,i3) )
      ttbgggppp = ttbgggppp + mtsq*s1345**(-1) * (  - 1/(za(i1,i5))/(
     &    za(i2,i3))/(za(i2,i5))*za(i4,i5)*za(i5,i7)*zb(i1,i4)*zb(i4,i3
     &    ) - 1/(za(i1,i5))/(za(i2,i3))/(za(i5,i3))*za(i4,i5)*za(i5,i7)
     &    *zb(i1,i4)*zb(i2,i4) - 1/(za(i2,i3))/(za(i2,i5))*za(i5,i7)*
     &    zb(i1,i3)*zb(i1,i4) + 1/(za(i2,i3))/(za(i5,i3))*za(i5,i7)*zb(
     &    i1,i2)*zb(i1,i4) )
      ttbgggppp = ttbgggppp + mtsq*s6789**(-1) * ( 1/(za(i1,i2))/(za(i1
     &    ,i5))/(za(i5,i3))*za(i5,i6)*za(i5,i7)*zb(i2,i4)*zb(i6,i3) +
     &    1/(za(i1,i2))/(za(i1,i5))/(za(i5,i3))*za(i5,i7)**2*zb(i2,i4)*
     &    zb(i7,i3) + 1/(za(i1,i2))/(za(i2,i5))/(za(i5,i3))*za(i5,i6)*
     &    za(i5,i7)*zb(i1,i4)*zb(i6,i3) + 1/(za(i1,i2))/(za(i2,i5))/(
     &    za(i5,i3))*za(i5,i7)**2*zb(i1,i4)*zb(i7,i3) )

      return
      end

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


program test_pslcal
  use simsphere_mod, only: thv, w2g, fs, zg, zp, b1, b2, rmin, vegheight, sc,&
                           xlai, raf, rmin, psice, vfl, qstf, rst, ps1, sol, &
                           steady, rcut, rscrit, beta, psisup, psig, psiwc,  &
                           psim, psie, mintemp, maxtemp, wpsi, rlpsi, thmax, &
                           rks, cosbyb, psis, rccap, zst, psist, zstini,     &
                           volrel, volini, rzcap, psix, frhgt, rlelf, frzp,  &
                           rkocap, eq
  use mod_testing, only: assert, initialize_tests, report_tests
  implicit none

  logical, dimension(:), allocatable :: tests
  logical :: test_failed
  integer :: n, ntests

  real, parameter :: arg1 = 1.0
  real, parameter :: arg2 = 5.0
  integer :: arg3

  ! Expected values
  real, parameter :: wpsi_unsteady_exp = 23.7398777
  real, parameter :: wpsi_steady_exp = 29.0206356
  real, parameter :: rlpsi_unsteady_exp = 24.9543571
  real, parameter :: rlpsi_steady_exp = 28.9508419

  n = 1
  ntests = 4
  call initialize_tests(tests,ntests)

  call pslcal_init
  arg3 = 1
  steady = 'N'
  call pslcal(arg1,arg2,arg3)
  tests(n) = assert(eq(wpsi,wpsi_unsteady_exp), 'wpsi unsteady')
  n = n + 1
  tests(n) = assert(eq(rlpsi,rlpsi_unsteady_exp), 'rlpsi unsteady')
  n = n + 1

  call pslcal_init
  arg3 = 1
  steady = 'Y'
  call pslcal(arg1,arg2,arg3)
  !write(6,*) 'pslcal steady: wpsi,rlpsi: ',wpsi,rlpsi
  tests(n) = assert(eq(wpsi,wpsi_steady_exp), 'wpsi steady')
  n = n + 1
  tests(n) = assert(eq(rlpsi,rlpsi_steady_exp), 'rlpsi steady')
  n = n + 1

  test_failed = .false.
  call report_tests(tests,test_failed)
  if (test_failed) stop 1

contains
  subroutine pslcal_init
    thv = 3.0
    thmax = 5.0
    w2g = 0.25
    fs = 1.0
    zg = 1.0
    zp = 1.7
    frzp = 1.0
    b1 = 5.0
    b2 = 2.0
    cosbyb = 3.0
    rmin = 1.0
    vegheight = 2.0
    xlai = 1.0
    raf = 3.0
    rmin = 0.1
    psice = 20.0
    vfl = 1.0
    qstf = 2.5
    rst = 1.3
    rks = 3.0
    ps1 = 2.7
    rcut = 1.0
    rscrit = 2.5
    beta = 0.3
    psisup = 0.5
    psig = 1.3
    psiwc = 2.7
    psist = 1.5
    psim = 1.3
    psie = 1.5
    psis = 1.0
    mintemp = 0.0
    maxtemp = 100.0
    wpsi = 20.0
    rlpsi = 14.0
    rccap = 1.0
    sc = 1.0
    sol = 1.0
    volrel = 1.0
    volini = 1.0
    psix = 1.0
    frhgt = 1.0
    rlelf = 1.0
    rzcap = 1.0
    rkocap = 1.0
    zstini = 0.001
    zst = 0.001
    return
  end subroutine pslcal_init

end program test_pslcal


text:
module broken_mod
    integer
end

program p1
!$claw ignore
    use broken_mod
    integer :: x = 1, y = 2, z
    z = x + y
!$claw end ignore
end program p1

max_stars_repo_name: clementval/claw-compiler

max_issues_repo_name: skosukhin/claw-compiler

max_forks_repo_name: skosukhin/claw-compiler


text: Warren R. Pugh served on the City Council from 20 April 1918 to 20 April 1940. He was born in Davis 15 October 1877 (WWII draft registration and the CA Death Index summary) to David H. Pugh, a farm laborer, and his wife Alida Hoffman Pugh. He was a barber with his own shop (1910, 1920 census) and farmer (1930 census). He died in Alameda County, CA, 30 July 1945.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


             program demo_pdec
             use M_verify, only : pdec
             call pdec(' ABCDEFG abcdefg    ')
             end program demo_pdec


text:
      SUBROUTINE slF5HZ (R5,D5,EPOCH,RH,DH)
*+
*     - - - - - -
*      F 5 H Z
*     - - - - - -
*
*  Transform an FK5 (J2000) star position into the frame of the
*  Hipparcos catalogue, assuming zero Hipparcos proper motion.
*
*  (double precision)
*
*  This routine converts a star position from the FK5 system to
*  the Hipparcos system, in such a way that the Hipparcos proper
*  motion is zero.  Because such a star has, in general, a non-zero
*  proper motion in the FK5 system, the routine requires the epoch
*  at which the position in the FK5 system was determined.
*
*  Given:
*     R5        d      FK5 RA (radians), equinox J2000, epoch EPOCH
*     D5        d      FK5 Dec (radians), equinox J2000, epoch EPOCH
*     EPOCH     d      Julian epoch (TDB)
*
*  Returned (all Hipparcos):
*     RH        d      RA (radians)
*     DH        d      Dec (radians)
*
*  Called:  slDS2C, slDAVM, slDIMV, slDMXV, slDC2S,
*           slDA2P
*
*  Notes:
*
*  1)  The FK5 to Hipparcos transformation consists of a pure
*      rotation and spin;  zonal errors in the FK5 catalogue are
*      not taken into account.
*
*  2)  The published orientation and spin components are interpreted
*      as "axial vectors".  An axial vector points at the pole of the
*      rotation and its length is the amount of rotation in radians.
*
*  3)  See also slFK5H, slHFK5, slHF5Z.
*
*  Reference:
*
*     M.Feissel & F.Mignard, Astron. Astrophys. 331, L33-L36 (1998).
*
*  P.T.Wallace   Starlink   22 June 1999
*
*  Copyright (C) 1999 Rutherford Appleton Laboratory
*
*  License:
*    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; either version 2 of the License, or
*    (at your option) any later version.
*
*    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 (see SLA_CONDITIONS); if not, write to the
*    Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
*    Boston, MA  02110-1301  USA
*
*  Copyright (C) 1995 Association of Universities for Research in Astronomy Inc.
*-

      IMPLICIT NONE

      DOUBLE PRECISION R5,D5,EPOCH,RH,DH

      DOUBLE PRECISION AS2R
      PARAMETER (AS2R=0.484813681109535994D-5)

*  FK5 to Hipparcos orientation and spin (radians, radians/year)
      DOUBLE PRECISION EPX,EPY,EPZ
      DOUBLE PRECISION OMX,OMY,OMZ

      PARAMETER ( EPX = -19.9D-3 * AS2R,
     :            EPY =  -9.1D-3 * AS2R,
     :            EPZ = +22.9D-3 * AS2R )

      PARAMETER ( OMX = -0.30D-3 * AS2R,
     :            OMY = +0.60D-3 * AS2R,
     :            OMZ = +0.70D-3 * AS2R )

      DOUBLE PRECISION P5E(3),ORTN(3),R5H(3,3),T,VST(3),RST(3,3),
     :                 P5(3),PH(3),W

      DOUBLE PRECISION slDA2P

*  FK5 barycentric position vector.
      CALL slDS2C(R5,D5,P5E)

*  FK5 to Hipparcos orientation matrix.
      ORTN(1) = EPX
      ORTN(2) = EPY
      ORTN(3) = EPZ
      CALL slDAVM(ORTN,R5H)

*  Time interval from epoch to J2000.
      T = 2000D0-EPOCH

*  Axial vector:  accumulated Hipparcos wrt FK5 spin over that interval.
      VST(1) = OMX*T
      VST(2) = OMY*T
      VST(3) = OMZ*T

*  Express the accumulated spin as a rotation matrix.
      CALL slDAVM(VST,RST)

*  Derotate the vector's FK5 axes back to epoch.
      CALL slDIMV(RST,P5E,P5)

*  Rotate the vector into the Hipparcos frame.
      CALL slDMXV(R5H,P5,PH)

*  Hipparcos vector to spherical.
      CALL slDC2S(PH,W,DH)
      RH = slDA2P(W)

      END

max_stars_repo_name: ysBach/irafdocgen

max_issues_repo_name: kirxkirx/iraf

max_forks_repo_name: kirxkirx/iraf


text:
      SUBROUTINE AUSPECF(N,LAGH1,CXX1,P1,P2,Q)
C

C
C     PROGRAM 5.2.1   POWER SPECTRUM
C-----------------------------------------------------------------------
C     ** DESIGNED BY H. AKAIKE, THE INSTITUTE OF STATISTICAL MATHEMATICS
C     ** PROGRAMMED BY E. ARAHATA, THE INSTITUTE OF STATISTICAL MATHEMAT
C        TOKYO
C     ** DATE OF THE LATEST REVISION: MARCH 25, 1977
C     ** THIS PROGRAM WAS ORIGINALLY PUBLISHED IN
C        "DAINAMIKKU SISTEMU NO TOKEI-TEKI KAISEKI TO SEIGYO (STATISTICA
C        ANALYSIS AND CONTROL OF DYNAMIC SYSTEMS)" BY H. AKAIKE AND
C        T. NAKAGAWA, SAIENSU-SHA, TOKYO, 1972 (IN JAPANESE)
C-----------------------------------------------------------------------
C     THIS PROGRAM COMPUTES POWER SPECTRUM ESTIMATES FOR TWO
C     TRIGONOMETRIC WINDOWS OF BLACKMAN-TUKEY TYPE BY GOERTZEL METHOD.
C     ONLY ONE CARD OF LAGH(MAXIMUM LAG OF COVARIANCES TO BE USED FOR
C     POWER SPECTRUM COMPUTATION) SHOULD BE ADDED ON TOP OF THE OUTPUT
C     OF PROGRAM 5.1.1 AUTCOR TO FORM INPUT TO THIS PROGRAM.
C     OUTPUTS ARE ESTIMATES P1(I),P2(I) FOR FREQUENCIES I/(2LAGH*DELTAT)
C     AND THE TEST STATISTICS Q(I) FOR THE DIFFERENCES BETWEEN P1(I) AND
C     P2(I).   Q(I) GREATER THAN 1 MEANS SIGNIFICANT DIFFERENCE.
C
cc      !DEC$ ATTRIBUTES DLLEXPORT :: AUSPECF
c      USE DFLIB
C
cxx      IMPLICIT REAL*8(A-H,O-Z)
c      DIMENSION CXX(1001),FC(1001),P1(1001),P2(1001),Q(1001)
c      DIMENSION A1(10),A2(10)
cxx      DIMENSION CXX1(LAGH1),CXX(LAGH1)
cxx      DIMENSION FC(LAGH1),P1(LAGH1),P2(LAGH1),Q(LAGH1)
cxx      DIMENSION A1(2),A2(3)
      INTEGER :: N, LAGH1
      REAL(8) :: CXX1(LAGH1), P1(LAGH1), P2(LAGH1), Q(LAGH1)
      REAL(8) :: CXX(LAGH1), FC(LAGH1), A1(2), A2(3)
C
C     INPUT / OUTPUT DATA FILE OPEN
c	CHARACTER(100) DFNAM
c	DFNAM='auspec.out'
c	CALL SETWND
c	CALL FLOPN3(DFNAM,NFL)
c	IF (NFL.EQ.0) GO TO 999
C
C     WINDOW W1 DEFINITION
      MLA1=2
      A1(1)=0.5D-00
      A1(2)=0.25D-00
C     WINDOW W2 DEFINITION
      MLA2=3
      A2(1)=0.625D-00
      A2(2)=0.25D-00
      A2(3)=-0.0625D-00
C
C     LAGH SPECIFICATION
c      READ(5,1) LAGH
c      LAGH1=LAGH+1
C     READING THE OUTPUTS OF PROGRAM 5.1.1 AUTCOR
c      READ(5,1) N,LAGH0
c      LAGH3=LAGH0+1
C     INITIAL CONDITION PRINT OUT
c      WRITE(6,60)
c      WRITE(6,61)
c      WRITE(6,62) N,LAGH
c      WRITE(6,63)
c      CALL PRCOL1(A1,1,MLA1,1)
c      WRITE(6,64)
c      CALL PRCOL1(A2,1,MLA2,1)
C     AUTO COVARIANCE INPUT
c      READ(5,2) (CXX(I),I=1,LAGH3)
C     AUTO COVARIANCE PRINT OUT
c      WRITE(6,172)
c      WRITE(6,162)
c      CALL PRCOL1(CXX,1,LAGH1,1)
c      DO 10 I=2,LAGH
c   10 CXX(I)=CXX(I)+CXX(I)
      CXX(1)=CXX1(1)
      DO 10 I=2,LAGH1-1
cxx   10 CXX(I)=CXX1(I)+CXX1(I)
      CXX(I)=CXX1(I)+CXX1(I)
   10 CONTINUE
      CXX(LAGH1)=CXX1(LAGH1)
C     F-COS TRANSFORMATION
      CALL FGERCO(CXX,LAGH1,FC,LAGH1)
C     SPECTRUM SMOOTHING BY WINDOW W1
      CALL AUSP(FC,P1,LAGH1,A1,MLA1)
C     SPECTRUM SMOOTHING BY WINDOW W2
      CALL AUSP(FC,P2,LAGH1,A2,MLA2)
C     TEST STATISTICS COMPUTATION
      CALL SIGNIF(P1,P2,Q,LAGH1,N)
C     AUTO SPECTRUM AND TEST STATISTICS PRINT OUT
c      WRITE(6,66) N,LAGH
c      WRITE(6,67)
c      CALL PRCOL3(P1,P2,Q,1,LAGH1,1)
C     AUTO SPECTRUM SMOOTHED BY WINDOW W1 PUNCH OUT
c      WRITE(7,1) N,LAGH
c      WRITE(7,2) (P1(I),I=1,LAGH1)
c	CALL FLCLS3(NFL)
c    1 FORMAT(10I5)
c    2 FORMAT(4D20.10)
c   60 FORMAT(1H ,30HPROGRAM 5.2.1   POWER SPECTRUM)
c   62 FORMAT(1H ,2HN=,I5,5X,5HLAGH=,I5)
c   61 FORMAT(1H ,17HINITIAL CONDITION)
c   63 FORMAT(1H ,12X,9HWINDOW W1/1H ,4X,1HI,11X,5HA1(I))
c   64 FORMAT(1H ,12X,9HWINDOW W2/1H ,4X,1HI,11X,5HA2(I))
c   66 FORMAT(//1H ,14HPOWER SPECTRUM,5X,2HN=,I5,5X,5HLAGH=,I5)
c   67 FORMAT(1H ,4X,1HI,8X,8HPOWER W1,6X,8HPOWER W2,2X,12HSIGNIFICANCE)
c  162 FORMAT(1H ,4X,1HI,6X,6HCXX(I))
c  172 FORMAT(1H ,5X,15HAUTO COVARIANCE)
      RETURN
      END SUBROUTINE


text: Address(Levee Road) is a road in South Campus South and West Campus that runs along Putah Creek. It is open to public vehicles between Hopkins Road and Brooks Road where the Putah Creek Picnic Area is located. The rest is generally open to the public for recreation and research access to the Putah Creek Riparian Reserve.


text:
      SUBROUTINE ZLAR1V( N, B1, BN, LAMBDA, D, L, LD, LLD,
     $           PIVMIN, GAPTOL, Z, WANTNC, NEGCNT, ZTZ, MINGMA,
     $           R, ISUPPZ, NRMINV, RESID, RQCORR, WORK )
*
*  -- LAPACK auxiliary routine (version 3.2) --
*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
*     November 2006
*
*     .. Scalar Arguments ..
      LOGICAL            WANTNC
      INTEGER   B1, BN, N, NEGCNT, R
      DOUBLE PRECISION   GAPTOL, LAMBDA, MINGMA, NRMINV, PIVMIN, RESID,
     $                   RQCORR, ZTZ
*     ..
*     .. Array Arguments ..
      INTEGER            ISUPPZ( * )
      DOUBLE PRECISION   D( * ), L( * ), LD( * ), LLD( * ),
     $                  WORK( * )
      COMPLEX*16       Z( * )
*     ..
*
*  Purpose
*  =======
*
*  ZLAR1V computes the (scaled) r-th column of the inverse of
*  the sumbmatrix in rows B1 through BN of the tridiagonal matrix
*  L D L^T - sigma I. When sigma is close to an eigenvalue, the
*  computed vector is an accurate eigenvector. Usually, r corresponds
*  to the index where the eigenvector is largest in magnitude.
*  The following steps accomplish this computation :
*  (a) Stationary qd transform,  L D L^T - sigma I = L(+) D(+) L(+)^T,
*  (b) Progressive qd transform, L D L^T - sigma I = U(-) D(-) U(-)^T,
*  (c) Computation of the diagonal elements of the inverse of
*      L D L^T - sigma I by combining the above transforms, and choosing
*      r as the index where the diagonal of the inverse is (one of the)
*      largest in magnitude.
*  (d) Computation of the (scaled) r-th column of the inverse using the
*      twisted factorization obtained by combining the top part of the
*      the stationary and the bottom part of the progressive transform.
*
*  Arguments
*  =========
*
*  N        (input) INTEGER
*           The order of the matrix L D L^T.
*
*  B1       (input) INTEGER
*           First index of the submatrix of L D L^T.
*
*  BN       (input) INTEGER
*           Last index of the submatrix of L D L^T.
*
*  LAMBDA    (input) DOUBLE PRECISION
*           The shift. In order to compute an accurate eigenvector,
*           LAMBDA should be a good approximation to an eigenvalue
*           of L D L^T.
*
*  L        (input) DOUBLE PRECISION array, dimension (N-1)
*           The (n-1) subdiagonal elements of the unit bidiagonal matrix
*           L, in elements 1 to N-1.
*
*  D        (input) DOUBLE PRECISION array, dimension (N)
*           The n diagonal elements of the diagonal matrix D.
*
*  LD       (input) DOUBLE PRECISION array, dimension (N-1)
*           The n-1 elements L(i)*D(i).
*
*  LLD      (input) DOUBLE PRECISION array, dimension (N-1)
*           The n-1 elements L(i)*L(i)*D(i).
*
*  PIVMIN   (input) DOUBLE PRECISION
*           The minimum pivot in the Sturm sequence.
*
*  GAPTOL   (input) DOUBLE PRECISION
*           Tolerance that indicates when eigenvector entries are negligible
*           w.r.t. their contribution to the residual.
*
*  Z        (input/output) COMPLEX*16       array, dimension (N)
*           On input, all entries of Z must be set to 0.
*           On output, Z contains the (scaled) r-th column of the
*           inverse. The scaling is such that Z(R) equals 1.
*
*  WANTNC   (input) LOGICAL
*           Specifies whether NEGCNT has to be computed.
*
*  NEGCNT   (output) INTEGER
*           If WANTNC is .TRUE. then NEGCNT = the number of pivots < pivmin
*           in the  matrix factorization L D L^T, and NEGCNT = -1 otherwise.
*
*  ZTZ      (output) DOUBLE PRECISION
*           The square of the 2-norm of Z.
*
*  MINGMA   (output) DOUBLE PRECISION
*           The reciprocal of the largest (in magnitude) diagonal
*           element of the inverse of L D L^T - sigma I.
*
*  R        (input/output) INTEGER
*           The twist index for the twisted factorization used to
*           compute Z.
*           On input, 0 <= R <= N. If R is input as 0, R is set to
*           the index where (L D L^T - sigma I)^{-1} is largest
*           in magnitude. If 1 <= R <= N, R is unchanged.
*           On output, R contains the twist index used to compute Z.
*           Ideally, R designates the position of the maximum entry in the
*           eigenvector.
*
*  ISUPPZ   (output) INTEGER array, dimension (2)
*           The support of the vector in Z, i.e., the vector Z is
*           nonzero only in elements ISUPPZ(1) through ISUPPZ( 2 ).
*
*  NRMINV   (output) DOUBLE PRECISION
*           NRMINV = 1/SQRT( ZTZ )
*
*  RESID    (output) DOUBLE PRECISION
*           The residual of the FP vector.
*           RESID = ABS( MINGMA )/SQRT( ZTZ )
*
*  RQCORR   (output) DOUBLE PRECISION
*           The Rayleigh Quotient correction to LAMBDA.
*           RQCORR = MINGMA*TMP
*
*  WORK     (workspace) DOUBLE PRECISION array, dimension (4*N)
*
*  Further Details
*  ===============
*
*  Based on contributions by
*     Beresford Parlett, University of California, Berkeley, USA
*     Jim Demmel, University of California, Berkeley, USA
*     Inderjit Dhillon, University of Texas, Austin, USA
*     Osni Marques, LBNL/NERSC, USA
*     Christof Voemel, University of California, Berkeley, USA
*
*  =====================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   ZERO, ONE
      PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
      COMPLEX*16         CONE
      PARAMETER          ( CONE = ( 1.0D0, 0.0D0 ) )

*     ..
*     .. Local Scalars ..
      LOGICAL            SAWNAN1, SAWNAN2
      INTEGER            I, INDLPL, INDP, INDS, INDUMN, NEG1, NEG2, R1,
     $                   R2
      DOUBLE PRECISION   DMINUS, DPLUS, EPS, S, TMP
*     ..
*     .. External Functions ..
      LOGICAL DISNAN
      DOUBLE PRECISION   DLAMCH
      EXTERNAL           DISNAN, DLAMCH
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          ABS, DBLE
*     ..
*     .. Executable Statements ..
*
      EPS = DLAMCH( 'Precision' )

      IF( R.EQ.0 ) THEN
         R1 = B1
         R2 = BN
      ELSE
         R1 = R
         R2 = R
      END IF

*     Storage for LPLUS
      INDLPL = 0
*     Storage for UMINUS
      INDUMN = N
      INDS = 2*N + 1
      INDP = 3*N + 1

      IF( B1.EQ.1 ) THEN
         WORK( INDS ) = ZERO
      ELSE
         WORK( INDS+B1-1 ) = LLD( B1-1 )
      END IF

*
*     Compute the stationary transform (using the differential form)
*     until the index R2.
*
      SAWNAN1 = .FALSE.
      NEG1 = 0
      S = WORK( INDS+B1-1 ) - LAMBDA
      DO 50 I = B1, R1 - 1
         DPLUS = D( I ) + S
         WORK( INDLPL+I ) = LD( I ) / DPLUS
         IF(DPLUS.LT.ZERO) NEG1 = NEG1 + 1
         WORK( INDS+I ) = S*WORK( INDLPL+I )*L( I )
         S = WORK( INDS+I ) - LAMBDA
 50   CONTINUE
      SAWNAN1 = DISNAN( S )
      IF( SAWNAN1 ) GOTO 60
      DO 51 I = R1, R2 - 1
         DPLUS = D( I ) + S
         WORK( INDLPL+I ) = LD( I ) / DPLUS
         WORK( INDS+I ) = S*WORK( INDLPL+I )*L( I )
         S = WORK( INDS+I ) - LAMBDA
 51   CONTINUE
      SAWNAN1 = DISNAN( S )
*
 60   CONTINUE
      IF( SAWNAN1 ) THEN
*        Runs a slower version of the above loop if a NaN is detected
         NEG1 = 0
         S = WORK( INDS+B1-1 ) - LAMBDA
         DO 70 I = B1, R1 - 1
            DPLUS = D( I ) + S
            IF(ABS(DPLUS).LT.PIVMIN) DPLUS = -PIVMIN
            WORK( INDLPL+I ) = LD( I ) / DPLUS
            IF(DPLUS.LT.ZERO) NEG1 = NEG1 + 1
            WORK( INDS+I ) = S*WORK( INDLPL+I )*L( I )
            IF( WORK( INDLPL+I ).EQ.ZERO )
     $                      WORK( INDS+I ) = LLD( I )
            S = WORK( INDS+I ) - LAMBDA
 70      CONTINUE
         DO 71 I = R1, R2 - 1
            DPLUS = D( I ) + S
            IF(ABS(DPLUS).LT.PIVMIN) DPLUS = -PIVMIN
            WORK( INDLPL+I ) = LD( I ) / DPLUS
            WORK( INDS+I ) = S*WORK( INDLPL+I )*L( I )
            IF( WORK( INDLPL+I ).EQ.ZERO )
     $                      WORK( INDS+I ) = LLD( I )
            S = WORK( INDS+I ) - LAMBDA
 71      CONTINUE
      END IF
*
*     Compute the progressive transform (using the differential form)
*     until the index R1
*
      SAWNAN2 = .FALSE.
      NEG2 = 0
      WORK( INDP+BN-1 ) = D( BN ) - LAMBDA
      DO 80 I = BN - 1, R1, -1
         DMINUS = LLD( I ) + WORK( INDP+I )
         TMP = D( I ) / DMINUS
         IF(DMINUS.LT.ZERO) NEG2 = NEG2 + 1
         WORK( INDUMN+I ) = L( I )*TMP
         WORK( INDP+I-1 ) = WORK( INDP+I )*TMP - LAMBDA
 80   CONTINUE
      TMP = WORK( INDP+R1-1 )
      SAWNAN2 = DISNAN( TMP )

      IF( SAWNAN2 ) THEN
*        Runs a slower version of the above loop if a NaN is detected
         NEG2 = 0
         DO 100 I = BN-1, R1, -1
            DMINUS = LLD( I ) + WORK( INDP+I )
            IF(ABS(DMINUS).LT.PIVMIN) DMINUS = -PIVMIN
            TMP = D( I ) / DMINUS
            IF(DMINUS.LT.ZERO) NEG2 = NEG2 + 1
            WORK( INDUMN+I ) = L( I )*TMP
            WORK( INDP+I-1 ) = WORK( INDP+I )*TMP - LAMBDA
            IF( TMP.EQ.ZERO )
     $          WORK( INDP+I-1 ) = D( I ) - LAMBDA
 100     CONTINUE
      END IF
*
*     Find the index (from R1 to R2) of the largest (in magnitude)
*     diagonal element of the inverse
*
      MINGMA = WORK( INDS+R1-1 ) + WORK( INDP+R1-1 )
      IF( MINGMA.LT.ZERO ) NEG1 = NEG1 + 1
      IF( WANTNC ) THEN
         NEGCNT = NEG1 + NEG2
      ELSE
         NEGCNT = -1
      ENDIF
      IF( ABS(MINGMA).EQ.ZERO )
     $   MINGMA = EPS*WORK( INDS+R1-1 )
      R = R1
      DO 110 I = R1, R2 - 1
         TMP = WORK( INDS+I ) + WORK( INDP+I )
         IF( TMP.EQ.ZERO )
     $      TMP = EPS*WORK( INDS+I )
         IF( ABS( TMP ).LE.ABS( MINGMA ) ) THEN
            MINGMA = TMP
            R = I + 1
         END IF
 110  CONTINUE
*
*     Compute the FP vector: solve N^T v = e_r
*
      ISUPPZ( 1 ) = B1
      ISUPPZ( 2 ) = BN
      Z( R ) = CONE
      ZTZ = ONE
*
*     Compute the FP vector upwards from R
*
      IF( .NOT.SAWNAN1 .AND. .NOT.SAWNAN2 ) THEN
         DO 210 I = R-1, B1, -1
            Z( I ) = -( WORK( INDLPL+I )*Z( I+1 ) )
            IF( (ABS(Z(I))+ABS(Z(I+1)))* ABS(LD(I)).LT.GAPTOL )
     $           THEN
               Z( I ) = ZERO
               ISUPPZ( 1 ) = I + 1
               GOTO 220
            ENDIF
            ZTZ = ZTZ + DBLE( Z( I )*Z( I ) )
 210     CONTINUE
 220     CONTINUE
      ELSE
*        Run slower loop if NaN occurred.
         DO 230 I = R - 1, B1, -1
            IF( Z( I+1 ).EQ.ZERO ) THEN
               Z( I ) = -( LD( I+1 ) / LD( I ) )*Z( I+2 )
            ELSE
               Z( I ) = -( WORK( INDLPL+I )*Z( I+1 ) )
            END IF
            IF( (ABS(Z(I))+ABS(Z(I+1)))* ABS(LD(I)).LT.GAPTOL )
     $           THEN
               Z( I ) = ZERO
               ISUPPZ( 1 ) = I + 1
               GO TO 240
            END IF
            ZTZ = ZTZ + DBLE( Z( I )*Z( I ) )
 230     CONTINUE
 240     CONTINUE
      ENDIF

*     Compute the FP vector downwards from R in blocks of size BLKSIZ
      IF( .NOT.SAWNAN1 .AND. .NOT.SAWNAN2 ) THEN
         DO 250 I = R, BN-1
            Z( I+1 ) = -( WORK( INDUMN+I )*Z( I ) )
            IF( (ABS(Z(I))+ABS(Z(I+1)))* ABS(LD(I)).LT.GAPTOL )
     $         THEN
               Z( I+1 ) = ZERO
               ISUPPZ( 2 ) = I
               GO TO 260
            END IF
            ZTZ = ZTZ + DBLE( Z( I+1 )*Z( I+1 ) )
 250     CONTINUE
 260     CONTINUE
      ELSE
*        Run slower loop if NaN occurred.
         DO 270 I = R, BN - 1
            IF( Z( I ).EQ.ZERO ) THEN
               Z( I+1 ) = -( LD( I-1 ) / LD( I ) )*Z( I-1 )
            ELSE
               Z( I+1 ) = -( WORK( INDUMN+I )*Z( I ) )
            END IF
            IF( (ABS(Z(I))+ABS(Z(I+1)))* ABS(LD(I)).LT.GAPTOL )
     $           THEN
               Z( I+1 ) = ZERO
               ISUPPZ( 2 ) = I
               GO TO 280
            END IF
            ZTZ = ZTZ + DBLE( Z( I+1 )*Z( I+1 ) )
 270     CONTINUE
 280     CONTINUE
      END IF
*
*     Compute quantities for convergence test
*
      TMP = ONE / ZTZ
      NRMINV = SQRT( TMP )
      RESID = ABS( MINGMA )*NRMINV
      RQCORR = MINGMA*TMP
*
*
      RETURN
*
*     End of ZLAR1V
*
      END

max_stars_repo_name: IRFM/ALOHA

max_issues_repo_name: IRFM/ALOHA

max_forks_repo_name: IRFM/ALOHA


module fsutils
implicit none
contains

function expanduser(indir)

character(:), allocatable :: expanduser
character(*), intent(in) :: indir

! -- resolve home directory as Fortran does not understand tilde
! works for Linux, Mac, Windows and more

if (len_trim(indir) < 1) then
  stop 'must provide path to expand'
elseif (indir(1:1) /= '~') then
  expanduser = trim(adjustl(indir))
  return
elseif (len_trim(indir) < 3) then
  expanduser = homedir()
else
  expanduser = homedir() // trim(adjustl(indir(3:)))
endif

end function expanduser

function homedir()

character(:), allocatable :: homedir
character :: filesep
character(256) :: buf

! assume MacOS/Linux/BSD/Cygwin/WSL
filesep = '/'
call get_environment_variable("HOME", buf)

if (len_trim(buf) == 0) then  ! Windows
  call get_environment_variable("USERPROFILE", buf)
  filesep = char(92)
endif

homedir = trim(buf) // filesep

end function homedir

end module fsutils

!------ demo

program home

use fsutils
implicit none
! explores what '~' means for paths in Fortran
! Note: when testing, enclose argument in '~/test.txt' quotes or
!  shell will expand '~' before it gets to Fortran!

character(:), allocatable :: expanded
character(256) :: buf

call get_command_argument(1, buf)

expanded = expanduser(trim(buf))

print '(A)', trim(buf)
print '(A)', expanded

end program


text:
  function setupTarget (dataObj) result(errorFlag)

! ======================================================================
!
!     Neutron density, proton density, nucleon density .
!     Same functional form, different normalizations.
!     Determines zone boundaries in nucleus.
!     Finds average: Fermi momentum, Coulomb potential, nucleon density
!     in each zone of nucleus.
!
!   Rsm(ia) is the value of the radius at which the Woods-Saxon
!   nucleon density reaches the value of a(ia) * (the central value).
!   Hin is the integral from rsm(ia-1) (or 0) to rsm(ia) of
!   the density * r**2.
!   Sumhin is the integral of r**2/(1 + exp((r-rn)/bn) from 0 to
!   rsm(n)  [1/(4*pi) * (unnormalized) total mass inside of rsm(n)].
!   Rhon0 is the density normalization factor so all nucleons are inside
!   rsm(n).
!   By an integration by parts, the Coulomb potential at r is
!   proportional to:  1/r *Int[0,r]{r'**2 * rhop(r')}dr' +
!   Int[r,rsm(n)]{r' * rhop(r')}dr' .
!   Fis is proportional to this quantity * r**2.
!   Rhon is the average nucleon density between rsm(ia-1) and rsm(ia);
!   af is the average Coulomb potential (for a unit charge) in zone ia.
!   Tfp and tfn are the average Fermi momenta for protons and neutrons
!   in the nuclear zone ia.
!   Rbig is the normalized [to rsm(n)] radius of the outer boundary of
!   each zone.
!
! NOTE: Modified labeling according to:
!    (in options data type)
!       n    => dataObj%options%numZones
!       r0n  => dataObj%options%r0
!       bn   => dataObj%options%expDenom
!       rmax => dataObj%options%maxRad
!       a    => dataObj%options%aveDen
!    (in target data type)
!       rsm  => zoneBoundR
!       rbig => zoneBRFrac
!       rhop => protonDensity
!       rhon => neutronDensity
!       af   => coulombPote
!       tfp  => protFermiMom
!       tfn  => neutFermiMom
!
!
!    CEM95 written by S. G. Mashnik
!    Edited by A. J. Sierk   LANL  T-2,  February-March, 1996.
!   "Last" change: 14-AUG-2003 by NVMokhov
!    Edited by A. J. Sierk, LANL T-16, October, 2003.
!    Edited by AJS, LANL T-2, December, 2011.
!    Edited by LMK, XCP-3, July 2013 (included error protection).
!
! ======================================================================

    use, intrinsic:: iso_fortran_env, only: int32, real64
    use standardDCMDataParams, only: zro, one, thr, five, twpi, fourpi, &
         & esq, ato3rd, twthrd

    implicit none
    class(StandardDCMData), intent(inout) :: dataObj
    integer(int32)                        :: errorFlag

    integer(int32) :: ia, k, outerBound1
    real(real64)   :: endn, ennuc, fac, fi0, pden, rhon0, rn, &
         & sumhin, temp, tempR
    real(real64), dimension(maxZonesAllowed) :: hin=zro, fi=zro

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

    procedure(IntegralInterface), pointer :: sfisPtr => sfis
    procedure(IntegralInterface), pointer :: ss2Ptr => ss2

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

    errorFlag = properlyConstructed

    ! Obtain parameters for the density equation
    ia = nint(dataObj%target%numBaryons)
    rn = dataObj%options%r0*ato3rd(ia)

    ! Obtain fractions of each nucleon type in the nucleus
    ennuc = dataObj%target%numBaryons - dataObj%target%numProtons   ! Num neutrons
    pden = dataObj%target%numProtons/dataObj%target%numBaryons      ! Proton  fraction in nucleus
    endn = ennuc/dataObj%target%numBaryons      ! Neutron fraction in nucleus

    ! Obtain radial boundaries of each zone [fm]
    ! NOTE: Protection for IEEE errors is provided for in the constructor.
    do ia = 1, dataObj%options%numZones
       dataObj%target%zoneBoundR(ia) = rn + dataObj%options%expDenom * &
            & log(one/dataObj%options%aveDen(ia) - one)
    end do
    ! Last zone bound extends to the max allowed
    dataObj%target%zoneBoundR(dataObj%options%numZones+1) = dataObj%options%maxRad

    ! Obtain integral values regarding radial bounds for neutron density (hin) and fermi energy (fi)
    do ia = 1,dataObj%options%numZones
       if (ia == 1) then
          call sfint (sfisPtr, dataObj%target%zoneBoundR(ia), zro, &
               & dataObj%options%expDenom, rn, fi(ia), dataObj)
          call sfint (ss2Ptr, dataObj%target%zoneBoundR(ia),  zro, &
               & dataObj%options%expDenom, rn, hin(ia), dataObj)
       else
          call sfint (sfisPtr, dataObj%target%zoneBoundR(ia), dataObj%target%zoneBoundR(ia-1), &
               & dataObj%options%expDenom, rn, fi(ia), dataObj)
          call sfint (ss2Ptr, dataObj%target%zoneBoundR(ia), dataObj%target%zoneBoundR(ia-1), &
               & dataObj%options%expDenom, rn, hin(ia), dataObj)
       endif
    end do

    ! Sum all portions of integral
    sumhin = zro
    do ia = 1,dataObj%options%numZones
       sumhin = sumhin + hin(ia)
    end do

    temp = fourpi*sumhin   ! Assumed isotropic (hence *4pi)
    if (temp < div0Lim .and. temp > -div0Lim) then
       temp = div0Lim
       write(dataObj%io%message,1000) "108"
       call dataObj%io%print(4, 3, dataObj%io%message)
    end if
    ! Obtain base line neutron density
    rhon0 = thr*dataObj%target%numBaryons/(temp)
    temp = sumhin
    if (temp < div0Lim .and. temp > -div0Lim) then
       temp = div0Lim
       write(dataObj%io%message,1000) "114"
       call dataObj%io%print(4, 3, dataObj%io%message)
    end if
    ! Obtain base line neutron fermi energy
    fi0 = thr*esq*dataObj%target%numProtons/temp

    ! For all zones, establish neutron densities and fermi energies
    tempR = dataObj%target%zoneBoundR(dataObj%options%numZones)
    if ( tempR < div0Lim .and. temp > -div0Lim ) then
       write(dataObj%io%message,1000) "137"
       call dataObj%io%print(4, 3, dataObj%io%message)
       tempR = div0Lim
    end if
    do ia = 1,dataObj%options%numZones

       ! Obtain a volume component [1 / (deltaR**3)]
       if (ia == 1) then
          temp = dataObj%target%zoneBoundR(ia)
          if (temp < div0Lim .and. temp > -div0Lim) then
             temp = div0Lim
             write(dataObj%io%message,1000) "122"
             call dataObj%io%print(4, 3, dataObj%io%message)
          end if
          fac = one/temp**3
       else
          temp = dataObj%target%zoneBoundR(ia)**3 - dataObj%target%zoneBoundR(ia-1)**3
          if (temp < div0Lim .and. temp > -div0Lim) then
             temp = div0Lim
             write(dataObj%io%message,1000) "129"
             call dataObj%io%print(4, 3, dataObj%io%message)
          end if
          fac = one/(temp)
       endif

       ! Obtain neutron density
       dataObj%target%neutronDensity(ia) = hin(ia)*rhon0*fac
       dataObj%target%coulombPote(ia) = fi0*fi(ia)*fac

       ! Obtain fermi energies based on the respective nucleon fraction
       dataObj%target%protFermiMom(ia) = 0.198595d0 * &
            & (dataObj%target%neutronDensity(ia)*pden)**twthrd
       dataObj%target%neutFermiMom(ia) = 0.198322d0 * &
            & (dataObj%target%neutronDensity(ia)*endn)**twthrd

       ! Obtain fraction of zone bound to that of the last zone bound (excluding max)
       dataObj%target%zoneBRFrac(ia) = dataObj%target%zoneBoundR(ia)/tempR

    end do
    ! Obtain fraction of maximum zone bound to that of the last zone bound (excluding max)
    dataObj%target%zoneBRFrac(dataObj%options%numZones+1) = dataObj%options%maxRad/tempR

    ! Set outermost boundary and remaining bins for nucleon densities and Fermi energies
    outerBound1 = dataObj%options%numZones + 1
    dataObj%target%coulombPote(outerBound1)    = &
         & dataObj%target%numProtons*esq/dataObj%options%maxRad
    dataObj%target%protFermiMom(outerBound1)  = zro
    dataObj%target%neutFermiMom(outerBound1)  = zro
    dataObj%target%neutronDensity(outerBound1) = zro
    dataObj%target%protonDensity(outerBound1)  = zro
    k = dataObj%options%numZones + 2
    do ia = k, maxZonesAllowed
       dataObj%target%zoneBoundR(ia)     = zro
       dataObj%target%zoneBRFrac(ia)     = zro
       dataObj%target%neutronDensity(ia) = zro
       dataObj%target%protonDensity(ia)  = zro
       dataObj%target%coulombPote(ia)    = zro
       dataObj%target%protFermiMom(ia)  = zro
       dataObj%target%neutFermiMom(ia)  = zro
    end do

    ! Scale neutron/proton zone densities according to the neutron/proton fraction
    do ia = 1, dataObj%options%numZones
       dataObj%target%protonDensity(ia) = dataObj%target%neutronDensity(ia)*pden
       dataObj%target%neutronDensity(ia) = dataObj%target%neutronDensity(ia)*endn
    end do

    ! Obtain geometric cross section:
    dataObj%target%geomCrossSection = twpi * five * &
         & ( dataObj%target%zoneBoundR( dataObj%options%numZones )**2 )

    return
! ======================================================================
1000 format("Divide by zero error prevented in 'setuptarget.f90', ", &
          & "line(s) ", A)
! ======================================================================
  end function setupTarget

max_stars_repo_name: lanl/generalized-spallation-model

max_stars_repo_licenses: Intel Unlicense

max_issues_repo_name: lanl/generalized-spallation-model

max_issues_repo_licenses: Intel Unlicense

max_forks_repo_name: lanl/generalized-spallation-model

max_forks_repo_licenses: Intel Unlicense


text:
MODULE precision_05
     EXTERNAL MPI_WTIME
     REAL MPI_WTIME
END MODULE precision_05

subroutine timer()
  use precision_05
  implicit none

  real t1wc

! This reference to the function MPI_WTIME() causes an error in ROSE.
  t1wc = MPI_WTIME()

end subroutine timer

max_stars_repo_name: maurizioabba/rose

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
!------------------------------------------------------------------------------
! NAME:
!       Create_ProcessControl_File
!
! PURPOSE:
!       Program to create a generic Process Control file (profile and molecule
!       set independent) for the transmittance production software given an
!       input list of netCDF SRF data files.
!
! CATEGORY:
!       Transmittance Production
!
! LANGUAGE:
!       Fortran-95
!
! MODULES:
!       Type_Kinds:                Module containing definitions for kinds
!                                  of variable types.
!
!       File_Utility:              Module containing generic file utility routines
!
!       Message_Handler:             Module to define simple error codes and
!                                  handle error conditions
!                                  USEs: FILE_UTILITY module
!
!       SensorInfo_Define:         Module defining the SensorInfo data
!                                  structure and containing routines to
!                                  manipulate it.
!                                  USEs: TYPE_KINDS module
!                                        ERROR_HANDLER module
!
!       SensorInfo_LinkedList:     Module defining the SensorInfo Linked
!                                  List data structure and containing
!                                  routines to manipulate it.
!                                  USEs: TYPE_KINDS module
!                                        ERROR_HANDLER module
!                                        SENSORINFO_DEFINE module
!
!       SensorInfo_IO:             Module continaing routines to read and
!                                  write ASCII SensorInfo format files.
!                                  USEs: TYPE_KINDS module
!                                        FILE_UTILITY module
!                                        ERROR_HANDLER module
!                                        SENSORINFO_DEFINE module
!
!       SRF_Define:                Module defining the SRF data structure and
!                                  containing routines to manipulate it.
!                                  USEs: TYPE_KINDS module
!                                        ERROR_HANDLER module
!
!       SRF_netCDF_IO:             Module containing routines to read and write
!                                  netCDF SRF format files.
!                                  USEs: TYPE_KINDS module
!                                        ERROR_HANDLER module
!                                        SRF_DEFINE module
!                                        NETCDF module
!                                        NETCDF_UTILITY module
!
!       ProcessControl_Define:     Module containing the ProcessControl
!                                  data type definition and routines to
!                                  manipulate the structure.
!                                  USEs: TYPE_KINDS module
!                                        ERROR_HANDLER module
!
!       ProcessControl_IO:         Module containing routines to read and write
!                                  ASCII Process Control files.
!                                  USEs: TYPE_KINDS module
!                                        FILE_UTILITY module
!                                        ERROR_HANDLER module
!                                        PROCESSCONTROL_DEFINE module
!
!       Tau_Production_Parameters: Module defining parameters used in the LBL
!                                  transmittance production runs
!                                  USEs: TYPE_KINDS module
!                                        LBLRTM_PARAMETERS module
!
!       Tau_Production_Utility:    Module continaing utility routines for the LBL
!                                  transmittance production runs.
!                                  USEs: TYPE_KINDS module
!                                        FILE_UTILITY module
!                                        ERROR_HANDLER module
!                                        COMPARE_FLOAT_NUMBERS module
!                                        TRANSMITTANCE_PRODUCTION_PARAMETERS module
!
! CONTAINS:
!       None.
!
! INCLUDE FILES:
!       None.
!
! EXTERNALS:
!       None.
!
! COMMON BLOCKS:
!       None.
!
! FILES ACCESSED:
!       Input:  - SensorInfo file.
!               - netCDF format SRF data file(s).
!
!       Output: ASCII ProcessControl file.
!
! SIDE EFFECTS:
!       If the ProcessControl file already exists, it is overwritten.
!
! RESTRICTIONS:
!       The maximum number of input SRF files that can be processed
!       is limited to 100.
!
! CREATION HISTORY:
!       Written by:     Paul van Delst, CIMSS/SSEC 26-Apr-2002
!                       paul.vandelst@ssec.wisc.edu
!
!  Copyright (C) 2002 Paul van Delst
!
!
!------------------------------------------------------------------------------

PROGRAM Create_ProcessControl_File

  ! ------------
  ! Module usage
  ! ------------

  USE Type_Kinds
  USE File_Utility
  USE Message_Handler

  USE SensorInfo_Define
  USE SensorInfo_LinkedList
  USE SensorInfo_IO

  USE oSRF_File_Define

  USE ProcessControl_Define
  USE ProcessControl_IO

  USE Tau_Production_Parameters
  USE Tau_Production_Utility
  USE Compare_Float_Numbers

  ! -----------------------
  ! Disable implicit typing
  ! -----------------------

  IMPLICIT NONE

  ! ----------
  ! Parameters
  ! ----------

  CHARACTER( * ), PARAMETER :: PROGRAM_NAME = 'Create_ProcessControl_File'
  CHARACTER( * ), PARAMETER :: PROGRAM_RCS_ID = &
  CHARACTER( * ), PARAMETER :: PROGRAM_HEADER = &
  '**********************************************************'

  ! -- Maximum number of sensors to allow
  INTEGER, PARAMETER :: MAX_N_SENSORS = 150

  ! -- Integer keyword set
  INTEGER, PARAMETER :: SET = 1

  ! ---------
  ! Variables
  ! ---------

  CHARACTER( 256 ) :: Message

  INTEGER         :: pn_pos
  CHARACTER( 80 ) :: pn_fmt

  CHARACTER( 512 ) :: SensorInfo_Filename
  CHARACTER( 512 ) :: SRF_FilePath
  CHARACTER( 512 ) :: SRF_Filename
  CHARACTER( 512 ) :: ProcessControl_Filename

  INTEGER :: IO_Status
  INTEGER :: Error_Status
  INTEGER :: Allocate_Status

  INTEGER,          DIMENSION( MAX_N_SENSORS ) :: Available_Sensor
  CHARACTER( 512 ), DIMENSION( MAX_N_SENSORS ) :: Available_SRF_Filename

  INTEGER :: i, l, ch, ll

  INTEGER :: n_Available_Sensors
  INTEGER :: n_Channels
  INTEGER :: Valid_Channel

  INTEGER,         DIMENSION( : ), ALLOCATABLE :: n_SRF_Frequencies
  REAL( fp_kind ), DIMENSION( : ), ALLOCATABLE :: Begin_SRF_Frequency
  REAL( fp_kind ), DIMENSION( : ), ALLOCATABLE :: End_SRF_Frequency
  TYPE(oSRF_File_type),  ALLOCATABLE :: oSRF_File(:)

  REAL( fp_kind ) :: dF, dF1
  INTEGER         :: Idx, dF_Index

  INTEGER :: n_Sensors, n
  TYPE( SensorInfo_type )      :: SensorInfo
  TYPE( SensorInfo_List_type ) :: SensorInfo_List

  TYPE( ProcessControl_type ) :: ProcessControl

  !#----------------------------------------------------------------------------#
  !#                   -- CREATE THE SensorInfo LINKED LIST --                  #
  !#----------------------------------------------------------------------------#

  SensorInfo_List = New_SensorInfo_List()

  !#----------------------------------------------------------------------------#
  !#                       -- OUTPUT DESCRIPTIVE HEADER --                      #
  !#----------------------------------------------------------------------------#

  pn_pos = ( LEN( PROGRAM_HEADER ) / 2 ) - &
           ( LEN( PROGRAM_NAME ) / 2 )
  pn_pos = MAX( pn_pos, 0 ) + 5
  WRITE( pn_fmt, '( "( ",i2,"x, a )" )' ) pn_pos

  WRITE( *, '(/5x, a )' ) PROGRAM_HEADER
  WRITE( *, FMT = TRIM( pn_fmt ) ) PROGRAM_NAME
  WRITE( *, '(/5x, " Program to create a transmittance production Process")' )
  WRITE( *, '( 5x, "   Control file from netCDF format SRF files.")' )
  WRITE( *, '(/5x, " $Revision: 2.3 $")' )
  WRITE( *, '( 5x, a, / )' ) PROGRAM_HEADER

  !#----------------------------------------------------------------------------#
  !#                       -- READ THE SensorInfo FILE --                       #
  !#----------------------------------------------------------------------------#

  ! ---------------------------
  ! Get the SensorInfo filename
  ! ---------------------------

  WRITE( *, FMT     = '( /5x, "Enter a SensorInfo filename: " )', &
            ADVANCE = 'NO' )
  READ( *, FMT = '( a )' ) SensorInfo_Filename
  SensorInfo_Filename = ADJUSTL( SensorInfo_Filename )

  ! ------------------------
  ! Read the SensorInfo file
  ! ------------------------

  Error_Status = Read_SensorInfo( SensorInfo_Filename, &
                                  SensorInfo_List, &
                                  Quiet = SET )

  IF ( Error_Status /= SUCCESS ) THEN
    CALL Display_Message( PROGRAM_NAME, &
                          'Error reading SensorInfo file '//&
                          TRIM( SensorInfo_Filename ), &
                          FAILURE )
    STOP
  END IF

  ! -- Count the number of sensors
  n_Sensors = Count_SensorInfo_Nodes( SensorInfo_List )

  IF ( n_Sensors < 1 ) THEN
    CALL Display_Message( PROGRAM_NAME, &
                          'SensorInfo_List is empty.', &
                          FAILURE )
    STOP
  END IF

  !#----------------------------------------------------------------------------#
  !#                       -- GET DF index --                                   #
  !#----------------------------------------------------------------------------#

  WRITE( *, FMT     = '( /5x, "Enter the Frequency intervals index: " )', &
            ADVANCE = 'NO' )
  READ( *, FMT = '( I1 )' )  dF_Index

  !#----------------------------------------------------------------------------#
  !#                       -- GET SRF DATA FILE PATH --                         #
  !#----------------------------------------------------------------------------#

  WRITE( *, FMT     = '( /5x, "Enter the SRF data file path [e.g. ./SRF_Data/]: " )', &
            ADVANCE = 'NO' )
  READ( *, FMT = '( a )' ) SRF_FilePath

  SRF_FilePath = ADJUSTL( SRF_FilePath )

  !#----------------------------------------------------------------------------#
  !#             -- LOOP OVER SENSORS TO COUNT THE TOTAL NUMBER --              #
  !#             --     OF AVAILABLE SENSORS AND CHANNELS       --              #
  !#----------------------------------------------------------------------------#

  ! ----------------------------------------------------
  ! Initialise the available sensor and channel counters
  ! ----------------------------------------------------

  n_Available_Sensors = 0
  n_Channels          = 0

  ! -----------------
  ! Begin sensor loop
  ! -----------------

  Sensor_Count_Loop: DO n = 1, n_Sensors

    ! ---------------------------------------------
    ! Get the current SensorInfo data from the list
    ! ---------------------------------------------

    Error_Status = GetFrom_SensorInfo_List( SensorInfo_List, &
                                            n, &
                                            SensorInfo )

    IF ( Error_Status /= SUCCESS ) THEN
      Error_Status = FAILURE
      WRITE( Message, '( "Error retrieving SensorInfo data for sensor # ", i5, &
                        &" in count loop." )' ) n
      CALL Display_Message( PROGRAM_NAME, &
                            TRIM( Message ), &
                            Error_Status )
      STOP
    END IF

    ! --------------------------
    ! Construct the SRF filename
    ! --------------------------

    SRF_Filename = TRIM( SRF_FilePath )//TRIM( SensorInfo%Sensor_Id )//'.osrf.nc'

    !#--------------------------------------------------------------------------#
    !#                   -- ONLY INCLUDE FILES THAT EXIST --                    #
    !#--------------------------------------------------------------------------#

    Count_Available_Sensors: IF ( File_Exists( TRIM( SRF_Filename ) ) ) THEN

      ! -----------------------------------------
      ! Increment the number of available sensors
      ! -----------------------------------------

      n_Available_Sensors = n_Available_Sensors + 1

      IF ( n_Available_Sensors > MAX_N_SENSORS ) THEN
        WRITE( Message, '( "Maximum number of sensors, ", i3, ", exceeded." )' ) MAX_N_SENSORS
        CALL Display_Message( PROGRAM_NAME, &
                              TRIM( Message ), &
                              FAILURE )
        STOP
      END IF

      ! --------------------------------
      ! Increment the number of channels
      ! --------------------------------

      n_Channels = n_Channels + SensorInfo%n_Channels

      ! -------------------------------
      ! Save the SensorInfo node number
      ! and SRF filename of this sensor
      ! -------------------------------

      Available_Sensor( n_Available_Sensors ) = n
      Available_SRF_Filename( n_Available_Sensors ) = TRIM( SRF_Filename )

    END IF Count_Available_Sensors

    !#--------------------------------------------------------------------------#
    !#              -- DESTROY THE CURRENT SensorInfo STRUCTURE --              #
    !#--------------------------------------------------------------------------#

    Error_Status = Destroy_SensorInfo( SensorInfo )

    IF ( Error_Status /= SUCCESS ) THEN
      WRITE( Message, '( "Error destroying SensorInfo for sensor # ", i5, &
                        &" in count loop." )' ) n
      CALL Display_Message( PROGRAM_NAME, &
                            TRIM( Message ), &
                            Error_Status )
      STOP
    END IF

  END DO Sensor_Count_Loop

  !#----------------------------------------------------------------------------#
  !#                      -- CHECK THAT WE HAVE SOMETHING --                    #
  !#----------------------------------------------------------------------------#

  IF ( n_Available_Sensors < 1 .OR. n_Channels < 1 ) THEN
    CALL Display_Message( PROGRAM_NAME, &
                          'No available sensors of channels!', &
                          FAILURE )
    STOP
  END IF

  !#----------------------------------------------------------------------------#
  !#                -- ALLOCATE THE ProcessControl STRUCTURE --                 #
  !#----------------------------------------------------------------------------#

  Error_Status = Allocate_ProcessControl( n_Available_Sensors, &
                                          n_Channels, &
                                          ProcessControl )

  IF ( Error_Status /= SUCCESS ) THEN
    CALL Display_Message( PROGRAM_NAME, &
                          'Error allocating ProcessControl structure.', &
                          FAILURE )
    STOP
  END IF

  ALLOCATE( oSRF_File( ProcessControl%n_Files ),  &
            STAT = Allocate_Status )

  IF ( Allocate_Status /= 0 ) THEN
    WRITE( Message, '( "Error allocating oSRF_File array. STAT = ", i5 )' ) &
                    Allocate_Status
    CALL Display_Message( PROGRAM_NAME, &
                          TRIM( Message ), &
                          FAILURE )
    STOP
  END IF

  !#----------------------------------------------------------------------------#
  !#              -- LOOP OVER AVAILABLE SENSORS TO READ DATA --                #
  !#              --  AND FILL THE ProcessControl STRUCTURE   --                #
  !#----------------------------------------------------------------------------#

  ! -------------------------------------
  ! Initialise main channel counter index
  ! -------------------------------------

  l = 0

  ! -----------------
  ! Begin sensor loop
  ! -----------------

  Sensor_Read_Loop: DO n = 1, n_Available_Sensors

    ! ----------------------------------------------
    ! Get the required SensorInfo data from the list
    ! ----------------------------------------------

    Error_Status = GetFrom_SensorInfo_List( SensorInfo_List, &
                                            Available_Sensor(n), &
                                            SensorInfo )

    IF ( Error_Status /= SUCCESS ) THEN
      Error_Status = FAILURE
      WRITE( Message, '( "Error retrieving SensorInfo data for sensor # ", i5, &
                        &" in read loop." )' ) Available_Sensor(n)
      CALL Display_Message( PROGRAM_NAME, &
                            TRIM( Message ), &
                            Error_Status )
      STOP
    END IF

    ! ---------------------------------
    ! Allocate the SRF frequency arrays
    ! ---------------------------------

    ALLOCATE( n_SRF_Frequencies(   SensorInfo%n_Channels ), &
              Begin_SRF_Frequency( SensorInfo%n_Channels ), &
              End_SRF_Frequency(   SensorInfo%n_Channels ), &
              STAT = Allocate_Status )

    IF ( Allocate_Status /= 0 ) THEN
      CALL Display_Message( PROGRAM_NAME, &
                            'Error allocating frequency arrays for SRF netCDF file '//&
                            TRIM( Available_SRF_Filename(n) )//' read.', &
                            FAILURE )
      STOP
    END IF

    ! -----------------------
    ! Read the frequency data
    ! -----------------------

    Error_Status = oSRF_File_Read( &
                     oSRF_File(n)    , &
                     Available_SRF_Filename(n)  )

    IF ( Error_Status /= SUCCESS ) THEN
      CALL Display_Message( PROGRAM_NAME, &
                            'Error reading SRF file '//&
                            TRIM( Available_SRF_Filename(n) ), &
                            Error_Status )
      STOP
    END IF

    ! For Infrared sensor, only have one band for each channel
    DO ll = 1, SensorInfo%n_Channels
      n_SRF_Frequencies(ll) = oSRF_File(n)%oSRF(ll)%n_Points(1)
      Begin_SRF_Frequency(ll) =oSRF_File(n)%oSRF(ll)%f1(1)
      End_SRF_Frequency(ll) = oSRF_File(n)%oSRF(ll)%f2(oSRF_File(n)%oSRF(ll)%n_Bands)
    ENDDO

    ! -----------------------------
    ! Assign the file prefix string
    ! -----------------------------

    ProcessControl%File_Prefix(n) = TRIM( SensorInfo%Sensor_Id )

    ! --------------------------------------
    ! Determine the frequency interval index
    ! --------------------------------------

    Frequency_Interval_Loop: DO ch = 1, SensorInfo%n_Channels

      ! -- Compute the channel frequency interval
      dF = ( End_SRF_Frequency(ch) - Begin_SRF_Frequency(ch) ) / &
      !    ---------------------------------------------------
               REAL( n_SRF_Frequencies(ch) - 1, fp_kind )

      ! -- Compute the frequency interval index

      ! -- Is it valid?
      IF ( .NOT. Compare_Float( dF, FREQUENCY_INTERVAL(dF_Index), ULP = 10000 ) ) THEN
        WRITE( Message, '( "Invalid frequency interval, ", es13.6, &
                          &"cm^-1, detected for channel ", i5, &
                          &" SRF from ", a, "." )' ) &
                        dF, &
                        SensorInfo%Sensor_Channel(ch), &
                        TRIM( Available_SRF_Filename(n) )
        CALL Display_Message( PROGRAM_NAME, &
                              TRIM( Message ), &
                              FAILURE )
        STOP
      END IF

      ! -- Save and compare the dF
      IF ( ch == 1 ) THEN
        dF1 = dF
      ELSE
        IF ( .NOT. Compare_Float( dF, FREQUENCY_INTERVAL(dF_Index), ULP = 10000 ) ) THEN
          WRITE( Message, '( "Frequency interval, ", es13.6, &
                            &"cm^-1, for channel ", i5, &
                            &" is different from that for the other channels, ", es13.6, &
                            &"cm^-1, for ", a, "." )' ) &
                          dF1, &
                          SensorInfo%Sensor_Channel(ch), &
                          FREQUENCY_INTERVAL( dF_Index ), &
                          TRIM( Available_SRF_Filename(n) )
          CALL Display_Message( PROGRAM_NAME, &
                                TRIM( Message ), &
                                FAILURE )
          STOP
        END IF
      END IF

    END DO Frequency_Interval_Loop

    ! -- Assign the frequency interval index
    ProcessControl%dF_Index(n) = dF_Index

    ! ----------------------------------------
    ! Begin loop over channels for this sensor
    ! ----------------------------------------

    Valid_Channel = 0

    Channel_Loop: DO ch = 1, SensorInfo%n_Channels

      ! -------------------------------------------
      ! Only process current channel if it is valid
      ! -------------------------------------------

      IF ( ( Begin_SRF_Frequency(ch) >= FREQUENCY_BEGIN ) .AND. &
           (   End_SRF_Frequency(ch) <= FREQUENCY_END   )       ) THEN

        ! -- Update main channel counter index for ALL available sensors
        l = l + 1

        ! -- Update valid channel counter for CURRENT sensor
        Valid_Channel = Valid_Channel + 1

        ! -- Assign the begin and end VALID channel indices for CURRENT sensor
        IF ( Valid_Channel == 1 ) THEN
          ProcessControl%Channel_Index( 1, n ) = l
        END IF
        ProcessControl%Channel_Index( 2, n ) = l

        ! -- Assign the channel list data
        ProcessControl%List(l)%Channel       = SensorInfo%Sensor_Channel(ch)
        ProcessControl%List(l)%Begin_LBLband = Compute_LBL_band( Begin_SRF_Frequency(ch), &
                                                                 ProcessControl%dF_Index(n)  )
        ProcessControl%List(l)%End_LBLband   = Compute_LBL_band( End_SRF_Frequency(ch), &
                                                                 ProcessControl%dF_Index(n)  )

        ! -- Initialise processing flag
        ProcessControl%List(l)%Processed = 0

      END IF

    END DO Channel_Loop

    ! -----------------------------------
    ! Deallocate the SRF frequency arrays
    ! -----------------------------------

    DEALLOCATE( n_SRF_Frequencies, &
                Begin_SRF_Frequency, &
                End_SRF_Frequency, &
                STAT = Allocate_Status )

    IF ( Allocate_Status /= 0 ) THEN
      CALL Display_Message( PROGRAM_NAME, &
                            'Error deallocating frequency arrays for SRF netCDF file '//&
                            TRIM( Available_SRF_Filename(n) )//' read.', &
                            FAILURE )
      STOP
    END IF

    ! ----------------------------------------
    ! Destroy the current SensorInfo structure
    ! ----------------------------------------

    Error_Status = Destroy_SensorInfo( SensorInfo )

    IF ( Error_Status /= SUCCESS ) THEN
      WRITE( Message, '( "Error destroying SensorInfo for sensor # ", i5, &
                        &" in count loop." )' ) n
      CALL Display_Message( PROGRAM_NAME, &
                            TRIM( Message ), &
                            Error_Status )
      STOP
    END IF

  END DO Sensor_Read_Loop

  ! -----------------------------------------------------
  ! Update the total channel count. Some channels may not
  ! lie completely within the processing frequency bounds
  ! -----------------------------------------------------

  ProcessControl%n_Channels = l

  !#----------------------------------------------------------------------------#
  !#                   -- WRITE THE ProcessControl FILE --                      #
  !#----------------------------------------------------------------------------#

  ! ----------------
  ! Get the filename
  ! ----------------

  WRITE( *, FMT     = '( /5x, "Enter the Process Control output filename : " )', &
            ADVANCE = 'NO' )
  READ( *, FMT = '( a )' ) ProcessControl_Filename
  ProcessControl_Filename = ADJUSTL( ProcessControl_Filename )

  ! --------------
  ! Write the data
  ! --------------

  Error_Status = Write_ProcessControl( TRIM( ProcessControl_Filename ), &
                                       ProcessControl )

  IF ( Error_Status /= SUCCESS ) THEN
    CALL Display_Message( PROGRAM_NAME, &
                          'Error writing Process Control data to '//&
                          TRIM( ProcessControl_Filename )//'.', &
                          Error_Status )
    STOP
  END IF

  !#----------------------------------------------------------------------------#
  !#                  -- DESTROY THE SensorInfo LINKED LIST --                  #
  !#----------------------------------------------------------------------------#

  Error_Status = Destroy_SensorInfo_List( SensorInfo_List )

  IF ( Error_Status /= SUCCESS ) THEN
    CALL Display_Message( PROGRAM_NAME, &
                          'Error destroying SensorInfo linked list.', &
                          WARNING )
  END IF

  !#----------------------------------------------------------------------------#
  !#                   -- TEST READ OF ProcessControl FILE --                   #
  !#----------------------------------------------------------------------------#

  ! ------------------
  ! Read the data file
  ! ------------------

  Error_Status = Read_ProcessControl( TRIM( ProcessControl_Filename ), &
                                      ProcessControl )

  IF ( Error_Status /= SUCCESS ) THEN
    CALL Display_Message( PROGRAM_NAME, &
                          'Error reading ProcessControl file '//&
                          TRIM( ProcessControl_Filename )//'.', &
                          FAILURE )
    STOP
  END IF

  ! ----------------
  ! Output some data
  ! ----------------

  CALL oSRF_File_Destroy( oSRF_File )
  DEALLOCATE(oSRF_File, STAT = Allocate_Status )

  IF ( Allocate_Status /= 0 ) THEN
    WRITE( Message, '( "Error deallocating oSRF_File array. STAT = ", i5 )' ) &
                    Allocate_Status
    CALL Display_Message( PROGRAM_NAME, &
                          TRIM( Message ), &
                          WARNING )
  END IF

  ! ------------------------------------
  ! Destroy the ProcessControl structure
  ! ------------------------------------

  Error_Status = Destroy_ProcessControl( ProcessControl )

  IF ( Error_Status /= SUCCESS ) THEN
    CALL Display_Message( PROGRAM_NAME, &
                          'Error destroying ProcessControl data structure.', &
                          Error_Status )
    STOP
  END IF

  !#----------------------------------------------------------------------------#
  !#         -- CREATE A SIGNAL FILE INDICATING SUCCESSFUL COMPLETION --        #
  !#----------------------------------------------------------------------------#

  Error_Status = Create_Signal_File( TRIM( ProcessControl_Filename ) )

END PROGRAM Create_ProcessControl_File

max_stars_repo_name: hsbadr/crtm

max_issues_repo_name: hsbadr/crtm

max_forks_repo_name: hsbadr/crtm


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

!> Implementation of the dimer projection method for extended tight binding methods.
module frag
   use mctc_env, only : error_type, fatal_error, get_argument, wp
   use mctc_io, only : structure_type, new, write_structure
   use mctc_io_math, only : matinv_3x3
   use mctc_io_convert, only : autoev
   use frag_bondorder, only : get_wiberg_bondorder
   use frag_fragment, only : get_wiberg_fragment, unfold_fragment
   use frag_output, only : format_list, to_string
   use frag_rad, only : guess_bonds
   use frag_xtb, only : get_calculator
   use tblite_basis_type, only : get_cutoff, basis_type
   use tblite_blas, only : dot, gemv, gemm
   use tblite_context_type, only : context_type
   use tblite_cutoff, only : get_lattice_points
   use tblite_integral_overlap, only : get_overlap
   use tblite_output_format, only : format_string
   use tblite_wavefunction_type, only : wavefunction_type, new_wavefunction, &
      & get_density_matrix
   use tblite_xtb_calculator, only : xtb_calculator
   use tblite_xtb_singlepoint, only : xtb_singlepoint
   implicit none
   private

   public :: get_jab, jab_input

   !> Configuration data for calculation
   type :: jab_input
      !> Name of the requested tight binding method
      character(len=:), allocatable :: method
      !> List of fragments, generated if not given here
      integer, allocatable :: fragment(:)
      !> Threshold for generating fragments from bond orders
      real(wp) :: thr = 0.1_wp
      !> Calculation accuracy
      real(wp) :: accuracy = 1.0_wp
      !> Output verbosity
      integer :: verbosity = 2
      !> Electronic temperature in Kelvin
      real(wp) :: etemp = 300.0_wp
   end type jab_input

   !> Conversion factor from temperature to energy (Boltzmann's constant in atomic units)
   real(wp), parameter :: ktoau = 3.166808578545117e-06_wp

contains

!> Entry point for calculation of dipole projection related properties
subroutine get_jab(input, mol, error)
   !> Input data for the calculation
   type(jab_input), intent(in) :: input
   !> Molecular structure data
   type(structure_type), intent(inout) :: mol
   !> Error handling
   type(error_type), allocatable, intent(out) :: error

   integer :: spin, charge, stat, unit, ifr, nfrag, nao
   logical :: exist
   real(wp) :: energy, cutoff, jab, sab, jeff, invlat(3, 3)
   real(wp), allocatable :: gradient(:, :), sigma(:, :)
   type(context_type) :: ctx
   type(xtb_calculator) :: calc, fcalc
   type(structure_type), allocatable :: mfrag(:)
   type(wavefunction_type) :: wfn
   type(wavefunction_type), allocatable :: wfx(:)
   real(wp), allocatable :: overlap(:, :), trans(:, :), wbo(:, :), abc(:, :)
   real(wp), allocatable :: orbital(:, :), scmat(:, :), fdim(:, :), scratch(:), efrag(:)
   integer, allocatable :: fragment(:)

   call get_calculator(calc, mol, input%method, error)
   if (allocated(error)) return

   call new_wavefunction(wfn, mol%nat, calc%bas%nsh, calc%bas%nao, &
      & input%etemp * ktoau)

   call xtb_singlepoint(ctx, mol, calc, wfn, input%accuracy, energy, &
      & verbosity=input%verbosity-1)
   if (ctx%failed()) then
      call ctx%get_error(error)
      return
   end if

   allocate(overlap(calc%bas%nao, calc%bas%nao))
   cutoff = get_cutoff(calc%bas)
   call get_lattice_points(mol%periodic, mol%lattice, cutoff, trans)
   call get_overlap(mol, trans, cutoff, calc%bas, overlap)

   if (allocated(input%fragment)) then
      fragment = input%fragment
   else
      allocate(wbo(mol%nat, mol%nat))
      call get_wiberg_bondorder(calc%bas, overlap, wfn%density, wbo)

      allocate(fragment(mol%nat))
      call get_wiberg_fragment(fragment, wbo, input%thr)
   end if

   nfrag = maxval(fragment)
   select case(nfrag)
   case(:1)
      !call fatal_error(error, "Found no fragments in the input structure")
      !return
   case(2:)
      call ctx%message("Found "//to_string(nfrag)//" fragments:")
      do ifr = 1, nfrag
         call ctx%message("Fragment "//to_string(ifr)//":  "//format_list(fragment == ifr))
      end do
      call ctx%message("")
   end select

   call guess_bonds(mol, wbo)

   call get_wiberg_fragment(fragment, wbo, input%thr)

   nfrag = maxval(fragment)
   select case(nfrag)
   case(:1)
      !call fatal_error(error, "Found no fragments in the input structure")
      !return
   case(2:)
      call ctx%message("Found "//to_string(nfrag)//" fragments:")
      do ifr = 1, nfrag
         call ctx%message("Fragment "//to_string(ifr)//":  "//format_list(fragment == ifr))
      end do
      call ctx%message("")
   end select

   if (any(mol%periodic)) then
      invlat = matinv_3x3(mol%lattice)
      abc = matmul(invlat, mol%xyz)
      call unfold_fragment(abc, wbo, input%thr)
      mol%xyz(:, :) = matmul(mol%lattice, abc)
   end if

   allocate(mfrag(nfrag))
   call get_structure_fragment(mfrag(1), mol, fragment == 1)
   call write_structure(mfrag(1), "final.xyz", error)

end subroutine get_jab

!> Unpack coefficients from a fragment orbital space in the full orbital space
subroutine unpack_coeff(full_bas, frag_bas, full, frag, mask)
   !> Basis set information for full system
   type(basis_type), intent(in) :: full_bas
   !> Basis set information for fragment
   type(basis_type), intent(in) :: frag_bas
   !> Quantity in full orbital space
   real(wp), intent(out) :: full(:)
   !> Quantity in fragment orbital space
   real(wp), intent(in) :: frag(:)
   !> Atom resolved mask for this fragment
   logical, intent(in) :: mask(:)

   integer :: iat, jat, ish, ii, jj, iao, jao, nao

   jat = 0
   do iat = 1, size(mask)
      if (mask(iat)) then
         jat = jat + 1
         ii = full_bas%ish_at(iat)
         jj = frag_bas%ish_at(jat)
         do ish = 1, full_bas%nsh_at(iat)
            iao = full_bas%iao_sh(ii+ish)
            jao = frag_bas%iao_sh(jj+ish)
            nao = full_bas%nao_sh(ii+ish)
            full(iao+1:iao+nao) = frag(jao+1:jao+nao)
         end do
      else
         ii = full_bas%ish_at(iat)
         do ish = 1, full_bas%nsh_at(iat)
            iao = full_bas%iao_sh(ii+ish)
            nao = full_bas%nao_sh(ii+ish)
            full(iao+1:iao+nao) = 0.0_wp
         end do
      end if
   end do
end subroutine unpack_coeff

!> Extract a the fragment structure from the full structure
!>
!> Todo: This routine can currently only create neutral fragments,
!>       charge and spin information is not partition or transferred to fragments
subroutine get_structure_fragment(frag, mol, mask)
   !> Molecular structure data of the fragment
   type(structure_type), intent(out) :: frag
   !> Molecular structure data of the full system
   type(structure_type), intent(in) :: mol
   !> Atom resolved mask for this fragment
   logical, intent(in) :: mask(:)

   integer :: nat
   integer, allocatable :: num(:)
   real(wp), allocatable :: xyz(:, :)

   nat = count(mask)
   num = pack(mol%num(mol%id), mask)
   xyz = reshape(pack(mol%xyz, spread(mask, 1, 3)), [3, nat])
   call new(frag, num, xyz)
end subroutine get_structure_fragment

end module frag

max_stars_repo_name: awvwgk/frag-tool

max_issues_repo_name: awvwgk/frag-tool

max_forks_repo_name: awvwgk/frag-tool