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, us_prime, vs_prime, psfc_prime, tg_prime
50 real, intent (in) :: hs, roughness, xland
51 real, intent (out) :: u10_prime, v10_prime, t2_prime, q2_prime
52
53 ! Maximum number of iterations in computing psim, psih
54
55 integer, parameter :: k_iteration = 10
56 ! integer, parameter :: k_iteration = 1
57
58 ! h10 is the height of 10m where the wind observed
59 ! h2 is the height of 2m where the temperature and
60 ! moisture observed.
61
62 real, parameter :: h10 = 10.0, h2 = 2.0
63 !
64 ! Default roughness over the land
65
66 real, parameter :: zint0 = 0.01
67 !
68 ! Von Karman constant
69
70 real, parameter :: k_kar = 0.4
71 !
72 ! Working variables
73
74 real :: Vc2, Va2, V2
75 real :: rib, rcp, xx, yy, cc, Pi
76 real :: psiw, psiz, mol, ust, hol, holz, hol2
77 real :: psim, psimz, psim2, psih, psihz, psih2
78 real :: psit, psit2, psiq, psiq2
79 real :: gzsoz0, gz10oz0, gz2oz0
80 real :: eg, qg, tvg, tvs
81 real :: ths, thg, thvs, thvg
82 real :: zq0, z0
83
84 real :: Vc2_prime, Va2_prime, V2_prime
85 real :: rib_prime, xx_prime, yy_prime
86 real :: psiw_prime, psiz_prime, mol_prime, ust_prime, &
87 hol_prime, holz_prime, hol2_prime
88 real :: psim_prime, psimz_prime, psim2_prime, &
89 psih_prime, psihz_prime, psih2_prime
90 real :: psit_prime, psit2_prime, &
91 psiq_prime, psiq2_prime
92 real :: qg_prime, tvg_prime, tvs_prime
93 real :: ths_prime, thg_prime, thvs_prime, thvg_prime
94
95 real, parameter :: ka = 2.4E-5
96
97 integer :: iregime
98
99 if (trace_use) call da_trace_entry("da_sfc_wtq_lin")
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 + 0.608 * qs) + 0.608 * ts * qs_prime
139 tvs = ts * (1.0 + 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 + 0.608 * qg) + 0.608 * tg * qg_prime
151 tvg = tg * (1.0 + 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.0 / 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.0 / 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.0*us*us_prime + 2.0*vs*vs_prime
191 Va2 = us*us + vs*vs
192
193 ! Convective velocity:
194
195 if (thvg >= thvs) then
196 Vc2_prime = 4.0 * (thvg_prime - thvs_prime)
197 Vc2 = 4.0 * (thvg - thvs)
198 else
199 Vc2_prime = 0.0
200 Vc2 = 0.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.0
228 psimz_prime = 0.0
229 psim2_prime = 0.0
230 psim = -10.0*gzsoz0
231 psimz = -10.0*gz10oz0
232 psim2 = -10.0*gz2oz0
233 psim = max(psim,-10.0)
234 psimz = max(psimz,-10.0)
235 psim2 = max(psim2,-10.0)
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.0*rib)*(1.1 - 5.0*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.0 * rib) / (1.1 - 5.0*rib)
255 psim = gzsoz0 * Pi
256 psimz = gz10oz0 * Pi
257 psim2 = gz2oz0 * Pi
258
259 if (psim >= -10.0) then
260 psim = psim
261 psim_prime = psim_prime
262 else
263 psim = -10.0
264 psim_prime = 0.0
265 end if
266 if (psimz >= -10.0) then
267 psimz = psimz
268 psimz_prime = psimz_prime
269 else
270 psimz = -10.0
271 psimz_prime = 0.0
272 end if
273 if (psim2 >= -10.0) then
274 psim2 = psim2
275 psim2_prime = psim2_prime
276 else
277 psim2 = -10.0
278 psim2_prime = 0.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.0
294 psimz_prime = 0.0
295 psim2_prime = 0.0
296
297 psim = 0.0
298 psimz = 0.0
299 psim2 = 0.0
300
301 psih_prime = psim_prime
302 psihz_prime = psimz_prime
303 psih2_prime = psim2_prime
304
305 psih = psim
306 psihz = psimz
307 psih2 = psim2
308
309
310 ! 5.4 Free convection (REGIME 4)
311 ! --------------------------
312
313 case (4);
314
315 ! Calculate psi m and pshi h using iteration method
316
317 psim_prime = 0.0
318 psih_prime = 0.0
319 psim = 0.0
320 psih = 0.0
321 cc = 2.0 * atan(1.0)
322
323 ! do k = 1 , k_iteration
324
325 ! 5.4.1 Calculate ust, m/L (mol), h/L (hol)
326 ! --------------------------
327
328 ! Friction speed
329
330 ust = k_kar * sqrt(v2) /(gzsoz0 - psim)
331 ust_prime = (0.5/V2 * v2_prime + psim_prime /(gzsoz0 - psim)) * ust
332
333 ! Heat fux factor
334
335 mol = k_kar * (ths - thg)/(gzsoz0 - psih)
336 mol_prime = ((ths_prime - thg_prime) /(ths - thg) + &
337 psih_prime /(gzsoz0 - psih)) * mol
338
339 ! Ratio of PBL height to Monin-Obukhov length
340
341 if (ust .LT. 0.01) then
342 hol_prime = rib_prime * gzsoz0
343 hol = rib * gzsoz0
344 else
345 hol = k_kar * gravity * hs * mol / (ths * ust * ust)
346 hol_prime = (mol_prime / mol - ths_prime / ths &
347 - 2.0* ust_prime / ust) * hol
348 end if
349
350 ! 5.4.2 Calculate n, nz, R, Rz
351 ! --------------------------
352
353 if (hol >= 0.0) then
354 hol_prime = 0.0
355 hol = 0.0
356 else
357 hol_prime = hol_prime
358 hol = hol
359 end if
360 if (hol >= -10.0) then
361 hol_prime = hol_prime
362 hol = hol
363 else
364 hol_prime = 0.0
365 hol = -10.0
366 end if
367
368 holz_prime = (h10 / hs) * hol_prime
369 holz = (h10 / hs) * hol
370 if (holz >= 0.0) then
371 holz_prime = 0.0
372 holz = 0.0
373 else
374 holz_prime = holz_prime
375 holz = holz
376 end if
377 if (holz >= -10.0) then
378 holz_prime = holz_prime
379 holz = holz
380 else
381 holz_prime = 0.0
382 holz = -10.0
383 end if
384
385 hol2_prime = (h2 / hs) * hol_prime
386 hol2 = (h2 / hs) * hol
387 if (hol2 >= 0.0) then
388 hol2_prime = 0.0
389 hol2 = 0.0
390 else
391 hol2_prime = hol2_prime
392 hol2 = hol2
393 end if
394 if (hol2 >= -10.0) then
395 hol2_prime = hol2_prime
396 hol2 = hol2
397 else
398 hol2_prime = 0.0
399 hol2 = -10.0
400 end if
401
402 ! 5.4.3 Calculate Psim & psih
403 ! --------------------------
404
405 ! Using the continuous function:
406 xx_prime = -4.0* hol_prime /((1.0 - 16.0 * hol) ** 0.75)
407 xx = (1.0 - 16.0 * hol) ** 0.25
408 yy_prime = 2.0* xx * xx_prime /(1.0+xx*xx)
409 yy = log((1.0+xx*xx)/2.0)
410 psim_prime = 2 * xx_prime *(1.0/(1.0+xx)-1.0/(1+xx*xx)) + yy_prime
411 psim = 2.0 * log((1.0+xx)/2.0) + yy - 2.0 * atan(xx) + cc
412 psih_prime = 2.0 * yy_prime
413 psih = 2.0 * yy
414
415 ! Using the continuous function:
416 xx_prime = -4.0* holz_prime /((1.0 - 16.0 * holz) ** 0.75)
417 xx = (1.0 - 16.0 * holz) ** 0.25
418 yy_prime = 2.0* xx * xx_prime /(1.0+xx*xx)
419 yy = log((1.0+xx*xx)/2.0)
420 psimz_prime = 2.0* xx_prime *(1.0/(1.0+xx)-1.0/(1+xx*xx)) + yy_prime
421 psimz = 2.0 * log((1.0+xx)/2.0) + yy - 2.0 * atan(xx) + cc
422 psihz_prime = 2.0 * yy_prime
423 psihz = 2.0 * yy
424
425 ! Using the continuous function:
426 xx_prime = -4.0* hol2_prime /((1.0 - 16.0 * hol2) ** 0.75)
427 xx = (1.0 - 16.0 * hol2) ** 0.25
428 yy_prime = 2.0* xx * xx_prime /(1.0+xx*xx)
429 yy = log((1.0+xx*xx)/2.0)
430 psim2_prime = 2.0* xx_prime *(1.0/(1.0+xx)-1.0/(1+xx*xx)) + yy_prime
431 psim2 = 2.0 * log((1.0+xx)/2.0) + yy - 2.0 * atan(xx) + cc
432 psih2_prime = 2.0 * yy_prime
433 psih2 = 2.0 * yy
434
435 ! end do
436
437 ! 5.4.4 Define the limit value for psim & psih
438 ! --------------------------
439
440 if (psim <= 0.9*gzsoz0) then
441 psim_prime = psim_prime
442 psim = psim
443 else
444 psim = 0.9*gzsoz0
445 psim_prime = 0.0
446 end if
447 if (psimz <= 0.9*gz10oz0) then
448 psimz_prime = psimz_prime
449 psimz = psimz
450 else
451 psimz_prime = 0.0
452 psimz = 0.9*gz10oz0
453 end if
454 if (psim2 <= 0.9*gz2oz0) then
455 psim2_prime = psim2_prime
456 psim2 = psim2
457 else
458 psim2_prime = 0.0
459 psim2 = 0.9*gz2oz0
460 end if
461 if (psih <= 0.9*gzsoz0) then
462 psih_prime = psih_prime
463 psih = psih
464 else
465 psih_prime = 0.0
466 psih = 0.9*gzsoz0
467 end if
468 if (psihz <= 0.9*gz10oz0) then
469 psihz_prime = psihz_prime
470 psihz = psihz
471 else
472 psihz_prime = 0.0
473 psihz = 0.9*gz10oz0
474 end if
475 if (psih2 <= 0.9*gz2oz0) then
476 psih2_prime = psih2_prime
477 psih2 = psih2
478 else
479 psih2_prime = 0.0
480 psih2 = 0.9*gz2oz0
481 end if
482
483 case default;
484 write(unit=message(1),fmt='(A,I2,A)') &
485 "Regime=",iregime," is invalid."
486 call da_error(__FILE__,__LINE__,message(1:1))
487
488 end select
489
490 ! 6. CALCULATE PSI FOR WinD, TEMPERATURE AND MOISTURE
491 ! =======================================
492
493 psiw_prime = - psim_prime
494 psiw = gzsoz0 - psim
495 psiz_prime = - psimz_prime
496 psiz = gz10oz0 - psimz
497 psit_prime = - psih_prime
498 psit = gzsoz0 - psih
499 psit2_prime = - psih2_prime
500 psit2 = gz2oz0 - psih2
501
502 ust = k_kar * sqrt(v2) /(gzsoz0 - psim)
503 ust_prime = (0.5/V2 * v2_prime + psim_prime /(gzsoz0 - psim)) * ust
504
505 psiq2_prime = k_kar*hs/(ka*(k_kar*ust*hs/ka + hs / zq0))*ust_prime
506 psiq_prime = psiq2_prime - psih_prime
507 psiq2_prime = psiq2_prime - psih2_prime
508
509 psiq = log(k_kar*ust*hs/ka + hs / zq0) - psih
510 psiq2 = log(k_kar*ust*h2/ka + h2 / zq0) - psih2
511
512 ! 7. CALCULATE THE PERTURBATIONS for 10M WinD, 2M TEMPERATURE AND MOISTURE
513 ! =======================================
514
515 Pi = psiz / psiw
516 u10_prime= (us_prime + us/psiz * psiz_prime - us/psiw * psiw_prime) * Pi
517 v10_prime= (vs_prime + vs/psiz * psiz_prime - vs/psiw * psiw_prime) * Pi
518
519 t2_prime = ((1.0-psit2/psit) * thg_prime + (ths_prime + &
520 (ths - thg)/psit2 * psit2_prime - &
521 (ths - thg)/psit * psit_prime) * psit2/psit &
522 + rcp*(thg + (ths - thg)*psit2/psit)/psfc * psfc_prime) &
523 * (psfc/100000.0)**rcp
524
525 q2_prime = (1.0-psiq2/psiq) * qg_prime + psiq2/psiq * qs_prime + &
526 (qs -qg)*(psiq2/psiq) * (psiq2_prime/psiq2 - psiq_prime/psiq)
527
528 if (trace_use) call da_trace_exit("da_sfc_wtq_lin")
529
530 end subroutine da_sfc_wtq_lin
531
532