Context Navigation

nlte_tcool.F @ 2325

Last change on this file since 2325 was 1524, checked in by emillour, 9 years ago

All GCMS:
More updates to enforce dynamics/physics separation:

get rid of references to "temps_mod" from physics packages;
make a "time_phylmdz_mod.F90" module to store that
information and fill it via "iniphysiq".

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

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

Download in other formats:

Original Format