source: trunk/LMDZ.VENUS/libf/phyvenus/aerave_new.F @ 3461

Last change on this file since 3461 was 2560, checked in by slebonnois, 3 years ago

SL: Implementation of SW computation based on generic model. Switch between this new SW module or old module that reads R. Haus tables implemented with a key (solarchoice)

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