Ticket #545: ropp_fm_roprof2obs.f90_bangle

File ropp_fm_roprof2obs.f90_bangle, 10.6 KB (added by Ian Culverwell, 6 years ago)

ropp_fm_roprof2obs.f90_bangle

Line 
1SUBROUTINE ropp_fm_roprof2obs1dbangle(ro_data, y)
2
3! 1.1 Declarations
4! ----------------
5
6 USE typesizes, ONLY: wp => EightByteReal
7 USE ropp_utils
8 USE ropp_io
9 USE ropp_io_types, ONLY: ROprof
10 USE ropp_fm
11 USE ropp_fm_types, ONLY: Obs1dBangle
12 USE ropp_fm_copy, not_this => ropp_fm_roprof2obs1dbangle
13 USE geodesy
14
15 IMPLICIT NONE
16
17 TYPE(ROprof), INTENT(in) :: ro_data ! RO data structure
18 TYPE(Obs1dBangle), INTENT(inout) :: y ! Bending angle structure
19
20 INTEGER :: i, n, ipn
21 INTEGER, DIMENSION(8) :: DT8
22
23 CHARACTER(len = 256) :: routine
24 CHARACTER(len = 10) :: err_val
25 CHARACTER(LEN=6) :: sr_leo
26
27! 1.2 Error handling
28! ------------------
29
30 CALL message_get_routine(routine)
31 CALL message_set_routine('ropp_fm_roprof2obs (1D bending angles)')
32
33 y%obs_ok = .TRUE.
34
35! 1.3 Check and copy geolocation and time
36! ---------------------------------------
37
38 IF (isinrange(ro_data%georef%lon, ro_data%georef%range%lon)) THEN
39 y%lon = ro_data%georef%lon
40 ELSE
41 WRITE(err_val, '(e8.1)') ro_data%georef%lon
42 CALL message(msg_warn, &
43 "Longitude data for observations out of range or missing. " // &
44 "(longitude value = " // TRIM(err_val) // ")")
45 CALL message(msg_warn, "Check input data and valid_range attributes")
46 CALL message(msg_warn, "Set status flag obs%ok to FALSE")
47 y%obs_ok = .FALSE.
48 ENDIF
49
50 IF (isinrange(ro_data%georef%lat, ro_data%georef%range%lat)) THEN
51 y%lat = ro_data%georef%lat
52 ELSE
53 WRITE(err_val, '(e8.1)') ro_data%georef%lat
54 CALL message(msg_warn, &
55 "Latitude data for observations out of range or missing. " // &
56 "(latitude value = " // TRIM(err_val) // ")")
57 CALL message(msg_warn, "Check input data and valid_range attributes")
58 CALL message(msg_warn, "Set status flag obs%ok to FALSE")
59 y%obs_ok = .FALSE.
60 ENDIF
61
62 IF (isinrange(ro_data%georef%azimuth, ro_data%georef%range%azimuth)) THEN
63 y%azimuth = ro_data%georef%azimuth
64 ELSE
65 WRITE(err_val, '(e8.1)') ro_data%georef%azimuth
66 CALL message(msg_warn, "Azimuth (of " // TRIM(err_val) // " deg) out of range or missing.")
67 CALL message(msg_warn, "Check input data and valid_range attributes.")
68 CALL message(msg_warn, "Continuing, as azimuth not needed for 1D bending angles.")
69 ENDIF
70
71 IF ( isinrange(ro_data%DTocc%year, ro_data%DTocc%range%year) .AND. &
72 isinrange(ro_data%DTocc%month, ro_data%DTocc%range%month) .AND. &
73 isinrange(ro_data%DTocc%day, ro_data%DTocc%range%day) .AND. &
74 isinrange(ro_data%DTocc%hour, ro_data%DTocc%range%hour) .AND. &
75 isinrange(ro_data%DTocc%minute, ro_data%DTocc%range%minute) .AND. &
76 isinrange(ro_data%DTocc%second, ro_data%DTocc%range%second) ) THEN
77 DT8 = (/ro_data%DTocc%year, ro_data%DTocc%month, &
78 ro_data%DTocc%day, 0, &
79 ro_data%DTocc%hour, ro_data%DTocc%minute, &
80 ro_data%DTocc%second, ro_data%DTocc%msec/)
81 CALL TimeSince ( DT8, y%time, 1, Base="JS2000" )
82 ELSE
83 CALL message(msg_warn, &
84 "Time data for observations out of range or missing.")
85 CALL message(msg_warn, "Check input data and valid_range attributes")
86 CALL message(msg_warn, "Set status flag obs%ok to FALSE")
87 y%obs_ok = .FALSE.
88 ENDIF
89
90! 1.4 Check that profiles are increasing in height - 1st element towards surface
91! ------------------------------------------------------------------------------
92
93 CALL ropp_io_ascend(ro_data)
94
95! 1.5 Check and copy observation data
96! -----------------------------------
97
98 IF (ro_data%Lev1b%Npoints == 0) THEN
99 CALL message(msg_warn, &
100 "RO data has no level 1b (bending angle) parameters.")
101 CALL message(msg_warn, "Check input data and valid_range attributes")
102 CALL message(msg_warn, "Set status flag obs%ok to FALSE")
103 y%obs_ok = .FALSE.
104 RETURN
105 ENDIF
106
107 n = ro_data%Lev1b%Npoints
108
109 IF (.NOT. ro_data%Lev2c%direct_ion) THEN
110
111 ALLOCATE(y%impact(n), y%bangle(n), y%weights(n))
112 y%impact = ro_data%Lev1b%impact(1:n)
113 y%bangle = ro_data%Lev1b%bangle(1:n)
114 y%weights = 1.0_wp
115 y%n_L1 = n
116
117 ELSE
118
119 ALLOCATE(y%impact(2*n), y%bangle(2*n), y%weights(2*n))
120
121 IF (ALL(ro_data%Lev1b%impact_L1 < ropp_MDTV)) THEN ! Use neutral IP if no L1 IPs
122 y%impact(1:n) = ro_data%Lev1b%impact(1:n)
123 ELSE
124 y%impact(1:n) = ro_data%Lev1b%impact_L1(1:n)
125 END IF
126
127 IF (ALL(ro_data%Lev1b%impact_L2 < ropp_MDTV)) THEN ! Use neutral IP if no L2 IPs
128 y%impact(n+1:2*n) = ro_data%Lev1b%impact(1:n)
129 ELSE
130 y%impact(n+1:2*n) = ro_data%Lev1b%impact_L2(1:n)
131 END IF
132
133 y%bangle(1:n) = ro_data%Lev1b%bangle_L1(1:n)
134 y%bangle(n+1:2*n) = ro_data%Lev1b%bangle_L2(1:n)
135 y%weights = 1.0_wp
136 y%n_L1 = n
137
138 END IF
139
140 y%g_sfc = gravity(ro_data%georef%lat, 0.0_wp)
141
142 y%r_earth = R_eff(ro_data%georef%lat)
143
144 IF (isinrange(ro_data%georef%roc, ro_data%georef%range%roc)) THEN
145 y%r_curve = ro_data%georef%roc
146 ELSE
147 CALL message(msg_warn, "Radius of curvature out of range or missing.")
148 CALL message(msg_warn, "Check input data and valid_range attributes")
149 CALL message(msg_warn, "Set status flag obs%ok to FALSE")
150 y%obs_ok = .FALSE.
151 ENDIF
152
153 IF ( ro_data%Lev2c%direct_ion ) THEN ! Need the radius of the LEO orbit for single frequency work
154
155 IF ( ANY ( ABS(ro_data%georef%leo_pod%pos - ropp_MDFV) < 1.0_wp ) .OR. &
156 ANY ( ABS(ro_data%georef%r_coc - ropp_MDFV) < 1.0_wp ) ) THEN
157
158 CALL message (msg_warn, 'Missing LEO position or centre of curvature ... ' // &
159 'will try to estimate r_LEO from satellite ID.')
160
161 SELECT CASE (ro_data%leo_id(1:2)) ! Data from https://www.wmo-sat.info/oscar/
162
163 CASE ('C0')
164 CALL message (msg_info, 'Assuming LEO height of 800 km (COSMIC-1 profile)')
165 y%r_leo = y%r_curve + 800.E3_wp
166
167 CASE ('C2')
168 CALL message (msg_info, 'Assuming LEO height of 520 km (COSMIC-2a profile)')
169 y%r_leo = y%r_curve + 520.E3_wp
170
171 CASE ('CH')
172 CALL message (msg_info, 'Assuming LEO height of 470 km (CHAMP profile)')
173 y%r_leo = y%r_curve + 470.E3_wp
174
175 CASE ('FY')
176 CALL message (msg_info, 'Assuming LEO height of 840 km (FY-3C/D/... profile)')
177 y%r_leo = y%r_curve + 840.E3_wp
178
179 CASE ('GR')
180 SELECT CASE (ro_data%leo_id(1:4))
181 CASE ('GRAA', 'GRAB')
182 CALL message (msg_info, 'Assuming LEO height of 340 km (GRACE-A/B profile)')
183 y%r_leo = y%r_curve + 340.E3_wp
184 CASE ('GRAC', 'GRAD')
185 CALL message (msg_info, 'Assuming LEO height of 490 km (GRACE-FO profile)')
186 y%r_leo = y%r_curve + 490.E3_wp
187 CASE DEFAULT
188 CALL message(msg_warn, 'Unrecognised LEO satellite ID (' // ro_data%leo_id // &
189 ') ... r_LEO cannot be estimated.')
190 CALL message(msg_warn, 'Check that LEO POD and CoC positions, or LEO ID, are valid.')
191 CALL message(msg_warn, 'Setting status flag obs%ok to FALSE')
192 y%r_leo = ropp_MDFV
193 y%obs_ok = .FALSE.
194 END SELECT
195
196 CASE ('KO')
197 CALL message (msg_info, 'Assuming LEO height of 550 km (KOMPSAT-5 profile)')
198 y%r_leo = y%r_curve + 550.E3_wp
199
200 CASE ('ME')
201 CALL message (msg_info, 'Assuming LEO height of 830 km (Metop-A/B/C profile)')
202 y%r_leo = y%r_curve + 830.E3_wp
203
204 CASE ('MG')
205 CALL message (msg_info, 'Assuming LEO height of 870 km (Megha-Tropiques profile)')
206 y%r_leo = y%r_curve + 870.E3_wp
207
208 CASE ('PA')
209 CALL message (msg_info, 'Assuming LEO height of 510 km (PAZ profile)')
210 y%r_leo = y%r_curve + 510.E3_wp
211
212 CASE DEFAULT
213 CALL message(msg_warn, 'Unrecognised LEO satellite ID (' // ro_data%leo_id // &
214 ') ... r_LEO cannot be estimated.')
215 CALL message(msg_warn, 'Check that LEO POD and CoC positions, or LEO ID, are valid.')
216 CALL message(msg_warn, 'Setting status flag obs%ok to FALSE')
217 y%r_leo = ropp_MDFV
218 y%obs_ok = .FALSE.
219
220 END SELECT
221
222 ELSE
223
224 y%r_leo = SQRT ( SUM ( (ro_data%georef%leo_pod%pos - ro_data%georef%r_coc)**2 ) )
225 WRITE (sr_leo, FMT='(F6.1)') (y%r_leo - y%r_curve) * 1.0E-3_wp
226 CALL message (msg_info, 'Using a LEO altitude of ' // sr_leo // ' km.')
227
228 END IF
229
230 ELSE
231
232 y%r_leo = ropp_MDFV ! For safety's sake
233
234 END IF ! Single frequency work
235
236 IF (isinrange(ro_data%georef%undulation, ro_data%georef%range%undulation)) THEN
237 y%undulation = ro_data%georef%undulation
238 ELSE
239 y%undulation = ropp_MDFV
240 CALL message(msg_warn, "Undulation out of range or missing.")
241 CALL message(msg_warn, "Check input data and valid_range attributes")
242 CALL message(msg_warn, "Set status flag obs%ok to FALSE")
243 y%obs_ok = .FALSE.
244 END IF
245
246! 1.6 Check and copy sigmas to diagonal error covariance
247! ------------------------------------------------------
248
249 y%cov_ok = .TRUE.
250
251 IF (ASSOCIATED(y%cov%d)) DEALLOCATE(y%cov%d)
252
253 IF (.NOT. ro_data%Lev2c%direct_ion) THEN
254
255 CALL callocate(y%cov%d, n*(n+1)/2)
256
257 DO i = 1, n
258 IF (ro_data%Lev1b%bangle(i) > 0.0_wp) THEN
259 IF (ro_data%Lev1b%bangle_sigma(i) > 0.0_wp) THEN
260 ! matrix_pp type, uplo = 'U'
261 y%cov%d(i + i*(i-1)/2) = ro_data%Lev1b%bangle_sigma(i)**2
262 ELSE
263 y%cov_ok = .FALSE.
264 END IF
265 ELSE
266 y%cov%d(i + i*(i-1)/2) = 0.0003_wp ! = 1 deg**2
267 ENDIF
268 END DO
269
270 ELSE
271
272 CALL callocate(y%cov%d, n*(2*n+1))
273
274 DO i = 1, n
275
276 IF (ro_data%Lev1b%bangle_L1(i) > 0.0_wp) THEN
277 IF (ro_data%Lev1b%bangle_L1_sigma(i) > 0.0_wp) THEN
278 ! matrix_pp type, uplo = 'U'
279 y%cov%d(i + i*(i-1)/2) = ro_data%Lev1b%bangle_L1_sigma(i)**2
280 ELSE
281 y%cov_ok = .FALSE.
282 END IF
283 ELSE
284 y%cov%d(i + i*(i-1)/2) = 0.0003_wp ! = 1 deg**2
285 END IF
286
287 ipn = i + n
288 IF (ro_data%Lev1b%bangle_L2(i) > 0.0_wp) THEN
289 IF (ro_data%Lev1b%bangle_L2_sigma(i) > 0.0_wp) THEN
290 ! matrix_pp type, uplo = 'U'
291 y%cov%d(ipn + ipn*(ipn-1)/2) = ro_data%Lev1b%bangle_L2_sigma(i)**2
292 ELSE
293 y%cov_ok = .FALSE.
294 END IF
295 ELSE
296 y%cov%d(ipn + ipn*(ipn-1)/2) = 0.0003_wp ! = 1 deg**2
297 END IF
298
299 END DO
300
301 END IF
302
303 IF (ASSOCIATED(y%cov%e)) DEALLOCATE(y%cov%e)
304 IF (ASSOCIATED(y%cov%f)) DEALLOCATE(y%cov%f)
305 IF (ASSOCIATED(y%cov%s)) DEALLOCATE(y%cov%s)
306
307 y%cov%fact_chol = .FALSE.
308 y%cov%equi_chol = 'N'
309
310! 1.7 Clean up
311! ------------
312
313 CALL message_set_routine(routine)
314
315END SUBROUTINE ropp_fm_roprof2obs1dbangle