Context Navigation

nlte_tcool.F @ 1009

Last change on this file since 1009 was 759, checked in by emillour, 12 years ago

Mars GCM:

commented out some messages in nlte_tcool.F (which were informations

useful and meaningful during code development only)

FGG

File size: 44.6 KB

Rev	Line
[757]	1	c**********************************************************************
	2	c
	3	c Contains the following old 1-d model subroutines:
	4	c
	5	c -NLTEdlvr11_TCOOL_03
	6	c -NLTEdlvr11_CZALU_03
	7	c -NLTEdlvr11_FB626CTS_03
	8	c
	9	c
	10	c
	11	c * Old NLTEdlvr11_TCOOL_02 *
	12	c
[498]	13	c***********************************************************************
	14
	15	c***********************************************************************
[757]	16
	17	subroutine nlte_tcool(ngridgcm,n_gcm,
	18	$ p_gcm, t_gcm, z_gcm,
	19	$ co2vmr_gcm, n2vmr_gcm, covmr_gcm, o3pvmr_gcm,
	20	$ q15umco2_gcm )
	21
	22	c***********************************************************************
	23
	24	implicit none
	25
[498]	26	include "dimensions.h"
	27	include "dimphys.h"
	28	include 'nlte_paramdef.h'
	29	include 'nlte_commons.h'
	30	include "chimiedata.h"
	31	include "conc.h"
	32
[757]	33
	34	c Arguments
[498]	35	integer n_gcm,ngridgcm
	36	real p_gcm(ngridgcm,n_gcm), t_gcm(ngridgcm,n_gcm)
[757]	37	real z_gcm(ngridgcm,n_gcm)
[498]	38	real co2vmr_gcm(ngridgcm,n_gcm), n2vmr_gcm(ngridgcm,n_gcm)
	39	real covmr_gcm(ngridgcm,n_gcm), o3pvmr_gcm(ngridgcm,n_gcm)
	40	real q15umco2_gcm(ngridgcm,n_gcm)
[757]	41	! real auxgcm(n_gcm)
	42	real*8 auxgcmd(n_gcm), aux2gcmd(n_gcm)
	43	real zmin_gcm
	44	integer ierr
	45	real*8 varerr
[498]	46
	47
[757]	48
	49	c local variables and constants
	50	integer i,ig,l, indice, nl_cts_real, nzy_cts_real
	51	real*8 q15umco2_nltot(nltot), zld(nltot)
	52	real*8 hr110CTS(nl_cts)
	53	real xx,factor
	54
[498]	55	real p_ig(n_gcm),z_ig(n_gcm)
	56	real t_ig(n_gcm)
	57	real co2_ig(n_gcm),n2_ig(n_gcm),co_ig(n_gcm),o3p_ig(n_gcm)
	58	real mmean_ig(n_gcm),cpnew_ig(n_gcm)
	59
	60
[757]	61	c***************
	62	c***************
	63
[498]	64	do ig=1,ngridgcm
[757]	65	ierr = 0
	66	nl_cts_real = 0
	67	nzy_cts_real = 0
[498]	68	do l=1,n_gcm
	69	p_ig(l)=p_gcm(ig,l)
	70	t_ig(l)=t_gcm(ig,l)
	71	co2_ig(l)=co2vmr_gcm(ig,l)
	72	n2_ig(l)=n2vmr_gcm(ig,l)
	73	o3p_ig(l)=o3pvmr_gcm(ig,l)
	74	co_ig(l)=covmr_gcm(ig,l)
	75	z_ig(l)=z_gcm(ig,l)/1000.
	76	mmean_ig(l)=mmean(ig,l)
	77	cpnew_ig(l)=cpnew(ig,l)
	78	enddo
	79
[757]	80	! From GCM's grid to NLTE's grid
	81	call NLTEdlvr11_ZGRID_02 (n_gcm,
	82	$ p_ig, t_ig, z_ig,
	83	$ co2_ig, n2_ig, co_ig, o3p_ig,
	84	$ mmean_ig,cpnew_ig,
	85	$ nl_cts_real, nzy_cts_real )
[498]	86
[757]	87
	88	! Isotopic Tstar & VC at the NLTE grid
	89	call interdp_ESCTVCISO
	90
	91	! Tstar para NLTE-CTS
	92	call MZESC110 ( nl_cts_real, nzy_cts_real )
	93
	94	! 626FB C.M.
[498]	95	call leetvt
[757]	96	c110(1:nl,1:nl)=0.d0
	97	! call zerom (c110, nl)
	98	call zero2v (vc110,taustar11, nl)
	99	call MZTUD110 ( ierr, varerr )
	100	if (ierr .gt. 0) goto 900
[498]	101
[757]	102	input_cza = 0
	103	call NLTEdlvr11_CZALU
	104
	105	input_cza = 1
	106	call NLTEdlvr11_CZALU
	107
	108	! call NLTEdlvr11_FB626CTS
	109	! Falta un merging del hr110CTS con el HR110
	110
	111
	112	! ! Interpolation of Tstar11(nl) to the GCM grid (será auxgcm)
	113	! ! solo entre jlowerboundary y jtopboundary (la extension del NLTE
	114	! ! model)
	115	! call interhuntlimits( auxgcm, p_gcm,n_gcm,
	116	! @ jlowerboundary,jtopboundary,
	117	! @ taustar11, pl, nl, 3 )
	118	! ! Mejor inter+extra polacion de Tstar11(nl) to the GCM grid
	119	! call TSTAR11_extension (n_gcm, p_gcm, auxgcm )
	120
	121	! NLTE-CTS
	122	call NLTEdlvr11_FB626CTS ( hr110CTS , nl_cts_real )
	123
	124
	125
	126	! total TCR
	127	do i = 1, nl
	128	q15umco2_nltot(i) =hr110(i) + hr210(i) + hr310(i) + hr410(i)
[498]	129	@ + hr121(i)
[757]	130	enddo
	131
[498]	132
[757]	133	! Merging con / actualizacion del HR_total
	134	! Eliminamos el ultimo pto de hrTotal, y en el penultimo
	135	! (que coincide con i=1 en el grid nl_cts)
	136	! hacemos la media entre hrTotal y hr110CTS :
	137	i=nl-1
	138	q15umco2_nltot(i) = 0.5*( q15umco2_nltot(i) + hr110CTS(1) )
	139	do i=2,nl_cts_real
	140	indice = (nl-2) + i
	141	q15umco2_nltot(indice) = hr110CTS(i)
[498]	142	enddo
[757]	143	do i=nl_cts_real+1,nl_cts
	144	indice = (nl-2) + i
	145	q15umco2_nltot(indice) = 0.0d0
	146	enddo
	147
	148	! Interpol to original Pgrid
	149	!
	150	! Primero, la parte conocida ([1,nl_cts_real])
	151	do i=1,nl
	152	zld(i) = - dble ( alog(pl(i)) )
	153	!write (,) i, zld(i), q15umco2_nltot(i)
	154	enddo
	155	do i=3,nl_cts_real
	156	indice = (nl-2) + i
	157	zld(indice) = - dble ( alog(pl_cts(i)) )
	158	!write (,) indice, zld(indice), q15umco2_nltot(indice)
	159	enddo
	160	! En caso que nl_cts_real < nl_cts , extrapolo el grid alegremente
	161	factor = pl_cts(nl_cts_real)/pl_cts(nl_cts_real-1)
	162	xx = pl_cts(nl_cts_real)
	163	do i=nl_cts_real+1,nl_cts
	164	indice = (nl-2) + i
	165	xx = xx * factor
	166	zld(indice) = - dble ( alog(xx) )
	167	enddo
	168
	169	do i=1,n_gcm
	170	auxgcmd(i) = - dble( alog(p_gcm(ig,i)) )
	171	enddo
	172	! call zerov( aux2gcmd, n_gcm )
	173	aux2gcmd(1:n_gcm)=0.d0
	174	call interdp_limits
	175	$ ( aux2gcmd, auxgcmd, n_gcm, jlowerboundary,jtopCTS,
	176	$ q15umco2_nltot, zld, nltot, 1, nltot,
	177	$ 1 )
	178
	179	! Smoothing
	180	call suaviza ( aux2gcmd, n_gcm, 1, auxgcmd )
	181
	182	do i=1,n_gcm
	183	q15umco2_gcm(ig,i) = sngl( aux2gcmd(i) )
	184	enddo
	185
[498]	186	enddo
[757]	187	c end subroutine
[498]	188	return
	189
[757]	190	c Error messages
	191	900 write (,) ' ERROR in MZTUD (banda 110). ierr=',ierr
	192	write (,) ' VAR available : ', varerr
	193	return
[498]	194
[757]	195	901 write (,) ' ERROR in MZTVC for vc210. ierr=',ierr
	196	write (,) ' VAR available : ', varerr
	197	return
[498]	198
[757]	199	902 write (,) ' ERROR in MZTVC for vc310. ierr=',ierr
	200	write (,) ' VAR available : ', varerr
	201	return
	202
	203	903 write (,) ' ERROR in MZTVC for vc410. ierr=',ierr
	204	write (,) ' VAR available : ', varerr
	205	return
	206
	207	904 write (,) ' ERROR in mzescape_fb ierr=',ierr
	208	write (,) ' VAR available : ', varerr
	209	return
	210
	211	930 write (,) ' ERROR in mztvc3iso. Temp is NaN'
	212	write (,) ' ierr , VAR =', ierr, varerr
	213	if (ierr.eq.30) write (,) ' During computation of VC210.'
	214	if (ierr.eq.31) write (,) ' During computation of VC310.'
	215	if (ierr.eq.32) write (,) ' During computation of VC410.'
	216	return
	217
	218	c end subroutine
	219	end
	220
	221
[498]	222	c***********************************************************************
	223
[757]	224	subroutine NLTEdlvr11_ZGRID_02 (n_gcm,
	225	$ p_gcm, t_gcm, z_gcm, co2vmr_gcm, n2vmr_gcm,
	226	$ covmr_gcm, o3pvmr_gcm, mmean_gcm,cpnew_gcm,
	227	$ nl_cts_real, nzy_cts_real )
[498]	228
	229	c***********************************************************************
[757]	230
	231	implicit none
	232
[498]	233	include 'nlte_paramdef.h'
	234	include 'nlte_commons.h'
	235
[757]	236	c Arguments
	237	integer n_gcm ! I
	238	real p_gcm(n_gcm), t_gcm(n_gcm) ! I
	239	real co2vmr_gcm(n_gcm), n2vmr_gcm(n_gcm) ! I
	240	real covmr_gcm(n_gcm), o3pvmr_gcm(n_gcm) ! I
	241	real z_gcm(n_gcm) ! I
	242	real mmean_gcm(n_gcm) ! I
	243	real cpnew_gcm(n_gcm) ! I
	244	integer nl_cts_real, nzy_cts_real ! O
	245
	246	c local variables
	247	integer i, iz
	248	real distancia, meanm, gz, Hkm
	249	real zmin, zmax
[498]	250	real mmean_nlte(n_gcm),cpnew_nlte(n_gcm)
	251
[757]	252	c functions
	253	external hrkday_convert
	254	real hrkday_convert
[498]	255
[757]	256	c***********************************************************************
[498]	257
[757]	258	! Define el working grid para MZ1D (NL, ZL, ZMIN)
	259	! y otro mas fino para M.Curtis (NZ, ZX, ZXMIN = ZMIN)
	260	! Tambien el working grid para MZESC110 (NL_cts, ZL_cts, ZMIN_cts=??)
	261	! Para ello hace falta una z de ref del GCM, que voy a suponer la inferior
[498]	262
[757]	263	! Primero, construimos escala z_gcm
[498]	264
[757]	265	! z_gcm(1) = zmin_gcm ! [km]
[498]	266
[757]	267	! do iz = 2, n_gcm
	268	! meanm = ( co2vmr_gcm(iz)44. + o3pvmr_gcm(iz)16.
	269	! @ + n2vmr_gcm(iz)28. + covmr_gcm(iz)28. )
	270	! meanm = meanm / n_avog
	271	! distancia = ( radio + z_gcm(iz-1) )*1.e5
	272	! gz = gg * masa / ( distancia * distancia )
	273	! Hkm = 0.5( t_gcm(iz)+t_gcm(iz-1) ) / ( meanm gz )
	274	! Hkm = kboltzman * Hkm *1e-5 ! [km]
	275	! z_gcm(iz) = z_gcm(iz-1) - Hkm * log( p_gcm(iz)/p_gcm(iz-1) )
	276	! enddo
[498]	277
[757]	278	! Segundo, definimos los límites de los 2 modelos de NLTE.
	279	! NLTE model completo: indices [jlowerboundary,jtopboundary]
	280	! NLTE CTS : indices [jbotCTS,jtopCTS] donde jbotCTS = jtopboundary-2
[498]	281
[757]	282	!!!!!!!!!Primero el NLTE completo !!!!!!!!
[498]	283
[757]	284	! Bottom boundary for NLTE model :
	285	! Pbot_atm = 2e-2 mb = 1.974e-5 atm , lnp(nb)=9.9 (see mz1d.par)
	286	jlowerboundary = 1
[498]	287	do while ( p_gcm(jlowerboundary) .gt. Pbottom_atm )
	288	jlowerboundary = jlowerboundary + 1
[757]	289	if (jlowerboundary .gt. n_gcm) then
	290	write (,) 'Error in lower boundary pressure.'
	291	write (,) ' p_gcm too low or wrong. '
	292	write (,) ' p_gcm, Pbottom_atm =',
	293	$ p_gcm(n_gcm), Pbottom_atm
	294	stop ' Check input value "p_gcm" or modify "Pbottom_atm" '
	295	endif
[498]	296	enddo
	297
[757]	298	! Top boundary for NLTE model :
	299	! Ptop_atm = 1e-9 atm (see mz1d.par)
	300	jtopboundary = jlowerboundary
	301	do while ( p_gcm(jtopboundary) .gt. Ptop_atm )
[498]	302	jtopboundary = jtopboundary + 1
[757]	303	if (jtopboundary .gt. n_gcm) then
	304	write (,) '!!!!!!!! Warning in top boundary pressure. '
	305	write (,) ' Ptop_atm too high for p_gcm. '
	306	write (,) ' p_gcm, Ptop_atm =',
	307	$ p_gcm(n_gcm), Ptop_atm
	308	write (,) '!!!!!!!! NLTE upper boundary modified '//
	309	$ ' to match p_gcm'
	310	jtopboundary=n_gcm
	311	goto 5000
	312	endif
[498]	313	enddo
[757]	314	5000 continue
	315
	316	! Grid steps
	317
	318	zmin = z_gcm(jlowerboundary)
[498]	319	zmax = z_gcm(jtopboundary)
[757]	320	deltaz = (zmax-zmin) / (nl-1)
	321	do i=1,nl
[498]	322	zl(i) = zmin + (i-1) * deltaz
	323	enddo
	324
[757]	325
	326	! Creamos el perfil del NLTE modelo completo interpolando
	327
	328	call interhunt ( pl,zl,nl, p_gcm,z_gcm,n_gcm, 2) ! [atm]
	329	call interhunt5veces
	330	$ ( t, co2vmr, n2vmr, covmr, o3pvmr,
	331	$ zl, nl,
	332	$ t_gcm, co2vmr_gcm, n2vmr_gcm, covmr_gcm, o3pvmr_gcm,
	333	$ z_gcm, n_gcm,
	334	$ 1 )
	335	call interhunt ( mmean_nlte,zl,nl,mmean_gcm,z_gcm,n_gcm,1)
	336	call interhunt ( cpnew_nlte,zl,nl,cpnew_gcm,z_gcm,n_gcm,1)
	337
[498]	338	do i = 1, nl
[757]	339	nt(i) = 7.339e+21 * pl(i) / t(i) ! --> [cm-3]
[498]	340	co2(i) = nt(i) * co2vmr(i)
	341	n2(i) = nt(i) * n2vmr(i)
	342	co(i) = nt(i) * covmr(i)
	343	o3p(i) = nt(i) * o3pvmr(i)
[757]	344	! hrkday_factor(i) = hrkday_convert( t(i),
	345	! $ co2vmr(i), o3pvmr(i), n2vmr(i), covmr(i) )
[695]	346	hrkday_factor(i) = hrkday_convert(mmean_nlte(i)
[757]	347	& ,cpnew_nlte(i))
	348	enddo
	349
	350	! Comprobar que las temps no se salen del grid del histograma
[498]	351
[757]	352	do i=1,nl
	353	if (t(i) .gt. 400.0) then
	354	write (,) '!!!! WARNING Temp higher than Histogram.'
	355	write (,) ' Histogram will be extrapolated. '
	356	write (,) ' i, t(i), pl(i) =', i, t(i), pl(i)
	357	endif
	358	if (t(i) .lt. 50.0) then
	359	write (,) '!!!! WARNING Temp lower than Histogram.'
	360	write (,) ' Histogram will be extrapolated. '
	361	write (,) ' i, t(i), pl(i) =', i, t(i), pl(i)
	362	endif
[498]	363	enddo
	364
[757]	365	! Fine grid for transmittance calculations
[498]	366
[757]	367	zmin = z_gcm(jlowerboundary)
	368	zmax = z_gcm(jtopboundary)
	369	deltazy = (zmax-zmin) / (nzy-1)
	370	do i=1,nzy
	371	zy(i) = zmin + (i-1) * deltazy
	372	enddo
	373	call interhunt ( py,zy,nzy, p_gcm,z_gcm,n_gcm, 2) ! [atm]
	374	call interhunt2veces ( ty,co2y, zy,nzy,
	375	$ t_gcm,co2vmr_gcm, z_gcm,n_gcm, 1)
[498]	376
[757]	377	do i=1,nzy
	378	nty(i) = 7.339e+21 * py(i) / ty(i) ! --> [cm-3]
[498]	379	co2y(i) = co2y(i) * nty(i)
	380	enddo
	381
	382
[757]	383	!!!!!!!!!Segundo, el NLTE - CTS !!!!!!!!
[498]	384
[757]	385	! Grid steps
	386	deltaz_cts = deltaz
	387	zl_cts(1) = zl(nl-1)
	388	nl_cts_real = 1
	389	do i=2,nl_cts
	390	zl_cts(i) = zl_cts(1) + (i-1)*deltaz_cts
	391	if (zl_cts(i) .gt. z_gcm(n_gcm)) then
[759]	392	! write (,) '!!!!!!!! Warning in top CTS layers. '
	393	! write (,) ' zl_Cts too high for z_gcm. '
	394	! write (,) ' z_gcm, zl_cts(i), i =',
	395	! $ z_gcm(n_gcm), zl_cts(i), i
	396	! write (,) '!!!!!!!! NLTE-CTS upper boundary modified '//
	397	! $ ' to match z_gcm'
[757]	398	nl_cts_real=i-1
[759]	399	! write (,) ' Original,Real NL_CTS=', nl_cts,nl_cts_real
[757]	400	goto 6000
	401	endif
	402	enddo
	403	nl_cts_real = nl_cts
	404	6000 continue
	405
	406	! Creamos perfil por interpolacion
[498]	407
[757]	408	call interhuntlimits ( pl_cts,zl_cts,nl_cts, 1,nl_cts_real,
	409	$ p_gcm,z_gcm,n_gcm, 2)
	410	call interhuntlimits5veces
	411	$ ( t_cts, co2vmr_cts, n2vmr_cts, covmr_cts, o3pvmr_cts,
	412	$ zl_cts, nl_cts,
	413	$ 1,nl_cts_real,
	414	$ t_gcm, co2vmr_gcm, n2vmr_gcm, covmr_gcm, o3pvmr_gcm,
	415	$ z_gcm, n_gcm,
	416	$ 1 )
	417	call interhuntlimits( cpnew_cts,zl_cts,nl_cts,1,nl_cts_real,
	418	$ cpnew_gcm,z_gcm,n_gcm, 1)
	419	call interhuntlimits( mmean_cts,zl_cts,nl_cts,1,nl_cts_real,
	420	$ mmean_gcm,z_gcm,n_gcm, 1)
[498]	421
[757]	422	do i = 1, nl_cts_real
	423	nt_cts(i) = 7.339e+21 * pl_cts(i) / t_cts(i) ! --> [cm-3]
	424	co2_cts(i) = nt_cts(i) * co2vmr_cts(i)
	425	n2_cts(i) = nt_cts(i) * n2vmr_cts(i)
	426	co_cts(i) = nt_cts(i) * covmr_cts(i)
	427	o3p_cts(i) = nt_cts(i) * o3pvmr_cts(i)
	428	hrkday_factor_cts(i) = hrkday_convert( mmean_cts(i)
	429	& ,cpnew_cts(i) )
	430	enddo
	431
	432	! Comprobar que las temps no se salen del grid del histograma
	433	do i=1,nl_cts_real
	434	if (t_cts(i) .gt. 400.0) then
	435	write (,) '!!!! WARNING Temp higher than Histogram.'
	436	write (,) ' ZGRID: Histogram will be extrapolated. '
	437	write (,) ' i, t(i), pl(i) =', i, t_cts(i), pl_cts(i)
	438	endif
	439	if (t_cts(i) .lt. 50.0) then
	440	write (,) '!!!! WARNING Temp lower than Histogram.'
	441	write (,) ' ZGRID: Histogram will be extrapolated. '
	442	write (,) ' i, t(i), pl(i) =', i, t_cts(i), pl_cts(i)
	443	endif
	444	enddo
	445
	446	! Calculo del indice maximo del GCM hasta donde llega el NLTE-CTS
	447	jtopCTS = jtopboundary
	448	do while ( p_gcm(jtopCTS) .gt. pl_cts(nl_cts_real) )
	449	jtopCTS = jtopCTS + 1
	450	if (jtopCTS .gt. n_gcm) then
	451	write (,) '!!!!!!!! Warning in top boundary pressure. '
	452	write (,) ' Ptop_NLTECTS too high for p_gcm. '
	453	write (,) ' p_gcm, Ptop_NLTECTS =',
	454	$ p_gcm(n_gcm), pl_cts(nl_cts_real)
	455	write (,) '!!!!!!!! NLTE-CTS upper boundary modified '//
	456	$ ' to match p_gcm'
	457	jtopCTS=n_gcm
	458	goto 7000
	459	endif
	460	enddo
	461	7000 continue
	462
	463	! Fine grid for transmittance calculations
	464
	465	deltazy_cts = 0.25*deltaz_cts ! Comprobar el factor 4 en mz1d.par
	466	do i=1,nzy_cts
	467	zy_cts(i) = zl_cts(1) + (i-1) * deltazy_cts
	468	enddo
	469	nzy_cts_real = (nl_cts_real - 1)*4 + 1
	470	call interhuntlimits ( py_cts,zy_cts,nzy_cts, 1,nzy_cts_real,
	471	$ p_gcm, z_gcm, n_gcm, 2) ! [atm]
	472	call interhuntlimits2veces
	473	$ ( ty_cts,co2y_cts, zy_cts,nzy_cts, 1,nzy_cts_real,
	474	$ t_gcm,co2vmr_gcm, z_gcm,n_gcm, 1)
	475
	476	do i=1,nzy_cts_real
	477	nty_cts(i) = 7.339e+21 * py_cts(i) / ty_cts(i) ! --> [cm-3]
	478	co2y_cts(i) = co2y_cts(i) * nty_cts(i)
	479	enddo
	480
	481	! write (,) ' NL = ', NL
	482	! write (,) ' Original,Real NL_CTS=', nl_cts,nl_cts_real
	483	! write (,) ' Original,Real NZY_CTS =', nzy_cts,nzy_cts_real
	484
	485
	486
	487	c end
	488	return
	489	end
	490
	491
	492	c * Old NLTEdlvr11_CZALU_03 *
	493
	494	c**********************************************************************
	495
	496
	497	subroutine NLTEdlvr11_CZALU
	498
[498]	499	c***********************************************************************
	500
[757]	501	implicit none
	502
	503	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!common variables and constants
	504
[498]	505	include 'nlte_paramdef.h'
	506	include 'nlte_commons.h'
	507
[757]	508
	509	c local variables
	510
	511	! matrixes and vectors
	512
	513	real*8 e110(nl), e210(nl), e310(nl), e410(nl)
	514	real*8 e121(nl)
[498]	515	real*8 f1(nl,nl)
[757]	516
	517	real*8 cax1(nl,nl), cax2(nl,nl), cax3(nl,nl)
[498]	518	real*8 v1(nl), v2(nl), v3(nl)
[757]	519	real*8 alf11(nl,nl), alf12(nl,nl)
	520	real*8 alf21(nl,nl), alf31(nl,nl), alf41(nl,nl)
	521	real*8 a11(nl), a1112(nl,nl)
	522	real*8 a1121(nl,nl), a1131(nl,nl), a1141(nl,nl)
	523	real*8 a21(nl), a2131(nl,nl), a2141(nl,nl)
	524	real*8 a2111(nl,nl), a2112(nl,nl)
	525	real*8 a31(nl), a3121(nl,nl), a3141(nl,nl)
	526	real*8 a3111(nl,nl), a3112(nl,nl)
	527	real*8 a41(nl), a4121(nl,nl), a4131(nl,nl)
	528	real*8 a4111(nl,nl), a4112(nl,nl)
	529	real*8 a12(nl), a1211(nl,nl)
	530	real*8 a1221(nl,nl), a1231(nl,nl), a1241(nl,nl)
[498]	531
[757]	532	real*8 aalf11(nl,nl),aalf21(nl,nl),
	533	@ aalf31(nl,nl),aalf41(nl,nl)
	534	real*8 aa11(nl), aa1121(nl,nl), aa1131(nl,nl), aa1141(nl,nl)
	535	real*8 aa21(nl), aa2111(nl,nl), aa2131(nl,nl), aa2141(nl,nl)
	536	real*8 aa31(nl), aa3111(nl,nl), aa3121(nl,nl), aa3141(nl,nl)
	537	real*8 aa41(nl), aa4111(nl,nl), aa4121(nl,nl), aa4131(nl,nl)
	538	real*8 aa1211(nl,nl),aa1221(nl,nl),
	539	@ aa1231(nl,nl),aa1241(nl,nl)
	540	real*8 aa1112(nl,nl),aa2112(nl,nl),
	541	@ aa3112(nl,nl),aa4112(nl,nl)
	542
	543	real*8 aaalf11(nl,nl), aaalf31(nl,nl), aaalf41(nl,nl)
	544	real*8 aaa11(nl),aaa1131(nl,nl),aaa1141(nl,nl)
	545	real*8 aaa31(nl),aaa3111(nl,nl),aaa3141(nl,nl)
	546	real*8 aaa41(nl),aaa4111(nl,nl),aaa4131(nl,nl)
	547
	548	real*8 aaaalf11(nl,nl),aaaalf41(nl,nl)
	549	real*8 aaaa11(nl),aaaa1141(nl,nl)
	550	real*8 aaaa41(nl),aaaa4111(nl,nl)
	551
	552
	553	! populations
	554	real*8 n10(nl), n11(nl), n12(nl)
[498]	555	real*8 n20(nl), n21(nl)
	556	real*8 n30(nl), n31(nl)
	557	real*8 n40(nl), n41(nl)
	558
[757]	559	! productions and loses
	560	real*8 d19b1,d19c1
	561	real*8 d19bp1,d19cp1
	562	real*8 d19c2
	563	real*8 d19cp2
	564	real*8 d19c3
	565	real*8 d19cp3
	566	real*8 d19c4
	567	real*8 d19cp4
[498]	568
[757]	569	real*8 l11, l12, l21, l31, l41
	570	real*8 p11, p12, p21, p31, p41
	571	real*8 p1112, p1211, p1221, p1231, p1241
	572	real*8 p1121, p1131, p1141
	573	real*8 p2111, p2112, p2131, p2141
	574	real*8 p3111, p3112, p3121, p3141
	575	real*8 p4111, p4112, p4121, p4131
[498]	576
[757]	577	real*8 pl11, pl12, pl21, pl31, pl41
[498]	578
[757]	579
	580	c local constants and indexes
	581
	582	real*8 co2t, o3pdbl, codble, n2dble
	583	real*8 a12_einst(nl)
	584	real*8 a21_einst(nl), a31_einst(nl), a41_einst(nl)
	585	real tsurf
	586
	587	integer i, isot
	588
	589	c external functions and subroutines
	590
	591	external planckdp
	592	real*8 planckdp
	593
	594
	595	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!start program
	596
	597
[498]	598	call zero4v( aa11, aa21, aa31, aa41, nl)
	599	call zero4m( aa1121, aa1131, aa1141, aalf11, nl)
	600	call zero4m( aa2111, aa2131, aa2141, aalf21, nl)
	601	call zero4m( aa3111, aa3121, aa3141, aalf31, nl)
	602	call zero4m( aa4111, aa4121, aa4131, aalf41, nl)
	603	call zero4m( aa1112, aa2112, aa3112, aa4112, nl)
	604	call zero4m( aa1211, aa1221, aa1231, aa1241, nl)
	605	call zero3v( aaa41, aaa31, aaa11, nl )
	606	call zero3m( aaa4111, aaa4131, aaalf41, nl)
	607	call zero3m( aaa3111, aaa3141, aaalf31, nl)
	608	call zero3m( aaa1131, aaa1141, aaalf11, nl)
	609	call zero2v( aaaa11, aaaa41, nl )
	610	call zero2m( aaaa1141, aaaalf11, nl)
	611	call zero2m( aaaa4111, aaaalf41, nl)
	612
[757]	613	call zero2v (vt11,vt12,nl)
	614	call zero3v (vt21,vt31,vt41,nl)
	615	call zero2v (hr110,hr121,nl)
	616	call zero3v (hr210,hr310,hr410,nl)
	617	call zero2v (sl110,sl121,nl)
	618	call zero3v (sl210,sl310,sl410,nl)
[498]	619
[757]	620	call zero4v (el11,el21,el31,el41,nl)
	621	call zero4v (e110,e210,e310,e410,nl)
	622	call zero2v (el12,e121,nl)
[498]	623
[757]	624	call zero3m (cax1,cax2,cax3,nl)
	625	f1(1:nl,1:nl)=0.d0
	626	! call zerom (f1,nl)
	627
	628	call zero3v (v1,v2,v3,nl)
	629
	630	call zero4m (alf11,alf21,alf31,alf41,nl)
	631	alf12(1:nl,1:nl)=0.d0
	632	! call zerom (alf12,nl)
	633	call zero2v (a11,a12,nl)
	634	call zero3v (a21,a31,a41,nl)
	635
	636	call zero3m (a1121,a1131,a1141,nl)
	637	a1112(1:nl,1:nl)=0.d0
	638	! call zerom (a1112,nl)
	639
	640	call zero3m (a1221,a1231,a1241,nl)
	641	a1211(1:nl,1:nl)=0.d0
	642	! call zerom (a1211,nl)
	643
	644	call zero2m (a2111,a2112,nl)
	645	call zero2m (a2131,a2141,nl)
	646	call zero2m (a3111,a3112,nl)
	647	call zero2m (a3121,a3141,nl)
	648	call zero2m (a4111,a4112,nl)
	649	call zero2m (a4121,a4131,nl)
	650
	651	call zero2v (n11,n12,nl)
	652	call zero3v (n21,n31,n41,nl)
	653
	654	nu11 = dble(nu(1,1))
	655	nu12 = dble(nu(1,2))
	656	nu121 = nu12-nu11
	657	nu21 = dble(nu(2,1))
	658	nu31 = dble(nu(3,1))
	659	nu41 = dble(nu(4,1))
	660
	661	c
	662	c
	663	do i=1,nl
	664	n10(i) = dble( co2(i) * imr(1) )
	665	n20(i) = dble( co2(i) * imr(2) )
	666	n30(i) = dble( co2(i) * imr(3) )
	667	n40(i) = dble( co2(i) * imr(4) )
	668	if ( input_cza.ge.1 ) then
	669	n11(i) = n10(i) 2.d0 exp( -ee*nu11/v626t1(i) )
	670	n21(i) = n20(i) 2.d0 exp( -ee*nu21/v628t1(i) )
	671	n31(i) = n30(i) 2.d0 exp( -ee*nu31/v636t1(i) )
	672	n41(i) = n40(i) 2.d0 exp( -ee*nu41/v627t1(i) )
	673	end if
	674	enddo
	675
	676	c
	677	c curtis matrix calculation
	678	c
	679	call zero3m (c210,c310,c410, nl)
	680
	681	if ( input_cza.ge.1 ) then
	682
	683	if (itt_cza.eq.15 ) then
	684
	685	call MZMC121
	686
	687	elseif (itt_cza.eq.13) then
	688
	689	! call zerom ( c121, nl )
	690	c121(1:nl,1:nl)=0.d0
	691	call MZESC121
	692	call MZTVC121
	693
[498]	694	endif
	695
	696	endif
	697
[757]	698	! Lower Boundary
	699	tsurf = t(1)
	700	do i=1,nl
	701	sl110(i) = vc110(i) * planckdp( tsurf, nu11 )
	702	sl210(i) = vc210(i) * planckdp( tsurf, nu21 )
	703	sl310(i) = vc310(i) * planckdp( tsurf, nu31 )
	704	sl410(i) = vc410(i) * planckdp( tsurf, nu41 )
	705	end do
	706	if (input_cza.ge.1) then
	707	do i=1,nl
	708	sl121(i) = vc121(i) * planckdp( tsurf, nu121 )
	709	end do
	710	endif
[498]	711
	712
[757]	713
[498]	714	do 4,i=nl,1,-1 !----------------------------------------------
	715
[757]	716	co2t = dble( co2(i) *(imr(1)+imr(3)+imr(2)+imr(4)) )
	717	o3pdbl = dble( o3p(i) )
	718	n2dble = dble( n2(i) )
	719	codble = dble ( co(i) )
	720
	721	call GETK_dlvr11 ( t(i) )
	722
	723	! V-T productions and losses V-T
	724
[498]	725	isot = 1
[757]	726	d19b1 = k19ba(isot)co2t + k19bb(isot)n2dble
	727	@ + k19bc(isot)*codble
	728	d19c1 = k19ca(isot)co2t + k19cb(isot)n2dble
	729	@ + k19cc(isot)*codble
	730	d19bp1 = k19bap(isot)co2t + k19bbp(isot)n2dble
	731	@ + k19bcp(isot)*codble
	732	d19cp1 = k19cap(isot)co2t + k19cbp(isot)n2dble
	733	@ + k19ccp(isot)*codble
[498]	734	isot = 2
[757]	735	d19c2 = k19ca(isot)co2t + k19cb(isot)n2dble
	736	@ + k19cc(isot)*codble
	737	d19cp2 = k19cap(isot)co2t + k19cbp(isot)n2dble
	738	@ + k19ccp(isot)*codble
[498]	739	isot = 3
[757]	740	d19c3 = k19ca(isot)co2t + k19cb(isot)n2dble
	741	@ + k19cc(isot)*codble
	742	d19cp3 = k19cap(isot)co2t + k19cbp(isot)n2dble
	743	@ + k19ccp(isot)*codble
[498]	744	isot = 4
[757]	745	d19c4 = k19ca(isot)co2t + k19cb(isot)n2dble
	746	@ + k19cc(isot)*codble
	747	d19cp4 = k19cap(isot)co2t + k19cbp(isot)n2dble
	748	@ + k19ccp(isot)*codble
	749	!
	750	l11 = d19c1 + k20c(1)*o3pdbl
	751	p11 = ( d19cp1 + k20cp(1)o3pdbl ) n10(i)
	752	l21 = d19c2 + k20c(2)*o3pdbl
	753	p21 = ( d19cp2 + k20cp(2)o3pdbl ) n20(i)
	754	l31 = d19c3 + k20c(3)*o3pdbl
	755	p31 = ( d19cp3 + k20cp(3)o3pdbl ) n30(i)
	756	l41 = d19c4 + k20c(4)*o3pdbl
	757	p41 = ( d19cp4 + k20cp(4)o3pdbl ) n40(i)
	758
	759	! Addition of V-V
	760
	761	l11 = l11 + k21cp(2)n20(i) + k21cp(3)n30(i)
	762	@ + k21cp(4)*n40(i)
	763	p1121 = k21c(2) * n10(i)
	764	p1131 = k21c(3) * n10(i)
	765	p1141 = k21c(4) * n10(i)
	766	!
	767	l21 = l21 + k21c(2)n10(i) + k23k21cn30(i) + k24k21c*n40(i)
	768	p2111 = k21cp(2) * n20(i)
	769	p2131 = k23k21cp * n20(i)
	770	p2141 = k24k21cp * n20(i)
	771	!
	772	l31 = l31 + k21c(3)n10(i) + k23k21cpn20(i) + k34k21c*n40(i)
	773	p3111 = k21cp(3)* n30(i)
	774	p3121 = k23k21c * n30(i)
	775	p3141 = k34k21cp* n30(i)
	776	!
	777	l41 = l41 + k21c(4)n10(i) + k24k21cpn20(i) + k34k21cp*n30(i)
	778	p4111 = k21cp(4)* n40(i)
	779	p4121 = k24k21c * n40(i)
	780	p4131 = k34k21c * n40(i)
	781
	782
	783	if ( input_cza.ge.1 ) then
	784
	785	l12 = d19b1
	786	@ + k20b(1)*o3pdbl
	787	@ + k21b(1)*n10(i)
	788	@ + k33c*( n20(i) + n30(i) + n40(i) )
	789	p12 = k21bp(1)n11(i) n11(i)
	790	p1211 = d19bp1 + k20bp(1)*o3pdbl
	791	p1221 = k33cp(2)*n11(i)
	792	p1231 = k33cp(3)*n11(i)
	793	p1241 = k33cp(4)*n11(i)
	794
	795	l11 = l11 + d19bp1
	796	@ + k20bp(1)*o3pdbl
	797	@ + 2.d0 * k21bp(1) * n11(i)
[498]	798	@ + k33cp(2)n21(i) + k33cp(3)n31(i) + k33cp(4)*n41(i)
[757]	799	p1112 = d19b1
	800	@ + k20b(1)*o3pdbl
	801	@ + 2.d0k21b(1)n10(i)
	802	@ + k33c*( n20(i) + n30(i) + n40(i) )
[498]	803
[757]	804	l21 = l21 + k33cp(2)*n11(i)
	805	p2112 = k33c*n20(i)
[498]	806
[757]	807	l31 = l31 + k33cp(3)*n11(i)
	808	p3112 = k33c*n30(i)
[498]	809
[757]	810	l41 = l41 + k33cp(4)*n11(i)
	811	p4112 = k33c*n40(i)
[498]	812
[757]	813	end if
	814
	815
	816	! For ITT=13,15
	817
	818	a21_einst(i) = a2_010_000 * 1.8d0 / 4.d0 * taustar21(i)
	819	a31_einst(i) = a3_010_000 * 1.8d0 / 4.d0 * taustar31(i)
	820	a41_einst(i) = a4_010_000 * 1.8d0 / 4.d0 * taustar41(i)
	821
	822	l21 = l21 + a21_einst(i)
	823	l31 = l31 + a31_einst(i)
	824	l41 = l41 + a41_einst(i)
	825
	826	! For ITT=13
	827	if (input_cza.ge.1 .and. itt_cza.eq.13) then
	828	a12_einst(i) = a1_020_010/3.d0 * 1.8d0/4.d0 * taustar12(i)
	829	l12=l12+a12_einst(i)
[498]	830	endif
	831
	832
[757]	833	!
[498]	834
[757]	835	a11(i) = gammanu113.d0 1.d0/2.d0 * (p11) /
	836	@ (n10(i)*l11)
	837	a1121(i,i) = (nu11/nu21)*3.d0 n20(i)/n10(i) * p1121/l11
	838	a1131(i,i) = (nu11/nu31)*3.d0 n30(i)/n10(i) * p1131/l11
	839	a1141(i,i) = (nu11/nu41)*3.d0 n40(i)/n10(i) * p1141/l11
	840	e110(i) = 2.d0* vlightnu112.d0 1.d0/2.d0 /
	841	@ ( n10(i) * l11 )
	842
	843	a21(i) = gammanu213.d0 1.d0/2.d0 *
	844	@ (p21)/(n20(i)*l21)
	845	a2111(i,i) = (nu21/nu11)*3.d0 n10(i)/n20(i) * p2111/l21
	846	a2131(i,i) = (nu21/nu31)*3.d0 n30(i)/n20(i) * p2131/l21
	847	a2141(i,i) = (nu21/nu41)*3.d0 n40(i)/n20(i) * p2141/l21
	848	e210(i) = 2.d0vlightnu21*2.d0 1.d0/2.d0 /
	849	@ ( n20(i) * l21 )
	850
	851	a31(i) = gammanu313.d0 1.d0/2.d0 * (p31) /
	852	@ (n30(i)*l31)
	853	a3111(i,i) = (nu31/nu11)*3.d0 n10(i)/n30(i) * p3111/l31
	854	a3121(i,i) = (nu31/nu21)*3.d0 n20(i)/n30(i) * p3121/l31
	855	a3141(i,i) = (nu31/nu41)*3.d0 n40(i)/n30(i) * p3141/l31
	856	e310(i) = 2.d0vlightnu31*2.d0 1.d0/2.d0 /
	857	@ ( n30(i) * l31 )
	858
	859	a41(i) = gammanu413.d0 1.d0/2.d0 * (p41) /
	860	@ (n40(i)*l41)
	861	a4111(i,i) = (nu41/nu11)*3.d0 n10(i)/n40(i) * p4111/l41
	862	a4121(i,i) = (nu41/nu21)*3.d0 n20(i)/n40(i) * p4121/l41
	863	a4131(i,i) = (nu41/nu31)*3.d0 n30(i)/n40(i) * p4131/l41
	864	e410(i) = 2.d0vlightnu41*2.d0 1.d0/2.d0 /
	865	@ ( n40(i) * l41 )
	866
	867	if (input_cza.ge.1) then
	868
	869	a1112(i,i) = (nu11/nu121)*3.d0 n11(i)/n10(i) *
	870	@ p1112/l11
	871	a2112(i,i) = (nu21/nu121)*3.d0 n11(i)/n20(i) *
	872	@ p2112/l21
	873	a3112(i,i) = (nu31/nu121)*3.d0 n11(i)/n30(i) *
	874	@ p3112/l31
	875	a4112(i,i) = (nu41/nu121)*3.d0 n11(i)/n40(i) *
	876	@ p4112/l41
	877	a12(i) = gammanu1213.d0 2.d0/4.d0* (p12)/
	878	@ (n11(i)*l12)
	879	a1211(i,i) = (nu121/nu11)*3.d0 n10(i)/n11(i) *
	880	@ p1211/l12
	881	a1221(i,i) = (nu121/nu21)*3.d0 n20(i)/n11(i) *
	882	@ p1221/l12
	883	a1231(i,i) = (nu121/nu31)*3.d0 n30(i)/n11(i) *
	884	@ p1231/l12
	885	a1241(i,i) = (nu121/nu41)*3.d0 n40(i)/n11(i) *
	886	@ p1241/l12
	887	e121(i) = 2.d0vlightnu121*2.d0 2.d0/4.d0 /
	888	@ ( n11(i) * l12 )
	889
	890	end if
	891
	892
	893	4 continue !-------------------------------------------------------
	894
	895
	896
	897	!!!!!!!!!!!! Solucion del sistema
	898
	899	!! Paso 0 : Calculo de los alphas alf11, alf21, alf31, alf41, alf12
	900
	901	call unit ( cax2, nl )
	902
[498]	903	call diago ( cax1, e110, nl )
[757]	904	call mulmmf90 ( cax3, cax1,c110, nl )
	905	call resmmf90 ( alf11, cax2,cax3, nl )
[498]	906
[757]	907	call diago ( cax1, e210, nl )
	908	call mulmmf90 ( cax3, cax1,c210, nl )
	909	call resmmf90 ( alf21, cax2,cax3, nl )
[498]	910
[757]	911	call diago ( cax1, e310, nl )
	912	call mulmmf90 ( cax3, cax1,c310, nl )
	913	call resmmf90 ( alf31, cax2,cax3, nl )
	914
	915	call diago ( cax1, e410, nl )
	916	call mulmmf90 ( cax3, cax1,c410, nl )
	917	call resmmf90 ( alf41, cax2,cax3, nl )
	918
	919	if (input_cza.ge.1) then
	920	call diago ( cax1, e121, nl )
	921	call mulmmf90 ( cax3, cax1,c121, nl )
	922	call resmmf90 ( alf12, cax2,cax3, nl )
[498]	923	endif
[757]	924
	925	!! Paso 1 : Calculo de vectores y matrices con 1 barra (aa***)
	926
[498]	927	if (input_cza.eq.0) then ! Skip paso 1, pues el12 no se calcula
	928
[757]	929	! el11
[498]	930	call sypvvv( aa11, a11,e110,sl110, nl )
	931	call samem( aa1121, a1121, nl )
	932	call samem( aa1131, a1131, nl )
	933	call samem( aa1141, a1141, nl )
	934	call samem( aalf11, alf11, nl )
[757]	935
	936	! el21
[498]	937	call sypvvv( aa21, a21,e210,sl210, nl )
	938	call samem( aa2111, a2111, nl )
	939	call samem( aa2131, a2131, nl )
	940	call samem( aa2141, a2141, nl )
	941	call samem( aalf21, alf21, nl )
	942
[757]	943	! el31
[498]	944	call sypvvv( aa31, a31,e310,sl310, nl )
	945	call samem( aa3111, a3111, nl )
	946	call samem( aa3121, a3121, nl )
	947	call samem( aa3141, a3141, nl )
	948	call samem( aalf31, alf31, nl )
	949
[757]	950	! el41
[498]	951	call sypvvv( aa41, a41,e410,sl410, nl )
	952	call samem( aa4111, a4111, nl )
	953	call samem( aa4121, a4121, nl )
	954	call samem( aa4131, a4131, nl )
	955	call samem( aalf41, alf41, nl )
	956
	957
	958	else ! (input_cza.ge.1) , FH !
	959
	960
	961	call sypvvv( v1, a12,e121,sl121, nl ) ! a12 + e121 * sl121
	962
[757]	963	! aa11
[498]	964	call sypvvv( v2, a11,e110,sl110, nl )
	965	call trucommvv( aa11 , alf12,a1112,v2, v1, nl )
[757]	966
	967	! aalf11
[498]	968	call invdiag( cax1, a1112, nl )
[757]	969	call mulmmf90( cax2, alf12, cax1, nl ) ! alf12 * (1/a1112)
	970	call mulmmf90( cax3, cax2, alf11, nl )
	971	call resmmf90( aalf11, cax3, a1211, nl )
	972	! aa1121
[498]	973	call trucodiag(aa1121, alf12,a1112,a1121, a1221, nl)
[757]	974	! aa1131
[498]	975	call trucodiag(aa1131, alf12,a1112,a1131, a1231, nl)
[757]	976	! aa1141
[498]	977	call trucodiag(aa1141, alf12,a1112,a1141, a1241, nl)
	978
[757]	979
	980	! aa21
[498]	981	call sypvvv( v2, a21,e210,sl210, nl )
	982	call trucommvv( aa21 , alf12,a2112,v2, v1, nl )
	983
[757]	984	! aalf21
[498]	985	call invdiag( cax1, a2112, nl )
[757]	986	call mulmmf90( cax2, alf12, cax1, nl ) ! alf12 * (1/a2112)
	987	call mulmmf90( cax3, cax2, alf21, nl )
	988	call resmmf90( aalf21, cax3, a1221, nl )
	989	! aa2111
[498]	990	call trucodiag(aa2111, alf12,a2112,a2111, a1211, nl)
[757]	991	! aa2131
[498]	992	call trucodiag(aa2131, alf12,a2112,a2131, a1231, nl)
[757]	993	! aa2141
[498]	994	call trucodiag(aa2141, alf12,a2112,a2141, a1241, nl)
	995
[757]	996
	997	! aa31
	998	call sypvvv ( v2, a31,e310,sl310, nl )
[498]	999	call trucommvv( aa31 , alf12,a3112,v2, v1, nl )
[757]	1000	! aalf31
[498]	1001	call invdiag( cax1, a3112, nl )
[757]	1002	call mulmmf90( cax2, alf12, cax1, nl ) ! alf12 * (1/a3112)
	1003	call mulmmf90( cax3, cax2, alf31, nl )
	1004	call resmmf90( aalf31, cax3, a1231, nl )
	1005	! aa3111
[498]	1006	call trucodiag(aa3111, alf12,a3112,a3111, a1211, nl)
[757]	1007	! aa3121
[498]	1008	call trucodiag(aa3121, alf12,a3112,a3121, a1221, nl)
[757]	1009	! aa3141
[498]	1010	call trucodiag(aa3141, alf12,a3112,a3141, a1241, nl)
	1011
[757]	1012
	1013	! aa41
[498]	1014	call sypvvv( v2, a41,e410,sl410, nl )
	1015	call trucommvv( aa41 , alf12,a4112,v2, v1, nl )
[757]	1016	! aalf41
[498]	1017	call invdiag( cax1, a4112, nl )
[757]	1018	call mulmmf90( cax2, alf12, cax1, nl ) ! alf12 * (1/a4112)
	1019	call mulmmf90( cax3, cax2, alf41, nl )
	1020	call resmmf90( aalf41, cax3, a1241, nl )
	1021	! aa4111
[498]	1022	call trucodiag(aa4111, alf12,a4112,a4111, a1211, nl)
[757]	1023	! aa4121
[498]	1024	call trucodiag(aa4121, alf12,a4112,a4121, a1221, nl)
[757]	1025	! aa4131
[498]	1026	call trucodiag(aa4131, alf12,a4112,a4131, a1231, nl)
	1027
	1028	endif ! Final caso input_cza.ge.1
	1029
	1030
[757]	1031	!! Paso 2 : Calculo de vectores y matrices con 2 barras (aaa***)
[498]	1032
[757]	1033	! aaalf41
[498]	1034	call invdiag( cax1, aa4121, nl )
[757]	1035	call mulmmf90( cax2, aalf21, cax1, nl ) ! alf21 * (1/a4121)
	1036	call mulmmf90( cax3, cax2, aalf41, nl )
	1037	call resmmf90( aaalf41, cax3, aa2141, nl )
	1038	! aaa41
[498]	1039	call trucommvv(aaa41, aalf21,aa4121,aa41, aa21, nl)
[757]	1040	! aaa4111
[498]	1041	call trucodiag(aaa4111, aalf21,aa4121,aa4111, aa2111, nl)
[757]	1042	! aaa4131
[498]	1043	call trucodiag(aaa4131, aalf21,aa4121,aa4131, aa2131, nl)
	1044
[757]	1045	! aaalf31
[498]	1046	call invdiag( cax1, aa3121, nl )
[757]	1047	call mulmmf90( cax2, aalf21, cax1, nl ) ! alf21 * (1/a3121)
	1048	call mulmmf90( cax3, cax2, aalf31, nl )
	1049	call resmmf90( aaalf31, cax3, aa2131, nl )
	1050	! aaa31
[498]	1051	call trucommvv(aaa31, aalf21,aa3121,aa31, aa21, nl)
[757]	1052	! aaa3111
[498]	1053	call trucodiag(aaa3111, aalf21,aa3121,aa3111, aa2111, nl)
[757]	1054	! aaa3141
[498]	1055	call trucodiag(aaa3141, aalf21,aa3121,aa3141, aa2141, nl)
	1056
[757]	1057	! aaalf11
[498]	1058	call invdiag( cax1, aa1121, nl )
[757]	1059	call mulmmf90( cax2, aalf21, cax1, nl ) ! alf21 * (1/a1121)
	1060	call mulmmf90( cax3, cax2, aalf11, nl )
	1061	call resmmf90( aaalf11, cax3, aa2111, nl )
	1062	! aaa11
[498]	1063	call trucommvv(aaa11, aalf21,aa1121,aa11, aa21, nl)
[757]	1064	! aaa1131
[498]	1065	call trucodiag(aaa1131, aalf21,aa1121,aa1131, aa2131, nl)
[757]	1066	! aaa1141
[498]	1067	call trucodiag(aaa1141, aalf21,aa1121,aa1141, aa2141, nl)
	1068
	1069
[757]	1070	!! Paso 3 : Calculo de vectores y matrices con 3 barras (aaaa***)
[498]	1071
[757]	1072	! aaaalf41
[498]	1073	call invdiag( cax1, aaa4131, nl )
[757]	1074	call mulmmf90( cax2, aaalf31, cax1, nl ) ! aaalf31 * (1/aaa4131)
	1075	call mulmmf90( cax3, cax2, aaalf41, nl )
	1076	call resmmf90( aaaalf41, cax3, aaa3141, nl )
	1077	! aaaa41
[498]	1078	call trucommvv(aaaa41, aaalf31,aaa4131,aaa41, aaa31, nl)
[757]	1079	! aaaa4111
[498]	1080	call trucodiag(aaaa4111, aaalf31,aaa4131,aaa4111,aaa3111, nl)
	1081
[757]	1082	! aaaalf11
[498]	1083	call invdiag( cax1, aaa1131, nl )
[757]	1084	call mulmmf90( cax2, aaalf31, cax1, nl ) ! aaalf31 * (1/aaa4131)
	1085	call mulmmf90( cax3, cax2, aaalf11, nl )
	1086	call resmmf90( aaaalf11, cax3, aaa3111, nl )
	1087	! aaaa11
[498]	1088	call trucommvv(aaaa11, aaalf31,aaa1131,aaa11, aaa31, nl)
[757]	1089	! aaaa1141
[498]	1090	call trucodiag(aaaa1141, aaalf31,aaa1131,aaa1141,aaa3141, nl)
	1091
	1092
[757]	1093	!! Paso 4 : Calculo de vectores y matrices finales y calculo de J1
[498]	1094
	1095	call trucommvv(v1, aaaalf41,aaaa1141,aaaa11, aaaa41, nl)
[757]	1096	!
[498]	1097	call invdiag( cax1, aaaa1141, nl )
[757]	1098	call mulmmf90( cax2, aaaalf41, cax1, nl ) ! aaaalf41 * (1/aaaa1141)
	1099	call mulmmf90( cax3, cax2, aaaalf11, nl )
	1100	call resmmf90( cax1, cax3, aaaa4111, nl )
	1101	!
[498]	1102	call LUdec ( el11, cax1, v1, nl, nl2 )
	1103
[757]	1104	! Solucion para el41
[498]	1105	call sypvmv( v1, aaaa41, aaaa4111,el11, nl )
	1106	call LUdec ( el41, aaaalf41, v1, nl, nl2 )
	1107
[757]	1108	! Solucion para el31
[498]	1109	call sypvmv( v2, aaa31, aaa3111,el11, nl )
	1110	call sypvmv( v1, v2, aaa3141,el41, nl )
	1111	call LUdec ( el31, aaalf31, v1, nl, nl2 )
	1112
[757]	1113	! Solucion para el21
[498]	1114	call sypvmv( v3, aa21, aa2111,el11, nl )
	1115	call sypvmv( v2, v3, aa2131,el31, nl )
	1116	call sypvmv( v1, v2, aa2141,el41, nl )
	1117	call LUdec ( el21, aalf21, v1, nl, nl2 )
	1118
[757]	1119	!!!
	1120	el11(1) = planckdp( t(1), nu11 )
	1121	el21(1) = planckdp( t(1), nu21 )
	1122	el31(1) = planckdp( t(1), nu31 )
	1123	el41(1) = planckdp( t(1), nu41 )
	1124	el11(nl) = 2.d0 * el11(nl-1) - el11(nl2)
	1125	el21(nl) = 2.d0 * el21(nl-1) - el21(nl2)
	1126	el31(nl) = 2.d0 * el31(nl-1) - el31(nl2)
	1127	el41(nl) = 2.d0 * el41(nl-1) - el41(nl2)
[498]	1128
[757]	1129	call mulmv ( v1, c110,el11, nl )
	1130	call sumvv ( hr110, v1,sl110, nl )
[498]	1131
[757]	1132	! Solucion para el12
	1133	if (input_cza.ge.1) then
	1134
[498]	1135	call sypvmv( v1, a12, a1211,el11, nl )
	1136	call sypvmv( v3, v1, a1221,el21, nl )
	1137	call sypvmv( v2, v3, a1231,el31, nl )
	1138	call sypvmv( v1, v2, a1241,el41, nl )
	1139	call LUdec ( el12, alf12, v1, nl, nl2 )
	1140
[757]	1141	el12(1) = planckdp( t(1), nu121 )
	1142	el12(nl) = 2.d0 * el12(nl-1) - el12(nl2)
[498]	1143
[757]	1144	if (itt_cza.eq.15) then
	1145	call mulmv ( v1, c121,el12, nl )
	1146	call sumvv ( hr121, v1,sl121, nl )
[498]	1147	endif
	1148
[757]	1149	end if
[498]	1150
	1151
[757]	1152
	1153	if (input_cza.lt.1) then
	1154
	1155	do i=1,nl
	1156	pl11 = el11(i)/( gamma * nu11*3.0d0 1.d0/2.d0 /n10(i) )
	1157	pl21 = el21(i)/( gamma * nu21*3.0d0 1.d0/2.d0 /n20(i) )
	1158	pl31 = el31(i)/( gamma * nu31*3.0d0 1.d0/2.d0 /n30(i) )
	1159	pl41 = el41(i)/( gamma * nu41*3.0d0 1.d0/2.d0 /n40(i) )
	1160	vt11(i) = -eenu11 / log( abs(pl11) / (2.0d0n10(i)) )
	1161	vt21(i) = -eenu21 / log( abs(pl21) / (2.0d0n20(i)) )
	1162	vt31(i) = -eenu31 / log( abs(pl31) / (2.0d0n30(i)) )
	1163	vt41(i) = -eenu41 / log( abs(pl41) / (2.0d0n40(i)) )
	1164	hr210(i) = sl210(i) -hplanckvlightnu21 a21_einst(i)pl21
	1165	hr310(i) = sl310(i) -hplanckvlightnu31 a31_einst(i)pl31
	1166	hr410(i) = sl410(i) -hplanckvlightnu41 a41_einst(i)pl41
	1167	enddo
	1168
	1169	v626t1(1:nl)=vt11(1:nl)
	1170	v628t1(1:nl)=vt21(1:nl)
	1171	v636t1(1:nl)=vt31(1:nl)
	1172	v627t1(1:nl)=vt41(1:nl)
	1173	! call dinterconnection( v626t1, vt11 )
	1174	! call dinterconnection ( v628t1, vt21 )
	1175	! call dinterconnection ( v636t1, vt31 )
	1176	! call dinterconnection ( v627t1, vt41 )
	1177
[498]	1178	else
[757]	1179
	1180	do i=1,nl
	1181	pl21 = el21(i)/( gamma * nu21*3.0d0 1.d0/2.d0 / n20(i) )
	1182	pl31 = el31(i)/( gamma * nu31*3.0d0 1.d0/2.d0 / n30(i) )
	1183	pl41 = el41(i)/( gamma * nu41*3.0d0 1.d0/2.d0 / n40(i) )
	1184	hr210(i) = sl210(i) -hplanckvlightnu21 a21_einst(i)pl21
	1185	hr310(i) = sl310(i) -hplanckvlightnu31 a31_einst(i)pl31
	1186	hr410(i) = sl410(i) -hplanckvlightnu41 a41_einst(i)pl41
	1187	if (itt_cza.eq.13) then
	1188	pl12 = el12(i)/( gammanu1213.0d0 2.d0/4.d0 /n11(i) )
	1189	hr121(i) = - hplanckvlight nu121 * a12_einst(i)*pl12
	1190	hr121(i) = hr121(i) + sl121(i)
	1191	endif
[498]	1192	enddo
	1193
	1194	endif
	1195
[757]	1196	! K/Dday
	1197	do i=1,nl
	1198	hr110(i)=hr110(i)*dble( hrkday_factor(i) / nt(i) )
	1199	hr210(i)=hr210(i)*dble( hrkday_factor(i) / nt(i) )
	1200	hr310(i)=hr310(i)*dble( hrkday_factor(i) / nt(i) )
	1201	hr410(i)=hr410(i)*dble( hrkday_factor(i) / nt(i) )
	1202	hr121(i)=hr121(i)*dble( hrkday_factor(i) / nt(i) )
	1203	end do
[498]	1204
	1205
[757]	1206	c final
	1207	return
	1208	c
	1209	end
[498]	1210
	1211
[757]	1212	c * Old NLTEdlvr11_FB626CTS_02 *
[498]	1213
	1214	c***********************************************************************
[757]	1215
	1216	subroutine NLTEdlvr11_FB626CTS ( hr110CTS, nl_cts_real )
[498]	1217
	1218	c***********************************************************************
	1219
[757]	1220	implicit none
[498]	1221
[757]	1222	!!!!!!!!!!!!!!!!!! common variables and constants
[498]	1223
[757]	1224	include 'nlte_paramdef.h'
	1225	include 'nlte_commons.h'
	1226
	1227
	1228	c Arguments
	1229	real*8 hr110CTS(nl_cts) ! output
	1230	integer nl_cts_real ! i
	1231
	1232	c local variables
	1233
	1234	real*8 n11CTS(nl_cts), slopeTstar110(nl_cts)
	1235	real*8 n10(nl_cts), co2t, codbl, n2dbl, o3pdbl
	1236	real*8 d19c1, d19cp1, l11, p11
	1237	real*8 a11_einst(nl_cts), hcv, maxslope
	1238	integer i, isot
	1239
	1240	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! start program
	1241
	1242	nu11 = dble(nu(1,1))
	1243	hcv = hplanckvlightnu11
	1244
	1245	call zero2v (hr110CTS,n11CTS,nl_cts)
	1246
	1247	do i=1,nl_cts_real
	1248
	1249	co2t = dble ( co2_cts(i) *(imr(1)+imr(3)+imr(2)+imr(4)) )
	1250	n10(i) = dble( co2_cts(i) * imr(1) )
	1251	codbl = dble(co_cts(i))
	1252	o3pdbl = dble(o3p_cts(i))
	1253	n2dbl = dble(n2_cts(i))
	1254
	1255	call GETK_dlvr11 ( t_cts(i) )
	1256	isot = 1
	1257	d19c1 = k19ca(isot)co2t + k19cb(isot)n2dbl
	1258	$ + k19cc(isot)*codbl
	1259	d19cp1 = k19cap(isot)co2t + k19cbp(isot)n2dbl
	1260	$ + k19ccp(isot)*codbl
	1261	l11 = d19c1 + k20c(1)*o3pdbl
	1262	p11 = ( d19cp1 + k20cp(1)o3pdbl ) n10(i)
	1263
	1264	a11_einst(i) = a1_010_000 * 1.8d0/4.d0 * taustar11_cts(i)
	1265
	1266	n11CTS(i) = p11 / (l11 + a11_einst(i))
	1267
	1268	hr110CTS(i) = - n11CTS(i) * a11_einst(i) * hcv
	1269	hr110CTS(i) = hr110CTS(i)*
	1270	$ dble( hrkday_factor_cts(i) / nt_cts(i) ) !K/Day
	1271
	1272	enddo
	1273
	1274
	1275	c calculo de la altura de transicion, a partir de Tstar
	1276	c y merging con el hr110(i), ya calculado con CZALU
	1277
	1278	slopeTstar110(1) = taustar11_cts(2)-taustar11_cts(1)
	1279	slopeTstar110(nl_cts_real) = taustar11_cts(nl_cts_real) -
	1280	$ taustar11_cts(nl_cts_real-1)
	1281	maxslope = max( slopeTstar110(1),slopeTstar110(nl_cts_real))
	1282	if (nl_cts_real .gt. 2) then
	1283	do i=2,nl_cts_real-1
	1284	slopeTstar110(i) = ( taustar11_cts(i+1) -
	1285	$ taustar11_cts(i-1) ) * 0.5d0
	1286	if ( slopeTstar110(i) .gt. maxslope ) then
	1287	!write (,) i, pl_cts(i), maxslope, slopeTstar110(i)
	1288	maxslope=slopeTstar110(i)
	1289	endif
	1290	enddo
	1291	endif
	1292
	1293	c
	1294	return
	1295	end
	1296
	1297
[498]	1298	c***********************************************************************
[757]	1299	c hrkday_convert.f
	1300	c
	1301	c fortran function that returns the factor for conversion from
	1302	c hr' [erg s-1 cm-3] to hr [ k day-1 ]
	1303	c
	1304	c mar 2010 fgg adapted to GCM
	1305	c jan 99 malv add o2 as major component.
	1306	c ago 98 malv also returns cp_avg,pm_avg
	1307	c jul 98 malv first version.
[498]	1308	c***********************************************************************
[757]	1309
	1310	function hrkday_convert
	1311	@ ( mmean_nlte,cpmean_nlte )
	1312
	1313	implicit none
	1314
	1315	include 'comcstfi.h'
	1316	include 'param.h'
	1317
	1318	c argumentos
	1319	real mmean_nlte,cpmean_nlte
	1320	real hrkday_convert
	1321
	1322	ccccccccccccccccccccccccccccccccccccc
	1323
	1324	hrkday_convert = daysec * n_avog /
	1325	& ( cpmean_nlte * 1.e4 * mmean_nlte )
	1326
	1327	c end
	1328	return
	1329	end

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

Download in other formats:

Original Format