da_tpq_to_slp_adj.inc
References to this file elsewhere.
1 subroutine da_tpq_to_slp_adj (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
17 implicit none
18
19 real, intent(in) :: TERR, PSFC
20 real, intent(in) :: T(kms:kme), Q(kms:kme), P(kms:kme)
21 real, intent(out) :: T9(kms:kme), Q9(kms:kme), P9(kms:kme)
22 real, intent(in) :: SLP9
23 real, intent(out) :: PSFC9
24
25 ! real :: SLP
26
27 integer :: K, KLO, KHI
28 real :: PL, T0, TS, XTERM,TLO, THI, TL
29 real :: PL9,T09,TS9,TLO9,THI9,TL9,COEF1,COEF2
30 real :: T08,TS8
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_adj")
38
39 ! sea level pressure
40
41 xterm=gamma* gas_constant / gravity
42
43 ! compute pressure at pconst mb above surface (pl)
44
45 T9(:) = 0.0
46 Q9(:) = 0.0
47 P9(:) = 0.0
48
49 if (terr <= 0.0) then
50 psfc9=slp9
51 if (trace_use) call da_trace_exit("da_tpq_to_slp_adj")
52 return
53 end if
54
55 PL = psfc - PCONST
56 klo = 0
57
58 ! find 2 levels on sigma surfaces surrounding pl at each i,j
59
60 do k=kts, kte-1
61 if ((p(k) >= pl) .and. (p(k+1) < pl)) then
62 khi = k+1
63 klo = k
64 exit
65 end if
66 end do
67
68 if (klo < 1) then
69 write(unit=message(1),fmt='(A,F11.3,A)') &
70 'Error finding pressure level ',PCONST,' Mb above the surface'
71 write(unit=message(2),fmt='(A,F11.3,2X,A,F11.3)') &
72 'PL=',PL,' PSFC=',psfc
73 call da_error(__FILE__,__LINE__,message(1:2))
74 end if
75
76 ! get temperature at pl (tl), extrapolate t at surface (ts)
77 ! and t at sea level (t0) with 6.5 k/km lapse rate
78
79 TLO = t(KLO) * (EPS+q(KLO))/(EPS*(1.0+q(KLO)))
80 THI = t(KHI) * (EPS+q(KHI))/(EPS*(1.0+q(KHI)))
81 COEF1 = ALOG(PL/p(KHI))
82 COEF2 = ALOG(p(KLO)/p(KHI))
83 TL = THI-(THI-TLO)*COEF1/COEF2
84 TS = TL*(psfc/PL)**XTERM
85 TS8 = TS
86 T0 = TS +GAMMA*terr
87 T08 = T0
88
89 ! Correct sea level temperature if too hot
90
91 if (t0 >= tc) then
92 if (ts <= tc) then
93 t0 = tc
94 else
95 t0 = tc-0.005*(ts-tc)**2
96 end if
97 end if
98
99 psfc9 = slp9*EXP(2.0*gravity*TERR/(gas_constant*(TS+T0)))
100 TS9 = -slp9*psfc*EXP(2.0*gravity*TERR/(gas_constant*(TS+T0)))* &
101 2.0*gravity*TERR/(gas_constant*(TS+T0)**2)
102 T09 = -slp9*psfc*EXP(2.0*gravity*TERR/(gas_constant*(TS+T0)))* &
103 2.0*gravity*TERR/(gas_constant*(TS+T0)**2)
104
105 if ( t08 >= tc ) then
106 if ( ts8 > tc ) then
107 ts9=ts9-0.01*(ts-tc)*t09
108 end if
109 t09=0.0
110 end if
111
112 TS9 = TS9+T09
113 TL9 = TS9*(psfc/PL)**XTERM
114 psfc9 = psfc9+TS9*XTERM*(TL/PL)*(psfc/PL)**(XTERM-1)
115 PL9 = -TS9*XTERM*(TL*psfc/(PL*PL))*(psfc/PL)**(XTERM-1)
116 THI9 = (1.0-COEF1/COEF2)*TL9
117 TLO9 = COEF1/COEF2*TL9
118 PL9 = PL9-(THI-TLO)/(COEF2*PL)*TL9
119 p9(KHI) = p9(KHI)+((THI-TLO)/(COEF2*p(KHI))*(1.0-COEF1/COEF2))*TL9
120 p9(KLO) = p9(KLO)+(THI-TLO)*COEF1/(COEF2*COEF2*p(KLO))*TL9
121 t9(KHI) = t9(KHI)+THI9* (EPS+q(KHI))/(EPS*(1.0+q(KHI)))
122 q9(KHI) = q9(KHI)+THI9*t(KHI)*(1.0-EPS)/(EPS*(1.0+q(KHI))**2)
123 t9(KLO) = t9(KLO)+TLO9* (EPS+q(KLO))/(EPS*(1.0+q(KLO)))
124 q9(KLO) = q9(KLO)+TLO9*t(KLO)*(1.0-EPS)/(EPS*(1.0+q(KLO))**2)
125
126 psfc9=psfc9+PL9
127
128 if (trace_use) call da_trace_exit("da_tpq_to_slp_adj")
129
130 end subroutine da_tpq_to_slp_adj
131
132