Context Navigation

optci_1pt.F @ 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: 12.4 KB

Line
1	SUBROUTINE optci_1pt(zqaer_1pt,iopti,
2	. COSBI_1pt,DTAUI_1pt,TAUHI_1pt,TAUHID_1pt,TAUCI_1pt,
3	. TAUGI_1pt,TAUGID_1pt,WBARI_1pt,TAUI_1pt,IPRINT)
4	#include "dimensions.h"
5	#include "dimphy.h"
6	#include "microtab.h"
7	#include "numchimrad.h"
8	#include "clesphys.h"
9
10	PARAMETER(NLAYER=llm,NLEVEL=NLAYER+1)
11	PARAMETER (NSPECI=46,NSPC1I=47,NSPECV=24,NSPC1V=25)
12
13	c Arguments:
14	c ---------
15	integer IPRINT,iopti
16	C iopti: premier appel, on ne calcule qu'une fois les QM et QF
17	* nrad dans microtab.h
18	real zqaer_1pt(NLAYER,nrad)
19	real TAUHID_1pt(NLAYER,NSPECI)
20	real TAUGID_1pt(NLAYER,NSPECI)
21	real TAUHI_1pt(NSPECI),TAUCI_1pt(NSPECI)
22	real TAUGI_1pt(NSPECI)
23	real DTAUI_1pt(NLAYER,NSPECI),TAUI_1pt(NLEVEL,NSPECI)
24	real WBARI_1pt(NLAYER,NSPECI)
25	real COSBI_1pt(NLAYER,NSPECI)
26	c ---------
27
28	COMMON /ATM/ Z(NLEVEL),PRESS(NLEVEL),DEN(NLEVEL),TEMP(NLEVEL)
29
30	COMMON /GASS/ CH4(NLEVEL),XN2(NLEVEL),H2(NLEVEL),AR(NLEVEL)
31	& ,XMU(NLEVEL),GAS1(NLAYER),COLDEN(NLAYER)
32
33	COMMON /STRATO/ C2H2(NLAYER),C2H6(NLAYER)
34	COMMON /STRAT2/ HCN(NLAYER)
35
36	COMMON /AERSOL/ RADIUS(NLAYER), XNUMB(NLAYER)
37	& , REALI(NSPECI), XIMGI(NSPECI), REALV(NSPECV), XIMGV(NSPECV)
38
39	COMMON /CLOUD/ RADCLD(NLAYER), XNCLD(NLAYER)
40	& , RCLDI(NSPECI), XICLDI(NSPECI), RCLDV(NSPECV), XICLDV(NSPECV)
41
42	COMMON /SPECTI/ BWNI(NSPC1I), WNOI(NSPECI),
43	& DWNI(NSPECI), WLNI(NSPECI)
44
45	COMMON /PLANT/ CSUBP,RSFI,RSFV,F0PI
46	COMMON /ADJUST/ RHCH4,FH2,FHAZE,FHVIS,FHIR,TAUFAC,RCLOUD,FARGON
47	COMMON /CONST/RGAS,RHOP,PI,SIGMA
48	COMMON /part/v,rayon,vrat,dr,dv
49
50	DIMENSION PROD(NLEVEL)
51	* nrad dans microtab.h
52	real v(nrad),rayon(nrad),vrat,dr(nrad),dv(nrad)
53	real xv1(klev,nspeci),xv2(klev,nspeci)
54	real xv3(klev,nspeci)
55	REAL QF1(nrad,NSPECI),QF2(nrad,NSPECI)
56	REAL QF3(nrad,NSPECI),QF4(nrad,NSPECI)
57	REAL QM1(nrad,NSPECI),QM2(nrad,NSPECI)
58	REAL QM3(nrad,NSPECI),QM4(nrad,NSPECI)
59	real emu
60	REAL TAEROSM1(NSPECI),TAEROSCATM1(NSPECI),DELTAZM1(NSPECI)
61
62	save qf1,qf2,qf3,qf4,qm1,qm2,qm3,qm4
63
64	C THE PRESSURE INDUCED TRANSITIONS ARE FROM REGIS
65	C THE LAST SEVENTEEN INTERVALS ARE THE BANDS FROM GNF.
66	C
67	C THIS SUBROUTINE SETS THE OPTICAL CONSTANTS IN THE INFRARED
68	C IT CALCUALTES FOR EACH LAYER, FOR EACH SPECRAL INTERVAL IN THE IR
69	C LAYER: WBAR, DTAU, COSBAR
70	C LEVEL: TAU
71	C
72
73	DO 80 K=1,NSPECI
74	TAUHI_1pt(K)=0.
75	TAUCI_1pt(K)=0.
76	TAUGI_1pt(K)=0.
77	80 CONTINUE
78
79	c************************************************************************
80	DO 100 J=1,NLAYER ! BOUCLE SUR L'ALTITUDE
81	c************************************************************************
82
83	c print*,'ig,k,j ',ig,k,j
84
85	C SET UP THE COEFFICIENT TO REDUCE MASS PATH TO STP ...SEE NOTES
86	C T0 =273.15 PO=1.01325 BAR
87
88	TBAR=0.5*(TEMP(J)+TEMP(J+1))
89	PBAR=SQRT(PRESS(J)*PRESS(J+1))
90	BMU=0.5*(XMU(J+1)+XMU(J))
91	COEF1=RGAS273.152.5E5* (PRESS(J+1)2 - PRESS(J)2)
92	& /(1.01325*2 EFFG(Z(J))TBARBMU)
93
94	IF (IPRINT .GT. 9) WRITE(6,21) J,EFFG(Z(J)),TBAR,BMU,COEF1
95	21 FORMAT(' J, EFFG, TBAR, BMU, COEF1,: ',I3,1P6E10.3)
96
97	c------------------------------------------------------------------------
98	DO 101 K=1,NSPECI ! BOUCLE SUR LES L.D'O
99	c------------------------------------------------------------------------
100
101
102	C #1: HAZE
103	C---------------------
104
105	C FIRST COMPUTE TAU AEROSOL
106
107
108	c
109	c /\
110	c / \
111	c / \
112	c / _O \
113	c / \|/ \
114	c / / \ \
115	c / \|\ \/\ \
116	c / \|\| / \ \
117	c ----------------
118	c \| WARNING \|
119	c \| SLOW DOWN \|
120	c ----------------
121
122
123
124
125	c********* EN TRAVAUX *************************
126
127	TAEROS=0.
128	TAEROSCAT=0.
129	CBAR=0.
130
131
132	DO inq=1,nrad !BOUCLE SUR LES TAILLE D"AEROSOLS
133
134
135	IF (WNOI(K).lt.wco) THEN ! lamda > 56 um
136
137	if (iopti.eq.0) then
138
139	c CALL XMIE(rayon(inq)*1.e6,REALI(K),XIMGI(K),
140	c & QEXT,QSCT,QABS,QBAR,WNOI(K))
141
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	CALL OPTFRAC(XMONO,10000./WNOI(K)
195	& ,QEXT,QSCT,QABS,QBAR)
196
197
198	c CALL CFFFV11(1.e-2/WNOI(K),REALI(K),XIMGI(K),RF(inq),2.
199	c & ,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
220
221
222	CBAR=CBAR/TAEROSCAT
223
224	DELTAZ=Z(J)-Z(J+1)
225
226	c --------------------------------------------------------------------
227	c profil brume Pascal: fit T (sauf tropopause) et albedo
228	c -------------------
229	if( cutoff.eq.1) then
230	IF(PRESS(J).gt.9.e-3) THEN
231	TAEROS=TAEROSM1(K)DELTAZ/DELTAZM1(K)0.85
232	TAEROSCAT=TAEROSCATM1(K)DELTAZ/DELTAZM1(K)0.85
233	c TAEROS=0.
234	c TAEROSCAT=0.
235	ENDIF
236
237	IF(PRESS(J).gt.1.e-1) THEN
238	TAEROS=TAEROSM1(K)DELTAZ/DELTAZM1(K)1.15
239	TAEROSCAT=TAEROSCATM1(K)DELTAZ/DELTAZM1(K)1.15
240	c TAEROS=0.
241	c TAEROSCAT=0.
242	ENDIF
243	endif !cutoff=1
244
245	c profil brume pour fit T (y compris tropopause), mais ne fit plus albedo...
246	c -----------------------
247	if( cutoff.eq.2) then
248	IF(PRESS(J).gt.1.e-1) THEN
249	TAEROS=0.
250	TAEROSCAT=0.
251	ENDIF
252	endif !cutoff=2
253	c --------------------------------------------------------------------
254
255	TAEROSM1(K)=TAEROS
256	TAEROSCATM1(K)=TAEROSCAT
257	DELTAZM1(K)=DELTAZ
258
259
260	IF(TAEROSCAT.le.0.) CBAR=0.
261
262	c print*,'HERE, MCKAY AEROSOLS IR'
263	c TAEROS=xv1(j,k)
264	c TAEROSCAT=xv2(j,k)
265	c CBAR=xv3(j,k)
266
267	c if (ig.eq.1) then
268	c if (k.eq.NSPECV/2) then
269	c print*,'@IR',K,J,TAEROS,TAEROSCAT,CBAR
270	c stop'Pour faire des comparaisons'
271	c endif
272	c endif
273
274
275	c********* EN TRAVAUX *************************
276
277	C #2: CLOUD
278	C------------------
279
280	C NEXT COMPUTE TAU CLOUD
281	TAUCLD=0.0
282	CBARC =0.0
283	IF ( XNCLD(J) .GT. 0..and .taufac.gt.0.) THEN
284	print*,'Appel a xmie avec radcld=',radcld(j)
285	CALL XMIE(RADCLD(J),RCLDI(K),XICLDI(K),
286	& QEXTC,QSCTC,QABSC,CBARC,WNOI(K))
287	TAUCLD=QEXTC*XNCLD(J)
288	ENDIF
289
290
291	C #3: GAZ
292	C------------------
293
294	C NOW COMPUTE TAUGAS DUE TO THE PIA TERM ONLY FOR LAMDA LT 940
295	TAUGAS=0.0
296	IF (WNOI(K) .LT. 940. ) THEN
297	c if(ig.eq.1.and.k.eq.nspecv/2) print*,'avant PIA'
298	CALL PIA(K,TBAR,PNN,PCC,PCN,PHN)
299	c if(ig.eq.1.and.k.eq.nspecv/2) print*,'apres PIA'
300	C HERE IS WHERE WE COULD SCALE THE PIA COEFFICEINTS TO FIT DATA
301	C BASED ON REGIS' NOTES. ---TGM HAS THIS ADJUST IN IT AS DEFAULT
302	PCN=PCN*MIN(1.75 , AMAX1(1.0,WNOI(K)/200.))
303	C**REPLACE ABOVE WITH: PCN=PCN1.25*MIN(1.75 , AMAX1(1.0,WNOI(K)/200.))
304	C 1.25 FACTOR (NOT FROM DATA) SUGGESTED BY TOON et al. (1988)
305	TAUGAS=COEF1*
306	& (XN2(J)XN2(J)PNN + CH4(J)CH4(J)PCC
307	& + XN2(J)CH4(J)PCN + XN2(J)H2(J)PHN)
308	IF (J .EQ. NLAYER .AND. IPRINT .GT. 9)
309	& WRITE (6,22) WNOI(K),TAUGAS,XN2(J),CH4(J),H2(J),
310	& TBAR, PNN,PCC,PCN, PHN,
311	& XN2(J)XN2(J)PNN , CH4(J)CH4(J)PCC ,
312	& XN2(J)CH4(J)PCN , XN2(J)H2(J)PHN
313	22 FORMAT(1X,1P8E10.2)
314	ENDIF
315
316	IF (K .GT. 28) THEN
317	KGAS=K-28
318	C ??FLAG? HERE MUST BE WATCHED CAREFULLY
319	U=COLDEN(J)*6.02204E23/BMU
320	if(ig.eq.1.and.k.eq.nspecv/2) print*,'Avant GAS2'
321	if((ylellouch).or.(.not.hcnrad)) then
322	CALL GAS2_NOHCN(J, KGAS,TBAR,PBAR,U,TAU2)
323	else
324	CALL GAS2(J, KGAS,TBAR,PBAR,U,TAU2)
325	endif
326	if(ig.eq.1.and.k.eq.nspecv/2) print*,'Apres GAS2'
327	TAUGAS=TAUGAS+TAU2
328	ENDIF
329	C
330
331	IF (TAEROS + TAUCLD .GT. 0.) THEN
332	COSBI_1pt(J,K)=(CBARTAEROSCAT + CBARCTAUCLD )
333	& /(TAEROSCAT + TAUCLD)
334	ELSE
335	COSBI_1pt(J,K)=0.0
336	ENDIF
337
338	DTAUI_1pt(J,K)=TAUGAS+TAEROS+TAUCLD
339
340	TAUHI_1pt(K)=TAUHI_1pt(K) + TAEROS
341	TAUHID_1pt(J,K)=TAUHI_1pt(K)
342
343	TAUGI_1pt(K)=TAUGI_1pt(K) + TAUGAS
344	TAUGID_1pt(J,K)=TAUGI_1pt(K)
345
346	TAUCI_1pt(K)=TAUCI_1pt(K) + TAUCLD
347
348	c if (ig.eq.1) then
349	c if (k.eq.NSPECI/2.or.k.eq.1.or.k.eq.NSPECI) then
350	c print*,'@IR',K,J,DTAUI_1pt(J,K),TAUGAS,TAEROS,TAUCLD
351	c endif
352	c endif
353
354
355	C ??FLAG? SERIOUS PROBLEM WITH THE CODE HERE!
356
357	TLIMIT=1.E-16
358
359	IF (DTAUI_1pt(J,K) .GT. TLIMIT) THEN
360
361	c***************** ECHANGE
362	c WBARI(J,K)=(QSCTXNUMB(J) + QSCTCXNCLD(J))
363	c****************
364	CFC WBARI_1pt(J,K)=(TAEROSCAT + QSCTC*XNCLD(J))
365	c****************
366	WBARI_1pt(J,K)=TAEROSCAT
367	& /DTAUI_1pt(J,K)
368	c if(iwarning.eq.0)
369	c s print*,'WARNING!!! dans optci xncld non initialise'
370	c iwarning=1
371
372	ELSE
373
374	WBARI_1pt(J,K)=0.0
375	c PRINT*,'gaz ',taugas,'aerosols',taeros,'nuages',taucld
376	c WRITE (6,999) J,K,DTAUI_1pt(J,K)
377	999 FORMAT (' WARNING TAU MINIMUM AT J,K,DTAU:',2I3,1PE10.3)
378	DTAUI_1pt(J,K)=TLIMIT
379
380	ENDIF
381
382	c IF (IPRINT .GT. 9)
383	c & WRITE(6,73)J,K,TAUGAS,TAEROS,QEXT,QSCT
384	73 FORMAT(2I3,1P8E10.3)
385
386
387	c------------------------------------------------------------------------
388	101 CONTINUE ! FIN BOUCLE L D'O
389	c------------------------------------------------------------------------
390
391	iopti=1
392
393	c************************************************************************
394	100 CONTINUE ! FIN BOUCLE ALTITUDE
395	c************************************************************************
396
397	DO 119 K=1,NSPECI
398	TAUI_1pt(1,K)=0.0
399	DO 119 J=1,NLAYER
400	TAUI_1pt(J+1,K)=TAUI_1pt(J,K)+DTAUI_1pt(J,K)
401	119 CONTINUE
402
403	c IF (IPRINT .GT. 2) THEN
404	c WRITE (6,120)
405	c 120 FORMAT(///' OPTICAL CONSTANTS IN THE INFRARED')
406
407	c DO 200 K=1,NSPECI ! #2
408	c WRITE (6,190)
409	c WRITE (6,210)K,WLNI(K),WNOI(K),BWNI(K)
410	c & ,BWNI(K)+DWNI(K),DWNI(K)
411	c WRITE (6,230)REALI(K),XIMGI(K)
412
413	c DO 195 J=1,NLAYER ! #3
414	c WRITE (6,220)XNUMB(J), WBARI_1pt(J,K),COSBI_1pt(J,K)
415	c & , DTAUI_1pt(J,K),TAUI_1pt(J,K)
416	c 195 CONTINUE
417
418	c IF(ig.eq.12) WRITE (6,240) TAUI_1pt(NLEVEL,K)
419	c 200 CONTINUE
420
421	c END IF
422
423
424	c 210 FORMAT(1X,I3,F10.3,F10.2,F10.2,'-',F8.2,F10.3)
425	c 190 FORMAT(1X//' SNUM MICRONS WAVENU INTERVAL DELTA-WN')
426	c 230 FORMAT(1X,'NREAL(LAYER)= ',1PE10.3,' NIMG(LAYER)= ',E10.3/
427	c &' #AEROSOLS WBAR COSBAR DTAU TAU')
428	c 220 FORMAT(5(1X,G9.3))
429	c 240 FORMAT(41X,G9.3)
430
431	RETURN
432	END

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

Download in other formats:

Original Format