Context Navigation

new_cloud_sedim.F @ 1704

Last change on this file since 1704 was 1687, checked in by slebonnois, 8 years ago
SL: bugs corrections outputs/clouds_AStol/upper_atmosphere/start2archive
Property svn:executable set to ``*
File size: 16.1 KB

Rev	Line
[1305]	1	SUBROUTINE new_cloud_sedim(n_lon,n_lev,ptimestep,
	2	& pmidlay,pbndlay,
[1442]	3	& pt,
	4	& pq, pdqsed,pdqs_sed,nq,F_sed)
[1305]	5
	6	USE ioipsl
	7	USE dimphy
	8	USE chemparam_mod
	9	IMPLICIT NONE
	10
	11	c=======================================================================
	12	c
	13	c=======================================================================
	14
	15	c-----------------------------------------------------------------------
	16	c declarations:
	17	c -------------
	18	#include "YOMCST.h"
	19	c#include "dimphys.h"
	20	c#include "comcstfi.h"
	21	c#include "tracer.h"
	22	c#include "callkeys.h"
	23
	24	c
	25	c arguments:
	26	c ----------
	27
	28	INTEGER n_lon ! number of horizontal grid points
	29	INTEGER n_lev ! number of atmospheric layers
	30	REAL ptimestep ! physics time step (s)
	31	REAL pmidlay(n_lon,n_lev) ! pressure at middle layers (Pa)
	32	REAL pt(n_lon,n_lev) ! temperature at mid-layer (l)
	33	REAL pbndlay(n_lon,n_lev+1) ! pressure at layer boundaries
	34
	35	c Traceurs :
	36	integer nq ! number of tracers
	37	real pq(n_lon,n_lev,nq) ! tracers (kg/kg)
	38	c real pdqfi(n_lon,n_lev,nq) ! tendency before sedimentation (kg/kg.s-1)
	39	real pdqsed(n_lon,n_lev,2) ! tendency due to sedimentation (kg/kg)
	40	real pdqs_sed(n_lon) ! surface density (Flux if /ptimestep) at surface due to sedimentation (kg.m-2)
	41
	42	c local:
	43	c ------
[1442]	44	integer imode
[1305]	45	integer ig
	46	integer iq
	47	integer l
	48
	49	real zlev(n_lon,n_lev+1) ! altitude at layer boundaries
	50	real zlay(n_lon,n_lev) ! altitude at the midlle layer
	51	real zqi_wv(n_lon,n_lev) ! to locally store H2O tracer
	52	real zqi_sa(n_lon,n_lev) ! to locally store H2SO4 tracer
[1442]	53	real m_lay (n_lon,n_lev) ! Layer Pressure over gravity (Dp/g == kg.m-2)
[1305]	54	real wq(n_lon,n_lev+1) ! displaced tracer mass (kg.m-2)
	55
	56	c Physical constant
	57	c ~~~~~~~~~~~~~~~~~
	58	c Gas molecular viscosity (N.s.m-2)
[1442]	59	c real,parameter :: visc=1.e-5 ! CO2
	60	REAL :: VISCOSITY_CO2
[1305]	61	c Effective gas molecular radius (m)
	62	real,parameter :: molrad=2.2e-10 ! CO2
	63
[1687]	64	c Ratio radius shell model du mode 3
	65	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	66	c Ce ratio correspond aux mesures effectuées par J. Cimino (1982), Icarus
	67	c Fixer ce parametre a 0 revient a une gouttelette pure en liquide acide sulfurique
	68	c ATTENTION ! DOIT ETRE COHERENT AVEC new_cloud_venus !
	69	REAL, PARAMETER :: qrad = 0.97
	70	REAL :: qmass
	71	c masse volumique du coeur (kg.m-3)
	72	c ATTENTION ! DOIT ETRE COHERENT AVEC new_cloud_venus !
	73	REAL, PARAMETER :: rho_core = 2500.0
[1305]	74
	75	REAL, DIMENSION(n_lon,n_lev+1) ::
	76	+ wgt_SA ! Fraction of H2SO4 in droplet local
	77
	78	c Stokes speed and sedimentation flux variable
	79	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	80
	81	REAL :: A1,A2,A3,A4, ! coeff du DL du Flux de sedimentation
	82	+ D_stokes, ! coeff de la vitesse de Stokes
	83	+ Rp_DL, ! "Point" du DL
	84	+ l_mean, ! libre parcours moyen (m)
	85	+ a,b_exp,c ! coeff du calcul du Flux de sedimentation
	86	REAL, DIMENSION(n_lon,n_lev+1) ::
	87	+ F_sed ! Flux de sedimentation (kg.m-2.s-1 puis en output kg.m-2)
	88
	89
	90	REAL :: R_mode0 ! Rayon mode 0 (m), rayon le plus frequent
	91
	92
	93
[1442]	94	! PRINT*,'RHO_DROPLET new_cloud_sedim.F'
	95	! PRINT*,'rho_droplet',rho_droplet(16,21)
	96	! PRINT*,'T',pt(16,21),'WSA',WH2SO4(16,21)
	97
[1305]	98	c-----------------------------------------------------------------------
	99	c 1. Initialization
	100	c -----------------
	101
	102	c Updating the droplet mass mixing ratio with the partition H2O/H2SO4
	103	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	104
	105	do l=1,n_lev
	106	do ig=1,n_lon
	107	zqi_wv(ig,l) = pq(ig,l,i_h2oliq)
	108	zqi_sa(ig,l) = pq(ig,l,i_h2so4liq)
[1442]	109	wgt_SA(ig,l) = WH2SO4(ig,l)
[1305]	110	enddo
	111	enddo
	112
[1442]	113	c Init F_sed
	114	F_sed(:,:) = 0.0E+0
	115
	116	c Au niveau top+1 , tout égal a 0
	117	wgt_SA(:,n_lev+1) = 0.0E+0
	118
[1305]	119	c Computing the different layer properties
	120	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[1442]	121	c m_lay (kg.m-2)
[1305]	122	c Ici g=8.87, conflit pour g entre #include "YOMCST.h"
	123	c et #include "comcstfi.h"
	124
	125	do l=1,n_lev
	126	do ig=1, n_lon
[1442]	127	m_lay(ig,l)=(pbndlay(ig,l) - pbndlay(ig,l+1)) /8.87E+0
[1305]	128	IF (m_lay(ig,l).LE.0.0) THEN
	129	PRINT*,'!!!! STOP PROBLEME SEDIMENTATION!!!!'
	130	PRINT*,'!!!! m_lay <= 0 !!!!'
	131	PRINT*,'!!!! STOP PROBLEME SEDIMENTATION!!!!'
	132	ENDIF
	133	end do
	134	end do
	135
	136	c Computing sedimentation for droplet "tracer"
	137	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	138	c pbndlay(:,51)=0 (en parallèle c'est sûr), ne pas l'utiliser pour Fse
	139
[1687]	140	c Sedimentation pour une gouttelette mode 3 de type J. Cimino, 1982, Icarus
	141	c c.a.d 97% radius due a solide 3% radius acide sulfurique
	142	DO imode=1, nbr_mode - 1
[1442]	143	DO l = cloudmin, cloudmax
	144	DO ig=1,n_lon
[1305]	145
	146	c RD=1000.RNAVORKBOL/RMD avec RMD=43.44 Masse molaire atm venus en g.mol-1
[1442]	147	D_stokes=((rho_droplet(ig,l)-pmidlay(ig,l)/(RD*pt(ig,l))))
	148	& (2./9.)(RG/VISCOSITY_CO2(pt(ig,l)))
[1305]	149
	150	l_mean=(pt(ig,l)/pmidlay(ig,l))*
	151	& (0.707R/(4.RPI* molradmolrad RNAVO))
	152
[1442]	153	R_mode0=R_MEDIAN(ig,l,imode)*
	154	& EXP(-LOG(STDDEV(ig,l,imode))**2.)
	155	IF ((l_mean/(R_mode0)).GT.10.) THEN
	156	Rp_DL=R_MEDIAN(ig,l,imode)*
	157	& EXP(3.LOG(STDDEV(ig,l,imode))*2.)
[1305]	158	ELSE
[1442]	159	Rp_DL=R_MEDIAN(ig,l,imode)*
	160	& EXP(4.LOG(STDDEV(ig,l,imode))*2.)
[1305]	161	ENDIF
	162
	163	a=1.246*l_mean
	164
	165	c=0.87/l_mean
	166
	167	b_exp=0.42l_meanEXP(-c*Rp_DL)
	168
	169	A1=a+b_exp(1.+cRp_DL
	170	& +0.5(Rp_DLc)**2
	171	& +1./6.(Rp_DLc)**3)
	172
	173	A2=1.-b_exp*(c
	174	& +Rp_DLc*2
[1442]	175	& +0.5(Rp_DL2)(c**3))
[1305]	176
	177	A3=0.5b_exp(c*2+Rp_DLc**3)
	178
	179	A4=-b_exp1./6.c**3
[1442]	180
	181	c Addition des Flux de tous les modes presents
	182	F_sed(ig,l)=F_sed(ig,l)+(rho_droplet(ig,l)4./3.RPI*
	183	& NBRTOT(ig,l,imode)1.0E6D_stokes*(
	184	& A1R_MEDIAN(ig,l,imode)*4
	185	& EXP(8.0LOG(STDDEV(ig,l,imode))**2.)
	186	& +A2R_MEDIAN(ig,l,imode)*5
	187	& EXP(12.5LOG(STDDEV(ig,l,imode))**2.)
	188	& +A3R_MEDIAN(ig,l,imode)*6
	189	& EXP(18.0LOG(STDDEV(ig,l,imode))**2.)
	190	& +A4R_MEDIAN(ig,l,imode)*7
	191	& EXP(24.5LOG(STDDEV(ig,l,imode))**2.)))
[1305]	192
	193	c PRINT*,' APRES dTime: F_sed=',F_sed(ig,l), ig, l
	194
[1442]	195	IF (F_sed(ig,l).GT.m_lay(ig,l)) THEN
	196	PRINT*,'==============================================='
	197	PRINT*,'WARNING On a epuise la couche', ig, l
	198	PRINT*,'On epuise pas une couche avec une espèce
	199	& minoritaire, c est pas bien maaaaaal'
	200	PRINT*,'Water',zqi_wv(ig,l),'Sulfuric Acid',zqi_sa(ig,l)
	201	PRINT*,'F_sed:',F_sed(ig,l),'m_lay:',m_lay(ig,l)
	202	PRINT*,'F_sed/dtphy',F_sed(ig,l)/ptimestep
	203	PRINT*,'Pbnd top',pbndlay(ig,l+1),'Temp',pt(ig,l),'Rho',
	204	& rho_droplet(ig,l)
	205	PRINT*,'Ntot',NBRTOT(ig,l,:)
	206	PRINT*,'StdDev',STDDEV(ig,l,:),'Rmed',R_MEDIAN(ig,l,:)
	207	PRINT*,'K_MASS',K_MASS(ig,l,:)
	208	PRINT*,'WSA',WH2SO4(ig,l),'RHO',rho_droplet(ig,l)
[1305]	209
	210	c ELSE
	211	c
	212	c PRINT*,'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
	213	c PRINT*,'WARNING On a PAS epuise la couche', ig, l
	214	c PRINT*,'F_sed:',F_sed(ig,l),'m_lay:',m_lay(ig,l)
	215	c PRINT*,'F_sed/dtphy',F_sed(ig,l)/ptimestep
	216	c PRINT*,'Pbnd top',pbndlay(ig,l+1),'Temp',pt(ig,l),'Rho',
[1442]	217	c & rho_droplet(ig,l)(ig,l)
	218	c PRINT,'Ntot',NBRTOT(ig,l),'Ntot m3',NBRTOT(ig,l)1.0e6
	219	c PRINT*,'StdDev',STDDEV(ig,l),'Rmed',R_MEDIAN(ig,l)
	220	STOP
	221	ENDIF
[1305]	222
[1442]	223	IF (F_sed(ig,l).LT.0.0d0) THEN
[1305]	224	PRINT*,"F_sed est négatif !!!"
	225	PRINT*,'F_sed:',F_sed(ig,l),'m_lay:',m_lay(ig,l)
	226	PRINT*,'F_sed/dtphy',F_sed(ig,l)/ptimestep
	227	PRINT*,'Pbnd top',pbndlay(ig,l+1),'Pmid',pmidlay(ig,l)
	228	PRINT*,'Temp',pt(ig,l),'Rho',
[1442]	229	& rho_droplet(ig,l)
	230	PRINT*,'Ntot',NBRTOT(ig,l,imode),'Ntot m3',
	231	& NBRTOT(ig,l,imode)*1.0e6
	232	PRINT*,'StdDev',STDDEV(ig,l,imode),'Rmed',
	233	& R_MEDIAN(ig,l,imode)
[1305]	234	PRINT*,'A1',A1,'A2',A2
	235	PRINT*,'A3',A1,'A4',A2
	236	PRINT*,'D_stokes',D_stokes
	237	STOP
	238	ENDIF
[1442]	239
	240	ENDDO
	241
	242	c ELSE
	243	c F_sed(:,l)=0.0d0
	244	c ENDIF
	245
[1305]	246	ENDDO
	247	ENDDO
	248
[1687]	249	c****************************************************************
	250	c On calcule le F_sed du mode 3 + coeff*(Fsed1 + Fsed2)
	251	c****************************************************************
	252	DO l = cloudmin, cloudmax
	253	DO ig=1,n_lon
	254
	255	c calcul de qmass
	256	qmass=(rho_coreqrad*3)/
	257	& (rho_coreqrad3+rho_droplet(ig,l)(1.-qrad**3))
	258
	259	c RD=1000.RNAVORKBOL/RMD avec RMD=43.44 Masse molaire atm venus en g.mol-1
	260	D_stokes=(((qmassrho_core+(1.-qmass)rho_droplet(ig,l))
	261	& -pmidlay(ig,l)/(RD*pt(ig,l))))
	262	& (2./9.)(RG/VISCOSITY_CO2(pt(ig,l)))
	263
	264	l_mean=(pt(ig,l)/pmidlay(ig,l))*
	265	& (0.707R/(4.RPI* molradmolrad RNAVO))
	266
	267	R_mode0=R_MEDIAN(ig,l,3)*
	268	& EXP(-LOG(STDDEV(ig,l,3))**2.)
	269	IF ((l_mean/(R_mode0)).GT.10.) THEN
	270	Rp_DL=R_MEDIAN(ig,l,3)*
	271	& EXP(3.LOG(STDDEV(ig,l,3))*2.)
	272	ELSE
	273	Rp_DL=R_MEDIAN(ig,l,3)*
	274	& EXP(4.LOG(STDDEV(ig,l,3))*2.)
	275	ENDIF
	276
	277	a=1.246*l_mean
	278
	279	c=0.87/l_mean
	280
	281	b_exp=0.42l_meanEXP(-c*Rp_DL)
	282
	283	A1=a+b_exp(1.+cRp_DL
	284	& +0.5(Rp_DLc)**2
	285	& +1./6.(Rp_DLc)**3)
	286
	287	A2=1.-b_exp*(c
	288	& +Rp_DLc*2
	289	& +0.5(Rp_DL2)(c**3))
	290
	291	A3=0.5b_exp(c*2+Rp_DLc**3)
	292
	293	A4=-b_exp1./6.c**3
	294
	295	c Addition des Flux de tous les modes presents
	296	F_sed(ig,l)=F_sed(ig,l)
	297	& +((1.-qmass)/(1.-qmassK_MASS(ig,l,3)))(
	298	& (qmassrho_core+(1.-qmass)rho_droplet(ig,l))4./3.RPI*
	299	& NBRTOT(ig,l,3)1.0E6D_stokes*(
	300	& A1R_MEDIAN(ig,l,3)*4
	301	& EXP(8.0LOG(STDDEV(ig,l,3))**2.)
	302	& +A2R_MEDIAN(ig,l,3)*5
	303	& EXP(12.5LOG(STDDEV(ig,l,3))**2.)
	304	& +A3R_MEDIAN(ig,l,3)*6
	305	& EXP(18.0LOG(STDDEV(ig,l,3))**2.)
	306	& +A4R_MEDIAN(ig,l,3)*7
	307	& EXP(24.5LOG(STDDEV(ig,l,3))**2.)))
	308
	309	c PRINT*,' APRES dTime: F_sed=',F_sed(ig,l), ig, l
	310
	311	IF (F_sed(ig,l).GT.m_lay(ig,l)) THEN
	312	PRINT*,'==============================================='
	313	PRINT*,'WARNING On a epuise la couche', ig, l
	314	PRINT*,'On epuise pas une couche avec une espèce
	315	& minoritaire, c est pas bien maaaaaal'
	316	PRINT*,'Water',zqi_wv(ig,l),'Sulfuric Acid',zqi_sa(ig,l)
	317	PRINT*,'F_sed:',F_sed(ig,l),'m_lay:',m_lay(ig,l)
	318	PRINT*,'F_sed/dtphy',F_sed(ig,l)/ptimestep
	319	PRINT*,'Pbnd top',pbndlay(ig,l+1),'Temp',pt(ig,l),'Rho',
	320	& rho_droplet(ig,l)
	321	PRINT*,'Ntot',NBRTOT(ig,l,:)
	322	PRINT*,'StdDev',STDDEV(ig,l,:),'Rmed',R_MEDIAN(ig,l,:)
	323	PRINT*,'K_MASS',K_MASS(ig,l,:)
	324	PRINT*,'WSA',WH2SO4(ig,l),'RHO',rho_droplet(ig,l)
	325
	326	c ELSE
	327	c
	328	c PRINT*,'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
	329	c PRINT*,'WARNING On a PAS epuise la couche', ig, l
	330	c PRINT*,'F_sed:',F_sed(ig,l),'m_lay:',m_lay(ig,l)
	331	c PRINT*,'F_sed/dtphy',F_sed(ig,l)/ptimestep
	332	c PRINT*,'Pbnd top',pbndlay(ig,l+1),'Temp',pt(ig,l),'Rho',
	333	c & rho_droplet(ig,l)(ig,l)
	334	c PRINT,'Ntot',NBRTOT(ig,l),'Ntot m3',NBRTOT(ig,l)1.0e6
	335	c PRINT*,'StdDev',STDDEV(ig,l),'Rmed',R_MEDIAN(ig,l)
	336	STOP
	337	ENDIF
	338
	339	IF (F_sed(ig,l).LT.0.0d0) THEN
	340	PRINT*,"F_sed est négatif !!!"
	341	PRINT*,'F_sed:',F_sed(ig,l),'m_lay:',m_lay(ig,l)
	342	PRINT*,'F_sed/dtphy',F_sed(ig,l)/ptimestep
	343	PRINT*,'Pbnd top',pbndlay(ig,l+1),'Pmid',pmidlay(ig,l)
	344	PRINT*,'Temp',pt(ig,l),'Rho',
	345	& rho_droplet(ig,l)
	346	PRINT*,'Ntot',NBRTOT(ig,l,imode),'Ntot m3',
	347	& NBRTOT(ig,l,imode)*1.0e6
	348	PRINT*,'StdDev',STDDEV(ig,l,imode),'Rmed',
	349	& R_MEDIAN(ig,l,imode)
	350	PRINT*,'A1',A1,'A2',A2
	351	PRINT*,'A3',A1,'A4',A2
	352	PRINT*,'D_stokes',D_stokes
	353	STOP
	354	ENDIF
	355
	356	ENDDO
	357
	358	c ELSE
	359	c F_sed(:,l)=0.0d0
	360	c ENDIF
	361
	362	ENDDO
	363	c****************************************************************
	364
[1442]	365	c Passage du Flux au Flux pour un pas de temps (== kg.m-2)
	366	F_sed(:,:)=F_sed(:,:)*ptimestep
[1305]	367
[1442]	368
[1305]	369	c VENUS: le flux à la surface est fixé à 0
	370	c les conditions P/T en surface ne permettent pas la condensation
	371	DO ig=1,n_lon
	372	pdqs_sed(ig) = 0.0d0
	373	ENDDO
	374
	375	c Compute the final tendency:
	376	c ---------------------------
	377
	378	c Partie H2SO4l
	379	c ~~~~~~~~~~~~
	380
	381	DO l = 1, n_lev
	382	DO ig=1,n_lon
	383	zqi_sa(ig,l) = zqi_sa(ig,l) + (
	384	& F_sed(ig,l+1)*wgt_SA(ig,l+1)
	385	& - F_sed(ig,l)*wgt_SA(ig,l))
[1442]	386	& / m_lay(ig,l)
	387	c On peut avoir theoriquement le cas ou on epuise tout le VMR present
	388	IF (zqi_sa(ig,l).LT.0.0D0) THEN
[1639]	389	c PRINT*,'STOP sedimentation on epuise tout le VMR present'
	390	c PRINT*,'couche',ig,'level',l
[1442]	391	c STOP
	392	c Ce n est pas juste mais il faudrait alors adapter les pas
	393	c de tps de la phys, microphys et chimie
	394	c car dans ce cas, c est comme si on epuisait la couche pour un pdtphys
	395	c mais en fait on l epuise pour un pdt<pdtphys
	396	zqi_sa(ig,l) = 0.0D0
	397	ENDIF
[1305]	398	pdqsed(ig,l,1) = zqi_sa(ig,l) - pq(ig,l,i_h2so4liq)
	399	ENDDO
	400	ENDDO
	401
	402	c Partie H2Ol
	403	c ~~~~~~~~~~~
	404
	405	DO l = 1, n_lev
	406	DO ig=1,n_lon
	407	zqi_wv(ig,l) = zqi_wv(ig,l) + (
	408	& F_sed(ig,l+1)*(1. - wgt_SA(ig,l+1))
	409	& - F_sed(ig,l)*(1. - wgt_SA(ig,l)))
	410	& / m_lay(ig,l)
[1442]	411	c On peut avoir theoriquement le cas ou on epuise tout le VMR present
	412	IF (zqi_wv(ig,l).LT.0.0D0) THEN
[1639]	413	c PRINT*,'STOP sedimentation on epuise tout le VMR present'
	414	c PRINT*,'couche',ig,'level',l
[1442]	415	c STOP
	416	c Ce n est pas juste mais il faudrait alors adapter les pas
	417	c de tps de la phys, microphys et chimie
	418	c car dans ce cas, c est comme si on epuisait la couche pour un pdtphys
	419	c mais en fait on l epuise pour un pdt<pdtphys
	420	zqi_wv(ig,l) = 0.0D0
	421	ENDIF
[1305]	422	pdqsed(ig,l,2) = zqi_wv(ig,l) - pq(ig,l,i_h2oliq)
	423	ENDDO
	424	ENDDO
	425
	426	c Save output file in 1D model
	427	c ============================
	428
	429	c IF (n_lon .EQ. 1) THEN
	430	c PRINT*,'Save output sedim'
	431	c DO l = 1, n_lev
	432	c DO ig=1,n_lon
	433	c WRITE(77,"(i4,','11(e15.8,','))") l,pdqsed(ig,l),zqi(ig,l),
[1442]	434	c & (WH2SO4(ig,l)*pq(ig,l,i_h2so4liq)+
	435	c & (1.-WH2SO4(ig,l))*pq(ig,l,i_h2oliq)),
[1305]	436	c & pq(ig,l,i_h2so4liq),pq(ig,l,i_h2oliq)
	437	c ENDDO
	438	c ENDDO
	439	c ENDIF
	440
	441	RETURN
	442	END
	443
[1442]	444	*******************************************************************************
	445	REAL FUNCTION VISCOSITY_CO2(temp)
	446	c Aurélien Stolzenbach 2015
	447	c Calcul de la viscosité dynamique du CO2 80°K -> 300°K
	448	c Viscosité dynamique en Pa.s
	449	c Source: Johnston & Grilly (1942)
	450
	451	c température en °K
	452	REAL, INTENT(IN) :: temp
	453
	454	REAL :: denom, numer
	455
	456	c Calcul de la viscosité dynamique grâce à la formule de Jones (Lennard-Jones (1924))
	457
	458	numer = 200.**(2.27/4.27)-0.435
	459	denom = temp**(2.27/4.27)-0.435
	460
	461	VISCOSITY_CO2 = (numer/denom)1015.(temp/200.)**(3./2.)
	462
	463	c convertion de Poises*1e7 -> Pa.s
	464	VISCOSITY_CO2 = VISCOSITY_CO2*1.e-8
	465
	466	END FUNCTION VISCOSITY_CO2
	467	*******************************************************************************
	468
	469

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

Download in other formats:

Original Format