source: trunk/LMDZ.TITAN/libf/phytitan/optcv_1pt_2.F @ 201

Last change on this file since 201 was 175, checked in by slebonnois, 14 years ago

S.LEBONNOIS:

  • Revision majeure de la physique Titan => ajout des nuages version 10 bins (Jeremie Burgalat) Cette version reste a tester mais avec clouds=0, on reste sur l'ancienne.
  • Quelques ajouts dans la doc.
File size: 12.1 KB
Line 
1      SUBROUTINE optcv_1pt2(zqaer_1pt,rcdb,xfrb,ioptv,IPRINT)
2
3      use dimphy
4#include "dimensions.h"
5#include "microtab.h"
6#include "clesphys.h"
7
8      PARAMETER(NLAYER=llm,NLEVEL=NLAYER+1)
9      PARAMETER (NSPECI=46,NSPC1I=47,NSPECV=24,NSPC1V=25)
10
11c   Arguments:
12c   ---------
13      integer IPRINT,ioptv
14C ioptv: premier appel, on ne calcule qu'une fois les QM et QF
15* nrad dans microtab.h
16      real   zqaer_1pt(NLAYER,2*nrad)
17#include "optcv_1pt.h"
18c   ---------
19
20      COMMON /ATM/ Z(NLEVEL),PRESS(NLEVEL),DEN(NLEVEL),TEMP(NLEVEL)
21
22      COMMON /GASS/ CH4(NLEVEL),XN2(NLEVEL),H2(NLEVEL),AR(NLEVEL)
23     & ,XMU(NLEVEL),GAS1(NLAYER),COLDEN(NLAYER)
24
25      COMMON /VISGAS/SOLARF(NSPECV),NTERM(NSPECV),PEXPON(NSPECV),
26     &         ATERM(4,NSPECV),BTERM(4,NSPECV)
27
28      COMMON /AERSOL/ RADIUS(NLAYER), XNUMB(NLAYER)
29     & , REALI(NSPECI), XIMGI(NSPECI), REALV(NSPECV), XIMGV(NSPECV)
30
31      COMMON /CLOUD/
32     &               RCLDI(NSPECI), XICLDI(NSPECI)
33     &             , RCLDV(NSPECV), XICLDV(NSPECV)
34     &             , RCLDI2(NSPECI), XICLDI2(NSPECI)
35     &             , RCLDV2(NSPECV), XICLDV2(NSPECV)
36
37      COMMON /SPECTV/ BWNV(NSPC1V),WNOV(NSPECV)
38     &               ,DWNV(NSPECV),WLNV(NSPECV)
39
40      COMMON /PLANT/ CSUBP,RSFI,RSFV,F0PI
41      COMMON /ADJUST/ RHCH4,FH2,FHAZE,FHVIS,FHIR,TAUFAC,RCLOUD,FARGON
42      COMMON /CONST/ RGAS,RHOP,PI,SIGMA
43* nrad dans microtab.h
44      COMMON /part/ v(nrad),rayon(nrad),vrat,dr(nrad),dv(nrad)
45
46      REAL QF1(nrad,NSPECV),QF2(nrad,NSPECV)
47      REAL QF3(nrad,NSPECV),QF4(nrad,NSPECV)
48      REAL QM1(nrad,NSPECV),QM2(nrad,NSPECV)
49      REAL QM3(nrad,NSPECV),QM4(nrad,NSPECV)
50      REAL QC1(nrad,NSPECV),QC2(nrad,NSPECV)
51      REAL QC3(nrad,NSPECV),QC4(nrad,NSPECV)
52
53c---- NUAGES
54      real TNUABS,TNUSCAT
55      real   rcdb(NLAYER), xfrb(NLAYER,4)
56
57      save qf1,qf2,qf3,qf4,qm1,qm2,qm3,qm4,qc1,qc2,qc3,qc4
58
59
60C*
61C THIS SUBROUTINE SETS THE OPTICAL CONSTANTS IN THE VISIBLE
62C IT CALCUALTES FOR EACH LAYER, FOR EACH SPECRAL INTERVAL IN THE VIS
63C LAYER: WBAR, DTAU, COSBAR
64C LEVEL: TAU
65C
66C ZERO THE COLUMN OPTICAL DEPTHS OF EACH TYPE
67C ??FLAG? THE OPTICAL DEPTH OF THE TOP OF THE MODEL
68C MAY NOT BE ZERO.
69
70c******* DEBUT DES BOUCLES ************************
71      DO 100 K=1,NSPECV         !b! BOUCLE SUR LAMBDA
72
73      TAURV_1pt(K)=0.
74      TAUHV_1pt(K)=0.            ! INTEGRATED TAU.......INITIALIZATION.
75      TAUCV_1pt(K)=0.            ! Rayleigh, Haze, Cloud, Gas
76      TAUGV_1pt(K)=0.            !   sca,    abs,  abs  , abs
77
78      DO 100 J=1,NLAYER         !a! BOUCLE SUR L"ALTITUDE
79
80C #1:                   HAZE
81c---------------------------
82
83c     CALL THE MIE CODE TO GIVE THE AEROSOL PROPERTIES
84c     USE XFRAC FOR FRACTAL AEROSOLS PROPERTIES AT LAMBDA < 2. um
85
86
87
88
89c                    /\
90c                   /  \
91c                  /    \
92c                 / _O   \
93c                / |/     \
94c               /  / \     \
95c              /   |\ \/\   \
96c             /    || /  \   \
97c             ----------------
98c            |     WARNING    |
99c            |    SLOW DOWN   |
100c             ----------------
101
102
103
104
105c*********** EN TRAVAUX ***************************
106 
107         TAEROS=0.
108         TAEROSCAT=0.
109         CBAR=0.
110
111c       print*,"rayon=",rayon
112c       print*,"RF=",RF
113
114      DO inq=1,nrad         !BOUCLE SUR LES TAILLE D"AEROSOLS
115
116
117            IF (rayon(inq).lt.RF(inq)) THEN    ! aerosols spheriques
118
119             
120            if(ioptv.eq.0.and.J.eq.1) then
121c                  CALL XMIE(rayon(inq)*1.e6,REALV(K),XIMGV(K),
122c    &             QEXT,QSCT,QABS,QBAR,WNOV(K))
123
124                CALL CMIE(1.E-2/WNOV(K),REALV(K),XIMGV(K),rayon(inq),
125     &          QEXT,QSCT,QABS,QBAR)
126
127c       print*,'inq=',inq,' QM1=',QM1(inq,K),' QEXT=',QEXT
128
129              QM1(inq,K)=QEXT
130              QM2(inq,K)=QSCT
131              QM3(inq,K)=QABS
132              QM4(inq,K)=QBAR
133            endif
134
135       TAEROS=QM1(inq,K)*zqaer_1pt(NLAYER+1-J,inq)*1.e-4+TAEROS
136       TAEROSCAT=QM2(inq,K)*zqaer_1pt(NLAYER+1-J,inq)*1.e-4+TAEROSCAT
137       CBAR=CBAR+QM4(inq,K)*QM2(inq,K)*zqaer_1pt(NLAYER+1-J,inq)*1.e-4
138
139            ELSE                        ! aerosols fractals
140
141               XMONO=(rayon(inq)/RF(inq))**3.
142               XRULE=1.
143
144            if(XMONO.gt.16384./1.5) then
145             XRULE=(XMONO/16384.)
146             XMONO=16384.
147            endif
148
149            if(ioptv.eq.0.and.J.eq.1) then
150
151c        CALL OPTFRAC(XMONO,10000./WNOV(K)
152c     &                        ,QEXT,QSCT,QABS,QBAR)
153
154        CALL CFFFV11(1.e-2/WNOV(K),REALV(K),XIMGV(K),RF(inq),2.
155     &   ,XMONO,QSCT,QEXT,QABS,QBAR)
156
157
158              QF1(inq,K)=QEXT*XRULE
159              QF2(inq,K)=QSCT*XRULE
160              QF3(inq,K)=QABS*XRULE
161              QF4(inq,K)=QBAR
162
163c       print*,'inq=',inq,' QF1=',QF1(inq,K),' QEXT=',QEXT,' XRULE=',XRULE
164                   
165            endif
166
167        TAEROS=QF1(inq,K)*zqaer_1pt(NLAYER+1-J,inq)+TAEROS
168        TAEROSCAT=QF2(inq,K)*zqaer_1pt(NLAYER+1-J,inq)+TAEROSCAT
169        CBAR=CBAR+QF4(inq,K)*QF2(inq,K)*zqaer_1pt(NLAYER+1-J,inq)
170
171           ENDIF
172
173       ENDDO    ! nrad
174
175       IF(TAEROSCAT.le.0.) then
176        CBAR=0.
177       ELSE
178        CBAR=CBAR/TAEROSCAT
179       ENDIF
180
181        DELTAZ=Z(J)-Z(J+1)
182
183c --------------------------------------------------------------------
184c profil brume Pascal: fit T (sauf tropopause) et albedo
185c -------------------
186        if( cutoff.eq.1) then
187         IF(PRESS(J).gt.9.e-3) THEN
188          TAEROS=TAEROSM1*DELTAZ/DELTAZM1*0.85
189          TAEROSCAT=TAEROSCATM1*DELTAZ/DELTAZM1*0.85
190c         TAEROS=0.
191c         TAEROSCAT=0.
192         ENDIF
193
194         IF(PRESS(J).gt.1.e-1) THEN
195          TAEROS=TAEROSM1*DELTAZ/DELTAZM1*1.15
196          TAEROSCAT=TAEROSCATM1*DELTAZ/DELTAZM1*1.15
197c         TAEROS=0.
198c         TAEROSCAT=0.
199         ENDIF
200        endif !cutoff=1
201
202c profil brume pour fit T (y compris tropopause), mais ne fit plus albedo...
203c -----------------------
204        if( cutoff.eq.2) then
205         IF(PRESS(J).gt.1.e-1) THEN
206          TAEROS=0.
207          TAEROSCAT=0.
208         ENDIF
209        endif !cutoff=2
210c --------------------------------------------------------------------
211
212         TAEROSM1=TAEROS
213         TAEROSCATM1=TAEROSCAT
214         DELTAZM1=DELTAZ
215
216
217       IF (TAEROSCAT.le.0.) CBAR=0.
218
2191699  FORMAT(a3,2I3,3(ES15.7,1X))
220
221c*********** EN TRAVAUX ***************************
222
223C #2:                   RAYLEIGH
224c-------------------------------
225
226C RAYLEIGH SCATTERING STRAIGHT FROM HANSEN AND TRAVIS...SEE NOTES
227C RATIOED BY THE LAYER COLUMN NUMBER TO THE TOTAL
228C COLUMN NUMBER ON EARTH. CM-2
229C THIS IS THE SCATTERING BY THE ATMOSPHERE
230
231      TAURAY=(COLDEN(J)*28.9/(XMU(J)*1013.25))*
232     &(.008569/WLNV(K)**4)*(1.+.0113/WLNV(K)**2+.00013/WLNV(K)**4)
233
234
235C #3:                   CLOUD
236c----------------------------
237C NEXT COMPUTE TAU CLOUD
238c
239c  Menu special :
240c  Afin d'eviter la surcharge de calcul on ne calcule les
241c  propriétes optiques des nuages qu'une seule fois
242c  avec un rayon de particule effectif de 3um et une composition
243c  de goutte : 90% CH4 / 10% NOYAUX
244c  Puis on ajute les section efficace par la surface reelle de
245c  la goutte.
246c
247c  ---> A TESTER !!!!
248c
249       IF (clouds.eq.0) THEN
250         CNBAR=0.
251         TNUSCAT=0.
252         TNUABS=0.
253       ELSE
254         IF (ioptv.eq.0.and.j.eq.1) THEN !--> au premier appel
255           QEXTC=0.
256           QSCTC=0.
257           QABSC=0.
258           CBARC=0.
259           DO inq=1,nrad         !BOUCLE SUR LES NQMX TAILLE D"AEROSOLS
260             QC1(inq,K)=0.
261             QC2(inq,K)=0.
262             QC3(inq,K)=0.
263             QC4(inq,K)=0.
264           ENDDO
265** OPTICAL CONSTANT : MIXING RULES
266** Fraction volumique fixe :
267** 10% noyaux.
268** 90% methane.
269           XNR = 0.5 * REALI(K)
270     &         + 0.5 * RCLDI(K)
271           XNI = 0.5 * XIMGI(K)
272     &         + 0.5 * XICLDI(K)
273**
274**   Efficacite : particule de 3um de rayon
275           CALL CMIE(1.E-2/WNOV(K),XNR,XNI,3.e-6,
276     &                QEXTC,QSCTC,QABSC,CBARC)
277
278           DO inq=1,nrad
279             QC1(inq,K)=QEXTC/xnuf
280             QC2(inq,K)=QSCTC/xnuf
281             QC3(inq,K)=QABSC/xnuf
282             QC4(inq,K)=CBARC
283           ENDDO
284         ENDIF   ! ioptv = 0
285         TNUABS=0.
286         TNUSCAT=0.
287         CNBAR=0.
288         IF (rcdb(nlayer+1-J).gt.1.1e-10) THEN
289           DO inq=1,nrad
290             TNUABS=QC1(inq,K)*(rcdb(nlayer+1-J)/3.e-6)**2.*1.e-4*
291     &              zqaer_1pt(NLAYER+1-J,inq+nrad) +
292     &              TNUABS
293             TNUSCAT=QC2(inq,K)*(rcdb(nlayer+1-J)/3.e-6)**2.*1.e-4*
294     &               zqaer_1pt(NLAYER+1-J,inq+nrad) +
295     &               TNUSCAT
296             CNBAR=QC4(inq,K)*QC2(inq,K)*(rcdb(nlayer+1-J)/3.e-6)**2.*
297     &             1.e-4*zqaer_1pt(NLAYER+1-J,inq+nrad) +
298     &             CNBAR
299           ENDDO
300         ENDIF
301
302         IF(TNUSCAT.EQ.0.) THEN
303           CNBAR=0.
304         ELSE
305           CNBAR=CNBAR/TNUSCAT
306         ENDIF
307       ENDIF  ! Cond. CLD
308
309       TAUCV_1pt(K)=TAUCV_1pt(K)+TNUABS
310       TAUCVD_1pt(J,K)=TAUCV_1pt(K)
311
312       TAURV_1pt(K)=TAURV_1pt(K)+TAURAY
313       TAUGVD_1pt(J,K)=TAURV_1pt(K)
314
315       TAUHV_1pt(K)=TAUHV_1pt(K)+TAEROS ! INTEGRATED Quant.
316       TAUHVD_1pt(J,K)=TAUHV_1pt(K)
317
318
319
320C #4:                  TAUGAS
321C----------------------------
322
323C LOOP OVER THE NTERMS
324C THIS IS THE ABSORPTION BY THE ATMOSPHERE (METHANE)
325
326
327       DO 909 NT=1,NTERM(K)
328         TAUGAS=COLDEN(J)*GAS1(J)*BTERM(NT,K)*
329     &   (   (PRESS(J+1) + PRESS(J))*.5  )**PEXPON(K)
330
331
332*  COSBV ET COSBVP
333*-----------------
334
335         IF(TAEROSCAT+TNUSCAT+TAURAY .ne. 0.) THEN
336           COSBV_1pt(J,K,NT)=(CBAR*TAEROSCAT + CNBAR*TNUSCAT)
337     &     /(TAEROSCAT+TNUSCAT+TAURAY) !CBAR_RAY=0.
338         ELSE
339           COSBV_1pt(J,K,NT)=0.
340         ENDIF
341
342         IF(TAEROSCAT+TAURAY .ne. 0.) THEN
343           COSBVP_1pt(J,K,NT)=(CBAR*TAEROSCAT)
344     &     /(TAEROSCAT+TAURAY) !CBAR_RAY=0.
345         ELSE
346           COSBVP_1pt(J,K,NT)=0.
347         ENDIF
348
349*  DTAUV ET DTAUVP
350*-----------------
351
352         DTAUV_1pt(J,K,NT) =TAUGAS+TAEROS+TAURAY+TNUABS !TAU_ABS_METH
353         DTAUVP_1pt(J,K,NT)=TAUGAS+TAEROS+TAURAY       !TAU_ABS_METH
354
355         TAUGV_1pt(K)=TAUGV_1pt(K)+TAUGAS*ATERM(NT,K) !INTEG.
356
357*  WBARV ET WBARVP
358*-----------------
359
360         IF(TAUGAS+TAEROS+TAURAY+TNUABS .ne.  0.) THEN
361           WBARV_1pt(J,K,NT)=(TAEROSCAT+TAURAY*0.9999999 + TNUSCAT)
362     &     /(TAUGAS+TAEROS+TAURAY+TNUABS)
363         ELSE
364           WBARV_1pt(J,K,NT)=0.
365         ENDIF
366
367         IF(TAUGAS+TAEROS+TAURAY .ne.  0.) THEN
368           WBARVP_1pt(J,K,NT)=(TAEROSCAT+TAURAY*0.9999999 )
369     &     /(TAUGAS+TAEROS+TAURAY)
370         ELSE
371           WBARVP_1pt(J,K,NT)=0.
372         ENDIF
373   
374 909   CONTINUE
375 
376       TAUGVD_1pt(J,K)=TAUGVD_1pt(J,K)+TAUGV_1pt(K)
377
378 100  CONTINUE
379
380       ioptv=1
381
382c HERE END OF THE LOOPS *******
383c******************************
384         
385C TOTAL EXTINCTION OPTICAL DEPTHS
386          DO 119 K=1,NSPECV
387C LOOP OVER NTERMS
388           DO 119 NT=1,NTERM(K)
389           TAUV_1pt(1,K,NT)=0.0
390           TAUVP_1pt(1,K,NT)=0.0
391             DO 119 J=1,NLAYER
392             TAUV_1pt(J+1,K,NT)=TAUV_1pt(J,K,NT)+DTAUV_1pt(J,K,NT)
393             TAUVP_1pt(J+1,K,NT)=TAUVP_1pt(J,K,NT)+DTAUVP_1pt(J,K,NT)
394 119     CONTINUE
395
396
397c       print*,'SETUP'
398c      do i=1,NSPECV
399c      print*,WLNV(i)
400c       do j=1,NLAYER+1
401c       print*,Z(j),TAUV(1,j,i,1),WBARV(1,j,i,1),COSBV(1,j,i,1)
402c       enddo
403c      enddo
404c
405c     IF (IPRINT .GT. 1) THEN
406c           NT=1
407c     IF (2 .GT. 1) THEN
408c          WRITE (6,120)
409c 120      FORMAT(///'  OPTICAL CONSTANTS IN THE VISIBLE (@EQUATOR) ')
410c          WRITE(6,*) 'latitude:',ig
411c          DO 200 K=1,NSPECV
412c          WRITE (6,190)
413c          WRITE (6,210)K,WLNV(K),WNOV(K),BWNV(K)
414c    &    ,BWNV(K)+DWNV(K),DWNV(K)
415c          WRITE (6,230)REALV(K),XIMGV(K)
416c          DO 195 J=1,NLAYER,NLAYER
417C RECALCULATE FOR PRINT OUT ONLY, ONLY FIRST NTERM AT ig=12 (EQUATOR)
418c          WRITE (6,220)XNUMB(J), WBARV_1pt(J,K,NT),COSBV_1pt(J,K,NT)
419c    &      ,DTAUV_1pt(J,K,NT),TAUV_1pt(J,K,NT)
420c 195      CONTINUE
421c          WRITE (6,240) TAUV_1pt(NLEVEL,K,NT)
422c 200      CONTINUE
423c     END IF
424
425c  210 FORMAT(1X,I3,F10.3,F10.2,F10.2,'-',F8.2,F10.3)
426c  190 FORMAT(1X//'  SNUM  MICRONS   WAVENU   INTERVAL    DELTA-WN')
427c  230 FORMAT(1X,'NREAL(LAYER)= ',1PE10.3,' NIMG(LAYER)= ',E10.3/
428c     &' #AEROSOLS   WBAR  COSBAR       DTAU     TAU'
429c     & ,9X,'RAY     GAS    AEROSOL')
430c  220 FORMAT(8(1X,F9.3))
431c  240 FORMAT(41X,F9.3)
432
433      RETURN
434      END
Note: See TracBrowser for help on using the repository browser.