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, czr,czs,czg, zrr,zds,zws,zg,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 jung et al 2008
  68    !------------------------
  69
  70    real    :: qvp,qra,qsn,qgr ! mixing ratio
  71    real    :: dqra,dqsn,dqgr,dtmk,dqvp
  72    real    :: dqnr,dqns,dqng
  73    real    :: zmm,zmm_ref
  74    real    :: qnr,qns,qng            ! number concentration of rain snow and graupel
  75    real    :: tmk,prs                ! temperature and pressure
  76    real    :: dbz                    ! reflectivity in dBZ
  77    real    :: rn0_r,rn0_s,rn0_g      ! intercept parameter of rain snow and graupel
  78    real    :: rhos,rhog              ! density of snow and graupel
  79
  80    !------------------------
  81
  82    alog_10 = alog(10.0)
  83
  84    ! Ze=zv*(ro*v)**1.75
  85    ! Zdb=10*log10(Ze)
  86    zrr = 3.63*1.00e+9  ! rainwater
  87    zds = 9.80*1.00e+8  ! dry snow
  88    zws = 4.26*1.00e+11 ! wet snow
  89    zg  = 4.33*1.00e+10 ! grauple
  90
  91    !------------------------
  92    !  for jung et al 2008
  93    !------------------------
  94    qnr=0
  95    qns=0
  96    qng=0
  97    rn0_r=8e6
  98    rn0_s=3e6
  99    rn0_g=4e6
 100    rhos=100.0
 101    rhog=400.0
 102    !------------------------
 103    if (trace_use) call da_trace_entry("da_get_innov_vector_radar")
 104
 105    irv = 0; irvf = 0; irf = 0; irff = 0
 106
 107
 108    ! No point in going through and allocating all these variables if we're just going to quit anyway
 109
 110    if ( use_radar_rf .and. use_radar_rhv ) then
 111       write(unit=message(1),fmt='(A)') "Both 'use_radar_rf' and 'use_radar_rhv' are set to true"
 112       write(unit=message(2),fmt='(A)') "You must only choose one of these options"
 113       call da_error(__FILE__,__LINE__,message(1:2))
 114    end if
 115
 116
 117    allocate (model_p(iv%info(radar)%max_lev,iv%info(radar)%n1:iv%info(radar)%n2))
 118    allocate (model_u(iv%info(radar)%max_lev,iv%info(radar)%n1:iv%info(radar)%n2))
 119    allocate (model_v(iv%info(radar)%max_lev,iv%info(radar)%n1:iv%info(radar)%n2))
 120    allocate (model_w(iv%info(radar)%max_lev,iv%info(radar)%n1:iv%info(radar)%n2))
 121
 122    if ( use_radar_rv ) allocate (model_rv(iv%info(radar)%max_lev,iv%info(radar)%n1:iv%info(radar)%n2))
 123    if ( use_radar_rf ) allocate (model_rf(iv%info(radar)%max_lev,iv%info(radar)%n1:iv%info(radar)%n2))
 124    allocate (model_ps(iv%info(radar)%n1:iv%info(radar)%n2))
 125
 126    allocate (model_qv(iv%info(radar)%max_lev,iv%info(radar)%n1:iv%info(radar)%n2))
 127    allocate (model_qs(iv%info(radar)%max_lev,iv%info(radar)%n1:iv%info(radar)%n2))
 128    allocate (model_tc(iv%info(radar)%max_lev,iv%info(radar)%n1:iv%info(radar)%n2))
 129
 130    allocate (model_rho(iv%info(radar)%max_lev,iv%info(radar)%n1:iv%info(radar)%n2))
 131    allocate (model_qrn(iv%info(radar)%max_lev,iv%info(radar)%n1:iv%info(radar)%n2))
 132    allocate (model_qcl(iv%info(radar)%max_lev,iv%info(radar)%n1:iv%info(radar)%n2))
 133    allocate (model_qci(iv%info(radar)%max_lev,iv%info(radar)%n1:iv%info(radar)%n2))
 134    allocate (model_qsn(iv%info(radar)%max_lev,iv%info(radar)%n1:iv%info(radar)%n2))
 135    allocate (model_qgr(iv%info(radar)%max_lev,iv%info(radar)%n1:iv%info(radar)%n2))
 136
 137    if ( use_radar_rqv ) then
 138       allocate (model_lcl(iv%info(radar)%n1:iv%info(radar)%n2))
 139       allocate (model_qs_ice(iv%info(radar)%max_lev,iv%info(radar)%n1:iv%info(radar)%n2))
 140       model_qs_ice(:,:) = 0. !initialize
 141    end if
 142
 143    model_p(:,:)    = 0.
 144    model_u(:,:)    = 0.
 145    model_v(:,:)    = 0.
 146    model_w(:,:)    = 0.
 147
 148    if ( use_radar_rv ) model_rv(:,:)   = 0.
 149    if ( use_radar_rf ) model_rf(:,:)   = 0.
 150    model_ps(:)     = 0.
 151
 152    model_qv(:,:)   = 0.
 153    model_qs(:,:)   = 0.
 154    model_tc(:,:)   = 0.
 155
 156    model_rho(:,:)  = 0.
 157    model_qrn(:,:)  = 0.
 158    model_qci(:,:)  = 0.
 159    model_qcl(:,:)  = 0.
 160    model_qsn(:,:)  = 0.
 161    model_qgr(:,:)  = 0.
 162
 163    if ( it > 1 ) then
 164      do n=iv%info(radar)%n1,iv%info(radar)%n2
 165         do k=1,iv%info(radar)%levels(n)
 166            if (iv%radar(n)%rv(k)%qc == fails_error_max) iv%radar(n)%rv(k)%qc = 0
 167            if (iv%radar(n)%rf(k)%qc == fails_error_max) iv%radar(n)%rf(k)%qc = 0
 168         end do
 169      end do
 170    end if
 171
 172    do n=iv%info(radar)%n1,iv%info(radar)%n2
 173       if (iv%info(radar)%levels(n) < 1) cycle
 174
 175       ! [1.0] Get horizontal interpolation weights:
 176
 177       i   = iv%info(radar)%i(1,n)
 178       j   = iv%info(radar)%j(1,n)
 179       dx  = iv%info(radar)%dx(1,n)
 180       dy  = iv%info(radar)%dy(1,n)
 181       dxm = iv%info(radar)%dxm(1,n)
 182       dym = iv%info(radar)%dym(1,n)
 183
 184       do k=kts,kte
 185          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))
 186       end do
 187
 188       do k=1, iv%info(radar)%levels(n)
 189          call da_to_zk(iv%radar(n)%height(k), v_h, v_interp_h, iv%info(radar)%zk(k,n))
 190
 191          if (iv%info(radar)%zk(k,n) < 1.0) then
 192             iv%radar(n)%height_qc(k) = below_model_surface
 193          else if (iv%info(radar)%zk(k,n) > mkz) then
 194             iv%radar(n)%height_qc(k) = above_model_lid
 195          end if
 196       end do
 197    end do
 198
 199    call da_convert_zk (iv%info(radar))
 200
 201    ! [2.0] Interpolate horizontally to ob points:
 202
 203    call da_interp_lin_3d (grid%xb % p,   iv%info(radar), model_p)
 204 #ifdef A2C
 205    call da_interp_lin_3d (grid%xb % u,   iv%info(radar), model_u,'u')
 206    call da_interp_lin_3d (grid%xb % v,   iv%info(radar), model_v,'v')
 207 #else
 208    call da_interp_lin_3d (grid%xb % u,   iv%info(radar), model_u)
 209    call da_interp_lin_3d (grid%xb % v,   iv%info(radar), model_v)
 210 #endif
 211    call da_interp_lin_3d (grid%xb % wh,  iv%info(radar), model_w)
 212    call da_interp_lin_3d (grid%xb % rho, iv%info(radar), model_rho)
 213    call da_interp_lin_3d (grid%xb % qrn, iv%info(radar), model_qrn)
 214    call da_interp_lin_3d (grid%xb % qcw, iv%info(radar), model_qcl)
 215    call da_interp_lin_3d (grid%xb % qci, iv%info(radar), model_qci)
 216    call da_interp_lin_3d (grid%xb % qsn, iv%info(radar), model_qsn)
 217 IF ( ASSOCIATED( grid%xb%qgr ) ) THEN
 218    call da_interp_lin_3d (grid%xb % qgr, iv%info(radar), model_qgr)
 219 END IF
 220    call da_interp_lin_3d (grid%xb % q,  iv%info(radar), model_qv)
 221    call da_interp_lin_3d (grid%xb % qs, iv%info(radar), model_qs)
 222    call da_interp_lin_3d (grid%xb % t,  iv%info(radar), model_tc)
 223
 224    model_tc = model_tc - 273.15 ! degree K to degree C
 225
 226    if ( use_radar_rqv ) then
 227       do n=iv%info(radar)%n1,iv%info(radar)%n2
 228          do k=1,iv%info(radar)%levels(n)
 229             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)))
 230          end do
 231       end do
 232    end if
 233
 234    ! calculate background/model LCL to be used by use_radar_rqv
 235    if ( use_radar_rqv .and. cloudbase_calc_opt == 2 ) then
 236       do j = jts, jte
 237          do i = its, ite
 238             zlcl(i,j) = 125.0*(grid%xb%t(i,j,1)-grid%xb%td(i,j,1))
 239          end do
 240       end do
 241    end if ! lcl for use_radar_rqv
 242
 243    do n=iv%info(radar)%n1,iv%info(radar)%n2
 244
 245       if ( use_radar_rf .and. radar_rf_opt==1) then
 246          ! for Xiao's default scheme
 247          ! Test 5.0E-8 (kg/kg) as critical value. It can not be smaller.
 248          do k=1,iv%info(radar)%levels(n)
 249             model_qrn(k,n)=amax1(5.0E-8,model_qrn(k,n))
 250          end do
 251       end if
 252
 253       i   = iv%info(radar)%i(1,n)
 254       j   = iv%info(radar)%j(1,n)
 255       dx  = iv%info(radar)%dx(1,n)
 256       dy  = iv%info(radar)%dy(1,n)
 257       dxm = iv%info(radar)%dxm(1,n)
 258       dym = iv%info(radar)%dym(1,n)
 259
 260
 261       model_ps(n) = dxm *(dym * grid%xb % psac(i,  j) + dy * grid%xb%psac(i+1,  j)) + &
 262                  dx  *(dym * grid%xb % psac(i,j+1) + dy * grid%xb%psac(i+1,j+1)) + &
 263                  grid%xb % ptop
 264
 265       ! calculate model LCL at ob locations
 266       if ( use_radar_rqv ) then
 267          select case ( cloudbase_calc_opt )
 268          case ( 1 )  !KNU scheme
 269             ! get surface variable to calculate cloud base
 270             h_sfc = dxm *(dym * grid%xb % terr(i,  j) + dy * grid%xb% terr(i+1,  j)) + &
 271                     dx  *(dym * grid%xb % terr(i,j+1) + dy * grid%xb% terr(i+1,j+1))
 272             p_sfc = dxm *(dym * grid%xb % psfc(i,  j) + dy * grid%xb% psfc(i+1,  j)) + &
 273                     dx  *(dym * grid%xb % psfc(i,j+1) + dy * grid%xb% psfc(i+1,j+1))
 274             t_sfc = dxm *(dym * grid%xb % t2(i,  j) + dy * grid%xb% t2(i+1,  j)) + &
 275                     dx  *(dym * grid%xb % t2(i,j+1) + dy * grid%xb% t2(i+1,j+1))
 276             q_sfc = dxm *(dym * grid%xb % q2(i,  j) + dy * grid%xb% q2(i+1,  j)) + &
 277                     dx  *(dym * grid%xb % q2(i,j+1) + dy * grid%xb% q2(i+1,j+1))
 278
 279             ! calculate cloud base height based on model background variables
 280             ev_sfc = q_sfc*p_sfc/0.622
 281             td_sfc = 1./(1./273.15 - 461.51/2.501e6*log(ev_sfc/611.2))
 282             model_lcl(n) = h_sfc + 125.*(t_sfc-td_sfc)
 283             model_lcl(n) = amax1(h_sfc+500.,model_lcl(n))
 284             model_lcl(n) = amin1(3000.     ,model_lcl(n))
 285
 286          case ( 2 )  !NCAR scheme
 287             model_lcl(n) = dxm *(dym * zlcl(i,j)   + dy * zlcl(i+1, j))  + &
 288                            dx  *(dym * zlcl(i,j+1) + dy * zlcl(i+1,j+1))
 289             ! zlcl_mean is model grid mean zlcl calculated in da_transform_wrftoxb.inc
 290             model_lcl(n) = max(model_lcl(n), zlcl_mean)
 291
 292             ! add radar elevation to the calculated model_lcl
 293             ! iv%radar(n)%height(k) (rather than zr) is later used for checking
 294             model_lcl(n) = model_lcl(n) + iv%radar(n)%stn_loc%elv
 295
 296          case default
 297             model_lcl(n) = 1500.0
 298
 299          end select
 300       end if ! lcl for use_radar_rqv
 301
 302       iv%radar(n)%model_p(1:iv%info(radar)%levels(n))   = model_p(1:iv%info(radar)%levels(n),n)
 303       iv%radar(n)%model_rho(1:iv%info(radar)%levels(n)) = model_rho(1:iv%info(radar)%levels(n),n)
 304       iv%radar(n)%model_qrn(1:iv%info(radar)%levels(n)) = model_qrn(1:iv%info(radar)%levels(n),n)
 305       iv%radar(n)%model_qcl(1:iv%info(radar)%levels(n)) = model_qcl(1:iv%info(radar)%levels(n),n)
 306       iv%radar(n)%model_qci(1:iv%info(radar)%levels(n)) = model_qci(1:iv%info(radar)%levels(n),n)
 307       iv%radar(n)%model_qsn(1:iv%info(radar)%levels(n)) = model_qsn(1:iv%info(radar)%levels(n),n)
 308       iv%radar(n)%model_qgr(1:iv%info(radar)%levels(n)) = model_qgr(1:iv%info(radar)%levels(n),n)
 309
 310       iv%radar(n)%model_ps     = model_ps(n)
 311
 312       ! [3.0] Calculate rv, rf at OBS location and initialise components of &
 313       ! innovation vector:
 314
 315       if (fg_format == fg_format_wrf_arw_regional .or. &
 316           fg_format == fg_format_wrf_arw_global ) then
 317          call da_llxy_wrf(map_info, &
 318             iv%radar(n)%stn_loc%lat, iv%radar(n)%stn_loc%lon, &
 319             iv%radar(n)%stn_loc%x,   iv%radar(n)%stn_loc%y)
 320       else
 321          call da_llxy_default( iv%radar(n)%stn_loc%lat, iv%radar(n)%stn_loc%lon, &
 322             iv%radar(n)%stn_loc%x,   iv%radar(n)%stn_loc%y)
 323       end if
 324
 325       xr = grid%xb%ds *(iv%info(radar)%x(1,n) - iv%radar(n)%stn_loc%x)
 326       yr = grid%xb%ds *(iv%info(radar)%y(1,n) - iv%radar(n)%stn_loc%y)
 327
 328       level_loop: do k=1, iv%info(radar)%levels(n)
 329          iv % radar(n) % rv(k) % inv = 0.0
 330          iv % radar(n) % rf(k) % inv = 0.0
 331
 332          ! initialize
 333          echo_non_precip = .false.
 334          echo_rf_good    = .false.
 335
 336          if (iv % radar(n) % height_qc(k) /= below_model_surface .and.  &
 337              iv % radar(n) % height_qc(k) /= above_model_lid) then
 338
 339             zr = iv%radar(n)%height(k) - iv%radar(n)%stn_loc%elv
 340
 341             ! identify if non-precip obs (rf = radar_non_precip_rf)
 342             echo_non_precip = abs(ob%radar(n)%rf(k) - radar_non_precip_rf) < 0.1
 343
 344             ! identify if valid rf obs to process
 345             ! this includes the echo_non_precip case
 346             echo_rf_good = (abs(ob % radar(n) % rf(k) - missing_r) > 1.0) &
 347                            .and. (iv % radar(n) % rf(k) % qc >= obs_qc_pointer)
 348
 349             if (use_radar_rv) then
 350                !set qc to missing_data for rv of -999.99 (radar_non_precip_rf)
 351                if ( abs(ob%radar(n)%rv(k) - radar_non_precip_rf) < 0.1 ) then
 352                    iv % radar(n) % rv(k) % qc = missing_data
 353                end if
 354                if (abs(iv % radar(n) % rv(k) % qc - missing_data) > 1) then
 355                   if (abs(ob % radar(n) % rv(k) - missing_r) > 1.0 .AND. &
 356                        iv % radar(n) % rv(k) % qc >= obs_qc_pointer) then
 357
 358                      !reference: Sun and Crook (1997)
 359                      call da_radial_velocity(model_rv(k,n), model_p(k,n),  &
 360                         model_u(k,n), model_v(k,n), model_w(k,n),          &
 361                         model_qrn(k,n), model_ps(n), xr, yr, zr,           &
 362                         model_rho(k,n))
 363
 364                      iv % radar(n) % rv(k) % inv = ob % radar(n) % rv(k) - model_rv(k,n)*radar_rv_rscl
 365                   end if
 366                end if
 367             end if
 368
 369             if (use_radar_rf .and. radar_rf_opt==2) then
 370               iv % radar(n) % zmm(k) % inv = 0
 371               iv % radar(n) % rf (k) % qc = -5 ! assume bad rf obs first
 372               echo_rf_good = ob % radar(n) % rf(k) >= rfmin
 373               if ( echo_rf_good ) then
 374                   tmk=model_tc(k,n)+273.15
 375                   prs=model_p(k,n)
 376                   zmm=0
 377                   zmm_ref=0
 378                   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,        &
 379                                        0,0,0,rn0_r,rn0_s,rn0_g,                                                                   &
 380                                        rhos,rhog,dtmk,dqvp,dqra,dqsn,dqgr,dqnr,dqns,dqng,zmm,0,                                   &
 381                                        0,zmm_ref)
 382                   model_rf(k,n)=dbz*radar_rf_rscl
 383                   iv % radar(n) % zmm(k) % inv = zmm
 384                   iv % radar(n) % rf(k) % inv = ob % radar(n) % rf(k) - model_rf(k,n)
 385                   if (model_rf(k,n) >= rfmin) iv % radar(n) % rf (k) % qc  = 0
 386               end if
 387             end if
 388
 389             if (use_radar_rf.and.radar_rf_opt==1) then
 390                if ( echo_rf_good ) then
 391                   model_rf(k,n) = leh1 + leh2 * alog10(model_rho(k,n) * model_qrn(k,n) * 1.0e+3)
 392                   iv % radar(n) % rf(k) % inv = ob % radar(n) % rf(k) - model_rf(k,n)
 393                end if
 394             end if
 395
 396             ! calculate background/model reflectivity
 397             if (use_radar_rhv .or. use_radar_rqv) then
 398                if ( echo_rf_good ) then
 399                   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)
 400                   bg_rf = 10.0*log10(bg_rze)  ! Z to dBZ
 401                   iv % radar(n) % zmm(k) % inv = bg_rf
 402                   iv % radar(n) % rf(k) % inv = ob % radar(n) % rf(k) - bg_rf
 403                end if
 404             end if
 405
 406             ! calculate retrieved hydrometeorological variables
 407             ! Jidong Gao JAS 2013
 408             if (use_radar_rhv) then
 409                if ( echo_rf_good ) then
 410
 411                   ! compute retrieved hydrometeors rhv
 412                   if (it.eq.1) then
 413
 414                      ob_radar_rf = ob % radar(n) % rf(k)
 415
 416                      if ( radar_non_precip_opt > 0 ) then ! assimilate non_precip echo
 417                         if ( echo_non_precip ) then ! ob is non-precip
 418                            if ( bg_rf > -15.0 ) then
 419                               ! when background/model is precip
 420                               ! retrieve hydrometeors from -15dBZ
 421                               ob_radar_rf = -15.0
 422                            else
 423                               ! when background/model is also non-precip
 424                               ! do not need to assimilate this ob
 425                               iv % radar(n) % rf (k) % qc = -5
 426                               iv % radar(n) % rrn(k) % qc = -5
 427                               iv % radar(n) % rsn(k) % qc = -5
 428                               iv % radar(n) % rgr(k) % qc = -5
 429                               cycle level_loop
 430                            end if
 431                         end if !if echo_non_precip
 432                      end if
 433
 434                      ob_radar_rf = min(ob_radar_rf, 55.0) ! if dBZ>55.0, set to 55.0
 435                      rze = 10.0**(ob_radar_rf*0.1) ! dBZ to Z
 436
 437                      if (model_tc(k,n).ge.5.0) then
 438                         ! contribution from rain only
 439                         ! Z_Qr = 3.63*1.0e9*(rho*Qr)**1.75
 440                         iv % radar(n) % rrno(k) = exp ( log(rze/zrr)/1.75 )/model_rho(k,n)
 441                         iv % radar(n) % rrn(k) % qc = 0
 442
 443                         ! rrn and rrno were assigned missing values in read_obs_radar_ascii.inc
 444                         ! maximum value check, use the data under threshold 15g/kg
 445                         iv % radar(n) % rrno(k) = min(iv%radar(n)%rrno(k), 0.015)
 446
 447                      else if (model_tc(k,n).lt.5.0 .and. model_tc(k,n).gt.-5.0 ) then
 448                         ! contribution from rain, snow and graupel
 449                         ! Ze = c * Z_Qr + (1-c) * (Z_Qs+Z_Qg)
 450                         ! the factor c varies linearly between 0 at t=-5C and 1 at t=5C
 451                         czr=(model_tc(k,n)+5)/10.0
 452                         if (model_tc(k,n).le.0.0) then
 453                            czs = (1.0-czr)*zds/(zds+zg) ! dry snow
 454                            czg = (1.0-czr)*zg/(zds+zg)
 455                            iv % radar(n) % rsno(k) = exp ( log(czs*rze/zds)/1.75 )/model_rho(k,n)
 456                         else
 457                            czs = (1.0-czr)*zws/(zws+zg) ! wet snow
 458                            czg = (1.0-czr)*zg/(zws+zg)
 459                            iv % radar(n) % rsno(k) = exp ( log(czs*rze/zws)/1.75 )/model_rho(k,n)
 460                         end if
 461                         iv % radar(n) % rrno(k) = exp ( log(czr*rze/zrr)/1.75 )/model_rho(k,n)
 462                         iv % radar(n) % rgro(k) = exp ( log(czg*rze/zg )/1.75 )/model_rho(k,n)
 463                         iv % radar(n) % rrn(k) % qc = 0
 464                         iv % radar(n) % rsn(k) % qc = 0
 465                         iv % radar(n) % rgr(k) % qc = 0
 466
 467                      else if (model_tc(k,n).le.-5.0) then
 468                         ! contribution from snow and graupel
 469                         czs = zds/(zds+zg)
 470                         czg = 1.0 - czs
 471                         iv % radar(n) % rsno(k) = exp ( log(czs*rze/zds)/1.75 )/model_rho(k,n)
 472                         iv % radar(n) % rgro(k) = exp ( log(czg*rze/zg )/1.75 )/model_rho(k,n)
 473                         iv % radar(n) % rsn(k) % qc = 0
 474                         iv % radar(n) % rgr(k) % qc = 0
 475                      end if  ! temp
 476
 477                      if ( radar_rhv_err_opt == 1 ) then
 478                         ! rainwater error
 479                         iv % radar(n) % rrn(k) % error = iv % radar(n) % rf(k) % error * iv % radar(n) % rrno(k) * alog_10/leh2
 480                         iv % radar(n) % rrn(k) % error = amax1(0.0005,iv % radar(n) % rrn(k) % error)
 481                         iv % radar(n) % rrn(k) % error = amin1( 0.001,iv % radar(n) % rrn(k) % error)
 482                         ! snow error
 483                         iv % radar(n) % rsn(k) % error = iv % radar(n) % rf(k) % error * iv % radar(n) % rsno(k) * alog_10/leh2
 484                         iv % radar(n) % rsn(k) % error = amax1(0.0005,iv % radar(n) % rsn(k) % error)
 485                         iv % radar(n) % rsn(k) % error = amin1( 0.001,iv % radar(n) % rsn(k) % error)
 486                         ! graupel error
 487                         iv % radar(n) % rgr(k) % error = iv % radar(n) % rf(k) % error * iv % radar(n) % rgro(k) * alog_10/leh2
 488                         iv % radar(n) % rgr(k) % error = amax1(0.0005,iv % radar(n) % rgr(k) % error)
 489                         iv % radar(n) % rgr(k) % error = amin1( 0.001,iv % radar(n) % rgr(k) % error)
 490                      else if ( radar_rhv_err_opt == 2 ) then
 491                         ! use settings in the namelist
 492                         ! rainwater error
 493                         iv % radar(n) % rrn(k) % error = radar_rhv_rrn_err * 0.001 ! g/kg to kg/kg
 494                         ! snow error
 495                         iv % radar(n) % rsn(k) % error = radar_rhv_rsn_err * 0.001 ! g/kg to kg/kg
 496                         ! graupel error
 497                         iv % radar(n) % rgr(k) % error = radar_rhv_rgr_err * 0.001 ! g/kg to kg/kg
 498                      else
 499                         write(unit=message(1),fmt='(A)') "radar_rhv_err_opt should be 1 or 2"
 500                         call da_error(__FILE__,__LINE__,message(1:1))
 501                      end if
 502
 503                   end if  ! it=1
 504
 505                   if (iv % radar(n) % rrn(k) % qc >= obs_qc_pointer) then
 506                      iv % radar(n) % rrn(k) % inv = iv % radar(n) % rrno(k) - model_qrn(k,n)
 507                   end if
 508                   if (iv % radar(n) % rsn(k) % qc >= obs_qc_pointer) then
 509                      iv % radar(n) % rsn(k) % inv = iv % radar(n) % rsno(k) - model_qsn(k,n)
 510                   end if
 511                   if (iv % radar(n) % rgr(k) % qc >= obs_qc_pointer) then
 512                      iv % radar(n) % rgr(k) % inv = iv % radar(n) % rgro(k) - model_qgr(k,n)
 513                   end if
 514
 515                end if ! echo_rf_good check
 516             end if ! rhv
 517
 518             ! retrieved water vapor
 519             if (use_radar_rqv) then
 520                !iv%%rqv and iv%%rqvo were assigned to missing values in read_obs_radar_ascii.inc
 521                !iter=1, rqv is missing; for second loop, dont change rqv value
 522                if (it .eq. 1) then
 523                   iv % radar(n) % rqvo(k) = 1.0*model_qs(k,n)
 524                   iv % radar(n) % rqv(k) % error = amax1(0.0005,0.1*iv % radar(n) % rqvo(k))
 525                end if
 526
 527                if ( echo_rf_good ) then
 528
 529                   ! initialize as not-assimilating rqv
 530                   iv % radar(n) % rqv(k) % qc = -5
 531
 532                   if ( echo_non_precip ) then ! ob is non-precip
 533                      if ( radar_non_precip_opt > 0 ) then ! assimilate non_precip echo
 534
 535                         if ( bg_rf >= 20.0 .and. iv%radar(n)%height(k) > model_lcl(n) ) then
 536                            iv % radar(n) % rqv(k) % qc = 0
 537
 538                            if ( model_tc(k,n) > 5.0 ) then
 539                               iv%radar(n)%rqvo(k) =  0.01*radar_non_precip_rh_w*model_qs(k,n)
 540                            else if ( model_tc(k,n) < -5.0 ) then
 541                               iv%radar(n)%rqvo(k) =  0.01*radar_non_precip_rh_i*model_qs_ice(k,n)
 542                            else
 543                               cwr = (model_tc(k,n)+5.0)/10.0
 544                               cws = 1.0 - cwr
 545                               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)
 546                            end if
 547
 548                            iv % radar(n) % rqv(k) % inv = iv % radar(n) % rqvo(k) - model_qv(k,n)
 549                            ! iv % radar(n) % rqv(k) % inv must be <= 0.0
 550                            iv % radar(n) % rqv(k) % inv = amin1(0.0,iv % radar(n) % rqv(k) % inv )
 551                            iv % radar(n) % rqv(k) % error = amax1(0.0005,0.10*model_qs(k,n))
 552                         end if
 553                      end if
 554
 555                   else  ! ob is precip
 556                      if ( ob % radar(n) % rf(k) >= radar_saturated_rf  .and. iv%radar(n)%height(k) >= model_lcl(n) )then
 557                         iv % radar(n) % rqv(k) % qc = 0
 558
 559                         if ( radar_rqv_h_lbound >= 0.0 .and. radar_rqv_h_ubound >= 0.0 ) then
 560                            ! height ranges for applying retrieved rqv are set in the namelist
 561                            if ( iv%radar(n)%height(k) < radar_rqv_h_lbound .or. &
 562                                 iv%radar(n)%height(k) > radar_rqv_h_ubound ) then
 563                               ! do not assimilate rqv outside the specified height ranges
 564                               iv % radar(n) % rqv(k) % qc = -5
 565                            end if
 566                         end if
 567                      end if
 568
 569                      if (iv % radar(n) % rqv(k) % qc >= obs_qc_pointer) then
 570                         ! reduce rqvo for certain rf ranges
 571                         if ( radar_rqv_thresh1 >= 0.0 ) then
 572                            if ( ob % radar(n) % rf(k) >= radar_saturated_rf .and. &
 573                                 ob % radar(n) % rf(k) < radar_rqv_thresh1 ) then
 574                               if ( radar_rqv_rh1 >= 0.0 ) then ! RH in percentage
 575                                  iv % radar(n) % rqvo(k)= 0.01*radar_rqv_rh1*model_qs(k,n)
 576                               end if
 577                            end if
 578                         end if
 579                         if ( radar_rqv_thresh1 >= 0.0 .and. radar_rqv_thresh2 >= 0.0 .and. &
 580                              radar_rqv_thresh2 >= radar_rqv_thresh1 ) then
 581                            if ( ob % radar(n) % rf(k) >= radar_rqv_thresh1 .and. &
 582                                 ob % radar(n) % rf(k) < radar_rqv_thresh2) then
 583                               if ( radar_rqv_rh2 >= 0.0 ) then ! RH in percentage
 584                                  iv % radar(n) % rqvo(k)= 0.01*radar_rqv_rh2*model_qs(k,n)
 585                               end if
 586                            end if
 587                         end if
 588                         iv % radar(n) % rqv(k) % inv = iv % radar(n) % rqvo(k) - model_qv(k,n)
 589                         ! iv % radar(n) % rqv(k) % inv must be >= 0.0
 590                         iv % radar(n) % rqv(k) % inv = amax1(0.0,iv % radar(n) % rqv(k) % inv )
 591                         iv % radar(n) % rqv(k) % error = amax1(0.001,0.20*iv % radar(n) % rqvo(k))
 592                      end if
 593
 594                   end if ! echo_non_precip check
 595
 596                end if  ! echo_rf_good check
 597             end if  ! use_radar_rqv
 598
 599          end if  ! not surface or model lid
 600       end do level_loop
 601    end do
 602
 603    !------------------------------------------------------------------------
 604    ! [4.0] Perform optional maximum error check:
 605    !------------------------------------------------------------------------
 606
 607    if (check_max_iv)  then
 608       call da_check_max_iv_radar(iv, it, irv, irf, irvf, irff)
 609    end if
 610
 611    if (rootproc .and. check_max_iv_print) then
 612       write(unit = check_max_iv_unit, fmt ='(/,A,i5,A)')   &
 613          'For outer iteration ', it, ', Total Rejections for radar follows:'
 614
 615       if (use_radar_rv) then
 616           write( unit = check_max_iv_unit, fmt = '(/,2(A,I6))') &
 617             'Number of failed rv observations:     ',irvf, ' on ',irv
 618       end if
 619
 620       if (use_radar_rf) then
 621          write( unit = check_max_iv_unit, fmt = '(/,2(A,I6))') &
 622             'Number of failed rf observations:     ',irff, ' on ',irf
 623       end if
 624    end if
 625
 626    deallocate (model_p)
 627    deallocate (model_u)
 628    deallocate (model_v)
 629    deallocate (model_w)
 630
 631    if ( allocated(model_rv) ) deallocate (model_rv)
 632    if ( allocated(model_rf) ) deallocate (model_rf)
 633    deallocate (model_ps)
 634
 635    deallocate (model_qv)
 636    deallocate (model_qs)
 637    deallocate (model_tc)
 638
 639    deallocate (model_qrn)
 640    deallocate (model_rho)
 641    deallocate (model_qcl)
 642    deallocate (model_qci)
 643    deallocate (model_qsn)
 644    deallocate (model_qgr)
 645
 646    if ( use_radar_rqv ) then
 647       deallocate (model_lcl)
 648       deallocate (model_qs_ice)
 649    end if
 650
 651    if (trace_use) call da_trace_exit("da_get_innov_vector_radar")
 652
 653 end subroutine da_get_innov_vector_radar
 654
 655