fortran的流函数 势函数问题

cycycycycy

想了三天的流函数势函数，结果还是不对{:cry:}才10^2,我知道肯定是10^7左右，思路好乱{:cry:}已经黔驴技穷了，站点里面现有的13年版的fortran程序，渣基础看不懂哇{:cry:}
我用的流函数第一边条件，和利布曼法，不知道为什么有问题{:cry:}
应该是思路的问题，有大神在线帮忙看一下，或者有让我借鉴的吗
十分感谢啊！真的脑子要炸啦
!流函数steam的边条件
dingzh=0.0
loop=0.0
loop_abs=0.0
do k=1,mz
  do t=1,mt
    do j=1,my
          do i=1,mx
           !速度的环路积分 loop(i,j,k,t)  loop_abs(i,j,k,t)
           loop(i,j,k,t)=loop(i,j,k,t)+a*delta*pi/180*(cos(lat(j))*(-q_v(i,1,k,t)+q_v(i,my,k,t))+(q_u(mx,j,k,t)-q_u(1,j,k,t)))
           !速度绝对值的环路积分
           loop_abs(i,j,k,t)=loop_abs(i,j,k,t)+a*delta*pi/180*(cos(lat(j))*(abs(-q_v(i,1,k,t))+abs(q_v(i,my,k,t)))+(abs(q_u(mx,j,k,t))+abs(-q_u(1,j,k,t))))
           enddo
    enddo
        !东南西北边界的订正系数
    dingzh(k,t)=-loop(mx,my,k,t)/loop_abs(mx,my,k,t)
  enddo
enddo
!订正流函数的u，v分量
do k=1,mz
  do t=1,mt
    do j=1,my
          do i=1,mx
            q_uu(i,j,k,t)=q_u(i,j,k,t)+dingzh(k,t)*abs(q_u(i,j,k,t))
                q_vv(i,j,k,t)=q_v(i,j,k,t)+dingzh(k,t)*abs(q_v(i,j,k,t))
      enddo
        enddo
  enddo
enddo
!边界上的steam值,采用中央差分
steam=0.0
do k=1,mz
  do t=1,mt
        do i=1+1,mx
        !j=my上的steam
    steam(i,my,k,t)=steam(i-1,my,k,t)+(q_vv(i-1,my,k,t)+q_vv(i,my,k,t))/2*a*delta*pi/180*cos(lat(my))
    enddo
  enddo
enddo
do k=1,mz
  do t=1,mt
    do j=my-1,1,-1
        !i=1.and.i=mx上的steam
          steam(1,j,k,t)=steam(1,j+1,k,t)+(q_uu(1,j+1,k,t)+q_uu(1,j,k,t))/2*a*delta*pi/180
          steam(mx,j,k,t)=steam(1,j+1,k,t)+(q_uu(1,j+1,k,t)+q_uu(1,j,k,t))/2*a*delta*pi/180
    enddo
  enddo
enddo
do k=1,mz
  do t=1,mt
    do i=1+1,mx
        !j=1上的steam
    steam(i,1,k,t)=steam(i-1,1,k,t)+(q_vv(i-1,1,k,t)+q_vv(i,1,k,t))/2*a*delta*pi/180*cos(lat(1))
    enddo
  enddo
enddo
print*,((steam(i,j,11,1),i=1,mx),j=1,my)
!不对啊不对啊，1.0e+7啊采用利布曼法，减少迭代次数，节省存储空间，不必要同时分配函数的n步和n+1步。
RR=0.0
RR1=1.e-6
do k=1,mz
do t=1,mt
  do while(RR1(1,1,k,t)>1.e-7)
  n=n+1
        do j=1+1,my-1
          do i=1+1,mx-1
                !初始估计值不是真值，与该点的散度存在偏差RR，称为残差或余差
             RR(i,j,k,t)=(steam(i+1,j,k,t)+steam(i-1,j,k,t))/((a*delta*pi/180*cos(lat(j)))**2)+(steam(i,j+1,k,t)+steam(i,j-1,k,t))/((a*delta*pi/180)**2)-(2.0/(a*delta*pi/180*cos(lat(j)))**2+2.0/(a*delta*pi/180)**2)*steam(i,j,k,t)-vapor_voro(i,j,k,t)
             steam(i,j,k,t)=steam(i,j,k,t)+RR(i,j,k,t)/(2*(1./(a*delta*pi/180*cos(lat(j))**2))+1./((a*delta*pi/180)**2))     
      enddo
        enddo
    RR1=abs(RR)
    do j=1,my
      do i=1,mx
        if(RR1(1,1,k,t)<=RR1(i,j,k,t))then
          RR1(1,1,k,t)=RR1(i,j,k,t)
             endif
      enddo
    enddo
  enddo
enddo
enddo

季报台风

我也不明白咋回事

cycycycycy

又重新写了一下，发现边值也不对啊啊啊啊啊啊啊
!------计算流函数steam function------
!流函数steam的边条件
dingzh=0.0
loop=0.0
loop_abs=0.0
do k=1,mz
  do t=1,mt
    do j=1,my
          do i=1,mx
           !速度的环路积分 loop(i,j,k,t)  loop_abs(i,j,k,t)
           loop(i,j,k,t)=loop(i,j,k,t)+a*delta*pi/180*(cos(lat(j))*(-q_v(i,1,k,t)+q_v(i,my,k,t))+(q_u(mx,j,k,t)-q_u(1,j,k,t)))
           !速度绝对值的环路积分
           loop_abs(i,j,k,t)=loop_abs(i,j,k,t)+a*delta*pi/180*(cos(lat(j))*(abs(-q_v(i,1,k,t))+abs(q_v(i,my,k,t)))+(abs(q_u(mx,j,k,t))+abs(-q_u(1,j,k,t))))
           enddo
    enddo
        !东南西北边界的订正系数
    dingzh(k,t)=-loop(mx,my,k,t)/loop_abs(mx,my,k,t)
  enddo
enddo
!订正流函数的u，v分量
do k=1,mz
  do t=1,mt
    do j=1,my
          do i=1,mx
            q_uu(i,j,k,t)=q_u(i,j,k,t)+dingzh(k,t)*abs(q_u(i,j,k,t))
                q_vv(i,j,k,t)=q_v(i,j,k,t)+dingzh(k,t)*abs(q_v(i,j,k,t))
      enddo
        enddo
  enddo
enddo
!边界上的steam值,采用两点差分
steam=0.0
do k=1,mz
  do t=1,mt
        do i=1+1,mx
        !j=my上的steam
    steam(i,my,k,t)=steam(i-1,my,k,t)+(q_vv(i-1,my,k,t)+q_vv(i,my,k,t))/2*a*delta*pi/180*cos(lat(my))
    enddo
  enddo
enddo
do k=1,mz
  do t=1,mt
    do j=my-1,1,-1
        !i=1.and.i=mx上的steam
          steam(1,j,k,t)=steam(1,j+1,k,t)+(q_uu(1,j+1,k,t)+q_uu(1,j,k,t))/2*a*delta*pi/180
          steam(mx,j,k,t)=steam(1,j+1,k,t)+(q_uu(1,j+1,k,t)+q_uu(1,j,k,t))/2*a*delta*pi/180
    enddo
  enddo
enddo
do k=1,mz
  do t=1,mt
    do i=1+1,mx
        !j=1上的steam
    steam(i,1,k,t)=steam(i-1,1,k,t)+(q_vv(i-1,1,k,t)+q_vv(i,1,k,t))/2*a*delta*pi/180*cos(lat(1))
    enddo
  enddo
enddo
!不对啊不对啊，1.0e+7啊采用利布曼法，减少迭代次数，节省存储空间，不必要同时分配函数的n步和n+1步。
RR=0.0
RR1=1.e-6
do k=1,mz
do t=1,mt
  do while(RR1(1,1,k,t)>1.e-7)
  n=n+1
        do j=1+1,my-1
          do i=1+1,mx-1
                !初始估计值不是真值，与该点的散度存在偏差RR，称为残差或余差
             RR(i,j,k,t)=(steam(i+1,j,k,t)+steam(i-1,j,k,t))/((a*delta*pi/180*cos(lat(j)))**2)+(steam(i,j+1,k,t)+steam(i,j-1,k,t))/((a*delta*pi/180)**2)-(2.0/(a*delta*pi/180*cos(lat(j)))**2+2.0/(a*delta*pi/180)**2)*steam(i,j,k,t)-vapor_voro(i,j,k,t)
             steam(i,j,k,t)=steam(i,j,k,t)+RR(i,j,k,t)/(2*(1./((a*delta*pi/180*cos(lat(j)))**2)+1./((a*delta*pi/180)**2)))     
      enddo
        enddo
    RR1=abs(RR)
    do j=1,my
      do i=1,mx
        if(RR1(1,1,k,t)<=RR1(i,j,k,t))then
          RR1(1,1,k,t)=RR1(i,j,k,t)
             endif
      enddo
    enddo
  enddo
enddo
enddo
!计算函数potential function
potential=0.0
RR=0.0
RR1=1.e-6
do k=1,mz
do t=1,mt
  do while(RR1(1,1,k,t)>1.e-7)
  n=n+1
        do j=1+1,my-1
          do i=1+1,mx-1
          !初始估计值不是真值，与该点的散度存在偏差RR，称为残差或余差
          RR(i,j,k,t)=(potential(i+1,j,k,t)+potential(i-1,j,k,t))/((a*delta*pi/180*cos(lat(j)))**2)+(potential(i,j+1,k,t)+potential(i,j-1,k,t))/((a*delta*pi/180)**2)-(2.0/(a*delta*pi/180*cos(lat(j)))**2+2.0/(a*delta*pi/180)**2)*potential(i,j,k,t)+div(i,j,k,t)
          enddo
        enddo
         do j=1+1,my-1
            do i=1+1,mx-1
                 potential(i,j,k,t)=potential(i,j,k,t)+RR(i,j,k,t)/(2*(1./((a*delta*pi/180*cos(lat(j)))**2)+1./((a*delta*pi/180)**2)))     
        enddo
         enddo
        do j=1,my
          do i=1,mx
        RR1(i,j,k,t)=abs(RR(i,j,k,t))
          enddo
        enddo
    do j=1,my
      do i=1,mx
        if(RR1(1,1,k,t)<=RR1(i,j,k,t))then
          RR1(1,1,k,t)=RR1(i,j,k,t)
             endif
      enddo
    enddo
  enddo
enddo
enddo

cycycycycy

下面这个应该是对的吧，如果有人的话就帮我看一下吧
!------计算流函数steam function------
!流函数steam的边条件
dingzh=0.0
loop=0.0
loop_abs=0.0
do k=1,mz
  do t=1,mt
    do j=1,my
          do i=1,mx
           !速度的环路积分 loop(i,j,k,t)  loop_abs(i,j,k,t)
           loop(i,j,k,t)=loop(i,j,k,t)+a*delta*pi/180*(cos(lat(j))*(-q_v(i,1,k,t)+q_v(i,my,k,t))+(q_u(mx,j,k,t)-q_u(1,j,k,t)))
           !速度绝对值的环路积分
           loop_abs(i,j,k,t)=loop_abs(i,j,k,t)+a*delta*pi/180*(cos(lat(j))*(abs(-q_v(i,1,k,t))+abs(q_v(i,my,k,t)))+(abs(q_u(mx,j,k,t))+abs(-q_u(1,j,k,t))))
           enddo
    enddo
        !东南西北边界的订正系数
    dingzh(k,t)=-loop(mx,my,k,t)/loop_abs(mx,my,k,t)
  enddo
enddo
!订正流函数的u，v分量
do k=1,mz
  do t=1,mt
    do j=1,my
          do i=1,mx
            q_uu(i,j,k,t)=q_u(i,j,k,t)+dingzh(k,t)*abs(q_u(i,j,k,t))
                q_vv(i,j,k,t)=q_v(i,j,k,t)+dingzh(k,t)*abs(q_v(i,j,k,t))
      enddo
        enddo
  enddo
enddo
!边界上的steam值
steam=0.0
do k=1,mz
  do t=1,mt
        do i=1+1,mx
        !j=my上的steam
    steam(i,my,k,t)=steam(i-1,my,k,t)+(q_vv(i-1,my,k,t)+q_vv(i,my,k,t))/2*a*delta*pi/180*cos(lat(my))
    enddo
  enddo
enddo
do k=1,mz
  do t=1,mt
    do j=my-1,1,-1
        !i=1.and.i=mx上的steam
          steam(1,j,k,t)=steam(1,j+1,k,t)+(q_uu(1,j+1,k,t)+q_uu(1,j,k,t))/2*a*delta*pi/180
          steam(mx,j,k,t)=steam(1,j+1,k,t)+(q_uu(1,j+1,k,t)+q_uu(1,j,k,t))/2*a*delta*pi/180
    enddo
  enddo
enddo
do k=1,mz
  do t=1,mt
    do i=1+1,mx
        !j=1上的steam
    steam(i,1,k,t)=steam(i-1,1,k,t)+(q_vv(i-1,1,k,t)+q_vv(i,1,k,t))/2*a*delta*pi/180*cos(lat(1))
    enddo
  enddo
enddo

!1.0e+3采用利布曼法，减少迭代次数，节省存储空间，不必要同时分配函数的n步和n+1步。
!先将所有RR(I,J)计算一遍，再计算steam,-->判断RR是否全部小于<1.e-7
!将每步计算后的残差RR赋值到RR1,寻找最大值，若MAX<1.e-7,则停止迭代
RR=0.0
RR1=1.e-6
do k=1,mz
do t=1,mt
  do while(RR1(1,1,k,t)>1.e-7)
  n=n+1 !迭代次数
        do j=1+1,my-1
          do i=1+1,mx-1
                !初始估计值不是真值，与该点的散度存在偏差RR，称为残差或余差
             RR(i,j,k,t)=(steam(i+1,j,k,t)+steam(i-1,j,k,t))/((a*delta*pi/180*cos(lat(j)))**2)+(steam(i,j+1,k,t)+steam(i,j-1,k,t))/((a*delta*pi/180)**2)-(2.0/(a*delta*pi/180*cos(lat(j)))**2+2.0/(a*delta*pi/180)**2)*steam(i,j,k,t)-vapor_voro(i,j,k,t)
          enddo
    enddo
        do j=1+1,my-1
          do i=1+1,mx-1   
                 steam(i,j,k,t)=steam(i,j,k,t)+RR(i,j,k,t)/(2*(1./((a*delta*pi/180*cos(lat(j)))**2)+1./((a*delta*pi/180)**2)))     
      enddo
        enddo
        !find MAX,令RR1(1,1,k,t)=max
    RR1=abs(RR)
    do j=1,my
      do i=1,mx
        if(RR1(1,1,k,t)<=RR1(i,j,k,t))then
          RR1(1,1,k,t)=RR1(i,j,k,t)
             endif
      enddo
    enddo
  enddo
enddo
enddo
!计算函数potential function
potential=0.0
RR=0.0
RR1=1.e-6
n=0
do k=1,mz
do t=1,mt
  do while(RR1(1,1,k,t)>1.e-7)
  n=n+1 !迭代次数
        do j=1+1,my-1
          do i=1+1,mx-1
          !初始估计值不是真值，与该点的散度存在偏差RR，称为残差或余差
          RR(i,j,k,t)=(potential(i+1,j,k,t)+potential(i-1,j,k,t))/((a*delta*pi/180*cos(lat(j)))**2)+(potential(i,j+1,k,t)+potential(i,j-1,k,t))/((a*delta*pi/180)**2)-(2.0/(a*delta*pi/180*cos(lat(j)))**2+2.0/(a*delta*pi/180)**2)*potential(i,j,k,t)+div(i,j,k,t)
          enddo
        enddo
         do j=1+1,my-1
            do i=1+1,mx-1
                 potential(i,j,k,t)=potential(i,j,k,t)+RR(i,j,k,t)/(2*(1./((a*delta*pi/180*cos(lat(j)))**2)+1./((a*delta*pi/180)**2)))     
        enddo
         enddo
        do j=1,my
          do i=1,mx
        RR1(i,j,k,t)=abs(RR(i,j,k,t))
          enddo
        enddo
    do j=1,my
      do i=1,mx
        if(RR1(1,1,k,t)<=RR1(i,j,k,t))then
          RR1(1,1,k,t)=RR1(i,j,k,t)
             endif
      enddo
    enddo
  enddo
enddo
enddo

登录登登录登登

感谢楼主分享~

yeah...

顶一个！！！！！！！！！！！我也想问

风儿飘飘

好复杂，看不懂

风儿飘飘

不知计算中缺测值如何处理的

porpoft

请问楼主计算的是哪个范围的势函数和流函数啊？时全球的吗？

yeah...

请问楼主改出来了吗