text:
! $Id$
! -----------------------------------------------------------
! Module spline_module
! -----------------------------------------------------------
! Code area 9: general
! -----------------------------------------------------------

!!****h* Conquest/spline_module
!!  NAME
!!
!!  PURPOSE
!!
!!  AUTHOR
!!   D.R.Bowler
!!  CREATION DATE
!!   23/05/01
!!  MODIFICATION HISTORY
!!   31/05/2001 dave
!!    Added RCS Id tag to header, increased ROBODoc documentation
!!   20/06/2001 dave
!!    Added RCS Log tag and cq_abort
!!   15:48, 08/04/2003 drb
!!    Bug fixes in splint and dsplint
!!***
module spline_module

  use global_module, only: area_general

  implicit none

  ! RCS tag for object file identification
  character(len=80), save, private :: RCSid = "$Id$"

contains

! -----------------------------------------------------------
! Subroutine spline
! -----------------------------------------------------------

!!****f* spline_module/spline *
!!
!!  NAME
!!   spline
!!  USAGE
!!
!!  PURPOSE
!!   Given an array y(1:n) containing a function tabulated at
!!   REGULAR intervals with spacing delta_x, and given values
!!   yp1 and ypn for the first derivative of the function at
!!   the end points, this subroutine returns an array d2y(1:n)
!!   containing the second derivatives of the interpolating
!!   function at the tabulated points. If yp1 and/or ypn are
!!   equal or larger than 1.0d30, the routine is signaled to set
!!   the corresponding boundary condition for a natural spline,
!!   with zero second derivative at that point. This subroutine
!!   is based on subroutine spline(), from Numerical Recipes (2nd edition).
!!  INPUTS
!!   integer :: n - size of table to be splined
!!   real(double) :: delta_x - spacing of values in table
!!   real(double) :: yp1, ypn - end points giving boundaries
!!   real(double), dimension(n) :: y, y2 - table and second derivative
!!  USES
!!   datatypes, numbers
!!  AUTHOR
!!   E.H.Hernandez
!!  CREATION DATE
!!   29/02/96
!!  MODIFICATION HISTORY
!!   23/05/2001 dave
!!    F90, ROBODoc, indentation
!!   31/05/2001 dave
!!    More ROBODoc comments
!!  SOURCE
!!
  subroutine spline( n, delta_x, y, yp1, ypn, y2 )

    use datatypes
    use numbers
    use GenComms, only: cq_abort
    use memory_module, only: reg_alloc_mem, reg_dealloc_mem, type_dbl

    implicit none

    ! Passed variables
    integer :: n

    real(double) :: delta_x, yp1, ypn
    real(double), dimension(n) :: y, y2

    ! Local variables
    integer :: i, k, stat
    real(double) :: p, qn, sig, un
    real(double), parameter ::huge = 1.0e30_double
    real(double), dimension(:), allocatable :: u, x

    allocate(u(n), x(n), STAT=stat)
    if (stat /= 0) call cq_abort("spline: Error alloc mem: ", n)
    call reg_alloc_mem(area_general, 2*n, type_dbl)

    ! generate the array of x entries
    do i=1, n
       x(i) = delta_x * real(i-1,double)
    end do
    if ( yp1 .gt. huge ) then ! lower boundary condition set to 'natural'
       y2(1) = zero
       u(1) = zero
    else ! or else to have a specified first derivative
       y2(1) = -half
       u(1) = ( three/(x(2) - x(1))) *  &
            ((y(2) - y(1))/(x(2) - x(1)) - yp1 )
    end if

    ! this is the decomposition loop for the tridiagonal algorithm; y2 and
    ! u are used for temporary storage of the decomposed factors
    do i=2, n-1
       sig = (x(i) - x(i-1)) / (x(i+1) - x(i-1))
       p = sig * y2(i-1) + two
       y2(i) = (sig - one) / p
       u(i) = (six * ((y(i+1) - y(i))/(x(i+1) - x(i)) -  &
            (y(i) - y(i-1))/ &
            (x(i) - x(i-1)))/(x(i+1) - x(i-1)) - sig * u(i-1))/p
    end do

    if (ypn>huge) then ! the upper boundary condition is set to be 'natural'...
       qn = zero
       un = zero
    else ! or else to have a specified first derivative
       qn = half
       un = ( three/(x(n) - x(n-1))) *  &
            (ypn - (y(n) - y(n-1))/(x(n) - x(n-1)))
    end if
    y2(n) = (un - qn * u(n-1)) / (qn * y2(n-1) + one )

    ! this is the backsubstitution loop of the tridiagonal algorithm
    do k=n-1, 1, -1
       y2(k) = y2(k) * y2(k+1) + u(k)
    end do

    deallocate(u, x, STAT=stat)
    if (stat /= 0) call cq_abort("spline: Error dealloc mem")
    call reg_dealloc_mem(area_general, 2*n, type_dbl)

    return
  end subroutine spline
!!***

! -----------------------------------------------------------
! Subroutine spline_nonU
! -----------------------------------------------------------

!!****f* spline_module/spline_nonU *
!!
!!  NAME
!!   spline_nonU - makes a spline table for non-uniform x array
!!  USAGE
!!   spline_nonU(number, x array, y array,start dy/dx,final dy/dx, d2y/dx2)
!!   spline_nonU(n, y, y, yp1, ypn, y2)
!!  PURPOSE
!!   Given an array y(1:n) containing a function tabulated at
!!   IRREGULAR intervals given by array x(1:n), and given values
!!   yp1 and ypn for the first derivative of the function at
!!   the end points, this subroutine returns an array d2y(1:n)
!!   containing the second derivatives of the interpolating
!!   function at the tabulated points. If yp1 and/or ypn are
!!   equal or larger than 1.0d30, the routine is signaled to set
!!   the corresponding boundary condition for a natural spline,
!!   with zero second derivative at that point. This subroutine
!!   is based on subroutine spline(), from Numerical Recipes (2nd edition).
!!  INPUTS
!!
!!
!!  USES
!!
!!  AUTHOR
!!   D.R.Bowler
!!  CREATION DATE
!!   20/11/97
!!  MODIFICATION HISTORY
!!   23/05/2001 dave
!!    F90, ROBODoc, indentation
!!  SOURCE
!!
  subroutine spline_nonU(n, x, y, yp1, ypn, y2)

    use datatypes
    use numbers
    use GenComms, only: cq_abort
    use memory_module, only: reg_alloc_mem, reg_dealloc_mem, type_dbl

    implicit none

    ! Passed variables

    integer :: n
    real(double) :: yp1, ypn
    real(double), dimension(n) :: x, y, y2

    !     Local variables
    integer i, k, stat
    real(double) :: p, qn, sig, un
    real(double), parameter :: huge = 1.0e30_double
    real(double), dimension(:), allocatable :: u

    allocate(u(n), STAT=stat)
    if (stat /= 0) call cq_abort("spline_nonU: Error alloc mem: ", n)
    call reg_alloc_mem(area_general, n, type_dbl)

    ! set up boundary conditions
    if (yp1>huge) then ! lower boundary condition set to 'natural'
       y2(1) = zero
       u(1) = zero
    else ! or else to have a specified first derivative
       y2(1) = -half
       u(1) = ( three/(x(2) - x(1))) * ((y(2) - y(1))/(x(2) - x(1)) - yp1 )
    end if

    ! this is the decomposition loop for the tridiagonal algorithm; y2 and
    ! u are used for temporary storage of the decomposed factors
    do i=2, n-1
       sig = (x(i) - x(i-1)) / (x(i+1) - x(i-1))
       p = sig * y2(i-1) + two
       y2(i) = (sig - one) / p
       u(i) = (six * ((y(i+1) - y(i))/(x(i+1) - x(i)) -  &
            (y(i) - y(i-1))/ &
            (x(i) - x(i-1)))/(x(i+1) - x(i-1)) - sig * u(i-1))/p
    end do

    if (ypn>huge) then ! the upper boundary condition is set to be 'natural'...
       qn = zero
       un = zero
    else ! or else to have a specified first derivative
       qn = half
       un = ( three/(x(n) - x(n-1))) *  &
            (ypn - (y(n) - y(n-1))/(x(n) - x(n-1)))
    end if
    y2(n) = (un - qn * u(n-1)) / (qn * y2(n-1) + one )

    ! this is the backsubstitution loop of the tridiagonal algorithm
    do k=n-1, 1, -1
       y2(k) = y2(k) * y2(k+1) + u(k)
    end do

    deallocate(u, STAT=stat)
    if (stat /= 0) call cq_abort("spline_nonU: Error dealloc mem")
    call reg_dealloc_mem(area_general, n, type_dbl)

    return
  end subroutine spline_nonU
!!***

! -----------------------------------------------------------
! Subroutine splint_nonU
! -----------------------------------------------------------

!!****f* spline_module/splint_nonU *
!!
!!  NAME
!!   splint_nonU
!!  USAGE
!!   splint_nonU(xarray, yarray, d2y/dx2 array, size of array, x value, y, flag)
!!   splint_nonU(xa,ya,y2a,n,x,y, flag)
!!  PURPOSE
!!   NR routine for interpolation using spline table set up by spline
!!   Takes a spline table with non-uniform x spacings and, for the
!!   specified value x, finds the y value.
!!  INPUTS
!!   integer :: n (size of table)
!!   real(double) :: x,y (x value specified, y value to be found)
!!   real(double), dimension(n) :: xa, ya, y2a (x, y, d2y/dx2 arrays)
!!   logical :: flag (flag to warn if given x is out of tabulated range)
!!  USES
!!   datatypes
!!  AUTHOR
!!   D.R.Bowler
!!  CREATION DATE
!!   20/11/97
!!  MODIFICATION HISTORY
!!   23/05/2001 dave
!!    F90, ROBODoc, indentation
!!   20/06/2001 dave
!!    Added cq_abort
!!   15/04/02 lkd
!!    Changed title from splint to splint_nonU for consistency with rest of module
!!   15/04/02 lkd
!!    Added "use numbers" statement
!!   17/04/2002 lkd
!!    Added logical flag to return false if all OK, true if the given x is out of range
!!
!!  SOURCE
!!
  subroutine splint_nonU(xa,ya,y2a,n,x,y,flag)

    use datatypes
    use numbers
    use GenComms, ONLY: cq_abort

    ! Passed variables
    integer, intent(in) :: n
    real(double), intent(in) :: x
    real(double), intent(out) :: y
    real(double), dimension(n), intent(in) :: xa,y2a,ya
    logical, intent(out) :: flag

    ! Local variables
    integer :: k,khi,klo
    real(double) :: a,b,h

    ! initialize warning flag
    flag = .false.
    if (x .gt. xa(n) .or. x .lt. xa(1)) then
       flag = .true. ! the given x lies outside the tabulated values, warn the calling routine
    end if

    ! Search through the table and bracket x value
    klo = 1
    khi = n
    do while((khi-klo)>1)
       k=(khi+klo)/2
       if(xa(k).gt.x) then ! Increment high point
          khi=k
       else                ! Increment low point
          klo=k
       end if
    end do
    h=xa(khi)-xa(klo)      ! Find x spacing
    if(h==0) then
       call cq_abort('splint_nonU: bisection failed, h=0')
    end if
    ! Create value of y
    a = (xa(khi)-x)/h
    b = (x-xa(klo))/h
    y=a*ya(klo)+b*ya(khi)+ &
         ( (a*a*a-a)*y2a(klo)+(b*b*b-b)*y2a(khi) )*(h*h)/six
    return
  end subroutine splint_nonU
!!***

! -----------------------------------------------------------
! Subroutine splint
! -----------------------------------------------------------

!!****f* spline_module/splint *
!!
!!  NAME
!!   splint
!!  USAGE
!!   splint(xarray, yarray, d2y/dx2 array, size of array, x value, y)
!!   splint(xa,ya,y2a,n,x,y)
!!  PURPOSE
!!   NR routine for interpolation using spline table set up by spline
!!   Takes a spline table with uniform x spacing delta_x and, for the
!!   specified value x, finds the y value.
!!  INPUTS
!!   integer :: n (size of table)
!!   real(double) :: x, y, delta_x (x value specified, y value to be found, x-value spacing)
!!   real(double), dimension(n) :: ya, y2a (y, d2y/dx2 arrays)
!!   logical :: flag (flag to warn if given x is out of tabulated range)
!!  USES
!!   datatypes, numbers, GenComms
!!  AUTHOR
!!   L.K.Dash
!!  CREATION DATE
!!   15/04/02
!!  MODIFICATION HISTORY
!!   15:47, 08/04/2003 drb
!!    Fixed major bug: array references were all one element off
!!   10:00, 28/08/2003 drb
!!    Fixed ANOTHER major bug: I'd forgotten that the array reference was used to create the x position so needs one taken off it !
!!   08:56, 2003/11/20 dave
!!    Added return to prevent array overrun
!!  SOURCE
!!
  subroutine splint(delta_x,ya,y2a,n,x,y,flag)

    use datatypes
    use numbers
    use GenComms, ONLY: cq_abort

    ! Passed variables
    integer, intent(in) :: n
    real(double), intent(in) :: x
    real(double), intent(out) :: y
    real(double), intent(in) :: delta_x
    real(double), dimension(n), intent(in) :: y2a,ya
    logical, intent(out) :: flag

    ! Local variables
    integer :: khi,klo
    real(double) :: a,b,xlo,xhi

    ! initialize warning flag
    flag = .false.
    if (x>real(n-1,double)*delta_x .or. x<0) then
       flag = .true. ! the given x lies outside the tabulated values, warn the calling routine
       y = 0.0_double ! Avoid array overruns
       return
    end if

    ! bracket x value
    klo = aint(x/delta_x) + 1 ! integer division gives lower bound, add 1 for array value
    khi = klo + 1
    if(khi>n) call cq_abort("Error in splint: index beyond range ",khi,n)
    ! Remember that we added to get array value and subtract it off again !
    xlo =delta_x*real(klo-1,double)
    xhi =delta_x*real(khi-1,double)

    ! Create value of y
    a = (xhi-x)/delta_x
    b = (x-xlo)/delta_x
    y=a*ya(klo)+b*ya(khi)+ &
         ( (a*a*a-a)*y2a(klo)+(b*b*b-b)*y2a(khi) )*(delta_x*delta_x)/six
    return
  end subroutine splint
!!***

! -----------------------------------------------------------
! Subroutine dsplint
! -----------------------------------------------------------

!!****f* spline_module/dsplint *
!!
!!  NAME
!!   dsplint
!!  USAGE
!!   dsplint(xarray, yarray, d2y/dx2 array, size of array, x value, dy)
!!   dsplint(xa,ya,y2a,n,x,dy)
!!  PURPOSE
!!   NR routine for interpolation using spline table set up by spline
!!   Takes a spline table with uniform x spacing delta_x and, for the
!!   specified value x, finds the derivative value.
!!  INPUTS
!!   integer :: n (size of table)
!!   real(double) :: x, y, delta_x (x value specified, y value to be found, x-value spacing)
!!   real(double), dimension(n) :: ya, y2a (y, d2y/dx2 arrays)
!!   logical :: flag (flag to warn if given x is out of tabulated range)
!!  USES
!!   datatypes, numbers, GenComms
!!  AUTHOR
!!   D.R.Bowler
!!  CREATION DATE
!!   07:50, 2003/03/19 dave
!!  MODIFICATION HISTORY
!!   15:47, 08/04/2003 drb
!!    Fixed major bug: array references were all one element off (carried over from splint)
!!   11:04, 03/09/2003 drb
!!    Added xlo/xhi correction that was added to splint (should have been done earlier !)
!!   08:56, 2003/11/20 dave
!!    Added return to prevent array overrun
!!   08:01, 2003/12/04 dave
!!    Added return of spline-interpolated value as well as derivative
!!  SOURCE
!!
  subroutine dsplint(delta_x,ya,y2a,n,x,y,dy,flag)

    use datatypes
    use numbers
    use GenComms, ONLY: cq_abort

    ! Passed variables
    integer, intent(in) :: n
    real(double), intent(in) :: x
    real(double), intent(out) :: y,dy
    real(double), intent(in) :: delta_x
    real(double), dimension(n), intent(in) :: y2a,ya
    logical, intent(out) :: flag

    ! Local variables
    integer :: khi,klo
    real(double) :: a,b,da,db,dc,dd, xlo, xhi

    ! initialize warning flag
    flag = .false.
    if (x>(n-1)*delta_x.OR.x<0) then
       flag = .true. ! the given x lies outside the tabulated values, warn the calling routine
       y = 0.0_double
       dy = 0.0_double ! Avoid array overruns
       return
    end if

    ! bracket x value
    klo = aint(x/delta_x) + 1 ! integer division gives lower bound, add 1 for array value
    khi = klo + 1
    ! Remember that we added to get array value and subtract it off again !
    xlo =delta_x*real(klo-1,double)
    xhi =delta_x*real(khi-1,double)

    ! Create value of y
    a = (xhi-x)/delta_x
    b = (x-xlo)/delta_x
    da = -one/delta_x
    db =  one/delta_x
    dc = -delta_x*(three*a*a-one)/six
    dd =  delta_x*(three*b*b-one)/six
    y=a*ya(klo)+b*ya(khi)+ &
         ( (a*a*a-a)*y2a(klo)+(b*b*b-b)*y2a(khi) )*(delta_x*delta_x)/six
    dy =  da*ya(klo)+db*ya(khi)+dc*y2a(klo)+dd*y2a(khi)
    return
  end subroutine dsplint
!!***

end module spline_module

max_stars_repo_name: negf/CONQUEST-release

max_issues_repo_name: negf/CONQUEST-release

max_forks_repo_name: negf/CONQUEST-release


text: The Otto and Doris Pleis Bench is one of the benches in the ValleyWise Garden near the Arboretum gazebo.

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
      subroutine err_label(where,what_type,what_label)

      implicit none
      include 'stdunit.h'

      character(len=*), intent(in) ::
     &     where,what_type, what_label

      call quit(1,'err_label for:'//trim(where),
     &     trim(what_type)//
     &     ' label not on list: "'//
     &     trim(what_label)//'"')

      end

max_stars_repo_name: ak-ustutt/GeCCo-public

max_issues_repo_name: ak-ustutt/GeCCo-public

max_forks_repo_name: ak-ustutt/GeCCo-public


program C20211219 !Solving definite integral by rectangular method
implicit none
    integer n , i
    real S,dx,a,b
    real h
    real, external :: rectangle,f,f1,f2,f3,f4,f5,f6

    read(*,*)a,b,n
        dx = (b - a) / n
        S = 0
    do i=1, n
        h = f(a+(i-0.5)*dx)
        S = S + rectangle(h,dx)
    end do
    write(*,"('y=',F8.1)") S

end program

real function rectangle(h,dx)
implicit none
    real,intent(in):: h,dx
    rectangle = h * dx
    return
end function

real function f(x)
implicit none
    real,intent(in):: x
    real, external :: f1,f2,f3,f4,f5,f6
    if(x<0) then
        f = f1(x)
    else if(x>=0 .and. x<1) then
        f = f2(x)
    else if(x>=1 .and. x<4) then
        f = f3(x)
    else if(x>=4 .and. x<16) then
        f = f4(x)
    else if(x>=16 .and. x<64) then
        f = f5(x)
    else if(x>=64) then
        f = f6(x)
    end if

    return
end function

real function f1(x)!x<0
implicit none
    real,intent(in):: x
    f1 = 1 - exp(x)
    return
end function

real function f2(x)!x>=0 .and. x<1
implicit none
    real,intent(in):: x
    f2 = x**4
    return
end function

real function f3(x)!x>=1 .and. x<4
implicit none
    real,intent(in):: x
    f3 = x**3
    return
end function

real function f4(x)!x>=4 .and. x<16
implicit none
    real,intent(in):: x
    f4 = x**2 + 48
    return
end function

real function f5(x)!x>=16 .and. x<64
implicit none
    real,intent(in):: x
    f5 = x + 288
    return
end function

real function f6(x)!x>=64
implicit none
    real,intent(in):: x
    f6 = x**0.5 + 344
    return
end function


module io
! use vars
use maths
implicit none
contains

!------------------------------------------------------------
! Read the main configuration file
!------------------------------------------------------------
	subroutine read_config
	character(len=100) :: config_file
	integer :: nargs, ios

		input_model_file = 'ATOMS/helium.mod'
		input_experiment = 'init_parameters.dat'
		verbose_mode = 0
		linear_solver = 0
		synthesis_mode = 5
		working_mode = 1

	end subroutine read_config

!------------------------------------------------------------
! Read the file with the data for the experiment
!------------------------------------------------------------
	subroutine read_experiment(file,in_params,in_fixed)
	character(len=120) :: file
	type(variable_parameters) :: in_params
	type(fixed_parameters) :: in_fixed
	real(kind=8) :: delta, height
	integer :: j

		open(unit=12,file=file,action='read',status='old')

! Read if we want to include the stimulated emission
		call lb(12,5)
		read(12,*) isti
		if (isti == 1 .and. verbose_mode == 1) then
			print *, 'Including stimulated emission...'
		endif

! Read if we want to include the magnetic field
		call lb(12,2)
		read(12,*) imag
		if (imag == 1 .and. verbose_mode == 1) then
			print *, 'Including magnetic field ...'
		endif

! Read if we want to include the depolarizing collisions
		call lb(12,2)
		read(12,*) idep
		if (idep == 1 .and. verbose_mode == 1) then
			print *, 'Including depolarizing collisions ...'
		endif

! Value of the delta parameter for the depolarizing collisions
		call lb(12,2)
		read(12,*) in_params%delta_collision

! Use Paschen-Back calculation or linear Zeeman effect
		call lb(12,2)
		read(12,*) use_paschen_back
		if (verbose_mode == 1) then
			if (use_paschen_back == 1) then
				print *, 'Taking into account the Paschen-Back effect...'
			else
				print *, 'Not taking into account the Paschen-Back effect...'
			endif
		endif

! Number of slabs
		call lb(12,2)
		read(12,*) in_params%nslabs

! Values characterizing the magnetic field
		call lb(12,2)
		if (in_params%nslabs == 3 .or. in_params%nslabs == -2) then
			read(12,*) in_params%bgauss, in_params%thetabd, in_params%chibd, in_params%bgauss2, in_params%thetabd2, in_params%chibd2
		else
			read(12,*) in_params%bgauss, in_params%thetabd, in_params%chibd
		endif

		if (imag == 1 .and. verbose_mode == 1) then
			print *, 'Field properties'
			write(*,FMT='(4X,A,F8.3,3X,A,F6.2,3X,A,F6.2)') 'B : ', in_params%bgauss, 'thetab : ', in_params%thetabd, 'chib : ',&
				in_params%chibd
		endif

! Values of the height for the calculation of the radiation anisotropy
		call lb(12,2)
		read(12,*) in_params%height

		if (verbose_mode == 1) then
			write(*,FMT='(A,I1,A)') 'Using ', abs(in_params%nslabs), ' slabs'
			if (in_params%nslabs == -2) then
				write(*,FMT='(A)') 'with filling factor'
			endif
		endif

! Value of the optical depth in the maximum (SLAB) or the strength of the line (ME)
		call lb(12,2)
		if (in_params%nslabs == 1) then
			read(12,*) in_params%dtau
		endif

		if (in_params%nslabs == 2 .or. in_params%nslabs == 3) then
			read(12,*) in_params%dtau, in_params%dtau2
		endif

! Two components with filling factor
		if (in_params%nslabs == -2) then
			read(12,*) in_params%dtau, in_params%dtau2, in_params%ff
		endif

! Value of the gradient of the source function (ME)
		call lb(12,2)
		read(12,*) in_params%beta

! Incident Stokes parameters
		call lb(12,2)
		read(12,*) (in_fixed%Stokes_incident(j),j=0,3)

! Transition where to compute the emission
		call lb(12,2)
		read(12,*) in_fixed%nemiss
		if (verbose_mode == 1) then
			write(*,FMT='(A,3X,I3)') ' Calculating emission from line : ', in_fixed%nemiss
		endif

! Use atomic polarization
		call lb(12,2)
		read(12,*) in_fixed%use_atomic_pol
		if (verbose_mode == 1) then
			if (in_fixed%use_atomic_pol == 1) then
				write(*,FMT='(A)') ' Using atomic polarization'
			endif
			if (in_fixed%use_atomic_pol == 0) then
				write(*,FMT='(A)') ' Neglecting atomic polarization'
			endif
			if (in_fixed%use_atomic_pol == 2) then
				write(*,FMT='(A)') ' Neglecting anisotropy'
			endif
		endif

! Observation angle with respect to the vertical
		call lb(12,2)
		read(12,*) in_fixed%thetad, in_fixed%chid, in_fixed%gammad
		if (verbose_mode == 1) then
			print *, 'LOS reference system'
			write(*,FMT='(4X,A,F6.2,3X,A,F6.2,3X,A,F6.2)') 'thetad : ', in_fixed%thetad, 'chid : ', in_fixed%chid, 'gammad : ', &
				in_fixed%gammad
		endif

! Wavelength axis
		call lb(12,2)
		read(12,*) in_fixed%omin, in_fixed%omax, in_fixed%no
		if (verbose_mode == 1) then
			print *, 'Wavelength axis'
			write(*,FMT='(4X,A,F8.3,3X,A,F8.3,3X,A,I5)') 'Min : ', in_fixed%omin, 'Max : ', in_fixed%omax, 'N. steps : ', &
				in_fixed%no
		endif

! Line wavelength and Doppler width
! If two components with different fields, read the two Doppler widths
		call lb(12,2)
		if (in_params%nslabs == 1 .or. in_params%nslabs == 2) then
			read(12,*) in_fixed%wl, in_params%vdopp, in_params%damping
		endif
		if (in_params%nslabs == 3 .or. in_params%nslabs == -2) then
			read(12,*) in_fixed%wl, in_params%vdopp, in_params%vdopp2, in_params%damping
		endif

! Test for the damping treatment. We can set if fixed to a number or estimated using the natural broadening
		in_fixed%damping_treatment = 0
		if (in_params%damping < 0.d0) then
			in_fixed%damping_treatment = 1
		endif

		if (verbose_mode == 1) then
			if (in_params%nslabs <= 2) then
				if (in_fixed%damping_treatment == 1) then
					write(*,FMT='(4X,A,F10.4,3X,A,F6.3,A,F6.3)') 'Wavelength [A]: ', in_fixed%wl, 'Doppler vel. [km/s] : ', &
						in_params%vdopp, 'Damping : automatic'
				else
					write(*,FMT='(4X,A,F10.4,3X,A,F6.3,A,F6.3)') 'Wavelength [A]: ', in_fixed%wl, 'Doppler vel. [km/s] : ', &
						in_params%vdopp, 'Damping : ', in_params%damping
				endif
			else
				if (in_fixed%damping_treatment == 1) then
					write(*,FMT='(4X,A,F10.4,3X,A,F6.3,A,F6.3,A,F6.3)') 'Wavelength [A]: ', in_fixed%wl, 'Doppler vel. [km/s] : ', &
						in_params%vdopp, 'Doppler vel 2. [km/s] : ', in_params%vdopp2, 'Damping : automatic'
				else
					write(*,FMT='(4X,A,F10.4,3X,A,F6.3,A,F6.3,A,F6.3)') 'Wavelength [A]: ', in_fixed%wl, 'Doppler vel. [km/s] : ', &
						in_params%vdopp, 'Doppler vel 2. [km/s] : ', in_params%vdopp2, 'Damping : ', in_params%damping
				endif
			endif
		endif

! Macroscopic velocity
		call lb(12,2)
		if (in_params%nslabs == 1) then
			read(12,*) in_params%vmacro
		else
			read(12,*) in_params%vmacro, in_params%vmacro2
		endif

! Use magneto-optical terms
		call lb(12,2)
		read(12,*) use_mag_opt_RT

! Use stimulated emission terms in the absorption coefficients
		call lb(12,2)
		read(12,*) use_stim_emission_RT

		close(12)

		in_fixed%total_forward_modeling = 0

! Take into account the observation angle (scattering angle) to calculate the height from the apparent height
		if (in_params%height < 0.d0) then
			delta = abs(90.d0 - in_fixed%thetad)
			height = (RSUN + in_params%height) / (cos(delta*PI/180.d0)) - RSUN
			write(*,FMT='(A,F7.3)') 'The height taking into account the scattering angle is : ', height
		endif

	end subroutine read_experiment

!------------------------------------------------------------
! Read the main configuration file
!------------------------------------------------------------
	subroutine read_model_file(file)
	character(len=120) :: file
	integer :: i, j, nJ, n, i1, i2, ir, jmax2, jmin2, jp2, j2, kmin, kmax, k, q, loop
	real(kind=8) :: rnujjp
	character(len=128) :: fileModel(23)

! Generate a scratch file that will be deleted later
		fileModel(1) = '2'
		fileModel(2) = '5'
		fileModel(3) = '1        0               '
		fileModel(4) = '                    0.00 '
		fileModel(5) = '2        2        '
		fileModel(6) = '                    0.00'
		fileModel(7) = '                   -0.987913'
		fileModel(8) = '                   -1.064340'
		fileModel(9) = '3        0       '
		fileModel(10) = '                    0.00'
		fileModel(11) = '4        2         '
		fileModel(12) = '                    0.00'
		fileModel(13) = '                   -0.270647'
		fileModel(14) = '                   -0.292616'
		fileModel(15) = ' 5        4'
		fileModel(16) = '                    0.00'
		fileModel(17) = '                   -0.044187'
		fileModel(18) = '                   -0.046722'
		fileModel(19) = ' 4'
		fileModel(20) = '1    1    2    1.022d7    10829.0911    1.0000000    1.0000000    0.0000000'
		fileModel(21) = '2    1    4    9.478d6    3888.6046    0.2000000    1.0000000    0.0000000'
		fileModel(22) = '3    2    3    2.780d7    7065.7085    1.0000000    1.0000000    0.0000000'
		fileModel(23) = '4    2    5    7.060d7    5875.9663    1.0000000    1.0000000    0.0000000'

		loop = 0

		loop = loop + 1
		is2 = 2

		loop = loop + 1
		n_terms = 5

		if (allocated(lsto2)) deallocate(lsto2)
		allocate(lsto2(n_terms))

		lsto2(1) = 0
		lsto2(2) = 2
		lsto2(3) = 0
		lsto2(4) = 2
		lsto2(5) = 4

		nrhos = 0

! 		open(unit=13,file='mterm.tab',action='write',status='replace')

		jlimit2 = 0
		do i = 1, n_terms
			loop = loop + 1
			jmin2 = abs(is2-lsto2(i))
			jmax2 = is2 + lsto2(i)
			! do j2 = jmin2, jmax2, 2
			! 	loop = loop + 1
			! 	print *, i, j2, fileModel(loop)

			! enddo
			if (verbose_mode == 1) then
				print *, 'Level ', i
				print *, '         Jmin : ', jmin2/2.d0
				print *, '         Jmax : ', jmax2/2.d0
			endif
			jlimit2 = max(jmax2,jlimit2)
		enddo

		if (allocated(energy)) deallocate(energy)
		allocate(energy(n_terms,0:jlimit2))

		energy(1,2) = 0.0
        energy(2,0) = 0.0
        energy(2,2) = -0.987913
        energy(2,4) = -1.064340
        energy(3,2) = 0.0
        energy(4,0) = 0.0
        energy(4,2) = -0.270647
        energy(4,4) = -0.292616
        energy(5,2) = 0.0
        energy(5,4) = -0.044187
        energy(5,6) = -0.046722

		loop = 2

		do i = 1, n_terms

			jmin2 = abs(is2-lsto2(i))
			jmax2 = is2 + lsto2(i)

			do j2 = jmin2, jmax2, 2
				do jp2 = j2, jmax2, 2
					rnujjp = PC * (energy(i,j2)-energy(i,jp2))
					kmin = (jp2-j2) / 2
					kmax = (jp2+j2) / 2
					do k = kmin, kmax
						do q = -k, k
							if (j2 == jp2 .and. q < 0) then
							else
								do ir = 1, 2
									if (j2 == jp2 .and. q == 0 .and. ir == 2) then
									else
										nrhos = nrhos + 1
									endif
								enddo
							endif
						enddo
					enddo
				enddo
			enddo
		enddo

		loop = 2

		if (allocated(ntab)) deallocate(ntab)
		if (allocated(j2tab)) deallocate(j2tab)
		if (allocated(jp2tab)) deallocate(jp2tab)
		if (allocated(ktab)) deallocate(ktab)
		if (allocated(qtab)) deallocate(qtab)
		if (allocated(irtab)) deallocate(irtab)
		if (allocated(rnutab)) deallocate(rnutab)

		allocate(ntab(nrhos))
		allocate(j2tab(nrhos))
		allocate(jp2tab(nrhos))
		allocate(ktab(nrhos))
		allocate(qtab(nrhos))
		allocate(irtab(nrhos))
		allocate(rnutab(nrhos))

		nrhos = 0

		do i = 1, n_terms
			n = i

			jmin2 = abs(is2-lsto2(i))
			jmax2 = is2 + lsto2(i)

			do j2 = jmin2, jmax2, 2
				do jp2 = j2, jmax2, 2
					rnujjp = PC * (energy(i,j2)-energy(i,jp2))
					kmin = (jp2-j2) / 2
					kmax = (jp2+j2) / 2
					do k = kmin, kmax
						do q = -k, k
							if (j2 == jp2 .and. q < 0) then
							else
								do ir = 1, 2
									if (j2 == jp2 .and. q == 0 .and. ir == 2) then
									else
										nrhos = nrhos + 1
										ntab(nrhos) = n
										j2tab(nrhos) = j2
										jp2tab(nrhos) = jp2
										ktab(nrhos) = k
										qtab(nrhos) = q
										irtab(nrhos) = ir
										rnutab(nrhos) = rnujjp
! 										write(13,FMT='(7(1X,I5),6X,E12.5)') nrhos, n, j2, jp2, k, q, ir, rnujjp
									endif
								enddo
							endif
						enddo
					enddo
				enddo
			enddo
		enddo

		! Now read transitions
		! loop = loop + 1
		! read(fileModel(loop),*) atom%ntran

		atom%ntran = 4

		allocate(atom%nterml(atom%ntran))
		allocate(atom%ntermu(atom%ntran))
		allocate(atom%ae(atom%ntran))
		allocate(atom%wavelength(atom%ntran))
		allocate(atom%reduction_factor(atom%ntran))
		allocate(atom%reduction_factor_omega(atom%ntran))
		allocate(atom%j10(atom%ntran))

		atom%nterml(1) = 1
		atom%ntermu(1) = 2
		atom%ae(1) = 1.022d7
		atom%wavelength(1) = 10829.0911d0
		atom%reduction_factor(1) = 1.d0
		atom%reduction_factor_omega(1) = 1.d0
		atom%j10(1) = 0.d0

		atom%nterml(2) = 1
		atom%ntermu(2) = 4
		atom%ae(2) = 9.478d6
		atom%wavelength(2) = 3888.6046d0
		atom%reduction_factor(2) = 0.2d0
		atom%reduction_factor_omega(2) = 1.0d0
		atom%j10(2) = 0.0d0

		atom%nterml(3) = 2
		atom%ntermu(3) = 3
		atom%ae(3) = 2.780d7
		atom%wavelength(3) = 7065.7085d0
		atom%reduction_factor(3) = 1.0d0
		atom%reduction_factor_omega(3) = 1.0d0
		atom%j10(3) = 0.0d0

		atom%nterml(4) = 2
		atom%ntermu(4) = 5
		atom%ae(4) = 7.060d7
		atom%wavelength(4) = 5875.9663d0
		atom%reduction_factor(4) = 1.0d0
		atom%reduction_factor_omega(4) = 1.0d0
		atom%j10(4) = 0.0d0

		! do i = 1, atom%ntran
		! 	loop = loop + 1
		! 	read(fileModel(loop),*) k, atom%nterml(i), atom%ntermu(i), atom%ae(i), atom%wavelength(i), &
		! 		atom%reduction_factor(i), atom%reduction_factor_omega(i), atom%j10(i)
		! enddo

		if (verbose_mode == 1) then
			print *, 'Number of unknowns : ', nrhos
		endif

		if (verbose_mode == 1) then
			print *, 'Number of unknowns : ', nrhos
		endif

	end subroutine read_model_file

!------------------------------------------------------------
! Clean all the allocated memory
!------------------------------------------------------------
	subroutine clean

		deallocate(lsto2)
		deallocate(energy)
		deallocate(ntab)
		deallocate(ktab)
		deallocate(qtab)
		deallocate(irtab)
		deallocate(j2tab)
		deallocate(jp2tab)
		deallocate(rnutab)

	end subroutine clean

!------------------------------------------------------------
! Number of columns in a line
!------------------------------------------------------------
	function ncolumns(string)
	integer :: ncolumns
	character(len=*) :: string
	character(len=1) :: ch, chprev
	integer :: i, l, blanks
		l = len(trim(string))
		blanks = 1
		chprev = string(1:1)
		do i = 2, l
			ch = string(i:i)
			if (ch /= ' ' .and. chprev == ' ') then
				blanks = blanks + 1
			endif
			chprev = ch
		enddo
		ncolumns = blanks - 1
	end function ncolumns
end module io


text:
      MODULE write_mgrid
      USE biotsavart
      USE mgrid_mod, ONLY: vn_nextcur, vn_mgmode, vn_ir,                       &
     &                     vn_jz, vn_kp, vn_nfp, vn_rmin, vn_rmax,             &
     &                     vn_zmin, vn_zmax, vn_coilgrp, vn_coilcur,           &
     &                     vn_br0, vn_bz0, vn_bp0, vn_ar0, vn_az0,             &
     &                     vn_ap0
      IMPLICIT NONE
C-----------------------------------------------
C   L o c a l   P a r a m e t e r s
C-----------------------------------------------
      REAL (rprec), PARAMETER      :: one = 1
      CHARACTER (LEN=*), PARAMETER ::                                          &
     &   coildim(2) = (/'stringsize          ', 'external_coil_groups'/)
      CHARACTER (LEN=*), PARAMETER :: groupdim(1)= (/'external_coils'/)
      CHARACTER (LEN=*), PARAMETER :: cylcoord(3)= (/'rad','zee','phi'/)
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER                                       :: ir = 121
      INTEGER                                       :: jz = 121
      INTEGER                                       :: iextc
      INTEGER                                       :: nfp
      INTEGER                                       :: nextcur
      INTEGER                                       :: istat
      INTEGER                                       :: kp
      INTEGER                                       :: kp2
      INTEGER                                       :: jz2
      INTEGER                                       :: kp_odd
      INTEGER                                       :: jz_odd
      REAL (rprec), ALLOCATABLE, DIMENSION(:, :, :) :: br
      REAL (rprec), ALLOCATABLE, DIMENSION(:, :, :) :: bz
      REAL (rprec), ALLOCATABLE, DIMENSION(:, :, :) :: bp
      REAL (rprec), ALLOCATABLE, DIMENSION(:, :, :) :: ar
      REAL (rprec), ALLOCATABLE, DIMENSION(:, :, :) :: az
      REAL (rprec), ALLOCATABLE, DIMENSION(:, :, :) :: ap
      REAL (rprec), ALLOCATABLE, DIMENSION(:)       :: extcur
      REAL (rprec)                                  :: rmin
      REAL (rprec)                                  :: rmax
      REAL (rprec)                                  :: zmin
      REAL (rprec)                                  :: zmax
      REAL (rprec)                                  :: fperiod
      REAL (rprec)                                  :: delr
      REAL (rprec)                                  :: delz
      REAL (rprec)                                  :: delp
      LOGICAL                                       ::                         &
     &   lstell_sym = .false.
      CHARACTER (LEN=1)                             :: mgrid_mode = 'S'
      CHARACTER (LEN=70)                            :: mgrid_file
      CHARACTER (LEN=70)                            :: coil_file
      CHARACTER (LEN=60)                            :: mgrid_ext
      LOGICAL                                       :: use_eddy
      CHARACTER (len=30), DIMENSION(:), ALLOCATABLE :: temp_sname

      PRIVATE                                       :: istat
C-----------------------------------------------

      CONTAINS
#if defined(NETCDF)
      SUBROUTINE write_mgrid_nc
      USE ezcdf
      IMPLICIT NONE
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER                                        :: ngrid
      INTEGER                                        :: ig
      CHARACTER (LEN=100)                            :: temp
      CHARACTER (LEN=100), ALLOCATABLE, DIMENSION(:) :: vn_br
      CHARACTER (LEN=100), ALLOCATABLE, DIMENSION(:) :: vn_bp
      CHARACTER (LEN=100), ALLOCATABLE, DIMENSION(:) :: vn_bz
      CHARACTER (LEN=100), ALLOCATABLE, DIMENSION(:) :: vn_ar
      CHARACTER (LEN=100), ALLOCATABLE, DIMENSION(:) :: vn_ap
      CHARACTER (LEN=100), ALLOCATABLE, DIMENSION(:) :: vn_az
C-----------------------------------------------
      mgrid_file = TRIM(mgrid_file) // '.nc'

      CALL cdf_open(ngrid,mgrid_file,'w',istat)
      IF (istat .ne. 0) THEN
         STOP 'Error opening WOUT.nc file VMEC WROUT'
      END IF

!
!     DEFINE DATA VARIABLES, DIMENSION NAMES
!
      CALL cdf_define(ngrid, vn_ir, ir)
      CALL cdf_define(ngrid, vn_jz, jz)
      CALL cdf_define(ngrid, vn_kp, kp)
      CALL cdf_define(ngrid, vn_nfp, nfp)
      CALL cdf_define(ngrid, vn_nextcur, nextcur)
      CALL cdf_define(ngrid, vn_rmin, rmin)
      CALL cdf_define(ngrid, vn_zmin, zmin)
      CALL cdf_define(ngrid, vn_rmax, rmax)
      CALL cdf_define(ngrid, vn_zmax, zmax)
      IF (nextcur .eq. 1) THEN
         CALL cdf_define(ngrid, vn_coilgrp,coil_group(1)%s_name)
      ELSE IF (use_eddy) THEN
         ALLOCATE(temp_sname(nextcur))
         temp_sname(1:nextcur - 1) = coil_group(1:nextcur - 1)%s_name
         temp_sname(nextcur) = 'Eddy currents'
         CALL cdf_define(ngrid, vn_coilgrp, temp_sname, dimname=coildim)
      ELSE
         CALL cdf_define(ngrid, vn_coilgrp,coil_group(1:nextcur)%s_name,
     &                   dimname=coildim)
      END IF
      CALL cdf_define(ngrid, vn_mgmode, mgrid_mode)
      CALL cdf_define(ngrid, vn_coilcur, extcur(1:nextcur),
     &                dimname=groupdim)
!
!     STORED AS 3D ARRAYS (ACTUALLY 4D, BUT CUT THROUGH IG)
!
      ALLOCATE (vn_br(nextcur), vn_bz(nextcur), vn_bp(nextcur))
      ALLOCATE (vn_ar(nextcur), vn_az(nextcur), vn_ap(nextcur))

      DO ig = 1, nextcur
         write (temp, '(a,i3.3)') "_",ig
         vn_br(ig) = vn_br0 // temp
         vn_bp(ig) = vn_bp0 // temp
         vn_bz(ig) = vn_bz0 // temp
         CALL cdf_define(ngrid, vn_br(ig), br, dimname=cylcoord)
         CALL cdf_define(ngrid, vn_bp(ig), bp, dimname=cylcoord)
         CALL cdf_define(ngrid, vn_bz(ig), bz, dimname=cylcoord)

         vn_ar(ig) = vn_ar0 // temp
         vn_ap(ig) = vn_ap0 // temp
         vn_az(ig) = vn_az0 // temp
         CALL cdf_define(ngrid, vn_ar(ig), ar, dimname=cylcoord)
         CALL cdf_define(ngrid, vn_ap(ig), ap, dimname=cylcoord)
         CALL cdf_define(ngrid, vn_az(ig), az, dimname=cylcoord)
      END DO

!
!     WRITE OUT DATA
!
      CALL cdf_write(ngrid, vn_ir, ir)
      CALL cdf_write(ngrid, vn_jz, jz)
      CALL cdf_write(ngrid, vn_kp, kp)
      CALL cdf_write(ngrid, vn_nfp, nfp)
      CALL cdf_write(ngrid, vn_nextcur, nextcur)
      CALL cdf_write(ngrid, vn_rmin, rmin)
      CALL cdf_write(ngrid, vn_zmin, zmin)
      CALL cdf_write(ngrid, vn_rmax, rmax)
      CALL cdf_write(ngrid, vn_zmax, zmax)
      IF (nextcur .eq. 1) THEN
         CALL cdf_write(ngrid, vn_coilgrp, coil_group(1)%s_name)
      ELSE IF (use_eddy) THEN
         CALL cdf_write(ngrid, vn_coilgrp, temp_sname)
         DEALLOCATE(temp_sname)
      ELSE
         CALL cdf_write(ngrid, vn_coilgrp, coil_group(1:nextcur)%s_name)
      END IF

!
!     SET UP CYLINDRICAL COMPONENTS OF MAGNETIC FIELD ON GRID
!     SUM OVER CURRENT GROUPS IG = 1,NEXTCUR
!     NOTE: USER MUST SUPPLY SUBROUTINE "	" TO COMPUTE THESE VALUES
!
      GROUPS: DO ig = 1, nextcur
         IF (ig .eq. nextcur .and. use_eddy) THEN
            br = 0
            bp = 0
            bz = 0
            ar = 0
            ap = 0
            az = 0
            extcur(nextcur) = 1
         ELSE
            CALL compute_bfield(ig)
         END IF

         CALL cdf_write(ngrid, vn_br(ig), br)
         CALL cdf_write(ngrid, vn_bp(ig), bp)
         CALL cdf_write(ngrid, vn_bz(ig), bz)

         CALL cdf_write(ngrid, vn_ar(ig), ar)
         CALL cdf_write(ngrid, vn_ap(ig), ap)
         CALL cdf_write(ngrid, vn_az(ig), az)

      END DO GROUPS

      CALL cdf_write(ngrid, vn_mgmode, mgrid_mode)
      CALL cdf_write(ngrid, vn_coilcur, extcur(1:nextcur))

      CALL cdf_close(ngrid)

      DEALLOCATE (vn_br, vn_bz, vn_bp)
      DEALLOCATE (vn_ar, vn_az, vn_ap)

      END SUBROUTINE write_mgrid_nc
#endif

      SUBROUTINE write_mgrid_bin
      USE safe_open_mod
      IMPLICIT NONE
C-----------------------------------------------
C   L o c a l   P a r a m e t e r s
C-----------------------------------------------
      INTEGER, PARAMETER :: igrid0 = 10
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER            :: igrid
      INTEGER            :: ig
      INTEGER            :: i
C-----------------------------------------------

      igrid = igrid0
      mgrid_file = TRIM(mgrid_file) // '.bin'
      CALL safe_open(igrid, istat, TRIM(mgrid_file),                           &
     &               'replace', 'unformatted')

      IF (istat .ne. 0) THEN
         WRITE (6,*) 'XGRID could not create ', TRIM(mgrid_file)
         WRITE (6,*) 'IOSTAT = ', istat,' IUNIT = ', igrid
         STOP 20
      END IF

!
!     MARK NEW-STYLE CODE WITH NEXTCUR < 0 SO READ_MGRID CAN DISTINGUISH
!     FROM OLD-STYLE MGRID FILE
!
      WRITE (igrid) ir, jz, kp, nfp, -nextcur
      WRITE (igrid) rmin, zmin, rmax, zmax
      WRITE (igrid) (coil_group(i)%s_name, i = 1, nextcur)

!
!     SET UP CYLINDRICAL COMPONENTS OF MAGNETIC FIELD ON GRID
!     SUM OVER CURRENT GROUPS IG = 1,NEXTCUR
!     NOTE: USER MUST SUPPLY SUBROUTINE "BFIELD" TO COMPUTE THESE VALUES
!
      GROUPS: DO ig = 1, nextcur

         CALL compute_bfield(ig)

         WRITE (igrid, iostat=istat) br, bp, bz
         WRITE (igrid, iostat=istat) ar, ap, az
!        WRITE(igrid,'(1p,9e12.4)') (((br(i,j,k), bz(i,j,k),              !OLD STYLE
!    1         bp(i,j,k), i=1,ir), j=1,jz),k=1,kp)
         IF (istat .ne. 0) THEN
            STOP ' Error writing bfield components'
         END IF

      END DO GROUPS

      WRITE (igrid) mgrid_mode
      WRITE (igrid) extcur(1:nextcur)

!
!     ADD RECONSTRUCTION/OBSERVATION STUFF HERE
!

      CLOSE (igrid)

      END SUBROUTINE write_mgrid_bin

      SUBROUTINE compute_bfield(ig)
      IMPLICIT NONE
      INTEGER, INTENT(IN) :: ig
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER             :: icoil
      INTEGER             :: numcoils
      INTEGER             :: k
      INTEGER             :: j
      INTEGER             :: i
      INTEGER             :: numfils
      INTEGER             :: curindex(1)
      REAL (rprec)        :: rad
      REAL (rprec)        :: zee
      REAL (rprec)        :: phi
      REAL (rprec)        :: extcur_ig
      REAL (rprec)        :: extcur_ig_inv
C-----------------------------------------------
      numcoils = coil_group(ig)%ncoil
      curindex = MAXLOC(abs(coil_group(ig)%coils(1:numcoils)%current))
      extcur_ig = coil_group(ig)%coils(curindex(1))%current

      extcur(ig) = extcur_ig
      IF (extcur_ig .eq. zero) THEN
         WRITE (6, 100) ig
      END IF
 100  FORMAT ('In COILS file, current(s) vanished for coil group ',i4)

      numfils = 0
      DO icoil = 1, numcoils
         numfils = numfils +                                                   &
     &      MAX(SIZE(coil_group(ig)%coils(icoil)%xnod, 2), 2) - 1
!
!        Compute field for unit current
!
!  JDH 2011-08-16. Comment out below - don't want to change bsc_coil data
!    Adjust B values late on with Raw test
!    Will behave differently if extcur_ig == 0.
!         IF (extcur_ig .ne. zero) THEN
!            coil_group(ig) % coils(icoil) % current =
!     1      coil_group(ig) % coils(icoil) % current /extcur_ig
!         ELSE
!            coil_group(ig) % coils(icoil) % current = 1
!         END IF
      END DO

      WRITE (6, '(/,2a)')      ' COIL GROUP          :  ',                     &
     &                         TRIM(coil_group(ig) % s_name)
      WRITE (6, '(a,i6,a,i6)') ' TOTAL COILS IN GROUP: ', numcoils,            &
     &                         ' TOTAL FILAMENTS: ', numfils

      k = 1                     ! this is always a symmetry plane
      phi = (k - 1)*delp

!$omp parallel do private(i, j, zee, rad)
      DO j = 1, jz2 + jz_odd
         zee = zmin + (j - 1)*delz
         DO i = 1, ir
            rad = rmin + (i - 1)*delr
            CALL bfield(rad, phi, zee, br(i,j,k),
     &                  bp(i,j,k), bz(i,j,k), ig)

            IF (lstell_sym) THEN
               br(i,jz + 1 - j,k) = -br(i,j,k)
               bz(i,jz + 1 - j,k) =  bz(i,j,k)
               bp(i,jz + 1 - j,k) =  bp(i,j,k)
            END IF

            CALL afield (rad, phi, zee, ar(i,j,k),                             &
     &                   ap(i,j,k), az(i,j,k), ig)

            IF (lstell_sym) THEN
               ar(i,jz + 1 - j,k) = -ar(i,j,k)
               az(i,jz + 1 - j,k) =  az(i,j,k)
               ap(i,jz + 1 - j,k) =  ap(i,j,k)
            END IF

         END DO
      END DO
!$omp end parallel do
      IF (kp .ne. 1) THEN
         WRITE (6, '(a,i4,a,i4,a)') ' K = ',k,' (OUT OF ',KP,')'
      END IF

!$omp parallel do private(i, j, k, zee, rad, phi)
      DO k = 2, kp2 + kp_odd
         phi = (k - 1)*delp
         DO j = 1, jz
            zee = zmin + (j - 1)*delz
            DO i = 1, ir
               rad = rmin + (i - 1)*delr
               CALL bfield(rad, phi, zee, br(i,j,k),                           &
     &                     bp(i,j,k), bz(i,j,k), ig)
               IF (lstell_sym) THEN
                  br(i,jz + 1 - j,kp + 2 - k) = -br(i,j,k)
                  bz(i,jz + 1 - j,kp + 2 - k) =  bz(i,j,k)
                  bp(i,jz + 1 - j,kp + 2 - k) =  bp(i,j,k)
               END IF

               CALL afield(rad, phi, zee, ar(i,j,k),                           &
     &                     ap(i,j,k), az(i,j,k), ig)
               IF (lstell_sym) THEN
                  ar(i,jz + 1 - j,kp + 2 - k) = -ar(i,j,k)
                  az(i,jz + 1 - j,kp + 2 - k) =  az(i,j,k)
                  ap(i,jz + 1 - j,kp + 2 - k) =  ap(i,j,k)
               END IF
            END DO
         END DO
         WRITE (6,'(a,i4)') ' K = ',k
      END DO
!$omp end parallel do

      IF ((kp_odd .eq. 0) .and. lstell_sym) THEN       ! another symmetry plane
         k = kp2 + 1
         phi = (k - 1)*delp
!$omp parallel do private(i, j, zee, rad)
         DO j = 1, jz2 + jz_odd
            zee = zmin + (j - 1)*delz
            DO i = 1, ir
               rad = rmin + (i - 1)*delr
               CALL bfield(rad, phi, zee, br(i,j,k),                           &
     &                     bp(i,j,k), bz(i,j,k), ig)
               br(i,jz + 1 - j,k) = -br(i,j,k)
               bz(i,jz + 1 - j,k) =  bz(i,j,k)
               bp(i,jz + 1 - j,k) =  bp(i,j,k)

               CALL afield (rad, phi, zee, ar(i,j,k),                          &
     &                      ap(i,j,k), az(i,j,k), ig)
               ar(i,jz + 1 - j,k) = -ar(i,j,k)
               az(i,jz + 1 - j,k) =  az(i,j,k)
               ap(i,jz + 1 - j,k) =  ap(i,j,k)
            END DO
         END DO
!$omp end parallel do
         WRITE (6,'(a,i4)') ' K = ',k
      END IF

      IF (mgrid_mode .eq. 'R') THEN
!  JDH 2011-080-16. Comment out below (1 line)
!         br = br*extcur_ig;  bp = bp*extcur_ig; bz = bz*extcur_ig    ! "Raw" fields
         IF (ig .lt. 10) THEN
            WRITE (iextc, 220) ig, extcur_ig, extcur_ig      ! currents for "raw" mode
         END IF
         IF (ig .ge. 10) THEN
            WRITE (iextc, 225) ig, extcur_ig, extcur_ig
         END IF
      ELSE
!  JDH 2011-08-16. Add below, scale b. 2011-08-24 - clean up extcur_ig_inv
         IF (ABS(extcur_ig) .gt. 1.E-100_rprec) THEN
            extcur_ig_inv = 1./extcur_ig
         ELSE
            extcur_ig_inv = 1.E100_rprec
         ENDIF
!         extcur_ig_inv = 1. / MAX(extcur_ig, 1.E-100)
         br = br*extcur_ig_inv
         bp = bp*extcur_ig_inv
         bz = bz*extcur_ig_inv                                 ! "Scaled" fields

         ar = ar*extcur_ig_inv
         ap = ap*extcur_ig_inv
         az = az*extcur_ig_inv                                 ! "Scaled" fields
!  JDH 2011-08-16. end addition
         IF (ig .lt. 10) THEN
            WRITE (iextc, 210) ig, extcur_ig      ! currents for "scaled" mode
         END IF
         IF (ig .ge. 10) THEN
            WRITE (iextc, 215) ig, extcur_ig
         END IF
      END IF

 210  FORMAT('EXTCUR(', i1,')  = ', 1p,e22.14)
 215  FORMAT('EXTCUR(', i2,') = ', 1p,e22.14)
 220  FORMAT('EXTCUR(', i1,')  = ', 1p,e22.14,'/',e22.14)
 225  FORMAT('EXTCUR(', i2,') = ', 1p,e22.14,'/',e22.14)

      END SUBROUTINE compute_bfield

      END MODULE write_mgrid

max_stars_repo_name: ORNL-Fusion/MAKEGRID

max_issues_repo_name: ORNL-Fusion/MAKEGRID

max_forks_repo_name: ORNL-Fusion/MAKEGRID


text:
C    Copyright(C) 1999-2020 National Technology & Engineering Solutions
C    of Sandia, LLC (NTESS).  Under the terms of Contract DE-NA0003525 with
C    NTESS, the U.S. Government retains certain rights in this software.
C
C    See packages/seacas/LICENSE for details

      SUBROUTINE SELTIM (TIMES, ITMSEL)
      DIMENSION TIMES(*)
      LOGICAL ITMSEL(*)
      CHARACTER*16 ENGNOT, STRA, STRB
      CHARACTER*80 STRTMP
      EXTERNAL ENGNOT

      include 'nu_ptim.blk'

C ... TOLER is the tolerance for matching a timestep.  If the difference
C     is less than TOLER, then a match occurs.  TOLER1 = TOLER + 1

      PARAMETER (TOLER = 1.0E-3)
      TOLERP1 = 1.0 + TOLER
      TOLERM1 = 1.0 - TOLER

      CALL INILOG (NSTEP, .FALSE., ITMSEL)

      NLAST  = 0
      LSTSEL = NSTEP
      TIMGET = STMIN
      NUMSEL = 0

      IF (STDEL .EQ. 0.0) THEN
   10    CONTINUE
         NLAST = NLAST + 1
         IF (NLAST .GT. NSTEP) GO TO 30
         IF (TIMES(NLAST) .LE. TOLERP1 * STMAX .AND.
     *       TIMES(NLAST) .GE. TOLERM1 * STMIN ) THEN
            ITMSEL(NLAST) = .TRUE.
            NUMSEL = NUMSEL + 1
            LSTSEL = NLAST
         ELSE IF (TIMES(NLAST) .GT. TOLERP1 * STMAX) THEN
            LSTSEL = NLAST - 1
            GO TO 30
         END IF
         GO TO 10
      ELSE
   20    CONTINUE
         NLAST = NLAST + 1
         IF (NLAST .GT. NSTEP) GO TO 30
         TDELT = ABS(TIMES(NLAST) - TIMGET)
         IF (TDELT .LT. (TOLER * TMAX) .AND.
     *      TIMES(NLAST) .LE. TOLERP1 * STMAX) THEN
            ITMSEL(NLAST) = .TRUE.
            NUMSEL = NUMSEL + 1
            TIMGET = MIN ( TIMES(NLAST) + STDEL, TMAX)
            LSTSEL = NLAST
         ELSE IF (TIMES(NLAST) .GT. TOLERP1 * STMAX) THEN
            LSTSEL = NLAST - 1
            GO TO 30
         END IF
         GO TO 20
      END IF
   30 CONTINUE

      IF (NUMSEL .EQ. 0) THEN
         CALL PRTERR ('WARNING', 'No time steps selected.')
      ELSE
        STRA = ENGNOT(STMIN,2)
        STRB = ENGNOT(STMAX,2)
        WRITE (STRTMP, 40) NUMSEL, STRA, STRB
        CALL SQZSTR(STRTMP, LSTR)
        WRITE (*, 50) STRTMP(:LSTR)
      END IF
      RETURN
   40 FORMAT (I5,' Steps Selected from ',A16,' to ',A16)
   50 FORMAT (/,5X,A)
      END

max_stars_repo_name: jschueller/seacas

max_issues_repo_name: jschueller/seacas

max_forks_repo_name: jschueller/seacas


text:
! $Id$
! -----------------------------------------------------------
! Module mult_init_module
! -----------------------------------------------------------
! Code area 2: matrices
! -----------------------------------------------------------

!!****h* Conquest/mult_init_module *
!!  NAME
!!   mult_init_module
!!  PURPOSE
!!   Routines to initialise matrix multiplications.
!!   Contains mult_ini and check_mm
!!  AUTHOR
!!   D.R.Bowler
!!  CREATION DATE
!!   25/11/99
!!  MODIFICATION HISTORY
!!   27/01/00 by D.R.Bowler to use new communications
!!   19/02/00 by DRB to implement transpose initialisation
!!   02/03/00 by DRB to finish transpose work
!!   12/04/00 by DRB to add itran_addr to trans_remote
!!   09/06/00 by DRB to adapt for Conquest
!!   20/06/2001 dave
!!    Added ROBODoc headers, RCS Id and Log tags and cq_abort
!!***
module mult_init_module

contains

!!****f* mult_init_module/mult_ini *
!!
!!  NAME
!!   mult_ini
!!  USAGE
!!
!!  PURPOSE
!!   Initiates a multiplication given by the
!!   matrices in a_b_c.
!!  INPUTS
!!
!!
!!  USES
!!
!!  AUTHOR
!!   D.R.Bowler
!!  CREATION DATE
!!   25/11/99
!!  MODIFICATION HISTORY
!!   20/06/2001 dave
!!    Added ROBODoc header and cq_abort
!!  SOURCE
!!
  subroutine mult_ini(a_b_c,ind,myid,n_prim,parts)

    ! Module usage
    use GenComms, ONLY: cq_abort
    use matrix_module
    use basic_types
    use matrix_comms_module

    ! Passed variables
    integer :: n_prim,myid,mx_part
    type(matrix_mult)::a_b_c
    type(group_set) :: parts
    integer, dimension(:), target :: ind

    ! Local variables
    integer :: ierr, indlen

    indlen = size(ind)
    a_b_c%bindex => ind(1:indlen)
    ! Check for mult_type
    if(a_b_c%mult_type/=1.AND.a_b_c%mult_type/=2) then
       call cq_abort('mult_ini: mult_type must be 1 or 2: ',a_b_c%mult_type)
    endif
    ! Check that the dimension are all correct
    if(n_prim.gt.0)then
       call check_mm( a_b_c%ahalo%np_in_halo,&
            a_b_c%ahalo%lab_hcell,a_b_c%ahalo%j_beg,a_b_c%ahalo%nh_part,&
            a_b_c%ahalo%j_seq,a_b_c%ahalo%i_halo,&
            a_b_c%ahalo%mx_part, &
            parts%mx_gcell,a_b_c%ahalo%mx_halo,parts%mx_mem_grp)
       call check_mm( a_b_c%chalo%np_in_halo,&
            a_b_c%chalo%lab_hcell,a_b_c%chalo%j_beg,a_b_c%chalo%nh_part,&
            a_b_c%chalo%j_seq,a_b_c%chalo%i_halo,&
            a_b_c%chalo%mx_part, &
            parts%mx_gcell,a_b_c%chalo%mx_halo,parts%mx_mem_grp)
    endif
    call init_mult_comms(parts,a_b_c,myid)
    return
  end subroutine mult_ini
!!***

!!****f* mult_init_module/check_mm *
!!
!!  NAME
!!   check_mm
!!  USAGE
!!
!!  PURPOSE
!!   Checks limits for a matrix multiplication
!!  INPUTS
!!
!!
!!  USES
!!
!!  AUTHOR
!!   D.R.Bowler
!!  CREATION DATE
!!   25/11/99
!!  MODIFICATION HISTORY
!!   20/06/2001 dave
!!    Added ROBODoc header and cq_abort
!!  SOURCE
!!
  subroutine check_mm(np_in_ahalo, lab_ahcell,kbeg,nah_part,kseq,k_halo,&
       mx_apart,mx_pcell, mx_ahalo,mx_part)

    ! Module usage
    use datatypes
    use GenComms, ONLY: cq_abort

    implicit none

    ! Passed variables
    integer :: mx_apart,mx_pcell
    integer :: mx_ahalo,mx_part
    integer :: np_in_ahalo
    integer :: lab_ahcell(mx_apart),kbeg(mx_apart),nah_part(mx_apart)
    integer :: kseq(mx_ahalo),k_halo(mx_ahalo)

    ! Local variables
    integer :: kpart,k

    if((np_in_ahalo<=0).or.(np_in_ahalo>mx_apart)) then
       call cq_abort('check_mm: np_in_halo wrong')
    endif
    do kpart=1,np_in_ahalo
       if(lab_ahcell(kpart)>mx_pcell) then
          call cq_abort('check_mm: lab_ahcell wrong')
       endif
       if(kbeg(kpart)+nah_part(kpart)-1>mx_ahalo) then
          call cq_abort('check_mm: kseq arg wrong')
       endif
       if(nah_part(kpart)<=0) then
          call cq_abort('check_mm: nahpart < 1')
       endif
       do k=1,nah_part(kpart)
          if(kseq(kbeg(kpart)+k-1)>mx_part) then
             call cq_abort('check_mm: ibind_rem wrong')
          endif
          if(k_halo(kbeg(kpart)+k-1)>mx_ahalo) then
             call cq_abort('check_mm: nahnab arg wrong')
          endif
       enddo
    enddo
    return
  end subroutine check_mm
!!***
end module mult_init_module


text:
! { dg-do compile }
! Tests the implementation of compile-time shape testing, required to fix
! PR19015.  The functionality of maxloc and friends is tested by existing
! testcases.
!
! Contributed by Thomas Koeing  <Thomas.Koenig@online.de>
!
  integer, dimension(0:1,0:1) :: n
  integer, dimension(1) :: i
  n = reshape((/1, 2, 3, 4/), shape(n))
  i = maxloc(n) ! { dg-error "Different shape for array assignment" }
  i = maxloc(n,dim=1) ! { dg-error "Different shape for array assignment" }
!  print *,i
end program

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
*
*  Originally derived from LAPACK 3.0 test suite failure.
*
*  David Billinghurst, (David.Billinghurst@riotinto.com.au)
*  23 February 2000
*
      INTEGER N, I, SLASQX
      N = 20
      I = SLASQX( N )
      IF ( I .NE. 2*N ) THEN
         WRITE(6,*) 'I = ', I, ' but should be ', 2*N
         CALL ABORT()
      END IF
      END

      INTEGER FUNCTION SLASQX( N )
      INTEGER  N, I0, I, K
      I0 = 1
      DO I = 4*I0, 2*( I0+N-1 ), 4
         K = I
      END DO
      SLASQX = K
      RETURN
      END


text:
!!  Copyright (C)  Stichting Deltares, 2012-2016.
!!
!!  This program is free software: you can redistribute it and/or modify
!!  it under the terms of the GNU General Public License version 3,
!!  as published by the Free Software Foundation.
!!
!!  This program is distributed in the hope that it will be useful,
!!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!!  GNU General Public License for more details.
!!
!!  You should have received a copy of the GNU General Public License
!!  along with this program. If not, see <http://www.gnu.org/licenses/>.
!!
!!  contact: delft3d.support@deltares.nl
!!  Stichting Deltares
!!  P.O. Box 177
!!  2600 MH Delft, The Netherlands
!!
!!  All indications and logos of, and references to registered trademarks
!!  of Stichting Deltares remain the property of Stichting Deltares. All
!!  rights reserved.

      subroutine caltau ( pmsa   , fl     , ipoint , increm , noseg  ,
     &                    noflux , iexpnt , iknmrk , noq1   , noq2   ,
     &                    noq3   , noq4   )
!>\file
!>       Calculation of bottom friction

!
!     Description of the module :
!
!        General water quality module for DELWAQ:
!        CALCULATE SHEER STRESS, CAUSED BY FLOW, WIND AND HUMAN ACIVITYS
!        (SHIPS, DREDGING). MODEL IS VERTICALLY AVERAGED.
!
!        THE SUBROUTINE ALSO PRODUCES A CALCULATED STREAM VELOCITY
!        BASED ON TAU-TOTAL FOR COMPARISON WITH CRITICAL STREAM VELOCITIES
!        IN STEAD OF CRITICAL TAUS (DELTA STUDY, T692).
!
!        AVERAGED MODELS
!
! Name    T   L I/O  Description                              Units
! ----    --- -  -   -------------------                      ----
! CHZ     R   1  L   Chezy coefficient                         [sqrt(m)/s]
! DEPTH   R   1  I   Water depth                                       [m]
! TOTDEP  R   1  I   Total water depth                                 [m]
! G       R   1  I   Acceleration of gravity                        [m/s2]
! RHOW    R   1  I   Density of water                              [kg/m3]
! TAUWIN  R   1  O   Shearstress by wind                  [kg/m/s2 = N/m2]
! TAUFLO  R   1  O   Shearstress by flow                            [N/m2]
! TAUSCH  R   1  O   Shearstress by ships and human activity        [N/m2]
! TAU     R   1  O   Total shearstress                              [N/m2]
! TAUVEL  R   1  O   Calculated velocity based on TAU                [m/s]
! VELOC   R   1  I   Velocity                                        [m/s]

!     Logical Units : -

!     Modules called : -

!     Name     Type   Library
!     ------   -----  ------------

      IMPLICIT REAL (A-H,J-Z)

      REAL     PMSA  ( * ) , FL    (*)
      INTEGER  IPOINT( * ) , INCREM(*) , NOSEG , NOFLUX,
     +         IEXPNT(4,*) , IKNMRK(*) , NOQ1, NOQ2, NOQ3, NOQ4
!
!     Local declarations, constants in source
!
      PARAMETER ( G      =     9.8    ,
     +            RHOW   =  1000.0    ,
     +            PI     = 3.14159265 ,
     +            KARMAN = 0.41       ,
     +            GRVITY = 9.811      )
      INTEGER     LUNREP
!
      IP1  = IPOINT( 1)
      IP2  = IPOINT( 2)
      IP3  = IPOINT( 3)
      IP4  = IPOINT( 4)
      IP5  = IPOINT( 5)
      IP6  = IPOINT( 6)
      IP7  = IPOINT( 7)
      IP8  = IPOINT( 8)
      IP9  = IPOINT( 9)
      IP10 = IPOINT(10)
      IP11 = IPOINT(11)
      IP12 = IPOINT(12)
      IP13 = IPOINT(13)
      IP14 = IPOINT(14)
      IP15 = IPOINT(15)
!
      KARMC1 = SQRT(GRVITY) / KARMAN
      KARMC2 = KARMAN / SQRT(GRVITY)
!
      IFLUX = 0
      DO 9000 ISEG = 1 , NOSEG
!!    CALL DHKMRK(1,IKNMRK(ISEG),IKMRK1)
!!    IF (IKMRK1.EQ.1) THEN
      IF (BTEST(IKNMRK(ISEG),0)) THEN
      CALL DHKMRK(2,IKNMRK(ISEG),IKMRK2)
      IF (IKMRK2.EQ.0 .OR. IKMRK2 .EQ. 3) THEN
!
      H       = PMSA(IP1 )
      RL      = PMSA(IP2 )
      T       = PMSA(IP3 )
      TAUSCH  = PMSA(IP4 )
      ISWTAUVELOC = NINT(PMSA(IP5 ))
      TAUFLO  = PMSA(IP6 )
      VELOC   = PMSA(IP7 )
      CHZ     = PMSA(IP8 )
      TOTDEP  = PMSA(IP9 )
      ISWTAU  = NINT(PMSA(IP10))
      DEPTH   = PMSA(IP11)
!
!     Nelson criteria
!
      IF (H .GT. 0.55*TOTDEP) H = 0.55*TOTDEP
!
      CHZ=MAX(1.0,CHZ)
      IF ((ISWTAU .NE. 1) .AND. (ISWTAU .NE. 2) .AND. (ISWTAU .NE. 3))
     &     CALL ERRSYS ('Invald switch for tau (ISWTAU) in CALTAU', 1 )
      IF ((ISWTAUVELOC .NE. 1) .AND. (ISWTAUVELOC .NE. 2))
     &     CALL ERRSYS ('Invald switch for tau (ISWTAUVELOC) in CALTAU', 1 )

!     Initialisation
      TAU     = 0.0
      TAUWIN  = 0.0
      TAUVEL  = 0.0

!     Calculate roughness from Chezy
      ROUGH = 12. * TOTDEP / (10**(CHZ/18.))

!     Calculate Chezy 3D
      IF (IKMRK2.EQ.3) THEN
! we are doing a bottom cell --> it must be 3d
         CHZ3D = KARMC1 *   LOG( 1 + ( ( 0.5 * DEPTH )
     +                  / (TOTDEP * EXP(-1.*(1+(KARMC2 * CHZ))))
     +                       )     )
         CHZ3D=MAX(1.0,CHZ3D)
      ELSE
         CHZ3D = CHZ
      ENDIF

!     Shear stress by flow, calculate if wanted otherwise from input

      IF ( ISWTAUVELOC .EQ. 1 ) THEN
         TAUFLO = RHOW * G * VELOC**2 / CHZ3D**2
      ENDIF

!     Shear stress by waves

      IF ( ABS(H) .LT. 1.E-20 .OR. ABS(RL) .LT. 1.E-20 .OR.
     +     ABS(T) .LT. 1.E-20                              ) THEN
         TAUWIN = 0
      ELSE
         A6     = 2.0 * PI*TOTDEP/RL
         IF ( A6 .LE. 10.0 ) THEN
            UBG    = PI  * H/(T * SINH(A6))
            ALM    = UBG * T / ( 2*PI )
            RLF    = ALM/ROUGH
            IF (ISWTAU .EQ. 1) THEN
!           option Tamminga
                FWG    = 0.16* SQRT( ROUGH / ALM )
            ELSEIF (ISWTAU .EQ. 2) THEN
!           option Swart
                IF ( RLF .GT. (PI/2) ) THEN
                  FWG = 0.00251 * EXP ( 5.213 * RLF**(-0.19) )
                ELSE
                  FWG = 0.32
                ENDIF
            ELSEIF (ISWTAU .EQ. 3) THEN
!           option Soulsby
                FWG    = 0.237* RLF**(-0.52)
            ENDIF

            TAUWIN = FWG * 0.25 * RHOW * UBG**2
         ELSE
            TAUWIN = 0.0
         ENDIF
      ENDIF

  150 CONTINUE

!     Total shear stress moet eigenlijk via windrichting opgeteld worden,
!     gebeurt hier niet

!     Als TAUSCH = -1.0 dan wordt tau -1.0
      IF (TAUSCH .EQ. -1.0) THEN
          TAU = -1.0
          TAUVEL = -1.0
      ELSE
         TAU = TAUFLO  + TAUWIN + TAUSCH
!PBO  Re-calculate total sheerstress (TAU) to a total stream velocity
         TAUVEL = SQRT ( TAU * CHZ3D**2 / (RHOW * G) )
      ENDIF

      PMSA (IP12) = TAU
      PMSA (IP13) = TAUFLO
      PMSA (IP14) = TAUWIN
      PMSA (IP15) = TAUVEL
!
      ENDIF
      ENDIF
!
      IFLUX = IFLUX + NOFLUX
      IP1   = IP1   + INCREM (  1 )
      IP2   = IP2   + INCREM (  2 )
      IP3   = IP3   + INCREM (  3 )
      IP4   = IP4   + INCREM (  4 )
      IP5   = IP5   + INCREM (  5 )
      IP6   = IP6   + INCREM (  6 )
      IP7   = IP7   + INCREM (  7 )
      IP8   = IP8   + INCREM (  8 )
      IP9   = IP9   + INCREM (  9 )
      IP10  = IP10  + INCREM ( 10 )
      IP11  = IP11  + INCREM ( 11 )
      IP12  = IP12  + INCREM ( 12 )
      IP13  = IP13  + INCREM ( 13 )
      IP14  = IP14  + INCREM ( 14 )
      IP15  = IP15  + INCREM ( 15 )
!
 9000 CONTINUE
!

      RETURN
      END

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


text:
great_rsd(k1,l1).
great_rsd(a1,t1).
great_rsd(f1,a1).
great_rsd(a1,h1).
great_rsd(p1,x1).
great_rsd(cc1,kk1).
great_rsd(f1,kk1).
great_rsd(e1,dd1).
great_rsd(jj1,dd1).
great_rsd(p1,kk1).
great_rsd(t1,h1).
great_rsd(v1,hh1).
great_rsd(c1,x1).
great_rsd(c1,jj1).
great_rsd(v1,g1).
great_rsd(p1,w1).
great_rsd(u1,g1).
great_rsd(u1,ee1).
great_rsd(e1,p1).
great_rsd(f1,ee1).
great_rsd(v1,aa1).
great_rsd(k1,dd1).
great_rsd(cc1,z1).
great_rsd(q1,ff1).
great_rsd(e1,k1).

max_stars_repo_name: manoelfranca/cilppp

max_issues_repo_name: manoelfranca/cilppp

max_forks_repo_name: manoelfranca/cilppp


text:
!###########################################################
!###########################################################
!###########################################################
!###########################################################
subroutine cmpdt(uu,gg,dx,dt,ncell)
  use amr_parameters
  use hydro_parameters
  use const
  use rt_cooling_module,only:rt_isoPress,iRTisoPressVar                     !RTPress
  use rt_parameters,only:rt_advect                                          !RTpress
  implicit none
  integer::ncell
  real(dp)::dx,dt
  real(dp),dimension(1:nvector,1:nvar)::uu
  real(dp),dimension(1:nvector,1:ndim)::gg

  real(dp)::dtcell,smallp
  integer::k,idim

  smallp = smallc**2/gamma

  ! Convert to primitive variables
  do k = 1,ncell
     uu(k,1)=max(uu(k,1),smallr)
  end do
  ! Velocity
  do idim = 1,ndim
     do k = 1, ncell
        uu(k,idim+1) = uu(k,idim+1)/uu(k,1)
     end do
  end do
  ! Internal energy
  do idim = 1,ndim
     do k = 1, ncell
        uu(k,ndim+2) = uu(k,ndim+2)-half*uu(k,1)*uu(k,idim+1)**2
     end do
  end do

  ! Debug
  if(debug)then
     do k = 1, ncell
        if(uu(k,ndim+2).le.0.or.uu(k,1).le.smallr)then
           write(*,*)'stop in cmpdt'
           write(*,*)'dx   =',dx
           write(*,*)'ncell=',ncell
           write(*,*)'rho  =',uu(k,1)
           write(*,*)'P    =',uu(k,ndim+2)
           write(*,*)'vel  =',uu(k,2:ndim+1)
           stop
        end if
     end do
  end if

  ! Compute pressure
  do k = 1, ncell
     uu(k,ndim+2) = max((gamma-one)*uu(k,ndim+2),uu(k,1)*smallp)
  end do

  ! Isotropic pressure from isotropic radiation:
  if(rt_isoPress .and. rt_advect) then                                !RTpress
     do k=1,ncell                                                     !RTpress
        uu(k,ndim+2) = uu(k,ndim+2) + uu(k,iRTisoPressVar) / uu(k,1)  !RTpress
     end do                                                           !RTpress
  endif                                                               !RTpress

  ! Compute sound speed
  do k = 1, ncell
     uu(k,ndim+2)=sqrt(gamma*uu(k,ndim+2)/uu(k,1))
  end do

  ! Compute wave speed
  do k = 1, ncell
     uu(k,ndim+2)=dble(ndim)*uu(k,ndim+2)
  end do
  do idim = 1,ndim
     do k = 1, ncell
        uu(k,ndim+2)=uu(k,ndim+2)+abs(uu(k,idim+1))
     end do
  end do

  ! Compute gravity strength ratio
  do k = 1, ncell
     uu(k,1)=zero
  end do
  do idim = 1,ndim
     do k = 1, ncell
        uu(k,1)=uu(k,1)+abs(gg(k,idim))
     end do
  end do
  do k = 1, ncell
     uu(k,1)=uu(k,1)*dx/uu(k,ndim+2)**2
     uu(k,1)=MAX(uu(k,1),0.0001_dp)
  end do

  ! Compute maximum time step for each authorized cell
  dt = courant_factor*dx/smallc
  do k = 1,ncell
     dtcell=dx/uu(k,ndim+2)*(sqrt(one+two*courant_factor*uu(k,1))-one)/uu(k,1)
     dt = min(dt,dtcell)
  end do

end subroutine cmpdt
!###########################################################
!###########################################################
!###########################################################
!###########################################################
subroutine hydro_refine(ug,um,ud,ok,nn)
  use amr_parameters
  use hydro_parameters
  use const
#ifdef RT
  use rt_parameters
#endif
  implicit none
  ! dummy arguments
  integer nn
  real(dp)::ug(1:nvector,1:nvar)
  real(dp)::um(1:nvector,1:nvar)
  real(dp)::ud(1:nvector,1:nvar)
  logical ::ok(1:nvector)

  integer::k,idim
  real(dp),dimension(1:nvector),save::eking,ekinm,ekind
  real(dp)::dg,dm,dd,pg,pm,pd,vg,vm,vd,cg,cm,cd,error

  ! Convert to primitive variables
  do k = 1,nn
     ug(k,1) = max(ug(k,1),smallr)
     um(k,1) = max(um(k,1),smallr)
     ud(k,1) = max(ud(k,1),smallr)
  end do
  ! Velocity
  do idim = 1,ndim
     do k = 1,nn
        ug(k,idim+1) = ug(k,idim+1)/ug(k,1)
        um(k,idim+1) = um(k,idim+1)/um(k,1)
        ud(k,idim+1) = ud(k,idim+1)/ud(k,1)
     end do
  end do
  ! Pressure
  do k = 1,nn
     eking(k) = zero
     ekinm(k) = zero
     ekind(k) = zero
  end do
  do idim = 1,ndim
     do k = 1,nn
        eking(k) = eking(k) + half*ug(k,1)*ug(k,idim+1)**2
        ekinm(k) = ekinm(k) + half*um(k,1)*um(k,idim+1)**2
        ekind(k) = ekind(k) + half*ud(k,1)*ud(k,idim+1)**2
     end do
  end do
  do k = 1,nn
     ug(k,ndim+2) = (gamma-one)*(ug(k,ndim+2)-eking(k))
     um(k,ndim+2) = (gamma-one)*(um(k,ndim+2)-ekinm(k))
     ud(k,ndim+2) = (gamma-one)*(ud(k,ndim+2)-ekind(k))
  end do
  ! Passive scalars
#if NVAR > NDIM + 2
  do idim = ndim+3,nvar
     do k = 1,nn
        ug(k,idim) = ug(k,idim)/ug(k,1)
        um(k,idim) = um(k,idim)/um(k,1)
        ud(k,idim) = ud(k,idim)/ud(k,1)
     end do
  end do
#endif

  ! Compute errors
  if(err_grad_d >= 0.)then
     do k=1,nn
        dg=ug(k,1); dm=um(k,1); dd=ud(k,1)
        error=2.0d0*MAX( &
             & ABS((dd-dm)/(dd+dm+floor_d)) , &
             & ABS((dm-dg)/(dm+dg+floor_d)) )
        ok(k) = ok(k) .or. error > err_grad_d
     end do
  end if

  if(err_grad_p >= 0.)then
     do k=1,nn
        pg=ug(k,ndim+2); pm=um(k,ndim+2); pd=ud(k,ndim+2)
        error=2.0d0*MAX( &
             & ABS((pd-pm)/(pd+pm+floor_p)), &
             & ABS((pm-pg)/(pm+pg+floor_p)) )
        ok(k) = ok(k) .or. error > err_grad_p
     end do
  end if

  if(err_grad_u >= 0.)then
     do idim = 1,ndim
        do k=1,nn
           vg=ug(k,idim+1); vm=um(k,idim+1); vd=ud(k,idim+1)
           cg=sqrt(max(gamma*ug(k,ndim+2)/ug(k,1),floor_u**2))
           cm=sqrt(max(gamma*um(k,ndim+2)/um(k,1),floor_u**2))
           cd=sqrt(max(gamma*ud(k,ndim+2)/ud(k,1),floor_u**2))
           error=2.0d0*MAX( &
                & ABS((vd-vm)/(cd+cm+ABS(vd)+ABS(vm)+floor_u)) , &
                & ABS((vm-vg)/(cm+cg+ABS(vm)+ABS(vg)+floor_u)) )
           ok(k) = ok(k) .or. error > err_grad_u
        end do
     end do
  end if

#ifdef RT
  ! Ionization state (only Hydrogen)
  if(rt_err_grad_xHII >= 0.) then !---------------------------------------
     do k=1,nn
        dg=min(1.d0,max(0.d0,ug(k,iIons)))
        dm=min(1.d0,max(0.d0,um(k,iIons)))
        dd=min(1.d0,max(0.d0,ud(k,iIons)))
        error=2.0d0*MAX( &
             & ABS((dd-dm)/(dd+dm+rt_floor_xHII)) , &
             & ABS((dm-dg)/(dm+dg+rt_floor_xHII)) )
        ok(k) = ok(k) .or. error > rt_err_grad_xHII
     end do
  end if

  ! Neutral state (only Hydrogen)
  if(rt_err_grad_xHI  >= 0.) then !---------------------------------------
     do k=1,nn
        dg=min(1.d0,max(0.d0,1.d0 - ug(k,iIons)))
        dm=min(1.d0,max(0.d0,1.d0 - um(k,iIons)))
        dd=min(1.d0,max(0.d0,1.d0 - ud(k,iIons)))
        error=2.0d0*MAX( &
             & ABS((dd-dm)/(dd+dm+rt_floor_xHI)) , &
             & ABS((dm-dg)/(dm+dg+rt_floor_xHI)) )
        ok(k) = ok(k) .or. error > rt_err_grad_xHI
     end do
  end if
#endif

end subroutine hydro_refine
!###########################################################
!###########################################################
!###########################################################
!###########################################################
subroutine riemann_approx(qleft,qright,qgdnv,fgdnv,ngrid)
  use amr_parameters
  use hydro_parameters
  use const
  implicit none

  ! dummy arguments
  integer::ngrid
  real(dp),dimension(1:nvector,1:nvar)::qleft,qright,qgdnv,fgdnv

  ! local arrays
  real(dp),dimension(1:nvector),save::rl   ,ul   ,pl   ,cl
  real(dp),dimension(1:nvector),save::rr   ,ur   ,pr   ,cr
  real(dp),dimension(1:nvector),save::ro   ,uo   ,po   ,co
  real(dp),dimension(1:nvector),save::rstar,ustar,pstar,cstar
  real(dp),dimension(1:nvector),save::wl   ,wr   ,wo
  real(dp),dimension(1:nvector),save::sgnm ,spin ,spout,ushock
  real(dp),dimension(1:nvector),save::frac ,delp ,pold
  integer ,dimension(1:nvector),save::ind  ,ind2

  ! local variables
  real(dp)::smallp, gamma6, ql, qr, usr, usl, wwl, wwr, smallpp, entho, etot
  integer ::i, j, n, iter, n_new

  ! constants
  smallp = smallc**2/gamma
  smallpp = smallr*smallp
  gamma6 = (gamma+one)/(two*gamma)
  entho = one/(gamma-one)

  ! Pressure, density and velocity
  do i=1,ngrid
     rl(i)=MAX(qleft (i,1),smallr)
     ul(i)=    qleft (i,2)
     pl(i)=MAX(qleft (i,3),rl(i)*smallp)
     rr(i)=MAX(qright(i,1),smallr)
     ur(i)=    qright(i,2)
     pr(i)=MAX(qright(i,3),rr(i)*smallp)
  end do

  ! Lagrangian sound speed
  do i=1,ngrid
     cl(i) = gamma*pl(i)*rl(i)
     cr(i) = gamma*pr(i)*rr(i)
  end do

  ! First guess
  do i=1,ngrid
     wl(i)= sqrt(cl(i)); wr(i) = sqrt(cr(i))
     pstar(i) = ((wr(i)*pl(i)+wl(i)*pr(i))+wl(i)*wr(i)*(ul(i)-ur(i)))/(wl(i)+wr(i))
     pstar(i) = MAX(pstar(i),0.0_dp)
     pold (i)= pstar(i)
  end do
  n = ngrid
  do i = 1,n
     ind(i)=i
  end do

  ! Newton-Raphson iterations to find pstar at the required accuracy
  ! for a two-shock Riemann problem
  do iter = 1,niter_riemann
     do i=1,n
        wwl=sqrt(cl(ind(i))*(one+gamma6*(pold(i)-pl(ind(i)))/pl(ind(i))))
        wwr=sqrt(cr(ind(i))*(one+gamma6*(pold(i)-pr(ind(i)))/pr(ind(i))))
        ql=two*wwl**3/(wwl**2+cl(ind(i)))
        qr=two*wwr**3/(wwr**2+cr(ind(i)))
        usl=ul(ind(i))-(pold(i)-pl(ind(i)))/wwl
        usr=ur(ind(i))+(pold(i)-pr(ind(i)))/wwr
        delp(i)=MAX(qr*ql/(qr+ql)*(usl-usr),-pold(i))
     end do
     do i=1,n
        pold(i)=pold(i)+delp(i)
     end do
     ! Convergence indicator
     do i=1,n
        uo(i)=ABS(delp(i)/(pold(i)+smallpp))
     end do
     n_new=0
     do i=1,n
        if(uo(i)>1.d-06)then
           n_new=n_new+1
           ind2(n_new)=ind (i)
           po  (n_new)=pold(i)
        end if
     end do
     j=n_new
     do i=1,n
        if(uo(i)<=1.d-06)then
           n_new=n_new+1
           ind2(n_new)=ind (i)
           po  (n_new)=pold(i)
        end if
     end do
     ind (1:n)=ind2(1:n)
     pold(1:n)=po  (1:n)
     n=j
  end do

  ! Star region pressure
  ! for a two-shock Riemann problem
  do i=1,ngrid
     pstar(ind(i))=pold(i)
  end do
  do i=1,ngrid
     wl(i) = sqrt(cl(i)*(one+gamma6*(pstar(i)-pl(i))/pl(i)))
     wr(i) = sqrt(cr(i)*(one+gamma6*(pstar(i)-pr(i))/pr(i)))
  end do

  ! Star region velocity
  ! for a two shock Riemann problem
  do i=1,ngrid
     ustar(i) = half*(ul(i) + (pl(i)-pstar(i))/wl(i) + &
          &           ur(i) - (pr(i)-pstar(i))/wr(i) )
  end do

  ! Left going or right going contact wave
  do i=1,ngrid
     sgnm(i) = sign(one,ustar(i))
  end do

  ! Left or right unperturbed state
  do i=1,ngrid
     if(sgnm(i)==one)then
        ro(i) = rl(i)
        uo(i) = ul(i)
        po(i) = pl(i)
        wo(i) = wl(i)
     else
        ro(i) = rr(i)
        uo(i) = ur(i)
        po(i) = pr(i)
        wo(i) = wr(i)
     end if
  end do
  do i=1,ngrid
     co(i) = max(smallc,sqrt(abs(gamma*po(i)/ro(i))))
  end do

  ! Star region density
  do i=1,ngrid
     if(pstar(i)>= po(i))then
        ! Shock
        rstar(i) = ro(i)/(one+ro(i)*(po(i)-pstar(i))/wo(i)**2)
     else
        ! Rarefaction
        rstar(i) = ro(i)*(pstar(i)/po(i))**(one/gamma)
     end if
  end do

  do i=1,ngrid
     ! Prevent vacuum formation in star region
     rstar(i) = max(rstar(i),smallr)
     ! Star region sound speed
     cstar(i) = sqrt(abs(gamma*pstar(i)/rstar(i)))
     cstar(i) = max(cstar(i),smallc)
     ! Compute rarefaction head and tail speed
     spout(i) = co   (i)-sgnm(i)*uo   (i)
     spin (i) = cstar(i)-sgnm(i)*ustar(i)
     ! Compute shock speed
     ushock(i) = wo(i)/ro(i)-sgnm(i)*uo(i)
  end do

  do i=1,ngrid
     if(pstar(i)>=po(i))then
        spout(i)=ushock(i)
        spin (i)=spout (i)
     end if
  end do

  ! Sample the solution at x/t=0
  do i=1,ngrid
     if(spout(i)<=zero)then
        qgdnv(i,1) = ro(i)
        qgdnv(i,2) = uo(i)
        qgdnv(i,3) = po(i)
     else if(spin(i)>=zero)then
        qgdnv(i,1) = rstar(i)
        qgdnv(i,2) = ustar(i)
        qgdnv(i,3) = pstar(i)
     else
        frac(i)=spout(i)/(spout(i)-spin(i))
        qgdnv(i,2) = frac(i)*ustar(i) + (one - frac(i))*uo(i)
        qgdnv(i,3) = frac(i)*pstar(i) + (one - frac(i))*po(i)
        qgdnv(i,1) = ro(i)*(qgdnv(i,3)/po(i))**(one/gamma)
     end if
  end do

  ! Passive scalars
#if NVAR > 3
  do n = 4,nvar
     do i=1,ngrid
        if(sgnm(i)==one)then
           qgdnv(i,n) = qleft(i,n)
        else
           qgdnv(i,n) = qright(i,n)
        end if
     end do
  end do
#endif

  ! Compute fluxes
  do i = 1, ngrid
     fgdnv(i,1) = qgdnv(i,1)*qgdnv(i,2)  ! Mass density
     fgdnv(i,2) = qgdnv(i,3)+qgdnv(i,1)*qgdnv(i,2)**2  ! Normal momentum
     etot = qgdnv(i,3)*entho + half*qgdnv(i,1)*qgdnv(i,2)**2
#if NDIM>1
     etot = etot + half*qgdnv(i,1)*qgdnv(i,4)**2
#endif
#if NDIM>2
     etot = etot + half*qgdnv(i,1)*qgdnv(i,5)**2
#endif
     fgdnv(i,3) = qgdnv(i,2)*(etot+qgdnv(i,3))     ! Total energy
  end do
  ! Other advected quantities
#if NVAR > 3
  do n = 4, nvar
     do i = 1, ngrid
        fgdnv(i,n) = fgdnv(i,1)*qgdnv(i,n)
     end do
  end do
#endif

end subroutine riemann_approx
!###########################################################
!###########################################################
!###########################################################
!###########################################################
subroutine riemann_acoustic(qleft,qright,qgdnv,fgdnv,ngrid)
  use amr_parameters
  use hydro_parameters
  use const
  implicit none

  ! dummy arguments
  integer::ngrid
  real(dp),dimension(1:nvector,1:nvar)::qleft,qright,qgdnv,fgdnv

  ! local variables
  integer::i,n
  real(dp)::smallp, entho, etot

  ! local arrays
  real(dp),dimension(1:nvector),save::rl   ,ul   ,pl   ,cl
  real(dp),dimension(1:nvector),save::rr   ,ur   ,pr   ,cr
  real(dp),dimension(1:nvector),save::ro   ,uo   ,po   ,co
  real(dp),dimension(1:nvector),save::rstar,ustar,pstar,cstar
  real(dp),dimension(1:nvector),save::wl   ,wr   ,wo
  real(dp),dimension(1:nvector),save::sgnm ,spin ,spout,ushock
  real(dp),dimension(1:nvector),save::frac

  ! constants
  smallp = smallc**2/gamma
  entho = one/(gamma-one)

  ! Initial states pressure, density and velocity
  do i=1,ngrid
     rl(i)=max(qleft (i,1),smallr)
     ul(i)=    qleft (i,2)
     pl(i)=max(qleft (i,3),rl(i)*smallp)
     rr(i)=max(qright(i,1),smallr)
     ur(i)=    qright(i,2)
     pr(i)=max(qright(i,3),rr(i)*smallp)
  end do

  ! Acoustic star state
  do i=1,ngrid
     cl(i) = sqrt(gamma*pl(i)/rl(i))
     cr(i) = sqrt(gamma*pr(i)/rr(i))
     wl(i) = cl(i)*rl(i)
     wr(i) = cr(i)*rr(i)
     pstar(i) = ((wr(i)*pl(i)+wl(i)*pr(i))+wl(i)*wr(i)*(ul(i)-ur(i))) &
          &   / (wl(i)+wr(i))
     ustar(i) = ((wr(i)*ur(i)+wl(i)*ul(i))+(pl(i)-pr(i))) &
          &   / (wl(i)+wr(i))
!!$     pstar(i) = MAX(pstar(i),zero)
  end do

  ! Left going or right going contact wave
  do i=1,ngrid
     sgnm(i) = sign(one,ustar(i))
  end do

  ! Left or right unperturbed state
  do i=1,ngrid
     if(sgnm(i)==one)then
        ro(i) = rl(i)
        uo(i) = ul(i)
        po(i) = pl(i)
        wo(i) = wl(i)
        co(i) = cl(i)
     else
        ro(i) = rr(i)
        uo(i) = ur(i)
        po(i) = pr(i)
        wo(i) = wr(i)
        co(i) = cr(i)
     end if
  end do

  ! Star region density and sound speed
  do i=1,ngrid
     rstar(i) = ro(i)+(pstar(i)-po(i))/co(i)**2
     rstar(i) = max(rstar(i),smallr)
     cstar(i) = sqrt(abs(gamma*pstar(i)/rstar(i)))
     cstar(i) = max(cstar(i),smallc)
  end do

  ! Head and tail speed of rarefaction
  do i=1,ngrid
     spout(i) = co   (i)-sgnm(i)*uo   (i)
     spin (i) = cstar(i)-sgnm(i)*ustar(i)
  end do

  ! Shock speed
  do i=1,ngrid
     ushock(i) = half*(spin(i)+spout(i))
     ushock(i) = max(ushock(i),-sgnm(i)*ustar(i))
  end do
  do i=1,ngrid
     if(pstar(i)>=po(i))then
        spout(i)=ushock(i)
        spin (i)=spout (i)
     end if
  end do

  ! Sample the solution at x/t=0
  do i=1,ngrid
     if(spout(i)<zero)then      ! Initial state
        qgdnv(i,1) = ro(i)
        qgdnv(i,2) = uo(i)
        qgdnv(i,3) = po(i)
     else if(spin(i)>=zero)then  ! Star region
        qgdnv(i,1) = rstar(i)
        qgdnv(i,2) = ustar(i)
        qgdnv(i,3) = pstar(i)
     else                        ! Rarefaction
        frac(i) = spout(i)/(spout(i)-spin(i))
        qgdnv(i,1) = frac(i)*rstar(i) + (one - frac(i))*ro(i)
        qgdnv(i,2) = frac(i)*ustar(i) + (one - frac(i))*uo(i)
        qgdnv(i,3) = frac(i)*pstar(i) + (one - frac(i))*po(i)
     end if
  end do

  ! Passive scalars
#if NVAR > 3
  do n = 4,nvar
     do i=1,ngrid
        if(sgnm(i)==one)then
           qgdnv(i,n) = qleft (i,n)
        else
           qgdnv(i,n) = qright(i,n)
        end if
     end do
  end do
#endif

  ! Compute fluxes
  do i = 1, ngrid
     fgdnv(i,1) = qgdnv(i,1)*qgdnv(i,2)  ! Mass density
     fgdnv(i,2) = qgdnv(i,3)+qgdnv(i,1)*qgdnv(i,2)**2  ! Normal momentum
     etot = qgdnv(i,3)*entho + half*qgdnv(i,1)*qgdnv(i,2)**2
#if NDIM>1
     etot = etot             + half*qgdnv(i,1)*qgdnv(i,4)**2
#endif
#if NDIM>2
     etot = etot             + half*qgdnv(i,1)*qgdnv(i,5)**2
#endif
     fgdnv(i,3) = qgdnv(i,2)*(etot+qgdnv(i,3))     ! Total energy
  end do
  ! Other advected quantities
#if NVAR > 3
  do n = 4, nvar
     do i = 1, ngrid
        fgdnv(i,n) = fgdnv(i,1)*qgdnv(i,n)
     end do
  end do
#endif

end subroutine riemann_acoustic
!###########################################################
!###########################################################
!###########################################################
!###########################################################
subroutine riemann_llf(qleft,qright,qgdnv,fgdnv,ngrid)
  use amr_parameters
  use hydro_parameters
  use const
  implicit none

  ! dummy arguments
  integer::ngrid
  real(dp),dimension(1:nvector,1:nvar)::qleft,qright,qgdnv,fgdnv

  ! local arrays
  real(dp),dimension(1:nvector,1:nvar),save::fleft,fright,uleft,uright
  real(dp),dimension(1:nvector),save::cmax

  ! local variables
  integer::i,n
  real(dp)::smallp, entho
  real(dp)::rl   ,ul   ,pl   ,cl
  real(dp)::rr   ,ur   ,pr   ,cr

  ! constants
  smallp = smallc**2/gamma
  entho = one/(gamma-one)

  ! Maximum wave speed
  do i=1,ngrid
     rl=max(qleft (i,1),smallr)
     ul=    qleft (i,2)
     pl=max(qleft (i,3),rl*smallp)
     rr=max(qright(i,1),smallr)
     ur=    qright(i,2)
     pr=max(qright(i,3),rr*smallp)
     cl= sqrt(gamma*pl/rl)
     cr= sqrt(gamma*pr/rr)
     cmax(i)=max(abs(ul)+cl,abs(ur)+cr)
  end do

  ! Compute average velocity
  do i=1,ngrid
     qgdnv(i,2) = half*(qleft(i,2)+qright(i,2))
  end do

  ! Compute conservative variables
  do i = 1, ngrid
     ! Mass density
     uleft (i,1) = qleft (i,1)
     uright(i,1) = qright(i,1)
     ! Normal momentum
     uleft (i,2) = qleft (i,1)*qleft (i,2)
     uright(i,2) = qright(i,1)*qright(i,2)
     ! Total energy
     uleft (i,3) = qleft (i,3)*entho + half*qleft (i,1)*qleft (i,2)**2
     uright(i,3) = qright(i,3)*entho + half*qright(i,1)*qright(i,2)**2
#if NDIM>1
     uleft (i,3) = uleft (i,3)       + half*qleft (i,1)*qleft (i,4)**2
     uright(i,3) = uright(i,3)       + half*qright(i,1)*qright(i,4)**2
#endif
#if NDIM>2
     uleft (i,3) = uleft (i,3)       + half*qleft (i,1)*qleft (i,5)**2
     uright(i,3) = uright(i,3)       + half*qright(i,1)*qright(i,5)**2
#endif
  end do
  ! Other advected quantities
#if NVAR > 3
  do n = 4, nvar
     do i = 1, ngrid
        uleft (i,n) = qleft (i,1)*qleft (i,n)
        uright(i,n) = qright(i,1)*qright(i,n)
     end do
  end do
#endif

  ! Compute left and right fluxes
  do i = 1, ngrid
     ! Mass density
     fleft (i,1) = uleft (i,2)
     fright(i,1) = uright(i,2)
     ! Normal momentum
     fleft (i,2) = qleft (i,3)+uleft (i,2)*qleft (i,2)
     fright(i,2) = qright(i,3)+uright(i,2)*qright(i,2)
     ! Total energy
     fleft (i,3) = qleft (i,2)*(uleft (i,3)+qleft (i,3))
     fright(i,3) = qright(i,2)*(uright(i,3)+qright(i,3))
  end do
  ! Other advected quantities
#if NVAR > 3
  do n = 4, nvar
     do i = 1, ngrid
        fleft (i,n) = fleft (i,1)*qleft (i,n)
        fright(i,n) = fright(i,1)*qright(i,n)
     end do
  end do
#endif

  ! Compute Lax-Friedrich fluxes
  do n = 1, nvar
     do i = 1, ngrid
        fgdnv(i,n) = half*(fleft(i,n)+fright(i,n)-cmax(i)*(uright(i,n)-uleft(i,n)))
     end do
  end do

end subroutine riemann_llf
!###########################################################
!###########################################################
!###########################################################
!###########################################################
!###########################################################
!###########################################################
!###########################################################
!###########################################################
subroutine riemann_hllc(qleft,qright,qgdnv,fgdnv,ngrid)
  use amr_parameters
  use hydro_parameters
  use const
  use rt_cooling_module,only:rt_isoPress,iRTisoPressVar                     !RTPress
  use rt_parameters,only:rt_advect                                          !RTpress
  implicit none

  ! HLLC Riemann solver (Toro)
  integer::ngrid
  real(dp),dimension(1:nvector,1:nvar)::qleft,qright,qgdnv,fgdnv

  REAL(dp)::SL,SR
  REAL(dp)::entho
  REAL(dp)::rl,pl,ul,ecinl,etotl,ptotl
  REAL(dp)::rr,pr,ur,ecinr,etotr,ptotr
  REAL(dp)::cfastl,rcl,rstarl
  REAL(dp)::cfastr,rcr,rstarr
  REAL(dp)::etotstarl,etotstarr
  REAL(dp)::ustar,ptotstar
  REAL(dp)::ro,uo,ptoto,etoto
  REAL(dp)::smallp

  INTEGER ::ivar, i

  ! constants
  smallp = smallc**2/gamma
  entho = one/(gamma-one)

  do i=1,ngrid

     ! Left variables
     rl=max(qleft (i,1),smallr)
     Pl=max(qleft (i,3),rl*smallp)
     ul=    qleft (i,2)

     ecinl = half*rl*ul*ul
#if NDIM>1
     ecinl=ecinl+half*rl*qleft(i,4)**2
#endif
#if NDIM>2
     ecinl=ecinl+half*rl*qleft(i,5)**2
#endif
     etotl = Pl*entho+ecinl
     ! Here update Ptot but not etot, according to O. Agertz !joki!!!!      !RTpress
     ! Add isotropic radiation pressure                                     !RTpress
     if(rt_isoPress .and. rt_advect) then                                   !RTpress
        Pl = Pl + MAX(qleft(i,iRTisoPressVar),0d0)                          !RTpress
     endif                                                                  !RTpress
     Ptotl = Pl

     ! Right variables
     rr=max(qright(i,1),smallr)
     Pr=max(qright(i,3),rr*smallp)
     ur=    qright(i,2)

     ecinr = half*rr*ur*ur
#if NDIM>1
     ecinr=ecinr+half*rr*qright(i,4)**2
#endif
#if NDIM>2
     ecinr=ecinr+half*rr*qright(i,5)**2
#endif
     etotr = Pr*entho+ecinr
     ! Add isotropic radiation pressure                                     !RTpress
     if(rt_isoPress .and. rt_advect) then                                   !RTpress
        Pr = Pr + MAX(qright(i,iRTisoPressVar),0d0)                         !RTpress
     endif                                                                  !RTpress
     Ptotr = Pr

     ! Find the largest eigenvalues in the normal direction to the interface
     cfastl=sqrt(max(gamma*Pl/rl,smallc**2))
     cfastr=sqrt(max(gamma*Pr/rr,smallc**2))

     ! Compute HLL wave speed
     SL=min(ul,ur)-max(cfastl,cfastr)
     SR=max(ul,ur)+max(cfastl,cfastr)

     ! Compute lagrangian sound speed
     rcl=rl*(ul-SL)
     rcr=rr*(SR-ur)

     ! Compute acoustic star state
     ustar   =(rcr*ur   +rcl*ul   +  (Ptotl-Ptotr))/(rcr+rcl)
     Ptotstar=(rcr*Ptotl+rcl*Ptotr+rcl*rcr*(ul-ur))/(rcr+rcl)

     ! Left star region variables
     rstarl=rl*(SL-ul)/(SL-ustar)
     etotstarl=((SL-ul)*etotl-Ptotl*ul+Ptotstar*ustar)/(SL-ustar)

     ! Right star region variables
     rstarr=rr*(SR-ur)/(SR-ustar)
     etotstarr=((SR-ur)*etotr-Ptotr*ur+Ptotstar*ustar)/(SR-ustar)

     ! Sample the solution at x/t=0
     if(SL>0d0)then
        ro=rl
        uo=ul
        Ptoto=Ptotl
        etoto=etotl
     else if(ustar>0d0)then
        ro=rstarl
        uo=ustar
        Ptoto=Ptotstar
        etoto=etotstarl
     else if (SR>0d0)then
        ro=rstarr
        uo=ustar
        Ptoto=Ptotstar
        etoto=etotstarr
     else
        ro=rr
        uo=ur
        Ptoto=Ptotr
        etoto=etotr
     end if

     ! Compute the Godunov flux
     fgdnv(i,1) = ro*uo
     fgdnv(i,2) = ro*uo*uo+Ptoto
     fgdnv(i,3) = (etoto+Ptoto)*uo
#if NVAR > 3
     do ivar = 4,nvar
        if(fgdnv(i,1)>0)then
           fgdnv(i,ivar) = fgdnv(i,1)*qleft (i,ivar)
        else
           fgdnv(i,ivar) = fgdnv(i,1)*qright(i,ivar)
        endif
     end do
#endif
     ! Compute the Godunov velocity
     qgdnv(i,2) = uo

  end do

end subroutine riemann_hllc
!###########################################################
!###########################################################
!###########################################################
!###########################################################
subroutine riemann_hll(qleft,qright,qgdnv,fgdnv,ngrid)
  USE amr_parameters
  USE const
  USE hydro_parameters
  ! 1D HLL Riemann solver
  IMPLICIT NONE
  integer::ngrid
  real(dp),dimension(1:nvector,1:nvar)::qleft,qright,qgdnv,fgdnv
  real(dp),dimension(1:nvector,1:nvar),save::fleft,fright,uleft,uright
  real(dp),dimension(1:nvector),save::SL,SR
  integer::i,n
  real(dp)::smallp, entho
  real(dp)::rl   ,ul   ,pl   ,cl
  real(dp)::rr   ,ur   ,pr   ,cr

  ! constants
  smallp = smallc**2/gamma
  entho = one/(gamma-one)

  ! Maximum wave speed
  do i=1,ngrid
     rl=max(qleft (i,1),smallr)
     ul=    qleft (i,2)
     pl=max(qleft (i,3),rl*smallp)
     rr=max(qright(i,1),smallr)
     ur=    qright(i,2)
     pr=max(qright(i,3),rr*smallp)
     cl= sqrt(gamma*pl/rl)
     cr= sqrt(gamma*pr/rr)
     SL(i)=min(min(ul,ur)-max(cl,cr),zero)
     SR(i)=max(max(ul,ur)+max(cl,cr),zero)
  end do

  ! Compute average velocity
  do i=1,ngrid
     qgdnv(i,2) = half*(qleft(i,2)+qright(i,2))
  end do

  ! Compute conservative variables
  do i = 1, ngrid
     uleft (i,1) = qleft (i,1)
     uright(i,1) = qright(i,1)
     uleft (i,2) = qleft (i,1)*qleft (i,2)
     uright(i,2) = qright(i,1)*qright(i,2)
     uleft (i,3) = qleft (i,3)*entho + half*qleft (i,1)*qleft (i,2)**2
     uright(i,3) = qright(i,3)*entho + half*qright(i,1)*qright(i,2)**2
#if NDIM>1
     uleft (i,3) = uleft (i,3)       + half*qleft (i,1)*qleft (i,4)**2
     uright(i,3) = uright(i,3)       + half*qright(i,1)*qright(i,4)**2
#endif
#if NDIM>2
     uleft (i,3) = uleft (i,3)       + half*qleft (i,1)*qleft (i,5)**2
     uright(i,3) = uright(i,3)       + half*qright(i,1)*qright(i,5)**2
#endif
  end do
  ! Other advected quantities
#if NVAR > 3
  do n = 4, nvar
     do i = 1, ngrid
        uleft (i,n) = qleft (i,1)*qleft (i,n)
        uright(i,n) = qright(i,1)*qright(i,n)
     end do
  end do
#endif

  ! Compute left and right fluxes
  do i = 1, ngrid
     fleft (i,1) = uleft (i,2)
     fright(i,1) = uright(i,2)
     fleft (i,2) = qleft (i,3)+uleft (i,2)*qleft (i,2)
     fright(i,2) = qright(i,3)+uright(i,2)*qright(i,2)
     fleft (i,3) = qleft (i,2)*(uleft (i,3)+qleft (i,3))
     fright(i,3) = qright(i,2)*(uright(i,3)+qright(i,3))
  end do
  ! Other advected quantities
#if NVAR > 3
  do n = 4, nvar
     do i = 1, ngrid
        fleft (i,n) = fleft (i,1)*qleft (i,n)
        fright(i,n) = fright(i,1)*qright(i,n)
     end do
  end do
#endif

  ! Compute HLL fluxes
  do n = 1, nvar
     do i = 1, ngrid
        fgdnv(i,n) = (SR(i)*fleft(i,n)-SL(i)*fright(i,n) &
             & + SR(i)*SL(i)*(uright(i,n)-uleft(i,n)))/(SR(i)-SL(i))
     end do
  end do

end subroutine riemann_hll
!###########################################################
!###########################################################
!###########################################################
!###########################################################

max_stars_repo_name: rknop/amuse

max_issues_repo_name: rknop/amuse

max_forks_repo_name: rknop/amuse


text:
! { dg-do run }
! PR 34887 - reverse tabs followed by a slash used to confuse I/O.
      program main
      character(len=2) :: b, a
      open(10,form="formatted")
      write (10,'(3X, A, T1, A,/)') 'aa', 'bb'
      rewind(10)
      read (10,'(A2,1X,A2)') b,a
      if (a /= 'aa' .or. b /= 'bb') STOP 1
      close(10,status="delete")
      end


MODULE trdmxl_oce
  USE par_oce
  IMPLICIT NONE
  PRIVATE
  PUBLIC :: trdmxl_oce_alloc
  INTEGER, PUBLIC, PARAMETER :: jpltrd = 12
  INTEGER, PUBLIC :: jpktrd
  INTEGER, PUBLIC, PARAMETER :: jpmxl_xad = 1
  INTEGER, PUBLIC, PARAMETER :: jpmxl_yad = 2
  INTEGER, PUBLIC, PARAMETER :: jpmxl_zad = 3
  INTEGER, PUBLIC, PARAMETER :: jpmxl_ldf = 4
  INTEGER, PUBLIC, PARAMETER :: jpmxl_zdf = 5
  INTEGER, PUBLIC, PARAMETER :: jpmxl_npc = 6
  INTEGER, PUBLIC, PARAMETER :: jpmxl_bbc = 7
  INTEGER, PUBLIC, PARAMETER :: jpmxl_bbl = 8
  INTEGER, PUBLIC, PARAMETER :: jpmxl_for = 9
  INTEGER, PUBLIC, PARAMETER :: jpmxl_dmp = 10
  INTEGER, PUBLIC, PARAMETER :: jpmxl_zdfp = 11
  INTEGER, PUBLIC, PARAMETER :: jpmxl_atf = 12
  INTEGER, PUBLIC :: nn_ctls = 0
  REAL(KIND = wp), PUBLIC :: rn_rho_c = 0.01
  REAL(KIND = wp), PUBLIC :: rn_ucf = 1.
  CHARACTER(LEN = 32), PUBLIC :: cn_trdrst_in = "restart_mxl"
  CHARACTER(LEN = 32), PUBLIC :: cn_trdrst_out = "restart_mxl"
  LOGICAL, PUBLIC :: ln_trdmxl_instant = .FALSE.
  LOGICAL, PUBLIC :: ln_trdmxl_restart = .FALSE.
  CHARACTER(LEN = 80), PUBLIC :: clname, ctrd(jpltrd + 1, 2)
  INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:, :) :: nmxl
  INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:, :) :: nbol
  REAL(KIND = wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:, :, :) :: wkx
  REAL(KIND = wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:, :) :: hmxl, tml, sml, tmlb, smlb, tmlbb, smlbb, tmlbn, smlbn, tmltrdm, smltrdm, tml_sum, tml_sumb, tmltrd_atf_sumb, sml_sum, sml_sumb, smltrd_atf_sumb, hmxl_sum, hmxlbn
  REAL(KIND = wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:, :) :: tmlatfb, tmlatfn, smlatfb, smlatfn, tmlatfm, smlatfm
  REAL(KIND = wp), PUBLIC, ALLOCATABLE, DIMENSION(:, :, :) :: tmltrd, smltrd, tmltrd_sum, tmltrd_csum_ln, tmltrd_csum_ub, smltrd_sum, smltrd_csum_ln, smltrd_csum_ub
  CONTAINS
  INTEGER FUNCTION trdmxl_oce_alloc()
    USE lib_mpp
    INTEGER :: ierr(5)
    jpktrd = jpk
    ierr(:) = 0
    ALLOCATE(nmxl(jpi, jpj), nbol(jpi, jpj), wkx(jpi, jpj, jpk), hmxl(jpi, jpj), tml(jpi, jpj), sml(jpi, jpj), tmlb(jpi, jpj), smlb(jpi, jpj), tmlbb(jpi, jpj), smlbb(jpi, jpj), STAT = ierr(1))
    ALLOCATE(tmlbn(jpi, jpj), smlbn(jpi, jpj), tmltrdm(jpi, jpj), smltrdm(jpi, jpj), tml_sum(jpi, jpj), tml_sumb(jpi, jpj), tmltrd_atf_sumb(jpi, jpj), STAT = ierr(2))
    ALLOCATE(sml_sum(jpi, jpj), sml_sumb(jpi, jpj), smltrd_atf_sumb(jpi, jpj), hmxl_sum(jpi, jpj), hmxlbn(jpi, jpj), tmlatfb(jpi, jpj), tmlatfn(jpi, jpj), STAT = ierr(3))
    ALLOCATE(smlatfb(jpi, jpj), smlatfn(jpi, jpj), tmlatfm(jpi, jpj), smlatfm(jpi, jpj), tmltrd(jpi, jpj, jpltrd), smltrd(jpi, jpj, jpltrd), STAT = ierr(4))
    ALLOCATE(tmltrd_sum(jpi, jpj, jpltrd), tmltrd_csum_ln(jpi, jpj, jpltrd), tmltrd_csum_ub(jpi, jpj, jpltrd), smltrd_sum(jpi, jpj, jpltrd), smltrd_csum_ln(jpi, jpj, jpltrd), smltrd_csum_ub(jpi, jpj, jpltrd), STAT = ierr(5))
    trdmxl_oce_alloc = MAXVAL(ierr)
    IF (lk_mpp) CALL mpp_sum(trdmxl_oce_alloc)
    IF (trdmxl_oce_alloc /= 0) CALL ctl_warn('trdmxl_oce_alloc: failed to allocate arrays')
  END FUNCTION trdmxl_oce_alloc
END MODULE trdmxl_oce