text: The account LPE was solely used to make several racial attacks on various users and Vent Magazine.  This account has been banned.  If you are the person behind LPE, I urge you to sleep on it and return when youre a bit more calm.

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
Image(shapiro.jpg, Stolen from dateline article, 300, right, thumbnail)Image(11305.jpg, Weather report from Jan 13, 2006., 200, right, thumbnail)Dr. Arthur Shapiro is a Faculty Professor of Evolution and Ecology at UC Davis. His wild black and Grey beard and mane of hair, unwashedlooking jeans, constantly faded and worn shirts, and leathery brown skin all make him look like a homeless guy  some even refer to him as the homeless professor.

Arthur Shapiro hangs out in the MU and can be seen walking all over Campus campus, sometimes catching Butterflies & Moths butterflies and moths.  He is one of the worlds foremost experts on butterflies, and sponsors the yearly Butterflies for Beer competition. He has the website http://butterfly.ucdavis.edu/ butterfly.ucdavis.edu.

Dr. Shapiro has been married since 1969 and has two children. He does not drive, and neither does anyone in his family. His wife is also brilliant and can be found at notable lectures whenever they happen around town. Their children, now grown and moved away from Davis, are also highly intelligent and knowledgeable. As for how he appears, he responds that I dress the way I do because I believe in the supremacy of content over display. Im also more comfortable. He considers appearance superficial and arbitrary; anyone can buy a suit at a department store and it doesnt change anything. Study hard and contribute to mankinds knowledge, now thats accomplishment worthy of respect.

Dr. Shapiro is also an amateur weather forecaster. His forecasts are hand written and posted on the door of the EVE department at 2320 Storer Hall. They are remarkable in detail and showcase a certain amount of weather wit. They may indeed be the best weather forecast for the Davis area.

 Art in the News
  http://wwwdateline.ucdavis.edu/011703/dl_shapiro.html Dateline
  http://www.knowledgeplex.org/news/78949.html From butterflies to beer Dont go by his looks: UC Davis profs an expert on lots more than insects (link broken, alternate source needed)
  http://www.news.ucdavis.edu/search/news_detail.lasso?id7324 On butterfly migration
  http://www.sacbee.com/content/news/story/12561209p13416298c.html Sac Bee (requires login)

Dr. Shapiro is definitely a genius.  I heard/watched him speak about ecological/evolutionary succession once; it was wonderful!  Users/CynthiaHe

20050923 22:35:24 nbsp Dr. Shapiro will be speaking at the Fall 2005 convocation Users/ArlenAbraham

Dr. Shapiro was the most animated speaker of the convocation, with an http://mypersonalcompass.ucdavis.edu/full_essay.php?id21 inspired essay on some bleeping thing. Users/TomHinds

20060820 11:13:57 nbsp I think he also goes to the Graduate on a fairly regular basis Users/MatthewTom

20110525 15:22:41 nbsp ask him about his spiritual beliefs :)
 Users/AnBan

20120413 13:43:34 nbsp If you like thinking about how small things in science tie into the greater everything, take one of his classes. If you like learning more than just whats written in the course guide, take a class. If you like good storytellers, take a class. Hes absolutely brilliant and is easily one of my favourite professors.  Users/JaniceW


text:
c	ELLIPTICITY CORRECTIONS - Polynomial interpolation
c
      subroutine ellip()
      real sc0
      real sc1
      real sc2
c SQUARE ROOT OF 3. TIMES .5
      real s3
c EPICENTAL CO-LATITUDE - RADIANS
      real ecolat
c EPICENTRAL DISTANCE - RADIANS
      real edist
c AZIMUTH FROM EPICENTER TO RECEIVER - RADIANS
      real az
c EPICENTRAL DEPTH (km)
      real edepth
c ELLIPTICITY CORRECTION -OUTPUT
      real tcor
      character*(*) phase
c TAU's
      real t0, t1, t2
c	CONSTANTS FOR POLYNOMIAL
      real t0con(8,10), t1con(8,10), t2con(8,10)
      integer ii, j
      real adepth
      data (t0con(1,j),j=1,10) /-0.01711,-1.7791,0.,0.,0.,-0.9630,-13.
     &2326,13.7390,0.,0./
      data (t0con(2,j),j=1,10) /-0.08291,-2.1455,2.4538,-0.7907,0.,2.
     &0258,-12.9357,2.1287,5.2668,-0.9229/
      data (t0con(3,j),j=1,10) /-1.5022,-0.0943,1.9655,-1.1661,.1393,3.
     &4920,-9.9051,-0.3875,5.3581,-0.0686/
      data (t0con(4,j),j=1,10) /2.9971,-2.9549,0.4082,0.,0.,28.1650,9.
     &2160,-17.9030,-5.2995,3.2029/
      data (t0con(5,j),j=1,10) /3.6775,-2.2221,0.,0.,0.,-1.3127,-6.2476,
     &1.6684,0.,0./
      data (t0con(6,j),j=1,10) /-10.6238,15.4993,-7.4840,1.0673,0.,3.
     &2763,-6.4596,-0.4923,0.,0./
      data (t0con(7,j),j=1,10) /-0.01332,-3.2777,-1.2243,7.5246,0.,-3.
     &4856,-10.3187,43.4834,-70.5341,-50.2287/
      data (t0con(8,j),j=1,10) /-0.07859,-4.0924,4.6116, -1.4760,0.,2.
     &9104,-17.8661, 4.6262,7.1486,-1.9154/
      data (t1con(1,j),j=1,10) /.0040,-0.7841,6.0441,-17.5535,0.,-0.
     &2549,2.0519,-19.0605,-37.8235,54.5110/
      data (t1con(2,j),j=1,10) /-.0048, .0839,-2.2705,2.4137,-0.5957,-2.
     &4241,-4.2792,1.9728,3.5644,-0.5285/
      data (t1con(3,j),j=1,10) /.0033,-1.3485,0.1735,1.1583,-0.4162,-0.
     &1096,0.2576,-0.5978,0.1888,0.1600/
      data (t1con(4,j),j=1,10) /2.6249,-.0025,-0.2086,-0.0184,0.,-1.
     &5077,0.9904,0.3513,0.,0./
      data (t1con(5,j),j=1,10) /3.4213,-0.9359,0.,0.,0.,0.,0.,0.,0.,0./
      data (t1con(6,j),j=1,10) /-8.0633,8.0238,-1.7407,0.,0.,0.,0.,0.,0.
     &,0./
      data (t1con(7,j),j=1,10) /0.0109,-1.2300,8.9145,-27.5847,0.,-0.
     &6951,5.6201,-33.0908,-83.8233,102.4333/
      data (t1con(8,j),j=1,10) /-0.0311,0.1896,-4.0694,4.2599,-1.0387,-
     &3.9368,-8.4379,2.6814,6.9535,-0.6086/
      data (t2con(1,j),j=1,10) /0.0107,0.0275,-0.6912,0.0347,0.1157,-0.
     &1836,0.,0.0296,0.,0./
      data (t2con(2,j),j=1,10) /0.0107,0.0275,-0.6912,0.0347,0.1157,-0.
     &1836,0.,0.0296,0.,0./
      data (t2con(3,j),j=1,10) /0.0005,-0.01231,-1.0156,0.4396,0.,0.,0.,
     &0.,0.,0./
      data (t2con(4,j),j=1,10) /-3.5838,2.1474,-0.3548,0.,0.,-1.3369,-5.
     &4889,0.6809,1.5096,-0.0763/
      data (t2con(5,j),j=1,10) /-2.9912,1.0313,0.,0.,0.,0.,0.,0.,0.,0./
      data (t2con(6,j),j=1,10) /3.2814,-7.1224,3.5418,-0.5115,0.,0.,0.,
     &0.,0.,0./
      data (t2con(7,j),j=1,10) /0.00025,0.1685,-2.2435,3.3433,0.,-0.
     &0503,0.5353,1.5362,-14.3118,-3.2938/
      data (t2con(8,j),j=1,10) /0.0843,-0.2917,-0.6767,-0.2934,0.2779,-
     &0.4336,0.0306,0.07113,0.,0./
      data s3 /.8660254/
      save sc0, sc1, sc2
c
c	INITIAL CALL TO SET UP CONSTANTS
      entry ellref(ecolat)
c
      sc0 = .25*(1.+3.*cos(2.*ecolat))
      sc1 = s3 * sin(2.*ecolat)
      sc2 = s3 * sin(ecolat)*sin(ecolat)
      return
c
c	CALLED ONCE FOR EACH STATION - RETURNS ELLIPTICITY CORRECTION IN tcor
      entry ellcor(edist, azim, edepth, phase,tcor)
      adepth = edepth/6371.
c DETERMINE INDEX FOR POLYNOMIAL
      if(.not.(phase .eq. 'P'))goto 23000
         if(.not.(edist .lt. (15.*(3.14159265/180.))))goto 23002
            ii = 1
            goto 23003
c        else
23002       continue
            ii = 2
23003    continue
         goto 23001
c     else
23000    continue
         if(.not.(phase .eq. 'PcP'))goto 23004
            ii = 3
            goto 23005
c        else
23004       continue
            if(.not.(phase .eq. 'PKPab'))goto 23006
               ii = 4
               goto 23007
c           else
23006          continue
               if(.not.(phase .eq. 'PKPbc'))goto 23008
                  ii = 5
                  goto 23009
c              else
23008             continue
                  if(.not.(phase .eq. 'PKIKP'))goto 23010
                     ii = 6
                     goto 23011
c                 else
23010                continue
                     if(.not.(phase .eq. 'S'))goto 23012
                        if(.not.(edist .lt. (15.*(3.14159265/180.))))
     &                    goto 23014
                           ii = 7
                           goto 23015
c                       else
23014                      continue
                           ii = 8
23015                   continue
                        goto 23013
c                    else
23012                   continue
                        tcor = 0.
                        return
23013                continue
23011             continue
23009          continue
23007       continue
23005    continue
23001 continue
c	COMPUTE TAU's
c      write(6,*) ii
      t0 = t0con(ii,1) + edist*(t0con(ii,2) + edist*(t0con(ii,3) +
     &edist*(t0con(ii,4) + edist*t0con(ii,5)))) +adepth*(t0con(ii,6) +
     &adepth*t0con(ii,7)) + adepth*edist*(t0con(ii,8) + t0con(ii,9)*
     &adepth + t0con(ii,10)*edist)
      t1 = t1con(ii,1) + edist*(t1con(ii,2) + edist*(t1con(ii,3) +
     &edist*(t1con(ii,4) + edist*t1con(ii,5)))) +adepth*(t1con(ii,6) +
     &adepth*t1con(ii,7)) + adepth*edist*(t1con(ii,8) + t1con(ii,9)*
     &adepth + t1con(ii,10)*edist)
      t2 = t2con(ii,1) + edist*(t2con(ii,2) + edist*(t2con(ii,3) +
     &edist*(t2con(ii,4) + edist*t2con(ii,5)))) +adepth*(t2con(ii,6) +
     &adepth*t2con(ii,7)) + adepth*edist*(t2con(ii,8) + t2con(ii,9)*
     &adepth + t2con(ii,10)*edist)
      tcor = sc0 * t0 + sc1 * cos(azim) * t1 + sc2 * cos(2.*azim) * t2
      return
      end

C ELLIPTICITY CORRECTIONS FOR AK135 MODEL (full set of phases)
C
       SUBROUTINE kellip()
C==========================================================================
C
C    Ellipticity correction for any given phase using
C    Dziewonski & Gilbert representation
C
C      The ellipticity corrections are found by linear interpolation
C    in terms of values calculated for the ak135 model for a wide
C    range of phases to match the output of the iasp software
C
C     first call:  ellref(ecolat)
C                        - to set up source dependent constants
C     2nd call  :  ellcor(phase,edist,depth,ecolat,azim,tcor,abrt)
C                        - to calculate correction for a station                                                                                C
C    Parameters:
C    character
C          phase : a  string specifying the PHASE,   -e.g P, ScP etc.
C
C    real
C          edist  :  epicentral distance to station (in degrees)
C          edepth :  depth of event
C          ecolat :  epicentral co-latitude of source (in radians)
C          azim   :  azimuth from source to station (in radians)
C
C          tcor   :  time correction for path to allow for ellipticity
C
C    logical
C          abrt   :  a logical variable -usally set to .FALSE.
C                    which is set to .TRUE. if a phase for
C                    which no data is available is chosen
C
C==========================================================================
C   B.L.N. Kennett RSES,ANU        May 1995, August 1996
C   (based on earlier routine by D.J. Brown)
C   with input from W. Spakman, Utrecht
C=========================================================================
      character *(*) phase
      character*8 phcod(57)
      integer phind(57),phspn(57),phnch(57)
      real edist,edepth,ecolat,azim,
     ^     sc0,sc1,sc2,s3,tcor,
     ^     tau0, a0,b0,h0,d0,e0,f0,g0,
     ^     tau1, a1,b1,h1,d1,e1,f1,g1,
     ^     tau2, a2,b2,h2,d2,e2,f2,g2
      real dpth(6),delta(50)
      real t0(50,6),t1(50,6),t2(50,6)
      integer Ne,Nd
      logical abrt
      data phcod/
     & "Pup   ","P     ","Pdiff ","PKPab ","PKPbc ","PKPdf ",
     & "PKiKP ","pP    ","pPKPab","pPKPbc","pPKPdf","pPKiKP",
     & "sP    ","sPKPab","sPKPbc","sPKPdf","sPKiKP","PcP   ",
     & "ScP   ","SKPab ","SKPbc ","SKPdf ","SKiKP ","PKKPab",
     & "PKKPbc","PKKPdf","SKKPab","SKKPbc","SKKPdf","PP    ",
     & "P'P'  ","Sup   ","S     ","Sdiff ","SKSac ","SKSdf ",
     & "pS    ","pSKSac","pSKSdf","sS    ","sSKSac","sSKSdf",
     & "ScS   ","PcS   ","PKSab ","PKSbc ","PKSdf ","PKKSab",
     & "PKKSbc","PKKSdf","SKKSac","SKKSdf","SS    ","S'S'  ",
     & "SP    ","PS    ","PnS   "/
      data phind/
     &        1,      14,      91,     136,     165,     178,
     &      235,     364,     433,     462,     475,     532,
     &      661,     742,     771,     784,     841,     970,
     &     1047,    1100,    1113,    1134,    1195,    1316,
     &     1337,    1382,    1507,    1516,    1573,    1702,
     &     1827,    1932,    1945,    2022,    2067,    2132,
     &     2197,    2234,    2295,    2356,    2425,    2490,
     &     2551,    2628,    2681,    2694,    2711,    2772,
     &     2781,    2838,    2967,    3140,    3273,    3398,
     &     3587,    3656,    3697/
      data phspn/
     &        3,      19,      11,       7,       3,      14,
     &       32,      17,       7,       3,      14,      32,
     &       20,       7,       3,      14,      32,      19,
     &       13,       3,       5,      15,      30,       5,
     &       11,      31,       2,      14,      32,      31,
     &       26,       3,      19,      11,      16,      16,
     &        9,      15,      15,      17,      16,      15,
     &       19,      13,       3,       4,      15,       2,
     &       14,      32,      43,      33,      31,      47,
     &       17,      10,       6/
      data phnch/
     &        3,       1,       5,       5,       5,       5,
     &        5,       2,       6,       6,       6,       6,
     &        2,       6,       6,       6,       6,       3,
     &        3,       5,       5,       5,       5,       6,
     &        6,       6,       6,       6,       6,       2,
     &        4,       3,       1,       5,       5,       5,
     &        2,       6,       6,       2,       6,       6,
     &        3,       3,       5,       5,       5,       6,
     &        6,       6,       6,       6,       2,       4,
     &        2,       2,       3/
      data dpth/ 0.0, 100.0, 200.0, 300.0, 500.0, 700.0 /
      save sc0,sc1,sc2
c...
c     In addition to the phase names listed above a number of phase aliases
c     are available in the routine phase_alias, e.g. Pn --> P etc
c     The input phase code is first checked against the phcod array
c     and next against the phase aliases.
c<sc>
c	                     initial call to set up source dependent constants
      entry kellref(ecolat)
c
      s3 = sqrt(3.0)/2.0
      sc0 = 0.25*(1.0+3.0*cos(2.0*ecolat))
      sc1 = s3*sin(2.0*ecolat)
      sc2 = s3*sin(ecolat)*sin(ecolat)
      return
c<sc>
c<ec>                                           phase identification
      entry kellcor(phase,edist,edepth,ecolat,azim,tcor,abrt)
*      write(6,*) "phase,edist,edepth,ecolat,azim"
*      write(6,*)  phase,edist,edepth,ecolat,azim
      Nd = 6
      NUMPH = 57
      deldst = 5.0
      abrt = .FALSE.
c                                             check on the length of phase
      l=len(phase)
      if(l.lt.8) then
       stop
     >    'character variable `phase` should have at least length 8'
      endif

c                                             select phase
      ip = -1
      nc=min(lnblk(phase),8)
      do 10 i=1,NUMPH
        if(nc.ne.phnch(i)) goto 10
        if (phase(1:nc) .eq. phcod(i)(1:nc)) then
          ip = i
          go to 11
        endif
 10   continue
 11   continue

      if(ip.eq.-1) then
c                                             check phase aliases
        call phase_alias(phase,edist,ip)
      endif
      Ne = phspn(ip)
*      write(6,*) "ip:",ip
c                                              phase not found
      if(ip.lt.0) then
        write(6,*) phase,'  is not available'
        abrt = .true.
        return
      endif
c                                               special case of upgoing waves
*
c
c                                                acquire phase information
       nr = phind(ip)
*       write(6,*) "nrec:",nr
       read(21,61,rec=nr) phcod(ip),np,d1,d2
*       write(6,*) "phcode,np,d1,d2: ", phcod(ip),np,d1,d2
       nr = nr+1
       if(np.ne.Ne) write(6,*) "HELP! - index wrong"
       do 15 i=1,np
         read(21,62,rec=nr) delta(i)
         nr = nr+1
         read(21,63,rec=nr) (t0(i,m),m=1,6)
         nr = nr+1
         read(21,63,rec=nr) (t1(i,m),m=1,6)
         nr = nr+1
         read(21,63,rec=nr) (t2(i,m),m=1,6)
         nr = nr+1
 15    continue
 61    format(a8,i10,2f10.0)
 62    format(f10.0)
 63    format(6f10.4)
c                                  distance index
       idist = 1 + int( (edist-d1)/ deldst )
       if(edist.lt.d1) idist =1
       if(edist.gt.d2) idist= np-1
c                                  depth index
       do 25 j = 1,Nd-1
         if ((dpth(j).le.edepth).and.(dpth(j+1).ge.edepth))then
            jdepth = j
            goto 26
         endif
 25    continue
 26    continue
*       write(6,*) "idist, jdepth;",idist,jdepth
c
*                      need to allow for zero entries (where phase
*                      description strongly depth dependent)
c tau0
         a0 = t0(idist,jdepth)
         b0 = t0(idist,jdepth+1)
         h0 = t0(idist+1,jdepth+1)
         d0 = t0(idist+1,jdepth)
         e0 = a0 +
     ^       (d0-a0)*(edist-delta(idist))/(delta(idist+1)-delta(idist))
         f0 = b0 +
     ^       (h0-b0)*(edist-delta(idist))/(delta(idist+1)-delta(idist))
         g0 = e0 +
     ^       (f0-e0)*(edepth-dpth(jdepth))/(dpth(jdepth+1)-dpth(jdepth))
         tau0 = g0
c tau1
         a1 = t1(idist,jdepth)
         b1 = t1(idist,jdepth+1)
         h1 = t1(idist+1,jdepth+1)
         d1 = t1(idist+1,jdepth)
         e1 = a1 +
     ^       (d1-a1)*(edist-delta(idist))/(delta(idist+1)-delta(idist))
         f1 = b1 +
     ^       (h1-b1)*(edist-delta(idist))/(delta(idist+1)-delta(idist))
         g1 = e1 +
     ^       (f1-e1)*(edepth-dpth(jdepth))/(dpth(jdepth+1)-dpth(jdepth))
         tau1 = g1
c tau2
         a2 = t2(idist,jdepth)
         b2 = t2(idist,jdepth+1)
         h2 = t2(idist+1,jdepth+1)
         d2 = t2(idist+1,jdepth)
         e2 = a2 +
     ^       (d2-a2)*(edist-delta(idist))/(delta(idist+1)-delta(idist))
         f2 = b2 +
     ^       (h2-b2)*(edist-delta(idist))/(delta(idist+1)-delta(idist))
         g2 = e2 +
     ^       (f2-e2)*(edepth-dpth(jdepth))/(dpth(jdepth+1)-dpth(jdepth))
         tau2 = g2
c
*         write(6,*) "tau0,tau1,tau2:",tau0,tau1,tau2
         caz = cos(azim)
         cbz = cos(2.0*azim)
*         write(6,*) "azim,caz,cbz",azim,caz,cbz
c
         tcor = sc0*tau0 + sc1*cos(azim)*tau1 + sc2*cos(2.0*azim)*tau2
c
      return
c<ec>
      end
      subroutine phase_alias(phase,delta,ip)

c-    check for alternative phase names
c     input phase, delta
c     output ip (index of phcod)

      character*(*) phase
      if(phase(1:3).eq.'Pg ') then
c       phase='P       '
        ip=2
      else if(phase(1:3).eq.'Sg ') then
c       phase='S       '
        ip=33
      else if(phase(1:4).eq.'pPg ') then
c       phase='pP      '
        ip=8
      else if(phase(1:4).eq.'sPg ') then
c       phase='sP      '
        ip=13
      else if(phase(1:4).eq.'pSg ') then
c       phase='pS      '
        ip=37
      else if(phase(1:4).eq.'sSg ') then
c       phase='sS      '
        ip=40
c
      elseif(phase(1:3).eq.'Pb ') then
c       phase='P       '
        ip=2
      else if(phase(1:3).eq.'Sb ') then
c       phase='S       '
        ip=33
      else if(phase(1:4).eq.'pPb ') then
c       phase='pP      '
        ip=8
      else if(phase(1:4).eq.'sPb ') then
c       phase='sP      '
        ip=13
      else if(phase(1:4).eq.'pSb ') then
c       phase='pS      '
        ip=37
      else if(phase(1:4).eq.'sSb ') then
c       phase='sS      '
c
      elseif(phase(1:3).eq.'Pn ') then
c       phase='P       '
        ip=2
      else if(phase(1:3).eq.'Sn ') then
c       phase='S       '
        ip=33
      else if(phase(1:4).eq.'pPn ') then
c       phase='pP      '
        ip=8
      else if(phase(1:4).eq.'sPn ') then
c       phase='sP      '
        ip=13
      else if(phase(1:4).eq.'pSn ') then
c       phase='pS      '
        ip=37
      else if(phase(1:4).eq.'sSn ') then
c       phase='sS      '
        ip=40
      else if(phase(1:4).eq.'SPn ') then
c       phase='SP      '
        ip=55
      else if(phase(1:4).eq.'SPb ') then
c       phase='SP      '
        ip=55
      else if(phase(1:4).eq.'SPg ') then
c       phase='SP      '
        ip=55
      else if(phase(1:4).eq.'SnP ') then
c       phase='SP      '
        ip=55
      else if(phase(1:4).eq.'PSn ') then
c       phase='PS      '
        ip=56
      else if(phase(1:5).eq.'PnPn ') then
c       phase='PP      '
        ip=30
      else if(phase(1:5).eq.'SnSn ') then
c       phase='SS      '
        ip=53
c                                       upgoing P, S
      else if(phase(1:2).eq.'p ') then
c       phase='Pup     '
        ip=1
      else if(phase(1:2).eq.'s ') then
c       phase='Sup     '
        ip=32
c
      else if(delta.le.100.0.and.phase.eq.'pPdiff  ') then
c       phase='pP      '
        ip=8
      else if(delta.le.100.0.and.phase.eq.'sPdiff  ') then
c       phase='sP      '
        ip=13
      else if(delta.le.100.0.and.phase.eq.'pSdiff  ') then
c       phase='pS      '
        ip=37
      else if(delta.le.100.0.and.phase.eq.'sSdiff  ') then
c       phase='sS      '
        ip=40
      else if(delta.ge.100.0.and.phase.eq.'pPdiff  ') then
c       phase='Pdiff   '
        ip=8
      else if(delta.ge.100.0.and.phase.eq.'sPdiff  ') then
c       phase='Pdiff   '
        ip=13
      else if(delta.ge.100.0.and.phase.eq.'pSdiff  ') then
c       phase='Sdiff   '
        ip=37
      else if(delta.ge.100.0.and.phase.eq.'sSdiff  ') then
c       phase='Sdiff    '
        ip=40
c
      else if(delta.le.165.0.and.phase.eq.'PKPdiff ') then
c       phase='PKPbc '
        ip=5
      else if(delta.le.165.0.and.phase.eq.'pPKPdiff') then
c       phase='pPKPbc '
        ip=10
      else if(delta.le.165.0.and.phase.eq.'sPKPdiff') then
c       phase='sPKPbc '
        ip=15
c
      else if(phase(1:4).eq."P'P'") then
c       phase="P'P'    "
        ip =31
      else if(phase(1:4).eq."S'S'") then
c       phase="P'P'    "
        ip =54
c                            diffractions (approx)
      else if(delta.gt.100.0.and.phase.eq.'pPdiff  ') then
c       phase='Pdiff   '
        ip=8
      else if(delta.gt.100.0.and.phase.eq.'sPdiff  ') then
c       phase='Pdiff   '
        ip=13
      else if(delta.gt.100.0.and.phase.eq.'pSdiff  ') then
c       phase='Sdiff   '
        ip=37
      else if(delta.gt.100.0.and.phase.eq.'sSdiff  ') then
c       phase='Sdiff    '
c
      else
        ip=-1
      endif
      return
      end
c-----
      integer function lnblk(s)
      character *(*) s
      l = len(s)
      do i=l,1,-1
        if (s(i:i).gt.' ') then
          lnblk=i
          return
        endif
      end do
      lnblk=0
      return
      end


text:
SUBROUTINE VEL_PREP
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
USE CLUSTER, only : natom, itype, ipas, vnu, tstep, tfin, irand_seed, g1,g2,g3!,thermostat_off
USE POTENTIAL, only : t2m
USE ENFORCE
!
IMPLICIT NONE
REAL :: dlaran
REAL :: qx(nsiz),qy(nsiz),qz(nsiz)
REAL :: ecinet,poguy,vargau,vel2
REAL :: vdfa
REAL :: vdfae
INTEGER :: nat3
INTEGER :: ir(4)
INTEGER :: i
!
vdfa=vnu*tstep+1.
IF(ipas==1)ir(:)=irand_seed(:)
!
     ecinet=0.d0
     poguy=0.d0
     nat3=natom
!
DO i=1,nat3

      qx(i)=(vx(i)+t2m(itype(i))*(fx(i)+dfx(i)))
      qy(i)=(vy(i)+t2m(itype(i))*(fy(i)+dfy(i)))
      qz(i)=(vz(i)+t2m(itype(i))*(fz(i)+dfz(i)))
!
!THERMALIZATION OF SYSTEM AT TINIT
!Through the Andersen Thermostat,
!Honeycutt, J. D., and H. C. Andersen, 1987, J. Phys, Chem. 91, 4950.
!Velocity is velocity per time step (in arete units!)
!
!if (thermostat_off) then !skip the thermostat if its turned off
!	!write(999,*) '################SKIPPING THERMOSTAT###############'
!	vx(i)=qx(i)
 !       vy(i)=qy(i)
  !      vz(i)=qz(i)
!else
  vdfae=dlaran(irand_seed)
     IF(vdfae > vdfa) THEN
!
           vx(i)=qx(i)
           vy(i)=qy(i)
           vz(i)=qz(i)
!
     ELSE
!
!CALL TO THE ANDERSEN THERMOSTAT
!
            vargau=SQRT(2.d0*t2m(itype(i))*cbol*tfin)
            CALL gauss

             vx(i)=g1*vargau
             vy(i)=g2*vargau
             vz(i)=g3*vargau
     ENDIF

!endif
!
!decision of thermostat call
!
! Kinetic energy calculation
! POGUY was introduced by Guy Treglia to speed up the system

         vel2=vx(i)*vx(i)+vy(i)*vy(i)+vz(i)*vz(i)
         poguy=vel2/(4.d0*t2m(itype(i)))
         ecinet=ecinet+poguy
!

ENDDO   !! LOOP on i

!Equipartition Theorem is applied
!
!
      ecin = ecinet
      temp=2.d0*ecin/(3.d0*(nat3)*cbol)
      etot=ener+ecin
!
END SUBROUTINE VEL_PREP

max_stars_repo_name: kcl-tscm/LoDiS

max_issues_repo_name: kcl-tscm/LoDiS

max_forks_repo_name: kcl-tscm/LoDiS


text:
      FUNCTION TWOT(OM,S,ZMS,TBM2)
C
C     'TWOT' IS A FUNCTION WHICH SUMMARIZES THE TOTAL KINETIC ENERGY
C     OF THE SYSTEM.
C
      IMPLICIT REAL*8 (A-H,O-Z)
      REAL*8 J1,J2,J3
C
      COMMON/RNRGY3/ J1(3,3),J2(3,3),J3(3,3),XICODT(3),W(3),ETA(3),
     .               ZETA(3)
C
      DIMENSION OM(3),S(3,3),OMT(3),ST(3,3),WT(3),OJ1(3),WJ2(3),WJ3(3)
C
C
      CALL MATRAN(OM,OMT,1,3)
      CALL MATRAN(S,ST,3,3)
      CALL MATRAN(W,WT,1,3)
C
      XSQR=0.0D0
      OMTE=0.0D0
      WTZ=0.0D0
C
      DO 10 I=1,3
      XSQR=XSQR + XICODT(I)**2
      OMTE=OMTE + OMT(I)*ETA(I)
   10 WTZ=WTZ + WT(I)*ZETA(I)
C
      XSQR=-ZMS*XSQR
      OMTE=2.D0*OMTE
      WTZ=2.D0*WTZ
C
      DO 20 I=1,3
      OJ1(I)=0.0D0
      WJ2(I)=0.0D0
      WJ3(I)=0.0D0
      DO 20 J=1,3
      OJ1(I)=OJ1(I) + OM(J)*J1(J,I)
      WJ2(I)=WJ2(I) + WT(J)*J2(J,I)
   20 WJ3(I)=WJ3(I) + WT(J)*J3(J,I)
C
      OJ1O=0.0D0
      WJ2W=0.0D0
      WJ3O=0.0D0
C
      DO 30 I=1,3
      OJ1O=OJ1O + OJ1(I)*OM(I)
      WJ2W=WJ2W + WJ2(I)*W(I)
   30 WJ3O=WJ3O + WJ3(I)*OM(I)
C
      TWOT=XSQR + OJ1O + WJ2W + TBM2 + 2.D0*WJ3O + OMTE + WTZ
      TWOT=TWOT/2.D0
C
      RETURN
      END

ext: for


! ==============================================================================
! Subroutine: TINT12 (LIST,EB)
!
! Purpose: Performs a table lookup for the internal energy for an asymmetric
!          mismatch pair between two basepairs in a helix.
!
! Method: Uses the MFOLD 3.0 energy function table for RNA @ T=37.
!
! Arguments:
!
!          LIST - Array of length 4 containing the nucleotides in
!                 numerical code (A=1,C=2,G=3,U=4) for the
!                 following locations:
!
!                         (3)
!                 5' (1) A .    X (6) 3'
!                 3' (2) U .  . Y (7) 5'
!                         (4)(5)
!
!                 where LIST(1) = letter code for position 1 etc.
!
!            EB - (OUTPUT) MFOLD 3.0 internal loop energy of the sequence
!                 provided in LIST.
!
! History:
!
! Version    Date         Comment
! --------   ----------   -----------------------
!            01/01/2015   Original Code
!
! Dependencies:
!
! Modules -
! Functions -
! Subroutines -
!
! Author(s): Eric Dykeman
!            Copyright (c) 2015 (Please Refer to LICENCE)
!
! ==============================================================================

      SUBROUTINE TINT12 (LIST,EB)

        IMPLICIT NONE

        !=== ARGUMENTS ===!

        INTEGER, INTENT(IN) :: list(7)
        REAL, INTENT(INOUT) :: eb

        !=== VARIABLES ===!

        INTEGER :: i,j,k,i1,i2,i3,i4,i5,i6,i7
        REAL :: au(16,24),cg(16,24),gc(16,24)
        REAL :: ua(16,24),gu(16,24),ug(16,24)

        DATA (au(1,i),i=1,24) / 3.90, 3.70, 3.10, 5.50, 3.20, 3.00, &
            & 2.40, 4.80, 3.20, 3.00, 2.40, 4.80, 3.90, 3.70, 3.10, &
            & 5.50, 3.90, 3.70, 3.10, 5.50, 3.90, 3.70, 3.10, 5.50 /
        DATA (au(2,i),i=1,24) / 3.80, 3.70, 5.50, 3.70, 3.10, 3.00, &
            & 4.80, 3.00, 3.10, 3.00, 4.80, 3.00, 3.80, 3.70, 5.50, &
            & 3.70, 3.80, 3.70, 5.50, 3.70, 3.80, 3.70, 5.50, 3.70 /
        DATA (au(3,i),i=1,24) / 3.20, 5.50, 2.30, 5.50, 2.50, 4.80, &
            & 1.60, 4.80, 2.50, 4.80, 1.60, 4.80, 3.20, 5.50, 2.30, &
            & 5.50, 3.20, 5.50, 2.30, 5.50, 3.20, 5.50, 2.30, 5.50 /
        DATA (au(4,i),i=1,24) / 5.50, 5.50, 5.50, 5.50, 4.80, 4.80, &
            & 4.80, 4.80, 4.80, 4.80, 4.80, 4.80, 5.50, 5.50, 5.50, &
            & 5.50, 5.50, 5.50, 5.50, 5.50, 5.50, 5.50, 5.50, 5.50 /
        DATA (au(5,i),i=1,24) / 3.60, 3.20, 3.10, 5.50, 2.90, 2.50, &
            & 2.40, 4.80, 2.90, 2.50, 2.40, 4.80, 3.60, 3.20, 3.10, &
            & 5.50, 3.60, 3.20, 3.10, 5.50, 3.60, 3.20, 3.10, 5.50 /
        DATA (au(6,i),i=1,24) / 3.70, 4.00, 5.50, 3.70, 3.00, 3.30, &
            & 4.80, 3.00, 3.00, 3.30, 4.80, 3.00, 3.70, 4.00, 5.50, &
            & 3.70, 3.70, 4.00, 5.50, 3.70, 3.70, 4.00, 5.50, 3.70 /
        DATA (au(7,i),i=1,24) / 5.50, 5.50, 5.50, 5.50, 4.80, 4.80, &
            & 4.80, 4.80, 4.80, 4.80, 4.80, 4.80, 5.50, 5.50, 5.50, &
            & 5.50, 5.50, 5.50, 5.50, 5.50, 5.50, 5.50, 5.50, 5.50 /
        DATA (au(8,i),i=1,24) / 5.50, 3.70, 5.50, 2.80, 4.80, 3.00, &
            & 4.80, 2.10, 4.80, 3.00, 4.80, 2.10, 5.50, 3.70, 5.50, &
            & 2.80, 5.50, 3.70, 5.50, 2.80, 5.50, 3.70, 5.50, 2.80 /
        DATA (au(9,i),i=1,24) / 2.50, 2.10, 1.90, 5.50, 1.80, 1.40, &
            & 1.20, 4.80, 1.80, 1.40, 1.20, 4.80, 2.50, 2.10, 1.90, &
            & 5.50, 2.50, 2.10, 1.90, 5.50, 2.50, 2.10, 1.90, 5.50 /
        DATA (au(10,i),i=1,24) / 5.50, 5.50, 5.50, 5.50, 4.80, 4.80, &
             & 4.80, 4.80, 4.80, 4.80, 4.80, 4.80, 5.50, 5.50, 5.50, &
             & 5.50, 5.50, 5.50, 5.50, 5.50, 5.50, 5.50, 5.50, 5.50 /
        DATA (au(11,i),i=1,24) / 2.30, 5.50, 3.70, 5.50, 1.60, 4.80, &
             & 3.00, 4.80, 1.60, 4.80, 3.00, 4.80, 2.30, 5.50, 3.70, &
             & 5.50, 2.30, 5.50, 3.70, 5.50, 2.30, 5.50, 3.70, 5.50 /
        DATA (au(12,i),i=1,24) / 5.50, 5.50, 5.50, 5.50, 4.80, 4.80, &
             & 4.80, 4.80, 4.80, 4.80, 4.80, 4.80, 5.50, 5.50, 5.50, &
             & 5.50, 5.50, 5.50, 5.50, 5.50, 5.50, 5.50, 5.50, 5.50 /
        DATA (au(13,i),i=1,24) / 5.50, 5.50, 5.50, 5.50, 4.80, 4.80, &
             & 4.80, 4.80, 4.80, 4.80, 4.80, 4.80, 5.50, 5.50, 5.50, &
             & 5.50, 5.50, 5.50, 5.50, 5.50, 5.50, 5.50, 5.50, 5.50 /
        DATA (au(14,i),i=1,24) / 4.00, 3.40, 5.50, 3.70, 3.30, 2.70, &
             & 4.80, 3.00, 3.30, 2.70, 4.80, 3.00, 4.00, 3.40, 5.50, &
             & 3.70, 4.00, 3.40, 5.50, 3.70, 4.00, 3.40, 5.50, 3.70 /
        DATA (au(15,i),i=1,24) / 5.50, 5.50, 5.50, 5.50, 4.80, 4.80, &
             & 4.80, 4.80, 4.80, 4.80, 4.80, 4.80, 5.50, 5.50, 5.50, &
             & 5.50, 5.50, 5.50, 5.50, 5.50, 5.50, 5.50, 5.50, 5.50 /
        DATA (au(16,i),i=1,24) / 5.50, 3.20, 5.50, 2.70, 4.80, 2.50, &
             & 4.80, 2.00, 4.80, 2.50, 4.80, 2.00, 5.50, 3.20, 5.50, &
             & 2.70, 5.50, 3.20, 5.50, 2.70, 5.50, 3.20, 5.50, 2.70 /

        DATA (cg(1,i),i=1,24) / 3.20, 3.00, 2.40, 4.80, 2.30, 2.20, &
            & 1.10, 4.00, 2.40, 2.20, 1.60, 4.00, 3.20, 3.00, 2.40, &
            & 4.80, 3.20, 3.00, 2.40, 4.80, 3.20, 3.00, 2.40, 4.80 /
        DATA (cg(2,i),i=1,24) / 3.10, 3.00, 4.80, 3.00, 2.30, 2.20, &
            & 4.00, 2.20, 2.30, 2.20, 4.00, 2.20, 3.10, 3.00, 4.80, &
            & 3.00, 3.10, 3.00, 4.80, 3.00, 3.10, 3.00, 4.80, 3.00 /
        DATA (cg(3,i),i=1,24) / 2.50, 4.80, 1.60, 4.80, 1.70, 4.00, &
            & 0.80, 4.00, 1.70, 4.00, 0.80, 4.00, 2.50, 4.80, 1.60, &
            & 4.80, 2.50, 4.80, 1.60, 4.80, 2.50, 4.80, 1.60, 4.80 /
        DATA (cg(4,i),i=1,24) / 4.80, 4.80, 4.80, 4.80, 4.00, 4.00, &
            & 4.00, 4.00, 4.00, 4.00, 4.00, 4.00, 4.80, 4.80, 4.80, &
            & 4.80, 4.80, 4.80, 4.80, 4.80, 4.80, 4.80, 4.80, 4.80 /
        DATA (cg(5,i),i=1,24) / 2.90, 2.50, 2.40, 4.80, 2.10, 1.70, &
            & 1.60, 4.00, 2.10, 1.70, 1.60, 4.00, 2.90, 2.50, 2.40, &
            & 4.80, 2.90, 2.50, 2.40, 4.80, 2.90, 2.50, 2.40, 4.80 /
        DATA (cg(6,i),i=1,24) / 3.00, 3.30, 4.80, 3.00, 2.20, 2.50, &
            & 4.00, 2.20, 2.20, 2.50, 4.00, 2.20, 3.00, 3.30, 4.80, &
            & 3.00, 3.00, 3.30, 4.80, 3.00, 3.00, 3.30, 4.80, 3.00 /
        DATA (cg(7,i),i=1,24) / 4.80, 4.80, 4.80, 4.80, 4.00, 4.00, &
            & 4.00, 4.00, 4.00, 4.00, 4.00, 4.00, 4.80, 4.80, 4.80, &
            & 4.80, 4.80, 4.80, 4.80, 4.80, 4.80, 4.80, 4.80, 4.80 /
        DATA (cg(8,i),i=1,24) / 4.80, 3.00, 4.80, 2.10, 4.00, 2.20, &
            & 4.00, 1.50, 4.00, 2.20, 4.00, 1.30, 4.80, 3.00, 4.80, &
            & 2.10, 4.80, 3.00, 4.80, 2.10, 4.80, 3.00, 4.80, 2.10 /
        DATA (cg(9,i),i=1,24) / 1.80, 1.40, 1.20, 4.80, 0.80, 0.60, &
            & 0.40, 4.00, 1.00, 0.60, 0.40, 4.00, 1.80, 1.40, 1.20, &
            & 4.80, 1.80, 1.40, 1.20, 4.80, 1.80, 1.40, 1.20, 4.80 /
        DATA (cg(10,i),i=1,24) / 4.80, 4.80, 4.80, 4.80, 4.00, 4.00, &
             & 4.00, 4.00, 4.00, 4.00, 4.00, 4.00, 4.80, 4.80, 4.80, &
             & 4.80, 4.80, 4.80, 4.80, 4.80, 4.80, 4.80, 4.80, 4.80 /
        DATA (cg(11,i),i=1,24) / 1.60, 4.80, 3.00, 4.80, 0.80, 4.00, &
             & 2.20, 4.00, 0.80, 4.00, 2.20, 4.00, 1.60, 4.80, 3.00, &
             & 4.80, 1.60, 4.80, 3.00, 4.80, 1.60, 4.80, 3.00, 4.80 /
        DATA (cg(12,i),i=1,24) / 4.80, 4.80, 4.80, 4.80, 4.00, 4.00, &
             & 4.00, 4.00, 4.00, 4.00, 4.00, 4.00, 4.80, 4.80, 4.80, &
             & 4.80, 4.80, 4.80, 4.80, 4.80, 4.80, 4.80, 4.80, 4.80 /
        DATA (cg(13,i),i=1,24) / 4.80, 4.80, 4.80, 4.80, 4.00, 4.00, &
             & 4.00, 4.00, 4.00, 4.00, 4.00, 4.00, 4.80, 4.80, 4.80, &
             & 4.80, 4.80, 4.80, 4.80, 4.80, 4.80, 4.80, 4.80, 4.80 /
        DATA (cg(14,i),i=1,24) / 3.30, 2.70, 4.80, 3.00, 2.50, 1.90, &
             & 4.00, 2.20, 2.50, 1.90, 4.00, 2.20, 3.30, 2.70, 4.80, &
             & 3.00, 3.30, 2.70, 4.80, 3.00, 3.30, 2.70, 4.80, 3.00 /
        DATA (cg(15,i),i=1,24) / 4.80, 4.80, 4.80, 4.80, 4.00, 4.00, &
             & 4.00, 4.00, 4.00, 4.00, 4.00, 4.00, 4.80, 4.80, 4.80, &
             & 4.80, 4.80, 4.80, 4.80, 4.80, 4.80, 4.80, 4.80, 4.80 /
        DATA (cg(16,i),i=1,24) / 4.80, 2.50, 4.80, 2.00, 4.00, 1.70, &
             & 4.00, 1.20, 4.00, 1.70, 4.00, 1.20, 4.80, 2.50, 4.80, &
             & 2.00, 4.80, 2.50, 4.80, 2.00, 4.80, 2.50, 4.80, 2.00 /

        DATA (gc(1,i),i=1,24) / 3.20, 3.00, 2.40, 4.80, 2.40, 2.20, &
            & 1.60, 4.00, 2.50, 2.20, 2.10, 4.00, 3.20, 3.00, 2.40, &
            & 4.80, 3.20, 3.00, 2.40, 4.80, 3.20, 3.00, 2.40, 4.80 /
        DATA (gc(2,i),i=1,24) / 3.10, 3.00, 4.80, 3.00, 2.30, 2.20, &
            & 4.00, 2.20, 2.30, 2.20, 4.00, 2.20, 3.10, 3.00, 4.80, &
            & 3.00, 3.10, 3.00, 4.80, 3.00, 3.10, 3.00, 4.80, 3.00 /
        DATA (gc(3,i),i=1,24) / 2.50, 4.80, 1.60, 4.80, 1.70, 4.00, &
            & 0.80, 4.00, 1.70, 4.00, 0.80, 4.00, 2.50, 4.80, 1.60, &
            & 4.80, 2.50, 4.80, 1.60, 4.80, 2.50, 4.80, 1.60, 4.80 /
        DATA (gc(4,i),i=1,24) / 4.80, 4.80, 4.80, 4.80, 4.00, 4.00, &
            & 4.00, 4.00, 4.00, 4.00, 4.00, 4.00, 4.80, 4.80, 4.80, &
            & 4.80, 4.80, 4.80, 4.80, 4.80, 4.80, 4.80, 4.80, 4.80 /
        DATA (gc(5,i),i=1,24) / 2.90, 2.50, 2.40, 4.80, 2.10, 1.70, &
            & 1.60, 4.00, 2.10, 1.70, 1.60, 4.00, 2.90, 2.50, 2.40, &
            & 4.80, 2.90, 2.50, 2.40, 4.80, 2.90, 2.50, 2.40, 4.80 /
        DATA (gc(6,i),i=1,24) / 3.00, 3.30, 4.80, 3.00, 2.20, 2.50, &
            & 4.00, 2.20, 2.20, 2.50, 4.00, 2.20, 3.00, 3.30, 4.80, &
            & 3.00, 3.00, 3.30, 4.80, 3.00, 3.00, 3.30, 4.80, 3.00 /
        DATA (gc(7,i),i=1,24) / 4.80, 4.80, 4.80, 4.80, 4.00, 4.00, &
            & 4.00, 4.00, 4.00, 4.00, 4.00, 4.00, 4.80, 4.80, 4.80, &
            & 4.80, 4.80, 4.80, 4.80, 4.80, 4.80, 4.80, 4.80, 4.80 /
        DATA (gc(8,i),i=1,24) / 4.80, 3.00, 4.80, 2.10, 4.00, 2.20, &
            & 4.00, 1.30, 4.00, 2.20, 4.00, 1.20, 4.80, 3.00, 4.80, &
            & 2.10, 4.80, 3.00, 4.80, 2.10, 4.80, 3.00, 4.80, 2.10 /
        DATA (gc(9,i),i=1,24) / 1.80, 1.40, 1.20, 4.80, 1.00, 0.60, &
            & 0.40, 4.00, 1.20, 0.60, 0.40, 4.00, 1.80, 1.40, 1.20, &
            & 4.80, 1.80, 1.40, 1.20, 4.80, 1.80, 1.40, 1.20, 4.80 /
        DATA (gc(10,i),i=1,24) / 4.80, 4.80, 4.80, 4.80, 4.00, 4.00, &
             & 4.00, 4.00, 4.00, 4.00, 4.00, 4.00, 4.80, 4.80, 4.80, &
             & 4.80, 4.80, 4.80, 4.80, 4.80, 4.80, 4.80, 4.80, 4.80 /
        DATA (gc(11,i),i=1,24) / 1.60, 4.80, 3.00, 4.80, 0.80, 4.00, &
             & 2.20, 4.00, 0.80, 4.00, 2.20, 4.00, 1.60, 4.80, 3.00, &
             & 4.80, 1.60, 4.80, 3.00, 4.80, 1.60, 4.80, 3.00, 4.80 /
        DATA (gc(12,i),i=1,24) / 4.80, 4.80, 4.80, 4.80, 4.00, 4.00, &
             & 4.00, 4.00, 4.00, 4.00, 4.00, 4.00, 4.80, 4.80, 4.80, &
             & 4.80, 4.80, 4.80, 4.80, 4.80, 4.80, 4.80, 4.80, 4.80 /
        DATA (gc(13,i),i=1,24) / 4.80, 4.80, 4.80, 4.80, 4.00, 4.00, &
             & 4.00, 4.00, 4.00, 4.00, 4.00, 4.00, 4.80, 4.80, 4.80, &
             & 4.80, 4.80, 4.80, 4.80, 4.80, 4.80, 4.80, 4.80, 4.80 /
        DATA (gc(14,i),i=1,24) / 3.30, 2.70, 4.80, 3.00, 2.50, 1.90, &
             & 4.00, 2.20, 2.50, 1.90, 4.00, 2.20, 3.30, 2.70, 4.80, &
             & 3.00, 3.30, 2.70, 4.80, 3.00, 3.30, 2.70, 4.80, 3.00 /
        DATA (gc(15,i),i=1,24) / 4.80, 4.80, 4.80, 4.80, 4.00, 4.00, &
             & 4.00, 4.00, 4.00, 4.00, 4.00, 4.00, 4.80, 4.80, 4.80, &
             & 4.80, 4.80, 4.80, 4.80, 4.80, 4.80, 4.80, 4.80, 4.80 /
        DATA (gc(16,i),i=1,24) / 4.80, 2.50, 4.80, 2.00, 4.00, 1.70, &
             & 4.00, 1.20, 4.00, 1.70, 4.00, 1.20, 4.80, 2.50, 4.80, &
             & 2.00, 4.80, 2.50, 4.80, 2.00, 4.80, 2.50, 4.80, 2.00 /

        DATA (ua(1,i),i=1,24) / 3.90, 3.70, 3.10, 5.50, 3.20, 3.00, &
            & 2.40, 4.80, 3.20, 3.00, 2.40, 4.80, 3.90, 3.70, 3.10, &
            & 5.50, 3.90, 3.70, 3.10, 5.50, 3.90, 3.70, 3.10, 5.50 /
        DATA (ua(2,i),i=1,24) / 3.80, 3.70, 5.50, 3.70, 3.10, 3.00, &
            & 4.80, 3.00, 3.10, 3.00, 4.80, 3.00, 3.80, 3.70, 5.50, &
            & 3.70, 3.80, 3.70, 5.50, 3.70, 3.80, 3.70, 5.50, 3.70 /
        DATA (ua(3,i),i=1,24) / 3.20, 5.50, 2.30, 5.50, 2.50, 4.80, &
            & 1.60, 4.80, 2.50, 4.80, 1.60, 4.80, 3.20, 5.50, 2.30, &
            & 5.50, 3.20, 5.50, 2.30, 5.50, 3.20, 5.50, 2.30, 5.50 /
        DATA (ua(4,i),i=1,24) / 5.50, 5.50, 5.50, 5.50, 4.80, 4.80, &
            & 4.80, 4.80, 4.80, 4.80, 4.80, 4.80, 5.50, 5.50, 5.50, &
            & 5.50, 5.50, 5.50, 5.50, 5.50, 5.50, 5.50, 5.50, 5.50 /
        DATA (ua(5,i),i=1,24) / 3.60, 3.20, 3.10, 5.50, 2.90, 2.50, &
            & 2.40, 4.80, 2.90, 2.50, 2.40, 4.80, 3.60, 3.20, 3.10, &
            & 5.50, 3.60, 3.20, 3.10, 5.50, 3.60, 3.20, 3.10, 5.50 /
        DATA (ua(6,i),i=1,24) / 3.70, 4.00, 5.50, 3.70, 3.00, 3.30, &
            & 4.80, 3.00, 3.00, 3.30, 4.80, 3.00, 3.70, 4.00, 5.50, &
            & 3.70, 3.70, 4.00, 5.50, 3.70, 3.70, 4.00, 5.50, 3.70 /
        DATA (ua(7,i),i=1,24) / 5.50, 5.50, 5.50, 5.50, 4.80, 4.80, &
            & 4.80, 4.80, 4.80, 4.80, 4.80, 4.80, 5.50, 5.50, 5.50, &
            & 5.50, 5.50, 5.50, 5.50, 5.50, 5.50, 5.50, 5.50, 5.50 /
        DATA (ua(8,i),i=1,24) / 5.50, 3.70, 5.50, 2.80, 4.80, 3.00, &
            & 4.80, 2.10, 4.80, 3.00, 4.80, 2.10, 5.50, 3.70, 5.50, &
            & 2.80, 5.50, 3.70, 5.50, 2.80, 5.50, 3.70, 5.50, 2.80 /
        DATA (ua(9,i),i=1,24) / 2.50, 2.10, 1.90, 5.50, 1.80, 1.40, &
            & 1.20, 4.80, 1.80, 1.40, 1.20, 4.80, 2.50, 2.10, 1.90, &
            & 5.50, 2.50, 2.10, 1.90, 5.50, 2.50, 2.10, 1.90, 5.50 /
        DATA (ua(10,i),i=1,24) / 5.50, 5.50, 5.50, 5.50, 4.80, 4.80, &
             & 4.80, 4.80, 4.80, 4.80, 4.80, 4.80, 5.50, 5.50, 5.50, &
             & 5.50, 5.50, 5.50, 5.50, 5.50, 5.50, 5.50, 5.50, 5.50 /
        DATA (ua(11,i),i=1,24) / 2.30, 5.50, 3.70, 5.50, 1.60, 4.80, &
             & 3.00, 4.80, 1.60, 4.80, 3.00, 4.80, 2.30, 5.50, 3.70, &
             & 5.50, 2.30, 5.50, 3.70, 5.50, 2.30, 5.50, 3.70, 5.50 /
        DATA (ua(12,i),i=1,24) / 5.50, 5.50, 5.50, 5.50, 4.80, 4.80, &
             & 4.80, 4.80, 4.80, 4.80, 4.80, 4.80, 5.50, 5.50, 5.50, &
             & 5.50, 5.50, 5.50, 5.50, 5.50, 5.50, 5.50, 5.50, 5.50 /
        DATA (ua(13,i),i=1,24) / 5.50, 5.50, 5.50, 5.50, 4.80, 4.80, &
             & 4.80, 4.80, 4.80, 4.80, 4.80, 4.80, 5.50, 5.50, 5.50, &
             & 5.50, 5.50, 5.50, 5.50, 5.50, 5.50, 5.50, 5.50, 5.50 /
        DATA (ua(14,i),i=1,24) / 4.00, 3.40, 5.50, 3.70, 3.30, 2.70, &
             & 4.80, 3.00, 3.30, 2.70, 4.80, 3.00, 4.00, 3.40, 5.50, &
             & 3.70, 4.00, 3.40, 5.50, 3.70, 4.00, 3.40, 5.50, 3.70 /
        DATA (ua(15,i),i=1,24) / 5.50, 5.50, 5.50, 5.50, 4.80, 4.80, &
             & 4.80, 4.80, 4.80, 4.80, 4.80, 4.80, 5.50, 5.50, 5.50, &
             & 5.50, 5.50, 5.50, 5.50, 5.50, 5.50, 5.50, 5.50, 5.50 /
        DATA (ua(16,i),i=1,24) / 5.50, 3.20, 5.50, 2.70, 4.80, 2.50, &
             & 4.80, 2.00, 4.80, 2.50, 4.80, 2.00, 5.50, 3.20, 5.50, &
             & 2.70, 5.50, 3.20, 5.50, 2.70, 5.50, 3.20, 5.50, 2.70 /

        DATA (gu(1,i),i=1,24) / 3.90, 3.70, 3.10, 5.50, 3.20, 3.00, &
            & 2.40, 4.80, 3.20, 3.00, 2.40, 4.80, 3.90, 3.70, 3.10, &
            & 5.50, 3.90, 3.70, 3.10, 5.50, 3.90, 3.70, 3.10, 5.50 /
        DATA (gu(2,i),i=1,24) / 3.80, 3.70, 5.50, 3.70, 3.10, 3.00, &
            & 4.80, 3.00, 3.10, 3.00, 4.80, 3.00, 3.80, 3.70, 5.50, &
            & 3.70, 3.80, 3.70, 5.50, 3.70, 3.80, 3.70, 5.50, 3.70 /
        DATA (gu(3,i),i=1,24) / 3.20, 5.50, 2.30, 5.50, 2.50, 4.80, &
            & 1.60, 4.80, 2.50, 4.80, 1.60, 4.80, 3.20, 5.50, 2.30, &
            & 5.50, 3.20, 5.50, 2.30, 5.50, 3.20, 5.50, 2.30, 5.50 /
        DATA (gu(4,i),i=1,24) / 5.50, 5.50, 5.50, 5.50, 4.80, 4.80, &
            & 4.80, 4.80, 4.80, 4.80, 4.80, 4.80, 5.50, 5.50, 5.50, &
            & 5.50, 5.50, 5.50, 5.50, 5.50, 5.50, 5.50, 5.50, 5.50 /
        DATA (gu(5,i),i=1,24) / 3.60, 3.20, 3.10, 5.50, 2.90, 2.50, &
            & 2.40, 4.80, 2.90, 2.50, 2.40, 4.80, 3.60, 3.20, 3.10, &
            & 5.50, 3.60, 3.20, 3.10, 5.50, 3.60, 3.20, 3.10, 5.50 /
        DATA (gu(6,i),i=1,24) / 3.70, 4.00, 5.50, 3.70, 3.00, 3.30, &
            & 4.80, 3.00, 3.00, 3.30, 4.80, 3.00, 3.70, 4.00, 5.50, &
            & 3.70, 3.70, 4.00, 5.50, 3.70, 3.70, 4.00, 5.50, 3.70 /
        DATA (gu(7,i),i=1,24) / 5.50, 5.50, 5.50, 5.50, 4.80, 4.80, &
            & 4.80, 4.80, 4.80, 4.80, 4.80, 4.80, 5.50, 5.50, 5.50, &
            & 5.50, 5.50, 5.50, 5.50, 5.50, 5.50, 5.50, 5.50, 5.50 /
        DATA (gu(8,i),i=1,24) / 5.50, 3.70, 5.50, 2.80, 4.80, 3.00, &
            & 4.80, 2.10, 4.80, 3.00, 4.80, 2.10, 5.50, 3.70, 5.50, &
            & 2.80, 5.50, 3.70, 5.50, 2.80, 5.50, 3.70, 5.50, 2.80 /
        DATA (gu(9,i),i=1,24) / 2.50, 2.10, 1.90, 5.50, 1.80, 1.40, &
            & 1.20, 4.80, 1.80, 1.40, 1.20, 4.80, 2.50, 2.10, 1.90, &
            & 5.50, 2.50, 2.10, 1.90, 5.50, 2.50, 2.10, 1.90, 5.50 /
        DATA (gu(10,i),i=1,24) / 5.50, 5.50, 5.50, 5.50, 4.80, 4.80, &
             & 4.80, 4.80, 4.80, 4.80, 4.80, 4.80, 5.50, 5.50, 5.50, &
             & 5.50, 5.50, 5.50, 5.50, 5.50, 5.50, 5.50, 5.50, 5.50 /
        DATA (gu(11,i),i=1,24) / 2.30, 5.50, 3.70, 5.50, 1.60, 4.80, &
             & 3.00, 4.80, 1.60, 4.80, 3.00, 4.80, 2.30, 5.50, 3.70, &
             & 5.50, 2.30, 5.50, 3.70, 5.50, 2.30, 5.50, 3.70, 5.50 /
        DATA (gu(12,i),i=1,24) / 5.50, 5.50, 5.50, 5.50, 4.80, 4.80, &
             & 4.80, 4.80, 4.80, 4.80, 4.80, 4.80, 5.50, 5.50, 5.50, &
             & 5.50, 5.50, 5.50, 5.50, 5.50, 5.50, 5.50, 5.50, 5.50 /
        DATA (gu(13,i),i=1,24) / 5.50, 5.50, 5.50, 5.50, 4.80, 4.80, &
             & 4.80, 4.80, 4.80, 4.80, 4.80, 4.80, 5.50, 5.50, 5.50, &
             & 5.50, 5.50, 5.50, 5.50, 5.50, 5.50, 5.50, 5.50, 5.50 /
        DATA (gu(14,i),i=1,24) / 4.00, 3.40, 5.50, 3.70, 3.30, 2.70, &
             & 4.80, 3.00, 3.30, 2.70, 4.80, 3.00, 4.00, 3.40, 5.50, &
             & 3.70, 4.00, 3.40, 5.50, 3.70, 4.00, 3.40, 5.50, 3.70 /
        DATA (gu(15,i),i=1,24) / 5.50, 5.50, 5.50, 5.50, 4.80, 4.80, &
             & 4.80, 4.80, 4.80, 4.80, 4.80, 4.80, 5.50, 5.50, 5.50, &
             & 5.50, 5.50, 5.50, 5.50, 5.50, 5.50, 5.50, 5.50, 5.50 /
        DATA (gu(16,i),i=1,24) / 5.50, 3.20, 5.50, 2.70, 4.80, 2.50, &
             & 4.80, 2.00, 4.80, 2.50, 4.80, 2.00, 5.50, 3.20, 5.50, &
             & 2.70, 5.50, 3.20, 5.50, 2.70, 5.50, 3.20, 5.50, 2.70 /

        DATA (ug(1,i),i=1,24) / 3.90, 3.70, 3.10, 5.50, 3.20, 3.00, &
            & 2.40, 4.80, 3.20, 3.00, 2.40, 4.80, 3.90, 3.70, 3.10, &
            & 5.50, 3.90, 3.70, 3.10, 5.50, 3.90, 3.70, 3.10, 5.50 /
        DATA (ug(2,i),i=1,24) / 3.80, 3.70, 5.50, 3.70, 3.10, 3.00, &
            & 4.80, 3.00, 3.10, 3.00, 4.80, 3.00, 3.80, 3.70, 5.50, &
            & 3.70, 3.80, 3.70, 5.50, 3.70, 3.80, 3.70, 5.50, 3.70 /
        DATA (ug(3,i),i=1,24) / 3.20, 5.50, 2.30, 5.50, 2.50, 4.80, &
            & 1.60, 4.80, 2.50, 4.80, 1.60, 4.80, 3.20, 5.50, 2.30, &
            & 5.50, 3.20, 5.50, 2.30, 5.50, 3.20, 5.50, 2.30, 5.50 /
        DATA (ug(4,i),i=1,24) / 5.50, 5.50, 5.50, 5.50, 4.80, 4.80, &
            & 4.80, 4.80, 4.80, 4.80, 4.80, 4.80, 5.50, 5.50, 5.50, &
            & 5.50, 5.50, 5.50, 5.50, 5.50, 5.50, 5.50, 5.50, 5.50 /
        DATA (ug(5,i),i=1,24) / 3.60, 3.20, 3.10, 5.50, 2.90, 2.50, &
            & 2.40, 4.80, 2.90, 2.50, 2.40, 4.80, 3.60, 3.20, 3.10, &
            & 5.50, 3.60, 3.20, 3.10, 5.50, 3.60, 3.20, 3.10, 5.50 /
        DATA (ug(6,i),i=1,24) / 3.70, 4.00, 5.50, 3.70, 3.00, 3.30, &
            & 4.80, 3.00, 3.00, 3.30, 4.80, 3.00, 3.70, 4.00, 5.50, &
            & 3.70, 3.70, 4.00, 5.50, 3.70, 3.70, 4.00, 5.50, 3.70 /
        DATA (ug(7,i),i=1,24) / 5.50, 5.50, 5.50, 5.50, 4.80, 4.80, &
            & 4.80, 4.80, 4.80, 4.80, 4.80, 4.80, 5.50, 5.50, 5.50, &
            & 5.50, 5.50, 5.50, 5.50, 5.50, 5.50, 5.50, 5.50, 5.50 /
        DATA (ug(8,i),i=1,24) / 5.50, 3.70, 5.50, 2.80, 4.80, 3.00, &
            & 4.80, 2.10, 4.80, 3.00, 4.80, 2.10, 5.50, 3.70, 5.50, &
            & 2.80, 5.50, 3.70, 5.50, 2.80, 5.50, 3.70, 5.50, 2.80 /
        DATA (ug(9,i),i=1,24) / 2.50, 2.10, 1.90, 5.50, 1.80, 1.40, &
            & 1.20, 4.80, 1.80, 1.40, 1.20, 4.80, 2.50, 2.10, 1.90, &
            & 5.50, 2.50, 2.10, 1.90, 5.50, 2.50, 2.10, 1.90, 5.50 /
        DATA (ug(10,i),i=1,24) / 5.50, 5.50, 5.50, 5.50, 4.80, 4.80, &
             & 4.80, 4.80, 4.80, 4.80, 4.80, 4.80, 5.50, 5.50, 5.50, &
             & 5.50, 5.50, 5.50, 5.50, 5.50, 5.50, 5.50, 5.50, 5.50 /
        DATA (ug(11,i),i=1,24) / 2.30, 5.50, 3.70, 5.50, 1.60, 4.80, &
             & 3.00, 4.80, 1.60, 4.80, 3.00, 4.80, 2.30, 5.50, 3.70, &
             & 5.50, 2.30, 5.50, 3.70, 5.50, 2.30, 5.50, 3.70, 5.50 /
        DATA (ug(12,i),i=1,24) / 5.50, 5.50, 5.50, 5.50, 4.80, 4.80, &
             & 4.80, 4.80, 4.80, 4.80, 4.80, 4.80, 5.50, 5.50, 5.50, &
             & 5.50, 5.50, 5.50, 5.50, 5.50, 5.50, 5.50, 5.50, 5.50 /
        DATA (ug(13,i),i=1,24) / 5.50, 5.50, 5.50, 5.50, 4.80, 4.80, &
             & 4.80, 4.80, 4.80, 4.80, 4.80, 4.80, 5.50, 5.50, 5.50, &
             & 5.50, 5.50, 5.50, 5.50, 5.50, 5.50, 5.50, 5.50, 5.50 /
        DATA (ug(14,i),i=1,24) / 4.00, 3.40, 5.50, 3.70, 3.30, 2.70, &
             & 4.80, 3.00, 3.30, 2.70, 4.80, 3.00, 4.00, 3.40, 5.50, &
             & 3.70, 4.00, 3.40, 5.50, 3.70, 4.00, 3.40, 5.50, 3.70 /
        DATA (ug(15,i),i=1,24) / 5.50, 5.50, 5.50, 5.50, 4.80, 4.80, &
             & 4.80, 4.80, 4.80, 4.80, 4.80, 4.80, 5.50, 5.50, 5.50, &
             & 5.50, 5.50, 5.50, 5.50, 5.50, 5.50, 5.50, 5.50, 5.50 /
        DATA (ug(16,i),i=1,24) / 5.50, 3.20, 5.50, 2.70, 4.80, 2.50, &
             & 4.80, 2.00, 4.80, 2.50, 4.80, 2.00, 5.50, 3.20, 5.50, &
             & 2.70, 5.50, 3.20, 5.50, 2.70, 5.50, 3.20, 5.50, 2.70 /

        !         (3)            !
        ! 5' (1) A .    X (6) 3' !
        ! 3' (2) U .  . Y (7) 5' !
        !         (4)(5)         !

        i1 = list(1)
        i2 = list(2)
        i3 = list(3)
        i4 = list(4)
        i5 = list(5)
        i6 = list(6)
        i7 = list(7)

        IF ( i5 == 1 ) j = 0 + i3
        IF ( i5 == 2 ) j = 4 + i3
        IF ( i5 == 3 ) j = 8 + i3
        IF ( i5 == 4 ) j =12 + i3

        IF ( i6 == 1 .and. i7 == 4 ) k = 0 + i4
        IF ( i6 == 2 .and. i7 == 3 ) k = 4 + i4
        IF ( i6 == 3 .and. i7 == 2 ) k = 8 + i4
        IF ( i6 == 4 .and. i7 == 1 ) k =12 + i4
        IF ( i6 == 3 .and. i7 == 4 ) k =16 + i4
        IF ( i6 == 4 .and. i7 == 3 ) k =20 + i4

        IF ( i1 == 1 .and. i2 == 4 ) eb = eb + au(j,k)
        IF ( i1 == 2 .and. i2 == 3 ) eb = eb + cg(j,k)
        IF ( i1 == 3 .and. i2 == 2 ) eb = eb + gc(j,k)
        IF ( i1 == 4 .and. i2 == 1 ) eb = eb + ua(j,k)
        IF ( i1 == 3 .and. i2 == 4 ) eb = eb + gu(j,k)
        IF ( i1 == 4 .and. i2 == 3 ) eb = eb + ug(j,k)

        RETURN

      END SUBROUTINE TINT12


text:
!> @ingroup
!>
!> @brief
!> Run  on a vector of  values.
!>
!> @details
!> Run  on a vector of  values to sort them in a specified order.
!>
!> @param
!> The vector to be sorted.
!>
!> @param ascending
!> Indicates if the ordering is ascending or decending.
!>
!> | Value | Meaning         |
!> |:-----:| :-------------- |
!> | true  | Ascending order |
!> | false | Decending order |
!>
!> @author Arin R. Bratt
!> @date 11/07/2014

pure subroutine (,ascending)
  ! Arguments
  , intent(inout) :: (:)
  logical, intent(in) :: ascending

  ! Main
  if (ascending) then
     call _asc()
  else
     call _dec()
  end if
end subroutine

max_stars_repo_name: arinrb/bsort

max_issues_repo_name: arinrb/bsort

max_forks_repo_name: arinrb/bsort


text:
C
C	$Id: gactm.f,v 1.5 2008-07-27 00:20:56 haley Exp $
C
C                Copyright (C)  2000
C        University Corporation for Atmospheric Research
C                All Rights Reserved
C
C The use of this Software is governed by a License Agreement.
C
      SUBROUTINE GACTM(MINP,X0,Y0,DX,DY,PHI,SX,SY,SW,MOUT)
C
C  ACCUMULATE TRANSFORMATION MATRIX
C
      INTEGER EACTM
      PARAMETER (EACTM=106)
C
      REAL X0,Y0,DX,DY,PHI,MINP(2,3),MOUT(2,3)
      INTEGER SW
C
C  Check if GKS is in the proper state.
C
      CALL GZCKST(8,EACTM,IER)
      IF (IER .NE. 0) RETURN
C
C  The transformation takes place in NDC space, so do the
C  appropriate conversions.
C
      IF (SW .EQ. 0) THEN
        CALL GZW2NX(1,X0,XN0)
        CALL GZW2NY(1,Y0,YN0)
        CALL GZW2NX(1,DX,DTX)
        CALL GZW2NX(1,DY,DTY)
        CALL GZW2NX(1,0.,D0X)
        CALL GZW2NX(1,0.,D0Y)
        DNX = DTX - D0X
        DNY = DTY - D0Y
      ELSE
        XN0 = X0
        YN0 = Y0
        DNX = DX
        DNY = DY
      ENDIF
C
C  Calculate the new transformtion matrix.
C
      IF (PHI .EQ. 0.) THEN
        MOUT(1,1) = SX*MINP(1,1)
        MOUT(1,2) = SX*MINP(1,2)
        MOUT(1,3) = SX*MINP(1,3)+DNX+XN0-XN0*SX
        MOUT(2,1) = SY*MINP(2,1)
        MOUT(2,2) = SY*MINP(2,2)
        MOUT(2,3) = SY*MINP(2,3)+DNY+YN0-YN0*SY
        GO TO 10
      ELSE
        CPHI = COS(PHI)
        SPHI = SIN(PHI)
C
        T11 = SX*CPHI
        T12 = -SY*SPHI
        T13 = DNX+XN0-XN0*SX*CPHI+YN0*SY*SPHI
        T21 = SX*SPHI
        T22 = SY*CPHI
        T23 = DNY+YN0-XN0*SX*SPHI-YN0*SY*CPHI
        MOUT(1,1) =  T11*MINP(1,1)+T12*MINP(2,1)
        MOUT(1,2) =  T11*MINP(1,2)+T12*MINP(2,2)
        MOUT(1,3) =  T11*MINP(1,3)+T12*MINP(2,3)+T13
        MOUT(2,1) =  T21*MINP(1,1)+T22*MINP(2,1)
        MOUT(2,2) =  T21*MINP(1,2)+T22*MINP(2,2)
        MOUT(2,3) =  T21*MINP(1,3)+T22*MINP(2,3)+T23
        GO TO 10
      ENDIF
C
   10 CONTINUE
      RETURN
      END

max_stars_repo_name: tenomoto/ncl

max_issues_repo_name: tenomoto/ncl

max_forks_repo_name: tenomoto/ncl


text:
! Associate creates an alias for expressions or variables.
program assoc
implicit none
integer, allocatable :: u(:)
u = [4,9]
associate (v => u)
v = v - 2
print*,v ! 2 7
print*,u ! 2 7
associate (w => [4,9])
print*,w
! w = w - 2 ! illegal since w set to a constant expression
end associate
end associate
end program assoc

max_stars_repo_name: awvwgk/FortranTip

max_issues_repo_name: awvwgk/FortranTip

max_forks_repo_name: awvwgk/FortranTip


text:
c       this file contains the following user-callable routines:
c
c
c       routine idzr_svd computes an approximation of specified rank
c       to a given matrix, in the usual SVD form U S V^*,
c       where U has orthonormal columns, V has orthonormal columns,
c       and S is diagonal.
c
c       routine idzp_svd computes an approximation of specified
c       precision to a given matrix, in the usual SVD form U S V^*,
c       where U has orthonormal columns, V has orthonormal columns,
c       and S is diagonal.
c
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c
c
c
        subroutine idzr_svd(m,n,a,krank,u,v,s,ier,r)
c
c       constructs a rank-krank SVD  u diag(s) v^*  approximating a,
c       where u is an m x krank matrix whose columns are orthonormal,
c       v is an n x krank matrix whose columns are orthonormal,
c       and diag(s) is a diagonal krank x krank matrix whose entries
c       are all nonnegative. This routine combines a QR code
c       (which is based on plane/Householder reflections)
c       with the LAPACK routine zgesdd.
c
c       input:
c       m -- first dimension of a and u
c       n -- second dimension of a, and first dimension of v
c       a -- matrix to be SVD'd
c       krank -- desired rank of the approximation to a
c
c       output:
c       u -- left singular vectors of a corresponding
c            to the k greatest singular values of a
c       v -- right singular vectors of a corresponding
c            to the k greatest singular values of a
c       s -- k greatest singular values of a
c       ier -- 0 when the routine terminates successfully;
c              nonzero when the routine encounters an error
c
c       work:
c       r -- must be at least
c            (krank+2)*n+8*min(m,n)+6*krank**2+8*krank
c            complex*16 elements long
c
c       _N.B._: This routine destroys a. Also, please beware that
c               the source code for this routine could be clearer.
c
        implicit none
        character*1 jobz
        integer m,n,k,krank,ifadjoint,ldr,ldu,ldvadj,lwork,
     1          info,j,ier,io
        real*8 s(krank)
        complex*16 a(m,n),u(m,krank),v(n*krank),r(*)
c
c
        io = 8*min(m,n)
c
c
        ier = 0
c
c
c       Compute a pivoted QR decomposition of a.
c
        call idzr_qrpiv(m,n,a,krank,r,r(io+1))
c
c
c       Extract R from the QR decomposition.
c
        call idz_retriever(m,n,a,krank,r(io+1))
c
c
c       Rearrange R according to ind.
c
        call idz_permuter(krank,r,krank,n,r(io+1))
c
c
c       Use LAPACK to SVD r,
c       storing the krank (krank x 1) left singular vectors
c       in r(io+krank*n+1 : io+krank*n+krank*krank).
c
        jobz = 'S'
        ldr = krank
        lwork = 2*(krank**2+2*krank+n)
        ldu = krank
        ldvadj = krank
c
        call zgesdd(jobz,krank,n,r(io+1),ldr,s,r(io+krank*n+1),ldu,
     1              v,ldvadj,r(io+krank*n+krank*krank+1),lwork,
     2              r(io+krank*n+krank*krank+lwork+1),r,info)
c
        if(info .ne. 0) then
          ier = info
          return
        endif
c
c
c       Multiply the U from R from the left by Q to obtain the U
c       for A.
c
        do k = 1,krank
c
          do j = 1,krank
            u(j,k) = r(io+krank*n+j+krank*(k-1))
          enddo ! j
c
          do j = krank+1,m
            u(j,k) = 0
          enddo ! j
c
        enddo ! k
c
        ifadjoint = 0
        call idz_qmatmat(ifadjoint,m,n,a,krank,krank,u,r)
c
c
c       Take the adjoint of v to obtain r.
c
        call idz_adjer(krank,n,v,r)
c
c
c       Copy r into v.
c
        do k = 1,n*krank
          v(k) = r(k)
        enddo ! k
c
c
        return
        end
c
c
c
c
        subroutine idzp_svd(lw,eps,m,n,a,krank,iu,iv,is,w,ier)
c
c       constructs a rank-krank SVD  U Sigma V^*  approximating a
c       to precision eps, where U is an m x krank matrix whose
c       columns are orthonormal, V is an n x krank matrix whose
c       columns are orthonormal, and Sigma is a diagonal krank x krank
c       matrix whose entries are all nonnegative.
c       The entries of U are stored in w, starting at w(iu);
c       the entries of V are stored in w, starting at w(iv).
c       The diagonal entries of Sigma are stored in w,
c       starting at w(is). This routine combines a QR code
c       (which is based on plane/Householder reflections)
c       with the LAPACK routine zgesdd.
c
c       input:
c       lw -- maximum usable length of w (in complex*16 elements)
c       eps -- precision to which the SVD approximates a
c       m -- first dimension of a and u
c       n -- second dimension of a, and first dimension of v
c       a -- matrix to be SVD'd
c
c       output:
c       krank -- rank of the approximation to a
c       iu -- index in w of the first entry of the matrix
c             of orthonormal left singular vectors of a
c       iv -- index in w of the first entry of the matrix
c             of orthonormal right singular vectors of a
c       is -- index in w of the first entry of the array
c             of singular values of a; the singular values are stored
c             as complex*16 numbers whose imaginary parts are zeros
c       w -- array containing the singular values and singular vectors
c            of a; w doubles as a work array, and so must be at least
c            (krank+1)*(m+2*n+9)+8*min(m,n)+6*krank**2
c            complex*16 elements long, where krank is the rank
c            output by the present routine
c       ier -- 0 when the routine terminates successfully;
c              -1000 when lw is too small;
c              other nonzero values when zgesdd bombs
c
c       _N.B._: This routine destroys a. Also, please beware that
c               the source code for this routine could be clearer.
c               w must be at least
c               (krank+1)*(m+2*n+9)+8*min(m,n)+6*krank**2
c               complex*16 elements long, where krank is the rank
c               output by the present routine.
c
        implicit none
        character*1 jobz
        integer m,n,k,krank,ifadjoint,ldr,ldu,ldvadj,lwork,
     1          info,j,ier,io,iu,iv,is,ivi,isi,lu,lv,ls,lw
        real*8 eps
        complex*16 a(m,n),w(*)
c
c
        io = 8*min(m,n)
c
c
        ier = 0
c
c
c       Compute a pivoted QR decomposition of a.
c
        call idzp_qrpiv(eps,m,n,a,krank,w,w(io+1))
c
c
        if(krank .gt. 0) then
c
c
c         Extract R from the QR decomposition.
c
          call idz_retriever(m,n,a,krank,w(io+1))
c
c
c         Rearrange R according to ind.
c
          call idz_permuter(krank,w,krank,n,w(io+1))
c
c
c         Use LAPACK to SVD R,
c         storing the krank (krank x 1) left singular vectors
c         in w(io+krank*n+1 : io+krank*n+krank*krank).
c
          jobz = 'S'
          ldr = krank
          lwork = 2*(krank**2+2*krank+n)
          ldu = krank
          ldvadj = krank
c
          ivi = io+krank*n+krank*krank+lwork+3*krank**2+4*krank+1
          lv = n*krank
c
          isi = ivi+lv
          ls = krank
c
          if(lw .lt. isi+ls+m*krank-1) then
            ier = -1000
            return
          endif
c
          call zgesdd(jobz,krank,n,w(io+1),ldr,w(isi),w(io+krank*n+1),
     1                ldu,w(ivi),ldvadj,w(io+krank*n+krank*krank+1),
     2                lwork,w(io+krank*n+krank*krank+lwork+1),w,info)
c
          if(info .ne. 0) then
            ier = info
            return
          endif
c
c
c         Take the adjoint of w(ivi:ivi+lv-1) to obtain V.
c
          iv = 1
          call idz_adjer(krank,n,w(ivi),w(iv))
c
c
c         Copy w(isi:isi+ls/2) into w(is:is+ls-1).
c
          is = iv+lv
c
          call idz_realcomp(ls,w(isi),w(is))
c
c
c         Multiply the U from R from the left by Q to obtain the U
c         for A.
c
          iu = is+ls
          lu = m*krank
c
          do k = 1,krank
c
            do j = 1,krank
              w(iu-1+j+krank*(k-1)) = w(io+krank*n+j+krank*(k-1))
            enddo ! j
c
          enddo ! k
c
          do k = krank,1,-1
c
            do j = m,krank+1,-1
              w(iu-1+j+m*(k-1)) = 0
            enddo ! j
c
            do j = krank,1,-1
              w(iu-1+j+m*(k-1)) = w(iu-1+j+krank*(k-1))
            enddo ! j
c
          enddo ! k
c
          ifadjoint = 0
          call idz_qmatmat(ifadjoint,m,n,a,krank,krank,w(iu),
     1                     w(iu+lu+1))
c
c
        endif ! krank .gt. 0
c
c
        return
        end
c
c
c
c
        subroutine idz_realcomp(n,a,b)
c
c       copies the real*8 array a into the complex*16 array b.
c
c       input:
c       n -- length of a and b
c       a -- real*8 array to be copied into b
c
c       output:
c       b -- complex*16 copy of a
c
        integer n,k
        real*8 a(n)
        complex*16 b(n)
c
c
        do k = 1,n
          b(k) = a(k)
        enddo ! k
c
c
        return
        end
c
c
c
c
        subroutine idz_permuter(krank,ind,m,n,a)
c
c       permutes the columns of a according to ind obtained
c       from routine idzr_qrpiv or idzp_qrpiv, assuming that
c       a = q r from idzr_qrpiv or idzp_qrpiv.
c
c       input:
c       krank -- rank specified to routine idzr_qrpiv
c                or obtained from routine idzp_qrpiv
c       ind -- indexing array obtained from routine idzr_qrpiv
c              or idzp_qrpiv
c       m -- first dimension of a
c       n -- second dimension of a
c       a -- matrix to be rearranged
c
c       output:
c       a -- rearranged matrix
c
        implicit none
        integer k,krank,m,n,j,ind(krank)
        complex*16 cswap,a(m,n)
c
c
        do k = krank,1,-1
          do j = 1,m
c
            cswap = a(j,k)
            a(j,k) = a(j,ind(k))
            a(j,ind(k)) = cswap
c
          enddo ! j
        enddo ! k
c
c
        return
        end
c
c
c
c
        subroutine idz_retriever(m,n,a,krank,r)
c
c       extracts R in the QR decomposition specified by the output a
c       of the routine idzr_qrpiv or idzp_qrpiv
c
c       input:
c       m -- first dimension of a
c       n -- second dimension of a and r
c       a -- output of routine idzr_qrpiv or idzp_qrpiv
c       krank -- rank specified to routine idzr_qrpiv,
c                or output by routine idzp_qrpiv
c
c       output:
c       r -- triangular factor in the QR decomposition specified
c            by the output a of the routine idzr_qrpiv or idzp_qrpiv
c
        implicit none
        integer m,n,j,k,krank
        complex*16 a(m,n),r(krank,n)
c
c
c       Copy a into r and zero out the appropriate
c       Householder vectors that are stored in one triangle of a.
c
        do k = 1,n
          do j = 1,krank
            r(j,k) = a(j,k)
          enddo ! j
        enddo ! k
c
        do k = 1,n
          if(k .lt. krank) then
            do j = k+1,krank
              r(j,k) = 0
            enddo ! j
          endif
        enddo ! k
c
c
        return
        end
c
c
c
c
        subroutine idz_adjer(m,n,a,aa)
c
c       forms the adjoint aa of a.
c
c       input:
c       m -- first dimension of a and second dimension of aa
c       n -- second dimension of a and first dimension of aa
c       a -- matrix whose adjoint is to be taken
c
c       output:
c       aa -- adjoint of a
c
        implicit none
        integer m,n,j,k
        complex*16 a(m,n),aa(n,m)
c
c
        do k = 1,n
          do j = 1,m
            aa(k,j) = conjg(a(j,k))
          enddo ! j
        enddo ! k
c
c
        return
        end

max_stars_repo_name: Ennosigaeon/scipy

max_issues_repo_name: Ennosigaeon/scipy

max_forks_repo_name: Ennosigaeon/scipy


!---------------------------------------------------------------------!
! OWNER: Ithaca Combustion Enterprise, LLC                            !
! COPYRIGHT: © 2012, Ithaca Combustion Enterprise, LLC                !
! LICENSE: BSD 3-Clause License (The complete text of the license can !
!  be found in the `LICENSE-ICE.txt' file included in the ISAT-CK7    !
!  source directory.)                                                 !
!---------------------------------------------------------------------!

subroutine ell_rad_lower( n, gg, r )

!  Determine the radius r of the ball inscribed in the ellipsoid E given
!  by { x | norm(G^T * x) <=1 ), where G is an n x n lower triangular
!  matrix.  The array gg contains the matrix G in packed format.

!  Method:  r = 1 / sqrt(lam_max)  where lam_max is the largest
!           eigenvalue of G * G^T

!  S.B. Pope 6/12/04

implicit none

integer, parameter      :: k_dp = kind(1.d0)
integer, intent(in)     :: n
real(k_dp), intent(in)  :: gg((n*(n+1))/2)
real(k_dp), intent(out) :: r

integer    :: info, i, j, k, nev, lwork, iwork(5*n), ifail(n)
real(k_dp) :: g(n,n), a(n,n), dum(1), ev(n), dumz(n,1), work(n*(n+8))

lwork = n*(n+8)

!  unpack gg into g
k = 0
g = 0.d0
do j = 1, n
   do i = j, n
      k = k + 1
	  g(i,j) = gg(k)
   end do
end do

!  form a = A = G * G^T
a = g
call dtrmm ( 'R', 'L', 'T', 'N', n, n, 1.d0, g, n, a, n )

!  find first (smallest) eigenvalues of A

call dsyevx( 'N', 'I', 'L', n, a, n, dum, dum, n, n, 0.d0, nev, ev, dumz, n, &
             work, lwork, iwork, ifail, info )

if( info /= 0 ) then  !  use SVD instead
   call ell_chol2eig( n, gg, a, work )
   r = 1.d0 / work(1)
   return
endif

if( ev(1) <= 0.d0 ) then
   write(0,*)'ell_rad_lower: lam_min <=0 ', ev(1)
   stop
endif

r = 1.d0 / sqrt( ev(1) )

return
end subroutine ell_rad_lower


program simple_simd

    use helpers
    use functions

    use, intrinsic :: iso_fortran_env

    implicit none

    ! General variables
    type(t_arguments)                                       :: args ! CLI arguments
    integer(INT8), dimension(:, :), allocatable, target     :: field_one, field_two ! Cell data array
    integer(INT8), dimension(:, :), pointer                 :: field_current, field_temp ! Cell data pointers
    real(REAL64)                                            :: time_start, time_finish, time_delta, &
                                                                time_sum, clock_delta, clock_sum ! Timing stamps
    integer(INT64)                                          :: clock_start, clock_finish, clock_rate = 0 ! Wallclock

    integer                                                 :: alloc_stat_one, alloc_stat_two ! Cell array allocation status
    integer                                                 :: k ! Step index

    ! Algorithmus specific variables
    integer(INT8)                                           :: cell_sum
    integer(INT32)                                          :: i, j

    write(*, "(A)") "Program: Simple optimized further: Vectorized sums using SIMD"

    ! Parse CLI arguments
    call arguments_get(args, .true.)

    write(*, "(A)") "Initializing..."
    call system_clock(clock_start, clock_rate)
    call cpu_time(time_start)

    ! Allocate cell data array
    ! Note that we allocate double the memory, this is used to compute the next state while still reading the old
    ! Note that we allocate a border of one cell, this is to capture any outflow
    ! By default, arrays are 1-indexed in fortran, so we use index 0 and n+1
    allocate(field_one(0:args%height + 1, 0:args%width + 1), stat = alloc_stat_one)
    if (alloc_stat_one .ne. 0) then
        write(*, "(A, I0)") "Error: Cannot allocate field_one memory: ", alloc_stat_one
        stop 1
    end if
    allocate(field_two(0:args%height + 1, 0:args%width + 1), stat = alloc_stat_two)
    if (alloc_stat_two .ne. 0) then
        write(*, "(A, I0)") "Error: Cannot allocate field_two memory: ", alloc_stat_two
        stop 1
    end if

    ! Initialize data
    field_one = 0
    field_two = 0
    field_current => field_one
    field_temp => field_two
    ! Initialize cells randomly
    write(*, "(A)") "Generating..."
    ! call field_pattern(field_current)
    call field_randomize(field_current, args%width, args%height)

    call cpu_time(time_finish)
    call system_clock(clock_finish)
    time_delta = time_finish - time_start
    clock_delta = real(clock_finish - clock_start) / real(clock_rate)
    ! Print initialization diagnostics
    call print_init_report(args, time_delta, clock_delta, field_current, .false.)

    ! Main computation loop
    write(*, "(A, I0, A)") "Computing ", args%steps, " steps..."
    do k = 1, args%steps
        ! We donst swap around anymore, we use the other field for intermediate storage

        call system_clock(clock_start)
        call cpu_time(time_start)

        ! Naive implementation with lookup and row sums
        ! We sum three colums next to each other row-wise and store them
        ! This loop can be vectorized using SIMD
        do i = 1, args%width
            field_temp(:, i) = sum(field_current(:, i-1:i+1), dim=2)
        end do

        ! We iterate column-wise to exploit CPU cache locality,
        ! because fortran lays out its array memory column-wise.
        do i = 1, args%width
            do j = 1, args%height
                ! Because we have a outflow border, we do not have to worry about edge cases
                ! Substract center cell, we only want neighbours
                ! We now only have to sum 3 elements in the current column, further reducing cache misses
                ! Also, this summation can be vectorized using SIMD
                cell_sum = sum(field_temp(j-1:j+1, i), dim=1) - field_current(j, i)

                ! We decide on the next state of this cell based on the count of neighbours
                ! Instead of explicit comparisions, we look up the new state in a lookup table
                ! Note that we write back to the same field
                field_current(j, i) = neighbour_lookup(cell_sum, field_current(j, i))
            end do
        end do

        call cpu_time(time_finish)
        call system_clock(clock_finish)
        time_delta = time_finish - time_start
        time_sum = time_sum + time_delta
        clock_delta = real(clock_finish - clock_start) / real(clock_rate)
        clock_sum = clock_sum + clock_delta
        ! Print step diagnostics
        call print_step_report(args, time_delta, clock_delta, k, field_current, .false.)
    end do

    ! Print concluding diagnostics
    call print_report(args, time_sum, clock_sum, "simple_simd")
    deallocate(field_one)
    deallocate(field_two)
end


text:
      subroutine MASKVALUE(
     &           mask
     &           ,imasklo0
     &           ,imaskhi0
     &           ,test
     &           ,itestlo0
     &           ,itesthi0
     &           ,ibxlo0
     &           ,ibxhi0
     &           ,val
     &           ,onoff
     &           )
      implicit none
      integer*8 ch_flops
      COMMON/ch_timer/ ch_flops
      integer imasklo0
      integer imaskhi0
      REAL*8 mask(
     &           imasklo0:imaskhi0)
      integer itestlo0
      integer itesthi0
      integer test(
     &           itestlo0:itesthi0)
      integer ibxlo0
      integer ibxhi0
      integer val
      integer onoff
      integer i0
      if (onoff .eq. 1) then
      do i0 = ibxlo0,ibxhi0
         if (test(i0) .eq. val) then
            mask(i0) = (1.0d0)
         else
            mask(i0) = 0
         endif
      enddo
      else
      do i0 = ibxlo0,ibxhi0
         if (test(i0) .eq. val) then
            mask(i0) = 0
         else
            mask(i0) = (1.0d0)
         endif
      enddo
      endif
      return
      end


text:
c retrieve a 1-d Berkeley format BINTABLE spectrum
c and write linearized version to ascii file.

      program writedeepspec

      include 'specarray.h'

      real wave(NWAVEMAX),spec(NWAVEMAX),serr(NWAVEMAX)
      character fname*160,oname*160

      call pgbeg(0,'?',1,1)
      call pgscf(2)
      call pgsch(1.2)

 100  continue
      write(*,'("1-d spec fits file [quit]: ",$)')
      read(*,'(a)') fname
      if (fname(1:5) .eq. '     ') go to 666

 200  continue
      write(*,'("output ascii file [writedeepspec.out]: ",$)')
      read(*,'(a)') oname
      if (oname(1:5) .eq. '     ') then
         oname='writedeepspec.out'
      end if
      open(2,file=oname,status='unknown',err=200)

      call get1ddeepspec(fname,nsp,spec,serr,w0,dw,ifer)
      if (ifer .ne. 0) then
         write(*,*) "Error reading spectrum!"
      end if

      do i=1,nsp
         wave(i) = w0 + i*dw
         write(*,'(i6,3x,f9.3,2(3x,f10.4))') i,wave(i),spec(i),serr(i)
      end do
      close(2)

      call showspecerr(nsp,wave,spec,serr)
      call pglabel("linearized wavelength","counts and error",fname)
      go to 100

 666  continue
      call pgend()

      end


text:
!> Contains wrapper for \c MPI_ABORT.
module mpifx_abort_module
  use mpifx_common_module
  implicit none
  private

  public :: mpifx_abort

contains

  !> Aborts MPI processes for the given communicator.
  !!
  !! \param mycomm  MPI handler.
  !! \param errorcode  Exit error code for the operating system. (default: 1)
  !! \param error  Optional error flag.
  !!
  !! \see MPI documentation (\c MPI_ABORT)
  !!
  !! Example:
  !!
  !!     program test_abort
  !!       use libmpifx_module
  !!       implicit none
  !!
  !!       type(mpifx_comm) :: mycomm
  !!
  !!       call mpifx_init()
  !!       call mycomm%init()
  !!       :
  !!       ! Stoping the program (e.g. due to error we can not handle)
  !!       call mpifx_abort(mycomm, 2)
  !!
  !!     end program test_abort
  !!
  subroutine mpifx_abort(mycomm, errorcode, error)
    type(mpifx_comm), intent(in) :: mycomm
    integer, intent(in), optional :: errorcode
    integer, intent(out), optional :: error

    integer :: error0, errorcode0

    if (present(errorcode)) then
      errorcode0 = errorcode
    else
      errorcode0 = -1
    end if

    call mpi_abort(mycomm%id, errorcode0, error0)
    call handle_errorflag(error0, "MPI_ABORT in mpifx_abort", error)

  end subroutine mpifx_abort

end module mpifx_abort_module

ext: fpp

max_stars_repo_name: terminationshock/mpifx

max_issues_repo_name: terminationshock/mpifx

max_forks_repo_name: terminationshock/mpifx


text:
subroutine plcd_ComputeForcesVectors(a)
   use typre
   use Mod_PLCD
   implicit none
   class(PLCDProblem) :: a

   interface
      subroutine plcd_EnditeElmope(a)
         use typre
         use Mod_PLCD
         implicit none
         class(PLCDProblem) :: a
      end subroutine

      subroutine plcd_Enditebouope(a)
         use typre
         use MOD_PLCD
         implicit none
         class(PLCDProblem) :: a
      end subroutine

      subroutine plcd_bouope(a)
         use typre
         use Mod_PLCD
         implicit none
         class(PLCDProblem) :: a
      end subroutine
   end interface

   integer(ip) :: ndime

   !Set Internal and External Forces to Zero to check convergence
   a%InternalForcesVector = 0.0_rp
   a%ExternalForcesVector = 0.0_rp
   call a%Mesh%GetNdime(ndime)
   a%ExternalForcesVector(1:ndime,:) = a%css%CurrentLoadFactor*a%cs%NodalForces

   !Add Tractions in External Forces Vector
   call plcd_Enditebouope(a)

   !So that we know that we are not inside the endite, but in the begste
   call plcd_EnditeElmope(a)

   call a%Mesh%ArrayCommunicator%GhostCommunicate(size(a%InternalForcesVector,1),a%InternalForcesVector)
   call a%Mesh%ArrayCommunicator%GhostCommunicate(size(a%ExternalForcesVector,1),a%ExternalForcesVector)

   !Computation of the residual vector
   a%ResidualForcesVector = a%ExternalForcesVector-a%InternalForcesVector
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 cliefu()
c
c     This subroutine clears the IEEE flag for floating-point underflow,
c     if such a flag is present, to avoid getting an unnecessary system
c     warning message. Underflow is a normal condition in EQ3/6.
c
c     The presence of an IEEE flag for underflow is platform and
c     compiler dependent. If the flag is present, the need to clear
c     it may also depend on the platform and compiler. Currently,
c     the flag is present and needs to be cleared on Sun SPARCstations
c     in the case of Sun's Fortran 77 compiler. This is not the case
c     with any Fortran 90 compiler that has been seen to date,
c     including Sun's.
c
c     This subroutine is called by:
c
c       EQPT/eqpt.f
c       EQ3NR/eq3nr.f
c       EQ6/eq6.f
c
c-----------------------------------------------------------------------
c
c     Input:
c
c       None
c
c     Output:
c
c       None
c
c-----------------------------------------------------------------------
c
      implicit none
c
c-----------------------------------------------------------------------
c
c     Calling sequence variable declarations.
c
c       None
c
c-----------------------------------------------------------------------
c
com   BEGIN_SPARC_DEPENDENT_CODE
com     BEGIN_F77_DEPENDENT_CODE
com
com       Here cieeef() is a C subroutine that is distributed with
com       EQ3/6. The object file is usually attached to the library
com       (.a) file for the EQLIBU library, for convenience in linking.
com
c         call cieeef()
c         go to 999
com
com     END_F77_DEPENDENT_CODE
com   END_SPARC_DEPENDENT_CODE
c
c* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
  999 continue
      end