Context Navigation

optci_1pt_2.F @ 201

Last change on this file since 201 was 175, checked in by slebonnois, 14 years ago

S.LEBONNOIS:

Revision majeure de la physique Titan => ajout des nuages version 10 bins (Jeremie Burgalat) Cette version reste a tester mais avec clouds=0, on reste sur l'ancienne.
Quelques ajouts dans la doc.

File size: 13.7 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format