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: | 104 | 107 | 97.2 %
Date: | 2024-12-17 22:39:59 |  | 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     2984544 : subroutine zm_conv_convtran_run(ncol, pver, &
          22     2984544 :                     doconvtran,q       ,ncnst   ,mu      ,md      , &
          23     2984544 :                     du      ,eu      ,ed      ,dp      ,dsubcld , &
          24     5969088 :                     jt      ,mx      ,ideep   ,il1g    ,il2g    , &
          25     8953632 :                     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     5969088 :    real(kind_phys) chat(ncol,pver)     ! Mix ratio in env at interfaces
          91     5969088 :    real(kind_phys) cond(ncol,pver)     ! Mix ratio in downdraft at interfaces
          92     5969088 :    real(kind_phys) const(ncol,pver)    ! Gathered tracer array
          93     5969088 :    real(kind_phys) fisg(ncol,pver)     ! gathered insoluble fraction of tracer
          94     5969088 :    real(kind_phys) conu(ncol,pver)     ! Mix ratio in updraft at interfaces
          95     5969088 :    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     5969088 :    real(kind_phys) dutmp(ncol,pver)       ! Mass detraining from updraft
         106     5969088 :    real(kind_phys) eutmp(ncol,pver)       ! Mass entraining from updraft
         107     5969088 :    real(kind_phys) edtmp(ncol,pver)       ! Mass entraining from downdraft
         108     5969088 :    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     2984544 :    small = 1.e-36_kind_phys
         115             : ! mbsth is the threshold below which we treat the mass fluxes as zero (in mb/s)
         116     2984544 :    mbsth = 1.e-15_kind_phys
         117             :
         118             : ! Find the highest level top and bottom levels of convection
         119     2984544 :    ktm = pver
         120     2984544 :    kbm = pver
         121    11988030 :    do i = il1g, il2g
         122     9003486 :       ktm = min(ktm,jt(i))
         123    11988030 :       kbm = min(kbm,mx(i))
         124             :    end do
         125             :
         126             : ! Loop ever each constituent
         127   122366304 :    do m = 2, ncnst
         128   122366304 :       if (doconvtran(m)) then
         129             :
         130    31334616 :          if (cnst_get_type_byind(m).eq.'dry') then
         131  1539807984 :             do k = 1,pver
         132  6135762336 :                do i =il1g,il2g
         133  4595954352 :                   dptmp(i,k) = dpdry(i,k)
         134  4595954352 :                   dutmp(i,k) = du(i,k)*dp(i,k)/dpdry(i,k)
         135  4595954352 :                   eutmp(i,k) = eu(i,k)*dp(i,k)/dpdry(i,k)
         136  6119381400 :                   edtmp(i,k) = ed(i,k)*dp(i,k)/dpdry(i,k)
         137             :                end do
         138             :             end do
         139             :          else
         140  1436980536 :             do k = 1,pver
         141  5600747580 :                do i =il1g,il2g
         142  4195101660 :                   dptmp(i,k) = dp(i,k)
         143  4195101660 :                   dutmp(i,k) = du(i,k)
         144  4195101660 :                   eutmp(i,k) = eu(i,k)
         145  5585793900 :                   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  2945453904 :          do k = 1,pver
         153 11736509916 :             do i =il1g,il2g
         154  8791056012 :                const(i,k) = q(ideep(i),k,m)
         155 11705175300 :                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  2945453904 :          do k = 1,pver
         163  2914119288 :             km1 = max(1,k-1)
         164 11736509916 :             do i = il1g, il2g
         165  8791056012 :                minc = min(const(i,km1),const(i,k))
         166  8791056012 :                maxc = max(const(i,km1),const(i,k))
         167  8791056012 :                if (minc < 0) then
         168             :                   cdifr = 0._kind_phys
         169             :                else
         170  8791056012 :                   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  8791056012 :                if (cdifr > 1.E-6_kind_phys) then
         176  6430731375 :                   cabv = max(const(i,km1),maxc*1.e-12_kind_phys)
         177  6430731375 :                   cbel = max(const(i,k),maxc*1.e-12_kind_phys)
         178  6430731375 :                   chat(i,k) = log(cabv/cbel)/(cabv-cbel)*cabv*cbel
         179             :
         180             :                else             ! Small diff, so just arithmetic mean
         181  2360324637 :                   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  8791056012 :                conu(i,k) = chat(i,k)
         186  8791056012 :                cond(i,k) = chat(i,k)
         187             :
         188             : !              provisional tends
         189 11705175300 :                dcondt(i,k) = 0._kind_phys
         190             :
         191             :             end do
         192             :          end do
         193             :
         194             : ! Do levels adjacent to top and bottom
         195   125862100 :          k = 2
         196   125862100 :          km1 = 1
         197   125862100 :          kk = pver
         198   125862100 :          do i = il1g,il2g
         199    94527484 :             mupdudp = mu(i,kk) + dutmp(i,kk)*dptmp(i,kk)
         200    94527484 :             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   125862100 :             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  2914119288 :          do kk = pver-1,1,-1
         210  2882784672 :             kkp1 = min(pver,kk+1)
         211 11610647816 :             do i = il1g,il2g
         212  8696528528 :                mupdudp = mu(i,kk) + dutmp(i,kk)*dptmp(i,kk)
         213 11579313200 :                if (mupdudp > mbsth) then
         214  1494239044 :                   conu(i,kk) = (  mu(i,kkp1)*conu(i,kkp1)+eutmp(i,kk)*fisg(i,kk)* &
         215  2988478088 :                                   const(i,kk)*dptmp(i,kk) )/mupdudp
         216             :                endif
         217             :             end do
         218             :          end do
         219             :
         220             : ! Downdraft from top to bottom
         221  2882784672 :          do k = 3,pver
         222  2851450056 :             km1 = max(1,k-1)
         223 11484785716 :             do i = il1g,il2g
         224 11453451100 :                if (md(i,k) < -mbsth) then
         225  2318136520 :                   cond(i,k) =  (  md(i,km1)*cond(i,km1)-edtmp(i,km1)*fisg(i,km1)*const(i,km1) &
         226  2318136520 :                                   *dptmp(i,km1) )/md(i,k)
         227             :                endif
         228             :             end do
         229             :          end do
         230             :
         231             :
         232   335152766 :          do k = ktm,pver
         233   303818150 :             km1 = max(1,k-1)
         234   303818150 :             kp1 = min(pver,k+1)
         235  2578978106 :             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  8975301360 :                fluxin =  mu(i,kp1)*conu(i,kp1)+ mu(i,k)*min(chat(i,k),const(i,km1)) &
         241 11219126700 :                          -(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  2243825340 :                          -(md(i,kp1)*cond(i,kp1) + md(i,k)*min(chat(i,k),const(i,k)))
         244             :
         245  2243825340 :                netflux = fluxin - fluxout
         246  2243825340 :                if (abs(netflux) < max(fluxin,fluxout)*1.e-12_kind_phys) then
         247   148281237 :                   netflux = 0._kind_phys
         248             :                endif
         249  2547643490 :                dcondt(i,k) = netflux/dptmp(i,k)
         250             :             end do
         251             :          end do
         252             : ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
         253             : !
         254   129131786 :          do k = kbm,pver
         255    97797170 :             km1 = max(1,k-1)
         256   734134289 :             do i = il1g,il2g
         257   702799673 :                if (k == mx(i)) then
         258             :
         259    94527484 :                   fluxin =  mu(i,k)*min(chat(i,k),const(i,km1)) - md(i,k)*cond(i,k)
         260    94527484 :                   fluxout = mu(i,k)*conu(i,k) - md(i,k)*min(chat(i,k),const(i,k))
         261             :
         262    94527484 :                   netflux = fluxin - fluxout
         263    94527484 :                   if (abs(netflux) < max(fluxin,fluxout)*1.e-12_kind_phys) then
         264           0 :                      netflux = 0._kind_phys
         265             :                   endif
         266    94527484 :                   dcondt(i,k) = netflux/dptmp(i,k)
         267   510475019 :                else if (k > mx(i)) then
         268   428131580 :                   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 48690349704 :          dqdt(:,:,m) = 0._kind_phys
         275  2945453904 :          do k = 1,pver
         276  2914119288 :             kp1 = min(pver,k+1)
         277 11736509916 :             do i = il1g,il2g
         278 11705175300 :                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     2984544 :    return
         287             : end subroutine zm_conv_convtran_run
         288             :
         289             :
         290             : end module zm_conv_convtran

Generated by: LCOV version 1.14