Context Navigation

mztud.F @ 455

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

LMDZ.MARS : NEW NLTE MODEL FROM GRANADA AMIGOS

23/11/11 == FGG + MALV

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

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

-inifis.F. Reading of nltemodel is added.

-callkeys.h Declaration of nltemodel is added.

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

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

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

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

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

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

-inifis: Reading new flag nircorr.

-callkeys.h: Declaration of nircorr.

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

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

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

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

File size: 24.9 KB

Rev	Line
[414]	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