MODULE COMMpreCan
      IMPLICIT NONE
      SAVE
						 
       INTEGER, PARAMETER :: nXPoint=355 
       INTEGER, PARAMETER :: nYPoint=390
       INTEGER, PARAMETER :: BASESYEAR=1986
       INTEGER, PARAMETER :: BASEEYEAR=2005
       INTEGER, PARAMETER :: SSYEAR=1961
       INTEGER, PARAMETER :: EEYEAR=2005
       							
      character(20) :: STNID
      integer(4)    :: STDSPAN, PRCPNN
      integer(4)    :: SYEAR, EYEAR, TOT, YRS, BYRS, WINSIZE, SS

      real          :: LATITUDE, PRCP(500*365)
      real          :: prin(nXPoint,nYPoint)
      real          :: txin(nXPoint,nYPoint)
      real          :: tnin(nXPoint,nYPoint)
      real          :: TMAX(500*365), TMIN(500*365), MISSING
      integer(4)    :: YMD(500*365,3), MNASTAT(500,12,3),YNASTAT(500,3)
      integer(4)    :: MON(12),MONLEAP(12)

      INTEGER       :: ierr,num_procs,mype,js_x,je_x

      data MON/31,28,31,30,31,30,31,31,30,31,30,31/
      data MONLEAP/31,29,31,30,31,30,31,31,30,31,30,31/
      data WINSIZE/5/

      END MODULE COMMpreCan


C########################
C#  Main program start  #
C########################

      use COMMpreCan
						
      INCLUDE "mpif.h" 
													
      character(140) :: ifile
      character(160) :: ifile_dates      
      character(160) :: ifile_pr, ifile_tx, ifile_tn
      character(160) :: omissf, otmpfile, ofile	      

      character(90) :: header
      integer(4)    :: stnnum
      integer       :: rno
						
      CHARACTER(LEN=255) :: PATH_out
      CHARACTER(LEN=3)   :: xpoint 
      CHARACTER(LEN=3)   :: ypoint 
						
      INTEGER            :: j 
      INTEGER            :: i
      INTEGER            :: TOTDIAS
      CALL MPI_INIT(ierr)
      call mpi_comm_size(mpi_comm_world,num_procs, ierr )
      call mpi_comm_rank(mpi_comm_world,mype, ierr )
				
*      print*,   num_procs, mype
 

      MISSING=-99.9
      SS=int(WINSIZE/2)

CDiego Calcula o numero de dias para ler o arquivao
      ANOBI=0
      TOTANOS=EEYEAR-SSYEAR+1
      do ano=SSYEAR,EEYEAR
        if(mod(ano,4.).eq.0) then
           ANOBI=ANOBI+1
        endif
      enddo
      
      TOTDIAS=(TOTANOS*365)+ANOBI

C     Algumas variaveis do para.txt						
      STDSPAN=4
      PRCPNN=25
						
      call para_range(1,nXPoint,num_procs,mype,js_x,je_x)  
      print *, ' task id, jsta, end = ',mype ,js_x,je_x

      stnnum=1				
      DO i=js_x,je_x
      WRITE(xpoint,'(I3.3)')i
						
      PATH_out='/scratchout/grupos/grpetamc/home/andre.lyra/'
     +//'CORDEX_Group1/Eta_SAM20/CanESM2/'
     +//'fclimdex_ascii/'//TRIM(xpoint)//'/'
      
      CALL SYSTEM ( "mkdir -p "//PATH_out ) 
      
      rno=1					
      ifile=TRIM(PATH_out)//TRIM(xpoint)//'.txt'  
      
      ifile_dates='/scratchout/grupos/grpetamc//home/andre.lyra/'//
     +'CORDEX_Group1/Eta_SAM20/CanESM2/dates.txt'
     
      ifile_pr='/scratchout/grupos/grpetamc//home/andre.lyra/'//
     +'CORDEX_Group1/Eta_SAM20/CanESM2/pr.dat'
      ifile_tx='/scratchout/grupos/grpetamc//home/andre.lyra/'//
     +'CORDEX_Group1/Eta_SAM20/CanESM2/tasmax.dat'
      ifile_tn='/scratchout/grupos/grpetamc//home/andre.lyra/'//
     +'CORDEX_Group1/Eta_SAM20/CanESM2/tasmin.dat'
      
      open(10, file=ifile_dates, STATUS="OLD", IOSTAT=ios10)
      
      open(11, file=ifile_pr, STATUS='OLD', FORM='UNFORMATTED',
     + ACCESS='direct', RECL=nXPoint*nYPoint*4, IOSTAT=ios11)
      open(12, file=ifile_tx, STATUS='OLD', FORM='UNFORMATTED',
     + ACCESS='direct', RECL=nXPoint*nYPoint*4, IOSTAT=ios12)
      open(13, file=ifile_tn, STATUS='OLD', FORM='UNFORMATTED',
     + ACCESS='direct', RECL=nXPoint*nYPoint*4, IOSTAT=ios13)

      if(ios10.ne.0) then
        write(6,*) "ERROR during opening file: ", trim(ifile_dates)
        write(6,*) "Program STOP!!"
        stop
      endif

      if(ios11.ne.0 .or. ios12.ne.0 .or. ios13.ne.0) then
        write(6,*) "ERROR during opening file: ", trim(ifile_pr)
        write(6,*) "ERROR during opening file: ", trim(ifile_pr)
        write(6,*) "ERROR during opening file: ", trim(ifile_pr)
        write(6,*) "Program STOP!!"
        stop
      endif

      
      otmpfile=trim(ifile)//'_prcpQC'
      open(81, file=otmpfile)	
      otmpfile=trim(ifile)//'_tempQC'
      open(82, file=otmpfile)

      omissf=trim(ifile)//"_NASTAT"
      open(20,file=trim(omissf))
      
C    data chrtmp/"FD","SU","ID","TR"/
      ofile=trim(ifile)//"_FD"
      open(31,file=ofile)
      ofile=trim(ifile)//"_SU"      
      open(32,file=ofile)
      ofile=trim(ifile)//"_ID"      
      open(33,file=ofile)
      ofile=trim(ifile)//"_TR"      
      open(34,file=ofile)
      
      ofile=trim(ifile)//"_GSL"
      open(35,file=ofile)

C    data chrtmp/"TXx","TXn","TNx","TNn"/      
      ofile=trim(ifile)//"_TXx"
      open(36,file=ofile)
      ofile=trim(ifile)//"_TXn"
      open(37,file=ofile)
      ofile=trim(ifile)//"_TNx"      
      open(38,file=ofile)
      ofile=trim(ifile)//"_TNn"      
      open(39,file=ofile)
      
      ofile=trim(ifile)//"_DTR"
      open(70,file=ofile)

C    data chrtmp/"R10mm","R20mm","Rnnmm"/
      ofile=trim(ifile)//"_R10mm"
      open(40,file=ofile)
      ofile=trim(ifile)//"_R20mm"
      open(41,file=ofile)
      ofile=trim(ifile)//"_Rnnmm"      
      open(42,file=ofile)
      ofile=trim(ifile)//"_SDII"  
      open(43,file=ofile)
   
      ofile=trim(ifile)//"_RX1day"      
      open(44,file=ofile)
      ofile=trim(ifile)//"_RX5day"  
      open(45,file=ofile)

      ofile=trim(ifile)//"_CDD"      
      open(46,file=ofile)
      ofile=trim(ifile)//"_CWD"  
      open(47,file=ofile)

      ofile=trim(ifile)//"_PRCPTOT"
      open(48,file=ofile)
      ofile=trim(ifile)//"_R95p"
      open(49,file=ofile)
      ofile=trim(ifile)//"_p95"      
      open(50,file=ofile)
      ofile=trim(ifile)//"_R99p"  
      open(51,file=ofile)
      ofile=trim(ifile)//"_p99"      
      open(52,file=ofile)               

      ofile=trim(ifile)//"_thresan10"
      open(53,file=ofile)
      ofile=trim(ifile)//"_thresan90"
      open(54,file=ofile)
      ofile=trim(ifile)//"_thresax10"      
      open(55,file=ofile)
      ofile=trim(ifile)//"_thresax90"  
      open(56,file=ofile)

      ofile=trim(ifile)//"_thresbn10"
      open(57,file=ofile)
      ofile=trim(ifile)//"_thresbn90"
      open(58,file=ofile)
      ofile=trim(ifile)//"_thresbx10"      
      open(59,file=ofile)
      ofile=trim(ifile)//"_thresbx90"  
      open(60,file=ofile)

      ofile=trim(ifile)//"_TX90p"
      open(61,file=ofile)
      ofile=trim(ifile)//"_TX50p"
      open(62,file=ofile)
      ofile=trim(ifile)//"_TX10p"      
      open(63,file=ofile)

      ofile=trim(ifile)//"_TN90p"
      open(64,file=ofile)
      ofile=trim(ifile)//"_TN50p"
      open(65,file=ofile)
      ofile=trim(ifile)//"_TN10p"      
      open(66,file=ofile)

      ofile=trim(ifile)//"_WSDI"
      open(67,file=ofile)
      ofile=trim(ifile)//"_CSDI"      
      open(68,file=ofile)    
        
      BYRS=BASEEYEAR-BASESYEAR+1

      DO j=1,nYPoint
      WRITE(ypoint,'(I3.3)')j
      
CAndre      IF (mype.eq.0) print*, "j, rno :", j, rno	
     
      LATITUDE = (0.2*(j-1)) - 50.0
						
      STNID=TRIM(xpoint)//TRIM(ypoint) 
CAndre      print *, "STNID", STNID
      
      call qc(ifile,TOTDIAS,mype,rno,i,j)
C      call FD(ifile)    ! FD, SU, ID, TR
C      call GSL(ifile)   ! GSL
      call TXX(ifile)   ! TXx, TXn, TNx, TNn, DTR
C      call Rnnmm(ifile) ! R10mm, R20mm, Rnnmm, SDII
      call RX5day(ifile)! Rx1day, Rx5day
      call CDD(ifile,TOTDIAS)   ! CDD, CWD
      call R95p(ifile)  ! R95p, R99p, PRCPTOT
      call TX10p(ifile) ! TX10p, TN10p, TX90p, TN90p
      
      stnnum=stnnum+1

      ENDDO !Y 
      
      ENDDO !X 
      
      CALL MPI_FINALIZE(ierr)
						
      stnnum=stnnum-1
      IF (mype.eq.0) write (*,*) ""
      IF (mype.eq.0) write(6,*) "Total ",stnnum,"stations be calculated"
      
						
      END PROGRAM
						
						
						
      integer function leapyear(iyear)
      integer iyear

      if(mod(iyear,400).eq.0) then
        leapyear=1
      else
        if(mod(iyear,100).eq.0) then
          leapyear=0
        else
          if(mod(iyear,4).eq.0) then
           leapyear=1
          else
            leapyear=0
        endif
       endif
      endif

      end

      real function percentile(x, length, per)
      use COMMpreCan
      integer length
      real x(length), per
      real xtos(length),bb,cc
      integer nn
      
      if(per.gt.1.or.per.lt.0) then
        print *, "Function percentile return error: parameter perc"
        stop
      endif

      nn=0
      do i=1, length
        if(x(i).ne.MISSING)then
          nn=nn+1
          xtos(nn)=x(i)
        endif
      enddo

      if(nn.eq.0) then
        percentile=MISSING
      else
        call sort(nn,xtos)
        bb=nn*per+per/3.+1/3.
        cc=real(int(bb))
        if(int(cc).ge.nn) then
          percentile=xtos(nn)
        else
          percentile=xtos(int(cc))+(bb-cc)*
     &          (xtos(int(cc)+1)-xtos(int(cc)))
        endif
      endif

      end

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE PARA_RANGE (N1,N2,NPROCS,IRANK,ISTA,IEND)

C USAGE:    CALL COLLECT(A)
C   INPUT ARGUMENT LIST:
C     N1 - FIRST INTERATE VALUE
C     N2 - LAST INTERATE VALUE
C     NPROCS - NUMBER OF MPI TASKS
C     IRANK - MY TAKS ID
C   OUTPUT ARGUMENT LIST:
C     ISTA - FIRST LOOP VALUE
C     IEND - LAST LOOP VALUE

      implicit none
      integer n1,n2,nprocs,irank,ista,iend
      integer iwork1, iwork2
      iwork1 = ( n2 - n1 + 1 ) / nprocs
      iwork2 = mod ( n2 - n1 + 1, nprocs )
C     print*, irank, iwork1, iwork2
      ista = irank * iwork1 + n1 + min ( irank, iwork2 )
      iend = ista + iwork1 - 1
      if ( iwork2 .gt. irank ) iend = iend + 1

      END SUBROUTINE PARA_RANGE

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
c---Sorts an array arr(1:n) into ascending numerical order using the Quicksort
c   algorithm. n is inpu; arr is replace on output by its sorted rearrangement.
c   Parameters: M is the size of subarrays sorted by straight insertion
c   and NSTACK is the required auxiliary.
      SUBROUTINE sort(n,arr)
      INTEGER n,M,NSTACK
      REAL arr(n)
      PARAMETER (M=7,NSTACK=50)
      INTEGER i,ir,j,jstack,k,l,istack(NSTACK)
      REAL a,temp
      jstack=0
      l=1
      ir=n
1     if(ir-l.lt.M)then
        do 12 j=l+1,ir
          a=arr(j)
          do 11 i=j-1,1,-1
            if(arr(i).le.a)goto 2
            arr(i+1)=arr(i)
11        continue
          i=0
2         arr(i+1)=a
12      continue
        if(jstack.eq.0)return
        ir=istack(jstack)
        l=istack(jstack-1)
        jstack=jstack-2
      else
        k=(l+ir)/2
        temp=arr(k)
        arr(k)=arr(l+1)
        arr(l+1)=temp
        if(arr(l+1).gt.arr(ir))then
          temp=arr(l+1)
          arr(l+1)=arr(ir)
          arr(ir)=temp
        endif
        if(arr(l).gt.arr(ir))then
          temp=arr(l)
          arr(l)=arr(ir)
          arr(ir)=temp
        endif
        if(arr(l+1).gt.arr(l))then
          temp=arr(l+1)
          arr(l+1)=arr(l)
          arr(l)=temp
        endif
        i=l+1
        j=ir
        a=arr(l)
3       continue
          i=i+1
        if(arr(i).lt.a)goto 3
4       continue
          j=j-1
        if(arr(j).gt.a)goto 4
        if(j.lt.i)goto 5
        temp=arr(i)
        arr(i)=arr(j)
        arr(j)=temp
        goto 3
5       arr(l)=arr(j)
        arr(j)=a
        jstack=jstack+2
        if(jstack.gt.NSTACK)pause 'NSTACK too small in sort'
        if(ir-i+1.ge.j-l)then
          istack(jstack)=ir
          istack(jstack-1)=i
          ir=j-1
        else
          istack(jstack)=j-1
          istack(jstack-1)=l
          l=i
        endif
      endif
      goto 1
      END

C  (C) Copr. 1986-92 Numerical Recipes Software &#5,.
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

CCC: OBS- 366*n. de anos no to ascii*390,3
      subroutine qc(ifile,TOTDIAS,mype,rno,xi,yj)
      use COMMpreCan
      character*140 ifile
      character*140 omissf, otmpfile
      character*140  title(3)
      integer ios, rno, tmpymd(366*45*390,3), i, ith
      real tmpdata(366*45*390,3),stddata(365,500,3)
      real stdval(365,3),m1(365,3)
      integer kth,month,k,trno,ymiss(3),mmiss(3),stdcnt(3),
     &        missout(500,13,3),tmpcnt
      integer TOTDIAS, mype
      integer xi, yj     

      data title/"PRCPMISS","TMAXMISS","TMINMISS"/
      
      omissf=trim(ifile)//"_NASTAT"

CAndre     print*, BASESYEAR, BASEEYEAR
CAndre      open(10, file=ifile, STATUS="OLD", IOSTAT=ios)
CAndre     print *, ifile, ios

CAndre    if(ios.ne.0) then
CAndre      write(6,*) "ERROR during opening file: ", trim(ifile)
CAndre      write(6,*) "Program STOP!!"
CAndre      stop
CAndre    endif

      otmpfile=trim(ifile)//'_prcpQC'
CAndre      open(81, file=otmpfile)	
CAndre      print*,otmpfile 
      write(81,*) "PRCP Quality Control Log File:"
      otmpfile=trim(ifile)//'_tempQC'
CAndre      open(82, file=otmpfile)
CAndre      print*,otmpfile 
      write(82,*) "TMAX and TMIN Quality Control Log File:"

CAndre      rno=1
CAndre 88   read(10,*,end=110) (tmpymd(rno,j),j=1,3),
CAndre     &                   (tmpdata(rno,j),j=1,3)

      do it=1,TOTDIAS

         read(10,*) (tmpymd(rno,j),j=1,3)

         read(11,rec=it) prin(:,:)
         read(12,rec=it) txin(:,:)	
         read(13,rec=it) tnin(:,:)	
	
	 tmpdata(rno,1)=prin(xi,yj)
	 tmpdata(rno,2)=txin(xi,yj)
	 tmpdata(rno,3)=tnin(xi,yj)
           
	 rno=rno+1
      enddo

CAndre      rno=rno-1      
CAndre      close(10)

      SYEAR=tmpymd(1,1)
      EYEAR=tmpymd(rno-1,1)
      YRS=EYEAR-SYEAR+1


      TOT=0
      do i=SYEAR,EYEAR
        do month=1,12
          if(leapyear(i)==1) then
                kth=MONLEAP(month)
          else
                kth=MON(month)
          endif
          do k=1,kth
            TOT=TOT+1
            YMD(tot,1)=i
            YMD(tot,2)=month
            YMD(tot,3)=k
          enddo
        enddo
      enddo

CAndre      j=1
      j=rno-TOTDIAS
      do i=1, TOT
111     if(YMD(i,1)*10000+YMD(i,2)*100+YMD(i,3).eq.
     &     tmpymd(j,1)*10000+tmpymd(j,2)*100+tmpymd(j,3)) then
          PRCP(i)=tmpdata(j,1)
          TMAX(i)=tmpdata(j,2)
          TMIN(i)=tmpdata(j,3)
          j=j+1
        elseif(YMD(i,1)*10000+YMD(i,2)*100+YMD(i,3).lt.
     &     tmpymd(j,1)*10000+tmpymd(j,2)*100+tmpymd(j,3)) then
          PRCP(i)=MISSING
          TMAX(i)=MISSING
          TMIN(i)=MISSING
        elseif(YMD(i,1)*10000+YMD(i,2)*100+YMD(i,3).gt.
     &  tmpymd(rno-1,1)*10000+tmpymd(rno-1,2)*100+tmpymd(rno-1,3)) then
          PRCP(i)=MISSING
          TMAX(i)=MISSING
          TMIN(i)=MISSING
        else
          j=j+1
          goto 111
        endif
      enddo

      trno=0
      MNASTAT=0
      YNASTAT=0
      do i=SYEAR,EYEAR
        ymiss=0
        stdcnt=0
        do month=1,12
          mmiss=0
          if(leapyear(i)==1) then
            kth=MONLEAP(month)
          else
            kth=MON(month)
          endif
          do k=1,kth
            trno=trno+1
            if(TMAX(trno).le.MISSING) TMAX(trno)=MISSING
            if(TMIN(trno).le.MISSING) TMIN(trno)=MISSING
            if(PRCP(trno).le.MISSING) PRCP(trno)=MISSING
	    
            if(TMAX(trno).lt.TMIN(trno).and.TMAX(trno).ne.MISSING
     &                   .and.TMIN(trno).ne.MISSING) then
              TMAX(trno)=MISSING
              TMIN(trno)=MISSING
              write(82, *) i*10000+month*100+k, "TMAX<TMIN!!"
            endif

            if(month.ne.2.or.k.ne.29) then
              stdcnt=stdcnt+1
              stddata(stdcnt, i-SYEAR+1, 1)=PRCP(trno)
              stddata(stdcnt, i-SYEAR+1, 2)=TMAX(trno)
              stddata(stdcnt, i-SYEAR+1, 3)=TMIN(trno)
            endif

            if((TMAX(trno).lt.-70..or.TMAX(trno).gt.70.).and.
     &         TMAX(trno).ne.MISSING) then
              TMAX(trno)=MISSING
              write(82, *) i*10000+month*100+k, "TMAX over bound!!"
            endif
            if((TMIN(trno).lt.-70..or.TMIN(trno).gt.70.).and.
     &         TMIN(trno).ne.MISSING) then
              TMIN(trno)=MISSING
              write(82, *) i*10000+month*100+k, "TMIN over bound!!"
            endif
            if(PRCP(trno).lt.0.and.PRCP(trno).ne.MISSING) then
              PRCP(trno)=MISSING
              write(81,*) i*10000+month*100+k, "PRCP less then 0!!"
            endif
            if(PRCP(trno).eq.MISSING) then
              mmiss(1)=mmiss(1)+1
              ymiss(1)=ymiss(1)+1
            endif
            if(TMAX(trno).eq.MISSING) then
              mmiss(2)=mmiss(2)+1
              ymiss(2)=ymiss(2)+1
            endif
            if(TMIN(trno).eq.MISSING) then
              mmiss(3)=mmiss(3)+1
              ymiss(3)=ymiss(3)+1
            endif
          enddo
          do k=1,3
            missout(i-SYEAR+1,month,k)=mmiss(k)
            if (mmiss(k).gt.3) then
              MNASTAT(i-SYEAR+1,month,k)=1
            endif
          enddo
        enddo
        do k=1,3
          missout(i-SYEAR+1,13,k)=ymiss(k)
          if(ymiss(k).gt.15) then
            YNASTAT(i-SYEAR+1,k)=1
          endif
        enddo
      enddo

c     do i=1,YRS
c     print *,(YNASTAT(i,k),k=1,3)
c     enddo

C Calculate STD for PRCP, TMAX and TMIN; then figure out outliers
c     stdval=0.
      do i=1,365
        do j=1,3
          stdval(i,j)=0.
        enddo
      enddo
      m1=0.
      do i=1,365
        stdcnt=0
        do j=1,YRS
          do k=2,3
            if(stddata(i,j,k).ne.MISSING) then
              stdcnt(k)=stdcnt(k)+1
              m1(i,k)=m1(i,k)+stddata(i,j,k)
            endif
          enddo
        enddo
        do k=2,3
          if(stdcnt(k).gt.0) then
            m1(i,k)=m1(i,k)/real(stdcnt(k))
          endif
        enddo
        stdtmp=0.
        do j=1,YRS
          do k=2,3
            if(stdcnt(k).gt.2.and.stddata(i,j,k).ne.MISSING) then
               stdval(i,k)=stdval(i,k)+
     &         (stddata(i,j,k)-m1(i,k))**2./(real(stdcnt(k))-1.)
            endif
          enddo
        enddo

        do k=2,3
          if(stdcnt(k).gt.2) then
            stdval(i,k)=stdval(i,k)**0.5
          else 
            stdval(i,k)=MISSING
          endif
        enddo
      enddo

      trno=0
      do i=SYEAR,EYEAR
        tmpcnt=0
        do month=1,12
          if(leapyear(i)==1) then
            kth=MONLEAP(month)
          else 
            kth=MON(month)
          endif
          do k=1, kth
            trno=trno+1
            if(month.ne.2.or.k.ne.29) tmpcnt=tmpcnt+1
c           if(abs(PRCP(trno)-m1(tmpcnt,1)).gt.
c    &          stdval(tmpcnt,1)*STDSPAN.and.PRCP(trno).ne.MISSING)
c    &      write(81, *) "Outlier: ", i, month, k, "PRCP: ", PRCP(trno),
c    &            "Lower limit:",m1(tmpcnt,1)-stdval(tmpcnt,1)*STDSPAN,
c    &            "Upper limit:",m1(tmpcnt,1)+stdval(tmpcnt,1)*STDSPAN
            if(stdval(tmpcnt,2).ne.MISSING)then
              if(abs(TMAX(trno)-m1(tmpcnt,2)).gt.
     &          stdval(tmpcnt,2)*STDSPAN.and.TMAX(trno).ne.MISSING)
     &      write(82, *) "Outlier: ", i, month, k, "TMAX: ", TMAX(trno),
     &            "Lower limit:",m1(tmpcnt,2)-stdval(tmpcnt,2)*STDSPAN,
     &            "Upper limit:",m1(tmpcnt,2)+stdval(tmpcnt,2)*STDSPAN
            endif
            if(stdval(tmpcnt,3).ne.MISSING)then
              if(abs(TMIN(trno)-m1(tmpcnt,3)).gt.
     &          stdval(tmpcnt,3)*STDSPAN.and.TMIN(trno).ne.MISSING)
     &      write(82, *) "Outlier: ", i, month, k, "TMIN: ", TMIN(trno),
     &            "Lower limit:",m1(tmpcnt,3)-stdval(tmpcnt,3)*STDSPAN,
     &            "Upper limit:",m1(tmpcnt,3)+stdval(tmpcnt,3)*STDSPAN
            endif
          enddo ! end do day
        enddo !end do month
      enddo ! end do year

C      open(20,file=trim(omissf))
      do i=SYEAR,EYEAR
        do k=1,3
          write(20,'(i4,2x,a8,13i4)')
     &    i,title(k),(missout(i+1-SYEAR,j,k),j=1,13)
        enddo
      enddo
c      close(20)
c      close(81)
c      close(82)

C  QC part finished, prepared data set: YMD(3), PRCP, TMAX & TMIN
C  and NASTAT dataset for missing values monthly and annual
      end

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      subroutine FD(ifile)
      use COMMpreCan
      character*140 ifile

      integer year, trno, kth, month 
      real oout(YRS,4)
      character*2 chrtmp(4)
      character*140 ofile
C oout(,1)--FD, oout(,2)--SU, oout(,3)--ID, oout(,4)--TR      
      data chrtmp/"FD","SU","ID","TR"/

      trno=0
      oout=0
      do i=1,YRS
        year=i+SYEAR-1
        do month=1,12
          if(leapyear(year)==1) then
                  kth=MONLEAP(month)
          else
                  kth=MON(month)
          endif
          do day=1,kth
            trno=trno+1
            if(YMD(trno,3).ne.day) then
              print *, 'ERROR1 at FD!!!'
              stop
            endif
            if(TMIN(trno).ne.MISSING.and.TMIN(trno).lt.0) 
     &          oout(i,1)=oout(i,1)+1
CCC: OBS- Limiar da variavel SU
            if(TMAX(trno).ne.MISSING.and.TMAX(trno).gt.25) 
     &          oout(i,2)=oout(i,2)+1
            if(TMAX(trno).ne.MISSING.and.TMAX(trno).lt.0) 
     &          oout(i,3)=oout(i,3)+1
            if(TMIN(trno).ne.MISSING.and.TMIN(trno).gt.20) 
     &          oout(i,4)=oout(i,4)+1
          enddo
        enddo
      enddo

      do i=1,YRS
        if(YNASTAT(i,2)==1) then
                oout(i,2)=MISSING  ! SU
                oout(i,3)=MISSING  ! ID
        endif
        if(YNASTAT(i,3)==1) then
                oout(i,1)=MISSING  ! FD
                oout(i,4)=MISSING  ! TR
        endif
      enddo

C      do j=1,4
C        ofile=trim(ifile)//"_"//chrtmp(j)
	
C        open(22,file=ofile)
          
	  write(31, *) "year    ", chrtmp(1)
          do i=1,YRS
            write(31, '(i8,f8.1)') i+SYEAR-1, oout(i,1)
          enddo
	  
          write(32, *) "year    ", chrtmp(2)
          do i=1,YRS
            write(32, '(i8,f8.1)') i+SYEAR-1, oout(i,2)
          enddo
          
	  write(33, *) "year    ", chrtmp(3)
          do i=1,YRS
            write(33, '(i8,f8.1)') i+SYEAR-1, oout(i,3)
          enddo
          
	  write(34, *) "year    ", chrtmp(4)
          do i=1,YRS
            write(34, '(i8,f8.1)') i+SYEAR-1, oout(i,4)
          enddo
	  	  	  
C        close(22)
C      enddo

      end

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      subroutine GSL(ifile)
      use COMMpreCan
      character*140 ifile
      character*140 ofile

      integer year,cnt,kth,month,day,marks,marke

      real TG,oout,strt(YRS),ee(YRS)


      strt=MISSING
      ee=MISSING
      cnt=0
      do i=1,YRS
        year=i+SYEAR-1
        marks=0
        marke=0
        do month=1,6
          if(leapyear(year).eq.1) then
            kth=MONLEAP(month)
          else
            kth=MON(month)
          endif
          do day=1,kth
            cnt=cnt+1
            if(YMD(cnt,1)*10000+YMD(cnt,2)*100+YMD(cnt,3).ne.
     &         year*10000+month*100+day) then
              print*, 'date count ERROR in GSL!'
              print*, YMD(cnt,1)*10000+YMD(cnt,2)*100+YMD(cnt,3),
     &                year*10000+month*100+day
              stop
            endif
            if(TMAX(cnt).ne.MISSING.and.TMIN(cnt).ne.MISSING) then
              TG=(TMAX(cnt)+TMIN(cnt))/2.
            else
              TG=MISSING
            endif
            if(LATITUDE.lt.0) then
              if(TG.ne.MISSING.and.TG.lt.5.)then
                marke=marke+1
              else
                marke=0
              endif
              if(marke.ge.6.and.i.gt.1.and.ee(i-1).eq.MISSING)then
                ee(i-1)=cnt
              endif
            else
              if(TG.ne.MISSING.and.TG.gt.5.)then
                marks=marks+1
              else
                marks=0
              endif
              if(marks.ge.6.and.strt(i).eq.MISSING)then
                strt(i)=cnt
              endif
            endif
          enddo
        enddo 
        marks=0
        marke=0
        do month=7,12
          do day=1,MON(month)
            cnt=cnt+1
            if(TMAX(cnt).ne.MISSING.and.TMIN(cnt).ne.MISSING) then
              TG=(TMAX(cnt)+TMIN(cnt))/2.
            else
              TG=MISSING
            endif
            if(LATITUDE.lt.0) then
              if(TG.ne.MISSING.and.TG.gt.5.)then
                marks=marks+1
              else
                marks=0
              endif
              if(marks.ge.6.and.strt(i).eq.MISSING)then
                strt(i)=cnt
              endif
            else
              if(TG.ne.MISSING.and.TG.lt.5.)then
                marke=marke+1
              else
                marke=0
              endif
              if(marke.ge.6.and.ee(i).eq.MISSING)then
                ee(i)=cnt
              endif
            endif
          enddo
        enddo
      enddo

CAndre      ofile=trim(ifile)//"_GSL"
CAndre      open(22,file=ofile)
      write(35,*) "  year    gsl  "
      do i=1,YRS
        year=i+SYEAR-1
        if(strt(i).ne.MISSING.and.ee(i).ne.MISSING.and.
     &     YNASTAT(i,2).ne.1.and.YNASTAT(i,3).ne.1)then
          oout=ee(i)-strt(i)
        elseif(strt(i).eq.MISSING.or.ee(i).eq.MISSING) then
          oout=0.
        endif
        if(YNASTAT(i,2).eq.1.or.YNASTAT(i,3).eq.1) oout=MISSING
c       if(year.eq.1923) print *, year, YNASTAT(i,2),YNASTAT(i,3)
        write(35,'(i6,f8.1)') year, oout
      enddo
CAndre      close(22)
      end

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      subroutine TXX(ifile)
      use COMMpreCan
      character*140 ifile
      character*140 ofile
      character*3 chrtmp(4)

      integer year,month,day,kth,cnt,nn
      real oout(YRS,12,4),yout(YRS,4), dtr(YRS,13)

      data chrtmp/"TXx","TXn","TNx","TNn"/

      oout=MISSING
      dtr=0.
      cnt=0
      do i=1,YRS
        year=i+SYEAR-1
        do month=1,12
          if(leapyear(year)==1) then
            kth=MONLEAP(month)
          else
            kth=MON(month)
          endif
          nn=0
          do day=1,kth
            cnt=cnt+1
            if(TMAX(cnt).ne.MISSING.and.TMIN(cnt).ne.MISSING) then
              dtr(i,month)=dtr(i,month)+(TMAX(cnt)-TMIN(cnt))
              nn=nn+1
            endif
            if(TMAX(cnt).ne.MISSING.and.(oout(i,month,1).eq.MISSING.or.
     &         TMAX(cnt).gt.oout(i,month,1))) then
              oout(i,month,1)=TMAX(cnt) ! TXX
            endif
            if(TMAX(cnt).ne.MISSING.and.(oout(i,month,2).eq.MISSING.or.
     &         TMAX(cnt).lt.oout(i,month,2))) then
              oout(i,month,2)=TMAX(cnt) ! TXN
            endif
            if(TMIN(cnt).ne.MISSING.and.(oout(i,month,3).eq.MISSING.or.
     &         TMIN(cnt).gt.oout(i,month,3))) then
              oout(i,month,3)=TMIN(cnt) ! TNX
            endif
            if(TMIN(cnt).ne.MISSING.and.(oout(i,month,4).eq.MISSING.or.
     &         TMIN(cnt).lt.oout(i,month,4))) then
              oout(i,month,4)=TMIN(cnt) ! TNN
            endif
          enddo 
          if(nn.gt.0.and.MNASTAT(i,month,2).eq.0.and.
     &          MNASTAT(i,month,3).eq.0) then
            dtr(i,month)=dtr(i,month)/nn
          else
            dtr(i,month)=MISSING
          endif
          if(MNASTAT(i,month,2).eq.1)then
            oout(i,month,1)=MISSING
            oout(i,month,2)=MISSING
          endif
          if(MNASTAT(i,month,3).eq.1)then
            oout(i,month,3)=MISSING
            oout(i,month,4)=MISSING
          endif
        enddo
      enddo

      yout=MISSING
      do i=1,YRS
        nn=0
        do month=1,12
          if(oout(i,month,1).ne.MISSING.and.(yout(i,1).eq.MISSING.or.
     &       oout(i,month,1).gt.yout(i,1))) then
            yout(i,1)=oout(i,month,1)
          endif
          if(oout(i,month,2).ne.MISSING.and.(yout(i,2).eq.MISSING.or.
     &       oout(i,month,2).lt.yout(i,2))) then
            yout(i,2)=oout(i,month,2)
          endif
          if(oout(i,month,3).ne.MISSING.and.(yout(i,3).eq.MISSING.or.
     &       oout(i,month,3).gt.yout(i,3))) then
            yout(i,3)=oout(i,month,3)
          endif
          if(oout(i,month,4).ne.MISSING.and.(yout(i,4).eq.MISSING.or.
     &       oout(i,month,4).lt.yout(i,4))) then
            yout(i,4)=oout(i,month,4)
          endif
          if(dtr(i,month).ne.MISSING) then
            dtr(i,13)=dtr(i,13)+dtr(i,month)
            nn=nn+1
          endif
        enddo
        if(nn.gt.0.and.YNASTAT(i,2).eq.0.and.YNASTAT(i,3).eq.0) then
          dtr(i,13)=dtr(i,13)/nn
        else
          dtr(i,13)=MISSING
        endif
        if(YNASTAT(i,2).eq.1) then
          yout(i,1)=MISSING
          yout(i,2)=MISSING
        endif
        if(YNASTAT(i,3).eq.1) then
          yout(i,3)=MISSING
          yout(i,4)=MISSING
        endif
      enddo

CAndre      do k=1,4
CAndre        ofile=trim(ifile)//"_"//chrtmp(k)
CANdre        open(22,file=ofile)
        write(36, *) "  year  jan   feb   mar   apr   may   jun  ",
     &               " jul   aug   sep   oct   nov   dec annual"
        do i=1,YRS
          write(36,'(i6,13f9.1)') i+SYEAR-1,(oout(i,j,1),j=1,12),
     &                            yout(i,1)
        enddo

        write(37, *) "  year  jan   feb   mar   apr   may   jun  ",
     &               " jul   aug   sep   oct   nov   dec annual"
        do i=1,YRS
          write(37,'(i6,13f9.1)') i+SYEAR-1,(oout(i,j,2),j=1,12),
     &                            yout(i,2)
        enddo
	
        write(38, *) "  year  jan   feb   mar   apr   may   jun  ",
     &               " jul   aug   sep   oct   nov   dec annual"
        do i=1,YRS
          write(38,'(i6,13f9.1)') i+SYEAR-1,(oout(i,j,3),j=1,12),
     &                            yout(i,3)
        enddo
	
        write(39, *) "  year  jan   feb   mar   apr   may   jun  ",
     &               " jul   aug   sep   oct   nov   dec annual"
        do i=1,YRS
          write(39,'(i6,13f9.1)') i+SYEAR-1,(oout(i,j,4),j=1,12),
     &                            yout(i,4)
        enddo		
CAndre        close(22)
CAndre      enddo

CAndre      ofile=trim(ifile)//"_DTR"
CAndre      open(22,file=ofile)
      write(70, *) "  year  jan   feb   mar   apr   may   jun  ",
     &             " jul   aug   sep   oct   nov   dec annual"
      do i=1,YRS
        write(70,'(i6,13f7.2)') i+SYEAR-1,(dtr(i,j),j=1,13)
      enddo
CAndre      close(22)

      end

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      subroutine Rnnmm(ifile)
      use COMMpreCan
      character*140 ifile
      character*140 ofile
      character*5 chrtmp(3)
      integer year,month,day,kth,cnt,nn

      real oout(YRS,3),sdii(YRS)

      data chrtmp/"R10mm","R20mm","Rnnmm"/
      cnt=0
      oout=0.
      sdii=0.
      do i=1,YRS
        nn=0
        year=i+SYEAR-1
        do month=1,12
          if(leapyear(year).eq.1) then
            kth=MONLEAP(month)
          else
            kth=MON(month)
          endif
          do day=1,kth
            cnt=cnt+1
            if(PRCP(cnt).ge.1.) then
              sdii(i)=sdii(i)+PRCP(cnt)
              nn=nn+1
            endif
            if(PRCP(cnt).ge.10.) oout(i,1)=oout(i,1)+1.
            if(PRCP(cnt).ge.20.) oout(i,2)=oout(i,2)+1.
            if(PRCP(cnt).ge.PRCPNN) oout(i,3)=oout(i,3)+1.
          enddo
        enddo
        if(nn.gt.0) then
          sdii(i)=sdii(i)/nn
        endif
      enddo

      do i=1,YRS
        if(YNASTAT(i,1).eq.1) then
          do k=1,3
            oout(i,k)=MISSING
          enddo
          sdii(i)=MISSING
        endif
      enddo

CAndre      do k=1,3
CAndre        ofile=trim(ifile)//"_"//chrtmp(k)
CAndre        open(22,file=ofile)
        write(40,*) "  year  ",chrtmp(1)
        do i=1,YRS
          write(40,'(i6,f8.1)') i+SYEAR-1, oout(i,1)
        enddo
        write(41,*) "  year  ",chrtmp(2)
        do i=1,YRS
          write(41,'(i6,f8.1)') i+SYEAR-1, oout(i,2)
        enddo
        write(42,*) "  year  ",chrtmp(3)
        do i=1,YRS
          write(42,'(i6,f8.1)') i+SYEAR-1, oout(i,3)
        enddo
CAndre        close(22)
CAndre      enddo

CAndre      ofile=trim(ifile)//"_SDII"
CAndre      open(22,file=ofile)
      write(43,*) "  year    sdii"
      do i=1,YRS
        write(43,'(i6,f8.1)') i+SYEAR-1, sdii(i)
      enddo
CAndre      close(22)

      end

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      subroutine RX5day(ifile)
      use COMMpreCan
      character*140 ifile
      character*140 ofile

      integer year, month, day,cnt

      real r1(YRS,13), r5(YRS,13), r5prcp

      cnt=0
      r1=MISSING
      r5=MISSING
      do i=1,YRS
        year=i+SYEAR-1
        do month=1,12
          if(leapyear(year).eq.1) then
            kth=MONLEAP(month)
          else
            kth=MON(month)
          endif
          do day=1,kth
            cnt=cnt+1
c           if(year.eq.1904.and.month.eq.1) then
c             print *, year, month, day, PRCP(cnt)
c           endif
            if(cnt.gt.5)then
              r5prcp=0.
              do k=cnt-4,cnt
                if(PRCP(k).ne.MISSING)then
                  r5prcp=r5prcp+PRCP(k)
                endif
              enddo
            else
              r5prcp=MISSING
            endif
            if(PRCP(cnt).ne.MISSING.and.(r1(i,month).eq.MISSING
     &         .or.PRCP(cnt).gt.r1(i,month))) then
              r1(i,month)=PRCP(cnt)
            endif
            if(PRCP(cnt).ne.MISSING.and.r5prcp.gt.r5(i,month)) then
              r5(i,month)=r5prcp
            endif
          enddo
          if(MNASTAT(i,month,1).eq.1) then
            r1(i,month)=MISSING
            r5(i,month)=MISSING
          endif
          if(r1(i,month).ne.MISSING.and.(r1(i,13).eq.MISSING
     &       .or.r1(i,month).gt.r1(i,13))) then
            r1(i,13)=r1(i,month)
          endif
          if(r5(i,month).ne.MISSING.and.(r5(i,13).eq.MISSING
     &       .or.r5(i,month).gt.r5(i,13))) then
            r5(i,13)=r5(i,month)
          endif
        enddo
        if(YNASTAT(i,1).eq.1) then
          r1(i,13)=MISSING
          r5(i,13)=MISSING
        endif
      enddo

CAndre      ofile=trim(ifile)//"_RX1day"
CAndre      open(22,file=ofile)
      write(44, *) "  year  jan    feb    mar    apr    may    jun  ",
     &             "  jul    aug    sep    oct    nov    dec  annual"
      do i=1,YRS
        write(44, '(i6,13f7.1)') i+SYEAR-1,(r1(i,j),j=1,13)
      enddo
CAndre      close(22)

CAndre      ofile=trim(ifile)//"_RX5day"
CAndre      open(22,file=ofile)
      write(45, *) "  year  jan    feb    mar    apr    may    jun  ",
     &             "  jul    aug    sep    oct    nov    dec  annual"
      do i=1,YRS
        write(45, '(i6,13f7.1)') i+SYEAR-1,(r5(i,j),j=1,13)
      enddo
CAndre      close(22)

      end

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      subroutine CDD(ifile,TOTDIAS)
      use COMMpreCan
      character*140 ifile
      character*140 ofile

      integer year, month, day, kth, cnt

      real ocdd(YRS), ocwd(YRS), nncdd, nncwd

      cnt=0
      ocdd=0.
      ocwd=0.
      do i=1,YRS
        if(i==1)nncdd=0.
        if(i==1)nncwd=0.
        year=i+SYEAR-1
        do month=1,12
          if(leapyear(year).eq.1) then
            kth=MONLEAP(month)
          else
            kth=MON(month)
          endif
          do day=1,kth
            cnt=cnt+1
            if(PRCP(cnt).eq.MISSING) then
              nncdd=0.
              nncwd=0.
            elseif(PRCP(cnt).lt.1) then
              nncdd=nncdd+1.
              if(nncwd.gt.ocwd(i)) ocwd(i)=nncwd
              nncwd=0.
            else
              nncwd=nncwd+1.
              if(nncdd.gt.ocdd(i)) ocdd(i)=nncdd
              nncdd=0.
            endif
c           if(year.eq.1959.and.month.eq.12) then 
c                   print *, month,day,nncdd, ocdd(i)
c           endif
          enddo
        enddo

        if(ocwd(i).lt.nncwd) then
          if(year.eq.EYEAR) then
                  ocwd(i)=nncwd
          elseif(PRCP(cnt+1).lt.1..or.PRCP(cnt+1).eq.MISSING)then
                  ocwd(i)=nncwd
          endif
        endif

        if(ocdd(i).lt.nncdd) then
          if(year.eq.EYEAR) then
                  ocdd(i)=nncdd
	  	          if(nncdd.eq.TOTDIAS) then
			           ocdd(i)=MISSING
		          endif				  
          elseif(PRCP(cnt+1).ge.1..or.PRCP(cnt+1).eq.MISSING)then
                  ocdd(i)=nncdd
          endif

          if(ocdd(i).eq.0) ocdd(i)=MISSING

        endif

	if(leapyear(year).eq.1) then
             if(ocdd(i).ge.366) ocdd(i)=366
	     if(ocwd(i).ge.366) ocwd(i)=366
        else
             if(ocdd(i).ge.365) ocdd(i)=365
	     if(ocwd(i).ge.365) ocwd(i)=365
        endif

        if(YNASTAT(i,1).eq.1) then
          ocdd(i)=MISSING
          ocwd(i)=MISSING
        endif
      enddo

CAndre      ofile=trim(ifile)//"_CDD"
CAndre      open(22,file=ofile)
      write(46,*) "  year   cdd  "
      do i=1,YRS
        write(46,'(i6,f8.1)') i+SYEAR-1, ocdd(i)
      enddo
CAndre      close(22)

CAndre      ofile=trim(ifile)//"_CWD"
CAndre      open(22,file=ofile)
      write(47,*) "  year   cwd  "
      do i=1,YRS
        write(47,'(i6,f8.1)') i+SYEAR-1, ocwd(i)
      enddo
CAndre      close(22)

      end

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      subroutine R95p(ifile)
      use COMMpreCan
      character*140 ifile
      character*140 ofile
      integer year, month, day, kth,cnt,leng

      real r95out(YRS), prcptmp(TOT),r99out(YRS), prcpout(YRS), p95, p99

      cnt=0
      leng=0
      prcptmp=MISSING
      do i=1,YRS
        year=i+SYEAR-1
        do month=1,12
          if(leapyear(year).eq.1)then
            kth=MONLEAP(month)
          else
            kth=MON(month)
          endif
          do day=1,kth
            cnt=cnt+1
            if(year.ge.BASESYEAR.and.year.le.BASEEYEAR.and.
     &         PRCP(cnt).ne.MISSING.and.PRCP(cnt).ge.1.)then
              leng=leng+1
              prcptmp(leng)=PRCP(cnt)
            endif
          enddo
        enddo
      enddo
      p95=percentile(prcptmp,leng,0.95)
      p99=percentile(prcptmp,leng,0.99)

      cnt=0
      r95out=0.
      r99out=0.
      prcpout=0.
      do i=1,YRS
        year=i+SYEAR-1
        do month=1,12
          if(leapyear(year).eq.1)then
            kth=MONLEAP(month)
          else
            kth=MON(month)
          endif
          do day=1,kth
            cnt=cnt+1
            if(PRCP(cnt).ge.1..and.PRCP(cnt).ne.MISSING)then
              prcpout(i)=prcpout(i)+PRCP(cnt)
              if(PRCP(cnt).gt.p95) r95out(i)=r95out(i)+PRCP(cnt)
              if(PRCP(cnt).gt.p99) r99out(i)=r99out(i)+PRCP(cnt)
            endif
          enddo
        enddo
        if(YNASTAT(i,1).eq.1) then
          prcpout(i)=MISSING
          r95out(i)=MISSING
          r99out(i)=MISSING
        endif
      enddo

CAndre      ofile=trim(ifile)//"_PRCPTOT"
CAndre      open(22,file=ofile)
      write(48,*) "  year prcptot "
      do i=1,YRS
        write(48,'(i6,f8.1)') i+SYEAR-1, prcpout(i)
      enddo
CAndre      close(22)
CAndre      ofile=trim(ifile)//"_R95p"
CAndre      open(22,file=ofile)
      write(49,*) "  year  r95p  "
      do i=1,YRS
        write(49,'(i6,f8.1)') i+SYEAR-1, r95out(i)
      enddo
CAndre      close(22)
CAndre      ofile=trim(ifile)//"_p95"
CAndre      open(22,file=ofile)
       write(50,*) "  year  p95  "
      do i=1,YRS
        write(50,'(i6,f8.1)') i+SYEAR-1, p95
      enddo
CAndre      close(22)
CAndre      ofile=trim(ifile)//"_R99p"
CAndre      open(22,file=ofile)
      write(51,*) "  year  r99p  "
      do i=1,YRS
        write(51,'(i6,f8.1)') i+SYEAR-1, r99out(i)
      enddo
CAndre      close(22)
CAndre      ofile=trim(ifile)//"_p99"
CAndre      open(22,file=ofile)
      write(52,*) "  year  p99  "
      do i=1,YRS
        write(52,'(i6,f8.1)') i+SYEAR-1, p99
      enddo
CAndre      close(22)

      end

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      subroutine TX10p(ifile)
      use COMMpreCan
      character*140 ifile
      character*140 ofile

      integer year, month, day, kth, cnt, nn,  missxcnt, missncnt,
     &        iter, cntx, cntn,i,byear,flgtn,flgtx,flg

      real tmaxbase(TOT),tminbase(TOT),txdata(BYRS,365+2*SS),
     &     tndata(BYRS,365+2*SS),thresan10(365),txdtmp(BYRS,365),
     &     tndtmp(BYRS,365),tnboot(BYRS,365+2*SS),
     &     txboot(BYRS,365+2*SS),thresan90(365),thresax10(365),
     &     thresax90(365),tx10out(YRS,13),tx90out(YRS,13),
     &     thresax50(365),thresan50(365),tx50out(YRS,13),
     &     tn50out(YRS,13),thresbx50(365,BYRS,BYRS-1),
     &     thresbn50(365,BYRS,BYRS-1),
     &     tn10out(YRS,13),tn90out(YRS,13),thresbn90(365,BYRS,BYRS-1),
     &     thresbn10(365,BYRS,BYRS-1),thresbx90(365,BYRS,BYRS-1),
     &     thresbx10(365,BYRS,BYRS-1),threstmp(365),wsdi(YRS),csdi(YRS)
      logical ismiss,nomiss
      
      cnt=0
      nn=0
      txdtmp=MISSING
      tndtmp=MISSING
      do i=1,YRS
        year=i+SYEAR-1
        nn=0
        do month=1,12
          if(leapyear(year)==1) then
            kth=MONLEAP(month)
          else
            kth=MON(month)
          endif
          do day=1,kth
            cnt=cnt+1
           if(year.ge.BASESYEAR.and.year.le.BASEEYEAR.and.(month.ne.2
     &         .or.day.ne.29))then
              nn=nn+1
              txdtmp(i+SYEAR-BASESYEAR,nn)=TMAX(cnt)
              tndtmp(i+SYEAR-BASESYEAR,nn)=TMIN(cnt)
            endif
          enddo
        enddo
        if(year.ge.BASESYEAR.and.year.le.BASEEYEAR.and.nn.ne.365)then
          print *,year, BASESYEAR, BASEEYEAR
          print *,"date count error in TX10p!", nn
          stop
        endif
      enddo
      
      do i=1,BYRS
        do j=1,SS
          if(i.eq.1) then
            tndata(i,j)=tndtmp(i,1)
            txdata(i,j)=txdtmp(i,1)
          else 
            tndata(i,j)=tndtmp(i-1,365+j-SS)
            txdata(i,j)=txdtmp(i-1,365+j-SS)
          endif
        enddo
        do j=1,365
          tndata(i,j+SS)=tndtmp(i,j)
          txdata(i,j+SS)=txdtmp(i,j)
        enddo
        do j=1,SS
          if(i.eq.BYRS)then
            tndata(i,j+365+SS)=tndtmp(i,365)
            txdata(i,j+365+SS)=txdtmp(i,365)
          else
            tndata(i,j+365+SS)=tndtmp(i+1,j)
            txdata(i,j+365+SS)=txdtmp(i+1,j)
          endif
        enddo
      enddo
      
      thresan10=MISSING
      flgtn=0
      flgtx=0
      call threshold(tndata,.1,thresan10, flgtn)
      if(flgtn.eq.1) then
        write(6,*) "TMIN Missing value overflow in exceedance rate"
C        tx10out=MISSING
C        tx50out=MISSING
C        tx90out=MISSING
C        tn10out=MISSING
C        tn50out=MISSING
C        tn90out=MISSING
c       goto 130
      else
        call threshold(tndata,.5,thresan50, flgtn)
        call threshold(tndata,.9,thresan90, flgtn)
      endif

      call threshold(txdata,.1,thresax10, flgtx)
      if(flgtx.eq.1) then
        write(6,*) "TMAX Missing value overflow in exceedance rate"
C        tx10out=MISSING
C        tx50out=MISSING
C        tx90out=MISSING
C        tn10out=MISSING
C        tn90out=MISSING
c       goto 130
      else
        call threshold(txdata,.5,thresax50, flgtx)
        call threshold(txdata,.9,thresax90, flgtx)
      endif

CAndre      ofile=trim(ifile)//"_thresan10"
CAndre      open(22,file=ofile)
      write(53,*) "  year  thresan10  "
        write(53,'(i6,365f8.2)') i+SYEAR-1, thresan10
CAndre      close(22)
						
CAndre      ofile=trim(ifile)//"_thresan90"
CAndre      open(22,file=ofile)
      write(54,*) "  year  thresan90  "
        write(54,'(i6,365f8.2)') i+SYEAR-1, thresan90
CAndre      close(22)
		
CAndre      ofile=trim(ifile)//"_thresax10"
CAndre      open(22,file=ofile)
      write(55,*) "  year  thresax10  "
        write(55,'(i6,365f8.2)') i+SYEAR-1, thresax10
CAndre      close(22)
						
CAndre      ofile=trim(ifile)//"_thresax90"
CAndre      open(22,file=ofile)
      write(56,*) "  year  thresax90  "
        write(56,'(i6,365f8.2)') i+SYEAR-1, thresax90
CAndre      close(22)

      do i=1,BYRS
        txboot=txdata
        tnboot=tndata
        nn=0
        do iter=1,BYRS
          if(iter.ne.i) then
            nn=nn+1
            do day=1,365+2*SS
              if(flgtx.eq.0) txboot(i,day)=txboot(iter,day)
              if(flgtn.eq.0) tnboot(i,day)=tnboot(iter,day)
            enddo
            if(flgtx.eq.0)then
              call threshold(txboot,.9,threstmp,flg)
              do day=1,365
                thresbx90(day,i,nn)=threstmp(day)
              enddo
              call threshold(txboot,.5,threstmp,flg)
              do day=1,365
                thresbx50(day,i,nn)=threstmp(day)
              enddo
              call threshold(txboot,.1,threstmp,flg)
              do day=1,365
                thresbx10(day,i,nn)=threstmp(day)
              enddo
            endif

            if(flgtn.eq.0) then
              call threshold(tnboot,.9,threstmp,flg)
              do day=1,365
                thresbn90(day,i,nn)=threstmp(day)
              enddo
              call threshold(tnboot,.5,threstmp,flg)
              do day=1,365
                thresbn50(day,i,nn)=threstmp(day)
              enddo
              call threshold(tnboot,.1,threstmp,flg)
              do day=1,365
                thresbn10(day,i,nn)=threstmp(day)
              enddo
            endif
          endif
        enddo
      enddo

CAndre      ofile=trim(ifile)//"_thresbn10"
CAndre      open(22,file=ofile)
      write(57,*) "  year  thresbn10  "
      write(57,*) i+SYEAR-1, thresbn10
CAndre      close(22)
						
CAndre	    ofile=trim(ifile)//"_thresbn90"
CAndre      open(22,file=ofile)
      write(58,*) "  year  thresbn90  "
	write(58,*) i+SYEAR-1, thresbn90
CAndre      close(22)
		
CAndre      ofile=trim(ifile)//"_thresbx10"
CAndre      open(22,file=ofile)
      write(59,*) "  year  thresbx10  "
	write(59,*) i+SYEAR-1, thresbx10
CAndre      close(22)
						
CAndre      ofile=trim(ifile)//"_thresbx90"
CAndre      open(22,file=ofile)
      write(60,*) "  year  thresbx90  "
	write(60,*) i+SYEAR-1, thresbx90
CAndre      close(22)

      if(flgtx.eq.0)then
        tx10out=0.
        tx50out=0.
        tx90out=0.
      endif
      if(flgtn.eq.0)then
        tn10out=0.
        tn50out=0.
        tn90out=0.
      endif
      cnt=0
      do i=1,YRS
        year=i+SYEAR-1
        byear=year-BASESYEAR+1
        nn=0
        do month=1,12
          missncnt=0
          missxcnt=0
          if(leapyear(year).eq.1)then
            kth=MONLEAP(month)
          else
            kth=MON(month)
          endif
          do day=1,kth
            if(month.ne.2.or.day.ne.29) nn=nn+1
            cnt=cnt+1
            if(TMAX(cnt).ne.MISSING)then
              if(year.lt.BASESYEAR.or.year.gt.BASEEYEAR) then
                if(TMAX(cnt).gt.thresax90(nn)) tx90out(i,month)=
     &                   tx90out(i,month)+1
                if(TMAX(cnt).gt.thresax50(nn)) tx50out(i,month)=
     &                   tx50out(i,month)+1
                if(TMAX(cnt).lt.thresax10(nn)) tx10out(i,month)=
     &                   tx10out(i,month)+1
              else
                do iter=1,BYRS-1
c                   if(byear.gt.30.or.iter.gt.29) then
c                     print *, i, year,month,day
c                     stop
c                   endif
                  if(TMAX(cnt).gt.thresbx90(nn,byear,iter))then
                    tx90out(i,month)=tx90out(i,month)+1
                  endif
                  if(TMAX(cnt).gt.thresbx50(nn,byear,iter))then
                    tx50out(i,month)=tx50out(i,month)+1
                  endif
                  if(TMAX(cnt).lt.thresbx10(nn,byear,iter))
     &              tx10out(i,month)=tx10out(i,month)+1
                enddo
              endif
            else
              missxcnt=missxcnt+1
            endif
            if(TMIN(cnt).ne.MISSING)then
              if(year.lt.BASESYEAR.or.year.gt.BASEEYEAR) then
                if(TMIN(cnt).gt.thresan90(nn)) tn90out(i,month)=
     &                 tn90out(i,month)+1
                if(TMIN(cnt).gt.thresan50(nn)) tn50out(i,month)=
     &                 tn50out(i,month)+1
                if(TMIN(cnt).lt.thresan10(nn)) tn10out(i,month)=
     &                   tn10out(i,month)+1
              else
                do iter=1,BYRS-1
                  if(TMIN(cnt).gt.thresbn90(nn,byear,iter))
     &              tn90out(i,month)=tn90out(i,month)+1
                  if(TMIN(cnt).gt.thresbn50(nn,byear,iter))
     &              tn50out(i,month)=tn50out(i,month)+1
                  if(TMIN(cnt).lt.thresbn10(nn,byear,iter))
     &              tn10out(i,month)=tn10out(i,month)+1
                enddo
              endif
            else
              missncnt=missncnt+1
            endif
          enddo ! do day=1,kth

c         if(year.ge.BASESYEAR.and.year.le.BASEEYEAR)then
c           print *, year,month,tx10out(i,month),tx90out(i,month),
c    &        tn10out(i,month),tn90out(i,month),missxcnt,missncnt
c         endif

          if(year.ge.BASESYEAR.and.year.le.BASEEYEAR)then
            tn90out(i,month)=tn90out(i,month)/(BYRS-1.)
            tn50out(i,month)=tn50out(i,month)/(BYRS-1.)
            tn10out(i,month)=tn10out(i,month)/(BYRS-1.)
            tx90out(i,month)=tx90out(i,month)/(BYRS-1.)
            tx50out(i,month)=tx50out(i,month)/(BYRS-1.)
            tx10out(i,month)=tx10out(i,month)/(BYRS-1.)
          endif

c         if(year.eq.1952) then
c           print *,year,month,tn10out(i,month),tn90out(i,month)
c         endif

          if(missxcnt.lt.kth.and.flgtx.eq.0)then
            tx90out(i,13)=tx90out(i,13)+tx90out(i,month)
            tx90out(i,month)=tx90out(i,month)*100./(kth-missxcnt)
            tx50out(i,13)=tx50out(i,13)+tx50out(i,month)
            tx50out(i,month)=tx50out(i,month)*100./(kth-missxcnt)
            tx10out(i,13)=tx10out(i,13)+tx10out(i,month)
            tx10out(i,month)=tx10out(i,month)*100./(kth-missxcnt)
          else
            tx90out(i,month)=MISSING
            tx50out(i,month)=MISSING
            tx10out(i,month)=MISSING
          endif
          if(missncnt.lt.kth.and.flgtn.eq.0)then
            tn90out(i,13)=tn90out(i,13)+tn90out(i,month)
            tn90out(i,month)=tn90out(i,month)*100./(kth-missncnt)
            tn50out(i,13)=tn50out(i,13)+tn50out(i,month)
            tn50out(i,month)=tn50out(i,month)*100./(kth-missncnt)
            tn10out(i,13)=tn10out(i,13)+tn10out(i,month)
            tn10out(i,month)=tn10out(i,month)*100./(kth-missncnt)
          else
            tn90out(i,month)=MISSING
            tn50out(i,month)=MISSING
            tn10out(i,month)=MISSING
          endif
        enddo ! do month=1,12
        if(YNASTAT(i,3).eq.1.or.flgtn.eq.1) then
          tn10out(i,13)=MISSING
          tn50out(i,13)=MISSING
          tn90out(i,13)=MISSING
        else
          tn10out(i,13)=tn10out(i,13)*100/365.
          tn50out(i,13)=tn50out(i,13)*100/365.
          tn90out(i,13)=tn90out(i,13)*100/365.
        endif
        if(YNASTAT(i,2).eq.1.or.flgtx.eq.1) then
          tx10out(i,13)=MISSING
          tx50out(i,13)=MISSING
          tx90out(i,13)=MISSING
        else
          tx10out(i,13)=tx10out(i,13)*100/365.
          tx50out(i,13)=tx50out(i,13)*100/365.
          tx90out(i,13)=tx90out(i,13)*100/365.
        endif
      enddo

130   continue
c      if(flgtx.eq.0)then
CAndre      ofile=trim(ifile)//"_TX90p"
CAndre      open(22,file=ofile)
      write(61, *) "  year  jan   feb   mar   apr   may   jun  ",
     &             " jul   aug   sep   oct   nov   dec annual"
      do i=1,YRS
        write(61,'(i6,13f6.1)') i+SYEAR-1,(tx90out(i,j),j=1,13)
      enddo
CAndre      close(22)

CAndre      ofile=trim(ifile)//"_TX50p"
CAndre      open(22,file=ofile)
      write(62, *) "  year  jan   feb   mar   apr   may   jun  ",
     &             " jul   aug   sep   oct   nov   dec annual"
      do i=1,YRS
        write(62,'(i6,13f6.1)') i+SYEAR-1,(tx50out(i,j),j=1,13)
      enddo
CAndre      close(22)

CAndre      ofile=trim(ifile)//"_TX10p"
CAndre      open(22,file=ofile)
      write(63, *) "  year  jan   feb   mar   apr   may   jun  ",
     &             " jul   aug   sep   oct   nov   dec annual"
      do i=1,YRS
        write(63,'(i6,13f6.1)') i+SYEAR-1,(tx10out(i,j),j=1,13)
      enddo
CAndre      close(22)
c      endif

C      if(flgtn.eq.0)then
CAndre      ofile=trim(ifile)//"_TN90p"
CAndre      open(22,file=ofile)
      write(64, *) "  year  jan   feb   mar   apr   may   jun  ",
     &             " jul   aug   sep   oct   nov   dec annual"
      do i=1,YRS
        write(64,'(i6,13f6.1)') i+SYEAR-1,(tn90out(i,j),j=1,13)
      enddo
CAndre      close(22)

CAndre      ofile=trim(ifile)//"_TN50p"
CAndre      open(22,file=ofile)
      write(65, *) "  year  jan   feb   mar   apr   may   jun  ",
     &             " jul   aug   sep   oct   nov   dec annual"
      do i=1,YRS
        write(65,'(i6,13f6.1)') i+SYEAR-1,(tn50out(i,j),j=1,13)
      enddo
CAndre      close(22)

CAndre      ofile=trim(ifile)//"_TN10p"
CAndre      open(22,file=ofile)
      write(66, *) "  year  jan   feb   mar   apr   may   jun  ",
     &             " jul   aug   sep   oct   nov   dec annual"
      do i=1,YRS
        write(66,'(i6,13f6.1)') i+SYEAR-1,(tn10out(i,j),j=1,13)
      enddo
CAndre      close(22)
c      endif
c      write (*,*) "de la soul hear"
      cnt=0
      wsdi=0.
      csdi=0.
c     if(flg.eq.1) then
c       wsdi=MISSING
c       csdi=MISSING
c       goto 140
c     endif

      do i=1,YRS
        cntx=0
        cntn=0
        nn=0
        year=i+SYEAR-1
        do month=1,12
          if(leapyear(year).eq.1)then
            kth=MONLEAP(month)
          else
            kth=MON(month)
          endif
          do day=1,kth
            if(month.ne.2.or.day.ne.29) nn=nn+1
            cnt=cnt+1
            if(TMAX(cnt).gt.thresax90(nn).and.TMAX(cnt).ne.MISSING)then
              cntx=cntx+1
              if(month.eq.12.and.day.eq.31.and.cntx.ge.6)
     &                  wsdi(i)=wsdi(i)+cntx
            elseif(cntx.ge.6)then
              wsdi(i)=wsdi(i)+cntx
              cntx=0
            else
              cntx=0
            endif
            if(TMIN(cnt).lt.thresan10(nn).and.TMIN(cnt).ne.MISSING)then
              cntn=cntn+1
              if(month.eq.12.and.day.eq.31.and.cntn.ge.6)
     &                  csdi(i)=csdi(i)+cntn
            elseif(cntn.ge.6)then
              csdi(i)=csdi(i)+cntn
              cntn=0
            else
              cntn=0
            endif
          enddo  ! day
        enddo    ! month
        if(YNASTAT(i,3).eq.1) csdi(i)=MISSING
        if(YNASTAT(i,2).eq.1) wsdi(i)=MISSING
      enddo      ! year

140   continue
      if(flgtn.eq.0) then
CAndre      ofile=trim(ifile)//"_WSDI"
CAndre      open(22,file=ofile)
      write(67,*) " year     wsdi"
      do i=1,YRS
        write(67,'(i6,f6.1)') i+SYEAR-1,wsdi(i)
      enddo
CAndre      close(22)
      endif

      if(flgtx.eq.0)then
CAndre      ofile=trim(ifile)//"_CSDI"
CAndre      open(22,file=ofile)
      write(68,*) " year     csdi"
      do i=1,YRS
        write(68,'(i6,f6.1)') i+SYEAR-1,csdi(i)
      enddo
CAndre      close(22)
      endif

      end

      subroutine threshold(idata, lev, odata, flg)
      use COMMpreCan
      real idata(BYRS,365+2*SS),odata(365), lev
      integer flg

      real tosort(BYRS*WINSIZE)
      integer nn

      do i=1,365
        nn=0
        do j=1,BYRS
          do k=i,i+2*SS
            if(idata(j,k).ne.MISSING) then
              nn=nn+1
              tosort(nn)=idata(j,k)
            endif
          enddo
        enddo
CAndre        print*,"nn= ",nn
        if (nn.lt.int(BYRS*WINSIZE*.85)) then
           flg=1
	   odata(i)=-99.9
	   else 
              odata(i)=percentile(tosort,nn,lev)
        endif

      enddo

      end