Context Navigation

nlte_tcool.F @ 3884

Last change on this file since 3884 was 3884, checked in by ikovalenko, 4 months ago

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

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

Download in other formats:

Original Format