program xspl
c  creates input files for anomaly interpolation from CRU 
c   monthly time series files
c  only years with all months missing are excluded
c  sets values gt defined sigma equal to missing
      parameter (xmiss=-9999)
      
      character*70 infl,mhfmt,mdfmt,splfmt,logfmt
      character name*36
      
      integer data(1701:2000,12),itot(1701:2000,12)
      integer lat1,lon1,lat2,lon2
      real sum(12),n(12),anom(12),xdata(12),var(12),sum2(12),sd(12)
      integer limit(4)
      
c      mhfmt='(i7,i5,i6,i5,a33,2i4)'
      mhfmt='(i7,i6,i7,i5,a36,i4,x,i4)'
      mdfmt='(i4,12i5)'
      splfmt='(3f8.2,12f8.1,12f8.4,,1x,i7,1x,a36,1x,i4)'
      logfmt='(3f8.2,12f8.1,12f8.4,x,i7,x,a36,x,i4)'
c
      open(99,file='log.dat',status='replace')
c
      write(*,*)'Enter start year, end year, start normal period,'
      write(*,*)' end normal period, nyears for normal'
      read(*,*)n1,n2,nor1,nor2,nyrs
      write(*,*)n1,n2,nor1,nor2,nyrs
      write(*,*) 'Enter min and max abs OK values (-999=no constraint)'
      read(*,*) minok,maxok
      write(*,*) 'Enter plus and minus sigma factors for extremes check'
      write(*,*) ' e.g. 5.0 -5.0'
      read(*,*) sigplus,sigminus
      write(*,*) 'Convert to percent anomalies (y=1 / no=0)?'
      read(*,*)iperc
c
      suspect=0.0
5     call openf(1,'Input time series file','old')
      write(*,*)'Enter missing value in time series data file'
      read(*,*)imiss
      write(*,*)imiss
      write(*,*)'Specify window (yes=1 /no=0)?'
      read(*,*)window
      if(window.eq.1)then
       write(*,*)'Enter latmin lonmin latmax lonmax'
       read(*,*)lat1,lon1,lat2,lon2
       lat1=lat1*100
       lat2=lat2*100
       lon1=lon1*100
       lon2=lon2*100
      else
       lat1=-9000
       lon1=-18000
       lat2=9000
       lon2=18000
      endif
      write(*,*)'Reading and writing time series file'
1     read(1,mhfmt,end=99)iwmo,lat,lon,ielv,name,iy1,iy2
      icount=icount+1
      if(500*(icount/500).eq.icount)write(*,*)icount,match,ival
      ncount=0
      do iy=iy1,iy2
        read(1,mdfmt)iyear,(data(iy,im),im=1,12)
        if(iy.ge.nor1.and.iy.le.nor2)ncount=ncount+1
      enddo
      if(lat.lt.lat1.or.lon.lt.lon1.or.lat.gt.lat2.or.lon.gt.lon2.or.
     &    ncount.lt.nyrs)goto 1
      match=match+1
      do im=1,12
       sum(im)=0
       sum2(im)=0
       n(im)=0
      enddo

      do iy=iy1,iy2
       do im=1,12
        if(iy.ge.nor1.and.iy.le.nor2.and.data(iy,im).gt.imiss)then
         sum(im)=sum(im)+real(data(iy,im))/10.0
	 sum2(im)=sum2(im)+(real(data(iy,im))/10.0)**2
         n(im)=n(im)+1
        endif
       enddo
      enddo
      do iy=iy1,iy2
       ipres=0 ; seqstatus=0.0
       do 55 im=1,12
        if(iy.eq.iy1.and.n(im).gt.0)then
	 sum(im)=sum(im)/n(im)
	 sum2(im)=sum2(im)/n(im)
	 diff=sum2(im)-sum(im)**2
	 if(diff.lt.0)diff=0.0
	 sd(im)=sqrt(diff)
        endif
        if(iy.ge.n1.and.iy.le.n2.and.data(iy,im).gt.imiss.and.
     &  	n(im).ge.nyrs)then
         valuestatus=0
         condmin=1 ; condmax=1
         if(minok.eq.-999.or.data(iy,im).ge.minok)condmin=0
         if(maxok.eq.-999.or.data(iy,im).le.maxok)condmax=0
         if(condmin.eq.0.and.condmax.eq.0)then
           anom(im)=real(data(iy,im))/10.0-sum(im)
	   var(im)=1/n(im)
           ipres=1
	   sdt=sdtest(n1,n2,iy,im,data,imiss,nyrs)
           condneg=1 ; condpos=1
           if(anom(im).lt.sigminus*sdt.and.sdt.gt.0)condneg=0
           if(anom(im).gt.sigplus*sdt.and.sdt.gt.0)condpos=0
           if(condneg.eq.0.or.condpos.eq.0)then
             if((anom(im)/sdt).gt.sigplus+1.or.anom(im)/sdt.lt.
     &			sigminus-1)then
               valuestatus=2 
               seqstatus=(seqstatus*2.0)+1.0
               if(seqstatus.gt.2)valuestatus=3
             else
               seqstatus=(seqstatus*2.0)+1.0
               if(seqstatus.gt.2)valuestatus=3
             endif             
           endif
         else
           valuestatus=3
           seqstatus=(seqstatus*2.0)+1.0
         endif
         if(seqstatus.gt.0)seqstatus=seqstatus-0.5
         if(seqstatus.lt.0)seqstatus=0.0
c
	 if(valuestatus.eq.3)then
	  write(99,'(i4,i3,i8,1x,a36)')iy,im,iwmo,name
	  iqy1=iy-10 ; if(iqy1.lt.iy1)iqy1=iy1
	  iqy2=iy+10 ; if(iqy2.gt.iy2)iqy2=iy2
	  do iqy=iqy1,iqy2
	    write(99,'(i4,12i5)')iqy,(data(iqy,iqm),iqm=1,12)
	  enddo
	  write(*,*)'Enter the first mon,yr & last mon,yr to reject:'
	  write(*,*)'  the earliest useful yr is the present=', iy
	  read(*,*)irejmon0,irejyr0,irejmon1,irejyr1
	  do irejyr=irejyr0,irejyr1
	    do irejmon=1,12
	      if(irejyr.eq.irejyr0.and.irejmon.lt.irejmon0)then
	      else if(irejyr.eq.irejyr1.and.irejmon.gt.irejmon1)then
	      else
	        data(irejyr,irejmon)=imiss
	        if(irejyr.eq.iy)then
	          anom(irejmon)=xmiss ; var(irejmon)=-9 ; ipres=0
	        endif
	      endif
	    enddo
	  enddo
	 else if(valuestatus.gt.0)then
	  suspect=suspect+1
	  anom(im)=xmiss
	  var(im)=-9    
	  ipres=0
	 endif
        else
	 anom(im)=xmiss
	 var(im)=-9    
        endif
55     enddo
       if(iperc.eq.1.and.ipres.eq.1)then
	do im=1,12
	 if(sum(im).eq.0)anom(im)=xmiss
	 if(anom(im).ne.xmiss)anom(im)=100*anom(im)/sum(im)
	 if(anom(im).gt.500)anom(im)=500
        enddo
       endif
c 			@@@@ added by TDM to avoid untraceable format errors in write to iy
       do im=1,12
         if(var(im).ne.-9)then
           if(var(im).lt.nyrs)var(im)=1.0/real(nyrs-1)
           if(var(im).gt.nyrs)var(im)=1.0/real(nyrs+1)
         endif
       enddo
c       write(99,'(i7)')iwmo
c       if(iwmo.eq.-965333)write(*,*)real(lat)/100,real(lon)/100,
c     &                   real(ielv)/1000,(anom(im),im=1,12),
c     &			 (var(im),im=1,12),iwmo,name,iy
c       if(ipres.eq.1)write(99,logfmt)real(lat)/100,real(lon)/100,
c     &                   real(ielv)/1000,(anom(im),im=1,12),
c     &			 (var(im),im=1,12),iwmo,name,iy
       if(ipres.eq.1)write(iy,splfmt)real(lat)/100,real(lon)/100,
     &                   real(ielv)/1000,(anom(im),im=1,12),
     &                   (var(im),im=1,12),iwmo,name,iy
c       if(ipres.eq.1)write(99,'(a7)')'written'
      enddo
      write(98,'(i7,a36,12f5.1)')iwmo,name,(sd(im),im=1,12)
      goto 1
99    close(1)
      write(*,*) 'Extract from another time series file (y=1 / n=0)?'
      read(*,*) repeat
      if(repeat.eq.1)goto 5
      if(suspect.gt.0)then
       write(*,*)suspects,' suspect values were found'
       write(*,*)'      see fort.99'
      endif
c
      close(99)
c
      end
      
      subroutine openf(iunit,prompt,oldnew)

      character*(*) prompt,oldnew
      character fname*80,yes*1
      logical fexist
      
1     write(*,*)prompt
      write(*,*)'or enter ''XX'' to quit'
      read(*,'(a80)')fname
      if(fname(1:2).eq.'XX')stop
      write(*,*)fname
      write(*,*)
      do i=1,75
       if(fname(i:i+5).eq.'     ')goto 5
      enddo
5     continue
      inquire(file=fname,exist=fexist)
      if(oldnew.eq.'new')then 
       if(fexist)then 
        write(*,*)'File already exists - open it anyway (y/n)'
        read(*,'(a1)')yes
        write(*,*)
        if(yes.eq.'y')then
	 open(iunit,file=fname,status='old')
        else
         goto 1
        endif
       else
        open(iunit,file=fname,status='new')
       endif
      endif
      if(oldnew.eq.'old')then 
       if(.not.fexist)then 
        write(*,*)'File does not exist - open it anyway (y/n)'
        read(*,'(a1)')yes
        write(*,*)
        if(yes.eq.'y')then
	 open(iunit,file=fname,status='new')
        else
         goto 1
        endif
       else
        open(iunit,file=fname,status='old')
       endif
      endif
      if(oldnew.eq.'unknown')open(iunit,file=fname,status='unknown')
      end

      function sdtest(n1,n2,iy,im,data,imiss,nyrs)
      integer data(1701:2000,12)
      sum=0.0
      sum2=0.0
      n=0
      do jy=n1,n2
       if(jy.ne.iy.and.data(jy,im).ne.imiss)then
	sum=sum+0.1*real(data(jy,im))
	sum2=sum2+(0.1*real(data(jy,im)))**2
	n=n+1
       endif
      enddo
      if(n.gt.nyrs)then
       sum=sum/n
       sum2=sum2/n
       diff=sum2-sum**2
       if(diff.lt.0)diff=0.0
       sdtest=sqrt(diff)
      else
       sdtest=imiss
      endif
      return
      end