Update gw_convect

components/eam/bld/build-namelist

@@ -830,6 +830,9 @@ add_default($nl,'use_hetfrz_classnuc');
 add_default($nl,'hist_hetfrz_classnuc');
 add_default($nl,'gw_convect_hcf') if (get_default_value('gw_convect_hcf'));
 add_default($nl,'hdepth_scaling_factor') if (get_default_value('hdepth_scaling_factor'));
+add_default($nl,'gw_convect_hdepth_min') if (get_default_value('gw_convect_hdepth_min'));
+add_default($nl,'gw_convect_storm_speed_min') if (get_default_value('gw_convect_storm_speed_min'));
+add_default($nl,'use_gw_convect_old', 'val'=>'.true.');
 add_default($nl,'linoz_psc_T');
 if ($cfg->get('microphys') =~ /^mg2/) {
     add_default($nl,'micro_mg_dcs_tdep');

components/eam/bld/namelist_files/namelist_defaults_eam.xml

@@ -1890,6 +1890,8 @@ with se_tstep, dt_remap_factor, dt_tracer_factor set to -1
 <gw_convect_hcf        phys="default" > 10.0  </gw_convect_hcf>
 <hdepth_scaling_factor phys="default" > 0.50  </hdepth_scaling_factor>
 <hdepth_scaling_factor phys="default" use_MMF="1" >  1.0  </hdepth_scaling_factor>
+<gw_convect_hdepth_min      phys="default">2.5</gw_convect_hdepth_min>
+<gw_convect_storm_speed_min phys="default">10.0</gw_convect_storm_speed_min>
 <effgw_oro phys="default"> 0.375      </effgw_oro>
 <use_gw_energy_fix phys="default"> .true.      </use_gw_energy_fix>
 <clubb_C14 phys="default"> 2.5D0      </clubb_C14>

components/eam/bld/namelist_files/namelist_definition.xml

@@ -1110,6 +1110,23 @@ Scaling factor for the heating depth
 Default: 1.0
 </entry>
 
+<entry id="gw_convect_hdepth_min" type="real" category="gw_drag"
+       group="gw_drag_nl" valid_values="" >
+minimum hdepth for for convective GWD spectrum lookup table
+Default: 2.5
+</entry>
+
+<entry id="gw_convect_storm_speed_min" type="real" category="gw_drag"
+       group="gw_drag_nl" valid_values="" >
+minimum convective storm speed for convective GWD
+Default: 10.0 m/s
+</entry>
+
+<entry id="use_gw_convect_old" type="logical" category="gw_drag"
+       group="gw_drag_nl" valid_values="" >
+switch to revert to old calculation of Beres scheme for heating depth and max
+</entry>
+
 <entry id="effgw_beres" type="real" category="gw_drag"
        group="gw_drag_nl" valid_values="" >
 Efficiency associated with convective gravity waves from the Beres

components/eam/src/physics/cam/gw_convect.F90

@@ -4,10 +4,9 @@ module gw_convect
 ! This module handles gravity waves from convection, and was extracted from
 ! gw_drag in May 2013.
 !
-
-use gw_utils, only: r8
-
-use gw_common, only: pver, pgwv
+use cam_logfile, only: iulog
+use gw_utils,    only: r8
+use gw_common,   only: pver, pgwv
 
 implicit none
 private
@@ -21,8 +20,8 @@ module gw_convect
 ! Dimension for mean wind in heating.
 integer :: maxuh
 
-! Index for level at 700 mb.
-integer :: k700
+! Index for level for storm/steering flow (usually 700 mb)
+integer :: k_src_wind
 
 ! Table of source spectra.
 real(r8), allocatable :: mfcc(:,:,:)
@@ -31,19 +30,20 @@ module gw_convect
 
 !==========================================================================
 
-subroutine gw_convect_init(k700_in, mfcc_in, errstring)
-  ! Index at 700 mb.
-  integer, intent(in) :: k700_in
-  ! Source spectra to keep as table.
-  real(r8), intent(in) :: mfcc_in(:,:,:)
-  ! Report any errors from this routine.
-  character(len=*), intent(out) :: errstring
-
+subroutine gw_convect_init( plev_src_wind, mfcc_in, errstring)
+  use ref_pres, only: pref_edge
+  real(r8), intent(in) :: plev_src_wind        ! previously hardcoded to 70000._r8
+  real(r8), intent(in) :: mfcc_in(:,:,:)       ! Source spectra to keep as table
+  character(len=*), intent(out) :: errstring   ! Report any errors from this routine
   integer :: ierr
 
   errstring = ""
 
-  k700 = k700_in
+  do k = 0, pver
+    if ( pref_edge(k+1) < plev_src_wind ) k_src_wind = k+1
+  end do
+
+  if (masterproc) write (iulog,*) 'gw_convect: steering flow level = ',k_src_wind
 
   ! First dimension is maxh.
   maxh = size(mfcc_in,1)
@@ -60,7 +60,9 @@ end subroutine gw_convect_init
 
 subroutine gw_beres_src(ncol, ngwv, lat, u, v, netdt, &
      zm, src_level, tend_level, tau, ubm, ubi, xv, yv, c, &
-     hdepth, maxq0, CF, hdepth_scaling_factor)
+     hdepth, maxq0_out, maxq0_conversion_factor, hdepth_scaling_factor, &
+     hdepth_min, storm_speed_min, &
+     use_gw_convect_old)
 !-----------------------------------------------------------------------
 ! Driver for multiple gravity wave drag parameterization.
 !
@@ -90,11 +92,20 @@ subroutine gw_beres_src(ncol, ngwv, lat, u, v, netdt, &
   real(r8), intent(in) :: zm(ncol,pver)
 
   ! Heating conversion factor
-  real(r8), intent(in) :: CF
+  real(r8), intent(in) :: maxq0_conversion_factor
 
   ! Scaling factor for the heating depth
   real(r8), intent(in) :: hdepth_scaling_factor
 
+  ! minimum hdepth for for spectrum lookup table
+  real(r8), intent(in) :: hdepth_min
+
+  ! minimum convective storm speed
+  real(r8), intent(in) :: storm_speed_min
+
+  ! switch for restoring legacy method
+  logical, intent(in) :: use_gw_convect_old
+
   ! Indices of top gravity wave source level and lowest level where wind
   ! tendencies are allowed.
   integer, intent(out) :: src_level(ncol)
@@ -110,19 +121,19 @@ subroutine gw_beres_src(ncol, ngwv, lat, u, v, netdt, &
   real(r8), intent(out) :: c(ncol,-pgwv:pgwv)
 
   ! Heating depth and maximum heating in each column.
-  real(r8), intent(out) :: hdepth(ncol), maxq0(ncol)
+  real(r8), intent(out) :: hdepth(ncol), maxq0_out(ncol)
 
 !---------------------------Local Storage-------------------------------
   ! Column and level indices.
   integer :: i, k
 
-  ! Zonal/meridional wind at 700mb.
-  real(r8) :: u700(ncol), v700(ncol)
+  ! Zonal/meridional source wind
+  real(r8) :: u_src(ncol), v_src(ncol)
   ! 3.14...
   real(r8), parameter :: pi = 4._r8*atan(1._r8)
 
   ! Maximum heating rate.
-  real(r8) :: q0(ncol)
+  real(r8) :: maxq0(ncol)
 
   ! Bottom/top heating range index.
   integer  :: mini(ncol), maxi(ncol)
@@ -135,36 +146,36 @@ subroutine gw_beres_src(ncol, ngwv, lat, u, v, netdt, &
   ! Source level tau for a column.
   real(r8) :: tau0(-PGWV:PGWV)
   ! Speed of convective cells relative to storm.
-  integer :: CS(ncol)
+  integer :: storm_speed(ncol)
   ! Index to shift spectra relative to ground.
   integer :: shift
 
-  ! Heating rate conversion factor. Change to take the value from caller and controllable by namelist (to tune QBO)
-  ! real(r8), parameter :: CF = 20._r8
-  ! Averaging length.
-  real(r8), parameter :: AL = 1.0e5_r8
+  ! fixed parameters (we may want to expose these in the namelist for tuning)
+  real(r8), parameter :: tau_avg_length       = 1.0e5_r8 ! spectrum averaging length
+  real(r8), parameter :: heating_altitude_max = 20e3     ! max altitude to check heating (probably don't need this)
+
+  integer :: ndepth_pos
+  integer :: ndepth_tot
 
   !----------------------------------------------------------------------
   ! Initialize tau array
   !----------------------------------------------------------------------
 
-  tau = 0.0_r8
+  tau    = 0.0_r8
   hdepth = 0.0_r8
-  q0 = 0.0_r8
-  tau0 = 0.0_r8
+  maxq0  = 0.0_r8
+  tau0   = 0.0_r8
 
   !------------------------------------------------------------------------
-  ! Determine 700 mb layer wind and unit vectors, then project winds.
+  ! Determine source layer wind and unit vectors, then project winds.
   !------------------------------------------------------------------------
 
-  ! Just use the 700 mb interface values for the source wind speed and
-  ! direction (unit vector).
-
-  u700 = u(:,k700)
-  v700 = v(:,k700)
+  ! source wind speed and direction
+  u_src = u(:,k_src_wind)
+  v_src = v(:,k_src_wind)
 
   ! Get the unit vector components and magnitude at the surface.
-  call get_unit_vector(u700, v700, xv, yv, ubi(:,k700))
+  call get_unit_vector(u_src, v_src, xv, yv, ubi(:,k_src_wind))
 
   ! Project the local wind at midpoints onto the source wind.
   do k = 1, pver
@@ -184,33 +195,62 @@ subroutine gw_beres_src(ncol, ngwv, lat, u, v, netdt, &
   ! which heating rate is continuously positive.
   !-----------------------------------------------------------------------
 
-  ! First find the indices for the top and bottom of the heating range.
+  ! Find indices for the top and bottom of the heating range.
   mini = 0
   maxi = 0
-  do k = pver, 1, -1
-     do i = 1, ncol
+
+  if (use_gw_convect_old) then
+    !---------------------------------------------------------------------
+    ! original version used in CAM4/5/6 and EAMv1/2/3
+    do k = pver, 1, -1
+      do i = 1, ncol
         if (mini(i) == 0) then
-           ! Detect if we are outside the maximum range (where z = 20 km).
-           if (zm(i,k) >= 20000._r8) then
-              mini(i) = k
-              maxi(i) = k
-           else
-              ! First spot where heating rate is positive.
-              if (netdt(i,k) > 0.0_r8) mini(i) = k
-           end if
+          ! Detect if we are outside the maximum range (where z = 20 km).
+          if (zm(i,k) >= heating_altitude_max) then
+            mini(i) = k
+            maxi(i) = k
+          else
+            ! First spot where heating rate is positive.
+            if (netdt(i,k) > 0.0_r8) mini(i) = k
+          end if
         else if (maxi(i) == 0) then
-           ! Detect if we are outside the maximum range (z = 20 km).
-           if (zm(i,k) >= 20000._r8) then
-              maxi(i) = k
-           else
-              ! First spot where heating rate is no longer positive.
-              if (.not. (netdt(i,k) > 0.0_r8)) maxi(i) = k
-           end if
+          ! Detect if we are outside the maximum range (z = 20 km).
+          if (zm(i,k) >= heating_altitude_max) then
+            maxi(i) = k
+          else
+            ! First spot where heating rate is no longer positive.
+            if (.not. (netdt(i,k) > 0.0_r8)) maxi(i) = k
+          end if
         end if
-     end do
-     ! When all done, exit
-     if (all(maxi /= 0)) exit
-  end do
+      end do
+      ! When all done, exit
+      if (all(maxi /= 0)) exit
+    end do
+    !---------------------------------------------------------------------
+  else
+    !---------------------------------------------------------------------
+    ! cleaner version that addresses bug in original where heating max and
+    ! depth were too low whenever heating <=0 occurred in the middle of
+    ! the heating profile (ex. at the melting level)
+    do i = 1, ncol
+      do k = pver, 1, -1
+        if ( zm(i,k) < heating_altitude_max ) then
+          if ( netdt(i,k) > 0.0_r8 ) then
+            ! Set mini as first spot where heating rate is positive
+            if ( mini(i)==0 ) mini(i) = k
+            ! Set maxi to current level
+            maxi(i) = k
+          end if
+        else
+          ! above the max check if indices were found
+          if ( mini(i)==0 ) mini(i) = k
+          if ( maxi(i)==0 ) maxi(i) = k
+        end if
+      end do
+    end do
+    !---------------------------------------------------------------------
+   end if
+
 
   ! Heating depth in km.
   hdepth = [ ( (zm(i,maxi(i))-zm(i,mini(i)))/1000._r8, i = 1, ncol ) ]
@@ -223,19 +263,19 @@ subroutine gw_beres_src(ncol, ngwv, lat, u, v, netdt, &
   ! Maximum heating rate.
   do k = minval(maxi), maxval(mini)
      where (k >= maxi .and. k <= mini)
-        q0 = max(q0, netdt(:ncol,k))
+        maxq0 = max(maxq0, netdt(:ncol,k))
      end where
   end do
 
   !output max heating rate in K/day
-  maxq0 = q0*24._r8*3600._r8
+  maxq0_out = maxq0*24._r8*3600._r8
 
   ! Multipy by conversion factor
-  q0 = q0 * CF
+  maxq0 = maxq0 * maxq0_conversion_factor
 
-  ! Taking ubm at 700 mb to be the storm speed, find the cell speed where
-  ! the storm speed is > 10 m/s.
-  CS = int(sign(max(abs(ubm(:,k700))-10._r8, 0._r8), ubm(:,k700)))
+  ! Taking ubm at assumed source level to be the storm speed, 
+  ! find the cell speed where the storm speed is > storm_speed_min
+  storm_speed = int(sign(max(abs(ubm(:,k_src_wind))-storm_speed_min, 0._r8), ubm(:,k_src_wind)))
 
   uh = 0._r8
   do k = minval(maxi), maxval(mini)
@@ -244,7 +284,7 @@ subroutine gw_beres_src(ncol, ngwv, lat, u, v, netdt, &
      end where
   end do
 
-  uh = uh - real(CS, r8)
+  uh = uh - real(storm_speed, r8)
 
   ! Limit uh to table range.
   uh = min(uh, real(maxuh, r8))
@@ -270,8 +310,19 @@ subroutine gw_beres_src(ncol, ngwv, lat, u, v, netdt, &
      !---------------------------------------------------------------------
      ! Look up spectrum only if depth >= 2.5 km, else set tau0 = 0.
      !---------------------------------------------------------------------
-
-     if ((hdepth(i) >= 2.5_r8) .and. (abs(lat(i)) < (pi/2._r8))) then
+      if ((hdepth(i) >= hdepth_min)) then
+         ndepth_pos = 0
+         ndepth_tot = 0
+         do k = 1,pver
+            if ( k>=maxi(i).and.k<=mini(i) )  then
+               ndepth_tot = ndepth_tot + 1
+               if ( netdt(i,k)>0 )  ndepth_pos = ndepth_pos + 1
+            end if
+         end do
+         ! write (iulog,*) 'WHDEBUG - i: ',i,' Hd: ',hdepth(i),' Hn: ',ndepth_pos,' N: ',ndepth_tot
+      end if
+
+     if ((hdepth(i) >= hdepth_min) .and. (abs(lat(i)) < (pi/2._r8))) then
 
         !------------------------------------------------------------------
         ! Look up the spectrum using depth and uh.
@@ -280,11 +331,11 @@ subroutine gw_beres_src(ncol, ngwv, lat, u, v, netdt, &
         tau0 = mfcc(NINT(hdepth(i)),NINT(uh(i)),:)
 
         ! Shift spectrum so that it is relative to the ground.
-        shift = -nint(real(CS(i), r8)/dc)
+        shift = -nint(storm_speed(i)/dc)
         tau0 = cshift(tau0,shift)
 
         ! Adjust magnitude.
-        tau0 = tau0*q0(i)*q0(i)/AL
+        tau0 = tau0*maxq0(i)*maxq0(i)/tau_avg_length
 
         ! Adjust for critical level filtering.
         Umini = max(nint(Umin(i)/dc),-PGWV)