Context Navigation

nlte_aux.F @ 757

Last change on this file since 757 was 757, checked in by emillour, 12 years ago

Mars GCM:

Improvement of the NLTE 15um scheme (for running with nltemodel = 2); now MUCH faster than previously (by a factor 5 or so):
Improvements included to the parameterization:
- Cool-to-space calculation included above P(atm)=1e-10, with a soft merging to the full result (without the CTS approximation) below that level
- exhaustive cleaning of the code, including FTNCHK and reordering of loops, subroutines and internal calls
- simplification of the precomputed tables of CO2 bands' atmospheric transmittances
- the two internal grids (the one used in the CTS region and the one below) have been , in order to reduce the CPU time consumption
- reading of the spectroscopic histograms is made only once, at the beginning of the GCM, to avoid repetitive readings of ASCII files
- F90 matrix operations (matmul,...) included.
Changes in routines:
- removed nlte_leedat.F
- updated nlte_calc.F, nlte_tcool.F, nlte_aux.F
- updated nlte_commons.h, nlte_paramdef.h
- added nlte_setup.F
Important: The input files (in the NLTEDAT directory) read as input by these routines have changed. the NLTEDAT directory should now on contain only the following files:

deltanu26.dat enelow27.dat hid26-3.dat parametp_Tstar_IAA1204.dat
deltanu27.dat enelow28.dat hid26-4.dat parametp_VC_IAA1204.dat
deltanu28.dat enelow36.dat hid27-1.dat
deltanu36.dat hid26-1.dat hid28-1.dat
enelow26.dat hid26-2.dat hid36-1.dat

FGG+MALV

File size: 62.9 KB

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

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

Download in other formats:

Original Format