!
!  Copyright 2019-2020 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.
!
!--------10--------20--------30--------40--------50--------60--------70--------80--------90--------100-------110-------120-------130
module em_field
  implicit none

contains

!===================================================================================================================================

subroutine calc_emfields(itt,nspin,curr_in,rt)
  use structures, only : s_rt
  use math_constants, only : pi
  use phys_constants, only: cspeed_au
  use salmon_global, only : dt,trans_longi,film_thickness,epsilon_em
  implicit none
  integer   ,intent(in)    :: itt,nspin
  real(8)   ,intent(in)    :: curr_in(3,nspin)
  type(s_rt),intent(inout) :: rt
  !
  integer :: j
  real(8) :: n1,n2
  integer,parameter :: m=100

! electric field
  rt%E_ext(:,itt) = -( rt%Ac_ext(:,itt) - rt%Ac_ext(:,itt-1) )/dt
  rt%E_ind(:,itt) = -( rt%Ac_ind(:,itt) - rt%Ac_ind(:,itt-1) )/dt
  rt%E_tot(:,itt) = -( rt%Ac_tot(:,itt) - rt%Ac_tot(:,itt-1) )/dt

! current density
  if(nspin==1) then
    rt%curr(1:3,itt) = curr_in(1:3,1)
  else if(nspin==2) then
    rt%curr(1:3,itt) = curr_in(1:3,1) + curr_in(1:3,2)
  end if

! vector potential for next step
  if(trans_longi=="lo")then
    rt%Ac_ind(:,itt+1) = 2d0*rt%Ac_ind(:,itt) -rt%Ac_ind(:,itt-1) -4d0*Pi*rt%curr(:,itt)*dt**2
  else if(trans_longi=="tr")then
    rt%Ac_ind(:,itt+1) = 0d0
  else if(trans_longi=="2d") then
    n1 = sqrt(epsilon_em(1)) ! refractive index of medium 1
    n2 = sqrt(epsilon_em(2)) ! refractive index of medium 2
    ! Ac_ind(:,itt+1) = A_ref(t+dt)/c, Ac_tot(:,itt+1) = A_trans(t+dt)/c
    do j=1,m
      rt%Ac_ind(:,itt+1) = rt%Ac_ind(:,0) + (dt/m)*( &
      & ((n1-n2)/(n1+n2))* ( rt%Ac_ext(:,itt+1) - rt%Ac_ext(:,itt) )/dt &
      & - (4d0*pi/(cspeed_au*(n1+n2)))* film_thickness* rt%curr(:,itt) )
      rt%Ac_ind(:,0) = rt%Ac_ind(:,itt+1)
    end do
  else
    stop "error: invalid trans_longi"
  end if
  rt%Ac_tot(:,itt+1) = rt%Ac_ext(:,itt+1) + rt%Ac_ind(:,itt+1)

end subroutine calc_emfields

!===================================================================================================================================

Subroutine calc_Ac_ext(t,Ac_ext)
  real(8),intent(in) :: t
  real(8)            :: Ac_ext(1:3)
  real(8)            :: Ac_ext_t(1:3, 0:0)

  call calc_Ac_ext_t(t, 0d0, 0, 0, Ac_ext_t)
  Ac_ext(1:3) = Ac_ext_t(1:3, 0)
  return
end subroutine calc_Ac_ext

!===================================================================================================================================

Subroutine calc_Ac_ext_t(t0, delta_t, is, ie, Ac_ext_t)
  use math_constants,only: zi,pi
  use salmon_global, only: I_wcm2_1,I_wcm2_2,E_amplitude1,E_amplitude2,ae_shape1,ae_shape2, &
                         & epdir_re1,epdir_re2,epdir_im1,epdir_im2,tw1,tw2,t1_start,omega1,omega2, &
                         & phi_CEP1,phi_CEP2,T1_T2,e_impulse,file_input1
  implicit none
  real(8),intent(in) :: t0
  real(8),intent(in) :: delta_t
  integer,intent(in) :: is
  integer,intent(in) :: ie
  real(8),intent(out) :: Ac_ext_t(1:3, is:ie)
  !
  integer :: i,npower
  real(8) :: t,f0_1,f0_2,tt,T1_T2_tmp

  Ac_ext_t(:,:) = 0d0
!  if(t0 < 0d0) return

  if(I_wcm2_1 < 0d0)then
    f0_1 = E_amplitude1
  else
    f0_1=5.338d-9*sqrt(I_wcm2_1)      ! electric field in a.u.
  end if
  if(I_wcm2_2 < 0d0)then
    f0_2 = E_amplitude2
  else
    f0_2=5.338d-9*sqrt(I_wcm2_2)      ! electric field in a.u.
  end if

  T1_T2_tmp = T1_T2

  select case(AE_shape1)
  case('impulse')
    do i=is, ie
      t=t0+i*delta_t
      if (0d0 <= t) then
        Ac_ext_t(1,i) = epdir_re1(1)*e_impulse
        Ac_ext_t(2,i) = epdir_re1(2)*e_impulse
        Ac_ext_t(3,i) = epdir_re1(3)*e_impulse
      end if
    end do

  case('Acos2','Acos3','Acos4','Acos6','Acos8')

    select case(ae_shape1)
    case('Acos2'); npower = 2
    case('Acos3'); npower = 3
    case('Acos4'); npower = 4
    case('Acos6'); npower = 6
    case('Acos8'); npower = 8
    case default
      stop 'Error in init_rt.f90'
    end select

    do i=is, ie
      t=t0+i*delta_t
      tt = t - 0.5d0*tw1 - t1_start
      if (abs(tt)<0.5d0*tw1) then
        Ac_ext_t(:,i) = -f0_1/omega1*(cos(pi*tt/tw1))**npower &
          *aimag( (epdir_re1(:) + zI*epdir_im1(:)) &
          *exp(zI*(omega1*tt+phi_CEP1*2d0*pi))  &
          )
      end if
    end do
    T1_T2_tmp = T1_T2 + t1_start

  case('Ecos2')

    if(phi_CEP1 /= 0.75d0)then
      stop "Error: phi_cep1 should be 0.75 when ae_shape1 is 'Ecos2'."
    end if
    if(sum(abs(epdir_im1(:)))>1.0d-8)then
      stop "Error: ae_shape1 should be 'Acos2' when epdir_im1 is used."
    end if
    do i=is, ie
      t=t0+i*delta_t
      tt = t - 0.5d0*tw1 - t1_start
      if (abs(tt)<0.5d0*tw1) then
        Ac_ext_t(:,i) = -epdir_re1(:)*f0_1/(8d0*pi**2*omega1 - 2d0*tw1**2*omega1**3) &
          *( &
          (-4d0*pi**2+tw1**2*omega1**2 + tw1**2*omega1**2*cos(2d0*pi*tt/tw1))*cos(omega1*tt) &
          +2d0*pi*(2d0*pi*cos(tw1*omega1/2d0) &
          +tw1*omega1*sin(2d0*pi*tt/tw1)*sin(omega1*tt)))
      end if
    end do
    T1_T2_tmp = T1_T2 + t1_start

  case('Esin2sin')

    stop "Esin2sin is not implemented"

  case('Asin2cos')

    ! pulse shape : A(t)=f0/omega*sin(Pi t/T)**2 *cos (omega t+phi_CEP*2d0*pi)
    ! pump laser
    do i=is, ie
      t=t0+i*delta_t
      tt = t
      if (tt<tw1) then
        Ac_ext_t(:,i) = -epdir_re1(:)*f0_1/omega1*(sin(pi*tt/tw1))**2*cos(omega1*tt+phi_CEP1*2d0*pi)
      end if
    end do

  case('input')

    call add_Ac_from_file(trim(file_input1), t0, delta_t, is, ie, Ac_ext_t)
    !  stop "ae_shape1='input' is not implemented"

  case('Asin2_cw')

    ! pulse shape : A(t)=f0/omega*sin(Pi t/T)**2 *cos (omega t+phi_CEP*2d0*pi)
    ! pump laser
    do i=is, ie
      t=t0+i*delta_t
      tt = t
      if (tt<tw1*0.5d0) then
        Ac_ext_t(:,i) = -Epdir_re1(:)*f0_1/omega1*(sin(pi*tt/tw1))**2*cos(omega1*tt+phi_CEP1*2d0*pi)
      else
        Ac_ext_t(:,i) = -Epdir_re1(:)*f0_1/omega1*cos(omega1*tt+phi_CEP1*2d0*pi)
      end if
    end do

  case('none')

    Ac_ext_t(:,:) = 0d0

  case default

    stop "Invalid pulse_shape_1 parameter!"

  end select

! Probe
  select case(ae_shape2)
  case('impulse')

    do i=is, ie
      t=t0+i*delta_t
      tt = t - 0.5d0*tw1 - T1_T2_tmp
      if(tt > 0d0)then
        Ac_ext_t(1,i) = Ac_ext_t(1,i) + epdir_re2(1)*e_impulse
        Ac_ext_t(2,i) = Ac_ext_t(2,i) + epdir_re2(2)*e_impulse
        Ac_ext_t(3,i) = Ac_ext_t(3,i) + epdir_re2(3)*e_impulse
      end if
    end do

  case('Acos2','Acos3','Acos4','Acos6','Acos8')

    select case(ae_shape2)
    case('Acos2'); npower = 2
    case('Acos3'); npower = 3
    case('Acos4'); npower = 4
    case('Acos6'); npower = 6
    case('Acos8'); npower = 8
    case default
      stop 'Error in init_Ac.f90'
    end select

    do i=is, ie
      t=t0+i*delta_t
      tt = t - 0.5d0*tw1 - T1_T2_tmp
      if (abs(tt)<0.5d0*tw2) then
        Ac_ext_t(:,i)=Ac_ext_t(:,i) &
          -f0_2/omega2*(cos(pi*tt/tw2))**npower &
          *aimag( (epdir_re2(:) + zI*epdir_im2(:)) &
          *exp(zI*(omega2*tt+phi_CEP2*2d0*pi))  &
          )
      end if
    end do

  case('Ecos2')

    if(phi_CEP2 /= 0.75d0)then
      stop "Error: phi_cep2 should be 0.75 when ae_shape2 is 'Ecos2'."
    end if
    if(sum(abs(epdir_im2(:)))>1.0d-8)then
      stop "Error: ae_shape2 should be 'Acos2' when epdir_im2 is used."
    end if

    do i=is, ie
      t=t0+i*delta_t
      tt = t - 0.5d0*tw1 - T1_T2_tmp
      if (abs(tt)<0.5d0*tw2) then
        Ac_ext_t(:,i)=Ac_ext_t(:,i) &
          -epdir_re2(:)*f0_2/(8d0*pi**2*omega2 - 2d0*tw2**2*omega2**3) &
          *( &
          (-4d0*pi**2+tw2**2*omega2**2 + tw2**2*omega2**2*cos(2d0*pi*tt/tw2))*cos(omega2*tt) &
          +2d0*pi*(2d0*pi*cos(tw2*omega2/2d0) &
          +tw2*omega2*sin(2d0*pi*tt/tw2)*sin(omega2*tt)))
      end if
    end do

  case('Esin2sin')

    stop "Esin2sin is not implemented"

  case('Asin2cos')

      ! pulse shape : A(t)=f0/omega*sin(Pi t/T)**2 *cos (omega t+phi_CEP*2d0*pi)
    ! probe laser

    do i=is, ie
      t=t0+i*delta_t
      tt = t
      if ( (tt-T1_T2_tmp>0d0) .and. (tt-T1_T2_tmp<tw2) ) then
        Ac_ext_t(:,i) = Ac_ext_t(:,i) &
          &-Epdir_re2(:)*f0_2/omega2*(sin(pi*(tt-T1_T2_tmp)/tw2))**2*cos(omega2*(tt-T1_T2_tmp)+phi_CEP2*2d0*pi)
      endif
    end do

  case('input')
  case('Asin2_cw')
  case('none')
  case default

    stop "Invalid pulse_shape_2 parameter!"

  end select

  return
End Subroutine calc_Ac_ext_t

!===================================================================================================================================

Subroutine calc_E_ext(ipulse,t,E_ext,amp_flag)
  use salmon_global,  only: ae_shape1,ae_shape2,tw1,tw2,omega1,omega2,&
                            phi_cep1,phi_cep2,E_amplitude1,E_amplitude2,t1_t2
  use math_constants, only: Pi
  implicit none
  integer,intent(in)      :: ipulse  ! 1: first pulse, 2: second pulse
  real(8),intent(in)      :: t
  real(8),intent(out)     :: E_ext
  character(1),intent(in) :: amp_flag
  real(8)                 :: alpha,beta,theta1,theta2,amp
  character(16)           :: tae_shape

  if(ipulse==1)then
    ! cos(theta1)**2
    theta1 = Pi/tw1*(t-0.5d0*tw1)
    alpha  = Pi/tw1
    ! cos(theta2)
    theta2 = omega1*(t-0.5d0*tw1)+phi_cep1*2d0*pi
    beta   = omega1
    !other
    tae_shape = ae_shape1
    amp       = E_amplitude1
  else if(ipulse==2)then
    ! cos(theta1)**2
    theta1 = Pi/tw2*(t-t1_t2-0.5d0*tw1)
    alpha  = Pi/tw2
    ! cos(theta2)
    theta2 = omega2*(t-t1_t2-0.5d0*tw1)+phi_cep2*2d0*pi
    beta   = omega2
    !other
    tae_shape = ae_shape2
    amp       = E_amplitude2
  end if

  if(tae_shape=='Ecos2')then
    E_ext = cos(theta1)**2*cos(theta2)
  else if(tae_shape=='Acos2')then
    E_ext = -( -alpha*sin(2d0*theta1)*cos(theta2)   &
               -beta*cos(theta1)**2*sin(theta2) )/beta
  end if

  if(amp_flag=='y') E_ext = amp*E_ext

  return
End Subroutine calc_E_ext

!===================================================================================================================================

subroutine add_Ac_from_file(filename, t0, delta_t, is, ie, Ac_ext_t)
  use read_rtdata_file, only: count_rows_from_rtdata_file, load_Ac_from_rtdata_file, load_Ac_from_rtdata_file_Ac_tot
  implicit none
  character(256), intent(in) :: filename
  real(8),intent(in) :: t0
  real(8),intent(in) :: delta_t
  integer,intent(in) :: is, ie
  real(8),intent(inout) :: Ac_ext_t(1:3, is:ie)

  integer :: n_dat, i
  real(8), allocatable :: t_dat(:), Ac_dat(:, :)
  real(8) :: Ac(1:3)   !,t,dt_dat

  n_dat = count_rows_from_rtdata_file(trim(filename))

  allocate(t_dat(1:n_dat), Ac_dat(1:3, 1:n_dat))
  !(choose Ac_ext or Ac_tot to load)
  call load_Ac_from_rtdata_file(trim(filename), n_dat, t_dat, Ac_dat)
!  call load_Ac_from_rtdata_file_Ac_tot(trim(filename), n_dat,t_dat,Ac_dat)

  do i = is, ie
      call interp_Ac(t0 + delta_t * i, Ac)
      Ac_ext_t(1:3, i) = Ac_ext_t(1:3, i) + Ac(1:3)
  end do

  deallocate(t_dat, Ac_dat)

  return

contains

  subroutine interp_Ac(tt, Ac)
    ! Linear interpolation of `Ac` at time `tt`.
    ! If a value outside the range of the data is requested,
    ! the value at the nearest end is returned.
    implicit none
    real(8), intent(in) :: tt
    real(8), intent(out) :: Ac(1:3)
    integer :: ii
    real(8) :: x, rii
    if (tt < t_dat(1)) then
      Ac(1:3) = 0d0 !Ac_dat(1:3, 1)
    else if (t_dat(n_dat) <= tt) then
      Ac(1:3) = 0d0 !Ac_dat(1:3, n_dat)
    else
      ! Linear estimation of row index: ii
      ii = int(floor((tt - t_dat(1)) / (t_dat(n_dat) - t_dat(1)) * (n_dat - 1))) + 1
      if (.not. (t_dat(ii) <= tt .and. tt < t_dat(ii+1))) then
        ! Bluteforce search of row index: ii
        do ii = 1, n_dat-1
          if ((t_dat(ii) <= tt) .and. (tt < t_dat(ii+1))) exit
        end do
      end if
      x = (tt - t_dat(ii)) / (t_dat(ii+1) - t_dat(ii))
      Ac(1:3) = (1d0 - x) * Ac_dat(1:3, ii) + x * Ac_dat(1:3, ii+1)
    end if
    return
  end subroutine interp_Ac

end subroutine add_Ac_from_file

!===================================================================================================================================

subroutine calcVbox(mg,lg,itt_t,system,rt,Vbox)
  use structures, only: s_rgrid, s_dft_system, s_rt, s_scalar
  use communication, only: comm_is_root
  use misc_routines, only: get_wtime
  use inputoutput
  implicit none
  type(s_rgrid),intent(in) :: mg,lg
  integer :: itt_t
  type(s_dft_system),intent(inout) :: system
  type(s_scalar),intent(inout)     :: Vbox
  type(s_rt),intent(in) :: rt
  integer :: ix,iy,iz,jj
  integer :: ix_sta_Vbox(3),ix_end_Vbox(3)
  integer :: ipulse
  real(8) :: env_trigon_1,env_trigon_2
  integer,parameter :: Nd = 4

  if(iperiodic==0)then
    if(yn_md=='y' .or. yn_out_rvf_rt=='y')then
      do jj=1,3
        if(lg%is(jj)==mg%is(jj))then
          ix_sta_Vbox(jj)=mg%is(jj)
        else
          ix_sta_Vbox(jj)=mg%is(jj)-Nd
        end if
        if(lg%ie(jj)==mg%ie(jj))then
          ix_end_Vbox(jj)=mg%ie(jj)
        else
          ix_end_Vbox(jj)=mg%ie(jj)+Nd
        end if
      end do
    else
      ix_sta_Vbox(1:3)=mg%is(1:3)
      ix_end_Vbox(1:3)=mg%ie(1:3)
    end if
  end if

  !$OMP parallel do collapse(2) private(ix,iy,iz)
  do iz=mg%is(3),mg%ie(3)
  do iy=mg%is(2),mg%ie(2)
  do ix=mg%is(1),mg%ie(1)
    Vbox%f(ix,iy,iz)=0.d0
  end do
  end do
  end do

  system%vec_E_ext(1:3)=0.d0

  if(iperiodic==0)then
    if(ae_shape1 == 'impulse')then
      continue
    else
        if(dt*dble(itt_t) <= tw1)then
          ipulse=1
          call calc_E_ext(ipulse,dt*dble(itt_t),env_trigon_1,'y')
        !$OMP parallel do collapse(2) private(ix,iy,iz)
          do iz=ix_sta_Vbox(3),ix_end_Vbox(3)
          do iy=ix_sta_Vbox(2),ix_end_Vbox(2)
          do ix=ix_sta_Vbox(1),ix_end_Vbox(1)
            Vbox%f(ix,iy,iz)=Vbox%f(ix,iy,iz)  &
                             +( epdir_re1(1)*lg%coordinate(ix,1)   &
                               +epdir_re1(2)*lg%coordinate(iy,2)   &
                               +epdir_re1(3)*lg%coordinate(iz,3) )*env_trigon_1  &
                             +( epdir_im1(1)*lg%coordinate(ix,1)   &
                               +epdir_im1(2)*lg%coordinate(iy,2)   &
                               +epdir_im1(3)*lg%coordinate(iz,3) )*env_trigon_1
          end do
          end do
          end do
          system%vec_E_ext(1:3)=system%vec_E_ext(1:3)+E_amplitude1*epdir_re1(1:3)*env_trigon_1   &
                                                       +E_amplitude1*epdir_im1(1:3)*env_trigon_1
        end if
        if(abs(dt*dble(itt_t)-0.5d0*tw1-t1_t2) < 0.5d0*tw2)then
          ipulse=2
          call calc_E_ext(ipulse,dt*dble(itt_t),env_trigon_2,'y')
          !$OMP parallel do collapse(2) private(ix,iy,iz)
          do iz=ix_sta_Vbox(3),ix_end_Vbox(3)
          do iy=ix_sta_Vbox(2),ix_end_Vbox(2)
          do ix=ix_sta_Vbox(1),ix_end_Vbox(1)
            Vbox%f(ix,iy,iz)=Vbox%f(ix,iy,iz)   &
                             +( epdir_re2(1)*lg%coordinate(ix,1)   &
                               +epdir_re2(2)*lg%coordinate(iy,2)   &
                               +epdir_re2(3)*lg%coordinate(iz,3) )*env_trigon_2  &
                             +( epdir_im2(1)*lg%coordinate(ix,1)   &
                               +epdir_im2(2)*lg%coordinate(iy,2)   &
                               +epdir_im2(3)*lg%coordinate(iz,3) )*env_trigon_2
          end do
          end do
          end do
          system%vec_E_ext(1:3)=system%vec_E_ext(1:3)+E_amplitude2*epdir_re2(1:3)*env_trigon_2   &
                                                       +E_amplitude2*epdir_im2(1:3)*env_trigon_2
        end if
    end if
  end if
  system%vec_E(1:3)=system%vec_E_ext(1:3)
  system%vec_Ac(1:3)=system%vec_Ac(1:3)-system%vec_E(1:3)*dt
  system%vec_Ac_ext(1:3)=system%vec_Ac(1:3)

  if(num_dipole_source>=1)then
    if(dt*dble(itt_t) <= tw1)then
      ipulse=1
      call calc_E_ext(ipulse,dt*dble(itt_t),env_trigon_1,'n')
!$OMP parallel do collapse(2) private(ix,iy,iz)
      do iz=mg%is(3),mg%ie(3)
      do iy=mg%is(2),mg%ie(2)
      do ix=mg%is(1),mg%ie(1)
        Vbox%f(ix,iy,iz)=Vbox%f(ix,iy,iz)+rt%vonf%f(ix,iy,iz)*env_trigon_1
      end do
      end do
      end do
    end if
    if(abs(dt*dble(itt_t)-0.5d0*tw1-t1_t2) < 0.5d0*tw2)then
      ipulse=2
      call calc_E_ext(ipulse,dt*dble(itt_t),env_trigon_2,'n')
!$OMP parallel do collapse(2) private(ix,iy,iz)
      do iz=mg%is(3),mg%ie(3)
      do iy=mg%is(2),mg%ie(2)
      do ix=mg%is(1),mg%ie(1)
        Vbox%f(ix,iy,iz)=Vbox%f(ix,iy,iz)+rt%vonf%f(ix,iy,iz)*env_trigon_2
      end do
      end do
      end do
    end if
  end if

  return
end subroutine calcVbox

!===================================================================================================================================

subroutine set_vonf(mg,lg,Hgs,rt)
!$ use omp_lib
  use structures, only: s_rgrid,s_rt
  use communication, only: comm_is_root, comm_summation
  use salmon_global
  implicit none
  real(8),intent(in) :: Hgs(3)
  type(s_rgrid) :: mg,lg
  type(s_rt) :: rt
  integer :: i
  integer :: ix,iy,iz,iix,iiy,iiz
  integer :: max_icell
  real(8) :: rr

  if(iperiodic==0)then
    max_icell=0
  else if(iperiodic==3)then
    max_icell=2
  end if

  rt%vonf%f=0.d0

  do i=1,num_dipole_source
    do iiz=-max_icell,max_icell
    do iiy=-max_icell,max_icell
    do iix=-max_icell,max_icell
      do iz=mg%is(3),mg%ie(3)
      do iy=mg%is(2),mg%ie(2)
      do ix=mg%is(1),mg%ie(1)
        rr=sqrt((lg%coordinate(ix,1)-(cood_dipole_source(1,i)+dble(iix*lg%num(1))*Hgs(1)))**2 &
               +(lg%coordinate(iy,2)-(cood_dipole_source(2,i)+dble(iiy*lg%num(2))*Hgs(2)))**2 &
               +(lg%coordinate(iz,3)-(cood_dipole_source(3,i)+dble(iiz*lg%num(3))*Hgs(3)))**2)
        if(rr>rad_dipole_source)then
          rt%vonf%f(ix,iy,iz)=rt%vonf%f(ix,iy,iz)  &
                        -(vec_dipole_source(1,i)*(lg%coordinate(ix,1)-(cood_dipole_source(1,i)+dble(iix*lg%num(1))*Hgs(1))) &
                         +vec_dipole_source(2,i)*(lg%coordinate(iy,2)-(cood_dipole_source(2,i)+dble(iiy*lg%num(2))*Hgs(2))) &
                         +vec_dipole_source(3,i)*(lg%coordinate(iz,3)-(cood_dipole_source(3,i)+dble(iiz*lg%num(3))*Hgs(3))))/rr**3
        else
          rt%vonf%f(ix,iy,iz)=rt%vonf%f(ix,iy,iz)  &
                        -(vec_dipole_source(1,i)*(lg%coordinate(ix,1)-(cood_dipole_source(1,i)+dble(iix*lg%num(1))*Hgs(1))) &
                         +vec_dipole_source(2,i)*(lg%coordinate(iy,2)-(cood_dipole_source(2,i)+dble(iiy*lg%num(2))*Hgs(2))) &
                         +vec_dipole_source(3,i)*(lg%coordinate(iz,3)-(cood_dipole_source(3,i)+dble(iiz*lg%num(3))*Hgs(3))))&
                                                                                /rad_dipole_source**3
        end if
      end do
      end do
      end do
    end do
    end do
    end do
  end do

  return
end subroutine set_vonf

end module em_field
!--------10--------20--------30--------40--------50--------60--------70--------80--------90--------100-------110-------120-------130


text:
! Sample solution to the problem posed in Lab 12.
! Version with coarse grain parallelism

program main

    use omp_lib
    use random_util, only: init_random_seed

    implicit none
    integer :: seed, nthreads
    integer, parameter :: ntests = 7
    real (kind=8), allocatable :: x(:)
    real (kind=8) :: sum, count1, count4
    real (kind=8), dimension(ntests)  :: mean,q1,q4
    integer :: i,n,k,kstart,kend,thread_num,points_per_thread

    ! Specify number of threads to use:
    nthreads = 1       ! need this value in serial mode
    !$ nthreads = 2
    !$ call omp_set_num_threads(nthreads)
    !$ print "('Using OpenMP with ',i3,' threads')", nthreads

    print *, "input seed"
    read *, seed

    call init_random_seed(seed)

    print *, "      n          mean        quartile 1    quartile 4"

    points_per_thread = (ntests + nthreads - 1) / nthreads
    print *, "points_per_thread = ",points_per_thread

    !$omp parallel private(k,x,count1,count4,sum,kstart,kend,thread_num,n)

    thread_num = 0     ! needed in serial mode
    !$ thread_num = omp_get_thread_num()    ! unique for each thread

    ! Determine start and end index for the set of points to be
    ! handled by this thread:
    kstart = thread_num * points_per_thread + 1
    kend = min((thread_num+1) * points_per_thread, ntests)

    !$omp critical
    print 201, thread_num, kstart, kend
    !$omp end critical
201 format("Thread ",i2," will take k = ",i6," through k = ",i6)

    do k=kstart,kend
        n = 10**(k+1)
        allocate(x(n))
        call random_number(x)
        sum = 0.d0
        count1 = 0
        count4 = 0
        do i=1,n
            sum = sum + x(i)
            if (x(i) < 0.25) then
                count1 = count1 + 1
                endif
            if (x(i) > 0.75) then
                count4 = count4 + 1
                endif
            enddo
        mean(k) = sum / float(n)
        q1(k) = count1 / float(n)
        q4(k) = count4 / float(n)
        deallocate(x)
        enddo
    !$omp end parallel

    do k=1,ntests
        n = 10**(k+1)
        print 601, n,mean(k),q1(k),q4(k)
 601    format(i10,3f15.8)
        enddo

end program main

max_stars_repo_name: philipwangdk/HPC

max_issues_repo_name: philipwangdk/HPC

max_forks_repo_name: philipwangdk/HPC


text:
great_ne(jj1,b1).
great_ne(m1,v1).
great_ne(p1,y1).
great_ne(dd1,ii1).
great_ne(q1,t1).
great_ne(t1,f1).
great_ne(hh1,a1).
great_ne(cc1,hh1).
great_ne(v1,ff1).
great_ne(dd1,m1).
great_ne(q1,ee1).
great_ne(u1,f1).
great_ne(q1,w1).
great_ne(v1,y1).
great_ne(w1,a1).
great_ne(q1,a1).
great_ne(z1,u1).

ext: f

max_stars_repo_name: manoelfranca/cilppp

max_issues_repo_name: manoelfranca/cilppp

max_forks_repo_name: manoelfranca/cilppp


text:
! one dimensional fourier transform, to use this subroutine
! 1.define {emin,de} and make them common in the main program
! 2.define an output array for {out}
	program main

	implicit real*8(a-h,o-z)
	integer*4 Nt,i,j
	real*8 dt,Emax,x,x0,p0,a0,Pi,t(nt)
	real*8 ymin,y,en,emin,de
	complex*16 :: data(Nt)
	complex*16,intent(OUT) :: out(Nt)
        complex*16 :: im

        open(100,file='fft.dat',status='new')

	open(10,file='corr')
        do i=1,nt
          read(10,*) t(i), data(i)
        enddo

	Emax = 0.1d0
	im = (0d0,1d0)
	Pi = 4.0d0*atan(1d0)

        de = 1d0/(Nt*dt)
        emin = -0.18d0
        out = (0d0,0d0)

	i = 1
	do while (Emin+(i-1)*de < Emax)

        en = Emin + (i-1)*de
!	if(en > Emax) goto 100
        do j=1,Nt
          t = 0d0 +j*dt
!	data(j) = exp(-a0*(x-x0)**2+im*p0*(x-x0))
          out(i) = data(j)*exp(im*en*t)*dt + out(i)
        enddo

        i = i+1

        write(100,1000) y,out(i)
        end do

1000    format(1000(e14.7,1x))

        return
        end program


text:
*     ---------------------------------------------
*     enter double precision
*     ---------------------------------------------
      SUBROUTINE ENTDOU (QUEST,X)
      IMPLICIT NONE

*     FORMAL_PARAMETERS:
      CHARACTER QUEST*(*)
      DOUBLE PRECISION X

*     integer codes for tokens
      INCLUDE 'rdndec.gin'

*     value info
      INCLUDE 'rddecinf.inc'

**    local variables
      CHARACTER INPUT*40
      INTEGER TTYPE, IB, IE, NERR, NWAR
      SAVE

*     ask the question and read the answer
10    CONTINUE
      CALL RDERRI ('Keyboard',0,0,.TRUE.,.FALSE.)
      CALL ENTCHA (QUEST,INPUT)
      CALL PARSWORD (INPUT,TTYPE,IB,IE,NERR,NWAR)

      IF (NERR.EQ.0) THEN
         IF (TTYPE .EQ. TFLOAT) THEN
*           floating point value
            X = VFLOAT
         ELSE
*           other data type
            CALL RDERR (2,'Not a floating point value')
            GOTO 10
         END IF
      ELSE
         GOTO 10
      END IF

      RETURN
      END

ext: for

max_issues_repo_name: zhumiansanren/FIRI-WHS


text:
Module shaderclass

    use types, only : vector
    use image, only : RGBAimage

    implicit none

    !abstract shader type
    type, abstract :: shader
        real :: varying_intensity(3)
        Contains
            procedure(generic_frag), deferred, pass :: fragment
            procedure(generic_vert), deferred, pass :: vertex
    end type shader

    abstract interface
        function generic_vert(this, vertex, light)
            use triangleclass
            use camera, only : viewport, projection, modelview, v2m, m2v
            import :: shader
            class(shader) :: this
            type(triangle), intent(IN) :: vertex
            type(vector),   intent(IN) :: light
            real :: generic_vert(4,3)

        end function generic_vert
    end interface

    abstract interface
        logical function generic_frag(this, bar_c, colour)
            use image, only : RGBA, operator(*)
            use types, only: operator(.dot.), vector
            import :: shader
            class(shader) :: this
            type(vector), intent(IN)    :: bar_c
            type(RGBA),   intent(INOUT) :: colour
        end function generic_frag
    end interface

    !default shader which is a gourand shader
    type, extends(shader) :: gourand
        Contains
            procedure, pass(this) :: fragment => fragment_gourand
            procedure, pass(this) :: vertex => vertex_gourand
    end type gourand
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    !shader which is a texture mapped shader
    type, extends(shader) :: tmap
        type(vector)    :: varying_uv(3)
        type(RGBAimage) :: texture
        Contains
            procedure, pass(this) :: fragment => fragment_tmap
            procedure, pass(this) :: vertex => vertex_tmap
    end type tmap

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    !shader which is a texture mapped shader
    type, extends(shader) :: wireframe
        real :: line_thickness = 0.0025
        type(vector) :: altitudes
        integer :: width, height
        Contains
            procedure, pass(this) :: fragment => fragment_wire
            procedure, pass(this) :: vertex => vertex_wire
    end type wireframe
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    private :: fragment_gourand, vertex_gourand, fragment_tmap, vertex_tmap, fragment_wire, vertex_wire
    public :: shader, gourand, tmap, wireframe

    Contains
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    logical function fragment_gourand(this, bar_c, colour)

        use image, only : RGBA, operator(*)
        use types, only: operator(.dot.)

        implicit none

        class(gourand) :: this
        type(vector), intent(IN)    :: bar_c
        type(RGBA),   intent(INOUT) :: colour

        type(vector) :: tmp
        real :: intensity

        tmp = vector(this%varying_intensity(1), this%varying_intensity(2), this%varying_intensity(3))

        intensity = abs(tmp .dot. bar_c)
        colour = RGBA(255,255,255,255) * intensity
        fragment_gourand = .false.

    end function fragment_gourand

    function vertex_gourand(this, vertex, light)

        use triangleclass
        use camera, only : viewport, projection, modelview, v2m, m2v

        implicit none

        class(gourand)             :: this
        type(triangle), intent(IN) :: vertex
        type(vector),   intent(IN) :: light
        real :: gl_vertex(4,1), vertex_gourand(4,3)
        integer :: j

        do j = 1, size(vertex%vert)
            this%varying_intensity(j) = max(0.01, (vertex%norms(j) .dot. light))
            gl_vertex = v2m(vertex%vert(j))
            gl_vertex = matmul(matmul(matmul(viewport,projection),modelview),gl_vertex)
            vertex_gourand(:, j) = gl_vertex(:,1)
        end do
    end function vertex_gourand

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    logical function fragment_tmap(this, bar_c, colour)

        use image, only : RGBA, operator(*), get_pixel
        use types, only: operator(.dot.), operator(*), operator(+)

        implicit none

        class(tmap) :: this
        type(vector), intent(IN)    :: bar_c
        type(RGBA),   intent(INOUT) :: colour

        type(vector) :: tmp, uv
        real :: intensity

        tmp = vector(this%varying_intensity(1), this%varying_intensity(2), this%varying_intensity(3))

        intensity = tmp .dot. bar_c
        uv = this%varying_uv(1)*bar_c%x + this%varying_uv(2)*bar_c%y + this%varying_uv(3)*bar_c%z
        call get_pixel(this%texture, int(uv%x*this%texture%width), int(uv%y*this%texture%height), colour)
        colour = colour * intensity
        fragment_tmap = .false.

    end function fragment_tmap

    function vertex_tmap(this, vertex, light)

        use triangleclass
        use camera, only : viewport, projection, modelview, v2m, m2v

        implicit none

        class(tmap)              :: this
        type(triangle), intent(IN) :: vertex
        type(vector),   intent(IN) :: light

        real    :: gl_vertex(4,1), vertex_tmap(4,3)
        integer :: j

        do j = 1, size(vertex%vert)
            this%varying_uv(j) = vertex%uvs(j)
            this%varying_intensity(j) = max(0., (vertex%norms(j) .dot. light))
            gl_vertex = v2m(vertex%vert(j))
            gl_vertex = matmul(matmul(matmul(viewport,projection),modelview),gl_vertex)
            vertex_tmap(:, j) = gl_vertex(:,1)
        end do
    end function vertex_tmap

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    logical function fragment_wire(this, bar_c, colour)

        use image, only : RGBA, operator(*), get_pixel, alpha_comp
        use types, only: operator(.dot.), operator(*), operator(+)

        implicit none

        class(wireframe) :: this
        type(vector), intent(IN)    :: bar_c
        type(RGBA),   intent(INOUT) :: colour

        type(vector) :: tmp
        real         :: eintensity, d, intensity

        tmp = bar_c%x * vector(this%altitudes%x,0.,0.) + bar_c%y * vector(0.,this%altitudes%y,0.) + &
              bar_c%z * vector(0.,0.,this%altitudes%z)

        d = min(tmp%x,tmp%y,tmp%z)
        ! if(d > 2.)then
        !     tmp = vector(this%varying_intensity(1), this%varying_intensity(2), this%varying_intensity(3))
        !     intensity = abs(tmp .dot. bar_c)
        !     colour = RGBA(255,0,0,255) * intensity
        ! else
            eintensity = 2.**(-.1*d*d)
            colour =  RGBA(255,255,255,255) * eintensity + RGBA(255,0,0,255) * (1. - eintensity)
        ! end if

        fragment_wire = .false.

    end function fragment_wire

    function vertex_wire(this, vertex, light)

        use triangleclass
        use camera, only : viewport, projection, modelview, v2m, m2v

        implicit none

        class(wireframe)           :: this
        type(triangle), intent(IN) :: vertex
        type(vector),   intent(IN) :: light

        type(vector) :: p0,p1,p2,v0,v1,v2

        real    :: gl_vertex(4,1), vertex_wire(4,3), area
        integer :: j

        do j = 1, size(vertex%vert)
            this%varying_intensity(j) = max(0., (vertex%norms(j) .dot. light))
            gl_vertex = v2m(vertex%vert(j))
            gl_vertex = matmul(matmul(matmul(viewport,projection),modelview),gl_vertex)
            vertex_wire(:,j) = gl_vertex(:,1)
        end do

        p0%z = 0.
        p1%z = 0.
        p2%z = 0.

        p0%x = this%width *  vertex%vert(1)%x/vertex_wire(4,1)
        p0%y = this%height * vertex%vert(1)%y/vertex_wire(4,1)

        p1%x = this%width *  vertex%vert(2)%x/vertex_wire(4,2)
        p1%y = this%height * vertex%vert(2)%y/vertex_wire(4,2)

        p2%x = this%width *  vertex%vert(3)%x/vertex_wire(4,3)
        p2%y = this%height * vertex%vert(3)%y/vertex_wire(4,3)

        v0 = p2 - p1
        v1 = p2 - p0
        v2 = p1 - p0

        area = abs(v1%x*v2%y - v1%y*v2%x)

        this%altitudes = vector(area/magnitude(v0), area/magnitude(v1), area/magnitude(v2))

        ! this%altitudes = vector(altitude(vertex%vert(1),vertex%vert(2),vertex%vert(3)), &
        !                         altitude(vertex%vert(2),vertex%vert(3),vertex%vert(1)), &
        !                         altitude(vertex%vert(3),vertex%vert(1),vertex%vert(2)))
        ! print*,this%altitudes
        ! stop
    end function vertex_wire

    real function altitude(a, b, c)

        use types

        implicit none

        type(vector), intent(IN) :: a, b, c
        type(vector) :: ba, bc, ba_on_bc

        ba = a - b
        bc = c - b
        ba_on_bc = (ba .dot. bc) * bc
        altitude = magnitude(ba - ba_on_bc)

    end function altitude
end Module shaderclass

max_stars_repo_name: lewisfish/FortGL

max_issues_repo_name: lewisfish/FortGL

max_forks_repo_name: lewisfish/FortGL


text:
!--------------------------------------------------------------------------------
! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
! This file is part of FLEUR and available as free software under the conditions
! of the MIT license as expressed in the LICENSE file in more detail.
!--------------------------------------------------------------------------------

PROGRAM inpgen
!----------------------------------------------------------------------------+
!   Set up a FLEUR inp-file from basic input data; for use and docu please   !
!   refer to inpgen.html (or see http://www.flapw.de/docs/inpgen.html)       !
!                                                                            |
!   The program is based on the input conventions of the FLAIR-code, so that !
!   some compatibility is ensured. The symmetry generator was written by     !
!   M.Weinert and implemented in the FLAIR-code by G.Schneider.              !
!                                                                    gb`02   |
!----------------------------------------------------------------------------+
  USE m_juDFT
  USE m_inpgen_help
  USE m_read_inpgen_input
  USE m_make_crystal
  USE m_make_atomic_defaults
  USE m_make_defaults
  USE m_make_kpoints
  USE m_winpxml
  USE m_xsf_io
  USE m_types_input
  USE m_types_atoms
  USE m_types_cell
  USE m_types_sym
  USE m_types_noco
  USE m_types_vacuum
  USE m_types_banddos
  USE m_types_hybinp
  USE m_types_xcpot_inbuild_nofunction
  USE m_types_forcetheo
  USE m_types_kpts
  USE m_types_gfinp
  USE m_types_enpara
  USE m_types_oneD
  USE m_types_sliceplot
  USE m_types_stars
  use m_read_old_inp
  use m_fleurinput_read_xml
  USE m_types_mpinp
  USE m_constants

      IMPLICIT NONE

      REAL,    ALLOCATABLE :: atompos(:, :),atomid(:)
      CHARACTER(len=20), ALLOCATABLE :: atomLabel(:)
      LOGICAL               :: l_fullinput,l_explicit,l_inpxml,l_include(4)

      TYPE(t_input)    :: input
      TYPE(t_atoms)    :: atoms
      TYPE(t_cell)     :: cell
      TYPE(t_sym)      :: sym
      TYPE(t_noco)     :: noco
      TYPE(t_vacuum)   :: vacuum
      TYPE(t_banddos)  :: banddos
      TYPE(t_mpinp)    :: mpinp
      TYPE(t_hybinp)   :: hybinp
      TYPE(t_xcpot_inbuild_nf)::xcpot
      TYPE(t_enpara)   :: enpara
      TYPE(t_forcetheo):: forcetheo
      TYPE(t_kpts)     :: kpts
      TYPE(t_oned)     :: oned
      TYPE(t_sliceplot):: sliceplot
      TYPE(t_stars)    :: stars
      TYPE(t_gfinp)    :: gfinp
      TYPE(t_enparaXML):: enparaxml

      CHARACTER(len=40):: kpts_str,filename
      LOGICAL          :: l_exist
      INTEGER          :: idum

      INTERFACE
       FUNCTION dropDefaultEConfig() BIND(C, name="dropDefaultEconfig")
         USE iso_c_binding
         INTEGER(c_int) dropDefaultEConfig
       END FUNCTION dropDefaultEConfig
      END INTERFACE

      CALL judft_init(oUnit,.FALSE.)

      kpts_str=""

      !Start program and greet user
      CALL inpgen_help()
      l_explicit=judft_was_argument("-explicit")

      INQUIRE(file='default.econfig',exist=l_exist)
      IF (.NOT.l_exist) idum=dropDefaultEconfig()

      OPEN(oUnit,file='out')

      INQUIRE(file='inp.xml',exist=l_inpxml)
      IF (l_inpxml.AND..NOT.(judft_was_argument("-inp.xml").or.judft_was_argument("-overwrite")))&
           CALL judft_error("inp.xml exists and can not be overwritten")

      IF (judft_was_argument("-inp")) THEN
         call read_old_inp(input,atoms,cell,stars,sym,noco,vacuum,forcetheo,&
              sliceplot,banddos,enpara,xcpot,kpts,hybinp, oneD)
         l_fullinput=.TRUE.
      ELSEIF (judft_was_argument("-inp.xml")) THEN
         !not yet
         l_fullinput=.true. !will be set to false if old inp.xml is read
         call Fleurinput_read_xml(cell,sym,atoms,input,noco,vacuum,&
         sliceplot=Sliceplot,banddos=Banddos,hybinp=hybinp,oned=Oned,xcpot=Xcpot,kpts=Kpts,enparaXML=enparaXML,old_version=l_fullinput)
         Call Cell%Init(Dot_product(Atoms%Volmts(:),Atoms%Neq(:)))
         call atoms%init(cell)
         Call Sym%Init(Cell,Input%Film)
         CALL xcpot%init(atoms%ntype)
         CALL enpara%init_enpara(atoms,input%jspins,input%film,enparaXML)
         l_fullinput=.TRUE.
      ELSEIF(judft_was_argument("-f")) THEN
         !read the input

         CALL read_inpgen_input(atompos,atomid,atomlabel,kpts_str,&
              input,sym,noco,vacuum,stars,xcpot,cell,hybinp)
         l_fullinput=.FALSE.
      ELSE
         CALL judft_error("You should either specify -inp,-inp.xml or -f command line options. Check -h if unsure")
      ENDIF
      IF (.NOT.l_fullinput) THEN
         !First we determine the spacegoup and map the atoms to groups
         CALL make_crystal(input%film,atomid,atompos,atomlabel,vacuum%dvac,noco,&
              cell,sym,atoms)

         !All atom related parameters are set here. Note that some parameters might
         !have been set in the read_input call before by adding defaults to the atompar module
         CALL make_atomic_defaults(input,vacuum,cell,oneD,atoms,enpara)

         !Set all defaults that have not been specified before or can not be specified in inpgen
         CALL make_defaults(atoms,sym,cell,vacuum,input,stars,&
                   xcpot,noco,mpinp,hybinp)
      ENDIF
      !
      ! k-points can also be modified here
      !
      call make_kpoints(kpts,cell,sym,hybinp,input%film,noco%l_ss.or.noco%l_soc,&
                        input%bz_integration,kpts_str)
      !
      !Now the IO-section
      !
      IF (.NOT.l_inpxml.or.judft_was_argument("-overwrite")) THEN
         call determine_includes(l_include)
         !the inp.xml file
         !CALL dump_FleurInputSchema()
         filename="inp.xml"
         if (judft_was_argument("-o")) filename=juDFT_string_for_argument("-o")
         CALL w_inpxml(&
              atoms,vacuum,input,stars,sliceplot,forcetheo,banddos,&
              cell,sym,xcpot,noco,oneD,mpinp,hybinp,kpts,enpara,gfinp,&
              l_explicit,l_include,filename)
         if (.not.l_include(1)) CALL sym%print_XML(99,"sym.xml")
      ENDIF
      IF (.NOT.l_include(2).OR.judft_was_argument("-k")) THEN
        !inquire(file="kpts.xml",exist=l_include(2))
        !if (l_include(2)) call system("rm kpts.xml")
        CALL kpts%print_XML(99,"kpts.xml")
      endif
      ! Structure in  xsf-format
      OPEN (55,file="struct.xsf")
      CALL xsf_WRITE_atoms(55,atoms,input%film,.FALSE.,cell%amat)
      CLOSE (55)
      CLOSE(6)

      CALL juDFT_end("All done")

    CONTAINS
      SUBROUTINE determine_includes(l_include)
        LOGICAL,INTENT(out)::l_include(4)  !kpts,operations,species,position

        CHARACTER(len=100)::str=''
        LOGICAL           ::incl

        l_include=[.FALSE.,.FALSE.,.TRUE.,.TRUE.]

        IF (judft_was_argument("-inc")) THEN
           str=judft_string_for_argument("-inc")

           DO WHILE(LEN_TRIM(str)>0)
              IF (str(1:1)=='-') THEN
                 incl=.FALSE.
                 str=str(2:)
              ELSE
                 incl=.TRUE.
                 IF (str(1:1)=='+') str=str(2:)
              ENDIF
              SELECT CASE(str(1:1))
              CASE ('k','K')
                 l_include(1)=incl
              CASE ('o','O')
                 l_include(2)=incl
              CASE ('s','S')
                 l_include(3)=incl
              CASE ('p','P')
                 l_include(4)=incl
              CASE ('a','A')
                 l_include(:)=incl
              END SELECT
              IF (INDEX(str,"'")>0) THEN
                 str=str(INDEX(str,"'")+1:)
              ELSE
                 str=""
              END IF
           END DO
        ENDIF

      END SUBROUTINE determine_includes

    END PROGRAM inpgen

max_stars_repo_name: MRedies/FLEUR

max_issues_repo_name: MRedies/FLEUR

max_forks_repo_name: MRedies/FLEUR


text:
C Copyright (c) 2008-2017 National Technology & Engineering Solutions
C of Sandia, LLC (NTESS).  Under the terms of Contract DE-NA0003525 with
C NTESS, the U.S. Government retains certain rights in this software.
C
C Redistribution and use in source and binary forms, with or without
C modification, are permitted provided that the following conditions are
C met:
C
C     * Redistributions of source code must retain the above copyright
C       notice, this list of conditions and the following disclaimer.
C
C     * Redistributions in binary form must reproduce the above
C       copyright notice, this list of conditions and the following
C       disclaimer in the documentation and/or other materials provided
C       with the distribution.
C
C     * Neither the name of NTESS nor the names of its
C       contributors may be used to endorse or promote products derived
C       from this software without specific prior written permission.
C
C THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
C "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
C LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
C A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
C OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
C SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
C LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
C DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
C THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
C (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
C OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
C

C=======================================================================
      SUBROUTINE MUNESS (NUMESS, ISTAT, LESSEL, LESSDL,
     &  IDESS, NEESS, NEDSS, IXEESS, IXEDSS,
     &  LTEESS, LTSSS, FACSS,
     &  LTEX, LTSX, TDX, IXESS, IXDSS, NEX, NDX, ISCR, USESDF)
C=======================================================================
C $Id: muness.f,v 1.2 2001/06/26 17:38:54 gdsjaar Exp $
C
C   --*** MUNESS *** (GJOIN) Compress and rearrange element side sets
C   --   Written by Amy Gilkey - revised 02/25/88
C   --
C   --MUNESS processes the element side sets according to the set status.
C   --Sets may be combined or deleted.
C   --
C   --Parameters:
C   --   NUMESS - IN/OUT - the number of element side sets
C   --   ISTAT - IN - the status of each set:
C   --      0 = same
C   --      - = delete
C   --      n = combine with set n
C   --   LESSEL - IN/OUT - the length of the element side sets element list
C   --   LESSDL - IN/OUT - the length of the element side sets dist-fact list
C   --   IDESS - IN/OUT - the element side set ID for each set
C   --   NEESS - IN/OUT - the number of elements for each set
C   --   NEDSS - IN/OUT - the number of dist-fact for each set
C   --   IXEESS - IN/OUT - the index of the first element for each set
C   --   IXEDSS - IN/OUT - the index of the first dist-fact for each set
C   --   LTEESS - IN/OUT - the elements for all sets
C   --   LTSSS - IN/OUT - the sides for all sets
C   --   FACSS - IN/OUT - the dist-fact for all sets
C   --   LTEX - SCRATCH - sized to hold the set elements
C   --   LTSX - SCRATCH - sized to hold the set sides
C   --   TDX - SCRATCH - sized to hold the set dist-fact
C   --   IXESS - SCRATCH - size = NUMESS
C   --   IXDSS - SCRATCH - size = NUMESS -- dist-fact
C   --   NEX - SCRATCH - size = NUMESS
C   --   NDX - SCRATCH - size = NUMESS -- dist-face
C   --   ISCR - SCRATCH - size = NUMESS

      INTEGER ISTAT(*)
      INTEGER IDESS(*)
      INTEGER NEESS(*), NEDSS(*)
      INTEGER IXEESS(*)
      INTEGER LTEESS(*), LTEX(*)
      INTEGER LTSSS(*), LTSX(*)
      INTEGER IXESS(*), IXEDSS(*)
      INTEGER NEX(*), NDX(*)
      INTEGER ISCR(*)
      REAL    FACSS(*), TDX(*)
      LOGICAL USESDF

      IF (NUMESS .LE. 0) RETURN

      JESS = 0
      JNE = 1
      JND = 1
      DO 110 IESS = 1, NUMESS

         IF (ISTAT(IESS) .EQ. 0) THEN
            NINSET = 1
            ISCR(NINSET) = IESS
         ELSE IF (ISTAT(IESS) .EQ. IESS) THEN
            CALL GETALL (IESS, NUMESS, ISTAT, NINSET, ISCR)
         ELSE
            NINSET = 0
         END IF

         IF (NINSET .GT. 0) THEN
            JESS = JESS + 1
            IXESS(JESS) = IESS
            NEX(JESS) = 0
            NDX(JESS) = 0
         END IF
         DO 100 ISET = 1, NINSET
            N = ISCR(ISET)
            CALL MOVINT (NEESS(N), LTEESS(IXEESS(N)), LTEX(JNE))
            CALL MOVINT (NEESS(N), LTSSS(IXEESS(N)),  LTSX(JNE))
            IF (USESDF) THEN
              CALL MOVREA (NEDSS(N), FACSS(IXEDSS(N)),  TDX(JND))
            ENDIF
            JNE = JNE + NEESS(N)
            JND = JND + NEDSS(N)
            NEX(JESS) = NEX(JESS) + NEESS(N)
            NDX(JESS) = NDX(JESS) + NEDSS(N)
  100    CONTINUE
  110 CONTINUE

      CALL ORDIX (JESS, IXESS, NUMESS, IDESS, ISCR, IDESS)
      CALL MOVINT (JESS, NEX, NEESS)
      CALL MOVINT (JESS, NDX, NEDSS)
      NUMESS = JESS
      JNE = 1
      JND = 1
      DO 120 IESS = 1, NUMESS
         IXEESS(IESS) = JNE
         IXEDSS(IESS) = JND
         JNE = JNE + NEESS(IESS)
         JND = JND + NEDSS(IESS)
  120 CONTINUE
      LESSEL = JNE - 1
      LESSDL = JND - 1
      CALL MOVINT (LESSEL, LTEX, LTEESS)
      CALL MOVINT (LESSEL, LTSX, LTSSS)
      if (usesdf) then
        CALL MOVREA (LESSDL, TDX,  FACSS)
      end if
      RETURN
      END

max_stars_repo_name: mathstuf/seacas

max_stars_repo_licenses: BSD-3-Clause NetCDF Zlib MIT

max_issues_repo_name: mathstuf/seacas

max_issues_repo_licenses: BSD-3-Clause NetCDF Zlib MIT

max_forks_repo_name: mathstuf/seacas

max_forks_repo_licenses: BSD-3-Clause NetCDF Zlib MIT


text:
      SUBROUTINE MOVAVR (ITASK,NAME,IP,IN,OUT)
      IMPLICIT NONE

*     FORMAL_PARAMETERS:
      INTEGER ITASK,IP
      REAL IN, OUT
      CHARACTER*(*) NAME

**    local variables
      INTEGER MNP,MNNAM,NNAM
      PARAMETER (MNP=100,MNNAM=10)
      REAL ARR(MNP,MNNAM)
      INTEGER NPA(MNNAM)
      CHARACTER*31 NAMESA(MNNAM),LNAME
      INTEGER I1,N,IFINDC,NP
      SAVE

      IF (ITASK.EQ.1) THEN

*        initialize storage arrays

         DO 10 I1=1,MNNAM
            NPA(I1)    = 0
            NAMESA(I1) = ' '
10       CONTINUE
         NNAM = 0
      ELSE

*        calculate running average

*        look up NAME in list of names
         IF (IP.GT.MNP) CALL FATALERR ('MOVAVR',
     &      'maximum number of points exceeded')
         LNAME = NAME
         N = IFINDC (NAMESA,MNNAM,1,NNAM,LNAME)
         IF (N.EQ.0) THEN
*           NAME was not found in list of names, add to list
            NNAM = NNAM+1
            IF (NNAM.GT.MNNAM) CALL FATALERR ('MOVAVR','too many names')
            N = NNAM
            NAMESA(N) = LNAME
         END IF

*        shift old values for NAME
         NP = MIN (NPA(N)+1, IP)

         DO 20 I1=NP,2,-1
            ARR(I1,N) = ARR(I1-1,N)
20       CONTINUE

*        add new value to array
         ARR(1,N) = IN

*        calculate running average by summing and dividing
         OUT = 0.
         DO 30 I1=1,NP
            OUT = OUT+ARR(I1,N)
30       CONTINUE
         OUT    = OUT/REAL (NP)
         NPA(N) = NP
      END IF

      RETURN
      END


text:
      subroutine ksando (obj,g,fscale,offset,ncon,nobj)
      implicit double precision (a-h,o-z)
      dimension obj(*),g(*),fscale(*),offset(*)
      data eps/1.0e-6/
c
c          routine to compute scale and offset for obj
c
c          author   - Gregory A. Wrenn
c          location - Lockheed Engineering and Sciences Co.
c                     144 Research Drive
c                     Hampton, Va. 23666
c
c          last modification - 17 July 1996
c
      if (ncon .le. 0) go to 20
      gmax = g(1)
      if (ncon .eq. 1) go to 20
      do 10 i = 2,ncon
        if (g(i) .gt. gmax) gmax = g(i)
   10 continue
c
   20 continue
      do 30 j = 1,nobj
        fscl = abs(obj(j))
        foff = abs(obj(j))
        if (fscl .lt. 1.0) fscl = 1.0
        if (foff .gt. 1.0) foff = 1.0
        fscale(j) = fscl
        offset(j) = 0.0
c
c          compute offset if ncon > 0
c
        if (ncon .le. 0) go to 30
        if (obj(j) .ge. 0.0) offset(j) = -foff - gmax - eps
        if (obj(j) .lt. 0.0) offset(j) =  foff - gmax + eps
   30 continue
c
      return
      end

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


text:
! { dg-lto-do link }
! { dg-lto-options {{ -O1 -flto }} }
! { dg-suppress-ld-options "-O1" }

      SUBROUTINE ylm4(ylm)
      COMPLEX, INTENT (OUT):: ylm(1)
      INTEGER l,m
      COMPLEX ylms
      REAL, ALLOCATABLE, SAVE  :: ynorm(:)
      ylms = 0
      DO m = 1, 1
         DO l = m, 1
            ylm(m) = conjg(ylms)*ynorm(m)
         ENDDO
      ENDDO
      END SUBROUTINE ylm4

      PROGRAM test
      END

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
!> NodelForceInLine
!!  @param iShiftU = nint(commonPar(1))
!!	@param iShiftDeltaU  = nint(CommonPar(2))
!!	@param iFemType = iFemType
!!	@param loadType = nint(CommonPar(4))
!!	@param LoadParam = CommonPar(5:10)
!! @author Rocha, Felipe Figueredo

!     ------------------------------------------------------------------
Subroutine NodalForceInLine(AE, BE, MaxLRows, XLL, NDim, iDofT, NodELT, Sol0, &
						Sol1, CommonPar, Param, JParam, DelT, DTm, Time)
!     ------------------------------------------------------------------
	use funcAux
	use loadingLib
	use finiteStrainLib
	use ptsGaussLib2

	IMPLICIT NONE

	!   ===== SUBROUTINE ARGUMENTS  =======
	integer :: MaxLRows,Ndim,iDofT,NodELT ! all integers
	Real*8 :: DelT, DTm,Time ! all reals
	! Reals Vectors and matrices
	Real*8  :: AE(MaxLRows,MaxLRows), BE(MaxLRows), XLL(*), Sol0(*), Sol1(*), &
		Param(*), JParam(*), CommonPar(*)

	!   =====   END ARGUMENTS  =======

	Integer :: nG, i,j,k,l,m,n,p,q, A, B, ApRow, BpCol , App,Bqp , ip, ipp ! counters
	Integer :: NodG,  pOrder , NGP, iShiftU, iShiftDeltaU, iBu, iSimplex, iFemType
	integer, parameter :: Ndim2D = 2
	Real*8 :: PhiA, dv , pm,  dPhiA(Ndim2D)
	real*8 , allocatable ::  SolU(:) , Xel(:) !! all have dimension NodG*NdimE
	Real*8 :: normal0(Ndim2D) , vecU(Ndim2D), lengthLine, lengthLine0
	type(ptGaussClass) :: PtG

	iShiftU = nint(CommonPar(1))
	iShiftDeltaU = nint(CommonPar(2))
	iFemType = nint(CommonPar(3))
	call chooseLoad(nint(CommonPar(4)))
	call pressureLoad(pm,Time,DelT,CommonPar(5:10))

	call setNodG(iFEMtype,NodG)

	call getSliceAllocate(SolU,Sol1,1,NodG,iShiftU + 1 ,iShiftU + Ndim2D, iDofT)
	call getSliceAllocate(Xel,XLL,1,NodG,1 ,Ndim2D, Ndim)

	VecU = Xel(Ndim2D + 1:2*Ndim2D) - Xel(1:Ndim2D)
	lengthLine0 = dsqrt(dot_product(vecU,vecU))
	vecU = vecU/lengthline0

	AE = 0.0d0
	BE = 0.0d0

	if(NodG == 2) then
!~ 		write(0,*) "===========================>>>>>>>>>> pm = " , pm
		BE(iShiftDeltaU + 1) = 0.5d0*pm*lengthLine0*VecU(1) !! in x
		BE(iShiftDeltaU + 2) = 0.5d0*pm*lengthLine0*VecU(2) !! in y
		BE(iDofT + iShiftDeltaU + 1) = 0.5d0*pm*lengthLine0*VecU(1) !! in x
		BE(iDofT + iShiftDeltaU + 2) = 0.5d0*pm*lengthLine0*VecU(2) !! in y
	end if

	if(NodG == 3) then
!~ 		write(0,*) "===========================>>>>>>>>>> pm = " , pm
		BE(iShiftDeltaU + 1) = 0.16666666667d0*pm*lengthLine0*VecU(1) !! in x
		BE(iShiftDeltaU + 2) = 0.16666666667d0*pm*lengthLine0*VecU(2) !! in y
		BE(iDofT + iShiftDeltaU + 1) = 0.66666666667d0*pm*lengthLine0*VecU(1) !! in x
		BE(iDofT + iShiftDeltaU + 2) = 0.66666666667d0*pm*lengthLine0*VecU(2) !! in y
		BE(2*iDofT + iShiftDeltaU + 1) = 0.16666666667d0*pm*lengthLine0*VecU(1) !! in x
		BE(2*iDofT + iShiftDeltaU + 2) = 0.16666666667d0*pm*lengthLine0*VecU(2) !! in y
	end if

end subroutine

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Subroutine NodalForceInLineS(Coupling,CommonPar,iDofT,Ndim,MaxLRows,iAdd)
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	use funcAux
	use ptsGaussLib2

	IMPLICIT NONE
	Integer :: MaxLRows,Ndim,iDofT, iAdd
	Integer   Coupling(MaxLRows,MaxLRows)
	Real*8 CommonPar(*)
	integer :: NodG, Ndim2D, iFemType
	Integer :: i,j,ipRow,ipCol,k,l, iShiftDeltaU

	iShiftDeltaU=nint(CommonPar(2))
	iFemType = nint(commonPar(3))
    call setNodG(iFemtype, NodG)

	Ndim2D = 2

	Coupling = 0
	do i = 1, NodG
		ipRow = (i-1)*iDofT + iShiftDeltaU
		do j = 1, NodG
			ipCol = (j-1)*iDofT + iShiftDeltaU
			do k = 1,Ndim2D
				do l = 1,Ndim2D
					Coupling(ipRow + k,ipCol + l) =  1
				end do
			end do
		end do
	end do

end Subroutine

max_stars_repo_name: felipefr/Piola

max_issues_repo_name: felipefr/Piola

max_forks_repo_name: felipefr/Piola


text:
program test_vegflx
  use simsphere_mod
  implicit none

! Call subroutine VEGFLX to determine variables changed, develop tests
! Will need at least 2 tests: XLEG >= 0 and XLEG < 0

  real :: EVAPV, rprime, XLEGN

  EVAPV = 0.3
  TAF = 10.0
  CHF = 10.0
  F = 100.0
  CHG = 0.1
  QAF = 0.5
  HG = 1.0
  CHG = 100.5
  XLEGN = 1.0
  XLEF = 1.2
  CHA = 100.0
  TA = 0.4
  CHF = 10.0
  TF = 1.0
  RST = 1.0
  QA = 100.0
  QSTF = 10.0

  call VEGFLX(EVAPV)

  write(*,*) QAF

end program test_vegflx

max_stars_repo_name: derekmorr/simsphere

max_issues_repo_name: derekmorr/simsphere

max_forks_repo_name: derekmorr/simsphere


subroutine trace_bend_2d(xzt,yzt,xst,yst,ips,	tout)

real xrmin(1000),yrmin(1000)
real t_segm(1000)
integer indexx(100)
real dxtmp(1000),dytmp(1000),xtmp(1000),ytmp(1000)

common/ray_param/ds_ini,ds_segm_min,bend_min0,bend_max0

common/ray/ nodes,xray(1000),yray(1000)
common/ray_part/ npart,xpart(1000),ypart(1000),spart(1000),kod_part(1000)
common/shift/ dxpart(1000),dypart(1000)

totdist=sqrt((xzt-xst)*(xzt-xst)+(yzt-yst)*(yzt-yst))
tout=0

if(totdist.lt.ds_segm_min) then
		x1=xzt
		x2=xst
		xm=(x1+x2)/2.
		y1=yzt
		y2=yst
		ym=(y1+y2)/2.

		vvv=velocity (xm,ym,ips)
		tout=ds/vvv
		return
end if
nsegm_max = int_best(totdist / ds_segm_min)

if(nsegm_max.eq.0) then
	return
end if

call streight_line(xzt,yzt,xst,yst,ips, tout)
!write(*,*)' streight line: tout=',tout
!write(*,*)ds_ini,ds_segm_min,bend_min0,bend_max0
!pause

tmin=tout

do nsegm=1,nsegm_max

!open(22,file='nodes.dat')

	indexx=0
	dpart=1./nsegm
	val_cur0=bend_max0/nsegm
	xtmp=xray
	ytmp=yray

	do iseg=1,nsegm

		part1=(iseg-1)*dpart
		part2=iseg*dpart

		call part_ray(xzt,yzt,xst,yst,part1,part2)
		!write(*,*)' iseg=',iseg,' part1=',part1,' part2=',part2

!write(22,*) xpart(1),ypart(1)
		ttt=0
		do inode=1,npart-1
			x1=xpart(inode)
			x2=xpart(inode+1)
			xm=(x1+x2)/2.
			y1=ypart(inode)
			y2=ypart(inode+1)
			ym=(y1+y2)/2.
			ds=sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1))
			vvv=velocity (xm,ym,ips)
			ttt=ttt+ds/vvv
			!write(*,*)inode,ypart(inode),dypart(inode),' v=',vvv
		end do
		tini=ttt
		!write(*,*)' tini=',tini,' npart=',npart

		dxpart=0
		dypart=0
		dxtmp=0
		dytmp=0
		valtot=0
		tmin=tini
		t_segm(iseg)=tmin
		val_cur=val_cur0

332		continue
		do icase=1,2
			ind=0
			val_bend = (-1)**icase * val_cur

331			continue
			val = val_bend

			sss2=spart(npart)/2.
			sss=0
			do inode=2,npart-1

				x1=xpart(inode-1)
				y1=ypart(inode-1)
				x2=xpart(inode)
				y2=ypart(inode)
				x3=xpart(inode+1)
				y3=ypart(inode+1)

				ds=sqrt((x2-x1)**2+(y2-y1)**2)
				sss=sss+ds

				dx=x3-x1
				dy=y3-y1
				ds2=sqrt(dx*dx+dy*dy)

				d_value = val * (1-abs(sss-sss2)/sss2)

				dxtmp(inode)= - d_value * dy / ds2
				dytmp(inode)=   d_value * dx / ds2
				!write(*,*)inode,d_value,dxpart(inode),dypart(inode)

			end do

			ttt=0
			sss=0
			do inode=1,npart-1
				x1=	xpart(inode)+	dxpart(inode)+		dxtmp(inode)
				x2=	xpart(inode+1)+	dxpart(inode+1)+	dxtmp(inode+1)
				xm=(x1+x2)/2.

				y1=	ypart(inode)+	dypart(inode)+		dytmp(inode)
				y2=	ypart(inode+1)+	dypart(inode+1)+	dytmp(inode+1)
				ym=(y1+y2)/2.

				ds=sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1))

				sss=sss+ds
				vvv=velocity (xm,ym,ips)
				ttt=ttt+ds/vvv
				!write(*,*)inode,ypart(inode),dypart(inode),' v=',vvv
			end do

			tout=ttt
			!write(*,*)' val_bend=',val_bend,' tout=',tout

			if(tout.lt.tmin) then
				ind=ind+1
				tmin=tout
				dxpart=dxpart+dxtmp
				dypart=dypart+dytmp
				valtot=valtot+val_bend
				!write(*,*)' valtot=',valtot,' val_bend=',val_bend,' tmin=',tmin
				goto 331
			else
				if(ind.ne.0) exit
			end if

		end do

		val_cur = val_cur / 2.
		if(abs(val_cur).gt.bend_min0) goto 332

		do ipp=1,npart
			kod=kod_part(ipp)
			if(kod.eq.0) cycle
			xtmp(kod)=xtmp(kod)+dxpart(ipp)
			ytmp(kod)=ytmp(kod)+dypart(ipp)
			xpart(ipp)=xpart(ipp)+dxpart(ipp)
			ypart(ipp)=ypart(ipp)+dypart(ipp)
		end do
		t_segm(iseg)=tmin

	end do

333	continue
	ttot=0
	do iseg=1,nsegm
		!write(*,*)' indexx(ip)=',indexx(ip)
		ttot=ttot+t_segm(iseg)
	end do

!write(22,*)xray(nodes),yray(nodes)
!close(22)

	xray=xtmp
	yray=ytmp
	!write(*,*)' val_cur0=',val_cur0,' ttot=',ttot

	nrefl=0
	call remeshing(nrefl)

	!open(21,file='ray_iter.bln')
	!write(21,*)nodes
	!do i=1,nodes
		!write(21,*)xray(i),yray(i)
	!end do
	!close(21)
	!pause

end do
tout=ttot

ttt=0
sss=0
do inode=1,nodes-1
	x1=	xray(inode)
	x2=	xray(inode+1)
	xm=(x1+x2)/2.

	y1=	yray(inode)
	y2=	yray(inode+1)
	ym=(y1+y2)/2.

	ds=sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1))

	sss=sss+ds

	vvv=velocity (xm,ym,ips)
	!dv= dv_mod_2d(xm,ym)

!	write(*,*)xm,ym,' vvv=',vvv
	!call vel_mod_2d(xm,ym, dv000,dv060,dv120)
	!write(*,*)' dv000=',dv000,' dv060=',dv060,' dv120=',dv120

	ttt=ttt+ds/vvv
	!write(*,*)inode,ypart(inode),dypart(inode),' v=',vvv
end do

tout=ttt

return
end


text:
! Program to test the stack variable size limit.
program stack
   call sub1
   call sub2 (1)
contains

   ! Local variables larger than 32768 in byte size shall be placed in static
   ! storage area, while others be put on stack by default.
   subroutine sub1
      real a, b(32768/4), c(32768/4+1)
      integer m, n(1024,4), k(1024,1024)
      a = 10.0
      b = 20.0
      c = 30.0
      m = 10
      n = 20
      k = 30
      if ((a .ne. 10.0).or.(b(1) .ne. 20.0).or.(c(1) .ne. 30.0)) call abort
      if ((m .ne. 10).or.(n(256,4) .ne. 20).or.(k(1,1024) .ne. 30)) call abort
   end subroutine

   ! Local variables defined in recursive subroutine are always put on stack.
   recursive subroutine sub2 (n)
      real a (32769)
      a (1) = 42
      if (n .ge. 1) call sub2 (n-1)
      if (a(1) .ne. 42) call abort
      a (1) = 0
   end subroutine
end

max_issues_repo_name: vidkidz/crossbridge

max_forks_repo_name: vidkidz/crossbridge


text:
These are reviews of The Graduate from 2005 and 2006.  For more current reviews, visit The Graduate the entry.

20050721 11:22:33 nbsp I personally try to avoid this location.  Last time I went there they were giving away free chewing tobacco at the door.  Users/AllanRae

20050905 14:13:00 nbsp I think the bouncers at this place are training to be cops.  You better not be in a pool game when they decide to close the doors.  Users/PatrickLynch

20060105 18:59:10 nbsp The Graduate was having 1418 Night on Graduation Night last quarter, and as such no graduates could attend The Graduate on their Graduation Night.  I think they should have their name revoked! Users/KrisFricke

20060205 23:43:08 nbsp If you live in a residential area nearby, you can (and should) complain to the Davis Police Department police when the Grad gets too noisy, as they often do. Users/GrahamFreeman

20060206 12:59:14 nbsp The food here is pretty decent, I ate here a few times and enjoyed myself.  They always have a bunch of sports games on in the day time which is also nice.  The one thing I dont like about the day time is the amount of little kids running around.  Lots of soccer or baseball teams seem to come for postgame meals and such so it can get a little crazy. Users/BenWebber

20060313 19:36:31 nbsp Swing nights are coming back, after being absent for a couple years. This time its better than ever. DJ Scottie is a personal friend of mine, and hes great. All of you should try swing dancing! Lessons are usually at 9 or 9:30, on certain Sundays. Just look around for flyers or ask me. Thursday nights at the Grad are Salsa nights. Thats also lots of fun. The bouncers, staff, etc. arent all that mean like people say...just dont do anything stupid or illegal. Users/LaurenWong

20060807 19:01:05 nbsp Couple of changes to the information about the Grad: Country nights are still Saturday and Monday, with an occasional Sunday thrown in, and it costs $11 if youre under 21, $8 for under 21 with a student ID, $8 if youre over 21, and $5 if youre over 21 and have a student ID.  If youre over 21 and you get there before 9:00pm, theyll let you in for free.  Dance lessons start at 9:30pm and general dancing starts at 10:00pm.  It usually goes until about 1:30am, with a break for hiphop and slow dances around midnight. Users/SteffiChristiana

20060912 20:49:10 nbsp Tuesday nights are now Rock It! which is going to be a weekly Indie Rock dance party. 18+, drink specials and fun music Users/JordanSmart

20060923 20:27:38 nbsp The pizza isnt very good, and its even worse reheated the next day.  Users/AnnaJones

20060928 17:36:14 nbsp Dude, all ages dancing was halarious! Also, I think its a law that all 1st and 2nd year students go to the Grad once you get back into Davis city limits. Users/SteveBlais

20060929 12:41:15 nbsp I like their fries, but their chicken tenders are WAY to spicey... Users/JoeRunnels

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


! ------------------------------------------------------------------
! Programmer(s): David J. Gardner, Cody J. Balos @ LLNL
!                Daniel R. Reynolds @ SMU
! ------------------------------------------------------------------
! SUNDIALS Copyright Start
! Copyright (c) 2002-2021, Lawrence Livermore National Security
! and Southern Methodist University.
! All rights reserved.
!
! See the top-level LICENSE and NOTICE files for details.
!
! SPDX-License-Identifier: BSD-3-Clause
! SUNDIALS Copyright End
! ------------------------------------------------------------------
! This test simulates a brusselator problem from chemical kinetics.
! This is an ODE system with 3 components, Y = [u,v,w], satisfying
! the equations,
!
!    du/dt = a - (w+1)*u + v*u^2
!    dv/dt = w*u - v*u^2
!    dw/dt = (b-w)/ep - w*u
!
! for t in the interval [0.0, 10.0], with initial conditions
! Y0 = [3.9, 1.1, 2.8], and parameter values a=1.2, b=2.5, and
! ep=1.0e-5
!
! Here, all three components exhibit a rapid transient change during
! the first 0.2 time units, followed by a slow and smooth evolution.
! ------------------------------------------------------------------

module ode_mod

  !======= Inclusions ===========
  use, intrinsic :: iso_c_binding

  !======= Declarations =========
  implicit none

  ! number of equations
  integer(c_long), parameter :: neq = 3

  ! ODE parameters
  double precision, parameter :: a  = 1.2d0
  double precision, parameter :: b  = 2.5d0
  double precision, parameter :: ep = 1.0d-5

contains

  ! ----------------------------------------------------------------
  ! RhsFn provides the right hand side function for the
  ! ODE: dy/dt = f(t,y)
  !
  ! Return values:
  !    0 = success,
  !    1 = recoverable error,
  !   -1 = non-recoverable error
  ! ----------------------------------------------------------------
  integer(c_int) function RhsFn(tn, sunvec_y, sunvec_f, user_data) &
       result(ierr) bind(C,name='RhsFn')

    !======= Inclusions ===========
    use, intrinsic :: iso_c_binding
    use fsundials_nvector_mod

    !======= Declarations =========
    implicit none

    ! calling variables
    real(c_double), value :: tn        ! current time
    type(N_Vector)        :: sunvec_y  ! solution N_Vector
    type(N_Vector)        :: sunvec_f  ! rhs N_Vector
    type(c_ptr),    value :: user_data ! user-defined data

    ! pointers to data in SUNDIALS vectors
    real(c_double), pointer :: yvec(:)
    real(c_double), pointer :: fvec(:)

    !======= Internals ============

    ! get data arrays from SUNDIALS vectors
    yvec => FN_VGetArrayPointer(sunvec_y)
    fvec => FN_VGetArrayPointer(sunvec_f)

    ! fill RHS vector
    fvec(1) = a  -  (yvec(3) + 1.0d0) * yvec(1)  +  yvec(2) * yvec(1) * yvec(1)
    fvec(2) = yvec(3) * yvec(1)  -  yvec(2) * yvec(1) * yvec(1)
    fvec(3) = (b-yvec(3))/ep - yvec(3) * yvec(1)

    ! return success
    ierr = 0
    return

  end function RhsFn

  ! ----------------------------------------------------------------
  ! JacFn: The Jacobian of the ODE hand side function J = df/dy
  !
  ! Return values:
  !    0 = success,
  !    1 = recoverable error,
  !   -1 = non-recoverable error
  ! ----------------------------------------------------------------
  integer(c_int) function JacFn(tn, sunvec_y, sunvec_f, sunmat_J, &
       user_data, tmp1, tmp2, tmp3) &
       result(ierr) bind(C,name='JacFn')

    !======= Inclusions ===========
    use, intrinsic :: iso_c_binding
    use fsundials_nvector_mod
    use fsunmatrix_dense_mod
    use fsundials_matrix_mod

    !======= Declarations =========
    implicit none

    ! calling variables
    real(c_double), value :: tn               ! current time
    type(N_Vector)        :: sunvec_y         ! current solution N_Vector
    type(N_Vector)        :: sunvec_f         ! current rhs N_Vector
    type(SUNMatrix)       :: sunmat_J         ! Jacobian SUNMatrix
    type(c_ptr), value    :: user_data        ! user-defined data
    type(N_Vector)        :: tmp1, tmp2, tmp3 ! workspace N_Vectors

    ! pointers to data in SUNDIALS vectors
    real(c_double), pointer :: yvec(:)

    ! pointer to data in SUNDIALS matrix
    real(c_double), pointer :: Jmat(:)

    !======= Internals ============

    ! get data array from SUNDIALS vector
    yvec => FN_VGetArrayPointer(sunvec_y)

    ! get data arrays from SUNDIALS vectors
    Jmat => FSUNDenseMatrix_Data(sunmat_J)

    ! fill Jacobian matrix
    Jmat = [-(yvec(3)+1.0d0) + 2.0d0*yvec(1)*yvec(2),&
            yvec(3) - 2.0d0*yvec(1)*yvec(2), -yvec(3),&
            yvec(1)*yvec(1), -yvec(1)*yvec(1), 0.0d0,&
            -yvec(1), yvec(1), -1.0d0/ep - yvec(1)]

    ! return success
    ierr = 0
    return

  end function JacFn

end module ode_mod

program main

  !======= Inclusions ===========
  use, intrinsic :: iso_c_binding

  use fsundials_context_mod
  use fcvode_mod                 ! Fortran interface to CVODE
  use fnvector_serial_mod        ! Fortran interface to serial N_Vector
  use fsunmatrix_dense_mod       ! Fortran interface to dense SUNMatrix
  use fsunlinsol_dense_mod       ! Fortran interface to dense SUNLinearSolver
  use fsundials_linearsolver_mod ! Fortran interface to generic SUNLinearSolver
  use fsundials_matrix_mod       ! Fortran interface to generic SUNMatrix
  use fsundials_nvector_mod      ! Fortran interface to generic N_Vector
  use ode_mod                    ! ODE functions

  !======= Declarations =========
  implicit none

                                                 ! local variables
  real(c_double)                 :: tstart       ! initial time
  real(c_double)                 :: tend         ! final time
  real(c_double)                 :: rtol, atol   ! relative and absolute tolerance
  real(c_double)                 :: dtout        ! output time interval
  real(c_double)                 :: tout         ! output time
  real(c_double)                 :: tcur(1)      ! current time

  integer(c_int)                 :: ierr         ! error flag from C functions
  integer(c_int)                 :: nout         ! number of outputs

  integer                        :: outstep      ! output loop counter

  type(c_ptr)                    :: ctx          ! SUNDIALS context
  type(c_ptr)                    :: cvode_mem    ! CVODE memory
  type(N_Vector),        pointer :: sunvec_y     ! sundials vector
  type(SUNMatrix),       pointer :: sunmat_A     ! sundials matrix
  type(SUNLinearSolver), pointer :: sunlinsol_LS ! sundials linear solver

  ! solution vector, neq is set in the ode_mod module
  real(c_double) :: yvec(neq)

  !======= Internals ============

  ! initialize ODE
  tstart = 0.0d0
  tend   = 10.0d0
  tcur   = tstart
  tout   = tstart
  dtout  = (tend-tstart)/10.d0
  nout   = ceiling(tend/dtout)

  ! initialize solution vector
  yvec(1) = 3.9d0
  yvec(2) = 1.1d0
  yvec(3) = 2.8d0

  ! create SUNDIALS context
  ierr = FSUNContext_Create(c_null_ptr, ctx)

  ! create SUNDIALS N_Vector
  sunvec_y => FN_VMake_Serial(neq, yvec, ctx)
  if (.not. associated(sunvec_y)) then
     print *, 'ERROR: sunvec = NULL'
     stop 1
  end if

  ! create a dense matrix
  sunmat_A => FSUNDenseMatrix(neq, neq, ctx)
  if (.not. associated(sunmat_A)) then
     print *, 'ERROR: sunmat = NULL'
     stop 1
  end if

  ! create a dense linear solver
  sunlinsol_LS => FSUNLinSol_Dense(sunvec_y, sunmat_A, ctx)
  if (.not. associated(sunlinsol_LS)) then
     print *, 'ERROR: sunlinsol = NULL'
     stop 1
  end if

  ! create CVode memory
  cvode_mem = FCVodeCreate(CV_BDF, ctx)
  if (.not. c_associated(cvode_mem)) then
     print *, 'ERROR: cvode_mem = NULL'
     stop 1
  end if

  ! initialize CVode
  ierr = FCVodeInit(cvode_mem, c_funloc(RhsFn), tstart, sunvec_y)
  if (ierr /= 0) then
     print *, 'Error in FCVodeInit, ierr = ', ierr, '; halting'
     stop 1
  end if

  ! set relative and absolute tolerances
  rtol = 1.0d-5
  atol = 1.0d-10

  ierr = FCVodeSStolerances(cvode_mem, rtol, atol)
  if (ierr /= 0) then
     print *, 'Error in FCVodeSStolerances, ierr = ', ierr, '; halting'
     stop 1
  end if

  ! attach linear solver
  ierr = FCVodeSetLinearSolver(cvode_mem, sunlinsol_LS, sunmat_A);
  if (ierr /= 0) then
     print *, 'Error in FCVodeSetLinearSolver, ierr = ', ierr, '; halting'
     stop 1
  end if

  ! set Jacobian routine
  ierr = FCVodeSetJacFn(cvode_mem, c_funloc(JacFn))
  if (ierr /= 0) then
     print *, 'Error in FCVodeSetJacFn, ierr = ', ierr, '; halting'
     stop 1
  end if

  ! start time stepping
  print *, '   '
  print *, 'Finished initialization, starting time steps'
  print *, '   '
  print *, '      t           u           v           w'
  print *, '----------------------------------------------------'
  print '(1x,4(es12.5,1x))', tcur, yvec(1), yvec(2), yvec(3)
  do outstep = 1,nout

     ! call CVode
     tout = min(tout + dtout, tend)
     ierr = FCVode(cvode_mem, tout, sunvec_y, tcur, CV_NORMAL)
     if (ierr /= 0) then
        print *, 'Error in FCVODE, ierr = ', ierr, '; halting'
        stop 1
     endif

     ! output current solution
     print '(1x,4(es12.5,1x))', tcur, yvec(1), yvec(2), yvec(3)

  enddo

  ! diagnostics output
  call CVodeStats(cvode_mem)

  ! clean up
  call FCVodeFree(cvode_mem)
  ierr = FSUNLinSolFree(sunlinsol_LS)
  call FSUNMatDestroy(sunmat_A)
  call FN_VDestroy(sunvec_y)
  ierr = FSUNContext_Free(ctx)

end program Main

! ----------------------------------------------------------------
! CVodeStats
!
! Print CVODE statstics to standard out
! ----------------------------------------------------------------
subroutine CVodeStats(cvode_mem)

  !======= Inclusions ===========
  use iso_c_binding
  use fcvode_mod

  !======= Declarations =========
  implicit none

  type(c_ptr), intent(in) :: cvode_mem ! solver memory structure

  integer(c_int)  :: ierr          ! error flag

  integer(c_long) :: nsteps(1)     ! num steps
  integer(c_long) :: nfevals(1)    ! num function evals
  integer(c_long) :: nlinsetups(1) ! num linear solver setups
  integer(c_long) :: netfails(1)   ! num error test fails

  integer(c_int)  :: qlast(1)      ! method order in last step
  integer(c_int)  :: qcur(1)       ! method order for next step

  real(c_double)  :: hinused(1)    ! initial step size
  real(c_double)  :: hlast(1)      ! last step size
  real(c_double)  :: hcur(1)       ! step size for next step
  real(c_double)  :: tcur(1)       ! internal time reached

  integer(c_long) :: nniters(1)    ! nonlinear solver iterations
  integer(c_long) :: nncfails(1)   ! nonlinear solver fails

  integer(c_long) :: njevals(1)    ! num Jacobian evaluations

  !======= Internals ============

  ! general solver statistics
  ierr = FCVodeGetIntegratorStats(cvode_mem, nsteps, nfevals, nlinsetups, &
       netfails, qlast, qcur, hinused, hlast, hcur, tcur)
  if (ierr /= 0) then
     print *, 'Error in FCVodeGetIntegratorStats, ierr = ', ierr, '; halting'
     stop 1
  end if

  ! nonlinear solver statistics
  ierr = FCVodeGetNonlinSolvStats(cvode_mem, nniters, nncfails)
  if (ierr /= 0) then
     print *, 'Error in FCVodeGetNonlinSolvStats, ierr = ', ierr, '; halting'
     stop 1
  end if

  ! nonlinear solver statistics
  ierr = FCVodeGetNumJacEvals(cvode_mem, njevals)
  if (ierr /= 0) then
     print *, 'Error in FCVodeGetNumJacEvals, ierr = ', ierr, '; halting'
     stop 1
  end if

  print *, ' '
  print *, ' General Solver Stats:'
  print '(4x,A,i9)'    ,'Total internal steps taken =',nsteps
  print '(4x,A,i9)'    ,'Total rhs function calls   =',nfevals
  print '(4x,A,i9)'    ,'Num lin solver setup calls =',nlinsetups
  print '(4x,A,i9)'    ,'Num error test failures    =',netfails
  print '(4x,A,i9)'    ,'Last method order          =',qlast
  print '(4x,A,i9)'    ,'Next method order          =',qcur
  print '(4x,A,es12.5)','First internal step size   =',hinused
  print '(4x,A,es12.5)','Last internal step size    =',hlast
  print '(4x,A,es12.5)','Next internal step size    =',hcur
  print '(4x,A,es12.5)','Current internal time      =',tcur
  print '(4x,A,i9)'    ,'Num nonlinear solver iters =',nniters
  print '(4x,A,i9)'    ,'Num nonlinear solver fails =',nncfails
  print '(4x,A,i9)'    ,'Num Jacobian evaluations   =',njevals
  print *, ' '

  return

end subroutine CVodeStats


text:
C     LAST UPDATE 16/03/89
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
      SUBROUTINE IMDISP (ITERM,IPRINT,BUFF,NPIX,NRAST)
      IMPLICIT   NONE
C
C Purpose: Check for image display.
C
      REAL    BUFF(1)
      INTEGER ITERM,IPRINT,NPIX,NRAST
C
C ITERM  : TERMINAL INPUT STREAM
C IPRINT : TERMINAL OUTPUT STREAM
C BUFF   : IMAGE DATA
C NPIX   : NOS. OF PIXELS
C NRAST  : NOS. OF RASTERS
C
C Calls   2: ASKNO  , PUTIMAGE
C Called by:
C
C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-
C LOCAL VARIABLES:
C
      INTEGER IRC
      LOGICAL ASKNO,NOPLOT
C
C IRC    : RETURN CODE
C PLOT   : PLOT REQUIRED FLAG
C
C-----------------------------------------------------------------------
      WRITE (IPRINT,'(A,$)') ' Do you want to display image [Y/N] [N]: '
      NOPLOT=ASKNO (ITERM,IRC)
      IF (IRC.EQ.0) THEN
         IF (NOPLOT) THEN
            RETURN
         ELSE
            CALL PUTIMAGE (BUFF,NPIX,NRAST)
         ENDIF
      ENDIF
      END