source: trunk/libf/phytitan/optci.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: 15.5 KB
Line 
1      SUBROUTINE OPTCI(ykim,nmicro,IPRINT)
2#include "dimensions.h"
3#include "dimphy.h"
4#include "microtab.h"
5#include "numchimrad.h"
6#include "clesphys.h"
7
8c   Arguments:
9c   ---------
10      REAL    ykim(klon,klev,nqmx)
11      integer nmicro
12c   ---------
13
14      PARAMETER(NLAYER=llm,NLEVEL=NLAYER+1)
15      PARAMETER (NSPECI=46,NSPC1I=47,NSPECV=24,NSPC1V=25)
16
17      COMMON /ATM/ Z(NLEVEL),PRESS(NLEVEL),DEN(NLEVEL),TEMP(NLEVEL)
18
19      COMMON /GASS/ CH4(NLEVEL),XN2(NLEVEL),H2(NLEVEL),AR(NLEVEL)
20     & ,XMU(NLEVEL),GAS1(NLAYER),COLDEN(NLAYER)
21
22      COMMON /STRATO/ C2H2(NLAYER),C2H6(NLAYER)
23      COMMON /STRAT2/ HCN(NLAYER)
24
25      COMMON /AERSOL/ RADIUS(NLAYER), XNUMB(NLAYER)
26     & , REALI(NSPECI), XIMGI(NSPECI), REALV(NSPECV), XIMGV(NSPECV)
27
28      COMMON /CLOUD/ RADCLD(NLAYER), XNCLD(NLAYER)
29     & , RCLDI(NSPECI), XICLDI(NSPECI), 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
42      COMMON /OPTICI/ DTAUI(klon,NLAYER,NSPECI)
43     &               ,TAUI (klon,NLEVEL,NSPECI)
44     &               ,WBARI(klon,NLAYER,NSPECI)
45     &               ,COSBI(klon,NLAYER,NSPECI)
46
47      COMMON /SPECTI/ BWNI(NSPC1I), WNOI(NSPECI),
48     &                DWNI(NSPECI), WLNI(NSPECI)
49
50      COMMON /PLANT/ CSUBP,RSFI,RSFV,F0PI
51      COMMON /ADJUST/ RHCH4,FH2,FHAZE,FHVIS,FHIR,TAUFAC,RCLOUD,FARGON
52      COMMON /CONST/RGAS,RHOP,PI,SIGMA
53      COMMON /traceurs/qaer
54      COMMON /part/v,r,vrat,dr,dv
55
56      DIMENSION PROD(NLEVEL)
57      real qaer(klon,nlayer,nqmx)
58      real v(nqmx),r(nqmx),vrat,dr(nqmx),dv(nqmx)
59      real xv1(klev,nspeci),xv2(klev,nspeci)
60      real xv3(klev,nspeci)
61      REAL QF1(nqmx,NSPECI),QF2(nqmx,NSPECI)
62      REAL QF3(nqmx,NSPECI),QF4(nqmx,NSPECI)
63      REAL QM1(nqmx,NSPECI),QM2(nqmx,NSPECI)
64      REAL QM3(nqmx,NSPECI),QM4(nqmx,NSPECI)
65      real emu
66      REAL TAEROSM1(NSPECI),TAEROSCATM1(NSPECI),DELTAZM1(NSPECI)
67     
68      save qf1,qf2,qf3,qf4,qm1,qm2,qm3,qm4
69
70      integer iopti,iwarning
71      integer ig_
72      save iopti,iwarning
73      data iopti,iwarning/0,0/
74
75C THE PRESSURE INDUCED TRANSITIONS ARE FROM REGIS
76C THE LAST SEVENTEEN INTERVALS ARE THE BANDS FROM GNF.
77C
78C THIS SUBROUTINE SETS THE OPTICAL CONSTANTS IN THE INFRARED
79C IT CALCUALTES FOR EACH LAYER, FOR EACH SPECRAL INTERVAL IN THE IR
80C LAYER: WBAR, DTAU, COSBAR
81C LEVEL: TAU
82C
83       print*,'START OPTCI'
84c      sum=0.
85c      do nng=2,klon
86c        do i=1,klev
87c         do j=1,nqmx
88c          print*,'j,rj',j,r(j)
89c          print*,'paer',qaer(nng,i,j)
90c           sum=sum+qaer(nng,i,j)*r(j)**3.*1.3333*3.1415*1000.
91c         enddo
92c        enddo
93c        enddo
94c      print*,sum/(klon-1),'SOMME COLONNE/OPTCI'
95
96
97c      do inq=1,nqmx
98c          print*,inq,r(inq),vrat,qaer(12,25,inq)
99c      enddo
100             
101      DO 420 K=1,NSPECI
102C LETS USE THE THOLIN OPTICAL CONSTANTS FOR THE HAZE.
103          CALL THOLIN(WLNI(K),TNR,TNI)
104          REALI(K)=TNR
105          XIMGI(K)=TNI*FHIR
106C SET UP THE OPTICAL CONSTANTS FOR THE CLOUD
107          RCLDI(K)=1.27
108          XICLDI(K)=REFLIQ(WNOI(K))
109 420  CONTINUE
110C
111C
112C ZERO ALL OPTICAL DEPTHS.
113C ??FLAG? FOR SOME APPLCIATIONS THE TOP OPACITY MAY NOT VANISH
114
115c      open  (unit=1,file='xsetupi')
116c      do i=1,klev
117c       read(1,*) a
118c        do j=1,nspeci
119c            read(1,*) xv1(i,j),xv2(i,j),xv3(i,j)
120c        enddo
121c       enddo
122c       close(1)
123
124c************************************************************************
125c************************************************************************
126      DO 79 ig=1,klon      ! BOUCLE SUR GRILLE HORIZONTALE   
127c      print*,'ig NEW optci',ig
128c************************************************************************
129c************************************************************************
130
131        DO 80 K=1,NSPECI
132           TAUHI(ig,K)=0.
133           TAUCI(ig,K)=0.
134           TAUGI(ig,K)=0.
135 80     CONTINUE
136
137        if (.not.ylellouch) then
138       
139            XN2(1) = ykim(ig,1,iradn2)
140            CH4(1) = ykim(ig,1,iradch4)
141             H2(1) = ykim(ig,1,iradh2)
142            do j=2,nlayer
143               XN2(j) = (ykim(ig,j,iradn2)+ykim(ig,j-1,iradn2))/2.
144               CH4(j) = (ykim(ig,j,iradch4)+ykim(ig,j-1,iradch4))/2.
145                H2(j) = (ykim(ig,j,iradh2)+ykim(ig,j-1,iradh2))/2.
146            enddo
147            XN2(nlevel) = ykim(ig,nlayer,iradn2)
148            CH4(nlevel) = ykim(ig,nlayer,iradch4)
149             H2(nlevel) = ykim(ig,nlayer,iradh2)     
150
151            do j=1,nlayer
152               emu = ( xmu(j) + xmu(j+1) )/2.
153               C2H2(j) = ykim(ig,j,iradc2h2) * 26./emu
154               C2H6(j) = ykim(ig,j,iradc2h6) * 30./emu
155                HCN(j) = ykim(ig,j,iradhcn ) * 27./emu
156            enddo
157                 
158        endif
159
160c************************************************************************
161          DO 100 J=1,NLAYER    ! BOUCLE SUR L'ALTITUDE
162c************************************************************************
163
164c      print*,'ig,k,j ',ig,k,j
165
166C SET UP THE COEFFICIENT TO REDUCE MASS PATH TO STP ...SEE NOTES
167C T0 =273.15   PO=1.01325 BAR
168
169        TBAR=0.5*(TEMP(J)+TEMP(J+1))
170        PBAR=SQRT(PRESS(J)*PRESS(J+1))
171        BMU=0.5*(XMU(J+1)+XMU(J))
172        COEF1=RGAS*273.15**2*.5E5* (PRESS(J+1)**2 - PRESS(J)**2)
173     & /(1.01325**2 *EFFG(Z(J))*TBAR*BMU)
174
175      IF (IPRINT .GT. 9) WRITE(6,21) J,EFFG(Z(J)),TBAR,BMU,COEF1
176 21   FORMAT(' J, EFFG, TBAR, BMU, COEF1,: ',I3,1P6E10.3)
177
178c------------------------------------------------------------------------
179         DO 101 K=1,NSPECI     ! BOUCLE SUR LES L.D'O
180c------------------------------------------------------------------------
181
182
183C #1:             HAZE
184C---------------------
185
186C FIRST COMPUTE TAU AEROSOL
187
188
189c
190c                    /\
191c                   /  \
192c                  /    \
193c                 / _O   \
194c                / |/     \
195c               /  / \     \
196c              /   |\ \/\   \
197c             /    || /  \   \
198c             ----------------
199c            |     WARNING    |
200c            |    SLOW DOWN   |
201c             ----------------
202
203
204
205
206c*********** EN TRAVAUX ***************************
207
208         TAEROS=0.
209         TAEROSCAT=0.
210         CBAR=0.
211
212      DO inq=1,nmicro         !BOUCLE SUR LES NMICRO TAILLE D"AEROSOLS
213
214
215      IF (WNOI(K).lt.wco) THEN    ! lamda > 56 um
216
217           if (iopti.eq.0) then
218
219c          CALL XMIE(R(inq)*1.e6,REALI(K),XIMGI(K),
220c    &     QEXT,QSCT,QABS,QBAR,WNOI(K))
221
222
223           CALL CMIE(1.E-2/WNOI(K),REALI(K),XIMGI(K),R(inq),
224     &     QEXT,QSCT,QABS,QBAR)
225
226
227           QM1(inq,K)=QEXT
228           QM2(inq,K)=QSCT
229           QM3(inq,K)=QABS
230           QM4(inq,K)=QBAR
231
232          endif         ! end iopti
233
234
235      if (microfi.eq.1) then
236         ig_=ig
237      else
238         ig_=12
239      endif
240
241
242      TAEROS=QM1(inq,K)*qaer(ig_,nlayer+1-J,inq)*1.e-4+TAEROS
243      TAEROSCAT=QM2(inq,K)*qaer(ig_,nlayer+1-J,inq)*1.e-4+TAEROSCAT
244      CBAR=CBAR+QM4(inq,K)*QM2(inq,K)*qaer(ig_,nlayer+1-J,inq)*1.e-4
245
246
247         ELSE                           ! 0.2 < lambda < 56 um
248
249
250            if(R(inq).lt.RF(inq)) THEN
251
252              if (iopti.eq.0) then
253
254                   CALL XMIE(R(inq)*1.e6,REALI(K),XIMGI(K),
255     &             QEXT,QSCT,QABS,QBAR,WNOI(K))
256
257              QM1(inq,K)=QEXT
258              QM2(inq,K)=QSCT
259              QM3(inq,K)=QABS
260              QM4(inq,K)=QBAR
261
262              endif
263
264      if (microfi.eq.1) then
265         ig_=ig
266      else
267         ig_=12
268      endif
269
270
271        TAEROS=QM1(inq,K)*qaer(ig_,nlayer+1-J,inq)*1.e-4+TAEROS
272        TAEROSCAT=QM2(inq,K)*qaer(ig_,nlayer+1-J,inq)*1.e-4+TAEROSCAT
273        CBAR=CBAR+QM4(inq,K)*QM2(inq,K)*qaer(ig_,nlayer+1-J,inq)*1.e-4
274
275           else
276
277               XMONO=(r(inq)/RF(inq))**3.
278               XRULE=1.
279
280            if(XMONO.gt.16384./1.5) then
281             XRULE=(XMONO/16384.)
282             XMONO=16384.
283            endif
284
285
286             if (iopti.eq.0) then
287
288       CALL OPTFRAC(XMONO,10000./WNOI(K)
289     &                         ,QEXT,QSCT,QABS,QBAR)
290
291
292c      CALL CFFFV11(1.e-2/WNOI(K),REALI(K),XIMGI(K),RF(inq),2.
293c    &                ,XMONO,QSCT,QEXT,QABS,QBAR)
294
295
296              QF1(inq,K)=QEXT*XRULE
297              QF2(inq,K)=QSCT*XRULE
298              QF3(inq,K)=QABS*XRULE
299              QF4(inq,K)=QBAR
300
301              endif
302
303      if (microfi.eq.1) then
304         ig_=ig
305      else
306         ig_=12
307      endif
308
309        TAEROS=QF1(inq,K)*qaer(ig_,nlayer+1-J,inq)+TAEROS
310        TAEROSCAT=QF2(inq,K)*qaer(ig_,nlayer+1-J,inq)+TAEROSCAT
311        CBAR=CBAR+QF4(inq,K)*QF2(inq,K)*qaer(ig_,nlayer+1-J,inq)
312
313           endif
314
315               IF(TAEROS.LT.1.e-10) TAEROS=1.e-10 
316
317         ENDIF
318       ENDDO             ! FIN DE LA BOUCLE SUR NMICRO
319
320
321
322
323        CBAR=CBAR/TAEROSCAT
324
325        DELTAZ=Z(J)-Z(J+1)
326
327c --------------------------------------------------------------------
328c profil brume Pascal: fit T (sauf tropopause) et albedo
329c -------------------
330        if( cutoff.eq.1) then
331         IF(PRESS(J).gt.9.e-3) THEN
332          TAEROS=TAEROSM1(K)*DELTAZ/DELTAZM1(K)*0.85
333          TAEROSCAT=TAEROSCATM1(K)*DELTAZ/DELTAZM1(K)*0.85
334c         TAEROS=0.
335c         TAEROSCAT=0.
336         ENDIF
337
338         IF(PRESS(J).gt.1.e-1) THEN
339          TAEROS=TAEROSM1(K)*DELTAZ/DELTAZM1(K)*1.15
340          TAEROSCAT=TAEROSCATM1(K)*DELTAZ/DELTAZM1(K)*1.15
341c         TAEROS=0.
342c         TAEROSCAT=0.
343         ENDIF
344        endif !cutoff=1
345
346c profil brume pour fit T (y compris tropopause), mais ne fit plus albedo...
347c -----------------------
348        if( cutoff.eq.2) then
349         IF(PRESS(J).gt.1.e-1) THEN
350          TAEROS=0.
351          TAEROSCAT=0.
352         ENDIF
353        endif !cutoff=2
354c --------------------------------------------------------------------
355
356         TAEROSM1(K)=TAEROS
357         TAEROSCATM1(K)=TAEROSCAT
358         DELTAZM1(K)=DELTAZ
359
360
361      IF(TAEROSCAT.le.0.) CBAR=0.
362
363c     print*,'HERE, MCKAY AEROSOLS IR'
364c     TAEROS=xv1(j,k)
365c     TAEROSCAT=xv2(j,k)
366c     CBAR=xv3(j,k)
367
368c     if (ig.eq.1) then
369c     if (k.eq.NSPECV/2) then
370c      print*,'@IR',K,J,TAEROS,TAEROSCAT,CBAR
371c     stop'Pour faire des comparaisons'
372c     endif
373c     endif
374
375
376c*********** EN TRAVAUX ***************************
377
378C #2:         CLOUD
379C------------------
380
381C NEXT COMPUTE TAU CLOUD
382      TAUCLD=0.0
383       IF ( XNCLD(J) .GT. 0..and .taufac.gt.0.) THEN
384         print*,'Appel a xmie avec radcld=',radcld(j)
385                CALL XMIE(RADCLD(J),RCLDI(K),XICLDI(K),
386     &                         QEXTC,QSCTC,QABSC,CBARC,WNOI(K))
387                TAUCLD=QEXTC*XNCLD(J)
388        ENDIF
389
390
391C #3:          GAZ
392C------------------
393
394C NOW COMPUTE TAUGAS DUE TO THE PIA TERM ONLY FOR LAMDA LT 940
395       TAUGAS=0.0
396       IF (WNOI(K) .LT. 940. ) THEN
397c        if(ig.eq.1.and.k.eq.nspecv/2) print*,'avant PIA'
398                 CALL PIA(K,TBAR,PNN,PCC,PCN,PHN)
399c        if(ig.eq.1.and.k.eq.nspecv/2) print*,'apres PIA'
400C HERE IS WHERE WE COULD SCALE THE PIA COEFFICEINTS TO FIT DATA
401C BASED ON REGIS' NOTES. ---TGM HAS THIS ADJUST IN IT AS DEFAULT
402                 PCN=PCN*MIN(1.75 , AMAX1(1.0,WNOI(K)/200.))
403C***REPLACE ABOVE WITH: PCN=PCN*1.25*MIN(1.75 , AMAX1(1.0,WNOI(K)/200.))
404C 1.25 FACTOR (NOT FROM DATA) SUGGESTED BY TOON et al. (1988)
405                 TAUGAS=COEF1*
406     &           (XN2(J)*XN2(J)*PNN + CH4(J)*CH4(J)*PCC
407     &           + XN2(J)*CH4(J)*PCN + XN2(J)*H2(J)*PHN)
408            IF (J .EQ. NLAYER .AND. IPRINT .GT. 9)
409     &          WRITE (6,22) WNOI(K),TAUGAS,XN2(J),CH4(J),H2(J),
410     &          TBAR, PNN,PCC,PCN, PHN,
411     &          XN2(J)*XN2(J)*PNN , CH4(J)*CH4(J)*PCC ,
412     &          XN2(J)*CH4(J)*PCN , XN2(J)*H2(J)*PHN
413 22             FORMAT(1X,1P8E10.2)
414       ENDIF
415
416       IF (K .GT. 28) THEN
417                KGAS=K-28
418C     ??FLAG? HERE MUST BE WATCHED CAREFULLY
419                     U=COLDEN(J)*6.02204E23/BMU
420          if(ig.eq.1.and.k.eq.nspecv/2) print*,'Avant GAS2'
421                     if((ylellouch).or.(.not.hcnrad)) then
422                       CALL GAS2_NOHCN(J, KGAS,TBAR,PBAR,U,TAU2)
423                     else
424                       CALL GAS2(J, KGAS,TBAR,PBAR,U,TAU2)
425                     endif
426          if(ig.eq.1.and.k.eq.nspecv/2) print*,'Apres GAS2'
427                     TAUGAS=TAUGAS+TAU2
428       ENDIF
429C
430
431      IF (TAEROS + TAUCLD .GT. 0.) THEN
432         COSBI(ig,J,K)=(CBAR*TAEROSCAT + CBARC*TAUCLD )
433     &                     /(TAEROSCAT + TAUCLD)
434      ELSE
435         COSBI(ig,J,K)=0.0
436      ENDIF
437
438      DTAUI(ig,J,K)=TAUGAS+TAEROS+TAUCLD
439
440      TAUHI(ig,K)=TAUHI(ig,K) + TAEROS
441      TAUHID(ig,J,K)=TAUHI(ig,K)
442
443      TAUGI(ig,K)=TAUGI(ig,K) + TAUGAS
444      TAUGID(ig,J,K)=TAUGI(ig,K)
445
446      TAUCI(ig,K)=TAUCI(ig,K) + TAUCLD
447 
448c     if (ig.eq.1) then
449c     if (k.eq.NSPECI/2.or.k.eq.1.or.k.eq.NSPECI) then
450c      print*,'@IR',K,J,DTAUI(ig,J,K),TAUGAS,TAEROS,TAUCLD         
451c     endif
452c     endif
453
454C ??FLAG? SERIOUS PROBLEM WITH THE CODE HERE!
455
456      TLIMIT=1.E-16
457
458      IF (DTAUI(ig,J,K) .GT. TLIMIT) THEN
459
460c***************** ECHANGE
461c        WBARI(J,K)=(QSCT*XNUMB(J) + QSCTC*XNCLD(J))
462c****************
463CFC        WBARI(ig,J,K)=(TAEROSCAT + QSCTC*XNCLD(J))
464c****************
465           WBARI(ig,J,K)=TAEROSCAT
466     &   /DTAUI(ig,J,K)
467        if(iwarning.eq.0)
468     s   print*,'WARNING!!! dans optci xncld non initialise'
469        iwarning=1
470
471      ELSE
472
473         WBARI(ig,J,K)=0.0
474c        PRINT*,'gaz ',taugas,'aerosols',taeros,'nuages',taucld
475c        WRITE (6,999) J,K,DTAUI(ig,J,K)
476 999     FORMAT ('  WARNING TAU MINIMUM AT J,K,DTAU:',2I3,1PE10.3)
477         DTAUI(ig,J,K)=TLIMIT
478
479      ENDIF
480
481c     IF (IPRINT .GT. 9)
482c    & WRITE(6,73)J,K,TAUGAS,TAEROS,QEXT,QSCT
483  73           FORMAT(2I3,1P8E10.3)
484 
485
486c------------------------------------------------------------------------
487 101  CONTINUE   ! FIN BOUCLE L D'O
488c------------------------------------------------------------------------
489 
490      iopti=1
491
492c************************************************************************
493 100  CONTINUE   ! FIN BOUCLE ALTITUDE
494c************************************************************************
495 
496c************************************************************************
497c************************************************************************
498  79  CONTINUE   ! FIN BOUCLE GRILLE HORIZONTALE
499c************************************************************************
500c************************************************************************
501
502
503C TOTAL EXTINCTION OPTICAL DEPTHS
504
505       DO 78 ig=1,klon           ! #1
506
507        DO 119 K=1,NSPECI
508           TAUI(ig,1,K)=0.0
509        DO 119 J=1,NLAYER
510           TAUI(ig,J+1,K)=TAUI(ig,J,K)+DTAUI(ig,J,K)
511 119    CONTINUE
512
513      IF (IPRINT .GT. 2) THEN
514c          WRITE (6,120)
515  120      FORMAT(///'  OPTICAL CONSTANTS IN THE INFRARED')
516
517c        DO 200 K=1,NSPECI           ! #2
518c          WRITE (6,190)
519c          WRITE (6,210)K,WLNI(K),WNOI(K),BWNI(K)
520c    &    ,BWNI(K)+DWNI(K),DWNI(K)
521c          WRITE (6,230)REALI(K),XIMGI(K)
522
523c        DO 195 J=1,NLAYER         !   #3
524c          WRITE (6,220)XNUMB(J), WBARI(ig,J,K),COSBI(ig,J,K)
525c    &                          , DTAUI(ig,J,K),TAUI(ig,J,K)
526c 195    CONTINUE
527
528c        IF(ig.eq.12) WRITE (6,240) TAUI(ig,NLEVEL,K)
529c 200    CONTINUE
530
531          END IF
532  78      CONTINUE
533
534           print*,"TAUI(1400,:,10)=",TAUI(1400,:,10)
535           print*,"DTAUI(1400,:,10)=",DTAUI(1400,:,10)
536
537
538  210 FORMAT(1X,I3,F10.3,F10.2,F10.2,'-',F8.2,F10.3)
539  190 FORMAT(1X//'  SNUM  MICRONS   WAVENU   INTERVAL    DELTA-WN')
540  230 FORMAT(1X,'NREAL(LAYER)= ',1PE10.3,' NIMG(LAYER)= ',E10.3/
541     &' #AEROSOLS   WBAR  COSBAR DTAU TAU')
542  220 FORMAT(5(1X,G9.3))
543  240 FORMAT(41X,G9.3)
544
545c     iopti=1
546      print*, 'FIN OPTCI'
547      stop
548      RETURN
549      END
Note: See TracBrowser for help on using the repository browser.