Context Navigation

nlte_tcool.F @ 3017

Last change on this file since 3017 was 3012, checked in by emillour, 3 years ago
Mars PCM: Code cleanup concerning chemistry. Turn chimiedata.h into module chemistrydata.F90 and integrate read_phototable.F in it. Also fix an OpenMP issue with read_phototable; different OpenMP threads should not simultaneously open/read a file. Let the master do the job and then broadcast the result to all cores. While at it, also turned nltecool.F and nlte_tcool.F into modules. EM
File size: 56.2 KB

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

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

Download in other formats:

Original Format