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 real, intent (inout) :: ps_prime, ts_prime, qs_prime, us_prime, vs_prime, psfc_prime, tg_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 real :: ust_s, hol_s, psim_s, psim2_s, psimz_s, psih_s, psih2_s, psihz_s
84
85 real :: Vc2_prime, Va2_prime, V2_prime
86 real :: rib_prime, xx_prime, yy_prime
87 real :: psiw_prime, psiz_prime, mol_prime, ust_prime, &
88 hol_prime, holz_prime, hol2_prime
89 real :: psim_prime, psimz_prime, psim2_prime, &
90 psih_prime, psihz_prime, psih2_prime
91 real :: psit_prime, psit2_prime, &
92 psiq_prime, psiq2_prime
93 real :: qg_prime, tvg_prime, tvs_prime
94 real :: ths_prime, thg_prime, thvs_prime, thvg_prime
95
96 real, parameter :: ka = 2.4E-5
97
98 integer :: iregime
99
100 if (trace_use) call da_trace_entry("da_sfc_wtq_adj")
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.0 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.0 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 + 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 + 0.608 * qg)
187
188 ! 3.0 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.0 / 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.0 / 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.0 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.0 * (thvg - thvs)
228 else
229 Vc2 = 0.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.0 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.0*gzsoz0
250 psimz = -10.0*gz10oz0
251 psim2 = -10.0*gz2oz0
252
253 psim_s = psim
254 psimz_s = psimz
255 psim2_s = psim2
256
257 psim = max(psim,-10.0)
258 psimz = max(psimz,-10.0)
259 psim2 = max(psim2,-10.0)
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.0 * rib) / (1.1 - 5.0*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.0) then
277 psim = psim
278 else
279 psim = -10.0
280 end if
281 if (psimz >= -10.0) then
282 psimz = psimz
283 else
284 psimz = -10.0
285 end if
286 if (psim2 >= -10.0) then
287 psim2 = psim2
288 else
289 psim2 = -10.0
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.0
302 psimz = 0.0
303 psim2 = 0.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.0
317 psih = 0.0
318 cc = 2.0 * 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.0) then
350 hol = 0.0
351 else
352 hol = hol
353 end if
354 if (hol >= -10.0) then
355 hol = hol
356 else
357 hol = -10.0
358 end if
359
360 holz = (h10 / hs) * hol
361 if (holz >= 0.0) then
362 holz = 0.0
363 else
364 holz = holz
365 end if
366 if (holz >= -10.0) then
367 holz = holz
368 else
369 holz = -10.0
370 end if
371
372 hol2 = (h2 / hs) * hol
373 if (hol2 >= 0.0) then
374 hol2 = 0.0
375 else
376 hol2 = hol2
377 end if
378 if (hol2 >= -10.0) then
379 hol2 = hol2
380 else
381 hol2 = -10.0
382 end if
383
384 ! 5.4.3 Calculate Psim & psih
385 ! --------------------------
386
387 ! Using the continuous function:
388 xx = (1.0 - 16.0 * hol) ** 0.25
389 yy = log((1.0+xx*xx)/2.0)
390 psim = 2.0 * log((1.0+xx)/2.0) + yy - 2.0 * atan(xx) + cc
391 psih = 2.0 * yy
392
393 ! Using the continuous function:
394 xx = (1.0 - 16.0 * holz) ** 0.25
395 yy = log((1.0+xx*xx)/2.0)
396 psimz = 2.0 * log((1.0+xx)/2.0) + yy - 2.0 * atan(xx) + cc
397 psihz = 2.0 * yy
398
399 ! Using the continuous function:
400 xx = (1.0 - 16.0 * hol2) ** 0.25
401 yy = log((1.0+xx*xx)/2.0)
402 psim2 = 2.0 * log((1.0+xx)/2.0) + yy - 2.0 * atan(xx) + cc
403 psih2 = 2.0 * 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.0 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.0 Calculate adjoint solution
469 ! +++++++++++++++++++++++++++++++++
470
471 ! 7.0 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.0
479
480 Pi = (psfc/100000.0)**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.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.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.0
498
499 ! 6.0 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.0 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.0
533 psimz_prime = 0.0
534 psim2_prime = 0.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.0) then
550 psim2_prime = psim2_prime
551 else
552 psim2_prime = 0.0
553 end if
554 if (psimz >= -10.0) then
555 psimz_prime = psimz_prime
556 else
557 psimz_prime = 0.0
558 end if
559 if (psim >= -10.0) then
560 psim_prime = psim_prime
561 else
562 psim_prime = 0.0
563 end if
564
565 Pi = -1.0 / ((1.1 - 5.*rib)*(1.1 - 5.0*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
579 psim2_prime = 0.0
580 psimz_prime = 0.0
581 psim_prime = 0.0
582
583 ! 5.4 Free convection (REGIME 4)
584 ! --------------------------
585
586 case (4);
587
588 ! 5.4.4 Define the limit value for psim & psih
589 ! -------------------------------------
590
591 ! Recover the values:
592
593 psim = psim_s
594 psim2 = psim2_s
595 psimz = psimz_s
596 psihz = psihz_s
597 psih = psih_s
598 psih2 = psih2_s
599
600 if (psih2 <= 0.9*gz2oz0) then
601 psih2_prime = psih2_prime
602 else
603 psih2_prime = 0.0
604 end if
605 if (psihz <= 0.9*gz10oz0) then
606 psihz_prime = psihz_prime
607 else
608 psihz_prime = 0.0
609 end if
610 if (psih <= 0.9*gzsoz0) then
611 psih_prime = psih_prime
612 else
613 psih_prime = 0.0
614 end if
615 if (psim2 <= 0.9*gz2oz0) then
616 psim2_prime = psim2_prime
617 else
618 psim2_prime = 0.0
619 end if
620 if (psimz <= 0.9*gz10oz0) then
621 psimz_prime = psimz_prime
622 else
623 psimz_prime = 0.0
624 end if
625 if (psim <= 0.9*gzsoz0) then
626 psim_prime = psim_prime
627 else
628 psim_prime = 0.0
629 end if
630
631 ! 5.4.3 Calculate Psim & psih
632 ! --------------------------
633
634 ! Recover ust:
635 ust = ust_s
636 psim = 0.0
637 psih = 0.0
638
639 xx = (1.0 - 16.0 * hol2) ** 0.25
640 yy = log((1.0+xx*xx)/2.0)
641 yy_prime = 2.0 * psih2_prime
642 psih2_prime = 0.0
643
644 xx_prime = 2.0*(1.0/(1.0+xx)-1.0/(1+xx*xx)) * psim2_prime
645 yy_prime = yy_prime + psim2_prime
646 psim2_prime = 0.0
647
648 xx_prime = xx_prime + 2.0* xx/(1.0+xx*xx) * yy_prime
649 yy_prime = 0.0
650 hol2_prime = -4.0/((1.0 - 16.0 * hol2) ** 0.75) * xx_prime
651 xx_prime = 0.0
652
653 xx = (1.0 - 16.0 * holz) ** 0.25
654 yy = log((1.0+xx*xx)/2.0)
655 yy_prime = 2.0 * psihz_prime
656 psihz_prime = 0.0
657
658 xx_prime = 2.0*(1.0/(1.0+xx)-1.0/(1+xx*xx)) * psimz_prime
659 yy_prime = yy_prime + psimz_prime
660 psimz_prime = 0.0
661
662 xx_prime = xx_prime + 2.0* xx/(1.0+xx*xx) * yy_prime
663 yy_prime = 0.0
664 holz_prime = -4.0/((1.0 - 16.0 * holz) ** 0.75) * xx_prime
665 xx_prime = 0.0
666
667 xx = (1.0 - 16.0 * hol) ** 0.25
668 yy = log((1.0+xx*xx)/2.0)
669 yy_prime = 2.0 * psih_prime
670 psih_prime = 0.0
671
672 xx_prime = 2.0*(1.0/(1.0+xx)-1.0/(1+xx*xx)) * psim_prime
673 yy_prime = yy_prime + psim_prime
674 psim_prime = 0.0
675
676 xx_prime = xx_prime + 2.0* xx/(1.0+xx*xx)*yy_prime
677 yy_prime = 0.0
678 hol_prime = -4.0/((1.0 - 16.0 * hol) ** 0.75)*xx_prime
679 xx_prime = 0.0
680
681 ! 5.4.2 Calculate n, nz, R, Rz
682 ! --------------------------
683
684 ! Recover the values:
685
686 hol = hol_s
687
688 hol2 = (h2 / hs) * hol
689 if (hol2 >= -10.0) then
690 hol2_prime = hol2_prime
691 else
692 hol2_prime = 0.0
693 end if
694 if (hol2 >= 0.0) then
695 hol2_prime = 0.0
696 else
697 hol2_prime = hol2_prime
698 end if
699
700 hol_prime = hol_prime + (h2 / hs) * hol2_prime
701 hol2_prime = 0.0
702
703 holz = (h10 / hs) * hol
704 if (holz >= -10.0) then
705 holz_prime = holz_prime
706 else
707 holz_prime = 0.0
708 end if
709 if (holz >= 0.0) then
710 holz_prime = 0.0
711 else
712 holz_prime = holz_prime
713 end if
714
715 hol_prime = hol_prime + (h10 / hs) * holz_prime
716 holz_prime = 0.0
717
718 if (hol >= -10.0) then
719 hol_prime = hol_prime
720 else
721 hol_prime = 0.0
722 end if
723 if (hol >= 0.0) then
724 hol_prime = 0.0
725 else
726 hol_prime = hol_prime
727 end if
728
729 ! 5.4.1 Calculate ust, m/L (mol), h/L (hol)
730 ! --------------------------
731
732 ! Ratio of PBL height to Monin-Obukhov length
733 if (ust .LT. 0.01) then
734 rib_prime = hol_prime * gzsoz0
735 hol_prime = 0.0
736 else
737 mol_prime = hol / mol * hol_prime
738 ths_prime = ths_prime - hol / ths * hol_prime
739 ust_prime = - 2.0 * hol / ust * hol_prime
740 hol_prime = 0.0
741 end if
742
743 ! Heat flux factor
744 ths_prime = ths_prime + mol/(ths - thg) * mol_prime
745 thg_prime = thg_prime - mol/(ths - thg) * mol_prime
746 psih_prime = psih_prime + mol/(gzsoz0 - psih) * mol_prime
747 mol_prime = 0.0
748
749 ! Friction speed
750
751 v2_prime = V2_prime + 0.5 * ust/V2 * ust_prime
752 psim_prime = psim_prime + ust/(gzsoz0 - psim) * ust_prime
753 ust_prime = 0.0
754
755 ! Calculate psi m and pshi h using iteration method
756
757 psim_prime = 0.0
758 psih_prime = 0.0
759
760 case default;
761 write(unit=message(1),fmt='(A,I2,A)') &
762 "Regime=",iregime," is invalid."
763 call da_error(__FILE__,__LINE__,message(1:1))
764
765 end select
766
767 ! 4.0 BULK RICHARDSON NUMBER AND MONI-OBUKOV LENGTH
768 ! =================================================
769
770 ! 4.2 Bulk richardson number
771 ! ----------------------
772
773 Pi = gravity * hs / (ths*V2)
774 ths_prime = ths_prime - Pi * (thvs-thvg)/ths * rib_prime
775 V2_prime = V2_prime - Pi * (thvs-thvg)/V2 * rib_prime
776 thvs_prime = thvs_prime + Pi * rib_prime
777 thvg_prime = thvg_prime - Pi * rib_prime
778 rib_prime = 0.0
779
780 ! 4.1 Velocity
781 ! --------
782
783 ! Convective velocity:
784
785 Va2_prime = V2_prime
786 Vc2_prime = V2_prime
787 V2_prime = 0.0
788
789 if (thvg >= thvs) then
790 thvg_prime = thvg_prime + 4.0 * Vc2_prime
791 thvs_prime = thvs_prime - 4.0 * Vc2_prime
792 Vc2_prime = 0.0
793 else
794 Vc2_prime = 0.0
795 end if
796
797 ! Wind speed:
798
799 us_prime = us_prime + 2.0 *us * Va2_prime
800 vs_prime = vs_prime + 2.0 *vs * Va2_prime
801 Va2_prime = 0.0
802
803 ! 3.0 Virtual potential temperature
804 ! ==================================
805
806 ! 3.4 Virtual potential temperature at ground
807 ! ---------------------------------------
808
809 Pi = (100000.0 / psfc) ** rcp
810 tvg_prime = tvg_prime + Pi * thvg_prime
811 psfc_prime = psfc_prime - thvg_prime * rcp * tvg/psfc * Pi
812 thvg_prime = 0.0
813
814 ! 3.3 Potential temperature at the ground
815 ! -----------------------------------
816
817 tg_prime = tg_prime + Pi * thg_prime
818 psfc_prime = psfc_prime - thg_prime *rcp * tg/psfc * Pi
819 thg_prime = 0.0
820
821 ! 3.2 Virtual potential temperature on the lowest half sigma level
822 ! ------------------------------------------------------------
823
824 Pi = (100000.0 / ps) ** rcp
825 tvs_prime = tvs_prime + PI * thvs_prime
826 ps_prime = ps_prime - thvs_prime *rcp * tvs/ps * Pi
827 thvs_prime = 0.0
828
829 ! 3.1 Potential temperature on the lowest half sigma level
830 ! ----------------------------------------------------
831 ts_prime = ts_prime + Pi * ths_prime
832 ps_prime = ps_prime - ths_prime * rcp * ts/ps * Pi
833 ths_prime = 0.0
834
835 ! 2.0 Calculate the virtual temperature
836 ! =====================================
837
838 ! 2.2 Compute the ground saturated mixing ratio and the ground virtual
839 ! temperature
840 ! ----------------------------------------------------------------
841
842 tg_prime = tg_prime + tvg_prime * (1.0 + 0.608 * qg)
843 qg_prime = qg_prime + tvg_prime * 0.608 * tg
844 tvg_prime = 0.0
845
846 qg_prime = qg_prime * qg
847 call da_tp_to_qs_adj1(tg, psfc, eg, tg_prime, psfc_prime, qg_prime)
848
849 ! 2.1 Compute Virtual temperature on the lowest half sigma level
850 ! ---------------------------------------------------------
851
852 ts_prime = ts_prime + tvs_prime * (1.0 + 0.608 * qs)
853 qs_prime = qs_prime + tvs_prime * 0.608 * ts
854 tvs_prime = 0.0
855
856 if (trace_use) call da_trace_exit("da_sfc_wtq_adj")
857
858 end subroutine da_sfc_wtq_adj
859
860