source: trunk/LMDZ.TITAN/libf/phytitan/optcv_1pt_3.F @ 1243

Last change on this file since 1243 was 1126, checked in by slebonnois, 11 years ago

SL: update of Titan photochemical module to include computation of chemistry up to 1300 km
File size: 11.6 KB
RevLine 
[808]1      SUBROUTINE optcv_1pt3(zqaer_1pt,rcdb,xfrb,ioptv,IPRINT)
2
3      use dimphy
4      IMPLICIT NONE
5#include "dimensions.h"
6#include "microtab.h"
7#include "clesphys.h"
8
9      integer nlayer, nlevel, nspeci, nspc1i, nspecv, nspc1v,nterm
10      real z,press, den, temp, ch4, xn2, h2, ar, xmu, gas1,
11     &     colden, c2h2, c2h6, hcn, radius, xnumb, reali
12     &     ximgi, realv, ximgv, rcldi, xicldi, rcldv, xicldv, rcldi2,
13     &     xicldi2, rcldv2, xicldv2,real bwni, wnoi, dwni, wlni, csubp,
[1056]14     &     f0pi, rhch4, fh2, fhaze, fhvis
[808]15     &     reali, ximgi, bwni, fhir, taufac, rcloud, fargon, rgas, rhop,
16     &     pi, sigma, prod,reali,fhvis
17
18      integer k, j,inq,nt
19
20      real tbar, pbar, bmu, coef1, effg, taeros, taeroscat, cbar,
21     &     qext, qsct, qabs, qbar, xmono, xrule, deltaz,
22     &     cnbar, qextc, qsctc, qabsc, qbarc, taugas, pnn,
23     &     pcc, pcn, phn, kgas, u, ig, tau2, tlimit,
24     &     solarf, pexpon, aterm, bterm, bwnv, wnov, dwnv,
25     &     wlnv, v, rayon, vrat, dr, dv, taerosm1, deltazm1,
26     &     taeroscatm1, tauray
27
28      PARAMETER(NLAYER=llm,NLEVEL=NLAYER+1)
29      PARAMETER (NSPECI=46,NSPC1I=47,NSPECV=24,NSPC1V=25)
30
31c   Arguments:
32c   ---------
33      integer IPRINT,ioptv
34C ioptv: premier appel, on ne calcule qu'une fois les QM et QF
35* nrad dans microtab.h
36      real   zqaer_1pt(NLAYER,2*nrad)
37#include "optcv_1pt.h"
38c   ---------
39
40      COMMON /ATM/ Z(NLEVEL),PRESS(NLEVEL),DEN(NLEVEL),TEMP(NLEVEL)
41
42      COMMON /GASS/ CH4(NLEVEL),XN2(NLEVEL),H2(NLEVEL),AR(NLEVEL)
43     & ,XMU(NLEVEL),GAS1(NLAYER),COLDEN(NLAYER)
44
45      COMMON /VISGAS/SOLARF(NSPECV),NTERM(NSPECV),PEXPON(NSPECV),
46     &         ATERM(4,NSPECV),BTERM(4,NSPECV)
47
48      COMMON /AERSOL/ RADIUS(NLAYER), XNUMB(NLAYER)
49     & , REALI(NSPECI), XIMGI(NSPECI), REALV(NSPECV), XIMGV(NSPECV)
50
51      COMMON /CLOUD/
52     &               RCLDI(NSPECI), XICLDI(NSPECI)
53     &             , RCLDV(NSPECV), XICLDV(NSPECV)
54     &             , RCLDI2(NSPECI), XICLDI2(NSPECI)
55     &             , RCLDV2(NSPECV), XICLDV2(NSPECV)
56
57      COMMON /SPECTV/ BWNV(NSPC1V),WNOV(NSPECV)
58     &               ,DWNV(NSPECV),WLNV(NSPECV)
59
[1056]60      COMMON /PLANT/ CSUBP,F0PI
[808]61      COMMON /ADJUST/ RHCH4,FH2,FHAZE,FHVIS,FHIR,TAUFAC,RCLOUD,FARGON
62      COMMON /CONST/ RGAS,RHOP,PI,SIGMA
63* nrad dans microtab.h
64      COMMON /part/ v(nrad),rayon(nrad),vrat,dr(nrad),dv(nrad)
65
66      REAL QF1(nrad,NSPECV),QF2(nrad,NSPECV)
67      REAL QF3(nrad,NSPECV),QF4(nrad,NSPECV)
68      REAL QM1(nrad,NSPECV),QM2(nrad,NSPECV)
69      REAL QM3(nrad,NSPECV),QM4(nrad,NSPECV)
70
71c---- NUAGES
72      real TNUEXT,TNUSCAT
73      real   rcdb(NLAYER), xfrb(NLAYER,4)
74
[814]75      save qf1,qf2,qf3,qf4,qm1,qm2,qm3,qm4
[808]76
77
78C*
79C THIS SUBROUTINE SETS THE OPTICAL CONSTANTS IN THE VISIBLE
80C IT CALCUALTES FOR EACH LAYER, FOR EACH SPECRAL INTERVAL IN THE VIS
81C LAYER: WBAR, DTAU, COSBAR
82C LEVEL: TAU
83C
84C ZERO THE COLUMN OPTICAL DEPTHS OF EACH TYPE
85C ??FLAG? THE OPTICAL DEPTH OF THE TOP OF THE MODEL
86C MAY NOT BE ZERO.
87
88c******* DEBUT DES BOUCLES ************************
89      DO 100 K=1,NSPECV         !b! BOUCLE SUR LAMBDA
90
91      TAURV_1pt(K)=0.
92      TAUHV_1pt(K)=0.            ! INTEGRATED TAU.......INITIALIZATION.
93      TAUCV_1pt(K)=0.            ! Rayleigh, Haze, Cloud, Gas
94      TAUGV_1pt(K)=0.            !   sca,    abs,  abs  , abs
95
96      DO 100 J=1,NLAYER         !a! BOUCLE SUR L"ALTITUDE
97
98C #1:                   HAZE
99c---------------------------
100
101c     CALL THE MIE CODE TO GIVE THE AEROSOL PROPERTIES
102c     USE XFRAC FOR FRACTAL AEROSOLS PROPERTIES AT LAMBDA < 2. um
103
104
105
106
107c                    /\
108c                   /  \
109c                  /    \
110c                 / _O   \
111c                / |/     \
112c               /  / \     \
113c              /   |\ \/\   \
114c             /    || /  \   \
115c             ----------------
116c            |     WARNING    |
117c            |    SLOW DOWN   |
118c             ----------------
119
120
121
122
123c*********** EN TRAVAUX ***************************
124 
125         TAEROS=0.
126         TAEROSCAT=0.
127         CBAR=0.
128
129c       print*,"rayon=",rayon
130c       print*,"RF=",RF
131
132      DO inq=1,nrad         !BOUCLE SUR LES TAILLE D"AEROSOLS
133
134
135            IF (rayon(inq).lt.RF(inq)) THEN    ! aerosols spheriques
136
137             
138            if(ioptv.eq.0.and.J.eq.1) then
139c                  CALL XMIE(rayon(inq)*1.e6,REALV(K),XIMGV(K),
140c    &             QEXT,QSCT,QABS,QBAR,WNOV(K))
141
142                CALL CMIE(1.E-2/WNOV(K),REALV(K),XIMGV(K),rayon(inq),
143     &          QEXT,QSCT,QABS,QBAR)
144
145c       print*,'inq=',inq,' QM1=',QM1(inq,K),' QEXT=',QEXT
146
147              QM1(inq,K)=QEXT
148              QM2(inq,K)=QSCT
149              QM3(inq,K)=QABS
150              QM4(inq,K)=QBAR
151            endif
152
153       TAEROS=QM1(inq,K)*zqaer_1pt(NLAYER+1-J,inq)*1.e-4+TAEROS
154       TAEROSCAT=QM2(inq,K)*zqaer_1pt(NLAYER+1-J,inq)*1.e-4+TAEROSCAT
155       CBAR=CBAR+QM4(inq,K)*QM2(inq,K)*zqaer_1pt(NLAYER+1-J,inq)*1.e-4
156
157            ELSE                        ! aerosols fractals
158
159               XMONO=(rayon(inq)/RF(inq))**3.
160               XRULE=1.
161
162            if(XMONO.gt.16384./1.5) then
163             XRULE=(XMONO/16384.)
164             XMONO=16384.
165            endif
166
167            if(ioptv.eq.0.and.J.eq.1) then
168
169c        CALL OPTFRAC(XMONO,10000./WNOV(K)
170c     &                        ,QEXT,QSCT,QABS,QBAR)
171
172        CALL CFFFV11(1.e-2/WNOV(K),REALV(K),XIMGV(K),RF(inq),2.
173     &   ,XMONO,QSCT,QEXT,QABS,QBAR)
174
175
176              QF1(inq,K)=QEXT*XRULE
177              QF2(inq,K)=QSCT*XRULE
178              QF3(inq,K)=QABS*XRULE
179              QF4(inq,K)=QBAR
180
181c       print*,'inq=',inq,' QF1=',QF1(inq,K),' QEXT=',QEXT,' XRULE=',XRULE
182                   
183            endif
184
185        TAEROS=QF1(inq,K)*zqaer_1pt(NLAYER+1-J,inq)+TAEROS
186        TAEROSCAT=QF2(inq,K)*zqaer_1pt(NLAYER+1-J,inq)+TAEROSCAT
187        CBAR=CBAR+QF4(inq,K)*QF2(inq,K)*zqaer_1pt(NLAYER+1-J,inq)
188
189           ENDIF
190
191       ENDDO    ! nrad
192
193
[1058]194       if (TAEROSCAT.ne.0.) CBAR=CBAR/TAEROSCAT
[808]195
196        DELTAZ=Z(J)-Z(J+1)
197
198c --------------------------------------------------------------------
199c profil brume Pascal: fit T (sauf tropopause) et albedo
200c -------------------
201        if( cutoff.eq.1) then
202         IF(PRESS(J).gt.9.e-3) THEN
203          TAEROS=TAEROSM1*DELTAZ/DELTAZM1*0.85
204          TAEROSCAT=TAEROSCATM1*DELTAZ/DELTAZM1*0.85
205c         TAEROS=0.
206c         TAEROSCAT=0.
207         ENDIF
208
209         IF(PRESS(J).gt.1.e-1) THEN
210          TAEROS=TAEROSM1*DELTAZ/DELTAZM1*1.15
211          TAEROSCAT=TAEROSCATM1*DELTAZ/DELTAZM1*1.15
212c         TAEROS=0.
213c         TAEROSCAT=0.
214         ENDIF
215        endif !cutoff=1
216
217c profil brume pour fit T (y compris tropopause), mais ne fit plus albedo...
218c -----------------------
219        if( cutoff.eq.2) then
220         IF(PRESS(J).gt.1.e-1) THEN
221          TAEROS=0.
222          TAEROSCAT=0.
223         ENDIF
224        endif !cutoff=2
225c --------------------------------------------------------------------
226
227         TAEROSM1=TAEROS
228         TAEROSCATM1=TAEROSCAT
229         DELTAZM1=DELTAZ
230
231
232       IF (TAEROSCAT.le.0.) CBAR=0.
233
[1126]234c     print*, 'HERE, CIRS AEROSOLS'
235c     call cirs_haze(PRESS(J),WNOV(K),TAEROS,TAEROSCAT,CBAR)
236
[808]2371699  FORMAT(a3,2I3,3(ES15.7,1X))
238
239c*********** EN TRAVAUX ***************************
240
241C #2:                   RAYLEIGH
242c-------------------------------
243
244C RAYLEIGH SCATTERING STRAIGHT FROM HANSEN AND TRAVIS...SEE NOTES
245C RATIOED BY THE LAYER COLUMN NUMBER TO THE TOTAL
246C COLUMN NUMBER ON EARTH. CM-2
247C THIS IS THE SCATTERING BY THE ATMOSPHERE
248
249      TAURAY=(COLDEN(J)*28.9/(XMU(J)*1013.25))*
250     &(.008569/WLNV(K)**4)*(1.+.0113/WLNV(K)**2+.00013/WLNV(K)**4)
251
252
253C #3:                   CLOUD
254c----------------------------
255C NEXT COMPUTE TAU CLOUD
256c
257c  Menu special :
258c  On utilise ici une look-up table afin de calculer
259c  les proprietes optique des nuages.
260c  Le principe est le suivant :
261c  La look-up table contient les proprietes optique d'une goutte
262c  de methane pur de 3 um.
263c  On approxime les proprietes optiques pour une goutte de rayon r a
264c  de la table.
265c
266       TNUEXT=0.
267       TNUSCAT=0.
268       CNBAR=0.
269       IF (clouds.eq.1) THEN
270
271         CALL getoptcld(1.E-2/WNOV(K),rcdb(nlayer+1-J),
272     &                  QEXTC,QSCTC,QABSC,QBARC)
273         TNUEXT=0.
274         TNUSCAT=0.
275         CNBAR=0.
276         IF (rcdb(nlayer+1-J).gt.1.1e-10) THEN
277           TNUEXT =QEXTC/xnuf*SUM(zqaer_1pt(NLAYER+1-J,nrad+1:2*nrad))
278           TNUSCAT=QSCTC/xnuf*SUM(zqaer_1pt(NLAYER+1-J,nrad+1:2*nrad))
279           CNBAR  =QBARC
280         ENDIF
281           IF(TNUSCAT.GE.0.8*TNUEXT) TNUSCAT=0.8*TNUEXT
282       ENDIF  ! Cond. CLD
283
284       TAUCV_1pt(K)=TAUCV_1pt(K)+TNUEXT
285       TAUCVD_1pt(J,K)=TAUCV_1pt(K)
286
287       TAURV_1pt(K)=TAURV_1pt(K)+TAURAY
288       TAUGVD_1pt(J,K)=TAURV_1pt(K)
289
290       TAUHV_1pt(K)=TAUHV_1pt(K)+TAEROS ! INTEGRATED Quant.
291       TAUHVD_1pt(J,K)=TAUHV_1pt(K)
292
293
294
295C #4:                  TAUGAS
296C----------------------------
297
298C LOOP OVER THE NTERMS
299C THIS IS THE ABSORPTION BY THE ATMOSPHERE (METHANE)
300
301
302       DO 909 NT=1,NTERM(K)
303         TAUGAS=COLDEN(J)*GAS1(J)*BTERM(NT,K)*
304     &   (   (PRESS(J+1) + PRESS(J))*.5  )**PEXPON(K)
305
306
307*  COSBV ET COSBVP
308*-----------------
309
310         IF(TAEROSCAT+TNUSCAT+TAURAY .ne. 0.) THEN
311           COSBV_1pt(J,K,NT)=(CBAR*TAEROSCAT + CNBAR*TNUSCAT)
312     &     /(TAEROSCAT+TNUSCAT+TAURAY) !CBAR_RAY=0.
313         ELSE
314           COSBV_1pt(J,K,NT)=0.
315         ENDIF
316
317         IF(TAEROSCAT+TAURAY .ne. 0.) THEN
318           COSBVP_1pt(J,K,NT)=(CBAR*TAEROSCAT)
319     &     /(TAEROSCAT+TAURAY) !CBAR_RAY=0.
320         ELSE
321           COSBVP_1pt(J,K,NT)=0.
322         ENDIF
323
324*  DTAUV ET DTAUVP
325*-----------------
326
327         DTAUV_1pt(J,K,NT) =TAUGAS+TAEROS+TAURAY+TNUEXT !TAU_ABS_METH
328         DTAUVP_1pt(J,K,NT)=TAUGAS+TAEROS+TAURAY       !TAU_ABS_METH
329
330         TAUGV_1pt(K)=TAUGV_1pt(K)+TAUGAS*ATERM(NT,K) !INTEG.
331
332*  WBARV ET WBARVP
333*-----------------
334
335         IF(TAUGAS+TAEROS+TAURAY+TNUEXT .ne.  0.) THEN
336           WBARV_1pt(J,K,NT)=(TAEROSCAT+TAURAY*0.9999999 + TNUSCAT)
337     &     /(TAUGAS+TAEROS+TAURAY+TNUEXT)
338         ELSE
339           WBARV_1pt(J,K,NT)=0.
340         ENDIF
341
342         IF(TAUGAS+TAEROS+TAURAY .ne.  0.) THEN
343           WBARVP_1pt(J,K,NT)=(TAEROSCAT+TAURAY*0.9999999 )
344     &     /(TAUGAS+TAEROS+TAURAY)
345         ELSE
346           WBARVP_1pt(J,K,NT)=0.
347         ENDIF
348   
349 909   CONTINUE
350 
351       TAUGVD_1pt(J,K)=TAUGVD_1pt(J,K)+TAUGV_1pt(K)
352
353 100  CONTINUE
354
355       ioptv=1
356
357c HERE END OF THE LOOPS *******
358c******************************
359         
360C TOTAL EXTINCTION OPTICAL DEPTHS
361          DO 119 K=1,NSPECV
362C LOOP OVER NTERMS
363           DO 119 NT=1,NTERM(K)
364           TAUV_1pt(1,K,NT)=0.0
365           TAUVP_1pt(1,K,NT)=0.0
366             DO 119 J=1,NLAYER
367             TAUV_1pt(J+1,K,NT)=TAUV_1pt(J,K,NT)+DTAUV_1pt(J,K,NT)
368             TAUVP_1pt(J+1,K,NT)=TAUVP_1pt(J,K,NT)+DTAUVP_1pt(J,K,NT)
369 119     CONTINUE
370
371
372c       print*,'SETUP'
373c      do i=1,NSPECV
374c      print*,WLNV(i)
375c       do j=1,NLAYER+1
376c       print*,Z(j),TAUV(1,j,i,1),WBARV(1,j,i,1),COSBV(1,j,i,1)
377c       enddo
378c      enddo
379c
380c     IF (IPRINT .GT. 1) THEN
381c           NT=1
382c     IF (2 .GT. 1) THEN
383c          WRITE (6,120)
384c 120      FORMAT(///'  OPTICAL CONSTANTS IN THE VISIBLE (@EQUATOR) ')
385c          DO 200 K=1,NSPECV
386c          WRITE (6,190)
387c          WRITE (6,210)K,WLNV(K),WNOV(K),BWNV(K)
388c    &    ,BWNV(K)+DWNV(K),DWNV(K)
389c          WRITE (6,230)REALV(K),XIMGV(K)
390c          DO 195 J=1,NLAYER,NLAYER
391c          WRITE (6,220)XNUMB(J), WBARV_1pt(J,K,NT),COSBV_1pt(J,K,NT)
392c    &      ,DTAUV_1pt(J,K,NT),TAUV_1pt(J,K,NT)
393c 195      CONTINUE
394c          WRITE (6,240) TAUV_1pt(NLEVEL,K,NT)
395c 200      CONTINUE
396c     END IF
397
398c  210 FORMAT(1X,I3,F10.3,F10.2,F10.2,'-',F8.2,F10.3)
399c  190 FORMAT(1X//'  SNUM  MICRONS   WAVENU   INTERVAL    DELTA-WN')
400c  230 FORMAT(1X,'NREAL(LAYER)= ',1PE10.3,' NIMG(LAYER)= ',E10.3/
401c     &' #AEROSOLS   WBAR  COSBAR       DTAU     TAU'
402c     & ,9X,'RAY     GAS    AEROSOL')
403c  220 FORMAT(8(1X,F9.3))
404c  240 FORMAT(41X,F9.3)
405
406      RETURN
407      END
Note: See TracBrowser for help on using the repository browser.