Context Navigation

nlte_tcool.F @ 3556

Last change on this file since 3556 was 2048, checked in by slebonnois, 6 years ago
SL: VENUS, autres details
File size: 56.5 KB

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

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

Download in other formats:

Original Format