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: | 98 | 107 | 91.6 %
Date: | 2024-12-17 17:57:11 |  | 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     5720188 :    do i = il1g, il2g
         122     4224820 :       ktm = min(ktm,jt(i))
         123     5720188 :       kbm = min(kbm,mx(i))
         124             :    end do
         125             :
         126             : ! Loop ever each constituent
         127    16449048 :    do m = 2, ncnst
         128    16449048 :       if (doconvtran(m)) then
         129             :
         130    14953680 :          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  1420599600 :             do k = 1,pver
         141  5334728520 :                do i =il1g,il2g
         142  3929082600 :                   dptmp(i,k) = dp(i,k)
         143  3929082600 :                   dutmp(i,k) = du(i,k)
         144  3929082600 :                   eutmp(i,k) = eu(i,k)
         145  5319774840 :                   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  1405645920 :          do k = 1,pver
         153  5334728520 :             do i =il1g,il2g
         154  3929082600 :                const(i,k) = q(ideep(i),k,m)
         155  5319774840 :                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  1405645920 :          do k = 1,pver
         163  1390692240 :             km1 = max(1,k-1)
         164  5334728520 :             do i = il1g, il2g
         165  3929082600 :                minc = min(const(i,km1),const(i,k))
         166  3929082600 :                maxc = max(const(i,km1),const(i,k))
         167  3929082600 :                if (minc < 0) then
         168             :                   cdifr = 0._kind_phys
         169             :                else
         170  3929082600 :                   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  3929082600 :                if (cdifr > 1.E-6_kind_phys) then
         176  2011955509 :                   cabv = max(const(i,km1),maxc*1.e-12_kind_phys)
         177  2011955509 :                   cbel = max(const(i,k),maxc*1.e-12_kind_phys)
         178  2011955509 :                   chat(i,k) = log(cabv/cbel)/(cabv-cbel)*cabv*cbel
         179             :
         180             :                else             ! Small diff, so just arithmetic mean
         181  1917127091 :                   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  3929082600 :                conu(i,k) = chat(i,k)
         186  3929082600 :                cond(i,k) = chat(i,k)
         187             :
         188             : !              provisional tends
         189  5319774840 :                dcondt(i,k) = 0._kind_phys
         190             :
         191             :             end do
         192             :          end do
         193             :
         194             : ! Do levels adjacent to top and bottom
         195    57201880 :          k = 2
         196    57201880 :          km1 = 1
         197    57201880 :          kk = pver
         198    57201880 :          do i = il1g,il2g
         199    42248200 :             mupdudp = mu(i,kk) + dutmp(i,kk)*dptmp(i,kk)
         200    42248200 :             if (mupdudp > mbsth) then
         201           0 :                conu(i,kk) = (+eutmp(i,kk)*fisg(i,kk)*const(i,kk)*dptmp(i,kk))/mupdudp
         202             :             endif
         203    57201880 :             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  1390692240 :          do kk = pver-1,1,-1
         210  1375738560 :             kkp1 = min(pver,kk+1)
         211  5277526640 :             do i = il1g,il2g
         212  3886834400 :                mupdudp = mu(i,kk) + dutmp(i,kk)*dptmp(i,kk)
         213  5262572960 :                if (mupdudp > mbsth) then
         214   675426140 :                   conu(i,kk) = (  mu(i,kkp1)*conu(i,kkp1)+eutmp(i,kk)*fisg(i,kk)* &
         215  1350852280 :                                   const(i,kk)*dptmp(i,kk) )/mupdudp
         216             :                endif
         217             :             end do
         218             :          end do
         219             :
         220             : ! Downdraft from top to bottom
         221  1375738560 :          do k = 3,pver
         222  1360784880 :             km1 = max(1,k-1)
         223  5220324760 :             do i = il1g,il2g
         224  5205371080 :                if (md(i,k) < -mbsth) then
         225  1045511660 :                   cond(i,k) =  (  md(i,km1)*cond(i,km1)-edtmp(i,km1)*fisg(i,km1)*const(i,km1) &
         226  1045511660 :                                   *dptmp(i,km1) )/md(i,k)
         227             :                endif
         228             :             end do
         229             :          end do
         230             :
         231             :
         232   153173160 :          do k = ktm,pver
         233   138219480 :             km1 = max(1,k-1)
         234   138219480 :             kp1 = min(pver,k+1)
         235  1166219970 :             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  4052187240 :                fluxin =  mu(i,kp1)*conu(i,kp1)+ mu(i,k)*min(chat(i,k),const(i,km1)) &
         241  5065234050 :                          -(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  1013046810 :                          -(md(i,kp1)*cond(i,kp1) + md(i,k)*min(chat(i,k),const(i,k)))
         244             :
         245  1013046810 :                netflux = fluxin - fluxout
         246  1013046810 :                if (abs(netflux) < max(fluxin,fluxout)*1.e-12_kind_phys) then
         247   146297763 :                   netflux = 0._kind_phys
         248             :                endif
         249  1151266290 :                dcondt(i,k) = netflux/dptmp(i,k)
         250             :             end do
         251             :          end do
         252             : ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
         253             : !
         254    59572740 :          do k = kbm,pver
         255    44619060 :             km1 = max(1,k-1)
         256   328335810 :             do i = il1g,il2g
         257   313382130 :                if (k == mx(i)) then
         258             :
         259    42248200 :                   fluxin =  mu(i,k)*min(chat(i,k),const(i,km1)) - md(i,k)*cond(i,k)
         260    42248200 :                   fluxout = mu(i,k)*conu(i,k) - md(i,k)*min(chat(i,k),const(i,k))
         261             :
         262    42248200 :                   netflux = fluxin - fluxout
         263    42248200 :                   if (abs(netflux) < max(fluxin,fluxout)*1.e-12_kind_phys) then
         264           0 :                      netflux = 0._kind_phys
         265             :                   endif
         266    42248200 :                   dcondt(i,k) = netflux/dptmp(i,k)
         267   226514870 :                else if (k > mx(i)) then
         268   189933590 :                   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 23236279920 :          dqdt(:,:,m) = 0._kind_phys
         275  1405645920 :          do k = 1,pver
         276  1390692240 :             kp1 = min(pver,k+1)
         277  5334728520 :             do i = il1g,il2g
         278  5319774840 :                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