da_wrf_tpq_2_slp.inc
References to this file elsewhere.
1 subroutine da_wrf_tpq_2_slp (grid)
2
3 !---------------------------------------------------------------------
4 ! Purpose: computes sea level pressure from the rule
5 ! t1/t2=(p1/p2)**(gamma*r/g).
6 !
7 ! input t temperature cross 3d
8 ! q mixing ratio cross 2d
9 ! p pressure perturbation cross 2d
10 ! terr terrain cross 2d
11 ! psfc surface pressure cross 2d
12 !
13 ! output slp sea level pressure cross 2d
14 !---------------------------------------------------------------------
15
16 implicit none
17
18
19 type (domain), intent(inout) :: grid
20
21 integer :: I, J, K, KLO, KHI
22 real :: PL, T0, TS, XTERM, &
23 TLO, THI, TL
24
25 real, parameter :: GAMMA = 6.5E-3, &
26 TC=t_kelvin+17.5, &
27 PCONST=10000.0 , &
28 EPS = 0.622
29
30 if (trace_use_dull) call da_trace_entry("da_wrf_tpq_2_slp")
31
32
33 ! sea level pressure
34
35 xterm=gamma* gas_constant / gravity
36
37 ! compute pressure at pconst mb above surface (pl)
38
39
40 j_loop: do j=jts, jte
41 i_loop: do i=its, ite
42 if (grid%xb%terr(i,j) <= 0.0) then
43 grid%xb%slp(i,j) = grid%xb%psfc(i,j)
44 cycle i_loop
45 end if
46
47 PL = grid%xb%psfc(i,j) - PCONST
48 klo = 0
49
50 ! FinD 2 LEVELS ON SIGMA SURFACES SURROUNDinG PL AT EACH I,J
51
52 k_loop: do k=kts, kte-1
53 if ((grid%xb%p(i,j,k) >= pl) .and. (grid%xb%p(i,j,k+1) < pl)) then
54 khi = k+1
55 klo = k
56 exit k_loop
57 end if
58 end do k_loop
59
60 if (klo < 1) then
61 write(unit=message(1),fmt='(A,F11.3,A)') &
62 'Cannot find pressure levelL ',PCONST,' mb above the surface'
63 write(unit=message(2),fmt='(A,F11.3,2X,A,F11.3)') &
64 'PL=',PL,' PSFC=',grid%xb%psfc(i,j)
65 call da_error(__FILE__,__LINE__,message(1:2))
66 end if
67
68 ! get temperature at pl (tl), extrapolate t at surface (ts)
69 ! and t at sea level (t0) with 6.5 k/km lapse rate
70
71 tlo=grid%xb%t(i,j,klo) * (eps+grid%xb%q(i,j,klo))/(eps*(1.0+grid%xb%q(i,j,klo)))
72 thi=grid%xb%t(i,j,khi) * (eps+grid%xb%q(i,j,khi))/(eps*(1.0+grid%xb%q(i,j,khi)))
73 tl=thi-(thi-tlo)*log(pl/grid%xb%p(i,j,khi)) &
74 /log(grid%xb%p(i,j,klo)/grid%xb%p(i,j,khi))
75 ts=tl*(grid%xb%psfc(i,j)/pl)**xterm
76 t0=ts +gamma*grid%xb%terr(i,j)
77
78 ! correct sea level temperature if too hot
79
80 if (t0 >= tc) then
81 if (ts <= tc) then
82 t0 = tc
83 else
84 t0 = tc-0.005*(ts-tc)**2
85 end if
86 end if
87
88 ! compute sea level pressure
89
90 grid%xb%slp(i,j)=grid%xb%psfc(i,j)*EXP(2.0*gravity*grid%xb%terr(i,j)/ &
91 (gas_constant*(TS+T0)))
92 end do i_loop
93 end do j_loop
94
95 if (trace_use_dull) call da_trace_exit("da_wrf_tpq_2_slp")
96
97 end subroutine da_wrf_tpq_2_slp
98
99