da_sfc_wtq_adj.inc
References to this file elsewhere.
1 subroutine da_sfc_wtq_adj (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) :: hs, roughness, xland
50 real , intent (in) :: u10_prime, v10_prime, t2_prime, q2_prime
51
52 real , intent (inout) :: ps_prime, ts_prime, qs_prime , &
53 us_prime, vs_prime, psfc_prime, tg_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 real :: ust_s, hol_s, psim_s, psim2_s, psimz_s, psih_s, psih2_s, psihz_s
86
87 real :: Vc2_prime, Va2_prime, V2_prime
88 real :: rib_prime, xx_prime, yy_prime
89 real :: psiw_prime, psiz_prime, mol_prime, ust_prime, &
90 hol_prime, holz_prime, hol2_prime
91 real :: psim_prime, psimz_prime, psim2_prime, &
92 psih_prime, psihz_prime, psih2_prime
93 real :: psit_prime, psit2_prime, &
94 psiq_prime, psiq2_prime
95 real :: qg_prime, tvg_prime, tvs_prime
96 real :: ths_prime, thg_prime, thvs_prime, thvg_prime
97
98 real, parameter :: ka = 2.4E-5
99
100 integer :: iregime
101
102 !-----------------------------------------------------------------------------!
103 ! initialize perturbation value
104 !------- ----------------------------------------------------------------------!
105 ! tg_prime = 0.0
106 ! us_prime = 0.0
107 ! vs_prime = 0.0
108 ! ts_prime = 0.0
109 ! ps_prime = 0.0
110 ! qs_prime = 0.0
111 ! psfc_prime = 0.0
112
113 psim_prime = 0.0
114 psimz_prime = 0.0
115 psim2_prime = 0.0
116 psih2_prime = 0.0
117 psihz_prime = 0.0
118 psih_prime = 0.0
119
120 psiw_prime = 0.0
121 psiz_prime = 0.0
122 psit_prime = 0.0
123 psit2_prime = 0.0
124 psiq_prime = 0.0
125 psiq2_prime = 0.0
126
127 qg_prime = 0.0
128 ths_prime = 0.0
129 thg_prime = 0.0
130
131 thvs_prime = 0.0
132 thvg_prime = 0.0
133 V2_prime = 0.0
134 rib_prime = 0.0
135 ust_prime = 0.0
136 tvg_prime = 0.0
137 tvs_prime = 0.0
138 va2_prime = 0.0
139 vc2_prime = 0.0
140
141 ! +++++++++++++++++++++++++++++++++
142 ! A. Calculate Basic state
143 ! +++++++++++++++++++++++++++++++++
144 rcp = gas_constant/cp
145
146 ! 1 Compute the roughness length based upon season and land use
147 ! =====================================
148
149 ! 1.1 Define the rouhness length
150 ! -----------------
151
152 z0 = roughness
153
154 if (z0 < 0.0001) z0 = 0.0001
155
156 ! 1.2 Define the rouhgness length for moisture
157 ! -----------------
158
159 if (xland .ge. 1.5) then
160 zq0 = z0
161 else
162 zq0 = zint0
163 end if
164
165 ! 1.3 Define the some constant variable for psi
166 ! -----------------
167
168 gzsoz0 = log(hs/z0)
169
170 gz10oz0 = log(h10/z0)
171
172 gz2oz0 = log(h2/z0)
173
174
175 ! 2. Calculate the virtual temperature
176 ! =====================================
177
178 ! 2.1 Compute Virtual temperature on the lowest half sigma level
179 ! ---------------------------------------------------------
180
181 tvs = ts * (1. + 0.608 * qs)
182
183 ! 2.2 Convert ground virtual temperature assuming it's saturated
184 ! ----------------------------------------------------------
185 call da_tp_to_qs(tg, psfc, eg, qg)
186 tvg = tg * (1. + 0.608 * qg)
187
188 ! 3. Compute the potential temperature and virtual potential temperature
189 ! =======================================================================
190
191 ! 3.1 Potential temperature on the lowest half sigma level
192 ! ----------------------------------------------------
193
194 Pi = (100000. / ps) ** rcp
195 ths = ts * Pi
196
197 ! 3.2 Virtual potential temperature on the lowest half sigma level
198 ! ------------------------------------------------------------
199
200 thvs = tvs * Pi
201
202 ! 3.3 Potential temperature at the ground
203 ! -----------------------------------
204
205 Pi = (100000. / psfc) ** rcp
206 thg = tg * Pi
207
208 ! 3.4 Virtual potential temperature at ground
209 ! ---------------------------------------
210
211 thvg = tvg * Pi
212
213
214 ! 4. BULK RICHARDSON NUMBER AND MONI-OBUKOV LENGTH
215 ! =================================================
216
217 ! 4.1 Velocity
218 ! --------
219
220 ! Wind speed:
221
222 Va2 = us*us + vs*vs
223
224 ! Convective velocity:
225
226 if (thvg >= thvs) then
227 Vc2 = 4. * (thvg - thvs)
228 else
229 Vc2 = 0.
230 end if
231
232 V2 = Va2 + Vc2
233
234 ! 4.2 Bulk richardson number
235 ! ----------------------
236
237 rib = (gravity * hs / ths) * (thvs - thvg) / V2
238 ! 5. CALCULATE PSI BASED UPON REGIME
239 ! =======================================
240
241 iregime = int(regime)
242 select case (iregime)
243
244 ! 5.1 Stable conditions (REGIME 1)
245 ! ---------------------------
246
247 case (1);
248
249 psim = -10.*gzsoz0
250 psimz = -10.*gz10oz0
251 psim2 = -10.*gz2oz0
252
253 psim_s = psim
254 psimz_s = psimz
255 psim2_s = psim2
256
257 psim = max(psim,-10.)
258 psimz = max(psimz,-10.)
259 psim2 = max(psim2,-10.)
260
261 psih = psim
262 psihz = psimz
263 psih2 = psim2
264
265 ! 5.2 Mechanically driven turbulence (REGIME 2)
266 ! ------------------------------------------
267
268 case (2);
269 Pi = (-5. * rib) / (1.1 - 5.*rib)
270 psim = gzsoz0 * Pi
271 psimz = gz10oz0 * Pi
272 psim2 = gz2oz0 * Pi
273 psim_s = psim
274 psimz_s = psimz
275 psim2_s = psim2
276 if (psim >= -10.) then
277 psim = psim
278 else
279 psim = -10.
280 end if
281 if (psimz >= -10.) then
282 psimz = psimz
283 else
284 psimz = -10.
285 end if
286 if (psim2 >= -10.) then
287 psim2 = psim2
288 else
289 psim2 = -10.
290 end if
291
292 psih = psim
293 psihz = psimz
294 psih2 = psim2
295
296 ! 5.3 Unstable Forced convection (REGIME 3)
297 ! -------------------------------------
298
299 case (3);
300
301 psim = 0.
302 psimz = 0.
303 psim2 = 0.
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 = 0.
317 psih = 0.
318 cc = 2. * atan(1.0)
319 !
320 ! do k = 1 , k_iteration
321
322 ! 5.4.1 Calculate ust, m/L (mol), h/L (hol)
323 ! --------------------------
324
325 ! Friction speed
326
327 ust = k_kar * sqrt(v2) /(gzsoz0 - psim)
328
329 ! save ust for adjoint:
330 ust_s = ust
331
332 ! Heat flux factor
333
334 mol = k_kar * (ths - thg)/(gzsoz0 - psih)
335
336 ! Ratio of PBL height to Monin-Obukhov length
337
338 if (ust .LT. 0.01) then
339 hol = rib * gzsoz0
340 else
341 hol = k_kar * gravity * hs * mol / (ths * ust * ust)
342 end if
343
344 ! 5.4.2 Calculate n, nz, R, Rz
345 ! --------------------------
346
347 hol_s = hol
348
349 if (hol >= 0.) then
350 hol = 0.
351 else
352 hol = hol
353 end if
354 if (hol >= -10.) then
355 hol = hol
356 else
357 hol = -10.
358 end if
359
360 holz = (h10 / hs) * hol
361 if (holz >= 0.) then
362 holz = 0.
363 else
364 holz = holz
365 end if
366 if (holz >= -10.) then
367 holz = holz
368 else
369 holz = -10.
370 end if
371
372 hol2 = (h2 / hs) * hol
373 if (hol2 >= 0.) then
374 hol2 = 0.
375 else
376 hol2 = hol2
377 end if
378 if (hol2 >= -10.) then
379 hol2 = hol2
380 else
381 hol2 = -10.
382 end if
383
384 ! 5.4.3 Calculate Psim & psih
385 ! --------------------------
386
387 ! Using the continuous function:
388 xx = (1. - 16. * hol) ** 0.25
389 yy = log((1.+xx*xx)/2.)
390 psim = 2. * log((1.+xx)/2.) + yy - 2. * atan(xx) + cc
391 psih = 2. * yy
392
393 ! Using the continuous function:
394 xx = (1. - 16. * holz) ** 0.25
395 yy = log((1.+xx*xx)/2.)
396 psimz = 2. * log((1.+xx)/2.) + yy - 2. * atan(xx) + cc
397 psihz = 2. * yy
398
399 ! Using the continuous function:
400 xx = (1. - 16. * hol2) ** 0.25
401 yy = log((1.+xx*xx)/2.)
402 psim2 = 2. * log((1.+xx)/2.) + yy - 2. * atan(xx) + cc
403 psih2 = 2. * yy
404
405 ! end do
406
407 ! 5.4.4 Define the limit value for psim & psih
408 ! --------------------------
409
410 ! Save the values for adjoint:
411
412 psim_s = psim
413 psimz_s = psimz
414 psim2_s = psim2
415 psih_s = psih
416 psihz_s = psihz
417 psih2_s = psih2
418
419 if (psim <= 0.9*gzsoz0) then
420 psim = psim
421 else
422 psim = 0.9*gzsoz0
423 end if
424 if (psimz <= 0.9*gz10oz0) then
425 psimz = psimz
426 else
427 psimz = 0.9*gz10oz0
428 end if
429 if (psim2 <= 0.9*gz2oz0) then
430 psim2 = psim2
431 else
432 psim2 = 0.9*gz2oz0
433 end if
434 if (psih <= 0.9*gzsoz0) then
435 psih = psih
436 else
437 psih = 0.9*gzsoz0
438 end if
439 if (psihz <= 0.9*gz10oz0) then
440 psihz = psihz
441 else
442 psihz = 0.9*gz10oz0
443 end if
444 if (psih2 <= 0.9*gz2oz0) then
445 psih2 = psih2
446 else
447 psih2 = 0.9*gz2oz0
448 end if
449
450 case default;
451 write(unit=message(1),fmt='(A,I2,A)') &
452 "Regime=",iregime," is invalid."
453 call da_error(__FILE__,__LINE__,message(1:1))
454
455 end select
456
457 ! 6. CALCULATE PSI FOR WinD, TEMPERATURE AND MOISTURE
458 ! =======================================
459
460 psiw = gzsoz0 - psim
461 psiz = gz10oz0 - psimz
462 psit = gzsoz0 - psih
463 psit2 = gz2oz0 - psih2
464 ust = k_kar * sqrt(v2) /(gzsoz0 - psim)
465 psiq = log(k_kar*ust*hs/ka + hs / zq0) - psih
466 psiq2 = log(k_kar*ust*h2/ka + h2 / zq0) - psih2
467 ! +++++++++++++++++++++++++++++++++
468 ! B. Calculate adjoint solution
469 ! +++++++++++++++++++++++++++++++++
470
471 ! 7. CALCULATE 10M WinD, 2M TEMPERATURE AND MOISTURE
472 ! =======================================
473
474 qg_prime = (1.0 - psiq2/psiq) * q2_prime
475 qs_prime = qs_prime + psiq2/psiq * q2_prime
476 psiq2_prime = (qs -qg)/psiq * q2_prime
477 psiq_prime = -(qs -qg)*psiq2/(psiq*psiq) * q2_prime
478 ! q2_prime = 0.
479
480 Pi = (psfc/100000.)**rcp
481 thg_prime = (1.0 - psit2/psit) * Pi * t2_prime
482 ths_prime = (psit2/psit) * Pi * t2_prime
483 psit2_prime = (ths - thg)/psit * Pi * t2_prime
484 psit_prime = - (ths - thg)*psit2/(psit*psit) * Pi * t2_prime
485 psfc_prime = psfc_prime + Pi * rcp*(thg + (ths - thg)*psit2/psit)/psfc * t2_prime
486 ! t2_prime = 0.
487
488 Pi = psiz / psiw
489 vs_prime = vs_prime + Pi * v10_prime
490 psiz_prime = vs / psiz * Pi * v10_prime
491 psiw_prime = - vs / psiw * Pi * v10_prime
492 ! v10_prime = 0.
493
494 us_prime = us_prime + Pi * u10_prime
495 psiz_prime = psiz_prime + us / psiz * Pi * u10_prime
496 psiw_prime = psiw_prime - us / psiw * Pi * u10_prime
497 ! u10_prime = 0.
498
499 ! 6. CALCULATE PSI FOR WinD, TEMPERATURE AND MOISTURE
500 ! =======================================
501
502 ! moisture part:
503 psih2_prime = - psiq2_prime
504 psih_prime = - psiq_prime
505 psiq2_prime = psiq2_prime + psiq_prime
506 ust_prime = k_kar*hs/(ka*(k_kar*ust*hs/ka + hs / zq0)) * psiq2_prime
507
508 v2_prime = 0.5/V2 * ust * ust_prime
509 psim_prime = ust / (gzsoz0 - psim) * ust_prime
510 ust_prime = 0.0
511
512 ! temperature part:
513 psih2_prime = psih2_prime - psit2_prime
514 psih_prime = psih_prime - psit_prime
515
516 ! wind part:
517 psimz_prime = psimz_prime - psiz_prime
518 psim_prime = psim_prime - psiw_prime
519
520 ! 5. CALCULATE PSI BASED UPON REGIME
521 ! =======================================
522 select case (iregime)
523
524 ! 5.1 Stable conditions (REGIME 1)
525 ! ---------------------------
526
527 case (1);
528
529 psim2_prime = psim2_prime + psih2_prime
530 psimz_prime = psimz_prime + psihz_prime
531 psim_prime = psim_prime + psih_prime
532 psim_prime = 0.
533 psimz_prime = 0.
534 psim2_prime = 0.
535
536 ! 5.2 Mechanically driven turbulence (REGIME 2)
537 ! ------------------------------------------
538
539 case (2);
540
541 psim2_prime = psim2_prime + psih2_prime
542 psimz_prime = psimz_prime + psihz_prime
543 psim_prime = psim_prime + psih_prime
544
545 psim = psim_s
546 psimz = psimz_s
547 psim2 = psim2_S
548
549 if (psim2 >= -10.) then
550 psim2_prime = psim2_prime
551 else
552 psim2_prime = 0.
553 end if
554 if (psimz >= -10.) then
555 psimz_prime = psimz_prime
556 else
557 psimz_prime = 0.
558 end if
559 if (psim >= -10.) then
560 psim_prime = psim_prime
561 else
562 psim_prime = 0.
563 end if
564
565 Pi = -1.0 / ((1.1 - 5.*rib)*(1.1 - 5.*rib))
566 rib_prime = 5.5 * gz2oz0 * psim2_prime * Pi
567 rib_prime = rib_prime + 5.5 * gz10oz0 * psimz_prime * Pi
568 rib_prime = rib_prime + 5.5 * gzsoz0 * psim_prime * Pi
569
570 ! 5.3 Unstable Forced convection (REGIME 3)
571 ! -------------------------------------
572
573 case (3);
574
575 psim2_prime = psih2_prime
576 psimz_prime = psihz_prime
577 psim_prime = psih_prime
578 psim2_prime = 0.
579 psimz_prime = 0.
580 psim_prime = 0.
581
582 ! 5.4 Free convection (REGIME 4)
583 ! --------------------------
584
585 case (4);
586
587 ! 5.4.4 Define the limit value for psim & psih
588 ! -------------------------------------
589
590 ! Recover the values:
591
592 psim = psim_s
593 psim2 = psim2_s
594 psimz = psimz_s
595 psihz = psihz_s
596 psih = psih_s
597 psih2 = psih2_s
598
599 if (psih2 <= 0.9*gz2oz0) then
600 psih2_prime = psih2_prime
601 else
602 psih2_prime = 0.
603 end if
604 if (psihz <= 0.9*gz10oz0) then
605 psihz_prime = psihz_prime
606 else
607 psihz_prime = 0.
608 end if
609 if (psih <= 0.9*gzsoz0) then
610 psih_prime = psih_prime
611 else
612 psih_prime = 0.
613 end if
614 if (psim2 <= 0.9*gz2oz0) then
615 psim2_prime = psim2_prime
616 else
617 psim2_prime = 0.
618 end if
619 if (psimz <= 0.9*gz10oz0) then
620 psimz_prime = psimz_prime
621 else
622 psimz_prime = 0.
623 end if
624 if (psim <= 0.9*gzsoz0) then
625 psim_prime = psim_prime
626 else
627 psim_prime = 0.
628 end if
629
630 ! 5.4.3 Calculate Psim & psih
631 ! --------------------------
632
633 ! Recover ust:
634 ust = ust_s
635 psim = 0.
636 psih = 0.
637
638 xx = (1. - 16. * hol2) ** 0.25
639 yy = log((1.+xx*xx)/2.)
640 yy_prime = 2. * psih2_prime
641 psih2_prime = 0.
642
643 xx_prime = 2.*(1.0/(1.+xx)-1.0/(1+xx*xx)) * psim2_prime
644 yy_prime = yy_prime + psim2_prime
645 psim2_prime = 0.
646
647 xx_prime = xx_prime + 2.* xx/(1.+xx*xx) * yy_prime
648 yy_prime = 0.
649 hol2_prime = -4./((1. - 16. * hol2) ** 0.75) * xx_prime
650 xx_prime = 0.
651
652 xx = (1. - 16. * holz) ** 0.25
653 yy = log((1.+xx*xx)/2.)
654 yy_prime = 2. * psihz_prime
655 psihz_prime = 0.
656
657 xx_prime = 2.*(1.0/(1.+xx)-1.0/(1+xx*xx)) * psimz_prime
658 yy_prime = yy_prime + psimz_prime
659 psimz_prime = 0.
660
661 xx_prime = xx_prime + 2.* xx/(1.+xx*xx) * yy_prime
662 yy_prime = 0.
663 holz_prime = -4./((1. - 16. * holz) ** 0.75) * xx_prime
664 xx_prime = 0.
665
666 xx = (1. - 16. * hol) ** 0.25
667 yy = log((1.+xx*xx)/2.)
668 yy_prime = 2. * psih_prime
669 psih_prime = 0.
670
671 xx_prime = 2.*(1.0/(1.+xx)-1.0/(1+xx*xx)) * psim_prime
672 yy_prime = yy_prime + psim_prime
673 psim_prime = 0.
674
675 xx_prime = xx_prime + 2.* xx/(1.+xx*xx)*yy_prime
676 yy_prime = 0.
677 hol_prime = -4./((1. - 16. * hol) ** 0.75)*xx_prime
678 xx_prime = 0.
679
680 ! 5.4.2 Calculate n, nz, R, Rz
681 ! --------------------------
682
683 ! Recover the values:
684
685 hol = hol_s
686
687 hol2 = (h2 / hs) * hol
688 if (hol2 >= -10.) then
689 hol2_prime = hol2_prime
690 else
691 hol2_prime = 0.
692 end if
693 if (hol2 >= 0.) then
694 hol2_prime = 0.
695 else
696 hol2_prime = hol2_prime
697 end if
698
699 hol_prime = hol_prime + (h2 / hs) * hol2_prime
700 hol2_prime = 0.
701
702 holz = (h10 / hs) * hol
703 if (holz >= -10.) then
704 holz_prime = holz_prime
705 else
706 holz_prime = 0.
707 end if
708 if (holz >= 0.) then
709 holz_prime = 0.
710 else
711 holz_prime = holz_prime
712 end if
713
714 hol_prime = hol_prime + (h10 / hs) * holz_prime
715 holz_prime = 0.
716
717 if (hol >= -10.) then
718 hol_prime = hol_prime
719 else
720 hol_prime = 0.
721 end if
722 if (hol >= 0.) then
723 hol_prime = 0.
724 else
725 hol_prime = hol_prime
726 end if
727
728 ! 5.4.1 Calculate ust, m/L (mol), h/L (hol)
729 ! --------------------------
730
731 ! Ratio of PBL height to Monin-Obukhov length
732 if (ust .LT. 0.01) then
733 rib_prime = hol_prime * gzsoz0
734 hol_prime = 0.
735 else
736 mol_prime = hol / mol * hol_prime
737 ths_prime = ths_prime - hol / ths * hol_prime
738 ust_prime = - 2.0 * hol / ust * hol_prime
739 hol_prime = 0.
740 end if
741
742 ! Heat flux factor
743 ths_prime = ths_prime + mol/(ths - thg) * mol_prime
744 thg_prime = thg_prime - mol/(ths - thg) * mol_prime
745 psih_prime = psih_prime + mol/(gzsoz0 - psih) * mol_prime
746 mol_prime = 0.
747
748 ! Friction speed
749
750 v2_prime = V2_prime + 0.5 * ust/V2 * ust_prime
751 psim_prime = psim_prime + ust/(gzsoz0 - psim) * ust_prime
752 ust_prime = 0.
753
754 ! Calculate psi m and pshi h using iteration method
755
756 psim_prime = 0.
757 psih_prime = 0.
758
759 case default;
760 write(unit=message(1),fmt='(A,I2,A)') &
761 "Regime=",iregime," is invalid."
762 call da_error(__FILE__,__LINE__,message(1:1))
763
764 end select
765
766 ! 4. BULK RICHARDSON NUMBER AND MONI-OBUKOV LENGTH
767 ! =================================================
768
769 ! 4.2 Bulk richardson number
770 ! ----------------------
771
772 Pi = gravity * hs / (ths*V2)
773 ths_prime = ths_prime - Pi * (thvs-thvg)/ths * rib_prime
774 V2_prime = V2_prime - Pi * (thvs-thvg)/V2 * rib_prime
775 thvs_prime = thvs_prime + Pi * rib_prime
776 thvg_prime = thvg_prime - Pi * rib_prime
777 rib_prime = 0.
778
779 ! 4.1 Velocity
780 ! --------
781
782 ! Convective velocity:
783
784 Va2_prime = V2_prime
785 Vc2_prime = V2_prime
786 V2_prime = 0.
787
788 if (thvg >= thvs) then
789 thvg_prime = thvg_prime + 4. * Vc2_prime
790 thvs_prime = thvs_prime - 4. * Vc2_prime
791 Vc2_prime = 0.
792 else
793 Vc2_prime = 0.
794 end if
795
796 ! Wind speed:
797
798 us_prime = us_prime + 2. *us * Va2_prime
799 vs_prime = vs_prime + 2. *vs * Va2_prime
800 Va2_prime = 0.
801
802 ! 3. Virtual potential temperature
803 ! ==================================
804
805 ! 3.4 Virtual potential temperature at ground
806 ! ---------------------------------------
807
808 Pi = (100000. / psfc) ** rcp
809 tvg_prime = tvg_prime + Pi * thvg_prime
810 psfc_prime = psfc_prime - thvg_prime * rcp * tvg/psfc * Pi
811 thvg_prime = 0.
812
813 ! 3.3 Potential temperature at the ground
814 ! -----------------------------------
815
816 tg_prime = tg_prime + Pi * thg_prime
817 psfc_prime = psfc_prime - thg_prime *rcp * tg/psfc * Pi
818 thg_prime = 0.
819
820 ! 3.2 Virtual potential temperature on the lowest half sigma level
821 ! ------------------------------------------------------------
822
823 Pi = (100000. / ps) ** rcp
824 tvs_prime = tvs_prime + PI * thvs_prime
825 ps_prime = ps_prime - thvs_prime *rcp * tvs/ps * Pi
826 thvs_prime = 0.
827
828 ! 3.1 Potential temperature on the lowest half sigma level
829 ! ----------------------------------------------------
830 ts_prime = ts_prime + Pi * ths_prime
831 ps_prime = ps_prime - ths_prime * rcp * ts/ps * Pi
832 ths_prime = 0.
833
834 ! 2. Calculate the virtual temperature
835 ! =====================================
836
837 ! 2.2 Compute the ground saturated mixing ratio and the ground virtual
838 ! temperature
839 ! ----------------------------------------------------------------
840
841 tg_prime = tg_prime + tvg_prime * (1. + 0.608 * qg)
842 qg_prime = qg_prime + tvg_prime * 0.608 * tg
843 tvg_prime = 0.
844
845 qg_prime = qg_prime * qg
846 call da_tp_to_qs_adj1(tg, psfc, eg, tg_prime, psfc_prime, qg_prime)
847
848 ! 2.1 Compute Virtual temperature on the lowest half sigma level
849 ! ---------------------------------------------------------
850
851 ts_prime = ts_prime + tvs_prime * (1. + 0.608 * qs)
852 qs_prime = qs_prime + tvs_prime * 0.608 * ts
853 tvs_prime = 0.
854
855 end subroutine da_sfc_wtq_adj
856
857