da_tpq_to_slp_lin.inc
References to this file elsewhere.
1 subroutine da_tpq_to_slp_lin ( T, Q, P, TERR, PSFC, T9, Q9, P9, PSFC9, SLP9)
2
3 !-----------------------------------------------------------------------
4 ! Purpose: computes sea level pressure from the rule
5 ! t1/t2=(p1/p2)**(gamma*r/g).
6 !
7 ! input t temperature
8 ! q mixing ratio
9 ! p pressure
10 ! terr terrain
11 ! psfc surface pressure
12 !
13 ! output slp sea level pressure
14 !-----------------------------------------------------------------------
15
16 implicit none
17
18 real, intent(in) :: terr, psfc, psfc9
19 real, intent(in) :: t(kms:kme)
20 real, intent(in) :: q(kms:kme)
21 real, intent(in) :: p(kms:kme)
22 real, intent(in) :: t9(kms:kme)
23 real, intent(in) :: q9(kms:kme)
24 real, intent(in) :: p9(kms:kme)
25 ! real :: slp
26 real, intent(out) :: slp9
27
28 integer :: k, klo, khi
29 real :: pl, t0, ts, xterm,tlo, thi, tl
30 real :: pl9,t09,ts9,tlo9,thi9,tl9,coef1,coef2
31
32 real, parameter :: gamma = 6.5e-3
33 real, parameter :: tc = t_kelvin+17.5
34 real, parameter :: pconst = 10000.0
35 real, parameter :: eps = 0.622
36
37 if (trace_use) call da_trace_entry("da_tpq_to_slp_lin")
38
39 ! ... sea level pressure
40
41 xterm=gamma* gas_constant / gravity
42
43 ! compute pressure at pconst mb above surface (pl)
44
45 if (terr <= 0.0) then
46 slp9 = psfc9
47 ! slp = psfc
48 if (trace_use) call da_trace_exit("da_tpq_to_slp_lin")
49 return
50 end if
51
52
53 pl9 = psfc9
54 pl = psfc - pconst
55 klo = 0
56
57 ! find 2 levels on sigma surfaces surrounding pl at each i,j
58
59 k_loop: do k=kts, kte-1
60 if ((p(k) >= pl) .and. (p(k+1) < pl)) then
61 khi = k+1
62 klo = k
63 exit k_loop
64 end if
65 end do k_loop
66
67 if (klo < 1) then
68 write(unit=message(1),fmt='(A,F11.3,A)') &
69 'error finding pressure level ',pconst,' mb above the surface'
70 write(unit=message(2),fmt='(A,F11.3,2X,A,F11.3)') 'PL=',PL,' PSFC=',psfc
71 call da_error(__FILE__,__LINE__,message(1:2))
72 end if
73
74 ! get temperature at pl (tl), extrapolate t at surface (ts)
75 ! and t at sea level (t0) with 6.5 k/km lapse rate
76
77 tlo9=t9(klo) * (eps+q(klo))/(eps*(1.0+q(klo))) + &
78 q9(klo)*t(klo)*(1.0-eps)/(eps*(1.0+q(klo))**2)
79 tlo=t(klo) * (eps+q(klo))/(eps*(1.0+q(klo)))
80 thi9=t9(khi) * (eps+q(khi))/(eps*(1.0+q(khi)))+ &
81 q9(khi)*t(khi)*(1.0-eps)/(eps*(1.0+q(khi))**2)
82 thi=t(khi) * (eps+q(khi))/(eps*(1.0+q(khi)))
83 coef1=alog(pl/p(khi))
84 coef2=alog(p(klo)/p(khi))
85 tl9=(1.0-coef1/coef2)*thi9+coef1/coef2*tlo9 &
86 -(thi-tlo)/(coef2*pl)*pl9 &
87 +((thi-tlo)/(coef2*p(khi))*(1.0-coef1/coef2))*p9(khi) &
88 +(thi-tlo)*coef1/(coef2*coef2*p(klo))*p9(klo)
89 tl=thi-(thi-tlo)*coef1/coef2
90 ts9=tl9*(psfc/pl)**xterm+psfc9*xterm*(tl/pl)*(psfc/pl)** &
91 (xterm-1)-pl9*xterm*(tl*psfc/(pl*pl))*(psfc/pl)**(xterm-1)
92 ts=tl*(psfc/pl)**xterm
93 t09=ts9
94 t0=ts +gamma*terr
95
96 ! correct sea level temperature if too hot
97
98 if ( t0 >= tc ) then
99 if ( ts <= tc ) then
100 t09 = 0.0
101 t0 = tc
102 else
103 t09 = -0.01*(ts-tc)*ts9
104 t0 = tc-0.005*(ts-tc)**2
105 end if
106 end if
107
108 ! compute sea level pressure
109
110 slp9=psfc9*exp(2.0*gravity*terr/(gas_constant*(ts+t0))) &
111 -psfc*exp(2.0*gravity*terr/(gas_constant*(ts+t0)))* &
112 2.0*gravity*terr/(gas_constant*(ts+t0)**2)*(ts9+t09)
113 ! slp=psfc*exp(2.0*gravity*terr/(gas_constant*(ts+t0)))
114
115 if (trace_use) call da_trace_exit("da_tpq_to_slp_lin")
116
117 end subroutine da_tpq_to_slp_lin
118
119