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

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

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