Context Navigation

optci_1pt_3.F @ 3557

Last change on this file since 3557 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.1 KB

Line
1	SUBROUTINE optci_1pt3(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	IMPLICIT NONE
7	#include "dimensions.h"
8	#include "microtab.h"
9	#include "numchimrad.h"
10	#include "clesphys.h"
11
12	integer nlayer, nlevel, nspeci, nspc1i, nspecv, nspc1v
13	real z,press, den, temp, ch4, xn2, h2, ar, xmu, gas1,
14	& colden, c2h2, c2h6, hcn, radius, xnumb,
15	& ximgi, realv, ximgv, rcldi, xicldi, rcldv, xicldv, rcldi2,
16	& xicldi2, rcldv2, xicldv2
17	real bwni, wnoi, dwni, wlni,
18	& prod,reali
19
20	integer k, j,inq,kgas
21
22	real tbar, pbar, bmu, coef1, effg, taeros, taeroscat, cbar,
23	& qext, qsct, qabs, qbar, xmono, xrule, deltaz, tnuext,
24	& tnuscat, cnbar, qextc, qsctc, qabsc, qbarc, taugas, pnn,
25	& pcc, pcn, phn, u, ig, tau2, tlimit
26
27	PARAMETER(NLAYER=llm,NLEVEL=NLAYER+1)
28	PARAMETER (NSPECI=46,NSPC1I=47,NSPECV=24,NSPC1V=25)
29
30	c Arguments:
31	c ---------
32	integer IPRINT,iopti
33	C iopti: premier appel, on ne calcule qu'une fois les QM et QF
34	* nrad dans microtab.h
35	real zqaer_1pt(NLAYER,2*nrad)
36	#include "optci_1pt.h"
37	c ---------
38
39	COMMON /ATM/ Z(NLEVEL),PRESS(NLEVEL),DEN(NLEVEL),TEMP(NLEVEL)
40
41	COMMON /GASS/ CH4(NLEVEL),XN2(NLEVEL),H2(NLEVEL),AR(NLEVEL)
42	& ,XMU(NLEVEL),GAS1(NLAYER),COLDEN(NLAYER)
43
44	COMMON /STRATO/ C2H2(NLAYER),C2H6(NLAYER)
45	COMMON /STRAT2/ HCN(NLAYER)
46
47	COMMON /AERSOL/ RADIUS(NLAYER), XNUMB(NLAYER)
48	& , REALI(NSPECI), XIMGI(NSPECI), REALV(NSPECV), XIMGV(NSPECV)
49
50	COMMON /CLOUD/
51	& RCLDI(NSPECI), XICLDI(NSPECI)
52	& , RCLDV(NSPECV), XICLDV(NSPECV)
53	& , RCLDI2(NSPECI), XICLDI2(NSPECI)
54	& , RCLDV2(NSPECV), XICLDV2(NSPECV)
55
56	COMMON /SPECTI/ BWNI(NSPC1I), WNOI(NSPECI),
57	& DWNI(NSPECI), WLNI(NSPECI)
58
59	COMMON /part/v,rayon,vrat,dr,dv
60
61	DIMENSION PROD(NLEVEL)
62	* nrad dans microtab.h
63	real v(nrad),rayon(nrad),vrat,dr(nrad),dv(nrad)
64	real xv1(klev,nspeci),xv2(klev,nspeci)
65	real xv3(klev,nspeci)
66	REAL QF1(nrad,NSPECI),QF2(nrad,NSPECI)
67	REAL QF3(nrad,NSPECI),QF4(nrad,NSPECI)
68	REAL QM1(nrad,NSPECI),QM2(nrad,NSPECI)
69	REAL QM3(nrad,NSPECI),QM4(nrad,NSPECI)
70	REAL QC1(nrad,NSPECI),QC2(nrad,NSPECI)
71	REAL QC3(nrad,NSPECI),QC4(nrad,NSPECI)
72	real emu
73	REAL TAEROSM1(NSPECI),TAEROSCATM1(NSPECI),DELTAZM1(NSPECI)
74
75	c ---- nuages
76	REAL TNUAGE,TNUAGESCAT
77	REAL rcdb(nlayer),xfrb(nlayer,4)
78
79	save qf1,qf2,qf3,qf4,qm1,qm2,qm3,qm4,qc1,qc2,qc3,qc4
80
81
82	C THE PRESSURE INDUCED TRANSITIONS ARE FROM REGIS
83	C THE LAST SEVENTEEN INTERVALS ARE THE BANDS FROM GNF.
84	C
85	C THIS SUBROUTINE SETS THE OPTICAL CONSTANTS IN THE INFRARED
86	C IT CALCUALTES FOR EACH LAYER, FOR EACH SPECRAL INTERVAL IN THE IR
87	C LAYER: WBAR, DTAU, COSBAR
88	C LEVEL: TAU
89	C
90
91	DO 80 K=1,NSPECI
92	TAUHI_1pt(K)=0.
93	TAUCI_1pt(K)=0.
94	TAUGI_1pt(K)=0.
95	80 CONTINUE
96
97	c************************************************************************
98	DO 100 J=1,NLAYER ! BOUCLE SUR L'ALTITUDE
99	c************************************************************************
100
101	C SET UP THE COEFFICIENT TO REDUCE MASS PATH TO STP ...SEE NOTES
102	C T0 =273.15 PO=1.01325 BAR
103
104	TBAR=0.5*(TEMP(J)+TEMP(J+1))
105	PBAR=SQRT(PRESS(J)*PRESS(J+1))
106	BMU=0.5*(XMU(J+1)+XMU(J))
107	c attention ici, Z en km doit etre passe en m
108	COEF1=RGAS273.152.5E5* (PRESS(J+1)2 - PRESS(J)2)
109	& /(1.01325*2 EFFG(Z(J)1000.)TBAR*BMU)
110
111	IF (IPRINT .GT. 9) WRITE(6,21) J,EFFG(Z(J)*1000.),TBAR,BMU,COEF1
112	21 FORMAT(' J, EFFG, TBAR, BMU, COEF1,: ',I3,1P6E10.3)
113
114	c------------------------------------------------------------------------
115	DO 101 K=1,NSPECI ! BOUCLE SUR LES L.D'O
116	c------------------------------------------------------------------------
117
118
119	C #1: HAZE
120	C---------------------
121
122	C FIRST COMPUTE TAU AEROSOL
123
124
125	c
126	c /\
127	c / \
128	c / \
129	c / _O \
130	c / \|/ \
131	c / / \ \
132	c / \|\ \/\ \
133	c / \|\| / \ \
134	c ----------------
135	c \| WARNING \|
136	c \| SLOW DOWN \|
137	c ----------------
138
139
140
141
142	c********* EN TRAVAUX *************************
143
144	TAEROS=0.
145	TAEROSCAT=0.
146	CBAR=0.
147
148
149	DO inq=1,nrad !BOUCLE SUR LES TAILLE D"AEROSOLS
150
151
152	IF (WNOI(K).lt.wco) THEN ! lamda > 56 um
153
154	if (iopti.eq.0) then
155
156	c CALL XMIE(rayon(inq)*1.e6,REALI(K),XIMGI(K),
157	c & QEXT,QSCT,QABS,QBAR,WNOI(K))
158
159	CALL CMIE(1.E-2/WNOI(K),REALI(K),XIMGI(K),rayon(inq),
160	& QEXT,QSCT,QABS,QBAR)
161
162
163	QM1(inq,K)=QEXT
164	QM2(inq,K)=QSCT
165	QM3(inq,K)=QABS
166	QM4(inq,K)=QBAR
167
168	endif ! end iopti
169
170
171	TAEROS=QM1(inq,K)zqaer_1pt(nlayer+1-J,inq)1.e-4+TAEROS
172	TAEROSCAT=QM2(inq,K)zqaer_1pt(nlayer+1-J,inq)1.e-4+TAEROSCAT
173	CBAR=CBAR+QM4(inq,K)QM2(inq,K)zqaer_1pt(nlayer+1-J,inq)*1.e-4
174
175
176	ELSE ! 0.2 < lambda < 56 um
177
178
179	if(rayon(inq).lt.RF(inq)) THEN
180
181	if (iopti.eq.0) then
182
183	CALL XMIE(rayon(inq)*1.e6,REALI(K),XIMGI(K),
184	& QEXT,QSCT,QABS,QBAR,WNOI(K))
185
186	QM1(inq,K)=QEXT
187	QM2(inq,K)=QSCT
188	QM3(inq,K)=QABS
189	QM4(inq,K)=QBAR
190	endif ! end iopti
191
192
193	TAEROS=QM1(inq,K)zqaer_1pt(nlayer+1-J,inq)1.e-4+TAEROS
194	TAEROSCAT=QM2(inq,K)zqaer_1pt(nlayer+1-J,inq)1.e-4+TAEROSCAT
195	CBAR=CBAR+QM4(inq,K)QM2(inq,K)zqaer_1pt(nlayer+1-J,inq)*1.e-4
196
197	else
198
199	XMONO=(rayon(inq)/RF(inq))**3.
200	XRULE=1.
201
202	if(XMONO.gt.16384./1.5) then
203	XRULE=(XMONO/16384.)
204	XMONO=16384.
205	endif
206
207
208	if (iopti.eq.0) then
209
210	c CALL OPTFRAC(XMONO,10000./WNOI(K)
211	c & ,QEXT,QSCT,QABS,QBAR)
212
213
214	CALL CFFFV11(1.e-2/WNOI(K),REALI(K),XIMGI(K),RF(inq),2.
215	& ,XMONO,QSCT,QEXT,QABS,QBAR)
216
217
218	QF1(inq,K)=QEXT*XRULE
219	QF2(inq,K)=QSCT*XRULE
220	QF3(inq,K)=QABS*XRULE
221	QF4(inq,K)=QBAR
222	endif ! end iopti
223
224	TAEROS=QF1(inq,K)*zqaer_1pt(nlayer+1-J,inq)+TAEROS
225	TAEROSCAT=QF2(inq,K)*zqaer_1pt(nlayer+1-J,inq)+TAEROSCAT
226	CBAR=CBAR+QF4(inq,K)QF2(inq,K)zqaer_1pt(nlayer+1-J,inq)
227
228	endif
229
230	IF(TAEROS.LT.1.e-10) TAEROS=1.e-10
231
232	ENDIF
233	ENDDO ! FIN DE LA BOUCLE SUR nrad
234
235
236
237
238	if (TAEROSCAT.ne.0.) CBAR=CBAR/TAEROSCAT
239
240	DELTAZ=Z(J)-Z(J+1)
241
242	c --------------------------------------------------------------------
243	c profil brume Pascal: fit T (sauf tropopause) et albedo
244	c -------------------
245	if( cutoff.eq.1) then
246	IF(PRESS(J).gt.9.e-3) THEN
247	TAEROS=TAEROSM1(K)DELTAZ/DELTAZM1(K)0.85
248	TAEROSCAT=TAEROSCATM1(K)DELTAZ/DELTAZM1(K)0.85
249	c TAEROS=0.
250	c TAEROSCAT=0.
251	ENDIF
252
253	IF(PRESS(J).gt.1.e-1) THEN
254	TAEROS=TAEROSM1(K)DELTAZ/DELTAZM1(K)1.15
255	TAEROSCAT=TAEROSCATM1(K)DELTAZ/DELTAZM1(K)1.15
256	c TAEROS=0.
257	c TAEROSCAT=0.
258	ENDIF
259	endif !cutoff=1
260
261	c profil brume pour fit T (y compris tropopause), mais ne fit plus albedo...
262	c -----------------------
263	if( cutoff.eq.2) then
264	IF(PRESS(J).gt.1.e-1) THEN
265	TAEROS=0.
266	TAEROSCAT=0.
267	ENDIF
268	endif !cutoff=2
269	c --------------------------------------------------------------------
270
271	TAEROSM1(K)=TAEROS
272	TAEROSCATM1(K)=TAEROSCAT
273	DELTAZM1(K)=DELTAZ
274
275
276	IF(TAEROSCAT.le.0.) CBAR=0.
277
278	c print*,'HERE, MCKAY AEROSOLS IR'
279	c TAEROS=xv1(j,k)
280	c TAEROSCAT=xv2(j,k)
281	c CBAR=xv3(j,k)
282
283	c print*, 'HERE, CIRS AEROSOLS'
284	c call cirs_haze(PRESS(J),WNOI(K),TAEROS,TAEROSCAT,CBAR)
285
286	c********* EN TRAVAUX *************************
287
288	C #2: CLOUD
289	C------------------
290	C NEXT COMPUTE TAU CLOUD
291	c
292	c Menu special :
293	c On utilise ici une look-up table afin de calculer
294	c les proprietes optique des nuages.
295	c Le principe est le suivant :
296	c La look-up table contient les proprietes optique d'une goutte
297	c de methane pur de 3 um.
298	c On approxime les proprietes optiques pour une goutte de rayon r a
299	c de la table.
300	c
301	c
302	TNUEXT=0.
303	TNUSCAT=0.
304	CNBAR=0.
305	IF (clouds.eq.1) THEN
306
307	CALL getoptcld(1.E-2/WNOI(K),rcdb(nlayer+1-J),
308	& QEXTC,QSCTC,QABSC,QBARC)
309
310
311	c ----- On ne calcule les constante optiques que si Rgoutte > 1e-10
312	IF (rcdb(nlayer+1-J).gt.1.1e-10) THEN
313	TNUEXT =QEXTC/xnufSUM(zqaer_1pt(NLAYER+1-J,nrad+1:2nrad))
314	TNUSCAT=QSCTC/xnufSUM(zqaer_1pt(NLAYER+1-J,nrad+1:2nrad))
315	CNBAR =QBARC
316	ENDIF
317
318	IF(TNUSCAT.EQ.0.) THEN
319	CNBAR=0.
320	ELSE
321	CNBAR=CNBAR/TNUSCAT
322	ENDIF
323
324	ENDIF ! Cond CLD
325	c
326	C #3: GAZ
327	C------------------
328
329	C NOW COMPUTE TAUGAS DUE TO THE PIA TERM ONLY FOR LAMDA LT 940
330	TAUGAS=0.0
331	IF (WNOI(K) .LT. 940. ) THEN
332	CALL PIA(K,TBAR,PNN,PCC,PCN,PHN)
333	C HERE IS WHERE WE COULD SCALE THE PIA COEFFICEINTS TO FIT DATA
334	C BASED ON REGIS' NOTES. ---TGM HAS THIS ADJUST IN IT AS DEFAULT
335	PCN=PCN*MIN(1.75 , AMAX1(1.0,WNOI(K)/200.))
336	C**REPLACE ABOVE WITH: PCN=PCN1.25*MIN(1.75 , AMAX1(1.0,WNOI(K)/200.))
337	C 1.25 FACTOR (NOT FROM DATA) SUGGESTED BY TOON et al. (1988)
338	TAUGAS=COEF1*
339	& (XN2(J)XN2(J)PNN + CH4(J)CH4(J)PCC
340	& + XN2(J)CH4(J)PCN + XN2(J)H2(J)PHN)
341	IF (J .EQ. NLAYER .AND. IPRINT .GT. 9)
342	& WRITE (6,22) WNOI(K),TAUGAS,XN2(J),CH4(J),H2(J),
343	& TBAR, PNN,PCC,PCN, PHN,
344	& XN2(J)XN2(J)PNN , CH4(J)CH4(J)PCC ,
345	& XN2(J)CH4(J)PCN , XN2(J)H2(J)PHN
346	22 FORMAT(1X,1P8E10.2)
347	ENDIF
348
349	IF (K .GT. 28) THEN
350	KGAS=K-28
351	C ??FLAG? HERE MUST BE WATCHED CAREFULLY
352	U=COLDEN(J)*6.02204E23/BMU
353	if((ylellouch).or.(.not.hcnrad)) then
354	CALL GAS2_NOHCN(J, KGAS,TBAR,PBAR,U,TAU2)
355	else
356	CALL GAS2(J, KGAS,TBAR,PBAR,U,TAU2)
357	endif
358	TAUGAS=TAUGAS+TAU2
359	ENDIF
360	C
361
362	DTAUI_1pt(J,K)=TAUGAS+TAEROS+TNUEXT
363	DTAUIP_1pt(J,K)=TAUGAS+TAEROS
364
365	TAUHI_1pt(K)=TAUHI_1pt(K) + TAEROS
366	TAUHID_1pt(J,K)=TAUHI_1pt(K)
367
368	TAUGI_1pt(K)=TAUGI_1pt(K) + TAUGAS
369	TAUGID_1pt(J,K)=TAUGI_1pt(K)
370
371	TAUCI_1pt(K)=TAUCI_1pt(K) + TNUEXT
372	TAUCID_1pt(J,K)=TAUCI_1pt(K)
373
374	C ??FLAG? SERIOUS PROBLEM WITH THE CODE HERE!
375
376	TLIMIT=1.E-16
377
378
379	IF (TAEROSCAT + TNUSCAT .GT. 0.) THEN
380	COSBI_1pt(J,K)=(CBARTAEROSCAT + CNBARTNUSCAT )
381	& /(TAEROSCAT + TNUSCAT)
382	ELSE
383	COSBI_1pt(J,K)=0.0
384	ENDIF
385
386	IF (TAEROSCAT .GT. 0.) THEN
387	COSBIP_1pt(J,K)=(CBAR*TAEROSCAT)
388	& /(TAEROSCAT)
389	ELSE
390	COSBIP_1pt(J,K)=0.0
391	ENDIF
392
393	*---------
394
395	IF (DTAUI_1pt(J,K) .GT. TLIMIT) THEN
396	WBARI_1pt(J,K)=(TAEROSCAT+TNUSCAT) /DTAUI_1pt(J,K)
397	ELSE
398	WBARI_1pt(J,K)=0.0
399	DTAUI_1pt(J,K)=TLIMIT
400	ENDIF
401
402	IF (DTAUIP_1pt(J,K) .GT. TLIMIT) THEN
403	WBARIP_1pt(J,K)=(TAEROSCAT) /DTAUIP_1pt(J,K)
404	ELSE
405	WBARIP_1pt(J,K)=0.0
406	DTAUIP_1pt(J,K)=TLIMIT
407	ENDIF
408
409
410	c IF (IPRINT .GT. 9)
411	c & WRITE(6,73)J,K,TAUGAS,TAEROS,QEXT,QSCT
412	73 FORMAT(2I3,1P8E10.3)
413
414
415	c------------------------------------------------------------------------
416	101 CONTINUE ! FIN BOUCLE L D'O
417	c------------------------------------------------------------------------
418
419	iopti=1
420
421	c************************************************************************
422	100 CONTINUE ! FIN BOUCLE ALTITUDE
423	c************************************************************************
424
425	DO 119 K=1,NSPECI
426	TAUI_1pt(1,K)=0.0
427	TAUIP_1pt(1,K)=0.0
428	DO 119 J=1,NLAYER
429	TAUI_1pt(J+1,K)=TAUI_1pt(J,K)+DTAUI_1pt(J,K)
430	TAUIP_1pt(J+1,K)=TAUIP_1pt(J,K)+DTAUIP_1pt(J,K)
431	119 CONTINUE
432
433	c IF (IPRINT .GT. 2) THEN
434	c WRITE (6,120)
435	c 120 FORMAT(///' OPTICAL CONSTANTS IN THE INFRARED')
436
437	c DO 200 K=1,NSPECI ! #2
438	c WRITE (6,190)
439	c WRITE (6,210)K,WLNI(K),WNOI(K),BWNI(K)
440	c & ,BWNI(K)+DWNI(K),DWNI(K)
441	c WRITE (6,230)REALI(K),XIMGI(K)
442
443	c DO 195 J=1,NLAYER ! #3
444	c WRITE (6,220)XNUMB(J), WBARI_1pt(J,K),COSBI_1pt(J,K)
445	c & , DTAUI_1pt(J,K),TAUI_1pt(J,K)
446	c 195 CONTINUE
447
448	c 200 CONTINUE
449
450	c END IF
451
452
453	c 210 FORMAT(1X,I3,F10.3,F10.2,F10.2,'-',F8.2,F10.3)
454	c 190 FORMAT(1X//' SNUM MICRONS WAVENU INTERVAL DELTA-WN')
455	c 230 FORMAT(1X,'NREAL(LAYER)= ',1PE10.3,' NIMG(LAYER)= ',E10.3/
456	c &' #AEROSOLS WBAR COSBAR DTAU TAU')
457	c 220 FORMAT(5(1X,G9.3))
458	c 240 FORMAT(41X,G9.3)
459
460	RETURN
461	END

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

Download in other formats:

Original Format