应《网络安全法》要求，自2017年10月1日起，未进行实名认证将不得使用互联网跟帖服务。为保障您的帐号能够正常使用，请尽快对帐号进行手机号验证，感谢您的理解与支持！

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 计算模拟区 » 程序语言 » Fortran » fortran中的数组越界问题，
北京石油化工学院2026年研究生招生接收调剂公告
41/1返回列表
查看: 3076  |  回复: 3
只看楼主 @他人 存档 新回复提醒 (忽略) 收藏    在APP中查看

dailjiao0706

新虫 (初入文坛)

[求助] 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
!================================================
回复此楼

» 猜你喜欢

» 本主题相关价值贴推荐，对您同样有帮助:
1楼 2013-10-16 20:39:09
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
回帖置顶 ( 共有1个 )

dailjiao0706

新虫 (初入文坛)
dailjiao0706: 回帖置顶 2013-10-17 13:54:48
引用回帖:
2楼: Originally posted by virtualzx at 2013-10-16 21:56:02
除非你手动指定，fortran不要求编译器如何初始化变量，没有默认情况这说，任何一个编译器都可以随意指定。很多都是不处理，初值是随机的，之前那块内存是什么值就出来什么值。第一次用i_inner之前要先手动赋值。

大侠，我覆初值了   i_inner=2，write(*,*) "i_inner=",i_inner，求大侠帮忙
赞一下 回复此楼
3楼2013-10-17 13:45:26
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
普通回帖

virtualzx

木虫 (著名写手)

【答案】应助回帖

感谢参与，应助指数 +1
除非你手动指定，fortran不要求编译器如何初始化变量，没有默认情况这说，任何一个编译器都可以随意指定。很多都是不处理，初值是随机的，之前那块内存是什么值就出来什么值。第一次用i_inner之前要先手动赋值。
赞一下 回复此楼
2楼2013-10-16 21:56:02
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

virtualzx

木虫 (著名写手)

【答案】应助回帖

内容已删除
赞一下 回复此楼
4楼2013-10-17 23:58:35
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 dailjiao0706 的主题更新
41/1返回列表
普通表情 龙 兔 虎 猫
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[考研] 求调剂 一志愿西南交通大学085701环境工程 282分 +7 多多爱吃汉堡 2026-04-04 7/350 2026-04-05 19:47 by 福农资环_环境基
[考研] 本科211，293分请求调剂 +8 莲菜就是藕吧 2026-04-03 9/450 2026-04-05 19:12 by 蓝云思雨
[考研] 化学0703-一志愿211-338分求调剂 +7 vants 2026-04-05 7/350 2026-04-05 18:17 by cql1109
[考研] 285求调剂 +11 哦呦呼o 2026-04-04 11/550 2026-04-05 17:59 by 猪会飞
[考研] 一志愿上海海洋大学083200食品学硕，求调剂，接受其他专业083200 +4 what张 2026-04-04 5/250 2026-04-05 14:07 by chw1980_0
[考研] 081200-11408-276学硕求调剂 +4 崔wj 2026-04-04 5/250 2026-04-05 14:06 by imissbao
[考研] 283分求调剂 +9 试试看呗 2026-04-04 9/450 2026-04-05 10:27 by 果冻大王
[考研] 271分求调剂学校 +12 zph158488！ 2026-04-02 13/650 2026-04-05 10:13 by lqwchd
[考研] 调剂求助 +10 想换手机不想解� 2026-04-02 13/650 2026-04-05 09:41 by sam3303
[考研] 求调剂 +7 xzghyuj 2026-04-04 7/350 2026-04-04 22:25 by oooqiao
[考研] 085400电子信息319求调剂（接受跨专业调剂） +5 星星不眨眼喽 2026-04-03 6/300 2026-04-04 21:50 by hemengdong
[考研] 11408 一志愿西电，277分求调剂 +4 zhouzhen654 2026-04-03 4/200 2026-04-04 18:10 by 猪会飞
[考研] 297求调剂 +11 ljy20040718！ 2026-04-03 13/650 2026-04-04 09:23 by 来看流星雨10
[考研] 265求调剂 +20 梁梁校校 2026-04-01 21/1050 2026-04-04 00:38 by userper
[考研] 350一志愿北京航空航天大学08500材料科学与工程求调剂 +5 kjnasfss 2026-04-03 5/250 2026-04-03 22:29 by 无际的草原
[考研] 数一英一285求调剂 +7 AZMK 2026-04-03 9/450 2026-04-03 13:03 by ms629
[考研] 农学考研求调剂 +3 dkdkxm 2026-04-01 3/150 2026-04-02 16:04 by wangjagri
[考博] 材料工程专业硕士申博 +3 麟正宇 2026-03-30 3/150 2026-04-02 15:04 by greychen00
[考研] 379求调剂 +3 ?苦瓜不苦 2026-04-01 3/150 2026-04-01 20:09 by 暮云清寒
[考研] 物理学调剂 +4 小羊36 2026-03-30 4/200 2026-03-31 16:16 by lishahe
信息提示
关闭
请填处理意见
关闭