Context Navigation

NLTEdlvr09_CZALU.F @ 496

Last change on this file since 496 was 414, checked in by aslmd, 14 years ago

LMDZ.MARS : NEW NLTE MODEL FROM GRANADA AMIGOS

23/11/11 == FGG + MALV

New parameterization of the NLTE 15 micron cooling. The old parameterization is kept as an option, including or not variable atomic oxygen concentration. A new flag is introduced in callphys.def, nltemodel, to select which parameterization wants to be used (new one, old one with variable [O], or old one with fixed [O], see below). Includes many new subroutines and commons in phymars. Some existing routines are also modified:

-physiq.F. Call to the new subroutine NLTE_leedat in first call. Call to nltecool modified to allow for variable atomic oxygen. Depending on the value of nltemodel, the new subroutine NLTEdlvr09_TCOOL is called instead of nltecool.

-inifis.F. Reading of nltemodel is added.

-callkeys.h Declaration of nltemodel is added.

The following lines should be added to callphys.def (ideally after setting callnlte):

# NLTE 15um scheme to use.
# 0-> Old scheme, static oxygen
# 1-> Old scheme, dynamic oxygen
# 2-> New scheme
nltemodel = 2

A new directory, NLTEDAT, has to be included in datagcm.

Improvements into NLTE NIR heating parameterization to take into account variability with O/CO2 ratio and SZA. A new subroutine, NIR_leedat.F, and a new common, NIRdata.h, are included. A new flag, nircorr, is added in callphys.def, to include or not these corrections. The following files are modified:

-nirco2abs.F: nq and pq are added in the arguments. The corrections factors are interpolated to the GCM grid and included in the clculation. A new subroutine, interpnir, is added at the end of the file.

-physiq.F: Call to NIR_leedat added at first call. Modified call to nirco2abs

-inifis: Reading new flag nircorr.

-callkeys.h: Declaration of nircorr.

The following lines have to be added to callphys.def (ideally after callnirco2):

# NIR NLTE correction for variable SZA and O/CO2?
# matters only if callnirco2=T
# 0-> no correction
# 1-> include correction
nircorr=1

A new data file, NIRcorrection_feb2011.dat, has to be included in datagcm.

Small changes to the molecular diffusion scheme to fix the number of species considered, to avoid problems when compiling with more than 15 tracers (for example, when CH4 is included). Modified subroutines: aeronomars/moldiff.F and aeronomars/moldiffcoeff.F

File size: 36.1 KB

Line
1	c***********************************************************************
2
3	subroutine NLTEdlvr09_CZALU(ig)
4
5	c jul 2011 malv+fgg
6	c***********************************************************************
7
8	implicit none
9
10	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! common variables and constants
11
12	include 'nltedefs.h'
13	include 'nlte_atm.h'
14	include 'nlte_data.h'
15	include 'nlte_results.h'
16	include 'tcr_15um.h'
17	include 'nlte_matrix.h'
18	include 'nlte_rates.h'
19	include 'nlte_curtis.h'
20
21	c arguments
22
23	integer ig !ADDED FOR TRACEBACK
24
25	c local variables
26
27	! matrixes and vectors
28
29	real*8 e110(nl), e210(nl), e310(nl), e410(nl)
30	real*8 e121(nl), e112(nl)
31
32	real*8 f1(nl,nl)
33
34	real*8 cax1(nl,nl), cax2(nl,nl), cax3(nl,nl)
35	real*8 v1(nl), v2(nl), v3(nl)
36
37	real*8 alf11(nl,nl), alf12(nl,nl)
38	real*8 alf21(nl,nl), alf31(nl,nl), alf41(nl,nl)
39	real*8 a11(nl), a1112(nl,nl)
40	real*8 a1121(nl,nl), a1131(nl,nl), a1141(nl,nl)
41	real*8 a21(nl), a2131(nl,nl), a2141(nl,nl)
42	real*8 a2111(nl,nl), a2112(nl,nl)
43	real*8 a31(nl), a3121(nl,nl), a3141(nl,nl)
44	real*8 a3111(nl,nl), a3112(nl,nl)
45	real*8 a41(nl), a4121(nl,nl), a4131(nl,nl)
46	real*8 a4111(nl,nl), a4112(nl,nl)
47	real*8 a12(nl), a1211(nl,nl)
48	real*8 a1221(nl,nl), a1231(nl,nl), a1241(nl,nl)
49
50	real*8 aalf11(nl,nl),aalf21(nl,nl),aalf31(nl,nl),aalf41(nl,nl)
51	real*8 aa11(nl), aa1121(nl,nl), aa1131(nl,nl), aa1141(nl,nl)
52	real*8 aa21(nl), aa2111(nl,nl), aa2131(nl,nl), aa2141(nl,nl)
53	real*8 aa31(nl), aa3111(nl,nl), aa3121(nl,nl), aa3141(nl,nl)
54	real*8 aa41(nl), aa4111(nl,nl), aa4121(nl,nl), aa4131(nl,nl)
55	real*8 aa12(nl)
56	real*8 aa1211(nl,nl), aa1221(nl,nl), aa1231(nl,nl), aa1241(nl,nl)
57	real*8 aa1112(nl,nl), aa2112(nl,nl), aa3112(nl,nl), aa4112(nl,nl)
58
59	real*8 aaalf11(nl,nl),aaalf21(nl,nl),aaalf31(nl,nl),aaalf41(nl,nl)
60	real*8 aaa11(nl),aaa1121(nl,nl),aaa1131(nl,nl),aaa1141(nl,nl)
61	real*8 aaa21(nl),aaa2111(nl,nl),aaa2131(nl,nl),aaa2141(nl,nl)
62	real*8 aaa31(nl),aaa3111(nl,nl),aaa3121(nl,nl),aaa3141(nl,nl)
63	real*8 aaa41(nl),aaa4111(nl,nl),aaa4121(nl,nl),aaa4131(nl,nl)
64
65	real*8 aaaalf11(nl,nl),aaaalf41(nl,nl)
66	real*8 aaaa11(nl),aaaa1141(nl,nl)
67	real*8 aaaa41(nl),aaaa4111(nl,nl)
68
69
70
71	! populations
72	real*8 n10(nl), n11(nl)
73	real*8 n20(nl), n21(nl)
74	real*8 n30(nl), n31(nl)
75	real*8 n40(nl), n41(nl)
76
77
78	! productions and loses
79	real*8 d19a1,d19b1,d19c1
80	real*8 d19ap1,d19bp1,d19cp1
81	real*8 d19a2,d19b2,d19c2
82	real*8 d19ap2,d19bp2,d19cp2
83	real*8 d19a3,d19b3,d19c3
84	real*8 d19ap3,d19bp3,d19cp3
85	real*8 d19a4,d19b4,d19c4
86	real*8 d19ap4,d19bp4,d19cp4
87
88	real*8 l11, l12, l21, l31, l41
89	real*8 p11, p12, p21, p31, p41
90	real*8 p1112, p1211, p1221, p1231, p1241
91	real*8 p1121, p1131, p1141
92	real*8 p2111, p2112, p2131, p2141
93	real*8 p3111, p3112, p3121, p3141
94	real*8 p4111, p4112, p4121, p4131
95
96
97	real*8 ps11, ps21, ps31, ps41, ps12
98
99	real*8 pl11, pl12, pl21, pl31, pl41
100
101	c local constants and indexes
102
103	integer ii ! decides if output of tv,hr
104	integer icurt ! decides if read/comp c.matrix
105
106	real*8 co2t
107	real*8 ftest
108
109	real*8 a11_einst(nl), a12_einst(nl)
110	real*8 a21_einst(nl), a31_einst(nl), a41_einst(nl)
111	real tsurf
112
113	real*8 nu11, nu12, nu121, nu21, nu31, nu41
114
115	integer i, j, ik, isot , icurtishb
116	integer i_by15sh, i_col020, i_col010636
117
118
119	c external functions and subroutines
120
121	external planckdp
122	real*8 planckdp
123
124	! subroutines called:
125	! mz4sub, dmzout, readc_mz4, mztf
126
127	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! start program
128
129
130	ii = 4
131	icurt = 1
132
133	call zero4v( aa11, aa21, aa31, aa41, nl)
134	call zero4m( aa1121, aa1131, aa1141, aalf11, nl)
135	call zero4m( aa2111, aa2131, aa2141, aalf21, nl)
136	call zero4m( aa3111, aa3121, aa3141, aalf31, nl)
137	call zero4m( aa4111, aa4121, aa4131, aalf41, nl)
138	call zero4m( aa1112, aa2112, aa3112, aa4112, nl)
139	call zero4m( aa1211, aa1221, aa1231, aa1241, nl)
140	call zero3v( aaa41, aaa31, aaa11, nl )
141	call zero3m( aaa4111, aaa4131, aaalf41, nl)
142	call zero3m( aaa3111, aaa3141, aaalf31, nl)
143	call zero3m( aaa1131, aaa1141, aaalf11, nl)
144	call zero2v( aaaa11, aaaa41, nl )
145	call zero2m( aaaa1141, aaaalf11, nl)
146	call zero2m( aaaa4111, aaaalf41, nl)
147
148	!write (,) ' --- c z a simple --- input_cza : ', input_cza
149
150
151	call zero3v (vt11,vt12,vt13,nl)
152	call zero3v (vt21,vt22,vt23,nl)
153	call zero3v (vt31,vt32,vt33,nl)
154	call zero3v (vt41,vt42,vt43,nl)
155
156	call zero3v (hr110,hr121,hr132,nl)
157	call zero3v (hr210,hr221,hr232,nl)
158	call zero3v (hr310,hr321,hr332,nl)
159	call zero3v (hr410,hr421,hr432,nl)
160	call zero3v (sl110,sl121,sl132,nl)
161	call zero3v (sl210,sl221,sl232,nl)
162	call zero3v (sl310,sl321,sl332,nl)
163	call zero3v (sl410,sl421,sl432,nl)
164
165	call zero4v (el11,el21,el31,el41,nl)
166	call zero4v (e110,e210,e310,e410,nl)
167	call zero3v (el12,e121,e112,nl)
168
169	call zero3m (cax1,cax2,cax3,nl)
170	call zerom (f1,nl)
171	call zero3v (v1,v2,v3,nl)
172
173	call zero4m (alf11,alf21,alf31,alf41,nl)
174	call zerom (alf12,nl)
175	call zero2v (a11,a12,nl)
176	call zero3v (a21,a31,a41,nl)
177
178	call zero3m (a1121,a1131,a1141,nl)
179	call zerom (a1112,nl)
180
181	call zero3m (a1221,a1231,a1241,nl)
182	call zerom (a1211,nl)
183
184	call zero2m (a2111,a2112,nl)
185	call zero2m (a2131,a2141,nl)
186	call zero2m (a3111,a3112,nl)
187	call zero2m (a3121,a3141,nl)
188	call zero2m (a4111,a4112,nl)
189	call zero2m (a4121,a4131,nl)
190
191
192	call zero4v (n11,n21,n31,n41,nl)
193
194	nu11 = nu(1,1)
195	nu12 = nu(1,2)
196	nu121 = nu12-nu11
197
198	nu21 = nu(2,1)
199
200	nu31 = nu(3,1)
201
202	nu41 = nu(4,1)
203
204	ftest = 1.d0
205	i_by15sh = 1
206	i_col020 = 1
207
208	i_col010636 = 1
209
210
211	101 format(a1)
212	180 format(a80)
213
214
215	c establishing molecular populations needed as input
216	do i=1,nl
217	n10(i) = dble( co2(i) * imr(1) )
218	n20(i) = dble( co2(i) * imr(2) )
219	n30(i) = dble( co2(i) * imr(3) )
220	n40(i) = dble( co2(i) * imr(4) )
221	if ( input_cza.ge.1 ) then
222	n11(i) = n10(i) 2.d0 exp( dble(-ee*nu(1,1))/v626t1(i) )
223	n21(i) = n20(i) 2.d0 exp( dble(-ee*nu(2,1))/v628t1(i) )
224	n31(i) = n30(i) 2.d0 exp( dble(-ee*nu(3,1))/v636t1(i) )
225	n41(i) = n40(i) 2.d0 exp( dble(-ee*nu(4,1))/v627t1(i) )
226	end if
227	enddo
228
229	cc
230	cc curtis matrix calculation
231	cc
232	if ( input_cza.ge.1 ) then
233
234	if (itt_cza.eq.15 ) then
235
236	call cm15um_hb_simple ( ig,icurt )
237
238	elseif (itt_cza.eq.13) then
239
240	call mztvc_626fh(ig)
241
242	endif
243
244	endif
245
246
247
248	do 4,i=nl,1,-1 !----------------------------------------------
249
250	co2t = dble ( co2(i) *(imr(1)+imr(3)+imr(2)+imr(4)) )
251
252	call getk ( t(i) )
253
254	ps11 = 0.d0
255	ps21 = 0.d0
256	ps31 = 0.d0
257	ps41 = 0.d0
258	ps12 = 0.d0
259
260	! V-T productions and losses V-T
261
262	isot = 1
263	d19b1 = dble(k19ba(isot)co2t+k19bb(isot)n2(i))
264	@ + dble(k19bc(isot)*co(i))
265	d19c1 = dble(k19ca(isot)co2t+k19cb(isot)n2(i))
266	@ + dble(k19cc(isot)*co(i))
267	d19bp1 = dble( k19bap(isot)co2t + k19bbp(isot)n2(i) )
268	@ + dble( k19bcp(isot)*co(i) )
269	d19cp1 = dble( k19cap(isot)co2t + k19cbp(isot)n2(i) )
270	@ + dble( k19ccp(isot)*co(i) )
271	isot = 2
272	d19c2 =dble(k19ca(isot)co2t+k19cb(isot)n2(i))
273	@ + dble(k19cc(isot)*co(i))
274	d19cp2 =dble( k19cap(isot)co2t + k19cbp(isot)n2(i) )
275	@ + dble( k19ccp(isot)*co(i) )
276	isot = 3
277	d19c3 =dble(k19ca(isot)co2t+k19cb(isot)n2(i))
278	@ + dble(k19cc(isot)*co(i))
279	d19cp3 =dble( k19cap(isot)co2t + k19cbp(isot)n2(i) )
280	@ + dble( k19ccp(isot)*co(i) )
281	isot = 4
282	d19c4 =dble(k19ca(isot)co2t+k19cb(isot)n2(i))
283	@ + dble(k19cc(isot)*co(i))
284	d19cp4 =dble( k19cap(isot)co2t + k19cbp(isot)n2(i) )
285	@ + dble(k19ccp(isot)*co(i) )
286	!
287	l11 = d19c1 + k20c(1)*dble(o3p(i))
288	p11 = ( d19cp1 + k20cp(1)dble(o3p(i)) ) n10(i)
289	l21 = d19c2 + k20c(2)*dble(o3p(i))
290	p21 = ( d19cp2 + k20cp(2)dble(o3p(i)) ) n20(i)
291	l31 = d19c3 + k20c(3)*dble(o3p(i))
292	p31 = ( d19cp3 + k20cp(3)dble(o3p(i)) ) n30(i)
293	l41 = d19c4 + k20c(4)*dble(o3p(i))
294	p41 = ( d19cp4 + k20cp(4)dble(o3p(i)) ) n40(i)
295
296	! Addition of V-V
297
298	l11 = l11 + k21cp(2)n20(i) + k21cp(3)n30(i) + k21cp(4)*n40(i)
299	p1121 = k21c(2) * n10(i)
300	p1131 = k21c(3) * n10(i)
301	p1141 = k21c(4) * n10(i)
302	!
303	l21 = l21 + k21c(2)n10(i) + k23k21cn30(i) + k24k21c*n40(i)
304	p2111 = k21cp(2) * n20(i)
305	p2131 = k23k21cp * n20(i)
306	p2141 = k24k21cp * n20(i)
307	!
308	l31 = l31 + k21c(3)n10(i) + k23k21cpn20(i) + k34k21c*n40(i)
309	p3111 = k21cp(3)* n30(i)
310	p3121 = k23k21c * n30(i)
311	p3141 = k34k21cp* n30(i)
312	!
313	l41 = l41 + k21c(4)n10(i) + k24k21cpn20(i) + k34k21cp*n30(i)
314	p4111 = k21cp(4)* n40(i)
315	p4121 = k24k21c * n40(i)
316	p4131 = k34k21c * n40(i)
317
318
319	if ( input_cza.ge.1 ) then
320
321	l12 = d19b1
322	@ + k20b(1)*dble(o3p(i))
323	@ + k21b(1)*n10(i)
324	@ + k33c*( n20(i) + n30(i) + n40(i) )
325	p12 = k21bp(1)n11(i) n11(i)
326	p1211 = d19bp1 + k20bp(1)*dble(o3p(i))
327	p1221 = k33cp(2)*n11(i)
328	p1231 = k33cp(3)*n11(i)
329	p1241 = k33cp(4)*n11(i)
330
331	l11 = l11 + d19bp1
332	@ + k20bp(1)*dble(o3p(i))
333	@ + 2.d0 * k21bp(1) * n11(i)
334	@ + k33cp(2)n21(i) + k33cp(3)n31(i) + k33cp(4)*n41(i)
335	p1112 = d19b1
336	@ + k20b(1)*dble(o3p(i))
337	@ + 2.d0k21b(1)n10(i)
338	@ + k33c*( n20(i) + n30(i) + n40(i) )
339
340	l21 = l21 + k33cp(2)*n11(i)
341	p2112 = k33c*n20(i)
342
343	l31 = l31 + k33cp(3)*n11(i)
344	p3112 = k33c*n30(i)
345
346	l41 = l41 + k33cp(4)*n11(i)
347	p4112 = k33c*n40(i)
348
349	end if
350
351
352	! Changes in local losses for ITT=13,15 cases
353
354	a21_einst(i) = 1.3452d00 * 1.8 / 4.0 * taustar21(i)
355	a31_einst(i) = 1.1878d00 * 1.8 / 4.0 * taustar31(i)
356	a41_einst(i) = 1.2455d00 * 1.8 / 4.0 * taustar41(i)
357
358	l21 = l21 + a21_einst(i)
359	l31 = l31 + a31_einst(i)
360	l41 = l41 + a41_einst(i)
361
362	if (input_cza.ge.1 .and. itt_cza.eq.13) then
363	a12_einst(i) = 4.35d00 / 3.0d0 * 1.8 / 4.0 * taustar12(i)
364	l12=l12+a12_einst(i)
365	endif
366
367	if (itt_cza.eq.24) then
368	a11_einst(i) = a11_einst(i) * 1.8 / 4.0 * taustar11(i)
369	l11 = l11 + a11_einst(i)
370	endif
371
372
373	! vectors and matrices for the formulation
374
375	a11(i) = dble(gammanu113.) 1.d0/2.d0 * (p11+ps11) /
376	@ (n10(i)*l11)
377	a1121(i,i) = dble((nu11/nu21))*3.d0 n20(i)/n10(i) * p1121/l11
378	a1131(i,i) = dble((nu11/nu31))*3.d0 n30(i)/n10(i) * p1131/l11
379	a1141(i,i) = dble((nu11/nu41))*3.d0 n40(i)/n10(i) * p1141/l11
380	e110(i) = 2.d0* dble(vlightnu112.) 1.d0/2.d0 /
381	@ ( n10(i) * l11 )
382
383	a21(i) = dble( gammanu213.) 1.d0/2.d0 *
384	@ (p21+ps21)/(n20(i)*l21)
385	a2111(i,i) = dble((nu21/nu11))*3.d0 n10(i)/n20(i) * p2111/l21
386	a2131(i,i) = dble((nu21/nu31))*3.d0 n30(i)/n20(i) * p2131/l21
387	a2141(i,i) = dble((nu21/nu41))*3.d0 n40(i)/n20(i) * p2141/l21
388	e210(i) = 2.d0dble(vlightnu21*2.) 1.d0/2.d0 /
389	@ ( n20(i) * l21 )
390
391	a31(i) = dble(gammanu313.) 1.d0/2.d0 * (p31+ps31) /
392	@ (n30(i)*l31)
393	a3111(i,i) = dble((nu31/nu11))*3.d0 n10(i)/n30(i) * p3111/l31
394	a3121(i,i) = dble((nu31/nu21))*3.d0 n20(i)/n30(i) * p3121/l31
395	a3141(i,i) = dble((nu31/nu41))*3.d0 n40(i)/n30(i) * p3141/l31
396	e310(i) = 2.d0dble(vlightnu31*2.) 1.d0/2.d0 /
397	@ ( n30(i) * l31 )
398
399	a41(i) = dble(gammanu413.) 1.d0/2.d0 * (p41+ps41) /
400	@ (n40(i)*l41)
401	a4111(i,i) = dble((nu41/nu11))*3.d0 n10(i)/n40(i) * p4111/l41
402	a4121(i,i) = dble((nu41/nu21))*3.d0 n20(i)/n40(i) * p4121/l41
403	a4131(i,i) = dble((nu41/nu31))*3.d0 n30(i)/n40(i) * p4131/l41
404	e410(i) = 2.d0dble(vlightnu41*2.) 1.d0/2.d0 /
405	@ ( n40(i) * l41 )
406
407	if (input_cza.ge.1) then
408
409	a1112(i,i) = dble((nu11/nu121))*3.d0 n11(i)/n10(i) *
410	@ p1112/l11
411	a2112(i,i) = dble((nu21/nu121))*3.d0 n11(i)/n20(i) *
412	@ p2112/l21
413	a3112(i,i) = dble((nu31/nu121))*3.d0 n11(i)/n30(i) *
414	@ p3112/l31
415	a4112(i,i) = dble((nu41/nu121))*3.d0 n11(i)/n40(i) *
416	@ p4112/l41
417	e112(i) = -2.d0dble(vlightnu11**3.)/nu121 /2.d0 /
418	@ ( n10(i)*l11 )
419	a12(i) = dble( gammanu1213.) 2.d0/4.d0* (p12+ps12)/
420	@ (n11(i)*l12)
421	a1211(i,i) = dble((nu121/nu11))*3.d0 n10(i)/n11(i) *
422	@ p1211/l12
423	a1221(i,i) = dble((nu121/nu21))*3.d0 n20(i)/n11(i) *
424	@ p1221/l12
425	a1231(i,i) = dble((nu121/nu31))*3.d0 n30(i)/n11(i) *
426	@ p1231/l12
427	a1241(i,i) = dble((nu121/nu41))*3.d0 n40(i)/n11(i) *
428	@ p1241/l12
429	e121(i) = 2.d0dble(vlightnu121*2.) 2.d0/4.d0 /
430	@ ( n11(i) * l12 )
431
432	end if
433
434
435	4 continue !-------------------------------------------------------
436
437
438	! Change C.M.
439
440	do i=1,nl
441	do j=1,nl
442	c210(i,j) = 0.0d0
443	c310(i,j) = 0.0d0
444	c410(i,j) = 0.0d0
445	end do
446	end do
447	if ( itt_cza.eq.13 ) then
448	do i=1,nl
449	do j=1,nl
450	c121(i,j) = 0.0d0
451	end do
452	end do
453	endif
454	!Añadido para hacer diagonal C121
455	! if ( itt_cza.eq.15 ) then
456	! do i=1,nl
457	! do j=1,nl
458	! if(abs(i-j).eq.1.or.abs(i-j).eq.2) c121(i,j) = 0.0d0
459	! end do
460	! end do
461	! endif
462	if ( itt_cza.eq.24 ) then
463	do i=1,nl
464	do j=1,nl
465	c110(i,j) = 0.0d0
466	end do
467	end do
468	endif
469
470	! Lower Boundary
471	tsurf = t(1) + tsurf_excess
472	do i=1,nl
473	sl110(i) = sl110(i) + vc110(i) * planckdp( tsurf, nu11 )
474	sl210(i) = sl210(i) + vc210(i) * planckdp( tsurf, nu21 )
475	sl310(i) = sl310(i) + vc310(i) * planckdp( tsurf, nu31 )
476	sl410(i) = sl410(i) + vc410(i) * planckdp( tsurf, nu41 )
477	end do
478	if (input_cza.ge.1) then
479	do i=1,nl
480	sl121(i) = sl121(i) + vc121(i) * planckdp( tsurf, nu121 )
481	end do
482	endif
483
484
485	!!!!!!!!!!!! Solucion del sistema
486
487	!! Paso 0 : Calculo de los alphas alf11, alf21, alf31, alf41, alf12
488
489	call unit ( cax2, nl )
490
491	call diago ( cax1, e110, nl )
492	call mulmm ( cax3, cax1,c110, nl )
493	! cax3=matmul(cax1,c110)
494	call resmm ( alf11, cax2,cax3, nl )
495
496	call diago ( cax1, e210, nl )
497	call mulmm ( cax3, cax1,c210, nl )
498	! cax3=matmul(cax1,c210)
499	call resmm ( alf21, cax2,cax3, nl )
500
501	call diago ( cax1, e310, nl )
502	call mulmm ( cax3, cax1,c310, nl )
503	! cax3=matmul(cax1,c310)
504	call resmm ( alf31, cax2,cax3, nl )
505	!
506	call diago ( cax1, e410, nl )
507	call mulmm ( cax3, cax1,c410, nl )
508	! cax3=matmul(cax1,c410)
509	call resmm ( alf41, cax2,cax3, nl )
510	!
511	! if(ig.eq.2223.and.input_cza.eq.1) then
512	! open(168,file='output_curtis_c121diagminus2.dat')
513	! do i=1,nl
514	! do j=1,nl
515	! write(168,*)i,j,c110(i,j),c121(i,j)
516	! enddo
517	! enddo
518	! close(168)
519	! open(178,file='output_taustar.dat')
520	! do i=1,nl
521	! write(178,*)i,taustar21(i),taustar31(i),taustar41(i)
522	! enddo
523	! close(178)
524	! endif
525	if (input_cza.ge.1) then
526	call diago ( cax1, e121, nl )
527	call mulmm ( cax3, cax1,c121, nl )
528	! cax3=matmul(cax1,c121)
529	call resmm ( alf12, cax2,cax3, nl )
530	endif
531
532	!! Paso 1 : Calculo de vectores y matrices con 1 barra (aa***)
533
534	if (input_cza.eq.0) then ! Skip paso 1, pues el12 no se calcula
535
536	! el11
537	call sypvvv( aa11, a11,e110,sl110, nl )
538	call samem( aa1121, a1121, nl )
539	call samem( aa1131, a1131, nl )
540	call samem( aa1141, a1141, nl )
541	call samem( aalf11, alf11, nl )
542
543	! el21
544	call sypvvv( aa21, a21,e210,sl210, nl )
545	call samem( aa2111, a2111, nl )
546	call samem( aa2131, a2131, nl )
547	call samem( aa2141, a2141, nl )
548	call samem( aalf21, alf21, nl )
549
550	! el31
551	call sypvvv( aa31, a31,e310,sl310, nl )
552	call samem( aa3111, a3111, nl )
553	call samem( aa3121, a3121, nl )
554	call samem( aa3141, a3141, nl )
555	call samem( aalf31, alf31, nl )
556
557	! el41
558	call sypvvv( aa41, a41,e410,sl410, nl )
559	call samem( aa4111, a4111, nl )
560	call samem( aa4121, a4121, nl )
561	call samem( aa4131, a4131, nl )
562	call samem( aalf41, alf41, nl )
563
564
565	else ! (input_cza.ge.1) , FH !
566
567
568	call sypvvv( v1, a12,e121,sl121, nl ) ! a12 + e121 * sl121
569
570	! aa11
571	call sypvvv( v2, a11,e110,sl110, nl )
572	call trucommvv( aa11 , alf12,a1112,v2, v1, nl )
573
574	! aalf11
575	call invdiag( cax1, a1112, nl )
576	call mulmm( cax2, alf12, cax1, nl ) ! alf12 * (1/a1112)
577	! cax2=matmul(alf12,cax1)
578	call mulmm( cax3, cax2, alf11, nl )
579	! cax3=matmul(cax2,alf11)
580
581	call resmm( aalf11, cax3, a1211, nl )
582	! aa1121
583	call trucodiag(aa1121, alf12,a1112,a1121, a1221, nl)
584	! aa1131
585	call trucodiag(aa1131, alf12,a1112,a1131, a1231, nl)
586	! aa1141
587	call trucodiag(aa1141, alf12,a1112,a1141, a1241, nl)
588
589
590	! aa21
591	call sypvvv( v2, a21,e210,sl210, nl )
592	call trucommvv( aa21 , alf12,a2112,v2, v1, nl )
593
594	! aalf21
595	call invdiag( cax1, a2112, nl )
596	call mulmm( cax2, alf12, cax1, nl ) ! alf12 * (1/a2112)
597	! cax2=matmul(alf12,cax1)
598	call mulmm( cax3, cax2, alf21, nl )
599	! cax3=matmul(cax2,alf21)
600	call resmm( aalf21, cax3, a1221, nl )
601	! aa2111
602	call trucodiag(aa2111, alf12,a2112,a2111, a1211, nl)
603	! aa2131
604	call trucodiag(aa2131, alf12,a2112,a2131, a1231, nl)
605	! aa2141
606	call trucodiag(aa2141, alf12,a2112,a2141, a1241, nl)
607
608
609	! aa31
610	call sypvvv( v2, a31,e310,sl310, nl )
611	call trucommvv( aa31 , alf12,a3112,v2, v1, nl )
612	! aalf31
613	call invdiag( cax1, a3112, nl )
614	call mulmm( cax2, alf12, cax1, nl ) ! alf12 * (1/a3112)
615	! cax2=matmul(alf12,cax1)
616	call mulmm( cax3, cax2, alf31, nl )
617	! cax3=matmul(cax2,alf31)
618	call resmm( aalf31, cax3, a1231, nl )
619	! aa3111
620	call trucodiag(aa3111, alf12,a3112,a3111, a1211, nl)
621	! aa3121
622	call trucodiag(aa3121, alf12,a3112,a3121, a1221, nl)
623	! aa3141
624	call trucodiag(aa3141, alf12,a3112,a3141, a1241, nl)
625
626
627	! aa41
628	call sypvvv( v2, a41,e410,sl410, nl )
629	call trucommvv( aa41 , alf12,a4112,v2, v1, nl )
630	! aalf41
631	call invdiag( cax1, a4112, nl )
632	call mulmm( cax2, alf12, cax1, nl ) ! alf12 * (1/a4112)
633	! cax2=matmul(alf12,cax1)
634	call mulmm( cax3, cax2, alf41, nl )
635	! cax3=matmul(cax2,alf41)
636	call resmm( aalf41, cax3, a1241, nl )
637	! aa4111
638	call trucodiag(aa4111, alf12,a4112,a4111, a1211, nl)
639	! aa4121
640	call trucodiag(aa4121, alf12,a4112,a4121, a1221, nl)
641	! aa4131
642	call trucodiag(aa4131, alf12,a4112,a4131, a1231, nl)
643
644	endif ! Final caso input_cza.ge.1
645
646
647	!! Paso 2 : Calculo de vectores y matrices con 2 barras (aaa***)
648
649	! aaalf41
650	call invdiag( cax1, aa4121, nl )
651	call mulmm( cax2, aalf21, cax1, nl ) ! alf21 * (1/a4121)
652	! cax2=matmul(aalf21,cax1)
653	call mulmm( cax3, cax2, aalf41, nl )
654	! cax3=matmul(cax2,aalf41)
655	call resmm( aaalf41, cax3, aa2141, nl )
656	! aaa41
657	call trucommvv(aaa41, aalf21,aa4121,aa41, aa21, nl)
658	! aaa4111
659	call trucodiag(aaa4111, aalf21,aa4121,aa4111, aa2111, nl)
660	! aaa4131
661	call trucodiag(aaa4131, aalf21,aa4121,aa4131, aa2131, nl)
662
663	! aaalf31
664	call invdiag( cax1, aa3121, nl )
665	call mulmm( cax2, aalf21, cax1, nl ) ! alf21 * (1/a3121)
666	! cax2=matmul(aalf21,cax1)
667	call mulmm( cax3, cax2, aalf31, nl )
668	! cax3=matmul(cax2,aalf31)
669	call resmm( aaalf31, cax3, aa2131, nl )
670	! aaa31
671	call trucommvv(aaa31, aalf21,aa3121,aa31, aa21, nl)
672	! aaa3111
673	call trucodiag(aaa3111, aalf21,aa3121,aa3111, aa2111, nl)
674	! aaa3141
675	call trucodiag(aaa3141, aalf21,aa3121,aa3141, aa2141, nl)
676
677	! aaalf11
678	call invdiag( cax1, aa1121, nl )
679	call mulmm( cax2, aalf21, cax1, nl ) ! alf21 * (1/a1121)
680	! cax2=matmul(aalf21,cax1)
681	call mulmm( cax3, cax2, aalf11, nl )
682	! cax3=matmul(cax2,aalf11)
683	call resmm( aaalf11, cax3, aa2111, nl )
684	! aaa11
685	call trucommvv(aaa11, aalf21,aa1121,aa11, aa21, nl)
686	! aaa1131
687	call trucodiag(aaa1131, aalf21,aa1121,aa1131, aa2131, nl)
688	! aaa1141
689	call trucodiag(aaa1141, aalf21,aa1121,aa1141, aa2141, nl)
690
691
692	!! Paso 3 : Calculo de vectores y matrices con 3 barras (aaaa***)
693
694	! aaaalf41
695	call invdiag( cax1, aaa4131, nl )
696	call mulmm( cax2, aaalf31, cax1, nl ) ! aaalf31 * (1/aaa4131)
697	! cax2=matmul(aaalf31,cax1)
698	call mulmm( cax3, cax2, aaalf41, nl )
699	! cax3=matmul(cax2,aaalf41)
700	call resmm( aaaalf41, cax3, aaa3141, nl )
701
702	! aaaa41
703	call trucommvv(aaaa41, aaalf31,aaa4131,aaa41, aaa31, nl)
704	! aaaa4111
705	call trucodiag(aaaa4111, aaalf31,aaa4131,aaa4111,aaa3111, nl)
706
707	! aaaalf11
708	call invdiag( cax1, aaa1131, nl )
709	call mulmm( cax2, aaalf31, cax1, nl ) ! aaalf31 * (1/aaa4131)
710	! cax2=matmul(aaalf31,cax1)
711	call mulmm( cax3, cax2, aaalf11, nl )
712	! cax3=matmul(cax2,aaalf11)
713	call resmm( aaaalf11, cax3, aaa3111, nl )
714	! aaaa11
715	call trucommvv(aaaa11, aaalf31,aaa1131,aaa11, aaa31, nl)
716	! aaaa1141
717	call trucodiag(aaaa1141, aaalf31,aaa1131,aaa1141,aaa3141, nl)
718
719
720	!! Paso 4 : Calculo de vectores y matrices finales y calculo de J1
721
722	call trucommvv(v1, aaaalf41,aaaa1141,aaaa11, aaaa41, nl)
723	!
724	call invdiag( cax1, aaaa1141, nl )
725	call mulmm( cax2, aaaalf41, cax1, nl ) ! aaaalf41 * (1/aaaa1141)
726	! cax2=matmul(aaaalf41,cax1)
727	call mulmm( cax3, cax2, aaaalf11, nl )
728	! cax3=matmul(cax2,aaaalf11)
729	call resmm( cax1, cax3, aaaa4111, nl )
730	!
731	call LUdec ( el11, cax1, v1, nl, nl2 )
732
733	! Solucion para el41
734	call sypvmv( v1, aaaa41, aaaa4111,el11, nl )
735	call LUdec ( el41, aaaalf41, v1, nl, nl2 )
736
737	! Solucion para el31
738	call sypvmv( v2, aaa31, aaa3111,el11, nl )
739	call sypvmv( v1, v2, aaa3141,el41, nl )
740	call LUdec ( el31, aaalf31, v1, nl, nl2 )
741
742	! Solucion para el21
743	call sypvmv( v3, aa21, aa2111,el11, nl )
744	call sypvmv( v2, v3, aa2131,el31, nl )
745	call sypvmv( v1, v2, aa2141,el41, nl )
746	call LUdec ( el21, aalf21, v1, nl, nl2 )
747
748	!!!
749	el11(1) = planckdp( t(1), nu11 )
750	el21(1) = planckdp( t(1), nu21 )
751	el31(1) = planckdp( t(1), nu31 )
752	el41(1) = planckdp( t(1), nu41 )
753	el11(nl) = 2.d0 * el11(nl-1) - el11(nl2)
754	el21(nl) = 2.d0 * el21(nl-1) - el21(nl2)
755	el31(nl) = 2.d0 * el31(nl-1) - el31(nl2)
756	el41(nl) = 2.d0 * el41(nl-1) - el41(nl2)
757
758	call mulmv ( v1, c110,el11, nl )
759	call sumvv ( hr110, v1,sl110, nl )
760
761	! Solucion para el12
762	if (input_cza.ge.1) then
763
764	call sypvmv( v1, a12, a1211,el11, nl )
765	call sypvmv( v3, v1, a1221,el21, nl )
766	call sypvmv( v2, v3, a1231,el31, nl )
767	call sypvmv( v1, v2, a1241,el41, nl )
768	call LUdec ( el12, alf12, v1, nl, nl2 )
769
770	el12(1) = planckdp( t(1), nu121 )
771	el12(nl) = 2.d0 * el12(nl-1) - el12(nl2)
772
773	if (itt_cza.eq.15) then
774	call mulmv ( v1, c121,el12, nl )
775	call sumvv ( hr121, v1,sl121, nl )
776	endif
777
778	end if
779
780
781
782	if (input_cza.lt.1) then
783
784	do i=1,nl
785	pl11 = el11(i)/dble( gamma * nu11*3.0d0 1./2. / n10(i) )
786	pl21 = el21(i)/dble( gamma * nu21*3.0d0 1./2. / n20(i) )
787	pl31 = el31(i)/dble( gamma * nu31*3.0d0 1./2. / n30(i) )
788	pl41 = el41(i)/dble( gamma * nu41*3.0d0 1./2. / n40(i) )
789	vt11(i) = dble(-eenu11) / log( abs(pl11) / (2.0d0n10(i)) )
790	vt21(i) = dble(-eenu21) / log( abs(pl21) / (2.0d0n20(i)) )
791	vt31(i) = dble(-eenu31) / log( abs(pl31) / (2.0d0n30(i)) )
792	vt41(i) = dble(-eenu41) / log( abs(pl41) / (2.0d0n40(i)) )
793	hr210(i) = sl210(i) - hplanckvlightnu21 * a21_einst(i)*pl21
794	hr310(i) = sl310(i) - hplanckvlightnu31 * a31_einst(i)*pl31
795	hr410(i) = sl410(i) - hplanckvlightnu41 * a41_einst(i)*pl41
796	! hr410(i) = 0.
797	enddo
798
799	call dinterconnection ( v626t1, vt11 )
800	call dinterconnection ( v628t1, vt21 )
801	call dinterconnection ( v636t1, vt31 )
802	call dinterconnection ( v627t1, vt41 )
803
804	else
805
806	do i=1,nl
807	pl21 = el21(i)/dble( gamma * nu21*3.0d0 1./2. / n20(i) )
808	pl31 = el31(i)/dble( gamma * nu31*3.0d0 1./2. / n30(i) )
809	pl41 = el41(i)/dble( gamma * nu41*3.0d0 1./2. / n40(i) )
810	hr210(i) = sl210(i) - hplanckvlightnu21 * a21_einst(i)*pl21
811	hr310(i) = sl310(i) - hplanckvlightnu31 * a31_einst(i)*pl31
812	hr410(i) = sl410(i) - hplanckvlightnu41 * a41_einst(i)*pl41
813	! hr410(i) = 0.
814	if (itt_cza.eq.13) then
815	pl12 = el12(i)/dble( gamma * nu121*3.0d0 2./4. / n11(i) )
816	hr121(i) = - hplanckvlight nu121 * a12_einst(i) * pl12
817	hr121(i) = hr121(i) + sl121(i)
818	endif
819	enddo
820
821	endif
822
823	! K/Dday
824	do i=1,nl
825	hr110(i)=hr110(i)*( hrkday_factor(i) / nt(i) )
826	hr210(i)=hr210(i)*( hrkday_factor(i) / nt(i) )
827	hr310(i)=hr310(i)*( hrkday_factor(i) / nt(i) )
828	hr410(i)=hr410(i)*( hrkday_factor(i) / nt(i) )
829	hr121(i)=hr121(i)*( hrkday_factor(i) / nt(i) )
830	end do
831
832
833
834	c output
835
836	!codigo = codeout
837	!call dmzout_tv ( 1 )
838	!call dmzout_hr ( 1 )
839
840	c final subrutina
841	return
842	end

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

Download in other formats:

Original Format