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

Last change on this file since 2639 was 2616, checked in by romain.vande, 4 years ago

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