da_get_innov_vector_satem.inc
References to this file elsewhere.
1 subroutine da_get_innov_vector_satem( it, xb, xp, ob, iv)
2
3 !-----------------------------------------------------------------------
4 ! Purpose: TBD
5 !-----------------------------------------------------------------------
6
7 implicit none
8
9 integer, intent(in) :: it ! External iteration.
10 type(xb_type), intent(in) :: xb ! first guess state.
11 type(xpose_type), intent(in) :: xp ! Domain decomposition vars.
12 type(y_type), intent(in) :: ob ! Observation structure.
13 type(ob_type), intent(inout) :: iv ! O-B structure.
14
15 integer :: n ! Loop counter.
16 integer :: i, j, k ! Index dimension.
17 integer :: num_levs ! Number of obs levels.
18 real :: dx, dxm ! Interpolation weights.
19 real :: dy, dym ! Interpolation weights.
20 real, dimension(1:max_ob_levels) :: model_thickness !Model thickness at ob loc
21 integer :: ithickness,ithicknessf
22
23 real, dimension(xp%kts-1:xp%kte+1) :: pre_ma,tv_ma
24 integer :: layer1,layer2,ks,ke
25 real :: tv1,tv2,pres2
26
27 if (iv % num_satem < 1) return
28
29 if (trace_use) call da_trace_entry("da_get_innov_vector_satem")
30
31 ithickness = 0 ; ithicknessf = 0
32
33 do n=iv%ob_numb(iv%current_ob_time-1)%satem + 1, iv%ob_numb(iv%current_ob_time)%satem
34 num_levs = iv % satem(n) % info % levels
35
36 if (num_levs < 1) cycle
37
38 model_thickness(:) = 0.0
39
40 ! [1.0] Get cross pt. horizontal interpolation weights:
41
42 i = iv%satem(n)%loc%i
43 dy = iv%satem(n)%loc%dy
44 dym = iv%satem(n)%loc%dym
45 j = iv%satem(n)%loc%j
46 dx = iv%satem(n)%loc%dx
47 dxm = iv%satem(n)%loc%dxm
48 ks = xp%kts; ke=xp%kte
49
50 !------------------------------------------------------------------------
51
52 ! [2.0] Calculate vertical profile of virtual temperature at obs pt.
53
54 call da_tv_profile(xp,xb,i,j,dx,dxm,dy,dym,pre_ma,tv_ma)
55
56 ! [3.0] Find model vertical position of pressure and do interp.
57
58 call da_find_layer(layer2,tv2,iv%satem(n)%ref_p,pre_ma,tv_ma,ks,ke)
59 pres2 = iv%satem(n)%ref_p
60
61 ! [4.0] Thickness innovations calculation
62
63 do k = 1, num_levs
64 iv % satem(n) % thickness(k) % inv = 0.0
65 call da_find_layer(layer1,tv1,iv%satem(n)%p(k),pre_ma,tv_ma,ks,ke)
66 call da_thickness(pre_ma,tv_ma,ks,ke,tv1,tv2,layer1,layer2, &
67 iv%satem(n)%p(k),pres2,model_thickness(k))
68
69 if (ABS(ob % satem(n) % thickness(k) - missing_r) > 1. .and. &
70 iv % satem(n) % thickness(k)%qc /= missing_data) then
71 iv % satem(n) % thickness(k) % inv = &
72 ob % satem(n) % thickness(k) - model_thickness(k)
73 ! write(unit=stdout,fmt="(A, 2F10.3,F10.0,A,F5.0,A)") &
74 ! "observed, model_thickness, layer = ", &
75 ! ob%satem(n)%thickness(k), &
76 ! model_thickness(k), 0.01*pres2, " -", &
77 ! iv%satem(n)%p(k)*0.01,'hPa'
78 end if
79
80 pres2 = iv%satem(n)%p(k)
81 layer2 = layer1
82 tv2 = tv1
83 end do
84
85 !------------------------------------------------------------------------
86 ! [5.0] Perform optional maximum error check:
87 !------------------------------------------------------------------------
88
89 if (check_max_iv) then
90 call da_check_max_iv_satem(it, iv % satem(n), ithickness,ithicknessf)
91 end if
92
93 !------------------------------------------------------------------------
94 ! [6.0] Perform land/ocean check
95 !------------------------------------------------------------------------
96
97 if (xb%landmask(i,j ) /= 0. .or. xb%landmask(i+1,j ) /= 0. .or. &
98 xb%landmask(i,j+1) /= 0. .or. xb%landmask(i+1,j+1) /= 0.) then
99 iv % satem(n) % thickness(1) % inv = 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) call da_trace_exit("da_get_innov_vector_satem")
113
114 end subroutine da_get_innov_vector_satem
115
116