source: trunk/LMDZ.TITAN/libf/phytitan/optci_1pt_3.F @ 1056

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

SL: Titan runs ! see DOC/chantiers/commit_importants.log

File size: 13.2 KB
Line 
1      SUBROUTINE optci_1pt3(zqaer_1pt,rcdb,xfrb,iopti,IPRINT)
2      use dimphy
3      IMPLICIT NONE
4#include "dimensions.h"
5#include "microtab.h"
6#include "numchimrad.h"
7#include "clesphys.h"
8
9      integer nlayer, nlevel, nspeci, nspc1i, nspecv, nspc1v
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,
14     &     f0pi, rhch4, fh2, fhaze, fhvis
15     &     reali, ximgi, bwni, fhir, taufac, rcloud, fargon, rgas, rhop,
16     &     pi, sigma, prod,reali,fhvis
17
18      integer k, j,inq,kgas
19
20      real tbar, pbar, bmu, coef1, effg, taeros, taeroscat, cbar,
21     &     qext, qsct, qabs, qbar, xmono, xrule, deltaz, tnuext,
22     &     tnuscat, cnbar, qextc, qsctc, qabsc, qbarc, taugas, pnn,
23     &     pcc, pcn, phn, u, ig, tau2, tlimit
24
25      PARAMETER(NLAYER=llm,NLEVEL=NLAYER+1)
26      PARAMETER (NSPECI=46,NSPC1I=47,NSPECV=24,NSPC1V=25)
27
28c   Arguments:
29c   ---------
30      integer IPRINT,iopti
31C iopti: premier appel, on ne calcule qu'une fois les QM et QF
32* nrad dans microtab.h
33      real   zqaer_1pt(NLAYER,2*nrad)
34#include "optci_1pt.h"
35c   ---------
36
37      COMMON /ATM/ Z(NLEVEL),PRESS(NLEVEL),DEN(NLEVEL),TEMP(NLEVEL)
38
39      COMMON /GASS/ CH4(NLEVEL),XN2(NLEVEL),H2(NLEVEL),AR(NLEVEL)
40     & ,XMU(NLEVEL),GAS1(NLAYER),COLDEN(NLAYER)
41
42      COMMON /STRATO/ C2H2(NLAYER),C2H6(NLAYER)
43      COMMON /STRAT2/ HCN(NLAYER)
44
45      COMMON /AERSOL/ RADIUS(NLAYER), XNUMB(NLAYER)
46     & , REALI(NSPECI), XIMGI(NSPECI), REALV(NSPECV), XIMGV(NSPECV)
47
48      COMMON /CLOUD/
49     &               RCLDI(NSPECI), XICLDI(NSPECI)
50     &             , RCLDV(NSPECV), XICLDV(NSPECV)
51     &             , RCLDI2(NSPECI), XICLDI2(NSPECI)
52     &             , RCLDV2(NSPECV), XICLDV2(NSPECV)
53
54      COMMON /SPECTI/ BWNI(NSPC1I), WNOI(NSPECI),
55     &                DWNI(NSPECI), WLNI(NSPECI)
56
57      COMMON /PLANT/ CSUBP,F0PI
58      COMMON /ADJUST/ RHCH4,FH2,FHAZE,FHVIS,FHIR,TAUFAC,RCLOUD,FARGON
59      COMMON /CONST/RGAS,RHOP,PI,SIGMA
60      COMMON /part/v,rayon,vrat,dr,dv
61
62      DIMENSION PROD(NLEVEL)
63* nrad dans microtab.h
64      real v(nrad),rayon(nrad),vrat,dr(nrad),dv(nrad)
65      real xv1(klev,nspeci),xv2(klev,nspeci)
66      real xv3(klev,nspeci)
67      REAL QF1(nrad,NSPECI),QF2(nrad,NSPECI)
68      REAL QF3(nrad,NSPECI),QF4(nrad,NSPECI)
69      REAL QM1(nrad,NSPECI),QM2(nrad,NSPECI)
70      REAL QM3(nrad,NSPECI),QM4(nrad,NSPECI)
71      REAL QC1(nrad,NSPECI),QC2(nrad,NSPECI)
72      REAL QC3(nrad,NSPECI),QC4(nrad,NSPECI)
73      real emu
74      REAL TAEROSM1(NSPECI),TAEROSCATM1(NSPECI),DELTAZM1(NSPECI)
75
76c ---- nuages     
77      REAL TNUAGE,TNUAGESCAT
78      REAL rcdb(nlayer),xfrb(nlayer,4)
79
80      save qf1,qf2,qf3,qf4,qm1,qm2,qm3,qm4,qc1,qc2,qc3,qc4
81
82
83C THE PRESSURE INDUCED TRANSITIONS ARE FROM REGIS
84C THE LAST SEVENTEEN INTERVALS ARE THE BANDS FROM GNF.
85C
86C THIS SUBROUTINE SETS THE OPTICAL CONSTANTS IN THE INFRARED
87C IT CALCUALTES FOR EACH LAYER, FOR EACH SPECRAL INTERVAL IN THE IR
88C LAYER: WBAR, DTAU, COSBAR
89C LEVEL: TAU
90C
91
92        DO 80 K=1,NSPECI
93           TAUHI_1pt(K)=0.
94           TAUCI_1pt(K)=0.
95           TAUGI_1pt(K)=0.
96 80     CONTINUE
97
98c************************************************************************
99          DO 100 J=1,NLAYER    ! BOUCLE SUR L'ALTITUDE
100c************************************************************************
101
102C SET UP THE COEFFICIENT TO REDUCE MASS PATH TO STP ...SEE NOTES
103C T0 =273.15   PO=1.01325 BAR
104
105        TBAR=0.5*(TEMP(J)+TEMP(J+1))
106        PBAR=SQRT(PRESS(J)*PRESS(J+1))
107        BMU=0.5*(XMU(J+1)+XMU(J))
108c attention ici, Z en km doit etre passe en m
109        COEF1=RGAS*273.15**2*.5E5* (PRESS(J+1)**2 - PRESS(J)**2)
110     & /(1.01325**2 *EFFG(Z(J)*1000.)*TBAR*BMU)
111
112      IF (IPRINT .GT. 9) WRITE(6,21) J,EFFG(Z(J)*1000.),TBAR,BMU,COEF1
113 21   FORMAT(' J, EFFG, TBAR, BMU, COEF1,: ',I3,1P6E10.3)
114
115c------------------------------------------------------------------------
116         DO 101 K=1,NSPECI     ! BOUCLE SUR LES L.D'O
117c------------------------------------------------------------------------
118
119
120C #1:             HAZE
121C---------------------
122
123C FIRST COMPUTE TAU AEROSOL
124
125
126c
127c                    /\
128c                   /  \
129c                  /    \
130c                 / _O   \
131c                / |/     \
132c               /  / \     \
133c              /   |\ \/\   \
134c             /    || /  \   \
135c             ----------------
136c            |     WARNING    |
137c            |    SLOW DOWN   |
138c             ----------------
139
140
141
142
143c*********** EN TRAVAUX ***************************
144
145         TAEROS=0.
146         TAEROSCAT=0.
147         CBAR=0.
148
149
150      DO inq=1,nrad         !BOUCLE SUR LES TAILLE D"AEROSOLS
151
152
153      IF (WNOI(K).lt.wco) THEN    ! lamda > 56 um
154
155           if (iopti.eq.0) then
156
157c          CALL XMIE(rayon(inq)*1.e6,REALI(K),XIMGI(K),
158c    &     QEXT,QSCT,QABS,QBAR,WNOI(K))
159
160           CALL CMIE(1.E-2/WNOI(K),REALI(K),XIMGI(K),rayon(inq),
161     &     QEXT,QSCT,QABS,QBAR)
162
163
164           QM1(inq,K)=QEXT
165           QM2(inq,K)=QSCT
166           QM3(inq,K)=QABS
167           QM4(inq,K)=QBAR
168
169          endif         ! end iopti
170
171
172      TAEROS=QM1(inq,K)*zqaer_1pt(nlayer+1-J,inq)*1.e-4+TAEROS
173      TAEROSCAT=QM2(inq,K)*zqaer_1pt(nlayer+1-J,inq)*1.e-4+TAEROSCAT
174      CBAR=CBAR+QM4(inq,K)*QM2(inq,K)*zqaer_1pt(nlayer+1-J,inq)*1.e-4
175
176
177         ELSE                           ! 0.2 < lambda < 56 um
178
179
180            if(rayon(inq).lt.RF(inq)) THEN
181
182              if (iopti.eq.0) then
183
184                   CALL XMIE(rayon(inq)*1.e6,REALI(K),XIMGI(K),
185     &             QEXT,QSCT,QABS,QBAR,WNOI(K))
186
187              QM1(inq,K)=QEXT
188              QM2(inq,K)=QSCT
189              QM3(inq,K)=QABS
190              QM4(inq,K)=QBAR
191              endif         ! end iopti
192
193
194        TAEROS=QM1(inq,K)*zqaer_1pt(nlayer+1-J,inq)*1.e-4+TAEROS
195        TAEROSCAT=QM2(inq,K)*zqaer_1pt(nlayer+1-J,inq)*1.e-4+TAEROSCAT
196        CBAR=CBAR+QM4(inq,K)*QM2(inq,K)*zqaer_1pt(nlayer+1-J,inq)*1.e-4
197
198           else
199
200               XMONO=(rayon(inq)/RF(inq))**3.
201               XRULE=1.
202
203            if(XMONO.gt.16384./1.5) then
204             XRULE=(XMONO/16384.)
205             XMONO=16384.
206            endif
207
208
209             if (iopti.eq.0) then
210
211c       CALL OPTFRAC(XMONO,10000./WNOI(K)
212c     &                         ,QEXT,QSCT,QABS,QBAR)
213
214
215       CALL CFFFV11(1.e-2/WNOI(K),REALI(K),XIMGI(K),RF(inq),2.
216     &                ,XMONO,QSCT,QEXT,QABS,QBAR)
217
218
219              QF1(inq,K)=QEXT*XRULE
220              QF2(inq,K)=QSCT*XRULE
221              QF3(inq,K)=QABS*XRULE
222              QF4(inq,K)=QBAR
223             endif         ! end iopti
224
225        TAEROS=QF1(inq,K)*zqaer_1pt(nlayer+1-J,inq)+TAEROS
226        TAEROSCAT=QF2(inq,K)*zqaer_1pt(nlayer+1-J,inq)+TAEROSCAT
227        CBAR=CBAR+QF4(inq,K)*QF2(inq,K)*zqaer_1pt(nlayer+1-J,inq)
228
229           endif
230
231               IF(TAEROS.LT.1.e-10) TAEROS=1.e-10
232
233         ENDIF
234       ENDDO             ! FIN DE LA BOUCLE SUR nrad
235
236
237
238
239        CBAR=CBAR/TAEROSCAT
240
241        DELTAZ=Z(J)-Z(J+1)
242
243c --------------------------------------------------------------------
244c profil brume Pascal: fit T (sauf tropopause) et albedo
245c -------------------
246        if( cutoff.eq.1) then
247         IF(PRESS(J).gt.9.e-3) THEN
248          TAEROS=TAEROSM1(K)*DELTAZ/DELTAZM1(K)*0.85
249          TAEROSCAT=TAEROSCATM1(K)*DELTAZ/DELTAZM1(K)*0.85
250c         TAEROS=0.
251c         TAEROSCAT=0.
252         ENDIF
253
254         IF(PRESS(J).gt.1.e-1) THEN
255          TAEROS=TAEROSM1(K)*DELTAZ/DELTAZM1(K)*1.15
256          TAEROSCAT=TAEROSCATM1(K)*DELTAZ/DELTAZM1(K)*1.15
257c         TAEROS=0.
258c         TAEROSCAT=0.
259         ENDIF
260        endif !cutoff=1
261
262c profil brume pour fit T (y compris tropopause), mais ne fit plus albedo...
263c -----------------------
264        if( cutoff.eq.2) then
265         IF(PRESS(J).gt.1.e-1) THEN
266          TAEROS=0.
267          TAEROSCAT=0.
268         ENDIF
269        endif !cutoff=2
270c --------------------------------------------------------------------
271
272         TAEROSM1(K)=TAEROS
273         TAEROSCATM1(K)=TAEROSCAT
274         DELTAZM1(K)=DELTAZ
275
276
277      IF(TAEROSCAT.le.0.) CBAR=0.
278
279c     print*,'HERE, MCKAY AEROSOLS IR'
280c     TAEROS=xv1(j,k)
281c     TAEROSCAT=xv2(j,k)
282c     CBAR=xv3(j,k)
283
284
285c*********** EN TRAVAUX ***************************
286
287C #2:         CLOUD
288C------------------
289C NEXT COMPUTE TAU CLOUD
290c
291c  Menu special :
292c  On utilise ici une look-up table afin de calculer
293c  les proprietes optique des nuages.
294c  Le principe est le suivant :
295c  La look-up table contient les proprietes optique d'une goutte
296c  de methane pur de 3 um.
297c  On approxime les proprietes optiques pour une goutte de rayon r a
298c  de la table.
299c
300c
301        TNUEXT=0.
302        TNUSCAT=0.
303        CNBAR=0.
304        IF (clouds.eq.1) THEN
305
306          CALL getoptcld(1.E-2/WNOI(K),rcdb(nlayer+1-J),
307     &                   QEXTC,QSCTC,QABSC,QBARC)
308
309
310c ----- On ne calcule les constante optiques que si Rgoutte > 1e-10
311          IF (rcdb(nlayer+1-J).gt.1.1e-10) THEN
312            TNUEXT =QEXTC/xnuf*SUM(zqaer_1pt(NLAYER+1-J,nrad+1:2*nrad))
313            TNUSCAT=QSCTC/xnuf*SUM(zqaer_1pt(NLAYER+1-J,nrad+1:2*nrad))
314            CNBAR  =QBARC
315          ENDIF
316
317          IF(TNUSCAT.EQ.0.) THEN
318            CNBAR=0.
319          ELSE
320            CNBAR=CNBAR/TNUSCAT
321          ENDIF
322
323        ENDIF    ! Cond CLD
324c
325C #3:          GAZ
326C------------------
327
328C NOW COMPUTE TAUGAS DUE TO THE PIA TERM ONLY FOR LAMDA LT 940
329       TAUGAS=0.0
330       IF (WNOI(K) .LT. 940. ) THEN
331                 CALL PIA(K,TBAR,PNN,PCC,PCN,PHN)
332C HERE IS WHERE WE COULD SCALE THE PIA COEFFICEINTS TO FIT DATA
333C BASED ON REGIS' NOTES. ---TGM HAS THIS ADJUST IN IT AS DEFAULT
334                 PCN=PCN*MIN(1.75 , AMAX1(1.0,WNOI(K)/200.))
335C***REPLACE ABOVE WITH: PCN=PCN*1.25*MIN(1.75 , AMAX1(1.0,WNOI(K)/200.))
336C 1.25 FACTOR (NOT FROM DATA) SUGGESTED BY TOON et al. (1988)
337                 TAUGAS=COEF1*
338     &           (XN2(J)*XN2(J)*PNN + CH4(J)*CH4(J)*PCC
339     &           + XN2(J)*CH4(J)*PCN + XN2(J)*H2(J)*PHN)
340            IF (J .EQ. NLAYER .AND. IPRINT .GT. 9)
341     &          WRITE (6,22) WNOI(K),TAUGAS,XN2(J),CH4(J),H2(J),
342     &          TBAR, PNN,PCC,PCN, PHN,
343     &          XN2(J)*XN2(J)*PNN , CH4(J)*CH4(J)*PCC ,
344     &          XN2(J)*CH4(J)*PCN , XN2(J)*H2(J)*PHN
345 22             FORMAT(1X,1P8E10.2)
346       ENDIF
347
348       IF (K .GT. 28) THEN
349                KGAS=K-28
350C     ??FLAG? HERE MUST BE WATCHED CAREFULLY
351                     U=COLDEN(J)*6.02204E23/BMU
352                     if((ylellouch).or.(.not.hcnrad)) then
353                       CALL GAS2_NOHCN(J, KGAS,TBAR,PBAR,U,TAU2)
354                     else
355                       CALL GAS2(J, KGAS,TBAR,PBAR,U,TAU2)
356                     endif
357                     TAUGAS=TAUGAS+TAU2
358       ENDIF
359C
360
361      DTAUI_1pt(J,K)=TAUGAS+TAEROS+TNUEXT
362      DTAUIP_1pt(J,K)=TAUGAS+TAEROS
363
364      TAUHI_1pt(K)=TAUHI_1pt(K) + TAEROS
365      TAUHID_1pt(J,K)=TAUHI_1pt(K)
366
367      TAUGI_1pt(K)=TAUGI_1pt(K) + TAUGAS
368      TAUGID_1pt(J,K)=TAUGI_1pt(K)
369
370      TAUCI_1pt(K)=TAUCI_1pt(K) + TNUEXT
371      TAUCID_1pt(J,K)=TAUCI_1pt(K)
372 
373C ??FLAG? SERIOUS PROBLEM WITH THE CODE HERE!
374
375      TLIMIT=1.E-16
376
377
378      IF (TAEROSCAT + TNUSCAT .GT. 0.) THEN
379         COSBI_1pt(J,K)=(CBAR*TAEROSCAT + CNBAR*TNUSCAT )
380     &                     /(TAEROSCAT + TNUSCAT)
381      ELSE
382         COSBI_1pt(J,K)=0.0
383      ENDIF
384
385      IF (TAEROSCAT  .GT. 0.) THEN
386         COSBIP_1pt(J,K)=(CBAR*TAEROSCAT)
387     &                     /(TAEROSCAT)
388      ELSE
389         COSBIP_1pt(J,K)=0.0
390      ENDIF
391
392*---------
393
394      IF (DTAUI_1pt(J,K) .GT.  TLIMIT) THEN
395          WBARI_1pt(J,K)=(TAEROSCAT+TNUSCAT) /DTAUI_1pt(J,K)
396      ELSE
397         WBARI_1pt(J,K)=0.0
398         DTAUI_1pt(J,K)=TLIMIT
399      ENDIF
400
401      IF (DTAUIP_1pt(J,K) .GT.  TLIMIT) THEN
402          WBARIP_1pt(J,K)=(TAEROSCAT) /DTAUIP_1pt(J,K)
403      ELSE
404         WBARIP_1pt(J,K)=0.0
405         DTAUIP_1pt(J,K)=TLIMIT
406      ENDIF
407
408
409c     IF (IPRINT .GT. 9)
410c    & WRITE(6,73)J,K,TAUGAS,TAEROS,QEXT,QSCT
411  73           FORMAT(2I3,1P8E10.3)
412 
413
414c------------------------------------------------------------------------
415 101  CONTINUE   ! FIN BOUCLE L D'O
416c------------------------------------------------------------------------
417 
418      iopti=1
419
420c************************************************************************
421 100  CONTINUE   ! FIN BOUCLE ALTITUDE
422c************************************************************************
423 
424        DO 119 K=1,NSPECI
425           TAUI_1pt(1,K)=0.0
426           TAUIP_1pt(1,K)=0.0
427        DO 119 J=1,NLAYER
428           TAUI_1pt(J+1,K)=TAUI_1pt(J,K)+DTAUI_1pt(J,K)
429           TAUIP_1pt(J+1,K)=TAUIP_1pt(J,K)+DTAUIP_1pt(J,K)
430 119    CONTINUE
431
432c      IF (IPRINT .GT. 2) THEN
433c          WRITE (6,120)
434c  120      FORMAT(///'  OPTICAL CONSTANTS IN THE INFRARED')
435
436c        DO 200 K=1,NSPECI           ! #2
437c          WRITE (6,190)
438c          WRITE (6,210)K,WLNI(K),WNOI(K),BWNI(K)
439c    &    ,BWNI(K)+DWNI(K),DWNI(K)
440c          WRITE (6,230)REALI(K),XIMGI(K)
441
442c        DO 195 J=1,NLAYER         !   #3
443c          WRITE (6,220)XNUMB(J), WBARI_1pt(J,K),COSBI_1pt(J,K)
444c    &                          , DTAUI_1pt(J,K),TAUI_1pt(J,K)
445c 195    CONTINUE
446
447c 200    CONTINUE
448
449c          END IF
450
451
452c  210 FORMAT(1X,I3,F10.3,F10.2,F10.2,'-',F8.2,F10.3)
453c  190 FORMAT(1X//'  SNUM  MICRONS   WAVENU   INTERVAL    DELTA-WN')
454c  230 FORMAT(1X,'NREAL(LAYER)= ',1PE10.3,' NIMG(LAYER)= ',E10.3/
455c     &' #AEROSOLS   WBAR  COSBAR DTAU TAU')
456c  220 FORMAT(5(1X,G9.3))
457c  240 FORMAT(41X,G9.3)
458
459      RETURN
460      END
Note: See TracBrowser for help on using the repository browser.