Context Navigation

mzescape.F @ 455

Last change on this file since 455 was 414, checked in by aslmd, 14 years ago

LMDZ.MARS : NEW NLTE MODEL FROM GRANADA AMIGOS

23/11/11 == FGG + MALV

New parameterization of the NLTE 15 micron cooling. The old parameterization is kept as an option, including or not variable atomic oxygen concentration. A new flag is introduced in callphys.def, nltemodel, to select which parameterization wants to be used (new one, old one with variable [O], or old one with fixed [O], see below). Includes many new subroutines and commons in phymars. Some existing routines are also modified:

-physiq.F. Call to the new subroutine NLTE_leedat in first call. Call to nltecool modified to allow for variable atomic oxygen. Depending on the value of nltemodel, the new subroutine NLTEdlvr09_TCOOL is called instead of nltecool.

-inifis.F. Reading of nltemodel is added.

-callkeys.h Declaration of nltemodel is added.

The following lines should be added to callphys.def (ideally after setting callnlte):

# NLTE 15um scheme to use.
# 0-> Old scheme, static oxygen
# 1-> Old scheme, dynamic oxygen
# 2-> New scheme
nltemodel = 2

A new directory, NLTEDAT, has to be included in datagcm.

Improvements into NLTE NIR heating parameterization to take into account variability with O/CO2 ratio and SZA. A new subroutine, NIR_leedat.F, and a new common, NIRdata.h, are included. A new flag, nircorr, is added in callphys.def, to include or not these corrections. The following files are modified:

-nirco2abs.F: nq and pq are added in the arguments. The corrections factors are interpolated to the GCM grid and included in the clculation. A new subroutine, interpnir, is added at the end of the file.

-physiq.F: Call to NIR_leedat added at first call. Modified call to nirco2abs

-inifis: Reading new flag nircorr.

-callkeys.h: Declaration of nircorr.

The following lines have to be added to callphys.def (ideally after callnirco2):

# NIR NLTE correction for variable SZA and O/CO2?
# matters only if callnirco2=T
# 0-> no correction
# 1-> include correction
nircorr=1

A new data file, NIRcorrection_feb2011.dat, has to be included in datagcm.

Small changes to the molecular diffusion scheme to fix the number of species considered, to avoid problems when compiling with more than 15 tracers (for example, when CH4 is included). Modified subroutines: aeronomars/moldiff.F and aeronomars/moldiffcoeff.F

File size: 25.9 KB

Rev	Line
[414]	1	c***********************************************************************
	2	c mzescape.f
	3	c***********************************************************************
	4	c
	5	c program for calculating atmospheric escape functions, from
	6	c a calculation of transmittances and derivatives of these ones
	7
	8	subroutine mzescape(ig,taustar,tauinf,tauii, ib,isot, iirw,iimu)
	9
	10	c jul 2011 malv+fgg adapted to LMD-MGCM
	11	c nov 99 malv adapt mztf to compute taustar (pg.23b-ma
	12	c nov 98 malv allow for overlaping in the lorentz line
	13	c jan 98 malv version for mz1d. based on curtis/mztf.for
	14	c 17-jul-96 mlp&crs change the calculation of mr.
	15	c evitar: divide por cero. anhadiendo: ff
	16	c oct-92 malv correct s(t) dependence for all histogr bands
	17	c june-92 malv proper lower levels for laser bands
	18	c may-92 malv new temperature dependence for laser bands
	19	c @ 991 malv boxing for the averaged absorber amount and t
	20	c ? malv extension up to 200 km altitude in mars
	21	c 13-nov-86 mlp include the temperature weighted to match
	22	c the eqw in the strong doppler limit.
	23	c***********************************************************************
	24
	25	implicit none
	26
	27	include 'nltedefs.h'
	28	include 'nlte_atm.h'
	29	include 'nlte_data.h'
	30	include 'nlte_curtis.h'
	31	include 'tcr_15um.h'
	32	include 'nlte_results.h'
	33
	34
	35	c arguments
	36	integer ig ! ADDED FOR TRACEBACK
	37	real*8 taustar(nl) ! o
	38	real*8 tauinf(nl) ! o
	39	real*8 tauii(nl) ! o
	40	integer ib ! i
	41	integer isot ! i
	42	integer iirw ! i
	43	integer iimu ! i
	44
	45
	46	c local variables and constants
	47	integer i, in, ir, im, k,j
	48	integer nmu
	49	parameter (nmu = 8)
	50	! real*8 tauinf(nl)
	51	real*8 con(nzy), coninf
	52	real*8 c1, c2, ccc
	53	real*8 t1, t2
	54	real*8 p1, p2
	55	real*8 mr1, mr2
	56	real*8 st1, st2
	57	real*8 c1box(70), c2box(70)
	58	real*8 ff ! to avoid too small numbers
	59	real*8 tvtbs(nzy)
	60	real*8 st, beta, ts, eqwmu
	61	real*8 mu(nmu), amu(nmu)
	62	real*8 zld(nl), zyd(nzy)
	63	real*8 correc
	64	real deltanux ! width of vib-rot band (cm-1)
	65	character isotcode*2
	66	real*8 maximum
	67	real*8 csL, psL, Desp, wsL ! for Strong Lorentz limit
	68
	69	c formats
	70	111 format(a1)
	71	112 format(a2)
	72	101 format(i1)
	73	202 format(i2)
	74	180 format(a80)
	75	181 format(a80)
	76	c***********************************************************************
	77
	78	c some needed values
	79	! rl=sqrt(log(2.d0))
	80	! pi2 = 3.14159265358989d0
	81	beta = 1.8d0
	82	! imrco = 0.9865
	83
	84	c esto es para que las subroutines de mztfsub calculen we
	85	c de la forma apropiada para mztf, no para fot
	86	icls=icls_mztf
	87
	88	c codigos para filenames
	89	! if (isot .eq. 1) isotcode = '26'
	90	! if (isot .eq. 2) isotcode = '28'
	91	! if (isot .eq. 3) isotcode = '36'
	92	! if (isot .eq. 4) isotcode = '27'
	93	! if (isot .eq. 5) isotcode = '62'
	94	! if(ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5
	95	! @ .or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9) then
	96	!! encode(2,101,ibcode1)ib
	97	! write ( ibcode1, 101) ib
	98	! else
	99	!! encode(2,202,ibcode2)ib
	100	! write (ibcode2, 202) ib
	101	! endif
	102	! write (*,'( 30h calculating curtis matrix : ,2x,
	103	! @ 8h band = ,i2,2x, 11h isotope = ,i2)') ib, isot
	104
	105	c integration in angle !!!!!!!!!!!!!!!!!!!!
	106
	107	c------- diffusivity approx.
	108	if (iimu.eq.1) then
	109	! write (,) ' diffusivity approx. beta = ',beta
	110	mu(1) = 1.0d0
	111	amu(1)= 1.0d0
	112	c-------data for 8 points integration
	113	elseif (iimu.eq.4) then
	114	write (,)' 4 points for the gauss-legendre angle quadrature.'
	115	mu(1)=(1.0d0+0.339981043584856)/2.0d0
	116	mu(2)=(1.0d0-0.339981043584856)/2.0d0
	117	mu(3)=(1.0d0+0.861136311594053)/2.0d0
	118	mu(4)=(1.0d0-0.861136311594053)/2.0d0
	119	amu(1)=0.652145154862546
	120	amu(2)=amu(1)
	121	amu(3)=0.347854845137454
	122	amu(4)=amu(3)
	123	beta=1.0d0
	124	c-------data for 8 points integration
	125	elseif(iimu.eq.8) then
	126	write (,)' 8 points for the gauss-legendre angle quadrature.'
	127	mu(1)=(1.0d0+0.183434642495650)/2.0d0
	128	mu(2)=(1.0d0-0.183434642495650)/2.0d0
	129	mu(3)=(1.0d0+0.525532409916329)/2.0d0
	130	mu(4)=(1.0d0-0.525532409916329)/2.0d0
	131	mu(5)=(1.0d0+0.796666477413627)/2.0d0
	132	mu(6)=(1.0d0-0.796666477413627)/2.0d0
	133	mu(7)=(1.0d0+0.960289856497536)/2.0d0
	134	mu(8)=(1.0d0-0.960289856497536)/2.0d0
	135	amu(1)=0.362683783378362
	136	amu(2)=amu(1)
	137	amu(3)=0.313706645877887
	138	amu(4)=amu(3)
	139	amu(5)=0.222381034453374
	140	amu(6)=amu(5)
	141	amu(7)=0.101228536290376
	142	amu(8)=amu(7)
	143	beta=1.0d0
	144	end if
	145	c!!!!!!!!!!!!!!!!!!!!!!!
	146
	147	ccc
	148	ccc determine abundances included in the absorber amount
	149	ccc
	150
	151	c first, set up the grid ready for interpolation.
	152	do i=1,nzy
	153	zyd(i) = dble(zy(i))
	154	enddo
	155	do i=1,nl
	156	zld(i) = dble(zl(i))
	157	enddo
	158
	159	c vibr. temp of the bending mode :
	160	if (isot.eq.1) call interdp ( tvtbs,zyd,nzy, v626t1,zld,nl, 1 )
	161	if (isot.eq.2) call interdp ( tvtbs,zyd,nzy, v628t1,zld,nl, 1 )
	162	if (isot.eq.3) call interdp ( tvtbs,zyd,nzy, v636t1,zld,nl, 1 )
	163	if (isot.eq.4) call interdp ( tvtbs,zyd,nzy, v627t1,zld,nl, 1 )
	164
	165	c 2nd: correccion a la n10(i) (cantidad de absorbente en el lower state)
	166	c por similitud a la que se hace en cza.for
	167
	168	do i=1,nzy
	169	if (isot.eq.5) then
	170	con(i) = dble( coy(i) * imrco )
	171	else
	172	con(i) = dble( co2y(i) * imr(isot) )
	173	correc = 2.d0 * dexp( dble(-ee*elow(isot,2))/tvtbs(i) )
	174	con(i) = con(i) * ( 1.d0 - correc )
	175	endif
	176	c-----------------------------------------------------------------------
	177	c mlp & cristina. 17 july 1996
	178	c change the calculation of mr. it is used for calculating partial press
	179	c alpha = alpha(self,co2)pp +alpha(n2)(pt-pp)
	180	c for an isotope, if mr is obtained by co2*imr(iso)/nt we are considerin
	181	c collisions with other co2 isotopes (including the major one, 626)
	182	c as if they were with n2. assuming mr as co2/nt, we consider collisions
	183	c of type 628-626 as of 626-626 instead of as 626-n2.
	184	c mrx(i)=con(i)/ntx(i) ! old malv
	185
	186	! mrx(i)= dble(co2x(i)/ntx(i)) ! mlp & crs
	187
	188	c jan 98:
	189	c esta modif de mlp implica anular el correc (deberia revisar esto)
	190	mr(i) = dble(co2y(i)/nty(i)) ! malv, jan 98
	191
	192	c-----------------------------------------------------------------------
	193
	194	end do
	195
	196	! como beta y 1.d5 son comunes a todas las weighted absorber amounts,
	197	! los simplificamos:
	198	! coninf = beta * 1.d5 * dble( con(n) / log( con(n-1) / con(n) ) )
	199	coninf = dble( con(nzy) / log( con(nzy-1) / con(nzy) ) )
	200
	201	! write (,) ' coninf =', coninf
	202
	203	ccc
	204	ccc temp dependence of the band strength and
	205	ccc nlte correction factor for the absorber amount
	206	ccc
	207	call mztf_correccion ( coninf, con, ib, isot, 0 )
	208
	209	ccc
	210	ccc reads histogrammed spectral data (strength for lte and vmr=1)
	211	ccc
	212	!hfile1 = dirspec//'hi'//dn ! Ya no distinguimos entre d/n
	213	!! hfile1 = dirspec//'hid' ! (see why in his.for)
	214	! hfile1='hid'
	215	!! if (ib.eq.13 .or. ib.eq.14 ) hfile1 = dirspec//'his'
	216	! if (ib.eq.13 .or. ib.eq.14 ) hfile1 = 'his'
	217
	218	! if(ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5
	219	! @ .or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9) then
	220	! if (isot.eq.1) hisfile = hfile1//'26-'//ibcode1//'.dat'
	221	! if (isot.eq.2) hisfile = hfile1//'28-'//ibcode1//'.dat'
	222	! if (isot.eq.3) hisfile = hfile1//'36-'//ibcode1//'.dat'
	223	! if (isot.eq.4) hisfile = hfile1//'27-'//ibcode1//'.dat'
	224	! if (isot.eq.5) hisfile = hfile1//'62-'//ibcode1//'.dat'
	225	! else
	226	! if (isot.eq.1) hisfile = hfile1//'26-'//ibcode2//'.dat'
	227	! if (isot.eq.2) hisfile = hfile1//'28-'//ibcode2//'.dat'
	228	! if (isot.eq.3) hisfile = hfile1//'36-'//ibcode2//'.dat'
	229	! if (isot.eq.4) hisfile = hfile1//'27-'//ibcode2//'.dat'
	230	! if (isot.eq.5) hisfile = hfile1//'62-'//ibcode2//'.dat'
	231	! endif
	232	!write (,) ' /MZESCAPE/ hisfile: ', hisfile
	233
	234	! the argument to rhist is to make this compatible with mztf_comp.f,
	235	! which is a useful modification of mztf.f (to change strengths of bands
	236	! call rhist (1.0)
	237	if(ib.eq.1) then
	238	if(isot.eq.1) then !Case 1
	239	mm=mm_c1
	240	nbox=nbox_c1
	241	tmin=tmin_c1
	242	tmax=tmax_c1
	243	do i=1,nbox_max
	244	no(i)=no_c1(i)
	245	dist(i)=dist_c1(i)
	246	do j=1,nhist
	247	sk1(j,i)=sk1_c1(j,i)
	248	xls1(j,i)=xls1_c1(j,i)
	249	xln1(j,i)=xln1_c1(j,i)
	250	xld1(j,i)=xld1_c1(j,i)
	251	enddo
	252	enddo
	253	do j=1,nhist
	254	thist(j)=thist_c1(j)
	255	enddo
	256	else if(isot.eq.2) then !Case 2
	257	mm=mm_c2
	258	nbox=nbox_c2
	259	tmin=tmin_c2
	260	tmax=tmax_c2
	261	do i=1,nbox_max
	262	no(i)=no_c2(i)
	263	dist(i)=dist_c2(i)
	264	do j=1,nhist
	265	sk1(j,i)=sk1_c2(j,i)
	266	xls1(j,i)=xls1_c2(j,i)
	267	xln1(j,i)=xln1_c2(j,i)
	268	xld1(j,i)=xld1_c2(j,i)
	269	enddo
	270	enddo
	271	do j=1,nhist
	272	thist(j)=thist_c2(j)
	273	enddo
	274	else if(isot.eq.3) then !Case 3
	275	mm=mm_c3
	276	nbox=nbox_c3
	277	tmin=tmin_c3
	278	tmax=tmax_c3
	279	do i=1,nbox_max
	280	no(i)=no_c3(i)
	281	dist(i)=dist_c3(i)
	282	do j=1,nhist
	283	sk1(j,i)=sk1_c3(j,i)
	284	xls1(j,i)=xls1_c3(j,i)
	285	xln1(j,i)=xln1_c3(j,i)
	286	xld1(j,i)=xld1_c3(j,i)
	287	enddo
	288	enddo
	289	do j=1,nhist
	290	thist(j)=thist_c3(j)
	291	enddo
	292	else if(isot.eq.4) then !Case 4
	293	mm=mm_c4
	294	nbox=nbox_c4
	295	tmin=tmin_c4
	296	tmax=tmax_c4
	297	do i=1,nbox_max
	298	no(i)=no_c4(i)
	299	dist(i)=dist_c4(i)
	300	do j=1,nhist
	301	sk1(j,i)=sk1_c4(j,i)
	302	xls1(j,i)=xls1_c4(j,i)
	303	xln1(j,i)=xln1_c4(j,i)
	304	xld1(j,i)=xld1_c4(j,i)
	305	enddo
	306	enddo
	307	do j=1,nhist
	308	thist(j)=thist_c4(j)
	309	enddo
	310	else
	311	write(,)'isot must be 2,3 or 4 for ib=1!!'
	312	write(,)'stop at mzescape/312'
	313	stop
	314	endif
	315	else if (ib.eq.2) then
	316	if(isot.eq.1) then !Case 5
	317	mm=mm_c5
	318	nbox=nbox_c5
	319	tmin=tmin_c5
	320	tmax=tmax_c5
	321	do i=1,nbox_max
	322	no(i)=no_c5(i)
	323	dist(i)=dist_c5(i)
	324	do j=1,nhist
	325	sk1(j,i)=sk1_c5(j,i)
	326	xls1(j,i)=xls1_c5(j,i)
	327	xln1(j,i)=xln1_c5(j,i)
	328	xld1(j,i)=xld1_c5(j,i)
	329	enddo
	330	enddo
	331	do j=1,nhist
	332	thist(j)=thist_c5(j)
	333	enddo
	334	else
	335	write(,)'isot must be 1 for ib=2!!'
	336	write(,)'stop at mzescape/336'
	337	stop
	338	endif
	339	else if (ib.eq.3) then
	340	if(isot.eq.1) then !Case 6
	341	mm=mm_c6
	342	nbox=nbox_c6
	343	tmin=tmin_c6
	344	tmax=tmax_c6
	345	do i=1,nbox_max
	346	no(i)=no_c6(i)
	347	dist(i)=dist_c6(i)
	348	do j=1,nhist
	349	sk1(j,i)=sk1_c6(j,i)
	350	xls1(j,i)=xls1_c6(j,i)
	351	xln1(j,i)=xln1_c6(j,i)
	352	xld1(j,i)=xld1_c6(j,i)
	353	enddo
	354	enddo
	355	do j=1,nhist
	356	thist(j)=thist_c6(j)
	357	enddo
	358	else
	359	write(,)'isot must be 1 for ib=3!!'
	360	write(,)'stop at mzescape/360'
	361	stop
	362	endif
	363	else if (ib.eq.4) then
	364	if(isot.eq.1) then !Case 7
	365	mm=mm_c7
	366	nbox=nbox_c7
	367	tmin=tmin_c7
	368	tmax=tmax_c7
	369	do i=1,nbox_max
	370	no(i)=no_c7(i)
	371	dist(i)=dist_c7(i)
	372	do j=1,nhist
	373	sk1(j,i)=sk1_c7(j,i)
	374	xls1(j,i)=xls1_c7(j,i)
	375	xln1(j,i)=xln1_c7(j,i)
	376	xld1(j,i)=xld1_c7(j,i)
	377	enddo
	378	enddo
	379	do j=1,nhist
	380	thist(j)=thist_c7(j)
	381	enddo
	382	else
	383	write(,)'isot must be 1 for ib=4!!'
	384	write(,)'stop at mzescape/384'
	385	stop
	386	endif
	387	else
	388	write(,)'ib must be 1,2,3 or 4!!'
	389	write(,)'stop at mzescape/389'
	390	endif
	391	if (isot.ne.5) deltanux = deltanu(isot,ib)
	392	if (isot.eq.5) deltanux = deltanuco
	393
	394	c******
	395	c****** calculation of tauinf(nl)
	396	c******
	397	call initial
	398
	399	ff=1.0e10
	400
	401	do i=nl,1,-1
	402
	403	if(i.eq.nl)then
	404
	405	call intz (zl(i),c2,p2,mr2,t2, con)
	406	do kr=1,nbox
	407	ta(kr)=t2
	408	end do
	409	! write (,) ' i, t2 =', i, t2
	410	call interstrength (st2,t2,ka,ta)
	411	aa = p2 * coninf * mr2 * (st2 * ff)
	412	bb = p2 * coninf * st2
	413	cc = coninf * st2
	414	dd = t2 * coninf * st2
	415	do kr=1,nbox
	416	ccbox(kr) = coninf * ka(kr)
	417	ddbox(kr) = t2 * ccbox(kr)
	418	! c2box(kr) = c2 * ka(kr) * beta * dble(deltaz) * 1.d5
	419	c2box(kr) = c2 * ka(kr) * dble(deltaz)
	420	end do
	421	! c2 = c2 * st2 * beta * dble(deltaz) * 1.d5
	422	c2 = c2 * st2 * dble(deltaz)
	423
	424	else
	425	call intz (zl(i),c1,p1,mr1,t1, con)
	426	do kr=1,nbox
	427	ta(kr)=t1
	428	end do
	429	! write (,) ' i, t1 =', i, t1
	430	call interstrength (st1,t1,ka,ta)
	431	do kr=1,nbox
	432	! c1box(kr) = c1 * ka(kr) * beta * dble(deltaz) * 1.d5
	433	c1box(kr) = c1 * ka(kr) * dble(deltaz)
	434	end do
	435	! c1 = c1 * st1 * beta * dble(deltaz) * 1.d5
	436	c1 = c1 * st1 * dble(deltaz)
	437	aa = aa + ( p1mr1(c1ff) + p2mr2(c2ff)) / 2.d0
	438	bb = bb + ( p1c1 + p2c2 ) / 2.d0
	439	cc = cc + ( c1 + c2 ) / 2.d0
	440	ccc = ( c1 + c2 ) / 2.d0
	441	dd = dd + ( t1c1 + t2c2 ) / 2.d0
	442	do kr=1,nbox
	443	ccbox(kr) = ccbox(kr) +
	444	@ ( c1box(kr) + c2box(kr) )/2.d0
	445	ddbox(kr) = ddbox(kr) +
	446	@ ( t1c1box(kr)+t2c2box(kr) )/2.d0
	447	end do
	448
	449	mr2 = mr1
	450	c2=c1
	451	do kr=1,nbox
	452	c2box(kr) = c1box(kr)
	453	end do
	454	t2=t1
	455	p2=p1
	456	end if
	457
	458	pt = bb / cc
	459	pp = aa / (cc*ff)
	460
	461	! ta=dd/cc
	462	! tdop = ta
	463	ts = dd/cc
	464	do kr=1,nbox
	465	ta(kr) = ddbox(kr) / ccbox(kr)
	466	end do
	467	! write (,) ' i, ts =', i, ts
	468	call interstrength(st,ts,ka,ta)
	469	! call intershape(alsa,alna,alda,tdop)
	470	call intershape(alsa,alna,alda,ta)
	471
	472	* ua = cc/st
	473
	474	c next loop calculates the eqw for an especified path ua,pp,pt,ta
	475
	476	eqwmu = 0.0d0
	477	do im = 1,iimu
	478	eqw=0.0d0
	479	do kr=1,nbox
	480	ua(kr) = ccbox(kr) / ka(kr) * beta * 1.0d5 / mu(im)
	481	if(ua(kr).lt.0.)write(,)'mzescape/480',ua(kr),
	482	$ ccbox(kr),ka(kr),beta,mu(im),kr,im,i,nl
	483
	484	call findw (ig,iirw, 0, csL,psL, Desp, wsL)
	485	if ( i_supersat .eq. 0 ) then
	486	eqw=eqw+no(kr)*w
	487	else
	488	eqw = w + (no(kr)-1) * ( asat_box*dist(kr) )
	489	endif
	490	end do
	491	eqwmu = eqwmu + eqw * mu(im)*amu(im)
	492	end do
	493
	494	! tauinf(i) = exp( - eqwmu / dble(deltanux) )
	495	tauinf(i) = 1.d0 - eqwmu / dble(deltanux)
	496	if (tauinf(i).lt.0.d0) tauinf(i) = 0.0d0
	497
	498	if (i.eq.nl) then
	499	taustar(i) = 0.0d0
	500	else
	501	taustar(i) = dble(deltanux) * (tauinf(i+1)-tauinf(i))
	502	! ~ / ( beta * cc * 1.d5 )
	503	~ / ( beta * ccc * 1.d5 )
	504	endif
	505
	506	end do ! i continue
	507
	508
	509	c******
	510	c****** calculation of tau(in,ir) for n<=r
	511	c******
	512
	513	do 1 in=1,nl-1
	514
	515	call initial
	516
	517	call intz (zl(in), c1,p1,mr1,t1, con)
	518	do kr=1,nbox
	519	ta(kr) = t1
	520	end do
	521	call interstrength (st1,t1,ka,ta)
	522	do kr=1,nbox
	523	c1box(kr) = c1 * ka(kr) * dble(deltaz)
	524	end do
	525	c1 = c1 * st1 * dble(deltaz)
	526
	527	call intz (zl(in+1), c2,p2,mr2,t2, con)
	528	do kr=1,nbox
	529	ta(kr) = t2
	530	end do
	531	call interstrength (st2,t2,ka,ta)
	532	do kr=1,nbox
	533	c2box(kr) = c2 * ka(kr) * dble(deltaz)
	534	end do
	535	c2 = c2 * st2 * dble(deltaz)
	536
	537	aa = aa + ( p1mr1(c1ff) + p2mr2(c2ff)) / 2.d0
	538	bb = bb + ( p1c1 + p2c2 ) / 2.d0
	539	cc = cc + ( c1 + c2 ) / 2.d0
	540	dd = dd + ( t1c1 + t2c2 ) / 2.d0
	541	do kr=1,nbox
	542	ccbox(kr) = ccbox(kr) + (c1box(kr)+c2box(kr))/2.d0
	543	ddbox(kr) = ddbox(kr) + (t1c1box(kr)+t2c2box(kr))/2.d0
	544	end do
	545
	546	mr1=mr2
	547	t1=t2
	548	c1=c2
	549	p1=p2
	550	do kr=1,nbox
	551	c1box(kr) = c2box(kr)
	552	end do
	553	pt = bb / cc
	554	pp = aa / (cc * ff)
	555	ts = dd/cc
	556	do kr=1,nbox
	557	ta(kr) = ddbox(kr) / ccbox(kr)
	558	end do
	559	call interstrength(st,ts,ka,ta)
	560	call intershape(alsa,alna,alda,ta)
	561
	562	eqwmu = 0.0d0
	563	do im = 1,iimu
	564	eqw=0.0d0
	565	do kr=1,nbox
	566	ua(kr) = ccbox(kr) / ka(kr) * beta * 1.0d5 / mu(im)
	567	if(ua(kr).lt.0.)write(,)'mzescape/566',ua(kr),
	568	$ ccbox(kr),ka(kr),beta,mu(im),kr,im,i,nl
	569
	570	call findw (ig,iirw, 0, csL,psL, Desp, wsL)
	571	if ( i_supersat .eq. 0 ) then
	572	eqw=eqw+no(kr)*w
	573	else
	574	eqw = w + (no(kr)-1) * ( asat_box*dist(kr) )
	575	endif
	576	end do
	577	eqwmu = eqwmu + eqw * mu(im)*amu(im)
	578	end do
	579
	580	tauii(in) = exp( - eqwmu / dble(deltanux) )
	581	!write (,) 'i,tauii=',in,tauii(in)
	582
	583	1 continue
	584	tauii(nl) = 1.0d0
	585
	586
	587	c end
	588	return
	589	end
	590
	591
	592
	593
	594
	595
	596
	597
	598
	599

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

Download in other formats:

Original Format