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

Last change on this file since 3303 was 1461, checked in by emillour, 9 years ago

Titan GCM:
Turned the common block "tgmdat.F" into a module "tgmdat_mod.F90".
This fixes issues in "debug" mode with common variables which seemed to not be correctly shared between routines.
EM
File size: 11.9 KB
Line 
1      SUBROUTINE optcv_1pt2(zqaer_1pt,rcdb,xfrb,ioptv,IPRINT)
2
3      use dimphy
4      USE TGMDAT_MOD, ONLY: RHCH4,FH2,FHAZE,FHVIS,FHIR,TAUFAC,
5     &                      RCLOUD,FARGON
6#include "dimensions.h"
7#include "microtab.h"
8#include "clesphys.h"
9
10      PARAMETER(NLAYER=llm,NLEVEL=NLAYER+1)
11      PARAMETER (NSPECI=46,NSPC1I=47,NSPECV=24,NSPC1V=25)
12
13c   Arguments:
14c   ---------
15      integer IPRINT,ioptv
16C ioptv: premier appel, on ne calcule qu'une fois les QM et QF
17* nrad dans microtab.h
18      real   zqaer_1pt(NLAYER,2*nrad)
19#include "optcv_1pt.h"
20c   ---------
21
22      COMMON /ATM/ Z(NLEVEL),PRESS(NLEVEL),DEN(NLEVEL),TEMP(NLEVEL)
23
24      COMMON /GASS/ CH4(NLEVEL),XN2(NLEVEL),H2(NLEVEL),AR(NLEVEL)
25     & ,XMU(NLEVEL),GAS1(NLAYER),COLDEN(NLAYER)
26
27      COMMON /VISGAS/SOLARF(NSPECV),NTERM(NSPECV),PEXPON(NSPECV),
28     &         ATERM(4,NSPECV),BTERM(4,NSPECV)
29
30      COMMON /AERSOL/ RADIUS(NLAYER), XNUMB(NLAYER)
31     & , REALI(NSPECI), XIMGI(NSPECI), REALV(NSPECV), XIMGV(NSPECV)
32
33      COMMON /CLOUD/
34     &               RCLDI(NSPECI), XICLDI(NSPECI)
35     &             , RCLDV(NSPECV), XICLDV(NSPECV)
36     &             , RCLDI2(NSPECI), XICLDI2(NSPECI)
37     &             , RCLDV2(NSPECV), XICLDV2(NSPECV)
38
39      COMMON /SPECTV/ BWNV(NSPC1V),WNOV(NSPECV)
40     &               ,DWNV(NSPECV),WLNV(NSPECV)
41
42* nrad dans microtab.h
43      COMMON /part/ v(nrad),rayon(nrad),vrat,dr(nrad),dv(nrad)
44
45      REAL QF1(nrad,NSPECV),QF2(nrad,NSPECV)
46      REAL QF3(nrad,NSPECV),QF4(nrad,NSPECV)
47      REAL QM1(nrad,NSPECV),QM2(nrad,NSPECV)
48      REAL QM3(nrad,NSPECV),QM4(nrad,NSPECV)
49      REAL QC1(nrad,NSPECV),QC2(nrad,NSPECV)
50      REAL QC3(nrad,NSPECV),QC4(nrad,NSPECV)
51
52c---- NUAGES
53      real TNUABS,TNUSCAT
54      real   rcdb(NLAYER), xfrb(NLAYER,4)
55
56      save qf1,qf2,qf3,qf4,qm1,qm2,qm3,qm4,qc1,qc2,qc3,qc4
57
58
59C*
60C THIS SUBROUTINE SETS THE OPTICAL CONSTANTS IN THE VISIBLE
61C IT CALCUALTES FOR EACH LAYER, FOR EACH SPECRAL INTERVAL IN THE VIS
62C LAYER: WBAR, DTAU, COSBAR
63C LEVEL: TAU
64C
65C ZERO THE COLUMN OPTICAL DEPTHS OF EACH TYPE
66C ??FLAG? THE OPTICAL DEPTH OF THE TOP OF THE MODEL
67C MAY NOT BE ZERO.
68
69c******* DEBUT DES BOUCLES ************************
70      DO 100 K=1,NSPECV         !b! BOUCLE SUR LAMBDA
71
72      TAURV_1pt(K)=0.
73      TAUHV_1pt(K)=0.            ! INTEGRATED TAU.......INITIALIZATION.
74      TAUCV_1pt(K)=0.            ! Rayleigh, Haze, Cloud, Gas
75      TAUGV_1pt(K)=0.            !   sca,    abs,  abs  , abs
76
77      DO 100 J=1,NLAYER         !a! BOUCLE SUR L"ALTITUDE
78
79C #1:                   HAZE
80c---------------------------
81
82c     CALL THE MIE CODE TO GIVE THE AEROSOL PROPERTIES
83c     USE XFRAC FOR FRACTAL AEROSOLS PROPERTIES AT LAMBDA < 2. um
84
85
86
87
88c                    /\
89c                   /  \
90c                  /    \
91c                 / _O   \
92c                / |/     \
93c               /  / \     \
94c              /   |\ \/\   \
95c             /    || /  \   \
96c             ----------------
97c            |     WARNING    |
98c            |    SLOW DOWN   |
99c             ----------------
100
101
102
103
104c*********** EN TRAVAUX ***************************
105 
106         TAEROS=0.
107         TAEROSCAT=0.
108         CBAR=0.
109
110c       print*,"rayon=",rayon
111c       print*,"RF=",RF
112
113      DO inq=1,nrad         !BOUCLE SUR LES TAILLE D"AEROSOLS
114
115
116            IF (rayon(inq).lt.RF(inq)) THEN    ! aerosols spheriques
117
118             
119            if(ioptv.eq.0.and.J.eq.1) then
120c                  CALL XMIE(rayon(inq)*1.e6,REALV(K),XIMGV(K),
121c    &             QEXT,QSCT,QABS,QBAR,WNOV(K))
122
123                CALL CMIE(1.E-2/WNOV(K),REALV(K),XIMGV(K),rayon(inq),
124     &          QEXT,QSCT,QABS,QBAR)
125
126c       print*,'inq=',inq,' QM1=',QM1(inq,K),' QEXT=',QEXT
127
128              QM1(inq,K)=QEXT
129              QM2(inq,K)=QSCT
130              QM3(inq,K)=QABS
131              QM4(inq,K)=QBAR
132            endif
133
134       TAEROS=QM1(inq,K)*zqaer_1pt(NLAYER+1-J,inq)*1.e-4+TAEROS
135       TAEROSCAT=QM2(inq,K)*zqaer_1pt(NLAYER+1-J,inq)*1.e-4+TAEROSCAT
136       CBAR=CBAR+QM4(inq,K)*QM2(inq,K)*zqaer_1pt(NLAYER+1-J,inq)*1.e-4
137
138            ELSE                        ! aerosols fractals
139
140               XMONO=(rayon(inq)/RF(inq))**3.
141               XRULE=1.
142
143            if(XMONO.gt.16384./1.5) then
144             XRULE=(XMONO/16384.)
145             XMONO=16384.
146            endif
147
148            if(ioptv.eq.0.and.J.eq.1) then
149
150        CALL CFFFV11(1.e-2/WNOV(K),REALV(K),XIMGV(K),RF(inq),2.
151     &   ,XMONO,QSCT,QEXT,QABS,QBAR)
152
153
154              QF1(inq,K)=QEXT*XRULE
155              QF2(inq,K)=QSCT*XRULE
156              QF3(inq,K)=QABS*XRULE
157              QF4(inq,K)=QBAR
158
159c       print*,'inq=',inq,' QF1=',QF1(inq,K),' QEXT=',QEXT,' XRULE=',XRULE
160                   
161            endif
162
163        TAEROS=QF1(inq,K)*zqaer_1pt(NLAYER+1-J,inq)+TAEROS
164        TAEROSCAT=QF2(inq,K)*zqaer_1pt(NLAYER+1-J,inq)+TAEROSCAT
165        CBAR=CBAR+QF4(inq,K)*QF2(inq,K)*zqaer_1pt(NLAYER+1-J,inq)
166
167           ENDIF
168
169       ENDDO    ! nrad
170
171       IF(TAEROSCAT.le.0.) then
172        CBAR=0.
173       ELSE
174        CBAR=CBAR/TAEROSCAT
175       ENDIF
176
177        DELTAZ=Z(J)-Z(J+1)
178
179c --------------------------------------------------------------------
180c profil brume Pascal: fit T (sauf tropopause) et albedo
181c -------------------
182        if( cutoff.eq.1) then
183         IF(PRESS(J).gt.9.e-3) THEN
184          TAEROS=TAEROSM1*DELTAZ/DELTAZM1*0.85
185          TAEROSCAT=TAEROSCATM1*DELTAZ/DELTAZM1*0.85
186c         TAEROS=0.
187c         TAEROSCAT=0.
188         ENDIF
189
190         IF(PRESS(J).gt.1.e-1) THEN
191          TAEROS=TAEROSM1*DELTAZ/DELTAZM1*1.15
192          TAEROSCAT=TAEROSCATM1*DELTAZ/DELTAZM1*1.15
193c         TAEROS=0.
194c         TAEROSCAT=0.
195         ENDIF
196        endif !cutoff=1
197
198c profil brume pour fit T (y compris tropopause), mais ne fit plus albedo...
199c -----------------------
200        if( cutoff.eq.2) then
201         IF(PRESS(J).gt.1.e-1) THEN
202          TAEROS=0.
203          TAEROSCAT=0.
204         ENDIF
205        endif !cutoff=2
206c --------------------------------------------------------------------
207
208         TAEROSM1=TAEROS
209         TAEROSCATM1=TAEROSCAT
210         DELTAZM1=DELTAZ
211
212
213       IF (TAEROSCAT.le.0.) CBAR=0.
214
2151699  FORMAT(a3,2I3,3(ES15.7,1X))
216
217c*********** EN TRAVAUX ***************************
218
219C #2:                   RAYLEIGH
220c-------------------------------
221
222C RAYLEIGH SCATTERING STRAIGHT FROM HANSEN AND TRAVIS...SEE NOTES
223C RATIOED BY THE LAYER COLUMN NUMBER TO THE TOTAL
224C COLUMN NUMBER ON EARTH. CM-2
225C THIS IS THE SCATTERING BY THE ATMOSPHERE
226
227      TAURAY=(COLDEN(J)*28.9/(XMU(J)*1013.25))*
228     &(.008569/WLNV(K)**4)*(1.+.0113/WLNV(K)**2+.00013/WLNV(K)**4)
229
230
231C #3:                   CLOUD
232c----------------------------
233C NEXT COMPUTE TAU CLOUD
234c
235c  Menu special :
236c  Afin d'eviter la surcharge de calcul on ne calcule les
237c  propriétes optiques des nuages qu'une seule fois
238c  avec un rayon de particule effectif de 3um et une composition
239c  de goutte : 90% CH4 / 10% NOYAUX
240c  Puis on ajute les section efficace par la surface reelle de
241c  la goutte.
242c
243c  ---> A TESTER !!!!
244c
245       IF (clouds.eq.0) THEN
246         CNBAR=0.
247         TNUSCAT=0.
248         TNUABS=0.
249       ELSE
250         IF (ioptv.eq.0.and.j.eq.1) THEN !--> au premier appel
251           QEXTC=0.
252           QSCTC=0.
253           QABSC=0.
254           CBARC=0.
255           DO inq=1,nrad         !BOUCLE SUR LES NQMX TAILLE D"AEROSOLS
256             QC1(inq,K)=0.
257             QC2(inq,K)=0.
258             QC3(inq,K)=0.
259             QC4(inq,K)=0.
260           ENDDO
261** OPTICAL CONSTANT : MIXING RULES
262** Fraction volumique fixe :
263** 10% noyaux.
264** 90% methane.
265           XNR = 0.5 * REALI(K)
266     &         + 0.5 * RCLDI(K)
267           XNI = 0.5 * XIMGI(K)
268     &         + 0.5 * XICLDI(K)
269**
270**   Efficacite : particule de 3um de rayon
271           CALL CMIE(1.E-2/WNOV(K),XNR,XNI,3.e-6,
272     &                QEXTC,QSCTC,QABSC,CBARC)
273
274           DO inq=1,nrad
275             QC1(inq,K)=QEXTC/xnuf
276             QC2(inq,K)=QSCTC/xnuf
277             QC3(inq,K)=QABSC/xnuf
278             QC4(inq,K)=CBARC
279           ENDDO
280         ENDIF   ! ioptv = 0
281         TNUABS=0.
282         TNUSCAT=0.
283         CNBAR=0.
284         IF (rcdb(nlayer+1-J).gt.1.1e-10) THEN
285           DO inq=1,nrad
286             TNUABS=QC1(inq,K)*(rcdb(nlayer+1-J)/3.e-6)**2.*1.e-4*
287     &              zqaer_1pt(NLAYER+1-J,inq+nrad) +
288     &              TNUABS
289             TNUSCAT=QC2(inq,K)*(rcdb(nlayer+1-J)/3.e-6)**2.*1.e-4*
290     &               zqaer_1pt(NLAYER+1-J,inq+nrad) +
291     &               TNUSCAT
292             CNBAR=QC4(inq,K)*QC2(inq,K)*(rcdb(nlayer+1-J)/3.e-6)**2.*
293     &             1.e-4*zqaer_1pt(NLAYER+1-J,inq+nrad) +
294     &             CNBAR
295           ENDDO
296         ENDIF
297
298         IF(TNUSCAT.EQ.0.) THEN
299           CNBAR=0.
300         ELSE
301           CNBAR=CNBAR/TNUSCAT
302         ENDIF
303       ENDIF  ! Cond. CLD
304
305       TAUCV_1pt(K)=TAUCV_1pt(K)+TNUABS
306       TAUCVD_1pt(J,K)=TAUCV_1pt(K)
307
308       TAURV_1pt(K)=TAURV_1pt(K)+TAURAY
309       TAUGVD_1pt(J,K)=TAURV_1pt(K)
310
311       TAUHV_1pt(K)=TAUHV_1pt(K)+TAEROS ! INTEGRATED Quant.
312       TAUHVD_1pt(J,K)=TAUHV_1pt(K)
313
314
315
316C #4:                  TAUGAS
317C----------------------------
318
319C LOOP OVER THE NTERMS
320C THIS IS THE ABSORPTION BY THE ATMOSPHERE (METHANE)
321
322
323       DO 909 NT=1,NTERM(K)
324         TAUGAS=COLDEN(J)*GAS1(J)*BTERM(NT,K)*
325     &   (   (PRESS(J+1) + PRESS(J))*.5  )**PEXPON(K)
326
327
328*  COSBV ET COSBVP
329*-----------------
330
331         IF(TAEROSCAT+TNUSCAT+TAURAY .ne. 0.) THEN
332           COSBV_1pt(J,K,NT)=(CBAR*TAEROSCAT + CNBAR*TNUSCAT)
333     &     /(TAEROSCAT+TNUSCAT+TAURAY) !CBAR_RAY=0.
334         ELSE
335           COSBV_1pt(J,K,NT)=0.
336         ENDIF
337
338         IF(TAEROSCAT+TAURAY .ne. 0.) THEN
339           COSBVP_1pt(J,K,NT)=(CBAR*TAEROSCAT)
340     &     /(TAEROSCAT+TAURAY) !CBAR_RAY=0.
341         ELSE
342           COSBVP_1pt(J,K,NT)=0.
343         ENDIF
344
345*  DTAUV ET DTAUVP
346*-----------------
347
348         DTAUV_1pt(J,K,NT) =TAUGAS+TAEROS+TAURAY+TNUABS !TAU_ABS_METH
349         DTAUVP_1pt(J,K,NT)=TAUGAS+TAEROS+TAURAY       !TAU_ABS_METH
350
351         TAUGV_1pt(K)=TAUGV_1pt(K)+TAUGAS*ATERM(NT,K) !INTEG.
352
353*  WBARV ET WBARVP
354*-----------------
355
356         IF(TAUGAS+TAEROS+TAURAY+TNUABS .ne.  0.) THEN
357           WBARV_1pt(J,K,NT)=(TAEROSCAT+TAURAY*0.9999999 + TNUSCAT)
358     &     /(TAUGAS+TAEROS+TAURAY+TNUABS)
359         ELSE
360           WBARV_1pt(J,K,NT)=0.
361         ENDIF
362
363         IF(TAUGAS+TAEROS+TAURAY .ne.  0.) THEN
364           WBARVP_1pt(J,K,NT)=(TAEROSCAT+TAURAY*0.9999999 )
365     &     /(TAUGAS+TAEROS+TAURAY)
366         ELSE
367           WBARVP_1pt(J,K,NT)=0.
368         ENDIF
369   
370 909   CONTINUE
371 
372       TAUGVD_1pt(J,K)=TAUGVD_1pt(J,K)+TAUGV_1pt(K)
373
374 100  CONTINUE
375
376       ioptv=1
377
378c HERE END OF THE LOOPS *******
379c******************************
380         
381C TOTAL EXTINCTION OPTICAL DEPTHS
382          DO 119 K=1,NSPECV
383C LOOP OVER NTERMS
384           DO 119 NT=1,NTERM(K)
385           TAUV_1pt(1,K,NT)=0.0
386           TAUVP_1pt(1,K,NT)=0.0
387             DO 119 J=1,NLAYER
388             TAUV_1pt(J+1,K,NT)=TAUV_1pt(J,K,NT)+DTAUV_1pt(J,K,NT)
389             TAUVP_1pt(J+1,K,NT)=TAUVP_1pt(J,K,NT)+DTAUVP_1pt(J,K,NT)
390 119     CONTINUE
391
392
393c       print*,'SETUP'
394c      do i=1,NSPECV
395c      print*,WLNV(i)
396c       do j=1,NLAYER+1
397c       print*,Z(j),TAUV(1,j,i,1),WBARV(1,j,i,1),COSBV(1,j,i,1)
398c       enddo
399c      enddo
400c
401c     IF (IPRINT .GT. 1) THEN
402c           NT=1
403c     IF (2 .GT. 1) THEN
404c          WRITE (6,120)
405c 120      FORMAT(///'  OPTICAL CONSTANTS IN THE VISIBLE (@EQUATOR) ')
406c          WRITE(6,*) 'latitude:',ig
407c          DO 200 K=1,NSPECV
408c          WRITE (6,190)
409c          WRITE (6,210)K,WLNV(K),WNOV(K),BWNV(K)
410c    &    ,BWNV(K)+DWNV(K),DWNV(K)
411c          WRITE (6,230)REALV(K),XIMGV(K)
412c          DO 195 J=1,NLAYER,NLAYER
413C RECALCULATE FOR PRINT OUT ONLY, ONLY FIRST NTERM AT ig=12 (EQUATOR)
414c          WRITE (6,220)XNUMB(J), WBARV_1pt(J,K,NT),COSBV_1pt(J,K,NT)
415c    &      ,DTAUV_1pt(J,K,NT),TAUV_1pt(J,K,NT)
416c 195      CONTINUE
417c          WRITE (6,240) TAUV_1pt(NLEVEL,K,NT)
418c 200      CONTINUE
419c     END IF
420
421c  210 FORMAT(1X,I3,F10.3,F10.2,F10.2,'-',F8.2,F10.3)
422c  190 FORMAT(1X//'  SNUM  MICRONS   WAVENU   INTERVAL    DELTA-WN')
423c  230 FORMAT(1X,'NREAL(LAYER)= ',1PE10.3,' NIMG(LAYER)= ',E10.3/
424c     &' #AEROSOLS   WBAR  COSBAR       DTAU     TAU'
425c     & ,9X,'RAY     GAS    AEROSOL')
426c  220 FORMAT(8(1X,F9.3))
427c  240 FORMAT(41X,F9.3)
428
429      RETURN
430      END
Note: See TracBrowser for help on using the repository browser.