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

Last change on this file since 3436 was 2803, checked in by aslmd, 2 years ago

Initialisation of Radiative Generic Condensable Aerosols

We can activate the scheme by putting aerogeneric = # of aerosols in callphys.def. This is the only needed thing for activating the radiative effects. They must be tracer in modern tracer.def

Commented out the abort if we use more than 4 aerosols

Added reading of optprop files for Radiative Generic Condensable Aerosols

We use the same file for IR and VI channel. For now, only MnS, Cr,Fe and Mg2SiO4 can be read. If you want to add another specie, check the code, it is explained how to quickly do that (right above the Initializations)

Added radii calculation for Radiative Generic Condensable aerosols

Changed the hardcoded size of the totalemit array

The hardcoded size is now 1900 instead of 100 so we don't exceed the array size when working at high spectral resolution (very rare case)

Added opacity computation for Radiative Generic Condensable aerosols

We do this computation in the same fashion as what's performed on water and dust.

switch iniaerosol and initracer order, to prepare for the RGCS scheme

Needed to switch the order of initialization so we can use the RGCS scheme without the assumption that ice and vap tracer of the same specie are following each other in the traceur.def file
File size: 9.9 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
58      INTEGER,PARAMETER :: nirmx=1900
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
69      INTEGER,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
85c check ordering of longdata
86      DO idata=1,ndata-1
87        IF (longdata(1).LT.longdata(ndata)) THEN
88          IF (.not.(longdata(idata).LT.longdata(idata+1))) THEN
89           call abort_physic("aerave_new",
90     &     "Non descending order in longdata",1)
91          ENDIF
92        ELSEIF (longdata(1).GT.longdata(ndata)) THEN
93          IF (.not.(longdata(idata).GT.longdata(idata+1))) THEN
94           call abort_physic("aerave_new",
95     &     "Non ascending order in longdata",1)
96          ENDIF
97        ENDIF
98      ENDDO
99c
100     
101       
102
103
104c********************************************************
105c interpolation
106c wavelengths (longdata) from data file in ascending order
107      IF (longdata(1).LT.longdata(ndata)) THEN
108        ilong=1
109        DO idata=2,ndata
110          IF (long.gt.longdata(idata)) ilong=idata
111        ENDDO
112        i1=ilong
113        i2=ilong+1
114        IF (i2.gt.ndata) i2=ndata
115        IF (long.lt.longdata(1)) i2=1
116        IF (i1.eq.i2) THEN
117          c1=1.E+0
118          c2=0.E+0
119        ELSE
120          c1=(longdata(i2)-long) / (longdata(i2)-longdata(i1))
121          c2=(longdata(i1)-long) / (longdata(i1)-longdata(i2))
122        ENDIF
123        qref=c1*epdata(i1)+c2*epdata(i2)
124        omegaref=c1*omegdata(i1)+c2*omegdata(i2)
125        factep=qref/epref
126        DO idata=1,ndata
127          qextcorrdata(idata)=epdata(idata)/factep
128        ENDDO
129c wavelengths (longdata) from data file in descending order
130      ELSEIF (longdata(1).GT.longdata(ndata)) THEN
131        ilong=1
132        DO idata=2,ndata
133          IF (long.lt.longdata(idata)) ilong=idata
134        ENDDO
135        i1=ilong+1
136        i2=ilong
137        IF (i1.gt.ndata) i1=ndata
138        IF (long.gt.longdata(1)) i1=1
139        IF (i1.eq.i2) THEN
140          c1=1.E+0
141          c2=0.E+0
142        ELSE
143          c1=(longdata(i2)-long) / (longdata(i2)-longdata(i1))
144          c2=(longdata(i1)-long) / (longdata(i1)-longdata(i2))
145        ENDIF
146        qref=c1*epdata(i1)+c2*epdata(i2)
147        omegaref=c1*omegdata(i1)+c2*omegdata(i2)
148        factep=qref/epref
149        DO idata=1,ndata
150          qextcorrdata(idata)=epdata(idata)/factep
151        ENDDO
152      ENDIF
153
154c********************************************************
155c.......................................................................
156c wavelengths (longdata) from data file in ascending order
157c.......................................................................
158      IF (longdata(1).LT.longdata(ndata)) THEN
159        DO iir=1,nir
160c
161c.......................................................................
162c
163          deltalong=(longir(iir+1)-longir(iir)) / nbande
164          totalemit(iir)=0.E+0
165          epir(iir)=0.E+0
166          omegir(iir)=0.E+0
167          gir(iir)=0.E+0
168c
169c.......................................................................
170c
171          DO ibande=1,nbande
172c
173c.......................................................................
174c
175            long=longir(iir) + (ibande-0.5E+0) * deltalong
176            CALL blackl(DBLE(long),DBLE(temp),tmp1)
177            emit=REAL(tmp1)
178c
179c.......................................................................
180c
181c interpolation
182            ilong=1
183            DO idata=2,ndata
184              IF (long.gt.longdata(idata)) ilong=idata
185            ENDDO
186            i1=ilong
187            i2=ilong+1
188            IF (i2.gt.ndata) i2=ndata
189            IF (long.lt.longdata(1)) i2=1
190            IF (i1.eq.i2) THEN
191              c1=1.E+0
192              c2=0.E+0
193            ELSE
194              c1=(longdata(i2)-long) / (longdata(i2)-longdata(i1))
195              c2=(longdata(i1)-long) / (longdata(i1)-longdata(i2))
196            ENDIF
197            qextcorr=c1*qextcorrdata(i1)+c2*qextcorrdata(i2)
198            omeg=c1*omegdata(i1)+c2*omegdata(i2)
199            g=c1*gdata(i1)+c2*gdata(i2)
200c
201c.......................................................................
202c
203            totalemit(iir)=totalemit(iir)+deltalong*emit
204            epir(iir)=epir(iir)+deltalong*emit*qextcorr
205            omegir(iir)=omegir(iir)+deltalong*emit*omeg*qextcorr
206            gir(iir)=gir(iir)+deltalong*emit*omeg*qextcorr*g
207c
208c.......................................................................
209c
210          ENDDO
211c
212c.......................................................................
213c
214          gir(iir)=gir(iir)/omegir(iir)
215          omegir(iir)=omegir(iir)/epir(iir)
216          epir(iir)=epir(iir)/totalemit(iir)
217c
218c.......................................................................
219c
220        ENDDO
221c.......................................................................
222c wavelengths (longdata) from data file in descending order
223c.......................................................................
224      ELSEIF (longdata(1).GT.longdata(ndata)) THEN
225        DO iir=1,nir
226c
227c.......................................................................
228c
229          deltalong=(longir(iir+1)-longir(iir)) / nbande
230          totalemit(iir)=0.E+0
231          epir(iir)=0.E+0
232          omegir(iir)=0.E+0
233          gir(iir)=0.E+0
234c
235c.......................................................................
236c
237          DO ibande=1,nbande
238c
239c.......................................................................
240c
241            long=longir(iir) + (ibande-0.5E+0) * deltalong
242            CALL blackl(DBLE(long),DBLE(temp),tmp1)
243            emit=REAL(tmp1)
244c
245c.......................................................................
246c
247c interpolation
248            ilong=1
249            DO idata=2,ndata
250              IF (long.lt.longdata(idata)) ilong=idata
251            ENDDO
252            i1=ilong+1
253            i2=ilong
254            IF (i1.gt.ndata) i1=ndata
255            IF (long.gt.longdata(1)) i1=1
256            IF (i1.eq.i2) THEN
257              c1=1.E+0
258              c2=0.E+0
259            ELSE
260              c1=(longdata(i2)-long) / (longdata(i2)-longdata(i1))
261              c2=(longdata(i1)-long) / (longdata(i1)-longdata(i2))
262            ENDIF
263            qextcorr=c1*qextcorrdata(i1)+c2*qextcorrdata(i2)
264            omeg=c1*omegdata(i1)+c2*omegdata(i2)
265            g=c1*gdata(i1)+c2*gdata(i2)
266c
267c.......................................................................
268c
269            totalemit(iir)=totalemit(iir)+deltalong*emit
270            epir(iir)=epir(iir)+deltalong*emit*qextcorr
271            omegir(iir)=omegir(iir)+deltalong*emit*omeg*qextcorr
272            gir(iir)=gir(iir)+deltalong*emit*omeg*qextcorr*g
273c
274c.......................................................................
275c
276          ENDDO
277c
278c.......................................................................
279c
280          gir(iir)=gir(iir)/omegir(iir)
281          omegir(iir)=omegir(iir)/epir(iir)
282          epir(iir)=epir(iir)/totalemit(iir)
283c
284c.......................................................................
285c
286        ENDDO
287      ENDIF
288c
289c********************************************************
290c
291c......................................................................
292c
293c     Diagnostic de controle si on moyenne sur tout le spectre vis ou IR :
294c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
295c     tmp2=0.E+0
296c     DO iir=1,nir
297c       tmp2=tmp2+totalemit(iir)
298c     ENDDO
299c     tmp3=5.67E-8 * temp**4
300c     IF (abs((tmp2-tmp3)/tmp3).gt.0.05E+0) THEN
301c       PRINT *,'!!!! <---> il manque du Planck (voir moyenne.F)'
302c       PRINT *,'somme des bandes :',tmp2,'--- Planck:',tmp3
303c     ENDIF
304c
305c......................................................................
306c
307      END
Note: See TracBrowser for help on using the repository browser.