Context Navigation

nlte_tcool.F @ 1543

Last change on this file since 1543 was 1530, checked in by emillour, 9 years ago

Venus and Titan GCMs:
Updates in the physics to keep up with updates in LMDZ5 (up to
LMDZ5 trunk, rev 2350) concerning dynamics/physics separation:

Adapted makelmdz and makelmdz_fcm script to stop if trying to compile 1d model or newstart or start2archive in parallel.
got rid of references to "dimensions.h" in physics. Within physics packages, use nbp_lon (=iim), nbp_lat (=jjmp1) and nbp_lev (=llm) from module mod_grid_phy_lmdz (in phy_common) instead. Only partially done for Titan, because of many hard-coded commons; a necessary first step will be to clean these up (using modules).

File size: 56.4 KB

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

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

Download in other formats:

Original Format