text:
*+CEL_CARTS
        subroutine CEL_CARTS( ra, dec, vector )

        implicit none

        double precision deg2rad
        parameter( deg2rad = 0.01745329252D+00 )

        double precision ra, dec
        double precision vector( 3 )

*       Description:
*         Given RA and Dec, returns the 3-d vector
*
*       Arguments:
*         ra        (i) : input RA (decimal hours)
*         dec       (i) : input Dec
*         vector    (o) : 3-d vector
*
*       Dependencies:
*         Assumes the presence of trig functions on degrees on the system
*
*       Origin:
*         Conceived as a part of the 3-d vector routines
*
*       Author:
*         Koji Mukai, 1993 Mar 15, Original version
*-CEL_CARTS

        double precision cos_dec, ra_rad, dec_rad

	dec_rad = dec * deg2rad
        cos_dec = cos( dec_rad )
        ra_rad = ra * 1.5d+01 * deg2rad
        vector( 1 ) = cos_dec * cos( ra_rad )
        vector( 2 ) = cos_dec * sin( ra_rad )
        vector( 3 ) = sin( dec_rad )

        end

*+ECL_CARTS
        subroutine ECL_CARTS( lambda, beta, vector )

        implicit none

        double precision deg2rad
        parameter( deg2rad = 0.01745329252D+00 )

        double precision lambda, beta
        double precision vector( 3 )

*       Description:
*         Given ecliptic longitude and latitude, returns the 3-d vector
*
*       Arguments:
*         ra        (i) : input RA (decimal hours)
*         dec       (i) : input Dec
*         vector    (o) : 3-d vector
*
*       Dependencies:
*         Assumes the presence of trig functions on degrees on the system
*
*       Origin:
*         Conceived as a part of the 3-d vector routines
*
*       Author:
*         Koji Mukai, 1993 Mar 15, Original version
*-CEL_CARTS

        double precision cos_beta, lambda_rad, beta_rad

	beta_rad = beta * deg2rad
        cos_beta = cos( beta_rad )
        lambda_rad = lambda * deg2rad
        vector( 1 ) = cos_beta * cos( lambda_rad )
        vector( 2 ) = cos_beta * sin( lambda_rad )
        vector( 3 ) = sin( beta_rad )

        end

*+VEC_ANGLE
        double precision function VEC_ANGLE( vector1, vector2 )

        implicit none

        double precision rad2deg
	parameter( rad2deg = 57.29577951D+00 )

        double precision vector1( 3 ), vector2( 3 )

*       Description:
*         Angle between two cartesian vectors
*
*       Arguments:
*         vector1, vector2 (i) : two input Cartesian vectors
*         <VEC_ANGLE>      (r) : Angle between the two
*
*       Dependencies:
*         Assumes the presence of trig functions on degrees on the system
*
*       Origin:
*         Conceived as a part of the 3-d vector routines
*
*       Author:
*         Koji Mukai, 1993 Mar 15, Original version
*-VEC_ANGLE

        double precision inner, temp, work
        integer j

        inner = 0.0d+00
        do j = 1, 3
          inner = inner + vector1( j ) * vector2( j )
        end do

        if( inner .gt. 0.71d+00 ) then
*            Use sine, which is more accurate in this regime
*            Length( vector1 - vector2 ) * 0.5 = sin( angle * 0.5 )
*            (Cosine of small x is ~1.0-x^2/2, whereas sin x ~ x)
          temp = 0.0
          do j = 1, 3
            work = vector1( j ) - vector2( j )
            temp = temp + work * work
          end do
          temp = sqrt( temp ) * 0.5d+00
          VEC_ANGLE = asin( temp ) * 2.0d+00 * rad2deg
        else if( inner .ge. -0.71d+00 ) then
*            Use the more familiar cosine formula
          VEC_ANGLE = acos( inner ) * rad2deg
        else
*            Back to sine, but two vectors are ~180 degrees apart
          temp = 0.0
          do j = 1, 3
            work = vector1( j ) + vector2( j )
            temp = temp + work * work
          end do
          temp = sqrt( temp ) * 0.5d+00
          VEC_ANGLE = 1.8d+02 - asin( temp ) * 2.0d+00 * rad2deg
        end if

        end

*+VEC_PRDCT
        subroutine VEC_PRDCT( vector1, vector2, vector3 )

        implicit none

        double precision vector1( 3 ), vector2( 3 ), vector3( 3 )

*       Description:
*         Calculates product of two vectors
*
*       Arguments:
*         vector1, vector2 (i) : two input Cartesian vectors
*         vector3          (o) : product of the two input vectors
*
*       Dependencies:
*         None
*
*       Origin:
*         Conceived as a part of the 3-d vector routines
*
*       Author:
*         Koji Mukai, 1993 Mar 15, Original version
*-VEC_PRDCT

        integer i, j, k

        do i = 1, 3
          j = i + 1
          if( j .eq. 4 ) j = 1
          k = j + 1
          if( k .eq. 4 ) k = 1
          vector3( i ) = vector1( j ) * vector2( k )
     &                                - vector1( k ) * vector2( j )
        end do

        end

*+VEC_RENRM
        subroutine VEC_RENRM( vector )

        implicit none

        double precision vector( 3 )

*       Description:
*         Renormalizes a Cartesian vector
*
*       Arguments:
*         vector  (i/o) : Cartesian vector to be renormalized
*
*       Dependencies:
*         None
*
*       Origin:
*         Conceived as a part of the 3-d vector routines
*
*       Author:
*         Koji Mukai, 1993 Mar 15, Original version
*-VEC_RENRM

        double precision temp
        integer j

        temp = 0.0d+00
        do j = 1, 3
          temp = temp + vector( j ) * vector( j )
        end do
        temp = 1.0d+00 / sqrt( temp )
        do j = 1, 3
          vector( j ) = vector( j ) * temp
        end do

        end

max_stars_repo_name: barentsen/tvguide

max_issues_repo_name: barentsen/tvguide

max_forks_repo_name: barentsen/tvguide


      module sd_channel_module

      implicit none

      integer :: maxint                           !number of intervals in hydrograph for degredation
      real :: peakrate, sed_reduc_t, no3_reduc_kg, tp_reduc_kg, tp_reduc, srp_reduc_kg
      real, dimension(:), allocatable :: hyd_rad  !m^2        |hydraulic radius for each hydrograph time step
      real, dimension(:), allocatable :: timeint  !days       |time spent in each hydrograph time step

      type swatdeg_hydsed_data
        character(len=16) :: name
        character(len=16) :: order
        real :: chw             !m         |channel width
        real :: chd             !m         |channel depth
        real :: chs             !m/m       |channel slope
        real :: chl             !km        |channel length
        real :: chn             !          |channel Manning's n
        real :: chk             !mm/h      |channel bottom conductivity
        real :: cherod          !          |channel erodibility
        real :: cov             !0-1       |channel cover factor
        real :: hc_cov          !0-1       |head cut cover factor
        real :: chseq           !m/m       |equilibrium channel slope
        real :: d50             !mm        |channel median sediment size
        real :: clay            !%         |clay percent of bank and bed
        real :: carbon          !%         |cabon percent of bank and bed
        real :: bd              !t/m3      |dry bulk density
        real :: chss            !          |channel side slope
        real :: bedldcoef       !          |percent of sediment entering the channel that is bed material
        real :: tc              !          |time of concentration
        real :: shear_bnk       !0-1       |bank shear coefficient - fraction of bottom shear
        real :: hc_kh           !          |headcut erodibility
        real :: hc_hgt          !m         |headcut height
        real :: hc_ini          !km        |initial channel length for gullies
      end type swatdeg_hydsed_data
      type (swatdeg_hydsed_data), dimension (:), allocatable :: sd_chd

      type swatdeg_init_datafiles
        integer :: init = 1                 !initial data-points to initial.cha
        integer :: org_min = 1              !points to initial organic-mineral input file
        integer :: pest = 1                 !points to initial pesticide input file
        integer :: path = 1                 !points to initial pathogen input file
        integer :: hmet = 1                 !points to initial heavy metals input file
        integer :: salt = 1                 !points to initial salt input file
      end type swatdeg_init_datafiles
      type (swatdeg_init_datafiles), dimension(:), allocatable :: sd_init

      type swatdeg_datafiles
        character(len=16) :: name = ""
        character(len=16) :: initc = ""
        character(len=16) :: hydc = ""
        character(len=16) :: sedc = ""
        character(len=16) :: nutc = ""
        integer :: init = 1
        integer :: hyd = 1
        integer :: sed = 1
        integer :: nut = 1
      end type swatdeg_datafiles
      type (swatdeg_datafiles), dimension(:),allocatable :: sd_dat

      type swatdeg_channel_dynamic
        character(len=13) :: name = "default"
        integer :: props
        integer :: obj_no
        integer :: aqu_link = 0             !aquifer the channel is linked to
        integer :: aqu_link_ch = 0          !sequential channel number in the aquifer
        character(len=16) :: region
        character(len=16) :: order
        real :: chw = 3.        !m          |channel width
        real :: chd = .5        !m          |channel depth
        real :: chs = .01       !m/m        |channel slope
        real :: chl = .1        !km         |channel length
        real :: chn             !           |channel Manning's n
        real :: cov             !0-1        |channel cover factor
        real :: cherod          !           |channel erodibility
        real :: shear_bnk       !0-1        |bank shear coefficient - fraction of bottom shear
        real :: hc_erod         !           |headcut erodibility
        real :: hc_co = 0.      !m/m        |proportionality coefficient for head cut
        real :: hc_len = 0.     !m          |length of head cut
        real :: hc_hgt          !m          |headcut height
        real, dimension (:), allocatable :: kd      !           |aquatic mixing velocity (diffusion/dispersion)-using mol_wt
        real, dimension (:), allocatable :: aq_mix  ! m/day     |aquatic mixing velocity (diffusion/dispersion)-using mol_wt
        character (len=2) :: overbank               !           |"ib"=in bank; "ob"=overbank flood
      end type swatdeg_channel_dynamic
      type (swatdeg_channel_dynamic), dimension (:), allocatable :: sd_ch
      type (swatdeg_channel_dynamic), dimension (:), allocatable :: sdch_init

      type sd_ch_output
        real :: flo_in = 0.             ! (m^3/s)      !ave flow rate
        real :: aqu_in = 0.             ! (m^3/s)      !ave flow rate
        real :: flo = 0.                ! (m^3/s)      !ave flow rate
        real :: peakr = 0.              ! (m^3/s)      |peak runoff rate
        real :: sed_in = 0.             ! (tons)       !total sed in
        real :: sed_out = 0.            ! (tons)       !total sed out
        real :: washld = 0.             ! (tons)       !wash load
        real :: bedld = 0.              ! (tons)       !bed load
        real :: dep = 0.                ! (tons)       !deposition
        real :: deg_btm = 0.            ! (tons)       !bottom erosion
        real :: deg_bank = 0.           ! (tons)       !bank erosion
        real :: hc_sed = 0.             ! (tons)       !headcut erosion
        real :: width = 0.              !
        real :: depth = 0.              !
        real :: slope = 0.              !
        real :: deg_btm_m = 0.          ! (m)          !downcutting
        real :: deg_bank_m = 0.         ! (m)          !widening
        real :: hc_m = 0.               ! (m)          !headcut retreat
      end type sd_ch_output

      type (sd_ch_output), dimension(:), allocatable, save :: chsd_d
      type (sd_ch_output), dimension(:), allocatable, save :: chsd_m
      type (sd_ch_output), dimension(:), allocatable, save :: chsd_y
      type (sd_ch_output), dimension(:), allocatable, save :: chsd_a
      type (sd_ch_output), dimension(:), allocatable, save :: schsd_d
      type (sd_ch_output), dimension(:), allocatable, save :: schsd_m
      type (sd_ch_output), dimension(:), allocatable, save :: schsd_y
      type (sd_ch_output), dimension(:), allocatable, save :: schsd_a
      type (sd_ch_output) :: bchsd_d
      type (sd_ch_output) :: bchsd_m
      type (sd_ch_output) :: bchsd_y
      type (sd_ch_output) :: bchsd_a
      type (sd_ch_output) :: chsdz

      type sdch_header
          character (len=6) :: day        =  "  jday"
          character (len=6) :: mo         =  "   mon"
          character (len=6) :: day_mo     =  "   day"
          character (len=6) :: yrc        =  "    yr"
          character (len=8) :: isd        =  "   unit "
          character (len=8) :: id         =  " gis_id "
          character (len=16) :: name      =  " name              "
          character(len=16) :: flo_in     =  "         flo_in"        ! (m^3/s)
          character(len=16) :: aqu_in     =  "         aqu_in"        ! (m^3/s)
          character(len=16) :: flo        =  "         flo_out"       ! (m^3/s)
          character(len=15) :: peakr      =  "          peakr"        ! (m^3/s)
          character(len=15) :: sed_in     =  "         sed_in"        ! (tons)
          character(len=15) :: sed_out    =  "        sed_out"        ! (tons)
          character(len=15) :: washld     =  "         washld"        ! (tons)
          character(len=15) :: bedld      =  "          bedld"        ! (tons)
          character(len=15) :: dep        =  "            dep"        ! (tons)
          character(len=15) :: deg_btm    =  "        deg_btm"        ! (tons)
          character(len=15) :: deg_bank   =  "       deg_bank"        ! (tons)
          character(len=15) :: hc_sed     =  "         hc_sed"        ! (tons)
          character(len=15) :: width      =  "          width"        ! (m)
          character(len=15) :: depth      =  "          depth"        ! (m)
          character(len=15) :: slope      =  "          slope"        ! (m/m)
          character(len=15) :: deg_btm_m  =  "        deg_btm"        ! (m)
          character(len=15) :: deg_bank_m =  "       deg_bank"        ! (m)
          character(len=15) :: hc_len     =  "         hc_len"        ! (m)
      end type sdch_header
      type (sdch_header) :: sdch_hdr

     type sdch_header_units
          character (len=6) :: day        =  "      "
          character (len=6) :: mo         =  "      "
          character (len=6) :: day_mo     =  "      "
          character (len=6) :: yrc        =  "      "
          character (len=8) :: isd        =  "        "
          character (len=8) :: id         =  "        "
          character (len=16) :: name      =  "                   "
          character(len=16) :: flo_in     =  "           m^3/s"       ! (m^3/s)
          character(len=16) :: aqu_in     =  "           m^3/s"       ! (m^3/s)
          character(len=16) :: flo        =  "           m^3/s"       ! (m^3/s)
          character(len=15) :: peakr      =  "          m^3/s"        ! (m^3/s)
          character(len=15) :: sed_in     =  "           tons"        ! (tons)
          character(len=15) :: sed_out    =  "           tons"        ! (tons)
          character(len=15) :: washld     =  "           tons"        ! (tons)
          character(len=15) :: bedld      =  "           tons"        ! (tons)
          character(len=15) :: dep        =  "           tons"        ! (tons)
          character(len=15) :: deg_btm    =  "           tons"        ! (tons)
          character(len=15) :: deg_bank   =  "           tons"        ! (tons)
          character(len=15) :: hc_sed     =  "           tons"        ! (tons)
          character(len=15) :: width      =  "              m"        ! (m)
          character(len=15) :: depth      =  "              m"        ! (m)
          character(len=15) :: slope      =  "            m/m"        ! (m/m)
          character(len=15) :: deg_btm_m  =  "              m"        ! (m)
          character(len=15) :: deg_bank_m =  "              m"        ! (m)
          character(len=15) :: hc_len     =  "              m"        ! (m)
      end type sdch_header_units
      type (sdch_header_units) :: sdch_hdr_units

      interface operator (+)
        module procedure chsd_add
      end interface

      interface operator (/)
        module procedure chsd_div
      end interface

      interface operator (//)
        module procedure chsd_ave
      end interface

      interface operator (*)
        module procedure chsd_mult
      end interface

      contains
!! routines for swatdeg_hru module

      function chsd_add(cho1,cho2) result (cho3)
      type (sd_ch_output),  intent (in) :: cho1
      type (sd_ch_output),  intent (in) :: cho2
      type (sd_ch_output) :: cho3
       cho3%flo_in = cho1%flo_in + cho2%flo_in
       cho3%aqu_in = cho1%aqu_in + cho2%aqu_in
       cho3%flo = cho1%flo + cho2%flo
       cho3%peakr = cho1%peakr + cho2%peakr
       cho3%sed_in = cho1%sed_in + cho2%sed_in
       cho3%sed_out = cho1%sed_out + cho2%sed_out
       cho3%washld = cho1%washld + cho2%washld
       cho3%bedld = cho1%bedld + cho2%bedld
       cho3%dep = cho1%dep + cho2%dep
       cho3%deg_btm = cho1%deg_btm + cho2%deg_btm
       cho3%deg_bank = cho1%deg_bank + cho2%deg_bank
       cho3%hc_sed = cho1%hc_sed + cho2%hc_sed
       cho3%width = cho1%width + cho2%width
       cho3%depth = cho1%depth + cho2%depth
       cho3%slope = cho1%slope + cho2%slope
       cho3%deg_btm_m = cho1%deg_btm_m + cho2%deg_btm_m
       cho3%deg_bank_m = cho1%deg_bank_m + cho2%deg_bank_m
       cho3%hc_m = cho1%hc_m + cho2%hc_m
      end function

      function chsd_div (ch1,const) result (ch2)
        type (sd_ch_output), intent (in) :: ch1
        real, intent (in) :: const
        type (sd_ch_output) :: ch2
        ch2%flo_in = ch1%flo_in / const
        ch2%aqu_in = ch1%aqu_in / const
        ch2%flo = ch1%flo / const
        ch2%peakr = ch1%peakr
        ch2%sed_in = ch1%sed_in / const
        ch2%sed_out = ch1%sed_out / const
        ch2%washld = ch1%washld / const
        ch2%bedld = ch1%bedld / const
        ch2%dep = ch1%dep / const
        ch2%deg_btm = ch1%deg_btm / const
        ch2%deg_bank = ch1%deg_bank / const
        ch2%hc_sed = ch1%hc_sed / const
        ch2%width = ch1%width
        ch2%depth = ch1%depth
        ch2%slope = ch1%slope
        ch2%deg_btm_m = ch1%deg_btm_m
        ch2%deg_bank_m = ch1%deg_bank_m
        ch2%hc_m = ch1%hc_m
      end function chsd_div

      function chsd_ave (ch1,const) result (ch2)
        type (sd_ch_output), intent (in) :: ch1
        real, intent (in) :: const
        type (sd_ch_output) :: ch2
        ch2%flo_in = ch1%flo_in
        ch2%aqu_in = ch1%aqu_in
        ch2%flo = ch1%flo
        ch2%peakr = ch1%peakr / const
        ch2%sed_in = ch1%sed_in
        ch2%sed_out = ch1%sed_out
        ch2%washld = ch1%washld
        ch2%bedld = ch1%bedld
        ch2%dep = ch1%dep
        ch2%deg_btm = ch1%deg_btm
        ch2%deg_bank = ch1%deg_bank
        ch2%hc_sed = ch1%hc_sed
        ch2%width = ch1%width / const
        ch2%depth = ch1%depth / const
        ch2%slope = ch1%slope / const
        ch2%deg_btm_m = ch1%deg_btm_m / const
        ch2%deg_bank_m = ch1%deg_bank_m / const
        ch2%hc_m = ch1%hc_m / const
      end function chsd_ave

      function chsd_mult (const, chn1) result (chn2)
        type (sd_ch_output), intent (in) :: chn1
        real, intent (in) :: const
        type (sd_ch_output) :: chn2
        chn2%flo_in = const * chn1%flo_in
        chn2%aqu_in = const * chn1%aqu_in
        chn2%flo = const * chn1%flo
        chn2%peakr = const * chn1%peakr
        chn2%sed_in = const * chn1%sed_in
        chn2%sed_out = const * chn1%sed_out
        chn2%washld = const * chn1%washld
        chn2%bedld = const * chn1%bedld
        chn2%dep = const * chn1%dep
        chn2%deg_btm = const * chn1%deg_btm
        chn2%deg_bank = const * chn1%deg_bank
        chn2%hc_sed = const * chn1%hc_sed
        chn2%width = const * chn1%width
        chn2%depth = const * chn1%depth
        chn2%slope = const * chn1%slope
        chn2%deg_btm_m = const * chn1%deg_btm_m
        chn2%deg_bank_m = const * chn1%deg_bank_m
        chn2%hc_m = const * chn1%hc_m
      end function chsd_mult

      end module sd_channel_module


text:
      SUBROUTINE POIDS(X,N,W)
      IMPLICIT REAL*8(A-H,O-Z)
      REAL*8 l
      DIMENSION X(N)

      l=360.0D0
      un=1.0D0
      rho=0.1D0

      W=l*(un/X(1)+un/X(2)+un/X(4))
      W=W+l*DSQRT(2.0D0)*(un/X(3)+un/X(5)+un/X(6))
      W=W*rho

      WRITE(*,*)'Poids actuel = ',W,' lbs'

      RETURN
      END

ext: for

max_stars_repo_name: rboman/progs

max_issues_repo_name: rboman/progs

max_forks_repo_name: rboman/progs


text:
!>
module UIVortexRings
    use PrecisionMod, only: MK
    implicit none

    real(MK),parameter :: MINNORM=1e-4
contains

    !     subroutine fUi_VortexRings(CPs,PRingCenter,R,GammaOriented, SmoothModel,SmoothParam,bComputeGrad,Uind,Grad)
    !         use MathConstants, only: fourpi_inv
    !         ! Arguments declarations
    !         integer, intent(in) :: nCPs
    !         integer, intent(in) :: nRings
    !         logical, intent(in) :: bComputeGrad !<
    !         real(MK), dimension(3,nCPs), intent(in)    :: CPs         !< side effects
    !         real(MK), dimension(3,nCPs), intent(inout) :: Uind        !< side effects
    !         real(MK), dimension(9,nCPs), intent(inout) :: Grad        !< side effects
    !         real(MK), dimension(3,nRings), intent(in)  :: PRingCenter !<
    !         real(MK), dimension(nSrc), intent(in)      :: Intensities !<
    !         integer, intent(in)                        :: SmoothModel !<
    !         real(MK), dimension(nSrc), intent(in)      :: SmoothParam !<
    ! !         ! Variable declarations
    !     end subroutine
    !
    subroutine fUi_VortexRing_11(CP,PRingCenter,R,GammaOriented, SmoothModel,SmoothParam,bComputeGrad,Uind,Grad)
        ! Arguments declarations
        real(MK), dimension(3), intent(in)    :: CP           !<
        real(MK), dimension(3), intent(in)    :: PRingCenter  !< PRingCenter : Position of ring center
        real(MK), intent(in)                  :: R            !< Ring radius
        real(MK), dimension(3), intent(in)    :: GammaOriented!< Vector norm is ring intensity, direction is normal to the ring
        integer, intent(in)                   :: SmoothModel  !<
        real(MK), intent(in)                  :: SmoothParam  !<
        logical, intent(in)                   :: bComputeGrad !<
        real(MK), dimension(3), intent(inout) :: Uind         !< side effects
        real(MK), dimension(9), intent(inout) :: Grad         !< side effects
        ! Variable declarations
        real(MK), dimension(3) :: ex  !<
        real(MK), dimension(3) :: ey  !<
        real(MK), dimension(3) :: ez  !<
        real(MK), dimension(3,3) :: T_r2g
        real(MK), dimension(3,3) :: T_g2r
        real(MK) :: Gamma  !<
        integer :: imin
        real(MK), dimension(3) :: CP_loc  !<
        real(MK) :: R_cp
        real(MK) :: Z_cp
        real(MK) :: phi
        real(MK), dimension(3) :: Uindtmp         !<
        real(MK), dimension(9) :: Gradtmp         !<
        real(MK) :: urad
        real(MK) :: upsi
        !         ! Note: theories return only vr (index 1) and vz (index 3), there is some cartesian handling though, but does it work ? Is it consistent with the rest of the codes. Does it make use of the inputed vpsi to reduce the error in the angle estimate.
        !
        !! Setting up an orthonornal cartesian system normal to the ring
        Gamma=norm2(GammaOriented)
        ez(1:3)=GammaOriented(1:3)/Gamma
        imin=minloc(abs(ez),1)
        ey=ez
        ey(imin)=0
        ey=ey([ 2,3,1 ] )
        ey(imin)=-ey(imin)
        ey=ey/norm2(ey)
        ex=cross(ey,ez)
        !! Now setting up rotation matrix
        T_r2g(1:3,1)=ex(1:3)
        T_r2g(1:3,2)=ey(1:3)
        T_r2g(1:3,3)=ez(1:3)
        T_g2r=transpose(T_r2g)

        !! Coordinate transformation to a cartesian frame normal to the ring
        CP_loc(1)=T_g2r(1,1)*(CP(1)-PRingCenter(1)) + T_g2r(1,2)*(CP(2)-PRingCenter(2)) + T_g2r(1,3)*(CP(3)-PRingCenter(3))
        CP_loc(2)=T_g2r(2,1)*(CP(1)-PRingCenter(1)) + T_g2r(2,2)*(CP(2)-PRingCenter(2)) + T_g2r(2,3)*(CP(3)-PRingCenter(3))
        CP_loc(3)=T_g2r(3,1)*(CP(1)-PRingCenter(1)) + T_g2r(3,2)*(CP(2)-PRingCenter(2)) + T_g2r(3,3)*(CP(3)-PRingCenter(3))

        !! Coordinate transform to a canonical polar coordinates
        Z_cp=CP_loc(3) ! Since we are already in the ring coordinate system
        R_cp=sqrt(CP_loc(1)**2+CP_loc(2)**2)
        !print*,'CP_loc',CP_loc

        !! Computing velocity in polar coordinates
        Uindtmp=0.0_MK
        Gradtmp=0.0_MK
        if(SmoothModel==0) then
            call fUi_Ring11_loc(R_cp, Z_cp, Gamma, R, Uindtmp)
            !print*,'Uind pol',Uindtmp
        else
            print*,'Smooth param',SmoothParam ! for unused var compiler warning
            write(*,*)'Smooth Ring not done'
            STOP -1
        endif

        !! Converting to cartesian
        phi=atan2(CP_loc(2),CP_loc(1))
        !print*,'phi',phi*180/3.141593
        urad=Uindtmp(1)
        upsi=0.0_MK
        Uindtmp(1)=+urad*cos(phi)
        Uindtmp(2)=+urad*sin(phi)
        !print*,'Uind cart',Uindtmp

        !! Going back to global frame
        Uind(1)=Uind(1) + T_r2g(1,1)*Uindtmp(1) + T_r2g(1,2)*Uindtmp(2) + T_r2g(1,3)*Uindtmp(3)
        Uind(2)=Uind(2) + T_r2g(2,1)*Uindtmp(1) + T_r2g(2,2)*Uindtmp(2) + T_r2g(2,3)*Uindtmp(3)
        Uind(3)=Uind(3) + T_r2g(3,1)*Uindtmp(1) + T_r2g(3,2)*Uindtmp(2) + T_r2g(3,3)*Uindtmp(3)

        !print*,'Uind glob',Uind
        !! TODO gradient
        if (bComputeGrad) then
            Grad = 0.0D0
        else
            Grad = 0.0D0
        end if
    contains
        !> Cross product of two vectors ! See crossprod in UTILS...
        function cross(a, b)
            real(MK), dimension(3) :: cross
            real(MK), dimension(3), intent(in) :: a, b
            cross(1) = a(2) * b(3) - a(3) * b(2)
            cross(2) = a(3) * b(1) - a(1) * b(3)
            cross(3) = a(1) * b(2) - a(2) * b(1)
        end function cross
    end subroutine

    subroutine test_ring()
        ! Arguments declarations
        real(MK):: Gamma !< Ring intensity
        real(MK):: R     !< Ring radius
        real(MK), dimension(3) :: Uind  !< side effects, in coord (r, theta, z)
        real(MK), dimension(3) :: CP           !<
        real(MK), dimension(3):: PRingCenter  !< PRingCenter : Position of ring center
        real(MK), dimension(3):: GammaOriented!< Vector norm is ring intensity, direction is normal to the ring
        integer:: SmoothModel  !<
        real(MK):: SmoothParam  !<
        logical:: bComputeGrad !<
        real(MK), dimension(9):: Grad         !< side effects

        R=3._MK
        Gamma=2._MK
        Uind=0_MK
        call fUi_Ring11_loc(1._MK, 2._MK, Gamma, R, Uind)
        GammaOriented=[2,0,0];
        PRingCenter=[3,2,1];
        CP=[5,2,1]
        CP=CP+[0.,cos(0.1),sin(0.1)];

        SmoothModel=0
        SmoothParam=0.0
        bComputeGrad=.false.
        Uind=[0,0,0]
        Grad(1:9)=0
        call fUi_VortexRing_11(CP,PRingCenter,R,GammaOriented, SmoothModel,SmoothParam,bComputeGrad,Uind,Grad)
    end subroutine

    !> Induction from a ring at a control point, within the referential of a ring: coordinates (r,theta, z), the ring is at z=0
    subroutine fUi_Ring11_loc(Rcp, Zcp, Gamma, R, Uind)
        use EllipticIntegrals, only: elliptic_ke
        use MathConstants, only: fourpi_inv
        ! Arguments declarations
        real(MK), intent(in) :: Rcp   !< Control point radius
        real(MK), intent(in) :: Zcp   !< Control point z location
        real(MK), intent(in) :: Gamma !< Ring intensity
        real(MK), intent(in) :: R     !< Ring radius
        real(MK), dimension(3), intent(inout) :: Uind  !< side effects, in coord (r, theta, z)
        ! Variable declarations
        real(MK), parameter :: MIN_DIST=1e-7 ! threshold for using axis formula
        real(MK) :: a  !<
        real(MK) :: AA !<
        real(MK) :: B  !<
        real(MK) :: E  !<
        real(MK) :: I1  !<
        real(MK) :: I2  !<
        real(MK) :: K  !<
        real(MK) :: m  !<
        !real(MK) :: aL,x,rr,cons,phi,a2  !<  Lewis notations

        if (Rcp<MIN_DIST*R) then
            !! Using Axis formula : v_z=-Gamma/(2R) *1 / (1+(z/R)^2)^(3/2) See Naca 1184 - Watsle, or Glauert Contraction of slipstream of an airscrew 1926
            Uind(1)=0
            Uind(3)=Gamma/(2_MK*R)*(1_MK/((1_MK+(Zcp/R)**2)**(1.5_MK))) ! no minus sign... it's a matter of convention
        else
            ! Formulation uses Formula from Yoon 2004.. it's not my favorite though..
            !Rcp=sqrt(Xcp.**2+Ycp.**2)
            !Zcp=Zcp-z ! !!!! In this case, rememeber to remove z to be in Ring coordinate system
            a=sqrt((Rcp+R)**2+(Zcp)**2)
            m=4_MK*Rcp*R/(a**2)
            AA=(Zcp)**2+Rcp**2+R**2
            B=-2_MK*Rcp*R
            call elliptic_ke(m, K, E)
            I1=4_MK/a*K
            I2=4_MK/a**3*E/(1_MK-m)
            Uind(1)=( (Zcp)/B )*(I1-AA*I2)      ! vr  !! Watchout GammaR/4pi removed , norm2(r_i) for yoon is R
            Uind(3)=((R+Rcp*AA/B)*I2-Rcp/B*I1 ) ! vz
            ! General scaling
            Uind(1:3)=Uind(1:3)*Gamma*fourpi_inv*R

            ! Lewis
!             x =Zcp/R;
!             rr=Rcp/R;
!             aL    = x**2+(rr-1._MK)**2
!             a2    = x**2+(rr+1._MK)**2
!             cons = 0.5_MK/(pi*R*sqrt(x**2+(rr+1._MK)**2))*Gamma
!             m    = 4._MK*rr/a2
!             phi=atan(sqrt(4._MK*rr/aL))
!             if(m>0.999) then
!                 K=log(4./cos(phi))
!                 E=1.+0.5*(K-1.0/1.2) * cos(phi)**2
!             else
!                 call elliptic_ke(m, K, E)
!             endif
!             Uind(3)= -cons*     (K-(1.+2.*(rr-1.)/aL)*E);
!             Uind(1)=  cons*x/rr*(K-(1.+2.*     rr/aL)*E);
            ! Oye with Grad
!                 zc=Zcp_prime(i)
!                 rc=Rcp_prime(i)
!                 a=sqrt((rc+R)**2+(zc)**2)
!                 m=4*rc*R/a**2
!                 A=(zc)**2+rc**2+R**2
!                 B=-2*rc*R
!                 call ellipke( m , K,E)
!                 I1=4/a*K
!                 I2=4/a**3*E/(1-m)
!                 uip(i,2)=0
!                 uip(i,1)=Gamma/(4*pi)*R*( (zc)/B )*(I1-A*I2) ! vr hopefully
!                 uip(i,3)=Gamma/(4*pi)*R*((R+rc*A/B)*I2-rc/B*I1 ) ! vz hopefully
!                 if (bComputeGrad) then
!                     a=(R + rc)**2 + zc**2
!                     b=(R + rc)**2 + zc**2
!                     ! X and Z are the same
!                     ! Y and R are the same
!                     Uzz=(1/2).*a.**(-3/2).*b.**(-2).*pi.**(-1).*zc.*(((-7).*R.**4+6.*R.**2.*rc.**2+rc.**4+2.*((-3).*R.**2+rc.**2).*zc.**2+zc.**4).*E+(-1).*b.*((-1).*R.**2+rc.**2+zc.**2).*K)
!
!                     Uzr=(1/2).*a.**(-3/2).*b.**(-2).*pi.**(-1).*rc.**(-1).*(((R.**2+(-1).*rc.**2).**2.*(R.**2+rc.**2)+2.*(R.**4+(-6).*R.**2.*rc.**2+rc.**4).*zc.**2+(R.**2+rc.**2).*zc.**4).*E+(-1).*b.*((R.**2+(-1).*rc.**2).**2+(R.**2+rc.**2).*zc.**2).*K)
!
!                     Urr= (-1/2).*a.**(-3/2).*b.**(-2).*pi.**(-1).*rc.**(-2).*zc.*((2.*rc.**6+5.* rc.**4.*(R.**2+zc.**2)+4.*rc.**2.*((-2).*R.**2+zc.**2).*(R.**2+zc.**2)+( R.**2+zc.**2).**3).*E+(-1).*b.*(2.*rc.**4+3.*rc.**2.*((-1).* R+zc).*(R+zc)+(R.**2+zc.**2).**2).*K)
!
!                     Urz=(1/2).*a.**(-3/2).*b.**(-2).*pi.**(-1).*rc.**(-1).*(((R.**2+(-1).* rc.**2).**2.*(R.**2+rc.**2)+2.*(R.**4+(-6).*R.**2.*rc.**2+rc.**4).*zc.**2+( R.**2+rc.**2).*zc.**4).*E+(-1).*b.*((R.**2+(-1).*rc.**2).**2+ (R.**2+rc.**2).*zc.**2).*K)
!
!                     grad(i,1)=Urr*Gamma
!                     grad(i,3)=Urz*Gamma
!                     grad(i,7)=Uzr*Gamma
!                     grad(i,9)=Uzz*Gamma
!                 end if

        endif
    end subroutine fUi_Ring11_loc

    subroutine fui_ring11_loc_c(Rcp, Zcp, Gamma, R, Uind) bind(C,name='fui_ring11_loc_c')
        use PrecisionMod, only: C_DOUBLE
        !DEC$ ATTRIBUTES DLLEXPORT :: fui_ring11_loc_c
        !GCC$ ATTRIBUTES DLLEXPORT :: fui_ring11_loc_c
        ! Arguments declarations
        real(C_DOUBLE), intent(in) :: Rcp   !< Control point radius
        real(C_DOUBLE), intent(in) :: Zcp   !< Control point z location
        real(C_DOUBLE), intent(in) :: Gamma !< Ring intensity
        real(C_DOUBLE), intent(in) :: R     !< Ring radius
        real(C_DOUBLE), dimension(3), intent(inout) :: Uind  !< side effects, in coord (r, theta, z)
        real(MK), dimension(3) :: Uind_MK

        Uind_MK=0.0_MK
        call fUi_Ring11_loc(Rcp, Zcp, Gamma, R, Uind_MK)
        Uind=real(Uind_MK,C_DOUBLE)

    end subroutine

end module

max_stars_repo_name: moonieann/welib

max_issues_repo_name: moonieann/welib

max_forks_repo_name: moonieann/welib


! { dg-do compile }
! PR 51502 - this was wrongly detected to be implicit pure.
module m
  integer :: i
contains
  subroutine foo(x)
    integer, intent(inout) :: x
    outer: block
      block
        i = 5
      end block
    end block outer
  end subroutine foo
end module m

! { dg-final { scan-module-absence "m" "IMPLICIT_PURE" } }


SUBROUTINE ReadSecondPart_of_in1(El,Elo,lapw,loor,ilo,nlo, Vr,ZZ,jri,mult,dx,r0,Efermi,aname,REL,fh_in1,fh_cout,fh_scf,nat,nrad,nloat,lmax,lomax,Qprint)
  IMPLICIT NONE
  REAL*8,  intent(out) :: El(1:lmax,nat), Elo(0:lomax,1:nloat,nat)
  LOGICAL, intent(out) :: lapw(0:lmax-1,nat), loor(0:lomax,nat)
  INTEGER, intent(out) :: ilo(0:lomax,nat)
  INTEGER, intent(out) :: nlo
  REAL*8,  intent(in)  :: Vr(nrad,nat), ZZ(nat)
  INTEGER, intent(in)  :: jri(nat), mult(nat)
  REAL*8,  intent(in)  :: dx(nat), r0(nat), Efermi
  CHARACTER*10, intent(in) :: aname(nat)
  LOGICAL, intent(in)  :: REL    ! relativistic Dirac or no
  INTEGER, intent(in)  :: fh_in1, fh_cout, fh_scf
  INTEGER, intent(in)  :: nat, nrad, nloat, lmax, lomax
  LOGICAL, intent(in)  :: Qprint
  ! locals
  INTEGER :: jatom, i, j, iapw, NL_read, L, L_previous
  REAL*8  :: RNET(nrad), Ei, dE
  CHARACTER*4 ::  EMAIN
  CHARACTER*67:: ERRMSG

  NLO = 0
  DO jatom = 1, NAT
     do i=1,jri(jatom)  ! Radial mesh created just once at the beginning
        RNET(i)=r0(jatom)*exp(dx(jatom)*(i-1.d0))
     enddo
     ! Reads the default setting from case.in1
     iapw=0
     READ (fh_in1,*,err=41) Ei, NL_read,iapw              ! EI==Global E par; NL_read==N cases; iapw=APW/LAPW
     if(abs(EI-0.3d0).lt.1.d-6) EI=Efermi-0.2d0  ! if EI is 0.3 than we set it to EF-0.2
     if(iapw.eq.1) then
        lapw(0:lmax-1,jatom)=.false.  ! all l's for this atom have apw+lo (if not changed below)
     else
        lapw(0:lmax-1,jatom)=.true.   ! all l's for this atom have lapw  (if not changed below)
     endif
41   continue

     if (Qprint) then
        if(iapw.eq.1) then
           WRITE (fh_scf,6000) ANAME(JATOM), jatom, Ei,'  APW'
        else
           WRITE (fh_scf,6000) ANAME(JATOM), jatom, Ei,' LAPW'
        endif
     endif
     ! initialization
     loor(0:lomax,jatom) = .false.
     ilo(0:lomax,jatom) = 0
     El(1:lmax,jatom) = EI
     if (nloat.gt.3) then ! LO type local orbitals possible
        ! for apw+lo (lapw+lo) we set to 99997.0 (99999.0)
        Elo(0:lomax,1:nloat,jatom) = 99997.0D+0    ! set to large number, so we know it is not present
        DO l=0,lomax
           if(lapw(l,jatom)) Elo(l,1:nloat,JATOM) = 99999.0D+0  ! large number, hence LO is not present
        ENDDO
     else                 ! LO type local orbitals not possible, just apw+lo might still be there
        ! for apw+lo (lapw+lo) we set to 997.0 (999.0)
        Elo(0:lomax,1:nloat,jatom) = 997.0D+0      ! set to large number, so we know it is not present
        DO l=0,lomax
           if(lapw(l,jatom)) Elo(l,1:nloat,JATOM) = 999.0D+0    ! large number, hence LO is not present
        ENDDO
     endif

     L = -1
     DO j=1,NL_read
        L_previous = L   ! previous L
        iapw=0
        READ (fh_in1,'(1X,I1,2F10.5,A4,i2)',err=42) L, EI, dE, EMAIN, iapw  !
        if (dE.lt.-1.0d-8) dE=-dE/10000.d0           ! if de<0 than de is really small positive, namely 1e4-times smaller than written
        if(abs(Ei-0.3d0).lt.1.d-6.and.L.ne.L_previous) Ei=Efermi - 0.2d0 ! if Ei=0.3 and not LO   => Ei=EF-0.2
        if(abs(Ei-0.3d0).lt.1.d-6.and.L.eq.L_previous) Ei=Efermi + 0.2d0 ! if Ei=0.3 and it is LO => Ei=EF+0.2
42      continue

        IF (dE .GT. 1.0D-5) THEN  ! finite de>1e-5 =>  we can search for the center of the band
           CALL SELECT(L,Ei,dE,EMAIN,REL,VR(1,JATOM),R0(jatom),dx(jatom),jri(jatom),ZZ(JATOM),jatom,Efermi)
        ELSE
           if (Qprint) then
              WRITE (fh_cout,6010)  L,jatom, L, Ei
              WRITE (fh_scf,6010) L,jatom, L, Ei
           endif
        ENDIF

        IF (L.NE.L_previous) THEN                ! the first n of certain l => must be apw or lapw
           El(l+1,jatom) = Ei                     ! Remember linearization energy for this apw/lapw
           if(iapw.eq.1) then                    ! this is apw+lo
              lapw(l,jatom)=.false.
              NLO = NLO + ((2*L+1))*MULT(JATOM)  ! we have apw+lo, constructued from u and dotu, hence many local orbitals in the ham.
              if (Qprint) then
                 write(fh_cout,*) '            APW+lo'
                 write(fh_scf,*) '            APW+lo'
              endif
              ilo(l,jatom)=1                     ! since this is the first for this l, we can set ilo to 1.
              Elo(l,1,jatom)=Ei                  ! Linearization energy for APW+lo is the same. We just date dotu.
              if(nloat.gt.3) then                ! real local orbitals are present
                 Elo(l,2,jatom)=99997.0D+0       ! next one is set to large number, so that we know when to stop
              else
                 Elo(l,2,jatom)=997.0D+0         ! next one is set to large number, hence it does not exist
              endif
           else                                  ! regular lapw
              if (Qprint) then
                 write(fh_cout,*) '            LAPW'
                 write(fh_scf,*) '            LAPW'
              endif
           endif
        ELSE ! not the first n of certain l => must be extra local orbital (not apw+lo but LO)
           ilo(l,jatom)=ilo(l,jatom)+1      ! yes, increase number of local orbitals
           if (ilo(l,jatom).gt.nloat) then  ! should not happen, as nloat should be larger than any number of LO's
              goto 902
           endif
           Elo(L,ilo(l,jatom),jatom) = Ei    ! remember this linearization energy for local orbital
           NLO = NLO + ((2*L+1))*MULT(JATOM) ! many local orbitals need to be added, namely, one for each equivalent atom and for each m.
           if (Qprint) then
              write(fh_cout,*)  '            LOCAL ORBITAL'
              write(fh_scf,*) '            LOCAL ORBITAL'
           endif
           loor(l,jatom)=.true.
        ENDIF
     ENDDO
  END DO

  !WRITE(6,*) 'loor,ilo='
  !DO jatom=1,NAT
  !   do l=0,2
  !      WRITE(6,*) l,jatom,loor(l,jatom), ilo(l,jatom)
  !   enddo
  !ENDDO

  return
902 CONTINUE
  WRITE (ERRMSG,9002) 'NLOAT too small', NLOAT,l,jatom
  CALL OUTERR('ATPAR',ERRMSG)
  STOP 'readSecondPart_of_in1 - Error'
  !
6000 FORMAT (/,10X,'ATOMIC SPHERE DEPENDENT PARAMETERS FOR ATOM  ',A10,/,':e__',i4.4,':',1X,'OVERALL ENERGY PARAMETER IS',F10.4/,10X,'OVERALL BASIS SET ON ATOM IS',A5)
6010 FORMAT (':E',i1,'_',i4.4,':',1X,'E(',I2,2H)=,F10.4)
9002 FORMAT (A,' (NLOAT =',I3,', L =',I3,' JATOM =',i3,')')
END SUBROUTINE ReadSecondPart_of_in1


text:
C    %W% %G%
C****************************************************************
C
C     File: clnup_ge.f
C
C     Purpose: Routine to complete PTI data processing from raw data
C              file.
C
c     Return code:  n = 0 : Success
c                   n = 1 : Error
c
C     Author: Walt Powell  Date: 21 May 1996
C     Called by: load_pti.f
C
C****************************************************************
      subroutine clnup_ge (error)
      integer error

      include 'ipfinc/parametr.inc'

      include 'ipfinc/prt.inc'
      include 'ipfinc/pti_data.inc'
      include 'ipfinc/blank.inc'
      include 'ipfinc/arcntl.inc'
      include 'ipfinc/bus.inc'
      include 'ipfinc/cbus.inc'
      include 'ipfinc/branch.inc'
      include 'ipfinc/area.inc'
      include 'ipfinc/alpha.inc'
      include 'ipfinc/qsdup.inc'
      include 'ipfinc/qksrt.inc'
      include 'ipfinc/alt_case.inc'
      include 'ipfinc/owner_cm.inc'
      include 'ipfinc/tx_misc.inc'

      common /is_batch / is_batch

      integer MAXWYEDELTA
      parameter (MAXWYEDELTA = 500)
      common /wye_delta/ num_wye_delta, lwye_delta(3,MAXWYEDELTA),
     &                   wye_delta(MAXWYEDELTA)

      integer MAXPTIRECORDS
      parameter (MAXPTIRECORDS = 16000)
      common /scratch/ count, array(4,MAXBUS), htable_2(MAXBUS),
     &                 nextptr_2(MAXBUS), count_newbus,
     &                 newbusno(MAXBUS), count_newzone,
     &                 newzoneno(MAXCZN), count_newown,
     &                 newownno(MAXOWN), tempc(MAXPTIRECORDS),
     &                 sort_tempc(MAXPTIRECORDS), txflag(MAXBUS),
     &                 txindex(2,MAXWYEDELTA), procnode(MAXBUS),
     &                 vangle(MAXBUS), itemp(MAXBUS)
      integer array, count, htable_2, count_newbus, newbusno,
     &        count_newzone, newzoneno, count_newown, newownno,
     &        sort_tempc, txflag, txindex, procnode
      character tempc*80

      character zn*2, xbuf*120, bus_name*8, id*2

      integer HASHSIZE
      parameter (HASHSIZE = MAXCZN - 1)

      external kmpzone, swpzon, kmparc, swparc, kmpari,
     &         swpari, kmp_ptiz, swp_ptiz, kmp_wydt, swp_wydt
      integer find_zon, find_ara, fnd_ptiz, fnd_ptib, fnd_ptic,
     &        status, rename_bus, rebldzon, h, p, pold,
     &        fnd_ptio, ptr, ftn_atoi, ptr_section(12), ptr9, ptr10,
     &        areanos(MAXCAR), find_hash, fnd_ptian, add_ptiz,
     &        fnd_ptiy
      logical finished
      character tempc1*10, tempc2*10, cbown*3,
     &          lastid*1, oldrcd(2)*120, newrcd*120

      error = 0
c
c     Rename any duplicate names
c
      do i = 1, num_znam
        num1 = fnd_ptiz (pti_znum(i))
        num2 = fnd_ptiy (pti_znam(i))
        if (num1 .ne. i .or. num2 .ne. i) then
          status = new_znnm (tempc1(1:2))
          if (status .eq. 0) then
c
c           Hash new name into old location
c
            h = 0
            do j = 1, 2
              h = h + h + ichar (tempc1(j:j))
            end do
            h = mod (h, MAXCZN-1) + 1
            p = htable_y(h)
            do while (p .gt. 0)
              if (tempc1(1:2) .ne. pti_znam(p)) then
                p = nextptr_y(p)
              else
                p = -p
              endif
            enddo
            if (p .eq. 0) then
              p = i
              pti_znam(p) = tempc1(1:2)
              nextptr_y(p) = htable_y(h)
              htable_y(h) = p
            endif
            write (errbuf(1), 10000) pti_znum(i), pti_znam(i),
     &        pti_znum(i), tempc1, pti_znum(p),
     &        pti_znam(p)
10000       format(' Duplicate zone names: ', i3, 1x, a, ' and ',
     &         i3, 1x, a, ' First pair renamed to ', i3, 1x, a)
            call prterx ('W', 1)
          endif
        endif
      enddo
c
c     Sort pti_anum() to generate missing zone, owner, and bus numbers
c
      do i = 1, MAXCAR
        areanos(i) = 0
        array(1,i) = i
        array(2,i) = 0
      enddo
c
c       key = sorts on the following fields:
c             1. bus(zone)
c             2. population
c             3. area
c             4. pti_znum
c             5. pti_anum
c             6. pti_znum, pti_anum
c
      key = 5
      call qiksrt (1, num_anam, kmp_ptiz, swp_ptiz)
      do i = 1, num_anam
        areanos(i) = pti_anum(array(1,i))
      enddo
      areanos(num_anam+1) = areanos(num_anam) + 10
c
c     Test zone-area uniqueness
c
      do i = 1, MAXBUS
         htable_2(i) = 0
         array(1,i) = 0
         array(2,i) = 0
         array(3,i) = 0
         array(4,i) = 0
         nextptr_2(i) = 0
      enddo
      count = 0
      do i = 1, num_hashn
         p = find_hash (pti_zone(i), pti_area(i))
         p = iabs (p)
         if (p .gt. 0) array(3,p) = array(3,p) + 1
      enddo

      dupsw = .false.
      key = 6                ! Set key to sort as follows:
c                            ! 1. Array(1,*) - PTI zone number
c                            ! 2. Array(3,*) - Zone-Area population
c                            ! 3. Array(2,*) - PTI area number

      call qiksrt (1, count, kmp_ptiz, swp_ptiz)
c
c     Rebuild zone-area hash tables
c
      do i = 1, MAXBUS
         htable_2(i) = 0
         nextptr_2(i) = 0
      enddo
      do i = 1, count
        h = mod (array(1,i), HASHSIZE) + 1
        p = htable_2(h)
        do while (p .gt. 0)         !search for existing entities
          if (array(1,i) .ne. array(1,p) .or.
     &        array(2,i) .ne. array(2,p)) then
            p = nextptr_2(p)
          else
            p = -p
          endif

        enddo
        if (p .eq. 0) then
          nextptr_2(i) = htable_2(h)
          htable_2(h) = i
        endif
      enddo
c
c     Find bus zone numbers missing in the zone data, create a
c     corresponding and unique name, and insert the number-name pair
c     into the number-name zone-hash tables.
c
      do i = 1, count
        num = fnd_ptiz (array(1,i)) ! Test if zone number exists
        if (num .eq. 0) then
          status = new_znnm (tempc1(1:2))
          if (status .eq. 0) then
            num = add_ptiz (array(1,i), tempc1(1:2))
            write (errbuf(1), 10010) array(1,i), tempc1(1:2)
10010       format(' Missing name for zone number ', i4,
     &        ' is given name ', a)
            call prterx ('W', 1)
          endif
        endif
      enddo

c     Flag incompatible zone-area designations
c
      last = 1
      next = last + 1
      ndup = 0
      do while (next .le. count + 1)
        if (next .le. count .and. array(1,next) .eq. array(1,last))
     &    then
          ndup = ndup + 1
        else
          if (ndup .gt. 0) then
c
c           First pass - flag duplicates zones (areas are unique)
c
            do i = 1, ndup+1
              j = next - (ndup + 2 - i)
              nptiz1 = fnd_ptiz (array(1,j))
              nptia1 = fnd_ptian (array(2,j))
              write (errbuf(1), 10020) array(1,j), pti_znam(nptiz1),
     &          array(2,j), pti_anam(nptia1), array(3,j)
10020         format(' Duplicate zone ', i4, 1x, a, ' area ', i4,
     &          1x, a, ' Population ', i4)
              call prterx ('W', 1)
            enddo
c
c           Second pass - rename duplicate zone of first, second, ...,
c           but not the last which is most populous.
c
            do i = 1, ndup
              j = next - (ndup + 2 - i)
              nptiz1 = fnd_ptiz (array(1,j))
              nptia1 = fnd_ptian (array(2,j))
c
c             Create a unique zone number higher than any existing
c             zone number
c
              numzone = 0
              do jx = 1, num_znam
                numzone = max0 (numzone, pti_znum(jx))
              enddo
              numzone = numzone + 1
c
c             Create a unique zone name
c
              status = new_znnm (tempc1(1:2))
c
c             Hash the number-name pair into the zone hash tables
c
              if (status .eq. 0) then
                num = add_ptiz (numzone, tempc1(1:2))
                if (num .gt. 0) then
                  write (errbuf(1), 10030) array(1,j),
     &              pti_znam(nptiz1), array(2,j), pti_anam(nptia1),
     &              pti_znum(num), pti_znam(num)
10030             format(' Zone ', i4, 1x, a, ' in area ', i4,
     &              1x, a, ' renamed to ', i4, 1x, a)
                  call prterx ('W', 1)
                  array(4,j) = num
                endif
              endif
            enddo
            ndup = 0
          endif
          last = next
        endif
        next = next + 1
      enddo

      do nb = 1, ntot
        if (jarzn(nb) .eq. 0) then
          tempc1 = ' '
          numarea = 0
        else
          nptia = fnd_ptic (jarzn(nb))
          jarzn(nb) = nptia
          if (nptia .gt. 0) then
            tempc1 = pti_anam(nptia)
            numarea = pti_anum(nptia)
          else
            tempc1 = ' '
            numarea = 0
          endif
        endif
        nptib = fnd_ptib (bus(nb), base(nb), tempc1)
        if (nptib .gt. 0) then
          numzone = pti_zone(nptib)
          alf2inp(nb) = numzone             ! Temporarily hold zone num
          nptiz = fnd_ptiz (numzone)
          if (nptiz .gt. 0) then
            p = find_hash (pti_znum(nptiz), numarea)
            p = iabs (p)
            if (array(4,p) .gt. 0) then
              write (errbuf(1), 10040) bus(nb), base(nb),
     &          pti_anam(nptia), pti_znam(nptiz), pti_znam(array(4,p))
10040         format(' Bus ', a, f6.1, ' in area ', a,
     &          ' has duplicate zone renamed from ', a, ' to ', a)
              call prterx ('W', 1)
              pti_zone(nptib) = pti_znum(array(4,p))
              zone(nb) = pti_znam(array(4,p))
            else
              zone(nb) = pti_znam(nptiz)
            endif
          else
            write (zone(nb), '(i2.2)') numarea
          endif

          nptio = fnd_ptio (pti_owner(nptib))
          if (nptio .gt. 0) then
            write (owner(nb), fmt='(i3.3)') pti_onum(nptio)
          else
            owner(nb) = '0'
          endif
        else
          alf2inp(nb) = 0
          jarzn(nb) = 0
          write (errbuf(1), 10050) bus(nb), base(nb)
10050     format(' Bus ', a, f6.1, ' has no assigned bus number ')
          call prterx ('W', 1)
        endif
      enddo

      status = rebldzon (error)
c
c     Sort zones (ACZNAM array)
c
      dupsw = .false.
      call qiksrt(1, nztot, kmpzone, swpzon)
      if (dupsw) then

c       Flag and remove duplicate records

        j = 1
        k = 2
        do while (k .le. nztot)
          if ( kmpzone(j,k) .lt. 0 ) then
            j = j + 1
            if (j .lt. k) call swpzon (j,k)
            k = k + 1
          else
            write (errbuf(1), 10060) acznam(k)
10060       format('0 Duplicate zones (', a, '). Case Aborted!')
            if (is_batch .eq. 0) then
              call prterx ('E',1)
            else
              call prterx ('F',1)
            endif
            k = k + 1
            error = 1
          endif
        enddo
      endif

c     Sort area interchange (ARCNAM array)

      dupsw = .false.
      call qiksrt (1, ntotc, kmparc, swparc)
      if (dupsw) then

c       Flag and remove duplicate records

        j = 1
        k = 2
        do while (k .le. ntotc)
          if ( kmparc(j,k) .lt. 0 ) then
            j = j + 1
            if (j .lt. k) call swparc (j,k)
            k = k + 1
          else
            write (errbuf(1), 10070)
10070       format('0 Duplicate area interchange records..',
     &             ' Case Aborted!')
            call bcdarc (j,xbuf)
            write (errbuf(2), 10080) xbuf(1:80)
10080       format(1x, '(', a, ')')
            call bcdarc (k,xbuf)
            write (errbuf(3), 10080) xbuf(1:80)
            if (is_batch .eq. 0) then
              call prterx ('E',3)
            else
              call prterx ('F',3)
            endif
            k = k + 1
            error = 1
          endif
        enddo
      endif

c     Build array(*,*) of zone-area cross-reference

      do nb = 1, ntot
        if (jarzn(nb) .gt. 0) then
          tempc1 = pti_anam(jarzn(nb))
        else
          tempc1 = ' '
        endif
        nptib = fnd_ptib (bus(nb), base(nb), tempc1)
        if (nptib .gt. 0) then
          array(1,nb) = find_zon (zone(nb))
        else
          array(1,nb) = 19999
          zone(nb) = ' '
        endif
        if (jarzn(nb) .gt. 0) then
          array(2,nb) = find_ara (pti_anam(jarzn(nb))(1:8))
        else
          array(2,nb) = 19999
        endif
      enddo

      count = ntot
      dupsw = .false.
      key = 2                ! Set key to sort as follows:
c                            ! 1. zone
c                            ! 2. area

      call qiksrt (1, count, kmp_ptiz, swp_ptiz)

      finished = .false.
      do while (count .gt. 0 .and. .not. finished)
        if (array(1,count) .ne. 19999) then
          finished = .true.
        else
          count = count - 1
        endif
      enddo

      if (dupsw) then
        j = 1
        k = 2
        do while (k .le. count)
          if (array(1,j) .eq. array(1,k) .and.
     &        array(2,j) .eq. array(2,k)) then
            k = k + 1
          else
            if (array(1,j) .eq. array(1,k) .and.
     &          array(2,j) .lt. array(2,k)) then
              if (array(2,j) .lt. 19999) then
                tempc1 = arcnam(array(2,j))
              else
                tempc1 = ' '
              endif
              if (array(2,k) .lt. 19999) then
                tempc2 = arcnam(array(2,k))
              else
                tempc2 = ' '
              endif
              write (errbuf(1), 10090) acznam(array(1,j)),
     &          tempc1, tempc2
10090         format(' Zone ', a, ' belongs to areas (', a,
     &           ') and (', a, ')')
              call prterx ('W', 1)
              error = 1
            endif
            j = j + 1
            if (j .lt. k) then
              array(1,j) = array(1,k)
              array(2,j) = array(2,k)
            endif
            k = k + 1
          endif
        enddo
        count = j
      endif

      do i = 1, nztot
        acznum(i) = 0
      enddo

      lastzone = 0
      lastarea = 0
      do i = 1, count
        izone = array(1,i)
        iarea = array(2,i)
        if (izone .eq. lastzone) then
          if (lastarea .lt. 19999) then
            tempc1 = arcnam(lastarea)
          else
            tempc1 = ' '
          endif
          if (iarea .lt. 19999) then
            tempc2 = arcnam(iarea)
          else
            tempc2 = ' '
          endif
          write (errbuf(1), 10100) acznam(izone), tempc1, tempc2
10100     format(' Zone ', a, ' belongs to areas (', a,
     &           ') and (', a, ')')
          call prterx ('W', 1)
          error = 1
        else
          acznum(izone) = iarea
        endif
        lastzone = izone
        lastarea = iarea
      enddo

      if (count .gt. nztot) then
        write (errbuf(1), 10110) count, nztot
10110   format(' Inconsistent PTI Zone-Area matching: ',i4,
     &          ' zone-areas encountered, ', i4,
     &          ' zones-area specified.')
        call prterx ('W',1)
        error = 1
      endif

      do i = 1, nztot
        iarea = acznum(i)
        if (iarea .gt. 0 .and. iarea .lt. 19999) then
          j = 1
          finished = .false.
          do while (j .le. MAXCZN .and. .not. finished)
            zn = arczns(j,iarea)
            if (zn .eq. '  ') then
              arczns(j,iarea) = acznam(i)
              finished = .true.
            else
              j = j + 1
            endif
          enddo
          if (.not. finished) then
            write (errbuf(1), 10120) MAXCZN, acznam(i), arcnam(iarea)
10120       format(' More than ', i4, ' zones in assigning zone (', a,
     &           ') to area (', a, ')')
            call prterx ('W', 1)
            error = 1
          endif
        endif
      enddo

      call space(1)
      write (outbuf,10130)
10130 format(t10,'AREA 1',t22,'AREA 2',t34,'SCHEDULED',t46,
     &           '   NET   ',t58,'  CONTROL BUS',t78,'ZONES(S)')
      call prtout(1)

      write (outbuf, 10140)
10140 format (t35,'EXPORT ',t46,' EXPORT ',t58,'NAME       BASE')
      call prtout(1)

      write (outbuf, 10150)
10150 format(t34,'  (MW)  ',t46,'  (MW)  ',t78,10('-- '))
      call prtout(1)

      do i=1,ntotc
         call space(1)
         call list_are (i)
      enddo

      outbuf= '0'
      call prtout(1)
c
c     Duplicate PTI owner data into /owner_cm/
c
      do i = 1, num_onam
	owner_number(i) = pti_onum(i)
        owner_code(i) = pti_onam(i)
      enddo
c
c     1. Delete buses without branches
c     2. Sort branches
c     3. Convert bus owners from ASCII numbers to names
c     4. Remove section numbers from non-sections
C     5. Add pi-equivalent for valid sections
C     6. Consolidate excessive sections (> 10)
c     7. Change brid() for "R" records to blank
c
      do nb = 1, ntot
        if (kbsdta(16,nb) .eq. 0 .and. bus(nb) .ne. srtlst) then
c
c         Remove isolated bus from bus hash table by renaming
c
          write (errbuf(1), 10160) bus(nb), base(nb), zone(nb)
10160     format (' Bus (', a8, f7.1, 1x, a2,
     &      ') has no branches and is deleted')
          call prterx ('W',1)

          bus_name = srtlst
          bus_base = 9999.0
          status = rename_bus (nb, bus_name, bus_base)
c
c         Flag nb as deleted
c
          do i = 1, 11
             kbsdta(i,nb) = 0
          enddo
          capcor(1,nb) = 0.0
          capcor(2,nb) = -9.0e10
        endif
      enddo

      call sortbus ()
      call srtbrnch ()

      do nb = 1, ntot
        if (kbsdta(16,nb) .eq. 0 .and. bus(nb) .ne. srtlst) then
        else
c
c         Convert bus owners from ASCII numbers to names
c
          numowner = ftn_atoi (owner(nb))
          nptio = fnd_ptio (numowner)
          if (nptio .gt. 0) then
            owner(nb) = pti_onam(nptio)
          else
            owner(nb) = ' '
          endif
c
c         Convert +bus owners from ASCII numbers to names
c
          ptr = kbsdta(15,nb)
          do while (ptr .gt. 0)
            call getchr (3, cbown, kbctbl(10,ptr))
            numowner = ftn_atoi (cbown)
            nptio = fnd_ptio (numowner)
            if (nptio .gt. 0) then
              cbown = pti_onam(nptio)
            else
              cbown = ' '
            endif
            call putchr (3, cbown, kbctbl(10,ptr))
            ptr = bctbl_nxt(ptr)
          enddo
c
c         3. Convert branch owners from ASCII numbers to names
c         4. Remove section numbers from non-sections
C         5. Add pi-equivalent for valid sections
c
          lastk2 = 0
          lastid = ' '
          numsec = 0
          ptr = kbsdta(16,nb)
          pold = 0
          do while (ptr .gt. 0)
            if (brtype(ptr) .eq. 4) then
              brid(ptr) = ' '
            else
              nbr = iabs (brnch_ptr(ptr))
              call getchr (3, cbown, kbrnch(3,nbr))
              numowner = ftn_atoi (cbown)
              nptio = fnd_ptio (numowner)
              if (nptio .gt. 0) then
                cbown = pti_onam(nptio)
              else
                cbown = ' '
              endif
              call putchr (3, cbown, kbrnch(3,nbr))
c
c             If non-zero section encountered, look ahead for additional
c             sections
c
              if (ky(ptr) .ne. lastk2 .or. brid(ptr) .ne. lastid) then
                numsec = 1
                ptr_section(numsec) = ptr
                p = brnch_nxt(ptr)
                do while (p .gt. 0 .and. ky(p) .eq. ky(ptr) .and.
     &                    brid(p) .eq. brid(ptr))
                  numsec = numsec + 1
                  ptr_section(numsec) = p
                  p = brnch_nxt(p)
                enddo
                if (numsec .gt. 1) then
                  ltot2 = ltot2 + 1
                  if (nb .lt. ky(ptr)) then
                    if (ltot+1 .ge. MAXBRN) then
                      write (errbuf(1), 10170) MAXBRN
10170                 format ('More than ',i5,
     &                 ' branch records. Overflow occurred adding pi-equ
     &ivalent section for branch:')
                      write (errbuf(2), 10180) bus(kx(ptr)),
     &                  base(kx(ptr)), bus(ky(ptr)), base(ky(ptr)),
     &                  brid(ptr), brsect(ptr)
10180                 format(11x, '(', a, f6.1, 1x, a8, f6.1, 1x, a, i2,
     &                  ')')
                      call prterx ('E',2)
                      ltot = 1
                      error = 1
                      go to 900
                    endif
                    ltot = ltot + 1
                    kbrnch(1,ltot) = 1
                    do i = 2, 18
                      kbrnch(i,ltot) = 0
                    enddo
                    nbr = ltot
                  else
                    p = numbrn (ky(ptr), kx(ptr), brid(ptr), 0)
                    if (p .le. 0) then
                      write (errbuf(1), 10190) bus(ky(ptr)),
     &                  base(ky(ptr)), bus(kx(ptr)), base(kx(ptr)),
     &                  brid(ptr), brsect(ptr)
10190                 format ('Error renaming section (', a8, f6.1, 1x,
     &                  a8, f6.1, 1x, a1, i2,
     &                  ') - could not locate branch')
                      call prterx ('E', 1)
                      error = 1
                      go to 900
                    else
                      nbr = -brnch_ptr(p)
                    endif
                  endif

c                 Link up branch.

                  if (pold .eq. 0) then
                    kbsdta(16,nb) = ltot2
                  else
                    brnch_nxt(pold) = ltot2
                  endif
                  brnch_nxt(ltot2) = ptr
                  brnch_ptr(ltot2) = nbr

                  kx(ltot2) = nb
                  ky(ltot2) = ky(ptr)
                  brid(ltot2) = brid(ptr)
                  brsect(ltot2) = 0
                  nbr = iabs (brnch_ptr(ltot2))
                  brtype(ltot2) = kbrnch(1,nbr)
c
c                 if more than 9 sections, consolidate 9 and 10
c
                  if (numsec .gt. 9) then
                    if (brnch_ptr(ltot2) .gt. 0) then
                      ptr9 = ptr_section(9)
                      ptr10 = ptr_section(10)
                      call bcdbrn (ptr9, oldrcd(1), ierror)
                      call bcdbrn (ptr10, oldrcd(2), ierror)
                      numtemp = 2
                      status = gtpieqiv (numtemp, oldrcd, newrcd)
                      if (status .eq. 0) then
                        write (outbuf, 10200)
10200                   format (' Old sections 9 and 10:')
                        call prtout (1)
                        write (outbuf, 10210) oldrcd(1)(1:80)
10210                   format (' (', a, ')')
                        call prtout (1)
                        write (outbuf, 10210) oldrcd(2)(1:80)
                        call prtout (1)
                        write (outbuf, 10220)
10220                   format (' New equivalent section 9:')
                        call prtout (1)
                        write (outbuf, 10210) newrcd(1:80)
                        call prtout (1)
                        nbr = iabs(brnch_ptr(ptr9))
                        read (newrcd, 10230, err=900) (brnch(k,nbr),
     &                    k=5,10)
10230                   format (bz, 38x, 6f6.5)
                        kbrnch(1,nbr) = 8
                        brtype(ptr9) = kbrnch(1,nbr)
                      endif
c
c                     De-link ptr10
c
                      brnch_nxt(ptr9) = brnch_nxt(ptr10)
                      write (errbuf(1), 10240) bus(kx(ptr10)),
     &                  base(kx(ptr10)), bus(ky(ptr10)),
     &                  base(ky(ptr10)), brid(ptr10), brsect(ptr10)
10240                 format(' Consolidating section (', a, f6.1,
     &                  1x, a8, f6.1, 1x, a, 1x, i2, ')')
                      call prterx ('W', 1)
                    else
c
c                     De-link ptr10
c
                      brnch_nxt(ltot2) = ptr_section(2)
                      ptr10 = ptr_section(1)
                      write (errbuf(1), 10240) bus(kx(ptr10)),
     &                  base(kx(ptr10)), bus(ky(ptr10)),
     &                  base(ky(ptr10)), brid(ptr10), brsect(ptr10)
                      call prterx ('W', 1)
                    endif
                  endif
                else
                  brsect(ptr) = 0
                endif
                lastk2 = ky(ptr)
                lastid = brid(ptr)

              endif
            endif
            pold = ptr
            ptr = brnch_nxt(ptr)
          enddo

        endif

      enddo
c
c     Check consistency of neutral buses in 3-terminal Tx's.
c
      do i = 1, num_3term
        k1 = tx_3term(1,i)                 ! from bus
        k2 = tx_3term(2,i)                 ! to bus
        k3 = tx_3term(3,i)                 ! reg bus
        k4 = tx_3term(4,i)                 ! 3-winding pt bus
        k5 = tx_3term(5,i)                 ! tertiary bus
        id = char (int (tx_3term(6,i)))
        ptr = kbsdta(16,k4)
        do while (ptr .gt. 0)
          if (ky(ptr) .eq. k1 .or. ky(ptr) .eq. k2 .or. ky(ptr) .eq. k5)
     &       then
          else
            write (errbuf(1), 10250) bus(k4), base(k4), bus(k1),
     &         base(k1), bus(k2), base(k2), bus(k5), base(k5)
10250       format (' Neutral bus ', a8, f6.1, ' in 3-winding Tx ',
     &              3(a8, f6.1, 1x))
            write (errbuf(2), 10260) bus(ky(ptr)), base(ky(ptr))
10260       format (' has an illegal branch to ', a8, f6.1)
            call prterx ('W', 2)
            error = 1
          endif
          ptr = brnch_nxt(ptr)
        enddo
      enddo
c
c     Correct wye-delta voltage transformation
c
      if (num_wye_delta .gt. 0) then
c
c       Create double entry entities, sort, and remove duplicates
c
        do i = 1, num_wye_delta
          lwye_delta(1,i+num_wye_delta) = lwye_delta(2,i)
          lwye_delta(2,i+num_wye_delta) = lwye_delta(1,i)
          wye_delta(i+num_wye_delta) = -wye_delta(i)
        enddo
        num_wye_delta = 2 * num_wye_delta
        call qiksrt (1, num_wye_delta, kmp_wydt, swp_wydt)
        i = 1
        do while (i .lt. num_wye_delta)
          j = i + 1
          do while (j .le. num_wye_delta)
            if (lwye_delta(1,i) .eq. lwye_delta(1,j) .and.
     &          lwye_delta(2,i) .eq. lwye_delta(2,j)) then
              do k = j, num_wye_delta - 1
                lwye_delta(1,k) = lwye_delta(1,k+1)
                lwye_delta(2,k) = lwye_delta(2,k+1)
                wye_delta(k) = wye_delta(k+1)
              enddo
              num_wye_delta = num_wye_delta - 1
            else
              j = num_wye_delta + 1
            endif
          enddo
          i = i + 1
        enddo
c
c       Initialize processing arrays
c
        do i = 1, ntot
          txflag(i) = 0
          procnode(i) = 0
          vangle(i) = 0.0
        enddo
        do i = 1, num_wye_delta
          txindex(1,i) = 0
          txindex(2,i) = 0
        enddo
        num_txflag = 0
        do i = 1, num_wye_delta
          num_txflag = num_txflag + 1
          k = lwye_delta(1,i)
          m = lwye_delta(2,i)
          next = txflag(k)
          last = 0
          do while (next .gt. 0)
            last = next
            next = txindex(2,next)
          enddo
          if (last .eq. 0) then
            txflag(k) = num_txflag
          else
            txindex(2,last) = num_txflag
          endif
          txindex(1,num_txflag) = i
        enddo
c
c       Propagate from highest base kV in lwye_delta array
c
        do i = 1, num_wye_delta
          k = lwye_delta(1,i)
          m = lwye_delta(2,i)
          if (procnode(k) .eq. 0) then
c
c           Flag all primary neighors (node "k") with same angle
c
            procnode(k) = i
            next = 1           ! "next" is next index to process
            last = 1           ! "last" is last available index in stack
            itemp(last) = k
            do while (next .le. last)
              nb = itemp(next) ! "nb" is kernel node
              ptr = kbsdta(16,nb)
              do while (ptr .gt. 0)
                mb = ky(ptr)
                if (procnode(mb) .eq. 0) then
c
c                 Test whether branch (nb,mb) is in wye_delta
c
                  jx = txflag(mb)
                  finished = (jx .eq. 0)
                  do while (.not. finished)
                    jy = txindex(1,jx)
                    if (lwye_delta(1,jy) .eq. mb .and.
     &                  lwye_delta(2,jy) .eq. nb) then
                      procnode(mb) = jy
                      vangle(mb) = vangle(nb) + wye_delta(jy)
                      finished = .true.
                      last = last + 1
                      itemp(last) = mb
                    else
                      jx = txindex(2,jx)
                      finished = (jx .eq. 0)
                    endif
                  enddo
                  if (procnode(mb) .eq. 0) then
                    procnode(mb) = procnode(nb)
                    vangle(mb) = vangle(nb)
                    last = last + 1
                    itemp(last) = mb
                  endif
                endif
                ptr = brnch_nxt(ptr)
              enddo
              next = next + 1
            enddo
          endif
          if (procnode(m) .eq. 0) then
            write (errbuf(1), 10270) bus(k), base(k), bus(m),
     &        base(m)
10270       format (' Wye-delta TX ', a8, f6.1, 1x, a8, f6.1,
     &        ' improperly processed')
            call prterx ('W', 1)
            error = 1
          endif
        enddo
c
c       Flag non-processed buses
c
        write (outbuf, 10280)
10280	format (' The following isolated buses could not have their vol
     &tages compensated')
        call prtout (1)
        do i = 1, ntot
          if (procnode(i) .ne. 0) then
            vmag = dsqrt (e(i) ** 2 + f(i) ** 2)
            angle = datan2 (f(i), e(i)) - 0.0174532 * vangle(i)
            e(i) = vmag * cos (angle)
            f(i) = vmag * sin (angle)
          else if (kbsdta(1,i) .gt. 0) then
            write (outbuf, 10290) bus(i), base(i)
10290       format (1x, a8, f7.1)
            call prtout (1)
          endif
        enddo
c
c       Test completedness of processing wye_delta
c
        do i = 1, num_wye_delta
          k = lwye_delta(1,i)
          m = lwye_delta(2,i)
          if (procnode(k) .eq. 0) then
            write (errbuf(1), 10300) 'Primary', bus(k), base(k)
10300       format (a, ' node of wye-delta Tx not processed ',
     &        a8, f7.1)
            call prterx ('W', 1)
          endif
          if (procnode(m) .eq. 0) then
            write (errbuf(1), 10300) 'Secondary', bus(m), base(m)
            call prterx ('W', 1)
          endif
          if (procnode(k) .ne. 0 .and. procnode(m) .ne. 0 .and.
     &        abs(vangle(k)-vangle(m)-wye_delta(i)) .gt. 0.01) then
            write (errbuf(1), 10310) bus(k), base(k), bus(m), base(m),
     &        vangle(k), vangle(m), wye_delta(i)
10310       format ('Angle error in wye-delta Tx ', a8, f7.1, 1x,
     &         a8, f7.1, ' Nodal bias', 2f7.1, ' Wye-delta ', f7.1)
            call prterx ('W', 1)
          endif
        enddo
      endif

  900 continue
      return
      end

ext: f

max_stars_repo_name: mbheinen/bpa-ipf-tsp

max_issues_repo_name: cuihantao/bpa-ipf-tsp

max_forks_repo_name: mbheinen/bpa-ipf-tsp


text:
!     .................................................
!             ____  _       _   ____  _____   _
!            |  _ \| |     |_| |  _ \|  ___| |_|
!            | |_) | |___   _  | |_) | |___   _
!            |  _ /|  _  | | | |  _ /|___  | | |
!            | |   | | | | | | | |    ___| | | |
!            |_|   |_| |_| |_| |_|   |_____| |_|
!     .................................................
!     PhiPsi:     a general-purpose computational
!                 mechanics program written in Fortran.
!     Website:    http://phipsi.top
!     Author:     Shi Fang from Huaiyin Institute of
!                 Technology, HuaiAn, JiangSu, China
!     Contact me: shifang@hyit.edu.cn
!     ------------------------------------------------
!     Please cite the following papers:
!     (1)Shi F, Wang X L, Liu C, Liu H, Wu H A. An
!        XFEM-based method with reduction technique
!        for modeling hydraulic fracture propagation
!        in formations containing frictional natural
!        fractures. Engineering Fracture Mechanics,
!        2017, 173: 64-90.
!     (2)Shi F, Wang X L, Liu C, Liu H, Wu H A. A
!        coupled extended finite element approach
!        for modeling hydraulic fracturing in
!        consideration of proppant. Journal of
!        Natural Gas Science and Engineering, 2016,
!        33: 885-897.
!     (3)Shi F, Wang X L, Liu C, Liu H, Wu H A. An
!        XFEM-based numerical model to calculate
!        conductivity of propped fracture considering
!        proppant transport, embedment and crushing.
!        Journal of Petroleum Science and Engineering,
!        2018, 167: 615-626..

      SUBROUTINE Boundary_Cond(Total_FD,isub,freeDOF,num_freeDOF)
c     Deal with the boudary condition.

      use Global_Common
      use Global_Model
      use Global_Elem_Area
      use Global_Crack

      implicit none
      integer,intent(in)::isub
      integer,intent(in)::Total_FD
      integer,intent(out)::freeDOF(Total_FD),num_freeDOF
      integer i,cur_Node,i_fd
      integer c_index,fixedDOF(Total_FD)

      print *,'    Dealing with the boundary condition......'

      !Initialize vector of fixed DOFs.
      freeDOF(1:Total_FD) = 0
      c_index = 0

      !Fix displacement in x-direction.
      do i = 1,Num_Bou_x
          cur_Node = Bou_x(i)
          c_index = c_index+1
	      fixedDOF(c_index) = 2*cur_Node-1
      end do

      !Fix displacement in y-direction.
      do i = 1,Num_Bou_y
          cur_Node = Bou_y(i)
          c_index = c_index+1
	      fixedDOF(c_index) = 2*cur_Node
      end do
      num_freeDOF = 0
      do i_fd =1,Total_FD
          if ( .not.(ANY( fixedDOF(1:c_index) .eq. i_fd)) ) then
              num_freeDOF = num_freeDOF +1
              freeDOF(num_freeDOF) = i_fd
          end if
      end do

      RETURN
      END SUBROUTINE Boundary_Cond


text:
SUBROUTINE RUNSPEC( FNAME, USEENV, SDATE, STIME, TSTEP, NRECS )

    !!***********************************************************************
    !! Version "$Id: runspec.f90 1 2017-06-10 18:05:20Z coats $"
    !! EDSS/Models-3 I/O API.
    !! (C) 2015-2016 UNC Institute for the Environment.
    !! Distributed under the GNU LESSER GENERAL PUBLIC LICENSE version 2.1
    !! See file "LGPL.txt" for conditions of use.
    !!.........................................................................
    !!  subroutine body starts at line  82
    !!
    !!  FUNCTION:
    !!      Return time step sequence SDATE:STIME:TSTEP:NRECS
    !!      compatible with FNAME as an input file.
    !!      If FNAME is blank, get default SDATE:STIME from wall-clock
    !!      and use
    !!      If USEENV, get SDATE, STIME, TSTEP, EDATE, ETIME from environment
    !!      Otherwise prompt the user for them.
    !!      NOTE:  May *NOT* assume that FNAME is already completed.
    !!      M3EXIT on failure.
    !!
    !!  PRECONDITIONS REQUIRED:
    !!      FNAME already opened.
    !!      If  USEENV:
    !!          setenv SDATE <starting date (YYYYDDD)>
    !!          setenv STIME <starting time  (HHMMSS)>
    !!          setenv TSTEP <time step      (H*MMSS)>
    !!          setenv EDATE <ending date   (YYYYDDD)>
    !!          setenv ETIME <ending date    (HHMMSS)>
    !!
    !!  REVISION  HISTORY:
    !!      Adapted 03/2015 by Carlie J. Coats, Jr., UNC IE, from algorithm
    !!      used in various "m3tools" programs
    !!      Bug-fix 05/2016 by CJC
    !!***********************************************************************

    IMPLICIT NONE

    INCLUDE 'PARMS3.EXT'
    INCLUDE 'FDESC3.EXT'
    INCLUDE 'IODECL3.EXT'

    !!...........   ARGUMENTS and their descriptions:

    CHARACTER(LEN=*), INTENT(IN   ) :: FNAME        !!  input file
    LOGICAL,          INTENT(IN   ) :: USEENV       !!  input file
    INTEGER,          INTENT(  OUT) :: SDATE        !!  starting date YYYYDDD
    INTEGER,          INTENT(  OUT) :: STIME        !!  starting time  H*MMSS
    INTEGER,          INTENT(  OUT) :: TSTEP        !!  time step      H*MMSS
    INTEGER,          INTENT(  OUT) :: NRECS        !!  Number of records

    !!...........   PARAMETERS and their descriptions:

    CHARACTER*1 , PARAMETER :: BLANK = ' '
    CHARACTER*16, PARAMETER :: PNAME = 'RUNSPEC'

    !!...........   EXTERNAL FUNCTIONS and their descriptions:

    LOGICAL, EXTERNAL :: ENVYN
    INTEGER, EXTERNAL :: ENVINT, GETNUM, JSTEP3, CURREC

    !!...........   SCRATCH LOCAL VARIABLES and their descriptions:

    INTEGER     SDATE1, JDATE1
    INTEGER     STIME1
    INTEGER     EDATE1
    INTEGER     ETIME1
    INTEGER     TSTEP1
    INTEGER     NRECS1
    INTEGER     JDATE, JTIME
    INTEGER     EDATE, ETIME
    INTEGER     ISTAT

    LOGICAL     EFLAG

    CHARACTER*256   MESG

    !!--------------------------------------------------------------
    !!   begin body of program RUNSPEC

    IF ( FNAME .EQ. BLANK ) THEN

        CALL GETDTTIME( JDATE1, STIME1 )
        SDATE1 = -1
        STIME1 = 10000 * ( STIME1 / 10000 ) !!  truncate to hour
        TSTEP1 = 10000
        EDATE1 = 9999999
        ETIME1 = 9999999

    ELSE IF ( .NOT.DESC3( FNAME ) ) THEN

        MESG = 'Could not get description for "'//TRIM( FNAME )//'"'
        CALL M3EXIT( PNAME, 0, 0, MESG, 2 )

    ELSE IF ( TSTEP3D .EQ. 0 ) THEN   !!  time independent case

        SDATE = SDATE3D
        STIME = STIME3D
        TSTEP = 0
        NRECS = 1
        RETURN

    ELSE

        SDATE1 = SDATE3D
        STIME1 = STIME3D
        TSTEP1 = TSTEP3D
        NRECS1 = MXREC3D
        CALL LASTTIME( SDATE1,STIME1,TSTEP1,NRECS1, EDATE1,ETIME1 )

    END IF

    EFLAG  = .FALSE.

    IF ( USEENV ) THEN

        SDATE = ENVINT( 'SDATE', 'Starting date (YYYYDDD)', SDATE1, ISTAT )
        IF ( ISTAT .GT. 0 ) THEN
            EFLAG = .TRUE.
            CALL M3MESG( 'Bad environment variable "SDATE"' )
        END IF

        STIME = ENVINT( 'STIME', 'Starting time (HHMMSS)', STIME1, ISTAT )
        IF ( ISTAT .GT. 0 ) THEN
            EFLAG = .TRUE.
            CALL M3MESG( 'Bad environment variable "STIME"' )
        END IF

        EDATE = ENVINT( 'EDATE', 'Ending date (YYYYDDD)', EDATE1, ISTAT )
        IF ( ISTAT .GT. 0 ) THEN
            EFLAG = .TRUE.
            CALL M3MESG( 'Bad environment variable "EDATE"' )
        END IF

        ETIME = ENVINT( 'ETIME', 'Ending time (HHMMSS)', ETIME1, ISTAT )
        IF ( ISTAT .GT. 0 ) THEN
            EFLAG = .TRUE.
            CALL M3MESG( 'Bad environment variable "ETIME"' )
        END IF

        TSTEP = ENVINT( 'TSTEP', 'Time step (H*MMSS)', TSTEP1, ISTAT )
        IF ( ISTAT .GT. 0 ) THEN
            EFLAG = .TRUE.
            CALL M3MESG( 'Bad environment variable "TSTEP"' )
        END IF

        IF ( EFLAG ) THEN
            CALL M3EXIT( PNAME,0,0, 'Bad environment', 2 )
        END IF

    ELSE

        CALL M3MESG( 'Now enter time step sequence parameters' )

        SDATE = GETNUM( SDATE1, EDATE1, SDATE1, 'Enter STARTING DATE (YYYYDDD)' )

        STIME = GETNUM( 0, 999999999,   STIME1, 'Enter STARTING TIME  (H*MMSS)' )

        TSTEP = GETNUM( TSTEP1, 999999999, TSTEP1, 'Enter     TIME STEP  (H*MMSS)' )

        EDATE = GETNUM( SDATE, EDATE1, EDATE1, 'Enter  ENDING DATE (YYYYDDD)' )

        ETIME = GETNUM( 0, 999999999,  ETIME1, 'Enter   ENDING TIME  (H*MMSS)' )

    END IF      !!  if useenv, or not

    IF ( FNAME .EQ. BLANK ) THEN
        CONTINUE
    ELSE IF ( JSTEP3(SDATE,STIME,SDATE1,STIME1,TSTEP1) .LT. 0 ) THEN
        EFLAG = .TRUE.
        WRITE( MESG, '(A, I9.7, A, I6.6, 2X, 3 A )' )       &
           'Starting date&time', SDATE, ':', STIME,  'not compatible with "', TRIM( FNAME ) , '"'
        CALL M3MESG( MESG )
    END IF

    JDATE = SDATE
    JTIME = STIME
    CALL NEXTIME( JDATE, JTIME, TSTEP )
    IF ( JSTEP3( JDATE,JTIME, SDATE,STIME,TSTEP ) .LT. 0 ) THEN
        EFLAG = .TRUE.
        WRITE( MESG, '(A, I10.6, 2X, 3 A )' ) 'Time step', TSTEP, 'not compatible with "', TRIM( FNAME ) , '"'
        CALL M3MESG( MESG )
    END IF

    NRECS = CURREC( EDATE, ETIME, SDATE, STIME, TSTEP, JDATE, JTIME )

    IF ( EFLAG ) THEN
        CALL M3EXIT( PNAME,0,0, 'Bad input time step sequence', 2 )
    END IF

    RETURN

END SUBROUTINE  RUNSPEC

max_stars_repo_name: mdhillmancmcl/TheWorldAvatar-CMCL-Fork

max_issues_repo_name: mdhillmancmcl/TheWorldAvatar-CMCL-Fork

max_forks_repo_name: mdhillmancmcl/TheWorldAvatar-CMCL-Fork