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

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 计算模拟区 » 程序语言 » 【求助】帮忙看下这个程序，实在调不出来啊
北京石油化工学院2026年研究生招生接收调剂公告
121/2返回列表12下一页
查看: 761  |  回复: 11
只看楼主 @他人 存档 新回复提醒 (忽略) 收藏    在APP中查看
当前主题已经存档。

yhwsmile

金虫 (小有名气)

[交流] 【求助】帮忙看下这个程序，实在调不出来啊

各位高手帮帮忙啊，程序一运行就会显示
run-time error M6201: math
- sqrt:domain error
这个应该是sqrt()中出现负数了吧，但我验证了我的程序不可能出现负数啊，而且我把sqrt语句注释了以后，运行还是出现上述错误提示，这是为什么啊?我实在调不出来了，程序是我们导师以前编的，我上传到附件里了，哪位大虾帮忙看下啊，谢谢了啊

由于程序中要读取的txt文件超出了附件上传限制我放到邮箱里了www.126.com
126：yhfortran
密码：123456

程序如下：      program ave !写均值及方差协方差文件程序
       parameter(n=1,n0=18000,n1=7,pi=3.14159265)
c          n为资料天数，n0为半小时资料记录数，n1为所使用的超声变量数,m为旋转矩阵的维数
c         integer i0(n1)
         real new_data(n0,n1)
         real av(n1),du(3),uv(12),ds(n1),star(3),uuvc(12),roll(2)
c      av:各变量半小时均值；du:平均水平风速及风向角；uv:风速与各变量协方差；ds:各变量方差；star:特征风、温、湿
c     aav:三次旋转订正后的平均风速；uuvc:三次旋转后的协方差；roll(2):二次旋转的角度！第一次v=0；第二次旋转=0；
c     dds:三次旋转后风速的方差,
     
         character c*27,cc1*11
         Character(100) cc
         integer error
         data  ba/99999.0/
       open(1,file='filename1.txt')

       open(3,file='e:\ave.txt')
      
      write(3,102) 'data_time   U(m/s)   V(m/s)   W(m/s)   C(mg/m3)
     *   Q(g/m3)  T(c)   P(hpa)    den(kg/m3)   UU     VV      WW     CC
     *    QQ     TT    PP         UV          UW     UC      UQ     UT
     *    VW     VC    VQ         VT          WC     WQ      WT     err'   


       do ii=1,n  
         read(1,*,iostat=error) c
       
         open(111,file='E:\'//c//'.txt')
         do j=1,48  !这里j为一天24小时，共48个半小时循环
       call dele_data1(n0,n0,n1,c,new_data,av,du,pi,ba,cc1)
         call uuv(n0,n1,new_data,av,ba,den,ds,uv,err)

ccc
       write(3,1001) cc1,(av(i),i=1,7),den,(ds(i),i=1,7),(uv(i),i=1,12),
     * err
         enddo
         write(*,*) cc1
       close(111)
       enddo



1001  format(a11,28f9.3)  
  102  format(a240)
      

       close(1)
c         close(2)
       close(3)
       end       


******************子程序1:用于去除奇异点，并计算均值，返回均值和去除奇异点后的新序列*****************                
         subroutine dele_data1(l,n0,n1,c,new_data,av2,du,pi,ba,cc1)
         character c*27,data_time(n0)*11,cc*100,cc1*11
         integer i0(n1),error1
         real sonic_data(n0,7),new_data(n0,n1)
         real aa(n1),av1(n1),av2(n1),cgm(n1),du(3),ds(n1)

         do i=1,n1
         i0(i)=0
         aa(i)=0.0
         cgm(i)=0.0
         enddo
         do j=1,l
         read(111,*,IOSTAT=error1) data_time(j),(sonic_data(j,jj),jj=1,7)
         if(error1/=0)exit
         do jj=1,n1
              if(abs(sonic_data(j,jj)).lt.ba) then
       new_data(j,jj)=sonic_data(j,jj)
         aa(jj)=aa(jj)+sonic_data(j,jj)
         i0(jj)=i0(jj)+1
         else
         new_data(j,jj)=ba
         endif
       enddo
         enddo
         cc1=data_time(l)
       do i=1,n1
         av1(i)=aa(i)/i0(i)  
         enddo
         do ii=1,l
         do i=1,n1
         if(abs(sonic_data(ii,i)).lt.ba) then
         cgm(i)=cgm(i)+((sonic_data(ii,i)-av1(i))**2/i0(i))
         endif
         enddo
         enddo
         do i=1,n1
         cgm(i)=sqrt(cgm(i))
         enddo
         do i=1,n1
         i0(i)=0
         aa(i)=0.0
         enddo
         do ii=1,l
         do i=1,n1
         if(new_data(ii,i).ne.ba.and.abs(new_data(ii,i)-av1(i))
     * .le.cgm(i)*4) then
         aa(i)=aa(i)+new_data(ii,i)
         i0(i)=i0(i)+1
         else
         new_data(ii,i)=ba
         endif
         enddo
         enddo
         do i=1,n1
         av2(i)=aa(i)/i0(i)
         enddo
         du(1)=SQRT(av2(1)*av2(1)+av2(2)*av2(2))
         du(2)=atan(abs(av2(1))/abs(av2(2)))      
ccc   处理合成风向

      if(av2(1).ge.0.)then
              if(av2(2).gt.0.) then
            du(2)=du(2)+pi
            else
            du(2)=2*pi-du(2)
            endif
        elseif(av2(1).lt.0.)then
          if(av2(2).gt.0.) then
           du(2)=pi-du(2)
          else
           du(2)=du(2)
          endif
       endif
         du(3)=du(2)*180/pi         

         if(l.lt.n0) then
         do i=l+1,n0
         do j=1,n1
         sonic_data(i,j)=ba
         new_data(i,j)=ba
         enddo
         sonic_data(i,7)=ba
         enddo
         endif
       end

***************子程序2：计算超声变量协方差程序********************
       subroutine uuv(n0,n1,new_data,av,ba,den,ds,uv,err)
         real new_data(n0,n1),av(n1),dd(n0,n1),ds(n1),bs(n1),uv(12)
       integer i0(n1),ii0(12)

        den=(av(7)/10)/(0.287*(273.15+av(6)))     
****注意：这里气压的单位为kpa，密度单位为Kg/m3

        do i=1,n1
        ds(i)=0.0
        i0(i)=0
        enddo

      dO i=1,n0
      do j=1,n1
        if(abs(new_data(i,j)).ne.ba)then
        dd(i,j)=new_data(i,j)-av(j)
        ds(j)=ds(j)+dd(i,j)**2
        write(50,*) ds(j)
        i0(j)=i0(j)+1
        else
        dd(i,j)=ba
        end if
        end do
        enddo
      
        err=0.0
        do i=1,n1
        ds(i)=ds(i)/i0(i)
        bs(i)=SQRT(ds(i))      
      if(bs(i).gt.2.0.and.i.ne.4) err=999   
      if(bs(i).gt.30.0.and.i.eq.4) err=999
      enddo

                     
        do ii=1,12
        uv(ii)=0.0
        ii0(ii)=0
        enddo
      
        DO i=1,n0
        if(dd(i,1).ne.ba.and.dd(i,2).ne.ba)then
        uv(1)=uv(1)+dd(i,1)*dd(i,2)
        ii0(1)=ii0(1)+1
        end if
        if(dd(i,1).ne.ba.and.dd(i,3).ne.ba)then
        uv(2)=uv(2)+dd(i,1)*dd(i,3)
        ii0(2)=ii0(2)+1
        end if
        if(dd(i,1).ne.ba.and.dd(i,4).ne.ba)then
        uv(3)=uv(3)+dd(i,1)*dd(i,4)
        ii0(3)=ii0(3)+1
        end if
        if(dd(i,1).ne.ba.and.dd(i,5).ne.ba)then
        uv(4)=uv(4)+dd(i,1)*dd(i,5)
        ii0(4)=ii0(4)+1
        end if
        if(dd(i,1).ne.ba.and.dd(i,6).ne.ba)then
        uv(5)=uv(5)+dd(i,1)*dd(i,6)
        ii0(5)=ii0(5)+1
        end if
      if(dd(i,2).ne.ba.and.dd(i,3).ne.ba)then
        uv(6)=uv(6)+dd(i,2)*dd(i,3)
        ii0(6)=ii0(6)+1
        end if
      if(dd(i,2).ne.ba.and.dd(i,4).ne.ba)then
        uv(7)=uv(7)+dd(i,2)*dd(i,4)
        ii0(7)=ii0(7)+1
        end if
      if(dd(i,2).ne.ba.and.dd(i,5).ne.ba)then
        uv(8)=uv(8)+dd(i,2)*dd(i,5)
        ii0(8)=ii0(8)+1
        end if
      if(dd(i,2).ne.ba.and.dd(i,6).ne.ba)then
        uv(9)=uv(9)+dd(i,2)*dd(i,6)
        ii0(9)=ii0(9)+1
        end if
      if(dd(i,3).ne.ba.and.dd(i,4).ne.ba)then
        uv(10)=uv(10)+dd(i,3)*dd(i,4)
        ii0(10)=ii0(10)+1
        end if
      if(dd(i,3).ne.ba.and.dd(i,5).ne.ba)then
        uv(11)=uv(11)+dd(i,3)*dd(i,5)
        ii0(11)=ii0(11)+1
        end if
      if(dd(i,3).ne.ba.and.dd(i,6).ne.ba)then
        uv(12)=uv(12)+dd(i,3)*dd(i,6)
        ii0(12)=ii0(12)+1
        end if
        enddo

        do i=1,12
      if(ii0(i).ne.0)then
        uv(i)=uv(i)/ii0(i)
        else
        uv(i)=ba
        endif
        enddo

      end


****一个二维矩阵于一个一维矩阵的相乘求解********************
        subroutine matrixmul2(a,b,m,n,c)
        dimension a(m,n),b(n),c(m)
      do i=1,m
         c(i)=0.0
          do l=1,n
          c(i)=c(i)+a(i,l)*b(l)
          end do
         end do
      end

***两个矩阵相乘运算求解程序子程序*****************************
      subroutine matrixmul(a,b,m,k,n,c)
        dimension a(m,n),b(n,k),c(m,k)
      do i=1,m
         do j=1,k
         c(i,j)=0.0
          do l=1,n
          c(i,j)=c(i,j)+a(i,l)*b(l,j)
          end do
        end do
         end do
      end

[ Last edited by yhwsmile on 2009-4-20 at 15:17 ]
回复此楼

» 猜你喜欢
1楼 2009-04-20 10:33:50
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

fspdlh

金虫 (正式写手)
★
kuhailangyu(金币+1,VIP+0):替楼主感谢！ 4-20 13:41
我不懂fortran,不过我想楼主可以让程序每运行一步都显示一下根号下的值，总会找到问题所在的
赞一下(10人) 回复此楼
2楼2009-04-20 12:15:58
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

chuanqingfu

金虫 (正式写手)
我在学习F过程中。看不出来。就当来学习了……
赞一下 回复此楼
3楼2009-04-21 08:33:09
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

西药

铜虫 (正式写手)
好久没用过FORTRAN 了
先仔细看下
赞一下 回复此楼
4楼2009-04-21 10:04:42
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

老虎大王

木虫 (著名写手)
★ ★
kuhailangyu(金币+2,VIP+0):替楼主感谢！ 4-22 16:05
你有这个子程序       
subroutine dele_data1(l,n0,n1,c,new_data,av2,du,pi,ba,cc1)
其中，有
                 do i=1,n1
         if(abs(sonic_data(ii,i)).lt.ba) then
         cgm(i)=cgm(i)+((sonic_data(ii,i)-av1(i))**2/i0(i))
         endif
         enddo

我初步监测到，当主程序第三次调用这个子程序的时候（即主程序do j=1,48，其中当j=3的时候），执行至此句之前，你的数组i0中各元素都是零。这样，上式中被零除，cgm(i)为NaN，这就造成后边的         cgm(i)=sqrt(cgm(i))  语句中，sqrt函数定义域错误。
同时，至少还导致了后边的：
         do i=1,n1
         av2(i)=aa(i)/i0(i)
         enddo
         du(1)=SQRT(av2(1)*av2(1)+av2(2)*av2(2))
                  du(2)=atan(abs(av2(1))/abs(av2(2)))  
也就有了同样的问题（av2为NaN，SQRT出错，atan也出错）......


请你检查你的算法中，数组i0的赋值方式。

Ps:  这只是初步发现的一个错误，我没有继续调试。先把它纠正了再说。
赞一下(10人) 回复此楼
5楼2009-04-21 13:27:28
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

snoopyzhao

至尊木虫 (职业作家)
★ ★
kuhailangyu(金币+2,VIP+0):替楼主向你致谢！ 4-22 16:05
俺找到问题了，其实很简单的问题。

在主文件中规定了 n 的取值，即：

parameter(n=1,n0=18000,n1=7,pi=3.14159265)

这样，主程序中的

       do ii=1,n
         read(1,*,iostat=error) c
      
         open(111,file='E:\'//c//'.txt')

只打开了第一个文件，而这个文件只有一个小时的数据，也就是当

         do j=1,48  !这里j为一天24小时，共48个半小时循环

中的 j 变为 2 时，就到了文件的结尾，但循环还在继续，所以后面的结果都是错的。

针对楼主提供的数据，把

        do ii=1,n  

中的 n 改为 3, 或者改 parameter 那行中的 n=3，然后将

do j=1,48 改为
do j=1,2

即可。

[ Last edited by snoopyzhao on 2009-4-22 at 15:45 ]
赞一下(10人) 回复此楼
6楼2009-04-22 15:43:41
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

老虎大王

木虫 (著名写手)
sunxiao(金币+0,VIP+0):呵呵，新虫，多包涵了 4-22 22:55
晕。楼主发了帖子就没有来过。咱们还研究个啥劲。
赞一下(10人) 回复此楼
7楼2009-04-22 16:06:32
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

yhwsmile

金虫 (小有名气)
引用回帖:
Originally posted by 老虎大王 at 2009-4-22 16:06:
晕。楼主发了帖子就没有来过。咱们还研究个啥劲。

最近去内蒙古出差了，所以没有及时看大家的回复，先谢谢各位了啊，我根据你们的好好看看啊！！
赞一下 回复此楼
8楼2009-04-25 21:23:08
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

yhwsmile

金虫 (小有名气)
上次我上传到邮箱中的数据为每小时的数据，后来经过研究，数据应该为日以上数据，所以我重新写了下，由于数据比较大，我只上传了一天的数据！！
txt文件在邮箱www.126.com
126：yhfortran
密码：123456
再次感谢以上虫友，在讨论啊！！！
赞一下 回复此楼
9楼2009-04-25 21:50:50
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

snoopyzhao

至尊木虫 (职业作家)
★ ★ ★
sunxiao(金币+3,VIP+0):谢谢参与，欢迎常来仿真编程版 4-26 23:35
给个补丁吧
CODE:
--- test6_old.for        Mon Apr 20 15:53:32 2009
+++ test6_new.for        Sat Apr 25 23:48:32 2009
@@ -8,7 +8,7 @@
c     aav:三次旋转订正后的平均风速；uuvc:三次旋转后的协方差；roll(2):二次旋转的角度！第一次v=0；第二次旋转=0；
c     dds:三次旋转后风速的方差,
      
-         character c*27,cc1*11
+         character c*25,cc1*11
         Character(100) cc
         integer error
         data  ba/99999.0/
@@ -23,7 +23,8 @@


        do ii=1,n  
-         read(1,*,iostat=error) c
+         read(1,'(a)',iostat=error) c
+C        write(*,*)c
       
         open(111,file='E:\'//c//'.txt')
         do j=1,48  !这里j为一天24小时，共48个半小时循环
@@ -62,7 +63,9 @@
         cgm(i)=0.0
         enddo
         do j=1,l
-         read(111,*,IOSTAT=error1) data_time(j),(sonic_data(j,jj),jj=1,7)
+         read(111,5,IOSTAT=error1)data_time(j),(sonic_data(j,jj),jj=1,7)
+5        format(A11,3F9.5,F8.3,2F8.4,F8.2)
+C        write(*,5)data_time(j),(sonic_data(j,jj),jj=1,7)
         if(error1/=0)exit
         do jj=1,n1
               if(abs(sonic_data(j,jj))

[ Last edited by snoopyzhao on 2009-4-25 at 23:57 ]
赞一下(10人) 回复此楼
10楼2009-04-25 23:51:03
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 yhwsmile 的主题更新
121/2返回列表12下一页
普通表情 龙 兔 虎 猫 高级回复 (可上传附件)
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[考研] 299求调剂 +8 15188958825 2026-03-25 8/400 2026-03-29 01:36 by fmesaito
[考研] 各位老师好，我的一志愿为北京科技大学085601材料专硕 +7 Koxui 2026-03-28 7/350 2026-03-29 00:58 by 我是小康
[考研] 315求调剂 +4 akie... 2026-03-28 5/250 2026-03-28 21:05 by zhq0425
[考研] 化学0703 调剂 306分 一志愿211 +4 26要上岸 2026-03-28 4/200 2026-03-28 15:30 by 1018329917
[考研] 复试调剂 +3 raojunqi0129 2026-03-28 3/150 2026-03-28 15:27 by 落睿可思
[考研] 085602 307分 求调剂 +7 不知道叫什么！ 2026-03-26 7/350 2026-03-28 09:57 by 神马都不懂
[考研] 265求调剂 +8 小木虫085600 2026-03-27 8/400 2026-03-27 22:16 by 无际的草原
[考研] 272求调剂 +7 脚滑的守法公民 2026-03-27 7/350 2026-03-27 17:23 by laoshidan
[考研] 0856调剂 +5 求求让我有书读� 2026-03-26 6/300 2026-03-27 15:12 by caszguilin
[考研] 279 分 求调剂 +4 睡个好觉_16 2026-03-24 4/200 2026-03-27 15:05 by 醉在风里
[考研] 314求调剂 +3 溪云珂 2026-03-26 3/150 2026-03-27 11:20 by sanrepian
[考研] 调剂求收留 +7 果然有我 2026-03-26 7/350 2026-03-27 00:26 by wxiongid
[考研] 321求调剂 +6 Ymlll 2026-03-24 6/300 2026-03-26 20:50 by 不吃魚的貓
[考研] 【双一流院校新能源、环境材料，材料加工与模拟招收大量调剂】 +4 Higraduate 2026-03-22 8/400 2026-03-26 20:34 by Higraduate
[考研] 生物学 296 求调剂 +4 朵朵- 2026-03-26 6/300 2026-03-26 19:01 by 不吃魚的貓
[考研] 085602 289分求调剂 +8 WWW西西弗斯 2026-03-24 8/400 2026-03-26 16:33 by 不吃魚的貓
[考研] 085601求调剂总分293英一数二 +4 钢铁大炮 2026-03-24 4/200 2026-03-26 16:28 by dick_runner
[考研] 263求调剂 +6 yqdszhdap－ 2026-03-22 10/500 2026-03-26 13:11 by 公瑾逍遥
[考研] 材料专硕 335 分求调剂 +4 拒绝冷暴力 2026-03-25 4/200 2026-03-25 18:45 by haxia
[考研] 336化工调剂 +4 王大坦1 2026-03-23 5/250 2026-03-23 18:32 by allen-yin
信息提示
关闭
请填处理意见
关闭