Ticket #658: ropp_fm_abel_refl.f90

File ropp_fm_abel_refl.f90, 12.6 KB (added by Ian Culverwell, 6 years ago)

ropp_fm_abel_refl.f90

Line 
1! $Id: ropp_fm_abel_refl.f90 $
2
3SUBROUTINE 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)
177write(*,*) 'in abel_refl values of r1 and r2 ',r1, r2
178write(*,*) 'in abel_refl values of refrac1 refrac2 ', refrac(1), refrac(2)
179write(*,*) '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)
190write(*,*) '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
355END SUBROUTINE ropp_fm_abel_refl
356