da_get_innov_vector_sound.inc

References to this file elsewhere.
1 subroutine da_get_innov_vector_sound (it, 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 
14    integer :: n, k        ! Loop counter.
15    integer :: i  (kms:kme)
16    integer :: j  (kms:kme)
17    real    :: dx (kms:kme)
18    real    :: dxm(kms:kme)  
19    real    :: dy (kms:kme)
20    real    :: dym(kms:kme)  
21 
22    real, allocatable :: model_u(:,:)  ! Model value u at ob location.
23    real, allocatable :: model_v(:,:)  ! Model value v at ob location.
24    real, allocatable :: model_t(:,:)  ! Model value t at ob location.
25    real, allocatable :: model_q(:,:)  ! Model value q at ob location.
26 
27    real    :: v_h(kts:kte)      ! Model value h at ob hor. location.
28    real    :: v_p(kts:kte)      ! Model value p at ob hor. location.
29 
30    integer :: itu,ituf,itvv,itvvf,itt,ittf,itqv,itqvf
31 
32    if (trace_use_dull) call da_trace_entry ("da_get_innov_vector_sound")
33 
34    itu   = 0; itvv    = 0;  itt  = 0; itqv  = 0;
35    ituf  = 0; itvvf   = 0;  ittf = 0; itqvf = 0;
36 
37    allocate (model_u(iv%info(sound)%max_lev,iv%info(sound)%n1:iv%info(sound)%n2))
38    allocate (model_v(iv%info(sound)%max_lev,iv%info(sound)%n1:iv%info(sound)%n2))
39    allocate (model_t(iv%info(sound)%max_lev,iv%info(sound)%n1:iv%info(sound)%n2))
40    allocate (model_q(iv%info(sound)%max_lev,iv%info(sound)%n1:iv%info(sound)%n2))
41 
42    model_u(:,:) = 0.0
43    model_v(:,:) = 0.0
44    model_t(:,:) = 0.0
45    model_q(:,:) = 0.0
46 
47    do n=iv%info(sound)%n1, iv%info(sound)%n2
48       if (iv%info(sound)%levels(n) < 1) cycle
49 
50       ! [1.1] Get horizontal interpolation weights:
51 
52       if (position_lev_dependant) then
53          i(:)   = iv%info(sound)%i(:,n)
54          j(:)   = iv%info(sound)%j(:,n)
55          dx(:)  = iv%info(sound)%dx(:,n)
56          dy(:)  = iv%info(sound)%dy(:,n)
57          dxm(:) = iv%info(sound)%dxm(:,n)
58          dym(:) = iv%info(sound)%dym(:,n)
59          do k=kts,kte
60             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)) &
61                + 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))
62             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)) &
63                + 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))
64          end do
65       else
66          i(1)   = iv%info(sound)%i(1,n)
67          j(1)   = iv%info(sound)%j(1,n)
68          dx(1)  = iv%info(sound)%dx(1,n)
69          dy(1)  = iv%info(sound)%dy(1,n)
70          dxm(1) = iv%info(sound)%dxm(1,n)
71          dym(1) = iv%info(sound)%dym(1,n)
72 
73          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)) &
74                        + 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))
75          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)) &
76                        + 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))
77       end if
78 
79       do k=1, iv%info(sound)%levels(n)
80          if (iv%sound(n)%p(k) > 1.0) then
81             call da_to_zk (iv%sound(n)%p(k), v_p, v_interp_p, iv%info(sound)%zk(k,n))
82          else if (iv%sound(n)%h(k) > 0.0) then
83             call da_to_zk (iv%sound(n)%h(k), v_h, v_interp_h, iv%info(sound)%zk(k,n))
84          end if
85       end do
86 
87    end do
88 
89    call da_convert_zk (iv%info(sound))
90 
91    if (.not. anal_type_verify) then
92       do n=iv%info(sound)%n1,iv%info(sound)%n2
93          do k=1, iv%info(sound)%levels(n)
94             if (iv%info(sound)%zk(k,n) < 0.0) then
95                iv%sound(n)%u(k)%qc = missing
96                iv%sound(n)%v(k)%qc = missing
97                iv%sound(n)%t(k)%qc = missing
98                iv%sound(n)%q(k)%qc = missing
99             end if
100          end do
101       end do
102    end if
103 
104    ! [1.2] Interpolate horizontally to ob:
105 
106    call da_interp_lin_3d (grid%xb%u, iv%info(sound), model_u)
107    call da_interp_lin_3d (grid%xb%v, iv%info(sound), model_v)
108    call da_interp_lin_3d (grid%xb%t, iv%info(sound), model_t)
109    call da_interp_lin_3d (grid%xb%q, iv%info(sound), model_q)
110 
111    do n=iv%info(sound)%n1, iv%info(sound)%n2
112       !----------------------------------------------------------------------
113       ! [2.0] Initialise components of innovation vector:
114       !----------------------------------------------------------------------
115 
116       do k = 1, iv%info(sound)%levels(n)
117          iv%sound(n)%u(k)%inv = 0.0
118          iv%sound(n)%v(k)%inv = 0.0
119          iv%sound(n)%t(k)%inv = 0.0
120          iv%sound(n)%q(k)%inv = 0.0
121 
122          !-------------------------------------------------------------------
123          ! [3.0] Interpolation:
124          !-------------------------------------------------------------------
125 
126          if (ob%sound(n)%u(k) > missing_r .AND. iv%sound(n)%u(k)%qc >= obs_qc_pointer) then
127             iv%sound(n)%u(k)%inv = ob%sound(n)%u(k) - model_u(k,n)
128          end if
129 
130          if (ob%sound(n)%v(k) > missing_r .AND. iv%sound(n)%v(k)%qc >= obs_qc_pointer) then
131             iv%sound(n)%v(k)%inv = ob%sound(n)%v(k) - model_v(k,n)
132          end if
133 
134          if (ob%sound(n)%t(k) > missing_r .AND. iv%sound(n)%t(k)%qc >= obs_qc_pointer) then
135             iv%sound(n)%t(k)%inv = ob%sound(n)%t(k) - model_t(k,n)
136          end if
137 
138          if (ob%sound(n)%q(k) > missing_r .AND. iv%sound(n)%q(k)%qc >= obs_qc_pointer) then
139             iv%sound(n)%q(k)%inv = ob%sound(n)%q(k) - model_q(k,n)
140          end if
141       end do
142    end do
143 
144    !----------------------------------------------------------------------
145    ! [5.0] Perform optional maximum error check:
146    !----------------------------------------------------------------------
147 
148    if (check_max_iv) call da_check_max_iv_sound (iv, it, itu,ituf,itvv,itvvf,itt,ittf,itqv,itqvf)
149 
150    if (rootproc .and. check_max_iv_print) then
151       write(unit = check_max_iv_unit, fmt ='(A,i5,A)')'For outer iteration ',it, &
152          ', Total Rejections for Sound follows:'
153 
154       write(unit = check_max_iv_unit, fmt = '(/,10(2(A,I6),/))') &
155          'Number of failed u-wind observations:     ',ituf, ' on ',itu,   &
156          'Number of failed v-wind observations:     ',itvvf,' on ',itvv,   &
157          'Number of failed temperature observations:',ittf, ' on ',itt,   &
158          'Number of failed mixing ratio observations:',itqvf,' on ',itqv, &
159          'Finally Total Sound rejections: ',ituf+itvvf+ittf+itqvf,' on ', itu+itvv+itt+itqv
160    end if
161 
162    deallocate (model_u)
163    deallocate (model_v)
164    deallocate (model_t)
165    deallocate (model_q)
166    
167    if (trace_use_dull) call da_trace_exit ("da_get_innov_vector_sound")
168 
169 end subroutine da_get_innov_vector_sound
170 
171