text:
c=======================================================================
c Soil nutrient uptake model with switches file processing.
c Modified from OurPrism2.f by Jonathan Rees.

c User documentation is in a separate file.  Don't get caught without it.

c Switch options: (default if not specified is "calculated")
c   Fixed <value>
c   Input
c   Calculated
c   Calculated <parameter>
c   Linear <initial>
c   Exponential <initial>
c   Output

c State variables and parameters:
c  Ruth     Carol
c   name     name    Units
c  ----     ------   -----
c  alpha    RtAbso   cm/s       Root absorbing power
c  b        SolBuf   1          Soil buffer capacity (cm3/cm3)
c  C0       CSurf    mol/cm3    Concentration at root surface
c  Cav      SChem    mol/cm3    Average nutrient concentration
c  De       DiffCo   cm2/s      Effective diffusion coefficient
c  dt       Interval s          Duration of a timestep
c  f        --       1          Fugacity (for computing v0) (NYI)
c  f_slope  --       1/s        Change in f (NYI)
c  growth   DlyGrwth 1/s        Growth rate (exponential; cm/cm/s)
c  growth   --       cm/s       Growth rate (linear; cm/s)
c  Km       Km       mol/cm3    Solution concentration at which net influx
c  L        cm       cm         Root length
c  r0       RtRad    cm         Root radius
c  rx       RadAve   cm         Average radial distance to the next root's
c                                zone of influence
c  source   --       mol/cm3/s  Constant nutrient supply/depletion (Cav pump)
c  t        --       s          Time
c  theta    --       1          Water content (cm3/cm3)   NYI
c  U        UPTK     mol        Uptake (this is the main output of the model)
c  Unew     Unew     mol        Uptake of new roots
c  steps    --       1          Number of timesteps to run
c  v0       VelIn    cm/s       Inward radial velocity of water at root surf
c  Vmax     Imax     mol/cm2/s  Maximum influx at high concentration is 1/2 Vmax
c  w0       RtWUp    cm3/s      Rate of water uptake
c
c Obsolescent (left over from OurPrism2):
c  --       Cmin     mol/cm3    Concentration at which alpha goes to 0
c  iter     Iteration --        Number of steps so far
c  oalpha   RtAbso   cm/s       Root absorbing power (calculated iteratively)

c To promote portability, I try to code to the FORTRAN 90 standard.  I
c have checked conformance using MS Fortran's "Adhere to Fortran 90
c Language Standard" option with array bounds and other checks enabled.
c
c Language features used:
c -- IF (...) THEN ... ELSEIF ... ENDIF
c -- STATUS= in OPEN statements
c -- Blank lines
c -- "!" comments
c -- Underscore as an alphabetic character.
c
c Language features not used:
c -- REAL*8 isn't official (this is an IBMish extension).  Replaced with
c    DOUBLE PRECISION.
c -- Statements not to go past column 72.
c -- I assume that only the first 6 characters of an identifier are
c    significant.
c -- Fortran 90 probably has some structured looping constructs, but I
c    have not bothered to learn them.
c
c Style notes:
c
c -- File extension ".fi" means "Fortran include" and is the extension
c    conventionally used by Microsoft Fortran Powerstation.
c
c -- I start statements in column 7, not in column 8 (as in many of
c    the TREGRO files).  Less whitespace (up to a point) means less
c    wasted space and therefore makes programs more readable, because
c    you can take in more information at a glance.
c
c -- Indentation: 2 characters for progressive indentation of IF's and
c    loops.
c
c -- Case:  Upper case is traditional for FORTRAN, but unnecessary
c    these days.  I am undecided which is more readable.
c    Use upper or lower case 'c' for comments?
c
c -- At present all the source code is in one file, Nut.f, with one
c    include file, Nut.fi.  This is better for development since it
c    makes global searches a little easier.  Also it's silly to repeat
c    the same long comments at the top of each file.  There is a
c    glossary of variables at the top, but it is not as detailed as
c    Carol's.

c *------------------------------------------------------------
c * Static information: Names of state variables.
c *------------------------------------------------------------

      BLOCK DATA MUMBLE
      INCLUDE 'Nut.fi'

c This list of variable name strings must be kept in sync with the
c corresponding PARAMETER definitions in Nut.fi.

      DATA VARNAMES /
     >  'alpha',
     >  'b',
     >  'C0',
     >  'Cav',
     >  'Cmin',
     >  'De',
     >  'Dl',
     >  'dt',
     >  'f',
     >  'f_slope',
     >  'growth',
     >  'iter',
     >  'Kd',
     >  'Km',
     >  'L',
     >  'oalpha',
     >  'r0',
     >  'rho',
     >  'rx',
     >  'source',
     >  'steps',
     >  't',
     >  'theta',
     >  'U',
     >  'Unew',
     >  'v0',
     >  'Vmax',
     >  'w0',
     >  'Lv' /
      DATA UNITNAMES /
     >  'cm/s',         ! alpha
     >  '1',            ! b
     >  'mol/cm3',      ! C0
     >  'mol/cm3',      ! Cav
     >  'mol/cm3',      ! Cmin
     >  'cm2/s',        ! De
     >  'cm3/g',        ! Dl
     >  's',            ! dt
     >  '1',            ! f
     >  '1/s',          ! f_slope
     >  '1/s',          ! growth
     >  '1',            ! iter
     >  'cm3/g',        ! Kd
     >  'mol/cm3',      ! Km
     >  'cm',           ! L
     >  'cm/s',         ! oalpha
     >  'cm',           ! r0
     >  'g/cm3',        ! rho
     >  'cm',           ! rx
     >  'mol/cm3/s',    ! source
     >  '1',            ! steps
     >  's',            ! t
     >  '1',            ! theta
     >  'mol',          ! U
     >  'mol',          ! Unew
     >  'cm/s',         ! v0
     >  'mol/cm2/s',    ! Vmax
     >  'cm3/s',        ! w0
     >  '1/cm2'/        ! Lv  (cm/cm3)
      END

c *------------------------------------------------------------
c * Nut: The main program.
c *------------------------------------------------------------

      PROGRAM NUT
      IMPLICIT NONE
      INCLUDE 'Nut.fi'

c Functions used:
      LOGICAL READIN
      INTEGER READFIELDS, RTRIM, VARNUM, IMEMQ

c Internal variables
      LOGICAL EXISTP
      CHARACTER*60 CFNAME, RFNAME
      INTEGER I, ITER, STEPS, J, NI
      CHARACTER*80 COMMENT
      DOUBLE PRECISION VAL
      DOUBLE PRECISION f_init, theta_init

c For input switches file processing: NF = number of fields
      INTEGER NF, MAXF
      PARAMETER (MAXF=NVARS+1)
      CHARACTER*15 FIELDS(MAXF)

      TAB = CHAR(9)

      WRITE (STDOUT, 901)
 901    FORMAT (' Soil Nutrient Uptake Model'/
     >          ' Original program by Carol Zollweg'/
     >          ' Modifications by Jonathan Rees, August 1997'/)

c     *------------------------------------------------
c     * Read name of the control file (which holds both switches and inputs).
c     *------------------------------------------------
 10   WRITE (STDOUT, *) 'Please enter the name of the control file: '
      READ (STDIN, *) CFNAME
      J = RTRIM(CFNAME)
      CFNAME = CFNAME(1:j)
      OPEN (1, FILE=CFNAME, STATUS='OLD', ERR=15)
      GOTO 18
 15   WRITE (STDOUT, *) '!! Couldn''t find control file ', CFNAME
      GOTO 10

 18   WRITE (STDOUT, *) '...Reading control information and data from ',
     >                  CFNAME

c     *------------------------------------
c     * Open the results file
c     *------------------------------------
 20   WRITE (STDOUT, *) 'Please enter the name of the results file: '
      READ (STDIN, *) RFNAME
      J = RTRIM(RFNAME)
      RFNAME = RFNAME(1:J)

      INQUIRE (FILE=RFNAME, EXIST=EXISTP)
      IF (EXISTP) THEN
        WRITE (STDOUT, *) 'A file by that name already exists.  ',
     >                    'Replace it? (y/n)'
        READ (STDIN, *) COMMENT
        IF (COMMENT .EQ. 'y' .OR. COMMENT .EQ. 'yes' .OR.
     >      COMMENT .EQ. 'Y' .OR. COMMENT .EQ. 'YES') GOTO 23
        GOTO 20
      ENDIF

 23   OPEN (3, FILE=RFNAME, ERR=25)
      REWIND (3)      !relictual. should be unnecessary
      GOTO 28
 25   WRITE (STDOUT, *) '!! Couldn''t open results file ', RFNAME,
     >  ' for writing'
      GOTO 20

 28   WRITE (STDOUT, *) '...Writing results to ', RFNAME

c     *---------------------------------------------
c     * Process control information
c     *---------------------------------------------

c The first few lines are an identifying comment.  They are simply copied
c to the result file.
c     LOOP
 30     READ (1, '(A80)', END=777) COMMENT
        J = RTRIM(COMMENT)
        WRITE (3, 33) COMMENT(1:MAX(J,1))
 33       FORMAT (A)
        IF (J .GT. 0) GOTO 30
c     END LOOP

c Initialize the variable status table, etc.
      DO 45 I = 1, NVARS
        STATUS(I) = CALCULATED
        VALS(I) = 0.0
        SCALE(I) = 1.0
 45   CONTINUE
c Oalpha is a special case.
      STATUS(Voalpha) = FIXED
      NUMINVARS = 0

c Default values
      VALS(Vdt) = 24 * 60 * 60    !seconds in a day

c Subsequent lines are parameter status specifiers:
c   name units status value...
c     LOOP
 50     NF = READFIELDS(1, FIELDS, MAXF)
        IF (NF .LT. 0) GOTO 777
        IF (NF .EQ. 0) GOTO 59
        IF (NF .LT. 3)
     >    WRITE(STDOUT, *) '!! Invalid status line syntax -- ',
     >                     FIELDS(1)

c Fields are:   <name> <units> <status> [<parameters>]
        I = VARNUM(FIELDS(1))
        IF (I .EQ. 0) THEN
          WRITE (STDOUT, *) '!! Unrecognized variable name ',
     >                      FIELDS(1), ' in switches table'
          GOTO 50
        ENDIF

c Perform unit normalization
        IF (FIELDS(2) .EQ. 'day') THEN
          FIELDS(2) = 's'
          SCALE(I) = 1.0/(24.0 * 60.0 * 60.0)
        ELSE IF (FIELDS(2) .EQ. '1/day') THEN
          FIELDS(2) = '1/s'
          SCALE(I) = (24.0 * 60.0 * 60.0)
        ELSE IF (FIELDS(2) .EQ. 'm') THEN
          FIELDS(2) = 'cm'
          SCALE(I) = 0.01
        ENDIF

c Check unit against permissible units
        IF (FIELDS(2) .NE. UNITNAMES(I)) THEN
          WRITE (STDOUT, *) '!! Unrecognized unit ', FIELDS(2),
     >                      ' for variable ', FIELDS(1)
          WRITE (STDOUT, *) '!! (Expected e.g. ', UNITNAMES(I), ')'
        ENDIF

c Read the initial value and scale it accordingly
        IF (NF .LT. 4 .OR. RTRIM(FIELDS(4)) .EQ. 0) THEN
          VALS(I) = 0.0
        ELSE
          READ(FIELDS(4), *) VAL
          VALS(I) = VAL / SCALE(I)
        ENDIF

c The status field is "input", "fixed", "calculated", etc.
        IF (FIELDS(3) .EQ. 'input') THEN
          STATUS(I) = INPUT
          NUMINVARS = NUMINVARS + 1
          IF (NF .GT. 3)
     >      WRITE (STDOUT, *) '!! Extraneous parameter for ',
     >        FIELDS(1)
        ELSE IF (FIELDS(3) .EQ. 'calculated') THEN
          STATUS(I) = CALCULATED
        ELSE
          IF (NF .LT. 4) WRITE (STDOUT, *) '!! Missing parameter for ',
     >      FIELDS(1)
          IF (FIELDS(3) .EQ. 'fixed') THEN
            STATUS(I) = FIXED
          ELSE IF (FIELDS(3) .EQ. 'linear') THEN
            STATUS(I) = LINEAR
          ELSE IF (FIELDS(3) .EQ. 'exponential') THEN
            STATUS(I) = EXPONENTIAL
          ELSE
            WRITE (STDOUT, *) '!! Invalid status: ', FIELDS(3)
            WRITE (STDOUT, *) '!! (Wanted input, fixed, calculated, ',
     >                        'linear, or exponential)'
          ENDIF
        ENDIF
        GOTO 50
c     END LOOP

c Blank line encountered.
 59   CONTINUE

c     *--------------------------------
c     * Process the column names line.  This determines what will be
c     * written to the results file, and also specifies the order of
c     * the fields in the input records.
c     *--------------------------------

      NCOLS = READFIELDS(1, FIELDS, MAXF)
      NI = 0
      DO 60 J = 1, NCOLS
        I = VARNUM(FIELDS(J))
        IF (I .EQ. 0) THEN
          WRITE (STDOUT, *) '!! Unrecognized variable name in inputs ',
     >                      'file: [', FIELDS(J), ']'
        ELSE IF (STATUS(I) .EQ. INPUT) THEN
          NI = NI + 1
          INVARS(NI) = I
        ENDIF
        COLVARS(J) = I
 60   CONTINUE

      IF (NI .NE. NUMINVARS) THEN
        WRITE (STDOUT, *) '!! Wrong number of inputs: switches list',
     >        ' has', NUMINVARS, ', but column headings list has ', NI
      ENDIF

c The results file should have at least one column, yes?
      IF (NCOLS .EQ. 0) THEN
        WRITE (STDOUT, *) '!! Warning: no output columns specified!'
        GOTO 98
      ENDIF

c The results file should have an uptake column, yes?
      IF (IMEMQ(VU, COLVARS, NCOLS) .EQ. 0) THEN
        WRITE (STDOUT, *) '!! Warning: no column specified for ',
     >    'uptake (U)'
      ENDIF

      if(debugp) then
        print *, '| columns: ',
     >    (varnames(colvars(j)), j = 1, ncols)
        if (numinvars .gt. 0) then
          print *, '| inputs:  ',
     >      (varnames(invars(j)), j = 1, numinvars)
        else
          print *, '| no inputs'
        endif
      endif

c     *---------------------------------
c     * Initialize, etc.
c     *---------------------------------

c Begin writing the results file (comment has already been written, above)

      DO 70 I = 1, NCOLS
        FIELDS(I) = VARNAMES(COLVARS(I))
 70   CONTINUE
      CALL WRITEFIELDS(3, FIELDS, NCOLS)

c Compute format to be used in output routine
      J = 1
      FORM(1:J) = '('
      DO 71 I = J, (NCOLS-1)*10 + J, 10
        FORM(I+1:I+10) = 'g12.6,1h-,'
        FORM(I+9:I+9) = TAB
   71 CONTINUE
      FORMLEN = NCOLS*10 + J + 1
      FORM(FORMLEN:FORMLEN) = ')'

      if(debugp)
     >  print *, '| Format = [', FORM(1:FORMLEN), ']'

c  Sanity check
      IF (NUMINVARS .LE. 0 .AND. STATUS(Vsteps) .NE. FIXED)
     >  WRITE (STDOUT, *) '!! No termination condition ',
     >                    '(no inputs, and steps isn''t fixed)'

c     *-----------------------------------------------------------------
c     * Begin iterating
c     *-----------------------------------------------------------------

c  The following are conversions to integer
      ITER = 1
      STEPS = VALS(Vsteps)
      if(debugp) print *, '| init done'

c     BEGIN LOOP
 90   CONTINUE

c  Read in the ith parameter set (if there is one)
        IF (NUMINVARS .GT. 0) THEN
          IF (.NOT. READIN(1)) GOTO 98
        ENDIF

        VALS(Viter) = ITER
        if(debugp) print *, 'Iteration #', ITER

c  Perform one time step (time t to time t+dt) of uptake model
        CALL RUNMODEL

c  Get initial values of f and theta to establish linear relationship
        IF (ITER .EQ. 1) THEN
          f_init = VALS(Vf)
          theta_init = VALS(Vtheta)
        ENDIF

c  Write results
        WRITE (3, FORM) (VALS(COLVARS(I)) * SCALE(COLVARS(I)),
     >                   I = 1, NCOLS)

c  Stop if it's time to stop
        IF (STATUS(Vsteps) .EQ. FIXED .AND. ITER .GE. STEPS)
     >    GOTO 98

c  Prepare for next iteration
        CALL UPDATE(f_init, theta_init)
        VALS(Vt) = VALS(Vt) + VALS(Vdt)
        ITER = ITER + 1

        GOTO 90
c     END LOOP

c     *-----------------------------------------------------------------
c     * End iterating
c     *-----------------------------------------------------------------

c  Close control and results files
 98   CLOSE(1)
      CLOSE(3)
      WRITE (STDOUT, *) 'Simulation completed.'
      WRITE (STDOUT, *) 'Ran', ITER, ' iterations.'

      GOTO 778

 777  WRITE (STDOUT,*) '!! Not enough stuff in control file'

 778  CONTINUE
      END

c *---------------------------------------------------------------------
c * RUNMODEL oversees execution of one iteration of the uptake model.
c *---------------------------------------------------------------------

      SUBROUTINE RUNMODEL
      IMPLICIT NONE
      INCLUDE 'Nut.fi'

c Set up inputs to the uptake model.
      alpha     = VALS(Valpha)
      b         = VALS(Vb)
      Cav       = VALS(VCav)
      De        = VALS(VDe)
      Dl        = VALS(VDl)
      dt        = VALS(Vdt)
      f         = VALS(Vf)
      growth    = VALS(Vgrowth)
      Km        = VALS(VKm)
      L         = VALS(VL)
      r0        = VALS(Vr0)
      rho       = VALS(Vrho)
      rx        = VALS(Vrx)
      source    = VALS(Vsource)
      theta     = VALS(Vtheta)
      Unew      = VALS(VUnew)
      v0        = VALS(Vv0)
      Vmax      = VALS(VVmax)
      w0        = VALS(Vw0)

c Obsolescent?
      Cmin      = VALS(VCmin)    ! Cmin (mol/cm3)

      if(debugp) print *, '| L:', L
      if(debugp) print *, '| r0:', r0
      if(debugp) print *, '| c0:', C0

c MODEL does the real work of computing values of dependent variables.
      CALL MODEL

c Copy calculated values back into the VALS vector.
      VALS(Valpha)  = alpha
      VALS(Vb)      = b
      VALS(VCav)    = Cav
      VALS(VC0)     = C0
      VALS(VDe)     = De
      VALS(VL)      = L
      VALS(Vrx)     = rx
      VALS(Vv0)     = v0
      VALS(VU)      = U
      VALS(VUnew)   = Unew

      RETURN
      END

c=======================================================================
c Derived from OurPrism2 file UpTk.f

c *------------------------------------------------------------
c * MODEL
c * Do all computations for one step of the uptake model.
c *
c * Calculate the root length for both new and established roots, the
c * root density, and the average radial distance to the next root's
c * zone of influence.  The type of root used depends on the soil
c * horizon.
c * The root length (cm) = m/gC * gC for the root * 100cm/m
c * The root density (m/m**3) = m/gC * gC for the root / soil volume
c * The average radial distance to the next root's zone of
c * influence, rx, (cm) =
c *      1/sqrt(Pi * RootDens (m/m**3)/10000cm**2/m**2)
c *---------------------------------------------------------------------

      SUBROUTINE MODEL
      IMPLICIT NONE

c Common variables
      INCLUDE 'Nut.fi'
c Function
      DOUBLE PRECISION CALCPc
c Internal variables
      DOUBLE PRECISION Deb, A, rav
      DOUBLE PRECISION GAMMA
      DOUBLE PRECISION DeltaL, Pc, RTRATIO
      DOUBLE PRECISION UEST
c      DOUBLE PRECISION Lv
c Temporaries
      DOUBLE PRECISION AA, AB, AC, AD, AE, AF, FOO

c EVENSMLR is important.  v0 goes down as low as 10^-11 (for Robinia
c  pseudoacacia), so the threshold should probably be smaller than that.
c  If v0 is smaller than that but nonzero, we should probably
c  generate a warning, because user should be informed that the
c  the diffusion-only formula will be invoked.
      DOUBLE PRECISION CM3_L, EVENSMLR
        PARAMETER (CM3_L = 1000.0, EVENSMLR = 1E-12)

c     *-----------------------------------------------------------*
c     * If there are no roots, there is no uptake!  Exit routine. *
c     *-----------------------------------------------------------*
      IF (L .LE. 0.0d0) THEN
        WRITE (STDOUT,*) '!! Root length <= 0'
        RETURN
      ENDIF

c     *-----------------------------------------------------------------
c     * Compute b if requested.
c     *-----------------------------------------------------------------
      IF (STATUS(Vb) .EQ. CALCULATED) THEN
        b = theta + rho*Kd
      ENDIF

c     *-----------------------------------------------------------------
c     * Compute De if requested.
c     *-----------------------------------------------------------------
      IF (STATUS(VDe) .EQ. CALCULATED) THEN
        Deb = Dl * theta * f
        De = Deb / b
      ELSE
        Deb = De * b
      ENDIF

c     *-----------------------------------------------------------------
c     * Compute rx (zone of influence), if required.
c     * Lv is root length per unit volume (cm of root per cm^3 of soil).
c     *-----------------------------------------------------------------
      IF (STATUS(Vrx) .EQ. CALCULATED) THEN
        rx = 1.0d0 / SQRT(PI * Lv)
      ENDIF

c     *---------------------------------
c     * Calculate the root ratio rx/r0.
c     *---------------------------------
      RTRATIO = rx / r0
      if(debugp) print *, '| RTRATIO:', RTRATIO

c     *------------------------------------------------------------
c     * Calculate v0 and gamma.  v0 equals
c     * the water uptake w0 (cm3/s) divided by the root surface
c     * area (cm2).  Don't allow the inward velocity to equal zero.
c     * Root surface area calculated with established roots
c     * only to make it consistent with root water uptake w0.
c     * Root water uptake is calculated with established roots.
c     *------------------------------------------------------------
      IF (STATUS(Vv0) .EQ. CALCULATED) THEN
        v0 = MAX(EVENSMLR, w0 / (2 * PI * r0 * L))
        if(debugp) print *, '| v0 (calculated):', v0
      ENDIF
      gamma = (r0 * v0)/Deb
      if(debugp) print *, '| gamma:', gamma

      IF (STATUS(Valpha) .EQ. CALCULATED) THEN
c     *-----------------------------------------------------------------
c     * Ruth's method for calculating root absorbing power (alpha).
c     *-----------------------------------------------------------------
c       -- CALCALPHA computes alpha and C0, but not Pc.
        CALL CALCALPHA(gamma, RTRATIO)
        if(debugp) print *, '| C0 (direct):', C0
        if(debugp) print *, '| alpha (direct):', alpha
        Pc = C0 / Cav
        if(debugp) print *, '| Pc (direct):', Pc
      ENDIF

      IF (STATUS(Voalpha) .EQ. CALCULATED) THEN
c     *-----------------------------------------------------------------
c     * Calculate the root absorbing power (alpha) as a function of
c     * Vmax and Km by iterating because alpha is a function of the
c     * concentration at the root surface and that concentration is a
c     * function of alpha.
c     *
c     * N.b. UPALPHA calculates C0 and Pc in addition to oalpha.
c     *-----------------------------------------------------------------
        CALL UPALPHA(EVENSMLR,
     >               0.1d0,      ! OLDC0 (OLDCSURF)
     >               Vmax, Km, Cmin, gamma,
     >               v0, RTRATIO, r0, Deb, b, Cav,
     >               .TRUE.,
     >               alpha,      ! This parameter appears to be ignored
     >               oalpha, C0, Pc, debugp)
        if(debugp) print *, '| C0 (iterative):', C0
        if(debugp) print *, '| alpha (iterative):', oalpha
        if(debugp) print *, '| Pc (iterative):', Pc
        alpha = oalpha
      ENDIF

c     *-----------------------------------------------------------------
c     * If alpha is calculated, then Pc and C0 are calculated as a side effect.
c     * If alpha is fixed, then we have to compute them.
c     *-----------------------------------------------------------------
      IF (STATUS(Valpha) .NE. CALCULATED .AND.
     >    STATUS(Voalpha) .NE. CALCULATED) THEN
c       -- If alpha is fixed, use it to compute Pc and C0.
        Pc = CALCPc(alpha, gamma, RTRATIO, Deb, EVENSMLR)
        C0 = Pc * Cav
      ENDIF

c     *-----------------------------------------------------------------
c     * Calculate rav = the radius where the concentration equals
c     * the average concentration. If v0 is close to zero, use a
c     * different equation.
c     *-----------------------------------------------------------------
      IF (v0 .LE. EVENSMLR) THEN
        AA = Deb/(alpha * r0)
        AB = (1 - Pc) / Pc
        FOO = EXP(AA * AB)
      ELSE
        AA = 2 / (2 - GAMMA)
        AC = -1 / GAMMA
        AD = RTRATIO**(2 - GAMMA) - 1
        AE = (RTRATIO * RTRATIO) - 1
        FOO = ((AA * (AD/AE))**AC)
      ENDIF
      rav = r0 * FOO
      if (debugp) print *, '| rav:', rav

c     *------------------------------------------------------------
c     * Calculate the uptake due to the established roots (U_est).
c     *------------------------------------------------------------
      Uest = L * 2 * PI * r0 * alpha * C0 * dt
      if (debugp) print *, '| Uest:', Uest

c     *-----------------------------------------------------------
c     * Calculate the uptake by the new roots (if any).
c     * To do this, first calculate A = the amount of solute removed
c     * from the depletion zone per unit length of new root.  If v0
c     * is very close to zero, use a different equation.
c     *-----------------------------------------------------
      IF (STATUS(VUnew) .EQ. CALCULATED) THEN
        IF (v0 .LE. EVENSMLR) THEN
          AB = r0 * r0
          AC = rav * rav
          AD = AC - AB
          AE = Pc * ALPHA * r0 / Deb
          AF = LOG(rav / r0)
          A = PI * b * Cav *
     >         (((1 - Pc) * AD) - (AE * ((AC * (AF - 0.5)) + AB*0.5)))
        ELSE
          AA = 1 - Pc*alpha/v0
          AB = rav*rav - r0*r0
          AC = 2*Pc/v0
          AD = (v0 - ALPHA)/(2 - GAMMA)
          AE = ((rav/r0) ** (2 - GAMMA)) - 1
          A = PI * b * Cav *
     >         (AA*AB - (AC * AD * r0*r0 * AE))
        ENDIF
        if (debugp) print *, '| A:', A

c     *-----------------------------------------------------------
c     * Length of roots that are brand new.
c     *-----------------------------------------------------------
        DeltaL = L * (growth * dt)

c     *-----------------------------------------------------------
c     * Don't allow uptake by new roots to be negative.
c     **** Ruth questions this.  Maybe it _should_ be allowed to go negative.
c     *-----------------------------------------------------------
        Unew = MAX(A * DeltaL, 0.0d0)
        if (debugp) print *, '| Unew:', Unew
      ENDIF

c     *---------------------------------------------------------
c     * The total amount of nutrient taken up by the root is
c     * the amount up by the established roots plus the amount
c     * taken up by the new roots.
c     *---------------------------------------------------------
      U = Uest + Unew
      if (debugp) print *, '| U:', U

      RETURN
      END

c=======================================================================
c Alpha computation copied mindlessly from Ruth's handwritten notes.
c
c Inputs:  gamma, RTRATIO
c          v0, Cav, Km, Vmax
c Outputs: alpha, C0

      SUBROUTINE CALCALPHA(gamma, RTRATIO)
      IMPLICIT NONE
      DOUBLE PRECISION gamma, RTRATIO

c Common variables.
      INCLUDE 'Nut.fi'

c Internal variables.
      DOUBLE PRECISION TEMP1, TEMP2
      DOUBLE PRECISION AA, BB
      DOUBLE PRECISION ROOT1, ROOT2

c Compute intermediate quantities.
        TEMP1 = RTRATIO**(2 - GAMMA) - 1
        TEMP2 = (RTRATIO * RTRATIO) - 1
      AA = (2 / (2 - GAMMA)) * (TEMP1 / TEMP2)
      if(debugp) print *, '| AA:', AA

      BB = v0 * Cav
      if(debugp) print *, '| BB:', BB

      CALL QUADRATIC(AA * v0,
     >               AA*Km*v0 + (1-AA)*Vmax - BB,
     >               BB * Km,
     >               root1, root2)
      if(debugp) print *, '| roots are', root1, root2

      C0 = MAX(root1, root2)
      IF (C0 .LT. 0) THEN
        WRITE (STDOUT, *) '!! negative roots in C0 computation'
      ENDIF

      alpha = Vmax / (Km + C0)

      RETURN
      END

c=======================================================================
c Update state variables for next iteration.

      SUBROUTINE UPDATE(f_init, theta_init)
      IMPLICIT NONE
      DOUBLE PRECISION f_init, theta_init

c Common and internal variables.
      INCLUDE 'Nut.fi'

c Formula for Cav copied from UpComp.f
      IF (STATUS(VCav) .EQ. CALCULATED) THEN
        Cav = Cav - U / (L * PI * (rx * rx) * b)
      ENDIF

      IF (STATUS(VCav) .EQ. LINEAR)
     >  Cav = Cav + (source * dt)

c Root growth can be linear or exponential
      IF (STATUS(VL) .EQ. LINEAR) THEN
        L = L + (growth * dt)
      ELSE IF (STATUS(VL) .EQ. EXPONENTIAL) THEN
        L = L * (1.0d0 + (growth * dt))
      ENDIF

c Update fugacity, which is linear in theta
      IF (STATUS(Vf) .EQ. CALCULATED) THEN
        f = f_init + f_slope*(theta - theta_init)
      ENDIF

      VALS(VCav) = Cav
      VALS(VL) = L
      VALS(Vf) = f

      END

c=======================================================================
c     *----------------------------------------------------------*
c     * Calculate the ratio Pc = C0/Cav.
c     * If v0 <= EvenSmlr a different calculation for Pc is
c     * used based on the diffusion-only case of Baldwin et al.
c     * 1973 (eqn III).  This equation is independent of  v0
c     * and gamma.
c     *----------------------------------------------------------*
      DOUBLE PRECISION FUNCTION
     >  CALCPc(myalpha, gamma, RTRATIO, Deb, EVENSMLR)
      IMPLICIT NONE
      DOUBLE PRECISION myalpha, gamma, RTRATIO, Deb, EVENSMLR
      INCLUDE 'Nut.fi'

      DOUBLE PRECISION Pc, AA, AB, AC, AD, AE

      IF (v0 .LE. EVENSMLR) THEN
c       rx = RTRATIO * r0      -- this has already been computed
        AA = myalpha * r0 / (2 * Deb)
        AB = rx * rx
        AC = AB * myalpha * r0 / Deb
        AD = AB - r0 * r0
        AE = LOG(rx/r0)
        Pc = 1/(1 - AA + AC/AD * AE)
        if(debugp) print *, '| Pc (small v0): ', Pc
      ELSE
        AA = v0 - myalpha
        AB = 2 / (2 - gamma)
        AC = RTRATIO ** (2 - gamma)
        AD = RTRATIO * RTRATIO
        AE = (AC - 1.0) / (AD - 1.0)
        Pc = v0 / (myalpha + AA * AB * AE)
        if(debugp) print *, '| Pc (normal): ', Pc
      ENDIF

      CALCPc = Pc
      RETURN
      END

c=======================================================================
c *-----------------------------------
c * Quadratic equation solver.
c * A numerical analyst would probably tell you that this is absolutely
c * the worst way to solve a quadratic equation, due to its
c * vulnerability to overflow, underflow, and roundoff errors.  A better
c * way might be found in a reference such as "Numerical Recipes," which I
c * haven't yet consulted.
c *-----------------------------------

      SUBROUTINE QUADRATIC(A, B, C, R1, R2)
      IMPLICIT NONE
      DOUBLE PRECISION A, B, C, R1, R2
      DOUBLE PRECISION DISCRIM, S

      DISCRIM = ((B * B) - (4 * A * C))
      IF (DISCRIM .LT. 0) THEN
        PRINT *, '!! Quadratic equation has complex roots'
        R1 = 0
        R2 = 0
        RETURN
      ENDIF
      S = SQRT(DISCRIM)
      R1 = (-B + S)/(2*A)
      R2 = (-B - S)/(2*A)
      RETURN
      END

c=======================================================================
c I/O utilities.

c *-----------------------------------
c * Read a line of input parameters.
c * Returns .true. normally, .false. at end of file.
c *-----------------------------------

      LOGICAL FUNCTION READIN(UNIT)
      IMPLICIT NONE
      INTEGER UNIT

      INCLUDE 'Nut.fi'
      INTEGER I
      DOUBLE PRECISION INPUTS(NVARS)

      READ (UNIT, *, END=777) (INPUTS(I), I=1,NUMINVARS)
      DO 20 I = 1, NUMINVARS
        VALS(INVARS(I)) = INPUTS(I)
  20  CONTINUE
      READIN = .TRUE.
      RETURN
c
 777  WRITE (STDOUT, *) 'End of input records'
      READIN = .FALSE.
      END

c *-------------------------------------------
c * Separate a line into tab-separated fields.
c * MAXF = maximum number of fields (i.e. size of output array).
c * Returns the number of fields read, or -1 on EOF.
c *-------------------------------------------
      INTEGER FUNCTION READFIELDS(UNIT, FIELDS, MAXF)
      IMPLICIT NONE
      INTEGER UNIT, MAXF
      CHARACTER*(*) FIELDS(MAXF)

      CHARACTER*320 BUFFER
      CHARACTER*1 TAB
      INTEGER POS, J, NF, FIN, RTRIM
c
      TAB = CHAR(9)
c
      POS = 1
      NF = 0
      READ (UNIT, '(A150)', END=9) BUFFER
      FIN = RTRIM(BUFFER)
 1    IF (POS .LE. FIN .AND. NF .LT. MAXF) THEN
        J = INDEX(BUFFER(POS:FIN), TAB)
        IF (J .GT. 0) THEN
           NF = NF + 1
           FIELDS(NF) = BUFFER(POS:POS+J-2)
           POS = POS + J
           GOTO 1
        ELSE
           NF = NF + 1
           FIELDS(NF) = BUFFER(POS:FIN)
        ENDIF
      ENDIF

      READFIELDS = NF
      RETURN

c end of file
 9    READFIELDS = -1
      RETURN
      END

c *---------------------------------------------------------------------
c * Write NF fields, separated by tabs.
c * [Trailing whitespace isn't written.  -- yes it is]
c *---------------------------------------------------------------------

      SUBROUTINE WRITEFIELDS(UNIT, FIELDS, NF)
      IMPLICIT NONE
      INTEGER UNIT, NF
      CHARACTER*15 FIELDS(NF)

      CHARACTER*320 BUFFER
      INTEGER POS, F, I
c      INTEGER RTRIM
      CHARACTER*1 TAB
c
      TAB = CHAR(9)
      POS = 1
      DO 10 F = 1, NF
c        I = RTRIM(FIELDS(F))
        I = 15     !TEMPORARY CHEAT, TRIAL PERIOD ONLY
        BUFFER(POS:POS+I-1) = FIELDS(F)(1:I)
           POS = POS + I
           IF (F < NF) THEN
              BUFFER(POS:POS) = TAB
              POS = POS + 1
           ENDIF
 10     CONTINUE
      WRITE (UNIT, 11) BUFFER(1:POS-1)
 11     FORMAT (A)
      RETURN
      END

c *---------------------------------------------------------------------
c * Find position of rightmost non-whitespace character in string
c *---------------------------------------------------------------------
      FUNCTION RTRIM(S)
      IMPLICIT NONE
      CHARACTER*(*) S
      INTEGER POS, RTRIM
      POS = LEN(S)
c     LOOP
 10     IF (S(POS:POS) .NE. ' ') THEN
          RTRIM = POS
          RETURN
        ENDIF
        POS = POS - 1
        IF (POS .LE. 0) THEN
          RTRIM = 0
          RETURN
        ENDIF
        GOTO 10
c     END LOOP

      END

c *---------------------------------------------------------------------
c * Determine whether a given integer occurs in a vector.
c * (MEMQ is a similar function in the Lisp programming language.)
c *---------------------------------------------------------------------
      INTEGER FUNCTION IMEMQ(N, V, LEN)
      INTEGER N, LEN, V(LEN)
      INTEGER I
      DO 1 I = 1, LEN
        IF (V(I) .EQ. N) THEN
          IMEMQ = I
          RETURN
        ENDIF
    1 CONTINUE
      IMEMQ = 0
      RETURN
      END

c *------------------------------------------------------------
c * Find the index of a variable given its name.
c *------------------------------------------------------------
      INTEGER FUNCTION VARNUM(NAME)
      IMPLICIT NONE
      INCLUDE 'Nut.fi'
      CHARACTER*(*) NAME
      INTEGER I

      DO 1 I = 1, NVARS
        IF (NAME .EQ. VARNAMES(I)) THEN
          VARNUM = I
          RETURN
        ENDIF
    1 CONTINUE
      VARNUM = 0
      RETURN
      END

c=======================================================================
c From Carol's file UpAlpha.f

c Old iterative method for computing alpha.

c INPUTS:  EVENSMLR, INC0, Vmax, KM, CMIN, GAMMA, v0,
c          RTRATIO, r0, Deb, B, Cav,
c          CALC_alpha,                  -- flag
c          INALPHA,      -- observed only if CALC_alpha is .FALSE.
c
c OUTPUTS:
c          ALPHA, C0, Pc

c Pc = The ratio between the concentration of the
c substance at the root surface and the
c average solution concentration.

      SUBROUTINE UPALPHA(EVENSMLR,   INC0,    Vmax,
     >                   Km,         Cmin,     GAMMA,
     >                   v0,         RTRATIO,  r0,
     >                   Deb,        b,        Cav,
     >                   CALC_ALPHA, INALPHA,  ALPHA,
     >                   C0,         Pc,       DEBUGP)
      IMPLICIT NONE
c Parameters
      DOUBLE PRECISION EVENSMLR, INC0, Vmax
      DOUBLE PRECISION Km, Cmin, GAMMA
      DOUBLE PRECISION RTRATIO, r0, Deb, b, Cav, v0
      LOGICAL   CALC_ALPHA
      DOUBLE PRECISION INALPHA
      DOUBLE PRECISION ALPHA, C0, Pc
      LOGICAL DEBUGP

c Function
      DOUBLE PRECISION CALCPc

c Internal variables
      DOUBLE PRECISION OLDC0
      INTEGER ITER
      LOGICAL DONE

c Constants
      DOUBLE PRECISION CRITERION
        PARAMETER (CRITERION = 0.001d0)
      INTEGER MAXITER
        PARAMETER (MAXITER = 100)

c Initialize
      ITER = 0
      OLDC0 = INC0

c     *--------------------------------------------------------------*
c     * If alpha is to be calculated, then check if concentration is *
c     * below minimum indicating that the root absorbing power is    *
c     * zero.  If it is, indicate that the routine converged in zero *
c     * iterations.  Otherwise, indicate that the routine has not    *
c     * converged.                                                   *
c     *--------------------------------------------------------------*
      IF (CALC_alpha) THEN
        IF (Cav .LE. Cmin) THEN
          ALPHA = 0.0
          DONE = .TRUE.
        ELSE
          DONE = .FALSE.
        ENDIF

c     *---------------------------------------------------------*
c     * If alpha is not to be calculated make sure that routine *
c     * does not iterate.                                       *
c     *---------------------------------------------------------*
      ELSE
        DONE = .TRUE.
      ENDIF

c     LOOP
 10     CONTINUE
c     *-----------------------------------------------*
c     * Increment the number of iterations completed. *
c     *-----------------------------------------------*
        ITER = ITER + 1
c     *-------------------------------------------------------*
c     * If alpha is to be calculated then calculate the root  *
c     * absorbing power based on the old value of             *
c     * concentration at the root surface. Otherwise, set it  *
c     * to the input value.   Vmax / (Km + c0)                *
c     *-------------------------------------------------------*
        IF (CALC_ALPHA) THEN
           ALPHA = Vmax / (Km + OLDC0)
        ELSE
           ALPHA = INALPHA
        ENDIF
        if(debugp) print *, '| alpha in loop: ', alpha

        Pc = CALCPc(alpha, gamma, RTRATIO, Deb, EVENSMLR)

c     *----------------------------------------------------------*
c     * Calculate the new concentration of the substance at the  *
c     * root surface: C_0 = P_c * C_av
c     *----------------------------------------------------------*
        C0 = Pc * Cav
c     *----------------------------------------------------*
c     * Is this new concentration close enough to the old  *
c     * concentration to consider the solution converged?  *
c     * (Only check if alpha is being calculated)          *
c     *----------------------------------------------------*
        IF (CALC_ALPHA) THEN
          if(debugp) print *, '| criterion: ', criterion
          IF (ABS((C0 - OLDC0) / (OLDC0 + Km)) .LE.
     >          CRITERION) THEN
            DONE = .TRUE.
          ELSE
            OLDC0 = C0
          ENDIF
        ENDIF

c     *---------------------------------------------------------*
c     * If the maximum number of iterations has been run, allow *
c     * routine to end.                                         *
c     *---------------------------------------------------------*
        IF (ITER .GE. MAXITER) DONE = .TRUE.

        IF (.NOT.DONE) GOTO 10
c     END LOOP

      RETURN
      END


text:
module nf95_inquire_dimension_m

  implicit none

contains

  subroutine nf95_inquire_dimension(ncid, dimid, name, nclen, ncerr)

    use handle_err_m, only: handle_err
    use netcdf, only: nf90_inquire_dimension

    integer,                       intent( in) :: ncid, dimid
    character (len = *), optional, intent(out) :: name
    integer,             optional, intent(out) :: nclen
    integer, intent(out), optional:: ncerr

    ! Variable local to the procedure:
    integer ncerr_not_opt

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

    ncerr_not_opt = nf90_inquire_dimension(ncid, dimid, name, nclen)
    if (present(ncerr)) then
       ncerr = ncerr_not_opt
    else
       call handle_err("nf95_inquire_dimension", ncerr_not_opt, ncid)
    end if

  end subroutine nf95_inquire_dimension

end module nf95_inquire_dimension_m

max_stars_repo_name: noaa-ocs-modeling/PaHM

max_issues_repo_name: noaa-ocs-modeling/PaHM


text:
+define+NETLIST+SDFSCAN
+notimingchecks
+nowarnTFNPC
+nowarnIWFA
+nowarnSVTL
src/timescale.v
src/UPA2_test.sv
//
// Uncomment for TSMC 180nm
-v /classes/eeee720/maieee/lib/tsmc-0.18/verilog/tsmc18.v
//
// Uncomment for TSMC 65nm
// -v /classes/ee620/maieee/lib/synopsys/TSMC_tcbc65/TSMCHOME/digital/Front_End/verilog/tcbn65lp_200a/tcbn65lp.v
//
// Uncomment for SAED 90nm
// -y /classes/ee620/maieee/lib/synopsys/SAED_EDK90nm/Digital_Standard_Cell_Library/verilog/
//
// Uncomment for SAED 32nm
// -v /classes/ee620/maieee/lib/synopsys/SAED_EDK32-28nm/SAED32_EDK/lib/stdcell_rvt/verilog/saed32nm.v
//
+libext+.vs
+librescan
//
// Uncomment for TSMC 180nm
-y netlist/UPA2_tsmc18_scan/
//
// Uncomment for TSMC 65nm
// -y netlist/UPA2_tsmc065_scan/
//
// Uncomment for SAED 90nm
// -y netlist/UPA2_saed90nm_scan/
//
// Uncomment for SAED 32nm
// -y netlist/UPA2_saed32nm_scan/

ext: f


text:
!
!     (c) 2019-2020 Guide Star Engineering, LLC
!     This Software was developed for the US Nuclear Regulatory Commission (US NRC) under contract
!     "Multi-Dimensional Physics Implementation into Fuel Analysis under Steady-state and Transients (FAST)",
!     contract # NRC-HQ-60-17-C-0007
!
module differentiable_field_interface
  !! author: Damian Rouson
  !! date: 12/19/2019
  !!
  !! A differentiable field abstraction, deferring function evaluation and differential operator implementations to child classes
  use grid_interface, only : grid
  implicit none

  private
  public :: differentiable_field

  type, abstract :: differentiable_field
  contains
    procedure(field_interface), deferred :: evaluate
    procedure(field_interface), deferred :: laplacian
  end type

  abstract interface
    function field_interface(this, grid_points) result(f)
      !! Compute function values sampled at the provided locations
      import differentiable_field, grid
      implicit none
      class(differentiable_field), intent(in)  :: this
      class(grid), intent(in)  :: grid_points
      class(grid), allocatable :: f
    end function
  end interface

end module differentiable_field_interface

max_stars_repo_name: sourceryinstitute/MORFEUS-Source

max_issues_repo_name: sourceryinstitute/MORFEUS-Source

max_forks_repo_name: sourceryinstitute/MORFEUS-Source


text:
      FUNCTION ei(x)
      INTEGER MAXIT
      REAL ei,x,EPS,EULER,FPMIN
      PARAMETER (EPS=6.e-8,EULER=.57721566,MAXIT=100,FPMIN=1.e-30)
      INTEGER k
      REAL fact,prev,sum,term
      if(x.le.0.) pause 'bad argument in ei'
      if(x.lt.FPMIN)then
        ei=log(x)+EULER
      else if(x.le.-log(EPS))then
        sum=0.
        fact=1.
        do 11 k=1,MAXIT
          fact=fact*x/k
          term=fact/k
          sum=sum+term
          if(term.lt.EPS*sum)goto 1
11      continue
        pause 'series failed in ei'
1       ei=sum+log(x)+EULER
      else
        sum=0.
        term=1.
        do 12 k=1,MAXIT
          prev=term
          term=term*k/x
          if(term.lt.EPS)goto 2
          if(term.lt.prev)then
            sum=sum+term
          else
            sum=sum-prev
            goto 2
          endif
12      continue
2       ei=exp(x)*(1.+sum)/x
      endif
      return
      END

ext: for

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


text:
version https://git-lfs.github.com/spec/v1
oid sha256:e5a2017a66d306346ab781b100253e52087b35cd8f4ec700e3244a0320370c0d
size 112251

max_stars_repo_name: initialz/tensorflow-yolo-face-ios

max_issues_repo_name: initialz/tensorflow-yolo-face-ios

max_forks_repo_name: initialz/tensorflow-yolo-face-ios


submodule(grid_interface) grid_implementation
  use assertions_interface, only : assert, assertions
  implicit none

contains

    module function get_dims(this) result(dims)
        class(grid_t) :: this
        integer :: dims(3)

        dims(1) = this%ime - this%ims + 1
        dims(2) = this%kme - this%kms + 1
        dims(3) = this%jme - this%jms + 1
    end function

    ! note currently this is supplied in the domain object, however it should be moved to the grid object along with associated my_ny, etc routines
    ! anything that defines the domain grid should be in here...

    !> -------------------------------
    !! Decompose the domain into as even a set of tiles as possible in two dimensions
    !!
    !! Searches through possible numbers of x and y tiles that multiple evenly to
    !! give the total number of images requested.
    !!
    !! For each x/y split compute the number of grid cells in both dimensions in each tile
    !! return the split that provides the closest match between the number of x and y grid cells
    !!
    !! -------------------------------
    ! module subroutine domain_decomposition(this, nx, ny, nimages)
    !     class(domain_t), intent(inout) :: this
    !     integer,         intent(in)    :: nx, ny, nimages
    !     integer :: ysplit, xsplit, xs, ys, i
    !     real :: best, current, x, y
    !
    !     xsplit = 1
    !     ysplit = nimages
    !     xs = xsplit
    !     ys = ysplit
    !
    !     x = (nx/real(xsplit))
    !     y = (ny/real(ysplit))
    !
    !     if (y > x) then
    !         best = abs(1 - ( y / x ))
    !     else
    !         best = abs(1 - ( x / y ))
    !     endif
    !
    !     do i=nimages,1,-1
    !         if (mod(nimages,i)==0) then
    !             ysplit = i
    !             xsplit = nimages / i
    !
    !             x = (nx/float(xsplit))
    !             y = (ny/float(ysplit))
    !
    !             if (y > x) then
    !                 current = abs(1 - ( y / x ))
    !             else
    !                 current = abs(1 - ( x / y ))
    !             endif
    !
    !             if (current < best) then
    !                 best = current
    !                 xs = xsplit
    !                 ys = ysplit
    !             endif
    !         endif
    !     enddo
    !
    !     this%ximages = xs
    !     this%yimages = ys
    !
    !     this%ximg = mod(this_image(),  this%ximages)
    !     this%yimg =     this_image() / this%ximages
    !
    !     if (assertions) call assert((xs*ys) == nimages, "Number of tiles does not sum to number of images")
    !
    ! end subroutine domain_decomposition

end submodule


text:
module bst_mod
use posting_mod
implicit none

type bst_node
	type(bst_node), pointer :: left, right
	character(len=:), allocatable :: key
	type(posting) :: post
end type bst_node

type bst
	type(bst_node), pointer :: root
	integer :: size
contains
	procedure :: initialize
	procedure :: insert
	procedure :: print
	procedure :: write
end type bst

private :: initialize, insert, print, write

contains

subroutine initialize(tree)
	class(bst) :: tree
	tree%size = 0
	nullify(tree%root)
end subroutine initialize

subroutine insert(tree, key, docid)
	class(bst) :: tree
	character(len=*), intent(in) :: key
	integer :: docid

	type(bst_node), pointer :: p

	tree%size = tree%size + 1

	if (.NOT. associated(tree%root)) then
		call make_node(tree%root, key, docid)
	end if

	p => tree%root

	do
!print *, "a"
		if (key < p%key) then
!print *, "b"
			if (.NOT. associated(p%left)) then
				call make_node(p%left, key, docid)
				return
			end if
			p => p%left
		else if (key > p%key) then
!print *, "b"
			if (.NOT. associated(p%right)) then
				call make_node(p%right, key, docid)
				return
			end if
			p => p%right
		else
!print *, "b"
			call p%post%append(docid)
			tree%size = tree%size - 1
			return
		end if
	end do
end subroutine insert

subroutine make_node(node, key, docid)
	type(bst_node), pointer, intent(inout) :: node
	character(len=*), intent(in) :: key
	integer, intent(in) :: docid

	allocate(node)
	nullify(node%left)
	nullify(node%right)
	node%key = key

	call node%post%initialize()
	call node%post%append(docid)
end subroutine make_node

subroutine print(tree)
	class(bst) :: tree
	call print_inorder(tree%root)
end subroutine print

recursive subroutine print_inorder(node)
	type(bst_node), pointer :: node
	if (.NOT. associated(node)) then
		return
	end if
	call print_inorder(node%left)
	print *, node%key
	call node%post%print()
	call print_inorder(node%right)
end subroutine print_inorder

subroutine write(tree)
	class(bst) :: tree
	write(10) tree%size
	call write_key_inorder(tree%root)
	call write_post_inorder(tree%root)
end subroutine write

recursive subroutine write_key_inorder(node)
	type(bst_node), pointer :: node
	if (.NOT. associated(node)) then
		return
	end if
	call write_key_inorder(node%left)
	write(10) len(node%key)
	write(10) node%key
	call write_key_inorder(node%right)
end subroutine write_key_inorder

recursive subroutine write_post_inorder(node)
	type(bst_node), pointer :: node
	if (.NOT. associated(node)) then
		return
	end if
	call write_post_inorder(node%left)
	call node%post%write()
	call write_post_inorder(node%right)
end subroutine write_post_inorder

end module bst_mod

max_stars_repo_name: vkitchen/fortsearch

max_issues_repo_name: vkitchen/fortsearch

max_forks_repo_name: vkitchen/fortsearch


text:
program modules_18
use iso_fortran_env, only: sp=>real32
use modules_18b, only: f, g
implicit none
integer :: i, a
real(sp) :: r32

a = 3
r32 = 5
i = f(a, r32)
print *, i
if (i /= 8) error stop
call g(a, r32, i)
print *, i
if (i /= 8) error stop

end

max_stars_repo_name: Thirumalai-Shaktivel/lfortran

max_issues_repo_name: Thirumalai-Shaktivel/lfortran

max_forks_repo_name: Thirumalai-Shaktivel/lfortran


text:
!
!     CalculiX - A 3-dimensional finite element program
!              Copyright (C) 1998-2021 Guido Dhondt
!
!     This program is free software; you can redistribute it and/or
!     modify it under the terms of the GNU General Public License as
!     published by the Free Software Foundation(version 2);
!
!
!     This program is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!     GNU General Public License for more details.
!
!     You should have received a copy of the GNU General Public License
!     along with this program; if not, write to the Free Software
!     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
!
      subroutine map3dtolayer(yn,ipkon,kon,lakon,nfield,
     &  ne,co,ielmat,mi)
!
!     interpolates 3d field nodal values to nodal values in the
!     layers of composite materials
!
      implicit none
!
      character*8 lakon(*)
!
      integer ipkon(*),kon(*),ne,nfield,i,j,k,l,mi(*),nelem,
     &  ielmat(mi(3),*),nlayer,iflag,nbot20(8),nmid20(8),ntop20(8),
     &  nbot15(6),nmid15(6),ntop15(6),ibot,itop,nope,nopes,nopeexp,
     &  indexe,konl(20),imid,node,nodebot,nodetop
!
      real*8 yn(nfield,*),shp(4,20),xsj(3),co(3,*),xl(3,20),
     &  xi,et,ze,dd,dt,xi20(8),et20(8),xi15(6),et15(6)
!
!
!
      data iflag /1/
!
      data nbot20 /1,2,3,4,9,10,11,12/
      data nmid20 /17,18,19,20,0,0,0,0/
      data ntop20 /5,6,7,8,13,14,15,16/
!
      data nbot15 /1,2,3,7,8,9/
      data nmid15 /13,14,15,0,0,0/
      data ntop15 /4,5,6,10,11,12/
!
      data xi20 /-1.d0,1.d0,1.d0,-1.d0,0.d0,1.d0,0.d0,-1.d0/
      data et20 /-1.d0,-1.d0,1.d0,1.d0,-1.d0,0.d0,1.d0,0.d0/
!
      data xi15 /0.d0,1.d0,0.d0,0.5d0,0.5d0,0.d0/
      data et15 /0.d0,0.d0,1.d0,0.d0,0.5d0,0.5d0/
!
      include "gauss.f"
!
      do nelem=1,ne
         if(lakon(nelem)(7:8).eq.'LC') then
!
!        composite materials
!
!        determining the number of layers
!
            nlayer=0
            do k=1,mi(3)
               if(ielmat(k,nelem).ne.0) then
                  nlayer=nlayer+1
               endif
            enddo
!
            indexe=ipkon(nelem)
!
            if(lakon(nelem)(4:5).eq.'20') then
               nope=20
               nopes=8
               nopeexp=28
            elseif(lakon(nelem)(4:5).eq.'15') then
               nope=15
               nopes=6
               nopeexp=21
            endif
!
            do i=1,nope
               konl(i)=kon(indexe+i)
               do j=1,3
                  xl(j,i)=co(j,konl(i))
               enddo
            enddo
!
            do i=1,nopes
               if(lakon(nelem)(4:5).eq.'20') then
                  ibot=nbot20(i)
                  imid=nmid20(i)
                  itop=ntop20(i)
                  xi=xi20(i)
                  et=et20(i)
               else
                  ibot=nbot15(i)
                  imid=nmid15(i)
                  itop=ntop15(i)
                  xi=xi15(i)
                  et=et15(i)
               endif
!
               nodebot=konl(ibot)
               nodetop=konl(itop)
               dd=sqrt((co(1,nodebot)-co(1,nodetop))**2+
     &                 (co(2,nodebot)-co(2,nodetop))**2+
     &                 (co(3,nodebot)-co(3,nodetop))**2)
               do j=0,nlayer-1
!
!                 bottom node
!
                  node=kon(indexe+nopeexp+j*nope+ibot)
                  dt=sqrt((co(1,nodebot)-co(1,node))**2+
     &                 (co(2,nodebot)-co(2,node))**2+
     &                 (co(3,nodebot)-co(3,node))**2)
                  ze=2.d0*dt/dd-1.d0
c                  write(*,*) 'map3dtolayer',node,xi,et,ze
!
!                 determining the value of the shape functions
!
                  if(lakon(nelem)(4:5).eq.'20') then
                     call shape20h(xi,et,ze,xl,xsj,shp,iflag)
                  elseif(lakon(nelem)(4:5).eq.'15') then
                     call shape15w(xi,et,ze,xl,xsj,shp,iflag)
                  endif
!
                  do k=1,nfield
                     yn(k,node)=0.d0
                     do l=1,nope
                        yn(k,node)=yn(k,node)+
     &                        shp(4,l)*yn(k,konl(l))
c                        write(*,*) 'xi',xi,et,ze,node
c                        write(*,*) l,shp(4,l),yn(k,konl(l))
                     enddo
                  enddo
!
!                 top node
!
                  node=kon(indexe+nopeexp+j*nope+itop)
                  dt=sqrt((co(1,nodebot)-co(1,node))**2+
     &                 (co(2,nodebot)-co(2,node))**2+
     &                 (co(3,nodebot)-co(3,node))**2)
                  ze=2.d0*dt/dd-1.d0
c                  write(*,*) 'map3dtolayer',node,xi,et,ze
!
!                 determining the value of the shape functions
!
                  if(lakon(nelem)(4:5).eq.'20') then
                     call shape20h(xi,et,ze,xl,xsj,shp,iflag)
                  elseif(lakon(nelem)(4:5).eq.'15') then
                     call shape15w(xi,et,ze,xl,xsj,shp,iflag)
                  endif
!
                  do k=1,nfield
                     yn(k,node)=0.d0
                     do l=1,nope
                        yn(k,node)=yn(k,node)+
     &                        shp(4,l)*yn(k,konl(l))
                     enddo
                  enddo
!
!                 middle node, if any
!
                  if(i.gt.nopes/2) cycle
!
                  node=kon(indexe+nopeexp+j*nope+imid)
                  dt=sqrt((co(1,nodebot)-co(1,node))**2+
     &                 (co(2,nodebot)-co(2,node))**2+
     &                 (co(3,nodebot)-co(3,node))**2)
                  ze=2.d0*dt/dd-1.d0
c                  write(*,*) 'map3dtolayer',node,xi,et,ze
!
!                 determining the value of the shape functions
!
                  if(lakon(nelem)(4:5).eq.'20') then
                     call shape20h(xi,et,ze,xl,xsj,shp,iflag)
                  elseif(lakon(nelem)(4:5).eq.'15') then
                     call shape15w(xi,et,ze,xl,xsj,shp,iflag)
                  endif
!
                  do k=1,nfield
                     yn(k,node)=0.d0
                     do l=1,nope
                        yn(k,node)=yn(k,node)+
     &                        shp(4,l)*yn(k,konl(l))
                     enddo
                  enddo
               enddo
            enddo
         endif
      enddo
!
      return
      end

max_stars_repo_name: alleindrach/calculix-desktop

max_issues_repo_name: alleindrach/calculix-desktop

max_forks_repo_name: alleindrach/calculix-desktop


text:
! MD5 of template: a9d1b24b38764cb05db2f7aeda80e7ca
module grid_io

  use core_lib
  use mpi_hdf5_io
  use type_grid

  implicit none
  save

  private
  public :: grid_exists
  public :: read_grid_3d
  public :: read_grid_4d
  public :: write_grid_3d
  public :: write_grid_4d

  interface read_grid_3d
     module procedure read_grid_3d_sp
     module procedure read_grid_3d_dp
     module procedure read_grid_3d_int
     module procedure read_grid_3d_int8
  end interface read_grid_3d

  interface read_grid_4d
     module procedure read_grid_4d_sp
     module procedure read_grid_4d_dp
     module procedure read_grid_4d_int
     module procedure read_grid_4d_int8
  end interface read_grid_4d

  interface write_grid_3d
     module procedure write_grid_3d_sp
     module procedure write_grid_3d_dp
     module procedure write_grid_3d_int
     module procedure write_grid_3d_int8
  end interface write_grid_3d

  interface write_grid_4d
     module procedure write_grid_4d_sp
     module procedure write_grid_4d_dp
     module procedure write_grid_4d_int
     module procedure write_grid_4d_int8
  end interface write_grid_4d

contains

  logical function grid_exists(group, name)
    implicit none
    integer(hid_t),intent(in) :: group
    character(len=*),intent(in) :: name
    integer(hid_t) :: g_level, g_grid
    if(mp_path_exists(group, 'level_00001')) then
       g_level = mp_open_group(group, 'level_00001')
       if(mp_path_exists(g_level, 'grid_00001')) then
          g_grid = mp_open_group(g_level, 'grid_00001')
          if(mp_path_exists(g_grid, name)) then
             grid_exists = .true.
          else
             grid_exists = .false.
          end if
       else
          grid_exists = .false.
       end if
    else
       grid_exists = .false.
    end if
  end function grid_exists

  subroutine read_grid_4d_int8(group, path, array, geo)

    implicit none

    integer(hid_t), intent(in) :: group
    character(len=*), intent(in) :: path
    integer(idp), intent(out) :: array(:,:)
    type(grid_geometry_desc),intent(in),target :: geo
    integer(idp), allocatable :: array4d(:,:,:,:)
    character(len=100) :: full_path
    integer :: ilevel, igrid, idust
    type(level_desc), pointer :: level
    type(grid_desc), pointer :: grid

    do ilevel=1,size(geo%levels)
       level => geo%levels(ilevel)
       do igrid=1,size(level%grids)
          grid => level%grids(igrid)
          write(full_path, '("level_", I5.5, "/grid_", I5.5,"/")') ilevel, igrid
          full_path = trim(full_path)//trim(path)
          call mp_read_array_auto(group, full_path, array4d)
          if(any(is_nan(array4d))) call error("read_grid_4d", "NaN values in 4D array")
          do idust=1,size(array4d, 4)
             where(grid%goto_grid(1:grid%n1,1:grid%n2,1:grid%n3) > 0)
                array4d(:,:,:,idust) = 0
             end where
          end do
          array(grid%start_id:grid%start_id + grid%n_cells - 1, :) = reshape(array4d, (/grid%n_cells, size(array, 2)/))
       end do
    end do

    ! The following three lines provide a workaround for the PGI Fortran
    ! compiler, which otherwise crashes with the following error:
    ! 0: RESHAPE: result type != SOURCE type
  contains
    subroutine test()
    end subroutine test

  end subroutine read_grid_4d_int8

  subroutine read_grid_3d_int8(group, path, array, geo)

    implicit none

    integer(hid_t), intent(in) :: group
    character(len=*), intent(in) :: path
    integer(idp), intent(out) :: array(:)
    type(grid_geometry_desc),intent(in),target :: geo
    integer(idp), allocatable :: array3d(:,:,:)
    character(len=100) :: full_path
    integer :: ilevel, igrid
    type(level_desc), pointer :: level
    type(grid_desc), pointer :: grid

    do ilevel=1,size(geo%levels)
       level => geo%levels(ilevel)
       do igrid=1,size(level%grids)
          grid => level%grids(igrid)
          write(full_path, '("level_", I5.5, "/grid_ ", I5.5,"/")') ilevel, igrid
          full_path = trim(full_path)//trim(path)
          call mp_read_array_auto(group, full_path, array3d)
          if(any(is_nan(array3d))) call error("read_grid_3d", "NaN values in 3D array")
          where(grid%goto_grid(1:grid%n1,1:grid%n2,1:grid%n3) > 0)
             array3d(:,:,:) = 0
          end where
          array(grid%start_id:grid%start_id + grid%n_cells - 1) = reshape(array3d, (/grid%n_cells/))
       end do
    end do

    ! The following three lines provide a workaround for the PGI Fortran
    ! compiler, which otherwise crashes with the following error:
    ! 0: RESHAPE: result type != SOURCE type
  contains
    subroutine test()
    end subroutine test

  end subroutine read_grid_3d_int8

  subroutine write_grid_4d_int8(group, path, array, geo)

    implicit none

    integer(hid_t), intent(in) :: group
    character(len=*), intent(in) :: path
    integer(idp), intent(in) :: array(:,:)
    type(grid_geometry_desc),intent(in),target :: geo
    integer(hid_t) :: g_level, g_grid
    character(len=100) :: name
    integer :: ilevel, igrid
    type(level_desc), pointer :: level
    type(grid_desc), pointer :: grid

    do ilevel=1,size(geo%levels)
       level => geo%levels(ilevel)
       write(name, '("level_", I5.5)') ilevel
       if(mp_path_exists(group, name)) then
          g_level = mp_open_group(group, name)
       else
          g_level = mp_create_group(group, name)
       end if
       do igrid=1,size(level%grids)
          grid => level%grids(igrid)
          write(name, '("grid_", I5.5)') igrid
          if(mp_path_exists(g_level, name)) then
             g_grid = mp_open_group(g_level, name)
          else
             g_grid = mp_create_group(g_level, name)
          end if
          call mp_write_array(g_grid, path, reshape(array(grid%start_id:grid%start_id + grid%n_cells - 1, :), &
               &                                     (/grid%n1, grid%n2, grid%n3, size(array,2)/)))
          call mp_close_group(g_grid)
       end do
       call mp_close_group(g_level)

    end do

  end subroutine write_grid_4d_int8

  subroutine write_grid_3d_int8(group, path, array, geo)

    implicit none

    integer(hid_t), intent(in) :: group
    character(len=*), intent(in) :: path
    integer(idp), intent(in) :: array(:)
    type(grid_geometry_desc),intent(in),target :: geo
    integer(hid_t) :: g_level, g_grid
    character(len=100) :: name
    integer :: ilevel, igrid
    type(level_desc), pointer :: level
    type(grid_desc), pointer :: grid

    do ilevel=1,size(geo%levels)
       level => geo%levels(ilevel)
       write(name, '("level_", I5.5)') ilevel
       if(mp_path_exists(group, name)) then
          g_level = mp_open_group(group, name)
       else
          g_level = mp_create_group(group, name)
       end if
       do igrid=1,size(level%grids)
          grid => level%grids(igrid)
          write(name, '("grid_", I5.5)') igrid
          if(mp_path_exists(g_level, name)) then
             g_grid = mp_open_group(g_level, name)
          else
             g_grid = mp_create_group(g_level, name)
          end if
          call mp_write_array(g_grid, path, reshape(array(grid%start_id:grid%start_id + grid%n_cells - 1), (/grid%n1, grid%n2, grid%n3/)))
          call mp_close_group(g_grid)
       end do
       call mp_close_group(g_level)

    end do

  end subroutine write_grid_3d_int8

  subroutine read_grid_4d_int(group, path, array, geo)

    implicit none

    integer(hid_t), intent(in) :: group
    character(len=*), intent(in) :: path
    integer, intent(out) :: array(:,:)
    type(grid_geometry_desc),intent(in),target :: geo
    integer, allocatable :: array4d(:,:,:,:)
    character(len=100) :: full_path
    integer :: ilevel, igrid, idust
    type(level_desc), pointer :: level
    type(grid_desc), pointer :: grid

    do ilevel=1,size(geo%levels)
       level => geo%levels(ilevel)
       do igrid=1,size(level%grids)
          grid => level%grids(igrid)
          write(full_path, '("level_", I5.5, "/grid_", I5.5,"/")') ilevel, igrid
          full_path = trim(full_path)//trim(path)
          call mp_read_array_auto(group, full_path, array4d)
          if(any(is_nan(array4d))) call error("read_grid_4d", "NaN values in 4D array")
          do idust=1,size(array4d, 4)
             where(grid%goto_grid(1:grid%n1,1:grid%n2,1:grid%n3) > 0)
                array4d(:,:,:,idust) = 0
             end where
          end do
          array(grid%start_id:grid%start_id + grid%n_cells - 1, :) = reshape(array4d, (/grid%n_cells, size(array, 2)/))
       end do
    end do

    ! The following three lines provide a workaround for the PGI Fortran
    ! compiler, which otherwise crashes with the following error:
    ! 0: RESHAPE: result type != SOURCE type
  contains
    subroutine test()
    end subroutine test

  end subroutine read_grid_4d_int

  subroutine read_grid_3d_int(group, path, array, geo)

    implicit none

    integer(hid_t), intent(in) :: group
    character(len=*), intent(in) :: path
    integer, intent(out) :: array(:)
    type(grid_geometry_desc),intent(in),target :: geo
    integer, allocatable :: array3d(:,:,:)
    character(len=100) :: full_path
    integer :: ilevel, igrid
    type(level_desc), pointer :: level
    type(grid_desc), pointer :: grid

    do ilevel=1,size(geo%levels)
       level => geo%levels(ilevel)
       do igrid=1,size(level%grids)
          grid => level%grids(igrid)
          write(full_path, '("level_", I5.5, "/grid_ ", I5.5,"/")') ilevel, igrid
          full_path = trim(full_path)//trim(path)
          call mp_read_array_auto(group, full_path, array3d)
          if(any(is_nan(array3d))) call error("read_grid_3d", "NaN values in 3D array")
          where(grid%goto_grid(1:grid%n1,1:grid%n2,1:grid%n3) > 0)
             array3d(:,:,:) = 0
          end where
          array(grid%start_id:grid%start_id + grid%n_cells - 1) = reshape(array3d, (/grid%n_cells/))
       end do
    end do

    ! The following three lines provide a workaround for the PGI Fortran
    ! compiler, which otherwise crashes with the following error:
    ! 0: RESHAPE: result type != SOURCE type
  contains
    subroutine test()
    end subroutine test

  end subroutine read_grid_3d_int

  subroutine write_grid_4d_int(group, path, array, geo)

    implicit none

    integer(hid_t), intent(in) :: group
    character(len=*), intent(in) :: path
    integer, intent(in) :: array(:,:)
    type(grid_geometry_desc),intent(in),target :: geo
    integer(hid_t) :: g_level, g_grid
    character(len=100) :: name
    integer :: ilevel, igrid
    type(level_desc), pointer :: level
    type(grid_desc), pointer :: grid

    do ilevel=1,size(geo%levels)
       level => geo%levels(ilevel)
       write(name, '("level_", I5.5)') ilevel
       if(mp_path_exists(group, name)) then
          g_level = mp_open_group(group, name)
       else
          g_level = mp_create_group(group, name)
       end if
       do igrid=1,size(level%grids)
          grid => level%grids(igrid)
          write(name, '("grid_", I5.5)') igrid
          if(mp_path_exists(g_level, name)) then
             g_grid = mp_open_group(g_level, name)
          else
             g_grid = mp_create_group(g_level, name)
          end if
          call mp_write_array(g_grid, path, reshape(array(grid%start_id:grid%start_id + grid%n_cells - 1, :), &
               &                                     (/grid%n1, grid%n2, grid%n3, size(array,2)/)))
          call mp_close_group(g_grid)
       end do
       call mp_close_group(g_level)

    end do

  end subroutine write_grid_4d_int

  subroutine write_grid_3d_int(group, path, array, geo)

    implicit none

    integer(hid_t), intent(in) :: group
    character(len=*), intent(in) :: path
    integer, intent(in) :: array(:)
    type(grid_geometry_desc),intent(in),target :: geo
    integer(hid_t) :: g_level, g_grid
    character(len=100) :: name
    integer :: ilevel, igrid
    type(level_desc), pointer :: level
    type(grid_desc), pointer :: grid

    do ilevel=1,size(geo%levels)
       level => geo%levels(ilevel)
       write(name, '("level_", I5.5)') ilevel
       if(mp_path_exists(group, name)) then
          g_level = mp_open_group(group, name)
       else
          g_level = mp_create_group(group, name)
       end if
       do igrid=1,size(level%grids)
          grid => level%grids(igrid)
          write(name, '("grid_", I5.5)') igrid
          if(mp_path_exists(g_level, name)) then
             g_grid = mp_open_group(g_level, name)
          else
             g_grid = mp_create_group(g_level, name)
          end if
          call mp_write_array(g_grid, path, reshape(array(grid%start_id:grid%start_id + grid%n_cells - 1), (/grid%n1, grid%n2, grid%n3/)))
          call mp_close_group(g_grid)
       end do
       call mp_close_group(g_level)

    end do

  end subroutine write_grid_3d_int

  subroutine read_grid_4d_dp(group, path, array, geo)

    implicit none

    integer(hid_t), intent(in) :: group
    character(len=*), intent(in) :: path
    real(dp), intent(out) :: array(:,:)
    type(grid_geometry_desc),intent(in),target :: geo
    real(dp), allocatable :: array4d(:,:,:,:)
    character(len=100) :: full_path
    integer :: ilevel, igrid, idust
    type(level_desc), pointer :: level
    type(grid_desc), pointer :: grid

    do ilevel=1,size(geo%levels)
       level => geo%levels(ilevel)
       do igrid=1,size(level%grids)
          grid => level%grids(igrid)
          write(full_path, '("level_", I5.5, "/grid_", I5.5,"/")') ilevel, igrid
          full_path = trim(full_path)//trim(path)
          call mp_read_array_auto(group, full_path, array4d)
          if(any(is_nan(array4d))) call error("read_grid_4d", "NaN values in 4D array")
          do idust=1,size(array4d, 4)
             where(grid%goto_grid(1:grid%n1,1:grid%n2,1:grid%n3) > 0)
                array4d(:,:,:,idust) = 0
             end where
          end do
          array(grid%start_id:grid%start_id + grid%n_cells - 1, :) = reshape(array4d, (/grid%n_cells, size(array, 2)/))
       end do
    end do

    ! The following three lines provide a workaround for the PGI Fortran
    ! compiler, which otherwise crashes with the following error:
    ! 0: RESHAPE: result type != SOURCE type
  contains
    subroutine test()
    end subroutine test

  end subroutine read_grid_4d_dp

  subroutine read_grid_3d_dp(group, path, array, geo)

    implicit none

    integer(hid_t), intent(in) :: group
    character(len=*), intent(in) :: path
    real(dp), intent(out) :: array(:)
    type(grid_geometry_desc),intent(in),target :: geo
    real(dp), allocatable :: array3d(:,:,:)
    character(len=100) :: full_path
    integer :: ilevel, igrid
    type(level_desc), pointer :: level
    type(grid_desc), pointer :: grid

    do ilevel=1,size(geo%levels)
       level => geo%levels(ilevel)
       do igrid=1,size(level%grids)
          grid => level%grids(igrid)
          write(full_path, '("level_", I5.5, "/grid_ ", I5.5,"/")') ilevel, igrid
          full_path = trim(full_path)//trim(path)
          call mp_read_array_auto(group, full_path, array3d)
          if(any(is_nan(array3d))) call error("read_grid_3d", "NaN values in 3D array")
          where(grid%goto_grid(1:grid%n1,1:grid%n2,1:grid%n3) > 0)
             array3d(:,:,:) = 0
          end where
          array(grid%start_id:grid%start_id + grid%n_cells - 1) = reshape(array3d, (/grid%n_cells/))
       end do
    end do

    ! The following three lines provide a workaround for the PGI Fortran
    ! compiler, which otherwise crashes with the following error:
    ! 0: RESHAPE: result type != SOURCE type
  contains
    subroutine test()
    end subroutine test

  end subroutine read_grid_3d_dp

  subroutine write_grid_4d_dp(group, path, array, geo)

    implicit none

    integer(hid_t), intent(in) :: group
    character(len=*), intent(in) :: path
    real(dp), intent(in) :: array(:,:)
    type(grid_geometry_desc),intent(in),target :: geo
    integer(hid_t) :: g_level, g_grid
    character(len=100) :: name
    integer :: ilevel, igrid
    type(level_desc), pointer :: level
    type(grid_desc), pointer :: grid

    do ilevel=1,size(geo%levels)
       level => geo%levels(ilevel)
       write(name, '("level_", I5.5)') ilevel
       if(mp_path_exists(group, name)) then
          g_level = mp_open_group(group, name)
       else
          g_level = mp_create_group(group, name)
       end if
       do igrid=1,size(level%grids)
          grid => level%grids(igrid)
          write(name, '("grid_", I5.5)') igrid
          if(mp_path_exists(g_level, name)) then
             g_grid = mp_open_group(g_level, name)
          else
             g_grid = mp_create_group(g_level, name)
          end if
          call mp_write_array(g_grid, path, reshape(array(grid%start_id:grid%start_id + grid%n_cells - 1, :), &
               &                                     (/grid%n1, grid%n2, grid%n3, size(array,2)/)))
          call mp_close_group(g_grid)
       end do
       call mp_close_group(g_level)

    end do

  end subroutine write_grid_4d_dp

  subroutine write_grid_3d_dp(group, path, array, geo)

    implicit none

    integer(hid_t), intent(in) :: group
    character(len=*), intent(in) :: path
    real(dp), intent(in) :: array(:)
    type(grid_geometry_desc),intent(in),target :: geo
    integer(hid_t) :: g_level, g_grid
    character(len=100) :: name
    integer :: ilevel, igrid
    type(level_desc), pointer :: level
    type(grid_desc), pointer :: grid

    do ilevel=1,size(geo%levels)
       level => geo%levels(ilevel)
       write(name, '("level_", I5.5)') ilevel
       if(mp_path_exists(group, name)) then
          g_level = mp_open_group(group, name)
       else
          g_level = mp_create_group(group, name)
       end if
       do igrid=1,size(level%grids)
          grid => level%grids(igrid)
          write(name, '("grid_", I5.5)') igrid
          if(mp_path_exists(g_level, name)) then
             g_grid = mp_open_group(g_level, name)
          else
             g_grid = mp_create_group(g_level, name)
          end if
          call mp_write_array(g_grid, path, reshape(array(grid%start_id:grid%start_id + grid%n_cells - 1), (/grid%n1, grid%n2, grid%n3/)))
          call mp_close_group(g_grid)
       end do
       call mp_close_group(g_level)

    end do

  end subroutine write_grid_3d_dp

  subroutine read_grid_4d_sp(group, path, array, geo)

    implicit none

    integer(hid_t), intent(in) :: group
    character(len=*), intent(in) :: path
    real(sp), intent(out) :: array(:,:)
    type(grid_geometry_desc),intent(in),target :: geo
    real(sp), allocatable :: array4d(:,:,:,:)
    character(len=100) :: full_path
    integer :: ilevel, igrid, idust
    type(level_desc), pointer :: level
    type(grid_desc), pointer :: grid

    do ilevel=1,size(geo%levels)
       level => geo%levels(ilevel)
       do igrid=1,size(level%grids)
          grid => level%grids(igrid)
          write(full_path, '("level_", I5.5, "/grid_", I5.5,"/")') ilevel, igrid
          full_path = trim(full_path)//trim(path)
          call mp_read_array_auto(group, full_path, array4d)
          if(any(is_nan(array4d))) call error("read_grid_4d", "NaN values in 4D array")
          do idust=1,size(array4d, 4)
             where(grid%goto_grid(1:grid%n1,1:grid%n2,1:grid%n3) > 0)
                array4d(:,:,:,idust) = 0
             end where
          end do
          array(grid%start_id:grid%start_id + grid%n_cells - 1, :) = reshape(array4d, (/grid%n_cells, size(array, 2)/))
       end do
    end do

    ! The following three lines provide a workaround for the PGI Fortran
    ! compiler, which otherwise crashes with the following error:
    ! 0: RESHAPE: result type != SOURCE type
  contains
    subroutine test()
    end subroutine test

  end subroutine read_grid_4d_sp

  subroutine read_grid_3d_sp(group, path, array, geo)

    implicit none

    integer(hid_t), intent(in) :: group
    character(len=*), intent(in) :: path
    real(sp), intent(out) :: array(:)
    type(grid_geometry_desc),intent(in),target :: geo
    real(sp), allocatable :: array3d(:,:,:)
    character(len=100) :: full_path
    integer :: ilevel, igrid
    type(level_desc), pointer :: level
    type(grid_desc), pointer :: grid

    do ilevel=1,size(geo%levels)
       level => geo%levels(ilevel)
       do igrid=1,size(level%grids)
          grid => level%grids(igrid)
          write(full_path, '("level_", I5.5, "/grid_ ", I5.5,"/")') ilevel, igrid
          full_path = trim(full_path)//trim(path)
          call mp_read_array_auto(group, full_path, array3d)
          if(any(is_nan(array3d))) call error("read_grid_3d", "NaN values in 3D array")
          where(grid%goto_grid(1:grid%n1,1:grid%n2,1:grid%n3) > 0)
             array3d(:,:,:) = 0
          end where
          array(grid%start_id:grid%start_id + grid%n_cells - 1) = reshape(array3d, (/grid%n_cells/))
       end do
    end do

    ! The following three lines provide a workaround for the PGI Fortran
    ! compiler, which otherwise crashes with the following error:
    ! 0: RESHAPE: result type != SOURCE type
  contains
    subroutine test()
    end subroutine test

  end subroutine read_grid_3d_sp

  subroutine write_grid_4d_sp(group, path, array, geo)

    implicit none

    integer(hid_t), intent(in) :: group
    character(len=*), intent(in) :: path
    real(sp), intent(in) :: array(:,:)
    type(grid_geometry_desc),intent(in),target :: geo
    integer(hid_t) :: g_level, g_grid
    character(len=100) :: name
    integer :: ilevel, igrid
    type(level_desc), pointer :: level
    type(grid_desc), pointer :: grid

    do ilevel=1,size(geo%levels)
       level => geo%levels(ilevel)
       write(name, '("level_", I5.5)') ilevel
       if(mp_path_exists(group, name)) then
          g_level = mp_open_group(group, name)
       else
          g_level = mp_create_group(group, name)
       end if
       do igrid=1,size(level%grids)
          grid => level%grids(igrid)
          write(name, '("grid_", I5.5)') igrid
          if(mp_path_exists(g_level, name)) then
             g_grid = mp_open_group(g_level, name)
          else
             g_grid = mp_create_group(g_level, name)
          end if
          call mp_write_array(g_grid, path, reshape(array(grid%start_id:grid%start_id + grid%n_cells - 1, :), &
               &                                     (/grid%n1, grid%n2, grid%n3, size(array,2)/)))
          call mp_close_group(g_grid)
       end do
       call mp_close_group(g_level)

    end do

  end subroutine write_grid_4d_sp

  subroutine write_grid_3d_sp(group, path, array, geo)

    implicit none

    integer(hid_t), intent(in) :: group
    character(len=*), intent(in) :: path
    real(sp), intent(in) :: array(:)
    type(grid_geometry_desc),intent(in),target :: geo
    integer(hid_t) :: g_level, g_grid
    character(len=100) :: name
    integer :: ilevel, igrid
    type(level_desc), pointer :: level
    type(grid_desc), pointer :: grid

    do ilevel=1,size(geo%levels)
       level => geo%levels(ilevel)
       write(name, '("level_", I5.5)') ilevel
       if(mp_path_exists(group, name)) then
          g_level = mp_open_group(group, name)
       else
          g_level = mp_create_group(group, name)
       end if
       do igrid=1,size(level%grids)
          grid => level%grids(igrid)
          write(name, '("grid_", I5.5)') igrid
          if(mp_path_exists(g_level, name)) then
             g_grid = mp_open_group(g_level, name)
          else
             g_grid = mp_create_group(g_level, name)
          end if
          call mp_write_array(g_grid, path, reshape(array(grid%start_id:grid%start_id + grid%n_cells - 1), (/grid%n1, grid%n2, grid%n3/)))
          call mp_close_group(g_grid)
       end do
       call mp_close_group(g_level)

    end do

  end subroutine write_grid_3d_sp

end module grid_io

max_stars_repo_name: bluescarni/hyperion

max_issues_repo_name: bluescarni/hyperion

max_forks_repo_name: bluescarni/hyperion


text:
program array5
implicit none

real, dimension(3) :: a, b, c
real :: e(0)

a = [1, 2, 3]
b = (/ 1, 2, 3 /)
c = [real:: 1, 2, 3]
e = [real ::]

end program


text:
!***********************************************************************
!                                                                      *
      SUBROUTINE ALCBUF(MODE)
!                                                                      *
!   The arrays in  MODULE BUFFER_C are allocated (MODE = 1), realloca-   *
!   ted (MODE = 2), and deallocated (MODE = 3) in this routine.        *
!                                                                      *
!   Call(s) to: [LIB92]: ALLOC, DALLOC, RALLOC.                        *
!                                                                      *
!   Written by Farid A. Parpia.             Last update: 06 Oct 1992   *
!   Editted after translation by C. Froese Fischer       29 Apr 2007   *
!                                                                      *
!***********************************************************************
!...Translated by Pacific-Sierra Research 77to90  4.3E  16:02:30   2/14/04
!...Modified by Charlotte Froese Fischer
!                     Gediminas Gaigalas  10/05/17
!====================================================================

!-----------------------------------------------
!   I n t e r f a c e   B l o c k s
!-----------------------------------------------
      USE memory_man
      use buffer_C,       ONLY: LABEL, COEFF, NBDIM
      IMPLICIT NONE
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      INTEGER, INTENT(IN) :: MODE
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      INTEGER :: NEWSIZ, IERR
!
!
      SELECT CASE (MODE)
      CASE (1)
!
!   Allocation
!
         NBDIM = 10
         CALL ALLOC (LABEL, 6,NBDIM, 'LABEL', 'ALCBUF')
         CALL ALLOC (COEFF, NBDIM, 'COEFF', 'ALCBUF')
!
      CASE (2)
!
!   Reallocation to double storage
!
         NEWSIZ = NBDIM + NBDIM
         CALL RALLOC(LABEL,6,NEWSIZ, 'LABEL', 'ALCBUF')
         CALL RALLOC(COEFF,NEWSIZ, 'COEFF', 'ALCBUF')
         NBDIM = NEWSIZ
!
      CASE (3)
!
!   Deallocation
!
          CALL DALLOC(label, 'LABEL', 'ALCBUF')
          CALL DALLOC(coeff, 'COEFF', 'ALCBUF')
!
      CASE DEFAULT
!
!   Argument error
!
         WRITE (6, *) 'ALCBUF: Invalid argument MODE = ', MODE
         STOP
!
      END SELECT
!
      RETURN
!
      END SUBROUTINE ALCBUF

max_stars_repo_name: sylas/grasp-continuum

max_issues_repo_name: sylas/grasp-continuum

max_forks_repo_name: sylas/grasp-continuum


program band_util

  implicit none

  integer,parameter :: unit1=5, unit2=1
  integer :: nband, nkpt, nseg
  integer :: iseg,iband,ikpt,idat,jdat
  real(8),allocatable :: esp(:,:),kpt(:,:)
  integer :: idummy(2)
  real(8) :: rdummy(3),d,p

  write(*,*) "[input the # of bands, k points/segment, and segments]"
  write(*,*) "nband, nkpt, nsegment = ? ? ?"
  read(*,*) nband, nkpt, nseg

  allocate( kpt(3,nkpt*nseg) ); kpt=0.0d0
  allocate( esp(nband,nkpt*nseg) ); esp=0.0d0

  open(unit1,file='band.dat',status='old')
  idat=0
  jdat=0
  do iseg=1,nseg
     do ikpt=1,nkpt
        idat=idat+1
        read(unit1,*) idummy(1),rdummy(1:3),kpt(1:3,idat)
     end do
     do ikpt=1,nkpt
        jdat=jdat+1
        do iband=1,nband
           read(unit1,*) idummy(1:2),esp(iband,jdat)
        end do
     end do
  end do
  close(unit1)

  p=0.0d0
  do idat=1,nkpt*nseg
     if( idat > 1 )then
        d=sqrt( sum( (kpt(1:3,idat)-kpt(1:3,idat-1))**2 ) )
        p=p+d
     end if
     do iband=1,nband
        write(1000+iband,*) p, esp(iband,idat)
     end do
  end do

end program band_util


text:
      program wtype9
      call foo(x,i)
      end
      subroutine foo(x,y)
      i=1
      end


! { dg-do run }
! Test mismatched type kinds in a select statement.
program select_5
  integer(kind=1) i          ! kind = 1, -128 <= i < 127
  do i = 1, 3
    select case (i)

    ! kind = 4, reachable
    case (1_4)
      if (i /=  1_4) STOP 1

    ! kind = 8, reachable
    case (2_8)
      if (i /= 2_8) STOP 2

    ! kind = 4, unreachable because of range of i
    case (200)                       ! { dg-warning "not in the range" }
      STOP 3

    case default
      if (i /= 3) STOP 4
    end select
  end do
end program select_5


text:
   subroutine set_bc_mask(lo, hi, domlo, domhi, &
                          bcMask, bcMask_l1, bcMask_l2, bcMask_h1, bcMask_h2) &
      bind(C, name="set_bc_mask")

      use prob_params_module, only : physbc_lo, physbc_hi, problo, probhi, &
                                   Interior, Inflow, Outflow, Symmetry, SlipWall, NoSlipWall

      implicit none

     integer, intent(in   ) :: lo(2), hi(2), domlo(2), domhi(2)
     integer, intent(in   ) :: bcMask_l1, bcMask_l2, bcMask_h1, bcMask_h2
     integer, intent(inout) :: bcMask(bcMask_l1:bcMask_h1,bcMask_l2:bcMask_h2,2)

     ! comp=1 holds physical bc type (Inflow, Outflow, NoSlipWall
     ! comp=2 hold flag for whether is adiabatic (=0) or isothermal (=1)   TODO: check this

     if (bcMask_l1 < domlo(1)) then
        bcMask(domlo(1),bcMask_l2:bcMask_h2  ,1) = physbc_lo(1)
     end if

     if (bcMask_h1 > domhi(1)) then
        bcMask(domhi(1)+1,bcMask_l2:bcMask_h2,1) = physbc_hi(1)
     end if

     if (bcMask_l2 < domlo(2)) then
        bcMask(bcMask_l1:bcMask_h1,domlo(2)  ,1) = physbc_lo(2)
     end if

     if (bcMask_h2 > domhi(2)) then
        bcMask(bcMask_l1:bcMask_h1,domhi(2)+1,1) = physbc_hi(2)
     end if

   end subroutine set_bc_mask

max_stars_repo_name: wjge/PeleC

max_issues_repo_name: marient/PeleC

max_forks_repo_name: marient/PeleC


text:
c El fichero de abundancias debe contener una cabecera
c el numero de lineas de la cabecera no debe exceder de 20 y
c cada linea no debe tener mas de 70 caracters ( si tiene
c mas no se escribiran al modificar e fichero)
c La cabecera debe acabar necesariamente con una linea que utilice alguno
c de los simbolos (ver el 'if' de lectura en el 'do while')

	subroutine leeabun(iesc,abu)           !iesc=1 escribo

        parameter (na=92)
	character mensaje*26
        integer in
        real aa,abu(na)
        character*100 fichabun
        common/fichabun/fichabun
	common/calerr/icalerr     !si calculo errores=1 else =0

	mensaje=' containing the abundances'
	ican=1
	call cabecera(ican,fichabun,mensaje,ifail)
	if(ifail.eq.1)goto 999

         do i=1,na
           read(ican,*,err=100)in,aa
           abu(in)=aa          !10.**(aa-12.)
         end do
100	 close(1)
         return

999	print*,' '
	print*,'STOP: The file containing the abundances does NOT exist:'
	print*,fichabun
	print*,' '
	print*,'________________________________________________________'
	stop

        end


! ==============================================================================
! Subroutine: HELX_REAC (R,I)
!
! Purpose: Recomputes the reaction rates for opening one of the internal
!          base-pairs within a continuous section of helix.
!
! Method:
!
!               A-U           A-U
!               x-x   --->   x   x
!               U-A           U-A
!
! Arguments:
!
!             R - Class structure containing information on the
!                 RNA secondary structure and possible reactions.
!             I - The indx number of one of the nucelotides which
!                 belongs to the helix.
!
! History:
!
! Version    Date         Comment
! --------   ----------   -----------------------
!            01/01/2015   Original Code
!
! Dependancies:
!
! Modules -
! Functions -
! Subroutines -
!
! Author(s): Eric Dykeman
!            Copyright (c) 2015 (Please Refer to LICENCE)
!
! ==============================================================================

      SUBROUTINE HELX_REAC (R,I)

        IMPLICIT NONE

        !=== ARGUMENTS ===!

        TYPE(RNA_STRUC), INTENT(INOUT) :: r

        INTEGER, INTENT(IN) :: i

        !=== VARIABLES ===!

        INTEGER :: ke,ip,jp,indx

        DOUBLE PRECISION :: x,dg,atot

        jp = i

        ip = r% ibsp(jp)
        jp = MIN(ip,jp)

        DO WHILE ( r% link(jp) == 0 )

          jp = jp + 1

        ENDDO

        indx = r% link(jp)

        !=== Open BP Inside Helix ===!

        ke = r% ibsp(jp)
        ip = r% ibsp(jp)

        atot = r% wrk1(jp)

        IF ( r% link(ke) == 0 ) THEN

          ip = ip + 1
          jp = jp - 1

          DO WHILE ( r% link(ip) == 0 )

            CALL DELTAG_HI (r,ip,jp,dg)

            dg = dg / 2.0d0

            x = beta * dg
            x = DEXP(-x) * rateh

            r% wrk1(ip) = x
            r% wrk1(jp) = 0.0d0

            r% wrk2(ip) = 0.0d0
            r% wrk2(jp) = 0.0d0

            atot = atot + x

            ip = ip + 1
            jp = jp - 1

          ENDDO

        ENDIF

        r% ptot(indx) = atot

        CALL LOOP_RESUM (r,indx)

        RETURN

      END SUBROUTINE HELX_REAC


! { dg-do run }
!
! Check handling of errmsg.
!
implicit none
integer, allocatable :: a[:], b(:)[:], c, d(:)
integer :: stat
character(len=300) :: str

allocate(a[*], b(1)[*], c, d(2), stat=stat)

str = repeat('X', len(str))
allocate(a[*], stat=stat, errmsg=str)
!print *, stat, trim(str)
if (stat == 0 .or. str /= "Attempt to allocate an allocated object") &
  call abort ()

str = repeat('Y', len(str))
allocate(b(2)[*], stat=stat, errmsg=str)
!print *, stat, trim(str)
if (stat == 0 .or. str /= "Attempt to allocate an allocated object") &
  call abort ()

str = repeat('Q', len(str))
allocate(c, stat=stat, errmsg=str)
!print *, stat, trim(str)
if (stat == 0 .or. str /= "Attempt to allocate an allocated object") &
  call abort ()

str = repeat('P', len(str))
allocate(d(3), stat=stat, errmsg=str)
!print *, stat, trim(str)
if (stat == 0 .or. str /= "Attempt to allocate an allocated object") &
  call abort ()

end