fortran扩散编程

【实例简介】2d_convection_diffusion.f90
【实例截图】
【核心代码】program convection_diffusion
  integer :: i,j,k,ni,nj,nt
  real    :: dx,dt
  real    :: KK
  real    :: temp1,temp2,temp3,temp4
  real,allocatable :: x(:,:),y(:,:),c(:,:),c_old(:,:),U(:,:),V(:,:)
  real    :: x_coner
  real    :: y_coner
  real    :: Lh
  real    :: Lv
  !
  open(10,file='parameter_setup.dat')
  read(10,*)ni,nj
  read(10,*)nt
  read(10,*)dt
  read(10,*)x_coner,y_coner
  read(10,*)Lh,Lv
  read(10,*)KK
  close(10)
  !
  allocate(x(ni,nj),y(ni,nj),U(ni,nj),V(ni,nj),c(ni,nj),c_old(ni,nj))
  !
  dx = Lh/(ni-1)
  dy = Lv/(nj-1)
  !
  do j=1,nj
  do i=1,ni
     x(i,j)=x_coner (i-1)*dx
     y(i,j)=y_coner (j-1)*dy
     !
     U(i,j)=  0.0005*y(i,j)
     V(i,j)= -0.0005*x(i,j)
     !
  enddo
  enddo
  !
  c= 0.0
  c_old = 0.0
  !
  do j=1,nj
  do i=1,ni
     if(x(i,j) >=-280.0 .and. x(i,j) <=-100.0  .and.y(i,j) >=-280.0 .and. y(i,j) <=-100.0  )then
        c(i,j) = 200.0
        c_old(i,j)=200.0
     endif
  enddo
  enddo
  !
  open(10,file='initial_condition.plt')
  write(10,*) 'TITLE ="HSVA-3D"'
  write(10,*) 'VARIABLES= X, Y, C, U, V'
  write(10,*) 'ZONE T="zone1", I=',ni,', J=',nj, ' F=block'
  write(10,*) ((x(i,j),i=1,ni),j=1,ni)
  write(10,*) ((y(i,j),i=1,ni),j=1,nj)
  write(10,*) ((c(i,j),i=1,ni),j=1,nj)
  write(10,*) ((U(i,j),i=1,ni),j=1,nj)
  write(10,*) ((V(i,j),i=1,ni),j=1,nj)
  close(10)
  !
  do k=1,nt
    !
      write(*,*) 'time_step= ',k
    !      
     do j=2,nj-1
     do i=2,ni-1
        temp1 = KK*( c_old(i 1,j) - 2.*c_old(i,j) c_old(i-1,j) )/dx/dx
        temp2 = KK*( c_old(i,j 1) - 2.*c_old(i,j) c_old(i,j-1) )/dy/dy
        !
        if(U(i,j) >0.0) then
          temp3 = -U(i,j)*( c_old(i,j) - c_old(i-1,j) )/dx
        else
          temp3 = -U(i,j)*( c_old(i 1,j) - c_old(i,j) )/dx
        endif
        !
        if(V(i,j) >0.0) then
          temp4 = -V(i,j)*( c_old(i,j) - c_old(i,j-1) )/dy
        else
          temp4 = -V(i,j)*( c_old(i,j 1) - c_old(i,j) )/dy
        endif
        !
        c(i,j) = c_old(i,j) dt* (temp1 temp2 temp3 temp4)
        !
     enddo
     enddo
     !
     !boundary treatment
     !
     do i=1,ni
        c(i,1) = 0.0
        c(i,nj) = 0.0
     enddo
     !
      do j=1,nj
        c(1,j) = 0.0
        c(ni,j) = 0.0
     enddo
     !
     c_old = c
     !
  enddo
  !
  open(10,file='rsults.plt')
  write(10,*) 'TITLE ="HSVA-3D"'
  write(10,*) 'VARIABLES= X, Y, C'
  write(10,*) 'ZONE T="zone1", I=',ni,', J=',nj, ' F=block'
  write(10,*) ((x(i,j),i=1,ni),j=1,ni)
  write(10,*) ((y(i,j),i=1,ni),j=1,nj)
  write(10,*) ((c(i,j),i=1,ni),j=1,nj)
  write(10,*) ((U(i,j),i=1,ni),j=1,nj)
  write(10,*) ((V(i,j),i=1,ni),j=1,nj)
  close(10)
  !
  pause 'Finish!'
  !
end