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

Last change on this file since 3118 was 3018, checked in by emillour, 16 months ago

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