Skip to content

Add rain QC and neutral wind operator for scatwind #1843

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 1 commit into
base: develop
Choose a base branch
from
Open

Add rain QC and neutral wind operator for scatwind #1843

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 2 additions & 2 deletions Registry/Registry.EM_COMMON
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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" ""
Expand Down
7 changes: 4 additions & 3 deletions Registry/Registry.EM_COMMON.var
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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"

Expand Down Expand Up @@ -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" ""

#
#---------------------------------------------------------------------------------------------------------------------------------------
Expand Down
3 changes: 3 additions & 0 deletions Registry/registry.var
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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" "" ""
Expand Down Expand Up @@ -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" "" ""
Expand Down
30 changes: 30 additions & 0 deletions var/da/da_qscat/da_get_innov_vector_qscat.inc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down Expand Up @@ -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

!------------------------------------------------------------------------
Expand Down Expand Up @@ -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)

Expand Down
69 changes: 69 additions & 0 deletions var/da/da_qscat/da_qc_qscat.inc
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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
5 changes: 3 additions & 2 deletions var/da/da_qscat/da_qscat.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
37 changes: 37 additions & 0 deletions var/da/da_qscat/da_transform_xtoy_qscat.inc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -11,13 +11,18 @@ 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.

real, allocatable :: u(:,:)
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")

Expand All @@ -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

Expand Down
34 changes: 34 additions & 0 deletions var/da/da_qscat/da_transform_xtoy_qscat_adj.inc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -11,13 +11,18 @@ 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.

real, allocatable :: u(:,:)
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")

Expand All @@ -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
Expand Down