Context Navigation

nlte_calc.F @ 501

Last change on this file since 501 was 498, checked in by emillour, 13 years ago
Mars GCM: Cleanup of the NLTE routines which have been packed together to limit the number of files. Also enforced that file names are non-capitalized (needed by the create_make_gcm scripts to better evaluate dependencies when building the makefile). FGG+EM
File size: 172.8 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 'nlte_paramdef.h'
28	include 'nlte_commons.h'
29
30
31	c arguments
32	integer ig ! ADDED FOR TRACEBACK
33	real*8 taustar(nl) ! o
34	real*8 tauinf(nl) ! o
35	real*8 tauii(nl) ! o
36	integer ib ! i
37	integer isot ! i
38	integer iirw ! i
39	integer iimu ! i
40
41
42	c local variables and constants
43	integer i, in, ir, im, k,j
44	integer nmu
45	parameter (nmu = 8)
46	! real*8 tauinf(nl)
47	real*8 con(nzy), coninf
48	real*8 c1, c2, ccc
49	real*8 t1, t2
50	real*8 p1, p2
51	real*8 mr1, mr2
52	real*8 st1, st2
53	real*8 c1box(70), c2box(70)
54	real*8 ff ! to avoid too small numbers
55	real*8 tvtbs(nzy)
56	real*8 st, beta, ts, eqwmu
57	real*8 mu(nmu), amu(nmu)
58	real*8 zld(nl), zyd(nzy)
59	real*8 correc
60	real deltanux ! width of vib-rot band (cm-1)
61	character isotcode*2
62	real*8 maximum
63	real*8 csL, psL, Desp, wsL ! for Strong Lorentz limit
64
65	c formats
66	111 format(a1)
67	112 format(a2)
68	101 format(i1)
69	202 format(i2)
70	180 format(a80)
71	181 format(a80)
72	c***********************************************************************
73
74	c some needed values
75	! rl=sqrt(log(2.d0))
76	! pi2 = 3.14159265358989d0
77	beta = 1.8d0
78	! imrco = 0.9865
79
80	c esto es para que las subroutines de mztfsub calculen we
81	c de la forma apropiada para mztf, no para fot
82	icls=icls_mztf
83
84	c codigos para filenames
85	! if (isot .eq. 1) isotcode = '26'
86	! if (isot .eq. 2) isotcode = '28'
87	! if (isot .eq. 3) isotcode = '36'
88	! if (isot .eq. 4) isotcode = '27'
89	! if (isot .eq. 5) isotcode = '62'
90	! if(ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5
91	! @ .or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9) then
92	!! encode(2,101,ibcode1)ib
93	! write ( ibcode1, 101) ib
94	! else
95	!! encode(2,202,ibcode2)ib
96	! write (ibcode2, 202) ib
97	! endif
98	! write (*,'( 30h calculating curtis matrix : ,2x,
99	! @ 8h band = ,i2,2x, 11h isotope = ,i2)') ib, isot
100
101	c integration in angle !!!!!!!!!!!!!!!!!!!!
102
103	c------- diffusivity approx.
104	if (iimu.eq.1) then
105	! write (,) ' diffusivity approx. beta = ',beta
106	mu(1) = 1.0d0
107	amu(1)= 1.0d0
108	c-------data for 8 points integration
109	elseif (iimu.eq.4) then
110	write (,)' 4 points for the gauss-legendre angle quadrature.'
111	mu(1)=(1.0d0+0.339981043584856)/2.0d0
112	mu(2)=(1.0d0-0.339981043584856)/2.0d0
113	mu(3)=(1.0d0+0.861136311594053)/2.0d0
114	mu(4)=(1.0d0-0.861136311594053)/2.0d0
115	amu(1)=0.652145154862546
116	amu(2)=amu(1)
117	amu(3)=0.347854845137454
118	amu(4)=amu(3)
119	beta=1.0d0
120	c-------data for 8 points integration
121	elseif(iimu.eq.8) then
122	write (,)' 8 points for the gauss-legendre angle quadrature.'
123	mu(1)=(1.0d0+0.183434642495650)/2.0d0
124	mu(2)=(1.0d0-0.183434642495650)/2.0d0
125	mu(3)=(1.0d0+0.525532409916329)/2.0d0
126	mu(4)=(1.0d0-0.525532409916329)/2.0d0
127	mu(5)=(1.0d0+0.796666477413627)/2.0d0
128	mu(6)=(1.0d0-0.796666477413627)/2.0d0
129	mu(7)=(1.0d0+0.960289856497536)/2.0d0
130	mu(8)=(1.0d0-0.960289856497536)/2.0d0
131	amu(1)=0.362683783378362
132	amu(2)=amu(1)
133	amu(3)=0.313706645877887
134	amu(4)=amu(3)
135	amu(5)=0.222381034453374
136	amu(6)=amu(5)
137	amu(7)=0.101228536290376
138	amu(8)=amu(7)
139	beta=1.0d0
140	end if
141	c!!!!!!!!!!!!!!!!!!!!!!!
142
143	ccc
144	ccc determine abundances included in the absorber amount
145	ccc
146
147	c first, set up the grid ready for interpolation.
148	do i=1,nzy
149	zyd(i) = dble(zy(i))
150	enddo
151	do i=1,nl
152	zld(i) = dble(zl(i))
153	enddo
154
155	c vibr. temp of the bending mode :
156	if (isot.eq.1) call interdp ( tvtbs,zyd,nzy, v626t1,zld,nl, 1 )
157	if (isot.eq.2) call interdp ( tvtbs,zyd,nzy, v628t1,zld,nl, 1 )
158	if (isot.eq.3) call interdp ( tvtbs,zyd,nzy, v636t1,zld,nl, 1 )
159	if (isot.eq.4) call interdp ( tvtbs,zyd,nzy, v627t1,zld,nl, 1 )
160
161	c 2nd: correccion a la n10(i) (cantidad de absorbente en el lower state)
162	c por similitud a la que se hace en cza.for
163
164	do i=1,nzy
165	if (isot.eq.5) then
166	con(i) = dble( coy(i) * imrco )
167	else
168	con(i) = dble( co2y(i) * imr(isot) )
169	correc = 2.d0 * dexp( dble(-ee*elow(isot,2))/tvtbs(i) )
170	con(i) = con(i) * ( 1.d0 - correc )
171	endif
172	c-----------------------------------------------------------------------
173	c mlp & cristina. 17 july 1996
174	c change the calculation of mr. it is used for calculating partial press
175	c alpha = alpha(self,co2)pp +alpha(n2)(pt-pp)
176	c for an isotope, if mr is obtained by co2*imr(iso)/nt we are considerin
177	c collisions with other co2 isotopes (including the major one, 626)
178	c as if they were with n2. assuming mr as co2/nt, we consider collisions
179	c of type 628-626 as of 626-626 instead of as 626-n2.
180	c mrx(i)=con(i)/ntx(i) ! old malv
181
182	! mrx(i)= dble(co2x(i)/ntx(i)) ! mlp & crs
183
184	c jan 98:
185	c esta modif de mlp implica anular el correc (deberia revisar esto)
186	mr(i) = dble(co2y(i)/nty(i)) ! malv, jan 98
187
188	c-----------------------------------------------------------------------
189
190	end do
191
192	! como beta y 1.d5 son comunes a todas las weighted absorber amounts,
193	! los simplificamos:
194	! coninf = beta * 1.d5 * dble( con(n) / log( con(n-1) / con(n) ) )
195	coninf = dble( con(nzy) / log( con(nzy-1) / con(nzy) ) )
196
197	! write (,) ' coninf =', coninf
198
199	ccc
200	ccc temp dependence of the band strength and
201	ccc nlte correction factor for the absorber amount
202	ccc
203	call mztf_correccion ( coninf, con, ib, isot, 0 )
204
205	ccc
206	ccc reads histogrammed spectral data (strength for lte and vmr=1)
207	ccc
208	!hfile1 = dirspec//'hi'//dn ! Ya no distinguimos entre d/n
209	!! hfile1 = dirspec//'hid' ! (see why in his.for)
210	! hfile1='hid'
211	!! if (ib.eq.13 .or. ib.eq.14 ) hfile1 = dirspec//'his'
212	! if (ib.eq.13 .or. ib.eq.14 ) hfile1 = 'his'
213
214	! if(ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5
215	! @ .or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9) then
216	! if (isot.eq.1) hisfile = hfile1//'26-'//ibcode1//'.dat'
217	! if (isot.eq.2) hisfile = hfile1//'28-'//ibcode1//'.dat'
218	! if (isot.eq.3) hisfile = hfile1//'36-'//ibcode1//'.dat'
219	! if (isot.eq.4) hisfile = hfile1//'27-'//ibcode1//'.dat'
220	! if (isot.eq.5) hisfile = hfile1//'62-'//ibcode1//'.dat'
221	! else
222	! if (isot.eq.1) hisfile = hfile1//'26-'//ibcode2//'.dat'
223	! if (isot.eq.2) hisfile = hfile1//'28-'//ibcode2//'.dat'
224	! if (isot.eq.3) hisfile = hfile1//'36-'//ibcode2//'.dat'
225	! if (isot.eq.4) hisfile = hfile1//'27-'//ibcode2//'.dat'
226	! if (isot.eq.5) hisfile = hfile1//'62-'//ibcode2//'.dat'
227	! endif
228	!write (,) ' /MZESCAPE/ hisfile: ', hisfile
229
230	! the argument to rhist is to make this compatible with mztf_comp.f,
231	! which is a useful modification of mztf.f (to change strengths of bands
232	! call rhist (1.0)
233	if(ib.eq.1) then
234	if(isot.eq.1) then !Case 1
235	mm=mm_c1
236	nbox=nbox_c1
237	tmin=tmin_c1
238	tmax=tmax_c1
239	do i=1,nbox_max
240	no(i)=no_c1(i)
241	dist(i)=dist_c1(i)
242	do j=1,nhist
243	sk1(j,i)=sk1_c1(j,i)
244	xls1(j,i)=xls1_c1(j,i)
245	xln1(j,i)=xln1_c1(j,i)
246	xld1(j,i)=xld1_c1(j,i)
247	enddo
248	enddo
249	do j=1,nhist
250	thist(j)=thist_c1(j)
251	enddo
252	else if(isot.eq.2) then !Case 2
253	mm=mm_c2
254	nbox=nbox_c2
255	tmin=tmin_c2
256	tmax=tmax_c2
257	do i=1,nbox_max
258	no(i)=no_c2(i)
259	dist(i)=dist_c2(i)
260	do j=1,nhist
261	sk1(j,i)=sk1_c2(j,i)
262	xls1(j,i)=xls1_c2(j,i)
263	xln1(j,i)=xln1_c2(j,i)
264	xld1(j,i)=xld1_c2(j,i)
265	enddo
266	enddo
267	do j=1,nhist
268	thist(j)=thist_c2(j)
269	enddo
270	else if(isot.eq.3) then !Case 3
271	mm=mm_c3
272	nbox=nbox_c3
273	tmin=tmin_c3
274	tmax=tmax_c3
275	do i=1,nbox_max
276	no(i)=no_c3(i)
277	dist(i)=dist_c3(i)
278	do j=1,nhist
279	sk1(j,i)=sk1_c3(j,i)
280	xls1(j,i)=xls1_c3(j,i)
281	xln1(j,i)=xln1_c3(j,i)
282	xld1(j,i)=xld1_c3(j,i)
283	enddo
284	enddo
285	do j=1,nhist
286	thist(j)=thist_c3(j)
287	enddo
288	else if(isot.eq.4) then !Case 4
289	mm=mm_c4
290	nbox=nbox_c4
291	tmin=tmin_c4
292	tmax=tmax_c4
293	do i=1,nbox_max
294	no(i)=no_c4(i)
295	dist(i)=dist_c4(i)
296	do j=1,nhist
297	sk1(j,i)=sk1_c4(j,i)
298	xls1(j,i)=xls1_c4(j,i)
299	xln1(j,i)=xln1_c4(j,i)
300	xld1(j,i)=xld1_c4(j,i)
301	enddo
302	enddo
303	do j=1,nhist
304	thist(j)=thist_c4(j)
305	enddo
306	else
307	write(,)'isot must be 2,3 or 4 for ib=1!!'
308	write(,)'stop at mzescape/312'
309	stop
310	endif
311	else if (ib.eq.2) then
312	if(isot.eq.1) then !Case 5
313	mm=mm_c5
314	nbox=nbox_c5
315	tmin=tmin_c5
316	tmax=tmax_c5
317	do i=1,nbox_max
318	no(i)=no_c5(i)
319	dist(i)=dist_c5(i)
320	do j=1,nhist
321	sk1(j,i)=sk1_c5(j,i)
322	xls1(j,i)=xls1_c5(j,i)
323	xln1(j,i)=xln1_c5(j,i)
324	xld1(j,i)=xld1_c5(j,i)
325	enddo
326	enddo
327	do j=1,nhist
328	thist(j)=thist_c5(j)
329	enddo
330	else
331	write(,)'isot must be 1 for ib=2!!'
332	write(,)'stop at mzescape/336'
333	stop
334	endif
335	else if (ib.eq.3) then
336	if(isot.eq.1) then !Case 6
337	mm=mm_c6
338	nbox=nbox_c6
339	tmin=tmin_c6
340	tmax=tmax_c6
341	do i=1,nbox_max
342	no(i)=no_c6(i)
343	dist(i)=dist_c6(i)
344	do j=1,nhist
345	sk1(j,i)=sk1_c6(j,i)
346	xls1(j,i)=xls1_c6(j,i)
347	xln1(j,i)=xln1_c6(j,i)
348	xld1(j,i)=xld1_c6(j,i)
349	enddo
350	enddo
351	do j=1,nhist
352	thist(j)=thist_c6(j)
353	enddo
354	else
355	write(,)'isot must be 1 for ib=3!!'
356	write(,)'stop at mzescape/360'
357	stop
358	endif
359	else if (ib.eq.4) then
360	if(isot.eq.1) then !Case 7
361	mm=mm_c7
362	nbox=nbox_c7
363	tmin=tmin_c7
364	tmax=tmax_c7
365	do i=1,nbox_max
366	no(i)=no_c7(i)
367	dist(i)=dist_c7(i)
368	do j=1,nhist
369	sk1(j,i)=sk1_c7(j,i)
370	xls1(j,i)=xls1_c7(j,i)
371	xln1(j,i)=xln1_c7(j,i)
372	xld1(j,i)=xld1_c7(j,i)
373	enddo
374	enddo
375	do j=1,nhist
376	thist(j)=thist_c7(j)
377	enddo
378	else
379	write(,)'isot must be 1 for ib=4!!'
380	write(,)'stop at mzescape/384'
381	stop
382	endif
383	else
384	write(,)'ib must be 1,2,3 or 4!!'
385	write(,)'stop at mzescape/389'
386	endif
387	if (isot.ne.5) deltanux = deltanu(isot,ib)
388	if (isot.eq.5) deltanux = deltanuco
389
390	c******
391	c****** calculation of tauinf(nl)
392	c******
393	call initial
394
395	ff=1.0e10
396
397	do i=nl,1,-1
398
399	if(i.eq.nl)then
400
401	call intz (zl(i),c2,p2,mr2,t2, con)
402	do kr=1,nbox
403	ta(kr)=t2
404	end do
405	! write (,) ' i, t2 =', i, t2
406	call interstrength (st2,t2,ka,ta)
407	aa = p2 * coninf * mr2 * (st2 * ff)
408	bb = p2 * coninf * st2
409	cc = coninf * st2
410	dd = t2 * coninf * st2
411	do kr=1,nbox
412	ccbox(kr) = coninf * ka(kr)
413	ddbox(kr) = t2 * ccbox(kr)
414	! c2box(kr) = c2 * ka(kr) * beta * dble(deltaz) * 1.d5
415	c2box(kr) = c2 * ka(kr) * dble(deltaz)
416	end do
417	! c2 = c2 * st2 * beta * dble(deltaz) * 1.d5
418	c2 = c2 * st2 * dble(deltaz)
419
420	else
421	call intz (zl(i),c1,p1,mr1,t1, con)
422	do kr=1,nbox
423	ta(kr)=t1
424	end do
425	! write (,) ' i, t1 =', i, t1
426	call interstrength (st1,t1,ka,ta)
427	do kr=1,nbox
428	! c1box(kr) = c1 * ka(kr) * beta * dble(deltaz) * 1.d5
429	c1box(kr) = c1 * ka(kr) * dble(deltaz)
430	end do
431	! c1 = c1 * st1 * beta * dble(deltaz) * 1.d5
432	c1 = c1 * st1 * dble(deltaz)
433	aa = aa + ( p1mr1(c1ff) + p2mr2(c2ff)) / 2.d0
434	bb = bb + ( p1c1 + p2c2 ) / 2.d0
435	cc = cc + ( c1 + c2 ) / 2.d0
436	ccc = ( c1 + c2 ) / 2.d0
437	dd = dd + ( t1c1 + t2c2 ) / 2.d0
438	do kr=1,nbox
439	ccbox(kr) = ccbox(kr) +
440	@ ( c1box(kr) + c2box(kr) )/2.d0
441	ddbox(kr) = ddbox(kr) +
442	@ ( t1c1box(kr)+t2c2box(kr) )/2.d0
443	end do
444
445	mr2 = mr1
446	c2=c1
447	do kr=1,nbox
448	c2box(kr) = c1box(kr)
449	end do
450	t2=t1
451	p2=p1
452	end if
453
454	pt = bb / cc
455	pp = aa / (cc*ff)
456
457	! ta=dd/cc
458	! tdop = ta
459	ts = dd/cc
460	do kr=1,nbox
461	ta(kr) = ddbox(kr) / ccbox(kr)
462	end do
463	! write (,) ' i, ts =', i, ts
464	call interstrength(st,ts,ka,ta)
465	! call intershape(alsa,alna,alda,tdop)
466	call intershape(alsa,alna,alda,ta)
467
468	* ua = cc/st
469
470	c next loop calculates the eqw for an especified path ua,pp,pt,ta
471
472	eqwmu = 0.0d0
473	do im = 1,iimu
474	eqw=0.0d0
475	do kr=1,nbox
476	ua(kr) = ccbox(kr) / ka(kr) * beta * 1.0d5 / mu(im)
477	if(ua(kr).lt.0.)write(,)'mzescape/480',ua(kr),
478	$ ccbox(kr),ka(kr),beta,mu(im),kr,im,i,nl
479
480	call findw (ig,iirw, 0, csL,psL, Desp, wsL)
481	if ( i_supersat .eq. 0 ) then
482	eqw=eqw+no(kr)*w
483	else
484	eqw = w + (no(kr)-1) * ( asat_box*dist(kr) )
485	endif
486	end do
487	eqwmu = eqwmu + eqw * mu(im)*amu(im)
488	end do
489
490	! tauinf(i) = exp( - eqwmu / dble(deltanux) )
491	tauinf(i) = 1.d0 - eqwmu / dble(deltanux)
492	if (tauinf(i).lt.0.d0) tauinf(i) = 0.0d0
493
494	if (i.eq.nl) then
495	taustar(i) = 0.0d0
496	else
497	taustar(i) = dble(deltanux) * (tauinf(i+1)-tauinf(i))
498	! ~ / ( beta * cc * 1.d5 )
499	~ / ( beta * ccc * 1.d5 )
500	endif
501
502	end do ! i continue
503
504
505	c******
506	c****** calculation of tau(in,ir) for n<=r
507	c******
508
509	do 1 in=1,nl-1
510
511	call initial
512
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) * dble(deltaz)
520	end do
521	c1 = c1 * st1 * dble(deltaz)
522
523	call intz (zl(in+1), c2,p2,mr2,t2, con)
524	do kr=1,nbox
525	ta(kr) = t2
526	end do
527	call interstrength (st2,t2,ka,ta)
528	do kr=1,nbox
529	c2box(kr) = c2 * ka(kr) * dble(deltaz)
530	end do
531	c2 = c2 * st2 * dble(deltaz)
532
533	aa = aa + ( p1mr1(c1ff) + p2mr2(c2ff)) / 2.d0
534	bb = bb + ( p1c1 + p2c2 ) / 2.d0
535	cc = cc + ( c1 + c2 ) / 2.d0
536	dd = dd + ( t1c1 + t2c2 ) / 2.d0
537	do kr=1,nbox
538	ccbox(kr) = ccbox(kr) + (c1box(kr)+c2box(kr))/2.d0
539	ddbox(kr) = ddbox(kr) + (t1c1box(kr)+t2c2box(kr))/2.d0
540	end do
541
542	mr1=mr2
543	t1=t2
544	c1=c2
545	p1=p2
546	do kr=1,nbox
547	c1box(kr) = c2box(kr)
548	end do
549	pt = bb / cc
550	pp = aa / (cc * ff)
551	ts = dd/cc
552	do kr=1,nbox
553	ta(kr) = ddbox(kr) / ccbox(kr)
554	end do
555	call interstrength(st,ts,ka,ta)
556	call intershape(alsa,alna,alda,ta)
557
558	eqwmu = 0.0d0
559	do im = 1,iimu
560	eqw=0.0d0
561	do kr=1,nbox
562	ua(kr) = ccbox(kr) / ka(kr) * beta * 1.0d5 / mu(im)
563	if(ua(kr).lt.0.)write(,)'mzescape/566',ua(kr),
564	$ ccbox(kr),ka(kr),beta,mu(im),kr,im,i,nl
565
566	call findw (ig,iirw, 0, csL,psL, Desp, wsL)
567	if ( i_supersat .eq. 0 ) then
568	eqw=eqw+no(kr)*w
569	else
570	eqw = w + (no(kr)-1) * ( asat_box*dist(kr) )
571	endif
572	end do
573	eqwmu = eqwmu + eqw * mu(im)*amu(im)
574	end do
575
576	tauii(in) = exp( - eqwmu / dble(deltanux) )
577	!write (,) 'i,tauii=',in,tauii(in)
578
579	1 continue
580	tauii(nl) = 1.0d0
581
582
583	c end
584	return
585	end
586
587
588
589	c***********************************************************************
590	c mzescape_normaliz.f
591	c***********************************************************************
592	c
593	c program for correcting some strange values and for normalizing
594	c the atmospheric escape functions computed by mzescape_15um.f
595	c possibilities according to istyle (see mzescape_15um.f).
596	c
597
598	subroutine mzescape_normaliz ( taustar, istyle )
599
600
601	c dic 99 malv first version
602	c jul 2011 malv+fgg Adapted to LMD-MGCM
603	c***********************************************************************
604
605	implicit none
606	include 'nlte_paramdef.h'
607	include 'nlte_commons.h'
608
609
610	c arguments
611	real*8 taustar(nl) ! o
612	integer istyle ! i
613
614	c local variables and constants
615	integer i, imaximum
616	real*8 maximum
617
618	c***********************************************************************
619
620	!
621	! correcting strange values at top, eliminating local maxima, etc...
622	!
623	taustar(nl) = taustar(nl-1)
624
625	if ( istyle .eq. 1 ) then
626	imaximum = nl
627	maximum = taustar(nl)
628	do i=1,nl-1
629	if (taustar(i).gt.maximum) taustar(i) = taustar(nl)
630	enddo
631	elseif ( istyle .eq. 2 ) then
632	imaximum = nl
633	maximum = taustar(nl)
634	do i=nl-1,1,-1
635	if (taustar(i).gt.maximum) then
636	maximum = taustar(i)
637	imaximum = i
638	endif
639	enddo
640	do i=imaximum,nl
641	if (taustar(i).lt.maximum) taustar(i) = maximum
642	enddo
643	endif
644
645	!
646	! normalizing
647	!
648	do i=1,nl
649	taustar(i) = taustar(i) / maximum
650	enddo
651
652
653	c end
654	return
655	end
656
657
658
659	c***********************************************************************
660	c mzescape_fb.f
661	c***********************************************************************
662	subroutine mzescape_fb(ig)
663
664	c computes the escape functions of the most important 15um bands
665	c this calls mzescape ( taustar,tauinf,tauii, ib,isot, iirw,iimu
666
667	c nov 99 malv based on cm15um_fb.f
668	c jul 2011 malv+fgg adapted to LMD-MGCM
669	c***********************************************************************
670
671	implicit none
672
673	include 'nlte_paramdef.h'
674	include 'nlte_commons.h'
675
676	c local variables
677	integer i, ib, ik, istyle
678	integer ig !ADDED FOR TRACEBACK
679	real*8 tau_factor
680	real*8 aux(nl), aux2(nl), aux3(nl)
681
682	c***********************************************************************
683
684	call mzescape (ig,taustar21,tauinf210,tauii210,1,2,irw_mztf,imu)
685	call mzescape (ig,taustar31,tauinf310,tauii310,1,3,irw_mztf,imu)
686	call mzescape (ig,taustar41,tauinf410,tauii410,1,4,irw_mztf,imu)
687
688	istyle = 2
689	call mzescape_normaliz ( taustar21, istyle )
690	call mzescape_normaliz ( taustar31, istyle )
691	call mzescape_normaliz ( taustar41, istyle )
692
693
694	c end
695	return
696	end
697
698
699
700	c***********************************************************************
701	c mzescape_fh.f
702	c***********************************************************************
703	subroutine mzescape_fh(ig)
704
705	c jul 2011 malv+fgg
706	c***********************************************************************
707
708	implicit none
709
710	include 'nlte_paramdef.h'
711	include 'nlte_commons.h'
712
713	c local variables
714	integer i, ib, ik, istyle
715	integer ig ! ADDED FOR TRACEBACK
716	real*8 tau_factor
717	real*8 aux(nl), aux2(nl), aux3(nl)
718
719	c***********************************************************************
720
721	call zero4v( aux, taustar12,tauinf121,tauii121, nl)
722	do ik=1,3
723	ib=ik+1
724	call mzescape ( ig,aux,aux2,aux3, ib, 1,irw_mztf,imu )
725	tau_factor = 1.d0
726	if (ik.eq.1) tau_factor = dble(667.75/618.03)
727	if (ik.eq.3) tau_factor = dble(667.75/720.806)
728	do i=1,nl
729	taustar12(i) = taustar12(i) + aux(i) * tau_factor
730	tauinf121(i) = tauinf121(i) + aux2(i) * tau_factor
731	tauii121(i) = tauii121(i) + aux3(i) * tau_factor
732	enddo
733	enddo
734
735	istyle = 2
736	call mzescape_normaliz ( taustar12, istyle )
737
738
739
740	c end
741	return
742	end
743
744
745
746
747
748	c***********************************************************************
749	c mztud.f
750	c***********************************************************************
751
752	subroutine mztud ( ig,cf,cfup,cfdw,vc,taugr, ib,isot,
753	@ iirw,iimu,itauout,icfout,itableout )
754
755	c program for calculating atmospheric transmittances
756	c to be used in the calculation of curtis matrix coefficients
757	c i*out = 1 output of data
758	c i*out = 0 no output
759	c itableout = 30 output de toda la C.M. y el VC y las poblaciones de los
760	c estados 626(020), esta opcion nueva se añade porque
761	c itableout=1 saca o bien solamente de 5 en 5 capas
762	c o bien los elementos de C.M. desde una cierta capa
763	c (consultese elimin_mz1d.f que es quien lo hace); lo
764	c de las poblaciones (020) lo hace mztf_correcion.f
765
766	c jul 2011 malv+fgg Adapted to LMD-MGCM
767	c jan 07 malv Add new vertical fine grid zy, similar to zx
768	c sep-oct 01 malv update for fluxes for hb and fb, adapt to Linux
769	c nov 98 mavl allow for overlaping in the lorentz line
770	c jan 98 malv version for mz1d. based on curtis/mztf.for
771	c 17-jul-96 mlp&crs change the calculation of mr.
772	c evitar: divide por cero. anhadiendo: ff
773	c oct-92 malv correct s(t) dependence for all histogr bands
774	c june-92 malv proper lower levels for laser bands
775	c may-92 malv new temperature dependence for laser bands
776	c @ 991 malv boxing for the averaged absorber amount and t
777	c ? malv extension up to 200 km altitude in mars
778	c 13-nov-86 mlp include the temperature weighted to match
779	c the eqw in the strong doppler limit.
780	c***********************************************************************
781
782	implicit none
783
784	include 'nlte_paramdef.h'
785	include 'nlte_commons.h'
786
787	c arguments
788	integer ig !ADDED FOR TRACEBACK
789	real*8 cf(nl,nl), cfup(nl,nl), cfdw(nl,nl) ! o
790	real*8 vc(nl), taugr(nl) ! o
791	integer ib ! i
792	integer isot ! i
793	integer iirw ! i
794	integer iimu ! i
795	integer itauout ! i
796	integer icfout ! i
797	integer itableout ! i
798
799	c local variables and constants
800	integer i, in, ir, im, k,j
801	integer nmu
802	parameter (nmu = 8)
803	real*8 tau(nl,nl)
804	real*8 tauinf(nl)
805	real*8 con(nzy), coninf
806	real*8 c1, c2
807	real*8 t1, t2
808	real*8 p1, p2
809	real*8 mr1, mr2
810	real*8 st1, st2
811	real*8 c1box(70), c2box(70)
812	real*8 ff ! to avoid too small numbers
813	real*8 tvtbs(nzy)
814	real*8 st, beta, ts, eqwmu
815	real*8 mu(nmu), amu(nmu)
816	real*8 zld(nl), zyd(nzy)
817	real*8 correc
818	real deltanux ! width of vib-rot band (cm-1)
819	character isotcode*2
820	integer idummy
821	real*8 Desp,wsL
822
823	c formats
824	111 format(a1)
825	112 format(a2)
826	101 format(i1)
827	202 format(i2)
828	180 format(a80)
829	181 format(a80)
830	c***********************************************************************
831
832	c some needed values
833	! rl=sqrt(log(2.d0))
834	! pi2 = 3.14159265358989d0
835	beta = 1.8d0
836	! beta = 1.0d0
837	idummy = 0
838	Desp = 0.0d0
839	wsL = 0.0d0
840
841	! write (,) ' MZTUD/ iirw = ', iirw
842
843
844	c esto es para que las subroutines de mztfsub calculen we
845	c de la forma apropiada para mztf, no para fot
846	icls=icls_mztf
847
848	c codigos para filenames
849	! if (isot .eq. 1) isotcode = '26'
850	! if (isot .eq. 2) isotcode = '28'
851	! if (isot .eq. 3) isotcode = '36'
852	! if (isot .eq. 4) isotcode = '27'
853	! if (isot .eq. 5) isotcode = '62'
854	! if(ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5
855	! @ .or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9) then
856	! write (ibcode1,101) ib
857	! else
858	! write (ibcode2,202) ib
859	! endif
860	! write (*,'( 30h calculating curtis matrix : ,2x,
861	! @ 8h band = ,i2,2x, 11h isotope = ,i2)') ib, isot
862
863	c integration in angle !!!!!!!!!!!!!!!!!!!!
864	c------- diffusivity approx.
865	if (iimu.eq.1) then
866	! write (,) ' diffusivity approx. beta = ',beta
867	mu(1) = 1.0d0
868	amu(1)= 1.0d0
869	c-------data for 8 points integration
870	elseif (iimu.eq.4) then
871	write (,)' 4 points for the gauss-legendre angle quadrature.'
872	mu(1)=(1.0d0+0.339981043584856)/2.0d0
873	mu(2)=(1.0d0-0.339981043584856)/2.0d0
874	mu(3)=(1.0d0+0.861136311594053)/2.0d0
875	mu(4)=(1.0d0-0.861136311594053)/2.0d0
876	amu(1)=0.652145154862546
877	amu(2)=amu(1)
878	amu(3)=0.347854845137454
879	amu(4)=amu(3)
880	beta=1.0d0
881	c-------data for 8 points integration
882	elseif(iimu.eq.8) then
883	write (,)' 8 points for the gauss-legendre angle quadrature.'
884	mu(1)=(1.0d0+0.183434642495650)/2.0d0
885	mu(2)=(1.0d0-0.183434642495650)/2.0d0
886	mu(3)=(1.0d0+0.525532409916329)/2.0d0
887	mu(4)=(1.0d0-0.525532409916329)/2.0d0
888	mu(5)=(1.0d0+0.796666477413627)/2.0d0
889	mu(6)=(1.0d0-0.796666477413627)/2.0d0
890	mu(7)=(1.0d0+0.960289856497536)/2.0d0
891	mu(8)=(1.0d0-0.960289856497536)/2.0d0
892	amu(1)=0.362683783378362
893	amu(2)=amu(1)
894	amu(3)=0.313706645877887
895	amu(4)=amu(3)
896	amu(5)=0.222381034453374
897	amu(6)=amu(5)
898	amu(7)=0.101228536290376
899	amu(8)=amu(7)
900	beta=1.0d0
901	end if
902	c!!!!!!!!!!!!!!!!!!!!!!!
903
904	ccc
905	ccc determine abundances included in the absorber amount
906	ccc
907
908	c first, set up the grid ready for interpolation.
909	do i=1,nzy
910	zyd(i) = dble(zy(i))
911	enddo
912	do i=1,nl
913	zld(i) = dble(zl(i))
914	enddo
915	c vibr. temp of the bending mode :
916	if (isot.eq.1) call interdp ( tvtbs,zyd,nzy, v626t1,zld,nl, 1 )
917	if (isot.eq.2) call interdp ( tvtbs,zyd,nzy, v628t1,zld,nl, 1 )
918	if (isot.eq.3) call interdp ( tvtbs,zyd,nzy, v636t1,zld,nl, 1 )
919	if (isot.eq.4) call interdp ( tvtbs,zyd,nzy, v627t1,zld,nl, 1 )
920	!if (isot.eq.5) call interdp ( tvtbs,zxd,nz, vcot1,zld,nl, 1 )
921
922	c 2nd: correccion a la n10(i) (cantidad de absorbente en el lower state)
923	c por similitud a la que se hace en cza.for ; esto solo se hace para CO2
924	!write (,) 'imr(isot) = ', isot, imr(isot)
925	do i=1,nzy
926	if (isot.eq.5) then
927	con(i) = dble( coy(i) * imrco )
928	else
929	con(i) = dble( co2y(i) * imr(isot) )
930	correc = 2.d0 * dexp( dble(-ee*elow(isot,2))/tvtbs(i) )
931	con(i) = con(i) * ( 1.d0 - correc )
932	! write (,) ' iz, correc, co2y(i), con(i) =',
933	! @ i,correc,co2y(i),con(i)
934	endif
935
936	!-----------------------------------------------------------------
937	! mlp & cristina. 17 july 1996 change the calculation of mr.
938	! it is used for calculating partial press
939	! alpha = alpha(self,co2)pp +alpha(n2)(pt-pp)
940	! for an isotope, if mr is obtained by
941	! co2*imr(iso)/nt
942	! we are considerin collisions with other co2 isotopes
943	! (including the major one, 626) as if they were with n2.
944	! assuming mr as co2/nt, we consider collisions
945	! of type 628-626 as of 626-626 instead of as 626-n2.
946	! mrx(i)=con(i)/ntx(i) ! old malv
947	! mrx(i)= dble(co2x(i)/ntx(i)) ! mlp & crs
948
949	! jan 98:
950	! esta modif de mlp implica anular el correc (deberia revisar esto)
951
952	mr(i) = dble(co2y(i)/nty(i)) ! malv, jan 98
953
954	!-----------------------------------------------------------------
955
956	end do
957
958	! como beta y 1.d5 son comunes a todas las weighted absorber amounts,
959	! los simplificamos:
960	! coninf = beta * 1.d5 * dble( con(n) / log( con(n-1) / con(n) ) )
961	!write (,) ' con(nz), con(nz-1) =', con(nz), con(nz-1)
962	coninf = dble( con(nzy) / log( con(nzy-1) / con(nzy) ) )
963	!write (,) ' coninf =', coninf
964
965	ccc
966	ccc temp dependence of the band strength and
967	ccc nlte correction factor for the absorber amount
968	ccc
969	call mztf_correccion ( coninf, con, ib, isot, itableout )
970	ccc
971	ccc reads histogrammed spectral data (strength for lte and vmr=1)
972	ccc
973	!hfile1 = dirspec//'hi'//dn !Ya no hacemos distincion d/n en esto
974	! hfile1 = dirspec//'hid' !(see why in his.for)
975	! hfile1='hid'
976	!! if (ib.eq.13 .or. ib.eq.14 ) hfile1 = dirspec//'his'
977	! if (ib.eq.13 .or. ib.eq.14 ) hfile1 = 'his'
978
979	! if(ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5
980	! @ .or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9) then
981	! if (isot.eq.1) hisfile = hfile1//'26-'//ibcode1//'.dat'
982	! if (isot.eq.2) hisfile = hfile1//'28-'//ibcode1//'.dat'
983	! if (isot.eq.3) hisfile = hfile1//'36-'//ibcode1//'.dat'
984	! if (isot.eq.4) hisfile = hfile1//'27-'//ibcode1//'.dat'
985	! if (isot.eq.5) hisfile = hfile1//'62-'//ibcode1//'.dat'
986	! else
987	! if (isot.eq.1) hisfile = hfile1//'26-'//ibcode2//'.dat'
988	! if (isot.eq.2) hisfile = hfile1//'28-'//ibcode2//'.dat'
989	! if (isot.eq.3) hisfile = hfile1//'36-'//ibcode2//'.dat'
990	! if (isot.eq.4) hisfile = hfile1//'27-'//ibcode2//'.dat'
991	! if (isot.eq.5) hisfile = hfile1//'62-'//ibcode2//'.dat'
992	! endif
993	if(ib.eq.1) then
994	if(isot.eq.1) then !Case 1
995	mm=mm_c1
996	nbox=nbox_c1
997	tmin=tmin_c1
998	tmax=tmax_c1
999	do i=1,nbox_max
1000	no(i)=no_c1(i)
1001	dist(i)=dist_c1(i)
1002	do j=1,nhist
1003	sk1(j,i)=sk1_c1(j,i)
1004	xls1(j,i)=xls1_c1(j,i)
1005	xln1(j,i)=xln1_c1(j,i)
1006	xld1(j,i)=xld1_c1(j,i)
1007	enddo
1008	enddo
1009	do j=1,nhist
1010	thist(j)=thist_c1(j)
1011	enddo
1012	else if(isot.eq.2) then !Case 2
1013	mm=mm_c2
1014	nbox=nbox_c2
1015	tmin=tmin_c2
1016	tmax=tmax_c2
1017	do i=1,nbox_max
1018	no(i)=no_c2(i)
1019	dist(i)=dist_c2(i)
1020	do j=1,nhist
1021	sk1(j,i)=sk1_c2(j,i)
1022	xls1(j,i)=xls1_c2(j,i)
1023	xln1(j,i)=xln1_c2(j,i)
1024	xld1(j,i)=xld1_c2(j,i)
1025	enddo
1026	enddo
1027	do j=1,nhist
1028	thist(j)=thist_c2(j)
1029	enddo
1030	else if(isot.eq.3) then !Case 3
1031	mm=mm_c3
1032	nbox=nbox_c3
1033	tmin=tmin_c3
1034	tmax=tmax_c3
1035	do i=1,nbox_max
1036	no(i)=no_c3(i)
1037	dist(i)=dist_c3(i)
1038	do j=1,nhist
1039	sk1(j,i)=sk1_c3(j,i)
1040	xls1(j,i)=xls1_c3(j,i)
1041	xln1(j,i)=xln1_c3(j,i)
1042	xld1(j,i)=xld1_c3(j,i)
1043	enddo
1044	enddo
1045	do j=1,nhist
1046	thist(j)=thist_c3(j)
1047	enddo
1048	else if(isot.eq.4) then !Case 4
1049	mm=mm_c4
1050	nbox=nbox_c4
1051	tmin=tmin_c4
1052	tmax=tmax_c4
1053	do i=1,nbox_max
1054	no(i)=no_c4(i)
1055	dist(i)=dist_c4(i)
1056	do j=1,nhist
1057	sk1(j,i)=sk1_c4(j,i)
1058	xls1(j,i)=xls1_c4(j,i)
1059	xln1(j,i)=xln1_c4(j,i)
1060	xld1(j,i)=xld1_c4(j,i)
1061	enddo
1062	enddo
1063	do j=1,nhist
1064	thist(j)=thist_c4(j)
1065	enddo
1066	else
1067	write(,)'isot must be 2,3 or 4 for ib=1!!'
1068	write(,)'stop at mztud/324'
1069	stop
1070	endif
1071	else if (ib.eq.2) then
1072	if(isot.eq.1) then !Case 5
1073	mm=mm_c5
1074	nbox=nbox_c5
1075	tmin=tmin_c5
1076	tmax=tmax_c5
1077	do i=1,nbox_max
1078	no(i)=no_c5(i)
1079	dist(i)=dist_c5(i)
1080	do j=1,nhist
1081	sk1(j,i)=sk1_c5(j,i)
1082	xls1(j,i)=xls1_c5(j,i)
1083	xln1(j,i)=xln1_c5(j,i)
1084	xld1(j,i)=xld1_c5(j,i)
1085	enddo
1086	enddo
1087	do j=1,nhist
1088	thist(j)=thist_c5(j)
1089	enddo
1090	else
1091	write(,)'isot must be 1 for ib=2!!'
1092	write(,)'stop at mztud/348'
1093	stop
1094	endif
1095	else if (ib.eq.3) then
1096	if(isot.eq.1) then !Case 6
1097	mm=mm_c6
1098	nbox=nbox_c6
1099	tmin=tmin_c6
1100	tmax=tmax_c6
1101	do i=1,nbox_max
1102	no(i)=no_c6(i)
1103	dist(i)=dist_c6(i)
1104	do j=1,nhist
1105	sk1(j,i)=sk1_c6(j,i)
1106	xls1(j,i)=xls1_c6(j,i)
1107	xln1(j,i)=xln1_c6(j,i)
1108	xld1(j,i)=xld1_c6(j,i)
1109	enddo
1110	enddo
1111	do j=1,nhist
1112	thist(j)=thist_c6(j)
1113	enddo
1114	else
1115	write(,)'isot must be 1 for ib=3!!'
1116	write(,)'stop at mztud/372'
1117	stop
1118	endif
1119	else if (ib.eq.4) then
1120	if(isot.eq.1) then !Case 7
1121	mm=mm_c7
1122	nbox=nbox_c7
1123	tmin=tmin_c7
1124	tmax=tmax_c7
1125	do i=1,nbox_max
1126	no(i)=no_c7(i)
1127	dist(i)=dist_c7(i)
1128	do j=1,nhist
1129	sk1(j,i)=sk1_c7(j,i)
1130	xls1(j,i)=xls1_c7(j,i)
1131	xln1(j,i)=xln1_c7(j,i)
1132	xld1(j,i)=xld1_c7(j,i)
1133	enddo
1134	enddo
1135	do j=1,nhist
1136	thist(j)=thist_c7(j)
1137	enddo
1138	else
1139	write(,)'isot must be 1 for ib=4!!'
1140	write(,)'stop at mztud/396'
1141	stop
1142	endif
1143	else
1144	write(,)'ib must be 1,2,3 or 4!!'
1145	write(,)'stop at mztud/401'
1146	endif
1147
1148
1149
1150
1151	! write (,) 'hisfile: ', hisfile
1152	! the argument to rhist is to make this compatible with mztf_comp.f,
1153	! which is a useful modification of mztf.f (to change strengths of bands
1154	! call rhist (1.0)
1155	if (isot.ne.5) deltanux = deltanu(isot,ib)
1156	if (isot.eq.5) deltanux = deltanuco
1157
1158	c******
1159	c****** calculation of tauinf(nl)
1160	c******
1161	call initial
1162	ff=1.0e10
1163
1164	do i=nl,1,-1
1165
1166	if(i.eq.nl)then
1167
1168	call intz (zl(i),c2,p2,mr2,t2, con)
1169	do kr=1,nbox
1170	ta(kr)=t2
1171	end do
1172	! write (,) ' i, t2 =', i, t2
1173	call interstrength (st2,t2,ka,ta)
1174	aa = p2 * coninf * mr2 * (st2 * ff)
1175	bb = p2 * coninf * st2
1176	cc = coninf * st2
1177	dd = t2 * coninf * st2
1178	do kr=1,nbox
1179	ccbox(kr) = coninf * ka(kr)
1180	ddbox(kr) = t2 * ccbox(kr)
1181	! c2box(kr) = c2 * ka(kr) * beta * dble(deltaz) * 1.d5
1182	c2box(kr) = c2 * ka(kr) * dble(deltaz)
1183	end do
1184	! c2 = c2 * st2 * beta * dble(deltaz) * 1.d5
1185	c2 = c2 * st2 * dble(deltaz)
1186
1187	else
1188	call intz (zl(i),c1,p1,mr1,t1, con)
1189	do kr=1,nbox
1190	ta(kr)=t1
1191	end do
1192	! write (,) ' i, t1 =', i, t1
1193	call interstrength (st1,t1,ka,ta)
1194	do kr=1,nbox
1195	! c1box(kr) = c1 * ka(kr) * beta * dble(deltaz) * 1.d5
1196	c1box(kr) = c1 * ka(kr) * dble(deltaz)
1197	end do
1198	! c1 = c1 * st1 * beta * dble(deltaz) * 1.d5
1199	c1 = c1 * st1 * dble(deltaz)
1200	aa = aa + ( p1mr1(c1ff) + p2mr2(c2ff)) / 2.d0
1201	bb = bb + ( p1c1 + p2c2 ) / 2.d0
1202	cc = cc + ( c1 + c2 ) / 2.d0
1203	dd = dd + ( t1c1 + t2c2 ) / 2.d0
1204	do kr=1,nbox
1205	ccbox(kr) = ccbox(kr) +
1206	@ ( c1box(kr) + c2box(kr) )/2.d0
1207	ddbox(kr) = ddbox(kr) +
1208	@ ( t1c1box(kr)+t2c2box(kr) )/2.d0
1209	end do
1210
1211	mr2 = mr1
1212	c2=c1
1213	do kr=1,nbox
1214	c2box(kr) = c1box(kr)
1215	end do
1216	t2=t1
1217	p2=p1
1218	end if
1219
1220	pt = bb / cc
1221	pp = aa / (cc*ff)
1222
1223	! ta=dd/cc
1224	! tdop = ta
1225	ts = dd/cc
1226	do kr=1,nbox
1227	ta(kr) = ddbox(kr) / ccbox(kr)
1228	end do
1229	! write (,) ' i, ts =', i, ts
1230	call interstrength(st,ts,ka,ta)
1231	! call intershape(alsa,alna,alda,tdop)
1232	call intershape(alsa,alna,alda,ta)
1233	* ua = cc/st
1234
1235	c next loop calculates the eqw for an especified path uapp,pt,ta
1236
1237	eqwmu = 0.0d0
1238	do im = 1,iimu
1239	eqw=0.0d0
1240	do kr=1,nbox
1241	ua(kr) = ccbox(kr) / ka(kr) * beta * 1.0d5 / mu(im)
1242	if(ua(kr).lt.0.)write(,)'mztud/504',ua(kr),ccbox(kr),
1243	$ ka(kr),beta,mu(im),kr,im,i,nl
1244	call findw(ig,iirw, idummy,c1,p1,Desp,wsL)
1245	if ( i_supersat .eq. 0 ) then
1246	eqw=eqw+no(kr)*w
1247	else
1248	eqw = w + (no(kr)-1) * ( asat_box*dist(kr) )
1249	endif
1250	end do
1251	eqwmu = eqwmu + eqw * mu(im)*amu(im)
1252	end do
1253
1254	tauinf(i) = exp( - eqwmu / dble(deltanux) )
1255
1256	end do
1257	! if ( isot.eq.1 .and. ib.eq.2 ) then
1258	! write (,) ' tauinf(nl) = ', tauinf(nl)
1259	! write (,) ' tauinf(1) = ', tauinf(1)
1260	! endif
1261
1262	c******
1263	c****** calculation of tau(in,ir) for n<=r
1264	c******
1265
1266	do 1 in=1,nl-1
1267	call initial
1268	call intz (zl(in), c1,p1,mr1,t1, con)
1269	do kr=1,nbox
1270	ta(kr) = t1
1271	end do
1272	call interstrength (st1,t1,ka,ta)
1273	do kr=1,nbox
1274	! c1box(kr) = c1 * ka(kr) * beta * dble(deltaz) * 1.d5
1275	c1box(kr) = c1 * ka(kr) * dble(deltaz)
1276	end do
1277	! c1 = c1 * st1 * beta * dble(deltaz) * 1.d5
1278	c1 = c1 * st1 * dble(deltaz)
1279
1280	do 2 ir=in,nl-1
1281
1282	if (ir.eq.in) then
1283	tau(in,ir) = 1.d0
1284	goto 2
1285	end if
1286
1287	call intz (zl(ir), c2,p2,mr2,t2, con)
1288	do kr=1,nbox
1289	ta(kr) = t2
1290	end do
1291	call interstrength (st2,t2,ka,ta)
1292	do kr=1,nbox
1293	! c2box(kr) = c2 * ka(kr) * beta * dble(deltaz) * 1.d5
1294	c2box(kr) = c2 * ka(kr) * dble(deltaz)
1295	end do
1296	! c2 = c2 * st2 * beta * dble(deltaz) * 1.e5
1297	c2 = c2 * st2 * dble(deltaz)
1298
1299	c aa = aa + ( p1mr1c1 + p2mr2c2 ) / 2.d0
1300	aa = aa + ( p1mr1(c1ff) + p2mr2(c2ff)) / 2.d0
1301	bb = bb + ( p1c1 + p2c2 ) / 2.d0
1302	cc = cc + ( c1 + c2 ) / 2.d0
1303	dd = dd + ( t1c1 + t2c2 ) / 2.d0
1304	do kr=1,nbox
1305	ccbox(kr) = ccbox(kr) + ( c1box(kr) + c2box(kr) ) /2.d0
1306	ddbox(kr) = ddbox(kr) +
1307	$ ( t1c1box(kr) + t2c2box(kr) ) /2.d0
1308	end do
1309
1310	mr1=mr2
1311	t1=t2
1312	c1=c2
1313	p1=p2
1314	do kr=1,nbox
1315	c1box(kr) = c2box(kr)
1316	end do
1317
1318	pt = bb / cc
1319	pp = aa / (cc * ff)
1320
1321	* ta=dd/cc
1322	* tdop = ta
1323	ts = dd/cc
1324	do kr=1,nbox
1325	ta(kr) = ddbox(kr) / ccbox(kr)
1326	end do
1327	call interstrength(st,ts,ka,ta)
1328	call intershape(alsa,alna,alda,ta)
1329	* ua = cc/st
1330
1331	c next loop calculates the eqw for an especified path ua,pp,pt,ta
1332
1333	eqwmu = 0.0d0
1334	do im = 1,iimu
1335	eqw=0.0d0
1336	do kr=1,nbox
1337	ua(kr) = ccbox(kr) / ka(kr) * beta * 1.0d5 / mu(im)
1338
1339	call findw(ig,iirw, idummy,c1,p1,Desp,wsL)
1340	if ( i_supersat .eq. 0 ) then
1341	eqw=eqw+no(kr)*w
1342	else
1343	eqw = w + (no(kr)-1) * ( asat_box*dist(kr) )
1344	endif
1345	end do
1346	eqwmu = eqwmu + eqw * mu(im)*amu(im)
1347	end do
1348
1349	tau(in,ir) = exp( - eqwmu / dble(deltanux) )
1350
1351	2 continue
1352
1353	1 continue
1354	! if ( isot.eq.1 .and. ib.eq.2 ) then
1355	! write (,) ' tau(1,) , =1,20 '
1356	! write (,) ( sngl(tau(1,k)), k=1,20 )
1357	! endif
1358
1359
1360	c**********
1361	c********** calculation of tau(in,ir) for n>r
1362	c**********
1363
1364	in=nl
1365
1366	call initial
1367	call intz (zl(in), c1,p1,mr1,t1, con)
1368	do kr=1,nbox
1369	ta(kr) = t1
1370	end do
1371	call interstrength (st1,t1,ka,ta)
1372	do kr=1,nbox
1373	! c1box(kr) = c1 * ka(kr) * beta * dble(deltaz) * 1.d5
1374	c1box(kr) = c1 * ka(kr) * dble(deltaz)
1375	end do
1376	! c1 = c1 * st1 * beta * dble(deltaz) * 1.d5
1377	c1 = c1 * st1 * dble(deltaz)
1378
1379	do 4 ir=in-1,1,-1
1380
1381	call intz (zl(ir), c2,p2,mr2,t2, con)
1382	do kr=1,nbox
1383	ta(kr) = t2
1384	end do
1385	call interstrength (st2,t2,ka,ta)
1386	do kr=1,nbox
1387	! c2box(kr) = c2 * ka(kr) * beta * dble(deltaz) * 1.d5
1388	c2box(kr) = c2 * ka(kr) * dble(deltaz)
1389	end do
1390	! c2 = c2 * st2 * beta * dble(deltaz) * 1.d5
1391	c2 = c2 * st2 * dble(deltaz)
1392
1393	aa = aa + ( p1mr1(c1ff) + p2mr2(c2ff)) / 2.d0
1394	bb = bb + ( p1c1 + p2c2 ) / 2.d0
1395	cc = cc + ( c1 + c2 ) / 2.d0
1396	dd = dd + ( t1c1 + t2c2 ) / 2.d0
1397	do kr=1,nbox
1398	ccbox(kr) = ccbox(kr) + ( c1box(kr) + c2box(kr) ) /2.d0
1399	ddbox(kr) = ddbox(kr) +
1400	$ ( t1c1box(kr) + t2c2box(kr) ) /2.d0
1401	end do
1402
1403	mr1=mr2
1404	c1=c2
1405	t1=t2
1406	p1=p2
1407	do kr=1,nbox
1408	c1box(kr) = c2box(kr)
1409	end do
1410
1411	pt = bb / cc
1412	pp = aa / (cc * ff)
1413	ts = dd / cc
1414	do kr=1,nbox
1415	ta(kr) = ddbox(kr) / ccbox(kr)
1416	end do
1417	call interstrength (st,ts,ka,ta)
1418	call intershape (alsa,alna,alda,ta)
1419
1420	* ua = cc/st
1421
1422	c next loop calculates the eqw for an especified path ua,pp,pt,ta
1423
1424	eqwmu = 0.0d0
1425	do im = 1,iimu
1426	eqw=0.0d0
1427	do kr=1,nbox
1428	ua(kr) = ccbox(kr) / ka(kr) * beta * 1.0d5 / mu(im)
1429	if(ua(kr).lt.0.)write(,)'mztud/691',ua(kr),ccbox(kr),
1430	$ ka(kr),beta,mu(im),kr,im,i,nl
1431
1432	call findw(ig,iirw, idummy,c1,p1,Desp,wsL)
1433	if ( i_supersat .eq. 0 ) then
1434	eqw=eqw+no(kr)*w
1435	else
1436	eqw = w + (no(kr)-1) * ( asat_box*dist(kr) )
1437	endif
1438	end do
1439	eqwmu = eqwmu + eqw * mu(im)*amu(im)
1440	end do
1441
1442	tau(in,ir) = exp( - eqwmu / dble(deltanux) )
1443
1444	4 continue
1445
1446	c
1447	c due to the simmetry of the transmittances
1448	c
1449	do in=nl-1,2,-1
1450	do ir=in-1,1,-1
1451	tau(in,ir) = tau(ir,in)
1452	end do
1453	end do
1454
1455
1456	ccc
1457	ccc writing out transmittances
1458	ccc
1459	if (itauout.eq.1) then
1460
1461	! if(ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5
1462	! @ .or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9) then
1463	! open( 1, file=
1464	! @ dircurtis//'taul'//isotcode//dn//ibcode1//'.dat',
1465	! @ access='sequential', form='unformatted' )
1466	! else
1467	! open( 1, file=
1468	! @ dircurtis//'taul'//isotcode//dn//ibcode2//'.dat',
1469	! @ access='sequential', form='unformatted' )
1470	! endif
1471
1472	! write(1) dummy
1473	! write(1)' format: (tauinf(n),(tau(n,r),r=1,nl),n=1,nl)'
1474	! do in=1,nl
1475	! write (1) tauinf(in), ( tau(in,ir), ir=1,nl )
1476	! end do
1477	! close(unit=1)
1478
1479	elseif (itauout.eq.2) then
1480
1481	! if(ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5
1482	! @ .or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9) then
1483	! open( 1, file=
1484	! @ dircurtis//'taul'//isotcode//dn//ibcode1//'.dat')
1485	! else
1486	! open( 1, file=
1487	! @ dircurtis//'taul'//isotcode//dn//ibcode2//'.dat')
1488	! endif
1489
1490	! !write(1,*) dummy
1491	! !write(1,*) 'tij for curtis matrix calculations '
1492	! !write(1,*)' cira mars model atmosphere '
1493	! !write(1,*)' beta= ',beta,'deltanu= ',deltanux
1494	! write(1,*) nl
1495	! write(1,*)
1496	! @ ' format: (tauinf(in),(tau(in,ir),ir=1,nl),in=1,nl)'
1497
1498	! do in=1,nl
1499	! write (1,*) tauinf(in)
1500	! write (1,*) (tau(in,ir), ir=1,nl)
1501	! end do
1502	! close(unit=1)
1503
1504	! if(ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5
1505	! @ .or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9) then
1506	! write (*,'(1x, 31htransmitances written out in: ,a22)')
1507	! @ 'taul'//isotcode//dn//ibcode1
1508	! else
1509	! write (*,'(1x, 31htransmitances written out in: ,a22)')
1510	! @ 'taul'//isotcode//dn//ibcode2
1511	! endif
1512
1513	end if
1514
1515	c cleaning of transmittances
1516	! call elimin_tau(tau,tauinf,nl,nan,itableout,nw,dummy,
1517	! @ isotcode,dn,ibcode2)
1518
1519	c construction of the curtis matrix
1520
1521	call mzcud ( tauinf,tau, cf,cfup,cfdw, vc,taugr,
1522	@ ib,isot,icfout,itableout )
1523
1524	c end
1525	return
1526	end
1527
1528
1529
1530
1531
1532	c***********************************************************************
1533	c mzcud.f
1534	c***********************************************************************
1535
1536	subroutine mzcud( tauinf,tau, c,cup,cdw,vc,taugr,
1537	@ ib,isot,icfout,itableout )
1538
1539	c old times mlp first version of mzcf
1540	c a.k.murphy method to avoid extrapolation in the curtis matrix
1541	c feb-89 malv AKM method to avoid extrapolation in C.M.
1542	c 25-sept-96 cristina dejar las matrices en doble precision
1543	c jan 98 malv version para mz1d
1544	c oct 01 malv update version for fluxes for hb and fb
1545	c jul 2011 malv+fgg Adapted to LMD-MGCM
1546	c***********************************************************************
1547
1548	implicit none
1549
1550	include 'comcstfi.h'
1551	include 'nlte_paramdef.h'
1552	include 'nlte_commons.h'
1553
1554	c arguments
1555	real*8 c(nl,nl), cup(nl,nl), cdw(nl,nl) ! o
1556	real*8 vc(nl), taugr(nl) ! o
1557	real*8 tau(nl,nl) ! i
1558	real*8 tauinf(nl) ! i
1559	integer ib ! i
1560	integer isot ! i
1561	integer icfout, itableout ! i
1562
1563	c external
1564	external bandid
1565	character*2 bandid
1566
1567	c local variables
1568	integer i, in, ir, iw, itblout
1569	real*8 cfup(nl,nl), cfdw(nl,nl)
1570	real*8 a(nl,nl), cf(nl,nl)
1571	character isotcode2, bcode2
1572
1573	c formats
1574	101 format(i1)
1575	202 format(i2)
1576	180 format(a80)
1577	181 format(a80)
1578	c***********************************************************************
1579
1580	if (isot.eq.1) isotcode = '26'
1581	if (isot.eq.2) isotcode = '28'
1582	if (isot.eq.3) isotcode = '36'
1583	if (isot.eq.4) isotcode = '27'
1584	if (isot.eq.5) isotcode = 'co'
1585	bcode = bandid( ib )
1586
1587	! write (,) ' '
1588
1589	do in=1,nl
1590
1591	do ir=1,nl
1592
1593	cf(in,ir) = 0.0d0
1594	cfup(in,ir) = 0.0d0
1595	cfdw(in,ir) = 0.0d0
1596	c(in,ir) = 0.0d0
1597	cup(in,ir) = 0.0d0
1598	cdw(in,ir) = 0.0d0
1599	a(in,ir) = 0.0d0
1600
1601	end do
1602
1603	vc(in) = 0.0d0
1604	taugr(in) = 0.0d0
1605
1606	end do
1607
1608
1609	c the next lines are a reduced and equivalent way of calculating
1610	c the c(in,ir) elements for n=2,nl1 and r=1,nl
1611
1612
1613	c do in=2,nl1
1614	c do ir=1,nl
1615	c if(ir.eq.1)then
1616	c c(in,ir)=tau(in-1,1)-tau(in+1,1)
1617	c elseif(ir.eq.nl)then
1618	c c(in,ir)=tau(in+1,nl1)-tauinf(in+1)-tau(in-1,nl1)+tauinf(in-1)
1619	c else
1620	c c(in,ir)=tau(in+1,ir-1)-tau(in+1,ir)-tau(in-1,ir-1)+tau(in-1,ir)
1621	c end if
1622	c c(in,ir)=c(in,ir)pideltanu(ib)/(2.deltaz1.0e5)
1623	c end do
1624	c end do
1625	c go to 1000
1626
1627	c calculation of the matrix cfup(nl,nl)
1628
1629	cfup(1,1) = 1.d0 - tau(1,1)
1630
1631	do in=2,nl
1632	do ir=1,in
1633
1634	if (ir.eq.1) then
1635	cfup(in,ir) = tau(in,ir) - tau(in,1)
1636	elseif (ir.eq.in) then
1637	cfup(in,ir) = 1.d0 - tau(in,ir-1)
1638	else
1639	cfup(in,ir) = tau(in,ir) - tau(in,ir-1)
1640	end if
1641
1642	end do
1643	end do
1644
1645	! contribution to upwards fluxes from bb at bottom :
1646	do in=1,nl
1647	taugr(in) = tau(in,1)
1648	enddo
1649
1650	c calculation of the matrix cfdw(nl,nl)
1651
1652	cfdw(nl,nl) = 1.d0 - tauinf(nl)
1653
1654	do in=1,nl-1
1655	do ir=in,nl
1656
1657	if (ir.eq.in) then
1658	cfdw(in,ir) = 1.d0 - tau(in,ir)
1659	elseif (ir.eq.nl) then
1660	cfdw(in,ir) = tau(in,ir-1) - tauinf(in)
1661	else
1662	cfdw(in,ir) = tau(in,ir-1) - tau(in,ir)
1663	end if
1664
1665	end do
1666	end do
1667
1668
1669	c calculation of the matrix cf(nl,nl)
1670
1671	do in=1,nl
1672	do ir=1,nl
1673
1674	if (ir.eq.1) then
1675	! version con l_bb(tg) = l_bb(t(1))=j(1) (see also vc below)
1676	! cf(in,ir) = tau(in,ir)
1677	! version con l_bb(tg) =/= l_bb(t(1))=j(1) (see also vc below)
1678	cf(in,ir) = tau(in,ir) - tau(in,1)
1679	elseif (ir.eq.nl) then
1680	cf(in,ir) = tauinf(in) - tau(in,ir-1)
1681	else
1682	cf(in,ir) = tau(in,ir) - tau(in,ir-1)
1683	end if
1684
1685	end do
1686	end do
1687
1688
1689	c definition of the a(nl,nl) matrix
1690
1691	do in=2,nl-1
1692	do ir=1,nl
1693	if (ir.eq.in+1) a(in,ir) = -1.d0
1694	if (ir.eq.in-1) a(in,ir) = +1.d0
1695	a(in,ir) = a(in,ir) / ( 2.d0deltaz1.d5 )
1696	end do
1697	end do
1698	! this is not needed anymore in the akm scheme
1699	! a(1,1) = +3.d0
1700	! a(1,2) = -4.d0
1701	! a(1,3) = +1.d0
1702	! a(nl,nl) = -3.d0
1703	! a(nl,nl1) = +4.d0
1704	! a(nl,nl2) = -1.d0
1705
1706	c calculation of the final curtis matrix ("reduced" by murphy's method)
1707
1708	if (isot.ne.5) then
1709	do in=1,nl
1710	do ir=1,nl
1711	cf(in,ir) = cf(in,ir) * pi*deltanu(isot,ib)
1712	cfup(in,ir) = cfup(in,ir) * pi*deltanu(isot,ib)
1713	cfdw(in,ir) = cfdw(in,ir) * pi*deltanu(isot,ib)
1714	end do
1715	taugr(in) = taugr(in) * pi*deltanu(isot,ib)
1716	end do
1717	else
1718	do in=1,nl
1719	do ir=1,nl
1720	cf(in,ir) = cf(in,ir) * pi*deltanuco
1721	enddo
1722	taugr(in) = taugr(in) * pi*deltanuco
1723	enddo
1724	endif
1725
1726	do in=2,nl-1
1727
1728	do ir=1,nl
1729
1730	do i=1,nl
1731	! only c contains the matrix a. matrixes cup,cdw dont because
1732	! these two will be used for flux calculations, not
1733	! only for flux divergencies
1734
1735	c(in,ir) = c(in,ir) + a(in,i) * cf(i,ir)
1736	! from this matrix we will extract (see below) the
1737	! nl2 x nl2 "core" for the "reduced" final curtis matrix.
1738
1739	end do
1740	cup(in,ir) = cfup(in,ir)
1741	cdw(in,ir) = cfdw(in,ir)
1742
1743	end do
1744	! version con l_bb(tg) = l_bb(t(1))=j(1) (see cf above)
1745	!vc(in) = c(in,1)
1746	! version con l_bb(tg) =/= l_bb(t(1))=j(1) (see cf above)
1747	if (isot.ne.5) then
1748	vc(in) = pideltanu(isot,ib)/( 2.d0deltaz1.d5 )
1749	@ ( tau(in-1,1) - tau(in+1,1) )
1750	else
1751	vc(in) = pideltanuco/( 2.d0deltaz1.d5 )
1752	@ ( tau(in-1,1) - tau(in+1,1) )
1753	endif
1754
1755	end do
1756
1757	5 continue
1758
1759	! write (,) 'mztf/1/ c(2,*) =', (c(2,i), i=1,nl)
1760
1761	! call elimin_dibuja(c,nl,itableout)
1762
1763	c ventana del smoothing de c es nw=3 y de vc es 5 (puesto en lisa):
1764	c subroutine elimin_mz4(c,vc,ilayer,nl,nan,iw, nw)
1765
1766	iw = nan
1767	if (isot.eq.4) iw = 5 ! eliminates values < 1.d-19
1768	if (itableout.eq.30) then
1769	itblout = 0
1770	else
1771	itblout = itableout
1772	endif
1773	call elimin_mz1d (c,vc,0,iw,itblout,nw)
1774
1775	! upper boundary condition
1776	! j'(nl) = j'(nl1) ==> j(nl) = 2j(nl1) - j(nl2) ==>
1777	do in=2,nl-1
1778	c(in,nl-2) = c(in,nl-2) - c(in,nl)
1779	c(in,nl-1) = c(in,nl-1) + 2.d0*c(in,nl)
1780	cup(in,nl-2) = cup(in,nl-2) - cup(in,nl)
1781	cup(in,nl-1) = cup(in,nl-1) + 2.d0*cup(in,nl)
1782	cdw(in,nl-2) = cdw(in,nl-2) - cdw(in,nl)
1783	cdw(in,nl-1) = cdw(in,nl-1) + 2.d0*cdw(in,nl)
1784	end do
1785	! j(nl) = j(nl1) ==>
1786	! do in=2,nl1
1787	! c(in,nl1) = c(in,nl1) + c(in,nl)
1788	! end do
1789
1790	! 1000 continue
1791
1792
1793	if (icfout.eq.1) then
1794
1795	! if (ib.eq.1 .or. ib.eq.12 .or. ib.eq.16 .or. ib.eq.18) then
1796	! codmatrx = codmatrx_fb
1797	! else
1798	! codmatrx = codmatrx_hot
1799	! end if
1800	! if (ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5
1801	! @ .or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9) then
1802	! ibcode2 = '0'//ibcode1
1803	! else
1804	! write ( ibcode2, 202) ib
1805	! endif
1806
1807	! open ( 1, access='sequential', form='unformatted', file=
1808	! @ dircurtis//'cfl'//isotcode//dn//ibcode2//codmatrx//'.dat')
1809	! open ( 2, access='sequential', form='unformatted', file=
1810	! @ dircurtis//'cflup'//isotcode//dn//ibcode2//codmatrx//'.dat')
1811	! open ( 3, access='sequential', form='unformatted', file=
1812	! @ dircurtis//'cfldw'//isotcode//dn//ibcode2//codmatrx//'.dat')
1813	! open ( 4, access='sequential', form='unformatted', file=
1814	! @ dircurtis//'cflgr'//isotcode//dn//ibcode2//codmatrx//'.dat')
1815
1816	! write(1) dummy
1817	! write(1) ' format: (vc(n),(ch(n,r),r=2,nl-1),n=2,nl-1)'
1818	! do in=2,nl-1
1819	! write(1) vc(in), (c(in,ir) , ir=2,nl-1 )
1820	!! write (,) in, vc(in)
1821	! end do
1822
1823	! write(2) dummy
1824	! write(2)' format: (cfup(n,r),r=1,nl), n=1,nl)'
1825	! do in=1,nl
1826	! write(2) ( cup(in,ir) , ir=1,nl )
1827	! end do
1828
1829	! write(3) dummy
1830	! write(3)' format: (cfdw(n,r),r=1,nl), n=1,nl)'
1831	! do in=1,nl
1832	! write(3) (cdw(in,ir) , ir=1,nl )
1833	! end do
1834
1835	! write(4) dummy
1836	! write(4)' format: (taugr(n), n=1,nl)'
1837	! do in=1,nl
1838	! write(4) (taugr(in), ir=1,nl )
1839	! end do
1840	! !write (,) ' Last value in file: ', taugr(nl)
1841
1842	! write (*,'(1x,30hcurtis matrix written out in: ,a,a,a,a)' )
1843	! @ dircurtis//'cfl'//isotcode//dn//ibcode2//codmatrx//'.dat',
1844	! @ dircurtis//'cflup'//isotcode//dn//ibcode2//codmatrx//'.dat',
1845	! @ dircurtis//'cfldw'//isotcode//dn//ibcode2//codmatrx//'.dat',
1846	! @ dircurtis//'cflgr'//isotcode//dn//ibcode2//codmatrx//'.dat'
1847
1848	! close (1)
1849	! close (2)
1850	! close (3)
1851	! close (4)
1852
1853	else
1854
1855	! write (,) ' no curtis matrix output file ', char(10)
1856
1857	end if
1858
1859	if (itableout.eq.30) then ! Force output of C.M. in ascii file
1860
1861	! if (ib.eq.1 .or. ib.eq.12 .or. ib.eq.16 .or. ib.eq.18) then
1862	! codmatrx = codmatrx_fb
1863	! else
1864	! codmatrx = codmatrx_hot
1865	! end if
1866	! if (ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5
1867	! @ .or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9) then
1868	! ibcode2 = '0'//ibcode1
1869	! else
1870	! write ( ibcode2, 202) ib
1871	! endif
1872
1873	! open (10, file=
1874	! & dircurtis//'table'//isotcode//dn//ibcode2//codmatrx//'.dat')
1875	! write(10,*) nl, ' = number of layers '
1876	! write(10,*) ' format: (vc(n),(ch(n,r),r=2,nl-1),n=2,nl-1)'
1877	! do in=2,nl-1
1878	! write(10,*) vc(in), (c(in,ir) , ir=2,nl-1 )
1879	! enddo
1880	! close (10)
1881	endif
1882
1883	c end
1884	return
1885	end
1886
1887
1888
1889
1890
1891	c***********************************************************************
1892	c mztvc
1893	c***********************************************************************
1894
1895	subroutine mztvc ( ig,vc, ib,isot,
1896	@ iirw,iimu,itauout,icfout,itableout )
1897
1898	c jul 2011 malv+fgg
1899	c***********************************************************************
1900
1901	implicit none
1902
1903	include 'comcstfi.h'
1904	include 'nlte_paramdef.h'
1905	include 'nlte_commons.h'
1906
1907	c arguments
1908	integer ig ! ADDED FOR TRACEBACK
1909	real*8 cf(nl,nl), cfup(nl,nl), cfdw(nl,nl) ! o
1910	real*8 vc(nl), taugr(nl) ! o
1911	integer ib ! i
1912	integer isot ! i
1913	integer iirw ! i
1914	integer iimu ! i
1915	integer itauout ! i
1916	integer icfout ! i
1917	integer itableout ! i
1918
1919	c local variables and constants
1920	integer i, in, ir, im, k ,j
1921	integer nmu
1922	parameter (nmu = 8)
1923	real*8 tau(nl,nl)
1924	real*8 tauinf(nl)
1925	real*8 con(nzy), coninf
1926	real*8 c1, c2
1927	real*8 t1, t2
1928	real*8 p1, p2
1929	real*8 mr1, mr2
1930	real*8 st1, st2
1931	real*8 c1box(70), c2box(70)
1932	real*8 ff ! to avoid too small numbers
1933	real*8 tvtbs(nzy)
1934	real*8 st, beta, ts, eqwmu
1935	real*8 mu(nmu), amu(nmu)
1936	real*8 zld(nl), zyd(nzy)
1937	real*8 correc
1938	real deltanux ! width of vib-rot band (cm-1)
1939	character isotcode*2
1940	integer idummy
1941	real*8 Desp,wsL
1942
1943	c formats
1944	111 format(a1)
1945	112 format(a2)
1946	101 format(i1)
1947	202 format(i2)
1948	180 format(a80)
1949	181 format(a80)
1950	c***********************************************************************
1951
1952	c some needed values
1953	! rl=sqrt(log(2.d0))
1954	! pi2 = 3.14159265358989d0
1955	beta = 1.8d0
1956	! beta = 1.0d0
1957	idummy = 0
1958	Desp = 0.0d0
1959	wsL = 0.0d0
1960
1961	!write (,) ' MZTUD/ iirw = ', iirw
1962
1963
1964	c esto es para que las subroutines de mztfsub calculen we
1965	c de la forma apropiada para mztf, no para fot
1966	icls=icls_mztf
1967
1968	c codigos para filenames
1969	! if (isot .eq. 1) isotcode = '26'
1970	! if (isot .eq. 2) isotcode = '28'
1971	! if (isot .eq. 3) isotcode = '36'
1972	! if (isot .eq. 4) isotcode = '27'
1973	! if (isot .eq. 5) isotcode = '62'
1974	! if(ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5
1975	! @ .or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9) then
1976	! write (ibcode1,101) ib
1977	! else
1978	! write (ibcode2,202) ib
1979	! endif
1980	! write (*,'( 30h calculating curtis matrix : ,2x,
1981	! @ 8h band = ,i2,2x, 11h isotope = ,i2)') ib, isot
1982
1983	c integration in angle !!!!!!!!!!!!!!!!!!!!
1984
1985	c------- diffusivity approx.
1986	if (iimu.eq.1) then
1987	! write (,) ' diffusivity approx. beta = ',beta
1988	mu(1) = 1.0d0
1989	amu(1)= 1.0d0
1990	c-------data for 8 points integration
1991	elseif (iimu.eq.4) then
1992	write (,)' 4 points for the gauss-legendre angle quadrature.'
1993	mu(1)=(1.0d0+0.339981043584856)/2.0d0
1994	mu(2)=(1.0d0-0.339981043584856)/2.0d0
1995	mu(3)=(1.0d0+0.861136311594053)/2.0d0
1996	mu(4)=(1.0d0-0.861136311594053)/2.0d0
1997	amu(1)=0.652145154862546
1998	amu(2)=amu(1)
1999	amu(3)=0.347854845137454
2000	amu(4)=amu(3)
2001	beta=1.0d0
2002	c-------data for 8 points integration
2003	elseif(iimu.eq.8) then
2004	write (,)' 8 points for the gauss-legendre angle quadrature.'
2005	mu(1)=(1.0d0+0.183434642495650)/2.0d0
2006	mu(2)=(1.0d0-0.183434642495650)/2.0d0
2007	mu(3)=(1.0d0+0.525532409916329)/2.0d0
2008	mu(4)=(1.0d0-0.525532409916329)/2.0d0
2009	mu(5)=(1.0d0+0.796666477413627)/2.0d0
2010	mu(6)=(1.0d0-0.796666477413627)/2.0d0
2011	mu(7)=(1.0d0+0.960289856497536)/2.0d0
2012	mu(8)=(1.0d0-0.960289856497536)/2.0d0
2013	amu(1)=0.362683783378362
2014	amu(2)=amu(1)
2015	amu(3)=0.313706645877887
2016	amu(4)=amu(3)
2017	amu(5)=0.222381034453374
2018	amu(6)=amu(5)
2019	amu(7)=0.101228536290376
2020	amu(8)=amu(7)
2021	beta=1.0d0
2022	end if
2023	c!!!!!!!!!!!!!!!!!!!!!!!
2024
2025	ccc
2026	ccc determine abundances included in the absorber amount
2027	ccc
2028
2029	c first, set up the grid ready for interpolation.
2030	do i=1,nzy
2031	zyd(i) = dble(zy(i))
2032	enddo
2033	do i=1,nl
2034	zld(i) = dble(zl(i))
2035	enddo
2036
2037	c vibr. temp of the bending mode :
2038	if (isot.eq.1) call interdp ( tvtbs,zyd,nzy, v626t1,zld,nl, 1 )
2039	if (isot.eq.2) call interdp ( tvtbs,zyd,nzy, v628t1,zld,nl, 1 )
2040	if (isot.eq.3) call interdp ( tvtbs,zyd,nzy, v636t1,zld,nl, 1 )
2041	if (isot.eq.4) call interdp ( tvtbs,zyd,nzy, v627t1,zld,nl, 1 )
2042	!if (isot.eq.5) call interdp ( tvtbs,zxd,nz, vcot1,zld,nl, 1 )
2043
2044	c 2nd: correccion a la n10(i) (cantidad de absorbente en el lower state)
2045	c por similitud a la que se hace en cza.for ; esto solo se hace para CO2
2046
2047	!write (,) 'imr(isot) = ', isot, imr(isot)
2048	do i=1,nzy
2049	if (isot.eq.5) then
2050	con(i) = dble( coy(i) * imrco )
2051	else
2052	con(i) = dble( co2y(i) * imr(isot) )
2053	correc = 2.d0 * dexp( dble(-ee*elow(isot,2))/tvtbs(i) )
2054	con(i) = con(i) * ( 1.d0 - correc )
2055	! write (,) ' iz, correc, co2y(i), con(i) =',
2056	! @ i,correc,co2y(i),con(i)
2057	endif
2058
2059	!-----------------------------------------------------------------
2060	! mlp & cristina. 17 july 1996 change the calculation of mr.
2061	! it is used for calculating partial press
2062	! alpha = alpha(self,co2)pp +alpha(n2)(pt-pp)
2063	! for an isotope, if mr is obtained by
2064	! co2*imr(iso)/nt
2065	! we are considerin collisions with other co2 isotopes
2066	! (including the major one, 626) as if they were with n2.
2067	! assuming mr as co2/nt, we consider collisions
2068	! of type 628-626 as of 626-626 instead of as 626-n2.
2069	! mrx(i)=con(i)/ntx(i) ! old malv
2070	! mrx(i)= dble(co2x(i)/ntx(i)) ! mlp & crs
2071
2072	! jan 98:
2073	! esta modif de mlp implica anular el correc (deberia revisar esto)
2074
2075	mr(i) = dble(co2y(i)/nty(i)) ! malv, jan 98
2076
2077	!-----------------------------------------------------------------
2078
2079	end do
2080
2081	! como beta y 1.d5 son comunes a todas las weighted absorber amounts,
2082	! los simplificamos:
2083	! coninf = beta * 1.d5 * dble( con(n) / log( con(n-1) / con(n) ) )
2084	!write (,) ' con(nz), con(nz-1) =', con(nz), con(nz-1)
2085	coninf = dble( con(nzy) / log( con(nzy-1) / con(nzy) ) )
2086	!write (,) ' coninf =', coninf
2087
2088	ccc
2089	ccc temp dependence of the band strength and
2090	ccc nlte correction factor for the absorber amount
2091	ccc
2092	call mztf_correccion ( coninf, con, ib, isot, itableout )
2093
2094	ccc
2095	ccc reads histogrammed spectral data (strength for lte and vmr=1)
2096	ccc
2097	!hfile1 = dirspec//'hi'//dn !Ya no hacemos distincion d/n en esto
2098	!! hfile1 = dirspec//'hid' !(see why in his.for)
2099	! hfile1='hid'
2100	!! if (ib.eq.13 .or. ib.eq.14 ) hfile1 = dirspec//'his'
2101	! if (ib.eq.13 .or. ib.eq.14 ) hfile1 = 'his'
2102
2103	! if(ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5
2104	! @ .or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9) then
2105	! if (isot.eq.1) hisfile = hfile1//'26-'//ibcode1//'.dat'
2106	! if (isot.eq.2) hisfile = hfile1//'28-'//ibcode1//'.dat'
2107	! if (isot.eq.3) hisfile = hfile1//'36-'//ibcode1//'.dat'
2108	! if (isot.eq.4) hisfile = hfile1//'27-'//ibcode1//'.dat'
2109	! if (isot.eq.5) hisfile = hfile1//'62-'//ibcode1//'.dat'
2110	! else
2111	! if (isot.eq.1) hisfile = hfile1//'26-'//ibcode2//'.dat'
2112	! if (isot.eq.2) hisfile = hfile1//'28-'//ibcode2//'.dat'
2113	! if (isot.eq.3) hisfile = hfile1//'36-'//ibcode2//'.dat'
2114	! if (isot.eq.4) hisfile = hfile1//'27-'//ibcode2//'.dat'
2115	! if (isot.eq.5) hisfile = hfile1//'62-'//ibcode2//'.dat'
2116	! endif
2117	! write (,) 'hisfile: ', hisfile
2118
2119	! the argument to rhist is to make this compatible with mztf_comp.f,
2120	! which is a useful modification of mztf.f (to change strengths of bands
2121	! call rhist (1.0)
2122	if(ib.eq.1) then
2123	if(isot.eq.1) then !Case 1
2124	mm=mm_c1
2125	nbox=nbox_c1
2126	tmin=tmin_c1
2127	tmax=tmax_c1
2128	do i=1,nbox_max
2129	no(i)=no_c1(i)
2130	dist(i)=dist_c1(i)
2131	do j=1,nhist
2132	sk1(j,i)=sk1_c1(j,i)
2133	xls1(j,i)=xls1_c1(j,i)
2134	xln1(j,i)=xln1_c1(j,i)
2135	xld1(j,i)=xld1_c1(j,i)
2136	enddo
2137	enddo
2138	do j=1,nhist
2139	thist(j)=thist_c1(j)
2140	enddo
2141	else if(isot.eq.2) then !Case 2
2142	mm=mm_c2
2143	nbox=nbox_c2
2144	tmin=tmin_c2
2145	tmax=tmax_c2
2146	do i=1,nbox_max
2147	no(i)=no_c2(i)
2148	dist(i)=dist_c2(i)
2149	do j=1,nhist
2150	sk1(j,i)=sk1_c2(j,i)
2151	xls1(j,i)=xls1_c2(j,i)
2152	xln1(j,i)=xln1_c2(j,i)
2153	xld1(j,i)=xld1_c2(j,i)
2154	enddo
2155	enddo
2156	do j=1,nhist
2157	thist(j)=thist_c2(j)
2158	enddo
2159	else if(isot.eq.3) then !Case 3
2160	mm=mm_c3
2161	nbox=nbox_c3
2162	tmin=tmin_c3
2163	tmax=tmax_c3
2164	do i=1,nbox_max
2165	no(i)=no_c3(i)
2166	dist(i)=dist_c3(i)
2167	do j=1,nhist
2168	sk1(j,i)=sk1_c3(j,i)
2169	xls1(j,i)=xls1_c3(j,i)
2170	xln1(j,i)=xln1_c3(j,i)
2171	xld1(j,i)=xld1_c3(j,i)
2172	enddo
2173	enddo
2174	do j=1,nhist
2175	thist(j)=thist_c3(j)
2176	enddo
2177	else if(isot.eq.4) then !Case 4
2178	mm=mm_c4
2179	nbox=nbox_c4
2180	tmin=tmin_c4
2181	tmax=tmax_c4
2182	do i=1,nbox_max
2183	no(i)=no_c4(i)
2184	dist(i)=dist_c4(i)
2185	do j=1,nhist
2186	sk1(j,i)=sk1_c4(j,i)
2187	xls1(j,i)=xls1_c4(j,i)
2188	xln1(j,i)=xln1_c4(j,i)
2189	xld1(j,i)=xld1_c4(j,i)
2190	enddo
2191	enddo
2192	do j=1,nhist
2193	thist(j)=thist_c4(j)
2194	enddo
2195	else
2196	write(,)'isot must be 2,3 or 4 for ib=1!!'
2197	write(,)'stop at mztvc/310'
2198	stop
2199	endif
2200	else if (ib.eq.2) then
2201	if(isot.eq.1) then !Case 5
2202	mm=mm_c5
2203	nbox=nbox_c5
2204	tmin=tmin_c5
2205	tmax=tmax_c5
2206	do i=1,nbox_max
2207	no(i)=no_c5(i)
2208	dist(i)=dist_c5(i)
2209	do j=1,nhist
2210	sk1(j,i)=sk1_c5(j,i)
2211	xls1(j,i)=xls1_c5(j,i)
2212	xln1(j,i)=xln1_c5(j,i)
2213	xld1(j,i)=xld1_c5(j,i)
2214	enddo
2215	enddo
2216	do j=1,nhist
2217	thist(j)=thist_c5(j)
2218	enddo
2219	else
2220	write(,)'isot must be 1 for ib=2!!'
2221	write(,)'stop at mztvc/334'
2222	stop
2223	endif
2224	else if (ib.eq.3) then
2225	if(isot.eq.1) then !Case 6
2226	mm=mm_c6
2227	nbox=nbox_c6
2228	tmin=tmin_c6
2229	tmax=tmax_c6
2230	do i=1,nbox_max
2231	no(i)=no_c6(i)
2232	dist(i)=dist_c6(i)
2233	do j=1,nhist
2234	sk1(j,i)=sk1_c6(j,i)
2235	xls1(j,i)=xls1_c6(j,i)
2236	xln1(j,i)=xln1_c6(j,i)
2237	xld1(j,i)=xld1_c6(j,i)
2238	enddo
2239	enddo
2240	do j=1,nhist
2241	thist(j)=thist_c6(j)
2242	enddo
2243	else
2244	write(,)'isot must be 1 for ib=3!!'
2245	write(,)'stop at mztvc/358'
2246	stop
2247	endif
2248	else if (ib.eq.4) then
2249	if(isot.eq.1) then !Case 7
2250	mm=mm_c7
2251	nbox=nbox_c7
2252	tmin=tmin_c7
2253	tmax=tmax_c7
2254	do i=1,nbox_max
2255	no(i)=no_c7(i)
2256	dist(i)=dist_c7(i)
2257	do j=1,nhist
2258	sk1(j,i)=sk1_c7(j,i)
2259	xls1(j,i)=xls1_c7(j,i)
2260	xln1(j,i)=xln1_c7(j,i)
2261	xld1(j,i)=xld1_c7(j,i)
2262	enddo
2263	enddo
2264	do j=1,nhist
2265	thist(j)=thist_c7(j)
2266	enddo
2267	else
2268	write(,)'isot must be 1 for ib=4!!'
2269	write(,)'stop at mztvc/382'
2270	stop
2271	endif
2272	else
2273	write(,)'ib must be 1,2,3 or 4!!'
2274	write(,)'stop at mztvc/387'
2275	endif
2276
2277
2278	c******
2279	c****** calculation of tau(1,ir) for 1<=r
2280	c******
2281	call initial
2282
2283	ff=1.0e10
2284
2285	in=1
2286
2287	tau(in,1) = 1.d0
2288
2289	call initial
2290	call intz (zl(in), c1,p1,mr1,t1, con)
2291	do kr=1,nbox
2292	ta(kr) = t1
2293	end do
2294	call interstrength (st1,t1,ka,ta)
2295	do kr=1,nbox
2296	c1box(kr) = c1 * ka(kr) * dble(deltaz)
2297	end do
2298	c1 = c1 * st1 * dble(deltaz)
2299
2300	do 2 ir=2,nl
2301
2302	call intz (zl(ir), c2,p2,mr2,t2, con)
2303	do kr=1,nbox
2304	ta(kr) = t2
2305	end do
2306	call interstrength (st2,t2,ka,ta)
2307	do kr=1,nbox
2308	c2box(kr) = c2 * ka(kr) * dble(deltaz)
2309	end do
2310	c2 = c2 * st2 * dble(deltaz)
2311
2312	aa = aa + ( p1mr1(c1ff) + p2mr2(c2ff)) / 2.d0
2313	bb = bb + ( p1c1 + p2c2 ) / 2.d0
2314	cc = cc + ( c1 + c2 ) / 2.d0
2315	dd = dd + ( t1c1 + t2c2 ) / 2.d0
2316	do kr=1,nbox
2317	ccbox(kr) = ccbox(kr) + ( c1box(kr) + c2box(kr) ) /2.d0
2318	ddbox(kr) = ddbox(kr) +
2319	$ ( t1c1box(kr) + t2c2box(kr) ) /2.d0
2320	end do
2321
2322	mr1=mr2
2323	t1=t2
2324	c1=c2
2325	p1=p2
2326	do kr=1,nbox
2327	c1box(kr) = c2box(kr)
2328	end do
2329
2330	pt = bb / cc
2331	pp = aa / (cc * ff)
2332
2333	ts = dd/cc
2334	do kr=1,nbox
2335	ta(kr) = ddbox(kr) / ccbox(kr)
2336	end do
2337	call interstrength(st,ts,ka,ta)
2338	call intershape(alsa,alna,alda,ta)
2339
2340
2341	eqwmu = 0.0d0
2342	do im = 1,iimu
2343	eqw=0.0d0
2344	do kr=1,nbox
2345	ua(kr) = ccbox(kr) / ka(kr) * beta * 1.0d5 / mu(im)
2346	call findw(ig,iirw, idummy,c1,p1,Desp,wsL)
2347	if ( i_supersat .eq. 0 ) then
2348	eqw=eqw+no(kr)*w
2349	else
2350	eqw = w + (no(kr)-1) * ( asat_box*dist(kr) )
2351	endif
2352	end do
2353	eqwmu = eqwmu + eqw * mu(im)*amu(im)
2354	end do
2355
2356	tau(in,ir) = exp( - eqwmu / dble(deltanu(isot,ib)) )
2357
2358	2 continue
2359
2360
2361
2362	c
2363	c due to the simmetry of the transmittances
2364	c
2365	do in=nl,2,-1
2366	tau(in,1) = tau(1,in)
2367	end do
2368
2369	vc(1) = 0.0d0
2370	vc(nl) = 0.0d0
2371	do in=2,nl-1 ! poner aqui nl-1 luego
2372	vc(in) = pideltanu(isot,ib)/( 2.d0deltaz1.d5 )
2373	@ ( tau(in-1,1) - tau(in+1,1) )
2374	end do
2375
2376
2377	c end
2378	return
2379	end
2380
2381
2382
2383
2384
2385	c***********************************************************************
2386	c mztvc_626fh.F
2387	c***********************************************************************
2388
2389	subroutine mztvc_626fh(ig)
2390
2391	c jul 2011 malv+fgg
2392	c***********************************************************************
2393
2394	implicit none
2395
2396	!!!!!!!!!!!!!!!!!!!!!!!
2397	! common variables & constants
2398
2399	include 'nlte_paramdef.h'
2400	include 'nlte_commons.h'
2401
2402	!!!!!!!!!!!!!!!!!!!!!!!
2403	! arguments
2404
2405	integer ig ! ADDED FOR TRACEBACK
2406
2407	!!!!!!!!!!!!!!!!!!!!!!!
2408	! local variables
2409
2410	real*4 cdummy(nl,nl), csngl(nl,nl)
2411
2412	real*8 cax1(nl,nl), cax2(nl,nl), cax3(nl,nl)
2413	real*8 v1(nl), v2(nl), v3(nl), cm_factor, vc_factor
2414
2415	integer itauout,icfout,itableout, interpol,ismooth, isngldble
2416	integer i,j,ik,ist,isot,ib,itt
2417
2418	!character bandcode*2
2419	character isotcode*2
2420	!character codmatrx_hot*5
2421
2422	!!!!!!!!!!!!!!!!!!!!!!!
2423	! external functions
2424
2425	external bandid
2426	character*2 bandid
2427
2428	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2429	! subroutines called:
2430	! mz4sub, dmzout, readc_mz4, readcupdw, mztf
2431
2432	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2433	! formatos
2434	132 format(i2)
2435
2436	************************************************************************
2437	************************************************************************
2438
2439	isngldble = 1 ! =1 --> dble precission
2440
2441	fileroot = 'cfl'
2442
2443	ist = 1
2444	isot = 26
2445	write (isotcode,132) isot
2446
2447	call zerov( vc121, nl )
2448
2449	do 11, ik=1,3
2450
2451	ib=ik+1
2452
2453	call mztvc (ig,v1, ib, 1, irw_mztf, imu, 0,0,0 )
2454
2455	do i=1,nl
2456
2457	if(ik.eq.1)then
2458	vc_factor = dble(667.75/618.03)
2459	elseif(ik.eq.2)then
2460	vc_factor = 1.d0
2461	elseif(ik.eq.3)then
2462	vc_factor = dble(667.75/720.806)
2463	end if
2464
2465	vc121(i) = vc121(i) + v1(i) * vc_factor
2466
2467	end do
2468
2469	11 continue
2470
2471
2472	return
2473	end
2474
2475
2476
2477
2478
2479	c***********************************************************************
2480	c mztf_correccion
2481	c***********************************************************************
2482
2483	subroutine mztf_correccion (coninf, con, ib, isot, icurt_pop)
2484
2485	c including the dependence of the absort. coeff. on temp., vibr. temp.,
2486	c function, etc.., when neccessary. imr is already corrected in his.for
2487	c we follow pg.39b-43a (l5):
2488	c tvt1 is the vibr temp of the upper level
2489	c tvt is the vibr temp of the transition itself
2490	c tvtbs is the vibr temp of the bending mode (used in qv)
2491	c for fundamental bands, they are not used at the moment.
2492	c for the 15 fh and sh bands, only tvt0 is used at the moment.
2493	c for the laser band, all of them are used following pg. 41a -l5- :
2494	c we need s(z) and we can read s(tk) from the histogram (also called
2495	c what we have to calculate now is the factor s(z)/s(tk) or following
2496	c l5 notebook notation, s_nlte/s_lte.
2497	c s_nlte/s_lte = xfactor = xlower * xqv * xes
2498
2499	c icurt_pop = 30 -> Output of populations of the 0200,0220,1000 states
2500	c = otro -> no output of these populations
2501
2502	c oct 92 malv
2503	c jan 98 malv version for mz1d
2504	c jul 2011 malv+fgg adapted to LMD-MGCM
2505	c***********************************************************************
2506
2507	implicit none
2508
2509	include 'nlte_paramdef.h'
2510	include 'nlte_commons.h'
2511
2512	c arguments
2513	integer ib, isot
2514	integer icurt_pop ! output of Fermi states population
2515	real*8 con(nzy), coninf
2516
2517	! local variables
2518	integer i
2519	real*8 tvt0(nzy),tvt1(nzy),tvtbs(nzy), zld(nl),zyd(nzy)
2520	real xalfa, xbeta, xtv1000, xtv0200, xtv0220, xfactor
2521	real xqv, xnu_trans, xtv_trans, xes, xlower
2522	c***********************************************************************
2523
2524	xfactor = 1.0
2525
2526	do i=1,nzy
2527	zyd(i) = dble(zy(i))
2528	enddo
2529	do i=1,nl
2530	zld(i) = dble( zl(i) )
2531	end do
2532
2533	! tvtbs is the bending mode of the molecule. used in xqv.
2534	if (isot.eq.1) call interdp (tvtbs,zyd,nzy, v626t1,zld,nl, 1 )
2535	if (isot.eq.2) call interdp (tvtbs,zyd,nzy, v628t1,zld,nl, 1 )
2536	if (isot.eq.3) call interdp (tvtbs,zyd,nzy, v636t1,zld,nl, 1 )
2537	if (isot.eq.4) call interdp (tvtbs,zyd,nzy, v627t1,zld,nl, 1 )
2538	if (isot.eq.5) call interdp (tvtbs,zyd,nzy, vcot1,zld,nl, 1 )
2539
2540	! tvt0 is the lower level of the transition. used in xlower.
2541	if (ib.eq.2 .or. ib.eq.3 .or. ib.eq.4 .or. ib.eq.15) then
2542	if (isot.eq.1) call interdp (tvt0,zyd,nzy, v626t1,zld,nl, 1 )
2543	if (isot.eq.2) call interdp (tvt0,zyd,nzy, v628t1,zld,nl, 1 )
2544	if (isot.eq.3) call interdp (tvt0,zyd,nzy, v636t1,zld,nl, 1 )
2545	if (isot.eq.4) call interdp (tvt0,zyd,nzy, v627t1,zld,nl, 1 )
2546	elseif (ib.eq.6 .or. ib.eq.8 .or. ib.eq.10
2547	@ .or. ib.eq.13 .or. ib.eq.14
2548	@ .or. ib.eq.17 .or. ib.eq.19 .or. ib.eq.20) then
2549	if (isot.eq.1) call interdp ( tvt0,zyd,nzy, v626t2,zld,nl, 1 )
2550	if (isot.eq.2) call interdp ( tvt0,zyd,nzy, v628t2,zld,nl, 1 )
2551	if (isot.eq.3) call interdp ( tvt0,zyd,nzy, v636t2,zld,nl, 1 )
2552	if (isot.eq.4) then
2553	call interdp ( tvt0,zyd,nzy, v627t2,zld,nl, 1 )
2554	endif
2555	else
2556	do i=1,nzy
2557	tvt0(i) = dble( ty(i) )
2558	end do
2559	end if
2560
2561	c tvt is the vt of the transition. used in xes.
2562	c since xes=1.0 except for the laser bands, tvt is only needed for them
2563	c but it is actually calculated from the tv of the upper and lower level
2564	c of the transition. hence, only tvt1 remains to be read for the laser b
2565	c tvt1 is the upper level of the transition.
2566	if (ib.eq.13 .or. ib.eq.14) then
2567	if (isot.eq.1) call interdp ( tvt1,zyd,nzy, v626t4,zld,nl, 1 )
2568	if (isot.eq.2) call interdp ( tvt1,zyd,nzy, v628t4,zld,nl, 1 )
2569	if (isot.eq.3) call interdp ( tvt1,zyd,nzy, v636t4,zld,nl, 1 )
2570	if (isot.eq.4) call interdp ( tvt1,zyd,nzy, v627t4,zld,nl, 1 )
2571	end if
2572
2573	c here we weight the absorber amount by a factor which compensate the l
2574	c value of the strength read from hitran. we use that factor in order t
2575	c correct the product s*m when we later multiply those two variables.
2576
2577	! if ( isot.eq.1 .and. icurt_pop.eq.30 ) then
2578	! open (30, file='020populations.dat')
2579	! write (30,*) ' z tv(020) tv0200 tv0220 tv1000 '
2580	! endif
2581
2582	do i=1,nzy
2583
2584	if (isot.eq.1) then
2585
2586	!!! vt of the 3 levels in (020) (see pag. 36a-sn1 for this)
2587	xalfa = 1.d0/2.d0* exp( dble(-ee*(nu12_1000-nu(1,2))/ty(i)) )
2588	xbeta = 1.d0/2.d0* exp( dble(-ee*(nu12_0200-nu(1,2))/ty(i)) )
2589	xtv0200 = dble( - ee * nu12_0200 ) /
2590	@ ( log( xbeta/(1.d0+xalfa+xbeta) ) -
2591	@ dble(ee*nu(1,2))/tvt0(i) )
2592	xtv0220 = dble( - ee * nu(1,2) ) /
2593	@ ( log( 1.d0/(1.d0+xalfa+xbeta) ) -
2594	@ dble(ee*nu(1,2))/tvt0(i) )
2595	xtv1000 = dble( - ee * nu12_1000 ) /
2596	@ ( log( xalfa/(1.d0+xalfa+xbeta) ) -
2597	@ dble(ee*nu(1,2))/tvt0(i) )
2598	!!! correccion 8-Nov-04 (see pag.9b-Marte4-)
2599	xtv0200 = dble( - ee * nu12_0200 /
2600	@ (log(4.xbeta/(1.d0+xalfa+xbeta))-eenu(1,2)/tvt0(i)))
2601	xtv0220 = dble( - ee * nu(1,2) /
2602	@ ( log(2./(1.d0+xalfa+xbeta)) - ee*nu(1,2)/tvt0(i) ) )
2603	xtv1000 = dble( - ee * nu12_1000 /
2604	@ ( log(4.xalfa/(1.d0+xalfa+xbeta))-eenu(1,2)/tvt0(i)))
2605
2606	! if ( icurt_pop.eq.30 ) then
2607	! write (30,'( 1x,f7.2, 3x,f8.3, 2x,3(1x,f8.3) )')
2608	! @ zx(i), tvt0(i), xtv0200, xtv0220, xtv1000
2609	! endif
2610
2611	!!! xlower and xes for the band
2612	if (ib.eq.19) then
2613	xlower = exp( dble(eeelow(isot,ib))
2614	@ ( 1.d0/dble(ty(i))-1.d0/xtv0200 ) )
2615	xes = 1.0d0
2616	elseif (ib.eq.17) then
2617	xlower = exp( dble(eeelow(isot,ib))
2618	@ ( 1.d0/dble(ty(i))-1.d0/xtv1000 ) )
2619	xes = 1.0d0
2620	elseif (ib.eq.20) then
2621	xlower = exp( dble(eeelow(isot,ib))
2622	@ ( 1.d0/dble(ty(i))-1.d0/xtv0220 ) )
2623	xes = 1.0d0
2624	elseif (ib.eq.14) then
2625	xlower = exp( dble(eenu12_1000)
2626	@ ( 1.d0/dble(ty(i))-1.d0/xtv1000 ) )
2627	xnu_trans = dble( nu(1,4)-nu12_1000 )
2628	xtv_trans = xnu_trans / dble(nu(1,4)/tvt1(i)-
2629	@ nu12_1000/xtv1000)
2630	xes = (1.d0-exp( dble(-ee*xnu_trans/xtv_trans) )) /
2631	@ (1.d0-exp( dble(-ee*xnu_trans/ty(i)) ))
2632	elseif (ib.eq.13) then
2633	xlower = exp( dble(eenu12_0200)
2634	@ ( 1.d0/dble(ty(i))-1.d0/xtv0200 ) )
2635	xnu_trans = dble(nu(1,4)-nu12_0200)
2636	xtv_trans = xnu_trans / dble(nu(1,4)/tvt1(i)-
2637	@ nu12_0200/xtv0200)
2638	xes = (1.d0-exp( dble(-ee*xnu_trans/xtv_trans) )) /
2639	@ (1.d0-exp( dble(-ee*xnu_trans/ty(i)) ))
2640	else
2641	xlower = exp( dble(eeelow(isot,ib))
2642	@ ( 1.d0/dble(ty(i))-1.d0/tvt0(i) ) )
2643	xes = 1.0d0
2644	end if
2645	xqv = (1.d0-exp( dble(-ee*667.3801/tvtbs(i)) )) /
2646	@ (1.d0-exp( dble(-ee*667.3801/ty(i)) ))
2647	xfactor = xlower * xqv*2.d0 xes
2648
2649	elseif (isot.eq.2) then
2650
2651	xalfa = 1.d0/2.d0* exp( dble(-ee*(nu22_1000-nu(2,2))/
2652	@ ty(i)) )
2653	xbeta = 1.d0/2.d0* exp( dble(-ee*(nu22_0200-nu(2,2))/
2654	@ ty(i)) )
2655	xtv0200 = dble( - ee * nu22_0200 ) /
2656	@ ( log( xbeta/(1.d0+xalfa+xbeta) ) - dble(ee*nu(2,2))/
2657	@ tvt0(i) )
2658	xtv1000 = dble( - ee * nu22_1000 ) /
2659	@ ( log( xalfa/(1.d0+xalfa+xbeta) ) - dble(ee*nu(2,2))/
2660	@ tvt0(i) )
2661
2662	if (ib.eq.14) then
2663	xlower = exp( dble(eenu22_1000)
2664	@ ( 1.d0/dble(ty(i))-1.d0/xtv1000 ) )
2665	xnu_trans = dble(nu(2,4)-nu22_1000)
2666	xtv_trans = xnu_trans / dble(nu(2,4)/tvt1(i)-nu22_1000/
2667	@ xtv1000)
2668	xes = (1.d0-exp( dble(-ee*xnu_trans/xtv_trans) )) /
2669	@ (1.d0-exp( dble(-ee*xnu_trans/ty(i)) ))
2670	elseif (ib.eq.13) then
2671	xlower = exp( dble(eenu22_0200)
2672	@ ( 1.d0/dble(ty(i))-1.d0/xtv0200 ) )
2673	xnu_trans = dble( nu(2,4)-nu22_0200 )
2674	xtv_trans = xnu_trans / dble(nu(2,4)/tvt1(i)-nu22_0200/
2675	@ xtv0200)
2676	xes = (1.d0-exp( dble(-ee*xnu_trans/xtv_trans) )) /
2677	@ (1.d0-exp( dble(-ee*xnu_trans/ty(i)) ))
2678	else
2679	xlower = exp( dble(eeelow(isot,ib))
2680	@ ( 1.d0/dble(ty(i))-1.d0/tvt0(i) ) )
2681	xes = 1.0d0
2682	end if
2683	xqv = (1.d0-exp( dble(-ee*662.3734/tvtbs(i)) )) /
2684	@ (1.d0-exp( dble(-ee*662.3734/ty(i)) ))
2685	xfactor = xlower * xqv*2.d0 xes
2686
2687	elseif (isot.eq.3) then
2688
2689	xalfa = 1.d0/2.d0* exp( dble(-ee*(nu32_1000-nu(3,2))/
2690	@ ty(i)) )
2691	xbeta = 1.d0/2.d0* exp( dble(-ee*(nu32_0200-nu(3,2))/
2692	@ ty(i)) )
2693	xtv0200 = dble( - ee * nu32_0200 ) /
2694	@ ( log( xbeta/(1.d0+xalfa+xbeta) ) - dble(ee*nu(3,2))/
2695	@ tvt0(i) )
2696	xtv1000 = dble( - ee * nu32_1000 ) /
2697	@ ( log( xalfa/(1.d0+xalfa+xbeta) ) - dble(ee*nu(3,2))/
2698	@ tvt0(i) )
2699
2700	if (ib.eq.14) then
2701	xlower = exp( dble(eenu32_1000)
2702	@ ( 1.d0/dble(ty(i))-1.d0/xtv1000 ) )
2703	xnu_trans = dble( nu(3,4)-nu32_1000 )
2704	xtv_trans = xnu_trans / dble(nu(3,4)/tvt1(i)-nu32_1000/
2705	@ xtv1000)
2706	xes = (1.d0-exp( dble(-ee*xnu_trans/xtv_trans) )) /
2707	@ (1.d0-exp( dble(-ee*xnu_trans/ty(i)) ))
2708	elseif (ib.eq.13) then
2709	xlower = exp( dble(eenu32_0200)
2710	@ ( 1.d0/dble(ty(i))-1.d0/xtv0200 ) )
2711	xnu_trans = dble( nu(3,4)-nu32_0200 )
2712	xtv_trans = xnu_trans / dble(nu(3,4)/tvt1(i)-nu32_0200/
2713	@ xtv0200)
2714	xes = (1.d0-exp( dble(-ee*xnu_trans/xtv_trans) )) /
2715	@ (1.d0-exp( dble(-ee*xnu_trans/ty(i)) ))
2716	else
2717	xlower = exp( dble(eeelow(isot,ib))
2718	@ ( 1.d0/dble(ty(i))-1.d0/tvt0(i) ) )
2719	xes = 1.0d0
2720	end if
2721	xqv = (1.d0-exp( dble(-ee*648.4784/tvtbs(i)) )) /
2722	@ (1.d0-exp( dble(-ee*648.4784/ty(i)) ))
2723	xfactor = xlower * xqv*2.d0 xes
2724
2725	elseif (isot.eq.4) then
2726
2727	xalfa = 1.d0/2.d0* exp( dble(-ee*(nu42_1000-nu(4,2))/
2728	@ ty(i)) )
2729	xbeta = 1.d0/2.d0* exp( dble(-ee*(nu42_0200-nu(4,2))/
2730	@ ty(i)) )
2731	xtv0200 = dble( - ee * nu42_0200 ) /
2732	@ ( log( xbeta/(1.d0+xalfa+xbeta) ) - dble(ee*nu(4,2))/
2733	@ tvt0(i) )
2734	xtv1000 = dble( - ee * nu42_1000 ) /
2735	@ ( log( xalfa/(1.d0+xalfa+xbeta) ) - dble(ee*nu(4,2))/
2736	@ tvt0(i) )
2737
2738	if (ib.eq.14) then
2739	xlower = exp( dble(eenu42_1000)
2740	@ ( 1.d0/dble(ty(i))-1.d0/xtv1000 ) )
2741	xnu_trans = dble( nu(4,4)-nu42_1000 )
2742	xtv_trans = xnu_trans / dble(nu(4,4)/tvt1(i)-nu42_1000/
2743	@ xtv1000)
2744	xes = (1.d0-exp( dble(-ee*xnu_trans/xtv_trans) )) /
2745	@ (1.d0-exp( dble(-ee*xnu_trans/ty(i)) ))
2746	elseif (ib.eq.13) then
2747	xlower = exp( dble(eenu42_0200)
2748	$ ( 1.d0/dble(ty(i))-1.d0/xtv0200 ) )
2749	xnu_trans = dble( nu(4,4)-nu42_0200 )
2750	xtv_trans = xnu_trans / dble(nu(4,4)/tvt1(i)-nu42_0200/
2751	@ xtv0200)
2752	xes = (1.d0-exp( dble(-ee*xnu_trans/xtv_trans) )) /
2753	@ (1.d0-exp( dble(-ee*xnu_trans/ty(i)) ))
2754	else
2755	xlower = exp( dble(eeelow(isot,ib))
2756	@ ( 1.d0/dble(ty(i))-1.d0/tvt0(i) ) )
2757	xes = 1.0d0
2758	end if
2759	xqv = (1.d0-exp( dble(-ee*664.7289/tvtbs(i)) )) /
2760	@ (1.d0-exp( dble(-ee*664.7289/ty(i)) ))
2761	xfactor = xlower * xqv*2.d0 xes
2762
2763	elseif (isot.eq.5 .and. ib.eq.1) then
2764
2765	xlower = 1.d0
2766	xes = 1.0d0
2767	xqv = (1.d0-exp( dble(-ee*nuco_10/tvtbs(i)) )) /
2768	@ (1.d0-exp( dble(-ee*nuco_10/ty(i)) ))
2769	xfactor = xlower * xqv * xes
2770
2771	end if
2772
2773	con(i) = con(i) * xfactor
2774	if (i.eq.nzy) coninf = coninf * xfactor
2775
2776	end do
2777
2778	! if ( isot.eq.1 .and. icurt_pop.eq.30 ) then
2779	! close (30)
2780	! endif
2781
2782	return
2783	end
2784
2785
2786
2787
2788
2789	c***********************************************************************
2790	c mztf.f
2791	c***********************************************************************
2792	c
2793	c program for calculating atmospheric transmittances
2794	c to be used in the calculation of curtis matrix coefficients
2795
2796	subroutine mztf ( ig,cf,cfup,cfdw,vc,taugr, ib,isot,
2797	@ iirw,iimu,itauout,icfout,itableout )
2798
2799	c i*out = 1 output of data
2800	c i*out = 0 no output
2801	c
2802	c jul 2011 malv+fgg adapted to LMD-MGCM
2803	c nov 98 mavl allow for overlaping in the lorentz line
2804	c jan 98 malv version for mz1d. based on curtis/mztf.for
2805	c 17-jul-96 mlp&crs change the calculation of mr.
2806	c evitar: divide por cero. anhadiendo: ff
2807	c oct-92 malv correct s(t) dependence for all histogr bands
2808	c june-92 malv proper lower levels for laser bands
2809	c may-92 malv new temperature dependence for laser bands
2810	c @ 991 malv boxing for the averaged absorber amount and t
2811	c ? malv extension up to 200 km altitude in mars
2812	c 13-nov-86 mlp include the temperature weighted to match
2813	c the eqw in the strong doppler limit.
2814	c***********************************************************************
2815
2816	implicit none
2817
2818	include 'nlte_paramdef.h'
2819	include 'nlte_commons.h'
2820
2821
2822	c arguments
2823	integer ig !ADDED FOR TRACEBACK
2824	real*8 cf(nl,nl), cfup(nl,nl), cfdw(nl,nl) ! o.
2825	real*8 vc(nl), taugr(nl) ! o
2826	integer ib ! i
2827	integer isot ! i
2828	integer iirw ! i
2829	integer iimu ! i
2830	integer itauout ! i
2831	integer icfout ! i
2832	integer itableout ! i
2833
2834	c local variables and constants
2835	integer i, in, ir, im, k ,j
2836	integer nmu
2837	parameter (nmu = 8)
2838	real*8 tau(nl,nl)
2839	real*8 tauinf(nl)
2840	real*8 con(nzy), coninf
2841	real*8 c1, c2
2842	real*8 t1, t2
2843	real*8 p1, p2
2844	real*8 mr1, mr2
2845	real*8 st1, st2
2846	real*8 c1box(70), c2box(70)
2847	real*8 ff ! to avoid too small numbers
2848	real*8 tvtbs(nzy)
2849	real*8 st, beta, ts, eqwmu
2850	real*8 mu(nmu), amu(nmu)
2851	real*8 zld(nl), zyd(nzy)
2852	real*8 correc
2853	real deltanux ! width of vib-rot band (cm-1)
2854	! character isotcode*2
2855	integer idummy
2856	real*8 Desp,wsL
2857
2858	c formats
2859	! 111 format(a1)
2860	! 112 format(a2)
2861	101 format(i1)
2862	202 format(i2)
2863	! 180 format(a80)
2864	! 181 format(a80)
2865	c***********************************************************************
2866
2867	c some needed values
2868	! rl=sqrt(log(2.d0))
2869	! pi2 = 3.14159265358989d0
2870	beta = 1.8d0
2871	idummy = 0
2872	Desp = 0.d0
2873	wsL = 0.d0
2874
2875	c esto es para que las subroutines de mztfsub calculen we
2876	c de la forma apropiada para mztf, no para fot
2877	icls=icls_mztf
2878
2879	c codigos para filenames
2880	! if (isot .eq. 1) isotcode = '26'
2881	! if (isot .eq. 2) isotcode = '28'
2882	! if (isot .eq. 3) isotcode = '36'
2883	! if (isot .eq. 4) isotcode = '27'
2884	! if (isot .eq. 5) isotcode = '62'
2885	! if(ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5
2886	! @ .or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9) then
2887	! write (ibcode1,101) ib
2888	! else
2889	! write (ibcode2,202) ib
2890	! endif
2891	! write (*,'( 30h calculating curtis matrix : ,2x,
2892	! @ 8h band = ,i2,2x, 11h isotope = ,i2)') ib, isot
2893
2894	c integration in angle !!!!!!!!!!!!!!!!!!!!
2895
2896	c------- diffusivity approx.
2897	if (iimu.eq.1) then
2898	! write (,) ' diffusivity approx. beta = ',beta
2899	mu(1) = 1.0d0
2900	amu(1)= 1.0d0
2901	c-------data for 8 points integration
2902	elseif (iimu.eq.4) then
2903	write (,)' 4 points for the gauss-legendre angle quadrature.'
2904	mu(1)=(1.0d0+0.339981043584856)/2.0d0
2905	mu(2)=(1.0d0-0.339981043584856)/2.0d0
2906	mu(3)=(1.0d0+0.861136311594053)/2.0d0
2907	mu(4)=(1.0d0-0.861136311594053)/2.0d0
2908	amu(1)=0.652145154862546
2909	amu(2)=amu(1)
2910	amu(3)=0.347854845137454
2911	amu(4)=amu(3)
2912	beta=1.0d0
2913	c-------data for 8 points integration
2914	elseif(iimu.eq.8) then
2915	write (,)' 8 points for the gauss-legendre angle quadrature.'
2916	mu(1)=(1.0d0+0.183434642495650)/2.0d0
2917	mu(2)=(1.0d0-0.183434642495650)/2.0d0
2918	mu(3)=(1.0d0+0.525532409916329)/2.0d0
2919	mu(4)=(1.0d0-0.525532409916329)/2.0d0
2920	mu(5)=(1.0d0+0.796666477413627)/2.0d0
2921	mu(6)=(1.0d0-0.796666477413627)/2.0d0
2922	mu(7)=(1.0d0+0.960289856497536)/2.0d0
2923	mu(8)=(1.0d0-0.960289856497536)/2.0d0
2924	amu(1)=0.362683783378362
2925	amu(2)=amu(1)
2926	amu(3)=0.313706645877887
2927	amu(4)=amu(3)
2928	amu(5)=0.222381034453374
2929	amu(6)=amu(5)
2930	amu(7)=0.101228536290376
2931	amu(8)=amu(7)
2932	beta=1.0d0
2933	end if
2934	c!!!!!!!!!!!!!!!!!!!!!!!
2935
2936	ccc
2937	ccc determine abundances included in the absorber amount
2938	ccc
2939
2940	c first, set up the grid ready for interpolation.
2941	do i=1,nzy
2942	zyd(i) = dble(zy(i))
2943	enddo
2944	do i=1,nl
2945	zld(i) = dble(zl(i))
2946	enddo
2947
2948
2949	c 2nd: correccion a la n10(i) (cantidad de absorbente en el lower state)
2950	c por similitud a la que se hace en cza.for
2951
2952	do i=1,nzy
2953	if (isot.eq.5) then
2954	con(i) = dble( coy(i) * imrco )
2955	else
2956	con(i) = dble( co2y(i) * imr(isot) )
2957	c vibr. temp of the bending mode :
2958	if(isot.eq.1) call interdp(tvtbs,zyd,nzy,v626t1,zld,nl,1)
2959	if(isot.eq.2) call interdp(tvtbs,zyd,nzy,v628t1,zld,nl,1)
2960	if(isot.eq.3) call interdp(tvtbs,zyd,nzy,v636t1,zld,nl,1)
2961	if(isot.eq.4) call interdp(tvtbs,zyd,nzy,v627t1,zld,nl,1)
2962	correc = 2.d0 * dexp( dble(-ee*elow(isot,2))/tvtbs(i) )
2963	con(i) = con(i) * ( 1.d0 - correc )
2964	endif
2965	c-----------------------------------------------------------------------
2966	c mlp & cristina. 17 july 1996
2967	c change the calculation of mr. it is used for calculating partial press
2968	c alpha = alpha(self,co2)pp +alpha(n2)(pt-pp)
2969	c for an isotope, if mr is obtained by co2*imr(iso)/nt we are considerin
2970	c collisions with other co2 isotopes (including the major one, 626)
2971	c as if they were with n2. assuming mr as co2/nt, we consider collisions
2972	c of type 628-626 as of 626-626 instead of as 626-n2.
2973	c mrx(i)=con(i)/ntx(i) ! old malv
2974
2975	! mrx(i)= dble(co2x(i)/ntx(i)) ! mlp & crs
2976
2977	c jan 98:
2978	c esta modif de mlp implica anular el correc (deberia revisar esto)
2979	mr(i) = dble(co2y(i)/nty(i)) ! malv, jan 98
2980
2981	c-----------------------------------------------------------------------
2982
2983	end do
2984
2985	! como beta y 1.d5 son comunes a todas las weighted absorber amounts,
2986	! los simplificamos:
2987	! coninf = beta * 1.d5 * dble( con(n) / log( con(n-1) / con(n) ) )
2988	coninf = dble( con(nzy) / log( con(nzy-1) / con(nzy) ) )
2989
2990	! write (,) ' coninf =', coninf
2991
2992	ccc
2993	ccc temp dependence of the band strength and
2994	ccc nlte correction factor for the absorber amount
2995	ccc
2996	call mztf_correccion ( coninf, con, ib, isot, itableout )
2997
2998	ccc
2999	ccc reads histogrammed spectral data (strength for lte and vmr=1)
3000	ccc
3001	!hfile1 = dirspec//'hi'//dn ! ya no distinguimos entre d/n
3002	!! hfile1 = dirspec//'hid' ! (see why in his.for)
3003	! hfile='hid'
3004	!! if (ib.eq.13 .or. ib.eq.14 ) hfile1 = dirspec//'his'
3005	! if (ib.eq.13 .or. ib.eq.14 ) hfile1 = 'his'
3006	!
3007	! if(ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5
3008	! @ .or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9) then
3009	! if (isot.eq.1) hisfile = hfile1//'26-'//ibcode1//'.dat'
3010	! if (isot.eq.2) hisfile = hfile1//'28-'//ibcode1//'.dat'
3011	! if (isot.eq.3) hisfile = hfile1//'36-'//ibcode1//'.dat'
3012	! if (isot.eq.4) hisfile = hfile1//'27-'//ibcode1//'.dat'
3013	! if (isot.eq.5) hisfile = hfile1//'62-'//ibcode1//'.dat'
3014	! else
3015	! if (isot.eq.1) hisfile = hfile1//'26-'//ibcode2//'.dat'
3016	! if (isot.eq.2) hisfile = hfile1//'28-'//ibcode2//'.dat'
3017	! if (isot.eq.3) hisfile = hfile1//'36-'//ibcode2//'.dat'
3018	! if (isot.eq.4) hisfile = hfile1//'27-'//ibcode2//'.dat'
3019	! if (isot.eq.5) hisfile = hfile1//'62-'//ibcode2//'.dat'
3020	! endif
3021	! write (,) 'hisfile: ', hisfile
3022
3023	! the argument to rhist is to make this compatible with mztf_comp.f,
3024	! which is a useful modification of mztf.f (to change strengths of bands
3025	! call rhist (1.0)
3026	if(ib.eq.1) then
3027	if(isot.eq.1) then !Case 1
3028	mm=mm_c1
3029	nbox=nbox_c1
3030	tmin=tmin_c1
3031	tmax=tmax_c1
3032	do i=1,nbox_max
3033	no(i)=no_c1(i)
3034	dist(i)=dist_c1(i)
3035	do j=1,nhist
3036	sk1(j,i)=sk1_c1(j,i)
3037	xls1(j,i)=xls1_c1(j,i)
3038	xln1(j,i)=xln1_c1(j,i)
3039	xld1(j,i)=xld1_c1(j,i)
3040	enddo
3041	enddo
3042	do j=1,nhist
3043	thist(j)=thist_c1(j)
3044	enddo
3045	else if(isot.eq.2) then !Case 2
3046	mm=mm_c2
3047	nbox=nbox_c2
3048	tmin=tmin_c2
3049	tmax=tmax_c2
3050	do i=1,nbox_max
3051	no(i)=no_c2(i)
3052	dist(i)=dist_c2(i)
3053	do j=1,nhist
3054	sk1(j,i)=sk1_c2(j,i)
3055	xls1(j,i)=xls1_c2(j,i)
3056	xln1(j,i)=xln1_c2(j,i)
3057	xld1(j,i)=xld1_c2(j,i)
3058	enddo
3059	enddo
3060	do j=1,nhist
3061	thist(j)=thist_c2(j)
3062	enddo
3063	else if(isot.eq.3) then !Case 3
3064	mm=mm_c3
3065	nbox=nbox_c3
3066	tmin=tmin_c3
3067	tmax=tmax_c3
3068	do i=1,nbox_max
3069	no(i)=no_c3(i)
3070	dist(i)=dist_c3(i)
3071	do j=1,nhist
3072	sk1(j,i)=sk1_c3(j,i)
3073	xls1(j,i)=xls1_c3(j,i)
3074	xln1(j,i)=xln1_c3(j,i)
3075	xld1(j,i)=xld1_c3(j,i)
3076	enddo
3077	enddo
3078	do j=1,nhist
3079	thist(j)=thist_c3(j)
3080	enddo
3081	else if(isot.eq.4) then !Case 4
3082	mm=mm_c4
3083	nbox=nbox_c4
3084	tmin=tmin_c4
3085	tmax=tmax_c4
3086	do i=1,nbox_max
3087	no(i)=no_c4(i)
3088	dist(i)=dist_c4(i)
3089	do j=1,nhist
3090	sk1(j,i)=sk1_c4(j,i)
3091	xls1(j,i)=xls1_c4(j,i)
3092	xln1(j,i)=xln1_c4(j,i)
3093	xld1(j,i)=xld1_c4(j,i)
3094	enddo
3095	enddo
3096	do j=1,nhist
3097	thist(j)=thist_c4(j)
3098	enddo
3099	else
3100	write(,)'isot must be 2,3 or 4 for ib=1!!'
3101	write(,)'stop at mztf_overlap/317'
3102	stop
3103	endif
3104	else if (ib.eq.2) then
3105	if(isot.eq.1) then !Case 5
3106	mm=mm_c5
3107	nbox=nbox_c5
3108	tmin=tmin_c5
3109	tmax=tmax_c5
3110	do i=1,nbox_max
3111	no(i)=no_c5(i)
3112	dist(i)=dist_c5(i)
3113	do j=1,nhist
3114	sk1(j,i)=sk1_c5(j,i)
3115	xls1(j,i)=xls1_c5(j,i)
3116	xln1(j,i)=xln1_c5(j,i)
3117	xld1(j,i)=xld1_c5(j,i)
3118	enddo
3119	enddo
3120	do j=1,nhist
3121	thist(j)=thist_c5(j)
3122	enddo
3123	else
3124	write(,)'isot must be 1 for ib=2!!'
3125	write(,)'stop at mztf_overlap/341'
3126	stop
3127	endif
3128	else if (ib.eq.3) then
3129	if(isot.eq.1) then !Case 6
3130	mm=mm_c6
3131	nbox=nbox_c6
3132	tmin=tmin_c6
3133	tmax=tmax_c6
3134	do i=1,nbox_max
3135	no(i)=no_c6(i)
3136	dist(i)=dist_c6(i)
3137	do j=1,nhist
3138	sk1(j,i)=sk1_c6(j,i)
3139	xls1(j,i)=xls1_c6(j,i)
3140	xln1(j,i)=xln1_c6(j,i)
3141	xld1(j,i)=xld1_c6(j,i)
3142	enddo
3143	enddo
3144	do j=1,nhist
3145	thist(j)=thist_c6(j)
3146	enddo
3147	else
3148	write(,)'isot must be 1 for ib=3!!'
3149	write(,)'stop at mztf_overlap/365'
3150	stop
3151	endif
3152	else if (ib.eq.4) then
3153	if(isot.eq.1) then !Case 7
3154	mm=mm_c7
3155	nbox=nbox_c7
3156	tmin=tmin_c7
3157	tmax=tmax_c7
3158	do i=1,nbox_max
3159	no(i)=no_c7(i)
3160	dist(i)=dist_c7(i)
3161	do j=1,nhist
3162	sk1(j,i)=sk1_c7(j,i)
3163	xls1(j,i)=xls1_c7(j,i)
3164	xln1(j,i)=xln1_c7(j,i)
3165	xld1(j,i)=xld1_c7(j,i)
3166	enddo
3167	enddo
3168	do j=1,nhist
3169	thist(j)=thist_c7(j)
3170	enddo
3171	else
3172	write(,)'isot must be 1 for ib=4!!'
3173	write(,)'stop at mztf_overlap/389'
3174	stop
3175	endif
3176	else
3177	write(,)'ib must be 1,2,3 or 4!!'
3178	write(,)'stop at mztf_overlap/394'
3179	endif
3180
3181	if (isot.ne.5) deltanux = deltanu(isot,ib)
3182	if (isot.eq.5) deltanux = deltanuco
3183
3184	c******
3185	c****** calculation of tauinf(nl)
3186	c******
3187	call initial
3188
3189	ff=1.0e10
3190
3191	do i=nl,1,-1
3192
3193	if(i.eq.nl)then
3194
3195	call intz (zl(i),c2,p2,mr2,t2, con)
3196	do kr=1,nbox
3197	ta(kr)=t2
3198	end do
3199	! write (,) ' i, t2 =', i, t2
3200	call interstrength (st2,t2,ka,ta)
3201	aa = p2 * coninf * mr2 * (st2 * ff)
3202	bb = p2 * coninf * st2
3203	cc = coninf * st2
3204	dd = t2 * coninf * st2
3205	do kr=1,nbox
3206	ccbox(kr) = coninf * ka(kr)
3207	ddbox(kr) = t2 * ccbox(kr)
3208	! c2box(kr) = c2 * ka(kr) * beta * dble(deltaz) * 1.d5
3209	c2box(kr) = c2 * ka(kr) * dble(deltaz)
3210	end do
3211	! c2 = c2 * st2 * beta * dble(deltaz) * 1.d5
3212	c2 = c2 * st2 * dble(deltaz)
3213
3214	else
3215	call intz (zl(i),c1,p1,mr1,t1, con)
3216	do kr=1,nbox
3217	ta(kr)=t1
3218	end do
3219	! write (,) ' i, t1 =', i, t1
3220	call interstrength (st1,t1,ka,ta)
3221	do kr=1,nbox
3222	! c1box(kr) = c1 * ka(kr) * beta * dble(deltaz) * 1.d5
3223	c1box(kr) = c1 * ka(kr) * dble(deltaz)
3224	end do
3225	! c1 = c1 * st1 * beta * dble(deltaz) * 1.d5
3226	c1 = c1 * st1 * dble(deltaz)
3227	aa = aa + ( p1mr1(c1ff) + p2mr2(c2ff)) / 2.d0
3228	bb = bb + ( p1c1 + p2c2 ) / 2.d0
3229	cc = cc + ( c1 + c2 ) / 2.d0
3230	dd = dd + ( t1c1 + t2c2 ) / 2.d0
3231	do kr=1,nbox
3232	ccbox(kr) = ccbox(kr) + ( c1box(kr) + c2box(kr) )/2.d0
3233	ddbox(kr) = ddbox(kr) + (t1c1box(kr)+t2c2box(kr))/2.d0
3234	end do
3235
3236	mr2 = mr1
3237	c2=c1
3238	do kr=1,nbox
3239	c2box(kr) = c1box(kr)
3240	end do
3241	t2=t1
3242	p2=p1
3243	end if
3244
3245	pt = bb / cc
3246	pp = aa / (cc*ff)
3247
3248	! ta=dd/cc
3249	! tdop = ta
3250	ts = dd/cc
3251	do kr=1,nbox
3252	ta(kr) = ddbox(kr) / ccbox(kr)
3253	end do
3254	! write (,) ' i, ts =', i, ts
3255	call interstrength(st,ts,ka,ta)
3256	! call intershape(alsa,alna,alda,tdop)
3257	call intershape(alsa,alna,alda,ta)
3258
3259	* ua = cc/st
3260
3261	c next loop calculates the eqw for an especified path ua,pp,pt,ta
3262
3263	eqwmu = 0.0d0
3264	do im = 1,iimu
3265	eqw=0.0d0
3266	do kr=1,nbox
3267	ua(kr) = ccbox(kr) / ka(kr) * beta * 1.0d5 / mu(im)
3268	if(ua(kr).lt.0.)write(,)'mztf_overlap/483',ua(kr),
3269	$ ccbox(kr),ka(kr),beta,mu(im),kr,im,i,nl
3270
3271	call findw(ig,iirw, idummy,c1,p1,Desp,wsL)
3272	if ( i_supersat .eq. 0 ) then
3273	eqw=eqw+no(kr)*w
3274	else
3275	eqw = w + (no(kr)-1) * ( asat_box*dist(kr) )
3276	endif
3277	end do
3278	eqwmu = eqwmu + eqw * mu(im)*amu(im)
3279	end do
3280
3281	tauinf(i) = exp( - eqwmu / dble(deltanux) )
3282
3283	end do ! i continue
3284
3285	! if ( isot.eq.1 .and. ib.eq.2 ) then
3286	! write (,) ' tauinf(nl) = ', tauinf(nl)
3287	! write (,) ' tauinf(1) = ', tauinf(1)
3288	! endif
3289
3290	c******
3291	c****** calculation of tau(in,ir) for n<=r
3292	c******
3293
3294	do 1 in=1,nl-1
3295
3296	call initial
3297	call intz (zl(in), c1,p1,mr1,t1, con)
3298	do kr=1,nbox
3299	ta(kr) = t1
3300	end do
3301	call interstrength (st1,t1,ka,ta)
3302	do kr=1,nbox
3303	! c1box(kr) = c1 * ka(kr) * beta * dble(deltaz) * 1.d5
3304	c1box(kr) = c1 * ka(kr) * dble(deltaz)
3305	end do
3306	! c1 = c1 * st1 * beta * dble(deltaz) * 1.d5
3307	c1 = c1 * st1 * dble(deltaz)
3308
3309	do 2 ir=in,nl-1
3310
3311	if (ir.eq.in) then
3312	tau(in,ir) = 1.d0
3313	goto 2
3314	end if
3315
3316	call intz (zl(ir), c2,p2,mr2,t2, con)
3317	do kr=1,nbox
3318	ta(kr) = t2
3319	end do
3320	call interstrength (st2,t2,ka,ta)
3321	do kr=1,nbox
3322	! c2box(kr) = c2 * ka(kr) * beta * dble(deltaz) * 1.d5
3323	c2box(kr) = c2 * ka(kr) * dble(deltaz)
3324	end do
3325	! c2 = c2 * st2 * beta * dble(deltaz) * 1.e5
3326	c2 = c2 * st2 * dble(deltaz)
3327
3328	c aa = aa + ( p1mr1c1 + p2mr2c2 ) / 2.d0
3329	aa = aa + ( p1mr1(c1ff) + p2mr2(c2ff)) / 2.d0
3330	bb = bb + ( p1c1 + p2c2 ) / 2.d0
3331	cc = cc + ( c1 + c2 ) / 2.d0
3332	dd = dd + ( t1c1 + t2c2 ) / 2.d0
3333	do kr=1,nbox
3334	ccbox(kr) = ccbox(kr) + ( c1box(kr) + c2box(kr) ) /2.d0
3335	ddbox(kr) = ddbox(kr) +
3336	$ ( t1c1box(kr) + t2c2box(kr) ) /2.d0
3337	end do
3338
3339	mr1=mr2
3340	t1=t2
3341	c1=c2
3342	p1=p2
3343	do kr=1,nbox
3344	c1box(kr) = c2box(kr)
3345	end do
3346
3347	pt = bb / cc
3348	pp = aa / (cc * ff)
3349
3350	* ta=dd/cc
3351	* tdop = ta
3352	ts = dd/cc
3353	do kr=1,nbox
3354	ta(kr) = ddbox(kr) / ccbox(kr)
3355	end do
3356	call interstrength(st,ts,ka,ta)
3357	call intershape(alsa,alna,alda,ta)
3358	* ua = cc/st
3359
3360	c next loop calculates the eqw for an especified path ua,pp,pt,ta
3361
3362	eqwmu = 0.0d0
3363	do im = 1,iimu
3364	eqw=0.0d0
3365	do kr=1,nbox
3366	ua(kr) = ccbox(kr) / ka(kr) * beta * 1.0d5 / mu(im)
3367	if(ua(kr).lt.0.)write(,)'mztf_overlap/581',ua(kr),
3368	$ ccbox(kr),ka(kr),beta,mu(im),kr,im,i,nl
3369
3370	call findw(ig,iirw, idummy,c1,p1,Desp,wsL)
3371	if ( i_supersat .eq. 0 ) then
3372	eqw=eqw+no(kr)*w
3373	else
3374	eqw = w + (no(kr)-1) * ( asat_box*dist(kr) )
3375	endif
3376	end do
3377	eqwmu = eqwmu + eqw * mu(im)*amu(im)
3378	end do
3379
3380	tau(in,ir) = exp( - eqwmu / dble(deltanux) )
3381
3382	2 continue
3383
3384	1 continue
3385
3386	! if ( isot.eq.1 .and. ib.eq.2 ) then
3387	! write (,) ' tau(1,) , =1,20 '
3388	! write (,) ( sngl(tau(1,k)), k=1,20 )
3389	! endif
3390
3391
3392	c**********
3393	c********** calculation of tau(in,ir) for n>r
3394	c**********
3395
3396	in=nl
3397
3398	call initial
3399	call intz (zl(in), c1,p1,mr1,t1, con)
3400	do kr=1,nbox
3401	ta(kr) = t1
3402	end do
3403	call interstrength (st1,t1,ka,ta)
3404	do kr=1,nbox
3405	! c1box(kr) = c1 * ka(kr) * beta * dble(deltaz) * 1.d5
3406	c1box(kr) = c1 * ka(kr) * dble(deltaz)
3407	end do
3408	! c1 = c1 * st1 * beta * dble(deltaz) * 1.d5
3409	c1 = c1 * st1 * dble(deltaz)
3410
3411	do 4 ir=in-1,1,-1
3412
3413	call intz (zl(ir), c2,p2,mr2,t2, con)
3414	do kr=1,nbox
3415	ta(kr) = t2
3416	end do
3417	call interstrength (st2,t2,ka,ta)
3418	do kr=1,nbox
3419	! c2box(kr) = c2 * ka(kr) * beta * dble(deltaz) * 1.d5
3420	c2box(kr) = c2 * ka(kr) * dble(deltaz)
3421	end do
3422	! c2 = c2 * st2 * beta * dble(deltaz) * 1.d5
3423	c2 = c2 * st2 * dble(deltaz)
3424
3425	aa = aa + ( p1mr1(c1ff) + p2mr2(c2ff)) / 2.d0
3426	bb = bb + ( p1c1 + p2c2 ) / 2.d0
3427	cc = cc + ( c1 + c2 ) / 2.d0
3428	dd = dd + ( t1c1 + t2c2 ) / 2.d0
3429	do kr=1,nbox
3430	ccbox(kr) = ccbox(kr) + ( c1box(kr) + c2box(kr) ) /2.d0
3431	ddbox(kr) = ddbox(kr) + ( t1c1box(kr) + t2c2box(kr) )/2.d0
3432	end do
3433
3434	mr1=mr2
3435	c1=c2
3436	t1=t2
3437	p1=p2
3438	do kr=1,nbox
3439	c1box(kr) = c2box(kr)
3440	end do
3441
3442	pt = bb / cc
3443	pp = aa / (cc * ff)
3444	ts = dd / cc
3445	do kr=1,nbox
3446	ta(kr) = ddbox(kr) / ccbox(kr)
3447	end do
3448	call interstrength (st,ts,ka,ta)
3449	call intershape (alsa,alna,alda,ta)
3450
3451	* ua = cc/st
3452
3453	c next loop calculates the eqw for an especified path ua,pp,pt,ta
3454
3455	eqwmu = 0.0d0
3456	do im = 1,iimu
3457	eqw=0.0d0
3458	do kr=1,nbox
3459	ua(kr) = ccbox(kr) / ka(kr) * beta * 1.0d5 / mu(im)
3460	if(ua(kr).lt.0.)write(,)'mztf_overlap/674',ua(kr),
3461	$ ccbox(kr),ka(kr),beta,mu(im),kr,im,i,nl
3462
3463	call findw(ig,iirw, idummy,c1,p1,Desp,wsL)
3464	if ( i_supersat .eq. 0 ) then
3465	eqw=eqw+no(kr)*w
3466	else
3467	eqw = w + (no(kr)-1) * ( asat_box*dist(kr) )
3468	endif
3469	end do
3470	eqwmu = eqwmu + eqw * mu(im)*amu(im)
3471	end do
3472
3473	tau(in,ir) = exp( - eqwmu / dble(deltanux) )
3474
3475	4 continue
3476
3477	c
3478	c due to the simmetry of the transmittances
3479	c
3480	do in=nl-1,2,-1
3481	do ir=in-1,1,-1
3482	tau(in,ir) = tau(ir,in)
3483	end do
3484	end do
3485
3486
3487	ccc
3488	ccc writing out transmittances
3489	ccc
3490	if (itauout.eq.1) then
3491
3492	! if(ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5
3493	! @ .or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9) then
3494	! open( 1, file=
3495	! @ dircurtis//'taul'//isotcode//dn//ibcode1//'.dat',
3496	! @ access='sequential', form='unformatted' )
3497	! else
3498	! open( 1, file=
3499	! @ dircurtis//'taul'//isotcode//dn//ibcode2//'.dat',
3500	! @ access='sequential', form='unformatted' )
3501	! endif
3502
3503	! write(1) dummy
3504	! write(1)' format: (tauinf(n),(tau(n,r),r=1,nl),n=1,nl)'
3505	! do in=1,nl
3506	! write (1) tauinf(in), ( tau(in,ir), ir=1,nl )
3507	! end do
3508	! close(unit=1)
3509
3510	elseif (itauout.eq.2) then
3511
3512	! if(ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5
3513	! @ .or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9) then
3514	! open( 1, file=
3515	! @ dircurtis//'taul'//isotcode//dn//ibcode1//'.dat')
3516	! else
3517	! open( 1, file=
3518	! @ dircurtis//'taul'//isotcode//dn//ibcode2//'.dat')
3519	! endif
3520
3521	! !write(1,*) dummy
3522	! !write(1,*) 'tij for curtis matrix calculations '
3523	! !write(1,*)' cira mars model atmosphere '
3524	! write(1,*)' beta= ',beta,'deltanu= ',deltanux
3525	! write(1,*)' number of elements (in,ir)= ',nl,nl
3526	! write(1,*)' format: (tauinf(in),(tau(in,ir),ir=1,nl),in=1,nl)'
3527
3528	! do in=1,nl
3529	! write (1,*) tauinf(in)
3530	! do ir=1,nl
3531	! write(1,*) tau(in,ir)
3532	! end do
3533	! end do
3534	! close(unit=1)
3535
3536	! if(ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5
3537	! @ .or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9) then
3538	! write (*,'(1x, 31htransmitances written out in: ,a22)')
3539	! @ 'taul'//isotcode//dn//ibcode1
3540	! else
3541	! write (*,'(1x, 31htransmitances written out in: ,a22)')
3542	! @ 'taul'//isotcode//dn//ibcode2
3543	! endif
3544
3545	end if
3546
3547	c cleaning of transmittances
3548	! call elimin_tau(tau,tauinf,nl,nan,itableout,nw,dummy,
3549	! @ isotcode,dn,ibcode2)
3550
3551	c construction of the curtis matrix
3552
3553	call mzcf ( tauinf,tau, cf,cfup,cfdw, vc,taugr,
3554	@ ib,isot,icfout,itableout )
3555
3556
3557	c end
3558	return
3559	end
3560
3561
3562
3563
3564	c***********************************************************************
3565	c mzcf
3566	c***********************************************************************
3567
3568	subroutine mzcf( tauinf,tau, c,cup,cdw,vc,taugr,
3569	@ ib,isot,icfout,itableout )
3570
3571	c a.k.murphy method to avoid extrapolation in the curtis matrix
3572	c feb-89 m. angel granada
3573	c 25-sept-96 cristina dejar las matrices en doble precision
3574	c jan 98 malv version para mz1d
3575	c jul 2011 malv+fgg adapted to LMD-MGCM
3576	c***********************************************************************
3577
3578	implicit none
3579
3580	include 'comcstfi.h'
3581	include 'nlte_paramdef.h'
3582	include 'nlte_commons.h'
3583
3584	c arguments
3585	real*8 c(nl,nl), cup(nl,nl), cdw(nl,nl) ! o
3586	real*8 vc(nl), taugr(nl) ! o
3587	real*8 tau(nl,nl) ! i
3588	real*8 tauinf(nl) ! i
3589	integer ib ! i
3590	integer isot ! i
3591	integer icfout, itableout ! i
3592
3593	c external
3594	external bandid
3595	character*2 bandid
3596
3597	c local variables
3598	integer i, in, ir, iw
3599	real*8 cfup(nl,nl), cfdw(nl,nl)
3600	real*8 a(nl,nl), cf(nl,nl)
3601	character isotcode2, bcode2
3602
3603	c formats
3604	101 format(i1)
3605	202 format(i2)
3606	180 format(a80)
3607	181 format(a80)
3608	c***********************************************************************
3609
3610	if (isot.eq.1) isotcode = '26'
3611	if (isot.eq.2) isotcode = '28'
3612	if (isot.eq.3) isotcode = '36'
3613	if (isot.eq.4) isotcode = '27'
3614	if (isot.eq.5) isotcode = 'co'
3615	bcode = bandid( ib )
3616
3617	! write (,) ' '
3618
3619	do in=1,nl
3620
3621	do ir=1,nl
3622
3623	cf(in,ir) = 0.0d0
3624	cfup(in,ir) = 0.0d0
3625	cfdw(in,ir) = 0.0d0
3626	c(in,ir) = 0.0d0
3627	cup(in,ir) = 0.0d0
3628	cdw(in,ir) = 0.0d0
3629	a(in,ir) = 0.0d0
3630
3631	end do
3632
3633	vc(in) = 0.0d0
3634	taugr(in) = 0.0d0
3635
3636	end do
3637
3638
3639	c the next lines are a reduced and equivalent way of calculating
3640	c the c(in,ir) elements for n=2,nl1 and r=1,nl
3641
3642
3643	c do in=2,nl1
3644	c do ir=1,nl
3645	c if(ir.eq.1)then
3646	c c(in,ir)=tau(in-1,1)-tau(in+1,1)
3647	c elseif(ir.eq.nl)then
3648	c c(in,ir)=tau(in+1,nl1)-tauinf(in+1)-tau(in-1,nl1)+tauinf(in-1)
3649	c else
3650	c c(in,ir)=tau(in+1,ir-1)-tau(in+1,ir)-tau(in-1,ir-1)+tau(in-1,ir)
3651	c end if
3652	c c(in,ir)=c(in,ir)pideltanu(ib)/(2.deltaz1.0e5)
3653	c end do
3654	c end do
3655	c go to 1000
3656
3657	c calculation of the matrix cfup(nl,nl)
3658
3659	cfup(1,1) = 1.d0 - tau(1,1)
3660
3661	do in=2,nl
3662	do ir=1,in
3663
3664	if (ir.eq.1) then
3665	cfup(in,ir) = tau(in,ir) - tau(in,1)
3666	elseif (ir.eq.in) then
3667	cfup(in,ir) = 1.d0 - tau(in,ir-1)
3668	else
3669	cfup(in,ir) = tau(in,ir) - tau(in,ir-1)
3670	end if
3671
3672	end do
3673	end do
3674
3675	! contribution to upwards fluxes from bb at bottom :
3676	do in=1,nl
3677	taugr(in) = tau(in,1)
3678	enddo
3679
3680	c calculation of the matrix cfdw(nl,nl)
3681
3682	cfdw(nl,nl) = 1.d0 - tauinf(nl)
3683
3684	do in=1,nl-1
3685	do ir=in,nl
3686
3687	if (ir.eq.in) then
3688	cfdw(in,ir) = 1.d0 - tau(in,ir)
3689	elseif (ir.eq.nl) then
3690	cfdw(in,ir) = tau(in,ir-1) - tauinf(in)
3691	else
3692	cfdw(in,ir) = tau(in,ir-1) - tau(in,ir)
3693	end if
3694
3695	end do
3696	end do
3697
3698
3699	c calculation of the matrix cf(nl,nl)
3700
3701	do in=1,nl
3702	do ir=1,nl
3703
3704	if (ir.eq.1) then
3705	! version con l_bb(tg) = l_bb(t(1))=j(1) (see also vc below)
3706	! cf(in,ir) = tau(in,ir)
3707	! version con l_bb(tg) =/= l_bb(t(1))=j(1) (see also vc below)
3708	cf(in,ir) = tau(in,ir) - tau(in,1)
3709	elseif (ir.eq.nl) then
3710	cf(in,ir) = tauinf(in) - tau(in,ir-1)
3711	else
3712	cf(in,ir) = tau(in,ir) - tau(in,ir-1)
3713	end if
3714
3715	end do
3716	end do
3717
3718
3719	c definition of the a(nl,nl) matrix
3720
3721	do in=2,nl-1
3722	do ir=1,nl
3723	if (ir.eq.in+1) a(in,ir) = -1.d0
3724	if (ir.eq.in-1) a(in,ir) = +1.d0
3725	a(in,ir) = a(in,ir) / ( 2.d0deltaz1.d5 )
3726	end do
3727	end do
3728	! this is not needed anymore in the akm scheme
3729	! a(1,1) = +3.d0
3730	! a(1,2) = -4.d0
3731	! a(1,3) = +1.d0
3732	! a(nl,nl) = -3.d0
3733	! a(nl,nl1) = +4.d0
3734	! a(nl,nl2) = -1.d0
3735
3736	c calculation of the final curtis matrix ("reduced" by murphy's method)
3737
3738	if (isot.ne.5) then
3739	do in=1,nl
3740	do ir=1,nl
3741	cf(in,ir) = cf(in,ir) * pi*deltanu(isot,ib)
3742	cfup(in,ir) = cfup(in,ir) * pi*deltanu(isot,ib)
3743	cfdw(in,ir) = cfdw(in,ir) * pi*deltanu(isot,ib)
3744	end do
3745	taugr(in) = taugr(in) * pi*deltanu(isot,ib)
3746	end do
3747	else
3748	do in=1,nl
3749	do ir=1,nl
3750	cf(in,ir) = cf(in,ir) * pi*deltanuco
3751	enddo
3752	taugr(in) = taugr(in) * pi*deltanuco
3753	enddo
3754	endif
3755
3756	do in=2,nl-1
3757
3758	do ir=1,nl
3759
3760	do i=1,nl
3761	! only c contains the matrix a. matrixes cup,cdw dont because
3762	! these two will be used for flux calculations, not
3763	! only for flux divergencies
3764
3765	c(in,ir) = c(in,ir) + a(in,i) * cf(i,ir)
3766	! from this matrix we will extract (see below) the
3767	! nl2 x nl2 "core" for the "reduced" final curtis matrix.
3768
3769	end do
3770	cup(in,ir) = cfup(in,ir)
3771	cdw(in,ir) = cfdw(in,ir)
3772
3773	end do
3774	! version con l_bb(tg) = l_bb(t(1))=j(1) (see cf above)
3775	!vc(in) = c(in,1)
3776	! version con l_bb(tg) =/= l_bb(t(1))=j(1) (see cf above)
3777	vc(in) = pideltanu(isot,ib)/( 2.d0deltaz1.d5 )
3778	@ ( tau(in-1,1) - tau(in+1,1) )
3779
3780	end do
3781
3782	5 continue
3783
3784	! write (,) 'mztf/1/ c(2,*) =', (c(2,i), i=1,nl)
3785
3786	! call elimin_dibuja(c,nl,itableout)
3787
3788	c ventana del smoothing de c es nw=3 y de vc es 5 (puesto en lisa):
3789	c subroutine elimin_mz4(c,vc,ilayer,nl,nan,iw, nw)
3790
3791	iw = nan
3792	if (isot.eq.4) iw = 5
3793	call elimin_mz1d (c,vc,0,iw,itableout,nw)
3794
3795	! upper boundary condition
3796	! j'(nl) = j'(nl1) ==> j(nl) = 2j(nl1) - j(nl2) ==>
3797	do in=2,nl-1
3798	c(in,nl-2) = c(in,nl-2) - c(in,nl)
3799	c(in,nl-1) = c(in,nl-1) + 2.d0*c(in,nl)
3800	cup(in,nl-2) = cup(in,nl-2) - cup(in,nl)
3801	cup(in,nl-1) = cup(in,nl-1) + 2.d0*cup(in,nl)
3802	cdw(in,nl-2) = cdw(in,nl-2) - cdw(in,nl)
3803	cdw(in,nl-1) = cdw(in,nl-1) + 2.d0*cdw(in,nl)
3804	end do
3805	! j(nl) = j(nl1) ==>
3806	! do in=2,nl1
3807	! c(in,nl1) = c(in,nl1) + c(in,nl)
3808	! end do
3809
3810	! 1000 continue
3811
3812	if (icfout.eq.1) then
3813
3814	! if (ib.eq.1 .or. ib.eq.12 .or. ib.eq.16 .or. ib.eq.18) then
3815	! codmatrx = codmatrx_fb
3816	! else
3817	! codmatrx = codmatrx_hot
3818	! end if
3819
3820	! if(ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5
3821	! @ .or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9) then
3822
3823	! open ( 1, access='sequential', form='unformatted', file=
3824	! @ dircurtis//'cfl'//isotcode//dn//ibcode1//codmatrx//'.dat')
3825	! open ( 2, access='sequential', form='unformatted', file=
3826	! @ dircurtis//'cflup'//isotcode//dn//ibcode1//codmatrx//'.dat')
3827	! open ( 3, access='sequential', form='unformatted', file=
3828	! @ dircurtis//'cfldw'//isotcode//dn//ibcode1//codmatrx//'.dat')
3829
3830	! else
3831
3832	! open ( 1, access='sequential', form='unformatted', file=
3833	! @ dircurtis//'cfl'//isotcode//dn//ibcode2//codmatrx//'.dat')
3834	! open ( 2, access='sequential', form='unformatted', file=
3835	! @ dircurtis//'cflup'//isotcode//dn//ibcode2//codmatrx//'.dat')
3836	! open ( 3, access='sequential', form='unformatted', file=
3837	! @ dircurtis//'cfldw'//isotcode//dn//ibcode2//codmatrx//'.dat')
3838
3839	! endif
3840
3841	! write(1) dummy
3842	! write(1)' format: (vc(n),(ch(n,r),r=2,nl-1),n=2,nl-1)'
3843	! do in=2,nl-1
3844	! write(1) vc(in), (c(in,ir) , ir=2,nl-1 )
3845	! es mas importante la precision que ocupar mucho espacio asi que
3846	! escribiremos las matrices en doble precision y por tanto en
3847	! [lib]readc_mz4.for no hay que reconvertirlas a doble precision.
3848	! ch is stored in single prec. to save storage space.
3849	! end do
3850
3851	! write(2) dummy
3852	! write(2)' format: (cfup(n,r),r=1,nl), n=1,nl)'
3853	! do in=1,nl
3854	! write(2) ( cup(in,ir) , ir=1,nl )
3855	! end do
3856
3857	! write(3) dummy
3858	! write(3)' format: (cfdw(n,r),r=1,nl), n=1,nl)'
3859	! do in=1,nl
3860	! write(3) (cdw(in,ir) , ir=1,nl )
3861	! end do
3862
3863	! if(ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5
3864	! @ .or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9) then
3865	! write (*,'(1x,30hcurtis matrix written out in: ,a50)' )
3866	! @ dircurtis//'cfl'//isotcode//dn//ibcode1//codmatrx//'.dat'
3867	! else
3868	! write (*,'(1x,30hcurtis matrix written out in: ,a50)' )
3869	! @ dircurtis//'cfl'//isotcode//dn//ibcode2//codmatrx//'.dat'
3870	! endif
3871
3872	else
3873
3874	! write (,) ' no curtis matrix output file ', char(10)
3875
3876	end if
3877
3878
3879	c end
3880	return
3881	end
3882
3883
3884
3885
3886
3887	c***********************************************************************
3888	c cm15um_hb_simple
3889	c***********************************************************************
3890
3891	subroutine cm15um_hb_simple (ig,icurt)
3892
3893	c computing the curtix matrixes for the 15 um hot bands
3894	c (las de las bandas fudnamentales las calcula cm15um_fb)
3895
3896	c jan 98 malv version de mod3/cm_15um.f para mz1d
3897	c jul 2011 malv+fgg adapted to LMD-MGCM
3898	c***********************************************************************
3899
3900	implicit none
3901
3902	!!!!!!!!!!!!!!!!!!!!!!!
3903	! common variables & constants
3904
3905	include 'nlte_paramdef.h'
3906	include 'nlte_commons.h'
3907
3908	!!!!!!!!!!!!!!!!!!!!!!!
3909	! arguments
3910
3911	integer ig ! ADDED FOR TRACEBACK
3912	integer icurt ! icurt=0,1,2
3913	! new calculations? (see caa.f heads)
3914
3915	!!!!!!!!!!!!!!!!!!!!!!!
3916	! local variables
3917
3918	real*4 cdummy(nl,nl), csngl(nl,nl)
3919
3920	real*8 cax1(nl,nl), cax2(nl,nl), cax3(nl,nl)
3921	real*8 v1(nl), v2(nl), v3(nl), cm_factor, vc_factor
3922
3923	integer itauout,icfout,itableout, interpol,ismooth, isngldble
3924	integer i,j,ik,ist,isot,ib,itt
3925
3926	!character bandcode*2
3927	character isotcode*2
3928	character codmatrx_hot*5
3929
3930	!!!!!!!!!!!!!!!!!!!!!!!
3931	! external functions
3932
3933	external bandid
3934	character*2 bandid
3935
3936	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3937	! subroutines called:
3938	! mz4sub, dmzout, readc_mz4, readcupdw, mztf
3939
3940	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3941	! formatos
3942	132 format(i2)
3943
3944	************************************************************************
3945	************************************************************************
3946
3947	call zerom (c121,nl)
3948
3949	call zerov (vc121,nl)
3950
3951	call zerom (cup121,nl)
3952	call zerom (cdw121,nl)
3953
3954	call zerov (taugr121,nl)
3955
3956
3957	itauout = 0 ! =1 --> with output of tau
3958	icfout = 0 ! =1 --> with output of cf
3959	itableout = 0 ! =1 --> with output of table of taus
3960	isngldble = 1 ! =1 --> dble precission
3961
3962	codmatrx_hot=' '
3963	if (icurt.eq.2) then
3964	icfout=1
3965	elseif (icurt.eq.0) then
3966	write (*,'(a,a$)')
3967	@ ' hot bands. code for old matrixes (5 chars): '
3968	read (*,'(a)') codmatrx_hot
3969	endif
3970
3971	fileroot = 'cfl'
3972
3973	! ====================== curtis matrixes for fh bands ==================
3974
3975
3976	! una piedra en el camino ...
3977	! write (,) ' cm15um_hb/1 '
3978
3979	ccc
3980	if ( input_cza.ge.1 ) then
3981	ccc
3982
3983	if (icurt.eq.2) then
3984	write (*,'(a,a$)')
3985	@ ' new calculation of curt. mat. for fh bands.',
3986	@ ' code for new matrixes : '
3987	read (*,'(a)') codmatrx_hot
3988	elseif (icurt.eq.0) then
3989	write (*,'(a,a$)')
3990	@ ' reading in curt. mat. for fh bands.',
3991	@ ' code for old matrixes : '
3992	read (*,'(a)') codmatrx_hot
3993	else
3994	! write (*,'(a)')
3995	! @ ' new calculation of curt. mat. for fh bands.'
3996	end if
3997
3998	! fh bands for the 626 isotope ================================-
3999
4000	ist = 1
4001	isot = 26
4002	! encode (2,132,isotcode) isot
4003	write (isotcode,132) isot
4004
4005	do 11, ik=1,3
4006
4007	ib=ik+1
4008
4009	if (icurt.gt.0) then
4010	call zero3m (cax1,cax2,cax3,nl)
4011	! una piedra en el camino ...
4012	!write (,) ' cm15um_hb/11 '
4013	!write (,) ' ib, ist, irw, imu =', ib, ist, irw_mztf, imu
4014	call mztf ( ig,cax1,cax2,cax3,v1,v2, ib,ist,irw_mztf,imu,
4015	@ itauout,icfout,itableout)
4016	! else
4017	! bandcode = bandid(ib)
4018	! filend=isotcode//dn//bandcode//codmatrx_hot
4019	!! write (,) char(9), fileroot//filend
4020	! call zero3m (cax1,cax2,cax3,nl)
4021	! call readcud_mz1d ( cax1,cax2,cax3, v1, v2,
4022	! @ fileroot,filend, csngl, nl,nan,0,isngldble)
4023	end if
4024
4025	c calculating the total c121(n,r) matrix for the first hot band
4026	do i=1,nl
4027
4028	if ( ib .eq. 4 ) then
4029	! write (,) ' '
4030	! write (,) i, ' ib,ist, altura :', ib,ist, zl(i)
4031	endif
4032
4033	! if ( v1(i) .le. 1.d-99 ) v1(i) = 0.0d0
4034	! if ( v2(i) .le. 1.d-99 ) v2(i) = 0.0d0
4035
4036
4037	if(ik.eq.1)then
4038	cm_factor = (dble(618.03/667.75))*2.d0
4039	@ exp( dble(ee*(667.75-618.03)/t(i)) )
4040	vc_factor = dble(667.75/618.03)
4041	elseif(ik.eq.2)then
4042	cm_factor = 1.d0
4043	vc_factor = 1.d0
4044	elseif(ik.eq.3)then
4045	cm_factor = ( dble(720.806/667.75) )*2.d0
4046	@ exp( dble(ee*(667.75-720.806)/t(i)) )
4047	vc_factor = dble(667.75/720.806)
4048	end if
4049	do j=1,nl
4050	! if (cax1(i,j) .le. 1.d-99 ) cax1(i,j) = 0.0d0
4051	! if (cax2(i,j) .le. 1.d-99 ) cax2(i,j) = 0.0d0
4052	! if (cax3(i,j) .le. 1.d-99 ) cax3(i,j) = 0.0d0
4053	c121(i,j) = c121(i,j) + cax1(i,j) * cm_factor
4054	cup121(i,j) = cup121(i,j) + cax2(i,j) * cm_factor
4055	cdw121(i,j) = cdw121(i,j) + cax3(i,j) * cm_factor
4056	end do
4057
4058	! write (,) ' i =', i
4059	! write (,) ' vc_factor =', vc_factor
4060	! write (,) ' v1 =', v1(i)
4061	! write (,) ' v2 =', v2(i)
4062	! write (,) vc121(i), taugr121(i)
4063	! write (,) v1(i) * vc_factor
4064	! write (,) vc121(i) + v1(i) * vc_factor
4065
4066	vc121(i) = vc121(i) + v1(i) * vc_factor
4067
4068
4069	! write (,) v2(i) * vc_factor
4070	! write (,) taugr121(i) + v2(i) * vc_factor
4071
4072	taugr121(i) = taugr121(i) + v2(i) * vc_factor
4073
4074	end do
4075	11 continue
4076
4077	ccc
4078	end if
4079	ccc
4080
4081
4082	return
4083	end
4084
4085
4086
4087
4088

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

Download in other formats:

Original Format