LCOV - code coverage report
Current view: | top level - physics/cam - gw_convect.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 81 | 84 | 96.4 %
Date: | 2024-12-17 22:39:59 |  | Functions: | 2 | 4 | 50.0 %

              Line data    Source code
           1             : module gw_convect
           2             :
           3             : !
           4             : ! This module handles gravity waves from convection, and was extracted from
           5             : ! gw_drag in May 2013.
           6             : !
           7             :
           8             : use gw_utils, only: r8
           9             :
          10             : implicit none
          11             : private
          12             : save
          13             :
          14             : public :: BeresSourceDesc
          15             : public :: gw_beres_src
          16             :
          17             : type :: BeresSourceDesc
          18             :    ! Whether wind speeds are shifted to be relative to storm cells.
          19             :    logical :: storm_shift
          20             :    ! Index for level where wind speed is used as the source speed.
          21             :    integer :: k
          22             :    ! Heating depths below this value [m] will be ignored.
          23             :    real(r8) :: min_hdepth
          24             :    ! Table bounds, for convenience. (Could be inferred from shape(mfcc).)
          25             :    integer :: maxh
          26             :    integer :: maxuh
          27             :    ! Heating depths [m].
          28             :    real(r8), allocatable :: hd(:)
          29             :    ! Table of source spectra.
          30             :    real(r8), allocatable :: mfcc(:,:,:)
          31             : end type BeresSourceDesc
          32             :
          33             : contains
          34             :
          35             : !==========================================================================
          36             :
          37     1489176 : subroutine gw_beres_src(ncol, band, desc, u, v, &
          38     2978352 :      netdt, zm, src_level, tend_level, tau, ubm, ubi, xv, yv, &
          39     1489176 :      c, hdepth, maxq0)
          40             : !-----------------------------------------------------------------------
          41             : ! Driver for multiple gravity wave drag parameterization.
          42             : !
          43             : ! The parameterization is assumed to operate only where water vapor
          44             : ! concentrations are negligible in determining the density.
          45             : !
          46             : ! Beres, J.H., M.J. Alexander, and J.R. Holton, 2004: "A method of
          47             : ! specifying the gravity wave spectrum above convection based on latent
          48             : ! heating properties and background wind". J. Atmos. Sci., Vol 61, No. 3,
          49             : ! pp. 324-337.
          50             : !
          51             : !-----------------------------------------------------------------------
          52             :   use gw_utils, only: get_unit_vector, dot_2d, midpoint_interp
          53             :   use gw_common, only: GWBand, pver, qbo_hdepth_scaling
          54             :
          55             : !------------------------------Arguments--------------------------------
          56             :   ! Column dimension.
          57             :   integer, intent(in) :: ncol
          58             :
          59             :   ! Wavelengths triggered by convection.
          60             :   type(GWBand), intent(in) :: band
          61             :
          62             :   ! Settings for convection type (e.g. deep vs shallow).
          63             :   type(BeresSourceDesc), intent(in) :: desc
          64             :
          65             :   ! Midpoint zonal/meridional winds.
          66             :   real(r8), intent(in) :: u(ncol,pver), v(ncol,pver)
          67             :   ! Heating rate due to convection.
          68             :   real(r8), intent(in) :: netdt(:,:)
          69             :   ! Midpoint altitudes.
          70             :   real(r8), intent(in) :: zm(ncol,pver)
          71             :
          72             :   ! Indices of top gravity wave source level and lowest level where wind
          73             :   ! tendencies are allowed.
          74             :   integer, intent(out) :: src_level(ncol)
          75             :   integer, intent(out) :: tend_level(ncol)
          76             :
          77             :   ! Wave Reynolds stress.
          78             :   real(r8), intent(out) :: tau(ncol,-band%ngwv:band%ngwv,pver+1)
          79             :   ! Projection of wind at midpoints and interfaces.
          80             :   real(r8), intent(out) :: ubm(ncol,pver), ubi(ncol,pver+1)
          81             :   ! Unit vectors of source wind (zonal and meridional components).
          82             :   real(r8), intent(out) :: xv(ncol), yv(ncol)
          83             :   ! Phase speeds.
          84             :   real(r8), intent(out) :: c(ncol,-band%ngwv:band%ngwv)
          85             :
          86             :   ! Heating depth [m] and maximum heating in each column.
          87             :   real(r8), intent(out) :: hdepth(ncol), maxq0(ncol)
          88             :
          89             : !---------------------------Local Storage-------------------------------
          90             :   ! Column and level indices.
          91             :   integer :: i, k
          92             :
          93             :   ! Zonal/meridional wind at roughly the level where the convection occurs.
          94     2978352 :   real(r8) :: uconv(ncol), vconv(ncol)
          95             :
          96             :   ! Maximum heating rate.
          97     2978352 :   real(r8) :: q0(ncol)
          98             :
          99             :   ! Bottom/top heating range index.
         100     2978352 :   integer  :: boti(ncol), topi(ncol)
         101             :   ! Index for looking up heating depth dimension in the table.
         102     2978352 :   integer  :: hd_idx(ncol)
         103             :   ! Mean wind in heating region.
         104     2978352 :   real(r8) :: uh(ncol)
         105             :   ! Min/max wavenumber for critical level filtering.
         106     2978352 :   integer :: Umini(ncol), Umaxi(ncol)
         107             :   ! Source level tau for a column.
         108     2978352 :   real(r8) :: tau0(-band%ngwv:band%ngwv)
         109             :   ! Speed of convective cells relative to storm.
         110     2978352 :   real(r8) :: CS(ncol)
         111             :   ! Index to shift spectra relative to ground.
         112             :   integer :: shift
         113             :
         114             :   ! Heating rate conversion factor.
         115             :   real(r8), parameter :: CF = 20._r8
         116             :   ! Averaging length.
         117             :   real(r8), parameter :: AL = 1.0e5_r8
         118             :
         119             :   !----------------------------------------------------------------------
         120             :   ! Initialize tau array
         121             :   !----------------------------------------------------------------------
         122             :
         123 >15207*10^7 :   tau = 0.0_r8
         124    24865776 :   hdepth = 0.0_r8
         125    24865776 :   q0 = 0.0_r8
         126    98285616 :   tau0 = 0.0_r8
         127             :
         128             :   !------------------------------------------------------------------------
         129             :   ! Determine wind and unit vectors approximately at the source level, then
         130             :   ! project winds.
         131             :   !------------------------------------------------------------------------
         132             :
         133             :   ! Source wind speed and direction.
         134    24865776 :   uconv = u(:,desc%k)
         135    24865776 :   vconv = v(:,desc%k)
         136             :
         137             :   ! Get the unit vector components and magnitude at the source level.
         138     1489176 :   call get_unit_vector(uconv, vconv, xv, yv, ubi(:,desc%k+1))
         139             :
         140             :   ! Project the local wind at midpoints onto the source wind.
         141   139982544 :   do k = 1, pver
         142   139982544 :      ubm(:,k) = dot_2d(u(:,k), v(:,k), xv, yv)
         143             :   end do
         144             :
         145             :   ! Compute the interface wind projection by averaging the midpoint winds.
         146             :   ! Use the top level wind at the top interface.
         147    24865776 :   ubi(:,1) = ubm(:,1)
         148             :
         149     1489176 :   ubi(:,2:pver) = midpoint_interp(ubm)
         150             :
         151             :   !-----------------------------------------------------------------------
         152             :   ! Calculate heating depth.
         153             :   !
         154             :   ! Heating depth is defined as the first height range from the bottom in
         155             :   ! which heating rate is continuously positive.
         156             :   !-----------------------------------------------------------------------
         157             :
         158             :   ! First find the indices for the top and bottom of the heating range.
         159    24865776 :   boti = 0
         160    24865776 :   topi = 0
         161    84061881 :   do k = pver, 1, -1
         162  1403715195 :      do i = 1, ncol
         163  1403715195 :         if (boti(i) == 0) then
         164             :            ! Detect if we are outside the maximum range (where z = 20 km).
         165  1115303792 :            if (zm(i,k) >= 20000._r8) then
         166    18921569 :               boti(i) = k
         167    18921569 :               topi(i) = k
         168             :            else
         169             :               ! First spot where heating rate is positive.
         170  1096382223 :               if (netdt(i,k) > 0.0_r8) boti(i) = k
         171             :            end if
         172   204349522 :         else if (topi(i) == 0) then
         173             :            ! Detect if we are outside the maximum range (z = 20 km).
         174    55645860 :            if (zm(i,k) >= 20000._r8) then
         175           0 :               topi(i) = k
         176             :            else
         177             :               ! First spot where heating rate is no longer positive.
         178    55645860 :               if (.not. (netdt(i,k) > 0.0_r8)) topi(i) = k
         179             :            end if
         180             :         end if
         181             :      end do
         182             :      ! When all done, exit.
         183   146927255 :      if (all(topi /= 0)) exit
         184             :   end do
         185             :
         186             :   ! Heating depth in m.
         187    48242376 :   hdepth = [ ( (zm(i,topi(i))-zm(i,boti(i))), i = 1, ncol ) ]
         188             :
         189             :   ! J. Richter: this is an effective reduction of the GW phase speeds (needed to drive the QBO)
         190    24865776 :   hdepth = hdepth*qbo_hdepth_scaling
         191             :
         192     1489176 :   hd_idx = index_of_nearest(hdepth, desc%hd)
         193             :
         194             :   ! hd_idx=0 signals that a heating depth is too shallow, i.e. that it is
         195             :   ! either not big enough for the lowest table entry, or it is below the
         196             :   ! minimum allowed for this convection type.
         197             :   ! Values above the max in the table still get the highest value, though.
         198    24865776 :   where (hdepth < max(desc%min_hdepth, desc%hd(1))) hd_idx = 0
         199             :
         200             :   ! Maximum heating rate.
         201    79895164 :   do k = minval(topi), maxval(boti)
         202   531288868 :      where (k >= topi .and. k <= boti)
         203    31652788 :         q0 = max(q0, netdt(:,k))
         204             :      end where
         205             :   end do
         206             :
         207             :   !output max heating rate in K/day
         208    24865776 :   maxq0 = q0*24._r8*3600._r8
         209             :
         210             :   ! Multipy by conversion factor
         211    24865776 :   q0 = q0 * CF
         212             :
         213     1489176 :   if (desc%storm_shift) then
         214             :
         215             :      ! Find the cell speed where the storm speed is > 10 m/s.
         216             :      ! Storm speed is taken to be the source wind speed.
         217    24865776 :      CS = sign(max(abs(ubm(:,desc%k))-10._r8, 0._r8), ubm(:,desc%k))
         218             :
         219             :      ! Average wind in heating region, relative to storm cells.
         220    24865776 :      uh = 0._r8
         221    79895164 :      do k = minval(topi), maxval(boti)
         222   531288868 :         where (k >= topi .and. k <= boti)
         223    31652788 :            uh = uh + ubm(:,k)/(boti-topi+1)
         224             :         end where
         225             :      end do
         226             :
         227    24865776 :      uh = uh - CS
         228             :
         229             :   else
         230             :
         231             :      ! For shallow convection, wind is relative to ground, and "heating
         232             :      ! region" wind is just the source level wind.
         233           0 :      uh = ubm(:,desc%k)
         234             :
         235             :   end if
         236             :
         237             :   ! Limit uh to table range.
         238    24865776 :   uh = min(uh, real(desc%maxuh, r8))
         239    24865776 :   uh = max(uh, -real(desc%maxuh, r8))
         240             :
         241             :   ! Speeds for critical level filtering.
         242    24865776 :   Umini =  band%ngwv
         243    24865776 :   Umaxi = -band%ngwv
         244    79895164 :   do k = minval(topi), maxval(boti)
         245  1559235464 :      where (k >= topi .and. k <= boti)
         246    31652788 :         Umini = min(Umini, nint(ubm(:,k)/band%dc))
         247             :         Umaxi = max(Umaxi, nint(ubm(:,k)/band%dc))
         248             :      end where
         249             :   end do
         250             :
         251    24865776 :   Umini = max(Umini, -band%ngwv)
         252    24865776 :   Umaxi = min(Umaxi, band%ngwv)
         253             :
         254             :   !-----------------------------------------------------------------------
         255             :   ! Gravity wave sources
         256             :   !-----------------------------------------------------------------------
         257             :   ! Start loop over all columns.
         258             :   !-----------------------------------------------------------------------
         259    24865776 :   do i=1,ncol
         260             :
         261             :      !---------------------------------------------------------------------
         262             :      ! Look up spectrum only if the heating depth is large enough, else set
         263             :      ! tau0 = 0.
         264             :      !---------------------------------------------------------------------
         265             :
         266    24865776 :      if (hd_idx(i) > 0) then
         267             :
         268             :         !------------------------------------------------------------------
         269             :         ! Look up the spectrum using depth and uh.
         270             :         !------------------------------------------------------------------
         271             :
         272    67079100 :         tau0 = desc%mfcc(hd_idx(i),nint(uh(i)),:)
         273             :
         274     1016350 :         if (desc%storm_shift) then
         275             :            ! For deep convection, the wind was relative to storm cells, so
         276             :            ! shift the spectrum so that it is now relative to the ground.
         277     1016350 :            shift = -nint(CS(i)/band%dc)
         278    67079100 :            tau0 = eoshift(tau0, shift)
         279             :         end if
         280             :
         281             :         ! Adjust magnitude.
         282    67079100 :         tau0 = tau0*q0(i)*q0(i)/AL
         283             :
         284             :         ! Adjust for critical level filtering.
         285     6346713 :         tau0(Umini(i):Umaxi(i)) = 0.0_r8
         286             :
         287    67079100 :         tau(i,:,topi(i)+1) = tau0
         288             :
         289             :      end if ! heating depth above min and not at the pole
         290             :
         291             :   enddo
         292             :
         293             :   !-----------------------------------------------------------------------
         294             :   ! End loop over all columns.
         295             :   !-----------------------------------------------------------------------
         296             :
         297             :   ! Output the source level.
         298    24865776 :   src_level = topi
         299    24865776 :   tend_level = topi
         300             :
         301             :   ! Set phase speeds; just use reference speeds.
         302     1489176 :   c = spread(band%cref, 1, ncol)
         303             :
         304     1489176 : end subroutine gw_beres_src
         305             :
         306             : ! Short routine to get the indices of a set of values rounded to their
         307             : ! nearest points on a grid.
         308     1489176 : function index_of_nearest(x, grid) result(idx)
         309             :   real(r8), intent(in) :: x(:)
         310             :   real(r8), intent(in) :: grid(:)
         311             :
         312             :   integer :: idx(size(x))
         313             :
         314     2978352 :   real(r8) :: interfaces(size(grid)-1)
         315             :   integer :: i, n
         316             :
         317     1489176 :   n = size(grid)
         318    31272696 :   interfaces = (grid(:n-1) + grid(2:))/2._r8
         319             :
         320    24865776 :   idx = 1
         321    29783520 :   do i = 1, n-1
         322   473938920 :      where (x > interfaces(i)) idx = i + 1
         323             :   end do
         324             :
         325     1489176 : end function index_of_nearest
         326             :
         327           0 : end module gw_convect

Generated by: LCOV version 1.14