在使用wrf时一不小心发现nssl_2mom方案中的一个错误

lxc0818

第一次发帖有点紧张

在使用wrf模式模拟一次飑线的过程中，我选取了编号18的微物理方案，即nssl 2mom ,含CCN prediction。最初模拟的结果感觉很好， 直到有一天想观赏一下这个微物理方案的代码，然后发现这个16000多行的代码中的一个小错误，应该错了吧，我想。认为没错的请看这一段代码：

!  Subroutine to calculate number concentrations from initial state that has only mixing ratio.
!  N will be in #/kg, NOT #/m^3, since sedimentation is done next.
!
      subroutine calcnfromq(nx,ny,nz,an,na,nor,norz,dn)
      
      implicit none
      integer nx,ny,nz,nor,norz,na,ngt,jgs,ixcol
      real an(-nor+1:nx+nor,-nor+1:ny+nor,-norz+1:nz+norz,na)  ! scalars (q, N, Z)
      real dn(nx,nz+1)  ! air density
      
      integer ixe,kze
      real    alpha
      real    qmin
      real    xvmn,xvmx
      integer ipconc
      integer lvol ! index for volume
      integer infall
      
      
      integer ix,jy,kz
      double precision vr,q,nrx,rd,g1h,g1r,g1s,zx,chw,z,znew,zt,zxt,n1,laminv1
      double precision :: zr, zs, zh, dninv
      real, parameter :: xn0s = 3.0e6, xn0r = 8.0e6, xn0h = 4.0e4
      real, parameter :: xdnr = 1000., xdns = 100. ,xdnh = 900.0
      real, parameter :: zhfac = 1./(pi*xdnh*xn0h)
      real, parameter :: zrfac = 1./(pi*xdnr*xn0r)
      real, parameter :: zsfac = 1./(pi*xdns*xn0s)
      real, parameter :: g0 = (6.0)*(5.0)*(4.0)/((3.0)*(2.0)*(1.0))
      real xv,xdn
      integer :: ndbz, nmwgt, nnwgt, nwlessthanz
! ------------------------------------------------------------------
      
      
      jy = 1
      
      
         g1h = (6.0 + alphah)*(5.0 + alphah)*(4.0 + alphah)/  &
     &        ((3.0 + alphah)*(2.0 + alphah)*(1.0 + alphah))
     
         IF ( imurain == 3 ) THEN
         g1r = (rnu+2.0)/(rnu+1.0)
         ELSE ! imurain == 1
         g1r = (6.0 + alphar)*(5.0 + alphar)*(4.0 + alphar)/  &
     &        ((3.0 + alphar)*(2.0 + alphar)*(1.0 + alphar))
         ENDIF
         g1s = (snu+2.0)/(snu+1.0)
      
      DO kz = 1,nz
       DO ix = 1,nx ! ixcol
         dninv = 1./dn(ix,kz)
         
   !  Cloud droplets
         
         IF ( lnc > 1 ) THEN
           IF ( an(ix,jy,kz,lnc) <= 0.0 .and. an(ix,jy,kz,lc) > qxmin(lc) ) THEN
             an(ix,jy,kz,lnc) = qccn
           ENDIF
         ENDIF
   !  rain
         
         IF ( lnr > 1 ) THEN
           IF ( an(ix,jy,kz,lnr) <= 0.0 .and. an(ix,jy,kz,lr) > qxmin(lr) ) THEN
             q = an(ix,jy,kz,lr)
            
             laminv1 = (dn(ix,kz) * q * zrfac)**(0.25)  ! inverse of slope
            
             n1 = laminv1*xn0h  ! number concentration for inv. exponential single moment input
            
             nrx =  n1*g1r/g0   ! number concentration for different shape parameter
             an(ix,jy,kz,lnr) = nrx ! *dninv ! convert to number mixing ratio
            
           ENDIF
         ENDIF
  ! snow
         IF ( lns > 1 ) THEN
           IF ( an(ix,jy,kz,lns) <= 0.0 .and. an(ix,jy,kz,ls) > qxmin(ls) ) THEN
             q = an(ix,jy,kz,ls)
            
             laminv1 = (dn(ix,kz) * q * zsfac)**(0.25)  ! inverse of slope
            
             n1 = laminv1*xn0s  ! number concentration for inv. exponential single moment input
            
             nrx =  n1*g1s/g0   ! number concentration for different shape parameter
             an(ix,jy,kz,lns) = nrx ! *dninv ! convert to number mixing ratio
            
           ENDIF
         ENDIF
         
    ! graupel
         IF ( lnh > 1 ) THEN
           IF ( an(ix,jy,kz,lnh) <= 0.0 .and. an(ix,jy,kz,lh) > qxmin(lh) ) THEN
             IF ( lvh > 1 ) THEN
               IF ( an(ix,jy,kz,lvh) <= 0.0 ) THEN
                 an(ix,jy,kz,lvh) = an(ix,jy,kz,lh)/xdnh
               ENDIF
             ENDIF
             q = an(ix,jy,kz,lh)
            
             laminv1 = (dn(ix,kz) * q * zhfac)**(0.25)  ! inverse of slope
            
             n1 = laminv1*xn0h  ! number concentration for inv. exponential single moment input
            
             nrx =  n1*g1h/g0   ! number concentration for different shape parameter
             an(ix,jy,kz,lnh) = nrx ! *dninv ! convert to number mixing ratio
           ENDIF
         ENDIF

      ENDDO ! ix
      ENDDO ! kz
      
      RETURN
      
      END subroutine calcnfromq


仔细一看，计算霰粒子和雨滴居然用的是同一个初始值。。。这明细是应该有问题的。。。也就是说我的毕业论文的模拟结果是靠不住的，wocao。
这直接差了2个量级额。

附个16000多行的代码。

lxc0818

自己水自己

lxc0818

Mansell回复说：You are correct!Thank you for pointing this out.....

李舒文

只能说，楼主很细心很较劲，简单粗暴已经，楼主很牛呀~！

李舒文

只能说，楼主很细心很较劲，简单粗暴一句，楼主很牛呀~！

freekiller

很多方案都有bug。多数时候影响不明显没有被发现，你可以改对了试试你之前的模拟，你可能会发现几乎没差别。

qtang1990

llllllllllll

lxc0818

freekiller 发表于 2018-1-24 21:40
很多方案都有bug。多数时候影响不明显没有被发现，你可以改对了试试你之前的模拟，你可能会发现几乎没差别 ...

雨滴数浓度有差几个量级

freekiller

lxc0818 发表于 2018-4-16 10:13
雨滴数浓度有差几个量级

模拟的降水差异大吗？强度和降水落区。谢谢

lxc0818

freekiller 发表于 2018-4-16 10:46
模拟的降水差异大吗？强度和降水落区。谢谢

没具体比较，但对天气系统的结构以及空间范围影响很大。但是在WRF3.7以后的版本中此贴中的错误已经得到纠正。