diff --git a/Registry/Registry.EM_COMMON b/Registry/Registry.EM_COMMON index 4c92f156b3..668332c6eb 100644 --- a/Registry/Registry.EM_COMMON +++ b/Registry/Registry.EM_COMMON @@ -1021,7 +1021,7 @@ state real FM ij misc 1 - - "FM" state real FH ij misc 1 - - "FH" "INTEGRATED FUNCTION FOR HEAT" "" state real WSPD ij misc 1 - - "WSPD" "Wind Speed At Lowest Model Level (may contain vconv)" "m s-1" i1 real GZ1OZ0 ij misc 1 - - "GZ1OZ0" "LOG OF Z1 over Z0" "" -state real BR ij misc 1 - - "BR" "Bulk Richardson" "" +state real BR ij misc 1 - h "BR" "Bulk Richardson" "" state real ZOL ij misc 1 - - "ZOL" "z/L" "" # ysupbl variables for grims shallow convection @@ -1915,7 +1915,7 @@ state integer NYEAR - misc 1 - r "NY state real NDAY - misc 1 - r "NDAY" "ACCUM TIMESTEPS IN A DAY" "" state real XLAND ij misc 1 - i02rhd=(interp_fcnm_imask)u=(copy_fcnm) "XLAND" "LAND MASK (1 FOR LAND, 2 FOR WATER)" "" state real cplmask i{ncpldom}j misc 1 z i0r "CPLMASK" "COUPLING MASK (0:VALUE FROM SST UPDATE; 1:VALUE FROM COUPLED OCEAN), vertical dim is number of external domains" "" -state real ZNT ij misc 1 - i3r "ZNT" "TIME-VARYING ROUGHNESS LENGTH" "m" +state real ZNT ij misc 1 - i3rh "ZNT" "TIME-VARYING ROUGHNESS LENGTH" "m" state real CK ij misc 1 - r "CK" "ENTHALPY EXCHANGE COEFF AT 10 m" "" state real CKA ij misc 1 - r "CKA" "ENTHALPY EXCHANGE COEFF AT LOWEST MODEL LVL" "" state real CD ij misc 1 - r "CD" "DRAG COEFF AT 10m" "" diff --git a/Registry/Registry.EM_COMMON.var b/Registry/Registry.EM_COMMON.var index 2802a80680..89a34efa75 100644 --- a/Registry/Registry.EM_COMMON.var +++ b/Registry/Registry.EM_COMMON.var @@ -148,7 +148,7 @@ state real cf3 - misc - - irh "cf3 # State for derived time quantities. state integer itimestep - - - - rh "itimestep" "" "" -state real xtime - - - - rh "xtime" "minutes since simulation start" "" +state real xtime - - - - irh "xtime" "minutes since simulation start" "" state real julian - - - - - "julian" "day of year, 0.0 at 0Z on 1 Jan." "days" # input file descriptor for lbcs on parent domain @@ -305,8 +305,8 @@ state real max_msftx - misc - - rh "m state real max_msfty - misc - - rh "max_mstfy" "Max map factor in domain" "" state logical v4_metgrid - misc - - - "v4_metgrid" "for real, T/F: identify if this is v4 metgrid data" "" -state real RAINC ij misc 1 - rhdu "RAINC" "ACCUMULATED TOTAL CUMULUS PRECIPITATION" "mm" -state real RAINNC ij misc 1 - rhdu "RAINNC" "ACCUMULATED TOTAL GRID SCALE PRECIPITATION" "mm" +state real RAINC ij misc 1 - irhdu "RAINC" "ACCUMULATED TOTAL CUMULUS PRECIPITATION" "mm" +state real RAINNC ij misc 1 - irhdu "RAINNC" "ACCUMULATED TOTAL GRID SCALE PRECIPITATION" "mm" state real RAINCV ij misc 1 - r "RAINCV" "TIME-STEP CUMULUS PRECIPITATION" "mm" state real RAINNCV ij misc 1 - r "RAINNCV" "TIME-STEP NONCONVECTIVE PRECIPITATION" "mm" @@ -346,6 +346,7 @@ state real HFX ij misc 1 - irh "H state real QFX ij misc 1 - irh "QFX" "UPWARD MOISTURE FLUX AT THE SURFACE" "kg m-2 s-1" state real REGIME ij misc 1 - irh "REGIME" "FLAGS: 1=Night/Stable, 2=Mechanical Turbulent, 3=Forced Conv, 4=Free Conv" "" state integer KPBL ij misc 1 - irh "KPBL" "LEVEL OF PBL TOP" "" +state real BR ij misc 1 - irh "BR" "Bulk Richardson" "" # #--------------------------------------------------------------------------------------------------------------------------------------- diff --git a/Registry/registry.var b/Registry/registry.var index b9a41e7c12..f42e0845b7 100644 --- a/Registry/registry.var +++ b/Registry/registry.var @@ -207,6 +207,7 @@ rconfig integer varbc_tamdar_bm namelist,wrfvar4 1 1 - "For rconfig integer varbc_tamdar_nbgerr namelist,wrfvar4 1 1000 - "Scaling B-beta" "" "" rconfig integer varbc_tamdar_nobsmin namelist,wrfvar4 1 6 - "OB count req to run VarBC""" "" rconfig real varbc_tamdar_pred0 namelist,wrfvar4 1 1.0 - "Predictor(1)" "" "" +rconfig logical qscat_neutral namelist,wrfvar4 1 .false. - "qscat_neutral" "Assimilation of scatterometer data as equivalent-neutral wind. (Hersbach H, 2010)" "" rconfig logical check_max_iv namelist,wrfvar5 1 .true. - "check_max_iv" "" "" rconfig real max_error_t namelist,wrfvar5 1 5.0 - "max_error_t" "" "" rconfig real max_error_uv namelist,wrfvar5 1 5.0 - "max_error_uv" "" "" @@ -268,6 +269,8 @@ rconfig integer report_end namelist,wrfvar5 1 10000000 - "rep rconfig integer tovs_start namelist,wrfvar5 1 1 - "tovs_start" "" "" rconfig integer tovs_end namelist,wrfvar5 1 10000000 - "tovs_end" "" "" rconfig logical gpsref_thinning namelist,wrfvar5 1 .false. - "gpsref_thinning" "" "" +rconfig logical qc_qscat namelist,wrfvar5 1 .false. - "qc_qscat" "real time rain quality control for scat" "" +rconfig real qc_qscat_rr namelist,wrfvar5 1 2.5 - "qc_qscat_rr" "rain rate threashold for real time rain quality control" "" rconfig integer max_ext_its namelist,wrfvar6 1 1 - "max_ext_its" "" "" rconfig integer ntmax namelist,wrfvar6 max_outer_iterations 75 - "ntmax" "" "" rconfig logical use_inverse_squarerootb namelist,wrfvar6 1 .false. - "use_inverse_squarerootb" "" "" diff --git a/var/da/da_qscat/da_get_innov_vector_qscat.inc b/var/da/da_qscat/da_get_innov_vector_qscat.inc index 5f4273f6bf..b11b75526a 100644 --- a/var/da/da_qscat/da_get_innov_vector_qscat.inc +++ b/var/da/da_qscat/da_get_innov_vector_qscat.inc @@ -19,6 +19,7 @@ subroutine da_get_innov_vector_qscat (it,num_qcstat_conv, grid, ob, iv) real :: dx, dxm ! Interpolation weights. real :: dy, dym ! Interpolation weights. real :: speed, direction + real :: znt, rib, ht, za, bn, fm, bd real, allocatable :: model_u(:,:) ! Model value u at ob location. real, allocatable :: model_v(:,:) ! Model value v at ob location. @@ -79,6 +80,33 @@ subroutine da_get_innov_vector_qscat (it,num_qcstat_conv, grid, ob, iv) call da_interp_lin_3d (grid%xb%v, iv%info(qscat), model_v) #endif + if (qscat_neutral) then + do n=iv%info(qscat)%n1,iv%info(qscat)%n2 + ! [1.1] Get horizontal interpolation weights: + i = iv%info(qscat)%i(1,n) + j = iv%info(qscat)%j(1,n) + dx = iv%info(qscat)%dx(1,n) + dy = iv%info(qscat)%dy(1,n) + dxm = iv%info(qscat)%dxm(1,n) + dym = iv%info(qscat)%dym(1,n) + + znt = dym*(dxm*grid%znt(i,j) + dx*grid%znt(i+1,j)) + dy*(dxm*grid%znt(i,j+1) + dx*grid%znt(i+1,j+1)) + rib = dym*(dxm*grid%br(i,j) + dx*grid%br(i+1,j)) + dy*(dxm*grid%br(i,j+1) + dx*grid%br(i+1,j+1)) + ht = dym*(dxm*grid%ht(i,j) + dx*grid%ht(i+1,j)) + dy*(dxm*grid%ht(i,j+1) + dx*grid%ht(i+1,j+1)) + + za = v_h(kts) - ht + + bn = log( 1 + za/znt ) + if (rib < 0) then + fm=1-10 * rib / ( 1-75 * rib * bn * sqrt( 1+10/znt ) ) + else + fm = 1 / (1 + 10 * rib / sqrt(1 + 5*rib ) ) + end if + bd = bn/sqrt(fm) + model_u(1, n) = model_u(1, n) * log( 1 + 10/znt ) / bd + model_v(1, n) = model_v(1, n) * log( 1 + 10/znt ) / bd + end do + end if do n=iv%info(qscat)%n1,iv%info(qscat)%n2 !------------------------------------------------------------------------ @@ -120,6 +148,8 @@ subroutine da_get_innov_vector_qscat (it,num_qcstat_conv, grid, ob, iv) if ( check_max_iv ) & call da_check_max_iv_qscat(iv, it, num_qcstat_conv) + if ( qc_qscat ) call da_qc_qscat(it, ob, grid, iv, num_qcstat_conv) + deallocate (model_u) deallocate (model_v) diff --git a/var/da/da_qscat/da_qc_qscat.inc b/var/da/da_qscat/da_qc_qscat.inc new file mode 100644 index 0000000000..81a0623723 --- /dev/null +++ b/var/da/da_qscat/da_qc_qscat.inc @@ -0,0 +1,69 @@ +subroutine da_qc_qscat (it, ob, grid, iv, num_qcstat_conv) + + !--------------------------------------------------------------------------- + ! Purpose: perform quality control for scatterometer data. + ! + ! METHOD: separated QC for each sensor + !--------------------------------------------------------------------------- + + implicit none + + integer , intent(in) :: it ! outer loop count + type (y_type), intent(in) :: ob ! Observation structure. + type(domain), intent(in) :: grid ! first guess state. + type (iv_type), intent(inout) :: iv ! O-B structure. + integer, intent(inout) :: num_qcstat_conv(:,:,:,:) + + integer :: n ! Loop counter. + real, allocatable :: model_rainc(:) ! interpolated rainc + real, allocatable :: model_rainnc(:) ! interpolated rainnc + real :: rr ! rain rate + real*8 :: obs_time, analysis_time + integer :: iyear, imonth, iday, ihour, imin + real :: ub, vb, sb, db, uo, vo, so, do ! + + if (trace_use_dull) call da_trace_entry("da_qc_qscat") + + allocate (model_rainc(iv%info(qscat)%n1:iv%info(qscat)%n2)) + allocate (model_rainnc(iv%info(qscat)%n1:iv%info(qscat)%n2)) + + call da_interp_lin_2d (grid%rainc, iv%info(qscat), 1, model_rainc) + call da_interp_lin_2d (grid%rainnc, iv%info(qscat), 1, model_rainnc) + + ! read (analysis_date,'(i4,4(1x,i2))') iyear, imonth, iday, ihour, imin + ! call da_get_julian_time (iyear,imonth,iday,ihour,imin,analysis_time) + + do n=iv%info(qscat)%n1, iv%info(qscat)%n2 + ! if ( TRIM(iv%info(qscat)%id(n)) .ne. "CFOSAT" ) cycle ! only for CFOSAT + if ( iv%qscat(n)%u%qc < 0 ) cycle + + ! read (iv%info(qscat)%date_char(n),'(i4,4(1x,i2))') iyear, imonth, iday, ihour, imin + ! call da_get_julian_time (iyear,imonth,iday,ihour,imin,obs_time) ! get obs time + ! if ( obs_time > analysis_time) cycle ! if obs time greater then ana time, pass + if (grid%xtime <= 0) then + rr = 0 + else + rr = (model_rainc(n) + model_rainnc(n)) / (grid%xtime / 60.0) ! rain rate (mm/h) = total/time + end if + ! print *, "rr=",rr + if ( rr > qc_qscat_rr ) then + iv%qscat(n)%u%qc = -16 + iv%qscat(n)%u%inv = 0.0 + iv%qscat(n)%v%qc = -16 + iv%qscat(n)%v%inv = 0.0 + num_qcstat_conv(2,qscat,1,1) = num_qcstat_conv(2,qscat,1,1) + 1 + num_qcstat_conv(2,qscat,2,1) = num_qcstat_conv(2,qscat,2,1) + 1 + + if ( write_rej_obs_conv ) write(qcstat_conv_unit,'(2x,a10,2x,a4,2f12.2,f12.2)')& + 'qscat',ob_vars(1),iv%info(qscat)%lat(1,n),iv%info(qscat)%lon(1,n), rr + if ( write_rej_obs_conv ) write(qcstat_conv_unit,'(2x,a10,2x,a4,2f12.2,f12.2)')& + 'qscat',ob_vars(2),iv%info(qscat)%lat(1,n),iv%info(qscat)%lon(1,n), rr + endif + end do + + deallocate (model_rainc) + deallocate (model_rainnc) + + if (trace_use_dull) call da_trace_exit("da_qc_qscat") + +end subroutine da_qc_qscat \ No newline at end of file diff --git a/var/da/da_qscat/da_qscat.f90 b/var/da/da_qscat/da_qscat.f90 index 619b422d9e..ed9d7d5f77 100644 --- a/var/da/da_qscat/da_qscat.f90 +++ b/var/da/da_qscat/da_qscat.f90 @@ -10,13 +10,13 @@ module da_qscat max_error_bt, max_error_buv, anal_type_verify, kms,kme,kts,kte,& ob_vars,qcstat_conv_unit, fails_error_max, & convert_fd2uv,convert_uv2fd,max_error_spd,max_error_dir,max_omb_spd,max_omb_dir,pi,qc_rej_both,& - wind_sd_qscat, wind_stats_sd, write_rej_obs_conv + wind_sd_qscat, wind_stats_sd, write_rej_obs_conv, qc_qscat, qc_qscat_rr, qscat_neutral use da_grid_definitions, only : da_ffdduv, da_ffdduv_model, da_ffdduv_diagnose use da_physics, only : da_uv_to_sd_lin, da_uv_to_sd_adj use da_define_structures, only : maxmin_type, iv_type, y_type, jo_type, & bad_data_type, x_type, number_type, bad_data_type use da_interpolation, only : da_to_zk, & - da_interp_lin_3d,da_interp_lin_3d_adj + da_interp_lin_3d,da_interp_lin_3d_adj,da_interp_lin_2d use da_par_util, only : da_proc_stats_combine use da_par_util1, only : da_proc_sum_int use da_statistics, only : da_stats_calculate @@ -51,5 +51,6 @@ module da_qscat #include "da_transform_xtoy_qscat.inc" #include "da_transform_xtoy_qscat_adj.inc" #include "da_calculate_grady_qscat.inc" +#include "da_qc_qscat.inc" end module da_qscat diff --git a/var/da/da_qscat/da_transform_xtoy_qscat.inc b/var/da/da_qscat/da_transform_xtoy_qscat.inc index 46893cbadd..60d37a650b 100644 --- a/var/da/da_qscat/da_transform_xtoy_qscat.inc +++ b/var/da/da_qscat/da_transform_xtoy_qscat.inc @@ -11,6 +11,10 @@ subroutine da_transform_xtoy_qscat(grid, iv, y) type (domain), intent(in) :: grid type (iv_type), intent(in) :: iv ! Innovation vector (O-B). type (y_type), intent(inout) :: y ! y = h (grid%xa) (linear) + integer :: i, j, k ! Index dimension. + real :: dx, dxm ! Interpolation weights. + real :: dy, dym ! Interpolation weights. + real :: znt, rib, ht, za, bn, fm, bd integer :: n ! Loop counter. @@ -18,6 +22,7 @@ subroutine da_transform_xtoy_qscat(grid, iv, y) real, allocatable :: v(:,:) real, allocatable :: ub(:,:) real, allocatable :: vb(:,:) + real :: v_h(kms:kme) ! Model value h at ob hor. location. if (trace_use_dull) call da_trace_entry("da_transform_xtoy_qscat") @@ -43,6 +48,38 @@ subroutine da_transform_xtoy_qscat(grid, iv, y) y%qscat(n)%u = u(1,n) y%qscat(n)%v = v(1,n) end if + if (qscat_neutral) then + i = iv%info(qscat)%i(1,n) + j = iv%info(qscat)%j(1,n) + dx = iv%info(qscat)%dx(1,n) + dy = iv%info(qscat)%dy(1,n) + dxm = iv%info(qscat)%dxm(1,n) + dym = iv%info(qscat)%dym(1,n) + do k=kts,kte + v_h(k) = dym*(dxm*grid%xb%h(i,j,k)+dx*grid%xb%h(i+1,j,k)) + dy*(dxm*grid%xb%h(i,j+1,k)+dx*grid%xb%h(i+1,j+1,k)) + end do + + znt = dym*(dxm*grid%znt(i,j) + dx*grid%znt(i+1,j)) + dy*(dxm*grid%znt(i,j+1) + dx*grid%znt(i+1,j+1)) + rib = dym*(dxm*grid%br(i,j) + dx*grid%br(i+1,j)) + dy*(dxm*grid%br(i,j+1) + dx*grid%br(i+1,j+1)) + ht = dym*(dxm*grid%ht(i,j) + dx*grid%ht(i+1,j)) + dy*(dxm*grid%ht(i,j+1) + dx*grid%ht(i+1,j+1)) + + za = v_h(kts) - ht + + bn = log( 1 + za/znt ) + if (rib < 0) then + fm=1-10 * rib / ( 1-75 * rib * bn * sqrt( 1+10/znt ) ) + else + fm = 1 / (1 + 10 * rib / sqrt(1 + 5*rib ) ) + end if + bd = bn/sqrt(fm) + + if (wind_sd_qscat) then + y%qscat(n)%u = y%qscat(n)%u * log( 1 + 10/znt ) / bd + else + y%qscat(n)%u = y%qscat(n)%u * log( 1 + 10/znt ) / bd + y%qscat(n)%v = y%qscat(n)%v * log( 1 + 10/znt ) / bd + end if + end if end do diff --git a/var/da/da_qscat/da_transform_xtoy_qscat_adj.inc b/var/da/da_qscat/da_transform_xtoy_qscat_adj.inc index e4bee88a15..f3e7c510f2 100644 --- a/var/da/da_qscat/da_transform_xtoy_qscat_adj.inc +++ b/var/da/da_qscat/da_transform_xtoy_qscat_adj.inc @@ -11,6 +11,10 @@ subroutine da_transform_xtoy_qscat_adj(grid, iv, jo_grad_y, jo_grad_x) type (iv_type), intent(in) :: iv ! obs. inc vector (o-b). type (y_type) , intent(in) :: jo_grad_y ! grad_y(jo) type (x_type) , intent(inout) :: jo_grad_x ! grad_x(jo) + integer :: i, j, k ! Index dimension. + real :: dx, dxm ! Interpolation weights. + real :: dy, dym ! Interpolation weights. + real :: znt, rib, ht, za, bn, fm, bd integer :: n ! Loop counter. @@ -18,6 +22,7 @@ subroutine da_transform_xtoy_qscat_adj(grid, iv, jo_grad_y, jo_grad_x) real, allocatable :: v(:,:) real, allocatable :: ub(:,:) real, allocatable :: vb(:,:) + real :: v_h(kms:kme) ! Model value h at ob hor. location. if (trace_use_dull) call da_trace_entry("da_transform_xtoy_qscat_adj") @@ -38,6 +43,35 @@ subroutine da_transform_xtoy_qscat_adj(grid, iv, jo_grad_y, jo_grad_x) u(1,n) = jo_grad_y%qscat(n)%u v(1,n) = jo_grad_y%qscat(n)%v end if + + if (qscat_neutral) then + i = iv%info(qscat)%i(1,n) + j = iv%info(qscat)%j(1,n) + dx = iv%info(qscat)%dx(1,n) + dy = iv%info(qscat)%dy(1,n) + dxm = iv%info(qscat)%dxm(1,n) + dym = iv%info(qscat)%dym(1,n) + do k=kts,kte + v_h(k) = dym*(dxm*grid%xb%h(i,j,k)+dx*grid%xb%h(i+1,j,k)) + dy*(dxm*grid%xb%h(i,j+1,k)+dx*grid%xb%h(i+1,j+1,k)) + end do + + znt = dym*(dxm*grid%znt(i,j) + dx*grid%znt(i+1,j)) + dy*(dxm*grid%znt(i,j+1) + dx*grid%znt(i+1,j+1)) + rib = dym*(dxm*grid%br(i,j) + dx*grid%br(i+1,j)) + dy*(dxm*grid%br(i,j+1) + dx*grid%br(i+1,j+1)) + ht = dym*(dxm*grid%ht(i,j) + dx*grid%ht(i+1,j)) + dy*(dxm*grid%ht(i,j+1) + dx*grid%ht(i+1,j+1)) + + za = v_h(kts) - ht + + bn = log( 1 + za/znt ) + if (rib < 0) then + fm=1-10 * rib / ( 1-75 * rib * bn * sqrt( 1+10/znt ) ) + else + fm = 1 / (1 + 10 * rib / sqrt(1 + 5*rib ) ) + end if + bd = bn/sqrt(fm) + + u(1,n) = u(1,n) * log( 1 + 10/znt ) / bd + v(1,n) = v(1,n) * log( 1 + 10/znt ) / bd + end if end do #ifdef A2C