text:
 implicit none

 integer :: int

 print*,'Give an integer between 1 and 99'; read*,int
 if (int<1.or.int>99) then
   print*,'Read the instructions more carefully next time! Good bye.'
 else if (int==8.or.int==88) then
   print*,'A lucky number; Congratulations!'
 else if (int==13.or.int==4) then
   print*,'Bad luck...not a good number; beware!'
 else
   print*,'Nothing special with this number, '
   if (mod(int,2)==0) then
     print*,'but I can tell you that it is an even number.'
   else
     print*,'but I can tell you that it is an odd number.'
   end if
 end if

 end

max_stars_repo_name: agforero/FTFramework

max_issues_repo_name: agforero/fortran-testing-framework

max_forks_repo_name: agforero/fortran-testing-framework


! How to remove ^M in source code on linux %s/^M//g (^M = Ctrl+v+m)
! Note on subroutine random_number
! Returns a single pseudorandom number or an array of pseudorandom numbers from the uniform distribution over the range [0, 1).
! The runtime-library implements the xorshift1024* random number generator(RNG). This generator has a period of 2^{1024} - 1,
! and when using multiple threads up to 2^{512} threads can each generate 2^{512} random numbers before any aliasing occurs.

module parameters
    integer, parameter :: latt_leng = 80
    integer, parameter :: EXCH_CONST = 1 ! exchange constant
    integer, dimension(-1 : latt_leng, -1 : latt_leng) :: spin_of ! (-1 ... latt_leng) X (-1 ... latt_leng)
    double precision, dimension(-2 : 2) :: accept_prob_chart
end module

program Ising_Model
    use parameters
    implicit none ! Require all variables to be explicitly declared
    integer(kind=8), parameter :: MC_max = 1E7 * latt_leng**2
    double precision, parameter :: T = 2.269d0 ! Unit = kT/J, Tc computed with each of the physical observables

    call create_lattice_random
    if(EXCH_CONST == -1) then ! Anti-ferromagnetic, EXCH_CONST = -1
        accept_prob_chart = (/ exp(-8d0 / T), exp(-4d0 / T), &
        1d0, exp(4d0 / T), exp(8d0 / T) /) ! Acceptance probability chart
    else ! Ferromagnetic, EXCH_CONST = 1
        accept_prob_chart = (/ exp(8d0 / T), exp(4d0 / T), &
        1d0, exp(-4d0 / T), exp(-8d0 / T) /)
    end if
    call flip_sites_metropolis(MC_max) ! Iterate MC simulation MC_max times per temperature
end program

subroutine create_lattice_random
    use parameters
    integer :: x, y
    double precision :: u

    do x = 0, latt_leng - 1
        do y = 0, latt_leng - 1
            call random_number(u) ! Return random real number in [0,1)
            if (u >= 0.5) then
                spin_of(x, y) = 1 ! Assign random spin values
            else
                spin_of(x, y) = -1
            end if
        end do
    end do

    do x = 0, latt_leng - 1 ! Assign spins to top and bottom line of padding sites to make it periodic
        spin_of(x, -1) = spin_of(x, latt_leng - 1)
        spin_of(x, latt_leng) = spin_of(x, 0)
    end do
    do y = 0, latt_leng -1  ! Assign spins to left and right line of padding sites to make it periodic
        spin_of(-1, y) = spin_of(latt_leng - 1, y)
        spin_of(latt_leng, y) = spin_of(0, y)
    end do
end subroutine

subroutine create_lattice_aligned
    use parameters
    integer :: x, y

    do x = 0, latt_leng - 1 ! Assign spins on each site
        do y = 0, latt_leng - 1
            spin_of(x, y) = 1 ! All up spins
        end do
    end do
    do x = 0, latt_leng - 1 ! Assign spins to top and bottom line of padding sites to make it periodic
        spin_of(x, -1) = spin_of(x, latt_leng - 1)
        spin_of(x, latt_leng) = spin_of(x, 0)
    end do
    do y = 0, latt_leng -1  ! Assign spins to left and right line of padding sites to make it periodic
        spin_of(-1, y) = spin_of(latt_leng - 1, y)
        spin_of(latt_leng, y) = spin_of(0, y)
    end do
end subroutine

subroutine flip_sites_metropolis(MC_max)
    use parameters
    integer :: first_open = 1
    integer :: x, y, sum_spin_nbr, sum_nbr_spin
    integer :: out_unit1 = 10, out_unit2 = 20 !,out_unit3 = 30
    integer(kind=8), intent(in) :: MC_max
    integer(kind=8) :: MC_step
    double precision :: accept_prob, tot_en, tot_mag, tot_en2, tot_E, tot_M, tot_E2
    double precision :: u, v, r

    n_sample = 0 ! Initialize variable to be used in averaging samples
    sum_tot_E = 0d0
    sum_tot_M = 0d0
    sum_tot_E2 = 0d0
    call compute_tot_E_M_E2(tot_en, tot_mag, tot_en2) ! Compute initial values of physical observables
    tot_E = tot_en
    tot_M = tot_mag
    tot_E2 = tot_en2
    do MC_step = 0, MC_max
        call random_number(u) ! Return random real number in [0,1)
        call random_number(v)
        x = FLOOR(u * latt_leng) ! Pick random site, [0, latt_leng - 1]
        y = FLOOR(v * latt_leng)
        sum_nbr_spin = spin_of(x - 1, y) + spin_of(x + 1, y) + spin_of(x, y + 1) + spin_of(x, y - 1)
        accept_prob = accept_prob_chart((spin_of(x, y) * sum_nbr_spin)/2) ! Sigma(SS) = spin_of(x, y) * sum_surrounding = -4, -2 , 0, 2, 4
        call random_number(r) ![0, 1)
        if (r < accept_prob) then
            spin_of(x, y) = - spin_of(x, y) ! If random number < accept_prob => flip site

            if (x == 0) then ! Update padding sites
                spin_of(latt_leng, y) = spin_of(x, y)
            else if (x == latt_leng - 1) then
                spin_of(-1, y) = spin_of(x, y)
            end if
            if (y == 0) then
                spin_of(x, latt_leng) = spin_of(x, y)
            else if (y == latt_leng - 1) then
                spin_of(x, -1) = spin_of(x, y)
            end if

            tot_M = tot_M + 2*spin_of(x, y) ! Update physical observables
            sum_spin_nbr = spin_of(x + 1, y) + spin_of(x - 1, y) + spin_of(x, y + 1) + spin_of(x, y - 1)
            tot_E = tot_E + (- EXCH_CONST) * spin_of(x, y) * sum_spin_nbr
            tot_E2 = tot_E**2
        end if

        if (MOD(MC_step, latt_leng**2) == 0) then
            if(first_open == 1) then
                open (out_unit1, FILE = 'output_ising_model_fixed_T_E_change_plot.txt', action = "write", status = 'replace')
                write (out_unit1,*) MC_step / latt_leng**2 , " ", tot_E
                close(out_unit1)
                open (out_unit2, FILE = 'output_ising_model_fixed_T_M_change_plot.txt', action = "write", status = 'replace')
                write (out_unit2,*) MC_step / latt_leng**2, " ", tot_M
                close(out_unit2)
                first_open = -1
            else
                open (out_unit1, FILE = 'output_ising_model_fixed_T_E_change_plot.txt', action = "write", position = 'append')
                write (out_unit1,*) MC_step / latt_leng**2, " ", tot_E
                close(out_unit1)
                open (out_unit2, FILE = 'output_ising_model_fixed_T_M_change_plot.txt', action = "write", position = 'append')
                write (out_unit2,*) MC_step / latt_leng**2, " ", tot_M
                close(out_unit2)
            end if
        end if
    end do
end subroutine

subroutine compute_tot_E_M_E2(tot_E, tot_M, tot_E2)
    use parameters
    integer :: x, y, sum_spin_nbr
    double precision, intent(out) :: tot_E, tot_M, tot_E2

    tot_E = 0d0
    tot_M = 0d0
    tot_E2 = 0d0
    do x = 0, latt_leng - 1
        do y = 0, latt_leng - 1
            tot_M = tot_M + spin_of(x, y)
            sum_spin_nbr = spin_of(x + 1, y) + spin_of(x - 1, y) + spin_of(x, y + 1) + spin_of(x, y - 1)
            tot_E = tot_E + (- EXCH_CONST) * spin_of(x, y) * sum_spin_nbr
        end do
    end do
    tot_E = tot_E / 2
    tot_E2 = tot_E**2
end subroutine


program teams_new_index
  !! author: Nathan Weeks
  !!
  !! Test FORM TEAM NEW_INDEX= specifier
  use iso_fortran_env, only : team_type
  use oc_assertions_interface, only : assert
  implicit none

  integer :: new_index
  type(team_type) :: team

  new_index = merge(num_images(), this_image()-1, this_image() == 1)

  form team (1, team, new_index=new_index)
  change team(team)
    call assert( this_image() == new_index, "this_image() /= new_index" )
  end team

  sync all

  if (this_image() == 1) print *,"Test passed."

end program


text:
The 4way intersection of Olive Drive and Richards Boulevard is widely considered the most inconvenient, inefficient, and generally heinous intersection in Davis.
  So is the worst intersection only on the olive drive side, or does it include the 1st and E part as well?  Users/ArlenAbraham arlen
  I think it should include both, and maybe even the OnRamp of Doom on Interstate 80 at Richards. Users/BrentLaabs The Labster

To be fair to the planners, there are a number of reasons to leave it the way it is.  First, anybody who is holding onto the concept of Davis as a small town forever, it provides a way to limit traffic into downtown.  That also means that a change to this intersection will affect downtowns traffic patterns, requiring significant changes and likely increasing pedestrian and biker accidents.

One of the most compelling hidden aspects that keep the intersection from being expanded is the fact that the Richards underpass is the last remaining rail span bridge (look up when you go through it).  Once, hundreds of these bridges gave access to the western half of the country and were common sights.  Now this last one remains, having originally allowed the Lincoln Highway to pass underneath the railway.

Also see The Worst Intersection on Campus.

 Description

 Overhead View

 Size
One lane each way, with lines sometimes stretching several blocks in each direction

 Features
   Interstate 80 dumps out onto Richards one block away
   South Davis residents must all funnel through to reach UCD or downtown
   InNOut, Redrum Burger and gas station attract more traffic
   Bridge for train makes expansion of Richards Blvd. Richards unfeasible  I heard it had something to do with the Boy Scout Cabin
   Incline of overpass makes bikes an unattractive alternative
   Close, skewed intersection of Richards Blvd. Richards with First and E streets disorients drivers
   Bizarre merging of Interstate 80 I80 offramp with eastbound Richards
   No leftturn signal for drivers on Olive
   New apartment complexes continually added along Olive
   Offramp from Interstate 80 I80 to far end of Olive, but no corresponding onramp
   No other outlet from Olive
   Traffic lights are set to a timer, not a sensor, guaranteeing a wait at any hour
   With the frequent pedestrians, timer allows roughly 4 cars on Olive to turn left
   Unitrans W Bus Line stops right at the intersection

 Events
   Gridlock
   Sometimes the train goes by overhead
   In evening, headlights gleam like strings of Christmas lights
   In daytime, glare from windshields resembles tinsel from a distance
   If you are turning left from Olive during the day, you will wait roughly 12 minutes
   If a vehicle breaks down, or a truck does not have enough clearance to make it under the train bridge, you will be stuck for hours
   The Davis Police Department occasionally has http://www.livejournal.com/community/ucdavis/2111235.html DUI checkpoints here, stopping every car to see if the people are drunk.
  Prone to flooding. When I was a student in Davis there were at least two times when the tunnel under the railroad tracks flooded, causing traffic to backup and only the people with tall vehicles (trucks, etc.) could go through. I tried once...in my Ford Escort (college car) came up midway to my door and I was lucky it didnt completely kill my car!

Truly a gem of urban planning.

 Tips to avoid this intersection

 By Bike
  Take the South Davis Bike Path, see map on Bike Paths this page.

 By Car
  Use Old Davis Road and the UC Davis I80 underpass.
  Use 2nd Street and the Mace Boulevard Mace Blvd I80 overpass.
  Use Pole Line Road.

 See also
  Olive Drive
  Richards Blvd.
  Richards Underpass

 Proposed Solutions

 Roundabout
David Zetland cites this Wiki entry in his paper http://www.kysq.org/pubs/roundabouts.pdf Roundabouts in Davis: A Comprehensive Policy Analysis, and proposes that a roundabout in place of the current form of the intersection would improve the experience of commuters and noncommuters alike.
   Theres an Environmental Impact Report the city commissioned in 1996 on this in the Shields Library. DOCSTA CA/DAV C6 R52

 Bike and Pedestrian Bridge
Mitagating solution, but needs willingness to spend the money:  This is the worse intersection in the city. Lets push the city council to build a bridge across Richards for bikes and pedestrians, exclusively. Then disallow bikes and pedestrians from crossing either Olive or Richards at that intersection. The bridge, and the trail leading up to it from either side, will have a zero grade from one side to the other and will easily connect to other bike trails and sidewalks. We can build it parallel to the existing Union Pacific rail bridge without disturbing it.  We would need to take some of the Design Houses parking lot, just along the edge, but could compensate them with some of the towns property behind the Design House.  I am quite sure the town owns the property on the other side of Richards that is use to pile up old cars.  I encourage all of you to take a look at this area and see what Im talking about.  Its pedestrian and bike traffic, which I respect more than auto traffic, that is the biggest factor making this intersection a nightmare.  So lets take them out of the equation with a safe and attractive alternative.  Sure it will cost, but federal money floods into this town all the time. Lets use some of it for this.  If you agree, start making noise about this in the paper and at the City Council meetings. Users/WayneSchiller Wayne Schiller

  Several years ago the Davis City Council voted on whether or not to expand the tunnel and the interesection so it would be more pedestrian and bike friendly. Unfortunately, they believe the tunnel has historical value and should not be changed. So the next time you almsot die trying to cross the intersection just remember that Davis history is worth more than you. It was meant as a joke, but it is also a half truth. If there is money to be made out of destroying history the Davis City Council does not have a problem with it. A few years before the underpass vote, the council pushed the Borders  shopping area....well in building the Borders area and the nearby housing development an ancient Indian village site had to be covered up. So I guess I will have to revise my statement in that City of Davis history is worth more than Native American history. Also, several years ago downtown residents wanted to fill in the many pothalls that dot the alleyways in central Davis. Residents were getting injured walking and their cars were getting beat up. Well the Council thought the potholes were historical, so they have never been covered up.  Users/JimSchwab

Nice bit of revisionist history here, but not true unfortunately. The widening of Richards was not rejected by the city council, it was rejected by votes of the people: Three Times. The most recent came after the city council voted to widen the intersection despite the previous two votes against doing so. In their weknowwhatsgoodforyou arrogance, the council went forward with the atrocious reconstruction of the south side of the intersection as though the widening would occur. When it was rejected by yet another overwhelming vote, we were stuck with a stupid traffic pattern that simply does not work.

There were lots of good reasons for not widening Richards. The historical significance of the railroad bridge is a relatively minor one. The primary reason is that downtown Davis cannot handle the increased traffic flow that would result. Every traffic study that was done was clear on this point.

That is so true, but my above bridge proposal doesnt even have to touch the UP bridge, it will just run parallel to the historic bridge.  I like the idea of preserving history, but not at the expense of safety or urban sprawl. Users/WayneSchiller Wayne Schiller

If the bridge runs parrell to the tracks, how do you get across the tracks? Can you sketch this? Im confused.  Users/ArlenAbraham arlen

The bridge doesnt cross the tracks it only crosses over Richards in the same direction as the rail bridge does.  To cross the tracks one would follow the new path from the new bridge on the south side (the NAPA Auto parts side) of Richards back towards the Olive Richards intersection that would connect with the existing bike path that goes under the tracks into downtown.  Others could use the same new path on the south side to go east on East Olive and unto the existing bike paths that runs to campus behind Borders, or towards south Davis.  The bridge is constructed on the south side of the rail bridge. Is that clear? I wouldnt know how to submit a sketch here. Users/WayneSchiller Wayne Schiller

Is this what you mean Wayne?Users/MattJojola

YES, THATS EXACTLY WHAT I MEAN! Very good, Matt.  Hey everyone, doesnt it make sense? Users/WayneSchiller Wayne Schiller

Wouldnt the money be better spent installing a bike crossing (or even a full crossing) at the L and 2nd Garden and improving the Olive Richards intersection? The intersection improvements could include a left hand turn lane off olive headed downtown that would separate out a portion of the traffic that makes riding a bike westbound on Olive through the intersection so funky. Fixing the intersection would also improve things for bikes traveling on Richards.  Users/JasonAller

Good rendition Matt!  Waynes idea was recommended to the Planning Commission and the City Council two years ago.  However, the search through university civil engineering courses for a solution are still conducted.  None of the student groups have come up with this idea yet!  The problem with putting a bike/pedestrian crossing at the L & 2nd St crossing is that the city has to purchase an unused railroad crossing arm from another site somewhere in California.  UP doesnt just give these out willynilly.  The city currently owns the empty property in front of NAPA Auto parts and the other side behind Design House.  Making the idea above a fabulous solution!  Since UP agreed to that undercrossing for the new bike/ped path on the West Olive Dr. side, why wouldnt they be amenable to a bike crossing next to the tracks???!!!  Lastly, I dont agree with this:

...just remember that Davis history is worth more than you. KellyT

I have informally studied this intersection for about three years and have concluded that the two main problems areas are turning left off of Olive onto Richards south and the back up under the tunnel.  Turning left onto Ruchards north, into downtown, hasnt been that big of a problem because not that many cars come from Olive east. I truly believe that its just a bad place to mix cars and bikes/pedestrians. Many cities in the Bay Area have prohibited pedestrians from certain problem intersections.  I also think if bikes ride with cars, instead of using the cross walks, that would help, and we certainly should never prohibit bikes from the roads if they use them properly. I am not against a bike crossing at L, but fat chance getting UP to agree. Users/WayneSchiller Wayne Schiller

Jason, thats a beautiful aerial photograph you got there (from where)? I compared it with the printed UC Davis campus map, and one can see significant distortions in their map. For instance, they show Richards as being straight, Olive without the curve, and the Putah Creek Bicycle Undercrossing as being 3 blocks long! Users/SteveDavison

I just wanted to tell of my history with this intersection.  Being new to Davis and seeing that this intersection was busy, I decided to actually follow the pedestrian walk signs on my way back from InNOut.  No sooner did I get into the opposite lane (I did check for traffic, mind you) did I come about six inches from being run over by some car going what felt like 40mph.  I vote for the pedestrian bridge (or at least let me walk on the other side of the road under the bridge). Users/TimMonk

Speaking of when trucks get stuck, this happened to me one summer day while riding the W line back from work. Normally this wouldnt be so bad, but I was already late for my sisters high school graduation in Roseville and it also happened to be about 100 degrees outside. I ended up having to walk all of the way from Borders to Avalon, but not before snapping this:

Users/LaineWalkerAvina
 That is a really good picture. Users/NickSchmalenberger

20080223 16:31:14 nbsp I dont know much about traffic controls and stuff, but wouldnt making the light for cars headed West on Olive stay green a little longer fix a lot of problems?  Even like 10 seconds longer would let more cars through.  When pedestrians are crossing sometimes only 1 car can get through, which makes more cars back up.  People from InNOut and Shell cant get out, but try anyway which tends to block Eastbound traffic, which means more cars backed up going North on Richards.  That in turn makes exiting InNOut / Shell from the Western entrance/exit seem impossible (at least to me) and makes more cars come out on Olive drive, which kinda starts the problem over from the beginning.  Users/BradBenedict

20081115 16:43:45 nbsp Almost every day at noon, on Saturdays and general days with more heavy traffic getting off the Freeway from Olive I am blocked from getting out of my residence by an almost mile long back up of rude impatient drivers.  People exit the freeway onto Olive so they can drive into Downtown, and get backed up on the worlds most obnoxious light.  Both sides are always backed up, but the left turn lane more so than the right.  Neither has enough time to clear traffic manageably.  ITS THE TOWN JOKE!

And RATHER than spend a little money to retime the lights to make the Olive drive thrutraffic longer, they hire an old crossing guard to make sure that NO CARS can get through, and were screwed all the way down to Olive Drive Market.

When I am coming into my residence I cannot even enter, but have to wait, blocking traffic behind me because of the ridiculous backup of cars.  Not to mention the fact that there is almost NO WHERE to uturn down Olive Drive.

Let me also point out that if a heavy influx of traffic INTO downtown were such a problem, then why cant they make the Olive Drive thru traffic light LONGER to prevent it.  There have only been a FEW very occasional times when Richards became backed up going into Downtown.  We are almost always backed up here on Olive.  It makes me so angry.  And as much as it is unfair, it makes me hate large crowds of pedestrians and people who walk their bikes.  If we had any time at all to turn left none of this would be a problem.  75% of the time I am FORCED to drive through downtown rather than take the freeway around because the left turn is so fing backed up.

Fire the Davis city planners.  PLEASE! Users/CarrieBishop

20091023 03:47:10 nbsp It seems like if they just added a green arrow for the Olive traffic, things would be so much nicer. Users/sonictrance

20110808 10:54:26 nbsp I want to add my rant here as well.  I think many of the comments here already prove my point, so forgive me if this seems redundant: There have been numerous attempts over the years to improve this horrendous intersection (specifically, the subway underpass into downtown), that hasnt changed since 1917. Yes, you read that right: 1917. As in almost 100 years ago WWWhttp://en.wikipedia.org/wiki/Davis,_California. The main artery into downtown has not changed since the population of Davis was less than 1,000. The last attempt (I believe), was done in the mid90s. A very tasteful, well crafted expansion was proposed (my opinion of course), and put on the ballot. If failed, due to the locals who played the population with scare tactics, spouting asinine statements such as big rigs will start coming downtown! traffic will increase downtown! Itll be like a local freeway! So, the measure failed. Sound familiar? Voters being frightened off by a few citizens drumming up questionable facts to scare off the populace? See Covell Village. Meanwhile, the facts remain that emergency vehicles could potentially not get through this bottleneck if/when it floods in the winter (this has happened). During peak traffic hours, anyone whos ever tried to get through this intersection can be stuck for 10 minutes or more. The bike tunnel is great, but in my humble opinion, this continues to be the single worst traffic spot in the whole city. The premise of slowing people down through the intersection by creating gridlock just isnt the right way to do things. I would also qualify my remarks by saying I am completely for most every conceivable effort to retain anything linking Davis to its history. However, to me, the argument against expanding this intersection rings more of a misplaced perspective on city history (see historic potholes: WWWhttp://cityofdavis.org/cdd/cultural/30years/chapt10.cfm) than of a valued piece of our past. I realize now that the site has been cited in the National Register of Historic Places (WWWhttp://www.cr.nps.gov/nr/), there is little chance of it ever being tampered with (WWWhttp://cityofdavis.org/bicycles/tour.cfm). This being said, I still truly believe that, done tastefully, CORRECTLY, that there could be a site that fits the city in both appearance and function. I sincerely invite any discussion/opposition to any of the above.  Users/OldDavis73

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
C
C	$Id$
C
      DOUBLE PRECISION FUNCTION CSSIG2 (N1,N2,N,X,Y,Z,H,LIST,
     .                                LPTR,LEND,GRAD,TOL,
     .                                IFLGS, SIGMA, IER)
      INTEGER N1, N2, N, LIST(*), LPTR(*), LEND(N), IFLGS,
     .        IER
      DOUBLE PRECISION X(N), Y(N), Z(N), H(N), GRAD(3,N),
     .                 TOL, SIGMA(*)
C
C***********************************************************
C
C                                              From SSRFPACK
C                                            Robert J. Renka
C                                  Dept. of Computer Science
C                                       Univ. of North Texas
C                                           renka@cs.unt.edu
C                                                   07/25/96
C
C   Given a triangulation of a set of nodes on the unit
C sphere, along with data values H and gradients GRAD at the
C nodes, this function determines the smallest tension fac-
C tor CSSIG2 such that the Hermite interpolatory tension
C spline H(A), defined by CSSIG2 and the endpoint values and
C directional derivatives associated with an arc N1-N2,
C preserves convexity (or concavity) of the data:
C
C   HP1 .LE. S .LE. HP2 implies HPP(A) .GE. 0, and
C   HP1 .GE. S .GE. HP2 implies HPP(A) .LE. 0
C
C for all A in the open interval (A1,A2) corresponding to
C the arc, where HP1 and HP2 are the derivative values of H
C at the endpoints, S is the slope of the linear interpolant
C of the endpoint data values, HPP denotes the second deriv-
C ative of H, and A is arc-length.  Note, however, that
C infinite tension is required if HP1 = S or HP2 = S (unless
C HP1 = HP2 = S).
C
C On input:
C
C       N1,N2 = Nodal indexes of the endpoints of an arc for
C               which the tension factor is to be computed.
C               The indexes must be distinct and lie in the
C               range 1 to N, and if IFLGS .GE. 1, they must
C               correspond to adjacent nodes in the triangu-
C               lation.
C
C       N = Number of nodes in the triangulation.  N .GE. 3.
C
C       X,Y,Z = Arrays of length N containing coordinates of
C               the nodes.  X(I)**2 + Y(I)**2 + Z(I)**2 = 1.
C
C       H = Array of length N containing data values at the
C           nodes.  H(I) is associated with (X(I),Y(I),Z(I))
C           for I = 1 to N.
C
C       LIST,LPTR,LEND = Data structure defining the trian-
C                        gulation.  Refer to STRIPACK
C                        Subroutine CSTRMESH.
C
C       GRAD = Array dimensioned 3 by N whose columns con-
C              gradients at the nodes.  GRAD( ,J) must be
C              orthogonal to node J:  GRAD(1,J)*X(J) +
C              GRAD(2,J)*Y(J) + GRAD(3,J)*Z(J) = 0.  Refer
C              to Subroutines CSGRADG, CSGRADL, and CSSMSURF.
C
C       TOL = Tolerance whose magnitude determines how close
C             CSSIG2 is to its optimal value when nonzero
C             finite tension is necessary and sufficient to
C             satisfy convexity or concavity.  In the case
C             convexity, CSSIG2 is chosen so that 0 .LE.
C             HPPMIN .LE. abs(TOL), where HPPMIN is the
C             minimum value of HPP on the arc.  In the case
C             of concavity, the maximum value of HPP satis-
C             fies -abs(TOL) .LE. HPPMAX .LE. 0.  Thus, the
C             constraint is satisfied but possibly with more
C             tension than necessary.
C
C       IFLGS = Tension array option indicator:
C               IFLGS .LE. 0 if SIGMA is not to be used.
C               IFLGS .GE. 1 if SIGMA is to be updated by
C                            storing CSSIG2 in the appropriate
C                            locations.
C
C The above parameters are not altered by this function.
C
C       SIGMA = Dummy parameter (IFLGS .LE. 0) or array
C               containing tension factors associated with
C               arcs in one-to-one correspondence with LIST
C               entries (IFLGS .GE. 1).  Refer to Subroutine
C               CSGETSIG.
C
C On output:
C
C       SIGMA = Tension factor array updated with the new
C               value if and only if IFLGS .GE. 1 and IER
C               .GE. 0.
C
C       IER = Error indicator:
C             IER = 0 if no errors were encountered and fin-
C                     ite tension is sufficient to satisfy
C                     convexity (or concavity).
C             IER = 1 if no errors were encountered but in-
C                     finite tension is required to satisfy
C                     convexity.
C             IER = 2 if the data does not satisfy convexity
C                     or concavity.
C             IER = -1 if N1, N2, or N is outside its valid
C                      range.
C             IER = -2 if nodes N1 and N2 coincide or IFLGS
C                      .GE. 1 and the nodes are not adja-
C                      cent.
C
C       CSSIG2 = Minimum tension factor defined above unless
C              IER < 0, in which case CSSIG2 = -1.  If IER
C              = 1, CSSIG2 is set to 85, resulting in an
C              approximation to the linear interpolant of
C              the endpoint values.  If IER = 2, CSSIG2 = 0,
C              resulting in the Hermite cubic interpolant.
C
C STRIPACK module required by CSSIG2:  CSSTORE
C
C SSRFPACK modules required by CSSIG2:  CSARCLEN, CSSNHCSH
C
C Intrinsic functions called by CSSIG2:  ABS, EXP, MAX, MIN,
C                                        SQRT
C
C***********************************************************
C
      DOUBLE PRECISION CSARCLEN, CSSTORE
      INTEGER LP1, LP2, LPL, LUN, NIT
      DOUBLE PRECISION AL, COSHM, D1, D1D2, D2, DSIG, DUMMY,
     .                 EMS, F, FP, FTOL, P1(3), P2(3), RTOL,
     .                 S, SBIG, SIG, SINHM, SSM, T, T1, TP1,
     .                 UN(3), UNORM
C
      DATA SBIG/85./,  LUN/-1/
C
C Print a heading.
C
      IF (LUN .GE. 0) WRITE (LUN,100) N1, N2
  100 FORMAT (//1X,'CSSIG2 -- N1 =',I4,', N2 =',I4)
C
C Test for errors and set local parameters.
C
      IER = -1
      IF (MIN(N1,N2) .LT. 1  .OR.  N1 .EQ. N2  .OR.
     .    MAX(N1,N2,3) .GT. N) GO TO 11
      IER = -2
      IF (IFLGS .GT. 0) THEN
C
C   Set LP1 and LP2 to the pointers to N2 as a neighbor of
C     N1 and N1 as a neighbor of N2, respectively.
C
        LPL = LEND(N1)
        LP1 = LPTR(LPL)
    1   IF (LIST(LP1) .EQ. N2) GO TO 2
          LP1 = LPTR(LP1)
          IF (LP1 .NE. LPL) GO TO 1
        IF (ABS(LIST(LP1)) .NE. N2) GO TO 11
C
    2   LPL = LEND(N2)
        LP2 = LPTR(LPL)
    3   IF (LIST(LP2) .EQ. N1) GO TO 4
          LP2 = LPTR(LP2)
          IF (LP2 .NE. LPL) GO TO 3
        IF (ABS(LIST(LP2)) .NE. N1) GO TO 11
      ENDIF
C
C Store nodal coordinates P1 and P2, compute arc-length AL
C   and unit normal UN = (P1 X P2)/UNORM, and test for
C   coincident nodes.
C
    4 P1(1) = X(N1)
      P1(2) = Y(N1)
      P1(3) = Z(N1)
      P2(1) = X(N2)
      P2(2) = Y(N2)
      P2(3) = Z(N2)
      AL = CSARCLEN(P1,P2)
      UN(1) = P1(2)*P2(3) - P1(3)*P2(2)
      UN(2) = P1(3)*P2(1) - P1(1)*P2(3)
      UN(3) = P1(1)*P2(2) - P1(2)*P2(1)
      UNORM = SQRT(UN(1)*UN(1) + UN(2)*UN(2) + UN(3)*UN(3))
      IF (UNORM .EQ. 0.  .OR.  AL .EQ. 0.) GO TO 11
C
C Compute first and second differences and test for infinite
C   tension required.
C
      S = H(N2) - H(N1)
      D1 = S - AL*(GRAD(1,N1)*P2(1) + GRAD(2,N1)*P2(2) +
     .             GRAD(3,N1)*P2(3))/UNORM
      D2 = -AL*(GRAD(1,N2)*P1(1) + GRAD(2,N2)*P1(2) +
     .          GRAD(3,N2)*P1(3))/UNORM - S
      D1D2 = D1*D2
      IER = 1
      SIG = SBIG
      IF (D1D2 .EQ. 0.  .AND.  D1 .NE. D2) GO TO 10
C
C Test for a valid constraint.
C
      IER = 2
      SIG = 0.
      IF (D1D2 .LT. 0.) GO TO 10
C
C Test for SIG = 0 sufficient.
C
      IER = 0
      IF (D1D2 .EQ. 0.) GO TO 10
      T = MAX(D1/D2,D2/D1)
      IF (T .LE. 2.) GO TO 10
C
C Find a zero of F(SIG) = SIG*COSHM(SIG)/SINHM(SIG) - (T+1).
C   Since the derivative of F vanishes at the origin, a
C   quadratic approximation is used to obtain an initial
C   estimate for the Newton method.
C
      TP1 = T + 1.
      SIG = SQRT(10.*T-20.)
      NIT = 0
C
C   Compute an absolute tolerance FTOL = abs(TOL) and a
C     relative tolerance RTOL = 100*Macheps.
C
      FTOL = ABS(TOL)
      RTOL = 1.
    5 RTOL = RTOL/2.
        IF (CSSTORE(RTOL+1.) .GT. 1.) GO TO 5
      RTOL = RTOL*200.
C
C Top of loop:  evaluate F and its derivative FP.
C
    6 IF (SIG .LE. .5) THEN
C
C   Use approximations designed to avoid cancellation error
C     in the hyperbolic functions.
C
        CALL CSSNHCSH (SIG, SINHM,COSHM,DUMMY)
        T1 = COSHM/SINHM
        FP = T1 + SIG*(SIG/SINHM - T1*T1 + 1.)
      ELSE
C
C   Scale SINHM and COSHM by 2*exp(-SIG) in order to avoid
C     overflow.
C
        EMS = EXP(-SIG)
        SSM = 1. - EMS*(EMS+SIG+SIG)
        T1 = (1.-EMS)*(1.-EMS)/SSM
        FP = T1 + SIG*(2.*SIG*EMS/SSM - T1*T1 + 1.)
      ENDIF
C
      F = SIG*T1 - TP1
C
C   Update the number of iterations NIT.
C
      NIT = NIT + 1
      IF (LUN .GE. 0) WRITE (LUN,110) NIT, SIG, F, FP
  110 FORMAT (1X,3X,I2,' -- SIG = ',E15.8,', F = ',
     .        E15.8/1X,31X,'FP = ',E15.8)
C
C   Test for convergence.
C
      IF (FP .LE. 0.) GO TO 10
      DSIG = -F/FP
      IF (ABS(DSIG) .LE. RTOL*SIG  .OR.  (F .GE. 0.  .AND.
     .    F .LE. FTOL)  .OR.  ABS(F) .LE. RTOL) GO TO 10
C
C   Bottom of loop:  update SIG.
C
      SIG = SIG + DSIG
      GO TO 6
C
C No errors encountered.
C
   10 CSSIG2 = SIG
      IF (IFLGS .LE. 0) RETURN
      SIGMA(LP1) = SIG
      SIGMA(LP2) = SIG
      RETURN
C
C Error termination.
C
   11 CSSIG2 = -1.
      RETURN
      END

max_stars_repo_name: xylar/cdat

max_issues_repo_name: xylar/cdat

max_forks_repo_name: CDAT/uvcdat


text:
subroutine GeneralRestart(a,array,itask)
   !------------------------------------------------------------------------
   !    This routine reads the initial values from the restart file,
   !    interpolates them if necessary and writes the new restart
   !------------------------------------------------------------------------
   use typre
   use def_parame
   use Mod_iofile
   use Mod_int2str
   use Mod_PhysicalProblem
   implicit none
   class(PhysicalProblem) :: a
   integer(ip), intent(in)    :: itask

   integer(ip) :: kfl_twost, kfl_twost_pre,kfl_trasg_pre, kfl_repro_pre
   integer(ip) :: ndime,npoin,oldNpoin,ipoin,idofn,ndofn
   real(rp)    :: array(:,:)
   real(rp), allocatable :: auxarray(:,:)

   call a%Timer%Total%Tic
   call a%Timer%Restar%Tic

   select case (itask)

   case (1)

      call a%Mesh%GetNdime(ndime)
      call a%Mesh%GetNpoin(npoin)
      ndofn = size(array,1)

      if (a%kfl_inter /= 0_ip) then
         call a%OldMesh%GetNpoin(oldNpoin)
      else
         oldNpoin = npoin
      endif

      call a%Memor%alloc(ndofn,oldNpoin,auxarray,'auxarray','php_GeneralRestart')

      ! Read the previous time scheme used: one integer number
      !read(a%lun_rsta2,*) kfl_twost, kfl_trasg_pre, kfl_repro_pre

      read(a%lun_rstar) ((auxarray(idofn,ipoin),idofn=1,ndofn),ipoin=1,oldNpoin)

      if (a%kfl_inter/=1_ip) then
         array = auxarray
      else
         call a%Int_Restart%Interpolate(ndofn,auxarray,array)
      endif

      call a%Memor%dealloc(ndofn,oldNpoin,auxarray,'auxarray','php_GeneralRestart')

   case (2)

      call a%Mesh%GetNpoin(npoin)
      call a%Mesh%GetNdime(ndime)
      ndofn = size(array,1)

      ! Write the previous time scheme used: one integer number
      !write(a%lun_rsta2,*) kfl_twost, kfl_trasg_pre, kfl_repro_pre

      write(a%lun_rstar) ((array(idofn,ipoin),idofn=1,ndofn),ipoin=1,npoin)

   end select

   call a%Timer%Total%Toc
   call a%Timer%Restar%Toc

100 format(5x,'ERROR:   ',a)
101 format(5x,'WARNING: ',a)

end subroutine GeneralRestart

max_stars_repo_name: ciaid-colombia/InsFEM

max_issues_repo_name: ciaid-colombia/InsFEM

max_forks_repo_name: ciaid-colombia/InsFEM


text:
      program capron03

      logical a(10)
      integer c

      read *, a

      c = 0
      do i = 1, n
         c = c + a(i)
      enddo

      print *, c
      end


text:
!	CHTEST.F   CALCULATION OF CHEBYSHEV COEFFICIENT AND APPROX.

implicit double precision (a-h, o-z)
parameter(ns = 200, nc = 100, ntst = 1001, mxi = 350, nn = 319, mfit = 10, nwcy = 10, &
        nd = 10, xlambd = 0.0005543d00, a = 0.d00, perc = 0.01d00, cht0 = 1.d-4)
dimension c(nc), tt(2, 0:ntst), d0(nd), r39(0:ns), xi(0:2 * nc, nd, mxi)
dimension xmat(nd, 0:ns, mxi), lc(mxi), vc(nd), tilab(ns), telab(ns)
dimension ys(0:ns), wt(ns), lista(nc), covar(nc, nc), alpha(nc, nc)
dimension dyda(nc), chq(nwcy), c0(nc)
common /cont/ ncyc, nmaxi(nwcy), nmaxo(nwcy)
common /var/ tti(nwcy, 2, ntst), tto(nwcy, 2, ntst), &
        agei(nwcy, 2, ns), ageo(nwcy, 2, ns)
double precision slope
external slope
character*9  tab1
tab1 = char(09)
nn1 = nn

!	OUTPUT FILES

open(unit = 20, file = 'param-mon.out')
open(unit = 21, file = 'mages-sd.samp')
open(unit = 32, file = 'mchisq.dat')
open(unit = 40, file = 'mchistall-inp.dat')
open(unit = 42, file = 'mchistall-out.dat')
open(unit = 44, file = 'magesall-inp.dat')
open(unit = 46, file = 'magesall-out.dat')

!	INPUT FILES
open(unit = 22, file = 'age.in', status = 'old')
open(unit = 24, file = 'sig.in', status = 'old')
open(unit = 30, file = 'arr-me.in', status = 'old')
open(unit = 17, file = 'temstep.in', status = 'old')

open(unit = 51, file = 'autoagemon.cl', status = 'old')

!     INITIALIZATION OF OUTPUT FILES FOR XVGR SETS
kc0 = 48
kc1 = 48
mcyc = 0
call wgraph(40, kc1)
call wgraph(42, kc1)
call wgraph(44, kc1)
call wgraph(46, kc1)

!      print *, 'Insert Number of runnings'
read (51, *) nrun
!       print *, 'Insert max. plateau age'
read (51, *) tt(1, 1)

agein = tt(1, 1)
idem = -tt(1, 1) * nrun
np = 3
tt(2, 1) = 500.d00
tt(2, 2) = 300.d00
tt(1, 2) = tt(1, 1) / 2.d0
tt(2, 3) = 0.d00
tt(1, 3) = 0.d0
ncyc = 1
chqmin = 1.0d20
read(17, *)ni
do 12 j = 1, ni
    r39(j) = 0.d0
    read(17, *)telab(j), tilab(j)
    if(telab(j).gt.1373.)goto 11
12      tilab(j) = tilab(j) / 5.256d11
11      ni = j - 1
close(17)

read(30, *)nst
do 14 js = 1, nst
    read(30, *)e, d0(js), vc(js)
    d0(js) = 10.d0**d0(js)
    d0(js)= d0(js)/4.d0 *(24.d0 *3600.d0 *365d06)
    14 end do

close(30)
!     CALCULATION OF LC
nchini = 0
nsq = 2
nsum = 0
m = -1
do 32 mi = 1, nn
    if(nsum.gt.20)then
        nsum = 1
        nsq = 2 * nsq
    endif
    m = m + nsq
    if(m.gt.200)nsum = nsum + 1
    lc(mi) = m
32    continue

call lab(r39, d0, e, vc, telab, tilab, lc, xmat, ni, nst, nn)
call agesamp(r39, ys, wt, ni)
call chebft(a, tt(1, 1), c, nc, slope, tt, np)
do 13 i = 1, mfit
    lista(i) = i
    c0(i) = c(i)
13       continue
34     call ranchist(c, c0, mfit, idem, tt(2, 1))
dch = 0d00
alamda = -1
call mrqmin(r39, ys, wt, ni, c, nc, lista, mfit, covar, alpha, &
        nc, chisq, alamda, tt, xmat, lc, vc, d0, E, nst, nn1)
k = 1
itst = 0
write(20, *)'  chi-squared:  ', ni * chisq
print *, '  chi-squared:  ', ni * chisq

alam = alamda
1       write(20, *)'it #=', k, '   alamda:', alamda
       print *,'it #=', k, '   alamda:', alamda
if(dch.gt.0.d0)then
    write(20, *)
    write(20, *)'  chi-squared:  ', ni * chisq
    print *,'  chi-squared:  ', ni * chisq
endif
k = k + 1
ochisq = chisq
call mrqmin(r39, ys, wt, ni, c, nc, lista, mfit, covar, alpha, &
        nc, chisq, alamda, tt, xmat, lc, vc, d0, E, nst, nn1)
dch = abs(ochisq - chisq)
cht = cht0 * chisq
if (dch.lt.cht.and.alamda.lt.alam) then
    itst = itst + 1
    write(20, *)'itst= ', itst
    print *,'itst= ', itst
endif
alam = alamda
if (itst.lt.3.) then
    goto 1
endif
write(20, *)'difference= ', (ochisq - chisq) * ni
print *,'difference= ', (ochisq - chisq) * ni
chq(ncyc) = chisq
aq = 0.5d0 * (ni - 2.d0)
q = gammq(aq, 0.5 * ni * chisq)
write(20, *)'chisq= ', ni * chisq, '  GOF= ', q
print *,'chisq= ', ni * chisq, '  GOF= ', q
alamda = 0.d0
c40 = 1.d0 - dexp(-xlambd * tt(1, 1))
call zita(c, nc, lista, mfit, lc, xi, tt(1, 1), d0, E, tt, nst&
        , nn, perc, alamda)

!         write(20,'(1x,t5,a,t12,a,t19,a,t26,a,t33,a,t40,a,t47,a,t54,
!     $    a,t61,a,t68,a)')'c(1)','c(2)','c(3)'
!     $	 ,'c(4)','c(5)','c(6)','c(7)','c(8)','c(9)','c(10)'
!           write(20,'(1x,10f7.3)') (c(i),i=1,10)
do 42 i = 1, ni
    call age(i, c, ymod, dyda, nc, lista, mfit, vc, xmat, r39, xi, &
            c40, nn, nst, ni, perc, alamda)
42    continue
print *, nchini
write(32, *) nchini, chisq
if (nchini.eq.0)chqmin = chisq

nchini = nchini + 1
write(20, *)
write(20, *) 'New cycle:', nchini
write(20, *)

if (ncyc.eq.nwcy.or.nchini.eq.nrun)then
    open(unit = 10, file = 'mchist-inp.dat')
    open(unit = 12, file = 'mchist-out.dat')
    open(unit = 18, file = 'mages-inp.dat')
    open(unit = 19, file = 'mages-out.dat')
    open(unit = 48, file = 'status.info')
    if(nchini.gt.nwcy)then
        call readfile(10)
        call readfile(12)
        call readfile(18)
        call readfile(19)
        call readfile(48)
    else
        call wgraph(10, kc0)
        call wgraph(12, kc0)
        call wgraph(18, kc0)
        call wgraph(19, kc0)
        write(48, *)'program initialization'
    endif
    if(ncyc.eq.nchini)then
        do 300 j = 1, ncyc
            if(chq(j).lt.chqmin)chqmin = chq(j)
        300            continue
    endif

    do 310 j = 1, ncyc
        kcyc = nchini - ncyc + j
        !   WRITE OF ALL SETS
        do 362 j1 = 1, ni
            write(44, 100)agei(j, 1, j1), tab1, agei(j, 2, j1)
            write(44, 100)agei(j, 1, j1 + 1), tab1, agei(j, 2, j1)
            write(46, 100)ageo(j, 1, j1), tab1, ageo(j, 2, j1)
            write(46, 100)ageo(j, 1, j1 + 1), tab1, ageo(j, 2, j1)
        362               continue
        do 322 j1 = 1, nmaxi(j)
            write(40, 100)tti(j, 1, j1), tab1, tti(j, 2, j1)
        322               continue
        do 342 j1 = 1, nmaxo(j)
            write(42, 100)tto(j, 1, j1), tab1, tto(j, 2, j1)
        342               continue
        if(kcyc.ne.nrun)then
            write(40, *)'&
                    &  cycle #', kcyc
            write(42, *)'&
                    &  cycle #', kcyc, ' chisq= ', chq(j) * ni
            write(44, *)'&
                    &  cycle #', kcyc
            write(46, *)'&
                    &  cycle #', kcyc
        endif
        if(kcyc / 30 * 30 - kcyc.eq.0.and.kcyc.ne.nrun)then
            kc1 = kc1 + 1
            call wgraph(40, kc1)
            call wgraph(42, kc1)
            call wgraph(44, kc1)
            call wgraph(46, kc1)
        endif

        !   WRITE OF SELECTED SETS
        if(chq(j).lt.chqmin * 1.5d0)then
            mcyc = mcyc + 1
            do 320 j1 = 1, nmaxi(j)
                write(10, 100)tti(j, 1, j1), tab1, tti(j, 2, j1)
            320               continue
            n1 = 0
            do 340 j1 = 1, nmaxo(j)
                write(12, 100)tto(j, 1, j1), tab1, tto(j, 2, j1)
                if(tto(j, 2, j1).lt.100.and.tto(j, 1, j1 - 1).lt.agein&
                        .and.n1.eq.0)then
                    agein = tto(j, 1, j1 - 1)
                    n1 = 1
                endif
            340          continue
            write(48, *)'total runs done= ', kcyc, ' &
                    sucessful runs= ', mcyc
            do 360 j1 = 1, ni
                write(18, 100)agei(j, 1, j1), tab1, agei(j, 2, j1)
                write(18, 100)agei(j, 1, j1 + 1), tab1, agei(j, 2, j1)
                write(19, 100)ageo(j, 1, j1), tab1, ageo(j, 2, j1)
                write(19, 100)ageo(j, 1, j1 + 1), tab1, ageo(j, 2, j1)
            360               continue
            write(10, *)'&  cycle #', kcyc
            write(12, *)'&  cycle #', kcyc, ' chisq= ', chq(j) * ni
            write(18, *)'&  cycle #', kcyc
            write(19, *)'&  cycle #', kcyc
            if(mcyc / 30 * 30 - mcyc.eq.0.and.kcyc.ne.nrun)then
                kc0 = kc0 + 1
                call wgraph(10, kc0)
                call wgraph(12, kc0)
                call wgraph(18, kc0)
                call wgraph(19, kc0)
            endif
        endif
    310            continue
    ncyc = 1
    close(10)
    close(12)
    close(18)
    close(19)
    close(48)
else
    ncyc = ncyc + 1
endif

if(nchini.ne.nrun)goto 34

!     REMOVE LAST '&' FROM OUTPUT FILES ???

write(20, *)'End of CH cycles, total= ', nrun, '  completed=', mcyc
close(12)
close(22)
close(24)
close(30)
call conf_int(agein, tt(1, 1))
stop 'End of CH cycles'
140   format(1x, 4(f12.4, a1))
120   format(1x, f5.2, 6(f12.4))
100   format(1x, 6(g20.8, a1))

end

!   SUBROUTINE RANCHIST

subroutine ranchist(c, c0, mc, idem, tini)
    implicit double precision (a-h, o-z)
    dimension c(mc), c0(mc)
    fsum2 = 0.d0
    do 30 l = 1, mc
        perc = 5.d0 * ran1(idem) * (0.6 - ran1(idem))
        c(l) = c0(l) + perc
    30       continue
    20    continue
    tini = (0.5 - ran1(idem)) * 200 + 400
    return
end

!   FUNCTION RAN1

double precision function ran1(idum)
    implicit double precision (a-h, o-z)
    save ix1, ix2, ix3, r
    dimension r(97)
    parameter (m1 = 259200, ia1 = 7141, ic1 = 54773, rm1 = 3.8580247e-6)
    parameter (m2 = 134456, ia2 = 8121, ic2 = 28411, rm2 = 7.4373773e-6)
    parameter (m3 = 243000, ia3 = 4561, ic3 = 51349)
    data iff /0/
    if (idum.lt.0.or.iff.eq.0) then
        iff = 1
        ix1 = mod(ic1 - idum, m1)
        ix1 = mod(ia1 * ix1 + ic1, m1)
        ix2 = mod(ix1, m2)
        ix1 = mod(ia1 * ix1 + ic1, m1)
        ix3 = mod(ix1, m3)
        do 11 j = 1, 97
            ix1 = mod(ia1 * ix1 + ic1, m1)
            ix2 = mod(ia2 * ix2 + ic2, m2)
            r(j) = (float(ix1) + float(ix2) * rm2) * rm1
        11      continue
        idum = 1
    endif
    ix1 = mod(ia1 * ix1 + ic1, m1)
    ix2 = mod(ia2 * ix2 + ic2, m2)
    ix3 = mod(ia3 * ix3 + ic3, m3)
    j = 1 + (97 * ix3) / m3
    if(j.gt.97.or.j.lt.1)then
        write(20, *)'ERROR(RAN1): (j.gt.97.or.j.lt.1)'
        stop 'ERROR(RAN1): (j.gt.97.or.j.lt.1)'
    endif
    ran1 = r(j)
    r(j) = (float(ix1) + float(ix2) * rm2) * rm1
    return
end

!   SUBROUTINE AGESAMP

subroutine agesamp(r39, ys, wt, ni)
    implicit double precision (a-h, o-z)
    parameter (ns = 200)
    dimension xs(0:ns), ys(0:ni), r39(0:ni + 1), ya(0:ns)
    dimension wt(ni + 1), sig(ns)
    j = 1
    1         read(24, *, END = 10)sig(j)
    if(sig(j).le.0.)then
        write(20, *)'ERROR(AGESAMP): (sig(j).le.0)'
        stop 'ERROR(AGESAMP): (sig(j).le.0)'
    endif
    read(22, *, END = 10)xs(j), ya(j)
    write(21, 100)xs(j - 1) * 100.d0, ya(j) + sig(j)
    write(21, 100)xs(j), ya(j) + sig(j)
    xs(j) = xs(j) / 100.d00
    wt(j) = 1.d00 / dsqrt(dabs(r39(j) - r39(j - 1))) * sig(j)
    j = j + 1
    goto 1
    10      continue
    do 90 k = 1, j - 1
        write(21, 100)xs(j - k) * 100., ya(j - k) - sig(j - k)
        write(21, 100)xs(j - k - 1) * 100., ya(j - k) - sig(j - k)
    90    continue
    write(21, 100)xs(0), ya(1) + sig(1)

    ! EXTRAPOLATION FROM NI FOR FITTING CALCULATION

    ya(0) = ya(1)
    xs(ni + 1) = 1.d00
    ys(ni + 1) = ys(ni)
    jold = 1
    do 30 k = 1, ni
        ys(k) = 0.d0
        do 20 j = jold, ni + 1
            if(xs(j - 1).lt.r39(k).and.xs(j).gt.r39(k))then
                if(jold.eq.j)then
                    ys(k) = ya(j)
                else
                    do 34 jm = jold + 1, j - 1
                        ys(k) = ya(jm) * (xs(jm) - xs(jm - 1)) + ys(k)
                    34          continue
                    ys(k) = ys(k) + ya(jold) * (xs(jold) - r39(k - 1)) + ya(j) * (r39(k) - xs(j - 1))
                    ys(k) = ys(k) / (r39(k) - r39(k - 1))
                    jold = j
                endif
                goto 30
            endif
        20      continue
    30      continue
    return
    100   format(1x, 4(f14.8))
end

!    SUBROUTINE LAB

subroutine lab(r39, d0, e, vc, telab, tilab, lc, xmat, ni, nst, nn)
    implicit double precision (a-h, o-z)
    parameter(pi = 3.14159265359d00, imp = 2, &
            b = 8.d0, ns = 200, R = 1.987d-03, nd = 10)
    dimension xmat(1:nd, 0:ns, 1:nn), lc(nn), zita(0:ns), r39(0:ni)
    dimension tilab(ni), telab(ni), d0(nst), vc(nst)
    !      character*9  tab1
    !      tab1=char(09)
    zita(0) = 0.d00
    nsq = 2
    nsq1 = 1
    do  60 j = 1, nn
        an2 = (lc(j) * pi)**2
        xmat(1, ni + 1, j) = b / an2
        if(j.ge.121)then
            if((j - 121) / 20 * 20 - (j - 121).eq.0)then
                nsq = nsq * 2
                nsq2 = nsq / imp
            endif
            do 64 mi = 1, nsq1 - 1
                an = ((lc(j) - mi * imp) * pi)**2
                xmat(1, ni + 1, j) = xmat(1, ni + 1, j) + b / an * (nsq1 - mi) / nsq1
            64         continue
            if(j.eq.nn)goto 60
            do 68 mi = 1, nsq2 - 1
                an = ((lc(j) + mi * imp) * pi)**2
                xmat(1, ni + 1, j) = xmat(1, ni + 1, j) + b / an * (nsq2 - mi) / nsq2
            68         continue
            nsq1 = nsq2
        endif
    60      continue
    do 10 k = 1, nst
        do 10 nt = 1, ni
            if(k.eq.1)r39(nt) = 0.d0
            zita(nt) = zita(nt - 1) + d0(k) * tilab(nt) * dexp(-e / (r * telab(nt)))
            nsq = 2
            nsq1 = 1
            do 48 j = 1, nn
                an2 = (lc(j) * pi)**2
                if(an2 * zita(nt).gt.60)then
                    do 70 jx = j, nn
                        xmat(k, nt, jx) = xmat(1, ni + 1, jx)
                    70                      continue
                    goto 72
                else
                    xmat(k, nt, j) = xmat(1, ni + 1, j) - b * dexp(-an2 * zita(nt)) / an2
                endif
                if(j.ge.121)then
                    if((j - 121) / 20 * 20 - (j - 121).eq.0)then
                        nsq = nsq * 2
                        nsq2 = nsq / imp
                    endif
                    do 42 mi = 1, nsq1 - 1
                        an = ((lc(j) - mi * imp) * pi)**2
                        xmat(k, nt, j) = xmat(k, nt, j) - b * dexp(-an * zita(nt))&
                                / an * (nsq1 - mi) / nsq1
                    42         continue
                    if(j.eq.nn)goto 50
                    do 44 mi = 1, nsq2 - 1
                        an = ((lc(j) + mi * imp) * pi)**2
                        xmat(k, nt, j) = xmat(k, nt, j) - b * dexp(-an * zita(nt))&
                                / an * (nsq2 - mi) / nsq2
                    44          continue
                    nsq1 = nsq2
                endif
                50           continue
            48      continue
            72             do 20 k1 = 1, 80000, imp
                an2 = (k1 * pi)**2
                r39(nt) = r39(nt) + b / an2 * vc(k)
                if (zita(nt) * (k1 * pi)**2.gt.100.) then
                    ab = 0.d00
                else
                    ab = dexp(-zita(nt) * (k1 * pi)**2) / (k1 * pi)**2
                    r39(nt) = r39(nt) - b * ab * vc(k)
                endif
            20       continue
        10    continue
    return
end

!   SUBROUTINE CHEBFT

subroutine chebft(a, b, c, n, func, tt, np)
    implicit double precision (a-h, o-z)
    parameter (nmax = 100, pi = 3.141592653589793d0)
    dimension c(n), f(nmax), tt(2, 0:np)
    bma = 0.5d00 * (b - a)
    bpa = 0.5d00 * (b + a)
    do 11 k = 1, n
        y = dcos(pi * (k - 0.5d0) / n)
        f(k) = func(tt, np, y * bma + bpa)
    11    continue
    fac = 2.d0 / n
    do 13 j = 1, n
        sum = 0.d0
        do 12 k = 1, n
            sum = sum + f(k) * dcos((pi * (j - 1)) * ((k - 0.5d0) / n))
        12      continue
        c(j) = fac * sum
    13    continue
    return
end

!    DOUBLE PRECISION FUNCTION SLOPE

double precision function slope(a, n, x)
    implicit double precision (a-h, o-z)
    dimension a(2, 0:n)
    do 10 j = 2, n
        if(a(1, j).lt.x.and.a(1, j - 1).gt.x)then
            slope = dsqrt(dabs((a(2, j - 1) - a(2, j)) / (a(1, j - 1) - a(1, j))))
            return
        endif
    10    continue
end

!    SUBROUTINE MRQMIN

subroutine mrqmin(x, y, sig, ndata, a, ma, lista, mfit, &
        covar, alpha, nca, chisq, alamda, tt, xmat, &
        lc, vc, d0, E, nst, nn)
    implicit double precision (a-h, o-z)
    parameter (mmax = 100, ntst = 1001, nd = 10, ns = 200)
    dimension x(0:ndata), y(0:ndata), sig(ndata), a(ma), lista(ma), &
            covar(nca, nca), alpha(nca, nca), atry(mmax), beta(mmax), da(mmax), &
            xmat(1:nd, 0:ns, 1:nn), d0(nst), vc(nst), tt(2, 0:ntst), lc(nn)&
            , covaux(mmax, mmax), daux(mmax)
    save beta, atry, da, ochisq
    if(alamda.lt.0.)then
        kk = mfit + 1
        do 12 j = 1, ma
            ihit = 0
            do 11 k = 1, mfit
                if(lista(k).eq.j)ihit = ihit + 1
            11        continue
            if (ihit.eq.0) then
                lista(kk) = j
                kk = kk + 1
            else if (ihit.gt.1) then
                write(20, *)'ERROR(MRQMIN): improper permutation in lista-1'
                stop 'ERROR(MRQMIN): improper permutation in lista-1'
            endif
        12      continue
        if (kk.ne.(ma + 1))then
            write(20, *)'ERROR(MRQMIN): improper perm. in lista-2'
            stop 'ERROR(MRQMIN): improper perm. in lista-2'
        endif
        call mrqcof(x, y, sig, ndata, a, ma, lista, mfit, alpha, beta, nca, chisq, &
                tt, xmat, lc, vc, d0, e, nst, nn, alamda)
        alamda = 0.001d0
        ochisq = chisq
        do 13 j = 1, ma
            atry(j) = a(j)
        13      continue
    endif
    do 15 j = 1, mfit
        do 14 k = 1, mfit
            covar(j, k) = alpha(j, k)
        14      continue
        covar(j, j) = alpha(j, j) * (1.d0 + alamda)
        da(j) = beta(j)
    15    continue
    ng = 1
    call gaussj(covar, mfit, nca, da, ng, 1)
    if(ng.eq.-1)stop 'Singular Matrix'
    if(alamda.eq.0.)then
        call covsrt(covar, nca, ma, lista, mfit)
        return
    endif
    do 16 j = 1, mfit
        if(alamda.ge.1.00d20)then
            da(j) = 0.d0
        endif
        atry(lista(j)) = a(lista(j)) + da(j)
    16    continue
    call mrqcof(x, y, sig, ndata, atry, ma, lista, mfit, covar, da, nca, chisq, &
            tt, xmat, lc, vc, d0, e, nst, nn, alamda)
    if(chisq.le.ochisq)then

        !  CHECK IF COVAR IS NOT SINGULAR
        do 19 j = 1, mfit
            do 20 k = 1, mfit
                covaux(j, k) = covar(j, k)
            20        continue
            daux(j) = da(j)
        19       continue
        call gaussj(covaux, mfit, nca, daux, ng, 1)
        if(ng.eq.-1)then
            alamda = 10.d0 * alamda
            chisq = ochisq
            ng = 1
            return
        endif
        !   END CHECKING FOR SINGULAR MATRIX

        alamda = 0.1d0 * alamda
        ochisq = chisq
        do 18 j = 1, mfit
            do 17 k = 1, mfit
                alpha(j, k) = covar(j, k)
            17        continue
            beta(j) = da(j)
            a(lista(j)) = atry(lista(j))
        18      continue
    else
        alamda = 10.d0 * alamda
        chisq = ochisq
    endif
    return
end

!     SUBROUTINE MRQCOF

subroutine mrqcof(x, y, sig, ndata, a, ma, lista, mfit, alpha, beta, nalp, &
        chisq, tt, xmat, lc, vc, d0, e, nst, nn, al)
    implicit double precision (a-h, o-z)
    parameter (mmax = 100, ntst = 1001, perc = 0.01d0, xlambd = 0.0005543d0, &
            nc = 100, mxi = 500, nd = 10, ns = 200)
    dimension x(0:ndata), y(0:ndata), sig(ndata), alpha(nalp, nalp), &
            dyda(mmax), lista(mfit), a(ma), xmat(1:nd, 0:ns, 1:nn), d0(nst), &
            vc(nst), tt(2, 0:ntst), lc(nn), beta(ma), xi(0:2 * nc, nd, mxi)
    t0 = tt(1, 1)
    c40 = 1.d0 - dexp(-xlambd * t0)
    do 12 j = 1, mfit
        do 11 k = 1, j
            alpha(j, k) = 0.d0
        11      continue
        beta(j) = 0.d0
    12    continue
    call zita(a, ma, lista, mfit, lc, xi, t0, d0, E, tt, nst, nn, perc, al)
    if(tt(2, ntst).eq.-1.)then
        tt(2, ntst) = 0.d0
        chisq = chisq * 2.d0
        return
    endif
    chisq = 0.d0
    do 15 i = 1, ndata
        call age(i, a, ymod, dyda, ma, lista, mfit, vc, xmat, x, xi, c40, nn, &
                nst, ndata, perc, al)
        sig2i = 1.d0 / (sig(i) * sig(i))
        dy = y(i) - ymod
        do 14 j = 1, mfit
            wt = dyda(lista(j)) * sig2i
            do 13 k = 1, j
                alpha(j, k) = alpha(j, k) + wt * dyda(lista(k))
            13        continue
            beta(j) = beta(j) + dy * wt
        14      continue
        chisq = chisq + dy * dy * sig2i
    15    continue
    do 17 j = 2, mfit
        do 16 k = 1, j - 1
            alpha(k, j) = alpha(j, k)
        16      continue
    17    continue
    return
end

!   SUBROUTINE COVSRT

subroutine covsrt(covar, ncvm, ma, lista, mfit)
    implicit double precision (a-h, o-z)
    dimension covar(ncvm, ncvm), lista(mfit)
    do 12 j = 1, ma - 1
        do 11 i = j + 1, ma
            covar(i, j) = 0.d0
        11      continue
    12    continue
    do 14 i = 1, mfit - 1
        do 13 j = i + 1, mfit
            if(lista(j).gt.lista(i)) then
                covar(lista(j), lista(i)) = covar(i, j)
            else
                covar(lista(i), lista(j)) = covar(i, j)
            endif
        13      continue
    14    continue
    swap = covar(1, 1)
    do 15 j = 1, ma
        covar(1, j) = covar(j, j)
        covar(j, j) = 0.d0
    15    continue
    covar(lista(1), lista(1)) = swap
    do 16 j = 2, mfit
        covar(lista(j), lista(j)) = covar(1, j)
    16    continue
    do 18 j = 2, ma
        do 17 i = 1, j - 1
            covar(i, j) = covar(j, i)
        17      continue
    18    continue
    return
end

!   SUBROUTINE GAUSSJ

subroutine gaussj(a, n, np, b, m, mp)
    implicit double precision (a-h, o-z)
    parameter (nmax = 100)
    dimension a(np, np), b(np, mp), ipiv(nmax), indxr(nmax), indxc(nmax)
    do 11 j = 1, n
        ipiv(j) = 0
    11    continue
    do 22 i = 1, n
        big = 0.d0
        do 13 j = 1, n
            if(ipiv(j).ne.1)then
                do 12 k = 1, n
                    if (ipiv(k).eq.0) then
                        if (dabs(a(j, k)).ge.big)then
                            big = dabs(a(j, k))
                            irow = j
                            icol = k
                        endif
                    else if (ipiv(k).gt.1) then
                        write(20, *)'ERROR(GAUSS): singular matrix-1'
                        m = -1
                        return
                    endif
                12          continue
            endif
        13      continue
        ipiv(icol) = ipiv(icol) + 1
        if (irow.ne.icol) then
            do 14 l = 1, n
                dum = a(irow, l)
                a(irow, l) = a(icol, l)
                a(icol, l) = dum
            14        continue
            do 15 l = 1, m
                dum = b(irow, l)
                b(irow, l) = b(icol, l)
                b(icol, l) = dum
            15        continue
        endif
        indxr(i) = irow
        indxc(i) = icol
        if (a(icol, icol).eq.0.) then
            write(20, *)'ERROR(GAUSS): singular matrix-2'
            m = -1
            return
        endif
        pivinv = 1.d0 / a(icol, icol)
        a(icol, icol) = 1.d0
        do 16 l = 1, n
            a(icol, l) = a(icol, l) * pivinv
        16      continue
        do 17 l = 1, m
            b(icol, l) = b(icol, l) * pivinv
        17      continue
        do 21 ll = 1, n
            if(ll.ne.icol)then
                dum = a(ll, icol)
                a(ll, icol) = 0.d0
                do 18 l = 1, n
                    a(ll, l) = a(ll, l) - a(icol, l) * dum
                18          continue
                do 19 l = 1, m
                    b(ll, l) = b(ll, l) - b(icol, l) * dum
                19          continue
            endif
        21      continue
    22    continue
    do 24 l = n, 1, -1
        if(indxr(l).ne.indxc(l))then
            do 23 k = 1, n
                dum = a(k, indxr(l))
                a(k, indxr(l)) = a(k, indxc(l))
                a(k, indxc(l)) = dum
            23        continue
        endif
    24    continue
    return
end

!   SUBROUTINE ZITA

subroutine zita(c, mc, lista, mfit, lc, xi, t0, d0, E, tt, nst, &
        nn, perc, al)
    implicit double precision (a-h, o-z)
    parameter(ntst = 1001, a = 0, nc = 100, nwcy = 10, ns = 200)
    dimension lc(nn), d0(nst), lista(mfit)
    dimension c(mc), tt(2, 0:ntst), xi(0:2 * nc, 1:nst, 1:nn)
    common /var/ tti(nwcy, 2, ntst), tto(nwcy, 2, ntst), &
            agei(nwcy, 2, ns), ageo(nwcy, 2, ns)
    common /cont/ ncyc, nmaxi(nwcy), nmaxo(nwcy)
    character*9  tab1
    tab1 = char(09)
    n = 0
    call sched(a, t0, c, mc, tt, imax)
    call geosec(imax, d0, E, tt, xi, nn, lc, nst, n)
    write(20, *)'# of CH points= ', imax
    if(al.le.0.)then
        if (al.lt.0) then
            nmaxi(ncyc) = imax
        else
            nmaxo(ncyc) = imax
        endif
        do 40 j = 1, imax
            if (al.lt.0) then
                tti(ncyc, 1, j) = tt(1, j)
                tti(ncyc, 2, j) = tt(2, j)
            else
                tto(ncyc, 1, j) = tt(1, j)
                tto(ncyc, 2, j) = tt(2, j)
            endif
        40      continue
    endif
    do 20 k = 1, mfit
        csk = c(lista(k))
        dc = perc * dabs(c(lista(k)))
        do 30 js = 1, 2
            n = n + 1
            c(lista(k)) = c(lista(k)) + (-1)**js * dc
            call sched2(a, t0, c, mc, tt, imax)
            call geosec(imax, d0, E, tt, xi, nn, lc, nst, n)
            c(lista(k)) = csk
        30        continue
    20    continue

    100   format(1x, 6(g20.8, a1))
    return
end

!    SUBROUTINE AGE

subroutine age(nt, c, y, dyda, mc, lista, mfit, vc, xmat, r39, xi, c40, &
        nn, nst, ndata, perc, al)
    implicit double precision (a-h, o-z)
    parameter(xlambd = 0.0005543d0, nc = 100, nd = 10, ns = 200, nwcy = 10, &
            ntst = 1001)
    dimension xmat(1:nd, 0:ns, 1:nn), r39(0:ndata), lista(mfit)
    dimension c(mc), dyda(mc), vc(nst), xi(0:2 * nc, 1:nst, 1:nn)
    common /var/ tti(nwcy, 2, ntst), tto(nwcy, 2, ntst), &
            agei(nwcy, 2, ns), ageo(nwcy, 2, ns)
    common /cont/ ncyc, nmaxi(nwcy), nmaxo(nwcy)
    character*9  tab1
    tab1 = char(09)
    c39wd = 1.d0 - c40
    s1 = 0.d0
    do 52 k = 1, nst
        do 52 m = 1, nn
            s1 = s1 + vc(k) * xi(0, k, m) * (xmat(k, nt, m) - xmat(k, nt - 1, m))
        52    continue
    dr39 = r39(nt) - r39(nt - 1)
    dr40 = dr39 * (c40 - 1) + s1
    y = 0.d0
    if(dr40 / dr39.gt.0.)then
        y = dlog(1.d0 + dr40 / dr39 / c39wd) / xlambd
        dya = 1.d0 / (dr39 * c39wd + dr40) / xlambd
    endif
    if(al.lt.0)then
        agei(ncyc, 1, nt + 1) = r39(nt) * 100.d0
        agei(ncyc, 2, nt) = y
    else if (al.eq.0.) then
        ageo(ncyc, 1, nt + 1) = r39(nt) * 100.d0
        ageo(ncyc, 2, nt) = y
    endif
    42    continue
    do 56 ls = 1, mfit
        l = ls * 2 - 1
        dc = perc * dabs(c(lista(ls)))
        s1 = 0.d0
        do 54 k = 1, nst
            do 54 m = 1, nn
                dxi = (xi(l + 1, k, m) - xi(l, k, m)) / (2 * dc)
                s1 = s1 + vc(k) * (xmat(k, nt, m) - xmat(k, nt - 1, m)) * dxi
            54    continue
        dyda(lista(ls)) = s1 * dya
    56    continue
    100   format(1x, 4(g20.8, a1))
    return
end

!     SUBROUTINE GEOSEC

subroutine geosec(imax, d0, E, tt, xi, nn, lc, nst, j1)
    implicit double precision (a-h, o-z)
    parameter(xlambd = 5.543d-04, R = 1.987d-03, max = 1002, &
            pi = 3.14159265359d0, nc = 100)
    dimension xi(0:2 * nc, 1:nst, 1:nn), tt(2, 0:imax), dzita(0:max)
    dimension lc(nn), d0(nst), d(0:max)
    !      character*9   tab1
    !      tab1=char(09)
    !      nn=319
    do 80 k = 1, nst
        nend = imax - 1
        do 10 j = 1, nend
            avtemp = (tt(2, j + 1) + tt(2, j)) / 2.d0 + 273.d0
            if(e / r / avtemp.gt.80)then
                d(j) = 0.d0
                goto 10
            endif
            d(j) = d0(k) * dexp(-e / r / avtemp)
        10    continue
        60    dzita(nend + 1) = 0.d0
        do 20 j = nend, 1, -1
        20    dzita(j) = dzita(j + 1) + dabs(tt(1, j + 1) - &
                tt(1, j)) * d(j)
        !     COMPUTO OF XI(M) - M<80000 --> nn=319
        do 30 mi = 1, nn
            xlogm = 2 * dlog(lc(mi) * pi)
            sum = 0.d0
            do 40 n = 1, nend
                if(d(n).eq.0.)goto 40
                uplus = (lc(mi) * pi)**2 * dzita(n + 1) + xlambd * (tt(1, 1) - tt(1, n + 1))
                if (uplus - xlogm.gt.25.)goto 40
                al = (lc(mi) * pi)**2 * d(n) - xlambd
                xal = al * (tt(1, n) - tt(1, n + 1))
                if (dabs(xal).gt.30)then
                    camal = 1.d0 / al
                else
                    if (dabs(xal).gt.0.001)then
                        camal = (1.d0 - dexp(-xal)) / al
                    else
                        camal = (tt(1, n) - tt(1, n + 1)) * (1 - xal / 2 + xal**2 / 6.d0)
                    endif
                endif
                sum = sum + d(n) * dexp(-uplus) * camal
            40        continue
            tzita = dzita(1) * (pi * lc(mi))**2
            if(tzita.lt.30.)then
                xfact = dexp(-tzita)
            else
                xfact = 0.
            endif
            xi(j1, k, mi) = sum * (pi * lc(mi))**2 + xfact
        30    continue
    80    continue
    return
end

!      SUBROUTINE SCHED

subroutine sched(a, t0, c, m, tt, i)
    implicit double precision (a-h, o-z)
    parameter(step = 4.d0, nt = 10000, ntst = 1001, tmax = 100.d0)
    dimension tt(2, 0:ntst), c(m)
    i = 1
    slp = (chebev(a, t0, c, m, t0))**2
    dt = t0 / nt
    sum = slp * dt / 2.d0
    do 10 j = 1, nt
        xj = t0 - j * dt
        if(j.eq.nt)xj = 0.d0
        slp = (chebev(a, t0, c, m, xj))**2
        sum = sum + slp * dt / 2
        if(sum.ge.step.or.xj.eq.0.)then
            tt(2, i + 1) = tt(2, i) - sum
            tt(1, i + 1) = xj
            sum = 0.d0
            i = i + 1
            if(tt(2, i).le.tmax)goto 20
        endif
        sum = sum + slp * dt / 2
    10    continue
    return
    20      tt(2, i) = 0.d0
    tt(1, i) = 0.d0
    return
end

subroutine sched2(a, t0, c, m, tt, i)
    implicit double precision (a-h, o-z)
    parameter(nt = 10000, ntst = 1001, tmax = 100.d0)
    dimension tt(2, 0:ntst), c(m)
    i = 1
    slp = (chebev(a, t0, c, m, t0))**2
    dt = t0 / nt
    sum = slp * dt / 2.d0
    do 10 j = 1, nt
        xj = t0 - j * dt
        if(j.eq.nt)xj = 0.d0
        slp = (chebev(a, t0, c, m, xj))**2
        sum = sum + slp * dt / 2
        if(xj.eq.tt(1, i + 1))then
            tt(2, i + 1) = tt(2, i) - sum
            sum = 0.d0
            i = i + 1
            if(tt(2, i).le.tmax)goto 20
        endif
        sum = sum + slp * dt / 2
    10    continue
    return
    20      tt(2, i) = 0.d0
    tt(1, i) = 0.d0
    return
end

!     SUBROUTINE CHEBEV

function chebev(a, b, c, m, x)
    implicit double precision (a-h, o-z)
    dimension c(m)
    if ((x - a) * (x - b).gt.0.) then
        write(20, *)'ERROR(CHEBEV): x not in range'
        stop 'ERROR(CHEBEV): x not in range'
    endif
    d = 0.d0
    dd = 0.d0
    y = (2.d0 * x - a - b) / (b - a)
    y2 = 2.d0 * y
    do 11 j = m, 2, -1
        sv = d
        d = y2 * d - dd + c(j)
        dd = sv
    11    continue
    chebev = y * d - dd + 0.5d0 * c(1)
    return
end

double precision function gammq(a, x)
    implicit double precision (a-h, o-z)
    if(x.lt.0..or.a.le.0.)then
        write(20, *)'ERROR(GAMMQ): (x.lt.0..or.a.le.0.)'
        stop 'ERROR(GAMMQ): (x.lt.0..or.a.le.0.)'
    endif
    if(x.lt.a + 1.)then
        call gser(gamser, a, x, gln)
        gammq = 1. - gamser
    else

        call gcf(gammcf, a, x, gln)
        gammq = gammcf
    endif
    return
end

subroutine gcf(gammcf, a, x, gln)
    implicit double precision (a-h, o-z)
    parameter (itmax = 100, eps = 3.e-7)
    gln = gammln(a)
    gold = 0.
    a0 = 1.
    a1 = x
    b0 = 0.
    b1 = 1.
    fac = 1.
    do 11 n = 1, itmax
        an = float(n)
        ana = an - a
        a0 = (a1 + a0 * ana) * fac
        b0 = (b1 + b0 * ana) * fac
        anf = an * fac
        a1 = x * a0 + anf * a1
        b1 = x * b0 + anf * b1
        if(a1.ne.0.)then
            fac = 1. / a1
            g = b1 * fac
            if(dabs((g - gold) / g).lt.eps)go to 1
            gold = g
        endif
    11    continue
    write(20, *)'ERROR(GCF): a too large, itmax too small'
    stop 'ERROR(GCF): a too large, itmax too small'
    1     gammcf = dexp(-x + a * dlog(x) - gln) * g
    return
end

subroutine gser(gamser, a, x, gln)
    implicit double precision (a-h, o-z)
    parameter (itmax = 100, eps = 3.e-7)
    gln = gammln(a)
    if(x.le.0.)then
        if(x.lt.0.)then
            write(20, *)'ERROR(GSER): (x.lt.0.)'
            stop 'ERROR(GSER): (x.lt.0.)'
        endif
        gamser = 0.
        return
    endif
    ap = a
    sum = 1. / a
    del = sum
    do 11 n = 1, itmax
        ap = ap + 1.
        del = del * x / ap
        sum = sum + del
        if(dabs(del).lt.dabs(sum) * eps)go to 1
    11    continue
    write(20, *)'ERROR(GSER): a too large, itmax too small'
    stop 'ERROR(GSER): a too large, itmax too small'
    1     gamser = sum * dexp(-x + a * dlog(x) - gln)
    return
end

!   SUBROUTINE WGRAPH
subroutine wgraph(n, k)
    character *2,graph
    if(k.lt.48.or.k.gt.57)return
    graph = 'g' // char(k)
    write(n, 120)'@WITH ' // graph
    write(n, 130)'@' // graph // ' ON'
    120    format(a8)
    130    format(a6)
    return
end

subroutine readfile(l)
    10      read(l, *, end = 20)
    goto 10
    20      return
end

!  CONF_INT.F  CALCULATES THE INTERVAL OF CONFIDENCE FROM
!  A SET OF THERMAL HISTORIES.
subroutine conf_int(agein, agend)
    implicit double precision (a-h, o-z)
    parameter(ns = 200, cn = 12, R = 1.987D-3, dt = 20.d0, tau1 = 3., tau2 = 1)
    parameter(tim = 3.15576d13, pi = 3.14159265, nt = 1000)
    dimension temp(nt, ns), time(nt, ns), xtem(ns), nch(ns)
    dimension xmed(ns), ave(ns), adev(ns), sdev(ns), xmed2(ns), &
            var(ns), skew(ns), curt(ns), xage(ns), xmed1(ns)
    external temfc
    open(unit = 22, file = 'mchist-out.dat', status = 'old')
    !      open(unit=24,file='statanal.dat')
    !      open(unit=26,file='average.dat')
    !      open(unit=28,file='confint.dat')
    open(unit = 29, file = 'confmed.dat')
    !      open(unit=30,file='input.dat')
    ncyc = 0
    ncurv = 0
    nsteps = 100
    dage = (agend - agein) / nsteps
    !  READ THERMAL HISTORIES
    ks = 1
    kend = 0
    1    n = 1
    2     read(22, *, END = 20, ERR = 30)time(n, ks), temp(n, ks)
    n = n + 1
    goto 2
    20    kend = -1
    30    nch(ks) = n - 1
    if (nch(ks).eq.0.and.kend.eq.0)goto 1
    if (nch(ks).eq.0.and.kend.eq.-1)goto 13
    ks = ks + 1
    if(ks.gt.ns)stop 'too many thermal histories'
    goto 1
    13    ks = ks - 1
    if(ks.lt.20)then
        write(29, *)'Not enough solutions to calculate &
                confidence intervals'
        return
    endif
    !  CALCULATES THE XTEM VECTOR FOR EACH TIME
    do 10 j = 1, nsteps + 1
        xage(j) = agein + (j - 1) * dage
        xmed(j) = 0
        do 22 k = 1, ks
            if (xage(j).gt.time(1, k))goto 10
            xtem(k) = temfc(time(1, k), temp(1, k), nch(k), xage(j))
        22        continue
        !	    STATISTICAL CALCULATIONS
        !	    MEAN OF THE DISTRIBUTION
        call moment(xtem, ks, ave(j), adev(j), sdev(j), &
                var(j), skew(j), curt(j))
        call mdian1(xtem, ks, xmed(j))
        nperc = ks * .05
        call binom(ks + 1, j1, j2)
        if(j.eq.1)then
            xmed0 = xtem(j2)
            xmed02 = xtem(ks - nperc)
        endif
        xmed1(j) = xtem(j1)
        xmed2(j) = xtem(nperc)
        !  WRITE INTERVAL OF CONFIDENCE OF MEDIAN and DISTRIBUTION
        !          write(24,110)xage(j),ave(j),xmed(j),sdev(j),adev(j),
        !     $      var(j),skew(j),curt(j)
        !          write(26,110)xage(j),ave(j),xmed(j),xmed(j)+sdev(j),
        !     $      xmed(j)-sdev(j)
        !          write(28,110)xage(j),xtem(j2),xmed(j)+sdev(j),
        !     $      xmed(j)+2.*sdev(j)
        write(29, 110)xage(j), xtem(ks - nperc), xtem(j2)
    10    continue
    do 12 j = nsteps + 1, 1, -1
        !          write(28,110)xage(j),xmed1(j),xmed(j)-sdev(j),
        !     $      xmed(j)-2.*sdev(j)
        write(29, 110)xage(j), xmed2(j), xmed1(j)
    12    continue
    !          write(28,110)xage(1),xmed0,xmed(1)+sdev(1),
    !     $      xmed(1)+2.*sdev(1)
    write(29, 110)xage(1), xmed02, xmed0
    !  WRITE INTERVAL OF CONFIDENCE
    14        return
    100    format(100(I5, 1x))
    110   format(8f8.2)
end

! FUNCTION TEMFC (CALCULATES THE TEMPERATURE AT ANY GIVEN AGE)
double precision function temfc(time, temp, ni, x)
    implicit double precision (a-h, o-z)
    dimension time(ni), temp(ni)
    temfc = 0.d0
    if(x.eq.time(1))then
        temfc = temp(1)
        return
    endif
    do 10 j = ni - 1, 1, -1
        if(time(j).gt.x.and.time(j + 1).le.x)then
            slope = (temp(j + 1) - temp(j)) / (time(j + 1) - time(j))
            temfc = temp(j) + slope * (x - time(j))
            return
        endif
    10    continue
end

!  SUBROUTINE MDIAN1
!  Given an array X of N numbers, returns their median value XMED. The
!  array X is modified and returned sorted into ascending order
!  (Num. Recipes, p.460)
subroutine mdian1(x, n, xmed)
    implicit double precision (a-h, o-z)
    dimension x(n)
    call sort(n, x)
    n2 = n / 2
    if(2 * n2.eq.n)then
        xmed = 0.5 * (x(n2) + x(n2 + 1))
    else
        xmed = x(n2 + 1)
    endif
    return
end

!  SUBROUTINE SORT
!  Sorts an array RA of length N into ascending numerical order using
!  the Heapsort algorithm. N is input; RA is replaced on output by its
!  sorted rearrangement. (Num. Recipes p.231)
subroutine sort(n, ra)
    implicit double precision (a-h, o-z)
    dimension ra(n)
    l = n / 2 + 1
    ir = n
    10    continue
    if(l.gt.1)then
        l = l - 1
        rra = ra(l)
    else
        rra = ra(ir)
        ra(ir) = ra(1)
        ir = ir - 1
        if(ir.eq.1)then
            ra(1) = rra
            return
        endif
    endif
    i = l
    j = l + l
    20      if(j.le.ir)then
        if(j.lt.ir)then
            if(ra(j).lt.ra(j + 1))j = j + 1
        endif
        if(rra.lt.ra(j))then
            ra(i) = ra(j)
            i = j
            j = j + j
        else
            j = ir + 1
        endif
        go to 20
    endif
    ra(i) = rra
    go to 10
end

!  SUBROUTINE MOMENT
!  Given an array of DATA of length N, this routine returns its mean AVE
!  average deviation ADEV, standard deviation SDEV,  variance VAR, skewness
!  SKEW, and kurtosis CURT. (Num. Recipes. p.458)
subroutine moment(data, n, ave, adev, sdev, var, skew, curt)
    implicit double precision (a-h, o-z)
    dimension data(n)
    if(n.le.1)pause 'n must be at least 2'
    s = 0.
    do 11 j = 1, n
        s = s + data(j)
    11    continue
    ave = s / n
    adev = 0.
    var = 0.
    skew = 0.
    curt = 0.
    do 12 j = 1, n
        s = data(j) - ave
        adev = adev + abs(s)
        p = s * s
        var = var + p
        p = p * s
        skew = skew + p
        p = p * s
        curt = curt + p
    12    continue
    adev = adev / n
    var = var / (n - 1)
    sdev = sqrt(var)
    if(var.ne.0.)then
        skew = skew / (n * sdev**3)
        curt = curt / (n * var**2) - 3.
    else
        print *, 'zero variance: no skew or kurtosis'
        skew = 0.
        curt = 0.
    endif
    return
end

!	BINOM.F
! ***********************************************************************
subroutine binom(n, j1, j2)
    implicit double precision (a-h, o-z)
    if(n / 2 * 2 - n.eq.0)then
        nmed = n / 2
    else
        nmed = (n + 1) / 2
    endif
    !  SEARCH FOR F(M)=0.5
    do 10 j = nmed, 1, -1
        faux = f
        f = 0.
        j1 = j
        do 12 k = 1, j
            f = bico(n, k) * 0.5**n + f
        12       continue
        if(f.le.0.05)then
            if(dabs(faux - .05).lt.dabs(f - .05))j1 = j1 + 1
            goto 14
        endif
    10    continue
    14    do 20 j = nmed + 1, n
        faux = f
        f = 0.
        j2 = j
        do 22 k = 1, j
            f = bico(n, k) * 0.5**n + f
        22       continue
        if(f.ge.0.95)then
            if(dabs(faux - .95).lt.dabs(f - .95))j2 = j2 - 1
            return
        endif
    20    continue
    stop 'error on binom sub'
end

double precision function bico(n, k)
    implicit double precision (a-h, o-z)
    bico = dnint(dexp(factln(n) - factln(k) - factln(n - k)))
    return
end

double precision function factln(n)
    implicit double precision (a-h, o-z)
    dimension a(100)
    data a/100*-1.d0/
    if (n.lt.0) pause 'negative factorial'
    y = 1.d0 + n
    if (n.le.99) then
        if (a(n + 1).lt.0.) a(n + 1) = gammln(y)
        factln = a(n + 1)
    else
        factln = gammln(y)
    endif
    return
end

function gammln(xx)
    implicit double precision (a-h, o-z)
    real*8 cof(6), stp, half, one, fpf, x, tmp, ser
    data cof, stp/76.18009173d0, -86.50532033d0, 24.01409822d0, &
            -1.231739516d0, .120858003d-2, -.536382d-5, 2.50662827465d0/
    data half, one, fpf/0.5d0, 1.0d0, 5.5d0/
    x = xx - one
    tmp = x + fpf
    tmp = (x + half) * dlog(tmp) - tmp
    ser = one
    do 11 j = 1, 6
        x = x + one
        ser = ser + cof(j) / x
    11    continue
    gammln = tmp + dlog(stp * ser)
    return
end

max_stars_repo_name: ASUPychron/pychron

max_issues_repo_name: ASUPychron/pychron

max_forks_repo_name: UIllinoisHALPychron/pychron


text:
! =====================================================================
! subroutine GRDR
!
! This subroutine generates the radial coordinate grid points.
! ---------------------------------------------------------------------
      subroutine GRDR(nR,rA,rB,dR,R,RH)
      implicit none

      integer :: nR
      double precision :: rA
      double precision :: rB
      double precision :: dR
      double precision :: R(nR)
      double precision :: RH(nR+1)

      integer :: i

      do i = 1, nR
         R(i) = rA + dfloat(i)*dR
         RH(i) = R(i) - 0.5d0*dR
      enddo
      RH(nR+1) = rB - 0.5d0*dR

      return
      end
! =====================================================================

max_stars_repo_name: mkandes/itpprp

max_issues_repo_name: mkandes/itpprp

max_forks_repo_name: mkandes/itpprp


text:
Citrus Circuits is a group of students at Davis High School who build robots.  Our most visible project around town is their RoboChair, This years robot is named Lemon Squeeze (Juice Master 1678).  Our main sponsor is Schilling Robotics.

For more robots in town, see the FIRST Robotics Competition.

Do you want to build a robot?! Just wait a year after you see this and hope that Jan (now retired) sponsors the robotics club. hahahahaha classic  Users/StevenDaubert

20070124 10:36:34 nbsp we are currently building a robot for FIRST, but Jan was the club supporter, now it is Steve Harvy Users/DevinCastellucci

20070128 17:57:55 nbsp Yeah I figured, the robotics club was always a standing joke in 03 with the tag line of Do you want to build a robot they finally got it formed and competing in what, 04 or 05? Users/StevenDaubert

20070128 18:42:06 nbsp it was offically formed in 2004, and the club competed in the 2005 FIRST Robotics Competition, we meet thursdays after school Users/DevinCastellucci

20071006 21:41:25 nbsp Hey guys, i need something to put on the flyers, got any ideas? Users/JessicaKobold

20071006 21:43:24 nbsp And no sarcasm unless youre actually going to contribute. I would love to here it, but not if I end up getting nothing done because of it. Users/JessicaKobold

20071006 21:44:29 nbsp by the way, Origami GoKart? Users/JessicaKobold

20071006 22:09:56 nbsp Also done with posters for booth and other such fundraising phenomena (literally), all i need is something to put on the flyer, so please try to come up with something. Thanks Users/JessicaKobold

20071006 22:12:46 nbsp And, yes i noticed the typo on one of my earlier messages, i am sorry, apologetic, and remorseful, so there. Users/JessicaKobold

20090402 08:23:48 nbsp JK microsystems might be willing to sponsor you and give you some hardware.  If you are interested, send me an email, jstewart@jkmicro.com Users/JimStewart


text:
!===============================================================================
! One of Andy Nowacki's Fortran utility modules for dealing with seismic
! anisotropy and other problems.
!
! Andy Nowacki <andy.nowacki@bristol.ac.uk>
!
! See the file LICENCE for licence details.
!===============================================================================
program gcp_from_point_azi_prog
!===============================================================================
! Wrapper program for sg_gcp_from_point_azi.  Given a point on a sphere and an
! azimuth through that point, return the pole to the great circle along that
! azimuth.  Reads lon, lat and azi from stdin, all in degrees

   use spherical_geometry, only: sg_gcp_from_point_azi

   implicit none

   real(8) :: lon,lat,azi,long,latg
   integer :: iostat

   if (command_argument_count() /= 0) call usage

   iostat = 0
   do while (iostat == 0)
      read(*,*,iostat=iostat) lon,lat,azi
      if (iostat > 0) then
         write(0,'(a)') 'gcp_from_point_azi: Error: problem reading lon,lat,azi (degrees) from stdin'
         stop
      endif
      if (iostat < 0) exit
      call sg_gcp_from_point_azi(lon,lat,azi,long,latg,degrees=.true.)
      write(*,'(f11.6,1x,f10.6)') long,latg
   enddo

contains
   !============================================================================
   subroutine usage()
   !============================================================================
      implicit none
      write(0,'(a)') 'Usage: gcp_from_point_azi < [lon] [lat] [azi]',&
                     '   Reads point on sphere and azimuth from stdin (degrees)',&
                     '   Writes lon and lat of pole to great circle defined by the input to stdout.'
      stop
   end subroutine usage
   !----------------------------------------------------------------------------

end program gcp_from_point_azi_prog
!-------------------------------------------------------------------------------

max_stars_repo_name: JackWalpole/seismo-fortran-fork

max_issues_repo_name: JackWalpole/seismo-fortran-fork

max_forks_repo_name: JackWalpole/seismo-fortran-fork