geometry.f90

Go to the documentation of this file.

!*********************************************************************
!
!*********************************************************************
MODULE geometry

  IMPLICIT NONE

  REAL*8, ALLOCATABLE :: z_comp(:)

  REAL*8, ALLOCATABLE :: z_stag(:)

  REAL*8, ALLOCATABLE :: radius_stag(:)

  REAL*8 :: z0
  REAL*8 :: zn
  REAL*8 :: zeta_exit
  REAL*8 :: dz

  REAL*8 :: pi

  REAL*8 :: radius

  CHARACTER*30 :: radius_model

  REAL*8 :: radius_fixed
  REAL*8 :: radius_min
  REAL*8 :: radius_max
  REAL*8 :: radius_z
  REAL*8 :: radius_z_sig

  INTEGER :: comp_cells
  INTEGER :: comp_interfaces

CONTAINS

  !*********************************************************************
  !
  !*********************************************************************

  SUBROUTINE init_grid

    IMPLICIT none

    INTEGER j

    pi = 4.d0 * datan(1.d0)

    comp_interfaces = comp_cells+1

    ALLOCATE( z_comp(comp_cells) )
    ALLOCATE( z_stag(comp_interfaces) )
    ALLOCATE( radius_stag(comp_interfaces) )
    
    dz = ( zn - z0 ) / comp_cells

    z_stag(1) = z0
    z_comp(1) = z0 + 0.5 * dz
    
    DO j=1,comp_cells

       z_stag(j+1) = min( z_stag(j) + dz , zn )

       z_comp(j) = 0.5 * ( z_stag(j) + z_stag(j+1) )

    END DO

    SELECT CASE ( radius_model )
 
    CASE DEFAULT

       radius_stag(1:comp_interfaces) = radius_fixed

      
    CASE ('fixed' )

       radius_stag(1:comp_interfaces) = radius_fixed

    CASE ( 'linear' )

       DO j=1,comp_interfaces
          
          radius_stag(j) = radius_min + (radius_max - radius_min) * z_stag(j) / (zn - z0)
          
       END DO
       
    CASE ( 'trans1' )

       DO j=1,comp_interfaces
          
          IF( z_stag(j) > (radius_z + radius_z_sig) ) THEN
             
             radius_stag(j) = radius_max
             
          ELSEIF( z_stag(j) < (radius_z - radius_z_sig) ) THEN
             
             radius_stag(j) = radius_min
             
          ELSE

             radius_stag(j) = radius_min + (radius_max - radius_min) * &
                  ( z_stag(j) - radius_z + radius_z_sig) / (2.d0 * radius_z_sig)
             
          END IF
          
       END DO
 
    CASE ( 'trans2' )
      
       DO j=1,comp_interfaces
          
          IF( z_stag(j) > (radius_z + radius_z_sig) ) THEN
             
             radius_stag(j) = radius_min
             
          ELSEIF( z_stag(j) > (radius_z) ) THEN
             
             radius_stag(j) = radius_max - (radius_max - radius_min) * &
                  ( z_stag(j) - radius_z ) / ( radius_z_sig )
             
          ELSEIF( z_stag(j) > (radius_z - radius_z_sig) ) THEN
             
             radius_stag(j) = radius_min + (radius_max - radius_min) * &
                  ( z_stag(j) - radius_z + radius_z_sig) / ( radius_z_sig )
             
          ELSE
             
             radius_stag(j) = radius_min
             
          END IF
          
       END DO

    CASE ( 'external' )
       
       OPEN( unit=10, file='DataRadius.txt' )
       
       DO j=1,comp_interfaces
          
          READ(10,*) radius_stag(j)
          
       END DO
       
       CLOSE(10)
       
    END SELECT

  END SUBROUTINE init_grid

  SUBROUTINE update_radius(zeta)

    IMPLICIT NONE

    REAL*8, INTENT(IN) :: zeta
    REAL*8 :: coeff_interp
    INTEGER :: j , z_idx

    z_idx = 1

    DO j=1,comp_interfaces-1

       IF ( z_stag(j) < zeta ) z_idx = j

    END DO    

    ! zeta is between z_stag(j) and z_stag(j+1)
    
    coeff_interp = ( zeta - z_stag(z_idx) ) / ( z_stag(z_idx+1) - z_stag(z_idx) )  

    radius = coeff_interp * radius_stag(z_idx+1) + ( 1.d0 - coeff_interp ) *        &
         radius_stag(z_idx)

  END SUBROUTINE update_radius

END MODULE geometry