求助模式地面温度因为凝结持续降温的问题

凌波微步！

运行的是一个云解析模式SAM，各种初值、边界、辐射通量都用ERA5提前给定了。
二维、32格点、4000米分辨率、10秒时间步长。因此1小时内积分360次。
主要关注的是这个模式用的微物理模块SAM1MOM。发现一个严重问题：
每次积分，地面比湿总是大于饱和比湿，持续凝结降温，最后竟然使地面低于0摄氏度，
夏季出现了降雪！除了地面温度，其他层温度没有问题。
这个微物理模块是用温度阈值来分离云水（云中液态水、冰水两类）、雨水（雨、雪、雹三类），
没有显式的雪融化过程，是根据温度来隐式表示的。

求助大佬指点迷津！

主要代码段如下：
        an = 1./(tbgmax-tbgmin)       
        bn = tbgmin * an
        ap = 1./(tprmax-tprmin)       
        bp = tprmin * ap
        fac1 = fac_cond+(1+bp)*fac_fus
        fac2 = fac_fus*ap
        ag = 1./(tgrmax-tgrmin)       

        kmax=0
        kmin=nzm+1
!       print *,'cloud1,tabs(16,1,1)',tabs(16,1,1)
!       print *,'cloud11,tabs(16,1,2)',tabs(16,1,2)
        do k = 1, nzm
         do j = 1, ny
          do i = 1, nx

             q(i,j,k)=max(0.,q(i,j,k))

c
c Initail guess for temperature assuming no cloud water/ice:
c

            tabs(i,j,k) = t(i,j,k)-gamaz(k)
            tabs1=(tabs(i,j,k)+fac1*qp(i,j,k))/(1.+fac2*qp(i,j,k))
!       print *,'initial tabs(16,1,1),tabs1(16,1,1)',tabs(16,1,1),tabs1

c      Warm cloud:

            if(tabs1.ge.tbgmax) then

              tabs1=tabs(i,j,k)+fac_cond*qp(i,j,k)
              qsat(i,j,k) = qsatw_crm(tabs1,pres(k))
              omegan(i,j,k) = 1.

c      Ice cloud:

            elseif(tabs1.le.tbgmin) then

              tabs1=tabs(i,j,k)+fac_sub*qp(i,j,k)             
              qsat(i,j,k) = qsati_crm(tabs1,pres(k))
              omegan(i,j,k) = 0.

c      Mixed-phase cloud:

            else

              om = an*tabs1-bn
              qsat(i,j,k) = om*qsatw_crm(tabs1,pres(k))+
     &                         (1.-om)*qsati_crm(tabs1,pres(k))
              omegan(i,j,k) = om

            endif

c
c  Test if condensation is possible:
c

            if(q(i,j,k).gt.qsat(i,j,k) ) then
!       print *,'(i,j,k)',i,j,k
!       print *,'q(i,j,k)',q(i,j,k)
              niter=0
              dtabs = 100.
              do while(abs(dtabs).gt.0.01.and.niter.lt.10)
                if(tabs1.ge.tbgmax) then
                   om=1.
                   lstarn=fac_cond
                   dlstarn=0.
                   qsat(i,j,k)=qsatw_crm(tabs1,pres(k))
                   dqsat=dtqsatw_crm(tabs1,pres(k))
                else if(tabs1.le.tbgmin) then
                   om=0.
                   lstarn=fac_sub
                   dlstarn=0.
                   qsat(i,j,k)=qsati_crm(tabs1,pres(k))
                   dqsat=dtqsati_crm(tabs1,pres(k))
                else
                   om=an*tabs1-bn
                   lstarn=fac_cond+(1.-om)*fac_fus


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! I changed acoording to SAM6.11.7/cloud.f90  2022.09.21
                   dlstarn=an*fac_fus
!!       print *,'an,an*fac_fus',an,dlstarn
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!111
!                   dlstarn=an
                   qsat(i,j,k)=om*qsatw_crm(tabs1,pres(k))+
     &                        (1.-om)*qsati_crm(tabs1,pres(k))
                   dqsat=om*dtqsatw_crm(tabs1,pres(k))+
     &                        (1.-om)*dtqsati_crm(tabs1,pres(k))
                endif
                if(tabs1.ge.tprmax) then
                   omp=1.
                   lstarp=fac_cond
                   dlstarp=0.
                else if(tabs1.le.tprmin) then
                   omp=0.
                   lstarp=fac_sub
                   dlstarp=0.
                else
                   omp=ap*tabs1-bp
                   lstarp=fac_cond+(1.-omp)*fac_fus
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! I changed acoording to SAM6.11.7/cloud.f90  2022.09.21
                   dlstarp=ap*fac_fus
!!       print *,'ap,ap*fac_fus',ap,dlstarp
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!111
!                   dlstarp=ap

                endif
                fff = tabs(i,j,k)-tabs1+
     &                lstarn*(q(i,j,k)-qsat(i,j,k))+lstarp*qp(i,j,k)
                dfff=dlstarn*(q(i,j,k)-qsat(i,j,k))+dlstarp*qp(i,j,k)-
     &                                 lstarn*dqsat-1.
                dtabs=-fff/dfff
                niter=niter+1
                tabs1=tabs1+dtabs
      if(i.eq.16. .and. j .eq.1 .and. k .eq. 1) then
!       print *,'i,j,k',i,j,k
!       print *,'tabs(i,j,k)',tabs(i,j,k)
!       print *,'tabs1(i,j,k)',tabs1
!       print *,'lstarn(i,j,k)',lstarn
!       print *,'q(i,j,k)',q(i,j,k)
!       print *,'qsat(i,j,k)',qsat(i,j,k)
!       print *,'qp(i,j,k)',qp(i,j,k)
!       print *,'fff(i,j,k)',fff
!       print *,'dfff(i,j,k)',dfff
!       print *,'dqsat(i,j,k)',dqsat
!       print *,'dtabs(i,j,k)',dtabs
      endif

              end do   

                   qsat(i,j,k) = qsat(i,j,k) + dqsat * dtabs
              qn(i,j,k) = max(0.,q(i,j,k)-qsat(i,j,k))
              omegan(i,j,k) = om
              omegap(i,j,k) = omp
              if(om.lt.1.) then
                 kmin = min(kmin,k)
                 kmax = max(kmax,k)
              end if
            else

              qn(i,j,k) = 0.
              omegap(i,j,k) = max(0.,min(1.,ap*(tabs(i,j,k)-tprmin)))

            endif

            tabs(i,j,k) = tabs1
            omegag(i,j,k) = max(0.,min(1.,ag*(tabs(i,j,k)-tgrmin)))
            qp(i,j,k) = max(0.,qp(i,j,k)) ! just in case

          end do
         end do
        end do

凌波微步！

联系过多次模式作者，因为用的是模式中比较简单低级的微物理方案，可简单方案代码容易看。后来他不回我了。