source: trunk/LMDZ.GENERIC/libf/phystd/aerave_new.F @ 1477

Last change on this file since 1477 was 135, checked in by aslmd, 14 years ago

CHANGEMENT ARBORESCENCE ETAPE 2 -- NON COMPLET

File size: 6.6 KB
Line 
1      SUBROUTINE aerave_new ( ndata,
2     & longdata,epdata,omegdata,gdata,         
3     &            longref,epref,temp,nir,longir
4     &            ,epir,omegir,gir,qref,omegaref        )
5
6
7      IMPLICIT NONE
8c.......................................................................
9c
10c R.Fournier 02/1996
11c (modif F.Forget 02/1996)
12c le spectre est decoupe en "nir" bandes et cette routine calcule
13c les donnees radiatives moyenne sur chaque bande : l'optimisation
14c est faite pour une temperature au sol "temp" et une epaisseur
15c optique de l'atmosphere "epref" a la longueur d'onde "longref"
16c
17c dans la version actuelle, les ponderations sont independantes de
18c l'epaisseur optique : c'est a dire que "omegir", "gir"
19c et "epir/epre" sont independants de "epref".
20c en effet les ponderations sont choisies pour une solution exacte
21c en couche mince et milieu isotherme.
22c
23c entree
24c
25c    ndata : taille des champs data
26c    longdata,epdata,omegdata,gdata : proprietes radiative de l'aerosol
27c                  (longdata longueur d'onde en METRES)
28c  * longref : longueur d'onde a laquelle l'epaisseur optique
29c              est connue
30c  * epref : epaisseur optique a longref
31c  * temp : temperature choisie pour la ponderation (Planck)
32c  * nir : nombre d'intervals dans la discretisation spectrale
33c           du GCM
34c  * longir : longueurs d'onde definissant ces intervals
35c
36c sortie
37c
38c  * epir : epaisseur optique moyenne pour chaque interval
39c  * omegir : "scattering albedo" moyen pour chaque interval
40c  * gir : "assymetry factor" moyen pour chaque interval
41c  * qref : extinction coefficient at reference wavelength
42c  * omegaref : single scat. albedo at reference wavelength
43c
44c.......................................................................
45c
46      REAL longref
47      REAL epref
48      REAL temp
49      INTEGER nir
50      REAL*8 longir(nir+1)
51      REAL epir(nir)
52      REAL omegir(nir)
53      REAL gir(nir)
54c
55c.......................................................................
56c
57      INTEGER iir,nirmx
58      PARAMETER (nirmx=100)
59      INTEGER idata,ndata
60c
61c.......................................................................
62c
63      REAL emit
64      REAL totalemit(nirmx)
65      REAL longdata(ndata),epdata(ndata)
66     &    ,omegdata(ndata),gdata(ndata)
67      REAL qextcorrdata(ndata)
68      INTEGER ibande,nbande
69      PARAMETER (nbande=1000)
70      REAL long,deltalong
71      INTEGER ilong
72      INTEGER i1,i2
73      REAL c1,c2
74      REAL factep,qextcorr,omeg,g
75      REAL qref,omegaref
76c
77c.......................................................................
78c
79      DOUBLE PRECISION tmp1
80      REAL tmp2,tmp3
81c
82c
83      long=longref
84
85
86      !if(nir.eq.27)then
87      !print*,'long',long
88      !print*,'longdata',longdata
89      !print*,'epdata',epdata
90      !print*,'omegdata',omegdata
91      !print*,'gdata',gdata
92      !print*,'data looks aok!'
93
94      !print*,'ndata=',ndata
95      !print*,'longdata',shape(longdata)
96      !print*,'epdata',shape(epdata)
97      !print*,'omegdata',shape(omegdata)
98      !print*,'gdata',shape(gdata)
99      ! print*,'longref',longref
100      !print*,'epref',epref
101      !print*,'temp',temp
102      !print*,'nir',nir
103      !print*,'longir',longir
104      !print*,'epir',epir
105      !print*,'omegir',gir
106      !print*,'qref',qref
107      !print*,'longir=',longir
108      !stop
109      !endif
110
111
112c********************************************************
113c interpolation
114      ilong=1
115      DO idata=2,ndata
116        IF (long.gt.longdata(idata)) ilong=idata
117      ENDDO
118      i1=ilong
119      i2=ilong+1
120      IF (i2.gt.ndata) i2=ndata
121      IF (long.lt.longdata(1)) i2=1
122      IF (i1.eq.i2) THEN
123        c1=1.E+0
124        c2=0.E+0
125      ELSE
126        c1=(longdata(i2)-long) / (longdata(i2)-longdata(i1))
127        c2=(longdata(i1)-long) / (longdata(i1)-longdata(i2))
128      ENDIF
129c********************************************************
130c
131      qref=c1*epdata(i1)+c2*epdata(i2)
132      omegaref=c1*omegdata(i1)+c2*omegdata(i2)
133      factep=qref/epref
134      DO idata=1,ndata
135        qextcorrdata(idata)=epdata(idata)/factep
136      ENDDO
137c
138c.......................................................................
139c
140      DO iir=1,nir
141c
142c.......................................................................
143c
144        deltalong=(longir(iir+1)-longir(iir)) / nbande
145        totalemit(iir)=0.E+0
146        epir(iir)=0.E+0
147        omegir(iir)=0.E+0
148        gir(iir)=0.E+0
149c
150c.......................................................................
151c
152        DO ibande=1,nbande
153c
154c.......................................................................
155c
156          long=longir(iir) + (ibande-0.5E+0) * deltalong
157          CALL blackl(DBLE(long),DBLE(temp),tmp1)
158          emit=REAL(tmp1)
159c
160c.......................................................................
161c
162c********************************************************
163c interpolation
164      ilong=1
165      DO idata=2,ndata
166        IF (long.gt.longdata(idata)) ilong=idata
167      ENDDO
168      i1=ilong
169      i2=ilong+1
170      IF (i2.gt.ndata) i2=ndata
171      IF (long.lt.longdata(1)) i2=1
172      IF (i1.eq.i2) THEN
173        c1=1.E+0
174        c2=0.E+0
175      ELSE
176        c1=(longdata(i2)-long) / (longdata(i2)-longdata(i1))
177        c2=(longdata(i1)-long) / (longdata(i1)-longdata(i2))
178      ENDIF
179c********************************************************
180c
181          qextcorr=c1*qextcorrdata(i1)+c2*qextcorrdata(i2)
182          omeg=c1*omegdata(i1)+c2*omegdata(i2)
183          g=c1*gdata(i1)+c2*gdata(i2)
184c
185c.......................................................................
186c
187          totalemit(iir)=totalemit(iir)+deltalong*emit
188          epir(iir)=epir(iir)+deltalong*emit*qextcorr
189          omegir(iir)=omegir(iir)+deltalong*emit*omeg*qextcorr
190          gir(iir)=gir(iir)+deltalong*emit*omeg*qextcorr*g
191c
192c.......................................................................
193c
194        ENDDO
195c
196c.......................................................................
197c
198        gir(iir)=gir(iir)/omegir(iir)
199        omegir(iir)=omegir(iir)/epir(iir)
200        epir(iir)=epir(iir)/totalemit(iir)
201c
202c.......................................................................
203c
204      ENDDO
205c
206c......................................................................
207c
208c     Diagnostic de controle si on moyenne sur tout le spectre vis ou IR :
209c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
210c     tmp2=0.E+0
211c     DO iir=1,nir
212c       tmp2=tmp2+totalemit(iir)
213c     ENDDO
214c     tmp3=5.67E-8 * temp**4
215c     IF (abs((tmp2-tmp3)/tmp3).gt.0.05E+0) THEN
216c       PRINT *,'!!!! <---> il manque du Planck (voir moyenne.F)'
217c       PRINT *,'somme des bandes :',tmp2,'--- Planck:',tmp3
218c     ENDIF
219c
220c......................................................................
221c
222      RETURN
223      END
Note: See TracBrowser for help on using the repository browser.