SUBROUTINE VISINI
C     =================
C
C     Auxiliary procedure for RESOLV
C     initialization of necessary quantities for treating the viscosity
C     in disks
C
      INCLUDE 'IMPLIC.FOR'
      INCLUDE 'BASICS.FOR'
      INCLUDE 'ATOMIC.FOR'
      INCLUDE 'MODELQ.FOR'
      INCLUDE 'ITERAT.FOR'
C
      AMUV0=DMVISC**(ZETA0+UN)
      AMUV1=UN-AMUV0
      GP=0.
      GN=UN
      IF(INMP.GT.0) THEN
         GP=UN
         GN=0.
      END IF
C
      IF(IVISC.LE.1) THEN 
         x=0.
         DO ID=1,ND
            DMD=DM(1)
            IF(ID.GT.1) DMD=(DM(ID)+DM(ID-1))*HALF
            IF(DM(ID).LE.DMVISC*DM(ND)) THEN
               VISCD(ID)=(UN-FRACTV)*(ZETA1+UN)/
     *                   DMVISC**(ZETA1+UN)*(DM(ID)/DM(ND))**ZETA1
               THETAV(ID)=(UN-FRACTV)*(DMD/DMVISC/DM(ND))**(ZETA1+UN)
             ELSE
               VISCD(ID)=FRACTV*(ZETA0+UN)/AMUV1*
     *                   (DM(ID)/DM(ND))**ZETA0
               THETAV(ID)=(UN-FRACTV)+FRACTV*((DMD/DM(ND))**(ZETA0+UN)-
     *                    AMUV0)/AMUV1
             END IF
             TVISC(ID)=EDISC*VISCD(ID)*DENS(ID)
             DTVIST(ID)=0.
             DTVISR(ID)=EDISC*VISCD(ID)*WMM(ID)
             DTVISN(ID)=0.
             if(id.gt.1) then
             X=X+HALF*(TVISC(ID)/DENS(ID)+TVISC(ID-1)/DENS(ID-1))*
     *         (DM(ID)-DM(ID-1))
             end if
          END DO
          vtot=x
          x=0.
          if(iter.eq.niter) write(6,600)
          do id=1,nd
             AN=DENS(ID)/WMM(ID)+ELEC(ID)
             PGS(ID)=BOLK*TEMP(ID)*AN
             if(id.gt.1) then
             X=X+HALF*(TVISC(ID)/DENS(ID)+TVISC(ID-1)/DENS(ID-1))*
     *         (DM(ID)-DM(ID-1))
             end if
             alp=tvisc(id)/omeg32/pgs(id)*12.5664
             if(iter.eq.niter)
     *       write(6,601) id,dm(id),tvisc(id),thetav(id),x/vtot,
     *                    viscd(id),alp
             write(96,601) id,dm(id),tvisc(id),thetav(id),x/vtot,
     *                    viscd(id),alp
          if(id.eq.nd) 
     *    write(96,601) id,edisc,viscd(id),dens(id),pgs(id),omeg32
          end do
  601     format(i5,1p6e12.4)
       ELSE IF(IVISC.EQ.2) THEN
          X=0.
          THETAV(1)=0.
          DO ID=1,ND
             AN=DENS(ID)/WMM(ID)+ELEC(ID)
             PGS(ID)=BOLK*TEMP(ID)*AN
             TVISC(ID)=OMEG32*ALPHAV*PGS(ID)/12.5664
             DTVIST(ID)=TVISC(ID)/TEMP(ID)
             DTVISN(ID)=TVISC(ID)/AN
             DTVISR(ID)=0.
             if(id.gt.1) then
c            X=X+HALF*(TVISC(ID)+TVISC(ID-1))*(DM(ID)-DM(ID-1))
             X=X+HALF*(TVISC(ID)/DENS(ID)+TVISC(ID-1)/DENS(ID-1))*
     *         (DM(ID)-DM(ID-1))
             end if
          END DO
          VTOT=X
          X=0.
          write(6,602)
  600 format(/'    ID    DM          TVISC   THETAV(orig)   THETAV',
     *        '      VISCD       ALPHA'/)
  602 format(/'    ID    DM          TVISC     THETAV        PGAS',
     *        '       VISCD       ALPHA'/)
          DO ID=1,ND
             if(id.gt.1) then
c            X=X+HALF*(TVISC(ID)+TVISC(ID-1))*(DM(ID)-DM(ID-1))
             X=X+HALF*(TVISC(ID)/DENS(ID)+TVISC(ID-1)/DENS(ID-1))*
     *         (DM(ID)-DM(ID-1))
             end if
             THETAV(ID)=X/VTOT
             viscd(id)=tvisc(id)/dens(id)/edisc
             write(6,601) id,dm(id),tvisc(id),thetav(id),pgs(id),
     *                    viscd(id),alphav
             write(96,601) id,dm(id),tvisc(id),thetav(id),pgs(id),
     *                    viscd(id),alphav
             if(id.eq.nd) 
     *       write(96,601) id,edisc,viscd(id),dens(id),pgs(id),omeg32
          END DO
      END IF
      RETURN
      END