Context Navigation

optcv_1pt.F @ 855

Last change on this file since 855 was 495, checked in by slebonnois, 14 years ago
Mise a jour physique Titan, ajout des forces de marees (dans la dynamique, sous flag titan). SL.
File size: 13.5 KB

Rev	Line
[175]	1	SUBROUTINE optcv_1pt(zqaer_1pt,rcdb,xfrb,ioptv,IPRINT)
[3]	2
	3
[102]	4	use dimphy
[3]	5	#include "dimensions.h"
	6	#include "microtab.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,ioptv
	15	C ioptv: premier appel, on ne calcule qu'une fois les QM et QF
	16	* nrad dans microtab.h
[175]	17	real zqaer_1pt(NLAYER,2*nrad)
	18	#include "optcv_1pt.h"
[3]	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 /VISGAS/SOLARF(NSPECV),NTERM(NSPECV),PEXPON(NSPECV),
	27	& ATERM(4,NSPECV),BTERM(4,NSPECV)
	28
	29	COMMON /AERSOL/ RADIUS(NLAYER), XNUMB(NLAYER)
	30	& , REALI(NSPECI), XIMGI(NSPECI), REALV(NSPECV), XIMGV(NSPECV)
	31
[175]	32	COMMON /CLOUD/
	33	& RCLDI(NSPECI), XICLDI(NSPECI)
[3]	34	& , RCLDV(NSPECV), XICLDV(NSPECV)
[175]	35	& , RCLDI2(NSPECI), XICLDI2(NSPECI)
	36	& , RCLDV2(NSPECV), XICLDV2(NSPECV)
[3]	37
	38	COMMON /SPECTV/ BWNV(NSPC1V),WNOV(NSPECV)
	39	& ,DWNV(NSPECV),WLNV(NSPECV)
	40
[495]	41	COMMON /PLANT/ CSUBP,F0PI
[3]	42	COMMON /ADJUST/ RHCH4,FH2,FHAZE,FHVIS,FHIR,TAUFAC,RCLOUD,FARGON
	43	COMMON /CONST/ RGAS,RHOP,PI,SIGMA
	44	* nrad dans microtab.h
	45	COMMON /part/ v(nrad),rayon(nrad),vrat,dr(nrad),dv(nrad)
	46
	47	REAL QF1(nrad,NSPECV),QF2(nrad,NSPECV)
	48	REAL QF3(nrad,NSPECV),QF4(nrad,NSPECV)
	49	REAL QM1(nrad,NSPECV),QM2(nrad,NSPECV)
	50	REAL QM3(nrad,NSPECV),QM4(nrad,NSPECV)
[175]	51	REAL QC1(nrad,NSPECV),QC2(nrad,NSPECV)
	52	REAL QC3(nrad,NSPECV),QC4(nrad,NSPECV)
[3]	53
[175]	54	c---- NUAGES
	55	real TNUABS,TNUSCAT
	56	real rcdb(NLAYER)
	57	real xfrb(NLAYER,4)
	58
[3]	59	save qf1,qf2,qf3,qf4,qm1,qm2,qm3,qm4
	60
[175]	61	integer ilat,jalt
	62	common/toto/ilat,jalt
	63
[3]	64	C*
	65	C THIS SUBROUTINE SETS THE OPTICAL CONSTANTS IN THE VISIBLE
	66	C IT CALCUALTES FOR EACH LAYER, FOR EACH SPECRAL INTERVAL IN THE VIS
	67	C LAYER: WBAR, DTAU, COSBAR
	68	C LEVEL: TAU
	69	C
	70	C ZERO THE COLUMN OPTICAL DEPTHS OF EACH TYPE
	71	C ??FLAG? THE OPTICAL DEPTH OF THE TOP OF THE MODEL
	72	C MAY NOT BE ZERO.
	73
	74	c***** DEBUT DES BOUCLES **********************
	75	DO 100 K=1,NSPECV !b! BOUCLE SUR LAMBDA
	76
	77	TAURV_1pt(K)=0.
	78	TAUHV_1pt(K)=0. ! INTEGRATED TAU.......INITIALIZATION.
	79	TAUCV_1pt(K)=0. ! Rayleigh, Haze, Cloud, Gas
	80	TAUGV_1pt(K)=0. ! sca, abs, abs , abs
	81
	82	DO 100 J=1,NLAYER !a! BOUCLE SUR L"ALTITUDE
[175]	83	jalt=j
[3]	84	C #1: HAZE
	85	c---------------------------
	86
	87	c CALL THE MIE CODE TO GIVE THE AEROSOL PROPERTIES
	88	c USE XFRAC FOR FRACTAL AEROSOLS PROPERTIES AT LAMBDA < 2. um
	89
	90
	91
	92
	93	c /\
	94	c / \
	95	c / \
	96	c / _O \
	97	c / \|/ \
	98	c / / \ \
	99	c / \|\ \/\ \
	100	c / \|\| / \ \
	101	c ----------------
	102	c \| WARNING \|
	103	c \| SLOW DOWN \|
	104	c ----------------
	105
	106
	107
	108
	109	c********* EN TRAVAUX *************************
	110
	111	TAEROS=0.
	112	TAEROSCAT=0.
	113	CBAR=0.
	114
	115	c print*,"rayon=",rayon
	116	c print*,"RF=",RF
	117
	118	DO inq=1,nrad !BOUCLE SUR LES TAILLE D"AEROSOLS
	119
	120
	121	IF (rayon(inq).lt.RF(inq)) THEN ! aerosols spheriques
	122
	123
	124	if(ioptv.eq.0.and.J.eq.1) then
	125	c CALL XMIE(rayon(inq)*1.e6,REALV(K),XIMGV(K),
	126	c & QEXT,QSCT,QABS,QBAR,WNOV(K))
	127
	128	CALL CMIE(1.E-2/WNOV(K),REALV(K),XIMGV(K),rayon(inq),
	129	& QEXT,QSCT,QABS,QBAR)
	130
	131	c print*,'inq=',inq,' QM1=',QM1(inq,K),' QEXT=',QEXT
	132
	133	QM1(inq,K)=QEXT
	134	QM2(inq,K)=QSCT
	135	QM3(inq,K)=QABS
	136	QM4(inq,K)=QBAR
	137	endif
	138
	139	TAEROS=QM1(inq,K)zqaer_1pt(NLAYER+1-J,inq)1.e-4+TAEROS
	140	TAEROSCAT=QM2(inq,K)zqaer_1pt(NLAYER+1-J,inq)1.e-4+TAEROSCAT
	141	CBAR=CBAR+QM4(inq,K)QM2(inq,K)zqaer_1pt(NLAYER+1-J,inq)*1.e-4
	142
	143	ELSE ! aerosols fractals
	144
	145	XMONO=(rayon(inq)/RF(inq))**3.
	146	XRULE=1.
	147
	148	if(XMONO.gt.16384./1.5) then
	149	XRULE=(XMONO/16384.)
	150	XMONO=16384.
	151	endif
	152
	153	if(ioptv.eq.0.and.J.eq.1) then
	154
[175]	155	c CALL OPTFRAC(XMONO,10000./WNOV(K)
	156	c & ,QEXT,QSCT,QABS,QBAR)
[3]	157
[175]	158	CALL CFFFV11(1.e-2/WNOV(K),REALV(K),XIMGV(K),RF(inq),2.
	159	& ,XMONO,QSCT,QEXT,QABS,QBAR)
[3]	160
	161
	162	QF1(inq,K)=QEXT*XRULE
	163	QF2(inq,K)=QSCT*XRULE
	164	QF3(inq,K)=QABS*XRULE
	165	QF4(inq,K)=QBAR
	166
	167	c print*,'inq=',inq,' QF1=',QF1(inq,K),' QEXT=',QEXT,' XRULE=',XRULE
	168
	169	endif
	170
	171	TAEROS=QF1(inq,K)*zqaer_1pt(NLAYER+1-J,inq)+TAEROS
	172	TAEROSCAT=QF2(inq,K)*zqaer_1pt(NLAYER+1-J,inq)+TAEROSCAT
	173	CBAR=CBAR+QF4(inq,K)QF2(inq,K)zqaer_1pt(NLAYER+1-J,inq)
	174
	175	ENDIF
	176
	177
	178	ENDDO ! nrad
	179
	180
	181	CBAR=CBAR/TAEROSCAT
	182
	183	DELTAZ=Z(J)-Z(J+1)
	184
	185	c --------------------------------------------------------------------
	186	c profil brume Pascal: fit T (sauf tropopause) et albedo
	187	c -------------------
	188	if( cutoff.eq.1) then
	189	IF(PRESS(J).gt.9.e-3) THEN
	190	TAEROS=TAEROSM1DELTAZ/DELTAZM10.85
	191	TAEROSCAT=TAEROSCATM1DELTAZ/DELTAZM10.85
	192	c TAEROS=0.
	193	c TAEROSCAT=0.
	194	ENDIF
	195
	196	IF(PRESS(J).gt.1.e-1) THEN
	197	TAEROS=TAEROSM1DELTAZ/DELTAZM11.15
	198	TAEROSCAT=TAEROSCATM1DELTAZ/DELTAZM11.15
	199	c TAEROS=0.
	200	c TAEROSCAT=0.
	201	ENDIF
	202	endif !cutoff=1
	203
	204	c profil brume pour fit T (y compris tropopause), mais ne fit plus albedo...
	205	c -----------------------
	206	if( cutoff.eq.2) then
	207	IF(PRESS(J).gt.1.e-1) THEN
	208	TAEROS=0.
	209	TAEROSCAT=0.
	210	ENDIF
	211	endif !cutoff=2
	212	c --------------------------------------------------------------------
	213
	214	TAEROSM1=TAEROS
	215	TAEROSCATM1=TAEROSCAT
	216	DELTAZM1=DELTAZ
	217
	218
	219	IF (TAEROSCAT.le.0.) CBAR=0.
	220
[175]	221	c if (IPRINT.eq.1) then
	222	c if (k.eq.NSPECV/2) then
	223	c write(*,1699) '@VI',K,J,TAEROS,TAEROSCAT,CBAR
	224	c write(*,1699) '@ ',K,J,QF1(1,K),QF2(1,K),zqaer_1pt(NLAYER+1-J,1)
	225	c write(*,1699) '@ ',K,J,QF1(3,K),QF2(3,K),zqaer_1pt(NLAYER+1-J,3)
	226	c write(*,1699) '@ ',K,J,QF1(5,K),QF2(5,K),zqaer_1pt(NLAYER+1-J,5)
	227	c write(*,1699) '@ ',K,J,QF1(7,K),QF2(7,K),zqaer_1pt(NLAYER+1-J,7)
	228	c write(*,1699) '@ ',K,J,QF1(9,K),QF2(9,K),zqaer_1pt(NLAYER+1-J,9)
	229	c print*
	230	c endif
	231	c endif
[3]	232
[175]	233	1699 FORMAT(a3,2I3,3(ES15.7,1X))
[3]	234
	235	c********* EN TRAVAUX *************************
	236
	237	C #2: RAYLEIGH
	238	c-------------------------------
	239
	240	C RAYLEIGH SCATTERING STRAIGHT FROM HANSEN AND TRAVIS...SEE NOTES
	241	C RATIOED BY THE LAYER COLUMN NUMBER TO THE TOTAL
	242	C COLUMN NUMBER ON EARTH. CM-2
	243	C THIS IS THE SCATTERING BY THE ATMOSPHERE
	244
	245	TAURAY=(COLDEN(J)28.9/(XMU(J)1013.25))*
	246	&(.008569/WLNV(K)*4)(1.+.0113/WLNV(K)2+.00013/WLNV(K)4)
	247
	248	c PRINT*,WLNV(K)
	249	c COLX=0.
	250	c COLP=0.
	251	c COLT=0.
	252	c DO IU=1,NLAYER
	253	c COLP=COLDEN(IU)1.e+11.35+COLP
	254	c TAURAY=(COLDEN(IU)28.9/(XMU(IU)1013.25))*
	255	c & (.008569/WLNV(K)*4)(1.+.0113/WLNV(K)**2
	256	c & +.00013/WLNV(K)**4)
	257	c COLT=COLT+TAURAY
	258	c COLX=COLDEN(IU)1.e+1/(1.E528./22.4E3)1.e-10.0933e-1+COLX
	259	c \|
	260	c \|
	261	c g/cm2->kg/m2 \| m2/kg
	262	c Print*,IU, tauray,
	263	c & COLDEN(IU)1.e+1/(1.E528./22.4E3)1.e-10.543e-1
	264	c ENDDO
	265	c PRINT*,COLP,' PRESSURE AT GROUND;'
	266	c PRINT*,COLX,' TAU_GAS AT GROUND;'
	267	c print*,colt,colx,' COLT, COLX'
	268	c STOP
	269
	270	c DZ=Z(J)-Z(J+1)
	271	c PRINT*, Z(J),WLNV(K),
	272	c &(28.9/(XMU(J)1013.25))(.008569/WLNV(K)*4)
	273	c &(1.+.0113/WLNV(K)2+.00013/WLNV(K)4)
	274	c & ,COLDEN(J)/DZ/100000.,
	275	c &(28.9/(XMU(J)1013.25))(.008569/WLNV(K)*4)
	276	c &(1.+.0113/WLNV(K)2+.00013/WLNV(K)4)
	277	c & *COLDEN(J)/DZ/100000.
	278
	279
	280
	281	C #3: CLOUD
	282	c----------------------------
	283
	284	C NEXT COMPUTE TAU CLOUD
	285
[175]	286	IF (clouds.eq.0) THEN
	287	CNBAR=0.
	288	TNUSCAT=0.
	289	TNUABS=0.
	290	TBNUABS=0.
	291	ELSE
	292	CNBAR=0.
	293	TNUSCAT=0.
	294	TNUABS=0.
	295	TBNUABS=0.
	296	QEXTC=0.
	297	QSCTC=0.
	298	QABSC=0.
	299	CBARC=0.
[3]	300
[175]	301	do inq=1,nrad
	302	QC1(INQ,k)=0.
	303	QC2(INQ,k)=0.
	304	QC3(INQ,k)=0.
	305	QC4(INQ,k)=0.
	306	enddo
[3]	307
[175]	308	IF (rcdb(nlayer+1-J).gt.1.1e-10) THEN
[3]	309
[175]	310	** OPTICAL CONSTANT : MIXING RULES
	311
	312	XNR=xfrb(nlayer+1-J,1)*REALV(K) !
	313	& +xfrb(nlayer+1-J,2)*RCLDV(K) !
	314	& +xfrb(nlayer+1-J,3)*RCLDV2(K) !
	315	& +xfrb(nlayer+1-J,4)*RCLDV2(K) !
	316
	317	XNI=xfrb(nlayer+1-J,1)*XIMGV(K)
	318	& +xfrb(nlayer+1-J,2)*XICLDV(K)
	319	& +xfrb(nlayer+1-J,3)*XICLDV2(K)
	320	& +xfrb(nlayer+1-J,4)*XICLDV2(K)
	321
	322	** OPTICAL CONSTANT : LIQUID DROP = THOLIN
	323	IF(xfrb(nlayer+1-J,1).ge.0.01) THEN
	324	XNI=XIMGV(K)
	325	XNR=REALV(K)
	326	ENDIF
	327
	328	IF (XNI.gt.1.e-10 .and. XNR.gt.1.00) THEN
	329	CALL CMIE(1.E-2/WNOV(K),XNR,XNI,
	330	& rcdb(nlayer+1-J),
	331	& QEXTC,QSCTC,QABSC,CBARC)
	332	ELSE
	333	PRINT*,' WARNING XNR/XNI in optcv: ',XNR,XNI
	334	QEXTC=0.
	335	QSCTC=0.
	336	QABSC=0.
	337	CBARC=0.
	338	STOP
	339	ENDIF
	340	ELSE
	341	QEXTC=0.
	342	QSCTC=0.
	343	QABSC=0.
	344	CBARC=0.
	345	ENDIF
	346
	347	DO inq=1,nrad
	348
	349	QC1(INQ,k)=QEXTC/xnuf
	350	QC2(INQ,k)=QSCTC/xnuf
	351	QC3(INQ,k)=QABSC/xnuf
	352	QC4(INQ,k)=CBARC
	353
	354	TNUABS=QC1(inq,K)zqaer_1pt(NLAYER+1-J,inq+nrad)1.e-4
	355	& +TNUABS
	356
	357	TNUSCAT=QC2(inq,K)zqaer_1pt(NLAYER+1-J,inq+nrad)1.e-4
	358	& +TNUSCAT
	359
	360	CNBAR=QC4(inq,K)QC2(inq,K)
	361	& zqaer_1pt(NLAYER+1-J,inq+nrad)*1.e-4 + CNBAR
	362
	363	ENDDO
	364
	365	IF(TNUSCAT.EQ.0) CNBAR=0.
	366	IF(TNUSCAT.NE.0.) CNBAR=CNBAR/TNUSCAT
	367
	368
	369	ENDIF ! Cond. CLD
	370
	371	TAURV_1pt(K)=TAURV_1pt(K)+TAURAY
	372	TAUGVD_1pt(J,K)=TAURV_1pt(K)
	373
	374	TAUHV_1pt(K)=TAUHV_1pt(K)+TAEROS ! INTEGRATED Quant.
	375	TAUHVD_1pt(J,K)=TAUHV_1pt(K)
	376
	377	TAUCV_1pt(K)=TAUCV_1pt(K)+TNUABS
	378	TAUCVD_1pt(J,K)=TAUCV_1pt(K)
	379
	380
[3]	381	C #4: TAUGAS
	382	C----------------------------
	383
	384	C LOOP OVER THE NTERMS
	385	C THIS IS THE ABSORPTION BY THE ATMOSPHERE (METHANE)
	386
	387
	388	DO 909 NT=1,NTERM(K)
	389	TAUGAS=COLDEN(J)GAS1(J)BTERM(NT,K)*
	390	& ( (PRESS(J+1) + PRESS(J)).5 )*PEXPON(K)
	391
	392
[175]	393	* COSBV ET COSBVP
	394	*-----------------
[3]	395
[175]	396	IF(TAEROSCAT+TNUSCAT+TAURAY .ne. 0.) THEN
	397	COSBV_1pt(J,K,NT)=(CBARTAEROSCAT + CNBARTNUSCAT)
	398	& /(TAEROSCAT+TNUSCAT+TAURAY) !CBAR_RAY=0.
	399	ELSE
	400	COSBV_1pt(J,K,NT)=0.
	401	ENDIF
[3]	402
[175]	403	IF(TAEROSCAT+TAURAY .ne. 0.) THEN
	404	COSBVP_1pt(J,K,NT)=(CBAR*TAEROSCAT)
	405	& /(TAEROSCAT+TAURAY) !CBAR_RAY=0.
	406	ELSE
	407	COSBVP_1pt(J,K,NT)=0.
	408	ENDIF
[3]	409
[175]	410	* DTAUV ET DTAUVP
	411	*-----------------
[3]	412
[175]	413	DTAUV_1pt(J,K,NT) =TAUGAS+TAEROS+TAURAY+TNUABS !TAU_ABS_METH
	414	DTAUVP_1pt(J,K,NT)=TAUGAS+TAEROS+TAURAY !TAU_ABS_METH
[3]	415
[175]	416	TAUGV_1pt(K)=TAUGV_1pt(K)+TAUGAS*ATERM(NT,K) !INTEG.
[3]	417
[175]	418	* WBARV ET WBARVP
	419	*-----------------
[3]	420
[175]	421	IF(TAUGAS+TAEROS+TAURAY+TNUABS .ne. 0.) THEN
	422	WBARV_1pt(J,K,NT)=(TAEROSCAT+TAURAY*0.9999999 + TNUSCAT)
	423	& /(TAUGAS+TAEROS+TAURAY+TNUABS)
	424	ELSE
	425	WBARV_1pt(J,K,NT)=0.
	426	ENDIF
	427
	428	IF(TAUGAS+TAEROS+TAURAY .ne. 0.) THEN
	429	WBARVP_1pt(J,K,NT)=(TAEROSCAT+TAURAY*0.9999999 )
	430	& /(TAUGAS+TAEROS+TAURAY)
	431	ELSE
	432	WBARVP_1pt(J,K,NT)=0.
	433	ENDIF
	434
[3]	435	909 CONTINUE
	436	TAUGVD_1pt(J,K)=TAUGVD_1pt(J,K)+TAUGV_1pt(K)
	437	100 CONTINUE
	438	ioptv=1
	439
	440	c HERE END OF THE LOOPS *******
	441	c******************************
	442
	443	C TOTAL EXTINCTION OPTICAL DEPTHS
	444	DO 119 K=1,NSPECV
	445	C LOOP OVER NTERMS
	446	DO 119 NT=1,NTERM(K)
	447	TAUV_1pt(1,K,NT)=0.0
[175]	448	TAUVP_1pt(1,K,NT)=0.0
[3]	449	DO 119 J=1,NLAYER
	450	TAUV_1pt(J+1,K,NT)=TAUV_1pt(J,K,NT)+DTAUV_1pt(J,K,NT)
[175]	451	TAUVP_1pt(J+1,K,NT)=TAUVP_1pt(J,K,NT)+DTAUVP_1pt(J,K,NT)
[3]	452	119 CONTINUE
	453
	454	c print*,'SETUP'
	455	c do i=1,NSPECV
	456	c print*,WLNV(i)
	457	c do j=1,NLAYER+1
	458	c print*,Z(j),TAUV(1,j,i,1),WBARV(1,j,i,1),COSBV(1,j,i,1)
	459	c enddo
	460	c enddo
	461	c
	462	c IF (IPRINT .GT. 1) THEN
	463	c NT=1
	464	c IF (2 .GT. 1) THEN
	465	c WRITE (6,120)
	466	c 120 FORMAT(///' OPTICAL CONSTANTS IN THE VISIBLE (@EQUATOR) ')
	467	c WRITE(6,*) 'latitude:',ig
	468	c DO 200 K=1,NSPECV
	469	c WRITE (6,190)
	470	c WRITE (6,210)K,WLNV(K),WNOV(K),BWNV(K)
	471	c & ,BWNV(K)+DWNV(K),DWNV(K)
	472	c WRITE (6,230)REALV(K),XIMGV(K)
	473	c DO 195 J=1,NLAYER,NLAYER
	474	C RECALCULATE FOR PRINT OUT ONLY, ONLY FIRST NTERM AT ig=12 (EQUATOR)
	475	c WRITE (6,220)XNUMB(J), WBARV_1pt(J,K,NT),COSBV_1pt(J,K,NT)
	476	c & ,DTAUV_1pt(J,K,NT),TAUV_1pt(J,K,NT)
	477	c 195 CONTINUE
	478	c WRITE (6,240) TAUV_1pt(NLEVEL,K,NT)
	479	c 200 CONTINUE
	480	c END IF
	481
	482	c 210 FORMAT(1X,I3,F10.3,F10.2,F10.2,'-',F8.2,F10.3)
	483	c 190 FORMAT(1X//' SNUM MICRONS WAVENU INTERVAL DELTA-WN')
	484	c 230 FORMAT(1X,'NREAL(LAYER)= ',1PE10.3,' NIMG(LAYER)= ',E10.3/
	485	c &' #AEROSOLS WBAR COSBAR DTAU TAU'
	486	c & ,9X,'RAY GAS AEROSOL')
	487	c 220 FORMAT(8(1X,F9.3))
	488	c 240 FORMAT(41X,F9.3)
	489
	490	if (IPRINT.eq.1) stop
	491
	492	RETURN
	493	END

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

Download in other formats:

Original Format