var/da/da_radar/da_transform_xtoy_radar.inc

   1 subroutine da_transform_xtoy_radar (grid, iv, y)
   2
   3    !-----------------------------------------------------------------------
   4    ! Purpose: calculate the Doppler radial velocity and
   5    ! reflectivity at the observation location from the first guess.
   6    ! It is linearized.
   7    ! History:
   8    !    Updated for Analysis on Arakawa-C grid
   9    !    Author: Syed RH Rizvi,  MMM/ESSL/NCAR,  Date: 10/22/2008
  10    !    08/2017 - bug fix for Vr operator (Siou-Ying Jiang, CWB, Taiwan)
  11    !---------------------------------------------------------------------
  12
  13    implicit none
  14
  15    type (domain),  intent(in)    :: grid
  16    type (iv_type), intent(in)    :: iv       ! Innovation vector (O-B).
  17    type (y_type),  intent(inout) :: y        ! y = h (grid%xa) (linear)
  18
  19    integer :: n, k
  20
  21    real, allocatable :: model_p(:,:)
  22    real, allocatable :: model_rho(:,:)
  23    real, allocatable :: model_u(:,:)
  24    real, allocatable :: model_v(:,:)
  25    real, allocatable :: model_w(:,:)
  26    real, allocatable :: model_qrn(:,:)
  27    real, allocatable :: model_qrnb(:,:)
  28    real, allocatable :: model_ps(:)
  29    real, allocatable :: model_qsn(:,:)
  30    real, allocatable :: model_qgr(:,:)
  31    real, allocatable :: model_qv(:,:)
  32    real, allocatable :: model_qvb(:,:)
  33    real, allocatable :: model_t(:,:)
  34    real, allocatable :: model_tb(:,:)
  35    !------------------------
  36    !  for jung et al 2008
  37    !------------------------
  38    real, allocatable :: model_qsnb(:,:)
  39    real, allocatable :: model_qgrb(:,:)
  40    real, allocatable :: model_qnrb(:,:)
  41    real, allocatable :: model_qnsb(:,:)
  42    real, allocatable :: model_qngb(:,:)
  43    real, allocatable :: model_qnr(:,:)
  44    real, allocatable :: model_qns(:,:)
  45    real, allocatable :: model_qng(:,:)
  46    !------------------------
  47
  48    real    :: xr,yr,zr
  49
  50    real    :: alog_10
  51
  52    !------------------------
  53    !  for jung et al 2008
  54    !------------------------
  55    real    :: qvp,qra,qsn,qgr ! mixing ratio
  56                               ! to match the type of argument in subroutine
  57    real    :: dqra,dqsn,dqgr,dtmk,dqvp
  58    real    :: dqnr,dqns,dqng
  59    real    :: zmm,zmm_ref
  60    real    :: qnr,qns,qng            ! number concentration of rain snow and graupel
  61    real    :: tmk,prs                ! temperature and pressure
  62    real    :: dbz                    ! reflectivity in dBZ
  63    real    :: rn0_r,rn0_s,rn0_g      ! intercept parameter of rain snow and graupel
  64    real    :: rhos,rhog              ! density of snow and graupel
  65    !------------------------
  66
  67    if (trace_use) call da_trace_entry("da_transform_xtoy_radar")
  68
  69    alog_10 = alog(10.0)
  70
  71    !------------------------
  72    !  for jung et al 2008
  73    !------------------------
  74    qnr=0
  75    qns=0
  76    qng=0
  77    rn0_r=8e6
  78    rn0_s=3e6
  79    rn0_g=4e6
  80    rhos=100.0
  81    rhog=400.0
  82    !------------------------
  83
  84    allocate (model_p(iv%info(radar)%max_lev,iv%info(radar)%n1:iv%info(radar)%n2))
  85    allocate (model_rho(iv%info(radar)%max_lev,iv%info(radar)%n1:iv%info(radar)%n2))
  86    allocate (model_u(iv%info(radar)%max_lev,iv%info(radar)%n1:iv%info(radar)%n2))
  87    allocate (model_v(iv%info(radar)%max_lev,iv%info(radar)%n1:iv%info(radar)%n2))
  88    allocate (model_w(iv%info(radar)%max_lev,iv%info(radar)%n1:iv%info(radar)%n2))
  89    allocate (model_qrn(iv%info(radar)%max_lev,iv%info(radar)%n1:iv%info(radar)%n2))
  90    allocate (model_qrnb(iv%info(radar)%max_lev,iv%info(radar)%n1:iv%info(radar)%n2))
  91    allocate (model_ps(iv%info(radar)%n1:iv%info(radar)%n2))
  92    allocate (model_qsn(iv%info(radar)%max_lev,iv%info(radar)%n1:iv%info(radar)%n2))
  93    allocate (model_qgr(iv%info(radar)%max_lev,iv%info(radar)%n1:iv%info(radar)%n2))
  94    allocate (model_qv(iv%info(radar)%max_lev,iv%info(radar)%n1:iv%info(radar)%n2))
  95    allocate (model_t(iv%info(radar)%max_lev,iv%info(radar)%n1:iv%info(radar)%n2))
  96
  97    !------------------------
  98    !  for jung et al 2008
  99    !------------------------
 100    if (use_radar_rf .and. radar_rf_opt==2) then
 101      allocate (model_qsnb(iv%info(radar)%max_lev,iv%info(radar)%n1:iv%info(radar)%n2))
 102      allocate (model_qgrb(iv%info(radar)%max_lev,iv%info(radar)%n1:iv%info(radar)%n2))
 103      allocate (model_qnrb(iv%info(radar)%max_lev,iv%info(radar)%n1:iv%info(radar)%n2))
 104      allocate (model_qnsb(iv%info(radar)%max_lev,iv%info(radar)%n1:iv%info(radar)%n2))
 105      allocate (model_qngb(iv%info(radar)%max_lev,iv%info(radar)%n1:iv%info(radar)%n2))
 106      allocate (model_qnr(iv%info(radar)%max_lev,iv%info(radar)%n1:iv%info(radar)%n2))
 107      allocate (model_qns(iv%info(radar)%max_lev,iv%info(radar)%n1:iv%info(radar)%n2))
 108      allocate (model_qng(iv%info(radar)%max_lev,iv%info(radar)%n1:iv%info(radar)%n2))
 109    end if
 110    !------------------------
 111
 112    do n=iv%info(radar)%n1,iv%info(radar)%n2
 113       do k = 1, iv%info(radar)%levels(n)
 114          model_qrnb(k,n) = iv%radar(n)%model_qrn(k)
 115          model_p(k,n)    = iv%radar(n)%model_p(k)
 116          model_rho(k,n)  = iv%radar(n)%model_rho(k)
 117          !------------------------
 118          !  for jung et al 2008
 119          !------------------------
 120         if (use_radar_rf .and. radar_rf_opt==2) then
 121          model_qsnb(k,n) = iv%radar(n)%model_qsn(k)
 122          model_qgrb(k,n) = iv%radar(n)%model_qgr(k)
 123          model_qnrb(k,n) = 0
 124          model_qnsb(k,n) = 0
 125          model_qngb(k,n) = 0
 126         end if
 127          !------------------------
 128       end do
 129
 130       model_ps(n) = iv%radar(n)%model_ps
 131    end do
 132
 133    ! [1.4] Interpolate horizontally from dot points:
 134 #ifdef A2C
 135    call da_interp_lin_3d (grid%xa%u,   iv%info(radar), model_u,'u')
 136    call da_interp_lin_3d (grid%xa%v,   iv%info(radar), model_v,'v')
 137 #else
 138    call da_interp_lin_3d (grid%xa%u,   iv%info(radar), model_u)
 139    call da_interp_lin_3d (grid%xa%v,   iv%info(radar), model_v)
 140 #endif
 141    call da_interp_lin_3d (grid%xa%wh,  iv%info(radar), model_w)
 142    model_qsn = 0.0
 143    model_qgr = 0.0
 144    model_qrn = 0.0
 145    if ( cloud_cv_options >= 1 ) then
 146       call da_interp_lin_3d (grid%xa%qrn, iv%info(radar), model_qrn)
 147       if ( cloud_cv_options >= 2 ) then
 148          call da_interp_lin_3d (grid%xa%qsn, iv%info(radar), model_qsn)
 149          call da_interp_lin_3d (grid%xa%qgr, iv%info(radar), model_qgr)
 150       end if
 151    end if
 152    call da_interp_lin_3d (grid%xa%q,   iv%info(radar), model_qv)
 153    call da_interp_lin_3d (grid%xa%t,   iv%info(radar), model_t)
 154
 155    !if ( use_radar_rqv ) then
 156       !basic states
 157       allocate (model_tb(iv%info(radar)%max_lev,iv%info(radar)%n1:iv%info(radar)%n2))
 158       allocate (model_qvb(iv%info(radar)%max_lev,iv%info(radar)%n1:iv%info(radar)%n2))
 159       call da_interp_lin_3d (grid%xb%t,   iv%info(radar), model_tb)
 160       call da_interp_lin_3d (grid%xb%q,   iv%info(radar), model_qvb)
 161    !end if
 162
 163    do n=iv%info(radar)%n1,iv%info(radar)%n2
 164
 165       ! [1.7] Calculate rv and rf at OBS location
 166
 167       xr = grid%xb%ds * (iv%info(radar)%x(1,n) - iv%radar(n)%stn_loc%x)
 168       yr = grid%xb%ds * (iv%info(radar)%y(1,n) - iv%radar(n)%stn_loc%y)
 169
 170       do k = 1, iv%info(radar)%levels(n)
 171          if (iv % radar(n) % height_qc(k) /= below_model_surface .and.  &
 172               iv % radar(n) % height_qc(k) /= above_model_lid) then
 173             if (use_radar_rv) then
 174                if (iv % radar(n) % rv(k) % qc >= obs_qc_pointer) then
 175                   zr=iv%radar(n)%height(k) - iv%radar(n)%stn_loc%elv
 176
 177                   call da_radial_velocity_lin(y%radar(n)%rv(k), &
 178                      model_p(k,n), &
 179                      model_u(k,n), model_v(k,n), model_w(k,n), model_qrn(k,n),    &
 180                      model_ps(n), xr, yr, zr, model_qrnb(k,n), model_rho(k,n))
 181                   y%radar(n)%rv(k)=y%radar(n)%rv(k)*radar_rv_rscl
 182                end if
 183             end if
 184
 185             if (use_radar_rf .and. radar_rf_opt==2) then
 186               if (iv % radar(n) % rf(k) % qc >= obs_qc_pointer) then
 187                    tmk=model_tb(k,n)
 188                    prs=model_p(k,n)
 189                    zmm=iv % radar(n) % zmm(k) % inv
 190                    zmm_ref=0
 191                    call da_radzicevar_tl(model_qvb(k,n),model_qrnb(k,n),model_qsnb(k,n),model_qgrb(k,n),qnr,qns,qng,tmk,prs,dbz,                   &
 192                                       0,0,0,rn0_r,rn0_s,rn0_g,                                                                                  &
 193                                       rhos,rhog,model_t(k,n),model_qv(k,n),model_qrn(k,n),model_qsn(k,n),model_qgr(k,n),dqnr,dqns,dqng,zmm,2,   &
 194                                       0,zmm_ref)
 195                    y%radar(n)%rf(k) =radar_rf_rscl*dbz
 196               end if
 197             end if
 198
 199             if (use_radar_rf .and. radar_rf_opt==1) then
 200                if (iv % radar(n) % rf(k) % qc >= obs_qc_pointer) then
 201                   y%radar(n)%rf(k) = leh2 * model_qrn(k,n) /(model_qrnb(k,n)*alog_10)
 202                end if
 203             end if
 204
 205             if (.not.use_radar_rf .and. use_radar_rhv) then
 206                if (iv % radar(n) % rrn(k) % qc >= obs_qc_pointer) then
 207                   y%radar(n)%rrn(k) = model_qrn(k,n)
 208                end if
 209                if (iv % radar(n) % rsn(k) % qc >= obs_qc_pointer) then
 210                   y%radar(n)%rsn(k) = model_qsn(k,n)
 211                end if
 212                if (iv % radar(n) % rgr(k) % qc >= obs_qc_pointer) then
 213                   y%radar(n)%rgr(k) = model_qgr(k,n)
 214                end if
 215             end if
 216
 217             if (use_radar_rqv) then
 218                !dqv=qs*drh+(c2*c3*/(T+c3)**2.0-c4/T)*qv**dT
 219                !c2=17.67 is es_beta
 220                !c3=243.5 is es_gamma
 221                !c4=0.622 is a_ew
 222                ! use qc from get_inv.
 223                if (iv % radar(n) % rqv(k) % qc >= obs_qc_pointer) then
 224                   y%radar(n)%rqv(k) = model_qv(k,n)
 225                   ! Wang JAMC
 226                   y%radar(n)%rqv(k) = y%radar(n)%rqv(k) + ( es_beta*es_gamma/(model_tb(k,n)+es_gamma)**2.0 )*model_qvb(k,n)*model_t(k,n)
 227                end if
 228             end if
 229
 230          end if
 231       end do
 232    end do
 233
 234    deallocate (model_p)
 235    deallocate (model_u)
 236    deallocate (model_v)
 237    deallocate (model_w)
 238    deallocate (model_qrn)
 239    deallocate (model_qrnb)
 240    deallocate (model_ps)
 241    deallocate (model_qsn)
 242    deallocate (model_qgr)
 243    deallocate (model_qv)
 244    deallocate (model_t)
 245    deallocate (model_rho)
 246    !if ( use_radar_rqv ) then
 247       deallocate (model_tb)
 248       deallocate (model_qvb)
 249    !end if
 250
 251    if (use_radar_rf .and. radar_rf_opt==2) then
 252      deallocate (model_qsnb)
 253      deallocate (model_qgrb)
 254      deallocate (model_qnr)
 255      deallocate (model_qns)
 256      deallocate (model_qng)
 257      deallocate (model_qnrb)
 258      deallocate (model_qnsb)
 259      deallocate (model_qngb)
 260    end if
 261
 262    if (trace_use) call da_trace_exit("da_transform_xtoy_radar")
 263
 264 end subroutine da_transform_xtoy_radar