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
C NOW WITH ADDED NORMALS !!!!
      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)'
      nrmfmt='(i7,i6,i7,i5,a36,i4,x,i4,x,i2,x,i3)'
      mdfmt='(i4,12i5)'
      splfmt='(3f8.2,12f8.1,12f8.4,x,i7,x,a20,x,i4)'
      logfmt='(3f8.2,12f8.1,12f8.4,x,i7,x,a20,x,i4)'
      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 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 ITERATION for each .cts file begins here
      suspect=0.0
5     call openf(1,'Input .cts 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 .cts .nrm and writing fort.YYYY'
C PROCESS NORMALS FILE
C ...
C READ .cts STATION HEADER
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
C READ .cts STATION DATA
      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
C CALC NORMAL FOR STATION
      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
C PROCESS time-series data for station
      do iy=iy1,iy2
       ipres=0
       do 55 im=1,12
C CALC NORMAL and STDEV
        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
C CALC ANOMALY
        if(iy.ge.n1.and.iy.le.n2.and.data(iy,im).gt.imiss.and.
     &  n(im).ge.nyrs)then
         anom(im)=real(data(iy,im))/10.0-sum(im)
	 var(im)=1/n(im)
         ipres=1
C PERFORM STDEV TEST
	 sdt=sdtest(n1,n2,iy,im,data,imiss,nyrs)
	 if(anom(im).lt.sigminus*sdt.and.sdt.gt.0)then
	  suspect=suspect+1
	  write(99,'(i4,i3,i8,x,a36,3f10.2)')iy,im,iwmo,name,anom(im),
     &     sum(im),sdt
	  anom(im)=xmiss
	  var(im)=-9    
	  ipres=0
         endif
	 if(anom(im).gt.sigplus*sdt.and.sdt.gt.0)then
	  suspect=suspect+1
	  write(99,'(i4,i3,i8,1x,a36,3f10.2)')iy,im,iwmo,name,anom(im),
     &     sum(im),sdt
	  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
       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(2:21),iy
       if(ipres.eq.1)write(99,'(a6)')'dumped'
      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(*,*)suspect,' suspect values were found'
       write(*,*)'      see fort.99'
      endif
      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