Use Qv or specific humidity in interpolation if they are available (#…

…1959)

TYPE: enhancement

KEYWORDS: Qv, specific humidity, vertical interpolation

SOURCE: internal

DESCRIPTION OF CHANGES:
Problem:
In current real program, if either water mixing ratio (Qv) or specific humidity (SH) is available, it first converts it to RH, perform vertical interpolation using RH, and then compute Qv for the model input. Assuming Qv and SH are coming from other model's native levels, which tends to be plenty and generally true, it should not be necessary to go through this Qv (SH) to RH, and RH to Qv steps.

Solution:
Use the flag_qv and flag_sh to avoid using interpolated RH (even though RH is still computed and vertically interpolated). If flag_sh is 1, SH is converted to Qv or mixing ratio. Then Qv is interpolated vertically. A namelist, use_sh_qv, is added for user to choose to do so if SH or Qv data is available. A test case shows that there is some differences in the final Qv in the lower levels (a bit less Qv), particularly near terrain slopes.

LIST OF MODIFIED FILES:
M Registry/Registry.EM_COMMON
M dyn_em/module_initialize_real.F
M run/README.namelist

TESTS CONDUCTED:
- It doesn't affect existing option if RH is the only incoming variable. Tested the code, and it gives bit-for-bit results before and after this change.
- The Jenkins tests are all passing.

RELEASE NOTE: If water mixing ratio or specific humidity is available in the metgrid output, they can be used directly in vertical interpolation. Previously these fields are first converted to RH, interpolated vertically, and the water vapor mixing ratio is computed at model levels. Now one can use the new namelist, use_sh_qv = T to choose to use SH or Qv data directly.

Registry/Registry.EM_COMMON

@@ -2295,6 +2295,7 @@ rconfig integer interp_method_type namelist,domains 1 2
 rconfig logical aggregate_lu namelist,domains 1 .false. irh "aggregate_lu" "T/F aggregate the grass, shrubs, trees in LU"
 rconfig logical rh2qv_wrt_liquid namelist,domains 1 .true. irh "rh2qv_wrt_liquid" "T = rh=>Qv assumes RH wrt liquid water, F = allows ice"
 rconfig integer rh2qv_method namelist,domains 1 1 irh "rh2qv_method" "1=old MM5 method, 2=new WMO method"
+rconfig logical use_sh_qv namelist,domains 1 .false. irh "use_sh_qv" "T/F whether to use SH or mixing ratio in input"
 rconfig real qv_max_p_safe namelist,domains 1 10000 irh "qv_max_p_safe" "Threshhold pressure, Qv > flag set to value" "Pa"
 rconfig real qv_max_flag namelist,domains 1 1.E-5 irh "qv_max_flag" "Qv flag for max" "kg kg{-1}"
 rconfig real qv_max_value namelist,domains 1 3.E-6 irh "qv_max_value" "Qv value for max" "kg kg{-1}"

dyn_em/module_initialize_real.F

@@ -1181,6 +1181,7 @@ SUBROUTINE init_domain_rk ( grid &
  END IF
 
  ! Some data sets do not provide a 3d geopotential height field. 
+ ! This calculation is more accurate if the data is bottom-up.
 
  IF ( grid%ght_gc(i_valid,grid%num_metgrid_levels/2,j_valid) .LT. 1 ) THEN
  DO j = jts, MIN(jte,jde-1)
@@ -1239,6 +1240,15 @@ SUBROUTINE init_domain_rk ( grid &
  END DO
  END IF
 
+ IF ( flag_sh .EQ. 1 ) THEN
+ DO j = jts, min(jde-1,jte)
+ DO i = its, min(ide-1,ite)
+ IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
+ grid%q2(i,j)=grid%qv_gc(i,1,j)
+ END DO
+ END DO
+ END IF
+
  ! The requested ptop for real data cases.
 
  p_top_requested = grid%p_top_requested
@@ -1708,6 +1718,23 @@ SUBROUTINE init_domain_rk ( grid &
  ims , ime , jms , jme , kms , kme , &
  its , ite , jts , jte , kts , kte )
 
+ ! when specific humidity is available, qv_gc is computed from sh_gc
+ IF (config_flags%use_sh_qv .and. (flag_sh .eq. 1 .or. flag_qv .eq. 1)) THEN
+ CALL vert_interp ( grid%qv_gc , grid%pd_gc , moist(:,:,:,P_QV) , grid%pb , &
+ grid%hgtmaxw , grid%hgttrop , grid%pmaxw , grid%ptrop , &
+ grid%pmaxwnn , grid%ptropnn , &
+ 0 , 0 , &
+ config_flags%maxw_horiz_pres_diff , config_flags%trop_horiz_pres_diff , &
+ config_flags%maxw_above_this_level , &
+ num_metgrid_levels , 'Q' , &
+ interp_type , lagrange_order , extrap_type , &
+ lowest_lev_from_sfc , use_levels_below_ground , use_surface , &
+ zap_close_levels , force_sfc_in_vinterp , grid%id , &
+ ids , ide , jds , jde , kds , kde , &
+ ims , ime , jms , jme , kms , kme , &
+ its , ite , jts , jte , kts , kte )
+ END IF
+
  ! If this is theta being interpolated, AND we have extra levels for temperature,
  ! convert those extra levels (trop and max wind) to potential temp.
 
@@ -1778,6 +1805,8 @@ SUBROUTINE init_domain_rk ( grid &
  its , ite , jts , jte , kts , kte )
  END IF
 
+ ! do not compute qv from RH if flag_sh or flag_qv = 1, or use_sh_qv = F
+ IF ( flag_sh .ne. 1 .or. flag_qv .ne. 1 .or. .not. config_flags%use_sh_qv ) THEN
  IF ( config_flags%rh2qv_method .eq. 1 ) THEN
  CALL rh_to_mxrat1(grid%u_1, grid%v_1, grid%p , moist(:,:,:,P_QV) , &
  config_flags%rh2qv_wrt_liquid , &
@@ -1799,6 +1828,7 @@ SUBROUTINE init_domain_rk ( grid &
  ims , ime , jms , jme , kms , kme , &
  its , ite , jts , jte , kts , kte-1 )
  END IF
+ END IF
 
  IF ( .NOT. config_flags%interp_theta ) THEN
  CALL t_to_theta ( grid%t_2 , grid%p , p00 , &
@@ -4042,6 +4072,7 @@ SUBROUTINE init_domain_rk ( grid &
  ims , ime , jms , jme , kms , kme , &
  its , ite , jts , jte , kts , kte )
 
+ IF ( flag_sh .ne. 1 .or. flag_qv .ne. 1 .or. .not. config_flags%use_sh_qv ) THEN
  IF ( config_flags%rh2qv_method .eq. 1 ) THEN
  CALL rh_to_mxrat1(grid%u_1, grid%v_1, grid%p_hyd , moist(:,:,:,P_QV) , &
  config_flags%rh2qv_wrt_liquid , &
@@ -4063,6 +4094,7 @@ SUBROUTINE init_domain_rk ( grid &
  ims , ime , jms , jme , kms , kme , &
  its , ite , jts , jte , kts , kte-1 )
  END IF
+ END IF
  
  ! Compute pressure similarly to how computed within model, with final Qv.
 

run/README.namelist

@@ -251,7 +251,8 @@ Namelist variables specifically for the WPS input for real:
  rh2qv_method = 1, ! which method to use to computer mixing ratio from RH:
  default is option 1, the old MM5 method; option 2 uses a WMO 
  recommended method (WMO-No. 49, corrigendum, August 2000) - 
- there is a difference between the two methods though small
+ use_sh_qv = .false., ! whether to use specific humidity or mixing ratio data from input
+ recommended if input data has high vertical resolution
  interp_theta = .false. ! If set to .false., it will vertically interpolate temperature 
  instead of potential temperature, which may reduce bias when 
  compared with input data