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

Last change on this file since 2266 was 2151, checked in by aslmd, 5 years ago

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