Context Navigation

nlte_tcool.F @ 2932

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

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

Download in other formats:

Original Format