source: trunk/LMDZ.TITAN/libf/phytitan/optcv_1pt.F @ 134

Last change on this file since 134 was 102, checked in by slebonnois, 14 years ago

SL : corrections et modifications dans phytitan correspondant a celles
faites apres compilation Venus. Titan pas encore compile.
File size: 11.9 KB
Line 
1      SUBROUTINE optcv_1pt(zqaer_1pt,ioptv,
2     .            COSBV_1pt,DTAUV_1pt,TAUHV_1pt,TAUHVD_1pt,TAUCV_1pt,
3     .       TAURV_1pt,TAUGV_1pt,TAUGVD_1pt,WBARV_1pt,TAUV_1pt,IPRINT)
4
5
6      use dimphy
7#include "dimensions.h"
8#include "microtab.h"
9#include "clesphys.h"
10
11      PARAMETER(NLAYER=llm,NLEVEL=NLAYER+1)
12      PARAMETER (NSPECI=46,NSPC1I=47,NSPECV=24,NSPC1V=25)
13
14c   Arguments:
15c   ---------
16      integer IPRINT,ioptv
17C ioptv: premier appel, on ne calcule qu'une fois les QM et QF
18* nrad dans microtab.h
19      real   zqaer_1pt(NLAYER,nrad)
20      real   TAUHVD_1pt(NLAYER,NSPECV)
21      real   TAUGVD_1pt(NLAYER,NSPECV)
22      real   TAUHV_1pt(NSPECV),TAUCV_1pt(NSPECV)
23      real   TAURV_1pt(NSPECV),TAUGV_1pt(NSPECV)
24      real   DTAUV_1pt(NLAYER,NSPECV,4),TAUV_1pt(NLEVEL,NSPECV,4)
25      real   WBARV_1pt(NLAYER,NSPECV,4)
26      real   COSBV_1pt(NLAYER,NSPECV,4)
27c   ---------
28
29      COMMON /ATM/ Z(NLEVEL),PRESS(NLEVEL),DEN(NLEVEL),TEMP(NLEVEL)
30
31      COMMON /GASS/ CH4(NLEVEL),XN2(NLEVEL),H2(NLEVEL),AR(NLEVEL)
32     & ,XMU(NLEVEL),GAS1(NLAYER),COLDEN(NLAYER)
33
34      COMMON /VISGAS/SOLARF(NSPECV),NTERM(NSPECV),PEXPON(NSPECV),
35     &         ATERM(4,NSPECV),BTERM(4,NSPECV)
36
37      COMMON /AERSOL/ RADIUS(NLAYER), XNUMB(NLAYER)
38     & , REALI(NSPECI), XIMGI(NSPECI), REALV(NSPECV), XIMGV(NSPECV)
39
40      COMMON /CLOUD/ RADCLD(NLAYER), XNCLD(NLAYER)
41     &             , RCLDI(NSPECI), XICLDI(NSPECI)
42     &             , RCLDV(NSPECV), XICLDV(NSPECV)
43
44      COMMON /SPECTV/ BWNV(NSPC1V),WNOV(NSPECV)
45     &               ,DWNV(NSPECV),WLNV(NSPECV)
46
47      COMMON /PLANT/ CSUBP,RSFI,RSFV,F0PI
48      COMMON /ADJUST/ RHCH4,FH2,FHAZE,FHVIS,FHIR,TAUFAC,RCLOUD,FARGON
49      COMMON /CONST/ RGAS,RHOP,PI,SIGMA
50* nrad dans microtab.h
51      COMMON /part/ v(nrad),rayon(nrad),vrat,dr(nrad),dv(nrad)
52
53      REAL QF1(nrad,NSPECV),QF2(nrad,NSPECV)
54      REAL QF3(nrad,NSPECV),QF4(nrad,NSPECV)
55      REAL QM1(nrad,NSPECV),QM2(nrad,NSPECV)
56      REAL QM3(nrad,NSPECV),QM4(nrad,NSPECV)
57
58      save qf1,qf2,qf3,qf4,qm1,qm2,qm3,qm4
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 OPTFRAC(XMONO,10000./WNOV(K)
152     &                        ,QEXT,QSCT,QABS,QBAR)
153
154c       CALL CFFFV11(1.e-2/WNOV(K),REALV(K),XIMGV(K),RF(inq),2.
155c    &   ,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
174       ENDDO    ! nrad
175
176
177       CBAR=CBAR/TAEROSCAT
178
179        DELTAZ=Z(J)-Z(J+1)
180
181c --------------------------------------------------------------------
182c profil brume Pascal: fit T (sauf tropopause) et albedo
183c -------------------
184        if( cutoff.eq.1) then
185         IF(PRESS(J).gt.9.e-3) THEN
186          TAEROS=TAEROSM1*DELTAZ/DELTAZM1*0.85
187          TAEROSCAT=TAEROSCATM1*DELTAZ/DELTAZM1*0.85
188c         TAEROS=0.
189c         TAEROSCAT=0.
190         ENDIF
191
192         IF(PRESS(J).gt.1.e-1) THEN
193          TAEROS=TAEROSM1*DELTAZ/DELTAZM1*1.15
194          TAEROSCAT=TAEROSCATM1*DELTAZ/DELTAZM1*1.15
195c         TAEROS=0.
196c         TAEROSCAT=0.
197         ENDIF
198        endif !cutoff=1
199
200c profil brume pour fit T (y compris tropopause), mais ne fit plus albedo...
201c -----------------------
202        if( cutoff.eq.2) then
203         IF(PRESS(J).gt.1.e-1) THEN
204          TAEROS=0.
205          TAEROSCAT=0.
206         ENDIF
207        endif !cutoff=2
208c --------------------------------------------------------------------
209
210         TAEROSM1=TAEROS
211         TAEROSCATM1=TAEROSCAT
212         DELTAZM1=DELTAZ
213
214
215       IF (TAEROSCAT.le.0.) CBAR=0.
216
217      if (IPRINT.eq.1) then
218      if (k.eq.NSPECV/2) then   
219       print*,'@VI',K,J,TAEROS,TAEROSCAT,CBAR
220       print*,'@  ',K,J,QF1(1,K),QF2(1,K),zqaer_1pt(NLAYER+1-J,1)
221       print*,'@  ',K,J,QF1(3,K),QF2(3,K),zqaer_1pt(NLAYER+1-J,3)
222       print*,'@  ',K,J,QF1(5,K),QF2(5,K),zqaer_1pt(NLAYER+1-J,5)
223       print*,'@  ',K,J,QF1(7,K),QF2(7,K),zqaer_1pt(NLAYER+1-J,7)
224       print*,'@  ',K,J,QF1(9,K),QF2(9,K),zqaer_1pt(NLAYER+1-J,9)
225       print*
226      endif
227      endif
228
229
230
231c*********** EN TRAVAUX ***************************
232
233C #2:                   RAYLEIGH
234c-------------------------------
235
236C RAYLEIGH SCATTERING STRAIGHT FROM HANSEN AND TRAVIS...SEE NOTES
237C RATIOED BY THE LAYER COLUMN NUMBER TO THE TOTAL
238C COLUMN NUMBER ON EARTH. CM-2
239C THIS IS THE SCATTERING BY THE ATMOSPHERE
240
241      TAURAY=(COLDEN(J)*28.9/(XMU(J)*1013.25))*
242     &(.008569/WLNV(K)**4)*(1.+.0113/WLNV(K)**2+.00013/WLNV(K)**4)
243
244c       PRINT*,WLNV(K)
245c      COLX=0.
246c      COLP=0.
247c      COLT=0.
248c     DO IU=1,NLAYER
249c      COLP=COLDEN(IU)*1.e+1*1.35+COLP
250c     TAURAY=(COLDEN(IU)*28.9/(XMU(IU)*1013.25))*
251c    & (.008569/WLNV(K)**4)*(1.+.0113/WLNV(K)**2
252c    & +.00013/WLNV(K)**4)
253c      COLT=COLT+TAURAY
254c      COLX=COLDEN(IU)*1.e+1/(1.E5*28./22.4E3)*1.e-1*0.0933e-1+COLX
255c                           |   
256c                           |   
257c           g/cm2->kg/m2    |  m2/kg   
258c      Print*,IU, tauray,
259c    &   COLDEN(IU)*1.e+1/(1.E5*28./22.4E3)*1.e-1*0.543e-1
260c     ENDDO
261c       PRINT*,COLP,' PRESSURE AT GROUND;'
262c       PRINT*,COLX,' TAU_GAS AT GROUND;'
263c       print*,colt,colx,' COLT, COLX'
264c      STOP
265                                                       
266c       DZ=Z(J)-Z(J+1)
267c     PRINT*, Z(J),WLNV(K),
268c    &(28.9/(XMU(J)*1013.25))*(.008569/WLNV(K)**4)*
269c    &(1.+.0113/WLNV(K)**2+.00013/WLNV(K)**4)
270c    & ,COLDEN(J)/DZ/100000.,
271c    &(28.9/(XMU(J)*1013.25))*(.008569/WLNV(K)**4)*
272c    &(1.+.0113/WLNV(K)**2+.00013/WLNV(K)**4)
273c    & *COLDEN(J)/DZ/100000.
274   
275 
276
277C #3:                   CLOUD
278c----------------------------
279
280C NEXT COMPUTE TAU CLOUD
281
282      TAUCLD=0.0
283      CBARC =0.0
284      QEXTC =0.0
285      QSCTC =0.0
286c             XNCLD(J)=0.
287      IF ( XNCLD(J) .GT. 0. .and .taufac.gt.0.) THEN
288                CALL XMIE(RADCLD(J),RCLDV(K),XICLDV(K),
289     &                         QEXTC,QSCTC,QABSC,CBARC,WNOV(K))
290                TAUCLD=QEXTC*XNCLD(J)         
291      ENDIF
292C
293      TAURV_1pt(K)=TAURV_1pt(K)+TAURAY
294      TAUGVD_1pt(J,K)=TAURV_1pt(K)
295
296      TAUHV_1pt(K)=TAUHV_1pt(K)+TAEROS          ! INTEGRATED Quant.
297      TAUHVD_1pt(J,K)=TAUHV_1pt(K)
298
299      TAUCV_1pt(K)=TAUCV_1pt(K)+TAUCLD
300
301C #4:                  TAUGAS
302C----------------------------
303
304C LOOP OVER THE NTERMS
305C THIS IS THE ABSORPTION BY THE ATMOSPHERE (METHANE)
306
307
308      DO 909 NT=1,NTERM(K)
309      TAUGAS=COLDEN(J)*GAS1(J)*BTERM(NT,K)*
310     &  (   (PRESS(J+1) + PRESS(J))*.5  )**PEXPON(K)
311
312
313C COMPUTE THE AVERAGE COSBAR AND WBAR
314C&&
315
316c     CBAR=MIN(1.0,1.05*CBAR)       ! THE HAZE FORWARD SCATTERING 5%(WHY?)
317      COSBV_1pt(J,K,NT)=(CBAR*TAEROSCAT + CBARC*TAUCLD)
318     &  /(TAEROSCAT+TAUCLD+TAURAY)     !CBAR_RAY=0.
319c        print*,'CBV',J,K,NT,CBAR,TAEROSCAT,CBARC,TAUCLD       
320
321      DTAUV_1pt(J,K,NT)=TAUGAS+TAEROS+TAURAY+TAUCLD       !TOTAL TAU_EXT
322      TAUGV_1pt(K)=TAUGV_1pt(K)+TAUGAS*ATERM(NT,K)         !TAU_ABS_METH INTEG.
323
324C WE LET W RAYLEIGH BE .999 OR W=1 WHEN ONLY RAYLEIGH=PROBLEM FOR TRID
325c WE HAVE ASSUMED ABOVE THAT COSBAR FOR RAYLEIGH IS ZERO.
326      if (IPRINT.eq.1) then
327      if (k.eq.NSPECV/2) then   
328       print*,'@VI',K,J,DTAUV_1pt(J,K,1),TAUGAS,TAEROS,TAUCLD
329      endif
330      endif
331
332
333c***************** ECHANGE
334c     WBARV(J,K,NT)=(QSCT*XNUMB(J)+TAURAY*0.9999999 + QSCTC*XNCLD(J) )
335c****************
336      WBARV_1pt(J,K,NT)=(TAEROSCAT+TAURAY*0.9999999 + QSCTC*XNCLD(J) )
337c     WBARV_1pt(J,K,NT)=(TAEROSCAT+TAURAY*0.9999999 )
338     &            /(TAUGAS+TAEROS+TAURAY+TAUCLD)
339c****************
340      IF((TAEROS+TAUCLD+TAURAY+TAUCLD).le.0.) WBARV_1pt(J,K,NT)=0.
341      IF((TAEROS+TAUCLD+TAURAY).le.0.) COSBV_1pt(J,K,NT)=0.
342
343c     print*,'WBV',J,K,NT,TAEROSCAT,TAURAY,QSCTC*XNCLD(J)
344c     print*,'WBV',J,K,NT,TAEROS,TAUGAS,TAURAY,TAUCLD
345c     print*,Z(j),J,K,NT,TAUV(1,j,K,NT),WBARV(1,j,K,NT),COSBV(1,j,K,NT)
346
347 909  CONTINUE
348      TAUGVD_1pt(J,K)=TAUGVD_1pt(J,K)+TAUGV_1pt(K)
349 100  CONTINUE
350       ioptv=1
351
352c HERE END OF THE LOOPS *******
353c******************************
354         
355C TOTAL EXTINCTION OPTICAL DEPTHS
356          DO 119 K=1,NSPECV
357C LOOP OVER NTERMS
358           DO 119 NT=1,NTERM(K)
359           TAUV_1pt(1,K,NT)=0.0
360             DO 119 J=1,NLAYER
361             TAUV_1pt(J+1,K,NT)=TAUV_1pt(J,K,NT)+DTAUV_1pt(J,K,NT)
362 119     CONTINUE
363
364c       print*,'SETUP'
365c      do i=1,NSPECV
366c      print*,WLNV(i)
367c       do j=1,NLAYER+1
368c       print*,Z(j),TAUV(1,j,i,1),WBARV(1,j,i,1),COSBV(1,j,i,1)
369c       enddo
370c      enddo
371c
372c     IF (IPRINT .GT. 1) THEN
373c           NT=1
374c     IF (2 .GT. 1) THEN
375c          WRITE (6,120)
376c 120      FORMAT(///'  OPTICAL CONSTANTS IN THE VISIBLE (@EQUATOR) ')
377c          WRITE(6,*) 'latitude:',ig
378c          DO 200 K=1,NSPECV
379c          WRITE (6,190)
380c          WRITE (6,210)K,WLNV(K),WNOV(K),BWNV(K)
381c    &    ,BWNV(K)+DWNV(K),DWNV(K)
382c          WRITE (6,230)REALV(K),XIMGV(K)
383c          DO 195 J=1,NLAYER,NLAYER
384C RECALCULATE FOR PRINT OUT ONLY, ONLY FIRST NTERM AT ig=12 (EQUATOR)
385c          WRITE (6,220)XNUMB(J), WBARV_1pt(J,K,NT),COSBV_1pt(J,K,NT)
386c    &      ,DTAUV_1pt(J,K,NT),TAUV_1pt(J,K,NT)
387c 195      CONTINUE
388c          WRITE (6,240) TAUV_1pt(NLEVEL,K,NT)
389c 200      CONTINUE
390c     END IF
391
392c  210 FORMAT(1X,I3,F10.3,F10.2,F10.2,'-',F8.2,F10.3)
393c  190 FORMAT(1X//'  SNUM  MICRONS   WAVENU   INTERVAL    DELTA-WN')
394c  230 FORMAT(1X,'NREAL(LAYER)= ',1PE10.3,' NIMG(LAYER)= ',E10.3/
395c     &' #AEROSOLS   WBAR  COSBAR       DTAU     TAU'
396c     & ,9X,'RAY     GAS    AEROSOL')
397c  220 FORMAT(8(1X,F9.3))
398c  240 FORMAT(41X,F9.3)
399
400       if (IPRINT.eq.1) stop
401
402      RETURN
403      END
Note: See TracBrowser for help on using the repository browser.