da_sfc_wtq_lin.inc
References to this file elsewhere.
1 subroutine da_sfc_wtq_lin(psfc, tg, ps, ts, qs, us, vs, regime, &
2 psfc_prime,tg_prime,ps_prime,ts_prime,qs_prime, &
3 us_prime,vs_prime, hs,roughness,xland, &
4 u10_prime,v10_prime,t2_prime,q2_prime)
5
6 !---------------------------------------------------------------------------
7 ! Purpose: Calculate the 10m wind, 2m temperature and moisture based on the
8 ! similarity theory/
9 !
10 ! Reference:
11 ! ---------
12 !
13 ! input Variables(basic state):
14 !
15 ! psfc, tg : surface pressure and ground temperature
16 ! ps, ts, qs, us, vs, hs : model variable at lowlest half sigma level
17 ! regime : PBL regime
18 !
19 ! input Variables(pertubation):
20 !
21 ! psfc_prime, tg_prime : Surface pressure and ground temperature
22 ! ps_prime, ts_prime, : Model variables at the lowest half sigma
23 ! qs_prime, us_prime, : level
24 ! vs_prime :
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(pertubation):
33 !
34 ! u10_prime, v10_prime : 10-m high observed wind components
35 ! t2_prime , q2_prime : 2-m high observed temperature and mixing ratio
36 !
37 !---------------------------------------------------------------------------
38 !
39 ! psim : mechanical psi at lowlest sigma level
40 ! psim2 : mechanical psi at 2m
41 ! psimz : mechanical psi at 10m
42 !
43 !---------------------------------------------------------------------------
44
45 implicit none
46
47 real , intent (in) :: regime
48 real , intent (in) :: ps , ts , qs , us, vs, psfc, tg
49 real , intent (in) :: ps_prime, ts_prime, qs_prime , &
50 us_prime, vs_prime, psfc_prime, tg_prime
51 real , intent (in) :: hs, roughness, xland
52
53 real , intent (out) :: u10_prime, v10_prime, t2_prime, q2_prime
54
55 ! Maximum number of iterations in computing psim, psih
56
57 integer, parameter :: k_iteration = 10
58 ! integer, parameter :: k_iteration = 1
59
60 ! h10 is the height of 10m where the wind observed
61 ! h2 is the height of 2m where the temperature and
62 ! moisture observed.
63
64 real, parameter :: h10 = 10., h2 = 2.
65 !
66 ! Default roughness over the land
67
68 real, parameter :: zint0 = 0.01
69 !
70 ! Von Karman constant
71
72 real, parameter :: k_kar = 0.4
73 !
74 ! Working variables
75
76 real :: Vc2, Va2, V2
77 real :: rib, rcp, xx, yy, cc, Pi
78 real :: psiw, psiz, mol, ust, hol, holz, hol2
79 real :: psim, psimz, psim2, psih, psihz, psih2
80 real :: psit, psit2, psiq, psiq2
81 real :: gzsoz0, gz10oz0, gz2oz0
82 real :: eg, qg, tvg, tvs
83 real :: ths, thg, thvs, thvg
84 real :: zq0, z0
85
86 real :: Vc2_prime, Va2_prime, V2_prime
87 real :: rib_prime, xx_prime, yy_prime
88 real :: psiw_prime, psiz_prime, mol_prime, ust_prime, &
89 hol_prime, holz_prime, hol2_prime
90 real :: psim_prime, psimz_prime, psim2_prime, &
91 psih_prime, psihz_prime, psih2_prime
92 real :: psit_prime, psit2_prime, &
93 psiq_prime, psiq2_prime
94 real :: qg_prime, tvg_prime, tvs_prime
95 real :: ths_prime, thg_prime, thvs_prime, thvg_prime
96
97 real, parameter :: ka = 2.4E-5
98
99 integer :: iregime
100
101 rcp = gas_constant/cp
102
103 ! 1 Compute the roughness length based upon season and land use
104 ! =====================================
105
106 ! 1.1 Define the rouhness length
107 ! -----------------
108
109 z0 = roughness
110
111 if (z0 < 0.0001) z0 = 0.0001
112
113 ! 1.2 Define the rouhgness length for moisture
114 ! -----------------
115
116 if (xland .ge. 1.5) then
117 zq0 = z0
118 else
119 zq0 = zint0
120 end if
121
122 ! 1.3 Define the some constant variable for psi
123 ! -----------------
124
125 gzsoz0 = log(hs/z0)
126
127 gz10oz0 = log(h10/z0)
128
129 gz2oz0 = log(h2/z0)
130
131
132 ! 2. Calculate the virtual temperature
133 ! =====================================
134
135 ! 2.1 Compute Virtual temperature on the lowest half sigma level
136 ! ---------------------------------------------------------
137
138 tvs_prime = ts_prime * (1. + 0.608 * qs) + 0.608 * ts * qs_prime
139 tvs = ts * (1. + 0.608 * qs)
140
141 ! 2.2 Compute the ground saturated mixing ratio and the ground virtual
142 ! temperature
143 ! ----------------------------------------------------------------
144
145 call da_tp_to_qs(tg, psfc, eg, qg)
146 call da_tp_to_qs_lin1(tg, psfc, eg, tg_prime, psfc_prime, qg_prime)
147
148 qg_prime = qg_prime * qg
149
150 tvg_prime = tg_prime * (1. + 0.608 * qg) + 0.608 * tg * qg_prime
151 tvg = tg * (1. + 0.608 * qg)
152
153 ! 3. Compute the potential temperature and virtual potential temperature
154 ! =======================================================================
155
156 ! 3.1 Potential temperature on the lowest half sigma level
157 ! ----------------------------------------------------
158
159 Pi = (100000. / ps) ** rcp
160 ths_prime = (ts_prime - ps_prime * rcp * ts/ps) * Pi
161 ths = ts * Pi
162
163 ! 3.2 Virtual potential temperature on the lowest half sigma level
164 ! ------------------------------------------------------------
165
166 thvs_prime = (tvs_prime - ps_prime * rcp * tvs/ps) * Pi
167 thvs = tvs * Pi
168
169 ! 3.3 Potential temperature at the ground
170 ! -----------------------------------
171
172 Pi = (100000. / psfc) ** rcp
173 thg_prime = (tg_prime - psfc_prime * rcp * tg/psfc) * Pi
174 thg = tg * Pi
175
176 ! 3.4 Virtual potential temperature at ground
177 ! ---------------------------------------
178
179 thvg_prime = (tvg_prime - psfc_prime * rcp * tvg/psfc) * Pi
180 thvg = tvg * Pi
181
182 ! 4. BULK RICHARDSON NUMBER AND MONI-OBUKOV LENGTH
183 ! =================================================
184
185 ! 4.1 Velocity
186 ! --------
187
188 ! Wind speed:
189
190 Va2_prime = 2.*us*us_prime + 2.*vs*vs_prime
191 Va2 = us*us + vs*vs
192
193 ! Convective velocity:
194
195 if (thvg >= thvs) then
196 Vc2_prime = 4. * (thvg_prime - thvs_prime)
197 Vc2 = 4. * (thvg - thvs)
198 else
199 Vc2_prime = 0.
200 Vc2 = 0.
201 end if
202
203 V2_prime = Va2_prime+ Vc2_prime
204 V2 = Va2 + Vc2
205
206 ! 4.2 Bulk richardson number
207 ! ----------------------
208
209 Pi = gravity * hs / (ths*V2)
210 rib_prime = (thvs_prime - thvg_prime &
211 - (thvs-thvg)/V2 * V2_prime &
212 - (thvs-thvg)/ths * ths_prime) * Pi
213 rib = (thvs - thvg) * Pi
214
215 ! 5. CALCULATE PSI BASED UPON REGIME
216 ! =======================================
217
218 iregime = int(regime)
219
220 select case (iregime)
221
222 ! 5.1 Stable conditions (REGIME 1)
223 ! ---------------------------
224
225 case (1);
226
227 psim_prime = 0.
228 psimz_prime = 0.
229 psim2_prime = 0.
230 psim = -10.*gzsoz0
231 psimz = -10.*gz10oz0
232 psim2 = -10.*gz2oz0
233 psim = max(psim,-10.)
234 psimz = max(psimz,-10.)
235 psim2 = max(psim2,-10.)
236
237 psih_prime = psim_prime
238 psihz_prime = psimz_prime
239 psih2_prime = psim2_prime
240 psih = psim
241 psihz = psimz
242 psih2 = psim2
243
244 ! 5.2 Mechanically driven turbulence (REGIME 2)
245 ! ------------------------------------------
246
247 case (2);
248
249 Pi = - 1.0 / ((1.1 - 5.*rib)*(1.1 - 5.*rib))
250 psim_prime = 5.5 * gzsoz0 * rib_prime * Pi
251 psimz_prime = 5.5 * gz10oz0 * rib_prime * Pi
252 psim2_prime = 5.5 * gz2oz0 * rib_prime * Pi
253
254 Pi = (-5. * rib) / (1.1 - 5.*rib)
255 psim = gzsoz0 * Pi
256 psimz = gz10oz0 * Pi
257 psim2 = gz2oz0 * Pi
258
259 if (psim >= -10.) then
260 psim = psim
261 psim_prime = psim_prime
262 else
263 psim = -10.
264 psim_prime = 0.
265 end if
266 if (psimz >= -10.) then
267 psimz = psimz
268 psimz_prime = psimz_prime
269 else
270 psimz = -10.
271 psimz_prime = 0.
272 end if
273 if (psim2 >= -10.) then
274 psim2 = psim2
275 psim2_prime = psim2_prime
276 else
277 psim2 = -10.
278 psim2_prime = 0.
279 end if
280
281 psih_prime = psim_prime
282 psihz_prime = psimz_prime
283 psih2_prime = psim2_prime
284 psih = psim
285 psihz = psimz
286 psih2 = psim2
287
288 ! 5.3 Unstable Forced convection (REGIME 3)
289 ! -------------------------------------
290
291 case (3);
292
293 psim_prime = 0.
294 psimz_prime = 0.
295 psim2_prime = 0.
296
297 psim = 0.
298 psimz = 0.
299 psim2 = 0.
300
301 psih_prime = psim_prime
302 psihz_prime = psimz_prime
303 psih2_prime = psim2_prime
304 psih = psim
305 psihz = psimz
306 psih2 = psim2
307
308
309 ! 5.4 Free convection (REGIME 4)
310 ! --------------------------
311
312 case (4);
313
314 ! Calculate psi m and pshi h using iteration method
315
316 psim_prime = 0.
317 psih_prime = 0.
318 psim = 0.
319 psih = 0.
320 cc = 2. * atan(1.0)
321
322 ! do k = 1 , k_iteration
323
324 ! 5.4.1 Calculate ust, m/L (mol), h/L (hol)
325 ! --------------------------
326
327 ! Friction speed
328
329 ust = k_kar * sqrt(v2) /(gzsoz0 - psim)
330 ust_prime = (0.5/V2 * v2_prime + psim_prime /(gzsoz0 - psim)) * ust
331
332 ! Heat fux factor
333
334 mol = k_kar * (ths - thg)/(gzsoz0 - psih)
335 mol_prime = ((ths_prime - thg_prime) /(ths - thg) + &
336 psih_prime /(gzsoz0 - psih)) * mol
337
338 ! Ratio of PBL height to Monin-Obukhov length
339
340 if (ust .LT. 0.01) then
341 hol_prime = rib_prime * gzsoz0
342 hol = rib * gzsoz0
343 else
344 hol = k_kar * gravity * hs * mol / (ths * ust * ust)
345 hol_prime = (mol_prime / mol - ths_prime / ths &
346 - 2.* ust_prime / ust) * hol
347 end if
348
349 ! 5.4.2 Calculate n, nz, R, Rz
350 ! --------------------------
351
352 if (hol >= 0.) then
353 hol_prime = 0.
354 hol = 0.
355 else
356 hol_prime = hol_prime
357 hol = hol
358 end if
359 if (hol >= -10.) then
360 hol_prime = hol_prime
361 hol = hol
362 else
363 hol_prime = 0.
364 hol = -10.
365 end if
366
367 holz_prime = (h10 / hs) * hol_prime
368 holz = (h10 / hs) * hol
369 if (holz >= 0.) then
370 holz_prime = 0.
371 holz = 0.
372 else
373 holz_prime = holz_prime
374 holz = holz
375 end if
376 if (holz >= -10.) then
377 holz_prime = holz_prime
378 holz = holz
379 else
380 holz_prime = 0.
381 holz = -10.
382 end if
383
384 hol2_prime = (h2 / hs) * hol_prime
385 hol2 = (h2 / hs) * hol
386 if (hol2 >= 0.) then
387 hol2_prime = 0.
388 hol2 = 0.
389 else
390 hol2_prime = hol2_prime
391 hol2 = hol2
392 end if
393 if (hol2 >= -10.) then
394 hol2_prime = hol2_prime
395 hol2 = hol2
396 else
397 hol2_prime = 0.
398 hol2 = -10.
399 end if
400
401 ! 5.4.3 Calculate Psim & psih
402 ! --------------------------
403
404 ! Using the continuous function:
405 xx_prime = -4.* hol_prime /((1. - 16. * hol) ** 0.75)
406 xx = (1. - 16. * hol) ** 0.25
407 yy_prime = 2.* xx * xx_prime /(1.+xx*xx)
408 yy = log((1.+xx*xx)/2.)
409 psim_prime = 2 * xx_prime *(1.0/(1.+xx)-1.0/(1+xx*xx)) + yy_prime
410 psim = 2. * log((1.+xx)/2.) + yy - 2. * atan(xx) + cc
411 psih_prime = 2. * yy_prime
412 psih = 2. * yy
413
414 ! Using the continuous function:
415 xx_prime = -4.* holz_prime /((1. - 16. * holz) ** 0.75)
416 xx = (1. - 16. * holz) ** 0.25
417 yy_prime = 2.* xx * xx_prime /(1.+xx*xx)
418 yy = log((1.+xx*xx)/2.)
419 psimz_prime = 2.* xx_prime *(1.0/(1.+xx)-1.0/(1+xx*xx)) + yy_prime
420 psimz = 2. * log((1.+xx)/2.) + yy - 2. * atan(xx) + cc
421 psihz_prime = 2. * yy_prime
422 psihz = 2. * yy
423
424 ! Using the continuous function:
425 xx_prime = -4.* hol2_prime /((1. - 16. * hol2) ** 0.75)
426 xx = (1. - 16. * hol2) ** 0.25
427 yy_prime = 2.* xx * xx_prime /(1.+xx*xx)
428 yy = log((1.+xx*xx)/2.)
429 psim2_prime = 2.* xx_prime *(1.0/(1.+xx)-1.0/(1+xx*xx)) + yy_prime
430 psim2 = 2. * log((1.+xx)/2.) + yy - 2. * atan(xx) + cc
431 psih2_prime = 2. * yy_prime
432 psih2 = 2. * yy
433
434 ! end do
435
436 ! 5.4.4 Define the limit value for psim & psih
437 ! --------------------------
438
439 if (psim <= 0.9*gzsoz0) then
440 psim_prime = psim_prime
441 psim = psim
442 else
443 psim = 0.9*gzsoz0
444 psim_prime = 0.
445 end if
446 if (psimz <= 0.9*gz10oz0) then
447 psimz_prime = psimz_prime
448 psimz = psimz
449 else
450 psimz_prime = 0.
451 psimz = 0.9*gz10oz0
452 end if
453 if (psim2 <= 0.9*gz2oz0) then
454 psim2_prime = psim2_prime
455 psim2 = psim2
456 else
457 psim2_prime = 0.
458 psim2 = 0.9*gz2oz0
459 end if
460 if (psih <= 0.9*gzsoz0) then
461 psih_prime = psih_prime
462 psih = psih
463 else
464 psih_prime = 0.
465 psih = 0.9*gzsoz0
466 end if
467 if (psihz <= 0.9*gz10oz0) then
468 psihz_prime = psihz_prime
469 psihz = psihz
470 else
471 psihz_prime = 0.
472 psihz = 0.9*gz10oz0
473 end if
474 if (psih2 <= 0.9*gz2oz0) then
475 psih2_prime = psih2_prime
476 psih2 = psih2
477 else
478 psih2_prime = 0.
479 psih2 = 0.9*gz2oz0
480 end if
481
482 case default;
483 write(unit=message(1),fmt='(A,I2,A)') &
484 "Regime=",iregime," is invalid."
485 call da_error(__FILE__,__LINE__,message(1:1))
486
487 end select
488
489 ! 6. CALCULATE PSI FOR WinD, TEMPERATURE AND MOISTURE
490 ! =======================================
491
492 psiw_prime = - psim_prime
493 psiw = gzsoz0 - psim
494 psiz_prime = - psimz_prime
495 psiz = gz10oz0 - psimz
496 psit_prime = - psih_prime
497 psit = gzsoz0 - psih
498 psit2_prime = - psih2_prime
499 psit2 = gz2oz0 - psih2
500
501 ust = k_kar * sqrt(v2) /(gzsoz0 - psim)
502 ust_prime = (0.5/V2 * v2_prime + psim_prime /(gzsoz0 - psim)) * ust
503
504 psiq2_prime = k_kar*hs/(ka*(k_kar*ust*hs/ka + hs / zq0))*ust_prime
505 psiq_prime = psiq2_prime - psih_prime
506 psiq2_prime = psiq2_prime - psih2_prime
507
508 psiq = log(k_kar*ust*hs/ka + hs / zq0) - psih
509 psiq2 = log(k_kar*ust*h2/ka + h2 / zq0) - psih2
510
511 ! 7. CALCULATE THE PERTURBATIONS for 10M WinD, 2M TEMPERATURE AND MOISTURE
512 ! =======================================
513
514 Pi = psiz / psiw
515 u10_prime= (us_prime + us/psiz * psiz_prime - us/psiw * psiw_prime) * Pi
516 v10_prime= (vs_prime + vs/psiz * psiz_prime - vs/psiw * psiw_prime) * Pi
517
518 t2_prime = ((1.0-psit2/psit) * thg_prime + (ths_prime + &
519 (ths - thg)/psit2 * psit2_prime - &
520 (ths - thg)/psit * psit_prime) * psit2/psit &
521 + rcp*(thg + (ths - thg)*psit2/psit)/psfc * psfc_prime) &
522 * (psfc/100000.)**rcp
523
524 q2_prime = (1.0-psiq2/psiq) * qg_prime + psiq2/psiq * qs_prime + &
525 (qs -qg)*(psiq2/psiq) * (psiq2_prime/psiq2 - psiq_prime/psiq)
526
527 end subroutine da_sfc_wtq_lin
528
529