Context Navigation

nlte_tcool.F @ 3325

Last change on this file since 3325 was 3018, checked in by emillour, 16 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

Line
1	MODULE nlte_tcool_mod
2
3	IMPLICIT NONE
4
5	CONTAINS
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	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
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
20	c -NLTEdlvr11_ERRORS
21	c
22	c
23	c
24	c * Old NLTEdlvr11_TCOOL_02 *
25	c
26	c***********************************************************************
27
28	c***********************************************************************
29
30	subroutine nlte_tcool(ngridgcm,n_gcm,
31	$ p_gcm, t_gcm, z_gcm,
32	$ co2vmr_gcm, n2vmr_gcm, covmr_gcm, o3pvmr_gcm,
33	$ q15umco2_gcm , ierr, varerr)
34
35	c***********************************************************************
36
37	use conc_mod, only: cpnew, mmean
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
44	implicit none
45
46	c Arguments
47	integer n_gcm,ngridgcm
48	real p_gcm(ngridgcm,n_gcm), t_gcm(ngridgcm,n_gcm)
49	real z_gcm(ngridgcm,n_gcm)
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)
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
58
59
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
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
73	c***************
74	c***************
75
76	do ig=1,ngridgcm
77	ierr = 0
78	nl_cts_real = 0
79	nzy_cts_real = 0
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
92	! From GCM's grid to NLTE's grid
93	call NLTEdlvr11_ZGRID (n_gcm,
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 )
98
99
100	! Isotopic Tstar & VC at the NLTE grid
101	call interdp_ESCTVCISO
102
103	! Tstar para NLTE-CTS
104	call MZESC110 ( ig,nl_cts_real, nzy_cts_real,ierr,varerr )
105	if (ierr .gt. 0) call ERRORS (ierr,varerr)
106
107	! 626FB C.M.
108	call leetvt
109	c110(1:nl,1:nl)=0.d0
110	! call zerom (c110, nl)
111	call zero2v (vc110,taustar11, nl)
112	call MZTUD110 ( ierr, varerr )
113	if (ierr .gt. 0) call ERRORS (ierr,varerr)
114
115	input_cza = 0
116	call NLTEdlvr11_CZALU(ierr,varerr)
117	if (ierr .gt. 0) call ERRORS (ierr,varerr)
118
119	input_cza = 1
120	call NLTEdlvr11_CZALU(ierr,varerr)
121	if (ierr .gt. 0) call ERRORS (ierr,varerr)
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)
144	@ + hr121(i)
145	enddo
146
147
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
153	q15umco2_nltot(i) = 0.5d0*( q15umco2_nltot(i) + hr110CTS(1) )
154	do i=2,nl_cts_real
155	indice = (nl-2) + i
156	q15umco2_nltot(indice) = hr110CTS(i)
157	enddo
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
201	enddo
202	c end subroutine
203
204	end subroutine nlte_tcool
205
206
207	c***********************************************************************
208
209	subroutine NLTEdlvr11_ZGRID (n_gcm,
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	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
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
240	real zaux_gcm(n_gcm)
241
242	c local variables
243	integer i, iz
244	real distancia, meanm, gz, Hkm
245	real zmin, zmax
246	real mmean_nlte(n_gcm),cpnew_nlte(n_gcm)
247	real gg,masa,radio,kboltzman
248
249	c functions
250	! external hrkday_convert
251	! real hrkday_convert
252
253	c***********************************************************************
254
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
259
260	! Primero, construimos escala z_gcm
261
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
267
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
278
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
282
283	!!!!!!!!!Primero el NLTE completo !!!!!!!!
284
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
288	do while ( p_gcm(jlowerboundary) .gt. Pbottom_atm )
289	jlowerboundary = jlowerboundary + 1
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
295	call abort_physic("nlte_tcool",
296	& 'Check input value "p_gcm" or modify "Pbottom_atm"',1)
297	endif
298	enddo
299
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 )
304	jtopboundary = jtopboundary + 1
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
315	enddo
316	5000 continue
317
318	! Grid steps
319
320	zmin = z_gcm(jlowerboundary)
321	zmax = z_gcm(jtopboundary)
322	deltaz = (zmax-zmin) / (nl-1)
323	do i=1,nl
324	zl(i) = zmin + (i-1) * deltaz
325	enddo
326
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
340	do i = 1, nl
341	nt(i) = 7.339e+21 * pl(i) / t(i) ! --> [cm-3]
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)
346	! hrkday_factor(i) = hrkday_convert( t(i),
347	! $ co2vmr(i), o3pvmr(i), n2vmr(i), covmr(i) )
348	hrkday_factor(i) = hrkday_convert(mmean_nlte(i)
349	& ,cpnew_nlte(i))
350	enddo
351
352	! Comprobar que las temps no se salen del grid del histograma
353
354	do i=1,nl
355	if (t(i) .gt. 500.0) then
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
365	enddo
366
367	! Fine grid for transmittance calculations
368
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)
378
379	do i=1,nzy
380	nty(i) = 7.339e+21 * py(i) / ty(i) ! --> [cm-3]
381	co2y(i) = co2y(i) * nty(i)
382	enddo
383
384
385	!!!!!!!!!Segundo, el NLTE - CTS !!!!!!!!
386
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
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'
400	nl_cts_real=i-1
401	! write (,) ' Original,Real NL_CTS=', nl_cts,nl_cts_real
402	goto 6000
403	endif
404	enddo
405	nl_cts_real = nl_cts
406	6000 continue
407
408	! Creamos perfil por interpolacion
409
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)
423
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
436	if (t_cts(i) .gt. thist_stored(1,mm_stored(1))) then
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
491	end subroutine NLTEdlvr11_ZGRID
492
493
494	c * Old NLTEdlvr11_CZALU_03 *
495
496	c**********************************************************************
497
498
499	subroutine NLTEdlvr11_CZALU(ierr,varerr)
500
501	c***********************************************************************
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
524	implicit none
525
526	c Arguments
527
528	integer ierr
529	real*8 varerr
530
531
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)
538	real*8 f1(nl,nl)
539
540	real*8 cax1(nl,nl), cax2(nl,nl), cax3(nl,nl)
541	real*8 v1(nl), v2(nl), v3(nl)
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)
554
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)
578	real*8 n20(nl), n21(nl)
579	real*8 n30(nl), n31(nl)
580	real*8 n40(nl), n41(nl)
581
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
591
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
599
600	real*8 pl11, pl12, pl21, pl31, pl41
601
602	real*8 minvt11, minvt21, minvt31, minvt41
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
621	ierr = 0
622	varerr = 0.d0
623
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
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)
645
646	call zero4v (el11,el21,el31,el41,nl)
647	call zero4v (e110,e210,e310,e410,nl)
648	call zero2v (el12,e121,nl)
649
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
718	call MZTVC121( ierr,varerr )
719	if (ierr .gt. 0) call ERRORS (ierr,varerr)
720
721	endif
722
723	endif
724
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
738
739
740
741	do 4,i=nl,1,-1 !----------------------------------------------
742
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
752	isot = 1
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
761	isot = 2
762	d19c2 = k19ca(isot)co2t + k19cb(isot)n2dble
763	@ + k19cc(isot)*codble
764	d19cp2 = k19cap(isot)co2t + k19cbp(isot)n2dble
765	@ + k19ccp(isot)*codble
766	isot = 3
767	d19c3 = k19ca(isot)co2t + k19cb(isot)n2dble
768	@ + k19cc(isot)*codble
769	d19cp3 = k19cap(isot)co2t + k19cbp(isot)n2dble
770	@ + k19ccp(isot)*codble
771	isot = 4
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)
825	@ + k33cp(2)n21(i) + k33cp(3)n31(i) + k33cp(4)*n41(i)
826	p1112 = d19b1
827	@ + k20b(1)*o3pdbl
828	@ + 2.d0k21b(1)n10(i)
829	@ + k33c*( n20(i) + n30(i) + n40(i) )
830
831	l21 = l21 + k33cp(2)*n11(i)
832	p2112 = k33c*n20(i)
833
834	l31 = l31 + k33cp(3)*n11(i)
835	p3112 = k33c*n30(i)
836
837	l41 = l41 + k33cp(4)*n11(i)
838	p4112 = k33c*n40(i)
839
840	end if
841
842
843	! For ITT=13,15
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
852	! For ITT=13
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)
856	endif
857
858
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	!
885
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
954	call diago ( cax1, e110, nl )
955	call mulmmf90 ( cax3, cax1,c110, nl )
956	call resmmf90 ( alf11, cax2,cax3, nl )
957
958	call diago ( cax1, e210, nl )
959	call mulmmf90 ( cax3, cax1,c210, nl )
960	call resmmf90 ( alf21, cax2,cax3, nl )
961
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 )
974	endif
975
976	!! Paso 1 : Calculo de vectores y matrices con 1 barra (aa***)
977
978	if (input_cza.eq.0) then ! Skip paso 1, pues el12 no se calcula
979
980	! el11
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 )
986
987	! el21
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
994	! el31
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
1001	! el41
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
1014	! aa11
1015	call sypvvv( v2, a11,e110,sl110, nl )
1016	call trucommvv( aa11 , alf12,a1112,v2, v1, nl )
1017
1018	! aalf11
1019	call invdiag( cax1, a1112, nl )
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
1024	call trucodiag(aa1121, alf12,a1112,a1121, a1221, nl)
1025	! aa1131
1026	call trucodiag(aa1131, alf12,a1112,a1131, a1231, nl)
1027	! aa1141
1028	call trucodiag(aa1141, alf12,a1112,a1141, a1241, nl)
1029
1030
1031	! aa21
1032	call sypvvv( v2, a21,e210,sl210, nl )
1033	call trucommvv( aa21 , alf12,a2112,v2, v1, nl )
1034
1035	! aalf21
1036	call invdiag( cax1, a2112, nl )
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
1041	call trucodiag(aa2111, alf12,a2112,a2111, a1211, nl)
1042	! aa2131
1043	call trucodiag(aa2131, alf12,a2112,a2131, a1231, nl)
1044	! aa2141
1045	call trucodiag(aa2141, alf12,a2112,a2141, a1241, nl)
1046
1047
1048	! aa31
1049	call sypvvv ( v2, a31,e310,sl310, nl )
1050	call trucommvv( aa31 , alf12,a3112,v2, v1, nl )
1051	! aalf31
1052	call invdiag( cax1, a3112, nl )
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
1057	call trucodiag(aa3111, alf12,a3112,a3111, a1211, nl)
1058	! aa3121
1059	call trucodiag(aa3121, alf12,a3112,a3121, a1221, nl)
1060	! aa3141
1061	call trucodiag(aa3141, alf12,a3112,a3141, a1241, nl)
1062
1063
1064	! aa41
1065	call sypvvv( v2, a41,e410,sl410, nl )
1066	call trucommvv( aa41 , alf12,a4112,v2, v1, nl )
1067	! aalf41
1068	call invdiag( cax1, a4112, nl )
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
1073	call trucodiag(aa4111, alf12,a4112,a4111, a1211, nl)
1074	! aa4121
1075	call trucodiag(aa4121, alf12,a4112,a4121, a1221, nl)
1076	! aa4131
1077	call trucodiag(aa4131, alf12,a4112,a4131, a1231, nl)
1078
1079	endif ! Final caso input_cza.ge.1
1080
1081
1082	!! Paso 2 : Calculo de vectores y matrices con 2 barras (aaa***)
1083
1084	! aaalf41
1085	call invdiag( cax1, aa4121, nl )
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
1090	call trucommvv(aaa41, aalf21,aa4121,aa41, aa21, nl)
1091	! aaa4111
1092	call trucodiag(aaa4111, aalf21,aa4121,aa4111, aa2111, nl)
1093	! aaa4131
1094	call trucodiag(aaa4131, aalf21,aa4121,aa4131, aa2131, nl)
1095
1096	! aaalf31
1097	call invdiag( cax1, aa3121, nl )
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
1102	call trucommvv(aaa31, aalf21,aa3121,aa31, aa21, nl)
1103	! aaa3111
1104	call trucodiag(aaa3111, aalf21,aa3121,aa3111, aa2111, nl)
1105	! aaa3141
1106	call trucodiag(aaa3141, aalf21,aa3121,aa3141, aa2141, nl)
1107
1108	! aaalf11
1109	call invdiag( cax1, aa1121, nl )
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
1114	call trucommvv(aaa11, aalf21,aa1121,aa11, aa21, nl)
1115	! aaa1131
1116	call trucodiag(aaa1131, aalf21,aa1121,aa1131, aa2131, nl)
1117	! aaa1141
1118	call trucodiag(aaa1141, aalf21,aa1121,aa1141, aa2141, nl)
1119
1120
1121	!! Paso 3 : Calculo de vectores y matrices con 3 barras (aaaa***)
1122
1123	! aaaalf41
1124	call invdiag( cax1, aaa4131, nl )
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
1129	call trucommvv(aaaa41, aaalf31,aaa4131,aaa41, aaa31, nl)
1130	! aaaa4111
1131	call trucodiag(aaaa4111, aaalf31,aaa4131,aaa4111,aaa3111, nl)
1132
1133	! aaaalf11
1134	call invdiag( cax1, aaa1131, nl )
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
1139	call trucommvv(aaaa11, aaalf31,aaa1131,aaa11, aaa31, nl)
1140	! aaaa1141
1141	call trucodiag(aaaa1141, aaalf31,aaa1131,aaa1141,aaa3141, nl)
1142
1143
1144	!! Paso 4 : Calculo de vectores y matrices finales y calculo de J1
1145
1146	call trucommvv(v1, aaaalf41,aaaa1141,aaaa11, aaaa41, nl)
1147	!
1148	call invdiag( cax1, aaaa1141, nl )
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	!
1153	call LUdec ( el11, cax1, v1, nl, nl2 )
1154
1155	! Solucion para el41
1156	call sypvmv( v1, aaaa41, aaaa4111,el11, nl )
1157	call LUdec ( el41, aaaalf41, v1, nl, nl2 )
1158
1159	! Solucion para el31
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
1164	! Solucion para el21
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
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)
1179
1180	call mulmv ( v1, c110,el11, nl )
1181	call sumvv ( hr110, v1,sl110, nl )
1182
1183	! Solucion para el12
1184	if (input_cza.ge.1) then
1185
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
1192	el12(1) = planckdp( t(1), nu121 )
1193	el12(nl) = 2.d0 * el12(nl-1) - el12(nl2)
1194
1195	if (itt_cza.eq.15) then
1196	call mulmv ( v1, c121,el12, nl )
1197	call sumvv ( hr121, v1,sl121, nl )
1198	endif
1199
1200	end if
1201
1202
1203
1204	if (input_cza.lt.1) then
1205
1206	minvt11 = 1.d6
1207	minvt21 = 1.d6
1208	minvt31 = 1.d6
1209	minvt41 = 1.d6
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
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) )
1227	enddo
1228
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
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
1257	else
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
1271	enddo
1272
1273	endif
1274
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
1283
1284
1285	c final
1286
1287	c
1288	end subroutine NLTEdlvr11_CZALU
1289
1290
1291	c * Old NLTEdlvr11_FB626CTS_02 *
1292
1293	c***********************************************************************
1294
1295	subroutine NLTEdlvr11_FB626CTS ( hr110CTS, nl_cts_real )
1296
1297	c***********************************************************************
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
1306	implicit none
1307
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
1374	end subroutine NLTEdlvr11_FB626CTS
1375
1376
1377	c***********************************************************************
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.
1387	c***********************************************************************
1388
1389	function hrkday_convert
1390	@ ( mmean_nlte,cpmean_nlte )
1391	use time_phylmdz_mod, only: daysec
1392	use param_v4_h, only: n_avog
1393	implicit none
1394
1395	c argumentos
1396	real mmean_nlte,cpmean_nlte
1397	real hrkday_convert
1398
1399	ccccccccccccccccccccccccccccccccccccc
1400
1401	hrkday_convert = daysec * n_avog /
1402	& ( cpmean_nlte * 1.e4 * mmean_nlte )
1403
1404	c end
1405
1406	end function hrkday_convert
1407
1408
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
1666	call abort_physic("nlte_tcool",
1667	& 'Stopped in NLTE scheme due to severe error.',1)
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