Context Navigation

nlte_aux.F @ 3026

Last change on this file since 3026 was 3018, checked in by emillour, 17 months ago
Mars PCM: Further code cleanup with NLTE routines; converted nlte_paramdef.h to module nlte_paramdef_h.F90 and nlte_commons.h to module nlte_commons_h.F90 (could not turn nlte_aux.F, nlte_setup.F and nlte_calc.F into modules due to circular dependencies; would require further code reorganization). EM
File size: 67.1 KB

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

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

Download in other formats:

Original Format