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 > 0) then
28 ithickness = 0 ; ithicknessf = 0
29
30 do n=iv%ob_numb(iv%current_ob_time-1)%satem + 1, iv%ob_numb(iv%current_ob_time)%satem
31 num_levs = iv % satem(n) % info % levels
32
33 if (num_levs < 1) cycle
34
35 model_thickness(:) = 0.0
36
37 ! [1.0] Get cross pt. horizontal interpolation weights:
38
39 i = iv%satem(n)%loc%i
40 dy = iv%satem(n)%loc%dy
41 dym = iv%satem(n)%loc%dym
42 j = iv%satem(n)%loc%j
43 dx = iv%satem(n)%loc%dx
44 dxm = iv%satem(n)%loc%dxm
45 ks = xp%kts; ke=xp%kte
46
47 !------------------------------------------------------------------------
48
49 ! [2.0] Calculate vertical profile of virtual temperature at obs pt.
50
51 call da_tv_profile(xp,xb,i,j,dx,dxm,dy,dym,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,ks,ke)
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,ks,ke)
63 call da_thickness(pre_ma,tv_ma,ks,ke,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. .and. &
67 iv % satem(n) % thickness(k)%qc /= missing_data) then
68 iv % satem(n) % thickness(k) % inv = &
69 ob % satem(n) % thickness(k) - model_thickness(k)
70 ! write(unit=stdout,fmt="(A, 2F10.3,F10.0,A,F5.0,A)") &
71 ! "observed, model_thickness, layer = ", &
72 ! ob%satem(n)%thickness(k), &
73 ! model_thickness(k), 0.01*pres2, " -", &
74 ! iv%satem(n)%p(k)*0.01,'hPa'
75 end if
76
77 pres2 = iv%satem(n)%p(k)
78 layer2 = layer1
79 tv2 = tv1
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(it, iv % satem(n), ithickness,ithicknessf)
88 end if
89
90 !------------------------------------------------------------------------
91 ! [6.0] Perform land/ocean check
92 !------------------------------------------------------------------------
93
94 if (xb%landmask(i,j ) /= 0. .or. xb%landmask(i+1,j ) /= 0. .or. &
95 xb%landmask(i,j+1) /= 0. .or. xb%landmask(i+1,j+1) /= 0.) then
96 iv % satem(n) % thickness(1) % inv = 0.
97 ithicknessf = ithicknessf + 1
98 end if
99 end do
100
101 if (rootproc .and. check_max_iv_print) then
102 write(unit = check_max_iv_unit, fmt ='(A,i5,A)')'For outer iteration ',it, &
103 ', Total Rejections for Satem follows:'
104 write(unit = check_max_iv_unit, fmt = '(/,10(2(A,I6),/))') &
105 'Number of failed thicknesses observations: ',&
106 ithicknessf, ' on ',ithickness
107 end if
108 end if
109
110 end subroutine da_get_innov_vector_satem
111
112