	program step3II 

	use karin

c	Karin Gleason, 7/2002

c	The sum of (a) percentage of the US in severe
c	drought based on the PDSI and (b) percentage of the
c	US with severe moisture surplus based on the PDSI.

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_step3.II.x iBegYr iEndYr iBegMo iEndMo rThresh

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

	parameter ( iBegYear = 1900 )
	parameter ( iEndYear = 2050 )

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

	integer daysinmth(1:12),iDivisions(1:48)
	integer iExist(1:48,1:10)
	integer iMo(1:5),iMoYr(1:5)
	real rPercentAreaMoist(iBegYear:iEndYear)
	real rPercentAreaDrought(iBegYear:iEndYear)
	real rTotalAreaMoist(iBegYear:iEndYear)
	real rTotalAreaDrought(iBegYear:iEndYear)
	real rTotalArea(iBegYear:iEndYear)
	real rData(1:48,1:10,iBegYear:iEndYear,1:12)
	real rLocalData(1:600)
	real rAnnual(iBegYear:iEndYear)
	real rDummyData(1:12)
	real rDumArea(1:48,1:11)
	real rVal(1:13)
	real rArea(1:11)
        character*2, cBegMo,cEndMo
        character*4, argv
        character*62, cFile
        character *2, cReg(1:9)
        character *8, cPeriod(1:5)
        character*80, cOutFile,cOutFile1,cOutFile2
        character*80, cOutFile3,cOutFile31,cOutFile32

	data daysinmth/31,28,31,30,31,30,31,31,30,31,30,31/
	data iDivisions/8,7,9,7,5,3,2,7,9,10,9,9,9,9,4,9,3,
     +                  8,3,10,9,10,6,7,9,4,2,3,8,10,8,9,
     +                 10,9,9,10,1,7,9,4,10,7,3,6,10,6,9,10/
       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,'(f4.1)') rHiPrcnt

	rLoPrcnt = 100.0 - 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

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

	do i=1,48
	  do j=1,10
	    do k=iBegYr,iCurYr
	      do l=1,12
	        rData(i,j,k,l) = -99.99
	      enddo
	    enddo
	  enddo
	enddo

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

c	Read in all data

	write (*,*) ' '
	write (*,*) 'Reading in pdsi05dv.'
	write (*,*) ' '
  
       cFile='../data/pdsi05dv.'
	open(unit=12,file=cFile)
5	read (12,10,end=50) ist,idv,iYear,(rVal(i),i=1,13)
10	format (3x,i2,i2,2x,i4,13(f7.2))

	if((iYear.lt.iBegYr).or.(iYear.gt.iCurYr).or.
     +   (ist.gt.48)) goto 5

	do iMonth=1,12
	  rData(ist,idv,iYear,iMonth) = rVal(iMonth)
	enddo
	 
	goto 5
	
50	close(12)

************************************************************
*	Find out which divisions are legitimate. 
************************************************************

	do i=1,48
	  do j=1,10
	    iExist(i,j)=0
	  enddo
	enddo
	do i=1,48
	  do j=1,iDivisions(i)
	    iExist(i,j)=1
	  enddo
	enddo
	iExist(25,4)=0    ! No Nebraska Div 4

***************************************************************
*	Store divisional areas in an array:
*	Loop through all divisions and sum the areas for a
*	total U.S. area.
***************************************************************
	open(73,file='divarea.')
 11     read(73,12,end=13) iState,(rArea(l),l=1,11)
 12     format(i2,11(f7.0))
   
 	do iDiv=1,11 
	  rDumArea(iState,iDiv) = rArea(iDiv)
	enddo
        goto 11
 13     continue
  	close(73)

**********************************************************************
*  	Analysis for current month, 3-month, 6-month, 12-month and YTD
**********************************************************************

      Do iLoop=1,5

	write (*,*) 'Calculating for ...',cPeriod(iLoop)
	write (*,*) ' '

       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
	
	do i=iBegYr,iEndYr
	   rTotalArea(k) = 0.
	   rTotalAreaMoist(k) = 0.
	   rTotalAreaDrought(k) = 0.
	enddo

c  File for saving extreme divisions for each year:
	open(14,file='../data/extr-pdsi-divs.'//cBegMo//'-'//cEndMo//
     +               '.txt')
***************************************************************
*	Create an array that has a mean annual pdsi value for
*	a division
***************************************************************

	do i=1,48
	  do j=1,10
	   if (iExist(i,j).eq.1) then
	    iCountYears=0	
	    do k=iBegYr,iEndYr
	      iCountYears=iCountYears+1
	      iCount=0
*******************************************************
*       For individual month calculations
*******************************************************
              if(iBegMo.eq.iEndMo) then
	        rLocalData(iCountYears) = rData(i,j,k,iBegMo)
	        rAnnual(k)=rData(i,j,k,iBegMo)
*******************************************************
*       For Winter or other year-spanning calculations
*******************************************************
              elseif(iBegMo.gt.iEndMo) then
                do l=iBegMo,12
	          iCount=iCount+1
	          rDummyData(iCount) = rData(i,j,k-1,l)
	        enddo
                do l=1,iEndMo
	          iCount=iCount+1
	          rDummyData(iCount) = rData(i,j,k,l)
	        enddo
	        call average(iCount,-99.99,rDummyData,rAve)
	        rLocalData(iCountYears) = rAve
	        rAnnual(k)=rAve
*******************************************************
*       For other seasonal/annual or month combinations
*******************************************************
              else
                do l=iBegMo,iEndMo
	          iCount=iCount+1
	          rDummyData(iCount) = rData(i,j,k,l)
	        enddo
	        call average(iCount,-99.99,rDummyData,rAve)
	        rLocalData(iCountYears) = rAve
	        rAnnual(k)=rAve
		endif
	    enddo

***************************************************************
            call sort(iCountYears,rLocalData)

            rSevereMoist=
     +          quickpile(rLocalData,iCountYears,rHiPrcnt,-99.99)
            rSevereDrought=
     +          quickpile(rLocalData,iCountYears,rLoPrcnt,-99.99)

***************************************************************
*	Loop through all years and find those years where values
*	are greater than or equal to 3 (severe moisture surplus) 
*	and less than or equal to -3 (severe drought). 
*	Sum up total area and partial areas meeting criteria.
***************************************************************

	do k=iBegYr,iEndYr
	   rTotalArea(k) = rTotalArea(k) + rDumArea(i,j)
	      if (rAnnual(k).ge.rSevereMoist) then
	        rTotalAreaMoist(k) = 
     +                   rTotalAreaMoist(k) + rDumArea(i,j)
	         write(14,444) i,j,k
444	         format(i2.2,1x,i2.2,1x,'90',1x,i4)
	      endif
	      if (rAnnual(k).le.rSevereDrought) then
	        rTotalAreaDrought(k) = 
     +                   rTotalAreaDrought(k) + rDumArea(i,j)
	         write(14,555) i,j,k
555	         format(i2.2,1x,i2.2,1x,'10',1x,i4)
	      endif
	    enddo
	   endif
	  enddo
	enddo

	close(14)
***************************************************************
c 	Calculate Area of Contiguous U.S. with meeting much above and
c	much below normal criteria.

	do i=iBegYr,iEndYr
	  rPercentAreaMoist(i) = 
     +           rTotalAreaMoist(i)/rTotalArea(i)
	enddo
	do i=iBegYr,iEndYr
	  rPercentAreaDrought(i) = 
     +           rTotalAreaDrought(i)/rTotalArea(i)
	enddo

*************************************************************** 
c	Write out final results and multiply by two to get twice the area

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

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

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

	if(iLoop.eq.5) then
         cOutFile='../data/dk-step3.ytd.results'
         cOutFile1='../data/dk-step3-lo.ytd.results'
         cOutFile2='../data/dk-step3-hi.ytd.results'

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

	 do i=iBegYr,iEndYr
	  rFinal=rPercentAreaMoist(i)+
     +           rPercentAreaDrought(i)
	  write (13,14) i,rFinal
	  write (14,14) i,rPercentAreaDrought(i)
	  write (15,14) i,rPercentAreaMoist(i)
	 enddo

	 close(13)	  	
	 close(14)	  	
	 close(15)	  	
	endif

c  Now, for regional results:
      Enddo

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

	end

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

	subroutine divarea(i,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 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
