1module bgckCSR_mod
2 ! MODULE bgckCSR_mod (prefix='csrbg' category='1. High-level functionality')
3 !
4 !:Purpose: To perform CSR (clear-sky radiance from geostationary
5 ! satellite platforms) data background check.
6 !
7 use midasMpi_mod
8 use utilities_mod
9 use obsSpaceData_mod
10 use tovsNL_mod
11 use obsErrors_mod
12
13 implicit none
14 save
15 private
16
17 real, parameter :: csrbg_realMissing=-99.
18
19 ! Public functions/subroutines
20
21 public :: csrbg_bgCheckCSR
22
23 real, parameter :: csrbg_ompThreshold = 4.2
24 integer, parameter :: maxNumsat = 20 ! nb max de satellites
25 integer, parameter :: maxNumchan = 15 ! nb max de canaux
26
27 ! namelist variables
28 character (len=9) :: burpSatName(maxNumSat) ! list of platform names to treat (BURP file station id)
29 integer :: satCloudCoverLimit(maxNumSat,maxNumchan) ! maximum limit of cloud cover (careful this is an integer!)
30
31 namelist /namcsr/burpSatName,satCloudCoverLimit
32
33contains
34
35 !----------------------------------------------------------------------------------------
36 ! csrbg_init
37 !----------------------------------------------------------------------------------------
38 subroutine csrbg_init()
39 !
40 !:Purpose: This subroutine reads the namelist section NAMCSR
41 ! for the module.
42 implicit none
43
44 ! Locals:
45 integer :: nulnam, ierr
46 integer, external :: fnom, fclos
47
48 nulnam = 0
49 ierr = fnom(nulnam, './flnml','FTN+SEQ+R/O', 0)
50 read(nulnam, nml=namcsr, iostat=ierr)
51 if (ierr /= 0) call utl_abort('csrbg_init: Error reading namelist')
52 if (mmpi_myid == 0) write(*, nml=namcsr)
53 ierr = fclos(nulnam)
54
55 end subroutine csrbg_init
56
57 !----------------------------------------------------------------------------------------
58 ! csrbg_bgCheckCSR
59 !----------------------------------------------------------------------------------------
60 subroutine csrbg_bgCheckCSR (obsSpaceData)
61 !
62 !: Purpose: Effectuer le controle que qualite sur les donnees CSR.
63 ! Modifier les marqueurs de donnees selon le type de rejet.
64 !
65 implicit none
66
67 ! Arguments:
68 type(struct_obs), intent(inout) :: obsSpaceData ! obspaceData Object
69
70 ! Locals:
71 real , allocatable :: obsTb(:) ! brightness temperature (btyp=9248/9264,ele=12163)
72 real , allocatable :: ompTb(:) ! OMP values
73 real , allocatable :: satZenithAngle(:) ! satellite zenith angle (btyp=3072,ele=7024)
74 integer, allocatable :: obsFlags(:) ! data flags
75 integer, allocatable :: cloudAmount(:) ! data flags
76 integer, allocatable :: obsChannels(:) ! Tb Channels
77 integer, allocatable :: obsDate(:) ! date YYYYMMDD
78 integer, allocatable :: obsHour(:) ! Hour HHMM
79 integer :: sensorIndex ! find tvs_sensor index corresponding to current obs
80 character(len=9) :: burpFileSatId ! Platform Name
81 ! Data derived from elements read in obspace data
82 logical, allocatable :: maxAngleReached(:) ! satellite angle exceed max angle at obs
83 logical, allocatable :: topographicData(:) ! data flagged as topo data
84 logical, allocatable :: nonCorrectedData(:) ! data non corrected by bias corr
85 logical, allocatable :: isTbPresent(:) ! non missing data
86 logical, allocatable :: isClearSky(:) ! clear sky obs
87 logical, allocatable :: strayLight(:) !
88 logical, allocatable :: goesMidi(:) ! goes noon
89 logical, allocatable :: isToAssim(:) ! is channel assimilable
90 logical, allocatable :: ompOutOfRange(:)
91 integer :: categorieRejet(7)
92 integer :: headerIndex
93 integer :: codtyp
94 logical :: csrDataPresent
95
96 categorieRejet(:) = 0
97 csrDataPresent = .false.
98 call obs_set_current_header_list(obsSpaceData,'TO')
99 HEADER0: do
100 headerIndex = obs_getHeaderIndex(obsSpaceData)
101 if (headerIndex < 0) exit HEADER0
102 codtyp = obs_headElem_i(obsSpaceData, OBS_ITY, headerIndex)
103 if (tvs_isIdBurpInst(codtyp,'radianceclear' )) csrDataPresent = .true.
104 end do HEADER0
105
106 if ( .not. csrDataPresent ) then
107 write(*,*) 'WARNING: WILL NOT RUN csrbg_bgCheckCSR since no CSR Data'
108 return
109 end if
110
111 call utl_tmg_start(116,'--BgckCSR')
112 write(*,*) ' CSRBG QC PROGRAM STARTS ....'
113 ! Read Namelist
114 call csrbg_init()
115
116 call obs_set_current_header_list(obsSpaceData,'TO')
117 HEADER: do
118 headerIndex = obs_getHeaderIndex(obsSpaceData)
119 if (headerIndex < 0) exit HEADER
120 codtyp = obs_headElem_i(obsSpaceData, OBS_ITY, headerIndex)
121 if ( .not. tvs_isIdBurpInst(codtyp,'radianceclear') ) then
122 write(*,*) 'WARNING: Observation with codtyp = ', codtyp, ' is not a CSR data'
123 cycle HEADER
124 end if
125
126 !###############################################################################
127 ! STEP 1) read obs from obsSpacedata to start QC !
128 !###############################################################################
129
130 call csrbg_readObsFromObsSpace(obsSpaceData, headerIndex, obsTb, ompTb, satZenithAngle, obsFlags, cloudAmount, &
131 obsChannels, sensorIndex, obsDate, obsHour, burpFileSatId, &
132 maxAngleReached, topographicData, nonCorrectedData, isTbPresent, &
133 isClearSky, strayLight, goesMidi, isToAssim, ompOutOfRange)
134
135 !###############################################################################
136 ! STEP 2) Controle de qualite des CSR. Data QC flags (obsFlags) are modified here!
137 !###############################################################################
138
139 call csrbg_csrCheckQc(obsFlags, categorierejet, maxAngleReached, topographicData, &
140 nonCorrectedData, isTbPresent, isClearSky, strayLight, goesMidi, isToAssim, &
141 ompOutOfRange)
142
143 !###############################################################################
144 ! STEP 3) Update Flags and obs in obsspace data
145 !###############################################################################
146
147 call csrbg_updateObsSpaceAfterQc(obsSpaceData, obsFlags, headerIndex, sensorIndex)
148
149 end do HEADER
150 write(*,*) "Nombre de donnees rejetees"
151 write(*,*) "Attention, une donnee peut etre rejetee pour plusieurs raisons"
152 write(*,*) "Maxangle, straylight ou goesmid " , categorieRejet(1)
153 write(*,*) "Topographie " , categorieRejet(2)
154 write(*,*) "TB sans correction de biais " , categorieRejet(3)
155 write(*,*) "Pas TB " , categorieRejet(4)
156 write(*,*) "std(O-P)*sigma trop grand " , categorieRejet(5)
157 write(*,*) "canal non assimile " , categorieRejet(6)
158 write(*,*) "Ciel non clair " , categorieRejet(7)
159 write(*,*) "*******"
160
161 call utl_tmg_stop(116)
162
163 end subroutine csrbg_bgCheckCSR
164
165 !--------------------------------------------------------------------------
166 ! csrbg_readObsFromObsSpace
167 !--------------------------------------------------------------------------
168 subroutine csrbg_readObsFromObsSpace(obsSpaceData, headerIndex, obsTb, ompTb, satZenithAngle, obsFlags, cloudAmount, &
169 obsChannels, sensorIndex, obsDate, obsHour, burpFileSatId, &
170 maxAngleReached, topographicData, nonCorrectedData, isTbPresent, &
171 isClearSky, strayLight, goesMidi, isToAssim, ompOutOfRange)
172 !
173 !:Purpose: copy headers and bodies from obsSpaceData object to arrays
174 ! compute some parameters from the read variables
175 !
176 implicit None
177
178 ! Arguments:
179 integer , intent(in) :: headerIndex ! current header index
180 real , allocatable, intent(out) :: obsTb(:) ! brightness temperature (btyp=9248/9264,ele=12163)
181 real , allocatable, intent(out) :: ompTb(:) ! OMP values
182 real , allocatable, intent(out) :: satZenithAngle(:) ! satellite zenith angle (btyp=3072,ele=7024)
183 integer, allocatable, intent(out) :: obsFlags(:) ! data flags
184 integer, allocatable, intent(out) :: cloudAmount(:) ! data flags
185 integer, allocatable, intent(out) :: obsChannels(:) ! Tb Channels
186 integer, allocatable, intent(out) :: obsDate(:) ! date YYYYMMDD
187 integer, allocatable, intent(out) :: obsHour(:) ! Hour HHMM
188 integer, intent(out) :: sensorIndex ! find tvs_sensor index corresponding to current obs
189 character(*), intent(out) :: burpFileSatId ! Platform Name
190 type(struct_obs), intent(inout) :: obsSpaceData ! obspaceData Object
191 logical, allocatable, intent(out) :: maxAngleReached(:) ! satellite angle exceed max angle at obs
192 logical, allocatable, intent(out) :: topographicData(:) ! data flagged as topo data
193 logical, allocatable, intent(out) :: nonCorrectedData(:) ! data non corrected by bias corr
194 logical, allocatable, intent(out) :: isTbPresent(:) ! non missing data
195 logical, allocatable, intent(out) :: isClearSky(:) ! clear sky obs
196 logical, allocatable, intent(out) :: strayLight(:) !
197 logical, allocatable, intent(out) :: goesMidi(:) ! goes noon
198 logical, allocatable, intent(out) :: isToAssim(:) ! is channel assimilable
199 logical, allocatable, intent(out) :: ompOutOfRange(:) !
200
201 ! Locals:
202 integer :: bodyIndex
203 integer :: bodyIndexbeg
204 integer :: headerCompt
205 integer :: currentChannelNumber
206 integer :: channelIndex
207 integer :: actualnumChannel
208 integer :: numObsToProcess
209 integer :: iplatform
210 integer :: instrum
211 integer :: isat, iplatf
212 integer :: instr
213 integer :: mois
214 integer :: jour
215 integer :: heure
216 integer :: midnight
217 integer :: dataIndex
218 integer :: indexSat
219 logical :: sensorIndexFound
220 logical :: indexSatFound
221 real :: errorThreshold
222
223 ! find tvs_sensor index corresponding to current obs
224
225 iplatf = obs_headElem_i( obsSpaceData, OBS_SAT, headerIndex )
226 instr = obs_headElem_i( obsSpaceData, OBS_INS, headerIndex )
227
228 call tvs_mapSat( iplatf, iplatform, isat )
229 call tvs_mapInstrum( instr, instrum )
230
231 sensorIndexFound = .false.
232 do sensorIndex =1, tvs_nsensors
233 if ( iplatform == tvs_platforms(sensorIndex) .and. &
234 isat == tvs_satellites(sensorIndex) .and. &
235 instrum == tvs_instruments(sensorIndex) ) then
236 sensorIndexFound = .true.
237 exit
238 end if
239 end do
240
241 if ( .not. sensorIndexFound ) call utl_abort('csrbg_readObsFromObsSpace: sensor Index not found')
242
243 ! find actual Number of channels
244 actualNumChannel = tvs_coefs(sensorIndex)%coef%fmv_ori_nchn
245
246 headerCompt = 1
247 numObsToProcess = 1
248 ! Allocate Header elements
249 call utl_reAllocate(obsDate, numObsToProcess)
250 call utl_reAllocate(obsHour, numObsToProcess)
251 call utl_reAllocate(satZenithAngle, numObsToProcess)
252
253 ! Allocate Body elements
254 call utl_reAllocate(obsTb, numObsToProcess*actualNumChannel)
255 call utl_reAllocate(ompTb, numObsToProcess*actualNumChannel)
256 call utl_reAllocate(obsFlags, numObsToProcess*actualNumChannel)
257 call utl_reAllocate(cloudAmount, numObsToProcess*actualNumChannel)
258 call utl_reAllocate(obsChannels, numObsToProcess*actualNumChannel)
259
260 !initialization
261 obsTb(:) = csrbg_realMissing
262 ompTb(:) = csrbg_realMissing
263 obsChannels(:) = csrbg_realMissing
264 cloudAmount(:) = csrbg_realMissing
265 !
266 ! Header elements
267 !
268 burpFileSatId = obs_elem_c ( obsSpaceData, 'STID' , headerIndex )
269 satZenithAngle(headerCompt) = obs_headElem_r( obsSpaceData, OBS_SZA, headerIndex )
270 obsDate(headerCompt) = obs_headElem_i( obsSpaceData, OBS_DAT, headerIndex )
271 obsHour(headerCompt) = obs_headElem_i( obsSpaceData, OBS_ETM, headerIndex )
272 !
273 ! Body elements
274 !
275 bodyIndexbeg = obs_headElem_i( obsSpaceData, OBS_RLN, headerIndex )
276
277 BODY: do bodyIndex = bodyIndexbeg, bodyIndexbeg + obs_headElem_i( obsSpaceData, OBS_NLV, headerIndex ) - 1
278 currentChannelNumber = nint(obs_bodyElem_r( obsSpaceData, OBS_PPP, bodyIndex )) - &
279 tvs_channelOffset(sensorIndex)
280 obsTb(currentChannelNumber) = obs_bodyElem_r( obsSpaceData, OBS_VAR, bodyIndex )
281 ompTb(currentChannelNumber) = obs_bodyElem_r( obsSpaceData, OBS_OMP, bodyIndex )
282 obsFlags(currentChannelNumber) = obs_bodyElem_i( obsSpaceData, OBS_FLG, bodyIndex )
283 cloudAmount(currentChannelNumber) = obs_bodyElem_i( obsSpaceData, OBS_CLA, bodyIndex)
284
285 end do BODY
286
287 do channelIndex=1, actualNumChannel
288 obsChannels(channelIndex) = channelIndex+tvs_channelOffset(sensorIndex)
289 end do
290
291
292 ! Get index for satCloudCoverLimit for current current sat
293 indexSatFound = .false.
294 do indexSat = 1, maxNumSat
295 if (trim(burpSatName(indexSat)) == trim(burpFileSatId)) then
296 indexSatFound = .true.
297 exit
298 end if
299 end do
300 if ( .not. indexSatFound ) call utl_abort('csrbg_readObsFromObsSpace: Cloud Cover Limit Not &
301 Found for ' // trim(burpFileSatId))
302 ! compute data for QC
303 call utl_reAllocate(topographicData, numObsToProcess)
304 call utl_reAllocate(maxAngleReached, numObsToProcess)
305 call utl_reAllocate(strayLight, numObsToProcess)
306 call utl_reAllocate(goesMidi, numObsToProcess)
307 call utl_reAllocate(isTbPresent, numObsToProcess*actualNumChannel)
308 call utl_reAllocate(nonCorrectedData, numObsToProcess*actualNumChannel)
309 call utl_reAllocate(isClearSky, numObsToProcess*actualNumChannel)
310 call utl_reAllocate(ompOutOfRange, numObsToProcess*actualNumChannel)
311 call utl_reAllocate(isToAssim, numObsToProcess*actualNumChannel)
312
313 isTbPresent(:) = .false.
314 isClearSky(:) = .false.
315 nonCorrectedData = .false.
316
317 do channelIndex=1, numObsToProcess*actualNumChannel
318 if (obsTb(channelIndex) /= csrbg_realMissing) isTbPresent(channelIndex) = .true.
319 if (cloudAmount(channelIndex) /= csrbg_realMissing .and. &
320 cloudAmount(channelIndex) < satCloudCoverlimit(indexSat,channelIndex)) then
321 isClearSky(channelIndex) = .true.
322 end if
323 if (.not. btest(obsFlags(channelIndex), 6)) nonCorrectedData(channelIndex) = .true.
324 end do
325
326 topographicData(:) = .false.
327 maxAngleReached(:) = .false.
328 do dataIndex=1, numObsToProcess
329 if (btest(obsFlags(dataIndex), 18)) topographicData(dataIndex) = .true.
330 if (satZenithAngle(dataIndex) > 15250) maxAngleReached(dataIndex) = .true.
331 end do
332
333 ! check O-P of Tb and set assim flag for each channel of satellite numsat using stat_iutilst
334 ! If channel not found in stats file, NO ASSIMILATE flag is set and O-P check is skipped
335 ompOutOfRange(:) = .false.
336 isToAssim(:) = .false.
337 do channelIndex=1, numObsToProcess*actualNumChannel
338 if ( .not. oer_tovutil(obsChannels(channelIndex), sensorIndex) == 0) then
339 isToAssim(channelIndex) = .true.
340 end if
341 errorThreshold = oer_toverrst(obsChannels(channelIndex), sensorIndex)*csrbg_ompThreshold
342 if (abs(ompTb(channelIndex)) > errorThreshold) ompOutOfRange(channelIndex) = .true.
343 end do
344
345 strayLight(:) = .false.
346 goesMidi(:) = .false.
347 do dataIndex = 1, numObsToProcess
348 mois = (obsDate(dataIndex)/100) - (obsDate(dataIndex)/10000)*100
349 jour = (obsDate(dataIndex)/100) - (obsDate(dataIndex)/100)*100
350 heure = obsHour(dataIndex)/100
351 if (burpFileSatId == '^METSAT7') then
352 if ( heure == 21 ) straylight(dataIndex) = 1
353 if ( (heure == 20) .or. (heure == 22) ) then
354 if ( (jour + (mois-1)*30) > 20 .and. (jour + (mois-1)*30) < 135 ) straylight(dataIndex) = .true.
355 if ( (jour + (mois-1)*30) > 210 .and. (jour + (mois-1)*30) < 300 ) straylight(dataIndex) = .true.
356 end if
357 end if
358 if (burpFileSatId == '^GOES11' .or. burpFileSatId == '^GOES15') midnight = 9
359 if (burpFileSatId == '^GOES13' .or. burpFileSatId == '^GOES14') midnight = 5
360 if (burpFileSatId == '^GOES11' .or. burpFileSatId == '^GOES15' .or. &
361 burpFileSatId == '^GOES13' .or. burpFileSatId == '^GOES14') then
362 if ( (heure == (midnight - 1)) .or. (heure == midnight) ) goesMidi(dataIndex) = .true.
363 end if
364 end do
365
366 end subroutine csrbg_readObsFromObsSpace
367
368 !--------------------------------------------------------------------------
369 ! csrbg_csrCheckQc
370 !--------------------------------------------------------------------------
371 subroutine csrbg_csrCheckQc(obsFlags, categorieRejet, maxAngleReached, topographicData, &
372 nonCorrectedData, isTbPresent, isClearSky, strayLight, goesMidi, isToAssim, &
373 ompOutOfRange)
374 !
375 !:Purpose: Modify obsFlags
376 !
377 implicit None
378
379 ! Arguments:
380 integer, intent(inout) :: obsFlags(:) ! obs Flags to update
381 integer, intent(inout) :: categorieRejet(7) ! the 7 categories of rejections
382 logical, intent(in) :: maxAngleReached(:) ! satellite angle exceed max angle at obs
383 logical, intent(in) :: topographicData(:) ! data flagged as topo data
384 logical, intent(in) :: nonCorrectedData(:) ! data non corrected by bias corr
385 logical, intent(in) :: isTbPresent(:) ! non missing data
386 logical, intent(in) :: isClearSky(:) ! clear sky obs
387 logical, intent(in) :: strayLight(:) !
388 logical, intent(in) :: goesMidi(:) ! goes noon
389 logical, intent(in) :: isToAssim(:) ! is channel assimilable
390 logical, intent(in) :: ompOutOfRange(:) ! abs of omp greater than threshold
391
392 ! Locals:
393 integer :: currentChannelNumber
394 integer :: channelIndex
395 integer :: numObsToProcess
396 integer :: numData
397 integer :: dataIndex
398
399 numObsToProcess = 1
400 currentChannelNumber = size(obsFlags)/numObsToProcess
401
402 ! tests de controle de qualite sur les profils et canaux
403 numData = 0
404 do dataIndex = 1, numObsToProcess
405 do channelIndex = 1, currentChannelNumber
406 numData = numData+1
407 !Maxangle & Straylight & Goesmid : Bit 7 et 9 pour tous les canaux
408 if (maxAngleReached(dataIndex) == .true. .or. &
409 strayLight(dataIndex) == .true. .or. &
410 goesMidi(dataIndex) == .true. ) then
411 obsFlags(numData) = ibset(obsFlags(numData),7)
412 obsFlags(numData) = ibset(obsFlags(numData),9)
413 categorieRejet(1) = categorieRejet(1) + 1
414 end if
415
416 !Topo : bit 9 et 18 sont déjà là provenant du EnVar
417 if (topographicData(dataIndex) == .true.) then
418 obsFlags(numData) = ibset(obsFlags(numData),9)
419 obsFlags(numData) = ibset(obsFlags(numData),18)
420 categorieRejet(2) = categorieRejet(2) + 1
421 end if
422
423 !Pas corrige : bit 11
424 if (nonCorrectedData(numData) == .true.) then
425 obsFlags(numData) = ibset(obsFlags(numData),11)
426 categorieRejet(3) = categorieRejet(3) + 1
427 end if
428
429 !Tbyela: bit 7 et 9 pour les canaux concernés
430 if (isTbPresent(numData) == .false.) then
431 obsFlags(numData) = ibset(obsFlags(numData),7)
432 obsFlags(numData) = ibset(obsFlags(numData),9)
433 categorieRejet(4) = categorieRejet(4) + 1
434 end if
435
436 !omp Out of range: Bit 9 et 16 pour les canaux concernés
437 if (ompOutOfRange(numData) == .true.) then
438 obsFlags(numData) = ibset(obsFlags(numData),9)
439 obsFlags(numData) = ibset(obsFlags(numData),16)
440 categorieRejet(5) = categorieRejet(5) + 1
441 end if
442
443 !Assim: Bit 11 pour les canaux concernés
444 if (isToAssim(numData) == .false.) then
445 obsFlags(numData) = ibset(obsFlags(numData),11)
446 categorieRejet(6) = categorieRejet(6) + 1
447 end if
448
449 !Clearsky : Bit 7 & Bit 9
450 if (isClearSky(numData) == .false.) then
451 obsFlags(numData) = ibset(obsFlags(numData),7)
452 obsFlags(numData) = ibset(obsFlags(numData),9)
453 categorieRejet(7) = categorieRejet(7) + 1
454 end if
455 end do
456 end do
457
458 end subroutine csrbg_csrCheckQc
459
460 !--------------------------------------------------------------------------
461 ! csrbg_updateObsSpaceAfterQc
462 !--------------------------------------------------------------------------
463 subroutine csrbg_updateObsSpaceAfterQc(obsSpaceData, obsFlags, headerIndex, sensorIndex)
464 !
465 !:Purpose: Update obspacedata variables (obstTB and obs flags) after QC
466 !
467 implicit None
468
469 ! Arguments:
470 type(struct_obs), intent(inout) :: obsSpaceData ! obspaceData Object
471 integer, intent(in) :: obsFlags(:) ! data flags
472 integer, intent(in) :: sensorIndex ! sensor Index
473 integer, intent(in) :: headerIndex ! sensor Index
474
475 ! Locals:
476 integer :: bodyIndex
477 integer :: bodyIndexbeg
478 integer :: currentChannelNumber
479
480 bodyIndexbeg = obs_headElem_i( obsSpaceData, OBS_RLN, headerIndex )
481
482 BODY: do bodyIndex = bodyIndexbeg, bodyIndexbeg + obs_headElem_i( obsSpaceData, OBS_NLV, headerIndex ) - 1
483 currentChannelNumber=nint(obs_bodyElem_r( obsSpaceData, OBS_PPP, bodyIndex ))-tvs_channelOffset(sensorIndex)
484 call obs_bodySet_i(obsSpaceData, OBS_FLG, bodyIndex, obsFlags(currentChannelNumber))
485 end do BODY
486
487 end subroutine csrbg_updateObsSpaceAfterQc
488
489end module bgckCSR_mod
