Context Navigation

optci_1pt.F @ 779

Last change on this file since 779 was 495, checked in by slebonnois, 13 years ago
Mise a jour physique Titan, ajout des forces de marees (dans la dynamique, sous flag titan). SL.
File size: 13.9 KB

Line
1	SUBROUTINE optci_1pt(zqaer_1pt,rcdb,xfrb,iopti,IPRINT)
2	use dimphy
3	#include "dimensions.h"
4	#include "microtab.h"
5	#include "numchimrad.h"
6	#include "clesphys.h"
7
8	PARAMETER(NLAYER=llm,NLEVEL=NLAYER+1)
9	PARAMETER (NSPECI=46,NSPC1I=47,NSPECV=24,NSPC1V=25)
10
11	c Arguments:
12	c ---------
13	integer IPRINT,iopti
14	C iopti: premier appel, on ne calcule qu'une fois les QM et QF
15	* nrad dans microtab.h
16	real zqaer_1pt(NLAYER,2*nrad)
17	#include "optci_1pt.h"
18	c ---------
19
20	COMMON /ATM/ Z(NLEVEL),PRESS(NLEVEL),DEN(NLEVEL),TEMP(NLEVEL)
21
22	COMMON /GASS/ CH4(NLEVEL),XN2(NLEVEL),H2(NLEVEL),AR(NLEVEL)
23	& ,XMU(NLEVEL),GAS1(NLAYER),COLDEN(NLAYER)
24
25	COMMON /STRATO/ C2H2(NLAYER),C2H6(NLAYER)
26	COMMON /STRAT2/ HCN(NLAYER)
27
28	COMMON /AERSOL/ RADIUS(NLAYER), XNUMB(NLAYER)
29	& , REALI(NSPECI), XIMGI(NSPECI), REALV(NSPECV), XIMGV(NSPECV)
30
31	COMMON /CLOUD/
32	& RCLDI(NSPECI), XICLDI(NSPECI)
33	& , RCLDV(NSPECV), XICLDV(NSPECV)
34	& , RCLDI2(NSPECI), XICLDI2(NSPECI)
35	& , RCLDV2(NSPECV), XICLDV2(NSPECV)
36
37	COMMON /SPECTI/ BWNI(NSPC1I), WNOI(NSPECI),
38	& DWNI(NSPECI), WLNI(NSPECI)
39
40	COMMON /PLANT/ CSUBP,F0PI
41	COMMON /ADJUST/ RHCH4,FH2,FHAZE,FHVIS,FHIR,TAUFAC,RCLOUD,FARGON
42	COMMON /CONST/RGAS,RHOP,PI,SIGMA
43	COMMON /part/v,rayon,vrat,dr,dv
44
45	DIMENSION PROD(NLEVEL)
46	* nrad dans microtab.h
47	real v(nrad),rayon(nrad),vrat,dr(nrad),dv(nrad)
48	real xv1(klev,nspeci),xv2(klev,nspeci)
49	real xv3(klev,nspeci)
50	REAL QF1(nrad,NSPECI),QF2(nrad,NSPECI)
51	REAL QF3(nrad,NSPECI),QF4(nrad,NSPECI)
52	REAL QM1(nrad,NSPECI),QM2(nrad,NSPECI)
53	REAL QM3(nrad,NSPECI),QM4(nrad,NSPECI)
54	REAL QC1(nrad,NSPECI),QC2(nrad,NSPECI)
55	REAL QC3(nrad,NSPECI),QC4(nrad,NSPECI)
56	real emu
57	REAL TAEROSM1(NSPECI),TAEROSCATM1(NSPECI),DELTAZM1(NSPECI)
58
59	c ---- nuages
60	REAL TNUAGE,TNUAGESCAT
61	REAL rcdb(nlayer),xfrb(nlayer,4)
62
63	save qf1,qf2,qf3,qf4,qm1,qm2,qm3,qm4
64
65	C THE PRESSURE INDUCED TRANSITIONS ARE FROM REGIS
66	C THE LAST SEVENTEEN INTERVALS ARE THE BANDS FROM GNF.
67	C
68	C THIS SUBROUTINE SETS THE OPTICAL CONSTANTS IN THE INFRARED
69	C IT CALCUALTES FOR EACH LAYER, FOR EACH SPECRAL INTERVAL IN THE IR
70	C LAYER: WBAR, DTAU, COSBAR
71	C LEVEL: TAU
72	C
73
74	DO 80 K=1,NSPECI
75	TAUHI_1pt(K)=0.
76	TAUCI_1pt(K)=0.
77	TAUGI_1pt(K)=0.
78	80 CONTINUE
79
80	c************************************************************************
81	DO 100 J=1,NLAYER ! BOUCLE SUR L'ALTITUDE
82	c************************************************************************
83
84	c print*,'ig,k,j ',ig,k,j
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
144	CALL CMIE(1.E-2/WNOI(K),REALI(K),XIMGI(K),rayon(inq),
145	& QEXT,QSCT,QABS,QBAR)
146
147
148	QM1(inq,K)=QEXT
149	QM2(inq,K)=QSCT
150	QM3(inq,K)=QABS
151	QM4(inq,K)=QBAR
152
153	endif ! end iopti
154
155
156	TAEROS=QM1(inq,K)zqaer_1pt(nlayer+1-J,inq)1.e-4+TAEROS
157	TAEROSCAT=QM2(inq,K)zqaer_1pt(nlayer+1-J,inq)1.e-4+TAEROSCAT
158	CBAR=CBAR+QM4(inq,K)QM2(inq,K)zqaer_1pt(nlayer+1-J,inq)*1.e-4
159
160
161	ELSE ! 0.2 < lambda < 56 um
162
163
164	if(rayon(inq).lt.RF(inq)) THEN
165
166	if (iopti.eq.0) then
167
168	CALL XMIE(rayon(inq)*1.e6,REALI(K),XIMGI(K),
169	& QEXT,QSCT,QABS,QBAR,WNOI(K))
170
171	QM1(inq,K)=QEXT
172	QM2(inq,K)=QSCT
173	QM3(inq,K)=QABS
174	QM4(inq,K)=QBAR
175	endif ! end iopti
176
177
178	TAEROS=QM1(inq,K)zqaer_1pt(nlayer+1-J,inq)1.e-4+TAEROS
179	TAEROSCAT=QM2(inq,K)zqaer_1pt(nlayer+1-J,inq)1.e-4+TAEROSCAT
180	CBAR=CBAR+QM4(inq,K)QM2(inq,K)zqaer_1pt(nlayer+1-J,inq)*1.e-4
181
182	else
183
184	XMONO=(rayon(inq)/RF(inq))**3.
185	XRULE=1.
186
187	if(XMONO.gt.16384./1.5) then
188	XRULE=(XMONO/16384.)
189	XMONO=16384.
190	endif
191
192
193	if (iopti.eq.0) then
194
195	c CALL OPTFRAC(XMONO,10000./WNOI(K)
196	c & ,QEXT,QSCT,QABS,QBAR)
197
198
199	CALL CFFFV11(1.e-2/WNOI(K),REALI(K),XIMGI(K),RF(inq),2.
200	& ,XMONO,QSCT,QEXT,QABS,QBAR)
201
202
203	QF1(inq,K)=QEXT*XRULE
204	QF2(inq,K)=QSCT*XRULE
205	QF3(inq,K)=QABS*XRULE
206	QF4(inq,K)=QBAR
207	endif ! end iopti
208
209	TAEROS=QF1(inq,K)*zqaer_1pt(nlayer+1-J,inq)+TAEROS
210	TAEROSCAT=QF2(inq,K)*zqaer_1pt(nlayer+1-J,inq)+TAEROSCAT
211	CBAR=CBAR+QF4(inq,K)QF2(inq,K)zqaer_1pt(nlayer+1-J,inq)
212
213	endif
214
215	IF(TAEROS.LT.1.e-10) TAEROS=1.e-10
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	IF (clouds.eq.0) THEN
282	TNUAGE=0.
283	TNUAGESCAT=0.
284	CNBAR=0.
285	ELSE
286	TNUAGE=0.
287	TNUAGESCAT=0.
288	CNBAR=0.
289	QEXTC=0.
290	QSCTC=0.
291	QABSC=0.
292	CBARC=0.
293
294	DO inq=1,nrad !BOUCLE SUR LES NQMX TAILLE D"AEROSOLS
295	QC1(inq,K)=0.
296	QC2(inq,K)=0.
297	QC3(inq,K)=0.
298	QC4(inq,K)=0.
299	ENDDO
300
301	** OPTICAL CONSTANT : MIXING RULES
302
303	IF (rcdb(nlayer+1-J).gt.1.1e-10) THEN
304
305	XNR=xfrb(nlayer+1-J,1) *REALI(K)
306	& +xfrb(nlayer+1-J,2) *RCLDI(K)
307	& +xfrb(nlayer+1-J,3) *RCLDI2(K)
308	& +xfrb(nlayer+1-J,4) *RCLDI2(K)
309
310	XNI=xfrb(nlayer+1-J,1) *XIMGI(K)
311	& +xfrb(nlayer+1-J,2) *XICLDI(K)
312	& +xfrb(nlayer+1-J,3) *XICLDI2(K)
313	& +xfrb(nlayer+1-J,4) *XICLDI2(K)
314
315	** OPTICAL CONSTANT : LIQUID DROP = THOLIN
316
317	IF(xfrb(nlayer+1-J,1).ge.0.1) THEN
318	XNI=XIMGI(K)
319	XNR=REALI(K)
320	ENDIF
321
322	IF (XNI.gt.1.e-10 .and. XNR.gt.1.00) THEN
323	CALL CMIE(1.E-2/WNOI(K),XNR,XNI,
324	& rcdb(nlayer+1-J),
325	& QEXTC,QSCTC,QABSC,CBARC)
326	ELSE
327	PRINT*,' WARNING XNR/XNI in optci: ',XNR,XNI
328	QEXTC=0.
329	QSCTC=0.
330	QABSC=0.
331	CBARC=0.
332	ENDIF
333	ELSE
334	QEXTC=0.
335	QSCTC=0.
336	QABSC=0.
337	CBARC=0.
338	ENDIF
339
340	DO inq=1,nrad !BOUCLE SUR LES NQMX TAILLE D"AEROSOLS
341	QC1(inq,K)=QEXTC/xnuf
342	QC2(inq,K)=QSCTC/xnuf
343	QC3(inq,K)=QABSC/xnuf
344	QC4(inq,K)=CBARC
345	TNUAGE=QC1(inq,K)zqaer_1pt(nlayer+1-J,inq+nrad)1.e-4
346	& +TNUAGE
347	TNUAGESCAT=QC2(inq,K)zqaer_1pt(nlayer+1-J,inq+nrad)1.e-4
348	& +TNUAGESCAT
349	CNBAR=QC4(inq,K)*QC2(inq,K)
350	& zqaer_1pt(nlayer+1-J,inq+nrad)1.e-4+CNBAR
351	ENDDO
352
353	IF(TNUAGESCAT.EQ.0.) THEN
354	CNBAR=0.
355	ELSE
356	CNBAR=CNBAR/TNUAGESCAT
357	ENDIF
358	ENDIF ! Cond CLD
359
360
361	C #3: GAZ
362	C------------------
363
364	C NOW COMPUTE TAUGAS DUE TO THE PIA TERM ONLY FOR LAMDA LT 940
365	TAUGAS=0.0
366	IF (WNOI(K) .LT. 940. ) THEN
367	c if(ig.eq.1.and.k.eq.nspecv/2) print*,'avant PIA'
368	CALL PIA(K,TBAR,PNN,PCC,PCN,PHN)
369	c if(ig.eq.1.and.k.eq.nspecv/2) print*,'apres PIA'
370	C HERE IS WHERE WE COULD SCALE THE PIA COEFFICEINTS TO FIT DATA
371	C BASED ON REGIS' NOTES. ---TGM HAS THIS ADJUST IN IT AS DEFAULT
372	PCN=PCN*MIN(1.75 , AMAX1(1.0,WNOI(K)/200.))
373	C**REPLACE ABOVE WITH: PCN=PCN1.25*MIN(1.75 , AMAX1(1.0,WNOI(K)/200.))
374	C 1.25 FACTOR (NOT FROM DATA) SUGGESTED BY TOON et al. (1988)
375	TAUGAS=COEF1*
376	& (XN2(J)XN2(J)PNN + CH4(J)CH4(J)PCC
377	& + XN2(J)CH4(J)PCN + XN2(J)H2(J)PHN)
378	IF (J .EQ. NLAYER .AND. IPRINT .GT. 9)
379	& WRITE (6,22) WNOI(K),TAUGAS,XN2(J),CH4(J),H2(J),
380	& TBAR, PNN,PCC,PCN, PHN,
381	& XN2(J)XN2(J)PNN , CH4(J)CH4(J)PCC ,
382	& XN2(J)CH4(J)PCN , XN2(J)H2(J)PHN
383	22 FORMAT(1X,1P8E10.2)
384	ENDIF
385
386	IF (K .GT. 28) THEN
387	KGAS=K-28
388	C ??FLAG? HERE MUST BE WATCHED CAREFULLY
389	U=COLDEN(J)*6.02204E23/BMU
390	c if(ig.eq.1.and.k.eq.nspecv/2) print*,'Avant GAS2'
391	if((ylellouch).or.(.not.hcnrad)) then
392	CALL GAS2_NOHCN(J, KGAS,TBAR,PBAR,U,TAU2)
393	else
394	CALL GAS2(J, KGAS,TBAR,PBAR,U,TAU2)
395	endif
396	c if(ig.eq.1.and.k.eq.nspecv/2) print*,'Apres GAS2'
397	TAUGAS=TAUGAS+TAU2
398	ENDIF
399	C
400
401	DTAUI_1pt(J,K)=TAUGAS+TAEROS+TNUAGE
402	DTAUIP_1pt(J,K)=TAUGAS+TAEROS
403
404	TAUHI_1pt(K)=TAUHI_1pt(K) + TAEROS
405	TAUHID_1pt(J,K)=TAUHI_1pt(K)
406
407	TAUGI_1pt(K)=TAUGI_1pt(K) + TAUGAS
408	TAUGID_1pt(J,K)=TAUGI_1pt(K)
409
410	TAUCI_1pt(K)=TAUCI_1pt(K) + TNUAGE
411	TAUCID_1pt(J,K)=TAUCI_1pt(K)
412
413	C ??FLAG? SERIOUS PROBLEM WITH THE CODE HERE!
414
415	TLIMIT=1.E-16
416
417
418	IF (TAEROSCAT + TNUAGESCAT .GT. 0.) THEN
419	COSBI_1pt(J,K)=(CBARTAEROSCAT + CNBARTNUAGESCAT )
420	& /(TAEROSCAT + TNUAGESCAT)
421	ELSE
422	COSBI_1pt(J,K)=0.0
423	ENDIF
424
425	IF (TAEROSCAT .GT. 0.) THEN
426	COSBIP_1pt(J,K)=(CBAR*TAEROSCAT)
427	& /(TAEROSCAT)
428	ELSE
429	COSBIP_1pt(J,K)=0.0
430	ENDIF
431
432	*---------
433
434	IF (DTAUI_1pt(J,K) .GT. TLIMIT) THEN
435	WBARI_1pt(J,K)=(TAEROSCAT+TNUAGESCAT) /DTAUI_1pt(J,K)
436	ELSE
437	WBARI_1pt(J,K)=0.0
438	DTAUI_1pt(J,K)=TLIMIT
439	ENDIF
440
441	IF (DTAUIP_1pt(J,K) .GT. TLIMIT) THEN
442	WBARIP_1pt(J,K)=(TAEROSCAT) /DTAUIP_1pt(J,K)
443	ELSE
444	WBARIP_1pt(J,K)=0.0
445	DTAUIP_1pt(J,K)=TLIMIT
446	ENDIF
447
448
449	c IF (IPRINT .GT. 9)
450	c & WRITE(6,73)J,K,TAUGAS,TAEROS,QEXT,QSCT
451	73 FORMAT(2I3,1P8E10.3)
452
453
454	c------------------------------------------------------------------------
455	101 CONTINUE ! FIN BOUCLE L D'O
456	c------------------------------------------------------------------------
457
458	iopti=1
459
460	c************************************************************************
461	100 CONTINUE ! FIN BOUCLE ALTITUDE
462	c************************************************************************
463
464	DO 119 K=1,NSPECI
465	TAUI_1pt(1,K)=0.0
466	TAUIP_1pt(1,K)=0.0
467	DO 119 J=1,NLAYER
468	TAUI_1pt(J+1,K)=TAUI_1pt(J,K)+DTAUI_1pt(J,K)
469	TAUIP_1pt(J+1,K)=TAUIP_1pt(J,K)+DTAUIP_1pt(J,K)
470	119 CONTINUE
471
472	c IF (IPRINT .GT. 2) THEN
473	c WRITE (6,120)
474	c 120 FORMAT(///' OPTICAL CONSTANTS IN THE INFRARED')
475
476	c DO 200 K=1,NSPECI ! #2
477	c WRITE (6,190)
478	c WRITE (6,210)K,WLNI(K),WNOI(K),BWNI(K)
479	c & ,BWNI(K)+DWNI(K),DWNI(K)
480	c WRITE (6,230)REALI(K),XIMGI(K)
481
482	c DO 195 J=1,NLAYER ! #3
483	c WRITE (6,220)XNUMB(J), WBARI_1pt(J,K),COSBI_1pt(J,K)
484	c & , DTAUI_1pt(J,K),TAUI_1pt(J,K)
485	c 195 CONTINUE
486
487	c IF(ig.eq.12) WRITE (6,240) TAUI_1pt(NLEVEL,K)
488	c 200 CONTINUE
489
490	c END IF
491
492
493	c 210 FORMAT(1X,I3,F10.3,F10.2,F10.2,'-',F8.2,F10.3)
494	c 190 FORMAT(1X//' SNUM MICRONS WAVENU INTERVAL DELTA-WN')
495	c 230 FORMAT(1X,'NREAL(LAYER)= ',1PE10.3,' NIMG(LAYER)= ',E10.3/
496	c &' #AEROSOLS WBAR COSBAR DTAU TAU')
497	c 220 FORMAT(5(1X,G9.3))
498	c 240 FORMAT(41X,G9.3)
499
500	RETURN
501	END

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

Download in other formats:

Original Format