c PROGRAM Filtrage

c   filtrage sur olr JAN-DEC de 1968

      parameter (nb=3538.,period1=10.,period2=30.)

       dimension olr(24,15,nb),vb(nb)
       dimension vvb(nb),vvvb(nb),olrf(24,15,nb)

      open(1,file='olr_ctl.dat' 
     *,form='unformatted',access='direct',recl=nb*24*15*4)
      read(1,rec=1) 
     *(((olr(i,j,k),i=1,24),j=1,15),k=1,nb)
      close(1)

c  l'ordre maximal du filtrage est tel que 2*MOR+1 < nb
c  veuillez changer dans la subroutine FILTRE le parameter MOR

      do lon0=1,24
      do lat0=1,15
      do i=1,nb
      vb(i)=olr(lon0,lat0,i)
      enddo
      call filtre(vb,vvb,period1,nb)
      call filtre(vb,vvvb,period2,nb)
      do i=1,nb
      olrf(lon0,lat0,i)=vvb(i)-vvvb(i)
      enddo
      enddo
      enddo
c      print*,(olrf(1,1,j),j=1,nb)

        open(2,file='olrf10-30_clim_sst.80.dat' 
     *,form='unformatted',access='direct',recl=24*15*nb*4)
        write(2,rec=1)(((olrf(lon0,lat0,l),lon0=1,24), 
     *                lat0=1,15),l=1,nb)
      close(2)

        END

      SUBROUTINE FILTRE(F,F1,PERIO,N)
      PARAMETER(KOR=4,JOR=4,MOR=50)
      DIMENSION F(N),W(-MOR:mor),G(-MOR:mor),F1(N)
      PI=ACOS(-1.)
      FC=1./PERIO
      CALL KISER(G,MOR)
      DO 1 I=-MOR,MOR
      IF (I.EQ.0) THEN
      W(I)=2.*FC
      ELSE
      W(I)=SIN(2.*PI*FC*FLOAT(I))/(PI*FLOAT(I))*G(I)
      ENDIF
   1  CONTINUE

      DO 2 I=1,N
      F1(I)=0.
      AT=0.

      IF(I.LE.KOR) THEN
      L1=-KOR
      L2=I-1
      ENDIF

      IF((I.GE.KOR+1).AND.(I.LE.MOR)) THEN
      L1=-I+1
      L2=I-1
      ENDIF

      IF((I.GE.MOR+1).AND.(I.LE.N-MOR)) THEN
      L1=-MOR
      L2=MOR
      ENDIF

      IF((I.GE.N-MOR+1).AND.(I.LE.N-JOR)) THEN
      L1=-N+I
      L2=N-I
      ENDIF

      IF(I.GE.N-JOR+1) THEN
      L1=-N+I
      L2=JOR
      ENDIF

      DO 3 K=L1,L2
      F1(I)=F1(I)+W(K)*F(I-K)
      AT=AT+W(K)
   3  CONTINUE
      F1(I)=F1(I)/AT
   2  CONTINUE

      RETURN
      END

      SUBROUTINE KISER(W,MOR)
      PARAMETER(LOR=100)
      DIMENSION W(-MOR:MOR),CO(-LOR:LOR)

      A=30.

      IF(A.LE.21.) THEN
      ALPHA=0.
      ENDIF

      IF((A.LT.50.).AND.(A.GT.21.)) THEN
      ALPHA=0.5842*(A-21.)**0.4+0.07886*(A-21.)
      ENDIF

      IF(A.GE.50.) THEN
      ALPHA=0.1102*(A-8.7)
      ENDIF

      DO 2 I=-MOR,MOR
      CO(I)=ALPHA*SQRT(1.-(FLOAT(I)/FLOAT(MOR))**2)

      N=0
      AS=1.
      AU=1.
      AS1=1.
      AU1=1.

      DO 10 K=1,200
      N=N+1

      AU=AU*((CO(I)/2.)/FLOAT(N))**2
      AS=AS+AU

      AU1=AU1*((ALPHA/2.)/FLOAT(N))**2
      AS1=AS1+AU1
  10  CONTINUE
      W(I)=AS/AS1
   2  CONTINUE
      RETURN
      END