text:
Youre in Davis. Its the weekend. OR maybe its not. But regardless youre bored and youre tired and want to get out of Davis.  But where can you go and how long will it take to get there?  The destinations are broken down by approximate driving time if you drive without stopping (or traffic slowdown) and at the posted speed limit.

Also see Outdoor Activities and follow some links or check out Hiking and Backpacking for lots of wilderness destination leads. The Northern California page contains a lot of information also. UC Davis UCD students interested in archaeology and outdoor activities might consider the field schools offered by the Anthropology department.  And dont forget about Beaches, although some particular ones are listed below.

If you require something more linear, see Interstate 80 I80 cartrips (Idiot proof trips).

For trips inside California, you can check road conditions on the http://www.dot.ca.gov/ Department of Transportation Website or visit the website http://www.trafficreport.com/ TrafficReport.com to obtain traffic information for areas in California.

 < 1 Hour Away
  Sacramento (wiki:sacramento:Front_Page Sacramento Wiki), Woodland (wiki:woodland:Front_Page Woodland Wiki), Winters (wiki:winters:Front_Page Winters Wiki) and Dixon (wiki:dixon:Front_Page Dixon Wiki)
  wiki:sacramento:Sutter%27s_Fort Sutters Fort SHP
  wiki:woodland:museums Museums in Woodland CA, wiki:sacramento:museums Museums in Sacramento CA
  A couple Indian Casinos (See Gambling.)
  http://www.parks.ca.gov/?page_id500 Folsom Lake (longer if you Bicycle Rides bike)If the water levels are especially low, you can see the remains of a gold rush town called http://files.myfolsom.com/photosmormonisland.shtml Mormon Island.
  Davis Paintball Center
  http://www.dirosapreserve.org di Rosa Preserve (46 mi., 55 min.) art and sculpture garden
  Jepson Prairie Reserve
  Stevenson Bridge
  Stebbins Cold Canyon Reserve
  Yolo County Grasslands Regional Park  archery, birdwatching, R/C glider soaring
  Jelly Belly Factory Tour with samples and goodies.
  http://www.budweisertours.com/docs/Fairfield/hours.htm Budweiser Factory Tour with free samples and free soda for the kids!
  The Baxter House
  Yolo County Recreation
  http://www.sixflags.com/discoverykingdom/index.aspx Six Flags Discovery Kingdom (formerly Marine World) Roller Coaster Theme Park in Vallejo.
  Woodland Corn Maze  Seasonal, in September and October.

 1 Hour Away
  Napa Napa Valley  (wine Wine tasting, sightseeing, gourmet food, Hot Air Balloons Hot Air Balloon Rides)
  Sonoma County Sonoma Valley  (Wine tasting, sightseeing, gourmet food)
  Various Indian Casinos (See Gambling.)
  Lake Berryessa (Fishing, water sports, sightseeing) and the Morning Glory SpillwayLow water levels means land access to some normally unreachable islands.
  Cache Creek (Hiking and Backpacking hiking, Camping camping, rafting, nature)
  Berkeley (Shopping, mischief, excellent ethnic food)
     Telegraph Avenue
     College Avenue
     Berkeley Outlets
     Zacharys Pizza (An adventure for your taste buds)

Its a shame UC Berkeley doesnt have a wiki for Berkeley because there are so many restaurants there which are ridiculously good: IN THE MOOD FOR... Indian? Naann Curry or TKs. Hot Dogs? Top Dog or Fat Dog. Sandwiches? IBs or Bongo Burger. Latin? Cafe de la Paz. (its a bit expensive, but good)  Emily Tung
    They used to have one, but it went down. Users/JaimeRaba jr
     Technically they still do, but no one does anything with it Users/MichelleAccurso
   Davis, California does indeed have Indian, Hot Dogs, Sandwiches, and Latin cuisine as well.  Why waste the gas/train/bike lube money? users/JasonDunne Jason Dunne

 1.5 Hours Away
  San Francisco (Gourmet food, museums, art galleries, clubbing, sightseeing, shopping)
     AT&T Park
     North Beach
     Chinatown
     Fishermans Wharf
     Castro District
     Haight and Ashbury
  http://www.applehill.com/ Apple Hill (located up HWY 50 past Placerville; offers a scenic drive and during the fall farms are open to guests to pick apples, pick pumpkins, eat apple pie, eat caramel apples, and drink apple cider)
  Lovers Leap (best climbing in the area; located on highway 50 before Lake Tahoe)
  http://www.fs.fed.us/r5/eldorado/recreation/wild/deso/ Desolation Wilderness: 63,960 acres of subalpine and alpine forest, granite peaks, and glaciallyformed valleys and lakes in the El Dorado National Forest.
  http://www.wilburhotsprings.com Wilbur Hot Springs A sanctuary for the self. Established 1865. 1h 32m away.
  Muir Beach and http://www.nps.gov/muwo/ Muir Woods (beach, hiking, fresh air, peace and quiet, history of parks, redwoods).  Just north of San Francisco.
  Mount Saint Helena  (Hiking and Backpacking Hiking and Climbing climbing after nice drive with great views.)
  Clear Lake  (Hiking and Backpacking hiking ,fishing fishing, water sports, and wine). The town of Clearlake Oaks offers access to the lake and is 90 miles from Davis.
  http://www.waterfallswest.com/waterfall.php?id210 University Falls Located near the town of Georgetown in El Dorado National Forest.  Its about a 45 minute hike from the road where you park.  There are 4 waterfalls, and you can slide down the first threenatural granite watersides!
  Yuba River Depending on the area, you could enjoy a calming river float or spend the day scrambling over boulders to reach swimming holes.

 2 Hours Away
  Reno, NV  (Gambling, Prostitution technically not in town but nearby)
  Lake Tahoe, CA (Water sports, Snow Sports snow sports, Hiking and Backpacking hiking)
  Lake Tahoe, NV (Gambling, Water sports, Snow Sports snow sports, Hiking and Backpacking hiking)
  Santa Cruz, CA  (Shopping, sightseeing, beaches, Hiking and Backpacking hiking, general awesomeness)
      Beach Boardwalk
      UC Santa Cruz
  Chico, CA, including wiki:Chico Sierra Nevada Brewing Company
  http://www1.cedarfair.com/greatamerica/index.cfm Great America amusement park
  http://www.hersheys.com/discover/oakdale.asp Hershey Plant in Oakdale. 99 south to CA120. They no longer do factory tours.
  http://www.harbin.org Harbin Hot Springs Soak, or go hiking, socialize, or attend workshops. 2h 10m away.

 2.5 Hours Away
  Monterey, CA (aquarium, ocean, beach, Hiking and Backpacking hiking)

 3 Hours Away
  Mendocino, CA (Hiking and Backpacking hiking, wine, ocean, hippies)

  Lake Shasta (Houseboats houseboating)
  Lassen Volcanic National Park Lassen National Park (wiki:WikiPedia:Lassen_Volcanic_National_Park Wikipedia entry) (hiking/backpacking etc.)
   Ive once gotten there in under 2 hours. I love I5  Users/JulienBiewerElstob

 4.5 Hours Away
  Yosemite National Park (http://www.nps.gov/yose/ website) (hiking/backpacking, Climbing climbing, sightseeing, nature, fresh air)
  Mono Lake (wiki:WikiPedia:Mono_lake Wikipedia entry) (sightseeing, nature, unique ecosystem)
  San Luis Obispo, CA (Ocean, shopping, yuppies)
  Eureka, CA (Hiking and Backpacking hiking, the Lost Coast)
  Black Rock City, NV for Burning Man (wiki:WikiPedia:Burning_Man Wikipedia entry)

 5 Hours Away
  Valencia, CA (Six Flags Magic Mountain!)
  Arcata, CA (http://www.kineticsculpturerace.org/ Kinetic Sculpture Race, http://www.humboldt.edu/ Humboldt State University)

 5.5 Hours Away
  Mammoth Lakes, CA (Snow Sports snow sports, water sports, mountain biking, Hiking and Backpacking hiking and backpacking, sightseeing, nature (volcanoes!)
  Los Angeles, CA  (Shopping, sightseeing, clubbing, smog, beach)
    Venice Beach
    Hollywood
    Garment District
  http://www.nps.gov/seki/ Sequoia & Kings Canyon National Parks (Hiking and Backpacking hiking and backpacking, camping, sightseeing, nature, fresh air)

 6.5 Hours Away
   Anaheim, CA (Disneyland!)
   Irvine, CA
      UC Irvine
      Irvine Spectrum
   http://www.nps.gov/redw/index.htm Redwood National and State Parks

 8 Hours Away
  San Diego, CA  (Shopping, sightseeing, beach)
    Gas Lamp District
  Tijuana, B.C.  Mexico (Practice Spanish, shopping, sightseeing, smog, beach)
  Death Valley National Park (http://www.nps.gov/deva/ website) (hiking/backpacking, camping, sightseeing, nature, fresh air)
     Note:  Southern Route (by way of Bakersfield) is safer and faster, but less fun. Northern Route (by way of Reno) is more scenic.  Twolane roads that go through towns make traveltime extremely unpredictable for both routes.
  Joshua Tree National Park (http://www.nps.gov/jotr/ website) (Hiking and Backpacking hiking and backpacking, Climbing climbing, camping, sightseeing, nature, fresh air)
  Las Vegas, NV  (Gambling, sightseeing, food, clubbing)
    The Strip (Las Vegas Blvd.)
    Fremont Street  (Classic Downtown casinos)
  Portland, OR  (Shopping, sightseeing, food, fresh air)
    Powells Bookstore

 As far away as possible
 District de Crozet

 Adventuring Tips

 If driving somewhere far like Las Vegas, consider car rental renting a car, especially if you have a gas guzzler.  With rental car rates as low as $1020/day, sometimes this is a net cost savings when you can get a car with 3545mpg.  Plus you spare your car the added mileage.
   Of course, if youre not at least 25, several companies will not rent to you, while others will charge you an outrageous fee.
 Be sure to bring one or many coadventurers.  If none of your friends share your wanderlust, perhaps the Ride Board can hook you up.
  AAA is your friend. In addition to towing, lockout and fuel delivery services, they also can make road trip calamities better. If your car breaks down, many AAA memberships will pay for your hotel and food costs while your vehicle is being fixed or, alternatively, your return transportation (either air, bus, train or rental car). Keep in mind, however, that the AAA frequently advocates against bicyclefriendly laws and infrastructure when there is no benefit to automobiles. A more bicyclefriendly alternative is http://www.betterworldclub.com/ The Better World Club.
 Dont forget your cell phone charger.  And dont expect your cell phone to work in the middle of Nevada.
 The Davis Travelaires lead trips that can help get you comfortable and excited about travel.
 Note: If you are driving through LA on the I5, remember that the http://en.wikipedia.org/wiki/Grapevine%2C_California Grapevine can be tricky for older cars to navigate smoothly, and that downtown LA traffic is not your friend.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


!!$ A C-program for MT19937: Real number version
!!$   genrand() generates one pseudorandom real number (double)
!!$ which is uniformly distributed on [0,1]-interval, for each
!!$ call. sgenrand(seed) set initial values to the working area
!!$ of 624 words. Before genrand(), sgenrand(seed) must be
!!$ called once. (seed is any 32-bit integer except for 0).
!!$ Integer generator is obtained by modifying two lines.
!!$   Coded by Takuji Nishimura, considering the suggestions by
!!$ Topher Cooper and Marc Rieffel in July-Aug. 1997.
!!$
!!$ This library is free software; you can redistribute it and/or
!!$ modify it under the terms of the GNU Library General Public
!!$ License as published by the Free Software Foundation; either
!!$ version 2 of the License, or (at your option) any later
!!$ version.
!!$ This library is distributed in the hope that it will be useful,
!!$ but WITHOUT ANY WARRANTY; without even the implied warranty of
!!$ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
!!$ See the GNU Library General Public License for more details.
!!$ You should have received a copy of the GNU Library General
!!$ Public License along with this library; if not, write to the
!!$ Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
!!$ 02111-1307  USA
!!$
!!$ Copyright (C) 1997 Makoto Matsumoto and Takuji Nishimura.
!!$ When you use this, send an email to: matumoto@math.keio.ac.jp
!!$ with an appropriate reference to your work.
!!$
!!$***********************************************************************
!!$ Fortran translation by Hiroshi Takano.  Jan. 13, 1999.
!!$
!!$ This program uses the following non-standard intrinsics.
!!$   ishft(i,n): If n>0, shifts bits in i by n positions to left.
!!$               If n<0, shifts bits in i by n positions to right.
!!$   iand (i,j): Performs logical AND on corresponding bits of i and j.
!!$   ior  (i,j): Performs inclusive OR on corresponding bits of i and j.
!!$   ieor (i,j): Performs exclusive OR on corresponding bits of i and j.
!!$
!!$************************************************************************
!!$ Fortran 95 translation and modularized to use in gyrokinetics project
!!$ by R. Numata, June 2, 2010.
!!$  * removed statement functions (TSHFTU,TSHFTS,TSHFTT,TSHFTL)
!!$  * uses [0,1) interval
!!$  * the above bit manipulation functions became standard.
!!$  * definition of UMASK is corrected (the original value cannot be
!!$    represented as a kind=4 integer.
!!$***********************************************************************

module mt19937

  implicit none

  private
  public :: sgrnd, grnd

  integer, parameter :: default_seed=4357

  ! Period parameters
  integer, parameter :: N = 624
  integer, parameter :: M = 397
  integer, parameter :: MATA = -1727483681 ! constant vector a
  integer, parameter :: LMASK = 2147483647 ! least significant r bits
  integer, parameter :: UMASK = -LMASK-1 ! most significant w-r bits

  ! Tempering parameters
  integer, parameter :: TMASKB= -1658038656
  integer, parameter :: TMASKC= -272236544

  integer, save :: mt(0:N-1) ! the array for the state vector
  integer, save :: mti=N+1 ! mti==N+1 means mt[N] is not initialized
  integer, save :: mag01(0:1)=(/ 0, MATA /) ! mag01(x) = x * MATA for x=0,1

contains

  subroutine sgrnd(seed)
    implicit none
    integer, intent(in) :: seed

!!$      setting initial seeds to mt[N] using
!!$      the generator Line 25 of Table 1 in
!!$      [KNUTH 1981, The Art of Computer Programming
!!$         Vol. 2 (2nd Ed.), pp102]

    mt(0)= iand(seed,-1)
    do mti=1,N-1
       mt(mti) = iand(69069 * mt(mti-1),-1)
    end do

    return
  end subroutine sgrnd

  function grnd ()
    implicit none
    real :: grnd
    real, parameter :: pow=4294967296.0 ! 2**32
    real, parameter :: div=1./pow       ! devided by 2**32 [0,1)-real-interval
    ! real, parameter :: div=1./(pow-1.)  ! devided by 2**32-1 [0,1]-real-interval
    integer :: y, kk

    if(mti >= N) then ! generate N words at one time
       if(mti == N+1) then ! if sgrnd() has not been called,
          call sgrnd(default_seed) ! a default initial seed is used
       endif

       do kk=0,N-M-1
          y=ior(iand(mt(kk),UMASK),iand(mt(kk+1),LMASK))
          mt(kk)=ieor(ieor(mt(kk+M),ishft(y,-1)),mag01(iand(y,1)))
       end do

       do kk=N-M,N-2
          y=ior(iand(mt(kk),UMASK),iand(mt(kk+1),LMASK))
          mt(kk)=ieor(ieor(mt(kk+(M-N)),ishft(y,-1)),mag01(iand(y,1)))
       end do

       y=ior(iand(mt(N-1),UMASK),iand(mt(0),LMASK))
       mt(N-1)=ieor(ieor(mt(M-1),ishft(y,-1)),mag01(iand(y,1)))
       mti = 0
    endif

    y=mt(mti)
    mti=mti+1
    y=ieor(y,ishft(y,-11))
    y=ieor(y,iand(ishft(y,7),TMASKB))
    y=ieor(y,iand(ishft(y,15),TMASKC))
    y=ieor(y,ishft(y,-18))

    grnd = real(y)
    if(grnd < 0.) grnd=grnd + pow
    grnd = grnd * div

    return
  end function grnd

end module mt19937

!!$ this main outputs first 1000 generated numbers
!!$program main
!!$  use mt19937, only: grnd
!!$  ! use mt19937, only: sgrnd
!!$  implicit none
!!$  integer, parameter :: no=1000
!!$  real :: r(0:7)
!!$  integer :: j,k
!!$  ! call sgrnd(4357) ! any nonzero integer can be used as a seed
!!$
!!$  do j=0,no-1
!!$     r(mod(j,8))=grnd()
!!$     if(mod(j,8) == 7) then
!!$        write(*,'(8(f8.6,'' ''))') (r(k),k=0,7)
!!$     else if(j == no-1) then
!!$        write(*,'(8(f8.6,'' ''))') (r(k),k=0,mod(no-1,8))
!!$     endif
!!$  end do
!!$
!!$  stop
!!$end program main


text:
Describe Users/KyleLanderhoffis here.

WTF? meaning, wheres the fun

20071001 19:26:47 nbsp Welcome to the Wiki Users/RyanJames

ext: f


text:
! imap.f90
!
! Example that shows how to examine the INBOX of an IMAP mailbox, using
! libcurl. Based on the following C implementation:
!
!     https://curl.haxx.se/libcurl/c/imap-examine.html
!
! Author:  Philipp Engel
! Licence: ISC
module callback_imap
    use, intrinsic :: iso_c_binding
    use :: curl, only: c_f_str_ptr
    implicit none
    private
    public :: write_callback
contains
    ! static size_t callback(void *ptr, size_t size, size_t nmemb, void *data)
    function write_callback(ptr, sze, nmemb, data) bind(c)
        !! Callback function for `CURLOPT_WRITEFUNCTION` that simply outputs
        !! the IMAP messages.
        type(c_ptr),            intent(in), value :: ptr
        integer(kind=c_size_t), intent(in), value :: sze
        integer(kind=c_size_t), intent(in), value :: nmemb
        type(c_ptr),            intent(in), value :: data
        integer(kind=c_size_t)                    :: write_callback
        character(len=:), allocatable             :: str

        write_callback = int(0, kind=c_size_t)

        write (*, '(a)') '---'
        write (*, '(a)') '--- CALLBACK INVOKED'
        write (*, '(a)') '---'

        if (.not. c_associated(ptr)) return
        call c_f_str_ptr(ptr, str, nmemb)
        write (*, '(a)', advance='no') str
        write_callback = nmemb
    end function write_callback
end module callback_imap

program main
    use, intrinsic :: iso_c_binding
    use :: curl
    use :: callback_imap
    implicit none

    character(len=*), parameter :: URL      = 'imaps://example.com' ! IMAP server (SSL).
    character(len=*), parameter :: USERNAME = 'user'                ! IMAP user name.
    character(len=*), parameter :: PASSWORD = 'secret'              ! IMAP password.

    type(c_ptr) :: curl_ptr
    integer     :: rc

    curl_ptr = curl_easy_init()

    if (c_associated(curl_ptr)) then
        ! Set curl options.
        rc = curl_easy_setopt(curl_ptr, CURLOPT_URL,            URL // c_null_char)
        rc = curl_easy_setopt(curl_ptr, CURLOPT_USERNAME,       USERNAME // c_null_char)
        rc = curl_easy_setopt(curl_ptr, CURLOPT_PASSWORD,       PASSWORD // c_null_char)
        rc = curl_easy_setopt(curl_ptr, CURLOPT_CUSTOMREQUEST,  'EXAMINE INBOX' // c_null_char)
        rc = curl_easy_setopt(curl_ptr, CURLOPT_SSL_VERIFYPEER, int(1, kind=8))
        rc = curl_easy_setopt(curl_ptr, CURLOPT_SSL_VERIFYHOST, int(1, kind=8))
        rc = curl_easy_setopt(curl_ptr, CURLOPT_WRITEFUNCTION,  c_funloc(write_callback))
        rc = curl_easy_setopt(curl_ptr, CURLOPT_VERBOSE,        int(1, kind=8))

        ! Perform request..
        if (curl_easy_perform(curl_ptr) /= CURLE_OK) then
            print '(a)', 'Error: curl_easy_perform() failed'
        end if

        call curl_easy_cleanup(curl_ptr)
    end if
end program main

max_stars_repo_name: awvwgk/fortran-curl

max_issues_repo_name: awvwgk/fortran-curl

max_forks_repo_name: interkosmos/fortran-curl


text:
c  Collection of fftpack obtained from netlib.org needed to computer
c   the inverse fft such that the result is 'real'
c
      SUBROUTINE ZFFTI (N,WSAVE)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION       WSAVE(1)
      IF (N .EQ. 1) RETURN
      IW1 = N+N+1
      IW2 = IW1+N+N
      CALL CFFTI1 (N,WSAVE(IW1),WSAVE(IW2))
      RETURN
      END
      SUBROUTINE DZFFTI (N,WSAVE)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION       WSAVE(1)
      IF (N .EQ. 1) RETURN
      CALL EZFFT1 (N,WSAVE(2*N+1),WSAVE(3*N+1))
      RETURN
      END
      SUBROUTINE CFFTI1 (N,WA,IFAC)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION       WA(*)      ,IFAC(*)    ,NTRYH(4)
      DATA NTRYH(1),NTRYH(2),NTRYH(3),NTRYH(4)/3,4,2,5/
      NL = N
      NF = 0
      J = 0
  101 J = J+1
      IF (J-4) 102,102,103
  102 NTRY = NTRYH(J)
      GO TO 104
  103 NTRY = NTRY+2
  104 NQ = NL/NTRY
      NR = NL-NTRY*NQ
      IF (NR) 101,105,101
  105 NF = NF+1
      IFAC(NF+2) = NTRY
      NL = NQ
      IF (NTRY .NE. 2) GO TO 107
      IF (NF .EQ. 1) GO TO 107
      DO 106 I=2,NF
         IB = NF-I+2
         IFAC(IB+2) = IFAC(IB+1)
  106 CONTINUE
      IFAC(3) = 2
  107 IF (NL .NE. 1) GO TO 104
      IFAC(1) = N
      IFAC(2) = NF
      TPI =  6.28318530717958647692D0
      ARGH = TPI/FLOAT(N)
      I = 2
      L1 = 1
      DO 110 K1=1,NF
         IP = IFAC(K1+2)
         LD = 0
         L2 = L1*IP
         IDO = N/L2
         IDOT = IDO+IDO+2
         IPM = IP-1
         DO 109 J=1,IPM
            I1 = I
            WA(I-1) = 1.0D0
            WA(I) = 0.0D0
            LD = LD+L1
            FI = 0.0D0
            ARGLD = FLOAT(LD)*ARGH
            DO 108 II=4,IDOT,2
               I = I+2
               FI = FI+1.D0
               ARG = FI*ARGLD
               WA(I-1) = COS(ARG)
               WA(I) = SIN(ARG)
  108       CONTINUE
            IF (IP .LE. 5) GO TO 109
            WA(I1-1) = WA(I-1)
            WA(I1) = WA(I)
  109    CONTINUE
         L1 = L2
  110 CONTINUE
      RETURN
      END
      SUBROUTINE EZFFT1 (N,WA,IFAC)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION       WA(*)      ,IFAC(*)    ,NTRYH(4)
      DATA NTRYH(1),NTRYH(2),NTRYH(3),NTRYH(4)/4,2,3,5/
     1    ,TPI/6.28318530717958647692D0/
      NL = N
      NF = 0
      J = 0
  101 J = J+1
      IF (J-4) 102,102,103
  102 NTRY = NTRYH(J)
      GO TO 104
  103 NTRY = NTRY+2
  104 NQ = NL/NTRY
      NR = NL-NTRY*NQ
      IF (NR) 101,105,101
  105 NF = NF+1
      IFAC(NF+2) = NTRY
      NL = NQ
      IF (NTRY .NE. 2) GO TO 107
      IF (NF .EQ. 1) GO TO 107
      DO 106 I=2,NF
         IB = NF-I+2
         IFAC(IB+2) = IFAC(IB+1)
  106 CONTINUE
      IFAC(3) = 2
  107 IF (NL .NE. 1) GO TO 104
      IFAC(1) = N
      IFAC(2) = NF
      ARGH = TPI/FLOAT(N)
      IS = 0
      NFM1 = NF-1
      L1 = 1
      IF (NFM1 .EQ. 0) RETURN
      DO 111 K1=1,NFM1
         IP = IFAC(K1+2)
         L2 = L1*IP
         IDO = N/L2
         IPM = IP-1
         ARG1 = FLOAT(L1)*ARGH
         CH1 = 1.0D0
         SH1 = 0.0D0
         DCH1 = COS(ARG1)
         DSH1 = SIN(ARG1)
         DO 110 J=1,IPM
            CH1H = DCH1*CH1-DSH1*SH1
            SH1 = DCH1*SH1+DSH1*CH1
            CH1 = CH1H
            I = IS+2
            WA(I-1) = CH1
            WA(I) = SH1
            IF (IDO .LT. 5) GO TO 109
            DO 108 II=5,IDO,2
               I = I+2
               WA(I-1) = CH1*WA(I-3)-SH1*WA(I-2)
               WA(I) = CH1*WA(I-2)+SH1*WA(I-3)
  108       CONTINUE
  109       IS = IS+IDO
  110    CONTINUE
         L1 = L2
  111 CONTINUE
      RETURN
      END
      SUBROUTINE DZFFTB (N,R,AZERO,A,B,WSAVE)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION       R(*)       ,A(*)       ,B(*)       ,WSAVE(*)
      IF (N-2) 101,102,103
  101 R(1) = AZERO
      RETURN
  102 R(1) = AZERO+A(1)
      R(2) = AZERO-A(1)
      RETURN
  103 NS2 = (N-1)/2
      DO 104 I=1,NS2
         R(2*I) = 0.5D0*A(I)
         R(2*I+1) = -0.5D0*B(I)
  104 CONTINUE
      R(1) = AZERO
      IF (MOD(N,2) .EQ. 0) R(N) = A(NS2+1)
      CALL DFFTB (N,R,WSAVE(N+1))
      RETURN
      END
      SUBROUTINE DFFTB (N,R,WSAVE)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION       R(1)       ,WSAVE(1)
      IF (N .EQ. 1) RETURN
      CALL RFFTB1 (N,R,WSAVE,WSAVE(N+1),WSAVE(2*N+1))
      RETURN
      END
      SUBROUTINE RFFTB1 (N,C,CH,WA,IFAC)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION       CH(*)      ,C(*)       ,WA(*)      ,IFAC(*)
      NF = IFAC(2)
      NA = 0
      L1 = 1
      IW = 1
      DO 116 K1=1,NF
         IP = IFAC(K1+2)
         L2 = IP*L1
         IDO = N/L2
         IDL1 = IDO*L1
         IF (IP .NE. 4) GO TO 103
         IX2 = IW+IDO
         IX3 = IX2+IDO
         IF (NA .NE. 0) GO TO 101
         CALL RADB4 (IDO,L1,C,CH,WA(IW),WA(IX2),WA(IX3))
         GO TO 102
  101    CALL RADB4 (IDO,L1,CH,C,WA(IW),WA(IX2),WA(IX3))
  102    NA = 1-NA
         GO TO 115
  103    IF (IP .NE. 2) GO TO 106
         IF (NA .NE. 0) GO TO 104
         CALL RADB2 (IDO,L1,C,CH,WA(IW))
         GO TO 105
  104    CALL RADB2 (IDO,L1,CH,C,WA(IW))
  105    NA = 1-NA
         GO TO 115
  106    IF (IP .NE. 3) GO TO 109
         IX2 = IW+IDO
         IF (NA .NE. 0) GO TO 107
         CALL RADB3 (IDO,L1,C,CH,WA(IW),WA(IX2))
         GO TO 108
  107    CALL RADB3 (IDO,L1,CH,C,WA(IW),WA(IX2))
  108    NA = 1-NA
         GO TO 115
  109    IF (IP .NE. 5) GO TO 112
         IX2 = IW+IDO
         IX3 = IX2+IDO
         IX4 = IX3+IDO
         IF (NA .NE. 0) GO TO 110
         CALL RADB5 (IDO,L1,C,CH,WA(IW),WA(IX2),WA(IX3),WA(IX4))
         GO TO 111
  110    CALL RADB5 (IDO,L1,CH,C,WA(IW),WA(IX2),WA(IX3),WA(IX4))
  111    NA = 1-NA
         GO TO 115
  112    IF (NA .NE. 0) GO TO 113
         CALL RADBG (IDO,IP,L1,IDL1,C,C,C,CH,CH,WA(IW))
         GO TO 114
  113    CALL RADBG (IDO,IP,L1,IDL1,CH,CH,CH,C,C,WA(IW))
  114    IF (IDO .EQ. 1) NA = 1-NA
  115    L1 = L2
         IW = IW+(IP-1)*IDO
  116 CONTINUE
      IF (NA .EQ. 0) RETURN
      DO 117 I=1,N
         C(I) = CH(I)
  117 CONTINUE
      RETURN
      END
      SUBROUTINE RADB2 (IDO,L1,CC,CH,WA1)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION       CC(IDO,2,L1)           ,CH(IDO,L1,2)           ,
     1                WA1(1)
      DO 101 K=1,L1
         CH(1,K,1) = CC(1,1,K)+CC(IDO,2,K)
         CH(1,K,2) = CC(1,1,K)-CC(IDO,2,K)
  101 CONTINUE
      IF (IDO-2) 107,105,102
  102 IDP2 = IDO+2
      DO 104 K=1,L1
         DO 103 I=3,IDO,2
            IC = IDP2-I
            CH(I-1,K,1) = CC(I-1,1,K)+CC(IC-1,2,K)
            TR2 = CC(I-1,1,K)-CC(IC-1,2,K)
            CH(I,K,1) = CC(I,1,K)-CC(IC,2,K)
            TI2 = CC(I,1,K)+CC(IC,2,K)
            CH(I-1,K,2) = WA1(I-2)*TR2-WA1(I-1)*TI2
            CH(I,K,2) = WA1(I-2)*TI2+WA1(I-1)*TR2
  103    CONTINUE
  104 CONTINUE
      IF (MOD(IDO,2) .EQ. 1) RETURN
  105 DO 106 K=1,L1
         CH(IDO,K,1) = CC(IDO,1,K)+CC(IDO,1,K)
         CH(IDO,K,2) = -(CC(1,2,K)+CC(1,2,K))
  106 CONTINUE
  107 RETURN
      END
      SUBROUTINE RADB3 (IDO,L1,CC,CH,WA1,WA2)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION       CC(IDO,3,L1)           ,CH(IDO,L1,3)           ,
     1                WA1(1)     ,WA2(1)
C     *** TAUI IS SQRT(3)/2 ***
      DATA TAUR,TAUI /-0.5D0,0.86602540378443864676D0/
      DO 101 K=1,L1
         TR2 = CC(IDO,2,K)+CC(IDO,2,K)
         CR2 = CC(1,1,K)+TAUR*TR2
         CH(1,K,1) = CC(1,1,K)+TR2
         CI3 = TAUI*(CC(1,3,K)+CC(1,3,K))
         CH(1,K,2) = CR2-CI3
         CH(1,K,3) = CR2+CI3
  101 CONTINUE
      IF (IDO .EQ. 1) RETURN
      IDP2 = IDO+2
      DO 103 K=1,L1
         DO 102 I=3,IDO,2
            IC = IDP2-I
            TR2 = CC(I-1,3,K)+CC(IC-1,2,K)
            CR2 = CC(I-1,1,K)+TAUR*TR2
            CH(I-1,K,1) = CC(I-1,1,K)+TR2
            TI2 = CC(I,3,K)-CC(IC,2,K)
            CI2 = CC(I,1,K)+TAUR*TI2
            CH(I,K,1) = CC(I,1,K)+TI2
            CR3 = TAUI*(CC(I-1,3,K)-CC(IC-1,2,K))
            CI3 = TAUI*(CC(I,3,K)+CC(IC,2,K))
            DR2 = CR2-CI3
            DR3 = CR2+CI3
            DI2 = CI2+CR3
            DI3 = CI2-CR3
            CH(I-1,K,2) = WA1(I-2)*DR2-WA1(I-1)*DI2
            CH(I,K,2) = WA1(I-2)*DI2+WA1(I-1)*DR2
            CH(I-1,K,3) = WA2(I-2)*DR3-WA2(I-1)*DI3
            CH(I,K,3) = WA2(I-2)*DI3+WA2(I-1)*DR3
  102    CONTINUE
  103 CONTINUE
      RETURN
      END
      SUBROUTINE RADB4 (IDO,L1,CC,CH,WA1,WA2,WA3)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION       CC(IDO,4,L1)           ,CH(IDO,L1,4)           ,
     1                WA1(1)     ,WA2(1)     ,WA3(1)
      DATA SQRT2 /1.41421356237309504880D0/
      DO 101 K=1,L1
         TR1 = CC(1,1,K)-CC(IDO,4,K)
         TR2 = CC(1,1,K)+CC(IDO,4,K)
         TR3 = CC(IDO,2,K)+CC(IDO,2,K)
         TR4 = CC(1,3,K)+CC(1,3,K)
         CH(1,K,1) = TR2+TR3
         CH(1,K,2) = TR1-TR4
         CH(1,K,3) = TR2-TR3
         CH(1,K,4) = TR1+TR4
  101 CONTINUE
      IF (IDO-2) 107,105,102
  102 IDP2 = IDO+2
      DO 104 K=1,L1
         DO 103 I=3,IDO,2
            IC = IDP2-I
            TI1 = CC(I,1,K)+CC(IC,4,K)
            TI2 = CC(I,1,K)-CC(IC,4,K)
            TI3 = CC(I,3,K)-CC(IC,2,K)
            TR4 = CC(I,3,K)+CC(IC,2,K)
            TR1 = CC(I-1,1,K)-CC(IC-1,4,K)
            TR2 = CC(I-1,1,K)+CC(IC-1,4,K)
            TI4 = CC(I-1,3,K)-CC(IC-1,2,K)
            TR3 = CC(I-1,3,K)+CC(IC-1,2,K)
            CH(I-1,K,1) = TR2+TR3
            CR3 = TR2-TR3
            CH(I,K,1) = TI2+TI3
            CI3 = TI2-TI3
            CR2 = TR1-TR4
            CR4 = TR1+TR4
            CI2 = TI1+TI4
            CI4 = TI1-TI4
            CH(I-1,K,2) = WA1(I-2)*CR2-WA1(I-1)*CI2
            CH(I,K,2) = WA1(I-2)*CI2+WA1(I-1)*CR2
            CH(I-1,K,3) = WA2(I-2)*CR3-WA2(I-1)*CI3
            CH(I,K,3) = WA2(I-2)*CI3+WA2(I-1)*CR3
            CH(I-1,K,4) = WA3(I-2)*CR4-WA3(I-1)*CI4
            CH(I,K,4) = WA3(I-2)*CI4+WA3(I-1)*CR4
  103    CONTINUE
  104 CONTINUE
      IF (MOD(IDO,2) .EQ. 1) RETURN
  105 CONTINUE
      DO 106 K=1,L1
         TI1 = CC(1,2,K)+CC(1,4,K)
         TI2 = CC(1,4,K)-CC(1,2,K)
         TR1 = CC(IDO,1,K)-CC(IDO,3,K)
         TR2 = CC(IDO,1,K)+CC(IDO,3,K)
         CH(IDO,K,1) = TR2+TR2
         CH(IDO,K,2) = SQRT2*(TR1-TI1)
         CH(IDO,K,3) = TI2+TI2
         CH(IDO,K,4) = -SQRT2*(TR1+TI1)
  106 CONTINUE
  107 RETURN
      END
      SUBROUTINE RADB5 (IDO,L1,CC,CH,WA1,WA2,WA3,WA4)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION       CC(IDO,5,L1)           ,CH(IDO,L1,5)           ,
     1                WA1(1)     ,WA2(1)     ,WA3(1)     ,WA4(1)
C     *** TR11=COS(2*PI/5), TI11=SIN(2*PI/5)
C     *** TR12=COS(4*PI/5), TI12=SIN(4*PI/5)
      DATA TR11,TI11,TR12,TI12 /0.3090169943749474241D0,
     +     0.95105651629515357212D0,
     +     -0.8090169943749474241D0,0.58778525229247312917D0/
      DO 101 K=1,L1
         TI5 = CC(1,3,K)+CC(1,3,K)
         TI4 = CC(1,5,K)+CC(1,5,K)
         TR2 = CC(IDO,2,K)+CC(IDO,2,K)
         TR3 = CC(IDO,4,K)+CC(IDO,4,K)
         CH(1,K,1) = CC(1,1,K)+TR2+TR3
         CR2 = CC(1,1,K)+TR11*TR2+TR12*TR3
         CR3 = CC(1,1,K)+TR12*TR2+TR11*TR3
         CI5 = TI11*TI5+TI12*TI4
         CI4 = TI12*TI5-TI11*TI4
         CH(1,K,2) = CR2-CI5
         CH(1,K,3) = CR3-CI4
         CH(1,K,4) = CR3+CI4
         CH(1,K,5) = CR2+CI5
  101 CONTINUE
      IF (IDO .EQ. 1) RETURN
      IDP2 = IDO+2
      DO 103 K=1,L1
         DO 102 I=3,IDO,2
            IC = IDP2-I
            TI5 = CC(I,3,K)+CC(IC,2,K)
            TI2 = CC(I,3,K)-CC(IC,2,K)
            TI4 = CC(I,5,K)+CC(IC,4,K)
            TI3 = CC(I,5,K)-CC(IC,4,K)
            TR5 = CC(I-1,3,K)-CC(IC-1,2,K)
            TR2 = CC(I-1,3,K)+CC(IC-1,2,K)
            TR4 = CC(I-1,5,K)-CC(IC-1,4,K)
            TR3 = CC(I-1,5,K)+CC(IC-1,4,K)
            CH(I-1,K,1) = CC(I-1,1,K)+TR2+TR3
            CH(I,K,1) = CC(I,1,K)+TI2+TI3
            CR2 = CC(I-1,1,K)+TR11*TR2+TR12*TR3
            CI2 = CC(I,1,K)+TR11*TI2+TR12*TI3
            CR3 = CC(I-1,1,K)+TR12*TR2+TR11*TR3
            CI3 = CC(I,1,K)+TR12*TI2+TR11*TI3
            CR5 = TI11*TR5+TI12*TR4
            CI5 = TI11*TI5+TI12*TI4
            CR4 = TI12*TR5-TI11*TR4
            CI4 = TI12*TI5-TI11*TI4
            DR3 = CR3-CI4
            DR4 = CR3+CI4
            DI3 = CI3+CR4
            DI4 = CI3-CR4
            DR5 = CR2+CI5
            DR2 = CR2-CI5
            DI5 = CI2-CR5
            DI2 = CI2+CR5
            CH(I-1,K,2) = WA1(I-2)*DR2-WA1(I-1)*DI2
            CH(I,K,2) = WA1(I-2)*DI2+WA1(I-1)*DR2
            CH(I-1,K,3) = WA2(I-2)*DR3-WA2(I-1)*DI3
            CH(I,K,3) = WA2(I-2)*DI3+WA2(I-1)*DR3
            CH(I-1,K,4) = WA3(I-2)*DR4-WA3(I-1)*DI4
            CH(I,K,4) = WA3(I-2)*DI4+WA3(I-1)*DR4
            CH(I-1,K,5) = WA4(I-2)*DR5-WA4(I-1)*DI5
            CH(I,K,5) = WA4(I-2)*DI5+WA4(I-1)*DR5
  102    CONTINUE
  103 CONTINUE
      RETURN
      END
      SUBROUTINE RADBG (IDO,IP,L1,IDL1,CC,C1,C2,CH,CH2,WA)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION       CH(IDO,L1,IP)          ,CC(IDO,IP,L1)          ,
     1                C1(IDO,L1,IP)          ,C2(IDL1,IP),
     2                CH2(IDL1,IP)           ,WA(1)
      DATA TPI/6.28318530717958647692D0/
      ARG = TPI/FLOAT(IP)
      DCP = COS(ARG)
      DSP = SIN(ARG)
      IDP2 = IDO+2
      NBD = (IDO-1)/2
      IPP2 = IP+2
      IPPH = (IP+1)/2
      IF (IDO .LT. L1) GO TO 103
      DO 102 K=1,L1
         DO 101 I=1,IDO
            CH(I,K,1) = CC(I,1,K)
  101    CONTINUE
  102 CONTINUE
      GO TO 106
  103 DO 105 I=1,IDO
         DO 104 K=1,L1
            CH(I,K,1) = CC(I,1,K)
  104    CONTINUE
  105 CONTINUE
  106 DO 108 J=2,IPPH
         JC = IPP2-J
         J2 = J+J
         DO 107 K=1,L1
            CH(1,K,J) = CC(IDO,J2-2,K)+CC(IDO,J2-2,K)
            CH(1,K,JC) = CC(1,J2-1,K)+CC(1,J2-1,K)
  107    CONTINUE
  108 CONTINUE
      IF (IDO .EQ. 1) GO TO 116
      IF (NBD .LT. L1) GO TO 112
      DO 111 J=2,IPPH
         JC = IPP2-J
         DO 110 K=1,L1
            DO 109 I=3,IDO,2
               IC = IDP2-I
               CH(I-1,K,J) = CC(I-1,2*J-1,K)+CC(IC-1,2*J-2,K)
               CH(I-1,K,JC) = CC(I-1,2*J-1,K)-CC(IC-1,2*J-2,K)
               CH(I,K,J) = CC(I,2*J-1,K)-CC(IC,2*J-2,K)
               CH(I,K,JC) = CC(I,2*J-1,K)+CC(IC,2*J-2,K)
  109       CONTINUE
  110    CONTINUE
  111 CONTINUE
      GO TO 116
  112 DO 115 J=2,IPPH
         JC = IPP2-J
         DO 114 I=3,IDO,2
            IC = IDP2-I
            DO 113 K=1,L1
               CH(I-1,K,J) = CC(I-1,2*J-1,K)+CC(IC-1,2*J-2,K)
               CH(I-1,K,JC) = CC(I-1,2*J-1,K)-CC(IC-1,2*J-2,K)
               CH(I,K,J) = CC(I,2*J-1,K)-CC(IC,2*J-2,K)
               CH(I,K,JC) = CC(I,2*J-1,K)+CC(IC,2*J-2,K)
  113       CONTINUE
  114    CONTINUE
  115 CONTINUE
  116 AR1 = 1.0D0
      AI1 = 0.0D0
      DO 120 L=2,IPPH
         LC = IPP2-L
         AR1H = DCP*AR1-DSP*AI1
         AI1 = DCP*AI1+DSP*AR1
         AR1 = AR1H
         DO 117 IK=1,IDL1
            C2(IK,L) = CH2(IK,1)+AR1*CH2(IK,2)
            C2(IK,LC) = AI1*CH2(IK,IP)
  117    CONTINUE
         DC2 = AR1
         DS2 = AI1
         AR2 = AR1
         AI2 = AI1
         DO 119 J=3,IPPH
            JC = IPP2-J
            AR2H = DC2*AR2-DS2*AI2
            AI2 = DC2*AI2+DS2*AR2
            AR2 = AR2H
            DO 118 IK=1,IDL1
               C2(IK,L) = C2(IK,L)+AR2*CH2(IK,J)
               C2(IK,LC) = C2(IK,LC)+AI2*CH2(IK,JC)
  118       CONTINUE
  119    CONTINUE
  120 CONTINUE
      DO 122 J=2,IPPH
         DO 121 IK=1,IDL1
            CH2(IK,1) = CH2(IK,1)+CH2(IK,J)
  121    CONTINUE
  122 CONTINUE
      DO 124 J=2,IPPH
         JC = IPP2-J
         DO 123 K=1,L1
            CH(1,K,J) = C1(1,K,J)-C1(1,K,JC)
            CH(1,K,JC) = C1(1,K,J)+C1(1,K,JC)
  123    CONTINUE
  124 CONTINUE
      IF (IDO .EQ. 1) GO TO 132
      IF (NBD .LT. L1) GO TO 128
      DO 127 J=2,IPPH
         JC = IPP2-J
         DO 126 K=1,L1
            DO 125 I=3,IDO,2
               CH(I-1,K,J) = C1(I-1,K,J)-C1(I,K,JC)
               CH(I-1,K,JC) = C1(I-1,K,J)+C1(I,K,JC)
               CH(I,K,J) = C1(I,K,J)+C1(I-1,K,JC)
               CH(I,K,JC) = C1(I,K,J)-C1(I-1,K,JC)
  125       CONTINUE
  126    CONTINUE
  127 CONTINUE
      GO TO 132
  128 DO 131 J=2,IPPH
         JC = IPP2-J
         DO 130 I=3,IDO,2
            DO 129 K=1,L1
               CH(I-1,K,J) = C1(I-1,K,J)-C1(I,K,JC)
               CH(I-1,K,JC) = C1(I-1,K,J)+C1(I,K,JC)
               CH(I,K,J) = C1(I,K,J)+C1(I-1,K,JC)
               CH(I,K,JC) = C1(I,K,J)-C1(I-1,K,JC)
  129       CONTINUE
  130    CONTINUE
  131 CONTINUE
  132 CONTINUE
      IF (IDO .EQ. 1) RETURN
      DO 133 IK=1,IDL1
         C2(IK,1) = CH2(IK,1)
  133 CONTINUE
      DO 135 J=2,IP
         DO 134 K=1,L1
            C1(1,K,J) = CH(1,K,J)
  134    CONTINUE
  135 CONTINUE
      IF (NBD .GT. L1) GO TO 139
      IS = -IDO
      DO 138 J=2,IP
         IS = IS+IDO
         IDIJ = IS
         DO 137 I=3,IDO,2
            IDIJ = IDIJ+2
            DO 136 K=1,L1
               C1(I-1,K,J) = WA(IDIJ-1)*CH(I-1,K,J)-WA(IDIJ)*CH(I,K,J)
               C1(I,K,J) = WA(IDIJ-1)*CH(I,K,J)+WA(IDIJ)*CH(I-1,K,J)
  136       CONTINUE
  137    CONTINUE
  138 CONTINUE
      GO TO 143
  139 IS = -IDO
      DO 142 J=2,IP
         IS = IS+IDO
         DO 141 K=1,L1
            IDIJ = IS
            DO 140 I=3,IDO,2
               IDIJ = IDIJ+2
               C1(I-1,K,J) = WA(IDIJ-1)*CH(I-1,K,J)-WA(IDIJ)*CH(I,K,J)
               C1(I,K,J) = WA(IDIJ-1)*CH(I,K,J)+WA(IDIJ)*CH(I-1,K,J)
  140       CONTINUE
  141    CONTINUE
  142 CONTINUE
  143 RETURN
      END
      SUBROUTINE DZFFTF (N,R,AZERO,A,B,WSAVE)
C
C                       VERSION 3  JUNE 1979
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION       R(*)       ,A(*)       ,B(*)       ,WSAVE(*)
      IF (N-2) 101,102,103
  101 AZERO = R(1)
      RETURN
  102 AZERO = 0.5D0*(R(1)+R(2))
      A(1) = 0.5D0*(R(1)-R(2))
      RETURN
  103 DO 104 I=1,N
         WSAVE(I) = R(I)
  104 CONTINUE
      CALL DFFTF (N,WSAVE,WSAVE(N+1))
      CF = 2.0D0/FLOAT(N)
      CFM = -CF
      AZERO = 0.5D0*CF*WSAVE(1)
      NS2 = (N+1)/2
      NS2M = NS2-1
      DO 105 I=1,NS2M
         A(I) = CF*WSAVE(2*I)
         B(I) = CFM*WSAVE(2*I+1)
  105 CONTINUE
      IF (MOD(N,2) .EQ. 1) RETURN
      A(NS2) = 0.5D0*CF*WSAVE(N)
      B(NS2) = 0.0D0
      RETURN
      END


text:
! Check that we correctly simplify IS_IOSTAT_END and IS_IOSTAT_EOR.
! Not very useful, but required by the standards
!
! This test relies on the error numbers for END and EOR being -1 and -2.
! This is good to actual
!
! { dg-do compile }
!

  use iso_fortran_env, only : iostat_end, iostat_eor
  implicit none

  integer(kind=merge(4, 0, is_iostat_end(-1))) :: a
  integer(kind=merge(4, 0, is_iostat_end(-1_1))) :: b
  integer(kind=merge(4, 0, is_iostat_end(-1_2))) :: c
  integer(kind=merge(4, 0, is_iostat_end(-1_4))) :: d
  integer(kind=merge(4, 0, is_iostat_end(-1_8))) :: e

  integer(kind=merge(4, 0, is_iostat_eor(-2))) :: f
  integer(kind=merge(4, 0, is_iostat_eor(-2_1))) :: g
  integer(kind=merge(4, 0, is_iostat_eor(-2_2))) :: h
  integer(kind=merge(4, 0, is_iostat_eor(-2_4))) :: i
  integer(kind=merge(4, 0, is_iostat_eor(-2_8))) :: j

  integer(kind=merge(0, 4, is_iostat_eor(-1))) :: k
  integer(kind=merge(0, 4, is_iostat_end(-2))) :: l

  integer(kind=merge(0, 4, is_iostat_eor(0))) :: m
  integer(kind=merge(0, 4, is_iostat_end(0))) :: n

  integer(kind=merge(4, 0, is_iostat_end(0))) :: o ! { dg-error "not supported for type" }
  integer(kind=merge(4, 0, is_iostat_eor(0))) :: p ! { dg-error "not supported for type" }

  integer(kind=merge(4, 0, is_iostat_eor(iostat_eor))) :: q
  integer(kind=merge(4, 0, is_iostat_end(iostat_end))) :: r
  integer(kind=merge(0, 4, is_iostat_end(iostat_eor))) :: s
  integer(kind=merge(0, 4, is_iostat_eor(iostat_end))) :: t

  end

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
      SUBROUTINE CAPIT(STR)

C$$$  SUBPROGRAM DOCUMENTATION BLOCK
C
C SUBPROGRAM:    CAPIT
C   PRGMMR: WOOLLEN          ORG: NP20       DATE: 2002-05-14
C
C ABSTRACT: THIS SUBROUTINE CAPITALIZES A STRING OF CHARACTERS. THIS
C   ENABLES THE USE OF MIXED CASE IN THE UNIT SECTION OF THE ASCII
C   BUFR TABLES.  AN EXAMPLE: A PROGRAM WHICH GENERATES AN ASCII BUFR
C   TABLE FROM THE "MASTER TABLE B" MIGHT END UP COPYING SOME UNITS
C   FIELDS IN MIXED OR LOWER CASE.  IF THE UNITS ARE 'CODE TABLE' OR
C   'FLAG TABLE' OR CERTAIN OTHER UNIT DESIGNATIONS, THE TABLE WILL BE
C   PARSED INCORRECTLY, AND THE DATA READ OR INCORRECTLY AS A RESULT.
C   THIS MAKES SURE ALL UNIT DESIGNATIONS ARE SEEN BY THE PARSER IN
C   UPPER CASE TO AVOID THESE TYPES OF PROBLEMS.
C
C PROGRAM HISTORY LOG:
C 2002-05-14  J. WOOLLEN -- ORIGINAL AUTHOR
C 2003-11-04  S. BENDER  -- ADDED REMARKS/BUFRLIB ROUTINE
C                           INTERDEPENDENCIES
C 2003-11-04  D. KEYSER  -- UNIFIED/PORTABLE FOR WRF; ADDED
C                           DOCUMENTATION (INCLUDING HISTORY)
C 2012-03-02  J. ATOR    -- CHANGED NAME OF UPS ARRAY TO UPCS TO AVOID
C                           NAMESPACE CONTENTION WITH NEW FUNCTION UPS
C
C USAGE:    CALL CAPIT (STR)
C   INPUT ARGUMENT LIST:
C     STR      - CHARACTER*(*): STRING POSSIBLY CONTAINING MIXED UPPER-
C                AND LOWER-CASE CHARACTERS
C
C   OUTPUT ARGUMENT LIST:
C     STR      - CHARACTER*(*): SAME STRING AS INPUT BUT NOW CONTAINING
C                ALL UPPER-CASE CHARACTERS
C
C REMARKS:
C    THIS ROUTINE CALLS:        None
C    THIS ROUTINE IS CALLED BY: CMPMSG   ELEMDX   STBFDX   STDMSG
C                               STRCPT
C                               Normally not called by any application
C                               programs but it could be.
C
C ATTRIBUTES:
C   LANGUAGE: FORTRAN 77
C   MACHINE:  PORTABLE TO ALL PLATFORMS
C
C$$$

      CHARACTER*(*) STR
      CHARACTER*26 UPCS,LWCS
      DATA UPCS/'ABCDEFGHIJKLMNOPQRSTUVWXYZ'/
      DATA LWCS/'abcdefghijklmnopqrstuvwxyz'/

      DO 20 I=1,LEN(STR)
      DO 10 J=1,26
      IF(STR(I:I).EQ.LWCS(J:J)) THEN
         STR(I:I) = UPCS(J:J)
         GOTO 20
      ENDIF
10    CONTINUE
20    CONTINUE

      RETURN
      END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
subroutine lev_reampi(a)
   use typre
   use MPI
   use Mod_Memor
   use Mod_Listen
   use Mod_PhysicalProblem
   use Mod_LevelSet
   implicit none

   class(LevelSetProblem) :: a
   integer :: ierr,igauge,ndime

   !Communicate lev_reaphy
   CALL MPI_BCAST(a%kfl_advec, 1, MPI_INTEGER4, a%MPIroot, a%MPIcomm, ierr)
   CALL MPI_BCAST(a%kfl_ReinitLevel, 1, MPI_INTEGER4, a%MPIroot, a%MPIcomm, ierr)
   CALL MPI_BCAST(a%kfl_ExactLevel, 1, MPI_INTEGER4, a%MPIroot, a%MPIcomm, ierr)

   !Communicate lev_reanut
   CALL MPI_BCAST(a%staco, size(a%staco), MPI_REAL8, a%MPIroot, a%MPIcomm, ierr)
   CALL MPI_BCAST(a%NstepsReinitLevel, 1, MPI_INTEGER4, a%MPIroot, a%MPIcomm, ierr)
   CALL MPI_BCAST(a%kfl_ForceEulerianAdvection, 1, MPI_INTEGER4, a%MPIroot, a%MPIcomm, ierr)

   !Adaptivity
   CALL MPI_BCAST(a%OutLayers, 1, MPI_REAL8, a%MPIroot, a%MPIcomm, ierr)
   CALL MPI_BCAST(a%GeneralRefinementLevels, 1, MPI_REAL8, a%MPIroot, a%MPIcomm, ierr)
   CALL MPI_BCAST(a%InterfaceRefinementLevels, 1, MPI_REAL8, a%MPIroot, a%MPIcomm, ierr)

   !Mass Correction
   CALL MPI_BCAST(a%kfl_MassCorrection, 1, MPI_INTEGER4, a%MPIroot, a%MPIcomm, ierr)

   !Initial Redistancing
   CALL MPI_BCAST(a%kfl_InitialRedistance, 1, MPI_INTEGER4, a%MPIroot, a%MPIcomm, ierr)

   !Communicate lev_reaous
   !HeightGauges
   CALL MPI_BCAST(a%nHeightGauges, 1, MPI_INTEGER4, a%MPIroot, a%MPIcomm, ierr)
   if (a%MPIrank /= a%MPIroot) then
      allocate(a%HeightGauges(a%nHeightGauges))
      call a%Memor%allocObj(0,'HeightGauges','lev_reaous',a%nHeightGauges)
   endif
   call a%Mesh%GetNdime(ndime)
   do igauge = 1,a%nHeightGauges
      CALL MPI_BCAST(a%HeightGauges(igauge)%Origin, size(a%HeightGauges(igauge)%Origin), MPI_REAL8, a%MPIroot, a%MPIcomm, ierr)
      CALL MPI_BCAST(a%HeightGauges(igauge)%DirectionVector, size(a%HeightGauges(igauge)%DirectionVector), MPI_REAL8, a%MPIroot, a%MPIcomm, ierr)
      call a%HeightGauges(igauge)%SetMPI(a%MPIcomm,a%MPIsize,a%MPIroot,a%Mpirank)
      call a%HeightGauges(igauge)%SetNdime(ndime)
   enddo

   !For fixing CN
   if (a%kfl_tsche_1st_datafile == 'CNOBS') a%kfl_tsche_1st_datafile = 'CN   '

end subroutine

max_stars_repo_name: ciaid-colombia/InsFEM

max_issues_repo_name: ciaid-colombia/InsFEM

max_forks_repo_name: ciaid-colombia/InsFEM


text:
      subroutine tprca(abeta,acphi,alamca,alpha,alphca,amua2c,amuc2a,
     $ icapr,ipbtmx,jpfcmx,natmax,na2ctr,ncapr,ncvca,nc2atr,nerr,
     $ noutpt,npx2mx,npx2t,nttyo,nwarn,pcvca,qpdca,uaqsp,upair,zaqsp)
c
c     Test and process the Pitzer data for ca (cation, anion) pairs
c     read from the DATA0 file. Find and flag errors, such as duplication
c     of data (e.g., two data blocks for the same ca pair). Calculate
c     the conventional primitive Pitzer parameters from the observable
c     compound parameters read from the data file:
c
c       beta(n)(ca) -> lambda(n)(ca)   (n = 0,2)
c       Cphi(ca)    -> mu(cca) and mu(aac)
c
c     Check the coverage of entered Pitzer data against all possible
c     ca pairs that can be composed of the aqueous cations and anions
c     present on the data file.
c
c     This subroutine is called by:
c
c       EQPT/eqpt.f
c
c-----------------------------------------------------------------------
c
c     Principal input:
c
c       nat    = the number of aqueous species
c       uaqsp  = array of names of aqueous species
c
c     Principal output:
c
c       nerr   = cumulative error counter
c       nwarn  = cumulative warning counter
c
c-----------------------------------------------------------------------
c
      implicit none
c
c-----------------------------------------------------------------------
c
c     Calling sequence variable declarations.
c
      integer ipbtmx,jpfcmx,natmax,na2ctr,ncapr,nc2atr,npx2mx
c
      integer noutpt,nttyo
c
      integer icapr(2,ncapr)
c
      integer ncvca,nerr,npx2t,nwarn
c
      logical qpdca(ncapr)
c
      character(len=24) uaqsp(natmax),upair(2,npx2mx)
c
      real*8 abeta(jpfcmx,0:ipbtmx,npx2mx),acphi(jpfcmx,npx2mx),
     $ alamca(jpfcmx,0:ipbtmx,ncapr),alphca(ipbtmx,ncapr),
     $ alpha(2,npx2mx),amua2c(jpfcmx,na2ctr),amuc2a(jpfcmx,nc2atr),
     $ zaqsp(natmax)
c
      real*8 pcvca
c
c-----------------------------------------------------------------------
c
c     Local variable declarations.
c
      integer i,j,jp,jpair,j2,j3,j4,j5,n,ncount,ndupl,nlistl,nn,
     $ nodatc
c
      integer ilnobl
c
      character*56 ustr56
      character*24 unam1,unam2
      character*8 ux8
c
      real*8 fmuc2a,fmua2c,z1,z2
c
c-----------------------------------------------------------------------
c
c     Limit on the list of pairs for which no data were found.
c
      data nlistl / 20 /
c
c-----------------------------------------------------------------------
c
c     Initialize the data arrays.
c
      do n = 1,ncapr
        qpdca(n) = .false.
        do i = 1,ipbtmx
          alphca(i,n) = 0.
        enddo
      enddo
c
      do n = 1,ncapr
        do i = 0,ipbtmx
          do j = 1,jpfcmx
            alamca(j,i,n) = 0.
          enddo
        enddo
      enddo
c
      do n = 1,ncapr
        do j = 1,jpfcmx
          amuc2a(j,n) = 0.
          amua2c(j,n) = 0.
        enddo
      enddo
c
c     Check the entered data for ca pairs.
c
      nodatc = 0
c
      do n = 1,ncapr
        i = icapr(1,n)
        j = icapr(2,n)
        unam1 = uaqsp(i)
        unam2 = uaqsp(j)
        z1 = zaqsp(i)
        z2 = zaqsp(j)
c
c       Search for unam1, unam2 in the upair array.
c       That array corresponds to the species pairs blocks.
c
        call srch22(jpair,unam1,unam2,upair,npx2mx,npx2t)
c
        if (jpair .gt. 0) then
c
c         Have found an entry.
c
          qpdca(n) = .true.
c
c         Store the data.
c
          do i = 1,ipbtmx
            alphca(i,n) = alpha(i,jpair)
          enddo
c
          do i = 0,ipbtmx
            do j = 1,jpfcmx
              alamca(j,i,n) = abeta(j,i,jpair)
            enddo
          enddo
c
          fmuc2a = sqrt(-z1/z2)/6.
          fmua2c = sqrt(-z2/z1)/6.
c
          do j = 1,jpfcmx
            amuc2a(j,n)   = acphi(j,jpair)*fmuc2a
            amua2c(j,n)   = acphi(j,jpair)*fmua2c
          enddo
c
c         Check for duplicate data sets.
c
          ndupl = 0
          do jp = jpair + 1,npx2t
            if (unam1(1:24) .eq. upair(1,jp)(1:24)) then
              if (unam2(1:24) .eq. upair(2,jp)(1:24)) then
                ndupl = ndupl + 1
              endif
            endif
          enddo
c
          if (ndupl .gt. 0) then
            j2 = ilnobl(unam1)
            j3 = ilnobl(unam2)
            if (ndupl .eq. 1) then
              write (noutpt,1010) unam1(1:j2),unam2(1:j3)
              write (nttyo,1010) unam1(1:j2),unam2(1:j3)
 1010         format(/' * Error - (EQPT/tprca) Have found a duplicate',
     $        ' data block on the DATA0 file',/7x,'for the ca pair ',
     $        a,', ',a,'.')
            else
              ux8 = ' '
              write (ux8,'(i5)') ndupl
              call lejust(ux8)
              j5 = ilnobl(ux8)
              write (noutpt,1020) ux8(1:j5),unam1(1:j2),unam2(1:j3)
              write (nttyo,1020) ux8(1:j5),unam1(1:j2),unam2(1:j3)
 1020         format(/' * Error - (EQPT/tprca) Have found ',a,
     $        ' duplicate data blocks on the DATA0 file',
     $        /7x,'for the ca pair ',a,', ',a,'.')
            endif
            nerr = nerr + ndupl
          endif
c
        else
c
c         No data block was found on the DATA0 file.
c         Note that qpdca(n) is left with a value of .false.
c
          nodatc = nodatc + 1
        endif
c
      enddo
c
c* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
      if (nodatc .gt. 0) then
        write (noutpt,1040)
        write (nttyo,1040)
 1040   format(/' * Warning - (EQPT/tprca) Did not find a data',
     $  ' block on the DATA0 file',/7x,'for any of the following',
     $  ' ca pairs:',/)
c
        ncount = 0
        do n = 1,ncapr
          if (.not.qpdca(n)) then
            ncount = ncount + 1
            i = icapr(1,n)
            j = icapr(2,n)
            unam1 = uaqsp(i)
            unam2 = uaqsp(j)
            j2 = ilnobl(unam1)
            j3 = ilnobl(unam2)
            ustr56 = unam1(1:j2) // ', ' // unam2(1:j3)
            j4 = ilnobl(ustr56)
            write (noutpt,1050) ustr56(1:j4)
            write (nttyo,1050) ustr56(1:j4)
 1050       format(9x,a)
            if (ncount .eq. nlistl) go to 200
          endif
        enddo
  200   continue
c
        nn = nodatc - ncount
        if (nn .gt. 0) then
          write (ux8,'(i5)') nn
          call lejust(ux8)
          j3 = ilnobl(ux8)
          write (noutpt,1060) ux8(1:j3)
          write (nttyo,1060) ux8(1:j3)
 1060     format(/9x,'plus ',a,' others')
        endif
        write (noutpt,1070)
        write (nttyo,1070)
 1070   format(1x)
        nwarn = nwarn + 1
      endif
c
      ncvca = ncapr - nodatc
      pcvca = (100.*ncvca)/float(ncapr)
c
c* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
      end


text:
subroutine one_phase_finite(s, E, mu, rho1, H1, km, fm, r, fn, a,layer,loc,st, m,n)
implicit none
complex*16, intent(in) :: s, fn
integer, intent(in) ::  m,N, layer, loc, st
complex*16,intent(out) :: a
complex*16 :: disp(4,4), T(4,4), GSK(n*2+2,n*2+2), s4inv(n*2+2,n*2+2), disp_st(4,4), T_st(4,4)
real*8, parameter :: pi = 4.D0*DATAN(1.D0)
integer ::  mode,NL, IS1, j, i,pa, layer_st
real*8, intent(in) ::  E(N), mu(N), rho1(N),H1(N), km(m),fm(m),r
real*8 ::  uniload(N*2+2), rho, H, unit1(4,4)
complex*16 :: CI, D3inv(2,2),uhat(N*2+2),kp, ks, stress(4),lamda,fnn1,fnn2, fnn3, stress1(4,4)
complex*16 :: D3(2,2),T3(2,2), SK3(2,2), lmu,cp2,cs2, S2INV(4,4),sk(4,4), a1(N*2+2)
real*8 :: td, ts

!f2py intent(in) s, fn, E, mu, rho1, H1, km,fm,r, layer,loc, st, m,n
!f2py intent(out) a
!f2py depend(m) km, fm
!f2py depend(n) E , mu , rho1 ,H1,GSK , s4inv, uniload, uhat

CI = (0.0 , 1.0); uniload = 0;  uniload(loc) = 1.0; a = 0; pa= n*2+2
unit1 = 0; unit1(1,1) = -1.0; unit1(2,2) = -1.0; unit1(3,3) = 1.0; unit1(4,4) = 1.0

do  mode = 1,m
    GSK = 0.0
    DO NL = 1,n     !for different layers
	H   = H1(NL)
    rho = rho1(NL)
    lamda = E(NL)*mu(NL) / ((1-2*mu(NL))*(1+mu(NL)))  !*(1+CI*0.1) artifical damping
	lmu  = (E(NL)/(2*(1+mu(NL)))) !/s*(1+0.001*s)
	CP2= (lamda+2.0*lmu )/rho
	CS2 =(lmu )/rho
    kp = sqrt( s**2 /CP2 + km(mode)**2)
	ks = sqrt( s**2 /CS2 + km(mode)**2)

    DISP(1,1) = -km(mode)
    DISP(1,2) = -km(mode)*EXP(-kp*H)
    DISP(1,3) =  ks
    DISP(1,4) = -ks*EXP(-ks*H)

    DISP(2,1) = -kp
    DISP(2,2) =  kp*EXP(-kp*H)
    DISP(2,3) =  km(mode)
    DISP(2,4) =  km(mode) *EXP(-ks*H)

    DISP(3,1) = -km(mode)*EXP(-kp*H)
    DISP(3,2) = -km(mode)
    DISP(3,3) =  ks*EXP(-ks*H)
    DISP(3,4) = - ks

    DISP(4,1) = -kp*EXP(-kp*H)
    DISP(4,2) =  kp
    DISP(4,3) = km(mode)*EXP(-ks*H)
    DISP(4,4) = km(mode)

!	*******************STRESS MATRIX**********************************
    T(1,1) =  2.0*kp*km(mode)*lmu
    T(1,2) =  -2.0*kp*km(mode)*lmu*EXP(-kp*H)
    T(1,3) = -(ks**2+km(mode)**2)*lmu
    T(1,4) = -(ks**2+km(mode)**2)*EXP(-ks*H)*lmu

    T(2,1) =  (kp**2*(lamda+2*lmu)-km(mode)**2*lamda)
    T(2,2) =  (kp**2*(lamda+2*lmu)-km(mode)**2*lamda)*EXP(-kp*H)
    T(2,3) = -2*km(mode)*ks*lmu
    T(2,4) = 2*km(mode)*ks*EXP(-ks*H)*lmu

    T(3,1) = 2.0*kp*EXP(-kp*H)*km(mode)*lmu
    T(3,2) = -2.0*kp*km(mode)*lmu
    T(3,3) = -(ks**2+km(mode)**2)*EXP(-ks*H)*lmu
    T(3,4) = -(ks**2+km(mode)**2)*lmu

    T(4,1) =  ((kp**2*(lamda+2*lmu)-km(mode)**2 *lamda))*EXP(-kp*H)
    T(4,2) =  kp**2*(lamda+2*lmu)-km(mode)**2*lamda
    T(4,3) = -2*km(mode)*ks*EXP(-ks*H)*lmu
    T(4,4) =  2*km(mode) *ks *lmu

!	******** FIND INVERSE OF DISPLACEMENT MATRIX **********************
    call zInverse(DISP, S2INV, 4)
    SK = matmul(T,S2INV)   !MULTIPLICATION TO GET STIFFNESS MATRIX OF AN ELEMENT
    if ((Int(layer/4)+1) == NL) then
        T_st = SK
    end if
    SK = MATMUL(UNIT1, SK)
    IS1=2*(NL-1)  !NL IS number

    DO I = 1,4
        DO J = 1,4
            GSK(I+IS1,J+IS1) = GSK(I+IS1,J+IS1) + SK(I,J)
        ENDDO
    ENDDO
    END DO		 !END NL LOOP

    call zInverse(GSK, s4inv,n*2+2)
    uhat =  matmul(s4inv,uniload)
    ! displacement
    if (st==0 .and. mod(layer,2)>0) then
    a  =  a  + uhat(layer) * Fm(mode)* Fn * DBESJ1(km(mode)*r)
    else if (st==0 .and. mod(layer,2)==0) then
    a  =  a  + uhat(layer) * Fm(mode)* Fn * DBESJ0(km(mode)*r)
    end if

    stress =  matmul(T_st, uhat((Int(layer/4)+1)*2+2-3:(Int(layer/4)+1)*2+2))

    if ((st==1) .and. (mod(layer+4,2)==1)) then
    a =  a + stress(mod(layer+4,4))* Fm(mode)* Fn * DBESJ1(km(mode)*r)
    else if ((st==1) .and. (mod(layer+4,2) .NE. 1)) then
    a =  a + stress((mod(layer+4,4)))* Fm(mode)* Fn * DBESJ0(km(mode)*r)
    end if
end do
end subroutine

subroutine one_phase_infinite(s, E, mu, rho1, H1, km, fm, r, fn, a,layer,loc,st,m,n)
implicit none
complex*16, intent(in) :: s, fn
integer, intent(in) ::  m,N, layer, loc,st
complex*16,intent(out) :: a
complex*16 :: disp(4,4), T(4,4), GSK(n*2,n*2), s4inv(n*2,n*2)
real*8, parameter :: pi = 4.D0*DATAN(1.D0)
integer ::  mode,NL, IS1, j, i,pa, layer_st
real*8, intent(in) ::  E(N), mu(N), rho1(N),H1(N), km(m),fm(m),r
real*8 ::  uniload(N*2), rho, H, unit1(4,4)
complex*16 :: CI, D3inv(2,2),uhat(N*2),kp, ks, stress(4), stress2(2), lamda,fnn1,fnn2, fnn3, T_st(4,4)
complex*16 :: D3(2,2),T3(2,2), SK3(2,2), lmu,cp2,cs2, S2INV(4,4),sk(4,4), a1(N*2), T_st2(2,2)
real*8 :: td, ts

!f2py intent(in) s, fn, E, mu, rho1, H1, km,fm,r, layer, loc,st, m,n
!f2py intent(out) a
!f2py depend(m) km, fm
!f2py depend(n) E , mu , rho1 ,H1,GSK , s4inv, uniload, uhat

CI = (0.0 , 1.0); uniload = 0;  uniload(loc) = 1.0; a = 0; pa= n*2
unit1 = 0; unit1(1,1) = -1.0; unit1(2,2) = -1.0; unit1(3,3) = 1.0; unit1(4,4) = 1.0

do mode = 1,m
    GSK = 0.0
    DO NL = 1,n                           !for different layers
		H   = H1(NL)
        rho = rho1(NL)
        lamda = E(NL)*mu(NL) / ((1-2*mu(NL))*(1+mu(NL)))  !*(1+CI*0.1) artifical damping
		lmu  = (E(NL)/(2*(1+mu(NL)))) !/s*(1+0.001*s)
		CP2= (lamda+2.0*lmu )/rho
		CS2 =(lmu )/rho
        kp = sqrt( s**2 /CP2 + km(mode)**2)
		ks = sqrt( s**2 /CS2 + km(mode)**2)

!*************************For Layers With Finite Thickness******************************
        IF (NL<=n-1) THEN
            DISP(1,1) = -km(mode)
            DISP(1,2) = -km(mode)*EXP(-kp*H)
            DISP(1,3) =  ks
            DISP(1,4) = -ks*EXP(-ks*H)

            DISP(2,1) = -kp
            DISP(2,2) =  kp*EXP(-kp*H)
            DISP(2,3) =  km(mode)
            DISP(2,4) =  km(mode) *EXP(-ks*H)

            DISP(3,1) = -km(mode)*EXP(-kp*H)
            DISP(3,2) = -km(mode)
            DISP(3,3) =  ks*EXP(-ks*H)
            DISP(3,4) = - ks

            DISP(4,1) = -kp*EXP(-kp*H)
            DISP(4,2) =  kp
            DISP(4,3) = km(mode)*EXP(-ks*H)
            DISP(4,4) = km(mode)

        !	*******************STRESS MATRIX**********************************
            T(1,1) =  2.0*kp*km(mode)*lmu
            T(1,2) =  -2.0*kp*km(mode)*lmu*EXP(-kp*H)
            T(1,3) = -(ks**2+km(mode)**2)*lmu
            T(1,4) = -(ks**2+km(mode)**2)*EXP(-ks*H)*lmu

            T(2,1) =  (kp**2*(lamda+2*lmu)-km(mode)**2*lamda)
            T(2,2) =  (kp**2*(lamda+2*lmu)-km(mode)**2*lamda)*EXP(-kp*H)
            T(2,3) = -2*km(mode)*ks*lmu
            T(2,4) = 2*km(mode)*ks*EXP(-ks*H)*lmu

            T(3,1) = 2.0*kp*EXP(-kp*H)*km(mode)*lmu
            T(3,2) = -2.0*kp*km(mode)*lmu
            T(3,3) = -(ks**2+km(mode)**2)*EXP(-ks*H)*lmu
            T(3,4) = -(ks**2+km(mode)**2)*lmu

            T(4,1) =  ((kp**2*(lamda+2*lmu)-km(mode)**2 *lamda))*EXP(-kp*H)
            T(4,2) =  kp**2*(lamda+2*lmu)-km(mode)**2*lamda
            T(4,3) = -2*km(mode)*ks*EXP(-ks*H)*lmu
            T(4,4) =  2*km(mode) *ks *lmu

!******** FIND INVERSE OF DISPLACEMENT MATRIX **********************

            call zInverse(DISP, S2INV, 4)
            SK = matmul(T,S2INV)   !	MULTIPLICATION TO GET STIFFNESS MATRIX OF AN ELEMENT
            if ((Int(layer/4)+2) == NL) then
                T_st = SK
            end if
            SK = MATMUL(UNIT1, SK)

            IS1=2*(NL-1)  ! NL IS number

            DO I = 1,4
                DO J = 1,4
                    GSK(I+IS1,J+IS1) = GSK(I+IS1,J+IS1) + SK(I,J)
                ENDDO
            ENDDO
!********************FOR LAST LAYER - 2DEGREE OF FREEDOM****************
        ELSE
            D3(1,1) = -km(mode)
            D3(1,2) = ks
            D3(2,1) = -kp
            D3(2,2) = km(mode)

            T3(1,1) =  2.0*lmu*km(mode)*kp
            T3(1,2) = -lmu *(ks**2+km(mode)**2)
            T3(2,1) = (kp**2*(lamda+2*lmu)-km(mode)**2*lamda)
            T3(2,2) = -2.0*lmu *km(mode)*ks

            call zInverse(D3, d3inv, 2)
            SK3 = matmul(T3,D3inv)
            if (n==1) then
            T_st2 = SK3
            end if
            SK3 = -SK3
            GSK(pa-1,pa-1) = GSK(pa-1,pa-1)+ SK3(1,1)
            GSK(pa-1,pa) = GSK(pa-1,pa)+ SK3(1,2)
            GSK(pa,pa-1) = GSK(pa,pa-1)+ SK3(2,1)
            GSK(pa,pa) = GSK(pa,pa)+ SK3(2,2)
        END IF
    END DO		 !END NL LOOP

    call zInverse(GSK, s4inv,n*2)
    uhat =  matmul(s4inv,uniload)

    ! displacement
    if (st==0 .and. mod(layer,2)>0) then
    a  =  a  + uhat(layer) * Fm(mode)* Fn * DBESJ1(km(mode)*r)
    else if (st==0 .and. mod(layer,2)==0) then
    a  =  a  + uhat(layer) * Fm(mode)* Fn * DBESJ0(km(mode)*r)
    end if

    if (n==1) then
    stress2 =  matmul(T_st2, uhat(1:2))
    else
    stress =  matmul(T_st, uhat((Int(layer/4)+2)*2-3:(Int(layer/4)+2)*2))
    end if

    if (n==1 .and. (st==1) .and. (mod(layer+4,2)==1)) then
    a =  a + stress2(layer)* Fm(mode)* Fn * DBESJ1(km(mode)*r)
    else if (n==1 .and. (st==1) .and. (mod(layer+4,2) .NE. 1)) then
    a =  a + stress2(layer)* Fm(mode)* Fn * DBESJ0(km(mode)*r)
    end if

    if (n>1 .and. (st==1) .and. (mod(layer+4,2)==1)) then
    a =  a + stress(mod(layer+4,4))* Fm(mode)* Fn * DBESJ1(km(mode)*r)
    else if (n>1 .and. (st==1) .and. (mod(layer+4,2) .NE. 1)) then
    a =  a + stress((mod(layer+4,4)))* Fm(mode)* Fn * DBESJ0(km(mode)*r)
    end if
end do
end subroutine

subroutine one_phase_dispersion(s, E, mu, rho1, H1, km, GSK,n)
implicit none
complex*16, intent(in) :: s
integer, intent(in) :: N
complex*16,intent(out) :: GSK(n*2,n*2)
complex*16 :: fn, disp(4,4), T(4,4)
real*8, parameter :: pi = 4.D0*DATAN(1.D0)
integer ::  mode, NL, IS1, j, i,pa
real*8, intent(in) ::  E(N), mu(N), rho1(N),H1(N), km
real*8 :: uniload(N*2), rho, H, unit1(4,4), d0, d1, mm
complex*16 :: CI, T3inv(2,2), D3inv(2,2),uhat(5),kp, ks, S3INV(4,4), stress(2),lamda
complex*16 :: D3(2,2),T3(2,2), SK3(2,2), lmu,cp2,cs2, S2INV(4,4),sk(4,4)

!f2py intent(in) s, E, mu, rho1, H1, km,fm, r,load,m,n
!f2py intent(out) GSK
!f2py depend(n) E , mu , rho1 ,H1,GSK, s4inv, uniload

CI = (0.0 , 1.0)
pa= n*2
unit1 = 0
unit1(1,1) = -1.0
unit1(2,2) = -1.0
unit1(3,3) = 1.0
unit1(4,4) = 1.0
GSK = 0

DO NL = 1,n                           !for different layers

H  = H1(NL)
rho = rho1(NL)

lamda = E(NL)*mu(NL)/ ((1-2*mu(NL))*(1+mu(NL)))  !*(1+damp*s*CI)
lmu  = (E(NL)/(2*(1+mu(NL))))                    !*(1+damp*s*CI)

CP2= (lamda+2.0*lmu)/rho
CS2 =(lmu)/rho

kp = sqrt( (s*ci)**2/CP2 + km**2)    !s represents the frequency in this case, not the laplace  variable
ks = sqrt( (s*ci)**2/CS2 + km**2)

if (kp .NE. 0 .and. ks .NE. 0 ) then

!*************************FOR 1-3 LAYER******************************
      IF (NL<=n-1) THEN

            DISP(1,1) = -km
            DISP(1,2) = -km*EXP(-kp*H)
            DISP(1,3) =  ks
            DISP(1,4) = -ks*EXP(-ks*H)

            DISP(2,1) = -kp
            DISP(2,2) =  kp*EXP(-kp*H)
            DISP(2,3) =  km
            DISP(2,4) =  km*EXP(-ks*H)

            DISP(3,1) = -km*EXP(-kp*H)
            DISP(3,2) = -km
            DISP(3,3) =  ks*EXP(-ks*H)
            DISP(3,4) = - ks

            DISP(4,1) = -kp*EXP(-kp*H)
            DISP(4,2) =  kp
            DISP(4,3) = km*EXP(-ks*H)
            DISP(4,4) = km

!	*******************STRESS MATRIX**********************************
            T(1,1) =  2.0*kp*km*lmu
            T(1,2) =  -2.0*kp*km*lmu*EXP(-kp*H)
            T(1,3) = -(ks**2+km**2)*lmu
            T(1,4) = -(ks**2+km**2)*EXP(-ks*H)*lmu

            T(2,1) =  (kp**2*(lamda+2*lmu)-km**2*lamda)
            T(2,2) =  (kp**2*(lamda+2*lmu)-km**2*lamda)*EXP(-kp*H)
            T(2,3) = -2*km*ks*lmu
            T(2,4) = 2*km*ks*EXP(-ks*H)*lmu

            T(3,1) = 2.0*kp*EXP(-kp*H)*km*lmu
            T(3,2) = -2.0*kp*km*lmu
            T(3,3) = -(ks**2+km**2)*EXP(-ks*H)*lmu
            T(3,4) = -(ks**2+km**2)*lmu

            T(4,1) =  ((kp**2*(lamda+2*lmu)-km**2 *lamda))*EXP(-kp*H)
            T(4,2) =  kp**2*(lamda+2*lmu)-km**2*lamda
            T(4,3) = -2*km*ks*EXP(-ks*H)*lmu
            T(4,4) =  2*km*ks *lmu
!	********FIND INVERSE OF DISPLACEMENT MATRIX**********************

            call zInverse(DISP, S2INV, 4)

            SK = matmul(T,S2INV)   !	MULTIPLICATION TO GET STIFFNESS MATRIX OF AN ELEMENT

            SK = MATMUL(UNIT1, SK)

            IS1=2*(NL-1)  ! NL IS number

            DO I = 1,4
            DO J = 1,4
                  GSK(I+IS1,J+IS1) = GSK(I+IS1,J+IS1) + SK(I,J)
            ENDDO
            ENDDO

!!********************FOR LAST LAYER - 2DEGREE OF FREEDOM****************
    ELSE
            D3(1,1) = -km
            D3(1,2) = ks
            D3(2,1) = -kp
            D3(2,2) = km

            T3(1,1) =  2.0*lmu*km*kp
            T3(1,2) = -lmu *(ks**2+km**2)
            T3(2,1) = (kp**2*(lamda+2*lmu)-km**2*lamda)
            T3(2,2) = -2.0*lmu *km*ks

            call zInverse(D3, d3inv, 2)
            SK3 = matmul(T3,D3inv)
            SK3 = -SK3
            GSK(pa-1,pa-1) = GSK(pa-1,pa-1)+ SK3(1,1)
            GSK(pa-1,pa) = GSK(pa-1,pa)+ SK3(1,2)
            GSK(pa,pa-1) = GSK(pa,pa-1)+ SK3(2,1)
            GSK(pa,pa) = GSK(pa,pa)+ SK3(2,2)
    END IF   !FOR LAYERS

    else
            GSK = 0

    END IF
    END DO	 !END NL LOOP
end subroutine

subroutine zInverse(a,ra,n)
      integer, intent(in)::n
      integer:: lda,ipiv(n),info,lwork, i
      complex*16, intent(in)::a(n,n)
	  complex*16, intent(out)::ra(n,n)
	  complex*16 :: work(n)
	  !f2py intent(in) n, a
      !f2py intent(out) ra
      !f2py depend(n) a, ra	, ipiv
      ra=a
      lwork=n
      lda=n
      call zgetrf(n, n, ra, lda, ipiv, info)
      if(info/=0) Then
    !   open(unit = 10, file = 'matrix.txt')
      write(0,*)  'This inversion has erros'
    !   do i = 1, n
    !     write(10,*) real(a(i,:))
    !   end do
    !   close(10)
      end if

      call zgetri(n, ra, lda, ipiv, work, lwork, info)
      if(info/=0) then
    !   open(unit = 20, file = 'matrix2.txt')
      write(0,*)  'This inversion has erros'
    !   do i = 1, n
    !     write(20,*) real(a(i,:))
    !   end do
    !   close(20)
      end if
end subroutine

subroutine eigen(a,w,vr,n)
    Implicit None
    integer, intent(in):: n
    integer :: lda,ldvl,ldvr,lwork(2*N),info
    complex*16, intent(in)::a(n,n)
    complex*16, intent(out):: W(N),vr(N,N)
    complex*16 ::work(2*N), vl(N,N)
    double  precision :: rwork(2*N)

    !f2py intent(in) a, n
    !f2py intent(out) w, vr
    !f2py depend(n) a,W, vr, work, vl, rwork
    lda=N
    ldvl=1
    ldvr=N
    lwork=2*N
    rwork=2*N
    call zgeev('N', 'v', N, A, lda, w, vl, ldvl, vr, ldvr, work, lwork,rwork, info)
end subroutine

max_stars_repo_name: Siglab-code/GeoNDT

max_issues_repo_name: Siglab-code/GeoNDT

max_forks_repo_name: Siglab-code/GeoNDT


text:
	subroutine put_gemgrid(iflno,parm_in,gribid,ivc,lev1,
     &     lev2,iyear,imonth,iday,ihour,ny,nx,data,nbits,
     &	   rlat1,rlon1,rlat2,rlon2)

	INCLUDE 'GEMPRM.PRM'
	INCLUDE 'nacmn.cmn'

        real	data(*)
	character*(*) parm_in
	character parm*12,lname*12,lname1*12
	integer	  time(5),level(2),ighdr(LLGDHD),gribid
	character dattim(2)*20,reftim*11
	real    fgrid(LLMXGD)

        parm = mparms(gribid)
        iscale = mpscal(gribid)
        rmissval = rmssvl(gribid)
        lname = mvcord(ivc)
        ivscale = mvscal(ivc)
        if(iscale .ne. 0) then
           factor = 10.**float(iscale)
        else
           factor = 1.0
        endif
        if(ivscale .ne. 0) then
           vfactor = 10.**float(ivscale)
        else
           vfactor = 1.0
        endif
        if(parm .eq. ' ') then
	   if(len(parm_in).gt.12) then
	      parm = parm_in(1:12)
           else
	      parm = parm_in(1:len(parm_in))
           endif
        endif

	CALL LV_CORD ( lname, lname1, ivcord, iret)

	level(1) = lev1 * vfactor
        if (lev2 .ne. -1) then
	   level(2) = lev2 * vfactor
        else
           level(2) = lev2
        endif

	time(1) = iyear
	time(2) = imonth
	time(3) = iday
	time(4) = ihour
	time(5) = 0

	reftim = ' '
	call ti_itoc(time, reftim, iret)
        if(iret.ne.0) then
           write(*,*) 'error with time ',time
        endif
	dattim(1) = reftim
	dattim(2) = ' '

	CALL GD_RDAT (iflno, dattim, level, ivcord, parm, fgrid,
     &        igx, igy, ighdr, iret)

        DO iy = 1, ny
           ipxy = (iy-1)*nx
           ip = (ny-iy)*nx
           DO ix = 1, nx
              ip = ip + 1
              fgrid(ipxy+ix) = data(ip) * factor
c             write(*,*) 'look ',ip,' fgrid ',ipxy+ix,' data ',data(ip),
c    &           fgrid(ipxy+ix)
           END DO
        END DO

        write(*,*) dattim(1),parm,lname,lev1,lev2
	CALL GD_WPGD (iflno, fgrid, nx, ny, ighdr, dattim, level,
     &     ivcord, parm, .TRUE., MDGGRB, nbits, iret)

	return
	end


text:
C
C
*+ mapcat_verify

       subroutine mapcat_verify(status)
C      --------------------------------
C
C Verify the map catalogue
C
C Returned:
C    Status
       integer         status
C
*-

       include 'mapcat_cat.inc'
       include '/mrao/include/chrlib_functions.inc'
       include '/mrao/include/iolib_constants.inc'

C counter
       integer       n
C local status variable
       integer       istat
C local string variable
       character     string*(iolen_file)

C check status on entry
       if (status.ne.0) return

C run through the list and test each entry in turn
       call mapcat_open(status)
       do n = 1,max_cat_entries
         call mapcat_read(n,status)
         if (status.ne.0) goto 999
         if (current_map_status(ip_data).eq.true) then
           istat = 0
           call io_namfil(
     *             current_filename(1:chr_lenb(current_filename)),
     *             string,0,istat)
           if (istat.ne.0) then
             call mapcat_remove(n,.false.,status)
           end if
         end if
       end do

999    call mapcat_err(status,'mapcat_verify',' ')

       end

max_stars_repo_name: CavendishAstrophysics/anmap

max_issues_repo_name: CavendishAstrophysics/anmap

max_forks_repo_name: CavendishAstrophysics/anmap


text:
-makelib xcelium_lib/xil_defaultlib \
  "../../../../eco_flow.gen/sources_1/ip/clk_wiz_0/clk_wiz_0_clk_wiz.v" \
  "../../../../eco_flow.gen/sources_1/ip/clk_wiz_0/clk_wiz_0.v" \
-endlib
-makelib xcelium_lib/xil_defaultlib \
  glbl.v
-endlib


! DART software - Copyright UCAR. This open source software is provided
! by UCAR, "as is", without charge, subject to all terms of use at
! http://www.image.ucar.edu/DAReS/DART/DART_download

program convert_amsu_L1

! Program to convert the AMSU/A 'AIRABRAD' brightness temperatures
! from netCDF to a DART observation sequence file.
!
! See the REAMDE.rst in this directory for more information, but the
! data citation information for the dataset is
!
! Title: AIRS/Aqua L1B AMSU (A1/A2) geolocated and calibrated brightness temperatures V005
! Version: 005
! Creator: AIRS project
! Publisher: Goddard Earth Sciences Data and Information Services Center (GES DISC)
! Release Date: 2007-07-26T00:00:00.000Z
! Linkage: https://disc.gsfc.nasa.gov/datacollection/AIRABRAD_005.html
!
! The data are originally distributed in HDF-EOS format and can be converted to netCDF
! by the h4tonccf_nc4 program available from http://hdfeos.org/software/h4cflib.php

use         types_mod, only : r8, deg2rad, PI

use  obs_sequence_mod, only : obs_sequence_type, write_obs_seq, &
                              static_init_obs_sequence, destroy_obs_sequence, &
                              print_obs_seq_summary, get_num_obs

use    utilities_mod, only : initialize_utilities, &
                             error_handler, finalize_utilities, E_ERR, E_MSG, &
                             find_namelist_in_file, check_namelist_read, &
                             do_nml_file, do_nml_term, set_filename_list, &
                             nmlfileunit, get_next_filename

use     amsua_bt_mod, only : amsua_bt_granule

use amsua_netCDF_support_mod, only : initialize_amsua_netcdf, &
                                     channel_list_to_indices, &
                                     read_amsua_bt_netCDF_granule, &
                                     make_obs_sequence, &
                                     combine_sequences, &
                                     max_possible_obs, &
                                     append_or_create

implicit none

integer, parameter :: MAXFILES = 512
integer, parameter :: AMSUA_BT_CHANNEL = 15

! ----------------------------------------------------------------------
! Declare namelist parameters
! ----------------------------------------------------------------------

character(len=256) :: l1_files(MAXFILES) = ''
character(len=256) :: l1_file_list       = ''
character(len=256) :: outputfile         = 'obs_seq.amsua'

real(r8) :: lon1 =   0.0_r8,  &   !  lower longitude bound
            lon2 = 360.0_r8,  &   !  upper longitude bound
            lat1 = -90.0_r8,  &   !  lower latitude bound
            lat2 =  90.0_r8       !  upper latitude bound

character(len=8) :: channel_list(AMSUA_BT_CHANNEL) = 'null'
integer  :: cross_track_thin = 0
integer  :: along_track_thin = 0
logical  :: append_output = .true.
integer  :: verbose = 0  ! 0 is quiet, 3 is lots of output

namelist /convert_amsu_L1_nml/ l1_files, l1_file_list, &
                               outputfile, &
                               lon1, lon2, lat1, lat2, &
                               channel_list, &
                               cross_track_thin, &
                               along_track_thin, &
                               append_output, &
                               verbose

! ----------------------------------------------------------------------
! Declare local parameters
! ----------------------------------------------------------------------

integer                  :: io, iunit, ifile
integer                  :: filecount
integer                  :: max_num
integer                  :: num_inserted
type(amsua_bt_granule)   :: granule
type(obs_sequence_type)  :: big_sequence, small_sequence
logical                  :: use_channels(AMSUA_BT_CHANNEL) = .false.

! version controlled file description for error handling, do not edit
character(len=*), parameter :: source   = 'convert_amsu_L1.f90'

character(len=512) :: string1

! ----------------------------------------------------------------------
! start of executable program code
! ----------------------------------------------------------------------

call initialize_utilities('convert_amsu_L1')
call static_init_obs_sequence()

!----------------------------------------------------------------------
! Read the namelist
!----------------------------------------------------------------------

call find_namelist_in_file('input.nml', 'convert_amsu_L1_nml', iunit)
read(iunit, nml = convert_amsu_L1_nml, iostat = io)
call check_namelist_read(iunit, io, 'convert_amsu_L1_nml')

! Record the namelist values used for the run ...
if (do_nml_file()) write(nmlfileunit, nml=convert_amsu_L1_nml)
if (do_nml_term()) write(    *      , nml=convert_amsu_L1_nml)

! when this routine returns, the l1_files variable will have
! all the filenames, regardless of which way they were specified.
filecount = set_filename_list(l1_files, l1_file_list, "convert_amsu_L1")

! use the first file to initialize the module
call initialize_amsua_netcdf(l1_files(1),verbose)

! FIXME ... maybe ... are all the files guaranteed to have the
! same channels - I think so.
call channel_list_to_indices(channel_list, use_channels)

! Figure out how many (more) observations there may be in the output.
! It is impossible to predict how the lat/lon subsetting may affect this
! number, but the thinning and channels can be taken into account.

max_num = max_possible_obs(filecount, cross_track_thin, along_track_thin, use_channels)

! either read existing obs_seq or create a new one

call append_or_create(append_output,outputfile,max_num,big_sequence)

! read from netCDF file into a derived type that holds all the information
! convert derived type information to DART sequence for that file
! combine that sequence to the output sequence
! release the file sequence memory in preparation for creating it new next iteration

do ifile = 1,filecount

   call read_amsua_bt_netCDF_granule(l1_files(ifile), granule)

   call make_obs_sequence(small_sequence, granule, lon1, lon2, lat1, lat2, &
                          use_channels, along_track_thin, cross_track_thin)

   num_inserted = combine_sequences(small_sequence, big_sequence, l1_files(ifile))

   if (verbose > 2) then
      call print_obs_seq_summary(small_sequence)
   elseif (verbose > 1) then
      write(string1,*)trim(l1_files(ifile))//' added ',num_inserted,' observations.'
      call error_handler(E_MSG,string1,'')
   endif

   call destroy_obs_sequence(small_sequence)

enddo

! Report on the current status of the entire observation sequence
if (verbose > 0) then
   call print_obs_seq_summary(big_sequence)
endif

! DISCUSSION
! There is some discussion in the group about whether this program should
! error out or continue if there are no observations in the sequence.
! If used in a batch job, maybe you're OK with not having an output observation
! sequence for this collection of files -OR- you might want to know that there
! are no viable observations for a particular collection of files.
! If you want to allow this program to exit without returning an error code,
! change E_ERR to E_MSG in the following call to error_handler().

! Write the sequence to a disk file
if ( get_num_obs(big_sequence) > 0 )  then
   call write_obs_seq(big_sequence, outputfile)
else
   call error_handler(E_ERR,'NO OBSERVATIONS TO WRITE',source, &
              text2='see the DISCUSSION comment in the source to exit without error.')
endif

! release the sequence memory
call destroy_obs_sequence(big_sequence)

call error_handler(E_MSG, source, 'Finished successfully.',source)
call finalize_utilities()

end program convert_amsu_L1


text:
C***********************************************************************
C*
C*	Copyright 1996, University Corporation for Atmospheric Research.
C*
C*	psurf.f
C*
C*	Plots surface data. Derived from the GEMPAK program SFMAP.
C*
C*	History:
C*
C*	11/96	COMET		Original copy
C*	 2/97 J. Cowie/COMET	Add day of the week to title.
C*	 4/97	COMET		Added gprintf to support logging.
C*	 8/97	COMET		Added pmarker, marker.
C*	10/97	COMET		Added nstns, call to get_data_area()
C*	11/97	COMET		Added ptitle to display clickable titles.
C*	11/97	COMET		Added call to dscolor, cleaned up.
C*
C************************************************************************

	SUBROUTINE psurf ( area, sfparm, dattim, sffile, colors, title,
     +                     filter, text, msize, mwidth, marker, pmarker,
     +                     skpmis, nstns, ititle, verbose, iret )
C************************************************************************
C************************************************************************
	INCLUDE		'GEMPRM.PRM'
C*
	CHARACTER	area*(*), sfparm*(*),
     +			dattim*(*), sffile*(*), colors*(*),
     +			title*(*), filter*(*),
     +			text*(*), msize*(*), mwidth*(*)
    	INTEGER		verbose, pmarker, mrkwid

	LOGICAL		skpmis
C*
	CHARACTER	sffcur*72, arecur*48, fmt*80, arg*80
	CHARACTER 	pmdset (MMPARM)*4, prmlst (MMPARM)*4
	CHARACTER	tstn*8, sta*8, ttlstr*48, ttt*72
        CHARACTER       prcons (MMPARM)*16, chd (MMPARM)*8
	CHARACTER	area1*48, area2*48
	CHARACTER	blank*2
	INTEGER		icolor (MMPARM), iscale (MMPARM)
	LOGICAL		proces, newfil, chrflg (MMPARM)
	LOGICAL		wndflg, plot
	REAL		offset (4), sxplt (LLSTFL), outd (MMPARM)
	REAL		syplt (LLSTFL), data (MMPARM), size
	REAL		ccvals (MMPARM*LLCLEV), tsize(MMPARM)
	INTEGER         numccc (MMPARM), icclrs (MMPARM*LLCLEV)
	INTEGER         icrprm (MMPARM), jwide(MMPARM)
	CHARACTER	endflg*1

C------------------------------------------------------------------------
	iret = 0
	nstns = 0
	blank = ' ' // char(0)
C
C	Set up the initial conditions
C
	proces = .true.
	iflno  = 0
	iopt   = 0
	imark  = 0
	ioldclr = 0
	sffcur = ' '
	read ( msize, '(f5.1)' ) size
	read ( mwidth, '(i5)' ) mrkwid

	if ( verbose .gt. 0 ) call gfprints (
     +     'psurf'//char(0), blank )
	if ( verbose .gt. 1 ) then
	    call gfprints ( '  area = ' // char(0), area )
	    call gfprints ( '  sfparm = ' // char(0), sfparm )
	    call gfprints ( '  dattim = ' // char(0), dattim )
	    call gfprints ( '  sffile = ' // char(0), sffile )
	    call gfprints ( '  colors = ' // char(0), colors )
	    call gfprints ( '  title = ' // char(0), title )
	    call gfprints ( '  filter = ' // char(0), filter )
	    call gfprints ( '  text = ' // char(0), text )
	    call gfprints ( '  msize = ' // char(0), msize )
	    call gfprints ( '  mwidth = ' // char(0), mwidth )
	    call gfprinti ( '  marker = ' // char(0), marker )
	    call gfprinti ( '  pmarker = ' // char(0), pmarker )
	    call gfprinti ( '  skpmis = ' // char(0), skpmis )
	    call gfprinti ( '  ititle = ' // char(0), ititle )
	    call gfprinti ( '  verbose = ' // char(0), verbose )
	    call gfprinti ( '  iret = ' // char(0), iret )
	end if

C
        CALL IN_TEXT ( text, ier )
C
C*	Process filename.
C
	CALL SFMFIL  (  sffile, sffcur, iflno, newfil, pmdset,
     +			npmdst, iret )
	IF  ( iret .ne. 0 )  proces = .false.
C
C*          Process text, marker, title, filter and parms.
C
	IF  ( proces)  THEN
C
	        CALL GSMRKR ( marker, 0, size, mrkwid, ier )
		CALL IN_FILT ( filter, filtfc, ier )
C
C*		Process parameter names.
C
C		CALL SFMPRM  ( newfil, sfparm, pmdset, npmdst,
C    +			       prmlst, chrflg, ncprm, prcons,
C    +			       wndflg, ier )
		CALL SFMPRM ( sfparm, pmdset, npmdst, colors,
     +                        prmlst, chrflg, ncprm, prcons,
     +                        wndflg, icolor,
     +                        ccvals, icclrs, numccc, icrprm,
     +			      tsize, jwide,
     +                        iaddcl, endflg, ier )
C
C*		Determine whether any data will be plotted.
C
		IF   ( (ncprm .eq. 0) .and. (imark .eq. 0) ) THEN
	             plot = .false.
		  ELSE
	             plot = .true.
		END IF
C
C*		Translate input colors.
C
c		IF ( ncprm .gt. 0 )
c     +              CALL SFMCLR ( ncprm, prmlst, colors, icolor, ier )
		IF ( pmarker .gt. 0 )  imark = icolor(1)
C
C*              Get offsets for filtering.
C
		IF  ( ( filtfc .ne. 0. ) .and. plot )
     +             CALL SFMCOF ( ncprm - iaddcl, prmlst, wndflg,
     +			    	 filtfc, offset, ier )
C
C*		Take care of the special case of plotting a list of
C*		stations before plotting an area with the filter on.
C
		ipos2 = INDEX ( area, '/' )
		IF  ( area(1:1) .eq. '@' .and. ( ipos2 .gt. 4 ) ) THEN
		    area1 = area(:ipos2-1)
		    area2 = area(ipos2+1:)
		    iloop = 1
		  ELSE
		    area1 = area
		    iloop = 2
		END IF
C
C*		Convert the garea to a broader data area
C
		CALL get_data_area ( area1, iret )
C
		CALL LC_UARE  ( area1, newfil, iflno, arecur, tstn,
     +		  	        ier )
		IF  ( ier .ne. 0 )  proces = .false.
	END IF
C
C*	    Begin processing if inputs are ok.
C
		nplot = 0
	        CALL SF_STIM  ( iflno, dattim, ier )
C
C*	        Process clear, define panel, set up filtering
C
		IF  ( proces )  THEN
C
C*		   Intialize coordinate arrays for filtering.
C
	           IF  ( ( filtfc .ne. 0. ) .and. plot )  THEN
		       DO  m = 1, LLSTFL
			   sxplt (m) = RMISSD
		           syplt (m) = RMISSD
		       END DO
		   END IF
C
C*		   For special plotting, change the area on the
C*		   second time through.
C
		   DO  lll = iloop, 2
		     IF  ( ( lll .eq. 2 ) .and. ( iloop .eq. 1 ) )
     +							 THEN
		       CALL LC_UARE  ( area2, newfil, iflno,
     +		      		       arecur, tstn, ier )
		       IF  ( ier .ne. 0 )  plot = .false.
		     END IF
C
C*		     Station loop.
C
		     iout = 0
		     DO  WHILE  ( plot .and. (iout .eq. 0) )
			 CALL SF_SNXT ( iflno, sta, id, slat,
     +					slon, selv, ispri, iout )
	                 IF  ( iout .eq. 0 )  THEN
C
C*			     Get the data.
C
			     CALL SF_RDAT ( iflno, data, ihhmm, ier )
C
C*			     Increment the station counter
C
			     if ( ier .eq. 0 ) nstns = nstns + 1
C
C*			     Check for missing data and filter.
C
		             IF  ( ier .eq. 0 ) THEN
				 CALL PC_SSTN ( sta, id, slat, slon,
     +					selv, ispri, ihhmm, 1, ier )
				 CALL PC_CMVS ( 0., 0, data,
     +						outd, chd, ier )
			     END IF
C*
		             IF  ( ier .eq. 0 ) THEN

C*				Convert to plot coordinates.
C
                                CALL GTRANS  ( 'M', 'P', 1, slat, slon,
     +					        sx, sy, ier )
C
C*			        Filter, if requested.
C
			        IF  ( ( filtfc .ne. 0. ) .and.
     +				      ( lll .eq. 2 ) ) THEN
				   CALL SFMOVR  ( sx, sy, sxplt, syplt,
     +						  nplot, offset, ier )
C
C*				   Save x/y for no overlap.
C
				   IF  ( ier .eq. 0 )  THEN
				      nplot = nplot + 1
				      sxplt (nplot) = sx
				      syplt (nplot) = sy
				   END IF
				 ELSE IF  ( ( filtfc .ne. 0. ) ) THEN
				   nplot = nplot + 1
				   sxplt (nplot) = sx
				   syplt (nplot) = sy
		                 END IF
			     END IF
C
C*    	                     Plot if we are ok to here.
C
	                     IF  ( ier .eq. 0 )  THEN
c			       CALL SFMPLT  ( icolor, imark, sx, sy,
c     +					      slat, slon, chrflg,
c     +					      prmlst, ncprm, outd, chd,
c     +					      skpmis, ier )
				CALL SFMPLT  ( icolor, sx, sy, slat,
     +                             slon, chrflg, prmlst, ncprm, outd,
     +                             chd, ccvals, icclrs, numccc,
     +                             icrprm, tsize, jwide, endflg, ier)
			       IF ( pmarker .ne. 0 )
     +                         CALL GMARK ( 'P', 1, sx, sy, ier )
C
	                     END IF
			 END IF
	             END DO
C
C*		     Create the title string.
C
		     CALL IN_TITL ( title, -3, ititl, linttl,
     +				    ttlstr, ier )
		     DO ii = 1, ncprm
			iscale (ii) = 0
		     END DO
		     IF  ( ititl .gt. 0 )  THEN
			CALL GSCOLR  ( ititl, ier )
			CALL DSCOLR  ( ititl, ioldclr, iret )
			lens = LEN ( ttlstr )
			ttlstr(lens:lens) = char(0)
			call ptitle ( ttlstr, ititle )
		     END IF
C
		   END DO
C
C	Not processing but at least plot a title.
C
		ELSE
		     CALL IN_TITL ( title, -3, ititl, linttl,
     +				    ttlstr, ier )
		     IF  ( ititl .gt. 0 )  THEN
			CALL GSCOLR  ( ititl, ier )
			CALL DSCOLR  ( ititl, ioldclr, iret )
		        call ptitle ( ttlstr, ititle )
		     END IF
	     	END IF

C
C	End the plot to flush everything out
C
        CALL SF_CLOS(iflno,iret)
C	CALL GFLUSH ( iret )
C*
	if ( verbose .gt. 0 ) call gfprinti (
     +  'returning from psurf - iret = ' // char(0), iret )

	RETURN
	END

max_stars_repo_name: Unidata/Garp

max_issues_repo_name: Unidata/Garp

max_forks_repo_name: Unidata/Garp


text:
Copies a file with 5 comma-separated values to a line, appending a column holding their sum.
      INTEGER N			!Instead of littering the source with "5"
      PARAMETER (N = 5)		!Provide some provenance.
      CHARACTER*6 HEAD(N)	!A perfect size?
      INTEGER X(N)		!Integers suffice.
      INTEGER LINPR,IN		!I/O unit numbers.
      LINPR = 6		!Standard output via this unit number.
      IN = 10		!Some unit number for the input file.
      OPEN (IN,FILE="CSVtest.csv",STATUS="OLD",ACTION="READ")	!For formatted input.

      READ (IN,*) HEAD	!The first line has texts as column headings.
      WRITE (LINPR,1) (TRIM(HEAD(I)), I = 1,N),"Sum"	!Append a "Sum" column.
    1 FORMAT (666(A:","))	!The : sez "stop if no list element awaits".
    2 READ (IN,*,END = 10) X	!Read a line's worth of numbers, separated by commas or spaces.
      WRITE (LINPR,3) X,SUM(X)	!Write, with a total appended.
    3 FORMAT (666(I0:","))	!I0 editing uses only as many columns as are needed.
      GO TO 2			!Do it again.

   10 CLOSE (IN)	!All done.
      END	!That's all.

max_stars_repo_name: mullikine/RosettaCodeData

max_issues_repo_name: seanwallawalla-forks/RosettaCodeData

max_forks_repo_name: seanwallawalla-forks/RosettaCodeData


text:
Address(Sandpiper Drive) is a residential street in North Davis.

 Intersecting Streets

  Anderson Road
  Tern Place
  Pintail Place
  Curlew Street


text:
	!--------------------------------------------------------------------
	!
	!
	! PURPOSE
	!
	! This program use a monte carlo method to calculate pi
	!
	! .. Parameters ..
	!  npts				= total number of Monte Carlo points
	!  xmin				= lower bound for integration region
	!  xmax             = upper bound for integration region
	! .. Scalars ..
	!  i                = loop counter
	!  f				= average value from summation
	!  sum              = total sum
	!  randnum          = random number generated from (0,1) uniform
	!                     distribution
	!  x                = current Monte Carlo location
	! .. Arrays ..
	!
	! .. Vectors ..
	!
	! REFERENCES
	! http://chpc.wustl.edu/mpi-fortran.html
	! Gropp, Lusk and Skjellum, "Using MPI" MIT press (1999)
	!
	! ACKNOWLEDGEMENTS
	! The program below was modified from one available at the internet
	! address in the references. This internet address was last checked
	! on 30 March 2012
	!
	! ACCURACY
	!
	! ERROR INDICATORS AND WARNINGS
	!
	! FURTHER COMMENTS
	!
	!--------------------------------------------------------------------
	! External routines required
	!
	! External libraries required
	! None
 	PROGRAM monte_carlo
    IMPLICIT NONE

    INTEGER(kind=8), PARAMETER 	    :: npts = 1e10
    REAL(kind=8), PARAMETER 		:: xmin=0.0d0,xmax=1.0d0
    INTEGER(kind=8) 				:: i
    REAL(kind=8) 					:: f,sum, randnum,x

    DO i=1,npts
      CALL random_number(randnum)
      x = (xmax-xmin)*randnum + xmin
      sum = sum + 4.0d0/(1.0d0 + x**2)
    END DO
    f = sum/npts
    PRINT*,'PI calculated with ',npts,' points = ',f

    STOP
    END

max_stars_repo_name: bcloutier/PSNM

max_issues_repo_name: bcloutier/PSNM

max_forks_repo_name: bcloutier/PSNM


text:
#include "rundeck_opts.h"

!@sum  TRACER_COM: Exists alone to minimize the number of dependencies
!@+    This version for simple trace gases/chemistry/isotopes
!@auth Jean Lerner/Gavin Schmidt
!@ver  1.0

      MODULE TRACER_COM
!@sum  TRACER_COM tracer variables
!@auth Jean Lerner
!@ver  1.0
#ifdef TRACERS_ON
      USE QUSDEF, only: nmom
      USE MODEL_COM, only: im,jm,lm
#endif
      IMPLICIT NONE
      SAVE

C**** Each tracer has a variable name and a unique index
!@param NTM number of tracers
!@var TRNAME: Name for each tracer >>> MUST BE LEFT-JUSTIFIED <<<
#ifdef TRACERS_SPECIAL_O18
      integer, parameter :: ntm=3
      character*8, parameter :: trname(ntm)=(/
     *     'Water   ','H2O18   ','HDO     '/)  !,'HTO     '/)
#else
#ifdef TRACERS_SPECIAL_Lerner
      integer, parameter :: ntm=9
      character*8, parameter :: trname(ntm)= (/
     *     'Air     ','SF6     ','Rn222   ','CO2     ','N2O     ',
     *     'CFC11   ','14CO2   ','CH4     ','O3      '/)
#else
#ifdef TRACERS_SPECIAL_Shindell
      integer, parameter :: ntm=15
      character*8, parameter :: trname(ntm)=(/
     *    'Ox      ','NOx     ','N2O5    ','HNO3    ','H2O2    ',
     *    'CH3OOH  ','HCHO    ','HO2NO2  ','CO      ','CH4     ',
     *    'PAN     ','Isoprene','AlkylNit','Alkenes ','Paraffin'/)
#else
#ifdef TRACERS_WATER
      integer, parameter :: ntm=2
      character*8, parameter :: trname(ntm)=(/'Air     ','Water   '/)
#else
#ifdef TRACERS_OCEAN
      integer, parameter :: ntm=1
      character*8, parameter :: trname(ntm)=(/'Water   '/)
#else ! default for TRACERS_ON
      integer, parameter :: ntm=1
      character*8, parameter :: trname(ntm)=(/'Air     '/)
#endif
#endif
#endif
#endif
#endif
!@var N_XXX: variable names of indices for tracers
      integer ::
     *     n_Air,    n_SF6,   n_Rn222, n_CO2,      n_N2O,
     *     n_CFC11,  n_14CO2, n_CH4,   n_O3,       n_water,
     *     n_H2O18,  n_HDO,   n_HTO,   n_Ox,       n_NOx,
     *     n_N2O5,   n_HNO3,  n_H2O2,  n_CH3OOH,   n_HCHO,
     *     n_HO2NO2, n_CO,    n_PAN,   n_Isoprene, n_AlkylNit,
     *     n_Alkenes,n_Paraffin

C****    The following are set in tracer_IC
!@var T_QLIMIT: if t_qlimit=.true. tracer is maintained as positive
      logical, dimension(ntm) :: t_qlimit
!@var trdecay radioactive decay constant (1/s) (=0 for stable tracers)
      real*8, dimension(ntm) :: trdecay
!@dbparam ITIME_TR0: start time for each tracer (hours)
      integer, dimension(ntm) :: itime_tr0
!@var MTRACE: timing index for tracers
      integer mtrace

#ifdef TRACERS_ON
!@var NTM_POWER: Power of 10 associated with each tracer (for printing)
      integer, dimension(ntm) :: ntm_power
!@var TR_MM: molecular mass of each tracer (g/mole)
      real*8, dimension(ntm) :: tr_mm
!@var needtrs: true if surface tracer value from PBL is required
      logical, dimension(ntm) :: needtrs

C**** Note units for these parameters!
C**** Example: clay dust; trpdens=2.5d3, trradius=0.73d-6
C****          silt dust; trpdens=2.65d3, trradius=6.1d-6
C****
!@var trpdens tracer particle density (kg/m^3)
!@+               (=0 for non-particle tracers)
      real*8, dimension(ntm) :: trpdens
!@var trradius tracer effective radius (m) (=0 for non particle tracers)
      real*8, dimension(ntm) :: trradius

!@var TRM: Tracer array (kg)
      real*8, dimension(im,jm,lm,ntm) :: trm
!@var TRMOM: Second order moments for tracers (kg)
      real*8, dimension(nmom,im,jm,lm,ntm) :: trmom

!@var ntsurfsrcmax maximum number of surface 2D sources/sinks
      integer, parameter :: ntsurfsrcmax=14
!@var ntsurfsrc no. of non-interactive surface sources for each tracer
      integer, dimension(ntm) :: ntsurfsrc
!@var nt3Dsrcmax maximum number of 3D tracer sources/sinks
#ifdef TRACERS_SPECIAL_Shindell
      integer, parameter :: nt3Dsrcmax=4
#else
      integer, parameter :: nt3Dsrcmax=3
#endif
#endif

#ifdef TRACERS_WATER
!@param nWD_TYPES number of tracer types for wetdep purposes
!@param nGAS   index for wetdep tracer type = gas
!@param nPART  index for wetdep tracer type = particle/aerosol
!@param nWATER index for wetdep tracer type = water
!@+       (original condense method)
      integer, parameter :: nWD_TYPES=3, nGAS=1, nPART=2, nWATER=3
!@param tr_evap_fact fraction of re-evaporation by tracer type
C note, tr_evap_fact is not dimensioned as NTM:
      REAL*8, parameter, dimension(nWD_TYPES) :: tr_evap_fact=
     *     (/1.d0, 0.5d0,  1.d0/)

!@var tr_wd_TYPE: tracer wet dep type (gas, particle, water)
      integer, dimension(ntm) :: tr_wd_TYPE
!@var tr_RKD: Henry's Law coefficient (in mole/Joule please !)
      real*8, dimension(ntm) :: tr_RKD
!@var tr_DHD: coefficient of temperature-dependence term of Henry's
!@+   Law coefficient (in Joule/mole please !)
      real*8, dimension(ntm) :: tr_DHD
!@var tr_H2ObyCH4 conc. of tracer in water from methane oxidation
      real*8, dimension(ntm) :: tr_H2ObyCH4
!@var TRWM tracer in cloud liquid water amount (kg)
      real*8, dimension(im,jm,lm,ntm) :: trwm
#ifdef TRACERS_SPECIAL_O18
!@dbparam supsatfac factor controlling super saturation for isotopes
      real*8 :: supsatfac = 2d-3
#endif
#endif

#if (defined TRACERS_WATER) || (defined TRACERS_OCEAN)
!@var TRW0 default tracer concentration in water (kg/kg)
      real*8, dimension(ntm) :: trw0
!@var NTROCN scaling power factor for ocean/ice tracer concentrations
      integer, dimension(ntm) :: ntrocn
#endif

#ifdef TRACERS_OCEAN
!@var TRGLAC tracer ratio in glacial runoff to ocean (kg/kg)
      real*8, dimension(ntm) :: trglac
#endif

      END MODULE TRACER_COM