Context Navigation

nlte_tcool.F @ 1233

Last change on this file since 1233 was 1226, checked in by aslmd, 11 years ago
LMDZ.MARS : Replaced comcstfi and planete includes by modules.
File size: 56.0 KB

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

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

Download in other formats:

Original Format