Context Navigation

nlte_tcool.F @ 2800

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

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

Download in other formats:

Original Format