	program step2II 

	use karin

c	Byron Gleason, 6/2000
c	modified by Karin Gleason, 07/2001-7/2002
c	modified by Karin Gleason, 06/2015
c	Reads in 5km gridded Tmin file & computes step 1 results.
c	/gcad/nclimdiv/us/por/tmin/pnt/YYYYMM.pnt (input files)

c	The sum of (a) percentage of the US with minimum
c	temperatures much below normal and (b) percentage
c	of the US with minimum temperatures much above normal.

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

c       SYNTAX:
c       dk_step2.II.x iEndYr iEndMo rThresh

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

	parameter ( iBegYear = 1895 )
	parameter ( iEndYear = 2050 )
	parameter ( rMsg = -99.99 )
	parameter ( iMaxGrid= 469758 )

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

	integer iExist(1:iMaxGrid,1:9)
	integer jCount(1:10)
        integer iKountGrids(iBegYear:iEndYear)
        integer iMo(1:5),iMoYr(1:5)
	real rTotalAreaMuchAboveNormal(iBegYear:iEndYear,1:9)
	real rTotalAreaMuchBelowNormal(iBegYear:iEndYear,1:9)
	real rUSArea(iBegYear:iEndYear,1:9)
	real rPercentAreaMuchAboveNormal(iBegYear:iEndYear)
	real rPercentAreaMuchBelowNormal(iBegYear:iEndYear)
        real rRegPercentAreaMuchAboveNormal(iBegYear:iEndYear,1:9)
        real rRegPercentAreaMuchBelowNormal(iBegYear:iEndYear,1:9)
c	real rData(1:iMaxGrid,iBegYear:iEndYear,1:12)
        real, dimension(:,:,:),allocatable :: rData
	real rMetaLat(1:iMaxGrid),rMetaLon(1:iMaxGrid)
	real rLocalData(1:600)
	real rAnnual(iBegYear:iEndYear)
	real rDummyData(1:12)
	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,cJ
        character *4, cI 
        character *5, argv 
        character *8, cPeriod(1:5)
        character *10, cFileName
        character *2, cReg(1:9)
        character *80, cOutFile,cOutFile1,cOutFile2

	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     '/
********************************************************

	allocate(rData(1:iMaxGrid,iBegYear:iEndYear,1:12))

********************************************************
********************************************************
	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)') rHiPrcnt
********************************************************
********************************************************

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
	
        do i=1,iMaxGrid
          do k=iBegYr,iCurYr
            do l=1,12
              rData(i,k,l) = rMsg
            enddo
          enddo
        enddo

************************************************************
c	Initialization
************************************************************

	rLoPrcnt = 100.0 - rHiPrcnt

        do m=1,9
          do i=1,iMaxGrid
            iExist(i,m)=0
          enddo
        enddo

**********************************************************************
c  Read in data and regions ONCE for all standard period calculations:
**********************************************************************
c	Read in all data
	write (*,*) ' '
	write (*,*) 'Reading in POR for 5km gridded Tmin data'
	write (*,*) ' '

	do i=iBegYr,iCurYr
	  do j=1,12

c  Kick out of read loop once current year/month have been stored 
c  into memory:
	    if((i.eq.iCurYr).and.(j.gt.iCurMo)) goto 60

	    write(cI,'(i4)') i
	    write(cJ,'(i2.2)') j
	    cFileName = cI//cJ//'.pnt'
	    write(*,*) 'Reading: ',cFileName

	    open (unit=12,file= cFileName)
	    iCount = 1
5	    read (12,10,end=50) rLat,rLon,rValue
10	    format (f8.4,1x,f9.4,1x,f8.2)

	    rData(iCount,i,j) = rValue
	    rMetaLat(iCount) = rLat 
	    rMetaLon(iCount) = rLon
	    iCount = iCount + 1
	    goto 5
50	    close(12)
 	    enddo
	  enddo	

60	  continue

	write(*,*) 'iCount = ',iCount

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

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

	call regread(jCount,xx,yy)

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

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

	do i=1,iCount
	  x=rMetaLon(i)
	  y=rMetaLat(i)
	    do l=1,10
	      iKount = jCount(l)
	      do m=1,iKount
	        rDummyX(m)=xx(l,m)
	        rDummyY(m)=yy(l,m)
	      enddo

	      iFlag = 0
	      call pnpoly(x,y,rDummyX,rDummyY,iKount,inout)
	      if (inout.ge.0) then
	        iFlag = 1
	        if(l.le.5) then
	          iExist(i,l)=1
c	          print*,'Assigned to ',cReg(l)
	        else
	          iExist(i,l-1)=1
c	          print*,'Assigned to ',cReg(l-1)
	        endif
	        goto 77
	      endif
	    enddo
 77       continue
	enddo

************************************************************
************************************************************
************************************************************
************************************************************
************************************************************
      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	Initialization for each processing period
************************************************************

        do k=iBegYr,iEndYr
          rPercentAreaMuchAboveNormal(k) = 0.0
          rPercentAreaMuchBelowNormal(k) = 0.0
	  iKountGrids(k) = 0
          do m=1,9
            rUSArea(k,m)=0.0
            rTotalAreaMuchAboveNormal(k,m) = 0.0
            rTotalAreaMuchBelowNormal(k,m) = 0.0
            rRegPercentAreaMuchAboveNormal(k,m) = 0.0
            rRegPercentAreaMuchBelowNormal(k,m) = 0.0
          enddo
        enddo

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

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

        open(14,file='extr-tmin-grids.'//cBegMo//'-'//cEndMo//
     +               '.txt')

	do m=1,9

	  do i=1,iCount
c	    Create an array that has average annual tmin for a grid point
	    iCountYears=0	
           do k=iBegYr,iEndYr
             rAnnual(k) = rMsg
           enddo
	   do k=iBegYr,iEndYr
	    if(iExist(i,m).eq.1) then
	      iKountGrids(k) = iKountGrids(k) + 1
	      iCountYears=iCountYears+1
	      iKount=0	
*******************************************************
*       For individual month calculations
*******************************************************
              if(iBegMo.eq.iEndMo) then
                rLocalData(iCountYears) = rData(i,k,iBegMo)
                rAnnual(k)=rData(i,k,iBegMo)
*******************************************************
*       For Winter or other year-spanning calculations
*******************************************************
              elseif(iBegMo.gt.iEndMo) then
                do l=iBegMo,12
                  iKount=iKount+1
                  rDummyData(iKount) = rData(i,k-1,l)
                enddo
                do l=1,iEndMo
                  iKount=iKount+1
                  rDummyData(iKount) = rData(i,k,l)
                enddo
                call average(iKount,rMsg,rDummyData,rAve)
                rLocalData(iCountYears) = rAve
                rAnnual(k)=rAve
*******************************************************
*       For other seasonal/annual or month combinations
*******************************************************
              else
	        do l=iBegMo,iEndMo
	          iKount=iKount+1
	          rDummyData(iKount) = rData(i,k,l)
	        enddo
	        call average(iKount,rMsg,rDummyData,rAve)
	        rLocalData(iCountYears) = rAve
	        rAnnual(k)=rAve
	      endif
	    endif
	   enddo

	    if(iCountYears.eq.0) goto 15
	    
c 	    Sort the rLocalData array in ascending order

	    call sort(iCountYears,rLocalData)

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

            r90thPtileValue=
     +		quickpile(rLocalData,iCountYears,rHiPrcnt,-99.9)


c	    Find 10th percentile value of sorted array.  Anything less
c	    than this value will be considered much below normal. 

            r10thPtileValue=
     +          quickpile(rLocalData,iCountYears,rLoPrcnt,-99.9)

*************************************************************
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.  THEN ...
*************************************************************
c       Go through all years and compare annual values to the 
c	90th and 10th percentile values to determine whether for 
c	that year the particular grid box in question is much
c	above or much below normal 
*************************************************************

	    call area(rMetaLat(i),rArea)
	    do k=iBegYr,iEndYr
	      if(iExist(i,m).eq.1) then
	      rUSArea(k,m)=rUSArea(k,m)+rArea
	      if (rAnnual(k).gt.r90thPtileValue) then
	        rTotalAreaMuchAboveNormal(k,m) = 
     +                   rTotalAreaMuchAboveNormal(k,m) + rArea
                write(14,144) rMetaLon(i),rMetaLat(i),k
144             format(f9.4,1x,f8.4,1x,'90 ',i4)
	      endif
	      if ((rAnnual(k).lt.r10thPtileValue).and.
     +            (rAnnual(k).ne.rMsg)) then
	        rTotalAreaMuchBelowNormal(k,m) = 
     +                   rTotalAreaMuchBelowNormal(k,m) + rArea
                write(14,145) rMetaLon(i),rMetaLat(i),k
145             format(f9.4,1x,f8.4,1x,'10 ',i4)
	      endif
	      endif
	    enddo
 15       continue
	  enddo

       	enddo
	close(14)

***************************************************************
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 minimum temperatures much below normal and above
c	normal.
***************************************************************

c  First, for CONUS:
	do m=1,9
	 do i=iBegYr,iEndYr
  	 rPercentAreaMuchAboveNormal(i) = 
     + 	 (rPercentAreaMuchAboveNormal(i) +
     + (rTotalAreaMuchAboveNormal(i,m)/rUSArea(i,m))*rRegArea(m))
	 enddo
	enddo
	do m=1,9
	 do i=iBegYr,iEndYr
	  rPercentAreaMuchBelowNormal(i) = 
     +	  (rPercentAreaMuchBelowNormal(i) + 
     + (rTotalAreaMuchBelowNormal(i,m)/rUSArea(i,m))*rRegArea(m))
	 enddo
	enddo

c  Now, sum areas for the 9 Regions separately:
        do m=1,9
         do i=iBegYr,iEndYr
         rRegPercentAreaMuchAboveNormal(i,m) =
     +   rRegPercentAreaMuchAboveNormal(i,m) +
     + (rTotalAreaMuchAboveNormal(i,m)/rUSArea(i,m))
         enddo
        enddo
        do m=1,9
         do i=iBegYr,iEndYr
          rRegPercentAreaMuchBelowNormal(i,m) =
     +    rRegPercentAreaMuchBelowNormal(i,m) +
     + (rTotalAreaMuchBelowNormal(i,m)/rUSArea(i,m))
         enddo
        enddo

*************************************************************** 
c	Add much above and below percentages to get the final
c	percentage.
*************************************************************** 

        cOutFile='../data/dk-step2.'//cBegMo//'-'//cEndMo//'.results'
        cOutFile1='../data/dk-step2-lo.'//cBegMo//'-'//cEndMo//
     +            '.results'
        cOutFile2='../data/dk-step2-hi.'//cBegMo//'-'//cEndMo//
     +            '.results'

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

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

c  Write out YTD output:
        if(iLoop.eq.5) then
         cOutFile='../data/dk-step2.ytd.results'
         cOutFile1='../data/dk-step2-lo.ytd.results'
         cOutFile2='../data/dk-step2-hi.ytd.results'
         open (unit=13,file=cOutFile)
         open (unit=14,file=cOutFile1)
         open (unit=15,file=cOutFile2)

         do i=iBegYr,iEndYr
          rFinal=rPercentAreaMuchBelowNormal(i)+
     +           rPercentAreaMuchAboveNormal(i)
          write (13,13) i,rFinal
          write (14,13) i,rPercentAreaMuchBelowNormal(i)
          write (15,13) i,rPercentAreaMuchAboveNormal(i)
         enddo
         close(13)              
         close(14)              
         close(15)              
        endif

c  Now, for Regions:

      do m=1,9
        cOutFile='../data/dk-step2'//
     +              cReg(m)//'.'//cBegMo//'-'//cEndMo//'.results'
        cOutFile1='../data/dk-step2'//
     +              cReg(m)//'-lo.'//cBegMo//'-'//cEndMo//'.results'
        cOutFile2='../data/dk-step2'//
     +              cReg(m)//'-hi.'//cBegMo//'-'//cEndMo//'.results'

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

        do i=iBegYr,iEndYr
          rFinal=rRegPercentAreaMuchBelowNormal(i,m)+
     +           rRegPercentAreaMuchAboveNormal(i,m)
          write (13,13) i,rFinal
          write (14,13) i,rRegPercentAreaMuchBelowNormal(i,m)
          write (15,13) i,rRegPercentAreaMuchAboveNormal(i,m)
        enddo	
        close(13)
        close(14)
        close(15)
      enddo

c  Now, for YTD regional output:
        if(iLoop.eq.5) then
         do m=1,9
          cOutFile = '../data/dk-step2'//cReg(m)//'.ytd.results'
          cOutFile1 = '../data/dk-step2'//cReg(m)//'-lo.ytd.results'
          cOutFile2 = '../data/dk-step2'//cReg(m)//'-hi.ytd.results'    

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

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

      enddo

	write (*,*) 'Done!'
	end

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

c  NOTE:  This is different since each grid box is 5km rather than 1 degree.
c         1 degree of latitude is equivalent to 111.2 km.
	subroutine area(rLat,rArea)
	
	PI = 3.14159
	EARTHRADIUS = 6371.229

	y = rLat * (PI/180.0)

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

	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)

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 Shells 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 Shells 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
