subroutine FindQuasiSteadyTdist() ! use particles use scattering use simulation implicit none integer nBin,TempQuasi_Bin_Lower,TempQuasi_Bin_Upper,step, err real ( kind = 8 ):: Quasi_OutStep,QuasiTimeAvg,QuasiTime,boundFactor,f_pA,exp_pA,ppInterUpRate,ppInterDnRate,TempVariation,fUnitVolAPEne,fTemperAP,UnitVolAPEne real ( kind = 8 ), dimension(:), allocatable :: QuasiPrev_TemperAP, QuasiPrev_TemperOP real ( kind = 8 ), dimension(:), allocatable :: TempAvg_BinQsupA,TempAvg_BinQsupO,TempAvg_BinQsupTotal,TempAvg_BinAPEne,TempAvg_BinOPEne,TempAvg_Bin_OccupOP real ( kind = 8 ), dimension(:), allocatable :: TempAvg_Bin_TemperAP,TempAvg_Bin_TemperOP,TempAvg_Quasi_peOWm2,TempAvg_Quasi_peTOWm2,TempAvg_Quasi_peLOWm2,TempAvg_Quasi_ppUpWm2 allocate ( QuasiPrev_TemperAP(num_Bins) ) allocate ( QuasiPrev_TemperOP(num_Bins) ) QuasiPrev_TemperAP(1:num_Bins) = 0.0d0 QuasiPrev_TemperOP(1:num_Bins) = 0.0d0 allocate ( TempAvg_BinQsupA(num_Bins) ) allocate ( TempAvg_BinQsupO(num_Bins) ) allocate ( TempAvg_BinQsupTotal(num_Bins) ) allocate ( TempAvg_BinAPEne(num_Bins) ) allocate ( TempAvg_BinOPEne(num_Bins) ) allocate ( TempAvg_Bin_OccupOP(num_Bins) ) allocate ( TempAvg_Bin_TemperAP(num_Bins) ) allocate ( TempAvg_Bin_TemperOP(num_Bins) ) allocate ( TempAvg_Quasi_peOWm2(num_Bins) ) allocate ( TempAvg_Quasi_peTOWm2(num_Bins) ) allocate ( TempAvg_Quasi_peLOWm2(num_Bins) ) allocate ( TempAvg_Quasi_ppUpWm2(num_Bins) ) TempAvg_BinQsupA(1:num_Bins) = 0.0d0 TempAvg_BinQsupO(1:num_Bins) = 0.0d0 TempAvg_BinQsupTotal(1:num_Bins) = 0.0d0 TempAvg_BinAPEne(1:num_Bins) = 0.0d0 TempAvg_BinOPEne(1:num_Bins) = 0.0d0 TempAvg_Bin_OccupOP(1:num_Bins) = 0.0d0 TempAvg_Bin_TemperAP(1:num_Bins) = 0.0d0 TempAvg_Bin_TemperOP(1:num_Bins) = 0.0d0 TempAvg_Quasi_peOWm2(1:num_Bins) = 0.0d0 TempAvg_Quasi_peTOWm2(1:num_Bins) = 0.0d0 TempAvg_Quasi_peLOWm2(1:num_Bins) = 0.0d0 TempAvg_Quasi_ppUpWm2(1:num_Bins) = 0.0d0 step = 0 Quasi_OutStep = 0.0 call InitQuasiSteadyTdist() ! Initialize T distribution if (QuasiBoundConsid.eq.1)then TempQuasi_Bin_Lower = num_boundBin+1 TempQuasi_Bin_Upper = num_Bins-num_boundBin else TempQuasi_Bin_Lower = 1 TempQuasi_Bin_Upper = num_Bins endif TempVariation = 10*MaxTempChange open (unit = 63, file = 'QuasiSteadyTempAT.txt', status = 'replace', action = 'write' ) open (unit = 64, file = 'QuasiSteadyTempOT.txt', status = 'replace', action = 'write' ) write(63,'(2x,a,a)'), 'Step Time ', 'QBin_TemperOP(1:num_Bins)' ! last Step... write(64,'(2x,a,a)'), 'Step Time ', 'QBin_TemperOP(1:num_Bins)' ! last Step.. do while((QuasiTime.le.QuasiTimeMax).or.(TempVariation.gt.MaxTempChange)) ! to find step = step + 1 QuasiTime = QuasiTimeStep*real(step,kind = 8) TempVariation = 0.0 ! smallest do nBin = TempQuasi_Bin_Lower, TempQuasi_Bin_Upper QuasiPrev_TemperAP(nBin) = Bin_TemperAP(nBin) QuasiPrev_TemperOP(nBin) = Bin_TemperOP(nBin) if ((nBin.eq.1).or.(nBin.eq.num_Bins)) then boundFactor = 0.5 else boundFactor = 1.0 endif exp_pA = exp((0.5*RepresentEpLO*e_c)/(kB*QuasiPrev_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 Quasi_ppUpWm2(nBin) = (-exp(-ppInterUpRate*QuasiTimeStep)+exp(-ppInterDnRate*QuasiTimeStep))*Bin_OccupOP(nBin)*ConstForEneOp/QuasiTimeStep ! Energy to Acoustic LO to AO for timestep if (QuasiSdotMode.eq.1) then ! changing Quasi_peLOWm2(nBin) = PreLoop_Bin_EpLOWm2(nBin)*(1+Bin_OccupOP(nBin))/(1+PreLoop_Bin_OccupOP(nBin)) Quasi_peTOWm2(nBin) = PreLoop_Bin_EpTOWm2(nBin)*(1+Bin_OccupOP(nBin))/(1+PreLoop_Bin_OccupOP(nBin)) Quasi_peOWm2(nBin) = Quasi_peLOWm2(nBin) + Quasi_peTOWm2(nBin) endif if (nBin.eq.(TempQuasi_Bin_Lower)) then ! first ConstForEneOp = RepresentEpLO*e_c*BinSize/primV if (LOmodeOccPhonSim.eq.1) then ! change in internal energy BinOPEne(nBin) = (qsupplyO-BinQsupO(nBin)-Quasi_peLOWm2(nBin)+Quasi_ppUpWm2(nBin)*boundFactor)*QuasiTimeStep + BinOPEne(nBin) ! energy J/m2 BinAPEne(nBin) = (qsupplyA-BinQsupA(nBin)-Quasi_peTOWm2(nBin)-Quasi_ppUpWm2(nBin)*boundFactor)*QuasiTimeStep + BinAPEne(nBin) else BinOPEne(nBin) = (qsupplyO-BinQsupO(nBin)-Quasi_peOWm2(nBin)+Quasi_ppUpWm2(nBin)*boundFactor)*QuasiTimeStep + BinOPEne(nBin) ! energy J/m2 BinAPEne(nBin) = (qsupplyA-BinQsupA(nBin)-Quasi_ppUpWm2(nBin)*boundFactor)*QuasiTimeStep + BinAPEne(nBin) endif else if (LOmodeOccPhonSim.eq.1) then ! change in internal energy BinOPEne(nBin) = (BinQsupO(nBin-1)-BinQsupO(nBin)-Quasi_peLOWm2(nBin)+Quasi_ppUpWm2(nBin)*boundFactor)*QuasiTimeStep + BinOPEne(nBin) ! energy J/m2 BinAPEne(nBin) = (BinQsupA(nBin-1)-BinQsupA(nBin)-Quasi_peTOWm2(nBin)-Quasi_ppUpWm2(nBin)*boundFactor)*QuasiTimeStep + BinAPEne(nBin) else BinOPEne(nBin) = (BinQsupO(nBin-1)-BinQsupO(nBin)-Quasi_peOWm2(nBin)+Quasi_ppUpWm2(nBin)*boundFactor)*QuasiTimeStep + BinOPEne(nBin) ! energy J/m2 BinAPEne(nBin) = (BinQsupA(nBin-1)-BinQsupA(nBin)-Quasi_ppUpWm2(nBin)*boundFactor)*QuasiTimeStep + BinAPEne(nBin) endif endif if ((BinOPEne(nBin).le.0).or.(BinAPEne(nBin).le.0)) then err = 0 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 = ',QuasiTime ! endif if (abs(QuasiPrev_TemperAP(nBin)-Bin_TemperAP(nBin)).gt.TempVariation) then TempVariation = abs(QuasiPrev_TemperAP(nBin)-Bin_TemperAP(nBin)) endif if (abs(QuasiPrev_TemperOP(nBin)-Bin_TemperOP(nBin)).gt.TempVariation) then TempVariation = abs(QuasiPrev_TemperOP(nBin)-Bin_TemperOP(nBin)) endif enddo ! do nBin = TempQuasi_Bin_Lower, TempQuasi_Bin_Upper if (BoundQSetting.eq.6) then ! no more heat source and sink BinQsupA(TempQuasi_Bin_Upper) = -(T_room - Bin_TemperAP(TempQuasi_Bin_Upper))/ExternalResAP if (ExternalResOP.ne.0) then BinQsupO(TempQuasi_Bin_Upper) = -(T_room - Bin_TemperOP(TempQuasi_Bin_Upper))/ExternalResOP else BinQsupO(TempQuasi_Bin_Upper) = 0 endif else BinQsupA(TempQuasi_Bin_Upper) = BinQsupA(TempQuasi_Bin_Upper-1) BinQsupO(TempQuasi_Bin_Upper) = BinQsupO(TempQuasi_Bin_Upper-1) endif do nBin = (TempQuasi_Bin_Lower), (TempQuasi_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 = TempQuasi_Bin_Lower,TempQuasi_Bin_Upper !col_len = RCutoff_SiIn BinQsupTotal(nBin) = BinQsupA(nBin) + BinQsupO(nBin) enddo Quasi_OutStep = Quasi_OutStep+QuasiTimeStep if ((QuasiTime.ge.(0.1*QuasiTimeMax)).and.(Quasi_OutStep.ge.(0.01*QuasiTimeMax-QuasiTimeStep))) then write(63,'(1x,i8,e17.5,210f17.6)'),step,QuasiTime,Bin_TemperAP(TempQuasi_Bin_Lower:TempQuasi_Bin_Upper) write(64,'(1x,i8,e17.5,210f17.6)'),step,QuasiTime,Bin_TemperOP(TempQuasi_Bin_Lower:TempQuasi_Bin_Upper) Quasi_OutStep = 0.0 elseif ((QuasiTime.lt.(0.1*QuasiTimeMax)).and.(Quasi_OutStep.ge.(0.001*QuasiTimeMax-QuasiTimeStep))) then write(63,'(1x,i8,e17.5,210f17.6)'),step,QuasiTime,Bin_TemperAP(TempQuasi_Bin_Lower:TempQuasi_Bin_Upper) write(64,'(1x,i8,e17.5,210f17.6)'),step,QuasiTime,Bin_TemperOP(TempQuasi_Bin_Lower:TempQuasi_Bin_Upper) Quasi_OutStep = 0.0 endif enddo ! do while.. QuasiTimeStep QuasiTimeMax QuasiTime ! PreLoop_Bin_EpTOWm2(j) ! Initial Setup ... outside bound writing and intial setup for main simulation step = 0 do while(QuasiTimeAvg.le.(QuasiTimeMax*0.05)) ! to find average step = step + 1 QuasiTimeAvg = QuasiTimeStep*real(step,kind = 8) TempVariation = 0.0 ! smallest do nBin = TempQuasi_Bin_Lower, TempQuasi_Bin_Upper if ((nBin.eq.1).or.(nBin.eq.num_Bins)) then boundFactor = 0.5 else boundFactor = 1.0 endif 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 Quasi_ppUpWm2(nBin) = (-exp(-ppInterUpRate*QuasiTimeStep)+exp(-ppInterDnRate*QuasiTimeStep))*Bin_OccupOP(nBin)*ConstForEneOp/QuasiTimeStep ! Energy to Acoustic LO to AO for timestep if (QuasiSdotMode.eq.1) then ! changing Quasi_peLOWm2(nBin) = PreLoop_Bin_EpLOWm2(nBin)*(1+Bin_OccupOP(nBin))/(1+PreLoop_Bin_OccupOP(nBin)) Quasi_peTOWm2(nBin) = PreLoop_Bin_EpTOWm2(nBin)*(1+Bin_OccupOP(nBin))/(1+PreLoop_Bin_OccupOP(nBin)) Quasi_peOWm2(nBin) = Quasi_peLOWm2(nBin) + Quasi_peTOWm2(nBin) endif if (nBin.eq.(TempQuasi_Bin_Lower)) then ! first ConstForEneOp = RepresentEpLO*e_c*BinSize/primV if (LOmodeOccPhonSim.eq.1) then ! change in internal energy BinOPEne(nBin) = (qsupplyO-BinQsupO(nBin)-Quasi_peLOWm2(nBin)+Quasi_ppUpWm2(nBin)*boundFactor)*QuasiTimeStep + BinOPEne(nBin) ! energy J/m2 BinAPEne(nBin) = (qsupplyA-BinQsupA(nBin)-Quasi_peTOWm2(nBin)-Quasi_ppUpWm2(nBin)*boundFactor)*QuasiTimeStep + BinAPEne(nBin) else BinOPEne(nBin) = (qsupplyO-BinQsupO(nBin)-Quasi_peOWm2(nBin)+Quasi_ppUpWm2(nBin)*boundFactor)*QuasiTimeStep + BinOPEne(nBin) ! energy J/m2 BinAPEne(nBin) = (qsupplyA-BinQsupA(nBin)-Quasi_ppUpWm2(nBin)*boundFactor)*QuasiTimeStep + BinAPEne(nBin) endif else if (LOmodeOccPhonSim.eq.1) then ! change in internal energy BinOPEne(nBin) = (BinQsupO(nBin-1)-BinQsupO(nBin)-Quasi_peLOWm2(nBin)+Quasi_ppUpWm2(nBin)*boundFactor)*QuasiTimeStep + BinOPEne(nBin) ! energy J/m2 BinAPEne(nBin) = (BinQsupA(nBin-1)-BinQsupA(nBin)-Quasi_peTOWm2(nBin)-Quasi_ppUpWm2(nBin)*boundFactor)*QuasiTimeStep + BinAPEne(nBin) else BinOPEne(nBin) = (BinQsupO(nBin-1)-BinQsupO(nBin)-Quasi_peOWm2(nBin)+Quasi_ppUpWm2(nBin)*boundFactor)*QuasiTimeStep + BinOPEne(nBin) ! energy J/m2 BinAPEne(nBin) = (BinQsupA(nBin-1)-BinQsupA(nBin)-Quasi_ppUpWm2(nBin))*QuasiTimeStep + BinAPEne(nBin) endif 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 enddo ! do nBin = TempQuasi_Bin_Lower, TempQuasi_Bin_Upper if (BoundQSetting.eq.6) then ! no more heat source and sink BinQsupA(TempQuasi_Bin_Upper) = -(T_room - Bin_TemperAP(TempQuasi_Bin_Upper))/ExternalResAP if (ExternalResOP.ne.0) then BinQsupO(TempQuasi_Bin_Upper) = -(T_room - Bin_TemperOP(TempQuasi_Bin_Upper))/ExternalResOP else BinQsupO(TempQuasi_Bin_Upper) = 0 endif else BinQsupA(TempQuasi_Bin_Upper) = BinQsupA(TempQuasi_Bin_Upper-1) BinQsupO(TempQuasi_Bin_Upper) = BinQsupO(TempQuasi_Bin_Upper-1) endif do nBin = (TempQuasi_Bin_Lower), (TempQuasi_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 = TempQuasi_Bin_Lower,TempQuasi_Bin_Upper !col_len = RCutoff_SiIn BinQsupTotal(nBin) = BinQsupA(nBin) + BinQsupO(nBin) enddo do nBin = TempQuasi_Bin_Lower,TempQuasi_Bin_Upper !col_len = RCutoff_SiIn TempAvg_BinQsupA(nBin) = TempAvg_BinQsupA(nBin) + BinQsupA(nBin) TempAvg_BinQsupO(nBin) = TempAvg_BinQsupO(nBin) + BinQsupO(nBin) TempAvg_BinQsupTotal(nBin) = TempAvg_BinQsupTotal(nBin) + BinQsupTotal(nBin) TempAvg_BinAPEne(nBin) = TempAvg_BinAPEne(nBin) + BinAPEne(nBin) TempAvg_BinOPEne(nBin) = TempAvg_BinOPEne(nBin) + BinOPEne(nBin) TempAvg_Bin_OccupOP(nBin) = TempAvg_Bin_OccupOP(nBin) + Bin_OccupOP(nBin) TempAvg_Bin_TemperAP(nBin) = TempAvg_Bin_TemperAP(nBin) + Bin_TemperAP(nBin) TempAvg_Bin_TemperOP(nBin) = TempAvg_Bin_TemperOP(nBin) + Bin_TemperOP(nBin) TempAvg_Quasi_peOWm2(nBin) = TempAvg_Quasi_peOWm2(nBin) + Quasi_peOWm2(nBin) TempAvg_Quasi_peTOWm2(nBin) = TempAvg_Quasi_peTOWm2(nBin) + Quasi_peTOWm2(nBin) TempAvg_Quasi_peLOWm2(nBin) = TempAvg_Quasi_peLOWm2(nBin) + Quasi_peLOWm2(nBin) TempAvg_Quasi_ppUpWm2(nBin) = TempAvg_Quasi_ppUpWm2(nBin) + Quasi_ppUpWm2(nBin) enddo enddo ! for average do nBin = TempQuasi_Bin_Lower,TempQuasi_Bin_Upper !col_len = RCutoff_SiIn BinQsupA(nBin) = TempAvg_BinQsupA(nBin)/real(step,kind = 8) BinQsupO(nBin) = TempAvg_BinQsupO(nBin)/real(step,kind = 8) BinQsupTotal(nBin) = TempAvg_BinQsupTotal(nBin)/real(step,kind = 8) BinAPEne(nBin) = TempAvg_BinAPEne(nBin)/real(step,kind = 8) BinOPEne(nBin) = TempAvg_BinOPEne(nBin)/real(step,kind = 8) Bin_OccupOP(nBin) = TempAvg_Bin_OccupOP(nBin)/real(step,kind = 8) Bin_TemperAP(nBin) = TempAvg_Bin_TemperAP(nBin)/real(step,kind = 8) Bin_TemperOP(nBin) = TempAvg_Bin_TemperOP(nBin)/real(step,kind = 8) Quasi_peOWm2(nBin) = TempAvg_Quasi_peOWm2(nBin)/real(step,kind = 8) Quasi_peTOWm2(nBin) = TempAvg_Quasi_peTOWm2(nBin)/real(step,kind = 8) Quasi_peLOWm2(nBin) = TempAvg_Quasi_peLOWm2(nBin)/real(step,kind = 8) Quasi_ppUpWm2(nBin) = TempAvg_Quasi_ppUpWm2(nBin)/real(step,kind = 8) enddo 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(TempQuasi_Bin_Lower) + qsupplyA*BinSize*(TempQuasi_Bin_Lower-nBin)/k_pA Bin_TemperOP(nBin) = Bin_TemperOP(TempQuasi_Bin_Lower) + qsupplyO*BinSize*(TempQuasi_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 Quasi_peOWm2(nBin) = 0 Quasi_peTOWm2(nBin) = 0 Quasi_peLOWm2(nBin) = 0 Quasi_ppUpWm2(nBin) = 0 enddo do nBin = TempQuasi_Bin_Upper+1,num_Bins !col_len = RCutoff_SiIn BinQsupA(nBin) = 0 BinQsupO(nBin) = 0 BinQsupTotal(nBin) = 0 BinAPEne(nBin) = BinAPEne(TempQuasi_Bin_Upper) BinOPEne(nBin) = BinOPEne(TempQuasi_Bin_Upper) Bin_OccupOP(nBin) = Bin_OccupOP(TempQuasi_Bin_Upper) Bin_TemperAP(nBin) = Bin_TemperAP(TempQuasi_Bin_Upper) Bin_TemperOP(nBin) = Bin_TemperOP(TempQuasi_Bin_Upper) Quasi_peOWm2(nBin) = 0 Quasi_peTOWm2(nBin) = 0 Quasi_peLOWm2(nBin) = 0 Quasi_ppUpWm2(nBin) = 0 enddo endif do nBin = 1,num_BinsComb ! PhononSim = 1 Hot Temp constant/ 2 Cold Temp Constant call TempVar_sc_table(nBin) ! constant temperature outside them enddo open (unit = 62, file = 'QuasiSteadyResult.txt', status = 'replace', action = 'write' ) write(62,*)'QuasiSteadyNumericalCalculation ','QuasiSteadySetting.txt' write(62,'(a,a)'), 'Setting: ', 'QuasiBoundConsid BoundQSetting QuasiTimeStep QuasiTimeMax QuasiTime MaxTempChange ExternalResOP ExternalResAP' write(62,'(2x,2i8,10e20.9)'),QuasiBoundConsid,BoundQSetting,QuasiTimeStep,QuasiTimeMax,QuasiTime,MaxTempChange,ExternalResOP,ExternalResAP write(62,'(a)'), '----------------------------------------------------------------------------------------------------------------------------------------' write(62,'(a,a)'), 'Pos ', ' BinQsupA BinQsupO BinQsupTotal Bin_TemperAP Bin_TemperOP Quasi_peOWm2 Quasi_peTOWm2 Quasi_peLOWm2 Quasi_ppUpWm2' write(62,'(a)'), '----------------------------------------------------------------------------------------------------------------------------------------' do nBin = 1, num_Bins write(62,'(1x,e20.9,10e20.9)'),loc_bin(nBin),BinQsupA(nBin),BinQsupO(nBin),BinQsupTotal(nBin),Bin_TemperAP(nBin),Bin_TemperOP(nBin),Quasi_peOWm2(nBin),Quasi_peTOWm2(nBin),Quasi_peLOWm2(nBin),Quasi_ppUpWm2(nBin) enddo close(62) close(63) close(64) return end subroutine InitQuasiSteadyTdist() ! use particles use scattering use simulation implicit none integer :: nBin,TempNumBins, TempQuasi_Bin_Upper, TempQuasi_Bin_Lower real ( kind = 8 ) :: UnitVolAPEne,avePPupWm2 UnitVolAPEne = 0.0d0 do nBin = 1, num_Bins Quasi_peLOWm2(nBin) = PreLoop_Bin_EpLOWm2(nBin) Quasi_peTOWm2(nBin) = PreLoop_Bin_EpTOWm2(nBin) Quasi_peOWm2(nBin) = Quasi_peLOWm2(nBin)+Quasi_peTOWm2(nBin) enddo ! parameter ! call InitializeAPTable() ! energy for AP etc if (LOmodeOccPhonSim.eq.1) then ! Included in initialize If we choose LOmode1.. kp should be smaller ConstForEneOp = RepresentEpLO*e_c*BinSize/Default_primV ! J/m^2 else ConstForEneOp = RepresentEpLO*e_c*BinSize*3.0/Default_primV ! J/m^2 endif if (GivenQsupply.ne.1) then ! 1 given q 2: updated q if (QuasiBoundConsid.eq.1)then ! total energy Hot Phonon factor qsupply = PreLoop_TP_EneWm2_eb else qsupply = PreLoop_TP_EneWm2 endif endif qsupplyO = HotPhononFactor*qsupply*k_pO/(k_pA+k_pO) ! A/k_GaAs of heat flows through qsupplyA = qsupply - qsupplyO if (QuasiBoundConsid.eq.1)then ! total energy Hot Phonon factor TempQuasi_Bin_Lower = num_boundBin+1 TempQuasi_Bin_Upper = num_Bins-num_boundBin else TempQuasi_Bin_Lower = 1 TempQuasi_Bin_Upper = num_Bins endif ! avePPupWm2 = ((HotPhononFactor-1)/HotPhononFactor)*qsupplyO/real((TempQuasi_Bin_Upper-TempQuasi_Bin_Lower+1),kind = 8) avePPupWm2 = qsupplyA/real((TempQuasi_Bin_Upper-TempQuasi_Bin_Lower+1),kind = 8) do nBin = TempQuasi_Bin_Lower,TempQuasi_Bin_Upper ! PhononSim = 1 Hot Temp constant/ 2 Cold Temp Constant if (nBin.eq.TempQuasi_Bin_Lower) then BinQsupO(nBin) = qsupplyO - PreLoop_Bin_EpLOWm2(nBin) + avePPupWm2 BinQsupA(nBin) = qsupplyA - PreLoop_Bin_EpTOWm2(nBin) - avePPupWm2 else BinQsupO(nBin) = BinQsupO(nBin-1) - PreLoop_Bin_EpLOWm2(nBin)+ avePPupWm2 BinQsupA(nBin) = BinQsupA(nBin-1) - PreLoop_Bin_EpTOWm2(nBin)- avePPupWm2 endif enddo TempNumBins = TempQuasi_Bin_Upper-TempQuasi_Bin_Lower + 1 do nBin = 1,TempNumBins ! BinSize if (nBin.eq.1) then Bin_TemperAP(TempNumBins+TempQuasi_Bin_Lower-nBin) = T_room + ExternalResAP*BinQsupA(TempNumBins+TempQuasi_Bin_Lower-nBin) if (k_pO.gt.0) then Bin_TemperOP(TempNumBins+TempQuasi_Bin_Lower-nBin) = T_room + ExternalResOP*BinQsupO(TempNumBins+TempQuasi_Bin_Lower-nBin) else Bin_TemperOP(TempNumBins+TempQuasi_Bin_Lower-nBin) = Bin_TemperAP(TempNumBins+TempQuasi_Bin_Lower-nBin) endif else Bin_TemperAP(TempNumBins+TempQuasi_Bin_Lower-nBin) = Bin_TemperAP(TempNumBins+TempQuasi_Bin_Lower-nBin+1) + BinSize*BinQsupA(TempNumBins+TempQuasi_Bin_Lower-nBin)/k_pA if (k_pO.gt.0) then Bin_TemperOP(TempNumBins+TempQuasi_Bin_Lower-nBin) = Bin_TemperOP(TempNumBins+TempQuasi_Bin_Lower-nBin+1) + BinSize*BinQsupO(TempNumBins+TempQuasi_Bin_Lower-nBin)/k_pO else Bin_TemperOP(TempNumBins+TempQuasi_Bin_Lower-nBin) = Bin_TemperOP(TempNumBins+TempQuasi_Bin_Lower-nBin+1) endif endif enddo do nBin = TempQuasi_Bin_Lower,TempQuasi_Bin_Upper ! PhononSim = 1 Hot Temp constant/ 2 Cold Temp Constant Bin_OccupOP(nBin) = 1/(exp((RepresentEpLO*e_c)/(kB*Bin_TemperOP(nBin))) - 1) BinQsupTotal(nBin) = BinQsupO(nBin) + BinQsupA(nBin) 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 if (QuasiBoundConsid.eq.1)then do nBin = 1,num_boundBin !col_len = RCutoff_SiIn BinQsupA(nBin) = 0 BinQsupO(nBin) = 0 BinQsupTotal(nBin) = 0 BinAPEne(nBin) = BinAPEne(TempQuasi_Bin_Lower) BinOPEne(nBin) = BinOPEne(TempQuasi_Bin_Lower) Bin_OccupOP(nBin) = Bin_OccupOP(TempQuasi_Bin_Lower) Bin_TemperAP(nBin) = Bin_TemperAP(TempQuasi_Bin_Lower) Bin_TemperOP(nBin) = Bin_TemperOP(TempQuasi_Bin_Lower) Quasi_peOWm2(nBin) = 0 Quasi_peTOWm2(nBin) = 0 Quasi_peLOWm2(nBin) = 0 Quasi_ppUpWm2(nBin) = 0 enddo do nBin = TempQuasi_Bin_Upper+1,num_Bins !col_len = RCutoff_SiIn BinQsupA(nBin) = 0 BinQsupO(nBin) = 0 BinQsupTotal(nBin) = 0 BinAPEne(nBin) = BinAPEne(TempQuasi_Bin_Upper) BinOPEne(nBin) = BinOPEne(TempQuasi_Bin_Upper) Bin_OccupOP(nBin) = Bin_OccupOP(TempQuasi_Bin_Upper) Bin_TemperAP(nBin) = Bin_TemperAP(TempQuasi_Bin_Upper) Bin_TemperOP(nBin) = Bin_TemperOP(TempQuasi_Bin_Upper) Quasi_peOWm2(nBin) = 0 Quasi_peTOWm2(nBin) = 0 Quasi_peLOWm2(nBin) = 0 Quasi_ppUpWm2(nBin) = 0 enddo endif ! ExternalResAP = (Bin_TemperAP(num_Bins) - T_room)/qsupplyA ! if (k_pO.ne.0) then ExternalResOP = ExternalResAP*(k_pA/k_pO) else ExternalResOP = 0 endif return end