Context Navigation

optci_1pt_2.F @ 696

Last change on this file since 696 was 495, checked in by slebonnois, 14 years ago
Mise a jour physique Titan, ajout des forces de marees (dans la dynamique, sous flag titan). SL.
File size: 13.6 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,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	CALL CFFFV11(1.e-2/WNOI(K),REALI(K),XIMGI(K),RF(inq),2.
195	& ,XMONO,QSCT,QEXT,QABS,QBAR)
196
197	QF1(inq,K)=QEXT*XRULE
198	QF2(inq,K)=QSCT*XRULE
199	QF3(inq,K)=QABS*XRULE
200	QF4(inq,K)=QBAR
201	endif ! end iopti
202
203	TAEROS=QF1(inq,K)*zqaer_1pt(nlayer+1-J,inq)+TAEROS
204	TAEROSCAT=QF2(inq,K)*zqaer_1pt(nlayer+1-J,inq)+TAEROSCAT
205	CBAR=CBAR+QF4(inq,K)QF2(inq,K)zqaer_1pt(nlayer+1-J,inq)
206
207	endif
208
209	IF(TAEROS.LT.1.e-10) TAEROS=1.e-10
210
211	ENDIF
212	ENDDO ! FIN DE LA BOUCLE SUR nrad
213
214	IF(TAEROSCAT.le.0.) then
215	CBAR=0.
216	ELSE
217	CBAR=CBAR/TAEROSCAT
218	ENDIF
219
220	DELTAZ=Z(J)-Z(J+1)
221
222	c --------------------------------------------------------------------
223	c profil brume Pascal: fit T (sauf tropopause) et albedo
224	c -------------------
225	if( cutoff.eq.1) then
226	IF(PRESS(J).gt.9.e-3) THEN
227	TAEROS=TAEROSM1(K)DELTAZ/DELTAZM1(K)0.85
228	TAEROSCAT=TAEROSCATM1(K)DELTAZ/DELTAZM1(K)0.85
229	c TAEROS=0.
230	c TAEROSCAT=0.
231	ENDIF
232
233	IF(PRESS(J).gt.1.e-1) THEN
234	TAEROS=TAEROSM1(K)DELTAZ/DELTAZM1(K)1.15
235	TAEROSCAT=TAEROSCATM1(K)DELTAZ/DELTAZM1(K)1.15
236	c TAEROS=0.
237	c TAEROSCAT=0.
238	ENDIF
239	endif !cutoff=1
240
241	c profil brume pour fit T (y compris tropopause), mais ne fit plus albedo...
242	c -----------------------
243	if( cutoff.eq.2) then
244	IF(PRESS(J).gt.1.e-1) THEN
245	TAEROS=0.
246	TAEROSCAT=0.
247	ENDIF
248	endif !cutoff=2
249	c --------------------------------------------------------------------
250
251	TAEROSM1(K)=TAEROS
252	TAEROSCATM1(K)=TAEROSCAT
253	DELTAZM1(K)=DELTAZ
254
255
256	IF(TAEROSCAT.le.0.) CBAR=0.
257
258	c print*,'HERE, MCKAY AEROSOLS IR'
259	c TAEROS=xv1(j,k)
260	c TAEROSCAT=xv2(j,k)
261	c CBAR=xv3(j,k)
262
263	c********* EN TRAVAUX *************************
264
265	C #2: CLOUD
266	C------------------
267	C NEXT COMPUTE TAU CLOUD
268	c
269	c Menu special :
270	c Afin d'eviter la surcharge de calcul on ne calcule les
271	c propriétes optiques des nuages qu'une seule fois
272	c avec un rayon de particule effectif de 3um et une composition
273	c de goutte : 90% CH4 / 10% NOYAUX
274	c Puis on ajute les section efficace par la surface reelle de
275	c la goutte.
276	c
277	c ---> A TESTER !!!!
278	c
279	TNUAGE=0.
280	TNUAGESCAT=0.
281	CNBAR=0.
282	IF (clouds.eq.1) THEN
283	IF (iopti.eq.0) THEN !--> au premier appel
284	QEXTC=0.
285	QSCTC=0.
286	QABSC=0.
287	CBARC=0.
288	DO inq=1,nrad !BOUCLE SUR LES NQMX TAILLE D"AEROSOLS
289	QC1(inq,K)=0.
290	QC2(inq,K)=0.
291	QC3(inq,K)=0.
292	QC4(inq,K)=0.
293	ENDDO
294	** OPTICAL CONSTANT : MIXING RULES
295	** Fraction volumique fixe :
296	** 10% noyaux.
297	** 90% methane.
298	XNR = 0.5 * REALI(K)
299	& + 0.5 * RCLDI(K)
300	XNI = 0.5 * XIMGI(K)
301	& + 0.5 * XICLDI(K)
302	**
303	** Efficacite : particule de 3um de rayon
304	CALL CMIE(1.E-2/WNOI(K),XNR,XNI,3.e-6,
305	& QEXTC,QSCTC,QABSC,CBARC)
306	**
307	** ATTENTION CE SONT DES EFFICACITES : il faut les x par la surface REELLE de la goutte.
308	DO inq=1,nrad
309	QC1(inq,K)=QEXTC/xnuf
310	QC2(inq,K)=QSCTC/xnuf
311	QC3(inq,K)=QABSC/xnuf
312	QC4(inq,K)=CBARC
313	ENDDO
314	ENDIF ! iopti = 0
315
316	c ----- On ne calcule les constante optiques que si Rgoutte > 1e-10
317	IF (rcdb(nlayer+1-J).gt.1.1e-10) THEN
318	DO inq=1,nrad
319	TNUAGE=QC1(inq,K)(rcdb(nlayer+1-J)/3.e-6)2.1.e-4*
320	& zqaer_1pt(nlayer+1-J,inq+nrad) +
321	& TNUAGE
322	TNUAGESCAT=QC2(inq,K)(rcdb(nlayer+1-J)/3.e-6)2.
323	& 1.e-4*zqaer_1pt(nlayer+1-J,inq+nrad) +
324	& TNUAGESCAT
325	CNBAR=QC4(inq,K)QC2(inq,K)
326	& (rcdb(nlayer+1-J)/3.e-6)*2.
327	& 1.e-4*zqaer_1pt(nlayer+1-J,inq+nrad) +
328	& CNBAR
329	ENDDO
330	ENDIF
331
332	IF(TNUAGESCAT.EQ.0.) THEN
333	CNBAR=0.
334	ELSE
335	CNBAR=CNBAR/TNUAGESCAT
336	ENDIF
337
338	ENDIF ! Cond CLD
339	c
340	C #3: GAZ
341	C------------------
342
343	C NOW COMPUTE TAUGAS DUE TO THE PIA TERM ONLY FOR LAMDA LT 940
344	TAUGAS=0.0
345	IF (WNOI(K) .LT. 940. ) THEN
346	CALL PIA(K,TBAR,PNN,PCC,PCN,PHN)
347	C HERE IS WHERE WE COULD SCALE THE PIA COEFFICEINTS TO FIT DATA
348	C BASED ON REGIS' NOTES. ---TGM HAS THIS ADJUST IN IT AS DEFAULT
349	PCN=PCN*MIN(1.75 , AMAX1(1.0,WNOI(K)/200.))
350	C**REPLACE ABOVE WITH: PCN=PCN1.25*MIN(1.75 , AMAX1(1.0,WNOI(K)/200.))
351	C 1.25 FACTOR (NOT FROM DATA) SUGGESTED BY TOON et al. (1988)
352	TAUGAS=COEF1*
353	& (XN2(J)XN2(J)PNN + CH4(J)CH4(J)PCC
354	& + XN2(J)CH4(J)PCN + XN2(J)H2(J)PHN)
355	IF (J .EQ. NLAYER .AND. IPRINT .GT. 9)
356	& WRITE (6,22) WNOI(K),TAUGAS,XN2(J),CH4(J),H2(J),
357	& TBAR, PNN,PCC,PCN, PHN,
358	& XN2(J)XN2(J)PNN , CH4(J)CH4(J)PCC ,
359	& XN2(J)CH4(J)PCN , XN2(J)H2(J)PHN
360	22 FORMAT(1X,1P8E10.2)
361	ENDIF
362
363	IF (K .GT. 28) THEN
364	KGAS=K-28
365	C ??FLAG? HERE MUST BE WATCHED CAREFULLY
366	U=COLDEN(J)*6.02204E23/BMU
367	if((ylellouch).or.(.not.hcnrad)) then
368	CALL GAS2_NOHCN(J, KGAS,TBAR,PBAR,U,TAU2)
369	else
370	CALL GAS2(J, KGAS,TBAR,PBAR,U,TAU2)
371	endif
372	TAUGAS=TAUGAS+TAU2
373	ENDIF
374	C
375
376	DTAUI_1pt(J,K)=TAUGAS+TAEROS+TNUAGE
377	DTAUIP_1pt(J,K)=TAUGAS+TAEROS
378
379	TAUHI_1pt(K)=TAUHI_1pt(K) + TAEROS
380	TAUHID_1pt(J,K)=TAUHI_1pt(K)
381
382	TAUGI_1pt(K)=TAUGI_1pt(K) + TAUGAS
383	TAUGID_1pt(J,K)=TAUGI_1pt(K)
384
385	TAUCI_1pt(K)=TAUCI_1pt(K) + TNUAGE
386	TAUCID_1pt(J,K)=TAUCI_1pt(K)
387
388	C ??FLAG? SERIOUS PROBLEM WITH THE CODE HERE!
389
390	TLIMIT=1.E-16
391
392
393	IF (TAEROSCAT + TNUAGESCAT .GT. 0.) THEN
394	COSBI_1pt(J,K)=(CBARTAEROSCAT + CNBARTNUAGESCAT )
395	& /(TAEROSCAT + TNUAGESCAT)
396	ELSE
397	COSBI_1pt(J,K)=0.0
398	ENDIF
399
400	IF (TAEROSCAT .GT. 0.) THEN
401	COSBIP_1pt(J,K)=(CBAR*TAEROSCAT)
402	& /(TAEROSCAT)
403	ELSE
404	COSBIP_1pt(J,K)=0.0
405	ENDIF
406
407	*---------
408
409	IF (DTAUI_1pt(J,K) .GT. TLIMIT) THEN
410	WBARI_1pt(J,K)=(TAEROSCAT+TNUAGESCAT) /DTAUI_1pt(J,K)
411	ELSE
412	WBARI_1pt(J,K)=0.0
413	DTAUI_1pt(J,K)=TLIMIT
414	ENDIF
415
416	IF (DTAUIP_1pt(J,K) .GT. TLIMIT) THEN
417	WBARIP_1pt(J,K)=(TAEROSCAT) /DTAUIP_1pt(J,K)
418	ELSE
419	WBARIP_1pt(J,K)=0.0
420	DTAUIP_1pt(J,K)=TLIMIT
421	ENDIF
422
423
424	c IF (IPRINT .GT. 9)
425	c & WRITE(6,73)J,K,TAUGAS,TAEROS,QEXT,QSCT
426	73 FORMAT(2I3,1P8E10.3)
427
428
429	c------------------------------------------------------------------------
430	101 CONTINUE ! FIN BOUCLE L D'O
431	c------------------------------------------------------------------------
432
433	iopti=1
434
435	c************************************************************************
436	100 CONTINUE ! FIN BOUCLE ALTITUDE
437	c************************************************************************
438
439	DO 119 K=1,NSPECI
440	TAUI_1pt(1,K)=0.0
441	TAUIP_1pt(1,K)=0.0
442	DO 119 J=1,NLAYER
443	TAUI_1pt(J+1,K)=TAUI_1pt(J,K)+DTAUI_1pt(J,K)
444	TAUIP_1pt(J+1,K)=TAUIP_1pt(J,K)+DTAUIP_1pt(J,K)
445	119 CONTINUE
446
447	c IF (IPRINT .GT. 2) THEN
448	c WRITE (6,120)
449	c 120 FORMAT(///' OPTICAL CONSTANTS IN THE INFRARED')
450
451	c DO 200 K=1,NSPECI ! #2
452	c WRITE (6,190)
453	c WRITE (6,210)K,WLNI(K),WNOI(K),BWNI(K)
454	c & ,BWNI(K)+DWNI(K),DWNI(K)
455	c WRITE (6,230)REALI(K),XIMGI(K)
456
457	c DO 195 J=1,NLAYER ! #3
458	c WRITE (6,220)XNUMB(J), WBARI_1pt(J,K),COSBI_1pt(J,K)
459	c & , DTAUI_1pt(J,K),TAUI_1pt(J,K)
460	c 195 CONTINUE
461
462	c 200 CONTINUE
463
464	c END IF
465
466
467	c 210 FORMAT(1X,I3,F10.3,F10.2,F10.2,'-',F8.2,F10.3)
468	c 190 FORMAT(1X//' SNUM MICRONS WAVENU INTERVAL DELTA-WN')
469	c 230 FORMAT(1X,'NREAL(LAYER)= ',1PE10.3,' NIMG(LAYER)= ',E10.3/
470	c &' #AEROSOLS WBAR COSBAR DTAU TAU')
471	c 220 FORMAT(5(1X,G9.3))
472	c 240 FORMAT(41X,G9.3)
473
474	RETURN
475	END

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

Download in other formats:

Original Format