da_sfc_wtq.inc
References to this file elsewhere.
1 subroutine da_sfc_wtq (psfc, tg, ps, ts, qs, us, vs, ps2, ts2, qs2, &
2 hs, roughness, xland, u10, v10, t2, q2, regime)
3
4 !---------------------------------------------------------------------------
5 ! Purpose: Calculate the 10m wind, 2m temperature and moisture based on the
6 ! similarity theory/
7 !
8 ! The unit for pressure : psfc, ps, ps2 is Pa.
9 ! The unit for temperature: tg, ts, ts2, t2 is K.
10 ! The unit for moisture : qs, qs2, q2 is kg/kg.
11 ! The unit for wind : us, vs, u10, v10 is m/s.
12 ! The unit for height : hs, roughness is m.
13 ! xland and regime are dimensionless.
14 !
15 ! Reference:
16 ! ---------
17 !
18 ! input Variables:
19 !
20 ! psfc, tg : surface pressure and ground temperature
21 ! ps, ts, qs, us, vs, hs : model variable at lowlest half sigma level
22 ! ps2, ts2, qs2 : model variables at the second lowest half
23 ! sigma level
24 !
25 !
26 ! Constants:
27 !
28 ! hs : height at the lowest half sigma level
29 ! roughness : roughness
30 ! xland : land-water-mask (=2 water, =1 land)
31 !
32 ! output Variables:
33 !
34 ! regime : PBL regime
35 ! u10, v10 : 10-m high observed wind components
36 ! t2 , q2 : 2-m high observed temperature and mixing ratio
37 !
38 !---------------------------------------------------------------------------
39 !
40 ! psim : mechanical psi at lowlest sigma level
41 ! psim2 : mechanical psi at 2m
42 ! psimz : mechanical psi at 10m
43 !
44 !---------------------------------------------------------------------------
45
46 implicit none
47
48 real, intent (in) :: ps , ts , qs , us, vs
49 real, intent (in) :: ps2, ts2, qs2, psfc, tg
50 real, intent (in) :: hs, roughness, xland
51 real, intent (out) :: regime
52 real, intent (out) :: u10, v10, t2, q2
53
54 ! Maximum number of iterations in computing psim, psih
55
56 integer, parameter :: k_iteration = 10
57 ! integer, parameter :: k_iteration = 1
58
59 ! h10 is the height of 10m where the wind observed
60 ! h2 is the height of 2m where the temperature and
61 ! moisture observed.
62
63 real, parameter :: h10 = 10.0, h2 = 2.0
64
65 ! Default roughness over the land
66
67 real, parameter :: zint0 = 0.01
68
69 ! Von Karman constant
70
71 real, parameter :: k_kar = 0.4
72
73 ! Working variables
74
75 real :: Vc2, Va2, V2
76 real :: rib, rcp, xx, yy, cc
77 real :: psiw, psiz, mol, ust, hol, holz, hol2
78 real :: psim, psimz, psim2, psih, psihz, psih2
79 real :: psit, psit2, psiq, psiq2
80 real :: gzsoz0, gz10oz0, gz2oz0
81 real :: eg, qg, tvg, tvs, tvs2
82 real :: ths, thg, thvs, thvg, thvs2
83 real :: zq0, z0
84
85 real, parameter :: ka = 2.4E-5
86
87 if (trace_use_dull) call da_trace_entry("da_sfc_wtq")
88
89 rcp = gas_constant/cp
90
91 ! 1 Compute the roughness length based upon season and land use
92
93 ! 1.1 Define the roughness length
94
95 z0 = roughness
96
97 if (z0 < 0.0001) z0 = 0.0001
98
99 ! 1.2 Define the rouhgness length for moisture
100
101 if (xland .ge. 1.5) then
102 zq0 = z0
103 else
104 zq0 = zint0
105 end if
106
107 ! 1.3 Define the some constant variable for psi
108
109 gzsoz0 = log(hs/z0)
110
111 gz10oz0 = log(h10/z0)
112
113 gz2oz0 = log(h2/z0)
114
115
116 ! 2. Calculate the virtual temperature
117
118 ! 2.1 Compute Virtual temperature on the lowest half sigma level
119
120 tvs = ts * (1.0 + 0.608 * qs)
121
122 ! 2.2 Compute Virtual temperature on the second lowest half sigma level
123
124 tvs2 = ts2 * (1.0 + 0.608 * qs2)
125
126 ! 2.3 Convert ground virtual temperature assuming it's saturated
127
128 call da_tp_to_qs(tg, psfc, eg, qg)
129 tvg = tg * (1.0 + 0.608 * qg)
130
131 ! 3. Compute the potential temperature
132
133 ! 3.1 Potential temperature on the lowest half sigma level
134
135 ths = ts * (1000.0 / (ps/100.0)) ** rcp
136
137 ! 3.2 Potential temperature at the ground
138
139 thg = tg * (1000.0 / (psfc/100.0)) ** rcp
140
141 ! 4. Virtual potential temperature
142
143 ! 4.1 Virtual potential temperature on the lowest half sigma level
144
145 thvs = tvs * (1000.0 / (ps/100.0)) ** rcp
146
147 ! 4.2 Virtual potential temperature on the second lowest half sigma level
148
149 thvs2 = tvs2 * (1000.0 / (ps2/100.0)) ** rcp
150
151 ! 4.3 Virtual potential temperature at ground
152
153 thvg = tvg * (1000.0 / (psfc/100.0)) ** rcp
154
155
156 ! 5. BULK RICHARDSON NUMBER AND MONI-OBUKOV LENGTH
157
158 ! 5.1 Velocity
159
160 ! Wind speed:
161
162 Va2 = us*us + vs*vs
163 !
164 ! Convective velocity:
165
166 if (thvg >= thvs) then
167 Vc2 = 4.0 * (thvg - thvs)
168 else
169 Vc2 = 0.0
170 end if
171
172 V2 = Va2 + Vc2
173
174 ! 5.2 Bulk richardson number
175
176 rib = (gravity * hs / ths) * (thvs - thvg) / V2
177
178 ! 6. CALCULATE PSI BASED UPON REGIME
179
180
181 if (rib .GE. 0.2) then
182 ! 6.1 Stable conditions (REGIME 1)
183 ! ---------------------------
184 regime = 1.1
185 psim = -10.0*gzsoz0
186 psimz = -10.0*gz10oz0
187 psim2 = -10.0*gz2oz0
188 psim = max(psim,-10.0)
189 psimz = max(psimz,-10.0)
190 psim2 = max(psim2,-10.0)
191 psih = psim
192 psihz = psimz
193 psih2 = psim2
194
195 else if ((rib .LT. 0.2) .AND. (rib .GT. 0.0)) then
196 ! 6.2 Mechanically driven turbulence (REGIME 2)
197
198 regime = 2.1
199 psim = (-5.0 * rib) * gzsoz0 / (1.1 - 5.0*rib)
200 psimz = (-5.0 * rib) * gz10oz0 / (1.1 - 5.0*rib)
201 psim2 = (-5.0 * rib) * gz2oz0 / (1.1 - 5.0*rib)
202
203 psim = max(psim,-10.0)
204 psimz = max(psimz,-10.0)
205 psim2 = max(psim2,-10.0)
206 psih = psim
207 psihz = psimz
208 psih2 = psim2
209
210 else if ((rib .EQ. 0.0) .or. (rib.LT.0.0 .and. thvs2.GT.thvs)) then
211 ! 6.3 Unstable Forced convection (REGIME 3)
212
213 regime = 3.1
214 psim = 0.0
215 psimz = 0.0
216 psim2 = 0.0
217 psih = psim
218 psihz = psimz
219 psih2 = psim2
220
221 else
222 ! 6.4 Free convection (REGIME 4)
223 regime = 4.1
224
225 ! Calculate psi m and pshi h using iteration method
226
227 psim = 0.0
228 psih = 0.0
229 cc = 2.0 * atan(1.0)
230
231 ! do k = 1 , k_iteration
232
233 ! 6.4.1 Calculate ust, m/L (mol), h/L (hol)
234
235 ! Friction speed
236
237 ust = k_kar * sqrt(v2) /(gzsoz0 - psim)
238
239 ! Heat flux factor
240
241 mol = k_kar * (ths - thg)/(gzsoz0 - psih)
242
243 ! Ratio of PBL height to Monin-Obukhov length
244
245 if (ust .LT. 0.01) then
246 hol = rib * gzsoz0
247 else
248 hol = k_kar * gravity * hs * mol / (ths * ust * ust)
249 end if
250
251 ! 6.4.2 Calculate n, nz, R, Rz
252
253 hol = min(hol,0.0)
254 hol = max(hol,-10.0)
255
256 holz = (h10 / hs) * hol
257 holz = min(holz,0.0)
258 holz = max(holz,-10.0)
259
260 hol2 = (h2 / hs) * hol
261 hol2 = min(hol2,0.0)
262 hol2 = max(hol2,-10.0)
263
264 ! 6.4.3 Calculate Psim & psih
265
266 ! Using the look-up table:
267 ! nn = int(-100.0 * hol)
268 ! rr = (-100.0 * hol) - nn
269 ! psim = psimtb(nn) + rr * (psimtb(nn+1) - psimtb(nn))
270 ! psih = psihtb(nn) + rr * (psihtb(nn+1) - psihtb(nn))
271 ! Using the continuous function:
272 xx = (1.0 - 16.0 * hol) ** 0.25
273 yy = log((1.0+xx*xx)/2.0)
274 psim = 2.0 * log((1.0+xx)/2.0) + yy - 2.0 * atan(xx) + cc
275 psih = 2.0 * yy
276
277 ! Using the look-up table:
278 ! nz = int(-100.0 * holz)
279 ! rz = (-100.0 * holz) - nz
280 ! psimz = psimtb(nz) + rz * (psimtb(nz+1) - psimtb(nz))
281 ! psihz = psihtb(nz) + rz * (psihtb(nz+1) - psihtb(nz))
282 ! Using the continuous function:
283 xx = (1.0 - 16.0 * holz) ** 0.25
284 yy = log((1.0+xx*xx)/2.0)
285 psimz = 2.0 * log((1.0+xx)/2.0) + yy - 2.0 * atan(xx) + cc
286 psihz = 2.0 * yy
287
288 ! Using the look-up table:
289 ! n2 = int(-100.0 * hol2)
290 ! r2 = (-100.0 * hol2) - n2
291 ! psim2 = psimtb(n2) + r2 * (psimtb(n2+1) - psimtb(n2))
292 ! psih2 = psihtb(n2) + r2 * (psihtb(n2+1) - psihtb(n2))
293 ! Using the continuous function:
294 xx = (1.0 - 16.0 * hol2) ** 0.25
295 yy = log((1.0+xx*xx)/2.0)
296 psim2 = 2.0 * log((1.0+xx)/2.0) + yy - 2.0 * atan(xx) + cc
297 psih2 = 2.0 * yy
298
299 ! end do
300
301 ! 6.4.4 Define the limit value for psim & psih
302
303 psim = min(psim,0.9*gzsoz0)
304 psimz = min(psimz,0.9*gz10oz0)
305 psim2 = min(psim2,0.9*gz2oz0)
306 psih = min(psih,0.9*gzsoz0)
307 psihz = min(psihz,0.9*gz10oz0)
308 psih2 = min(psih2,0.9*gz2oz0)
309 end if ! Regime
310
311 ! 7. Calculate psi for wind, temperature and moisture
312
313 psiw = gzsoz0 - psim
314 psiz = gz10oz0 - psimz
315 psit = gzsoz0 - psih
316 psit2 = gz2oz0 - psih2
317
318 ! Friction speed
319 ust = k_kar * sqrt(v2) /(gzsoz0 - psim)
320
321 psiq = log(k_kar*ust*hs/ka + hs / zq0) - psih
322 psiq2 = log(k_kar*ust*h2/ka + h2 / zq0) - psih2
323
324 ! 8. Calculate 10m wind, 2m temperature and moisture
325
326 u10 = us * psiz / psiw
327 v10 = vs * psiz / psiw
328 t2 = (thg + (ths - thg)*psit2/psit)*((psfc/100.0)/1000.0)**rcp
329 q2 = qg + (qs - qg)*psiq2/psiq
330
331 if (trace_use_dull) call da_trace_exit("da_sfc_wtq")
332
333 end subroutine da_sfc_wtq
334
335