From 16c9287e629a9c918630b54b93880130fdb1997a Mon Sep 17 00:00:00 2001 From: "Zhiquan (Jake) Liu" Date: Mon, 17 Jan 2022 09:32:59 -0700 Subject: [PATCH] bug fixes for radar_rf_opt=2 (#1642) Bug fix for direct assimilation of radar reflectivity (radar_rf_opt=2) TYPE: bug fix KEYWORDS: reflectivity operator, quality control SOURCE: Zhiquan (Jake) Liu (NCAR/MMM) DESCRIPTION OF CHANGES: Problem: The namelist parameter "qthres", which is used for hydrometeor base state of TL/AD for no-rain background, should not be used in reflectivity forward operator. QC for radar_rf_opt=2 should not be confused with that of retrieval radar DA. Solution: 1. Remove the use of "qthres" in forward operator. 2. Simplify QC to keep obs with both observed and background reflectivity >=rfmin (a namelist parameter). LIST OF MODIFIED FILES: M var/da/da_radar/da_get_innov_vector_radar.inc M var/da/da_radar/da_radzicevar.inc M var/da/da_radar/da_radzicevar_upper_f.inc M var/da/da_radar/da_transform_xtoy_radar.inc M var/da/da_radar/da_transform_xtoy_radar_adj.inc M var/da/da_setup_structures/da_setup_obs_structures_radar.inc M var/da/da_setup_structures/da_setup_structures.f90 TESTS CONDUCTED: 1. WFRDA Regtest passed. 2. Qualified yes. The parallel netcdf4 mods introduced an error in DA to the develop branch that was not caught. RELEASE NOTE: Bug fixes for direct assimilation of radar reflectivity. The namelist parameter "qthres", which is used for hydrometeor base state of TL/AD for no-rain background, should not be used in reflectivity forward operator. QC for radar_rf_opt=2 should not be confused with that of retrieval radar DA. --- var/da/da_radar/da_get_innov_vector_radar.inc | 29 ++++++------------- var/da/da_radar/da_radzicevar.inc | 14 +-------- var/da/da_radar/da_radzicevar_upper_f.inc | 2 +- var/da/da_radar/da_transform_xtoy_radar.inc | 2 +- .../da_radar/da_transform_xtoy_radar_adj.inc | 2 +- .../da_setup_obs_structures_radar.inc | 2 +- .../da_setup_structures.f90 | 2 +- 7 files changed, 15 insertions(+), 38 deletions(-) diff --git a/var/da/da_radar/da_get_innov_vector_radar.inc b/var/da/da_radar/da_get_innov_vector_radar.inc index 427f186618..44fcc0c83f 100644 --- a/var/da/da_radar/da_get_innov_vector_radar.inc +++ b/var/da/da_radar/da_get_innov_vector_radar.inc @@ -96,7 +96,7 @@ subroutine da_get_innov_vector_radar (it, grid, ob, iv) qng=0 rn0_r=8e6 rn0_s=3e6 - rn0_g=4e5 + rn0_g=4e6 rhos=100.0 rhog=400.0 !------------------------ @@ -211,8 +211,8 @@ subroutine da_get_innov_vector_radar (it, grid, ob, iv) call da_interp_lin_3d (grid%xb % wh, iv%info(radar), model_w) call da_interp_lin_3d (grid%xb % rho, iv%info(radar), model_rho) call da_interp_lin_3d (grid%xb % qrn, iv%info(radar), model_qrn) - call da_interp_lin_3d (grid%xb % qrn, iv%info(radar), model_qcl) - call da_interp_lin_3d (grid%xb % qrn, iv%info(radar), model_qci) + call da_interp_lin_3d (grid%xb % qcw, iv%info(radar), model_qcl) + call da_interp_lin_3d (grid%xb % qci, iv%info(radar), model_qci) call da_interp_lin_3d (grid%xb % qsn, iv%info(radar), model_qsn) IF ( ASSOCIATED( grid%xb%qgr ) ) THEN call da_interp_lin_3d (grid%xb % qgr, iv%info(radar), model_qgr) @@ -242,7 +242,7 @@ END IF do n=iv%info(radar)%n1,iv%info(radar)%n2 - if ( use_radar_rf ) then + if ( use_radar_rf .and. radar_rf_opt==1) then ! for Xiao's default scheme ! Test 5.0E-8 (kg/kg) as critical value. It can not be smaller. do k=1,iv%info(radar)%levels(n) @@ -367,9 +367,10 @@ END IF end if if (use_radar_rf .and. radar_rf_opt==2) then + iv % radar(n) % zmm(k) % inv = 0 + iv % radar(n) % rf (k) % qc = -5 ! assume bad rf obs first + echo_rf_good = ob % radar(n) % rf(k) >= rfmin if ( echo_rf_good ) then - !select case (radar_rf_opt) - !case(2) tmk=model_tc(k,n)+273.15 prs=model_p(k,n) zmm=0 @@ -380,20 +381,8 @@ END IF 0,zmm_ref) model_rf(k,n)=dbz*radar_rf_rscl iv % radar(n) % zmm(k) % inv = zmm - iv % radar(n) % rf(k) % inv = 0 - if(ob % radar(n) % rf(k) < rfmin .and. ( model_qrn(k,n)<=2.0*rf_qthres.and. & - model_qsn(k,n)<=2.0*rf_qthres.and. & - model_qgr(k,n)<=2.0*rf_qthres)) then - iv % radar(n) % zmm(k) % inv = 0 - iv % radar(n) % rf (k) % qc = -5 - else if( ob % radar(n) % rf(k) < rfmin .and. ( model_qrn(k,n)>2.0*rf_qthres.or. & - model_qsn(k,n)>2.0*rf_qthres.or. & - model_qgr(k,n)>2.0*rf_qthres) ) then - iv % radar(n) % rf(k) % inv = rfmin - model_rf(k,n) - else if( ob % radar(n) % rf(k) >= rfmin ) then - iv % radar(n) % rf(k) % inv = ob % radar(n) % rf(k) - model_rf(k,n) - endif - !end select + iv % radar(n) % rf(k) % inv = ob % radar(n) % rf(k) - model_rf(k,n) + if (model_rf(k,n) >= rfmin) iv % radar(n) % rf (k) % qc = 0 end if end if diff --git a/var/da/da_radar/da_radzicevar.inc b/var/da/da_radar/da_radzicevar.inc index 9c974f1a25..79263236a0 100644 --- a/var/da/da_radar/da_radzicevar.inc +++ b/var/da/da_radar/da_radzicevar.inc @@ -133,7 +133,6 @@ subroutine da_radzicevar(qvp0,qra0,qsn0,qgr0,qnr0,qns0,qng0,tmk0,prs,dbz, integer :: temi01,temi02,temi03 ! temporary variable for integer real :: para_pr,para_pdx_dq,para_pdx_df ! for - real :: tc !--------------------------------------------------------------------------------- @@ -162,17 +161,6 @@ subroutine da_radzicevar(qvp0,qra0,qsn0,qgr0,qnr0,qns0,qng0,tmk0,prs,dbz, ! if(qsn0=1.and.gropt>=1) then dqra=0 dqgr=0 @@ -408,7 +396,7 @@ subroutine da_radzicevar(qvp0,qra0,qsn0,qgr0,qnr0,qns0,qng0,tmk0,prs,dbz, ! Convert to dBZ !------------------------------------- if(tlopt==0) then - dbz = max(-15.0,10. * log10(z_e)) + dbz = 10. * log10(z_e) endif ! save z_e mm^6 mm^-3 zmm=z_e diff --git a/var/da/da_radar/da_radzicevar_upper_f.inc b/var/da/da_radar/da_radzicevar_upper_f.inc index 40d9403dd6..4e3636db6c 100644 --- a/var/da/da_radar/da_radzicevar_upper_f.inc +++ b/var/da/da_radar/da_radzicevar_upper_f.inc @@ -12,7 +12,7 @@ if(flg==2) fmax=0.3 ! for hail/graupel upper_f=0 - if(qr>=qthres.and.qice>=qthres) then + if(qr>1.0e-12 .and. qice>1.0e-12) then upper_f=fmax*min(qice/qr,qr/qice)**0.3 endif diff --git a/var/da/da_radar/da_transform_xtoy_radar.inc b/var/da/da_radar/da_transform_xtoy_radar.inc index 24115ef06f..b99c0b826a 100644 --- a/var/da/da_radar/da_transform_xtoy_radar.inc +++ b/var/da/da_radar/da_transform_xtoy_radar.inc @@ -76,7 +76,7 @@ subroutine da_transform_xtoy_radar (grid, iv, y) qng=0 rn0_r=8e6 rn0_s=3e6 - rn0_g=4e5 + rn0_g=4e6 rhos=100.0 rhog=400.0 !------------------------ diff --git a/var/da/da_radar/da_transform_xtoy_radar_adj.inc b/var/da/da_radar/da_transform_xtoy_radar_adj.inc index b65cb2a985..a54fd218de 100644 --- a/var/da/da_radar/da_transform_xtoy_radar_adj.inc +++ b/var/da/da_radar/da_transform_xtoy_radar_adj.inc @@ -85,7 +85,7 @@ subroutine da_transform_xtoy_radar_adj(grid, iv, jo_grad_y, jo_grad_x) qng=0 rn0_r=8e6 rn0_s=3e6 - rn0_g=4e5 + rn0_g=4e6 rhos=100.0 rhog=400.0 !------------------------ diff --git a/var/da/da_setup_structures/da_setup_obs_structures_radar.inc b/var/da/da_setup_structures/da_setup_obs_structures_radar.inc index efc981a3df..94991e6664 100644 --- a/var/da/da_setup_structures/da_setup_obs_structures_radar.inc +++ b/var/da/da_setup_structures/da_setup_obs_structures_radar.inc @@ -100,7 +100,7 @@ subroutine da_setup_obs_structures_radar( grid, ob, iv ) ! Calculate DT for RF DA - if (use_radar_rf) then + if (use_radar_rf .and. radar_rf_opt==1) then if (.not. DT_cloud_model) then do j = jts,jte do i = its, ite diff --git a/var/da/da_setup_structures/da_setup_structures.f90 b/var/da/da_setup_structures/da_setup_structures.f90 index 7b19ce967f..582a14a112 100644 --- a/var/da/da_setup_structures/da_setup_structures.f90 +++ b/var/da/da_setup_structures/da_setup_structures.f90 @@ -21,7 +21,7 @@ module da_setup_structures analysis_date,coarse_ix,coarse_ds,map_projection,coarse_jy, c2,dsm,phic, & pole, cone_factor, start_x,base_pres,ptop,psi1,start_y, base_lapse,base_temp,truelat2_3dv, & truelat1_3dv,xlonc,t0,num_fft_factors,pi,print_detail_spectral, global, print_detail_obs, & - use_radar_rf, use_radar_rhv, use_radar_rqv, & + use_radar_rf, use_radar_rhv, use_radar_rqv, radar_rf_opt, & num_ob_indexes,kts, kte, time_window_max, time_window_min, & max_fgat_time, num_fgat_time, dt_cloud_model, & use_ssmiretrievalobs,use_radarobs,use_ssmitbobs,use_qscatobs, num_procs, use_rainobs, &