Ticket #255: ropp_pp_ionospheric_correction.f90

File ropp_pp_ionospheric_correction.f90, 19.0 KB (added by Ian Culverwell, 14 years ago)
Line 
1! $Id: ropp_pp_ionospheric_correction.f90 2511 2010-05-27 15:22:27Z frhl $
2
3SUBROUTINE ropp_pp_ionospheric_correction(impact_L1, ba_L1, impact_L2, ba_L2, &
4 impact_LM, ba_LM, config, impact_LC, ba_LC, diag_out, WLC) !HGL!
5
6!****s* IonosphericCorrection/ropp_pp_ionospheric_correction *
7!
8! NAME
9! ropp_pp_ionospheric_correction -
10! Calculate neutral and ionospheric bending angle
11! profile on L1 impact heights from L1/L2 bending angles
12!
13! SYNOPSIS
14! call ropp_pp_ionospheric_correction(impact_L1, ba_L1,
15! impact_L2, ba_L2,
16! impact_LM, ba_LM,
17! config, impact_LC, ba_LC, diag_out, WLC) !HGL!
18!
19! DESCRIPTION
20! This routine calculates bending angles at a given set of impact parameters
21! from vertical profiles of bending angles at the two measurement
22! frequencies (channels) L1 and L2.
23!
24! INPUTS
25! real(wp), dimension(:) :: impact_L1 ! Impact parameters of channel L1 (m)
26! real(wp), dimension(:) :: ba_L1 ! Bending angles for channel L1 (rad)
27! real(wp), dimension(:) :: impact_L2 ! Impact parameters of channel L2 (m)
28! real(wp), dimension(:) :: ba_L2 ! Bending angles for channel L2 (rad)
29! real(wp), dimension(:) :: impact_LM ! Model impact parameters (m)
30! real(wp), dimension(:) :: ba_LM ! Model bending angles (rad)
31! type(ppConfig) :: config ! Configuration parameters
32!
33! OUTPUT
34! real(wp), dimension(:) :: impact_LC ! Impact parameters of channel L1
35! real(wp), dimension(:) :: ba_LC ! Corrected bending angles
36! real(wp), dimension(:) :: WLC ! Weight of the linear combination of L1 and L2 !HGL!
37! type(ppDiag) :: diag_out ! Output diagnostics
38!
39! NOTES
40! Method:
41! 1. Calculation of strongly smoothed ionospheric signal
42! (using the external scale). Further deviations from
43! this signal are calculated.
44! 2. Estimation of ionospheric signal and noise covariances
45! using the highest part (> 50 km) of the occultation.
46! 3. Calculation of relative mean deviation of neutral refraction
47! from the model refraction using signal at heights 12-35 km.
48! 4. Optimal linear combination for the same impact parameters using the
49! covariances.
50!
51! REFERENCES
52! M.E. Gorbunov
53! Ionospheric correction and statistical optimization of radio
54! occultation data
55! Radio Science, 37(5), 1084
56!
57! AUTHOR
58! M Gorbunov, Russian Academy of Sciences, Russia.
59! Any comments on this software should be given via the GRAS SAF
60! Helpdesk at http://www.grassaf.org
61!
62! COPYRIGHT
63! Copyright (c) 1998-2010 Michael Gorbunov <michael.gorbunov@zmaw.de>
64! For further details please refer to the file COPYRIGHT
65! which you should have received as part of this distribution.
66!
67!****
68
69!-------------------------------------------------------------------------------
70! 1. Declarations
71!-------------------------------------------------------------------------------
72
73 USE typesizes, ONLY: wp => EightByteReal
74 USE ropp_pp, not_this => ropp_pp_ionospheric_correction
75 USE ropp_pp_types, ONLY: PPConfig, PPDiag
76
77 IMPLICIT NONE
78
79 REAL(wp), DIMENSION(:), INTENT(in) :: impact_L1 ! L1 impact parameters (m)
80 REAL(wp), DIMENSION(:), INTENT(in) :: ba_L1 ! L1 bending angles (rad)
81 REAL(wp), DIMENSION(:), INTENT(in) :: impact_L2 ! L2 impact parameters (m)
82 REAL(wp), DIMENSION(:), INTENT(in) :: ba_L2 ! L2 bending angles (rad)
83 REAL(wp), DIMENSION(:), INTENT(in) :: impact_LM ! Model impact parameters
84 REAL(wp), DIMENSION(:), INTENT(in) :: ba_LM ! Model bending angles
85 TYPE(ppConfig), INTENT(inout) :: config ! Configuration parameters
86 REAL(wp), DIMENSION(:), INTENT(out) :: impact_LC ! Output impact parameters
87 REAL(wp), DIMENSION(:), INTENT(out) :: ba_LC ! Corrected bending angles
88 TYPE(ppDiag), INTENT(inout) :: diag_out ! Output diagnostics
89 REAL(wp), DIMENSION(:), INTENT(inout), OPTIONAL :: WLC ! Weight of the linear combination of L1 and L2 !HGL!
90
91 REAL(wp), DIMENSION(:), ALLOCATABLE :: impact_LH ! Hi-res impact grid
92 REAL(wp), DIMENSION(:), ALLOCATABLE :: ba_L1H ! L1 bangle on impact_LH
93 REAL(wp), DIMENSION(:), ALLOCATABLE :: ba_L2H ! L2 bangle on impact_LH
94 REAL(wp), DIMENSION(:), ALLOCATABLE :: ba_LMH ! LM bangle on impact_LH
95 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ba_diff ! bangle differences
96 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ba_low ! low-freq ba_ob-ba_model
97 REAL(wp), DIMENSION(:), ALLOCATABLE :: ba_is ! smoothed ionospheric ba
98 REAL(wp), DIMENSION(:), ALLOCATABLE :: ba_ion ! ionospheric bangle
99 REAL(wp), DIMENSION(:), ALLOCATABLE :: ba_neut ! neutral bangle
100 REAL(wp), DIMENSION(:), ALLOCATABLE :: ba_nfilt ! filtered neutral bangle
101 REAL(wp), DIMENSION(:), ALLOCATABLE :: ba_ifilt ! filtered ionospheric ba
102 REAL(wp), DIMENSION(:), ALLOCATABLE :: c_nfilt ! error covariance ba_nfilt
103 REAL(wp), DIMENSION(:), ALLOCATABLE :: c_ifilt ! error covariance ba_ifilt
104 REAL(wp), DIMENSION(:), ALLOCATABLE :: ba_LI ! ionospheric bangle ba_LI
105 REAL(wp), DIMENSION(:), ALLOCATABLE :: err_LI ! error covariance ba_LI
106 REAL(wp), DIMENSION(:), ALLOCATABLE :: err_LC ! error covariance ba_LC
107 REAL(wp), DIMENSION(:), ALLOCATABLE :: WLCH ! weight of the linear combination !HGL!
108 LOGICAL, DIMENSION(:), ALLOCATABLE :: m_diff ! array mask
109
110 REAL(wp), DIMENSION(2) :: cin ! covariance of ionospheric noise
111 REAL(wp), DIMENSION(2) :: cis ! covariance of ionospheric signal
112 REAL(wp) :: sih ! covariance of L1 ionospheric signal
113 REAL(wp) :: sem ! mean square of relative model error
114 REAL(wp), DIMENSION(2,2) :: K ! system matrix
115 REAL(wp), DIMENSION(2,2) :: cs ! inverse signal covariance
116 REAL(wp), DIMENSION(2,2) :: cn ! inverse noie covariance
117 REAL(wp), DIMENSION(2,2) :: KC ! K^T*CN
118 REAL(wp), DIMENSION(2,2) :: KCK ! K^T*CN*K
119 REAL(wp), DIMENSION(2,2) :: KK ! K^T*CN*K + CS
120 REAL(wp), DIMENSION(2,2) :: KI ! (K^T*CN*K + CS)^(-1)
121 REAL(wp), DIMENSION(2,2) :: KQI(2,2) ! (K^T*CN*K + CS)^(-1)*K^T*CN !HGL!
122 REAL(wp), DIMENSION(2) :: ba_12 ! L1 and L2 bending angle
123 REAL(wp), DIMENSION(2) :: ba_ni ! neut and ion bending angle
124
125 REAL(wp) :: ba_diff0, ba_diffS
126 REAL(wp) :: impact_lt
127 REAL(wp) :: pmin, pmax
128 INTEGER :: iil, iiu, i_str, i_ltr
129 INTEGER :: i, n_obs, nh
130 INTEGER :: whi, wei
131 INTEGER :: w_smooth
132
133!-------------------------------------------------------------------------------
134! 2. Useful variables
135!-------------------------------------------------------------------------------
136
137 n_obs = SIZE(impact_L1)
138 pmax = MAXVAL(impact_L1)
139 pmin = MINVAL(impact_L1)
140
141 ALLOCATE(ba_LI(n_obs))
142 ALLOCATE(err_LI(n_obs))
143 ALLOCATE(err_LC(n_obs))
144
145! IMPLEMENTATION FROM 'INVERT'
146 IF (config%dpi < 50.) THEN
147 w_smooth = Ceiling(config%f_width*(config%npoints-1) / &
148 Abs(config%Pmax-config%Pmin))
149 ELSE
150! IMPLEMENTATION FROM 'OCC'
151 w_smooth = CEILING(config%f_width/config%dpi)
152 ENDIF
153
154!-------------------------------------------------------------------------------
155! 3. Array allocation
156!-------------------------------------------------------------------------------
157
158 nh = CEILING( ( pmax - pmin ) / config%delta_p ) + 1
159
160 ALLOCATE(impact_LH(nh))
161 ALLOCATE(ba_LMH(nh))
162 ALLOCATE(ba_is(nh))
163 ALLOCATE(ba_ion(nh))
164 ALLOCATE(ba_neut(nh))
165 ALLOCATE(ba_ifilt(nh))
166 ALLOCATE(ba_nfilt(nh))
167 ALLOCATE(c_ifilt(nh))
168 ALLOCATE(c_nfilt(nh))
169 ALLOCATE(ba_diff(2,nh))
170 ALLOCATE(ba_low(2,nh))
171 ALLOCATE(WLCH(nh)) !HGL!
172
173!-------------------------------------------------------------------------------
174! 4. Interpolation to homogeneous grid
175!-------------------------------------------------------------------------------
176
177 DO i=1,nh
178 impact_LH(i) = pmin + (i-1.0_wp)*( pmax - pmin )/(nh-1.0_wp)
179 ENDDO
180
181 ALLOCATE(ba_L1H(nh))
182 ALLOCATE(ba_L2H(nh))
183
184 CALL ropp_pp_interpol(impact_L1, impact_LH, ba_L1, ba_L1H)
185 CALL ropp_pp_interpol(impact_L2, impact_LH, ba_L2, ba_L2H)
186 CALL ropp_pp_interpol(impact_LM, impact_LH, ba_LM, ba_LMH)
187
188 ba_diff(1,:) = ba_L1H(:) - ba_LMH(:)
189 ba_diff(2,:) = ba_L2H(:) - ba_LMH(:)
190
191 DEALLOCATE(ba_L1H)
192 DEALLOCATE(ba_L2H)
193
194!-------------------------------------------------------------------------------
195! 5. Calculate smoothing windows
196!-------------------------------------------------------------------------------
197
198 whi = CEILING(w_smooth * (nh-1.0_wp)/( n_obs - 1.0_wp))
199 wei = CEILING(config%s_smooth * (nh-1.0_wp)/( pmax - pmin))
200
201!-------------------------------------------------------------------------------
202! 6. Region separation
203!-------------------------------------------------------------------------------
204
205 ! 6.1 Stratospheric upper limit
206
207 i_str = SUM(MINLOC(impact_LH(:), impact_LH(:)-config%r_curve > config%z_str))
208
209 ! 6.2 Initial estimation of upper limit of lower-tropospheric L2 perturbations
210
211 i_ltr = SUM(MINLOC(impact_LH(:), impact_LH(:)-config%r_curve > config%z_ltr))
212
213 ! 6.3 Lower limit of ionospheric signal/noise
214
215 iil = SUM(MINLOC(impact_LH(:), impact_LH(:)-config%r_curve > config%z_ion))
216 IF(iil < 1 .OR. iil > nh) iil = 0
217
218 ! 6.4 Dynamic upper limit for ionospheric noise estimate
219
220 IF(iil == 0) THEN
221 iiu = nh
222 ELSE
223 ba_diff0 = SUM(ba_diff(2,i_ltr:iil)-ba_diff(1,i_ltr:iil))/(iil-i_ltr+1.0_wp)
224 ba_diffS = SQRT(SUM((ba_diff(2,i_ltr:iil) - ba_diff(1,i_ltr:iil) - &
225 ba_diff0)**2) / (iil-i_ltr+1.0_wp))
226
227 ALLOCATE(m_diff(nh))
228 m_diff = ABS( ba_diff(2,:) - ba_diff(1,:) - ba_diff0) > 6.0_wp*ba_diffS
229 IF (ANY(m_diff(iil:nh) .AND. &
230 impact_LH(iil:nh)-config%r_curve > 75000.0_wp)) THEN
231 iiu = SUM(MaxLoc(impact_LH(:), impact_LH(:)-config%r_curve < 70000.0_wp))
232 ELSE IF (ANY(m_diff(iil:nh) .AND. &
233 impact_LH(iil:nh)-config%r_curve > 70000.0_wp)) THEN
234 iiu = SUM(MaxLoc(impact_LH(:), impact_LH(:)-config%r_curve < 65000.0_wp))
235 ELSE IF (ANY(m_diff(iil:nh) .AND. &
236 impact_LH(iil:nh)-config%r_curve > 65000.0_wp)) THEN
237 iiu = SUM(MaxLoc(impact_LH(:), impact_LH(:)-config%r_curve < 60000.0_wp))
238 ELSE
239 iiu = SUM(MaxLoc(impact_LH(:), impact_LH(:)-config%r_curve < 80000.0_wp))
240 ENDIF
241 DEALLOCATE(m_diff)
242 ENDIF
243
244!-------------------------------------------------------------------------------
245! 7. Covariance estimation
246!-------------------------------------------------------------------------------
247
248 ! 7.1 Smoothing with external scale
249
250 SELECT CASE (config%filter_method)
251 CASE('optest')
252 CALL ropp_pp_filter(impact_LH(2) - impact_LH(1), ba_diff(:, i_ltr:nh), wei, &
253 config%n_smooth, ba_low(:, i_ltr:nh))
254 CASE('slpoly')
255 CALL ropp_pp_sliding_polynomial(impact_LH(:), ba_diff(:, i_ltr:nh), wei, &
256 config%np_smooth, ba_low(:, i_ltr:nh))
257 END SELECT
258
259 ! Calculate low-frequency component of ionospheric refraction for L1
260
261 ba_is(i_ltr:nh) = (ba_low(1,i_ltr:nh) - ba_low(2,i_ltr:nh)) * &
262 f_L2**2/(f_L2**2 - f_L1**2)
263
264 ! 7.2 Smoothing with ionospheric scale
265
266 SELECT CASE (config%filter_method)
267 CASE('optest')
268 CALL ropp_pp_filter(impact_LH(2) - impact_LH(1), ba_diff(:, i_ltr:nh), whi, &
269 config%n_smooth, ba_low(:, i_ltr:nh))
270 CASE('slpoly')
271 CALL ropp_pp_sliding_polynomial(impact_LH(:), ba_diff(:, i_ltr:nh), whi, &
272 config%np_smooth, ba_low(:, i_ltr:nh))
273 END SELECT
274
275 ! 7.3 Estimation of ionospheric noise covariance
276
277 ba_neut(i_ltr:nh) = ba_low(1,i_ltr:nh) - &
278 ba_low(2,i_ltr:nh) * (f_L2/f_L1)**2
279
280 IF (iil > 0) THEN
281 cin(1) = (SUM(ba_neut(iil:iiu)**2)/(iiu-iil+1))/2.0_wp
282 cin(2) = ((f_L1/f_L2)**4)*(SUM(ba_neut(iil:iiu)**2)/(iiu-iil+1))/2.0_wp
283 ELSE
284 cin(:) = 0.0_wp
285 ENDIF
286
287 ! 7.4 Estimation of ionospheric signal covariance
288
289 ba_ion(i_ltr:nh) = (ba_low(1,i_ltr:nh) + &
290 ba_low(2,i_ltr:nh) * (f_L2/f_L1)**2)/2.0_wp - &
291 ba_is(i_ltr:nh)
292
293 IF (iil > 0) THEN
294 cis(1) = SUM(ba_ion(iil:iiu)**2)/(iiu-iil+1)
295 cis(2) = ((f_L1/f_L2)**4)*(SUM(ba_ion(iil:iiu)**2)/(iiu-iil+1))
296 cis(:) = MAX(0.0_wp, cis(:) - cin(:)/2)
297 ELSE
298 cis(:) = 0.0_wp
299 ENDIF
300
301 ! 7.5 Estimation of neutral signal covariance
302
303 sem = SUM((ba_diff(1, i_ltr:i_str)/ba_LMH(i_ltr:i_str))**2)/(i_str-i_ltr+1)
304
305 ! 7.6 Final estimation of lower tropospheric height
306
307 ba_ion(:) = (ba_diff(1, :) - ba_diff(2, :)) * &
308 f_L2**2/(f_L2**2 - f_L1**2)
309
310 sih = SUM((ba_ion(i_ltr:i_str)-ba_is(i_ltr:i_str))**2)/(i_str-i_ltr+1)
311
312 i = SUM(MAXLOC(impact_LH(1:i_ltr), &
313 (ba_ion(1:i_ltr) - ba_is(i_ltr+wei))**2 > 100*sih))
314 IF ( i < 1 .OR. i > i_ltr ) THEN
315 i = 0
316 ENDIF
317
318 i = MAX(1, i)
319 impact_lt = impact_LH(i) + 1000.0_wp
320 i_ltr = SUM(MINLOC(impact_LH(:), impact_LH(:) >= impact_lt))
321
322!-------------------------------------------------------------------------------
323! 8. Statistical optimization
324!-------------------------------------------------------------------------------
325
326 ! 8.1 Smoothing with external scale
327
328 SELECT CASE(config%filter_method)
329 CASE('optest')
330 CALL ropp_pp_filter(impact_LH(2) - impact_LH(1), ba_diff(:, i_ltr:nh), wei, &
331 config%n_smooth, ba_low(:, i_ltr:nh))
332 CASE('slpoly')
333 CALL ropp_pp_sliding_polynomial(impact_LH(:), ba_diff(:, i_ltr:nh), wei, &
334 config%np_smooth, ba_low(:, i_ltr:nh))
335 END SELECT
336
337 ba_is(i_ltr:nh) = (ba_low(1,i_ltr:nh) - ba_low(2,i_ltr:nh)) * &
338 f_L2**2/(f_L2**2 - f_L1**2)
339
340 ! 8.2 Smoothing with ionospheric scale
341
342 SELECT CASE(config%filter_method)
343 CASE('optest')
344 CALL ropp_pp_filter(impact_LH(2) - impact_LH(1), ba_diff(:, i_ltr:nh), whi, &
345 config%n_smooth, ba_low(:, i_ltr:nh))
346 CASE('slpoly')
347 CALL ropp_pp_sliding_polynomial(impact_LH(:), ba_diff(:, i_ltr:nh), whi, &
348 config%np_smooth, ba_low(:, i_ltr:nh))
349 END SELECT
350
351 ! 8.3 Border elimination
352
353 i_ltr = i_ltr + wei
354 impact_lt = impact_LH(i_ltr)
355
356 ! 8.4 Matrix calculation
357
358 K(:,:) = 1.0_wp
359 K(2,2) = (f_L1/f_L2)**2
360
361 IF (iil > 0) THEN
362 cn(:,:) = Diag(1.0_wp/cin(:))
363 ELSE
364 cn(:,:) = Diag((/1.0_wp, 1.0_wp/))
365 ENDIF
366
367 KC(:,:) = MATMUL(TRANSPOSE(K), CN)
368 KCK(:,:) = MATMUL(KC, K)
369
370 ! 8.5 Optimal linear combination of L1 and L2 bending angles
371
372 DO i=i_ltr, nh
373
374 IF (iil > 0) THEN
375 cs(:,:) = Diag((/1.0_wp/(sem*ba_LMH(i)**2), 1.0_wp/cis(1) /))
376 ELSE
377 cs(:,:) = 0.0_wp
378 ENDIF
379
380 KK(:,:) = KCK(:,:) + CS(:,:)
381
382 CALL ropp_pp_invert_matrix(KK, KI)
383
384 ba_12 = (/ ba_low(1,i) - ba_is(i), &
385 ba_low(2,i) - ba_is(i)*(f_L1/f_L2)**2 /)
386
387 ba_ni = MATMUL(KI, MATMUL(KC, ba_12))
388
389 KQI = MATMUL(KI, KC) !HGL!
390 WLCH(i) = KQI(1,1) + KQI(1,2) !HGL!
391
392 ba_nfilt(i) = ba_ni(1) + ba_LMH(i)
393 ba_ifilt(i) = ba_ni(2) + ba_is(i)
394
395 c_nfilt(i) = KI(1,1)
396 c_ifilt(i) = KI(2,2)
397
398 ENDDO
399
400 ! 8.6 Interpolation
401
402 DO i=1,n_obs
403
404 IF (impact_L1(i) > impact_lt) THEN
405
406 CALL ropp_pp_interpol(impact_LH(i_ltr:nh),impact_L1(i), &
407 ba_ifilt(i_ltr:nh),ba_LI(i))
408 CALL ropp_pp_interpol(impact_LH(i_ltr:nh),impact_L1(i), &
409 ba_nfilt(i_ltr:nh),ba_LC(i))
410 CALL ropp_pp_interpol(impact_LH(i_ltr:nh),impact_L1(i), &
411 c_ifilt(i_ltr:nh),err_LI(i))
412 CALL ropp_pp_interpol(impact_LH(i_ltr:nh),impact_L1(i), &
413 c_nfilt(i_ltr:nh),err_LC(i))
414 IF (PRESENT(WLC)) THEN !HGL!
415 CALL ropp_pp_interpol(impact_LH(i_ltr:nh), impact_L1(i), WLCH(i_ltr:nh), WLC(i)) !HGL!
416 ENDIF !HGL!
417
418 ELSE
419
420 ba_LI(i) = ba_is(i_ltr)
421 ba_LC(i) = ba_L1(i) - ba_LI(i)
422
423 err_LI(i) = c_ifilt(i_ltr)
424 err_LC(i) = c_nfilt(i_ltr)
425
426 IF (PRESENT(WLC)) THEN !HGL!
427 WLC(i) = 1.0 !HGL!
428 ENDIF !HGL!
429
430 ENDIF
431
432 ENDDO
433
434 impact_LC(:) = impact_L1(:)
435
436!-------------------------------------------------------------------------------
437! 9. Linear combination plus statistical optimization (if required)
438!-------------------------------------------------------------------------------
439
440 SELECT CASE(config%so_method)
441 CASE ( 'lcso' )
442 cin(1) = cin(1)*2.0_wp*f_L1**4 / (f_L1**2 - f_L2**2)**2
443 ba_LC(:) = ba_L1(:) - ba_LI(:)
444 ba_LC(:) = ba_LM(:) + (ba_LC(:)-ba_LM(:))*sem*ba_LM(:)**2 / &
445 (sem*ba_LM(:)**2 + cin(1))
446
447 END SELECT
448
449!-------------------------------------------------------------------------------
450! 10. Output diagnostics
451!-------------------------------------------------------------------------------
452
453 diag_out%sq = 100.0_wp * MAXVAL(SQRT(err_LC(:))/ba_LC(:)) ! 'Badness score'
454 ALLOCATE(diag_out%ba_ion(n_obs))
455 diag_out%ba_ion(:) = ba_LI(:) ! Ionospheric bending angle
456 ALLOCATE(diag_out%err_ion(n_obs))
457 diag_out%err_ion(:) = err_LI(:) ! Error covariance ionospheric bending
458 ALLOCATE(diag_out%err_neut(n_obs))
459 diag_out%err_neut(:) = err_LC(:) ! Error covariance neutral bending
460
461!-------------------------------------------------------------------------------
462! 11. Clean up
463!-------------------------------------------------------------------------------
464
465 DEALLOCATE(impact_LH)
466 DEALLOCATE(ba_LMH)
467 DEALLOCATE(ba_is)
468 DEALLOCATE(ba_ion)
469 DEALLOCATE(ba_neut)
470 DEALLOCATE(ba_ifilt)
471 DEALLOCATE(ba_nfilt)
472 DEALLOCATE(c_ifilt)
473 DEALLOCATE(c_nfilt)
474 DEALLOCATE(ba_diff)
475 DEALLOCATE(ba_low)
476
477 DEALLOCATE(ba_LI)
478 DEALLOCATE(err_LI)
479 DEALLOCATE(err_LC)
480
481 DEALLOCATE(WLCH) !HGL!
482
483CONTAINS
484
485!-------------------------------------------------------------------------------
486! 11. Generation of diagonal matrix
487!-------------------------------------------------------------------------------
488
489 FUNCTION diag(A)
490
491 IMPLICIT NONE
492
493 REAL(wp), DIMENSION(:), INTENT(in) :: A ! Array of diagonal elements
494 REAL(wp), DIMENSION(SIZE(A),SIZE(A)) :: diag ! Diagonal matrix
495
496 INTEGER :: i ! Diagonal index
497
498 diag(:,:) = 0.0_wp
499
500 DO i=1,SIZE(A)
501 diag(i,i) = A(i)
502 ENDDO
503
504 END FUNCTION diag
505
506END SUBROUTINE ropp_pp_ionospheric_correction