Context Navigation

nlte_calc.F @ 690

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

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

Download in other formats:

Original Format