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

              Line data    Source code
           1             : module zm_conv_momtran
           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             :   public zm_conv_momtran_run      ! convective momentum transport
          10             :
          11             :
          12             : contains
          13             :
          14             : !===============================================================================
          15             : !> \section arg_table_zm_conv_momtran_run Argument Table
          16             : !! \htmlinclude zm_conv_momtran_run.html
          17             : !!
          18     1495368 : subroutine zm_conv_momtran_run(ncol, pver, pverp, &
          19     1495368 :                     domomtran,windu, windv,num_winds, mu, md, &
          20             :                     momcu, momcd, &
          21     1495368 :                     du, eu, ed, dp, dsubcld , &
          22     1495368 :                     jt, mx, ideep , il1g, il2g, &
          23     1495368 :                     nstep, windu_tend, windv_tend, pguallu, pguallv, pgdallu, pgdallv, &
          24     2990736 :                     icwuu, icwuv, icwdu, icwdv, dt, seten)
          25             : !-----------------------------------------------------------------------
          26             : !
          27             : ! Purpose:
          28             : ! Convective transport of momentum
          29             : !
          30             : ! Mixing ratios may be with respect to either dry or moist air
          31             : !
          32             : ! Method:
          33             : ! Based on the convtran subroutine by P. Rasch
          34             : ! <Also include any applicable external references.>
          35             : !
          36             : ! Author: J. Richter and P. Rasch
          37             : !
          38             : !-----------------------------------------------------------------------
          39             : ! CACNOTE - use CCPP constituents object
          40             :    use constituents,    only: cnst_get_type_byind
          41             :
          42             :    implicit none
          43             : !-----------------------------------------------------------------------
          44             : !
          45             : ! Input arguments
          46             : !
          47             :    integer, intent(in) :: ncol                  ! number of atmospheric columns
          48             :    integer, intent(in) :: num_winds             ! number of wind directions
          49             :    integer, intent(in) :: pver, pverp
          50             :    logical, intent(in) :: domomtran(:)      ! flag for doing convective transport    (num_winds)
          51             :    real(kind_phys), intent(in) :: windu(:,:)  ! U Wind array                                    (ncol,pver)
          52             :    real(kind_phys), intent(in) :: windv(:,:)  ! V Wind array                                    (ncol,pver)
          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) :: momcu
          56             :    real(kind_phys), intent(in) :: momcd
          57             :    real(kind_phys), intent(in) :: du(:,:)       ! Mass detraining from updraft              (ncol,pver)
          58             :    real(kind_phys), intent(in) :: eu(:,:)       ! Mass entraining from updraft              (ncol,pver)
          59             :    real(kind_phys), intent(in) :: ed(:,:)       ! Mass entraining from downdraft            (ncol,pver)
          60             :    real(kind_phys), intent(in) :: dp(:,:)       ! Delta pressure between interfaces         (ncol,pver)
          61             :    real(kind_phys), intent(in) :: dsubcld(:)       ! Delta pressure from cloud base to sfc  (ncol)
          62             :    real(kind_phys), intent(in) :: dt            ! time step in seconds
          63             :
          64             :    integer, intent(in) :: jt(:)         ! Index of cloud top for each column         (ncol)
          65             :    integer, intent(in) :: mx(:)         ! Index of cloud top for each column         (ncol)
          66             :    integer, intent(in) :: ideep(:)      ! Gathering array                            (ncol)
          67             :    integer, intent(in) :: il1g              ! Gathered min lon indices over which to operate
          68             :    integer, intent(in) :: il2g              ! Gathered max lon indices over which to operate
          69             :    integer, intent(in) :: nstep             ! Time step index
          70             :
          71             :
          72             :
          73             : ! input/output
          74             :
          75             :    real(kind_phys), intent(out) :: windu_tend(:,:)  ! U wind tendency
          76             :    real(kind_phys), intent(out) :: windv_tend(:,:)  ! V wind tendency
          77             :
          78             : !--------------------------Local Variables------------------------------
          79             :
          80             :    integer i                 ! Work index
          81             :    integer k                 ! Work index
          82             :    integer kbm               ! Highest altitude index of cloud base
          83             :    integer kk                ! Work index
          84             :    integer kkp1              ! Work index
          85             :    integer kkm1              ! Work index
          86             :    integer km1               ! Work index
          87             :    integer kp1               ! Work index
          88             :    integer ktm               ! Highest altitude index of cloud top
          89             :    integer m                 ! Work index
          90             :    integer ii                 ! Work index
          91             :
          92             :    real(kind_phys) cabv                 ! Mix ratio of constituent above
          93             :    real(kind_phys) cbel                 ! Mix ratio of constituent below
          94             :    real(kind_phys) cdifr                ! Normalized diff between cabv and cbel
          95     2990736 :    real(kind_phys) chat(ncol,pver)     ! Mix ratio in env at interfaces
          96     2990736 :    real(kind_phys) cond(ncol,pver)     ! Mix ratio in downdraft at interfaces
          97     2990736 :    real(kind_phys) const(ncol,pver)    ! Gathered wind array
          98     2990736 :    real(kind_phys) conu(ncol,pver)     ! Mix ratio in updraft at interfaces
          99     2990736 :    real(kind_phys) dcondt(ncol,pver)   ! Gathered tend array
         100             :    real(kind_phys) mbsth                ! Threshold for mass fluxes
         101             :    real(kind_phys) mupdudp              ! A work variable
         102             :    real(kind_phys) minc                 ! A work variable
         103             :    real(kind_phys) maxc                 ! A work variable
         104             :    real(kind_phys) fluxin               ! A work variable
         105             :    real(kind_phys) fluxout              ! A work variable
         106             :    real(kind_phys) netflux              ! A work variable
         107             :
         108             :
         109             :    real(kind_phys) sum                  ! sum
         110             :    real(kind_phys) sum2                  ! sum2
         111             :
         112     2990736 :    real(kind_phys) mududp(ncol,pver) ! working variable
         113     2990736 :    real(kind_phys) mddudp(ncol,pver)     ! working variable
         114             :
         115     2990736 :    real(kind_phys) pgu(ncol,pver)      ! Pressure gradient term for updraft
         116     2990736 :    real(kind_phys) pgd(ncol,pver)      ! Pressure gradient term for downdraft
         117             :
         118             :    real(kind_phys),intent(out) ::  pguallu(:,:)      ! Apparent force from  updraft PG on U winds  ! (ncol,pver)
         119             :    real(kind_phys),intent(out) ::  pguallv(:,:)      ! Apparent force from  updraft PG on V winds  ! (ncol,pver)
         120             :    real(kind_phys),intent(out) ::  pgdallu(:,:)      ! Apparent force from  downdraft PG on U winds! (ncol,pver)
         121             :    real(kind_phys),intent(out) ::  pgdallv(:,:)      ! Apparent force from  downdraft PG on V winds! (ncol,pver)
         122             :
         123             :    real(kind_phys),intent(out) ::  icwuu(:,:)      ! In-cloud U winds in updraft           ! (ncol,pver)
         124             :    real(kind_phys),intent(out) ::  icwuv(:,:)      ! In-cloud V winds in updraft           ! (ncol,pver)
         125             :    real(kind_phys),intent(out) ::  icwdu(:,:)      ! In-cloud U winds in downdraft         ! (ncol,pver)
         126             :    real(kind_phys),intent(out) ::  icwdv(:,:)      ! In-cloud V winds in downdraft         ! (ncol,pver)
         127             :
         128             :    real(kind_phys),intent(out) ::  seten(:,:) ! Dry static energy tendency                ! (ncol,pver)
         129     2990736 :    real(kind_phys)                 gseten(ncol,pver) ! Gathered dry static energy tendency
         130             :
         131     2990736 :    real(kind_phys) :: winds(ncol,pver,num_winds)       ! combined winds array
         132     2990736 :    real(kind_phys) :: wind_tends(ncol,pver,num_winds)  ! combined tendency array
         133     2990736 :    real(kind_phys) :: pguall(ncol,pver,num_winds)      ! Combined apparent force from  updraft PG on U winds
         134     2990736 :    real(kind_phys) :: pgdall(ncol,pver,num_winds)      ! Combined apparent force from  downdraft PG on U winds
         135     2990736 :    real(kind_phys) :: icwu(ncol,pver,num_winds)        ! Combined In-cloud winds in updraft
         136     2990736 :    real(kind_phys) :: icwd(ncol,pver,num_winds)        ! Combined In-cloud winds in downdraft
         137             :
         138     2990736 :    real(kind_phys)  mflux(ncol,pverp,num_winds)   ! Gathered momentum flux
         139             :
         140     2990736 :    real(kind_phys)  wind0(ncol,pver,num_winds)       !  gathered  wind before time step
         141     2990736 :    real(kind_phys)  windf(ncol,pver,num_winds)       !  gathered  wind after time step
         142             :    real(kind_phys) fkeb, fket, ketend_cons, ketend, utop, ubot, vtop, vbot, gset2
         143             :
         144             :
         145             : !-----------------------------------------------------------------------
         146             : !
         147             : ! Combine winds in single array
         148   652189104 :    winds(:,:,1) = windu(:,:)
         149   650693736 :    winds(:,:,2) = windv(:,:)
         150             :
         151             : ! Initialize outgoing fields
         152  1302882840 :    pguall(:,:,:)     = 0.0_kind_phys
         153  1302882840 :    pgdall(:,:,:)     = 0.0_kind_phys
         154             : ! Initialize in-cloud winds to environmental wind
         155  1302882840 :    icwu(:ncol,:,:)       = winds(:ncol,:,:)
         156  1302882840 :    icwd(:ncol,:,:)       = winds(:ncol,:,:)
         157             :
         158             : ! Initialize momentum flux and  final winds
         159  1352821176 :    mflux(:,:,:)       = 0.0_kind_phys
         160  1302882840 :    wind0(:,:,:)         = 0.0_kind_phys
         161  1302882840 :    windf(:,:,:)         = 0.0_kind_phys
         162             :
         163             : ! Initialize dry static energy
         164             :
         165   650693736 :    seten(:,:)         = 0.0_kind_phys
         166   650693736 :    gseten(:,:)         = 0.0_kind_phys
         167             :
         168             : ! mbsth is the threshold below which we treat the mass fluxes as zero (in mb/s)
         169     1495368 :    mbsth = 1.e-15_kind_phys
         170             :
         171             : ! Find the highest level top and bottom levels of convection
         172     1495368 :    ktm = pver
         173     1495368 :    kbm = pver
         174    13077368 :    do i = il1g, il2g
         175    11582000 :       ktm = min(ktm,jt(i))
         176    13077368 :       kbm = min(kbm,mx(i))
         177             :    end do
         178             :
         179             : ! Loop ever each wind component
         180     4486104 :    do m = 1, num_winds                    !start at m = 1 to transport momentum
         181     4486104 :       if (domomtran(m)) then
         182             :
         183             : ! Gather up the winds and set tend to zero
         184    80749872 :          do k = 1,pver
         185   683013872 :             do i =il1g,il2g
         186   602264000 :                const(i,k) = winds(ideep(i),k,m)
         187   680023136 :                 wind0(i,k,m) = const(i,k)
         188             :             end do
         189             :          end do
         190             :
         191             :
         192             : ! From now on work only with gathered data
         193             :
         194             : ! Interpolate winds to interfaces
         195             :
         196    80749872 :          do k = 1,pver
         197    77759136 :             km1 = max(1,k-1)
         198   683013872 :             do i = il1g, il2g
         199             :
         200             :                ! use arithmetic mean
         201   602264000 :                chat(i,k) = 0.5_kind_phys* (const(i,k)+const(i,km1))
         202             :
         203             : ! Provisional up and down draft values
         204   602264000 :                conu(i,k) = chat(i,k)
         205   602264000 :                cond(i,k) = chat(i,k)
         206             :
         207             : !              provisional tends
         208   680023136 :                dcondt(i,k) = 0._kind_phys
         209             :
         210             :             end do
         211             :          end do
         212             :
         213             :
         214             : !
         215             : ! Pressure Perturbation Term
         216             : !
         217             :
         218             :       !Top boundary:  assume mu is zero
         219             :
         220     2990736 :          k=1
         221    29145472 :          pgu(:il2g,k) = 0.0_kind_phys
         222    26154736 :          pgd(:il2g,k) = 0.0_kind_phys
         223             :
         224    74768400 :          do k=2,pver-1
         225    71777664 :             km1 = max(1,k-1)
         226    71777664 :             kp1 = min(pver,k+1)
         227   630704400 :             do i = il1g,il2g
         228             :
         229             :                !interior points
         230             :
         231  1667808000 :                mududp(i,k) =  ( mu(i,k) * (const(i,k)- const(i,km1))/dp(i,km1) &
         232  2223744000 :                            +  mu(i,kp1) * (const(i,kp1) - const(i,k))/dp(i,k))
         233             :
         234   555936000 :                pgu(i,k) = - momcu * 0.5_kind_phys * mududp(i,k)
         235             :
         236             :
         237   555936000 :                mddudp(i,k) =  ( md(i,k) * (const(i,k)- const(i,km1))/dp(i,km1) &
         238   555936000 :                            +  md(i,kp1) * (const(i,kp1) - const(i,k))/dp(i,k))
         239             :
         240   627713664 :                pgd(i,k) = - momcd * 0.5_kind_phys * mddudp(i,k)
         241             :
         242             :
         243             :             end do
         244             :          end do
         245             :
         246             :        ! bottom boundary
         247     2990736 :        k = pver
         248     2990736 :        km1 = max(1,k-1)
         249    26154736 :        do i=il1g,il2g
         250             :
         251    23164000 :           mududp(i,k) =   mu(i,k) * (const(i,k)- const(i,km1))/dp(i,km1)
         252    23164000 :           pgu(i,k) = - momcu *  mududp(i,k)
         253             :
         254    23164000 :           mddudp(i,k) =   md(i,k) * (const(i,k)- const(i,km1))/dp(i,km1)
         255             :
         256    26154736 :           pgd(i,k) = - momcd * mddudp(i,k)
         257             :
         258             :        end do
         259             :
         260             :
         261             : !
         262             : ! In-cloud velocity calculations
         263             : !
         264             :
         265             : ! Do levels adjacent to top and bottom
         266    26154736 :          k = 2
         267    26154736 :          km1 = 1
         268    26154736 :          kk = pver
         269    26154736 :          kkm1 = max(1,kk-1)
         270    26154736 :          do i = il1g,il2g
         271    23164000 :             mupdudp = mu(i,kk) + du(i,kk)*dp(i,kk)
         272    23164000 :             if (mupdudp > mbsth) then
         273             :
         274    22223084 :                conu(i,kk) = (+eu(i,kk)*const(i,kk)*dp(i,kk)+pgu(i,kk)*dp(i,kk))/mupdudp
         275             :             endif
         276    26154736 :             if (md(i,k) < -mbsth) then
         277           0 :                cond(i,k) =  (-ed(i,km1)*const(i,km1)*dp(i,km1))-pgd(i,km1)*dp(i,km1)/md(i,k)
         278             :             endif
         279             :
         280             :
         281             :          end do
         282             :
         283             :
         284             :
         285             : ! Updraft from bottom to top
         286    77759136 :          do kk = pver-1,1,-1
         287    74768400 :             kkm1 = max(1,kk-1)
         288    74768400 :             kkp1 = min(pver,kk+1)
         289   656859136 :             do i = il1g,il2g
         290   579100000 :                mupdudp = mu(i,kk) + du(i,kk)*dp(i,kk)
         291   653868400 :                if (mupdudp > mbsth) then
         292             :
         293   236049804 :                   conu(i,kk) = (  mu(i,kkp1)*conu(i,kkp1)+eu(i,kk)* &
         294   236049804 :                                   const(i,kk)*dp(i,kk)+pgu(i,kk)*dp(i,kk))/mupdudp
         295             :                endif
         296             :             end do
         297             :
         298             :          end do
         299             :
         300             :
         301             : ! Downdraft from top to bottom
         302    74768400 :          do k = 3,pver
         303    71777664 :             km1 = max(1,k-1)
         304   630704400 :             do i = il1g,il2g
         305   627713664 :                if (md(i,k) < -mbsth) then
         306             :
         307   196768404 :                   cond(i,k) =  (  md(i,km1)*cond(i,km1)-ed(i,km1)*const(i,km1) &
         308   196768404 :                                   *dp(i,km1)-pgd(i,km1)*dp(i,km1) )/md(i,k)
         309             :
         310             :                endif
         311             :             end do
         312             :          end do
         313             :
         314             :
         315             :          sum = 0._kind_phys
         316             :          sum2 = 0._kind_phys
         317             :
         318             :
         319    18397672 :          do k = ktm,pver
         320    15406936 :             km1 = max(1,k-1)
         321    15406936 :             kp1 = min(pver,k+1)
         322   186633220 :             do i = il1g,il2g
         323   168235548 :                ii = ideep(i)
         324             :
         325             : ! version 1 hard to check for roundoff errors
         326   168235548 :                dcondt(i,k) =  &
         327   168235548 :                            +(mu(i,kp1)* (conu(i,kp1)-chat(i,kp1)) &
         328   168235548 :                            -mu(i,k)*   (conu(i,k)-chat(i,k))      &
         329   168235548 :                            +md(i,kp1)* (cond(i,kp1)-chat(i,kp1)) &
         330   168235548 :                            -md(i,k)*   (cond(i,k)-chat(i,k)) &
         331   520113580 :                           )/dp(i,k)
         332             :
         333             :             end do
         334             :          end do
         335             :
         336             :   ! dcont for bottom layer
         337             :           !
         338     6246608 :           do k = kbm,pver
         339    28967488 :              km1 = max(1,k-1)
         340    31958224 :              do i = il1g,il2g
         341    28967488 :                 if (k == mx(i)) then
         342             :
         343             :                    ! version 1
         344    23164000 :                    dcondt(i,k) = (1._kind_phys/dp(i,k))*   &
         345    23164000 :                         (-mu(i,k)*(conu(i,k)-chat(i,k)) &
         346    23164000 :                         -md(i,k)*(cond(i,k)-chat(i,k)) &
         347    46328000 :                         )
         348             :                 end if
         349             :              end do
         350             :           end do
         351             :
         352             : ! Initialize to zero everywhere, then scatter tendency back to full array
         353  1301387472 :          wind_tends(:,:,m) = 0._kind_phys
         354             :
         355    80749872 :          do k = 1,pver
         356   683013872 :             do i = il1g,il2g
         357   602264000 :                ii = ideep(i)
         358   602264000 :                wind_tends(ii,k,m) = dcondt(i,k)
         359             :     ! Output apparent force on the mean flow from pressure gradient
         360   602264000 :                pguall(ii,k,m) = -pgu(i,k)
         361   602264000 :                pgdall(ii,k,m) = -pgd(i,k)
         362   602264000 :                icwu(ii,k,m)   =  conu(i,k)
         363   680023136 :                icwd(ii,k,m)   =  cond(i,k)
         364             :             end do
         365             :          end do
         366             :
         367             :           ! Calculate momentum flux in units of mb*m/s2
         368             :
         369    18397672 :           do k = ktm,pver
         370   186633220 :              do i = il1g,il2g
         371   168235548 :                 ii = ideep(i)
         372   168235548 :                 mflux(i,k,m) = &
         373   168235548 :                      -mu(i,k)*   (conu(i,k)-chat(i,k))      &
         374   351878032 :                      -md(i,k)*   (cond(i,k)-chat(i,k))
         375             :              end do
         376             :           end do
         377             :
         378             :
         379             :           ! Calculate winds at the end of the time step
         380             :
         381    18397672 :           do k = ktm,pver
         382   186633220 :              do i = il1g,il2g
         383   168235548 :                 ii = ideep(i)
         384   168235548 :                 km1 = max(1,k-1)
         385   168235548 :                 kp1 = k+1
         386   183642484 :                 windf(i,k,m) = const(i,k)    -   (mflux(i,kp1,m) - mflux(i,k,m)) * dt /dp(i,k)
         387             :
         388             :              end do
         389             :           end do
         390             :
         391             :        end if      ! for domomtran
         392             :    end do
         393             :
         394             :  ! Need to add an energy fix to account for the dissipation of kinetic energy
         395             :     ! Formulation follows from Boville and Bretherton (2003)
         396             :     ! formulation by PJR
         397             :
         398     9198836 :     do k = ktm,pver
         399     7703468 :        km1 = max(1,k-1)
         400     7703468 :        kp1 = min(pver,k+1)
         401    93316610 :        do i = il1g,il2g
         402             :
         403    84117774 :           ii = ideep(i)
         404             :
         405             :           ! calculate the KE fluxes at top and bot of layer
         406             :           ! based on a discrete approximation to b&b eq(35) F_KE = u*F_u + v*F_v at interface
         407    84117774 :           utop = (wind0(i,k,1)+wind0(i,km1,1))/2._kind_phys
         408    84117774 :           vtop = (wind0(i,k,2)+wind0(i,km1,2))/2._kind_phys
         409    84117774 :           ubot = (wind0(i,kp1,1)+wind0(i,k,1))/2._kind_phys
         410    84117774 :           vbot = (wind0(i,kp1,2)+wind0(i,k,2))/2._kind_phys
         411    84117774 :           fket = utop*mflux(i,k,1)   + vtop*mflux(i,k,2)    ! top of layer
         412    84117774 :           fkeb = ubot*mflux(i,k+1,1) + vbot*mflux(i,k+1,2)  ! bot of layer
         413             :
         414             :           ! divergence of these fluxes should give a conservative redistribution of KE
         415    84117774 :           ketend_cons = (fket-fkeb)/dp(i,k)
         416             :
         417             :           ! tendency in kinetic energy resulting from the momentum transport
         418    84117774 :           ketend = ((windf(i,k,1)**2 + windf(i,k,2)**2) - (wind0(i,k,1)**2 + wind0(i,k,2)**2))/dt
         419             :
         420             :           ! the difference should be the dissipation
         421    84117774 :           gset2 = ketend_cons - ketend
         422    91821242 :           gseten(i,k) = gset2
         423             :
         424             :        end do
         425             :
         426             :     end do
         427             :
         428             :     ! Scatter dry static energy to full array
         429    40374936 :     do k = 1,pver
         430   341506936 :        do i = il1g,il2g
         431   301132000 :           ii = ideep(i)
         432   340011568 :           seten(ii,k) = gseten(i,k)
         433             :
         434             :        end do
         435             :     end do
         436             :
         437             : ! Split out the wind tendencies
         438   652189104 :    windu_tend(:,:) = wind_tends(:,:,1)
         439   652189104 :    windv_tend(:,:) = wind_tends(:,:,2)
         440             :
         441   652189104 :    pguallu(:,:)     = pguall(:,:,1)
         442   652189104 :    pguallv(:,:)     = pguall(:,:,2)
         443   652189104 :    pgdallu(:,:)     = pgdall(:,:,1)
         444   652189104 :    pgdallv(:,:)     = pgdall(:,:,2)
         445   652189104 :    icwuu(:ncol,:)       = icwu(:,:,1)
         446   652189104 :    icwuv(:ncol,:)       = icwu(:,:,2)
         447   652189104 :    icwdu(:ncol,:)       = icwd(:,:,1)
         448   652189104 :    icwdv(:ncol,:)       = icwd(:,:,2)
         449             :
         450     1495368 :    return
         451             : end subroutine zm_conv_momtran_run
         452             :
         453             :
         454             : end module zm_conv_momtran

Generated by: LCOV version 1.14