24小时热门版块排行榜

返回列表

当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖

dailjiao0706

新虫 (初入文坛)

应助: 0 (幼儿园)
金币: 4
帖子: 10
在线: 4.9小时
虫号: 2059068
注册: 2012-10-13
专业: 原子和分子物理

[求助] fortran中的数组越界问题，

MODULE  my_Constants
!------------------------------------------------------------------------------------
IMPLICIT NONE

integer, parameter :: r8=selected_real_kind(p=13) !double precision
real(r8), PARAMETER :: pi = 3.1415926535897932 !pi
real(r8), PARAMETER :: c = 2.99792458e08       !Speed of light (m/s)
real(r8), PARAMETER :: e0=1.0e07/(4*pi*c**2)    !vacuum permittivity
real(r8), PARAMETER :: er=1             !relative permittivity
real(r8), PARAMETER :: u0=pi*4e-07       !vacuum permeability
real(r8), PARAMETER :: mp=1.6726231e-27          !Proton rest mass(kg)
real(r8), PARAMETER :: me=9.1093897e-31          !Electron rest mass(kg)
real(r8), PARAMETER :: e=1.60217733e-19          !Elementary charge (C)
END MODULE  my_Constants

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

program CE
!------------------------------------------------------
use my_Constants

implicit none
!------------- parameters -----------------------------------------
real(r8), PARAMETER :: e2_4pie0=e*e/(4*pi*e0)
real(r8), PARAMETER :: Iinner=4.3247737e19       !intensity for full inner ionization(W/m2)
real(r8), PARAMETER :: lambda=8e-7          !electric field's wavelength(m)
real(r8), PARAMETER :: we=2*pi*c/lambda    !electric field's frequency
real(r8), PARAMETER :: th=50e-15          !FWHM: tp=th/(2*sqrt(log(2.0)))
real(r8), PARAMETER :: tstart=0.0
INTEGER, PARAMETER :: M=5005          !M----t
INTEGER, PARAMETER :: iduation=5004
real(r8), dimension(M) :: RH_bord,RI_bord,Ne,Inten,pou_e,ttime
real(r8), PARAMETER :: IeV=1.6e19          !intensity for ionic energy ~300eV(W/m2)
real(r8), PARAMETER :: er0=1e-4
real(r8), PARAMETER :: KeV=1000*e    ! change energy unit from J to KeV
real(r8), PARAMETER :: dt=10.0e-17
real(r8), PARAMETER :: dt2=dt*dt
real(r8), PARAMETER :: dr=5.0e-11
real(r8), PARAMETER :: F0=2.0e21             !peak intensity of electric field(W/cm2)
INTEGER :: i,j,b,y,y0,y1,i_inner,yHmax,yImax,y2,dy2,jmax
real(r8) :: t,tp,Nemin,RH_y,RI_y,EH,EI
real(r8) :: yminH,yminI,RIHmin,RHImin,yHbord,yIbord
real(r8) :: EkHmax,EkImax,EaverH,EaverI,dNH,dNI,EH_term,EI_term,sumH,sumI
!-----------the following parameters relates to the size distribution ---------------------
real(r8), PARAMETER :: R0=5.0e-9    !4e-9 for 2w, 2e-9 for 3w
real(r8), PARAMETER :: n0=1.6e28             !cluster density  (/m**3)
real(r8), PARAMETER :: NI=(4.0/3)*pi*(R0**3)*n0
real(r8), PARAMETER :: NH=NI*4
real(r8), PARAMETER :: pou0_I=n0
real(r8), PARAMETER :: pou0_H=n0*4
INTEGER, PARAMETER :: QI=4
real(r8), PARAMETER :: mI=12*mp
!----------------------------------------------
real(r8), PARAMETER :: dr_alongr=5.0e-11
integer :: Ndr       ! Ndr=int(R0(j)/dr_alongr)+1
INTEGER :: ierr
real(r8),dimension(

, allocatable :: RH_0,RI_0
real(r8),dimension(:,

, allocatable :: RH,RI,pout_H,pout_I,NHH_y,NIH_y,NHI_y,NII_y,F_HH,F_HI,F_IH,F_II,F_He,F_Ie ! 原程序为RH,RI
write(*,*) "i_inner=",i_inner  ! 输出i_inner 值为0，说明在默认情况下不赋值则系统给整数参数赋值为0，这就导致i_inner-1=-1,显然数组的下标不能为负值。再有，新版fortran数组下标从1开始,所以给i_inner赋初值为2
i_inner=2
!--------------------------------------------------------
write(*,*) "i_inner=",i_inner
!------------------------------------------------------------------
open(unit=9,file='F_t.dat')  !
open(unit=15,file='pout_H.dat')  !
open(unit=16,file='RH_r_t_R0_5.dat')  !
open(unit=17,file='RI_r_t_R0_5.dat')  !
open(unit=18,file='pout_I.dat')  !
open(unit=19,file='NII.dat') !
open(unit=20,file='NIH.dat')  !
open(unit=21,file='NHH.dat') !
open(unit=22,file='NHI.dat')  !


!************************************************************************************

!========== loop for the cluster-size distribution ====================================
Ndr=int(R0/dr_alongr)+1
write(*,*) "Ndr=",Ndr
!======= Allocate arrays according to the exact different cluster-size ===============
allocate (RH_0(Ndr),stat=ierr)
if(ierr/=0)then
write(*,*)"Unable to allocate memory for array RH_0 "
stop "Stopping"
end if
!
allocate (RI_0(Ndr),stat=ierr)
if(ierr/=0)then
write(*,*)"Unable to allocate memory for array RI_0 "
stop "Stopping"
end if
!
allocate (RH(Ndr,M),stat=ierr)
if(ierr/=0)then
write(*,*)"Unable to allocate memory for array RH "
stop "Stopping"
end if
!
allocate (RI(Ndr,M),stat=ierr)
if(ierr/=0)then
write(*,*)"Unable to allocate memory for array RI "
stop "Stopping"
end if
!RH_0=0
!pause

!======= divide many layers along the radius ========================================
layer1_loop: do y=1,Ndr
RH_0(y)=(y-1)*dr_alongr
RI_0(y)=(y-1)*dr_alongr
end do layer1_loop
!======== initialize the cluster radius & charge =========================
t=tstart
layer2_loop: do i=1,iduation-1
   ttime(i)=t
         RH_bord(i)=R0
         RI_bord(i)=R0
   Ne(i)=0.0
layer3_loop: do y=1,Ndr
   RH(y,i)=RH_0(y)
   RI(y,i)=RI_0(y)
end do layer3_loop
t=t+dt
end do layer2_loop
!------ set the connecting condition between inner and outer ionization ----------------
      Ne(i_inner-1)=0.0
      Ne(i_inner)=0.0
      Ne(i_inner)=0.0
pou_e(i_inner)=0.0
!==== compute temporal evolution of layers along  radius during outer ionization period ===
!
outer1_loop: do i=2,iduation
!------ judge if the cluster is outer-ionized fully ----------
full_check1: if (Ne(i)<=0.0) then
            Ne(i)=0.0
!    pou_e(i)=0.0
end if full_check1
!---- Ne(t) at the next time point must be smaller than that of the present time ----
reasonable_check: if (Ne(i)>Ne(i-1)) then
   Ne(i)=0.0
end if reasonable_check
   pou_e(i)=0.0
!==== compute temporal evolution of layers along  radius==========
layer4_loop: do y0=2,Ndr
RH_bord(i+1)=RH_bord(i)
RI_bord(i+1)=RI_bord(i)
      !------ set the RH_y/RC_y as the layer subjeIt -----------
RH_y=RH(y0,i)
RI_y=RI(y0,i)
!----- assume the electrons distribute uniformly within the cluster -----------
F_He=0.0
F_Ie=0.0
!------------------------
NHH_y=0.0
NIH_y=0.0
NHI_y=0.0
NII_y=0.0
RIHmin=0.0
RHImin=0.0
!==== count the number of C/H ions included in a certain radius RH_y/RI_y ==========
layer5_loop: do y1=2,Ndr
!------------------------
NHH_count: if (RH(y1,i)<=RH_y) then
NHH_y(y1,1)=NHH_y(y1-1,1)+pou0_H*4*pi/3*          (RH(y1,1)**3-RH(y1-1,1)**3)
end if NHH_count
!------------------------
NII_count: if (RI(y1,i)<=RI_y) then
NII_y(y1,1)=NII_y(y1-1,1)+pou0_I*4*pi/3*(RI(y1,1)**3-RI(y1-1,1)**3)
end if NII_count
!------------------------
NIH_count: if (RH(y1,i)<=RI_y) then
NIH_y(y1,1)=NIH_y(y1-1,1)+pou0_H*4*pi/3*(RH(y1,1)**3-RH(y1-1,1)**3)
if (RH(y1,i)>RIHmin) then
RIHmin=RH(y1,i)
yminH=y1
end if
end if NIH_count
!------------------------
NHI_count: if (RI(y1,i)<=RH_y) then
NHI_y(y1,1)=NHI_y(y1-1,1)+pou0_I*4*pi/3*(RI(y1,1)**3-RI(y1-1,1)**3)
if (RI(y1,i)>RHImin) then
RHImin=RI(y1,i)
yminI=y1
end if
end if NHI_count
!------------------------
end do  layer5_loop
!===== make some repairs for perfect =============================================
if (NIH_y(y0,i)==0.0) then
pout_H(y0,i)=pou0_H*(RH(2,1)/RH(y0,i))**3
NIH_y(y0,i)=4*pi/3*RI_y**3*pout_H(2,i)
      else if (RI_y<=RH_bord(i)) then
! pout_H(y0,i)=NH/(4*pi/3*RH_bord(i)**3)    !--little effected--!
pout_H(y0,i)=pou0_H*(RH(yminH,1)**3-RH(yminH-1,1)**3)/abs(RH(yminH,i)**3-RH(yminH-1,i)**3)
NIH_y(y0,i)=NIH_y(y0-1,i)+pout_H(y0,i)*4*pi/3*(RI_y**3-RIHmin**3)
end if
!---------------------------------------------
if (NHI_y(y0,i)==0.0) then
pout_I(y0,i)=pou0_I*(RI(2,1)/RI(y0,i))**3                            !
NHI_y(y0,i)=4*pi/3*RH_y**3*pout_I(y0,i)
else if(RH_y<=RI_bord(i)) then
! pout_I=NI/(4*pi/3*RI_bord(i)**3)    !--little effected--!
pout_I(y0,i)=pou0_I*(RI(yminI,1)**3-RI(yminI-1,1)**3)/abs(RI(yminI,i)**3-RI(yminI-1,i)**3)
NHI_y(y0,i)=NHI_y(y0-1,i)+pout_I(y0,i)*4*pi/3*(RH_y**3-RHImin**3)
end if
!------------------------
F_HH(y0,i)=e2_4pie0*NHH_y(y0,i)/(RH_y**2)
F_HI(y0,i)=e2_4pie0*QI*NHI_y(y0,i)/(RH_y**2)
F_IH(y0,i)=QI*e2_4pie0*NIH_y(y0,i)/(RI_y**2)
F_II(y0,i)=QI*QI*e2_4pie0*NII_y(y0,i)/(RI_y**2)
!------------------------------------------------------------------------
RH(y0,i+1)=dt2/mp*(F_HH(y0,i+1)+F_HI(y0,i+1)+F_He(y0,i+1))+2*RH(y0,i)-RH(y0,i-1)
RI(y0,i+1)=dt2/mI*(F_IH(y0,i+1)+F_II(y0,i+1)+F_Ie(y0,i+1))+2*RI(y0,i)-RI(y0,i-1)
!---------- find the new maximum borders --------------------------------
RH_border: if (RH(y0,i+1)>RH_bord(i+1)) then
yHbord=y0
RH_bord(i+1)=RH(y0,i+1)
end if RH_border
RI_border: if (RI(y0,i+1)>RI_bord(i+1)) then
yIbord=y0
RI_bord(i+1)=RI(y0,i+1)
end if RI_border
!------------------------------------------------------------------------
end do layer4_loop
!--------- set Ne(t) at the next time point ---------------------------------
  full_check2: if (Ne(i)<=0.0) then       ! if Ne(i)<=0  "it's full outer ionization."
Ne(i+1)=0.0
! pou_e(i+1)=0.0
  else full_check2
         Ne(i+1)=0.0
! pou_e(i+1)=Ne(i+1)/(4.0*pi*RH_bord(i+1)**3/3)
  end if full_check2
!--------------------------------------------------------------------
end do outer1_loop
!******************************************************************************************
!======== output the necessary data =====================================
write(16,*) R0
      write(16,*) Ndr
out_time: do i=3,iduation,50
write(16,*) ttime(i)/1e-15
write(17,*) ttime(i)/1e-15
out_layer: do y=1,Ndr
write(16,*) RH(y,i)
write(16,*) RH(y,i-1)
write(15,*) pout_H(y,i)
      write(18,*) pout_I(y,i)
      write(19,*) NII_y(y,i)
      write(20,*) NIH_y(y,i)
      write(21,*) NHH_y(y,i)
      write(22,*) NHI_y(y,i)
end do out_layer
end do out_time
!======It is important to deallocate the arrays ====================================
deallocate(RH_0,RI_0,RH,RI)
!*********************************************************************
close(15)
close(16)
close(17)
!----------------------------------------------------------------------------
end program CE
!================================================

回复此楼