Context Navigation

mztf_overlap.F @ 469

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

Line
1	c***********************************************************************
2	c mztf.f
3	c***********************************************************************
4	c
5	c program for calculating atmospheric transmittances
6	c to be used in the calculation of curtis matrix coefficients
7
8	subroutine mztf ( ig,cf,cfup,cfdw,vc,taugr, ib,isot,
9	@ iirw,iimu,itauout,icfout,itableout )
10
11	c i*out = 1 output of data
12	c i*out = 0 no output
13	c
14	c jul 2011 malv+fgg adapted to LMD-MGCM
15	c nov 98 mavl allow for overlaping in the lorentz line
16	c jan 98 malv version for mz1d. based on curtis/mztf.for
17	c 17-jul-96 mlp&crs change the calculation of mr.
18	c evitar: divide por cero. anhadiendo: ff
19	c oct-92 malv correct s(t) dependence for all histogr bands
20	c june-92 malv proper lower levels for laser bands
21	c may-92 malv new temperature dependence for laser bands
22	c @ 991 malv boxing for the averaged absorber amount and t
23	c ? malv extension up to 200 km altitude in mars
24	c 13-nov-86 mlp include the temperature weighted to match
25	c the eqw in the strong doppler limit.
26	c***********************************************************************
27
28	implicit none
29
30	include 'nltedefs.h'
31	include 'nlte_atm.h'
32	include 'nlte_data.h'
33	include 'nlte_curtis.h'
34	include 'tcr_15um.h'
35	include 'nlte_results.h'
36
37
38	c arguments
39	integer ig !ADDED FOR TRACEBACK
40	real*8 cf(nl,nl), cfup(nl,nl), cfdw(nl,nl) ! o.
41	real*8 vc(nl), taugr(nl) ! o
42	integer ib ! i
43	integer isot ! i
44	integer iirw ! i
45	integer iimu ! i
46	integer itauout ! i
47	integer icfout ! i
48	integer itableout ! i
49
50	c local variables and constants
51	integer i, in, ir, im, k ,j
52	integer nmu
53	parameter (nmu = 8)
54	real*8 tau(nl,nl)
55	real*8 tauinf(nl)
56	real*8 con(nzy), coninf
57	real*8 c1, c2
58	real*8 t1, t2
59	real*8 p1, p2
60	real*8 mr1, mr2
61	real*8 st1, st2
62	real*8 c1box(70), c2box(70)
63	real*8 ff ! to avoid too small numbers
64	real*8 tvtbs(nzy)
65	real*8 st, beta, ts, eqwmu
66	real*8 mu(nmu), amu(nmu)
67	real*8 zld(nl), zyd(nzy)
68	real*8 correc
69	real deltanux ! width of vib-rot band (cm-1)
70	! character isotcode*2
71	integer idummy
72	real*8 Desp,wsL
73
74	c formats
75	! 111 format(a1)
76	! 112 format(a2)
77	101 format(i1)
78	202 format(i2)
79	! 180 format(a80)
80	! 181 format(a80)
81	c***********************************************************************
82
83	c some needed values
84	! rl=sqrt(log(2.d0))
85	! pi2 = 3.14159265358989d0
86	beta = 1.8d0
87	idummy = 0
88	Desp = 0.d0
89	wsL = 0.d0
90
91	c esto es para que las subroutines de mztfsub calculen we
92	c de la forma apropiada para mztf, no para fot
93	icls=icls_mztf
94
95	c codigos para filenames
96	! if (isot .eq. 1) isotcode = '26'
97	! if (isot .eq. 2) isotcode = '28'
98	! if (isot .eq. 3) isotcode = '36'
99	! if (isot .eq. 4) isotcode = '27'
100	! if (isot .eq. 5) isotcode = '62'
101	! if(ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5
102	! @ .or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9) then
103	! write (ibcode1,101) ib
104	! else
105	! write (ibcode2,202) ib
106	! endif
107	! write (*,'( 30h calculating curtis matrix : ,2x,
108	! @ 8h band = ,i2,2x, 11h isotope = ,i2)') ib, isot
109
110	c integration in angle !!!!!!!!!!!!!!!!!!!!
111
112	c------- diffusivity approx.
113	if (iimu.eq.1) then
114	! write (,) ' diffusivity approx. beta = ',beta
115	mu(1) = 1.0d0
116	amu(1)= 1.0d0
117	c-------data for 8 points integration
118	elseif (iimu.eq.4) then
119	write (,)' 4 points for the gauss-legendre angle quadrature.'
120	mu(1)=(1.0d0+0.339981043584856)/2.0d0
121	mu(2)=(1.0d0-0.339981043584856)/2.0d0
122	mu(3)=(1.0d0+0.861136311594053)/2.0d0
123	mu(4)=(1.0d0-0.861136311594053)/2.0d0
124	amu(1)=0.652145154862546
125	amu(2)=amu(1)
126	amu(3)=0.347854845137454
127	amu(4)=amu(3)
128	beta=1.0d0
129	c-------data for 8 points integration
130	elseif(iimu.eq.8) then
131	write (,)' 8 points for the gauss-legendre angle quadrature.'
132	mu(1)=(1.0d0+0.183434642495650)/2.0d0
133	mu(2)=(1.0d0-0.183434642495650)/2.0d0
134	mu(3)=(1.0d0+0.525532409916329)/2.0d0
135	mu(4)=(1.0d0-0.525532409916329)/2.0d0
136	mu(5)=(1.0d0+0.796666477413627)/2.0d0
137	mu(6)=(1.0d0-0.796666477413627)/2.0d0
138	mu(7)=(1.0d0+0.960289856497536)/2.0d0
139	mu(8)=(1.0d0-0.960289856497536)/2.0d0
140	amu(1)=0.362683783378362
141	amu(2)=amu(1)
142	amu(3)=0.313706645877887
143	amu(4)=amu(3)
144	amu(5)=0.222381034453374
145	amu(6)=amu(5)
146	amu(7)=0.101228536290376
147	amu(8)=amu(7)
148	beta=1.0d0
149	end if
150	c!!!!!!!!!!!!!!!!!!!!!!!
151
152	ccc
153	ccc determine abundances included in the absorber amount
154	ccc
155
156	c first, set up the grid ready for interpolation.
157	do i=1,nzy
158	zyd(i) = dble(zy(i))
159	enddo
160	do i=1,nl
161	zld(i) = dble(zl(i))
162	enddo
163
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	c vibr. temp of the bending mode :
174	if (isot.eq.1) call interdp ( tvtbs,zyd,nzy, v626t1,zld,nl, 1 )
175	if (isot.eq.2) call interdp ( tvtbs,zyd,nzy, v628t1,zld,nl, 1 )
176	if (isot.eq.3) call interdp ( tvtbs,zyd,nzy, v636t1,zld,nl, 1 )
177	if (isot.eq.4) call interdp ( tvtbs,zyd,nzy, v627t1,zld,nl, 1 )
178	correc = 2.d0 * dexp( dble(-ee*elow(isot,2))/tvtbs(i) )
179	con(i) = con(i) * ( 1.d0 - correc )
180	endif
181	c-----------------------------------------------------------------------
182	c mlp & cristina. 17 july 1996
183	c change the calculation of mr. it is used for calculating partial press
184	c alpha = alpha(self,co2)pp +alpha(n2)(pt-pp)
185	c for an isotope, if mr is obtained by co2*imr(iso)/nt we are considerin
186	c collisions with other co2 isotopes (including the major one, 626)
187	c as if they were with n2. assuming mr as co2/nt, we consider collisions
188	c of type 628-626 as of 626-626 instead of as 626-n2.
189	c mrx(i)=con(i)/ntx(i) ! old malv
190
191	! mrx(i)= dble(co2x(i)/ntx(i)) ! mlp & crs
192
193	c jan 98:
194	c esta modif de mlp implica anular el correc (deberia revisar esto)
195	mr(i) = dble(co2y(i)/nty(i)) ! malv, jan 98
196
197	c-----------------------------------------------------------------------
198
199	end do
200
201	! como beta y 1.d5 son comunes a todas las weighted absorber amounts,
202	! los simplificamos:
203	! coninf = beta * 1.d5 * dble( con(n) / log( con(n-1) / con(n) ) )
204	coninf = dble( con(nzy) / log( con(nzy-1) / con(nzy) ) )
205
206	! write (,) ' coninf =', coninf
207
208	ccc
209	ccc temp dependence of the band strength and
210	ccc nlte correction factor for the absorber amount
211	ccc
212	call mztf_correccion ( coninf, con, ib, isot, itableout )
213
214	ccc
215	ccc reads histogrammed spectral data (strength for lte and vmr=1)
216	ccc
217	!hfile1 = dirspec//'hi'//dn ! ya no distinguimos entre d/n
218	!! hfile1 = dirspec//'hid' ! (see why in his.for)
219	! hfile='hid'
220	!! if (ib.eq.13 .or. ib.eq.14 ) hfile1 = dirspec//'his'
221	! if (ib.eq.13 .or. ib.eq.14 ) hfile1 = 'his'
222	!
223	! if(ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5
224	! @ .or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9) then
225	! if (isot.eq.1) hisfile = hfile1//'26-'//ibcode1//'.dat'
226	! if (isot.eq.2) hisfile = hfile1//'28-'//ibcode1//'.dat'
227	! if (isot.eq.3) hisfile = hfile1//'36-'//ibcode1//'.dat'
228	! if (isot.eq.4) hisfile = hfile1//'27-'//ibcode1//'.dat'
229	! if (isot.eq.5) hisfile = hfile1//'62-'//ibcode1//'.dat'
230	! else
231	! if (isot.eq.1) hisfile = hfile1//'26-'//ibcode2//'.dat'
232	! if (isot.eq.2) hisfile = hfile1//'28-'//ibcode2//'.dat'
233	! if (isot.eq.3) hisfile = hfile1//'36-'//ibcode2//'.dat'
234	! if (isot.eq.4) hisfile = hfile1//'27-'//ibcode2//'.dat'
235	! if (isot.eq.5) hisfile = hfile1//'62-'//ibcode2//'.dat'
236	! endif
237	! write (,) 'hisfile: ', hisfile
238
239	! the argument to rhist is to make this compatible with mztf_comp.f,
240	! which is a useful modification of mztf.f (to change strengths of bands
241	! call rhist (1.0)
242	if(ib.eq.1) then
243	if(isot.eq.1) then !Case 1
244	mm=mm_c1
245	nbox=nbox_c1
246	tmin=tmin_c1
247	tmax=tmax_c1
248	do i=1,nbox_max
249	no(i)=no_c1(i)
250	dist(i)=dist_c1(i)
251	do j=1,nhist
252	sk1(j,i)=sk1_c1(j,i)
253	xls1(j,i)=xls1_c1(j,i)
254	xln1(j,i)=xln1_c1(j,i)
255	xld1(j,i)=xld1_c1(j,i)
256	enddo
257	enddo
258	do j=1,nhist
259	thist(j)=thist_c1(j)
260	enddo
261	else if(isot.eq.2) then !Case 2
262	mm=mm_c2
263	nbox=nbox_c2
264	tmin=tmin_c2
265	tmax=tmax_c2
266	do i=1,nbox_max
267	no(i)=no_c2(i)
268	dist(i)=dist_c2(i)
269	do j=1,nhist
270	sk1(j,i)=sk1_c2(j,i)
271	xls1(j,i)=xls1_c2(j,i)
272	xln1(j,i)=xln1_c2(j,i)
273	xld1(j,i)=xld1_c2(j,i)
274	enddo
275	enddo
276	do j=1,nhist
277	thist(j)=thist_c2(j)
278	enddo
279	else if(isot.eq.3) then !Case 3
280	mm=mm_c3
281	nbox=nbox_c3
282	tmin=tmin_c3
283	tmax=tmax_c3
284	do i=1,nbox_max
285	no(i)=no_c3(i)
286	dist(i)=dist_c3(i)
287	do j=1,nhist
288	sk1(j,i)=sk1_c3(j,i)
289	xls1(j,i)=xls1_c3(j,i)
290	xln1(j,i)=xln1_c3(j,i)
291	xld1(j,i)=xld1_c3(j,i)
292	enddo
293	enddo
294	do j=1,nhist
295	thist(j)=thist_c3(j)
296	enddo
297	else if(isot.eq.4) then !Case 4
298	mm=mm_c4
299	nbox=nbox_c4
300	tmin=tmin_c4
301	tmax=tmax_c4
302	do i=1,nbox_max
303	no(i)=no_c4(i)
304	dist(i)=dist_c4(i)
305	do j=1,nhist
306	sk1(j,i)=sk1_c4(j,i)
307	xls1(j,i)=xls1_c4(j,i)
308	xln1(j,i)=xln1_c4(j,i)
309	xld1(j,i)=xld1_c4(j,i)
310	enddo
311	enddo
312	do j=1,nhist
313	thist(j)=thist_c4(j)
314	enddo
315	else
316	write(,)'isot must be 2,3 or 4 for ib=1!!'
317	write(,)'stop at mztf_overlap/317'
318	stop
319	endif
320	else if (ib.eq.2) then
321	if(isot.eq.1) then !Case 5
322	mm=mm_c5
323	nbox=nbox_c5
324	tmin=tmin_c5
325	tmax=tmax_c5
326	do i=1,nbox_max
327	no(i)=no_c5(i)
328	dist(i)=dist_c5(i)
329	do j=1,nhist
330	sk1(j,i)=sk1_c5(j,i)
331	xls1(j,i)=xls1_c5(j,i)
332	xln1(j,i)=xln1_c5(j,i)
333	xld1(j,i)=xld1_c5(j,i)
334	enddo
335	enddo
336	do j=1,nhist
337	thist(j)=thist_c5(j)
338	enddo
339	else
340	write(,)'isot must be 1 for ib=2!!'
341	write(,)'stop at mztf_overlap/341'
342	stop
343	endif
344	else if (ib.eq.3) then
345	if(isot.eq.1) then !Case 6
346	mm=mm_c6
347	nbox=nbox_c6
348	tmin=tmin_c6
349	tmax=tmax_c6
350	do i=1,nbox_max
351	no(i)=no_c6(i)
352	dist(i)=dist_c6(i)
353	do j=1,nhist
354	sk1(j,i)=sk1_c6(j,i)
355	xls1(j,i)=xls1_c6(j,i)
356	xln1(j,i)=xln1_c6(j,i)
357	xld1(j,i)=xld1_c6(j,i)
358	enddo
359	enddo
360	do j=1,nhist
361	thist(j)=thist_c6(j)
362	enddo
363	else
364	write(,)'isot must be 1 for ib=3!!'
365	write(,)'stop at mztf_overlap/365'
366	stop
367	endif
368	else if (ib.eq.4) then
369	if(isot.eq.1) then !Case 7
370	mm=mm_c7
371	nbox=nbox_c7
372	tmin=tmin_c7
373	tmax=tmax_c7
374	do i=1,nbox_max
375	no(i)=no_c7(i)
376	dist(i)=dist_c7(i)
377	do j=1,nhist
378	sk1(j,i)=sk1_c7(j,i)
379	xls1(j,i)=xls1_c7(j,i)
380	xln1(j,i)=xln1_c7(j,i)
381	xld1(j,i)=xld1_c7(j,i)
382	enddo
383	enddo
384	do j=1,nhist
385	thist(j)=thist_c7(j)
386	enddo
387	else
388	write(,)'isot must be 1 for ib=4!!'
389	write(,)'stop at mztf_overlap/389'
390	stop
391	endif
392	else
393	write(,)'ib must be 1,2,3 or 4!!'
394	write(,)'stop at mztf_overlap/394'
395	endif
396
397	if (isot.ne.5) deltanux = deltanu(isot,ib)
398	if (isot.eq.5) deltanux = deltanuco
399
400	c******
401	c****** calculation of tauinf(nl)
402	c******
403	call initial
404
405	ff=1.0e10
406
407	do i=nl,1,-1
408
409	if(i.eq.nl)then
410
411	call intz (zl(i),c2,p2,mr2,t2, con)
412	do kr=1,nbox
413	ta(kr)=t2
414	end do
415	! write (,) ' i, t2 =', i, t2
416	call interstrength (st2,t2,ka,ta)
417	aa = p2 * coninf * mr2 * (st2 * ff)
418	bb = p2 * coninf * st2
419	cc = coninf * st2
420	dd = t2 * coninf * st2
421	do kr=1,nbox
422	ccbox(kr) = coninf * ka(kr)
423	ddbox(kr) = t2 * ccbox(kr)
424	! c2box(kr) = c2 * ka(kr) * beta * dble(deltaz) * 1.d5
425	c2box(kr) = c2 * ka(kr) * dble(deltaz)
426	end do
427	! c2 = c2 * st2 * beta * dble(deltaz) * 1.d5
428	c2 = c2 * st2 * dble(deltaz)
429
430	else
431	call intz (zl(i),c1,p1,mr1,t1, con)
432	do kr=1,nbox
433	ta(kr)=t1
434	end do
435	! write (,) ' i, t1 =', i, t1
436	call interstrength (st1,t1,ka,ta)
437	do kr=1,nbox
438	! c1box(kr) = c1 * ka(kr) * beta * dble(deltaz) * 1.d5
439	c1box(kr) = c1 * ka(kr) * dble(deltaz)
440	end do
441	! c1 = c1 * st1 * beta * dble(deltaz) * 1.d5
442	c1 = c1 * st1 * dble(deltaz)
443	aa = aa + ( p1mr1(c1ff) + p2mr2(c2ff)) / 2.d0
444	bb = bb + ( p1c1 + p2c2 ) / 2.d0
445	cc = cc + ( c1 + c2 ) / 2.d0
446	dd = dd + ( t1c1 + t2c2 ) / 2.d0
447	do kr=1,nbox
448	ccbox(kr) = ccbox(kr) + ( c1box(kr) + c2box(kr) )/2.d0
449	ddbox(kr) = ddbox(kr) + ( t1c1box(kr)+t2c2box(kr) )/2.d0
450	end do
451
452	mr2 = mr1
453	c2=c1
454	do kr=1,nbox
455	c2box(kr) = c1box(kr)
456	end do
457	t2=t1
458	p2=p1
459	end if
460
461	pt = bb / cc
462	pp = aa / (cc*ff)
463
464	! ta=dd/cc
465	! tdop = ta
466	ts = dd/cc
467	do kr=1,nbox
468	ta(kr) = ddbox(kr) / ccbox(kr)
469	end do
470	! write (,) ' i, ts =', i, ts
471	call interstrength(st,ts,ka,ta)
472	! call intershape(alsa,alna,alda,tdop)
473	call intershape(alsa,alna,alda,ta)
474
475	* ua = cc/st
476
477	c next loop calculates the eqw for an especified path ua,pp,pt,ta
478
479	eqwmu = 0.0d0
480	do im = 1,iimu
481	eqw=0.0d0
482	do kr=1,nbox
483	ua(kr) = ccbox(kr) / ka(kr) * beta * 1.0d5 / mu(im)
484	if(ua(kr).lt.0.)write(,)'mztf_overlap/483',ua(kr),
485	$ ccbox(kr),ka(kr),beta,mu(im),kr,im,i,nl
486
487	call findw(ig,iirw, idummy,c1,p1,Desp,wsL)
488	if ( i_supersat .eq. 0 ) then
489	eqw=eqw+no(kr)*w
490	else
491	eqw = w + (no(kr)-1) * ( asat_box*dist(kr) )
492	endif
493	end do
494	eqwmu = eqwmu + eqw * mu(im)*amu(im)
495	end do
496
497	tauinf(i) = exp( - eqwmu / dble(deltanux) )
498
499	end do ! i continue
500
501	! if ( isot.eq.1 .and. ib.eq.2 ) then
502	! write (,) ' tauinf(nl) = ', tauinf(nl)
503	! write (,) ' tauinf(1) = ', tauinf(1)
504	! endif
505
506	c******
507	c****** calculation of tau(in,ir) for n<=r
508	c******
509
510	do 1 in=1,nl-1
511
512	call initial
513	call intz (zl(in), c1,p1,mr1,t1, con)
514	do kr=1,nbox
515	ta(kr) = t1
516	end do
517	call interstrength (st1,t1,ka,ta)
518	do kr=1,nbox
519	! c1box(kr) = c1 * ka(kr) * beta * dble(deltaz) * 1.d5
520	c1box(kr) = c1 * ka(kr) * dble(deltaz)
521	end do
522	! c1 = c1 * st1 * beta * dble(deltaz) * 1.d5
523	c1 = c1 * st1 * dble(deltaz)
524
525	do 2 ir=in,nl-1
526
527	if (ir.eq.in) then
528	tau(in,ir) = 1.d0
529	goto 2
530	end if
531
532	call intz (zl(ir), c2,p2,mr2,t2, con)
533	do kr=1,nbox
534	ta(kr) = t2
535	end do
536	call interstrength (st2,t2,ka,ta)
537	do kr=1,nbox
538	! c2box(kr) = c2 * ka(kr) * beta * dble(deltaz) * 1.d5
539	c2box(kr) = c2 * ka(kr) * dble(deltaz)
540	end do
541	! c2 = c2 * st2 * beta * dble(deltaz) * 1.e5
542	c2 = c2 * st2 * dble(deltaz)
543
544	c aa = aa + ( p1mr1c1 + p2mr2c2 ) / 2.d0
545	aa = aa + ( p1mr1(c1ff) + p2mr2(c2ff)) / 2.d0
546	bb = bb + ( p1c1 + p2c2 ) / 2.d0
547	cc = cc + ( c1 + c2 ) / 2.d0
548	dd = dd + ( t1c1 + t2c2 ) / 2.d0
549	do kr=1,nbox
550	ccbox(kr) = ccbox(kr) + ( c1box(kr) + c2box(kr) ) /2.d0
551	ddbox(kr) = ddbox(kr) + ( t1c1box(kr) + t2c2box(kr) ) /2.d0
552	end do
553
554	mr1=mr2
555	t1=t2
556	c1=c2
557	p1=p2
558	do kr=1,nbox
559	c1box(kr) = c2box(kr)
560	end do
561
562	pt = bb / cc
563	pp = aa / (cc * ff)
564
565	* ta=dd/cc
566	* tdop = ta
567	ts = dd/cc
568	do kr=1,nbox
569	ta(kr) = ddbox(kr) / ccbox(kr)
570	end do
571	call interstrength(st,ts,ka,ta)
572	call intershape(alsa,alna,alda,ta)
573	* ua = cc/st
574
575	c next loop calculates the eqw for an especified path ua,pp,pt,ta
576
577	eqwmu = 0.0d0
578	do im = 1,iimu
579	eqw=0.0d0
580	do kr=1,nbox
581	ua(kr) = ccbox(kr) / ka(kr) * beta * 1.0d5 / mu(im)
582	if(ua(kr).lt.0.)write(,)'mztf_overlap/581',ua(kr),
583	$ ccbox(kr),ka(kr),beta,mu(im),kr,im,i,nl
584
585	call findw(ig,iirw, idummy,c1,p1,Desp,wsL)
586	if ( i_supersat .eq. 0 ) then
587	eqw=eqw+no(kr)*w
588	else
589	eqw = w + (no(kr)-1) * ( asat_box*dist(kr) )
590	endif
591	end do
592	eqwmu = eqwmu + eqw * mu(im)*amu(im)
593	end do
594
595	tau(in,ir) = exp( - eqwmu / dble(deltanux) )
596
597	2 continue
598
599	1 continue
600
601	! if ( isot.eq.1 .and. ib.eq.2 ) then
602	! write (,) ' tau(1,) , =1,20 '
603	! write (,) ( sngl(tau(1,k)), k=1,20 )
604	! endif
605
606
607	c**********
608	c********** calculation of tau(in,ir) for n>r
609	c**********
610
611	in=nl
612
613	call initial
614	call intz (zl(in), c1,p1,mr1,t1, con)
615	do kr=1,nbox
616	ta(kr) = t1
617	end do
618	call interstrength (st1,t1,ka,ta)
619	do kr=1,nbox
620	! c1box(kr) = c1 * ka(kr) * beta * dble(deltaz) * 1.d5
621	c1box(kr) = c1 * ka(kr) * dble(deltaz)
622	end do
623	! c1 = c1 * st1 * beta * dble(deltaz) * 1.d5
624	c1 = c1 * st1 * dble(deltaz)
625
626	do 4 ir=in-1,1,-1
627
628	call intz (zl(ir), c2,p2,mr2,t2, con)
629	do kr=1,nbox
630	ta(kr) = t2
631	end do
632	call interstrength (st2,t2,ka,ta)
633	do kr=1,nbox
634	! c2box(kr) = c2 * ka(kr) * beta * dble(deltaz) * 1.d5
635	c2box(kr) = c2 * ka(kr) * dble(deltaz)
636	end do
637	! c2 = c2 * st2 * beta * dble(deltaz) * 1.d5
638	c2 = c2 * st2 * dble(deltaz)
639
640	aa = aa + ( p1mr1(c1ff) + p2mr2(c2ff)) / 2.d0
641	bb = bb + ( p1c1 + p2c2 ) / 2.d0
642	cc = cc + ( c1 + c2 ) / 2.d0
643	dd = dd + ( t1c1 + t2c2 ) / 2.d0
644	do kr=1,nbox
645	ccbox(kr) = ccbox(kr) + ( c1box(kr) + c2box(kr) ) /2.d0
646	ddbox(kr) = ddbox(kr) + ( t1c1box(kr) + t2c2box(kr) ) /2.d0
647	end do
648
649	mr1=mr2
650	c1=c2
651	t1=t2
652	p1=p2
653	do kr=1,nbox
654	c1box(kr) = c2box(kr)
655	end do
656
657	pt = bb / cc
658	pp = aa / (cc * ff)
659	ts = dd / cc
660	do kr=1,nbox
661	ta(kr) = ddbox(kr) / ccbox(kr)
662	end do
663	call interstrength (st,ts,ka,ta)
664	call intershape (alsa,alna,alda,ta)
665
666	* ua = cc/st
667
668	c next loop calculates the eqw for an especified path ua,pp,pt,ta
669
670	eqwmu = 0.0d0
671	do im = 1,iimu
672	eqw=0.0d0
673	do kr=1,nbox
674	ua(kr) = ccbox(kr) / ka(kr) * beta * 1.0d5 / mu(im)
675	if(ua(kr).lt.0.)write(,)'mztf_overlap/674',ua(kr),
676	$ ccbox(kr),ka(kr),beta,mu(im),kr,im,i,nl
677
678	call findw(ig,iirw, idummy,c1,p1,Desp,wsL)
679	if ( i_supersat .eq. 0 ) then
680	eqw=eqw+no(kr)*w
681	else
682	eqw = w + (no(kr)-1) * ( asat_box*dist(kr) )
683	endif
684	end do
685	eqwmu = eqwmu + eqw * mu(im)*amu(im)
686	end do
687
688	tau(in,ir) = exp( - eqwmu / dble(deltanux) )
689
690	4 continue
691
692	c
693	c due to the simmetry of the transmittances
694	c
695	do in=nl-1,2,-1
696	do ir=in-1,1,-1
697	tau(in,ir) = tau(ir,in)
698	end do
699	end do
700
701
702	ccc
703	ccc writing out transmittances
704	ccc
705	if (itauout.eq.1) then
706
707	! if(ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5
708	! @ .or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9) then
709	! open( 1, file=
710	! @ dircurtis//'taul'//isotcode//dn//ibcode1//'.dat',
711	! @ access='sequential', form='unformatted' )
712	! else
713	! open( 1, file=
714	! @ dircurtis//'taul'//isotcode//dn//ibcode2//'.dat',
715	! @ access='sequential', form='unformatted' )
716	! endif
717
718	! write(1) dummy
719	! write(1)' format: (tauinf(n),(tau(n,r),r=1,nl),n=1,nl)'
720	! do in=1,nl
721	! write (1) tauinf(in), ( tau(in,ir), ir=1,nl )
722	! end do
723	! close(unit=1)
724
725	elseif (itauout.eq.2) then
726
727	! if(ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5
728	! @ .or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9) then
729	! open( 1, file=
730	! @ dircurtis//'taul'//isotcode//dn//ibcode1//'.dat')
731	! else
732	! open( 1, file=
733	! @ dircurtis//'taul'//isotcode//dn//ibcode2//'.dat')
734	! endif
735
736	! !write(1,*) dummy
737	! !write(1,*) 'tij for curtis matrix calculations '
738	! !write(1,*)' cira mars model atmosphere '
739	! write(1,*)' beta= ',beta,'deltanu= ',deltanux
740	! write(1,*)' number of elements (in,ir)= ',nl,nl
741	! write(1,*)' format: (tauinf(in),(tau(in,ir),ir=1,nl),in=1,nl)'
742
743	! do in=1,nl
744	! write (1,*) tauinf(in)
745	! do ir=1,nl
746	! write(1,*) tau(in,ir)
747	! end do
748	! end do
749	! close(unit=1)
750
751	! if(ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5
752	! @ .or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9) then
753	! write (*,'(1x, 31htransmitances written out in: ,a22)')
754	! @ 'taul'//isotcode//dn//ibcode1
755	! else
756	! write (*,'(1x, 31htransmitances written out in: ,a22)')
757	! @ 'taul'//isotcode//dn//ibcode2
758	! endif
759
760	end if
761
762	c cleaning of transmittances
763	! call elimin_tau(tau,tauinf,nl,nan,itableout,nw,dummy,
764	! @ isotcode,dn,ibcode2)
765
766	c construction of the curtis matrix
767
768	call mzcf ( tauinf,tau, cf,cfup,cfdw, vc,taugr,
769	@ ib,isot,icfout,itableout )
770
771
772	c end
773	return
774	end

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

Download in other formats:

Original Format