da_obs_sfc_correction.inc
References to this file elsewhere.
1 subroutine da_obs_sfc_correction(info, sfc_obs, n, xb)
2
3 !--------------------------------------------------------------------
4 ! Purpose: correct the surface measurements (wind,
5 ! temperature, and pressure) from the observed height to the WRF
6 ! model's lowest half-zeta level before going to the minimization.
7 !
8 ! Wind : based on the similarity theory
9 ! Temperature: Frank Ruggiero's (1996) method
10 ! Pressure : Hydrostatic correction
11 !
12 ! The order of the vertical index is "kts=1(bottom) and kte(top)".
13 ! With cv_options=2 and sfc_assi_option=1, this procedure must be
14 ! gone through, otherwise unrealistic results may be obtained.
15 !--------------------------------------------------------------------
16
17 implicit none
18
19 type(infa_type), intent(in) :: info
20 type(synop_type), intent(inout) :: sfc_obs
21 integer, intent(in) :: n
22 type(xb_type), intent(in) :: xb
23
24 real :: roughness, psfc, mslp, dx, dxm, dy, dym, ho, po, to, qo
25 real :: hm, pm, tm, qm, um, vm, correc, val
26
27 integer :: i, j, k
28 real :: t_mdl(kts:kte)
29 real :: q_mdl(kts:kte)
30 real :: u_mdl(kts:kte)
31 real :: v_mdl(kts:kte)
32 real :: height(kts:kte)
33 real :: pressure(kts:kte)
34
35 if (trace_use_dull) call da_trace_entry("da_obs_sfc_correction")
36
37 ! 1. Check if it needs to do the surface correction at the first level
38
39 ! 1.1 Surface reports located at far below the lowest model level
40
41 ! 2. Model profile at OBS site for surface correction
42
43 i = info%i(1,n)
44 j = info%j(1,n)
45 dx = info%dx(1,n)
46 dy = info%dy(1,n)
47 dxm = info%dxm(1,n)
48 dym = info%dym(1,n)
49
50 ! Model roughness at the obs site
51
52 roughness = dym*(dxm*xb%rough(i,j) + dx*xb%rough(i+1,j)) &
53 + dy *(dxm*xb%rough(i,j+1) + dx*xb%rough(i+1,j+1))
54
55 do k = kts, kte
56 pressure(k) = dym*(dxm*xb%p(i,j,k) + dx*xb%p(i+1,j,k)) + dy*(dxm*xb%p(i,j+1,k) + dx*xb%p(i+1,j+1,k))
57 height(k) = dym*(dxm*xb%h(i,j,k) + dx*xb%h(i+1,j,k)) + dy*(dxm*xb%h(i,j+1,k) + dx*xb%h(i+1,j+1,k))
58 t_mdl(k) = dym*(dxm*xb%t(i,j,k) + dx*xb%t(i+1,j,k)) + dy*(dxm*xb%t(i,j+1,k) + dx*xb%t(i+1,j+1,k))
59 q_mdl(k) = dym*(dxm*xb%q(i,j,k) + dx*xb%q(i+1,j,k)) + dy*(dxm*xb%q(i,j+1,k) + dx*xb%q(i+1,j+1,k))
60 u_mdl(k) = dym*(dxm*xb%u(i,j,k) + dx*xb%u(i+1,j,k)) + dy*(dxm*xb%u(i,j+1,k) + dx*xb%u(i+1,j+1,k))
61 v_mdl(k) = dym*(dxm*xb%v(i,j,k) + dx*xb%v(i+1,j,k)) + dy*(dxm*xb%v(i,j+1,k) + dx*xb%v(i+1,j+1,k))
62 end do
63
64 ! 3. OBS data and recover the surface pressure from the
65 ! mean sea level pressure (mslp)
66
67 ho = sfc_obs % h
68 po = sfc_obs % p % inv
69 to = sfc_obs % t % inv
70 qo = sfc_obs % q % inv
71
72 ! 3.1 Compute the surface OBS pressure from mean sea level pressure
73
74 mslp = info%slp(n)%inv
75 if (abs(mslp - missing_r) > 1.0) then
76 psfc = missing_r
77 if (abs(ho - missing_r) > 1.0) then
78 if (abs(to - missing_r) > 1.0) then
79 call da_sfcprs (kts, kte, pressure, t_mdl, q_mdl, height, psfc, mslp, ho, to)
80 else
81 call da_sfcprs (kts, kte, pressure, t_mdl, q_mdl, height, psfc, mslp, ho)
82 end if
83 end if
84 sfc_obs % p % inv = psfc
85 ! YRG: to allow assmilate the Psfc from mslp:
86 sfc_obs % p % qc = 0
87 end if
88
89 if (sfc_obs % p % inv < 1.0) then
90 sfc_obs % p % qc = missing
91 end if
92
93 po = sfc_obs % p % inv
94
95 ! 3.2 Check that obs pressure and height are present
96 ! ----------------------------------------------
97
98 if (abs(po - missing_r) < 1.0 .OR. abs(ho - missing_r) < 1.0) then
99
100 return
101
102 ! write(unit=message(1), fmt='(/3(1x,a))') &
103 ! 'MISSinG HEIGHT OR PRESSURE OBSERVATION ID ', &
104 ! trim (sfc_obs%info % id), trim (sfc_obs%info % name)
105
106 ! write(unit=message(2), fmt='(2(A,F12.3,/))') &
107 ! ' height = ',ho,&
108 ! ' pressure = ',po
109 ! call da_error(__FILE__,__LINE__,message(1:2))
110
111 end if
112
113 ! 4. Bring surface observation below model levels with good height quality
114 ! ================================================
115
116 if (sfc_obs % h < height(kts)) then
117
118 ! 2.3 Make use of local variables for model
119 ! -------------------------------------
120
121 um = u_mdl(kts)
122 vm = v_mdl(kts)
123 tm = t_mdl(kts)
124 pm = pressure (kts)
125 qm = q_mdl(kts)
126 hm = height(kts)
127
128 ! 3.2 Correction wind based on similarity laws
129 ! -------------------------------------------
130
131 if ((abs(sfc_obs%u%inv - missing_r) > 1.0) .AND. (abs(sfc_obs%v%inv - missing_r) > 1.0)) then
132 ! 3.2.1 Correction factor
133
134 ! temperature and moisture at sigma level:
135
136 if (abs(to - missing_r) < 1.0) then
137 correc = da_mo_correction(ho, po, tm, qo, hm, pm, tm, qm, um ,vm, roughness)
138 else
139 correc = da_mo_correction(ho, po, to, qo, hm, pm, tm, qm, um ,vm, roughness)
140 end if
141
142 ! 3.2.2 Wind correction
143 ! ---------------
144
145 ! Correct data and replace any previous wind qcs
146 ! with surface correction flag
147
148 sfc_obs % u % inv = correc * sfc_obs % u % inv
149 sfc_obs % u % qc = surface_correction
150
151 sfc_obs % v % inv = correc * sfc_obs % v % inv
152 sfc_obs % v % qc = surface_correction
153 end if
154
155 ! 3.4 Correct pressure
156 ! ----------------
157
158 if (sfc_obs % p % qc >= 0) then
159 ! Correct data
160 if (abs(to - missing_r) > 1.0 .and. abs(qo - missing_r) > 1.0) then
161 call da_intpsfc_prs (val, ho, po, hm, tm, qm, to, qo)
162 else if (abs(to - missing_r) > 1.0) then
163 call da_intpsfc_prs (val, ho, po, hm, tm, qm, to)
164 else
165 call da_intpsfc_prs (val, ho, po, hm, tm, qm)
166 end if
167
168 ! Replace any previous pressure qc by the surface correction
169
170 sfc_obs % p % inv = val
171 sfc_obs % p % qc = surface_correction
172 end if
173
174 ! 3.5 Correct temperature
175 ! -------------------
176
177 if (abs(sfc_obs % t % inv - missing_r) > 1.0) then
178 ! Correct data
179 call da_intpsfc_tem (val, ho, po, to, height, pressure, t_mdl, kts, kte)
180
181 sfc_obs % t % inv = val
182
183 ! Replace any previous temperature qc by the surface correction
184 sfc_obs % t % qc = surface_correction
185 end if
186
187 ! 3.6 Assign model lowest level height + 1m to observation
188 ! -----------------------------------------------------
189 ! sfc_obs % h = height(kts) + 1.0
190
191 ! 3.7 Height QC
192 ! ---------
193 ! sfc_obs % height_qc = surface_correction
194 end if
195
196 if (trace_use_dull) call da_trace_exit("da_obs_sfc_correction")
197
198 end subroutine da_obs_sfc_correction
199
200