source: trunk/libf/phytitan/optcv.v1 @ 6

Last change on this file since 6 was 3, checked in by slebonnois, 15 years ago

Creation de repertoires:

  • chantiers : pour communiquer sur nos projets de modifs
  • documentation : pour stocker les docs

Ajout de:

  • libf/phytitan : physique de Titan
  • libf/chimtitan: chimie de Titan
  • libf/phyvenus : physique de Venus
File size: 13.6 KB
RevLine 
[3]1      SUBROUTINE OPTCV(nmicro,IPRINT)
2
3
4#include "dimensions.h"
5#include "dimphy.h"
6#include "microtab.h"
7#include "clesphys.h"
8
9c   Argument:
10c   ---------
11      integer nmicro
12c   ---------
13
14      PARAMETER(NLAYER=llm,NLEVEL=NLAYER+1)
15      PARAMETER (NSPECI=46,NSPC1I=47,NSPECV=24,NSPC1V=25)
16      COMMON /ATM/ Z(NLEVEL),PRESS(NLEVEL),DEN(NLEVEL),TEMP(NLEVEL)
17
18      COMMON /GASS/ CH4(NLEVEL),XN2(NLEVEL),H2(NLEVEL),AR(NLEVEL)
19     & ,XMU(NLEVEL),GAS1(NLAYER),COLDEN(NLAYER)
20
21      COMMON /VISGAS/SOLARF(NSPECV),NTERM(NSPECV),PEXPON(NSPECV),
22     &         ATERM(4,NSPECV),BTERM(4,NSPECV)
23
24      COMMON /AERSOL/ RADIUS(NLAYER), XNUMB(NLAYER)
25     & , REALI(NSPECI), XIMGI(NSPECI), REALV(NSPECV), XIMGV(NSPECV)
26
27      COMMON /CLOUD/ RADCLD(NLAYER), XNCLD(NLAYER)
28     &             , RCLDI(NSPECI), XICLDI(NSPECI)
29     &             , RCLDV(NSPECV), XICLDV(NSPECV)
30
31      COMMON /TAUS/   TAUHI(klon,NSPECI), TAUCI(klon,NSPECI)
32     &               ,TAUGI(klon,NSPECI), TAURV(klon,NSPECV)
33     &               ,TAUHV(klon,NSPECV) ,TAUCV(klon,NSPECV)
34     &               ,TAUGV(klon,NSPECV)
35
36      COMMON /TAUD/   TAUHID(klon,NLAYER,NSPECI)
37     &               ,TAUGID(klon,NLAYER,NSPECI)
38     &               ,TAUHVD(klon,NLAYER,NSPECV)
39     &               ,TAUGVD(klon,NLAYER,NSPECV)
40
41      COMMON /OPTICV/ DTAUV(klon,NLAYER,NSPECV,4)
42     &               ,TAUV(klon,NLEVEL,NSPECV,4)
43     &               ,WBARV(klon,NLAYER,NSPECV,4)
44     &               ,COSBV(klon,NLAYER,NSPECV,4)
45
46      COMMON /SPECTV/ BWNV(NSPC1V),WNOV(NSPECV)
47     &               ,DWNV(NSPECV),WLNV(NSPECV)
48
49      COMMON /PLANT/ CSUBP,RSFI,RSFV,F0PI
50      COMMON /ADJUST/ RHCH4,FH2,FHAZE,FHVIS,FHIR,TAUFAC,RCLOUD,FARGON
51      COMMON /CONST/ RGAS,RHOP,PI,SIGMA
52      COMMON /traceurs/qaer(klon,nlayer,nqmx)
53      COMMON /part/ v(nqmx),r(nqmx),vrat,dr(nqmx),dv(nqmx)
54
55      REAL xv1(klev,NSPECV)
56      REAL xv2(klev,NSPECV)
57      REAL xv3(klev,NSPECV)
58
59      REAL QF1(nqmx,NSPECV),QF2(nqmx,NSPECV)
60      REAL QF3(nqmx,NSPECV),QF4(nqmx,NSPECV)
61      REAL QM1(nqmx,NSPECV),QM2(nqmx,NSPECV)
62      REAL QM3(nqmx,NSPECV),QM4(nqmx,NSPECV)
63
64      save qf1,qf2,qf3,qf4,qm1,qm2,qm3,qm4
65 
66      integer ioptv,iwarning
67      integer ig_
68      save ioptv,iwarning
69      data ioptv,iwarning/0,0/
70
71C*
72C*
73C THIS SUBROUTINE SETS THE OPTICAL CONSTANTS IN THE VISIBLE
74C IT CALCUALTES FOR EACH LAYER, FOR EACH SPECRAL INTERVAL IN THE VIS
75C LAYER: WBAR, DTAU, COSBAR
76C LEVEL: TAU
77C
78       sum=0.
79       PRINT*,'OPTCV'
80       print*,'ATTENTION, TAU UNIFORME DANS OPTCV'
81
82c      do nng=2,klon
83c       do i=1,klev           
84c        do j=1,nqmx
85c          sum=sum+qaer(nng,i,j)*r(j)**3.*1.3333*3.1415*1000.
86c        enddo
87c       enddo
88c       enddo
89c       print*,sum/(klon-1),'SOMME COLONNE/OPTCV'
90
91c      open (unit=1,file='xsetupv')
92c       do j=1,nspecv
93c        read(1,*) a
94c        do i=1,klev
95c            read(1,*) xv1(i,j),xv2(i,j),xv3(i,j)
96c        enddo
97c       enddo
98c       close(1)
99
100      DO 130 K=1,NSPECV
101C LETS USE THE OPTICAL CONSTANTS FOR THOLIN
102      CALL THOLIN(WLNV(K),TNR,TNI)
103      REALV(K)=TNR
104      XIMGV(K)=TNI*FHVIS
105C BUT WE NOW USE THE GEOMETRIC ALBEDO FITTED RESULTS
106C      XIMGV(K)=FITEDT(WLNV(K))
107C      XIMGV(K)=FITEDN(WLNV(K))
108C THE CLOUD IS CLEAR IN THE VISIBLE
109      RCLDV(K)=1.27
110      XICLDV(K)=1.E-7
111 130  CONTINUE
112C
113c******* DEBUT DES BOUCLES ************************
114c     PRINT*, 'AEROSOLS EN VISIBLE'
115
116      DO 101 ig=1,klon       !c! BOUCLE SUR GRILLE HORIZONTALE
117
118      DO 100 K=1,NSPECV         !b! BOUCLE SUR LAMBDA
119
120C ZERO THE COLUMN OPTICAL DEPTHS OF EACH TYPE
121C ??FLAG? THE OPTICAL DEPTH OF THE TOP OF THE MODEL
122C MAY NOT BE ZERO.
123      TAURV(ig,K)=0.
124      TAUHV(ig,K)=0.            ! INTEGRATED TAU.......INITIALIZATION.
125      TAUCV(ig,K)=0.            ! Rayleigh, Haze, Cloud, Gas
126      TAUGV(ig,K)=0.            !   sca,    abs,  abs  , abs
127
128      DO 100 J=1,NLAYER         !a! BOUCLE SUR L"ALTITUDE
129
130C #1:                   HAZE
131c---------------------------
132
133c     CALL THE MIE CODE TO GIVE THE AEROSOL PROPERTIES
134c     USE XFRAC FOR FRACTAL AEROSOLS PROPERTIES AT LAMBDA < 2. um
135
136
137
138
139c                    /\
140c                   /  \
141c                  /    \
142c                 / _O   \
143c                / |/     \
144c               /  / \     \
145c              /   |\ \/\   \
146c             /    || /  \   \
147c             ----------------
148c            |     WARNING    |
149c            |    SLOW DOWN   |
150c             ----------------
151
152
153
154
155c*********** EN TRAVAUX ***************************
156 
157         TAEROS=0.
158         TAEROSCAT=0.
159         CBAR=0.
160
161
162      DO inq=1,nmicro         !BOUCLE SUR LES nmicro TAILLE D"AEROSOLS
163
164
165            IF (R(inq).lt.RF(inq)) THEN    ! aerosols spheriques
166
167             
168            if(ioptv.eq.0.and.J.eq.1) then
169c                  CALL XMIE(R(inq)*1.e6,REALV(K),XIMGV(K),
170c    &             QEXT,QSCT,QABS,QBAR,WNOV(K))
171
172                   CALL CMIE(1.E-2/WNOV(K),REALV(K),XIMGV(K),R(inq),
173     &             QEXT,QSCT,QABS,QBAR)
174
175c       print*,'inq=',inq,' QM1=',QM1(inq,K),' QEXT=',QEXT
176
177              QM1(inq,K)=QEXT
178              QM2(inq,K)=QSCT
179              QM3(inq,K)=QABS
180              QM4(inq,K)=QBAR
181            endif
182
183
184      if (microfi.eq.1) then
185         ig_=ig
186      else
187         ig_=12
188      endif
189
190       TAEROS=QM1(inq,K)*qaer(ig_,NLAYER+1-J,inq)*1.e-4+TAEROS
191       TAEROSCAT=QM2(inq,K)*qaer(ig_,NLAYER+1-J,inq)*1.e-4+TAEROSCAT
192       CBAR=CBAR+QM4(inq,K)*QM2(inq,K)*qaer(ig_,NLAYER+1-J,inq)*1.e-4
193
194            ELSE                        ! aerosols fractals
195
196               XMONO=(R(inq)/RF(inq))**3.
197               XRULE=1.
198
199            if(XMONO.gt.16384./1.5) then
200             XRULE=(XMONO/16384.)
201             XMONO=16384.
202            endif
203
204            if(ioptv.eq.0.and.J.eq.1) then
205
206        CALL OPTFRAC(XMONO,10000./WNOV(K)
207     &                        ,QEXT,QSCT,QABS,QBAR)
208
209c       CALL CFFFV11(1.e-2/WNOV(K),REALV(K),XIMGV(K),RF(inq),2.
210c    &   ,XMONO,QSCT,QEXT,QABS,QBAR)
211
212
213              QF1(inq,K)=QEXT*XRULE
214              QF2(inq,K)=QSCT*XRULE
215              QF3(inq,K)=QABS*XRULE
216              QF4(inq,K)=QBAR
217
218c       print*,'inq=',inq,' QF1=',QF1(inq,K),' QEXT=',QEXT,' XRULE=',XRULE
219                   
220            endif
221
222
223      if (microfi.eq.1) then
224         ig_=ig
225      else
226         ig_=12
227      endif
228
229        TAEROS=QF1(inq,K)*qaer(ig_,NLAYER+1-J,inq)+TAEROS
230        TAEROSCAT=QF2(inq,K)*qaer(ig_,NLAYER+1-J,inq)+TAEROSCAT
231        CBAR=CBAR+QF4(inq,K)*QF2(inq,K)*qaer(ig_,NLAYER+1-J,inq)
232
233
234           ENDIF
235
236
237       ENDDO    ! nmicro
238
239
240       CBAR=CBAR/TAEROSCAT
241
242        DELTAZ=Z(J)-Z(J+1)
243
244c --------------------------------------------------------------------
245c profil brume Pascal: fit T (sauf tropopause) et albedo
246c -------------------
247        if( cutoff.eq.1) then
248         IF(PRESS(J).gt.9.e-3) THEN
249          TAEROS=TAEROSM1*DELTAZ/DELTAZM1*0.85
250          TAEROSCAT=TAEROSCATM1*DELTAZ/DELTAZM1*0.85
251c         TAEROS=0.
252c         TAEROSCAT=0.
253         ENDIF
254
255         IF(PRESS(J).gt.1.e-1) THEN
256          TAEROS=TAEROSM1*DELTAZ/DELTAZM1*1.15
257          TAEROSCAT=TAEROSCATM1*DELTAZ/DELTAZM1*1.15
258c         TAEROS=0.
259c         TAEROSCAT=0.
260         ENDIF
261        endif !cutoff=1
262
263c profil brume pour fit T (y compris tropopause), mais ne fit plus albedo...
264c -----------------------
265        if( cutoff.eq.2) then
266         IF(PRESS(J).gt.1.e-1) THEN
267          TAEROS=0.
268          TAEROSCAT=0.
269         ENDIF
270        endif !cutoff=2
271c --------------------------------------------------------------------
272
273         TAEROSM1=TAEROS
274         TAEROSCATM1=TAEROSCAT
275         DELTAZM1=DELTAZ
276
277
278       IF (TAEROSCAT.le.0.) CBAR=0.
279
280c          if(ig.eq.12) then
281c          if(j.eq.1) print*,'NEWK',wlnv(k)
282c          print*,j,TAEROS,xv1(j,k),'   ', TAEROSCAT/TAEROS,
283c    &     xv2(j,k)/xv1(j,k),'  ',CBAR,xv3(j,k)
284c          print*,' '
285c          endif
286           
287
288c      print*,'HERE; MCKAY AEROSOLS VIS'
289
290c      TAEROSCAT=xv2(j,k)
291c      TAEROS=xv1(j,k)
292c      CBAR=xv3(j,k)
293
294c     if (ig.eq.1) then
295c     if (k.eq.NSPECV/2) then   
296c      print*,'@VI',K,J,TAEROS,TAEROSCAT,CBAR
297c      print*,'@  ',K,J,QF1(1,K),QF2(1,K),qaer(12,NLAYER+1-J,1)
298c      print*,'@  ',K,J,QF1(3,K),QF2(3,K),qaer(12,NLAYER+1-J,3)
299c      print*,'@  ',K,J,QF1(5,K),QF2(5,K),qaer(12,NLAYER+1-J,5)
300c      print*,'@  ',K,J,QF1(7,K),QF2(7,K),qaer(12,NLAYER+1-J,7)
301c      print*,'@  ',K,J,QF1(9,K),QF2(9,K),qaer(12,NLAYER+1-J,9)
302c      print*
303c     endif
304c     endif
305
306
307
308c*********** EN TRAVAUX ***************************
309
310C #2:                   RAYLEIGH
311c-------------------------------
312
313C RAYLEIGH SCATTERING STRAIGHT FROM HANSEN AND TRAVIS...SEE NOTES
314C RATIOED BY THE LAYER COLUMN NUMBER TO THE TOTAL
315C COLUMN NUMBER ON EARTH. CM-2
316C THIS IS THE SCATTERING BY THE ATMOSPHERE
317
318      TAURAY=(COLDEN(J)*28.9/(XMU(J)*1013.25))*
319     &(.008569/WLNV(K)**4)*(1.+.0113/WLNV(K)**2+.00013/WLNV(K)**4)
320
321c       PRINT*,WLNV(K)
322c      COLX=0.
323c      COLP=0.
324c      COLT=0.
325c     DO IU=1,NLAYER
326c      COLP=COLDEN(IU)*1.e+1*1.35+COLP
327c     TAURAY=(COLDEN(IU)*28.9/(XMU(IU)*1013.25))*
328c    & (.008569/WLNV(K)**4)*(1.+.0113/WLNV(K)**2
329c    & +.00013/WLNV(K)**4)
330c      COLT=COLT+TAURAY
331c      COLX=COLDEN(IU)*1.e+1/(1.E5*28./22.4E3)*1.e-1*0.0933e-1+COLX
332c                           |   
333c                           |   
334c           g/cm2->kg/m2    |  m2/kg   
335c      Print*,IU, tauray,
336c    &   COLDEN(IU)*1.e+1/(1.E5*28./22.4E3)*1.e-1*0.543e-1
337c     ENDDO
338c       PRINT*,COLP,' PRESSURE AT GROUND;'
339c       PRINT*,COLX,' TAU_GAS AT GROUND;'
340c       print*,colt,colx,' COLT, COLX'
341c      STOP
342                                                       
343c       DZ=Z(J)-Z(J+1)
344c     PRINT*, Z(J),WLNV(K),
345c    &(28.9/(XMU(J)*1013.25))*(.008569/WLNV(K)**4)*
346c    &(1.+.0113/WLNV(K)**2+.00013/WLNV(K)**4)
347c    & ,COLDEN(J)/DZ/100000.,
348c    &(28.9/(XMU(J)*1013.25))*(.008569/WLNV(K)**4)*
349c    &(1.+.0113/WLNV(K)**2+.00013/WLNV(K)**4)
350c    & *COLDEN(J)/DZ/100000.
351   
352 
353
354C #3:                   CLOUD
355c----------------------------
356
357C NEXT COMPUTE TAU CLOUD
358
359      TAUCLD=0.0
360c             XNCLD(J)=0.
361      IF ( XNCLD(J) .GT. 0. .and .taufac.gt.0.) THEN
362                CALL XMIE(RADCLD(J),RCLDV(K),XICLDV(K),
363     &                         QEXTC,QSCTC,QABSC,CBARC,WNOV(K))
364                TAUCLD=QEXTC*XNCLD(J)         
365      ENDIF
366C
367      TAURV(ig,K)=TAURV(ig,K)+TAURAY
368      TAUGVD(ig,J,K)=TAURV(ig,K)
369
370      TAUHV(ig,K)=TAUHV(ig,K)+TAEROS          ! INTEGRATED Quant.
371      TAUHVD(ig,J,K)=TAUHV(ig,K)
372
373      TAUCV(ig,K)=TAUCV(ig,K)+TAUCLD
374
375C #4:                  TAUGAS
376C----------------------------
377
378C LOOP OVER THE NTERMS
379C THIS IS THE ABSORPTION BY THE ATMOSPHERE (METHANE)
380
381
382      DO 909 NT=1,NTERM(K)
383      TAUGAS=COLDEN(J)*GAS1(J)*BTERM(NT,K)*
384     &  (   (PRESS(J+1) + PRESS(J))*.5  )**PEXPON(K)
385
386
387C COMPUTE THE AVERAGE COSBAR AND WBAR
388C&&
389
390c     CBAR=MIN(1.0,1.05*CBAR)       ! THE HAZE FORWARD SCATTERING 5%(WHY?)
391      COSBV(ig,J,K,NT)=(CBAR*TAEROSCAT + CBARC*TAUCLD)
392     &  /(TAEROSCAT+TAUCLD+TAURAY)     !CBAR_RAY=0.
393c        print*,'CBV',J,K,NT,CBAR,TAEROSCAT,CBARC,TAUCLD       
394
395      DTAUV(ig,J,K,NT)=TAUGAS+TAEROS+TAURAY+TAUCLD       !TOTAL TAU_EXT
396      TAUGV(ig,K)=TAUGV(ig,K)+TAUGAS*ATERM(NT,K)         !TAU_ABS_METH INTEG.
397
398C WE LET W RAYLEIGH BE .999 OR W=1 WHEN ONLY RAYLEIGH=PROBLEM FOR TRID
399c WE HAVE ASSUMED ABOVE THAT COSBAR FOR RAYLEIGH IS ZERO.
400c     if (ig.eq.1) then
401c     if (k.eq.NSPECV/2) then   
402c      print*,'@VI',K,J,DTAUV(ig,J,K,1),TAUGAS,TAEROS,TAUCLD
403c     endif
404c     endif
405
406
407c***************** ECHANGE
408c     WBARV(J,K,NT)=(QSCT*XNUMB(J)+TAURAY*0.9999999 + QSCTC*XNCLD(J) )
409c****************
410      WBARV(ig,J,K,NT)=(TAEROSCAT+TAURAY*0.9999999 + QSCTC*XNCLD(J) )
411c     WBARV(ig,J,K,NT)=(TAEROSCAT+TAURAY*0.9999999 )
412     &            /(TAUGAS+TAEROS+TAURAY+TAUCLD)
413c****************
414      IF((TAEROS+TAUCLD+TAURAY+TAUCLD).le.0.) WBARV(ig,J,K,NT)=0.
415      IF((TAEROS+TAUCLD+TAURAY).le.0.) COSBV(ig,J,K,NT)=0.
416
417c     print*,'WBV',J,K,NT,TAEROSCAT,TAURAY,QSCTC*XNCLD(J)
418c     print*,'WBV',J,K,NT,TAEROS,TAUGAS,TAURAY,TAUCLD
419c     print*,Z(j),J,K,NT,TAUV(1,j,K,NT),WBARV(1,j,K,NT),COSBV(1,j,K,NT)
420
421 909  CONTINUE
422      TAUGVD(ig,J,K)=TAUGVD(ig,J,K)+TAUGV(ig,K)
423 100  CONTINUE
424       ioptv=1
425 101  CONTINUE
426c HERE END OF THE LOOPS *******
427c******************************
428         
429        DO 102 ig=1,klon
430
431C TOTAL EXTINCTION OPTICAL DEPTHS
432          DO 119 K=1,NSPECV
433C LOOP OVER NTERMS
434           DO 119 NT=1,NTERM(K)
435           TAUV(ig,1,K,NT)=0.0
436             DO 119 J=1,NLAYER
437             TAUV(ig,J+1,K,NT)=TAUV(ig,J,K,NT)+DTAUV(ig,J,K,NT)
438 119     CONTINUE
439
440c       print*,'SETUP'
441c      do i=1,NSPECV
442c      print*,WLNV(i)
443c       do j=1,NLAYER+1
444c       print*,Z(j),TAUV(1,j,i,1),WBARV(1,j,i,1),COSBV(1,j,i,1)
445c       enddo
446c      enddo
447c
448c     IF (IPRINT .GT. 1) THEN
449c           NT=1
450c     IF (2 .GT. 1) THEN
451c          WRITE (6,120)
452c 120      FORMAT(///'  OPTICAL CONSTANTS IN THE VISIBLE (@EQUATOR) ')
453c          WRITE(6,*) 'latitude:',ig
454c          DO 200 K=1,NSPECV
455c          WRITE (6,190)
456c          WRITE (6,210)K,WLNV(K),WNOV(K),BWNV(K)
457c    &    ,BWNV(K)+DWNV(K),DWNV(K)
458c          WRITE (6,230)REALV(K),XIMGV(K)
459c          DO 195 J=1,NLAYER,NLAYER
460C RECALCULATE FOR PRINT OUT ONLY, ONLY FIRST NTERM AT ig=12 (EQUATOR)
461c          WRITE (6,220)XNUMB(J), WBARV(ig,J,K,NT),COSBV(ig,J,K,NT)
462c    &      ,DTAUV(ig,J,K,NT),TAUV(ig,J,K,NT)
463c 195      CONTINUE
464c          WRITE (6,240) TAUV(IG,NLEVEL,K,NT)
465c 200      CONTINUE
466c     END IF
467  102 CONTINUE
468  210 FORMAT(1X,I3,F10.3,F10.2,F10.2,'-',F8.2,F10.3)
469  190 FORMAT(1X//'  SNUM  MICRONS   WAVENU   INTERVAL    DELTA-WN')
470  230 FORMAT(1X,'NREAL(LAYER)= ',1PE10.3,' NIMG(LAYER)= ',E10.3/
471     &' #AEROSOLS   WBAR  COSBAR       DTAU     TAU'
472     & ,9X,'RAY     GAS    AEROSOL')
473  220 FORMAT(8(1X,F9.3))
474  240 FORMAT(41X,F9.3)
475
476           print*,"TAUV(1400,:,10,2)=",TAUV(1400,:,10,2)
477           print*,"DTAUV(1400,:,10,2)=",DTAUV(1400,:,10,2)
478c      ioptv=1
479       PRINT*, 'FIN OPTCV'
480       stop
481      RETURN
482      END
Note: See TracBrowser for help on using the repository browser.