Context Navigation

nlte_aux.F @ 2800

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

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

Download in other formats:

Original Format