      program monthly_adjustments

c     originally used to test many options, this program now (8/95)
c     takes the annual adjustment from output from E-P program
c     reads in monthly data & reference series & difference series,
c     calculates the proper monthly adjustment, and creates a finished,
c     homogeneity adjusted data set.  (actually QC should be rerun on
c     the adjusted data).
c     We have changed how we eliminate stations with too little data
c     so this program also accesses a list of stations to be included
c     in the final adjusted data set (the E-P results don't mention
c     anything if the station passes without adjustments).
c                                          t. peterson

c     current list of programs needed:
c     1.  pe.f reference series creation program both annual and monthly
c     2.  ep.f adjustment program to determine when and how large the
c                   adjustments should be.
c     3.  sta.use.f  to determine which of the annual time series stations
c                   we should use
c     4.  adj.mon.f   to spread the adjustment amounts around the year
c       plus QC before and after, nearest neighbor, and filter canadian ...



c     This program is needed because we discovered that annual reference
c     series are much more reliable than monthly.  So we adjust based on
c     the annual year that was determined and the annual amount that was
c     determined, but we let the monthly reference series determine how
c     the annual is distributed (e.g. more in winter and less in summer).
c     But since the monthly has considerable noise, we smooth the resulting
c     adjustments.





      integer ista(7000),iyr(3,7000)
      integer iuse(9000) !the homogeneous stations to use
      integer kyr(3,100)
      real adj(7000),aadj(100)
      real dif(150,12),dat(150,12)
      real d(2,12)
      real r(12),ddd(12),d4(12)
      real finadj(150,12) !the final adjustment value to add
      real preadj(150,12) ! a preliminary before summed up
      real rnumc(13) !how many missing per year

      zindef=-99.99 !missing
      zin=zindef+1
      id=1
      ng=1
      bmin=9999

	write (*,*) ' '
	write (*,*) '1 = Max. Temperature'
	write (*,*) ' '
	write (*,*) '2 = Min. Temperature'
	write (*,*) ' '
	read (*,*) iChoice
	write (*,*) ' '


	if (iChoice.eq.1) then
          open(unit=1,file='max.raw_month_and_ann_ref')! input file
          open(unit=11,file='tmax.mon.adj') ! output file
	else
          open(unit=1,file='min.raw_month_and_ann_ref') ! input file
          open(unit=11,file='tmin.mon.adj') ! output file
	endif
	  
cccccccccccccccccccccccccccccccccccccccc
c  the stations to use file:  (ep results don't list these)
      if (iChoice.eq.1) then
        open(unit=9,file='max.inv')
      else 
	open(unit=9,file='min.inv')
      endif
      do i=1,99999
         read(9,*,end=50)iuse(i)
      enddo
 50   continue
      nuse=i-1  !the number of stations to use
      write(*,*)'number of stations to use: ',nuse,' last: ',iuse(nuse)
      close(unit=9)

ccccccccccccccccccccccccccccccccccccccccc
c  the results from ep:  years of offsets
      if (iChoice.eq.1) then
	open(unit=3,file='tmax.stat.HII.F')
      else
	open(unit=3,file='tmin.stat.HII.F')
      endif

      do i=1,7000
        read(3,*,end=5)j1,ista(i),ikj,iyr(1,i),adj(i),iyr(2,i),iyr(3,i)
c        write(*,*)i,ista(i)
      enddo
 5    continue
      nadj=i-1  !the number of adjustments
      write(*,*)'last adj: ',nadj,ista(nadj),adj(nadj)
      close(unit=3)


      



cccccccccccccccccccccccccccccccccccc
c     read cand-reference series
 10   continue
      ksta=jsta
      do i=1,150
         do j=1,12
            dif(i,j)=zindef
            dat(i,j)=zindef
         enddo
      enddo

      if(istart.ne.0)then
         iystart=jyr-1850
         dif(jyr-1850,jmon)=dd
         dat(jyr-1850,jmon)=f1
      endif

 100  continue
      if(iend.eq.1)go to 500
      read(1,*,end=101)jsta,jyr,jmon,f1,f2,dd
      if(f1.lt.zin)f1=zindef
      if(f2.lt.zin)f2=zindef
      if(dd.lt.zin)dd=zindef  !make sure all have same missing value
      go to 102
 101  continue
      iend=1
      go to 200
 102  continue
      if(istart.eq.0)then
         iystart=jyr-1850
         istart=1
         ksta=jsta
      endif
      if(jyr.lt.1851)go to 100 !useless before 1850 so skip
      if(ksta.eq.jsta)then
         dif(jyr-1850,jmon)=dd
         dat(jyr-1850,jmon)=f1
         iyend=jyr-1850
         go to 100
      endif

 200  continue !all monthly data for that station read in

      do i=ng,nuse
         if(iuse(i).eq.ksta)then !a station to use
            ng=i
            go to 222
         endif
      enddo
c  if not a station to use then
      ng=1
      go to 10

 222  continue
      


c find the adjustment times & years
      do i=1,100
         aadj(i)=zindef
         do j=1,3
            kyr(j,i)=zindef
         enddo
      enddo

      iloop1=0
      num=0
      do i=id,nadj
         iloop2=0
         if(iloop1.eq.1)iloop2=1
         iloop1=0
         if(ista(i).eq.ksta)then
            num=num+1
            do j=1,3
               kyr(j,num)=iyr(j,i)
            enddo
            aadj(num)=adj(i)
            iloop1=1
         endif
         if(iloop1.eq.0.and.iloop2.eq.1)then !time to leave
            id=i-1
            go to 300
         endif
      enddo
c   if no adjustments for that station then
      id=1
      go to 450 ! to write out the data
 300  continue


      do i=1,150
         do j=1,12
            preadj(i,j)=0
            finadj(i,j)=0
         enddo
      enddo



      do i=1,num !loop through all the adjustment times
         
         j1=kyr(2,i)-1850
         j2=kyr(1,i)-1850 !the year of the discontinuity
         j3=kyr(3,i)-1850 !chronological order


c  set the j's so that if the window is wide, we use a smaller window
c  as defined by Dave in the EP paper.  -for monthly, we must use the 
c  same period as the annual adjustments used
         if(j2-j1.gt.12)j1=j2-12
         if(j3-(j2+1).gt.12)j3=j2+12
         if(j1.lt.1)then
            j1=1
            go to 400
         endif
         if(j2.lt.1)then
         endif
         do k=1,12
            r(k)=0
            do j=1,2
               d(j,k)=0
            enddo
         enddo

c calculate mean of part prior to the discontinuity
            do k=1,12
               r(k)=0
            enddo
         do j=j1,j2
            do k=1,12
               if(dif(j,k).gt.zin)then
                  d(1,k)=d(1,k)+dif(j,k)
                  r(k)=r(k)+1
               endif
            enddo
         enddo
         do k=1,12
            if(r(k).gt.0.5)then
               d(1,k)=d(1,k)/r(k)
            else 
               d(1,k)=zindef
            endif
         enddo
            do k=1,12
               r(k)=0
            enddo

c  calculate the means after the discontinuity
         do j=j2+1,j3
            do k=1,12
               if(dif(j,k).gt.zin)then
                  d(2,k)=d(2,k)+dif(j,k)
                  r(k)=r(k)+1
               endif
            enddo
         enddo
         rnum=0
         do k=1,12
            if(r(k).gt.0.5)then
               d(2,k)=d(2,k)/r(k)
            else 
               d(2,k)=zindef
            endif
         enddo


c now calculate subtract early from late
         ds1=0
         ds2=0
         do k=1,12
            if(d(1,k).gt.zin)ds1=ds1+d(1,k)
            if(d(2,k).gt.zin)ds2=ds2+d(2,k)
         enddo
         do k=1,12
            if(d(1,k).gt.zin.and.d(2,k).gt.zin)then
               ddd(k)=d(2,k)-d(1,k)
            else 
               ddd(k)=zindef
            endif
         enddo

c adjust up to meet same mean:  adding a value 
c  mean of monthly should equal the annual value
         sa=0
         r12=0
         do k=1,12
            if(ddd(k).gt.zin)then
               sa=sa+ddd(k)
               r12=r12+1
            endif
         enddo
         if(r12.lt.0.1)then  !divide by 0 problem
            write(*,*)'divide by zero problem r12'
            aval=0
         else
            aval=aadj(i)-(sa/r12)
         endif
         do k=1,12
            if(ddd(k).gt.zin)ddd(k)=ddd(k)+aval
         enddo



c     9 pt binomial filter
         do k=1,12
            k1=k-4
            k2=k-3
            k3=k-2
            k4=k-1
            k5=k
            k6=k+1
            k7=k+2
            k8=k+3
            k9=k+4
            if(k1.eq.0)k1=12
            if(k1.eq.-1)k1=11
            if(k1.eq.-2)k1=10
            if(k1.eq.-3)k1=9
            if(k2.eq.0)k2=12
            if(k2.eq.-1)k2=11
            if(k2.eq.-2)k2=10
            if(k3.eq.0)k3=12
            if(k3.eq.-1)k3=11
            if(k4.eq.0)k4=12
            if(k6.eq.13)k6=1
            if(k7.eq.13)k7=1
            if(k7.eq.14)k7=2
            if(k8.eq.13)k8=1
            if(k8.eq.14)k8=2
            if(k8.eq.15)k8=3
            if(k9.eq.13)k9=1
            if(k9.eq.14)k9=2
            if(k9.eq.15)k9=3
            if(k9.eq.16)k9=4

            rnum=0
            d4(k)=0
            ddiv=0 ! weighting to divide by
            if(ddd(k1).gt.zin)then
               d4(k)=d4(k)+ddd(k1)*1
               ddiv=ddiv+1
            endif
            if(ddd(k2).gt.zin)then
               d4(k)=d4(k)+ddd(k2)*8
               ddiv=ddiv+8
            endif
            if(ddd(k3).gt.zin)then
               d4(k)=d4(k)+ddd(k3)*28
               ddiv=ddiv+28
            endif
            if(ddd(k4).gt.zin)then
               d4(k)=d4(k)+ddd(k4)*56
               ddiv=ddiv+56
            endif
            if(ddd(k5).gt.zin)then
               d4(k)=d4(k)+ddd(k5)*70
               ddiv=ddiv+70
            endif
            if(ddd(k6).gt.zin)then
               d4(k)=d4(k)+ddd(k6)*56
               ddiv=ddiv+56
            endif
            if(ddd(k7).gt.zin)then
               d4(k)=d4(k)+ddd(k7)*28
               ddiv=ddiv+28
            endif
            if(ddd(k8).gt.zin)then
               d4(k)=d4(k)+ddd(k8)*8
               ddiv=ddiv+8
            endif
            if(ddd(k9).gt.zin)then
               d4(k)=d4(k)+ddd(k9)*1
               ddiv=ddiv+1
            endif
c  temporary fix for problems when dates wrong
            if(ddiv.lt.0.1)then !no data
               d4(k)=0
            else
               d4(k)=d4(k)/ddiv
            endif
         enddo
   

 400     continue

         rnum=0
         do ij=1,12
            if(ddd(ij).lt.zin)rnum=rnum+1
         enddo

         if(rnum.eq.5)then
            write(*,*)
            write(*,*)ksta,j2+1850
            write(*,*)ddd
         endif


         if(rnum.eq.12)write(*,*)'problem with: ',ksta,rnum,j2
         rnumc(rnum+1)=rnumc(rnum+1)+1
         
c protect against monthly reference series not being present
c by earlier counts this is very rare.  Occasionally one month isn't
c there - no problem, just average in without it
c no case of gt 6 missing, but this would catch it.
c if this happened frequently - say 2 or 3 months missing, we might
c need to address it slightly differently, but 3 out of 3,000 we'll
c just average through the missing essentially treating them as if
c they were similar to their neighbor months which is probably true.
c Here is the solution for stations when an adjustment window doesn't
c have adequate months to use:  right now we set that for .gt.5
c so we can have 5 of the 9 months missing and still get a value from
c the other 4:  
c in max&min, 5 happened 3 times (each), 4 none.  But looking at the
c values, they were pretty smooth, which is what the purpose of this
c program is, so it is OK to use 5 for these rare occasions, I guess.

         if(rnum.gt.5)then
            isk=isk+1
            write(*,*)
            write(*,*)'inadequate monthly ref series:::'
            write(*,*)'station: ',ksta,j2+1850,' months miss: ',rnum
            write(*,*)'we will skip station: ',ksta, 'our ',isk,'skip'
            write(*,*)
            do k=1,12
               d4(k)=0
            enddo
            go to 10
         endif


c  record the adjustment values in the discontinuity year
c  adjustment is applied to that year and all previous years
         do k=1,12
            preadj(j2,k)=d4(k)
         enddo

      enddo !from 1,num

c   now sum up all the adjustments

      do i=149,1,-1
         do j=1,12
            finadj(i,j)=finadj(i+1,j)+preadj(i,j)
         enddo
      enddo

c  create final adjusted data
      do i=iystart,iyend
         do j=1,12
            if(dat(i,j).gt.zin)dat(i,j)=dat(i,j)+finadj(i,j)
         enddo
      enddo

 450  continue
c  write out the final data
      do i=iystart,iyend
         write(11,451)ksta,i+1850,(dat(i,j),j=1,12)
      enddo
 451  format(i7,i5,12f7.2)

      go to 10
 500   continue

       do i=1,13
        write(*,*)'num of yrs with ',i-1,' months missing: ',rnumc(i)
       enddo
 
 989   continue
       write(*,*)'DONE'
       end








