LCOV - code coverage report
Current view: | top level - physics/cam - gw_movmtn.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 105 | 126 | 83.3 %
Date: | 2024-12-17 17:57:11 |  | Functions: | 3 | 5 | 60.0 %

              Line data    Source code
           1             : module gw_movmtn
           2             :
           3             : !
           4             : ! This module parameterizes gravity waves generated by the obstacle effect produced by
           5             : ! boundary layer turbulence for convection.
           6             : !
           7             :
           8             : use gw_utils, only: r8
           9             :
          10             : implicit none
          11             : private
          12             : save
          13             :
          14             : public :: MovMtnSourceDesc
          15             : public :: gw_movmtn_src
          16             :
          17             : type :: MovMtnSourceDesc
          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 !-bounds of the lookup table heating depths
          26             :    integer :: maxuh ! bounds of the lookup table wind
          27             :    ! Heating depths [m].
          28             :    real(r8), allocatable :: hd(:), uh(:)
          29             :    ! Table of source spectra.
          30             :    real(r8), allocatable :: mfcc(:,:,:)  !is the lookup table f(depth, wind, phase speed)
          31             : end type MovMtnSourceDesc
          32             :
          33             : contains
          34             :
          35             : !==========================================================================
          36             :
          37     2978352 : subroutine gw_movmtn_src(ncol,lchnk, band, desc, u, v, &
          38     2978352 :      netdt, netdt_shcu, xpwp_shcu, &
          39     2978352 :      zm, alpha_gw_movmtn, src_level, tend_level, tau, ubm, ubi, xv, yv, &
          40     1489176 :      c, hdepth)
          41             : !-----------------------------------------------------------------------
          42             : ! Flexible driver for gravity wave source from obstacle effects produced
          43             : ! by boundary layer turbulence or deep convection
          44             : !-----------------------------------------------------------------------
          45             :   use gw_utils, only: get_unit_vector, dot_2d, midpoint_interp
          46             :   use gw_common, only: GWBand, pver, qbo_hdepth_scaling
          47             :   use cam_history, only: outfld
          48             :   use phys_control,   only: use_gw_movmtn_pbl
          49             :   use physconst,      only: rair, gravit
          50             : !------------------------------Arguments--------------------------------
          51             :   ! Column dimension.
          52             :   integer, intent(in) :: ncol , lchnk
          53             :
          54             :   ! Wavelengths triggered by convection.
          55             :   type(GWBand), intent(in) :: band
          56             :
          57             :   ! Settings for convection type (e.g. deep vs shallow).
          58             :   type(MovMtnSourceDesc), intent(in) :: desc
          59             :
          60             :   ! Midpoint zonal/meridional winds.
          61             :   real(r8), intent(in) :: u(ncol,pver), v(ncol,pver)
          62             :   ! Heating rate due to convection.
          63             :   real(r8), intent(in) :: netdt(:,:)  !from deep scheme
          64             :   ! Heating rate due to shallow convection and PBL turbulence.
          65             :   real(r8), intent(in) :: netdt_shcu(:,:)
          66             :   ! Higher order flux from ShCu/PBL.
          67             :   real(r8), intent(in) :: xpwp_shcu(ncol,pver+1)
          68             :   ! Midpoint altitudes.
          69             :   real(r8), intent(in) :: zm(ncol,pver)
          70             :   ! tunable parameter controlling proportion of PBL momentum flux emitted as GW
          71             :   real(r8), intent(in) :: alpha_gw_movmtn
          72             :
          73             :   ! Indices of top gravity wave source level and lowest level where wind
          74             :   ! tendencies are allowed.
          75             :   integer, intent(out) :: src_level(ncol)
          76             :   integer, intent(out) :: tend_level(ncol)
          77             :
          78             :   ! Wave Reynolds stress.
          79             :   real(r8), intent(out) :: tau(ncol,-band%ngwv:band%ngwv,pver+1) !tau = momentum flux (m2/s2) at interface level ngwv = band of phase speeds
          80             :   ! Projection of wind at midpoints and interfaces.
          81             :   real(r8), intent(out) :: ubm(ncol,pver), ubi(ncol,pver+1)
          82             :   ! Unit vectors of source wind (zonal and meridional components).
          83             :   real(r8), intent(out) :: xv(ncol), yv(ncol) !determined by vector direction of wind at source
          84             :   ! Phase speeds.
          85             :   real(r8), intent(out) :: c(ncol,-band%ngwv:band%ngwv)
          86             :
          87             :   ! Heating depth [m] and maximum heating in each column.
          88             :   real(r8), intent(out) :: hdepth(ncol)    !calculated here in this code
          89             :
          90             : !---------------------------Local Storage-------------------------------
          91             :   ! Column and (vertical) level indices.
          92             :   integer :: i, k
          93             :
          94             :   ! Zonal/meridional wind at steering level, i.e., 'cell speed'.
          95             :   ! May be later modified by retrograde motion ....
          96     2978352 :   real(r8) :: usteer(ncol), vsteer(ncol)
          97     2978352 :   real(r8) :: uwavef(ncol,pver),vwavef(ncol,pver)
          98             :   ! Steering level (integer converted to real*8)
          99     2978352 :   real(r8) :: steer_level(ncol)
         100             :   ! Retrograde motion of Cell
         101     2978352 :   real(r8) :: Cell_Retro_Speed(ncol)
         102             :
         103             :   ! Maximum heating rate.
         104     2978352 :   real(r8) :: q0(ncol), qj(ncol)
         105             :   ! unit vector components at steering level and mag
         106     2978352 :   real(r8) :: xv_steer(ncol), yv_steer(ncol), umag_steer(ncol)
         107             :   ! Bottom/top heating range index.
         108     2978352 :   integer  :: boti(ncol), topi(ncol)
         109             :   ! Index for looking up heating depth dimension in the table.
         110     2978352 :   integer  :: hd_idx(ncol)
         111             :   ! Mean wind in heating region.
         112     2978352 :   real(r8) :: uh(ncol)
         113             :   ! Min/max wavenumber for critical level filtering.
         114             :   integer :: Umini(ncol), Umaxi(ncol)
         115             :   ! Source level tau for a column.
         116     2978352 :   real(r8) :: tau0(-band%ngwv:band%ngwv)
         117             :   ! Speed of convective cells relative to storm.
         118     2978352 :   real(r8) :: CS(ncol),CS1(ncol)
         119             :   ! Wind speeds in wave direction
         120     2978352 :   real(r8) :: udiff(ncol),vdiff(ncol)
         121             :   ! "on-crest" source level wind
         122     2978352 :   real(r8) :: ubmsrc(ncol),ubisrc(ncol)
         123             :
         124             :   ! Index to shift spectra relative to ground.
         125             :   integer :: shift
         126             :   ! Other wind quantities
         127     2978352 :   real(r8) :: ut(ncol),uc(ncol),umm(ncol)
         128             :   ! Tau from moving mountain lookup table
         129     2978352 :   real(r8) :: taumm(ncol)
         130             :   ! Heating rate conversion factor.  -> tuning factors
         131             :   real(r8), parameter :: CF = 20._r8  !(1/ (5%))  -> 5% of grid cell is covered with convection
         132             :   ! Averaging length.
         133             :   real(r8), parameter :: AL = 1.0e5_r8
         134             :   ! Index for moving mountain lookuptable
         135     2978352 :   integer :: hdmm_idx(ncol), uhmm_idx(ncol)
         136             :   ! Index for ground based phase speed bin
         137     2978352 :   real(r8) :: c0(ncol,-band%ngwv:band%ngwv)
         138     2978352 :   integer :: c_idx(ncol,-band%ngwv:band%ngwv)
         139             :   ! Flux source from ShCu/PBL
         140     1489176 :   real(r8) :: xpwp_src(ncol)
         141             :   ! Manual steering level set
         142             :   integer :: Steer_k
         143             :
         144             :   !----------------------------------------------------------------------
         145             :   ! Initialize tau array
         146             :   !----------------------------------------------------------------------
         147  2478854664 :   tau = 0.0_r8
         148    24865776 :   hdepth = 0.0_r8
         149    24865776 :   q0 = 0.0_r8
         150     2978352 :   tau0 = 0.0_r8
         151             :
         152             :   !----------------------------------------------------------------------
         153             :   ! Calculate flux source from ShCu/PBL
         154             :   !----------------------------------------------------------------------
         155     1489176 :   xpwp_src = shcu_flux_src( xpwp_shcu, ncol, pver+1, alpha_gw_movmtn )
         156             :
         157             :   !------------------------------------------------------------------------
         158             :   ! Determine wind and unit vectors approximately at the source (steering level), then
         159             :   ! project winds.
         160             :   !------------------------------------------------------------------------
         161             :
         162             :   ! Winds at 'steering level'
         163     1489176 :   Steer_k = pver-1
         164    24865776 :   usteer = u(:,Steer_k)  !k defined in line21 (at specified altitude)
         165    24865776 :   vsteer = v(:,Steer_k)
         166    24865776 :   steer_level = real(Steer_k,r8)
         167             :
         168             :   ! all GW calculations on a plane, which in our case is the wind at source level -> ubi is wind in this plane
         169             :   ! Get the unit vector components and magnitude at the source level.
         170     1489176 :   call get_unit_vector(usteer, vsteer, xv_steer, yv_steer, umag_steer)
         171             :
         172             :   !-------------------------------------------------------------------------
         173             :   ! If we want to account for some retorgrade cell motion,
         174             :   ! it should be done by vector subtraction from (usteer,vsteer).
         175             :   ! We assume the retrograde motion is in the same direction as
         176             :   ! (usteer,vsteer) or the unit vector (xv_steer,yv_steer). Then, the
         177             :   ! vector retrograde motion is just:
         178             :   !      = -Cell_Retrograde_Speed * (xv_steer,yv_steer)
         179             :   ! and we would modify usteer and vsteer
         180             :   !     usteer = usteer - Cell_Retrograde_Speed * xv_steer
         181             :   !     vsteer = vsteer - Cell_Retrograde_Speed * yv_steer
         182             :   !-----------------------------------------------------------------------
         183             :   ! Cell_Retro_Speed is always =0 for now
         184             :   !-----------------------------------------------------------------------
         185    24865776 :   do i=1,ncol
         186    24865776 :      Cell_Retro_Speed(i) = min( sqrt(usteer(i)**2 + vsteer(i)**2), 0._r8)
         187             :   end do
         188    24865776 :   do i=1,ncol
         189    23376600 :      usteer(i) = usteer(i) - xv_steer(i)*Cell_Retro_Speed(i)
         190    24865776 :      vsteer(i) = vsteer(i) - yv_steer(i)*Cell_Retro_Speed(i)
         191             :   end do
         192             :   !-------------------------------------------------------------------------
         193             :   ! At this point (usteer,vsteer) is the cell-speed, or equivalently, the 2D
         194             :   ! ground based wave phase speed for moving mountain GW
         195             :   !-------------------------------------------------------------------------
         196             :
         197             :
         198             :   ! Calculate heating depth.
         199             :   !
         200             :   ! Heating depth is defined as the first height range from the bottom in
         201             :   ! which heating rate is continuously positive.
         202             :   !-----------------------------------------------------------------------
         203             :
         204             :   ! First find the indices for the top and bottom of the heating range.
         205             :   !nedt is heating profile from Zhang McFarlane (it's pressure coordinates, therefore k=0 is the top)
         206             :
         207    24865776 :   boti = 0 !bottom
         208    24865776 :   topi = 0  !top
         209             :
         210     1489176 :   if (use_gw_movmtn_pbl) then
         211    24865776 :      boti=pver
         212    24865776 :      topi=Steer_k-10 ! desc%k-5
         213             :   else
         214           0 :     do k = pver, 1, -1 !start at surface
         215           0 :        do i = 1, ncol
         216           0 :           if (boti(i) == 0) then
         217             :              ! Detect if we are outside the maximum range (where z = 20 km).
         218           0 :              if (zm(i,k) >= 20000._r8) then
         219           0 :                 boti(i) = k
         220           0 :                 topi(i) = k
         221             :              else
         222             :                 ! First spot where heating rate is positive.
         223           0 :                 if (netdt(i,k) > 0.0_r8) boti(i) = k
         224             :              end if
         225           0 :           else if (topi(i) == 0) then
         226             :              ! Detect if we are outside the maximum range (z = 20 km).
         227           0 :              if (zm(i,k) >= 20000._r8) then
         228           0 :                 topi(i) = k
         229             :              else
         230             :                 ! First spot where heating rate is no longer positive.
         231           0 :                 if (.not. (netdt(i,k) > 0.0_r8)) topi(i) = k
         232             :              end if
         233             :           end if
         234             :        end do
         235             :        ! When all done, exit.
         236           0 :        if (all(topi /= 0)) exit
         237             :     end do
         238             :   end if
         239             :   ! Heating depth in m.  (top-bottom altitudes)
         240    48242376 :   hdepth = [ ( (zm(i,topi(i))-zm(i,boti(i))), i = 1, ncol ) ]
         241     1489176 :   hd_idx = index_of_nearest(hdepth, desc%hd)
         242             :
         243             :   ! hd_idx=0 signals that a heating depth is too shallow, i.e. that it is
         244             :   ! either not big enough for the lowest table entry, or it is below the
         245             :   ! minimum allowed for this convection type.
         246             :   ! Values above the max in the table still get the highest value, though.
         247             :
         248    24865776 :   where (hdepth < max(desc%min_hdepth, desc%hd(1))) hd_idx = 0
         249             :
         250             :   ! Maximum heating rate.
         251    66112488 :   do k = minval(topi), maxval(boti)
         252   299878488 :      where (k >= topi .and. k <= boti)
         253    17870112 :         q0 = max(q0, netdt(:,k))
         254             :      end where
         255             :   end do
         256             :
         257             :   ! Multiply by conversion factor
         258             :   ! (now 20* larger than what Zhang McFarlane said as they try to describe heating over 100km grid cell)
         259    24865776 :   q0 = q0 * CF
         260    24865776 :   qj = gravit/rair*q0 ! unit conversion to m/s3
         261             :
         262             :   !-------------------------------------------------
         263             :   ! CS1 and CS should be equal in current implemen-
         264             :   ! tation.
         265             :   !-------------------------------------------------
         266    24865776 :   CS1 = sqrt( usteer**2._r8 + vsteer**2._r8 )
         267    24865776 :   CS = CS1*xv_steer + CS1*yv_steer
         268             :
         269             :   ! -----------------------------------------------------------
         270             :   ! Calculate winds in reference frame of wave (uwavef,vwavef).
         271             :   ! This is like "(U-c)" in GW literature, where U and c are in
         272             :   ! ground-based speeds in a plane perpendicular to wave fronts.
         273             :   !------------------------------------------------------------
         274    24865776 :   do i=1,ncol
         275    23376600 :      udiff(i) = u(i,topi(i)) - usteer(i)
         276    23376600 :      vdiff(i) = v(i,topi(i)) - vsteer(i)
         277  2198889576 :      do k=1,pver
         278  2174023800 :         uwavef(i, k ) = u(i, k ) - usteer(i)
         279  2197400400 :         vwavef(i, k ) = v(i, k ) - vsteer(i)
         280             :      end do
         281             :   end do
         282             :   !----------------------------------------------------------
         283             :   ! Wave relative wind at source level. This determines
         284             :   ! orientation of wave in the XY plane, and therefore the
         285             :   ! direction in which force from dissipating GW will be
         286             :   ! applied.
         287             :   !----------------------------------------------------------
         288    24865776 :   do i=1,ncol
         289    23376600 :      udiff(i) = uwavef( i, topi(i) )
         290    24865776 :      vdiff(i) = vwavef( i, topi(i) )
         291             :   end do
         292             :   !-----------------------------------------------------------
         293             :   ! Unit vector components (xv,yv) in direction of wavevector
         294             :   ! i.e., in which force will be applied
         295             :   !-----------------------------------------------------------
         296     1489176 :   call get_unit_vector(udiff , vdiff , xv, yv, ubisrc )
         297             :
         298     1489176 :   call outfld('UCELL_MOVMTN', usteer, ncol, lchnk)
         299     1489176 :   call outfld('VCELL_MOVMTN', vsteer, ncol, lchnk)
         300     1489176 :   call outfld('CS_MOVMTN', CS, ncol, lchnk)
         301     1489176 :   call outfld('STEER_LEVEL_MOVMTN',steer_level, ncol, lchnk )
         302     1489176 :   call outfld('XPWP_SRC_MOVMTN', xpwp_src , ncol, lchnk )
         303             :
         304             :   !----------------------------------------------------------
         305             :   ! Project the local wave relative wind at midpoints onto the
         306             :   !  direction of the wavevector.
         307             :   !----------------------------------------------------------
         308   139982544 :   do k = 1, pver
         309   139982544 :      ubm(:,k) = dot_2d(uwavef(:,k), vwavef(:,k), xv, yv)
         310             :   end do
         311             :   ! Source level on-crest wind
         312    24865776 :   do i=1,ncol
         313    24865776 :      ubmsrc(i) = ubm(i,topi(i))
         314             :   end do
         315             :
         316             :   !---------------------------------------------------------------
         317             :   ! adjust everything so that source level wave relative on-crest
         318             :   ! wind is always positive. Also adjust unit vector comps xv,yv
         319             :   !--------------------------------------------------------------
         320   139982544 :   do k=1,pver
         321  2314006344 :      do i=1,ncol
         322  2312517168 :         ubm(i,k) = sign( 1._r8 , ubmsrc(i) )* ubm(i,k)
         323             :      end do
         324             :   end do
         325             :   !
         326    24865776 :   do i=1,ncol
         327    23376600 :      xv(i) = sign( 1._r8 , ubmsrc(i) ) * xv(i)
         328    24865776 :      yv(i) = sign( 1._r8 , ubmsrc(i) ) * yv(i)
         329             :   end do
         330             :
         331             :
         332             :
         333             :   ! Compute the interface wind projection by averaging the midpoint winds. (both same wind profile,
         334             :   ! just at different points of the grid)
         335             :
         336             :   ! Use the top level wind at the top interface.
         337    24865776 :   ubi(:,1) = ubm(:,1)
         338             :
         339     1489176 :   ubi(:,2:pver) = midpoint_interp(ubm)
         340             :
         341             :   !-----------------------------------------------------------------------
         342             :   ! determine wind for lookup table
         343             :   ! need wind speed at the top of the convecitve cell and at the steering level
         344    24865776 :   uh = 0._r8
         345    24865776 :   do i=1,ncol
         346    23376600 :      ut(i) = ubm(i,topi(i))
         347    24865776 :      uh(i) = ut(i) - CS(i) ! wind at top in the frame moving with the cell
         348             :   end do
         349             :
         350             :   ! Set phase speeds; just use reference speeds.
         351    24865776 :   c(:,0) = 0._r8
         352             :
         353             :   !-----------------------------------------------------------------------
         354             :   ! Gravity wave sources
         355             :   !-----------------------------------------------------------------------
         356             :   ! Start loop over all columns.
         357             :   !-----------------------------------------------------------------------
         358    24865776 :   do i=1,ncol
         359             :
         360             :      !---------------------------------------------------------------------
         361             :      ! Look up spectrum only if the heating depth is large enough, else leave
         362             :      ! tau = 0.
         363             :      !---------------------------------------------------------------------
         364    24865776 :      if (.not. use_gw_movmtn_pbl) then
         365           0 :         if (hd_idx(i) > 0) then
         366             :            !------------------------------------------------------------------
         367             :            ! Look up the spectrum using depth and uh.
         368             :            !------------------------------------------------------------------
         369             :            !hdmm_idx = index_of_nearest(hdepth, desc%hd)
         370           0 :            uhmm_idx = index_of_nearest(uh, desc%uh)
         371           0 :            taumm(i) = abs(desc%mfcc(uhmm_idx(i),hd_idx(i),0))
         372           0 :            taumm(i) = taumm(i)*qj(i)*qj(i)/AL/1000._r8
         373             :            ! assign sign to MF based on the ground based phase speed, ground based phase speed = CS
         374           0 :            taumm(i) = -1._r8*sign(taumm(i),CS(i))
         375             :            !find the right phase speed bin
         376           0 :            c0(i,:) = CS(i)
         377           0 :            c_idx(i,:) = index_of_nearest(c0(i,:),c(i,:))
         378             :
         379             :            !input tau to top +1 level, interface level just below top of heating, remember it's in pressure
         380             :            ! everything is upside down (source level of GWs, level where GWs are launched)
         381           0 :            tau(i,c_idx(i,:),topi(i):topi(i)+1) = taumm(i)
         382             :
         383             :         end if ! heating depth above min and not at the pole
         384             :      else
         385   327272400 :         tau(i,0,topi(i):pver+1 ) = xpwp_src(i) ! 0.1_r8/10000._r8
         386             :      endif
         387             :
         388             :   enddo
         389             :   !-----------------------------------------------------------------------
         390             :   ! End loop over all columns.
         391             :   !-----------------------------------------------------------------------
         392             :
         393             :   ! Output the source level.
         394    24865776 :   src_level = topi
         395    24865776 :   tend_level = topi
         396             :
         397             :
         398     1489176 : end subroutine gw_movmtn_src
         399             :
         400             : ! Short routine to get the indices of a set of values rounded to their
         401             : ! nearest points on a grid.
         402     1489176 : pure function index_of_nearest(x, grid) result(idx)
         403             :   real(r8), intent(in) :: x(:)
         404             :   real(r8), intent(in) :: grid(:)
         405             :
         406             :   integer :: idx(size(x))
         407             :
         408     2978352 :   real(r8) :: interfaces(size(grid)-1)
         409             :   integer :: i, n
         410             :
         411     1489176 :   n = size(grid)
         412    23826816 :   interfaces = (grid(:n-1) + grid(2:))/2._r8
         413             :
         414    24865776 :   idx = 1
         415    22337640 :   do i = 1, n-1
         416   349610040 :      where (x > interfaces(i)) idx = i + 1
         417             :   end do
         418             :
         419     2978352 : end function index_of_nearest
         420             :
         421             : !!!!!!!!!!!!!!!!!!!!!!!!!!!
         422     1489176 : pure function shcu_flux_src (xpwp_shcu , ncol, pverx, alpha_gw_movmtn ) result(xpwp_src)
         423             :   integer, intent(in) :: ncol,pverx
         424             :   real(r8), intent(in) :: xpwp_shcu (ncol,pverx)
         425             :   real(r8), intent(in) :: alpha_gw_movmtn
         426             :
         427             :   real(r8) :: xpwp_src(ncol)
         428             :
         429             :   integer :: k, nlayers
         430             :
         431             :   !-----------------------------------
         432             :   ! Simple average over layers.
         433             :   ! Probably can do better
         434             :   !-----------------------------------
         435             :   nlayers=5
         436    24865776 :   xpwp_src(:) =0._r8
         437     8935056 :   do k = 0, nlayers-1
         438   125818056 :      xpwp_src(:) = xpwp_src(:) + xpwp_shcu(:,pverx-k)
         439             :   end do
         440    24865776 :   xpwp_src(:) = alpha_gw_movmtn * xpwp_src(:)/(1.0_r8*nlayers)
         441             :
         442     1489176 : end function shcu_flux_src
         443             :
         444           0 : end module gw_movmtn

Generated by: LCOV version 1.14