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

Last change on this file since 2541 was 2398, checked in by emillour, 5 years ago

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