subroutine Update_Tdist() use particles use scattering use simulation ! e-p from MC simulation integer nBin,Temp_Bin_Lower,Temp_Bin_Upper ! initial and final valley real ( kind = 8 ) :: fTemperAP,fUnitVolAPEne,UnitVolAPEne,Time_step_PhononUpdate, exp_pA, f_pA, ppInterUpRate, ppInterDnRate, epInterEne, epInterEneTO, epInterEneWm3, epInterEneTOWm2 ! local variouble if (QuasiBoundConsid.eq.1)then Temp_Bin_Lower = num_boundBin+1 Temp_Bin_Upper = num_Bins-num_boundBin else Temp_Bin_Lower = 1 Temp_Bin_Upper = num_Bins endif do nBin = Temp_Bin_Lower,Temp_Bin_Upper Time_step_PhononUpdate = time_step*real(TempUpdateFreq,kind = 8) epInterEne = 0 ! pp interaction rate.. ! New_NpLO = 1/(exp((Ep_LO*e_c)/(kB*Initial_T_pLO)) - 1) ! exp_pA = exp((0.5*Ep_LO*e_c)/(kB*Initial_T_pAO)) ! f_pa(t) = 1/(exp_pA - 1) Old_NpLO*(-exp(-gammadotUp*StepTps)+exp(-gammadotDn*StepTps)) ! gammadotDn = hbar*Gru*Gru*R*(omega_AP^2)*(omega_oLO^3)*(f_pa(t-1)+1)*(f_pa(t-1)+1)*Old_NpLO/(2*pi**(u_pA^5)); gammadotUp = hbar*Gru*Gru*R*(omega_AP^2)*(omega_oLO^3)*f_pa(t-1)*f_pa(t-1)*(Old_NpLO+1)/(2*pi*rho_GaAs*(u_pA^5)); exp_pA = exp((0.5*RepresentEpLO*e_c)/(kB*Bin_TemperAP(nBin))) f_pA = 1/(exp_pA - 1) ppInterUpRate = (Bin_OccupOP(nBin)+1)*f_pA*f_pA*MatElement ! Energy to Acoustic LO to AO for timestep ppInterDnRate = Bin_OccupOP(nBin)*(f_pA+1)*(f_pA+1)*MatElement if (LOmodeOccPhonSim.eq.1) then ! OP mode energy.. Occupation*number density of phonon* energy J/m^2 epInterEne = epInterEne+real((TempVar_BinScatCount(nBin,2,1)-TempVar_BinScatCount(nBin,3,1)),kind = 8)*polar_Ep_G+real((TempVar_BinScatCount(nBin,2,2)-TempVar_BinScatCount(nBin,3,2)),kind = 8)*polar_Ep_L+real((TempVar_BinScatCount(nBin,2,3)-TempVar_BinScatCount(nBin,3,3)),kind = 8)*polar_Ep_X ! absorption - emission net absorption .. phonon to electron.. optical energy decrease epInterEneTO = 0 ! LOmodeOccPhonSim.eq.1 epInterEneTOWm2 = 0 epInterEneTO = epInterEneTO+real((TempVar_BinScatCount(nBin,4,1)-TempVar_BinScatCount(nBin,5,1)),kind = 8)*IntVal_Ep_G_L+real((TempVar_BinScatCount(nBin,6,1)-TempVar_BinScatCount(nBin,7,1)),kind = 8)*IntVal_Ep_G_X epInterEneTO = epInterEneTO+real((TempVar_BinScatCount(nBin,4,2)-TempVar_BinScatCount(nBin,5,2)),kind = 8)*IntVal_Ep_L_G+real((TempVar_BinScatCount(nBin,6,2)-TempVar_BinScatCount(nBin,7,2)),kind = 8)*IntVal_Ep_L_L+real((TempVar_BinScatCount(nBin,8,2)-TempVar_BinScatCount(nBin,9,2)),kind = 8)*IntVal_Ep_L_X epInterEneTO = epInterEneTO+real((TempVar_BinScatCount(nBin,4,3)-TempVar_BinScatCount(nBin,5,3)),kind = 8)*IntVal_Ep_X_G+real((TempVar_BinScatCount(nBin,6,3)-TempVar_BinScatCount(nBin,7,3)),kind = 8)*IntVal_Ep_X_L+real((TempVar_BinScatCount(nBin,8,3)-TempVar_BinScatCount(nBin,9,3)),kind = 8)*IntVal_Ep_X_X epInterEneTOWm2 = (epInterEneTO*e_c)*n_d*BinSize/real(num_InitNode,kind = 8)/Time_step_PhononUpdate else epInterEne = epInterEne+real((TempVar_BinScatCount(nBin,2,1)-TempVar_BinScatCount(nBin,3,1)),kind = 8)*polar_Ep_G+real((TempVar_BinScatCount(nBin,4,1)-TempVar_BinScatCount(nBin,5,1)),kind = 8)*IntVal_Ep_G_L+real((TempVar_BinScatCount(nBin,6,1)-TempVar_BinScatCount(nBin,7,1)),kind = 8)*IntVal_Ep_G_X epInterEne = epInterEne+real((TempVar_BinScatCount(nBin,2,2)-TempVar_BinScatCount(nBin,3,2)),kind = 8)*polar_Ep_L+real((TempVar_BinScatCount(nBin,4,2)-TempVar_BinScatCount(nBin,5,2)),kind = 8)*IntVal_Ep_L_G+real((TempVar_BinScatCount(nBin,6,2)-TempVar_BinScatCount(nBin,7,2)),kind = 8)*IntVal_Ep_L_L+real((TempVar_BinScatCount(nBin,8,2)-TempVar_BinScatCount(nBin,9,2)),kind = 8)*IntVal_Ep_L_X epInterEne = epInterEne+real((TempVar_BinScatCount(nBin,2,3)-TempVar_BinScatCount(nBin,3,3)),kind = 8)*polar_Ep_X+real((TempVar_BinScatCount(nBin,4,3)-TempVar_BinScatCount(nBin,5,3)),kind = 8)*IntVal_Ep_X_G+real((TempVar_BinScatCount(nBin,6,3)-TempVar_BinScatCount(nBin,7,3)),kind = 8)*IntVal_Ep_X_L+real((TempVar_BinScatCount(nBin,8,3)-TempVar_BinScatCount(nBin,9,3)),kind = 8)*IntVal_Ep_X_X epInterEneTOWm2 = 0 endif if ((nBin.eq.1).or.(nBin.eq.num_Bins)) then boundFactor = 0.5 else boundFactor = 1.0 endif ! p-p interaction from the previous step Up-Dn .. Thus OP increses ppInterEneWm2(nBin) = (-exp(-ppInterUpRate*Time_step_PhononUpdate)+exp(-ppInterDnRate*Time_step_PhononUpdate))*Bin_OccupOP(nBin)*boundFactor*ConstForEneOp/Time_step_PhononUpdate ! Energy to Acoustic LO to AO for timestep epInterEneWm3 = (epInterEne*e_c)*n_d/real(num_InitNode,kind = 8)/Time_step_PhononUpdate ! W/m^3 - BinSize*TA_bin_ParNumber(j) RepresentEpLO*e_c*BinSize/primV epInterEneWm2(nBin) = epInterEneWm3*BinSize ! Total energy absorption in bin W/m^2 ! prevOPEne(nBin) = Bin_OccupOP(nBin)*RepresentEpLO*e_c*BinSize/primV ! only LO phonon considered if (nBin.eq.Temp_Bin_Lower) then ! first ConstForEneOp = RepresentEpLO*e_c*BinSize/primV BinOPEne(nBin) = (qsupplyO-BinQsupO(nBin)-epInterEneWm2(nBin)+ppInterEneWm2(nBin))*Time_step_PhononUpdate + BinOPEne(nBin) ! energy J/m2 BinAPEne(nBin) = (qsupplyA-BinQsupA(nBin)-epInterEneTOWm2-ppInterEneWm2(nBin))*Time_step_PhononUpdate + BinAPEne(nBin) else BinOPEne(nBin) = (BinQsupO(nBin-1)-BinQsupO(nBin)-epInterEneWm2(nBin)+ ppInterEneWm2(nBin))*Time_step_PhononUpdate + BinOPEne(nBin)! BinAPEne(nBin) = (BinQsupA(nBin-1)-BinQsupA(nBin)-epInterEneTOWm2-ppInterEneWm2(nBin))*Time_step_PhononUpdate + BinAPEne(nBin) endif Bin_OccupOP(nBin) = BinOPEne(nBin)/(boundFactor*ConstForEneOp) ! occupancy from total energy Bin_TemperOP(nBin) = RepresentEpLO*e_c/(kB*(log(Bin_OccupOP(nBin)+1)-log(Bin_OccupOP(nBin)))) fUnitVolAPEne = BinAPEne(nBin)/(boundFactor*BinSize) call APEneToTempCal(fUnitVolAPEne,fTemperAP) Bin_TemperAP(nBin) = fTemperAP if ((Bin_TemperAP(nBin).lt.APTableTempLow).or.(Bin_TemperAP(nBin).gt.APTableTempHigh)) then ! for error check print*,' T = ', Bin_TemperAP(nBin) print*,' Bin = ',nBin print*,' Time = ',time endif enddo ! Acoustic PhononEneAP = PhononEneAP() - ppInterEne(nBin) + HeatFlowInA - HeatFlowOutA ! call PhononEneToOccTemp() BinQsupA(nBin) = BinQsupA(nBin-1) + ppInterEne(nBin) if (BoundQSetting.eq.6) then ! no more heat source and sink BinQsupA(Temp_Bin_Upper) = -(T_room - Bin_TemperAP(Temp_Bin_Upper))/ExternalResAP if (ExternalResOP.ne.0) then BinQsupO(Temp_Bin_Upper) = -(T_room - Bin_TemperOP(Temp_Bin_Upper))/ExternalResOP else BinQsupO(Temp_Bin_Upper) = 0 endif else BinQsupA(Temp_Bin_Upper) = BinQsupA(Temp_Bin_Upper-1) BinQsupO(Temp_Bin_Upper) = BinQsupO(Temp_Bin_Upper-1) endif do nBin = (Temp_Bin_Lower), (Temp_Bin_Upper-1) BinQsupA(nBin) = -k_pA*(Bin_TemperAP(nBin+1) - Bin_TemperAP(nBin))/BinSize BinQsupO(nBin) = -k_pO*(Bin_TemperOP(nBin+1) - Bin_TemperOP(nBin))/BinSize enddo do nBin = Temp_Bin_Lower,Temp_Bin_Upper !col_len = RCutoff_SiIn BinQsupTotal(nBin) = BinQsupA(nBin) + BinQsupO(nBin) enddo ! making new table for each position if (QuasiBoundConsid.eq.1)then do nBin = 1,num_boundBin !col_len = RCutoff_SiIn BinQsupA(nBin) = qsupplyA BinQsupO(nBin) = qsupplyO BinQsupTotal(nBin) = qsupply Bin_TemperAP(nBin) = Bin_TemperAP(Temp_Bin_Lower)+ qsupplyA*BinSize*(Temp_Bin_Lower-nBin)/k_pA Bin_TemperOP(nBin) = Bin_TemperOP(Temp_Bin_Lower)+ qsupplyO*BinSize*(Temp_Bin_Lower-nBin)/k_pO Bin_OccupOP(nBin) = 1/(exp((RepresentEpLO*e_c)/(kB*Bin_TemperOP(nBin))) - 1) call APTempToEneCal(Bin_TemperAP(nBin),UnitVolAPEne) if ((nBin.eq.1).or.(nBin.eq.num_Bins)) then ! first BinAPEne(nBin) = UnitVolAPEne*(0.5*BinSize) BinOPEne(nBin) = RepresentEpLO*e_c*Bin_OccupOP(nBin)*(0.5*BinSize)/Default_primV else ! middle of first last bin BinAPEne(nBin) = UnitVolAPEne*BinSize BinOPEne(nBin) = RepresentEpLO*e_c*Bin_OccupOP(nBin)*BinSize/Default_primV endif enddo do nBin = Temp_Bin_Upper+1,num_Bins !col_len = RCutoff_SiIn BinQsupA(nBin) = 0 BinQsupO(nBin) = 0 BinQsupTotal(nBin) = 0 BinAPEne(nBin) = BinAPEne(Temp_Bin_Upper) BinOPEne(nBin) = BinOPEne(Temp_Bin_Upper) Bin_OccupOP(nBin) = Bin_OccupOP(Temp_Bin_Upper) Bin_TemperAP(nBin) = Bin_TemperAP(Temp_Bin_Upper) Bin_TemperOP(nBin) = Bin_TemperOP(Temp_Bin_Upper) enddo endif do nBin = 1,num_BinsComb !col_len = RCutoff_SiIn call TempVar_sc_table(nBin) enddo return end ! if (PhononSim.eq.1) then ! left constant fixed ! if (nBin.eq.1) then ! first ! Bin_TemperAP(nBin) = TpA_Hot ! else ! Bin_TemperAP(nBin) = Bin_TemperAP(nBin-1) - BinSize*BinQsupA(nBin-1)/k_pA ! endif ! elseif (PhononSim.eq.2) then ! from RightTemp Fixed. ! if (nBin.eq.1) then ! first ! Bin_TemperAP(num_Bins+1-nBin) = TpA_Cold ! else ! Bin_TemperAP(num_Bins+1-nBin) = Bin_TemperAP(num_Bins+2-nBin) + BinSize*BinQsupA(num_Bins+2-nBin)/k_pA ! endif ! endif subroutine APEneToTempCal(fUnitVolAPEne,fTemperAP) use particles use scattering use simulation real ( kind = 8 ) :: fUnitVolAPEne real ( kind = 8 ) :: RefdiffEne,fTemperAP ! local variouble integer:: repeat, preInx, checkGap, prevGap, j checkGap = num_APTable ! first check - within range or not if (fUnitVolAPEne.ge.APTableEnergy(num_APTable+1)) then ! not accurate regime RefdiffEne = (fUnitVolAPEne - APTableEnergy(num_APTable+1))/(APTableEnergy(num_APTable+1) - APTableEnergy(num_APTable)) fTemperAP = APTableTempHigh + step_APTemp*RefdiffEne elseif (fUnitVolAPEne.lt.APTableEnergy(1)) then RefdiffEne = (APTableEnergy(1)-fUnitVolAPEne)/(APTableEnergy(2)- APTableEnergy(1)) ! not accurate regime fTemperAP = APTableTempLow - step_APTemp*RefdiffEne else preInx = 0 do repeat = 1, num_APmidCheck prevGap = checkGap checkGap = int(checkGap/10) if (checkGap.le.1) then preInx = preInx + (SelectedIdx(repeat-1)-1)*prevGap do j = 1, prevGap ! lowerNum = preInx + j ! upperNum = preInx + j+1 ! if ((fUnitVolAPEne.ge.APTableEnergy(lowerNum)).and.(fUnitVolAPEne.lt.APTableEnergy(upperNum))) then SelectedIdx(repeat) = j endif enddo else if (repeat.gt.1) then preInx = preInx + (SelectedIdx(repeat-1)-1)*prevGap endif do j = 1, 10 ! if (repeat.eq.1) then lowerNum = (j-1)*checkGap+1 ! upperNum = j*checkGap+1 ! else lowerNum = preInx + (j-1)*checkGap+1 ! upperNum = preInx + j*checkGap+1 ! endif if ((fUnitVolAPEne.ge.APTableEnergy(lowerNum)).and.(fUnitVolAPEne.lt.APTableEnergy(upperNum))) then SelectedIdx(repeat) = j endif enddo ! do j = 1, 10 ! endif ! enddo ! ! interpolation preInx = preInx + SelectedIdx(repeat-1) RefdiffEne = (fUnitVolAPEne- APTableEnergy(preInx))/(APTableEnergy(preInx+1)- APTableEnergy(preInx)) fTemperAP = APTableTemper(preInx+1)*RefdiffEne+APTableTemper(preInx)*(1.0-RefdiffEne) endif return end subroutine APTempToEneCal(TemperAP,UnitVolAPEne) use particles use scattering use simulation real ( kind = 8 ) :: Refdif, UnitVolAPEne, TemperAP ! local variouble integer:: TempGuess,upperInx,lowerInx, i, foundIdx if (TemperAP.ge.APTableTempHigh) then Refdif = (TemperAP - APTableTempHigh)/step_APTemp UnitVolAPEne = APTableEnergy(num_APTable+1) + Refdif*(APTableEnergy(num_APTable+1)-APTableEnergy(num_APTable)) elseif (TemperAP.lt.APTableTempLow) then Refdif = (APTableTempLow - TemperAP)/step_APTemp UnitVolAPEne = APTableEnergy(1) - Refdif*(APTableEnergy(2)-APTableEnergy(1)) else TempGuess = int((TemperAP-APTableTempLow)/step_APTemp) +1 if (TempGuess.ge.(num_APTable-2)) then upperInx = num_APTable else upperInx = TempGuess + 2 endif if (TempGuess.lt.3) then lowerInx = 1 else lowerInx = TempGuess - 2 endif do i = lowerInx,upperInx if ((TemperAP.ge.APTableTemper(i)).and.(TemperAP.lt.APTableTemper(i+1))) then foundIdx = i endif enddo Refdif = (TemperAP- APTableTemper(foundIdx))/(APTableTemper(foundIdx+1)- APTableTemper(foundIdx)) UnitVolAPEne = Refdif*APTableEnergy(foundIdx+1) + (1.0-Refdif)*APTableEnergy(foundIdx) endif return end subroutine InitializeAPTable() use particles use scattering use simulation integer:: DataPartMin,DataPartMax, numT, i real ( kind = 8 ) :: integTPhononEne,tempAPTemp,tempTstar,tempTotalHeatCap, expLO,fpLO,tempLOphononEne ! in input.txt determine.. APTableTempLow = 100 APTableTempHigh = 600 num_APTable = 10000 num_APmidCheck = 4 APParaCapRT = 322 APParaCapZ = 50 APParaBeta = 1.6 ! Global step_APTemp *input *allo APTableTemper, APTableHeatCp, APTableEnergy integer, dimension(:), allocatable :: ! Local DataPartMin DataPartMax tempAPTemp tempTstar step_APTemp = (APTableTempHigh-APTableTempLow)/real(num_APTable,kind = 8) ! calculate step allocate ( SelectedIdx(num_APmidCheck) ) allocate ( APTableTemper(num_APTable+1) ) allocate ( APTableHeatCp(num_APTable+1) ) allocate ( APTableEnergy(num_APTable+1) ) SelectedIdx(:) = 0 APTableTemper(:) = 0.0 APTableHeatCp(:) = 0.0 APTableEnergy(:) = 0.0 DataPartMin = int(APTableTempLow/step_APTemp) DataPartMax = int(APTableTempHigh/step_APTemp) numT = 0 integTPhononEne = 0 do i = 1, DataPartMax tempAPTemp = real(i,kind = 8)*step_APTemp tempTstar = (tempAPTemp/300.0)**APParaBeta tempTotalHeatCap = APParaCapRT + APParaCapZ*(tempTstar-1.0)/(tempTstar+APParaCapZ/APParaCapRT) ! J/kg-K integTPhononEne = integTPhononEne + tempTotalHeatCap*step_APTemp ! J/kg * rho_GaAs ! dfpLOdT = RepresentEpLO*expLO*e_c/((expLO-1)*(expLO-1)*kB*tempAPTemp*tempAPTemp) dUdT(i) = dfpLOdT(i)*EpLO*e_c/primM integForU = integForU + dUdT(i)*step_APTemp if (i.ge.DataPartMin) then numT = numT + 1 expLO = exp(RepresentEpLO*e_c/(kB*tempAPTemp)) fpLO = 1/(expLO-1) if (LOmodeOccPhonSim.eq.1) then tempLOphononEne = fpLO*RepresentEpLO*e_c/Default_primV ! J/m3 else tempLOphononEne = fpLO*RepresentEpLO*3.0*e_c/Default_primV ! J/m3 endif APTableTemper(numT) = tempAPTemp APTableHeatCp(numT) = tempTotalHeatCap APTableEnergy(numT) = integTPhononEne*Default_rho_GaAs - tempLOphononEne endif enddo return end !if (BoundQSetting.gt.3) then ! if (BoundQSetting.eq.4) then ! balance q ! if (LOmodeOccPhonSim.eq.1) then ! BinQsupA(num_Bins) = BinQsupA(num_Bins-1)-epInterEneTOWm2-ppInterEneWm2(num_Bins)*0.5 ! else ! BinQsupA(num_Bins) = BinQsupA(num_Bins-1)-ppInterEneWm2(num_Bins)*0.5 ! endif ! BinQsupO(num_Bins) = BinQsupO(num_Bins-1)-epInterEneWm2(num_Bins)+ppInterEneWm2(num_Bins)*0.5 ! elseif (BoundQSetting.eq.5) then ! no more heat source and sink ! BinQsupA(num_Bins) = BinQsupA(num_Bins-1) ! BinQsupO(num_Bins) = BinQsupO(num_Bins-1) ! elseif (BoundQSetting.eq.6) then ! no more heat source and sink ! BinQsupA(num_Bins) = -(T_room - Bin_TemperAP(num_Bins))/ExternalResAP ! if (ExternalResOP.ne.0) then ! BinQsupO(num_Bins) = -(T_room - Bin_TemperOP(num_Bins))/ExternalResOP ! else ! BinQsupO(num_Bins) = 0 !!! endif ! endif ! do nBin = 1,(num_Bins-1) ! BinSize ! BinQsupA(nBin) = -k_pA*(Bin_TemperAP(nBin+1) - Bin_TemperAP(nBin))/BinSize ! BinQsupO(nBin) = -k_pO*(Bin_TemperOP(nBin+1) - Bin_TemperOP(nBin))/BinSize ! enddo ! else ! do nBin = 1,num_Bins ! BinSize ! if (nBin.ne.num_Bins) then ! first ! BinQsupA(nBin) = -k_pA*(Bin_TemperAP(nBin+1) - Bin_TemperAP(nBin))/BinSize ! BinQsupO(nBin) = -k_pO*(Bin_TemperOP(nBin+1) - Bin_TemperOP(nBin))/BinSize ! else ! consider 3 options... 1. constant T 2. Consistent Q 3. balance Q ! if (BoundQSetting.eq.1) then ! balance q ! if (LOmodeOccPhonSim.eq.1) then ! BinQsupA(nBin) = BinQsupA(nBin-1)-epInterEneTOWm2-ppInterEneWm2(nBin)*0.5 ! else ! BinQsupA(nBin) = BinQsupA(nBin-1)-ppInterEneWm2(nBin)*0.5 ! endif ! BinQsupO(nBin) = BinQsupO(nBin-1)-epInterEneWm2(nBin)+ppInterEneWm2(nBin)*0.5 ! elseif (BoundQSetting.eq.2) then ! no more heat source and sink ! BinQsupA(nBin) = BinQsupA(nBin-1) ! BinQsupO(nBin) = BinQsupO(nBin-1) ! else ! constant T ! BinQsupA(nBin) = -k_pA*(RightBoundTempAP - Bin_TemperAP(nBin))/BinSize ! BinQsupO(nBin) = -k_pO*(RightBoundTempOP - Bin_TemperOP(nBin))/BinSize ! endif ! endif ! enddo ! endif