da_get_innov_vector_geoamv.inc
References to this file elsewhere.
1 subroutine da_get_innov_vector_geoamv( it, grid, ob, iv)
2
3 !-------------------------------------------------------------------------
4 ! Purpose: Calculates innovation vector does QC for geoamv
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(in) :: ob ! Observation structure.
12 type(iv_type), intent(inout) :: iv ! O-B structure.
13
14 integer :: n ! Loop counter.
15 integer :: k ! Index dimension.
16 integer :: num_levs ! Number of obs levels.
17
18 integer :: i (kms:kme)
19 integer :: j (kms:kme)
20 real :: dx (kms:kme)
21 real :: dxm(kms:kme)
22 real :: dy (kms:kme)
23 real :: dym(kms:kme)
24 real,allocatable :: model_u(:,:)
25 real,allocatable :: model_v(:,:)
26
27 real :: v_p(kts:kte) ! Model value p at ob hor. location.
28
29 integer :: itu,ituf,itvv,itvvf
30
31 if (trace_use_dull) call da_trace_entry ("da_get_innov_vector_geoamv")
32
33 itu = 0
34 itvv = 0
35 ituf = 0
36 itvvf = 0
37
38 allocate (model_u(iv%info(geoamv)%max_lev,iv%info(geoamv)%n1:iv%info(geoamv)%n2))
39 allocate (model_v(iv%info(geoamv)%max_lev,iv%info(geoamv)%n1:iv%info(geoamv)%n2))
40
41 model_u(:,:) = 0.0
42 model_v(:,:) = 0.0
43
44 do n = iv%info(geoamv)%n1, iv%info(geoamv)%n2
45 ! [1.3] Get horizontal interpolation weights:
46
47 num_levs = iv%info(geoamv)%levels(n)
48 if (num_levs < 1) cycle
49
50 ! slower
51 ! i(:) = iv%info(geoamv)%i(:,n)
52 ! j(:) = iv%info(geoamv)%j(:,n)
53 ! dx(:) = iv%info(geoamv)%dx(:,n)
54 ! dy(:) = iv%info(geoamv)%dy(:,n)
55 ! dxm(:) = iv%info(geoamv)%dxm(:,n)
56 ! dym(:) = iv%info(geoamv)%dym(:,n)
57
58 ! faster
59 i(1) = iv%info(geoamv)%i(1,n)
60 j(1) = iv%info(geoamv)%j(1,n)
61 dx(1) = iv%info(geoamv)%dx(1,n)
62 dy(1) = iv%info(geoamv)%dy(1,n)
63 dxm(1) = iv%info(geoamv)%dxm(1,n)
64 dym(1) = iv%info(geoamv)%dym(1,n)
65
66 ! if position varies with height, slower
67 ! do k=kts,kte
68 ! 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)) &
69 ! + 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))
70 ! end do
71
72 ! If position does not, faster
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
76 do k=1, iv%info(geoamv)%levels(n)
77 if (iv%geoamv(n)%p(k) > 1.0) then
78 call da_to_zk (iv%geoamv(n)%p(k), v_p, v_interp_p, iv%info(geoamv)%zk(k,n))
79 end if
80 end do
81
82 end do
83
84 call da_convert_zk (iv%info(geoamv))
85
86 if (.not. anal_type_verify) then
87 do n = iv%info(geoamv)%n1, iv%info(geoamv)%n2
88 do k=1, iv%info(geoamv)%levels(n)
89 if (iv%info(geoamv)%zk(k,n) < 0.0) then
90 iv%geoamv(n)%u(k)% qc = missing
91 iv%geoamv(n)%v(k)% qc = missing
92 end if
93 end do
94 end do
95 end if
96
97 call da_interp_lin_3d (grid%xb%u, iv%info(geoamv), model_u)
98 call da_interp_lin_3d (grid%xb%v, iv%info(geoamv), model_v)
99
100 do n = iv%info(geoamv)%n1, iv%info(geoamv)%n2
101 do k = 1, iv%info(geoamv)%levels(n)
102 iv%geoamv(n)%u(k)%inv = 0.0
103 iv%geoamv(n)%v(k)%inv = 0.0
104 if (ob%geoamv(n)%u(k) > missing_r .AND. iv%geoamv(n)%u(k)%qc >= obs_qc_pointer) then
105 iv%geoamv(n)%u(k)%inv = ob%geoamv(n)%u(k) - model_u(k,n)
106 end if
107
108 if (ob%geoamv(n)%v(k) > missing_r .AND. iv%geoamv(n)%v(k)%qc >= obs_qc_pointer) then
109 iv%geoamv(n)%v(k)%inv = ob%geoamv(n)%v(k) - model_v(k,n)
110 end if
111 end do
112 end do
113
114 deallocate (model_u)
115 deallocate (model_v)
116
117 !------------------------------------------------------------------------
118 ! Perform optional maximum error check:
119 !------------------------------------------------------------------------
120
121 if (check_max_iv) then
122 call da_check_max_iv_geoamv (iv,it,itu,ituf,itvv,itvvf)
123 end if
124
125 if (rootproc .and. check_max_iv_print) then
126 write(unit = check_max_iv_unit, fmt ='(A,i5,A)')'For outer iteration ',it, &
127 ', Total Rejections for Geo. AMVs follows:'
128 write(unit = check_max_iv_unit, fmt = '(/,10(2(A,I6),/))') &
129 'Number of failed u-wind observations: ',ituf, ' on ',itu, &
130 'Number of failed v-wind observations: ',itvvf,' on ',itvv, &
131 'Finally Total Geo. AMVs rejections ',ituf+itvvf,' on ',itu+itvv
132 end if
133
134
135 if (trace_use_dull) call da_trace_exit ("da_get_innov_vector_geoamv")
136
137 end subroutine da_get_innov_vector_geoamv
138
139