Context Navigation

optci_1pt_3.F @ 1242

Last change on this file since 1242 was 1126, checked in by slebonnois, 11 years ago
SL: update of Titan photochemical module to include computation of chemistry up to 1300 km
File size: 13.3 KB

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

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

Download in other formats:

Original Format