text:
! <mod_util.f90 - A component of the City-scale
!                 Chemistry Transport Model EPISODE-CityChem>
!*****************************************************************************!
!*
!* EPISODE - An urban-scale air quality model
!* ==========================================
!* Copyright (C) 2018  NILU - Norwegian Institute for Air Research
!*                     Instituttveien 18
!*                     PO Box 100
!*                     NO-2027 Kjeller
!*                     Norway
!*
!*                     Contact persons: Gabriela Sousa Santos - gss@nilu.no
!*                                      Paul Hamer - pdh@nilu.no
!*
!* Unless explicitly acquired and licensed from Licensor under another license,
!* the contents of this file are subject to the Reciprocal Public License ("RPL")
!* Version 1.5, https://opensource.org/licenses/RPL-1.5 or subsequent versions as
!* allowed by the RPL, and You may not copy or use this file in either source code
!* or executable form, except in compliance with the terms and conditions of the RPL.
!*
!* All software distributed under the RPL is provided strictly on an "AS IS" basis,
!* WITHOUT WARRANTY OF ANY KIND, EITHER EXPRESS OR IMPLIED, AND LICENSOR HEREBY
!* DISCLAIMS ALL SUCH WARRANTIES, INCLUDING WITHOUT LIMITATION, ANY WARRANTIES OF
!* MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, QUIET ENJOYMENT, OR NON-INFRINGEMENT.
!* See the RPL for specific language governing rights and limitations under the RPL.
!*
!* ==========================================
!* The dispersion model EPISODE (Grønskei et. al., 1993; Larssen et al., 1994;
!* Walker et al., 1992, 1999; Slørdal et al., 2003, 2008) is an Eulerian grid model
!* with embedded subgrid models for calculations of pollutant concentrations resulting
!* from different types of sources (area-, line- and point sources). EPISODE solves the
!* time dependent advection/-diffusion equation on a 3 dimensional grid.
!* Finite difference numerical methods are applied to integrate the solution forward in time.
!* It also includes extensions as the implementation of a simplified EMEP photochemistry
!* scheme for urban areas (Walker et al. 2004) and a scheme for Secondary Organic Aerosol
!* implemented by Håvard Slørdal
!*
!* Grønskei, K.E., Walker, S.E., Gram, F. (1993) Evaluation of a model for hourly spatial
!*    concentrations distributions. Atmos. Environ., 27B, 105-120.
!* Larssen, S., Grønskei, K.E., Gram, F., Hagen, L.O., Walker, S.E. (1994) Verification of
!*    urban scale time-dependent dispersion model with sub-grid elements in Oslo, Norway.
!*    In: Air poll. modelling and its appl. X. New York, Plenum Press.
!* Slørdal, L.H., McInnes, H., Krognes, T. (2008): The Air Quality Information System AirQUIS.
!*    Info. Techn. Environ. Eng., 1, 40-47, 2008.
!* Slørdal, L.H., Walker, S.-E., Solberg, S. (2003) The Urban Air Dispersion Model EPISODE
!*    applied in AirQUIS. Technical Description. NILU TR 12/2003. ISBN 82-425-1522-0.
!* Walker, S.E., Grønskei, K.E. (1992) Spredningsberegninger for on-line overvåking i Grenland.
!*    Programbeskrivelse og brukerveiledning. Lillestrøm,
!*    Norwegian Institute for Air Research (NILU OR 55/92).
!* Walker, S.E., Slørdal, L.H., Guerreiro, C., Gram, F., Grønskei, K.E. (1999) Air pollution
!*    exposure monitoring and estimation. Part II. Model evaluation and population exposure.
!*    J. Environ. Monit, 1, 321-326.
!* Walker, S.-E., Solberg, S., Denby, B. (2003) Development and implementation of a
!*    simplified EMEP photochemistry scheme for urban areas in EPISODE. NILU TR 13/2013.
!*    ISBN 82-425-1524-7
!*****************************************************************************!

      module mod_util

! ----------------------------------------------------------------------------------
! Based on:
! Version Episode 5.5 (May29, 2012) prepared for BB-Stand-Alone
! Original source code of EPISODE by Sam-Erik Walker (NILU)
!
! Sam-Erik Walker
! Norwegian Institute for Air Research (NILU)
! Instituttveien 18 P.O. Box 100
! N-2027 Kjeller, NORWAY
! Tel: +47 63898000 Fax: +47 63898050
! E-mail: sam-erik.walker@nilu.no
!
! ----------------------------------------------------------------------------------
! ----------------------------------------------------------------------------------
! REVISION HISTORY
!
!           2017  M. Karl: copied here the functions: EXP4, XLINE, YLINE (from hwylne.f90).
!                          Added functions func_canyon, qromb, trapzd, polint
!                          from Numerical Recipes
!
! ----------------------------------------------------------------------------------

      double precision :: avogad
      real :: busi
      real :: ctok
      real :: g
      real :: gamma
      real :: kappa
      real :: miss
      double precision :: pi
      real :: rad
      real :: rgas_j
      real :: charn

      parameter (avogad = 6.022045e23)
      parameter (busi   = 4.7)
      parameter (ctok   = 273.15)
      parameter (g      = 9.80665)
      parameter (gamma  = 0.27)
      parameter (kappa  = 0.41)
      parameter (miss   = -9900.)
      parameter (pi     = 3.141592654)
      parameter (rad    = pi/180.)
      parameter (rgas_j = 8.31441)
!MSK start
      parameter (charn  = 0.013)
!MSK end

!     avogad - Avogadros constant
!     busi   - Businger's surface layer theory constant
!     ctok   - Conversion from degrees Celcius to degrees Kelvin
!     g      - Acceleration of gravity (m/s^2)
!     gamma  - Venkatram constant (1984)
!     kappa  - Von Karman's constant
!     miss   - Missing data value
!     pi     - The mathematical constant pi
!     rad    - Conversion from degrees to radians
!     rgas_j - The molar gas constant (J/(mol*K))
!            - When pressure is measured in hPa = mb then rgas_j must be
!            - multiplied by 10^-2 since 1 Pa = 1 N/m^2
!     charn  - Charnock constant = 0.013; Units : none

      contains

!MSK copied function EXP4 into mod_util from EXP4.for
!   ******************************************************************

      REAL FUNCTION EXP4(X)

! The function computes exp(x) for all ranges of x

! Scalar arguments

      REAL X

! Parameters

      REAL XMIN,XMAX
      PARAMETER (XMIN = -87.3,XMAX = 88.7)

! Calculate exp(x)

      IF (X .GT. XMIN .AND. X .LT. XMAX) THEN
          EXP4 = EXP(X)
      ELSEIF (X .LE. XMIN) THEN
          EXP4 = 0.
      ELSE
          EXP4 = EXP(XMAX)
      ENDIF

      RETURN

! End of real function EXP4

      END FUNCTION EXP4
!   ******************************************************************

! MSK functions used in line source module hwylne.f90
!   ******************************************************************
!
! *** MODIFIED OCT. 1979 TO ADD ROMBERG INTEGRATION ENHANCEMENTS.
! *** MODIFIED SEP. 1989 TO IMPROVE LOCATION SEARCH
! *** XLINE is UPWIND DISTANCE OF R,S FROM RREC,SREC

      REAL FUNCTION XLINE(R,S,SINT,COST,RREC,SREC)

      REAL R,S,SINT,COST,RREC,SREC

!   -----------------------------------------------------------------
!
      XLINE = (R-RREC) *SINT + (S-SREC) * COST

      RETURN
      END FUNCTION XLINE

!   ******************************************************************
!
! *** MODIFIED OCT. 1979 TO ADD ROMBERG INTEGRATION ENHANCEMENTS.
! *** MODIFIED SEP. 1989 TO IMPROVE LOCATION SEARCH
! *** YLINE is CROSSWIND DISTANCE OF R,S FROM RREC,SREC

      REAL FUNCTION YLINE(R,S,SINT,COST,RREC,SREC)

      REAL R,S,SINT,COST,RREC,SREC

!   -----------------------------------------------------------------
!
      YLINE = (S-SREC) * SINT - (R-RREC) * COST

      RETURN
      END FUNCTION YLINE

!   ******************************************************************

!MSK functions used in line source module csubl.f90

!   ******************************************************************
!
! *** STREET CANYON
! *** Function in integral(a,b)=func(x)dx
! *** The function is not used because the integrals in OSPM can be
! *** solved analytically.

      real function func_canyon(x,const)

      REAL x,const

!   -----------------------------------------------------------------
!
      func_canyon = 1.0 / (x + const)

      return
      end function func_canyon

!   ******************************************************************

!MSK numerical recipe subroutines

      subroutine qromb(const, a, b, ssv)
    !----------------------------------------------------------------------
    !
    !****
    !
    !      purpose
    !      -------
    !      Romberg Integration of order 2K using polynomial extrapolation
    !
    !      interface
    !      ---------
    !
    !
    !      method
    !      ------
    !      NUMERICAL RECIPES
    !      Returns as ssv the integratl of the function func from a to b.
    !      Integration is performed by Romberg's method of order 2K, where
    !      e.g. K=2 is Simpson's rule.
    !      Parameters:
    !      EPS is the fractional accuracy desired, as determined by the
    !      extrapolation error estimate;
    !      JMAX is the total number of steps
    !      K is the number of points used in the extrapolation
    !      Uses:
    !      subroutine polint
    !      subroutine trapzd
    !
    !      external
    !      --------
    !      none
    !
    !      reference
    !      ---------
    !      Book Numerical Recipes, Ch. 4.3
    !
    !
    !------------------------------------------------------------------
    implicit none

    real, intent(in)                     :: a,b
    real, intent(in)                     :: const
    real, intent(out)                    :: ssv

    real, parameter                      :: EPS   = 1.e-6
    integer, parameter                   :: JMAX  = 20
    integer, parameter                   :: JMAXP = JMAX+1
    integer, parameter                   :: K     = 5
    integer, parameter                   :: KM    = K-1

    real, dimension(JMAXP)               :: h
    real, dimension(JMAXP)               :: s
    integer                              :: j
    real                                 :: dss

    ! Arrays h and s store the succesive trapezoidal approximation
    ! and their relative step sizes.

          h(1) = 1
          !print *,'qromb const',const
          do j=1, JMAX
            call trapzd(const,a,b,s(j),j)
            if (j.ge.K) then

    ! Call polint(xa,ya,PMAX,x, y,dy), polynomial extrapolation
              call polint(h(j-KM),s(j-KM),K,0., ssv,dss)

              if (ABS(dss).le.EPS*ABS(ssv)) return

            endif
            s(j+1) = s(j)
    ! This is a key step: The factor is 0.25 even though the step size
    ! is decreased by only 0.5. This makes the extrapolation a polynomial
    ! in h^2 as allowed in equation (4.2.1) [in NR], not just a polynomial in h.
            h(j+1) = 0.25 * h(j)
          enddo
    ! Should not get here, just for testing
          print *,'Warning: Too many steps in qromb (mod_util)'

      end subroutine qromb

      subroutine trapzd(const, a, b, s, n)
    !----------------------------------------------------------------------
    !
    !****
    !
    !      purpose
    !      -------
    !      Extended trapezoidal rule
    !
    !      interface
    !      ---------
    !
    !
    !      method
    !      ------
    !      NUMERICAL RECIPES
    !      Routine computes the nth stage of refinement of an extended
    !      trapezoidal rule.
    !      Input func is the name of the function to be integrated with
    !      limits a and b, also input. When called with n=1 the routine
    !      returns as s, the crudest estimate of integral(a,b)f(x)dx.
    !      Subsequent calls with n=2, 3,... (in that sequential order) will
    !      improve the accuracy of s, by adding 2^(n-2) additional interior
    !      points. s should NOT be modified between the subsequent calls
    !
    !      external
    !      --------
    !      none
    !
    !      reference
    !      ---------
    !      Book Numerical Recipes, Ch. 4.2
    !
    !
    !------------------------------------------------------------------
    implicit none

    real, intent(in)                     :: a,b
    real, intent(in)                     :: const
    integer, intent(in)                  :: n
    real, intent(out)                    :: s

    integer                              :: it
    integer                              :: j
    real                                 :: del
    real                                 :: summe
    real                                 :: tnm
    real                                 :: x

          if (n.eq.1) then
            s = 0.5 * (b-a) * ( func_canyon(a,const) + func_canyon(b,const) )
          else
            it=2**(n-2)
            tnm=REAL(it)
    ! This is the spacing of the points to be added
            del=(b-a)/tnm
            x=a+0.5*del
            summe=0.
            do j = 1, it
              summe=summe+func_canyon(x,const)
              x=x+del
            enddo
    ! This replaces s by its refined value
            s = 0.5 * ( s+(b-a)*summe/tnm )
          endif

          return

      end subroutine trapzd

      subroutine polint(xa,ya,PMAX,x, y,dy)
    !----------------------------------------------------------------------
    !
    !****
    !
    !      purpose
    !      -------
    !      Polynomial Interpolation and Extrapolation
    !
    !      interface
    !      ---------
    !
    !
    !      method
    !      ------
    !      NUMERICAL RECIPES
    !      Given arrays xa and ya, each of length n, and given a value x,
    !      this routine returns a value y, and an error estimate dy.
    !      If P(x) is the polynomial of degree N-1 such that
    !      P(xai) = yai, i = 1, ... , n, then the returned value y = P(x).
    !
    !      external
    !      --------
    !      none
    !
    !      reference
    !      ---------
    !      Book Numerical Recipes, Ch. 3.1
    !
    !
    !------------------------------------------------------------------
    implicit none

    REAL, intent(in)                     :: x
    integer, intent(in)                  :: PMAX

    REAL, dimension(PMAX), intent(in)    :: xa,ya
    REAL, intent(out)                    :: y
    REAL, intent(out)                    :: dy

    REAL                                 :: den,dif,dift,ho,hp,w

    integer, parameter                   :: NMAX=10
    real, dimension(NMAX)                :: c
    real, dimension(NMAX)                :: d
    integer                              :: i,m,ns

          ns=1
          dif=abs(x-xa(1))

          ! Here we find the index ns of the closest table entry
          do i=1,PMAX
            dift=abs(x-xa(i))
            if (dift .lt. dif) then
              ns=i
              dif=dift
            endif
            ! initialize the tableau of c's and d's.
            c(i)=ya(i)
            d(i)=ya(i)
          end do

          ! This is the initial approximation to y.
          y=ya(ns)
          ns=ns-1

          ! For each column of the tableau
          ! we loop over the current c's and d's and update them.
          ! After each column in the tableau is completed, we decide
          ! which correction, c or d, we want to add to our accumulating           CRV(IC) = CSUM
          ! value of y, i.e., which path to take through
          ! the tableau forking up or down. We do this in such a
          ! way as to take the most "straight line" route through the
          ! tableau to its apex, updating ns accordingly to keep track
          ! of where we are. This route keeps the partial approximations
          ! centered (insofar as possible) on the target x.
          ! The last dy added is thus the error indication.
          do m=1,PMAX-1
            do i=1,PMAX-m
              ho=xa(i)-x
              hp=xa(i+m)-x
              w=c(i+1)-d(i)
              den=ho-hp
              ! if(den.eq.0.)pause 'failure in polint'
              ! This error can occur only if two input xa's are (to within roundoff) identical.
              if (den .eq. 0.0) then
                den= den+1.e-32
              endif
              den=w/den
              d(i)=hp*den
              c(i)=ho*den
            end do
            if (2*ns .lt. PMAX-m) then
              dy=c(ns+1)
            else
              dy=d(ns)
              ns=ns-1
            endif
            y=y+dy
          end do

      end subroutine polint

!   ******************************************************************

      subroutine FreeUtilMemory()

!DEC$ ATTRIBUTES DLLEXPORT,STDCALL,ALIAS:'FreeUtilMemory' :: FreeUtilMemory

      implicit none

!     End of subroutine FreeUtilMemory
      end subroutine FreeUtilMemory

      end module mod_util

max_stars_repo_name: mdhillmancmcl/TheWorldAvatar-CMCL-Fork

max_issues_repo_name: mdhillmancmcl/TheWorldAvatar-CMCL-Fork

max_forks_repo_name: mdhillmancmcl/TheWorldAvatar-CMCL-Fork


text:
program strings

  implicit none

  character(16) :: mystring1
  character(len=16) :: mystring2
  character :: mystring3(16)

end program strings

max_stars_repo_name: EraYaN/transpyle

max_issues_repo_name: EraYaN/transpyle

max_forks_repo_name: EraYaN/transpyle


text:
Poot & Scoot is a licensed Plumbers plumber Lic.# 695378 that does plumbing repairs, sewer and drain cleaning, root control and video inspections of sewage lines. USE Credit Union members receive a 25% discount.  The comments left for him overwhelmingly praise his honesty, prices, and skills as a plumber. He is always on call!

They have no known relationship to Pootleberry Dinglehopper this poot.

Other Plumbers

20071115 20:20:38 nbsp David is a super nice guy, and Davis resident. Prior to repairing our clogged sewer line, he advised me to call the City of Davis Dept. of Public Works to locate the cleanout pipe in our front yard. This ended up saving us a bundle of money. The actual work he did on our house was professional and his prices were reasonable.  Users/JanieBooth

20080203 22:19:05 nbsp A good friend from Lyon Real Estate recommended Poot & Scoot to me. Dave, from Poot & Scoot stopped by to give me some advice on how to fix my sewer line after Plumbing M.D. damaged it. He went out of his way to help me out. He was very pleasant to deal with.  Users/CarolBlake

20080307 11:36:31 nbsp David is wonderful!!! I had a broken garbage disposal and he kindly helped me through it. Im telling all of my friends about Poot and Scoot. Hes kind, honest, and patient! Call him! Users/TomLee

20080307 12:23:26 nbsp Im sorry, but this may be one of the most amusing names ever for a business. Users/JabberWokky

20080320 19:39:13 nbsp David is great.  Did a great job, good prices, very helpful and knowledgeable. Users/HarrySpanglet

20080401 13:19:11 nbsp I called David after having problems trying to solder a pipe in the crawl space under my house.  The price he quoted me to come over was EXTREMELY reasonable.    Then, he asked me questions about the job and gave me all sorts of tips on how to get it done right.  He said, try them out, then call me if you still have a problem.  I did what he said, and it worked; I ended up with free advice!  He is a really nice guy, very knowledgeable and clearly very generous too.  Next job I have, I will definitely call him.   Users/JosephBiello

20080624 22:38:18 nbsp David is the best plumber in Davis.  Work is very high quality, reasonably priced and stayed within quoted value.  He even found and fixed an unnoticed leak behind the same wall that was probably the reason for termite & ant problems in that part of the house. Users/VincentFox

20080715 10:22:06 nbsp David has got to be the best plumber Ive ever come across.  We had a clogged sewer line, and he came by immediately to check things out.  He mentioned I should try to call the City of Davis Public Works, to get them to come check on their part of the line, which they would fix, if something was wrong.  The problem was fixed in an hour, for FREE by Public Works.  David told me, I hate to talk myself out of business, but I hate it even more when people arent informed abouthow things work.  Absolutely amazing.  He is the only person Ill call from here on out.  Thank You again David! Users/CarDogg13

20081008 10:26:58 nbsp Dave, this link is for you:  http://www.nytimes.com/2008/10/09/us/09water.html?hp Users/JosephBiello

20081008 11:04:26 nbsp the poot and scoot truck is pretty nifty looking, I see it around town all the time Users/StevenDaubert

20081024 14:22:08 nbsp I was extremely happy with the service provided by Poot & Scoot.  David was very helpful and flexible with scheduling.  He recognized the problem right away, bought the parts, and fixed it.  I would highly recommend Poot and Scoot for plumbing services.  Very reliable and trusting company. Users/jessicakk

20081216 12:53:07 nbsp David is excellent. Twice this year (once at the place we were renting in Davis, once at our new place in Woodland) we called him and had problems fixed in less than 3 hours after our initial call. He is friendly and extremely reasonable with his pricing! Use him for any plumbing problem. Users/Adie

20090116 09:38:17 nbsp I would recommend Poot & Scoot.  I had a good experience with David taking care of a clogged toilet.  He is better than the big plumbing companies. Users/jakedavis

20090608 14:37:44 nbsp This guy is a Davis legend; came over within 10 mins of my call and figured out the problem on the spot (something that had previously eluded the city of Davis public works dept, incidentally). Look no further.    Users/drhanszarkov

20090926 00:24:26 nbsp Great experience. My sink wouldnt drain, so I looked up plumbers on Daviswiki. After reading reviews, I decided to call Poot & Scoot. David, the owner, answered and first walked me through a series of techniques to unclog the sink myself. Turns out the clog was deeper than the sink pipes, so it didnt work. But I learned a lot. He fixed it that night, in 1.5 hrs. (He charges no more than 2 hrs for sink problems.) Im a poor student, but he was flexible with payment too. Users/Devonian

20090929 17:14:10 nbsp Dave is easy to find, comes promptly if he can, and is always happy to talk you through the problem so you can fix some of it yourself, or make adjustments after he leaves.  Very good to work with.   Users/TedBuehler

20100125 11:10:05 nbsp MUST USE!
    My roommates and I just had the sort of experience that makes me feel the world can be a fair place, and honor is still important to some people. David White, owner of Poot & Scoot plumbing in Davis, CA, 5307567883, is an honest man, and he fixed our plumbing problem at a very reasonable cost. He explained everything, recommended a $10.00 tool and a $1.99 bottle of Palmolive for next time, and grew incensed that another plumber had just given us a $400.00 quote for the job he accomplished for $95.00.  Users/Leah

20100125 17:52:06 nbsp Like Leah, Dave has restored my faith in the goodness and honesty of some people. I had a clogged drain since Saturday and tried everything in the book to clear it to no avail. Dave came over, spent about an hour here, explained possible causes, and drain unclogged. Much to my surprise he only charged 1/2 hours time! If you need a plumber and want an honest and fair person, Daves your man! Users/Kiefer

20100628 20:32:26 nbsp Look no further for a plumber.  Hes a gem. Users/stephthrasher

20100628 23:01:14 nbsp Dave is a great guy, I called him for toilet questions and he did some research to find the best toilet for the price that would fit my unusually small bathroom.  He told me where to buy it so I could install it myself.  I really appreciate the helpfulness without a sales pitch.  Users/DagonJones

20100716 18:31:54 nbsp Best plumber I ever had!  Dave installed a water line for our refrigerator and did a beautiful job.  Also once came to my house to fix my faucet handle and didnt charge for it!  He is super friendly and absolutely competent. Users/EstherWright

20100924 17:24:27 nbsp After reading all the good reviews here, I decided to call Poot & Scoot. It was Sunday morning, David picked up the phone and said he would come over sometime in the evening and he would call me to confirm. Sure enough, David came over in the evening and fixed the problem we had with our dishwasher. He charged us for a half hours time. We were extremely happy with the service provided. Users/EastDavisResident

20101129 18:08:15 nbsp I called David the AM to see if he could come and fix my water heater... He asked me some pertinent questions, and told me what the issue likely was,  could come late afternoon..., and that it would take 1.25 hours to fix andplus parts, and what his rate was.  He came on time, fixed the water heater, and finished quickly, in the time he estimated!  Couldnt ask for a better experience unless it was free!  I wish!  Great plumber for small repairs! Users/Carsonsig

20101225 10:20:22 nbsp Incredible service & very reasonable price! We first tried to repair a leaking valve on the toilet, misunderstood the connection, and mangled our pipe. He fixed everything. A Davis hero!  Users/LynnEpstein

20101231 16:17:03 nbsp Highest possible recommendation. I had a water leak in my main line on Christmas Eve and based on recommendations from this Wiki, I called Dave. He managed to get over and employ some clever diagnostic techniques to help me validate where my leak was coming from. He worked with me to minimize my cost (I did a lot of the digging), and I probably wound up talking with him at least a half a dozen times, and he made probably half a dozen trips to my house to get the thing finally fixed. He let me borrow several of his tools, which were extremely timesaving. He gave me lots of tips. He was up front on all his estimates (didnt charge for the diagnostics work or all the consults) and did a real nice job on the actual pipe repair. As honest, humble, helpful, straighforward a guy as I can imagine. He doesnt do certain kinds of work (acoustics work to find underground leaks, gas lines, water heaters) but for line repairs, toilets, etc., I give Dave the highest possible recommendation. In the end, I would have been satisfied paying twice or thrice what I wound up paying him for his services (and Im picky). The stereotypical fears about fees, etc. may apply to many plumbers, but Daves service was exceptional and his cost was far below what I had steeled myself to expect. And I always like helping out a local, small biz operation... Users/PaulCrawford

20110126 11:08:19 nbsp We are lucky to have someone like David living here in Davis.

My experience: Our house had a backed up drain. After spending hours snaking it myself to no avail, I called David Smith. Even though he was busy, he talked to me over the phone about the problem in detail to see if there was some way I could fix it myself. He then told me that he was too busy to come over, but recommended another plumber else who could do it. When that plumber could not resolve the problem, he came over at 6PM that night to fix the problem. (He called us to tell us that he felt bad that the other guy couldnt do it.)

He told us it would probably take two hours, but the problems was a real bear, and it took him almost four hours of snaking through various cleanouts before he finally found the clog and got it out. He charged us only for 2 hours. He even included the replacement of the kitchen drain trap at no charge.

This was definitely service well above and beyond the call. I had used him once before for a much more routine problem, and also had a great experience with him.

David is also an incredibly nice man. I will definitely call him for any other plumbing issues I have. Users/MarkusLuty

20110301 20:16:10 nbsp You cant do any better than Dave.  His prices are reasonable, he did an excellent and thorough job, and most importantly, he is a true professional who takes pride in his work  even for a relatively small job.  Plus, hes got the best business name of all the plumbers.   Users/CovertProfessor

20110301 21:32:15 nbsp Dave was great. He replaced our kitchen sink, installed new faucet and garbage disposal for 1/3 the cost of another plumbing company in Davis. He was straightforward with the costs, and was stunned with the itemized receipt the other company wanted to charge me. He explained how to maintain everything, and was very easy to talk to. Its nice to be able to find a trustworthy person who gets the job done right and at a fair cost. Users/AMoore

20110328 18:05:11 nbsp David is as honest & straightfoward a guy as they come. I call Poot & Scoot when in need of a plumber because David is skilled as a plumber, his prices are reasonable, and hes very trustworthy...he gets my highest recommendation! Users/RickHolmes

20110526 22:28:41 nbsp I called David today at 4:30 and left him a voicemail; at 8pm he had me, the happy customer, gladly handing him a check for my fixed toilet. When he called me back and left me a voicemail, he first suggested that if I was having a problem with any other drains in the house that I should call Public Works and they would fix it for me for free. He left me detailed directions and provided me with the phone numbers, too! But, since that wasnt the case, he came right out, spent about 2 hours at my place (exactly as he had quoted), and pulled a tree root out of my drain pipe! Now it flushes beautifully. As he was wrapping up, he mentioned to us that the city of Davis is offering rebates on lowflow toilets and recommended some good models that he had seen at a good price in case we were interested. All in all, I was really impressed with his timeliness, honesty and eagerness to keep us informed. Users/StephanieStudebaker

20110601 12:24:19 nbsp David is awesome. I called to see if he could come fix an issue I was having with my sink faucet. Rather than take my money, he guided me through the repair over the phone! I cant imagine calling any other plumber. Users/DanLawyer

20110805 16:32:30 nbsp As others have observed, David is great.  Hes the third plumber weve used, and hell be our goto guy now.  He fixed two complicated leaky faucets that required ordering parts.  He actually ordered extra parts and showed me how to fix it myself in case the other faucet started leaking.  Compared to the constant upsell of other plumbers and inconsistent service, David is honest, friendly and local. Users/ScottLay

20110905 17:01:26 nbsp I searched online for a plumber in Davis and I read all these wonderful five stars reviews about David on Daviswiki and Yelp.com.  So I left David a message on Sunday night around 10pm.  Then David returned my call 8am on Monday (labor day)morning.  He brought all the parts and fixed our problem in the afternoon.  Hes really as nice and honest as the reviews said.  I really appreciate him for coming out even on a holiday.  He only charged us the labor.  We paid for the parts but he didnt charge extra for the parts which he brought from the store.  I will definitely recommend him to all my friends. Users/AmyC

20110914 11:09:20 nbsp Dave was wonderful!!!  Came the same day (within 2 hrs) to examine our leaking shower faucet, figured out the part that was needed, gave an estimate, went to Roseville to pick up the part, came back first thing the next morning and repaired it very quickly.  Wonderful plumber, and will be using him again! Users/ChristineHaas

20110916 08:56:15 nbsp David is as great as all the reviews say he is. I called him midday on Wednesday and he came at 6pm that evening to work on our leaky toilet. He didnt want us to have to take time off of work, and was more than willing to come after hours. Even though he spent extra time with us talking about the toilet and giving us tips (and just chatting! such a nice fellow!) he roundeddown on the time he spent working when he wroteup the bill. Hes very honest and good at what he does. Users/EmiliaC

20111011 10:53:14 nbsp Highly recommended.  Called David at 6PM about a clogged drain, and he arrived 15 minutes later.  Worked for nearly 4 hours on a difficult problem, did a great job and charged exactly the price he quoted over the phone despite the problems that cropped up during repair.   Users/SteveNoctor

20111011 14:16:02 nbsp David came out to fix the clogged bathtub  we tried fixing it ourselves but the clog was really deep down and we didnt have the equipment to go that far  I assumed Id have to pay for the full hour but he only charged me for the half hour it took for him to fix it! That was a nice surprise! Users/Ravyn

20111023 18:40:05 nbsp Another satisfied customer. David did a great job fixing a challenging leak. His work was high quality. His rate was reasonable. And he completed the work on the weekend. He’s a true professional. Highly recommended! Users/AngeloDeSantis

20111130 09:36:17 nbsp I completely agree with the positive reviews of Poot & Scoot, and the generally negative reviews of Plumbing M.D. on the wiki.

When my drain was blocked my wife pulled Plumbing MDs number out of a hat. They responded quickly, however they claimed they ran 100 feet of drain snake and did not encounter the blockage.  They charged us $125 for an hour of work without any noticeable effect, then tried to sell my wife on a $400 service plan whereby we would receive a %10 discount on future work done and get priority service.  In other words, in the future it would be easier and faster for me to pay them to not accomplish anything for a reasonable $400.  They also tried to sell a $300 video inspection of our pipes  however I understand from other plumbers that when the drain is blocked you are unlikely to see anything on the camera.  It definitely felt like they were trying to squeeze us for as much money as they could.

In contrast, once we took the time to ask friends and read some reviews about Poot and Scoot, we wish we had called him in the first place.  Dave lived up to the accolades by responding cheerfully on a Saturday morning, letting himself in as I was at work, clearing the blockage in the drain (which was only 70 feet in  unknown why Plumming MD couldnt fix it with 100 feet of snake), and leaving a clean work area and a reasonable bill on the step on his way out.

Also in stark contrast with Pumming MD: Dave will not charge if he cant fix your problem, and guarantees his work for 6 months (he will come back out at no cost if the drain gets plugged again). This shows his confidence in his own ability (he says he has only been stumped by a drain twice in 6 years), and IMHO is all you need to know when selecting a new plumber.

Please support excellent business practices and local small business and pick Poot & Scoot as your next plumber.  Users/jamesrbronson

20111212 18:54:27 nbsp We decided to go with Poot & Scoot because of the rave reviews on the Wiki, and we were not disappointed! David started by suggesting we contact the city, as the problem might be theirs  allowing us to fix the problem with no cost at all. Unfortunately, that was not the problem. He then came out the very next day — apologizing he could not make it out at 8pm that very (Sunday) night, fixed the problem at an extremely affordable and surprisingly low cost. He then gave detailed instructions on how to install a new pipe ourselves so we would not have to pay for any additional work — basically giving free advice. For any plumbing needs we ever have again, we will absolutely be calling Poot & Scoot!  Users/HeatherMcKibben

20120110 15:39:58 nbsp Excellent, prompt service; normal industry rates. As noted above, he first recommended having the city do a sewer line check on their side. They did; it wasnt their problem, so he was out within the hour. Explained the work in detail, gave an itemized bill, and discussed future issues. Highly recommended.  Users/DonShor

20120110 15:43:17 nbsp Every time I see the Poot & Scoot vehicle around town I laugh.  In a good way!  I havent needed a plumber in Davis, but if I do I know who Ill call.  Talk about a memorable name. Users/TomGarberson
 I spoke with David over the phone once at like 2 AM, he was super chill ! Daubert

20120129 11:02:41 nbsp The garbage disposal stopped working this morning (motor humming but no action) and despite finding instructions online on how to fix the problem I was still a bit too timid to follow the instructions since Im just a lowly tenant. I called Poot and Scoot (on a Sunday morning mind you) and within 10 minutes David called me back and walked me through how to fix the problem myself (free for me, but took up Davids time). Ten thumbs up!  Users/StellaCopeland

20120830 14:41:40 nbsp I cannot speak highly enough of Dave! Weve used him several times to help us out! We have old pipes and roots that interfere with our cleanout and he always knows exactly how to fix it. He only ever charges us for an hour and a half when sometimes he is working on our clogs for over 2 hours! He is friendly and generous and I will never use another plumber! Use Poot and Scoot!! Users/AllisonOlson

20121122 09:37:08 nbsp   Exemplary local business.  Called with a drain problem.  Rather than make a weekend/holiday high dollar trip to take care of it, he gave me instructions on how I might be able to do the work myself.  I was.  Thanks to Poot & Scoot and to Davis Wiki for putting me on to him. Users/ButchySchwank

20130107 10:08:02 nbsp We found David to be very helpful, fair, and honest. I think he gave us a free quote, didnt charge us for a cheap part that he had on hand, etc. However, after fixing our problem with our shower, he left us with a big hole in our bathroom wall because he said he didnt do sheetrock repair, and left us to find someone to fix that. That was very frustrating. I also just called him to see if he could install a new hot water heater for me, and he said that he doesnt do that. So it seems like hes great at what he does, but what he does is somewhat limited, so if I were to hire him again Id want to be sure that he would be able to complete the job or would take care of outsourcing the stuff he cant do. Users/harvardslacker

20130311 12:26:41 nbsp Great service for a fair price.  Called him this morning, we made an appointment for a few hours later, he arrived on time and fixed my sad kitchen drain quickly and without a mess.  Charged me for only 40 minutes of time total and gave me some tips to avoid clogs in the future. Users/SarahmDavis

20130423 14:17:59 nbsp We have enjoyed Davids service on multiple occasions.  Davis is very polite and always willing to help.  He is very honest about what he can do and what he cannot.  For the jobs that he was willing to take on, he was great in terms of quality of work and fees charged.  Yesterday we had a water hammer issue in the house and called him in the evening.  He offered some suggestions over the phone and, after knowing that we were not able to perform certain tasks, came over (after 8pm) and helped us resolve the issue.  Afterwards he did not want to take any payment for his help.  After my insistence, he asked for 20 dollars.  We are really impressed by his kindness and character.   Users/LW

20140410 11:36:38 nbsp I saw David last year and he cleaned my sewage blockage so that it worked again when we flushed the toilets. I was satisfied, but he thought the problem was only a temporary fix so he said No charge even though he spent an hour and a half on it. I wrote a check anyway and had to keep offering it to him before he finally accepted it. I rarely run into people with such strong morals and superior work ethics. Users/ChrisOkusako

20140416 22:41:46 nbsp We had an after hours emergency and called David.  He got back to us fairly quickly and once we explained the problem, he advised us to call the Department of Public Works.  Turns out, it was a city problem….not a problem with our home plumbing, so we didnt have to pay for it.  He saved us from having to pay for an after hours visit which was an honest thing to do.  The best part was that he assured us that if the Department of Public Works didnt fix the problem, he would come out and fix it even though it was after 8pm.  We would call him again without a doubt!   Users/darcid4

20140613 15:30:57 nbsp The Best, period. Dave came out late evening to fix a problem I couldnt, and had it solved in 1/2 hour. He was only going to charge $45 or so but I gave him all the cash I had on hand ($60$70) and I would have paid double what he asked for if I had it.  Users/RalphFinch

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text: #redirect Users/FreddieOakley


text:
!
!  Written by Leandro Martínez, 2009-2011.
!  Copyright (c) 2009-2018, Leandro Martínez, Jose Mario Martinez,
!  Ernesto G. Birgin.
!
!-----------------------------------------------------------------------------
!
! http://m3g.iqm.unicamp.br/packmol
!
! Usage (see the page above for further information):
!
! ./packmol < inputfile.inp
!
! References:
!
! L. Martinez, R. Andrade, E. G. Birgin, J. M. Martinez,
! PACKMOL: A package for building initial configurations for
! molecular dynamics simulations, J. Comp. Chem. 30:2157-2164, 2009.
!
! J. M. Martinez and L. Martinez,
! Packing optimization for the automated generation of complex
! system's initial configurations for molcular dynamics and
! docking. J. Comp. Chem. 24:819-825, 2003.
!
! This version of Packmol uses the optimization method GENCAN which
! is a part of the TANGO (Trustable Algorithms for Nonlinear General
! Optimization) project.
! Reference:
! E. G. Birgin, J. M. Martinez, Comp. Opt. Appl. 23:101-125, 2002.
! http://www.ime.usp.br/~egbirgin/tango
!
!

program packmol

  use sizes
  use compute_data
  use input
  use usegencan
  use flashsort
  use swaptypemod
  use ahestetic
  implicit none

  integer :: itype, irest, idatom, iatom
  integer :: idtemp, nmtemp, natemp, input_itypetemp
  integer :: linesttmp1, linesttmp2, jtype
  integer :: ntmol, n, iftype, icart, imol, iicart, iline_atoms
  integer :: i, iline, iiatom, iat, iirest, iratcount, ival
  integer :: loop
  integer :: resntemp, nloop_tmp
  integer :: strlength, ioerr

  double precision, allocatable :: x(:), xprint(:) ! (nn)
  double precision :: v1(3),v2(3),v3(3)
  double precision :: radscale, value
  double precision :: cmx, cmy, cmz, beta, gama, teta
  double precision :: xtemp, ytemp, ztemp
  double precision :: fx, bestf, flast, fprint, all_type_fx
  double precision :: fimp, fimprov
  double precision, parameter :: pi=4.d0*datan(1.d0)

  real :: etime, tarray(2), time0

  character(len=200) :: record, restart_from_temp, restart_to_temp
  character(len=80) :: xyzfile
  character(len=1) :: chain_tmp

  logical :: fixtmp
  logical :: rests
  logical :: movebadprint
  logical :: changechains_tmp

  logical, allocatable :: fixed(:) ! ntype

  ! Printing title

  call title()

  ! Set dimensions of all arrays

  call setsizes()

  ! Allocate local array

  allocate(fixed(ntype),x(nn),xprint(nn),xfull(nn))

  ! Start time computation

  time0 = etime(tarray)

  ! Reading input file

  call getinp()

  ! Put molecules in their center of mass

  call cenmass()

  ! Writting some input data

  write(*,*) ' Total number of atoms: ', ntotat

  ! Put fixed molecules in the specified position

  do itype = 1, ntype
    fixed(itype) = .false.
  end do

  do irest = 1, nrest
    if(ityperest(irest).eq.1) then
      do itype = 1, ntype
        if(irestline(irest).gt.linestrut(itype,1).and.&
           irestline(irest).lt.linestrut(itype,2)) then
          cmx = restpars(irest,1)
          cmy = restpars(irest,2)
          cmz = restpars(irest,3)
          beta = restpars(irest,4)
          gama = restpars(irest,5)
          teta = restpars(irest,6)

          ! Compute rotation matrix from euler angles

          call eulerfixed(beta,gama,teta,v1,v2,v3)

          idatom = idfirst(itype) - 1
          do iatom = 1, natoms(itype)
            idatom = idatom + 1
            xtemp =   coor(idatom,1)*v1(1) &
                    + coor(idatom,2)*v2(1) &
                    + coor(idatom,3)*v3(1)
            ytemp =   coor(idatom,1)*v1(2) &
                    + coor(idatom,2)*v2(2) &
                    + coor(idatom,3)*v3(2)
            ztemp =   coor(idatom,1)*v1(3) &
                    + coor(idatom,2)*v2(3) &
                    + coor(idatom,3)*v3(3)
            coor(idatom, 1) = xtemp + cmx
            coor(idatom, 2) = ytemp + cmy
            coor(idatom, 3) = ztemp + cmz
          end do
          record = name(itype)
          write(*,*) ' Molecule ',record(1:strlength(record)),'(',itype,') will be fixed.'
          fixed(itype) = .true.
          if(nmols(itype).gt.1) then
            write(*,*)' ERROR: Cannot set number > 1',' for fixed molecules. '
            write(*,*) '       Structure: ', itype,': ', trim(adjustl(record))
            stop
          end if
          if ( restart_from(itype) /= 'none' .or. &
               restart_to(itype) /= 'none' ) then
            write(*,*) ' ERROR: Restart files cannot be used for fixed molecules. '
            write(*,*) '        Structure: ', itype,': ', trim(adjustl(record))
            stop
          end if
        end if
      end do
    end if
  end do

  ! Reseting parameters for removing the fixed molecules
  ! fix is the logical variable that informs that there are fixed molecules

  fix = .false.
  ntemp = 0
  do itype = 1, ntype

    ! input_itype and fixedoninput vectors are used only to preserve the
    ! order of input in the output files

    input_itype(itype) = itype
    if(fixed(itype)) then
      fix = .true.
      fixedoninput(itype) = .true.
    else
      ntemp = ntemp + 1
      fixedoninput(itype) = .false.
    end if
  end do
  ntfix = ntype
  ntype = ntemp

  do i = 1, ntfix - ntype
    do itype = 1, ntfix - 1
      if(fixed(itype)) then
        record = name(itype)
        restart_to_temp = restart_to(itype)
        restart_from_temp = restart_from(itype)
        fixtmp = fixed(itype)
        idtemp = idfirst(itype)
        input_itypetemp = input_itype(itype)
        nmtemp = nmols(itype)
        natemp = natoms(itype)
        resntemp = resnumbers(itype)
        if(pdb) xyzfile = pdbfile(itype)
        linesttmp1 = linestrut(itype,1)
        linesttmp2 = linestrut(itype,2)
        changechains_tmp = changechains(itype)
        chain_tmp = chain(itype)
        nloop_tmp = nloop_type(itype)
        jtype = itype + 1
        if(.not.fixed(jtype)) then
          name(itype) = name(jtype)
          name(jtype) = record(1:10)
          restart_to(itype) = restart_to(jtype)
          restart_to(jtype) = restart_to_temp
          restart_from(itype) = restart_from(jtype)
          restart_from(jtype) = restart_from_temp
          idfirst(itype) = idfirst(jtype)
          idfirst(jtype) = idtemp
          input_itype(itype) = input_itype(jtype)
          input_itype(jtype) = input_itypetemp
          fixed(itype) = fixed(jtype)
          fixed(jtype) = fixtmp
          nmols(itype) = nmols(jtype)
          nmols(jtype) = nmtemp
          natoms(itype) = natoms(jtype)
          natoms(jtype) = natemp
          resnumbers(itype) = resnumbers(jtype)
          resnumbers(jtype) = resntemp
          changechains(itype) = changechains(jtype)
          changechains(jtype) = changechains_tmp
          chain(itype) = chain(jtype)
          chain(jtype) = chain_tmp
          nloop_type(itype) = nloop_type(jtype)
          nloop_type(jtype) = nloop_tmp
          if(pdb) then
            pdbfile(itype) = pdbfile(jtype)
            pdbfile(jtype) = xyzfile
          end if
          linestrut(itype,1) = linestrut(jtype,1)
          linestrut(itype,2) = linestrut(jtype,2)
          linestrut(jtype,1) = linesttmp1
          linestrut(jtype,2) = linesttmp2
        end if
      end if
    end do
  end do

  ! Computing the number of variables
  !
  ! ntype: 1...ntype (counter for the number of free structures)
  !
  ! ntfix: 1...ntype...ntfix (counter for the total number of structures)
  !

  ntmol = 0
  do itype = 1, ntfix
    ntmol = ntmol + nmols(itype)
  end do
  ntotmol = 0
  do itype = 1, ntype
    ntotmol = ntotmol + nmols(itype)
  end do
  n = ntotmol * 6
  write(*,*) ' Total number of molecules: ', ntmol
  write(*,*) ' Number of fixed molecules: ', ntmol - ntotmol
  write(*,*) ' Number of free molecules: ', ntotmol
  write(*,*) ' Number of variables: ', n

  ! Computing the total number of fixed atoms

  natfix = 0
  if(fix) then
    do iftype = ntype + 1, ntfix
      natfix = natfix + natoms(iftype)
    end do
  end if
  write(*,*) ' Total number of fixed atoms: ', natfix

  ! Setting the array that contains the restrictions per atom

  icart = 0
  do itype = 1, ntype
    rests = .false.
    do imol = 1, nmols(itype)
      idatom = idfirst(itype) - 1
      do iatom = 1, natoms(itype)
        icart = icart + 1
        idatom = idatom + 1
        nratom(icart) = 0
        iratcount = 0
        do i = 1, mrperatom
          iratom(icart,i) = 0
        end do
        iline = linestrut(itype,1)
        do while(iline.lt.linestrut(itype,2))
          iline = iline + 1
          if(keyword(iline,1).eq.'atoms') then
            iiatom = -1
            do iat = 2, maxkeywords
              read(keyword(iline,iat),*,iostat=ioerr) iiatom
              if ( ioerr /= 0 ) then
                if ( iiatom == -1 ) then
                  write(*,*) ' ERROR: Could not read atom selection for type: ', itype
                  stop
                else
                  exit
                end if
              end if
              if ( iiatom > natoms(itype) ) then
                write(*,*) ' ERROR: atom selection with index greater than number of '
                write(*,*) '        atoms in structure ', itype
                stop
              end if
              if(iatom.eq.iiatom) exit
            end do
            do while(keyword(iline,1).ne.'end'.and.&
                     keyword(iline,2).ne.'atoms')
              iline = iline + 1
              if(iatom.eq.iiatom) then
                if(keyword(iline,1).eq.'inside'.or.&
                   keyword(iline,1).eq.'outside'.or.&
                   keyword(iline,1).eq.'over'.or.&
                   keyword(iline,1).eq.'below') then
                  nratom(icart) = nratom(icart) + 1
                  iratcount = iratcount + 1
                  do irest = 1, nrest
                    if(irestline(irest).eq.iline) iirest = irest
                  end do
                  iratom(icart,iratcount) = iirest
                end if
              end if
            end do
            iline = iline - 1
          else if(keyword(iline,1).eq.'inside'.or.&
                  keyword(iline,1).eq.'outside'.or.&
                  keyword(iline,1).eq.'over'.or.&
                  keyword(iline,1).eq.'below') then
            nratom(icart) = nratom(icart) + 1
            iratcount = iratcount + 1
            do irest = 1, nrest
              if(irestline(irest).eq.iline) iirest = irest
            end do
            iratom(icart,iratcount) = iirest
          end if
        end do
        if(nratom(icart).gt.0) rests = .true.
      end do
      if(.not.rests) then
        write(*,*) ' ERROR: Some molecule has no geometrical',&
                   ' restriction defined: nothing to do.'
        stop
      end if
    end do
  end do

  ! Read the constraints to rotations about axis, if set

  do itype = 1, ntype
    constrain_rot(itype,1) = .false.
    constrain_rot(itype,2) = .false.
    constrain_rot(itype,3) = .false.
    iline = linestrut(itype,1)
    do while(iline.lt.linestrut(itype,2))
      iline = iline + 1
      if(keyword(iline,1).eq.'constrain_rotation') then
        if(iline.gt.linestrut(itype,1).and.&
           iline.lt.linestrut(itype,2)) then

           ! Note that for movable molecules, teta is a rotation on the x-axis,
           !                                  gama is a rotation on the z-axis,
           !                                  beta is a rotation on the y-axis
           !                                  (see eulerrmat routine)

          if(keyword(iline,2).eq.'x') then
            constrain_rot(itype,3) = .true.
            read(keyword(iline,3),*) rot_bound(itype,3,1)
            read(keyword(iline,4),*) rot_bound(itype,3,2)
            rot_bound(itype,3,1) = rot_bound(itype,3,1)*pi/180.d0
            rot_bound(itype,3,2) = rot_bound(itype,3,2)*pi/180.d0

            write(*,*) ' Rotations about x axis of molecules of ',&
                       ' type ', itype, ' will be constrained. '
          end if
          if(keyword(iline,2).eq.'y') then
            constrain_rot(itype,1) = .true.
            read(keyword(iline,3),*) rot_bound(itype,1,1)
            read(keyword(iline,4),*) rot_bound(itype,1,2)
            rot_bound(itype,1,1) = rot_bound(itype,1,1)*pi/180.d0
            rot_bound(itype,1,2) = rot_bound(itype,1,2)*pi/180.d0

            write(*,*) ' Rotations about y axis of molecules of ',&
                       ' type ', itype, ' will be constrained. '
          end if
          if(keyword(iline,2).eq.'z') then
            constrain_rot(itype,2) = .true.
            read(keyword(iline,3),*) rot_bound(itype,2,1)
            read(keyword(iline,4),*) rot_bound(itype,2,2)
            rot_bound(itype,2,1) = rot_bound(itype,2,1)*pi/180.d0
            rot_bound(itype,2,2) = rot_bound(itype,2,2)*pi/180.d0

            write(*,*) ' Rotations about z axis of molecules of ',&
                       ' type ', itype, ' will be constrained. '
          end if
          if ( keyword(iline,2) /= 'x' .and. &
               keyword(iline,2) /= 'y' .and. &
               keyword(iline,2) /= 'z' ) then
            write(*,*) ' ERROR: constrain_rotation option not properly defined (not x, y, or z) '
            stop
          end if
        end if
      end if
    end do
  end do

  ! Setting the vector that contains the default tolerances

  do i = 1, ntotat
    radius(i) = dism/2.d0
    fscale(i) = 1.d0
    if ( use_short_tol ) then
      use_short_radius(i) = .true.
    else
      use_short_radius(i) = .false.
    end if
    short_radius(i) = short_tol_dist/2.d0
    short_radius_scale(i) = short_tol_scale
  end do

  ! Setting the radius defined for atoms of each molecule,
  ! but not atom-specific, first

  icart = 0
  do itype = 1, ntfix
    iline = linestrut(itype,1)
    iline_atoms = 0
    do while( iline <= linestrut(itype,2) )
      if ( keyword(iline,1) == "atoms" ) then
        iline_atoms = iline
        iline = iline + 1
        cycle
      end if
      if ( keyword(iline,1) == "end" .and. &
           keyword(iline,2) == "atoms" ) then
        iline_atoms = 0
        iline = iline + 1
        cycle
      end if
      if ( iline_atoms == 0 ) then
        !
        ! Read radius
        !
        if ( keyword(iline,1) == "radius" ) then
          read(keyword(iline,2),*,iostat=ioerr) value
          if ( ioerr /= 0 ) then
            write(*,*) ' ERROR: Could not read radius from keyword. '
            stop
          end if
          iicart = icart
          do imol = 1, nmols(itype)
            do iatom = 1, natoms(itype)
              iicart = iicart + 1
              radius(iicart) = value
            end do
          end do
        end if
        !
        ! Read minimum-distance function scale
        !
        if ( keyword(iline,1) == "fscale" ) then
          read(keyword(iline,2),*,iostat=ioerr) value
          if ( ioerr /= 0 ) then
            write(*,*) ' ERROR: Could not read fscale value from keyword. '
            stop
          end if
          iicart = icart
          do imol = 1, nmols(itype)
            do iatom = 1, natoms(itype)
              iicart = iicart + 1
              fscale(iicart) = value
            end do
          end do
        end if
        !
        ! Read short_radius
        !
        if ( keyword(iline,1) == "short_radius" ) then
          read(keyword(iline,2),*,iostat=ioerr) value
          if ( ioerr /= 0 ) then
            write(*,*) ' ERROR: Could not read short_radius value from keyword. '
            stop
          end if
          iicart = icart
          do imol = 1, nmols(itype)
            do iatom = 1, natoms(itype)
              iicart = iicart + 1
              short_radius(iicart) = value
              use_short_radius(iicart) = .true.
            end do
          end do
        end if
        !
        ! Read short_radius scale
        !
        if ( keyword(iline,1) == "short_radius_scale" ) then
          read(keyword(iline,2),*,iostat=ioerr) value
          if ( ioerr /= 0 ) then
            write(*,*) ' ERROR: Could not read short_radius_scale value from keyword. '
            stop
          end if
          iicart = icart
          do imol = 1, nmols(itype)
            do iatom = 1, natoms(itype)
              iicart = iicart + 1
              short_radius_scale(iicart) = value
              use_short_radius(iicart) = .true.
            end do
          end do
        end if
      end if
      iline = iline + 1
    end do
    icart = icart + nmols(itype)*natoms(itype)
  end do

  ! If some radius was defined using atom-specific definitions, overwrite
  ! the general radius defined for the molecule

  icart = 0
  do itype = 1, ntfix
    iline = linestrut(itype,1)
    iline_atoms = 0
    do while( iline <= linestrut(itype,2) )
      if ( keyword(iline,1) == "atoms" ) then
        iline_atoms = iline
        iline = iline + 1
        cycle
      end if
      if ( keyword(iline,1) == "end" .and. &
           keyword(iline,2) == "atoms" ) then
        iline_atoms = 0
        iline = iline + 1
        cycle
      end if
      if ( iline_atoms /= 0 ) then
        !
        ! Read atom specific radius
        !
        if ( keyword(iline,1) == "radius" ) then
          read(keyword(iline,2),*,iostat=ioerr) value
          if ( ioerr /= 0 ) then
            write(*,*) ' ERROR: Could not read radius from keyword. '
            stop
          end if
          ival = 2
          do
            read(keyword(iline_atoms,ival),*,iostat=ioerr) iat
            if ( ioerr /= 0 ) exit
            if ( iat > natoms(itype) ) then
              write(*,*) ' ERROR: atom selection with index greater than number of '
              write(*,*) '        atoms in structure ', itype
              stop
            end if
            radius(icart+iat) = value
            ival = ival + 1
          end do
        end if
        !
        ! Read atom specific function scale
        !
        if ( keyword(iline,1) == "fscale" ) then
          read(keyword(iline,2),*,iostat=ioerr) value
          if ( ioerr /= 0 ) then
            write(*,*) ' ERROR: Could not read fscale value from keyword. '
            stop
          end if
          ival = 2
          do
            read(keyword(iline_atoms,ival),*,iostat=ioerr) iat
            if ( ioerr /= 0 ) exit
            if ( iat > natoms(itype) ) then
              write(*,*) ' ERROR: atom selection with index greater than number of '
              write(*,*) '        atoms in structure ', itype
              stop
            end if
            fscale(icart+iat) = value
            ival = ival + 1
          end do
        end if
        !
        ! Read atom specific short radius
        !
        if ( keyword(iline,1) == "short_radius" ) then
          read(keyword(iline,2),*,iostat=ioerr) value
          if ( ioerr /= 0 ) then
            write(*,*) ' ERROR: Could not read short_radius value from keyword. '
            stop
          end if
          ival = 2
          do
            read(keyword(iline_atoms,ival),*,iostat=ioerr) iat
            if ( ioerr /= 0 ) exit
            if ( iat > natoms(itype) ) then
              write(*,*) ' ERROR: atom selection with index greater than number of '
              write(*,*) '        atoms in structure ', itype
              stop
            end if
            short_radius(icart+iat) = value
            use_short_radius(icart+iat) = .true.
            ival = ival + 1
          end do
        end if
        !
        ! Read atom specific short radius function scale
        !
        if ( keyword(iline,1) == "short_radius_scale" ) then
          read(keyword(iline,2),*,iostat=ioerr) value
          if ( ioerr /= 0 ) then
            write(*,*) ' ERROR: Could not read short_radius_scale value from keyword. '
            stop
          end if
          ival = 2
          do
            read(keyword(iline_atoms,ival),*,iostat=ioerr) iat
            if ( ioerr /= 0 ) exit
            if ( iat > natoms(itype) ) then
              write(*,*) ' ERROR: atom selection with index greater than number of '
              write(*,*) '        atoms in structure ', itype
              stop
            end if
            short_radius_scale(icart+iat) = value
            use_short_radius(icart+iat) = .true.
            ival = ival + 1
          end do
        end if
      end if
      iline = iline + 1
    end do
    iicart = icart
    icart = icart + natoms(itype)
    do imol = 2, nmols(itype)
      do iatom = 1, natoms(itype)
        icart = icart + 1
        radius(icart) = radius(iicart+iatom)
        fscale(icart) = fscale(iicart+iatom)
        short_radius(icart) = short_radius(iicart+iatom)
        short_radius_scale(icart) = short_radius_scale(iicart+iatom)
        use_short_radius(icart) = use_short_radius(iicart+iatom)
      end do
    end do
  end do

  ! Check if the short radii were set correctly, if the case

  ioerr = 0
  do i = 1, ntotat
    if ( use_short_radius(i) ) then
     if ( short_radius(i) >= radius(i) ) then
       write(*,*) ' ERROR: The short radius must be smaller than the default radius. '
       write(*,*) '        (the default radius is one half of the default tolerance).'
       stop
     end if
   end if
  end do

  ! If there are no variables (only fixed molecules, stop)

  if(n.eq.0) then
    call output(n,x)
    write(*,dash1_line)
    write(*,*) ' There are only fixed molecules, therefore there is nothing to do. '
    write(*,*) ' The output file contains the fixed molecules in the desired positions. '
    write(*,dash1_line)
    write(*,*) ' Wrote output file: ', trim(adjustl(xyzout))
    if ( crd ) write(*,*) ' ... and to CRD file: ', trim(adjustl(crdfile))
    write(*,dash1_line)
    stop
  end if

  !
  ! (Re)setting parameters and building initial point
  !

  call initial(n,x)

  ! Computing the energy at the initial point

  radscale = 1.d0
  do i = 1, ntotat
    radius_ini(i) = radius(i)
  end do
  call computef(n,x,all_type_fx)
  write(*,*) ' Objective function at initial point: ', all_type_fx
  fprint = all_type_fx

  ! Stop if only checking the initial approximation

  if(check) then
    call output(n,x)
    write(*,*) ' Wrote initial point to output file: ', xyzout(1:strlength(xyzout))
    if ( crd ) write(*,*) ' ... and to CRD file: ', trim(adjustl(crdfile))
    stop
  end if

  !
  ! Main loop: first pack types of molecules separately, then
  ! pack all molecules together
  !

  call swaptype(n,x,itype,0) ! Save all-molecule vector data
  itype = 0
  main : do while(itype <= ntype)
    itype = itype + 1
    if ( packall ) itype = ntype + 1

    ! Use larger tolerance than required to improve separation

    radscale = discale
    do i = 1, ntotat
      radius(i) = discale*radius_ini(i)
    end do

    ! Set vectors for specific or all-molecule packing

    if ( itype <= ntype ) then
      call swaptype(n,x,itype,1) ! Set vectors to pack only this type of molecule
    else
      call swaptype(n,x,itype,3) ! Restore all-molecule vectors
    end if

    ! Print titles

    write(*,hash3_line)
    if ( itype <= ntype ) then
      write(*,*) ' Packing molecules of type: ', input_itype(itype)
    else
      write(*,*) ' Packing all molecules together '
    end if
    write(*,hash3_line)

    ! Checking if first approximation is a solution

    call computef(n,x,fx)

    if ( fdist < precision .and. frest < precision ) then

      write(*,*)
      write(*,*) ' Initial approximation is a solution. Nothing to do. '
      write(*,*)
      call swaptype(n,x,itype,3) ! Restore all-molecule vectors
      call output(n,x)
      if( itype == ntype + 1 ) then
        write(*,*) ' Solution written to file: ', trim(adjustl(xyzout))
        if ( crd ) write(*,*) ' ... and to CRD file: ', trim(adjustl(crdfile))
      else
        write(*,*) ' Current point written to file: ', trim(adjustl(xyzout))
        if ( crd ) write(*,*) ' ... and to CRD file: ', trim(adjustl(crdfile))
      end if
      call writesuccess(itype,fdist,frest,fx)

    ! Otherwise, pack the molecules

    else

      loop = -1
      gencanloop : do while(loop.lt.nloop)
        loop = loop + 1

        ! Reseting the parameters relative to the improvement of the function

        if(loop.eq.0) then
          fimp = 1.d99
          fimprov = fimp
          do i = 1, ntotat
            radiuswork(i) = radius(i)
            radius(i) = radius_ini(i)
          end do
          call computef(n,x,fx)
          do i = 1, ntotat
            radius(i) = radiuswork(i)
          end do
          bestf = fx
          flast = fx
        end if

        ! Moving bad molecules

        if(radscale == 1.d0 .and. fimp.le.10.d0) then
          movebadprint = .true.
          call movebad(n,x,fx,movebadprint)
          flast = fx
        end if

        write(*,dash3_line)
        write(*,*) ' Starting GENCAN loop: ', loop
        write(*,*) ' Scaling radii by: ', radscale
        write(*,*)

        ! CALL GENCAN

        write(*,prog1_line)
        call pgencan(n,x,fx)

        !
        ! Compute the statistics of the last optimization loop
        !

        ! Use the user-specified radii for statistics

        do i = 1, ntotat
          radiuswork(i) = radius(i)
          radius(i) = radius_ini(i)
        end do
        call computef(n,x,fx)

        if(bestf.gt.0.d0) fimprov = -100.d0 * (fx - bestf) / bestf
        if(bestf.eq.0.d0) fimprov = 100.d0
        if(flast.gt.0.d0) fimp = -100.d0 * (fx - flast) / flast
        if(flast.eq.0.d0) fimp = 100.d0
        fimp = dmin1(99.99d0,dmax1(-99.99d0,fimp))
        fimprov = dmin1(99.99d0,dmax1(-99.99d0,fimprov))

        write(*,"(/&
                  &'  Function value from last GENCAN loop: f = ', e10.5, /&
                  &'  Best function value before: f = ', e10.5,           /&
                  &'  Improvement from best function value: ', f8.2, ' %',/&
                  &'  Improvement from last loop: ', f8.2, ' %',          /&
                  &'  Maximum violation of target distance: ', f12.6,     /&
                  &'  Maximum violation of the constraints: ', e10.5       &
                  &)") fx, bestf, fimprov, fimp, fdist, frest
        flast = fx

        !
        ! Analysis of final loop packing and output data
        !

        if ( itype <= ntype ) then

          ! Save best function value for this packing

          if ( fx < bestf ) bestf = fx

          ! Check if this point is a solution

          call swaptype(n,x,itype,2) ! Save this type current point
          ! If the solution was found for this type
          if( fdist < precision .and. frest < precision ) then
            call swaptype(n,x,itype,3) ! Restore all molecule vectors
            call output(n,x)
            write(*,*) ' Current structure written to file: ', trim(adjustl(xyzout))
            if ( crd ) write(*,*) ' ... and to CRD file: ', trim(adjustl(crdfile))
            call writesuccess(itype,fdist,frest,fx)
            exit gencanloop
          end if

          ! Compute and report function value for all-type packing

          call swaptype(n,x,itype,3) ! Restore all molecule vectors
          call computef(n,x,all_type_fx)
          write(*,"('  All-type function value: ', e10.5 )") all_type_fx

        else

          call computef(n,x,fx)
          all_type_fx = fx
          if ( fx < bestf ) bestf = fx
          ! If solution was found for all system
          if ( fdist < precision .and. frest < precision ) then
            call output(n,x)
            call writesuccess(itype,fdist,frest,fx)
            write(*,*) ' Solution written to file: ', xyzout(1:strlength(xyzout))
            if ( crd ) write(*,*) ' ... and to CRD file: ', trim(adjustl(crdfile))
            write(*,dash3_line)
            exit main
          end if

        end if
        write(*,dash3_line)

        ! If this is the best structure so far
        if( mod(loop+1,writeout) == 0 .and. all_type_fx < fprint ) then
          call output(n,x)
          write(*,*) ' Current solution written to file: ', trim(adjustl(xyzout))
          if ( crd ) write(*,*) ' ... and to CRD file: ', trim(adjustl(crdfile))
          fprint = all_type_fx
          do i = 1, n
            xprint(i) = x(i)
          end do

        ! If the user required printing even bad structures
        else if ( mod(loop+1,writeout) == 0 .and. writebad ) then
          call output(n,x)
          write(*,*) ' Writing current (perhaps bad) structure to file: ', trim(adjustl(xyzout))
          if ( crd ) write(*,*) ' ... and to CRD file: ', trim(adjustl(crdfile))
        end if

        ! Restore vector for packing this type of molecule, if the case

        if ( itype <= ntype ) then
          call swaptype(n,x,itype,0) ! Reset type vectors
          call swaptype(n,x,itype,1) ! Set vector for molecules of this type
          call computef(n,x,fx)
        end if

        ! Restore the working radii

        do i = 1, ntotat
          radius(i) = radiuswork(i)
        end do
        if ( radscale > 1.d0 ) then
          if( ( fdist < precision .and. fimp < 10.d0 ) .or. &
                fimp < 2.d0 ) then
            radscale = dmax1(0.9*radscale,1.d0)
            do i = 1, ntotat
              radius(i) = dmax1(radius_ini(i),0.9d0*radius(i))
            end do
          end if
        end if

        if(loop.eq.nloop) then
          if ( itype .eq. ntype+1 ) then
            write(*,*)' STOP: Maximum number of GENCAN loops achieved.'
            call checkpoint(n,xprint)
            exit main
          else
            write(*,*)' Maximum number of GENCAN loops achieved.'
          end if
        end if

      end do gencanloop

    end if

  end do main

  write(*,*) '  Running time: ', etime(tarray) - time0,' seconds. '
  write(*,dash3_line)
  write(*,*)

end program packmol

max_stars_repo_name: SSchott/packmol

max_issues_repo_name: SSchott/packmol

max_forks_repo_name: SSchott/packmol


text:
Address(Condor Court) is a residential street in North Davis.

 Intersecting Streets

  Eagle Place
  Valencia Avenue


text:
subroutine sub_one(ix, oy, ioz)
   real, intent(in) :: ix ! input
   real, intent(out)  :: oy ! output
   real, intent(inout) :: ioz ! input and output

   oy = sqrt(real(ix)) + ioz
   ioz = max(ix, oy)

end subroutine

function func(i) result(j)
   integer, intent(in) :: i ! input
   integer              :: j ! output

   j = sqrt(real(i)) + log(real(i))
end function

program main

   implicit none

   integer :: i
   integer :: func
   real :: ix, oy, ioz

   ix = 10.0
   ioz = 20.0

   i = huge(1)
   print *, "i=", i, char(10), " sqrt(i) + log(i) =", func(i)

   print *, 'Before:', ix, oy, ioz
   call sub_one(ix, oy, ioz)
   print *, 'After :', ix, oy, ioz

end program

max_stars_repo_name: WVUHPC/Modern-Fortran

max_issues_repo_name: WVUHPC/Modern-Fortran

max_forks_repo_name: WVUHPC/Modern-Fortran


text:
	SUBROUTINE MT_ALTM ( stralt, idecd, iret )
C************************************************************************
C* MT_ALTM                                                              *
C*                                                                      *
C* This subroutine decodes the altimeter group in a METAR report.	*
C* If the group exists, the altimeter value is stored in common using	*
C* the units in which it was reported.                                  *
C* 								        *
C* MT_ALTM ( STRALT, IDECD, IRET )                                      *
C*								        *
C* Input parameters: 						        *
C*      STRALT		CHAR*		Possible altimeter field        *
C*								        *
C* Output parameters:						        *
C*	RIVALS(IRALTI)	REAL		Altimeter setting in inches     *
C*	RIVALS(IRALTM)	REAL		Altimeter setting in hPa (mb)   *
C*	IDECD		INTEGER		Decode decision flag            *
C*					     1 = field decoded       	*
C*					     0 = field not decoded    	*
C*	IRET		INTEGER		Return code                     *
C*	      		    		     0 = normal return		*
C*	                                    13 = miscoded field		*
C*	                                    -1 = field not found       	*
C*					                                *
C**								        *
C* Log:							        	*
C* D. Kidwell/NCEP 	 4/95                                           *
C* D. Kidwell/NCEP 	 2/96   Changed altim upper bound to 32.48      *
C* K. Tyle/GSC		 1/97	Reorganized header and comments; 	*
C*				changed call to DC_WLOG			*
C* K. Tyle/GSC		 2/97	Changed error processing		*
C* D. Kidwell/NCEP 	 6/97   ST_LSTR -> INDEX and ST_CRNM -> ST_INTG *
C* D. Kidwell/NCEP 	 4/98   New interface                           *
C************************************************************************
	INCLUDE 	'GEMPRM.PRM'
	INCLUDE 	'mtcmn.cmn'
C*
	CHARACTER*(*) 	stralt
C*
	LOGICAL 	alterr
C------------------------------------------------------------------------
	alterr = .false.
	iret   = 0
	idecd  = 0
C
C*	Look for group starting with 'A' (inches) or 'Q' (hectopascals)
C
	IF ( ( stralt ( 1:1 ) .eq. 'A' ) .or.
     +       ( stralt ( 1:1 ) .eq. 'Q' ) ) THEN
	    lens = INDEX ( stralt, ' ' ) - 1
	    IF ( lens .ge. 5 ) THEN
	        CALL ST_INTG ( stralt ( 2:5 ), ialtim, iret )
	        IF ( iret .eq. 0 ) THEN
		    idecd = 1
		    IF ( stralt ( 1:1 ) .eq. 'A' ) THEN
C
C*		        Check bounds, then convert to inches if good.
C
			IF ( ialtim .lt. 1000 .or. ialtim .gt. 3350 )
     +			     THEN
			    alterr = .true.
			  ELSE
			    rivals ( iralti ) = FLOAT ( ialtim ) *.01
			END IF
		      ELSE
C
C*			Check bounds in hPa.
C
			IF ( ialtim .lt.  400 .or. ialtim .gt. 1100 )
     +			     THEN
			    alterr = .true.
			  ELSE
			    rivals ( iraltm ) = FLOAT ( ialtim )
			END IF
		    END IF
		  ELSE
		    alterr = .true.
		END IF
	      ELSE
		alterr = .true.
	    END IF
	  ELSE
C
C*	    No 'A' or 'Q' group was found.
C
	    iret = -1
	END IF
C
	IF ( alterr ) THEN
	    iret = 13
	    idecd = 1
	END IF
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


! { dg-do run }
! PR77707 formatted direct access: nextrec off by one
program directaccess_formatted
  integer nextrec
  open(10, status='scratch', form='formatted', access='direct', recl=10*4)
  write(10,'(10i4)',rec=9) 1,2,3,4,5,6,7,8,9,10
  inquire(unit=10,nextrec=nextrec)
  if (nextrec.ne.10) call abort
  close(10)
end


! { dg-do compile }
! PR fortran/54223
subroutine r(d)
    implicit none
    integer, optional :: d
    integer :: h, q
    q(d) = d + 1     ! statement function statement
    h = q(d)
end subroutine r

subroutine s(x)
    implicit none
    integer, optional :: x
    integer :: g, z
    g(x) = x + 1     ! statement function statement
    z = g()          ! { dg-error "Missing actual argument" }
end subroutine s

subroutine t(a)
    implicit none
    integer :: a
    integer :: f, y
    f(a) = a + 1     ! statement function statement
    y = f()          ! { dg-error "Missing actual argument" }
end subroutine t
! { dg-prune-output " Obsolescent feature" }


text:
subroutine nrmaxl(rfin,maxk,rinit,kopt,v_dic)
  implicit none
  logical            :: kNN
  integer,intent(in) :: maxk
  real*8,intent(in)  :: rinit
  integer,intent(in) :: kopt
  real*8,intent(in)  :: v_dic(maxk)
  !real*8,intent(out)  :: vi(maxk)
  real*8,intent(out)  :: rfin
  real*8 :: b,L0,a,stepmax,ga,gb,Cov2(2,2),Covinv2(2,2)
  real*8 :: jf,t,s,tt,func,sigma,sa,sb,vi(maxk)
  integer :: j,niter
  !write (6,*) maxk,rinit,kopt,(v_dic(j),j=1,kopt)
  !Integer, Parameter :: kdp = selected_real_kind(15)
  !real (kind=kdp) :: dimset                    ! dimension in which is embedded the dataset (real)
  !real (kind=kdp), Parameter :: nearzero = tiny(dimset)  ! Minimum positive value
  b=rinit
  L0=0.
  a=0.
  stepmax=0.1*abs(b)
  gb=float(kopt)
  ga=float(kopt+1)*float(kopt)/2.
  Cov2(1,1)=0. !gbb
  Cov2(1,2)=0. !gab
  Cov2(2,2)=0. !gaa
  kNN=.false.
  vi(1)=v_dic(1)
  do j=2,kopt
    vi(j)=v_dic(j)-v_dic(j-1)
    if ((vi(j)<1D-300).and.(.not.kNN)) then
       !write (6,*) rinit,kopt,kNN,vi(j),v_dic(j),v_dic(j-1),j
       kNN=.true.
    endif
  enddo
  if (.not.kNN) then
    do j=1,kopt
      jf=float(j)
      t=b+a*jf
      s=exp(t)
      tt=vi(j)*s
      L0=L0+t-tt
      gb=gb-tt
      ga=ga-jf*tt
      Cov2(1,1)=Cov2(1,1)-tt
      Cov2(1,2)=Cov2(1,2)-jf*tt
      Cov2(2,2)=Cov2(2,2)-jf*jf*tt
    enddo
    Cov2(2,1)=Cov2(1,2)
    Covinv2=matinv2(Cov2)
    func=100.
    niter=0
    do while ((func>1D-3).and.(niter.lt.1000))
      sb=(Covinv2(1,1)*gb+Covinv2(1,2)*ga)
      sa=(Covinv2(2,1)*gb+Covinv2(2,2)*ga)
      niter=niter+1
      sigma=0.1
      if (abs(sigma*sb).gt.stepmax) then
        sigma=abs(stepmax/sb)
      endif
      b=b-sigma*sb
      a=a-sigma*sa
      L0=0.
      gb=float(kopt)
      ga=float(kopt+1)*float(kopt)/2.
      Cov2(1,1)=0. !gbb
      Cov2(1,2)=0. !gab
      Cov2(2,2)=0. !gaa
      do j=1,kopt
        jf=float(j)
        t=b+a*jf
        s=exp(t)
        tt=vi(j)*s
        L0=L0+t-tt
        gb=gb-tt
        ga=ga-jf*tt
        Cov2(1,1)=Cov2(1,1)-tt
        Cov2(1,2)=Cov2(1,2)-jf*tt
        Cov2(2,2)=Cov2(2,2)-jf*jf*tt
      enddo
      Cov2(2,1)=Cov2(1,2)
      Covinv2=matinv2(Cov2)
      if ((abs(a).le.tiny(a)).or.(abs(b).le.tiny(b))) then
        func=max(gb,ga)
      else
        func=max(abs(gb/b),abs(ga/a))
      endif
    enddo
  endif
  Cov2(:,:)=-Cov2(:,:)
  Covinv2=matinv2(Cov2)
  rfin=b
contains
  pure function matinv2(A) result(B)
!! Performs a direct calculation of the inverse of a 2×2 matrix.
    real*8, intent(in) :: A(2,2)   !! Matrix
    real*8             :: B(2,2)   !! Inverse matrix
    real*8             :: detinv
    detinv = 1./(A(1,1)*A(2,2) - A(1,2)*A(2,1))
    B(1,1) = +detinv * A(2,2)
    B(2,1) = -detinv * A(2,1)
    B(1,2) = -detinv * A(1,2)
    B(2,2) = +detinv * A(1,1)
  end function
end subroutine NRmaxL

max_stars_repo_name: giovannidoni/DPA

max_issues_repo_name: giovannidoni/DPA

max_forks_repo_name: giovannidoni/DPA


text:
C    %W% %G%
      subroutine int3fh (a)
C
C  THIS SUBROUTINE CALCULATES INITIAL STATE VECTORS AND PAST VALUES
C  FOR EXCITER TYPE FH
C
      include 'tspinc/param.inc'
      include 'tspinc/prt.inc'
      include 'tspinc/params.inc'
      include 'tspinc/gentbla.inc'
      include 'tspinc/pointr.inc'
      include 'tspinc/gentblb.inc'
      include 'tspinc/int3.inc'
      dimension array(40),a(40)
       equivalence         (array(1),  vamax), (array(2),  vamin),
     1 (array(3),   vref), (array(4),  esatx), (array(5),  csatx),
     2 (array(6),    hbr), (array(7),    cke), (array(8),    cka),
     3 (array(9),     dr), (array(10),   cfd), (array(11),   ckf),
     4 (array(12),    df), (array(13),   vlv), (array(14),   ckn),
     5 (array(15),   ckr), (array(16),   ckc), (array(17),    db),
     6 (array(18),    dc), (array(19),   efdn), (array(20),    de),
     7 (array(21),  rcex), (array(22),    dk), (array(23),  xcex),
     8 (array(24),    vaio), (array(25),    da), (array(26),   hba),
     9 (array(27),   vco), (array(28),     hbf), (array(29),   efd),
     * (array(30),   ckd), (array(31),  dk1), (array(32),efdmax),
     * (array(33),   vepo), (array(34),   hcb), (array(35),   hbe),
     * (array(36),    dep), (array(37),bprim), (array(38),   veo),
     * (array(39),  ckap), (array(40),  ckep)

C
C     TRANSFER CITER DATA TO ARRAY
C
      l6 = 6
      do 100 itr = 1,40
100   array(itr) = a(itr)
          efd = efdo
C
C         CALCULATE VOLTAGE TRANSDUCER CIRCUIT
C
          if(xt .eq. 0.0)  then
             vtr = vhd
             vti = vhq
          else
             vtr = vhd + rt*agend - xt*agenq
             vti = vhq + rt*agenq + xt*agend
          endif
          crt = agend
          cit = agenq
          zcr = rcex*crt - xcex*cit
          zci = rcex*cit + xcex*crt
          vcr = zcr + vtr
          vci = zci + vti
          vc = sqrt(vcr**2 + vci**2)
          vco = vc
          vc1 = vc
C
C         CALCULATE R-BLOCK
C
          dr = 2.0*dr/dt + 1.0
          hbr = (dr - 2.0)*vc1 + vc
C
C         CALCULATE FEX-BLOCK (INPUT: EFDO,CFD,CKC;
C                              OUTPUT: FEX,VE)
C
C         Note: Initial field current CFD equals initial field EFDO.
C
          cfd = efdo
C         CFDO = EFDO                           ! removed extraneous var
          call getfex (fex,efdo,ve,cfd,ckc)
C
C         CALCULATE E BLOCK
C
          ckep = csatx * exp(esatx*abs(ve))
          ck1 = ckep + ckep*ve*esatx
          ck2 = ve**2 * esatx*ckep
          vfe = (cke+ck1)*ve + ckd*cfd - ck2
          vepo = 0.0
          de = 2.*de/dt
          dep = de+ck1 +cke
          hbe = ve*(dep-(ck1+cke)) +ck2
          veo = ve
C
C         CALCULATE A-BLOCK (REGULATOR)
C
          vr = vfe
          va = vr/(efdo*ckr)
          x3 = va
          vai = va/cka
          da = 2.0*da/dt + 1.0
          hba = (da - 2.0)*va + cka*vai
C
C         CHECK REGULATOR LIMITS
C
          if (va .ge. vamax .or. va .le. vamin) then
              write (errbuf, 2036) name, base, id, vr
              call prterr ('E',1)
 2036         format('0 INITIAL VA VIOLATES LIMITS FOR ', a8,
     1           2x, f5.1, 2x, a1, 5x, 'INITIAL VA = ', f6.3)
              iabort = 1
              go to 2300
          endif
C
C         CALCULATE F-BLOCK
C
          df = 2.0*df/dt + 1.0
          if(efdo .lt. efdn)then
          fconst = 0.0
          gconst = ckf
          else
          fconst = efdn*(ckf-ckn)
          gconst = ckn
          endif
          dk = 2.0*gconst/dt
          dk1 = 2.*fconst/dt
          hbf = -efdo*dk - dk1
C
C         CALCULATE B-C BLOCK
C
          vaii = vai
          db = 2.0*db/dt + 1.0
          dc = 2.0*dc/dt + 1.0
          hcb = (db - 2.0)*vai - (dc - 2.0)*vaii
          verr = vai
          vref = verr + vc1
          ckep = 1.0
          ckap = cka
C
C   CALCULATE APPROXIMATION EFD = ECONST*VE -DCONST
C
          cin = abs(ckc*cfd/ve)
          if(cin .gt. 0.715)then
          econst = 1.68
          dconst = 1.714*ckc*cfd
          else
          if(cin .lt. 0.51)then
          econst = 1.0
          dconst = 0.58*ckc*cfd
          else
          econst = 0.93227 + 0.865*ckc*ckc*cfd*cfd/(veo*veo)
          dconst = -1.730*ckc*ckc*cfd/veo -0.01429*ckc
          endif
          endif
C
C      CALCULATE REDUCED EQUATION COEFFICIENTS
C
          aconst =-(dk*dc*cka*ckr)/(df*dep*da*db)
          bprim = aconst*hbf/dk + (hcb*cka*ckr)/(dep*da*db)
     1 + (hba*ckr)/(dep*da)
          bconst = bprim -1.0/econst- aconst*df*verr/dk
          cconst = hbe/dep - dconst/econst - ckd*cfd/dep
          efdtst = (-bconst - sqrt(bconst*bconst - 4.0*aconst*cconst))
     1 /(2.0*aconst)
C
C         TRANSFER DATA TO CITER AND RETURN TO INITAL3
C
 2300     do 2400 itr = 1,40
 2400     a(itr) = array(itr)
          return
          end

max_stars_repo_name: mbheinen/bpa-ipf-tsp

max_issues_repo_name: cuihantao/bpa-ipf-tsp

max_forks_repo_name: mbheinen/bpa-ipf-tsp