Context Navigation

nlte_tcool.F @ 1396

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

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

Download in other formats:

Original Format