source: trunk/LMDZ.VENUS/libf/phyvenus/nlte_calc.F @ 3567

Last change on this file since 3567 was 1310, checked in by slebonnois, 10 years ago

SL: VENUS VERTICAL EXTENSION. NLTE and thermospheric processes, to be run with 78 levels and specific inputs.
File size: 43.6 KB
RevLine 
[1310]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 (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
39c     old arguments
40      integer   ierr      ! o
41      real*8    varerr    ! o
42
43c     local variables and constants
44      integer   i, iaquiHIST , iaquiZ
45      integer   isot
46      real*8    argumento
47      real*8    tauinf(nl_cts)
48      real*8    con(nzy_cts), coninf
49      real*8    c1, c2 , ccc
50      real*8    t1, t2
51      real*8    p1, p2
52      real*8    mr1, mr2
53      real*8    st1, st2
54      real*8    c1box(nbox_max), c2box(nbox_max)
55      real*8    ff      ! to avoid too small numbers
56      real*8    st, beta, ts
57      real*8    tyd(nzy_cts)
58      real*8    correc
59      real*8    deltanudbl, deltazdbl
60      real*8    yy
61 
62c     external function
63      external  we_clean
64      real*8    we_clean
65       
66c***********************************************************************
67      ierr = 0
68      varerr = 0.d0
69
70      beta = 1.8d5
71      ibcode1 = '1'
72      isot = 1
73      deltanudbl = dble(deltanu(1,1))
74      deltazdbl = dble(deltaz_cts)
75      ff=1.0d10
76
77ccc   
78      do i=1,nzy_cts_real
79         tyd(i) = dble(ty_cts(i))
80         con(i) =  dble( co2y_cts(i) * imr(isot) )
81         correc = 2.d0 * dexp( -ee*dble(elow(isot,2))/tyd(i) )
82         con(i) = con(i) * ( 1.d0 - correc )
83         mr_cts(i) = dble(co2y_cts(i)/nty_cts(i))
84      end do
85      if ( con(nzy_cts_real) .le. 0.0d0 ) then
86         ierr = 33
87         varerr = con(nzy_cts_real)
88         return
89      elseif ( con(nzy_cts_real-1) .le. con(nzy_cts_real) ) then
90         write (*,*) ' WARNING in MZESC110 '
91         write (*,*) '    [CO2] growing with altitude at TOA.'
92         write (*,*) '    [CO2] @ TOA = ', con(nzy_cts_real)
93         coninf = dble( con(nzy_cts_real) )
94      else
95         coninf = dble( con(nzy_cts_real) /
96     @        log( con(nzy_cts_real-1) / con(nzy_cts_real) ) )
97      endif
98ccc   
99      call gethist_03 ( 1 )
100
101c     
102c     tauinf
103c     
104      call initial
105
106      iaquiHIST = nhist/2
107      iaquiZ = nzy_cts_real - 2
108
109      do i=nl_cts_real,1,-1
110
111         if(i.eq.nl_cts_real)then
112
113            call intzhunt_cts (iaquiZ, zl_cts(i), nzy_cts_real,
114     @           c2,p2,mr2,t2, con)
115            do kr=1,nbox
116               ta(kr)=t2
117            end do
118            call interstrhunt (iaquiHIST, st2,t2,ka,ta)
119                                ! Check interpolation errors :
120            if (c2.le.0.0d0) then
121               ierr=15
122               varerr=c2
123               return
124            elseif (p2.le.0.0d0) then
125               ierr=16
126               varerr=p2
127               return
128            elseif (mr2.le.0.0d0) then
129               ierr=17
130               varerr=mr2
131               return
132            elseif (t2.le.0.0d0) then
133               ierr=18
134               varerr=t2
135               return
136            elseif (st2.le.0.0d0) then
137               ierr=19
138               varerr=st2
139               return
140            endif
141                                !
142            aa = p2 * coninf * mr2 * (st2 * ff)
143            cc = coninf * st2
144            dd = t2 * coninf * st2
145            do kr=1,nbox
146               ccbox(kr) = coninf * ka(kr)
147               ddbox(kr) = t2 * ccbox(kr)
148               c2box(kr) = c2 * ka(kr) * deltazdbl
149            end do
150            c2 = c2 * st2 * deltazdbl
151
152         else
153
154            call intzhunt_cts (iaquiZ, zl_cts(i), nzy_cts_real,
155     @           c1,p1,mr1,t1, con)
156            do kr=1,nbox
157               ta(kr)=t1
158            end do
159            call interstrhunt (iaquiHIST, st1,t1,ka,ta)
160            do kr=1,nbox
161               c1box(kr) = c1 * ka(kr) * deltazdbl
162            end do
163            c1 = c1 * st1 * deltazdbl
164            aa = aa + ( p1*mr1*(c1*ff) + p2*mr2*(c2*ff)) / 2.d0
165            cc = cc + ( c1 + c2 ) / 2.d0
166            ccc = ( c1 + c2 ) / 2.d0
167            dd = dd + ( t1*c1 + t2*c2 ) / 2.d0
168            do kr=1,nbox
169               ccbox(kr) = ccbox(kr) +
170     @              ( c1box(kr) + c2box(kr) )/2.d0
171               ddbox(kr) = ddbox(kr) +
172     @              ( t1*c1box(kr)+t2*c2box(kr) )/2.d0
173            end do
174
175            mr2 = mr1
176            c2=c1
177            do kr=1,nbox
178               c2box(kr) = c1box(kr)
179            end do
180            t2=t1
181            p2=p1
182         end if
183
184         pp = aa / (cc*ff)
185
186         ts = dd/cc
187         do kr=1,nbox
188            ta(kr) = ddbox(kr) / ccbox(kr)
189         end do
190         call interstrhunt(iaquiHIST, st,ts,ka,ta)
191         call intershphunt(iaquiHIST, alsa,alda,ta)
192
193c     
194         eqw=0.0d0
195         do  kr=1,nbox
196            yy = ccbox(kr) * beta
197            w = we_clean ( yy, pp, alsa(kr),alda(kr) )
198            eqw = eqw + no(kr)*w
199         end do
200
201         argumento = eqw / deltanudbl
202         tauinf(i) = dexp( - argumento )
203
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
265
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 * dexp( -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) = dexp( - 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               stop ' ik can only be = 2,3,4.   Check needed.'
770            end if
771            do j=1,nl
772               c121(i,j) = c121(i,j) + cax1(i,j) * cm_factor
773            end do
774
775            vc121(i) = vc121(i) + v1(i) * vc_factor
776
777         end do
778
779 11   continue
780
781      return
782      end
783
784
785c     *** Old MZTUD121_dlvr11_03.f ***
786
787c***********************************************************************
788      subroutine MZTUD121 ( cf,vc, ib, ierr, varerr )
789c***********************************************************************
790
791      implicit none
792
793#include "nlte_paramdef.h"
794#include "nlte_commons.h"
795     
796
797c     arguments
798      real*8    cf(nl,nl) ! o.
799      real*8    vc(nl)    ! o
800      integer   ib        ! i
801      integer   ierr      ! o
802      real*8    varerr    ! o
803
804c     local variables and constants
805      integer   i, in, ir, iaquiHIST, iaquiZ
806      integer   nmu
807      parameter (nmu = 8)
808      real*8    tau(nl,nl), argumento, deltazdbl
809      real*8    tauinf(nl)
810      real*8    con(nzy), coninf
811      real*8    c1, c2
812      real*8    t1, t2
813      real*8    p1, p2
814      real*8    mr1, mr2
815      real*8    st1, st2
816      real*8    c1box(70), c2box(70)
817      real*8    ff      ! to avoid too small numbers
818      real*8    tvtbs(nzy)
819      real*8    st, beta, ts
820      real*8    zld(nl), zyd(nzy)
821      real*8    correc
822      real*8    deltanudbl
823      integer   isot
824      real*8    yy
825
826c     external function
827      external  we_clean
828      real*8    we_clean
829
830c     formats
831 101  format(i1)
832c***********************************************************************
833
834      ierr = 0
835      varerr = 0.d0
836
837c     some values
838      beta = 1.8d5
839      isot = 1
840      write (ibcode1,101) ib
841      deltanudbl = dble( deltanu(isot,ib) )
842      ff=1.0d10
843      deltazdbl = dble(deltaz)
844
845ccc   
846ccc   
847ccc   
848      do i=1,nl
849         zld(i) = dble(zl(i))
850      enddo
851      do i=1,nzy
852         zyd(i) = dble(zy(i))
853      enddo
854
855      call interhuntdp ( tvtbs,zyd,nzy, v626t1,zld,nl, 1 )
856
857      do i=1,nzy
858         con(i) =  dble( co2y(i) * imr(isot) )
859         correc = 2.d0 * exp( -ee*dble(elow(isot,2))/tvtbs(i) )
860         con(i) = con(i) * ( 1.d0 - correc )
861         mr(i) = dble( co2y(i) / nty(i) )
862      end do
863
864      if ( con(nzy) .le. 0.0d0 ) then
865         ierr = 83
866         varerr = con(nzy)
867         return
868      elseif ( con(nzy-1) .le. con(nzy) ) then
869         write (*,*) ' WARNING in MZTUD121 '
870         write (*,*) '    [CO2] grows with height at CurtisMatrix top.'
871         write (*,*) '    [CO2] @ top = ', con(nzy)
872         coninf = dble( con(nzy) )
873      else
874         coninf = dble( con(nzy) / log( con(nzy-1) / con(nzy) ) )
875      endif
876      call mztf_correccion ( coninf, con, ib )
877
878ccc   
879      call gethist_03 ( ib )
880
881
882c     
883c     tauinf(nl)
884c     
885      call initial
886
887      iaquiZ = nzy - 2
888      iaquiHIST = nhist / 2
889
890      do i=nl,1,-1
891
892         if(i.eq.nl)then
893
894            call intzhunt ( iaquiZ, zl(i),c2,p2,mr2,t2, con)
895            do kr=1,nbox
896               ta(kr)=t2
897            end do
898            call interstrhunt (iaquiHIST, st2,t2,ka,ta)
899            aa = p2 * coninf * mr2 * (st2 * ff)
900            cc = coninf * st2
901            dd = t2 * coninf * st2
902            do kr=1,nbox
903               ccbox(kr) = coninf * ka(kr)
904               ddbox(kr) = t2 * ccbox(kr)
905               c2box(kr) = c2 * ka(kr) * deltazdbl
906            end do
907            c2 = c2 * st2 * deltazdbl
908
909         else
910            call intzhunt ( iaquiZ, zl(i),c1,p1,mr1,t1, con)
911            do kr=1,nbox
912               ta(kr)=t1
913            end do
914            call interstrhunt (iaquiHIST, st1,t1,ka,ta)
915            do kr=1,nbox
916               c1box(kr) = c1 * ka(kr) * deltazdbl
917            end do
918            ! Check interpolation errors :
919            if (c1.le.0.0d0) then
920               ierr=85
921               varerr=c1
922               return
923            elseif (p1.le.0.0d0) then
924               ierr=86
925               varerr=p1
926               return
927            elseif (mr1.le.0.0d0) then
928               ierr=87
929               varerr=mr1
930               return
931            elseif (t1.le.0.0d0) then
932               ierr=88
933               varerr=t1
934               return
935            elseif (st1.le.0.0d0) then
936               ierr=89
937               varerr=st1
938               return
939            endif
940            !
941            c1 = c1 * st1 * deltazdbl
942            aa = aa + ( p1*mr1*(c1*ff) + p2*mr2*(c2*ff)) / 2.d0
943            cc = cc + ( c1 + c2 ) / 2.d0
944            dd = dd + ( t1*c1 + t2*c2 ) / 2.d0
945            do kr=1,nbox
946               ccbox(kr) = ccbox(kr) +
947     @              ( c1box(kr) + c2box(kr) )/2.d0
948               ddbox(kr) = ddbox(kr) +
949     @              ( t1*c1box(kr)+t2*c2box(kr) )/2.d0
950            end do
951
952            mr2 = mr1
953            c2=c1
954            do kr=1,nbox
955               c2box(kr) = c1box(kr)
956            end do
957            t2=t1
958            p2=p1
959         end if
960
961         pp = aa / (cc*ff)
962
963         ts = dd/cc
964         do kr=1,nbox
965            ta(kr) = ddbox(kr) / ccbox(kr)
966         end do
967         call interstrhunt(iaquiHIST, st,ts,ka,ta)
968         call intershphunt(iaquiHIST, alsa,alda,ta)
969
970c     
971
972         eqw = 0.0d0
973         do  kr=1,nbox
974            yy = ccbox(kr) * beta
975            w = we_clean ( yy, pp, alsa(kr),alda(kr) )
976            eqw = eqw + no(kr)*w
977         end do
978
979         argumento = eqw / deltanudbl
980         tauinf(i) = dexp( - argumento )
981
982
983      end do                    ! i continue
984
985
986c     
987c     tau(in,ir) for n<=r
988c     
989
990      iaquiZ = 2
991      do 1 in=1,nl-1
992
993         call initial
994         call intzhunt ( iaquiZ, zl(in), c1,p1,mr1,t1, con)
995         do kr=1,nbox
996            ta(kr) = t1
997         end do
998         call interstrhunt (iaquiHIST, st1,t1,ka,ta)
999         do kr=1,nbox
1000            c1box(kr) = c1 * ka(kr) * deltazdbl
1001         end do
1002         c1 = c1 * st1 * deltazdbl
1003
1004         do 2 ir=in,nl-1
1005
1006            if (ir.eq.in) then
1007               tau(in,ir) = 1.d0
1008               goto 2
1009            end if
1010
1011            call intzhunt ( iaquiZ, zl(ir), c2,p2,mr2,t2, con)
1012            do kr=1,nbox
1013               ta(kr) = t2
1014            end do
1015            call interstrhunt (iaquiHIST, st2,t2,ka,ta)
1016            do kr=1,nbox
1017               c2box(kr) = c2 * ka(kr) * deltazdbl
1018            end do
1019            c2 = c2 * st2 * deltazdbl
1020
1021            aa = aa + ( p1*mr1*(c1*ff) + p2*mr2*(c2*ff)) / 2.d0
1022            cc = cc + ( c1 + c2 ) / 2.d0
1023            dd = dd + ( t1*c1 + t2*c2 ) / 2.d0
1024            do kr=1,nbox
1025               ccbox(kr) = ccbox(kr) +
1026     $              ( c1box(kr) + c2box(kr) ) / 2.d0
1027               ddbox(kr) = ddbox(kr) +
1028     $              ( t1*c1box(kr) + t2*c2box(kr) ) / 2.d0
1029            end do
1030
1031            mr1=mr2
1032            t1=t2
1033            c1=c2
1034            p1=p2
1035            do kr=1,nbox
1036               c1box(kr) = c2box(kr)
1037            end do
1038
1039            pp = aa / (cc * ff)
1040
1041            ts = dd/cc
1042            do kr=1,nbox
1043               ta(kr) = ddbox(kr) / ccbox(kr)
1044            end do
1045            call interstrhunt(iaquiHIST, st,ts,ka,ta)
1046            call intershphunt(iaquiHIST, alsa,alda,ta)
1047
1048c     
1049
1050            eqw = 0.0d0
1051            do kr=1,nbox
1052               yy = ccbox(kr) * beta
1053               w = we_clean ( yy, pp, alsa(kr),alda(kr) )
1054               eqw = eqw + no(kr)*w
1055            end do
1056
1057            argumento = eqw / deltanudbl
1058            tau(in,ir) = dexp( - argumento )
1059
1060 2       continue
1061
1062 1    continue
1063
1064c     
1065c     tau(in,ir) for n>r
1066c     
1067
1068      in=nl
1069
1070      call initial
1071      iaquiZ = nzy - 2
1072      call intzhunt ( iaquiZ, zl(in), c1,p1,mr1,t1, con)
1073      do kr=1,nbox
1074         ta(kr) = t1
1075      end do
1076      call interstrhunt (iaquiHIST, st1,t1,ka,ta)
1077      do kr=1,nbox
1078         c1box(kr) = c1 * ka(kr) * deltazdbl
1079      end do
1080      c1 = c1 * st1 * deltazdbl
1081
1082      do 4 ir=in-1,1,-1
1083
1084         call intzhunt ( iaquiZ, zl(ir), c2,p2,mr2,t2, con)
1085         do kr=1,nbox
1086            ta(kr) = t2
1087         end do
1088         call interstrhunt (iaquiHIST, st2,t2,ka,ta)
1089         do kr=1,nbox
1090            c2box(kr) = c2 * ka(kr) * deltazdbl
1091         end do
1092         c2 = c2 * st2 * deltazdbl
1093
1094         aa = aa + ( p1*mr1*(c1*ff) + p2*mr2*(c2*ff)) / 2.d0
1095         cc = cc + ( c1 + c2 ) / 2.d0
1096         dd = dd + ( t1*c1 + t2*c2 ) / 2.d0
1097         do kr=1,nbox
1098            ccbox(kr) = ccbox(kr) +
1099     $           ( c1box(kr) + c2box(kr) ) / 2.d0
1100            ddbox(kr) = ddbox(kr) +
1101     $           ( t1*c1box(kr) + t2*c2box(kr) ) / 2.d0
1102         end do
1103
1104         mr1=mr2
1105         c1=c2
1106         t1=t2
1107         p1=p2
1108         do kr=1,nbox
1109            c1box(kr) = c2box(kr)
1110         end do
1111
1112         pp = aa / (cc * ff)
1113         ts = dd / cc
1114         do kr=1,nbox
1115            ta(kr) = ddbox(kr) / ccbox(kr)
1116         end do
1117         call interstrhunt (iaquiHIST, st,ts,ka,ta)
1118         call intershphunt (iaquiHIST, alsa,alda,ta)
1119
1120c     
1121         eqw=0.0d0
1122         do kr=1,nbox
1123            yy = ccbox(kr) * beta
1124            w = we_clean ( yy, pp, alsa(kr),alda(kr) )
1125            eqw = eqw + no(kr)*w
1126         end do
1127
1128         argumento = eqw / deltanudbl
1129         tau(in,ir) = dexp( - argumento )
1130
1131 4    continue
1132
1133c     
1134c     
1135c     
1136      do in=nl-1,2,-1
1137         do ir=in-1,1,-1
1138            tau(in,ir) = tau(ir,in)
1139         end do
1140      end do
1141
1142c     
1143      call MZCUD121 ( tauinf,tau, cf, vc, ib )
1144
1145
1146c     end
1147      return
1148      end
1149
1150
1151
1152c     *** Old MZCUD121_dlvr11.f ***
1153
1154c***********************************************************************
1155
1156      subroutine MZCUD121 ( tauinf,tau, c,vc, ib )
1157
1158c***********************************************************************
1159
1160      implicit none
1161
1162#include "nlte_paramdef.h"
1163#include "nlte_commons.h"
1164
1165
1166c     arguments
1167      real*8            c(nl,nl) ! o
1168      real*8            vc(nl)  ! o
1169      real*8            tau(nl,nl) ! i
1170      real*8            tauinf(nl) ! i
1171      integer           ib      ! i
1172
1173c     local variables
1174      integer   i, in, ir, isot
1175      real*8            a(nl,nl), cf(nl,nl), pideltanu, deltazdbl,pi
1176
1177c***********************************************************************
1178
1179      pi=3.141592
1180      isot = 1
1181      pideltanu = pi*dble(deltanu(isot,ib))
1182      deltazdbl = dble(deltaz)
1183c     
1184      do in=1,nl
1185
1186         do ir=1,nl
1187
1188            cf(in,ir) = 0.0d0
1189            c(in,ir) = 0.0d0
1190            a(in,ir) = 0.0d0
1191
1192         end do
1193
1194         vc(in) = 0.0d0
1195
1196      end do
1197
1198
1199c     
1200      do in=1,nl
1201         do ir=1,nl
1202
1203            if (ir.eq.1) then
1204               cf(in,ir) = tau(in,ir) - tau(in,1)
1205            elseif (ir.eq.nl) then
1206               cf(in,ir) = tauinf(in) - tau(in,ir-1)
1207            else
1208               cf(in,ir) = tau(in,ir) - tau(in,ir-1)
1209            end if
1210
1211         end do
1212      end do
1213
1214
1215c     
1216      do in=2,nl-1
1217         do ir=1,nl
1218            if (ir.eq.in+1) a(in,ir) = -1.d0
1219            if (ir.eq.in-1) a(in,ir) = +1.d0
1220            a(in,ir) = a(in,ir) / ( 2.d5*deltazdbl )
1221         end do
1222      end do
1223
1224c     
1225      do in=1,nl
1226         do ir=1,nl
1227            cf(in,ir) = cf(in,ir) * pideltanu
1228         end do
1229      end do
1230
1231
1232      do in=2,nl-1
1233         do ir=1,nl
1234            do i=1,nl
1235               c(in,ir) = c(in,ir) + a(in,i) * cf(i,ir)
1236            end do
1237         end do
1238         vc(in) =  pideltanu /( 2.d5*deltazdbl ) *
1239     @        ( tau(in-1,1) - tau(in+1,1) )
1240      end do
1241
1242c     
1243      do in=2,nl-1
1244         c(in,nl-2) = c(in,nl-2) - c(in,nl)
1245         c(in,nl-1) = c(in,nl-1) + 2.d0*c(in,nl)
1246      end do
1247
1248
1249c     end
1250      return
1251      end
1252
1253
1254
1255c     *** Old MZESC121_dlvr11_03.f ***
1256
1257c***********************************************************************
1258      subroutine MZESC121
1259c***********************************************************************
1260
1261      implicit none
1262
1263#include "nlte_paramdef.h"
1264#include "nlte_commons.h"
1265
1266
1267c     local variables
1268      integer         i,ierr
1269      real*8          factor0200, factor0220, factor1000
1270      real*8          aux_0200(nl), aux2_0200(nl)
1271      real*8          aux_0220(nl), aux2_0220(nl)
1272      real*8          aux_1000(nl), aux2_1000(nl)
1273      real*8          varerr
1274
1275c***********************************************************************
1276
1277!      call zerov (taustar12, nl)
1278      taustar12(1:nl)=0.d0
1279      call zero2v(aux_0200,aux2_0200, nl)
1280      call zero2v(aux_0220,aux2_0220, nl)
1281      call zero2v(aux_1000,aux2_1000, nl)
1282
1283      call MZESC121sub (aux_0200,aux2_0200, 2 , ierr, varerr)
1284      if (ierr .gt. 0) call ERRORS (ierr,varerr)
1285      call MZESC121sub (aux_0220,aux2_0220, 3 , ierr, varerr)
1286      if (ierr .gt. 0) call ERRORS (ierr,varerr)
1287      call MZESC121sub (aux_1000,aux2_1000, 4 , ierr, varerr)
1288      if (ierr .gt. 0) call ERRORS (ierr,varerr)
1289
1290      factor0220 = 1.d0
1291      factor0200 = dble( (nu(1,2)-nu(1,1)) / (nu12_0200-nu(1,1)) )
1292      factor1000 = dble( (nu(1,2)-nu(1,1)) / (nu12_1000-nu(1,1)) )
1293      do i=1,nl
1294         taustar12(i) = taustar12(i)
1295     @        + aux_0200(i) * factor0200
1296     @        + aux_0220(i) * factor0220
1297     @        + aux_1000(i) * factor1000
1298      enddo
1299
1300      call mzescape_normaliz ( taustar12, 2 )
1301
1302c     end
1303      return
1304      end
1305
1306
1307c     *** Old MZESC121sub_dlvr11_03.f ***
1308
1309c***********************************************************************
1310
1311      subroutine MZESC121sub (taustar,tauinf, ib, ierr, varerr )
1312
1313c***********************************************************************
1314
1315      implicit none
1316
1317#include "nlte_paramdef.h"
1318#include "nlte_commons.h"
1319
1320
1321c     arguments
1322      real*8          taustar(nl) ! o
1323      real*8          tauinf(nl) ! o
1324      integer         ib      ! i
1325      integer         ierr        ! o
1326      real*8          varerr      ! o
1327
1328
1329c     local variables and constants
1330      integer         i, iaquiHIST, iaquiZ, isot
1331      real*8          con(nzy), coninf
1332      real*8          c1, c2, ccc
1333      real*8          t1, t2
1334      real*8          p1, p2
1335      real*8          mr1, mr2
1336      real*8          st1, st2
1337      real*8          c1box(70), c2box(70)
1338      real*8          ff      ! to avoid too small numbers
1339      real*8          tvtbs(nzy)
1340      real*8          st, beta, ts
1341      real*8          zld(nl), zyd(nzy)
1342      real*8          correc
1343      real*8          deltanudbl, deltazdbl
1344      real*8          yy
1345
1346c     external function
1347      external        we_clean
1348      real*8          we_clean
1349
1350c     formats
1351 101  format(i1)
1352
1353c***********************************************************************
1354
1355      ierr = 0
1356      varerr = 0.d0
1357
1358      beta = 1.8d5
1359      isot = 1
1360      write ( ibcode1, 101) ib
1361      deltanudbl = dble( deltanu(isot,ib) )
1362      ff=1.0d10
1363      deltazdbl = dble(deltaz)
1364
1365c     
1366      do i=1,nzy
1367         zyd(i) = dble(zy(i))
1368      enddo
1369      do i=1,nl
1370         zld(i) = dble(zl(i))
1371      enddo
1372
1373      call interhuntdp ( tvtbs,zyd,nzy, v626t1,zld,nl, 1 )
1374
1375      do i=1,nzy
1376         con(i) =  dble( co2y(i) * imr(isot) )
1377         correc = 2.d0 * dexp( -ee*dble(elow(isot,2))/tvtbs(i) )
1378         con(i) = con(i) * ( 1.d0 - correc )
1379         mr(i) = dble(co2y(i)/nty(i))
1380      end do
1381      if ( con(nzy) .le. 0.0d0 ) then
1382         ierr = 63
1383         varerr = con(nzy)
1384         return
1385      elseif ( con(nzy-1) .le. con(nzy) ) then
1386         write (*,*) ' WARNING in MZESC121sub '
1387         write (*,*) '    [CO2] grows with height at CurtisMatrix top.'
1388         write (*,*) '    [CO2] @ top = ', con(nzy)
1389         coninf = dble( con(nzy) )
1390      else
1391         coninf = dble( con(nzy) / log( con(nzy-1) / con(nzy) ) )
1392      endif
1393      call mztf_correccion ( coninf, con, ib )
1394
1395c     
1396      call gethist_03 ( ib )
1397
1398c     
1399c     tauinf
1400c     
1401      call initial
1402
1403      iaquiHIST = nhist/2
1404      iaquiZ = nzy - 2
1405
1406      do i=nl,1,-1
1407
1408         if(i.eq.nl)then
1409
1410            call intzhunt (iaquiZ, zl(i),c2,p2,mr2,t2, con)
1411            do kr=1,nbox
1412               ta(kr)=t2
1413            end do
1414            call interstrhunt (iaquiHIST, st2,t2,ka,ta)
1415            ! Check interpolation errors :
1416            if (c2.le.0.0d0) then
1417               ierr=65
1418               varerr=c2
1419               return
1420            elseif (p2.le.0.0d0) then
1421               ierr=66
1422               varerr=p2
1423               return
1424            elseif (mr2.le.0.0d0) then
1425               ierr=67
1426               varerr=mr2
1427               return
1428            elseif (t2.le.0.0d0) then
1429               ierr=68
1430               varerr=t2
1431               return
1432            elseif (st2.le.0.0d0) then
1433               ierr=69
1434               varerr=st2
1435               return
1436            endif
1437            !
1438            aa = p2 * coninf * mr2 * (st2 * ff)
1439            cc = coninf * st2
1440            dd = t2 * coninf * st2
1441            do kr=1,nbox
1442               ccbox(kr) = coninf * ka(kr)
1443               ddbox(kr) = t2 * ccbox(kr)
1444               c2box(kr) = c2 * ka(kr) * deltazdbl
1445            end do
1446            c2 = c2 * st2 * deltazdbl
1447
1448         else
1449            call intzhunt (iaquiZ, zl(i),c1,p1,mr1,t1, con)
1450            do kr=1,nbox
1451               ta(kr)=t1
1452            end do
1453            call interstrhunt (iaquiHIST,st1,t1,ka,ta)
1454            do kr=1,nbox
1455               c1box(kr) = c1 * ka(kr) * deltazdbl
1456            end do
1457            c1 = c1 * st1 * deltazdbl
1458            aa = aa + ( p1*mr1*(c1*ff) + p2*mr2*(c2*ff)) / 2.d0
1459            cc = cc + ( c1 + c2 ) / 2.d0
1460            ccc = ( c1 + c2 ) / 2.d0
1461            dd = dd + ( t1*c1 + t2*c2 ) / 2.d0
1462            do kr=1,nbox
1463               ccbox(kr) = ccbox(kr) +
1464     @              ( c1box(kr) + c2box(kr) )/2.d0
1465               ddbox(kr) = ddbox(kr) +
1466     @              ( t1*c1box(kr)+t2*c2box(kr) )/2.d0
1467            end do
1468
1469            mr2 = mr1
1470            c2=c1
1471            do kr=1,nbox
1472               c2box(kr) = c1box(kr)
1473            end do
1474            t2=t1
1475            p2=p1
1476         end if
1477
1478         pp = aa / (cc*ff)
1479
1480         ts = dd/cc
1481         do kr=1,nbox
1482            ta(kr) = ddbox(kr) / ccbox(kr)
1483         end do
1484         call interstrhunt(iaquiHIST,st,ts,ka,ta)
1485         call intershphunt(iaquiHIST,alsa,alda,ta)
1486
1487c     
1488         eqw=0.0d0
1489         do  kr=1,nbox
1490            yy = ccbox(kr) * beta
1491            w = we_clean ( yy, pp, alsa(kr),alda(kr) )
1492            eqw = eqw + no(kr)*w
1493         end do
1494         tauinf(i) = dexp( - eqw / deltanudbl )
1495         if (tauinf(i).lt.0.d0) tauinf(i) = 0.0d0
1496
1497         if (i.eq.nl) then
1498            taustar(i) = 0.0d0
1499         else
1500            taustar(i) = deltanudbl * (tauinf(i+1)-tauinf(i))
1501     @           / ( beta * ccc  )
1502         endif
1503
1504      end do
1505
1506
1507
1508c     end
1509      return
1510      end
1511
1512
1513c     *** Old MZTVC121_dlvr11.f ***
1514
1515c***********************************************************************
1516
1517      subroutine MZTVC121
1518
1519c***********************************************************************
1520
1521      implicit none
1522
1523!!!!!!!!!!!!!!!!!!!!!!!
1524!     common variables & constants
1525
1526#include "nlte_paramdef.h"
1527#include "nlte_commons.h"
1528
1529
1530!     arguments
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      integer       nmu
1597      parameter     (nmu = 8)
1598      real*8        tau(nl,nl), argumento
1599      real*8        con(nzy), coninf
1600      real*8        c1, c2
1601      real*8        t1, t2
1602      real*8        p1, p2
1603      real*8        mr1, mr2
1604      real*8        st1, st2
1605      real*8        c1box(70), c2box(70)
1606      real*8        ff      ! to avoid too small numbers
1607      real*8        tvtbs(nzy)
1608      real*8        st, beta, ts
1609      real*8        zld(nl), zyd(nzy), deltazdbl
1610      real*8        correc
1611      real*8        deltanudbl, pideltanu,pi
1612      real*8        yy
1613      real*8        minvc, maxtau
1614
1615c     external function
1616      external      we_clean
1617      real*8        we_clean
1618
1619c     formats
1620 101  format(i1)
1621
1622c***********************************************************************
1623
1624      ierr = 0
1625      varerr = 0.d0
1626c     
1627      pi=3.141592
1628      isot = 1
1629      beta = 1.8d5
1630      write (ibcode1,101) ib
1631      deltanudbl = dble( deltanu(isot,ib) )
1632      pideltanu = pi*deltanudbl
1633      ff=1.0d10
1634      deltazdbl = dble(deltaz)
1635c     
1636c     
1637c     
1638
1639      do i=1,nzy
1640         zyd(i) = dble(zy(i))
1641      enddo
1642      do i=1,nl
1643         zld(i) = dble(zl(i))
1644      enddo
1645
1646      call interhuntdp ( tvtbs,zyd,nzy, v626t1,zld,nl, 1 )
1647
1648      do i=1,nzy
1649         con(i) =  dble( co2y(i) * imr(isot) )
1650         correc = 2.d0 * dexp( -ee*dble(elow(isot,2))/tvtbs(i) )
1651         con(i) = con(i) * ( 1.d0 - correc )
1652         mr(i) = dble(co2y(i)/nty(i))
1653      end do
1654
1655      if ( con(nzy) .le. 0.0d0 ) then
1656         ierr = 53
1657         varerr = con(nzy)
1658         return
1659      elseif ( con(nzy-1) .le. con(nzy) ) then
1660         write (*,*) ' WARNING in MZTVC121sub '
1661         write (*,*) '    [CO2] grows with height at CurtisMatrix top.'
1662         write (*,*) '    [CO2] @ top = ', con(nzy)
1663         coninf = dble( con(nzy) )
1664      else
1665         coninf = dble( con(nzy) / log( con(nzy-1) / con(nzy) ) )
1666      endif
1667      call mztf_correccion ( coninf, con, ib )
1668
1669ccc   
1670      call gethist_03 ( ib )
1671
1672c     
1673c     tau(1,ir)
1674c     
1675      call initial
1676
1677      iaquiHIST = nhist/2
1678
1679      in=1
1680
1681      tau(in,1) = 1.d0
1682
1683      call initial
1684      iaquiZ = 2
1685      call intzhunt ( iaquiZ, zl(in), c1,p1,mr1,t1, con)
1686      do kr=1,nbox
1687         ta(kr) = t1
1688      end do
1689      call interstrhunt (iaquiHIST, st1,t1,ka,ta)
1690      do kr=1,nbox
1691         c1box(kr) = c1 * ka(kr) * deltazdbl
1692      end do
1693      c1 = c1 * st1 * deltazdbl
1694                                ! Check interpolation errors :
1695      if (c1.le.0.0d0) then
1696         ierr=55
1697         varerr=c1
1698         return
1699      elseif (p1.le.0.0d0) then
1700         ierr=56
1701         varerr=p1
1702         return
1703      elseif (mr1.le.0.0d0) then
1704         ierr=57
1705         varerr=mr1
1706         return
1707      elseif (t1.le.0.0d0) then
1708         ierr=58
1709         varerr=t1
1710         return
1711      elseif (st1.le.0.0d0) then
1712         ierr=59
1713         varerr=st1
1714         return
1715      endif
1716                                !
1717
1718      do 2 ir=2,nl
1719
1720         call intzhunt (iaquiZ, zl(ir), c2,p2,mr2,t2, con)
1721         do kr=1,nbox
1722            ta(kr) = t2
1723         end do
1724         call interstrhunt (iaquiHIST, st2,t2,ka,ta)
1725         do kr=1,nbox
1726            c2box(kr) = c2 * ka(kr) * deltazdbl
1727         end do
1728         c2 = c2 * st2 * deltazdbl
1729
1730         aa = aa + ( p1*mr1*(c1*ff) + p2*mr2*(c2*ff)) / 2.d0
1731         cc = cc + ( c1 + c2 ) / 2.d0
1732         dd = dd + ( t1*c1 + t2*c2 ) / 2.d0
1733         do kr=1,nbox
1734            ccbox(kr) = ccbox(kr) + ( c1box(kr) + c2box(kr) ) /2.d0
1735            ddbox(kr) = ddbox(kr) +
1736     $           ( t1*c1box(kr) + t2*c2box(kr) ) / 2.d0
1737         end do
1738
1739         mr1=mr2
1740         t1=t2
1741         c1=c2
1742         p1=p2
1743         do kr=1,nbox
1744            c1box(kr) = c2box(kr)
1745         end do
1746
1747         pp = aa / (cc * ff)
1748
1749         ts = dd/cc
1750         do kr=1,nbox
1751            ta(kr) = ddbox(kr) / ccbox(kr)
1752         end do
1753         call interstrhunt(iaquiHIST, st,ts,ka,ta)
1754         call intershphunt(iaquiHIST, alsa,alda,ta)
1755
1756         eqw=0.0d0
1757         do kr=1,nbox
1758            yy = ccbox(kr) * beta
1759            w = we_clean ( yy, pp, alsa(kr),alda(kr) )
1760            eqw = eqw + no(kr)*w
1761         end do
1762
1763         argumento = eqw / deltanudbl
1764         tau(in,ir) = dexp( - argumento )
1765
1766 2    continue
1767
1768
1769c     
1770c     
1771c     
1772      do in=nl,2,-1
1773         tau(in,1) = tau(1,in)
1774      end do
1775
1776c     
1777      vc(1) = 0.0d0
1778      vc(nl) = 0.0d0
1779      do in=2,nl-1
1780         vc(in) =  pideltanu /( 2.d5*deltazdbl ) *
1781     @        ( tau(in-1,1) - tau(in+1,1) )
1782         if (vc(in) .lt. 0.0d0) vc(in) = vc(in-1)
1783      end do
1784
1785c     
1786c     Tracking potential numerical errors
1787c     
1788      minvc = 1.d6
1789      maxtau = tau(nl,1)
1790      do in=2,nl-1
1791         minvc = min( minvc, vc(in) )
1792         maxtau = max( maxtau, tau(in,1) )
1793      end do
1794      if (maxtau .gt. 1.0d0) then
1795         ierr = 52
1796         varerr = maxtau
1797         return
1798      else if (minvc .lt. 0.0d0) then
1799         ierr = 51
1800         varerr = minvc
1801         return
1802      endif
1803
1804c     end
1805      return
1806      end
1807
1808
1809
1810
1811
1812
1813
1814
1815
Note: See TracBrowser for help on using the repository browser.