Context Navigation

mzescape.F @ 426

Last change on this file since 426 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

Line
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