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