Context Navigation

optci_1pt.F @ 3094

Last change on this file since 3094 was 1461, checked in by emillour, 10 years ago
Titan GCM: Turned the common block "tgmdat.F" into a module "tgmdat_mod.F90". This fixes issues in "debug" mode with common variables which seemed to not be correctly shared between routines. EM
File size: 14.0 KB

Line
1	SUBROUTINE optci_1pt(zqaer_1pt,rcdb,xfrb,iopti,IPRINT)
2	use dimphy
3	USE TGMDAT_MOD, ONLY: RHCH4,FH2,FHAZE,FHVIS,FHIR,TAUFAC,
4	& RCLOUD,FARGON
5	USE TGMDAT_MOD, ONLY: RGAS
6	#include "dimensions.h"
7	#include "microtab.h"
8	#include "numchimrad.h"
9	#include "clesphys.h"
10
11	PARAMETER(NLAYER=llm,NLEVEL=NLAYER+1)
12	PARAMETER (NSPECI=46,NSPC1I=47,NSPECV=24,NSPC1V=25)
13
14	c Arguments:
15	c ---------
16	integer IPRINT,iopti
17	C iopti: premier appel, on ne calcule qu'une fois les QM et QF
18	* nrad dans microtab.h
19	real zqaer_1pt(NLAYER,2*nrad)
20	#include "optci_1pt.h"
21	c ---------
22
23	COMMON /ATM/ Z(NLEVEL),PRESS(NLEVEL),DEN(NLEVEL),TEMP(NLEVEL)
24
25	COMMON /GASS/ CH4(NLEVEL),XN2(NLEVEL),H2(NLEVEL),AR(NLEVEL)
26	& ,XMU(NLEVEL),GAS1(NLAYER),COLDEN(NLAYER)
27
28	COMMON /STRATO/ C2H2(NLAYER),C2H6(NLAYER)
29	COMMON /STRAT2/ HCN(NLAYER)
30
31	COMMON /AERSOL/ RADIUS(NLAYER), XNUMB(NLAYER)
32	& , REALI(NSPECI), XIMGI(NSPECI), REALV(NSPECV), XIMGV(NSPECV)
33
34	COMMON /CLOUD/
35	& RCLDI(NSPECI), XICLDI(NSPECI)
36	& , RCLDV(NSPECV), XICLDV(NSPECV)
37	& , RCLDI2(NSPECI), XICLDI2(NSPECI)
38	& , RCLDV2(NSPECV), XICLDV2(NSPECV)
39
40	COMMON /SPECTI/ BWNI(NSPC1I), WNOI(NSPECI),
41	& DWNI(NSPECI), WLNI(NSPECI)
42
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	c attention ici, Z en km doit etre passe en m
93	COEF1=RGAS273.152.5E5* (PRESS(J+1)2 - PRESS(J)2)
94	& /(1.01325*2 EFFG(Z(J)1000.)TBAR*BMU)
95
96	IF (IPRINT .GT. 9) WRITE(6,21) J,EFFG(Z(J)*1000.),TBAR,BMU,COEF1
97	21 FORMAT(' J, EFFG, TBAR, BMU, COEF1,: ',I3,1P6E10.3)
98
99	c------------------------------------------------------------------------
100	DO 101 K=1,NSPECI ! BOUCLE SUR LES L.D'O
101	c------------------------------------------------------------------------
102
103
104	C #1: HAZE
105	C---------------------
106
107	C FIRST COMPUTE TAU AEROSOL
108
109
110	c
111	c /\
112	c / \
113	c / \
114	c / _O \
115	c / \|/ \
116	c / / \ \
117	c / \|\ \/\ \
118	c / \|\| / \ \
119	c ----------------
120	c \| WARNING \|
121	c \| SLOW DOWN \|
122	c ----------------
123
124
125
126
127	c********* EN TRAVAUX *************************
128
129	TAEROS=0.
130	TAEROSCAT=0.
131	CBAR=0.
132
133
134	DO inq=1,nrad !BOUCLE SUR LES TAILLE D"AEROSOLS
135
136
137	IF (WNOI(K).lt.wco) THEN ! lamda > 56 um
138
139	if (iopti.eq.0) then
140
141	c CALL XMIE(rayon(inq)*1.e6,REALI(K),XIMGI(K),
142	c & QEXT,QSCT,QABS,QBAR,WNOI(K))
143
144
145	CALL CMIE(1.E-2/WNOI(K),REALI(K),XIMGI(K),rayon(inq),
146	& QEXT,QSCT,QABS,QBAR)
147
148
149	QM1(inq,K)=QEXT
150	QM2(inq,K)=QSCT
151	QM3(inq,K)=QABS
152	QM4(inq,K)=QBAR
153
154	endif ! end iopti
155
156
157	TAEROS=QM1(inq,K)zqaer_1pt(nlayer+1-J,inq)1.e-4+TAEROS
158	TAEROSCAT=QM2(inq,K)zqaer_1pt(nlayer+1-J,inq)1.e-4+TAEROSCAT
159	CBAR=CBAR+QM4(inq,K)QM2(inq,K)zqaer_1pt(nlayer+1-J,inq)*1.e-4
160
161
162	ELSE ! 0.2 < lambda < 56 um
163
164
165	if(rayon(inq).lt.RF(inq)) THEN
166
167	if (iopti.eq.0) then
168
169	CALL XMIE(rayon(inq)*1.e6,REALI(K),XIMGI(K),
170	& QEXT,QSCT,QABS,QBAR,WNOI(K))
171
172	QM1(inq,K)=QEXT
173	QM2(inq,K)=QSCT
174	QM3(inq,K)=QABS
175	QM4(inq,K)=QBAR
176	endif ! end iopti
177
178
179	TAEROS=QM1(inq,K)zqaer_1pt(nlayer+1-J,inq)1.e-4+TAEROS
180	TAEROSCAT=QM2(inq,K)zqaer_1pt(nlayer+1-J,inq)1.e-4+TAEROSCAT
181	CBAR=CBAR+QM4(inq,K)QM2(inq,K)zqaer_1pt(nlayer+1-J,inq)*1.e-4
182
183	else
184
185	XMONO=(rayon(inq)/RF(inq))**3.
186	XRULE=1.
187
188	if(XMONO.gt.16384./1.5) then
189	XRULE=(XMONO/16384.)
190	XMONO=16384.
191	endif
192
193
194	if (iopti.eq.0) then
195
196	c CALL OPTFRAC(XMONO,10000./WNOI(K)
197	c & ,QEXT,QSCT,QABS,QBAR)
198
199
200	CALL CFFFV11(1.e-2/WNOI(K),REALI(K),XIMGI(K),RF(inq),2.
201	& ,XMONO,QSCT,QEXT,QABS,QBAR)
202
203
204	QF1(inq,K)=QEXT*XRULE
205	QF2(inq,K)=QSCT*XRULE
206	QF3(inq,K)=QABS*XRULE
207	QF4(inq,K)=QBAR
208	endif ! end iopti
209
210	TAEROS=QF1(inq,K)*zqaer_1pt(nlayer+1-J,inq)+TAEROS
211	TAEROSCAT=QF2(inq,K)*zqaer_1pt(nlayer+1-J,inq)+TAEROSCAT
212	CBAR=CBAR+QF4(inq,K)QF2(inq,K)zqaer_1pt(nlayer+1-J,inq)
213
214	endif
215
216	IF(TAEROS.LT.1.e-10) TAEROS=1.e-10
217	ENDIF
218	ENDDO ! FIN DE LA BOUCLE SUR nrad
219
220
221
222
223	CBAR=CBAR/TAEROSCAT
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 if (ig.eq.1) then
269	c if (k.eq.NSPECV/2) then
270	c print*,'@IR',K,J,TAEROS,TAEROSCAT,CBAR
271	c stop'Pour faire des comparaisons'
272	c endif
273	c endif
274
275
276	c********* EN TRAVAUX *************************
277
278	C #2: CLOUD
279	C------------------
280
281	C NEXT COMPUTE TAU CLOUD
282	IF (clouds.eq.0) THEN
283	TNUAGE=0.
284	TNUAGESCAT=0.
285	CNBAR=0.
286	ELSE
287	TNUAGE=0.
288	TNUAGESCAT=0.
289	CNBAR=0.
290	QEXTC=0.
291	QSCTC=0.
292	QABSC=0.
293	CBARC=0.
294
295	DO inq=1,nrad !BOUCLE SUR LES NQMX TAILLE D"AEROSOLS
296	QC1(inq,K)=0.
297	QC2(inq,K)=0.
298	QC3(inq,K)=0.
299	QC4(inq,K)=0.
300	ENDDO
301
302	** OPTICAL CONSTANT : MIXING RULES
303
304	IF (rcdb(nlayer+1-J).gt.1.1e-10) THEN
305
306	XNR=xfrb(nlayer+1-J,1) *REALI(K)
307	& +xfrb(nlayer+1-J,2) *RCLDI(K)
308	& +xfrb(nlayer+1-J,3) *RCLDI2(K)
309	& +xfrb(nlayer+1-J,4) *RCLDI2(K)
310
311	XNI=xfrb(nlayer+1-J,1) *XIMGI(K)
312	& +xfrb(nlayer+1-J,2) *XICLDI(K)
313	& +xfrb(nlayer+1-J,3) *XICLDI2(K)
314	& +xfrb(nlayer+1-J,4) *XICLDI2(K)
315
316	** OPTICAL CONSTANT : LIQUID DROP = THOLIN
317
318	IF(xfrb(nlayer+1-J,1).ge.0.1) THEN
319	XNI=XIMGI(K)
320	XNR=REALI(K)
321	ENDIF
322
323	IF (XNI.gt.1.e-10 .and. XNR.gt.1.00) THEN
324	CALL CMIE(1.E-2/WNOI(K),XNR,XNI,
325	& rcdb(nlayer+1-J),
326	& QEXTC,QSCTC,QABSC,CBARC)
327	ELSE
328	PRINT*,' WARNING XNR/XNI in optci: ',XNR,XNI
329	QEXTC=0.
330	QSCTC=0.
331	QABSC=0.
332	CBARC=0.
333	ENDIF
334	ELSE
335	QEXTC=0.
336	QSCTC=0.
337	QABSC=0.
338	CBARC=0.
339	ENDIF
340
341	DO inq=1,nrad !BOUCLE SUR LES NQMX TAILLE D"AEROSOLS
342	QC1(inq,K)=QEXTC/xnuf
343	QC2(inq,K)=QSCTC/xnuf
344	QC3(inq,K)=QABSC/xnuf
345	QC4(inq,K)=CBARC
346	TNUAGE=QC1(inq,K)zqaer_1pt(nlayer+1-J,inq+nrad)1.e-4
347	& +TNUAGE
348	TNUAGESCAT=QC2(inq,K)zqaer_1pt(nlayer+1-J,inq+nrad)1.e-4
349	& +TNUAGESCAT
350	CNBAR=QC4(inq,K)*QC2(inq,K)
351	& zqaer_1pt(nlayer+1-J,inq+nrad)1.e-4+CNBAR
352	ENDDO
353
354	IF(TNUAGESCAT.EQ.0.) THEN
355	CNBAR=0.
356	ELSE
357	CNBAR=CNBAR/TNUAGESCAT
358	ENDIF
359	ENDIF ! Cond CLD
360
361
362	C #3: GAZ
363	C------------------
364
365	C NOW COMPUTE TAUGAS DUE TO THE PIA TERM ONLY FOR LAMDA LT 940
366	TAUGAS=0.0
367	IF (WNOI(K) .LT. 940. ) THEN
368	c if(ig.eq.1.and.k.eq.nspecv/2) print*,'avant PIA'
369	CALL PIA(K,TBAR,PNN,PCC,PCN,PHN)
370	c if(ig.eq.1.and.k.eq.nspecv/2) print*,'apres PIA'
371	C HERE IS WHERE WE COULD SCALE THE PIA COEFFICEINTS TO FIT DATA
372	C BASED ON REGIS' NOTES. ---TGM HAS THIS ADJUST IN IT AS DEFAULT
373	PCN=PCN*MIN(1.75 , AMAX1(1.0,WNOI(K)/200.))
374	C**REPLACE ABOVE WITH: PCN=PCN1.25*MIN(1.75 , AMAX1(1.0,WNOI(K)/200.))
375	C 1.25 FACTOR (NOT FROM DATA) SUGGESTED BY TOON et al. (1988)
376	TAUGAS=COEF1*
377	& (XN2(J)XN2(J)PNN + CH4(J)CH4(J)PCC
378	& + XN2(J)CH4(J)PCN + XN2(J)H2(J)PHN)
379	IF (J .EQ. NLAYER .AND. IPRINT .GT. 9)
380	& WRITE (6,22) WNOI(K),TAUGAS,XN2(J),CH4(J),H2(J),
381	& TBAR, PNN,PCC,PCN, PHN,
382	& XN2(J)XN2(J)PNN , CH4(J)CH4(J)PCC ,
383	& XN2(J)CH4(J)PCN , XN2(J)H2(J)PHN
384	22 FORMAT(1X,1P8E10.2)
385	ENDIF
386
387	IF (K .GT. 28) THEN
388	KGAS=K-28
389	C ??FLAG? HERE MUST BE WATCHED CAREFULLY
390	U=COLDEN(J)*6.02204E23/BMU
391	c if(ig.eq.1.and.k.eq.nspecv/2) print*,'Avant GAS2'
392	if((ylellouch).or.(.not.hcnrad)) then
393	CALL GAS2_NOHCN(J, KGAS,TBAR,PBAR,U,TAU2)
394	else
395	CALL GAS2(J, KGAS,TBAR,PBAR,U,TAU2)
396	endif
397	c if(ig.eq.1.and.k.eq.nspecv/2) print*,'Apres GAS2'
398	TAUGAS=TAUGAS+TAU2
399	ENDIF
400	C
401
402	DTAUI_1pt(J,K)=TAUGAS+TAEROS+TNUAGE
403	DTAUIP_1pt(J,K)=TAUGAS+TAEROS
404
405	TAUHI_1pt(K)=TAUHI_1pt(K) + TAEROS
406	TAUHID_1pt(J,K)=TAUHI_1pt(K)
407
408	TAUGI_1pt(K)=TAUGI_1pt(K) + TAUGAS
409	TAUGID_1pt(J,K)=TAUGI_1pt(K)
410
411	TAUCI_1pt(K)=TAUCI_1pt(K) + TNUAGE
412	TAUCID_1pt(J,K)=TAUCI_1pt(K)
413
414	C ??FLAG? SERIOUS PROBLEM WITH THE CODE HERE!
415
416	TLIMIT=1.E-16
417
418
419	IF (TAEROSCAT + TNUAGESCAT .GT. 0.) THEN
420	COSBI_1pt(J,K)=(CBARTAEROSCAT + CNBARTNUAGESCAT )
421	& /(TAEROSCAT + TNUAGESCAT)
422	ELSE
423	COSBI_1pt(J,K)=0.0
424	ENDIF
425
426	IF (TAEROSCAT .GT. 0.) THEN
427	COSBIP_1pt(J,K)=(CBAR*TAEROSCAT)
428	& /(TAEROSCAT)
429	ELSE
430	COSBIP_1pt(J,K)=0.0
431	ENDIF
432
433	*---------
434
435	IF (DTAUI_1pt(J,K) .GT. TLIMIT) THEN
436	WBARI_1pt(J,K)=(TAEROSCAT+TNUAGESCAT) /DTAUI_1pt(J,K)
437	ELSE
438	WBARI_1pt(J,K)=0.0
439	DTAUI_1pt(J,K)=TLIMIT
440	ENDIF
441
442	IF (DTAUIP_1pt(J,K) .GT. TLIMIT) THEN
443	WBARIP_1pt(J,K)=(TAEROSCAT) /DTAUIP_1pt(J,K)
444	ELSE
445	WBARIP_1pt(J,K)=0.0
446	DTAUIP_1pt(J,K)=TLIMIT
447	ENDIF
448
449
450	c IF (IPRINT .GT. 9)
451	c & WRITE(6,73)J,K,TAUGAS,TAEROS,QEXT,QSCT
452	73 FORMAT(2I3,1P8E10.3)
453
454
455	c------------------------------------------------------------------------
456	101 CONTINUE ! FIN BOUCLE L D'O
457	c------------------------------------------------------------------------
458
459	iopti=1
460
461	c************************************************************************
462	100 CONTINUE ! FIN BOUCLE ALTITUDE
463	c************************************************************************
464
465	DO 119 K=1,NSPECI
466	TAUI_1pt(1,K)=0.0
467	TAUIP_1pt(1,K)=0.0
468	DO 119 J=1,NLAYER
469	TAUI_1pt(J+1,K)=TAUI_1pt(J,K)+DTAUI_1pt(J,K)
470	TAUIP_1pt(J+1,K)=TAUIP_1pt(J,K)+DTAUIP_1pt(J,K)
471	119 CONTINUE
472
473	c IF (IPRINT .GT. 2) THEN
474	c WRITE (6,120)
475	c 120 FORMAT(///' OPTICAL CONSTANTS IN THE INFRARED')
476
477	c DO 200 K=1,NSPECI ! #2
478	c WRITE (6,190)
479	c WRITE (6,210)K,WLNI(K),WNOI(K),BWNI(K)
480	c & ,BWNI(K)+DWNI(K),DWNI(K)
481	c WRITE (6,230)REALI(K),XIMGI(K)
482
483	c DO 195 J=1,NLAYER ! #3
484	c WRITE (6,220)XNUMB(J), WBARI_1pt(J,K),COSBI_1pt(J,K)
485	c & , DTAUI_1pt(J,K),TAUI_1pt(J,K)
486	c 195 CONTINUE
487
488	c IF(ig.eq.12) WRITE (6,240) TAUI_1pt(NLEVEL,K)
489	c 200 CONTINUE
490
491	c END IF
492
493
494	c 210 FORMAT(1X,I3,F10.3,F10.2,F10.2,'-',F8.2,F10.3)
495	c 190 FORMAT(1X//' SNUM MICRONS WAVENU INTERVAL DELTA-WN')
496	c 230 FORMAT(1X,'NREAL(LAYER)= ',1PE10.3,' NIMG(LAYER)= ',E10.3/
497	c &' #AEROSOLS WBAR COSBAR DTAU TAU')
498	c 220 FORMAT(5(1X,G9.3))
499	c 240 FORMAT(41X,G9.3)
500
501	RETURN
502	END

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

Download in other formats:

Original Format