#! /bin/sh
#
#    Last change: Norel Rimbu 2004
#
if [ $# -lt 4 ]
then
  echo " "
  echo "           *****     COMPOZCONT     *****      "
  echo " "
  echo "   The program calculates continuous composite maps"
  echo " "
  echo " Syntax: compozcont <file1> <file2> <file3> <win_length>"
  echo " "
  echo "    <file1>      - input series (one record)"
  echo "    <file2>      - input field"
  echo "    <file3>      - output file with composites"
  echo "    <win_length> - window length"
  echo " "
  exit 0
else
 echo " "
 echo "          *****     COMPOZCONT     *****      "
 echo " "
 echo "file1="$1
 echo "file2="$2
 echo "file3="$3
 echo "Window length="$4
 cp $1 fort.1
 cp $2 fort.2
fi
cat > compozcont.f << EOF
C
      program compozcont
      parameter(undef=0.9e+10)
      real  x(1000),z(10000)
      real  cpp(10000),cpn(10000),cpd(10000)
      integer nt, ng, LF
      integer ia(4),ib(4)


C LENGTH OF WINDOW 
      read(*,'(i10)')LF

C START THE ALGORITHM

      open(1,form='unformatted') 
      open(2,form='unformatted') 
      open(3,form='unformatted')

   55 format(a64)

C  READ SERIES

          
      read(1)(ia(l),l=1,4)
      n=ia(4)
      read(1)(x(l),l=1,n)
              
C CALCULATE STANDARD DEV OF THE FIRST INDEX

      nf=nt-LF

      do  888 nn=1,nf
        s1=0    
        a=0.0
        na=LF+nn                   

        do 7 l=nn,na
          a=a+1.0
          s1=s1+x(l)
    7   continue
        xm=s1/a          
        s1=0
        a=0.0
        do 8 l=nn,na
          a=a+1.0
          s1=s1+(x(l)-xm)*(x(l)-xm)
    8   continue
        xd=sqrt(s1/a) 
        xpp=xm+0.75*xd
        xpn=xm-0.75*xd  
        kkp=0
        kkn=0

        do 100 l=1,ng
           cpd(l)=0.0
           cpp(l)=0.0
           cpn(l)=0.0
  100   continue
                  
        do 33 k=nn,na             
              
C  READ MAPS


          read(2)(ib(l),l=1,4)
          ng=ib(4)
          read(2)(z(l),l=1,ng)
          if( (x(k).ge.xpp)) then
            kkp=kkp+1 
            do 77 l=1,ng
              cpp(l)=cpp(l)+z(l)
   77       continue
          endif
                    
          if( (x(k).le.xpn)) then
            kkn=kkn+1  
            do 777 l=1,ng

              cpn(l)=cpn(l)+z(l)
  777       continue
          endif
                    
   33   continue                  

        do 88 l=1,ng
          a=kkp+0.001
          b=kkn+0.001
          aa=cpp(l)/a
          cpp(l)=aa
          bb=cpn(l)/b
          cpn(l)=bb 
          if((abs(aa).lt.100000.0).and.
     *(abs(bb).lt.100000.0)) then                  
            cpd(l)=cpp(l)-cpn(l)
          else
            cpd(l)=undef
          endif
   88   continue

        do 567 l=1,ng
          aa=cpp(l)
          bb=cpn(l)
          cc=cpd(l)
          if(abs(aa).ge.100000.0) then
            cpp(l)=undef
          endif
          if(abs(bb).ge.100000.0) then
            cpn(l)=undef
          endif
          if(abs(cc).ge.100000.0) then
            cpd(l)=undef
          endif
  567   continue


C start to write 


        write(3)(ib(l),l=1,4)
        write(3)(cpd(l),l=1,ng)

  888 continue

      close(1)
      close(2)
      close(3)
      print *,'REGULAR END!'
      
      end
EOF
f77 compozcont.f -o compozcont.x
compozcont.x $1 $2 $3 << M
$4
M
cp fort.3 $3
rm fort.1 fort.2 fort.3
rm compozcont.x compozcont.f
exit