Context Navigation

nlte_tcool.F @ 3461

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

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

Download in other formats:

Original Format