1 | ! $Id: ropp_fm_abel_refl.f90 $
|
---|
2 |
|
---|
3 | SUBROUTINE ropp_fm_abel_refl(nr, refrac, temp, roc, und, Tgrad_oper, impact, bangle)
|
---|
4 |
|
---|
5 | !****s* BendingAngle/ropp_fm_abel_refl *
|
---|
6 | !
|
---|
7 | ! NAME
|
---|
8 | ! ropp_fm_abel_refl - Forward model calculating a one dimensional bending
|
---|
9 | ! angle profile from refractivity / impact parameter profile
|
---|
10 | ! at state vector levels using a modified Fast Abel Transform
|
---|
11 | ! valid for impact parameters under the surface, which
|
---|
12 | ! correspond to signals reflected off the surface.
|
---|
13 | !
|
---|
14 | ! SYNOPSIS
|
---|
15 | ! call ropp_fm_abel_refl(nr, refrac, temp, roc, und, Tgrad_oper, impact, bangle)
|
---|
16 | !
|
---|
17 | ! DESCRIPTION
|
---|
18 | ! This routine calculates bending angles for reflected signals, evaluated
|
---|
19 | ! at a given set of impact parameters below the surface's impact parameter.
|
---|
20 | ! It uses a vertical profile of rafractivity given at
|
---|
21 | ! the state vector's set of
|
---|
22 | ! x = nr levels.
|
---|
23 | !
|
---|
24 | ! INPUTS
|
---|
25 | ! real(wp), dimension(:) :: nr ! x = nr product
|
---|
26 | ! real(wp), dimension(:) :: refrac ! Refractivity values
|
---|
27 | ! real(wp), dimension(:) :: temp ! temperature values
|
---|
28 | ! real(wp), :: roc ! radius of curvature
|
---|
29 | ! real(wp), :: und ! undulation
|
---|
30 | ! LOGICAL :: Tgrad_oper ! Use temp. gradient oper.
|
---|
31 | ! real(wp), dimension(:) :: impact ! Impact parameters
|
---|
32 | !
|
---|
33 | ! OUTPUT
|
---|
34 | ! real(wp), dimension(:) :: bangle ! Forward modelled bending angles
|
---|
35 | !
|
---|
36 | ! NOTES
|
---|
37 | ! The model's equation can be found in Aparicio et al., 2017.
|
---|
38 | ! The interpolation of bending angles calculated at the state vector's
|
---|
39 | ! geopotential levels to the observation vector's impact parameters is
|
---|
40 | ! carried out assuming that bending angle varies exponentially with
|
---|
41 | ! impact parameter.
|
---|
42 | !
|
---|
43 | ! SEE ALSO
|
---|
44 | ! ropp_fm_types
|
---|
45 | ! ropp_fm_bangle_1d_refl
|
---|
46 | ! ropp_fm_bangle_1d
|
---|
47 | ! ropp_fm_abel_ad
|
---|
48 | ! ropp_fm_abel_tl
|
---|
49 | !
|
---|
50 | ! AUTHOR
|
---|
51 | ! Institute of Space Studies of Catalonia (IEEC), Barcelona, Spain, based on
|
---|
52 | ! original code (no-reflection) by
|
---|
53 | ! Met Office, Exeter, UK.
|
---|
54 | ! Any comments on this software should be given via the ROM SAF
|
---|
55 | ! Helpdesk at http://www.romsaf.org
|
---|
56 | !
|
---|
57 | ! COPYRIGHT
|
---|
58 | ! (c) EUMETSAT. All rights reserved.
|
---|
59 | ! For further details please refer to the file COPYRIGHT
|
---|
60 | ! which you should have received as part of this distribution.
|
---|
61 | !
|
---|
62 | !****
|
---|
63 |
|
---|
64 | !-------------------------------------------------------------------------------
|
---|
65 | ! 1. Declarations
|
---|
66 | !-------------------------------------------------------------------------------
|
---|
67 | USE typesizes, ONLY: wp => EightByteReal
|
---|
68 | USE ropp_utils, ONLY: ropp_MDFV, ropp_ZERO, ropp_ZDTV
|
---|
69 | USE ropp_fm_constants, ONLY: pi, imp_ht_min
|
---|
70 |
|
---|
71 | IMPLICIT NONE
|
---|
72 |
|
---|
73 | REAL(wp), DIMENSION(:), INTENT(in) :: nr ! x = nr product
|
---|
74 | REAL(wp), DIMENSION(:), INTENT(in) :: refrac ! Refractivity
|
---|
75 | REAL(wp), DIMENSION(:), INTENT(in) :: temp ! Temperature
|
---|
76 | REAL(wp) , INTENT(in) :: roc ! Radius of curvature
|
---|
77 | REAL(wp) , INTENT(in) :: und ! Undulation
|
---|
78 | LOGICAL , INTENT(in) :: Tgrad_oper ! Use temp. gradient oper.
|
---|
79 |
|
---|
80 | REAL(wp), DIMENSION(:), INTENT(in) :: impact ! Impact parameter
|
---|
81 | REAL(wp), DIMENSION(:), INTENT(out) :: bangle ! Bending angle
|
---|
82 |
|
---|
83 | REAL(wp), DIMENSION(:), ALLOCATABLE :: kval ! Exponential decay rate
|
---|
84 | REAL(wp), DIMENSION(:), ALLOCATABLE :: beta ! Temperature gradient
|
---|
85 | REAL(wp), DIMENSION(:), ALLOCATABLE :: nr_mid ! Average nr product
|
---|
86 | REAL(wp), DIMENSION(:), ALLOCATABLE :: temp_mid ! Average temp. of two levels
|
---|
87 | REAL(wp), DIMENSION(:), ALLOCATABLE :: dval ! Useful in BA computation
|
---|
88 | REAL(wp) :: t_upper ! Upper bound integral
|
---|
89 | REAL(wp) :: t_lower ! Lower bound integral
|
---|
90 | REAL(wp) :: refrac_low ! Refrac at lower level
|
---|
91 | REAL(wp) :: nr_low ! x=nr at lower level
|
---|
92 | REAL(wp) :: zed ! Impact height
|
---|
93 |
|
---|
94 |
|
---|
95 | REAL(wp) :: integral_diff ! Integral approximation
|
---|
96 | REAL(wp) :: erf_up, erf_low
|
---|
97 | REAL(wp) :: int_up, int_low
|
---|
98 | REAL(wp) :: zt
|
---|
99 | REAL(wp) :: dn_dx,p1,p2,p3
|
---|
100 | REAL(wp) :: as ! surf. impact parameter
|
---|
101 | REAL(wp) :: r1,r2,AA,BB ! aux. variables
|
---|
102 |
|
---|
103 | REAL(wp), PARAMETER :: a = 0.3480242_wp
|
---|
104 | REAL(wp), PARAMETER :: b = 0.0958798_wp
|
---|
105 | REAL(wp), PARAMETER :: c = 0.7478556_wp
|
---|
106 |
|
---|
107 | INTEGER :: n_lev, n_lower, n_impact
|
---|
108 | INTEGER :: i, i_bot, l, n_top
|
---|
109 | !-------------------------------------------------------------------------------
|
---|
110 | ! 2. Useful variables
|
---|
111 | !-------------------------------------------------------------------------------
|
---|
112 |
|
---|
113 | n_lev = SIZE(nr)
|
---|
114 | n_impact = SIZE(impact)
|
---|
115 |
|
---|
116 | ALLOCATE(kval(n_lev-1))
|
---|
117 | ALLOCATE(beta(n_lev-1))
|
---|
118 | ALLOCATE(nr_mid(n_lev-1))
|
---|
119 | ALLOCATE(temp_mid(n_lev-1))
|
---|
120 | ALLOCATE(dval(n_lev-1))
|
---|
121 |
|
---|
122 | !-------------------------------------------------------------------------------
|
---|
123 | ! 3. Calculate lowest usable level (because of superrefraction)
|
---|
124 | !-------------------------------------------------------------------------------
|
---|
125 | n_lower = 1
|
---|
126 | DO i = n_lev, 2, -1
|
---|
127 | IF (nr(i) - nr(i-1) < 10.0_wp) THEN
|
---|
128 | n_lower = i
|
---|
129 | EXIT
|
---|
130 | ENDIF
|
---|
131 | ENDDO
|
---|
132 | IF (n_lower > 1) THEN
|
---|
133 | print *, "WARNING SUPERREFRACTION: current implementation"
|
---|
134 | print *, "of ropp_fm_abel_refl might not work properly "
|
---|
135 | print *, "in superrefractive profiles."
|
---|
136 | ENDIF
|
---|
137 |
|
---|
138 | !-------------------------------------------------------------------------------
|
---|
139 | ! 4. Calculate exponential decay rate and temp. gradient between levels
|
---|
140 | !-------------------------------------------------------------------------------
|
---|
141 |
|
---|
142 | DO i = 1, n_lev - 1
|
---|
143 |
|
---|
144 | kval(i) = LOG(refrac(i)/refrac(i+1)) / MAX(1.0_wp, (nr(i+1)-nr(i)))
|
---|
145 | kval(i) = MAX(1.0e-6_wp, kval(i))
|
---|
146 |
|
---|
147 | ! limit the maximum kval so that refractivity gradient is ~ half critical value
|
---|
148 |
|
---|
149 | kval(i) = MIN(kval(i),0.157_wp/refrac(i))
|
---|
150 |
|
---|
151 | ! beta is the temperature gradient
|
---|
152 |
|
---|
153 | beta(i) = (temp(i+1) - temp(i)) / MAX(1.0_wp, (nr(i+1)-nr(i)))
|
---|
154 |
|
---|
155 | ! mean nr
|
---|
156 |
|
---|
157 | nr_mid(i) = 0.5_wp*(nr(i) + nr(i+1))
|
---|
158 |
|
---|
159 | ! mean temp.
|
---|
160 |
|
---|
161 | temp_mid(i) = 0.5_wp*(temp(i) + temp(i+1))
|
---|
162 |
|
---|
163 | dval(i) = (nr(i) -nr_mid(i))**2
|
---|
164 | ENDDO
|
---|
165 |
|
---|
166 | !-------------------------------------------------------------------------------
|
---|
167 | ! 5. Calculate bending angles for observational heights
|
---|
168 | !-------------------------------------------------------------------------------
|
---|
169 |
|
---|
170 | bangle(:) = ropp_MDFV
|
---|
171 |
|
---|
172 | ! 5.1 Find top output impact value (surface level or closest below)
|
---|
173 | ! Extrapolating from bottom level of the state vector to surface
|
---|
174 |
|
---|
175 | r1 = nr(1)/(refrac(1)*1E-6_wp+1.0_wp)
|
---|
176 | r2 = nr(2)/(refrac(2)*1E-6_wp+1.0_wp)
|
---|
177 | write(*,*) 'in abel_refl values of r1 and r2 ',r1, r2
|
---|
178 | write(*,*) 'in abel_refl values of refrac1 refrac2 ', refrac(1), refrac(2)
|
---|
179 | write(*,*) 'in abel_refl values of roc und', roc, und, roc+und
|
---|
180 | BB = (LOG(refrac(1))-LOG(refrac(2)))/(r2-r1)
|
---|
181 | !AA = EXP(LOG(refrac(1))+BB*(r2-roc-und))
|
---|
182 | AA = refrac(1) / (exp(-BB*(r1-roc-und)))
|
---|
183 | write(*,*) 'in abel_refl values of BB AA',BB,AA
|
---|
184 | ! extrapolation of refractivity towards the surface:
|
---|
185 | ! N=AA*exp(-BB*(r-rsurface))
|
---|
186 | ! where AA is refractivity at the surface level.
|
---|
187 | ! Then impact parameter at the surface is:
|
---|
188 |
|
---|
189 | as = (AA*1.0E-6_wp+1.0_wp)*(roc+und)
|
---|
190 | write(*,*) 'in abel_refl value of impact at the surface ',as
|
---|
191 |
|
---|
192 | ! which limits the uppermost reflected impact parameter
|
---|
193 | ! to be evaluated within the given impact vector
|
---|
194 |
|
---|
195 | DO l = 1, n_lev
|
---|
196 | IF (nr(l) > as ) THEN
|
---|
197 | n_top = l-1
|
---|
198 | EXIT
|
---|
199 | ENDIF
|
---|
200 | ENDDO
|
---|
201 | write(*,*) 'n_top in abel and n_lev',n_top, n_lev
|
---|
202 | IF ( n_top < 1 ) THEN
|
---|
203 | print *, "WARNING: impact parameters to be evaluated as reflected signals must be below the surface one."
|
---|
204 | ENDIF
|
---|
205 |
|
---|
206 |
|
---|
207 | ! 5.1 Bottom state vector level is always the lowest one
|
---|
208 | ! In fact, it should be the surface level, but here
|
---|
209 | ! restricted to cases free of superrefractivity.
|
---|
210 | ! i_bot = n_lower which should be 1 except for
|
---|
211 | ! superrefractive profiles
|
---|
212 | ! ------------------------------------------------------
|
---|
213 |
|
---|
214 | !i_bot = n_lower ! NO, now set to 1 always
|
---|
215 | i_bot = 1
|
---|
216 |
|
---|
217 | ! 5.2 Loop over all levels above
|
---|
218 | ! ------------------------------
|
---|
219 | ! DO l = 1, n\_impact -- now only evaluating impact parameters up to surf
|
---|
220 | IF (n_top >= 1) THEN
|
---|
221 | DO l = 1, n_top
|
---|
222 |
|
---|
223 | bangle(l) = ropp_ZERO
|
---|
224 |
|
---|
225 |
|
---|
226 | DO i = i_bot, n_lev - 1
|
---|
227 |
|
---|
228 | ! 5.2.1 Values of refractivity and impact parameter at lower level
|
---|
229 | ! ----------------------------------------------------------------
|
---|
230 |
|
---|
231 | refrac_low = refrac(i)
|
---|
232 | nr_low = nr(i)
|
---|
233 |
|
---|
234 |
|
---|
235 | IF (refrac(i+1)-refrac(i) > - ropp_ZDTV ) THEN
|
---|
236 |
|
---|
237 | ! 5.2.2 If the refractivity gradient is +ve with height
|
---|
238 | ! -------------------------------------------------------------------
|
---|
239 |
|
---|
240 | ! This will handle the cases where the refractivity goes up with height
|
---|
241 |
|
---|
242 | dn_dx = (refrac(i+1)-refrac(i))/(nr(i+1)-nr(i))
|
---|
243 |
|
---|
244 | t_upper = SQRT( nr(i+1)-impact(l))
|
---|
245 |
|
---|
246 | t_lower = 0.0_wp
|
---|
247 |
|
---|
248 | IF (i > i_bot) t_lower = SQRT( nr(i)-impact(l))
|
---|
249 |
|
---|
250 | bangle(l) = bangle(l) - &
|
---|
251 | 2.0E-6_wp*SQRT(2.0_wp*impact(l))* dn_dx*(t_upper-t_lower)
|
---|
252 |
|
---|
253 | ELSE
|
---|
254 |
|
---|
255 | ! 5.2.3 Upper and lower bounds of the "normal" integral
|
---|
256 | ! -------------------------------------------------------------------
|
---|
257 |
|
---|
258 | t_upper = SQRT(kval(i) * (nr(i+1) - impact(l)))
|
---|
259 |
|
---|
260 | IF (i == i_bot) THEN
|
---|
261 | t_lower = 0.0_wp
|
---|
262 | ELSE
|
---|
263 | t_lower = SQRT(kval(i) * (nr(i) - impact(l)))
|
---|
264 | ENDIF
|
---|
265 |
|
---|
266 | ! 5.2.4 Error functions
|
---|
267 | ! --------------
|
---|
268 |
|
---|
269 | ! Approximate error function with polynomial
|
---|
270 |
|
---|
271 | zt = 1.0_wp / (1.0_wp + 0.47047_wp * t_lower)
|
---|
272 | erf_low = 1.0_wp - (a-(b-c*zt)*zt) * zt * EXP(-(t_lower*t_lower))
|
---|
273 | zt = 1.0_wp / (1.0_wp + 0.47047_wp*t_upper)
|
---|
274 | erf_up = 1.0_wp - (a-(b-c*zt)*zt) * zt * EXP(-(t_upper*t_upper))
|
---|
275 |
|
---|
276 | IF (i == n_lev-1) erf_up = 1.0_wp ! limit at infinity
|
---|
277 |
|
---|
278 |
|
---|
279 | ! 5.2.5 New terms for integral that now allows kval to vary within layer
|
---|
280 | ! ------------------------------------------------------------------
|
---|
281 |
|
---|
282 | ! These p1, p2,p3 values correspond the case where dT/dx = beta = 0, i.e kval is constant!
|
---|
283 |
|
---|
284 | p1 = kval(i)
|
---|
285 | p2 = 0.0_wp
|
---|
286 | p3 = 0.0_wp
|
---|
287 |
|
---|
288 | ! impact height of level
|
---|
289 |
|
---|
290 | zed = nr(i) - roc
|
---|
291 |
|
---|
292 | IF ( i < n_lev-1 .AND. zed > imp_ht_min .AND. Tgrad_oper) THEN
|
---|
293 |
|
---|
294 | ! compute the "p" values for temp. gradient beta
|
---|
295 | p1 = kval(i)*(1.0_wp + beta(i)/temp_mid(i)* &
|
---|
296 | (0.5_wp*kval(i)*((impact(l)-nr_mid(i))**2 - dval(i)) - &
|
---|
297 | (impact(l)-nr_mid(i))))
|
---|
298 |
|
---|
299 | p2 = kval(i)*beta(i)/temp_mid(i)* &
|
---|
300 | (kval(i)*(impact(l)-nr_mid(i))-1.0_wp)
|
---|
301 |
|
---|
302 | p3 = 0.5_wp*kval(i)**2*beta(i)/temp_mid(i)
|
---|
303 |
|
---|
304 | ENDIF
|
---|
305 |
|
---|
306 | int_up = SQRT(pi/kval(i))* &
|
---|
307 | (p1 + 0.5_wp/kval(i)*(p2+1.5_wp*p3/kval(i)))*erf_up
|
---|
308 |
|
---|
309 | int_up = int_up - EXP(-kval(i)*(nr(i+1)-impact(l)))* &
|
---|
310 | SQRT(nr(i+1)-impact(l))/kval(i)*(p2 + p3* ( &
|
---|
311 | (nr(i+1)-impact(l))+1.5_wp/kval(i)))
|
---|
312 |
|
---|
313 | ! lower limit of integral
|
---|
314 |
|
---|
315 | int_low = 0.0_wp
|
---|
316 |
|
---|
317 | IF (i > i_bot) THEN
|
---|
318 |
|
---|
319 | int_low = SQRT(pi/kval(i))* &
|
---|
320 | (p1 + 0.5_wp/kval(i)*(p2+ 1.5_wp*p3/kval(i)))*erf_low
|
---|
321 |
|
---|
322 | int_low = int_low - EXP(-kval(i)*(nr(i)-impact(l)))* &
|
---|
323 | SQRT(nr(i)-impact(l))/kval(i)*(p2 + p3* ( &
|
---|
324 | (nr(i)-impact(l))+1.5_wp/kval(i)))
|
---|
325 |
|
---|
326 | ENDIF
|
---|
327 |
|
---|
328 | integral_diff = int_up - int_low
|
---|
329 |
|
---|
330 | ! 5.2.6 Bending angle value
|
---|
331 | ! -------------------------
|
---|
332 |
|
---|
333 | bangle(l) = bangle(l) &
|
---|
334 | + 1.0e-6_wp * SQRT(2.0_wp*impact(l)) &
|
---|
335 | * refrac_low * EXP(kval(i) * (nr_low - impact(l))) &
|
---|
336 | * integral_diff
|
---|
337 |
|
---|
338 | ENDIF
|
---|
339 |
|
---|
340 | ENDDO
|
---|
341 |
|
---|
342 | ! 5.3 Adding the term due to the reflection angle
|
---|
343 | ! -----------------------------------------------
|
---|
344 | bangle(l) = bangle(l) - 2.0_wp * ACOS(impact(l)/as)
|
---|
345 |
|
---|
346 | ENDDO
|
---|
347 | ENDIF
|
---|
348 |
|
---|
349 | DEALLOCATE(beta)
|
---|
350 | DEALLOCATE(nr_mid)
|
---|
351 | DEALLOCATE(temp_mid)
|
---|
352 | DEALLOCATE(dval)
|
---|
353 | DEALLOCATE(kval)
|
---|
354 |
|
---|
355 | END SUBROUTINE ropp_fm_abel_refl
|
---|
356 |
|
---|