Context Navigation

nlte_tcool.F @ 1242

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

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

Download in other formats:

Original Format