Context Navigation

nlte_calc.F @ 706

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

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

Download in other formats:

Original Format