Context Navigation

NLTEdlvr09_CZALU.F @ 473

Last change on this file since 473 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

Rev	Line
[414]	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