var/da/da_radar/da_get_innov_vector_radar.inc

   1 subroutine da_get_innov_vector_radar (it, grid, ob, iv)
   2
   3    !-----------------------------------------------------------------------
   4    ! Purpose: Calculate innovation vectors for radar obs
   5    ! History:
   6    !    10/22/2008 - Updated for Analysis on Arakawa-C grid (Syed RH Rizvi, NCAR)
   7    !    ......
   8    !    08/08/2016 - Updated to include null-echo assimilation
   9    !                 (Yu-Shin Kim and Ki-Hong Min, School of Earth System
  10    !                  Sciences/Kyungpook National University, Daegu, S.Korea)
  11    !    08/2017 - bug fix for Vr operator (Siou-Ying Jiang, CWB, Taiwan)
  12    !-----------------------------------------------------------------------
  13
  14    implicit none
  15
  16    integer,          intent(in)    :: it       ! External iteration.
  17    type(domain),     intent(in)    :: grid     ! first guess state.
  18    type(y_type),     intent(inout) :: ob       ! Observation structure.
  19    type(iv_type),    intent(inout) :: iv       ! O-B structure.
  20
  21    integer :: n        ! Loop counter.
  22    integer :: i, j, k  ! Index dimension.
  23    real    :: dx, dxm  ! Interpolation weights.
  24    real    :: dy, dym  ! Interpolation weights.
  25
  26    real, allocatable :: model_p(:,:)
  27    real, allocatable :: model_u(:,:)
  28    real, allocatable :: model_v(:,:)
  29    real, allocatable :: model_w(:,:)
  30
  31    real, allocatable :: model_rv(:,:)
  32    real, allocatable :: model_rf(:,:)
  33    real, allocatable :: model_ps(:)
  34
  35    real, allocatable :: model_qv(:,:)
  36    real, allocatable :: model_qs(:,:)
  37    real, allocatable :: model_qs_ice(:,:)
  38    real, allocatable :: model_tc(:,:)
  39
  40    real, allocatable :: model_rho(:,:)
  41    real, allocatable :: model_qrn(:,:)
  42    real, allocatable :: model_qcl(:,:)
  43    real, allocatable :: model_qci(:,:)
  44    real, allocatable :: model_qsn(:,:)
  45    real, allocatable :: model_qgr(:,:)
  46
  47    real    :: v_h(kms:kme)      ! Model value h at ob hor. location.
  48
  49    real    :: xr,yr,zr
  50    integer :: irv, irvf
  51    integer :: irf, irff
  52
  53    real    :: alog_10, czrn, czds, czws, czgr, zrn, zds, zws, zgr, rze
  54    real    :: ob_radar_rf, bg_rze, bg_rf
  55    real    :: cwr, cws ! weighting coefficient for mixing ratio
  56
  57    ! variables for calculating cloud base
  58    real, allocatable :: model_lcl(:)           !model LCL at ob locations
  59    ! for opt=1
  60    real              :: h_sfc, p_sfc, t_sfc, q_sfc, ev_sfc, td_sfc
  61    ! for opt=2
  62    real              :: zlcl(ims:ime,jms:jme)
  63
  64    logical           :: echo_non_precip, echo_rf_good
  65
  66    !--------------------------------------------------------
  67    !  for background-dependent hydrmeteor retrieval scheme
  68    !--------------------------------------------------------
  69    character(len=filename_len) :: hydro_weight_file
  70    integer              :: hydro_weight_unit
  71    integer              :: tot_h_index, tot_z_index
  72    integer              :: ii, jj, kk, nk
  73    integer              :: h_index, z_index
  74    logical              :: file_exist, qg_exist
  75    real                 :: zern_ratio, zews_ratio, zeds_ratio, zegr_ratio  ! contributions of each hydrometeor to total reflectivity
  76    real, allocatable    :: num_sample(:,:)  ! number of samples from the background
  77    real, allocatable    :: avg_zern(:,:)    ! ze contributed by bin-averaged rainwater
  78    real, allocatable    :: avg_zeds(:,:)    ! ze contributed by bin-averaged dry snow
  79    real, allocatable    :: avg_zews(:,:)    ! ze contributed by bin-averaged wet snow
  80    real, allocatable    :: avg_zegr(:,:)    ! ze contributed by bin-averaged graupel
  81    real, allocatable    :: avg_qrn(:,:)     ! bin-averaged rainwater
  82    real, allocatable    :: avg_qds(:,:)     ! bin-averaged dry snow
  83    real, allocatable    :: avg_qws(:,:)     ! bin-averaged wet snow
  84    real, allocatable    :: avg_qgr(:,:)     ! bin-averaged graupel
  85    real, allocatable    :: ave_rho(:,:)     ! bin-averaged air density
  86
  87    !------------------------
  88    !  for jung et al 2008
  89    !------------------------
  90    real    :: qvp,qra,qsn,qgr ! mixing ratio
  91    real    :: dqra,dqsn,dqgr,dtmk,dqvp
  92    real    :: dqnr,dqns,dqng
  93    real    :: zmm,zmm_ref
  94    real    :: qnr,qns,qng            ! number concentration of rain snow and graupel
  95    real    :: tmk,prs                ! temperature and pressure
  96    real    :: dbz                    ! reflectivity in dBZ
  97    real    :: rn0_r,rn0_s,rn0_g      ! intercept parameter of rain snow and graupel
  98    real    :: rhos,rhog              ! density of snow and graupel
  99
 100    !------------------------
 101
 102    alog_10 = alog(10.0)
 103
 104    ! Ze=zv*(ro*v)**1.75
 105    ! Zdb=10*log10(Ze)
 106    zrn = 3.63*1.00e+9  ! rainwater
 107    zds = 9.80*1.00e+8  ! dry snow
 108    zws = 4.26*1.00e+11 ! wet snow
 109    zgr = 4.33*1.00e+10 ! grauple
 110
 111    !------------------------
 112    !  for jung et al 2008
 113    !------------------------
 114    qnr=0
 115    qns=0
 116    qng=0
 117    rn0_r=8e6
 118    rn0_s=3e6
 119    rn0_g=4e6
 120    rhos=100.0
 121    rhog=400.0
 122    !------------------------
 123    if (trace_use) call da_trace_entry("da_get_innov_vector_radar")
 124
 125    irv = 0; irvf = 0; irf = 0; irff = 0
 126
 127
 128    ! No point in going through and allocating all these variables if we're just going to quit anyway
 129
 130    if ( use_radar_rf .and. use_radar_rhv ) then
 131       write(unit=message(1),fmt='(A)') "Both 'use_radar_rf' and 'use_radar_rhv' are set to true"
 132       write(unit=message(2),fmt='(A)') "You must only choose one of these options"
 133       call da_error(__FILE__,__LINE__,message(1:2))
 134    end if
 135
 136
 137    allocate (model_p(iv%info(radar)%max_lev,iv%info(radar)%n1:iv%info(radar)%n2))
 138    allocate (model_u(iv%info(radar)%max_lev,iv%info(radar)%n1:iv%info(radar)%n2))
 139    allocate (model_v(iv%info(radar)%max_lev,iv%info(radar)%n1:iv%info(radar)%n2))
 140    allocate (model_w(iv%info(radar)%max_lev,iv%info(radar)%n1:iv%info(radar)%n2))
 141
 142    if ( use_radar_rv ) allocate (model_rv(iv%info(radar)%max_lev,iv%info(radar)%n1:iv%info(radar)%n2))
 143    if ( use_radar_rf ) allocate (model_rf(iv%info(radar)%max_lev,iv%info(radar)%n1:iv%info(radar)%n2))
 144    allocate (model_ps(iv%info(radar)%n1:iv%info(radar)%n2))
 145
 146    allocate (model_qv(iv%info(radar)%max_lev,iv%info(radar)%n1:iv%info(radar)%n2))
 147    allocate (model_qs(iv%info(radar)%max_lev,iv%info(radar)%n1:iv%info(radar)%n2))
 148    allocate (model_tc(iv%info(radar)%max_lev,iv%info(radar)%n1:iv%info(radar)%n2))
 149
 150    allocate (model_rho(iv%info(radar)%max_lev,iv%info(radar)%n1:iv%info(radar)%n2))
 151    allocate (model_qrn(iv%info(radar)%max_lev,iv%info(radar)%n1:iv%info(radar)%n2))
 152    allocate (model_qcl(iv%info(radar)%max_lev,iv%info(radar)%n1:iv%info(radar)%n2))
 153    allocate (model_qci(iv%info(radar)%max_lev,iv%info(radar)%n1:iv%info(radar)%n2))
 154    allocate (model_qsn(iv%info(radar)%max_lev,iv%info(radar)%n1:iv%info(radar)%n2))
 155    allocate (model_qgr(iv%info(radar)%max_lev,iv%info(radar)%n1:iv%info(radar)%n2))
 156
 157    if ( use_radar_rqv ) then
 158       allocate (model_lcl(iv%info(radar)%n1:iv%info(radar)%n2))
 159       allocate (model_qs_ice(iv%info(radar)%max_lev,iv%info(radar)%n1:iv%info(radar)%n2))
 160       model_qs_ice(:,:) = 0. !initialize
 161    end if
 162
 163    model_p(:,:)    = 0.
 164    model_u(:,:)    = 0.
 165    model_v(:,:)    = 0.
 166    model_w(:,:)    = 0.
 167
 168    if ( use_radar_rv ) model_rv(:,:)   = 0.
 169    if ( use_radar_rf ) model_rf(:,:)   = 0.
 170    model_ps(:)     = 0.
 171
 172    model_qv(:,:)   = 0.
 173    model_qs(:,:)   = 0.
 174    model_tc(:,:)   = 0.
 175
 176    model_rho(:,:)  = 0.
 177    model_qrn(:,:)  = 0.
 178    model_qci(:,:)  = 0.
 179    model_qcl(:,:)  = 0.
 180    model_qsn(:,:)  = 0.
 181    model_qgr(:,:)  = 0.
 182
 183    if ( it > 1 ) then
 184      do n=iv%info(radar)%n1,iv%info(radar)%n2
 185         do k=1,iv%info(radar)%levels(n)
 186            if (iv%radar(n)%rv(k)%qc == fails_error_max) iv%radar(n)%rv(k)%qc = 0
 187            if (iv%radar(n)%rf(k)%qc == fails_error_max) iv%radar(n)%rf(k)%qc = 0
 188         end do
 189      end do
 190    end if
 191
 192    do n=iv%info(radar)%n1,iv%info(radar)%n2
 193       if (iv%info(radar)%levels(n) < 1) cycle
 194
 195       ! [1.0] Get horizontal interpolation weights:
 196
 197       i   = iv%info(radar)%i(1,n)
 198       j   = iv%info(radar)%j(1,n)
 199       dx  = iv%info(radar)%dx(1,n)
 200       dy  = iv%info(radar)%dy(1,n)
 201       dxm = iv%info(radar)%dxm(1,n)
 202       dym = iv%info(radar)%dym(1,n)
 203
 204       do k=kts,kte
 205          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))
 206       end do
 207
 208       do k=1, iv%info(radar)%levels(n)
 209          call da_to_zk(iv%radar(n)%height(k), v_h, v_interp_h, iv%info(radar)%zk(k,n))
 210
 211          if (iv%info(radar)%zk(k,n) < 1.0) then
 212             iv%radar(n)%height_qc(k) = below_model_surface
 213          else if (iv%info(radar)%zk(k,n) > mkz) then
 214             iv%radar(n)%height_qc(k) = above_model_lid
 215          end if
 216       end do
 217    end do
 218
 219    call da_convert_zk (iv%info(radar))
 220
 221    ! [2.0] Interpolate horizontally to ob points:
 222
 223    call da_interp_lin_3d (grid%xb % p,   iv%info(radar), model_p)
 224 #ifdef A2C
 225    call da_interp_lin_3d (grid%xb % u,   iv%info(radar), model_u,'u')
 226    call da_interp_lin_3d (grid%xb % v,   iv%info(radar), model_v,'v')
 227 #else
 228    call da_interp_lin_3d (grid%xb % u,   iv%info(radar), model_u)
 229    call da_interp_lin_3d (grid%xb % v,   iv%info(radar), model_v)
 230 #endif
 231    call da_interp_lin_3d (grid%xb % wh,  iv%info(radar), model_w)
 232    call da_interp_lin_3d (grid%xb % rho, iv%info(radar), model_rho)
 233    call da_interp_lin_3d (grid%xb % qrn, iv%info(radar), model_qrn)
 234    call da_interp_lin_3d (grid%xb % qcw, iv%info(radar), model_qcl)
 235    call da_interp_lin_3d (grid%xb % qci, iv%info(radar), model_qci)
 236    call da_interp_lin_3d (grid%xb % qsn, iv%info(radar), model_qsn)
 237 IF ( ASSOCIATED( grid%xb%qgr ) ) THEN
 238    call da_interp_lin_3d (grid%xb % qgr, iv%info(radar), model_qgr)
 239 END IF
 240    call da_interp_lin_3d (grid%xb % q,  iv%info(radar), model_qv)
 241    call da_interp_lin_3d (grid%xb % qs, iv%info(radar), model_qs)
 242    call da_interp_lin_3d (grid%xb % t,  iv%info(radar), model_tc)
 243
 244    model_tc = model_tc - 273.15 ! degree K to degree C
 245
 246    if ( use_radar_rqv ) then
 247       do n=iv%info(radar)%n1,iv%info(radar)%n2
 248          do k=1,iv%info(radar)%levels(n)
 249             model_qs_ice(k,n) = model_qs(k,n)*exp((2.837E6 - 2.501E6)/461.5*(1./273.15 - 1./(model_tc(k,n)+273.15)))
 250          end do
 251       end do
 252    end if
 253
 254    ! calculate background/model LCL to be used by use_radar_rqv
 255    if ( use_radar_rqv .and. cloudbase_calc_opt == 2 ) then
 256       do j = jts, jte
 257          do i = its, ite
 258             zlcl(i,j) = 125.0*(grid%xb%t(i,j,1)-grid%xb%td(i,j,1))
 259          end do
 260       end do
 261    end if ! lcl for use_radar_rqv
 262
 263    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 264    ! background-dependent hydrometer retrieval scheme !
 265    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 266    if (use_radar_rhv .and. radar_rhv_opt == 2 ) then
 267       !! allocate variables
 268       tot_h_index = 40  ! from 500 m to 20 km, at an interval of 500 m
 269       tot_z_index = 7   ! from -5 dBZ to 65 dBZ, at an interval of 10 dBZ
 270       allocate (num_sample(tot_h_index,tot_z_index))
 271       allocate (avg_zern(tot_h_index,tot_z_index))
 272       allocate (avg_zeds(tot_h_index,tot_z_index))
 273       allocate (avg_zews(tot_h_index,tot_z_index))
 274       allocate (avg_zegr(tot_h_index,tot_z_index))
 275       allocate (avg_qrn(tot_h_index,tot_z_index))
 276       allocate (avg_qds(tot_h_index,tot_z_index))
 277       allocate (avg_qws(tot_h_index,tot_z_index))
 278       allocate (avg_qgr(tot_h_index,tot_z_index))
 279       allocate (ave_rho(tot_h_index,tot_z_index))
 280
 281       !! variable initialization
 282       num_sample = 0.
 283       avg_qrn    = 0.
 284       avg_qws    = 0.
 285       avg_qds    = 0.
 286       avg_qgr    = 0.
 287       ave_rho    = 0.
 288       avg_zern   = 0.
 289       avg_zews   = 0.
 290       avg_zeds   = 0.
 291       avg_zegr   = 0.
 292
 293       !! read historical statistics from hydro_mean.dat if available
 294       hydro_weight_file='hydro_mean.dat'
 295       inquire(file=trim(hydro_weight_file), exist=file_exist)
 296       if (file_exist) then
 297          call da_get_unit(hydro_weight_unit)
 298          open(unit=hydro_weight_unit, file=trim(hydro_weight_file), form='FORMATTED')
 299          read(unit=hydro_weight_unit, fmt='(A)')
 300          do z_index=1,tot_z_index
 301             do h_index=1,tot_h_index
 302                read(hydro_weight_unit, fmt='(2(10x), 4(f19.9,2x))') &
 303                     avg_zern(h_index,z_index), avg_zews(h_index,z_index), avg_zeds(h_index,z_index), avg_zegr(h_index,z_index)
 304             end do
 305          end do
 306          close(hydro_weight_unit)
 307          call da_free_unit(hydro_weight_unit)
 308       end if
 309
 310       !! calculate sum of background states in the current processor
 311       do kk=kds, kde
 312          do jj=jps, jpe
 313             do ii=ips, ipe
 314                !! calculate background reflectivity
 315                call da_radar_rf(grid%xb%qrn(ii,jj,kk), grid%xb%qsn(ii,jj,kk), grid%xb%qgr(ii,jj,kk), &
 316                                 grid%xb%t(ii,jj,kk)-273.15, grid%xb%rho(ii,jj,kk), bg_rze)
 317                bg_rf = 10.*log10(bg_rze)
 318                !! get the index of reflectvity
 319                z_index = nint(bg_rf/10.)+1
 320                z_index = max(z_index, 0) ! set to non-precip if below -5 dBZ
 321                z_index = min(z_index, tot_z_index) ! set to 65 dBZ if above
 322                !! get the height index
 323                h_index = nint(grid%xb%h(ii,jj,kk)/500.)
 324                h_index = max(h_index, 1) ! set to 500 m if below
 325                h_index = min(h_index,tot_h_index) ! set to 20 km if above
 326
 327                !! Sum of the model states of different model levels and reflectivity thresholds
 328                if (z_index .ne. 0 ) then
 329                   avg_qrn(h_index,z_index)  = avg_qrn(h_index,z_index)  + grid%xb%qrn(ii,jj,kk)
 330                   if ( grid%xb%t(ii,jj,kk) > 273.15 ) then
 331                      avg_qws(h_index,z_index) = avg_qws(h_index,z_index) + grid%xb%qsn(ii,jj,kk)
 332                   else
 333                      avg_qds(h_index,z_index) = avg_qds(h_index,z_index) + grid%xb%qsn(ii,jj,kk)
 334                   end if
 335                   avg_qgr(h_index,z_index) = avg_qgr(h_index,z_index) + grid%xb%qgr(ii,jj,kk)
 336                   ave_rho(h_index,z_index) = ave_rho(h_index,z_index) + grid%xb%rho(ii,jj,kk)
 337                   num_sample(h_index,z_index) = num_sample(h_index,z_index) + 1.
 338                end if
 339             end do ! west-east
 340          end do  ! south-north
 341       end do ! bottom-top
 342
 343       !! sum of all processors and get the averaged background states
 344       do z_index=1,tot_z_index
 345          do h_index=1,tot_h_index
 346             num_sample(h_index,z_index) = wrf_dm_sum_real(num_sample(h_index,z_index))
 347             if (num_sample(h_index,z_index) .gt. 0) then
 348                ave_rho(h_index,z_index) = wrf_dm_sum_real(ave_rho(h_index,z_index)) / num_sample(h_index,z_index)
 349                avg_qrn(h_index,z_index) = wrf_dm_sum_real(avg_qrn(h_index,z_index)) / num_sample(h_index,z_index)
 350                avg_qws(h_index,z_index) = wrf_dm_sum_real(avg_qws(h_index,z_index)) / num_sample(h_index,z_index)
 351                avg_qds(h_index,z_index) = wrf_dm_sum_real(avg_qds(h_index,z_index)) / num_sample(h_index,z_index)
 352                avg_qgr(h_index,z_index) = wrf_dm_sum_real(avg_qgr(h_index,z_index)) / num_sample(h_index,z_index)
 353             end if
 354          end do
 355       end do
 356
 357       !! calculate the contributions of each hydrometeor to total reflectivity and save them to hydro_mean.dat.update
 358       hydro_weight_file='hydro_mean.dat.update'
 359       if (rootproc) call da_get_unit(hydro_weight_unit)
 360       if (rootproc) open(unit=hydro_weight_unit, file=trim(hydro_weight_file), form='FORMATTED')
 361       if (rootproc) write(unit=hydro_weight_unit, fmt='(2(a8,2x), 4(a19,2x))') &
 362                           "z_index:", "h_index:", "===Rainwater===", "===Wet snow===", "===Dry snow===", "===Graupel==="
 363       do z_index=1,tot_z_index
 364          do h_index=1,tot_h_index
 365             if (num_sample(h_index,z_index) .gt. 10.) then
 366                if (avg_qrn(h_index,z_index) > 0.)  & !! rain water
 367                   avg_zern(h_index,z_index) =  zrn*(ave_rho(h_index,z_index)*avg_qrn(h_index,z_index))**1.75
 368                if (avg_qws(h_index,z_index) > 0.)  & !! wet snow
 369                   avg_zews(h_index,z_index) =  zws*(ave_rho(h_index,z_index)*avg_qws(h_index,z_index))**1.75
 370                if (avg_qds(h_index,z_index) > 0.)  & !! dry snow
 371                   avg_zeds(h_index,z_index) =  zds*(ave_rho(h_index,z_index)*avg_qds(h_index,z_index))**1.75
 372                if (avg_qgr(h_index,z_index) > 0.)  & !! graupel
 373                   avg_zegr(h_index,z_index) =  zgr*(ave_rho(h_index,z_index)*avg_qgr(h_index,z_index))**1.75
 374             end if
 375             if (rootproc) &
 376                write(unit=hydro_weight_unit, fmt='(2(i8, 2x), 4(f19.9,2x))')  z_index, h_index, &
 377                      avg_zern(h_index,z_index),avg_zews(h_index,z_index), avg_zeds(h_index,z_index), avg_zegr(h_index,z_index)
 378          end do
 379       end do   !bottom-top
 380       if (rootproc) close(hydro_weight_unit)
 381       if (rootproc) call da_free_unit(hydro_weight_unit)
 382    end if   !! use_radar_rhv .and. radar_rhv_opt == 2
 383
 384    do n=iv%info(radar)%n1,iv%info(radar)%n2
 385
 386       if ( use_radar_rf .and. radar_rf_opt==1) then
 387          ! for Xiao's default scheme
 388          ! Test 5.0E-8 (kg/kg) as critical value. It can not be smaller.
 389          do k=1,iv%info(radar)%levels(n)
 390             model_qrn(k,n)=amax1(5.0E-8,model_qrn(k,n))
 391          end do
 392       end if
 393
 394       i   = iv%info(radar)%i(1,n)
 395       j   = iv%info(radar)%j(1,n)
 396       dx  = iv%info(radar)%dx(1,n)
 397       dy  = iv%info(radar)%dy(1,n)
 398       dxm = iv%info(radar)%dxm(1,n)
 399       dym = iv%info(radar)%dym(1,n)
 400
 401       model_ps(n) = dxm *(dym * grid%xb % psac(i,  j) + dy * grid%xb%psac(i+1,  j)) + &
 402                  dx  *(dym * grid%xb % psac(i,j+1) + dy * grid%xb%psac(i+1,j+1)) + &
 403                  grid%xb % ptop
 404
 405       ! calculate model LCL at ob locations
 406       if ( use_radar_rqv ) then
 407          select case ( cloudbase_calc_opt )
 408          case ( 1 )  !KNU scheme
 409             ! get surface variable to calculate cloud base
 410             h_sfc = dxm *(dym * grid%xb % terr(i,  j) + dy * grid%xb% terr(i+1,  j)) + &
 411                     dx  *(dym * grid%xb % terr(i,j+1) + dy * grid%xb% terr(i+1,j+1))
 412             p_sfc = dxm *(dym * grid%xb % psfc(i,  j) + dy * grid%xb% psfc(i+1,  j)) + &
 413                     dx  *(dym * grid%xb % psfc(i,j+1) + dy * grid%xb% psfc(i+1,j+1))
 414             t_sfc = dxm *(dym * grid%xb % t2(i,  j) + dy * grid%xb% t2(i+1,  j)) + &
 415                     dx  *(dym * grid%xb % t2(i,j+1) + dy * grid%xb% t2(i+1,j+1))
 416             q_sfc = dxm *(dym * grid%xb % q2(i,  j) + dy * grid%xb% q2(i+1,  j)) + &
 417                     dx  *(dym * grid%xb % q2(i,j+1) + dy * grid%xb% q2(i+1,j+1))
 418
 419             ! calculate cloud base height based on model background variables
 420             ev_sfc = q_sfc*p_sfc/0.622
 421             td_sfc = 1./(1./273.15 - 461.51/2.501e6*log(ev_sfc/611.2))
 422             model_lcl(n) = h_sfc + 125.*(t_sfc-td_sfc)
 423             model_lcl(n) = amax1(h_sfc+500.,model_lcl(n))
 424             model_lcl(n) = amin1(3000.     ,model_lcl(n))
 425
 426          case ( 2 )  !NCAR scheme
 427             model_lcl(n) = dxm *(dym * zlcl(i,j)   + dy * zlcl(i+1, j))  + &
 428                            dx  *(dym * zlcl(i,j+1) + dy * zlcl(i+1,j+1))
 429             ! zlcl_mean is model grid mean zlcl calculated in da_transform_wrftoxb.inc
 430             model_lcl(n) = max(model_lcl(n), zlcl_mean)
 431
 432             ! add radar elevation to the calculated model_lcl
 433             ! iv%radar(n)%height(k) (rather than zr) is later used for checking
 434             model_lcl(n) = model_lcl(n) + iv%radar(n)%stn_loc%elv
 435
 436          case default
 437             model_lcl(n) = 1500.0
 438
 439          end select
 440       end if ! lcl for use_radar_rqv
 441
 442       iv%radar(n)%model_p(1:iv%info(radar)%levels(n))   = model_p(1:iv%info(radar)%levels(n),n)
 443       iv%radar(n)%model_rho(1:iv%info(radar)%levels(n)) = model_rho(1:iv%info(radar)%levels(n),n)
 444       iv%radar(n)%model_qrn(1:iv%info(radar)%levels(n)) = model_qrn(1:iv%info(radar)%levels(n),n)
 445       iv%radar(n)%model_qcl(1:iv%info(radar)%levels(n)) = model_qcl(1:iv%info(radar)%levels(n),n)
 446       iv%radar(n)%model_qci(1:iv%info(radar)%levels(n)) = model_qci(1:iv%info(radar)%levels(n),n)
 447       iv%radar(n)%model_qsn(1:iv%info(radar)%levels(n)) = model_qsn(1:iv%info(radar)%levels(n),n)
 448       iv%radar(n)%model_qgr(1:iv%info(radar)%levels(n)) = model_qgr(1:iv%info(radar)%levels(n),n)
 449
 450       iv%radar(n)%model_ps     = model_ps(n)
 451
 452       ! [3.0] Calculate rv, rf at OBS location and initialise components of &
 453       ! innovation vector:
 454
 455       if (fg_format == fg_format_wrf_arw_regional .or. &
 456           fg_format == fg_format_wrf_arw_global ) then
 457          call da_llxy_wrf(map_info, &
 458             iv%radar(n)%stn_loc%lat, iv%radar(n)%stn_loc%lon, &
 459             iv%radar(n)%stn_loc%x,   iv%radar(n)%stn_loc%y)
 460       else
 461          call da_llxy_default( iv%radar(n)%stn_loc%lat, iv%radar(n)%stn_loc%lon, &
 462             iv%radar(n)%stn_loc%x,   iv%radar(n)%stn_loc%y)
 463       end if
 464
 465       xr = grid%xb%ds *(iv%info(radar)%x(1,n) - iv%radar(n)%stn_loc%x)
 466       yr = grid%xb%ds *(iv%info(radar)%y(1,n) - iv%radar(n)%stn_loc%y)
 467
 468       level_loop: do k=1, iv%info(radar)%levels(n)
 469          iv % radar(n) % rv(k) % inv = 0.0
 470          iv % radar(n) % rf(k) % inv = 0.0
 471
 472          ! initialize
 473          echo_non_precip = .false.
 474          echo_rf_good    = .false.
 475
 476          if (iv % radar(n) % height_qc(k) /= below_model_surface .and.  &
 477              iv % radar(n) % height_qc(k) /= above_model_lid) then
 478
 479             zr = iv%radar(n)%height(k) - iv%radar(n)%stn_loc%elv
 480
 481             ! identify if non-precip obs (rf = radar_non_precip_rf)
 482             echo_non_precip = abs(ob%radar(n)%rf(k) - radar_non_precip_rf) < 0.1
 483
 484             ! identify if valid rf obs to process
 485             ! this includes the echo_non_precip case
 486             echo_rf_good = (abs(ob % radar(n) % rf(k) - missing_r) > 1.0) &
 487                            .and. (iv % radar(n) % rf(k) % qc >= obs_qc_pointer)
 488
 489             if (use_radar_rv) then
 490                !set qc to missing_data for rv of -999.99 (radar_non_precip_rf)
 491                if ( abs(ob%radar(n)%rv(k) - radar_non_precip_rf) < 0.1 ) then
 492                    iv % radar(n) % rv(k) % qc = missing_data
 493                end if
 494                if (abs(iv % radar(n) % rv(k) % qc - missing_data) > 1) then
 495                   if (abs(ob % radar(n) % rv(k) - missing_r) > 1.0 .AND. &
 496                        iv % radar(n) % rv(k) % qc >= obs_qc_pointer) then
 497
 498                      !reference: Sun and Crook (1997)
 499                      call da_radial_velocity(model_rv(k,n), model_p(k,n),  &
 500                         model_u(k,n), model_v(k,n), model_w(k,n),          &
 501                         model_qrn(k,n), model_ps(n), xr, yr, zr,           &
 502                         model_rho(k,n))
 503
 504                      iv % radar(n) % rv(k) % inv = ob % radar(n) % rv(k) - model_rv(k,n)*radar_rv_rscl
 505                   end if
 506                end if
 507             end if
 508
 509             if (use_radar_rf .and. radar_rf_opt==2) then
 510               iv % radar(n) % zmm(k) % inv = 0
 511               iv % radar(n) % rf (k) % qc = -5 ! assume bad rf obs first
 512               echo_rf_good = ob % radar(n) % rf(k) >= rfmin
 513               if ( echo_rf_good ) then
 514                   tmk=model_tc(k,n)+273.15
 515                   prs=model_p(k,n)
 516                   zmm=0
 517                   zmm_ref=0
 518                   call da_radzicevar(model_qv(k,n),model_qrn(k,n),model_qsn(k,n),model_qgr(k,n),qnr,qns,qng,tmk,prs,dbz,        &
 519                                        0,0,0,rn0_r,rn0_s,rn0_g,                                                                   &
 520                                        rhos,rhog,dtmk,dqvp,dqra,dqsn,dqgr,dqnr,dqns,dqng,zmm,0,                                   &
 521                                        0,zmm_ref)
 522                   model_rf(k,n)=dbz*radar_rf_rscl
 523                   iv % radar(n) % zmm(k) % inv = zmm
 524                   iv % radar(n) % rf(k) % inv = ob % radar(n) % rf(k) - model_rf(k,n)
 525                   if (model_rf(k,n) >= rfmin) iv % radar(n) % rf (k) % qc  = 0
 526               end if
 527             end if
 528
 529             if (use_radar_rf.and.radar_rf_opt==1) then
 530                if ( echo_rf_good ) then
 531                   model_rf(k,n) = leh1 + leh2 * alog10(model_rho(k,n) * model_qrn(k,n) * 1.0e+3)
 532                   iv % radar(n) % rf(k) % inv = ob % radar(n) % rf(k) - model_rf(k,n)
 533                end if
 534             end if
 535
 536             ! Calculate background/model reflectivity
 537             if (use_radar_rhv .or. use_radar_rqv) then
 538                if ( echo_rf_good ) then
 539                   call da_radar_rf (model_qrn(k,n),model_qsn(k,n),model_qgr(k,n),model_tc(k,n),model_rho(k,n),bg_rze)
 540                   bg_rf = 10.0*log10(bg_rze)  ! Z to dBZ
 541                   iv % radar(n) % zmm(k) % inv = bg_rf
 542                   iv % radar(n) % rf(k) % inv = ob % radar(n) % rf(k) - bg_rf
 543                end if
 544             end if
 545
 546             ! Calculate retrieved hydrometeorological variables
 547             ! Background-dependent retrieval scheme (Chen et al. 2020 AR; Chen et al. 2021 QJRMS)
 548             if (use_radar_rhv) then
 549                if ( echo_rf_good ) then
 550
 551                   ! compute retrieved hydrometeors rhv
 552                   if (it.eq.1) then
 553
 554                      ob_radar_rf = ob % radar(n) % rf(k)
 555
 556                      if ( radar_non_precip_opt > 0 ) then ! assimilate non_precip echo
 557                         if ( echo_non_precip ) then ! ob is non-precip
 558                            if ( bg_rf > -15.0 ) then
 559                               ! when background/model is precip
 560                               ! retrieve hydrometeors from -15dBZ
 561                               ob_radar_rf = -15.0
 562                            else
 563                               ! when background/model is also non-precip
 564                               ! do not need to assimilate this ob
 565                               iv % radar(n) % rf (k) % qc = -5
 566                               iv % radar(n) % rrn(k) % qc = -5
 567                               iv % radar(n) % rsn(k) % qc = -5
 568                               iv % radar(n) % rgr(k) % qc = -5
 569                               cycle level_loop
 570                            end if
 571                         end if !if echo_non_precip
 572                      end if
 573
 574                      ! The original WRFDA hydrometeor retrieval scheme
 575                      if (model_tc(k,n) .ge. 5.0) then
 576                         czrn = 1.0
 577                         czws = 0.0
 578                         czds = 0.0
 579                         czgr = 0.0
 580                      else if (model_tc(k,n) .ge.  0.0) then
 581                         czrn = (model_tc(k,n)+5.0)/10.0
 582                         czws = (1.0-czrn)*zws/(zws+zgr)
 583                         czds = 0.0
 584                         czgr = (1.0-czrn)*zgr/(zws+zgr)
 585                      else if (model_tc(k,n) .ge. -5.0) then
 586                         czrn = (model_tc(k,n)+5.0)/10.0
 587                         czws = 0.0
 588                         czds = (1.0-czrn)*zds/(zds+zgr)
 589                         czgr = (1.0-czrn)*zgr/(zds+zgr)
 590                      else if (model_tc(k,n) .lt. -5.0) then
 591                         czrn = 0.0
 592                         czws = 0.0
 593                         czds = zds/(zds+zgr)
 594                         czgr = zgr/(zds+zgr)
 595                      end if
 596
 597                      if (radar_rhv_opt == 2) then
 598                        ! backgound-dependent reflectivity retrival scheme (Chen et al. 2020, AR; Chen et al. 2021, QJRMS)
 599                        !! get the index of reflectvity
 600                        z_index = nint(ob_radar_rf/10.+1)
 601                        z_index = max(z_index, 0)
 602                        z_index = min(z_index, tot_z_index)
 603                        !! get the height index
 604                        h_index = nint(iv%radar(n)%height(k)/500.)
 605                        h_index = max(h_index, 1)
 606                        h_index = min(h_index, tot_h_index)
 607
 608                        if (z_index > 0) then
 609                           zern_ratio = avg_zern(h_index, z_index)
 610                           zews_ratio = avg_zews(h_index, z_index)
 611                           zeds_ratio = avg_zeds(h_index, z_index)
 612                           zegr_ratio = avg_zegr(h_index, z_index)
 613                           ! detect whether rain/snow/graupel exists in certain temperatures.
 614                           qg_exist = .true.
 615                           ! when T < 273.15K
 616                           if (model_tc(k,n) .lt. -5.0) zern_ratio = 0.
 617                           if (model_tc(k,n) .lt.  0.0) zews_ratio = 0.
 618                           ! when T >= 273.15K
 619                           if (model_tc(k,n) .ge. 0.0) then
 620                              zeds_ratio = 0.
 621                              qg_exist = .false.
 622                              do nk = k, iv%info(radar)%levels(n)
 623                                 if (model_tc(nk,n) .lt. -5.0 .and. ob % radar(n) % rf(nk) .ge. 40.) qg_exist = .true.
 624                              end do
 625                           end if
 626                           if (model_tc(k,n) .ge. 5.0) zews_ratio = 0.
 627                           if (.not. qg_exist .or. model_tc(k,n) .ge. 10.0) zegr_ratio = 0.
 628
 629                           ! determine the contributions of each hydrometeor to reflectivity
 630                           if ((zern_ratio+zews_ratio+zeds_ratio+zegr_ratio) .gt. 0.) then
 631                              czrn = zern_ratio/(zern_ratio+zews_ratio+zeds_ratio+zegr_ratio)
 632                              czws = zews_ratio/(zern_ratio+zews_ratio+zeds_ratio+zegr_ratio)
 633                              czds = zeds_ratio/(zern_ratio+zews_ratio+zeds_ratio+zegr_ratio)
 634                              czgr = zegr_ratio/(zern_ratio+zews_ratio+zeds_ratio+zegr_ratio)
 635                           end if
 636                         else
 637                           ob_radar_rf = -15.0 !! Assign reflectivity below -5.0 dBZ to -15.0 dbZ for suppression
 638                                               !! No need to tune the weights because of very small impacts
 639                         end if
 640                      end if
 641
 642                      ! convert dBZ to Z
 643                      ob_radar_rf = min(ob_radar_rf, 65.0) ! if dBZ>65.0, set to 65.0
 644                      rze = 10.0**(ob_radar_rf*0.1) ! dBZ to Z
 645
 646                      ! Rainwater mixing ratio
 647                      if (czrn .gt. 0.) then
 648                         iv % radar(n) % rrno(k) = exp ( log(czrn*rze/zrn)/1.75 )/model_rho(k,n)
 649                         iv % radar(n) % rrn(k) % qc = 0
 650                      end if
 651
 652                      ! Snow mixing ratio
 653                      if ((czws+czds) .gt. 0.) then
 654                         if (model_tc(k,n) .gt. 0.) then
 655                           iv % radar(n) % rsno(k) = exp ( log(czws*rze/zws)/1.75 )/model_rho(k,n)
 656                           iv % radar(n) % rsn(k) % qc = 0
 657                         else
 658                           iv % radar(n) % rsno(k) = exp ( log(czds*rze/zds)/1.75 )/model_rho(k,n)
 659                           iv % radar(n) % rsn(k) % qc = 0
 660                         end if
 661                      end if
 662
 663                      ! Graupel mixing ratio
 664                      if (czgr .gt. 0.) then
 665                         iv % radar(n) % rgro(k) = exp ( log(czgr*rze/zgr)/1.75 )/model_rho(k,n)
 666                         iv % radar(n) % rgr(k) % qc = 0
 667                      end if
 668
 669                      if ( radar_rhv_err_opt == 1 ) then
 670                         ! rainwater error
 671                         iv % radar(n) % rrn(k) % error = iv % radar(n) % rf(k) % error * iv % radar(n) % rrno(k) * alog_10/leh2
 672                         iv % radar(n) % rrn(k) % error = amax1(0.0005,iv % radar(n) % rrn(k) % error)
 673                         iv % radar(n) % rrn(k) % error = amin1( 0.001,iv % radar(n) % rrn(k) % error)
 674                         ! snow error
 675                         iv % radar(n) % rsn(k) % error = iv % radar(n) % rf(k) % error * iv % radar(n) % rsno(k) * alog_10/leh2
 676                         iv % radar(n) % rsn(k) % error = amax1(0.0005,iv % radar(n) % rsn(k) % error)
 677                         iv % radar(n) % rsn(k) % error = amin1( 0.001,iv % radar(n) % rsn(k) % error)
 678                         ! graupel error
 679                         iv % radar(n) % rgr(k) % error = iv % radar(n) % rf(k) % error * iv % radar(n) % rgro(k) * alog_10/leh2
 680                         iv % radar(n) % rgr(k) % error = amax1(0.0005,iv % radar(n) % rgr(k) % error)
 681                         iv % radar(n) % rgr(k) % error = amin1( 0.001,iv % radar(n) % rgr(k) % error)
 682                      else if ( radar_rhv_err_opt == 2 ) then
 683                         ! use settings in the namelist
 684                         ! rainwater error
 685                         iv % radar(n) % rrn(k) % error = radar_rhv_rrn_err * 0.001 ! g/kg to kg/kg
 686                         ! snow error
 687                         iv % radar(n) % rsn(k) % error = radar_rhv_rsn_err * 0.001 ! g/kg to kg/kg
 688                         ! graupel error
 689                         iv % radar(n) % rgr(k) % error = radar_rhv_rgr_err * 0.001 ! g/kg to kg/kg
 690                      else
 691                         write(unit=message(1),fmt='(A)') "radar_rhv_err_opt should be 1 or 2"
 692                         call da_error(__FILE__,__LINE__,message(1:1))
 693                      end if
 694
 695                   end if  ! it=1
 696
 697                   if (iv % radar(n) % rrn(k) % qc >= obs_qc_pointer) then
 698                      iv % radar(n) % rrn(k) % inv = iv % radar(n) % rrno(k) - model_qrn(k,n)
 699                   end if
 700                   if (iv % radar(n) % rsn(k) % qc >= obs_qc_pointer) then
 701                      iv % radar(n) % rsn(k) % inv = iv % radar(n) % rsno(k) - model_qsn(k,n)
 702                   end if
 703                   if (iv % radar(n) % rgr(k) % qc >= obs_qc_pointer) then
 704                      iv % radar(n) % rgr(k) % inv = iv % radar(n) % rgro(k) - model_qgr(k,n)
 705                   end if
 706
 707                end if ! echo_rf_good check
 708             end if ! rhv
 709
 710             ! retrieved water vapor
 711             if (use_radar_rqv) then
 712                !iv%%rqv and iv%%rqvo were assigned to missing values in read_obs_radar_ascii.inc
 713                !iter=1, rqv is missing; for second loop, dont change rqv value
 714                if (it .eq. 1) then
 715                   iv % radar(n) % rqvo(k) = 1.0*model_qs(k,n)
 716                   iv % radar(n) % rqv(k) % error = amax1(0.0005,0.1*iv % radar(n) % rqvo(k))
 717                end if
 718
 719                if ( echo_rf_good ) then
 720
 721                   ! initialize as not-assimilating rqv
 722                   iv % radar(n) % rqv(k) % qc = -5
 723
 724                   if ( echo_non_precip ) then ! ob is non-precip
 725                      if ( radar_non_precip_opt > 0 ) then ! assimilate non_precip echo
 726
 727                         if ( bg_rf >= 20.0 .and. iv%radar(n)%height(k) > model_lcl(n) ) then
 728                            iv % radar(n) % rqv(k) % qc = 0
 729
 730                            if ( model_tc(k,n) > 5.0 ) then
 731                               iv%radar(n)%rqvo(k) =  0.01*radar_non_precip_rh_w*model_qs(k,n)
 732                            else if ( model_tc(k,n) < -5.0 ) then
 733                               iv%radar(n)%rqvo(k) =  0.01*radar_non_precip_rh_i*model_qs_ice(k,n)
 734                            else
 735                               cwr = (model_tc(k,n)+5.0)/10.0
 736                               cws = 1.0 - cwr
 737                               iv%radar(n)%rqvo(k) =  cwr*0.01*radar_non_precip_rh_w*model_qs(k,n) + cws*0.01*radar_non_precip_rh_i*model_qs_ice(k,n)
 738                            end if
 739
 740                            iv % radar(n) % rqv(k) % inv = iv % radar(n) % rqvo(k) - model_qv(k,n)
 741                            ! iv % radar(n) % rqv(k) % inv must be <= 0.0
 742                            iv % radar(n) % rqv(k) % inv = amin1(0.0,iv % radar(n) % rqv(k) % inv )
 743                            iv % radar(n) % rqv(k) % error = amax1(0.0005,0.10*model_qs(k,n))
 744                         end if
 745                      end if
 746
 747                   else  ! ob is precip
 748                      if ( ob % radar(n) % rf(k) >= radar_saturated_rf  .and. iv%radar(n)%height(k) >= model_lcl(n) )then
 749                         iv % radar(n) % rqv(k) % qc = 0
 750
 751                         if ( radar_rqv_h_lbound >= 0.0 .and. radar_rqv_h_ubound >= 0.0 ) then
 752                            ! height ranges for applying retrieved rqv are set in the namelist
 753                            if ( iv%radar(n)%height(k) < radar_rqv_h_lbound .or. &
 754                                 iv%radar(n)%height(k) > radar_rqv_h_ubound ) then
 755                               ! do not assimilate rqv outside the specified height ranges
 756                               iv % radar(n) % rqv(k) % qc = -5
 757                            end if
 758                         end if
 759                      end if
 760
 761                      if (iv % radar(n) % rqv(k) % qc >= obs_qc_pointer) then
 762                         ! reduce rqvo for certain rf ranges
 763                         if ( radar_rqv_thresh1 >= 0.0 ) then
 764                            if ( ob % radar(n) % rf(k) >= radar_saturated_rf .and. &
 765                                 ob % radar(n) % rf(k) < radar_rqv_thresh1 ) then
 766                               if ( radar_rqv_rh1 >= 0.0 ) then ! RH in percentage
 767                                  iv % radar(n) % rqvo(k)= 0.01*radar_rqv_rh1*model_qs(k,n)
 768                               end if
 769                            end if
 770                         end if
 771                         if ( radar_rqv_thresh1 >= 0.0 .and. radar_rqv_thresh2 >= 0.0 .and. &
 772                              radar_rqv_thresh2 >= radar_rqv_thresh1 ) then
 773                            if ( ob % radar(n) % rf(k) >= radar_rqv_thresh1 .and. &
 774                                 ob % radar(n) % rf(k) < radar_rqv_thresh2) then
 775                               if ( radar_rqv_rh2 >= 0.0 ) then ! RH in percentage
 776                                  iv % radar(n) % rqvo(k)= 0.01*radar_rqv_rh2*model_qs(k,n)
 777                               end if
 778                            end if
 779                         end if
 780                         iv % radar(n) % rqv(k) % inv = iv % radar(n) % rqvo(k) - model_qv(k,n)
 781                         ! iv % radar(n) % rqv(k) % inv must be >= 0.0
 782                         iv % radar(n) % rqv(k) % inv = amax1(0.0,iv % radar(n) % rqv(k) % inv )
 783                         iv % radar(n) % rqv(k) % error = amax1(0.001,0.20*iv % radar(n) % rqvo(k))
 784                      end if
 785
 786                   end if ! echo_non_precip check
 787
 788                end if  ! echo_rf_good check
 789             end if  ! use_radar_rqv
 790
 791          end if  ! not surface or model lid
 792       end do level_loop
 793    end do
 794
 795    !------------------------------------------------------------------------
 796    ! [4.0] Perform optional maximum error check:
 797    !------------------------------------------------------------------------
 798
 799    if (check_max_iv)  then
 800       call da_check_max_iv_radar(iv, it, irv, irf, irvf, irff)
 801    end if
 802
 803    if (rootproc .and. check_max_iv_print) then
 804       write(unit = check_max_iv_unit, fmt ='(/,A,i5,A)')   &
 805          'For outer iteration ', it, ', Total Rejections for radar follows:'
 806
 807       if (use_radar_rv) then
 808           write( unit = check_max_iv_unit, fmt = '(/,2(A,I6))') &
 809             'Number of failed rv observations:     ',irvf, ' on ',irv
 810       end if
 811
 812       if (use_radar_rf) then
 813          write( unit = check_max_iv_unit, fmt = '(/,2(A,I6))') &
 814             'Number of failed rf observations:     ',irff, ' on ',irf
 815       end if
 816    end if
 817
 818    deallocate (model_p)
 819    deallocate (model_u)
 820    deallocate (model_v)
 821    deallocate (model_w)
 822
 823    if ( allocated(model_rv) ) deallocate (model_rv)
 824    if ( allocated(model_rf) ) deallocate (model_rf)
 825    deallocate (model_ps)
 826
 827    deallocate (model_qv)
 828    deallocate (model_qs)
 829    deallocate (model_tc)
 830
 831    deallocate (model_qrn)
 832    deallocate (model_rho)
 833    deallocate (model_qcl)
 834    deallocate (model_qci)
 835    deallocate (model_qsn)
 836    deallocate (model_qgr)
 837
 838    if ( allocated(num_sample) )  deallocate (num_sample)
 839    if ( allocated(avg_zern) )    deallocate (avg_zern)
 840    if ( allocated(avg_zeds) )    deallocate (avg_zeds)
 841    if ( allocated(avg_zews) )    deallocate (avg_zews)
 842    if ( allocated(avg_zegr) )    deallocate (avg_zegr)
 843    if ( allocated(avg_qrn) )     deallocate (avg_qrn)
 844    if ( allocated(avg_qds) )     deallocate (avg_qds)
 845    if ( allocated(avg_qws) )     deallocate (avg_qws)
 846    if ( allocated(avg_qgr) )     deallocate (avg_qgr)
 847    if ( allocated(ave_rho) )     deallocate (ave_rho)
 848
 849    if ( use_radar_rqv ) then
 850       deallocate (model_lcl)
 851       deallocate (model_qs_ice)
 852    end if
 853
 854    if (trace_use) call da_trace_exit("da_get_innov_vector_radar")
 855
 856 end subroutine da_get_innov_vector_radar
 857
 858