var/da/da_radar/da_transform_xtoy_radar_adj.inc

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