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