text:
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!    SUBROUTINES FOR BOUNDARY CONDITIONS    !!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine apply_BC()
	use constants

	implicit none
	integer :: i

	!boundary conditions for reflection at piston side
	do i = 1, nspec
		U1D(1,3*(i-1)+1) = rho(1,i)!U1D(2,3*(i-1)+1) - lm * F1D(2,3*(i-1)+1) !continuity
		U1D(1,3*(i-1)+2) = rho(1,i) * u(1,i) !momentum
		U1D(1,3*(i-1)+3) = p(1,i) / (g-1) + 0.5 * rho(1,i) * u(1,i)**2!U1D(2,3*(i-1)+3) - lm * F1D(2,3*(i-1)+3) !energy
	enddo
	!now electrons
	U1D(1,neqi+1) =  U1D(2,neqi+1) - lm * F1D(2,neqi+1)

	if (.not.(geom=="spherical")) then
		!boundary condition for reflection at wall side
		do i = 1, nspec
			U1D(nz,3*(i-1)+1) = U1D(nz-1,3*(i-1)+1) - lm * F1D(nz-1,3*(i-1)+1) !continuity
			U1D(nz,3*(i-1)+2) = 0. !momentum
			U1D(nz,3*(i-1)+3) = U1D(nz-1,3*(i-1)+3) - lm * F1D(nz-1,3*(i-1)+3) !energy
		enddo
		!now electrons
		U1D(nz,neqi+1) =  U1D(nz-1,neqi+1)
	else !spherical geometry
		!boundary condition for reflection at wall side
		do i = 1, nspec
			U1D(nz,3*(i-1)+1) = U1D(nz-1,3*(i-1)+1) - lm * F1D(nz-1,3*(i-1)+1) 	&
						- 2. * dtm_ / r(nz) * phi * F1D(nz-1,3*(i-1)+1)
			U1D(nz,3*(i-1)+2) = 0. !momentum
			U1D(nz,3*(i-1)+3) = U1D(nz-1,3*(i-1)+3) - lm * F1D(nz-1,3*(i-1)+3)  &
						- 2. * dtm_ / r(nz) * phi * F1D(nz-1,3*(i-1)+3)
		enddo
		!now electrons
		U1D(nz,neqi+1) =  U1D(nz-1,neqi+1)
	endif

	!BC for predictor and corrector
	U1D_p(1,:) = U1D(1,:)
	U1D_c(1,:) = U1D(1,:)
	U1D_p(nz,:) = U1D(nz,:)
	U1D_c(nz,:) = U1D(nz,:)

end subroutine apply_BC

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!    SUBROUTINE ART_VISCOSITY    !!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine artviscosity()

	use constants
	implicit none
	integer :: m,j,k
	real*8 :: a, b, y1, y2, x1, x2

	if (geom=="slab") then
		do k = 1, nz
		   eps_visc(k,1:neqi+1) = eps_visc_max
		enddo
	else !spherical
		y1 = log10(eps_visc_max)
		y2 = log10(eps_visc_max / eps_compress)
		x1 = log10(r(nz))
		x2 = log10(r(1))
		b  = ( y1*x2 - y2*x1 ) / ( x2 - x1 )
		a  = ( y1 - b ) / x1
		do k = 1,nz
			eps_visc(k,1:neqi+1) = r(k)**a * 10**b
		enddo
		write(*,*) "eps_visc(1,1) = ", eps_visc(1,1), "eps_visc(nz,1) = ", eps_visc(nz,1)
	endif

end subroutine artviscosity

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!    SUBROUTINE  INITSPACEGRID    !!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine initspacegrid()
	use constants

	implicit none
	integer :: m,j,k
	real*8 :: a, b, y1, y2, x1, x2

	r(1) = L + rmin
	r0(:,nregions+1) = r0(:,nregions+1) + rmin !correct also r0 array

	do k = 2, nz
		r(k) = r(k-1) + dr !note that dr is negative
	enddo
	r(nz) = rmin

	!find out position of shell in r
	k = 1
	do while (r(k)>=r0(2,2)) !remember that r starts at Rmax
	  k = k + 1
	enddo
	nz0 = k

end subroutine initspacegrid

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!    SUBROUTINE  INITVARIABLES    !!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine initvariables()

	use constants

	implicit none
	integer :: m,j,k

	do j = 1, nspec
		do m = 1, nregions
			do k = 1, nz
				if ( r(k) > r0(j,m) .and. r(k) <= r0(j,m+1)) then
					T0(k,j) = temp0(j,m)
					N0(k,j) = den0(j,m)
					V0(k,j) = vel0(j,m)
				endif
			enddo
		enddo
	enddo

	rho = 0.
	u   = 0.
	p   = 0.
	T   = 0.
	nz00 = nz0 - dnz

	!ions
    do k = 1, nz
    	do j = 1, nspec
			rho(k,j) = N0(k,j) * mi(j)
			u(k,j) = V0(k,j)
			p(k,j) = rho(k,j) / mi(j) * T0(k,j)
			T(k,j) = T0(k,j)
		enddo
    enddo

end subroutine initvariables

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!    SUBROUTINE  DO_SMOOTHING    !!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine do_smoothing()

	use constants

	implicit none
	integer :: m, j, k, kgrad(2)=0
	real*8 :: a, b, y1, y2, x1, x2
	real*8 :: grad, mingrad, maxgrad
	real*8 :: nA(2)=0., nB(2)=0., uA(2)=0., uB(2)=0., TA(2)=0., TB(2)=0.
	integer :: sgn1(2) = (/0,1/), sgn2(2) = (/1,0/), ok(2) = 0
	real*8, parameter :: epsilon = 1.e-10
	real*8 :: vel(nz)

	if (smoothing) then
		!now smooth out interface between nz0 and nz0+nsmooth
		!------------------------
		write(*,*) "density smoothing function"
		write(*,*) "--------------------------"
		!for each ion species, find position of max and min gradients, then smooth out density
		do j = 1, nspec
			maxgrad = 0.
			mingrad = 1.e20
			do k = 1, nz - 1
				grad = rho(k,j) / rho(k+1,j)
				if (grad > maxgrad) then
					maxgrad = grad
					kgrad(1) = k
					nA(1) = rho(k,j) / mi(j)
					nB(1) = rho(k+1,j) / mi(j)
				endif
				if (grad < mingrad) then
					mingrad = grad
					kgrad(2) = k
					nA(2) = rho(k,j) / mi(j)
					nB(2) = rho(k+1,j) / mi(j)
				endif
			enddo

			do m = 1,2
				do k = kgrad(m)-sgn1(m)*nsmooth, kgrad(m) + sgn2(m)*nsmooth
						!now smooth out interface
						y1 = log10( nA(m)  )
						y2 = log10( nB(m)  )
						x1 = log10( r(kgrad(m) - sgn1(m)*nsmooth) ) !note r is in micron
						x2 = log10( r(kgrad(m) + sgn2(m)*nsmooth) ) !note r is in micron
						b  = ( y1*x2 - y2*x1 ) / ( x2 - x1 )
						a  = ( y1 - b ) / x1

						rho(k,j) = (  10**( a*log10(r(k)) + b )   ) * mi(j)
				enddo

			enddo
		enddo

		write(*,*) "velocity smoothing function"
		write(*,*) "--------------------------"
		!for each ion species, find position of max and min gradients, then smooth out density
		ok = 0
		do j = 1, nspec
			maxgrad = 0.
			mingrad = 1.e20
			vel = max(epsilon,abs(u(:,j)))
			do k = 1, nz - 1
				grad = vel(k) / vel(k+1)
				if (grad > maxgrad) then
					maxgrad = grad
					kgrad(1) = k
					uA(1) = vel(k)
					uB(1) = vel(k+1)
					if ( k+tsmooth < nz .and. k-tsmooth>0 ) then
						ok(1) = 1
					endif
				endif
				if (grad < mingrad) then
					mingrad = grad
					kgrad(2) = k
					uA(2) = vel(k)
					uB(2) = vel(k+1)
					if ( k+tsmooth < nz .and. k-tsmooth>0 ) then
						ok(2) = 1
					endif
				endif
			enddo

			do m = 1,2
				if(ok(m).ne.1) cycle !because it is a bad point...
				do k = kgrad(m)-sgn1(m)*vsmooth, kgrad(m) + sgn2(m)*vsmooth
						!now smooth out interface
						y1 = log10( uA(m) )
						y2 = log10( uB(m) )
						x1 = log10( r(kgrad(m) - sgn1(m)*vsmooth) ) !note r is in micron
						x2 = log10( r(kgrad(m) + sgn2(m)*vsmooth) ) !note r is in micron
						b  = ( y1*x2 - y2*x1 ) / ( x2 - x1 )
						a  = ( y1 - b ) / x1

						u(k,j) = - (  10**( a*log10(r(k)) + b )   )
				enddo

			enddo
		enddo

		write(*,*) "temperature smoothing function"
		write(*,*) "------------------------------"
		!for each ion species, find position of max and min gradients, then smooth out density
		ok = 0
		do j = 1, nspec
			maxgrad = 0.
			mingrad = 1.e20
			do k = 1, nz - 1
				grad = T(k,j) / T(k+1,j)
				if (grad > maxgrad) then
					maxgrad = grad
					kgrad(1) = k
					TA(1) = T(k,j)
					TB(1) = T(k+1,j)
					if ( k+tsmooth < nz .and. k-tsmooth>0 ) then
						ok(1) = 1
					endif
				endif
				if (grad < mingrad) then
					mingrad = grad
					kgrad(2) = k
					TA(2) = T(k,j)
					TB(2) = T(k+1,j)
					if ( k+tsmooth < nz .and. k-tsmooth>0 ) then
						ok(2) = 1
					endif
				endif
			enddo

			do m = 1,2
				if(ok(m).ne.1) cycle !because it is a bad point...
				do k = kgrad(m)-sgn1(m)*tsmooth, kgrad(m) + sgn2(m)*tsmooth
						!now smooth out interface
						y1 = log10( TA(m) )
						y2 = log10( TB(m) )
						x1 = log10( r(kgrad(m) - sgn1(m)*tsmooth) ) !note r is in micron
						x2 = log10( r(kgrad(m) + sgn2(m)*tsmooth) ) !note r is in micron
						b  = ( y1*x2 - y2*x1 ) / ( x2 - x1 )
						a  = ( y1 - b ) / x1

						T(k,j) = (  10**( a*log10(r(k)) + b )   )
				enddo

			enddo

		enddo

		!now update pressure after all this smoothing...
		do j = 1, nspec
			p(:,j) = rho(:,j)/mi(j) * T(:,j)
		enddo

	endif
end subroutine do_smoothing

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!------------- SUBROUTINE READ_FILES -----------------------
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine read_files()

	use constants

	implicit none

	real*8, dimension(nz) :: var
	real*8 :: Tmean
	integer :: i,k,c1,c2,c3,ierr
	character(len=20) :: filename, x1
	character(len=8) :: fmt ! format descriptor

	fmt = '(I1)' ! an integer of width 5 with zeros at the left
	c1 = 0
	c2 = 0
	c3 = 0

	write(*,*) "--------------------------------"
	do i = 1,(nspec+1)*3+2
		if (i==1) then
			filename = 'r.dat'
			write(*,*) "reading ", filename
			write(*,*)
			open(unit=unt, file = filename, action = 'read', iostat = ierr)
				if(ierr/=0) then
					write(*,*) 'problems in opening file r.dat'
					stop
				else
					do k = 1, nz
						read(unt,'(E20.10E3)') r(k)
					enddo
				endif
			close(unt)
		else if (i==(nspec+1)*3+2) then
			filename = 'efield.dat'
			write(*,*) "reading ", filename
			open(unit=unt, file = filename, action = 'read', iostat = ierr)
				if(ierr/=0) then
					write(*,*) '  warning - efield.dat does not exist'
				else
					do k = 1, nz
						read(unt,'(E20.10E3)') Efield(k)
					enddo
				endif
			close(unt)
		else if (i>1 .and. i<=(nspec+1)+1) then
			c1 = c1 +1
			write (x1,fmt) c1
			filename = trim('vel'//trim(x1)//'.dat')
			write(*,*) "reading ", filename
			open(unit=unt, file = filename, action = 'read', iostat = ierr)
				if(ierr/=0) then
					write(*,*) 'problems in opening file', filename
					stop
				else
					do k = 1, nz
						read(unt,'(E20.10E3)') u(k,c1)  !m/s
					enddo
				endif
			close(unt)
			if (i==nspec+2) write(*,*)
		else if (i>(nspec+1)+1 .and. i<=2*(nspec+1)+1) then
			c2 = c2 +1
			write (x1,fmt) c2
			filename = trim('rho'//trim(x1)//'.dat')
			write(*,*) "reading ", filename
			open(unit=unt, file = filename, action = 'read', iostat = ierr)
				if(ierr/=0) then
					write(*,*) 'problems in opening file', filename
					stop
				else
					do k = 1, nz
						read(unt,'(E20.10E3)') rho(k,c2)  !kg/cm3
					enddo
				endif
			close(unt)
			if (i==2*(nspec+1)+1) write(*,*)
		else if (i>2*(nspec+1)+1 .and. i<=3*(nspec+1)+1) then
			c3 = c3 +1
			write (x1,fmt) c3
			filename = trim('temp'//trim(x1)//'.dat')
			write(*,*) "reading ", filename
			open(unit=unt, file = filename, action = 'read', iostat = ierr)
				if(ierr/=0) then
					write(*,*) 'problems in opening file', filename
					stop
				else
					do k = 1, nz
						read(unt,'(E20.10E3)') T(k,c3)  !Joule
					enddo
				endif
			close(unt)
			if (i==3*(nspec+1)+1) write(*,*)
		endif
	enddo
	write(*,*) "--------------------------------"

	!define position of boundary
	L = r(1)

	!make temperature equal at border
	Tmean = sum(T(1,:))/(nspec+1)
	T(1,:) = Tmean

end subroutine read_files

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!---------------------- SUBROUTINE DEFINE FLUXES ----------------------------
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine define_fluxes()

	use constants

	implicit none
	integer :: k,j

	U1D = 0.
	F1D = 0.
	R1D = 0.

	write(*,*)
	write(*,*) " WARNING - forcing quasi-neutrality and zero-current conditions"
	write(*,*)

	rho(:,nspec+1) = 0.
	u(:,nspec+1) = 0.

	!electrons
    do k = 1, nz
    	do j = 1, nspec
			rho(k,nspec+1) = rho(k,nspec+1) + me * Zi(j) * rho(k,j) / mi(j) !quasi-neutrality
			u(k,nspec+1) = u(k,nspec+1) + Zi(j) * rho(k,j) * u(k,j) / mi(j)   !zero-current condition
    	enddo
    	if(.not.(restart)) then !take species 1 as reference
			p(k,nspec+1) = rho(k,nspec+1) / me * T(k,1) !electron pressure
			T(k,nspec+1) = T(k,1) !electron temperature
		else !we have read the electron temperature from file
			p(k,nspec+1) = rho(k,nspec+1) / me * T(k,nspec+1) !electron pressure
		endif
    enddo

!   correct for electron velocity
    u(:,nspec+1) = u(:,nspec+1) / ( rho(:,nspec+1) / me )
!	apply BC at piston side
!    do j = 1, nspec+1
!        u(1,j) = V0(1,j)!0.
!    enddo

!	finally, calculate U1D and F1D for all species
	do j = 1, nspec
		U1D(:,3*(j-1)+1) = rho(:,j)
		U1D(:,3*(j-1)+2) = rho(:,j) * u(:,j)
		U1D(:,3*(j-1)+3) = p(:,j) / (g-1) + 0.5 * rho(:,j) * u(:,j)**2
		F1D(:,3*(j-1)+1) = rho(:,j) * u(:,j)
		F1D(:,3*(j-1)+2) = rho(:,j) * u(:,j)**2 + p(:,j)
		F1D(:,3*(j-1)+3) = u(:,j) * ( g / (g-1) * p(:,j) + 0.5 * rho(:,j) * u(:,j)**2 )
	enddo
	!electrons
	U1D(:,neqi+1) = p(:,nspec+1) / (g-1) + 0.5 * rho(:,nspec+1) * u(:,nspec+1)**2
	F1D(:,neqi+1) = u(:,nspec+1) * ( g / (g-1) * p(:,nspec+1) + 0.5 * rho(:,nspec+1) * u(:,nspec+1)**2 )

end subroutine define_fluxes

max_stars_repo_name: cbellei/mION

max_issues_repo_name: cbellei/mION

max_forks_repo_name: cbellei/mION


module geometry

  use constants,  only: dp
  use cell,       only: cell_type
  implicit none

  private

  ! default geom opt parameters:
  integer,        parameter ::  default_max_iter = 100
  real(kind=dp),  parameter ::  default_tpsd_initial_lambda = 1.0E-05_dp
  real(kind=dp),  parameter ::  default_max_disp = 5.0E-01_dp             ! Ang
  real(kind=dp),  parameter ::  default_force_tol = 1.0E-05_dp            ! eV/Ang
  real(kind=dp),  parameter ::  default_energy_tol = 1.0E-08_dp           ! eV
  real(kind=dp),  parameter ::  default_disp_tol = 1.0E-03_dp             ! Ang
  real(kind=dp),  parameter ::  default_stress_tol = 1.0E-01_dp           ! GPa
  logical,        parameter ::  default_fix_cell = .false.
  logical,        parameter ::  default_fix_atoms = .false.

  ! hardcoded for now - possibly get from input?
  real(kind=dp),  parameter   ::  d_epsilon = 5.0E-07_dp

  type, public  ::  geom_type
    integer                         ::  output_unit_num

    real(kind=dp),  dimension(3,3)              ::  orig_lattice_vectors

    real(kind=dp),  dimension(3,3)              ::  strain_tensor
    real(kind=dp),  dimension(3,3)              ::  stress_tensor
    real(kind=dp),  dimension(:,:), allocatable ::  inv_hessian
    real(kind=dp),  dimension(:,:), allocatable ::  atom_forces
    real(kind=dp),  dimension(:),   allocatable ::  x_vec
    real(kind=dp),  dimension(:),   allocatable ::  f_vec
    real(kind=dp),  dimension(:),   allocatable ::  x_vec_old
    real(kind=dp),  dimension(:),   allocatable ::  f_vec_old
    real(kind=dp),  dimension(:),   allocatable ::  delta_vec
    real(kind=dp)                               ::  lambda
    real(kind=dp)                               ::  pressure
    real(kind=dp)                               ::  enthalpy
    real(kind=dp)                               ::  potential_energy

    ! set from input
    logical                         ::  fix_cell
    logical                         ::  fix_atoms
    integer                         ::  max_iter
    real(kind=dp)                   ::  force_tol
    real(kind=dp)                   ::  energy_tol
    real(kind=dp)                   ::  disp_tol
    real(kind=dp)                   ::  stress_tol
    real(kind=dp)                   ::  max_disp
    real(kind=dp)                   ::  tpsd_initial_lambda
    real(kind=dp),  dimension(3,3)  ::  ext_pressure_tensor
  end type geom_type

  public  ::  geom_init
  public  ::  do_geom_opt
  public  ::  geom_deallocate
  public  ::  geom_calc_stress_tensor

contains

subroutine geom_init(geom, cell)
  use io,         only: io_open_file, seedname
  use constants,  only: units_natural_to_atomic
  implicit none
  type(geom_type),  intent(inout) ::  geom
  type(cell_type),  intent(inout) ::  cell

  ! get parameters from input file
  call geom_read_input(geom)

  ! convert parameters into atomic units:
  geom%force_tol = units_natural_to_atomic(force=geom%force_tol)
  geom%energy_tol = units_natural_to_atomic(energy=geom%energy_tol)
  geom%disp_tol = units_natural_to_atomic(length=geom%disp_tol)
  geom%stress_tol = units_natural_to_atomic(pressure=geom%stress_tol)
  geom%max_disp = units_natural_to_atomic(length=geom%max_disp)

  geom%ext_pressure_tensor(:,:) = units_natural_to_atomic(pressure=geom%ext_pressure_tensor(:,:))

  ! allocate everything
  call geom_allocate(geom, cell)

  geom%orig_lattice_vectors(:,:) = cell%lattice_vectors(:,:)

  ! pressure = Tr(pressure_tensor)/3
  geom%pressure = (geom%ext_pressure_tensor(1,1) + geom%ext_pressure_tensor(2,2) + geom%ext_pressure_tensor(3,3))/3.0_dp

  geom%strain_tensor(:,:) = 0.0_dp

  call geom_inv_hessian_init(geom)

  call io_open_file(trim(seedname)//'.geom', geom%output_unit_num, 'replace')
end subroutine geom_init

subroutine geom_allocate(geom, cell)
  use io, only: io_err
  implicit none
  type(geom_type),  intent(inout) ::  geom
  type(cell_type),  intent(in)    ::  cell
  integer ::  array_size, istat

  array_size = 9+3*cell%natoms

  allocate(geom%x_vec(array_size), stat=istat)
  if (istat .ne. 0) call io_err("geom_allocate: Could not allocate x_vec array")

  allocate(geom%x_vec_old(array_size), stat=istat)
  if (istat .ne. 0) call io_err("geom_allocate: Could not allocate x_vec_old array")

  allocate(geom%f_vec(array_size), stat=istat)
  if (istat .ne. 0) call io_err("geom_allocate: Could not allocate f_vec array")

  allocate(geom%f_vec_old(array_size), stat=istat)
  if (istat .ne. 0) call io_err("geom_allocate: Could not allocate f_vec_old array")

  allocate(geom%inv_hessian(array_size, array_size), stat=istat)
  if (istat .ne. 0) call io_err("geom_allocate: Could not allocate inv_hessian array")

  allocate(geom%delta_vec(array_size), stat=istat)
  if (istat .ne. 0) call io_err("geom_allocate: Could not allocate delta_vec array")

  allocate(geom%atom_forces(3, cell%natoms), stat=istat)
  if (istat .ne. 0) call io_err("geom_allocate: Could not allocate atom_forces array")
end subroutine geom_allocate

subroutine geom_deallocate(geom)
  use io, only: io_err
  implicit none
  type(geom_type),  intent(inout) ::  geom
  integer ::  istat

  deallocate(geom%x_vec, stat=istat)
  if (istat .ne. 0) call io_err("geom_deallocate: Could not deallocate x_vec array")

  deallocate(geom%x_vec_old, stat=istat)
  if (istat .ne. 0) call io_err("geom_deallocate: Could not deallocate x_vec_old array")

  deallocate(geom%f_vec, stat=istat)
  if (istat .ne. 0) call io_err("geom_deallocate: Could not deallocate f_vec array")

  deallocate(geom%f_vec_old, stat=istat)
  if (istat .ne. 0) call io_err("geom_deallocate: Could not deallocate f_vec_old array")

  deallocate(geom%inv_hessian, stat=istat)
  if (istat .ne. 0) call io_err("geom_deallocate: Could not deallocate inv_hessian array")

  deallocate(geom%delta_vec, stat=istat)
  if (istat .ne. 0) call io_err("geom_deallocate: Could not deallocate delta_vec array")

  deallocate(geom%atom_forces, stat=istat)
  if (istat .ne. 0) call io_err("geom_deallocate: Could not deallocate atom_forces array")
end subroutine geom_deallocate

subroutine geom_inv_hessian_init(geom)
  implicit none
  type(geom_type),  intent(inout) ::  geom
  integer ::  i

  ! for now, let's just initialize to the identity.. can eventually use Pfrommer method
  geom%inv_hessian(:,:) = 0.0_dp

  do i = 1, size(geom%inv_hessian, 1)
    geom%inv_hessian(i,i) = 1.0_dp
  end do
end subroutine geom_inv_hessian_init

subroutine geom_read_input(geom)
  use io, only: io_input_get_single_value, max_line_len, io_input_get_data, io_str_get_num_tokens, io_err
  implicit none
  type(geom_type),  intent(inout) ::  geom

  character(len=max_line_len) ::  keyword
  character(len=max_line_len),  allocatable,  dimension(:)  ::  input_data
  logical ::  kw_found
  integer ::  iline, istat

  ! single value parameters:
  call io_input_get_single_value('geom_max_iter', geom%max_iter, default_max_iter)
  if (geom%max_iter .le. 0) call io_err("geom_read_input: geom_max_iter must be > 0")

  call io_input_get_single_value('geom_tpsd_initial_lambda', geom%tpsd_initial_lambda, default_tpsd_initial_lambda)
  if (geom%tpsd_initial_lambda .le. 0.0_dp) call io_err("geom_read_input: geom_tpsd_initial_lambda must be > 0.0")

  call io_input_get_single_value('geom_force_tol', geom%force_tol, default_force_tol)
  if (geom%force_tol .le. 0.0_dp) call io_err("geom_read_input: geom_force_tol must be > 0.0")

  call io_input_get_single_value('geom_energy_tol', geom%energy_tol, default_energy_tol)
  if (geom%energy_tol .le. 0.0_dp) call io_err("geom_read_input: geom_energy_tol must be > 0.0")

  call io_input_get_single_value('geom_disp_tol', geom%disp_tol, default_disp_tol)
  if (geom%disp_tol .le. 0.0_dp) call io_err("geom_read_input: geom_disp_tol must be > 0.0")

  call io_input_get_single_value('geom_stress_tol', geom%stress_tol, default_stress_tol)
  if (geom%stress_tol .le. 0.0_dp) call io_err("geom_read_input: geom_stress_tol must be > 0.0")

  call io_input_get_single_value('geom_max_disp', geom%max_disp,  default_max_disp)
  if (geom%max_disp .le. 0.0_dp) call io_err("geom_read_input: geom_max_disp must be > 0.0")

  call io_input_get_single_value('geom_fix_cell', geom%fix_cell,  default_fix_cell)

  call io_input_get_single_value('geom_fix_atoms', geom%fix_atoms,  default_fix_atoms)

  ! block parameters:
  keyword = 'external_pressure'
  call io_input_get_data(keyword, input_data, kw_found)
  if (.not. kw_found) then
    ! zero external pressure:
    geom%ext_pressure_tensor(:,:) = 0.0_dp
  else
    ! input_data must be allocated
    if (size(input_data, 1) .ne. 3) call io_err("geom_read_input: "//trim(keyword)//" block must contain 3 lines")

    ! don't do anything smart - require whole 3x3 array - don't check that it's valid..
    do iline = 1, 3
      if (io_str_get_num_tokens(input_data(iline)) .ne. 3) &
        & call io_err("geom_read_input: "//trim(keyword)//" block must contain 3 values on each line")

      read(input_data(iline), *, iostat=istat) geom%ext_pressure_tensor(iline, :)
      if (istat .ne. 0) call io_err("geom_read_input: "//trim(keyword)//": Error reading line")
    end do
  end if

  ! clean up
  if (allocated(input_data)) then
    deallocate(input_data, stat=istat)
    if (istat .ne. 0) call io_err("geom_read_input: Could not deallocate input_data array")
  end if
end subroutine geom_read_input

subroutine do_geom_opt(geom, cell)
  use constants,  only: units_atomic_to_natural, element_symbol
  use cell,       only: cell_min_img_vec
  use potential,  only: pot_get_forces, pot_get_potential
  use io,         only: io_err, stdout
  implicit none
  type(geom_type),  intent(inout) ::  geom
  type(cell_type),  intent(inout) ::  cell

  real(kind=dp),  allocatable,  dimension(:,:)  ::  old_cart_pos
  real(kind=dp),  dimension(3)  ::  disp_vec
  real(kind=dp)   ::  old_enthalpy, delta_enthalpy
  real(kind=dp)   ::  max_force, max_stress_comp, max_disp
  real(kind=dp)   ::  current_force_sq, current_disp
  integer         ::  i, iter

  ! routine does TPSD optimization, fixed/variable cell, but no other constraints..

  ! need array to store old cart positions for calculating displacements if fixed cell
  if (geom%fix_cell) then
    allocate(old_cart_pos(3, cell%natoms), stat=i)
    if (i .ne. 0) call io_err("do_geom_opt: Could not allocate old_cart_pos array")

    old_cart_pos(:,:) = cell%atom_cart_pos(:,:)
  end if

  call geom_build_x_vec(geom, cell)

  if (geom%fix_atoms) then
    call pot_get_potential(cell, geom%potential_energy)
    geom%atom_forces = 0.0_dp
  else
    call pot_get_forces(cell, geom%atom_forces, geom%potential_energy)
  end if

  if (geom%fix_cell) then
    geom%stress_tensor(:,:) = 0.0_dp
  else
    call geom_calc_stress_tensor(geom%stress_tensor, cell)
  end if

  write(stdout,*) "Starting geometry optimisation,"
  write(stdout,*) "  Default energy_tol = ", units_atomic_to_natural(energy=geom%energy_tol)  , "eV"
  write(stdout,*) "  Default force_tol  = ", units_atomic_to_natural(force=geom%force_tol)    , "eV/Ang"
  write(stdout,*) "  Default disp_tol   = ", units_atomic_to_natural(length=geom%disp_tol)    , "Ang"
  if(.not. geom%fix_cell) write(stdout,*) "  Default stress_tol = ", units_atomic_to_natural(pressure=geom%stress_tol), "GPa"

  call geom_build_f_vec(geom, cell)
  ! geom%atom_forces are not in cartesian space at this point

  geom%enthalpy = geom%potential_energy + geom%pressure*cell%volume

  geom%lambda = geom%tpsd_initial_lambda

  write(stdout,*)
  write(stdout,*) "  Initial Lambda     =", geom%lambda
  write(stdout,*) "  Initial enthalpy   =", units_atomic_to_natural(energy  =geom%enthalpy), "eV"
  write(stdout,*) "  Initial max_force  =", units_atomic_to_natural(force   =maxval(geom%atom_forces)), "eV/Ang"
if(.not. geom%fix_cell) write(stdout,*)"  Initial max_stress =",units_atomic_to_natural(pressure=maxval(geom%stress_tensor)), "GPa"

  ! just to start the whole thing, make copies of x_vec and f_vec
  geom%x_vec_old = geom%x_vec
  geom%f_vec_old = geom%f_vec
  old_enthalpy = geom%enthalpy

  iter = 0
  call geom_write_output(geom, cell, iter)
  do !i = 1, geom%max_iter

    iter = iter + 1
    write(stdout,*)
    write(stdout,*) "Beginning TPSD iteration", iter

    ! X_(i+1) = X_i + lambda*H_i*F_i
    geom%x_vec = geom%x_vec + geom%lambda*matmul(geom%inv_hessian, geom%f_vec)

    ! with new positions, can now update everything else
    call geom_x_vec_to_cell(geom, cell)
    ! geom_x_vec_to_cell also updates geom%strain_tensor..

    ! should strain tensor be reinitialized to zero at every step?

    if (geom%fix_atoms) then
      call pot_get_potential(cell, geom%potential_energy)
      geom%atom_forces = 0.0_dp
    else
      call pot_get_forces(cell, geom%atom_forces, geom%potential_energy)
    end if

    ! not sure which "pressure" to use here.. External pressure?
    geom%enthalpy = geom%potential_energy + geom%pressure*cell%volume

    if (geom%fix_cell) then
      geom%stress_tensor(:,:) = 0.0_dp
    else
      call geom_calc_stress_tensor(geom%stress_tensor, cell)
    end if

    ! get max abs(force) before we build new f_vec
    max_force = 0.0_dp
    do i = 1, cell%natoms
      current_force_sq = dot_product(geom%atom_forces(:,i), geom%atom_forces(:,i))
      if (current_force_sq .gt. max_force) max_force = current_force_sq
    end do
    max_force = sqrt(max_force)

    ! max abs(stress) component - only calc if variable cell
    if (.not. geom%fix_cell) then
      max_stress_comp = maxval(abs(geom%stress_tensor/cell%volume+geom%ext_pressure_tensor))
    end if

    ! calc max abs(disp).. Always use same cell, but must use cell_min_img
    if (geom%fix_cell) then
      max_disp = 0.0_dp
      do i = 1, cell%natoms
        call cell_min_img_vec(cell, old_cart_pos(:,i), cell%atom_cart_pos(:,i), disp_vec, current_disp)
        if (current_disp .gt. max_disp) max_disp = current_disp
      end do
    end if

    ! enforce max disp? Don't want to move too far in a single step..

    ! calc change in enthalpy/energy:
    delta_enthalpy = old_enthalpy - geom%enthalpy

!!!!!!!!!!! to do
    ! can now write output...
    call geom_write_output(geom, cell, iter)
    write(stdout,*) "  Lambda     =", geom%lambda
    write(stdout,*) "  Enthalpy   =", units_atomic_to_natural(energy  =geom%enthalpy), 'eV'
    if(abs(delta_enthalpy) <= geom%energy_tol) then
      write(stdout,*) "  Y Max_dH     =", units_atomic_to_natural(energy  =delta_enthalpy), 'eV'
    else
      write(stdout,*) "  N Max_dH     =", units_atomic_to_natural(energy  =delta_enthalpy), 'eV'
    end if

    if(max_force <= geom%force_tol) then
      write(stdout,*) "  Y Max_force  =", units_atomic_to_natural(force   =max_force), 'eV/Ang'
    else
      write(stdout,*) "  N Max_force  =", units_atomic_to_natural(force   =max_force), 'eV/Ang'
    end if
    if(geom%fix_cell) then
      if(max_disp <= geom%disp_tol) then
        write(stdout,*) "  Y  Max_disp   =", units_atomic_to_natural(length  =max_disp), "Ang"
      else
        write(stdout,*) "  N  Max_disp   =", units_atomic_to_natural(length  =max_disp), "Ang"
      end if
    else
      if(max_stress_comp <= geom%stress_tol) then
        write(stdout,*) "  Y Max_stress =", units_atomic_to_natural(pressure=max_stress_comp), 'GPa'
      else
        write(stdout,*) "  N Max_stress =", units_atomic_to_natural(pressure=max_stress_comp), 'GPa'
      end if
    end if

!    write(stdout,*) 'Stress:                                   |  Strain'
!    do i=1,3
!      write(stdout,*) geom%stress_tensor(i,:), '|', geom%strain_tensor(i,:)
!    end do
    ! exit conditions
    if( (max_force <= geom%force_tol) .and. (abs(delta_enthalpy) <= geom%energy_tol) &
      .and. (max_disp <= geom%disp_tol) .and. (geom%fix_cell .or. (max_stress_comp <= geom%stress_tol))) then
      write(stdout,*)
      write(stdout,*) "Converged in ", iter, "iterations"
      write(stdout,*)
      write(stdout,*) "====== Final cell ======"
      write(stdout,*) "begin lattice_vectors"
      do i = 1, 3
        write(stdout,*) units_atomic_to_natural(length=cell%lattice_vectors(:,i))
      end do
      write(stdout,*) "end lattice_vectors"
      write(stdout,*)
      write(stdout,*) "begin atomic_positions"
      do i = 1, cell%natoms
        write(stdout,*) element_symbol(cell%atom_species(i)), cell%atom_frac_pos(:,i)
      end do
      write(stdout,*) "end atomic_positions"
      exit
    else if (iter .eq. geom%max_iter) then
      write(stdout,*)
      write(stdout,*) 'Failed to converge in', iter, 'iterations'
      exit
    end if
!!!!!!!!!! end to do

    ! below here necessary for next step only:
    call geom_build_f_vec(geom, cell)
    ! geom%atom_forces are not in cartesian space at this point

    ! for TPSD, let's not update geom%inv_hessian - this will need to be changed for other methods (eg: BFGS)
    ! calculate lambda for next time
    call geom_tpsd_calc_step_size(geom)

    ! store things for next time..
    if (geom%fix_cell) old_cart_pos(:,:) = cell%atom_cart_pos(:,:)
    geom%x_vec_old = geom%x_vec
    geom%f_vec_old = geom%f_vec
    old_enthalpy = geom%enthalpy
  end do
end subroutine do_geom_opt

subroutine geom_tpsd_calc_step_size(geom)
  implicit none
  type(geom_type),  intent(inout) ::  geom

  real(kind=dp) ::  dx_dot_dg, dg_dot_dg

  ! TPSD:
  ! lambda = <delta_x, delta_g>/<delta_g, delta_g>
  ! (or lambda = <delta_x, delta_x>/<delta_x, delta_g>)

  ! and:
  ! delta_g = -delta_f = -( F_i - F_(i-1) )

  ! in Pfrommer this is: -(H_i*F_i - H_(i-1)*F_(i-1)) because we multiply by H
  ! if H_i = H_(i-1) = H_0     =>     delta_g = -H_0*( F_i - F_(i-1) )

  ! use geom%delta_vec as a work array..
  geom%delta_vec = geom%f_vec - geom%f_vec_old
  geom%delta_vec = -matmul(geom%inv_hessian, geom%delta_vec)

  dx_dot_dg = dot_product((geom%x_vec - geom%x_vec_old), geom%delta_vec)
  dg_dot_dg = dot_product(geom%delta_vec, geom%delta_vec)

  geom%lambda = dx_dot_dg/dg_dot_dg
end subroutine geom_tpsd_calc_step_size

subroutine geom_build_x_vec(geom, cell)
  implicit none
  type(geom_type),  intent(inout) ::  geom
  type(cell_type),  intent(in)    ::  cell

  integer ::  i, j, k

  do i = 1, 3 ! row
    do j = 1, 3 ! col
      geom%x_vec(3*(i-1)+j) = geom%strain_tensor(i,j)
    end do
  end do

  ! s is the fractional coordinates..
  ! start at element 10..
  k = 9
  do i = 1, cell%natoms
    do j = 1, 3
      k = k + 1
      geom%x_vec(k) = cell%atom_frac_pos(j,i)
    end do
  end do
end subroutine geom_build_x_vec

subroutine geom_x_vec_to_cell(geom, cell)
  implicit none
  type(geom_type),  intent(inout) ::  geom
  type(cell_type),  intent(inout) ::  cell

  integer :: i, j, k

  do i = 1, 3 ! row
    do j = 1, 3 ! col
      geom%strain_tensor(i,j) = geom%x_vec(3*(i-1)+j)
    end do
  end do

  ! s is the fractional coordinates..
  ! start at element 10..
  k = 9
  do i = 1, cell%natoms
    do j = 1, 3
      k = k + 1
      cell%atom_frac_pos(j,i) = geom%x_vec(k)
    end do
  end do

  call geom_apply_strain(geom%strain_tensor, geom%orig_lattice_vectors, cell)
end subroutine geom_x_vec_to_cell

subroutine geom_apply_strain(strain, orig_lattice, cell)
  implicit none
  real(kind=dp),  dimension(3,3), intent(in)    ::  strain
  real(kind=dp),  dimension(3,3), intent(in)    ::  orig_lattice
  type(cell_type),                intent(inout) ::  cell

  real(kind=dp),  dimension(3,3)  ::  identity
  integer ::  i

  ! apply strain to change lattice vectors - resync everything else (cart positions etc)
  identity = 0.0_dp
  do i = 1, 3
    identity(i,i) = 1.0_dp
  end do

  cell%lattice_vectors = matmul((identity + strain), orig_lattice)

  call geom_sync_cell(cell)
end subroutine geom_apply_strain

subroutine geom_sync_cell(cell)
  use cell, only: cell_calc_vol, cell_get_recip, cell_shift_to_unit_cell, cell_frac_to_cart, cell_min_img_init
  implicit none
  type(cell_type),  intent(inout) ::  cell

  ! assuming that we have either changed cell vectors (either directly or applying a strain)
  ! and/or changed fractional positions

  call cell_calc_vol(cell)
  call cell_get_recip(cell)             ! need vol first
  call cell_shift_to_unit_cell(cell)
  call cell_frac_to_cart(cell)
  call cell_min_img_init(cell)          ! cell size/shape may have changed
end subroutine geom_sync_cell

subroutine geom_build_f_vec(geom, cell)
  use constants,  only: two_pi
  use algor,      only: algor_calc_inverse_real
  implicit none

  type(geom_type),  intent(inout)   ::  geom
  type(cell_type),  intent(in)      ::  cell

  real(kind=dp),    dimension(3,3)  ::  metric_tensor
  real(kind=dp),    dimension(3,3)  ::  f_strain
  integer                           ::  i, j, k

  ! strain components of f_vec:
  ! f^(strain) = -(stress + pressure*vol)*(1 + strain^T)^-1

  ! first calculate rightmost term:
  ! f_strain is identity first..
  f_strain = 0.0_dp
  do i = 1, 3
    f_strain(i,i) = 1.0_dp
  end do

  f_strain = f_strain + transpose(geom%strain_tensor)
  call algor_calc_inverse_real(f_strain)

  f_strain = matmul(-(geom%stress_tensor + geom%ext_pressure_tensor*cell%volume), f_strain)

  ! now have strain components.. can put into F_vec
  ! same sort of form as the X_vec...
  do j = 1, 3
    do i = 1, 3
      geom%f_vec(3*(i-1)+j) = f_strain(i,j)
    end do
  end do

  ! force terms:
  ! "obtained by multiplying the forces on the atoms in the lattice coordinates with the metric tensor"

  ! forces are currently in cartesian space -> convert to fractional..
  geom%atom_forces = matmul(cell%recip_lattice/two_pi, geom%atom_forces)

  metric_tensor = matmul(transpose(cell%lattice_vectors), cell%lattice_vectors)

  ! multiply by metric tensor
  geom%atom_forces = matmul(metric_tensor, geom%atom_forces)

  ! start at element 10..
  k = 9
  do j = 1, cell%natoms
    do i = 1, 3
      k = k + 1
      geom%f_vec(k) = geom%atom_forces(i,j)
    end do
  end do
end subroutine geom_build_f_vec

subroutine geom_calc_stress_tensor(stress_tensor, cell)
  use cell,       only: cell_copy, cell_deallocate
  use potential,  only: pot_get_potential
  implicit none
  real(kind=dp),    dimension(3,3), intent(out) ::  stress_tensor
  type(cell_type),                  intent(in)  ::  cell

  type(cell_type)                 ::  tmp_cell
  real(kind=dp),  dimension(3,3)  ::  d_strain
  real(kind=dp)                   ::  energy
  integer ::  i, j, k

  stress_tensor(:,:) = 0.0_dp

  do i = 1, 3
    do j = i, 3
      ! centered difference requires perturbations in each direction:
      do k = -1, 1, 2 ! -1 = back,  +1 = forward
        call cell_copy(cell, tmp_cell)

        ! d_strain is epsilon' tensor
        d_strain = 0.0_dp
        d_strain(i,j) = real(k,dp)*d_epsilon

        call geom_apply_strain(d_strain, tmp_cell%lattice_vectors, tmp_cell)

        ! calc new energy - don't need forces
        call pot_get_potential(tmp_cell, energy)

        ! sum energies..
        stress_tensor(i,j) = stress_tensor(i,j) + real(k,dp)*energy
      end do

      ! use centered difference
      ! dy/dx = (y(x+h) - y(x-h)) / 2*h
      stress_tensor(i,j) = stress_tensor(i,j)/(2.0_dp*d_epsilon*cell%volume)

      if (i .ne. j) stress_tensor(j,i) = stress_tensor(i,j)
    end do
  end do

  call cell_deallocate(tmp_cell)
end subroutine geom_calc_stress_tensor

subroutine geom_write_output(geom, cell, iter)
  use constants,  only: element_symbol
  implicit none
  type(geom_type),  intent(in)  ::  geom
  type(cell_type),  intent(in)  ::  cell
  integer,          intent(in)  ::  iter
  ! local vars
  real(kind=dp),  dimension(3,3)  ::  stress
  integer :: i

  if (iter .eq. 0) then
    write (geom%output_unit_num, 99) 'BEGIN header'
    write (geom%output_unit_num, *)  cell%natoms
    write (geom%output_unit_num, 99) 'END header'
    write (geom%output_unit_num, *)  ' '
  end if

  write (geom%output_unit_num, 1) iter!, dE_str, Fmax_str, dRmax_str, Smax_str
  write (geom%output_unit_num, 2) geom%potential_energy, geom%enthalpy

  ! lattice vectors (lines are a, b, c)
  do i = 1, 3
     write (geom%output_unit_num, 4) cell%lattice_vectors(:,i)
  end do

  if (.not. geom%fix_cell) then
     !cell stresses
     stress(:,:) = geom%stress_tensor(:,:)/cell%volume + geom%ext_pressure_tensor(:,:)
     do i=1,3
        write (geom%output_unit_num,6) stress(i,:)
     end do
  end if

  ! positions
  do i = 1, cell%natoms
    write (geom%output_unit_num, 7) element_symbol(cell%atom_species(i)), i, cell%atom_cart_pos(:,i)
  end do

  ! forces
  do i = 1, cell%natoms
    write (geom%output_unit_num, 9) element_symbol(cell%atom_species(i)), i, geom%atom_forces(:,i)
  end do

  write (geom%output_unit_num, *) ' '

  1   format(21x,i18,34x,4(3x,a),10x,'  <-- c')
  2   format(18x,2(3x,es24.16e3),27x,'  <-- E')
  4   format(18x,3(3x,es24.16e3),    '  <-- h')
  6   format(18x,3(3x,es24.16e3),    '  <-- S')
  7   format(1x,a8,1x,i8,3(3x,es24.16e3),'  <-- R')
  9   format(1x,a8,1x,i8,3(3x,es24.16e3),'  <-- F')
  99  format(1x,a)
end subroutine geom_write_output

end module geometry


text:
module sst_advection_boundary_average
#include <messenger.h>
    use mod_kinds,              only: rk,ik
    use mod_constants,          only: ZERO,ONE,TWO,THREE,HALF
    use mod_fluid,              only: omega, gam, Rgas
    use type_operator,          only: operator_t
    use type_chidg_worker,      only: chidg_worker_t
    use type_properties,        only: properties_t
    use DNAD_D
    implicit none

    !>
    !!
    !!
    !! @author  Eric M. Wolf
    !! @date    01/26/2018
    !!
    !--------------------------------------------------------------------------------
    type, extends(operator_t), public :: sst_advection_boundary_average_operator_t

    contains

        procedure   :: init
        procedure   :: compute

    end type sst_advection_boundary_average_operator_t
    !*****************************************************************************************

contains

    !>
    !!
    !!
    !! @author  Eric M. Wolf
    !! @date    01/26/2018
    !!
    !--------------------------------------------------------------------------------
    subroutine init(self)
        class(sst_advection_boundary_average_operator_t),   intent(inout)  :: self

        ! Set operator name
        call self%set_name('SST Advection Boundary Average Operator')

        ! Set operator type
        call self%set_operator_type('Boundary Advective Operator')

        ! Set operator equations
        call self%add_primary_field('Density * Omega')
        call self%add_primary_field('Density * k')
        call self%add_primary_field('Momentum-1')
        call self%add_primary_field('Momentum-2')
        call self%add_primary_field('Momentum-3')
        call self%add_primary_field('Energy')

    end subroutine init
    !*****************************************************************************************

    !>
    !!
    !!
    !! @author  Eric M. Wolf
    !! @date    01/26/2018
    !!
    !--------------------------------------------------------------------------------
    subroutine compute(self,worker,prop)
        class(sst_advection_boundary_average_operator_t),  intent(inout)   :: self
        type(chidg_worker_t),                                           intent(inout)   :: worker
        class(properties_t),                                            intent(inout)   :: prop

        type(AD_D), dimension(:), allocatable   ::  &
            density_m, density_p,                   &
            mom1_m, mom2_m, mom3_m,                 &
            mom1_p, mom2_p, mom3_p,                 &
            density_omega_m, density_omega_p,       &
            density_k_m, density_k_p,       &
            invdensity_m, invdensity_p,             &
            u_m, v_m, w_m,                          &
            u_p, v_p, w_p,                          &
            flux_1_m, flux_2_m, flux_3_m,           &
            flux_1_p, flux_2_p, flux_3_p, r

        !
        ! Interpolate solution to quadrature nodes
        !
        density_m         = worker%get_field('Density',           'value', 'face interior')
        density_p         = worker%get_field('Density',           'value', 'face exterior')

        mom1_m            = worker%get_field('Momentum-1',        'value', 'face interior')
        mom1_p            = worker%get_field('Momentum-1',        'value', 'face exterior')

        mom2_m            = worker%get_field('Momentum-2',        'value', 'face interior')
        mom2_p            = worker%get_field('Momentum-2',        'value', 'face exterior')

        mom3_m            = worker%get_field('Momentum-3',        'value', 'face interior')
        mom3_p            = worker%get_field('Momentum-3',        'value', 'face exterior')

        density_omega_m = worker%get_field('Density * Omega', 'value', 'face interior')
        density_omega_p = worker%get_field('Density * Omega', 'value', 'face exterior')

        !density_omega_m = density_m*density_omega_m
        !density_omega_p = density_p*density_omega_p

        density_k_m = worker%get_field('Density * k', 'value', 'face interior')
        density_k_p = worker%get_field('Density * k', 'value', 'face exterior')

        !density_k_m = density_m*density_k_m
        !density_k_p = density_p*density_k_p
        !
        ! Account for cylindrical. Get tangential momentum from angular momentum.
        !
        if (worker%coordinate_system() == 'Cylindrical') then
            r = worker%coordinate('1','face interior')
            mom2_m = mom2_m / r
            mom2_p = mom2_p / r
        end if

        !
        ! Get fluid advection velocity
        !
        invdensity_m = ONE/density_m
        u_m = mom1_m*invdensity_m
        v_m = mom2_m*invdensity_m
        w_m = mom3_m*invdensity_m

        invdensity_p = ONE/density_p
        u_p = mom1_p*invdensity_p
        v_p = mom2_p*invdensity_p
        w_p = mom3_p*invdensity_p

        !
        ! Compute average flux
        !
        flux_1_m = density_omega_m * u_m
        flux_2_m = density_omega_m * v_m
        flux_3_m = density_omega_m * w_m

        flux_1_p = density_omega_p * u_p
        flux_2_p = density_omega_p * v_p
        flux_3_p = density_omega_p * w_p

        !
        ! Integrate flux
        !
        call worker%integrate_boundary_average('Density * Omega','Advection', &
                                                flux_1_m, flux_2_m, flux_3_m,   &
                                                flux_1_p, flux_2_p, flux_3_p)

         !
        ! Compute average flux
        !
        flux_1_m = density_k_m * u_m
        flux_2_m = density_k_m * v_m
        flux_3_m = density_k_m * w_m

        flux_1_p = density_k_p * u_p
        flux_2_p = density_k_p * v_p
        flux_3_p = density_k_p * w_p

        !
        ! Integrate flux
        !
        call worker%integrate_boundary_average('Density * k','Advection', &
                                                flux_1_m, flux_2_m, flux_3_m,   &
                                                flux_1_p, flux_2_p, flux_3_p)

         !
        ! Compute average flux
        !
        flux_1_m = (TWO/THREE)*density_k_m
        flux_2_m = ZERO*(TWO/THREE)*density_k_m
        flux_3_m = ZERO*(TWO/THREE)*density_k_m

        flux_1_p = (TWO/THREE)*density_k_p
        flux_2_p = ZERO*(TWO/THREE)*density_k_p
        flux_3_p = ZERO*(TWO/THREE)*density_k_p

        !
        ! Integrate flux
        !
        call worker%integrate_boundary_average('Momentum-1','Advection', &
                                                flux_1_m, flux_2_m, flux_3_m,   &
                                                flux_1_p, flux_2_p, flux_3_p)
!
        ! Compute average flux
        !
        flux_1_m = ZERO*(TWO/THREE)*density_k_m
        flux_2_m = (TWO/THREE)*density_k_m
        flux_3_m = ZERO*(TWO/THREE)*density_k_m

        flux_1_p = ZERO*(TWO/THREE)*density_k_p
        flux_2_p = (TWO/THREE)*density_k_p
        flux_3_p = ZERO*(TWO/THREE)*density_k_p

        !
        ! Integrate flux
        !
        call worker%integrate_boundary_average('Momentum-2','Advection', &
                                                flux_1_m, flux_2_m, flux_3_m,   &
                                                flux_1_p, flux_2_p, flux_3_p)

!
        ! Compute average flux
        !
        flux_1_m = ZERO*(TWO/THREE)*density_k_m
        flux_2_m = ZERO*(TWO/THREE)*density_k_m
        flux_3_m = (TWO/THREE)*density_k_m

        flux_1_p = ZERO*(TWO/THREE)*density_k_p
        flux_2_p = ZERO*(TWO/THREE)*density_k_p
        flux_3_p = (TWO/THREE)*density_k_p

        !
        ! Integrate flux
        !
        call worker%integrate_boundary_average('Momentum-3','Advection', &
                                                flux_1_m, flux_2_m, flux_3_m,   &
                                                flux_1_p, flux_2_p, flux_3_p)

 !
        ! Compute average flux
        !
        flux_1_m = (TWO/THREE)*density_k_m * u_m
        flux_2_m = (TWO/THREE)*density_k_m * v_m
        flux_3_m = (TWO/THREE)*density_k_m * w_m

        flux_1_p = (TWO/THREE)*density_k_p * u_p
        flux_2_p = (TWO/THREE)*density_k_p * v_p
        flux_3_p = (TWO/THREE)*density_k_p * w_p

        !
        ! Integrate flux
        !
        call worker%integrate_boundary_average('Energy','Advection', &
                                                flux_1_m, flux_2_m, flux_3_m,   &
                                                flux_1_p, flux_2_p, flux_3_p)

    end subroutine compute
    !******************************************************************************************

end module sst_advection_boundary_average

max_stars_repo_name: nwukie/ChiDG

max_issues_repo_name: nwukie/ChiDG

max_forks_repo_name: nwukie/ChiDG


text:
Subroutine ro2et2(i1, i2, lbb1, lbb2, ei1, ei2, iblock, iseed)
  Parameter (maxstr=150001)
  Parameter (etam=0.5475, arho=0.77)
  Common /ee/id(maxstr), lb(maxstr)
  Common /ppb1/ene, factr2(6), fsum, ppinnb, s, wtot
  Common /ppmm/pprr, ppee, pppe, rpre, xopoe, rree
  Common /rndf77/nseed
  Save
  If (lb(i1)>=25 .And. lb(i1)<=27 .And. lb(i2)>=25 .And. lb(i2)<=27) Then
    iblock = 1895
    ei1 = etam
    ei2 = etam
    lbb1 = 0
    lbb2 = 0
  Else If (lb(i1)==0 .And. lb(i2)==0) Then
    iblock = 1896
    lbb1 = 25 + int(3*ranart(nseed))
    lbb2 = 25 + int(3*ranart(nseed))
    ei1 = arho
    ei2 = arho
  End If
  Return
End Subroutine ro2et2

max_stars_repo_name: xiaohaijin/AMPT

max_issues_repo_name: xiaohaijin/AMPT

max_forks_repo_name: xiaohaijin/AMPT


text:
program sbverify

  use globals
  use gdal
  use mapio
  use slezanbear
  use geo

  implicit none

  type dataentry
    real(dp) :: lat, lng
    real(dp) :: sky
  end type

  real(dp) :: llatsrc, llngsrc, ulatsrc, ulngsrc
  type(source), dimension(:,:), allocatable :: src

  type(dataentry), dimension(1024), target :: dat_all
  type(dataentry), dimension(:), pointer :: dat
  integer :: ndata, i, j
  type(modelpar) :: mpar

  character(128), dimension(9) :: fndem

  integer, parameter :: npar = 5
  real(dp), dimension(npar), parameter :: parlo &
  &   = [ 1e-6,  0.0, 0.5e3, 0.0, 2.8e-7]
  real(dp), dimension(npar), parameter :: parhi &
  &   = [10e-6, 24.0,  18e3, 0.8, 4.0e-7]
  real(dp), dimension(npar, npar + 1) :: par
  real(dp), dimension(npar + 1) :: modelerr
  real(dp), dimension(:,:), allocatable :: modelsky

  real(sp), dimension(:,:), allocatable :: mapi, maph
  type(geotransform) :: gti, gth
  integer :: errno

  character(*), parameter :: fni = 'SVDNB_npp_20150101-20151231_75N060W_vcm-orm-ntl_v10_c201701311200.avg_rade9.tif'
  character(*), parameter :: datafmt = '(F9.4,1X,F9.4,3X,F5.2,1X,F4.2,3X,F5.2)'
  character(128) :: buf

  call get_command_argument(1, buf)
  read (buf,*) mpar % alpha
  mpar % alpha = mpar % alpha * 1e-6

  call get_command_argument(2, buf)
  read (buf,*) mpar % relabs

  call get_command_argument(3, buf)
  read (buf,*) mpar % hscale

  call get_command_argument(4, buf)
  read (buf,*) mpar % beta

  call get_command_argument(5, buf)
  read (buf,*) mpar % skybg
  mpar % skybg = mpar % skybg * 1e-7

  write (*,'("parameters: ",F8.3,F7.3,F8.1,F6.3,F7.3)') &
        & mpar % alpha * 1e6, mpar % relabs, &
        & mpar % hscale, mpar % beta, (mpar % skybg) * 1e7

  terrain_attenuation = .false.

  ndata = 0
  open(33, file = 'data.txt', action = 'read')
  do i = 1,size(dat_all)
    read(33, *, iostat = errno) dat_all(i) % lat, dat_all(i) % lng, dat_all(i) % sky
    if (errno .ne. 0) exit
    ndata = ndata + 1
  end do
  close(33)

  dat => dat_all(1:ndata)

  llat = minval(dat % lat)
  ulat = maxval(dat % lat)
  llng = minval(dat % lng)
  ulng = maxval(dat % lng)

  fndem(1) = "eudem_dem_5deg_n45e010.tif"
  fndem(2) = "eudem_dem_5deg_n45e015.tif"
  fndem(3) = "eudem_dem_5deg_n45e020.tif"
  fndem(4) = "eudem_dem_5deg_n50e010.tif"
  fndem(5) = "eudem_dem_5deg_n50e015.tif"
  fndem(6) = "eudem_dem_5deg_n50e020.tif"
  fndem(7) = "eudem_dem_5deg_n55e010.tif"
  fndem(8) = "eudem_dem_5deg_n55e015.tif"
  fndem(9) = "eudem_dem_5deg_n55e020.tif"

  print '("dmax = ",F0.2," km")', dmax / 1e3
  call estimate_data_limits(llat, llng, ulat, ulng, dmax * 1.05, &
  &       llatsrc, llngsrc, ulatsrc, ulngsrc)
  print '("data limits = ",F6.2,1X,F7.2)', llatsrc, llngsrc, ulatsrc, ulngsrc

  call GDALAllRegister

  call map_read_section(fni, llatsrc, llngsrc, ulatsrc, ulngsrc, mapi, gti)
  write (0, '("satellite raster size = ", 2I6)') size(mapi,1), size(mapi,2)

  read_dem: block
    integer nox, noy

    call map_gen_size(llatsrc, llngsrc, ulatsrc, ulngsrc, elev_grid_meters, &
    &              nox, noy, gth)

    allocate( maph(nox,noy) )
    write (*,'("dem raster size = ",2I6)') size(maph,1), size(maph,2)

    call map_read_tiles(fndem, maph, gth)
  end block read_dem

  call random_seed()
  allocate(modelsky(ndata, npar + 1))

  printcalib: block
    real(dp) :: sky(4), hobs

    allocate( src(size(mapi,1), size(mapi,2)) )
    call mksources(mapi, gti, maph, gth, src)

    write (*,'(5A10)') 'LAT', 'LNG', 'DATA', 'MODEL', 'ERR'

    iterate_datpoints: do i = 1, ndata
      call onepoint(mpar, dat(i) % lat, dat(i) % lng, &
      &       src, maph, gth, sky, hobs)
      write (*,'(2F10.4,3F10.3)') &
      &       dat(i) % lat, dat(i) % lng, dat(i) % sky, ity2mag(sky(4)), &
      &       ity2mag(sky(4)) - dat(i) % sky
    end do iterate_datpoints

    deallocate(src)
  end block printcalib

  deallocate(mapi, maph)

end program sbverify

max_stars_repo_name: gronki/slezanbear

max_issues_repo_name: gronki/slezanbear

max_forks_repo_name: gronki/slezanbear


text:
!***********************************************************************
!                                                                      *
      INTEGER FUNCTION INDTPI (ITYPE, I)
!                                                                      *
!   Written by Farid A Parpia             Last revision: 09 Oct 1992   *
!                                                                      *
!***********************************************************************
!...Translated by Pacific-Sierra Research 77to90  4.3E  14:04:58   1/ 3/07
!...Modified by Charlotte Froese Fischer
!                     Gediminas Gaigalas  10/05/17
      USE bilst_C
!-----------------------------------------------
!   I n t e r f a c e   B l o c k s
!-----------------------------------------------
      IMPLICIT NONE
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      INTEGER, INTENT(IN) :: ITYPE
      INTEGER  :: I
!-----------------------------------------------
!
      SELECT CASE (ITYPE)
      CASE (1)
         INDTPI = INDTP1(I)
      CASE (2)
         INDTPI = INDTP2(I)
      CASE (3)
         INDTPI = INDTP3(I)
      CASE (4)
         INDTPI = INDTP4(I)
      CASE (5)
         INDTPI = INDTP5(I)
      CASE (6)
         INDTPI = INDTP6(I)
      END SELECT
!
      RETURN
      END FUNCTION INDTPI

max_stars_repo_name: sylas/grasp-continuum

max_issues_repo_name: compas/grasp-rci-qed

max_forks_repo_name: compas/grasp-rci-qed


text:
      INTEGER FUNCTION fortranwrite2 (IUNIT, STRING1, STRING2)
C
C     Interface routine for C/C++ functions to write to a file
C     under FORTRAN
!      use ifcore
C
      INTEGER IUNIT
      CHARACTER STRING1*(*), STRING2*(*)
C
      INTEGER ISTAT
C
      WRITE (IUNIT, 20, IOSTAT=ISTAT) STRING1, STRING2
 20   FORMAT(A,A)
      fortranwrite2 = ISTAT
!      ISTAT = COMMITQQ(IUNIT)
C
      RETURN
      END

ext: f

max_stars_repo_name: HydrologicEngineeringCenter/heclib

max_issues_repo_name: HydrologicEngineeringCenter/heclib

max_forks_repo_name: HydrologicEngineeringCenter/heclib


module readtree
  implicit none
  public
  integer :: nsteps !> Number of steps in the simulation

  integer, dimension(:), allocatable       :: nb_of_halos    !> halos in the tree
  integer, dimension(:), allocatable       :: nb_of_subhalos !> subhalos in the tree
  real (kind=4), dimension(:), allocatable :: aexp     !> expansion factor at each step
  real (kind=4), dimension(:), allocatable :: age_univ !> age of the universe at each step
  real (kind=4), dimension(:), allocatable :: omega_t  !> omega_t at at each step

  ! Internal variables
  integer :: istep = 1
  integer :: j = 1
  integer :: unit

  ! Some temp variable
  integer, allocatable       :: tmp_id_fathers(:), tmp_id_sons(:)
  real(kind=4), allocatable  :: tmp_m_fathers(:)

  private :: istep, j, unit, tmp_id_fathers, tmp_id_sons, tmp_m_fathers
contains
  !> Read a merger tree, storing the data about it in the module.
  subroutine read_tree(filename)

    character(len=*), intent(in) :: filename !> Tree filename

    open(newunit=unit, file=trim(filename), status='old', form='unformatted')
    read(unit) nsteps

    ! Deallocate all previously allocated data
    call free()

    allocate(nb_of_halos(nsteps), nb_of_subhalos(nsteps), aexp(nsteps), &
         omega_t(nsteps), age_univ(nsteps))
    read(unit) nb_of_halos, nb_of_subhalos
    read(unit) aexp
    read(unit) omega_t
    read(unit) age_univ

  end subroutine read_tree

  !> Read the next part of the tree. If the end is reached, set stop_now to True.
  !! To get the father / son of the halo, call directly after 'get_father' or 'get_son'
  subroutine iter_tree(stop_now, halo_id, BushID, step, level, hosthalo, hostsub, &
       nbsub, nextsub, nb_of_fathers, nb_of_sons, &
       m, macc, px, py, pz, vx, vy, vz, Lx, Ly, Lz, &
       r, ra, rb, rc, ek, ep, et, spin)

    logical, intent(out), optional         :: stop_now      !> true if end is reached
    integer(kind=4), intent(out), optional :: halo_id       !> id if the halo
    integer(kind=4), intent(out), optional :: BushID        !> ???
    integer(kind=4), intent(out), optional :: step          !> current step
    integer(kind=4), intent(out), optional :: level         !> level of the halo (1 is halo, 2+ are subhalos)
    integer(kind=4), intent(out), optional :: hosthalo      !> id of hosting halo
    integer(kind=4), intent(out), optional :: hostsub       !> ???
    integer(kind=4), intent(out), optional :: nbsub         !> number of subhalos
    integer(kind=4), intent(out), optional :: nextsub       !> ???
    integer(kind=4), intent(out), optional :: nb_of_fathers
    integer(kind=4), intent(out), optional :: nb_of_sons
    real(kind=4), intent(out), optional    :: m             !> Mass of the halo
    real(kind=4), intent(out), optional    :: macc          !> Accreted mass
    real(kind=4), intent(out), optional    :: px, py, pz    !> position of the halo
    real(kind=4), intent(out), optional    :: vx, vy, vz    !> velocitiy of the halo
    real(kind=4), intent(out), optional    :: Lx, Ly, Lz    !> angular momentum of the halo
    real(kind=4), intent(out), optional    :: r, ra, rb, rc !> ellispoid description
    real(kind=4), intent(out), optional    :: ek, ep, et    !> Kinetic, potential and thermic energy
    real(kind=4), intent(out), optional    :: spin          !> spin of the halo

    integer :: ntot

    stop_now = .false.
    ! Loop
    j = j + 1
    ntot = nb_of_subhalos(istep) + nb_of_halos(istep)

    if (j > ntot) then
       j = 1
       istep = istep + 1
       if (istep > nsteps) then
          call close_file()
          stop_now = .true.
          return
       end if
    end if

    ! Check that the halo_id is non negative (if it is, stop immediately)
    read(unit) halo_id
    if (halo_id == -1) then
       call close_file()
       stop_now = .true.
       return
    end if
    read(unit) BushID
    read(unit) step
    read(unit) level, hosthalo, hostsub, nbsub, nextsub
    read(unit) m
    read(unit) macc
    read(unit) px, py, pz
    read(unit) vx, vy, vz
    read(unit) Lx, Ly, Lz
    read(unit) r, ra, rb, rc
    read(unit) ek, ep, et
    read(unit) spin

    read(unit) nb_of_fathers

    if (nb_of_fathers > 0) then
       if (allocated(tmp_id_fathers)) deallocate(tmp_id_fathers)
       if (allocated(tmp_m_fathers )) deallocate(tmp_m_fathers)

       allocate(tmp_id_fathers(nb_of_fathers), tmp_m_fathers(nb_of_fathers))

       read(unit) tmp_id_fathers
       read(unit) tmp_m_fathers
    endif

    read(unit) nb_of_sons
    if(nb_of_sons>0)then
       if (allocated(tmp_id_sons)) then
          deallocate(tmp_id_sons)
       end if
       allocate(tmp_id_sons(nb_of_sons))

       read(unit) tmp_id_sons
    endif

    read(unit)
    read(unit)
    ! if(star) read(unit)
  end subroutine iter_tree

  !> Get the father of the last red halo
  subroutine get_fathers(nb_of_fathers, id_fathers, m_fathers)
    integer, intent(in)                                 :: nb_of_fathers !> number of son of the halo
    integer, dimension(nb_of_fathers), intent(out)      :: id_fathers     !> id of the father
    real(kind=4), dimension(nb_of_fathers), intent(out) :: m_fathers      !> mass of the father

    id_fathers = tmp_id_fathers
    m_fathers = tmp_m_fathers
  end subroutine get_fathers

  !> Get the son of the last red halo
  subroutine get_sons(nb_of_sons, id_sons)
    integer, intent(in)                         :: nb_of_sons!> number of son of the halo
    integer, dimension(nb_of_sons), intent(out) :: id_sons    !> ids of the son of the halo

    id_sons = tmp_id_sons
  end subroutine get_sons

  !> Close the tree file
  subroutine close_file()
    logical :: op
    ! TODO: inquire it is opened
    inquire(unit, opened=op)
    if (op) close(unit)
  end subroutine close_file

  !> Free allocated variables
  subroutine free()
    if (allocated(nb_of_halos    )) deallocate(nb_of_halos)
    if (allocated(nb_of_subhalos )) deallocate(nb_of_subhalos)
    if (allocated(aexp           )) deallocate(aexp)
    if (allocated(omega_t        )) deallocate(omega_t)
    if (allocated(age_univ       )) deallocate(age_univ)
  end subroutine free

  !> Revert the tree (from end to begin)
  subroutine revert_steps(filetree, outfich, star)
    character(len=*), intent(in) :: outfich, filetree
    logical, intent(in)            :: star

    integer                                    :: unit_in, unit_out
    integer                                    :: ntot, i, j, istep, idone
    integer(kind=4)                            :: nsteps, nbodies, nb_of_fathers, nb_of_sons, ncont
    integer(kind=4)                            :: mynumber, BushID, mystep, level, hosthalo, hostsub, nbsub, nextsub
    real(kind=4)                               ::m, macc, px, py, pz, vx, vy, vz, Lx, Ly, Lz, r, ra, rb, rc, ek, ep, et, spin
    real(kind=4)                               ::rvir, mvir, tvir, cvel, rho0, rc2
    integer(kind=4), allocatable, dimension(:) :: nb_of_halos, nb_of_subhalos
    integer(kind=4), allocatable, dimension(:) :: idfather, idson
    real(kind=4),    allocatable, dimension(:) :: aexp, age_univ, omega_t, mfather
    integer(kind=4), allocatable, dimension(:) :: nb_of_halos_rev, nb_of_subhalos_rev
    real(kind=4),    allocatable, dimension(:) :: aexp_rev, age_univ_rev, omega_t_rev

    !-----------------------------------------------
    ! Lecture du fichier contenant les galaxies
    !-----------------------------------------------
    write(*, *) filetree

    open(newunit=unit_in, file=filetree, status='old', form='unformatted')
    read(unit_in) nsteps
    allocate(nb_of_halos(1:nsteps), nb_of_subhalos(1:nsteps), aexp(1:nsteps), omega_t(1:nsteps), age_univ(1:nsteps))

    read(unit_in) nb_of_halos, nb_of_subhalos
    read(unit_in) aexp
    read(unit_in) omega_t
    read(unit_in) age_univ
    write(*, *) nsteps, ' steps'
    write(*, *) 'list of timesteps:'
    do i=1, nsteps
       write(*, *) 1d0/aexp(i)-1d0, age_univ(i)
    enddo
    write(*, *) MINVAL(nb_of_halos), MAXVAL(nb_of_halos), ' galaxies'
    write(*, *) MINVAL(nb_of_subhalos), MAXVAL(nb_of_subhalos), ' sub-galaxies'
    close(unit_in)

    allocate(nb_of_halos_rev(1:nsteps), nb_of_subhalos_rev(1:nsteps))
    allocate(aexp_rev(1:nsteps), omega_t_rev(1:nsteps), age_univ_rev(1:nsteps))
    do i=1, nsteps
       nb_of_halos_rev(nsteps-i+1)    = nb_of_halos(i)
       nb_of_subhalos_rev(nsteps-i+1) = nb_of_subhalos(i)
       aexp_rev(nsteps-i+1)     = aexp(i)
       omega_t_rev(nsteps-i+1)  = omega_t(i)
       age_univ_rev(nsteps-i+1) = age_univ(i)
    enddo

    open(newunit=unit_out, file=trim(outfich), form='unformatted')
    write(unit_out) nsteps
    write(unit_out) nb_of_halos_rev, nb_of_subhalos_rev
    write(unit_out) aexp_rev
    write(unit_out) omega_t_rev
    write(unit_out) age_univ_rev
    deallocate(nb_of_halos_rev, nb_of_subhalos_rev, aexp_rev, omega_t_rev, age_univ_rev)

    idone = nsteps

    do while(idone > 0)

       open(newunit=unit_in, file=trim(filetree), status='old', form='unformatted')
       read(unit_in) nsteps
       read(unit_in) nb_of_halos, nb_of_subhalos
       read(unit_in) aexp
       read(unit_in) omega_t
       read(unit_in) age_univ

       do istep=1, idone

          ntot = nb_of_halos(istep) + nb_of_subhalos(istep)

          if (istep == idone) then

             do j = 1, ntot
                read(unit_in) mynumber
                read(unit_in) BushID
                read(unit_in) mystep
                read(unit_in) level, hosthalo, hostsub, nbsub, nextsub
                read(unit_in) m
                read(unit_in) macc
                read(unit_in) px, py, pz
                read(unit_in) vx, vy, vz
                read(unit_in) Lx, Ly, Lz
                read(unit_in) r, ra, rb, rc
                read(unit_in) ek, ep, et
                read(unit_in) spin

                write(unit_out) mynumber
                write(unit_out) BushID
                write(unit_out) mystep
                write(unit_out) level, hosthalo, hostsub, nbsub, nextsub
                write(unit_out) m
                write(unit_out) macc
                write(unit_out) px, py, pz
                write(unit_out) vx, vy, vz
                write(unit_out) Lx, Ly, Lz
                write(unit_out) r, ra, rb, rc
                write(unit_out) ek, ep, et
                write(unit_out) spin

                read(unit_in) nb_of_fathers
                write(unit_out) nb_of_fathers
                if (nb_of_fathers > 0) then
                   allocate(idfather(nb_of_fathers), mfather(nb_of_fathers))
                   read(unit_in)   idfather
                   read(unit_in)   mfather

                   write(unit_out) idfather
                   write(unit_out) mfather
                   deallocate(idfather, mfather)
                endif

                read(unit_in) nb_of_sons
                write(unit_out) nb_of_sons

                if (nb_of_sons > 0) then
                   allocate(idson(nb_of_sons))
                   read(unit_in)   idson

                   write(unit_out) idson
                   deallocate(idson)
                endif

                read(unit_in) rvir, mvir, tvir, cvel
                read(unit_in) rho0, rc2
                if (star) read(unit_in) ncont

                write(unit_out) rvir, mvir, tvir, cvel
                write(unit_out) rho0, rc2
                if (star) write(unit_out) ncont
             enddo

             write(*, *)'wrote ', idone, 1d0/aexp(idone)-1d0
             idone = idone-1

          else if (istep < idone) then

             do j=1, ntot
                read(unit_in)
                read(unit_in)
                read(unit_in)
                read(unit_in)
                read(unit_in)
                read(unit_in)
                read(unit_in)
                read(unit_in)
                read(unit_in)
                read(unit_in)
                read(unit_in)
                read(unit_in)

                read(unit_in) nb_of_fathers
                if (nb_of_fathers > 0) then
                   allocate(idfather(1:nb_of_fathers), mfather(1:nb_of_fathers))
                   read(unit_in) idfather
                   read(unit_in) mfather
                   deallocate(idfather, mfather)
                endif

                read(unit_in) nb_of_sons
                if (nb_of_sons > 0) then
                   allocate(idson(1:nb_of_sons))
                   read(unit_in) idson
                   deallocate(idson)
                endif

                read(unit_in)
                read(unit_in)
                if (star) read(unit_in)
             enddo

          endif

       enddo
       close(unit_in)

    enddo
    close(unit_out)

  end subroutine revert_steps
end module readtree


text:
      program hycom_month_day
      implicit none
c
c     input:  ddd (ordinal day)
c     output: ddd, months that span the day, and their weights (linear interpolation)
c
c     example: echo 360 | hycom_month_day
c     result:  360 12 01 0.6967211 0.3032789
c
      integer jday,mon0,mon1
      real*4  w0,w1,x
c
      read(5,*) jday
      x    = 1.0 + mod( jday+366.0-15.25,366.0)/30.5
      mon0 = x
      mon1 = mod(mon0,12)+1
      w1   = mod(x,1.0) !fraction
      w0   = 1.0-w1
      write(6,'(i3.3,2i3.2,2f10.7)') jday,mon0,mon1,w0,w1
      call exit(0)
      end

max_stars_repo_name: TillRasmussen/HYCOM-tools

max_issues_repo_name: TillRasmussen/HYCOM-tools

max_forks_repo_name: TillRasmussen/HYCOM-tools


text:
  subroutine clearsky_optical_properties( &
       wavelength_ns, nsmax,      ns_elas,    &
       maxlay,        nstokes,    maxleg,     &
       relN2,         relO2,      relAr,      &
       csabs_O3,      csabs_O4,   csabs_NO2,  &
       csRay_N2,      csRay_O2,   csRay_Ar,   &
       csram_N2,      csram_O2,   csram_air,  &
       eps_N2,        eps_O2,     eps_Ar,     &
       dvO4,          dvO3,       dvNO2,      &
       dvair,         dz,         tautot,     &
       omega_Ray,     omega_Cab,  om_air,     &
       pRay_N2,       pCab_N2,    pRay_O2,    &
       pCab_O2,       pRay_Ar,    pCab_Ar,    &
       pRay_air,      pCab_air,   pRam_air,   &
       maxcoeffs)

!-------------------------------------------------------------------
!
! Subroutine:
!   clearsky_optical_properties
!
! Purpose:
!   Return optical properties of clear sky atmosphere
!   (optical thickness; single scattering albedo; phase functions;
!    expansion coefficient matrices)
!!
! Assumption: phase matrix does not change over wavelength range
! wavelength_ns(0 ... nsmax)
!
! omega = (csram+cscab) / csray
!
! For phase matrix and phase function see Table A.1 in
! appendix A of "Rotational Raman scattering  of polarized light in
! the Earth's atmosphere..." - Landgraf et al., jqsrt 2004
!-------------------------------------------------------------------
!
! The Rayleigh, Cabannes and Raman phase matrix P has the shape:
!
!        ( P11  P12   0    0  )
!   P =  ( P21  P22   0    0  )   where P11 = phase function
!        (  0    0   P33  P34 )
!        (  0    0   P43  P44 )
!
!        P11 = [ A [1 + (costheta)**2] + C] / N
!        P12 =  -A [1 - (costheta)**2] + C] / N
!        P22 =   A [1 + (costheta)**2] / N
!        P33 =  2A [costheta] / N
!        P44 =   B [costheta] / N
!        P21 =  P12
!        P43 =  -P34 = 0
!
!  Rayleigh:
!
!    A = 3(45 + eps)
!    B = 30(9 - eps)
!    C = 36 eps
!    N = 180 + 40 eps
!
!  Cabannes:
!
!    A = 3(180 + eps)
!    B = 30(36 - eps)
!    C = 36 eps
!    N = 40(18 + eps)
!
!  Raman:
!    A =   3
!    B =  -30 (erratum Landgraf et al. -10 -> -30)
!    C =  36
!    N =  40
!
!  where eps is the molecular anistropy,
!  related to the king factor (Bates 1984) as Fking = 1 + (2/9)eps
!
!  The phase matrix is expanded in generalized spherical function as
!  follows:
!                 2
!    P11       = sum [alpha1(ns)] * p(l,0,0)
!                l=0
!
!                 2
!    P22 + P33 = sum [alpha2(ns) + alpha3(ns)] * p(l,2,+2)
!                l=2
!
!                 2
!    P22 - P33 = sum [alpha2(ns) - alpha3(ns)] * p(l,2,-2)
!                l=2
!
!                 2
!    P44       = sum [alpha4(ns)] * p(l,0,0)
!                l=0
!
!                 2
!    P12 = P21 = sum [beta1(ns)] * p(l,0,+2)
!                l=2
!
!                 2
!    P34       = sum [beta2(ns)] * p(l,0,+2)
!                l=2
!
!    P43       = -P34
!
!  where generalized spherical functions p(l,m,n) are given by
!  [de Haan et al 1987, A&A 183, 371-391]:
!
!    l    p(l,m,n)
!   ---------------------------------------------------------
!    0    p(0,0,0)  = 1
!    1    p(1,0,0)  = costheta
!    2    p(2,0,0)  = 0.50*[3(costheta)**2 - 1]
!    0    p(0,0,+2) = p(0,0,-2) =  0
!    1    p(1,0,+2) = p(1,0,-2) =  0
!    2    p(2,0,+2) = p(2,0,-2) = -0.25*sqrt(6)*[1 - (costheta)**2]
!    0    p(0,2,+2) = p(0,2,-2) = 0
!    1    p(1,2,+2) = p(0,2,-2) = 0
!    2    p(2,2,+2) = 0.25*[1 + 2costheta + costheta**2]
!    2    p(2,2,-2) = 0.25*[1 - 2costheta + costheta**2]
!
!  see also subroutine "gsf.f90"
!
!  Expansion coefficient matrix for Rayleigh, Cabannes and Raman
!  scattering is returned as phaseRay, phaseCab and phaseRam :
!
!    ( alpha1  beta1     0       0     )
!    ( beta1   alpha2    0       0     )
!    (   0       0     alpha3   beta2  )
!    (   0       0     -beta2   alpha4 )
!
!-------------------------------------------------------------------

  implicit none

  integer,parameter::dp=selected_real_kind(8) !double precision

  integer,intent(in):: &
    nsmax,    & !nr of wavenumbers (incident light)
    ns_elas,  & !wavelength_ns elastic scattering
    maxlay,   &
    nstokes,  &
    maxleg

  real(kind=dp),intent(in):: &
    wavelength_ns(0:nsmax)

  real(kind=dp),intent(in):: &
     relN2, & !N2 volume mixing ratio
     relO2, & !O2 volume mixing ratio
     relAr    !Ar volume mixing ratio

  real(kind=dp),intent(in):: &
    dz(maxlay),              & !layer tickness [m]
    csabs_O3(maxlay,0:nsmax),  & !O3 csabs section [cm^2]
    csabs_O4(maxlay,0:nsmax),  & !O4 csabs section [cm^5 molec-2]
    csabs_NO2(maxlay,0:nsmax), & !NO2 csabs section [cm^2]
    csRay_N2(0:nsmax),         & !Rayleigh csabs section [cm2] for N2
    csRay_O2(0:nsmax),         & !Rayleigh csabs section [cm2] for O2
    csRay_Ar(0:nsmax),         & !Rayleigh csabs section [cm2] for Ar
    eps_N2(0:nsmax),           & !epsilon for N2
    eps_O2(0:nsmax),           & !epsilon for O2
    eps_Ar(0:nsmax),           & !epsilon for Ar = 0
    dvO4(maxlay),            & !diff. O4 column density [part^2./cm^5]
    dvO3(maxlay),            & !diff. O3 column density [part./cm^2]
    dvNO2(maxlay),           & !diff. NO2 column density [part./cm^2]
    dvair(maxlay)              !diff. O2 column density [part./cm^2]

  real(kind=dp),intent(in):: &
    csram_N2(maxlay,0:nsmax,0:nsmax), &
    csram_O2(maxlay,0:nsmax,0:nsmax), &
    csram_air(maxlay,0:nsmax,0:nsmax) !Cab+Ram

  real(kind=dp)::  &
    elasfac_N2,       & !elastic fraction
    elasfac_O2,       &
    elasfac_Ar,       &
    weight_N2,     &
    weight_O2,     &
    weight_Ar,     &
    tausca_Ray_N2, &
    tausca_Ray_O2, &
    tausca_Ray_Ar

! phase matrix for vector radiative transfer,
! elements are divided by 2m+1
  real(kind=dp):: &
    A_Ray_N2, & !Rayleigh
    A_Ray_O2, &
    A_Ray_Ar, &
    B_Ray_N2, &
    B_Ray_O2, &
    B_Ray_Ar, &
    C_Ray_N2, &
    C_Ray_O2, &
    C_Ray_Ar, &
    N_Ray_N2, &
    N_Ray_O2, &
    N_Ray_Ar, &
    A_Cab_N2, & !Cabannes
    A_Cab_O2, &
    A_Cab_Ar, &
    B_Cab_N2, &
    B_Cab_O2, &
    B_Cab_Ar, &
    C_Cab_N2, &
    C_Cab_O2, &
    C_Cab_Ar, &
    N_Cab_N2, &
    N_Cab_O2, &
    N_Cab_Ar, &
    A_Ram,    & !Raman
    B_Ram,    &
    C_Ram,    &
    N_Ram

  real(kind=dp):: &
    phaseRay_air(4,4,0:2), & !Rayleigh phase fct. moments
    phaseRay_N2(4,4,0:2),  & !Rayleigh phase fct. moments
    phaseRay_O2(4,4,0:2),  & !Rayleigh phase fct. moments
    phaseRay_Ar(4,4,0:2),  & !Rayleigh phase fct. moments
    phaseCab_air(4,4,0:2), & !Cabannes phase fct. moments
    phaseCab_N2(4,4,0:2),  & !Cabannes phase fct. moments
    phaseCab_O2(4,4,0:2),  & !Cabannes phase fct. moments
    phaseCab_Ar(4,4,0:2),  & !Cabannes phase fct. moments
    phaseRam(4,4,0:2)        !Raman phase fct. moments

  integer::ns,k,m,lg,i1,i2,i_st,j_st

  real(kind=dp),dimension(maxlay,0:nsmax):: &
    tauabs, & !total abs. optical depth
    tausca,&    !total scat. optical depth
    bscat_Ray, & !beta_scat [m-1]
    bscat_Cab  !beta_scat [m-1]

  real(kind=dp),dimension(maxlay,0:nsmax),intent(out):: &
    tautot,    & !total optical depth
    omega_Ray, & !Rayleigh single scattering albedo
    omega_Cab    !Cabannes single scattering albedo

  real(kind=dp),intent(out):: &
    om_air(maxlay,0:nsmax,0:nsmax) !watch out, different for each layer k

 real(kind=dp), &
  dimension(nstokes,nstokes,0:maxleg,maxlay),intent(out):: &
    pCab_N2,  & !Cabannes phase fct. moments
    pRay_N2,  & !Rayleigh phase fct. moments
    pCab_O2,  & !Cabannes phase fct. moments
    pRay_O2,  & !Rayleigh phase fct. moments
    pCab_Ar,  & !Cabannes phase fct. moments
    pRay_Ar,  & !Rayleigh phase fct. moments
    pCab_air, & !Cabannes phase fct. moments
    pRay_air, & !Rayleigh phase fct. moments
    pRam_air    !Raman    phase fct. moments

 integer::maxcoeffs(maxlay) !number of phase fct. coeff.

 logical,parameter::scalar =.true.

 logical,parameter:: &
   verbose_input  = .false., & !show setup expansion coefficients
   verbose_output = .false.    !show selected expans.coefficients

 do ns=0,nsmax

!-------------------------------------------------------------------
!    Optical depths and single scattering albedos
!-------------------------------------------------------------------

     elasfac_N2 = (18.d0+eps_N2(ns))/(18.d0+4.d0*eps_N2(ns))
     elasfac_O2 = (18.d0+eps_O2(ns))/(18.d0+4.d0*eps_O2(ns))
     elasfac_Ar =  1.d0

!    Total optical depth (Rayleigh + absorption) and
!    single scattering albedo for Cabannes and Rayleigh scattering

     do k = 1,maxlay

        tausca_Ray_N2 = csRay_N2(ns) * relN2 * dvair(k)
        tausca_Ray_O2 = csRay_O2(ns) * relO2 * dvair(k)
        tausca_Ray_Ar = csRay_Ar(ns) * relAr * dvair(k)

        tausca(k,ns) = tausca_Ray_N2 + &
                       tausca_Ray_O2 + &
                       tausca_Ray_Ar

        tauabs(k,ns) = csabs_O3(k,ns)  * dvO3(k)  + &
                       csabs_O4(k,ns)  * dvO4(k)  + &
                       csabs_NO2(k,ns) * dvNO2(k)

        tautot(k,ns) = tauabs(k,ns) + tausca(k,ns)

!        if(ns.eq.ns_elas) &
!           write(*,'(i4,f10.4,4d12.4)') &
!             k,wavelength_ns(ns),tausca(k,ns_elas),tauabs(k,ns_elas),tautot(k,ns_elas)

        omega_Ray(k,ns)    = 0.d0
        omega_Cab(k,ns)    = 0.d0
        om_air(k,0:nsmax,ns) = 0.d0

        if(tautot(k,ns).ne.0.d0)then

           omega_Ray(k,ns) = tausca(k,ns)/tautot(k,ns)
           omega_Cab(k,ns) =(elasfac_N2*tausca_Ray_N2 + &
                             elasfac_O2*tausca_Ray_O2 + &
                             elasfac_Ar*tausca_Ray_Ar)/tautot(k,ns)

           om_air(k,0:nsmax,ns) = &
             csram_air(k,0:nsmax,ns)*dvair(k)/tautot(k,ns)
        endif

        bscat_Ray(k,ns) = omega_Ray(k,ns)*tautot(k,ns)/dz(k) ![m-1]
        bscat_Cab(k,ns) = omega_Cab(k,ns)*tautot(k,ns)/dz(k) ![m-1]

     enddo

!--------------------------------------------------------------------
!    Setup Rayleigh phase matrix,
!    assume independent of ns (calculate for ns_elas)
!-------------------------------------------------------------------

     if(ns.eq.ns_elas)then

!       Save omega, tautot and pRay for selected wavelength
!       weight = omega_Ray_species(k)/omega_Ray_air(k,ns),
!       which does not depend on altitude

        weight_N2 = csRay_N2(ns)*relN2/(csRay_N2(ns)*relN2 + &
                                        csRay_O2(ns)*relO2 + &
                                        csRay_Ar(ns)*relAr)
        weight_O2 = csRay_O2(ns)*relO2/(csRay_N2(ns)*relN2 + &
                                        csRay_O2(ns)*relO2 + &
                                        csRay_Ar(ns)*relAr)
        weight_Ar = csRay_Ar(ns)*relAr/(csRay_N2(ns)*relN2 + &
                                        csRay_O2(ns)*relO2 + &
                                        csRay_Ar(ns)*relAr)

!       parameters A,B,C and normalization constant N in
!       Table A.1 landgraf et al.

        A_Ray_N2 = 3.d0 * (45.d0 + eps_N2(ns))
        A_Ray_O2 = 3.d0 * (45.d0 + eps_O2(ns))
        A_Ray_Ar = 3.d0 *  45.d0

        B_Ray_N2 = 30.d0 * (9.d0 - eps_N2(ns))
        B_Ray_O2 = 30.d0 * (9.d0 - eps_O2(ns))
        B_Ray_Ar = 30.d0 *  9.d0

        C_Ray_N2 = 36.d0 * eps_N2(ns)
        C_Ray_O2 = 36.d0 * eps_O2(ns)
        C_Ray_Ar =  0.d0

        N_Ray_N2 = 180.d0 + 40.d0 * eps_N2(ns)
        N_Ray_O2 = 180.d0 + 40.d0 * eps_O2(ns)
        N_Ray_Ar = 180.d0 + 40.d0 * eps_Ar(ns)

!       Rayleigh phase matrix elements
!       in generalized spherical function expansion

        call RayCaB_Ram_EXPANSION_coeff( &
             A_Ray_N2, B_Ray_N2, C_Ray_N2, N_Ray_N2, phaseRay_N2)

        call RayCaB_Ram_EXPANSION_coeff( &
             A_Ray_O2, B_Ray_O2, C_Ray_O2, N_Ray_O2, phaseRay_O2)

        call RayCaB_Ram_EXPANSION_coeff( &
             A_Ray_Ar, B_Ray_Ar, C_Ray_Ar, N_Ray_Ar, phaseRay_Ar)

        do m=0,2
           do i2=1,4
              do i1=1,4
                 phaseRay_air(i1,i2,m) = &
                      weight_N2 * phaseRay_N2(i1,i2,m) + &
                      weight_O2 * phaseRay_O2(i1,i2,m) + &
                      weight_Ar * phaseRay_Ar(i1,i2,m)
              enddo
           enddo
        enddo

     endif !ns_elas

!--------------------------------------------------------------------
!    Setup Cabannes phase matrix, assume independent of ns
!    (Calculate for ns_elas)
!--------------------------------------------------------------------

     if(ns.eq.ns_elas)then

!       weight = omega_Cab_species(k)/omega_Cab(k,ns)
!       which does not depend on altitude

        weight_N2 = elasfac_N2 * csRay_N2(ns) * relN2 / &
                   (elasfac_N2 * csRay_N2(ns) * relN2 + &
                    elasfac_O2 * csRay_O2(ns) * relO2 + &
                    elasfac_Ar * csRay_Ar(ns) * relAr)
        weight_O2 = elasfac_O2 * csRay_O2(ns) * relO2 / &
                   (elasfac_N2 * csRay_N2(ns) * relN2 + &
                    elasfac_O2 * csRay_O2(ns) * relO2 + &
                    elasfac_Ar * csRay_Ar(ns) * relAr)
        weight_Ar = elasfac_Ar * csRay_Ar(ns) * relAr / &
                   (elasfac_N2 * csRay_N2(ns) * relN2 + &
                    elasfac_O2 * csRay_O2(ns) * relO2 + &
                    elasfac_Ar * csRay_Ar(ns) * relAr)

!       parameters A,B,C and normalization constant N in
!       Table A.1 landgraf et al.

        A_Cab_N2 = 3.d0 * (180.d0 + eps_N2(ns))
        A_Cab_O2 = 3.d0 * (180.d0 + eps_O2(ns))
        A_Cab_Ar = 3.d0 *  180.d0

        B_Cab_N2 = 30.d0 * (36.d0 - eps_N2(ns))
        B_Cab_O2 = 30.d0 * (36.d0 - eps_O2(ns))
        B_Cab_Ar = 30.d0 *  36.d0

        C_Cab_N2 = 36.d0 * eps_N2(ns)
        C_Cab_O2 = 36.d0 * eps_O2(ns)
        C_Cab_Ar =  0.d0

        N_Cab_N2 = 40.d0 * (18.d0 + eps_N2(ns))
        N_Cab_O2 = 40.d0 * (18.d0 + eps_O2(ns))
        N_Cab_Ar = 40.d0 *  18.d0

!       Cabannes phase matrix elements
!       in generalized spherical function expansion

        call RayCaB_Ram_EXPANSION_coeff( &
             A_Cab_N2, B_Cab_N2, C_Cab_N2, N_Cab_N2, phaseCab_N2)

        call RayCaB_Ram_EXPANSION_coeff( &
             A_Cab_O2, B_Cab_O2, C_Cab_O2, N_Cab_O2, phaseCab_O2)

        call RayCaB_Ram_EXPANSION_coeff( &
             A_Cab_Ar, B_Cab_Ar, C_Cab_Ar, N_Cab_Ar, phaseCab_Ar)

        do m=0,2
           do i2=1,4
              do i1=1,4
                 phaseCab_air(i1,i2,m) = &
                      weight_N2 * phaseCab_N2(i1,i2,m) + &
                      weight_O2 * phaseCab_O2(i1,i2,m) + &
                      weight_Ar * phaseCab_Ar(i1,i2,m)
              enddo
           enddo
        enddo

     endif !ns_elas only

  enddo !ns

!--------------------------------------------------------------------
! Raman phase matrix, phase function
! independent of wavelength and species
!--------------------------------------------------------------------

! parameters A,B,C and normalization constant N in
! Table A.1 landgraf et al.

  A_Ram =   3.d0
  B_Ram = -30.d0 !correction, used to be -10.d0 in paper Landgraf 2004
  C_Ram =  36.d0
  N_Ram =  40.d0

  call RayCaB_Ram_EXPANSION_coeff( &
       A_Ram, B_Ram, C_Ram, N_Ram, phaseRam)

  if(verbose_input)then
     write(*,*) 'Setup Expansion coefficients for Rayleigh scattering'
     write(*,'(A4,6A9)') &
       'lg','alpha1','alpha2','alpha3','alpha4','beta1','beta2'
     do lg=0,2
        write(*,'(I4,6F9.5)')   &
          lg,                   &
          phaseRay_air(1,1,lg), &
          phaseRay_air(2,2,lg), &
          phaseRay_air(3,3,lg), &
          phaseRay_air(4,4,lg), &
          phaseRay_air(1,2,lg), &
          phaseRay_air(3,4,lg)
      enddo
      write(*,*) &
        'Setup Expansion coefficients for Cabannes scattering'
      write(*,'(A4,6A9)') &
        'lg','alpha1','alpha2','alpha3','alpha4','beta1','beta2'
      do lg=0,2
         write(*,'(I4,6F9.5)')   &
           lg,                   &
           phaseCab_air(1,1,lg), &
           phaseCab_air(2,2,lg), &
           phaseCab_air(3,3,lg), &
           phaseCab_air(4,4,lg), &
           phaseCab_air(1,2,lg), &
           phaseCab_air(3,4,lg)
       enddo
       write(*,*) &
         'Setup Expansion coefficients for Raman scattering'
       write(*,'(A4,6A9)') &
         'lg','alpha1','alpha2','alpha3','alpha4','beta1','beta2'
       do lg=0,2
          write(*,'(I4,6F9.5)')  &
            lg,                  &
            phaseRam(1,1,lg),    &
            phaseRam(2,2,lg),    &
            phaseRam(3,3,lg),    &
            phaseRam(4,4,lg),    &
            phaseRam(1,2,lg),    &
            phaseRam(3,4,lg)
       enddo
  endif

!--------------------------------------------------------------------
! Select phase matrix elements
! for (4,4,lgmax) -> (nstokes,nstokes,maxleg)
!--------------------------------------------------------------------

  do k = 1,maxlay
     maxcoeffs(k) = 2
  enddo

  do k=1,maxlay
   do lg = 0,maxleg
    do i_st = 1,nstokes
     do j_st = 1,nstokes
      if(lg.lE.2)then
         pRay_N2( j_st,i_st,lg,k) = phaseRay_N2( j_st,i_st,lg)
         pCab_N2( j_st,i_st,lg,k) = phaseCab_N2( j_st,i_st,lg)
         pRay_O2( j_st,i_st,lg,k) = phaseRay_O2( j_st,i_st,lg)
         pCab_O2( j_st,i_st,lg,k) = phaseCab_O2( j_st,i_st,lg)
         pRay_Ar( j_st,i_st,lg,k) = phaseRay_Ar( j_st,i_st,lg)
         pCab_Ar( j_st,i_st,lg,k) = phaseCab_Ar( j_st,i_st,lg)
         pRay_air(j_st,i_st,lg,k) = phaseRay_air(j_st,i_st,lg)
         pCab_air(j_st,i_st,lg,k) = phaseCab_air(j_st,i_st,lg)
         pRam_air(j_st,i_st,lg,k) = phaseRam(j_st,i_st,lg)
         if(scalar)then
            pRay_N2( j_st,i_st,lg,k) = 0.d0
            pCab_N2( j_st,i_st,lg,k) = 0.d0
            pRay_O2( j_st,i_st,lg,k) = 0.d0
            pCab_O2( j_st,i_st,lg,k) = 0.d0
            pRay_Ar( j_st,i_st,lg,k) = 0.d0
            pCab_Ar( j_st,i_st,lg,k) = 0.d0
            pRay_air(j_st,i_st,lg,k) = 0.d0
            pCab_air(j_st,i_st,lg,k) = 0.d0
            pRam_air(j_st,i_st,lg,k) = 0.d0
            if((j_st.eq.1).AND.(i_st.eq.1))then
               pRay_N2( j_st,i_st,lg,k) = phaseRay_N2( j_st,i_st,lg)
               pCab_N2( j_st,i_st,lg,k) = phaseCab_N2( j_st,i_st,lg)
               pRay_O2( j_st,i_st,lg,k) = phaseRay_O2( j_st,i_st,lg)
               pCab_O2( j_st,i_st,lg,k) = phaseCab_O2( j_st,i_st,lg)
               pRay_Ar( j_st,i_st,lg,k) = phaseRay_Ar( j_st,i_st,lg)
               pCab_Ar( j_st,i_st,lg,k) = phaseCab_Ar( j_st,i_st,lg)
               pRay_air(j_st,i_st,lg,k) = phaseRay_air(j_st,i_st,lg)
               pCab_air(j_st,i_st,lg,k) = phaseCab_air(j_st,i_st,lg)
               pRam_air(j_st,i_st,lg,k) = phaseRam(j_st,i_st,lg)
            endif
         endif
      else
          pRay_N2( j_st,i_st,lg,k) = 0.d0
          pCab_N2( j_st,i_st,lg,k) = 0.d0
          pRay_O2( j_st,i_st,lg,k) = 0.d0
          pCab_O2( j_st,i_st,lg,k) = 0.d0
          pRay_Ar( j_st,i_st,lg,k) = 0.d0
          pCab_Ar( j_st,i_st,lg,k) = 0.d0
          pRay_air(j_st,i_st,lg,k) = 0.d0
          pCab_air(j_st,i_st,lg,k) = 0.d0
          pRam_air(j_st,i_st,lg,k) = 0.d0
      endif
     enddo
    enddo
   enddo
  enddo

 goto 44
  write(*,*)'Isotropic phasefunctions'
  do lg = 0,maxleg
   do i_st = 1,nstokes
    do j_st = 1,nstokes
     do k = 1,maxlay
        if((i_st.eq.1).AND.(j_st.eq.1).AND.(lg.eq.0))then
           pRay_N2( j_st,i_st,lg,k) = 1.d0
           pCab_N2( j_st,i_st,lg,k) = 1.d0
           pRay_O2( j_st,i_st,lg,k) = 1.d0
           pCab_O2( j_st,i_st,lg,k) = 1.d0
           pRay_Ar( j_st,i_st,lg,k) = 1.d0
           pCab_Ar( j_st,i_st,lg,k) = 1.d0
           pRay_air(j_st,i_st,lg,k) = 1.d0
           pCab_air(j_st,i_st,lg,k) = 1.d0
        else
           pRay_N2( j_st,i_st,lg,k) = 0.d0
           pCab_N2( j_st,i_st,lg,k) = 0.d0
           pRay_O2( j_st,i_st,lg,k) = 0.d0
           pCab_O2( j_st,i_st,lg,k) = 0.d0
           pRay_Ar( j_st,i_st,lg,k) = 0.d0
           pCab_Ar( j_st,i_st,lg,k) = 0.d0
           pRay_air(j_st,i_st,lg,k) = 0.d0
           pCab_air(j_st,i_st,lg,k) = 0.d0
        endif
        if((i_st.eq.1).AND.(j_st.eq.1).AND.(lg.eq.0))then
           pRam_air(j_st,i_st,lg,k) = 1.d0
        else
          pRam_air(j_st,i_st,lg,k) = 0.d0
        endif
      enddo
     enddo
    enddo
   enddo
44 continue

  if(verbose_output)then
     write(*,*) &
       'Output Expansion coefficients for Rayleigh scattering'
     write(*,*) 'lambda [nm] = ',wavelength_ns(ns_elas)
     write(*,'(A4,6A9)') &
       'lg','alpha1','alpha2','alpha3','alpha4','beta1','beta2'
     k = maxlay
     do lg=0,maxleg
        if(nstokes.eq.1)then
           write(*,'(I4,6F9.5)')  &
             lg,                  &
             pRay_air(1,1,lg,k),  &
             0.d0,                &
             0.d0,                &
             0.d0,                &
             0.d0,                &
             0.d0
        endif
        if(nstokes.eq.3)then
           write(*,'(I4,6F9.5)')  &
             lg,                  &
             pRay_air(1,1,lg,k),  &
             pRay_air(2,2,lg,k),  &
             pRay_air(3,3,lg,k),  &
             0.d0,                &
             pRay_air(1,2,lg,k),  &
             0.d0
        endif
        if(nstokes.eq.4)then
           write(*,'(I4,6F9.5)')  &
             lg,                  &
             pRay_air(1,1,lg,k),  &
             pRay_air(2,2,lg,k),  &
             pRay_air(3,3,lg,k),  &
             pRay_air(4,4,lg,k),  &
             pRay_air(1,2,lg,k),  &
             pRay_air(3,4,lg,k)
        endif
     enddo
     write(*,*) &
       'Output Expansion coefficients for Cabannes scattering'
     write(*,'(A4,6A9)') &
       'lg','alpha1','alpha2','alpha3','alpha4','beta1','beta2'
     k = maxlay
     do lg=0,maxleg
        if(nstokes.eq.1)then
           write(*,'(I4,6F9.5)') &
             lg,                 &
             pCab_air(1,1,lg,k), &
             0.d0,               &
             0.d0,               &
             0.d0,               &
             0.d0,               &
             0.d0
         endif
         if(nstokes.eq.3)then
            write(*,'(I4,6F9.5)') &
              lg,                 &
              pCab_air(1,1,lg,k), &
              pCab_air(2,2,lg,k), &
              pCab_air(3,3,lg,k), &
              0.d0,               &
              pCab_air(1,2,lg,k), &
              0.d0
         endif
         if(nstokes.eq.4)then
            write(*,'(I4,6F9.5)') &
              lg,                 &
              pCab_air(1,1,lg,k), &
              pCab_air(2,2,lg,k), &
              pCab_air(3,3,lg,k), &
              pCab_air(4,4,lg,k), &
              pCab_air(1,2,lg,k), &
              pCab_air(3,4,lg,k)
         endif
      enddo
      write(*,*) &
        'Output Expansion coefficients for Raman scattering'
      write(*,*) 'lambda [nm] = ',wavelength_ns(ns_elas)
      write(*,'(A4,6A9)') &
        'lg','alpha1','alpha2','alpha3','alpha4','beta1','beta2'
      k = maxlay
      do lg=0,maxleg
         if(nstokes.eq.1)then
            write(*,'(I4,6F9.5)') &
              lg,                 &
              pRam_air(1,1,lg,k), &
              0.d0,               &
              0.d0,               &
              0.d0,               &
              0.d0,               &
              0.d0
         endif
         if(nstokes.eq.3)then
            write(*,'(I4,6F9.5)') &
              lg,                 &
              pRam_air(1,1,lg,k), &
              pRam_air(2,2,lg,k), &
              pRam_air(3,3,lg,k), &
              0.d0,               &
              pRam_air(1,2,lg,k), &
              0.d0
         endif
         if(nstokes.eq.4)then
            write(*,'(I4,6F9.5)') &
              lg,                 &
              pRam_air(1,1,lg,k), &
              pRam_air(2,2,lg,k), &
              pRam_air(3,3,lg,k), &
              pRam_air(4,4,lg,k), &
              pRam_air(1,2,lg,k), &
              pRam_air(3,4,lg,k)
         endif
      enddo
  endif

  return
end subroutine clearsky_optical_properties

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

subroutine RayCaB_Ram_expansion_coeff(A, B, C, N, coeff)

!--------------------------------------------------------------------
! Purpose:
!   Returns expansion coefficient matrix for Rayleigh, Cabannes
!   scattering in the form of weights,
!
!             ( alpha1  beta1     0       0     )
!     coeff = ( beta1   alpha2    0       0     )
!             (   0       0     alpha3   beta2  )
!             (   0       0     -beta2   alpha4 )
!
!   where expansion coeffcicients depend on Fourier moment lg=0,..,2
!   and wavelength [see Stam et al. JGR vol 107, D20, 2002]
!
!    lg   alpha1  alpha2  alpha3  alpha4      beta1      beta2
!   ------------------------------------------------------------
!     0      1       0       0       0          0           0
!     1      0       0       0      B/N         0           0
!     2    2A/3N    4A/N     0       0    4A/(N*(6**0.5))   0
!
!   the rest is zero. coeff=coeff(i1,i2,lg) where i1,i2 are index
!   for Stokes component, lg is Fourier moment index
!--------------------------------------------------------------------

  implicit none
  integer,parameter::dp=selected_real_kind(8) !double precision
  integer::lg,i1,i2
  logical,parameter::verbose = .false.
    real(kind=dp)::A,B,C,N,coeff(4,4,0:2)

  if(verbose) write(*,*) A,B,C,N

! Set all to zero
  do lg=0,2
     do i2=1,4
        do i1=1,4
           coeff(i1,i2,lg) = 0.d0
        enddo
     enddo
  enddo

! alpha1
  coeff(1,1,0) = 1.d0
  coeff(1,1,1) = 0.d0
  coeff(1,1,2) = 2.d0*A/(3.d0*N)

! alpha2
  coeff(2,2,0) = 0.d0
  coeff(2,2,1) = 0.d0
  coeff(2,2,2) = 4.d0*A/N

! alpha3
  coeff(3,3,0) = 0.d0
  coeff(3,3,1) = 0.d0
  coeff(3,3,2) = 0.d0

! alpha4
  coeff(4,4,0) = 0.d0
  coeff(4,4,1) = B/N
  coeff(4,4,2) = 0.d0

! beta1
  coeff(1,2,0) = 0.d0
  coeff(1,2,1) = 0.d0
  coeff(1,2,2) = 4.d0*A/(dsqrt(6.d0)*N)

  coeff(2,1,0) = coeff(1,2,0)
  coeff(2,1,1) = coeff(1,2,1)
  coeff(2,1,2) = coeff(1,2,2)

! beta2
  coeff(3,4,0) = 0.d0
  coeff(3,4,1) = 0.d0
  coeff(3,4,2) = 0.d0

  coeff(4,3,0) = -coeff(3,4,0)
  coeff(4,3,1) = -coeff(3,4,1)
  coeff(4,3,2) = -coeff(3,4,2)

  return
end subroutine RayCaB_Ram_expansion_coeff

max_stars_repo_name: rutgervandeelen/das

max_issues_repo_name: rutgervandeelen/das

max_forks_repo_name: rutgervandeelen/das


text:
MODULE GLOBAL
use assert, only: wp

IMPLICIT NONE

CHARACTER(*), PARAMETER :: VERSION = '1.1.7'

INTEGER :: STACK_SIZE                                         ! stack size
INTEGER, PARAMETER :: REG_SIZE = 10                                           ! number of storage registers
real(wp), allocatable :: STACK(:)                           ! real stack
real(wp), DIMENSION (0:REG_SIZE-1) :: REG                             ! real storage registers
COMPLEX(wp), allocatable :: CSTACK(:)                           ! complex stack
COMPLEX(wp), DIMENSION (0:REG_SIZE-1) :: CREG                             ! complex storage registers
INTEGER, allocatable :: RNSTACK(:), RDSTACK(:)                      ! rational stack
INTEGER, DIMENSION (0:REG_SIZE-1) :: RNREG, RDREG                             ! rational storage registers
real(wp) :: LASTX                                                     ! real LAST X register
COMPLEX(wp) :: CLASTX                                                     ! complex LAST X register
INTEGER :: RNLASTX, RDLASTX                                                   ! rational LAST X register
real(wp) :: NN, SUMX, SUMX2, SUMY, SUMY2, SUMXY                       ! real summation registers
COMPLEX(wp) :: CNN, CSUMX, CSUMX2, CSUMY, CSUMY2, CSUMXY                  ! complex summation registers
INTEGER :: RNNN, RNSUMX, RNSUMX2, RNSUMY, RNSUMY2, RNSUMXY                    ! rational summation registers (num.)
INTEGER :: RDNN, RDSUMX, RDSUMX2, RDSUMY, RDSUMY2, RDSUMXY                    ! rational summation registers (den.)
real(wp) :: ANGLE_FACTOR, FRACTOL
INTEGER :: ANGLE_MODE, DISP_MODE, DISP_DIGITS, BASE_MODE, DOMAIN_MODE, &
   FRACTION_MODE

INTEGER, PARAMETER :: INITIAL_ANGLE_MODE = 1                                  ! 1=deg 2=rad 3=grad 4=rev
INTEGER, PARAMETER :: INITIAL_DISP_MODE = 4                                   ! 1=fix, 2=sci, 3=eng, 4=all
INTEGER, PARAMETER :: INITIAL_DISP_DIGITS = 4
INTEGER, PARAMETER :: INITIAL_DOMAIN_MODE = 1                                 ! 1=real, 2=complex, 3=rational
INTEGER, PARAMETER :: INITIAL_BASE_MODE = 10                                  ! 2=bin, 8=oct, 10=dec, 16=hex
INTEGER, PARAMETER :: INITIAL_FRACTION_MODE = 1                               ! 1=improper, 2=mixed
real(wp), PARAMETER :: INITIAL_FRACTOL = 1.0D-4                       ! tolerance for decimal to fraction conversion

contains

subroutine init_stack()

integer :: i
character(16) :: argv

stack_size = 4  ! default
call get_command_argument(1, argv, status=i)
if (i==0) read(argv,'(I2)') stack_size

if (stack_size < 2) stack_size=2  ! doesn't make sense to have only 1 level

allocate(stack(stack_size), cstack(stack_size), rdstack(stack_size), rnstack(stack_size))

end subroutine init_stack

END MODULE GLOBAL

max_stars_repo_name: majacQ/rpn-calc-fortran

max_issues_repo_name: majacQ/rpn-calc-fortran

max_forks_repo_name: majacQ/rpn-calc-fortran