c     program bda/progf/iso.f
c
c     calcul des isochrones
c     avec deviations evolutives et gravites
c     avec deviations dans d/b2-v1  et dans  c1/beta
c
      dimension dm(60),dlm(60),dl(60,60),dt(60,60),da(60,60)
      dimension basl(60),bast(60),ym(2)
      dimension kl(2),dlog(60,60),va(60),yl(2),yt(2),xmas(60,60)
      dimension dmas(60,60),lex(60),dda(60,60),ddl(60,60),ddt(60,60)
      dimension ddmas(60,60),gg(60,60),zlage(60)
      dimension glog(5000),xm(5000),xlm(5000)
      dimension xdes(1000),ydes(1000)
      character stat*3,cage*4,cagf*5,filnam*45
      character titf*30,agedes*6,atype*1,path*35,paty(6)*4
      real masset
      real indix,indix7,indix8,indix9
      logical lx
      data xmax,xmin,ymax,ymin/-100.,100.,-100.,100./
c*****
      data paty/'theo','isu_','iso_','isq_','is2_','isg_'/
      data agl,agt,paglog/0.15,0.05,2.5/
      data xminl,xmint/0.05,0.01/
      data dseg/0.01/
      data iprim/0/
c*****    
c
      call getarg(1,atype)
      call getarg(2,agedes)
      read(atype, 504) itype
      read(agedes,505) aglog
  504 format(i1)
  505 format(f5.2)
c
      do 76 i=1,60
      lex(i)=6
   76 continue
c
      call getenv('BDA',path)
      do 74 m=1,35
      mm=36-m
      if(path(mm:mm).ne.' ') goto 75
   74 continue
   75 open(10,file=path(1:mm)//'/sequence/modele.dat',status='old')
c
c     lecture des donnes des modeles
c
      read(10,400) titf
  400 format(a30)
      read(10,402) nma,jmp
  402 format(i3,1x,i3)
      read(10,401)
  401 format(1x)
      do 77 ki=1,nma
      read(10,403) masset
  403 format(f7.3)
      dm(ki)=masset
      read(10,401)
      do 78 km=1,jmp
      read(10,404) ddl(km,ki),ddt(km,ki),dda(km,ki),ddmas(km,ki)
  404 format(f7.3,f6.3,1x,f10.4,f8.3)
   78 continue
      do 79 i=1,5
      read(10,401)
   79 continue
   77 continue
c
      do 80 m=1,nma
      dlm(m)=log10(dm(m))
      pex=lex(m)
      pex=10.**pex
      do 1 j=1,jmp
      dl(m,j)=ddl(j,m)
      dt(m,j)=ddt(j,m)
      xmas(m,j)=ddmas(j,m)
      da(m,j)=dda(j,m)*pex
      dmas(m,j)=log10(xmas(m,j))
      r2=3.1472+33.583+ddl(j,m)-4.*ddt(j,m)
      gg(m,j)=-7.176+33.299+log10(ddmas(j,m))-r2
      if(j.ne.1) goto 1
      basl(m)=dl(m,j)
      bast(m)=dt(m,j)
    1 dlog(m,j)=log10(da(m,j))
c
   80 continue
c
      if(aglog.gt.9.5) then
      agt=0.003
      xmint=0.001
      endif
c
c     choix de la masse superieure
c
      do 600 m=1,nma 
      zlage(m)=dlog(m,jmp)
  600 continue
      deumas=10**(fipoi(aglog,nma,zlage,dlm))
c
c     choix de la masse inferieure
c
      xaglog=aglog+paglog
      premas=10**(fipoi(xaglog,nma,zlage,dlm))
  733 xmpre=premas
      xmvpre=0.
      bmvpre=0.
      age=10.**aglog
c
c     ouverture des fichiers de sortie
c
      if(aglog.lt.10.0) then
        write(cage,500) aglog
  500   format(f4.2)
        filnam=path(1:mm)//'/isochrone/'//paty(itype)//cage
      else
        write(cagf,501) aglog
  501   format(f5.2)
        filnam=path(1:mm)//'/isochrone/'//paty(itype)//cagf
      endif
      inquire(file=filnam,exist=lx)
      if(.not.lx) then
        stat='new'
        else
        stat='old'
      endif
      open(15,file=filnam,status=stat)
c
c     calculs
c    
      i=1
      delm=0.
      xm2w=premas
   65 xm2w=xm2w+delm
      xmm=xm2w  
      if(delm.lt.0) print*,' delm.lt.0'
      if(xmm.gt.deumas) go to 40
      xlm(i)=log10(xmm)
      xm(i)=xmm     
      m=indice(xmm,dm,nma)
      do 8 k=1,jmp
      va(k)=dlog(m,k)+(dlog(m+1,k)-dlog(m,k))*(xlm(i)-dlm(m))/
     1(dlm(m+1)-dlm(m))
      if(aglog-va(1))11,11,12
   11 yt(1)=dt(m,1)
      yt(2)=dt(m+1,1)
      yl(1)=dl(m,1)
      yl(2)=dl(m+1,1)
      ym(1)=dmas(m,1)
      ym(2)=dmas(m+1,1)
      kl(1)=1
      kl(2)=1
      goto 18
   12 if(k-jmp)14,13,13
   13 if(aglog-va(jmp))14,16,40
   16 diff=1.
      kl(1)=jmp-1
      kl(2)=jmp
      goto 15
   14 if(aglog-va(k))19,8,8
   19 vak=10.**va(k)
      vakm1=10.**va(k-1)
      diff=(age-vakm1)/(vak-vakm1)
      kl(1)=k-1
      kl(2)=k
      goto 15
    8 continue
   15 do 17 l=1,2
      ll=m+l-1
      ym(l)=xmas(ll,k-1)+((xmas(ll,k)-xmas(ll,k-1))*diff)
      ym(l)=log10(ym(l))
      yt(l)=dt(ll,k-1)+((dt(ll,k)-dt(ll,k-1))*diff)
   17 yl(l)=dl(ll,k-1)+((dl(ll,k)-dl(ll,k-1))*diff)
   18 ff=(xlm(i)-dlm(m))/(dlm(m+1)-dlm(m))
      xl=yl(1)+ (yl(2)-yl(1))*ff
      xt=yt(1)+ (yt(2)-yt(1))*ff
      xmlog=ym(1)+(ym(2)-ym(1))*ff
      xmm=10.**xmlog
      glog(i)=-10.6113+xmlog+4.*xt-xl
c
c***********
c
      if(aglog.gt.8.8) then
        if(kl(2).le.13) then
        call bvmv(xl,xt,xmv,indix,itype)
        goto 51
      endif
       if(kl(2).eq.14) then
         xt7=dt(m,13)+(dt(m+1,13)-dt(m,13))*ff
         xl7=dl(m,13)+(dl(m+1,13)-dl(m,13))*ff
         xt8=dt(m,14)+(dt(m+1,14)-dt(m,14))*ff
         xl8=dl(m,14)+(dl(m+1,14)-dl(m,14))*ff
         call bvmv(xl7,xt7,xmv7,indix7,itype)
         call bvmiii(xl8,xt8,xmv8,indix8,itype)
         xmv=xmv7+(xmv8-xmv7)*diff
         indix=indix7+(indix8-indix7)*diff
         go to 51
       endif
        if(kl(2).gt.14) then
         call bvmiii(xl,xt,xmv,indix,itype)
         go to 51
        endif
      endif
c
c**********
c
      if(aglog.le.8.8) then
        if(kl(2).le.13) then
        call bvmv(xl,xt,xmv,indix,itype)
        goto 51
      endif
       if(kl(2).eq.14) then
         xt7=dt(m,13)+(dt(m+1,13)-dt(m,13))*ff
         xl7=dl(m,13)+(dl(m+1,13)-dl(m,13))*ff
         xt8=dt(m,14)+(dt(m+1,14)-dt(m,14))*ff
         xl8=dl(m,14)+(dl(m+1,14)-dl(m,14))*ff
         call bvmv(xl7,xt7,xmv7,indix7,itype)
         call bvmiii(xl8,xt8,xmv8,indix8,itype)
         xmv=xmv7+(xmv8-xmv7)*diff
         indix=indix7+(indix8-indix7)*diff
         go to 51
       endif
        if(kl(2).lt.20) then
         call bvmiii(xl,xt,xmv,indix,itype)
         go to 51
        endif
       if(kl(2).eq.20) then
         xt8=dt(m,19)+(dt(m+1,19)-dt(m,19))*ff
         xl8=dl(m,19)+(dl(m+1,19)-dl(m,19))*ff
         xt9=dt(m,20)+(dt(m+1,20)-dt(m,20))*ff
         xl9=dl(m,20)+(dl(m+1,20)-dl(m,20))*ff
         call bvmiii(xl8,xt8,xmv8,indix8,itype)
         call bvmi(xl9,xt9,xmv9,indix9,itype)
         xmv=xmv8+(xmv9-xmv8)*diff
         indix=indix8+(indix9-indix8)*diff
         go to 51
       endif
        if(kl(2).gt.20) then
        call bvmi(xl,xt,xmv,indix,itype)
        go to 51
        endif 
      endif
c**********
   51 xmbol=-2.5*xl+4.75
      base=fipoi(xt,nma,bast,basl)
      dev=2.5*(xl-base)
      x2logr=3.1472+33.583+xl-4.*xt
      xlogr=0.5*x2logr
      rsol=(10.**xlogr)/6.9599e+10
      xlogg=-7.176+33.299+xmlog-x2logr
      write(26,800) xt,xlogg
  800 format(1x,f8.3, 3x,f8.3)
c
c     variables pour le dessin
c
      if(itype.eq.1) then
        xdes(i)=xt
        ydes(i)=xmbol
      else
        ydes(i)=xmv
        xdes(i)=indix
      endif
c*****
      n1des=i
      nombre=i
c*****
      if(i.eq.1) then
       delm=(dm(m)-dm(m+1))/10.
      else
       distl=abs(ydes(i)-ydes(i-1))
       distt=abs(xdes(i)-xdes(i-1))
       if(distl.gt.agl.or.distt.gt.agt) then
        xm2w=xm2w-delm
        delm=delm/1.2
        go to 65
       endif
       if(distl.lt.xminl.and.distt.lt.xmint) then
        xm2w=xm2w-delm
        delm=10.*delm
        if((xm2w+delm).gt.deumas) delm=2.*delm/10. 
        go to 65
       endif
      endif
      i=i+1
      go to 65
c***** 
    7 continue
   40 continue
c*****
      do 601 jj=1,nombre
        write(15,118) ydes(jj),xdes(jj)
  601 continue
      close(15,status='keep')
c
  730 stop 
c
  100 format(e11.4,1x,f6.2,1x,f6.2,1x,i5,i3)
  101 format(f6.2,2x,i3)
  102 format(i2,1x,f6.4,2f7.3,f8.3)
c 103 format(1x,i4,2f10.3,2f8.3,1x,2f8.3,2x,2f9.3,2x,i3,' -',i3,
c    13x,f8.3,3x,f8.3,1x,e11.4)
  104 format(1x,'1')
  105 format(1x,' age',2x,e14.7,2x,' nombre',2x,i5,2x,' masse no. 1',2x,
     1f7.3,2x,' masse no. 2 ',f7.3,//)
  106 format(1x,///)
  107 format(8x,' Masse',12x,' log L',3x,' log Te',3x,' M Bol',13x,
     1' log g',4x,' r',8x,' dev'/)
c    1 ' Teff',8x,' xp',6x,' U-B',15x,' Mbol',7x,' dev',3x,' densli'/)
  108 format(1x,i4,f10.3,4f9.3,2i5,3f9.3)
  109 format(1x,2i4)
  110 format(1x,i4,f10.3,7f9.3)
  111 format(1x,' log age ',f8.3,//)
  114 format(1x,' bin ', 9x,f10.3,f8.3)
  118 format(2x,f7.3,2x,f7.3)
  209 format(1h1,' Donnees des traces evolutifs'/2i4)
  301 format(/,1x,f8.3,/)
  302 format(f8.3,f6.3,e15.7,2f8.3)
  303 format(1h1)
c
      end
c
c
      subroutine sortw ( a,n,j )
c
c
c     *******************************************************
c
c          sous - programme de mise en ordre de nombres
c
c     *******************************************************
c
c
c     ce sous - programme est une traduction fortran
c     de l algorithme 271  * quickersort * de l acm
c
c          paul bartholdi        -        mai 1969
c                  observatoire de geneve
c
c        ******************************************
c
      dimension a(j),n(j),lt(20),ut(20)
      integer p,q,ut
c
      jj = j
      i=1
      m=1
    1 if ( j-i-1 ) 12,12,2
    2 p = (i+j)/2
      t = a(p)
      a(p) = a(i)
      it = n(p)
      n(p) = n(i)
      q = j
      ip = i+1
      do 3 k=ip,q
      if ( a(k) - t ) 3,3,4
    4 if ( q-k ) 20,19,19
   19 if ( a(q) - t ) 6,5,5
    6 x = a(k)
      a(k) = a(q)
      a(q) = x
      is = n(k)
      n(k) = n(q)
      n(q) = is
      q = q-1
      go to 3
    5 q = q-1
      go to 4
   20 q = k-1
      go to 8
    3 continue
    8 a(i) = a(q)
      a(q) = t
      n(i) = n(q)
      n(q) = it
      if ( 2*q-i-j ) 10,10,9
    9 lt(m) = i
      ut(m) = q-1
      i = q+1
      go to 11
   10 lt(m) = q+1
      ut(m) = j
      j = q-1
   11 m = m+1
      go to 1
   12 if ( i-j ) 17,14,14
   17 if ( a(i) - a(j) ) 14,14,13
   13 x = a(i)
      a(i) = a(j)
      a(j) = x
      is = n(i)
      n(i) = n(j)
      n(j) = is
   14 m = m-1
      if ( m ) 16,16,15
   15 i = lt(m)
      j = ut(m)
      go to 1
   16 j = jj
      return
      end
c
      function fipoi(x,n,a,b)
c
c     utilise le sous-programme search.
c     sous-programme d interpolation.
c
      dimension a(n),b(n)
      k=indice(x,a,n)
c
c     interpolation lineaire.
c
      fipoi=b(k)+(b(k+1)-b(k))*(x-a(k))/(a(k+1)-a(k))
      return
      end
c
c
      function indice(x0,x,m)
c
c     recherche rapide de la position d une valeur x0 dans une table
c     monotone croissante ou decroissante de m nombres x(i)
c
c     si  k = indice(x0,x,m)  on aura   x0 compris entre x(k) et x(k+1)
c             ou x0 = x(k).
c     si x0 exterieur a la table indice=1 si x0 du cote de x(1)
c                                indice = m-1 si x0 du cote de x(m)
c
      dimension x(m)
c
      n=m
      k=1
c
    1 if(n-k-1) 2,5,2
    2 i=(k+n)/2
      if((x(i)-x0)*(x(n)-x(1))) 4,4,3
    3 n=i
      go to 1
    4 k=i
      go to 1
    5 indice = k
c
      return
c
      end 
c
c
c
      function ran1(idum)
c
      dimension r(97)
      parameter (m1=259200,ia1=7141,ic1=54773,rm1=3.8580247e-6)
      parameter (m2=134456,ia2=8121,ic2=28411,rm2=7.4373773e-6)
      parameter (m3=243000,ia3=4561,ic3=51349)
      data iff /0/
c
      if (idum.lt.0.or.iff.eq.0) then
        iff=1
        ix1=mod(ic1-idum,m1)
        ix1=mod(ia1*ix1+ic1,m1)
        ix2=mod(ix1,m2)
        ix1=mod(ia1*ix1+ic1,m1)
        ix3=mod(ix1,m3)
        do 11 j=1,97
          ix1=mod(ia1*ix1+ic1,m1)
          ix2=mod(ia2*ix2+ic2,m2)
          r(j)=(float(ix1)+float(ix2)*rm2)*rm1
11      continue
        idum=1
      endif
      ix1=mod(ia1*ix1+ic1,m1)
      ix2=mod(ia2*ix2+ic2,m2)
      ix3=mod(ia3*ix3+ic3,m3)
      j=1+(97*ix3)/m3
      if(j.gt.97.or.j.lt.1)pause
      ran1=r(j)
      r(j)=(float(ix1)+float(ix2)*rm2)*rm1
      return
      end
c
c
c
C
C
      SUBROUTINE BVMV(WL,WT,WMV,WIND,ITYPE)
C
C     RELATION Teff vs B-V D'APRES ARRIBAS AND MARTINEZ ROGER:1988,
C     A&A,206,63,POUR Teff DANS L'INTERVALLE [5000,8000],ET D'APRES
C     BOHM-VITENSE:1981,ANN. REV. A&A,19,295
C     CORR. BOLO. D'APRES MALAGNINI ET AL:1986,162,140
C     RELATION UBV SCHMIDT-KALER:1982,LANDOLT & BORNSTEIN
C     verifiee 17-4-89
C  
      DIMENSION ZT(28),ZBC(28),ZBV(28),ZUB(28),ZB2V1(28),ZBVG(28)
C
      DATA ZT/ 3.740,3.778,3.813,3.845,3.875,3.903,3.929,3.954,
     1         3.978,4.000,4.041,4.079,4.114,4.146,4.176,4.204,
     2         4.230,4.255,4.301,4.342,4.352,4.380,4.398,4.415,
     3         4.447,4.477,4.653,4.688/   
      DATA ZBC/-0.21,-0.08, 0.00, 0.04, 0.04, 0.02,-0.02,-0.09,
     1         -0.16,-0.25,-0.45,-0.66,-0.87,-1.07,-1.26,-1.43,
     2         -1.58,-1.72,-1.96,-2.17,-2.22,-2.37,-2.46,-2.56,
     3         -2.75,-2.91,-4.00,-4.30/
      DATA ZBV/ 0.70, 0.55, 0.420, 0.32, 0.22, 0.14, 0.10, 0.05,
     1          0.00,-0.03,-0.075,-0.11,-0.13,-0.15,-0.17,-0.18,
     2         -0.19,-0.20,-0.220,-0.24,-0.24,-0.26,-0.26,-0.27,
     3         -0.29,-0.30,-0.360,-0.37/
      DATA ZUB/ 0.23, 0.04,-0.016, 0.02, 0.10, 0.10, 0.09, 0.05,
     1          0.01,-0.06,-0.218,-0.34,-0.43,-0.50,-0.58,-0.62,
     2         -0.67,-0.71,-0.775,-0.84,-0.87,-0.92,-0.95,-0.98,
     3         -1.05,-1.08,-1.300,-1.34/             
      DATA ZB2V1 /+0.426,+0.321,+0.203,+0.118,+0.033,-0.044,-0.073,
     1            -0.108,-0.144,-0.165,-0.197,-0.222,-0.237,-0.251,
     2            -0.265,-0.272,-0.279,-0.286,-0.301,-0.315,-0.315,
     3            -0.329,-0.329,-0.336,-0.350,-0.358,-0.400,-0.407/
      DATA ZBVG  /-0.097,-0.286,-0.451,-0.569,-0.687,-0.781,-0.828,
     1            -0.890,-0.951,-0.987,-1.042,-1.084,-1.109,-1.133,
     2            -1.157,-1.169,-1.181,-1.193,-1.218,-1.242,-1.242,
     3            -1.266,-1.266,-1.278,-1.302,-1.315,-1.387,-1.399/
C
      CB=FIPOI(WT,28,ZT,ZBC)
      WMV=-2.5*WL+4.75-CB
C
      GOTO(1,2,3,4,5,6), ITYPE
    1 GOTO 10
    2 WIND=FIPOI(WT,28,ZT,ZUB)
      GOTO 10
    3 WIND=FIPOI(WT,28,ZT,ZBV)
      GOTO 10
    4 WUB=FIPOI(WT,28,ZT,ZUB)
      WBV=FIPOI(WT,28,ZT,ZBV)
      WIND=WUB-0.72*WBV
      GOTO 10
    5 WIND=FIPOI(WT,28,ZT,ZB2V1)
      GOTO 10
    6 WIND=FIPOI(WT,28,ZT,ZBVG)
C
   10 RETURN
      END                
C
C
C
      SUBROUTINE BVMIII(WL,WT,WMV,WIND,ITYPE)
C
C     RELATION Teff vs B-V ET
C     CORR BOL D'APRES FLOWER A A 54,31 (1977)
C     UBV (CLASSE III) SCHMIDT KALER,1982,LANDOLT ET B.
C
      DIMENSION ZT(38),ZBC(38),ZBV(38),ZUB(38),ZB2V1(38),ZBVG(38)
C
      DATA ZT/ 3.544,3.574,3.592,3.622,3.650,3.674,3.698,3.706,
     1         3.724,3.747,3.769,3.793,3.813,3.845,3.892,3.914,
     2         3.930,3.944,3.954,3.981,4.000,4.023,4.057,4.099,
     3         4.121,4.211,4.234,4.256,4.301,4.324,4.369,4.414,
     4         4.437,4.482,4.504,4.527,4.549,4.594/
      DATA ZBC/-1.66,-1.19,-0.92,-0.66,-0.49,-0.37,-0.26,-0.22,
     1         -0.16,-0.10,-0.06,-0.03,-0.01, 0.01, 0.00,-0.02,
     1         -0.03,-0.06,-0.09,-0.19,-0.24,-0.32,-0.46,-0.72,
     2         -0.84,-1.34,-1.46,-1.58,-1.84,-1.97,-2.22,-2.48,
     3         -2.60,-2.87,-2.99,-3.12,-3.25,-3.50/
      DATA ZBV/ 1.61, 1.54, 1.45, 1.30, 1.160, 1.04, 0.92, 0.88,
     1          0.80, 0.70, 0.60, 0.50, 0.430, 0.33, 0.20, 0.14,
     1          0.10, 0.07, 0.05, 0.00,-0.025,-0.05,-0.08,-0.11,
     2         -0.12,-0.16,-0.17,-0.18,-0.200,-0.21,-0.23,-0.25,
     3         -0.26,-0.28,-0.29,-0.30,-0.310,-0.33/
      DATA ZUB/ 1.88, 1.84, 1.72, 1.44, 1.160, 0.94, 0.66, 0.59,
     1          0.45, 0.28, 0.16, 0.10, 0.090, 0.08, 0.11, 0.11,
     2          0.10, 0.09, 0.06, 0.03,-0.053,-0.13,-0.24,-0.37,
     3         -0.40,-0.54,-0.58,-0.63,-0.740,-0.78,-0.87,-0.94,
     4         -0.97,-1.04,-1.08,-1.10,-1.120,-1.22/             
      DATA ZB2V1/+1.216,+1.148,+1.061,+0.926,+0.808,+0.707,+0.614,
     1           +0.582,+0.520,+0.442,+0.364,+0.273,+0.213,+0.126,
     2           +0.014,-0.042,-0.070,-0.090,-0.104,-0.139,-0.156,
     3           -0.174,-0.194,-0.215,-0.222,-0.250,-0.257,-0.264,
     4           -0.278,-0.285,-0.298,-0.312,-0.319,-0.333,-0.340,
     5           -0.347,-0.354,-0.368/
      DATA ZBVG /+1.075,+0.992,+0.886,+0.690,+0.506,+0.349,+0.191,
     1           +0.138,+0.033,-0.099,-0.230,-0.350,-0.435,-0.557,
     2           -0.715,-0.788,-0.836,-0.877,-0.901,-0.961,-0.991,
     3           -1.021,-1.057,-1.093,-1.105,-1.154,-1.166,-1.178,
     4           -1.202,-1.214,-1.238,-1.262,-1.274,-1.298,-1.310,
     5           -1.322,-1.334,-1.358/
C
      CB=FIPOI(WT,38,ZT,ZBC)
      WMV=-2.5*WL+4.75-CB
C
      GOTO(1,2,3,4,5,6), ITYPE
    1 GOTO 10
    2 WIND=FIPOI(WT,38,ZT,ZUB)
      GOTO 10
    3 WIND=FIPOI(WT,38,ZT,ZBV)
      GOTO 10
    4 WUB=FIPOI(WT,38,ZT,ZUB)
      WBV=FIPOI(WT,38,ZT,ZBV)
      WIND=WUB-0.72*WBV
      GOTO 10
    5 WIND=FIPOI(WT,38,ZT,ZB2V1)
      GOTO 10
    6 WIND=FIPOI(WT,38,ZT,ZBVG)
C
   10 RETURN
      END
C
      SUBROUTINE BVMI(WL,WT,WMV,WIND,ITYPE) 
C
C     RELATION Teff vs B-V ET
C     CORR BOL D'APRES FLOWER A A 54,31 (1977)
C     UBV (CLASSE Ia) SCHMIDT KALER,1982,LANDOLT ET B.
C
      DIMENSION ZT(38),ZBC(38),ZBV(38),ZUB(38),ZB2V1(38),ZBVG(38)
C
      DATA ZT/ 3.447,3.477,3.512,3.544,3.574,3.588,3.605,3.632,
     1         3.656,3.677,3.698,3.705,3.743,3.766,3.793,3.845,
     1         3.892,3.914,3.930,3.944,3.960,4.011,4.038,4.081,
     2         4.137,4.194,4.213,4.288,4.307,4.325,4.363,4.419,
     3         4.457,4.476,4.513,4.532,4.570,4.607/
      DATA ZBC/-3.36,-2.50,-1.72,-1.43,-1.00,-0.84,-0.67,-0.46,
     1         -0.35,-0.22,-0.14,-0.12,-0.01, 0.04, 0.08, 0.13,
     1          0.14, 0.09, 0.00,-0.10,-0.17,-0.38,-0.51,-0.64, 
     2         -0.82,-1.05,-1.16,-1.56,-1.67,-1.76,-2.02,-2.40,
     3         -2.68,-2.79,-3.06,-3.20,-3.46,-3.73/
      DATA ZBV/ 1.80, 1.76, 1.72, 1.70, 1.61, 1.54, 1.450, 1.30,
     1          1.16, 1.04, 0.92, 0.88, 0.70, 0.60, 0.500, 0.33,
     1          0.20, 0.14, 0.10, 0.07, 0.05, 0.00,-0.025,-0.05,  
     2         -0.08,-0.11,-0.12,-0.16,-0.17,-0.18,-0.200,-0.23, 
     3         -0.25,-0.26,-0.28,-0.29,-0.31,-0.33/
      DATA ZUB/ 2.11, 2.05, 1.99, 1.92, 1.81, 1.71, 1.572, 1.25,
     1          1.06, 0.84, 0.70, 0.66, 0.50, 0.45, 0.392, 0.28,
     2          0.17, 0.12,-0.04,-0.15,-0.20,-0.58,-0.635,-0.70,
     3         -0.76,-0.83,-0.85,-0.96,-0.97,-0.99,-1.013,-1.05,
     4         -1.09,-1.11,-1.13,-1.13,-1.16,-1.19/    
      DATA ZB2V1/ 1.435, 1.394, 1.355, 1.335, 1.246, 1.178, 1.092,
     1            0.951, 0.823, 0.715, 0.610, 0.576, 0.423, 0.341,
     2            0.260, 0.126, 0.027,-0.018,-0.048,-0.070,-0.085,
     3           -0.121,-0.139,-0.157,-0.178,-0.199,-0.207,-0.234,
     4           -0.241,-0.248,-0.262,-0.283,-0.297,-0.303,-0.317,
     5           -0.324,-0.337,-0.351/
      DATA ZBVG /+1.335,+1.282,+1.229,+1.202,+1.083,+0.991,+0.872,
     1           +0.674,+0.489,+0.330,+0.171,+0.118,-0.104,-0.228,
     2           -0.352,-0.564,-0.726,-0.801,-0.850,-0.888,-0.913,
     3           -0.975,-1.006,-1.037,-1.075,-1.112,-1.124,-1.174,
     4           -1.187,-1.199,-1.224,-1.261,-1.286,-1.299,-1.324,
     5           -1.336,-1.361,-1.386/
C
      CB=FIPOI(WT,38,ZT,ZBC)
      WMV=-2.5*WL+4.75-CB
C
      GOTO(1,2,3,4,5,6), ITYPE
    1 GOTO 10
    2 WIND=FIPOI(WT,38,ZT,ZUB)
      GOTO 10
    3 WIND=FIPOI(WT,38,ZT,ZBV)
      GOTO 10
    4 WUB=FIPOI(WT,38,ZT,ZUB)
      WBV=FIPOI(WT,38,ZT,ZBV)
      WIND=UB-0.72*WBV
      GOTO 10
    5 WIND=FIPOI(WT,38,ZT,ZB2V1)
      GOTO 10
    6 WIND=FIPOI(WT,38,ZT,ZBVG)
C
   10 RETURN
      END
