!-------------------------------------------------------------------------------
!> module ATMOSPHERE / Physics Aerosol Microphysics
!!
!! @par Description
!!          Container for mod_atmos_phy_ae
!!
!! @author Team SCALE
!!
!! @par History
!! @li      2014-05-04 (H.Yashiro)    [new]
!<
!-------------------------------------------------------------------------------
#include "inc_openmp.h"
module mod_atmos_phy_ae_vars
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision
  use scale_stdio
  use scale_prof
  use scale_grid_index
  use scale_tracer
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public procedure
  !
  public :: ATMOS_PHY_AE_vars_setup
  public :: ATMOS_PHY_AE_vars_fillhalo
  public :: ATMOS_PHY_AE_vars_restart_read
  public :: ATMOS_PHY_AE_vars_restart_write

  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  logical,               public :: ATMOS_PHY_AE_RESTART_OUTPUT       = .false.                !< output restart file?

  character(len=H_LONG), public :: ATMOS_PHY_AE_RESTART_IN_BASENAME  = ''                     !< basename of the restart file
  character(len=H_LONG), public :: ATMOS_PHY_AE_RESTART_OUT_BASENAME = ''                     !< basename of the output file
  character(len=H_MID),  public :: ATMOS_PHY_AE_RESTART_OUT_TITLE    = 'ATMOS_PHY_AE restart' !< title    of the output file
  character(len=H_MID),  public :: ATMOS_PHY_AE_RESTART_OUT_DTYPE    = 'DEFAULT'              !< REAL4 or REAL8

  real(RP), public, allocatable :: ATMOS_PHY_AE_RHOQ_t(:,:,:,:)                               ! tendency rho*QTRC [kg/kg/s]

  real(RP), public, allocatable :: ATMOS_PHY_AE_CCN(:,:,:)                                    ! cloud condensation nuclei [/m3]
  real(RP), public, allocatable :: ATMOS_PHY_AE_CCN_t(:,:,:)                                  ! tendency CCN [/m3/s]
  real(RP), public, allocatable :: ATMOS_PHY_AE_EMIT(:,:,:,:)                                 ! emission of aerosol and gas
  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
  integer,                private, parameter :: VMAX = 1       !< number of the variables
  integer,                private, parameter :: I_CCN = 1

  character(len=H_SHORT), private            :: VAR_NAME(VMAX) !< name  of the variables
  character(len=H_MID),   private            :: VAR_DESC(VMAX) !< desc. of the variables
  character(len=H_SHORT), private            :: VAR_UNIT(VMAX) !< unit  of the variables

  data VAR_NAME / 'CCN' /
  data VAR_DESC / 'cloud condensation nuclei' /
  data VAR_UNIT / '/m3' /

  !-----------------------------------------------------------------------------
contains
  !-----------------------------------------------------------------------------
  !> Setup
  subroutine ATMOS_PHY_AE_vars_setup
    use scale_process, only: &
       PRC_MPIstop
    use scale_const, only: &
       UNDEF => CONST_UNDEF
    implicit none

    NAMELIST / PARAM_ATMOS_PHY_AE_VARS / &
       ATMOS_PHY_AE_RESTART_IN_BASENAME,  &
       ATMOS_PHY_AE_RESTART_OUTPUT,       &
       ATMOS_PHY_AE_RESTART_OUT_BASENAME, &
       ATMOS_PHY_AE_RESTART_OUT_TITLE,    &
       ATMOS_PHY_AE_RESTART_OUT_DTYPE

    integer :: ierr
    integer :: iv
    !---------------------------------------------------------------------------

    if( IO_L ) write(IO_FID_LOG,*)
    if( IO_L ) write(IO_FID_LOG,*) '++++++ Module[VARS] / Categ[ATMOS_PHY_AE] / Origin[SCALE-LES]'

    allocate( ATMOS_PHY_AE_RHOQ_t(KA,IA,JA,QA) )
    ATMOS_PHY_AE_RHOQ_t(:,:,:, :) = UNDEF

    allocate( ATMOS_PHY_AE_CCN(KA,IA,JA) )
    ATMOS_PHY_AE_CCN(:,:,:) = UNDEF

    allocate( ATMOS_PHY_AE_CCN_t(KA,IA,JA) )
    ATMOS_PHY_AE_CCN_t(:,:,:) = UNDEF

    allocate( ATMOS_PHY_AE_EMIT(KA,IA,JA,QA_AE) )
    ATMOS_PHY_AE_EMIT(:,:,:,:) = 0.0_RP

    !--- read namelist
    rewind(IO_FID_CONF)
    read(IO_FID_CONF,nml=PARAM_ATMOS_PHY_AE_VARS,iostat=ierr)
    if( ierr < 0 ) then !--- missing
       if( IO_L ) write(IO_FID_LOG,*) '*** Not found namelist. Default used.'
    elseif( ierr > 0 ) then !--- fatal error
       write(*,*) 'xxx Not appropriate names in namelist PARAM_ATMOS_PHY_AE_VARS. Check!'
       call PRC_MPIstop
    endif
    if( IO_LNML ) write(IO_FID_LOG,nml=PARAM_ATMOS_PHY_AE_VARS)

    if( IO_L ) write(IO_FID_LOG,*)
    if( IO_L ) write(IO_FID_LOG,*) '*** [ATMOS_PHY_AE] prognostic/diagnostic variables'
    if( IO_L ) write(IO_FID_LOG,'(1x,A,A15,A,A32,3(A))') &
               '***       |','VARNAME        ','|', 'DESCRIPTION                     ','[', 'UNIT            ',']'
    do iv = 1, VMAX
       if( IO_L ) write(IO_FID_LOG,'(1x,A,i3,A,A15,A,A32,3(A))') &
                  '*** NO.',iv,'|',VAR_NAME(iv),'|',VAR_DESC(iv),'[',VAR_UNIT(iv),']'
    enddo

    if( IO_L ) write(IO_FID_LOG,*)
    if ( ATMOS_PHY_AE_RESTART_IN_BASENAME /= '' ) then
       if( IO_L ) write(IO_FID_LOG,*) '*** Restart input?  : ', trim(ATMOS_PHY_AE_RESTART_IN_BASENAME)
    else
       if( IO_L ) write(IO_FID_LOG,*) '*** Restart input?  : NO'
    endif
    if (       ATMOS_PHY_AE_RESTART_OUTPUT             &
         .AND. ATMOS_PHY_AE_RESTART_OUT_BASENAME /= '' ) then
       if( IO_L ) write(IO_FID_LOG,*) '*** Restart output? : ', trim(ATMOS_PHY_AE_RESTART_OUT_BASENAME)
    else
       if( IO_L ) write(IO_FID_LOG,*) '*** Restart output? : NO'
       ATMOS_PHY_AE_RESTART_OUTPUT = .false.
    endif

    return
  end subroutine ATMOS_PHY_AE_vars_setup

  !-----------------------------------------------------------------------------
  !> HALO Communication
  subroutine ATMOS_PHY_AE_vars_fillhalo
    use scale_comm, only: &
       COMM_vars8, &
       COMM_wait
    implicit none

    integer :: i, j
    !---------------------------------------------------------------------------

    do j  = JS, JE
    do i  = IS, IE
       ATMOS_PHY_AE_CCN(   1:KS-1,i,j) = ATMOS_PHY_AE_CCN(KS,i,j)
       ATMOS_PHY_AE_CCN(KE+1:KA,  i,j) = ATMOS_PHY_AE_CCN(KE,i,j)
    enddo
    enddo

    call COMM_vars8( ATMOS_PHY_AE_CCN(:,:,:), 1 )
    call COMM_wait ( ATMOS_PHY_AE_CCN(:,:,:), 1 )

    return
  end subroutine ATMOS_PHY_AE_vars_fillhalo

  !-----------------------------------------------------------------------------
  !> Read restart
  subroutine ATMOS_PHY_AE_vars_restart_read
    use scale_fileio, only: &
       FILEIO_read
    use scale_les_statistics, only: &
       STAT_total
    implicit none

    real(RP) :: total
    !---------------------------------------------------------------------------

    if( IO_L ) write(IO_FID_LOG,*)
    if( IO_L ) write(IO_FID_LOG,*) '*** Input restart file (ATMOS_PHY_AE) ***'

    if ( ATMOS_PHY_AE_RESTART_IN_BASENAME /= '' ) then
       if( IO_L ) write(IO_FID_LOG,*) '*** basename: ', trim(ATMOS_PHY_AE_RESTART_IN_BASENAME)

       call FILEIO_read( ATMOS_PHY_AE_CCN(:,:,:),                                     & ! [OUT]
                         ATMOS_PHY_AE_RESTART_IN_BASENAME, VAR_NAME(1), 'ZXY', step=1 ) ! [IN]

       call ATMOS_PHY_AE_vars_fillhalo

       call STAT_total( total, ATMOS_PHY_AE_CCN(:,:,:), VAR_NAME(1) )
    else
       if( IO_L ) write(IO_FID_LOG,*) '*** restart file for ATMOS_PHY_AE is not specified.'
    endif

    return
  end subroutine ATMOS_PHY_AE_vars_restart_read

  !-----------------------------------------------------------------------------
  !> Write restart
  subroutine ATMOS_PHY_AE_vars_restart_write
    use scale_time, only: &
       TIME_gettimelabel
    use scale_fileio, only: &
       FILEIO_write
    implicit none

    character(len=15)     :: timelabel
    character(len=H_LONG) :: basename
    !---------------------------------------------------------------------------

    if ( ATMOS_PHY_AE_RESTART_OUT_BASENAME /= '' ) then

       call TIME_gettimelabel( timelabel )
       write(basename,'(A,A,A)') trim(ATMOS_PHY_AE_RESTART_OUT_BASENAME), '_', trim(timelabel)

       if( IO_L ) write(IO_FID_LOG,*)
       if( IO_L ) write(IO_FID_LOG,*) '*** Output restart file (ATMOS_PHY_AE) ***'
       if( IO_L ) write(IO_FID_LOG,*) '*** basename: ', trim(basename)

       call FILEIO_write( ATMOS_PHY_AE_CCN(:,:,:), basename,            ATMOS_PHY_AE_RESTART_OUT_TITLE, & ! [IN]
                          VAR_NAME(1), VAR_DESC(1), VAR_UNIT(1), 'ZXY', ATMOS_PHY_AE_RESTART_OUT_DTYPE  ) ! [IN]

    endif

    return
  end subroutine ATMOS_PHY_AE_vars_restart_write

end module mod_atmos_phy_ae_vars


text: #redirect wiki:Sacramento:Sacramento Galleries

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
module bc_state_outlet_nrbc_lindblad
#include <messenger.h>
    use mod_kinds,              only: rk,ik
    use mod_constants,          only: ZERO, ONE, TWO, HALF, ME, CYLINDRICAL,    &
                                      XI_MIN, XI_MAX, ETA_MIN, ETA_MAX, ZETA_MIN, ZETA_MAX, PI
    use mod_fluid,              only: gam
    use mod_interpolation,      only: interpolate_linear, interpolate_linear_ad
    use mod_gridspace,          only: linspace
    use mod_dft,                only: idft_eval
    use mod_chimera,            only: find_gq_donor, find_gq_donor_parallel

    use type_point,             only: point_t
    use type_mesh,              only: mesh_t
    use type_bc_state,          only: bc_state_t
    use bc_nonlocal_nrbc_lindblad_base,       only: nonlocal_nrbc_lindblad_base_t
    use type_bc_patch,          only: bc_patch_t
    use type_chidg_worker,      only: chidg_worker_t
    use type_properties,        only: properties_t
    use type_face_info,         only: face_info_t, face_info_constructor
    use type_element_info,      only: element_info_t
    use mod_chidg_mpi,          only: IRANK
    use mod_interpolate,        only: interpolate_face_autodiff
    use mpi_f08,                only: MPI_REAL8, MPI_AllReduce, mpi_comm, MPI_INTEGER, MPI_BCast, MPI_MIN, MPI_MAX
    use ieee_arithmetic,        only: ieee_is_nan
    use DNAD_D
    implicit none

    !>  Name: Outlet - 3D Giles
    !!
    !!  Options:
    !!      : Average Pressure
    !!
    !!  Behavior:
    !!
    !!  References:
    !!
    !!  @author Nathan A. Wukie
    !!  @date   2/8/2018
    !!
    !---------------------------------------------------------------------------------
    type, public, extends(nonlocal_nrbc_lindblad_base_t) :: outlet_nrbc_lindblad_t

    contains

        procedure   :: init                 ! Set-up bc state with options/name etc.
        procedure   :: compute_bc_state     ! boundary condition function implementation
        procedure   :: get_q_exterior
        procedure   :: compute_absorbing_outlet
        procedure   :: apply_nonreflecting_condition
        procedure   :: compute_boundary_global_average

    end type outlet_nrbc_lindblad_t
    !*********************************************************************************

contains

    !>
    !!
    !!  @author Nathan A. average_pressure
    !!  @date   2/8/2017
    !!
    !--------------------------------------------------------------------------------
    subroutine init(self)
        class(outlet_nrbc_lindblad_t),   intent(inout) :: self

        ! Set name, family
        call self%set_name('Outlet - Nonlocal NRBC Lindblad')
        call self%set_family('Outlet')

        ! Add functions
        call self%bcproperties%add('Average Pressure',    'Required')
        call self%bcproperties%add('Pitch A',             'Required')
        call self%bcproperties%add('Pitch B',             'Required')
        call self%bcproperties%add('Spatial Periodicity', 'Required')

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

    !>
    !!
    !!  @author Nathan A. Wukie
    !!  @date   2/8/2018
    !!
    !!  @param[in]      worker  Interface for geometry, cache, integration, etc.
    !!  @param[inout]   prop    properties_t object containing equations and material_t objects
    !!
    !-------------------------------------------------------------------------------------------
    subroutine compute_bc_state(self,worker,prop,bc_comm)
        class(outlet_nrbc_lindblad_t),  intent(inout)   :: self
        type(chidg_worker_t),                       intent(inout)   :: worker
        class(properties_t),                        intent(inout)   :: prop
        type(mpi_comm),                             intent(in)      :: bc_comm

        ! Storage at quadrature nodes
        type(AD_D), allocatable, dimension(:)   ::                                                      &
            density_bc, mom1_bc, mom2_bc, mom3_bc, energy_bc, pressure_bc, vel1_bc, vel2_bc, vel3_bc,   &
            density_m, mom1_m, mom2_m, mom3_m, energy_m, pressure_m, vel1_m, vel2_m, vel3_m,   &
            density_reflection, vel1_reflection, vel2_reflection, vel3_reflection, pressure_reflection, &
            !c1_1d, c2_1d, c3_1d, c4_1d, c5_1d,  &
            ddensity, dvel1, dvel2, dvel3, dpressure, &
            c1_3d, c2_3d, c3_3d, c4_3d, &
            density_bar, vel1_bar, vel2_bar, vel3_bar, pressure_bar, c_bar, &
            ddensity_1d, dvel1_1d, dvel2_1d, dvel3_1d, dpressure_1d, &
            ddensity_3d, dvel1_3d, dvel2_3d, dvel3_3d, dpressure_3d, &
            grad1_density_m, grad1_mom1_m, grad1_mom2_m, grad1_mom3_m, grad1_energy_m,                  &
            grad2_density_m, grad2_mom1_m, grad2_mom2_m, grad2_mom3_m, grad2_energy_m,                  &
            grad3_density_m, grad3_mom1_m, grad3_mom2_m, grad3_mom3_m, grad3_energy_m

        type(AD_D), allocatable, dimension(:,:) ::                                  &
            density_check_real_gq, vel1_check_real_gq, vel2_check_real_gq, vel3_check_real_gq, pressure_check_real_gq, &
            density_check_imag_gq, vel1_check_imag_gq, vel2_check_imag_gq, vel3_check_imag_gq, pressure_check_imag_gq

        type(AD_D), allocatable, dimension(:,:,:) ::                                                                                &
            density_check_real_m, vel1_check_real_m, vel2_check_real_m, vel3_check_real_m, pressure_check_real_m, c_check_real_m,   &
            density_check_imag_m, vel1_check_imag_m, vel2_check_imag_m, vel3_check_imag_m, pressure_check_imag_m, c_check_imag_m,   &
            density_hat_real_m, vel1_hat_real_m, vel2_hat_real_m, vel3_hat_real_m, pressure_hat_real_m, c_hat_real_m,               &
            density_hat_imag_m, vel1_hat_imag_m, vel2_hat_imag_m, vel3_hat_imag_m, pressure_hat_imag_m, c_hat_imag_m,               &
            density_check_real_p, vel1_check_real_p, vel2_check_real_p, vel3_check_real_p, pressure_check_real_p, c_check_real_p,   &
            density_check_imag_p, vel1_check_imag_p, vel2_check_imag_p, vel3_check_imag_p, pressure_check_imag_p, c_check_imag_p,   &
            density_hat_real_p, vel1_hat_real_p, vel2_hat_real_p, vel3_hat_real_p, pressure_hat_real_p, c_hat_real_p,               &
            density_hat_imag_p, vel1_hat_imag_p, vel2_hat_imag_p, vel3_hat_imag_p, pressure_hat_imag_p, c_hat_imag_p,               &
            density_hat_real_abs, vel1_hat_real_abs, vel2_hat_real_abs, vel3_hat_real_abs, pressure_hat_real_abs,                   &
            density_hat_imag_abs, vel1_hat_imag_abs, vel2_hat_imag_abs, vel3_hat_imag_abs, pressure_hat_imag_abs,                   &
            density_hat_real_gq, vel1_hat_real_gq, vel2_hat_real_gq, vel3_hat_real_gq, pressure_hat_real_gq,    &
            density_hat_imag_gq, vel1_hat_imag_gq, vel2_hat_imag_gq, vel3_hat_imag_gq, pressure_hat_imag_gq,    &
            density_grid_m, vel1_grid_m, vel2_grid_m, vel3_grid_m, pressure_grid_m, c_grid_m,                   &
            density_grid_p, vel1_grid_p, vel2_grid_p, vel3_grid_p, pressure_grid_p, c_grid_p

        type(AD_D)  :: density_avg, vel1_avg, vel2_avg, vel3_avg, pressure_avg, c_avg, c5_1d

        type(AD_D),     allocatable, dimension(:)   :: p_user, r
        integer(ik) :: igq

        ! Get back pressure from function.
        p_user = self%bcproperties%compute('Average Pressure',worker%time(),worker%coords())

        ! Interpolate interior solution to face quadrature nodes
        grad1_density_m = worker%get_field('Density'   , 'grad1', 'face interior')
        grad2_density_m = worker%get_field('Density'   , 'grad2', 'face interior')
        grad3_density_m = worker%get_field('Density'   , 'grad3', 'face interior')

        grad1_mom1_m    = worker%get_field('Momentum-1', 'grad1', 'face interior')
        grad2_mom1_m    = worker%get_field('Momentum-1', 'grad2', 'face interior')
        grad3_mom1_m    = worker%get_field('Momentum-1', 'grad3', 'face interior')

        grad1_mom2_m    = worker%get_field('Momentum-2', 'grad1', 'face interior')
        grad2_mom2_m    = worker%get_field('Momentum-2', 'grad2', 'face interior')
        grad3_mom2_m    = worker%get_field('Momentum-2', 'grad3', 'face interior')

        grad1_mom3_m    = worker%get_field('Momentum-3', 'grad1', 'face interior')
        grad2_mom3_m    = worker%get_field('Momentum-3', 'grad2', 'face interior')
        grad3_mom3_m    = worker%get_field('Momentum-3', 'grad3', 'face interior')

        grad1_energy_m  = worker%get_field('Energy'    , 'grad1', 'face interior')
        grad2_energy_m  = worker%get_field('Energy'    , 'grad2', 'face interior')
        grad3_energy_m  = worker%get_field('Energy'    , 'grad3', 'face interior')

        ! Store boundary gradient state. Grad(Q_bc). Do this here, before we
        ! compute any transformations for cylindrical.
        call worker%store_bc_state('Density'   , grad1_density_m, 'grad1')
        call worker%store_bc_state('Density'   , grad2_density_m, 'grad2')
        call worker%store_bc_state('Density'   , grad3_density_m, 'grad3')

        call worker%store_bc_state('Momentum-1', grad1_mom1_m,    'grad1')
        call worker%store_bc_state('Momentum-1', grad2_mom1_m,    'grad2')
        call worker%store_bc_state('Momentum-1', grad3_mom1_m,    'grad3')

        call worker%store_bc_state('Momentum-2', grad1_mom2_m,    'grad1')
        call worker%store_bc_state('Momentum-2', grad2_mom2_m,    'grad2')
        call worker%store_bc_state('Momentum-2', grad3_mom2_m,    'grad3')

        call worker%store_bc_state('Momentum-3', grad1_mom3_m,    'grad1')
        call worker%store_bc_state('Momentum-3', grad2_mom3_m,    'grad2')
        call worker%store_bc_state('Momentum-3', grad3_mom3_m,    'grad3')

        call worker%store_bc_state('Energy'    , grad1_energy_m,  'grad1')
        call worker%store_bc_state('Energy'    , grad2_energy_m,  'grad2')
        call worker%store_bc_state('Energy'    , grad3_energy_m,  'grad3')

        ! Get primitive variables at (radius,theta,time) grid.
        call self%get_q_side(worker,bc_comm,'A', &
                             density_grid_m,     &
                             vel1_grid_m,        &
                             vel2_grid_m,        &
                             vel3_grid_m,        &
                             pressure_grid_m)
        c_grid_m = sqrt(gam*pressure_grid_m/density_grid_m)

        ! Compute Fourier decomposition of temporal data at points
        ! on the spatial transform grid.
        !   q_check(r,theta,omega) = DFT(q)[time]
        call self%compute_temporal_dft(worker,bc_comm,                                  &
                                       density_grid_m,                                  &
                                       vel1_grid_m,                                     &
                                       vel2_grid_m,                                     &
                                       vel3_grid_m,                                     &
                                       pressure_grid_m,                                 &
                                       c_grid_m,                                        &
                                       density_check_real_m,  density_check_imag_m,     &
                                       vel1_check_real_m,     vel1_check_imag_m,        &
                                       vel2_check_real_m,     vel2_check_imag_m,        &
                                       vel3_check_real_m,     vel3_check_imag_m,        &
                                       pressure_check_real_m, pressure_check_imag_m,    &
                                       c_check_real_m,        c_check_imag_m)

        ! Compute Fourier decomposition in theta at set of radial
        ! stations for each temporal mode:
        !   q_hat(r,m,omega) = DFT(q_check)[theta]
        call self%compute_spatial_dft(worker,bc_comm,'A',                           &
                                      density_check_real_m,  density_check_imag_m,  &
                                      vel1_check_real_m,     vel1_check_imag_m,     &
                                      vel2_check_real_m,     vel2_check_imag_m,     &
                                      vel3_check_real_m,     vel3_check_imag_m,     &
                                      pressure_check_real_m, pressure_check_imag_m, &
                                      c_check_real_m,        c_check_imag_m,        &
                                      density_hat_real_m,    density_hat_imag_m,    &
                                      vel1_hat_real_m,       vel1_hat_imag_m,       &
                                      vel2_hat_real_m,       vel2_hat_imag_m,       &
                                      vel3_hat_real_m,       vel3_hat_imag_m,       &
                                      pressure_hat_real_m,   pressure_hat_imag_m,   &
                                      c_hat_real_m,          c_hat_imag_m)

        ! Get primitive variables at (radius,theta,time) grid.
        call self%get_q_exterior(worker,bc_comm,    &
                                 density_grid_p,    &
                                 vel1_grid_p,       &
                                 vel2_grid_p,       &
                                 vel3_grid_p,       &
                                 pressure_grid_p)
        c_grid_p = sqrt(gam*pressure_grid_p/density_grid_p)

        ! Compute Fourier decomposition of temporal data at points
        ! on the spatial transform grid.
        !   q_check(r,theta,omega) = DFT(q)[time]
        call self%compute_temporal_dft(worker,bc_comm,                                  &
                                       density_grid_p,                                  &
                                       vel1_grid_p,                                     &
                                       vel2_grid_p,                                     &
                                       vel3_grid_p,                                     &
                                       pressure_grid_p,                                 &
                                       c_grid_p,                                        &
                                       density_check_real_p,  density_check_imag_p,     &
                                       vel1_check_real_p,     vel1_check_imag_p,        &
                                       vel2_check_real_p,     vel2_check_imag_p,        &
                                       vel3_check_real_p,     vel3_check_imag_p,        &
                                       pressure_check_real_p, pressure_check_imag_p,    &
                                       c_check_real_p,        c_check_imag_p)

        ! Compute Fourier decomposition in theta at set of radial
        ! stations for each temporal mode:
        !   q_hat(r,m,omega) = DFT(q_check)[theta]
        call self%compute_spatial_dft(worker,bc_comm,'B',                           &
                                      density_check_real_p,  density_check_imag_p,  &
                                      vel1_check_real_p,     vel1_check_imag_p,     &
                                      vel2_check_real_p,     vel2_check_imag_p,     &
                                      vel3_check_real_p,     vel3_check_imag_p,     &
                                      pressure_check_real_p, pressure_check_imag_p, &
                                      c_check_real_p,        c_check_imag_p,        &
                                      density_hat_real_p,    density_hat_imag_p,    &
                                      vel1_hat_real_p,       vel1_hat_imag_p,       &
                                      vel2_hat_real_p,       vel2_hat_imag_p,       &
                                      vel3_hat_real_p,       vel3_hat_imag_p,       &
                                      pressure_hat_real_p,   pressure_hat_imag_p,   &
                                      c_hat_real_p,          c_hat_imag_p)

        ! Compute absorbing state using 2D eigenmode decomposition.
        !
        !   q_abs = f(q_hat_m,q_hat_p)
        !
        call self%compute_absorbing_outlet(worker,bc_comm,                              &
                                           density_hat_real_m,    density_hat_imag_m,   &
                                           vel1_hat_real_m,       vel1_hat_imag_m,      &
                                           vel2_hat_real_m,       vel2_hat_imag_m,      &
                                           vel3_hat_real_m,       vel3_hat_imag_m,      &
                                           pressure_hat_real_m,   pressure_hat_imag_m,  &
                                           c_hat_real_m,          c_hat_imag_m,         &
                                           density_hat_real_p,    density_hat_imag_p,   &
                                           vel1_hat_real_p,       vel1_hat_imag_p,      &
                                           vel2_hat_real_p,       vel2_hat_imag_p,      &
                                           vel3_hat_real_p,       vel3_hat_imag_p,      &
                                           pressure_hat_real_p,   pressure_hat_imag_p,  &
                                           c_hat_real_p,          c_hat_imag_p,         &
                                           density_hat_real_abs,  density_hat_imag_abs, &
                                           vel1_hat_real_abs,     vel1_hat_imag_abs,    &
                                           vel2_hat_real_abs,     vel2_hat_imag_abs,    &
                                           vel3_hat_real_abs,     vel3_hat_imag_abs,    &
                                           pressure_hat_real_abs, pressure_hat_imag_abs)

        ! q_abs(r_gq) = I(q_abs(r_aux))
        call self%interpolate_raux_to_rgq(worker,bc_comm,                               &
                                          density_hat_real_abs,  density_hat_imag_abs,  &
                                          vel1_hat_real_abs,     vel1_hat_imag_abs,     &
                                          vel2_hat_real_abs,     vel2_hat_imag_abs,     &
                                          vel3_hat_real_abs,     vel3_hat_imag_abs,     &
                                          pressure_hat_real_abs, pressure_hat_imag_abs, &
                                          density_hat_real_gq,   density_hat_imag_gq,   &
                                          vel1_hat_real_gq,      vel1_hat_imag_gq,      &
                                          vel2_hat_real_gq,      vel2_hat_imag_gq,      &
                                          vel3_hat_real_gq,      vel3_hat_imag_gq,      &
                                          pressure_hat_real_gq,  pressure_hat_imag_gq)

        ! Reconstruct primitive variables at quadrature nodes from absorbing Fourier modes
        ! via inverse transform.
        !   q_check(rgq,theta,omega) = IDFT(q_hat)[m]
        call self%compute_spatial_idft_gq(worker,bc_comm,'A',                               &
                                          density_hat_real_gq,      density_hat_imag_gq,    &
                                          vel1_hat_real_gq,         vel1_hat_imag_gq,       &
                                          vel2_hat_real_gq,         vel2_hat_imag_gq,       &
                                          vel3_hat_real_gq,         vel3_hat_imag_gq,       &
                                          pressure_hat_real_gq,     pressure_hat_imag_gq,   &
                                          density_check_real_gq,    density_check_imag_gq,  &
                                          vel1_check_real_gq,       vel1_check_imag_gq,     &
                                          vel2_check_real_gq,       vel2_check_imag_gq,     &
                                          vel3_check_real_gq,       vel3_check_imag_gq,     &
                                          pressure_check_real_gq,   pressure_check_imag_gq)

        ! q(rgq,theta,t) = IDFT(q_check)[omega]
        call self%compute_temporal_idft_gq(worker,bc_comm,                                      &
                                           density_check_real_gq,    density_check_imag_gq,     &
                                           vel1_check_real_gq,       vel1_check_imag_gq,        &
                                           vel2_check_real_gq,       vel2_check_imag_gq,        &
                                           vel3_check_real_gq,       vel3_check_imag_gq,        &
                                           pressure_check_real_gq,   pressure_check_imag_gq,    &
                                           density_bc, vel1_bc, vel2_bc, vel3_bc, pressure_bc)

        !
        ! Form conserved variables
        !
        density_bc = density_bc
        mom1_bc    = density_bc*vel1_bc
        mom2_bc    = density_bc*vel2_bc
        mom3_bc    = density_bc*vel3_bc
        energy_bc  = pressure_bc/(gam - ONE)  + HALF*(mom1_bc*mom1_bc + mom2_bc*mom2_bc + mom3_bc*mom3_bc)/density_bc

        !
        ! Account for cylindrical. Convert tangential momentum back to angular momentum.
        !
        if (worker%coordinate_system() == 'Cylindrical') then
            r = worker%coordinate('1','boundary')
            mom2_bc = mom2_bc * r
        end if

        !
        ! Store boundary condition state. Q_bc
        !
        call worker%store_bc_state('Density'   , density_bc, 'value')
        call worker%store_bc_state('Momentum-1', mom1_bc,    'value')
        call worker%store_bc_state('Momentum-2', mom2_bc,    'value')
        call worker%store_bc_state('Momentum-3', mom3_bc,    'value')
        call worker%store_bc_state('Energy'    , energy_bc,  'value')

    end subroutine compute_bc_state
    !*********************************************************************************

    !>
    !!
    !!  @author Nathan A. Wukie
    !!  @date   6/18/2018
    !!
    !------------------------------------------------------------------------------------
    subroutine get_q_exterior(self,worker,bc_comm, density, vel1, vel2, vel3, pressure)
        class(outlet_nrbc_lindblad_t),  intent(inout)   :: self
        type(chidg_worker_t),                       intent(inout)   :: worker
        type(mpi_comm),                             intent(in)      :: bc_comm
        type(AD_D),     allocatable,                intent(inout)   :: density(:,:,:)
        type(AD_D),     allocatable,                intent(inout)   :: vel1(:,:,:)
        type(AD_D),     allocatable,                intent(inout)   :: vel2(:,:,:)
        type(AD_D),     allocatable,                intent(inout)   :: vel3(:,:,:)
        type(AD_D),     allocatable,                intent(inout)   :: pressure(:,:,:)

        type(AD_D), allocatable, dimension(:,:,:) ::  &
            mom1, mom2, mom3, energy

        real(rk),   allocatable, dimension(:)   :: p_user

        integer(ik) :: iradius, itime, nradius, ntheta, ntime, ierr

        ! Define Fourier space discretization to determine
        ! number of theta-samples are being taken
        nradius = size(self%r)
        ntheta  = size(self%theta_a,2)
        ntime   = worker%time_manager%ntime

        ! Allocate storage for discrete time instances
        allocate(density( nradius,ntheta,ntime),  &
                 mom1(    nradius,ntheta,ntime),  &
                 mom2(    nradius,ntheta,ntime),  &
                 mom3(    nradius,ntheta,ntime),  &
                 vel1(    nradius,ntheta,ntime),  &
                 vel2(    nradius,ntheta,ntime),  &
                 vel3(    nradius,ntheta,ntime),  &
                 energy(  nradius,ntheta,ntime),  &
                 pressure(nradius,ntheta,ntime), stat=ierr)
        if (ierr /= 0) call AllocationError

        ! Perform Fourier decomposition at each radial station.
        do iradius = 1,nradius
            do itime = 1,ntime
                ! Interpolate solution to physical_nodes at current radial station: [ntheta]
                density(iradius,:,itime) = worker%interpolate_field_general('Density',    donors=self%donor_a(iradius,:), donor_nodes=self%donor_node_a(iradius,:,:), itime=itime)
                mom1(iradius,:,itime)    = worker%interpolate_field_general('Momentum-1', donors=self%donor_a(iradius,:), donor_nodes=self%donor_node_a(iradius,:,:), itime=itime)
                mom2(iradius,:,itime)    = worker%interpolate_field_general('Momentum-2', donors=self%donor_a(iradius,:), donor_nodes=self%donor_node_a(iradius,:,:), itime=itime)
                mom3(iradius,:,itime)    = worker%interpolate_field_general('Momentum-3', donors=self%donor_a(iradius,:), donor_nodes=self%donor_node_a(iradius,:,:), itime=itime)
                energy(iradius,:,itime)  = worker%interpolate_field_general('Energy',     donors=self%donor_a(iradius,:), donor_nodes=self%donor_node_a(iradius,:,:), itime=itime)

                if (worker%coordinate_system() == 'Cylindrical') then
                    mom2(iradius,:,itime) = mom2(iradius,:,itime)/self%r(iradius)  ! convert to tangential momentum
                end if
            end do
        end do

        ! Compute velocities and pressure at each time
        density  = density
        vel1     = mom1/density
        vel2     = mom2/density
        vel3     = mom3/density

        ! Pressure from user
        pressure = (gam-ONE)*(energy - HALF*(mom1*mom1 + mom2*mom2 + mom3*mom3)/density)
        p_user = self%bcproperties%compute('Average Pressure',worker%time(),[point_t(ZERO,ZERO,ZERO)])
        pressure = p_user(1)

    end subroutine get_q_exterior
    !************************************************************************************

    !> DANIEL'S FORMULATION
    !!
    !!
    !!
    !!
    !--------------------------------------------------------------------------------
    subroutine compute_absorbing_outlet(self,worker,bc_comm,                    &
                                     density_real_m,    density_imag_m,      &
                                     vel1_real_m,       vel1_imag_m,         &
                                     vel2_real_m,       vel2_imag_m,         &
                                     vel3_real_m,       vel3_imag_m,         &
                                     pressure_real_m,   pressure_imag_m,     &
                                     c_real_m,          c_imag_m,            &
                                     density_real_p,    density_imag_p,      &
                                     vel1_real_p,       vel1_imag_p,         &
                                     vel2_real_p,       vel2_imag_p,         &
                                     vel3_real_p,       vel3_imag_p,         &
                                     pressure_real_p,   pressure_imag_p,     &
                                     c_real_p,          c_imag_p,            &
                                     density_real_abs,  density_imag_abs,    &
                                     vel1_real_abs,     vel1_imag_abs,       &
                                     vel2_real_abs,     vel2_imag_abs,       &
                                     vel3_real_abs,     vel3_imag_abs,       &
                                     pressure_real_abs, pressure_imag_abs)
        class(outlet_nrbc_lindblad_t),     intent(inout)   :: self
        type(chidg_worker_t),       intent(inout)   :: worker
        type(mpi_comm),             intent(in)      :: bc_comm
        type(AD_D),                 intent(inout)   :: density_real_m(:,:,:)
        type(AD_D),                 intent(inout)   :: density_imag_m(:,:,:)
        type(AD_D),                 intent(inout)   :: vel1_real_m(:,:,:)
        type(AD_D),                 intent(inout)   :: vel1_imag_m(:,:,:)
        type(AD_D),                 intent(inout)   :: vel2_real_m(:,:,:)
        type(AD_D),                 intent(inout)   :: vel2_imag_m(:,:,:)
        type(AD_D),                 intent(inout)   :: vel3_real_m(:,:,:)
        type(AD_D),                 intent(inout)   :: vel3_imag_m(:,:,:)
        type(AD_D),                 intent(inout)   :: pressure_real_m(:,:,:)
        type(AD_D),                 intent(inout)   :: pressure_imag_m(:,:,:)
        type(AD_D),                 intent(inout)   :: c_real_m(:,:,:)
        type(AD_D),                 intent(inout)   :: c_imag_m(:,:,:)
        type(AD_D),                 intent(inout)   :: density_real_p(:,:,:)
        type(AD_D),                 intent(inout)   :: density_imag_p(:,:,:)
        type(AD_D),                 intent(inout)   :: vel1_real_p(:,:,:)
        type(AD_D),                 intent(inout)   :: vel1_imag_p(:,:,:)
        type(AD_D),                 intent(inout)   :: vel2_real_p(:,:,:)
        type(AD_D),                 intent(inout)   :: vel2_imag_p(:,:,:)
        type(AD_D),                 intent(inout)   :: vel3_real_p(:,:,:)
        type(AD_D),                 intent(inout)   :: vel3_imag_p(:,:,:)
        type(AD_D),                 intent(inout)   :: pressure_real_p(:,:,:)
        type(AD_D),                 intent(inout)   :: pressure_imag_p(:,:,:)
        type(AD_D),                 intent(inout)   :: c_real_p(:,:,:)
        type(AD_D),                 intent(inout)   :: c_imag_p(:,:,:)
        type(AD_D), allocatable,    intent(inout)   :: density_real_abs(:,:,:)
        type(AD_D), allocatable,    intent(inout)   :: density_imag_abs(:,:,:)
        type(AD_D), allocatable,    intent(inout)   :: vel1_real_abs(:,:,:)
        type(AD_D), allocatable,    intent(inout)   :: vel1_imag_abs(:,:,:)
        type(AD_D), allocatable,    intent(inout)   :: vel2_real_abs(:,:,:)
        type(AD_D), allocatable,    intent(inout)   :: vel2_imag_abs(:,:,:)
        type(AD_D), allocatable,    intent(inout)   :: vel3_real_abs(:,:,:)
        type(AD_D), allocatable,    intent(inout)   :: vel3_imag_abs(:,:,:)
        type(AD_D), allocatable,    intent(inout)   :: pressure_real_abs(:,:,:)
        type(AD_D), allocatable,    intent(inout)   :: pressure_imag_abs(:,:,:)

        type(AD_D), allocatable, dimension(:,:,:) ::                &
            a1_real, a2_real, a3_real, a4_real, a5_real,            &
            a1_imag, a2_imag, a3_imag, a4_imag, a5_imag,            &
            a1_real_m, a2_real_m, a3_real_m, a4_real_m, a5_real_m,  &
            a1_imag_m, a2_imag_m, a3_imag_m, a4_imag_m, a5_imag_m,  &
            a1_real_p, a2_real_p, a3_real_p, a4_real_p, a5_real_p,  &
            a1_imag_p, a2_imag_p, a3_imag_p, a4_imag_p, a5_imag_p,  &
            c1_real_m, c2_real_m, c3_real_m, c4_real_m, c5_real_m,  &
            c1_imag_m, c2_imag_m, c3_imag_m, c4_imag_m, c5_imag_m

        type(AD_D)  :: density_avg, vel1_avg, vel2_avg, vel3_avg, pressure_avg, c_avg, &
                       density_bar, vel1_bar, vel2_bar, vel3_bar, pressure_bar, c_bar, &
                       B3_real, B3_imag, B4_real, B4_imag, beta

        type(AD_D), allocatable, dimension(:)   :: ddensity, dvel1, dvel2, dvel3, dpressure, c1_1d, c2_1d, c3_1d, c4_1d, c5_1d

        type(AD_D), allocatable, dimension(:)   :: p_user
        integer(ik) :: iradius, itheta, itime, ntheta

        ! Retrieve target average pressure
        p_user = self%bcproperties%compute('Average Pressure',worker%time(),worker%coords())

        !
        ! Step 1: handle higher-order modes: m /= 0
        !
        ! Project interior to eigenmodes
        call self%primitive_to_eigenmodes(worker,bc_comm,                   &
                                          density_real_m(:,1,1),            &
                                          vel1_real_m(:,1,1),               &
                                          vel2_real_m(:,1,1),               &
                                          vel3_real_m(:,1,1),               &
                                          pressure_real_m(:,1,1),           &
                                          c_real_m(:,1,1),                  &
                                          density_real_m,  density_imag_m,  &
                                          vel1_real_m,     vel1_imag_m,     &
                                          vel2_real_m,     vel2_imag_m,     &
                                          vel3_real_m,     vel3_imag_m,     &
                                          pressure_real_m, pressure_imag_m, &
                                          a1_real_m,       a1_imag_m,       &
                                          a2_real_m,       a2_imag_m,       &
                                          a3_real_m,       a3_imag_m,       &
                                          a4_real_m,       a4_imag_m,       &
                                          a5_real_m,       a5_imag_m)

        ! Project exterior to eigenmodes
        call self%primitive_to_eigenmodes(worker,bc_comm,                   &
                                          density_real_m(:,1,1),            &
                                          vel1_real_m(:,1,1),               &
                                          vel2_real_m(:,1,1),               &
                                          vel3_real_m(:,1,1),               &
                                          pressure_real_m(:,1,1),           &
                                          c_real_m(:,1,1),                  &
                                          density_real_p,  density_imag_p,  &
                                          vel1_real_p,     vel1_imag_p,     &
                                          vel2_real_p,     vel2_imag_p,     &
                                          vel3_real_p,     vel3_imag_p,     &
                                          pressure_real_p, pressure_imag_p, &
                                          a1_real_p,       a1_imag_p,       &
                                          a2_real_p,       a2_imag_p,       &
                                          a3_real_p,       a3_imag_p,       &
                                          a4_real_p,       a4_imag_p,       &
                                          a5_real_p,       a5_imag_p)

        call self%apply_nonreflecting_condition(worker,bc_comm,             &
                                                density_real_m(:,1,1),      &
                                                vel1_real_m(:,1,1),         &
                                                vel2_real_m(:,1,1),         &
                                                vel3_real_m(:,1,1),         &
                                                pressure_real_m(:,1,1),     &
                                                c_real_m(:,1,1),            &
                                                a1_real_m,    a1_imag_m,    &
                                                a2_real_m,    a2_imag_m,    &
                                                a3_real_m,    a3_imag_m,    &
                                                a4_real_m,    a4_imag_m,    &
                                                a5_real_m,    a5_imag_m,    &
                                                a1_real_p,    a1_imag_p,    &
                                                a2_real_p,    a2_imag_p,    &
                                                a3_real_p,    a3_imag_p,    &
                                                a4_real_p,    a4_imag_p,    &
                                                a5_real_p,    a5_imag_p,    &
                                                a1_real,      a1_imag,      &
                                                a2_real,      a2_imag,      &
                                                a3_real,      a3_imag,      &
                                                a4_real,      a4_imag,      &
                                                a5_real,      a5_imag)

        ! Convert back to primitive variables
        call self%eigenmodes_to_primitive(worker,bc_comm,                       &
                                          density_real_m(:,1,1),                &
                                          vel1_real_m(:,1,1),                   &
                                          vel2_real_m(:,1,1),                   &
                                          vel3_real_m(:,1,1),                   &
                                          pressure_real_m(:,1,1),               &
                                          c_real_m(:,1,1),                      &
                                          a1_real,           a1_imag,           &
                                          a2_real,           a2_imag,           &
                                          a3_real,           a3_imag,           &
                                          a4_real,           a4_imag,           &
                                          a5_real,           a5_imag,           &
                                          density_real_abs,  density_imag_abs,  &
                                          vel1_real_abs,     vel1_imag_abs,     &
                                          vel2_real_abs,     vel2_imag_abs,     &
                                          vel3_real_abs,     vel3_imag_abs,     &
                                          pressure_real_abs, pressure_imag_abs)

        !
        ! Step 2: handle m=0(radius-local space-time average) using 1D characteristics
        !         based on the perturbation of the local m=0 quantity from the boundary
        !         global space-time average.
        !

        ! Handle radius-local mean variation from global boundary average using 1D characteristics
        call self%compute_boundary_global_average(worker,bc_comm,density_real_m(:,1,1),     &
                                                                 vel1_real_m(:,1,1),        &
                                                                 vel2_real_m(:,1,1),        &
                                                                 vel3_real_m(:,1,1),        &
                                                                 pressure_real_m(:,1,1),    &
                                                                 density_avg,vel1_avg,vel2_avg,vel3_avg,pressure_avg)
        c_avg = sqrt(gam*pressure_avg/density_avg)

        ! Compute perturbation of radius local-avg from boundary-global avg
        ddensity  = density_real_m(:,1,1)  - density_avg
        dvel1     = vel1_real_m(:,1,1)     - vel1_avg
        dvel2     = vel2_real_m(:,1,1)     - vel2_avg
        dvel3     = vel3_real_m(:,1,1)     - vel3_avg
        dpressure = pressure_real_m(:,1,1) - pressure_avg

        ! Allocate/compute 1d characteristics
        c1_1d = ZERO*ddensity
        c2_1d = ZERO*ddensity
        c3_1d = ZERO*ddensity
        c4_1d = ZERO*ddensity
        c5_1d = ZERO*ddensity
        do iradius = 1,size(ddensity)
            c1_1d(iradius) = -c_avg*c_avg*ddensity(iradius)    +  dpressure(iradius)
            c2_1d(iradius) = density_avg*c_avg*dvel1(iradius)
            c3_1d(iradius) = density_avg*c_avg*dvel2(iradius)
            c4_1d(iradius) = density_avg*c_avg*dvel3(iradius)  +  dpressure(iradius)
            c5_1d(iradius) = -TWO*(pressure_avg - p_user(1))
        end do

        ! Begin composition of radius-local absorbing avg from boundary-global average
        density_real_abs(:,1,1)  = density_avg
        vel1_real_abs(:,1,1)     = vel1_avg
        vel2_real_abs(:,1,1)     = vel2_avg
        vel3_real_abs(:,1,1)     = vel3_avg
        pressure_real_abs(:,1,1) = pressure_avg

        ! Now add local perturbation from the average
        do iradius = 1,size(c1_1d)
            density_real_abs(iradius,1,1)  = density_real_abs(iradius,1,1)  + (-ONE/(c_avg*c_avg))*c1_1d(iradius)  +  (ONE/(TWO*c_avg*c_avg))*c4_1d(iradius)  +  (ONE/(TWO*c_avg*c_avg))*c5_1d(iradius)
            vel1_real_abs(iradius,1,1)     = vel1_real_abs(iradius,1,1)     + (ONE/(density_avg*c_avg))*c2_1d(iradius)
            vel2_real_abs(iradius,1,1)     = vel2_real_abs(iradius,1,1)     + (ONE/(density_avg*c_avg))*c3_1d(iradius)
            vel3_real_abs(iradius,1,1)     = vel3_real_abs(iradius,1,1)     + (ONE/(TWO*density_avg*c_avg))*c4_1d(iradius)  -  (ONE/(TWO*density_avg*c_avg))*c5_1d(iradius)
            pressure_real_abs(iradius,1,1) = pressure_real_abs(iradius,1,1) + HALF*c4_1d(iradius)  +  HALF*c5_1d(iradius)
        end do

    end subroutine compute_absorbing_outlet
    !********************************************************************************

    !>
    !!
    !!
    !!
    !-----------------------------------------------------------------------------------
    subroutine apply_nonreflecting_condition(self,worker,bc_comm, &
                                             density_bar_r, vel1_bar_r, vel2_bar_r, vel3_bar_r, pressure_bar_r, c_bar_r, &
                                             a1_real_m,    a1_imag_m,    &
                                             a2_real_m,    a2_imag_m,    &
                                             a3_real_m,    a3_imag_m,    &
                                             a4_real_m,    a4_imag_m,    &
                                             a5_real_m,    a5_imag_m,    &
                                             a1_real_p,    a1_imag_p,    &
                                             a2_real_p,    a2_imag_p,    &
                                             a3_real_p,    a3_imag_p,    &
                                             a4_real_p,    a4_imag_p,    &
                                             a5_real_p,    a5_imag_p,    &
                                             a1_real,      a1_imag,      &
                                             a2_real,      a2_imag,      &
                                             a3_real,      a3_imag,      &
                                             a4_real,      a4_imag,      &
                                             a5_real,      a5_imag)
        class(outlet_nrbc_lindblad_t),     intent(inout)   :: self
        type(chidg_worker_t),       intent(inout)   :: worker
        type(mpi_comm),             intent(in)      :: bc_comm
        type(AD_D),                 intent(in)      :: density_bar_r(:)
        type(AD_D),                 intent(in)      :: vel1_bar_r(:)
        type(AD_D),                 intent(in)      :: vel2_bar_r(:)
        type(AD_D),                 intent(in)      :: vel3_bar_r(:)
        type(AD_D),                 intent(in)      :: pressure_bar_r(:)
        type(AD_D),                 intent(in)      :: c_bar_r(:)
        type(AD_D),                 intent(inout)   :: a1_real_m(:,:,:)
        type(AD_D),                 intent(inout)   :: a1_imag_m(:,:,:)
        type(AD_D),                 intent(inout)   :: a2_real_m(:,:,:)
        type(AD_D),                 intent(inout)   :: a2_imag_m(:,:,:)
        type(AD_D),                 intent(inout)   :: a3_real_m(:,:,:)
        type(AD_D),                 intent(inout)   :: a3_imag_m(:,:,:)
        type(AD_D),                 intent(inout)   :: a4_real_m(:,:,:)
        type(AD_D),                 intent(inout)   :: a4_imag_m(:,:,:)
        type(AD_D),                 intent(inout)   :: a5_real_m(:,:,:)
        type(AD_D),                 intent(inout)   :: a5_imag_m(:,:,:)
        type(AD_D),                 intent(inout)   :: a1_real_p(:,:,:)
        type(AD_D),                 intent(inout)   :: a1_imag_p(:,:,:)
        type(AD_D),                 intent(inout)   :: a2_real_p(:,:,:)
        type(AD_D),                 intent(inout)   :: a2_imag_p(:,:,:)
        type(AD_D),                 intent(inout)   :: a3_real_p(:,:,:)
        type(AD_D),                 intent(inout)   :: a3_imag_p(:,:,:)
        type(AD_D),                 intent(inout)   :: a4_real_p(:,:,:)
        type(AD_D),                 intent(inout)   :: a4_imag_p(:,:,:)
        type(AD_D),                 intent(inout)   :: a5_real_p(:,:,:)
        type(AD_D),                 intent(inout)   :: a5_imag_p(:,:,:)
        type(AD_D), allocatable,    intent(inout)   :: a1_real(:,:,:)
        type(AD_D), allocatable,    intent(inout)   :: a1_imag(:,:,:)
        type(AD_D), allocatable,    intent(inout)   :: a2_real(:,:,:)
        type(AD_D), allocatable,    intent(inout)   :: a2_imag(:,:,:)
        type(AD_D), allocatable,    intent(inout)   :: a3_real(:,:,:)
        type(AD_D), allocatable,    intent(inout)   :: a3_imag(:,:,:)
        type(AD_D), allocatable,    intent(inout)   :: a4_real(:,:,:)
        type(AD_D), allocatable,    intent(inout)   :: a4_imag(:,:,:)
        type(AD_D), allocatable,    intent(inout)   :: a5_real(:,:,:)
        type(AD_D), allocatable,    intent(inout)   :: a5_imag(:,:,:)

        integer(ik) :: iradius, itheta, ntheta, itime

        type(AD_D)  :: beta, B3_real, B3_imag, B4_real, B4_imag, &
                       density_bar, vel1_bar, vel2_bar, vel3_bar, pressure_bar, c_bar

        ! Outoing amplitudes from interior
        a1_real = a1_real_m
        a1_imag = a1_imag_m
        a2_real = a2_real_m
        a2_imag = a2_imag_m
        a3_real = a3_real_m
        a3_imag = a3_imag_m
        a4_real = a4_real_m
        a4_imag = a4_imag_m

        ! Incoming amplitudes from exterior
        a5_real = ZERO*a5_real_p
        a5_imag = ZERO*a5_imag_p

!        itime = 1   ! Time-constant
!        do iradius = 1,size(a1_real_m,1)
!            ! Get average parts
!            density_bar  = density_bar_r(iradius)
!            vel1_bar     = vel1_bar_r(iradius)
!            vel2_bar     = vel2_bar_r(iradius)
!            vel3_bar     = vel3_bar_r(iradius)
!            pressure_bar = pressure_bar_r(iradius)
!            c_bar        = sqrt(gam*pressure_bar/density_bar)
!
!            ! starting with 2 here because the first mode is treated with 1D characteristics
!            ntheta = size(a1_real_m,2)
!            do itheta = 2,ntheta
!                ! Account for sign(mode) in the calculation of beta. The second half of the
!                ! modes are negative frequencies.
!                if (itheta <= (ntheta-1)/2 + 1) then
!                    beta = sqrt(c_bar*c_bar  -  (vel3_bar*vel3_bar + vel2_bar*vel2_bar))
!                else if (itheta > (ntheta-1)/2 + 1) then
!                    beta = -sqrt(c_bar*c_bar  -  (vel3_bar*vel3_bar + vel2_bar*vel2_bar))
!                end if
!
!                ! The imaginary part of beta has already been accounted for in
!                ! the expressions for B2 and B3
!                B3_real = -TWO*vel3_bar*vel2_bar/(vel2_bar*vel2_bar + beta*beta)
!                B3_imag = -TWO*beta*vel3_bar/(vel2_bar*vel2_bar + beta*beta)
!
!                B4_real = (beta*beta - vel2_bar*vel2_bar)/(beta*beta + vel2_bar*vel2_bar)
!                B4_imag = -TWO*beta*vel2_bar/(beta*beta + vel2_bar*vel2_bar)
!
!                a5_real(iradius,itheta,itime) = (B3_real*a3_real_m(iradius,itheta,itime) - B3_imag*a3_imag_m(iradius,itheta,itime))  &   ! A3*c3 (real)
!                                              - (B4_real*a4_real_m(iradius,itheta,itime) - B4_imag*a4_imag_m(iradius,itheta,itime))      ! A4*c4 (real)
!                a5_imag(iradius,itheta,itime) = (B3_imag*a3_real_m(iradius,itheta,itime) + B3_real*a3_imag_m(iradius,itheta,itime))  &   ! A3*c3 (imag)
!                                              - (B4_imag*a4_real_m(iradius,itheta,itime) + B4_real*a4_imag_m(iradius,itheta,itime))      ! A4*c4 (imag)
!            end do !itheta
!        end do !iradius

    end subroutine apply_nonreflecting_condition
    !********************************************************************************

    !>  Compute boundary global average by averaging spatio-temporal average over
    !!  radial stations.
    !!
    !!  @author Nathan A. Wukie
    !!  @date   4/16/2018
    !!
    !--------------------------------------------------------------------------------
    subroutine compute_boundary_global_average(self,worker,bc_comm,density_bar,vel1_bar,vel2_bar,vel3_bar,pressure_bar, &
                                                                   density_avg,vel1_avg,vel2_avg,vel3_avg,pressure_avg)
        class(outlet_nrbc_lindblad_t),  intent(inout)   :: self
        type(chidg_worker_t),                       intent(in)      :: worker
        type(mpi_comm),                             intent(in)      :: bc_comm
        type(AD_D),                                 intent(in)      :: density_bar(:)
        type(AD_D),                                 intent(in)      :: vel1_bar(:)
        type(AD_D),                                 intent(in)      :: vel2_bar(:)
        type(AD_D),                                 intent(in)      :: vel3_bar(:)
        type(AD_D),                                 intent(in)      :: pressure_bar(:)
        type(AD_D),                                 intent(inout)   :: density_avg
        type(AD_D),                                 intent(inout)   :: vel1_avg
        type(AD_D),                                 intent(inout)   :: vel2_avg
        type(AD_D),                                 intent(inout)   :: vel3_avg
        type(AD_D),                                 intent(inout)   :: pressure_avg

        real(rk)    :: area, dr
        integer(ik) :: irad

        density_avg  = ZERO*density_bar(1)
        vel1_avg     = ZERO*density_bar(1)
        vel2_avg     = ZERO*density_bar(1)
        vel3_avg     = ZERO*density_bar(1)
        pressure_avg = ZERO*density_bar(1)

        if (worker%coordinate_system() == 'Cartesian') then
            dr = self%r(2) - self%r(1)
            area = ZERO
            do irad = 1,size(density_bar)-1
                density_avg  = density_avg  + dr*(density_bar(irad+1)  + density_bar(irad))/TWO
                vel1_avg     = vel1_avg     + dr*(vel1_bar(irad+1)     + vel1_bar(irad))/TWO
                vel2_avg     = vel2_avg     + dr*(vel2_bar(irad+1)     + vel2_bar(irad))/TWO
                vel3_avg     = vel3_avg     + dr*(vel3_bar(irad+1)     + vel3_bar(irad))/TWO
                pressure_avg = pressure_avg + dr*(pressure_bar(irad+1) + pressure_bar(irad))/TWO
                area = area + dr
            end do

        else if (worker%coordinate_system() == 'Cylindrical') then
            dr = self%r(2) - self%r(1)
            area = ZERO
            do irad = 1,size(density_bar)-1
                density_avg  = density_avg  + dr*(self%r(irad+1)*density_bar(irad+1)  + self%r(irad)*density_bar(irad))/TWO
                vel1_avg     = vel1_avg     + dr*(self%r(irad+1)*vel1_bar(irad+1)     + self%r(irad)*vel1_bar(irad))/TWO
                vel2_avg     = vel2_avg     + dr*(self%r(irad+1)*vel2_bar(irad+1)     + self%r(irad)*vel2_bar(irad))/TWO
                vel3_avg     = vel3_avg     + dr*(self%r(irad+1)*vel3_bar(irad+1)     + self%r(irad)*vel3_bar(irad))/TWO
                pressure_avg = pressure_avg + dr*(self%r(irad+1)*pressure_bar(irad+1) + self%r(irad)*pressure_bar(irad))/TWO
                area = area + dr*(self%r(irad+1)+self%r(irad))/TWO
            end do

        end if

        density_avg  = density_avg  / area
        vel1_avg     = vel1_avg     / area
        vel2_avg     = vel2_avg     / area
        vel3_avg     = vel3_avg     / area
        pressure_avg = pressure_avg / area

    end subroutine compute_boundary_global_average
    !********************************************************************************

end module bc_state_outlet_nrbc_lindblad

max_stars_repo_name: nwukie/ChiDG

max_issues_repo_name: nwukie/ChiDG

max_forks_repo_name: nwukie/ChiDG


text:
!! The following example returns the INCORRECT results with version
!! 18.0.0 20170811 of the Intel compiler:
!!
!! $ ifort intel-bug-20171023.f90
!! $ ./a.out
!!   1065353216 (expect 1)
!!
!! CORRECT results are obtained with versions 17.0.1 and 16.0.3, as well
!! as the gfortran and NAG compilers.

module any_scalar_type

  type any_scalar
    class(*), allocatable :: value
  end type

contains

  function value_ptr(this) result(uptr)
    class(any_scalar), target, intent(in) :: this
    class(*), pointer :: uptr
    uptr => this%value
  end function

end module

program main

  use any_scalar_type

  type(any_scalar), target :: foo
  class(*), allocatable :: scalar
  class(*), pointer :: uptr

  type point; real x, y; end type
  allocate(scalar, source=point(1.0,2.0))

  allocate(foo%value, source=1)
  uptr => value_ptr(foo)
  deallocate(scalar)
  allocate(scalar, source=uptr) ! should hold the value 1

  select type (scalar)
  type is (integer)
    if (scalar == 1) then
      print *, 'pass:', scalar, '(expect 1)'
    else
      print *, 'fail:', scalar, '(expect 1)'
    end if
  end select

end program

max_stars_repo_name: dhnza/fortran-compiler-tests

max_issues_repo_name: dhnza/fortran-compiler-tests

max_forks_repo_name: dhnza/fortran-compiler-tests


	module pH2
	real*8, parameter :: a2bohr = 0.52917721092, hartree_wavenumber = 219474.63
	real*8, parameter :: re = 3.47005, De = 24.2288
	real*8, parameter :: r_ref = 4.60, Vmin = -24.2288

	real*8, parameter :: C6 = 5.820364d04, C8 = 2.87052154d05, C10 = 1.80757343d06
	real*8, parameter, dimension(0:10) :: b = (/-6.631d-02, 1.346d-01, -3.300d-02, 6d0, -1.4d01,  &
				-1.193d02, 2.290d02, 1.110d03, -1.850d03, -3.5d03, 6.0d03/)

	contains

	double precision function damp(r,n)
		implicit real*8(a-h,o-z)
		den = 1.10

		damp = (1.0 - exp(-3.30 * den * r/n - 0.423 * (den * r)**2/sqrt(dble(n))))**(n-1)

	end function

	end module

	program main

	use pH2, only : hartree_wavenumber

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

	real*8,allocatable :: ke(:),V(:),x(:),y(:),px(:),py(:),w(:),c(:)
	real*8,allocatable :: qfx(:),qfy(:),qp(:),fx(:),rx(:),      &
					   fy(:),s(:),dpx(:),dpy(:),s0(:),dr(:),ddr(:), &
					   ap(:),ar(:),ddp(:)
	real*8 :: f(3),err(2)

	real :: gasdev
	complex*16 :: cor,im
	common/smaple/xmin,xmax,dx

	open(100,file='energy.dat')
	open(102,file='xav.out')
	open(103,file='pes.out')
	open(101,file='p.out')
	open(104,file='ke.out')
	open(105,file='weight.out')
	open(106,file='traj.out')
	open(107,file='cor.out')
	open(119,file='err.out')
	open(120,file='px.out')

	open(5,file='IN')
	read(5,*) Ntraj
	read(5,*) kmax,dt
	read(5,*) am
	read(5,*) idum1
	read(5,*) ax
	read(5,*) qx0,p0
	read(5,*) fc
	close(5)

	allocate(ke(Ntraj),V(Ntraj),x(Ntraj))
	allocate(y(Ntraj),px(Ntraj),py(Ntraj),s(Ntraj),w(Ntraj),qfx(Ntraj))
	allocate(qfy(Ntraj),qp(Ntraj),fx(Ntraj),fy(Ntraj),dpx(Ntraj),dpy(Ntraj),s0(Ntraj))
	allocate(c(ntraj),rx(ntraj),dr(ntraj),ddr(ntraj),ap(ntraj),ar(ntraj),ddp(ntraj))

	dt2 = dt/2d0
	t = 0d0
	s = 0d0
	pi = 4.d0*atan(1.d0)
	pow = 6d0

	im = (0d0,1d0)

    ww = 0d0
      do i=1,Ntraj
        x(i) = gasdev(idum1)
        x(i)=x(i)/sqrt(4d0*ax)+qx0
!        y(i)=y(i)/sqrt(4d0*ay)+qy0
!        x(i) = xmin+dx*(i-1)
!        c(i) = -ax*(x(i)-qx0)**2+log(dsqrt(2d0*ax/pi))/2d0
!        w(i) = exp(2d0*c(i))*dx
        w(i) = 1d0/Ntraj
        rx(i) = - 2d0*ax*(x(i)-qx0)
        ww = ww+w(i)
      enddo

      write(*,*) 'Weights = ', ww

       px = p0

! ---- print out the initial conditions
	write(*,*) 'Initial Conditions'
	print *,'ax =',ax
	write(*,*) 'Number of trajectories =', Ntraj
	write(*,*) 'Time steps =',kmax,dt
	write(*,100) am, p0
100	format('Mass = ',f14.7, 'Initial mementum',f14.7)

	s = 0d0
	! --- initial value for action function s
	do i=1,Ntraj
		s(i) = px(i)*(x(i)-x0)
	enddo

! --- print intial trajectories and potential
	xmin = 4
	dx = 0.02
    do i=1,300
		x_grid = xmin + (i-1) * dx
		write(103,10000) x_grid,vpot(x_grid)
    enddo

      write(*,*) 'friction coeff =',fc

    print *,'initial trajectory', minval(x),maxval(x)

          call derivs(x,ntraj,fx,v)
          call fitp(am,x,px,rx,w,Ntraj,dpx,ddp,dr,ddr,qp,ap,ar,err)

    print *,'fx',fx(1:10)
    print *,'dr',dr(1:10)
    print *,'dp',dpx(1:10)

!     begin the time step loop
        do 10 k=1,kmax

          t = t + dt   ! increase t by dt

          do i=1,Ntraj

            px(i) = px(i)+ ( fx(i) + 1d0/2d0/am * (2d0 * rx(i)*dr(i) + ddr(i)))*dt2 - fc * px(i) * dt2
            rx(i) = rx(i) - 1d0/2d0/am * (2.0 * dpx(i) * rx(i) + ddp(i))*dt2

          enddo

        print *,'px',px(1:10)

          do i=1,ntraj
            x(i) = x(i) + px(i) / am * dt
          enddo

        print *,'x',x(1:10)

          call derivs(x,ntraj,fx,v)
          call fitp(am,x,px,rx,w,Ntraj,dpx,ddp,dr,ddr,qp,ap,ar,err)

          do i=1,Ntraj

            px(i) = px(i)+ ( fx(i) - fc * px(i) + 1d0/2d0/am * (2d0 * rx(i)*dr(i) + ddr(i)))*dt2
            rx(i) = rx(i) - 1d0/2d0/am * (2.0 * dpx(i) * rx(i) + ddp(i))*dt2

          enddo

          do i=1,Ntraj
            s(i) = s(i)+(- px(i)**2/2d0/am - v(i)-qp(i))*dt
          enddo

!
!	do  i=1,Ntraj
!          px(i)=px(i)+(-fx(i)-qfx(i))*dt2-px(i)*dpx(i)/m1*dt2
!          py(i)=py(i)+(-fy(i)-qfy(i))*dt2-py(i)*dpy(i)/m2*dt2
!        enddo
! update weights for guassian trajectories
!	do i=1,Ntraj
!	  w(i) = exp(-2d0*(s(i)-s0(i)))*1d0/Ntraj
!	enddo

         write(119,10000) t,(err(i),i=1,2)

! update potential, kinetic, and total energy for each trajectory
         Ek = 0d0
         do i=1,Ntraj
           Ek =  Ek + px(i)**2/(2d0*am) * w(i)
         enddo
! print out a random trajectory

         write(106,10000) t,(x(i),i=1,20)

! calculate the total energy, the sum of all the trajectories
          call aver(Ntraj,w,V,VAve)
          call aver(Ntraj,w,qp,UAVe)
        if(UAve < 0) then
			print *,'Quantum potential = ',UAve, ' Cannot be a negative number'
			stop
		endif

	print *, 'p',px(1:10)

! --- autocorrelation
      cor = (0d0,0d0)
      do i=1,ntraj
        cor = cor+w(i)*exp(2d0*s(i)*im)
      enddo
        write(107,10000) t,cor,abs(cor)

		Ek = Ek * hartree_wavenumber
		VAve = VAve * hartree_wavenumber
		UAve = UAVe * hartree_wavenumber
        Etot = (Ek + VAve + UAve)
        write(100,10000) t, Ek, VAve, UAve, Etot

10	end do

      x1 = 0d0
      x2 = 0d0

      do i=1,ntraj
        x1 = x1 + x(i)*w(i)
        x2 = x2 + x(i)**2*w(i)
      enddo

      write(*,*) 'Expectation of x, var',x1,x2 - x1**2
      write(*,*) 'Total Energy = ', Etot  , 'cm-1.'

      do i=1,ntraj
        write(120,10000) x(i),px(i)
      enddo

10000 format(100(e14.7,1x))
1001  format(A10,20(f10.6,1x))

      close(101)
      close(102)
      close(103)
      close(104)
      close(105)
      close(106)

      return
      end program main

!----------------------------------------------
! average over trajectories y = sum(w(i)**x(i))
!----------------------------------------------
        subroutine aver(Ntraj,w,x,y)
        implicit real*8(a-h,o-z)
        real*8 :: x(Ntraj),w(Ntraj)
        y = 0d0

        do i=1,Ntraj
          y = y+x(i)*w(i)
        enddo

        return
        end subroutine

!     -------------------------------------
!     potential and derivatives
!     -------------------------------------
      subroutine derivs(x,Ntraj,fx,v)
      use pH2
      implicit real*8(a-h,o-z)
      integer*4,intent(in) :: ntraj
      real*8,   intent(in) :: x(Ntraj)
      real*8,   intent(out) :: v(ntraj),fx(ntraj)
      real*8 :: dv(ntraj)
      character*20 :: PES

      PES = 'pH2'

      if(PES == 'Morse') then

      a0 = -0.180516
      a1 = 0.196123
      a2 = 0.95861
      a3 = 1.33249
      a4 = 0.0322998

      do i=1,ntraj
        z = 1d0-exp(-a2*(x(i)-a3))
        v(i)= a0 + a1*z**2 + a4/x(i)**4
        dv(i) = 2d0*a1*z*a2*x(i)*exp(-a2*(x(i)-a3)) - 4d0*a4/x(i)**5
        fx(i) = -dv(i)
      enddo

      elseif(PES == 'pH2') then

        dr = 1.0d-4

        do i = 1, Ntraj

           r = x(i) + re/a2bohr

           v(i) = vpot(r)
           dv(i) = (vpot(r+dr)-vpot(r))/dr    ! finite difference for force
         enddo

         fx = -dv

      endif

      return
      end subroutine

	double precision function vpot(r)
	use pH2
	implicit real*8(a-h,o-z)

	r = r * a2bohr

        beta_inf = log(2.0 * De / u_LR(re))
           s = 0.0
          do j = 0,10
            s = s + b(j) * y_ref(r,1)**j
          enddo

          beta = y_ref(r,6) * beta_inf + (1.0-y_ref(r,6))  * s

    vpot = De*(1.0 - u_LR(r)/u_LR(re) * exp(-beta * y_eq(r,6)))**2

	vpot = vpot + Vmin

    vpot = vpot / hartree_wavenumber

	end

    double precision function y_eq(r,n)
    use pH2, only : re
    implicit real*8(a-h,o-z)

    y_eq = (r**n - re**n)/(r**n + re**n)

    end function

    double precision function y_ref(r,n)
    use pH2, only : r_ref
    implicit real*8(a-h,o-z)

    y_ref = (r**n - r_ref**n)/(r**n + r_ref**n)

    end function

    double precision function u_LR(r)
    use pH2, only : C6, C8, C10, damp
    implicit real*8(a-h,o-z)

	u_LR = damp(r,6) * C6/r**6 + damp(r,8) * C8/r**8 + damp(r,10) * C10 / r**10

    end function


text:
      SUBROUTINE GOUT(SA,OMEG,TIME,IPLOT,IPLTPE,PHI,TESTY)
C
C
C     SUBROUTINE GOUT DETERMINES EULER ANGLES AND RATES WRT THE
C     LOCAL VERTICAL FRAME FOR GRAVITY GRADIENT SATELLITES.
C     IT CALLS ARANGE TO EXTRACT ELEMENT DATA ; DMOMNT TO CALCULATE
C     ELEMENT BENDING MOMENTS AND TENSIONS ; AND WRTPLT TO WRITE
C     TEMPORARY DATA SET FOR PLOTS .
C     MODIFICATION AS OF SEPT. 1977.
C
C
      IMPLICIT REAL*8 (A-H,O-Z)
C
      COMMON/CANTNA/ A(10,3),ADOT(10,3),B(10,3),BDOT(10,3),DIN(10,3),
     .               DINDOT(10,3),DOUT(10,3),DOUTDT(10,3)
C
C
      COMMON/COMBCR/ OMBC1(3)
C
      COMMON/CONSTS/ PI,TWOPI,RADIAN
C
      COMMON/CSTAT /XDEB(20)
C
      COMMON/DMMNT1/ZKBM(6),EMAK(10),EMBK(10),DUMM1(6),
     *    ITORK,IBENDM,ITENSE,ITNS1
C
      COMMON/HAMOUT/HAMILT,IHAMLT
C
      COMMON/HOUTPT/IHCALC
C
      COMMON/HVCOUT/HBODY(3),HINERT(3),HMAG
C
      COMMON/ICNTRL/KNTRL(10)
C
      COMMON/IPOOL1/ IGRAV,IDAMP,IK,K1,ITIM,IAB,IAPS,IBB,IBPS,NK(10),
     .               LK(10),LLK(10)
C
      COMMON/JBTEST/ IBTEST
C
      COMMON/KNERGY/ S,OMC,BT
C
      COMMON/NUTANG/ANUT
C
C
      COMMON/RMAIN1/ DELTAT,FACTOR,FREQ,TSTOP,DELMIT,
     .               UPBND(150),DNBND(150)
C
      COMMON/RPOOL2/ PO,SD(3),DAN(3,10),DBN(3,10),CFMT(3,3),DIY1(3),
     .               SD1(3),DT1,P1,AERO,DTO1,YIZK(3),PO1
C
      COMMON/SUNANG/ DSL2
C
      COMMON/VECTRS/ XSAT(3),XSATDT(3),AD(3)
C
      COMMON/XIN2  / ALFAE,BETAE,GAMAE,OMBC(3),ITEST2
C
C
      DIMENSION SA(3,3),OMEG(3),SAT(3,3),S(3,3),BS(3,3),BT(3,3),
     .          OMC(3),SXOMC(3)
C
C
      DATA IB /0/
C
      DATA ALPOLD,BETOLD,GAMOLD/ 3*0.D0/
C
C
      JTEST=1
      ITEST2=1
C
C
      RADIUS=0.D0
      DO 5 I=1,3
    5 RADIUS=RADIUS + XSAT(I)*XSAT(I)
      RADIUS=DSQRT(RADIUS)
C
      DO 10 I=1,3
   10 BS(3,I)=XSAT(I)/RADIUS
C
      YZ=XSAT(2)*XSATDT(3) - XSAT(3)*XSATDT(2)
      ZX=XSAT(3)*XSATDT(1) - XSAT(1)*XSATDT(3)
      XY=XSAT(1)*XSATDT(2) - XSAT(2)*XSATDT(1)
C
      RXVMAG=DSQRT(YZ*YZ + ZX*ZX + XY*XY)
C
      BS(2,1)=YZ/RXVMAG
      BS(2,2)=ZX/RXVMAG
      BS(2,3)=XY/RXVMAG
C
      BS(1,1)=BS(2,2)*BS(3,3)-BS(2,3)*BS(3,2)
      BS(1,2)=BS(2,3)*BS(3,1)-BS(2,1)*BS(3,3)
      BS(1,3)=BS(2,1)*BS(3,2)-BS(2,2)*BS(3,1)
C
      DO 20 I=1,3
      DO 20 L=1,3
   20 BT(I,L)=BS(L,I)
      SDLV2=BT(1,2)*SD(1)+BT(2,2)*SD(2)+BT(3,2)*SD(3)
      DELSLV=DARCOS(-1.0D0*SDLV2)
      DSL2=DELSLV/RADIAN
      DO 30 I=1,3
      DO 30 J=1,3
   30 SAT(I,J)=SA(J,I)
C
      DO 40 I=1,3
      DO 40 J=1,3
      S(I,J)=0.D0
      DO 40 K=1,3
   40 S(I,J)=S(I,J) + SAT(I,K)*BT(K,J)
C
C
      IF(S(3,1).EQ.0.D0.AND.S(3,3).EQ.0.D0) GO TO 50
      W=-S(3,2)
      ALFAE=DARSIN(W)
      IF(DABS(ALFAE/RADIAN).GT.89.D0.AND.DABS(ALFAE/RADIAN).LT.91.D0)
     .           GO TO 100
      SN=DCOS(ALFAE)
      TEST1=S(3,3)/SN
      W=S(3,1)/SN
      BETAE=DARSIN(W)
      IF(TEST1.LT.0.D0) BETAE=-DARSIN(W) + PI
      TEST2=S(2,2)/SN
      W=S(1,2)/SN
      GAMAE=DARSIN(W)
      IF(TEST2.LT.0.D0) GAMAE=-DARSIN(W) + PI
      ALFAE=ALFAE/RADIAN
      BETAE=BETAE/RADIAN
      GAMAE=GAMAE/RADIAN
      ALPOLD=ALFAE
      BETOLD=BETAE
      GAMOLD=GAMAE
      IF (IBTEST .EQ. 0)  GO TO 110
      IF (IB .EQ. 0)  BPAST=BETAE
      IB=1
      IF (BPAST*BETAE .GE. 0.D0)  GO TO 101
      IF (BETAE .LT. 0.D0)  GO TO 101
      TESTY=TIME
      IBTEST=0
      IB=0
      GO TO 110
  101 BPAST=BETAE
      GO TO 110
   50 IF(S(3,2)-1.D0) 60,80,60
   60 IF(S(3,2)+1.D0) 100,70,100
   70 ALFAE=90.D0
      ALPOLD=ALFAE
      GO TO 100
   80 ALFAE=-90.D0
      ALPOLD=ALFAE
  100 ITEST2=0
  110 OMC(1)=0.D0
      OMC(2)=RXVMAG/RADIUS**2
C
      DOTX=AD(1)*YZ
      DOTY=AD(2)*ZX
      DOTZ=AD(3)*XY
C
      OMC(3)=((DOTX + DOTY + DOTZ)*RADIUS)/RXVMAG**2
C
      DO 120 I=1,3
      SXOMC(I)=0.D0
      DO 120 J=1,3
  120 SXOMC(I)=SXOMC(I) + S(I,J)*OMC(J)
C
      DO 130 I=1,3
      OMBC(I)=(OMEG(I) - SXOMC(I))/RADIAN
      OMBC1(I)=OMC(I)
  130 OMC(I)=OMBC(I)*RADIAN
C
      ALFAE=ALPOLD
      BETAE=BETOLD
      GAMAE=GAMOLD
      PHILD=0.0D0
      IF(IDAMP.NE.0) PHILD=PHI/RADIAN
      AL=ALFAE*RADIAN
      BE=BETAE*RADIAN
      GA=GAMAE*RADIAN
      RRAT=OMBC(1)*DCOS(GA)-OMBC(2)*DSIN(GA)
      PRAT=(OMBC(1)*DSIN(GA)+OMBC(2)*DCOS(GA))/DCOS(AL)
      YRAT=OMBC(3)+DTAN(AL)*(OMBC(1)*DSIN(GA)+OMBC(2)*DCOS(GA))
C
      IF(KNTRL(1).EQ.0) GO TO 140
C
      IF(IHCALC.EQ.0) GO TO 140
C
      WS=HBODY(2)/HMAG
      IF(WS.GT.1.0D0) WS=1.0D0
      IF(WS.LT.-1.0D0) WS=-1.0D0
      ANUT=DARCOS(WS)/RADIAN
      XDEB(5)=ANUT
C
  140 CONTINUE
C
C     CALL ARANGE
C
      IF(ITENSE.NE.0.OR.IBENDM.NE.0) CALL DMOMNT
C
      IF(IHAMLT.NE.0) CALL ENERGY(2,JTEST)
C
      IF(IPLOT.NE.0) CALL WRTPLT(TIME,ALFAE,BETAE,GAMAE,OMBC,PHILD
     1                          ,RRAT,PRAT,YRAT,2)
C
      RETURN
      END

ext: for


text:
      subroutine setunits(numfil)
!***********************************************************************
!  Copyright, 1993, 2004,  The  Regents of the University of California.
!  This program was prepared by the Regents of the University of
!  California at Los Alamos National Laboratory (the University) under
!  contract No. W-7405-ENG-36 with the U.S. Department of Energy (DOE).
!  All rights in the program are reserved by the DOE and the University.
!  Permission is granted to the public to copy and use this software
!  without charge, provided that this Notice and any statement of
!  authorship are reproduced on all copies. Neither the U.S. Government
!  nor the University makes any warranty, express or implied, or
!  assumes any liability or responsibility for the use of this software.
C***********************************************************************
CD1
CD1 PURPOSE
CD1
CD1 Set file unit numbers.
CD1
C***********************************************************************
CD2
CD2 REVISION HISTORY
CD2
CD2 Revision                    ECD
CD2 Date         Programmer     Number  Comments
CD2
CD2 25-OCT-93    Z. Dash        22      Add prolog.
CD2              G. Zyvoloski           Initial implementation.
CD2
CD2 $Log:   /pvcs.config/fehm90/src/setunits.f_a  $
!D2
!D2    Rev 2.5   06 Jan 2004 10:43:54   pvcs
!D2 FEHM Version 2.21, STN 10086-2.21-00, Qualified October 2003
!D2
!D2    Rev 2.4   29 Jan 2003 09:15:18   pvcs
!D2 FEHM Version 2.20, STN 10086-2.20-00
!D2
!D2    Rev 2.3   14 Nov 2001 13:12:34   pvcs
!D2 FEHM Version 2.12, STN 10086-2.12-00
!D2
!D2    Rev 2.2   06 Jun 2001 13:26:36   pvcs
!D2 FEHM Version 2.11, STN 10086-2.11-00
!D2
!D2    Rev 2.1   30 Nov 2000 12:07:36   pvcs
!D2 FEHM Version 2.10, STN 10086-2.10-00
!D2
!D2    Rev 2.0   Fri May 07 14:45:24 1999   pvcs
!D2 FEHM Version 2.0, SC-194 (Fortran 90)
CD2
CD2    Rev 1.2   Thu Feb 01 16:08:14 1996   hend
CD2 Updated Requirements Traceability
CD2
CD2    Rev 1.1   03/18/94 15:55:44   gaz
CD2 Added solve_new and cleaned up memory management.
CD2
CD2    Rev 1.0   01/20/94 10:27:32   pvcs
CD2 original version in process of being certified
CD2
C***********************************************************************
CD3
CD3 INTERFACES
CD3
CD3 None
CD3
C***********************************************************************
CD4
CD4 GLOBAL OBJECTS
CD4
CD4 Global Constants
CD4
CD4   None
CD4
CD4 Global Types
CD4
CD4   None
CD4
CD4 Global Variables
CD4
CD4                            COMMON
CD4   Identifier      Type     Block  Description
CD4
CD4   incoor          INT      faai   Unit number for coordinate input file
CD4   inpt            INT      faai   Unit number for input file
CD4   inzone          INT      faai   Unit number for zone input file
CD4   iout            INT      faai   Unit number for output file
CD4   isave           INT      faai   Unit number for restart file (to write)
CD4   ischk           INT      faai   Unit number for input data check file
CD4   iscon           INT      faai   Unit number for contour data file
CD4   iscon1          INT      faai   Unit number for dual porosity or dpdp
CD4                                     contour data file
CD4   ishis           INT      faai   Unit number for history data file
CD4   isstor          INT      faai   Unit number for element coefficient file
CD4   istrc           INT      faai   Unit number for tracer history data file
CD4   iowork          CHAR     faay   File use/type designator
CD4   isave           INT      faai   Unit number for restart file (to write)
CD4
CD4 Global Subprograms
CD4
CD4   None
CD4
C***********************************************************************
CD5
CD5 LOCAL IDENTIFIERS
CD5
CD5 Local Constants
CD5
CD5   None
CD5
CD5 Local Types
CD5
CD5   None
CD5
CD5 Local variables
CD5
CD5   Identifier      Type     Description
CD5
CD5   numfil          INT ARY  I/O file unit numbers
CD5
CD5 Local Subprograms
CD5
CD5   None
CD5
C***********************************************************************
CD6
CD6 FUNCTIONAL DESCRIPTION
CD6
C***********************************************************************
CD7
CD7 ASSUMPTIONS AND LIMITATIONS
CD7
CD7 None
CD7
C***********************************************************************
CD8
CD8 SPECIAL COMMENTS
CD8
CD8  Requirements from SDN: 10086-RD-2.20-00
CD8    SOFTWARE REQUIREMENTS DOCUMENT (RD) for the
CD8    FEHM Application Version 2.20
CD8
C***********************************************************************
CD9
CD9 REQUIREMENTS TRACEABILITY
CD9
CD9 2.6 Provide Input/Output Data Files
CD9 3.0 INPUT AND OUTPUT REQUIREMENTS
CD9
C***********************************************************************
CDA
CDA REFERENCES
CDA
CDA None
CDA
C***********************************************************************
CPS
CPS PSEUDOCODE
CPS
CPS BEGIN setunits
CPS
CPS   Assign file unit number for each I/O file [iocntl, inpt, incoor, . . .
CPS   . . . inzone, iout, iread, isave, ishis, istrc, iscon, iscon1, . . .
CPS   . . . isstor, ischk]
CPS
CPS END setunits
CPS
C***********************************************************************

      use comai
      implicit none

      integer numfil(*)

C **** Set file unit numbers ****
      iocntl = numfil(1)
      inpt = numfil(2)
      incoor = numfil(3)
      inzone = numfil(4)
      iout = numfil(5)
      iread = numfil(6)
      isave = numfil(7)
      ishis = numfil(8)
      istrc = numfil(9)
      iscon = numfil(10)
      iscon1 = numfil(11)
      isstor = numfil(12)
      ischk = numfil(13)

      end

max_stars_repo_name: satkarra/FEHM

max_issues_repo_name: satkarra/FEHM

max_forks_repo_name: satkarra/FEHM


!
!  Copyright 2017 SALMON developers
!
!  Licensed under the Apache License, Version 2.0 (the "License");
!  you may not use this file except in compliance with the License.
!  You may obtain a copy of the License at
!
!      http://www.apache.org/licenses/LICENSE-2.0
!
!  Unless required by applicable law or agreed to in writing, software
!  distributed under the License is distributed on an "AS IS" BASIS,
!  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
!  See the License for the specific language governing permissions and
!  limitations under the License.
!
subroutine calcJxyz_all
use salmon_communication, only: comm_is_root, comm_summation
use scf_data
use allocate_psl_sub
implicit none
integer :: iatom,ix,iy,iz
  integer :: i,j
  integer :: mmx(mg_num(1)*mg_num(2)*mg_num(3))
  integer :: mmy(mg_num(1)*mg_num(2)*mg_num(3))
  integer :: mmz(mg_num(1)*mg_num(2)*mg_num(3))
  integer :: lx(lg_num(1)*lg_num(2)*lg_num(3))
  integer :: ly(lg_num(1)*lg_num(2)*lg_num(3))
  integer :: lz(lg_num(1)*lg_num(2)*lg_num(3))
  real(8) :: alx,aly,alz
  real(8) :: hx,hy,hz

  hx=Hgs(1)
  hy=Hgs(2)
  hz=Hgs(3)
  alx=Hgs(1)*dble(lg_num(1))
  aly=Hgs(2)*dble(lg_num(2))
  alz=Hgs(3)*dble(lg_num(3))

  do iz=mg_sta(3),mg_end(3)
  do iy=mg_sta(2),mg_end(2)
  do ix=mg_sta(1),mg_end(1)
    i=(iz-mg_sta(3))*mg_num(1)*mg_num(2)+(iy-mg_sta(2))*mg_num(1)+ix-mg_sta(1)+1
    mmx(i)=ix
    mmy(i)=iy
    mmz(i)=iz
  end do
  end do
  end do

  do iz=lg_sta(3),lg_end(3)
  do iy=lg_sta(2),lg_end(2)
  do ix=lg_sta(1),lg_end(1)
    i=(iz-lg_sta(3))*lg_num(1)*lg_num(2)+(iy-lg_sta(2))*lg_num(1)+ix-lg_sta(1)+1
    lx(i)=ix
    ly(i)=iy
    lz(i)=iz
  end do
  end do
  end do

  call calc_mps(pp,ppg,alx,aly,alz,lx,ly,lz,lg_num(1)*lg_num(2)*lg_num(3),   &
                                   mmx,mmy,mmz,mg_num(1)*mg_num(2)*mg_num(3),   &
                                   hx,hy,hz)
  call calc_mps(pp,ppg_all,alx,aly,alz,lx,ly,lz,lg_num(1)*lg_num(2)*lg_num(3),   &
                                       lx,ly,lz,lg_num(1)*lg_num(2)*lg_num(3),   &
                                       hx,hy,hz)
  Mps(1:MI)=ppg%mps(1:MI)
  Mps_all(1:MI)=ppg_all%mps(1:MI)

  call init_jxyz(ppg)
  call init_jxyz(ppg_all)

  call calc_jxyz(pp,ppg,alx,aly,alz,lx,ly,lz,lg_num(1)*lg_num(2)*lg_num(3),   &
                                    mmx,mmy,mmz,mg_num(1)*mg_num(2)*mg_num(3),   &
                                    hx,hy,hz)
  call calc_jxyz(pp,ppg_all,alx,aly,alz,lx,ly,lz,lg_num(1)*lg_num(2)*lg_num(3),   &
                                    lx,ly,lz,lg_num(1)*lg_num(2)*lg_num(3),   &
                                    hx,hy,hz)

  Jxyz=0
  Jxyz_all=0
  Jxxyyzz=0
  Jxxyyzz_all=0

  do iatom=1,MI
    do j=1,Mps(iatom)
      Jxyz(1,j,iatom)=ppg%jxyz(1,j,iatom)
      Jxyz(2,j,iatom)=ppg%jxyz(2,j,iatom)
      Jxyz(3,j,iatom)=ppg%jxyz(3,j,iatom)
    end do
    do j=1,Mps_all(iatom)
      Jxyz_all(1,j,iatom)=ppg_all%jxyz(1,j,iatom)
      Jxyz_all(2,j,iatom)=ppg_all%jxyz(2,j,iatom)
      Jxyz_all(3,j,iatom)=ppg_all%jxyz(3,j,iatom)
    end do
  end do

return

end subroutine calcJxyz_all


text:
module init_vel_module

  use bl_types
  use bl_constants_module
  use bc_module
  use multifab_physbc_module
  use define_bc_module
  use multifab_module
  use fill_3d_module
  use eos_module
  use variables
  use network
  use geometry
  use ml_layout_module
  use ml_cc_restriction_module
  use multifab_fill_ghost_module

  implicit none

  private
  public :: initveldata

contains

  subroutine initveldata(u,s0_init,p0_init,dx,bc,mla)

    use mt19937_module

    type(multifab) , intent(inout) :: u(:)
    real(kind=dp_t), intent(in   ) :: s0_init(:,0:,:)
    real(kind=dp_t), intent(in   ) :: p0_init(:,0:)
    real(kind=dp_t), intent(in   ) :: dx(:,:)
    type(bc_level) , intent(in   ) :: bc(:)
    type(ml_layout), intent(inout) :: mla

    real(kind=dp_t), pointer:: uop(:,:,:,:)
    integer :: lo(mla%dim),hi(mla%dim),ng,dm
    integer :: i,j,k,n,nlevs

    ! random numbers between -1 and 1
    real(kind=dp_t) :: alpha(3,3,3), beta(3,3,3), gamma(3,3,3)

    ! random numbers between 0 and 2*pi
    real(kind=dp_t) :: phix(3,3,3), phiy(3,3,3), phiz(3,3,3)

    ! L2 norm of k
    real(kind=dp_t) :: normk(3,3,3)

    ! random number
    real(kind=dp_t) :: rand

    dm = mla%dim
    nlevs = mla%nlevel

    ng = nghost(u(1))

    ! load in random numbers alpha, beta, gamma, phix, phiy, and phiz
    if (dm .eq. 3 .and. spherical .eq. 1) then
       call init_genrand(20908)
       do i=1,3
          do j=1,3
             do k=1,3
                rand = genrand_real1()
                rand = 2.0d0*rand - 1.0d0
                alpha(i,j,k) = rand
                rand = genrand_real1()
                rand = 2.0d0*rand - 1.0d0
                beta(i,j,k) = rand
                rand = genrand_real1()
                rand = 2.0d0*rand - 1.0d0
                gamma(i,j,k) = rand
                rand = genrand_real1()
                rand = 2.0d0*M_PI*rand
                phix(i,j,k) = rand
                rand = genrand_real1()
                rand = 2.0d0*M_PI*rand
                phiy(i,j,k) = rand
                rand = genrand_real1()
                rand = 2.0d0*M_PI*rand
                phiz(i,j,k) = rand
             enddo
          enddo
       enddo

       ! compute the norm of k
       do i=1,3
          do j=1,3
             do k=1,3
                normk(i,j,k) = sqrt(dble(i)**2+dble(j)**2+dble(k)**2)
             enddo
          enddo
       enddo
    end if

    do n=1,nlevs

       do i = 1, nfabs(u(n))
          uop => dataptr(u(n),i)
          lo =  lwb(get_box(u(n),i))
          hi =  upb(get_box(u(n),i))
          select case (dm)
          case (2)
             call bl_error('initveldata_2d not written')
          case (3)
             if (spherical .eq. 1) then
                call initveldata_3d_sphr(uop(:,:,:,:), lo, hi, ng, dx(n,:), &
                                         s0_init(1,:,:), p0_init(1,:), &
                                         alpha, beta, gamma, phix, phiy, phiz, normk)
             else
                call bl_error('initveldata_3d not written')
             end if
          end select
       end do

    enddo

    if (nlevs .eq. 1) then

       ! fill ghost cells for two adjacent grids at the same level
       ! this includes periodic domain boundary ghost cells
       call multifab_fill_boundary(u(nlevs))

       ! fill non-periodic domain boundary ghost cells
       call multifab_physbc(u(nlevs),1,1,dm,bc(nlevs))

    else

       ! the loop over nlevs must count backwards to make sure the finer grids are done first
       do n=nlevs,2,-1

          ! set level n-1 data to be the average of the level n data covering it
          call ml_cc_restriction(u(n-1),u(n),mla%mba%rr(n-1,:))

          ! fill level n ghost cells using interpolation from level n-1 data
          ! note that multifab_fill_boundary and multifab_physbc are called for
          ! both levels n-1 and n
          call multifab_fill_ghost_cells(u(n),u(n-1),ng,mla%mba%rr(n-1,:), &
                                         bc(n-1),bc(n),1,1,dm)
       enddo

    end if

  end subroutine initveldata

  ! the velocity is initialized to zero plus a perturbation which is a
  ! summation of 27 fourier modes with random amplitudes and phase
  ! shifts over a square of length "velpert_scale".  The parameter
  ! "velpert_radius" is the cutoff radius, such that the actual
  ! perturbation that gets added to the initial velocity decays to
  ! zero quickly if r > "velpert_radius", via the tanh function.
  ! The steepness of the cutoff is controlled by "velpert_steep".  The
  ! relative amplitude of the modes is controlled by
  ! "velpert_amplitude".
  subroutine initveldata_3d_sphr(u,lo,hi,ng,dx,s0_init,p0_init, &
                                 alpha,beta,gamma,phix,phiy,phiz,normk)

    use probin_module, only: prob_lo, prob_hi, &
         velpert_amplitude, velpert_radius, velpert_steep, velpert_scale

    integer, intent(in) :: lo(:), hi(:), ng
    real (kind = dp_t), intent(out) :: u(lo(1)-ng:,lo(2)-ng:,lo(3)-ng:,:)
    real (kind = dp_t), intent(in ) :: dx(:)
    real(kind=dp_t), intent(in   ) ::    s0_init(0:,:)
    real(kind=dp_t), intent(in   ) ::    p0_init(0:)

    ! random numbers between -1 and 1
    real(kind=dp_t), intent(in) :: alpha(3,3,3), beta(3,3,3), gamma(3,3,3)

    ! random numbers between 0 and 2*pi
    real(kind=dp_t), intent(in) :: phix(3,3,3), phiy(3,3,3), phiz(3,3,3)

    ! L2 norm of k
    real(kind=dp_t), intent(in) :: normk(3,3,3)

    ! Local variables
    integer :: i, j, k
    integer :: iloc, jloc, kloc

    ! cos and sin of (2*pi*kx/L + phix), etc
    real(kind=dp_t) :: cx(3,3,3), cy(3,3,3), cz(3,3,3)
    real(kind=dp_t) :: sx(3,3,3), sy(3,3,3), sz(3,3,3)

    ! location of center of star
    real(kind=dp_t) :: xc(3)

    ! radius, or distance, to center of star
    real(kind=dp_t) :: rloc

    ! the point we're at
    real(kind=dp_t) :: xloc(3)

    ! perturbational velocity to add
    real(kind=dp_t) :: upert(3)

    ! initialize the velocity to zero everywhere
    u = ZERO

    ! define where center of star is
    ! this currently assumes the star is at the center of the domain
    xc(1) = 0.5d0*(prob_lo(1)+prob_hi(1))
    xc(2) = 0.5d0*(prob_lo(2)+prob_hi(2))
    xc(3) = 0.5d0*(prob_lo(3)+prob_hi(3))

    ! now do the big loop over all points in the domain
    do iloc = lo(1),hi(1)
       do jloc = lo(2),hi(2)
          do kloc = lo(3),hi(3)

             ! set perturbational velocity to zero
             upert = ZERO

             ! compute where we physically are
             xloc(1) = prob_lo(1) + (dble(iloc)+0.5d0)*dx(1)
             xloc(2) = prob_lo(2) + (dble(jloc)+0.5d0)*dx(2)
             xloc(3) = prob_lo(3) + (dble(kloc)+0.5d0)*dx(3)

             ! compute distance to the center of the star
             rloc = ZERO
             do i=1,3
                rloc = rloc + (xloc(i) - xc(i))**2
             enddo
             rloc = sqrt(rloc)

             ! loop over the 27 combinations of fourier components
             do i=1,3
                do j=1,3
                   do k=1,3
                      ! compute cosines and sines
                      cx(i,j,k) = cos(2.0d0*M_PI*dble(i)*xloc(1)/velpert_scale + phix(i,j,k))
                      cy(i,j,k) = cos(2.0d0*M_PI*dble(j)*xloc(2)/velpert_scale + phiy(i,j,k))
                      cz(i,j,k) = cos(2.0d0*M_PI*dble(k)*xloc(3)/velpert_scale + phiz(i,j,k))
                      sx(i,j,k) = sin(2.0d0*M_PI*dble(i)*xloc(1)/velpert_scale + phix(i,j,k))
                      sy(i,j,k) = sin(2.0d0*M_PI*dble(j)*xloc(2)/velpert_scale + phiy(i,j,k))
                      sz(i,j,k) = sin(2.0d0*M_PI*dble(k)*xloc(3)/velpert_scale + phiz(i,j,k))
                   enddo
                enddo
             enddo

             ! loop over the 27 combinations of fourier components
             do i=1,3
                do j=1,3
                   do k=1,3
                      ! compute contribution from perturbation velocity from each mode
                      upert(1) = upert(1) + &
                           (-gamma(i,j,k)*dble(j)*cx(i,j,k)*cz(i,j,k)*sy(i,j,k) &
                             +beta(i,j,k)*dble(k)*cx(i,j,k)*cy(i,j,k)*sz(i,j,k)) &
                            / normk(i,j,k)

                      upert(2) = upert(2) + &
                           (gamma(i,j,k)*dble(i)*cy(i,j,k)*cz(i,j,k)*sx(i,j,k) &
                           -alpha(i,j,k)*dble(k)*cx(i,j,k)*cy(i,j,k)*sz(i,j,k)) &
                            / normk(i,j,k)

                      upert(3) = upert(3) + &
                           ( -beta(i,j,k)*dble(i)*cy(i,j,k)*cz(i,j,k)*sx(i,j,k) &
                            +alpha(i,j,k)*dble(j)*cx(i,j,k)*cz(i,j,k)*sy(i,j,k)) &
                            / normk(i,j,k)
                   enddo
                enddo
             enddo

             ! apply the cutoff function to the perturbational velocity
             do i=1,3
                upert(i) = velpert_amplitude*upert(i) &
                     *(0.5d0+0.5d0*tanh((velpert_radius-rloc)/velpert_steep))
             enddo

             ! add perturbational velocity to background velocity
             do i=1,3
                u(iloc,jloc,kloc,i) = u(iloc,jloc,kloc,i) + upert(i)
             enddo

          enddo
       enddo
    enddo

  end subroutine initveldata_3d_sphr

end module init_vel_module

max_stars_repo_name: sailoridy/MAESTRO

max_issues_repo_name: sailoridy/MAESTRO

max_forks_repo_name: sailoridy/MAESTRO


text:
The City of Davis had a program where they would loan you an electric vehicles electric car for a week after taking a class and meeting various requirements (e.g. insurance, parking). The cars were part of the citys GEM car fleet. The loaner program was discontinued because of declining interest and vehicle maintenance issues.

It is discussed in an old (Nov 3, 2003) California Aggie http://www.californiaaggie.com/news/2003/11/03/FrontPage/CityPromotes.Use.Of.Electric.Vehicles.Through.Lending.Program1315681.shtml article and on the citys http://cityofdavis.org/pcs/about/fleetsvcs.cfm fleet services page.

In the spring of 2008 the city http://www.cityofdavis.org/pdfs/GEM_auction_Release1.pdf auctioned off some of these vehicles.