Context Navigation

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
8	c Arguments:
9	c ---------
10	REAL ykim(klon,klev,nqmx)
11	integer nmicro
12	c ---------
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
75	C THE PRESSURE INDUCED TRANSITIONS ARE FROM REGIS
76	C THE LAST SEVENTEEN INTERVALS ARE THE BANDS FROM GNF.
77	C
78	C THIS SUBROUTINE SETS THE OPTICAL CONSTANTS IN THE INFRARED
79	C IT CALCUALTES FOR EACH LAYER, FOR EACH SPECRAL INTERVAL IN THE IR
80	C LAYER: WBAR, DTAU, COSBAR
81	C LEVEL: TAU
82	C
83	print*,'START OPTCI'
84	c sum=0.
85	c do nng=2,klon
86	c do i=1,klev
87	c do j=1,nqmx
88	c print*,'j,rj',j,r(j)
89	c print*,'paer',qaer(nng,i,j)
90	c sum=sum+qaer(nng,i,j)r(j)3.1.33333.14151000.
91	c enddo
92	c enddo
93	c enddo
94	c print*,sum/(klon-1),'SOMME COLONNE/OPTCI'
95
96
97	c do inq=1,nqmx
98	c print*,inq,r(inq),vrat,qaer(12,25,inq)
99	c enddo
100
101	DO 420 K=1,NSPECI
102	C 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
106	C SET UP THE OPTICAL CONSTANTS FOR THE CLOUD
107	RCLDI(K)=1.27
108	XICLDI(K)=REFLIQ(WNOI(K))
109	420 CONTINUE
110	C
111	C
112	C ZERO ALL OPTICAL DEPTHS.
113	C ??FLAG? FOR SOME APPLCIATIONS THE TOP OPACITY MAY NOT VANISH
114
115	c open (unit=1,file='xsetupi')
116	c do i=1,klev
117	c read(1,*) a
118	c do j=1,nspeci
119	c read(1,*) xv1(i,j),xv2(i,j),xv3(i,j)
120	c enddo
121	c enddo
122	c close(1)
123
124	c************************************************************************
125	c************************************************************************
126	DO 79 ig=1,klon ! BOUCLE SUR GRILLE HORIZONTALE
127	c print*,'ig NEW optci',ig
128	c************************************************************************
129	c************************************************************************
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
160	c************************************************************************
161	DO 100 J=1,NLAYER ! BOUCLE SUR L'ALTITUDE
162	c************************************************************************
163
164	c print*,'ig,k,j ',ig,k,j
165
166	C SET UP THE COEFFICIENT TO REDUCE MASS PATH TO STP ...SEE NOTES
167	C 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=RGAS273.152.5E5* (PRESS(J+1)2 - PRESS(J)2)
173	& /(1.01325*2 EFFG(Z(J))TBARBMU)
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
178	c------------------------------------------------------------------------
179	DO 101 K=1,NSPECI ! BOUCLE SUR LES L.D'O
180	c------------------------------------------------------------------------
181
182
183	C #1: HAZE
184	C---------------------
185
186	C FIRST COMPUTE TAU AEROSOL
187
188
189	c
190	c /\
191	c / \
192	c / \
193	c / _O \
194	c / \|/ \
195	c / / \ \
196	c / \|\ \/\ \
197	c / \|\| / \ \
198	c ----------------
199	c \| WARNING \|
200	c \| SLOW DOWN \|
201	c ----------------
202
203
204
205
206	c********* 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
219	c CALL XMIE(R(inq)*1.e6,REALI(K),XIMGI(K),
220	c & 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
292	c CALL CFFFV11(1.e-2/WNOI(K),REALI(K),XIMGI(K),RF(inq),2.
293	c & ,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
327	c --------------------------------------------------------------------
328	c profil brume Pascal: fit T (sauf tropopause) et albedo
329	c -------------------
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
334	c TAEROS=0.
335	c 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
341	c TAEROS=0.
342	c TAEROSCAT=0.
343	ENDIF
344	endif !cutoff=1
345
346	c profil brume pour fit T (y compris tropopause), mais ne fit plus albedo...
347	c -----------------------
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
354	c --------------------------------------------------------------------
355
356	TAEROSM1(K)=TAEROS
357	TAEROSCATM1(K)=TAEROSCAT
358	DELTAZM1(K)=DELTAZ
359
360
361	IF(TAEROSCAT.le.0.) CBAR=0.
362
363	c print*,'HERE, MCKAY AEROSOLS IR'
364	c TAEROS=xv1(j,k)
365	c TAEROSCAT=xv2(j,k)
366	c CBAR=xv3(j,k)
367
368	c if (ig.eq.1) then
369	c if (k.eq.NSPECV/2) then
370	c print*,'@IR',K,J,TAEROS,TAEROSCAT,CBAR
371	c stop'Pour faire des comparaisons'
372	c endif
373	c endif
374
375
376	c********* EN TRAVAUX *************************
377
378	C #2: CLOUD
379	C------------------
380
381	C 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
391	C #3: GAZ
392	C------------------
393
394	C NOW COMPUTE TAUGAS DUE TO THE PIA TERM ONLY FOR LAMDA LT 940
395	TAUGAS=0.0
396	IF (WNOI(K) .LT. 940. ) THEN
397	c if(ig.eq.1.and.k.eq.nspecv/2) print*,'avant PIA'
398	CALL PIA(K,TBAR,PNN,PCC,PCN,PHN)
399	c if(ig.eq.1.and.k.eq.nspecv/2) print*,'apres PIA'
400	C HERE IS WHERE WE COULD SCALE THE PIA COEFFICEINTS TO FIT DATA
401	C BASED ON REGIS' NOTES. ---TGM HAS THIS ADJUST IN IT AS DEFAULT
402	PCN=PCN*MIN(1.75 , AMAX1(1.0,WNOI(K)/200.))
403	C**REPLACE ABOVE WITH: PCN=PCN1.25*MIN(1.75 , AMAX1(1.0,WNOI(K)/200.))
404	C 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
418	C ??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
429	C
430
431	IF (TAEROS + TAUCLD .GT. 0.) THEN
432	COSBI(ig,J,K)=(CBARTAEROSCAT + CBARCTAUCLD )
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
448	c if (ig.eq.1) then
449	c if (k.eq.NSPECI/2.or.k.eq.1.or.k.eq.NSPECI) then
450	c print*,'@IR',K,J,DTAUI(ig,J,K),TAUGAS,TAEROS,TAUCLD
451	c endif
452	c endif
453
454	C ??FLAG? SERIOUS PROBLEM WITH THE CODE HERE!
455
456	TLIMIT=1.E-16
457
458	IF (DTAUI(ig,J,K) .GT. TLIMIT) THEN
459
460	c***************** ECHANGE
461	c WBARI(J,K)=(QSCTXNUMB(J) + QSCTCXNCLD(J))
462	c****************
463	CFC WBARI(ig,J,K)=(TAEROSCAT + QSCTC*XNCLD(J))
464	c****************
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
474	c PRINT*,'gaz ',taugas,'aerosols',taeros,'nuages',taucld
475	c 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
481	c IF (IPRINT .GT. 9)
482	c & WRITE(6,73)J,K,TAUGAS,TAEROS,QEXT,QSCT
483	73 FORMAT(2I3,1P8E10.3)
484
485
486	c------------------------------------------------------------------------
487	101 CONTINUE ! FIN BOUCLE L D'O
488	c------------------------------------------------------------------------
489
490	iopti=1
491
492	c************************************************************************
493	100 CONTINUE ! FIN BOUCLE ALTITUDE
494	c************************************************************************
495
496	c************************************************************************
497	c************************************************************************
498	79 CONTINUE ! FIN BOUCLE GRILLE HORIZONTALE
499	c************************************************************************
500	c************************************************************************
501
502
503	C 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
514	c WRITE (6,120)
515	120 FORMAT(///' OPTICAL CONSTANTS IN THE INFRARED')
516
517	c DO 200 K=1,NSPECI ! #2
518	c WRITE (6,190)
519	c WRITE (6,210)K,WLNI(K),WNOI(K),BWNI(K)
520	c & ,BWNI(K)+DWNI(K),DWNI(K)
521	c WRITE (6,230)REALI(K),XIMGI(K)
522
523	c DO 195 J=1,NLAYER ! #3
524	c WRITE (6,220)XNUMB(J), WBARI(ig,J,K),COSBI(ig,J,K)
525	c & , DTAUI(ig,J,K),TAUI(ig,J,K)
526	c 195 CONTINUE
527
528	c IF(ig.eq.12) WRITE (6,240) TAUI(ig,NLEVEL,K)
529	c 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
545	c iopti=1
546	print*, 'FIN OPTCI'
547	stop
548	RETURN
549	END

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

Download in other formats:

Original Format