source: trunk/src/mode_sahelien/progfiltrage_10_25.f @ 37

Last change on this file since 37 was 10, checked in by pinsard, 15 years ago

first commit with original work of Sebastien Gervois

File size: 2.5 KB
Line 
1
2c PROGRAM Filtrage
3
4c   filtrage sur olr JAN-DEC de 1968
5
6      parameter (nb=10500.,period1=10.,period2=25.)
7
8       dimension olr(25,17,nb),vb(nb)
9       dimension vvb(nb),vvvb(nb),olrf(25,17,nb)
10
11      open(1,file='olrint.80.dat' 
12     *,form='unformatted',access='direct',recl=nb*25*17*4)
13      read(1,rec=1) 
14     *(((olr(i,j,k),i=1,25),j=1,17),k=1,nb)
15      close(1)
16
17c  l'ordre maximal du filtrage est tel que 2*MOR+1 < nb
18c  veuillez changer dans la subroutine FILTRE le parameter MOR
19
20      do lon0=1,25
21      do lat0=1,17
22      do i=1,nb
23      vb(i)=olr(lon0,lat0,i)
24      enddo
25      call filtre(vb,vvb,period1,nb)
26      call filtre(vb,vvvb,period2,nb)
27      do i=1,nb
28      olrf(lon0,lat0,i)=vvb(i)-vvvb(i)
29      enddo
30      enddo
31      enddo
32c      print*,(olrf(1,1,j),j=1,nb)
33
34        open(2,file='olr_1025.dat'
35     *,form='unformatted',access='direct',recl=25*17*nb*4)
36        write(2,rec=1)(((olrf(lon0,lat0,l),lon0=1,25), 
37     *                lat0=1,17),l=1,nb)
38      close(2)
39
40        END
41
42      SUBROUTINE FILTRE(F,F1,PERIO,N)
43      PARAMETER(KOR=4,JOR=4,MOR=50)
44      DIMENSION F(N),W(-MOR:mor),G(-MOR:mor),F1(N)
45      PI=ACOS(-1.)
46      FC=1./PERIO
47      CALL KISER(G,MOR)
48      DO 1 I=-MOR,MOR
49      IF (I.EQ.0) THEN
50      W(I)=2.*FC
51      ELSE
52      W(I)=SIN(2.*PI*FC*FLOAT(I))/(PI*FLOAT(I))*G(I)
53      ENDIF
54   1  CONTINUE
55
56      DO 2 I=1,N
57      F1(I)=0.
58      AT=0.
59
60      IF(I.LE.KOR) THEN
61      L1=-KOR
62      L2=I-1
63      ENDIF
64
65      IF((I.GE.KOR+1).AND.(I.LE.MOR)) THEN
66      L1=-I+1
67      L2=I-1
68      ENDIF
69
70      IF((I.GE.MOR+1).AND.(I.LE.N-MOR)) THEN
71      L1=-MOR
72      L2=MOR
73      ENDIF
74
75      IF((I.GE.N-MOR+1).AND.(I.LE.N-JOR)) THEN
76      L1=-N+I
77      L2=N-I
78      ENDIF
79
80      IF(I.GE.N-JOR+1) THEN
81      L1=-N+I
82      L2=JOR
83      ENDIF
84
85      DO 3 K=L1,L2
86      F1(I)=F1(I)+W(K)*F(I-K)
87      AT=AT+W(K)
88   3  CONTINUE
89      F1(I)=F1(I)/AT
90   2  CONTINUE
91
92      RETURN
93      END
94
95      SUBROUTINE KISER(W,MOR)
96      PARAMETER(LOR=100)
97      DIMENSION W(-MOR:MOR),CO(-LOR:LOR)
98
99      A=30.
100
101      IF(A.LE.21.) THEN
102      ALPHA=0.
103      ENDIF
104
105      IF((A.LT.50.).AND.(A.GT.21.)) THEN
106      ALPHA=0.5842*(A-21.)**0.4+0.07886*(A-21.)
107      ENDIF
108
109      IF(A.GE.50.) THEN
110      ALPHA=0.1102*(A-8.7)
111      ENDIF
112
113      DO 2 I=-MOR,MOR
114      CO(I)=ALPHA*SQRT(1.-(FLOAT(I)/FLOAT(MOR))**2)
115
116      N=0
117      AS=1.
118      AU=1.
119      AS1=1.
120      AU1=1.
121
122      DO 10 K=1,200
123      N=N+1
124
125      AU=AU*((CO(I)/2.)/FLOAT(N))**2
126      AS=AS+AU
127
128      AU1=AU1*((ALPHA/2.)/FLOAT(N))**2
129      AS1=AS1+AU1
130  10  CONTINUE
131      W(I)=AS/AS1
132   2  CONTINUE
133      RETURN
134      END
135
136
Note: See TracBrowser for help on using the repository browser.