LCOV - code coverage report
Current view: | top level - atmos_phys/schemes/zhang_mcfarlane - zm_conv_convtran.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 99 | 107 | 92.5 %
Date: | 2025-01-13 21:54:50 |  | Functions: | 1 | 1 | 100.0 %

              Line data    Source code
           1             : module zm_conv_convtran
           2             :
           3             :   use ccpp_kinds, only:  kind_phys
           4             :
           5             :   implicit none
           6             :
           7             :   save
           8             :   private                         ! Make default type private to the module
           9             : !
          10             : ! PUBLIC: interfaces
          11             : !
          12             :   public zm_conv_convtran_run     ! convective transport
          13             :
          14             :
          15             : contains
          16             :
          17             : !===============================================================================
          18             : !> \section arg_table_zm_conv_convtran_run Argument Table
          19             : !! \htmlinclude zm_conv_convtran_run.html
          20             : !!
          21     1495368 : subroutine zm_conv_convtran_run(ncol, pver, &
          22     1495368 :                     doconvtran,q       ,ncnst   ,mu      ,md      , &
          23     1495368 :                     du      ,eu      ,ed      ,dp      ,dsubcld , &
          24     2990736 :                     jt      ,mx      ,ideep   ,il1g    ,il2g    , &
          25     4486104 :                     nstep   ,fracis  ,dqdt    ,dpdry)
          26             : !-----------------------------------------------------------------------
          27             : !
          28             : ! Purpose:
          29             : ! Convective transport of trace species
          30             : !
          31             : ! Mixing ratios may be with respect to either dry or moist air
          32             : !
          33             : ! Method:
          34             : ! <Describe the algorithm(s) used in the routine.>
          35             : ! <Also include any applicable external references.>
          36             : !
          37             : ! Author: P. Rasch
          38             : !
          39             : !-----------------------------------------------------------------------
          40             : ! CACNOTE - replace with CCPP constituents
          41             :    use constituents,    only: cnst_get_type_byind
          42             :
          43             :    implicit none
          44             : !-----------------------------------------------------------------------
          45             : !
          46             : ! Input arguments
          47             : !
          48             :    integer, intent(in) :: ncol
          49             :    integer, intent(in) :: pver
          50             :    integer, intent(in) :: ncnst          ! number of tracers to transport
          51             :    logical, intent(in) :: doconvtran(:)  ! flag for doing convective transport      (ncnst)
          52             :    real(kind_phys), intent(in) :: q(:,:,:)      ! Tracer array including moisture          (ncol,pver,ncnst)
          53             :    real(kind_phys), intent(in) :: mu(:,:)       ! Mass flux up                             (ncol,pver)
          54             :    real(kind_phys), intent(in) :: md(:,:)       ! Mass flux down                           (ncol,pver)
          55             :    real(kind_phys), intent(in) :: du(:,:)       ! Mass detraining from updraft             (ncol,pver)
          56             :    real(kind_phys), intent(in) :: eu(:,:)       ! Mass entraining from updraft             (ncol,pver)
          57             :    real(kind_phys), intent(in) :: ed(:,:)       ! Mass entraining from downdraft           (ncol,pver)
          58             :    real(kind_phys), intent(in) :: dp(:,:)       ! Delta pressure between interfaces        (ncol,pver)
          59             :    real(kind_phys), intent(in) :: dsubcld(:)    ! Delta pressure from cloud base to sfc    (ncol)
          60             :    real(kind_phys), intent(in) :: fracis(:,:,:) ! fraction of tracer that is insoluble     (ncol,pver,ncnst)
          61             :
          62             :    integer, intent(in) :: jt(:)          ! Index of cloud top for each column       (ncol)
          63             :    integer, intent(in) :: mx(:)          ! Index of cloud top for each column       (ncol)
          64             :    integer, intent(in) :: ideep(:)       ! Gathering array                          (ncol)
          65             :    integer, intent(in) :: il1g           ! Gathered min lon indices over which to operate
          66             :    integer, intent(in) :: il2g           ! Gathered max lon indices over which to operate
          67             :    integer, intent(in) :: nstep          ! Time step index
          68             :
          69             :    real(kind_phys), intent(in) :: dpdry(:,:)    ! Delta pressure between interfaces        (ncol,pver)
          70             :
          71             : ! input/output
          72             :
          73             :    real(kind_phys), intent(out) :: dqdt(:,:,:)  ! Tracer tendency array  (ncol,pver,ncnst)
          74             :
          75             : !--------------------------Local Variables------------------------------
          76             :
          77             :    integer i                 ! Work index
          78             :    integer k                 ! Work index
          79             :    integer kbm               ! Highest altitude index of cloud base
          80             :    integer kk                ! Work index
          81             :    integer kkp1              ! Work index
          82             :    integer km1               ! Work index
          83             :    integer kp1               ! Work index
          84             :    integer ktm               ! Highest altitude index of cloud top
          85             :    integer m                 ! Work index
          86             :
          87             :    real(kind_phys) cabv                 ! Mix ratio of constituent above
          88             :    real(kind_phys) cbel                 ! Mix ratio of constituent below
          89             :    real(kind_phys) cdifr                ! Normalized diff between cabv and cbel
          90     2990736 :    real(kind_phys) chat(ncol,pver)     ! Mix ratio in env at interfaces
          91     2990736 :    real(kind_phys) cond(ncol,pver)     ! Mix ratio in downdraft at interfaces
          92     2990736 :    real(kind_phys) const(ncol,pver)    ! Gathered tracer array
          93     2990736 :    real(kind_phys) fisg(ncol,pver)     ! gathered insoluble fraction of tracer
          94     2990736 :    real(kind_phys) conu(ncol,pver)     ! Mix ratio in updraft at interfaces
          95     2990736 :    real(kind_phys) dcondt(ncol,pver)   ! Gathered tend array
          96             :    real(kind_phys) small                ! A small number
          97             :    real(kind_phys) mbsth                ! Threshold for mass fluxes
          98             :    real(kind_phys) mupdudp              ! A work variable
          99             :    real(kind_phys) minc                 ! A work variable
         100             :    real(kind_phys) maxc                 ! A work variable
         101             :    real(kind_phys) fluxin               ! A work variable
         102             :    real(kind_phys) fluxout              ! A work variable
         103             :    real(kind_phys) netflux              ! A work variable
         104             :
         105     2990736 :    real(kind_phys) dutmp(ncol,pver)       ! Mass detraining from updraft
         106     2990736 :    real(kind_phys) eutmp(ncol,pver)       ! Mass entraining from updraft
         107     2990736 :    real(kind_phys) edtmp(ncol,pver)       ! Mass entraining from downdraft
         108     2990736 :    real(kind_phys) dptmp(ncol,pver)    ! Delta pressure between interfaces
         109             :    real(kind_phys) total(ncol)
         110             :    real(kind_phys) negadt,qtmp
         111             :
         112             : !-----------------------------------------------------------------------
         113             : !
         114     1495368 :    small = 1.e-36_kind_phys
         115             : ! mbsth is the threshold below which we treat the mass fluxes as zero (in mb/s)
         116     1495368 :    mbsth = 1.e-15_kind_phys
         117             :
         118             : ! Find the highest level top and bottom levels of convection
         119     1495368 :    ktm = pver
         120     1495368 :    kbm = pver
         121    13077368 :    do i = il1g, il2g
         122    11582000 :       ktm = min(ktm,jt(i))
         123    13077368 :       kbm = min(kbm,mx(i))
         124             :    end do
         125             :
         126             : ! Loop ever each constituent
         127     4486104 :    do m = 2, ncnst
         128     4486104 :       if (doconvtran(m)) then
         129             :
         130     2990736 :          if (cnst_get_type_byind(m).eq.'dry') then
         131           0 :             do k = 1,pver
         132           0 :                do i =il1g,il2g
         133           0 :                   dptmp(i,k) = dpdry(i,k)
         134           0 :                   dutmp(i,k) = du(i,k)*dp(i,k)/dpdry(i,k)
         135           0 :                   eutmp(i,k) = eu(i,k)*dp(i,k)/dpdry(i,k)
         136           0 :                   edtmp(i,k) = ed(i,k)*dp(i,k)/dpdry(i,k)
         137             :                end do
         138             :             end do
         139             :          else
         140    83740608 :             do k = 1,pver
         141   683013872 :                do i =il1g,il2g
         142   602264000 :                   dptmp(i,k) = dp(i,k)
         143   602264000 :                   dutmp(i,k) = du(i,k)
         144   602264000 :                   eutmp(i,k) = eu(i,k)
         145   680023136 :                   edtmp(i,k) = ed(i,k)
         146             :                end do
         147             :             end do
         148             :          endif
         149             : !        dptmp = dp
         150             :
         151             : ! Gather up the constituent and set tend to zero
         152    80749872 :          do k = 1,pver
         153   683013872 :             do i =il1g,il2g
         154   602264000 :                const(i,k) = q(ideep(i),k,m)
         155   680023136 :                fisg(i,k) = fracis(ideep(i),k,m)
         156             :             end do
         157             :          end do
         158             :
         159             : ! From now on work only with gathered data
         160             :
         161             : ! Interpolate environment tracer values to interfaces
         162    80749872 :          do k = 1,pver
         163    77759136 :             km1 = max(1,k-1)
         164   683013872 :             do i = il1g, il2g
         165   602264000 :                minc = min(const(i,km1),const(i,k))
         166   602264000 :                maxc = max(const(i,km1),const(i,k))
         167   602264000 :                if (minc < 0) then
         168             :                   cdifr = 0._kind_phys
         169             :                else
         170   602264000 :                   cdifr = abs(const(i,k)-const(i,km1))/max(maxc,small)
         171             :                endif
         172             :
         173             : ! If the two layers differ significantly use a geometric averaging
         174             : ! procedure
         175   602264000 :                if (cdifr > 1.E-6_kind_phys) then
         176   414024028 :                   cabv = max(const(i,km1),maxc*1.e-12_kind_phys)
         177   414024028 :                   cbel = max(const(i,k),maxc*1.e-12_kind_phys)
         178   414024028 :                   chat(i,k) = log(cabv/cbel)/(cabv-cbel)*cabv*cbel
         179             :
         180             :                else             ! Small diff, so just arithmetic mean
         181   188239972 :                   chat(i,k) = 0.5_kind_phys* (const(i,k)+const(i,km1))
         182             :                end if
         183             :
         184             : ! Provisional up and down draft values
         185   602264000 :                conu(i,k) = chat(i,k)
         186   602264000 :                cond(i,k) = chat(i,k)
         187             :
         188             : !              provisional tends
         189   680023136 :                dcondt(i,k) = 0._kind_phys
         190             :
         191             :             end do
         192             :          end do
         193             :
         194             : ! Do levels adjacent to top and bottom
         195    26154736 :          k = 2
         196    26154736 :          km1 = 1
         197    26154736 :          kk = pver
         198    26154736 :          do i = il1g,il2g
         199    23164000 :             mupdudp = mu(i,kk) + dutmp(i,kk)*dptmp(i,kk)
         200    23164000 :             if (mupdudp > mbsth) then
         201    22223084 :                conu(i,kk) = (+eutmp(i,kk)*fisg(i,kk)*const(i,kk)*dptmp(i,kk))/mupdudp
         202             :             endif
         203    26154736 :             if (md(i,k) < -mbsth) then
         204           0 :                cond(i,k) =  (-edtmp(i,km1)*fisg(i,km1)*const(i,km1)*dptmp(i,km1))/md(i,k)
         205             :             endif
         206             :          end do
         207             :
         208             : ! Updraft from bottom to top
         209    77759136 :          do kk = pver-1,1,-1
         210    74768400 :             kkp1 = min(pver,kk+1)
         211   656859136 :             do i = il1g,il2g
         212   579100000 :                mupdudp = mu(i,kk) + dutmp(i,kk)*dptmp(i,kk)
         213   653868400 :                if (mupdudp > mbsth) then
         214   118024902 :                   conu(i,kk) = (  mu(i,kkp1)*conu(i,kkp1)+eutmp(i,kk)*fisg(i,kk)* &
         215   236049804 :                                   const(i,kk)*dptmp(i,kk) )/mupdudp
         216             :                endif
         217             :             end do
         218             :          end do
         219             :
         220             : ! Downdraft from top to bottom
         221    74768400 :          do k = 3,pver
         222    71777664 :             km1 = max(1,k-1)
         223   630704400 :             do i = il1g,il2g
         224   627713664 :                if (md(i,k) < -mbsth) then
         225   196768404 :                   cond(i,k) =  (  md(i,km1)*cond(i,km1)-edtmp(i,km1)*fisg(i,km1)*const(i,km1) &
         226   196768404 :                                   *dptmp(i,km1) )/md(i,k)
         227             :                endif
         228             :             end do
         229             :          end do
         230             :
         231             :
         232    18397672 :          do k = ktm,pver
         233    15406936 :             km1 = max(1,k-1)
         234    15406936 :             kp1 = min(pver,k+1)
         235   186633220 :             do i = il1g,il2g
         236             :
         237             : ! limit fluxes outside convection to mass in appropriate layer
         238             : ! these limiters are probably only safe for positive definite quantitities
         239             : ! it assumes that mu and md already satify a courant number limit of 1
         240   672942192 :                fluxin =  mu(i,kp1)*conu(i,kp1)+ mu(i,k)*min(chat(i,k),const(i,km1)) &
         241   841177740 :                          -(md(i,k)  *cond(i,k) + md(i,kp1)*min(chat(i,kp1),const(i,kp1)))
         242             :                fluxout = mu(i,k)*conu(i,k) + mu(i,kp1)*min(chat(i,kp1),const(i,k)) &
         243   168235548 :                          -(md(i,kp1)*cond(i,kp1) + md(i,k)*min(chat(i,k),const(i,k)))
         244             :
         245   168235548 :                netflux = fluxin - fluxout
         246   168235548 :                if (abs(netflux) < max(fluxin,fluxout)*1.e-12_kind_phys) then
         247    12612158 :                   netflux = 0._kind_phys
         248             :                endif
         249   183642484 :                dcondt(i,k) = netflux/dptmp(i,k)
         250             :             end do
         251             :          end do
         252             : ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
         253             : !
         254     6246608 :          do k = kbm,pver
         255     3255872 :             km1 = max(1,k-1)
         256    31958224 :             do i = il1g,il2g
         257    28967488 :                if (k == mx(i)) then
         258             :
         259    23164000 :                   fluxin =  mu(i,k)*min(chat(i,k),const(i,km1)) - md(i,k)*cond(i,k)
         260    23164000 :                   fluxout = mu(i,k)*conu(i,k) - md(i,k)*min(chat(i,k),const(i,k))
         261             :
         262    23164000 :                   netflux = fluxin - fluxout
         263    23164000 :                   if (abs(netflux) < max(fluxin,fluxout)*1.e-12_kind_phys) then
         264           0 :                      netflux = 0._kind_phys
         265             :                   endif
         266    23164000 :                   dcondt(i,k) = netflux/dptmp(i,k)
         267     2547616 :                else if (k > mx(i)) then
         268      758318 :                   dcondt(i,k) = 0._kind_phys
         269             :                end if
         270             :             end do
         271             :          end do
         272             :
         273             : ! Initialize to zero everywhere, then scatter tendency back to full array
         274  1301387472 :          dqdt(:,:,m) = 0._kind_phys
         275    80749872 :          do k = 1,pver
         276    77759136 :             kp1 = min(pver,k+1)
         277   683013872 :             do i = il1g,il2g
         278   680023136 :                dqdt(ideep(i),k,m) = dcondt(i,k)
         279             :             end do
         280             :          end do
         281             :
         282             :       end if      ! for doconvtran
         283             :
         284             :    end do
         285             :
         286     1495368 :    return
         287             : end subroutine zm_conv_convtran_run
         288             :
         289             :
         290             : end module zm_conv_convtran

Generated by: LCOV version 1.14