text:
c----------------------------------------------------------------
      function ichash(atxt,rtxt,ntxt,ctxt)
c
c produce hash number from atom and residue name,rsidue number and chain
c name for charge assignment
c
	include 'qdiffpar4.h'
c	include 'qdiffpar1.h'

      character*6 atxt
      character*3 rtxt
      character*4 ntxt
      character*1 ctxt
c
c	character*1 chn
c	character*3 res
c	character*4 rnum
c	character*6 atm
c
      character*38 string
	integer n,irhash,j
      data string /'* 0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ'/
      n = 1
      do 9000 i = 1,3
        j = index(string,rtxt(i:i))
        n = 5*n + j
c      end do
9000	continue
      do 9001 i = 1,6
        j = index(string,atxt(i:i))
        n = 5*n + j
c      end do
9001	continue
      do 9002 i = 1,4
        j = index(string,ntxt(i:i))
        n = 5*n + j
c      end do
9002	continue
      do 9003 i = 1,1
        j = index(string,ctxt(i:i))
        n = 5*n + j
c      end do
9003	continue
	n = iabs(n)
      ichash = mod(n,nclist) + 1
      return
      end

max_stars_repo_name: caixiuhong/Stable-MCCE

max_issues_repo_name: caixiuhong/Stable-MCCE

max_forks_repo_name: caixiuhong/Stable-MCCE


text:
	SUBROUTINE IP_GSAV  ( pvar, iframe, nframe, iret )
C************************************************************************
C* IP_GSAV								*
C*									*
C* This subroutine exports the current X window into an output file	*
C* specified by the file extension in pvar.				*
C*									*
C* IP_GSAV  ( PVAR, IFRAME, NFRAME, IRET )				*
C*									*
C* Input parameters:							*
C*	PVAR		CHAR*		Input parameter name		*
C*	IFRAME		INT		Frame index number		*
C*	NFRAME		INT		Number of frames		*
C*									*
C* Output parameters:							*
C*	IRET		INTEGER		Return code			*
C*				   	  0 = normal return 		*
C*					-13 = Invalid window name	*
C**									*
C* Log:									*
C* C. Bailey/HPC	02/05						*
C************************************************************************
	INCLUDE		'ipcmn.cmn'
C*
	CHARACTER*(*)	pvar
	INTEGER		iframe, nframe, lenw
C-----------------------------------------------------------------------
	iret = 0
C
C* 	Add Null character to the end of the filnam string
C
	CALL ST_NULL(pvar, pvar, lenw, iret)
C
C*	Export the current window to a file.
C
	CALL GGSAVE ( pvar, iframe, nframe, ierr )
C
C*	Check for an invalid window name and write an error message.
C
	IF  ( ierr .ne. 0 ) THEN
	    CALL ER_WMSG ( 'GEMPLT', ierr, ' ', ier )
	END IF
C*
	CALL GEPLOT ( ier )
C*
	RETURN
	END

ext: f

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
Name: Lois Wolk

Office: California State Senate, 5th District

Website: http://sd03.senate.ca.gov/

Email: MailTo(senator.wolk@senate.ca.gov)

Phone Number: (916) 6514005

Other Information: Has been a community activist for many years.  Has served as Mayor of Davis, Yolo County Board of Supervisors county supervisor, assemblymember, and state senator.

As a wife and mother, Lois Wolk raised two children, one of whom is current Davis City Council member, Dan Wolk, and served as PTA president, AYSO and Little League coach.  She was a charter member of the board of the Explorit Science Center.  She served on the citys Department of Parks and Recreation following her successful involvement with a campaign to expand the Davis Branch Library.

She was elected to the California State Legislature state legislature in 2002 in the 8th Assembly District and to the 5th Senate District in 2008. In the Assembly she was succeeded by Mariko Yamada.

Lois Wolk served on the City Council from 6/20/1990 to 4/10/2002.  Served as Mayor Pro Tem from 1990 to 1991, Mayor from 1992 to 1993, Mayor Pro Tem in 1994 to 1995, and again Mayor from 1996 to 1998.

She spoke on the current assembly legislation and took questions from members of Davis College Democrats on October 17th 2007.

 Endorsements

Assemblymember Wolk has received ASUCD endorsements in the past.  A constitutional amendment had been passed that prevented ASUCD officials from endorsing political candidates, as it was a violation of ASUCDs taxexempt status.  Her endorsements from ASUCD officials used to be available at http://www.loiswolk.net/endorse.html  but now you have to go to the http://web.archive.org/web/20050312075356/http://www.loiswolk.net/endorse.html internet archive.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
!> Pimpact
!! \author huppd
!! \date 2018

!>  Modul: cmod_VectorField
!!
!! Impact functions \c Pimpact::VectorField e.g. scales, norms ...
module cmod_VectorField

  use iso_c_binding

  implicit none

contains

  !> \brief calculate distance to immersed boundary
  !!
  !! \param[in] x
  !! \param[in] y
  !! \param[in] z
  !! \param[in] x0
  !! \param[in] y0
  !! \param[in] R
  !! \param[in] dr
  FUNCTION distance2ib( x, y, z, x0, y0, R, dr ) RESULT(dis)

    IMPLICIT NONE

    REAL(c_double), INTENT(IN)  ::  x,y,z

    REAL(c_double), intent(in)  ::  x0
    REAL(c_double), intent(in)  ::  y0

    REAL(c_double), intent(in)  ::  R

    REAL(c_double), intent(in)  ::  dr

    REAL(c_double)              ::  z0

    REAL(c_double)              ::  dis

    ! geometric properties of disc
    !        x0 = L1/2. + L1*amp*SIN(2.*pi*freq*subtime)
    !        y0 = L2/4.
    z0 = 0.

    !        R = L1/10.
    !  write(*,*) 'dr=',dr

    dis = SQRT( (x0-x)**2 + (y0-y)**2 + (z0-z)**2 )

    if( dis <= R) then
      dis = 1.
    else if( dis <= R+dr) then
      dis = ABS(1./( 1. + EXP(1./(-ABS(dis-R)/dr) + 1./(1.-Abs(dis-R)/dr) ) ));
    else
      !IF( Abs(dis-R) < dr) THEN
      !dis = ABS(1./( 1. + EXP(1./(-ABS(dis-R)/dr) + 1./(1.-Abs(dis-R)/dr) ) ));
      !!    dis = (1-(dis-R)/dr)
      !ELSE
      dis = 0.
    end if
    !    dis = 0.

  END FUNCTION distance2ib

  !> \brief get initial slope/curvature for boundary layer equations as a
  !! function of sweep angle and suction velocity
  !!
  !! \param[in] kappa non-dimensional boundary layer suction velocity
  !! \param[in] sweep_angle_degrees sweep angle (in degrees for ease of interpolation)
  !! \param[out] fxxGuess final guess values (return values)
  !! \param[out] gxGuess final guess values (return values)
  SUBROUTINE init_value_bl_equation( kappa, sweep_angle_degrees, fxxGuess, gxGuess )

    implicit none

    !----------- mjohn 101111 ------------------------------------------------------------------
    ! - INPUT VARIABLES
    real(c_double), intent(in)  :: kappa
    real(c_double), intent(in)  :: sweep_angle_degrees

    ! - OUTPUT VARIABLES
    real(c_double), intent(out) :: fxxGuess, gxGuess

    real(c_double)   :: fxx(1:10,7), gx(1:10,7)    ! initial value matrices for f and g (structure see where declared)
    real(c_double)   :: fxx1, fxx2, fxx3, fxx4     ! interpolation stencil (2D, since kappa and sweep_angle vary) for f
    real(c_double)   :: gx1, gx2, gx3, gx4         ! interpolation stencil (2D, since kappa and sweep_angle vary) for g
    real(c_double)   :: dK, dPhi                   ! interpolation increment remainders (interpolation weights)
    real(c_double)   :: fxxLo, fxxHi, gxLo, gxHi   ! more dummy variables...
    !------------------------------------------------------------------------------------------

    ! v''(0) for kappa = (0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0) and
    ! phi =       (     0        10        20        30        40        50        60        70        80      90   )
    fxx(1:10,1) = (/ 1.232578, 1.204596, 1.122776, 0.993369, 0.826410, 0.635208, 0.435783, 0.246543, 0.089191, 0.0 /)
    fxx(1:10,2) = (/ 1.541729, 1.509206, 1.413891, 1.262421, 1.065455, 0.854631, 0.594496, 0.357035, 0.147418, 0.0 /)
    fxx(1:10,3) = (/ 1.889298, 1.851881, 1.742007, 1.566680, 1.337148, 1.068325, 0.778115, 0.486851, 0.217475, 0.0 /)
    fxx(1:10,4) = (/ 2.279667, 2.236852, 2.110903, 1.909177, 1.643485, 1.329463, 0.985815, 0.633578, 0.295836, 0.0 /) ! DUMMY VALUES (mere interpolation of line above and below)
    fxx(1:10,5) = (/ 2.670036, 2.621822, 2.479798, 2.251675, 1.949821, 1.590600, 1.193515, 0.780304, 0.374196, 0.0 /)
    fxx(1:10,6) = (/ 3.098327, 3.044144, 2.884339, 2.626989, 2.285053, 1.884300, 1.419229, 0.938387, 0.457070, 0.0 /) ! DUMMY VALUES (mere interpolation of line above and below)
    fxx(1:10,7) = (/ 3.526617, 3.466466, 3.288879, 3.002302, 2.620284, 2.178000, 1.644944, 1.096470, 0.539944, 0.0 /)
    ! w'(0):
    gx(1:10,1) = (/ 0.570454, 0.566105, 0.552986, 0.530868, 0.499286, 0.457358, 0.403374, 0.333619, 0.237717, 0.0 /)
    gx(1:10,2) = (/ 0.921733, 0.917607, 0.905180, 0.884295, 0.854631, 0.815588, 0.766028, 0.703608, 0.622597, 0.5 /)
    gx(1:10,3) = (/ 1.323674, 1.319995, 1.308946, 1.290479, 1.264492, 1.230778, 1.188932, 1.138139, 1.076642, 1.0 /)
    gx(1:10,4) = (/ 1.767516, 1.764282, 1.754585, 1.738445, 1.715877, 1.686874, 1.651353, 1.609055, 1.559273, 1.5 /) ! DUMMY VALUES (mere interpolation of line above and below)
    gx(1:10,5) = (/ 2.211359, 2.208569, 2.200225, 2.186411, 2.167263, 2.142970, 2.113775, 2.079971, 2.041903, 2.0 /)
    gx(1:10,6) = (/ 2.682235, 2.679764, 2.672382, 2.660189, 2.643351, 2.622097, 2.596728, 2.567613, 2.535191, 2.5 /) ! DUMMY VALUES (mere interpolation of line above and below)
    gx(1:10,7) = (/ 3.153111, 3.150959, 3.144539, 3.133967, 3.119438, 3.101224, 3.079682, 3.055254, 3.028479, 3.0 /)

    ! get interpolation stencil for v'' between starting values in sweep_angle and kappa direction
    fxx1 = fxx(floor(  sweep_angle_degrees/10.)+1,floor(  kappa/0.5)+1)
    fxx2 = fxx(floor(  sweep_angle_degrees/10.)+1,ceiling(kappa/0.5)+1)
    fxx3 = fxx(ceiling(sweep_angle_degrees/10.)+1,floor(  kappa/0.5)+1)
    fxx4 = fxx(ceiling(sweep_angle_degrees/10.)+1,ceiling(kappa/0.5)+1)
    ! get interpolation stencil for w' between starting values in sweep_angle and kappa direction
    gx1 = gx(floor(  sweep_angle_degrees/10.)+1,floor(  kappa/0.5)+1)
    gx2 = gx(floor(  sweep_angle_degrees/10.)+1,ceiling(kappa/0.5)+1)
    gx3 = gx(ceiling(sweep_angle_degrees/10.)+1,floor(  kappa/0.5)+1)
    gx4 = gx(ceiling(sweep_angle_degrees/10.)+1,ceiling(kappa/0.5)+1)

    ! get remainders in kappa and sweep_angle dimension
    dK = mod(kappa,0.5)
    dPhi = mod(sweep_angle_degrees,10.0)

    ! interpolate inside 4 point stencil for v''
    fxxLo = ((0.5-dK)*fxx1 + dK*fxx2)/0.5
    fxxHi = ((0.5-dK)*fxx3 + dK*fxx4)/0.5
    fxxGuess = ((10-dPhi)*fxxLo + dPhi*fxxHi)/10.0
    ! interpolate inside 4 point stencil for w'
    gxLo = ((0.5-dK)*gx1 + dK*gx2)/0.5
    gxHi = ((0.5-dK)*gx3 + dK*gx4)/0.5
    gxGuess = ((10-dPhi)*gxLo + dPhi*gxHi)/10.0

    !   ! --- FOR DEBUGGING PURPOSES ONLY
    !    write(*,*) 'Kappa and Phi are', kappa, ', ', sweep_angle_degrees, '.'
    !    write(*,*) 'Interpolation for f between the four values'
    !    write(*,*) fxx1, ', ', fxx2, ', ', fxx3, ', ', fxx4
    !    write(*,*) 'resulted in initial guess', fxxGuess, '.'
    !    write(*,*)
    !    write(*,*) 'Interpolation for g between the four values'
    !    write(*,*) gx1, ', ', gx2, ', ', gx3, ', ', gx4
    !    write(*,*) 'resulted in initial guess', gxGuess, '.'

  END SUBROUTINE init_value_bl_equation

  !> \brief shooting method for SHBL
  !!
  !! \param[out] v solution vector
  !! \param[in] grid array of grid points
  !! \param[in] n number of grid points
  !! \param[in] sweep_angle sweep angle PHI
  !! \param[in] kappa boundary suction KAPPA
  !! \param[in] s starting value for shooting
  !! \param[in] r end value for shooting
  SUBROUTINE shoot_v(v,grid,n,sweep_angle_degrees,kappa,s,r)

    implicit none

    !-----------------------------------------------------------------------
    integer(c_int), intent(in)    :: n
    real(c_double), intent(out)   :: v(n,3)
    real(c_double), intent(in)    :: grid(n)
    real(c_double), intent(inout) :: s,r
    real(c_double), intent(in)    :: sweep_angle_degrees
    real(c_double), intent(in)    :: kappa
    !-------------------------------------------------------------------------
    real(c_double) :: rhs(3)        !value of right-hand side
    real(c_double) :: g(3)          !intermediate value in Runge-Kutta integration
    real(c_double) :: t(3)          !intermediate value of v at x
    real(c_double) :: rk1(3), rk2(3)!Runge-Kutta coefficients
    real(c_double) :: dx            !integration step
    real(c_double) :: x             !independent variable
    integer(c_int) :: i, k          !counters
    real(c_double) :: cosPhi        !cos(sweep_angle)
    real(c_double) :: upperBound    !upper interval boundary for integration
    real(c_double) :: pi
    !=====================================================================

    pi = 4.*atan(1.)    !!set constants

    ! set Runge-Kutta coefficients
    rk1(1) = 0.
    rk1(2) = -5./9.
    rk1(3) = -153./128.
    rk2(1) = 1./3.
    rk2(2) = 15./16.
    rk2(3) = 8./15.

    !!set initial condition: v(0) = kappa, v'(0) = 0, initial guess depends on kappa, phi, alpha
    i = 1
    upperBound = 50.0

    cosPhi   = cos(sweep_angle_degrees*pi/180.)

    t(1) = -kappa
    t(2) = 0.
    t(3) = s
    x = 0.

    do
      ! -- mjohn 111111
      !! catch values in order to prevent divergence collapse (floating point overflow)
      if( abs(t(1)) >= 1e6 ) then
        r = t(2) + cosPhi;
        exit
      end if
      ! -- mjohn 111111

      !! save error and integration value
      r = t(2) + cosPhi;
      IF (x .EQ. grid(i)) THEN
        v(i,1:3) = t
        i = i + 1
      end if

      !! exit condition (compute until full grid is integrated AND 30 is reached)
      IF (i > n .and. x >= upperBound) EXIT

      !! standard step size
      dx = 1.e-5

      !! reduce step size accordingly if grid point is approached
      IF (i<=n .and. x+dx > grid(i)) THEN
        dx = grid(i)-x
        x = grid(i)
      else
        x = x+dx
      END IF

      !!three step runge-kutta
      g = 0.

      do k = 1, 3
        if( x <= upperBound ) then
          rhs(1) = t(2)
          rhs(2) = t(3)
          rhs(3) = t(3)*t(1) + cosPhi**2 - t(2)**2
        else
          rhs(1) = t(2)
          rhs(2) = 0.
          rhs(3) = 0.
        end if

        g = rk1(k)*g + rhs
        t = t + rk2(k)*dx*g

      end do

    end do

  END SUBROUTINE shoot_v

  SUBROUTINE shoot_w(w,v,grid,n,sweep_angle_degrees,kappa,sv,sw,r,blThick)

    implicit none

    !!!-----------------------------------------------------------------------
    integer(c_int), intent(in) :: n                     !> number of grid points
    real(c_double), intent(out):: w(n,2)                !> OUTPUT: solution vector
    real(c_double), intent(out):: v(n,3)                !> OUTPUT: solution vector
    real(c_double), intent(in) :: grid(n)               !> array of grid points
    real(c_double), intent(inout):: sv,sw,r             !> starting and end value for shooting
    real(c_double), intent(in)   :: sweep_angle_degrees !> sweep angle PHI
    real(c_double), intent(in)   :: kappa               !> boundary suction KAPPA
    real(c_double), intent(out)  :: blThick             !> blThickness (mathematical, non-dimensionalized)
    !!!-----------------------------------------------------------------------
    real(c_double) :: rhs(5)        !value of right-hand side
    real(c_double) :: g(5)          !intermediate value in Runge-Kutta integration
    real(c_double) :: t(5)          !intermediate value of v at x
    real(c_double) :: rk1(3), rk2(3)!Runge-Kutta coefficients
    real(c_double) :: dx            !integration step
    real(c_double) :: x             !independent variable
    integer(c_int) :: i, k          !counters
    real(c_double) :: cosPhi        !cos(sweep_angle)
    real(c_double) :: upperBound    !upper interval boundary for integration
    real(c_double) :: pi
    !!!=====================================================================

    pi = 4.*atan(1.)    !!set constants

    !!set Runge-Kutta coefficients
    rk1(1) = 0.
    rk1(2) = -5./9.
    rk1(3) = -153./128.
    rk2(1) = 1./3.
    rk2(2) = 15./16.
    rk2(3) = 8./15.

    !!set initial condition: v(0) = kappa, v'(0) = 0, initial guess depends on kappa, phi, alpha
    cosPhi = cos(sweep_angle_degrees*pi/180.)

    i = 1
    t(1) = -kappa
    t(2) = 0.
    t(3) = sv
    t(4) = 0.
    t(5) = sw
    x = 0.

    !! w needs special treatment for phi = 90° and kappa = 0, since solution becomes trivial
    !! activate VERY large integration domain for large phi and small kappa
    if (kappa < 0.01 .AND. sweep_angle_degrees > 85.0) then
      upperBound = min(2/(kappa+epsilon(upperbound)),200.0);
    else
      upperBound = 50.0;
    end if

    ! initialize boundary layer thickness to 0
    blThick = 0;

    !! numerical integration loop (let coordinate x run from 0 to infinity)
    do

      !!save values
      r = t(4) - 1.;
      IF (x .EQ. grid(i)) THEN
        v(i,1:3) = t(1:3)
        w(i,1:2) = t(4:5)
        i = i + 1
      END IF

      !! exit condition
      IF (i > n .and. x >= upperBound) EXIT
      dx = 1.e-5

      IF (i<=n .and. x+dx > grid(i)) THEN
        dx = grid(i)-x
        x = grid(i)
      ELSE
        x = x+dx
      END IF

      !!three step runge-kutta
      g = 0.
      do k = 1, 3
        if (x < upperBound) then
          rhs(1) = t(2)
          rhs(2) = t(3)
          rhs(3) = t(3)*t(1) + cosPhi**2 - t(2)**2
          rhs(4) = t(5)
          rhs(5) = t(5)*t(1)
        else
          rhs(1) = t(2)
          rhs(2) = 0.
          rhs(3) = 0.
          rhs(4) = 0.
          rhs(5) = 0.
        end if
        g = rk1(k)*g + rhs
        t = t + rk2(k)*dx*g
      end do

      !! carry out first order integral: int_0^\inf (g * (1-g)) dx
      blThick = blThick + (t(4)*(1-t(4))) *dx;

      if(blThick >= 0.25 * upperBound) then
        write(*,*) 'WARNING. B.L. THICKNESS GREATER THAN 25 % OF INTEGRATION INTERVAL. BASE FLOW PROFILE MIGHT BE WRONG.'
      end if

    end do

  END SUBROUTINE shoot_w

  !> \brief calculate a laminar swept attachment-line flow
  !!
  !! based on calcbasicflow by obristd, 1999/10/11
  !! extended by mjohn, 2011/11/08
  !! Purpose: calculate a generalized 3D boundary layer flow
  !! with SWEEP ANGLE and boundary suction (wall-normal) KAPPA
  !!
  !! 11/11/28: modified: baseflow now contains "real" profile, not only
  !! integration routine, i.e. multiplication with 1/Re, sin(phi), x3
  !! is performed before ub, vb, wb are returned
  !!
  !! 12/10/05: modified: baseflow may now be non-dimensionalized
  !! according to SHBL formalism (1 d.o.f.), wing formalism (2 d.o.f.)
  !! or novel scaling (1 d.o.f.)
  !!
  !! \param[in] rank
  !! \param[in] kappa boundary suction KAPPA
  !! \param[in] sweep_angle_degrees sweep angle PHI
  !! \param[in] n number of grid points
  !! \param[in] grid array of grid points
  !! \param[inout] ub chordwise baseflow profile u
  !! \param[inout] vb wall-normal baseflow profile v
  !! \param[inout] wb spanwise baseflow profile w
  !! \param[out] blthick
  subroutine calcbasicflow( &
      rank,                 &
      kappa,                &
      sweep_angle_degrees,  &
      n,                    &
      grid,                 &
      ub,                   &
      vb,                   &
      wb,                   &
      blthick )

    implicit none

    integer(c_int), intent(in)    :: rank
    real(c_double), intent(in)    :: kappa
    real(c_double), intent(in)    :: sweep_angle_degrees
    integer(c_int), intent(in)    :: n
    real(c_double), intent(in)    :: grid(n)
    real(c_double), intent(inout) :: ub(n)
    real(c_double), intent(inout) :: vb(n)
    real(c_double), intent(inout) :: wb(n)

    real(c_double),allocatable :: v(:,:)       ! dependent variable v and its derivatives
    real(c_double),allocatable :: w(:,:)       ! dependent variable w and its derivatives

    real(c_double), intent(out) :: blThick

    real(c_double) :: initv1, initv2 !starting values for v equation
    real(c_double) :: initw1, initw2 !starting values for w equation
    real(c_double) :: sv, s0, s1, s1old  !starting values
    real(c_double) :: r0, r1         !values of object function
    real(c_double) :: vxxGuess, wxGuess                ! final guess values (return values)  !-- mjohn 101111

    integer(c_int) :: i

    real(c_double) :: pi

    allocate(v(n,3),w(n,2))   !! allocate local memory
    pi = 4.*atan(1.)    !!set constants

    !!! ---------------------------------------------------------- INITIALIZE SHOOTING ROUTINE (GET STARTING VALUES) ----------------------------------------------------
    !----------- mjohn 101111
    ! -- get initial slope/curvature for boundary layer equations as a function of sweep angle and suction velocity
    call init_value_bl_equation( kappa, sweep_angle_degrees, vxxguess, wxguess )

    !!!*** solve equation for v
    ! move to VectorField::initField
    if( kappa>3.0 .or. kappa<0.0 .or. sweep_angle_degrees>90.0  .or. sweep_angle_degrees<0.0 ) then
      WRITE(*,*) 'WARNING: kappa or sweep angle outside allowed interval!! Shooting integration might fail. rank: ', rank
    end if

    initv1 = -vxxGuess*(1+1e-3) !value of ddv for first shot
    initv2 = -vxxGuess*(1-1e-3) !value of ddv for second shot
    initw1 =   wxGuess*(1+1e-3) !value of dw for first shot
    initw2 =   wxGuess*(1-1e-3) !value of dw for second shot
    !----------- mjohn 101111

    !!! ---------------------------------------------------------- SOLVE SHOOTING INTEGRATION PROBLEM FOR VELOCITY V ----------------------------------------------------
    if( rank == 0 ) then
      WRITE(*,*) ' Solving equation v'''''' = cos^2(sweep_angle) - v''^2  + v v'''' for v'
    end if

    !!set initial values
    s0 = initv1
    s1 = initv2

    !! two initial shots
    call shoot_v(v,grid,n,sweep_angle_degrees,kappa,s0,r0)
    call shoot_v(v,grid,n,sweep_angle_degrees,kappa,s1,r1)

    i = 1
    !! controlled shooting loop
    do

      !!calculate new initial value (only if residuum is large enough!)
      if( abs(r1) > 1E-13 .and. abs(r0-r1) > 1E-13) then
        s1old = s1
        s1 = (r0*s1 - r1*s0)/(r0-r1)
        s0 = s1old
        r0 = r1
      end if

      call shoot_v(v,grid,n,sweep_angle_degrees,kappa,s1,r1)

      if(abs(r1) <= 1e-10) then
        !!goal reached
        exit
      elseif ((s1 == s0) .or. (r1 == r0)) then
        !!no more progress
        exit
      end if
    end do

    !! save starting value
    sv = s1

    !! print output
    if( rank==0 ) then
      !      write(*,*) '   Ended shoot_v successfully. ddy(v)|0 = ', v(1,3)
      write(*,*) ' residuum =', r1
      !      if ( ((initv1+initv2)/2.0) .NE. 0.0) then
      !         write(*,*) 'Relative deviation from initial guess:', (v(1,3)-((initv1+initv2)/2))/((initv1+initv2)/2)
      !      else
      !         write(*,*) 'Total deviation from initial guess:   ', v(1,3)
      !      end if
      !      write(*,*)
    end if

    !!! ---------------------------------------------------------- SOLVE SHOOTING INTEGRATION PROBLEM FOR VELOCITY W ----------------------------------------------------
    if( rank==0 ) then
      WRITE(*,*) ' Solving equation w'''' = v w'' for w'
    end if

    !! set initial values
    s0 = initw1
    s1 = initw2

    !! two initial shots
    call shoot_w(w,v,grid,n,sweep_angle_degrees,kappa,sv,s0,r0,blThick)
    call shoot_w(w,v,grid,n,sweep_angle_degrees,kappa,sv,s1,r1,blThick)

    !! controlled shooting loop
    do

      !!calculate new initial value
      s1old = s1
      s1 = (r0*s1 - r1*s0)/(r0-r1)
      s0 = s1old
      r0 = r1
      !!shoot
      call shoot_w(w,v,grid,n,sweep_angle_degrees,kappa,sv,s1,r1,blThick)

      !!check result
      IF( abs(r1) <= 1E-10 ) THEN
        !!goal reached
        EXIT
      ELSEIF ((s1 .EQ. s0) .OR. (r1 .EQ. r0)) THEN
        !!no more progress
        EXIT
      END IF
    END DO

    !! print output
    if( rank==0 ) then
      !      write(*,*) '   Ended shoot_w successfully. dy(w)|0 = ', w(1,2)
      write(*,*) ' residuum =', r1
      write(*,*)
      !      if (((initw1+initw2)/2).NE.0) then
      !         write(*,*) 'Relative deviation from initial guess:', (w(1,2)-((initw1+initw2)/2))/((initw1+initw2)/2)
      !      else
      !         write(*,*) 'Total deviation from initial guess:   ', w(1,2)
      !      end if
      !      write(*,*)
    end if

    !!! ---------------------------------------------------------- CONVERT INTEGRATION RESULTS TO VERITABLE VELOCITIES --------------------------------------------------
    ub = -v(:,2)
    vb = v(:,1)
    wb = w(:,1)

    !! deallocate local memory
    deallocate(v,w)

  END SUBROUTINE calcbasicflow

  !> \brief init vector field with 2d pulsatile flow in x-direction
  !! \f[ u(y,t) = \hat{u}^+ \exp(i \omega t) + \hat{u}^- \exp(- i \omega t) \f]
  !! \f[ \hat{u}^+(y) = c^+ \left( \exp(+ \lambda_1 y ) + \exp(-\lambda\right) + \frac{p_x i}{\omega \f]
  !! \f[ \hat{u}^-(y) = c^- \left( \exp(+ \lambda_{-1} y ) \right) + \frac{p_x i}{\omega \f]
  subroutine VF_init_2DPulsatileXC(   &
      N,                              &
      bL,bU,                          &
      SU,NU,                          &
      SV,NV,                          &
      SW,NW,                          &
      L2,                             &
      x2p,                            &
      re_, om, px,                    &
      phiU, phiV, phiW ) bind ( c, name='VF_init_2DPulsatileXC' )

    implicit none

    integer(c_int), intent(in)    ::  N(3)

    integer(c_int), intent(in)     :: bL(3)
    integer(c_int), intent(in)     :: bU(3)

    integer(c_int), intent(in)     :: SU(3)
    integer(c_int), intent(in)     :: NU(3)

    integer(c_int), intent(in)     :: SV(3)
    integer(c_int), intent(in)     :: NV(3)

    integer(c_int), intent(in)     :: SW(3)
    integer(c_int), intent(in)     :: NW(3)

    real(c_double), intent(in)     :: L2

    real(c_double), intent(in)     :: x2p( bL(2):(N(2)+bU(2)) )

    real(c_double), intent(in)    ::  re_
    real(c_double), intent(in)    ::  om
    real(c_double), intent(in)    ::  px

    real(c_double),  intent(inout) :: phiU(bL(1):(N(1)+bU(1)),bL(2):(N(2)+bU(2)),bL(3):(N(3)+bU(3)))
    real(c_double),  intent(inout) :: phiV(bL(1):(N(1)+bU(1)),bL(2):(N(2)+bU(2)),bL(3):(N(3)+bU(3)))
    real(c_double),  intent(inout) :: phiW(bL(1):(N(1)+bU(1)),bL(2):(N(2)+bU(2)),bL(3):(N(3)+bU(3)))

    real(c_double) :: pi
    real(c_double) :: Lh
    real(c_double) :: mu
    real(c_double) :: ny
    real(c_double) :: c

    integer                ::  i, j, k

    !--- initial conditions for velocity ---
    ! note: - cf. sketch in file "usr_geometry.f90"
    !
    !         grid points in the domain and on the boundary
    !         |         |         |     velocity component
    !         |         |         |     |
    ! vel(S1U:N1U,S2U:N2U,S3U:N3U,1)
    ! vel(S1V:N1V,S2V:N2V,S3V:N3V,2)
    ! vel(S1W:N1W,S2W:N2W,S3W:N3W,3)
    !
    !  phiU = 0.
    !  phiV = 0.
    !  phiW = 0.
    pi = 4.*atan(1.)
    Lh  = L2/2.
    mu = sqrt( om*re_/2. )*Lh
    c  = px/( om*(cos(mu)**2*cosh(mu)**2 + sin(mu)**2*sinh(mu)**2) )

    do k = SU(3), NU(3)
      do j = SU(2), NU(2)
        do i = SU(1), NU(1)
          ny = sqrt( om*re_/2. )*( x2p(j)-Lh )
          phiU(i,j,k) = -c*( -cos(ny)*cosh(ny)*sin(mu)*sinh(mu) +sin(ny)*sinh(ny)*cos(mu)*cosh(mu) )
        end do
      end do
    end do

    do k = SV(3), NV(3)
      do j = SV(2), NV(2)
        do i = SV(1), NV(1)
          phiV(i,j,k) = 0
        end do
      end do
    end do

    do k = SW(3), NW(3)
      do j = SW(2), NW(2)
        do i = SW(1), NW(1)
          phiW(i,j,k) = 0.0
        end do
      end do
    end do

  end subroutine VF_init_2DPulsatileXC

  !> \brief init vector field with 2d pulsatile flow in x-direction
  !! \f[ u(y,t) = \hat{u}^+ \exp(i \omega t) + \hat{u}^- \exp(- i \omega t) \f]
  !! \f[ \hat{u}^+(y) = c^+ \left( \exp(+ \lambda_1 y ) + \exp(-\lambda\right) + \frac{p_x i}{\omega \f]
  !! \f[ \hat{u}^-(y) = c^- \left( \exp(+ \lambda_{-1} y ) \right) + \frac{p_x i}{\omega \f]
  subroutine VF_init_2DPulsatileYC(   &
      N,                              &
      bL,bU,                          &
      SU,NU,                          &
      SV,NV,                          &
      SW,NW,                          &
      L1,                             &
      x1,                             &
      re_, om_, px,                   &
      phiU, phiV, phiW ) bind ( c, name='VF_init_2DPulsatileYC' )
    ! (basic subroutine)

    implicit none

    integer(c_int), intent(in)    ::  N(3)

    integer(c_int), intent(in)    :: bL(3)
    integer(c_int), intent(in)    :: bU(3)

    integer(c_int), intent(in)    :: SU(3)
    integer(c_int), intent(in)    :: NU(3)

    integer(c_int), intent(in)    :: SV(3)
    integer(c_int), intent(in)    :: NV(3)

    integer(c_int), intent(in)    :: SW(3)
    integer(c_int), intent(in)    :: NW(3)

    real(c_double), intent(in)     :: L1
    real(c_double), intent(in)     :: x1( bL(1):(N(1)+bU(1)) )

    real(c_double), intent(in)    :: re_
    real(c_double), intent(in)    :: om_
    real(c_double), intent(in)    :: px

    real(c_double),  intent(inout) :: phiU(bL(1):(N(1)+bU(1)),bL(2):(N(2)+bU(2)),bL(3):(N(3)+bU(3)))
    real(c_double),  intent(inout) :: phiV(bL(1):(N(1)+bU(1)),bL(2):(N(2)+bU(2)),bL(3):(N(3)+bU(3)))
    real(c_double),  intent(inout) :: phiW(bL(1):(N(1)+bU(1)),bL(2):(N(2)+bU(2)),bL(3):(N(3)+bU(3)))

    real(c_double):: pi
    real(c_double):: Lh
    real(c_double):: mu
    real(c_double):: ny
    real(c_double):: c

    integer(c_int)                ::  i, j, k

    !--- initial conditions for velocity ---
    ! note: - cf. sketch in file "usr_geometry.f90"
    !
    !         grid points in the domain and on the boundary
    !         |         |         |     velocity component
    !         |         |         |     |
    ! vel(S1U:N1U,S2U:N2U,S3U:N3U,1)
    ! vel(S1V:N1V,S2V:N2V,S3V:N3V,2)
    ! vel(S1W:N1W,S2W:N2W,S3W:N3W,3)
    !
    !  phiU = 0.
    !  phiV = 0.
    !  phiW = 0.
    pi = 4.*atan(1.)
    Lh  = L1/2.
    mu = sqrt( om_*re_/2. )*Lh
    c  = px/( om_*(cos(mu)**2*cosh(mu)**2 + sin(mu)**2*sinh(mu)**2) )

    do k = SU(3), NU(3)
      do j = SU(2), NU(2)
        do i = SU(1), NU(1)
          phiU(i,j,k) = 0
        end do
      end do
    end do

    do k = SV(3), NV(3)
      do j = SV(2), NV(2)
        do i = SV(1), NV(1)
          ny = sqrt( om_*re_/2. )*( x1(i)-Lh )
          phiV(i,j,k) = -c*( -cos(ny)*cosh(ny)*sin(mu)*sinh(mu) +sin(ny)*sinh(ny)*cos(mu)*cosh(mu) )
        end do
      end do
    end do

    do k = SW(3), NW(3)
      do j = SW(2), NW(2)
        do i = SW(1), NW(1)
          phiW(i,j,k) = 0.0
        end do
      end do
    end do

  end subroutine VF_init_2DPulsatileYC

  !> \brief init vector field with 2d pulsatile flow in x-direction
  !! \f[ u(y,t) = \hat{u}^+ \exp(i \omega t) + \hat{u}^- \exp(- i \omega t) \f]
  !! \f[ \hat{u}^+(y) = c^+ \left( \exp(+ \lambda_1 y ) + \exp(-\lambda\right) + \frac{p_x i}{\omega \f]
  !! \f[ \hat{u}^-(y) = c^- \left( \exp(+ \lambda_{-1} y ) \right) + \frac{p_x i}{\omega \f]
  subroutine VF_init_2DPulsatileXS(   &
      N,                      &
      bL,bU,                  &
      SU,NU,                  &
      SV,NV,                  &
      SW,NW,                  &
      L2,                     &
      x2,                     &
      re_, om_, px,           &
      phiU, phiV, phiW ) bind ( c, name='VF_init_2DPulsatileXS' )
    ! (basic subroutine)

    implicit none

    integer(c_int), intent(in)    ::  N(3)

    integer(c_int), intent(in)     :: bL(3)
    integer(c_int), intent(in)     :: bU(3)

    integer(c_int), intent(in)     :: SU(3)
    integer(c_int), intent(in)     :: NU(3)

    integer(c_int), intent(in)     :: SV(3)
    integer(c_int), intent(in)     :: NV(3)

    integer(c_int), intent(in)     :: SW(3)
    integer(c_int), intent(in)     :: NW(3)

    real(c_double), intent(in)     :: L2
    real(c_double), intent(in)     :: x2( bL(2):(N(2)+bU(2)) )

    real(c_double), intent(in)    ::  re_
    real(c_double), intent(in)    ::  om_
    real(c_double), intent(in)    ::  px

    real(c_double),  intent(inout) :: phiU(bL(1):(N(1)+bU(1)),bL(2):(N(2)+bU(2)),bL(3):(N(3)+bU(3)))
    real(c_double),  intent(inout) :: phiV(bL(1):(N(1)+bU(1)),bL(2):(N(2)+bU(2)),bL(3):(N(3)+bU(3)))
    real(c_double),  intent(inout) :: phiW(bL(1):(N(1)+bU(1)),bL(2):(N(2)+bU(2)),bL(3):(N(3)+bU(3)))

    real(c_double):: pi
    real(c_double):: Lh
    real(c_double):: mu
    real(c_double):: ny
    real(c_double):: c

    integer(c_int)               ::  i, j, k

    !--- initial conditions for velocity ---
    ! note: - cf. sketch in file "usr_geometry.f90"
    !
    !         grid points in the domain and on the boundary
    !         |         |         |     velocity component
    !         |         |         |     |
    ! vel(S1U:N1U,S2U:N2U,S3U:N3U,1)
    ! vel(S1V:N1V,S2V:N2V,S3V:N3V,2)
    ! vel(S1W:N1W,S2W:N2W,S3W:N3W,3)
    !
    !  phiU = 0.
    !  phiV = 0.
    !  phiW = 0.
    pi = 4.*atan(1.)
    Lh  = L2/2.
    mu = sqrt( om_*re_/2. )*Lh
    c  = px/( om_*(cos(mu)**2*cosh(mu)**2 + sin(mu)**2*sinh(mu)**2) )

    do k = SU(3), NU(3)
      do j = SU(2), NU(2)
        do i = SU(1), NU(1)
          ny = sqrt( om_*Re_/2. )*( x2(j)-Lh )
          phiU(i,j,k) = -c*( cos(ny)*cosh(ny)*cos(mu)*cosh(mu) +sin(ny)*sinh(ny)*sin(mu)*sinh(mu) ) + px/om_
        end do
      end do
    end do

    do k = SV(3), NV(3)
      do j = SV(2), NV(2)
        do i = SV(1), NV(1)
          phiV(i,j,k) = 0
        end do
      end do
    end do

    do k = SW(3), NW(3)
      do j = SW(2), NW(2)
        do i = SW(1), NW(1)
          phiW(i,j,k) = 0.0
        end do
      end do
    end do

  end subroutine VF_init_2DPulsatileXS

  !> \brief init vector field with 2d pulsatile flow in x-direction
  !! \f[ u(y,t) = \hat{u}^+ \exp(i \omega t) + \hat{u}^- \exp(- i \omega t) \f]
  !! \f[ \hat{u}^+(y) = c^+ \left( \exp(+ \lambda_1 y ) + \exp(-\lambda\right) + \frac{p_x i}{\omega \f]
  !! \f[ \hat{u}^-(y) = c^- \left( \exp(+ \lambda_{-1} y ) \right) + \frac{p_x i}{\omega \f]
  subroutine VF_init_2DPulsatileYS(   &
      N,                      &
      bL,bU,                  &
      SU,NU,                  &
      SV,NV,                  &
      SW,NW,                  &
      L1,                     &
      x1,                     &
      re_, om_, px,           &
      phiU,phiV,phiW ) bind ( c, name='VF_init_2DPulsatileYS' )

    implicit none

    integer(c_int), intent(in)    ::  N(3)

    integer(c_int), intent(in)     :: bL(3)
    integer(c_int), intent(in)     :: bU(3)

    integer(c_int), intent(in)     :: SU(3)
    integer(c_int), intent(in)     :: NU(3)

    integer(c_int), intent(in)     :: SV(3)
    integer(c_int), intent(in)     :: NV(3)

    integer(c_int), intent(in)     :: SW(3)
    integer(c_int), intent(in)     :: NW(3)

    real(c_double), intent(in)     :: L1
    real(c_double), intent(in)     :: x1( bL(1):(N(1)+bU(1)) )

    real(c_double), intent(in)    ::  re_
    real(c_double), intent(in)    ::  om_
    real(c_double), intent(in)    ::  px

    real(c_double),  intent(inout) :: phiU(bL(1):(N(1)+bU(1)),bL(2):(N(2)+bU(2)),bL(3):(N(3)+bU(3)))
    real(c_double),  intent(inout) :: phiV(bL(1):(N(1)+bU(1)),bL(2):(N(2)+bU(2)),bL(3):(N(3)+bU(3)))
    real(c_double),  intent(inout) :: phiW(bL(1):(N(1)+bU(1)),bL(2):(N(2)+bU(2)),bL(3):(N(3)+bU(3)))

    real(c_double) :: pi
    real(c_double) :: Lh
    real(c_double) :: mu
    real(c_double) :: ny
    real(c_double) :: c

    integer(c_int)                ::  i, j, k

    !--- initial conditions for velocity ---
    ! note: - cf. sketch in file "usr_geometry.f90"
    !
    !         grid points in the domain and on the boundary
    !         |         |         |     velocity component
    !         |         |         |     |
    ! vel(S1U:N1U,S2U:N2U,S3U:N3U,1)
    ! vel(S1V:N1V,S2V:N2V,S3V:N3V,2)
    ! vel(S1W:N1W,S2W:N2W,S3W:N3W,3)
    !
    !  phiU = 0.
    !  phiV = 0.
    !  phiW = 0.
    pi = 4.*atan(1.)
    Lh  = L1/2.
    mu = sqrt( om_*re_/2. )*Lh
    c  = px/( om_*(cos(mu)**2*cosh(mu)**2 + sin(mu)**2*sinh(mu)**2) )

    do k = SU(3), NU(3)
      do j = SU(2), NU(2)
        do i = SU(1), NU(1)
          phiU(i,j,k) = 0
        end do
      end do
    end do

    do k = SV(3), NV(3)
      do j = SV(2), NV(2)
        do i = SV(1), NV(1)
          ny = sqrt( om_*re_/2. )*( x1(i)-Lh )
          phiV(i,j,k) = -c*( cos(ny)*cosh(ny)*cos(mu)*cosh(mu) +sin(ny)*sinh(ny)*sin(mu)*sinh(mu) ) + px/om_
        end do
      end do
    end do

    do k = SW(3), NW(3)
      do j = SW(2), NW(2)
        do i = SW(1), NW(1)
          phiW(i,j,k) = 0.0
        end do
      end do
    end do

  end subroutine VF_init_2DPulsatileYS

  subroutine VF_init_Vpoint(  &
      N,                      &
      bL,bU,                  &
      SU,NU,                  &
      SV,NV,                  &
      SW,NW,                  &
      L,                      &
      x1u,                    &
      x2p,                    &
      sig,                    &
      phiU,phiV,phiW ) bind ( c, name='VF_init_Vpoint' )
    ! (basic subroutine)

    implicit none

    integer(c_int), intent(in)    ::  N(3)

    integer(c_int), intent(in)     :: bL(3)
    integer(c_int), intent(in)     :: bU(3)

    integer(c_int), intent(in)     :: SU(3)
    integer(c_int), intent(in)     :: NU(3)

    integer(c_int), intent(in)     :: SV(3)
    integer(c_int), intent(in)     :: NV(3)

    integer(c_int), intent(in)     :: SW(3)
    integer(c_int), intent(in)     :: NW(3)

    real(c_double), intent(in)     :: L(3)
    real(c_double), intent(in)     :: x1u( bL(1):(N(1)+bU(1)) )
    real(c_double), intent(in)     :: x2p( bL(2):(N(2)+bU(2)) )

    real(c_double), intent(in)    :: sig

    real(c_double),  intent(inout) :: phiU(bL(1):(N(1)+bU(1)),bL(2):(N(2)+bU(2)),bL(3):(N(3)+bU(3)))
    real(c_double),  intent(inout) :: phiV(bL(1):(N(1)+bU(1)),bL(2):(N(2)+bU(2)),bL(3):(N(3)+bU(3)))
    real(c_double),  intent(inout) :: phiW(bL(1):(N(1)+bU(1)),bL(2):(N(2)+bU(2)),bL(3):(N(3)+bU(3)))

    integer(c_int)               ::  i, j, k

    !--- initial conditions for velocity ---
    ! note: - cf. sketch in file "usr_geometry.f90"
    !
    !         grid points in the domain and on the boundary
    !         |         |         |     velocity component
    !         |         |         |     |
    ! vel(S1U:N1U,S2U:N2U,S3U:N3U,1)
    ! vel(S1V:N1V,S2V:N2V,S3V:N3V,2)
    ! vel(S1W:N1W,S2W:N2W,S3W:N3W,3)
    !

    do k = SU(3), NU(3)
      do j = SU(2), NU(2)
        do i = SU(1), NU(1)
          phiU(i,j,k) = exp( -((x1u(i)-L(1)/2.)/sig)**2 -((x2p(j)-L(2)/2.)/sig/sig)**2 )
        end do
      end do
    end do

    do k = SV(3), NV(3)
      do j = SV(2), NV(2)
        do i = SV(1), NV(1)
          phiV(i,j,k) = 0.
        end do
      end do
    end do

    do k = SW(3), NW(3)
      do j = SW(2), NW(2)
        do i = SW(1), NW(1)
          phiW(i,j,k) = 0.0
        end do
      end do
    end do

  end subroutine VF_init_Vpoint

  subroutine VF_init_RankineVortex(   &
      N,                              &
      bL,bU,                          &
      SU,NU,                          &
      SV,NV,                          &
      SW,NW,                          &
      L,                              &
      x1p,x2p,                        &
      x1u,x2v,                        &
      phiU,phiV,phiW ) bind ( c, name='VF_init_RankineVortex' )
    ! (basic subroutine)

    implicit none

    integer(c_int), intent(in)    ::  N(3)

    integer(c_int), intent(in)     :: bL(3)
    integer(c_int), intent(in)     :: bU(3)

    integer(c_int), intent(in)     :: SU(3)
    integer(c_int), intent(in)     :: NU(3)

    integer(c_int), intent(in)     :: SV(3)
    integer(c_int), intent(in)     :: NV(3)

    integer(c_int), intent(in)     :: SW(3)
    integer(c_int), intent(in)     :: NW(3)

    real(c_double), intent(in)     :: L(3)

    real(c_double), intent(in)     :: x1p( bL(1):(N(1)+bU(1)) )
    real(c_double), intent(in)     :: x2p( bL(2):(N(2)+bU(2)) )

    real(c_double), intent(in)     :: x1u( bL(1):(N(1)+bU(1)) )
    real(c_double), intent(in)     :: x2v( bL(2):(N(2)+bU(2)) )

    real(c_double),  intent(inout) :: phiU(bL(1):(N(1)+bU(1)),bL(2):(N(2)+bU(2)),bL(3):(N(3)+bU(3)))
    real(c_double),  intent(inout) :: phiV(bL(1):(N(1)+bU(1)),bL(2):(N(2)+bU(2)),bL(3):(N(3)+bU(3)))
    real(c_double),  intent(inout) :: phiW(bL(1):(N(1)+bU(1)),bL(2):(N(2)+bU(2)),bL(3):(N(3)+bU(3)))

    real(c_double):: circ
    real(c_double):: rad
    real(c_double):: r
    real(c_double):: pi

    integer(c_int)               ::  i, j, k

    !--- initial conditions for velocity ---
    ! note: - cf. sketch in file "usr_geometry.f90"
    !
    !         grid points in the domain and on the boundary
    !         |         |         |     velocity component
    !         |         |         |     |
    ! vel(S1U:N1U,S2U:N2U,S3U:N3U,1)
    ! vel(S1V:N1V,S2V:N2V,S3V:N3V,2)
    ! vel(S1W:N1W,S2W:N2W,S3W:N3W,3)
    !
    pi = 4.*atan(1.)
    rad = L(1)/2./2.
    circ = 2.*pi*rad

    do k = SU(3), NU(3)
      do j = SU(2), NU(2)
        do i = SU(1), NU(1)
          r = sqrt( (x1u(i)-L(1)/2)**2 + (x2p(j)-L(2)/2)**2 )
          if( r<= rad ) then
            phiU(i,j,k) = -(x2p(j)-L(2)/2)/rad
          else
            phiU(i,j,k) = -(x2p(j)-L(2)/2)*rad/r/r
          endif
        end do
      end do
    end do

    do k = SV(3), NV(3)
      do j = SV(2), NV(2)
        do i = SV(1), NV(1)
          r = sqrt( (x1p(i)-L(1)/2)**2 + (x2v(j)-L(2)/2)**2 )
          if( r<=rad ) then
            phiV(i,j,k) = (x1p(i)-L(1)/2)/rad
          else
            phiV(i,j,k) = (x1p(i)-L(1)/2)*rad/r/r
          endif
        end do
      end do
    end do

    do k = SW(3), NW(3)
      do j = SW(2), NW(2)
        do i = SW(1), NW(1)
          phiW(i,j,k) = 0.0
        end do
      end do
    end do

  end subroutine VF_init_RankineVortex

  subroutine VF_init_GaussianForcing1D(   &
      N,                                  &
      bL,bU,                              &
      SU,NU,                              &
      SV,NV,                              &
      SW,NW,                              &
      L1,                                 &
      x1u,                                &
      phiU,phiV,phiW ) bind ( c, name='VF_init_GaussianForcing1D' )
    ! (basic subroutine)

    implicit none

    integer(c_int), intent(in)    ::  N(3)

    integer(c_int), intent(in)     :: bL(3)
    integer(c_int), intent(in)     :: bU(3)

    integer(c_int), intent(in)     :: SU(3)
    integer(c_int), intent(in)     :: NU(3)

    integer(c_int), intent(in)     :: SV(3)
    integer(c_int), intent(in)     :: NV(3)

    integer(c_int), intent(in)     :: SW(3)
    integer(c_int), intent(in)     :: NW(3)

    real(c_double), intent(in)     :: L1
    real(c_double), intent(in)     :: x1u( bL(1):(N(1)+bU(1)) )

    real(c_double),  intent(inout) :: phiU(bL(1):(N(1)+bU(1)),bL(2):(N(2)+bU(2)),bL(3):(N(3)+bU(3)))
    real(c_double),  intent(inout) :: phiV(bL(1):(N(1)+bU(1)),bL(2):(N(2)+bU(2)),bL(3):(N(3)+bU(3)))
    real(c_double),  intent(inout) :: phiW(bL(1):(N(1)+bU(1)),bL(2):(N(2)+bU(2)),bL(3):(N(3)+bU(3)))

    real(c_double):: sig

    integer(c_int)               ::  i, j, k

    !--- initial conditions for velocity ---
    ! note: - cf. sketch in file "usr_geometry.f90"
    !
    !         grid points in the domain and on the boundary
    !         |         |         |     velocity component
    !         |         |         |     |
    ! vel(S1U:N1U,S2U:N2U,S3U:N3U,1)
    ! vel(S1V:N1V,S2V:N2V,S3V:N3V,2)
    ! vel(S1W:N1W,S2W:N2W,S3W:N3W,3)
    !

    sig = 0.2

    do k = SU(3), NU(3)
      do j = SU(2), NU(2)
        do i = SU(1), NU(1)
          phiU(i,j,k) = exp( -((x1u(i)-L1/2)/sig)**2  )
        end do
      end do
    end do

    do k = SV(3), NV(3)
      do j = SV(2), NV(2)
        do i = SV(1), NV(1)
          phiV(i,j,k) = 0
        end do
      end do
    end do

    do k = SW(3), NW(3)
      do j = SW(2), NW(2)
        do i = SW(1), NW(1)
          phiW(i,j,k) = 0.0
        end do
      end do
    end do

  end subroutine VF_init_GaussianForcing1D

  subroutine VF_init_GaussianForcing2D(   &
      N,                                  &
      bL,bU,                              &
      SU,NU,                              &
      SV,NV,                              &
      SW,NW,                              &
      L,                                  &
      x1p,x2p,                            &
      x1u,x2v,                            &
      phiU,phiV,phiW ) bind ( c, name='VF_init_GaussianForcing2D' )
    ! (basic subroutine)

    implicit none

    integer(c_int), intent(in)    ::  N(3)

    integer(c_int), intent(in)     :: bL(3)
    integer(c_int), intent(in)     :: bU(3)

    integer(c_int), intent(in)     :: SU(3)
    integer(c_int), intent(in)     :: NU(3)

    integer(c_int), intent(in)     :: SV(3)
    integer(c_int), intent(in)     :: NV(3)

    integer(c_int), intent(in)     :: SW(3)
    integer(c_int), intent(in)     :: NW(3)

    real(c_double), intent(in)     :: L(3)

    real(c_double), intent(in)     :: x1p( bL(1):(N(1)+bU(1)) )
    real(c_double), intent(in)     :: x2p( bL(2):(N(2)+bU(2)) )

    real(c_double), intent(in)     :: x1u( bL(1):(N(1)+bU(1)) )
    real(c_double), intent(in)     :: x2v( bL(2):(N(2)+bU(2)) )

    real(c_double),  intent(inout) :: phiU(bL(1):(N(1)+bU(1)),bL(2):(N(2)+bU(2)),bL(3):(N(3)+bU(3)))
    real(c_double),  intent(inout) :: phiV(bL(1):(N(1)+bU(1)),bL(2):(N(2)+bU(2)),bL(3):(N(3)+bU(3)))
    real(c_double),  intent(inout) :: phiW(bL(1):(N(1)+bU(1)),bL(2):(N(2)+bU(2)),bL(3):(N(3)+bU(3)))
    real(c_double):: sig

    integer(c_int)               ::  i, j, k

    !--- initial conditions for velocity ---
    ! note: - cf. sketch in file "usr_geometry.f90"
    !
    !         grid points in the domain and on the boundary
    !         |         |         |     velocity component
    !         |         |         |     |
    ! vel(S1U:N1U,S2U:N2U,S3U:N3U,1)
    ! vel(S1V:N1V,S2V:N2V,S3V:N3V,2)
    ! vel(S1W:N1W,S2W:N2W,S3W:N3W,3)
    !

    sig = 0.2

    do k = SU(3), NU(3)
      do j = SU(2), NU(2)
        do i = SU(1), NU(1)
          phiU(i,j,k) = exp( -((x1u(i)   )/sig)**2 - ((x2p(j)   )/sig)**2 ) / sqrt(2.)    &
            + exp( -((x1u(i)-L(1))/sig)**2 - ((x2p(j)-L(2))/sig)**2 ) / sqrt(2.)    &
            + exp( -((x1u(i)-L(1))/sig)**2 - ((x2p(j)   )/sig)**2 ) / sqrt(2.)    &
            + exp( -((x1u(i)   )/sig)**2 - ((x2p(j)-L(2))/sig)**2 ) / sqrt(2.)
        end do
      end do
    end do

    do k = SV(3), NV(3)
      do j = SV(2), NV(2)
        do i = SV(1), NV(1)
          phiV(i,j,k) = exp( -((x1p(i)   )/sig)**2 -((x2v(j)   )/sig)**2 ) / sqrt(2.) &
            + exp( -((x1p(i)-L(1))/sig)**2 -((x2v(j)-L(2))/sig)**2 ) / sqrt(2.) &
            + exp( -((x1p(i)-L(1))/sig)**2 -((x2v(j)   )/sig)**2 ) / sqrt(2.) &
            + exp( -((x1p(i)   )/sig)**2 -((x2v(j)-L(2))/sig)**2 ) / sqrt(2.)
        end do
      end do
    end do

    do k = SW(3), NW(3)
      do j = SW(2), NW(2)
        do i = SW(1), NW(1)
          phiW(i,j,k) = 0.0
        end do
      end do
    end do

  end subroutine VF_init_GaussianForcing2D

  subroutine VF_init_BoundaryFilter1D(    &
      N,                                  &
      bL,bU,                              &
      SU,NU,                              &
      SV,NV,                              &
      SW,NW,                              &
      L1,                                 &
      x1u,                                &
      phiU,phiV,phiW ) bind ( c, name='VF_init_BoundaryFilter1D' )
    ! (basic subroutine)

    implicit none

    integer(c_int), intent(in)    ::  N(3)

    integer(c_int), intent(in)     :: bL(3)
    integer(c_int), intent(in)     :: bU(3)

    integer(c_int), intent(in)     :: SU(3)
    integer(c_int), intent(in)     :: NU(3)

    integer(c_int), intent(in)     :: SV(3)
    integer(c_int), intent(in)     :: NV(3)

    integer(c_int), intent(in)     :: SW(3)
    integer(c_int), intent(in)     :: NW(3)

    real(c_double), intent(in)     :: L1

    real(c_double), intent(in)     :: x1u( bL(1):(N(1)+bU(1)) )

    real(c_double),  intent(inout) :: phiU(bL(1):(N(1)+bU(1)),bL(2):(N(2)+bU(2)),bL(3):(N(3)+bU(3)))
    real(c_double),  intent(inout) :: phiV(bL(1):(N(1)+bU(1)),bL(2):(N(2)+bU(2)),bL(3):(N(3)+bU(3)))
    real(c_double),  intent(inout) :: phiW(bL(1):(N(1)+bU(1)),bL(2):(N(2)+bU(2)),bL(3):(N(3)+bU(3)))

    real(c_double):: h
    real(c_double):: h2

    integer(c_int)                ::  i, j, k

    !--- initial conditions for velocity ---
    ! note: - cf. sketch in file "usr_geometry.f90"
    !
    !         grid points in the domain and on the boundary
    !         |         |         |     velocity component
    !         |         |         |     |
    ! vel(S1U:N1U,S2U:N2U,S3U:N3U,1)
    ! vel(S1V:N1V,S2V:N2V,S3V:N3V,2)
    ! vel(S1W:N1W,S2W:N2W,S3W:N3W,3)
    !

    h = 0.1
    h2 = h*h

    do k = SU(3), NU(3)
      do j = SU(2), NU(2)
        do i = SU(1), NU(1)
          phiU(i,j,k) = 10*( exp( -(x1u(i)**2)/h2  ) + exp( -((x1u(i)-L1)**2)/h2  ) )
        end do
      end do
    end do

    do k = SV(3), NV(3)
      do j = SV(2), NV(2)
        do i = SV(1), NV(1)
          phiV(i,j,k) = 0
        end do
      end do
    end do

    do k = SW(3), NW(3)
      do j = SW(2), NW(2)
        do i = SW(1), NW(1)
          phiW(i,j,k) = 0.0
        end do
      end do
    end do

  end subroutine VF_init_BoundaryFilter1D

  subroutine VF_init_BoundaryFilter2D(    &
      N,                                  &
      bL,bU,                              &
      SU,NU,                              &
      SV,NV,                              &
      SW,NW,                              &
      L,                                  &
      x1p,x2p,                            &
      x1u,x2v,                            &
      phiU,phiV,phiW ) bind ( c, name='VF_init_BoundaryFilter2D' )
    ! (basic subroutine)

    implicit none

    integer(c_int), intent(in)    ::  N(3)

    integer(c_int), intent(in)     :: bL(3)
    integer(c_int), intent(in)     :: bU(3)

    integer(c_int), intent(in)     :: SU(3)
    integer(c_int), intent(in)     :: NU(3)

    integer(c_int), intent(in)     :: SV(3)
    integer(c_int), intent(in)     :: NV(3)

    integer(c_int), intent(in)     :: SW(3)
    integer(c_int), intent(in)     :: NW(3)

    real(c_double), intent(in)     :: L(3)

    real(c_double), intent(in)     :: x1p( bL(1):(N(1)+bU(1)) )
    real(c_double), intent(in)     :: x2p( bL(2):(N(2)+bU(2)) )

    real(c_double), intent(in)     :: x1u( bL(1):(N(1)+bU(1)) )
    real(c_double), intent(in)     :: x2v( bL(2):(N(2)+bU(2)) )

    real(c_double),  intent(inout) :: phiU(bL(1):(N(1)+bU(1)),bL(2):(N(2)+bU(2)),bL(3):(N(3)+bU(3)))
    real(c_double),  intent(inout) :: phiV(bL(1):(N(1)+bU(1)),bL(2):(N(2)+bU(2)),bL(3):(N(3)+bU(3)))
    real(c_double),  intent(inout) :: phiW(bL(1):(N(1)+bU(1)),bL(2):(N(2)+bU(2)),bL(3):(N(3)+bU(3)))

    real(c_double):: h
    real(c_double):: h2

    integer(c_int)               ::  i, j, k

    !--- initial conditions for velocity ---
    ! note: - cf. sketch in file "usr_geometry.f90"
    !
    !         grid points in the domain and on the boundary
    !         |         |         |     velocity component
    !         |         |         |     |
    ! vel(S1U:N1U,S2U:N2U,S3U:N3U,1)
    ! vel(S1V:N1V,S2V:N2V,S3V:N3V,2)
    ! vel(S1W:N1W,S2W:N2W,S3W:N3W,3)
    !

    h = 0.1
    h2 = h*h

    do k = SU(3), NU(3)
      do j = SU(2), NU(2)
        do i = SU(1), NU(1)
          phiU(i,j,k) = max( 10*( exp( -(x1u(i)**2)/h2  ) + exp( -((x1u(i)-L(1))**2)/h2  ) ) , &
            10*( exp( -(x2p(j)**2)/h2  ) + exp( -((x2p(j)-L(2))**2)/h2  ) ) )
        end do
      end do
    end do

    do k = SV(3), NV(3)
      do j = SV(2), NV(2)
        do i = SV(1), NV(1)
          phiV(i,j,k) = max( 10*( exp( -(x1p(i)**2)/h2  ) + exp( -((x1p(i)-L(1))**2)/h2  ) ) , &
            10*( exp( -(x2v(j)**2)/h2  ) + exp( -((x2v(j)-L(2))**2)/h2  ) ) )
        end do
      end do
    end do

    do k = SW(3), NW(3)
      do j = SW(2), NW(2)
        do i = SW(1), NW(1)
          phiW(i,j,k) = 0.0
        end do
      end do
    end do

  end subroutine VF_init_BoundaryFilter2D

  !> \f$ u = \min( 4*\exp( -((x-xm)/rad)**2 -((y-ym)/rad)**2 ),1.) \f$
  subroutine VF_init_Disc(    &
      N,                      &
      bL,bU,                  &
      SU,NU,                  &
      SV,NV,                  &
      SW,NW,                  &
      x1p,x2p,x3p,            &
      x1u,x2v,                &
      xm,ym, rad,sca,         &
      phiU,phiV,phiW ) bind ( c, name='VF_init_Disc' )

    implicit none

    integer(c_int), intent(in)    ::  N(3)

    integer(c_int), intent(in)     :: bL(3)
    integer(c_int), intent(in)     :: bU(3)

    integer(c_int), intent(in)     :: SU(3)
    integer(c_int), intent(in)     :: NU(3)

    integer(c_int), intent(in)     :: SV(3)
    integer(c_int), intent(in)     :: NV(3)

    integer(c_int), intent(in)     :: SW(3)
    integer(c_int), intent(in)     :: NW(3)

    real(c_double), intent(in)     :: x1p( bL(1):(N(1)+bU(1)) )
    real(c_double), intent(in)     :: x2p( bL(2):(N(2)+bU(2)) )
    real(c_double), intent(in)     :: x3p( bL(3):(N(3)+bU(3)) )

    real(c_double), intent(in)     :: x1u( bL(1):(N(1)+bU(1)) )
    real(c_double), intent(in)     :: x2v( bL(2):(N(2)+bU(2)) )

    real(c_double),  intent(inout) :: phiU(bL(1):(N(1)+bU(1)),bL(2):(N(2)+bU(2)),bL(3):(N(3)+bU(3)))
    real(c_double),  intent(inout) :: phiV(bL(1):(N(1)+bU(1)),bL(2):(N(2)+bU(2)),bL(3):(N(3)+bU(3)))
    real(c_double),  intent(inout) :: phiW(bL(1):(N(1)+bU(1)),bL(2):(N(2)+bU(2)),bL(3):(N(3)+bU(3)))

    real(c_double), intent(in)    :: xm
    real(c_double), intent(in)    :: ym
    real(c_double), intent(in)    :: rad
    real(c_double), intent(in)    :: sca

    real(c_double)   :: dr

    integer(c_int) ::  i, j, k

    !--- initial conditions for velocity ---
    ! note: - cf. sketch in file "usr_geometry.f90"
    !
    !         grid points in the domain and on the boundary
    !         |         |         |     velocity component
    !         |         |         |     |
    ! vel(S1U:N1U,S2U:N2U,S3U:N3U,1)
    ! vel(S1V:N1V,S2V:N2V,S3V:N3V,2)
    ! vel(S1W:N1W,S2W:N2W,S3W:N3W,3)
    !

    !        dr = SQRT( ( 2*L1/REAL(M1-1) )**2 + ( 2*L2/REAL(M2-1) )**2 )
    dr = SQRT( ( 2*(x1p(2)-x1p(1)) )**2 + ( 2*(x2p(2)-x2p(1)) )**2 )

    do k = SU(3), NU(3)
      do j = SU(2), NU(2)
        do i = SU(1), NU(1)
          phiU(i,j,k) = sca*distance2ib( x1u(i),x2p(j),x3p(k),xm,ym,rad,dr )
        end do
      end do
    end do

    do k = SV(3), NV(3)
      do j = SV(2), NV(2)
        do i = SV(1), NV(1)
          phiV(i,j,k) = sca*distance2ib( x1p(i),x2v(j),x3p(k),xm,ym,rad,dr )
        end do
      end do
    end do

    do k = SW(3), NW(3)
      do j = SW(2), NW(2)
        do i = SW(1), NW(1)
          phiW(i,j,k) = 0.0
        end do
      end do
    end do

  end subroutine VF_init_Disc

  subroutine VF_init_RotatingDisc(    &
      N,                              &
      bL,bU,                          &
      SU,NU,                          &
      SV,NV,                          &
      SW,NW,                          &
      x1p,x2p,                        &
      xm,ym, omega,                   &
      phiU,                           &
      phiV,                           &
      phiW ) bind ( c, name='VF_init_RotatingDisc' )

    implicit none

    integer(c_int), intent(in)    ::  N(3)

    integer(c_int), intent(in)     :: bL(3)
    integer(c_int), intent(in)     :: bU(3)

    integer(c_int), intent(in)     :: SU(3)
    integer(c_int), intent(in)     :: NU(3)

    integer(c_int), intent(in)     :: SV(3)
    integer(c_int), intent(in)     :: NV(3)

    integer(c_int), intent(in)     :: SW(3)
    integer(c_int), intent(in)     :: NW(3)

    real(c_double), intent(in)     :: x1p( bL(1):(N(1)+bU(1)) )
    real(c_double), intent(in)     :: x2p( bL(2):(N(2)+bU(2)) )

    real(c_double), intent(in)    :: xm
    real(c_double), intent(in)    :: ym
    real(c_double), intent(in)    :: omega

    real(c_double),  intent(inout) :: phiU(bL(1):(N(1)+bU(1)),bL(2):(N(2)+bU(2)),bL(3):(N(3)+bU(3)))
    real(c_double),  intent(inout) :: phiV(bL(1):(N(1)+bU(1)),bL(2):(N(2)+bU(2)),bL(3):(N(3)+bU(3)))
    real(c_double),  intent(inout) :: phiW(bL(1):(N(1)+bU(1)),bL(2):(N(2)+bU(2)),bL(3):(N(3)+bU(3)))

    integer(c_int)               ::  i, j, k

    !--- initial conditions for velocity ---
    ! note: - cf. sketch in file "usr_geometry.f90"
    !
    !         grid points in the domain and on the boundary
    !         |         |         |     velocity component
    !         |         |         |     |
    ! vel(S1U:N1U,S2U:N2U,S3U:N3U,1)
    ! vel(S1V:N1V,S2V:N2V,S3V:N3V,2)
    ! vel(S1W:N1W,S2W:N2W,S3W:N3W,3)
    !

    do k = SU(3), NU(3)
      do j = SU(2), NU(2)
        do i = SU(1), NU(1)
          phiU(i,j,k) = -omega*(x2p(j)-ym)
        end do
      end do
    end do

    do k = SV(3), NV(3)
      do j = SV(2), NV(2)
        do i = SV(1), NV(1)
          phiV(i,j,k) = omega*(x1p(i)-xm)
        end do
      end do
    end do

    do k = SW(3), NW(3)
      do j = SW(2), NW(2)
        do i = SW(1), NW(1)
          phiW(i,j,k) = 0.0
        end do
      end do
    end do

  end subroutine VF_init_RotatingDisc

  !> \brief init swept Hiemenz boundary layer
  !!
  !! \param rank
  !! \param iShift
  !! \param IB1
  !! \param M
  !! \param N
  !! \param bL
  !! \param bU
  !! \param dL
  !! \param dU
  !! \param SU
  !! \param NU
  !! \param SV
  !! \param NV
  !! \param SW
  !! \param NW
  !! \param y1p
  !! \param y1u
  !! \param x3w
  !! \param cIup
  !! \param Re
  !! \param nonDim
  !! \param kappa Properties of swept Hiemenz flow
  !! \param sweep_angle_degrees Properties of swept Hiemenz flow
  !! \param sweep_angle Properties of swept Hiemenz flow
  !! \param velu
  !! \param velv
  !! \param velw
  !!
  !! \note - baseflow_global(*,1) need to run from 0 to M1, i.e. full axis, on u
  !!         grid, since interpolated to p grid by init_BC
  !!       - baseflow_global(*,2) and (*,3) are computed on p grid, since they
  !!         start exactly on p(1) wall !
  !! difficulty: - full u grid serves also negative values below the wall, where shooting integration is invalid
  !! solution:   - compute velocity for y = 0 and y > 0 on u grid, then extrapolate from 0 to first (negative) u value
  !!             - perform extrapolation not on all ranks, but only those which touch the ground
  !!             - since only one value below the wall is relevant for
  !!               interpolation, extrapolation formula reads:
  !!               base(0,1)|_u = ( base(0,1)|_p - sum_0^d1U{cIup(j,1)*base(1+j,1)|_u} ) /
  !!               cIup(-1,1)
  subroutine VF_init_SHBF(  &
      rank,                 &
      iShift,               &
      IB1,                  &
      M,                    &
      N,                    &
      bL,bU,                &
      dL,dU,                &
      SU,NU,                &
      SV,NV,                &
      SW,NW,                &
      y1p,                  &
      y1u,                  &
      x3w,                  &
      cIup,                 &
      Re,                   &
      nonDim,               &
      kappa,                &
      sweep_angle_degrees,  &
      velU,                 &
      velV,                 &
      velW ) bind ( c, name='VF_init_SHBF' )

    implicit none

    integer(c_int), intent(in)     :: rank
    integer(c_int), intent(in)     :: iShift
    integer(c_int), intent(in)     :: IB1

    integer(c_int), intent(in)     :: M(3)
    integer(c_int), intent(in)     :: N(3)

    integer(c_int), intent(in)     :: bL(3)
    integer(c_int), intent(in)     :: bU(3)

    integer(c_int), intent(in)     :: dL(3)
    integer(c_int), intent(in)     :: dU(3)

    integer(c_int), intent(in)     :: SU(3)
    integer(c_int), intent(in)     :: NU(3)

    integer(c_int), intent(in)     :: SV(3)
    integer(c_int), intent(in)     :: NV(3)

    integer(c_int), intent(in)     :: SW(3)
    integer(c_int), intent(in)     :: NW(3)

    real(c_double), intent(in)     :: y1p( 1:M(1) )

    real(c_double), intent(in)     :: y1u( 0:M(1) )

    real(c_double), intent(in)     :: x3w( bL(3):(N(3)+bU(3)) )

    !real(c_double), intent(in)     :: cIup( dU(1):dL(1), 0:N(1) )
    real(c_double), intent(in)     :: cIup( dU(1):dL(1) )

    real(c_double), intent(in)     :: Re
    integer(c_int), intent(in)     :: nonDim
    real(c_double), intent(in)     :: kappa
    real(c_double), intent(in)     :: sweep_angle_degrees

    real(c_double),  intent(inout) :: velU(bL(1):(N(1)+bU(1)),bL(2):(N(2)+bU(2)),bL(3):(N(3)+bU(3)))
    real(c_double),  intent(inout) :: velV(bL(1):(N(1)+bU(1)),bL(2):(N(2)+bU(2)),bL(3):(N(3)+bU(3)))
    real(c_double),  intent(inout) :: velW(bL(1):(N(1)+bU(1)),bL(2):(N(2)+bU(2)),bL(3):(N(3)+bU(3)))

    integer(c_int)                ::  i, ii, j, k

    !! mjohn 051012 - variables required for changing the nondimensionalization
    real(c_double)                ::  blThick
    !! mjohn 051012

    real(c_double)                ::  y1u_temp(0:M(1))

    real(c_double)                ::  baseflow_global(bL(1):(M(1)+bU(1)),1:3)
    real(c_double)                ::  baseflow       (bL(1):(N(1)+bU(1)),1:3)
    !real(c_double)                ::  baseflow1p     (bL(1):(N(1)+bU(1)))

    real(c_double) :: pi

    pi = 4.*atan(1.)    !!set constants

    !---- mjohn 120207 --------------------------- SWEPT HIEMENZ BASE FLOW ---------------------

    ! initialize variables
    velU = 0.
    velV = 0.
    velW = 0.

    baseflow_global = 0.
    baseflow        = 0.

    ! ------------------------- get two tangential base flow components (staggered grid! need information on p grid)
    ! get global basic flow (exactly FULL y axis, over all ranks) on yp grid, store in baseflow_global, for two tangent coordinates
    call calcbasicflow( rank, kappa, sweep_angle_degrees, M(1),y1p(1),baseflow_global(1,3),baseflow_global(1,1),baseflow_global(1,2), blThick)
    ! --- mjohn 051012
    ! set different velocity fields if integration with Hiemenz is performed
    ! integration routines however are unaltered (see below), because variables are overwritten (see usr_config.f90)
    select case (nonDim)
    case(0)
      if( rank == 0 ) then
        write(*,'(a)') 'classical non-dimensionalization applied'
        write(*,*) 'Kappa and Phi are', kappa, ', ', sweep_angle_degrees, '.'
      end if
      ! traditinoal SHBL case
      baseflow(bL(1):(N(1)+bU(1)),2) = baseflow_global((bL(1)+iShift):(N(1)+bU(1)+iShift),2)
    case default
      if( rank == 0 ) then
        write(*,'(a)') 'novel non-dimensionalization applied'
        write(*,*) 'Kappa and Phi are', kappa, ', ', sweep_angle_degrees, '.'
      end if
      ! 2 degrees of freedom and novel formalism (cases 1 and 2)
      baseflow(bL(1):(N(1)+bU(1)),2) = baseflow_global((bL(1)+iShift):(N(1)+bU(1)+iShift),2) * sin(sweep_angle_degrees*pi/180.)
    end select
    baseflow(bL(1):(N(1)+bU(1)),3) = baseflow_global((bL(1)+iShift):(N(1)+bU(1)+iShift),3) / Re

    ! ------------------------- get wall-normal base flow component (staggered grid! need information on u grid)
    ! get global basic flow (ENTIRE y axis, excluding points beyond boundaries, over all ranks) on u grid, store in baseflow_global, for wall-normal coordinate
    y1u_temp(1:M(1))  = y1u(1:M(1))
    y1u_temp(0     )  = 0.

    call calcbasicflow(rank, kappa, sweep_angle_degrees, M(1)+1,y1u_temp(0),baseflow_global(0,3),baseflow_global(0,1),baseflow_global(0,2), blThick)

    !write(*,*)  'ii, cIup(ii,i)'
    !do ii = dL(1), dU(1)
    !write(*,*)  ii, cIup(ii)
    !end do

    if( 1==IB1 ) then ! change to BCL(1)
      do ii = 0, dU(1)
        !write(*,*) ii, cIup(ii),baseflow_global(1+ii,1)
        baseflow_global(0,1) = baseflow_global(0,1) - cIup(ii)*baseflow_global(1+ii,1)
      end do
      baseflow_global(0,1) = baseflow_global(0,1) / cIup(-1)

    end if
    !write(*,*) cIup(-1),baseflow_global(0,1)

    baseflow(bL(1):(N(1)+bU(1)),1) = baseflow_global((bL(1)+iShift):(N(1)+bU(1)+iShift),1) / Re

    !! TEST!!! - debugging purposes only
    if( rank == 0 ) then
      !write(*,*) "distance, velocity (both in x1-direction)"
      !DO i = SU(1), NU(1)
        !write(*,*)  i+iShift, y1u(i+iShift), baseflow(i,1)
      !end do
      !write(*,*)
      !write(*,*) "distance, velocity (both in x2-direction)"
      !DO i = SV(1), NV(1)
        !write(*,*) i+iShift, y1u(i+iShift), baseflow(i,2)
      !end do
      !write(*,*)
      !write(*,*) "distance, velocity (both in x3-direction)"
      !DO i = SW(1), NW(1)
        !write(*,*) i+iShift, y1u(i+iShift), baseflow(i,3)
      !end do
      write(*,*) 'blThick', blThick
    endif

    do k = SU(3), NU(3)
      do j = SU(2), NU(2)
        do i = SU(1), NU(1)
          velU(i,j,k) = velU(i,j,k) + baseflow(i,1)
        END DO
      END DO
    END DO
    do k = SV(3), NV(3)
      do j = SV(2), NV(2)
        do i = SV(1), NV(1)
          velV(i,j,k) = velV(i,j,k) + baseflow(i,2)
        END DO
      END DO
    END DO
    do k = SW(3), NW(3)
      do j = SW(2), NW(2)
        do i = SW(1), NW(1)
          velW(i,j,k) = velW(i,j,k) + baseflow(i,3)*x3w(k)
        END DO
      END DO
    END DO

  end subroutine VF_init_SHBF

  function eddy(x1,x2,dir) result(fn_val)

    implicit none

    real(c_double), intent(in   ) ::  x1
    real(c_double), intent(in   ) ::  x2
    integer(c_int), intent(in   ) ::  dir

    real(c_double)                ::  length, orient, pi
    real(c_double)                ::  fn_val

    pi = 4.*atan(1.)    !!set constants

    fn_val = 0.

    if (1 == 2) then

      length = sqrt(x1**2+x2**2)
      orient = atan2(x2,x1)

      fn_val = sin(pi*erf(1.*length))

      !if (dir == 1) fn_val = -fn_val*sin(orient)
      !if (dir == 2) fn_val =  fn_val*cos(orient)
      !if (dir == 3) fn_val =  fn_val*cos(orient)
      if (dir == 1) fn_val =  fn_val*sin(orient)
      if (dir == 2) fn_val = -fn_val*cos(orient)
      if (dir == 3) fn_val = -fn_val*cos(orient)

    else

      if (x1 >= -1. .and. x1 <= 1. .and. x2 >= -1. .and. x2 <= 1.) then
        !! if (dir == 1) fn_val = -sin(pi*x1)*(1.+cos(pi*x2))/2. ! "quatropol" mit sich teilender instabilitaet (eine verbleibt in der symmetriefläche)
        !! if (dir == 2) fn_val =  sin(pi*x2)*(1.+cos(pi*x1))/2.
        !! if (dir == 3) fn_val =  sin(pi*x2)*(1.+cos(pi*x1))/2.
        ! if (dir == 1) fn_val = -sin(pi*x2)*(1.+cos(pi*x1))/2.
        ! if (dir == 2) fn_val =  sin(pi*x1)*(1.+cos(pi*x2))/2.
        ! if (dir == 3) fn_val =  sin(pi*x1)*(1.+cos(pi*x2))/2.
        if (dir == 1) fn_val =  sin(pi*x2)*(1.+cos(pi*x1))/2.
        if (dir == 2) fn_val = -sin(pi*x1)*(1.+cos(pi*x2))/2.
        if (dir == 3) fn_val = -sin(pi*x1)*(1.+cos(pi*x2))/2.
      else
        fn_val = 0.
      end if

      !if (x1 < -1. .or. x1 > 1.) fn_val = 0.
      !if (x2 < -1. .or. x2 > 1.) fn_val = 0.

    end if

    return

  end function eddy

  !> \brief init disturbance
  !!
  !!  0 generic GH
  !!  1 counter-rotating primary vortices, even secondary vortices
  !!  2 counter-rotating primary vortices, odd secondary vortices
  !!  3 pair of antisymmetric streaks
  !!  4 array of periodic vortices
  !!  5 line of exaggerated spanwise velocity
  !!  6 Goertler-Haemmerlin disturbance in w component (spanwise)
  !!  7 generic noise (OU process)
  !!  8 pair of counter-rotating vortices with d/dz (init cond)
  !!  9 array of counter-rotating vortices with d/dz (init cond)
  subroutine VF_init_Dist(  &
      rank,                 &
      N,                    &
      bL,bU,                &
      SU,NU,                &
      SV,NV,                &
      SW,NW,                &
      BC_3L_global,         &
      x1u,                  &
      x1p,                  &
      x2p,                  &
      x3w,                  &
      x3p,                  &
      dist_type,            &
      vortex_ampli_prim,    &
      vortex_x1pos,         &
      vortex_x3pos,         &
      vortex_radius,        &
      vortex_band,          &
      velU,velV,velW ) bind ( c, name='VF_init_Dist' )

    implicit none

    integer(c_int), intent(in)     :: rank

    integer(c_int), intent(in)     :: N(3)

    integer(c_int), intent(in)     :: bL(3)
    integer(c_int), intent(in)     :: bU(3)

    integer(c_int), intent(in)     :: SU(3)
    integer(c_int), intent(in)     :: NU(3)

    integer(c_int), intent(in)     :: SV(3)
    integer(c_int), intent(in)     :: NV(3)

    integer(c_int), intent(in)     :: SW(3)
    integer(c_int), intent(in)     :: NW(3)

    integer(c_int), intent(in)     :: BC_3L_global

    real(c_double), intent(in)     :: x1u( bL(1):(N(1)+bU(1)) )
    real(c_double), intent(in)     :: x1p( bL(1):(N(1)+bU(1)) )
    real(c_double), intent(in)     :: x2p( bL(2):(N(2)+bU(2)) )

    real(c_double), intent(in)     :: x3p( bL(3):(N(3)+bU(3)) )
    real(c_double), intent(in)     :: x3w( bL(3):(N(3)+bU(3)) )

    integer(c_int), intent(in)     :: dist_type

    real(c_double), intent(in)     :: vortex_ampli_prim
    real(c_double), intent(in)     :: vortex_x1pos, vortex_x3pos
    real(c_double), intent(in)     :: vortex_radius, vortex_band

    real(c_double),  intent(inout) :: velU(bL(1):(N(1)+bU(1)),bL(2):(N(2)+bU(2)),bL(3):(N(3)+bU(3)))
    real(c_double),  intent(inout) :: velV(bL(1):(N(1)+bU(1)),bL(2):(N(2)+bU(2)),bL(3):(N(3)+bU(3)))
    real(c_double),  intent(inout) :: velW(bL(1):(N(1)+bU(1)),bL(2):(N(2)+bU(2)),bL(3):(N(3)+bU(3)))

    integer(c_int)                ::  i, ii, j, k

    !! mjohn 051012 - variables required for changing the nondimensionalization
    integer(c_int)                ::  No, Ne, aborted, merror
    real(c_double)                ::  pi
    !! mjohn 051012

    velU( SU(1):NU(1),SU(2):NU(2),SU(3):NU(3) ) = 0.
    velV( SV(1):NV(1),SV(2):NV(2),SV(3):NV(3) ) = 0.
    velW( SW(1):NW(1),SW(2):NW(2),SW(3):NW(3) ) = 0.

    aborted = 0

    pi = 4.*atan(1.)    !!set constants

    ! wall-normal v component
    do k = SU(3), NU(3)
      do j = SU(2), NU(2)
        do i = SU(1), NU(1)
          select case (dist_type)
          case (0)                 !----- mjohn 301111 - generic inflow
            aborted = 1
            ! do nothing
          case (1)                 !----- mjohn 301111 - pair of primary vortices, even secondary forcing
            if( BC_3L_global == -2 ) then
              velU(i,j,k) = velU(i,j,k) + vortex_ampli_prim*eddy((x1u(i)-vortex_x1pos)/vortex_radius,(x3p(k)-vortex_x3pos)/vortex_radius,1)
            ELSE
              velU(i,j,k) = velU(i,j,k) + vortex_ampli_prim*eddy((x1u(i)-vortex_x1pos)/vortex_radius,(x3p(k)-vortex_x3pos)/vortex_radius,1)
              velU(i,j,k) = velU(i,j,k) - vortex_ampli_prim*eddy((x1u(i)-vortex_x1pos)/vortex_radius,(x3p(k)+vortex_x3pos)/vortex_radius,1)
            END IF
          CASE (2)                 !----- mjohn 301111 - odd secondary forcing
            IF (BC_3L_global == -2) THEN
              WRITE(*,*) "ERROR! Wrong disturbance selected. Symmetric setup does not allow for odd secondary forcing. Simulation aborted."
              CALL MPI_FINALIZE(merror)
            ELSE
              velU(i,j,k) = velU(i,j,k) + vortex_ampli_prim*eddy((x1u(i)-vortex_x1pos)/vortex_radius,(x3p(k)-vortex_x3pos)/vortex_radius,1)
              velU(i,j,k) = velU(i,j,k) - vortex_ampli_prim*eddy((x1u(i)-vortex_x1pos)/vortex_radius,(x3p(k)+vortex_x3pos)/vortex_radius,1)
            END IF
          CASE (3)                 !----- mjohn 301111 - mixed secondary forcing
            IF (BC_3L_global == -2) THEN
              WRITE(*,*) "ERROR! Wrong disturbance selected. Symmetric setup does not allow for mixed secondary forcing. Simulation aborted."
              CALL MPI_FINALIZE(merror)
            ELSE
              velU(i,j,k) = velU(i,j,k) + vortex_ampli_prim*eddy((x1u(i)-vortex_x1pos)/vortex_radius,(x3p(k)-vortex_x3pos)/vortex_radius,1)
              velU(i,j,k) = velU(i,j,k) + vortex_ampli_prim*eddy((x1u(i)-vortex_x1pos)/vortex_radius,(x3p(k)+vortex_x3pos)/vortex_radius,1)
              velU(i,j,k) = velU(i,j,k) - vortex_ampli_prim*eddy((x1u(i)-vortex_x1pos)/vortex_radius,(x3p(k)             )/vortex_radius,1)
            END IF
          CASE (4)                 !----- mjohn 260912 - array of primary vortices, no secondary forcing (without distinction between symmetric and non-symmetric case / without IF(BC_3L_..) statement)
            ! get number of even and odd vortices to be initialized within domain
            Ne =   floor((vortex_band/vortex_x3pos +1)/4)
            No = ceiling((vortex_band/vortex_x3pos -1)/4)
            Ne = min(Ne,No)
            No = min(Ne,No)
            ! initialize No odd vortices at positions (1*x3pos, 5*x3pos, ...) and Ne even vortices (with opposite sign) at positions (3*x3pos, 7*x3pos, ...)
            do ii = -ceiling((No-1)/2.0), floor((No-1)/2.0)
              velU(i,j,k) = velU(i,j,k) + vortex_ampli_prim*eddy((x1u(i)-vortex_x1pos)/vortex_radius,(x3p(k) - vortex_x3pos - 4*ii*vortex_x3pos)/vortex_radius,1)
              velU(i,j,k) = velU(i,j,k) - vortex_ampli_prim*eddy((x1u(i)-vortex_x1pos)/vortex_radius,(x3p(k) + vortex_x3pos + 4*ii*vortex_x3pos)/vortex_radius,1)
            end do
          CASE (5)
            aborted = 1
          CASE (6)
            aborted = 1
          CASE (7)
            aborted = 1
          CASE (8)                 !----- mjohn 300513 - pair of primary vortices with modulation in 2 (z) direction
            velU(i,j,k) = velU(i,j,k) + vortex_ampli_prim*eddy((x1u(i)-vortex_x1pos)/vortex_radius,(x3p(k)-vortex_x3pos)/vortex_radius,1)*(0.75+0.25*sin(2*pi*x2p(j)/(4*vortex_radius)))
            velU(i,j,k) = velU(i,j,k) - vortex_ampli_prim*eddy((x1u(i)-vortex_x1pos)/vortex_radius,(x3p(k)+vortex_x3pos)/vortex_radius,1)*(0.75+0.25*sin(2*pi*x2p(j)/(4*vortex_radius)))
          CASE (9)                 !----- mjohn 310513 - array of primary vortices with modulation in 2 (z) direction
            ! get number of even and odd vortices to be initialized within domain
            Ne =   floor((vortex_band/vortex_x3pos +1)/4)
            No = ceiling((vortex_band/vortex_x3pos -1)/4)
            Ne = min(Ne,No)
            No = min(Ne,No)
            ! initialize No odd vortices at positions (1*x3pos, 5*x3pos, ...) and Ne even vortices (with opposite sign) at positions (3*x3pos, 7*x3pos, ...)
            do ii = -ceiling((No-1)/2.0), floor((No-1)/2.0)
              velU(i,j,k) = velU(i,j,k) + vortex_ampli_prim*eddy((x1u(i)-vortex_x1pos)/vortex_radius,(x3p(k) - vortex_x3pos - 4*ii*vortex_x3pos)/vortex_radius,1)*(0.75+0.25*sin(2*pi*x2p(j)/(4*vortex_radius)))
              velU(i,j,k) = velU(i,j,k) - vortex_ampli_prim*eddy((x1u(i)-vortex_x1pos)/vortex_radius,(x3p(k) + vortex_x3pos + 4*ii*vortex_x3pos)/vortex_radius,1)*(0.75+0.25*sin(2*pi*x2p(j)/(4*vortex_radius)))
            end do

          case default
            aborted = 1
            ! do nothing
          end select
        end do
      end do
    end do

    ! chordwise u component
    do k = SW(3), NW(3)
      do j = SW(2), NW(2)
        do i = SW(1), NW(1)
          select case (dist_type)
          case (0)                 !----- mjohn 301111 - generic inflow conditions of gh-type
            aborted = 1
            ! do nothing
          case (1)                 !----- mjohn 301111 - even secondary forcing
            if (BC_3L_global == -2) then
              velW(i,j,k) = velW(i,j,k) + vortex_ampli_prim*eddy((x1p(i)-vortex_x1pos)/vortex_radius,(x3w(k)-vortex_x3pos)/vortex_radius,3)
            else
              velW(i,j,k) = velW(i,j,k) + vortex_ampli_prim*eddy((x1p(i)-vortex_x1pos)/vortex_radius,(x3w(k)-vortex_x3pos)/vortex_radius,3)
              velW(i,j,k) = velW(i,j,k) - vortex_ampli_prim*eddy((x1p(i)-vortex_x1pos)/vortex_radius,(x3w(k)+vortex_x3pos)/vortex_radius,3)
            end if
          case (2)                 !----- mjohn 301111 - odd secondary forcing
            IF (BC_3L_global == -2) THEN
              WRITE(*,*) "ERROR! Wrong disturbance selected. Symmetric setup does not allow for odd secondary forcing. Simulation aborted."
              CALL MPI_FINALIZE(merror)
            else
              velW(i,j,k) = velW(i,j,k) + vortex_ampli_prim*eddy((x1p(i)-vortex_x1pos)/vortex_radius,(x3w(k)-vortex_x3pos)/vortex_radius,3)
              velW(i,j,k) = velW(i,j,k) - vortex_ampli_prim*eddy((x1p(i)-vortex_x1pos)/vortex_radius,(x3w(k)+vortex_x3pos)/vortex_radius,3)
            end if
          case (3)                 !----- mjohn 301111 - mixed secondary forcing
            IF (BC_3L_global == -2) THEN
              WRITE(*,*) "ERROR! Wrong disturbance selected. Symmetric setup does not allow for mixed secondary forcing. Simulation aborted."
              CALL MPI_FINALIZE(merror)
            else
              velW(i,j,k) = velW(i,j,k) + vortex_ampli_prim*eddy((x1p(i)-vortex_x1pos)/vortex_radius,(x3w(k)-vortex_x3pos)/vortex_radius,3)
              velW(i,j,k) = velW(i,j,k) + vortex_ampli_prim*eddy((x1p(i)-vortex_x1pos)/vortex_radius,(x3w(k)+vortex_x3pos)/vortex_radius,3)
              velW(i,j,k) = velW(i,j,k) - vortex_ampli_prim*eddy((x1p(i)-vortex_x1pos)/vortex_radius,(x3w(k)             )/vortex_radius,3)
            END IF
          case (4)                 !----- mjohn 260912 - array of primary vortices, no secondary forcing (without distinction between symmetric and non-symmetric case / without IF(BC_3L_..) statement)
            ! get number of even and odd vortices to be initialized within domain
            Ne =   floor((vortex_band/vortex_x3pos +1)/4)
            No = ceiling((vortex_band/vortex_x3pos -1)/4)
            Ne = min(Ne,No)
            No = min(Ne,No)
            ! initialize No odd vortices at positions (1*x3pos, 5*x3pos, ...) and Ne even (with opposite sign) vortices at positions (3*x3pos, 7*x3pos, ...)
            do ii = -ceiling((No-1)/2.0), floor((No-1)/2.0)
              velW(i,j,k) = velW(i,j,k) + vortex_ampli_prim*eddy((x1p(i)-vortex_x1pos)/vortex_radius,(x3w(k) - vortex_x3pos - 4*ii*vortex_x3pos)/vortex_radius,3)
              velW(i,j,k) = velW(i,j,k) - vortex_ampli_prim*eddy((x1p(i)-vortex_x1pos)/vortex_radius,(x3w(k) + vortex_x3pos + 4*ii*vortex_x3pos)/vortex_radius,3)
            end do
          case (5)
            aborted = 1
          case (6)
            aborted = 1
          case (7)
            aborted = 1
          case (8)                 !----- mjohn 300513 - pair of primary vortices with modulation in 2 (z) direction
            velU(i,j,k) = velU(i,j,k) + vortex_ampli_prim*eddy((x1u(i)-vortex_x1pos)/vortex_radius,(x3p(k)-vortex_x3pos)/vortex_radius,1)*(0.75+0.25*sin(2*pi*x2p(j)/(4*vortex_radius)))
            velU(i,j,k) = velU(i,j,k) - vortex_ampli_prim*eddy((x1u(i)-vortex_x1pos)/vortex_radius,(x3p(k)+vortex_x3pos)/vortex_radius,1)*(0.75+0.25*sin(2*pi*x2p(j)/(4*vortex_radius)))
          case (9)                 !----- mjohn 310513 - array of primary vortices with modulation in 2 (z) direction
            ! get number of even and odd vortices to be initialized within domain
            Ne =   floor((vortex_band/vortex_x3pos +1)/4)
            No = ceiling((vortex_band/vortex_x3pos -1)/4)
            Ne = min(Ne,No)
            No = min(Ne,No)
            ! initialize No odd vortices at positions (1*x3pos, 5*x3pos, ...) and Ne even (with opposite sign) vortices at positions (3*x3pos, 7*x3pos, ...)
            do ii = -ceiling((No-1)/2.0), floor((No-1)/2.0)
              velW(i,j,k) = velW(i,j,k) + vortex_ampli_prim*eddy((x1p(i)-vortex_x1pos)/vortex_radius,(x3w(k) - vortex_x3pos - 4*ii*vortex_x3pos)/vortex_radius,3)*(0.75+0.25*sin(2*pi*x2p(j)/(4*vortex_radius)))
              velW(i,j,k) = velW(i,j,k) - vortex_ampli_prim*eddy((x1p(i)-vortex_x1pos)/vortex_radius,(x3w(k) + vortex_x3pos + 4*ii*vortex_x3pos)/vortex_radius,3)*(0.75+0.25*sin(2*pi*x2p(j)/(4*vortex_radius)))
            end do
          case default
            aborted = 1
            ! do nothing
          end select
        end do
      end do
    end do

    !----- mjohn 260912
    if( rank == 0 .and. aborted == 1) then
      WRITE(*,*) "WARNING! No valid initial disturbance defined. Initial conditions equal to baseflow only."
    end if

  end subroutine VF_init_Dist

end module cmod_VectorField

max_stars_repo_name: huppd/PINTimpact

max_issues_repo_name: huppd/PINTimpact

max_forks_repo_name: huppd/PINTimpact


module global_module

  implicit none

  ! Time integration

  real(kind=8), parameter :: t_end      = 0.2d0 ! [s] total simulated time (from 0 to t_end)
  real(kind=8), parameter :: CFL_target = 0.2d0

  real(kind=8), parameter :: x_min = -0.5d0 ! [m]
  real(kind=8), parameter :: x_max = 0.5d0  ! [m]

  ! Total number of cells (including 4 ghost cells for each side)
  ! (using 4 ghost cells makes it easy to reach higher order accuracy)

  integer, parameter :: Nx = 200

  logical, parameter :: bool_PERIODIC = .false.

  ! Note: we use 4 ghost cells per side. This allows to reach up to 5th order
  ! (with the WENO method)
  real(kind=8), parameter :: dx = (x_max-x_min)/(Nx-8)

  ! Order in space
  ! 1: uniform sol in cell, first order in space
  ! 2: MUSCL with TVD slope limiter, second order in space
  ! 3: WENO 3rd order
  ! 5: WENO 5th order
  integer, parameter :: space_order = 3

  ! Order in time
  ! 1: forward Euler, first order in time
  ! 2: midpoint Euler, second order in time
  ! 3: explicit TVD Runge-Kutta, third order in time
  integer, parameter :: time_order = 3

  ! Space fluxes: choose ONLY one
  logical, parameter :: flux_type_Rusanov = .false.
  logical, parameter :: flux_type_HLL     = .true.

  ! Utilities
  real(kind=8) :: ws_maxabs ! Global variable, maximum wave speed

end module


text:
      SUBROUTINE MFFTA7(C,FAC)
*
*   PURPOSE:
*       ELEMENTARY COOLEY-TUKEY RADIX 2 STEP APPLIED TO A VECTOR-OF
*       2-VECTORS-OF-COMPLEX C[IMS,NM [IVS,NV [IES,NE]]].
*       SEE REF.[1] FOR NOTATIONS.
*       THIS ROUTINE CAN BE USED ONLY BY ROUTINE MFFTIM, WHICH CONTROLS
*       ITS OPERATION THROUGH THE MFFTPA COMMON
*
*   DUMMY ARGUMENTS :
*
*   C   ARRAY BEING FOURIER  TRANSFORMED
*   FAC PHASE FACTORS, PREPARED BY MFFTP; NOT MODIFIED IN OUTPUT
*
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON /MFFTPA/  IMS,IVS,IES,NM,NV,NE,MX,LX,MLIM,MSTEP,LLIM,LSTEP,
     $ NUSTEP,IVLIM,ILIM,MD2LIM,LD2LIM
      INTEGER NUSTEP
      DOUBLE COMPLEX C(0:NUSTEP-1,0:1)
      COMPLEX FAC(0:*)
      DOUBLE COMPLEX T0,F

        IF(2*MX.GT.LX)THEN
          IF(LX.EQ.1)GOTO 1000
*   ELSE HERE
          LAMF=MX
          DO 100 LAM=LSTEP,LD2LIM,LSTEP
            F=FAC(LAMF)
            DO 90 MU=LAM,LAM+MLIM,MSTEP
              DO 80 IV=MU,MU+IVLIM,IVS
              DO 80 I=IV,IV+ILIM,IES
                T0=C(I,1)*F
                C(I,1)=C(I,0)-T0
                C(I,0)=C(I,0)+T0
80            CONTINUE
90          CONTINUE
            LAMF=LAMF+MX
100       CONTINUE
          LAMF=LAMF-MX
          DO 101 LAM=LD2LIM+2*LSTEP,LLIM,LSTEP
            F=-CONJG(FAC(LAMF))
            DO 91 MU=LAM,LAM+MLIM,MSTEP
              DO 81 IV=MU,MU+IVLIM,IVS
              DO 81 I=IV,IV+ILIM,IES
                T0=C(I,1)*F
                C(I,1)=C(I,0)-T0
                C(I,0)=C(I,0)+T0
81            CONTINUE
91          CONTINUE
            LAMF=LAMF-MX
101       CONTINUE
          DO 93 MU=LD2LIM+LSTEP,LD2LIM+LSTEP+MLIM,MSTEP

            DO 83 IV=MU,MU+IVLIM,IVS
            DO 83 I=IV,IV+ILIM,IES
              T0=DCMPLX(-DIMAG(C(I,1)),DREAL(C(I,1)))
              C(I,1)=C(I,0)-T0
              C(I,0)=C(I,0)+T0
83          CONTINUE
93        CONTINUE
1000      DO 92 MU=0,MLIM,MSTEP

            DO 82 IV=MU,MU+IVLIM,IVS
            DO 82 I=IV,IV+ILIM,IES
              T0=C(I,1)
              C(I,1)=C(I,0)-T0
              C(I,0)=C(I,0)+T0
82          CONTINUE
92        CONTINUE
        ELSE
          DO 200 MU=0,MLIM,MSTEP
            LAMF=MX
            DO 190 LAM=MU+LSTEP,MU+LD2LIM,LSTEP
              F=FAC(LAMF)
              DO 180 IV=LAM,LAM+IVLIM,IVS
              DO 180 I=IV,IV+ILIM,IES
                T0=C(I,1)*F
                C(I,1)=C(I,0)-T0
                C(I,0)=C(I,0)+T0
180           CONTINUE
              LAMF=LAMF+MX
190         CONTINUE
            LAMF=LAMF-MX
            DO 191 LAM=MU+LD2LIM+2*LSTEP,MU+LLIM,LSTEP
              F=-CONJG(FAC(LAMF))
              DO 181 IV=LAM,LAM+IVLIM,IVS
              DO 181 I=IV,IV+ILIM,IES
                T0=C(I,1)*F
                C(I,1)=C(I,0)-T0
                C(I,0)=C(I,0)+T0
181           CONTINUE
              LAMF=LAMF-MX
191         CONTINUE
            DO 182 IV=MU,MU+IVLIM,IVS
            DO 182 I=IV,IV+ILIM,IES
              T0=C(I,1)
              C(I,1)=C(I,0)-T0
              C(I,0)=C(I,0)+T0
182         CONTINUE
            DO 183 IV=MU+LD2LIM+LSTEP,MU+LD2LIM+LSTEP+IVLIM,IVS
            DO 183 I=IV,IV+ILIM,IES
              T0=DCMPLX(-DIMAG(C(I,1)),DREAL(C(I,1)))
              C(I,1)=C(I,0)-T0
              C(I,0)=C(I,0)+T0
183         CONTINUE
200       CONTINUE
        ENDIF
      END

max_stars_repo_name: cea-trust-platform/trust-code

max_issues_repo_name: pledac/trust-code

max_forks_repo_name: pledac/trust-code


text:
      SUBROUTINE WF_R_MTRX
     &          (index)

c     Computes reflection coefficient matrix from p matrix and
c     returns it in r.

      implicit real*8 (a-h,o-z)

c     Wave Fields
      complex  * 16 c,s,p,dpdh,
     &              r11,r22,r12,r21,rbar11,rbar22

      common/wf_pmtx/
     &              c(2),s(2),p(4,2,2),dpdh(4,2,2)
     &      /wf_rmtx/
     &              r11,r22,r12,r21,rbar11,rbar22

c     Local
      complex  * 16 g12,g13,g14,g23,g24,g34,
     &              d00,d11,d22,d12,d21

      g12=p(1,1,index)*p(2,2,index)-p(1,2,index)*p(2,1,index)
      g13=p(1,1,index)*p(3,2,index)-p(1,2,index)*p(3,1,index)
      g14=p(1,1,index)*p(4,2,index)-p(1,2,index)*p(4,1,index)
      g23=p(2,1,index)*p(3,2,index)-p(2,2,index)*p(3,1,index)
      g24=p(2,1,index)*p(4,2,index)-p(2,2,index)*p(4,1,index)
      g34=p(3,1,index)*p(4,2,index)-p(3,2,index)*p(4,1,index)

      d00=-g13+c(index)*( g34-g12+c(index)*g24)
      d11= g13+c(index)*( g34+g12+c(index)*g24)
      d22= g13+c(index)*(-g34-g12+c(index)*g24)
      d12=2.d0*c(index)*g14
      d21=2.d0*c(index)*g23

      r11=d11/d00
      r22=d22/d00
      r12=d12/d00
      r21=d21/d00
      RETURN
      END      ! WF_R_MTRX

ext: for

max_stars_repo_name: spinorkit/LWPC

max_issues_repo_name: spinorkit/LWPC

max_forks_repo_name: spinorkit/LWPC


text:
! Additions by Luke Madaus, Nov 2015, to generate simple
! Surface observations for CM1 OSSE experiments

! BEGIN DART PREPROCESS KIND LIST
!TEMPERATURE_2M,             QTY_2M_TEMPERATURE,       COMMON_CODE
!U_WIND_10M,                 QTY_10M_U_WIND_COMPONENT, COMMON_CODE
!V_WIND_10M,                 QTY_10M_V_WIND_COMPONENT, COMMON_CODE
!SURFACE_PRESSURE,           QTY_SURFACE_PRESSURE,     COMMON_CODE
!SPECIFIC_HUMIDITY_2M,       QTY_SPECIFIC_HUMIDITY,    COMMON_CODE
! END DART PREPROCESS KIND LIST

! BEGIN DART PREPROCESS USE OF SPECIAL OBS_DEF MODULE
! END DART PREPROCESS USE OF SPECIAL OBS_DEF MODULE

! BEGIN DART PREPROCESS GET_EXPECTED_OBS_FROM_DEF
! END DART PREPROCESS GET_EXPECTED_OBS_FROM_DEF

! BEGIN DART PREPROCESS READ_OBS_DEF
! END DART PREPROCESS READ_OBS_DEF

! BEGIN DART PREPROCESS WRITE_OBS_DEF
! END DART PREPROCESS WRITE_OBS_DEF

! BEGIN DART PREPROCESS INTERACTIVE_OBS_DEF
! END DART PREPROCESS INTERACTIVE_OBS_DEF

! BEGIN DART PREPROCESS MODULE CODE
! END DART PREPROCESS MODULE CODE

! !!! Note about Specific Humidity observations:
! !!! UNITS in original BUFR are g/kg; This is converted to kg/kg by
! !!! the BUFR to obs_sequence conversion programs making it unnecessary
! !!! to multiply by 1000 at assimilation time.
! !!! PLEASE pay attention to units for specific humidity in models.

! <next few lines under version control, do not edit>
! $URL$
! $Id$
! $Revision$
! $Date$

max_stars_repo_name: hkershaw-brown/feature-preprocess

max_issues_repo_name: hkershaw-brown/feature-preprocess

max_forks_repo_name: hkershaw-brown/feature-preprocess


text:
C
C $Id: idsfft.f,v 1.6 2008-07-27 00:17:30 haley Exp $
C
C                Copyright (C)  2000
C        University Corporation for Atmospheric Research
C                All Rights Reserved
C
C The use of this Software is governed by a License Agreement.
C
      SUBROUTINE IDSFFT (MD,NDP,XD,YD,ZD,NXI,NYI,NZI,XI,YI,ZI,IWK,WK)
C
C DIMENSION OF           XD(NDP), YD(NDP), ZD(NDP), XI(NXI),
C ARGUMENTS              YI(NYI), ZI(NZI,NYI), WK(6*NDP),
C                        IWK(31*NDP + NXI*NYI)
C
C PURPOSE                This subroutine performs smooth surface
C                        fitting when the projections of the data
C                        points in the X-Y plane are irregularly
C                        distributed in the plane.
C
C USAGE                  CALL IDSFFT (MD,NDP,XD,YD,ZD,NXI,NYI,NZI,
C                                     XI,YI,ZI,IWK,WK)
C
C ARGUMENTS
C
C ON INPUT               MD
C                          Mode of computation (must be 1, 2, or 3,
C                          else an error return will occur).
C                           = 1 if this is the first call to this
C                               subroutine, or if the value of NDP
C                               has been changed from the previous
C                               call, or if the contents of the XD
C                               or YD arrays have been changed from
C                               the previous call.
C                           = 2 if the values of NDP and the XD and
C                               YD arrays are unchanged from the
C                               previous call, but new values for
C                               XI and YI are being used.  If MD = 2
C                               and NDP has been changed since the
C                               previous call to IDSFFT, an error
C                               return occurs.
C                           = 3 if the values of NDP, NXI, NYI, XD,
C                               YD, XI, and YI are unchanged from the
C                               previous call, i.e., if the only change
C                               on input to IDSFFT is in the ZD array.
C                               If MD = 3 and NDP, NXI, or NYI has been
C                               changed since the previous call to
C                               IDSFFT, an error return occurs.
C
C                           Between the call with MD = 2 or MD = 3 and
C                           the preceding call, the IWK and WK work
C                           arrays should not be disturbed.
C
C                        NDP
C                          Number of data points (must be 4 or
C                          greater, else an error return will occur).
C
C                        XD
C                          Array of dimension NDP containing the X
C                          coordinates of the data points.
C
C                        YD
C                          Array of dimension NDP containing the Y
C                          coordinates of the data points.
C
C                        ZD
C                          Array of dimension NDP containing the Z
C                          coordinates of the data points.
C
C                        NXI
C                          Number of output grid points in the X
C                          direction (must be 1 or greater, else
C                          an error return will occur).
C
C                        NYI
C                          Number of output grid points in the Y
C                          direction (must be 1 or greater, else
C                          an error return will occur).
C
C                        NZI
C                          First dimension of ZI as declared in the
C                          calling program.  NZI must be greater than
C                          or equal to NXI, else an error return will
C                          occur.
C
C                        XI
C                          Array of dimension NXI containing the
C                          X coordinates of the output grid points.
C
C                        YI
C                         Array of dimension NYI containing the
C                         Y coordinates of the output grid points.
C
C                        IWK
C                          Integer work array of dimension at
C                          least 31*NDP + NXI*NYI.
C
C                        WK
C                          Real work array of dimension at least 6*NDP.
C
C ON OUTPUT              ZI
C                           Real, two-dimensional array of dimension
C                           (NZI,NYI), storing the interpolated Z
C                           values at the output grid points.
C
C SPECIAL CONDITIONS     Inadequate work space IWK and WK may cause
C                        incorrect results.
C
C                        The data points must be distinct and their
C                        projections in the X-Y plane must not be
C                        collinear, else an error return occurs.
C
C IDSFFT calls the subroutines IDGRID, IDPDRV, IDPTIP, and IDTANG.
C
C Declaration statements.
C
      DIMENSION XD(NDP),YD(NDP),ZD(NDP),XI(NXI),YI(NYI),ZI(NZI,NYI),
     +          IWK(31*NDP+NXI*NYI),WK(6*NDP)
C
      COMMON /IDPT/ ITPV,DMMY(27)
      SAVE   /IDPT/
C
      COMMON /IDCOMN/ INTY,ITTY,ALSP,BLSP,CLSP,XAVG,YAVG
      SAVE   /IDCOMN/
C
C Check for an uncleared prior error.
C
        IF (ICFELL('IDSFFT (BIVAR) - UNCLEARED PRIOR ERROR',1).NE.0)
     +                                                        RETURN
C
C Check for input errors.
C
      IF (MD.LT.1.OR.MD.GT.3) THEN
        CALL SETER ('IDSFFT (BIVAR) - INPUT VARIABLE MD IS OUT OF RANGE'
     +,2,1)
        RETURN
      END IF
C
      IF (NDP.LT.4) THEN
        CALL SETER ('IDSFFT (BIVAR) - INPUT VARIABLE NDP IS OUT OF RANGE
     +',3,1)
        RETURN
      END IF
C
      IF (NXI.LT.1.OR.NYI.LT.1) THEN
        CALL SETER ('IDSFFT (BIVAR) - INPUT VARIABLE NXI OR NYI IS OUT O
     +F RANGE',4,1)
        RETURN
      END IF
C
      IF (NXI.GT.NZI) THEN
        CALL SETER ('IDSFFT (BIVAR) - INPUT VARIABLE NZI IS LESS THAN NX
     +I',5,1)
        RETURN
      END IF
C
      IF (MD.LE.1) THEN
C
        IWK(1)=NDP
C
      ELSE
C
        NDPPV=IWK(1)
C
        IF (NDP.NE.NDPPV) THEN
          CALL SETER  ('IDSFFT (BIVAR) - MD = 2 OR 3 BUT NDP WAS CHANGED
     + SINCE LAST CALL',6,1)
          RETURN
        END IF
C
      END IF
C
      IF (MD.LE.2) THEN
C
        IWK(2)=INTY
        IWK(3)=NXI
        IWK(4)=NYI
C
      ELSE
C
        INTYPV=IWK(2)
C
        IF (INTY.NE.INTYPV) THEN
          CALL SETER ('IDSFFT (BIVAR) - MD = 3 BUT ITY WAS CHANGED SINCE
     + LAST CALL',7,1)
           RETURN
        END IF
C
        NXIPV=IWK(3)
C
        IF (NXI.NE.NXIPV) THEN
          CALL SETER ('IDSFFT (BIVAR) - MD = 3 BUT NXI WAS CHANGED SINCE
     + LAST CALL',8,1)
           RETURN
        END IF
C
        NYIPV=IWK(4)
C
        IF (NYI.NE.NYIPV) THEN
          CALL SETER ('IDSFFT (BIVAR) - MD = 3 BUT NYI WAS CHANGED SINCE
     + LAST CALL',9,1)
          RETURN
        END IF
C
      END IF
C
C Allocate storage areas in the array IWK.
C
      JWIPT=16
      JWIWL=6*NDP+1
      JWNGP0=JWIWL-1
      JWIPL=24*NDP+1
      JWIWP=30*NDP+1
      JWIGP0=31*NDP
      JWWPD=5*NDP+1
C
C If MD = 1, triangulate the X/Y plane.  (If MD = 2 or 3, this has
C already been done and need not been redone.)
C
      IF (MD.EQ.1) THEN
        CALL IDTANG (NDP,XD,YD,NT,IWK(JWIPT),NL,IWK(JWIPL),IWK(JWIWL),
     +                                                  IWK(JWIWP),WK)
        IF (ICFELL('IDSFFT',10).NE.0) RETURN
        IWK(5)=NT
        IWK(6)=NL
      ELSE
        NT=IWK(5)
        NL=IWK(6)
      END IF
C
      IF (NT.EQ.0) RETURN
C
C If linear interpolation is activated, compute the coefficients of the
C least-squares-fit plane and the mean values of X and Y.
C
      IF (INTY.NE.0) THEN
        CALL IDLSQF (XD,YD,ZD,NDP,ALSP,BLSP,CLSP,XAVG,YAVG)
      END IF
C
C If MD = 1 or 2, sort output grid points in ascending order by the
C number of the containing region (a triangle number or a combined
C pair of border line segment numbers).  (If MD = 3, this has already
C been done and need not be redone.)
C
      IF (MD.LE.2) THEN
        CALL IDGRID (XD,YD,NT,IWK(JWIPT),NL,IWK(JWIPL),NXI,
     +               NYI,XI,YI,IWK(JWNGP0+1),IWK(JWIGP0+1))
        IF (ICFELL('IDSFFT',11).NE.0) RETURN
      END IF
C
C If quintic interpolation is activated, estimate partial derivatives.
C
      IF (INTY.EQ.0) THEN
        CALL IDPDRV (NDP,XD,YD,ZD,NT,IWK(JWIPT),WK,WK(JWWPD))
      END IF
C
C Interpolate to get ZI values.
C
      ITPV=0
      JIG0MX=0
      JIG1MN=NXI*NYI+1
      NNGP=NT+2*NL
C
      DO 105 JNGP=1,NNGP
C
        ITI=JNGP
        IF (JNGP.LE.NT) GO TO 101
        IL1=(JNGP-NT+1)/2
        IL2=(JNGP-NT+2)/2
        IF (IL2.GT.NL) IL2=1
        ITI=IL1*(NT+NL)+IL2
C
  101   JWNGP=JWNGP0+JNGP
        NGP0=IWK(JWNGP)
        IF (NGP0.EQ.0) GO TO 103
        JIG0MN=JIG0MX+1
        JIG0MX=JIG0MX+NGP0
C
        DO 102 JIGP=JIG0MN,JIG0MX
          JWIGP=JWIGP0+JIGP
          IZI=IWK(JWIGP)
          IYI=(IZI-1)/NXI+1
          IXI=IZI-NXI*(IYI-1)
          CALL IDPTIP(XD,YD,ZD,NT,IWK(JWIPT),NL,IWK(JWIPL),WK,
     +                        ITI,XI(IXI),YI(IYI),ZI(IXI,IYI))
  102   CONTINUE
C
  103   JWNGP=JWNGP0+2*NNGP+1-JNGP
        NGP1=IWK(JWNGP)
        IF (NGP1.EQ.0) GO TO 105
        JIG1MX=JIG1MN-1
        JIG1MN=JIG1MN-NGP1
C
        DO 104 JIGP=JIG1MN,JIG1MX
          JWIGP=JWIGP0+JIGP
          IZI=IWK(JWIGP)
          IYI=(IZI-1)/NXI+1
          IXI=IZI-NXI*(IYI-1)
          CALL IDPTIP(XD,YD,ZD,NT,IWK(JWIPT),NL,IWK(JWIPL),WK,
     +                        ITI,XI(IXI),YI(IYI),ZI(IXI,IYI))
  104   CONTINUE
C
  105 CONTINUE
C
C Done.
C
        RETURN
C
      END

max_stars_repo_name: tenomoto/ncl

max_issues_repo_name: tenomoto/ncl

max_forks_repo_name: tenomoto/ncl


text:
program basic_stats
  implicit none

  integer, parameter :: i64 = selected_int_kind(18)
  integer, parameter :: r64 = selected_real_kind(15)
  integer(i64), parameter :: samples = 1000000000_i64

  real(r64) :: r
  real(r64) :: mean, stddev
  real(r64) :: sumn = 0, sumnsq = 0
  integer(i64) :: n = 0
  integer(i64) :: bin(10) = 0
  integer :: i, ind

  call random_seed

  n = 0
  do while(n <= samples)
    call random_number(r)
    ind = r * 10 + 1
    bin(ind) = bin(ind) + 1_i64
    sumn = sumn + r
    sumnsq = sumnsq + r*r
    n = n + 1_i64
  end do

  mean = sumn / n
  stddev = sqrt(sumnsq/n - mean*mean)
  write(*, "(a, i0)") "sample size = ", samples
  write(*, "(a, f17.15)") "Mean :   ", mean,
  write(*, "(a, f17.15)") "Stddev : ", stddev
  do i = 1, 10
    write(*, "(f3.1, a, a)") real(i)/10.0, ": ", repeat("=", int(bin(i)*500/samples))
  end do

end program

max_issues_repo_name: seanwallawalla-forks/RosettaCodeData

max_forks_repo_name: seanwallawalla-forks/RosettaCodeData