module plant

use parameters_site
use parameters_plant
use environment
use brent

implicit none

! Plant variables
integer :: NOHARV ! simon removed NOHARV switch, so not used
real :: CRESMX,DAYLGE,FRACTV,GLVSI,GSTSI,LERG,LERV,LUEMXQ,NELLVG,PHENRF,PHOT,RESMOB
real :: RDLVD, ALLOTOT,GRESSI,GSHSI,GLAISI,SOURCE,SINK1T,CSTAV,TGE
real :: RDRFROST,RDRT,RDRL,RDRTOX,RESPGRT,RESPGSH,RESPHARD,RESPHARDSI,RESNOR,RLEAF,RplantAer,SLANEW
real :: RATEH,reHardPeriod,RDRTIL,RDRS,RDRW ! Simon renamed TV2TIL to RDRTIL
real :: CRESMN,DAYLGEMX
real :: ALLOSH, ALLORT, ALLOLV, ALLOST, FS, ALLOFRAC

! define variables for the the harvest fraction optimisation, needs to be global

real :: clv_cres_ect, fhageer, HAGRE_stuff, goal

contains

  real function f(x)
    IMPLICIT NONE
    real, intent(in):: x
  f = clv_cres_ect * x * 10 + x ** (1 - fhageer) * HAGRE_stuff * 10 - goal
  end function f

! Calculate Harvest GSTUB,HARVLA,HARVLV,HARVPH,HARVRE,HARVST,HARVTILG2,HARVFR
! Simon plant processes are now calculated as if harvest did not happen
Subroutine Harvest(day, NDAYS, NHARVCOL, BASAL, CLV,CRES,CST,CSTUB,CLVD,DAYS_HARVEST,LAI,PHEN,TILG2,TILG1,TILV, &
                             GSTUB,HARVLA,HARVLV,HARVLVD,HARVPH,HARVRE,HARVST, &
                             HARVTILG2,HARVFR,HARVFRIN,HARV,RDRHARV, WEED_HARV_FR, &
                    DM_RYE_RM, DM_WEED_RM, DMH_RYE, DMH_WEED)
  integer :: day
  integer :: NDAYS, NHARVCOL
  real, dimension(NDAYS, NHARVCOL) :: DAYS_HARVEST     ! major re-structure by Matt Hanson
  real    :: BASAL, CLV, CRES, CST, CSTUB, CLVD, LAI, PHEN, TILG2, TILG1, TILV
  real    :: GSTUB, HARVLV, HARVLVD, HARVLA, HARVRE, HARVTILG2, HARVST, HARVPH
  real    :: CLAI, HARVFR, TV1, HARVFRIN, RDRHARV, HARVFRST, DIESFRST, DMH_RYE, DMH_WEED
  real    :: WEED_HARV_FR ! fraction of harvest yeild from weed species, outputs to calc weed yeild
  integer :: HARV
!  integer :: i

  real ::  FRAC_HARV
  real ::  HARV_TRIG
  real ::  HARV_TARG
  real ::  WEED_DM_FRAC, DM_RM, DM_RYE_RM, DM_WEED_RM
  logical :: temp_opt_harvfrin

  ! set parameters from days_harvest
  FRAC_HARV = DAYS_HARVEST(day, 3)
  HARV_TRIG = DAYS_HARVEST(day, 4)
  HARV_TARG = DAYS_HARVEST(day, 5)
  WEED_DM_FRAC = DAYS_HARVEST(day, 6)
  temp_opt_harvfrin = opt_harvfrin

  ! calculate dry matter of ryegrass + weeds, include the harvestable fraction of dry matter
  DMH_RYE        = ((CLV+CST+CSTUB)/0.45 + CRES/0.40 + (CLVD * HARVFRD / 0.45)) * 10.0
  DMH_WEED =  WEED_DM_FRAC*DMH_RYE/BASAL*(1-BASAL)

  ! if above trigger and trigger >=0 (HARV_TRIG<0, flag for no harvest) then harvest
  if (((DMH_RYE + DMH_WEED) >= HARV_TRIG) .and. (HARV_TRIG>=0)) then
      HARV = 1
      if (FIXED_REMOVAL) then
          ! harvest assuming that the target is a fixed volume to harvest and
          ! harvest weeds and rye propotionally according to dry matter content
          DM_RM = HARV_TARG * FRAC_HARV ! amount of total dry matter to remove

          DM_RYE_RM = DM_RM * (DMH_RYE/(DMH_RYE+DMH_WEED)) ! amount of rye to harvest
          DM_WEED_RM = DM_RM * (DMH_WEED/(DMH_RYE+DMH_WEED)) ! amount of weed to harvest

          HARVFRIN = DM_RYE_RM/DMH_RYE
          WEED_HARV_FR = DM_WEED_RM/DM_RYE_RM

          ! set values for the optimisation function
          goal = DM_RYE_RM
          clv_cres_ect = (CLV / 0.45 + CLVD * HARVFRD / 0.45) + (CRES * CLV / (CLV + CST + CSTUB) / 0.40)
          fhageer= HAGERE
          HAGRE_stuff = (CST / 0.45 + CRES * CST / (CLV + CST + CSTUB) / 0.40)

      else
          ! harvest assuming that the goal is to harvest to the target dry matter
          ! harvest weeds and rye propotionally according to dry matter content
          DM_RM = ((DMH_RYE + DMH_WEED) - HARV_TARG) * FRAC_HARV ! amount of total dry matter to remove

          DM_RYE_RM = DM_RM * (DMH_RYE/(DMH_RYE+DMH_WEED)) ! amount of rye to harvest
          DM_WEED_RM = DM_RM * (DMH_WEED/(DMH_RYE+DMH_WEED)) ! amount of weed to harvest

          HARVFRIN = DM_RYE_RM/DMH_RYE
          WEED_HARV_FR = DM_WEED_RM/DM_RYE_RM

          ! set values for the optimisation function
          goal = DM_RYE_RM
          clv_cres_ect = (CLV / 0.45 + CLVD * HARVFRD / 0.45) + (CRES * CLV / (CLV + CST + CSTUB) / 0.40)
          fhageer= HAGERE
          HAGRE_stuff = (CST / 0.45 + CRES * CST / (CLV + CST + CSTUB) / 0.40)
      end if

      if (HARVFRIN<=0) then
        temp_opt_harvfrin = .false. ! to ensure that we cannot have a optimisation where there is zero harvest,
        ! this could break harvfrin
      end if
      if (temp_opt_harvfrin) then
     ! estimate the fraction of harvest to undertake using brent zero
     HARVFRIN = zero(0.0,1.0,& ! bounds
       1e-5, & ! machine tolerance
      1e-5, & ! tolerance
      f) ! function to minimize
      end if

  else
      ! do not harvest
      HARV = 0
      HARVFRIN = 0.0 ! no harvest
      WEED_HARV_FR = 0.0
      DM_RYE_RM = 0 ! amount of rye to harvest
      DM_WEED_RM = 0 ! amount of weed to harvest
      ! difference between HARVFRIN and HARVFR is unclear, but I belive that HARVFR is simply a legacy variable
  end if

  HARVFR = HARVFRIN  ! Simon just use entered fraction for leaf. Otherwise dead can have too much effect.

  FRACTV = (TILV + TILG1)/(TILG2 + TILG1 + TILV) ! Fraction of non-elongating tillers (Simon included TILG1)
  HARVFRST  = HARVFR ** (1-HAGERE)                                             ! Simon proportion of CST harvested
  DIESFRST  = 1.0 - HARVFRST                                                   ! Simon proportion of CST that dies
  TV1       = (HARVFR * CLV + HARVFRST * CST + 0 * CSTUB)/(CLV + CST + CSTUB)  ! Simon proportion of CRES harvested
  HARVFR    = HARVFR * HARV                                                    ! Simon only return HARVFR on HARV days
  RDRHARV   = RDRHARVMAX * HARVFR                                              ! Simon relative death rate due to harvest

! HARVFR = Fraction of leaf                                    that is harvested
! 1.0    = Fraction of leaf              in elongating tillers that is harvested (we assume)
! HAGERE = Fraction of stem and reserves in elongating tillers that is harvested (parameter)
! TV1    = Fraction of reserves                                that is harvested

  HARVLA    = (HARV   * LAI * HARVFR) / DELT
  HARVLV    = (HARV   * CLV * HARVFR) / DELT
  HARVLVD   = (HARV   * CLVD * HARVFR * HARVFRD) / DELT
  HARVPH    = (HARV   * PHEN        ) / DELT           ! PHEN zeroed after each harvest
  HARVST    = (HARV   * CST * HARVFRST) / DELT         ! Simon separated out GSTUB from HARVST
  GSTUB     = (HARV   * CST * DIESFRST) / DELT         ! Simon allowed stem survival when HARVFRST + DIESFRST < 1
  HARVRE    = (HARV   * CRES * TV1  ) / DELT
  HARVTILG2 = (HARV   * TILG2       ) / DELT           ! TILG2 zeroed after each harvest
end Subroutine Harvest

! Calculate RESNOR (relative amount of CRES)
Subroutine Biomass(AGE,CLV,CRES,CST,CSTUB)
  real :: AGE, CLV, CRES, CST, CSTUB
!  CRESMX = COCRESMX * (CLV + CRES + CST)     ! Maximum reserves in aboveground biomass (not stubble) in terms of C (not DM)
  CRESMX = COCRESMX * (CLV + CST)            ! Maximum reserves in aboveground biomass (not stubble) in terms of C (not DM)
  CRESMN = FCOCRESMN * CRESMX                ! Minimum reserves in aboveground biomass (not stubble) in terms of C (not DM)
!  RESNOR = max(0.0, min(1.0, (CRES-CRESMN)/(CRESMX-CRESMN) )) ! Simon revised normalisation of CRES relative to upper and lower "bounds" (seems to break the balance)
  RESNOR = max(0.0, min(1.0, CRES/CRESMX )) ! CRES normalised as a proportion of maximum
end Subroutine Biomass

! Calculate phenological changes
Subroutine Phenology(DAYL,TILG2,PHEN, DPHEN,GPHEN,HARVPH)
  real :: DAYL,TILG2,PHEN
  real :: DPHEN,GPHEN,HARVPH
  if (TILG2 > 0.0) then                                                       ! Simon PHEN only refers to elongating tillers
    GPHEN = max(0., (DAVTMP-0.01)*0.000144*24. * (min(DAYLP,DAYL)-0.24) ) ! Basically degree days * day length
    DPHEN = 0.
  else
	GPHEN = 0.
    DPHEN  = PHEN / DELT
  end if
!  if (DAYL < DAYLB) then                                       ! Simon adjusted resetting of PHEN whenever DAYL < DAYLB
!  if (DAYL < DAYLRV) then                                       ! Simon adjusted resetting of PHEN whenever DAYL < DAYLRV
!    GPHEN  = 0.0
!    DPHEN  = PHEN / DELT
!  end if
  PHENRF = max(0.0, min(1.0, (1 - PHEN)/(1 - PHENCR) ))        ! Phenological stage decreases leaf number and appearance on elongating tillers
!  DAYLGE = max(0.0, min(1.0, (DAYL - DAYLB)/(DLMXGE - DAYLB) ))! Day length increases tillering, leaf appearance, leaf elongation (very crude)
  if (DLMXGE /= DAYLB) then
    DAYLGE = DAYLGEMN + (1-DAYLGEMN) * max(0.0, min(1.0, (DAYL - DAYLB)/(DLMXGE - DAYLB) ))! Simon added DAYLGEMN following STICS model
  else if (DAYL >= DAYLB) then
    DAYLGE = 1.0
  else
    DAYLGE = DAYLGEMN
  end if
end Subroutine Phenology

! Simon added vernalisation function, based on STICS model (Brisson et al 2009)
! Calculate vernalisation VERN, which allows RGRTVG1 = relative growth rate of generative tillers
Subroutine Vernalisation(DAYL,PHEN,YDAYL,TMMN,TMMX,DAVTMP,Tsurf,VERN,VERND, DVERND)
  real :: DAYL, PHEN, YDAYL, TMMN, TMMX, DAVTMP,Tsurf
!  integer :: VERN
  real :: VERN
  real :: VERND, DVERND
!  real :: X, Y
!  ! assume no change in vernalisation
!  DVERND = 0.
!  ! accumulate sum of cold temperatures
!  if (VERN==0) then
!    if (TVERN.le.TMMN) then
!      DVERND = 0.0
!    else if (TVERN.ge.TMMX) then
!      DVERND = 1.0
!	else
!      DVERND = 1.0
!!	  Y      = (TVERN-TMMN)/(TMMX-TMMN) * 2.0 - 1.0   ! TVERN relative to max and min
!!	  X      = acos(Y)                                ! position on first half of cos wave
!!      DVERND = 1.0 - X / pi                           ! proportion of day below TVERN
!	end if
!  end if
!  ! does vernalisation occur?
!  if ((VERN==0).and.(VERND .ge. TVERND)) then
!	VERN = 1
! 	VERND = 0.0
!    DVERND = 0.0
!  end if
  if ((DAYL<YDAYL).and.(DAYL<=DAYLRV).and.(DAYLRV<=YDAYL)) then ! Reset vernalisation when daylength shortens after Solstice
	VERN = 0.0
	VERND  = 0.0
    DVERND = 0.0
  end if
  if (DAYL<=DAYLRV) then ! Vernalisation rate based on STICS and Streck models
    DVERND  = max(0.0, 1.0 - ((Tsurf - TVERN) / 7.5)**2)
  else
    DVERND  = 0.0
  end if
end Subroutine Vernalisation

! Simon renamed Foliage1() to CalcSLA()
! Calculate leaf elongation rates LERV, LERG and SLANEW of new leaves
Subroutine CalcSLA
  real :: EFFTMP, SLAMIN
  EFFTMP = max(TBASE, DAVTMP)
  ! Linear relationship based on Peacock 1976 (who did not include daylength effect)
  ! See also Hoglind et al 2001 - different eqn for LERG
  ! See also Hogling et al 2016 - DAYLGE applied to LERG (paper has typo)
!  LERV   =          max(0., (-0.76 + 0.52*EFFTMP)/1000. ) ! m d-1 leaf elongation rate on vegetative tillers (for timothy, Peacock 1976)
!  LERG   = DAYLGE * max(0., (-5.46 + 2.80*EFFTMP)/1000. ) ! Why is DAYLGE applied here and not to LERV? Bug when DAYLGE is always small?
!  LERV   =          max(0., (-1.13 + 0.75*EFFTMP)/1000. ) ! m d-1 leaf elongation rate on vegetative tillers (Simon, for ryegrass, Peacock 1976)
!  LERG   = DAYLGE * max(0., (-8.21 + 1.75*EFFTMP)/1000. ) ! m d-1 leaf elongation rate on generative tillers (Simon, for ryegrass, Peacock 1976)
  LERV   =          max(0., (LERVA + LERVB*EFFTMP)/1000. ) ! m d-1 leaf elongation rate on vegetative tillers (Hjelkrem et al EM 2017)
  LERG   =  max(0., (LERGA + LERGB*EFFTMP)/1000. ) ! m d-1 leaf elongation rate on generative tillers (Hjelkrem et al EM 2017)
  SLAMIN = SLAMAX * FSLAMIN
  SLANEW = SLAMAX - RESNOR * ( SLAMAX - SLAMIN )          ! m2 leaf gC-1 SLA of new leaves (depends on CRES) note unusual units!
end Subroutine CalcSLA

! Calculate light use efficiency LUEMXQ
Subroutine LUECO2TM(PARAV,BASAL) ! also uses KLUETILG, FRACTV, KLAI
!=============================================================================
! Calculate LUEMXQ (mol CO2 mol-1 PAR quanta)
! Inputs : PARAV (micromol PAR quanta m-2 s-1)
! See equations in M. van Oijen et al. / Ecological Modelling 179 (2004) 39-60 (for spring wheat)
! See aldo Rodriguez et al 1999
!=============================================================================
  real :: PARAV,BASAL
  real :: CO2I, EA, EAKMC, EAKMO, EAVCMX, EFF, GAMMAX, KC25, KMC, KMC25
  real :: KMO, KMO25, KOKC, O2, PMAX, R, RUBISCN, T, TMPFAC, VCMAX
  T      = DAVTMP                                            ! degC
  RUBISCN = RUBISC * (1.E6/550000.)                          ! mumol m-2 leaf Rubisco content of upper leaves
  EAVCMX =  68000                                            ! J mol-1 Activation energy for VCMAX
  EAKMC  =  65800                                            ! J mol-1 Activation energy for KMC
  EAKMO  =   1400                                            ! J mol-1 Activation energy for KMO
  KC25   =     20                                            ! mol CO2 mol-1 Rubisco s-1 Catalytic efficiency of Rubisco at 25 degC
  KMC25  =    460                                            ! ppm CO2 Km-value Rubisco for carboxylation at 25 degC
  KMO25  =     33                                            ! % O2	Km-value Rubisco for oxygenation at 25 degC
  KOKC   =      0.21                                         ! Catalytic efficiency ratio Rubisco oxygenation/carboxylation
  O2     =     21                                            ! % O2	Oxygen concentration in chloroplasts
  R      =      8.314                                        ! J K-1 mol-1 Universal gas constant
  CO2I   = 0.7 * CO2A                                        ! ppm CO2 concentration in chloroplasts (Eqn 8)
  VCMAX  = RUBISCN * KC25 * exp((1/298.-1/(T+273))*EAVCMX/R) ! micromol CO2 m-2 leaf s-1 Maximum carboxylation rate in upper leaves (Eqn 7a)
  KMC    =         KMC25 * exp((1/298.-1/(T+273))*EAKMC /R)  ! ppm CO2 Km-value Rubisco for carboxylation (Eqn 7b)
  KMO    =         KMO25 * exp((1/298.-1/(T+273))*EAKMO /R)  ! % O2	Km-value Rubisco for oxygenation (Eqn 7c)
  GAMMAX = 0.5 * KOKC * KMC * O2 / KMO                       ! ppm CO2)	CO2 compensation point at no mitochondrial respiration (Eqn 7d)
  PMAX   = VCMAX * (CO2I-GAMMAX) / (CO2I + KMC * (1+O2/KMO)) ! micromol CO2 m-2 s-1	Photosynthesis rate of upper leaves at light saturation (Eqn 6a)
  TMPFAC = max( 0., min( 1., (T+4.)/5. ) )                   ! Linear decrease of photosynthetic quantum yield at low temperature (below 1 degC)
  EFF    = TMPFAC * (1/2.1) * (CO2I-GAMMAX) / (4.5*CO2I+10.5*GAMMAX) ! mol CO2 mol-1 PAR quanta	Quantum yield of photosynthesis (Eqn 6b)
  LUEMXQ = EFF*PMAX*(1+KLUETILG*(1-FRACTV)) / (EFF*KLAI/BASAL*PARAV + PMAX)   ! mol CO2 mol-1 PAR Light-use efficiency (Eqn 5)
end Subroutine LUECO2TM

! Calculate RESPHARDSI respiration for use in Growth()
Subroutine HardeningSink(CLV,DAYL,doy,LT50,Tsurf)
  integer :: doy
  real :: CLV,DAYL,LT50,Tsurf
  real :: doySinceStart, reHardRedStart
  if ( LAT > 0 ) then ! correct for hemisphere
    reHardRedStart = modulo( reHardRedEnd - reHardRedDay, 365. ) ! Rehardening reduction start
  else
    reHardRedStart = modulo( reHardRedEnd + 183 - reHardRedDay, 365. ) ! Rehardening reduction adjusted for hemisphere
  end if
  doySinceStart  = modulo( doy-reHardRedStart       , 365. )
  if ( doySinceStart < (reHardRedDay+0.5*(365.-reHardRedDay)) ) then
    reHardPeriod = max( 0., 1.-doySinceStart/reHardRedDay )
  else
    reHardPeriod = 1.
  end if
  if ( (Tsurf>THARDMX) .or. (LT50<LT50MN) ) then
    RATEH = 0.
  else
    RATEH = reHardPeriod * Hparam * (THARDMX-Tsurf) * (LT50-LT50MN)
  end if
  RESPHARDSI = RATEH * CLV * KRESPHARD * max(0.,min(1., RESNOR*5. )) ! gC m-2 d-1 Sink strength from carbohydrate demand of hardening
end Subroutine HardeningSink

! Calculate all the growth rates
Subroutine Growth(CLV,CRES,CST,PARINT,TILG2,TILG1,TILV,TRANRF,AGE,LAI, GLV,GRES,GRT,GST)
  real :: CLV,CRES,CST,PARINT,TILG2,TILG1,TILV,TRANRF,AGE,LAI
  real :: GLV,GRES,GRT,GST
!  PHOT     = PARINT * TRANRF * 12. * LUEMXQ * NOHARV               ! gC m-2 d-1 Photosynthesis (12. = gC mol-1)
  PHOT     = PARINT * TRANRF * 12. * LUEMXQ                        ! gC m-2 d-1 Photosynthesis (12. = gC mol-1), Simon removed NOHARV
!  RESMOB   = (CRES * NOHARV / TCRES) * max(0.,min( 1.,DAVTMP/5. )) ! gC m-2 d-1	Mobilisation of reserves
  RESMOB   = max(0.0, CRES - CRESMN) / TCRES * max(0.0, min(1.0, DAVTMP/5.0)) ! gC m-2 d-1	Mobilisation of reserves, Simon removed NOHARV
  SOURCE   = RESMOB + PHOT                                         ! gC m-2 d-1	Source strength from photsynthesis and reserve mobilisation
  RESPHARD = min(SOURCE,RESPHARDSI)                                ! gC m-2 d-1	Plant hardening respiration
  ALLOTOT  = SOURCE - RESPHARD                                     ! gC m-2 d-1	Allocation of carbohydrates to sinks other than hardening
!  GRESSI   = 0.5 * (RESMOB + max(0., CRESMX-CRES) / DELT)         ! gC m-2 d-1 Sink strength of reserve pool (a fraction of CRESMX-(CRES-RESMOB))
  GRESSI   = FGRESSI * max(0., CRESMX-(CRES-RESMOB)) / DELT        ! gC m-2 d-1 Sink strength of reserve pool (a fraction of CRESMX-(CRES-RESMOB)), Simon parameterised
  if (TILG2 > 0.0) then
    CSTAV  = CST/TILG2                                             ! gC tiller-1 Average stem mass of elongating tillers
  else
    CSTAV  = 0.
  end if
  SINK1T   = max(0., 1 - (CSTAV/CSTAVM)) * SIMAX1T                 ! gC tiller-1 d-1 Sink strength of individual elongating tillers
  NELLVG   = PHENRF * NELLVM                                       ! leaves tiller-1 Growing leaves per elongating tiller.
!  GLAISI   = ((LERV*TILV*NELLVM*LFWIDV) + (LERG*TILG2*NELLVG*LFWIDG)) * LSHAPE * TRANRF ! m2 leaf m-2 d-1 Potential growth rate of leaf area
  GLAISI   = ((LERV*(TILV+TILG1)*NELLVM*LFWIDV) + (LERG*TILG2*NELLVG*LFWIDG)) * LSHAPE * TRANRF ! m2 leaf m-2 d-1 Potential growth rate of leaf area (Simon added TILG1)
!  GLVSI    = max(0.0, (GLAISI * NOHARV / SLANEW) / YG)              ! gC m-2 d-1 Potential growth rate of leaf mass
!  GSTSI    = max(0.0, (SINK1T * TILG2 * TRANRF * NOHARV) / YG)      ! gC m-2 d-1 Potential growth rate of stems
  GLVSI    = max(0.0, (GLAISI / SLANEW) / YG)              ! gC m-2 d-1 Potential growth rate of leaf mass, Simon removed NOHARV
  GSTSI    = max(0.0, (SINK1T * TILG2 * TRANRF) / YG)      ! gC m-2 d-1 Potential growth rate of stems, Simon removed NOHARV
  call Allocation(GRES,GRT,GLV,GST)
end Subroutine Growth

   ! Calculate allocation of CRES to GRES,GRT,GLV,GST
   Subroutine Allocation(GRES,GRT,GLV,GST)
     real :: GRES, GRT, GLV, GST
     ! Sinks RESPHARDSI, GLVSI, GSTSI, GRESSI,
     GSHSI = GLVSI + GSTSI
!     if (DAYLGE >= 0.1) then   ! Simon thinks maybe this value should be a parameter
     if (DAYL >= DAYLA) then   ! Simon modified since DAYLGE could remain high
     ! Situation 1: Growth has priority over storage (spring and growth period)
       ! Calculate amount of assimilates allocated to shoot
       ALLOSH = min( ALLOTOT, GSHSI )
       ! Calculate amount of assimilates allocated to reserves
       GRES   = min( ALLOTOT - ALLOSH, GRESSI)
     else
     ! Situation 2: Storage has priority over shoot (autumn)
       ! Calculate amount of assimilates allocated to reserves
       GRES   = min( ALLOTOT, GRESSI )
       ! Calculate amount of assimilates allocated to shoot
       ALLOSH = min( ALLOTOT - GRES, GSHSI )
     end if
     ! All surplus carbohydrate goes to roots
     ALLORT  = ALLOTOT - ALLOSH - GRES
     if (GSHSI == 0.) GSHSI = 1           ! avoid divide by zero error when GSHSI==0.
     ALLOLV  = GLVSI * (ALLOSH / GSHSI)
     ALLOFRAC = ALLOLV / GLVSI            ! Simon fraction of allocation to leaves (non-stem shoot)
     ALLOST  = GSTSI * (ALLOSH / GSHSI)
     GLV     = ALLOLV * YG
     GST     = ALLOST * YG
     GRT     = ALLORT * YG
     RESPGSH = (ALLOLV + ALLOST) * (1-YG) ! gC m-2 d-1 Respiration associated with shoot growth
     RESPGRT =  ALLORT           * (1-YG) ! gC m-2 d-1 Respiration associated with root growth
   end Subroutine Allocation

! Calculate RplantAer = gC m-2 d-1 Aerobic plant respiration
Subroutine PlantRespiration(FO2,RESPHARD)
  real :: FO2,RESPHARD
  real :: fAer
  fAer      = max(0.,min(1., FO2/FO2MX ))
  RplantAer = fAer * ( RESPGRT + RESPGSH + RESPHARD )
end Subroutine PlantRespiration

! Calculate death rates
Subroutine Senescence(CLV,CRT,CSTUB,doy,LAI,PARBASE,BASAL,LT50,PERMgas,TRANRF,TANAER,TILV,Tsurf,AGE, &
                                 DeHardRate,DLAI,DLV,DRT,DSTUB,dTANAER,DTILV,HardRate,RDRS,RDRW)
  integer :: doy
  real :: CLV,CRT,CSTUB,DAYL,LAI,PARBASE,BASAL,LT50,PERMgas,TRANRF,TANAER,TILV,Tsurf,AGE
  real :: DeHardRate,DLAI,DLV,DRT,DSTUB,dTANAER,DTILV,HardRate
  real :: RDRS, TV1, TV2, RDRW
  call AnaerobicDamage(LT50,PERMgas,TANAER, dTANAER)
  call Hardening(CLV,LT50,Tsurf, DeHardRate,HardRate)
!  if (LAI/BASAL < LAICR) then
!    TV1 = 0.0
!  else
!    TV1 = RDRSCO*(LAI/BASAL-LAICR)/LAICR            ! RDRSCO/LAICR is the slope of RDRS past LAICR
!  end if
!  RDRS   = min(TV1, RDRSMX)                         ! d-1 Relative leaf and tiller death rate due to shading, see Gastal & Lemaire 2015
!  RDRS   = max(0.0, min(RDRSCO*(LAI/BASAL-LAICR)/LAICR, RDRSMX)) ! d-1 Relative leaf and tiller death rate due to shading, Original rewritten on one line
!  RDRS   = max(0.0, min(RDRSCO*(LAI/BASAL-LAICR), RDRSMX)) ! d-1 Relative leaf and tiller death rate due to shading, Simon simplified original
  RDRS   = max(0.0, RDRSMX*(1 - exp(-KLAI*LAI/BASAL)/exp(-KLAI*LAICR/BASAL)))  ! d-1 Relative leaf and tiller death rate due to shading, Simon shading method
  RDRT   = max(RDRTMIN, RDRTEM * Tsurf)             ! d-1 Leaf turnover temperature dependent
  RDRW   = RDRWMAX * (1 - TRANRF / TRANRFCR)        ! d-1 Simon Relative death rate due to water stress
! Original
!  TV2    = NOHARV * max(RDRS,RDRT,RDRFROST,RDRTOX) ! d-1 Relative leaf death rate
!  RDRTIL = NOHARV * max(RDRS,     RDRFROST,RDRTOX) ! d-1 Relative death rate of non-elongating tillers
! Simon try different ways to combine death rates to make parameters more responsive
! Maximum stress
  TV2    = max(RDRS,RDRFROST,RDRTOX,RDRT,RDRW,RDRTILMIN)   ! d-1 Relative leaf death rate (can't be smaller than RDRTIL)
  RDRTIL = max(RDRS,RDRFROST,RDRTOX     ,RDRW,RDRTILMIN)   ! d-1 Relative death rate of non-elongating tillers
! Euclidean combination
!  TV2    = sqrt(RDRS*RDRS+RDRW*RDRW+RDRFROST*RDRFROST+RDRTOX*RDRTOX+RDRT*RDRT)           ! d-1 Relative leaf death rate
!  RDRTIL = sqrt(RDRS*RDRS+RDRW*RDRW+RDRFROST*RDRFROST+RDRTOX*RDRTOX+RDRTILMIN*RDRTILMIN) ! d-1 Relative death rate of non-elongating tillers, Simon added background death rate
! Joint survival probability
!  TV2    = 1 - (1-RDRS)*(1-RDRW)*(1-RDRFROST)*(1-RDRTOX)*(1-RDRT)           ! d-1 Relative leaf death rate
!  RDRTIL = 1 - (1-RDRS)*(1-RDRW)*(1-RDRFROST)*(1-RDRTOX)*(1-RDRTILMIN)      ! d-1 Relative death rate of non-elongating tillers, Simon added background death rate
! Additive stress
!  TV2    = RDRS+RDRW+RDRFROST+RDRTOX+RDRT      ! d-1 Relative leaf death rate
!  RDRTIL = RDRS+RDRW+RDRFROST+RDRTOX+RDRTILMIN ! d-1 Relative death rate of non-elongating tillers, Simon added background death rate
  RDRL   = TV2
  DLAI   = LAI    * TV2
  DLV    = CLV    * TV2
  DSTUB  = CSTUB  * RDRSTUB
  DTILV  = TILV   * RDRTIL
  DRT    = CRT    * RDRROOT

end Subroutine Senescence

   ! Calculate RDRTOX = d-1	Relative death rate of tillers due to anaerobic conditions
   Subroutine AnaerobicDamage(LT50,PERMgas,TANAER, dTANAER)
     real :: LT50,PERMgas,TANAER
     real :: dTANAER,LD50
     if (PERMgas==0.) then      ! d-1 Permeability of soil surface to gas exchange
       dTANAER = 1.             ! d d-1	Change in days since start anaerobic conditions
     else
       dTANAER = -TANAER / DELT ! d d-1	Change in days since start anaerobic conditions
     end if
     LD50 = LDT50A + LDT50B * LT50 ! d Duration of anaerobic conditions at which death rate is half the maximum
     if (TANAER > 0.) then      ! d	Time since start anaerobic conditions
       RDRTOX = KRDRANAER / (1.+exp(-KRDRANAER*(TANAER-LD50))) ! d-1 Relative death rate of tillers due to anaerobic conditions
     else
       RDRTOX = 0.
     end if
     end Subroutine AnaerobicDamage

   ! Calculate RDRFROST, DeHardRate, HardRate
   Subroutine Hardening(CLV,LT50,Tsurf, DeHardRate,HardRate)
     real :: CLV,LT50,Tsurf
     real :: DeHardRate,HardRate
     real :: RATED,RSR3H,RSRDAY
     RSR3H      = 1. / (1.+exp(-KRSR3H*(Tsurf-LT50))) ! d-1	Relative frost survival rate
     ! RDRFROST should be less than 1 to avoid numerical problems
     ! (loss of all biomass but keeping positive reserves). We cap it at 0.5.
     RSRDAY     = RSR3H ! In previous versions we had RSRDAY = RSR3H^8 which understimated survival
     RDRFROST   = min( 0.5, 1. - RSRDAY )             ! d-1 Relative death rate due to frost
     RATED      = min( Dparam*(LT50MX-LT50)*(Tsurf+TsurfDiff), (LT50MX-LT50)/DELT ) ! ?C d-1 Potential rate of dehardening, if below limit set by RATEDMX
     DeHardRate = max(0.,min( RATEDMX, RATED ))
     if ( CLV > 0.0 ) then
       HardRate   = RESPHARD / (CLV * KRESPHARD)
     else
       HardRate   = 0.0
     end if
   end Subroutine Hardening

! Simon added decomposition function
! Calculate decompositon of dead leaf
Subroutine Decomposition(CLVD,DAVTMP,WCLM, DLVD,RDLVD)
  real :: CLVD,DAVTMP,WCLM
  real :: DLVD
  real :: PSIA,PSIB,SWCS,PSIS!,DELD,DELE
  real :: EBIOMASS,CT,CP,WORMS
  real :: DTEMP,DWATER,DECOMP,RDLVD
!  EBIOMASSMAX = 131.0              ! g m-2
!  PSIA    = 3.0e-3                 ! Te Kowhai silt loam
!  PSIB    = 7.75                   ! Te Kowhai silt loam
!  BD      = 1.1                    ! Bulk density (Singleton pers comm)
!  DELD    = 0.0148                 ! Decomposition disappearance
!  DELE    = 0.0005                 ! Earthworm disappearance
  SWCS    = WCLM                    ! Volumetric soil water content near surface (WCL = in non-frozen root zone)
  ! PSIFC = -1500 kPa = -PSIA * (WCFC ** (-PSIB))
  ! PSIWP =   -20 kPa = -PSIA * (WCWP ** (-PSIB))
  PSIB    = -log(1500.0/20.0) / log(WCWP/WCFC)
  PSIA    = 20.0 / (WCFC ** (-PSIB))
  PSIS    =  -PSIA * (SWCS ** (-PSIB)) ! Soil water tension near surface
  ! Calculate number of worms and their grazing of dead matter
  ! Numbers at surface based on Baker et al., driven by GWCS
  ! Activity based on Daniels
  EBIOMASS= max(0.0, min(1.0, 5.0*SWCS/BD-1.0)) * EBIOMAX ! EBIOMASSMAX
  if (DAVTMP > 20.0) then
    CT    = 0.0
  else
    CT    = 0.515 * (20.0-DAVTMP) ** 1.84 * exp(-0.297*(20.0-DAVTMP))/2.345 ! Daniels
  end if
  if (PSIS > -12.3) then
    CP    = 1.0
  else
    CP    = 0.549 * (-PSIS) ** 0.793 * exp(0.113 * PSIS) ! Daniels
  end if
  WORMS   = DELE * EBIOMASS * CT * CP
  ! Calculate decomposition, based on Andren paper
  if (DAVTMP > 0.0) then
    DTEMP = 2.0 ** ((DAVTMP - 20.0)/10.0)
  else
    DTEMP = 0.0
  end if
  DWATER  = max(0.0, min(1.0, log(-7580.0 / PSIS) / log(-7580.0 / (-10.0))))
  if (RAIN > 0.0) DWATER = 1.0       ! decomp on rain days even if dry soil, McCall 1984
  DECOMP  = DELD * DTEMP * DWATER  ! total relative decomposition rate
  ! Total relative dead matter disappearance rate
  RDLVD   = DECOMP + WORMS
  DLVD    = CLVD    * RDLVD
end Subroutine Decomposition

! Simon renamed Foliage2() to Tillering()
! Calculate GLAI,GTILV,TILVG1,TILG1G2
Subroutine Tillering(DAYL,GLV,LAI,BASAL,TILV,TILG1,TRANRF,Tsurf,VERN,AGE, GLAI,RGRTV,GTILV,TILVG1,TILG1G2)
  real    :: DAYL,GLV,LAI,BASAL,TILV,TILG1,TRANRF,Tsurf,AGE
!  integer :: VERN
  real :: VERN
  real    :: GLAI,GTILV,TILVG1,TILG1G2
  real    :: RGRTV,RGRTVG1,TV1,TV2
  GLAI    = SLANEW * GLV                                                      ! Note SLANEW is in m2 leaf gC-1
  if (Tsurf < TBASE) then
    TV1   = 0.
  else
    TV1   = Tsurf/PHY                                                         ! d-1 Potential leaf appearance rate
  end if
!  RLEAF   = TV1 * TRANRF * DAYLGE * ( FRACTV + PHENRF * (1-FRACTV) )          ! d-1 Leaf appearance rate, Original
  RLEAF   = TV1 * TRANRF * ( FRACTV + PHENRF * (1-FRACTV) )                   ! d-1 Leaf appearance rate. Simon removed DAYLGE effect (Pararajasingham and Hunt 1995)
!  TV2     = max( 0.0, min(FSMAX, LAITIL - LAIEFT*LAI/BASAL ))                 ! tillers site-1 Ratio of tiller appearance and leaf apearance rates, Original
!  TV2     = min(FSMAX, FSMAX * exp( - LAIEFT * (LAI/BASAL-LAITIL) ))           ! tillers site-1 Ratio of tiller appearance and leaf apearance rates, Simon modifed
  TV2     = min(FSMAX, FSMAX * exp(-KLAI*LAI/BASAL) / exp(-KLAI*LAITIL/BASAL) )! tillers site-1 Ratio of tiller appearance and leaf apearance rates, Simon shading method
  FS      = TV2                                                               ! Simon record site filling fraction
  RGRTV   = max( 0.0       , TV2 * RESNOR * RLEAF )                           ! d-1 Relative rate of vegetative tiller appearance
  GTILV   = TILV  * RGRTV                                                     ! Simon deleted NOHARV switch
  TGE     = max( 0.0       , 1.0 - (abs(DAVTMP - TOPTGE))/(TOPTGE-TBASE))     ! Temperature effect on initiation of elongation in tillers
  RGRTVG1 = DAYLGE * TGE * RGENMX * VERN                                      ! d-1 Relative rate of vegetative tiller conversion to generative, Simon removed NOHARV
  TILVG1  = TILV  * RGRTVG1
  if (DAYL > DAYLG1G2) then                                                   ! Generative tiller elongation controlled by DAYL
    TILG1G2 = TILG1 * RGRTG1G2
!    TILG1G2 = TILG1 * RGRTG1G2 * TGE                                          ! Simon added temperature response
  else if (YDAYL < DAYL) then
    TILG1G2 = 0.                                                              ! no conversion yet
  else
    TILG1G2 = 0.                                                              ! Simon remaining generative tillers remain generative
!    TILG1G2 = TILG1                                                           ! Simon remaining generative tillers elongate
    TILVG1  = TILVG1 - TILG1                                                  ! Simon remaining generative tillers revert to vegetative
  end if
end Subroutine Tillering

  Subroutine Reseed(day, NDAYS, NHARVCOL, DAYS_HARVEST, BASAL, LAI, PHEN, TILG1, TILG2, TILV, & ! inputs
                    CLV, CRES, CST, CSTUB, &
                    RESEEDED) ! outputs
  ! add a re-seed option matt hanson
    integer :: day
    integer :: NDAYS, NHARVCOL
    real, dimension(NDAYS, NHARVCOL) :: DAYS_HARVEST     ! major re-structure by Matt Hanson
    real    :: BASAL, LAI, PHEN, TILG2, TILG1, TILV, CLV, CRES, CST, CSTUB ! values that may be modified.
    real    :: reseed_trig, reseed_basal, RESEEDED

    reseed_trig =  DAYS_HARVEST(day, 7)
    reseed_basal = DAYS_HARVEST(day, 8)
    RESEEDED = 0
    if ((reseed_trig>=0) .and. (BASAL<=reseed_trig)) then ! reseed_trig < 0 is a flag for do not re-seed
      RESEEDED = 1

      ! set parameters:
      BASAL = reseed_basal ! coverage fraction=

      if (reseed_LAI >=0) then
        LAI = reseed_LAI  ! Leaf area index
      end if

      PHEN =0 ! Phenological stage zeroed after re-seed
      if (reseed_TILG2>=0) then
        TILG2 = reseed_TILG2  ! Non-elongating generative tiller density
      end if
      if (reseed_TILG1>=0) then
        TILG1 = reseed_TILG1  ! Elongating generative tiller density
      end if
      if (reseed_TILV>=0) then
        TILV = reseed_TILV  ! Non-elongating tiller density
      end if
      ! set harvest delay by setting harv_trig to -1 for the day and the following days
      DAYS_HARVEST(day: day + reseed_harv_delay, 4) = -1

      ! add the carbon stores! on simon's reccomendations
      if (reseed_CLV>=0) then
        CLV = reseed_CLV ! Weight of leaves

      end if
      if (reseed_CRES>=0) then
        CRES = reseed_CRES  ! Weight of reserves

      end if
      if (reseed_CST>=0) then
        CST = reseed_CST  ! Weight of stems

      end if
      if (reseed_CSTUB>=0) then
        CSTUB = reseed_CSTUB  ! Weight of stubble

      end if

    end if

  End Subroutine Reseed

end module plant


text:
! RUN: %S/test_errors.sh %s %t %f18

! When a module subprogram has the MODULE prefix the following must match
! with the corresponding separate module procedure interface body:
! - C1549: characteristics and dummy argument names
! - C1550: binding label
! - C1551: NON_RECURSIVE prefix

module m1
  interface
    module subroutine s4(x)
      real, intent(in) :: x
    end
    module subroutine s5(x, y)
      real, pointer :: x
      real, value :: y
    end
    module subroutine s6(x, y)
      real :: x
      real :: y
    end
    module subroutine s7(x, y, z)
      real :: x(8)
      real :: y(8)
      real :: z(8)
    end
    module subroutine s8(x, y, z)
      real :: x(8)
      real :: y(*)
      real :: z(*)
    end
    module subroutine s9(x, y, z, w)
      character(len=4) :: x
      character(len=4) :: y
      character(len=*) :: z
      character(len=*) :: w
    end
  end interface
end

submodule(m1) sm1
contains
  module subroutine s4(x)
    !ERROR: The intent of dummy argument 'x' does not match the intent of the corresponding argument in the interface body
    real, intent(out) :: x
  end
  module subroutine s5(x, y)
    !ERROR: Dummy argument 'x' has the OPTIONAL attribute; the corresponding argument in the interface body does not
    real, pointer, optional :: x
    !ERROR: Dummy argument 'y' does not have the VALUE attribute; the corresponding argument in the interface body does
    real :: y
  end
  module subroutine s6(x, y)
    !ERROR: Dummy argument 'x' has type INTEGER(4); the corresponding argument in the interface body has type REAL(4)
    integer :: x
    !ERROR: Dummy argument 'y' has type REAL(8); the corresponding argument in the interface body has type REAL(4)
    real(8) :: y
  end
  module subroutine s7(x, y, z)
    integer, parameter :: n = 8
    real :: x(n)
    real :: y(2:n+1)
    !ERROR: The shape of dummy argument 'z' does not match the shape of the corresponding argument in the interface body
    real :: z(n+1)
  end
  module subroutine s8(x, y, z)
    !ERROR: The shape of dummy argument 'x' does not match the shape of the corresponding argument in the interface body
    real :: x(*)
    real :: y(*)
    !ERROR: The shape of dummy argument 'z' does not match the shape of the corresponding argument in the interface body
    real :: z(8)
  end
  module subroutine s9(x, y, z, w)
    character(len=4) :: x
    !ERROR: Dummy argument 'y' has type CHARACTER(KIND=1,LEN=5_4); the corresponding argument in the interface body has type CHARACTER(KIND=1,LEN=4_4)
    character(len=5) :: y
    character(len=*) :: z
    !ERROR: Dummy argument 'w' has type CHARACTER(KIND=1,LEN=4_4); the corresponding argument in the interface body has type CHARACTER(KIND=1,LEN=*)
    character(len=4) :: w
  end
end

module m2
  interface
    module subroutine s1(x, y)
      real, intent(in) :: x
      real, intent(out) :: y
    end
    module subroutine s2(x, y)
      real, intent(in) :: x
      real, intent(out) :: y
    end
    module subroutine s3(x, y)
      real(4) :: x
      procedure(real) :: y
    end
    module subroutine s4()
    end
    non_recursive module subroutine s5()
    end
  end interface
end

submodule(m2) sm2
contains
  !ERROR: Module subprogram 's1' has 3 args but the corresponding interface body has 2
  module subroutine s1(x, y, z)
    real, intent(in) :: x
    real, intent(out) :: y
    real :: z
  end
  !ERROR: Dummy argument name 'z' does not match corresponding name 'y' in interface body
  module subroutine s2(x, z)
    real, intent(in) :: x
    real, intent(out) :: y
  end
  module subroutine s3(x, y)
    !ERROR: Dummy argument 'x' is a procedure; the corresponding argument in the interface body is not
    procedure(real) :: x
    !ERROR: Dummy argument 'y' is a data object; the corresponding argument in the interface body is not
    real :: y
  end
  !ERROR: Module subprogram 's4' has NON_RECURSIVE prefix but the corresponding interface body does not
  non_recursive module subroutine s4()
  end
  !ERROR: Module subprogram 's5' does not have NON_RECURSIVE prefix but the corresponding interface body does
  module subroutine s5()
  end
end

module m2b
  interface
    module subroutine s1()
    end
    module subroutine s2() bind(c, name="s2")
    end
    module subroutine s3() bind(c, name="s3")
    end
  end interface
end

submodule(m2b) sm2b
  character(*), parameter :: suffix = "_xxx"
contains
  !ERROR: Module subprogram 's1' has a binding label but the corresponding interface body does not
  module subroutine s1() bind(c, name="s1")
  end
  !ERROR: Module subprogram 's2' does not have a binding label but the corresponding interface body does
  module subroutine s2()
  end
  !ERROR: Module subprogram 's3' has binding label "s3_xxx" but the corresponding interface body has "s3"
  module subroutine s3() bind(c, name="s3" // suffix)
  end
end

module m3
  interface
    module subroutine s1(x, y, z)
      procedure(real), pointer, intent(in) :: x
      procedure(real), pointer, intent(out) :: y
      procedure(real), pointer, intent(out) :: z
    end
    module subroutine s2(x, y)
      procedure(real), pointer :: x
      procedure(real) :: y
    end
  end interface
end

submodule(m3) sm3
contains
  module subroutine s1(x, y, z)
    procedure(real), pointer, intent(in) :: x
    !ERROR: The intent of dummy argument 'y' does not match the intent of the corresponding argument in the interface body
    procedure(real), pointer, intent(inout) :: y
    !ERROR: The intent of dummy argument 'z' does not match the intent of the corresponding argument in the interface body
    procedure(real), pointer :: z
  end
  module subroutine s2(x, y)
    !ERROR: Dummy argument 'x' has the OPTIONAL attribute; the corresponding argument in the interface body does not
    !ERROR: Dummy argument 'x' does not have the POINTER attribute; the corresponding argument in the interface body does
    procedure(real), optional :: x
    !ERROR: Dummy argument 'y' has the POINTER attribute; the corresponding argument in the interface body does not
    procedure(real), pointer :: y
  end
end

module m4
  interface
    subroutine s_real(x)
      real :: x
    end
    subroutine s_real2(x)
      real :: x
    end
    subroutine s_integer(x)
      integer :: x
    end
    module subroutine s1(x)
      procedure(s_real) :: x
    end
    module subroutine s2(x)
      procedure(s_real) :: x
    end
  end interface
end

submodule(m4) sm4
contains
  module subroutine s1(x)
    !OK
    procedure(s_real2) :: x
  end
  module subroutine s2(x)
    !ERROR: Dummy procedure 'x' does not match the corresponding argument in the interface body
    procedure(s_integer) :: x
  end
end

module m5
  interface
    module function f1()
      real :: f1
    end
    module subroutine s2()
    end
  end interface
end

submodule(m5) sm5
contains
  !ERROR: Module subroutine 'f1' was declared as a function in the corresponding interface body
  module subroutine f1()
  end
  !ERROR: Module function 's2' was declared as a subroutine in the corresponding interface body
  module function s2()
  end
end

module m6
  interface
    module function f1()
      real :: f1
    end
    module function f2()
      real :: f2
    end
    module function f3()
      real :: f3
    end
  end interface
end

submodule(m6) ms6
contains
  !OK
  real module function f1()
  end
  !ERROR: Return type of function 'f2' does not match return type of the corresponding interface body
  integer module function f2()
  end
  !ERROR: Return type of function 'f3' does not match return type of the corresponding interface body
  module function f3()
    real :: f3
    pointer :: f3
  end
end

module m7
  interface
    module subroutine s1(x, *)
      real :: x
    end
  end interface
end

submodule(m7) sm7
contains
  !ERROR: Dummy argument 1 of 's1' is an alternate return indicator but the corresponding argument in the interface body is not
  !ERROR: Dummy argument 2 of 's1' is not an alternate return indicator but the corresponding argument in the interface body is
  module subroutine s1(*, x)
    real :: x
  end
end

max_stars_repo_name: Machiry/checkedc-clang


text:
type scimage
   integer, dimension(:,:), pointer :: channel
   integer :: width, height
end type scimage

ext: f

max_issues_repo_name: seanwallawalla-forks/RosettaCodeData

max_forks_repo_name: seanwallawalla-forks/RosettaCodeData


text:
      PROGRAM M
      IMPLICIT NONE
      INTEGER Z, A

      Z = 12
      A = 567

      WRITE(*,10) Z, A
 10   FORMAT (I2, I3)
      STOP
      END PROGRAM M

ext: for


text:
      SUBROUTINE genmul(n,p,ncat,ix)
C**********************************************************************
C
C            SUBROUTINE GENMUL( N, P, NCAT, IX )
C     GENerate an observation from the MULtinomial distribution
C
C
C                              Arguments
C
C
C     N --> Number of events that will be classified into one of
C           the categories 1..NCAT
C                         INTEGER N
C
C     P --> Vector of probabilities.  P(i) is the probability that
C           an event will be classified into category i.  Thus, P(i)
C           must be [0,1]. Only the first NCAT-1 P(i) must be defined
C           since P(NCAT) is 1.0 minus the sum of the first
C           NCAT-1 P(i).
C                         REAL P(NCAT-1)
C
C     NCAT --> Number of categories.  Length of P and IX.
C                         INTEGER NCAT
C
C     IX <-- Observation from multinomial distribution.  All IX(i)
C            will be nonnegative and their sum will be N.
C                         INTEGER IX(NCAT)
C
C
C                              Method
C
C
C     Algorithm from page 559 of
C
C     Devroye, Luc
C
C     Non-Uniform Random Variate Generation.  Springer-Verlag,
C     New York, 1986.
C
C**********************************************************************
C     .. Scalar Arguments ..
      INTEGER n,ncat
C     ..
C     .. Array Arguments ..
      REAL p(*)
      INTEGER ix(*)
C     ..
C     .. Local Scalars ..
      REAL prob,ptot,sum
      INTEGER i,icat,ntot
C     ..
C     .. External Functions ..
      INTEGER ignbin
      EXTERNAL ignbin
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC abs
C     ..
C     .. Executable Statements ..

C     Check Arguments
      IF (n.LT.0) STOP 'N < 0 in GENMUL'
      IF (ncat.LE.1) STOP 'NCAT <= 1 in GENMUL'
      ptot = 0.0
      DO 10,i = 1,ncat - 1
          IF (p(i).LT.0.0) STOP 'Some P(i) < 0 in GENMUL'
          IF (p(i).GT.1.0) STOP 'Some P(i) > 1 in GENMUL'
          ptot = ptot + p(i)
   10 CONTINUE
      IF (ptot.GT.0.99999) STOP 'Sum of P(i) > 1 in GENMUL'

C     Initialize variables
      ntot = n
      sum = 1.0
      DO 20,i = 1,ncat
          ix(i) = 0
   20 CONTINUE

C     Generate the observation
      DO 30,icat = 1,ncat - 1
          prob = p(icat)/sum
          ix(icat) = ignbin(ntot,prob)
          ntot = ntot - ix(icat)
          IF (ntot.LE.0) RETURN
          sum = sum - p(icat)
   30 CONTINUE
      ix(ncat) = ntot

C     Finished
      RETURN

      END

max_stars_repo_name: ppernot/MC-ChemDB

max_issues_repo_name: ppernot/MC-ChemDB

max_forks_repo_name: ppernot/MC-ChemDB


text:
module multifab_coefbc_module

  use multifab_module
  use define_bc_module
  use bc_module

  implicit none

  private

  public :: multifab_coefbc

contains

  subroutine multifab_coefbc(s,start_scomp,ncomp,the_bc_level)

    ! this fills ghost cells for transport coefficients at walls by
    ! averaging ghost and interior into ghost

    type(multifab) , intent(inout) :: s
    integer        , intent(in   ) :: start_scomp,ncomp
    type(bc_level) , intent(in   ) :: the_bc_level

    ! Local
    integer                  :: lo(get_dim(s)),hi(get_dim(s))
    integer                  :: i,ng,dm,comp
    real(kind=dp_t), pointer :: sp(:,:,:,:)

    type(bl_prof_timer),save :: bpt

    call build(bpt,"multifab_coefbc")

    ng = nghost(s)
    dm = get_dim(s)

    do i=1,nfabs(s)
       sp => dataptr(s,i)
       lo = lwb(get_box(s,i))
       hi = upb(get_box(s,i))
       do comp=start_scomp,start_scomp+ncomp-1
          select case (dm)
          case (2)
             call coefbc_2d(sp(:,:,1,comp), lo, hi, ng, &
                            the_bc_level%phys_bc_level_array(i,:,:))
          case (3)
             call coefbc_3d(sp(:,:,:,comp), lo, hi, ng, &
                            the_bc_level%phys_bc_level_array(i,:,:))
          end select
       end do
    end do

    call destroy(bpt)

  end subroutine multifab_coefbc

  subroutine coefbc_2d(s,lo,hi,ng,bc)

    integer        , intent(in   ) :: lo(:),hi(:),ng
    real(kind=dp_t), intent(inout) ::    s(lo(1)-ng:,lo(2)-ng:)
    integer        , intent(in   ) :: bc(:,:)

    ! Local variables
    integer :: i,j

!!!!!!!!!!!!!!!!!!
! lo-x boundary
!!!!!!!!!!!!!!!!!!

    if (bc(1,1) .eq. SLIP_WALL .or. bc(1,1) .eq. NO_SLIP_WALL) then
       ! average interior and ghost into ghost
       do j=lo(2)-ng,hi(2)+ng
          s(lo(1)-ng:lo(1)-1,j) = 0.5d0*(s(lo(1),j)+s(lo(1)-1,j))
       end do
    end if

!!!!!!!!!!!!!!!!!!
! hi-x boundary
!!!!!!!!!!!!!!!!!!

    if (bc(1,2) .eq. SLIP_WALL .or. bc(1,2) .eq. NO_SLIP_WALL) then
       ! average interior and ghost into ghost
       do j=lo(2)-ng,hi(2)+ng
          s(hi(1)+1:hi(1)+ng,j) = 0.5d0*(s(hi(1),j)+s(hi(1)+1,j))
       end do
    end if

!!!!!!!!!!!!!!!!!!
! lo-y boundary
!!!!!!!!!!!!!!!!!!

    if (bc(2,1) .eq. SLIP_WALL .or. bc(2,1) .eq. NO_SLIP_WALL) then
       ! average interior and ghost into ghost
       do i=lo(1)-ng,hi(1)+ng
          s(i,lo(2)-ng:lo(2)-1) = 0.5d0*(s(i,lo(2))+s(i,lo(2)-1))
       end do
    end if

!!!!!!!!!!!!!!!!!!
! hi-y boundary
!!!!!!!!!!!!!!!!!!

    if (bc(2,2) .eq. SLIP_WALL .or. bc(2,2) .eq. NO_SLIP_WALL) then
       ! average interior and ghost into ghost
       do i=lo(1)-ng,hi(1)+ng
          s(i,hi(2)+1:hi(2)+ng) = 0.5d0*(s(i,hi(2))+s(i,hi(2)+1))
       end do
    end if

  end subroutine coefbc_2d

  subroutine coefbc_3d(s,lo,hi,ng,bc)

    integer        , intent(in   ) :: lo(:),hi(:),ng
    real(kind=dp_t), intent(inout) ::    s(lo(1)-ng:,lo(2)-ng:,lo(3)-ng:)
    integer        , intent(in   ) :: bc(:,:)

    ! Local variables
    integer :: i,j,k

!!!!!!!!!!!!!!!!!!
! lo-x boundary
!!!!!!!!!!!!!!!!!!

    if (bc(1,1) .eq. SLIP_WALL .or. bc(1,1) .eq. NO_SLIP_WALL) then
       ! average interior and ghost into ghost
       do k=lo(3)-ng,hi(3)+ng
          do j=lo(2)-ng,hi(2)+ng
             s(lo(1)-ng:lo(1)-1,j,k) = 0.5d0*(s(lo(1),j,k)+s(lo(1)-1,j,k))
          end do
       end do
    end if

!!!!!!!!!!!!!!!!!!
! hi-x boundary
!!!!!!!!!!!!!!!!!!

    if (bc(1,2) .eq. SLIP_WALL .or. bc(1,2) .eq. NO_SLIP_WALL) then
       ! average interior and ghost into ghost
       do k=lo(3)-ng,hi(3)+ng
          do j=lo(2)-ng,hi(2)+ng
             s(hi(1)+1:hi(1)+ng,j,k) = 0.5d0*(s(hi(1),j,k)+s(hi(1)+1,j,k))
          end do
       end do
    end if

!!!!!!!!!!!!!!!!!!
! lo-y boundary
!!!!!!!!!!!!!!!!!!

    if (bc(2,1) .eq. SLIP_WALL .or. bc(2,1) .eq. NO_SLIP_WALL) then
       ! average interior and ghost into ghost
       do k=lo(3)-ng,hi(3)+ng
          do i=lo(1)-ng,hi(1)+ng
             s(i,lo(2)-ng:lo(2)-1,k) = 0.5d0*(s(i,lo(2),k)+s(i,lo(2)-1,k))
          end do
       end do
    end if

!!!!!!!!!!!!!!!!!!
! hi-y boundary
!!!!!!!!!!!!!!!!!!

    if (bc(2,2) .eq. SLIP_WALL .or. bc(2,2) .eq. NO_SLIP_WALL) then
       ! average interior and ghost into ghost
       do k=lo(3)-ng,hi(3)+ng
          do i=lo(1)-ng,hi(1)+ng
             s(i,hi(2)+1:hi(2)+ng,k) = 0.5d0*(s(i,hi(2),k)+s(i,hi(2)+1,k))
          end do
       end do
    end if

!!!!!!!!!!!!!!!!!!
! lo-z boundary
!!!!!!!!!!!!!!!!!!

    if (bc(3,1) .eq. SLIP_WALL .or. bc(3,1) .eq. NO_SLIP_WALL) then
       ! average interior and ghost into ghost
       do j=lo(2)-ng,hi(2)+ng
          do i=lo(1)-ng,hi(1)+ng
             s(i,j,lo(3)-ng:lo(3)-1) = 0.5d0*(s(i,j,lo(3))+s(i,j,lo(3)-1))
          end do
       end do
    end if

!!!!!!!!!!!!!!!!!!
! hi-z boundary
!!!!!!!!!!!!!!!!!!

    if (bc(3,2) .eq. SLIP_WALL .or. bc(3,2) .eq. NO_SLIP_WALL) then
       ! average interior and ghost into ghost
       do j=lo(2)-ng,hi(2)+ng
          do i=lo(1)-ng,hi(1)+ng
             s(i,j,hi(3)+1:hi(3)+ng) = 0.5d0*(s(i,j,hi(3))+s(i,j,hi(3)+1))
          end do
       end do
    end if

  end subroutine coefbc_3d

end module multifab_coefbc_module

max_stars_repo_name: BoxLib-Codes/LowMachMixtures

max_issues_repo_name: BoxLib-Codes/LowMachMixtures

max_forks_repo_name: BoxLib-Codes/LowMachMixtures


text:
MODULE spline_module

  implicit none

  PRIVATE
  PUBLIC :: spline, splint

CONTAINS

  SUBROUTINE spline( x,y,n,yp1,ypn,y2 )
    implicit none
    integer,intent(IN)  :: n
    real(8),intent(IN)  :: x(n),y(n),yp1,ypn
    real(8),intent(OUT) :: y2(n)
    real(8),allocatable :: a(:),b(:),c(:),f(:)
    integer :: i
    allocate( a(n), b(n), c(n), f(n) )
    a=0.0d0 ; b=0.0d0 ; c=0.0d0 ; f=0.0d0
    do i=2,n-1
       a(i) = x(i) - x(i-1)
    end do
    do i=2,n-1
       b(i) = ( x(i+1) - x(i-1) )*2.0d0
    end do
    do i=2,n-2
       c(i) = x(i+1) - x(i)
    end do
    do i=2,n-1
       f(i) = 6.0d0*( (y(i+1)-y(i))/(x(i+1)-x(i)) &
                    - (y(i)-y(i-1))/(x(i)-x(i-1)) )
    end do
    do i=3,n-1
       b(i) = b(i) - a(i)*c(i-1)/b(i-1)
       f(i) = f(i) - a(i)*f(i-1)/b(i-1)
    end do
    y2(:)=0.0d0
    y2(n-1) = f(n-1)/b(n-1)
    do i=n-2,2,-1
       y2(i) = ( f(i) - c(i)*y2(i+1) )/b(i)
    end do
    deallocate( f, c, b, a )
  END SUBROUTINE spline

  SUBROUTINE splint( xa,ya,y2a,n,x,y )
    implicit none
    integer,intent(IN) :: n
    real(8),intent(IN) :: xa(n),ya(n),y2a(n),x
    real(8),intent(OUT) :: y
    integer :: k,khi,klo
    real(8) :: a,b,c,d,h
    klo = 1
    khi = n
1   if ( khi - klo > 1 ) then
       k = ( khi + klo )/2
       if ( xa(k) > x ) then
          khi = k
       else
          klo = k
       end if
       goto 1
    end if
    h = xa(khi) - xa(klo)
    if ( h == 0.0d0 ) stop "splint"
    a = ya(klo)
    b = ( ya(khi)-ya(klo) )/h - ( y2a(khi) + 2.0d0*y2a(klo) )*h/6.0d0
    c = y2a(klo)/2.0d0
    d = ( y2a(khi)-y2a(klo) )/(6.0d0*h)
    h = x - xa(klo)
    y = a + b*h + c*h**2 + d*h**3
  END SUBROUTINE splint

!
! from NUMERICAL RECIPES
!
  SUBROUTINE spline_nr(x,y,n,yp1,ypn,y2)
    implicit none
    integer :: n
    real(8) :: x(n),y(n),yp1,ypn,y2(n)
    integer,parameter :: nmax=50000
    integer :: i,k
    real(8) :: p,qn,sig,un,u(nmax)
    integer,save :: unit=300

    if ( yp1 > 1.d30 ) then
       y2(1) = 0.d0
       u(1) = 0.d0
    else
       y2(1) = -0.5d0
       u(1) = ( 3.d0/(x(2)-x(1)) )*( (y(2)-y(1))/(x(2)-x(1)) - yp1 )
    end if

    do i=2,n-1
       sig = ( x(i)-x(i-1) )/( x(i+1)-x(i-1) )
       p = sig*y2(i-1)+2.d0
       y2(i) = ( sig-1.d0 )/p
       u(i) = ( 6.d0*( (y(i+1)-y(i))/(x(i+1)-x(i)) &
            -(y(i)-y(i-1))/(x(i)-x(i-1)) )/(x(i+1)-x(i-1))-sig*u(i-1) )/p
    end do

    if ( ypn > 1.d30 ) then
       qn = 0.d0
       un = 0.d0
    else
       qn = 0.5d0
       un = ( 3.d0/(x(n)-x(n-1)) )*( ypn - (y(n)-y(n-1))/(x(n)-x(n-1)) )
    end if
    y2(n) = ( un - qn*u(n-1) )/( qn*y2(n-1) + 1.d0 )
    do k=n-1,1,-1
       y2(k) = y2(k)*y2(k+1) + u(k)
    end do

  END SUBROUTINE spline_nr

  SUBROUTINE splint_nr(xa,ya,y2a,n,x,y)
    implicit none
    integer :: n
    real(8) :: x,y,xa(n),y2a(n),ya(n)
    integer :: k,khi,klo
    real(8) :: a,b,h
    klo=1
    khi=n
1   if ( khi-klo > 1 ) then
       k=(khi+klo)/2
       if ( xa(k) > x ) then
          khi=k
       else
          klo=k
       end if
       goto 1
    end if
    h=xa(khi)-xa(klo)
    if ( h == 0.d0 ) stop "splint"
    a=(xa(khi)-x)/h
    b=(x-xa(klo))/h
    y=a*ya(klo)+b*ya(khi)+( (a**3-a)*y2a(klo)+(b**3-b)*y2a(khi) )*(h**2)/6.d0
  END SUBROUTINE splint_nr

END MODULE spline_module

max_stars_repo_name: j-iwata/RSDFT

max_issues_repo_name: j-iwata/RSDFT

max_forks_repo_name: j-iwata/RSDFT


text:
      subroutine dbdqr(n, D, E, c1, c2)
      implicit none
      integer n
      double precision D(*),E(*),c1,c2

      integer i
      double precision cs,sn,r

      if (n.lt.2) return
      do i=1,n-1
         call dlartg(d(i),e(i),cs,sn,r)
         d(i) = r
         e(i) = sn*d(i+1)
         d(i+1) = cs*d(i+1)
      enddo
      call dlartg(d(n),e(n),cs,sn,r)
      d(n) = r
      e(n) = 0.0
      c1 = sn
      c2 = cs
      end


text:
ALLSEL
!/NOPR
NSEL, NONE

startX_coord = 0
*DO, SleeperIndex, 1, AantalSleepers
	X_coord_min = startX_coord - S_Halve_Breedte
	X_coord_max = startX_coord + S_Halve_Breedte

	*USE, '%Scriptfolder(1)%/Contact/createSleeperRughellingplaatContact.MAC', X_coord_min, X_coord_max

	*IF, Railpad_thickness, GT, 0.000001, THEN
		*USE, '%Scriptfolder(1)%/Contact/createRughellingplaatRailpadContact.MAC', X_coord_min, X_coord_max
	*ELSE !Geen railpad
		*USE, '%Scriptfolder(1)%/Contact/createRughellingplaatRailContact.MAC', X_coord_min, X_coord_max
	*ENDIF

	*IF, SleeperIndex, GT, 1, AND, MakeUseOfSymmetry, NE, 'true', THEN

		*USE, '%Scriptfolder(1)%/Contact/createSleeperRughellingplaatContact.MAC', -X_coord_max, -X_coord_min

		*IF, Railpad_thickness, GT, 0.000001, THEN
			*USE, '%Scriptfolder(1)%/Contact/createRughellingplaatRailpadContact.MAC', -X_coord_max, -X_coord_min
		*ELSE !Geen railpad
			*USE, '%Scriptfolder(1)%/Contact/createRughellingplaatRailContact.MAC', -X_coord_max, -X_coord_min
		*ENDIF
	*ENDIF

	startX_coord = SleeperDistance*SleeperIndex !update start coordinate
*ENDDO

CSYS, 0

*DEL, SleeperIndex
*DEL, X_coord_min
*DEL, X_coord_max
*DEL, startX_coord

max_stars_repo_name: dryjognos/AnsysSleeperModel

max_issues_repo_name: dryjognos/AnsysSleeperModel

max_forks_repo_name: dryjognos/AnsysSleeperModel


text:
      SUBROUTINE GRDFCN(N,X,G,NPROB)                                    00000010
      INTEGER N,NPROB
      DOUBLE PRECISION X(N),G(N)
C     **********
C
C     SUBROUTINE GRDFCN
C
C     THIS SUBROUTINE DEFINES THE GRADIENT VECTORS OF EIGHTEEN
C     NONLINEAR UNCONSTRAINED MINIMIZATION PROBLEMS. THE PROBLEM
C     DIMENSIONS ARE AS DESCRIBED IN THE PROLOGUE COMMENTS OF OBJFCN.
C
C     THE SUBROUTINE STATEMENT IS
C
C       SUBROUTINE GRDFCN(N,X,G,NPROB)
C
C     WHERE
C
C       N IS A POSITIVE INTEGER INPUT VARIABLE.
C
C       X IS AN INPUT ARRAY OF LENGTH N.
C
C       G IS AN OUTPUT ARRAY OF LENGTH N WHICH CONTAINS THE COMPONENTS
C         OF THE GRADIENT VECTOR OF THE NPROB OBJECTIVE FUNCTION
C         EVALUATED AT X.
C
C       NPROB IS A POSITIVE INTEGER INPUT VARIABLE WHICH DEFINES THE
C         NUMBER OF THE PROBLEM. NPROB MUST NOT EXCEED 18.
C
C     SUBPROGRAMS CALLED
C
C       FORTRAN-SUPPLIED ... DABS,DATAN,DCOS,DEXP,DLOG,DSIGN,DSIN,
C                            DSQRT
C
C     MINPACK. VERSION OF JULY 1978.
C     BURTON S. GARBOW, KENNETH E. HILLSTROM, JORGE J. MORE
C
C     **********
      INTEGER I,IEV,IVAR,J
      DOUBLE PRECISION AP,ARG,BP,C2PDM6,CP0001,CP1,CP2,CP25,CP5,C1P5,
     1                 C2P25,C2P625,C3P5,C19P8,C20P2,C25,C29,C100,
     2                 C180,C200,C10000,C1PD6,D1,D2,EIGHT,FIFTY,FIVE,
     3                 FOUR,ONE,R,S1,S2,S3,T,T1,T2,T3,TEN,TH,THREE,
     4                 TPI,TWENTY,TWO,ZERO
      DOUBLE PRECISION FVEC(50),Y(15)
      DOUBLE PRECISION DFLOAT
      DATA ZERO,ONE,TWO,THREE,FOUR,FIVE,EIGHT,TEN,TWENTY,FIFTY
     1     /0.0D0,1.0D0,2.0D0,3.0D0,4.0D0,5.0D0,8.0D0,1.0D1,2.0D1,
     2      5.0D1/
      DATA C2PDM6,CP0001,CP1,CP2,CP25,CP5,C1P5,C2P25,C2P625,C3P5,
     1     C19P8,C20P2,C25,C29,C100,C180,C200,C10000,C1PD6
     2     /2.0D-6,1.0D-4,1.0D-1,2.0D-1,2.5D-1,5.0D-1,1.5D0,2.25D0,
     3      2.625D0,3.5D0,1.98D1,2.02D1,2.5D1,2.9D1,1.0D2,1.8D2,2.0D2,
     4      1.0D4,1.0D6/
      DATA AP,BP /1.0D-5,1.0D0/
      DATA Y(1),Y(2),Y(3),Y(4),Y(5),Y(6),Y(7),Y(8),Y(9),Y(10),Y(11),
     1     Y(12),Y(13),Y(14),Y(15)
     2     /9.0D-4,4.4D-3,1.75D-2,5.4D-2,1.295D-1,2.42D-1,3.521D-1,
     3      3.989D-1,3.521D-1,2.42D-1,1.295D-1,5.4D-2,1.75D-2,4.4D-3,
     4      9.0D-4/
      DFLOAT(IVAR) = IVAR
C
C     GRADIENT ROUTINE SELECTOR.
C
      GO TO (10,20,50,70,80,100,130,190,220,260,270,290,310,350,370,
     1       390,400,410), NPROB
C
C     HELICAL VALLEY FUNCTION.
C
   10 CONTINUE
      TPI = EIGHT*DATAN(ONE)
      TH = DSIGN(CP25,X(2))
      IF (X(1) .GT. ZERO) TH = DATAN(X(2)/X(1))/TPI
      IF (X(1) .LT. ZERO) TH = DATAN(X(2)/X(1))/TPI + CP5
      ARG = X(1)**2 + X(2)**2
      R = DSQRT(ARG)
      T = X(3) - TEN*TH
      S1 = TEN*T/(TPI*ARG)
      G(1) = C200*(X(1) - X(1)/R + X(2)*S1)
      G(2) = C200*(X(2) - X(2)/R - X(1)*S1)
      G(3) = TWO*(C100*T + X(3))
      GO TO 490
C
C     BIGGS EXP6 FUNCTION.
C
   20 CONTINUE
      DO 30 J = 1, N
         G(J) = ZERO
   30    CONTINUE
      DO 40 I = 1, 13
         D1 = DFLOAT(I)/TEN
         D2 = DEXP(-D1) - FIVE*DEXP(-TEN*D1) + THREE*DEXP(-FOUR*D1)
         S1 = DEXP(-D1*X(1))
         S2 = DEXP(-D1*X(2))
         S3 = DEXP(-D1*X(5))
         T = X(3)*S1 - X(4)*S2 + X(6)*S3 - D2
         TH = D1*T
         G(1) = G(1) - S1*TH
         G(2) = G(2) + S2*TH
         G(3) = G(3) + S1*T
         G(4) = G(4) - S2*T
         G(5) = G(5) - S3*TH
         G(6) = G(6) + S3*T
   40    CONTINUE
      G(1) = TWO*X(3)*G(1)
      G(2) = TWO*X(4)*G(2)
      G(3) = TWO*G(3)
      G(4) = TWO*G(4)
      G(5) = TWO*X(6)*G(5)
      G(6) = TWO*G(6)
      GO TO 490
C
C     GAUSSIAN FUNCTION.
C
   50 CONTINUE
      G(1) = ZERO
      G(2) = ZERO
      G(3) = ZERO
      DO 60 I = 1, 15
         D1 = CP5*DFLOAT(I-1)
         D2 = C3P5 - D1 - X(3)
         ARG = -CP5*X(2)*D2**2
         R = DEXP(ARG)
         T = X(1)*R - Y(I)
         S1 = R*T
         S2 = D2*S1
         G(1) = G(1) + S1
         G(2) = G(2) - D2*S2
         G(3) = G(3) + S2
   60    CONTINUE
      G(1) = TWO*G(1)
      G(2) = X(1)*G(2)
      G(3) = TWO*X(1)*X(2)*G(3)
      GO TO 490
C
C     POWELL BADLY SCALED FUNCTION.
C
   70 CONTINUE
      T1 = C10000*X(1)*X(2) - ONE
      S1 = DEXP(-X(1))
      S2 = DEXP(-X(2))
      T2 = S1 + S2 - ONE - CP0001
      G(1) = TWO*(C10000*X(2)*T1 - S1*T2)
      G(2) = TWO*(C10000*X(1)*T1 - S2*T2)
      GO TO 490
C
C     BOX 3-DIMENSIONAL FUNCTION.
C
   80 CONTINUE
      G(1) = ZERO
      G(2) = ZERO
      G(3) = ZERO
      DO 90 I = 1, 10
         D1 = DFLOAT(I)
         D2 = D1/TEN
         S1 = DEXP(-D2*X(1))
         S2 = DEXP(-D2*X(2))
         S3 = DEXP(-D2) - DEXP(-D1)
         T = S1 - S2 - S3*X(3)
         TH = D2*T
         G(1) = G(1) - S1*TH
         G(2) = G(2) + S2*TH
         G(3) = G(3) - S3*T
   90    CONTINUE
      G(1) = TWO*G(1)
      G(2) = TWO*G(2)
      G(3) = TWO*G(3)
      GO TO 490
C
C     VARIABLY DIMENSIONED FUNCTION.
C
  100 CONTINUE
      T1 = ZERO
      DO 110 J = 1, N
         T1 = T1 + DFLOAT(J)*(X(J) - ONE)
  110    CONTINUE
      T = T1*(ONE + TWO*T1**2)
      DO 120 J = 1, N
         G(J) = TWO*(X(J) - ONE + DFLOAT(J)*T)
  120    CONTINUE
      GO TO 490
C
C     WATSON FUNCTION.
C
  130 CONTINUE
      DO 140 J = 1, N
         G(J) = ZERO
  140    CONTINUE
      DO 180 I = 1, 29
         D1 = DFLOAT(I)/C29
         S1 = ZERO
         D2 = ONE
         DO 150 J = 2, N
            S1 = S1 + DFLOAT(J-1)*D2*X(J)
            D2 = D1*D2
  150       CONTINUE
         S2 = ZERO
         D2 = ONE
         DO 160 J = 1, N
            S2 = S2 + D2*X(J)
            D2 = D1*D2
  160       CONTINUE
         T = S1 - S2**2 - ONE
         S3 = TWO*D1*S2
         D2 = TWO/D1
         DO 170 J = 1, N
            G(J) = G(J) + D2*(DFLOAT(J-1) - S3)*T
            D2 = D1*D2
  170       CONTINUE
  180    CONTINUE
      T1 = X(2) - X(1)**2 - ONE
      G(1) = G(1) + X(1)*(TWO - FOUR*T1)
      G(2) = G(2) + TWO*T1
      GO TO 490
C
C     PENALTY FUNCTION I.
C
  190 CONTINUE
      T1 = -CP25
      DO 200 J = 1, N
         T1 = T1 + X(J)**2
  200    CONTINUE
      D1 = TWO*AP
      TH = FOUR*BP*T1
      DO 210 J = 1, N
         G(J) = D1*(X(J) - ONE) + X(J)*TH
  210    CONTINUE
      GO TO 490
C
C     PENALTY FUNCTION II.
C
  220 CONTINUE
      T1 = -ONE
      DO 230 J = 1, N
         T1 = T1 + DFLOAT(N-J+1)*X(J)**2
  230    CONTINUE
      D1 = DEXP(CP1)
      D2 = ONE
      TH = FOUR*BP*T1
      DO 250 J = 1, N
         G(J) = DFLOAT(N-J+1)*X(J)*TH
         S1 = DEXP(X(J)/TEN)
         IF (J .EQ. 1) GO TO 240
         S3 = S1 + S2 - D2*(D1 + ONE)
         G(J) = G(J) + AP*S1*(S3 + S1 - ONE/D1)/FIVE
         G(J-1) = G(J-1) + AP*S2*S3/FIVE
  240    CONTINUE
         S2 = S1
         D2 = D1*D2
  250    CONTINUE
      G(1) = G(1) + TWO*BP*(X(1) - CP2)
      GO TO 490
C
C     BROWN BADLY SCALED FUNCTION.
C
  260 CONTINUE
      T1 = X(1) - C1PD6
      T2 = X(2) - C2PDM6
      T3 = X(1)*X(2) - TWO
      G(1) = TWO*(T1 + X(2)*T3)
      G(2) = TWO*(T2 + X(1)*T3)
      GO TO 490
C
C     BROWN AND DENNIS FUNCTION.
C
  270 CONTINUE
      G(1) = ZERO
      G(2) = ZERO
      G(3) = ZERO
      G(4) = ZERO
      DO 280 I = 1, 20
         D1 = DFLOAT(I)/FIVE
         D2 = DSIN(D1)
         T1 = X(1) + D1*X(2) - DEXP(D1)
         T2 = X(3) + D2*X(4) - DCOS(D1)
         T = T1**2 + T2**2
         S1 = T1*T
         S2 = T2*T
         G(1) = G(1) + S1
         G(2) = G(2) + D1*S1
         G(3) = G(3) + S2
         G(4) = G(4) + D2*S2
  280    CONTINUE
      G(1) = FOUR*G(1)
      G(2) = FOUR*G(2)
      G(3) = FOUR*G(3)
      G(4) = FOUR*G(4)
      GO TO 490
C
C     GULF RESEARCH AND DEVELOPMENT FUNCTION.
C
  290 CONTINUE
      G(1) = ZERO
      G(2) = ZERO
      G(3) = ZERO
      D1 = TWO/THREE
      DO 300 I = 1, 99
         ARG = DFLOAT(I)/C100
         R = DABS((-FIFTY*DLOG(ARG))**D1+C25-X(2))
         T1 = R**X(3)/X(1)
         T2 = DEXP(-T1)
         T = T2 - ARG
         S1 = T1*T2*T
         G(1) = G(1) + S1
         G(2) = G(2) + S1/R
         G(3) = G(3) - S1*DLOG(R)
  300    CONTINUE
      G(1) = TWO*G(1)/X(1)
      G(2) = TWO*X(3)*G(2)
      G(3) = TWO*G(3)
      GO TO 490
C
C     TRIGONOMETRIC FUNCTION.
C
  310 CONTINUE
      S1 = ZERO
      DO 320 J = 1, N
         G(J) = DCOS(X(J))
         S1 = S1 + G(J)
  320    CONTINUE
      S2 = ZERO
      DO 330 J = 1, N
         TH = DSIN(X(J))
         T = DFLOAT(N+J) - TH - S1 - DFLOAT(J)*G(J)
         S2 = S2 + T
         G(J) = (DFLOAT(J)*TH - G(J))*T
  330    CONTINUE
      DO 340 J = 1, N
         G(J) = TWO*(G(J) + DSIN(X(J))*S2)
  340    CONTINUE
      GO TO 490
C
C     EXTENDED ROSENBROCK FUNCTION.
C
  350 CONTINUE
      DO 360 J = 1, N, 2
         T1 = ONE - X(J)
         G(J+1) = C200*(X(J+1) - X(J)**2)
         G(J) = -TWO*(X(J)*G(J+1) + T1)
  360    CONTINUE
      GO TO 490
C
C     EXTENDED POWELL FUNCTION.
C
  370 CONTINUE
      DO 380 J = 1, N, 4
         T = X(J) + TEN*X(J+1)
         T1 = X(J+2) - X(J+3)
         S1 = FIVE*T1
         T2 = X(J+1) - TWO*X(J+2)
         S2 = FOUR*T2**3
         T3 = X(J) - X(J+3)
         S3 = TWENTY*T3**3
         G(J) = TWO*(T + S3)
         G(J+1) = TWENTY*T + S2
         G(J+2) = TWO*(S1 - S2)
         G(J+3) = -TWO*(S1 + S3)
  380    CONTINUE
      GO TO 490
C
C     BEALE FUNCTION.
C
  390 CONTINUE
      S1 = ONE - X(2)
      T1 = C1P5 - X(1)*S1
      S2 = ONE - X(2)**2
      T2 = C2P25 - X(1)*S2
      S3 = ONE - X(2)**3
      T3 = C2P625 - X(1)*S3
      G(1) = -TWO*(S1*T1 + S2*T2 + S3*T3)
      G(2) = TWO*X(1)*(T1 + X(2)*(TWO*T2 + THREE*X(2)*T3))
      GO TO 490
C
C     WOOD FUNCTION.
C
  400 CONTINUE
      S1 = X(2) - X(1)**2
      S2 = ONE - X(1)
      S3 = X(2) - ONE
      T1 = X(4) - X(3)**2
      T2 = ONE - X(3)
      T3 = X(4) - ONE
      G(1) = -TWO*(C200*X(1)*S1 + S2)
      G(2) = C200*S1 + C20P2*S3 + C19P8*T3
      G(3) = -TWO*(C180*X(3)*T1 + T2)
      G(4) = C180*T1 + C20P2*T3 + C19P8*S3
      GO TO 490
C
C     CHEBYQUAD FUNCTION.
C
  410 CONTINUE
      DO 420 I = 1, N
         FVEC(I) = ZERO
  420    CONTINUE
      DO 440 J = 1, N
         T1 = ONE
         T2 = TWO*X(J) - ONE
         T = TWO*T2
         DO 430 I = 1, N
            FVEC(I) = FVEC(I) + T2
            TH = T*T2 - T1
            T1 = T2
            T2 = TH
  430       CONTINUE
  440    CONTINUE
      D1 = ONE/DFLOAT(N)
      IEV = -1
      DO 450 I = 1, N
         FVEC(I) = D1*FVEC(I)
         IF (IEV .GT. 0) FVEC(I) = FVEC(I) + ONE/(DFLOAT(I)**2 - ONE)
         IEV = -IEV
  450    CONTINUE
      DO 470 J = 1, N
         G(J) = ZERO
         T1 = ONE
         T2 = TWO*X(J) - ONE
         T = TWO*T2
         S1 = ZERO
         S2 = TWO
         DO 460 I = 1, N
            G(J) = G(J) + FVEC(I)*S2
            TH = FOUR*T2 + T*S2 - S1
            S1 = S2
            S2 = TH
            TH = T*T2 - T1
            T1 = T2
            T2 = TH
  460       CONTINUE
  470    CONTINUE
      D2 = TWO*D1
      DO 480 J = 1, N
         G(J) = D2*G(J)
  480    CONTINUE
  490 CONTINUE
      RETURN
C
C     LAST CARD OF SUBROUTINE GRDFCN.
C
      END


text:
SUBROUTINE SUB(N)
  CHARACTER(N) A(5)
  INQUIRE (IOLENGTH=IOL) A(1)  ! Inquire by output list
  OPEN (7,RECL=IOL)

END SUBROUTINE

max_stars_repo_name: maurizioabba/rose

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
!
!  Copyright 2017-2019 SALMON developers
!
!  Licensed under the Apache License, Version 2.0 (the "License");
!  you may not use this file except in compliance with the License.
!  You may obtain a copy of the License at
!
!      http://www.apache.org/licenses/LICENSE-2.0
!
!  Unless required by applicable law or agreed to in writing, software
!  distributed under the License is distributed on an "AS IS" BASIS,
!  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
!  See the License for the specific language governing permissions and
!  limitations under the License.
!
!--------10--------20--------30--------40--------50--------60--------70--------80--------90--------100-------110-------120-------130
subroutine Ion_Force_omp(Rion_update,GS_RT,ixy_m)
  use Global_Variables, only: zu_t,zu_m,zu_GS,NB,NBoccmax,calc_mode_gs,calc_mode_rt,Rion,Rion_m,force,force_m
  implicit none
  integer,intent(in) :: GS_RT
  logical,intent(in) :: Rion_update
  integer,intent(in),optional :: ixy_m

  select case(GS_RT)
    case(calc_mode_gs)
      call impl(Rion_update,zu_GS,NB)
    case(calc_mode_rt)
      if (present(ixy_m)) then
        Rion(:,:) = Rion_m(:,:,ixy_m)
        call impl(Rion_update,zu_m(:,:,:,ixy_m),NBoccmax)
        force_m(:,:,ixy_m) = force(:,:)
      else
        call impl(Rion_update,zu_t,NBoccmax)
      end if
    case default
      call err_finalize('ion_force_omp: gs_rt flag')
  end select

contains
  subroutine impl(Rion_update,zutmp,zu_NB)
    use Global_Variables
    use salmon_parallel, only: nproc_group_tdks,get_thread_id
    use salmon_communication, only: comm_summation, comm_is_root
    use timer
    use salmon_math
    use projector
    use opt_variables, only: NUMBER_THREADS_POW2
    implicit none
    logical,intent(in)       :: Rion_update
    integer,intent(in)       :: zu_NB
    complex(8),intent(inout) :: zutmp(NL,zu_NB,NK_s:NK_e)

    integer      :: ia,ib,ilma,ik,ix,iy,iz,n,j,i
    integer      :: ix_s,ix_e,iy_s,iy_e,iz_s,iz_e
    integer,allocatable:: idx(:),idy(:),idz(:)
    real(8)      :: rab(3),rab2,Gvec(3),G2,Gd
    real(8)      :: ftmp_l(3,NI)
    real(8)      :: FionR(3,NI), FionG(3,NI),nabt_wrk(4,3)
    complex(8)   :: uVpsi,duVpsi(3)
    complex(8)   :: dzudrzu(3),dzuekrdr(3)

    complex(8), allocatable :: dzudr(:,:,:,:)

    integer :: tid
    real(8) :: ftmp_t(3,NI,0:NUMBER_THREADS_POW2-1)

    allocate(dzudr(3,NL,zu_NB,NK_s:NK_e))
    ftmp_t(:,:,:) = 0.d0

    !flag_use_grad_wf_on_force is given in Gloval_Variable
    !if .true.   use gradient of wave-function (better accuracy)
    !if .false., old way which use gradient of potential is used (less accurate)

    call timer_begin(LOG_ION_FORCE)

    !ion-ion interaction (point charge Ewald)
    if (Rion_update) then

      !ion-ion: Ewald short interaction term (real space)
!$omp parallel private(tid)
      tid = get_thread_id()
      ftmp_t(:,:,tid)=0.d0
!$omp do private(ia, ik,ix,iy,iz,ib,rab,rab2) collapse(5)
      do ia=1,NI
      do ix=-NEwald,NEwald
      do iy=-NEwald,NEwald
      do iz=-NEwald,NEwald
      do ib=1,NI
        ik=Kion(ia)
        if(ix**2+iy**2+iz**2 == 0 .and. ia == ib) cycle
        rab(1)=Rion(1,ia)-ix*aLx-Rion(1,ib)
        rab(2)=Rion(2,ia)-iy*aLy-Rion(2,ib)
        rab(3)=Rion(3,ia)-iz*aLz-Rion(3,ib)
        rab2=sum(rab(:)**2)
        ftmp_t(:,ia,tid)=ftmp_t(:,ia,tid)&
             &-Zps(Kion(ia))*Zps(Kion(ib))*rab(:)/sqrt(rab2)*(-erfc_salmon(sqrt(aEwald*rab2))/rab2&
             &-2*sqrt(aEwald/(rab2*Pi))*exp(-aEwald*rab2))
      end do
      end do
      end do
      end do
      end do
!$omp end do
      call mreduce_omp(ftmp_l,ftmp_t,3*NI,tid)
!$omp end parallel
      FionR = ftmp_l

      !ion-ion: Ewald long interaction term (wave-number space)
!$omp parallel private(ia,tid)
      tid = get_thread_id()
      ftmp_t(:,:,tid)=0.d0
      do ia=1,NI
!$omp do private(ik,n,Gvec,G2,Gd)
      do n=NG_s,NG_e
        if(n == nGzero) cycle
        ik=Kion(ia)
        Gvec(1)=Gx(n); Gvec(2)=Gy(n); Gvec(3)=Gz(n)
        G2=sum(Gvec(:)**2)
        Gd=sum(Gvec(:)*Rion(:,ia))
        ftmp_t(:,ia,tid) = ftmp_t(:,ia,tid) &
        &   + Gvec(:)*(4*Pi/G2)*exp(-G2/(4*aEwald))*Zps(ik) &
        &     *zI*0.5d0*(conjg(rhoion_G(n))*exp(-zI*Gd)-rhoion_G(n)*exp(zI*Gd))
      end do
!$omp end do
      end do
      call mreduce_omp(ftmp_l,ftmp_t,3*NI,tid)
!$omp end parallel

      call comm_summation(ftmp_l,FionG,3*NI,nproc_group_tdks)
      Fion = FionR + FionG

    end if

    call update_projector_omp(kac)

    ! ion-electron
    if(flag_use_grad_wf_on_force)then

    ! Use gradient of wave-func for calculating force on ions

    !(prepare gradient of w.f.)
    ix_s = 0  ;  ix_e = NLx-1
    iy_s = 0  ;  iy_e = NLy-1
    iz_s = 0  ;  iz_e = NLz-1

    allocate(idx(ix_s-4:ix_e+4),idy(iy_s-4:iy_e+4),idz(iz_s-4:iz_e+4))
    do i=ix_s-4,ix_e+4
      idx(i) = mod(NLx+i,NLx)
    end do
    do i=iy_s-4,iy_e+4
      idy(i) = mod(NLy+i,NLy)
    end do
    do i=iz_s-4,iz_e+4
      idz(i) = mod(NLz+i,NLz)
    end do

    nabt_wrk(1:4,1) = nabx(1:4)
    nabt_wrk(1:4,2) = naby(1:4)
    nabt_wrk(1:4,3) = nabz(1:4)

!$omp parallel do collapse(2) default(none) &
!$omp          private(ik,ib) &
!$omp          shared(zutmp,dzudr,idx,idy,idz,nabt_wrk) &
!$omp          firstprivate(ix_s,ix_e,iy_s,iy_e,iz_s,iz_e,nk_s,nk_e,nboccmax)
    do ik=NK_s,NK_e
    do ib=1,NBoccmax
       call stencil_C_zu(zutmp(:,ib,ik) &
       &                ,dzudr(:,:,ib,ik) &
       &                ,ix_s,ix_e,iy_s,iy_e,iz_s,iz_e &
       &                ,idx,idy,idz,nabt_wrk)
    enddo
    enddo
!$omp end parallel do

    !Force from Vlocal with wave-func gradient --
!$omp parallel private(ia,tid)
    tid = get_thread_id()
    ftmp_t(:,:,tid) = 0.d0
!$omp do private(ik,ib,i,ia,dzudrzu) collapse(3)
    do ik=NK_s,NK_e
    do ib=1,NBoccmax
    do i=1,NL

       dzudrzu(:)=conjg(dzudr(:,i,ib,ik))*zutmp(i,ib,ik)

       do ia=1,NI
          ftmp_t(:,ia,tid) = ftmp_t(:,ia,tid) &
          &  -2d0* dble(dzudrzu(:)*Vpsl_ia(i,ia))*occ(ib,ik)*Hxyz
       enddo

    enddo
    enddo
    enddo
!$omp end do
    call mreduce_omp(ftmp_l,ftmp_t,3*NI,tid)
!$omp end parallel

    call timer_begin(LOG_ALLREDUCE)
    call comm_summation(ftmp_l,Floc,3*NI,nproc_group_tdks)
    call timer_end(LOG_ALLREDUCE)

    !Non-Local pseudopotential term using gradient of w.f.
!$omp parallel private(ia,tid)
    tid = get_thread_id()
    ftmp_t(:,:,tid) = 0.d0
!$omp do private(ik,j,i,ib,ilma,uVpsi,duVpsi,dzuekrdr) collapse(2)
    do ik=NK_s,NK_e
    do ib=1,NBoccmax
       do ilma=1,Nlma
          ia=a_tbl(ilma)
           uVpsi   =0.d0
          duVpsi(:)=0.d0
          do j=1,Mps(ia)
             i=Jxyz(j,ia)
            uVpsi      =uVpsi + uV(j,ilma)*ekr_omp(j,ia,ik)*zutmp(i,ib,ik)
            dzuekrdr(:)=(dzudr(:,i,ib,ik)+zI*kAc(ik,:)*zutmp(i,ib,ik))*ekr_omp(j,ia,ik)
            duVpsi(:)  =duVpsi(:) + conjg(dzuekrdr(:))*uV(j,ilma)
          end do
          uVpsi    =uVpsi    *Hxyz
          duVpsi(:)=duVpsi(:)*Hxyz
          ftmp_t(:,ia,tid) = ftmp_t(:,ia,tid) &
          &                  -2d0* dble(uVpsi*duVpsi(:))*iuV(ilma)*occ(ib,ik)
       end do
    end do
    end do
!$omp end do
    call mreduce_omp(ftmp_l,ftmp_t,3*NI,tid)
!$omp end parallel

    call timer_begin(LOG_ALLREDUCE)
    call comm_summation(ftmp_l,Fnl,3*NI,nproc_group_tdks)
    call timer_end(LOG_ALLREDUCE)

    else

    ! Use gradient of potential for calculating force on ions
    ! (older method, less accurate for larger grid size)

!$omp parallel private(ia,tid)
    tid = get_thread_id()
    ftmp_t(:,:,tid) = 0.d0
    do ia=1,NI
!$omp do private(ik,n,Gvec,G2,Gd)
    do n=NG_s,NG_e
      if(n == nGzero) cycle
      ik=Kion(ia)
      Gvec(1)=Gx(n); Gvec(2)=Gy(n); Gvec(3)=Gz(n)
      G2=sum(Gvec(:)**2)
      Gd=sum(Gvec(:)*Rion(:,ia))
      ftmp_t(:,ia,tid) = ftmp_t(:,ia,tid) &
      &    + zI*Gvec(:)*(4*Pi/G2)*Zps(ik)*exp(zI*Gd)*rhoe_G(n) &
      &    + conjg(rhoe_G(n))*dVloc_G(n,ik)*zI*Gvec(:)*exp(-zI*Gd)
    end do
!$omp end do
    end do
    call mreduce_omp(ftmp_l,ftmp_t,3*NI,tid)
!$omp end parallel

    call comm_summation(ftmp_l,Floc,3*NI,nproc_group_tdks)

    !non-local pseudopotential term
!$omp parallel private(ia,tid)
    tid = get_thread_id()
    ftmp_t(:,:,tid) = 0.d0
!$omp do private(ik,j,i,ib,ilma,uVpsi,duVpsi) collapse(2)
    do ik=NK_s,NK_e
    do ib=1,NBoccmax
      do ilma=1,Nlma
        ia=a_tbl(ilma)
        uVpsi=0.d0
        duVpsi(:)=0.d0
        do j=1,Mps(ia)
          i=Jxyz(j,ia)
          uVpsi=uVpsi+uV(j,ilma)*ekr_omp(j,ia,ik)*zutmp(i,ib,ik)
          duVpsi(:)=duVpsi(:)+duV(j,ilma,:)*ekr_omp(j,ia,ik)*zutmp(i,ib,ik)
        end do
        uVpsi=uVpsi*Hxyz; duVpsi(:)=duVpsi(:)*Hxyz
        ftmp_t(:,ia,tid)=ftmp_t(:,ia,tid)+(conjg(uVpsi)*duVpsi(:)+uVpsi*conjg(duVpsi(:)))*iuV(ilma)*occ(ib,ik)
      end do
    end do
    end do
!$omp end do
    call mreduce_omp(ftmp_l,ftmp_t,3*NI,tid)
!$omp end parallel

    call timer_begin(LOG_ALLREDUCE)
    call comm_summation(ftmp_l,fnl,3*NI,nproc_group_tdks)
    call timer_end(LOG_ALLREDUCE)

    endif ! flag_use_grad_wf_on_force

    call timer_begin(LOG_ALLREDUCE)
    force = Fion + Floc + Fnl
    call timer_end(LOG_ALLREDUCE)

    call timer_end(LOG_ION_FORCE)
  end subroutine

!OpenMP manual reduction routine
  subroutine mreduce_omp(vout,vtmp,vsize,tid)
    use opt_variables, only: NUMBER_THREADS_POW2
    implicit none
    integer,intent(in)  :: vsize,tid
    real(8),intent(out) :: vout(vsize)
    real(8)             :: vtmp(vsize,0:NUMBER_THREADS_POW2-1)

    integer :: i

    i = NUMBER_THREADS_POW2/2
    do while(i > 0)
      if(tid < i) then
        vtmp(:,tid) = vtmp(:,tid) + vtmp(:,tid + i)
      end if
      i = i/2
!$omp barrier
    end do

    if (tid == 0) then
      vout(:) = vtmp(:,0)
    end if
  end subroutine

!Gradient of wave function (du/dr) with nine points formura
# define DX(dt) iz,iy,idx(ix+(dt))
# define DY(dt) iz,idy(iy+(dt)),ix
# define DZ(dt) idz(iz+(dt)),iy,ix
  subroutine stencil_C_zu(zu0,dzu0dr &
  &                      ,ix_s,ix_e,iy_s,iy_e,iz_s,iz_e &
  &                      ,idx,idy,idz,nabt)
  implicit none
  integer   ,intent(in)  :: ix_s,ix_e,iy_s,iy_e,iz_s,iz_e
  integer   ,intent(in)  :: idx(ix_s-4:ix_e+4),idy(iy_s-4:iy_e+4),idz(iz_s-4:iz_e+4)
  real(8)   ,intent(in)  :: nabt(4,3)
  complex(8),intent(in)  :: zu0(iz_s:iz_e,iy_s:iy_e,ix_s:ix_e)
  complex(8),intent(out) :: dzu0dr(3,iz_s:iz_e,iy_s:iy_e,ix_s:ix_e)
  !
  integer :: iz,iy,ix
  complex(8) :: w(3)

  do ix=ix_s,ix_e
  do iy=iy_s,iy_e
!dir$ vector nontemporal(dzu0dr)
  do iz=iz_s,iz_e

    w(1) =  nabt(1,1)*(zu0(DX(1)) - zu0(DX(-1))) &
           +nabt(2,1)*(zu0(DX(2)) - zu0(DX(-2))) &
           +nabt(3,1)*(zu0(DX(3)) - zu0(DX(-3))) &
           +nabt(4,1)*(zu0(DX(4)) - zu0(DX(-4)))

    w(2) =  nabt(1,2)*(zu0(DY(1)) - zu0(DY(-1))) &
           +nabt(2,2)*(zu0(DY(2)) - zu0(DY(-2))) &
           +nabt(3,2)*(zu0(DY(3)) - zu0(DY(-3))) &
           +nabt(4,2)*(zu0(DY(4)) - zu0(DY(-4)))

    w(3) =  nabt(1,3)*(zu0(DZ(1)) - zu0(DZ(-1))) &
           +nabt(2,3)*(zu0(DZ(2)) - zu0(DZ(-2))) &
           +nabt(3,3)*(zu0(DZ(3)) - zu0(DZ(-3))) &
           +nabt(4,3)*(zu0(DZ(4)) - zu0(DZ(-4)))

    dzu0dr(:,iz,iy,ix) = w(:)
  end do
  end do
  end do

  return
  end subroutine stencil_C_zu

end subroutine Ion_Force_omp
!--------10--------20--------30--------40--------50--------60--------70--------80--------90--------100-------110-------120-------130

max_stars_repo_name: shunsuke-sato/SALMON

max_issues_repo_name: shunsuke-sato/SALMON

max_forks_repo_name: shunsuke-sato/SALMON


text:
c     Defines the functions h1(z,v) and h2(z,v) according to
c     the definition in eq.(14) of hep-ph/0002078

      function h1(z)
      implicit none
      complex *16 h1,z,zz
      real *8 t,v,acc,tmp1,tmp2
      real *8 err,aa
      integer i,j
      common/ij/i,j
      common/v/v
      common/c/zz
      external aa
      include 'const.F'
      zz=z
      acc=1d-3
c     Real part
      h1=(0d0,0d0)
      j=1
      do i=1,3
      call dadapt(aa,0d0,1d0,5,1d-15,1d-8,tmp1,err)
      h1=h1+tmp1
      enddo
c     Imaginary part
      j=2
      do i=1,3
      call dadapt(aa,0d0,1d0,5,1d-15,1d-8,tmp2,err)
      h1=h1+II*tmp2
      enddo
      h1=-h1/pi
      return
      end

      function h2(z)
      implicit none
      complex *16 h2,z,zz
      common/c/zz
      real *8 t,v,acc,tmp1,tmp2
      real *8 err,aa
      integer i,j
      common/ij/i,j
      common/v/v
      include 'const.F'
      external aa
      zz=z
      acc=1d-3
c     Real part
      h2=(0d0,0d0)
      j=1
      do i=4,6
      call dadapt(aa,0d0,1d0,5,1d-15,1d-8,tmp1,err)
      h2=h2+tmp1
      enddo
c     Imaginary part
      j=2
      do i=4,6
      call dadapt(aa,0d0,1d0,5,1d-15,1d-8,tmp2,err)
      h2=h2+II*tmp2
      enddo
      h2=-h2/pi
      return
      end

c     The function aa(t) give the parametrization of the real (j=1)
c     or imaginary (j=2) part of the function Exp[-i z sin(teta)]
c     and includes the jacobian
c     i=1,2,3 needed to compute h1
c     i=4,5,6 needed to compute h2

      function aa(t)
      implicit none
      complex *16 zz,a,teta,jac
      real *8 v,aa,t
      integer i,j
      common/c/zz
      common/ij/i,j
      common/v/v
      include 'const.F'
      if(i.eq.1) then
      teta=-II*pi*v*t
      jac=II*pi*v
      elseif(i.eq.2) then
      teta=-pi*t
      jac=-pi
      elseif(i.eq.3) then
      teta=-Pi+II*v*pi*t
      jac=II*v*pi
      elseif(i.eq.4) then
      teta=pi+II*pi*v*t
      jac=-II*pi*v
      elseif(i.eq.5) then
      teta=pi*t
      jac=-pi
      elseif(i.eq.6) then
      teta=-II*v*pi*t
      jac=-II*v*pi
      endif
      a=jac*exp(-II*zz*sin(teta))
      if(j.eq.1) then
      aa=dreal(a)
      elseif(j.eq.2) then
      aa=dimag(a)
      endif
      return
      end

max_stars_repo_name: gvita/RHEHpt

max_issues_repo_name: gvita/RHEHpt

max_forks_repo_name: gvita/RHEHpt


             program demo_vinit
             use M_draw
             use ISO_C_BINDING
             integer :: ios
             character(len=50) :: device

             ! read in device name and start graphics mode
             print*,'Enter output device:'
             read(*,'(a)',iostat=ios)device
             if(ios.ne.0)device=' '
             call vinit(device)
             ! by default the drawing surface is
             ! a square ranging from -1 to 1 in both
             ! the X and Y axis

             ! set font to large hardware font
             call font('large')

             ! set current color to black
             call color(D_BLACK)

             ! clear to current color
             call clear()

             ! we want to draw in green
             call color(D_GREEN)

             ! draw a horizontal line at y = 0
             call move2(-1.0, 0.0)
             call draw2(1.0, 0.0)

             ! pause for some input
             idum=getkey()

             ! draw a line along x = 0
             call move2(0.0, 0.0)
             call draw2(0.0, 1.0)

             ! move to the middle of the screen
             call move2(0.0, 0.0)

             ! draw 'Hello' starting at the origin
             call drawstr('Hello')

             ! pause again
             idum=getkey()

             !  wrap up and exit graphics mode
             call vexit()

             end program demo_vinit


text:
Interstate 505 connects Interstate 80 in north Vacaville with Interstate 5 in northern Yolo County. I505 goes mainly through rural areas, although it goes through the eastern edge of Winters. It is two lanes each way and is mainly used by truckers and other vehicles transporting goods. I505 was built before Highway 113 was a real highway and was meant to be a short cut through Yolo County into Solano County so trucks would not have to go all the way into Sacramento to reach I80.

From CR14 to I5 it is bike accessible.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
module neural

  ! Explicit typing only
  implicit none

  ! Everything is private unless specified
  private

  public :: neural_init_network
  public :: neural_use_network
  public :: neural_train_network_file
  public :: neural_train_network_direct
  public :: neural_load_network
  public :: neural_load_network_from_string
  public :: neural_save_network
  public :: neural_save_network_to_string
  public :: neural_convert_input
  public :: neural_convert_output
  public :: neural_sigmoid

  public :: network

  ! The neural network type, stores info about lauer sizes, weights and biases
  type network
      integer                                         :: nl            ! Number of layers
      integer                                         :: max_size
      integer, dimension(:), allocatable              :: ls            ! Layer sizes
      double precision, dimension(:, :, :), allocatable :: weights
      double precision, dimension(:, :), allocatable    :: biases
  end type

  interface neural_sigmoid
    module procedure neural_sigmoid_vector
    module procedure neural_sigmoid_array
  end interface

contains

  ! Train a certain network from a file, repeating a number of times
  subroutine neural_train_network_file(net, file_name, repeat, learn_rate)

    ! The network to train
    type(network), intent(inout) :: net

    ! The file to train on
    character(*), intent(in) :: file_name

    ! How many times to repeat the training
    integer, optional, intent(in) :: repeat

    ! How fast should it learn? [0,1]
    double precision, optional, intent(in) :: learn_rate

    ! Internal vars
    double precision, dimension(:, :), allocatable :: input, output
    integer :: file_num
    integer :: ierr, input_size, output_size, num_lines, i, j, num_data

    open(newunit=file_num, file=file_name, status="old", iostat=ierr)
    if (ierr /= 0) stop "Couldn't load file for training"

    ! Count the number of lines in the file
    num_lines = 0
    do

      read(file_num, *, iostat=ierr)
      if (ierr /= 0) exit
      num_lines = num_lines + 1

    end do

    ! Restart
    rewind(file_num)

    ! Get the size of the input and output arrays
    read(file_num, *, iostat=ierr) input_size
    read(file_num, *, iostat=ierr) output_size

    ! Check if these are right for the network and exit otherwise
    if (input_size /= net%ls(1) .or. output_size /= net%ls(net%nl)) then
      print *, "Error training from file, wrong number of inputs/outputs"
      return
    end if

    ! Determine how many data sets there are
    num_data = int((num_lines - 2) / (input_size + output_size))

    ! Allocate the arrays for storing the input/output data
    allocate(input(num_data, input_size), stat=ierr)
    if (ierr /= 0) stop "Error allocating input array in neural_train_network_file"
    allocate(output(num_data, output_size), stat=ierr)
    if (ierr /= 0) stop "Error allocating output array in neural_train_network_file"

    ! Load the data into the input/output arrays
    do i = 1, num_data

      do j = 1, input_size
        read(file_num, *, iostat=ierr) input(i, j)
      end do
      do j = 1, output_size
        read(file_num, *, iostat=ierr) output(i, j)
      end do

    end do

    close(file_num)

    ! Train the network using this data now in memory
    call neural_train_network_direct(net, input, output, repeat, learn_rate)

  end subroutine neural_train_network_file

  ! Train a network directly on some input/output data
  subroutine neural_train_network_direct(net, input, output, repeat, learn_rate)

    type(network), intent(inout) :: net
    double precision, dimension(:, :), intent(in) :: input
    double precision, dimension(:, :), intent(in) :: output
    integer, optional, intent(in) :: repeat
    double precision, optional, intent(in) :: learn_rate

    integer :: i,j,k,l, max, ierr
    double precision :: lr

    double precision, dimension(net%nl, net%max_size, net%max_size) :: delta_weights
    double precision, dimension(net%nl-1, net%max_size) :: delta_biases
    double precision, dimension(net%nl, size(input, 1), net%max_size) :: layer_activations, layer_outputs
    double precision, dimension(net%nl, size(input, 1), net%max_size) :: layer_errors, d_layer

    ! Check if these sizes are right for the network and exit otherwise
    if (size(input, 2) /= net%ls(1) .or. size(output, 2) /= net%ls(net%nl)) then
      stop 'Error in neural_train_network_direct: wrong number of inputs/outputs'
      return
    end if

    ! If given a learn_rate, use that instead of the default
    if (present(learn_rate)) then
      lr = learn_rate
    else
      lr = 1.0
    end if

    ! If given a repeat number, use that instead of the default
    if (present(repeat)) then
      max = repeat
    else
      max = 1
    end if

    ! Repeat as many times as requested
    do l=1, max

      ! Reset things just in case
      layer_activations = 0
      layer_outputs = 0
      layer_errors = 0
      d_layer = 0

      ! Treat the input like the output of the first layer
      layer_outputs(1, :, 1:net%ls(1)) = input

      ! Loop over the layers
      do j=2, net%nl

        ! Multiply the previous layer results by the weights
        layer_activations(j, :, 1:net%ls(j)) = matmul(layer_outputs(j-1,:,1:net%ls(j-1)), &
        & net%weights(j-1,1:net%ls(j-1),1:net%ls(j)))

        ! Add the biases
        do i=1, size(input, 1)
          layer_activations(j, i, 1:net%ls(j)) = layer_activations(j, i, 1:net%ls(j)) + net%biases(j-1, 1:net%ls(j))
        end do

        ! Sigmoid the activations to get the outputs
        layer_outputs(j, :, 1:net%ls(j)) = neural_sigmoid_array(layer_activations(j, :, 1:net%ls(j)))

      end do

      ! Do output layer seperate since slightly different
      layer_errors(net%nl, :, 1:net%ls(net%nl)) = output - layer_outputs(net%nl, :, 1:net%ls(net%nl))
      d_layer(net%nl, :, 1:net%ls(net%nl)) = layer_errors(net%nl, :, 1:net%ls(net%nl)) * &
      & neural_sigmoid_prime_array(layer_outputs(net%nl, :, 1:net%ls(net%nl)))

      ! Backpropagate the rest of the layers
      do j=net%nl-1, 2, -1

        layer_errors(j, :, 1:net%ls(j)) = matmul(d_layer(j+1, :, 1:net%ls(j+1)), &
        & transpose(net%weights(j, 1:net%ls(j), 1:net%ls(j+1))))
        d_layer(j, :, 1:net%ls(j)) = layer_errors(j, :, 1:net%ls(j)) * neural_sigmoid_prime_array(layer_outputs(j, :, 1:net%ls(j)))

      end do

      ! Determine the change in weights and biases
      do j=net%nl-1, 1, -1

        delta_weights(j,1:net%ls(j),1:net%ls(j+1)) = matmul(transpose(layer_outputs(j, :, 1:net%ls(j))), &
        & d_layer(j+1, :, 1:net%ls(j+1))) * lr
        delta_biases(j,1:net%ls(j+1)) = sum(d_layer(j+1, :, net%ls(j+1)), 1) * lr

      end do

      ! Update the network weights
      do i = 1, net%nl-1
        do j = 1, net%max_size
          do k = 1, net%max_size
             net%weights(i,j,k) = net%weights(i,j,k) + delta_weights(i,j,k)
          end do
        end do
      end do

      ! Update the network biases
      do i = 1, net%nl-1
        do j = 1, net%max_size
            net%biases(i,j) = net%biases(i,j) + delta_biases(i,j)
        end do
      end do

    end do

  end subroutine neural_train_network_direct

  ! Forward propogate some inputs through the network
  function neural_use_network(net, input)

    type(network), intent(in) :: net
    double precision, dimension(net%ls(1)), intent(in) :: input
    double precision, dimension(net%ls(net%nl)) :: neural_use_network
    double precision, dimension(net%nl, 1, net%max_size) :: layer_activations, layer_outputs
    integer :: j, ierr

    ! Treat the input like the output of the first layer
    layer_outputs(1, 1, 1:net%ls(1)) = input

    ! Loop over the layers
    do j=2, net%nl

      ! Multiply the previous layer results by the weights
      layer_activations(j, 1, 1:net%ls(j)) = matmul(layer_outputs(j-1,1,1:net%ls(j-1)), net%weights(j-1,1:net%ls(j-1),1:net%ls(j)))

      ! Add the biases
      layer_activations(j, 1, 1:net%ls(j)) = layer_activations(j, 1, 1:net%ls(j)) + net%biases(j-1, 1:net%ls(j))

      ! Sigmoid the activations to get the outputs
      layer_outputs(j, 1, 1:net%ls(j)) = neural_sigmoid_vector(layer_activations(j, 1, 1:net%ls(j)))

    end do

    ! The results are the last layer outputs
    neural_use_network = layer_outputs(net%nl, 1, 1:net%ls(net%nl))

  end function neural_use_network

  ! Initialise the network with random weights
  subroutine neural_init_network(net, sizes)

    type(network), intent(inout) :: net
    integer, dimension(:), intent(in) :: sizes
    integer :: i, j, k, ierr

    ! Dellocate if already initialized
    if (allocated(net%ls)) then
      deallocate(net%ls, stat=ierr)
      if (ierr /= 0) stop 'Error in deallocating net%ls in neural_init_network'
    end if
    if (allocated(net%weights)) then
      deallocate(net%weights, stat=ierr)
      if (ierr /= 0) stop 'Error in deallocating net%weights in neural_init_network'
    end if
    if (allocated(net%biases)) then
      deallocate(net%biases, stat=ierr)
      if (ierr /= 0) stop 'Error in deallocating net%biases in neural_init_network'
    end if

    net%nl = size(sizes(:))

    ! Allocate the layer sizes list
    allocate(net%ls(net%nl), stat=ierr)
    if (ierr /= 0) stop "failed to allocate array"

    net%ls = sizes
    net%max_size = maxval(sizes(:))

    ! Allocate the weights and biases arrays
    allocate(net%weights(net%nl-1, net%max_size, net%max_size), stat=ierr)
    if (ierr /= 0) stop "failed to allocate array"
    allocate(net%biases(net%nl-1, net%max_size), stat=ierr)
    if (ierr /= 0) stop "failed to allocate array"

    ! Init the random generator
    call RANDOM_seed()

    ! Randomise the initial weights
    net%weights = 0
    do i=1, net%nl-1
      do j=1, net%ls(i)
        do k=1, net%ls(i+1)
          net%weights(i,j,k) = neural_normal()
        end do
      end do
    end do

    ! Randomise the initial biases
    net%biases = 0
    do i=1, net%nl-1
      do j=1, net%max_size
          net%biases(i,j) = neural_normal()
      end do
    end do

  end subroutine neural_init_network

  ! Load the network from a file
  subroutine neural_load_network(net, file)

    type(network), intent(inout) :: net
    character(*), intent(in) :: file
    integer :: file_num
    integer :: i, j, k, ierr

    ! Deallocate if already allocated
    if (allocated(net%ls)) then
      deallocate(net%ls, stat=ierr)
      if (ierr /= 0) stop 'Error in deallocating net%ls in neural_load_network'
    end if
    if (allocated(net%weights)) then
      deallocate(net%weights, stat=ierr)
      if (ierr /= 0) stop 'Error in deallocating net%weights in neural_load_network'
    end if
    if (allocated(net%biases)) then
      deallocate(net%biases, stat=ierr)
      if (ierr /= 0) stop 'Error in deallocating net%biases in neural_load_network'
    end if

    ! Open the file containing network info
    open(unit=file_num, file=file, status="old")
    if (ierr /= 0) stop 'Failed to open file'

    ! Sizes first
    read(file_num, *) net%nl
    allocate(net%ls(net%nl), stat=ierr)
    if (ierr /= 0) stop "Failed to allocate array"

    do i = 1, net%nl
      read(file_num, *) net%ls(i)
    end do

    ! Recalculate the max size
    net%max_size = maxval(net%ls(:))

    ! Allocate the arrays
    allocate(net%weights(net%nl-1, net%max_size, net%max_size), stat=ierr)
    if (ierr /= 0) stop "failed to allocate array"
    allocate(net%biases(net%nl-1, net%max_size), stat=ierr)
    if (ierr /= 0) stop "failed to allocate array"

    ! Then weights
    do i = 1, net%nl-1
      do j = 1, net%ls(i)
        do k = 1, net%ls(i+1)
          read(file_num, *) net%weights(i,j,k)
        end do
      end do
    end do

    ! Then biases
    do i = 1, net%nl-1
      do j = 1, net%ls(i)
        read(file_num, *) net%biases(i,j)
      end do
    end do

    close(file_num)

  end subroutine neural_load_network

  ! Save the network to a file
  subroutine neural_save_network(net, file)

    type(network), intent(in) :: net
    character(*), intent(in) :: file
    integer :: file_num
    integer :: i, j, k, ierr

    ! Open the file to save the network to
    open(newunit=file_num, file=file, status="new")
    if (ierr /= 0) stop "Failed to open file"

    ! Sizes first
    write(file_num, *) net%nl
    do i = 1, net%nl
      write(file_num, *) net%ls(i)
    end do

    ! Then weights
    do i = 1, net%nl-1
      do j = 1, net%ls(i)
        do k = 1, net%ls(i+1)
          write(file_num, *) net%weights(i,j,k)
        end do
      end do
    end do

    ! Then biases
    do i = 1, net%nl-1
      do j = 1, net%ls(i)
        write(file_num, *) net%biases(i,j)
      end do
    end do

    close(file_num)

  end subroutine neural_save_network

  ! Save the network to a string and return it
  function neural_save_network_to_string(net)

    type(network), intent(in) :: net
    character(300) :: temp_result
    character(:), allocatable :: neural_save_network_to_string
    character(10) :: temp
    integer :: i, j, k, ierr

    ! Start with a blank string
    temp_result = ""
    temp = ""

    ! Sizes first
    write(temp, "(I3)") net%nl
    temp_result = trim(temp_result) // trim(temp) // ","
    do i = 1, net%nl
      write(temp, "(I3)") net%ls(i)
      temp_result = trim(temp_result) // trim(temp) // ","
    end do

    ! Then weights
    do i = 1, net%nl-1
      do j = 1, net%ls(i)
        do k = 1, net%ls(i+1)
          write(temp, "(f10.5)") net%weights(i,j,k)
          temp_result = trim(temp_result) // trim(temp) // ","
        end do
      end do
    end do

    ! Then biases
    do i = 1, net%nl-1
      do j = 1, net%ls(i)
        write(temp, "(f10.5)") net%biases(i,j)
        temp_result = trim(temp_result) // trim(temp) // ","
      end do
    end do

    ! Compress the result by removing whitespace
    neural_save_network_to_string = ""
    do i=1, len_trim(temp_result)
      if (temp_result(i:i) /= " ") then
        neural_save_network_to_string = trim(neural_save_network_to_string) // temp_result(i:i)
      end if
    end do

  end function neural_save_network_to_string

  ! Load the network from a string
  subroutine neural_load_network_from_string(net, string)

    type(network), intent(inout) :: net
    character(*), intent(in) :: string
    integer :: i, j, k, ierr, l, prev_l

    ! Deallocate if already allocated
    if (allocated(net%ls)) then
      deallocate(net%ls, stat=ierr)
      if (ierr /= 0) stop 'Error in deallocating net%ls in neural_load_network'
    end if
    if (allocated(net%weights)) then
      deallocate(net%weights, stat=ierr)
      if (ierr /= 0) stop 'Error in deallocating net%weights in neural_load_network'
    end if
    if (allocated(net%biases)) then
      deallocate(net%biases, stat=ierr)
      if (ierr /= 0) stop 'Error in deallocating net%biases in neural_load_network'
    end if

    ! Read the number of layers
    prev_l = 0
    l = index(string(:), ",")
    read(string(prev_l+1:l-1), *) net%nl
    allocate(net%ls(net%nl), stat=ierr)
    if (ierr /= 0) stop "failed to allocate array"

    ! Read the size of each layer
    do i = 1, net%nl
      prev_l = l
      l = l+index(string(l+1:), ",")
      read(string(prev_l+1:l-1), *) net%ls(i)
    end do

    ! Recalculate the max size
    net%max_size = maxval(net%ls(:))

    ! Allocate the arrays
    allocate(net%weights(net%nl-1, net%max_size, net%max_size), stat=ierr)
    if (ierr /= 0) stop "failed to allocate array"
    allocate(net%biases(net%nl-1, net%max_size), stat=ierr)
    if (ierr /= 0) stop "failed to allocate array"

    ! Then weights
    net%weights = 0
    do i = 1, net%nl-1
      do j = 1, net%ls(i)
        do k = 1, net%ls(i+1)
          prev_l = l
          l = l+index(string(l+1:), ",")
          read(string(prev_l+1:l-1), *) net%weights(i,j,k)
        end do
      end do
    end do

    ! Then biases
    net%biases = 0
    do i = 1, net%nl-1
      do j = 1, net%ls(i)
        prev_l = l
        l = l+index(string(l+1:), ",")
        read(string(prev_l+1:l-1), *) net%biases(i,j)
      end do
    end do

  end subroutine neural_load_network_from_string

  ! Apply the sigmoid to a 1D array
  function neural_sigmoid_vector(x)

    double precision, dimension(:), intent(in) :: x
    double precision, dimension(size(x)) :: neural_sigmoid_vector
    integer :: i

    do i=1, size(x)
      neural_sigmoid_vector(i) = 1.0 / (1.0 + exp(-x(i)))
    end do

  end function neural_sigmoid_vector

  ! Apply the sigmoid to a 2D array
  function neural_sigmoid_array(x)

    double precision, dimension(:,:), intent(in) :: x
    double precision, dimension(size(x,1), size(x,2)) :: neural_sigmoid_array
    integer :: i, j

    do i=1, size(x,1)
      do j=1, size(x,2)
        neural_sigmoid_array(i, j) = 1.0 / (1.0 + exp(-x(i, j)))
      end do
    end do

  end function neural_sigmoid_array

  ! Apply the inverse of the sigmoid to a 1D array
  function neural_sigmoid_inv_vector(x)

    double precision, dimension(:), intent(in) :: x
    double precision, dimension(size(x)) :: neural_sigmoid_inv_vector
    integer :: i

    do i=1, size(x)
      neural_sigmoid_inv_vector(i) = -log((1.0 / x(i)) - 1.0)
    end do

  end function neural_sigmoid_inv_vector

  ! Apply the derivative of the sigmoid to a 2D array
  function neural_sigmoid_prime_array(x)

    double precision, dimension(:,:), intent(in) :: x
    double precision, dimension(size(x, 1), size(x, 2)) :: neural_sigmoid_prime_array

    neural_sigmoid_prime_array(:, :) = x(:, :)*(1.0 - x(:, :))

  end function neural_sigmoid_prime_array

  ! Return a number in a normal (ish) distribution
  function neural_normal()

    double precision :: neural_normal
    double precision :: temp_rand1, temp_rand2, temp_rand3

    call random_number(temp_rand1)
    call random_number(temp_rand2)
    call random_number(temp_rand3)

    ! Simple approximation, tends to stick around 0.5 as desired
    neural_normal = (temp_rand1 + temp_rand2 + temp_rand3) / 3.0

  end function neural_normal

end module neural

max_stars_repo_name: Lumorti/fortran-ann

max_issues_repo_name: Lumorti/fortran-ann

max_forks_repo_name: Lumorti/fortran-ann