Context Navigation

nlte_aux.F @ 1543

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

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

Download in other formats:

Original Format