Context Navigation

optci_1pt.F @ 1242

Last change on this file since 1242 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: 14.0 KB

Rev	Line
[175]	1	SUBROUTINE optci_1pt(zqaer_1pt,rcdb,xfrb,iopti,IPRINT)
[102]	2	use dimphy
[3]	3	#include "dimensions.h"
	4	#include "microtab.h"
	5	#include "numchimrad.h"
	6	#include "clesphys.h"
	7
	8	PARAMETER(NLAYER=llm,NLEVEL=NLAYER+1)
	9	PARAMETER (NSPECI=46,NSPC1I=47,NSPECV=24,NSPC1V=25)
	10
	11	c Arguments:
	12	c ---------
	13	integer IPRINT,iopti
	14	C iopti: premier appel, on ne calcule qu'une fois les QM et QF
	15	* nrad dans microtab.h
[175]	16	real zqaer_1pt(NLAYER,2*nrad)
	17	#include "optci_1pt.h"
[3]	18	c ---------
	19
	20	COMMON /ATM/ Z(NLEVEL),PRESS(NLEVEL),DEN(NLEVEL),TEMP(NLEVEL)
	21
	22	COMMON /GASS/ CH4(NLEVEL),XN2(NLEVEL),H2(NLEVEL),AR(NLEVEL)
	23	& ,XMU(NLEVEL),GAS1(NLAYER),COLDEN(NLAYER)
	24
	25	COMMON /STRATO/ C2H2(NLAYER),C2H6(NLAYER)
	26	COMMON /STRAT2/ HCN(NLAYER)
	27
	28	COMMON /AERSOL/ RADIUS(NLAYER), XNUMB(NLAYER)
	29	& , REALI(NSPECI), XIMGI(NSPECI), REALV(NSPECV), XIMGV(NSPECV)
	30
[175]	31	COMMON /CLOUD/
	32	& RCLDI(NSPECI), XICLDI(NSPECI)
	33	& , RCLDV(NSPECV), XICLDV(NSPECV)
	34	& , RCLDI2(NSPECI), XICLDI2(NSPECI)
	35	& , RCLDV2(NSPECV), XICLDV2(NSPECV)
[3]	36
	37	COMMON /SPECTI/ BWNI(NSPC1I), WNOI(NSPECI),
	38	& DWNI(NSPECI), WLNI(NSPECI)
	39
[495]	40	COMMON /PLANT/ CSUBP,F0PI
[3]	41	COMMON /ADJUST/ RHCH4,FH2,FHAZE,FHVIS,FHIR,TAUFAC,RCLOUD,FARGON
	42	COMMON /CONST/RGAS,RHOP,PI,SIGMA
	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)
[175]	54	REAL QC1(nrad,NSPECI),QC2(nrad,NSPECI)
	55	REAL QC3(nrad,NSPECI),QC4(nrad,NSPECI)
[3]	56	real emu
	57	REAL TAEROSM1(NSPECI),TAEROSCATM1(NSPECI),DELTAZM1(NSPECI)
	58
[175]	59	c ---- nuages
	60	REAL TNUAGE,TNUAGESCAT
	61	REAL rcdb(nlayer),xfrb(nlayer,4)
	62
[3]	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))
[808]	92	c attention ici, Z en km doit etre passe en m
[3]	93	COEF1=RGAS273.152.5E5* (PRESS(J+1)2 - PRESS(J)2)
[808]	94	& /(1.01325*2 EFFG(Z(J)1000.)TBAR*BMU)
[3]	95
[808]	96	IF (IPRINT .GT. 9) WRITE(6,21) J,EFFG(Z(J)*1000.),TBAR,BMU,COEF1
[3]	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
[175]	196	c CALL OPTFRAC(XMONO,10000./WNOI(K)
	197	c & ,QEXT,QSCT,QABS,QBAR)
[3]	198
	199
[175]	200	CALL CFFFV11(1.e-2/WNOI(K),REALI(K),XIMGI(K),RF(inq),2.
	201	& ,XMONO,QSCT,QEXT,QABS,QBAR)
[3]	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
[175]	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.
[3]	339	ENDIF
	340
[175]	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
[3]	353
[175]	354	IF(TNUAGESCAT.EQ.0.) THEN
	355	CNBAR=0.
	356	ELSE
	357	CNBAR=CNBAR/TNUAGESCAT
	358	ENDIF
	359	ENDIF ! Cond CLD
	360
	361
[3]	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
[130]	391	c if(ig.eq.1.and.k.eq.nspecv/2) print*,'Avant GAS2'
[3]	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
[130]	397	c if(ig.eq.1.and.k.eq.nspecv/2) print*,'Apres GAS2'
[3]	398	TAUGAS=TAUGAS+TAU2
	399	ENDIF
	400	C
	401
[175]	402	DTAUI_1pt(J,K)=TAUGAS+TAEROS+TNUAGE
	403	DTAUIP_1pt(J,K)=TAUGAS+TAEROS
[3]	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
[175]	411	TAUCI_1pt(K)=TAUCI_1pt(K) + TNUAGE
	412	TAUCID_1pt(J,K)=TAUCI_1pt(K)
[3]	413
	414	C ??FLAG? SERIOUS PROBLEM WITH THE CODE HERE!
	415
	416	TLIMIT=1.E-16
	417
	418
[175]	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
[3]	425
[175]	426	IF (TAEROSCAT .GT. 0.) THEN
	427	COSBIP_1pt(J,K)=(CBAR*TAEROSCAT)
	428	& /(TAEROSCAT)
[3]	429	ELSE
[175]	430	COSBIP_1pt(J,K)=0.0
	431	ENDIF
[3]	432
[175]	433	*---------
	434
	435	IF (DTAUI_1pt(J,K) .GT. TLIMIT) THEN
	436	WBARI_1pt(J,K)=(TAEROSCAT+TNUAGESCAT) /DTAUI_1pt(J,K)
	437	ELSE
[3]	438	WBARI_1pt(J,K)=0.0
	439	DTAUI_1pt(J,K)=TLIMIT
[175]	440	ENDIF
[3]	441
[175]	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
[3]	447	ENDIF
	448
[175]	449
[3]	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
[175]	467	TAUIP_1pt(1,K)=0.0
[3]	468	DO 119 J=1,NLAYER
	469	TAUI_1pt(J+1,K)=TAUI_1pt(J,K)+DTAUI_1pt(J,K)
[175]	470	TAUIP_1pt(J+1,K)=TAUIP_1pt(J,K)+DTAUIP_1pt(J,K)
[3]	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