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

Last change on this file since 306 was 306, checked in by slebonnois, 13 years ago

SLebonnois: correction de bugs dans la physique Titan:

  • effg.F : Z doit etre en km, donc conversion
  • optc*_1pt_2.F : On utilise cfffv11 et plus optfrac Du coup, les fichiers input testag* ne sont plus necessaires.
  • phytrac.F : passage de la tendance aerosols en intensif dans tous les cas


File size: 12.0 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
151        CALL CFFFV11(1.e-2/WNOV(K),REALV(K),XIMGV(K),RF(inq),2.
152     &   ,XMONO,QSCT,QEXT,QABS,QBAR)
153
154
155              QF1(inq,K)=QEXT*XRULE
156              QF2(inq,K)=QSCT*XRULE
157              QF3(inq,K)=QABS*XRULE
158              QF4(inq,K)=QBAR
159
160c       print*,'inq=',inq,' QF1=',QF1(inq,K),' QEXT=',QEXT,' XRULE=',XRULE
161                   
162            endif
163
164        TAEROS=QF1(inq,K)*zqaer_1pt(NLAYER+1-J,inq)+TAEROS
165        TAEROSCAT=QF2(inq,K)*zqaer_1pt(NLAYER+1-J,inq)+TAEROSCAT
166        CBAR=CBAR+QF4(inq,K)*QF2(inq,K)*zqaer_1pt(NLAYER+1-J,inq)
167
168           ENDIF
169
170       ENDDO    ! nrad
171
172       IF(TAEROSCAT.le.0.) then
173        CBAR=0.
174       ELSE
175        CBAR=CBAR/TAEROSCAT
176       ENDIF
177
178        DELTAZ=Z(J)-Z(J+1)
179
180c --------------------------------------------------------------------
181c profil brume Pascal: fit T (sauf tropopause) et albedo
182c -------------------
183        if( cutoff.eq.1) then
184         IF(PRESS(J).gt.9.e-3) THEN
185          TAEROS=TAEROSM1*DELTAZ/DELTAZM1*0.85
186          TAEROSCAT=TAEROSCATM1*DELTAZ/DELTAZM1*0.85
187c         TAEROS=0.
188c         TAEROSCAT=0.
189         ENDIF
190
191         IF(PRESS(J).gt.1.e-1) THEN
192          TAEROS=TAEROSM1*DELTAZ/DELTAZM1*1.15
193          TAEROSCAT=TAEROSCATM1*DELTAZ/DELTAZM1*1.15
194c         TAEROS=0.
195c         TAEROSCAT=0.
196         ENDIF
197        endif !cutoff=1
198
199c profil brume pour fit T (y compris tropopause), mais ne fit plus albedo...
200c -----------------------
201        if( cutoff.eq.2) then
202         IF(PRESS(J).gt.1.e-1) THEN
203          TAEROS=0.
204          TAEROSCAT=0.
205         ENDIF
206        endif !cutoff=2
207c --------------------------------------------------------------------
208
209         TAEROSM1=TAEROS
210         TAEROSCATM1=TAEROSCAT
211         DELTAZM1=DELTAZ
212
213
214       IF (TAEROSCAT.le.0.) CBAR=0.
215
2161699  FORMAT(a3,2I3,3(ES15.7,1X))
217
218c*********** EN TRAVAUX ***************************
219
220C #2:                   RAYLEIGH
221c-------------------------------
222
223C RAYLEIGH SCATTERING STRAIGHT FROM HANSEN AND TRAVIS...SEE NOTES
224C RATIOED BY THE LAYER COLUMN NUMBER TO THE TOTAL
225C COLUMN NUMBER ON EARTH. CM-2
226C THIS IS THE SCATTERING BY THE ATMOSPHERE
227
228      TAURAY=(COLDEN(J)*28.9/(XMU(J)*1013.25))*
229     &(.008569/WLNV(K)**4)*(1.+.0113/WLNV(K)**2+.00013/WLNV(K)**4)
230
231
232C #3:                   CLOUD
233c----------------------------
234C NEXT COMPUTE TAU CLOUD
235c
236c  Menu special :
237c  Afin d'eviter la surcharge de calcul on ne calcule les
238c  propriétes optiques des nuages qu'une seule fois
239c  avec un rayon de particule effectif de 3um et une composition
240c  de goutte : 90% CH4 / 10% NOYAUX
241c  Puis on ajute les section efficace par la surface reelle de
242c  la goutte.
243c
244c  ---> A TESTER !!!!
245c
246       IF (clouds.eq.0) THEN
247         CNBAR=0.
248         TNUSCAT=0.
249         TNUABS=0.
250       ELSE
251         IF (ioptv.eq.0.and.j.eq.1) THEN !--> au premier appel
252           QEXTC=0.
253           QSCTC=0.
254           QABSC=0.
255           CBARC=0.
256           DO inq=1,nrad         !BOUCLE SUR LES NQMX TAILLE D"AEROSOLS
257             QC1(inq,K)=0.
258             QC2(inq,K)=0.
259             QC3(inq,K)=0.
260             QC4(inq,K)=0.
261           ENDDO
262** OPTICAL CONSTANT : MIXING RULES
263** Fraction volumique fixe :
264** 10% noyaux.
265** 90% methane.
266           XNR = 0.5 * REALI(K)
267     &         + 0.5 * RCLDI(K)
268           XNI = 0.5 * XIMGI(K)
269     &         + 0.5 * XICLDI(K)
270**
271**   Efficacite : particule de 3um de rayon
272           CALL CMIE(1.E-2/WNOV(K),XNR,XNI,3.e-6,
273     &                QEXTC,QSCTC,QABSC,CBARC)
274
275           DO inq=1,nrad
276             QC1(inq,K)=QEXTC/xnuf
277             QC2(inq,K)=QSCTC/xnuf
278             QC3(inq,K)=QABSC/xnuf
279             QC4(inq,K)=CBARC
280           ENDDO
281         ENDIF   ! ioptv = 0
282         TNUABS=0.
283         TNUSCAT=0.
284         CNBAR=0.
285         IF (rcdb(nlayer+1-J).gt.1.1e-10) THEN
286           DO inq=1,nrad
287             TNUABS=QC1(inq,K)*(rcdb(nlayer+1-J)/3.e-6)**2.*1.e-4*
288     &              zqaer_1pt(NLAYER+1-J,inq+nrad) +
289     &              TNUABS
290             TNUSCAT=QC2(inq,K)*(rcdb(nlayer+1-J)/3.e-6)**2.*1.e-4*
291     &               zqaer_1pt(NLAYER+1-J,inq+nrad) +
292     &               TNUSCAT
293             CNBAR=QC4(inq,K)*QC2(inq,K)*(rcdb(nlayer+1-J)/3.e-6)**2.*
294     &             1.e-4*zqaer_1pt(NLAYER+1-J,inq+nrad) +
295     &             CNBAR
296           ENDDO
297         ENDIF
298
299         IF(TNUSCAT.EQ.0.) THEN
300           CNBAR=0.
301         ELSE
302           CNBAR=CNBAR/TNUSCAT
303         ENDIF
304       ENDIF  ! Cond. CLD
305
306       TAUCV_1pt(K)=TAUCV_1pt(K)+TNUABS
307       TAUCVD_1pt(J,K)=TAUCV_1pt(K)
308
309       TAURV_1pt(K)=TAURV_1pt(K)+TAURAY
310       TAUGVD_1pt(J,K)=TAURV_1pt(K)
311
312       TAUHV_1pt(K)=TAUHV_1pt(K)+TAEROS ! INTEGRATED Quant.
313       TAUHVD_1pt(J,K)=TAUHV_1pt(K)
314
315
316
317C #4:                  TAUGAS
318C----------------------------
319
320C LOOP OVER THE NTERMS
321C THIS IS THE ABSORPTION BY THE ATMOSPHERE (METHANE)
322
323
324       DO 909 NT=1,NTERM(K)
325         TAUGAS=COLDEN(J)*GAS1(J)*BTERM(NT,K)*
326     &   (   (PRESS(J+1) + PRESS(J))*.5  )**PEXPON(K)
327
328
329*  COSBV ET COSBVP
330*-----------------
331
332         IF(TAEROSCAT+TNUSCAT+TAURAY .ne. 0.) THEN
333           COSBV_1pt(J,K,NT)=(CBAR*TAEROSCAT + CNBAR*TNUSCAT)
334     &     /(TAEROSCAT+TNUSCAT+TAURAY) !CBAR_RAY=0.
335         ELSE
336           COSBV_1pt(J,K,NT)=0.
337         ENDIF
338
339         IF(TAEROSCAT+TAURAY .ne. 0.) THEN
340           COSBVP_1pt(J,K,NT)=(CBAR*TAEROSCAT)
341     &     /(TAEROSCAT+TAURAY) !CBAR_RAY=0.
342         ELSE
343           COSBVP_1pt(J,K,NT)=0.
344         ENDIF
345
346*  DTAUV ET DTAUVP
347*-----------------
348
349         DTAUV_1pt(J,K,NT) =TAUGAS+TAEROS+TAURAY+TNUABS !TAU_ABS_METH
350         DTAUVP_1pt(J,K,NT)=TAUGAS+TAEROS+TAURAY       !TAU_ABS_METH
351
352         TAUGV_1pt(K)=TAUGV_1pt(K)+TAUGAS*ATERM(NT,K) !INTEG.
353
354*  WBARV ET WBARVP
355*-----------------
356
357         IF(TAUGAS+TAEROS+TAURAY+TNUABS .ne.  0.) THEN
358           WBARV_1pt(J,K,NT)=(TAEROSCAT+TAURAY*0.9999999 + TNUSCAT)
359     &     /(TAUGAS+TAEROS+TAURAY+TNUABS)
360         ELSE
361           WBARV_1pt(J,K,NT)=0.
362         ENDIF
363
364         IF(TAUGAS+TAEROS+TAURAY .ne.  0.) THEN
365           WBARVP_1pt(J,K,NT)=(TAEROSCAT+TAURAY*0.9999999 )
366     &     /(TAUGAS+TAEROS+TAURAY)
367         ELSE
368           WBARVP_1pt(J,K,NT)=0.
369         ENDIF
370   
371 909   CONTINUE
372 
373       TAUGVD_1pt(J,K)=TAUGVD_1pt(J,K)+TAUGV_1pt(K)
374
375 100  CONTINUE
376
377       ioptv=1
378
379c HERE END OF THE LOOPS *******
380c******************************
381         
382C TOTAL EXTINCTION OPTICAL DEPTHS
383          DO 119 K=1,NSPECV
384C LOOP OVER NTERMS
385           DO 119 NT=1,NTERM(K)
386           TAUV_1pt(1,K,NT)=0.0
387           TAUVP_1pt(1,K,NT)=0.0
388             DO 119 J=1,NLAYER
389             TAUV_1pt(J+1,K,NT)=TAUV_1pt(J,K,NT)+DTAUV_1pt(J,K,NT)
390             TAUVP_1pt(J+1,K,NT)=TAUVP_1pt(J,K,NT)+DTAUVP_1pt(J,K,NT)
391 119     CONTINUE
392
393
394c       print*,'SETUP'
395c      do i=1,NSPECV
396c      print*,WLNV(i)
397c       do j=1,NLAYER+1
398c       print*,Z(j),TAUV(1,j,i,1),WBARV(1,j,i,1),COSBV(1,j,i,1)
399c       enddo
400c      enddo
401c
402c     IF (IPRINT .GT. 1) THEN
403c           NT=1
404c     IF (2 .GT. 1) THEN
405c          WRITE (6,120)
406c 120      FORMAT(///'  OPTICAL CONSTANTS IN THE VISIBLE (@EQUATOR) ')
407c          WRITE(6,*) 'latitude:',ig
408c          DO 200 K=1,NSPECV
409c          WRITE (6,190)
410c          WRITE (6,210)K,WLNV(K),WNOV(K),BWNV(K)
411c    &    ,BWNV(K)+DWNV(K),DWNV(K)
412c          WRITE (6,230)REALV(K),XIMGV(K)
413c          DO 195 J=1,NLAYER,NLAYER
414C RECALCULATE FOR PRINT OUT ONLY, ONLY FIRST NTERM AT ig=12 (EQUATOR)
415c          WRITE (6,220)XNUMB(J), WBARV_1pt(J,K,NT),COSBV_1pt(J,K,NT)
416c    &      ,DTAUV_1pt(J,K,NT),TAUV_1pt(J,K,NT)
417c 195      CONTINUE
418c          WRITE (6,240) TAUV_1pt(NLEVEL,K,NT)
419c 200      CONTINUE
420c     END IF
421
422c  210 FORMAT(1X,I3,F10.3,F10.2,F10.2,'-',F8.2,F10.3)
423c  190 FORMAT(1X//'  SNUM  MICRONS   WAVENU   INTERVAL    DELTA-WN')
424c  230 FORMAT(1X,'NREAL(LAYER)= ',1PE10.3,' NIMG(LAYER)= ',E10.3/
425c     &' #AEROSOLS   WBAR  COSBAR       DTAU     TAU'
426c     & ,9X,'RAY     GAS    AEROSOL')
427c  220 FORMAT(8(1X,F9.3))
428c  240 FORMAT(41X,F9.3)
429
430      RETURN
431      END
Note: See TracBrowser for help on using the repository browser.