      module channel_module

      implicit none

      integer :: jhyd    !units         |description
      integer :: jsed    !units         |description
      integer :: jnut    !units         |description
      real :: rttime     !hr            |reach travel time
      real :: ben_area   !m2            |benthic area (bottom sediments)
      real :: rchdep     !m             |depth of flow on day
      real :: rtevp      !m^3 H2O       |evaporation from reach on day
      real :: rttlc      !m^3 H2O       |transmission losses from reach on day
      real :: pet_ch     ! mm           |potential evaporation from reach on day
      real, dimension (:), allocatable :: hrtwtr     !m^3 H2O       |water leaving reach
      real, dimension (:), allocatable :: hharea     !m^2           |cross-sectional area of flow
      real, dimension (:), allocatable :: hdepth     !m             |depth of flow
      real, dimension (:), allocatable :: rhy        !m H2O         |main channel hydraulic radius
      real, dimension (:), allocatable :: hsdti      !m^3/s         |flow rate in reach for hour
      real, dimension (:), allocatable :: hhtime     !hr            |flow travel time for hour
      real, dimension (:), allocatable :: hrttlc     !m^3 H2O       |transmission losses from reach during time step
      real, dimension (:), allocatable :: hrtevp     !m^3 H2O       |evaporation from reach during time step
      real, dimension (:), allocatable :: hhstor     !m^3 H2O       |water stored in reach at end of hour
      real, dimension (:), allocatable :: hrchwtr    !m^3 H2O       |water stored at beginning of day
      real, dimension (:), allocatable :: halgae     !mg alg/L      |algal biomass concentration in reach
      real, dimension (:), allocatable :: hbactlp    !# cfu/100mL   |less persistent bacteria in reach/outflow during hour
      real, dimension (:), allocatable :: hbactp     !# cfu/100mL   |persistent bacteria in reach/outflow during hour
      real, dimension (:), allocatable :: hbod       !mg O2/L       |carbonaceous biochemical oxygen demand inreach at end of hour
      real, dimension (:), allocatable :: hchla      !mg chl-a/L    |chlorophyll-a concentration in reach at end of hour
      real, dimension (:), allocatable :: hdisox     !mg O2/L       |dissolved oxygen concentration in reach at end of hour
      real, dimension (:), allocatable :: hnh4       !mg N/L        |ammonia concentration in reach at end of hour
      real, dimension (:), allocatable :: hno2       !mg N/L        |nitrite concentration in reach at end of hour
      real, dimension (:), allocatable :: hno3       !mg N/L        |nitrate concentration in reach at end of hour
      real, dimension (:), allocatable :: horgn      !mg N/L        |organic nitrogen concentration in reach at end of hour
      real, dimension (:), allocatable :: horgp      !mg P/L        |organic phosphorus concentration in reach at end of hour
      real, dimension (:), allocatable :: hsedst     !metric tons   |amount of sediment stored in reach at the end of hour
      real, dimension (:), allocatable :: hsedyld    !metric tons   |sediment transported out of reach during hour
      real, dimension (:), allocatable :: hsolp      !mg P/L        |dissolved phosphorus concentration in reach at end of hour
      real, dimension (:), allocatable :: hsolpst    !mg pst/m^3    |soluble pesticide concentration in outflow on day
      real, dimension (:), allocatable :: hsorpst    !mg pst/m^3    |sorbed pesticide concentration in outflow on day
      real, dimension (:), allocatable :: rchsep     !

      real :: peakr, rcharea, sdti
      real :: bnkrte              !
      real :: degrte              !
      real :: sedrch              !metric tons       |sediment transported out of reach on day
      real :: rch_san             !
      real :: rch_sil             !
      real :: rch_cla             !
      real :: rch_sag
      real :: rtwtr_d             !m^3 H2O           |water leaving reach during day
      real :: rt_delt             ! calculation time step in days
      real :: rch_lag             !
      real :: rch_gra             !
      real :: rtwtr               !m^3 H2O           |water leaving reach on day
      real :: wtrin               !m^3               |water entering reach during day
      integer:: sed_ch

      type channel
          real :: algae = 0.     ! mg alg/L      |algal biomass concentration in reach
          real :: ammonian = 0.  ! mg N/L        |ammonia concentration in reach
          real :: bankst = 0.    ! m^3 H2O       |bank storage
          real :: li = 0.        ! km            |initial length of main channel
          real :: orgn = 0.      !               |organic nitrogen contribution from channel erosion
          real :: orgp = 0.      !               |organic phosphorus contribution from channel erosion
          real :: si = 0.        !(m/n)          |slope of main channel
          real :: wi = 0.        !(m)            |width of main channel at top of bank
          real :: di = 0.        !(m)            |depth of main channel from top of bank to bottom
          real :: chlora = 0.    ! mg chl-a/L    |chlorophyll-a concentration in reach
          real :: pst_conc =0.   ! mg/(m**3)     |initial pesticide concentration in reach
          real :: dep_chan =0.   ! m             |average daily water depth in channel
          real :: disolvp = 0.   ! mg P/L        |dissolved P concentration in reach
          real :: drift = 0.     ! kg            |amount of pesticide drifting onto main channel in subbasin
          real :: flwin = 0.     ! m^3 H2O       |flow into reach on previous day
          real :: flwout = 0.    ! m^3 H2O       |flow out of reach on previous day
          real :: nitraten = 0.  ! mg N/L        |nitrate concentration in reach
          real :: nitriten = 0.  ! mg N/L        |nitrite concentration in reach
          real :: organicn = 0.  ! mg N/L        |organic nitrogen concentration in reach
          real :: organicp = 0.  ! mg P/L        |organic phosphorus concentration in reach
          real :: rch_bactlp= 0. ! # cfu/100ml   |less persistent bacteria stored in reach
          real :: rch_bactp = 0. ! # cfu/100ml   |persistent bacteria stored in reach
          real :: rch_cbod = 0.  ! mg O2/L       |carbonaceous biochemical oxygen demand in reach
          real :: rch_dox = 0.   ! mg O2/L       |dissolved oxygen concentration in reach
          real :: rchstor = 0.   ! m^3 H2O       |water stored in reach
          real :: sedst = 0.     ! metric tons   |amount of sediment stored in reach
          real :: vel_chan = 0.  ! m/s           |average flow velocity in channel
          real :: bed_san = 0.
          real :: bed_sil = 0.
          real :: bed_cla = 0.
          real :: bed_gra = 0.
          real :: bnk_san = 0.
          real :: bnk_sil = 0.
          real :: bnk_cla = 0.
          real :: bnk_gra = 0.
          real :: depfp = 0.
          real :: depprfp = 0.
          real :: depsilfp = 0.
          real :: depclafp = 0.
          real :: depch = 0.
          real :: depprch = 0.
          real :: depsanch = 0.
          real :: depsilch = 0.
          real :: depclach = 0.
          real :: depsagch= 0.
          real :: deplagch = 0.
          real :: depgrach = 0.
          real :: sanst = 0.
          real :: silst = 0.
          real :: clast = 0.
          real :: sagst = 0.
          real :: lagst = 0.
          real :: grast = 0.
          real :: wattemp = 0.
          real :: bactp = 0.
          real :: chfloodvol = 0.
          real :: bactlp = 0.
      end type channel
      type (channel), dimension(:), allocatable :: ch

      type ch_output
          real :: flo_in = 0.                  ! (ha-m)     |streamflow into reach during time step
          real :: flo_out = 0.                 ! (ha-m)     |streamflow out of reach during time step
          real :: evap = 0.                    ! (m^3/s)    |daily rate of water loss from reach by evaporation
          real :: tloss = 0.                   ! (m^3/s)    |rate of water loss from reach by transmission through the streambed
          real :: sed_in = 0.                  ! (tons)     |sediment transported with water into reach
          real :: sed_out = 0.                 ! (tons)     |sediment transported with water out of reach
          real :: sed_conc = 0.                ! (mg/L)     |concentration of sediment in reach
          real :: orgn_in = 0.                 ! (kg N)     |organic nitrogen transported with water into reach
          real :: orgn_out = 0.                ! (kg N)     |organic nitrogen transported with water out of reach
          real :: orgp_in = 0.                 ! (kg P)     |organic phosphorus transported with water into reach
          real :: orgp_out = 0.                ! (kg P)     |organic phosphorus transported with water out of reach
          real :: no3_in = 0.                  ! (kg N)     |nitrate transported with water into reach
          real :: no3_out = 0.                 ! (kg N)     |nitrate transported with water out of reach
          real :: nh4_in = 0.                  ! (kg)       |ammonium transported with water into reach
          real :: nh4_out = 0.                 ! (kg)       |ammonium transported with water out of reach
          real :: no2_in = 0.                  ! (kg)       |nitrite transported with water into reach
          real :: no2_out = 0.                 ! (kg)       |nitrite transported with water out of reach
          real :: solp_in = 0.                 ! (kg P)     |soluble pesticide transported with water into reach
          real :: solp_out = 0.                ! (kg P)     |soluble pesticide transported with water out of reach
          real :: chla_in = 0.                 ! (kg)       |amount of chlorophyll a transported into reach
          real :: chla_out = 0.                ! (kg)       |amount of chlorophyll a transported out of reach
          real :: cbod_in = 0.                 ! (kg)       |carbonaceous biochemical oxygen demand of material transported into reach
          real :: cbod_out = 0.                ! (kg)       |carbonaceous biochemical oxygen demand of material transported out of reach
          real :: dis_in = 0.                  ! (kg)       |amount of dissolved oxygen transported into reach
          real :: dis_out = 0.                 ! (kg)       |amount of dissolved oxygen transported out of reach
          real :: solpst_in = 0.               ! (mg pst)   |soluble pesticide transported with water into reach
          real :: solpst_out = 0.              ! (mg pst)   |soluble pesticide transported with water out of reach
          real :: sorbpst_in = 0.              ! (mg pst)   |pesticide sorbed to sediment transported with water into reach
          real :: sorbpst_out = 0.             ! (mg pst)   |pesticide sorbed to sediment transported with water out of reach
          real :: react = 0.                   ! (mg pst)   |loss of pesticide from water from reaction
          real :: volat = 0.                   ! (mg)       |loss of pesticide from water by volatilization
          real :: setlpst = 0.                 ! (mg pst)   |transfer of pesticide from water to river bed sediment by settling
          real :: resuspst = 0.                ! (mg)       |transfer of pesticide from river bed sediment to water by resuspension
          real :: difus = 0.                   ! mg         |transfer of pesticide from water to river bed sediment by diffusion
          real :: reactb = 0.                  ! (mg)       |loss of pesticide from river bed sediment by reaction
          real :: bury = 0.                    ! (mg)       |loss of pesticide from river bed sediment by burial
          real :: sedpest = 0.                 ! mg         |pesticide in river bed sediment
          real :: bacp = 0.                    ! # cfu/100mL  |number of persistent bacteria transported out of reach
          real :: baclp = 0.                   ! # cfu/100mL  |number of less persistent bacteria transported out of reach
          real :: met1 = 0.                    ! kg         |conservative metal #1 transported out of reach
          real :: met2 = 0.                    ! kg         |conservative metal #2 transported out of reach
          real :: met3 = 0.                    ! kg         |conservative metal #3 transported out of reach
          real :: sand_in = 0.                 ! tons       |sand in
          real :: sand_out = 0.                ! tons       |sand out
          real :: silt_in = 0.                 ! tons       |silt_in
          real :: silt_out = 0.                ! tons       |silt_out
          real :: clay_in = 0.                 ! tons       |clay_in
          real :: clay_out = 0.                ! tons       |clay_out
          real :: smag_in = 0.                 ! tons       |small aggregates transported into reach
          real :: smag_out = 0.                ! tons       |small aggregates transported out of reach
          real :: lag_in = 0.                  ! tons       |large aggregates transported into reachlg ag in
          real :: lag_out = 0.                 ! tons       |large aggregates transported out of reach
          real :: grvl_in = 0.                 ! tons       |gravel in
          real :: grvl_out = 0.                ! tons       |gravel out
          real :: bnk_ero = 0.                 ! tons       |bank erosion
          real :: ch_deg = 0.                  ! tons       |channel degradation
          real :: ch_dep = 0.                  ! tons       |channel deposition
          real :: fp_dep = 0.                  ! tons       |flood deposition
          real :: tot_ssed = 0.                ! mg/L       |total suspended sediments
      end type ch_output

      type regional_output_channel
        type (ch_output), dimension (:), allocatable :: ord
      end type regional_output_channel
      type (regional_output_channel), dimension (:), allocatable :: rch_d
      type (regional_output_channel), dimension (:), allocatable :: rch_m
      type (regional_output_channel), dimension (:), allocatable :: rch_y
      type (regional_output_channel), dimension (:), allocatable :: rch_a

      type (ch_output), dimension(:), allocatable, save :: ch_d
      type (ch_output), dimension(:), allocatable, save :: ch_m
      type (ch_output), dimension(:), allocatable, save :: ch_y
      type (ch_output), dimension(:), allocatable, save :: ch_a
      type (ch_output) :: bch_d
      type (ch_output) :: bch_m
      type (ch_output) :: bch_y
      type (ch_output) :: bch_a
      type (ch_output) :: chz

      type ch_header
          character (len=6) :: day =          " jday"
          character (len=6) :: mo =           "   mon"
          character (len=6) :: day_mo =       "   day"
          character (len=5) :: yrc =          "   yr"
          character (len=9) :: isd =          "    unit "
          character (len=8) :: id =           " gis_id "
          character (len=16) :: name =        " name              "
          character(len=16) :: flo_in =       "      flo_in    "        ! (ha-m)
          character(len=15) :: flo_out =      "    flo_out    "         ! (ha-m)
          character(len=15) :: evap =         "        evap     "       ! (ha-m)
          character(len=15) :: tloss =        "     tloss     "         ! (ha-m)
          character(len=15) :: sed_in =       "      sed_in    "        ! (tons)
          character(len=15) :: sed_out=       "     sed_out    "        ! (tons)
          character(len=15) :: sed_conc =     "    sed_conc    "        ! (mg/L)
          character(len=15) :: orgn_in =      "    orgn_in   "          ! (kg N)
          character(len=15) :: orgn_out =     "    orgn_out    "        ! (kg N)
          character(len=15) :: orgp_in =      "      orgp_in      "     ! (kg P)
          character(len=15) :: orgp_out =     "      orgp_out  "        ! (kg P)
          character(len=15) :: no3_in =       "       no3_in   "        ! (kg N)
          character(len=15) :: no3_out =      "      no3_out   "        ! (kg N)
          character(len=15) :: nh4_in =       "        nh4_in  "        ! (kg)
          character(len=15) :: nh4_out=       "       nh4_out  "        ! (kg)
          character(len=15) :: no2_in =       "        no2_in  "        ! (kg)
          character(len=15) :: no2_out =      "       no2_out  "        ! (kg)
          character(len=15) :: solp_in =      "      solp_in   "        ! (kg P)
          character(len=15) :: solp_out =     "     solp_out   "        ! (kg P)
          character(len=15) :: chla_in =      "       chla_in  "        ! (kg)
          character(len=15) :: chla_out =     "      chla_out  "        ! (kg)
          character(len=15) :: cbod_in =      "       cbod_in  "        ! (kg)
          character(len=15) :: cbod_out =     "      cbod_out  "        ! (kg)
          character(len=15) :: dis_in =       "        dis_in  "        ! (kg)
          character(len=15) :: dis_out =      "       dis_out  "        ! (kg)
          character(len=15) :: solpst_in =    "   solpst_in    "        ! (mg pst)
          character(len=15) :: solpst_out =   "   solpst_out   "        ! (mg pst)
          character(len=15) :: sorbpst_in =   "  sorbpst_in    "        ! (mg pst)
          character(len=15) :: sorbpst_out=   "  sorbpst_out   "        ! (mg pst)
          character(len=15) :: react =        "     react      "        ! (mg pst)
          character(len=15) :: volat =        "        volat   "        ! (mg)
          character(len=15) :: setlpst =      "   setlpst      "        ! (mg pst)
          character(len=15) :: resuspst =     "     resuspst   "        ! (mg)
          character(len=15) :: difus =        "     difus      "        ! (mg pst)
          character(len=15) :: reactb =       "    reactb      "        ! pst/sed (mg)
          character(len=15) :: bury =         "     bury       "        ! pst bury (mg)
          character(len=15) :: sedpest =      "      sedpest   "        ! pst in rivbed sed mg
          character(len=15) :: bacp =         "           bacp"        ! persistent bact out
          character(len=15) :: baclp =        "          baclp"        ! lpersistent bact out
          character(len=15) :: met1 =         "           met1"        ! cmetal #1
          character(len=15) :: met2 =         "           met2"        ! cmetal #2
          character(len=15) :: met3 =         "           met3"        ! cmetal #3
          character(len=15) :: sand_in =      "        sand_in"        ! sand in
          character(len=15) :: sand_out =     "       sand_out"        ! sand out
          character(len=15) :: silt_in =      "        silt_in"        ! silt_in
          character(len=15) :: silt_out =     "       silt_out"        ! silt_out
          character(len=15) :: clay_in =      "        clay_in"        ! clay_in
          character(len=15) :: clay_out =     "       clay_out"        ! clay_out
          character(len=15) :: smag_in =      "        smag_in"        ! sm ag in
          character(len=15) :: smag_out =     "       smag_out"        ! sm ag out
          character(len=15) :: lag_in =       "         lag_in"        ! lg ag in
          character(len=15) :: lag_out =      "        lag_out"        ! lg ag out
          character(len=15) :: grvl_in =      "        grvl_in"        ! gravel in
          character(len=15) :: grvl_out =     "       grvl_out"        ! gravel out
          character(len=15) :: bnk_ero =      "        bnk_ero"        ! bank erosion
          character(len=15) :: ch_deg =       "         ch_deg"        ! channel degradation
          character(len=15) :: ch_dep =       "         ch_dep"        ! channel deposition
          character(len=15) :: fp_dep =       "         fp_dep"        ! flood deposition
          character(len=15) :: tot_ssed =     "       tot_ssed"        ! total suspended sediments
      end type ch_header
      type (ch_header) :: ch_hdr

      type ch_header_units
          character (len=6) :: day         =   "     "
          character (len=6) :: mo          =   "     "
          character (len=6) :: day_mo      =   "      "
          character (len=5) :: yrc         =   "     "
          character (len=9) :: isd         =   "         "
          character (len=8) :: id          =   "        "
          character (len=16) :: name       =   "               "
          character(len=16) :: flo_in      =   "        ha-m"            ! (ha-m)
          character(len=15) :: flo_out     =   "       ha-m"             ! (ha-m)
          character(len=15) :: evap        =   "        ha-m"            ! (ha-m)
          character(len=15) :: tloss       =   "      ha-m"              ! (ha-m)
          character(len=15) :: sed_in      =   "        tons"            ! (tons)
          character(len=15) :: sed_out     =   "        tons"            ! (tons)
          character(len=15) :: sed_conc    =   "        mg/L"            ! (mg/L)
          character(len=15) :: orgn_in     =   "        kgN"             ! (kg N)
          character(len=15) :: orgn_out    =   "         kgN"            ! (kg N)
          character(len=15) :: orgp_in     =   "          kgP"           ! (kg P)
          character(len=15) :: orgp_out    =   "           kgP"          ! (kg P)
          character(len=15) :: no3_in      =   "          kgN"           ! (kg N)
          character(len=15) :: no3_out     =   "          kgN"           ! (kg N)
          character(len=15) :: nh4_in      =   "            kg"          ! (kg)
          character(len=15) :: nh4_out     =   "            kg"          ! (kg)
          character(len=15) :: no2_in      =   "            kg"          ! (kg)
          character(len=15) :: no2_out     =   "            kg"          ! (kg)
          character(len=15) :: solp_in     =   "          kgP"           ! (kg P)
          character(len=15) :: solp_out    =   "          kgP"           ! (kg P)
          character(len=15) :: chla_in     =   "            kg"          ! (kg)
          character(len=15) :: chla_out    =   "            kg"          ! (kg)
          character(len=15) :: cbod_in     =   "            kg"           ! (kg)
          character(len=15) :: cbod_out    =   "            kg"           ! (kg)
          character(len=15) :: dis_in      =   "            kg"           ! (kg)
          character(len=15) :: dis_out     =   "            kg"           ! (kg)
          character(len=15) :: solpst_in   =   "      mg_pst"           ! (mg pst)
          character(len=15) :: solpst_out  =   "       mg_pst"           ! (mg pst)
          character(len=15) :: sorbpst_in  =   "      mg_pst"           ! (mg pst)
          character(len=15) :: sorbpst_out =   "       mg_pst"           ! (mg pst)
          character(len=15) :: react       =   "    mg_pst"           ! (mg pst)
          character(len=15) :: volat       =   "           mg"           ! (mg)
          character(len=15) :: setlpst     =   "    mg_pst"           ! (mg pst)
          character(len=15) :: resuspst    =   "           mg"           ! (mg)
          character(len=15) :: difus       =   "    mg_pst"           ! (mg pst)
          character(len=15) :: reactb      =   "        mg"           ! pst/sed (mg)
          character(len=15) :: bury        =   "       mg"           ! pst bury (mg)
          character(len=15) :: sedpest     =   "           mg"           ! pst in rivbed sed mg
          character(len=15) :: bacp        =   "            ----"        ! persistent bact out
          character(len=15) :: baclp       =   "            ----"        ! lpersistent bact out
          character(len=15) :: met1        =   "            ----"        ! cmetal #1
          character(len=15) :: met2        =   "            ----"        ! cmetal #2
          character(len=15) :: met3        =   "            ----"        ! cmetal #3
          character(len=15) :: sand_in     =   "            ----"        ! sand in
          character(len=15) :: sand_out    =   "            ----"        ! sand out
          character(len=15) :: silt_in     =   "            ----"        ! silt_in
          character(len=15) :: silt_out    =   "            ----"        ! silt_out
          character(len=15) :: clay_in     =   "            ----"        ! clay_in
          character(len=15) :: clay_out    =   "            ----"        ! clay_out
          character(len=15) :: smag_in     =   "            ----"        ! sm ag in
          character(len=15) :: smag_out    =   "            ----"        ! sm ag out
          character(len=15) :: lag_in      =   "            ----"        ! lg ag in
          character(len=15) :: lag_out     =   "            ----"        ! lg ag out
          character(len=15) :: grvl_in     =   "            ----"        ! gravel in
          character(len=15) :: grvl_out    =   "            ----"        ! gravel out
          character(len=15) :: bnk_ero     =   "            ----"        ! bank erosion
          character(len=15) :: ch_deg      =   "            ----"        ! channel degradation
          character(len=15) :: ch_dep      =   "            ----"        ! channel deposition
          character(len=15) :: fp_dep      =   "            ----"        ! flood deposition
          character(len=15) :: tot_ssed    =   "            ----"        ! total suspended sediments
      end type ch_header_units
      type (ch_header_units) :: ch_hdr_units

      interface operator (+)
        module procedure ch_add
      end interface

      interface operator (/)
        module procedure ch_div
      end interface

      interface operator (*)
        module procedure ch_mult
      end interface

      contains

      function ch_add(cho1,cho2) result (cho3)
      type (ch_output),  intent (in) :: cho1
      type (ch_output),  intent (in) :: cho2
      type (ch_output) :: cho3
       cho3%flo_in = cho1%flo_in + cho2%flo_in
       cho3%flo_out = cho1%flo_out + cho2%flo_out
       cho3%evap = cho1%evap + cho2%evap
       cho3%tloss = cho1%tloss + cho2%tloss
       cho3%sed_in = cho1%sed_in + cho2%sed_in
       cho3%sed_out = cho1%sed_out + cho2%sed_out
       cho3%sed_conc = cho1%sed_conc + cho2%sed_conc
       cho3%orgn_in =cho1%orgn_in + cho2%orgn_in
       cho3%orgn_out = cho1%orgn_out + cho2%orgn_out
       cho3%orgp_in = cho1%orgp_in + cho2%orgp_in
       cho3%orgp_out = cho1%orgp_out + cho2%orgp_out
       cho3%no3_in = cho1%no3_in + cho2%no3_in
       cho3%no3_out = cho1%no3_out + cho2%no3_out
       cho3%nh4_in = cho1%nh4_in + cho2%nh4_in
       cho3%nh4_out = cho1%nh4_out + cho2%nh4_out
       cho3%no2_in = cho1%no2_in + cho2%no2_in
       cho3%no2_out = cho1%no2_out + cho2%no2_out
       cho3%solp_in = cho1%solp_in + cho2%solp_in
       cho3%solp_out = cho1%solp_out + cho2%solp_out
       cho3%chla_in = cho1%chla_in + cho2%chla_in
       cho3%chla_out = cho1%chla_out + cho2%chla_out
       cho3%cbod_in = cho1%cbod_in + cho2%cbod_in
       cho3%cbod_out = cho1%cbod_out + cho2%cbod_out
       cho3%dis_in = cho1%dis_in + cho2%dis_in
       cho3%dis_out = cho1%dis_out + cho2%dis_out
       cho3%solpst_in = cho1%solpst_in + cho2%solpst_in
       cho3%solpst_out = cho1%solpst_out + cho2%solpst_out
       cho3%sorbpst_in = cho1%sorbpst_in + cho2%sorbpst_in
       cho3%sorbpst_out = cho1%sorbpst_out + cho2%sorbpst_out
       cho3%react = cho1%react + cho2%react
       cho3%volat = cho1%volat + cho2%volat
       cho3%setlpst = cho1%setlpst + cho2%setlpst
       cho3%resuspst = cho1%resuspst + cho2%resuspst
       cho3%difus = cho1%difus + cho2%difus
       cho3%reactb = cho1%reactb + cho2%reactb
       cho3%bury = cho1%bury + cho2%bury
       cho3%sedpest = cho1%sedpest + cho2%sedpest
       cho3%bacp = cho1%bacp + cho2%bacp
       cho3%baclp = cho1%baclp + cho2%baclp
       cho3%met1 = cho1%met1 + cho2%met1
       cho3%met2 = cho1%met2 + cho2%met2
       cho3%met3 = cho1%met3 + cho2%met3
       cho3%sand_in = cho1%sand_in + cho2%sand_in
       cho3%sand_out = cho1%sand_out + cho2%sand_out
       cho3%silt_in = cho1%silt_in + cho2%silt_in
       cho3%silt_out = cho1%silt_out + cho2%silt_out
       cho3%clay_in = cho1%clay_in + cho2%clay_in
       cho3%clay_out = cho1%clay_out + cho2%clay_out
       cho3%smag_in = cho1%smag_in + cho2%smag_in
       cho3%smag_out = cho1%smag_out + cho2%smag_out
       cho3%lag_in = cho1%lag_in + cho2%lag_in
       cho3%lag_out = cho1%lag_out + cho2%lag_out
       cho3%grvl_in = cho1%grvl_in + cho2%grvl_in
       cho3%grvl_out = cho1%grvl_out + cho2%grvl_out
       cho3%bnk_ero = cho1%bnk_ero + cho2%bnk_ero
       cho3%ch_deg = cho1%ch_deg + cho2%ch_deg
       cho3%fp_dep = cho1%fp_dep + cho2%fp_dep
       cho3%tot_ssed = cho1%tot_ssed + cho2%tot_ssed
      end function

      function ch_div (ch1,const) result (ch2)
        type (ch_output), intent (in) :: ch1
        real, intent (in) :: const
        type (ch_output) :: ch2
        ch2%flo_in = ch1%flo_in / const
        ch2%flo_out = ch1%flo_out / const
        ch2%evap = ch1%evap / const
        ch2%tloss = ch1%tloss / const
        ch2%sed_in = ch1%sed_in / const
        ch2%sed_out = ch1%sed_out / const
        ch2%sed_conc = ch1%sed_conc / const
        ch2%orgn_in = ch1%orgn_in / const
        ch2%orgn_out = ch1%orgn_out / const
        ch2%orgp_in = ch1%orgp_in / const
        ch2%orgp_out = ch1%orgp_out / const
        ch2%no3_in = ch1%no3_in / const
        ch2%no3_out = ch1%no3_out / const
        ch2%nh4_in = ch1%nh4_in / const
        ch2%nh4_out = ch1%nh4_out / const
        ch2%no2_in = ch1%no2_in / const
        ch2%no2_out = ch1%no2_out / const
        ch2%solp_in = ch1%solp_in / const
        ch2%solp_out = ch1%solp_out / const
        ch2%chla_in = ch1%chla_in / const
        ch2%chla_out = ch1%chla_out / const
        ch2%cbod_in = ch1%cbod_in / const
        ch2%cbod_out = ch1%cbod_out / const
        ch2%dis_in = ch1%dis_in / const
        ch2%dis_out = ch1%dis_out / const
        ch2%solpst_in = ch1%solpst_in / const
        ch2%solpst_out = ch1%solpst_out / const
        ch2%sorbpst_in = ch1%sorbpst_in / const
        ch2%sorbpst_out = ch1%sorbpst_out / const
        ch2%react = ch1%react / const
        ch2%volat = ch1%volat / const
        ch2%setlpst = ch1%setlpst / const
        ch2%resuspst = ch1%resuspst / const
        ch2%difus = ch1%difus / const
        ch2%reactb = ch1%reactb / const
        ch2%bury = ch1%bury / const
        ch2%sedpest = ch1%sedpest / const
        ch2%bacp = ch1%bacp / const
        ch2%baclp = ch1%baclp / const
        ch2%met1 = ch1%met1 / const
        ch2%met2 = ch1%met2 / const
        ch2%met3 = ch1%met3 / const
        ch2%sand_in = ch1%sand_in / const
        ch2%sand_out = ch1%sand_out / const
        ch2%silt_in = ch1%silt_in / const
        ch2%silt_out = ch1%silt_out / const
        ch2%clay_in = ch1%clay_in / const
        ch2%clay_out = ch1%clay_out / const
        ch2%smag_in = ch1%smag_in / const
        ch2%smag_out = ch1%smag_out / const
        ch2%lag_in = ch1%lag_in / const
        ch2%lag_out = ch1%lag_out / const
        ch2%grvl_in = ch1%grvl_in / const
        ch2%grvl_out = ch1%grvl_out / const
        ch2%bnk_ero = ch1%bnk_ero / const
        ch2%ch_deg = ch1%ch_deg / const
        ch2%fp_dep = ch1%fp_dep / const
        ch2%tot_ssed = ch1%tot_ssed / const
      end function ch_div

      function ch_mult (const, chn1) result (chn2)
        type (ch_output), intent (in) :: chn1
        real, intent (in) :: const
        type (ch_output) :: chn2
        chn2%flo_in = const * chn1%flo_in
        chn2%flo_out = const * chn1%flo_out
        chn2%evap = const * chn1%evap
        chn2%tloss = const * chn1%tloss
        chn2%sed_in = const * chn1%sed_in
        chn2%sed_out = const * chn1%sed_out
        chn2%sed_conc = const * chn1%sed_conc
        chn2%orgn_in = const * chn1%orgn_in
        chn2%orgn_out = const * chn1%orgn_out
        chn2%orgp_in = const * chn1%orgp_in
        chn2%orgp_out = const * chn1%orgp_out
        chn2%no3_in = const * chn1%no3_in
        chn2%no3_out = const * chn1%no3_out
        chn2%nh4_in = const * chn1%nh4_in
        chn2%nh4_out = const * chn1%nh4_out
        chn2%no2_in = const * chn1%no2_in
        chn2%no2_out = const * chn1%no2_out
        chn2%solp_in = const * chn1%solp_in
        chn2%solp_out = const * chn1%solp_out
        chn2%chla_in = const * chn1%chla_in
        chn2%chla_out = const * chn1%chla_out
        chn2%cbod_in = const * chn1%cbod_in
        chn2%cbod_out = const * chn1%cbod_out
        chn2%dis_in = const * chn1%dis_in
        chn2%dis_out = const * chn1%dis_out
        chn2%solpst_in = const * chn1%solpst_in
        chn2%solpst_out = const * chn1%solpst_out
        chn2%sorbpst_in = const * chn1%sorbpst_in
        chn2%sorbpst_out = const * chn1%sorbpst_out
        chn2%react = const * chn1%react
        chn2%bury = const * chn1%bury
        chn2%sedpest = const * chn1%sedpest
        chn2%bacp = const * chn1%bacp
        chn2%baclp = const * chn1%baclp
        chn2%met1 = const * chn1%met1
        chn2%met2 = const * chn1%met2
        chn2%met3 = const * chn1%met3
        chn2%sand_in = const * chn1%sand_in
        chn2%sand_out = const * chn1%sand_out
        chn2%silt_in = const * chn1%silt_in
        chn2%silt_out = const * chn1%silt_out
        chn2%clay_in = const * chn1%clay_in
        chn2%clay_out = const * chn1%clay_out
        chn2%smag_in = const * chn1%smag_in
        chn2%smag_out = const * chn1%smag_out
        chn2%lag_in = const * chn1%lag_in
        chn2%lag_out = const * chn1%lag_out
        chn2%grvl_in = const * chn1%grvl_in
        chn2%grvl_out = const * chn1%grvl_out
        chn2%bnk_ero = const * chn1%bnk_ero
        chn2%ch_deg = const * chn1%ch_deg
        chn2%fp_dep = const * chn1%fp_dep
        chn2%tot_ssed = const * chn1%tot_ssed
      end function ch_mult

      end module channel_module


text:
!****h* root/fortran/test/tf_F90.f90
!
! NAME
!  tf_F90.f90
!
! FUNCTION
!  Module for when the compiler is not F2003 or F2008 compliant.
!  Needed by tf.f90 for the test programs.
!
! COPYRIGHT
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!   Copyright by The HDF Group.                                               *
!   Copyright by the Board of Trustees of the University of Illinois.         *
!   All rights reserved.                                                      *
!                                                                             *
!   This file is part of HDF5.  The full HDF5 copyright notice, including     *
!   terms governing use, modification, and redistribution, is contained in    *
!   the COPYING file, which can be found at the root of the source code       *
!   distribution tree, or in https://support.hdfgroup.org/ftp/HDF5/releases.  *
!   If you do not have access to either file, you may request a copy from     *
!   help@hdfgroup.org.                                                        *
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!
!*****
MODULE TH5_MISC_PROVISIONAL
  IMPLICIT NONE

  INTEGER, PARAMETER :: sp = KIND(0.0)
  INTEGER, PARAMETER :: dp = KIND(0.D0)

END MODULE TH5_MISC_PROVISIONAL

max_stars_repo_name: GWU-CFD/flash-distro

max_issues_repo_name: GWU-CFD/flash-distro

max_forks_repo_name: GWU-CFD/flash-distro


text:
MODULE vtk_attributes
    USE Precision, ONLY : i4k, r8k
    USE Misc,      ONLY : def_len
    USE VTK_DataArray_element, ONLY : DataArray_dt
    IMPLICIT NONE
    !! author: Ian Porter
    !! date: 12/13/2017
    !!
    !! This module contains the dataset attributes for vtk format
    !!
    !! The following dataset attributes are available:
    !! 1) scalars
    !! 2) vectors
    !! 3) normals
    !! 4) texture coordinates (1D, 2D & 3D)
    !! 5) 3x3 tensors
    !! 6) field data
    !!
    PRIVATE
    PUBLIC :: attribute, attributes, scalar, vector, normal, texture, tensor, field, field_data_array

    CHARACTER(LEN=*), PARAMETER :: default = 'default'     !! Default table name

    TYPE, ABSTRACT :: attribute
        !! Abstract DT of attribute information
        INTEGER(i4k) :: nvals = 0_i4k
        CHARACTER(LEN=:), ALLOCATABLE :: dataname
        CHARACTER(LEN=:), ALLOCATABLE :: datatype
    CONTAINS
        PROCEDURE(abs_read),  DEFERRED, PUBLIC :: read
        PROCEDURE(abs_write), DEFERRED, PUBLIC :: write
        PROCEDURE, NON_OVERRIDABLE,     PUBLIC :: init => initialize  !! Initialize the attribute
        PROCEDURE, PRIVATE :: check_for_diffs
        GENERIC :: OPERATOR(.diff.) => check_for_diffs
        PROCEDURE :: convert_to_dataarray
        !PROCEDURE(abs_get_name), DEFERRED, PUBLIC :: get_name
    END TYPE attribute

    TYPE, EXTENDS(attribute) :: scalar
        !! Scalar attribute DT
        INTEGER(i4k) :: numcomp = 0
        CHARACTER(LEN=:), ALLOCATABLE :: tablename
        INTEGER(i4k), DIMENSION(:), ALLOCATABLE :: ints
        REAL(r8k),    DIMENSION(:), ALLOCATABLE :: reals
    CONTAINS
        PROCEDURE :: read  => scalar_read
        PROCEDURE :: write => scalar_write
        PROCEDURE :: setup => scalar_setup
        PROCEDURE, PRIVATE :: check_for_diffs => scalar_check_for_diffs
        PROCEDURE :: convert_to_dataarray => scalar_convert_to_dataarray
    END TYPE scalar

    TYPE, EXTENDS(attribute) :: vector
        !! Vector attribute DT
        INTEGER(i4k), DIMENSION(:,:), ALLOCATABLE :: i_vector
        REAL(r8k),    DIMENSION(:,:), ALLOCATABLE :: r_vector
    CONTAINS
        PROCEDURE :: read  => vector_read
        PROCEDURE :: write => vector_write
        PROCEDURE :: setup => vector_setup
        PROCEDURE, PRIVATE :: check_for_diffs => vector_check_for_diffs
    END TYPE vector

    TYPE, EXTENDS(attribute) :: normal
        !! Normal attribute DT
        INTEGER(i4k), DIMENSION(:,:), ALLOCATABLE :: i_normal
        REAL(r8k),    DIMENSION(:,:), ALLOCATABLE :: r_normal
    CONTAINS
        PROCEDURE :: read  => normal_read
        PROCEDURE :: write => normal_write
        PROCEDURE :: setup => normal_setup
        PROCEDURE, PRIVATE :: check_for_diffs => normal_check_for_diffs
    END TYPE normal

    TYPE, EXTENDS(attribute) :: texture
        !! Texture attribute DT
        INTEGER(i4k), DIMENSION(:,:), ALLOCATABLE :: i_texture
        REAL(r8k),    DIMENSION(:,:), ALLOCATABLE :: r_texture
    CONTAINS
        PROCEDURE :: read  => texture_read
        PROCEDURE :: write => texture_write
        PROCEDURE :: setup => texture_setup
        PROCEDURE, PRIVATE :: check_for_diffs => texture_check_for_diffs
    END TYPE texture

    TYPE :: i_tensor_array
        !! Tensor integer data DT
        INTEGER(i4k), DIMENSION(3,3) :: val = 0_i4k
    END TYPE i_tensor_array

    TYPE :: r_tensor_array
        !! Tensor real data DT
        REAL(r8k),    DIMENSION(3,3) :: val = 0.0_r8k
    END TYPE r_tensor_array

    TYPE, EXTENDS(attribute) :: tensor
        !! Tensor attribute DT
        TYPE(i_tensor_array), DIMENSION(:), ALLOCATABLE :: i_tensor
        TYPE(r_tensor_array), DIMENSION(:), ALLOCATABLE :: r_tensor
    CONTAINS
        PROCEDURE :: read  => tensor_read
        PROCEDURE :: write => tensor_write
        PROCEDURE :: setup => tensor_setup
        PROCEDURE, PRIVATE :: check_for_diffs => tensor_check_for_diffs
    END TYPE tensor

    TYPE :: field_data_array
        !! Field data DT
        CHARACTER(LEN=:), ALLOCATABLE :: name
        INTEGER(i4k) :: numComponents = 0
        INTEGER(i4k) :: numTuples     = 0
        CHARACTER(LEN=:), ALLOCATABLE :: datatype
        REAL(r8k), DIMENSION(:,:), ALLOCATABLE :: data
    END TYPE field_data_array

    TYPE, EXTENDS(attribute) :: field
        !! Field attribute DT
        TYPE(field_data_array), DIMENSION(:), ALLOCATABLE :: array
    CONTAINS
        PROCEDURE :: read  => field_read
        PROCEDURE :: write => field_write
        PROCEDURE :: setup => field_setup
        PROCEDURE, PRIVATE :: check_for_diffs => field_check_for_diffs
    END TYPE field

    TYPE :: attributes
        CLASS(attribute), ALLOCATABLE :: attribute
    END TYPE attributes

    INTERFACE

        MODULE SUBROUTINE abs_read (me, unit)
        IMPLICIT NONE
        !! author: Ian Porter
        !! date: 12/13/2017
        !!
        !! Abstract for reading an attribute
        !!
        CLASS(attribute), INTENT(OUT) :: me
        INTEGER(i4k),     INTENT(IN)  :: unit

        END SUBROUTINE abs_read

        MODULE SUBROUTINE abs_write (me, unit)
        IMPLICIT NONE
        !! author: Ian Porter
        !! date: 12/13/2017
        !!
        !! Abstract for writing an attribute
        !!
        CLASS(attribute), INTENT(IN) :: me
        INTEGER(i4k),     INTENT(IN) :: unit

        END SUBROUTINE abs_write

        MODULE SUBROUTINE initialize (me, dataname, datatype, numcomp, tablename, int1d, int2d, int3d, &
          &                           real1d, real2d, real3d, field_arrays)
        IMPLICIT NONE
        !! author: Ian Porter
        !! date: 12/13/2017
        !!
        !! Abstract for performing the set-up of an attribute
        !!
        CLASS(attribute), INTENT(OUT) :: me
        CHARACTER(LEN=*), INTENT(IN)  :: dataname
        INTEGER(i4k),     INTENT(IN), OPTIONAL :: numcomp
        CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: datatype, tablename
        INTEGER(i4k), DIMENSION(:),     INTENT(IN), OPTIONAL :: int1d
        INTEGER(i4k), DIMENSION(:,:),   INTENT(IN), OPTIONAL :: int2d
        INTEGER(i4k), DIMENSION(:,:,:), INTENT(IN), OPTIONAL :: int3d
        REAL(r8k),    DIMENSION(:),     INTENT(IN), OPTIONAL :: real1d
        REAL(r8k),    DIMENSION(:,:),   INTENT(IN), OPTIONAL :: real2d
        REAL(r8k),    DIMENSION(:,:,:), INTENT(IN), OPTIONAL :: real3d
        TYPE(field_data_array), DIMENSION(:), INTENT(IN), OPTIONAL :: field_arrays

        END SUBROUTINE initialize

        MODULE FUNCTION check_for_diffs (me, you) RESULT (diffs)
        IMPLICIT NONE
        !! author: Ian Porter
        !! date: 12/13/2017
        !!
        !! Function checks for differences in an attribute
        !!
        CLASS(attribute), INTENT(IN) :: me, you
        LOGICAL                      :: diffs

        END FUNCTION check_for_diffs

        MODULE FUNCTION convert_to_dataarray (me) RESULT (array)
        IMPLICIT NONE
        !! author: Ian Porter
        !! date: 07/20/2019
        !!
        !! Function converts an attribute to a dataarray
        !!
        CLASS(attribute), INTENT(IN) :: me
        TYPE(DataArray_dt)           :: array

        END FUNCTION convert_to_dataarray
!********
! Scalars
!********
        MODULE SUBROUTINE scalar_read (me, unit)
        IMPLICIT NONE
        !! author: Ian Porter
        !! date: 12/13/2017
        !!
        !! Subroutine performs the read for a scalar attribute
        !!
        CLASS(scalar), INTENT(OUT) :: me
        INTEGER(i4k),  INTENT(IN)  :: unit

        END SUBROUTINE scalar_read

        MODULE SUBROUTINE scalar_write (me, unit)
        IMPLICIT NONE
        !! author: Ian Porter
        !! date: 12/13/2017
        !!
        !! Subroutine performs the write for a scalar attribute
        !!
        CLASS(scalar), INTENT(IN) :: me
        INTEGER(i4k),  INTENT(IN) :: unit

        END SUBROUTINE scalar_write

        MODULE SUBROUTINE scalar_setup (me, dataname, datatype, numcomp, tablename, int1d, real1d)
        IMPLICIT NONE
        !! author: Ian Porter
        !! date: 12/13/2017
        !!
        !! Subroutine performs the set-up for a scalar attribute
        !!
        CLASS(scalar),    INTENT(OUT) :: me
        CHARACTER(LEN=*), INTENT(IN)  :: dataname
        INTEGER(i4k),     INTENT(IN), OPTIONAL :: numcomp
        CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: datatype
        CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: tablename
        INTEGER(i4k), DIMENSION(:), INTENT(IN), OPTIONAL :: int1d
        REAL(r8k),    DIMENSION(:), INTENT(IN), OPTIONAL :: real1d

        END SUBROUTINE scalar_setup

        MODULE FUNCTION scalar_check_for_diffs (me, you) RESULT (diffs)
        IMPLICIT NONE
        !! author: Ian Porter
        !! date: 12/13/2017
        !!
        !! Function checks for differences in a scalar attribute
        !!
        CLASS(scalar),    INTENT(IN) :: me
        CLASS(attribute), INTENT(IN) :: you
        LOGICAL                      :: diffs

        END FUNCTION scalar_check_for_diffs

        MODULE FUNCTION scalar_convert_to_dataarray (me) RESULT (array)
        IMPLICIT NONE
        !! author: Ian Porter
        !! date: 07/20/2019
        !!
        !! Function converts an attribute to a dataarray
        !!
        CLASS(scalar), INTENT(IN) :: me
        TYPE(DataArray_dt)        :: array

        END FUNCTION scalar_convert_to_dataarray
!********
! Vectors
!********
        MODULE SUBROUTINE vector_read (me, unit)
        IMPLICIT NONE
        !! author: Ian Porter
        !! date: 12/14/2017
        !!
        !! Subroutine performs the read for a vector attribute
        !!
        CLASS(vector), INTENT(OUT) :: me
        INTEGER(i4k),  INTENT(IN)  :: unit

        END SUBROUTINE vector_read

        MODULE SUBROUTINE vector_write (me, unit)
        IMPLICIT NONE
        !! author: Ian Porter
        !! date: 12/13/2017
        !!
        !! Subroutine performs the write for a vector attribute
        !!
        CLASS(vector), INTENT(IN) :: me
        INTEGER(i4k),  INTENT(IN) :: unit

        END SUBROUTINE vector_write

        MODULE SUBROUTINE vector_setup (me, dataname, datatype, int2d, real2d)
        IMPLICIT NONE
        !! author: Ian Porter
        !! date: 12/14/2017
        !!
        !! Subroutine performs the set-up for a vector attribute
        !!
        CLASS(vector),    INTENT(OUT) :: me
        CHARACTER(LEN=*), INTENT(IN)  :: dataname
        CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: datatype
        INTEGER(i4k), DIMENSION(:,:), INTENT(IN), OPTIONAL :: int2d
        REAL(r8k),    DIMENSION(:,:), INTENT(IN), OPTIONAL :: real2d

        END SUBROUTINE vector_setup

        MODULE FUNCTION vector_check_for_diffs (me, you) RESULT (diffs)
        IMPLICIT NONE
        !! author: Ian Porter
        !! date: 12/14/2017
        !!
        !! Function checks for differences in a vector attribute
        !!
        CLASS(vector),    INTENT(IN) :: me
        CLASS(attribute), INTENT(IN) :: you
        LOGICAL                      :: diffs

        END FUNCTION vector_check_for_diffs
!********
! Normals
!********
        MODULE SUBROUTINE normal_read (me, unit)
        IMPLICIT NONE
        !! author: Ian Porter
        !! date: 12/14/2017
        !!
        !! Subroutine performs the read for a normal attribute
        !!
        CLASS(normal), INTENT(OUT) :: me
        INTEGER(i4k),  INTENT(IN)  :: unit

        END SUBROUTINE normal_read

        MODULE SUBROUTINE normal_write (me, unit)
        IMPLICIT NONE
        !! author: Ian Porter
        !! date: 12/13/2017
        !!
        !! Subroutine performs the write for a normal attribute
        !!
        CLASS(normal), INTENT(IN) :: me
        INTEGER(i4k),  INTENT(IN) :: unit

        END SUBROUTINE normal_write

        MODULE SUBROUTINE normal_setup (me, dataname, datatype, int2d, real2d)
        IMPLICIT NONE
        !! author: Ian Porter
        !! date: 12/14/2017
        !!
        !! Subroutine performs the set-up for a normal attribute
        !!
        CLASS(normal),    INTENT(OUT) :: me
        CHARACTER(LEN=*), INTENT(IN)  :: dataname
        CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: datatype
        INTEGER(i4k), DIMENSION(:,:), INTENT(IN), OPTIONAL :: int2d
        REAL(r8k),    DIMENSION(:,:), INTENT(IN), OPTIONAL :: real2d

        END SUBROUTINE normal_setup

        MODULE FUNCTION normal_check_for_diffs (me, you) RESULT (diffs)
        IMPLICIT NONE
        !! author: Ian Porter
        !! date: 12/14/2017
        !!
        !! Function checks for differences in a normal attribute
        !!
        CLASS(normal),    INTENT(IN) :: me
        CLASS(attribute), INTENT(IN) :: you
        LOGICAL                      :: diffs

        END FUNCTION normal_check_for_diffs
!********
! Textures
!********
        MODULE SUBROUTINE texture_read (me, unit)
        IMPLICIT NONE
        !! author: Ian Porter
        !! date: 12/14/2017
        !!
        !! Subroutine performs the read for a texture attribute
        !!
        CLASS(texture), INTENT(OUT) :: me
        INTEGER(i4k),   INTENT(IN)  :: unit

        END SUBROUTINE texture_read

        MODULE SUBROUTINE texture_write (me, unit)
        IMPLICIT NONE
        !! author: Ian Porter
        !! date: 12/13/2017
        !!
        !! Subroutine performs the write for a texture attribute
        !!
        CLASS(texture), INTENT(IN) :: me
        INTEGER(i4k),   INTENT(IN) :: unit

        END SUBROUTINE texture_write

        MODULE SUBROUTINE texture_setup (me, dataname, datatype, int2d, real2d)
        IMPLICIT NONE
        !! author: Ian Porter
        !! date: 12/14/2017
        !!
        !! Subroutine performs the set-up for a texture attribute
        !!
        CLASS(texture),   INTENT(OUT) :: me
        CHARACTER(LEN=*), INTENT(IN)  :: dataname
        CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: datatype
        INTEGER(i4k), DIMENSION(:,:), INTENT(IN), OPTIONAL :: int2d
        REAL(r8k),    DIMENSION(:,:), INTENT(IN), OPTIONAL :: real2d

        END SUBROUTINE texture_setup

        MODULE FUNCTION texture_check_for_diffs (me, you) RESULT (diffs)
        IMPLICIT NONE
        !! author: Ian Porter
        !! date: 12/14/2017
        !!
        !! Function checks for differences in a texture attribute
        !!
        CLASS(texture),   INTENT(IN) :: me
        CLASS(attribute), INTENT(IN) :: you
        LOGICAL                      :: diffs

        END FUNCTION texture_check_for_diffs
!********
! Tensors
!********
        MODULE SUBROUTINE tensor_read (me, unit)
        IMPLICIT NONE
        !! author: Ian Porter
        !! date: 12/14/2017
        !!
        !! Subroutine performs the read for a tensor attribute
        !!
        CLASS(tensor), INTENT(OUT) :: me
        INTEGER(i4k),  INTENT(IN)  :: unit

        END SUBROUTINE tensor_read

        MODULE SUBROUTINE tensor_write (me, unit)
        IMPLICIT NONE
        !! author: Ian Porter
        !! date: 12/13/2017
        !!
        !! Subroutine performs the write for a tensor attribute
        !!
        CLASS(tensor), INTENT(IN) :: me
        INTEGER(i4k),  INTENT(IN) :: unit

        END SUBROUTINE tensor_write

        MODULE SUBROUTINE tensor_setup (me, dataname, datatype, int3d, real3d)
        IMPLICIT NONE
        !! author: Ian Porter
        !! date: 12/14/2017
        !!
        !! Subroutine performs the set-up for a tensor attribute
        !!
        CLASS(tensor),    INTENT(OUT) :: me
        CHARACTER(LEN=*), INTENT(IN)  :: dataname
        CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: datatype
        INTEGER(i4k), DIMENSION(:,:,:), INTENT(IN), OPTIONAL :: int3d
        REAL(r8k),    DIMENSION(:,:,:), INTENT(IN), OPTIONAL :: real3d

        END SUBROUTINE tensor_setup

        MODULE FUNCTION tensor_check_for_diffs (me, you) RESULT (diffs)
        IMPLICIT NONE
        !! author: Ian Porter
        !! date: 12/14/2017
        !!
        !! Function checks for differences in a tensor attribute
        !!
        CLASS(tensor),    INTENT(IN) :: me
        CLASS(attribute), INTENT(IN) :: you
        LOGICAL                      :: diffs

        END FUNCTION tensor_check_for_diffs
!********
! Fields
!********
        MODULE SUBROUTINE field_read (me, unit)
        IMPLICIT NONE
        !! author: Ian Porter
        !! date: 12/14/2017
        !!
        !! Subroutine performs the read for a field attribute
        !!
        CLASS(field), INTENT(OUT) :: me
        INTEGER(i4k), INTENT(IN)  :: unit

        END SUBROUTINE field_read

        MODULE SUBROUTINE field_write (me, unit)
        IMPLICIT NONE
        !! author: Ian Porter
        !! date: 12/13/2017
        !!
        !! Subroutine performs the write for a field attribute
        !!
        CLASS(field), INTENT(IN) :: me
        INTEGER(i4k), INTENT(IN) :: unit

        END SUBROUTINE field_write

        MODULE SUBROUTINE field_setup (me, dataname, datatype, field_arrays)
        IMPLICIT NONE
        !! author: Ian Porter
        !! date: 12/14/2017
        !!
        !! Subroutine performs the set-up for a field attribute
        !!
        CLASS(field),     INTENT(OUT) :: me
        CHARACTER(LEN=*), INTENT(IN)  :: dataname
        CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: datatype
        TYPE(field_data_array), DIMENSION(:), INTENT(IN) :: field_arrays

        END SUBROUTINE field_setup

        MODULE FUNCTION field_check_for_diffs (me, you) RESULT (diffs)
        IMPLICIT NONE
        !! author: Ian Porter
        !! date: 12/14/2017
        !!
        !! Function checks for differences in a field attribute
        !!
        CLASS(field),     INTENT(IN) :: me
        CLASS(attribute), INTENT(IN) :: you
        LOGICAL                      :: diffs

        END FUNCTION field_check_for_diffs

    END INTERFACE

END MODULE vtk_attributes

max_stars_repo_name: zbeekman/vtkmofo

max_issues_repo_name: zbeekman/vtkmofo

max_forks_repo_name: zbeekman/vtkmofo


text:
function gtcpnt(pntnam    ,gdp       )
!----- GPL ---------------------------------------------------------------------
!
!  Copyright (C)  Stichting Deltares, 2011-2016.
!
!  This program is free software: you can redistribute it and/or modify
!  it under the terms of the GNU General Public License as published by
!  the Free Software Foundation version 3.
!
!  This program is distributed in the hope that it will be useful,
!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!  GNU General Public License for more details.
!
!  You should have received a copy of the GNU General Public License
!  along with this program.  If not, see <http://www.gnu.org/licenses/>.
!
!  contact: delft3d.support@deltares.nl
!  Stichting Deltares
!  P.O. Box 177
!  2600 MH Delft, The Netherlands
!
!  All indications and logos of, and references to, "Delft3D" and "Deltares"
!  are registered trademarks of Stichting Deltares, and remain the property of
!  Stichting Deltares. All rights reserved.
!
!-------------------------------------------------------------------------------
!  $Id: gtcpnt.f90 5717 2016-01-12 11:35:24Z mourits $
!  $HeadURL: https://svn.oss.deltares.nl/repos/delft3d/tags/6686/src/engines_gpl/flow2d3d/packages/data/src/allocation/gtcpnt.f90 $
!!--description-----------------------------------------------------------------
!
!   Look up a character pointer name using FSM
!   and return an index into the appropriate array.
!
!!--pseudo code and references--------------------------------------------------
!
!   - Determine the length of the pointer name.
!   - Call the FSM GETPTR routine.
!   - If this returns zero, an error occured: an error message will be printed
!     (by ERRPNT) and the program terminates (in ERRPNT)
!
!!--declarations----------------------------------------------------------------
    use precision
    use globaldata
    !
    implicit none
    !
    type(globdat),target :: gdp
    !
    ! The following list of pointer parameters is used to point inside the gdp structure
    !
    include 'fsm.i'
!
! Global variables
!
    integer(pntrsize) :: gtcpnt
    character(*)    :: pntnam
                                   !!  Character string containing array
                                   !!  name (hence max 6 characters).
!
! Local variables
!
    integer                        :: ind   ! Actual length of character string
!
!! executable statements -------------------------------------------------------
!
    ! Define length of pointer name (array name)
    ind = index(pntnam, ' ')
    if (ind==0) ind = len(pntnam) + 1
    !
    ! Call dynamic array declaration function GETPTR to get requested
    ! array pointer
    !
    gtcpnt = getptr(pntnam(:ind - 1))
    if (gtcpnt==0) then
       call errpnt(pntnam, 'character', 'lookup', gdp)
       ! program is terminated in errpnt
    endif
end function gtcpnt

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


! triangle - ASS configuration
! Copyright 2021 by Eric Conrad.
! License - MIT License - see: https://opensource.org/licenses/MIT
!
!   input - the measure of an angle the length of a pair of sides
!       that don't contain the angle (NB: one of the two sides will
!       be opposite the given angle)
!   remarks - there may be 0, 1 or 2 solutions
!   output -
!       the sides and angles
!       the semiperimeter and the area
!       type of triangle
!
! Corrections:
!  14 Jan 2021 -
!		- added copyright and license statement
!		- corrected text of second error message under "Remaining error checks",
!		  changing "given angles" (plural) to "given angle" (singular).

! When I took geometry in high school, this was called the ASS
! configuration.  In today's more puritannical times, it's called
! the SSA configuration.

      program triangle
      print*,"ASS Configuration"

! Notes:
!   100 grad = 90 degrees = pi/2 (rad) = 1 rt angle = 1/4 rev
!   200 grad = 180 degrees = pi (rad) = 2 rt angles = 1/2 rev

      print*,"Angle measure units (1-radian, 2-degree, 3-grad, 4-rev):"
      read(*,*) imeas
      if(imeas .lt. 1 .or. imeas .gt. 4) then
          print*,"ERROR: unsupported angle measure"
          stop 1
      end if

      print*,"Enter the measure of the given angle (A):"
      read(*,*) alpha
      print*,"Enter the length of the side opposite A (a):"
      read(*,*) a
      print*,"Enter the length of a side adjacent to A (b):"
      read(*,*) b

! First error condition:
! the given sides must be positive in length
      if(a.le.0.0 .or. b.le.0.0) then
          print*,"ERROR: The given sides must be positive in length"
          stop 1
      end if

! Convert the angle measures
      pi = 4*atan(1.d0)
      if(imeas.eq.1) then       ! input in radians
          dalpha = 180 * alpha / pi
          galpha = 200 * alpha / pi
          ralpha = alpha / (2 * pi)
      else if(imeas.eq.2) then  ! input in degrees
          dalpha = alpha
          alpha = pi * dalpha / 180
          galpha = 200 * dalpha / 180
          ralpha = dalpha / 360
      else if(imeas.eq.3) then  ! input in grads
          galpha = alpha
          alpha = pi * galpha / 200
          dalpha = 180 * galpha / 200
          ralpha = dalpha / 400
      else                      ! input in revolutions
          ralpha = alpha
          alpha = 2 * pi * galpha
          dalpha = 360 * galpha
          galpha = 400 * ralpha
      end if

! Remaining error checks:
!   (1) the given angle should positive
      if(alpha.le.0) then
          print*, "ERROR: the given angle must be positive"
          stop 1
      end if
!   (2) the given angle must be less than a straight angle
      if(alpha.ge.pi) then
          print*, "ERROR: the given angle must be less" // &
                "than a straight angle"
          stop 1
      end if

! law of sines
! it's safe for finding sides (apart from numerical errors), but
! must be used carefully when used to find angles...
! there will be either no solutions or two solutions (the principal
! acute angle and its obtuse counterpart) or a double solution (a
! right angle)...

! Step 1) Find angle B using the law of sines
!     sin(B) / B = sin(A) / A

      sinB = b * sin(alpha) / a
      print*, "sin(B) =", sinB
      if(sinB .gt. 1) then
          print*, "There is no triangle satifying these constraints."
          stop 1
      end if

! we have reduced the confirguration to AAS (A, B, a), but there
! are two values of B to consider... and there are still some
! landmines...

! the principal solution
      print*, "Principal solution"
      beta = asin(sinB)
      print*, "a =", a, ", b =", b
      print*, "A =", alpha, ", B =", beta
      if((alpha + beta) .ge. pi) then
          print*, "The two known angles exceed a straight angle."
          print*, "No solution."
          stop 1
      end if

      dbeta = 180 * beta / pi
      gbeta = 200 * beta / pi
      rbeta = beta / (2 * pi)

! the principal solution B gives a real triangle...

      gamm = pi - alpha - beta          ! third angle
      dgamma = 180 * gamm / pi
      ggamma = 200 * gamm / pi
      rgamma = gamm / (2 * pi)

      c = sin(gamm) * a / sin(alpha)

! output solution of triangle
      print*, "------- PRINCIPAL SOLUTION --------------------------"
      print*, "a =", a, "  (given)"
      print*, "b =", b, "  (given)"
      print*, "c =", c
      print*, "A =", alpha, "degrees:", dalpha, "  (given)"
      print*, "B =", beta, "degrees:", dbeta, "  (principal)"
      print*, "C =", gamm, "degrees:", dgamma
      print*, "-----------------------------------------------------"
      print*, "Angles in grads:"
      print*, "     ", galpha, gbeta, ggamma
      print*, "     sum:", galpha + gbeta + ggamma, "  (exp: 200)"
      print*, "Angles in revolutions:"
      print*, "     ", ralpha, rbeta, rgamma
      print*, "     sum:", ralpha + rbeta + rgamma, "  (exp: 0.5)"
      print*, "-----------------------------------------------------"

! output semiperimeter and area
      s = (a + b + c) / 2
      print*, "semiperimeter:", s
      area = sqrt(s * (s-a) * (s-b) * (s-c))    ! Heron's formula
      print*, "area:", area

! angle sum
      sum1 = alpha + beta + gamm
      err1 = sum1 - pi
      sum2 = dalpha + dbeta + dgamma
      err2 = sum2 - 180.0
      print*, "angle sum (radians):", sum1, "   error:", err1
      print*, "angle sum (degrees):", sum2, "   error:", err2

! triangle inequality check:
!           a + b < c; b + c < a; c + a < c
! A sufficient and necessary condition is that no side may
! exceed the semiperimeter.
      if (max(a,b,c) .ge. s) then
          print*,"ERROR: (triangle inequality) max(a,b,c) < s"
      end if

! Nesbitt's inequality
      symsum = a/(b+c) + b/(a+c) + c/(a+b)
      print*,"Nesbitt's symmetric sum:", symsum
      if (symsum.lt.1.5) then
          print*,"   Error: sum < 1.5"
      end if
      if (symsum.ge.2) then
          print*,"   Error: sum > 2"
      end if

! the secondary solution
      print*, "====================================================="
      print*, "Secondary solution"
      beta2 = pi - beta
      print*, "a =", a, ", b =", b
      print*, "A =", alpha, ", B =", beta2
      if(beta2 .eq. beta) then
          print*, "Angle B is a right angle."
          print*, "The solution is unique."
          stop 0
      end if
      beta = beta2
      if((alpha + beta) .ge. pi) then
          print*, "The two known angles exceed a straight angle."
          print*, "No secondary solution."
          stop 0
      end if

      dbeta = 180 * beta / pi
      gbeta = 200 * beta / pi
      rbeta = beta / (2 * pi)

! the secondary (obtuse) solution B gives a real triangle...

      gamm = pi - alpha - beta          ! third angle
      dgamma = 180 * gamm / pi
      ggamma = 200 * gamm / pi
      rgamma = gamm / (2 * pi)

      c = sin(gamm) * a / sin(alpha)

! output solution of triangle
      print*, "------- SECONDARY SOLUTION --------------------------"
      print*, "a =", a, "  (given)"
      print*, "b =", b, "  (given)"
      print*, "c =", c
      print*, "A =", alpha, "degrees:", dalpha, "  (given)"
      print*, "B =", beta, "degrees:", dbeta, "  (secondary)"
      print*, "C =", gamm, "degrees:", dgamma
      print*, "-----------------------------------------------------"
      print*, "Angles in grads:"
      print*, "     ", galpha, gbeta, ggamma
      print*, "     sum:", galpha + gbeta + ggamma, "  (exp: 200)"
      print*, "Angles in revolutions:"
      print*, "     ", ralpha, rbeta, rgamma
      print*, "     sum:", ralpha + rbeta + rgamma, "  (exp: 0.5)"
      print*, "-----------------------------------------------------"

! output semiperimeter and area
      s = (a + b + c) / 2
      print*, "semiperimeter:", s
      area = sqrt(s * (s-a) * (s-b) * (s-c))    ! Heron's formula
      print*, "area:", area

! angle sum
      sum1 = alpha + beta + gamm
      err1 = sum1 - pi
      sum2 = dalpha + dbeta + dgamma
      err2 = sum2 - 180.0
      print*, "angle sum (radians):", sum1, "   error:", err1
      print*, "angle sum (degrees):", sum2, "   error:", err2

! triangle inequality check:
!           a + b < c; b + c < a; c + a < c
! A sufficient and necessary condition is that no side may
! exceed the semiperimeter.
      if (max(a,b,c) .ge. s) then
          print*,"ERROR: (triangle inequality) max(a,b,c) < s"
      end if

! Nesbitt's inequality
      symsum = a/(b+c) + b/(a+c) + c/(a+b)
      print*,"Nesbitt's symmetric sum:", symsum
      if (symsum.lt.1.5) then
          print*,"   Error: sum < 1.5"
      end if
      if (symsum.ge.2) then
          print*,"   Error: sum > 2"
      end if

      end program


text:
C: Add the constrains to the first and second derivatives
C: Therefore, there are three constrains: Cp, DCp, and D2Cp
C: --------------------------------------------------------------
C: For NIST
C: Use a single matrix to calculate all 10 coefficients
C: [A]: M x 10
C: [X]: 10 x 1 (alpha1_low,alpha2_low,alpha3_low,alpha4_low,alpha5_low,
C:              alpha1_hig,alpha2_hig,alpha3_hig,alpha4_hig,alpha5_hig)
C: [C]: M x 1
C: Contrains: CP, dCP/dT, and d2CP/d2T
C:            at the intermediate temperature Tint
C: [B]: 3 x 10
C: [D]: 3 x 1
C----------------------------------------------------------------
C FOR NIST CHEMTHERMO DATA FORMAT
C number of coefficients: N=5 x 2
C number of discrete data point: M>=N, i.e., M=102
C
C use the DGGLSE to solve the linear least squares equality-constrained
C (LlSE) problem.
!***************************A*************************************
! Programmed by John Z. Wen, MIT, June 2005
!****************************************************************
C
      SUBROUTINE NIST(Tint,TNEW,THERM1,THERM2,TEMP1,C1,CH1,CS1)
C
      implicit none
C
      DOUBLE PRECISION Tint,TNEW,THERM1,THERM2
      DOUBLE PRECISION H_298,S_298,DLTH
C
      INTEGER M,N,P,LDA,LDB,LWORK,INFO
      INTEGER MH,NH,PH,LDAH,LDBH,LWORKH
      INTEGER MS,NS,PS,LDAS,LDBS,LWORKS
      DOUBLE PRECISION A1,B1,Bk,C1,D1,X1,WORK1
      DOUBLE PRECISION TEMP1,PHI_NIST1,PHI_NIST2,PHI_NIST3,
     &                   PHI_NIST4,PHI_NIST5,ROL
      DOUBLE PRECISION PHI_NIST_DCP1,PHI_NIST_DCP2,PHI_NIST_DCP3,
     &           PHI_NIST_DCP4,PHI_NIST_DCP5,PHI_NIST_D2CP1,
     &             PHI_NIST_D2CP2,PHI_NIST_D2CP3,PHI_NIST_D2CP4,
     &         PHI_NIST_D2CP5
      DOUBLE PRECISION A_H1,A_S1,A_H2,A_S2
      DOUBLE PRECISION AH1,BH1,CH1,DH1,XH1,WORKH,BHk,CH,CH2
      DOUBLE PRECISION AS1,BS1,CS1,DS1,XS1,WORKS,BSk,CS,CS2
      DOUBLE PRECISION PHI_NIST_H,PHI_NIST_S
      INTEGER MARK1
      DOUBLE PRECISION TEMP,C
C
C Required by LSEC solver:
! LDA >= MAX(1,M)
! LDB >= MAX(1,P)
! LWORK >= max(1,M+N+P)
! P=1 for only one constrain at the intermediate temp.
! P=3 for three constrains at the intermediate temp.
C
C CHANGE 'LWORK' if NECESSARY
      PARAMETER (M=102)
      PARAMETER (LDB=3,P=3)
      PARAMETER (N=10)
      DIMENSION A1(M,N),B1(LDB,N),C1(M-1),D1(P),X1(N),WORK1(M+N+P),
     &            Bk(LDB,N)
      PARAMETER (LDBH=2,PH=2)
      PARAMETER (NH=2)
      PARAMETER (LDBS=2,PS=2)
      PARAMETER (NS=2)
      DIMENSION AH1(M,NH),BH1(LDBH,NH),CH1(M-1),DH1(PH),XH1(NH),
     &            WORKH(M+NH+PH),BHk(LDBH,NH),CH(M),CH2(M)
      DIMENSION AS1(M,NS),BS1(LDBS,NS),CS1(M-1),DS1(PS),XS1(NS),
     &            WORKS(M+NS+PS),BSk(LDBS,NS),CS(M),CS2(M)
      DIMENSION TEMP1(M-1),THERM1(7),THERM2(7)
      DIMENSION TEMP(M),C(M)
      INTEGER I,J
C
! GAS CONSTANT
      COMMON /GAS/ROL,H_298,S_298,DLTH
C
      EXTERNAL PHI_NIST1,PHI_NIST2,PHI_NIST3,PHI_NIST4,PHI_NIST5,
     &         PHI_NIST_DCP1,PHI_NIST_DCP2,PHI_NIST_DCP3,
     &         PHI_NIST_DCP4,PHI_NIST_DCP5,PHI_NIST_D2CP1,
     &             PHI_NIST_D2CP2,PHI_NIST_D2CP3,PHI_NIST_D2CP4,
     &         PHI_NIST_D2CP5,PHI_NIST_H,PHI_NIST_S
C
!     WRITE(*,*) 'CALLED NIST'
C
      LDA   = M
      LWORK = M+N+P
      MH    = M
      LDAH  = MH
      LWORKH= MH+NH+PH
      MS    = M
      LDAS  = MS
      LWORKS= MS+NS+PS
C
C FIND 'Tint' and 'Tint1' in sample data
C
      MARK1 = 0
      DO I=1,M-1
        IF (TEMP1(I) .GT. Tint) THEN
           MARK1=I-1
           IF (Tint .NE. TEMP1(I-1)) THEN
                Tint=TEMP1(I-1)
              TNEW=Tint
!             WRITE(*,*) 'Tint has been replaced!'
           ENDIF
           GOTO 995
        ENDIF
      ENDDO
 995  CONTINUE
!     WRITE(*,*) 'MARK1=',MARK1
C
C Rearrange the matrics
C
      DO I=1,M
         TEMP(I)=0.D0
         C   (I)=0.D0
         CH  (I)=0.D0
         CS  (I)=0.D0
         CH2 (I)=0.D0
         CS2 (I)=0.D0
      ENDDO
C
      DO I=1,MARK1
         TEMP(I)=TEMP1(I)
         C (I)=C1(I)
         CH(I)=CH1(I)
         CS(I)=CS1(I)
      ENDDO
      TEMP(MARK1+1)=TEMP(MARK1)
      C (MARK1+1)=C (MARK1)
      CH(MARK1+1)=CH(MARK1)
      CS(MARK1+1)=CS(MARK1)
      DO I=MARK1+2,M
         TEMP(I)=TEMP1(I-1)
         C (I)=C1 (I-1)
         CH(I)=CH1(I-1)
         CS(I)=CS1(I-1)
      ENDDO
C
!     WRITE(*,*) 'Tint=', Tint
!     WRITE(*,*) TEMP
!     WRITE(*,*) C
!     WRITE(*,*) CH
!     WRITE(*,*) CS
!     PAUSE
C
      DO I=1,N
         X1(I)=0.D0
      ENDDO
      DO I=1,N
         DO J=1,LDA
           A1(J,I)=0.D0
         ENDDO
         DO J=1,LDB
           B1(J,I)=0.D0
           Bk(J,I)=0.D0
         ENDDO
      ENDDO
      DO I=1,LWORK
         WORK1(I)=0.D0
      ENDDO
C
! For PHI matrix, first temperature range
        DO I=1,MARK1
           A1(I,1)=PHI_NIST1(TEMP(I))
           A1(I,2)=PHI_NIST2(TEMP(I))
           A1(I,3)=PHI_NIST3(TEMP(I))
           A1(I,4)=PHI_NIST4(TEMP(I))
           A1(I,5)=PHI_NIST5(TEMP(I))
        ENDDO
! For PHI matrix, second temperature range
        DO I=MARK1+1,M
           A1(I,6) =PHI_NIST1(TEMP(I))
           A1(I,7) =PHI_NIST2(TEMP(I))
           A1(I,8) =PHI_NIST3(TEMP(I))
           A1(I,9) =PHI_NIST4(TEMP(I))
           A1(I,10)=PHI_NIST5(TEMP(I))
        ENDDO
C
!     WRITE(*,*) TEMP1
!     WRITE(*,*) C1
!     PAUSE
C
C For Constrains at Tint ------------------------
C
! For CP values
        Bk(1,1) = PHI_NIST1(Tint)
        Bk(1,2) = PHI_NIST2(Tint)
        Bk(1,3) = PHI_NIST3(Tint)
        Bk(1,4) = PHI_NIST4(Tint)
        Bk(1,5) = PHI_NIST5(Tint)
        Bk(1,6) =-PHI_NIST1(Tint)
        Bk(1,7) =-PHI_NIST2(Tint)
        Bk(1,8) =-PHI_NIST3(Tint)
        Bk(1,9) =-PHI_NIST4(Tint)
        Bk(1,10)=-PHI_NIST5(Tint)
        D1(1)=0.D0
C For DCp values
        Bk(2,1) = PHI_NIST_DCP1(Tint)
        Bk(2,2) = PHI_NIST_DCP2(Tint)
        Bk(2,3) = PHI_NIST_DCP3(Tint)
        Bk(2,4) = PHI_NIST_DCP4(Tint)
        Bk(2,5) = PHI_NIST_DCP5(Tint)
        Bk(2,6) =-PHI_NIST_DCP1(Tint)
        Bk(2,7) =-PHI_NIST_DCP2(Tint)
        Bk(2,8) =-PHI_NIST_DCP3(Tint)
        Bk(2,9) =-PHI_NIST_DCP4(Tint)
        Bk(2,10)=-PHI_NIST_DCP5(Tint)
        D1(2)=0.D0
C For D2Cp values
        Bk(3,1) = PHI_NIST_D2CP1(Tint)
        Bk(3,2) = PHI_NIST_D2CP2(Tint)
        Bk(3,3) = PHI_NIST_D2CP3(Tint)
        Bk(3,4) = PHI_NIST_D2CP4(Tint)
        Bk(3,5) = PHI_NIST_D2CP5(Tint)
        Bk(3,6) =-PHI_NIST_D2CP1(Tint)
        Bk(3,7) =-PHI_NIST_D2CP2(Tint)
        Bk(3,8) =-PHI_NIST_D2CP3(Tint)
        Bk(3,9) =-PHI_NIST_D2CP4(Tint)
        Bk(3,10)=-PHI_NIST_D2CP5(Tint)
        D1(3)=0.D0
C
      DO I=1,N
         DO J=1,LDB
            B1(J,I)=Bk(J,I)
         ENDDO
      ENDDO
C
! NIST FORMAT-----------------------------------------
C
      call DGGLSE(M,N,P,A1,LDA,B1,LDB,C,D1,X1,WORK1,LWORK,INFO)
!     WRITE(*,*) 'X1=',X1(1:10)
!     pause
C
C For enthalpy, NEED define 2 coeffs, N=2
C
      DO I=1,NH
         XH1(I)=0.D0
      ENDDO
      DO I=1,NH
         DO J=1,LDAH
           AH1(J,I)=0.D0
         ENDDO
         DO J=1,LDBH
           BH1(J,I)=0.D0
           BHk(J,I)=0.D0
         ENDDO
      ENDDO
      DO I=1,LWORKH
         WORKH(I)=0.D0
      ENDDO
C
      DO I=1,MARK1
         CH2(I)=CH(I)+ DLTH - H_298 - PHI_NIST_H(TEMP(I),X1(1:5))
      ENDDO
      DO I=MARK1+1,MH
         CH2(I)=CH(I)+ DLTH - H_298 - PHI_NIST_H(TEMP(I),X1(6:10))
      ENDDO
C
      DO I=1,MARK1
         AH1(I,1)=1.0D0
      ENDDO
      DO I=MARK1+1,MH
         AH1(I,2)=1.0D0
      ENDDO
C
C Contrain for H at Tint
C
      BHk(1,1)=1.0D0
      BHk(1,2)=-1.0D0
      DH1(1)=PHI_NIST_H(Tint,X1(6:10))
     &      -PHI_NIST_H(Tint,X1(1:5))
C
C Contrain for H at 298.15K
C
      BHk(2,1)=1.0D0
      BHk(2,2)=0.D0
      DH1(2)=H_298 + DLTH - H_298 - PHI_NIST_H(298.15D0,X1(1:5))
C
      DO I=1,NH
         DO J=1,LDBH
            BH1(J,I)=BHk(J,I)
         ENDDO
      ENDDO
C
      call DGGLSE(MH,NH,PH,AH1,LDAH,BH1,LDBH,CH2,DH1,XH1,WORKH,
     &              LWORKH,INFO)
!     WRITE(*,*) 'H_CONST1=',XH1(1)
!     WRITE(*,*) 'H_CONST2=',XH1(2)
C
      A_H1=XH1(1)
      A_H2=XH1(2)
C
C For entropy, NEED define 2 coeffs, N=2
      DO I=1,NS
         XS1(I)=0.D0
      ENDDO
      DO I=1,NS
         DO J=1,LDAS
           AS1(J,I)=0.D0
         ENDDO
         DO J=1,LDBS
           BS1(J,I)=0.D0
           BSk(J,I)=0.D0
         ENDDO
      ENDDO
      DO I=1,LWORKS
         WORKS(I)=0.D0
      ENDDO
C
      DO I=1,MARK1
         CS2(I)=CS(I) - PHI_NIST_S(TEMP(I),X1(1:5))
      ENDDO
      DO I=MARK1+1,MS
         CS2(I)=CS(I) - PHI_NIST_S(TEMP(I),X1(6:10))
      ENDDO
      DO I=1,MARK1
         AS1(I,1)=1.0D0
      ENDDO
      DO I=MARK1+1,MS
         AS1(I,2)=1.0D0
      ENDDO
C
C Contrain for S at Tint
C
      BSk(1,1)=1.0D0
      BSk(1,2)=-1.0D0
      DS1(1)=PHI_NIST_S(Tint,X1(6:10))
     &      -PHI_NIST_S(Tint,X1(1:5))
C
C Contrain for S at 298.15K
C
      BSk(2,1)=1.0D0
      BSk(2,2)=0.0D0
!     DS1(2)=CS(1) - PHI_NIST_S(298.D0,X1(1:5))
      DS1(2)=S_298 - PHI_NIST_S(298.15D0,X1(1:5))
C
      DO I=1,NS
         DO J=1,LDBS
            BS1(J,I)=BSk(J,I)
         ENDDO
      ENDDO
C
      call DGGLSE(MS,NS,PS,AS1,LDAS,BS1,LDBS,CS2,DS1,XS1,WORKS,
     &              LWORKS,INFO)
!     WRITE(*,*) XS1
!     WRITE(*,*) 'S_CONST1=',XS1(1)
!     WRITE(*,*) 'S_CONST2=',XS1(2)
C
      A_S1=XS1(1)
      A_S2=XS1(2)
C
C
        DO I=1,5
           THERM1(I)=X1(I)
        ENDDO
        THERM1(5+1)=A_H1
        THERM1(5+2)=A_S1
        DO I=1,5
           THERM2(I)=X1(I+5)
        ENDDO
        THERM2(5+1)=A_H2
        THERM2(5+2)=A_S2
C
!     WRITE(*,*) THERM1(1:7)
!     WRITE(*,*) THERM2(1:7)
C
      RETURN
C
      END
C
      DOUBLE PRECISION FUNCTION PHI_NIST1(TEM)
C
      implicit none
C
      DOUBLE PRECISION TEM,ROL,H_298,S_298,DLTH
      COMMON /GAS/ROL,H_298,S_298,DLTH
C
      PHI_NIST1=1.D0
C
      RETURN
C
      END
C
      DOUBLE PRECISION FUNCTION PHI_NIST2(TEM)
C
      implicit none
C
      DOUBLE PRECISION TEM,ROL,H_298,S_298,DLTH
      COMMON /GAS/ROL,H_298,S_298,DLTH
C
      PHI_NIST2=TEM/1000.D0
C
      RETURN
C
      END
C
      DOUBLE PRECISION FUNCTION PHI_NIST3(TEM)
C
      implicit none
C
      DOUBLE PRECISION TEM,ROL,H_298,S_298,DLTH
      COMMON /GAS/ROL,H_298,S_298,DLTH
C
      PHI_NIST3=(TEM/1000.D0)*(TEM/1000.D0)
C
      RETURN
C
      END
C
C
      DOUBLE PRECISION FUNCTION PHI_NIST4(TEM)
C
      implicit none
C
      DOUBLE PRECISION TEM,ROL,H_298,S_298,DLTH
      COMMON /GAS/ROL,H_298,S_298,DLTH
C
      PHI_NIST4=(TEM/1000.D0)*(TEM/1000.D0)*(TEM/1000.D0)
C
      RETURN
C
      END
C
      DOUBLE PRECISION FUNCTION PHI_NIST5(TEM)
C
      implicit none
C
      DOUBLE PRECISION TEM,ROL,H_298,S_298,DLTH
      COMMON /GAS/ROL,H_298,S_298,DLTH
C
      PHI_NIST5=1.0D0/(TEM/1000.D0)/(TEM/1000.D0)
C
      RETURN
C
      END
C
      DOUBLE PRECISION FUNCTION PHI_NIST_DCP1(TEM)
C
      implicit none
C
      DOUBLE PRECISION TEM,ROL,H_298,S_298,DLTH
      COMMON /GAS/ROL,H_298,S_298,DLTH
C
      PHI_NIST_DCP1=0.D0
C
      RETURN
C
      END
C
      DOUBLE PRECISION FUNCTION PHI_NIST_DCP2(TEM)
C
      implicit none
C
      DOUBLE PRECISION TEM,ROL,H_298,S_298,DLTH
      COMMON /GAS/ROL,H_298,S_298,DLTH
C
      PHI_NIST_DCP2=1.D0
C
      RETURN
C
      END
C
      DOUBLE PRECISION FUNCTION PHI_NIST_DCP3(TEM)
C
      implicit none
C
      DOUBLE PRECISION TEM,ROL,H_298,S_298,DLTH
      COMMON /GAS/ROL,H_298,S_298,DLTH
C
      PHI_NIST_DCP3=2.D0*(TEM/1000.D0)
C
      RETURN
C
      END
C
      DOUBLE PRECISION FUNCTION PHI_NIST_DCP4(TEM)
C
      implicit none
C
      DOUBLE PRECISION TEM,ROL,H_298,S_298,DLTH
      COMMON /GAS/ROL,H_298,S_298,DLTH
C
      PHI_NIST_DCP4=3.D0*(TEM/1000.D0)*(TEM/1000.D0)
C
      RETURN
C
      END
C
      DOUBLE PRECISION FUNCTION PHI_NIST_DCP5(TEM)
C
      implicit none
C
      DOUBLE PRECISION TEM,ROL,H_298,S_298,DLTH
      COMMON /GAS/ROL,H_298,S_298,DLTH
C
      PHI_NIST_DCP5=-2.D0/(TEM/1000.D0)/(TEM/1000.D0)/(TEM/1000.D0)
C
      RETURN
C
      END
C
      DOUBLE PRECISION FUNCTION PHI_NIST_D2CP1(TEM)
C
      implicit none
C
      DOUBLE PRECISION TEM,ROL,H_298,S_298,DLTH
      COMMON /GAS/ROL,H_298,S_298,DLTH
C
      PHI_NIST_D2CP1=0.D0
C
      RETURN
C
      END
C
      DOUBLE PRECISION FUNCTION PHI_NIST_D2CP2(TEM)
C
      implicit none
C
      DOUBLE PRECISION TEM,ROL,H_298,S_298,DLTH
      COMMON /GAS/ROL,H_298,S_298,DLTH
C
      PHI_NIST_D2CP2=0.D0
C
      RETURN
C
      END
C
      DOUBLE PRECISION FUNCTION PHI_NIST_D2CP3(TEM)
C
      implicit none
C
      DOUBLE PRECISION TEM,ROL,H_298,S_298,DLTH
      COMMON /GAS/ROL,H_298,S_298,DLTH
C
      PHI_NIST_D2CP3=2.D0
C
      RETURN
C
      END
C
      DOUBLE PRECISION FUNCTION PHI_NIST_D2CP4(TEM)
C
      implicit none
C
      DOUBLE PRECISION TEM,ROL,H_298,S_298,DLTH
      COMMON /GAS/ROL,H_298,S_298,DLTH
C
      PHI_NIST_D2CP4=6.D0*(TEM/1000.D0)
C
      RETURN
C
      END
C
      DOUBLE PRECISION FUNCTION PHI_NIST_D2CP5(TEM)
C
      implicit none
C
      DOUBLE PRECISION TEM,ROL,H_298,S_298,DLTH
      COMMON /GAS/ROL,H_298,S_298,DLTH
C
      PHI_NIST_D2CP5=6.D0/(TEM/1000.D0)/(TEM/1000.D0)
     &              /(TEM/1000.D0)/(TEM/1000.D0)
C
      RETURN
C
      END
C
      DOUBLE PRECISION FUNCTION PHI_NIST_H(TEM,alpha)
C
      implicit none
C
      DOUBLE PRECISION TEM,ROL,alpha,H_298,S_298,DLTH,T1000
      DIMENSION alpha(5)
      COMMON /GAS/ROL,H_298,S_298,DLTH
C
      T1000=TEM/1000.D0
      PHI_NIST_H=alpha(1)*T1000+alpha(2)*T1000*T1000/2.D0
     &          +alpha(3)*T1000*T1000*T1000/3.D0
     &          +alpha(4)*T1000*T1000*T1000*T1000/4.D0
     &          -alpha(5)/T1000
!      WRITE(*,*) 'PHI_NIST_H=',PHI_NIST_H
      WRITE(*,*) 'H calculated by Shomate is too small!'
      WRITE(*,*) 'warning! warning! warning!'
      PAUSE
C
      RETURN
C
      END
C
      DOUBLE PRECISION FUNCTION PHI_NIST_S(TEM,alpha)
C
      implicit none
C
      DOUBLE PRECISION TEM,ROL,alpha,H_298,S_298,DLTH,T1000
      DIMENSION alpha(5)
      COMMON /GAS/ROL,H_298,S_298,DLTH
C
!     WRITE(*,*) ALPHA(1:5)
!     WRITE(*,*) TEM
      T1000=TEM/1000.D0
      PHI_NIST_S=alpha(1)*DLOG(T1000)+alpha(2)*T1000
     &          +alpha(3)*T1000*T1000/2.D0
     &          +alpha(4)*T1000*T1000*T1000/3.D0
     &          -alpha(5)/T1000/T1000/2.D0
C
      RETURN
C
      END

ext: f

max_issues_repo_name: bslakman/RMG-Java

max_forks_repo_name: bslakman/RMG-Java


text:
c $Header: /data/petsun4/data1/src_solaris/imgblur_4dfp/RCS/fimgblur.f,v 1.3 2008/01/01 22:39:55 avi Exp $
c $Log: fimgblur.f,v $
c Revision 1.3  2008/01/01  22:39:55  avi
c gcc v4 compliant
c
c Revision 1.2  1999/03/10  02:23:51  avi
c subroutine imgblur3d
c
c Revision 1.1  1998/10/02  19:24:42  mcavoy
c Initial revision
c
c Revision 1.1  1997/10/30  15:38:21  tscull
c Initial revision
c
c Revision 1.1  1997/01/16  04:47:58  avi
c Initial revision
c
      subroutine fimgblur_rcs
      write(*,"('$Id: fimgblur.f,v 1.3 2008/01/01 22:39:55 avi Exp $')")
      return
      end

      subroutine imgblur3d(fwhm,alpha,voxsiz,img0,imgdim,img1)
      real*4 img0(1),img1(1),voxsiz(3)
      integer*4 imgdim(3)

      call fimgblur_rcs
      call imgblurx(fwhm,alpha,voxsiz(1),img0,imgdim(1),imgdim(2),imgdim(3),img1)
      call imgblury(fwhm,alpha,voxsiz(2),img1,imgdim(1),imgdim(2),imgdim(3),img0)
      call imgblurz(fwhm,alpha,voxsiz(3),img0,imgdim(1),imgdim(2),imgdim(3),img1)
      return
      end

      subroutine imgblurx(fwhm,alpha,voxsiz,imag,nx,ny,nz,imgb)
      real*4 imag(nx,ny,nz),imgb(nx,ny,nz)
      parameter (nkmax=256)
      real*4 h(-nkmax:nkmax)
      logical*4 ldebug/.false./

      if(fwhm.le.0..or.alpha.le.0.)then
        write(*,"('imgblurx: fwhm and alpha must be greater than 0')")
        return
      endif
      sigma=0.5*fwhm/sqrt(2.*alog(2.))
      nk=nint(sigma*sqrt(-2.*alog(alpha))/abs(voxsiz))
      if(nk.gt.nkmax)then
        write(*,"('imgblurx: nk reduced to',i5)")nkmax
        nk=nkmax
      endif
      if(ldebug) write(*,"('imgblurx: nk=',i5)")nk

      s=0.
      do 5 k=-nk,nk
      h(k)=exp(-0.5*(float(k)*voxsiz/sigma)**2)
    5 s=s+h(k)
      do 6 k=-nk,nk
    6 h(k)=h(k)/s

      do 1 iz=1,nz
      do 1 iy=1,ny
      do 1 ix=1,nx
      s=0.
      do 2 k=-nk,nk
      kx=min0(max0(ix+k,1),nx)
    2 s=s+h(k)*imag(kx,iy,iz)
    1 imgb(ix,iy,iz)=s

      return
      end

      subroutine imgblury(fwhm,alpha,voxsiz,imag,nx,ny,nz,imgb)
      real*4 imag(nx,ny,nz),imgb(nx,ny,nz)
      parameter (nkmax=256)
      real*4 h(-nkmax:nkmax)
      logical*4 ldebug/.false./

      if(fwhm.le.0..or.alpha.le.0.)then
        write(*,"('imgblurx: fwhm and alpha must be greater than 0')")
        return
      endif
      sigma=0.5*fwhm/sqrt(2.*alog(2.))
      nk=nint(sigma*sqrt(-2.*alog(alpha))/abs(voxsiz))
      if(nk.gt.nkmax)then
        write(*,"('imgblury: nk reduced to',i5)")nkmax
        nk=nkmax
      endif
      if(ldebug) write(*,"('imgblury: nk=',i5)")nk

      s=0.
      do 5 k=-nk,nk
      h(k)=exp(-0.5*(float(k)*voxsiz/sigma)**2)
    5 s=s+h(k)
      do 6 k=-nk,nk
    6 h(k)=h(k)/s

      do 1 iz=1,nz
      do 1 ix=1,nx
      do 1 iy=1,ny
      s=0.
      do 2 k=-nk,nk
      ky=min0(max0(iy+k,1),ny)
    2 s=s+h(k)*imag(ix,ky,iz)
    1 imgb(ix,iy,iz)=s

      return
      end

      subroutine imgblurz(fwhm,alpha,voxsiz,imag,nx,ny,nz,imgb)
      real*4 imag(nx,ny,nz),imgb(nx,ny,nz)
      parameter (nkmax=256)
      real*4 h(-nkmax:nkmax)
      logical*4 ldebug/.false./

      if(fwhm.le.0..or.alpha.le.0.)then
        write(*,"('imgblurz: fwhm and alpha must be greater than 0')")
        return
      endif
      sigma=0.5*fwhm/sqrt(2.*alog(2.))
      nk=nint(sigma*sqrt(-2.*alog(alpha))/abs(voxsiz))
      if(nk.gt.nkmax)then
        write(*,"('imgblurx: nk reduced to',i5)")nkmax
        nk=nkmax
      endif
      if(ldebug) write(*,"('imgblurz: nk=',i5)")nk

      s=0.
      do 5 k=-nk,nk
      h(k)=exp(-0.5*(float(k)*voxsiz/sigma)**2)
    5 s=s+h(k)
      do 6 k=-nk,nk
    6 h(k)=h(k)/s

      do 1 iy=1,ny
      do 1 ix=1,nx
      do 1 iz=1,nz
      s=0.
      do 2 k=-nk,nk
      kz=min0(max0(iz+k,1),nz)
    2 s=s+h(k)*imag(ix,iy,kz)
    1 imgb(ix,iy,iz)=s

      return
      end

max_stars_repo_name: llevitis/PUP

max_issues_repo_name: llevitis/PUP

max_forks_repo_name: llevitis/PUP


text:
!
!     CalculiX - A 3-dimensional finite element program
!              Copyright (C) 1998-2019 Guido Dhondt
!
!     This program is free software; you can redistribute it and/or
!     modify it under the terms of the GNU General Public License as
!     published by the Free Software Foundation(version 2);
!
!
!     This program is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!     GNU General Public License for more details.
!
!     You should have received a copy of the GNU General Public License
!     along with this program; if not, write to the Free Software
!     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
!
      subroutine changesolidsections(inpc,textpart,set,istartset,
     &  iendset,ialset,nset,ielmat,matname,nmat,ielorien,orname,norien,
     &  lakon,thicke,kon,ipkon,irstrt,istep,istat,n,iline,ipol,inl,
     &  ipoinp,inp,cs,mcs,iaxial,ipoinpc,mi,nelcon,ier)
!
!     reading the input deck: *CHANGE SOLID SECTION
!
      implicit none
!
      character*1 inpc(*)
      character*8 lakon(*)
      character*80 matname(*),orname(*),material,orientation
      character*81 set(*),elset
      character*132 textpart(16)
!
      integer mi(*),istartset(*),iendset(*),ialset(*),ielmat(mi(3),*),
     &  ielorien(mi(3),*),kon(*),ipkon(*),indexe,irstrt(*),nset,nmat,
     &  norien,nelcon(2,*),ier,
     &  istep,istat,n,key,i,j,k,l,imaterial,iorientation,ipos,
     &  iline,ipol,inl,ipoinp(2,*),inp(3,*),mcs,iaxial,ipoinpc(0:*)
!
      real*8 thicke(mi(3),*),thickness,pi,cs(17,*)
!
      if(istep.eq.0) then
         write(*,*) '*ERROR reading *CHANGE SOLID SECTION:'
         write(*,*) '       *CHANGE SOLID SECTION should'
         write(*,*) '       be placed within a step definition'
         ier=1
         return
      endif
!
      pi=4.d0*datan(1.d0)
!
      orientation='
     &                           '
      elset='
     &                      '
      ipos=0
!
      do i=2,n
         if(textpart(i)(1:9).eq.'MATERIAL=') then
            material=textpart(i)(10:89)
         elseif(textpart(i)(1:12).eq.'ORIENTATION=') then
            orientation=textpart(i)(13:92)
         elseif(textpart(i)(1:6).eq.'ELSET=') then
            elset=textpart(i)(7:86)
            elset(81:81)=' '
            ipos=index(elset,' ')
            elset(ipos:ipos)='E'
         else
            write(*,*) '*WARNING reading *CHANGE SOLID SECTION:'
            write(*,*) '         parameter not recognized:'
            write(*,*) '         ',
     &                 textpart(i)(1:index(textpart(i),' ')-1)
            call inputwarning(inpc,ipoinpc,iline,
     &"*CHANGE SOLID SECTIONQ%")
         endif
      enddo
!
!     check for the existence of the set,the material and orientation
!
      do i=1,nmat
         if(matname(i).eq.material) exit
      enddo
      if(i.gt.nmat) then
         do i=1,nmat
            if(matname(i)(1:11).eq.'ANISO_CREEP') then
               if(matname(i)(12:20).eq.material(1:9)) exit
            elseif(matname(i)(1:10).eq.'ANISO_PLAS') then
               if(matname(i)(11:20).eq.material(1:10)) exit
            endif
         enddo
      endif
      if(i.gt.nmat) then
         write(*,*)
     &     '*ERROR reading *CHANGE SOLID SECTION: nonexistent material'
         write(*,*) '  '
         call inputerror(inpc,ipoinpc,iline,
     &        "*CHANGE SOLID SECTION%",ier)
         return
      endif
      imaterial=i
!
      if(orientation.eq.'                    ') then
         iorientation=0
c      elseif(nelcon(1,i).eq.2) then
c         write(*,*)
c     &      '*INFO reading *CHANGE SOLID SECTION: an orientation'
c         write(*,*) '      is for isotropic materials irrelevant'
c         call inputinfo(inpc,ipoinpc,iline,
c     &"*CHANGE SOLID SECTION%")
c         iorientation=0
      else
         do i=1,norien
            if(orname(i).eq.orientation) exit
         enddo
         if(i.gt.norien) then
            write(*,*) '*ERROR reading *CHANGE SOLID SECTION:'
            write(*,*) '       nonexistent orientation'
            write(*,*) '  '
            call inputerror(inpc,ipoinpc,iline,
     &           "*CHANGE SOLID SECTION%",ier)
            return
         endif
         iorientation=i
      endif
!
      if(ipos.eq.0) then
         write(*,*)
     &    '*ERROR reading *CHANGE SOLID SECTION: no element set ',elset
         write(*,*) '       was been defined. '
         call inputerror(inpc,ipoinpc,iline,
     &        "*CHANGE SOLID SECTION%",ier)
         return
      endif
      do i=1,nset
         if(set(i).eq.elset) exit
      enddo
      if(i.gt.nset) then
         elset(ipos:ipos)=' '
         write(*,*) '*ERROR reading *CHANGE SOLID SECTION:'
         write(*,*) '       element set ',elset
         write(*,*) '  has not yet been defined. '
         call inputerror(inpc,ipoinpc,iline,
     &        "*CHANGE SOLID SECTION%",ier)
         return
      endif
!
!     assigning the elements of the set the appropriate material
!     and orientation number
!
      do j=istartset(i),iendset(i)
         if(ialset(j).gt.0) then
            if((lakon(ialset(j))(1:1).eq.'B').or.
     &         (lakon(ialset(j))(1:1).eq.'S')) then
               write(*,*) '*ERROR reading *CHANGE SOLID SECTION:'
               write(*,*) '       *CHANGE SOLID SECTION can'
               write(*,*) '       not be used for beam or shell elements
     &'
               write(*,*) '       Faulty element: ',ialset(j)
               ier=1
               return
            endif
            ielmat(1,ialset(j))=imaterial
            if(norien.gt.0) ielorien(1,ialset(j))=iorientation
         else
            k=ialset(j-2)
            do
               k=k-ialset(j)
               if(k.ge.ialset(j-1)) exit
               if((lakon(k)(1:1).eq.'B').or.
     &              (lakon(k)(1:1).eq.'S')) then
                  write(*,*) '*ERROR reading *CHANGE SOLID SECTION:'
                  write(*,*) '       *CHANGE SOLID SECTION can'
                  write(*,*) '       not be used for beam or shell eleme
     &nts'
                  write(*,*) '       Faulty element: ',k
                  ier=1
                  return
               endif
               ielmat(1,k)=imaterial
               if(norien.gt.0) ielorien(1,k)=iorientation
            enddo
         endif
      enddo
!
      call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl,
     &     ipoinp,inp,ipoinpc)
!
!     assigning a thickness to plane stress elements and an angle to
!     axisymmetric elements
!
      if(key.eq.0) then
         write(*,*) '*ERROR reading *CHANGE SOLID SECTION'
         write(*,*) '       no second line allowed'
         call inputerror(inpc,ipoinpc,iline,
     &        "*CHANGE SOLID SECTION%",ier)
         return
      endif
!
      return
      end

max_stars_repo_name: alleindrach/calculix-desktop

max_issues_repo_name: alleindrach/calculix-desktop

max_forks_repo_name: alleindrach/calculix-desktop


! { dg-do run }
!
! Test the behavior of lbound, ubound of shape with assumed rank arguments
! in an array context (without DIM argument).
!

program test

  integer              :: a(2:4,-2:5)
  integer, allocatable :: b(:,:)
  integer, allocatable :: c(:,:)
  integer, pointer     :: d(:,:)
  character(52)        :: buffer

  b = foo(a)
  !print *,b(:,1)
  if (any(b(:,1) /= [11, 101])) call abort
  buffer = 'ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz'
  write(buffer,*) b(:,1)
  if (buffer /= '          11         101') call abort

  !print *,b(:,2)
  if (any(b(:,2) /= [3, 8])) call abort
  buffer = 'ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz'
  write(buffer,*) b(:,2)
  if (buffer /= '           3           8') call abort

  !print *,b(:,3)
  if (any(b(:,3) /= [13, 108])) call abort
  buffer = 'ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz'
  write(buffer,*) b(:,3)
  if (buffer /= '          13         108') call abort

  allocate(c(1:2,-3:6))
  b = bar(c)
  !print *,b(:,1)
  if (any(b(:,1) /= [11, 97])) call abort
  buffer = 'ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz'
  write(buffer,*) b(:,1)
  if (buffer /= '          11          97') call abort

  !print *,b(:,2)
  if (any(b(:,2) /= [12, 106])) call abort
  buffer = 'ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz'
  write(buffer,*) b(:,2)
  if (buffer /= '          12         106') call abort

  !print *,b(:,3)
  if (any(b(:,3) /= [2, 10])) call abort
  buffer = 'ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz'
  write(buffer,*) b(:,3)
  if (buffer /= '           2          10') call abort

  allocate(d(3:5,-1:10))
  b = baz(d)
  !print *,b(:,1)
  if (any(b(:,1) /= [3, -1])) call abort
  buffer = 'ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz'
  write(buffer,*) b(:,1)
  if (buffer /= '           3          -1') call abort

  !print *,b(:,2)
  if (any(b(:,2) /= [15, 110])) call abort
  buffer = 'ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz'
  write(buffer,*) b(:,2)
  if (buffer /= '          15         110') call abort

  !print *,b(:,3)
  if (any(b(:,3) /= [13, 112])) call abort
  buffer = 'ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz'
  write(buffer,*) b(:,3)
  if (buffer /= '          13         112') call abort

contains
  function foo(arg) result(res)
    integer :: arg(..)
    integer, allocatable :: res(:,:)

    allocate(res(rank(arg), 3))

    res(:,1) = lbound(arg) + (/ 10, 100 /)
    res(:,2) = ubound(arg)
    res(:,3) = (/ 10, 100 /) + shape(arg)

  end function foo
  function bar(arg) result(res)
    integer, allocatable :: arg(..)
    integer, allocatable :: res(:,:)

    allocate(res(-1:rank(arg)-2, 3))

    res(:,1) = lbound(arg) + (/ 10, 100 /)
    res(:,2) = (/ 10, 100 /) + ubound(arg)
    res(:,3) = shape(arg)

  end function bar
  function baz(arg) result(res)
    integer, pointer     :: arg(..)
    integer, allocatable :: res(:,:)

    allocate(res(2:rank(arg)+1, 3))

    res(:,1) = lbound(arg)
    res(:,2) = (/ 10, 100 /) + ubound(arg)
    res(:,3) = shape(arg) + (/ 10, 100 /)

  end function baz
end program test