c   $Id: readprintpar.f,v 1.1 1998/10/20 13:23:03 zaleski Exp $
c   (c) S. Zaleski, CNRS, UPMC
c
c   $Author: zaleski $
c
c   New initialisation subroutine
c   All input must be in SI units 
c
c
c
      SUBROUTINE readglob(tmax,tprint)
      include 'undefined.h'
      INTEGER tmax,tprint,i
      DOUBLE PRECISION vel1,vel2, tau1
      INCLUDE 'size.h'
      INCLUDE 'global.f'
      DOUBLE PRECISION lambdasi,tausi,sigmasi,rholiqsi,rhogassi,
     *  viscgassi, viscliqsi, xvelsi, yvelsi,
     *  rholiq, gSI,duSI,circulSI,circul
      INTEGER GP_ifetch, ifetch, dfetch

c circul passed directly to inivort

      common /COMVORT/ circul

c parameters passed to in_and_out

      common /INOUT/ vel1,vel2, tau1

c  First read SI values
c
c 1) compulsory parameters

      if(ifetch("n1",nx).ne.1) then
         if(ifetch("nx",nx).ne.1) stop 'missing n1'
      endif
      if(ifetch("n2",ny).ne.1) then
         if(ifetch("ny",ny).ne.1) stop 'missing n2'
      endif
      nz=1

      if(dfetch("tau-si",tauSI).ne.1)   stop 'missing tau-si'
      if(dfetch("lambda-si",lambdaSI).ne.1)  then
        if(dfetch("L-si",lambdaSI).ne.1) 
     *      stop 'missing lambda-si or L-si'
      else
         if(dfetch("L-si",lambdaSI).eq.1) 
     *      stop 'Both L-si and lambda-si defined'
      endif

      if(dfetch("rhogas-si",rhogasSI).ne.1) stop 'missing rhogas-si'
      if(dfetch("rholiq-si",rholiqSI).ne.1) stop 'missing rholiq-si'
c
      if(ifetch("tmax",tmax).ne.1) stop 'missing tmax'
      if(ifetch("tprint",tprint).ne.1)   stop 'missing tprint'
c
c Optional parameters -- default values. 
c
      xvelSI=0.d0
      yvelSI=0.d0
      sigmaSI=0.d0
      i=dfetch("xvel-si",xvelSI)
      i=dfetch("yvel-si",yvelSI)
      i=dfetch("sigma-si",sigmaSI)
      if(dfetch("viscgas-si",viscgasSI).ne.1)  
     *  stop 'missing viscgas-si'
      if(dfetch("viscliq-si",viscliqSI).ne.1)  
     *  stop 'missing viscliq-si'
c
c    Now we have a scale of mass ( rholiqSI*L^3 ) and  of length (l-SI)
c
c  We determine the time scale from 
c   A) the velocity scale if present 
c   B) Otherwise, from the gravity scale gSI if present (I.E. we are not
c        in outer space)
c   C) If none of the above, timescale=1second.
c
c
      gSI=0.
      duSI=0.
      circulSI=0.

      i=dfetch("gravity-si",gSI)
      i=dfetch("deltau-si",duSI)
      i=dfetch("circul-si",circulSI)
      

c ifdrop

      if((yvelSI.gt.0.).or.(xvelSI.gt.0.)) then
         duSI = dmax1(dabs(yvelSI),dabs(xvelSI))
         if(GP_ifetch("acceleration",0,"WARNING").eq.1) then
            yvelSI=0.d0
         endif
      endif

c  Determine the timescale from yvelSI, if not zero.  Then duSI, if duSI 
c  is zero try other methods. - Changed 3/7/2007 - RDS
       if(yvelSI.ne.0.d0) then
          write(6,*) "determining the time scale from yvelSI"
          timescaleSI = lambdaSI/abs(yvelSI)
          du = duSI/abs(yvelSI)
          g = gSI / (lambdaSI/(timescaleSI**2))
c	 Note that we need to non-dimensionalize yvel no matter what, so
c	 that we get the sign right! - RDS 3/7/2007
      else if(duSI.ne.0.d0) then
         write(6,*) "determining the time scale from duSI"
         timescaleSI = lambdaSI/duSI
         du = 1.d0
         g=gSI / (lambdaSI/(timescaleSI**2))
      else if(gSI.ne.0.d0) then
         write(6,*) "determining the time scale from gSI"
         timescaleSI= dsqrt(lambdaSI/dabs(gSI))
         g=1.d0*gSI/dabs(gSI)
         du=0.d0
      else 
         write(6,*) "setting the time scale at 1s"
         timescaleSI=1
         g=0.d0
         du=0.d0
      endif

c    We now convert everything into dimensionless code units 

      tau = tauSI/timescaleSI
      rholiq=1
      rhogas=rhogasSI/rholiqSI
      mugas=viscgasSI*timescaleSI/(rholiqSI*lambdaSI**2)
      muliq=viscliqSI*timescaleSI/(rholiqSI*lambdaSI**2)
      sigma=sigmaSI*timescaleSI**2/(rholiqSI*lambdaSI**3)
      xvel=xvelSI*timescaleSI/lambdaSI
      yvel=yvelSI*timescaleSI/lambdaSI
      circul = circulSI/(lambdaSI**2/timescaleSI)

c Compute values for common block /INOUT/

      vel1=xvel
      vel2=xvel+du
      tau1=tau
      
      write(6,*) "deltau-si",duSI
      write(6,*) "lambda-si",lambdaSI
      write(6,*) "tau-si=",tausi
      write(6,*) "sigma-si=",sigmasi
      write(6,*) "rhogas-si=",rhogassi
      write(6,*) "rholiq-si=",rholiqsi
      write(6,*) "viscgas-si=",viscgassi
      write(6,*) "viscliq-si=",viscliqsi
      write(6,*) "timescale-si=",timescaleSI
c      write(6,*) "inner nozzle velocity",vel1
c      write(6,*) "outer nozzle velocity",vel2

      return
      end
c
      SUBROUTINE printhist(tprint)
      include 'undefined.h'
      DOUBLE PRECISION vmax,cfl,h,tauh
      INTEGER tprint
      INCLUDE 'size.h'
      INCLUDE 'global.f'

      h = 1.d0/(nx-2)
      tauh=tau/h
      vmax = dmax1(dabs(xvel),dabs(yvel))
      cfl = tau*vmax/h
      if(cfl.gt.MAXCFL) write(6,*) "cfl=",cfl
      write(6,*) "d3=", tprint*tau
      write(6,*) "tau=",tau
      write(6,*) "sigma=",sigma
      write(6,*) "rhogas=",rhogas
      write(6,*) "viscgas=",mugas
      write(6,*) "viscliq=",muliq
      write(6,*) "xvel=",xvel
      write(6,*) "yvel=",yvel
      write(6,*) "g=",g
      return
      end