1 | !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
2 | ! MODULE optcld : |
---|
3 | ! |
---|
4 | ! regroupe les variables communes au routines/fonctions pour l'extrapolation/interpolation |
---|
5 | ! des proprietes optiques des gouttes. |
---|
6 | ! |
---|
7 | ! Contient les routines : |
---|
8 | ! INITIALISATION - lecture de la look-up table et initialisation des variables communes. |
---|
9 | ! LOCATE - Recherche de valeur dans une table. |
---|
10 | ! INTERPOLEMOI - interpolation d'une valeur a partir d'un jeu de donnees. |
---|
11 | ! EXTRAPOLEMOI - extrapolation des valeurs en dehors de la look-up table. |
---|
12 | !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
13 | MODULE optcld |
---|
14 | IMPLICIT NONE |
---|
15 | REAL,SAVE :: A_ex,A_ab,A_gg, B_ex,B_ab,B_gg |
---|
16 | REAL,SAVE :: fmin_ex,fmin_ab,fmin_gg |
---|
17 | REAL,SAVE :: r0cld |
---|
18 | REAL,SAVE :: lmin,lmax |
---|
19 | INTEGER,SAVE :: npts |
---|
20 | REAL,ALLOCATABLE,SAVE :: tq_ex(:),tq_ab(:),tq_gg(:),tq_wln(:) |
---|
21 | REAL,ALLOCATABLE,SAVE :: ltq_ex(:),ltq_ab(:),ltq_gg(:), & |
---|
22 | ltq_wln(:) |
---|
23 | REAL,SAVE :: frac_c(3),fhvi_c,fhir_c,lseuil_c |
---|
24 | !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
25 | |
---|
26 | CONTAINS |
---|
27 | |
---|
28 | !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
29 | ! SUBROUTINE iniqcld : |
---|
30 | ! |
---|
31 | ! Initialisation des variables commune et lecture de la look-up table. |
---|
32 | ! |
---|
33 | !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
34 | SUBROUTINE iniqcld() |
---|
35 | IMPLICIT NONE |
---|
36 | ! ---------- LOCAL |
---|
37 | INTEGER :: i,iun |
---|
38 | LOGICAL :: ok |
---|
39 | CHARACTER*100 :: tmp |
---|
40 | iun=-1 |
---|
41 | ! rechercher une unite logique libre sur les 200 premieres. |
---|
42 | DO i=1,200 |
---|
43 | INQUIRE (UNIT=i, OPENED=ok) |
---|
44 | IF (.not.ok) THEN |
---|
45 | iun=i |
---|
46 | exit |
---|
47 | ENDIF |
---|
48 | ENDDO |
---|
49 | ! Une belle condition d'arret ==> aucun unite dispo. |
---|
50 | IF (iun.eq.-1) THEN |
---|
51 | PRINT*,"CATASTROPHE !" |
---|
52 | PRINT*,"Impossible de trouver une unite logique libre", & |
---|
53 | " sur les 200 premieres." |
---|
54 | PRINT*,"Je ne peux pas lire optcld.table." |
---|
55 | STOP "Je m'arrete (et salement en plus)..." |
---|
56 | ENDIF |
---|
57 | ! Lecture de la table !!! |
---|
58 | OPEN(iun,file='optcld.table') |
---|
59 | DO i=1,4 |
---|
60 | READ(iun,*) tmp |
---|
61 | ENDDO |
---|
62 | READ(iun,*) npts |
---|
63 | ! Petite pause dans la lecture pour allouer les tableaux |
---|
64 | ALLOCATE(tq_ex(npts)) |
---|
65 | ALLOCATE(tq_ab(npts)) |
---|
66 | ALLOCATE(tq_gg(npts)) |
---|
67 | ALLOCATE(tq_wln(npts)) |
---|
68 | ALLOCATE(ltq_ex(npts)) |
---|
69 | ALLOCATE(ltq_ab(npts)) |
---|
70 | ALLOCATE(ltq_gg(npts)) |
---|
71 | ALLOCATE(ltq_wln(npts)) |
---|
72 | ! Reprise de la |
---|
73 | READ(iun,*) tmp |
---|
74 | READ(iun,'(ES14.7)') r0cld |
---|
75 | READ(iun,*) tmp |
---|
76 | READ(iun,'(2(ES14.7,2X))') lmin,lmax |
---|
77 | DO i=1,3 |
---|
78 | READ(iun,*) tmp |
---|
79 | ENDDO |
---|
80 | READ(iun,'(2(ES14.7,2X))') A_ex,B_ex |
---|
81 | READ(iun,*) tmp |
---|
82 | READ(iun,'(2(ES14.7,2X))') A_ab,B_ab |
---|
83 | READ(iun,*) tmp |
---|
84 | READ(iun,'(2(ES14.7,2X))') A_gg,B_gg |
---|
85 | DO i=1,3 |
---|
86 | READ(iun,*) tmp |
---|
87 | ENDDO |
---|
88 | READ(iun,'(ES14.7)') fmin_ex |
---|
89 | READ(iun,*) tmp |
---|
90 | READ(iun,'(ES14.7)') fmin_ab |
---|
91 | READ(iun,*) tmp |
---|
92 | READ(iun,'(ES14.7)') fmin_gg |
---|
93 | DO i=1,3 |
---|
94 | READ(iun,*) tmp |
---|
95 | ENDDO |
---|
96 | READ(iun,*) (frac_c(i),i=1,3) |
---|
97 | DO i=1,2 |
---|
98 | READ(iun,*) tmp |
---|
99 | ENDDO |
---|
100 | READ(iun,*) fhvi_c,fhir_c,lseuil_c |
---|
101 | READ(iun,*) tmp |
---|
102 | |
---|
103 | DO i=1,npts |
---|
104 | READ(iun,'(4(ES23.15,1X))') & |
---|
105 | tq_wln(i),tq_ex(i),tq_ab(i),tq_gg(i) |
---|
106 | ! ------------ on passe tout en log pour les interpolations |
---|
107 | ltq_wln(i) = alog(tq_wln(i)) |
---|
108 | ltq_ex(i) = alog(tq_ex(i)) |
---|
109 | ltq_ab(i) = alog(tq_ab(i)) |
---|
110 | ltq_gg(i) = alog(tq_gg(i)) |
---|
111 | ENDDO |
---|
112 | CLOSE(iun) |
---|
113 | |
---|
114 | WRITE(*,*) & |
---|
115 | "LECTURE LOOK-UP optcld.table... TERMINEE :)" |
---|
116 | |
---|
117 | END SUBROUTINE iniqcld |
---|
118 | |
---|
119 | !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
120 | ! SUBROUTINE locate(xx,n,x,j) : |
---|
121 | ! |
---|
122 | ! Recherche de l'indice j de la valeur la plus proche par defaut de x dans le tableau xx de |
---|
123 | ! dimension n |
---|
124 | ! |
---|
125 | ! ARGUMENTS D'ENTREE : |
---|
126 | ! xx : tableau dans lequel rechercher l'indice. |
---|
127 | ! x : valeur recherchee |
---|
128 | ! n : dimension de xx |
---|
129 | ! |
---|
130 | ! ARGUMENT DE SORTIE : |
---|
131 | ! j : indice de la valeur la plus proche PAR DEFAUT de x |
---|
132 | !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
133 | SUBROUTINE locate(xx,n,x,j) |
---|
134 | IMPLICIT NONE |
---|
135 | ! ---------- INPUT |
---|
136 | INTEGER,INTENT(in) :: n |
---|
137 | REAL ,INTENT(in) :: x,xx(n) |
---|
138 | INTEGER,INTENT(out) :: j |
---|
139 | ! ---------- LOCAL |
---|
140 | INTEGER jl,jm,ju |
---|
141 | jl=0 |
---|
142 | ju=n+1 |
---|
143 | DO WHILE (ju-jl.gt.1) |
---|
144 | jm=(ju+jl)/2 |
---|
145 | IF (jm.eq.0) STOP "ALERTE jm=0 !!" |
---|
146 | IF((xx(n).ge.xx(1)).eqv.(x.ge.xx(jm))) THEN |
---|
147 | jl=jm |
---|
148 | ELSE |
---|
149 | ju=jm |
---|
150 | ENDIF |
---|
151 | ENDDO |
---|
152 | IF (x.eq.xx(1))THEN |
---|
153 | j=1 |
---|
154 | ELSE IF (x.eq.xx(n)) THEN |
---|
155 | j=n-1 |
---|
156 | ELSE |
---|
157 | j=jl |
---|
158 | ENDIF |
---|
159 | RETURN |
---|
160 | END SUBROUTINE locate |
---|
161 | |
---|
162 | !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
163 | ! SUBROUTINE interpolemoi(ii,x,xin,yin,npts,yout,ver,loga) : |
---|
164 | ! |
---|
165 | ! Interpolation d'une valeur dans un jeu de donnees xin(npts),yin(npts). |
---|
166 | ! |
---|
167 | ! ARGUMENTS D'ENTREE : |
---|
168 | ! ii : ind dans xin de la valeur la plus proche de x par defaut (locate est ton ami) |
---|
169 | ! x : abscisse de la valeur a interpoler. |
---|
170 | ! xin,yin : jeu de valeurs pour l'interpolation. |
---|
171 | ! npts : nombre de points de xin,yin. |
---|
172 | ! ver : type d'interpolation (0 = lineaire / 1 = quadratique) |
---|
173 | ! loga : interpolation en espace log |
---|
174 | ! |
---|
175 | ! ARGUMENT DE SORTIE : |
---|
176 | ! yout : valeur interpolee en x. ! |
---|
177 | ! |
---|
178 | ! NOTES : |
---|
179 | ! - Si loga est utilisee alors xin et yin doivent alors representer les logarithmes du |
---|
180 | ! jeu de donnees. |
---|
181 | ! - Quelque soit la valeur de loga, yout est la valeur interpolee dans l'espace normal |
---|
182 | !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
183 | SUBROUTINE interpolemoi(ii,x,xin,yin,npts,yout,ver,loga) |
---|
184 | IMPLICIT NONE |
---|
185 | ! ---------- INPUT |
---|
186 | INTEGER,INTENT(in) :: npts |
---|
187 | INTEGER,INTENT(inout) :: ii,ver ! ces 2 variables sont suceptible de changer |
---|
188 | LOGICAL,INTENT(in) :: loga |
---|
189 | REAL ,INTENT(in) :: x,xin(npts),yin(npts) |
---|
190 | ! ---------- OUTPUT |
---|
191 | REAL ,INTENT(out) :: yout |
---|
192 | ! ---------- LOCAL |
---|
193 | REAL :: myx,ytmp,denom |
---|
194 | |
---|
195 | ! ---------- on pourrait ameliorer la condition ici et tester le cas |
---|
196 | ! x plus proche de x(i+1) et utiliser dans ce cas l'interpolation |
---|
197 | ! quadratique... mais est ce vraiment nécessaire ??? |
---|
198 | ! Si on est sur les 2 premiers ou les 2 derniers points : |
---|
199 | ! ===> interpolation lineaire |
---|
200 | IF (ii.eq.1.or.ii.eq.npts-1) ver = 0 |
---|
201 | |
---|
202 | ! Interpolation lineaire |
---|
203 | IF (ver.eq.0) THEN |
---|
204 | myx = x |
---|
205 | IF (loga) myx=alog(myx) |
---|
206 | denom = (xin(ii+1)-xin(ii)) |
---|
207 | ytmp=((yin(ii+1)-yin(ii))*(myx-xin(ii)))/denom+yin(ii) |
---|
208 | IF (loga) THEN |
---|
209 | yout=exp(ytmp) |
---|
210 | ELSE |
---|
211 | yout=ytmp |
---|
212 | ENDIF |
---|
213 | ELSE |
---|
214 | ! ------------ Recherche de l'indice le plus proche de la valeur. |
---|
215 | ! Permet de choisir si l'on interpole avec : |
---|
216 | ! i-1;i;i+1 OU i;i+1;i+2 |
---|
217 | IF (x-xin(ii).gt.xin(ii+1)-x) ii = ii+1 |
---|
218 | myx=x |
---|
219 | IF (loga) myx=alog(myx) |
---|
220 | ytmp = (myx-xin(ii))*(myx-xin(ii+1)) / & |
---|
221 | ((xin(ii-1)-xin(ii)) * & |
---|
222 | (xin(ii-1)-xin(ii+1))) * & |
---|
223 | yin(ii-1) & |
---|
224 | + & |
---|
225 | (myx-xin(ii-1))*(myx-xin(ii+1)) / & |
---|
226 | ((xin(ii)-xin(ii-1)) * & |
---|
227 | (xin(ii)-xin(ii+1))) * & |
---|
228 | yin(ii) & |
---|
229 | + & |
---|
230 | (myx-xin(ii-1))*(myx-xin(ii)) / & |
---|
231 | ((xin(ii+1)-xin(ii-1)) * & |
---|
232 | (xin(ii+1)-xin(ii))) * & |
---|
233 | yin(ii) |
---|
234 | IF (loga) THEN |
---|
235 | yout=exp(ytmp) |
---|
236 | ELSE |
---|
237 | yout=ytmp |
---|
238 | ENDIF |
---|
239 | ENDIF |
---|
240 | RETURN |
---|
241 | END SUBROUTINE interpolemoi |
---|
242 | |
---|
243 | !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
244 | ! SUBROUTINE extrapolemoi(x,A,B,yout,loga) : |
---|
245 | ! |
---|
246 | ! Extrapolation lineaire d'une valeur a partir de coeff A et B issus d'un jeu de donnees |
---|
247 | ! xin(npts), yin(npts). [ f(x) = A*x + B OU exp(A*log(x)+B) ] |
---|
248 | ! |
---|
249 | ! ARGUMENTS D'ENTREE : |
---|
250 | ! x : abscisse de la valeur a extrapoler. |
---|
251 | ! A : coeff directeur de la droite. |
---|
252 | ! B : ordonnee a l'origine de la droite. |
---|
253 | ! loga : interpolation en espace log |
---|
254 | ! |
---|
255 | ! ARGUMENT DE SORTIE : |
---|
256 | ! yout : valeur extrapolee en x. |
---|
257 | ! |
---|
258 | ! NOTES : |
---|
259 | ! - Si loga est utilisee alors A et B n'ont pas la meme signification (voir forme f(x)) |
---|
260 | ! - Quelque soit la valeur de loga, yout est la valeur extrapolee dans l'espace normal |
---|
261 | !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
262 | SUBROUTINE extrapolemoi(x,A,B,yout,loga) |
---|
263 | IMPLICIT NONE |
---|
264 | ! ---------- INPUT |
---|
265 | REAL ,INTENT(in) :: x,A,B |
---|
266 | LOGICAL,INTENT(in) :: loga |
---|
267 | ! ---------- OUTPUT |
---|
268 | REAL ,INTENT(out) :: yout |
---|
269 | |
---|
270 | IF (loga) THEN |
---|
271 | yout = exp(A*alog(x)+B) |
---|
272 | ELSE |
---|
273 | yout = A*x+B |
---|
274 | ENDIF |
---|
275 | |
---|
276 | END SUBROUTINE extrapolemoi |
---|
277 | |
---|
278 | END MODULE optcld |
---|