	program step5II 

	use karin

c       This version of step 5 looks at the gridded files which are
c       made up of one station per grid box.  The station chosen is
c       the one which has the largest POR of precip within that grid.

c	Byron Gleason, 6/2000
c	Karin Gleason, 7/2001-7/2002

c	The sum of (a) percentage of the US with a much greater
c	than normal number of days with precipitation and (b)
c	percentage of the US with a much greater than normal
c	number of days without precipitation.

c	User must select:
c	1)   beginning and ending year
c	2)   beginning and ending month
c	3)   threshold percentile for these extreme events
c	     (usually 90.0)

c	SYNTAX:  
c	dk_step5.II.x iBegYr iEndYr iBegMo iEndMo rThresh
***********************************************************

	parameter ( iBegYear = 1900 )
	parameter ( iEndYear = 2050 )
	parameter ( rMsg = -99.9 )

***********************************************************

	integer daysinmth(1:12)
	integer iExist(66:126,23:51,iBegYear:iEndYear,1:9)
        integer jCount(1:10),iMo(1:5), iMoYr(1:5)
	real rRegPrcntAreaMuchAbvNormalWith(iBegYear:iEndYear,1:9)
	real rRegPrcntAreaMuchAbvNormalWithOut(iBegYear:iEndYear,1:9)
	real rPrcntAreaMuchAbvNormalWith(iBegYear:iEndYear)
	real rPrcntAreaMuchAbvNormalWithOut(iBegYear:iEndYear)
	real rTotAreaMuchAbvNormalWith(iBegYear:iEndYear,1:9)
	real rTotAreaMuchAbvNormalWithOut(iBegYear:iEndYear,1:9)
	real rData(66:126,23:51,iBegYear:iEndYear,1:12,1:31)
        real rCount(66:126,23:51,iBegYear:iEndYear)
        real rCountAllYears(66:126,23:51)
        real rCountYears(66:126,23:51)
        real rRatioYears(66:126,23:51)
        real rRatio(66:126,23:51,iBegYear:iEndYear)
        real rCountAll(66:126,23:51,iBegYear:iEndYear)
	real rLocalWith(1:150)
	real rLocalWithOut(1:150)
	real rAnnualWith(iBegYear:iEndYear)
	real rAnnualWithOut(iBegYear:iEndYear)
	real rUSArea(iBegYear:iEndYear,1:9)
        real xx(1:10,1:4000), yy(1:10,1:4000)
        real rDummyX(1:4000), rDummyY(1:4000)
        real rRegArea(1:9)

        character *2, cBegMo,cEndMo
	character *4, cYear
	character *5, argv 
	character *8, cPeriod(1:5)
	character *2, cReg(1:9)
	character *60, cInfile
        character *80, cOutFile,cOutFile1,cOutFile2,cExtrFile

	data daysinmth/31,28,31,30,31,30,31,31,30,31,30,31/
        data rRegArea/0.08230,0.08908,0.15551,0.14053,
     +                0.08428,0.18822,0.10271,0.06021,
     +                0.09715/

        data cReg/'nw','we','wn','sw','en','so','ce',
     +            'ne','se'/

        data cPeriod/'Current ','3-Month ','6-Month ',
     +               '12-Month','YTD     '/
********************************************************

        call GETARG (1,argv)
        read(argv,'(i4.4)') iBegYr
        call GETARG (2,argv)
        read(argv,'(i4.4)') iCurYr
	call GETARG (3,argv)
        read(argv,'(i2.2)') iCurMo
        call GETARG (4,argv)
        read(argv,'(f5.1)') rPercentile

********************************************************
********************************************************

c  Determine all month and year variables for production run:
c       current month, 3-mo, 6-mo, 12-mo and YTD
        iMo(1) = iCurMo
        iMoYr(1) = iBegYr

        if(iCurMo.ge.3) then
          iMo(2) = iCurMo-2
          iMoYr(2) = iBegYr
        else
          iMo(2) = iCurMo + 10
          iMoYr(2) = iBegYr + 1
        endif
        if(iCurMo.ge.6) then
          iMo(3) = iCurMo-5
          iMoYr(3) = iBegYr
        else
          iMo(3) = iCurMo + 7
          iMoYr(3) = iBegYr + 1
        endif
        if(icurMo.eq.12) then
          iMo(4) = 1
          iMoYr(4) = iBegYr
        else
          iMo(4) = iCurMo + 1
          iMoYr(4) = iBegYr + 1
        endif

c	Initialization of Data Read

	do i=66,126
	  do j=23,51
	    do k=iBegYr,iCurYr
	      do l=1,12
	        do m=1,31
	          rData(i,j,k,l,m) = rMsg
	        enddo
	      enddo
	    enddo
	  enddo
	enddo

************************************************************

c	Read in all data

	write (*,*) ' '
	write (*,*) 'Reading in Data'
	write (*,*) ' '

	do iYear=iBegYr,iCurYr
	  write (cYear,'(i4.4)') iYear
	  print*,'Reading in :',iYear
	  cInfile=cYear//'.1stn.ghcn.dat'
	  open (unit=12,file=cInfile)
5	  read (12,10,end=50) iMonth,iDay,rLon,rLat,rValue
10	  format (i2,1x,i2,1x,f8.2,1x,f8.2,1x,f8.2)

	  igp=int(abs(rLon))
	  jgp=int(abs(rLat)) 
	  if(rValue.eq.-888.) then
	    rData(igp,jgp,iYear,iMonth,iDay) = 0.005
	  else
	    rData(igp,jgp,iYear,iMonth,iDay) = rValue
	  endif
	 
	  goto 5

50	  close(12)
	enddo

************************************************************
c       Read in 9 region polygon files
************************************************************

        print*,'Reading in 9 regional polygon files'
        print*,' '

        call regread(jCount,xx,yy)

************************************************************
***Compute for all periods**********************************
************************************************************
************************************************************
************************************************************
      Do iLoop=1,5

        write(*,*) 'Working on: ', cPeriod(iLoop)

        if(iLoop.ne.5) then
          iBegMo = iMo(iLoop)
          iBegYr = iMoYr(iLoop)
        else
          iBegMo = 1
          iBegYr = iMoYr(1)
        endif
          iEndMo = iCurMo
          iEndYr = iCurYr

        write(cBegMo,'(i2.2)') iBegMo
        write(cEndMo,'(i2.2)') iEndMo

c  Initialize for each loop:

        do k=iBegYr,iEndYr
	  rPrcntAreaMuchAbvNormalWith(k)=0.0
	  rPrcntAreaMuchAbvNormalWithOut(k)=0.0
          do n=1,9
            rUSArea(k,n) = 0.0
	    rTotAreaMuchAbvNormalWith(k,n)=0.0
	    rTotAreaMuchAbvNormalWithOut(k,n)=0.0
	    rRegPrcntAreaMuchAbvNormalWith(k,n)=0.0
	    rRegPrcntAreaMuchAbvNormalWithOut(k,n)=0.0
	    do i=66,126
	      do j=23,51
                iExist(i,j,k,n)=0
              enddo
            enddo
          enddo
        enddo

************************************************************
c       Determine completeness of data for each year and POR
************************************************************
        write (*,*) ' '
        write (*,*) 'Determining valid data... '
        write (*,*) ' '

        do i=66,126
          do j=23,51
            do k=iBegYr,iEndYr
              rCountAll(i,j,k) = 0.
              rCount(i,j,k) = 0.
              if(iBegMo.gt.iEndMo) then
                do l=iBegMo,12
                  do m=1,daysinmth(l)
                    rCountAll(i,j,k) = rCountAll(i,j,k) + 1.
                    if(rData(i,j,k-1,l,m).ne.rMsg) then
                      rCount(i,j,k) = rCount(i,j,k) + 1.
                    endif
                  enddo
                enddo
                do l=1,iEndMo
                  do m=1,daysinmth(l)
                    rCountAll(i,j,k) = rCountAll(i,j,k) + 1.
                    if(rData(i,j,k,l,m).ne.rMsg) then
                      rCount(i,j,k) = rCount(i,j,k) + 1.
                    endif
                  enddo
                enddo
                if(rCountAll(i,j,k).ne.0.0) then
                  rRatio(i,j,k) = rCount(i,j,k)/rCountAll(i,j,k)
                else
                  rRatio(i,j,k) = 0.0
                endif
              else
                do l=iBegMo,iEndMo
                  do m=1,daysinmth(l)
                    rCountAll(i,j,k) = rCountAll(i,j,k) + 1.
                    if(rData(i,j,k,l,m).ne.rMsg) then
                      rCount(i,j,k) = rCount(i,j,k) + 1.
                    endif
                  enddo
                enddo
                if(rCountAll(i,j,k).ne.0.0) then
                  rRatio(i,j,k) = rCount(i,j,k)/rCountAll(i,j,k)
                else
                  rRatio(i,j,k) = 0.0
                endif
              endif
            enddo
            rCountAllYears(i,j) = 0.
            rCountYears(i,j) = 0.
            do k=iBegYr,iEndYr
              rCountAllYears(i,j) = rCountAllYears(i,j) + 1.
              if(rRatio(i,j,k).ge.0.90) then
               rCountYears(i,j) = rCountYears(i,j) + 1.
c  If last year is msg due to NRT report not available,
c  do not include: 4/2005 klg
              elseif(k.eq.iEndYr) then
               rCountAllYears(i,j) = rCountAllYears(i,j) - 1.
              endif
            enddo
            if(rCountAllYears(i,j).ne.0.0) then
              rRatioYears(i,j) = rCountYears(i,j)/rCountAllYears(i,j)
            else
              rRatioYears(i,j) = 0.0
            endif
          enddo
        enddo

**************************************************************
c       Find out which grid points are in each of the 9
c       regions.
**************************************************************

        print*,'Assigning grid points to each region'
        print*,' '

        do i=126,66, -1
          do j=51,23, -1
	    iFlag = 0
            x=real(i)*(-1.0)
            y=real(j)
            if(rRatioYears(i,j).ge.0.90) then
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c       Want 90% completeness throughout each year and POR
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
              do l=1,10
                iKount = jCount(l)
                do m=1,iKount
                  rDummyX(m)=xx(l,m)
                  rDummyY(m)=yy(l,m)
                enddo

                call pnpoly(x,y,rDummyX,rDummyY,iKount,inout)
                if (inout.ge.0) then
	 	  iFlag = 1
                  do k=iBegYr,iEndYr
                   if(rRatio(i,j,k).ge.0.90) then
                    if(l.le.5) then
                      iExist(i,j,k,l)=1
c                      print*,'Assigned to ',cReg(l)
                    else
                      iExist(i,j,k,l-1)=1
c                      print*,'Assigned to ',cReg(l-1)
                    endif
                   endif
                  enddo
                  goto 21
                endif
              enddo
*********** For grid points which are just outside the regional **********
*********** boundaries, assign to appropriate regions:          **********
              if(iFlag.eq.0) then
                if((i.ge.116).and.(j.ge.42)) then
                  iNewGrid=1
                elseif((i.ge.115).and.(j.lt.42)) then
                  iNewGrid=2
                elseif((i.ge.95).and.(j.ge.40)) then
                  iNewGrid=3
                elseif((i.ge.107).and.(j.lt.40)) then
                  iNewGrid=4
                elseif((i.ge.81).and.(j.ge.40)) then
                  iNewGrid=5
                elseif((i.ge.88).and.(j.lt.40)) then
                  iNewGrid=6
                elseif((i.lt.81).and.(j.ge.38)) then
                  iNewGrid=8
                elseif((i.lt.88).and.(j.lt.38)) then
                  iNewGrid=9
                endif
                do k=iBegYr,iEndYr
                  if(rRatio(i,j,k).ge.0.90) then
                    iExist(i,j,k,iNewGrid)=1
                  endif
                enddo
              endif
            endif
21          continue
          enddo
        enddo

**************************************************************
c  	Analysis
**************************************************************

	cExtrFile ='extr-days-grids.'//cBegMo//
     +              '-'//cEndMo//'.txt'
	open(19,file=cExtrFile)

	write (*,*) 'Calculating ...'
	write (*,*) ' '

	do n=1,9

	 do iYr=iBegYr,iEndYr
	  rPrcntAreaMuchAbvNormalWith(iYr)=0.0
	  rPrcntAreaMuchAbvNormalWithOut(iYr)=0.0
	 enddo

	do i=66,126
	  do j=23,51
	   iCountYears=0	

	    do k=iBegYr,iEndYr
	     rAnnualWith(k) = rMsg
	     rAnnualWithout(k) = rMsg
	     if(iExist(i,j,k,n).eq.1) then
	      iCountYears=iCountYears+1
	      rDaysWithPrecip=0.0
	      rDaysWithOutPrecip=0.0 
*******************************************************
*       For individual month calculations
*******************************************************
              if(iBegMo.eq.iEndMo) then
                do m=1,daysinmth(iBegMo)
                  if (rData(i,j,k,iBegMo,m).eq.0.0) then
                    rDaysWithOutPrecip=rDaysWithOutPrecip+1.0
                  elseif ((rData(i,j,k,iBegMo,m).gt.0.0).or.
     +                    (rData(i,j,k,iBegMo,m).eq.-888.)) then
                    rDaysWithPrecip=rDaysWithPrecip+1.0
                  endif
                enddo
                rAnnualWith(k) = rDaysWithPrecip
                rAnnualWithOut(k) = rDaysWithOutPrecip
                rLocalWith(iCountYears) = rDaysWithPrecip
                rLocalWithOut(iCountYears) = rDaysWithOutPrecip
*******************************************************
*       For Winter or other year-spanning calculations
*******************************************************
              elseif(iBegMo.gt.iEndMo) then
              do l=iBegMo,12
                do m=1,daysinmth(l)
                  if (rData(i,j,k-1,l,m).eq.0.0) then
                    rDaysWithOutPrecip=rDaysWithOutPrecip+1.0
                  elseif (rData(i,j,k-1,l,m).ne.rMsg) then
                    rDaysWithPrecip=rDaysWithPrecip+1.0
                  endif
                enddo
              enddo
              do l=1,iEndMo
                do m=1,daysinmth(l)
                  if (rData(i,j,k,l,m).eq.0.0) then
                    rDaysWithOutPrecip=rDaysWithOutPrecip+1.0
                  elseif (rData(i,j,k,l,m).ne.rMsg) then
                    rDaysWithPrecip=rDaysWithPrecip+1.0
                  endif
                enddo
              enddo
              rAnnualWith(k) = rDaysWithPrecip
              rAnnualWithOut(k) = rDaysWithOutPrecip
              rLocalWith(iCountYears) = rDaysWithPrecip
              rLocalWithOut(iCountYears) = rDaysWithOutPrecip
*******************************************************
*       For other seasonal/annual or month combinations
*******************************************************
              else
	       do l=iBegMo,iEndMo
	        do m=1,daysinmth(l)
	          if (rData(i,j,k,l,m).eq.0.0) then
	            rDaysWithOutPrecip=rDaysWithOutPrecip+1.0
                  elseif(rData(i,j,k,l,m).ne.rMsg) then
	            rDaysWithPrecip=rDaysWithPrecip+1.0
		  endif
	        enddo
	       enddo
               rAnnualWith(k) = rDaysWithPrecip
               rAnnualWithOut(k) = rDaysWithOutPrecip
               rLocalWith(iCountYears) = rDaysWithPrecip
               rLocalWithOut(iCountYears) = rDaysWithOutPrecip

	      endif

	      if((rAnnualWith(k).eq.rMsg).and.
     +          (rAnnualWithOut(k).eq.rMsg)) then
	         iCountYears = iCountYears - 1
	 	 goto 12
	      endif
	     endif
 12 	    continue
	    enddo

	    if(iCountYears.eq.0) goto 129

c 	    Sort the rLocal array in ascending order

	    call sort(iCountYears,rLocalWith)
	    call sort(iCountYears,rLocalWithOut)

c	    Find 90th percentile value of sorted array.  Anything greater
c	    than this value will be considered much above normal. 

            r90thPtileValueWith=
     +		quickpile(rLocalWith,iCountYears,rPercentile,-99.9)
            r90thPtileValueWithOut=
     +		quickpile(rLocalWithOut,iCountYears,rPercentile,-99.9)

c	    Go through all years and compare annual values to the 
c	    90th percentile values to determine whether for 
c	    that year the particular grid box in question is much above 
c	    normal 

	    call area(j,rArea)
	    do k=iBegYr,iEndYr
	      if(iExist(i,j,k,n).eq.1) then
	      if ((rAnnualWith(k).gt.r90thPtileValueWith).and.
     +            (rAnnualWith(k).ne.rMsg)) then
	        rTotAreaMuchAbvNormalWith(k,n) = 
     +                   rTotAreaMuchAbvNormalWith(k,n) + rArea
		write(19,777) -i,j,k
777		format(i5,1x,i2,1x,'90 ',i4)
	      endif
	      if ((rAnnualWithOut(k).gt.r90thPtileValueWithOut).and.
     +            (rAnnualWithOut(k).ne.rMsg)) then
	        rTotAreaMuchAbvNormalWithOut(k,n) = 
     +                   rTotAreaMuchAbvNormalWithOut(k,n) + rArea
		write(19,778) -i,j,k
778		format(i5,1x,i2,1x,'10 ',i4)
	      endif
	
	      endif
	    enddo
 129 	   continue
	  enddo
	enddo

	enddo

	close(19)

*************************************************************
c	Loop through all grid points and sum the areas of those
c	grid points with valid data for each year for each of
c	the 9 regions.
*************************************************************

	do n=1,9
	 do k=iBegYr,iEndYr
	  do i=66,126
	   do j=23,51
	    if (iExist(i,j,k,n).eq.1) then
	      call area(j,rArea)
	      rUSArea(k,n)=rUSArea(k,n)+rArea
	    endif
	   enddo
	  enddo
	 enddo
	enddo

***************************************************************
c       Area weight each region before computing a national
c       value.
c
c 	Calculate Area of Contiguous U.S. with meeting much
c	above and much below normal criteria.
***************************************************************

	do n=1,9
	 do i=iBegYr,iEndYr
	  if(rUSArea(i,n).ne.0.0) then
	    rPrcntAreaMuchAbvNormalWith(i) = 
     +	    (rPrcntAreaMuchAbvNormalWith(i) + 
     +      (rTotAreaMuchAbvNormalWith(i,n)/rUSArea(i,n))
     +      *rRegArea(n))
	    rPrcntAreaMuchAbvNormalWithOut(i) = 
     +	    (rPrcntAreaMuchAbvNormalWithOut(i) + 
     +      (rTotAreaMuchAbvNormalWithOut(i,n)/rUSArea(i,n))
     +      *rRegArea(n))
c  Do also for each region:
	    rRegPrcntAreaMuchAbvNormalWith(i,n) = 
     +	    (rRegPrcntAreaMuchAbvNormalWith(i,n) + 
     +      (rTotAreaMuchAbvNormalWith(i,n)/rUSArea(i,n)))
	    rRegPrcntAreaMuchAbvNormalWithOut(i,n) = 
     +	    (rRegPrcntAreaMuchAbvNormalWithOut(i,n) + 
     +      (rTotAreaMuchAbvNormalWithOut(i,n)/rUSArea(i,n)))
	  endif
	 enddo
	enddo

*************************************************************** 

c	Write out final results 

        cOutFile = '../data/dk-step5.'//
     +             cBegMo//'-'//cEndMo//'.results'
        cOutFile1 = '../data/dk-step5-wi.'//
     +             cBegMo//'-'//cEndMo//'.results'
        cOutFile2 = '../data/dk-step5-wo.'//
     +             cBegMo//'-'//cEndMo//'.results'

	open (unit=13,file=cOutFile) 
	open (unit=14,file=cOutFile1) 
	open (unit=15,file=cOutFile2) 

	do i=iBegYr,iEndYr
	  rFinal=rPrcntAreaMuchAbvNormalWith(i)+
     +           rPrcntAreaMuchAbvNormalWithOut(i)
	  write (13,13) i,rFinal
	  write (14,13) i,rPrcntAreaMuchAbvNormalWith(i)
	  write (15,13) i,rPrcntAreaMuchAbvNormalWithOut(i)
13	  format(i4,1x,f5.3)
	enddo
	  
	close(13)	  	
	close(14)	  	
	close(15)	  	

	if(iLoop.eq.5 ) then
         cOutFile = '../data/dk-step5.ytd.results'
         cOutFile1 = '../data/dk-step5-wi.ytd.results'
         cOutFile2 = '../data/dk-step5-wo.ytd.results'

	 open (unit=13,file=cOutFile) 
	 open (unit=14,file=cOutFile1) 
	 open (unit=15,file=cOutFile2) 

	 do i=iBegYr,iEndYr
	  rFinal=rPrcntAreaMuchAbvNormalWith(i)+
     +           rPrcntAreaMuchAbvNormalWithOut(i)
	  write (13,13) i,rFinal
	  write (14,13) i,rPrcntAreaMuchAbvNormalWith(i)
	  write (15,13) i,rPrcntAreaMuchAbvNormalWithOut(i)
	 enddo
	  
	 close(13)	  	
	 close(14)	  	
	 close(15)	  	

	endif

c  Now, for Regions;

      do m=1,9
        cOutFile = '../data/dk-step5'//
     +              cReg(m)//'.'//cBegMo//'-'//cEndMo//'.results'
        cOutFile1 = '../data/dk-step5'//
     +              cReg(m)//'-wi.'//cBegMo//'-'//cEndMo//'.results'
        cOutFile2 = '../data/dk-step5'//
     +              cReg(m)//'-wo.'//cBegMo//'-'//cEndMo//'.results'

        open (unit=13,file=cOutFile)
        open (unit=14,file=cOutFile1)
        open (unit=15,file=cOutFile2)

        do i=iBegYr,iEndYr
          rFinal=rRegPrcntAreaMuchAbvNormalWith(i,m)+
     +           rRegPrcntAreaMuchAbvNormalWithOut(i,m)
          write (13,13) i,rFinal
          write (14,13) i,rRegPrcntAreaMuchAbvNormalWith(i,m)
          write (15,13) i,rRegPrcntAreaMuchAbvNormalWithOut(i,m)
         enddo
       enddo
	close(13)
	close(14)
	close(15)
	
	if(iLoop.eq.5) then
         do m=1,9
          cOutFile = '../data/dk-step5'//cReg(m)//'.ytd.results'
          cOutFile1 = '../data/dk-step5'//
     +              cReg(m)//'-wi.ytd.results'
          cOutFile2 = '../data/dk-step5'//
     +              cReg(m)//'-wo.ytd.results'

          open (unit=13,file=cOutFile)
          open (unit=14,file=cOutFile1)
          open (unit=15,file=cOutFile2)

          do i=iBegYr,iEndYr
            rFinal=rRegPrcntAreaMuchAbvNormalWith(i,m)+
     +           rRegPrcntAreaMuchAbvNormalWithOut(i,m)
            write (13,13) i,rFinal
            write (14,13) i,rRegPrcntAreaMuchAbvNormalWith(i,m)
            write (15,13) i,rRegPrcntAreaMuchAbvNormalWithOut(i,m)
           enddo
          enddo
	  close(13)
	  close(14)
	  close(15)
	endif

      Enddo
	
	write (*,*) 'Done!'
	write (*,*) ' '

	end

***********************************************************

	subroutine area(j,rArea)
	
	PI = 3.14159
	EARTHRADIUS = 6371.229

	y=real(j)
	y = y * (PI/180.0)

	rArea = ((PI**2)*(EARTHRADIUS**2)*cos(y))/(180**2)

	return
	end

*********************************************************** 

c       this subroutine reads regional polygon files

        subroutine regread(jCount,xx,yy)

        integer jCount(1:10)
        real xx(1:10,1:4000), yy(1:10,1:4000)
        character*2, cRegion(1:4)
        character*3, cRegion2(1:8)
        character*4, cRegion3(1:10)
        character*60, cInFile

        data cRegion /'  ','  ','  ','sw'/
        data cRegion2 /'nw3','we4','wnc','  ','   ','   ',
     +                'so9','mw2'/
        data cRegion3 /'    ','    ','    ','    ','enc5',
     +                'enc6','    ','    ','ne10','se14'/

        do i=1,10
          if(i.eq.4) then
            cInFile = cRegion(i)//'.poly'
          elseif((i.lt.4).or.(i.eq.7).or.(i.eq.8)) then
            cInFile=cRegion2(i)//'.poly'
          else
            cInFile= cRegion3(i)//'.poly'
          endif
          open(unit=14,file=cInFile)

          print*,i,cInFile

          jCount(i) = 0
88        read(14,*,end=99) rLon,rLat
          jCount(i) = jCount(i) + 1
          xx(i,jCount(i)) = rLon
          yy(i,jCount(i)) = rLat
          goto 88
99        close(14)
        enddo

        return
        end

*********************************************************** 

c       this subroutine sorts the data

	subroutine sort(n,arr)

c	this is two different sort routines depending on length
c	the straight sort (part 1) is  only for n<50, as n
c	gets large, the time gets too long (goes with n**2)
c	the second is Shell's method.
c	both came from Numerical Recipes.  There was a warning
c	in there not to use a bubble sort


	real arr(n)


	if(n.gt.50)go to 500

	do j=2,n
	 a=arr(j)
	 do i=j-1,1,-1
	  if(arr(i).le.a)go to 100
	  arr(i+1)=arr(i)
	 enddo
	 i=0
 100	 continue
	 arr(i+1)=a
	enddo

	go to 900
 500    continue
c     this is Shell's routine for n>50

	aln2i=1./0.69314718
	tiny=1.E-5
	lognb2=int(alog(float(n))*aln2i+tiny)
	m=n
	do nn=1,lognb2
	 m=m/2
	 k=n-m
	 do j=1,k
	  i=j
 553	  continue
	  l=i+m
	  if(arr(l).lt.arr(i))then
	   t=arr(i)
	   arr(i)=arr(l)
	   arr(l)=t
	   i=i-m
	   if(i.ge.1)go to 553
	  endif
	 enddo
	enddo

 900    continue
	return

	end
