	program step4II 

	use karin

c	This version of step 4 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-10/2002

c	Twice the value of the percentage of the US with a
c	much greater than normal proportion of precipitation
c	derived from extreme 1-day precipitation events.

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_step4.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 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 rTallyExtr(66:126,23:51)
	real rCountYears(66:126,23:51)
	real rRatioYears(66:126,23:51)
	real rCountAll(66:126,23:51,iBegYear:iEndYear)
	real rRatio(66:126,23:51,iBegYear:iEndYear)
	real rTotal(iBegYear:iEndYear)
	real rExtremeTotal(iBegYear:iEndYear)
	real rProportion(iBegYear:iEndYear)
	real rAllProportions(66:126,23:51,iBegYear:iEndYear)
	real rExtrProp(66:126,23:51,iBegYear:iEndYear)
	real rLocalData(1:60000)
        real xx(1:10,1:4000), yy(1:10,1:4000)
        real rDummyX(1:4000), rDummyY(1:4000)
	real rTotalArea(iBegYear:iEndYear,1:9)
	real rPercentArea(iBegYear:iEndYear)
	real rRegPercentArea(iBegYear:iEndYear,1:9)
	real rUSArea(iBegYear:iEndYear,1:9)
        real rRegArea(1:9)

        character *2, cBegMo,cEndMo
	character *4, cYear
	character *5, argv
	character *2, cReg(1:9)
        character *8, cPeriod(1:5)
	character *60, cInfile,cInFile2
        character *80, cOutFile,cExtrGrids,cAllGrids

	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

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
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 (*,*) ' '
	do iYear=iBegYr,iCurYr
	  write (*,*) 'Reading in data for year: ',iYear
	  write (cYear,'(i4.4)') 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.2,1x,i2.2,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)

************************************************************
      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(*,*) iBegYr
        write(*,*) iEndYr
        write(*,*) iBegMo
        write(*,*) iEndMo

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

c Initialize for each loop:
	do i=66,126
	  do j=23,51
	    rTallyExtr(i,j) = rMsg
	    do k=iBegYr,iEndYr
	      rAllProportions(i,j,k) = rMsg
	      rExtrProp(i,j,k) = rMsg
	    enddo
	  enddo
	enddo

	do k=iBegYr,iEndYr
	  rPercentArea(k) = 0.0
	  do n=1,9
	    rUSArea(k,n) = 0.0
	    rTotalArea(k,n) = 0.0
	    rRegPercentArea(k,n) = 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.
              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)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c	Want 90% completeness throughout each year and POR
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	    if(rRatioYears(i,j).ge.0.90) then
              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
************************************************************

	cExtrGrids = '../data/extr-prcp-grids.'//cBegMo//'-'//
     +               cEndMo//'.txt'
	cAllGrids = '../data/all-prcp-grids.'//cBegMo//'-'//
     +               cEndMo//'.txt'
	iiBad = 0
	open(unit=19,file=cExtrGrids)
	open(unit=20,file=cAllGrids)
	do n=1,9

         do iYr=iBegYr,iEndYr
          rPercentArea(iYr)=0.0
          rPercentArea(iYr)=0.0
         enddo

	do i=66,126
	  do j=23,51
c	    Create an array with all daily values for all years
	    iCount=0	
	    do k=iBegYr,iEndYr
	     if (iExist(i,j,k,n).eq.1) then
	     rTallyExtr(i,j) = 0.
*******************************************************
*       For individual month calculations
*******************************************************
              if(iBegMo.eq.iEndMo) then
                do m=1,daysinmth(iBegMo)
	          if(rData(i,j,k,l,m).gt.0.) then
                    iCount=iCount+1
                    rLocalData(iCount) = rData(i,j,k,iBegMo,m)
		  endif
                enddo
*******************************************************
*       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,l,m).gt.0.) then
                    iCount=iCount+1
                    rLocalData(iCount) = rData(i,j,k-1,l,m)
		  endif
                 enddo
                enddo
                do l=1,iEndMo
                 do m=1,daysinmth(l)
	          if(rData(i,j,k,l,m).gt.0.) then
                    iCount=iCount+1
                    rLocalData(iCount) = rData(i,j,k,l,m)
	          endif
                 enddo
                enddo
*******************************************************
*       For other seasonal/annual or month combinations
*******************************************************
              else
	       do l=iBegMo,iEndMo
	        do m=1,daysinmth(l)
	          if(rData(i,j,k,l,m).gt.0.) then
	            iCount=iCount+1
	            rLocalData(iCount) = rData(i,j,k,l,m)
		  endif
	        enddo
	       enddo
	      endif
	      endif
	     enddo

	    if(iCount.eq.0) goto 12

c 	    Sort the rLocalData array in ascending order

	    call sort(iCount,rLocalData)

c	    Find 90th percentile value of sorted array.  Anything greater
c	    than this value will be considered 'extreme' precipitation. 

	    r90thPtileValue=quickpile(rLocalData,iCount,rPercentile,rMsg)

	    if(r90thPtileValue.eq.0.0) then
	      iiBad = iiBad + 1
	      do k=iBegYr,iEndYr
	        iExist(i,j,k,n) = 0
	      enddo
	      goto 12 
            else
              iFixI = i*(-1)
              iFixJ = j+1
	    endif

c	    Go through all days for a given year and calculate total   
c	    precipitation, and extreme total precipitation.  Then
c	    calculate the proportion of extreme to total.

	    do k=iBegYr,iEndYr
	    rProportion(k)=rMsg
	    if(iExist(i,j,k,n).eq.1) then
	    rExtrProp(i,j,k) = 1.5
	    rTotal(k)=0.0
	    rExtremeTotal(k)=0.0
*******************************************************
*       For individual month calculations
*******************************************************
            if(iBegMo.eq.iEndMo) then
              do m=1,daysinmth(iBegMo)
              if(rData(i,j,k,iBegMo,m).ge.0.0) then
                rTotal(k)=rTotal(k)+rData(i,j,k,iBegMo,m)
                if (rData(i,j,k,iBegMo,m).gt.r90thPtileValue) then
                  rExtremeTotal(k)=
     +                    rExtremeTotal(k)+rData(i,j,k,iBegMo,m)
                endif
              endif
              enddo
*******************************************************
*       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).ge.0.0) then
                rTotal(k)=rTotal(k)+rData(i,j,k-1,l,m)
                if (rData(i,j,k-1,l,m).gt.r90thPtileValue) then
                  rExtremeTotal(k)=rExtremeTotal(k)+rData(i,j,k-1,l,m)
                endif
              endif
              enddo
              enddo
              do l=1,iEndMo
              do m=1,daysinmth(l)
              if(rData(i,j,k,l,m).ge.0.0) then
                rTotal(k)=rTotal(k)+rData(i,j,k,l,m)
                if (rData(i,j,k,l,m).gt.r90thPtileValue) then
                  rExtremeTotal(k)=rExtremeTotal(k)+rData(i,j,k,l,m)
                endif
              endif
              enddo
              enddo
*******************************************************
*       For other seasonal/annual or month combinations
*******************************************************
              else
	       do l=iBegMo,iEndMo
	        do m=1,daysinmth(l)
                if(rData(i,j,k,l,m).ge.0.0) then
	          rTotal(k)=rTotal(k)+rData(i,j,k,l,m)
	          if (rData(i,j,k,l,m).gt.r90thPtileValue) then
                    rExtremeTotal(k)=rExtremeTotal(k)+rData(i,j,k,l,m)
	   	  endif
	   	endif
	        enddo
	       enddo
	      endif

	    if (rTotal(k).eq.0.0) then
	      rProportion(k)=0.0
	    else
	      rProportion(k)=rExtremeTotal(k)/rTotal(k)
	    endif
	    endif
	    enddo
                     
c	    Sort the values in the rProportion array in ascending order

	    iCount=0	  
	    do k=iBegYr,iEndYr
	      if((rProportion(k).ne.0.).and.
     +           (rProportion(k).ne.rMsg)) then
	        iCount=iCount+1
	        rLocalData(iCount)=rProportion(k)
	 	rAllProportions(i,j,k) = rProportion(k)
	      endif
	    enddo

	    call sort(iCount,rLocalData)

c	    find 90th percentile value in rLocalData array

	    r90thPtileValue=quickpile(rLocalData,iCount,rPercentile,rMsg)

cccccccccccccccccccccccc
c	Write to extreme-prop.txt file 
cccccccccccccccccccccccc
      	    iFixI = i*(-1)

c	    Loop through each year to determine if proportion of extreme
c	    precip. to total precip. is greater than the 90th percentile
c	    value for the period of record proportions.

********************************************************
	    do k=iBegYr,iEndYr
	      if (rProportion(k).gt.r90thPtileValue) then
c  Write out only extreme grids to output file:
              write(19,202) iFixI,j,k
202           format(i4,1x,i2,1x,'90 ',i4)
********************************************************
	        call area(j,rArea)
	        rTotalArea(k,n)=rTotalArea(k,n)+rArea
		rExtrProp(i,j,k) = rProportion(k)
		rTallyExtr(i,j) = rTallyExtr(i,j) + 1.
	      endif
	    enddo

 12       continue
	  enddo
	enddo

	enddo

***************************************************************
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
      	    iFixI = i*(-1)
	   do j=23,51
	    if (iExist(i,j,k,n).eq.1) then
              write(20,203) iFixI,j,k
203           format(i4,1x,i2,1x,'150 ',i4)
	      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 by year the percentage of the contiguous U.S. 
c	with much greater than normal proportion of precipitation
c	derived from extreme 1-day precipitation events.
***************************************************************

	do n=1,9
	 do i=iBegYr,iEndYr
c if statement needed for NRT data calculations where grid points
c do not exist in a particular region:
	  if(rUSArea(i,n).ne.0.0) then
	    rPercentArea(i) = (rPercentArea(i) + 
     +         (rTotalArea(i,n)/rUSArea(i,n))*rRegArea(n))
c  Do also for 9 regions separately:
	    rRegPercentArea(i,n) = (rRegPercentArea(i,n) + 
     +         (rTotalArea(i,n)/rUSArea(i,n)))
	  endif
	 enddo
	enddo

	close(19)
	close(20)

***************************************************************
***************************************************************
c	Write out final results 
***************************************************************
***************************************************************

        cOutFile='../data/dk-step4.'//cBegMo//'-'//cEndMo//'.results'

	open (unit=13,file=cOutFile) 

	do i=iBegYr,iEndYr
	  write (13,13) i,rPercentArea(i)
13	  format(i4,1x,f5.3)
	enddo
	  
	close(13)	  	

c  Now, write output to ytd file:
	
	if(iLoop.eq.5) then
         cOutFile='../data/dk-step4.ytd.results'

	 open (unit=13,file=cOutFile) 

	  do i=iBegYr,iEndYr
	    write (13,13) i,rPercentArea(i)
	  enddo
	  
	  close(13)	  	

	 endif

c  Now, for Regions:

      do m=1,9
        cOutFile = '../data/dk-step4'//
     +              cReg(m)//'.'//cBegMo//'-'//cEndMo//'.results'
        open (unit=13,file=cOutFile)

        do i=iBegYr,iEndYr
          write (13,13) i,rRegPercentArea(i,m)
        enddo
        close(13)

	if(iLoop.eq.5) then
          cOutFile = '../data/dk-step4'//cReg(m)//'.ytd.results'
          open (unit=13,file=cOutFile)

           do i=iBegYr,iEndYr
            write (13,13) i,rRegPercentArea(i,m)
           enddo
          close(13)

	endif
      enddo

      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*80, 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 = '/monitoring-shared/workspace/resources/'//
     +      'polygons/usa/9_climate_regions/'//cRegion(i)//'.poly'
          elseif((i.lt.4).or.(i.eq.7).or.(i.eq.8)) then
            cInFile='/monitoring-shared/workspace/resources/'//
     +      'polygons/usa/9_climate_regions/'//cRegion2(i)//'.poly'
          else
            cInFile='/monitoring-shared/workspace/resources/'//
     +      'polygons/usa/9_climate_regions/'//cRegion3(i)//'.poly'
          endif
          open(unit=14,file=cInFile)

c          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
