Context Navigation

optci_1pt_2.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: 13.7 KB

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

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

Download in other formats:

Original Format