Context Navigation

nlte_tcool.F @ 1009

Last change on this file since 1009 was 759, checked in by emillour, 12 years ago

Mars GCM:

commented out some messages in nlte_tcool.F (which were informations

useful and meaningful during code development only)

FGG

File size: 44.6 KB

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

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

Download in other formats:

Original Format