da_intpsfc_tem.inc
References to this file elsewhere.
1 subroutine da_intpsfc_tem (val, ho, po, to, height, pressure, temp, kts, kte)
2
3 !-----------------------------------------------------------------------
4 ! Purpose: Correct temperature between two levels.
5 !
6 ! Reference:
7 ! ---------
8 ! The use of surface observations in four dimensional data assmiilation
9 ! Using a mesoscale model.
10 ! Ruggiero et al., 1996, Monthly Weather Review, Volume 124, 1018-1033
11 !
12 !----------------------------------------------------------------------------
13
14 implicit none
15
16 real, intent (out) :: val
17 integer, intent (in) :: kts, kte
18 real, intent (in) :: ho, po, to
19 real, intent (in) :: height(kts:kte)
20 real, intent (in) :: pressure(kts:kte)
21 real, intent (in) :: temp(kts:kte)
22
23 real :: prs_mb(kts:kte)
24 ! calculated but never used
25 ! real :: dth_12, dth_21, dth_sfc, dth_obs
26 ! real :: dhe_12, dhe_21, dhe_sfc1, dhe_obs1, dhe_sfc2, dhe_obs2
27 real :: dth_21, dth_sfc, dth_obs
28 real :: dhe_12, dhe_sfc1, dhe_obs1, dhe_sfc2, dhe_obs2
29 real :: th_100mb, th_200mb, th_obs, th_sfc
30 real :: th_obs_int, th_sfc_int
31 real :: pdif, rcp
32 integer :: k_100mb, k_200mb, jk
33 integer :: inc_100mb, inc_200mb
34
35 if (trace_use_dull) call da_trace_entry("da_intpsfc_tem")
36
37 rcp = gas_constant/cp
38
39 ! 1.Find levels: model surface + 100hpa and model surface + 200hpa ar obs loc
40 ! ===========================================================================
41
42 ! 1.1 Convert model pressure profile from Pa to hPa
43
44 prs_mb = pressure / 100.0
45
46 ! 1.2 Find levels surface + 100hPA
47
48 inc_100mb = 100.0
49 k_100mb = kts
50
51 do jk = kts+1, kte
52 pdif = prs_mb (kts) - prs_mb (jk)
53 if (pdif .GE. inc_100mb) then
54 k_100mb = jk
55 exit
56 end if
57 end do
58
59 ! 1.2 Find levels surface + 200hPA
60
61 inc_200mb = 200.0
62 k_200mb = kts
63
64 do jk = kts+1, kte
65 pdif = prs_mb (kts) - prs_mb (jk)
66 if (pdif .GE. inc_200mb) then
67 k_200mb = jk
68 exit
69 end if
70 end do
71
72 ! 1.3 Check consistency
73
74 if ((k_100mb .LE. kts) .OR. (k_200mb .LE. kts) .OR. &
75 (k_200mb .LT. k_100mb)) then
76 write (unit=message(1),fmt='(A)') ' Cannot find sfc + 200hPa and sfc + 100hPa'
77 write (unit=message(2),fmt='(A,I2,A,F10.3)') ' P (',k_200mb,') = ',prs_mb (k_200mb)
78 write (unit=message(3),fmt='(A,I2,A,F10.3)') ' P (',k_100mb,') = ',prs_mb (k_100mb)
79 write (unit=message(4),fmt='(A,F10.3)') ' P_SFC = ', prs_mb (kts)
80 call da_warning(__FILE__,__LINE__,message(1:4))
81 val = missing_r
82 return
83 end if
84
85 ! 2. potential temperature
86 ! =========================
87
88 ! 2.1 Potential temperature at 100hPa above model surface
89
90 th_100mb = temp (k_100mb) * (1000.0 / prs_mb (k_100mb))**rcp
91
92 ! 2.2 Potential temperature at 200hPa above model surface
93
94 th_200mb = temp (K_200mb) * (1000.0 / prs_mb (k_200mb))**rcp
95
96 ! 2.3 Potential temperature at observation location
97
98 th_obs = to * (1000.0 / (po/100.0)) ** rcp
99
100
101 ! 3. lapse rate between surface+200hpa and surface+100hpa
102 ! =========================================================
103
104 ! 3.1 Potential temperature increment
105
106 dth_21 = th_100mb - th_200mb
107 ! never used
108 ! dth_12 = th_200mb - th_100mb
109
110 ! 3.1 Height increments
111
112 ! never used
113 ! dhe_21 = height (k_100mb)- height (k_200mb)
114 dhe_sfc1 = height (k_100mb)- height (kts)
115 dhe_obs1 = height (k_100mb)- ho
116
117 dhe_12 = height (k_200mb)- height (k_100mb)
118 dhe_sfc2 = height (k_200mb)- height (kts)
119 dhe_obs2 = height (k_200mb)- ho
120
121 ! 3.2 Extrapolated potential temperature at model surface and observation loc
122
123 th_sfc_int = th_100mb + (dth_21/dhe_12) * dhe_sfc1
124 th_obs_int = th_100mb + (dth_21/dhe_12) * dhe_obs1
125
126
127 ! 4. Bring temperature onto model surface
128 ! ========================================
129
130 ! 4.1 Difference at observation locations
131
132 dth_obs = th_obs_int - th_obs
133
134 ! 4.2 Difference at model surface
135
136 dth_sfc = (dhe_sfc2/dhe_obs2) * dth_obs
137
138 ! 4.3 Potentiel temperature brought to model surface
139
140 th_sfc = th_sfc_int - dth_sfc
141
142 ! 4.3 Corresponding Temperature
143
144 val = th_sfc * (prs_mb (kts) / 1000.0)**rcp
145
146 if (trace_use_dull) call da_trace_exit("da_intpsfc_tem")
147
148 end subroutine da_intpsfc_tem
149
150