text:
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      PROGRAM PROJ1
C     Project 1: Scattering by a central 6-12 potential
C  COMPUTATIONAL PHYSICS (FORTRAN VERSION)
C  by Steven E. Koonin and Dawn C. Meredith
C  Copyright 1989, Addison-Wesley Publishing Company Inc.
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      CALL INIT           !display header screen, setup parameters
5     CONTINUE            !main loop; execute once for each set of param
         CALL PARAM       !get input parameters
         CALL ARCHON      !calculate scattering angles
      GOTO 5
      END
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE ARCHON
C calculate scattering angles as a function of the impact parameter
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C Global variables
      INCLUDE 'PARAM.P1'
      INCLUDE 'IO.ALL'
C Local variables:
      REAL RMIN                 !radius of closest approach
      REAL ANGLE(0:MAX)         !final scattering angle
      REAL IMPACT(0:MAX)        !impact parameter
      REAL INT1, INT2           !values of the two integrals
      INTEGER IB                !index of impact param
      INTEGER NLINES            !number of lines printed so far
      INTEGER SCREEN            !send to terminal
      INTEGER PAPER             !make a hardcopy
      INTEGER FILE              !send to a file
      DATA SCREEN,PAPER,FILE/1,2,3/
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     output summary of parameters
      IF (TTERM) CALL PRMOUT(OUNIT,NLINES)
      IF (TFILE) CALL PRMOUT(TUNIT,NLINES)
      IF (GFILE) CALL PRMOUT(GUNIT,NLINES)
C
      DO 25 IB=0,(NB-1)                      !loop over impact param
C
        IMPACT(IB)=BMIN+(IB*DB)              !current value of B
        CALL THETA(IMPACT(IB),INT1,INT2,RMIN)!calculate THETA for each B
        ANGLE(IB)=(2*IMPACT(IB)*(INT1-INT2)*180)/PI !convert to degrees
C
C       text output
        IF (TTERM) CALL TXTOUT(OUNIT,NLINES,IMPACT(IB),RMIN,ANGLE(IB))
        IF (TFILE) CALL TXTOUT(TUNIT,NLINES,IMPACT(IB),RMIN,ANGLE(IB))
C
   25 CONTINUE
C
      IF (TTERM) CALL PAUSE('to continue...',1)
      IF (TTERM) CALL CLEAR
C
      IF (GTERM) CALL GRFOUT(SCREEN,ANGLE,IMPACT)   !graphics output
      IF (GFILE) CALL GRFOUT(FILE,ANGLE,IMPACT)
      IF (GHRDCP) CALL GRFOUT(PAPER,ANGLE,IMPACT)
C
      RETURN
      END
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE THETA(B,INT1,INT2,RMIN)
C  searches to find RMIN (closest approach)
C  integrates FNI1 from B    to RMAX to find INT1
C  integrates FN12 from RMIN to RMAX to find INT2
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C Global variables:
      INCLUDE 'PARAM.P1'
C Input/Output variables:
      REAL B                     !impact parameter (input)
      REAL INT1, INT2            !values for the two integrals (output)
      REAL RMIN                  !radius of closes approach (output)
C Local variables:
      REAL R                     !current radius
      REAL DR                    !current search step size
      REAL U                     !U=SQRT(R-B) or SQRT(R-RMIN)
      REAL UMAX                  !UMAX=SQRT(RMAX-B) or SQRT(RMAX-RMIN)
      REAL H                     !radial step for quadrature
      REAL SUM                   !sum for quadrature
      INTEGER IU                 !index on U
      REAL FNV                   !scattering potential (function)
      REAL FNI1,FNI2             !first and second integrands (function)
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     define functions
      FNV(R)      = 4*(R**(-12)-R**(-6))
      FNI1(B,R)   = 1/R**2/SQRT(1-(B/R)**2)
      FNI2(B,R,E) = 1/R**2/SQRT(1-(B/R)**2-FNV(R)/E)
C
      RMIN=RMAX                      !inward search for turning point
      DR = 0.2
1     IF (DR.GT.TOLR) THEN
         RMIN=RMIN-DR                !simple search
         IF ((1-((B/RMIN)**2)-(FNV(RMIN)/E)).LT.0.0) THEN
           RMIN=RMIN +DR             !RMIN is where FNI2 is infinite
           DR=DR/2                   !i.e., the denominator=0
         END IF
      GO TO 1
      END IF
C
      SUM=0                          !first integral by rectangle rule
      UMAX=SQRT(RMAX-B)              !change variable to U=SQRT(R-B)
      H=UMAX/NPTS                    !to remove singularity at R=B
      DO 10 IU=1,NPTS
            U=H*(IU-0.5)
            R=U**2+B                 !change back to R to eval integrand
            SUM=SUM+U*FNI1(B,R)
   10 CONTINUE
      INT1=2*H*SUM
C
      SUM=0                          !second integral by rectangle rule
      UMAX=SQRT(RMAX-RMIN)           !change of variable to
      H=UMAX/NPTS                    !to remove singularity at R=RMIN
      DO 22 IU=1,NPTS
            U=H*(IU-0.5)
            R=U**2+RMIN              !change back to R to eval integrand
            SUM=SUM+U*FNI2(B,R,E)
   22 CONTINUE
      INT2=2*H*SUM
C
      RETURN
      END
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE INIT
C initializes constants, displays header screen,
C initializes arrays for input parameters
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C Global variables:
      INCLUDE 'IO.ALL'
      INCLUDE 'MENU.ALL'
      INCLUDE 'PARAM.P1'
C Local parameters:
      CHARACTER*80 DESCRP        !program description
      DIMENSION DESCRP(22)
      INTEGER NHEAD,NTEXT,NGRF   !number of lines for each description
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      CALL SETUP                !get environment parameters
C
      !display header screen
      DESCRP(1)= 'PROJECT 1'
      DESCRP(2)= 'Scattering by a 6-12 potential'
      NHEAD=2
C
C     text output description
      DESCRP(3)= 'impact parameter, radius of closest approach, '
      DESCRP(4)= 'and the angle of deflection'
      NTEXT=2
C
C     graphics output description
      DESCRP(5)= 'deflection function'
     +   //' (scattering angle vs. impact parameter)'
      NGRF=1
C
      CALL HEADER(DESCRP,NHEAD,NTEXT,NGRF)
C
C     calculate constants
      PI=4*ATAN(1.0)
C
      CALL MENU                          !setup constant part of menu
C
      MTYPE(13)=FLOAT
      MPRMPT(13)= 'Enter incident energy E (scaled units)'
      MTAG(13)= 'Incident energy (scaled units)'
      MLOLIM(13)=.01
      MHILIM(13)=1000.
      MREALS(13)=1.0
C
      MTYPE(14)=FLOAT
      MPRMPT(14)= 'Enter the minimum impact parameter (scaled units)'
      MTAG(14)= 'Minimum impact parameter Bmin (scaled units)'
      MLOLIM(14)=.0
      MHILIM(14)=2.5
      MREALS(14)=0.1
C
      MTYPE(15)=FLOAT
      MPRMPT(15)= 'Enter the maximum impact parameter (scaled units)'
      MTAG(15)= 'Maximum impact parameter Bmax (scaled units)'
      MLOLIM(15)=.0
      MHILIM(15)=7.
      MREALS(15)=2.4
C
      MTYPE(16)=NUM
      MPRMPT(16)= 'Enter number of values for impact parameter'
      MTAG(16)= 'Number of values for impact parameter'
      MLOLIM(16)=1.
      MHILIM(16)=MAX
      MINTS(16)=20
C
      MTYPE(17)=SKIP
      MREALS(17)=35.
C
      MTYPE(38)=FLOAT
      MPRMPT(38)= 'Enter tolerance for turning point search'
      MTAG(38)= 'Turning point search tolerance'
      MLOLIM(38)=.000005
      MHILIM(38)=.1
      MREALS(38)=.0001
C
      MTYPE(39)=FLOAT
      MPRMPT(39)= 'Enter radius at which V(r)<<E (scaled units)'
      MTAG(39)= 'Maximum radius RMAX (scaled units)'
      MLOLIM(39)=0.
      MHILIM(39)=10.0
      MREALS(39)=2.5
C
      MTYPE(40)=NUM
      MPRMPT(40)= 'Enter number of quadrature points'
      MTAG(40)= 'Number of quadrature points'
      MLOLIM(40)=20.
      MHILIM(40)=300.
      MINTS(40)=40
C
      MTYPE(41)=SKIP
      MREALS(41)=60.
C
      MSTRNG(MINTS(75))= 'proj1.txt'
C
      MTYPE(76)=SKIP
      MREALS(76)=80.
C
      MSTRNG(MINTS(86))= 'proj1.grf'
C
      MTYPE(87)=SKIP
      MREALS(87)=90.
C
      RETURN
      END
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE PARAM
C gets parameters from screen
C ends program on request
C closes old files
C maps menu variables to program variables
C opens new files
C calculates all derivative parameters
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C Global variables:
      INCLUDE 'MENU.ALL'
      INCLUDE 'IO.ALL'
      INCLUDE 'PARAM.P1'
      INCLUDE 'MAP.P1'
C Function:
      LOGICAL LOGCVT          !converts 1 and 0 to true and false
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     get input from terminal
      CALL CLEAR
      CALL ASK(1,ISTOP)
C
C     stop program if requested
      IF (MREALS(IMAIN) .EQ. STOP) CALL FINISH
C
C     close files if necessary
      IF (TNAME .NE. MSTRNG(MINTS(ITNAME)))
     +     CALL FLCLOS(TNAME,TUNIT)
      IF (GNAME .NE. MSTRNG(MINTS(IGNAME)))
     +     CALL FLCLOS(GNAME,GUNIT)
C
C     physical and numerical parameters
      E=MREALS(IE)
      BMIN=MREALS(IBMIN)
      BMAX=MREALS(IBMAX)
      NB=MINTS(INB)
      RMAX=MREALS(IRMAX)
      TOLR=MREALS(ITOLR)
      NPTS=MINTS(INPTS)
C
C     text output
      TTERM=LOGCVT(MINTS(ITTERM))
      TFILE=LOGCVT(MINTS(ITFILE))
      TNAME=MSTRNG(MINTS(ITNAME))
C
C     graphics output
      GTERM=LOGCVT(MINTS(IGTERM))
      GHRDCP=LOGCVT(MINTS(IGHRD))
      GFILE=LOGCVT(MINTS(IGFILE))
      GNAME=MSTRNG(MINTS(IGNAME))
C
C     open files
      IF (TFILE) CALL FLOPEN(TNAME,TUNIT)
      IF (GFILE) CALL FLOPEN(GNAME,GUNIT)
C     files may have been renamed
      MSTRNG(MINTS(ITNAME))=TNAME
      MSTRNG(MINTS(IGNAME))=GNAME
C
      CALL PCHECK            !check 0<BMIN<BMAX<=RMAX
      DB=(BMAX-BMIN)/(NB-1)  !calculate step in B
      CALL CLEAR
C
      RETURN
      END
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE PCHECK
C to ensure that 0<BMIN<BMAX<=RMAX
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C Global variables:
      INCLUDE 'MENU.ALL'
      INCLUDE 'IO.ALL'
      INCLUDE 'PARAM.P1'
      INCLUDE 'MAP.P1'
C Function:
      REAL GETFLT                      !gets floating point number
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
40    IF (BMAX .LE. BMIN) THEN
         WRITE (OUNIT,10) BMAX,BMIN
         WRITE (OUNIT,20)
         MREALS(IBMIN)=GETFLT(MREALS(IBMIN),MLOLIM(IBMIN),MHILIM(IBMIN),
     +                 'Enter BMIN')
         MREALS(IBMAX)=GETFLT(MREALS(IBMAX),BMIN,MHILIM(IBMAX),
     +                 'Enter BMAX')
         BMAX=MREALS(IBMAX)
         BMIN=MREALS(IBMIN)
      GO TO 40
      END IF
C
50    IF (RMAX .LE. BMAX) THEN
         WRITE (OUNIT,30) RMAX,BMAX
         WRITE (OUNIT,20)
         !prompt for new values
         MREALS(IBMAX)=GETFLT(MREALS(IBMAX),BMIN,MHILIM(IBMAX),
     +                 'Enter BMAX')
         BMAX=MREALS(IBMAX)
         MREALS(IRMAX)=GETFLT(MREALS(IRMAX),MLOLIM(IRMAX),MHILIM(IRMAX),
     +                 'Enter RMAX')
         RMAX=MREALS(IRMAX)
      GO TO 50
      END IF
C
10    FORMAT (/,' BMAX (=',F6.3,')',' is less than BMIN (=',
     +        F6.3,'),')
20    FORMAT ('  but it should be the other way round; '
     +         ' reenter parameter values')
30    FORMAT (/,' RMAX (=',F6.3,')',' is less than BMAX (=',F6.3,')')
C
      RETURN
      END
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE PRMOUT(MUNIT,NLINES)
C outputs parameter summary to the specified unit
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C Global variables
      INCLUDE 'IO.ALL'
      INCLUDE 'PARAM.P1'
C Passed variables
      INTEGER MUNIT             !unit number for output (input)
      INTEGER NLINES            !number of lines written so far (output)
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      IF (MUNIT .EQ. OUNIT) CALL CLEAR
C
      WRITE (MUNIT,19)
      WRITE (MUNIT,21)
      WRITE (MUNIT,27)  E
      WRITE (MUNIT,23)  RMAX,NPTS
      WRITE (MUNIT,25)  TOLR
      WRITE (MUNIT,29)  BMIN,BMAX,NB
      WRITE (MUNIT,30)
      WRITE (MUNIT,19)
C
C     different data headers for graphics and text output
      IF (MUNIT .NE. GUNIT) THEN
        WRITE (MUNIT,31)
        WRITE (MUNIT,32)
        WRITE (MUNIT,33)
      ELSE
        WRITE (MUNIT,35)
      END IF
C
      NLINES = 11
C
19    FORMAT  (' ')
21    FORMAT  (' Output from project 1:',
     +         '  Scattering by a central 6-12 potential ')
23    FORMAT  (' Rmax=',F6.3,5X,'number of quadrature points =',I4)
25    FORMAT  (' Turning point tolerance = ',E12.5)
27    FORMAT  (' Incident energy = 'F8.3)
29    FORMAT  (' Bmin = ',F6.3,5X,'Bmax = ',F6.3,
     +         5X,'number of B values = ',I3)
30    FORMAT  (' Lengths are in scaled units, angles in degrees')
31    FORMAT (15X,'Impact',14X,'Closest',16X,'Scattering')
32    FORMAT (13X,'Parameter',13X,'Approach'17X,'Angle')
33    FORMAT (13X,'---------',13X,'--------'17X,'-----')
35    FORMAT  (2X,'Impact parameter',2X,'Scattering Angle')
C
      RETURN
      END
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE TXTOUT(MUNIT,NLINES,B,RMIN,ANGLE)
C writes results for one impact parameter to the requested unit
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C Global variables:
      INCLUDE 'IO.ALL'
C Passed variables:
      INTEGER MUNIT                !output unit specifier (input)
      INTEGER NLINES               !number of lines printed so far (I/0)
      REAL B                       !impact parameter (input)
      REAL RMIN                    !radius of closest approach (input)
      REAL ANGLE                   !scattering angle(input)
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     if screen is full, clear screen and retype headings
      IF ((MOD(NLINES,TRMLIN-4) .EQ. 0)
     +                          .AND. (MUNIT .EQ. OUNIT)) THEN
         CALL PAUSE('to continue...',1)
         CALL CLEAR
         WRITE (MUNIT,31)
         WRITE (MUNIT,32)
         WRITE (MUNIT,33)
         NLINES=3
      END IF
C
      WRITE (MUNIT,35) B,RMIN,ANGLE
C
C     keep track of printed lines only for terminal output
      IF (MUNIT .EQ. OUNIT) NLINES=NLINES+1
C
31    FORMAT (15X,'Impact',14X,'Closest',16X,'Scattering')
32    FORMAT (13X,'Parameter',13X,'Approach'17X,'Angle')
33    FORMAT (13X,'---------',13X,'--------'17X,'-----')
35    FORMAT  (15X,F6.3,15X,F6.3,15X,F9.3)
C
      RETURN
      END
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE GRFOUT(DEVICE,ANGLE,IMPACT)
C graphs scattering angle vs. impact parameter
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C Global variables
      INCLUDE 'IO.ALL'
      INCLUDE 'PARAM.P1'
      INCLUDE 'GRFDAT.ALL'
C Input variables:
      INTEGER DEVICE            !which device is being used?
      REAL ANGLE(0:MAX)         !final scattering angle
      REAL IMPACT(0:MAX)        !impact parameter
C Local variables
      INTEGER IB                !index on impact parameter
      CHARACTER*9 CE,CR         !energy, RMAX as a character string
      INTEGER LEN,RLEN          !length of string
      INTEGER SCREEN            !send to terminal
      INTEGER PAPER             !make a hardcopy
      INTEGER FILE              !send to a file
      DATA SCREEN,PAPER,FILE/1,2,3/
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     messages for the impatient
      IF (DEVICE .NE. SCREEN) WRITE (OUNIT,100)
C
C     calculate parameters for graphing
      IF (DEVICE .NE. FILE) THEN
C
          NPLOT=1                        !how many plots?
          IPLOT=1
C
          YMAX=0.                        !find limits on data points
          YMIN=0.
          DO 20 IB=0,NB-1
             IF (ANGLE(IB) .GT. YMAX) YMAX=ANGLE(IB)
             IF (ANGLE(IB) .LT. YMIN) YMIN=ANGLE(IB)
20        CONTINUE
          XMIN=BMIN
          XMAX=BMAX
          Y0VAL=BMIN
          X0VAL=0.
C
          NPOINT=NB
C
          ILINE=1                      !line and symbol styles
          ISYM=1
          IFREQ=1
          NXTICK=5
          NYTICK=5
C
          CALL CONVRT(E,CE,LEN)             !titles and labels
          CALL CONVRT(RMAX,CR,RLEN)
          INFO='RMAX = '//CR(1:RLEN)
          TITLE = 'Scattering from 6-12 potential, Energy='//CE
          LABEL(1)= 'Impact parameter (scaled units)'
          LABEL(2)= 'Scattering angle (degrees)'
C
          CALL GTDEV(DEVICE)                   !device nomination
          IF (DEVICE .EQ. SCREEN) CALL GMODE   !change to graphics mode
          CALL LNLNAX                          !draw axes
      END IF
C
C     output results
      IF (DEVICE .EQ. FILE) THEN
          WRITE (GUNIT,70) (IMPACT(IB),ANGLE(IB),IB=1,NB)
      ELSE
          CALL XYPLOT (IMPACT,ANGLE)
      END IF
C
C     end graphing session
      IF (DEVICE .NE. FILE) CALL GPAGE(DEVICE)  !close graphing package
      IF (DEVICE .EQ. SCREEN) CALL TMODE        !switch to text mode
C
70    FORMAT (2(5X,E11.3))
100   FORMAT (/,'Patience, please; output going to a file.')
C
      RETURN
      END

ext: for


! Have to enable optimizations, as otherwise builtins won't be expanded.
! { dg-additional-options "-O -fdump-rtl-expand" }

module openacc_kinds
  implicit none

  integer, parameter :: acc_device_kind = 4

end module openacc_kinds

module openacc
  use openacc_kinds
  implicit none

  integer (acc_device_kind), parameter :: acc_device_host = 2

  interface
     function acc_on_device (dev)
       use openacc_kinds
       logical (4) :: acc_on_device
       integer (acc_device_kind), intent (in) :: dev
     end function acc_on_device
  end interface
end module openacc

logical (4) function f ()
  use openacc
  implicit none

  integer (4), parameter :: dev = 2

  f = acc_on_device (dev)
  return
end function f

! With -fopenacc, we're expecting the builtin to be expanded, so no calls.
! TODO: not working.
! { dg-final { scan-rtl-dump-times "\\\(call \[^\\n\]* acc_on_device" 0 "expand" { xfail *-*-* } } }


text:
!------------------------------------------------------------------------------!
!  The Community Multiscale Air Quality (CMAQ) system software is in           !
!  continuous development by various groups and is based on information        !
!  from these groups: Federal Government employees, contractors working        !
!  within a United States Government contract, and non-Federal sources         !
!  including research institutions.  These groups give the Government          !
!  permission to use, prepare derivative works of, and distribute copies       !
!  of their work in the CMAQ system to the public and to permit others         !
!  to do so.  The United States Environmental Protection Agency                !
!  therefore grants similar permission to use the CMAQ system software,        !
!  but users are requested to provide copies of derivative works or            !
!  products designed to operate in the CMAQ system to the United States        !
!  Government without restrictions as to use by others.  Software              !
!  that is used with the CMAQ system but distributed under the GNU             !
!  General Public License or the GNU Lesser General Public License is          !
!  subject to their copyright restrictions.                                    !
!------------------------------------------------------------------------------!

SUBROUTINE gridout (sdate, stime)

!-------------------------------------------------------------------------------
! Name:     GRID Output -- Create "grid" (or time-independent) output for CTM
! Purpose:  Output time-independent fields.
! Revised:  18 Dec 2018  Original version in MCIPv5.0.  Subsumes part of
!                        gridout.f90 from MCIPv4.5.  (T. Spero)
!-------------------------------------------------------------------------------

  USE mcipparm, ONLY: ioform

  IMPLICIT NONE

  INTEGER,            INTENT(IN)    :: sdate
  INTEGER,            INTENT(IN)    :: stime

!-------------------------------------------------------------------------------
! Write time-independent output fields.
!-------------------------------------------------------------------------------

  SELECT CASE ( ioform )

    CASE ( 1 )  ! Models-3 I/O API
      CALL outgm3io (sdate, stime)

  END SELECT

!-------------------------------------------------------------------------------
! Write GRIDDESC file.
!-------------------------------------------------------------------------------

  CALL wrgdesc

!-------------------------------------------------------------------------------
! Print sample output to log file.
!-------------------------------------------------------------------------------

  CALL outglog

END SUBROUTINE gridout

max_stars_repo_name: Simeng-unique/CMAQ-changed

max_issues_repo_name: Simeng-unique/CMAQ-changed

max_forks_repo_name: Simeng-unique/CMAQ-changed


text:
SUBROUTINE clawpack5_rpt2adv(ixy,imp,maxm,meqn,mwaves,maux,&
     mbc,mx,ql,qr,aux1,aux2,aux3,asdq,bmasdq,bpasdq)
  IMPLICIT NONE

  INTEGER ixy,imp,maxm,meqn,mwaves,maux,mbc,mx

  DOUBLE PRECISION     ql(meqn, 1-mbc:maxm+mbc)
  DOUBLE PRECISION     qr(meqn, 1-mbc:maxm+mbc)
  DOUBLE PRECISION   asdq(meqn, 1-mbc:maxm+mbc)
  DOUBLE PRECISION bmasdq(meqn, 1-mbc:maxm+mbc)
  DOUBLE PRECISION bpasdq(meqn, 1-mbc:maxm+mbc)
  DOUBLE PRECISION   aux1(maux, 1-mbc:maxm+mbc)
  DOUBLE PRECISION   aux2(maux, 1-mbc:maxm+mbc)
  DOUBLE PRECISION   aux3(maux, 1-mbc:maxm+mbc)

  INTEGER iface, i, i1, mq

  iface = 3-ixy  !#  = 1 if ixy=2  or  = 2 if ixy=1
  DO i = 2-mbc,mx+mbc
     i1 = i-2+imp    !#  =  i-1 for amdq,  i for apdq
     do mq = 1,meqn
        bmasdq(mq,i) = min(aux2(iface,i1), 0.d0) * asdq(mq,i)
        bpasdq(mq,i) = max(aux3(iface,i1), 0.d0) * asdq(mq,i)
     end do
  END DO

  RETURN
END SUBROUTINE clawpack5_rpt2adv

max_stars_repo_name: ECLAIRWaveS/ForestClaw

max_issues_repo_name: ECLAIRWaveS/ForestClaw

max_forks_repo_name: ECLAIRWaveS/ForestClaw


text:
      SUBROUTINE allocate_modular_coils
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE boundary, ONLY: nfp
      USE modular_coils
      USE Vcoilpts
      USE Vwire
      USE coils
      USE mpi_params
      IMPLICIT NONE
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      INTEGER :: nc, i, status, modes
      LOGICAL :: lmaster
!-----------------------------------------------
      lmaster = (myid .EQ. master)
!     ALLOCATE modular coil arrays

      nc = nmod_coils_per_period
      IF (lmaster .AND. nc.LE.0)
     1  STOP 'nmod_coils_per_period must > 0!'

      ALLOCATE (modular(nc), stat=status)
      ALLOCATE (x_mod(nwdim1,3,ncdim), y_mod(nwdim1,3,ncdim),
     1          z_mod(nwdim1,3,ncdim), rho(ncdim,nwdim1) ,
     2          phi(ncdim,nwdim1), rcoil(ncdim,nwdim),
     3          zcoil(ncdim,nwdim), stat=status)

      IF (lmaster .AND. status.NE.0)
     1   STOP "Cannot ALLOCATE modular_coils"

      nfper = nfp
      nmod_coils = nc * nfper
      IF (lmaster .AND. nmod_coils.GT.ncdim)
     1   STOP 'nmod_coils > ncdim'

      nmid = (nc+1)/2
      nodd = MOD(nc,2)
      IF (lmaster) THEN
        PRINT *, 'no. modular coils (per period) = ', nc
      END IF

      IF (nmod_coils .LE. 0) RETURN

!     Initialize arrays to values of unique coil parameters

!     First consider the CASE with coils on both symmetry planes at
!     phi = 0 and phi = pi/nfp (lsymm = F). This implies that the number
!     of coils per field period must be even (nodd = 0).

      IF ((nodd.EQ.0) .AND. (.NOT.lsymm)) THEN
        nmid = nmid + 1

        IF (lmaster) THEN
          PRINT *, 'nmid = ', nmid, 'nodd = ', nodd
        END IF

!     Allocate only those coil components which are needed. There
!     is no need to store the Fourier coefficients for coils whose
!     locations in real space will be computed from stellarator
!     symmetry.

!     ALLOCATE coil components for symmetry coil at phi = 0.
        i = 1
        ALLOCATE(modular(i)%rhoc(0:nf_rho), modular(i)%rhos(0:nf_rho),
     1           modular(i)%phic(0:nf_phi), modular(i)%phis(0:nf_phi),
     2           modular(i)%phi(nwire1),    modular(i)%rho(nwire1),
     3           stat=status)
        IF (lmaster .AND. status.NE.0)
     1     STOP "Cannot ALLOCATE center coil components"

!     ALLOCATE coil components for coils 2 through nmid-1.
        DO i = 2, nmid-1
          ALLOCATE(modular(i)%rhoc(0:nf_rho), modular(i)%rhos(0:nf_rho),
     1             modular(i)%phic(0:nf_phi), modular(i)%phis(0:nf_phi),
     2             modular(i)%phi(nwire1),    modular(i)%rho(nwire1),
     3             stat = status)
        IF (lmaster .AND. status .NE. 0)
     1     STOP "Cannot ALLOCATE modular components"

        END DO               !! DO i = ... loop over half of coils

!     ALLOCATE coil components for symmetry coil at phi = pi/nfp.
        i = nmid

        ALLOCATE(modular(i)%rhoc(0:nf_rho), modular(i)%rhos(0:nf_rho),
     1           modular(i)%phic(0:nf_phi), modular(i)%phis(0:nf_phi),
     2           modular(i)%phi(nwire1),    modular(i)%rho(nwire1),
     3           stat=status)
        IF (lmaster .AND. status.NE.0)
     1     STOP "Cannot ALLOCATE center coil components"

!     ALLOCATE coil components phi, rhi for coils nmid + 1 through nc.
        DO i = 1, nmid-2

          ALLOCATE(modular(nc-i+1)%phi(nwire1),
     1             modular(nc-i+1)%rho(nwire1),
     2             stat=status)
          IF (lmaster .AND. status.NE.0)
     1       STOP "Cannot ALLOCATE modular components"

        END DO               !! DO i = ... loop over half of coils

!     Symmetry coil at phi = 0.
        i = 1
        modes = 0
        phic(i,modes) = 0
        modular(i)%phic(modes) = 0
        modular(i)%phis(modes) = phis(i,modes)
        DO modes = 1, nf_phi
          phic(i,modes) = 0
          modular(i)%phic(modes) = 0
          modular(i)%phis(modes) = phis(i,modes)
        END DO

        modes = 0
        DO modes = 0, nf_rho
          modular(i)%rhos(modes) = rhos(i,modes)
          modular(i)%rhoc(modes) = rhoc(i,modes)
        END DO
        modular(i)%rhos(modes) = 0

!     Coils 2 through nmid-1.
        DO i = 2, nmid-1
          modes = 0
          modular(i)%phic(modes) = phic(i,modes)
          modular(i)%phis(modes) = 0
          DO modes = 1,nf_phi
            modular(i)%phic(modes) = phic(i,modes)
            modular(i)%phis(modes) = phis(i,modes)
          END DO

          modes = 0
          DO modes = 0,nf_rho
            modular(i)%rhos(modes) = rhos(i,modes)
            modular(i)%rhoc(modes) = rhoc(i,modes)
          END DO
          modular(i)%rhos(modes) = 0
        END DO

!     Symmetry coil at phi = pi/nfp.
        i = nmid
        modes = 0
        phic(i,modes) = 0
        modular(i)%phic(modes) = 0
        modular(i)%phis(modes) = phis(i,modes)
        DO modes = 1, nf_phi
          phic(i,modes) = 0
          modular(i)%phic(modes) = 0
          modular(i)%phis(modes) = phis(i,modes)
        END DO

        modes = 0
        DO modes = 1, nf_rho
          modular(i)%rhos(modes) = rhos(i,modes)
          modular(i)%rhoc(modes) = rhoc(i,modes)
        END DO
        modular(i)%rhos(modes) = 0

!     Set currents from input

        DO i = 1, nmid
          modular(i)%current = curmod(i)
        END DO

!     Impose symmetry on coil-coil penalty weights, exponents

        DO i = 1, nmid-2
          dcc_wgt(nc+1-i) = dcc_wgt(i+1)
          dcc_exp(nc+1-i) = dcc_exp(i+1)
          dcc_tgt(nc+1-i) = dcc_tgt(i+1)
          rc_wgt(nc+1-i) = rc_wgt(i+1)
          rc_exp(nc+1-i) = rc_exp(i+1)
          rc_tgt(nc+1-i) = rc_tgt(i+1)
          r_ext(nc+1-i) = r_ext(i+1)
        END DO

!     Next consider the cases with a coil on phi = 0 (lsymm = t), or
!     on phi = pi/nfp (lsymm = F), but not both. Then there may be an
!     even number of coils per period (nodd = 0) or an odd number of
!     coils per period (nodd = 1).

      ELSE

        IF (lmaster) THEN
          PRINT *, 'nmid = ', nmid
          PRINT *, 'nodd = ', nodd
        END IF

        DO i = 1,nmid-nodd

!     Allocate only those coil components which are needed. There
!     is no need to store the Fourier coefficients for coils whose
!     locations in real space will be computed from stellarator
!     symmetry.

          ALLOCATE(modular(i)%rhoc(0:nf_rho), modular(i)%rhos(0:nf_rho),
     1             modular(i)%phic(0:nf_phi), modular(i)%phis(0:nf_phi),
     2             modular(i)%phi(nwire1),    modular(i)%rho(nwire1),
     3             modular(nmod_coils_per_period-i+1)%phi(nwire1),
     4             modular(nmod_coils_per_period-i+1)%rho(nwire1),
     5             stat=status)
          IF (lmaster .AND. status.NE.0)
     1       STOP "Cannot ALLOCATE modular components"

        END DO               !! DO i = ... loop over half of coils

        IF (nodd .eq. 1) THEN
          i = nmid

          ALLOCATE(modular(i)%rhoc(0:nf_rho), modular(i)%rhos(0:nf_rho),
     1             modular(i)%phic(0:nf_phi), modular(i)%phis(0:nf_phi),
     2             modular(i)%phi(nwire1),    modular(i)%rho(nwire1),
     3             stat = status)
          IF (lmaster .AND. status.NE.0)
     1       STOP "Cannot ALLOCATE center coil components"
        ENDIF

        DO i = 1, nmid-nodd
          modes = 0
          modular(i)%phic(modes) = phic(i,modes)
          modular(i)%phis(modes) = 0
          DO modes = 1,nf_phi
            modular(i)%phic(modes) = phic(i,modes)
            modular(i)%phis(modes) = phis(i,modes)
          END DO

          modes = 0
          modular(i)%rhos(modes) = 0
          DO modes = 1,nf_rho
            modular(i)%rhos(modes) = rhos(i,modes)
          END DO
        END DO

        IF (nodd .EQ. 1) THEN
          i = nmid
          modes = 0
          phic(i,modes) = 0
          modular(i)%phic(modes) = 0
          modular(i)%phis(modes) = phis(i,modes)
          DO modes = 1, nf_phi
            phic(i,modes) = 0
            modular(i)%phic(modes) = 0
            modular(i)%phis(modes) = phis(i,modes)
          END DO

          modes = 0
          modular(i)%rhos(modes) = 0
          DO modes = 1, nf_rho
            modular(i)%rhos(modes) = rhos(i,modes)
          END DO
        END IF

!     Set currents from input

        DO i = 1, nmid
          modular(i)%current = curmod(i)
        END DO

!     Impose symmetry on coil-coil penalty weights, exponents

        DO i = 1,nmid-nodd
          dcc_wgt(nc+1-i) = dcc_wgt(i)
          dcc_exp(nc+1-i) = dcc_exp(i)
          dcc_tgt(nc+1-i) = dcc_tgt(i)
          rc_wgt(nc+1-i) = rc_wgt(i)
          rc_exp(nc+1-i) = rc_exp(i)
          rc_tgt(nc+1-i) = rc_tgt(i)
          r_ext(nc+1-i) = r_ext(i)
        END DO

      END IF  ! END IF ((nodd .eq. 0) .and. (lsymm .eqv. .false.))

      nmod_unique_coils = nmid

      END SUBROUTINE allocate_modular_coils

ext: f

max_stars_repo_name: joseluisvelasco/STELLOPT

max_issues_repo_name: joseluisvelasco/STELLOPT

max_forks_repo_name: joseluisvelasco/STELLOPT


text:
   module thermoData
      ! This thermodynamic data is from Alberty2003 and his Mathematica
      ! BiochemThermo`BasicBiochemData3.m data.
      ! For bacterial structure, have used data on yeast from Battley1998.
      ! Energy units are in kJ/mole
      ! V3: This version is identical to files labeled with V2.  Not sure why this was given that different name
      !    This file was copied from: C:\Users\jvallino\Projects\PIE-LTER\Modeling\PIE-1D-MEP\MEP-phototroph-1D-OC
      ! V4, 17-Aug-2016: Giving this a new version, just to distinguish it from V3.  A few changes were made
      !    Made R a global module parameter, and changed name to RJ to distinguish it from other R values.
      !    Used NIST values for R (also used Wikipedia's calculations for others.)
      !    Replaced data statements with parameters statements
      ! V4.1, 28-Aug-2016
      !     Added hydrogen (aq) species to parameter list
      ! V4.2, 20-Apr-2017: Nothing changed, but did add some new comments and want to keep those as versions go forward.
      !     No, now it has changed.  There is a slight descrepancy in Alberty2003 tables 8.1 and 8.2.  The values for co2aqsp have
      !     been updated here.  Also, species for H2CO3, HCO3- and CO3-- have been added so that the fraction of h2co3, hco3 and co3 can be extracted
      !     from co2tot.  See freeHCO3 and freeCO3 functions.  It has been varified that freeCO2 + freeHCO3 + freeCO3 = 1 (see CarbonateChemistry2.nb in Mma directory)
      !     Also, added free energy of formation for aveAA, aveNA, SiO2 and Si(OH)4
      ! V4.3, 21-Jul-2017: Added function to calculate H2S in equilibrium between air and water.
      ! V4.4, 17-May-2017: Added
      !     glucosesp, avePLsp, aveChlsp
      ! V5.0, 24-Jul-2018:
      !     Replaced parameter attribute of GHdata with target attribute
      !     Introduced derived type GHdb and array GHznData that points string names to sp data
      !     Introduced new function dGf0n(name,T,is,pH) that can use a string name for compound
      !     Introduced new routine initThermo(), that must be called to populate GHznData needed by dGf0n
      !     Introduced new parameters maxSpec, maxName and maxDim (see below).
      !     added compounds: acetatesp, ethanolsp, xylosesp
      !     Added gas constant for bar, Rbar
      ! V5.1, 22-Jan-2019:
      !     Added Mg(OH)2 (aq) and Mg2+ (aq) to the database for Chl a synthesis reaction.
      implicit none
      save ! Probably not necessary
      integer, parameter:: maxSpec = 5 ! This is the maximum number of species allowed for compound
      integer, parameter:: maxName = 10 ! This is the maximum length allowed for species names
      integer, parameter:: maxDim = 38  ! maximum allowed lenght of GHznData that holds the character names and pointer assignment
      type GHdata
         real(8) dGf ! Standar free energy of formation, kJ/mol (25 C, 1bar)
         real(8) dHf ! Standard enthalpy of formation kJ/mol (25 C, 1 bar)
         real(8) zi  ! Charge on molecule
         real(8) nH  ! Number of hydrogen atoms at the given charge
      end type GHdata
      type GHdb
         ! This is introduced so that character names can be used to refer to species
         character(maxName) name
         type(GHdata), pointer:: sp(:)
      end type GHdb

      type(GHdb) GHznData(maxDim) ! This is needed to use string names in function dGf0n

      real(8), parameter:: Rp   = 0.082057338d0  ! gas constant (atm L/(g-mol K))
      real(8), parameter:: Rpm3 = 0.082057338d-6 ! gas constant (atm m^3/(g-mmol K)) (note, mmol not mol)
      real(8), parameter:: RkJ  = 8.3144598d-3   ! gas constant (kJ/(g-mol K) or J/(g-mmol K))
      real(8), parameter:: RJ   = 8.3144598d0    ! (J/g-mol/K) or (Pa m^3/(K g-mol))
      real(8), parameter:: Rbar = 8.3144598d-5   ! (bar m^3/K/g-mol) or (mbar m^3/K/g-mmol)

      ! Note, in order for all quantities to be treated as arrays, it is necessary to dimension species
      ! even if they only have one speces, such as h2osp(1).
      type(GHdata), target:: atpsp(3)		   = [GHdata(-2768.1,-3619.21,-4,12),GHdata(-2811.48,-3612.91,-3,13),GHdata(-2838.18,-3627.91,-2,14)]
      type(GHdata), target:: methaneaqsp(1)  = [GHdata(-34.33,-89.04,0,4)]
      type(GHdata), target:: methanegsp(1)	= [GHdata(-50.72,-74.81,0,4)]
      type(GHdata), target:: methanolsp(1)	= [GHdata(-175.31,-245.93,0,4)]
      type(GHdata), target:: methanolgsp(1)  = [GHdata(-162.3,-201.0,0,4)] ! gas phase data from CRC Handbook of Chemistry and Physics, Vol 91.
      type(GHdata), target:: adpsp(3)		   = [GHdata(-1906.13,-2626.54,-3,12),GHdata(-1947.1,-2620.94,-2,13),GHdata(-1971.98,-2638.54,-1,14)]
      type(GHdata), target:: pisp(2)		   = [GHdata(-1096.1,-1299.,-2,1),GHdata(-1137.3,-1302.6,-1,2)] ! phosphate h3po4
      type(GHdata), target:: h2osp(1)		   = [GHdata(-237.19,-285.83,0,2)]
      type(GHdata), target:: o2gsp(1)		   = [GHdata(0.,0.,0.,0.)]
      type(GHdata), target:: o2aqsp(1)		   = [GHdata(16.4,-11.7,0,0)]
      type(GHdata), target:: yeastsp(1)      = [GHdata(-87.96,-133.09,0,1.613)] !Battley1998: CH1.613 O0.557 N0.158 P0.012 S0.003 K0.022 Mg0.003 Ca0.001
      type(GHdata), target:: ch2osp(1)		   = [GHdata(-152.65,-210.365,0,2)] ! this is glucose, but on unit C basis (it will be 6 times less than glucose)
      type(GHdata), target:: nadoxsp(1)		= [GHdata(0,0,-1,26)]
      type(GHdata), target:: nadredsp(1)	   = [GHdata(22.65,-31.94,-2,27)]
      type(GHdata), target:: sulfatesp(2)	   = [GHdata(-744.53,-909.27,-2,0),GHdata(-755.91,-887.34,-1,1)]
      type(GHdata), target:: ammoniasp(2)    = [GHdata(-26.5,-80.29,0,3),GHdata(-79.31,-132.51,1,4)]
      type(GHdata), target:: co2gsp(1)		   = [GHdata(-394.36,-393.5,0,0)]
      type(GHdata), target:: co2aqsp(1)      = [GHdata(-385.92,-413.80,0,0)] ! [GHdata(-385.97,-413.81,0,0)] ! this is not in the alberty database, but is given in Table 8.2, pg 151, alberty2003
                                                          ! This values are slighty different so that h2co3(sp) equals co2(sp) + h2o
      type(GHdata), target:: h2co3sp(1)      = [GHdata(-623.11,-699.63,0,2)]
      type(GHdata), target:: hco3sp(1)	      = [GHdata(-586.77,-691.99,-1,1)] ! This is only the bicarbonate species (HCO3-)
      type(GHdata), target:: co3sp(1)	      = [GHdata(-527.81,-677.14,-2,0)] ! this is only the carbonate species (CO3--)
      type(GHdata), target:: co2totsp(3)	   = [GHdata(-527.81,-677.14,-2,0),GHdata(-586.77,-691.99,-1,1),GHdata(-623.11,-699.63,0,2)]
      type(GHdata), target:: nitratesp(2)    = [GHdata(-108.74,-205.,-1,0),GHdata(-111.25,-207.36,0,1)]
      type(GHdata), target:: nitritesp(2)	   = [GHdata(-32.2,-104.6,-1,0), GHdata(-50.6,-119.2,0,1)]
      type(GHdata), target:: n2aqsp(1)		   = [GHdata(18.7,-10.54,0,0)]
      type(GHdata), target:: h2saqsp(3)      = [GHdata(85.8,33.1,-2,0),GHdata(12.08,-17.6,-1,1),GHdata(-27.83,-39.7,0,2)]
      ! These data for H2S in g are from Thauer1977 (for the dGf), and also
      ! from the book "Chemical Principles" by P. Atkins and L. Jones, gotten from google books for dHf.
      ! These data should also be in "The NBS tables of chemical thermodynamic properties :
      ! selected values for inorganic and C1 and C2 organic substances in SI units" by Donald D Wagman; et al.
      type(GHdata), target:: h2sgsp(1)		   = [GHdata(-33.56,-20.63,0.,2.)]
      ! New additions with V4.1
      type(GHdata), target:: h2aqsp(1)		   = [GHdata(17.6,-4.2,0.,2.)]
      ! New with version 4.2
      ! See excel file BIOMASS-2017.xlsx in C:\Users\jvallino\models\BIOENERG for calculations to derived aveAA and aveNA
      type(GHdata), target:: aveAAsp(1)      = [GHdata(-297.2, 0., 0.031, 7.06)] ! Aver AA in E. coli, with stoichiometry: C3.5 H7.06 O1.7 N Note, don't know delHf, but does not matter very much
      type(GHdata), target:: aveNAsp(1)      = [GHdata(-1332.7, 0., 0., 13.8)] ! Aver NA in E. coli, with stoichiometry: C9.6 H13.8 O7.8 N3.9 P Note, don't know delHf, but does not matter very muc
      type(GHdata), target:: sio2sp(1)       = [GHdata(-856.281, -910.700, 0., 0.)] ! From Gunnarsson2000
      type(GHdata), target:: h4sio4aqsp(1)   = [GHdata(-1309.181, -1461.722, 0., 4)] ! From Gunnarsson2000

      ! New with 4.4
      type(GHdata), target:: glucosesp(1)    = [GHdata(-915.9, -1262.19, 0, 12)]
      ! Here use Dihydrogeranylgeranyl diphosphate C20H38O7P2 (Kegg C17439) for phospholipids.
      ! These data are from the Equilibrator site (http://equilibrator.weizmann.ac.il/compound?compoundId=C17439)
      ! and they determined these values using Noor2013 group contribution methods.  No enthalphy of formation was available, so just set to zero.
      type(GHdata), target:: avePLsp(4)      = [GHdata(-1595.0, 0., 0., 38.), GHdata(-1584.9, 0., -1., 37.), GHdata(-1566.6, 0., -2.,36.), GHdata(-1524.3, 0., -3., 35.)]
      ! Chlorophyll a (C55H57O5N4Mg) is not in Equilibrator, but strangely enough, Chl b is.
      ! Consequently, just using the thermo for Chl b here, as it can't be too different.
      type(GHdata), target:: aveChlsp(3)     = [GHdata(-289.1, 0, -1, 71), GHdata(-340.2, 0, 0, 72), GHdata(-341.9, 0, 1,73)]

      ! New with 5.0.  These are from {alberty2003}
      type(GHdata), target:: acetatesp(2)    = [GHdata(-369.31, -486.01, -1, 3), GHdata(-396.45, -485.76, 0, 4)]
      type(GHdata), target:: ethanolsp(1)    = [GHdata(-181.64, -288.3, 0, 6)]
      type(GHdata), target:: xylosesp(1)     = [GHdata(-750.49, -1045.94, 0, 10)]

      ! New with 5.1.
      ! Mg(OH)2 (aq) from wikapedia, which  provided the delG, and the delH determined by addition from OH- and Mg2+
      type(GHdata), target:: mgoh2aqsp(1)    = [GHdata(-769.4, -926.8, 0, 2)]
      ! Mg2+ (aq) from Sawyer, McCarty, Parkin, Chemistry for Environmental Engineering and Science, 5th Edition, McGraw-Hill, 2003.
      type(GHdata), target:: mg2aqsp(1)      = [GHdata(-454.8, -466.85, 2, 0)]

contains
      subroutine initThermo()
         ! This must be called before a call to the function dGf0n, as this populates GHznData array
         GHznData( 1)%name = 'atp';         GHznData( 1)%sp => atpsp
         GHznData( 2)%name = 'methaneaq';   GHznData( 2)%sp => methaneaqsp
         GHznData( 3)%name = 'methaneg';    GHznData( 3)%sp => methanegsp
         GHznData( 4)%name = 'methanol';    GHznData( 4)%sp => methanolsp
         GHznData( 5)%name = 'methanolg';   GHznData( 5)%sp => methanolgsp
         GHznData( 6)%name = 'adp';         GHznData( 6)%sp => adpsp
         GHznData( 7)%name = 'pi';          GHznData( 7)%sp => pisp
         GHznData( 8)%name = 'h2o';         GHznData( 8)%sp => h2osp
         GHznData( 9)%name = 'o2g';         GHznData( 9)%sp => o2gsp
         GHznData(10)%name = 'o2aq';        GHznData(10)%sp => o2aqsp
         GHznData(11)%name = 'yeast';       GHznData(11)%sp => yeastsp
         GHznData(12)%name = 'ch2o';        GHznData(12)%sp => ch2osp
         GHznData(13)%name = 'nadox';       GHznData(13)%sp => nadoxsp
         GHznData(14)%name = 'nadred';      GHznData(14)%sp => nadredsp
         GHznData(15)%name = 'sulfate';     GHznData(15)%sp => sulfatesp
         GHznData(16)%name = 'ammonia';     GHznData(16)%sp => ammoniasp
         GHznData(17)%name = 'co2g';        GHznData(17)%sp => co2gsp
         GHznData(18)%name = 'co2aq';       GHznData(18)%sp => co2aqsp
         GHznData(19)%name = 'h2co3';       GHznData(19)%sp => h2co3sp
         GHznData(20)%name = 'hco3';        GHznData(20)%sp => hco3sp
         GHznData(21)%name = 'co3';         GHznData(21)%sp => co3sp
         GHznData(22)%name = 'co2tot';      GHznData(22)%sp => co2totsp
         GHznData(23)%name = 'nitrate';     GHznData(23)%sp => nitratesp
         GHznData(24)%name = 'nitrite';     GHznData(24)%sp => nitritesp
         GHznData(25)%name = 'n2aq';        GHznData(25)%sp => n2aqsp
         GHznData(26)%name = 'h2saq';       GHznData(26)%sp => h2saqsp
         GHznData(27)%name = 'h2sg';        GHznData(27)%sp => h2sgsp
         GHznData(28)%name = 'aveAA';       GHznData(28)%sp => aveAAsp
         GHznData(29)%name = 'aveNA';       GHznData(29)%sp => aveNAsp
         GHznData(30)%name = 'h4sio4aq';    GHznData(30)%sp => h4sio4aqsp
         GHznData(31)%name = 'glucose';     GHznData(31)%sp => glucosesp
         GHznData(32)%name = 'avePL';       GHznData(32)%sp => avePLsp
         GHznData(33)%name = 'aveChl';      GHznData(33)%sp => aveChlsp
         GHznData(34)%name = 'acetate';     GHznData(34)%sp => acetatesp
         GHznData(35)%name = 'ethanol';     GHznData(35)%sp => ethanolsp
         GHznData(36)%name = 'xylose';      GHznData(36)%sp => xylosesp
         GHznData(37)%name = 'mg(oh)2aq';   GHznData(37)%sp => mgoh2aqsp
         GHznData(38)%name = 'mg2+aq';      GHznData(38)%sp => mg2aqsp
         return
      end subroutine initThermo

      real(8) function dGf0n(name,T,is,pH)
         character*(*), intent(in):: name        ! Name of compound, lowercase
         real(8),   intent(in):: T, is, pH   ! temperature (K), ionic strenght (M) and pH

         ! local declarations
         integer lloc

         ! Get location of GHzn data in GHznData based on string name
         lloc = findloc(GHznData(:)%name, value=name, dim=1)
         if (lloc == 0) then
            write(*,'(a)') 'Error, compound "'//trim(name)//'" does not exisit in thermodynamic GHznData data'
            dGf0n = 0.d0
            return
         end if
         !dGf0n = dGf0(GHznData(lloc)%GHzn(1:GHznData(lloc)%nSpecies),T,is,pH)
         dGf0n = dGf0(GHznData(lloc)%sp,T,is,pH)
         return
      end function dGf0n

      real(8) function dGf0(sp,T,is,pH)
         ! This function returns the Gibbs energy of formation for the compound sp at temperature T (K),
         ! ionic strength, is (M), and pH in kJ/mol.  This is based on Alberty2003.
         !   sp:  contains, dGf, dHf, zi and nH at 298.15 and zero is and pH.
         !
         ! 30 Mar 2006: First written; Joe Vallino
         ! 21-Apr-2017: Allocating and deallocating eterms seems to chew up a bunch of CPU time. Instead, just use fix dimentsion
         !     But, if new items are added to GHDATA that are larger than 5, then eterms will need to be made larger
         implicit none
         type(GHdata) sp(:)
         real(8) T, is, pH

         ! Local declarations
         integer i, nsp
         real(8) gibbscoeff, dGzeroT, pHterm, istermG, gpfnsp, Tref
         !real(8), allocatable:: eterms(:)
         real(8) eterms(5)
         real(8) maxet
         real(8), parameter:: ln10 = log(10.d0)

         !**** Begin Program *****
         nsp = size(sp)
         !allocate ( eterms(nsp) )
         Tref = 298.15d0
         do i=1,nsp  ! sum for each species of the compound sp
            gibbscoeff=(9.20483d0*T)/1.D3-(1.284668d0*T**2)/1.D5+(4.95199d0*T**3)/1.D8
            dGzeroT=(sp(i)%dGf*T)/Tref+sp(i)%dHf*(1-T/Tref)
            pHterm=-sp(i)%nH*pH*RJ*T*ln10/1000.d0
            istermG=gibbscoeff*(sp(i)%zi**2-sp(i)%nH)*sqrt(is)/(1.d0+1.6d0*sqrt(is));
            eterms(i)=-1.d3*(dGzeroT-pHterm-istermG)/(RJ*T)
         end do
         ! Need to bring the exponents down to values that can be calculated
         ! without overflow (underflow ok), so subtract the largest one.
         maxet = maxval(eterms(1:nsp))
         dGf0 = 0.d0
         do i=1,nsp
            dGf0 = dGf0 + exp(eterms(i) - maxet)
         end do
         dGf0 = -RJ*T*(log(dGf0) + maxet)/1.d3
         !deallocate ( eterms )
         return
      end function dGf0

      real(8) function solCH4(T, is, pH)
        ! This function returns the solubiility of methane in water (uM/atm) at temperature T (K),
        ! pH and ionic strength is (M).  It is based on Alberty's values for CH4 in gas and aquious phase.
        ! However, these reactions are not affected by pH nor is, so could still have used NIST data,
        ! but this function provides consistency with Alberty's data and approach, so will keep.
        implicit none
        real(8) T, pH, is
        ! Local declarations
        real(8) dGrxn
        ! Use the following reaction:  CH4(g) -> CH4(aq)
        dGrxn = dGf0(methaneaqsp,T,is,pH) - dGf0(methanegsp,T,is,pH)

        solCH4 = exp(-1.d3*dGrxn/(T*RJ))*1.0d6*1.01325d0 ! Last number converts Bar to ATM

        return
      end function solCH4

      real(8) function solCH3OH(T, is, pH)
         ! This function returns the solubiility of methanol in water (uM/atm) at temperature T (K),
         ! pH and ionic strength is (M).  It is based on Alberty's values for CH4 in aquious phase
         ! and CH4 in the gas phase from CRC data.
         ! However, these reactions are not affected by pH nor is, so could still have used NIST data,
         ! but this function provides consistency with Alberty's data and approach, so will keep.
         implicit none
         real(8) T, pH, is
         ! Local declarations
         real(8) dGrxn
         ! Use the following reaction:  CH3OH(g) -> CH3OH(aq)
         dGrxn = dGf0(methanolsp,T,is,pH) - dGf0(methanolgsp,T,is,pH)

         solCH3OH = exp(-1.d3*dGrxn/(T*RJ))*1.0d6*1.01325d0 ! Last number converts Bar to ATM

         return
      end function solCH3OH

      real(8) function solH2CO3(T, is, pH)
         ! This function returns the solubiility of carbon dioxide in water (uM/atm) at temperature T (K)
         ! pH, and ionic strength is (M).  That is, given the CO2 partial pressure in ATM, this returns the
         ! total concentration of all carbonate species (ie., DIC = solH2CO3(T,is,pH)*pCO2 ).  This is from Alberty2003, pg 151.
         ! Rxn is:  H2CO3(aq) -> CO2(g) + H2O
         ! *** Note, since mass transfer coef are based on free CO2, use co2g2l to get mass transfer driver, that is, use
         ! flux = (pCO2*co2g2l(T, is, pH) - freeCO2(T, is, pH)*[DIC])
         ! NOT (pCO2*solH2CO3(T, is, pH) - [DIC])
         implicit none
         real(8) T, pH, is
         ! Local declarations
         real(8) dGrxn
         dGrxn = dGf0(co2gsp,T,is,pH) + dGf0(h2osp,T,is,pH) - dGf0(co2totsp,T,is,pH)

         solH2CO3 = exp(1.d3*dGrxn/(T*RJ))*1.0d6*1.01325d0 ! Last number converts Bar to ATM

         return
      end function solH2CO3

      real(8) function freeCO2(T, is, pH)
         ! This function returns the ratio of CO2 to totalCO2 (i.e. DIC) at temperature T (K)
         ! pH, and ionic strength is (M).  freeCO2 is unitless.  That is, given the DIC concentration, this returns the
         ! free CO2 in water ([CO2] = freeCO2*[DIC]).
         ! Rxn is:  co2totsp -> h2co3
         ! Note, this use to be written as H2CO3(aq) -> CO2(aq) + H2O and co2aqsp was used, but instead this has been replaced with that above.
         ! see below for freeHCO3 and freeCO3 for similar calculations
         implicit none
         real(8) T, pH, is
         ! Local declarations
         real(8) dG0rxn
         dG0rxn = dGf0(h2co3sp,T,is,pH) - dGf0(co2totsp,T,is,pH)
         ! now dGrxn = dG0rxn + RT ln ([h2co3]/[co2tot]), but a equilibrium, dGrxn = 0, so dG0rxn = - RT ln ([h2co3]/[co2tot])
         ! so freeCO2  = [h2co3]/[co2tot] = Exp( -dG0rxn/(RT) ), but dG0rxn is kJ/mol, so must use RkJ = RJ/1000 to get units to cancel
         freeCO2 = exp( -dG0rxn/(T*RkJ) ) ! This is unitless ([CO2]/[DIC])

         return
      end function freeCO2

      real(8) function freeHCO3(T, is, pH)
         ! This function returns the ratio of HCO3 to totalCO2 (i.e. DIC) at temperature T (K)
         ! pH, and ionic strength is (M).  freeHCO3 is unitless.  That is, given the DIC concentration, this returns the
         ! free HCO3 in water ([HCO3] = freeHCO3*[DIC]).
         ! Rxn is:  co2tot(aq) -> HCO3(-1) + H(+); however, the H(+) does not have to be added to the free energy calculations as described by Alberty
         implicit none
         real(8) T, pH, is
         ! Local declarations
         real(8) dG0rxn
         dG0rxn = dGf0(hco3sp,T,is,pH) - dGf0(co2totsp,T,is,pH)
         freeHCO3 = exp( -dG0rxn/(T*RkJ) ) ! This is unitless ([HCO3]/[DIC])
         return
      end function freeHCO3

      real(8) function freeCO3(T, is, pH)
         ! This function returns the ratio of CO3 to totalCO2 (i.e. DIC) at temperature T (K)
         ! pH, and ionic strength is (M).  freeCO3 is unitless.  That is, given the DIC concentration, this returns the
         ! free CO3 in water ([CO3] = freeHO3*[DIC]).
         ! Rxn is:  co2tot(aq) -> CO3(2-) + 2H(+); however, the H(+) does not have to be added to the free energy calculations as described by Alberty
         implicit none
         real(8) T, pH, is
         ! Local declarations
         real(8) dG0rxn
         dG0rxn = dGf0(co3sp,T,is,pH) - dGf0(co2totsp,T,is,pH)
         freeCO3 = exp( -dG0rxn/(T*RkJ) ) ! This is unitless ([HCO3]/[DIC])
         return
      end function freeCO3

      real(8) function co2g2l(T, is, pH)
         ! This function returns the ratio of free CO2 (uM) to pCO2 (atm) at temperature T (K)
         ! pH, and ionic strength is (M).  units are (uM/atm).  That is, given the pCO2 (atm) , this returns the
         ! free CO2 in water ([CO2] = co2gas2l*pCO2).
         ! Rxn is:  CO2(aq) -> CO2(g)
         implicit none
         real(8) T, pH, is
         ! Local declarations
         real(8) dGrxn
         dGrxn = dGf0(co2gsp,T,is,pH) - dGf0(co2aqsp,T,is,pH)
         co2g2l = exp(1.d3*dGrxn/(T*RJ))*1.0d6*1.01325d0 ! Last number converts Bar to ATM
         return
      end function co2g2l

      real(8) function solO2(T, is, pH)
         ! This function returns the solubiility of oxygen in water (uM/atm) at temperature T (K),
         ! pH and ionic strength is (M).  It is based on Alberty's values for O2 in gas and aquious phase.
         ! However, these reactions are not affected by pH nor is, so could still have used NIST data,
         ! but this function provides consistency with Alberty's data and approach, so will keep.
         implicit none
         real(8) T, pH, is
         ! Local declarations
         real(8) dGrxn
         ! Use the following reaction:  O2(g) -> O2(aq)
         dGrxn = dGf0(o2aqsp,T,is,pH) - dGf0(o2gsp,T,is,pH)

         solO2 = exp(-1.d3*dGrxn/(T*RJ))*1.0d6*1.01325d0 ! Last number converts Bar to ATM

         return
      end function solO2

      real(8) function solH2S(T, is, pH)
         ! This function returns the solubiility of hydrogen sulfide in water (uM/atm) at temperature T (K),
         ! pH and ionic strength is (M).  It is based on Alberty's values for H2S aquious phase and H2S in the gas phase mentioned above.
         implicit none
         real(8) T, pH, is
         ! Local declarations
         real(8) dGrxn
         ! Use the following reaction:  H2S(g) -> H2S(aq)
         dGrxn = dGf0(h2saqsp,T,is,pH) - dGf0(h2sgsp,T,is,pH)

         solH2S = exp(-1.d3*dGrxn/(T*RJ))*1.0d6*1.01325d0 ! Last number converts Bar to ATM

         return
      end function solH2S

      !    real(8) function solCH4_old(T)
      !      ! This function returns the solubiility of methane in water (uM/atm) at temperature T (K).
      !      ! Data from http://webbook.nist.gov/chemistry/
      !        implicit none
      !        real(8) T
      !
      !        solCH4_old = 1.4d-3*dexp( 1700.*(1.0/T-1.0/298.0) )*1.0d6
      !
      !        return
      !    end function solCH4_old
      !    real(8) function solO2(T)
      !      ! This function returns the solubiility of oxygen in water (uM/atm) at temperature T (K).
      !      ! Data from http://webbook.nist.gov/chemistry/
      !        implicit none
      !        real(8) T
      !
      !        solO2 = 1.3d-3*dexp( 1700.*(1.0/T-1.0/298.0) )*1.0d6
      !
      !        return
      !    end function solO2
      !
      !    real(8) function solCO2(T)
      !        ! This function returns the solubiility of carbon dioxide in water (uM/atm) at temperature T (K).
      !        ! Data from http://webbook.nist.gov/chemistry/
      !        implicit none
      !        real(8) T
      !
      !        solCO2 = 3.4d-2*dexp( 2400.*(1.0/T-1.0/298.0) )*1.0d6
      !
      !        return
      !    end function solCO2

   end module thermoData

max_stars_repo_name: maxEntropyProd/MEP-AutoHetDet

max_issues_repo_name: maxEntropyProd/MEP-AutoHetDet

max_forks_repo_name: maxEntropyProd/MEP-AutoHetDet


text:
C=======================================================================
C  IPSOIL, Subroutine
C
C  Read soil parameters from RESCH???.SDA
C-----------------------------------------------------------------------
C  Revision history
C  06/15/1994 PWW Written
C  02/07/1993 PWW Header revision and minor changes
C  02/07/1993 PWW Added switch common block, restructured
C  06/09/1999 CHP Modular format
C  03/16/2000 GH  Incorporated in CROPGRO
C                 Note, file name should be dynamically created based on
C                 model name
C  03/26/2003 GH  Modified file name and location (SOL directory)
C  08/12/2003 CHP Added I/O error checking
C  08/22/2003 CHP Changed to read file once only and save arrays of data.
C  03/07/2006 CHP Added AM, RCP. Made output arguments optional
!  04/13/2006 CHP Revised format for input file.  Most arguments are
!                 no longer needed, but defaults will be provided
!                 rather than values read from file.  Note that the
!                 values which were previously included in the file
!                 were the same for each residue, so this will not
!                 change results from previous runs.
!  06/28/2007 CHP Added AM, WATFAC to SOILN045.SOL file
!  10/25/2007 CHP/GH Changed name of residue characteristics file to
!                     RESCH045.SDA
!  04/30/2008 CHP Changed units for SCN, SCP, RCN, RCP to %
!  04/30/2008 CHP Path for SDA files set in DSSATPRO file
!  06/25/2008 CHP Changed units for AM from ha/kg to cm2/g
!                               WATFAC from mm-ha/kg to kg[H2O]/kg[DM]
!                 Added EXTFAC - Light extinction coef. for mulch
!                 Check for file version number, use defaults if wrong
C-----------------------------------------------------------------------
C  Called : CROPGRO
C=======================================================================

      SUBROUTINE IPSOIL (CONTROL, RESTYPE, CROP,          !Input
     &    AM, DMINR, DSNC, EXTFAC, PRCEL, PRCHO,          !Output
     &    PRLIG, PSLIG, RCN, RCP, RDCEL, RDCHO, RDLIG,    !Output
     &    SCN, SCP, WATFAC)                               !Output

C-----------------------------------------------------------------------
      USE ModuleDefs     !Definitions of constructed variable types,
                         ! which contain control information, soil
                         ! parameters, hourly weather data.

      IMPLICIT     NONE
      SAVE

      TYPE (ControlType), INTENT(IN) :: CONTROL
      CHARACTER(len=5), INTENT(IN), OPTIONAL :: RESTYPE
      CHARACTER(len=2), INTENT(IN), OPTIONAL :: CROP

      REAL, INTENT(OUT), OPTIONAL :: PRLIG, PSLIG, RCN, RCP, SCN, SCP
      REAL, INTENT(OUT), OPTIONAL :: AM, DMINR, DSNC, EXTFAC, PRCEL
      REAL, INTENT(OUT), OPTIONAL :: PRCHO, RDCEL, RDCHO, RDLIG, WATFAC

      CHARACTER(len=1), PARAMETER :: BLANK = ' '
      CHARACTER(len=6), PARAMETER :: ERRKEY = 'IPSOIL'
      CHARACTER(len=3)   ModelVersionDecimal
      CHARACTER(len=9)   VERSION_TXT
      CHARACTER(len=12)  FILESS, NAMEF
      CHARACTER(len=30)  FILEIO
      CHARACTER(len=78)  MSG(4)
      CHARACTER(len=80)  PATHSD
      CHARACTER(len=92)  SOILNF
      CHARACTER(len=180) CHAR

      INTEGER FOUND, ERR, ISECT, LUNSOL, LNUM
      INTEGER I, J, JJ, LUNIO, NRES, PFLAG
      INTEGER, PARAMETER :: NR=50    !Maximum # of residue types
      REAL REDUCE_FRAC

!     Keep default values for legacy variables
      REAL ADMINR, ADSNC, APRCEL, APRCHO
      REAL ARDCEL, ARDCHO, ARDLIG

!     Saved values from file
      CHARACTER*2, DIMENSION(NR) :: ACROP
      CHARACTER*5, DIMENSION(NR) :: ARESTYPE
      REAL, DIMENSION(NR) :: AAM, AEXTFAC, ARCN, ARCP, ASCN, ASCP
      REAL, DIMENSION(NR) :: APRLIG, APSLIG, AWATFAC

      LOGICAL FIRST, FEXIST, EROR

      DATA FIRST /.TRUE./

C-----------------------------------------------------------------------
!     On first call to routine, read file and save values in arrays for
!       subsequent retrieval.
      IF (FIRST) THEN
        FIRST = .FALSE.
        EROR = .FALSE.

!       Assign default values
        ARESTYPE='RE001'
        AAM    = 32.0     !Area covered / unit dry mass of res (cm2/g)
        AWATFAC= 3.8      !Sat. water content of res. (kg[H2O]/kg[DM])
        AEXTFAC= 0.80     !Light extinction coef. for mulch
        APSLIG = 0.10     !Proportion of lignin is surface res. (%)
        ASCN   = 1.0      !N content of init surface residue (%)
        ASCP   = 0.1      !P content of init surface residue (%)
        APRLIG = 0.10     !Proportion of lignin in root residue (%)
        ARCN   = 1.0      !N content of init root residue (%)
        ARCP   = 0.02     !P content of init root residue (%)

        ADMINR = 8.3E-05
        ADSNC  = 20.0
        ARDCHO = 0.2
        ARDCEL = 0.05
        ARDLIG = 0.0095
        APRCHO = 0.20  !PRCHO, PRCEL, and PRLIG are proportions of
        APRCEL = 0.70  !carbohydrate, cellulose and lignin in residue
        APRLIG = 0.10  !and must sum to 1.0.

        FILEIO  = CONTROL % FILEIO
        LUNIO   = CONTROL % LUNIO

        FILESS = 'RESCH' // ModelVerTxt // '.SDA'
        WRITE(ModelVersionDecimal,'(I1,".",I1)')
     &     Version%Major, Version%Minor

        SOILNF = FILESS
        INQUIRE (FILE = SOILNF, EXIST = FEXIST)
        IF (.NOT. FEXIST) THEN
          CALL PATH('STD',CONTROL%DSSATP,PATHSD,PFLAG,NAMEF)
          SOILNF = TRIM(PATHSD) // FILESS
        ENDIF

        INQUIRE (FILE=SOILNF, EXIST=FEXIST)
        IF (FEXIST) THEN
C-----------------------------------------------------------------------
C       Open the RESCH???.SDA file
C-----------------------------------------------------------------------
          CALL GETLUN('FILESS', LUNSOL)
          OPEN (LUNSOL,FILE = SOILNF, STATUS = 'OLD', IOSTAT=ERR)
          IF (ERR .EQ. 0) THEN

            CALL FIND2(LUNSOL,'@VERSION',LNUM,FOUND)
            IF (FOUND == 0) THEN
!             Version not found in file
              EROR = .TRUE.
              ERR = 50
              GOTO 500
            ELSE
!             Read version number
              CALL IGNORE(LUNSOL,LNUM,ISECT,CHAR)
              READ(CHAR,'(A9)',IOSTAT=ERR) VERSION_TXT
              IF (INDEX(VERSION_TXT,ModelVersionDecimal) < 1 .OR.
     &                ERR /= 0) THEN
!               Version not found in file, or wrong version
                EROR = .TRUE.
                ERR = 50
                GOTO 500
              ENDIF
            ENDIF
C-----------------------------------------------------------------------
C     Read Soil Parameters from RESCH???.SDA
C-----------------------------------------------------------------------
            I = 0; ISECT = 1
!           EOF not portable. CHP 7/24/2007
!           ReadLoop: DO WHILE (.NOT. EOF (LUNSOL))
            ReadLoop: DO WHILE (ISECT > 0)
              CALL IGNORE(LUNSOL,LNUM,ISECT,CHAR)
              IF (ISECT == 0) EXIT ReadLoop
              IF (CHAR(1:1) .EQ. '*' .OR. CHAR(1:1) .EQ. '@') THEN
                CYCLE ReadLoop
              ENDIF
              I = I + 1
              READ(CHAR,'(1X,A5,1X,A2,3F8.0,2(F8.0,2F6.0))',IOSTAT=ERR)
     &            ARESTYPE(I), ACROP(I), AAM(I), AWATFAC(I), AEXTFAC(I),
     &            APSLIG(I), ASCN(I), ASCP(I),
     &            APRLIG(I), ARCN(I), ARCP(I)
              IF (ERR .NE. 0) THEN
                I = I - 1
                EROR = .TRUE.
              ENDIF
            ENDDO ReadLoop
          ELSE
            !File could not be opened
            EROR = .TRUE.
            ERR = 20
          ENDIF
        ELSE
          !File does not exist
          EROR = .TRUE.
          ERR = 40
        ENDIF

        NRES = I

        CONTINUE
        CLOSE (LUNSOL)

C-----------------------------------------------------------------------
  500   CONTINUE
!       Use default values if problems with reading file
        IF (EROR) THEN
          SELECT CASE(ERR)
            CASE (20)
              MSG(1) = 'Error opening file: '
            CASE (40)
              MSG(1) = 'File does not exist: '
            CASE (50)
              MSG(1) = 'Incorrect version for file:'
            CASE DEFAULT
              MSG(1) = 'Error reading file: '
          END SELECT
          MSG(2) = SOILNF(1:78)
          MSG(3) = 'Default residue characteristics will be used.'
          CALL WARNING(3,ERRKEY,MSG)
        ENDIF
      ENDIF       !End of FIRST block

C-----------------------------------------------------------------------
!     Assign values independant of restype (legacy values)
      IF (PRESENT(DMINR)) DMINR = ADMINR
      IF (PRESENT(DSNC )) DSNC  = ADSNC
      IF (PRESENT(RDCHO)) RDCHO = ARDCHO
      IF (PRESENT(RDCEL)) RDCEL = ARDCEL
      IF (PRESENT(RDLIG)) RDLIG = ARDLIG
      IF (PRESENT(PRCHO)) PRCHO = APRCHO
      IF (PRESENT(PRCEL)) PRCEL = APRCEL

!     Get appropriate residue or crop type, if included
      IF (PRESENT(RESTYPE) .OR. PRESENT(CROP)) THEN
        RLoop: DO J = 1, NR
          IF (PRESENT(RESTYPE)) THEN
            IF (RESTYPE == '     ') THEN
              JJ = 1
              EXIT RLoop
            ENDIF
            IF (RESTYPE .EQ. ARESTYPE(J)) THEN
              JJ = J
              EXIT RLoop
            ENDIF
          ELSEIF (PRESENT(CROP)) THEN
            IF (CROP == '  ') THEN
              JJ = 1
              EXIT RLoop
            ENDIF
            IF (CROP .EQ. ACROP(J)) THEN
              JJ = J
              EXIT RLoop
            ENDIF
          ENDIF
        ENDDO RLoop

        IF (J > NR) THEN
          JJ = 1
          IF (PRESENT(RESTYPE)) THEN
            WRITE(MSG(1),"('Residue type: ',A5,' not valid.')") RESTYPE
            MSG(2) = 'Default residue characteristics will be used.'
            CALL WARNING(2,ERRKEY,MSG)
          ENDIF
        ENDIF

      ELSE !Default to residue type 1
        JJ = 1
      ENDIF

      IF (PRESENT(AM    )) THEN
        AM    = AAM(JJ)
!       Check for incorrect units
        IF (AM < 0.01) THEN
          AM = AM * 1.E5
          MSG(1) = "Units of mulch aerial coverage appear to be wrong."
          WRITE(MSG(2),'(A,F8.2,A)') "Value changed to ",AM," cm2/g"
          MSG(3) = "May need to update file: "
          MSG(4) = SOILNF(1:78)
          CALL WARNING(4,ERRKEY,MSG)
        ENDIF
      ENDIF

      IF (PRESENT(WATFAC)) THEN
        WATFAC= AWATFAC(JJ)
!       Check for incorrect units
        IF (AM < 0.001) THEN
          WATFAC = WATFAC * 1.E4
          MSG(1) ="Units of mulch saturation factor appear to be wrong."
          WRITE(MSG(2),'(A,F8.3,A)')
     &      "Value changed to ",WATFAC," kg[H2O]/KG[DM]"
          MSG(3) = "May need to update file: "
          MSG(4) = SOILNF(1:78)
          CALL WARNING(4,ERRKEY,MSG)
        ENDIF
      ENDIF

      IF (PRESENT(EXTFAC)) EXTFAC= AEXTFAC(JJ)
      IF (PRESENT(PSLIG )) PSLIG = APSLIG(JJ)
      IF (PRESENT(SCN   )) SCN   = ASCN  (JJ)
      IF (PRESENT(SCP   )) SCP   = ASCP  (JJ)

!     PRLIG must sum with PRCHO and PRCEL to 1.0
      IF (PRESENT(PRLIG )) THEN
        PRLIG = APRLIG(JJ)
        IF (PRESENT(PRCHO) .AND. PRESENT(PRCEL)) THEN
          REDUCE_FRAC = (1.0 - PRLIG) / (PRCHO + PRCEL)
          PRCHO = PRCHO * REDUCE_FRAC
          PRCEL = PRCEL * REDUCE_FRAC
        ENDIF
      ENDIF

      IF (PRESENT(RCN   )) RCN   = ARCN  (JJ)
      IF (PRESENT(RCP   )) RCP   = ARCP  (JJ)

C-----------------------------------------------------------------------
      RETURN
      END SUBROUTINE IPSOIL

C=======================================================================

! IPSOIL and SCREATE Variables
!
! ACRO(I) List of variables names used in RESCH???.SDA file which contains
!           coefficients needed for simulating the decomposition of soil
!           organic matter and organic matter
! CHAR    Contains the contents of last record read
! DMINR   Maximum decomposition rate constant of stable organic matter
!           (d-1)
! DSNC    Depth to which C and N are integrated across all soil layers for
!           output in CARBON.OUT (cm)
! ERR     Error code for file operation
! FILESS  Path plus filename for species file (*.spe)
! I       Loop counter
! J       Loop counter
! LUNSOL  Logical unit number for RESCH???.SDA
! PRCEL   Cellulose fraction of the residue (fraction)
! PRCHO   Carbohydrate fraction of the residue (fraction)
! PRLIG   Lignin fraction of the residue (fraction)
! RCN     N content of initial root residue (%)
! RCP     P content of initial root residue (%)
! SCN     N content of initial shoot residue (%)
! SCP     P content of initial shoot residue (%)
! RDCEL   Maximum decomposition rate of cellulose (fraction / day)
! RDCHO   Maximum decomposition rate of carbohydrates (fraction / day)
! RDLIG   Maximum decomposition rate of lignin (fraction / day)
! RESTYPE Residue application method as sent from calling routine
! T(I)    Text array with explanatory text of the RESCH???.SDA file
!=======================================================================

!=======================================================================
      MODULE Interface_IpSoil
!     Interface needed for optional arguments with IPSOIL
      INTERFACE
        SUBROUTINE IPSOIL (CONTROL, RESTYPE, CROP,        !Input
     &    AM, DMINR, DSNC, EXTFAC, PRCEL, PRCHO,          !Output
     &    PRLIG, PSLIG, RCN, RCP, RDCEL, RDCHO, RDLIG,    !Output
     &    SCN, SCP, WATFAC)                               !Output
          USE ModuleDefs
          TYPE (ControlType), INTENT(IN) :: CONTROL
          CHARACTER(len=5), INTENT(IN), OPTIONAL :: RESTYPE
          CHARACTER(len=2), INTENT(IN), OPTIONAL :: CROP
          REAL, INTENT(OUT), OPTIONAL :: AM, DMINR, DSNC, EXTFAC
          REAL, INTENT(OUT), OPTIONAL :: PRCEL, PRCHO, PRLIG, PSLIG
          REAL, INTENT(OUT), OPTIONAL :: RCN, RCP, RDCEL, RDCHO, RDLIG
          REAL, INTENT(OUT), OPTIONAL :: SCN, SCP, WATFAC
        END SUBROUTINE
      END INTERFACE
      END MODULE Interface_IpSoil

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

max_stars_repo_name: vstocca/dssat-csm-os

max_issues_repo_name: vstocca/dssat-csm-os

max_forks_repo_name: vstocca/dssat-csm-os


text:
c test pure directive
      program pure_1
      external cube
      real cube
cfcd$ pure cube
      real A(100)
chpf$ processors P(4)
chpf$ template T(100)
chpf$ align A with T
chpf$ distribute T(block) onto P
      integer i
chpf$ independent
      do i=1, 100
        A(i) = 1.0/real(i)
      enddo
chpf$ independent
      do i=1, 100
        A(i) = cube(A(i))
      enddo
      print *, A(51)
      end
      real function cube(x)
cfcd$ pure
      real x
      cube=x*x*x
      end


text:
cdis
cdis    Open Source License/Disclaimer, Forecast Systems Laboratory
cdis    NOAA/OAR/FSL, 325 Broadway Boulder, CO 80305
cdis
cdis    This software is distributed under the Open Source Definition,
cdis    which may be found at http://www.opensource.org/osd.html.
cdis
cdis    In particular, redistribution and use in source and binary forms,
cdis    with or without modification, are permitted provided that the
cdis    following conditions are met:
cdis
cdis    - Redistributions of source code must retain this notice, this
cdis    list of conditions and the following disclaimer.
cdis
cdis    - Redistributions in binary form must provide access to this
cdis    notice, this list of conditions and the following disclaimer, and
cdis    the underlying source code.
cdis
cdis    - All modifications to this software must be clearly documented,
cdis    and are solely the responsibility of the agent making the
cdis    modifications.
cdis
cdis    - If significant modifications or enhancements are made to this
cdis    software, the FSL Software Policy Manager
cdis    (softwaremgr@fsl.noaa.gov) should be notified.
cdis
cdis    THIS SOFTWARE AND ITS DOCUMENTATION ARE IN THE PUBLIC DOMAIN
cdis    AND ARE FURNISHED "AS IS."  THE AUTHORS, THE UNITED STATES
cdis    GOVERNMENT, ITS INSTRUMENTALITIES, OFFICERS, EMPLOYEES, AND
cdis    AGENTS MAKE NO WARRANTY, EXPRESS OR IMPLIED, AS TO THE USEFULNESS
cdis    OF THE SOFTWARE AND DOCUMENTATION FOR ANY PURPOSE.  THEY ASSUME
cdis    NO RESPONSIBILITY (1) FOR THE USE OF THE SOFTWARE AND
cdis    DOCUMENTATION; OR (2) TO PROVIDE TECHNICAL SUPPORT TO USERS.
cdis
cdis
cdis
cdis
cdis

!       1997 Jun        Ken Dritz     Added calls to get_grid_dim_xy and
!                                     get_laps_dimensions to get values of
!                                     NX_L, NY_L, NZ_L.
!       1997 Jun        Ken Dritz     Now pass NX_L, NY_L, NZ_L to laps_accum.

c       integer j_status(20),iprod_number(20),i4time_array(20)
        integer j_status(20),iprod_number(20)
        character*9 a9_time

        call get_systime(i4time,a9_time,istatus)
        if(istatus .ne. 1)go to 999

        call get_grid_dim_xy(NX_L,NY_L,istatus)
        if (istatus .ne. 1) then
           write (6,*) 'Error getting horizontal domain dimensions'
           go to 999
        endif

        call get_laps_dimensions(NZ_L,istatus)
        if (istatus .ne. 1) then
           write (6,*) 'Error getting vertical domain dimension'
           go to 999
        endif

        call get_max_radar_files(MAX_RADAR_FILES, istatus)
        if (istatus .ne. 1) then
           write (6,*) 'Error getting max_radar_files'
           go to 999
        endif

        call laps_accum   (i4time,
     1                     NX_L,NY_L,NZ_L,MAX_RADAR_FILES,
     1                     i_diag,
     1                     n_prods,
     1                     iprod_number,
     1                     j_status)

999     continue

        end

max_stars_repo_name: maxinye/laps-mirror

max_stars_repo_licenses: Intel Unlicense OLDAP-2.2.1 NetCDF

max_issues_repo_name: longwosion/laps-mirror

max_issues_repo_licenses: Intel NetCDF OLDAP-2.2.1 Unlicense

max_forks_repo_name: longwosion/laps-mirror

max_forks_repo_licenses: Intel NetCDF OLDAP-2.2.1 Unlicense


text:
    !-------------------------------------------------------------------------------------------------------------
    !
    !> \file    ParseCSDataBlockSUBI.f90
    !> \brief   Parse the data block section corresponding to a SUBL phase of a ChemSage data-file.
    !> \author  M. Poschmann
    !> \date    Jan. 19, 2021
    !> \todo    Add capability to parse magnetic mixing terms in a SUBLM phase.
    !> \sa      ParseCSDataFile.f90
    !> \sa      ParseCSDataBlock.f90
    !> \sa      ParseCSDataBlockGibbs.f90
    !
    !
    ! DISCLAIMER
    ! ==========
    !
    ! All of the programming herein is original unless otherwise specified and is completely
    ! independent of ChemApp and related products, including Solgas, Solgasmix, Fact, FactSage
    ! and ChemSage.
    !
    !
    ! Revisions:
    ! ==========
    !
    !   Date            Programmer      Description of change
    !   ----            ----------      ---------------------
    !   01/19/2021      M. POschmann     Original code
    !
    !
    ! Purpose:
    ! ========
    !
    !> \details The purpose of this subroutine is to parse the "data block" section of a ChemSage data-file
    !! containing a "SUBI" phase (Ionic Liquid Model).
    !!
    !!
    !! The coefficients for the standard molar Gibbs energy equations for pure species originate from a ChemSage
    !! data-file that was parsed from the ParseCSDataFile program.  The format for the coefficients follow:
    !!
    !! \f$ g_i^{\circ} = A + BT + CTln(T) + DT^2 + ET^3 + F/T + (GT^U + HT^V + IT^W + JT^X + KT^Y + LT^Z) \f$
    !!
    !
    !
    !
    ! Pertinent variables:
    ! ====================
    !
    !> \param[in] i                 An integer scalar representing the solution phase index.
    !> \param[in] nCountSublatticeCS  An integer scalar representing a cummulative count of the number of phases
    !                            containing a sublattice.
    !
    ! INFO                      A scalar integer that indicates a successful exit or identifies an error.
    ! cSolnPhaseNameCS          A character vector representing the name of each solution phase.
    ! cSolnPhaseTypeCS          A character vector representing the type of each solution phase.
    ! nSublatticePhaseCS        An integer vector representing the number of sublattices per phase.
    ! dStoichSublatticeCS       A double real matrix representing the stoichiometry coefficient for each
    !                            sublattice of a phase.
    ! nConstituentSublatticeCS  An integer matrix representing the number of constituents per sublattice of a
    !                            particular phase.
    ! cConstituentNameSUBCS   A character matrix representing the names of each constituent on each sublattice
    !                            for a particular phase.
    ! iConstituentSublatticeCS  An integer matrix representing the constituent indices on each sublattice.
    ! iRegularParamCS           An integer matrix representing the index of the mixing parameter and important
    !                            information for mixing.  For iRegularParamCS(i,j) for a SUBL phase, i refers to
    !                            the mixing parameter index and j refers to a vector specific to that mixing term.
    !                            The first coefficient along j refers to the number of constituents (n) involved for
    !                            that parameter, coefficients 2:n+1 represent the indices of constituents,
    !                            coefficient n+2 represents the order of the mixing parameter.
    !
    !-------------------------------------------------------------------------------------------------------------

subroutine ParseCSDataBlockSUBI(i)

    USE ModuleParseCS

    implicit none

    integer                   :: c, i, j, k, n, s, v, iDummy, iSubLat
    logical                   :: lTripleTerm
    character(8)              :: cDummy
    ! Counter to determine if there are two or more SUBI
    ! phases present for a miscibility gap.
    iMiscSUBI = iMiscSUBI + 1

    if (iMiscSUBI == 1) then
        if (allocated(iSUBIMixTypeCS)) deallocate(iSUBIMixTypeCS)
        if (allocated(iSUBIParamDataCS)) deallocate(iSUBIParamDataCS)
        allocate(iSUBIMixTypeCS(1000),iSUBIParamDataCS(1000,6))
        iSUBIMixTypeCS = 0
        iSUBIParamDataCS = 0
    end if

    ! Initialize variables:
    n = 0

    ! Two sublattices per phase:
    nSublatticePhaseCS(nCountSublatticeCS) = 2

    ! Report an error if the number of sublattices is out of range:
    if ((nSublatticePhaseCS(nCountSublatticeCS) > nMaxSublatticeCS).OR. &
        (nSublatticePhaseCS(nCountSublatticeCS) <= 0)) then
        INFO = 33
        return
    end if

    ! Read in the number of constituents for each sublattice:
    read (1,*,IOSTAT = INFO) nConstituentSublatticeCS(nCountSublatticeCS,1:nSublatticePhaseCS(nCountSublatticeCS))

    ! Read in the name of each constituent for each sublattice:
    do s = 1, nSublatticePhaseCS(nCountSublatticeCS)
        read (1,*,IOSTAT = INFO) cConstituentNameSUBCS(nCountSublatticeCS,s,1:nConstituentSublatticeCS(nCountSublatticeCS,s))
    end do

    ! Convert Va, va, Va, and vA cases to Va
    do v = 1, nConstituentSublatticeCS(nCountSublatticeCS,2)
      cDummy = cConstituentNameSUBCS(nCountSublatticeCS,2,v)

      if ((cDummy == 'VA') .OR. &
          (cDummy == 'va') .OR. &
          (cDummy == 'Va') .OR. &
          (cDummy == 'vA')) then

          cConstituentNameSUBCS(nCountSublatticeCS,2,v) = 'Va'

      end if
    end do

    ! Record an error if necessary:
    if (INFO /= 0) then
        INFO = 2100 + i
        return
    end if

    ! I think this is # of sublattices (always 2), but need that earlier and no reason to overwrite
    read (1,*,IOSTAT = INFO) iDummy

    ! Read charges of constituents
    do s = 1, nSublatticePhaseCS(nCountSublatticeCS)
        read (1,*,IOSTAT = INFO) dSublatticeChargeCS(nCountSublatticeCS,s,1:nConstituentSublatticeCS(nCountSublatticeCS,s))
    end do

    ! Read in the constituent indices for each component on each sublattice:
    do s = 1, nSublatticePhaseCS(nCountSublatticeCS)
        ! Number of components for this phase:
        n = nSpeciesPhaseCS(i) - nSpeciesPhaseCS(i-1)

        read (1,*,IOSTAT = INFO) iConstituentSublatticeCS(nCountSublatticeCS, s, 1:n)
    end do

    ! Record an error and return if necessary:
    if (INFO /= 0) then
        INFO = 2200 + i
        return
    end if

    ! Loop through excess mixing parameters:
    LOOP_ExcessMixingSUBI: do
        ! Read in number of constituents involved in parameter:
        read (1,*,IOSTAT = INFO) iRegularParamCS(nParamCS+1,1)

        ! The end of the section of mixing terms is labelled "0".
        if (iRegularParamCS(nParamCS+1,1) == 0) exit LOOP_ExcessMixingSUBI

        ! Update counter of the number of parameters:
        nParamCS = nParamCS + 1

        ! Read in the list of constituent indices involved in this parameter:
        j = iRegularParamCS(nParamCS,1)

        read (1,*,IOSTAT = INFO) iRegularParamCS(nParamCS,2:j+2)

        ! Read in the mixing parameter:
        read (1,*,IOSTAT = INFO) dRegularParamCS(nParamCS,1:6)

        ! Store number of mixing parameters per component array:
        n = iRegularParamCS(nParamCS, j+2)

        ! Correct the indexing scheme (order of mixing) for the first parameter:
        iRegularParamCS(nParamCS, j+2) = 0

        LOOP_Mixing_Cases: do k = 2,iRegularParamCS(nParamCS,1) + 1

            ! Number of parameters in mixing case
            iSUBIParamDataCS(nParamCS,1) = iRegularParamCS(nParamCS,1)

            iSubLat = iRegularParamCS(nParamCS,k) - nConstituentSublatticeCS(nCountSublatticeCS,1)

            if (iSubLat <= 0) then
                s = 1
                c = iRegularParamCS(nParamCS,k)
            else
                s = 2
                c = iSubLat
            end if

            ! Distinguishing between SUBI mixing Case types 1 through 12
            ! Cation
            if (dSublatticeChargeCS(nCountSublatticeCS,s,c) > 0) then
                iSUBIParamDataCS(nParamCS,2) = iSUBIParamDataCS(nParamCS,2) + 1
            ! Nuetral
            else if (dSublatticeChargeCS(nCountSublatticeCS,s,c) == 0) then
                iSUBIParamDataCS(nParamCS,3) = iSUBIParamDataCS(nParamCS,3) + 1
            ! Vacancy
            else if (cConstituentNameSUBCS(nCountSublatticeCS,s,c) == 'Va') then
                iSUBIParamDataCS(nParamCS,4) = iSUBIParamDataCS(nParamCS,4) + 1
            ! Anion
            else
                iSUBIParamDataCS(nParamCS,5) = iSUBIParamDataCS(nParamCS,5) + 1
                ! Checking for cases with Dl
                if ((iSUBIParamDataCS(nParamCS,5) > 0) .AND. &
                    ((iSUBIParamDataCS(nParamCS,1) - iSUBIParamDataCS(nParamCS,2)) > 1)) then
                    ! if Dl present
                    iSUBIParamDataCS(nParamCS,6) = 1
                end if
            end if

        end do LOOP_Mixing_Cases

        ! Case 1: Determine the mixing parameter type: L_Ci,Cj:Ak
        if ((iSUBIParamDataCS(nParamCS,1) == 3) .AND. &
            (iSUBIParamDataCS(nParamCS,2) == 2) .AND. &
            (iSUBIParamDataCS(nParamCS,5) == 1)) then

            iSUBIMixTypeCS(nParamCS) = 1

        ! Case 2: Determine the mixing parameter type: L_Ci:Aj,Dk
        else if ((iSUBIParamDataCS(nParamCS,1) == 3) .AND. &
                 (iSUBIParamDataCS(nParamCS,2) == 1) .AND. &
                 (iSUBIParamDataCS(nParamCS,6) == 1)) then

            iSUBIMixTypeCS(nParamCS) = 2

        ! Case 3: Determine the mixing parameter type: L_Ci,Cj:Va
        else if ((iSUBIParamDataCS(nParamCS,1) == 3) .AND. &
                 (iSUBIParamDataCS(nParamCS,2) == 2) .AND. &
                 (iSUBIParamDataCS(nParamCS,4) == 1)) then

            iSUBIMixTypeCS(nParamCS) = 3

        ! Case 4: Determine the mixing parameter type: L_Ci:Va,Bj
        else if ((iSUBIParamDataCS(nParamCS,1) == 3) .AND. &
                 (iSUBIParamDataCS(nParamCS,2) == 1) .AND. &
                 (iSUBIParamDataCS(nParamCS,3) == 1) .AND. &
                 (iSUBIParamDataCS(nParamCS,4) == 1)) then

            iSUBIMixTypeCS(nParamCS) = 4

        ! Case 5: Determine the mixing parameter type: L_Ci:Bj,Bk
        else if ((iSUBIParamDataCS(nParamCS,1) == 3) .AND. &
                 (iSUBIParamDataCS(nParamCS,2) == 1) .AND. &
                 (iSUBIParamDataCS(nParamCS,3) == 2)) then

            iSUBIMixTypeCS(nParamCS) = 5

        ! Case 6: Determine the mixing parameter type: L_Ci,Cj,Ck:Al
        else if ((iSUBIParamDataCS(nParamCS,1) == 4) .AND. &
                 (iSUBIParamDataCS(nParamCS,2) == 3) .AND. &
                 (iSUBIParamDataCS(nParamCS,5) == 1)) then

            iSUBIMixTypeCS(nParamCS) = 6

        ! Case 7: Determine the mixing parameter type: L_Ci:Aj,Dk,Dl
        else if ((iSUBIParamDataCS(nParamCS,1) == 4) .AND. &
                 (iSUBIParamDataCS(nParamCS,2) == 1) .AND. &
                 (iSUBIParamDataCS(nParamCS,6) == 1)) then

            iSUBIMixTypeCS(nParamCS) = 7

        ! Case 8: Determine the mixing parameter type: L_Ci,Cj,Ck:Va
        else if ((iSUBIParamDataCS(nParamCS,1) == 4) .AND. &
                 (iSUBIParamDataCS(nParamCS,2) == 3) .AND. &
                 (iSUBIParamDataCS(nParamCS,4) == 1)) then

            iSUBIMixTypeCS(nParamCS) = 8

        ! Case 9: Determine the mixing parameter type: L_Ci:Va,Bj,Bk
        else if ((iSUBIParamDataCS(nParamCS,1) == 4) .AND. &
                 (iSUBIParamDataCS(nParamCS,2) == 1) .AND. &
                 (iSUBIParamDataCS(nParamCS,3) == 2) .AND. &
                 (iSUBIParamDataCS(nParamCS,4) == 1)) then

            iSUBIMixTypeCS(nParamCS) = 9

        ! Case 10: Determine the mixing parameter type: L_Ci:Bj,Bk,Bl
        else if ((iSUBIParamDataCS(nParamCS,1) == 4) .AND. &
                 (iSUBIParamDataCS(nParamCS,2) == 1) .AND. &
                 (iSUBIParamDataCS(nParamCS,3) == 3)) then

            iSUBIMixTypeCS(nParamCS) = 10

        ! Case 11: Determine the mixing parameter type: L_Ci,Cj:Ak,Dl
        else if ((iSUBIParamDataCS(nParamCS,1) == 4) .AND. &
                 (iSUBIParamDataCS(nParamCS,2) == 2) .AND. &
                 (iSUBIParamDataCS(nParamCS,6) == 1)) then

            iSUBIMixTypeCS(nParamCS) = 11

        ! Case 12: Determine the mixing parameter type: L_Ci,Cj:Va,Bk
        else if ((iSUBIParamDataCS(nParamCS,1) == 4) .AND. &
                 (iSUBIParamDataCS(nParamCS,2) == 2) .AND. &
                 (iSUBIParamDataCS(nParamCS,3) == 1) .AND. &
                 (iSUBIParamDataCS(nParamCS,4) == 1)) then

            iSUBIMixTypeCS(nParamCS) = 12
        end if

        ! lTripleTerm is false unless the specified conditions are met
        lTripleTerm = .false.
        !Some mixing cases require three interaction parameters to be created
        ! if the term is L_0
        !If one mixing term create three mixing terms
        if (n == 1) then
            if ((iSUBIMixTypeCS(nParamCS) == 6) .OR. &
                (iSUBIMixTypeCS(nParamCS) == 7) .OR. &
                (iSUBIMixTypeCS(nParamCS) == 8) .OR. &
                (iSUBIMixTypeCS(nParamCS) == 9) .OR. &
                (iSUBIMixTypeCS(nParamCS) == 10) .OR. &
                (iSUBIMixTypeCS(nParamCS) == 12))then

                ! Adds 2 mixing terms
                n = 3
                ! If lTripleTerm is true then interaction parameters will be added
                lTripleTerm = .TRUE.
            end if
        end if

        ! Loop through number of mixing terms per component array:
        do k = 2, n
            ! Update counter of the number of parameters:
            nParamCS = nParamCS + 1

            ! Add additional terms:
            iRegularParamCS(nParamCS,1:nParamMax*2+1) = iRegularParamCS(nParamCS-1,1:nParamMax*2+1)

            ! Correct the indexing scheme (order of mixing parameter):
            iRegularParamCS(nParamCS, j+2) = k - 1

            ! Duplicate SUBI mixing type for the same case
            iSUBIMixTypeCS(nParamCS) = iSUBIMixTypeCS(nParamCS - 1)

            ! Re-read in mixing parameters
            if (lTripleTerm) backspace(UNIT = 1)

            ! Read mixing terms:
            read (1,*,IOSTAT = INFO) dRegularParamCS(nParamCS,1:6)

        end do

    end do LOOP_ExcessMixingSUBI

    ! Report an error and return if necessary:
    if (INFO /= 0) then
        INFO = 2400 + i
        return
    end if

    return

end subroutine ParseCSDataBlockSUBI

max_stars_repo_name: ORNL-CEES/thermochimica

max_issues_repo_name: ORNL-CEES/thermochimica

max_forks_repo_name: ORNL-CEES/thermochimica


text:
Address(East 8th Street) is one half of 8th Street which changes names and address numbers at A Street.  At its eastern end East 8th Street curves / turns into Tulip Lane.

The addresses are not continuous which is a bane to new Pizza pizza delivery drivers. So if you live on 8th and are ordering a delivery, be very specific about which part of 8th you live on. Some have said in the past 8th is a good alternative to Russell Boulevard to avoid police; however, as of 2011 speed traps on 8th are rather common.  Not to mention that 8th is a residential street and if you speed on it the residents are likely to get really pissed off at you.

While the Speed limits speed limit is only 25mph for the length of 8th St, the traffic lights are timed quite well, making it an (often quicker) alternate to Russell Boulevard Russell Blvd.

West 8th Street, with Sycamore Lane, was added in 1967 as one of the few new postWWII streets with bike lanes.  Then UCD Professor Bob Sommer brought the idea to city planners in the mid 60s.  The initial design was a success, and became the standard bike lane design for the rest of Davis, California.  At first, bike lanes only went from Sycamore Lane to B Street, then later from J Street to Pole Line Road.  The gap from B Street to L Street was too narrow for bike lanes and parking.  Eventually, daytime parking was banned, then parking was permanently banned, and now 8th Street is easily the best street to go across town.  Had the parkers won a victory in 1970, it would be much more dangerous for bicyclists.  8th Street now has bike lines for its entire length, and specific bicyclistoperated traffic signals at the intersections with L Street and B Street.

Nowadays, the E. 8th St. corridor is a vital corridor for cyclists.  However, the stretch between B and F St. is especially narrow and residents put out toters for waste, recycling and yard clippings.  There should be no yard waste piles on this section.  High volumes of commuters, businesstruck owners/contractors, delivery vehicles and UNITRANS use this section.  Cycling lanes are marked but the wellposted 25 MPH speed limit is largely ignored by those hurrying between lights.

The intersection of E. 8th and B Streets is dicey for several reasons.  It is rather problematic and has had several collisions and there are many nearmisses.  First, five different elementary and secondary schools are within a twoblock radius, plus the Public Library.   That means TONS of cyclists traverse through here each day.  Secondly, as commute traffic lines up in the morning, lanes back up since there is but only one leftturn lane (a free leftturn on westbound E. 8th to southbound B).  Some motorists dodge around leftturners who are blocking their passage and use the bike lanes.  City officials claim that the traffic signal control unit is one of the oldest in town and cannot handle a controlled leftturn phase.  If it helps, this is the corner where the Davis Lutheran Church is located.

Be careful bicycling through the F Street intersection on 8th Street.  There are no left turn lanes for cars, so a bicyclist can be hidden behind a samedirectioncar while approaching the light, and then hit by a left turner in the opposite direction after the samedirection car has passed.  Or, if you are passing a car that is waiting to turn left, you are hidden from view and may have opposing left turners turn in front of you or hit you.  Several serious, nearly fatal accidents have occurred this way in recent years.

Before East 8th St. was built, Cemetery Road extended from G Street east across the railroad tracks to the Davis Cemetery.  Sweet Briar Road and Pennsylvania Place are remnants of Cemetery Road, and the section from J Street to the cemetery is now part of 8th Street.  Part of old Cemetery Road can be seen if you walk into the bushes from the corner of Sweet Briar and H Street  there is some concrete and a yellow center line on an abandoned piece of street where it once crossed the railroad tracks.  8th Street ended just east of the cemetery until the mid 1990s when it was connected to Tulip Lane.

Notable Locations along 8th Street from include:

 Food and Shopping
   Davis Manor Center
   Quick Shop Market

 Landmarks
  Davis Cemetery
  Train Tracks Railroad Tracks
  Davis Fire Crew parking
  Sunrise Farm

 Services
  A Fair Go
  Davis Lutheran Church at the corner of 8th and B Street
  Valley Oak Elementary School
  Pole Line Road Baptist Church

 Businesses
  Sherman Homes

 Neighborhoods
  Davis Manor

 Apartments
  Redwood Tree Apartments
  Regency Arms Townhouses Apartments
  Villa Calabria

 Streets
  A Street
  B Street
  D Street
  F Street
  G Street
  H Street & Sweet Briar Road  Sweet Briar Road may not reach all the way to 8th, it may already be H by then.
  Railroad Tracks
  J Street
  Hunt Way
  K Street
  L Street
  Simmons Way
  M Street
  Chestnut Place
  Pole Line Road
  North Diameter Drive
  Mesquite Drive
  Tulip Lane

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search