Context Navigation

nlte_tcool.F @ 1720

Last change on this file since 1720 was 1524, checked in by emillour, 9 years ago

All GCMS:
More updates to enforce dynamics/physics separation:

get rid of references to "temps_mod" from physics packages;
make a "time_phylmdz_mod.F90" module to store that
information and fill it via "iniphysiq".

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

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

Download in other formats:

Original Format