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 assimilation
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, dimension (kts:kte), intent (in) :: height, pressure, temp
20
21 real, dimension (kts:kte) :: prs_mb
22 ! calculated but never used
23 ! real :: dth_12, dth_21, dth_sfc, dth_obs
24 ! real :: dhe_12, dhe_21, dhe_sfc1, dhe_obs1, dhe_sfc2, dhe_obs2
25 real :: dth_21, dth_sfc, dth_obs
26 real :: dhe_12, dhe_sfc1, dhe_obs1, dhe_sfc2, dhe_obs2
27 real :: th_100mb, th_200mb, th_obs, th_sfc
28 real :: th_obs_int, th_sfc_int
29 real :: pdif, rcp
30 integer :: k_100mb, k_200mb, jk
31 integer :: inc_100mb, inc_200mb
32
33 rcp = gas_constant/cp
34
35 ! 1.Find levels: model surface + 100hpa and model surface + 200hpa ar obs loc
36 ! ===========================================================================
37
38 ! 1.1 Convert model pressure profile from Pa to hPa
39
40 prs_mb = pressure / 100.
41
42 ! 1.2 Find levels surface + 100hPA
43
44 inc_100mb = 100.0
45 k_100mb = kts
46
47 do jk = kts+1, kte
48 pdif = prs_mb (kts) - prs_mb (jk)
49 if (pdif .GE. inc_100mb) then
50 k_100mb = jk
51 exit
52 end if
53 end do
54
55 ! 1.2 Find levels surface + 200hPA
56
57 inc_200mb = 200.0
58 k_200mb = kts
59
60 do jk = kts+1, kte
61 pdif = prs_mb (kts) - prs_mb (jk)
62 if (pdif .GE. inc_200mb) then
63 k_200mb = jk
64 exit
65 end if
66 end do
67
68 ! 1.3 Check consistency
69
70 if ((k_100mb .LE. kts) .OR. (k_200mb .LE. kts) .OR. &
71 (k_200mb .LT. k_100mb)) then
72 write (unit=message(1),fmt='(A)') ' CANNOT FinD SFC + 200hPA AND SFC + 100hPA'
73 write (unit=message(2),fmt='(A,I2,A,F10.3)') ' P (',k_200mb,') = ',prs_mb (k_200mb)
74 write (unit=message(3),fmt='(A,I2,A,F10.3)') ' P (',k_100mb,') = ',prs_mb (k_100mb)
75 write (unit=message(4),fmt='(A,F10.3)') ' P_SFC = ', prs_mb (kts)
76 call da_warning(__FILE__,__LINE__,message(1:4))
77 val = missing_r
78 return
79 end if
80
81 ! 2. POTENTIAL TEMPERATURE
82 ! =========================
83
84 ! 2.1 Potential temperature at 100hPa above model surface
85
86 th_100mb = temp (k_100mb) * (1000. / prs_mb (k_100mb))**rcp
87
88 ! 2.2 Potential temperature at 200hPa above model surface
89
90 th_200mb = temp (K_200mb) * (1000. / prs_mb (k_200mb))**rcp
91
92 ! 2.3 Potential temperature at observation location
93
94 th_obs = to * (1000. / (po/100.)) ** rcp
95
96
97 ! 3. LAPSE RATE BETWEEN SURFACE+200hPA AND SURFACE+100hPa
98 ! =========================================================
99
100 ! 3.1 Potential temperature increment
101
102 dth_21 = th_100mb - th_200mb
103 ! never used
104 ! dth_12 = th_200mb - th_100mb
105
106 ! 3.1 Height increments
107
108 ! never used
109 ! dhe_21 = height (k_100mb)- height (k_200mb)
110 dhe_sfc1 = height (k_100mb)- height (kts)
111 dhe_obs1 = height (k_100mb)- ho
112
113 dhe_12 = height (k_200mb)- height (k_100mb)
114 dhe_sfc2 = height (k_200mb)- height (kts)
115 dhe_obs2 = height (k_200mb)- ho
116
117 ! 3.2 Extrapolated potential temperature at model surface and observation loc
118
119 th_sfc_int = th_100mb + (dth_21/dhe_12) * dhe_sfc1
120 th_obs_int = th_100mb + (dth_21/dhe_12) * dhe_obs1
121
122
123 ! 4. BRING TEMPERATURE ONTO MODEL SURFACE
124 ! ========================================
125
126 ! 4.1 Difference at observation locations
127
128 dth_obs = th_obs_int - th_obs
129
130 ! 4.2 Difference at model surface
131
132 dth_sfc = (dhe_sfc2/dhe_obs2) * dth_obs
133
134 ! 4.3 Potentiel temperature brought to model surface
135
136 th_sfc = th_sfc_int - dth_sfc
137
138 ! 4.3 Corresponding Temperature
139
140 val = th_sfc * (prs_mb (kts) / 1000.)**rcp
141
142 end subroutine da_intpsfc_tem
143
144