Context Navigation

nlte_tcool.F @ 1453

Last change on this file since 1453 was 1442, checked in by slebonnois, 9 years ago
SL: update of the Venus GCM, + corrections on routines used for newstart/start2archive for Titan and Venus, + some modifications on tools
File size: 56.4 KB

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

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

Download in other formats:

Original Format