Context Navigation

mztud.F @ 414

Last change on this file since 414 was 414, checked in by aslmd, 13 years ago

LMDZ.MARS : NEW NLTE MODEL FROM GRANADA AMIGOS

23/11/11 == FGG + MALV

New parameterization of the NLTE 15 micron cooling. The old parameterization is kept as an option, including or not variable atomic oxygen concentration. A new flag is introduced in callphys.def, nltemodel, to select which parameterization wants to be used (new one, old one with variable [O], or old one with fixed [O], see below). Includes many new subroutines and commons in phymars. Some existing routines are also modified:

-physiq.F. Call to the new subroutine NLTE_leedat in first call. Call to nltecool modified to allow for variable atomic oxygen. Depending on the value of nltemodel, the new subroutine NLTEdlvr09_TCOOL is called instead of nltecool.

-inifis.F. Reading of nltemodel is added.

-callkeys.h Declaration of nltemodel is added.

The following lines should be added to callphys.def (ideally after setting callnlte):

# NLTE 15um scheme to use.
# 0-> Old scheme, static oxygen
# 1-> Old scheme, dynamic oxygen
# 2-> New scheme
nltemodel = 2

A new directory, NLTEDAT, has to be included in datagcm.

Improvements into NLTE NIR heating parameterization to take into account variability with O/CO2 ratio and SZA. A new subroutine, NIR_leedat.F, and a new common, NIRdata.h, are included. A new flag, nircorr, is added in callphys.def, to include or not these corrections. The following files are modified:

-nirco2abs.F: nq and pq are added in the arguments. The corrections factors are interpolated to the GCM grid and included in the clculation. A new subroutine, interpnir, is added at the end of the file.

-physiq.F: Call to NIR_leedat added at first call. Modified call to nirco2abs

-inifis: Reading new flag nircorr.

-callkeys.h: Declaration of nircorr.

The following lines have to be added to callphys.def (ideally after callnirco2):

# NIR NLTE correction for variable SZA and O/CO2?
# matters only if callnirco2=T
# 0-> no correction
# 1-> include correction
nircorr=1

A new data file, NIRcorrection_feb2011.dat, has to be included in datagcm.

Small changes to the molecular diffusion scheme to fix the number of species considered, to avoid problems when compiling with more than 15 tracers (for example, when CH4 is included). Modified subroutines: aeronomars/moldiff.F and aeronomars/moldiffcoeff.F

File size: 24.9 KB

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

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

Download in other formats:

Original Format