da_get_innov_vector_satem.inc
References to this file elsewhere.
1 subroutine da_get_innov_vector_satem(grid, it, ob, iv)
2
3 !-----------------------------------------------------------------------
4 ! Purpose: TBD
5 !-----------------------------------------------------------------------
6
7 implicit none
8
9 type(domain), intent(in) :: grid
10 integer, intent(in) :: it ! External iteration.
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 :: i, j, k ! Index dimension.
16 integer :: num_levs ! Number of obs levels.
17 real :: dx, dxm ! Interpolation weights.
18 real :: dy, dym ! Interpolation weights.
19 real :: model_thickness(1:max_ob_levels) !Model thickness at ob loc
20 integer :: ithickness,ithicknessf
21
22 real :: pre_ma(kts-1:kte+1)
23 real :: tv_ma(kts-1:kte+1)
24 integer :: layer1,layer2
25 real :: tv1,tv2,pres2
26
27 if (trace_use_dull) call da_trace_entry("da_get_innov_vector_satem")
28
29 ithickness = 0 ; ithicknessf = 0
30
31 do n=iv%info(satem)%n1,iv%info(satem)%n2
32 num_levs = iv%info(satem)%levels(n)
33
34 if (num_levs < 1) cycle
35
36 model_thickness(:) = 0.0
37
38 ! [1.0] Get cross pt. horizontal interpolation weights:
39
40 i = iv%info(satem)%i(1,n)
41 dy = iv%info(satem)%dy(1,n)
42 dym = iv%info(satem)%dym(1,n)
43 j = iv%info(satem)%j(1,n)
44 dx = iv%info(satem)%dx(1,n)
45 dxm = iv%info(satem)%dxm(1,n)
46
47 !------------------------------------------------------------------------
48
49 ! [2.0] Calculate vertical profile of virtual temperature at obs pt.
50
51 call da_tv_profile(grid,iv%info(satem),n,pre_ma,tv_ma)
52
53 ! [3.0] Find model vertical position of pressure and do interp.
54
55 call da_find_layer(layer2,tv2,iv%satem(n)%ref_p,pre_ma,tv_ma,kts,kte)
56 pres2 = iv%satem(n)%ref_p
57
58 ! [4.0] Thickness innovations calculation
59
60 do k = 1, num_levs
61 iv % satem(n) % thickness(k) % inv = 0.0
62 call da_find_layer(layer1,tv1,iv%satem(n)%p(k),pre_ma,tv_ma,kts,kte)
63 call da_thickness(pre_ma,tv_ma,kts,kte,tv1,tv2,layer1,layer2, &
64 iv%satem(n)%p(k),pres2,model_thickness(k))
65
66 if (ABS(ob % satem(n) % thickness(k) - missing_r) > 1.0 .and. &
67 iv % satem(n) % thickness(k)%qc /= missing_data) then
68 iv % satem(n) % thickness(k) % inv = ob % satem(n) % thickness(k) - model_thickness(k)
69 ! write(unit=stdout,fmt="(A, 2F10.3,F10.0,A,F5.0,A)") &
70 ! "observed, model_thickness, layer = ", &
71 ! ob%satem(n)%thickness(k), &
72 ! model_thickness(k), 0.01*pres2, " -", &
73 ! iv%satem(n)%p(k)*0.01,'hPa'
74 end if
75
76 pres2 = iv%satem(n)%p(k)
77 layer2 = layer1
78 tv2 = tv1
79 end do
80 end do
81
82 !------------------------------------------------------------------------
83 ! [5.0] Perform optional maximum error check:
84 !------------------------------------------------------------------------
85
86 if (check_max_iv) then
87 call da_check_max_iv_satem(iv, it, ithickness,ithicknessf)
88 end if
89
90 !------------------------------------------------------------------------
91 ! [6.0] Perform land/ocean check
92 !------------------------------------------------------------------------
93
94 do n=iv%info(satem)%n1,iv%info(satem)%n2
95 i = iv%info(satem)%i(1,n)
96 j = iv%info(satem)%j(1,n)
97 if (grid%xb%landmask(i,j ) /= 0.0 .or. grid%xb%landmask(i+1,j ) /= 0. .or. &
98 grid%xb%landmask(i,j+1) /= 0.0 .or. grid%xb%landmask(i+1,j+1) /= 0.0) then
99 iv % satem(n) % thickness(1) % inv = 0.0
100 ithicknessf = ithicknessf + 1
101 end if
102 end do
103
104 if (rootproc .and. check_max_iv_print) then
105 write(unit = check_max_iv_unit, fmt ='(A,i5,A)')'For outer iteration ',it, &
106 ', Total Rejections for Satem follows:'
107 write(unit = check_max_iv_unit, fmt = '(/,10(2(A,I6),/))') &
108 'Number of failed thicknesses observations: ',&
109 ithicknessf, ' on ',ithickness
110 end if
111
112 if (trace_use_dull) call da_trace_exit("da_get_innov_vector_satem")
113
114 end subroutine da_get_innov_vector_satem
115
116