Context Navigation

nlte_aux.F @ 3026

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

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

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

Download in other formats:

Original Format