subroutine BoundaryTreat() use particles use scattering use simulation integer :: nAlBin,deficiency, deficiencyLeft, nonUseParNum, i, j, addInx, ebin, valley, InjectComplete, TempNumCarriers integer, dimension(:), allocatable :: nonUseList real ( kind = 8 ):: k,kx,ky,kz,e,tc,binFreq,fai,ct,st !! local real ( kind = 8 ):: sq_wvk_sum,abs_wvk,FactorEne, FactorVel, movingDist,LowAlBound,UpAlBound real ( kind = 8 ), dimension(DIM) :: Temp_wvk, Temp_vel ! carrier position and wave vector k nonUseParNum = 0 TempNumCarriers = 0 do i = 1, num_maxCarriers ! count if (Par_use(i).eq.1) then TempNumCarriers = TempNumCarriers + 1 endif enddo num_Carriers = TempNumCarriers allocate (nonUseList(num_maxCarriers-num_Carriers)) do i = 1, num_maxCarriers if (Par_use(i).ne.1) then nonUseParNum = nonUseParNum + 1 nonUseList(nonUseParNum) = i endif enddo Deficiency_bin(:) = 0 deficiency = 0 deficiencyLeft = 0 if (Algo_bcCheck.eq.1) then ! left concentration check do i=1,num_bcCheckBin Deficiency_bin(i) = num_InitBin - bin_ParNumber(i) if ((i.eq.1).and.(PropertyNode.eq.1)) then Deficiency_bin(i) = Deficiency_bin(i)*0.5 endif deficiency = deficiency + Deficiency_bin(i) enddo if (deficiency.lt.0) then deficiency = 0 endif deficiencyLeft = deficiency if (BothSideCheck.gt.2) then ! both axes do i=1,num_bcCheckBin Deficiency_bin(num_bcCheckBin+i) = num_InitBin - bin_ParNumber(num_Bins+1-i) if ((i.eq.1).and.(PropertyNode.eq.1)) then Deficiency_bin(num_bcCheckBin+i) = Deficiency_bin(num_bcCheckBin+i)*0.5 endif deficiency = deficiency + Deficiency_bin(num_bcCheckBin+i) enddo endif if (deficiency.lt.0) then deficiency = 0 endif elseif ((Algo_bcCheck.eq.2).or.(Algo_bcCheck.eq.3)) then ! 2 constant net injection + Out - only first bin 224 if (AccumInjParL.ge.1) then ! 3 constant injection deficiency = deficiency + 1 AccumInjParL = AccumInjParL - 1 goto 224 endif deficiencyLeft = deficiency if (BothSideCheck.gt.2) then ! both axes 227 if (AccumInjParR.ge.1) then ! 3 compensate constant current injection - only first bin deficiency = deficiency + 1 AccumInjParR = AccumInjParR - 1 goto 227 endif endif else ! constant injection endif if (deficiency.gt.nonUseParNum) then goto 226 elseif (deficiency.ge.1) then do i = 1, deficiency if ((BothSideCheck.ge.3).and.(i.gt.deficiencyLeft)) then call DecideAlBin(box_coor(2,1),nAlBin,LowAlBound,UpAlBound) else call DecideAlBin(box_coor(1,1),nAlBin,LowAlBound,UpAlBound) endif addInx = nonUseList(i) 228 call random_number(harvest = RandomNum) do while(RandomNum.eq.0) call random_number(harvest = RandomNum) enddo ebin = 1 do j = 2, num_EnergyLevel if ((RandomNum.gt.RelCumFreqDist(nAlBin,j-1)).and.(RandomNum.le.RelCumFreqDist(nAlBin,j))) then ebin = j endif enddo if (ebin.eq.1) then binFreq = RelCumFreqDist(nAlBin,ebin) e = StepE_eV*(RandomNum/RelCumFreqDist(nAlBin,ebin)) !(log(RelCumFreqDist(ebin))/log(RandomNum)) else binFreq = RelCumFreqDist(nAlBin,ebin) - RelCumFreqDist(nAlBin,ebin-1) e = StepE_eV*real((ebin-1),kind = 8) + StepE_eV *(1.-(log(RelCumFreqDist(nAlBin,ebin))-log(RandomNum))/(log(RelCumFreqDist(nAlBin,ebin))-log(RelCumFreqDist(nAlBin,ebin-1)))) endif e = BCEne_factor*e ! boundary Test BCPos_factor, BCEne_factor ! call random_number(harvest = RandomNum) call random_number(harvest = RandomNum) if (Use_prerun.ne.0) then valley = 1 if (RandomNum.le.PRValley(1)) then valley = 1 elseif (RandomNum.le.PRValley(2)) then valley = 2 elseif (RandomNum.le.PRValley(3)) then valley = 3 endif else valley = 1 ! this is for GaAs materials endif k=semh(nAlBin,valley)*sqrt(e*(1.+ nonParabolPara(nAlBin,valley)*e)) ! Initial wavevector call random_number(harvest = RandomNum) fai=2.*pi*RandomNum call random_number(harvest = RandomNum) ct=1.-2.*RandomNum st=sqrt(1.-ct*ct) kx=k*st*cos(fai) ky=k*st*sin(fai) kz=k*ct 225 call random_number(harvest = RandomNum) if(RandomNum.le.1.e-5) goto 225 tc = -(log(RandomNum))*tau_max(nAlBin,valley) ! random free flight time.. ! vel_InitBias velocity bias... for initial current... 1. v 2.add 3. new k vector, new e. .. FactorEne, FactorVel Temp_wvk(1) = kx Temp_wvk(2) = ky Temp_wvk(3) = kz FactorVel = 1./(1.+nonParabolPara2(nAlBin,valley)*e) do j = 1, DIM Temp_vel(j) = hbar*Temp_wvk(j)*FactorVel/eff_mass(nAlBin,valley) enddo if ((BothSideCheck.eq.1).or.(BothSideCheck.eq.3)) then Temp_vel(1) = Temp_vel(1)+ vel_InitBias elseif ((BothSideCheck.eq.4).and.(i.le.deficiencyLeft)) then Temp_vel(1) = Temp_vel(1)+ vel_InitBias endif if ((BothSideCheck.ge.3).and.(i.gt.deficiencyLeft)) then ! 0529 if (Temp_vel(1).ge.0) goto 228 ! velocity check - right injected else if (Temp_vel(1).le.0) goto 228 ! velocity check - injected through left boundary endif ! new wvk Temp_wvk(1) = eff_mass(nAlBin,valley)*Temp_vel(1)/hbar/FactorVel sq_wvk_sum = Temp_wvk(1)*Temp_wvk(1) + Temp_wvk(2)*Temp_wvk(2) + Temp_wvk(3)*Temp_wvk(3) abs_wvk = sqrt(sq_wvk_sum) FactorEne = hhem(nAlBin,valley)*sq_wvk_sum ! new energy e = 2*FactorEne/(1.+sqrt(1.+nonParabolPara4(nAlBin,valley)*FactorEne)) ! Map particle atributes Par_wvk(addInx,1) = Temp_wvk(1) Par_wvk(addInx,2) = Temp_wvk(2) Par_wvk(addInx,3) = Temp_wvk(3) Par_tau(addInx) = tc Par_valley(addInx) = valley Par_ene(addInx) = e Par_use(addInx) = 1 Par_Al(addInx) = nAlBin call random_number(harvest = RandomNum) movingDist = Temp_vel(1)*(time_step) if ((BCoption.eq.1).or.(BCoption.eq.3)) then RandomNum = 1.0 elseif ((BCoption.eq.5).or.(BCoption.eq.7)) then RandomNum = 1.0 endif if ((BothSideCheck.ge.3).and.(i.gt.deficiencyLeft)) then Par_pos(addInx,1) = RandomNum*movingDist + box_coor(2,1) ! or equal test 1 ! Par_pos(addInx,1) = BCPos_factor*movingDist + box_coor(2,1) ! or equal test 2 ! Par_pos(addInx,1) = 0.5*movingDist + box_coor(2,1) ! or equal test 3 EneBoundR = EneBoundR - Par_ene(addInx) else Par_pos(addInx,1) = RandomNum*movingDist + box_coor(1,1) ! or equal !Par_pos(addInx,1) = BCPos_factor*movingDist + box_coor(1,1) ! or equal ! Par_pos(addInx,1) = 0.5*movingDist + box_coor(1,1) ! or equal EneBoundL = EneBoundL - Par_ene(addInx) endif call random_number(harvest = RandomNum) Par_pos(addInx,2) = box_len(2)*RandomNum + box_coor(1,2) call random_number(harvest = RandomNum) Par_pos(addInx,3) = box_len(3)*RandomNum + box_coor(1,3) num_Carriers = num_Carriers + 1 EneBound = EneBound - Par_ene(addInx) AccumEneBound = AccumEneBound - Par_ene(addInx) if ((BothSideCheck.ge.3).and.(i.gt.deficiencyLeft)) then ParInNumBoundR = ParInNumBoundR + 1 AccumParInNumBoundR = AccumParInNumBoundR + 1 ParNumBoundPlus = ParNumBoundPlus - 1 AccumParNumBoundPlus = AccumParNumBoundPlus - 1 else ParInNumBoundL = ParInNumBoundL + 1 AccumParInNumBoundL = AccumParInNumBoundL + 1 ParNumBoundPlus = ParNumBoundPlus + 1 AccumParNumBoundPlus = AccumParNumBoundPlus + 1 endif enddo endif return 226 continue print*,'Reach To Max Carrier number!!' stop end