Context Navigation

nlte_aux.F @ 1009

Last change on this file since 1009 was 769, checked in by acolaitis, 12 years ago
MESOSCALE. minor modifications to allow compilation with new physics + nesting with latest svn revision
File size: 62.9 KB

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

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

Download in other formats:

Original Format