Context Navigation

optci_1pt_2.F @ 1243

Last change on this file since 1243 was 808, checked in by slebonnois, 12 years ago
SL: Many changes for VENUS (related to newstart) and TITAN (related to clouds). Please read DOC/chantiers/commit_importants.log (cf v808).
File size: 13.7 KB

Rev	Line
[175]	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
[495]	41	COMMON /PLANT/ CSUBP,F0PI
[175]	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))
[808]	92	c attention ici, Z en km doit etre passe en m
[175]	93	COEF1=RGAS273.152.5E5* (PRESS(J+1)2 - PRESS(J)2)
[808]	94	& /(1.01325*2 EFFG(Z(J)1000.)TBAR*BMU)
[175]	95
[808]	96	IF (IPRINT .GT. 9) WRITE(6,21) J,EFFG(Z(J)*1000.),TBAR,BMU,COEF1
[175]	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