source: trunk/LMDZ.MARS/libf/phymars/nlte_calc.F @ 1009

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

Mars GCM:

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

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

FGG+MALV
File size: 39.3 KB
Line 
1c**********************************************************************
2c     
3c     Includes the following 1-d model subroutines:
4c     
5c     -MZESC110_dlvr11_03.f
6c     -MZTUD110_dlvr11_03.f
7c     -MZMC121_dlvr11_03.f
8c     -MZTUD121_dlvr11_03.f
9c     -MZESC121_dlvr11_03.f
10c     -MZESC121sub_dlvr11_03.f
11c     -MZTVC121_dlvr11.f
12c     -MZTVC121sub_dlvr11_03.f
13
14
15
16c     *** Old MZESC110_dlvr11_03.f
17
18c**********************************************************************
19
20c***********************************************************************
21      subroutine MZESC110 (nl_cts_real, nzy_cts_real)
22c***********************************************************************
23
24      implicit none
25
26      include 'datafile.h'
27      include 'nlte_paramdef.h'
28      include 'nlte_commons.h'
29
30c     arguments
31      integer     nl_cts_real, nzy_cts_real ! i
32
33c     old arguments
34      integer         ierr      ! o
35      real*8          varerr    ! o
36
37c     local variables and constants
38      integer   i, in, ir, iaquiHIST , iaquiZ
39      integer         ib, isot
40      real*8            argumento
41      real*8            tauinf(nl_cts)
42      real*8            con(nzy_cts), coninf
43      real*8            c1, c2 , ccc
44      real*8            t1, t2
45      real*8            p1, p2
46      real*8            mr1, mr2
47      real*8            st1, st2
48      real*8            c1box(nbox_max), c2box(nbox_max)
49      real*8            ff      ! to avoid too small numbers
50      real*8            st, beta, ts
51      real*8    tyd(nzy_cts)
52      real*8            correc
53      real*8            deltanudbl, deltazdbl
54      real*8          yy
55
56c     external function
57      external        we_clean
58      real*8          we_clean
59
60c***********************************************************************
61
62c     
63      ib = 1
64      beta = 1.8d5
65      ibcode1 = '1'
66      isot = 1
67      deltanudbl = dble(deltanu(1,1))
68      deltazdbl = dble(deltaz_cts)
69      ff=1.0d10
70
71ccc   
72      do i=1,nzy_cts_real
73         tyd(i) = dble(ty_cts(i))
74         con(i) =  dble( co2y_cts(i) * imr(isot) )
75         correc = 2.d0 * dexp( -ee*dble(elow(isot,2))/tyd(i) )
76         con(i) = con(i) * ( 1.d0 - correc )
77         mr_cts(i) = dble(co2y_cts(i)/nty_cts(i))
78      end do
79      coninf = dble( con(nzy_cts_real) /
80     @     log( con(nzy_cts_real-1) / con(nzy_cts_real) ) )
81                                ! Correccion pequeña para la FB, nos la saltamos
82                                !call mztf_correccion_cts ( coninf, con, ib )
83
84ccc   
85      call gethist_03 ( 1 )
86
87c     
88c     tauinf
89c     
90      call initial
91
92      iaquiHIST = nhist/2
93      iaquiZ = nzy_cts_real - 2
94
95      do i=nl_cts_real,1,-1
96
97         if(i.eq.nl_cts_real)then
98
99            call intzhunt_cts (iaquiZ, zl_cts(i), nzy_cts_real,
100     @           c2,p2,mr2,t2, con)
101            do kr=1,nbox
102               ta(kr)=t2
103            end do
104            call interstrhunt (iaquiHIST, st2,t2,ka,ta)
105                                ! Check interpolation errors :
106            if (c2.le.0.0d0) then
107               ierr=15
108               varerr=c2
109               return
110            elseif (p2.le.0.0d0) then
111               ierr=16
112               varerr=p2
113               return
114            elseif (mr2.le.0.0d0) then
115               ierr=17
116               varerr=mr2
117               return
118            elseif (t2.le.0.0d0) then
119               ierr=18
120               varerr=t2
121               return
122            elseif (st2.le.0.0d0) then
123               ierr=19
124               varerr=st2
125               return
126            endif
127                                !
128            aa = p2 * coninf * mr2 * (st2 * ff)
129            cc = coninf * st2
130            dd = t2 * coninf * st2
131            do kr=1,nbox
132               ccbox(kr) = coninf * ka(kr)
133               ddbox(kr) = t2 * ccbox(kr)
134               c2box(kr) = c2 * ka(kr) * deltazdbl
135            end do
136            c2 = c2 * st2 * deltazdbl
137
138         else
139
140            call intzhunt_cts (iaquiZ, zl_cts(i), nzy_cts_real,
141     @           c1,p1,mr1,t1, con)
142            do kr=1,nbox
143               ta(kr)=t1
144            end do
145            call interstrhunt (iaquiHIST, st1,t1,ka,ta)
146            do kr=1,nbox
147               c1box(kr) = c1 * ka(kr) * deltazdbl
148            end do
149            c1 = c1 * st1 * deltazdbl
150            aa = aa + ( p1*mr1*(c1*ff) + p2*mr2*(c2*ff)) / 2.d0
151            cc = cc + ( c1 + c2 ) / 2.d0
152            ccc = ( c1 + c2 ) / 2.d0
153            dd = dd + ( t1*c1 + t2*c2 ) / 2.d0
154            do kr=1,nbox
155               ccbox(kr) = ccbox(kr) +
156     @              ( c1box(kr) + c2box(kr) )/2.d0
157               ddbox(kr) = ddbox(kr) +
158     @              ( t1*c1box(kr)+t2*c2box(kr) )/2.d0
159            end do
160
161            mr2 = mr1
162            c2=c1
163            do kr=1,nbox
164               c2box(kr) = c1box(kr)
165            end do
166            t2=t1
167            p2=p1
168         end if
169
170         pp = aa / (cc*ff)
171
172         ts = dd/cc
173         do kr=1,nbox
174            ta(kr) = ddbox(kr) / ccbox(kr)
175         end do
176         call interstrhunt(iaquiHIST, st,ts,ka,ta)
177         call intershphunt(iaquiHIST, alsa,alda,ta)
178
179c     
180         eqw=0.0d0
181         do  kr=1,nbox
182            yy = ccbox(kr) * beta
183            w = we_clean ( yy, pp, alsa(kr),alda(kr) )
184            eqw = eqw + no(kr)*w
185         end do
186
187         argumento = eqw / deltanudbl
188         tauinf(i) = dexp( - argumento )
189
190         if (i.eq.nl_cts_real) then
191            taustar11_cts(i) = 0.0d0
192         else
193            taustar11_cts(i) = deltanudbl * (tauinf(i+1)-tauinf(i))
194     @           / ( beta * ccc )
195         endif
196
197      end do
198
199
200      call mzescape_normaliz_02 ( taustar11_cts, nl_cts_real, 2 )
201
202c     end
203      return
204      end
205
206
207c     *** Old MZTUD110_dlvr11_03.f
208
209c***********************************************************************
210      subroutine MZTUD110( ierr, varerr )
211c***********************************************************************
212
213      implicit none
214
215      include 'datafile.h'
216      include 'nlte_paramdef.h'
217      include 'nlte_commons.h'
218
219
220c     arguments
221      integer         ierr      ! o
222      real*8          varerr    ! o
223
224c     local variables and constants
225      integer   i, in, ir, iaquiHIST , iaquiZ
226      integer         ib, isot
227      real*8            tau(nl,nl), argumento
228      real*8            tauinf(nl)
229      real*8            con(nzy), coninf
230      real*8            c1, c2
231      real*8            t1, t2
232      real*8            p1, p2
233      real*8            mr1, mr2
234      real*8            st1, st2
235      real*8            c1box(nbox_max), c2box(nbox_max)
236      real*8            ff      ! to avoid too small numbers
237      real*8            tvtbs(nzy)
238      real*8            st, beta, ts
239      real*8    zld(nl), zyd(nzy), deltazdbl
240      real*8            correc
241      real*8            deltanudbl
242      real*8          maxtau, yy
243
244c     external function
245      external        we_clean
246      real*8          we_clean
247
248c***********************************************************************
249
250c     
251      ib = 1
252      beta = 1.8d5
253      ibcode1 = '1'
254      isot = 1
255      deltanudbl = dble(deltanu(1,1))
256      deltazdbl = dble(deltaz)
257      ff=1.0d10
258
259ccc   
260      do i=1,nzy
261         zyd(i) = dble(zy(i))
262      enddo
263      do i=1,nl
264         zld(i) = dble(zl(i))
265      enddo
266      call interhuntdp ( tvtbs,zyd,nzy, v626t1,zld,nl, 1 )
267      do i=1,nzy
268         con(i) =  dble( co2y(i) * imr(isot) )
269         correc = 2.d0 * dexp( -ee*dble(elow(isot,2))/tvtbs(i) )
270         con(i) = con(i) * ( 1.d0 - correc )
271         mr(i) = dble(co2y(i)/nty(i))
272      end do
273      coninf = dble( con(nzy) / log( con(nzy-1) / con(nzy) ) )
274      call mztf_correccion ( coninf, con, ib )
275
276ccc   
277      call gethist_03 ( 1 )
278
279c     
280c     tauinf
281c     
282      call initial
283
284      iaquiHIST = nhist/2
285      iaquiZ = nzy - 2
286
287      do i=nl,1,-1
288
289         if(i.eq.nl)then
290
291            call intzhunt (iaquiZ, zl(i),c2,p2,mr2,t2, con)
292            do kr=1,nbox
293               ta(kr)=t2
294            end do
295            call interstrhunt (iaquiHIST, st2,t2,ka,ta)
296                                ! Check interpolation errors :
297            if (c2.le.0.0d0) then
298               ierr=15
299               varerr=c2
300               return
301            elseif (p2.le.0.0d0) then
302               ierr=16
303               varerr=p2
304               return
305            elseif (mr2.le.0.0d0) then
306               ierr=17
307               varerr=mr2
308               return
309            elseif (t2.le.0.0d0) then
310               ierr=18
311               varerr=t2
312               return
313            elseif (st2.le.0.0d0) then
314               ierr=19
315               varerr=st2
316               return
317            endif
318                                !
319            aa = p2 * coninf * mr2 * (st2 * ff)
320            cc = coninf * st2
321            dd = t2 * coninf * st2
322            do kr=1,nbox
323               ccbox(kr) = coninf * ka(kr)
324               ddbox(kr) = t2 * ccbox(kr)
325               c2box(kr) = c2 * ka(kr) * deltazdbl
326            end do
327            c2 = c2 * st2 * deltazdbl
328
329         else
330
331            call intzhunt (iaquiZ, zl(i),c1,p1,mr1,t1, con)
332            do kr=1,nbox
333               ta(kr)=t1
334            end do
335            call interstrhunt (iaquiHIST, st1,t1,ka,ta)
336            do kr=1,nbox
337               c1box(kr) = c1 * ka(kr) * deltazdbl
338            end do
339            c1 = c1 * st1 * deltazdbl
340            aa = aa + ( p1*mr1*(c1*ff) + p2*mr2*(c2*ff)) / 2.d0
341            cc = cc + ( c1 + c2 ) / 2.d0
342            dd = dd + ( t1*c1 + t2*c2 ) / 2.d0
343            do kr=1,nbox
344               ccbox(kr) = ccbox(kr) +
345     @              ( c1box(kr) + c2box(kr) )/2.d0
346               ddbox(kr) = ddbox(kr) +
347     @              ( t1*c1box(kr)+t2*c2box(kr) )/2.d0
348            end do
349
350            mr2 = mr1
351            c2=c1
352            do kr=1,nbox
353               c2box(kr) = c1box(kr)
354            end do
355            t2=t1
356            p2=p1
357         end if
358
359         pp = aa / (cc*ff)
360
361         ts = dd/cc
362         do kr=1,nbox
363            ta(kr) = ddbox(kr) / ccbox(kr)
364         end do
365         call interstrhunt(iaquiHIST, st,ts,ka,ta)
366         call intershphunt(iaquiHIST, alsa,alda,ta)
367
368c     
369         eqw=0.0d0
370         do  kr=1,nbox
371            yy = ccbox(kr) * beta
372            w = we_clean ( yy, pp, alsa(kr),alda(kr) )
373            eqw = eqw + no(kr)*w
374         end do
375
376         argumento = eqw / deltanudbl
377         tauinf(i) = dexp( - argumento )
378
379      end do
380
381
382c     
383c     tau
384c     
385
386      iaquiZ = 2
387      do 1 in=1,nl-1
388
389         call initial
390         call intzhunt (iaquiZ, zl(in), c1,p1,mr1,t1, con)
391         do kr=1,nbox
392            ta(kr) = t1
393         end do
394         call interstrhunt (iaquiHIST, st1,t1,ka,ta)
395         do kr=1,nbox
396            c1box(kr) = c1 * ka(kr) * deltazdbl
397         end do
398         c1 = c1 * st1 * deltazdbl
399
400         do 2 ir=in,nl-1
401
402            if (ir.eq.in) then
403               tau(in,ir) = 1.d0
404               goto 2
405            end if
406
407            call intzhunt (iaquiZ, zl(ir), c2,p2,mr2,t2, con)
408            do kr=1,nbox
409               ta(kr) = t2
410            end do
411            call interstrhunt (iaquiHIST, st2,t2,ka,ta)
412            do kr=1,nbox
413               c2box(kr) = c2 * ka(kr) * deltazdbl
414            end do
415            c2 = c2 * st2 * deltazdbl
416
417            aa = aa + ( p1*mr1*(c1*ff) + p2*mr2*(c2*ff)) / 2.d0
418            cc = cc + ( c1 + c2 ) / 2.d0
419            dd = dd + ( t1*c1 + t2*c2 ) / 2.d0
420            do kr=1,nbox
421               ccbox(kr) = ccbox(kr) +
422     $              ( c1box(kr) + c2box(kr) ) / 2.d0
423               ddbox(kr) = ddbox(kr) +
424     $              ( t1*c1box(kr) + t2*c2box(kr) ) / 2.d0
425            end do
426
427            mr1=mr2
428            t1=t2
429            c1=c2
430            p1=p2
431            do kr=1,nbox
432               c1box(kr) = c2box(kr)
433            end do
434
435            pp = aa / (cc * ff)
436
437            ts = dd/cc
438            do kr=1,nbox
439               ta(kr) = ddbox(kr) / ccbox(kr)
440            end do
441            call interstrhunt(iaquiHIST, st,ts,ka,ta)
442            call intershphunt(iaquiHIST, alsa,alda,ta)
443c     
444            eqw=0.0d0
445            do kr=1,nbox
446               yy = ccbox(kr) * beta
447               w = we_clean ( yy, pp, alsa(kr),alda(kr) )
448               eqw = eqw + no(kr)*w
449            end do
450
451            argumento = eqw / deltanudbl
452            tau(in,ir) = exp( - argumento )
453
454
455 2       continue
456
457 1    continue
458
459
460c     
461c     tau(in,ir) for n>r
462c     
463
464      in=nl
465
466      call initial
467
468      iaquiZ = nzy - 2
469      call intzhunt (iaquiZ, zl(in), c1,p1,mr1,t1, con)
470      do kr=1,nbox
471         ta(kr) = t1
472      end do
473      call interstrhunt (iaquiHIST,st1,t1,ka,ta)
474      do kr=1,nbox
475         c1box(kr) = c1 * ka(kr) * deltazdbl
476      end do
477      c1 = c1 * st1 * deltazdbl
478
479      do 4 ir=in-1,1,-1
480
481         call intzhunt (iaquiZ, zl(ir), c2,p2,mr2,t2, con)
482         do kr=1,nbox
483            ta(kr) = t2
484         end do
485         call interstrhunt (iaquiHIST, st2,t2,ka,ta)
486         do kr=1,nbox
487            c2box(kr) = c2 * ka(kr) * deltazdbl
488         end do
489         c2 = c2 * st2 * deltazdbl
490
491         aa = aa + ( p1*mr1*(c1*ff) + p2*mr2*(c2*ff)) / 2.d0
492         cc = cc + ( c1 + c2 ) / 2.d0
493         dd = dd + ( t1*c1 + t2*c2 ) / 2.d0
494         do kr=1,nbox
495            ccbox(kr) = ccbox(kr) +
496     $           ( c1box(kr) + c2box(kr) ) / 2.d0
497            ddbox(kr) = ddbox(kr) +
498     $           ( t1*c1box(kr) + t2*c2box(kr) ) / 2.d0
499         end do
500
501         mr1=mr2
502         c1=c2
503         t1=t2
504         p1=p2
505         do kr=1,nbox
506            c1box(kr) = c2box(kr)
507         end do
508
509         pp = aa / (cc * ff)
510         ts = dd / cc
511         do kr=1,nbox
512            ta(kr) = ddbox(kr) / ccbox(kr)
513         end do
514         call interstrhunt (iaquiHIST, st,ts,ka,ta)
515         call intershphunt (iaquiHIST, alsa,alda,ta)
516
517c     
518
519         eqw=0.0d0
520         do kr=1,nbox
521            yy = ccbox(kr) * beta
522            w = we_clean ( yy, pp, alsa(kr),alda(kr) )
523            eqw = eqw + no(kr)*w
524         end do
525
526         argumento = eqw / deltanudbl
527         tau(in,ir) = exp( - argumento )
528
529
530 4    continue
531
532c     
533c     
534c     
535      do in=nl-1,2,-1
536         do ir=in-1,1,-1
537            tau(in,ir) = tau(ir,in)
538         end do
539      end do
540
541c     
542c     Tracking potential numerical errors
543c     
544      maxtau = 0.0d0
545      do in=nl-1,2,-1
546         do ir=in-1,1,-1
547            maxtau = max( maxtau, tau(in,ir) )
548         end do
549      end do
550      if (maxtau .gt. 1.0d0) then
551         ierr = 13
552         varerr = maxtau
553         return
554      endif
555
556
557c     
558      call MZCUD110 ( tauinf,tau )
559
560c     end
561      return
562      end
563
564
565c     *** Old file MZCUD_dlvr11.f ***
566
567c***********************************************************************
568
569      subroutine MZCUD110 ( tauinf,tau )
570
571c***********************************************************************
572
573      implicit none
574
575      include 'nlte_paramdef.h'
576      include 'nlte_commons.h'
577
578c     arguments
579      real*8            tau(nl,nl) ! i
580      real*8            tauinf(nl) ! i
581
582
583c     local variables
584      integer   i, in, ir
585      real*8            a(nl,nl), cf(nl,nl), pideltanu, deltazdp, pi
586
587c***********************************************************************
588
589      pi = 3.141592
590      pideltanu = pi * dble(deltanu(1,1))
591      deltazdp = 2.0d5 * dble(deltaz)
592
593      do in=1,nl
594         do ir=1,nl
595            cf(in,ir) = 0.0d0
596            c110(in,ir) = 0.0d0
597            a(in,ir) = 0.0d0
598         end do
599         vc110(in) = 0.0d0
600      end do
601
602c     
603      do in=1,nl
604         do ir=1,nl
605
606            if (ir.eq.1) then
607               cf(in,ir) = tau(in,ir) - tau(in,1)
608            elseif (ir.eq.nl) then
609               cf(in,ir) = tauinf(in) - tau(in,ir-1)
610            else
611               cf(in,ir) = tau(in,ir) - tau(in,ir-1)
612            end if
613
614         end do
615      end do
616
617c     
618      do in=2,nl-1
619         do ir=1,nl
620            if (ir.eq.in+1) a(in,ir) = -1.d0
621            if (ir.eq.in-1) a(in,ir) = +1.d0
622            a(in,ir) = a(in,ir) / deltazdp
623         end do
624      end do
625
626c     
627      do in=1,nl
628         do ir=1,nl
629            cf(in,ir) = cf(in,ir) * pideltanu
630         end do
631      end do
632
633      do in=2,nl-1
634         do ir=1,nl
635            do i=1,nl
636               c110(in,ir) = c110(in,ir) + a(in,i) * cf(i,ir)
637            end do
638         end do
639      end do
640
641      do in=2,nl-1
642         vc110(in) =  pideltanu/deltazdp *
643     @        ( tau(in-1,1) - tau(in+1,1) )
644      end do
645
646
647c     
648      do in=2,nl-1
649         c110(in,nl-2) = c110(in,nl-2) - c110(in,nl)
650         c110(in,nl-1) = c110(in,nl-1) + 2.d0*c110(in,nl)
651      end do
652
653c     end
654      return
655      end
656
657
658c     *** Old MZMC121_dlvr11_03.f ***
659
660c***********************************************************************
661
662      subroutine MZMC121
663
664c***********************************************************************
665
666      implicit none
667
668                                ! common variables & constants
669     
670      include 'nlte_paramdef.h'
671      include 'nlte_commons.h'
672
673                                ! local variables
674
675      real*8 cax1(nl,nl)
676      real*8 v1(nl), cm_factor, vc_factor
677      real nuaux1, nuaux2, nuaux3
678      real*8 faux2,faux3, daux2,daux3
679
680      integer i,j,ik,ib
681
682************************************************************************
683
684      c121(1:nl,1:nl)=0.d0
685!      call zerom (c121,nl)
686      vc121(1:nl)=0.d0
687!      call zerov (vc121,nl)
688
689      nuaux1 = nu(1,2) - nu(1,1) ! 667.75
690      nuaux2 = nu12_0200-nu(1,1) ! 618.03
691      nuaux3 = nu12_1000-nu(1,1) ! 720.81
692      faux2 = dble(nuaux2/nuaux1)
693      faux3 = dble(nuaux3/nuaux1)
694      daux2 = dble(nuaux1-nuaux2)
695      daux3 = dble(nuaux1-nuaux3)
696
697      do 11, ik=1,3
698
699         ib=ik+1
700         cax1(1:nl,1:nl)=0.d0
701!         call zerom (cax1,nl)
702         call MZTUD121 ( cax1,v1, ib )
703
704         do i=1,nl
705
706            if(ik.eq.1)then
707               cm_factor = faux2**2.d0 * exp( daux2*ee/dble(t(i)) )
708               vc_factor = 1.d0/faux2
709            elseif(ik.eq.2)then
710               cm_factor = 1.d0
711               vc_factor = 1.d0
712            elseif(ik.eq.3)then
713               cm_factor = faux3**2.d0 * exp( daux3*ee/dble(t(i)) )
714               vc_factor = 1.d0 / faux3
715            else
716               write (*,*) ' Error in 626 hot band index  ik =', ik
717               stop ' ik can only be = 2,3,4.   Check needed.'
718            end if
719            do j=1,nl
720               c121(i,j) = c121(i,j) + cax1(i,j) * cm_factor
721            end do
722
723            vc121(i) = vc121(i) + v1(i) * vc_factor
724
725         end do
726
727 11   continue
728
729      return
730      end
731
732
733c     *** Old MZTUD121_dlvr11_03.f ***
734
735c***********************************************************************
736      subroutine MZTUD121 ( cf,vc, ib )
737c***********************************************************************
738
739      implicit none
740
741      include 'datafile.h'
742      include 'nlte_paramdef.h'
743      include 'nlte_commons.h'
744     
745
746c     arguments
747      real*8    cf(nl,nl)       ! o.
748      real*8            vc(nl)  ! o
749      integer           ib      ! i
750
751
752c     local variables and constants
753      integer   i, in, ir, iaquiHIST, iaquiZ
754      integer   nmu
755      parameter         (nmu = 8)
756      real*8            tau(nl,nl), argumento, deltazdbl
757      real*8            tauinf(nl)
758      real*8            con(nzy), coninf
759      real*8            c1, c2
760      real*8            t1, t2
761      real*8            p1, p2
762      real*8            mr1, mr2
763      real*8            st1, st2
764      real*8            c1box(70), c2box(70)
765      real*8            ff      ! to avoid too small numbers
766      real*8            tvtbs(nzy)
767      real*8            st, beta, ts
768
769      real*8    zld(nl), zyd(nzy)
770      real*8            correc
771      real*8    deltanudbl
772      integer         isot
773      real*8          yy
774
775c     external function
776      external        we_clean
777      real*8          we_clean
778
779
780c     formats
781 101  format(i1)
782c***********************************************************************
783
784c     some values
785      beta = 1.8d5
786      isot = 1
787      write (ibcode1,101) ib
788      deltanudbl = dble( deltanu(isot,ib) )
789      ff=1.0d10
790      deltazdbl = dble(deltaz)
791
792ccc   
793ccc   
794ccc   
795      do i=1,nl
796         zld(i) = dble(zl(i))
797      enddo
798      do i=1,nzy
799         zyd(i) = dble(zy(i))
800      enddo
801
802      call interhuntdp ( tvtbs,zyd,nzy, v626t1,zld,nl, 1 )
803
804      do i=1,nzy
805         con(i) =  dble( co2y(i) * imr(isot) )
806         correc = 2.d0 * exp( -ee*dble(elow(isot,2))/tvtbs(i) )
807         con(i) = con(i) * ( 1.d0 - correc )
808         mr(i) = dble( co2y(i) / nty(i) )
809      end do
810
811      coninf = dble( con(nzy) / log( con(nzy-1) / con(nzy) ) )
812      call mztf_correccion ( coninf, con, ib )
813
814ccc   
815      call gethist_03 ( ib )
816
817
818c     
819c     tauinf(nl)
820c     
821      call initial
822
823      iaquiZ = nzy - 2
824      iaquiHIST = nhist / 2
825
826      do i=nl,1,-1
827
828         if(i.eq.nl)then
829
830            call intzhunt ( iaquiZ, zl(i),c2,p2,mr2,t2, con)
831            do kr=1,nbox
832               ta(kr)=t2
833            end do
834            call interstrhunt (iaquiHIST, st2,t2,ka,ta)
835            aa = p2 * coninf * mr2 * (st2 * ff)
836            cc = coninf * st2
837            dd = t2 * coninf * st2
838            do kr=1,nbox
839               ccbox(kr) = coninf * ka(kr)
840               ddbox(kr) = t2 * ccbox(kr)
841               c2box(kr) = c2 * ka(kr) * deltazdbl
842            end do
843            c2 = c2 * st2 * deltazdbl
844
845         else
846            call intzhunt ( iaquiZ, zl(i),c1,p1,mr1,t1, con)
847            do kr=1,nbox
848               ta(kr)=t1
849            end do
850            call interstrhunt (iaquiHIST, st1,t1,ka,ta)
851            do kr=1,nbox
852               c1box(kr) = c1 * ka(kr) * deltazdbl
853            end do
854            c1 = c1 * st1 * deltazdbl
855            aa = aa + ( p1*mr1*(c1*ff) + p2*mr2*(c2*ff)) / 2.d0
856            cc = cc + ( c1 + c2 ) / 2.d0
857            dd = dd + ( t1*c1 + t2*c2 ) / 2.d0
858            do kr=1,nbox
859               ccbox(kr) = ccbox(kr) +
860     @              ( c1box(kr) + c2box(kr) )/2.d0
861               ddbox(kr) = ddbox(kr) +
862     @              ( t1*c1box(kr)+t2*c2box(kr) )/2.d0
863            end do
864
865            mr2 = mr1
866            c2=c1
867            do kr=1,nbox
868               c2box(kr) = c1box(kr)
869            end do
870            t2=t1
871            p2=p1
872         end if
873
874         pp = aa / (cc*ff)
875
876         ts = dd/cc
877         do kr=1,nbox
878            ta(kr) = ddbox(kr) / ccbox(kr)
879         end do
880         call interstrhunt(iaquiHIST, st,ts,ka,ta)
881         call intershphunt(iaquiHIST, alsa,alda,ta)
882
883c     
884
885         eqw = 0.0d0
886         do  kr=1,nbox
887            yy = ccbox(kr) * beta
888            w = we_clean ( yy, pp, alsa(kr),alda(kr) )
889            eqw = eqw + no(kr)*w
890         end do
891
892         argumento = eqw / deltanudbl
893         tauinf(i) = dexp( - argumento )
894
895
896      end do                    ! i continue
897
898
899c     
900c     tau(in,ir) for n<=r
901c     
902
903      iaquiZ = 2
904      do 1 in=1,nl-1
905
906         call initial
907         call intzhunt ( iaquiZ, zl(in), c1,p1,mr1,t1, con)
908         do kr=1,nbox
909            ta(kr) = t1
910         end do
911         call interstrhunt (iaquiHIST, st1,t1,ka,ta)
912         do kr=1,nbox
913            c1box(kr) = c1 * ka(kr) * deltazdbl
914         end do
915         c1 = c1 * st1 * deltazdbl
916
917         do 2 ir=in,nl-1
918
919            if (ir.eq.in) then
920               tau(in,ir) = 1.d0
921               goto 2
922            end if
923
924            call intzhunt ( iaquiZ, zl(ir), c2,p2,mr2,t2, con)
925            do kr=1,nbox
926               ta(kr) = t2
927            end do
928            call interstrhunt (iaquiHIST, st2,t2,ka,ta)
929            do kr=1,nbox
930               c2box(kr) = c2 * ka(kr) * deltazdbl
931            end do
932            c2 = c2 * st2 * deltazdbl
933
934            aa = aa + ( p1*mr1*(c1*ff) + p2*mr2*(c2*ff)) / 2.d0
935            cc = cc + ( c1 + c2 ) / 2.d0
936            dd = dd + ( t1*c1 + t2*c2 ) / 2.d0
937            do kr=1,nbox
938               ccbox(kr) = ccbox(kr) +
939     $              ( c1box(kr) + c2box(kr) ) / 2.d0
940               ddbox(kr) = ddbox(kr) +
941     $              ( t1*c1box(kr) + t2*c2box(kr) ) / 2.d0
942            end do
943
944            mr1=mr2
945            t1=t2
946            c1=c2
947            p1=p2
948            do kr=1,nbox
949               c1box(kr) = c2box(kr)
950            end do
951
952            pp = aa / (cc * ff)
953
954            ts = dd/cc
955            do kr=1,nbox
956               ta(kr) = ddbox(kr) / ccbox(kr)
957            end do
958            call interstrhunt(iaquiHIST, st,ts,ka,ta)
959            call intershphunt(iaquiHIST, alsa,alda,ta)
960
961c     
962
963            eqw = 0.0d0
964            do kr=1,nbox
965               yy = ccbox(kr) * beta
966               w = we_clean ( yy, pp, alsa(kr),alda(kr) )
967               eqw = eqw + no(kr)*w
968            end do
969
970            argumento = eqw / deltanudbl
971            tau(in,ir) = dexp( - argumento )
972
973 2       continue
974
975 1    continue
976
977c     
978c     tau(in,ir) for n>r
979c     
980
981      in=nl
982
983      call initial
984      iaquiZ = nzy - 2
985      call intzhunt ( iaquiZ, zl(in), c1,p1,mr1,t1, con)
986      do kr=1,nbox
987         ta(kr) = t1
988      end do
989      call interstrhunt (iaquiHIST, st1,t1,ka,ta)
990      do kr=1,nbox
991         c1box(kr) = c1 * ka(kr) * deltazdbl
992      end do
993      c1 = c1 * st1 * deltazdbl
994
995      do 4 ir=in-1,1,-1
996
997         call intzhunt ( iaquiZ, zl(ir), c2,p2,mr2,t2, con)
998         do kr=1,nbox
999            ta(kr) = t2
1000         end do
1001         call interstrhunt (iaquiHIST, st2,t2,ka,ta)
1002         do kr=1,nbox
1003            c2box(kr) = c2 * ka(kr) * deltazdbl
1004         end do
1005         c2 = c2 * st2 * deltazdbl
1006
1007         aa = aa + ( p1*mr1*(c1*ff) + p2*mr2*(c2*ff)) / 2.d0
1008         cc = cc + ( c1 + c2 ) / 2.d0
1009         dd = dd + ( t1*c1 + t2*c2 ) / 2.d0
1010         do kr=1,nbox
1011            ccbox(kr) = ccbox(kr) +
1012     $           ( c1box(kr) + c2box(kr) ) / 2.d0
1013            ddbox(kr) = ddbox(kr) +
1014     $           ( t1*c1box(kr) + t2*c2box(kr) ) / 2.d0
1015         end do
1016
1017         mr1=mr2
1018         c1=c2
1019         t1=t2
1020         p1=p2
1021         do kr=1,nbox
1022            c1box(kr) = c2box(kr)
1023         end do
1024
1025         pp = aa / (cc * ff)
1026         ts = dd / cc
1027         do kr=1,nbox
1028            ta(kr) = ddbox(kr) / ccbox(kr)
1029         end do
1030         call interstrhunt (iaquiHIST, st,ts,ka,ta)
1031         call intershphunt (iaquiHIST, alsa,alda,ta)
1032
1033c     
1034         eqw=0.0d0
1035         do kr=1,nbox
1036            yy = ccbox(kr) * beta
1037            w = we_clean ( yy, pp, alsa(kr),alda(kr) )
1038            eqw = eqw + no(kr)*w
1039         end do
1040
1041         argumento = eqw / deltanudbl
1042         tau(in,ir) = dexp( - argumento )
1043
1044 4    continue
1045
1046c     
1047c     
1048c     
1049      do in=nl-1,2,-1
1050         do ir=in-1,1,-1
1051            tau(in,ir) = tau(ir,in)
1052         end do
1053      end do
1054
1055c     
1056      call MZCUD121 ( tauinf,tau, cf, vc, ib )
1057
1058
1059c     end
1060      return
1061      end
1062
1063
1064
1065c     *** Old MZCUD121_dlvr11.f ***
1066
1067c***********************************************************************
1068
1069      subroutine MZCUD121 ( tauinf,tau, c,vc, ib )
1070
1071c***********************************************************************
1072
1073      implicit none
1074
1075      include 'nlte_paramdef.h'
1076      include 'nlte_commons.h'
1077
1078
1079c     arguments
1080      real*8            c(nl,nl) ! o
1081      real*8            vc(nl)  ! o
1082      real*8            tau(nl,nl) ! i
1083      real*8            tauinf(nl) ! i
1084      integer           ib      ! i
1085
1086c     local variables
1087      integer   i, in, ir, isot
1088      real*8            a(nl,nl), cf(nl,nl), pideltanu, deltazdbl,pi
1089
1090c***********************************************************************
1091
1092      pi=3.141592
1093      isot = 1
1094      pideltanu = pi*dble(deltanu(isot,ib))
1095      deltazdbl = dble(deltaz)
1096c     
1097      do in=1,nl
1098
1099         do ir=1,nl
1100
1101            cf(in,ir) = 0.0d0
1102            c(in,ir) = 0.0d0
1103            a(in,ir) = 0.0d0
1104
1105         end do
1106
1107         vc(in) = 0.0d0
1108
1109      end do
1110
1111
1112c     
1113      do in=1,nl
1114         do ir=1,nl
1115
1116            if (ir.eq.1) then
1117               cf(in,ir) = tau(in,ir) - tau(in,1)
1118            elseif (ir.eq.nl) then
1119               cf(in,ir) = tauinf(in) - tau(in,ir-1)
1120            else
1121               cf(in,ir) = tau(in,ir) - tau(in,ir-1)
1122            end if
1123
1124         end do
1125      end do
1126
1127
1128c     
1129      do in=2,nl-1
1130         do ir=1,nl
1131            if (ir.eq.in+1) a(in,ir) = -1.d0
1132            if (ir.eq.in-1) a(in,ir) = +1.d0
1133            a(in,ir) = a(in,ir) / ( 2.d5*deltazdbl )
1134         end do
1135      end do
1136
1137c     
1138      do in=1,nl
1139         do ir=1,nl
1140            cf(in,ir) = cf(in,ir) * pideltanu
1141         end do
1142      end do
1143
1144
1145      do in=2,nl-1
1146         do ir=1,nl
1147            do i=1,nl
1148               c(in,ir) = c(in,ir) + a(in,i) * cf(i,ir)
1149            end do
1150         end do
1151         vc(in) =  pideltanu /( 2.d5*deltazdbl ) *
1152     @        ( tau(in-1,1) - tau(in+1,1) )
1153      end do
1154
1155c     
1156      do in=2,nl-1
1157         c(in,nl-2) = c(in,nl-2) - c(in,nl)
1158         c(in,nl-1) = c(in,nl-1) + 2.d0*c(in,nl)
1159      end do
1160
1161
1162c     end
1163      return
1164      end
1165
1166
1167
1168c     *** Old MZESC121_dlvr11_03.f ***
1169
1170c***********************************************************************
1171      subroutine MZESC121
1172c***********************************************************************
1173
1174      implicit none
1175
1176      include 'nlte_paramdef.h'
1177      include 'nlte_commons.h'
1178
1179
1180c     local variables
1181      integer   i
1182      real*8          factor0200, factor0220, factor1000
1183      real*8          aux_0200(nl), aux2_0200(nl)
1184      real*8          aux_0220(nl), aux2_0220(nl)
1185      real*8          aux_1000(nl), aux2_1000(nl)
1186
1187c***********************************************************************
1188
1189!      call zerov (taustar12, nl)
1190      taustar12(1:nl)=0.d0
1191      call zero2v(aux_0200,aux2_0200, nl)
1192      call zero2v(aux_0220,aux2_0220, nl)
1193      call zero2v(aux_1000,aux2_1000, nl)
1194
1195      call MZESC121sub (aux_0200,aux2_0200, 2 )
1196      call MZESC121sub (aux_0220,aux2_0220, 3 )
1197      call MZESC121sub (aux_1000,aux2_1000, 4 )
1198
1199      factor0220 = 1.d0
1200      factor0200 = dble( (nu(1,2)-nu(1,1)) / (nu12_0200-nu(1,1)) )
1201      factor1000 = dble( (nu(1,2)-nu(1,1)) / (nu12_1000-nu(1,1)) )
1202      do i=1,nl
1203         taustar12(i) = taustar12(i)
1204     @        + aux_0200(i) * factor0200
1205     @        + aux_0220(i) * factor0220
1206     @        + aux_1000(i) * factor1000
1207      enddo
1208
1209      call mzescape_normaliz ( taustar12, 2 )
1210
1211c     end
1212      return
1213      end
1214
1215
1216c     *** Old MZESC121sub_dlvr11_03.f ***
1217
1218c***********************************************************************
1219
1220      subroutine MZESC121sub (taustar,tauinf, ib )
1221
1222c***********************************************************************
1223
1224      implicit none
1225
1226      include 'datafile.h'
1227      include 'nlte_paramdef.h'
1228      include 'nlte_commons.h'
1229
1230
1231c     arguments
1232      real*8            taustar(nl) ! o
1233      real*8            tauinf(nl) ! o
1234      integer           ib      ! i
1235
1236
1237c     local variables and constants
1238      integer   i, iaquiHIST, iaquiZ, isot
1239      real*8            con(nzy), coninf
1240      real*8            c1, c2, ccc
1241      real*8            t1, t2
1242      real*8            p1, p2
1243      real*8            mr1, mr2
1244      real*8            st1, st2
1245      real*8            c1box(70), c2box(70)
1246      real*8            ff      ! to avoid too small numbers
1247      real*8            tvtbs(nzy)
1248      real*8            st, beta, ts
1249      real*8    zld(nl), zyd(nzy)
1250      real*8            correc
1251      real*8            deltanudbl, deltazdbl
1252      real*8          yy
1253
1254c     external function
1255      external        we_clean
1256      real*8  we_clean
1257
1258c     formats
1259 101  format(i1)
1260
1261c***********************************************************************
1262
1263c     
1264      beta = 1.8d5
1265      isot = 1
1266      write ( ibcode1, 101) ib
1267      deltanudbl = dble( deltanu(isot,ib) )
1268      ff=1.0d10
1269      deltazdbl = dble(deltaz)
1270
1271c     
1272      do i=1,nzy
1273         zyd(i) = dble(zy(i))
1274      enddo
1275      do i=1,nl
1276         zld(i) = dble(zl(i))
1277      enddo
1278
1279      call interhuntdp ( tvtbs,zyd,nzy, v626t1,zld,nl, 1 )
1280
1281      do i=1,nzy
1282         con(i) =  dble( co2y(i) * imr(isot) )
1283         correc = 2.d0 * dexp( -ee*dble(elow(isot,2))/tvtbs(i) )
1284         con(i) = con(i) * ( 1.d0 - correc )
1285         mr(i) = dble(co2y(i)/nty(i))
1286      end do
1287
1288      coninf = dble( con(nzy) / log( con(nzy-1) / con(nzy) ) )
1289      call mztf_correccion ( coninf, con, ib )
1290
1291c     
1292      call gethist_03 ( ib )
1293
1294c     
1295c     tauinf
1296c     
1297      call initial
1298
1299      iaquiHIST = nhist/2
1300      iaquiZ = nzy - 2
1301
1302      do i=nl,1,-1
1303
1304         if(i.eq.nl)then
1305
1306            call intzhunt (iaquiZ, zl(i),c2,p2,mr2,t2, con)
1307            do kr=1,nbox
1308               ta(kr)=t2
1309            end do
1310            call interstrhunt (iaquiHIST, st2,t2,ka,ta)
1311            aa = p2 * coninf * mr2 * (st2 * ff)
1312            cc = coninf * st2
1313            dd = t2 * coninf * st2
1314            do kr=1,nbox
1315               ccbox(kr) = coninf * ka(kr)
1316               ddbox(kr) = t2 * ccbox(kr)
1317               c2box(kr) = c2 * ka(kr) * deltazdbl
1318            end do
1319            c2 = c2 * st2 * deltazdbl
1320
1321         else
1322            call intzhunt (iaquiZ, zl(i),c1,p1,mr1,t1, con)
1323            do kr=1,nbox
1324               ta(kr)=t1
1325            end do
1326            call interstrhunt (iaquiHIST,st1,t1,ka,ta)
1327            do kr=1,nbox
1328               c1box(kr) = c1 * ka(kr) * deltazdbl
1329            end do
1330            c1 = c1 * st1 * deltazdbl
1331            aa = aa + ( p1*mr1*(c1*ff) + p2*mr2*(c2*ff)) / 2.d0
1332            cc = cc + ( c1 + c2 ) / 2.d0
1333            ccc = ( c1 + c2 ) / 2.d0
1334            dd = dd + ( t1*c1 + t2*c2 ) / 2.d0
1335            do kr=1,nbox
1336               ccbox(kr) = ccbox(kr) +
1337     @              ( c1box(kr) + c2box(kr) )/2.d0
1338               ddbox(kr) = ddbox(kr) +
1339     @              ( t1*c1box(kr)+t2*c2box(kr) )/2.d0
1340            end do
1341
1342            mr2 = mr1
1343            c2=c1
1344            do kr=1,nbox
1345               c2box(kr) = c1box(kr)
1346            end do
1347            t2=t1
1348            p2=p1
1349         end if
1350
1351         pp = aa / (cc*ff)
1352
1353         ts = dd/cc
1354         do kr=1,nbox
1355            ta(kr) = ddbox(kr) / ccbox(kr)
1356         end do
1357         call interstrhunt(iaquiHIST,st,ts,ka,ta)
1358         call intershphunt(iaquiHIST,alsa,alda,ta)
1359
1360c     
1361         eqw=0.0d0
1362         do  kr=1,nbox
1363            yy = ccbox(kr) * beta
1364            w = we_clean ( yy, pp, alsa(kr),alda(kr) )
1365            eqw = eqw + no(kr)*w
1366         end do
1367         tauinf(i) = dexp( - eqw / deltanudbl )
1368         if (tauinf(i).lt.0.d0) tauinf(i) = 0.0d0
1369
1370         if (i.eq.nl) then
1371            taustar(i) = 0.0d0
1372         else
1373            taustar(i) = deltanudbl * (tauinf(i+1)-tauinf(i))
1374     @           / ( beta * ccc  )
1375         endif
1376
1377      end do
1378
1379
1380
1381c     end
1382      return
1383      end
1384
1385
1386c     *** Old MZTVC121_dlvr11.f ***
1387
1388c***********************************************************************
1389
1390      subroutine MZTVC121
1391
1392c***********************************************************************
1393
1394      implicit none
1395
1396!!!!!!!!!!!!!!!!!!!!!!!
1397!     common variables & constants
1398
1399      include 'nlte_paramdef.h'
1400      include 'nlte_commons.h'
1401
1402
1403!     arguments
1404      integer ierr
1405      real*8 varerr
1406
1407
1408!     local variables
1409
1410      real*8 v1(nl), vc_factor
1411      integer i,ik,ib
1412
1413************************************************************************
1414
1415!      call zerov( vc121, nl )
1416      vc121(1:nl)=0.d0
1417
1418      do 11, ik=1,3
1419
1420         ib=ik+1
1421
1422         call MZTVC121sub (v1, ib, ierr,varerr )
1423
1424         do i=1,nl
1425
1426            if(ik.eq.1)then
1427               vc_factor =
1428     @              dble( (nu(1,2)-nu(1,1)) / (nu12_0200-nu(1,1)) )
1429            elseif(ik.eq.2)then
1430               vc_factor = 1.d0
1431            elseif(ik.eq.3)then
1432               vc_factor =
1433     @              dble( (nu(1,2)-nu(1,1)) / (nu12_1000-nu(1,1)) )
1434            end if
1435
1436            vc121(i) = vc121(i) + v1(i) * vc_factor
1437
1438         end do
1439
1440 11   continue
1441
1442
1443      return
1444      end
1445
1446
1447c     *** Old MZTVC121sub_dlvr11_03.f ***
1448
1449c***********************************************************************
1450c     mztf.f
1451c***********************************************************************
1452
1453      subroutine MZTVC121sub  ( vc, ib,  ierr, varerr )
1454
1455c***********************************************************************
1456
1457      implicit none
1458
1459      include 'datafile.h'
1460      include 'nlte_paramdef.h'
1461      include 'nlte_commons.h'
1462
1463
1464c     arguments
1465      real*8            vc(nl)  ! o
1466      integer           ib      ! i
1467      integer ierr              ! o
1468      real*8  varerr            ! o
1469
1470c     local variables and constants
1471      integer   i, in, ir, iaquiHIST , iaquiZ, isot
1472      integer   nmu
1473      parameter         (nmu = 8)
1474      real*8            tau(nl,nl), argumento
1475      real*8            con(nzy), coninf
1476      real*8            c1, c2
1477      real*8            t1, t2
1478      real*8            p1, p2
1479      real*8            mr1, mr2
1480      real*8            st1, st2
1481      real*8            c1box(70), c2box(70)
1482      real*8            ff      ! to avoid too small numbers
1483      real*8            tvtbs(nzy)
1484      real*8            st, beta, ts
1485      real*8    zld(nl), zyd(nzy), deltazdbl
1486      real*8            correc
1487      real*8            deltanudbl, pideltanu,pi
1488      real*8          yy
1489      real*8          minvc, maxtau
1490
1491c     external function
1492      external        we_clean
1493      real*8          we_clean
1494
1495c     formats
1496 101  format(i1)
1497
1498c***********************************************************************
1499
1500c     
1501      pi=3.141592
1502      isot = 1
1503      beta = 1.8d5
1504      write (ibcode1,101) ib
1505      deltanudbl = dble( deltanu(isot,ib) )
1506      pideltanu = pi*deltanudbl
1507      ff=1.0d10
1508      deltazdbl = dble(deltaz)
1509c     
1510c     
1511c     
1512
1513      do i=1,nzy
1514         zyd(i) = dble(zy(i))
1515      enddo
1516      do i=1,nl
1517         zld(i) = dble(zl(i))
1518      enddo
1519
1520      call interhuntdp ( tvtbs,zyd,nzy, v626t1,zld,nl, 1 )
1521
1522      do i=1,nzy
1523         con(i) =  dble( co2y(i) * imr(isot) )
1524         correc = 2.d0 * dexp( -ee*dble(elow(isot,2))/tvtbs(i) )
1525         con(i) = con(i) * ( 1.d0 - correc )
1526         mr(i) = dble(co2y(i)/nty(i))
1527      end do
1528
1529      coninf = dble( con(nzy) / log( con(nzy-1) / con(nzy) ) )
1530      call mztf_correccion ( coninf, con, ib )
1531
1532ccc   
1533      call gethist_03 ( ib )
1534
1535c     
1536c     tau(1,ir)
1537c     
1538      call initial
1539
1540      iaquiHIST = nhist/2
1541
1542      in=1
1543
1544      tau(in,1) = 1.d0
1545
1546      call initial
1547      iaquiZ = 2
1548      call intzhunt ( iaquiZ, zl(in), c1,p1,mr1,t1, con)
1549      do kr=1,nbox
1550         ta(kr) = t1
1551      end do
1552      call interstrhunt (iaquiHIST, st1,t1,ka,ta)
1553      do kr=1,nbox
1554         c1box(kr) = c1 * ka(kr) * deltazdbl
1555      end do
1556      c1 = c1 * st1 * deltazdbl
1557                                ! Check interpolation errors :
1558      if (c1.le.0.0d0) then
1559         ierr=15
1560         varerr=c1
1561         return
1562      elseif (p1.le.0.0d0) then
1563         ierr=16
1564         varerr=p1
1565         return
1566      elseif (mr1.le.0.0d0) then
1567         ierr=17
1568         varerr=mr1
1569         return
1570      elseif (t1.le.0.0d0) then
1571         ierr=18
1572         varerr=t1
1573         return
1574      elseif (st1.le.0.0d0) then
1575         ierr=19
1576         varerr=st1
1577         return
1578      endif
1579                                !
1580
1581      do 2 ir=2,nl
1582
1583         call intzhunt (iaquiZ, zl(ir), c2,p2,mr2,t2, con)
1584         do kr=1,nbox
1585            ta(kr) = t2
1586         end do
1587         call interstrhunt (iaquiHIST, st2,t2,ka,ta)
1588         do kr=1,nbox
1589            c2box(kr) = c2 * ka(kr) * deltazdbl
1590         end do
1591         c2 = c2 * st2 * deltazdbl
1592
1593         aa = aa + ( p1*mr1*(c1*ff) + p2*mr2*(c2*ff)) / 2.d0
1594         cc = cc + ( c1 + c2 ) / 2.d0
1595         dd = dd + ( t1*c1 + t2*c2 ) / 2.d0
1596         do kr=1,nbox
1597            ccbox(kr) = ccbox(kr) + ( c1box(kr) + c2box(kr) ) /2.d0
1598            ddbox(kr) = ddbox(kr) +
1599     $           ( t1*c1box(kr) + t2*c2box(kr) ) / 2.d0
1600         end do
1601
1602         mr1=mr2
1603         t1=t2
1604         c1=c2
1605         p1=p2
1606         do kr=1,nbox
1607            c1box(kr) = c2box(kr)
1608         end do
1609
1610         pp = aa / (cc * ff)
1611
1612         ts = dd/cc
1613         do kr=1,nbox
1614            ta(kr) = ddbox(kr) / ccbox(kr)
1615         end do
1616         call interstrhunt(iaquiHIST, st,ts,ka,ta)
1617         call intershphunt(iaquiHIST, alsa,alda,ta)
1618
1619         eqw=0.0d0
1620         do kr=1,nbox
1621            yy = ccbox(kr) * beta
1622            w = we_clean ( yy, pp, alsa(kr),alda(kr) )
1623            eqw = eqw + no(kr)*w
1624         end do
1625
1626         argumento = eqw / deltanudbl
1627         tau(in,ir) = dexp( - argumento )
1628
1629 2    continue
1630
1631
1632c     
1633c     
1634c     
1635      do in=nl,2,-1
1636         tau(in,1) = tau(1,in)
1637      end do
1638
1639c     
1640      vc(1) = 0.0d0
1641      vc(nl) = 0.0d0
1642      do in=2,nl-1
1643         vc(in) =  pideltanu /( 2.d5*deltazdbl ) *
1644     @        ( tau(in-1,1) - tau(in+1,1) )
1645         if (vc(in) .lt. 0.0d0) vc(in) = vc(in-1)
1646      end do
1647
1648c     
1649c     Tracking potential numerical errors
1650c     
1651      minvc = 1.d6
1652      maxtau = tau(nl,1)
1653      do in=2,nl-1
1654         minvc = min( minvc, vc(in) )
1655         maxtau = max( maxtau, tau(in,1) )
1656      end do
1657      if (maxtau .gt. 1.0d0) then
1658         ierr = 13
1659         varerr = maxtau
1660         return
1661      else if (minvc .lt. 0.0d0) then
1662         ierr = 14
1663         varerr = minvc
1664         return
1665      endif
1666
1667c     end
1668      return
1669      end
1670
1671
1672
1673
1674
1675
1676
1677
1678
Note: See TracBrowser for help on using the repository browser.