wrfv2_fire/var/da/da_tamdar/da_get_innov_vector_tamdar.inc

   1 subroutine da_get_innov_vector_tamdar (it, num_qcstat_conv, grid, ob, iv)
   2
   3    !-----------------------------------------------------------------------
   4    ! Purpose: TBD
   5    !-----------------------------------------------------------------------
   6
   7    implicit none
   8
   9    integer,          intent(in)    :: it       ! External iteration.
  10    type(domain),     intent(in)    :: grid     ! first guess state.
  11    type(y_type),     intent(inout) :: ob       ! Observation structure.
  12    type(iv_type),    intent(inout) :: iv       ! O-B structure.
  13    integer,          intent(inout) :: num_qcstat_conv(:,:,:,:)
  14
  15    integer :: n, k        ! Loop counter.
  16    integer :: i  (kms:kme)
  17    integer :: j  (kms:kme)
  18    real    :: dx (kms:kme)
  19    real    :: dxm(kms:kme)
  20    real    :: dy (kms:kme)
  21    real    :: dym(kms:kme)
  22
  23    real, allocatable :: model_u(:,:)  ! Model value u at ob location.
  24    real, allocatable :: model_v(:,:)  ! Model value v at ob location.
  25    real, allocatable :: model_t(:,:)  ! Model value t at ob location.
  26    real, allocatable :: model_q(:,:)  ! Model value q at ob location.
  27
  28    real    :: v_h(kts:kte)      ! Model value h at ob hor. location.
  29    real    :: v_p(kts:kte)      ! Model value p at ob hor. location.
  30
  31    if (trace_use_dull) call da_trace_entry ("da_get_innov_vector_tamdar")
  32
  33    allocate (model_u(iv%info(tamdar)%max_lev,iv%info(tamdar)%n1:iv%info(tamdar)%n2))
  34    allocate (model_v(iv%info(tamdar)%max_lev,iv%info(tamdar)%n1:iv%info(tamdar)%n2))
  35    allocate (model_t(iv%info(tamdar)%max_lev,iv%info(tamdar)%n1:iv%info(tamdar)%n2))
  36    allocate (model_q(iv%info(tamdar)%max_lev,iv%info(tamdar)%n1:iv%info(tamdar)%n2))
  37
  38    model_u(:,:) = 0.0
  39    model_v(:,:) = 0.0
  40    model_t(:,:) = 0.0
  41    model_q(:,:) = 0.0
  42
  43    if ( it > 1 ) then
  44       do n=iv%info(tamdar)%n1,iv%info(tamdar)%n2
  45          do k=1, iv%info(tamdar)%levels(n)
  46             if (iv%tamdar(n)%u(k)%qc == fails_error_max) iv%tamdar(n)%u(k)%qc = 0
  47             if (iv%tamdar(n)%v(k)%qc == fails_error_max) iv%tamdar(n)%v(k)%qc = 0
  48             if (iv%tamdar(n)%t(k)%qc == fails_error_max) iv%tamdar(n)%t(k)%qc = 0
  49             if (iv%tamdar(n)%q(k)%qc == fails_error_max) iv%tamdar(n)%q(k)%qc = 0
  50          end do
  51       end do
  52    end if
  53
  54    do n=iv%info(tamdar)%n1, iv%info(tamdar)%n2
  55       if (iv%info(tamdar)%levels(n) < 1) cycle
  56
  57       ! [1.1] Get horizontal interpolation weights:
  58
  59       if (position_lev_dependant) then
  60          i(:)   = iv%info(tamdar)%i(:,n)
  61          j(:)   = iv%info(tamdar)%j(:,n)
  62          dx(:)  = iv%info(tamdar)%dx(:,n)
  63          dy(:)  = iv%info(tamdar)%dy(:,n)
  64          dxm(:) = iv%info(tamdar)%dxm(:,n)
  65          dym(:) = iv%info(tamdar)%dym(:,n)
  66          do k=kts,kte
  67             v_h(k) = dym(k)*(dxm(k)*grid%xb%h(i(k),j(k),k) + dx(k)*grid%xb%h(i(k)+1,j(k),k)) &
  68                + dy(k) *(dxm(k)*grid%xb%h(i(k),j(k)+1,k) + dx(k)*grid%xb%h(i(k)+1,j(k)+1,k))
  69             v_p(k) = dym(k)*(dxm(k)*grid%xb%p(i(k),j(k),k) + dx(k)*grid%xb%p(i(k)+1,j(k),k)) &
  70                + dy(k) *(dxm(k)*grid%xb%p(i(k),j(k)+1,k) + dx(k)*grid%xb%p(i(k)+1,j(k)+1,k))
  71          end do
  72       else
  73          i(1)   = iv%info(tamdar)%i(1,n)
  74          j(1)   = iv%info(tamdar)%j(1,n)
  75          dx(1)  = iv%info(tamdar)%dx(1,n)
  76          dy(1)  = iv%info(tamdar)%dy(1,n)
  77          dxm(1) = iv%info(tamdar)%dxm(1,n)
  78          dym(1) = iv%info(tamdar)%dym(1,n)
  79
  80          v_h(kts:kte) = dym(1) * (dxm(1)*grid%xb%h(i(1),j(1),kts:kte)   + dx(1)*grid%xb%h(i(1)+1,j(1),kts:kte)) &
  81                        + dy(1) * (dxm(1)*grid%xb%h(i(1),j(1)+1,kts:kte) + dx(1)*grid%xb%h(i(1)+1,j(1)+1,kts:kte))
  82          v_p(kts:kte) = dym(1) * (dxm(1)*grid%xb%p(i(1),j(1),kts:kte)   + dx(1)*grid%xb%p(i(1)+1,j(1),kts:kte)) &
  83                        + dy(1) * (dxm(1)*grid%xb%p(i(1),j(1)+1,kts:kte) + dx(1)*grid%xb%p(i(1)+1,j(1)+1,kts:kte))
  84       end if
  85
  86       do k=1, iv%info(tamdar)%levels(n)
  87          if (iv%tamdar(n)%p(k) > 1.0) then
  88             call da_to_zk (iv%tamdar(n)%p(k), v_p, v_interp_p, iv%info(tamdar)%zk(k,n))
  89          else if (iv%tamdar(n)%h(k) > 0.0) then
  90             call da_to_zk (iv%tamdar(n)%h(k), v_h, v_interp_h, iv%info(tamdar)%zk(k,n))
  91          end if
  92       end do
  93
  94    end do
  95
  96    call da_convert_zk (iv%info(tamdar))
  97
  98    if (.not. anal_type_verify) then
  99       do n=iv%info(tamdar)%n1,iv%info(tamdar)%n2
 100          do k=1, iv%info(tamdar)%levels(n)
 101             if (iv%info(tamdar)%zk(k,n) < 0.0) then
 102                iv%tamdar(n)%u(k)%qc = missing_data
 103                iv%tamdar(n)%v(k)%qc = missing_data
 104                iv%tamdar(n)%t(k)%qc = missing_data
 105                iv%tamdar(n)%q(k)%qc = missing_data
 106             end if
 107          end do
 108       end do
 109    end if
 110
 111    ! [1.2] Interpolate horizontally to ob:
 112 #ifdef A2C
 113    call da_interp_lin_3d (grid%xb%u, iv%info(tamdar), model_u,'u')
 114    call da_interp_lin_3d (grid%xb%v, iv%info(tamdar), model_v,'v')
 115 #else
 116    call da_interp_lin_3d (grid%xb%u, iv%info(tamdar), model_u)
 117    call da_interp_lin_3d (grid%xb%v, iv%info(tamdar), model_v)
 118 #endif
 119    call da_interp_lin_3d (grid%xb%t, iv%info(tamdar), model_t)
 120    call da_interp_lin_3d (grid%xb%q, iv%info(tamdar), model_q)
 121
 122    do n=iv%info(tamdar)%n1, iv%info(tamdar)%n2
 123       !----------------------------------------------------------------------
 124       ! [2.0] Initialise components of innovation vector:
 125       !----------------------------------------------------------------------
 126
 127       do k = 1, iv%info(tamdar)%levels(n)
 128          iv%tamdar(n)%u(k)%inv = 0.0
 129          iv%tamdar(n)%v(k)%inv = 0.0
 130          iv%tamdar(n)%t(k)%inv = 0.0
 131          iv%tamdar(n)%q(k)%inv = 0.0
 132
 133          !-------------------------------------------------------------------
 134          ! [3.0] Interpolation:
 135          !-------------------------------------------------------------------
 136
 137          if (ob%tamdar(n)%u(k) > missing_r .AND. iv%tamdar(n)%u(k)%qc >= obs_qc_pointer) then
 138             iv%tamdar(n)%u(k)%inv = ob%tamdar(n)%u(k) - model_u(k,n)
 139          end if
 140
 141          if (ob%tamdar(n)%v(k) > missing_r .AND. iv%tamdar(n)%v(k)%qc >= obs_qc_pointer) then
 142             iv%tamdar(n)%v(k)%inv = ob%tamdar(n)%v(k) - model_v(k,n)
 143          end if
 144
 145          if (ob%tamdar(n)%t(k) > missing_r .AND. iv%tamdar(n)%t(k)%qc >= obs_qc_pointer) then
 146             iv%tamdar(n)%t(k)%inv = ob%tamdar(n)%t(k) - model_t(k,n)
 147          end if
 148
 149          if (ob%tamdar(n)%q(k) > missing_r .AND. iv%tamdar(n)%q(k)%qc >= obs_qc_pointer) then
 150             iv%tamdar(n)%q(k)%inv = ob%tamdar(n)%q(k) - model_q(k,n)
 151          end if
 152       end do
 153    end do
 154
 155    !----------------------------------------------------------------------
 156    ! [5.0] Perform optional maximum error check:
 157    !----------------------------------------------------------------------
 158
 159    if (check_max_iv) call da_check_max_iv_tamdar (iv, it, num_qcstat_conv)
 160
 161    deallocate (model_u)
 162    deallocate (model_v)
 163    deallocate (model_t)
 164    deallocate (model_q)
 165    if (trace_use_dull) call da_trace_exit ("da_get_innov_vector_tamdar")
 166
 167 end subroutine da_get_innov_vector_tamdar
 168
 169