Merge remote-tracking branch 'origin/release-v4.6.1'
[WRF.git] / var / da / da_tamdar / da_get_innov_vector_tamdar_sfc.inc
blob2b1ecee9a5f75f6b4489bf48d36b1b4d354ab8da
1 subroutine da_get_innov_vector_tamdar_sfc( it, num_qcstat_conv, grid, ob, iv)
3    !-----------------------------------------------------------------------
4    ! Purpose: TBD    
5    !-----------------------------------------------------------------------
7    implicit none
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(:,:,:,:)
15    integer :: n        ! Loop counter.
16    integer :: i, j, k  ! Index dimension.
17    real    :: dx, dxm  ! Interpolation weights.
18    real    :: dy, dym  ! Interpolation weights.
19    real, allocatable :: model_u(:,:)  ! Model value u at oblocation.
20    real, allocatable :: model_v(:,:)  ! Model value v at oblocation.
21    real, allocatable :: model_t(:,:)  ! Model value t at oblocation.
22    real, allocatable :: model_p(:,:)  ! Model value p at oblocation.
23    real, allocatable :: model_q(:,:)  ! Model value q at oblocation.
24    real, allocatable :: model_hsm(:,:)
26    real  :: speed, direction
28    real    :: v_h(kms:kme)      ! Model value h at ob hor. location.
29    real    :: v_p(kms:kme)      ! Model value p at ob hor. location.
31    real    :: hd, psfcm
32    real    :: ho, to, qo
33    
34    if (trace_use_dull) call da_trace_entry("da_get_innov_vector_tamdar_sfc")
37    allocate (model_u(1,iv%info(tamdar_sfc)%n1:iv%info(tamdar_sfc)%n2))
38    allocate (model_v(1,iv%info(tamdar_sfc)%n1:iv%info(tamdar_sfc)%n2))
39    allocate (model_t(1,iv%info(tamdar_sfc)%n1:iv%info(tamdar_sfc)%n2))
40    allocate (model_p(1,iv%info(tamdar_sfc)%n1:iv%info(tamdar_sfc)%n2))
41    allocate (model_q(1,iv%info(tamdar_sfc)%n1:iv%info(tamdar_sfc)%n2))
42    allocate (model_hsm(1,iv%info(tamdar_sfc)%n1:iv%info(tamdar_sfc)%n2))
44    if ( it > 1 ) then
45       do n=iv%info(tamdar_sfc)%n1,iv%info(tamdar_sfc)%n2
46          if (iv%tamdar_sfc(n)%u%qc == fails_error_max) iv%tamdar_sfc(n)%u%qc = 0
47          if (iv%tamdar_sfc(n)%v%qc == fails_error_max) iv%tamdar_sfc(n)%v%qc = 0
48          if (iv%tamdar_sfc(n)%t%qc == fails_error_max) iv%tamdar_sfc(n)%t%qc = 0
49          if (iv%tamdar_sfc(n)%p%qc == fails_error_max) iv%tamdar_sfc(n)%p%qc = 0
50          if (iv%tamdar_sfc(n)%q%qc == fails_error_max) iv%tamdar_sfc(n)%q%qc = 0
51       end do
52    end if
54    if (sfc_assi_options == sfc_assi_options_1) then
55       do n=iv%info(tamdar_sfc)%n1,iv%info(tamdar_sfc)%n2
56          ! [1.1] Get horizontal interpolation weights:
58          i   = iv%info(tamdar_sfc)%i(1,n)
59          j   = iv%info(tamdar_sfc)%j(1,n)
60          dx  = iv%info(tamdar_sfc)%dx(1,n)
61          dy  = iv%info(tamdar_sfc)%dy(1,n)
62          dxm = iv%info(tamdar_sfc)%dxm(1,n)
63          dym = iv%info(tamdar_sfc)%dym(1,n)
65          ! Surface correction
67          iv%tamdar_sfc(n)%p%inv = ob%tamdar_sfc(n)%p
68          iv%tamdar_sfc(n)%t%inv = ob%tamdar_sfc(n)%t
69          iv%tamdar_sfc(n)%q%inv = ob%tamdar_sfc(n)%q
70          iv%tamdar_sfc(n)%u%inv = ob%tamdar_sfc(n)%u
71          iv%tamdar_sfc(n)%v%inv = ob%tamdar_sfc(n)%v
73          if (iv % tamdar_sfc(n) % h > missing_r) then
74             do k=kts,kte
75                v_h(k) = dym*(dxm*grid%xb%h(i,j  ,k) + dx*grid%xb%h(i+1,j  ,k)) &
76                   + dy *(dxm*grid%xb%h(i,j+1,k) + dx*grid%xb%h(i+1,j+1,k))
77             end do
79             hd = v_h(kts) - iv % tamdar_sfc(n) % h
81             if (abs(hd) <= Max_StHeight_Diff .or. anal_type_verify ) then
82                if (iv % tamdar_sfc(n) % h < v_h(kts)) then
83                   iv%info(tamdar_sfc)%zk(:,n) = 1.0+1.0e-6
84                   call da_obs_sfc_correction(iv%info(tamdar_sfc), iv%tamdar_sfc(n), n, grid%xb)
85                else
86                   call da_to_zk(iv % tamdar_sfc(n) % h, v_h, v_interp_h, iv%info(tamdar_sfc)%zk(1,n))
87                end if
88             end if
89          else if (ob % tamdar_sfc(n) % p > 1.0) then
90             do k=kts,kte
91                v_p(k) = dym*(dxm*grid%xb%p(i,j  ,k) + dx*grid%xb%p(i+1,j  ,k)) &
92                       + dy *(dxm*grid%xb%p(i,j+1,k) + dx*grid%xb%p(i+1,j+1,k))
93             end do
95             call da_to_zk(ob % tamdar_sfc(n) % p, v_p, v_interp_p, iv%info(tamdar_sfc)%zk(1,n))
97             if (iv%info(tamdar_sfc)%zk(1,n) < 0.0 .and. .not. anal_type_verify) then
98                iv % tamdar_sfc(n) % p % inv = missing_r
99                iv % tamdar_sfc(n) % p % qc  = missing_data
100                iv%info(tamdar_sfc)%zk(:,n) = 1.0+1.0e-6
101             end if
102          end if
103       end do
105       call da_convert_zk (iv%info(tamdar_sfc))
107       if (.not.anal_type_verify ) then
108          do n=iv%info(tamdar_sfc)%n1,iv%info(tamdar_sfc)%n2
109             if (iv%info(tamdar_sfc)%zk(1,n) < 0.0) then
110                iv % tamdar_sfc(n) % u % qc = missing_data
111                iv % tamdar_sfc(n) % v % qc = missing_data
112                iv % tamdar_sfc(n) % t % qc = missing_data
113                iv % tamdar_sfc(n) % q % qc = missing_data
114                iv % tamdar_sfc(n) % p % qc = missing_data
115             end if
116          end do
117       end if
119       ! [1.2] Interpolate horizontally:
120 #ifdef A2C
121       call da_interp_lin_3d (grid%xb%u, iv%info(tamdar_sfc), model_u,'u')
122       call da_interp_lin_3d (grid%xb%v, iv%info(tamdar_sfc), model_v,'v')
123 #else
124       call da_interp_lin_3d (grid%xb%u, iv%info(tamdar_sfc), model_u)
125       call da_interp_lin_3d (grid%xb%v, iv%info(tamdar_sfc), model_v)
126 #endif
127       call da_interp_lin_3d (grid%xb%t, iv%info(tamdar_sfc), model_t)
128       call da_interp_lin_3d (grid%xb%q, iv%info(tamdar_sfc), model_q)
129       call da_interp_lin_3d (grid%xb%p, iv%info(tamdar_sfc), model_p)
131    else if (sfc_assi_options == sfc_assi_options_2) then
132       ! 1.2.1 Surface assimilation approach 2(10-m u, v, 2-m t, q, and sfc_p)
133 #ifdef A2C
134       call da_interp_lin_3d (grid%xb%u, iv%info(tamdar_sfc), model_u,'u')
135       call da_interp_lin_3d (grid%xb%v, iv%info(tamdar_sfc), model_v,'v')
136 #else
137       call da_interp_lin_3d (grid%xb%u, iv%info(tamdar_sfc), model_u)
138       call da_interp_lin_3d (grid%xb%v, iv%info(tamdar_sfc), model_v)
139 #endif
140       call da_interp_lin_2d (grid%xb%t2,   iv%info(tamdar_sfc), 1,model_t)
141       call da_interp_lin_2d (grid%xb%q2,   iv%info(tamdar_sfc), 1,model_q)
142       call da_interp_lin_2d (grid%xb%psfc, iv%info(tamdar_sfc), 1,model_p)
144       do n=iv%info(tamdar_sfc)%n1,iv%info(tamdar_sfc)%n2
146          iv%tamdar_sfc(n)%p%inv = ob%tamdar_sfc(n)%p
147          iv%tamdar_sfc(n)%t%inv = ob%tamdar_sfc(n)%t
148          iv%tamdar_sfc(n)%q%inv = ob%tamdar_sfc(n)%q
149          iv%tamdar_sfc(n)%u%inv = ob%tamdar_sfc(n)%u
150          iv%tamdar_sfc(n)%v%inv = ob%tamdar_sfc(n)%v
151          if (iv%tamdar_sfc(n)%p%qc >= 0) then
152             ! model surface p, t, q, h at observed site:
154             call da_interp_lin_2d_partial (grid%xb%terr, iv%info(tamdar_sfc), 1, n, n, model_hsm(:,n))
156             ho = iv%tamdar_sfc(n)%h
157             to = -888888.0
158             qo = -888888.0
160             if (iv%tamdar_sfc(n)%t%qc >= 0 .and. iv%tamdar_sfc(n)%q%qc >= 0) then
161                to = ob%tamdar_sfc(n)%t
162                qo = ob%tamdar_sfc(n)%q
163                call da_sfc_pre(psfcm, model_p(1,n), model_t(1,n), model_q(1,n), model_hsm(1,n), ho, to, qo)
164             else if (iv%tamdar_sfc(n)%t%qc >= 0 .and. iv%tamdar_sfc(n)%q%qc < 0) then
165                to = ob%tamdar_sfc(n)%t
166                call da_sfc_pre(psfcm, model_p(1,n), model_t(1,n), model_q(1,n), model_hsm(1,n), ho,to)
167             else
168                call da_sfc_pre(psfcm, model_p(1,n), model_t(1,n), model_q(1,n), model_hsm(1,n), ho)
169             end if
171             ! Pressure at the observed height:
172             model_p(1,n) = psfcm
173          end if
174       end do
175    end if
177    do n=iv%info(tamdar_sfc)%n1,iv%info(tamdar_sfc)%n2
178       !-----------------------------------------------------------------------
179       ! [3.0] Fast interpolation:
180       !-----------------------------------------------------------------------
182        if(wind_sd_tamdar)then
183          call da_ffdduv_model (speed,direction,model_u(1,n), model_v(1,n), convert_uv2fd)
185          if (ob%tamdar_sfc(n)%u > missing_r .AND. iv%tamdar_sfc(n)%u%qc >= obs_qc_pointer) then
186              iv%tamdar_sfc(n)%u%inv = iv%tamdar_sfc(n)%u%inv - speed
187          else
188              iv % tamdar_sfc(n) % u % inv = 0.0
189          end if
191          if (ob%tamdar_sfc(n)%v > missing_r .AND. iv%tamdar_sfc(n)%v%qc >= obs_qc_pointer) then
192              iv%tamdar_sfc(n)%v%inv = iv%tamdar_sfc(n)%v%inv - direction
193              if (iv%tamdar_sfc(n)%v%inv > 180.0 ) iv%tamdar_sfc(n)%v%inv = iv%tamdar_sfc(n)%v%inv - 360.0
194              if (iv%tamdar_sfc(n)%v%inv < -180.0 ) iv%tamdar_sfc(n)%v%inv = iv%tamdar_sfc(n)%v%inv + 360.0
195          else
196              iv % tamdar_sfc(n) % v % inv = 0.0
197          end if
198        else
199           if (ob % tamdar_sfc(n) % u > missing_r .AND. iv % tamdar_sfc(n) % u % qc >= obs_qc_pointer) then
200               iv % tamdar_sfc(n) % u % inv = iv%tamdar_sfc(n)%u%inv - model_u(1,n)
201           else
202               iv % tamdar_sfc(n) % u % inv = 0.0
203           end if
205           if (ob % tamdar_sfc(n) % v > missing_r .AND. iv % tamdar_sfc(n) % v % qc >= obs_qc_pointer) then
206               iv % tamdar_sfc(n) % v % inv = iv%tamdar_sfc(n)%v%inv - model_v(1,n)
207           else
208               iv % tamdar_sfc(n) % v % inv = 0.0
209           end if
210        end if
213       if (ob % tamdar_sfc(n) % p > 0.0 .AND. iv % tamdar_sfc(n) % p % qc >= obs_qc_pointer) then
214          iv % tamdar_sfc(n) % p % inv = iv%tamdar_sfc(n)%p%inv - model_p(1,n)
215       else
216          iv % tamdar_sfc(n) % p % inv = 0.0
217       end if
219       if (ob % tamdar_sfc(n) % t > 0.0 .AND. iv % tamdar_sfc(n) % t % qc >= obs_qc_pointer) then
220          iv % tamdar_sfc(n) % t % inv = iv%tamdar_sfc(n)%t%inv - model_t(1,n)
221       else
222          iv % tamdar_sfc(n) % t % inv = 0.0
223       end if
225       if (ob % tamdar_sfc(n) % q > 0.0 .AND. iv % tamdar_sfc(n) % q % qc >= obs_qc_pointer) then
226          iv % tamdar_sfc(n) % q % inv = iv%tamdar_sfc(n)%q%inv - model_q(1,n)
227       else
228          iv % tamdar_sfc(n) % q % inv = 0.0
229       end if
230    end do
232    !-----------------------------------------------------------------------
233    !     [5.0] Perform optional maximum error check:
234    !-----------------------------------------------------------------------
236    if ( check_max_iv ) &
237       call da_check_max_iv_tamdar_sfc(iv,ob, it, num_qcstat_conv) 
239    deallocate (model_u)
240    deallocate (model_v)
241    deallocate (model_t)
242    deallocate (model_p)
243    deallocate (model_q)
244    deallocate (model_hsm)
245    if (trace_use_dull) call da_trace_exit("da_get_innov_vector_tamdar_sfc")
247 end subroutine da_get_innov_vector_tamdar_sfc