da_sfcprs.inc
References to this file elsewhere.
1 subroutine da_sfcprs (kts, kte, p, t, q, height, psfc, pslv, h_obs, t_obs)
2
3 ! Purpose: to compute the pressure at obs height (psfc) from the sea
4 ! level pressure (pslv), obs height (h_obs), temperature
5 ! (t_obs, optional), and the model profile (p,t,q).
6
7 implicit none
8
9 integer, intent(in) :: kts, kte
10 real, intent(in) :: height(kts:kte)
11 real, intent(in) :: p(kts:kte)
12 real, intent(in) :: t(kts:kte)
13 real, intent(in) :: q(kts:kte)
14 real, intent(in) :: h_obs, pslv
15 real, intent(out) :: psfc
16 real, optional, intent(in) :: t_obs
17
18 real, parameter :: g = gravity
19 real, parameter :: r = gas_constant
20 real, parameter :: rov2 = r / 2.0
21 real, parameter :: gamma = 6.5e-3 ! temperature lapse rate
22 real, parameter :: gammarg= gamma*gas_constant/g
23
24 real :: plmb(3)
25
26 integer :: k
27 integer :: k500
28 integer :: k700
29 integer :: k850
30
31 logical :: l1
32 logical :: l2
33 logical :: l3
34
35 real :: gamma78, gamma57
36 real :: ht
37 real :: p1
38 real :: t1
39 real :: t500
40 real :: t700
41 real :: t850
42 real :: tc
43 real :: tfixed
44 real :: tslv, tsfc
45
46 if (trace_use_dull) call da_trace_entry("da_sfcprs")
47
48 TC = t_kelvin + 17.5
49 PLMB = (/85000.0, 70000.0, 50000.0/)
50
51 ! Find the locations of the 850, 700 and 500 mb levels.
52
53 101 continue
54 K850 = 0 ! FinD K AT: P=850
55 K700 = 0 ! P=700
56 K500 = 0 ! P=500
57
58 do K = kte-1, kts, -1
59 if (p(k) > PLMB(1) .and. K850==0) then
60 K850 = K + 1
61 else if (p(k) > PLMB(2) .and. K700==0) then
62 K700 = K + 1
63 else if (p(k) > PLMB(3) .and. K500==0) then
64 K500 = K + 1
65 end if
66
67 end do
68
69 if ((K850 == 0) .OR. (K700 == 0) .OR. (K500 == 0)) then
70 ! write(unit=message(1),fmt='(A,3f8.0,A)') &
71 ! 'Error in finding p level for ',PLMB(1), PLMB(2), PLMB(3),' Pa.'
72 ! do K = 1, KX
73 ! write(unit=message(k+1),fmt='(A,I3,A,F10.2)') 'K = ',K,' PRESSURE = ',P(K)
74 ! end do
75 ! write(unit=message(kx+2),fmt='(A,2f8.0,A,f8.0,A)) ','Expected ', &
76 ! PLMB(1), PLMB(2),' and ',PLMB(3),' Pa. values, at least.'
77 ! message(kx+3)="not enough levels"
78 ! message(kx+4)="Change the pressure levels by -25mb"
79 ! call da_error(__FILE__,__LINE__,message(1:kx+4))
80 PLMB(1) = PLMB(1) - 2500.0
81 PLMB(2) = PLMB(2) - 2500.0
82 PLMB(3) = PLMB(3) - 2500.0
83 goto 101
84 end if
85
86 ! The 850 hPa level is called something special, and interpolated
87 ! to cross points. And then, we fill those last rows and columns.
88
89 ht = height(k850)
90
91 ! The variable ht is now -h_obs/ht(850 hPa). The plot thickens.
92
93 ht = -h_obs / ht
94
95 ! Make an isothermal assumption to get a first guess at the surface
96 ! pressure. This is to tell us which levels to use for the lapse
97 ! rates in a bit.
98
99 psfc = pslv * (pslv / p(k850)) ** ht
100
101 ! Get a pressure more than 100 hPa above the surface - P1. The
102 ! P1 is the top of the level that we will use for our lapse rate
103 ! computations.
104
105 if ((PSFC - p(k850) - 10000.0) >= 0.0) then
106 P1 = p(k850)
107 else if ((PSFC - p(k700)) >= 0.0) then
108 P1 = PSFC - 10000.0
109 else
110 P1 = p(k500)
111 end if
112
113 ! Compute virtual temperatures for k850, k700, and k500 layers. Now
114 ! you see WHY? we wanted Q on pressure levels, it all is beginning
115 ! to make sense.
116
117 t850 = t(k850) * (1.0 + 0.608 * q(k850))
118 t700 = t(k700) * (1.0 + 0.608 * q(k700))
119 t500 = t(k500) * (1.0 + 0.608 * q(k500))
120
121 ! Compute two lapse rates between these three levels. These are
122 ! environmental values for each (i,j).
123
124 gamma78 = LOG(t850 / t700) / LOG (p(k850) / p(k700))
125 gamma57 = LOG(t700 / t500) / LOG (p(k700) / p(k500))
126
127 if ((psfc - p(k850) - 10000.0) >= 0.0) then
128 t1 = t850
129 else if ((psfc - p(k850)) >= 0.0) then
130 t1 = t700 * (p1 / p(k700)) ** gamma78
131 else if ((psfc - p(k700)) >= 0.0) then
132 t1 = t500 * (p1 / p(k500)) ** gamma57
133 else
134 t1 = t500
135 end if
136
137 ! From our temperature way up in the air, we extrapolate down to
138 ! the sea level to get a guess at the sea level temperature.
139
140 tslv = t1 * (pslv / p1) ** (gammarg)
141
142 ! The new surface temperature is computed from the with new sea level
143 ! temperature, just using the elevation and a lapse rate. This lapse
144 ! rate is -6.5 K/km.
145
146 if (present (t_obs)) then
147 tsfc = t_obs
148 else
149 tsfc = tslv - gamma * h_obs
150 end if
151
152 ! A correction to the sea-level temperature, in case it is too warm.
153
154 TFIXED = TC - 0.005 * (TSFC - TC) ** 2
155
156 L1 = tslv .LT. tc
157 L2 = tsfc .LE. tc
158 L3 = .NOT. L1
159 if (L2 .AND. L3) then
160 tslv = tc
161 else if ((.NOT. L2) .AND. L3) then
162 tslv = tfixed
163 end if
164
165 ! Finally, we can get to the surface pressure.
166
167 p1 = -h_obs * g / (rov2 * (tsfc + tslv))
168 psfc = pslv * EXP(p1)
169
170 ! Surface pressure and sea-level pressure are the same at sea level.
171
172 if (ABS (h_obs) < 0.1) psfc = pslv
173
174 if (trace_use_dull) call da_trace_exit("da_sfcprs")
175
176 end subroutine da_sfcprs
177
178