Context Navigation

nlte_tcool.F @ 3026

Last change on this file since 3026 was 3018, checked in by emillour, 17 months ago
Mars PCM: Further code cleanup with NLTE routines; converted nlte_paramdef.h to module nlte_paramdef_h.F90 and nlte_commons.h to module nlte_commons_h.F90 (could not turn nlte_aux.F, nlte_setup.F and nlte_calc.F into modules due to circular dependencies; would require further code reorganization). EM
File size: 58.8 KB

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

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

Download in other formats:

Original Format