source: trunk/LMDZ.MARS/libf/phymars/nlte_aux.F @ 1635

Last change on this file since 1635 was 1124, checked in by emillour, 11 years ago

Mars GCM:

  • Improved 15um cooling routines, with also better handling of errors.

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