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