Context Navigation

nlte_aux.F @ 2203

Last change on this file since 2203 was 2048, checked in by slebonnois, 6 years ago
SL: VENUS, autres details
File size: 65.8 KB

Rev	Line
[1310]	1	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	2	! Fast scheme for NLTE cooling rates at 15um by CO2 in a Martian GCM !
	3	! Version dlvr11_03. 2012. !
	4	! Software written and provided by IAA/CSIC, Granada, Spain, !
	5	! under ESA contract "Mars Climate Database and Physical Models" !
	6	! Person of contact: Miguel Angel Lopez Valverde valverde@iaa.es !
	7	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	8	c**********************************************************************
	9
	10	c Includes the following old 1-D model files/subroutines
	11
	12	c -MZTCRSUB_dlvr11.f
	13	c *dinterconnection
	14	c *planckd
	15	c *leetvt
	16	c -MZTFSUB_dlvr11_02.f
	17	c *initial
	18	c *intershphunt
	19	c *interstrhunt
	20	c *intzhunt
	21	c *intzhunt_cts
	22	c *rhist
	23	c *we_clean
	24	c *mztf_correccion
	25	c *mzescape_normaliz
	26	c *mzescape_normaliz_02
	27	c -interdpESCTVCISO_dlvr11.f
	28	c -hunt_cts.f
	29	c -huntdp.f
	30	c -hunt.f
	31	c -interdp_limits.f
	32	c -interhunt2veces.f
	33	c -interhunt5veces.f
	34	c -interhuntdp3veces.f
	35	c -interhuntdp4veces.f
	36	c -interhuntdp.f
	37	c -interhunt.f
	38	c -interhuntlimits2veces.f
	39	c -interhuntlimits5veces.f
	40	c -interhuntlimits.f
	41	c -lubksb_dp.f
	42	c -ludcmp_dp.f
	43	c -LUdec.f
	44	c -mat_oper.f
	45	c *unit
	46	c *diago
	47	c *invdiag
	48	c *samem
	49	c *mulmv
	50	c *trucodiag
	51	c *trucommvv
	52	c *sypvmv
	53	c *mulmm
	54	c *resmm
	55	c *sumvv
	56	c *sypvvv
	57	c *zerom
	58	c *zero4m
	59	c *zero3m
	60	c *zero2m
	61	c *zerov
	62	c *zero4v
	63	c *zero3v
	64	c *zero2v
	65	c -suaviza.f
	66
	67	c**********************************************************************
	68
	69
	70	c * Old MZTCRSUB_dlvr11.f *
	71
	72	!************************************************************************
	73
	74	! subroutine dinterconnection ( v, vt )
	75
	76
	77	************************************************************************
	78
	79	! implicit none
	80	!#include "nlte_paramdef.h"
	81
	82	c argumentos
	83	! real*8 vt(nl), v(nl)
	84
	85	c local variables
	86	! integer i
	87
	88	c *************
	89	!
	90	! do i=1,nl
	91	! v(i) = vt(i)
	92	! end do
	93
	94	! return
	95	! end
	96
	97	c***********************************************************************
	98	function planckdp(tp,xnu)
	99	c***********************************************************************
	100
	101	implicit none
	102
	103	#include "nlte_paramdef.h"
	104
	105	real*8 planckdp
	106	real*8 xnu
	107	real tp
	108
	109	planckdp = gammaxnu3.0d0 / exp( eexnu/dble(tp) )
	110	!erg cm-2.sr-1/cm-1.
	111
	112	c end
	113	return
	114	end
	115
	116	c***********************************************************************
	117	subroutine leetvt
	118
	119	c***********************************************************************
	120
	121	implicit none
	122
	123	#include "nlte_paramdef.h"
	124	#include "nlte_commons.h"
	125
	126	c local variables
	127	integer i
	128	real*8 zld(nl), zyd(nzy)
	129	real*8 xvt11(nzy), xvt21(nzy), xvt31(nzy), xvt41(nzy)
	130
	131	c***********************************************************************
	132
	133	do i=1,nzy
	134	zyd(i) = dble(zy(i))
	135	xvt11(i)= dble( ty(i) )
	136	xvt21(i)= dble( ty(i) )
	137	xvt31(i)= dble( ty(i) )
	138	xvt41(i)= dble( ty(i) )
	139	end do
	140
	141	do i=1,nl
	142	zld(i) = dble( zl(i) )
	143	enddo
	144	call interhuntdp4veces ( v626t1,v628t1,v636t1,v627t1, zld,nl,
	145	$ xvt11, xvt21, xvt31, xvt41, zyd,nzy, 1 )
	146
	147
	148	c end
	149	return
	150	end
	151
	152
	153	c * MZTFSUB_dlvr11_02.f *
	154
	155
	156	c ****************************************************************
	157	subroutine initial
	158
	159	c ****************************************************************
	160
	161	implicit none
	162
	163	#include "nlte_paramdef.h"
	164	#include "nlte_commons.h"
	165
	166	c local variables
	167	integer i
	168
	169	c ***************
	170
	171	eqw = 0.0d00
	172	aa = 0.0d00
	173	cc = 0.0d00
	174	dd = 0.0d00
	175
	176	do i=1,nbox
	177	ccbox(i) = 0.0d0
	178	ddbox(i) = 0.0d0
	179	end do
	180
	181	return
	182	end
	183
	184	c **********************************************************************
	185
	186	subroutine intershphunt (i, alsx,adx,xtemp)
	187
	188	c **********************************************************************
	189
	190	implicit none
	191
	192	#include "nlte_paramdef.h"
	193	#include "nlte_commons.h"
	194
	195	c arguments
	196	real*8 alsx(nbox_max),adx(nbox_max) ! Output
	197	real*8 xtemp(nbox_max) ! Input
	198	integer i ! I , O
	199
	200	c local variables
	201	integer k
	202	real*8 factor
	203	real*8 temperatura ! para evitar valores ligeramnt out of limits
[2048]	204	character(len=100) text
[1310]	205
	206	c ***********
	207
	208	do 1, k=1,nbox
	209	temperatura = xtemp(k)
	210	if (abs(xtemp(k)-thist(1)).le.0.01d0) then
	211	temperatura=thist(1)
	212	elseif (abs(xtemp(k)-thist(nhist)).le.0.01d0) then
	213	temperatura=thist(nhist)
	214	elseif (xtemp(k).lt.thist(1)) then
	215	temperatura=thist(1)
	216	c write (,) ' WARNING intershphunt/ Too low atmosph Tk:'
	217	c write (,) ' WARNING k,xtemp = ', k,xtemp(k)
	218	c write (,) ' Minimum Tk in histogram used : ', thist(1)
	219	elseif (xtemp(k).gt.thist(nhist)) then
	220	temperatura=thist(nhist)
	221	c write (,) ' WARNING intershphunt/ Very high atmosph Tk:'
	222	c write (,) ' WARNING k,xtemp = ', k,xtemp(k)
	223	c write (,) ' Max Tk in histogram used : ', thist(nhist)
	224	endif
	225	call huntdp ( thist,nhist, temperatura, i )
	226	if ( i.eq.0 .or. i.eq.nhist ) then
	227	c write (,) ' HUNT/ Limits input grid:',
	228	c @ thist(1),thist(nhist)
	229	c write (,) ' HUNT/ location in grid:', xtemp(k)
[2048]	230	text='INTERSHP/ Interpolation error. T out of Histogram.'
	231	CALL abort_physic("nlte_aux",text, 1)
[1310]	232	endif
	233	factor = 1.d0 / (thist(i+1)-thist(i))
	234	alsx(k) = (( xls1(i,k)*(thist(i+1)-xtemp(k)) +
	235	@ xls1(i+1,k)(xtemp(k)-thist(i)) )) factor
	236	adx(k) = (( xld1(i,k)*(thist(i+1)-xtemp(k)) +
	237	@ xld1(i+1,k)(xtemp(k)-thist(i)) )) factor
	238	1 continue
	239
	240	return
	241	end
	242
	243	c **********************************************************************
	244
	245	subroutine interstrhunt (i, stx, ts, sx, xtemp )
	246
	247	c **********************************************************************
	248
	249	implicit none
	250
	251	#include "nlte_paramdef.h"
	252	#include "nlte_commons.h"
	253
	254	c arguments
	255	real*8 stx ! output, total band strength
	256	real*8 ts ! input, temp for stx
	257	real*8 sx(nbox_max) ! output, strength for each box
	258	real*8 xtemp(nbox_max) ! input, temp for sx
	259	integer i
	260
	261	c local variables
	262	integer k
	263	real*8 factor
	264	real*8 temperatura
[2048]	265	character(len=100) text
[1310]	266
	267	c ***********
	268
	269	do 1, k=1,nbox
	270	temperatura = xtemp(k)
	271	if (abs(xtemp(k)-thist(1)).le.0.01d0) then
	272	temperatura=thist(1)
	273	elseif (abs(xtemp(k)-thist(nhist)).le.0.01d0) then
	274	temperatura=thist(nhist)
	275	elseif (xtemp(k).lt.thist(1)) then
	276	temperatura=thist(1)
	277	c write (,) ' WARNING interstrhunt/ Too low atmosph Tk:'
	278	c write (,) ' WARNING k,xtemp(k) = ', k,xtemp(k)
	279	c write (,) ' Minimum Tk in histogram used : ', thist(1)
	280	elseif (xtemp(k).gt.thist(nhist)) then
	281	temperatura=thist(nhist)
	282	c write (,) ' WARNING interstrhunt/ Very high atmosph Tk:'
	283	c write (,) ' WARNING k,xtemp(k) = ', k,xtemp(k)
	284	c write (,) ' Max Tk in histogram used : ', thist(nhist)
	285	endif
	286	call huntdp ( thist,nhist, temperatura, i )
	287	if ( i.eq.0 .or. i.eq.nhist ) then
	288	c write(,)'HUNT/ Limits input grid:',
	289	c $ thist(1),thist(nhist)
	290	c write(,)'HUNT/ location in grid:',xtemp(k)
[2048]	291	text='INTERSTR/1/ Interpolation error. T out of Histogram.'
	292	CALL abort_physic("nlte_aux",text, 1)
[1310]	293	endif
	294	factor = 1.d0 / (thist(i+1)-thist(i))
	295	sx(k) = ( sk1(i,k) * (thist(i+1)-xtemp(k))
	296	@ + sk1(i+1,k) * (xtemp(k)-thist(i)) ) * factor
	297	1 continue
	298
	299
	300	temperatura = ts
	301	if (abs(ts-thist(1)).le.0.01d0) then
	302	temperatura=thist(1)
	303	elseif (abs(ts-thist(nhist)).le.0.01d0) then
	304	temperatura=thist(nhist)
	305	elseif (ts.lt.thist(1)) then
	306	temperatura=thist(1)
	307	c write (,) ' WARNING interstrhunt/ Too low atmosph Tk:'
	308	c write (,) ' WARNING ts = ', temperatura
	309	c write (,) ' Minimum Tk in histogram used : ', thist(1)
	310	elseif (ts.gt.thist(nhist)) then
	311	temperatura=thist(nhist)
	312	c write (,) ' WARNING interstrhunt/ Very high atmosph Tk:'
	313	c write (,) ' WARNING ts = ', temperatura
	314	c write (,) ' Max Tk in histogram used : ', thist(nhist)
	315	endif
	316	call huntdp ( thist,nhist, temperatura, i )
	317	if ( i.eq.0 .or. i.eq.nhist ) then
	318	c write (,) ' HUNT/ Limits input grid:',
	319	c @ thist(1),thist(nhist)
	320	c write (,) ' HUNT/ location in grid:', ts
[2048]	321	text='INTERSTR/2/ Interpolat error. T out of Histogram.'
	322	CALL abort_physic("nlte_aux",text, 1)
[1310]	323	endif
	324	factor = 1.d0 / (thist(i+1)-thist(i))
	325	stx = 0.d0
	326	do k=1,nbox
	327	stx = stx + no(k) * ( sk1(i,k)*(thist(i+1)-ts) +
	328	@ sk1(i+1,k)(ts-thist(i)) ) factor
	329	end do
	330
	331
	332	return
	333	end
	334
	335	c **********************************************************************
	336
	337	subroutine intzhunt (k, h, aco2,ap,amr,at, con)
	338
	339	c k lleva la posicion de la ultima llamada a intz , necesario para
	340	c que esto represente una aceleracion real.
	341	c **********************************************************************
	342
	343	implicit none
	344	#include "nlte_paramdef.h"
	345	#include "nlte_commons.h"
	346
	347	c arguments
	348	real h ! i
	349	real*8 con(nzy) ! i
	350	real*8 aco2, ap, at, amr ! o
	351	integer k ! i
	352
	353	c local variables
	354	real factor
	355
	356	c ************
	357
	358	call hunt ( zy,nzy, h, k )
	359	factor = (h-zy(k)) / (zy(k+1)-zy(k))
	360	ap = dble( exp( log(py(k)) + log(py(k+1)/py(k)) * factor ) )
	361	aco2 = dlog(con(k)) + dlog( con(k+1)/con(k) ) * dble(factor)
	362	aco2 = exp( aco2 )
	363	at = dble( ty(k) + (ty(k+1)-ty(k)) * factor )
	364	amr = dble( mr(k) + (mr(k+1)-mr(k)) * factor )
	365
	366
	367	return
	368	end
	369
	370	c **********************************************************************
	371
	372	subroutine intzhunt_cts (k, h, nzy_cts_real,
	373	@ aco2,ap,amr,at, con)
	374
	375	c k lleva la posicion de la ultima llamada a intz , necesario para
	376	c que esto represente una aceleracion real.
	377	c **********************************************************************
	378
	379	implicit none
	380	#include "nlte_paramdef.h"
	381	#include "nlte_commons.h"
	382
	383	c arguments
	384	real h ! i
	385	real*8 con(nzy_cts) ! i
	386	real*8 aco2, ap, at, amr ! o
	387	integer k ! i
	388	integer nzy_cts_real ! i
	389
	390	c local variables
	391	real factor
	392
	393	c ************
	394
	395	call hunt_cts ( zy_cts,nzy_cts, nzy_cts_real, h, k )
	396	factor = (h-zy_cts(k)) / (zy_cts(k+1)-zy_cts(k))
	397	ap = dble( exp( log(py_cts(k)) +
	398	@ log(py_cts(k+1)/py_cts(k)) * factor ) )
	399	aco2 = dlog(con(k)) + dlog( con(k+1)/con(k) ) * dble(factor)
	400	aco2 = exp( aco2 )
	401	at = dble( ty_cts(k) + (ty_cts(k+1)-ty_cts(k)) * factor )
	402	amr = dble( mr_cts(k) + (mr_cts(k+1)-mr_cts(k)) * factor )
	403
	404
	405	return
	406	end
	407
	408
	409	c **********************************************************************
	410
	411	real*8 function we_clean ( y,pl, xalsa, xalda )
	412
	413	c **********************************************************************
	414
	415	implicit none
	416
	417	#include "nlte_paramdef.h"
	418
	419	c arguments
	420	real8 y ! I. path's absorber amount strength
	421	real*8 pl ! I. path's partial pressure of CO2
	422	real*8 xalsa ! I. Self lorentz linewidth for 1 isot & 1 box
	423	real*8 xalda ! I. Doppler linewidth " "
	424
	425	c local variables
	426	integer i
	427	real*8 x,wl,wd,wvoigt
	428	real*8 cn(0:7),dn(0:7)
	429	real*8 factor, denom
	430	real*8 pi, pi2, sqrtpi
	431
	432	c data blocks
	433	data cn/9.99998291698d-1,-3.53508187098d-1,9.60267807976d-2,
	434	@ -2.04969011013d-2,3.43927368627d-3,-4.27593051557d-4,
	435	@ 3.42209457833d-5,-1.28380804108d-6/
	436	data dn/1.99999898289,5.774919878d-1,-5.05367549898d-1,
	437	@ 8.21896973657d-1,-2.5222672453,6.1007027481,
	438	@ -8.51001627836,4.6535116765/
	439
	440	c ***********
	441
	442	pi = 3.141592
	443	pi2= 6.28318531
	444	sqrtpi = 1.77245385
	445
	446	x=y / ( pi2 * xalsa*pl )
	447
	448
	449	c Lorentz
	450	wl=y/sqrt(1.0d0+pi*x/2.0d0)
	451
	452	c Doppler
	453	x = y / (xalda*sqrtpi)
	454	if (x .lt. 5.0d0) then
	455	wd = cn(0)
	456	factor = 1.d0
	457	do i=1,7
	458	factor = factor * x
	459	wd = wd + cn(i) * factor
	460	end do
	461	wd = xalda * x * sqrtpi * wd
	462	else
	463	wd = dn(0)
	464	factor = 1.d0 / log(x)
	465	denom = 1.d0
	466	do i=1,7
	467	denom = denom * factor
	468	wd = wd + dn(i) * denom
	469	end do
	470	wd = xalda * sqrt(log(x)) * wd
	471	end if
	472
	473	c Voigt
	474	wvoigt = wlwl + wdwd - (wdwl/y)(wd*wl/y)
	475
	476	if ( wvoigt .lt. 0.0d0 ) then
	477	write (,) ' Subroutine WE/ Error in Voift EQS calculation'
	478	write (,) ' WL, WD, X, Y = ', wl, wd, x, y
	479	stop ' ERROR : Imaginary EQW. Revise spectral data. '
	480	endif
	481
	482	we_clean = sqrt( wvoigt )
	483
	484
	485	return
	486	end
	487
	488
	489	c ***********************************************************************
	490
	491	subroutine mztf_correccion (coninf, con, ib )
	492
	493	c ***********************************************************************
	494
	495	implicit none
	496
	497	#include "nlte_paramdef.h"
	498	#include "nlte_commons.h"
	499
	500	c arguments
	501	integer ib
	502	real*8 con(nzy), coninf
	503
	504	! local variables
	505	integer i, isot
	506	real*8 tvt0(nzy), tvtbs(nzy), zld(nl),zyd(nzy)
	507	real*8 xqv, xes, xlower, xfactor
	508
	509	c *********
	510
	511	isot = 1
	512	nu11 = dble( nu(1,1) )
	513
	514	do i=1,nzy
	515	zyd(i) = dble(zy(i))
	516	enddo
	517	do i=1,nl
	518	zld(i) = dble( zl(i) )
	519	end do
	520
	521	! tvtbs
	522	call interhuntdp (tvtbs,zyd,nzy, v626t1,zld,nl, 1 )
	523
	524	! tvt0
	525	if (ib.eq.2 .or. ib.eq.3 .or. ib.eq.4) then
	526	call interhuntdp (tvt0,zyd,nzy, v626t1,zld,nl, 1 )
	527	else
	528	do i=1,nzy
	529	tvt0(i) = dble( ty(i) )
	530	end do
	531	end if
	532
	533	c factor
	534	do i=1,nzy
	535
	536	xlower = exp( eedble(elow(isot,ib))
	537	@ ( 1.d0/dble(ty(i))-1.d0/tvt0(i) ) )
	538	xes = 1.0d0
	539	xqv = ( 1.d0-exp( -ee*nu11/tvtbs(i) ) ) /
	540	@ (1.d0-exp( -ee*nu11/dble(ty(i)) ))
	541	xfactor = xlower * xqv*2.d0 xes
	542
	543	con(i) = con(i) * xfactor
	544	if (i.eq.nzy) coninf = coninf * xfactor
	545
	546	end do
	547
	548
	549	return
	550	end
	551
	552
	553	c ***********************************************************************
	554
	555	subroutine mzescape_normaliz ( taustar, istyle )
	556
	557	c ***********************************************************************
	558
	559	implicit none
	560	#include "nlte_paramdef.h"
	561
	562	c arguments
	563	real*8 taustar(nl) ! o
	564	integer istyle ! i
	565
	566	c local variables and constants
	567	integer i, imaximum
	568	real*8 maximum
	569
	570	c ***************
	571
	572	taustar(nl) = taustar(nl-1)
	573
	574	if ( istyle .eq. 1 ) then
	575	imaximum = nl
	576	maximum = taustar(nl)
	577	do i=1,nl-1
	578	if (taustar(i).gt.maximum) taustar(i) = taustar(nl)
	579	enddo
	580	elseif ( istyle .eq. 2 ) then
	581	imaximum = nl
	582	maximum = taustar(nl)
	583	do i=nl-1,1,-1
	584	if (taustar(i).gt.maximum) then
	585	maximum = taustar(i)
	586	imaximum = i
	587	endif
	588	enddo
	589	do i=imaximum,nl
	590	if (taustar(i).lt.maximum) taustar(i) = maximum
	591	enddo
	592	endif
	593
	594	do i=1,nl
	595	taustar(i) = taustar(i) / maximum
	596	enddo
	597
	598
	599	c end
	600	return
	601	end
	602
	603	c ***********************************************************************
	604
	605	subroutine mzescape_normaliz_02 ( taustar, nn, istyle )
	606
	607	c ***********************************************************************
	608
	609	implicit none
	610
	611	c arguments
	612	real*8 taustar(nn) ! i,o
	613	integer istyle ! i
	614	integer nn ! i
	615
	616	c local variables and constants
	617	integer i, imaximum
	618	real*8 maximum
	619
	620	c ***************
	621
	622	taustar(nn) = taustar(nn-1)
	623
	624	if ( istyle .eq. 1 ) then
	625	imaximum = nn
	626	maximum = taustar(nn)
	627	do i=1,nn-1
	628	if (taustar(i).gt.maximum) taustar(i) = taustar(nn)
	629	enddo
	630	elseif ( istyle .eq. 2 ) then
	631	imaximum = nn
	632	maximum = taustar(nn)
	633	do i=nn-1,1,-1
	634	if (taustar(i).gt.maximum) then
	635	maximum = taustar(i)
	636	imaximum = i
	637	endif
	638	enddo
	639	do i=imaximum,nn
	640	if (taustar(i).lt.maximum) taustar(i) = maximum
	641	enddo
	642	endif
	643
	644	do i=1,nn
	645	taustar(i) = taustar(i) / maximum
	646	enddo
	647
	648
	649	c end
	650	return
	651	end
	652
	653
	654	c * interdp_ESCTVCISO_dlvr11.f *
	655
	656	c***********************************************************************
	657
	658	subroutine interdp_ESCTVCISO
	659
	660	c***********************************************************************
	661
	662	implicit none
	663
	664	#include "nlte_paramdef.h"
	665	#include "nlte_commons.h"
	666
	667	c local variables
	668	integer i
	669	real*8 lnpnb(nl)
	670
	671
	672	c***********************************************************************
	673
	674	c Use pressure in the NLTE grid but in log and in nb
	675	do i=1,nl
	676	lnpnb(i) = log( dble( pl(i) * 1013.25 * 1.e6) )
	677	enddo
	678
	679	c Interpolations
	680
	681	call interhuntdp3veces
	682	@ ( taustar21,taustar31,taustar41, lnpnb, nl,
	683	@ tstar21tab,tstar31tab,tstar41tab, lnpnbtab, nztabul,
	684	@ 1 )
	685
	686	call interhuntdp3veces ( vc210,vc310,vc410, lnpnb, nl,
	687	@ vc210tab,vc310tab,vc410tab, lnpnbtab, nztabul, 2 )
	688
	689	c end
	690	return
	691	end
	692
	693
	694	c * hunt_cts.f *
	695
	696	cccc
	697	SUBROUTINE hunt_cts(xx,n,n_cts,x,jlo)
	698	c
	699	c La dif con hunt es el uso de un indice superior (n_cts) mas bajito que (n)
	700	c
	701	c Arguments
	702	INTEGER jlo ! O
	703	INTEGER n ! I
	704	INTEGER n_cts ! I
	705	REAL xx(n) ! I
	706	REAL x ! I
	707
	708	c Local variables
	709	INTEGER inc,jhi,jm
	710	LOGICAL ascnd
	711	c
	712	cccc
	713	c
	714	ascnd=xx(n_cts).ge.xx(1)
	715	if(jlo.le.0.or.jlo.gt.n_cts)then
	716	jlo=0
	717	jhi=n_cts+1
	718	goto 3
	719	endif
	720	inc=1
	721	if(x.ge.xx(jlo).eqv.ascnd)then
	722	1 jhi=jlo+inc
	723	! write (,) jlo
	724	if(jhi.gt.n_cts)then
	725	jhi=n_cts+1
	726	! write (,) jhi-1
	727	else if(x.ge.xx(jhi).eqv.ascnd)then
	728	jlo=jhi
	729	inc=inc+inc
	730	! write (,) jlo
	731	goto 1
	732	endif
	733	else
	734	jhi=jlo
	735	2 jlo=jhi-inc
	736	! write (,) jlo
	737	if(jlo.lt.1)then
	738	jlo=0
	739	else if(x.lt.xx(jlo).eqv.ascnd)then
	740	jhi=jlo
	741	inc=inc+inc
	742	goto 2
	743	endif
	744	endif
	745	3 if(jhi-jlo.eq.1)then
	746	if(x.eq.xx(n_cts))jlo=n_cts-1
	747	if(x.eq.xx(1))jlo=1
	748	! write (,) jlo
	749	return
	750	endif
	751	jm=(jhi+jlo)/2
	752	if(x.ge.xx(jm).eqv.ascnd)then
	753	jlo=jm
	754	else
	755	jhi=jm
	756	endif
	757	! write (,) jhi-1
	758	goto 3
	759	c
	760	END
	761
	762
	763	c * huntdp.f *
	764
	765	cccc
	766	SUBROUTINE huntdp(xx,n,x,jlo)
	767	c
	768	c Arguments
	769	INTEGER jlo ! O
	770	INTEGER n ! I
	771	REAL*8 xx(n) ! I
	772	REAL*8 x ! I
	773
	774	c Local variables
	775	INTEGER inc,jhi,jm
	776	LOGICAL ascnd
	777	c
	778	cccc
	779	c
	780	ascnd=xx(n).ge.xx(1)
	781	if(jlo.le.0.or.jlo.gt.n)then
	782	jlo=0
	783	jhi=n+1
	784	goto 3
	785	endif
	786	inc=1
	787	if(x.ge.xx(jlo).eqv.ascnd)then
	788	1 jhi=jlo+inc
	789	if(jhi.gt.n)then
	790	jhi=n+1
	791	else if(x.ge.xx(jhi).eqv.ascnd)then
	792	jlo=jhi
	793	inc=inc+inc
	794	goto 1
	795	endif
	796	else
	797	jhi=jlo
	798	2 jlo=jhi-inc
	799	if(jlo.lt.1)then
	800	jlo=0
	801	else if(x.lt.xx(jlo).eqv.ascnd)then
	802	jhi=jlo
	803	inc=inc+inc
	804	goto 2
	805	endif
	806	endif
	807	3 if(jhi-jlo.eq.1)then
	808	if(x.eq.xx(n))jlo=n-1
	809	if(x.eq.xx(1))jlo=1
	810	return
	811	endif
	812	jm=(jhi+jlo)/2
	813	if(x.ge.xx(jm).eqv.ascnd)then
	814	jlo=jm
	815	else
	816	jhi=jm
	817	endif
	818	goto 3
	819	c
	820	END
	821
	822
	823	c * hunt.f *
	824
	825	cccc
	826	SUBROUTINE hunt(xx,n,x,jlo)
	827	c
	828	c Arguments
	829	INTEGER jlo ! O
	830	INTEGER n ! I
	831	REAL xx(n) ! I
	832	REAL x ! I
	833
	834	c Local variables
	835	INTEGER inc,jhi,jm
	836	LOGICAL ascnd
	837	c
	838	cccc
	839	c
	840	ascnd=xx(n).ge.xx(1)
	841	if(jlo.le.0.or.jlo.gt.n)then
	842	jlo=0
	843	jhi=n+1
	844	goto 3
	845	endif
	846	inc=1
	847	if(x.ge.xx(jlo).eqv.ascnd)then
	848	1 jhi=jlo+inc
	849	! write (,) jlo
	850	if(jhi.gt.n)then
	851	jhi=n+1
	852	! write (,) jhi-1
	853	else if(x.ge.xx(jhi).eqv.ascnd)then
	854	jlo=jhi
	855	inc=inc+inc
	856	! write (,) jlo
	857	goto 1
	858	endif
	859	else
	860	jhi=jlo
	861	2 jlo=jhi-inc
	862	! write (,) jlo
	863	if(jlo.lt.1)then
	864	jlo=0
	865	else if(x.lt.xx(jlo).eqv.ascnd)then
	866	jhi=jlo
	867	inc=inc+inc
	868	goto 2
	869	endif
	870	endif
	871	3 if(jhi-jlo.eq.1)then
	872	if(x.eq.xx(n))jlo=n-1
	873	if(x.eq.xx(1))jlo=1
	874	! write (,) jlo
	875	return
	876	endif
	877	jm=(jhi+jlo)/2
	878	if(x.ge.xx(jm).eqv.ascnd)then
	879	jlo=jm
	880	else
	881	jhi=jm
	882	endif
	883	! write (,) jhi-1
	884	goto 3
	885	c
	886	END
	887
	888
	889	c * interdp_limits.f *
	890
	891	c ***********************************************************************
	892
	893	subroutine interdp_limits ( yy, zz, m, i1,i2,
	894	@ y, z, n, j1,j2, opt)
	895
	896	c Interpolation soubroutine.
	897	c Returns values between indexes i1 & i2, donde 1 =< i1 =< i2 =< m
	898	c Solo usan los indices de los inputs entre j1,j2, 1 =< j1 =< j2 =< n
	899	c Input values: y(n) , z(n) (solo se usarann los valores entre j1,j2)
	900	c zz(m) (solo se necesita entre i1,i2)
	901	c Output values: yy(m) (solo se calculan entre i1,i2)
	902	c Options: opt=1 -> lineal ,, opt=2 -> logarithmic
	903	c Difference with interdp:
	904	c here interpolation proceeds between indexes i1,i2 only
	905	c if i1=1 & i2=m, both subroutines are exactly the same
	906	c thus previous calls to interdp or interdp2 could be easily replaced
	907
	908	c JAN 98 MALV Version for mz1d
	909	c ***********************************************************************
	910
	911	implicit none
	912
	913	! Arguments
	914	integer n,m ! I. Dimensions
	915	integer i1, i2, j1, j2, opt ! I
	916	real*8 zz(m) ! I
	917	real*8 yy(m) ! O
	918	real*8 z(n),y(n) ! I
	919
	920	! Local variables
	921	integer i,j
	922	real*8 zmin,zzmin,zmax,zzmax
	923
	924	c *******************************
	925
	926	! write (,) ' d interpolating '
	927	! call mindp_limits (z,n,zmin, j1,j2)
	928	! call mindp_limits (zz,m,zzmin, i1,i2)
	929	! call maxdp_limits (z,n,zmax, j1,j2)
	930	! call maxdp_limits (zz,m,zzmax, i1,i2)
	931	zmin=minval(z(j1:j2))
	932	zzmin=minval(zz(i1:i2))
	933	zmax=maxval(z(j1:j2))
	934	zzmax=maxval(zz(i1:i2))
	935
	936	if(zzmin.lt.zmin)then
	937	write (,) 'from d interp: new variable out of limits'
	938	write (,) zzmin,'must be .ge. ',zmin
	939	stop
	940	! elseif(zzmax.gt.zmax)then
	941	! type *,'from interp: new variable out of limits'
	942	! type *,zzmax, 'must be .le. ',zmax
	943	! stop
	944	end if
	945
	946	do 1,i=i1,i2
	947
	948	do 2,j=j1,j2-1
	949	if(zz(i).ge.z(j).and.zz(i).lt.z(j+1)) goto 3
	950	2 continue
	951	c in this case (zz(i2).eq.z(j2)) and j leaves the loop with j=j2-1+1=j2
	952	if(opt.eq.1)then
	953	yy(i)=y(j2-1)+(y(j2)-y(j2-1))*(zz(i)-z(j2-1))/
	954	$ (z(j2)-z(j2-1))
	955	elseif(opt.eq.2)then
	956	if(y(j2).eq.0.0d0.or.y(j2-1).eq.0.0d0)then
	957	yy(i)=0.0d0
	958	else
	959	yy(i)=exp(log(y(j2-1))+log(y(j2)/y(j2-1))*
	960	@ (zz(i)-z(j2-1))/(z(j2)-z(j2-1)))
	961	end if
	962	else
	963	write (,) ' d interp : opt must be 1 or 2, opt= ',opt
	964	end if
	965	goto 1
	966	3 continue
	967	if(opt.eq.1)then
	968	yy(i)=y(j)+(y(j+1)-y(j))*(zz(i)-z(j))/(z(j+1)-z(j))
	969	! type *, ' '
	970	! type *, ' z(j),z(j+1) =', z(j),z(j+1)
	971	! type *, ' t(j),t(j+1) =', y(j),y(j+1)
	972	! type *, ' zz, tt = ', zz(i), yy(i)
	973	elseif(opt.eq.2)then
	974	if(y(j+1).eq.0.0d0.or.y(j).eq.0.0d0)then
	975	yy(i)=0.0d0
	976	else
	977	yy(i)=exp(log(y(j))+log(y(j+1)/y(j))*
	978	@ (zz(i)-z(j))/(z(j+1)-z(j)))
	979	end if
	980	else
	981	write (,) ' interp : opt must be 1 or 2, opt= ',opt
	982	end if
	983	1 continue
	984	return
	985	end
	986
	987
	988
	989	c * interhunt2veces.f *
	990
	991	c ***********************************************************************
	992
	993	subroutine interhunt2veces ( y1,y2, zz,m,
	994	@ x1,x2, z,n, opt)
	995
	996	c interpolation soubroutine basada en Numerical Recipes HUNT.FOR
	997	c input values: y(n) at z(n)
	998	c output values: yy(m) at zz(m)
	999	c options: 1 -> lineal
	1000	c 2 -> logarithmic
	1001	c ***********************************************************************
	1002
	1003	implicit none
	1004
	1005	! Arguments
	1006	integer n,m,opt ! I
	1007	real zz(m),z(n) ! I
	1008	real y1(m),y2(m) ! O
	1009	real x1(n),x2(n) ! I
	1010
	1011
	1012	! Local variables
	1013	integer i, j
	1014	real factor
	1015	real zaux
[2048]	1016	character(len=100) text
[1310]	1017
	1018	!!!!
	1019
	1020	j = 1 ! initial first guess (=n/2 is anothr pssblty)
	1021
	1022	do 1,i=1,m !
	1023
	1024	! Busca indice j donde ocurre q zz(i) esta entre [z(j),z(j+1)]
	1025	zaux = zz(i)
	1026	if (abs(zaux-z(1)).le.0.01) then
	1027	zaux=z(1)
	1028	elseif (abs(zaux-z(n)).le.0.01) then
	1029	zaux=z(n)
	1030	endif
	1031	call hunt ( z,n, zaux, j )
	1032	if ( j.eq.0 .or. j.eq.n ) then
	1033	write (,) ' HUNT/ Limits input grid:', z(1),z(n)
	1034	write (,) ' HUNT/ location in new grid:', zz(i)
[2048]	1035	text='interhunt2/ Interpolat error. zz out of limits.'
	1036	CALL abort_physic("nlte_aux",text, 1)
[1310]	1037	endif
	1038
	1039	! Perform interpolation
	1040	factor = (zz(i)-z(j))/(z(j+1)-z(j))
	1041	if (opt.eq.1) then
	1042	y1(i) = x1(j) + (x1(j+1)-x1(j)) * factor
	1043	y2(i) = x2(j) + (x2(j+1)-x2(j)) * factor
	1044	else
	1045	y1(i) = exp( log(x1(j)) + log(x1(j+1)/x1(j)) * factor )
	1046	y2(i) = exp( log(x2(j)) + log(x2(j+1)/x2(j)) * factor )
	1047	end if
	1048
	1049	1 continue
	1050
	1051	return
	1052	end
	1053
	1054
	1055	c * interhunt5veces.f *
	1056
	1057	c ***********************************************************************
	1058
	1059	subroutine interhunt5veces ( y1,y2,y3,y4,y5, zz,m,
	1060	@ x1,x2,x3,x4,x5, z,n, opt)
	1061
	1062	c interpolation soubroutine basada en Numerical Recipes HUNT.FOR
	1063	c input values: y(n) at z(n)
	1064	c output values: yy(m) at zz(m)
	1065	c options: 1 -> lineal
	1066	c 2 -> logarithmic
	1067	c ***********************************************************************
	1068
	1069	implicit none
	1070	! Arguments
	1071	integer n,m,opt ! I
	1072	real zz(m),z(n) ! I
	1073	real y1(m),y2(m),y3(m),y4(m),y5(m) ! O
	1074	real x1(n),x2(n),x3(n),x4(n),x5(n) ! I
	1075
	1076
	1077	! Local variables
	1078	integer i, j
	1079	real factor
	1080	real zaux
[2048]	1081	character(len=100) text
[1310]	1082
	1083	!!!!
	1084
	1085	j = 1 ! initial first guess (=n/2 is anothr pssblty)
	1086
	1087	do 1,i=1,m !
	1088
	1089	! Busca indice j donde ocurre q zz(i) esta entre [z(j),z(j+1)]
	1090	zaux = zz(i)
	1091	if (abs(zaux-z(1)).le.0.01) then
	1092	zaux=z(1)
	1093	elseif (abs(zaux-z(n)).le.0.01) then
	1094	zaux=z(n)
	1095	endif
	1096	call hunt ( z,n, zaux, j )
	1097	if ( j.eq.0 .or. j.eq.n ) then
	1098	write (,) ' HUNT/ Limits input grid:', z(1),z(n)
	1099	write (,) ' HUNT/ location in new grid:', zz(i)
[2048]	1100	text='interhunt5/ Interpolat error. zz out of limits.'
	1101	CALL abort_physic("nlte_aux",text, 1)
[1310]	1102	endif
	1103
	1104	! Perform interpolation
	1105	factor = (zz(i)-z(j))/(z(j+1)-z(j))
	1106	if (opt.eq.1) then
	1107	y1(i) = x1(j) + (x1(j+1)-x1(j)) * factor
	1108	y2(i) = x2(j) + (x2(j+1)-x2(j)) * factor
	1109	y3(i) = x3(j) + (x3(j+1)-x3(j)) * factor
	1110	y4(i) = x4(j) + (x4(j+1)-x4(j)) * factor
	1111	y5(i) = x5(j) + (x5(j+1)-x5(j)) * factor
	1112	else
	1113	y1(i) = exp( log(x1(j)) + log(x1(j+1)/x1(j)) * factor )
	1114	y2(i) = exp( log(x2(j)) + log(x2(j+1)/x2(j)) * factor )
	1115	y3(i) = exp( log(x3(j)) + log(x3(j+1)/x3(j)) * factor )
	1116	y4(i) = exp( log(x4(j)) + log(x4(j+1)/x4(j)) * factor )
	1117	y5(i) = exp( log(x5(j)) + log(x5(j+1)/x5(j)) * factor )
	1118	end if
	1119
	1120	1 continue
	1121
	1122	return
	1123	end
	1124
	1125
	1126
	1127	c * interhuntdp3veces.f *
	1128
	1129	c ***********************************************************************
	1130
	1131	subroutine interhuntdp3veces ( y1,y2,y3, zz,m,
	1132	@ x1,x2,x3, z,n, opt)
	1133
	1134	c interpolation soubroutine basada en Numerical Recipes HUNT.FOR
	1135	c input values: x(n) at z(n)
	1136	c output values: y(m) at zz(m)
	1137	c options: opt = 1 -> lineal
	1138	c opt=/=1 -> logarithmic
	1139	c ***********************************************************************
	1140	! Arguments
	1141	integer n,m,opt ! I
	1142	real*8 zz(m),z(n) ! I
	1143	real*8 y1(m),y2(m),y3(m) ! O
	1144	real*8 x1(n),x2(n),x3(n) ! I
	1145
	1146
	1147	! Local variables
	1148	integer i, j
	1149	real*8 factor
	1150	real*8 zaux
[2048]	1151	character(len=100) text
[1310]	1152
	1153	!!!!
	1154
	1155	j = 1 ! initial first guess (=n/2 is anothr pssblty)
	1156
	1157	do 1,i=1,m !
	1158
	1159	! Busca indice j donde ocurre q zz(i) esta entre [z(j),z(j+1)]
	1160	zaux = zz(i)
	1161	if (abs(zaux-z(1)).le.0.01d0) then
	1162	zaux=z(1)
	1163	elseif (abs(zaux-z(n)).le.0.01d0) then
	1164	zaux=z(n)
	1165	endif
	1166	call huntdp ( z,n, zaux, j )
	1167	if ( j.eq.0 .or. j.eq.n ) then
	1168	write (,) ' HUNT/ Limits input grid:', z(1),z(n)
	1169	write (,) ' HUNT/ location in new grid:', zz(i)
[2048]	1170	text='INTERHUNTDP3/ Interpolat error. zz out of limits.'
	1171	CALL abort_physic("nlte_aux",text, 1)
[1310]	1172	endif
	1173
	1174	! Perform interpolation
	1175	factor = (zz(i)-z(j))/(z(j+1)-z(j))
	1176	if (opt.eq.1) then
	1177	y1(i) = x1(j) + (x1(j+1)-x1(j)) * factor
	1178	y2(i) = x2(j) + (x2(j+1)-x2(j)) * factor
	1179	y3(i) = x3(j) + (x3(j+1)-x3(j)) * factor
	1180	else
	1181	y1(i) = dexp( dlog(x1(j)) + dlog(x1(j+1)/x1(j)) * factor )
	1182	y2(i) = dexp( dlog(x2(j)) + dlog(x2(j+1)/x2(j)) * factor )
	1183	y3(i) = dexp( dlog(x3(j)) + dlog(x3(j+1)/x3(j)) * factor )
	1184	end if
	1185
	1186	1 continue
	1187
	1188	return
	1189	end
	1190
	1191
	1192	c * interhuntdp4veces.f *
	1193
	1194	c ***********************************************************************
	1195
	1196	subroutine interhuntdp4veces ( y1,y2,y3,y4, zz,m,
	1197	@ x1,x2,x3,x4, z,n, opt)
	1198
	1199	c interpolation soubroutine basada en Numerical Recipes HUNT.FOR
	1200	c input values: x1(n),x2(n),x3(n),x4(n) at z(n)
	1201	c output values: y1(m),y2(m),y3(m),y4(m) at zz(m)
	1202	c options: 1 -> lineal
	1203	c 2 -> logarithmic
	1204	c ***********************************************************************
	1205
	1206	implicit none
	1207
	1208	! Arguments
	1209	integer n,m,opt ! I
	1210	real*8 zz(m),z(n) ! I
	1211	real*8 y1(m),y2(m),y3(m),y4(m) ! O
	1212	real*8 x1(n),x2(n),x3(n),x4(n) ! I
	1213
	1214
	1215	! Local variables
	1216	integer i, j
	1217	real*8 factor
	1218	real*8 zaux
[2048]	1219	character(len=100) text
[1310]	1220
	1221	!!!!
	1222
	1223	j = 1 ! initial first guess (=n/2 is anothr pssblty)
	1224
	1225	do 1,i=1,m !
	1226
	1227	! Caza del indice j donde ocurre que zz(i) esta entre [z(j),z(j+1)]
	1228	zaux = zz(i)
	1229	if (abs(zaux-z(1)).le.0.01d0) then
	1230	zaux=z(1)
	1231	elseif (abs(zaux-z(n)).le.0.01d0) then
	1232	zaux=z(n)
	1233	endif
	1234	call huntdp ( z,n, zaux, j )
	1235	if ( j.eq.0 .or. j.eq.n ) then
	1236	write (,) ' HUNT/ Limits input grid:', z(1),z(n)
	1237	write (,) ' HUNT/ location in new grid:', zz(i)
[2048]	1238	text='INTERHUNTDP4/ Interpolat error. zz out of limits.'
	1239	CALL abort_physic("nlte_aux",text, 1)
[1310]	1240	endif
	1241
	1242	! Perform interpolation
	1243	factor = (zz(i)-z(j))/(z(j+1)-z(j))
	1244	if (opt.eq.1) then
	1245	y1(i) = x1(j) + (x1(j+1)-x1(j)) * factor
	1246	y2(i) = x2(j) + (x2(j+1)-x2(j)) * factor
	1247	y3(i) = x3(j) + (x3(j+1)-x3(j)) * factor
	1248	y4(i) = x4(j) + (x4(j+1)-x4(j)) * factor
	1249	else
	1250	y1(i) = dexp( dlog(x1(j)) + dlog(x1(j+1)/x1(j)) * factor )
	1251	y2(i) = dexp( dlog(x2(j)) + dlog(x2(j+1)/x2(j)) * factor )
	1252	y3(i) = dexp( dlog(x3(j)) + dlog(x3(j+1)/x3(j)) * factor )
	1253	y4(i) = dexp( dlog(x4(j)) + dlog(x4(j+1)/x4(j)) * factor )
	1254	end if
	1255
	1256	1 continue
	1257
	1258	return
	1259	end
	1260
	1261
	1262	c * interhuntdp.f *
	1263
	1264	c ***********************************************************************
	1265
	1266	subroutine interhuntdp ( y1, zz,m,
	1267	@ x1, z,n, opt)
	1268
	1269	c interpolation soubroutine basada en Numerical Recipes HUNT.FOR
	1270	c input values: x1(n) at z(n)
	1271	c output values: y1(m) at zz(m)
	1272	c options: 1 -> lineal
	1273	c 2 -> logarithmic
	1274	c ***********************************************************************
	1275
	1276	implicit none
	1277
	1278	! Arguments
	1279	integer n,m,opt ! I
	1280	real*8 zz(m),z(n) ! I
	1281	real*8 y1(m) ! O
	1282	real*8 x1(n) ! I
	1283
	1284
	1285	! Local variables
	1286	integer i, j
	1287	real*8 factor
	1288	real*8 zaux
[2048]	1289	character(len=100) text
[1310]	1290
	1291	!!!!
	1292
	1293	j = 1 ! initial first guess (=n/2 is anothr pssblty)
	1294
	1295	do 1,i=1,m !
	1296
	1297	! Caza del indice j donde ocurre que zz(i) esta entre [z(j),z(j+1)]
	1298	zaux = zz(i)
	1299	if (abs(zaux-z(1)).le.0.01d0) then
	1300	zaux=z(1)
	1301	elseif (abs(zaux-z(n)).le.0.01d0) then
	1302	zaux=z(n)
	1303	endif
	1304	call huntdp ( z,n, zaux, j )
	1305	if ( j.eq.0 .or. j.eq.n ) then
	1306	write (,) ' HUNT/ Limits input grid:', z(1),z(n)
	1307	write (,) ' HUNT/ location in new grid:', zz(i)
[2048]	1308	text='INTERHUNT/ Interpolat error. zz out of limits.'
	1309	CALL abort_physic("nlte_aux",text, 1)
[1310]	1310	endif
	1311
	1312	! Perform interpolation
	1313	factor = (zz(i)-z(j))/(z(j+1)-z(j))
	1314	if (opt.eq.1) then
	1315	y1(i) = x1(j) + (x1(j+1)-x1(j)) * factor
	1316	else
	1317	y1(i) = dexp( dlog(x1(j)) + dlog(x1(j+1)/x1(j)) * factor )
	1318	end if
	1319
	1320	1 continue
	1321
	1322	return
	1323	end
	1324
	1325
	1326	c * interhunt.f *
	1327
	1328	c***********************************************************************
	1329
	1330	subroutine interhunt ( y1, zz,m,
	1331	@ x1, z,n, opt)
	1332
	1333	c interpolation soubroutine basada en Numerical Recipes HUNT.FOR
	1334	c input values: x1(n) at z(n)
	1335	c output values: y1(m) at zz(m)
	1336	c options: 1 -> lineal
	1337	c 2 -> logarithmic
	1338	c***********************************************************************
	1339
	1340	implicit none
	1341
	1342	! Arguments
	1343	integer n,m,opt ! I
	1344	real zz(m),z(n) ! I
	1345	real y1(m) ! O
	1346	real x1(n) ! I
	1347
	1348
	1349	! Local variables
	1350	integer i, j
	1351	real factor
	1352	real zaux
[2048]	1353	character(len=100) text
[1310]	1354
	1355	!!!!
	1356
	1357	j = 1 ! initial first guess (=n/2 is anothr pssblty)
	1358
	1359	do 1,i=1,m !
	1360
	1361	! Busca indice j donde ocurre q zz(i) esta entre [z(j),z(j+1)]
	1362	zaux = zz(i)
	1363	if (abs(zaux-z(1)).le.0.01) then
	1364	zaux=z(1)
	1365	elseif (abs(zaux-z(n)).le.0.01) then
	1366	zaux=z(n)
	1367	endif
	1368	call hunt ( z,n, zaux, j )
	1369	if ( j.eq.0 .or. j.eq.n ) then
	1370	write (,) ' HUNT/ Limits input grid:', z(1),z(n)
	1371	write (,) ' HUNT/ location in new grid:', zz(i)
[2048]	1372	text='interhunt/ Interpolat error. z out of limits.'
	1373	CALL abort_physic("nlte_aux",text, 1)
[1310]	1374	endif
	1375
	1376	! Perform interpolation
	1377	factor = (zz(i)-z(j))/(z(j+1)-z(j))
	1378	if (opt.eq.1) then
	1379	y1(i) = x1(j) + (x1(j+1)-x1(j)) * factor
	1380	else
	1381	y1(i) = exp( log(x1(j)) + log(x1(j+1)/x1(j)) * factor )
	1382	end if
	1383
	1384
	1385	1 continue
	1386
	1387	return
	1388	end
	1389
	1390
	1391	c * interhuntlimits2veces.f *
	1392
	1393	c***********************************************************************
	1394
	1395	subroutine interhuntlimits2veces
	1396	@ ( y1,y2, zz,m, limite1,limite2,
	1397	@ x1,x2, z,n, opt)
	1398
	1399	c Interpolation soubroutine basada en Numerical Recipes HUNT.FOR
	1400	c Input values: x1,x2(n) at z(n)
	1401	c Output values:
	1402	c y1,y2(m) at zz(m) pero solo entre los indices de zz
	1403	c siguientes: [limite1,limite2]
	1404	c Options: 1 -> linear in z and linear in x
	1405	c 2 -> linear in z and logarithmic in x
	1406	c 3 -> logarithmic in z and linear in x
	1407	c 4 -> logarithmic in z and logaritmic in x
	1408	c
	1409	c NOTAS: Esta subrutina extiende y generaliza la usual
	1410	c "interhunt5veces" en 2 direcciones:
	1411	c - la condicion en los limites es que zz(limite1:limite2)
	1412	c esté dentro de los limites de z (pero quizas no todo zz)
	1413	c - se han añadido 3 opciones mas al caso de interpolacion
	1414	c logaritmica, ahora se hace en log de z, de x o de ambos.
	1415	c Notese que esta subrutina engloba a la interhunt5veces
	1416	c ( esta es reproducible haciendo limite1=1 y limite2=m
	1417	c y usando una de las 2 primeras opciones opt=1,2 )
	1418	c***********************************************************************
	1419
	1420	implicit none
	1421
	1422	! Arguments
	1423	integer n,m,opt, limite1,limite2 ! I
	1424	real zz(m),z(n) ! I
	1425	real y1(m),y2(m) ! O
	1426	real x1(n),x2(n) ! I
	1427
	1428
	1429	! Local variables
	1430	integer i, j
	1431	real factor
	1432	real zaux
[2048]	1433	character(len=100) text
[1310]	1434
	1435	!!!!
	1436
	1437	j = 1 ! initial first guess (=n/2 is anothr pssblty)
	1438
	1439	do 1,i=limite1,limite2
	1440
	1441	! Busca indice j donde ocurre q zz(i) esta entre [z(j),z(j+1)]
	1442	zaux = zz(i)
	1443	if (abs(zaux-z(1)).le.0.01) then
	1444	zaux=z(1)
	1445	elseif (abs(zaux-z(n)).le.0.01) then
	1446	zaux=z(n)
	1447	endif
	1448	call hunt ( z,n, zaux, j )
	1449	if ( j.eq.0 .or. j.eq.n ) then
	1450	write (,) ' HUNT/ Limits input grid:', z(1),z(n)
	1451	write (,) ' HUNT/ location in new grid:', zz(i)
[2048]	1452	text='interhuntlimits/ Interpolat error. z out of limits.'
	1453	CALL abort_physic("nlte_aux",text, 1)
[1310]	1454	endif
	1455
	1456	! Perform interpolation
	1457	if (opt.eq.1) then
	1458	factor = (zz(i)-z(j))/(z(j+1)-z(j))
	1459	y1(i) = x1(j) + (x1(j+1)-x1(j)) * factor
	1460	y2(i) = x2(j) + (x2(j+1)-x2(j)) * factor
	1461
	1462	elseif (opt.eq.2) then
	1463	factor = (zz(i)-z(j))/(z(j+1)-z(j))
	1464	y1(i) = exp( log(x1(j)) + log(x1(j+1)/x1(j)) * factor )
	1465	y2(i) = exp( log(x2(j)) + log(x2(j+1)/x2(j)) * factor )
	1466	elseif (opt.eq.3) then
	1467	factor = (log(zz(i))-log(z(j)))/(log(z(j+1))-log(z(j)))
	1468	y1(i) = x1(j) + (x1(j+1)-x1(j)) * factor
	1469	y2(i) = x2(j) + (x2(j+1)-x2(j)) * factor
	1470	elseif (opt.eq.4) then
	1471	factor = (log(zz(i))-log(z(j)))/(log(z(j+1))-log(z(j)))
	1472	y1(i) = exp( log(x1(j)) + log(x1(j+1)/x1(j)) * factor )
	1473	y2(i) = exp( log(x2(j)) + log(x2(j+1)/x2(j)) * factor )
	1474	end if
	1475
	1476
	1477	1 continue
	1478
	1479	return
	1480	end
	1481
	1482
	1483	c * interhuntlimits5veces.f *
	1484
	1485	c***********************************************************************
	1486
	1487	subroutine interhuntlimits5veces
	1488	@ ( y1,y2,y3,y4,y5, zz,m, limite1,limite2,
	1489	@ x1,x2,x3,x4,x5, z,n, opt)
	1490
	1491	c Interpolation soubroutine basada en Numerical Recipes HUNT.FOR
	1492	c Input values: x1,x2,..,x5(n) at z(n)
	1493	c Output values:
	1494	c y1,y2,...,y5(m) at zz(m) pero solo entre los indices de zz
	1495	c siguientes: [limite1,limite2]
	1496	c Options: 1 -> linear in z and linear in x
	1497	c 2 -> linear in z and logarithmic in x
	1498	c 3 -> logarithmic in z and linear in x
	1499	c 4 -> logarithmic in z and logaritmic in x
	1500	c
	1501	c NOTAS: Esta subrutina extiende y generaliza la usual
	1502	c "interhunt5veces" en 2 direcciones:
	1503	c - la condicion en los limites es que zz(limite1:limite2)
	1504	c esté dentro de los limites de z (pero quizas no todo zz)
	1505	c - se han añadido 3 opciones mas al caso de interpolacion
	1506	c logaritmica, ahora se hace en log de z, de x o de ambos.
	1507	c Notese que esta subrutina engloba a la interhunt5veces
	1508	c ( esta es reproducible haciendo limite1=1 y limite2=m
	1509	c y usando una de las 2 primeras opciones opt=1,2 )
	1510	c***********************************************************************
	1511
	1512	implicit none
	1513
	1514	! Arguments
	1515	integer n,m,opt, limite1,limite2 ! I
	1516	real zz(m),z(n) ! I
	1517	real y1(m),y2(m),y3(m),y4(m),y5(m) ! O
	1518	real x1(n),x2(n),x3(n),x4(n),x5(n) ! I
	1519
	1520
	1521	! Local variables
	1522	integer i, j
	1523	real factor
	1524	real zaux
[2048]	1525	character(len=100) text
[1310]	1526
	1527	!!!!
	1528
	1529	j = 1 ! initial first guess (=n/2 is anothr pssblty)
	1530
	1531	do 1,i=limite1,limite2
	1532
	1533	! Busca indice j donde ocurre q zz(i) esta entre [z(j),z(j+1)]
	1534	zaux = zz(i)
	1535	if (abs(zaux-z(1)).le.0.01) then
	1536	zaux=z(1)
	1537	elseif (abs(zaux-z(n)).le.0.01) then
	1538	zaux=z(n)
	1539	endif
	1540	call hunt ( z,n, zaux, j )
	1541	if ( j.eq.0 .or. j.eq.n ) then
	1542	write (,) ' HUNT/ Limits input grid:', z(1),z(n)
	1543	write (,) ' HUNT/ location in new grid:', zz(i)
[2048]	1544	text='interhuntlimits/ Interpolat error. z out of limits.'
	1545	CALL abort_physic("nlte_aux",text, 1)
[1310]	1546	endif
	1547
	1548	! Perform interpolation
	1549	if (opt.eq.1) then
	1550	factor = (zz(i)-z(j))/(z(j+1)-z(j))
	1551	y1(i) = x1(j) + (x1(j+1)-x1(j)) * factor
	1552	y2(i) = x2(j) + (x2(j+1)-x2(j)) * factor
	1553	y3(i) = x3(j) + (x3(j+1)-x3(j)) * factor
	1554	y4(i) = x4(j) + (x4(j+1)-x4(j)) * factor
	1555	y5(i) = x5(j) + (x5(j+1)-x5(j)) * factor
	1556	elseif (opt.eq.2) then
	1557	factor = (zz(i)-z(j))/(z(j+1)-z(j))
	1558	y1(i) = exp( log(x1(j)) + log(x1(j+1)/x1(j)) * factor )
	1559	y2(i) = exp( log(x2(j)) + log(x2(j+1)/x2(j)) * factor )
	1560	y3(i) = exp( log(x3(j)) + log(x3(j+1)/x3(j)) * factor )
	1561	y4(i) = exp( log(x4(j)) + log(x4(j+1)/x4(j)) * factor )
	1562	y5(i) = exp( log(x5(j)) + log(x5(j+1)/x5(j)) * factor )
	1563	elseif (opt.eq.3) then
	1564	factor = (log(zz(i))-log(z(j)))/(log(z(j+1))-log(z(j)))
	1565	y1(i) = x1(j) + (x1(j+1)-x1(j)) * factor
	1566	y2(i) = x2(j) + (x2(j+1)-x2(j)) * factor
	1567	y3(i) = x3(j) + (x3(j+1)-x3(j)) * factor
	1568	y4(i) = x4(j) + (x4(j+1)-x4(j)) * factor
	1569	y5(i) = x5(j) + (x5(j+1)-x5(j)) * factor
	1570	elseif (opt.eq.4) then
	1571	factor = (log(zz(i))-log(z(j)))/(log(z(j+1))-log(z(j)))
	1572	y1(i) = exp( log(x1(j)) + log(x1(j+1)/x1(j)) * factor )
	1573	y2(i) = exp( log(x2(j)) + log(x2(j+1)/x2(j)) * factor )
	1574	y3(i) = exp( log(x3(j)) + log(x3(j+1)/x3(j)) * factor )
	1575	y4(i) = exp( log(x4(j)) + log(x4(j+1)/x4(j)) * factor )
	1576	y5(i) = exp( log(x5(j)) + log(x5(j+1)/x5(j)) * factor )
	1577	end if
	1578
	1579
	1580	1 continue
	1581
	1582	return
	1583	end
	1584
	1585
	1586
	1587	c * interhuntlimits.f *
	1588
	1589	c***********************************************************************
	1590
	1591	subroutine interhuntlimits ( y, zz,m, limite1,limite2,
	1592	@ x, z,n, opt)
	1593
	1594	c Interpolation soubroutine basada en Numerical Recipes HUNT.FOR
	1595	c Input values: x(n) at z(n)
	1596	c Output values: y(m) at zz(m) pero solo entre los indices de zz
	1597	c siguientes: [limite1,limite2]
	1598	c Options: 1 -> linear in z and linear in x
	1599	c 2 -> linear in z and logarithmic in x
	1600	c 3 -> logarithmic in z and linear in x
	1601	c 4 -> logarithmic in z and logaritmic in x
	1602	c
	1603	c NOTAS: Esta subrutina extiende y generaliza la usual "interhunt"
	1604	c en 2 direcciones:
	1605	c - la condicion en los limites es que zz(limite1:limite2)
	1606	c esté dentro de los limites de z (pero quizas no todo zz)
	1607	c - se han añadido 3 opciones mas al caso de interpolacion
	1608	c logaritmica, ahora se hace en log de z, de x o de ambos.
	1609	c Notese que esta subrutina engloba a la usual interhunt
	1610	c ( esta es reproducible haciendo limite1=1 y limite2=m
	1611	c y usando una de las 2 primeras opciones opt=1,2 )
	1612	c***********************************************************************
	1613
	1614	implicit none
	1615
	1616	! Arguments
	1617	integer n,m,opt, limite1,limite2 ! I
	1618	real zz(m),z(n) ! I
	1619	real y(m) ! O
	1620	real x(n) ! I
	1621
	1622
	1623	! Local variables
	1624	integer i, j
	1625	real factor
	1626	real zaux
[2048]	1627	character(len=100) text
[1310]	1628
	1629	!!!!
	1630
	1631	j = 1 ! initial first guess (=n/2 is anothr pssblty)
	1632
	1633	do 1,i=limite1,limite2
	1634
	1635	! Busca indice j donde ocurre q zz(i) esta entre [z(j),z(j+1)]
	1636	zaux = zz(i)
	1637	if (abs(zaux-z(1)).le.0.01) then
	1638	zaux=z(1)
	1639	elseif (abs(zaux-z(n)).le.0.01) then
	1640	zaux=z(n)
	1641	endif
	1642	call hunt ( z,n, zaux, j )
	1643	if ( j.eq.0 .or. j.eq.n ) then
	1644	write (,) ' HUNT/ Limits input grid:', z(1),z(n)
	1645	write (,) ' HUNT/ location in new grid:', zz(i)
[2048]	1646	text='interhuntlimits/ Interpolat error. z out of limits.'
	1647	CALL abort_physic("nlte_aux",text, 1)
[1310]	1648	endif
	1649
	1650	! Perform interpolation
	1651	if (opt.eq.1) then
	1652	factor = (zz(i)-z(j))/(z(j+1)-z(j))
	1653	y(i) = x(j) + (x(j+1)-x(j)) * factor
	1654	elseif (opt.eq.2) then
	1655	factor = (zz(i)-z(j))/(z(j+1)-z(j))
	1656	y(i) = exp( log(x(j)) + log(x(j+1)/x(j)) * factor )
	1657	elseif (opt.eq.3) then
	1658	factor = (log(zz(i))-log(z(j)))/(log(z(j+1))-log(z(j)))
	1659	y(i) = x(j) + (x(j+1)-x(j)) * factor
	1660	elseif (opt.eq.4) then
	1661	factor = (log(zz(i))-log(z(j)))/(log(z(j+1))-log(z(j)))
	1662	y(i) = exp( log(x(j)) + log(x(j+1)/x(j)) * factor )
	1663	end if
	1664
	1665
	1666	1 continue
	1667
	1668	return
	1669	end
	1670
	1671
	1672	c * lubksb_dp.f *
	1673
	1674	subroutine lubksb_dp(a,n,np,indx,b)
	1675
	1676	implicit none
	1677
	1678	integer,intent(in) :: n,np
	1679	real*8,intent(in) :: a(np,np)
	1680	integer,intent(in) :: indx(n)
	1681	real*8,intent(out) :: b(n)
	1682
	1683	real*8 sum
	1684	integer ii, ll, i, j
	1685
	1686	ii=0
	1687	do 12 i=1,n
	1688	ll=indx(i)
	1689	sum=b(ll)
	1690	b(ll)=b(i)
	1691	if (ii.ne.0)then
	1692	do 11 j=ii,i-1
	1693	sum=sum-a(i,j)*b(j)
	1694	11 continue
	1695	else if (sum.ne.0.0) then
	1696	ii=i
	1697	endif
	1698	b(i)=sum
	1699	12 continue
	1700	do 14 i=n,1,-1
	1701	sum=b(i)
	1702	if(i.lt.n)then
	1703	do 13 j=i+1,n
	1704	sum=sum-a(i,j)*b(j)
	1705	13 continue
	1706	endif
	1707	b(i)=sum/a(i,i)
	1708	14 continue
	1709	return
	1710	end
	1711
	1712
	1713	c * ludcmp_dp.f *
	1714
	1715	subroutine ludcmp_dp(a,n,np,indx,d)
	1716
	1717	implicit none
	1718
	1719	integer,intent(in) :: n, np
	1720	real*8,intent(inout) :: a(np,np)
	1721	real*8,intent(out) :: d
	1722	integer,intent(out) :: indx(n)
	1723
	1724	integer nmax, i, j, k, imax
	1725	real*8 tiny
	1726	parameter (nmax=100,tiny=1.0d-20)
	1727	real*8 vv(nmax), aamax, sum, dum
	1728
	1729
	1730	d=1.0d0
	1731	do 12 i=1,n
	1732	aamax=0.0d0
	1733	do 11 j=1,n
	1734	if (abs(a(i,j)).gt.aamax) aamax=abs(a(i,j))
	1735	11 continue
	1736	if (aamax.eq.0.0) then
	1737	write(,) 'ludcmp_dp: singular matrix!'
	1738	stop
	1739	endif
	1740	vv(i)=1.0d0/aamax
	1741	12 continue
	1742	do 19 j=1,n
	1743	if (j.gt.1) then
	1744	do 14 i=1,j-1
	1745	sum=a(i,j)
	1746	if (i.gt.1)then
	1747	do 13 k=1,i-1
	1748	sum=sum-a(i,k)*a(k,j)
	1749	13 continue
	1750	a(i,j)=sum
	1751	endif
	1752	14 continue
	1753	endif
	1754	aamax=0.0d0
	1755	do 16 i=j,n
	1756	sum=a(i,j)
	1757	if (j.gt.1)then
	1758	do 15 k=1,j-1
	1759	sum=sum-a(i,k)*a(k,j)
	1760	15 continue
	1761	a(i,j)=sum
	1762	endif
	1763	dum=vv(i)*abs(sum)
	1764	if (dum.ge.aamax) then
	1765	imax=i
	1766	aamax=dum
	1767	endif
	1768	16 continue
	1769	if (j.ne.imax)then
	1770	do 17 k=1,n
	1771	dum=a(imax,k)
	1772	a(imax,k)=a(j,k)
	1773	a(j,k)=dum
	1774	17 continue
	1775	d=-d
	1776	vv(imax)=vv(j)
	1777	endif
	1778	indx(j)=imax
	1779	if(j.ne.n)then
	1780	if(a(j,j).eq.0.0)a(j,j)=tiny
	1781	dum=1.0d0/a(j,j)
	1782	do 18 i=j+1,n
	1783	a(i,j)=a(i,j)*dum
	1784	18 continue
	1785	endif
	1786	19 continue
	1787	if(a(n,n).eq.0.0)a(n,n)=tiny
	1788	return
	1789	end
	1790
	1791
	1792	c * LUdec.f *
	1793
	1794	ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	1795	c
	1796	c Solution of linear equation without inverting matrix
	1797	c using LU decomposition:
	1798	c AA * xx = bb AA, bb: known
	1799	c xx: to be found
	1800	c AA and bb are not modified in this subroutine
	1801	c
	1802	c MALV , Sep 2007
	1803	ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	1804
	1805	subroutine LUdec(xx,aa,bb,m,n)
	1806
	1807	implicit none
	1808
	1809	! Arguments
	1810	integer,intent(in) :: m, n
	1811	real*8,intent(in) :: aa(m,m), bb(m)
	1812	real*8,intent(out) :: xx(m)
	1813
	1814
	1815	! Local variables
	1816	real*8 a(n,n), b(n), x(n), d
	1817	integer i, j, indx(n)
	1818
	1819
	1820	! Subrutinas utilizadas
	1821	! ludcmp_dp, lubksb_dp
	1822
	1823	!!!!!!!!!!!!!!!Comienza el programa !!!!!!!!!!!!!!
	1824
	1825	do i=1,n
	1826	b(i) = bb(i+1)
	1827	do j=1,n
	1828	a(i,j) = aa(i+1,j+1)
	1829	enddo
	1830	enddo
	1831
	1832	! Descomposicion de auxm1
	1833	call ludcmp_dp ( a, n, n, indx, d)
	1834
	1835	! Sustituciones foward y backwards para hallar la solucion
	1836	do i=1,n
	1837	x(i) = b(i)
	1838	enddo
	1839	call lubksb_dp( a, n, n, indx, x )
	1840
	1841	do i=1,n
	1842	xx(i+1) = x(i)
	1843	enddo
	1844
	1845	return
	1846	end
	1847
	1848
	1849	c * mat_oper.f *
	1850
	1851	ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	1852
	1853	c ***********************************************************************
	1854	subroutine unit(a,n)
	1855	c store the unit value in the diagonal of a
	1856	c ***********************************************************************
	1857	implicit none
	1858	real*8 a(n,n)
	1859	integer n,i,j,k
	1860	do 1,i=2,n-1
	1861	do 2,j=2,n-1
	1862	if(i.eq.j) then
	1863	a(i,j) = 1.d0
	1864	else
	1865	a(i,j)=0.0d0
	1866	end if
	1867	2 continue
	1868	1 continue
	1869	do k=1,n
	1870	a(n,k) = 0.0d0
	1871	a(1,k) = 0.0d0
	1872	a(k,1) = 0.0d0
	1873	a(k,n) = 0.0d0
	1874	end do
	1875	return
	1876	end
	1877
	1878	c ***********************************************************************
	1879	subroutine diago(a,v,n)
	1880	c store the vector v in the diagonal elements of the square matrix a
	1881	c ***********************************************************************
	1882	implicit none
	1883
	1884	integer n,i,j,k
	1885	real*8 a(n,n),v(n)
	1886
	1887	do 1,i=2,n-1
	1888	do 2,j=2,n-1
	1889	if(i.eq.j) then
	1890	a(i,j) = v(i)
	1891	else
	1892	a(i,j)=0.0d0
	1893	end if
	1894	2 continue
	1895	1 continue
	1896	do k=1,n
	1897	a(n,k) = 0.0d0
	1898	a(1,k) = 0.0d0
	1899	a(k,1) = 0.0d0
	1900	a(k,n) = 0.0d0
	1901	end do
	1902	return
	1903	end
	1904
	1905	c ***********************************************************************
	1906	subroutine invdiag(a,b,n)
	1907	c inverse of a diagonal matrix
	1908	c ***********************************************************************
	1909	implicit none
	1910
	1911	integer n,i,j,k
	1912	real*8 a(n,n),b(n,n)
	1913
	1914	do 1,i=2,n-1
	1915	do 2,j=2,n-1
	1916	if (i.eq.j) then
	1917	a(i,j) = 1.d0/b(i,i)
	1918	else
	1919	a(i,j)=0.0d0
	1920	end if
	1921	2 continue
	1922	1 continue
	1923	do k=1,n
	1924	a(n,k) = 0.0d0
	1925	a(1,k) = 0.0d0
	1926	a(k,1) = 0.0d0
	1927	a(k,n) = 0.0d0
	1928	end do
	1929	return
	1930	end
	1931
	1932
	1933	c ***********************************************************************
	1934	subroutine samem (a,m,n)
	1935	c store the matrix m in the matrix a
	1936	c ***********************************************************************
	1937	implicit none
	1938	real*8 a(n,n),m(n,n)
	1939	integer n,i,j,k
	1940	do 1,i=2,n-1
	1941	do 2,j=2,n-1
	1942	a(i,j) = m(i,j)
	1943	2 continue
	1944	1 continue
	1945	do k=1,n
	1946	a(n,k) = 0.0d0
	1947	a(1,k) = 0.0d0
	1948	a(k,1) = 0.0d0
	1949	a(k,n) = 0.0d0
	1950	end do
	1951	return
	1952	end
	1953
	1954
	1955	c ***********************************************************************
	1956	subroutine mulmv(a,b,c,n)
	1957	c do a(i)=b(i,j)*c(j). a, b, and c must be distint
	1958	c ***********************************************************************
	1959	implicit none
	1960
	1961	integer n,i,j
	1962	real*8 a(n),b(n,n),c(n),sum
	1963
	1964	do 1,i=2,n-1
	1965	sum=0.0d0
	1966	do 2,j=2,n-1
	1967	sum = sum + b(i,j) * c(j)
	1968	2 continue
	1969	a(i)=sum
	1970	1 continue
	1971	a(1) = 0.0d0
	1972	a(n) = 0.0d0
	1973	return
	1974	end
	1975
	1976
	1977	c ***********************************************************************
	1978	subroutine trucodiag(a,b,c,d,e,n)
	1979	c inputs: matrices b,c,d,e
	1980	c output: matriz diagonal a
	1981	c Operacion a realizar: a = b * c^(-1) * d + e
	1982	c La matriz c va a ser invertida
	1983	c Todas las matrices de entrada son diagonales excepto b
	1984	c Aprovechamos esa condicion para invertir c, acelerar el calculo, y
	1985	c ademas, para forzar que a sea diagonal
	1986	c ***********************************************************************
	1987	implicit none
	1988	real*8 a(n,n),b(n,n),c(n,n),d(n,n),e(n,n), sum
	1989	integer n,i,j,k
	1990	do 1,i=2,n-1
	1991	sum=0.0d0
	1992	do 2,j=2,n-1
	1993	sum=sum+ b(i,j) * d(j,j)/c(j,j)
	1994	2 continue
	1995	a(i,i) = sum + e(i,i)
	1996	1 continue
	1997	do k=1,n
	1998	a(n,k) = 0.0d0
	1999	a(1,k) = 0.0d0
	2000	a(k,1) = 0.0d0
	2001	a(k,n) = 0.0d0
	2002	end do
	2003	return
	2004	end
	2005
	2006
	2007	c ***********************************************************************
	2008	subroutine trucommvv(v,b,c,u,w,n)
	2009	c inputs: matrices b,c , vectores u,w
	2010	c output: vector v
	2011	c Operacion a realizar: v = b * c^(-1) * u + w
	2012	c La matriz c va a ser invertida
	2013	c c es diagonal, b no
	2014	c Aprovechamos esa condicion para invertir c, y acelerar el calculo
	2015	c ***********************************************************************
	2016	implicit none
	2017	real*8 v(n),b(n,n),c(n,n),u(n),w(n), sum
	2018	integer n,i,j
	2019	do 1,i=2,n-1
	2020	sum=0.0d0
	2021	do 2,j=2,n-1
	2022	sum=sum+ b(i,j) * u(j)/c(j,j)
	2023	2 continue
	2024	v(i) = sum + w(i)
	2025	1 continue
	2026	v(1) = 0.d0
	2027	v(n) = 0.d0
	2028	return
	2029	end
	2030
	2031
	2032	c ***********************************************************************
	2033	subroutine sypvmv(v,u,c,w,n)
	2034	c inputs: matriz diagonal c , vectores u,w
	2035	c output: vector v
	2036	c Operacion a realizar: v = u + c * w
	2037	c ***********************************************************************
	2038	implicit none
	2039	real*8 v(n),u(n),c(n,n),w(n)
	2040	integer n,i
	2041	do 1,i=2,n-1
	2042	v(i)= u(i) + c(i,i) * w(i)
	2043	1 continue
	2044	v(1) = 0.0d0
	2045	v(n) = 0.0d0
	2046	return
	2047	end
	2048
	2049
	2050	c ***********************************************************************
	2051	subroutine sumvv(a,b,c,n)
	2052	c a(i)=b(i)+c(i)
	2053	c ***********************************************************************
	2054	implicit none
	2055
	2056	integer n,i
	2057	real*8 a(n),b(n),c(n)
	2058
	2059	do 1,i=2,n-1
	2060	a(i)= b(i) + c(i)
	2061	1 continue
	2062	a(1) = 0.0d0
	2063	a(n) = 0.0d0
	2064	return
	2065	end
	2066
	2067
	2068	c ***********************************************************************
	2069	subroutine sypvvv(a,b,c,d,n)
	2070	c a(i)=b(i)+c(i)*d(i)
	2071	c ***********************************************************************
	2072	implicit none
	2073	real*8 a(n),b(n),c(n),d(n)
	2074	integer n,i
	2075	do 1,i=2,n-1
	2076	a(i)= b(i) + c(i) * d(i)
	2077	1 continue
	2078	a(1) = 0.0d0
	2079	a(n) = 0.0d0
	2080	return
	2081	end
	2082
	2083
	2084	c ***********************************************************************
	2085	! subroutine zerom(a,n)
	2086	c a(i,j)= 0.0
	2087	c ***********************************************************************
	2088	! implicit none
	2089	! integer n,i,j
	2090	! real*8 a(n,n)
	2091
	2092	! do 1,i=1,n
	2093	! do 2,j=1,n
	2094	! a(i,j) = 0.0d0
	2095	! 2 continue
	2096	! 1 continue
	2097	! return
	2098	! end
	2099
	2100
	2101	c ***********************************************************************
	2102	subroutine zero4m(a,b,c,d,n)
	2103	c a(i,j) = b(i,j) = c(i,j) = d(i,j) = 0.0
	2104	c ***********************************************************************
	2105	implicit none
	2106	real*8 a(n,n), b(n,n), c(n,n), d(n,n)
	2107	integer n
	2108	a(1:n,1:n)=0.d0
	2109	b(1:n,1:n)=0.d0
	2110	c(1:n,1:n)=0.d0
	2111	d(1:n,1:n)=0.d0
	2112	! do 1,i=1,n
	2113	! do 2,j=1,n
	2114	! a(i,j) = 0.0d0
	2115	! b(i,j) = 0.0d0
	2116	! c(i,j) = 0.0d0
	2117	! d(i,j) = 0.0d0
	2118	! 2 continue
	2119	! 1 continue
	2120	return
	2121	end
	2122
	2123
	2124	c ***********************************************************************
	2125	subroutine zero3m(a,b,c,n)
	2126	c a(i,j) = b(i,j) = c(i,j) = 0.0
	2127	c **********************************************************************
	2128	implicit none
	2129	real*8 a(n,n), b(n,n), c(n,n)
	2130	integer n
	2131	a(1:n,1:n)=0.d0
	2132	b(1:n,1:n)=0.d0
	2133	c(1:n,1:n)=0.d0
	2134	! do 1,i=1,n
	2135	! do 2,j=1,n
	2136	! a(i,j) = 0.0d0
	2137	! b(i,j) = 0.0d0
	2138	! c(i,j) = 0.0d0
	2139	! 2 continue
	2140	! 1 continue
	2141	return
	2142	end
	2143
	2144
	2145	c ***********************************************************************
	2146	subroutine zero2m(a,b,n)
	2147	c a(i,j) = b(i,j) = 0.0
	2148	c ***********************************************************************
	2149	implicit none
	2150	real*8 a(n,n), b(n,n)
	2151	integer n
	2152	a(1:n,1:n)=0.d0
	2153	b(1:n,1:n)=0.d0
	2154	! do 1,i=1,n
	2155	! do 2,j=1,n
	2156	! a(i,j) = 0.0d0
	2157	! b(i,j) = 0.0d0
	2158	! 2 continue
	2159	! 1 continue
	2160	return
	2161	end
	2162
	2163
	2164	c ***********************************************************************
	2165	! subroutine zerov(a,n)
	2166	c a(i)= 0.0
	2167	c ***********************************************************************
	2168	! implicit none
	2169	! real*8 a(n)
	2170	! integer n,i
	2171	! do 1,i=1,n
	2172	! a(i) = 0.0d0
	2173	! 1 continue
	2174	! return
	2175	! end
	2176
	2177
	2178	c ***********************************************************************
	2179	subroutine zero4v(a,b,c,d,n)
	2180	c a(i) = b(i) = c(i) = d(i,j) = 0.0
	2181	c ***********************************************************************
	2182	implicit none
	2183	real*8 a(n), b(n), c(n), d(n)
	2184	integer n
	2185	a(1:n)=0.d0
	2186	b(1:n)=0.d0
	2187	c(1:n)=0.d0
	2188	d(1:n)=0.d0
	2189	! do 1,i=1,n
	2190	! a(i) = 0.0d0
	2191	! b(i) = 0.0d0
	2192	! c(i) = 0.0d0
	2193	! d(i) = 0.0d0
	2194	! 1 continue
	2195	return
	2196	end
	2197
	2198
	2199	c ***********************************************************************
	2200	subroutine zero3v(a,b,c,n)
	2201	c a(i) = b(i) = c(i) = 0.0
	2202	c ***********************************************************************
	2203	implicit none
	2204	real*8 a(n), b(n), c(n)
	2205	integer n
	2206	a(1:n)=0.d0
	2207	b(1:n)=0.d0
	2208	c(1:n)=0.d0
	2209	! do 1,i=1,n
	2210	! a(i) = 0.0d0
	2211	! b(i) = 0.0d0
	2212	! c(i) = 0.0d0
	2213	! 1 continue
	2214	return
	2215	end
	2216
	2217
	2218	c ***********************************************************************
	2219	subroutine zero2v(a,b,n)
	2220	c a(i) = b(i) = 0.0
	2221	c ***********************************************************************
	2222	implicit none
	2223	real*8 a(n), b(n)
	2224	integer n
	2225	a(1:n)=0.d0
	2226	b(1:n)=0.d0
	2227	! do 1,i=1,n
	2228	! a(i) = 0.0d0
	2229	! b(i) = 0.0d0
	2230	! 1 continue
	2231	return
	2232	end
	2233
	2234	c ***********************************************************************
	2235
	2236
	2237	c****************************************************************************
	2238
	2239	c * suaviza.f *
	2240
	2241	c*****************************************************************************
	2242	c
	2243	subroutine suaviza ( x, n, ismooth, y )
	2244	c
	2245	c x - input and return values
	2246	c y - auxiliary vector
	2247	c ismooth = 0 --> no smoothing is performed
	2248	c ismooth = 1 --> weak smoothing (5 points, centred weighted)
	2249	c ismooth = 2 --> normal smoothing (3 points, evenly weighted)
	2250	c ismooth = 3 --> strong smoothing (5 points, evenly weighted)
	2251
	2252
	2253	c august 1991
	2254	c*****************************************************************************
	2255
	2256	implicit none
	2257
	2258	integer n, imax, imin, i, ismooth
	2259	real*8 x(n), y(n)
	2260	c*****************************************************************************
	2261
	2262	imin=1
	2263	imax=n
	2264
	2265	if (ismooth.eq.0) then
	2266
	2267	return
	2268
	2269	elseif (ismooth.eq.1) then ! 5 points, with central weighting
	2270
	2271	do i=imin,imax
	2272	if(i.eq.imin)then
	2273	y(i)=x(imin)
	2274	elseif(i.eq.imax)then
	2275	y(i)=x(imax-1)+(x(imax-1)-x(imax-3))/2.d0
	2276	elseif(i.gt.(imin+1) .and. i.lt.(imax-1) )then
	2277	y(i) = ( x(i+2)/4.d0 + x(i+1)/2.d0 + 2.d0*x(i)/3.d0 +
	2278	@ x(i-1)/2.d0 + x(i-2)/4.d0 )* 6.d0/13.d0
	2279	else
	2280	y(i)=(x(i+1)/2.d0+x(i)+x(i-1)/2.d0)/2.d0
	2281	end if
	2282	end do
	2283
	2284	elseif (ismooth.eq.2) then ! 3 points, evenly spaced
	2285
	2286	do i=imin,imax
	2287	if(i.eq.imin)then
	2288	y(i)=x(imin)
	2289	elseif(i.eq.imax)then
	2290	y(i)=x(imax-1)+(x(imax-1)-x(imax-3))/2.d0
	2291	else
	2292	y(i) = ( x(i+1)+x(i)+x(i-1) )/3.d0
	2293	end if
	2294	end do
	2295
	2296	elseif (ismooth.eq.3) then ! 5 points, evenly spaced
	2297
	2298	do i=imin,imax
	2299	if(i.eq.imin)then
	2300	y(i) = x(imin)
	2301	elseif(i.eq.(imin+1) .or. i.eq.(imax-1))then
	2302	y(i) = ( x(i+1)+x(i)+x(i-1) )/3.d0
	2303	elseif(i.eq.imax)then
	2304	y(i) = ( x(imax-1) + x(imax-1) + x(imax-2) ) / 3.d0
	2305	else
	2306	y(i) = ( x(i+2)+x(i+1)+x(i)+x(i-1)+x(i-2) )/5.d0
	2307	end if
	2308	end do
	2309
	2310	else
	2311
	2312	write (,) ' Error in suaviza.f Wrong ismooth value.'
	2313	stop
	2314
	2315	endif
	2316
	2317	c rehago el cambio, para devolver x(i)
	2318	do i=imin,imax
	2319	x(i)=y(i)
	2320	end do
	2321
	2322	return
	2323	end
	2324
	2325
	2326	c ***********************************************************************
	2327	subroutine mulmmf90(a,b,c,n)
	2328	c ***********************************************************************
	2329	implicit none
	2330	real*8 a(n,n), b(n,n), c(n,n)
	2331	integer n
	2332
	2333	a=matmul(b,c)
	2334	a(1,:)=0.d0
	2335	a(:,1)=0.d0
	2336	a(n,:)=0.d0
	2337	a(:,n)=0.d0
	2338
	2339	return
	2340	end
	2341
	2342
	2343	c ***********************************************************************
	2344	subroutine resmmf90(a,b,c,n)
	2345	c ***********************************************************************
	2346	implicit none
	2347	real*8 a(n,n), b(n,n), c(n,n)
	2348	integer n
	2349
	2350	a=b-c
	2351	a(1,:)=0.d0
	2352	a(:,1)=0.d0
	2353	a(n,:)=0.d0
	2354	a(:,n)=0.d0
	2355
	2356	return
	2357	end
	2358
	2359
	2360	c*******************************************************************
	2361
	2362	subroutine gethist_03 (ihist)
	2363
	2364	c*******************************************************************
	2365
	2366	implicit none
	2367
	2368	#include "nlte_paramdef.h"
	2369	#include "nlte_commons.h"
	2370
	2371
	2372	c arguments
	2373	integer ihist
	2374
	2375	c local variables
	2376	integer j, r
	2377
	2378	c ***************
	2379
	2380	nbox = nbox_stored(ihist)
	2381	do j=1,mm_stored(ihist)
	2382	thist(j) = thist_stored(ihist,j)
	2383	do r=1,nbox_stored(ihist)
	2384	no(r) = no_stored(ihist,r)
	2385	sk1(j,r) = sk1_stored(ihist,j,r)
	2386	xls1(j,r) = xls1_stored(ihist,j,r)
	2387	xld1(j,r) = xld1_stored(ihist,j,r)
	2388	enddo
	2389	enddo
	2390
	2391
	2392	return
	2393	end
	2394
	2395
	2396	c *******************************************************************
	2397
	2398	subroutine rhist_03 (ihist)
	2399
	2400	c *******************************************************************
	2401
	2402	implicit none
	2403
	2404	#include "nlte_paramdef.h"
	2405	#include "nlte_commons.h"
	2406
	2407
	2408	c arguments
	2409	integer ihist
	2410
	2411	c local variables
	2412	integer j, r
	2413	real*8 xx
	2414
	2415	c ***************
	2416
	2417	open(unit=3,file=hisfile,status='old')
	2418
	2419	read(3,*)
	2420	read(3,*)
	2421	read(3,*) mm_stored(ihist)
	2422	read(3,*)
	2423	read(3,*) nbox_stored(ihist)
	2424	read(3,*)
	2425
	2426	if ( nbox_stored(ihist) .gt. nbox_max ) then
	2427	write (,) ' nbox too large in input file ', hisfile
	2428	stop ' Check maximum number nbox_max in mz1d.par '
	2429	endif
	2430
	2431	do j=mm_stored(ihist),1,-1
	2432	read(3,*) thist_stored(ihist,j)
	2433	do r=1,nbox_stored(ihist)
	2434	read(3,*) no_stored(ihist,r),
	2435	& sk1_stored(ihist,j,r),
	2436	& xls1_stored(ihist,j,r),
	2437	& xx,
	2438	& xld1_stored(ihist,j,r)
	2439	enddo
	2440
	2441	enddo
	2442
	2443	close(unit=3)
	2444
	2445
	2446	return
	2447	end

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

Download in other formats:

Original Format