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