sigmoid diffusivity, command line options, removing comments

code/test.f90

@@ -6,5 +6,6 @@
     PROGRAM test
      USE ysolve
      IMPLICIT NONE
+     CALL parse
      CALL solve
     END PROGRAM test

code/ysolve.f90

@@ -25,8 +25,38 @@
                                 2114624349019./3568978502595.,   &
                                 5198255086312./14908931495163. /)
 
-     PUBLIC :: solve
+     LOGICAL :: read_itape, read_stdin
+     CHARACTER(100) :: itape_name
+     PUBLIC :: solve, parse
     CONTAINS
+!------------------------------------------------------------------------
+     SUBROUTINE parse
+        CHARACTER(100) :: num1char
+
+        !First, make sure the right number of inputs have been provided
+        IF(COMMAND_ARGUMENT_COUNT().EQ.0)THEN
+          ! read from itape
+          read_itape = .true.
+        ELSE IF(COMMAND_ARGUMENT_COUNT().EQ.1)THEN
+          CALL GET_COMMAND_ARGUMENT(1,num1char)   !first, read in the two values
+        ELSE
+          WRITE(*,*)'ERROR, TOO MANY COMMAND-LINE ARGUMENTS(MORE THAN ONE). STOPPING'
+          STOP
+        ENDIF
+
+        read_stdin = (num1char=="-")
+
+        IF(read_itape)THEN
+          itape_name = "itape"
+          WRITE(*,*) "Read from itape"
+        ELSE IF(read_stdin)THEN
+          WRITE(*,*) "Read from STDIN"
+        ELSE 
+          itape_name = num1char
+          WRITE(*,*) "Read from " // num1char
+        END IF
+
+     END SUBROUTINE parse
 !------------------------------------------------------------------------
      SUBROUTINE solve
       INTEGER :: i,j,k,savenum
@@ -104,7 +134,7 @@
        IF (ux(1,i).gt.1.) ux(1,i)=1.                                  
       ENDDO
           
-      dm = (/ (diff * diffusivity_T076(1.-y1(i,1,1)), i=1,nx) /)
+      dm = (/ (diff * diffusivity_sigmoid(1.-y1(i,1,1)), i=1,nx) /)
 
       nd=1
       CALL dfnonp(nx,hx,ux(1,:),dux(1,:),nd,1)  
@@ -255,10 +285,15 @@
 
      SUBROUTINE READ_INTRO
        CHARACTER(LEN=8) :: cdum
-       INTEGER :: itape=5,otape=301,ierr
+       INTEGER :: itape=300,otape=301,ierr
+       INTEGER :: i
 
+       IF (read_stdin) THEN
+         itape=5
+       ELSE
+         OPEN(itape,FILE=itape_name)
+       END IF
 
-       ! OPEN(itape,FILE='itape')
        OPEN(otape,FILE='otape')
         READ(itape,*) cdum,nx
         WRITE(otape,*) cdum,nx
@@ -303,7 +338,13 @@
         READ(itape,*) cdum,cstar
         WRITE(otape,*) cdum,cstar
 
-       ! CLOSE(itape)
+       DO i = 0,100
+        WRITE(*,*) diffusivity_sigmoid(i / 100.)
+       END DO
+
+       IF (.not.read_stdin) THEN
+         CLOSE(itape)
+       END IF
 
 !       hx=l_0*pi/REAL(nx)
        hx=l_0*pi/REAL(nx-1)
@@ -389,7 +430,7 @@
          ux(3,i)=r1(i,j,k)                              ! Y
         ENDDO
 
-      dm = (/ (diff * diffusivity_T076(1.-r1(i,1,1)), i=1,nx) /)
+      dm = (/ (diff * diffusivity_sigmoid(1.-r1(i,1,1)), i=1,nx) /)
 
         CALL dfnonp(nx,hx,ux(1:3,:),dux(1:3,:),3,1)
         CALL dfnonp(nx,hx,dm,duy(1,:),1,1)
@@ -418,50 +459,7 @@
          ENDIF
         ENDDO
        ENDDO
-!! y-direction
-!       DO i=1,nx
-!        DO j=1,ny
-!         uy(1,j)=v(i,j,k)*r1(i,j,k)                     ! v*Y
-!         uy(2,j)=v(i,j,k)                               ! v
-!         uy(3,j)=r1(i,j,k)                              ! Y
-!        ENDDO
-!        CALL dfyz3(ny,hy,uy(1:3,:),duy)
-!        CALL d2fyz1(ny,hy,uy(3:3,:),d2uy)
-!        DO j=1,ny
-!        f(i,j,k)=f(i,j,k) + &
-!        -0.5*( duy(1,j) + uy(2,j)*duy(3,j) + uy(3,j)*duy(2,j) ) &
-!                   + diff*d2uy(j)
-!        ENDDO
-!       ENDDO
       ENDDO
-!      f(1,1,1)=1.0
-!      f(nx,1,1)=-0.5*( dux(1,nx) + ux(2,nx)*dux(3,nx) + ux(3,nx)*dux(2,nx) )
-!      f(nx,1,1)=0.0
-!! z-direction
-!      DO j=1,ny
-!       DO i=1,nx
-!        DO k=1,nz
-!         uy(1,k)=w(i,j,k)*r1(i,j,k)                     ! w*Y
-!         uy(2,k)=w(i,j,k)                               ! w
-!         uy(3,k)=r1(i,j,k)                              ! Y
-!        ENDDO
-!        CALL dfyz3(ny,hy,uy(1:3,:),duy)
-!        CALL d2fyz1(ny,hy,uy(3:3,:),d2uy)
-!        DO k=1,nz
-!        f(i,j,k)=f(i,j,k) + &
-!        -0.5*( duy(1,k) + uy(2,k)*duy(3,k) + uy(3,k)*duy(2,k) ) &
-!                        + diff*d2uy(k)
-!        ENDDO
-!       ENDDO
-!      ENDDO
-
-!      DO j=1,ny
-!       DO k=1,nz
-!         f(1,j,k)=0.0
-!         f(nx,j,k)=-0.5*(dux(1,nx)+ux(2,nx)*dux(3,nx)+ux(3,nx)*&
-!           dux(2,nx))
-!       ENDDO
-!      ENDDO
 
 !    Boundary conditionS
       f(1,1,1)=0.
@@ -478,7 +476,7 @@
 !------------------------------------------------------------------------
      REAL FUNCTION diffusivity_sigmoid(c)
          REAL, INTENT(IN) :: c
-         REAL :: a = 200.
+         REAL :: a = 100.
          diffusivity_sigmoid = 1. + (rvis - 1.) * sigmoid(a*(c-cstar))
      END FUNCTION diffusivity_sigmoid
 !------------------------------------------------------------------------