1 | |
---|
2 | c PROGRAM Filtrage |
---|
3 | |
---|
4 | c filtrage sur olr JAN-DEC de 1968 |
---|
5 | |
---|
6 | parameter (nb=3416.,period1=10.,period2=30.) |
---|
7 | |
---|
8 | dimension olr(144,25,nb),vb(nb) |
---|
9 | dimension vvb(nb),vvvb(nb),olrf(144,25,nb) |
---|
10 | |
---|
11 | open(1,file='olr_filtre1030.dat' |
---|
12 | *,form='unformatted',access='direct',recl=nb*144*25*4) |
---|
13 | read(1,rec=1) |
---|
14 | *(((olr(i,j,k),i=1,144),j=1,25),k=1,nb) |
---|
15 | close(1) |
---|
16 | |
---|
17 | c l'ordre maximal du filtrage est tel que 2*MOR+1 < nb |
---|
18 | c veuillez changer dans la subroutine FILTRE le parameter MOR |
---|
19 | |
---|
20 | do lon0=1,144 |
---|
21 | do lat0=1,25 |
---|
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 |
---|
32 | c print*,(olrf(1,1,j),j=1,nb) |
---|
33 | |
---|
34 | open(2,file='olr_filtre1030_bis.dat' |
---|
35 | *,form='unformatted',access='direct',recl=144*25*nb*4) |
---|
36 | write(2,rec=1)(((olrf(lon0,lat0,l),lon0=1,144), |
---|
37 | * lat0=1,25),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 | |
---|