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, &
3 u10, v10, t2, q2, regime)
4
5 !---------------------------------------------------------------------------
6 ! Purpose: Calculate the 10m wind, 2m temperature and moisture based on the
7 ! similarity theory/
8 !
9 ! The unit for pressure : psfc, ps, ps2 is Pa.
10 ! The unit for temperature: tg, ts, ts2, t2 is K.
11 ! The unit for moisture : qs, qs2, q2 is kg/kg.
12 ! The unit for wind : us, vs, u10, v10 is m/s.
13 ! The unit for height : hs, roughness is m.
14 ! xland and regime are dimensionless.
15 !
16 ! Reference:
17 ! ---------
18 !
19 ! input Variables:
20 !
21 ! psfc, tg : surface pressure and ground temperature
22 ! ps, ts, qs, us, vs, hs : model variable at lowlest half sigma level
23 ! ps2, ts2, qs2 : model variables at the second lowest half
24 ! sigma level
25 !
26 !
27 ! Constants:
28 !
29 ! hs : height at the lowest half sigma level
30 ! roughness : roughness
31 ! xland : land-water-mask (=2 water, =1 land)
32 !
33 ! output Variables:
34 !
35 ! regime : PBL regime
36 ! u10, v10 : 10-m high observed wind components
37 ! t2 , q2 : 2-m high observed temperature and mixing ratio
38 !
39 !---------------------------------------------------------------------------
40 !
41 ! psim : mechanical psi at lowlest sigma level
42 ! psim2 : mechanical psi at 2m
43 ! psimz : mechanical psi at 10m
44 !
45 !---------------------------------------------------------------------------
46
47 implicit none
48
49 real, intent (in) :: ps , ts , qs , us, vs
50 real, intent (in) :: ps2, ts2, qs2, psfc, tg
51 real, intent (in) :: hs, roughness, xland
52
53 real, intent (out) :: regime
54 real , intent (out) :: u10, v10, t2, q2
55
56 ! Maximum number of iterations in computing psim, psih
57
58 integer, parameter :: k_iteration = 10
59 ! integer, parameter :: k_iteration = 1
60
61 ! h10 is the height of 10m where the wind observed
62 ! h2 is the height of 2m where the temperature and
63 ! moisture observed.
64
65 real, parameter :: h10 = 10., h2 = 2.
66
67 ! Default roughness over the land
68
69 real, parameter :: zint0 = 0.01
70
71 ! Von Karman constant
72
73 real, parameter :: k_kar = 0.4
74
75 ! Working variables
76
77 real :: Vc2, Va2, V2
78 real :: rib, rcp, xx, yy, cc
79 real :: psiw, psiz, mol, ust, hol, holz, hol2
80 real :: psim, psimz, psim2, psih, psihz, psih2
81 real :: psit, psit2, psiq, psiq2
82 real :: gzsoz0, gz10oz0, gz2oz0
83 real :: eg, qg, tvg, tvs, tvs2
84 real :: ths, thg, thvs, thvg, thvs2
85 real :: zq0, z0
86
87 real, parameter :: ka = 2.4E-5
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.608 * qs)
121
122 ! 2.2 Compute Virtual temperature on the second lowest half sigma level
123
124 tvs2 = ts2 * (1. + 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.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. / (ps/100.)) ** rcp
136
137 ! 3.2 Potential temperature at the ground
138
139 thg = tg * (1000. / (psfc/100.)) ** 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. / (ps/100.)) ** rcp
146
147 ! 4.2 Virtual potential temperature on the second lowest half sigma level
148
149 thvs2 = tvs2 * (1000. / (ps2/100.)) ** rcp
150
151 ! 4.3 Virtual potential temperature at ground
152
153 thvg = tvg * (1000. / (psfc/100.)) ** 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. * (thvg - thvs)
168 else
169 Vc2 = 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.*gzsoz0
186 psimz = -10.*gz10oz0
187 psim2 = -10.*gz2oz0
188 psim = max(psim,-10.)
189 psimz = max(psimz,-10.)
190 psim2 = max(psim2,-10.)
191 psih = psim
192 psihz = psimz
193 psih2 = psim2
194
195 else if ((rib .LT. 0.2) .AND. (rib .GT. 0.)) then
196 ! 6.2 Mechanically driven turbulence (REGIME 2)
197
198 regime = 2.1
199 psim = (-5. * rib) * gzsoz0 / (1.1 - 5.*rib)
200 psimz = (-5. * rib) * gz10oz0 / (1.1 - 5.*rib)
201 psim2 = (-5. * rib) * gz2oz0 / (1.1 - 5.*rib)
202
203 psim = max(psim,-10.)
204 psimz = max(psimz,-10.)
205 psim2 = max(psim2,-10.)
206 psih = psim
207 psihz = psimz
208 psih2 = psim2
209
210 else if ((rib .EQ. 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.
215 psimz = 0.
216 psim2 = 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.
228 psih = 0.
229 cc = 2. * 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.)
254 hol = max(hol,-10.)
255
256 holz = (h10 / hs) * hol
257 holz = min(holz,0.)
258 holz = max(holz,-10.)
259
260 hol2 = (h2 / hs) * hol
261 hol2 = min(hol2,0.)
262 hol2 = max(hol2,-10.)
263
264 ! 6.4.3 Calculate Psim & psih
265
266 ! Using the look-up table:
267 ! nn = int(-100. * hol)
268 ! rr = (-100. * 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. - 16. * hol) ** 0.25
273 yy = log((1.+xx*xx)/2.)
274 psim = 2. * log((1.+xx)/2.) + yy - 2. * atan(xx) + cc
275 psih = 2. * yy
276
277 ! Using the look-up table:
278 ! nz = int(-100. * holz)
279 ! rz = (-100. * 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. - 16. * holz) ** 0.25
284 yy = log((1.+xx*xx)/2.)
285 psimz = 2. * log((1.+xx)/2.) + yy - 2. * atan(xx) + cc
286 psihz = 2. * yy
287
288 ! Using the look-up table:
289 ! n2 = int(-100. * hol2)
290 ! r2 = (-100. * 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. - 16. * hol2) ** 0.25
295 yy = log((1.+xx*xx)/2.)
296 psim2 = 2. * log((1.+xx)/2.) + yy - 2. * atan(xx) + cc
297 psih2 = 2. * 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.)/1000.)**rcp
329 q2 = qg + (qs - qg)*psiq2/psiq
330
331 end subroutine da_sfc_wtq
332
333