fc2k_check_bc.f90

Go to the documentation of this file.

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! routine to check boundary condition consistency
! in ETS workflow
! the following checks are done:
! 1- nion>0 (number of ions in composition should be positive)
! 2 -ne bc is valid (>0)
! 3 - if electron density is predictive, one of ions should be interpretative (to force quasineutrality)
! 4 - ion bc are valid (number of bc should be the same as number of ions in composition)
!
!
! author:D.Yadykin
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 


subroutine interface_check_bc(nion_in,array_in,array_out,eflag_out,emessage_out)

implicit none

integer, intent(in),dimension(:) :: array_in
integer, intent(in) :: nion_in
integer, intent(out),dimension(:),pointer :: array_out
integer, intent(out) :: eflag_out
character(len=:), pointer, intent(out) :: emessage_out
!internal
character(len=:), allocatable :: emessage
!character(len=1000) :: emessage
!character(len=:), allocatable,target :: emessage_target 
character(len=1000) ::emessage_fix
integer :: eflag
integer :: nmes,narr_out
 integer,allocatable,dimension(:) ::local_array_out

interface
  subroutine check_bc (intarr_in,nion,intarr_out,error_flag,error_message)
   integer, intent(in), dimension(:) :: intarr_in
   integer, intent(in) :: nion
   integer, intent(out), allocatable,dimension(:) ::intarr_out
   integer, intent(out) :: error_flag
   character(len=:), allocatable :: error_message
   !character(len=1000) :: error_message
 end subroutine check_bc
end interface

eflag=0

call check_bc(array_in,nion_in,local_array_out,eflag,emessage)



if (eflag.lt.0) then
   nmes=len_trim(emessage)
   write(*,*) nmes,emessage
   allocate(character(len=nmes) :: emessage_out)
   emessage_out=emessage(1:nmes)
   eflag_out=eflag
   write(*,*) eflag_out,emessage_out
   !deallocate(emessage)
   !allocate(array_out(1))
   !array_out(1)=-1
end if
if (allocated(local_array_out)) then
  narr_out=size(local_array_out)
  allocate(array_out(narr_out))
  array_out(1:narr_out)=local_array_out(1:narr_out)
  !deallocate (local_array_out)
else
  allocate(array_out(1))
  array_out(1)=-1
  if (eflag_out.ge.0) then
    eflag_out=-1
    allocate(character(len=30) ::emessage_out)
    emessage_out='no output array allocated'
  end if
end if

return
end subroutine interface_check_bc


subroutine check_bc (intarr_in,nion,intarr_out,error_flag,error_message)

implicit none

integer, intent(in), dimension(:) :: intarr_in
integer, intent(out), allocatable, dimension(:) :: intarr_out
integer, allocatable, dimension(:) :: intarr
integer, allocatable, dimension (:) :: ni_bt,ti_bt,vi_bt

integer:: n_in,ti_shift,vi_shift,ni_shift,ne_shift,ne_bt,nbc_max
integer :: iion,ncount,ncount2,ncount3,ncount4
logical :: mask1,mask2,mask3,mask4

integer, intent(in)  :: nion
!this is set presently
integer, parameter :: nion_max=7

!error handling
integer, intent(out) :: error_flag
character(len=:), intent(out), allocatable :: error_message
!character(len=1000) :: error_message

!check 1
if (nion.lt.1) then
  error_flag=-1
  error_message='no ions found'
  return
end if

!this should be always, ne,te,ve
ne_shift=3
!this is set for nion_max ions 
ti_shift=ne_shift+1
ni_shift=ne_shift+nion_max+1
vi_shift=ne_shift+2*nion_max+1

!write(*,*) 'ti_shift,ni_shift,vi_shift',ti_shift,ni_shift,vi_shift

!for 7 ions this should be 24
n_in=size(intarr_in)


allocate(intarr(n_in))

!set otput
allocate(intarr_out(n_in))
intarr_out=intarr_in

!set work array
intarr=intarr_in

!set otput
intarr_out=intarr_in

!write(*,*) 'intarr', intarr
!get ne bc type
ne_bt=intarr(ne_shift)
!write(*,*) 'ne_bt',ne_bt

!allocate and get ni bc type array
allocate(ni_bt(nion_max))
ni_bt=intarr(ni_shift:ni_shift+nion_max-1)
!write(*,*) 'ni_bt',ni_bt

!check if ne is reliable (check 2)
if (ne_bt.lt.0) then
  error_flag=-1
  error_message='boundary type for ne less then zero, error'
  return
end if
!check if ne is interpretative (check 3)
!if not then
!check how many ni are predictive
!and if more than nion raise exception
if ((ne_bt.ge.1).and.(ne_bt.le.6)) then
  ncount=0
  do iion=1,nion_max
    if((ni_bt(iion).ge.1).and.(ni_bt(iion).le.6)) ncount=ncount+1
  enddo
  !write(*,*) 'ncount',ncount
  if (ncount.gt.nion-1) then
     error_flag=-1
     error_message='not possible to do predictive run, leave density of one of the ions or electrons free (OFF)'
     return
  end if   

end if

!get max number of bc set for ions (among t,v,n)
allocate(ti_bt(nion_max),vi_bt(nion_max))
ti_bt=intarr(ti_shift:ti_shift+nion_max-1)
vi_bt=intarr(vi_shift:vi_shift+nion_max-1)

!mask2=ni_bt.gt.0
!mask3=ti_bt.gt.0
!mask4=vi_bt.gt.0
ncount2=0
ncount3=0
ncount4=0
do iion=1,nion_max
  if (ni_bt(iion).ge.0) ncount2=ncount2+1
  if (ti_bt(iion).ge.0) ncount3=ncount3+1
  if (vi_bt(iion).ge.0) ncount4=ncount4+1
end do
deallocate(intarr)
deallocate(ni_bt,ti_bt,vi_bt)

!ncount2=all(ni_bt.gt.0)
!ncount3=all(ti_bt.gt.0)
!ncount4=all(ti_bt.gt.0)

!write(*,*) 'ncount2,ncount3,ncount4', ncount2,ncount3,ncount4

nbc_max=max(ncount2,ncount3,ncount4)

!write(*,*) 'nbc_max',nbc_max
!check 4
 if (nbc_max.lt.nion) then
   error_flag=-1
   error_message='error in composition, number of boundary conditions set is less than number of ions'
   return
 end if

return
end subroutine check_bc