Fortran迭代

faih

想通过wu_sp和eu_sp的两个公式进行迭代，但是自己写出来的迭代前后好像是一样的。有没有大神可以请教一下呀（哭泣


program main
implicit none

integer :: i,k,j,kk,it
integer :: msg = 0                    !  ic number of missing moisture levels at the top of model.
integer :: ncld = 1
integer :: il2g = 3                   !CORE GROUP REMOVE
integer :: jb(3) = (/21,4,3/)                  ! updraft base level

real,allocatable :: q(:,:)                    ! spec. humidity of env
real :: rl = 2500                             ! latent heat of vap
real,allocatable :: t(:,:)                    ! temp of env
real,allocatable :: s(:,:)                    ! normalized dry static energy of env
real,allocatable :: zf(:,:)


real,allocatable :: th(:,:)
real,allocatable :: qh(:,:)
real,allocatable :: dz(:,:)

real,allocatable :: wu_sp(:,:,:)              !vertical velocity of updraft
real,allocatable :: zkine_sp(:,:,:)
real,allocatable :: tu_sp(:,:,:)
real,allocatable :: zbuo_sp(:,:,:)
real,allocatable :: qu_sp(:,:,:)              ! spec hum of updraft
real,allocatable :: su_sp(:,:,:)              ! normalized dry stat energy of updraft
real,allocatable :: hu_sp(:,:,:)
real,allocatable :: mu_sp(:,:,:)               ! updraft mass flux
real,allocatable :: eu_sp(:,:,:)               ! entrainment rate of updraft
real,allocatable :: eu_spo(:,:,:)

real,parameter ::  retv = 0.608
real,parameter ::  grav = 9.8                  ! = gravit
real,parameter ::  cp = 1004                 !1.004                  ! heat capacity of dry air
real,parameter ::  arb = 0.167
real,parameter ::  ceu = 0.225

real mum_sp
real zbuom_sp
real wum_sp
real wu_sp_0

allocate(q(3,20))
allocate(t(3,20))
allocate(s(3,20))
allocate(zf(3,20))
allocate(th(3,21))
allocate(qh(3,21))
allocate(dz(3,21))
allocate(wu_sp(3,21,1))
allocate(zkine_sp(3,21,1))
allocate(tu_sp(3,21,1))
allocate(zbuo_sp(3,21,1))
allocate(qu_sp(3,21,1))
allocate(su_sp(3,21,1))
allocate(hu_sp(3,21,1))
allocate(mu_sp(3,21,1))
allocate(eu_sp(3,21,1))
allocate(eu_spo(3,21,1))

do i = 1,il2g
kk = jb(i)
!write(*,*) kk
     do j = 1, ncld
            mu_sp(i,kk,j) = 1
        wu_sp(i,kk,j) = 1
                zkine_sp(i,kk,j) = 0.5
                qu_sp(i,kk,j) = 0.02 !0.025
        hu_sp(i,kk,j) = 350000 !350
        tu_sp(i,kk,j) = 280  !290
        th(i,kk) = 278       !280
        qh(i,kk) = 0.03      !water vapor mixing ratio
        dz(i,kk) = 10  !10
        su_sp(i,kk,j) = (hu_sp(i,kk,j)-rl*qu_sp(i,kk,j))/cp
                zbuo_sp(i,kk,j) = grav*(tu_sp(i,kk,j)*(1+retv*qu_sp(i,kk,j))-    &
                    th(i,kk)*(1+retv*qh(i,kk)))/   &
                    (th(i,kk)*(1+retv*qh(i,kk)))
!write(*,*) su_sp(i,kk,j)      
!write(*,*) zbuo_sp(i,kk,j)
                eu_sp(i,kk,j) = mu_sp(i,kk,j)/dz(i,kk)
                do k= jb(i)-1,msg + 1,-1
              do it =1,2
                 if (it.eq.1) then
                    mum_sp =  mu_sp(i,k+1,j)
                    zbuom_sp= zbuo_sp(i,k+1,j)
                    wum_sp = wu_sp(i,k+1,j)
                 end if
!write(*,*) wum_sp
                s(i,k) = 285
                dz(i,k) = 10 !10
                q(i,k) = 0.03  !0.035
                zf(i,k) = 1000-10*k ! 100*(jb(i)-k) !0.5  
                t(i,k) = 278   !275
                                !eu_sp(i,k,j)= -1.2+ 3.8/wum_sp-0.36/(wum_sp**2)
                !eu_sp(i,k,j)= arb*ceu*zbuom_sp/(wum_sp**2)


                do while (.true.)
                eu_sp(i,k,j)= 0.15*(wum_sp**(-0.14))*(zbuom_sp**(-0.10))
                !eu_sp(i,k,j)= -1.2+ 3.8/wum_sp-0.36/(wum_sp**2)
                !eu_sp(i,k,j)= 0.23-1.22/wum_sp+10/(wum_sp**2)

!write(*,*) eu_sp(i,k,j)


                                eu_sp(i,k,j)= max(eu_sp(i,k,j),1.0e-5)
                eu_sp(i,k,j)= min(eu_sp(i,k,j),2.0e-3)  
                eu_spo(i,k,j)= eu_sp(i,k,j)
!write(*,*) eu_sp(i,k,j)
                eu_sp(i,k,j)= eu_sp(i,k,j)* mum_sp
!write(*,*) eu_sp(i,k,j)

                                mu_sp(i,k,j)= mu_sp(i,k+1,j)+ eu_sp(i,k,j)*dz(i,k)
!write(*,*) mu_sp(i,k,j)
                                if (mu_sp(i,k,j)==0) write(*,*) "mu_sp=0","i=",i,"j=",j,"k=",k
                                su_sp(i,k,j)= mu_sp(i,k+1,j)/mu_sp(i,k,j)*su_sp(i,k+1,j) +  &
                               eu_sp(i,k,j)*dz(i,k)*s(i,k)/mu_sp(i,k,j)
!write(*,*) su_sp(i,k,j)
                                qu_sp(i,k,j)= mu_sp(i,k+1,j)/mu_sp(i,k,j)*qu_sp(i,k+1,j) +  &
                               eu_sp(i,k,j)*dz(i,k)*q(i,k)/mu_sp(i,k,j)
!write(*,*) qu_sp(i,k,j)
                tu_sp(i,k,j)= su_sp(i,k,j) - grav/cp*zf(i,k)
!write(*,*) tu_sp(i,k,j)
                                zbuo_sp(i,k,j) = grav*(0.5*(tu_sp(i,k,j)+tu_sp(i,k+1,j))*    &
                              (1+retv*0.5*(qu_sp(i,k,j)+qu_sp(i,k+1,j)))-    &
                              t(i,k)*(1+retv*q(i,k)))/   &
                              (t(i,k)*(1+retv*q(i,k)))
!write(*,*) zbuo_sp(i,k,j)
!write(*,*) wum_sp

                !zkine_sp(i,k,j) = zkine_sp(i,k+1,j) + (arb*zbuo_sp(i,k,j)+1.2*wum_sp**(2)-3.8*wum_sp+0.36)*dz(i,k)
                !zkine_sp(i,k,j) = zkine_sp(i,k+1,j) + (arb*zbuo_sp(i,k,j)-eu_sp(i,k,j)*wu_sp_0**(2))*dz(i,k)
                zkine_sp(i,k,j) = (zkine_sp(i,k+1,j) + arb*zbuo_sp(i,k,j)*dz(i,k))/(1+2*eu_sp(i,k,j)*dz(i,k))
                                wu_sp(i,k,j) = min(15.0,sqrt(2.0*max(0.1,zkine_sp(i,k,j)) ))
                if (abs(wu_sp(i,k,j)- wum_sp)<0.001) then
                    exit
                else   
                 wum_sp = wu_sp(i,k,j)               
                end if
!!!!write(*,*) zkine_sp(i,k,j)
                end do

!write(*,*) zkine_sp(i,k,j)

!write(*,*) wu_sp(i,k,j)
!write(*,*) zkine_sp(i,k,j)
                                mum_sp = 0.5*(mu_sp(i,k+1,j)+mu_sp(i,k,j))
                zbuom_sp= 0.5*(zbuo_sp(i,k+1,j)+zbuo_sp(i,k,j))
                wum_sp = 0.5*(wu_sp(i,k+1,j)+wu_sp(i,k,j))

!write(*,*)  mum_sp,zbuom_sp,c
!write(*,*)  wu_sp(i,k,j),eu_sp(i,k,j)
open(unit=10,file='D:\test1\3.txt')
write(10,*)  i,k,wu_sp(i,k,j),eu_sp(i,k,j)
                          end do
         end do  !k
     end do  !j
end do
end