zm_conv.f90

来自「CCSM Research Tools: Community Atmospher」· F90 代码 · 共 1,872 行 · 第 1/5 页

         cme(ideep(i),k) = cmeg(i,k)
!++bee
         cmfdqr(ideep(i),k) = cmfdqrg(i,k)
!--bee
         zdu(ideep(i),k) = du2(i,k)
         mcon(ideep(i),k) = mc(i,k)
!**BAB not this tendency
!!$         qtg(ideep(i),k) = dqdt(i,k)
!**BAB report back heating in J/kg/s
!!$         ttg(ideep(i),k) = dsdt(i,k)
         ttg(ideep(i),k) = dsdt(i,k)*cpres
!            utg(ideep(i),k) = dudt(i,k)
!            vtg(ideep(i),k) = dvdt(i,k)
         dlf(ideep(i),k) = dlg(i,k)
         pflx(ideep(i),k) = pflxg(i,k)
         ql(ideep(i),k) = qlg(i,k)
      end do
   end do
!
   do i = 1,lengath
      jctop(ideep(i)) = jt(i)
!++bee
      jcbot(ideep(i)) = maxg(i)
!--bee
      pflx(ideep(i),pverp) = pflxg(i,pverp)
      psevap = psevap + totevp(i)
      psrain = psrain + totpcp(i)
   end do
!
! convert back to specific humidity from mixing ratio.
! take into account any moisture added to ensure positiveness
! of specific humidity at start of routine.
!
   do k = msg + 1,pver
      do i = 1,ncol
         q(i,k) = q(i,k) - max((qmin-qh(i,k,1)),0._r8)
      end do
   end do
   do k = pver,msg + 1,-1
      do i = 1,ncol
         sumq(i) = sumq(i) - dpp(i,k)* (q(i,k)-qh(i,k,1))
         qtg(i,k) = (q(i,k)-qh(i,k,1)) / (2.*delt)
!
! account for the detraining cloud water in the precip
!
         sumq(i) = sumq(i) - dpp(i,k)*dlf(i,k)*2*delt
         pcpck(i,k) = max(0._r8,sumq(i))
      end do
   end do
!
! obtain final precipitation rate.
!
   do i = 1,ncol
!         llo1 = ts(i) .ge. tfreez
!
! here pcpr and pcps are in units of kg/m^2, ie. precip per
! time step
!
!         pcpr(i) = cvmgt(pcpdh(i)*max(sumq(i),0._r8),0.,llo1)
!         pcps(i) = cvmgt(0.,pcpdh(i)*max(sumq(i),0._r8),llo1)
      precip = pcpdh(i)*max(sumq(i),0._r8)
      if (ts(i) >= tfreez) then
         pcpr(i) = precip
         pcps(i) = 0.
      else
         pcpr(i) = 0.
         pcps(i) = precip
      end if
   end do

!
! accumulate precipitation, the 1000. is the density of water, so
! paprc and paprs are now in units of meters.
!
   do i = 1,ncol
      paprc(i) = paprc(i) + (pcpr(i)+pcps(i))/1000.
      paprs(i) = paprs(i) + pcps(i)/1000.
   end do
!
! convert precipitation to m/s, ie, precip rate.
!
   do i = 1,ncol
      pcpr(i) = pcpr(i)/ (2.*delt)/1000.
      pcps(i) = pcps(i)/ (2.*delt)/1000.
      pcpc(i) = pcpr(i) + pcps(i)
      psheat = psheat + (pcps(i)+pcpr(i))*rl
   end do
   do k = msg + 1,pver
      do i = 1,ncol
         pcpck(i,k) = pcpdh(i)*pcpck(i,k)/ (2.*delt)
      end do
   end do
!
! calculate conservation of quantities.
!

!       if(lat(1).eq.jlatpr)then
!        do l=msg+1,pver
!        do i=1,lengath
!          sumdq(i) = sumdq(i) + 2.*delt*(rl/cpres)*dpp(ideep(i),l)*
!     1                          dqdt(i,l)
!          sumdt(i) = sumdt(i) + 2.*delt*dpp(ideep(i),l)*dsdt(i,l)
!        end do
!        end do
!
!        write(6,*)'sumdq,sumdt,sumde in convection subroutine########'
!        do i=1,lengath
!          sumde(i) = sumdt(i) + sumdq(i)
!          write(6, 901) sumdq(i), sumdt(i),sumde(i), i, ideep(i)
!        end do
!c
!        write(6,*)'sumdq,sumdt,sumde ... all points'
!      do i=1,ncol
!         sumdq(i) = 0.0
!         sumdt(i) = 0.0
!      end do
!c
!      do l=msg+1,pver
!      do i=1,ncol
!        deltaq(i,l) = qh(i,l,1) - deltaq(i,l)
!        deltat(i,l) = t (i,l) - deltat(i,l)
!        sumdq(i) = sumdq(i) + (rl/cpres)*dpp(i,l)*deltaq(i,l)
!        sumdt(i) = sumdt(i) + dpp(i,l)*deltat(i,l)
!      end do
!      end do
!      do i=1,ncol
!        sumde(i) = sumdt(i) + sumdq(i)
!      end do
!        write(6, 902) (i,sumdq(i),sumdt(i),sumde(i),i=1,ncol)
!
!  901   format(1x,3e20.12, i10, i10)
!  902   format(1x,i10, 3e20.12)
!
!      endif

   return
end subroutine zm_convr

!===============================================================================
  subroutine zm_conv_evap(state, ptend, pflx, precc, cldo, deltat, evappct)

!-----------------------------------------------------------------------
! compute tendencies due to evaporation of rain from ZM scheme
!-----------------------------------------------------------------------

    use dycore,         only: dycore_is
    use physics_types,  only: physics_state, physics_ptend
    use wv_saturation,  only: aqsat

#include <guang.h>
!------------------------------Arguments--------------------------------
    real(r8), intent(in   ) :: cldo(pcols,pver)   ! old cloud fraction
    real(r8), intent(in   ) :: deltat             ! time step
    type(physics_state), intent(in)    :: state   ! Physics state variables
    type(physics_ptend), intent(inout) :: ptend   ! indivdual parameterization tendencies

    real(r8), intent(inout) :: precc(pcols)       ! Convective-scale preciptn rate
    real(r8), intent(inout) :: pflx (pcols,pverp) ! Conv rain flux thru out btm of lev
    real(r8), intent(out)   :: evappct(pcols)     ! Convective-scale preciptn rate
!
!---------------------------Local storage-------------------------------

    real(r8) :: dpovrg                 ! Pressure depth over gravity
    real(r8) :: envevap                ! Evaporation rate from precitation
    real(r8) :: est (pcols,pver)       ! Saturation vapor pressure
    real(r8) :: ke                     ! Tunable evaporation efficiency
    real(r8) :: pflxtmp(pcols,pver)    ! Conv rain flux thru out btm of lev (temp)
    real(r8) :: qsat(pcols,pver)       ! Specific humidity at Saturation
    real(r8) :: rh                     ! Relative humidity
    real(r8) :: sumflx(pcols)          ! flux integral 

    integer :: i,k                     ! longitude,level indices

!-----------------------------------------------------------------------
!
! Evaporate some of the precip directly into the environment (Sundqvist)
! pflx = kg/m^2/s
!
    call aqsat (state%t    ,state%pmid  ,est    ,qsat    ,pcols   , &
         state%ncol ,pver  ,1       ,pver    )

    sumflx(:) = 0.

    if ( dycore_is ('LR') ) then
!====== finite-volume dynamics =============================
       ke = 7.5E-6     ! larger value to compensate for cloud fraction effect 
       sumflx = 0.
       pflx(:,1) = 0.

       do k=1,pver
          do i=1,state%ncol
             if(precc(i) > 0.001*sumflx(i) .and. state%t(i,k) > tfreez) then
                dpovrg  =  state%pdel(i,k)*rgrav
                pflxtmp(i,k) = max(pflx(i,k)-sumflx(i),0._r8)
                envevap = max(ke*(1.-state%q(i,k,1)/qsat(i,k))*sqrt(pflxtmp(i,k)),0._r8)

! Sundqvist's Cloud fraction effect

                envevap = envevap * (1. - cldo(i,k))
                envevap = max(min(envevap, 0.75*(qsat(i,k)-state%q(i,k,1))/deltat),0._r8)
                envevap = min(pflxtmp(i,k)/dpovrg,envevap)

! The factor 0.999 is to prevent rounding-level negative precc

                envevap = min(envevap,0.999*(1000.*precc(i)-sumflx(i))/dpovrg)
                envevap = min(envevap,max(0._r8,0.999*(pflx(i,k)-sumflx(i))/dpovrg))
                pflx(i,k) = pflx(i,k) - envevap*dpovrg - sumflx(i)
                sumflx(i) = sumflx(i) + envevap*dpovrg
                ptend%q(i,k,1) =  envevap
                ptend%s(i,k)   = -envevap*rl
             else
                ptend%q(i,k,1) =  0.
                ptend%s(i,k)   =  0.
              endif
          end do
       end do
    else
!====== EUL/SLD dynamics =============================
       ke = 3.0E-6
       do k=1,pver
          do i=1,state%ncol
             dpovrg  =  state%pdel(i,k)*rgrav
             rh      = state%q(i,k,1)/qsat(i,k)
             pflxtmp(i,k)= max(pflx(i,k)-sumflx(i),0._r8)

! Determine evap amount based on rh can't be negative                                

             envevap   = max(ke*(1.0 - rh)*sqrt(pflxtmp(i,k)),0._r8)                        

! Don't let envevap supersaturate layer (approx)                                     

             envevap   = max(min(envevap, (qsat(i,k)-state%q(i,k,1))/deltat),0._r8)          

! Can't evap more than what we have at this level                                    

             envevap   = min(pflxtmp(i,k)/dpovrg,envevap)                                

! Can't evap more than what we are told is in bottom level                           
! precc and pflx(,pver) are supposed to be identical but they differ by              
! roundoff, precc is used for this check since that is what is written out           

             envevap  = min((precc(i)*1.E3-sumflx(i))/dpovrg,envevap)                   
             pflx(i,k) = pflx(i,k) - envevap*dpovrg - sumflx(i)
             sumflx(i) = sumflx(i)+envevap*dpovrg
             ptend%q(i,k,1) =  envevap
             ptend%s(i,k)   = -envevap*rl
          end do
       end do
    endif
!
! adjust precc by the amount of precip evaporated
!
    evappct(:)=0.
    do i = 1, state%ncol
       precc(i)=precc(i)-(sumflx(i)/1.E3)
       if (precc(i) > 0) evappct(i)=(sumflx(i)/1.E3)/precc(i)
    end do

  end subroutine zm_conv_evap

subroutine buoyan(lchnk   ,ncol    , &
                  q       ,t       ,p       ,z       ,pf      , &
                  tp      ,qstp    ,tl      ,rl      ,cape    , &
                  pblt    ,lcl     ,lel     ,lon     ,mx      , &
                  rd      ,grav    ,cp      ,msg     ,nstep   , &
                  tpert   )
!----------------------------------------------------------------------- 
! 
! Purpose: 
! <Say what the routine does> 
! 
! Method: 
! <Describe the algorithm(s) used in the routine.> 
! <Also include any applicable external references.> 
! 
! Author:
! This is contributed code not fully standardized by the CCM core group.
! The documentation has been enhanced to the degree that we are able.
! Reviewed:          P. Rasch, April 1996
! 
!-----------------------------------------------------------------------
   implicit none
#include <guang.h>
!-----------------------------------------------------------------------
!
! input arguments
!
   integer, intent(in) :: lchnk                 ! chunk identifier
   integer, intent(in) :: ncol                  ! number of atmospheric columns

   real(r8), intent(in) :: q(pcols,pver)        ! spec. humidity
   real(r8), intent(in) :: t(pcols,pver)        ! temperature
   real(r8), intent(in) :: p(pcols,pver)        ! pressure
   real(r8), intent(in) :: z(pcols,pver)        ! height
   real(r8), intent(in) :: pf(pcols,pver+1)     ! pressure at interfaces
   real(r8), intent(in) :: pblt(pcols)          ! index of pbl depth
   real(r8), intent(in) :: tpert(pcols)         ! perturbation temperature by pbl processes

!
! output arguments
!
   real(r8), intent(out) :: tp(pcols,pver)       ! parcel temperature
   real(r8), intent(out) :: qstp(pcols,pver)     ! saturation mixing ratio of parcel
   real(r8), intent(out) :: tl(pcols)            ! parcel temperature at lcl
   real(r8), intent(out) :: cape(pcols)          ! convective aval. pot. energy.
   integer lcl(pcols)        !
   integer lel(pcols)        !
   integer lon(pcols)        ! level of onset of deep convection
   integer mx(pcols)         ! level of max moist static energy
!
!--------------------------Local Variables------------------------------
!
   real(r8) capeten(pcols,5)     ! provisional value of cape
   real(r8) tv(pcols,pver)       !
   real(r8) tpv(pcols,pver)      !
   real(r8) buoy(pcols,pver)

   real(r8) a1(pcols)
   real(r8) a2(pcols)
   real(r8) estp(pcols)
   real(r8) pl(pcols)
   real(r8) plexp(pcols)
   real(r8) hmax(pcols)
   real(r8) hmn(pcols)
   real(r8) y(pcols)

   logical plge600(pcols)
   integer knt(pcols)
   integer lelten(pcols,5)

   real(r8) cp
   real(r8) e
   real(r8) grav

   integer i
   integer k
   integer msg
   integer n
   integer nstep

   real(r8) rd
   real(r8) rl
#ifdef PERGRO
   real(r8) rhd
#endif
!
!-----------------------------------------------------------------------
!
   do n = 1,5
      do i = 1,ncol
         lelten(i,n) = pver
         capeten(i,n) = 0.
      end do
   end do
!
   do i = 1,ncol
      lon(i) = pver
      knt(i) = 0
      lel(i) = pver
      mx(i) = lon(i)
      cape(i) = 0.
      hmax(i) = 0.
   end do

   tp(:ncol,:) = t(:ncol,:)
   qstp(:ncol,:) = q(:ncol,:)
!
! set "launching" level(mx) to be at maximum moist static energy.
! search for this level stops at planetary boundary layer top.
!
#ifdef PERGRO
   do k = pver,msg + 1,-1
      do i = 1,ncol
         hmn(i) = cp*t(i,k) + grav*z(i,k) + rl*q(i,k)
!