source: trunk/LMDZ.TITAN.old/libf/phytitan/optcld.F90 @ 3557

Last change on this file since 3557 was 808, checked in by slebonnois, 12 years ago

SL: Many changes for VENUS (related to newstart) and TITAN (related to clouds). Please read DOC/chantiers/commit_importants.log (cf v808).

File size: 11.0 KB
Line 
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
Note: See TracBrowser for help on using the repository browser.