Context Navigation

nlte_aux.F @ 3556

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

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

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

Download in other formats:

Original Format