module cantera_thermo

  use fct
  use cantera_xml

  type phase_t
    integer :: thermo_id
    integer :: kin_id
    integer :: tran_id
    integer :: err
    integer :: nel
    integer :: nsp
    integer :: nrxn
  end type phase_t

! these definitions are for use with the equilibrate function.
  integer, parameter :: TV = 100
  integer, parameter :: HP = 101
  integer, parameter :: SP = 102
  integer, parameter :: PV = 103
  integer, parameter :: TP = 104
  integer, parameter :: UV = 105
  integer, parameter :: SV = 107

  integer, parameter :: VT = -100
  integer, parameter :: PH = -101
  integer, parameter :: PS = -102
  integer, parameter :: VP = -103
  integer, parameter :: PT = -104
  integer, parameter :: VU = -105
  integer, parameter :: VS = -107

contains

    type(phase_t) function newThermoPhase(xml_phase, index)
      implicit none
      type(XML_Node), intent(inout), optional :: xml_phase
      integer, intent(in), optional :: index
      type(phase_t) :: self
      if (present(index)) then
         self%thermo_id = index
         self%err = 0
      else if (present(xml_phase)) then
         self%thermo_id = newthermofromxml(xml_phase%xml_id)
         self%nel = phase_nelements(self%thermo_id)
         self%nsp = phase_nspecies(self%thermo_id)
         self%nrxn = 0
         self%err = 0
      else
         call cantera_error('newThermoPhase','xml_phase or id must be specified.')
      end if
      newThermoPhase = self
    end function newThermoPhase

    integer function ctthermo_nElements(self)
      implicit none
      type(phase_t), intent(inout) :: self
      ctthermo_nelements = phase_nelements(self%thermo_id)
    end function ctthermo_nElements

    integer function ctthermo_nSpecies(self)
      implicit none
      type(phase_t), intent(inout) :: self
      ctthermo_nspecies = phase_nspecies(self%thermo_id)
    end function ctthermo_nSpecies

    double precision function ctthermo_temperature(self)
      implicit none
      type(phase_t), intent(inout) :: self
      ctthermo_temperature = phase_temperature(self%thermo_id)
    end function ctthermo_temperature

    subroutine ctthermo_setTemperature(self, t)
      implicit none
      type(phase_t), intent(inout) :: self
      double precision, intent(in) :: t
      self%err = phase_settemperature(self%thermo_id, t)
    end subroutine ctthermo_setTemperature

    double precision function ctthermo_density(self)
      implicit none
      type(phase_t), intent(inout) :: self
      ctthermo_density = phase_density(self%thermo_id)
    end function ctthermo_density

    subroutine ctthermo_setDensity(self, rho)
      implicit none
      type(phase_t), intent(inout) :: self
      double precision, intent(in) :: rho
      self%err = phase_setdensity(self%thermo_id, rho)
    end subroutine ctthermo_setDensity

    double precision function ctthermo_molarDensity(self)
      implicit none
      type(phase_t), intent(inout) :: self
      ctthermo_molardensity = phase_molardensity(self%thermo_id)
    end function ctthermo_molarDensity

    double precision function ctthermo_meanMolecularWeight(self)
      implicit none
      type(phase_t), intent(inout) :: self
      ctthermo_meanmolecularweight = phase_meanmolecularweight(self%thermo_id)
    end function ctthermo_meanMolecularWeight

    integer function ctthermo_elementIndex(self, nm)
      implicit none
      type(phase_t), intent(inout) :: self
      character*(*), intent(in) :: nm
      ctthermo_elementindex = phase_elementindex(self%thermo_id, nm)
    end function ctthermo_elementIndex

    integer function ctthermo_speciesIndex(self, nm)
      implicit none
      type(phase_t), intent(inout) :: self
      character*(*), intent(in) :: nm
      ctthermo_speciesindex = phase_speciesindex(self%thermo_id, nm)
    end function ctthermo_speciesIndex

    subroutine ctthermo_getMoleFractions(self, x)
      implicit none
      type(phase_t), intent(inout) :: self
      double precision, intent(out) :: x(self%nsp)
      self%err = phase_getmolefractions(self%thermo_id, x)
    end subroutine ctthermo_getMoleFractions

    double precision function ctthermo_moleFraction(self, k)
      implicit none
      type(phase_t), intent(inout) :: self
      integer, intent(in) :: k
      ctthermo_molefraction = phase_molefraction(self%thermo_id, k)
    end function ctthermo_moleFraction

    subroutine ctthermo_getMassFractions(self, y)
      implicit none
      type(phase_t), intent(inout) :: self
      double precision, intent(out) :: y(self%nsp)
      self%err = phase_getmassfractions(self%thermo_id, y)
    end subroutine ctthermo_getMassFractions

    double precision function ctthermo_massFraction(self, k)
      implicit none
      type(phase_t), intent(inout) :: self
      integer, intent(in) :: k
      ctthermo_massfraction = phase_massfraction(self%thermo_id, k)
    end function ctthermo_massFraction

    subroutine ctthermo_setMoleFractions(self, x, norm)
      implicit none
      type(phase_t), intent(inout) :: self
      double precision, intent(in) :: x(self%nsp)
      integer, intent(in), optional :: norm
      integer :: n
      if (present(norm)) then
         n = norm
      else
         n = 1
      end if
      self%err = phase_setmolefractions(self%thermo_id, x, n)
    end subroutine ctthermo_setmolefractions

    subroutine ctthermo_setMoleFractionsByName(self, x)
      implicit none
      type(phase_t), intent(inout) :: self
      character*(*), intent(in) :: x
      self%err = phase_setmolefractionsbyname(self%thermo_id, x)
    end subroutine ctthermo_setmolefractionsbyname

    subroutine ctthermo_setMassFractions(self, y, norm)
      implicit none
      type(phase_t), intent(inout) :: self
      double precision, intent(in) :: y(self%nsp)
      integer, intent(in), optional :: norm
      integer :: n
      if (present(norm)) then
         n = norm
      else
         n = 1
      end if
      self%err = phase_setmassfractions(self%thermo_id, y, n)
    end subroutine ctthermo_setmassfractions

    subroutine ctthermo_setMassFractionsByName(self, y)
      implicit none
      type(phase_t), intent(inout) :: self
      character*(*), intent(in) :: y
      self%err = phase_setmassfractionsbyname(self%thermo_id, y)
    end subroutine ctthermo_setmassfractionsbyname

    subroutine ctthermo_getAtomicWeights(self, atw)
      implicit none
      type(phase_t), intent(inout) :: self
      double precision, intent(out) :: atw(self%nel)
      self%err = phase_getatomicweights(self%thermo_id, atw)
    end subroutine ctthermo_getatomicweights

    subroutine ctthermo_getMolecularWeights(self, mw)
      implicit none
      type(phase_t), intent(inout) :: self
      double precision, intent(out) :: mw(self%nsp)
      self%err = phase_getmolecularweights(self%thermo_id, mw)
    end subroutine ctthermo_getmolecularweights

    subroutine ctthermo_getSpeciesName(self, k, nm)
      implicit none
      type(phase_t), intent(inout) :: self
      integer, intent(in) :: k
      character*(*), intent(out) :: nm
      self%err = phase_getspeciesname(self%thermo_id, k, nm)
    end subroutine ctthermo_getspeciesname

    subroutine ctthermo_getElementName(self, m, nm)
      implicit none
      type(phase_t), intent(inout) :: self
      integer, intent(in) :: m
      character*(*), intent(out) :: nm
      self%err = phase_getelementname(self%thermo_id, m, nm)
    end subroutine ctthermo_getelementname

    double precision function ctthermo_nAtoms(self, k, m)
      implicit none
      type(phase_t), intent(inout) :: self
      integer, intent(in) :: k
      integer, intent(in) :: m
      ctthermo_natoms = phase_natoms(self%thermo_id, k, m)
    end function ctthermo_natoms

    integer function ctthermo_eosType(self)
      implicit none
      type(phase_t), intent(inout) :: self
      ctthermo_eostype = th_eostype(self%thermo_id)
    end function ctthermo_eostype

    double precision function ctthermo_enthalpy_mole(self)
      implicit none
      type(phase_t), intent(inout) :: self
      ctthermo_enthalpy_mole = th_enthalpy_mole(self%thermo_id)
    end function ctthermo_enthalpy_mole

    double precision function ctthermo_intEnergy_mole(self)
      implicit none
      type(phase_t), intent(inout) :: self
      ctthermo_intenergy_mole = th_intenergy_mole(self%thermo_id)
    end function ctthermo_intenergy_mole

    double precision function ctthermo_entropy_mole(self)
      implicit none
      type(phase_t), intent(inout) :: self
      ctthermo_entropy_mole = th_entropy_mole(self%thermo_id)
    end function ctthermo_entropy_mole

    double precision function ctthermo_gibbs_mole(self)
      implicit none
      type(phase_t), intent(inout) :: self
      ctthermo_gibbs_mole = th_gibbs_mole(self%thermo_id)
    end function ctthermo_gibbs_mole

    double precision function ctthermo_cp_mole(self)
      implicit none
      type(phase_t), intent(inout) :: self
      ctthermo_cp_mole = th_cp_mole(self%thermo_id)
    end function ctthermo_cp_mole

    double precision function ctthermo_cv_mole(self)
      implicit none
      type(phase_t), intent(inout) :: self
      ctthermo_cv_mole = th_cv_mole(self%thermo_id)
    end function ctthermo_cv_mole

    double precision function ctthermo_pressure(self)
      implicit none
      type(phase_t), intent(inout) :: self
      ctthermo_pressure = th_pressure(self%thermo_id)
    end function ctthermo_pressure

    double precision function ctthermo_enthalpy_mass(self)
      implicit none
      type(phase_t), intent(inout) :: self
      ctthermo_enthalpy_mass = th_enthalpy_mass(self%thermo_id)
    end function ctthermo_enthalpy_mass

    double precision function ctthermo_intEnergy_mass(self)
      implicit none
      type(phase_t), intent(inout) :: self
      ctthermo_intEnergy_mass = th_intenergy_mass(self%thermo_id)
    end function ctthermo_intEnergy_mass

    double precision function ctthermo_entropy_mass(self)
      implicit none
      type(phase_t), intent(inout) :: self
      ctthermo_entropy_mass = th_entropy_mass(self%thermo_id)
    end function ctthermo_entropy_mass

    double precision function ctthermo_gibbs_mass(self)
      implicit none
      type(phase_t), intent(inout) :: self
      ctthermo_gibbs_mass = th_gibbs_mass(self%thermo_id)
    end function ctthermo_gibbs_mass

    double precision function ctthermo_cp_mass(self)
      implicit none
      type(phase_t), intent(inout) :: self
      ctthermo_cp_mass = th_cp_mass(self%thermo_id)
    end function ctthermo_cp_mass

    double precision function ctthermo_cv_mass(self)
      implicit none
      type(phase_t), intent(inout) :: self
      ctthermo_cv_mass = th_cv_mass(self%thermo_id)
    end function ctthermo_cv_mass

    subroutine ctthermo_chemPotentials(self, mu)
      implicit none
      type(phase_t), intent(inout) :: self
      double precision, intent(out) :: mu(self%nsp)
      self%err = th_chempotentials(self%thermo_id, mu)
    end subroutine ctthermo_chempotentials

    subroutine ctthermo_setPressure(self, p)
      implicit none
      type(phase_t), intent(inout) :: self
      double precision, intent(in) :: p
      self%err = th_setpressure(self%thermo_id, p)
    end subroutine ctthermo_setpressure

    subroutine ctthermo_setState_TPX(self, t, p, x)
      implicit none
      type(phase_t), intent(inout) :: self
      double precision, intent(in) :: t
      double precision, intent(in) :: p
      double precision, intent(in) :: x(*)
      call ctthermo_setTemperature(self, t)
      call ctthermo_setMoleFractions(self, x)
      call ctthermo_setPressure(self, p)
    end subroutine ctthermo_setState_TPX

    subroutine ctstring_setState_TPX(self, t, p, x)
      implicit none
      type(phase_t), intent(inout) :: self
      double precision, intent(in) :: t
      double precision, intent(in) :: p
      character*(*), intent(in) :: x
      call ctthermo_setTemperature(self, t)
      call ctthermo_setMoleFractionsByName(self, x)
      call ctthermo_setPressure(self, p)
    end subroutine ctstring_setState_TPX

    subroutine ctthermo_setState_TRY(self, t, rho, y)
      implicit none
      type(phase_t), intent(inout) :: self
      double precision, intent(in) :: t
      double precision, intent(in) :: rho
      double precision, intent(in) :: y(*)
      call ctthermo_setTemperature(self, t)
      call ctthermo_setMassFractions(self, y)
      call ctthermo_setDensity(self, rho)
    end subroutine ctthermo_setState_TRY

    subroutine ctstring_setState_TRY(self, t, rho, y)
      implicit none
      type(phase_t), intent(inout) :: self
      double precision, intent(in) :: t
      double precision, intent(in) :: rho
      character*(*), intent(in) :: y
      call ctthermo_setTemperature(self, t)
      call ctthermo_setMassFractionsByName(self, y)
      call ctthermo_setDensity(self, rho)
    end subroutine ctstring_setState_TRY

    subroutine ctthermo_setState_HP(self, h, p)
      implicit none
      type(phase_t), intent(inout) :: self
      double precision, intent(in) :: h
      double precision, intent(in) :: p
      self%err = th_set_hp(self%thermo_id, h, p)
    end subroutine ctthermo_setstate_hp

    subroutine ctthermo_setState_UV(self, u, v)
      implicit none
      type(phase_t), intent(inout) :: self
      double precision, intent(in) :: u
      double precision, intent(in) :: v
      self%err = th_set_uv(self%thermo_id, u, v)
    end subroutine ctthermo_setstate_uv

    subroutine ctthermo_setState_SV(self, s, v)
      implicit none
      type(phase_t), intent(inout) :: self
      double precision, intent(in) :: s
      double precision, intent(in) :: v
      self%err = th_set_sv(self%thermo_id, s, v)
    end subroutine ctthermo_setstate_sv

    subroutine ctthermo_setState_SP(self, s, p)
      implicit none
      type(phase_t), intent(inout) :: self
      double precision, intent(in) :: s
      double precision, intent(in) :: p
      self%err = th_set_sp(self%thermo_id, s, p)
    end subroutine ctthermo_setstate_sp

    subroutine ctthermo_equilibrate(self, XY)
      implicit none
      type(phase_t), intent(inout) :: self
      integer, intent(in) :: XY
      self%err = th_equil(self%thermo_id, XY)
    end subroutine ctthermo_equilibrate

    double precision function ctthermo_refPressure(self)
      implicit none
      type(phase_t), intent(inout) :: self
      ctthermo_refpressure = th_refpressure(self%thermo_id)
    end function ctthermo_refpressure

    double precision function ctthermo_minTemp(self, k)
      implicit none
      type(phase_t), intent(inout) :: self
      integer, intent(in) :: k
      ctthermo_mintemp = th_mintemp(self%thermo_id, k)
    end function ctthermo_mintemp

    double precision function ctthermo_maxTemp(self, k)
      implicit none
      type(phase_t), intent(inout) :: self
      integer, intent(in) :: k
      ctthermo_maxtemp = th_maxtemp(self%thermo_id, k)
    end function ctthermo_maxtemp

    subroutine ctthermo_getEnthalpies_RT(self, h_rt)
      implicit none
      type(phase_t), intent(inout) :: self
      double precision, intent(out) :: h_rt(self%nsp)
      self%err = th_getenthalpies_rt(self%thermo_id, h_rt)
    end subroutine ctthermo_getenthalpies_rt

    subroutine ctthermo_getEntropies_R(self, s_r)
      implicit none
      type(phase_t), intent(inout) :: self
      double precision, intent(out) :: s_r(self%nsp)
      self%err = th_getentropies_r(self%thermo_id, s_r)
    end subroutine ctthermo_getentropies_r

    subroutine ctthermo_getCp_R(self, lenm, cp_r)
      implicit none
      type(phase_t), intent(inout) :: self
      integer, intent(out) :: lenm
      double precision, intent(out) :: cp_r(self%nsp)
      self%err = th_getcp_r(self%thermo_id, lenm, cp_r)
    end subroutine ctthermo_getcp_r

end module cantera_thermo


text:
!RUN: %f18 -fdebug-dump-symbols %s | FileCheck %s

program p
  ! CHECK: a size=4 offset=0: ObjectEntity type: LOGICAL(4)
  ! CHECK: b size=4 offset=4: ObjectEntity type: REAL(4)
  logical :: a = .false.
  real :: b = 9.73
  ! CHECK: a: AssocEntity type: REAL(4) expr:b
  ! CHECK: b: AssocEntity type: LOGICAL(4) expr:a
  associate (b => a, a => b)
    print*, a, b
  end associate
  print*, a, b
end

max_stars_repo_name: NSCoder/llvm-project

max_issues_repo_name: NSCoder/llvm-project

max_forks_repo_name: NSCoder/llvm-project


text:
!
!     CalculiX 3-dimensional finite element program
!     Copyright (C) 1998-2020 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 e_c3d_trhs(co,nk,konl,lakonl,p1,p2,omx,bodyfx,
     &     nbody,ff,nelem,nmethod,rhcon,nrhcon,
     &     ielmat,ntmat_,vold,vcon,nelemload,
     &     sideload,xload,nload,idist,dtime,matname,mi,
     &     ttime,time,istep,iinc,xloadold,reltimef,shcon,nshcon,cocon,
     &     ncocon,physcon,nelemface,sideface,nface,
     &     ipompc,nodempc,coefmpc,nmpc,ikmpc,ilmpc,compressible,
     &     vcontu,yy,turbulent,ipvar,var,ipvarf,varf,dt)
!
!     computation of the energy element matrix and rhs for the element with
!     element with the topology in konl: step 4
!
!     ff: rhs
!
      implicit none
!
      character*1 sideface(*)
      character*8 lakonl
      character*20 sideload(*)
      character*80 matname(*),amat
!
      integer konl(8),ifaceq(8,6),nelemload(2,*),nk,nbody,nelem,
     &     nload,idist,i,j,k,i1,i2,j1,ncocon(2,*),k1,node,nfield,mi(*),
     &     nmethod,ii,jj,id,ipointer,ig,kk,nrhcon(*),ielmat(mi(3),*),
     &     nshcon(*),flux,compressible,
     &     ntmat_,nope,nopes,imat,mint2d,mint3d,ifacet(6,4),nopev,
     &     ifacew(8,5),istep,iinc,layer,kspt,jltyp,iflag,nelemface(*),
     &     nface,igl,idf,ipompc(*),nodempc(3,*),nmpc,ikmpc(*),ilmpc(*),
     &     iscale,turbulent,ipvar(*),index,ipvarf(*),iemchange
!
      real*8 co(3,*),shp(4,8),xs2(3,7),dvi,dt(*),
     &     p1(3),p2(3),bodyf(3),bodyfx(3),ff(78),cond,enthalpy,
     &     bf(3),q(3),xsjmod,dtem(3),vkl(3,3),corio(3),sinktemp,
     &     rhcon(0:1,ntmat_,*),reltimef,t(3,3),tv(3),bfv,press,
     &     vel(3),div,shcon(0:3,ntmat_,*),pgauss(3),dxsj2,areaj,
     &     voldl(0:mi(2),8),xloadold(2,*),cocon(0:6,ntmat_,*),
     &     xl2(3,8),xsj2(3),shp2(7,8),vold(0:mi(2),*),xload(2,*),
     &     om,omx,xi,et,ze,const,xsj,field,physcon(*),tvn,
     &     temp,vcon(0:4,*),vconl(0:4,8),rho,xi3d,et3d,ze3d,
     &     weight,shpv(8),xlocal20(3,9,6),coefmpc(*),
     &     xlocal4(3,1,4),xlocal10(3,3,4),xlocal6(3,1,5),vl(0:mi(2),8),
     &     xlocal15(3,4,5),xlocal8(3,4,6),xlocal8r(3,1,6),omcor,
     &     shpvnithi(8),vcontu(2,*),vcontul(2,8),yy(*),yyl(8),
     &     y,xtuf,xkin,vort,un,unt,umt,f2,arg2,c1,c2,a1,var(*),varf(*)
!
      real*8 dtime,ttime,time,tvar(2),coords(3)
!
!
!
      include "gauss.f"
      include "xlocal.f"
!
      ifaceq=reshape((/4,3,2,1,11,10,9,12,
     &     5,6,7,8,13,14,15,16,
     &     1,2,6,5,9,18,13,17,
     &     2,3,7,6,10,19,14,18,
     &     3,4,8,7,11,20,15,19,
     &     4,1,5,8,12,17,16,20/),(/8,6/))
      ifacet=reshape((/1,3,2,7,6,5,
     &     1,2,4,5,9,8,
     &     2,3,4,6,10,9,
     &     1,4,3,8,10,7/),(/6,4/))
      ifacew=reshape((/1,3,2,9,8,7,0,0,
     &     4,5,6,10,11,12,0,0,
     &     1,2,5,4,7,14,10,13,
     &     2,3,6,5,8,15,11,14,
     &     4,6,3,1,12,15,9,13/),(/8,5/))
      iflag=3
      a1=0.31d0
!
      tvar(1)=time
      tvar(2)=ttime+dtime
!
      imat=ielmat(1,nelem)
      amat=matname(imat)
!
      if(lakonl(4:4).eq.'4') then
        nope=4
        mint3d=1
      elseif(lakonl(4:4).eq.'6') then
        nope=6
        mint3d=2
      elseif(lakonl(4:5).eq.'8R') then
        nope=8
        mint3d=1
      elseif(lakonl(4:4).eq.'8') then
        nope=8
        mint3d=8
      endif
!
!     initialisation for distributed forces
!
      do i=1,nope
        ff(i)=0.d0
      enddo
!
!     temperature, velocity and conservative variables
!     (rho*energy density, rho*velocity and rho)
!
      do i1=1,nope
        voldl(0,i1)=vold(0,konl(i1))
        voldl(4,i1)=vold(4,konl(i1))
        vconl(0,i1)=vcon(0,konl(i1))
      enddo
!
!     computation of the matrix: loop over the Gauss points
!
      index=ipvar(nelem)
      do kk=1,mint3d
        if(lakonl(4:5).eq.'8R') then
          weight=weight3d1(kk)
        elseif(lakonl(4:4).eq.'8') then
          weight=weight3d2(kk)
        elseif(lakonl(4:4).eq.'4') then
          weight=weight3d4(kk)
        elseif(lakonl(4:5).eq.'6 ') then
          weight=weight3d7(kk)
        endif
!
!     copying the shape functions, their derivatives and the
!     Jacobian determinant from field var
!
        do jj=1,nope
          do ii=1,4
            index=index+1
            shp(ii,jj)=var(index)
          enddo
        enddo
        index=index+1
        xsj=var(index)
!
        xsjmod=dtime*xsj*weight
!
!     calculating of the temperature gradient dtem
!     in the integration point
!
        enthalpy=0.d0
        do i1=1,3
          dtem(i1)=0.d0
        enddo
        do i1=1,nope
          do j1=1,3
            dtem(j1)=dtem(j1)+shp(j1,i1)*voldl(0,i1)
          enddo
        enddo
!
!     retrieving shpv and temp
!
        do i1=1,nope
          index=index+1
          shpv(i1)=var(index)
        enddo
        index=index+1
        temp=var(index)
!
!     calculating the total enthalpy
!     in the integration point
!
        do i1=1,nope
          enthalpy=enthalpy+shpv(i1)*(vconl(0,i1)+voldl(4,i1))
        enddo
!
!     material data (conductivity)
!
        call materialdata_cond(imat,ntmat_,temp,cocon,ncocon,cond)
!
!     retrieving stress x velocity
!
        do i1=1,3
          index=index+1
          tv(i1)=var(index)
        enddo
!
!     determining stress x velocity + conductivity x
!     temperature gradient
!
        do i1=1,3
          tv(i1)=tv(i1)+cond*dtem(i1)
        enddo
!
!     determination of the rhs of the energy equations
!
        do jj=1,nope
          ff(jj)=ff(jj)-xsjmod*(
     &         (shp(4,jj)+dt(konl(jj))*shpv(jj)/2.d0)*enthalpy+
     &         shp(1,jj)*tv(1)+shp(2,jj)*tv(2)+shp(3,jj)*tv(3))
        enddo
!
!     computation of contribution due to body forces
!
        if(nbody.ne.0) then
!
!     retrieving bfv (scalar product of the body force
!     with the velocity)
!
          index=index+1
          bfv=var(index)
          index=index+9
!
          do jj=1,nope
            ff(jj)=ff(jj)+xsjmod*(shp(4,jj)+
     &           dt(konl(jj))*shpv(jj)/2.d0)*bfv
          enddo
        else
          index=index+10
        endif
!
!     distributed heat flux
!
        if(nload.gt.0) then
          call nident2(nelemload,nelem,nload,id)
          areaj=xsj*weight
          do
            if((id.eq.0).or.(nelemload(1,id).ne.nelem)) exit
            if(sideload(id)(1:2).ne.'BF') then
              id=id-1
              cycle
            endif
            if(sideload(id)(3:4).eq.'NU') then
              do j=1,3
                pgauss(j)=0.d0
                do i1=1,nope
                  pgauss(j)=pgauss(j)+
     &                 shp(4,i1)*co(j,konl(i1))
                enddo
              enddo
              jltyp=1
              iscale=1
              call dflux(xload(1,id),temp,istep,iinc,tvar,
     &             nelem,kk,pgauss,jltyp,temp,press,sideload(id),
     &             areaj,vold,co,lakonl,konl,ipompc,nodempc,coefmpc,
     &             nmpc,ikmpc,ilmpc,iscale,mi)
              if((nmethod.eq.20).and.(iscale.ne.0))
     &             xload(1,id)=xloadold(1,id)+
     &             (xload(1,id)-xloadold(1,id))*reltimef
            endif
            do jj=1,nope
              ff(jj)=ff(jj)+xsjmod*(shp(4,jj)+
     &             dt(konl(jj))*shpv(jj)/2.d0)*xload(1,id)
            enddo
            exit
          enddo
        endif
!
      enddo
!
      if(nface.ne.0) then
        index=ipvarf(nelem)
!
!     free stream or solid surface boundaries
!
        nopes=0
        call nident(nelemface,nelem,nface,idf)
        do
          if((idf.eq.0).or.(nelemface(idf).ne.nelem)) exit
          ig=ichar(sideface(idf)(1:1))-48
!
!     check for distributed flux
!     an adiabatic face must be declared as a face with
!     distributed flux zero!
!
          flux=0
          call nident2(nelemload,nelem,nload,id)
          do
            if((id.eq.0).or.(nelemload(1,id).ne.nelem)) exit
            if((sideload(id)(1:1).ne.'F').and.
     &           (sideload(id)(1:1).ne.'R').and.
     &           (sideload(id)(1:1).ne.'S')) then
              id=id-1
              cycle
            endif
            igl=ichar(sideload(id)(2:2))-48
            if(igl.ne.ig) then
              id=id-1
              cycle
            endif
            flux=1
            exit
          enddo
!
          if(nopes.eq.0) then
            if(lakonl(4:4).eq.'4') then
              nopes=3
              mint2d=1
            elseif(lakonl(4:4).eq.'6') then
              mint2d=1
            elseif(lakonl(4:5).eq.'8R') then
              nopes=4
              mint2d=1
            elseif(lakonl(4:4).eq.'8') then
              nopes=4
              mint2d=4
            endif
          endif
!
          if(lakonl(4:4).eq.'6') then
            if(ig.le.2) then
              nopes=3
            else
              nopes=4
            endif
          endif
!
          do i=1,mint2d
!
!     local coordinates of the surface integration
!     point within the surface local coordinate system
!
            if((lakonl(4:5).eq.'8R').or.
     &           ((lakonl(4:4).eq.'6').and.(nopes.eq.4))) then
              weight=weight2d1(i)
            elseif(lakonl(4:4).eq.'8') then
              weight=weight2d2(i)
            elseif((lakonl(4:4).eq.'4').or.
     &             ((lakonl(4:4).eq.'6').and.(nopes.eq.3))) then
              weight=weight2d4(i)
            endif
!
!     local surface normal
!
            do i1=1,nopes
              index=index+1
              shp2(4,i1)=varf(index)
            enddo
            do i1=1,3
              index=index+1
              xsj2(i1)=varf(index)
            enddo
!
!     derivative of the volumetric shape functions
!     needed for the dissipative stress
!
            do i1=1,nope
              do j1=1,4
                index=index+1
                shp(j1,i1)=varf(index)
              enddo
            enddo
!
!     calculating of the temperature gradient dtem
!     in the integration point
!
            do i1=1,3
              dtem(i1)=0.d0
            enddo
            do i1=1,nope
              do j1=1,3
                dtem(j1)=dtem(j1)+shp(j1,i1)*voldl(0,i1)
              enddo
            enddo
!
!     retrieving the temperature
!
            index=index+1
            varf(index)=temp
!
!     material data (conductivity)
!
            call materialdata_cond(imat,ntmat_,temp,cocon,ncocon,
     &           cond)
!
!     determining  stress x velocity + conductivity x
!     temperature gradient
!
            do i1=1,3
              index=index+1
              tv(i1)=varf(index)
              if(flux.eq.0) then
                tv(i1)=tv(i1)+cond*dtem(i1)
              endif
            enddo
            index=index+5
!
            tvn=tv(1)*xsj2(1)+tv(2)*xsj2(2)+tv(3)*xsj2(3)
!
            if(flux.eq.1) then
              dxsj2=dsqrt(xsj2(1)*xsj2(1)+xsj2(2)*xsj2(2)+
     &             xsj2(3)*xsj2(3))
              areaj=dxsj2*weight
              sinktemp=xload(2,id)
!
!     for nonuniform load: determine the coordinates of the
!     point (transferred into the user subroutine)
!
              if((sideload(id)(3:4).eq.'NU').or.
     &             (sideload(id)(5:6).eq.'NU')) then
                if(nope.eq.8) then
                  do k=1,3
                    coords(k)=0.d0
                    do j=1,nopes
                      coords(k)=coords(k)+
     &                     co(k,konl(ifaceq(j,ig)))*shp2(4,j)
                    enddo
                  enddo
                elseif(nope.eq.4) then
                  do k=1,3
                    coords(k)=0.d0
                    do j=1,nopes
                      coords(k)=coords(k)+
     &                     co(k,konl(ifacet(j,ig)))*shp2(4,j)
                    enddo
                  enddo
                else
                  do k=1,3
                    coords(k)=0.d0
                    do j=1,nopes
                      coords(k)=coords(k)+
     &                     co(k,konl(ifacew(j,ig)))*shp2(4,j)
                    enddo
                  enddo
                endif
                jltyp=ichar(sideload(id)(2:2))-48
                jltyp=jltyp+10
                if(sideload(id)(1:1).eq.'S') then
                  iscale=1
                  call dflux(xload(1,id),temp,istep,iinc,tvar,
     &                 nelem,i,coords,jltyp,temp,press,
     &                 sideload(id),areaj,vold,co,lakonl,konl,
     &                 ipompc,nodempc,coefmpc,nmpc,ikmpc,ilmpc,
     &                 iscale,mi)
                  if((nmethod.eq.1).and.(iscale.ne.0))
     &                 xload(1,id)=xloadold(1,id)+
     &                 (xload(1,id)-xloadold(1,id))*reltimef
                elseif(sideload(id)(1:1).eq.'F') then
                  call film(xload(1,id),sinktemp,temp,istep,
     &                 iinc,tvar,nelem,i,coords,jltyp,field,
     &                 nfield,sideload(id),node,areaj,vold,mi)
                  if(nmethod.eq.1) xload(1,id)=xloadold(1,id)+
     &                 (xload(1,id)-xloadold(1,id))*reltimef
                elseif(sideload(id)(1:1).eq.'R') then
                  call radiate(xload(1,id),xload(2,id),temp,istep,
     &                 iinc,tvar,nelem,i,coords,jltyp,field,
     &                 nfield,sideload(id),node,areaj,vold,mi,
     &                 iemchange)
                  if(nmethod.eq.1) xload(1,id)=xloadold(1,id)+
     &                 (xload(1,id)-xloadold(1,id))*reltimef
                endif
              endif
!
              if(sideload(id)(1:1).eq.'S') then
!
!     flux INTO the face is positive (input deck convention)
!     this is different from the convention in the theory
!
                tvn=tvn+xload(1,id)*dxsj2
              elseif(sideload(id)(1:1).eq.'F') then
                tvn=tvn-xload(1,id)*(temp-sinktemp)*dxsj2
              elseif(sideload(id)(1:1).eq.'R') then
                tvn=tvn-physcon(2)*
     &               xload(1,id)*((temp-physcon(1))**4-
     &               (xload(2,id)-physcon(1))**4)*dxsj2
              endif
            endif
!
            xsjmod=tvn*weight*dtime
            do k=1,nopes
              if(nope.eq.8) then
                ipointer=ifaceq(k,ig)
              elseif(nope.eq.4) then
                ipointer=ifacet(k,ig)
              else
                ipointer=ifacew(k,ig)
              endif
              ff(ipointer)=ff(ipointer)+
     &             shp2(4,k)*xsjmod
            enddo
          enddo
          idf=idf-1
        enddo
      endif
!
      return
      end

ext: f

max_stars_repo_name: alleindrach/calculix-desktop

max_issues_repo_name: alleindrach/calculix-desktop

max_forks_repo_name: alleindrach/calculix-desktop


text:
      SUBROUTINE AB01OD( STAGES, JOBU, JOBV, N, M, A, LDA, B, LDB, U,
     $                   LDU, V, LDV, NCONT, INDCON, KSTAIR, TOL, IWORK,
     $                   DWORK, LDWORK, INFO )
C
C     SLICOT RELEASE 5.7.
C
C     Copyright (c) 2002-2020 NICONET e.V.
C
C     PURPOSE
C
C     To reduce the matrices A and B using (and optionally accumulating)
C     state-space and input-space transformations U and V respectively,
C     such that the pair of matrices
C
C        Ac = U' * A * U,    Bc = U' * B * V
C
C     are in upper "staircase" form. Specifically,
C
C             [ Acont     *    ]         [ Bcont ]
C        Ac = [                ],   Bc = [       ],
C             [   0    Auncont ]         [   0   ]
C
C        and
C
C                [ A11 A12  . . .  A1,p-1 A1p ]         [ B1 ]
C                [ A21 A22  . . .  A2,p-1 A2p ]         [ 0  ]
C                [  0  A32  . . .  A3,p-1 A3p ]         [ 0  ]
C        Acont = [  .   .   . . .    .     .  ],   Bc = [ .  ],
C                [  .   .     . .    .     .  ]         [ .  ]
C                [  .   .       .    .     .  ]         [ .  ]
C                [  0   0   . . .  Ap,p-1 App ]         [ 0  ]
C
C     where the blocks  B1, A21, ..., Ap,p-1  have full row ranks and
C     p is the controllability index of the pair.  The size of the
C     block Auncont is equal to the dimension of the uncontrollable
C     subspace of the pair (A, B).  The first stage of the reduction,
C     the "forward" stage, accomplishes the reduction to the orthogonal
C     canonical form (see SLICOT library routine AB01ND). The blocks
C     B1, A21, ..., Ap,p-1 are further reduced in a second, "backward"
C     stage to upper triangular form using RQ factorization. Each of
C     these stages is optional.
C
C     ARGUMENTS
C
C     Mode Parameters
C
C     STAGES  CHARACTER*1
C             Specifies the reduction stages to be performed as follows:
C             = 'F':  Perform the forward stage only;
C             = 'B':  Perform the backward stage only;
C             = 'A':  Perform both (all) stages.
C
C     JOBU    CHARACTER*1
C             Indicates whether the user wishes to accumulate in a
C             matrix U the state-space transformations as follows:
C             = 'N':  Do not form U;
C             = 'I':  U is internally initialized to the unit matrix (if
C                     STAGES <> 'B'), or updated (if STAGES = 'B'), and
C                     the orthogonal transformation matrix U is
C                     returned.
C
C     JOBV    CHARACTER*1
C             Indicates whether the user wishes to accumulate in a
C             matrix V the input-space transformations as follows:
C             = 'N':  Do not form V;
C             = 'I':  V is initialized to the unit matrix and the
C                     orthogonal transformation matrix V is returned.
C             JOBV is not referenced if STAGES = 'F'.
C
C     Input/Output Parameters
C
C     N       (input) INTEGER
C             The actual state dimension, i.e. the order of the
C             matrix A.  N >= 0.
C
C     M       (input) INTEGER
C             The actual input dimension.  M >= 0.
C
C     A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
C             On entry, the leading N-by-N part of this array must
C             contain the state transition matrix A to be transformed.
C             If STAGES = 'B', A should be in the orthogonal canonical
C             form, as returned by SLICOT library routine AB01ND.
C             On exit, the leading N-by-N part of this array contains
C             the transformed state transition matrix U' * A * U.
C             The leading NCONT-by-NCONT part contains the upper block
C             Hessenberg state matrix Acont in Ac, given by U' * A * U,
C             of a controllable realization for the original system.
C             The elements below the first block-subdiagonal are set to
C             zero.  If STAGES <> 'F', the subdiagonal blocks of A are
C             triangularized by RQ factorization, and the annihilated
C             elements are explicitly zeroed.
C
C     LDA     INTEGER
C             The leading dimension of array A.  LDA >= MAX(1,N).
C
C     B       (input/output) DOUBLE PRECISION array, dimension (LDB,M)
C             On entry, the leading N-by-M part of this array must
C             contain the input matrix B to be transformed.
C             If STAGES = 'B', B should be in the orthogonal canonical
C             form, as returned by SLICOT library routine AB01ND.
C             On exit with STAGES = 'F', the leading N-by-M part of
C             this array contains the transformed input matrix U' * B,
C             with all elements but the first block set to zero.
C             On exit with STAGES <> 'F', the leading N-by-M part of
C             this array contains the transformed input matrix
C             U' * B * V, with all elements but the first block set to
C             zero and the first block in upper triangular form.
C
C     LDB     INTEGER
C             The leading dimension of array B.  LDB >= MAX(1,N).
C
C     U       (input/output) DOUBLE PRECISION array, dimension (LDU,N)
C             If STAGES <> 'B' or JOBU = 'N', then U need not be set
C             on entry.
C             If STAGES = 'B' and JOBU = 'I', then, on entry, the
C             leading N-by-N part of this array must contain the
C             transformation matrix U that reduced the pair to the
C             orthogonal canonical form.
C             On exit, if JOBU = 'I', the leading N-by-N part of this
C             array contains the transformation matrix U that performed
C             the specified reduction.
C             If JOBU = 'N', the array U is not referenced and can be
C             supplied as a dummy array (i.e. set parameter LDU = 1 and
C             declare this array to be U(1,1) in the calling program).
C
C     LDU     INTEGER
C             The leading dimension of array U.
C             If JOBU = 'I', LDU >= MAX(1,N);  if JOBU = 'N', LDU >= 1.
C
C     V       (output) DOUBLE PRECISION array, dimension (LDV,M)
C             If JOBV = 'I', then the leading M-by-M part of this array
C             contains the transformation matrix V.
C             If STAGES = 'F', or JOBV = 'N', the array V is not
C             referenced and can be supplied as a dummy array (i.e. set
C             parameter  LDV = 1 and declare this array to be V(1,1) in
C             the calling program).
C
C     LDV     INTEGER
C             The leading dimension of array V.
C             If STAGES <> 'F' and JOBV = 'I', LDV >= MAX(1,M);
C             if STAGES = 'F' or JOBV = 'N', LDV >= 1.
C
C     NCONT   (input/output) INTEGER
C             The order of the controllable state-space representation.
C             NCONT is input only if STAGES = 'B'.
C
C     INDCON  (input/output) INTEGER
C             The number of stairs in the staircase form (also, the
C             controllability index of the controllable part of the
C             system representation).
C             INDCON is input only if STAGES = 'B'.
C
C     KSTAIR  (input/output) INTEGER array, dimension (N)
C             The leading INDCON elements of this array contain the
C             dimensions of the stairs, or, also, the orders of the
C             diagonal blocks of Acont.
C             KSTAIR is input if STAGES = 'B', and output otherwise.
C
C     Tolerances
C
C     TOL     DOUBLE PRECISION
C             The tolerance to be used in rank determination when
C             transforming (A, B). If the user sets TOL > 0, then
C             the given value of TOL is used as a lower bound for the
C             reciprocal condition number (see the description of the
C             argument RCOND in the SLICOT routine MB03OD);  a
C             (sub)matrix whose estimated condition number is less than
C             1/TOL is considered to be of full rank.  If the user sets
C             TOL <= 0, then an implicitly computed, default tolerance,
C             defined by  TOLDEF = N*N*EPS,  is used instead, where EPS
C             is the machine precision (see LAPACK Library routine
C             DLAMCH).
C             TOL is not referenced if STAGES = 'B'.
C
C     Workspace
C
C     IWORK   INTEGER array, dimension (M)
C             IWORK is not referenced if STAGES = 'B'.
C
C     DWORK   DOUBLE PRECISION array, dimension (LDWORK)
C             On exit, if INFO = 0, DWORK(1) returns the optimal value
C             of LDWORK.
C
C     LDWORK  INTEGER
C             The length of the array DWORK.
C             If STAGES <> 'B', LDWORK >= MAX(1, N + MAX(N,3*M));
C             If STAGES =  'B', LDWORK >= MAX(1, M + MAX(N,M)).
C             For optimum performance LDWORK should be larger.
C
C     Error Indicator
C
C     INFO    INTEGER
C             = 0:  successful exit;
C             < 0:  if INFO = -i, the i-th argument had an illegal
C                   value.
C
C     METHOD
C
C     Staircase reduction of the pencil [B|sI - A] is used. Orthogonal
C     transformations U and V are constructed such that
C
C
C                        |B |sI-A      *  . . .  *      *       |
C                        | 1|    11       .      .      .       |
C                        |  |  A    sI-A    .    .      .       |
C                        |  |   21      22    .  .      .       |
C                        |  |        .     .     *      *       |
C     [U'BV|sI - U'AU] = |0 |     0    .     .                  |
C                        |  |            A     sI-A     *       |
C                        |  |             p,p-1    pp           |
C                        |  |                                   |
C                        |0 |         0          0   sI-A       |
C                        |  |                            p+1,p+1|
C
C
C     where the i-th diagonal block of U'AU has dimension KSTAIR(i),
C     for i = 1,...,p. The value of p is returned in INDCON. The last
C     block contains the uncontrollable modes of the (A,B)-pair which
C     are also the generalized eigenvalues of the above pencil.
C
C     The complete reduction is performed in two stages. The first,
C     forward stage accomplishes the reduction to the orthogonal
C     canonical form. The second, backward stage consists in further
C     reduction to triangular form by applying left and right orthogonal
C     transformations.
C
C     REFERENCES
C
C     [1] Van Dooren, P.
C         The generalized eigenvalue problem in linear system theory.
C         IEEE Trans. Auto. Contr., AC-26, pp. 111-129, 1981.
C
C     [2] Miminis, G. and Paige, C.
C         An algorithm for pole assignment of time-invariant multi-input
C         linear systems.
C         Proc. 21st IEEE CDC, Orlando, Florida, 1, pp. 62-67, 1982.
C
C     NUMERICAL ASPECTS
C
C     The algorithm requires O((N + M) x N**2) operations and is
C     backward stable (see [1]).
C
C     FURTHER COMMENTS
C
C     If the system matrices A and B are badly scaled, it would be
C     useful to scale them with SLICOT routine TB01ID, before calling
C     the routine.
C
C     CONTRIBUTOR
C
C     Release 3.0: V. Sima, Katholieke Univ. Leuven, Belgium, Nov. 1996.
C     Supersedes Release 2.0 routine AB01CD by M. Vanbegin, and
C     P. Van Dooren, Philips Research Laboratory, Brussels, Belgium.
C
C     REVISIONS
C
C     January 14, 1997, February 12, 1998, September 22, 2003.
C
C     KEYWORDS
C
C     Controllability, generalized eigenvalue problem, orthogonal
C     transformation, staircase form.
C
C     ******************************************************************
C
C     .. Parameters ..
      DOUBLE PRECISION  ZERO, ONE
      PARAMETER         ( ZERO = 0.0D0, ONE = 1.0D0 )
C     .. Scalar Arguments ..
      CHARACTER         JOBU, JOBV, STAGES
      INTEGER           INDCON, INFO, LDA, LDB, LDU, LDV, LDWORK, M, N,
     $                  NCONT
      DOUBLE PRECISION  TOL
C     .. Array Arguments ..
      INTEGER           IWORK(*), KSTAIR(*)
      DOUBLE PRECISION  A(LDA,*), B(LDB,*), DWORK(*), U(LDU,*), V(LDV,*)
C     .. Local Scalars ..
      LOGICAL           LJOBUI, LJOBVI, LSTAGB, LSTGAB
      INTEGER           I, I0, IBSTEP, ITAU, J0, JINI, JWORK, MCRT, MM,
     $                  NCRT, WRKOPT
C     .. External Functions ..
      LOGICAL           LSAME
      EXTERNAL          LSAME
C     .. External Subroutines ..
      EXTERNAL          AB01ND, DGERQF, DLACPY, DLASET, DORGRQ, DORMRQ,
     $                  DSWAP, XERBLA
C     .. Intrinsic Functions ..
      INTRINSIC         INT, MAX, MIN
C     .. Executable Statements ..
C
      INFO = 0
      LJOBUI = LSAME( JOBU, 'I' )
C
      LSTAGB = LSAME( STAGES, 'B' )
      LSTGAB = LSAME( STAGES, 'A' ).OR.LSTAGB
C
      IF ( LSTGAB ) THEN
         LJOBVI = LSAME( JOBV, 'I' )
      END IF
C
C     Test the input scalar arguments.
C
      IF( .NOT.LSTGAB .AND. .NOT.LSAME( STAGES, 'F' ) ) THEN
         INFO = -1
      ELSE IF( .NOT.LJOBUI .AND. .NOT.LSAME( JOBU, 'N' ) ) THEN
         INFO = -2
      ELSE IF( N.LT.0 ) THEN
         INFO = -4
      ELSE IF( M.LT.0 ) THEN
         INFO = -5
      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
         INFO = -7
      ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
         INFO = -9
      ELSE IF( LDU.LT.1 .OR. ( LJOBUI .AND. LDU.LT.N ) ) THEN
         INFO = -11
      ELSE IF( .NOT.LSTAGB .AND. LDWORK.LT.MAX( 1, N + MAX( N, 3*M ) )
     $         .OR. LSTAGB .AND. LDWORK.LT.MAX( 1, M + MAX( N, M ) ) )
     $         THEN
         INFO = -20
      ELSE IF( LSTAGB .AND. NCONT.GT.N ) THEN
         INFO = -14
      ELSE IF( LSTAGB .AND. INDCON.GT.N ) THEN
         INFO = -15
      ELSE IF( LSTGAB ) THEN
         IF( .NOT.LJOBVI .AND. .NOT.LSAME( JOBV, 'N' ) ) THEN
            INFO = -3
         ELSE IF( LDV.LT.1 .OR. ( LJOBVI .AND. LDV.LT.M ) ) THEN
            INFO = -13
         END IF
      END IF
C
      IF ( INFO.NE.0 ) THEN
C
C        Error return.
C
         CALL XERBLA( 'AB01OD', -INFO )
         RETURN
      END IF
C
C     Quick return if possible.
C
      IF ( MIN( N, M ).EQ.0 ) THEN
         NCONT  = 0
         INDCON = 0
         IF( N.GT.0 .AND. LJOBUI )
     $      CALL DLASET( 'F', N, N, ZERO, ONE, U, LDU )
         IF( LSTGAB ) THEN
            IF( M.GT.0 .AND. LJOBVI )
     $         CALL DLASET( 'F', M, M, ZERO, ONE, V, LDV )
         END IF
         DWORK(1) = ONE
         RETURN
      END IF
C
C     (Note: Comments in the code beginning "Workspace:" describe the
C     minimal amount of real workspace needed at that point in the
C     code, as well as the preferred amount for good performance.
C     NB refers to the optimal block size for the immediately
C     following subroutine, as returned by ILAENV.)
C
      ITAU = 1
      WRKOPT = 1
C
      IF ( .NOT.LSTAGB ) THEN
C
C        Perform the forward stage computations of the staircase
C        algorithm on B and A: reduce the (A, B) pair to orthogonal
C        canonical form.
C
C        Workspace: N + MAX(N,3*M).
C
         JWORK = N + 1
         CALL AB01ND( JOBU, N, M, A, LDA, B, LDB, NCONT, INDCON,
     $                KSTAIR, U, LDU, DWORK(ITAU), TOL, IWORK,
     $                DWORK(JWORK), LDWORK-JWORK+1, INFO )
C
         WRKOPT = INT( DWORK(JWORK) ) + JWORK - 1
      END IF
C
C     Exit if no further reduction to triangularize B1 and subdiagonal
C     blocks of A is required, or if the order of the controllable part
C     is 0.
C
      IF ( .NOT.LSTGAB ) THEN
         DWORK(1) = WRKOPT
         RETURN
      ELSE IF ( NCONT.EQ.0 .OR. INDCON.EQ.0 ) THEN
         IF( LJOBVI )
     $      CALL DLASET( 'F', M, M, ZERO, ONE, V, LDV )
         DWORK(1) = WRKOPT
         RETURN
      END IF
C
C     Now perform the backward steps except the last one.
C
      MCRT = KSTAIR(INDCON)
      I0 = NCONT - MCRT + 1
      JWORK = M + 1
C
      DO 10 IBSTEP = INDCON, 2, -1
         NCRT = KSTAIR(IBSTEP-1)
         J0 = I0 - NCRT
         MM = MIN( NCRT, MCRT )
C
C        Compute the RQ factorization of the current subdiagonal block
C        of A, Ai,i-1 = R*Q (where i is IBSTEP), of dimension
C        MCRT-by-NCRT, starting in position (I0,J0).
C        The matrix Q' should postmultiply U, if required.
C        Workspace: need   M + MCRT;
C                   prefer M + MCRT*NB.
C
         CALL DGERQF( MCRT, NCRT, A(I0,J0), LDA, DWORK(ITAU),
     $                DWORK(JWORK), LDWORK-JWORK+1, INFO )
         WRKOPT = MAX( WRKOPT, INT( DWORK(JWORK) )+JWORK-1 )
C
C        Set JINI to the first column number in A where the current
C        transformation Q is to be applied, taking the block Hessenberg
C        form into account.
C
         IF ( IBSTEP.GT.2 ) THEN
            JINI = J0 - KSTAIR(IBSTEP-2)
         ELSE
            JINI = 1
C
C           Premultiply the first block row (B1) of B by Q.
C           Workspace: need   2*M;
C                      prefer M + M*NB.
C
            CALL DORMRQ( 'Left', 'No transpose', NCRT, M, MM, A(I0,J0),
     $                   LDA, DWORK(ITAU), B, LDB, DWORK(JWORK),
     $                   LDWORK-JWORK+1, INFO )
            WRKOPT = MAX( WRKOPT, INT( DWORK(JWORK) )+JWORK-1 )
         END IF
C
C        Premultiply the appropriate block row of A by Q.
C        Workspace: need   M + N;
C                   prefer M + N*NB.
C
         CALL DORMRQ( 'Left', 'No transpose', NCRT, N-JINI+1, MM,
     $                A(I0,J0), LDA, DWORK(ITAU), A(J0,JINI), LDA,
     $                DWORK(JWORK), LDWORK-JWORK+1, INFO )
         WRKOPT = MAX( WRKOPT, INT( DWORK(JWORK) )+JWORK-1 )
C
C        Postmultiply the appropriate block column of A by Q'.
C        Workspace: need   M +  I0-1;
C                   prefer M + (I0-1)*NB.
C
         CALL DORMRQ( 'Right', 'Transpose', I0-1, NCRT, MM, A(I0,J0),
     $                LDA, DWORK(ITAU), A(1,J0), LDA, DWORK(JWORK),
     $                LDWORK-JWORK+1, INFO )
         WRKOPT = MAX( WRKOPT, INT( DWORK(JWORK) )+JWORK-1 )
C
         IF ( LJOBUI ) THEN
C
C           Update U, postmultiplying it by Q'.
C           Workspace: need   M + N;
C                      prefer M + N*NB.
C
            CALL DORMRQ( 'Right', 'Transpose', N, NCRT, MM, A(I0,J0),
     $                   LDA, DWORK(ITAU), U(1,J0), LDU, DWORK(JWORK),
     $                   LDWORK-JWORK+1, INFO )
            WRKOPT = MAX( WRKOPT, INT( DWORK(JWORK) )+JWORK-1 )
         END IF
C
C        Zero the subdiagonal elements of the current subdiagonal block
C        of A.
C
         CALL DLASET( 'F', MCRT, NCRT-MCRT, ZERO, ZERO, A(I0,J0), LDA )
         IF ( I0.LT.N )
     $      CALL DLASET( 'L', MCRT-1, MCRT-1, ZERO, ZERO,
     $                   A(I0+1,I0-MCRT), LDA )
C
         MCRT = NCRT
         I0 = J0
C
   10 CONTINUE
C
C     Now perform the last backward step on B, V = Qb'.
C
C     Compute the RQ factorization of the first block of B, B1 = R*Qb.
C     Workspace: need   M + MCRT;
C                prefer M + MCRT*NB.
C
      CALL DGERQF( MCRT, M, B, LDB, DWORK(ITAU), DWORK(JWORK),
     $             LDWORK-JWORK+1, INFO )
      WRKOPT = MAX( WRKOPT, INT( DWORK(JWORK) )+JWORK-1 )
C
      IF ( LJOBVI ) THEN
C
C        Accumulate the input-space transformations V.
C        Workspace: need 2*M;  prefer M + M*NB.
C
         CALL DLACPY( 'F', MCRT, M-MCRT, B, LDB, V(M-MCRT+1,1), LDV )
         IF ( MCRT.GT.1 )
     $      CALL DLACPY( 'L', MCRT-1, MCRT-1, B(2,M-MCRT+1), LDB,
     $                   V(M-MCRT+2,M-MCRT+1), LDV )
         CALL DORGRQ( M, M, MCRT, V, LDV, DWORK(ITAU), DWORK(JWORK),
     $                LDWORK-JWORK+1, INFO )
C
         DO 20 I = 2, M
            CALL DSWAP( I-1, V(I,1), LDV, V(1,I), 1 )
   20    CONTINUE
C
         WRKOPT = MAX( WRKOPT, INT( DWORK(JWORK) )+JWORK-1 )
      END IF
C
C     Zero the subdiagonal elements of the submatrix B1.
C
      CALL DLASET( 'F', MCRT, M-MCRT, ZERO, ZERO, B, LDB )
      IF ( MCRT.GT.1 )
     $   CALL DLASET( 'L', MCRT-1, MCRT-1, ZERO, ZERO, B(2,M-MCRT+1),
     $                LDB )
C
C     Set optimal workspace dimension.
C
      DWORK(1) = WRKOPT
      RETURN
C *** Last line of AB01OD ***
      END

max_stars_repo_name: bnavigator/SLICOT-Reference

max_issues_repo_name: RJHKnight/slicotr

max_forks_repo_name: RJHKnight/slicotr


text:
      SUBROUTINE BOUND
C
C     TO OBTAIN THE CORRECTION DUE TO THE TURBULENT BOUNDARY LAYER
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMMON /GG/ GAM,GM,G1,G2,G3,G4,G5,G6,G7,G8,G9,GA,RGA,QT
      COMMON /CORR/ DLA(200),RCO(200),DAX(200),DRX(200),SL(200),DR2
      COMMON /COORD/ S(200),FS(200),WALTAN(200),SD(200),WMN(200),TTR(200
     1),DMDX(200),SPR(200),BTA(200),SREF(200),XBIN,XCIN,GMA,GMB,GMC,GMD
      COMMON /PROP/ AR,ZO,RO,VISC,VISM,SFOA,SBL,CONV
      COMMON /PARAM/ ETAD,RC,AMACH,BMACH,CMACH,EMACH,GMACH,FRC,SF,WWO,WW
     1OP,QM,WE,CBET,XE,ETA,EPSI,BPSI,XO,YO,RRC,SDO,XB,XC,AH,PP,SE,TYE,XA
      COMMON /HTTR/ HAIR,TAW,TWQ,TW,TWAT,QFUN,QFUNW,IPQ,IJ,IV,IW
      COMMON /CONTR/ ITLE(3),IE,LR,IT,JB,JQ,JX,KAT,KBL,KING,KO,LV,NOCON,
     1IN,MC,MCP
      DIMENSION Z(16), D(16), SCV(200), SK(200), CDS(200), RW(200)
      DATA ZRO/0.0D+0/,ONE/1.D+0/,TWO/2.D+0/,SIX/6.D+0/,HALF/5.D-1/
      DATA THR/3.D+0/,FOUR/4.D+0/,TEN/1.D+1/,TLV/1.2D+1/
      DATA CF1/3.865D-2/,CF2/4.561D+0/,CF3/5.46D-1/,FSI/3.17897971D+0/
      DATA LY/4H   Y/,LS/4H   S/,DD/8HD2Y/DX2 /,DK/8H CURV.  /
      DATA Z(1)/.052995325D-1/,Z(4)/.1222977958D+0/,Z(7)/.3591982246D+0/
      DATA Z(2)/.277124885D-1/,Z(5)/.1910618778D+0/,Z(8)/.4524937451D+0/
      DATA Z(3)/.671843988D-1/,Z(6)/.2709916112D+0/
      DATA D(1)/.135762297D-1/,D(2)/.31126762D-1/,D(3)/.475792558D-1/
      DATA D(4)/.623144856D-1/,D(5)/.747979944D-1/,D(6)/.845782597D-1/
      DATA D(7)/.913017075D-1/,D(8)/.947253052D-1/
      DO 1 J=9,16
      D(J)=D(17-J)
1     Z(J)=ONE-Z(17-J)
      DO 2 J=1,KAT
2     SREF(J)=S(J)
      SBIN=XBIN
      SCIN=XCIN
      TRPI=CONV/90.D+0
      FCC=2.05D+0+DLOG(.41D+0)
      CHAIR=GAM*G1*AR/RO/RO/777.64885D+0
      IF (IT .EQ. 0) XBL=SBL
C
3     READ (1,66,END=65) PPQ,T0,TWT,TWAT,QFUN,ALPH,IHT,IR,LD,LV
C
      PPS=PPQ
      RHO=144.D+0*PPS/ZO/AR/T0
      ID=IABS(LD)
      KOR=KO
      IF (IABS(IN) .EQ. 10) KOR=1
      IF (MCP .LT. 0) KOR=KING
      ROY=ONE
      IF (IE .EQ. 0) HW=AH
      IF ((ID .EQ. 0) .OR. (IE .EQ. 1)) HW=ZRO
      IF (HW .EQ. ZRO) YOH=ZRO
      IF (HW .EQ. ZRO) YOHA=ZRO
      ALF=DABS(ALPH)
      ARC=FRC
      IF (IHT .LT. 0) ARC=FRC**(IE+1)
      IPQ=0
      IW=1
      IF (LV .NE. 0) IW=IABS(LV)
      DO 4 J=1,KAT
      S(J)=SREF(J)
      SL(J)=S(J)
      RW(J)=FS(J)
      RCO(J)=FS(J)
      SCW=DSQRT(ONE+WALTAN(J)**2)
      SK(J)=SD(J)/SCW**3
      IF (KAT .EQ. KING) GO TO 4
      IF (S(J) .LT. SBL) KBL=J+2
4     DRX(J)=WALTAN(J)
      IF (KBL .GT. KAT) KBL=KAT+4
5     DO 58 IV=1,IW
      IF ((IV .GT. 1) .AND. (IV .LT. IW)) GO TO 15
      IF (LD .GE. 0) WRITE (2,80) ITLE,PPS,T0
      IF (ALPH .GT. ZRO) GO TO 6
      ALPHA=ZRO
      IF ((LD .GE. 0) .OR. (PPQ .EQ. ZRO)) WRITE (2,71)
      GO TO 7
6     ALPHA=ALPH
      IF ((LD .GE. 0) .OR. (PPQ .EQ. ZRO)) WRITE (2,70)
7     IF (IR .EQ. 2) GO TO 13
      IF (ALF .EQ. ONE) GO TO 8
      IF ((LD .GE. 0) .OR. (PPQ .EQ. ZRO)) WRITE (2,75)
      GO TO 9
8     IF ((LD .GE. 0) .OR. (PPQ .EQ. ZRO)) WRITE (2,72)
9     IF (IR) 10,11,12
10    IF ((LD .GE. 0) .OR. (PPQ .EQ. ZRO)) WRITE (2,74)
      GO TO 14
11    IF ((LD .GE. 0) .OR. (PPQ .EQ. ZRO)) WRITE (2,73)
      GO TO 14
12    IF ((LD .GE. 0) .OR. (PPQ .EQ. ZRO)) WRITE (2,76)
      GO TO 14
13    IF ((LD .GE. 0) .OR. (PPQ .EQ. ZRO)) WRITE (2,77)
14    IF (PPQ .EQ. ZRO) GO TO 60
15    CAPI=.55D+0
      IPP=0
      IJ=1
      DO 56 J=1,KAT
      BET=TTR(J)-ONE
      STR=ONE/TTR(J)
      TE=T0*STR
      RAJ=WMN(J)*(G7*STR)**GA
      IF (IHT .GE. 0) RAJ=RAJ**QT
      SCW=DSQRT(ONE+DRX(J)**2)
      EMU=VISC*TE*DSQRT(TE)/(TE+VISM)
      IF (TE .LT. VISM) EMU=HALF*VISC*TE/DSQRT(VISM)
      IF (VISM .LE. ONE) EMU=VISC*TE**VISM
      TAW=TE*(ONE+RO*BET)
      RHOE=RHO*STR**G1
      VE=WMN(J)*DSQRT(GAM*AR*TE)
      REO=RHOE*VE/EMU/TLV
      IF (HW .GT. ZRO) YOH=FS(J)/HW
      IF (IE .EQ. 0 .AND. HW .GT. ZRO) ROY=(HW/FS(J)+ONE)*TRPI
      K=J
      IF (J .EQ. 1) GO TO 19
      IF (J .GT. KOR) K=J-KOR+1
      IF (K-3) 16,17,18
16    DS=S(J)-S(J-1)
      SMD=HALF*DS
      GO TO 19
17    DT=S(J)-S(J-1)
      DST=DS+DT
      SMA=DST*(TWO-DT/DS)/SIX
      SMC=DST*(TWO-DS/DT)/SIX
      SMB=DST-SMA-SMC
      HB=H
      IF (IV .GT. 1) GO TO 19
      BMA=TWO/DS/DST
      BMB=-TWO/DS/DT
      BMC=TWO/DT/DST
      GO TO 19
18    DU=S(J)-S(J-1)
      DT=S(J-1)-S(J-2)
      DS=S(J-2)-S(J-3)
      DST=DS+DT
      DSTU=DST+DU
      DTU=DT+DU
      DUT=DU-DT
      DTS=DS-DT
      DTUS=DT+TWO*(DU-DS)
      DTSU=DT+TWO*(DS-DU)
      DSTTU=TWO*(DST+DTU)
      HA=HB
      HB=H
      QMA=HALF*DS*(ONE-DS*(THR+(DTU+DU)/DST)/DSTU/SIX)
      QMB=HALF*DS*(ONE+DS*(TWO+(DST+DT)/DTU)/DT/SIX)
      QMC=-DS**3*(ONE+(DTU+DU)/DST)/DT/DU/TLV
      QMD=DS**3*(DST+DT)/DU/DTU/DSTU/TLV
      SMA=HALF*DS+(DUT*DTU**3/DS-DS*DS*(DS+DSTTU))/DST/DSTU/TLV
      SMB=HALF*DST+(DS*DS*(DSTTU-DS)/DT+DT*DT*DTUS/DS-DU**3*(DSTU+DST)/D
     1S/DT)/DTU/TLV
      SMC=HALF*DTU+(DT*DT*DTSU/DU+DU*DU*(DSTTU-DU)/DT-DS**3*(DSTU+DTU)/D
     1T/DU)/DST/TLV
      SMD=HALF*DU+(DTS*DST**3/DU-DU*DU*(DU+DSTTU))/DTU/DSTU/TLV
19    IF (TWT .NE. ZRO) GO TO 20
      TW=TAW
      GO TO 21
20    TWD=(ARC*RAJ-ONE)*(TWT-TWAT)/(ARC-ONE)
      IF (TWD .LT. ZRO) TWD=ZRO
      TW=TWD+TWAT
21    WMU=VISC*TW*DSQRT(TW)/(TW+VISM)
      IF (VISM .LE. ONE) WMU=VISC*TW**VISM
      DL=TW/TE
      DM=ALPHA*(TAW-TW)/TE
      DN=ONE-DL-DM
      DA=ALF*(TAW-TW)
      DB=DA+TW-TE
      IF (DB) 22,23,24
22    DG=DSQRT(-DB*TE)
      DH=DSQRT(-DB*TW)
      DI=(TWO*(DG+TE-TW)-DA)/(TWO*DH+DA)
      DJ=DLOG(DI)
      TP=-DB/DJ/DJ
      GO TO 25
23    TP=(DSQRT(TE)+DSQRT(TW))**2/FOUR
      GO TO 25
24    DC=DSQRT(DA*DA+FOUR*TW*DB)
      DF=DASIN((DB+TW-TE)/DC)
      DE=DASIN(DA/DC)
      TP=DB/(DF+DE)/(DF+DE)
25    IF (IR) 26,27,28
26    FRD=TW*EMU/WMU/TP
      GO TO 29
27    FRD=EMU/WMU
      GO TO 29
28    FRD=TE*EMU/WMU/DSQRT(TP*TW)
29    IF (IPP .GT. 0) GO TO 31
      RTHI=1.D-2*REO*FS(1)
      RTII=RTHI
      RDLI=TEN*RTHI
      IF (IR .EQ. 1) GO TO 32
30    RTHG=DLOG10(RTHI)
      CFI=CF1/(RTHG+CF2)/(RTHG-CF3)
31    IF (IR .NE. 2) GO TO 33
      SCFI=DSQRT(CFI)
      TC=TW+17.2D+0*SCFI*DA-305.D+0*CFI*DB
      CMU=VISC*TC*DSQRT(TC)/(TC+VISM)
      IF (VISM .LE. ONE) CMU=VISC*TC**VISM
      TP=TW*CMU/WMU
      FRD=EMU/CMU
      GO TO 33
32    RDLG=DLOG10(RDLI)
      CFI=0.0444D+0/(RDLG+4.6221D+0)/(RDLG-1.4402D+0)
33    CF=CFI*TE/TP
      CFS=CF*SCW
      RTIG=DLOG10(RTII)
      XCF=.41D+0*DSQRT((RTIG+CF2)*(RTIG-CF3)/CF1)
34    C3=TWO+CAPI*(FSI+1.5D+0*CAPI)
      C2=ONE+CAPI
      C1=C2-C3/XCF
      FXCF=XCF+DLOG(C1/RTII)-FCC-TWO*CAPI
      FPCP=(XCF-FSI-THR*CAPI)/XCF/C1-TWO
      CAPI=CAPI-FXCF/FPCP
      IF (DABS(FXCF) .GT. 1.D-8) GO TO 34
      DOTI=XCF/C1
      XN=HALF*(DOTI+DSQRT(DOTI*(DOTI-SIX)+ONE)-THR)
      HI=ONE+TWO/XN
      SUMA=ZRO
      SUMB=ZRO
      SUMC=ZRO
      SUMD=ZRO
      DO 35 L=1,16
      UN=Z(L)**XN
      TR=DL+Z(L)*(DM+Z(L)*DN)
      ADD=D(L)*XN*UN/TR
      BDD=ADD*Z(L)
      CDD=ADD*UN
      DDD=BDD*UN
      SUMA=SUMA+ADD
      SUMB=SUMB+BDD
      SUMC=SUMC+CDD
35    SUMD=SUMD+DDD
      DOT=ONE/(SUMA-SUMB)
      DSOD=ONE-SUMA
      DSM=HALF-SUMC
      THM=SUMC-SUMD
      HU=DSOD*DOT
      IF (IPP .GT. 0) GO TO 36
      H=HU
      DOTR=DOT
36    FMY=(H+TWO-G9*BET)*DMDX(J)*STR/WMN(J)+ID*DRX(J)/(RW(J)+HW)
      IF (J.EQ.1) TH=CFS/FMY
      IF (K.EQ.2) TH=(THA+SMD*(DTHA+CFS))/(ONE+SMD*FMY)
      IF (K.EQ.3) TH=(THA+SMA*DTHA+SMB*DTHB+SMC*CFS)/(ONE+SMC*FMY)
      IF (K.GT.3) TH=(THA+SMA*DTHA+SMB*DTHB+SMC*DTHC+SMD*CFS)/(ONE+SMD*F
     1MY)
      DELST=H*TH
      ASEC=DELST+DSQRT(ID*DELST**2+(FS(J)*SCW*ROY)**2)
      DOR=ID*DOTR*TH/ASEC
      DSROD=DSOD-DOR*DSM
      IPP=1
      DOTR=ONE/(ONE/DOT-THM*DOR)
      HR=DSROD*DOTR
      IF (DABS(H-HR) .LT. 5.D-7) GO TO 37
      H=HR
      GO TO 36
37    DELTA=DOTR*TH
      THU=DELTA/DOT
      DSU=DELTA*DSOD
      RDEL=REO*DELTA
      RTII=RDEL/DOTI
      RDLX=FRD*RDEL
      RTHX=RDLX/DOT
      IF (RTHX .LT. 100.D+0) GO TO 38
      IF (IR .EQ. 1) GO TO 39
      IF (DABS(ONE-RTHX/RTHI) .LT. 1.D-6) GO TO 41
      RTHI=RTHX
      GO TO 30
38    WRITE (2,88) RTHX,REO,FRD,TH,DELTA,DOT
      RETURN
39    IF (DABS(ONE-RDLX/RDLI) .LT. 1.D-6) GO TO 40
      RDLI=RDLX
      GO TO 32
40    RTHG=HALF*(DSQRT((CF2+CF3)**2+FOUR*CF1/CF1)-CF2+CF3)
      RTHX=TEN**RTHG
41    IF (J .GT. 1) GO TO 42
      DTH=ZRO
      HAIR=RHOE*VE*CF*CHAIR
      TAIR=HAIR
      IF (TWAT.EQ.TWT .OR. QFUN.EQ.ZRO) GO TO 46
      TWQ=(HAIR*TAW+QFUN*(TWAT-15.D+0))/(HAIR+QFUN)
      CALL HEAT
      IF (IPQ .GT. 100) GO TO 65
      IF (DABS(TW-TWQ).LT.1.D-2.AND.DABS(QFUN-QFUNW).LT.1.D-5) GO TO 46
      TWT=TWAT+(TWQ-TWAT)*(ARC-ONE)/(ARC*RAJ-ONE)
      QFUN=QFUNW
      GO TO 20
42    DTH=CFS-TH*FMY
      IF (DTH .LT. ZRO) DTH=ZRO
      IF (J .EQ. KOR) GO TO 46
      IF (K-3) 43,45,44
43    DTHB=DTH
      GO TO 47
44    THA=THA+QMA*DTHA+QMB*DTHB+QMC*DTHC+QMD*DTH
      DTHA=DTHB
      DTHB=DTHC
      IF (K .GT. 5) GO TO 45
      SCU=DSQRT(ONE+DRX(J-2)**2)
      DELA=HA*THA
      IF ((IE .EQ. 1).OR.(ID .EQ. 0)) YSEC=FS(J-2)*SCU
      IF ((IE .EQ. 0).AND.(HW .GT. ZRO)) YSEC=SCU*(FS(J-2)+HW)*TRPI
      IF (HW .GT. ZRO) YOHA=FS(J-2)/HW
      ASCA=DELA+DSQRT(ID*DELA**2+YSEC**2)
      RW(J-2)=ASCA/SCU
      DLA(J-2)=SCU*(ASCA-YSEC)*(ONE+YOHA)
      RCO(J-2)=FS(J-2)+DLA(J-2)
45    DTHC=DTH
      GO TO 47
46    THA=TH
      DTHA=DTH
      IF ((IV .GT. 1) .AND. (IV .LT. IW)) GO TO 47
      IF (J .EQ. 1 .AND. LD .GE. 0) WRITE (2,82)
47    CDS(J)=ASEC-SCW*FS(J)*ROY
      DLA(J)=SCW*CDS(J)*(ONE+YOH)
      RCO(J)=FS(J)+DLA(J)
      RW(J)=ASEC/SCW
      IF (IV .LT. IW) GO TO 48
      BTA(J)=-DMDX(J)*DSU/WMN(J)/TTR(J)/SCW/CFI
      IF (J.EQ.1 .OR. J.GT.KO .OR. IHT.EQ.0) GO TO 48
      IF (MOD(J,IHT) .NE. 1) GO TO 48
      IJ=J
      HAIR=RHOE*VE*CF*CHAIR
      CALL HEAT
48    IF (LD.LT.0) GO TO 56
      IF ((IV.GT.1).AND.(IV.LT.IW)) GO TO 56
      CFIK=2000.D+0*CFI
      CFK=2000.D+0*CF
      CFSK=2000.D+0*CFS
      DTHK=1000.D+0*DTH
      CTH=TWO*TH/(ONE+DSQRT(ONE-TWO*TH*ID/ASEC))
      CH=CDS(J)/CTH
      IEO=REO+HALF
      ITHX=RTHX+HALF
      WRITE (2,83) J,TW,TE,TAW,TP,IEO,ITHX,FRD,CFIK,CFK,CFSK,H,HI,FMY,DT
     1HK,TH,DELTA,DELST
      IF (J.LT.KBL-3) GO TO 54
      IF (J-KBL+2) 49,50,51
49    CTHA=CTH
      XNA=XN
      DLTA=DELTA
      REOA=REO
      GO TO 55
50    CTHB=CTH
      XNB=XN
      DLTB=DELTA
      REOB=REO
      GO TO 55
51    IF (J-KBL) 52,53,54
52    CTHC=CTH
      XNC=XN
      DLTC=DELTA
      REOC=REO
      GO TO 55
53    IF (IT.GT.0) GO TO 55
      DLST=GMA*CDS(J-3)+GMB*CDS(J-2)+GMC*CDS(J-1)+GMD*CDS(J)
      THBL=GMA*CTHA+GMB*CTHB+GMC*CTHC+GMD*CTH
      HBL=DLST/THBL
      DLTBL=GMA*DLTA+GMB*DLTB+GMC*DLTC+GMD*DELTA
      REOBL=GMA*REOA+GMB*REOB+GMC*REOC+GMD*REO
      REOFT=TLV*REOBL
      RETH=THBL*REOBL
      REDL=DLTBL*REOBL
      RETHG=DLOG10(RETH)
      REDLG=DLOG10(REDL)
      XNBL=GMA*XNA+GMB*XNB+GMC*XNC+GMD*XN
      GO TO 55
54    IF ((J.GT.3) .AND. (MOD(J,10).NE.0)) GO TO 56
55    WRITE (2,86) S(J),DSU,THU,CTH,HU,H,CH,XN
56    CONTINUE
      RW(1)=RCO(1)
      CALL SCOND(S,DLA,DAX,KAT)
      DO 57 J=1,KAT
57    DRX(J)=WALTAN(J)+DAX(J)
      IF ((IT.GT.0) .OR. (LD.LT.0)) GO TO 58
      IF ((IV.GT.1) .AND. (IV.LT.IW)) GO TO 58
      IF (KBL.LE.KAT) WRITE (2,85) XBL,DLST,THBL,HBL,XNBL,DLTBL,REOFT,RE
     1TH,RETHG,REDL,REDLG
      IF (KBL.LE.KAT) GO TO 58
      HBL=CDS(KAT)/CTH
      REOFT=TLV*REO
      RETH=CTH*REO
      REDL=DELTA*REO
      RETHG=DLOG10(RETH)
      REDLG=DLOG10(REDL)
      IF (KBL.GT.KAT) WRITE (2,85) S(KAT),CDS(KAT),CTH,HBL,XN,DELTA,REOF
     1T,RETH,RETHG,REDL,REDLG
58    CONTINUE
      DD2=BMA*DLA(1)+BMB*DLA(2)+BMC*DLA(3)
      DR2=SD(1)+DD2
      DXS=DAX(1)/DR2
      XST=S(1)-DXS
      YST=RCO(1)-HALF*DAX(1)**2/DR2
      SCW=DSQRT(ONE+DAX(1)**2)
      DR2=DR2/SCW**3
      RCV=ONE/DR2/YST
      IF (IT.GT.0) XBIN=SBIN-XST
      IF (IT.GT.0) XCIN=SCIN-XST
      WRITE (2,78) ITLE,XBIN,XCIN,SF
      PPQ=ZRO
      WRITE (2,67) RC,ETAD,AMACH,BMACH,CMACH,EMACH,MC,AH
      IF (TWT.NE.ZRO) GO TO 59
      WRITE (2,81) PPS,T0
      GO TO 5
59    WRITE (2,79) PPS,T0,TWT,TWAT,TAIR
      GO TO 5
60    IF (IT.EQ.0) GO TO 63
      DO 61 K=1,KAT
      S(K)=SREF(K)-XST
61    SCV(K)=DSQRT(ONE+DRX(K)**2)
      SCV(1)=ONE
      SL(1)=ZRO
      IM=(KAT-1)/2
      DO 62 I=1,IM
      J=2*I
      SS=S(J)-S(J-1)
      IF (I.EQ.1) SS=S(2)
      TT=S(J+1)-S(J)
      ST=SS+TT
      S1=(TWO-TT/SS)*ST/SIX
      S3=(TWO-SS/TT)*ST/SIX
      S2=ST-S1-S3
      SA=(TWO+TT/ST)*SS/SIX
      SB=(TWO+ST/TT)*SS/SIX
      SC=SS-SA-SB
      SL(J)=SL(J-1)+SA*SCV(J-1)+SB*SCV(J)+SC*SCV(J+1)
62    SL(J+1)=SL(J-1)+S1*SCV(J-1)+S2*SCV(J)+S3*SCV(J+1)
      XST=ZRO
      WRITE (2,68) LS,DK
      WRITE (2,69) (K,S(K),SL(K),DLA(K),RCO(K),WALTAN(K),SK(K),DAX(K),DR
     1X(K),WMN(K),DMDX(K),SPR(K),BTA(K),K=1,KAT)
      IF (KBL.GT.KAT) GO TO 64
      CALL TWIXT (SL,GMA,GMB,GMC,GMD,SBL,KAT,KBL)
      XBL=GMA*S(KBL-3)+GMB*S(KBL-2)+GMC*S(KBL-1)+GMD*S(KBL)
      DLAB=GMA*DLA(KBL-3)+GMB*DLA(KBL-2)+GMC*DLA(KBL-1)+GMD*DLA(KBL)
      RCOB=GMA*RCO(KBL-3)+GMB*RCO(KBL-2)+GMC*RCO(KBL-1)+GMD*RCO(KBL)
      WRITE (2,89) XBL,SBL,DLAB,RCOB,GMA,GMB,GMC,GMD
      GO TO 64
63    WRITE (2,68) LY,DD
      WRITE (2,69) (K,S(K),FS(K),DLA(K),RCO(K),WALTAN(K),SD(K),DAX(K),DR
     1X(K),WMN(K),DMDX(K),SPR(K),BTA(K),K=1,KAT)
      IF (KBL.GT.KAT) GO TO 64
      CALL TWIXT (S,GMA,GMB,GMC,GMD,XBL,KAT,KBL)
      DLAB=GMA*DLA(KBL-3)+GMB*DLA(KBL-2)+GMC*DLA(KBL-1)+GMD*DLA(KBL)
      RCOB=GMA*RCO(KBL-3)+GMB*RCO(KBL-2)+GMC*RCO(KBL-1)+GMD*RCO(KBL)
      YBL=RCOB-DLAB
      WRITE (2,84) XBL,YBL,DLAB,RCOB,GMA,GMB,GMC,GMD
64    WRITE (2,87) XST,YST,DD2,DR2,RCV
      S(1)=XST
      RCO(1)=YST
      DRX(1)=ZRO
      IF (SBL .EQ. 1.D+3) RETURN
      IF (LV .GT. 0) GO TO 3
65    CONTINUE
      IF (J .EQ. 1) WRITE (2,90) IPQ,QFUNW,TWT
      RETURN
C
66    FORMAT (6E10.3,4I5)
67    FORMAT (1H ,4H RC=,F11.6,3X,5HETAD=,F8.4,4H DEG,3X,6HAMACH=,F10.7,
     &3X,6HBMACH=,F10.7,3X,6HCMACH=,F10.7,3X,6HEMACH=,F10.7,3X,A4,2HH=,F
     &11.7/)
68    FORMAT (1H ,7X,9HSTA(IN)  ,A4,40H(IN)    DELR(IN)    R(IN)    DY/D
     &X     ,A8,50H     DA/DX     DR/DX    MACH NO.    DM/DX    PE/PO,7X
     &,4HBETA/)
69    FORMAT (10(I4,0P2F11.6,2F11.7,4F10.7,F11.7,F10.7,1P2E12.4/))
70    FORMAT (1H+,5X,34HQUADRATIC TEMPERATURE DISTRIBUTION)
71    FORMAT (1H+,5X,34HPARABOLIC TEMPERATURE DISTRIBUTION)
72    FORMAT (1H+,44X,34HSPALDING-CHI REFERENCE TEMPERATURE)
73    FORMAT (1H+,83X,36HVAN DRIEST REFERENCE REYNOLDS NUMBER /)
74    FORMAT (1H+,83X,35HCOLES LAW REFERENCE REYNOLDS NUMBER /)
75    FORMAT (1H+,44X,34HMODIF. SPALDING-CHI REFERENCE TEMP)
76    FORMAT (1H+,83X,40HREFERENCE REYNOLDS NUMBER BASED ON DELTA /)
77    FORMAT (1H+,44X,29HMODIFIED COLES TRANSFORMATION /)
78    FORMAT (1H1,3A4,39HNOZZLE CONTOUR, RADIAL FLOW ENDS AT STA,F12.7,2
     &5H, TEST CONE BEGINS AT STA,F12.7,16H, SCALE FACTOR =,F13.8/)
79    FORMAT (1H ,1X,15HSTAG. PRESSURE=,F5.0,24H PSI, STAG. TEMPERATURE=
     &,F5.0,21H DEG R, THROAT TEMP.=,F5.0,19H DEG R, WALL TEMP.=,F4.0,24
     &H DEG R, THROAT HT COEF.=,F8.5//)
80    FORMAT (1H1,3A4,49HBOUNDARY LAYER CALCULATIONS, STAGNATION PRESSUR
     &E=,F5.0,28HPSI, STAGNATION TEMPERATURE=,F5.0,27H DEG R, N BASED ON
     & RE,DELTA //)
81    FORMAT (1H ,5X,15HSTAG. PRESSURE=,F5.0,24H PSI STAG. TEMPERATURE=,
     &F5.0,34H DEG R ADIABATIC WALL TEMPERATURE//)
82    FORMAT (1H ,5X,38HTW    TE    TAW    TP    RE/IN    RTHI,4X,3HFRD,
     &5X,4HKCF1,4X,3HKCF,5X,4HRCFS,5X,1HH,6X,2HHI,5X,38HFMY     KTHP THE
     &TA-1  DELTA  DELTA*-1 /)
83    FORMAT (1H ,I3,2F6.1,F7.1,F6.1,I9,I7,4F8.5,F8.4,F7.4,2F8.5,F9.6,F7
     &.4,F9.6)
84    FORMAT (1H ,3HSTA,2F11.6,2F11.7,7X,27HINTERPOLATION COEFFICIENTS,,
     & F12.8,1H,,F11.8,1H,,F11.8,1H,,F12.8 /)
85    FORMAT (1H0,5H   X=,F7.3,11H,   DELTA*=,F10.7,10H,   THETA=,F9.7,6
     &H,   H=,F10.6,6H,   N=,F10.7,10H,   DELTA=,F11.7,10H,   RE/FT=,F11
     &.0//35X,9HRE,THETA=,F9.0,8H,   LOG=,F8.5,1H,,16X,9HRE,DELTA=,F11.0
     &,8H,   LOG=,F8.5)
86    FORMAT (1H ,3X,2HX=,F7.3,8H,   DSU=,F8.5,8H,   THU=,F9.7,8H,   CTH
     &=,F9.7,7H,   HU=,F10.6,6H,   H=,F10.6,7H,   CH=,F10.6,6H,   N=,F8.
     &5)
87    FORMAT (1H ,3HSTA,F11.6,9H      Y*=,F11.7,14H,     D2A/DX2=,F12.9,
     &14H,     D2R/DX2=,F12.9,16H,     VISCID RC=,F14.8)
88    FORMAT (1H ,5HRTHX=,1PE12.5,6H, REO=,E12.5,6H, FRO=,0PF8.5,5H, TH=
     &,F8.5,8H, DELTA=,F8.5,6H, DOT=,F9.5)
89    FORMAT (1H ,3HSTA,2F11.6,2F11.7,7X,27HINTERPOLATION COEFFICIENTS,,
     &F12.8,1H,,F11.8,1H,,F11.8,1H,,F12.8 /)
90    FORMAT (1H0,10H ITERATION,I4,11H,    QFUN =,F8.5,18H,   THROAT TEM
     &P  =,F6.1 /)
      END

max_stars_repo_name: aldorona/contur

max_issues_repo_name: aldorona/contur

max_forks_repo_name: aldorona/contur


text:
      MODULE LINE
          USE DISPLAY

          IMPLICIT NONE
          PRIVATE
          PUBLIC :: DRAWLINE

          INTEGER :: DX, DY, X, Y
      CONTAINS
          RECURSIVE SUBROUTINE DRAWLINE(DISPLAY, P1, P2, COLOR)
              INTEGER, INTENT(INOUT), DIMENSION(:, :, :, :) :: DISPLAY
              INTEGER, INTENT(IN), DIMENSION(3) :: P1, P2
              INTEGER, INTENT(IN), DIMENSION(3) :: COLOR

              DX = P2(1) - P1(1)
              DY = P2(2) - P1(2)
              X = P1(1)
              Y = P1(2)

              IF (ABS(DY) < ABS(DX)) THEN
                  IF (P1(1) > P2(1)) THEN
                      CALL DRAWLINE(DISPLAY, P2, P1, COLOR)
                  ELSE
                      CALL DRAWLINELOW(DISPLAY, P2(1), COLOR)
                  ENDIF
              ELSE
                  IF (P1(2) > P2(2)) THEN
                      CALL DRAWLINE(DISPLAY, P2, P1, COLOR)
                  ELSE
                      CALL DRAWLINEHIGH(DISPLAY, P2(2), COLOR)
                  ENDIF
              ENDIF
          END SUBROUTINE DRAWLINE

          SUBROUTINE DRAWLINELOW(DISPLAY, MX, COLOR)
              INTEGER, INTENT(INOUT), DIMENSION(:, :, :, :) :: DISPLAY
              INTEGER, INTENT(IN)               :: MX
              INTEGER, INTENT(IN), DIMENSION(3) :: COLOR
              INTEGER                           :: D, I, YI

              YI = 1
              IF (DY < 0) THEN
                  YI = -1
                  DY = -DY
              ENDIF

              D = 2 * DY - DX
              DO I = X, MX
                  CALL PLOT(DISPLAY, I, Y, 1, COLOR)
                  D = D + 2 * DY
                  IF (D > 0) THEN
                      Y = Y + YI
                      D = D - 2 * DX
                  END IF
              END DO
          END SUBROUTINE DRAWLINELOW

          SUBROUTINE DRAWLINEHIGH(DISPLAY, MY, COLOR)
              INTEGER, INTENT(INOUT), DIMENSION(:, :, :, :) :: DISPLAY
              INTEGER, INTENT(IN)               :: MY
              INTEGER, INTENT(IN), DIMENSION(3) :: COLOR
              INTEGER                           :: D, I, XI

              XI = 1
              IF (DX < 0) THEN
                  XI = -1
                  DX = -DX
              ENDIF

              D = 2 * DX - DY
              DO I = Y, MY
                  CALL PLOT(DISPLAY, X, I, 1, COLOR)
                  D = D + 2 * DX
                  IF (D > 0) THEN
                      X = X + XI
                      D = D - 2 * DY
                  END IF
              END DO
          END SUBROUTINE DRAWLINEHIGH
      END MODULE LINE

ext: for


text:
Name: Myra

20051107 16:45:56 nbsp Im impressed that you dont own a cell phone. How cool! Users/JohnDudek

20051109 16:42:45 nbsp Thanks... moral highground over all the cell phone snobs! Users/MyaBrn

20051116 22:09:46 nbsp you could be my dream girl except your married damn. i have a cell only for emergencies and no regular phone i communicate almost exclusivly in person. Users/AnonEmuss

20051121 15:04:27 nbsp Im a fulltime student, only work on campus part time, but I hear a lot about the lack of a pay raise from coworkers and service workers in the VMTH. And yeah I did get a little carried away, eh? heh.  Sorry! Didnt mean to make it feel like I was trying to come down on you or anything. Fees were hard enough before, I just got frustrated and ranted randomly. Users/EdWins ES

20060130 16:22:44 nbsp With an acorn is a mighty oak grown... that Trailer Park page has led to at least three other pages (since it also begged the question of what is that complex between the Arbors and the Trailer Park on Olive).  Cool first entry! Users/JabberWokky

20070119 19:04:57 nbsp I added as much as I can about Grannys store (I found her eBay sellers site), but if you could add a sentence or three about Granny herself it would probably finish off the entry.  You noted that she was nice.  Any other impressions? Users/JabberWokky

20070215 09:18:41 nbsp I rear ended someone once and it wasnt my fault either.  Despite that I always hear if you rear end somebody it is automatically your fault.  In my case the guy had swerved and then slammed on his brakes.  When the police officer asked why he braked, he apparantly said for the stop sign.  The officer looked at him, asked what stop sign?, and pointed down the busy, but sign free street. Users/JabberWokky

20070215 09:21:44 nbsp Yeah it was definitly my fault, I did not break in time. Jabberwokky I will also add more comments about Grannys. What happened to your page? I cant leave a comment Users/MyaBrn

20070416 14:37:10 nbsp Go to The UCD Annual Summer Job Fair Wednesday April 18 at Freeborn 10am2pm! Users/MyaBrn

20070508 11:59:53 nbsp Shouldnt you be on Bike Commute Month this list somewhere?  :) Users/KevinChin

20080215 06:15:36 nbsp everything wants to think they had at least a small effect on their successors Users/fredchen

20080216 11:49:40 nbsp did you read it by itself?

it was in response to your comment Users/fredchen

20080805 18:30:25 nbsp Youre into education? Me too! Your Nestware comment has some rather banal grammar errors. That is, its is a  pronoun, but its means it is; and alot is two words, as in a lot of people cant speak English. Users/onnaca

20090509 17:19:30 nbsp Thanks :) Users/EdWins

20110412 22:40:27 nbsp back on davis wiki Users/GQuincy

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


module trapezoidal_channel

    use define_settings, only: setting
    use define_globals
    use define_indexes
    use define_keys

    implicit none

    !%-----------------------------------------------------------------------------
    !% Description:
    !% Trapezoidal channel geometry
    !%

    private

    public :: trapezoidal_depth_from_volume
    public :: trapezoidal_area_from_depth
    public :: trapezoidal_area_from_depth_singular
    public :: trapezoidal_topwidth_from_depth
    public :: trapezoidal_topwidth_from_depth_singular
    public :: trapezoidal_perimeter_from_depth
    public :: trapezoidal_perimeter_from_depth_singular
    public :: trapezoidal_hyddepth_from_depth
    public :: trapezoidal_hyddepth_from_depth_singular
    public :: trapezoidal_hydradius_from_depth_singular

    contains
!%==========================================================================
!% PUBLIC
!%==========================================================================
!%
    subroutine trapezoidal_depth_from_volume (elemPGx, Npack, thisCol)
        !%-----------------------------------------------------------------------------
        !% Description:
        !% Only applies on open channels (or non-surcharged trapezoidal conduits)
        !% Input elemPGx is pointer (already assigned) for elemPGalltm, elemPGetm or elemPGac
        !% Assumes that volume > 0 is enforced in volume computations.
        !% NOTE: this does NOT limit the depth by surcharge height at this point
        !% This will be done after the head is computed.
        !%-----------------------------------------------------------------------------
        integer, target, intent(in) :: elemPGx(:,:), Npack, thisCol
        integer, pointer :: thisP(:)
        real(8), pointer :: depth(:), volume(:), length(:), breadth(:)
        real(8), pointer :: lslope(:), rslope(:)
        !%-----------------------------------------------------------------------------
        thisP   => elemPGx(1:Npack,thisCol)
        depth   => elemR(:,er_Depth)
        volume  => elemR(:,er_Volume)
        length  => elemR(:,er_Length)
        breadth => elemSGR(:,esgr_Trapezoidal_Breadth)
        lslope  => elemSGR(:,esgr_Trapezoidal_LeftSlope)
        rslope  => elemSGR(:,esgr_Trapezoidal_RightSlope)
        !%-----------------------------------------------------------------------------

        depth(thisP)       = - onehalfR * (breadth(thisP)/(onehalfR*(lslope(thisP) + rslope(thisP))) &
                - sqrt((breadth(thisP)/(onehalfR*(lslope(thisP) + rslope(thisP)))) ** twoR &
                + fourR * volume(thisP)/(onehalfR*length(thisP)*(lslope(thisP) + rslope(thisP)))))

    end subroutine trapezoidal_depth_from_volume
!%
!%==========================================================================
!%==========================================================================
!%
    elemental real(8) function trapezoidal_area_from_depth (indx) result (outvalue)
        !%-----------------------------------------------------------------------------
        !% Description:
        !% Computes area from known depth for trapezoidal cross section
        !% input may be a scalar or a packed array of indexes
        !%-----------------------------------------------------------------------------
        integer, intent(in) :: indx
        !%-----------------------------------------------------------------------------
        outvalue = (elemSGR(indx,esgr_Trapezoidal_Breadth)                        &
                        + onehalfR * (  elemSGR(indx,esgr_Trapezoidal_LeftSlope)  &
                                      + elemSGR(indx,esgr_Trapezoidal_RightSlope) &
                                     ) * elemR(indx,er_Depth)                     &
                    ) * elemR(indx,er_Depth)

    end function trapezoidal_area_from_depth
!%
!%==========================================================================
!%==========================================================================
!%
    real(8) function trapezoidal_area_from_depth_singular (indx) result (outvalue)
        !%-----------------------------------------------------------------------------
        !% Description:
        !% Computes area from known depth for trapezoidal cross section of a single element
        !% The input indx is the row index in full data 2D array.
        !%-----------------------------------------------------------------------------
        integer, intent(in) :: indx
        real(8), pointer :: depth(:), breadth(:), lslope(:), rslope(:)
        !%-----------------------------------------------------------------------------
        depth   => elemR(:,er_Depth)
        breadth => elemSGR(:,esgr_Trapezoidal_Breadth)
        lslope  => elemSGR(:,esgr_Trapezoidal_LeftSlope)
        rslope  => elemSGR(:,esgr_Trapezoidal_RightSlope)
        !%-----------------------------------------------------------------------------
        outvalue = (breadth(indx) + onehalfR * (lslope(indx) + rslope(indx)) * depth(indx)) * depth(indx)

    end function trapezoidal_area_from_depth_singular
!%
!%==========================================================================
!%==========================================================================
!%
    subroutine trapezoidal_topwidth_from_depth (elemPGx, Npack, thisCol)
        !%
        !%-----------------------------------------------------------------------------
        !% Description:
        !% Computes the topwidth from a known depth in a rectangular channel
        !%-----------------------------------------------------------------------------
        integer, target, intent(in) :: elemPGx(:,:)
        integer, intent(in) ::  Npack, thisCol
        integer, pointer :: thisP(:)
        real(8), pointer :: topwidth(:), depth(:), breadth(:), lslope(:), rslope(:)
        !%-----------------------------------------------------------------------------
        thisP    => elemPGx(1:Npack,thisCol)
        topwidth => elemR(:,er_Topwidth)
        depth    => elemR(:,er_Depth)
        breadth  => elemSGR(:,esgr_Trapezoidal_Breadth)
        lslope   => elemSGR(:,esgr_Trapezoidal_LeftSlope)
        rslope   => elemSGR(:,esgr_Trapezoidal_RightSlope)
        !%-----------------------------------------------------------------------------

        topwidth(thisP) = breadth(thisP) + depth(thisP) * (lslope(thisP) + rslope(thisP))

    end subroutine trapezoidal_topwidth_from_depth
!%
!%==========================================================================
!%==========================================================================
!%
    real(8) function trapezoidal_topwidth_from_depth_singular (indx) result (outvalue)
        !%-----------------------------------------------------------------------------
        !% Description:
        !% Computes the topwidth for a trapezoidal cross section of a single element
        !%-----------------------------------------------------------------------------
        integer, intent(in) :: indx
        real(8), pointer :: depth(:), breadth(:), lslope(:), rslope(:)
        !%-----------------------------------------------------------------------------
        depth   => elemR(:,er_Depth)
        breadth => elemSGR(:,esgr_Trapezoidal_Breadth)
        lslope  => elemSGR(:,esgr_Trapezoidal_LeftSlope)
        rslope  => elemSGR(:,esgr_Trapezoidal_RightSlope)
        !%-----------------------------------------------------------------------------
        outvalue = breadth(indx) + depth(indx) * (lslope(indx) + rslope(indx))

    end function trapezoidal_topwidth_from_depth_singular
!%
!%==========================================================================
!%==========================================================================
!%
    subroutine trapezoidal_perimeter_from_depth (elemPGx, Npack, thisCol)
        !%
        !%-----------------------------------------------------------------------------
        !% Description:
        !% Computes the perimeter from a known depth in a trapezoidal channel
        !%-----------------------------------------------------------------------------
        integer, target, intent(in) :: elemPGx(:,:)
        integer, intent(in) ::  Npack, thisCol
        integer, pointer :: thisP(:)
        real(8), pointer :: perimeter(:), depth(:), breadth(:), lslope(:), rslope(:)
        !%-----------------------------------------------------------------------------
        thisP     => elemPGx(1:Npack,thisCol)
        perimeter => elemR(:,er_Perimeter)
        depth     => elemR(:,er_Depth)
        breadth   => elemSGR(:,esgr_Trapezoidal_Breadth)
        lslope    => elemSGR(:,esgr_Trapezoidal_LeftSlope)
        rslope    => elemSGR(:,esgr_Trapezoidal_RightSlope)
        !%-----------------------------------------------------------------------------

        perimeter(thisP) = breadth(thisP) + depth(thisP) * (sqrt(oneR + lslope(thisP)**twoR) &
                        + sqrt(oneR + rslope(thisP)**twoR))

    end subroutine trapezoidal_perimeter_from_depth
!%
!%==========================================================================
!%==========================================================================
!%
    real(8) function trapezoidal_perimeter_from_depth_singular (indx) result (outvalue)
        !%
        !%-----------------------------------------------------------------------------
        !% Description:
        !% Computes wetted perimeter from known depth for a trapezoidal cross section of
        !% a single element
        !%-----------------------------------------------------------------------------
        !%-----------------------------------------------------------------------------
        integer, intent(in) :: indx
        real(8), pointer :: depth(:), breadth(:), lslope(:), rslope(:)
        !%-----------------------------------------------------------------------------
        depth   => elemR(:,er_Depth)
        breadth => elemSGR(:,esgr_Trapezoidal_Breadth)
        lslope  => elemSGR(:,esgr_Trapezoidal_LeftSlope)
        rslope  => elemSGR(:,esgr_Trapezoidal_RightSlope)
        !%-----------------------------------------------------------------------------

        outvalue =  breadth(indx) + depth(indx) * (sqrt(oneR + lslope(indx)**twoR) + &
                    sqrt(oneR + rslope(indx)**twoR))

    end function trapezoidal_perimeter_from_depth_singular
!%
!%==========================================================================
!%==========================================================================
!%
    subroutine trapezoidal_hyddepth_from_depth (elemPGx, Npack, thisCol)
        !%
        !%-----------------------------------------------------------------------------
        !% Description:
        !% Computes the hydraulic (average) depth from a known depth in a rectangular channel
        !%-----------------------------------------------------------------------------
        integer, target, intent(in) :: elemPGx(:,:)
        integer, intent(in) ::  Npack, thisCol
        integer, pointer :: thisP(:)
        real(8), pointer :: hyddepth(:), depth(:), breadth(:), lslope(:), rslope(:)
        !%-----------------------------------------------------------------------------
        thisP     => elemPGx(1:Npack,thisCol)
        hyddepth  => elemR(:,er_HydDepth)
        depth   => elemR(:,er_Depth)
        breadth => elemSGR(:,esgr_Trapezoidal_Breadth)
        lslope  => elemSGR(:,esgr_Trapezoidal_LeftSlope)
        rslope  => elemSGR(:,esgr_Trapezoidal_RightSlope)
        !%-----------------------------------------------------------------------------

        hyddepth(thisP) = ((breadth(thisP) + onehalfR * (lslope(thisP) + rslope(thisP)) * &
                    depth(thisP)) * depth(thisP)) / (breadth(thisP) + depth(thisP) * &
                    (lslope(thisP) + rslope(thisP)))

    end subroutine trapezoidal_hyddepth_from_depth
!%
!%==========================================================================
!%==========================================================================
!%
    real(8) function trapezoidal_hyddepth_from_depth_singular (indx) result (outvalue)
        !%
        !%-----------------------------------------------------------------------------
        !% Description:
        !% Computes hydraulic depth from known depth for trapezoidal cross section of
        !% a single element
        !%-----------------------------------------------------------------------------
        integer, intent(in) :: indx
        real(8), pointer :: depth(:), breadth(:), lslope(:), rslope(:)
        !%-----------------------------------------------------------------------------
        depth   => elemR(:,er_Depth)
        breadth => elemSGR(:,esgr_Trapezoidal_Breadth)
        lslope  => elemSGR(:,esgr_Trapezoidal_LeftSlope)
        rslope  => elemSGR(:,esgr_Trapezoidal_RightSlope)
        !%-----------------------------------------------------------------------------

        outvalue = ((breadth(indx) + onehalfR * (lslope(indx) + rslope(indx)) * &
                    depth(indx)) * depth(indx)) / (breadth(indx) + depth(indx) * &
                    (lslope(indx) + rslope(indx)))

    end function trapezoidal_hyddepth_from_depth_singular
!%
!%==========================================================================
!%==========================================================================
!%
    real(8) function trapezoidal_hydradius_from_depth_singular (indx) result (outvalue)
        !%
        !%-----------------------------------------------------------------------------
        !% Description:
        !% Computes hydraulic radius from known depth for a trapezoidal cross section of
        !% a single element
        !%-----------------------------------------------------------------------------
        integer, intent(in) :: indx
        real(8), pointer :: depth(:), breadth(:), lslope(:), rslope(:)
        !%-----------------------------------------------------------------------------
        depth   => elemR(:,er_Depth)
        breadth => elemSGR(:,esgr_Trapezoidal_Breadth)
        lslope  => elemSGR(:,esgr_Trapezoidal_LeftSlope)
        rslope  => elemSGR(:,esgr_Trapezoidal_RightSlope)
        !%-----------------------------------------------------------------------------

        outvalue = ((breadth(indx) + onehalfR * (lslope(indx) + rslope(indx)) *  &
                    depth(indx)) * depth(indx)) / (breadth(indx) + depth(indx) * &
                    (sqrt(oneR + lslope(indx)**twoR) + sqrt(oneR + rslope(indx)**twoR)))

    end function trapezoidal_hydradius_from_depth_singular
!%
!%==========================================================================

!%
!%
!%==========================================================================
!%==========================================================================
!%
        !%
        !%-----------------------------------------------------------------------------
        !% Description:
        !%
        !%-----------------------------------------------------------------------------

        !%-----------------------------------------------------------------------------
        !%
!%
!%
!%==========================================================================
!%==========================================================================
!%
        !%
        !%-----------------------------------------------------------------------------
        !% Description:
        !%
        !%-----------------------------------------------------------------------------

        !%-----------------------------------------------------------------------------
        !%
!%
!%
!%==========================================================================
!%==========================================================================
!%
        !%
        !%------------------------------------------------------------------
        !% Description:
        !%
        !%------------------------------------------------------------------

        !%------------------------------------------------------------------
        !%
end module trapezoidal_channel