1MODULE biasCorrectionConv_mod
2 ! MODULE biasCorrectionConv_mod (prefix='bcc' category='1. High-level functionality')
3 !
4 !:Purpose: Performs bias correction for conventional observations.
5 !
6 use utilities_mod
7 use obsSpaceData_mod
8 use MathPhysConstants_mod
9 use midasMpi_mod
10 use bufr_mod
11 use codtyp_mod
12 use timeCoord_mod
13
14 implicit none
15 save
16 private
17 public :: bcc_applyAIBcor, bcc_applyGPBcor, bcc_applyUABcor
18 public :: bcc_biasActive
19
20 ! This variable is set to .true. when bcc_readConfig() [the routine that reads &NAMBIASCONV namelist]
21 ! is called for the first time to read/initialize the bias correction namelist variables.
22 logical :: initialized = .false.
23
24 integer, parameter :: nPhases=3, nLevels=5, nAircraftMax=100000
25 integer, parameter :: nStationMaxGP = 10000
26 integer, parameter :: phaseLevel = 3
27 integer, parameter :: phaseAscent = 5
28 integer, parameter :: phaseDescent = 6
29 integer, parameter :: phaseLevelIndex = 1
30 integer, parameter :: phaseAscentIndex = 2
31 integer, parameter :: phaseDescentIndex = 3
32
33 integer, parameter :: nSondesMax = 18, nMandLevs = 16
34 integer, parameter :: nStationMaxUA = 1000
35 integer, parameter :: Index500mb = 5
36
37 ! Missing values in the input bias correction files for AI, GP and UA families
38 real(8), parameter :: aiMissingValue = 99.d0
39 real(8), parameter :: gpMissingValue = -999.00d0
40 real(8), parameter :: uaMissingValue = -99.0d0
41
42 integer :: nbAircrafts, nbGpStations
43 real(8), allocatable :: AIttCorrections(:,:,:)
44 real(8), allocatable :: ztdCorrections(:)
45 real(8), allocatable :: ttCorrections(:,:,:,:), tdCorrections(:,:,:,:)
46 real(8), allocatable :: ttCorrectionsStn(:,:,:,:), tdCorrectionsStn(:,:,:,:)
47
48 character(len=9), allocatable :: aircraftIds(:), gpsStations(:), uaStations(:)
49 character(len=8) :: sondeTypes(nSondesMax)
50
51 real(8) :: mandLevs(nMandLevs), tolPress(nMandLevs)
52
53 logical :: bcc_aiBiasActive, bcc_gpBiasActive, bcc_uaBiasActive
54 logical, allocatable :: biasCorrPresentStype(:,:,:), biasCorrPresentStn(:,:,:)
55
56 ! Bias correction files (must be in program working directory)
57 character(len=8), parameter :: aiBcFile = "ai_bcors", gpBcFile = "gp_bcors"
58 character(len=14), parameter :: uaBcFileStype = "ua_bcors_stype", uaBcFileStn = "ua_bcors_stn"
59
60 integer, external :: fnom, fclos, newdate
61
62 ! Namelist variables
63 logical :: aiBiasActive ! Control if bias correction is applied to aircraft data
64 logical :: gpBiasActive ! Control if bias correction is applied to ground-based GPS data
65 logical :: uaBiasActive ! Control if bias correction is applied to radiosonde data
66 logical :: aiRevOnly ! Don't apply new correction but simply reverse any old corrections for AI
67 logical :: gpRevOnly ! Don't apply new correction but simply reverse any old corrections for GP
68 logical :: uaRevOnly ! Don't apply new correction but simply reverse any old corrections for UA
69 logical :: uaRejUnBcor ! Set DATA QC flag bit 11 on to exclude uncorrected UA observations from assimilation
70 integer :: uaNbiasCat ! Number of bias profile categories in UA bcor files, e.g. 1, or 2 for "asc" and "desc" phase categories
71 integer :: uaNlatBands ! Number of latitude bands in ua_bcors_stype bcor file (= 5 or 1). Set to 1 if there are no latitude bands in file
72 integer :: uaNprofsMin ! Min number of bias profiles required for a station/stype/time-of-day to use biases 'ua_bcors_stn' as corrections
73 character(len=8) :: nlSondeTypes(nSondesMax) ! List of radiosonde type names
74 integer :: nlSondeCodes(nSondesMax,20) ! List of radiosonde type codes
75 integer :: nlNbSondes ! Number of radiosonde types in lists
76
77 ! Structure to hold dictionary containing the BUFR sonde type codes associated with each sonde type
78 type sondeType
79 character(len=8) :: name
80 integer :: codes(20)
81 end type sondeType
82
83 type(sondeType), allocatable :: rsTypes(:)
84
85 namelist /nambiasconv/ aiBiasActive,gpBiasActive,aiRevOnly,gpRevOnly,uaBiasActive,uaRevOnly,uaNprofsMin,uaRejUnBcor,uaNbiasCat,uaNlatBands
86 namelist /namsondetypes/ nlNbSondes, nlSondeTypes, nlSondeCodes
87
88 ! 16 mandatory pressure levels (mb) on which radiosonde bias profiles are defined
89 data mandLevs /1000.d0,925.d0,850.d0,700.d0,500.d0,400.d0,300.d0,250.d0,200.d0,150.d0,100.d0,70.d0,50.d0,30.d0,20.d0,10.d0/
90 ! +/- tolerance (mb) for matching a radiosonde observation pressure to one of the mandatory levels
91 data tolPress / 10.d0, 10.d0, 10.d0, 10.d0, 10.d0, 10.d0, 10.d0, 5.d0, 5.d0, 5.d0, 5.d0, 2.d0, 2.d0, 2.d0, 1.d0, 1.d0/
92
93CONTAINS
94
95 !-----------------------------------------------------------------------
96 ! bcc_readConfig
97 !-----------------------------------------------------------------------
98 subroutine bcc_readConfig()
99 !
100 !:Purpose: Read NAMBIASCONV namelist section and NAMSONDETYPES section if uaBiasActive=true
101 !
102 implicit none
103
104 ! Locals:
105 integer :: ierr, nulnam, sondeIndex
106
107 if ( initialized ) then
108 write(*,*) "bcc_readConfig has already been called. Returning..."
109 return
110 end if
111
112 ! set default values for namelist variables
113 aiBiasActive = .false. ! bias correct AI data (TT)
114 gpBiasActive = .false. ! bias correct GP data (ZTD)
115 uaBiasActive = .false. ! bias correct UA data (TT,ES)
116 aiRevOnly = .false. ! AI: don't apply new correction but simply reverse any old corrections
117 gpRevOnly = .false. ! GP: don't apply new correction but simply reverse any old corrections
118 uaRevOnly = .false. ! UA: don't apply new correction but simply reverse any old corrections
119 uaRejUnBcor = .false. ! UA: Set DATA QC flag bit 11 on to exclude uncorrected UA observations from assimilation
120 uaNprofsMin = 100 ! UA: If Nprofs for a given stn/stype < uaNprofsMin, flag the bcors as unusable
121 uaNbiasCat = 1 ! UA: Number of bias profile categories in UA bcor files: 1, or 2 for "asc" and "desc" categories
122 uaNlatBands = 1 ! UA: Number of latitude bands in ua_bcors_stype file (1 or 5): 1 = no bands (global biases).
123 ! read in the namelist NAMBIASCONV
124 if ( utl_isNamelistPresent('nambiasconv','./flnml') ) then
125 nulnam = 0
126 ierr = fnom(nulnam,'./flnml','FTN+SEQ+R/O',0)
127 read(nulnam,nml=nambiasconv,iostat=ierr)
128 if ( ierr /= 0 ) call utl_abort('bcc_readConfig: Error reading namelist section NAMBIASCONV')
129 if ( mmpi_myid == 0 ) write(*,nml=nambiasconv)
130 ierr = fclos(nulnam)
131 else
132 write(*,*)
133 write(*,*) 'bcc_readconfig: NAMBIASCONV section is missing in the namelist. The default values will be used.'
134 end if
135
136 bcc_aiBiasActive = aiBiasActive
137 bcc_gpBiasActive = gpBiasActive
138 bcc_uaBiasActive = uaBiasActive
139
140 ! read in the namelist NAMSONDETYPES
141 if ( uaBiasActive .and. .not.uaRevOnly ) then
142 if ( utl_isNamelistPresent('namsondetypes','./flnml') ) then
143 nlNbSondes = MPC_missingValue_INT
144 nlSondeTypes(:) = 'empty'
145 nlSondeCodes(:,:) = MPC_missingValue_INT
146 nulnam = 0
147 ierr = fnom(nulnam,'./flnml','FTN+SEQ+R/O',0)
148 read(nulnam,nml=namsondetypes,iostat=ierr)
149 if ( ierr /= 0 ) call utl_abort('bcc_readConfig: Error reading namelist section NAMSONDETYPES')
150 if (nlNbSondes /= MPC_missingValue_INT) then
151 call utl_abort('bcc_readConfig: check namsondetypes namelist section, you should remove nlNbSondes')
152 end if
153 nlNbSondes = 0
154 do sondeIndex = 1, nSondesMax
155 if (trim(nlSondeTypes(sondeIndex)) == 'empty') exit
156 nlNbSondes = nlNbSondes + 1
157 end do
158 if ( mmpi_myid == 0 ) write(*,nml=namsondetypes)
159 ierr = fclos(nulnam)
160 else
161 write(*,*)
162 call utl_abort('bcc_readconfig: NAMSONDETYPES section is missing in the namelist!')
163 end if
164 if ( nlNbSondes > nSondesMax ) call utl_abort('bcc_readconfig: Number of sonde types in namelist exceeds nSondesMax!')
165 if ( nlNbSondes <= 0 ) call utl_abort('bcc_readconfig: Number of sonde types in namelist <= 0!')
166 allocate( rsTypes(nlNbSondes) )
167 do sondeIndex = 1, nlNbSondes
168 rsTypes(sondeIndex)%name = nlSondeTypes(sondeIndex)
169 rsTypes(sondeIndex)%codes(:) = nlSondeCodes(sondeIndex,:)
170 end do
171 end if
172
173 initialized = .true.
174
175 end subroutine bcc_readConfig
176
177 !-----------------------------------------------------------------------
178 ! bcc_UACorrection
179 !-----------------------------------------------------------------------
180 function bcc_UACorrection(timeOfDayX,corrNight,corrDay) result(uaCorrection)
181 !
182 !:Purpose: Returns a UA bias correction from day and night corrections
183 !
184 implicit none
185
186 ! Arguments:
187 real(8), intent(in) :: timeOfDayX ! 0.0 (night) <= timeOfDayX <= 1.0 (day), depends on solar_elev
188 real(8), intent(in) :: corrNight ! night bias correction
189 real(8), intent(in) :: corrDay ! day bias correction
190 ! Result:
191 real(8) :: uaCorrection
192
193 uaCorrection = MPC_missingValue_R8
194
195 if ( timeOfDayX == 0.0d0 ) then
196 uaCorrection = corrNight
197 else if ( timeOfDayX == 1.0d0 ) then
198 uaCorrection = corrDay
199 else
200 if ( corrNight /= MPC_missingValue_R8 .and. corrDay /= MPC_missingValue_R8 ) then
201 uaCorrection = timeOfDayX*corrDay + (1.0d0-timeOfDayX)*corrNight
202 end if
203 end if
204
205 end function bcc_UACorrection
206
207 !-----------------------------------------------------------------------
208 ! bcc_GetUACorrection
209 !-----------------------------------------------------------------------
210 subroutine bcc_GetUACorrection(varName,stnIndex,sondeTypeIndex,sondeType,biasProfileCategory,timeOfDayX,latband,obsPressure,corr,sourceCorr)
211 !
212 !:Purpose: Return a TT or TD bias correction (corr)
213
214 !
215 ! varName = (str) 'TT' or 'TD'
216 ! stnIndex = (int) station index (in uaStations) = -1 if station was not found in bcor file
217 ! sondeTypeIndex = (int) sonde type index (in rsTypes) = 0 if sonde-type code is not associated with any type
218 ! sondeType = (str) sonde-type from rsTypes
219 ! biasProfileCategory = (int) bias profile category = 1 (ascent/none) or 3 (descent)
220 ! timeOfDayX = (float) 0.0 (night) <= TimeOfDayX <= 1.0 (day), depends on solar_elev
221 ! latband = (int) 1-5
222 ! obsPressure = (float) pressure level (mb) of observation
223 !
224 ! Also returns info regarding the source of the correction returned:
225 ! sourceCorr = (str) 'none', 'stn', or 'stype'
226 !
227 ! Requires TT and TD correction profiles read from UA bcor file ua_bcors_stn and stored in
228 ! ttCorrectionsStn(stnIndex,sondeTypeIndex,j,MandLev)
229 ! tdCorrectionsStn(stnIndex,sondeTypeIndex,j,MandLev)
230 ! with backup corrections by sonde-type read from UA bcor file ua_bcors_stype and stored in
231 ! ttCorrections(sondeTypeIndex,latband,j,MandLev)
232 ! tdCorrections(sondeTypeIndex,latband,j,MandLev)
233 ! where MandLev = 1 (1000 mb) to 16 (10 mb)
234 ! If uaNbiasCat = 1 (biasProfileCategory = 1)
235 ! j = 1 (night)
236 ! 2 (day)
237 ! If uaNbiasCat = 2 (biasProfileCategory = 1 [ascent] or 3 [descent])
238 ! j = 1 (night-ascent)
239 ! 2 (day-ascent)
240 ! 3 (night-descent)
241 ! 4 (day-descent)
242 !
243 ! Interpolation of the correction profiles on mandatory levels is used to get the
244 ! correction at the observation level (obsPressure).
245 ! Persistence is applied for observations outside the range of the mandatory levels.
246 !
247 implicit none
248
249 ! Arguments:
250 integer, intent(in) :: stnIndex
251 integer, intent(in) :: sondeTypeIndex
252 integer, intent(in) :: biasProfileCategory
253 integer, intent(in) :: latband
254 real(8), intent(in) :: timeOfDayX
255 real(8), intent(in) :: obsPressure
256 character(len=*), intent(in) :: varName
257 character(len=*), intent(in) :: sondeType
258 real(8), intent(out) :: corr
259 character(len=*), intent(out) :: sourceCorr
260
261 ! Locals:
262 real(8) :: corrProfileStnDay(16), corrProfileStnNight(16), corrProfileStypeDay(16), corrProfileStypeNight(16)
263 real(8) :: corrDay, corrNight
264 real(8) :: weight1, weight2, deltaPressure, pressureAbove, pressureBelow, corrAbove, corrBelow
265 real(8) :: corrAboveNight, corrBelowNight, corrAboveDay, corrBelowDay
266 integer :: level, levelIndex
267 logical :: profileExistsStn, profileExistsStype, doInterp
268
269 sourceCorr = "none"
270
271 ! Bias correction by station and sonde-type
272 if ( stnIndex == -1 ) then
273 profileExistsStn = .false.
274 else
275 if ( sondeTypeIndex > 0 ) then
276 ! Check that correction profile for this station, sonde-type and TimeOfDay is available for use
277 profileExistsStn = biasCorrPresentStn(stnIndex,sondeTypeIndex,biasProfileCategory)
278 else ! There are no bias profiles for this station/stype combination
279 profileExistsStn = .false.
280 end if
281 end if
282
283 ! BACKUP: Bias correction by sonde-type
284 profileExistsStype = .false.
285 if ( trim(sondeType) /= 'unknown' .and. trim(sondeType) /= 'Others' .and. trim(sondeType) /= 'None' ) then
286 ! Check if correction profile for this sonde-type,latband,and TimeOfDay is available for use
287 profileExistsStype = biasCorrPresentStype(sondeTypeIndex,latband,biasProfileCategory)
288 end if
289
290 if ( .not.profileExistsStn .and. .not.profileExistsStype ) then
291 corr = MPC_missingValue_R8
292 return
293 end if
294
295 ! Fill the night and day bias correction profiles
296 if ( varName == 'TT' ) then
297 if ( profileExistsStn ) then
298 corrProfileStnNight(:) = ttCorrectionsStn(stnIndex,sondeTypeIndex,biasProfileCategory,:)
299 corrProfileStnDay(:) = ttCorrectionsStn(stnIndex,sondeTypeIndex,biasProfileCategory+1,:)
300 end if
301 if ( profileExistsStype ) then
302 corrProfileStypeNight(:) = ttCorrections(sondeTypeIndex,latband,biasProfileCategory,:)
303 corrProfileStypeDay(:) = ttCorrections(sondeTypeIndex,latband,biasProfileCategory+1,:)
304 end if
305 else if ( varName == 'TD' ) then
306 if ( profileExistsStn ) then
307 corrProfileStnNight(:) = tdCorrectionsStn(stnIndex,sondeTypeIndex,biasProfileCategory,:)
308 corrProfileStnDay(:) = tdCorrectionsStn(stnIndex,sondeTypeIndex,biasProfileCategory+1,:)
309 end if
310 if ( profileExistsStype ) then
311 corrProfileStypeNight(:) = tdCorrections(sondeTypeIndex,latband,biasProfileCategory,:)
312 corrProfileStypeDay(:) = tdCorrections(sondeTypeIndex,latband,biasProfileCategory+1,:)
313 end if
314 else
315 call utl_abort('bcc_GetUACorrection: Unsupported varName '//varName)
316 end if
317
318 corr = MPC_missingValue_R8
319 doInterp = .true.
320
321 !--------------------------------------------------------------------------------------
322 ! Get the correction at the observation level (obsPressure)
323 !--------------------------------------------------------------------------------------
324
325 ! Check if obsPressure is outside range of levels (no interpolation possible)
326 if ( obsPressure >= mandLevs(1) ) then
327 levelIndex = 1
328 else if ( obsPressure <= mandLevs(nMandLevs) ) then
329 levelIndex = nMandLevs
330 else
331 levelIndex = 0
332 end if
333
334 ! Check if obsPressure is close to one of the 16 mandatory levels (no interpolation needed)
335 if ( levelIndex == 0 ) then
336 do level = 1, nMandLevs
337 if ( abs(obsPressure-mandLevs(level)) <= tolPress(level) ) then
338 levelIndex = level
339 exit
340 end if
341 end do
342 end if
343
344 ! If obsPressure is close to one of the 16 mandatory levels get the correction for the mandatory level:
345 ! Use correction by station with correction by stype as backup
346 if ( levelIndex > 0 ) then
347 if ( profileExistsStn ) then
348 corrNight = corrProfileStnNight(levelIndex)
349 corrDay = corrProfileStnDay(levelIndex)
350 sourceCorr = "stn"
351 corr = bcc_UACorrection(timeOfDayX,corrNight,corrDay)
352 if ( corr == MPC_missingValue_R8 ) then
353 sourceCorr = "none"
354 if ( profileExistsStype ) then
355 corrNight = corrProfileStypeNight(levelIndex)
356 corrDay = corrProfileStypeDay(levelIndex)
357 sourceCorr = "stype"
358 corr = bcc_UACorrection(timeOfDayX,corrNight,corrDay)
359 if ( corr == MPC_missingValue_R8 ) sourceCorr = "none"
360 end if
361 end if
362 else
363 if ( profileExistsStype ) then
364 corrNight = corrProfileStypeNight(levelIndex)
365 corrDay = corrProfileStypeDay(levelIndex)
366 sourceCorr = "stype"
367 corr = bcc_UACorrection(timeOfDayX,corrNight,corrDay)
368 if ( corr == MPC_missingValue_R8 ) sourceCorr = "none"
369 end if
370 end if
371 doInterp = .false.
372 end if
373
374 ! or interpolate to get correction for observation level obsPressure:
375 ! Use corrections by station with corrections by stype as backup
376 if ( doInterp ) then
377 corrAbove = MPC_missingValue_R8
378 corrBelow = MPC_missingValue_R8
379
380 if ( profileExistsStn ) then
381 do level = 1, nMandLevs-1
382 if ( obsPressure <= mandLevs(level) .and. obsPressure > mandLevs(level+1) ) then
383 pressureBelow = mandLevs(level)
384 pressureAbove = mandLevs(level+1)
385 corrBelowNight = corrProfileStnNight(level)
386 corrAboveNight = corrProfileStnNight(level+1)
387 corrBelowDay = corrProfileStnDay(level)
388 corrAboveDay = corrProfileStnDay(level+1)
389 exit
390 end if
391 end do
392 sourceCorr = "stn"
393 if ( timeOfDayX == 0.0d0 ) then
394 corrBelow = corrBelowNight
395 corrAbove = corrAboveNight
396 else if ( timeOfDayX == 1.0d0 ) then
397 corrBelow = corrBelowDay
398 corrAbove = corrAboveDay
399 else
400 if ( corrBelowNight /= MPC_missingValue_R8 .and. corrBelowDay /= MPC_missingValue_R8 ) then
401 corrBelow = timeOfDayX*corrBelowDay + (1.0d0-timeOfDayX)*corrBelowNight
402 end if
403 if ( corrAboveNight /= MPC_missingValue_R8 .and. corrAboveDay /= MPC_missingValue_R8 ) then
404 corrAbove = timeOfDayX*corrAboveDay + (1.0d0-timeOfDayX)*corrAboveNight
405 end if
406 end if
407 end if
408
409 if ( corrAbove == MPC_missingValue_R8 .or. corrBelow == MPC_missingValue_R8 ) then
410 if ( profileExistsStype ) then
411 do level = 1, nMandLevs-1
412 if ( obsPressure <= mandLevs(level) .and. obsPressure > mandLevs(level+1) ) then
413 pressureBelow = mandLevs(level)
414 pressureAbove = mandLevs(level+1)
415 corrBelowNight = corrProfileStypeNight(level)
416 corrAboveNight = corrProfileStypeNight(level+1)
417 corrBelowDay = corrProfileStypeDay(level)
418 corrAboveDay = corrProfileStypeDay(level+1)
419 exit
420 end if
421 end do
422 sourceCorr = "stype"
423 if ( timeOfDayX == 0.0d0 ) then
424 corrBelow = corrBelowNight
425 corrAbove = corrAboveNight
426 else if ( timeOfDayX == 1.0d0 ) then
427 corrBelow = corrBelowDay
428 corrAbove = corrAboveDay
429 else
430 if ( corrBelowNight /= MPC_missingValue_R8 .and. corrBelowDay /= MPC_missingValue_R8 ) then
431 corrBelow = timeOfDayX*corrBelowDay + (1.0d0-timeOfDayX)*corrBelowNight
432 end if
433 if ( corrAboveNight /= MPC_missingValue_R8 .and. corrAboveDay /= MPC_missingValue_R8 ) then
434 corrAbove = timeOfDayX*corrAboveDay + (1.0d0-timeOfDayX)*corrAboveNight
435 end if
436 end if
437 end if
438 end if
439 deltaPressure = log10(pressureBelow)-log10(pressureAbove)
440 weight1 = (log10(pressureBelow) - log10(obsPressure)) / deltaPressure
441 weight2 = (log10(obsPressure) - log10(pressureAbove)) / deltaPressure
442 if ( corrAbove /= MPC_missingValue_R8 .and. corrBelow /= MPC_missingValue_R8 ) then
443 corr = weight1*corrAbove + weight2*corrBelow
444 else if ( corrAbove /= MPC_missingValue_R8 .and. max(weight1,weight2) == weight1 ) then
445 corr = corrAbove
446 else if ( corrBelow /= MPC_missingValue_R8 .and. max(weight1,weight2) == weight2 ) then
447 corr = corrBelow
448 else
449 corr = MPC_missingValue_R8
450 sourceCorr = "none"
451 end if
452 end if
453
454 end subroutine bcc_GetUACorrection
455
456 !-----------------------------------------------------------------------
457 ! bcc_StationIndex
458 !-----------------------------------------------------------------------
459 function bcc_StationIndex(station) result(stationIndexOut)
460 !
461 !:Purpose: Return the station index (order in array uaStations) corresponding to station
462 ! Returns -1 if station is not found in uaStations.
463 !
464 implicit none
465
466 ! Arguments:
467 character(len=*), intent(in) :: station
468 ! Result:
469 integer :: stationIndexOut
470
471 ! Locals:
472 integer :: stationIndex
473
474 if ( allocated(uaStations) ) then
475 stationIndexOut = -1
476 do stationIndex = 1, nStationMaxUA
477 if ( trim(uaStations(stationIndex)) == trim(station) ) then
478 stationIndexOut = stationIndex
479 exit
480 end if
481 end do
482 else
483 call utl_abort('bcc_StationIndex: array uaStations not allocated!')
484 end if
485
486 end function bcc_StationIndex
487
488 !-----------------------------------------------------------------------
489 ! bcc_SondeIndex
490 !-----------------------------------------------------------------------
491 function bcc_SondeIndex(sondeType) result(sondeIndex)
492 !
493 !:Purpose: Return the Sonde Type index (order in array rsTypes) corresponding to the SondeType
494 ! Requires array of sondeType structures (rsTypes) to be allocated and filled.
495 ! Returns -1 if SondeType is not found in rsTypes.
496 !
497 implicit none
498
499 ! Arguments:
500 character(len=*), intent(in) :: sondeType
501 ! Result:
502 integer :: sondeIndex
503
504 ! Locals:
505 integer :: typeIndex, ntypes
506
507 if ( allocated(rsTypes) ) then
508 ntypes = nlNbSondes
509 else
510 call utl_abort('bcc_SondeIndex: array rsTypes not allocated!')
511 end if
512
513 SondeIndex = -1
514 do typeIndex = 1, ntypes
515 if ( trim(rsTypes(typeIndex)%name) == trim(sondeType) ) then
516 sondeIndex = typeIndex
517 exit
518 end if
519 end do
520
521 end function bcc_SondeIndex
522
523 !-----------------------------------------------------------------------
524 ! bcc_GetSondeType
525 !-----------------------------------------------------------------------
526 subroutine bcc_GetSondeType(sondeTypeCode,sondeType,sondeTypeIndex)
527 !
528 !:Purpose: Returns the sonde type and index given a BUFR table sondeTypeCode (BUFR elem 002011)
529 ! Returns sondeType='unknown', sondeTypeIndex=0 if sondeTypeCode is not found.
530 !
531 !
532 ! Requires array of sonde_type structures (rsTypes) to be allocated and filled with the
533 ! sonde type codes associated with each sonde-type (read from namelist).
534 !
535 implicit none
536
537 ! Arguments:
538 integer, intent(in) :: sondeTypeCode
539 character(len=*), intent(out) :: sondeType
540 integer, intent(out) :: sondeTypeIndex
541
542 ! Locals:
543 integer :: typeIndex, ntypes, sondeCode
544
545 if ( allocated(rsTypes) ) then
546 ntypes = nlNbSondes
547 else
548 call utl_abort('bcc_GetSondeType: array rsTypes not allocated!')
549 end if
550
551 if ( sondeTypeCode == 190 .or. sondeTypeCode == 192 ) then ! NCAR dropsonde (BUFR only)
552 sondeCode = 13 ! RS92 code
553 else if ( sondeTypeCode == 191 .or. sondeTypeCode == 193 ) then ! NCAR dropsonde (BUFR only)
554 sondeCode = 41 ! RS41 code
555 else if ( sondeTypeCode >=100 .and. sondeTypeCode < 200 ) then
556 sondeCode = sondeTypeCode - 100
557 else
558 sondeCode = sondeTypeCode
559 end if
560
561 sondeType = 'unknown'
562 sondeTypeIndex = 0
563 do typeIndex = 1, ntypes
564 if ( ANY(rsTypes(typeIndex)%codes == sondeCode) ) then
565 sondeType = rsTypes(typeIndex)%name
566 sondeTypeIndex = typeIndex
567 exit
568 end if
569 end do
570
571 end subroutine bcc_GetSondeType
572
573 !-----------------------------------------------------------------------
574 ! bcc_GetSolarElevation
575 !-----------------------------------------------------------------------
576 subroutine bcc_GetSolarElevation(lat,lon,date,time,solarElev)
577 !
578 !:Purpose: Returns the solar elevation angle (degrees) given lat,lon,date(yyyymmdd),time(hhmm)
579 !
580
581 ! lat, lon = obsSpaceData header column OBS_LAT, OBS_LON (radians)
582 ! date, time = obsSpaceData header column OBS_DAT (yyyymmdd), OBS_ETM (hhmm) -or-
583 ! datestamp = tim_getDatestamp() -- date stamp for central (analysis) time
584 ! ier = newdate(datestamp,date,time,-3)
585 ! time = time/10000
586 !
587 implicit none
588
589 ! Arguments:
590 integer, intent(in) :: date ! yyyymmdd
591 integer, intent(in) :: time ! hhmm
592 real(8), intent(in) :: lat ! radians
593 real(8), intent(in) :: lon ! radians
594 real(8), intent(out) :: solarElev ! degrees
595
596 ! Locals:
597 integer :: days(13) = (/0,31,28,31,30,31,30,31,31,30,31,30,31/)
598 integer :: leap_years(7) = (/2016,2020,2024,2028,2032,2036,2040/)
599 integer :: yy, mmdd, mm, dd, hh, nn, doy
600 real(8) :: timeUTC, timeLCL, solarDec, hourAngle, csz, sza
601
602 yy = date/10000
603 mmdd = date - yy*10000
604 mm = mmdd/100
605 dd = mmdd - mm*100
606 hh = time/100
607 nn = time - hh*100
608
609 if ( ANY(leap_years == yy) ) days(3) = 29
610 doy = SUM(days(1:mm)) + dd
611 timeUTC = float(hh) + float(nn)/60.0d0
612 timeLCL = timeUTC + (lon*MPC_DEGREES_PER_RADIAN_R8)/15.0d0
613 if ( timeLCL < 0.0d0 ) timeLCL = 24.0d0 + timeLCL
614 if ( timeLCL > 24.0d0 ) timeLCL = timeLCL - 24.0d0
615 solarDec = 0.4093d0 * sin(((2.0d0*MPC_PI_R8)/365.0d0)*(284.0d0 + float(doy)))
616 hourAngle = 15.0d0*(timeLCL-12.0d0) * MPC_RADIANS_PER_DEGREE_R8
617 csz = sin(solarDec)*sin(lat) + cos(solarDec)*cos(lat)*cos(hourAngle)
618 sza = MPC_DEGREES_PER_RADIAN_R8 * acos(csz)
619 solarElev = 90.0d0 - sza
620
621 end subroutine bcc_GetSolarElevation
622
623 !-----------------------------------------------------------------------
624 ! bcc_GetTimeOfDay
625 !-----------------------------------------------------------------------
626 subroutine bcc_GetTimeOfDay(solarElev,timeOfDayX)
627 !
628 !:Purpose: Returns the time-of-day x value (0.0(night) <= x <= 1.0(day))
629 !
630 implicit none
631
632 ! Arguments:
633 real(8), intent(in) :: solarElev ! degrees
634 real(8), intent(out) :: timeOfDayX
635
636 if (solarElev < -7.5d0) then
637 timeOfDayX = 0.0d0
638 else if (solarElev < 22.5d0) then
639 timeOfDayX = (solarElev+7.5d0)/(22.5d0+7.5d0)
640 else
641 timeOfDayX = 1.0d0
642 end if
643
644 end subroutine bcc_GetTimeOfDay
645
646 !----------------------------------------------------------------------------------------
647 ! bcc_uaPhase
648 !----------------------------------------------------------------------------------------
649 function bcc_uaPhase(codeType) result(uaPhase)
650 !
651 !:Purpose: Returns the radiosonde phase (1=ascent, 2=descent) given a header code type
652 !
653 implicit none
654
655 ! Arguments:
656 integer, intent(in) :: codeType
657 ! Result:
658 integer :: uaPhase
659
660 if (codeType == codtyp_get_codtyp('tempdrop')) then
661 uaPhase = 2
662 else
663 uaPhase = 1
664 end if
665
666 end function bcc_uaPhase
667
668 !-----------------------------------------------------------------------
669 ! bcc_LatBand
670 !-----------------------------------------------------------------------
671 function bcc_LatBand(latInRadians) result(latBand)
672 !
673 !:Purpose: Returns latitude band number given latitude (radians)
674 !
675 implicit none
676
677 ! Arguments:
678 real(8), intent(in) :: latInRadians
679 ! Result:
680 integer :: latBand
681
682 ! Locals:
683 real(8) :: latInDegrees
684
685 if ( uaNlatBands /= 5 ) then
686 latBand = -1
687 else
688 latInDegrees = MPC_DEGREES_PER_RADIAN_R8 * latInRadians
689 if (latInDegrees < -90.0d0 ) latBand = 1 ! Should never be the case!
690 if (latInDegrees >= -90.0d0 .and. latInDegrees < -60.0d0) then
691 latBand = 1
692 else if (latInDegrees >= -60.0d0 .and. latInDegrees < -20.0d0) then
693 latBand = 2
694 else if (latInDegrees >= -20.0d0 .and. latInDegrees < 20.0d0) then
695 latBand = 3
696 else if (latInDegrees >= 20.0d0 .and. latInDegrees < 60.0d0) then
697 latBand = 4
698 else
699 latBand = 5
700 end if
701 end if
702
703 end function bcc_LatBand
704
705 !-----------------------------------------------------------------------
706 ! bcc_readAIBiases
707 !-----------------------------------------------------------------------
708 subroutine bcc_readAIBiases(biasEstimateFile)
709 !
710 !:Purpose: Read aircraft (AI) TT bias estimates from bias file and fill bias correction array ttCorrections.
711 ! The first line of the file is the number of aircraft plus one.
712 ! The rest of the file gives 15 values of Mean O-A for each aircraft, with each (AC,value) line written in format "a9,1x,f6.2".
713 ! The order is the same as what is written by genbiascorr script genbc.aircraft_bcor.py.
714 ! The first "aircraft" (AC name = BULKBCORS) values are the bulk biases by layer for All-AC (first 5 values), AIREP/ADS
715 ! (second 5 values) and AMDAR/BUFR (last 5 values).
716 ! Missing value = 99.0.
717 !
718 implicit none
719
720 ! Arguments:
721 character(len=*), intent(in) :: biasEstimateFile
722
723 ! Locals:
724 integer :: ierr, nulcoeff
725 integer :: stationIndex, phaseIndex, levelIndex
726 real(8) :: biasEstimate, correctionValue
727 character(len=9) :: stationId
728
729 if (.not.initialized) call bcc_readConfig()
730
731 if ( aiRevOnly ) return
732
733 nulcoeff = 0
734 ierr = fnom(nulcoeff, biasEstimateFile, 'FMT+R/O', 0)
735 if ( ierr /= 0 ) then
736 call utl_abort('bcc_readAIBiases: unable to open airplanes bias correction file ' // biasEstimateFile )
737 end if
738 read (nulcoeff, '(i5)', iostat=ierr ) nbAircrafts
739 if ( ierr /= 0 ) then
740 call utl_abort('bcc_readAIBiases: error 1 while reading airplanes bias correction file ' // biasEstimateFile )
741 end if
742
743 allocate( AIttCorrections(nAircraftMax,nPhases,nLevels) )
744 allocate( aircraftIds(nAircraftMax) )
745
746 AIttCorrections(:,:,:) = MPC_missingValue_R8
747
748 do stationIndex=1,nbAircrafts
749 do phaseIndex=1,3
750 do levelIndex=1,5
751 read (nulcoeff, *, iostat=ierr) stationId, biasEstimate
752 if ( ierr /= 0 ) then
753 call utl_abort('bcc_readAIBiases: error 2 while reading airplanes bias correction file ' // biasEstimateFile )
754 end if
755 if ( biasEstimate == aiMissingValue ) then
756 correctionValue = MPC_missingValue_R8
757 else
758 correctionValue = -1.0D0*biasEstimate
759 end if
760 AIttCorrections(stationIndex,phaseIndex,levelIndex) = correctionValue
761 aircraftIds(stationIndex) = stationId
762 !print*, stationIndex, phaseIndex, levelIndex, aircraftIds(stationIndex), AIttCorrections(stationIndex,phaseIndex,levelIndex)
763 end do
764 end do
765 end do
766 ierr = fclos(nulcoeff)
767
768 ! Check for bulk bias corrections at start of file
769 if ( aircraftIds(1) /= "BULKBCORS" ) then
770 call utl_abort('bcc_readAIBiases: Bulk bias corrections are missing in bias correction file!' )
771 end if
772
773 end subroutine bcc_readAIBiases
774
775 !-----------------------------------------------------------------------
776 ! bcc_applyAIBcor
777 !-----------------------------------------------------------------------
778 subroutine bcc_applyAIBcor(obsSpaceData)
779 !
780 !:Purpose: to apply aircraft (AI) bias corrections to observations in ObsSpaceData
781 !
782 implicit none
783
784 ! Arguments:
785 type(struct_obs), intent(inout) :: obsSpaceData
786
787 ! Locals:
788 integer :: headerIndex, bodyIndex, codtyp
789 integer :: flag, phase, bufrCode
790 integer :: phaseIndex, levelIndex, stationIndex, stationNumber
791 integer :: countTailCorrections, countBulkCorrections
792 integer :: headerFlag
793 real(8) :: corr, tt, oldCorr, pressure
794 character(len=9) :: stnid, stnId1, stnId2
795
796 if (.not.initialized) call bcc_readConfig()
797
798 if ( .not. aiBiasActive ) return
799
800 write(*,*) "bcc_applyAIBcor: start"
801
802 if ( .not.aiRevOnly ) call bcc_readAIBiases(aiBcFile)
803
804 countTailCorrections = 0
805 countBulkCorrections = 0
806
807 call obs_set_current_header_list(obsSpaceData,'AI')
808 HEADER: do
809 headerIndex = obs_getHeaderIndex(obsSpaceData)
810 if ( headerIndex < 0 ) exit HEADER
811
812 headerFlag = obs_headElem_i(obsSpaceData, OBS_ST1, headerIndex )
813
814 call obs_set_current_body_list(obsSpaceData, headerIndex)
815
816 BODY: do
817
818 bodyIndex = obs_getBodyIndex(obsSpaceData)
819 if ( bodyIndex < 0 ) exit BODY
820
821 if ( obs_bodyElem_i(obsSpaceData,OBS_ASS,bodyIndex) /= obs_assimilated ) cycle BODY
822
823 bufrCode = obs_bodyElem_i(obsSpaceData, OBS_VNM, bodyIndex )
824
825 if ( bufrCode == BUFR_NETT ) then
826 tt = obs_bodyElem_r(obsSpaceData, OBS_VAR, bodyIndex )
827 flag = obs_bodyElem_i(obsSpaceData, OBS_FLG, bodyIndex )
828 oldCorr = obs_bodyElem_r(obsSpaceData, OBS_BCOR, bodyIndex )
829 corr = MPC_missingValue_R8
830
831 if ( tt /= MPC_missingValue_R8 ) then
832
833 if ( btest(flag, 6) .and. oldCorr /= MPC_missingValue_R8 ) then
834 if ( btest(headerFlag, 15) ) then
835 tt = tt - oldCorr
836 else
837 tt = tt + oldCorr
838 end if
839 flag = ibclr(flag, 6)
840 end if
841 if ( aiRevOnly ) corr = 0.0D0
842
843 if ( .not.aiRevOnly ) then
844
845 pressure = obs_bodyElem_r(obsSpaceData, OBS_PPP, bodyIndex ) * MPC_MBAR_PER_PA_R8
846
847 ! Get level index and current (mar 2020) static bulk corrections applied to AI TT data at derivate stage
848 if ( (pressure <= 1100.d0) .and. (pressure > 700.d0) ) then
849 corr = 0.0D0
850 levelIndex = 5
851 else if ( (pressure <= 700.d0) .and. (pressure > 500.d0) ) then
852 corr = -0.1D0
853 levelIndex = 4
854 else if ( (pressure <= 500.d0) .and. (pressure > 400.d0) ) then
855 corr = -0.2D0
856 levelIndex = 3
857 else if ( (pressure <= 400.d0) .and. (pressure > 300.d0) ) then
858 corr = -0.3D0
859 levelIndex = 2
860 else if ( (pressure <= 300.d0) .and. (pressure > 100.d0) ) then
861 corr = -0.5D0
862 levelIndex = 1
863 else
864 levelIndex = 0
865 corr = 0.0D0
866 end if
867
868 codtyp = obs_headElem_i(obsSpaceData, OBS_ITY, headerIndex)
869
870 ! Default bulk corrections read from bcor file (applied if dynamic corrections are not availble for the aircraft)
871 select case( trim( codtyp_get_name(codtyp) ) )
872 case('airep','ads')
873 phaseIndex = phaseAscentIndex
874 case('amdar','acars')
875 phaseIndex = phaseDescentIndex
876 case default
877 write(*,*) 'bcc_applyAIBcor: codtyp=', codtyp
878 call utl_abort('bcc_applyAIBcor: unknown codtyp')
879 end select
880
881 if ( levelIndex /= 0 ) then
882 if ( AIttCorrections(1,phaseIndex,levelIndex) /= MPC_missingValue_R8 ) corr = AIttCorrections(1,phaseIndex,levelIndex)
883 countBulkCorrections = countBulkCorrections + 1
884 end if
885
886 headerIndex = obs_bodyElem_i(obsSpaceData, OBS_HIND, bodyIndex )
887 stnid = trim( obs_elem_c(obsSpaceData,'STID',headerIndex) )
888
889 ! on verifie si la station est dans le dictionnaire du fichier de correction de biais
890 !---------------------------------------------------------------------------------
891 stationNumber = 0
892 stnId2 = trim(stnid)
893 do stationIndex = 1, nbAircrafts
894 stnId1 = trim(aircraftIds(stationIndex))
895 if ( stnId2(2:9) == stnId1(1:8) ) stationNumber = stationIndex
896 end do
897
898 phase = obs_headElem_i(obsSpaceData, OBS_PHAS, headerIndex )
899
900 ! If the aircraft is in the bias correction file, get the new correction
901 ! AIttCorrections(stationNumber,phaseIndex,levelIndex) where
902 ! stationNumber = index for this AC in bias correction file (0 if not found)
903 ! phaseIndex = index for the 3 phases of flight (level, asc, desc)
904 ! levelIndex = index for the 5 layers (100-300, 300-400,400-500,500-700,700-1100)
905 ! and use it instead of bulk value (if it is not missing value).
906 if ( stationNumber /= 0 ) then
907 phaseIndex = 0
908 if ( phase == phaseLevel ) phaseIndex = phaseLevelIndex
909 if ( phase == phaseAscent ) phaseIndex = phaseAscentIndex
910 if ( phase == phaseDescent ) phaseIndex = phaseDescentIndex
911 if ( levelIndex /= 0 .and. phaseIndex /= 0 ) then
912 if ( AIttCorrections(stationNumber,phaseIndex,levelIndex) /= MPC_missingValue_R8 ) then
913 corr = AIttCorrections(stationNumber,phaseIndex,levelIndex)
914 countTailCorrections = countTailCorrections + 1
915 countBulkCorrections = countBulkCorrections - 1
916 end if
917 end if
918 end if
919
920 ! Apply the bias correction (bulk or new) and set the "bias corrected" bit in TT data flag ON
921 tt = tt + corr
922 flag = ibset(flag, 6)
923
924 end if
925
926 end if
927
928 call obs_bodySet_r( obsSpaceData, OBS_BCOR, bodyIndex, corr )
929 call obs_bodySet_r( obsSpaceData, OBS_VAR , bodyIndex, tt )
930 call obs_bodySet_i( obsSpaceData, OBS_FLG , bodyIndex, flag )
931
932 end if
933
934 end do BODY
935
936 ! Set header flag bit 15 on to indicate that bias correction has been ADDED to raw value.
937 ! In older versions of this routine, such as the version used in IC-3/GDPS v8.0, the correction was SUBTRACTED.
938 headerFlag = ibset(headerFlag, 15)
939 call obs_headSet_i( obsSpaceData, OBS_ST1, headerIndex, headerFlag )
940
941 end do HEADER
942
943 if ( countBulkCorrections + countTailCorrections /= 0 ) then
944 write (*, '(a58, i10)' ) "bcc_applyAIBcor: Number of obs with TT bulk correction = ", countBulkCorrections
945 write (*, '(a58, i10)' ) "bcc_applyAIBcor: Number of obs with TT tail correction = ", countTailCorrections
946 else
947 write(*,*) "bcc_applyAIBcor: No AI data found"
948 end if
949
950 if ( allocated(AIttCorrections) ) deallocate(AIttCorrections)
951 if ( allocated(aircraftIds) ) deallocate(aircraftIds)
952
953 write(*,*) "bcc_applyAIBcor: end"
954
955 end subroutine bcc_applyAIBcor
956
957 !-----------------------------------------------------------------------
958 ! bcc_readGPBiases
959 !-----------------------------------------------------------------------
960 subroutine bcc_readGPBiases(biasEstimateFile)
961 !
962 !:Purpose: Read GB-GPS bias estimates (mean ZTD O-A [mm] by station) and fill bias correction array ztdCorrections.
963 ! Missing value = -999.00
964 !
965 implicit none
966
967 ! Arguments:
968 character(len=*), intent(in) :: biasEstimateFile
969
970 ! Locals:
971 integer :: ierr, nulcoeff
972 integer :: stationIndex
973 real(8) :: biasEstimate
974 character(len=9) :: stationId
975
976 if (.not.initialized) call bcc_readConfig()
977
978 if ( gpRevOnly ) return
979
980 nulcoeff = 0
981 ierr = fnom(nulcoeff, biasEstimateFile, 'FMT+R/O', 0)
982 if ( ierr /= 0 ) then
983 call utl_abort('bcc_readGPBiases: unable to open GB-GPS bias correction file ' // biasEstimateFile )
984 end if
985 read (nulcoeff, '(i5)', iostat=ierr ) nbGpStations
986 if ( ierr /= 0 ) then
987 call utl_abort('bcc_readGPBiases: error 1 while reading GB-GPS bias correction file ' // biasEstimateFile )
988 end if
989
990 allocate( ztdCorrections(nStationMaxGP) )
991 allocate( gpsStations(nStationMaxGP) )
992
993 ztdCorrections(:) = MPC_missingValue_R8
994
995 do stationIndex=1,nbGpStations
996 read (nulcoeff, *, iostat=ierr) stationId, biasEstimate
997 if ( ierr /= 0 ) then
998 call utl_abort('bcc_readGPBiases: error 2 while reading GB-GPS bias correction file ' // biasEstimateFile )
999 end if
1000 if ( biasEstimate /= gpMissingValue ) ztdCorrections(stationIndex) = -1.0D0*(biasEstimate/1000.0D0) ! mm to m
1001 gpsStations(stationIndex) = stationId
1002 end do
1003 ierr = fclos(nulcoeff)
1004
1005 end subroutine bcc_readGPBiases
1006
1007 !-----------------------------------------------------------------------
1008 ! bcc_applyGPBcor
1009 !-----------------------------------------------------------------------
1010 subroutine bcc_applyGPBcor(obsSpaceData)
1011 !
1012 !:Purpose: to apply GB-GPS (GP) ZTD bias corrections to ZTD observations in ObsSpaceData
1013 !
1014 implicit none
1015
1016 ! Arguments:
1017 type(struct_obs), intent(inout) :: obsSpaceData
1018
1019 ! Locals:
1020 integer :: headerIndex, bodyIndex
1021 integer :: flag, bufrCode
1022 integer :: stationIndex, stationNumber
1023 integer :: nbCorrected
1024 real(8) :: corr, ztd, oldCorr
1025 character(len=9) :: stnid, stnId1, stnId2
1026
1027 if (.not.initialized) call bcc_readConfig()
1028
1029 if ( .not. gpBiasActive ) return
1030
1031 write(*,*) "bcc_applyGPBcor: start"
1032
1033 if ( .not.gpRevOnly ) call bcc_readGPBiases(gpBcFile)
1034
1035 nbCorrected = 0
1036
1037 call obs_set_current_header_list(obsSpaceData,'GP')
1038
1039 HEADER: do
1040 headerIndex = obs_getHeaderIndex(obsSpaceData)
1041 if ( headerIndex < 0 ) exit HEADER
1042
1043 call obs_set_current_body_list(obsSpaceData, headerIndex)
1044
1045 BODY: do
1046
1047 bodyIndex = obs_getBodyIndex(obsSpaceData)
1048 if ( bodyIndex < 0 ) exit BODY
1049
1050 if ( obs_bodyElem_i(obsSpaceData,OBS_ASS,bodyIndex) /= obs_assimilated ) cycle BODY
1051
1052 bufrCode = obs_bodyElem_i(obsSpaceData, OBS_VNM, bodyIndex )
1053
1054 if ( bufrCode == BUFR_NEZD ) then
1055
1056 ztd = obs_bodyElem_r(obsSpaceData, OBS_VAR, bodyIndex )
1057 oldCorr = obs_bodyElem_r(obsSpaceData, OBS_BCOR, bodyIndex )
1058 flag = obs_bodyElem_i(obsSpaceData, OBS_FLG, bodyIndex )
1059
1060 corr = MPC_missingValue_R8
1061
1062 if ( ztd /= MPC_missingValue_R8 ) then
1063
1064 ! Remove any previous bias correction
1065 if ( btest(flag, 6) .and. oldCorr /= MPC_missingValue_R8 ) then
1066 ztd = ztd - oldCorr
1067 flag = ibclr(flag, 6)
1068 end if
1069
1070 if ( .not. gpRevOnly ) then
1071 headerIndex = obs_bodyElem_i(obsSpaceData, OBS_HIND, bodyIndex )
1072 stnid = trim( obs_elem_c(obsSpaceData,'STID',headerIndex) )
1073
1074 ! on verifie si la station est dans le dictionnaire du fichier de correction de biais
1075 ! ---------------------------------------------------------------------------------
1076 stationNumber = 0
1077 stnId2 = trim(stnid)
1078 do stationIndex = 1, nbGpStations
1079 stnId1 = trim(gpsStations(stationIndex))
1080 if ( stnId2 == stnId1 ) stationNumber = stationIndex
1081 end do
1082
1083 if (stationNumber /= 0) then
1084 corr = ztdCorrections(stationNumber)
1085 end if
1086
1087 ! Apply the bias correction and set the "bias corrected" bit in ZTD data flag ON
1088 if ( corr /= MPC_missingValue_R8 ) then
1089 ztd = ztd + corr
1090 nbCorrected = nbCorrected + 1
1091 flag = ibset(flag, 6)
1092 else
1093 corr = 0.0D0
1094 end if
1095 else
1096 corr = 0.0D0
1097 end if
1098
1099 end if
1100
1101 call obs_bodySet_r( obsSpaceData, OBS_BCOR, bodyIndex, corr )
1102 call obs_bodySet_r( obsSpaceData, OBS_VAR , bodyIndex, ztd )
1103 call obs_bodySet_i( obsSpaceData, OBS_FLG , bodyIndex, flag )
1104
1105 end if
1106
1107 end do BODY
1108 end do HEADER
1109
1110 if ( nbCorrected /= 0 ) then
1111 write (*, '(a50, i10)' ) "bcc_applyGPBcor: Number of ZTD observations corrected = ", nbCorrected
1112 else
1113 write(*,*) "bcc_applyGPBcor: No GP data bias corrections made"
1114 end if
1115
1116 if ( allocated(ztdCorrections) ) deallocate( ztdCorrections )
1117 if ( allocated(gpsStations) ) deallocate( gpsStations )
1118
1119 write(*,*) "bcc_applyGPBcor: end"
1120
1121 end subroutine bcc_applyGPBcor
1122
1123
1124 !-----------------------------------------------------------------------
1125 ! bcc_readUABcorStype
1126 !-----------------------------------------------------------------------
1127 subroutine bcc_readUABcorStype(biasCorrectionFileName,nGroups)
1128 !
1129 !:Purpose: Read night and day TT, TD biases by SONDE TYPE and latitude band on 16 mandatory levels for UA family.
1130
1131 ! The first line is the sonde-type (and latitude band if uaNlatBands=5). Sonde type = 'END' for end of file.
1132 ! nGroups groups of of 16 lines follow, e.g., nGroups = 2: one group for ascending sonde observations and
1133 ! one for descending sonde observations
1134 ! -- 16 lines are 16 mandatory levels from 100000 Pa to 1000 Pa
1135 ! -- Lines contain four values: TT night-bias, TT day-bias, TD night-bias, TD day-bias
1136 ! Missing value = -99.0.
1137 !
1138 implicit none
1139
1140 ! Arguments:
1141 character(len=*), intent(in) :: biasCorrectionFileName
1142 integer, intent(in) :: nGroups
1143
1144 ! Locals:
1145 integer :: ierr, nulcoeff
1146 integer :: sondeTypeIndex, group, latBand, levelIndex, groupIndex, maxGroups
1147 real(8) :: ttBiasNight, ttBiasDay, tdBiasNight, tdBiasDay
1148 character(len=8) :: sondeType
1149
1150 if ( uaRevOnly ) return
1151
1152 nulcoeff = 0
1153 ierr = fnom(nulcoeff, biasCorrectionFileName, 'FMT+R/O', 0)
1154 if ( ierr /= 0 ) then
1155 call utl_abort('bcc_readUABcorStype: unable to open radiosonde bias correction file ' // biasCorrectionFileName )
1156 end if
1157
1158 maxGroups = nGroups*2
1159
1160 allocate( ttCorrections(nSondesMax,uaNlatBands,maxGroups,nMandLevs) )
1161 allocate( tdCorrections(nSondesMax,uaNlatBands,maxGroups,nMandLevs) )
1162 allocate( biasCorrPresentStype(nSondesMax,uaNlatBands,maxGroups) )
1163
1164 ttCorrections(:,:,:,:) = MPC_missingValue_R8
1165 tdCorrections(:,:,:,:) = MPC_missingValue_R8
1166 sondeTypes(:) = 'empty'
1167
1168 biasCorrPresentStype(:,:,:) = .false.
1169
1170 main_loop: do
1171
1172 if ( uaNlatBands == 1 ) then
1173 read (nulcoeff, *, iostat=ierr) sondeType
1174 latBand = 1
1175 else
1176 read (nulcoeff, *, iostat=ierr) sondeType, latBand
1177 end if
1178 if ( ierr /= 0 ) then
1179 call utl_abort('bcc_readUABcorStype: error reading sondeType, latBand in radiosonde bias correction file ' // biasCorrectionFileName )
1180 end if
1181 if ( trim(sondeType) == 'END' ) exit main_loop
1182 sondeTypeIndex = bcc_SondeIndex(sondeType)
1183 if ( sondeTypeIndex < 0 ) call utl_abort('bcc_readUABcorStype: Unknown sonde type '//sondeType)
1184 if ( latBand < 0 .or. latBand > uaNlatBands ) then
1185 write(*,*) 'sondeType, latband = ', sondeType, latBand
1186 call utl_abort('bcc_readUABcorStype: Latitude band out of range 0-5!')
1187 end if
1188 ! Skip sondeType/latband 0 (globe) if present
1189 if ( latBand == 0 ) then
1190 do group = 1, nGroups
1191 do levelIndex = 1, nMandLevs
1192 read (nulcoeff, *, iostat=ierr) ttBiasNight, ttBiasDay, tdBiasNight, tdBiasDay
1193 end do
1194 end do
1195 cycle main_loop
1196 end if
1197 sondeTypes(sondeTypeIndex) = sondeType
1198 do group = 1, nGroups
1199 groupIndex = (group-1)*2 + 1
1200 do levelIndex = 1, nMandLevs
1201 read (nulcoeff, *, iostat=ierr) ttBiasNight, ttBiasDay, tdBiasNight, tdBiasDay
1202 if ( ierr /= 0 ) then
1203 call utl_abort('bcc_readUABcorStype: error reading corrections in radiosonde bias correction file ' // biasCorrectionFileName )
1204 end if
1205 if ( ttBiasNight == uaMissingValue ) ttBiasNight = MPC_missingValue_R8
1206 if ( ttBiasDay == uaMissingValue ) ttBiasDay = MPC_missingValue_R8
1207 if ( tdBiasNight == uaMissingValue ) tdBiasNight = MPC_missingValue_R8
1208 if ( tdBiasDay == uaMissingValue ) tdBiasDay = MPC_missingValue_R8
1209 if ( ttBiasNight /= MPC_missingValue_R8 ) ttCorrections(sondeTypeIndex,latBand,groupIndex,levelIndex) = -1.0d0*ttBiasNight
1210 if ( ttBiasDay /= MPC_missingValue_R8 ) ttCorrections(sondeTypeIndex,latBand,groupIndex+1,levelIndex) = -1.0d0*ttBiasDay
1211 if ( tdBiasNight /= MPC_missingValue_R8 ) tdCorrections(sondeTypeIndex,latBand,groupIndex,levelIndex) = -1.0d0*tdBiasNight
1212 if ( tdBiasDay /= MPC_missingValue_R8 ) tdCorrections(sondeTypeIndex,latBand,groupIndex+1,levelIndex) = -1.0d0*tdBiasDay
1213 end do
1214 if ( ttCorrections(sondeTypeIndex,latBand,groupIndex,Index500mb) /= MPC_missingValue_R8 ) biasCorrPresentStype(sondeTypeIndex,latBand,groupIndex) = .true.
1215 if ( ttCorrections(sondeTypeIndex,latBand,groupIndex+1,Index500mb) /= MPC_missingValue_R8 ) biasCorrPresentStype(sondeTypeIndex,latBand,groupIndex+1) = .true.
1216 end do
1217
1218 end do main_loop
1219
1220 ierr = fclos(nulcoeff)
1221
1222 end subroutine bcc_readUABcorStype
1223
1224 !-----------------------------------------------------------------------
1225 ! bcc_readUABcorStn
1226 !-----------------------------------------------------------------------
1227 subroutine bcc_readUABcorStn(biasCorrectionFileName,nProfsMin,nGroups)
1228 !
1229 !:Purpose: Read TT, TD biases by STATION/sonde-type on 16 mandatory levels for UA family.
1230 !
1231 ! The first line is the station and sonde-type followed by number of profiles.
1232 !
1233 ! Station = 'END' for end of file.
1234 !
1235 ! nGroups groups of of 16 lines follow, e.g., nGroups = 2: one group for ascending sonde observations and
1236 ! one for descending sonde observations
1237 !
1238 ! -- 16 lines are 16 mandatory levels from 100000 Pa to 1000 Pa
1239 ! -- Lines contain 4 values: TT night-bias, TT day-bias, TD night-bias, TD day-bias
1240 !
1241 ! Missing value (file) = -99.0.
1242 !
1243 ! Sonde type in file is "None" if type was unknown/missing for the reports from a station.
1244 !
1245 ! biasCorrPresentStn(iStn,iType,TOD) = .true. for all TOD if total Nprofs >= nProfsMin (e.g. 100)
1246 !
1247 ! ttCorrectionsStn(iStn,iType,TOD,level) = MPC_missingValue_R8 if TTcorrectionValue == -99.0
1248 !
1249 implicit none
1250
1251 ! Arguments:
1252 character(len=*), intent(in) :: biasCorrectionFileName
1253 integer, intent(in) :: nProfsMin
1254 integer, intent(in) :: nGroups
1255
1256 ! Locals:
1257 integer :: ierr, nulcoeff, stationIndex, typeIndex
1258 integer :: group, levelIndex, nProfs, groupIndex, maxGroups
1259 real(8) :: ttBiasNight, ttBiasDay, tdBiasNight, tdBiasDay
1260 character(len=8) :: sondeType
1261 character(len=9) :: station, stnPrev
1262
1263 if ( uaRevOnly ) return
1264
1265 nulcoeff = 0
1266 ierr = fnom(nulcoeff, biasCorrectionFileName, 'FMT+R/O', 0)
1267 if ( ierr /= 0 ) then
1268 call utl_abort('bcc_readUABcorStn: unable to open radiosonde bias correction file ' // biasCorrectionFileName )
1269 end if
1270
1271 maxGroups = nGroups*2
1272
1273 allocate( ttCorrectionsStn(nStationMaxUA,nSondesMax,maxGroups,nMandLevs) )
1274 allocate( tdCorrectionsStn(nStationMaxUA,nSondesMax,maxGroups,nMandLevs) )
1275 allocate( uaStations(nStationMaxUA) )
1276 allocate( biasCorrPresentStn(nStationMaxUA,nSondesMax,maxGroups) )
1277
1278 ttCorrectionsStn(:,:,:,:) = MPC_missingValue_R8
1279 tdCorrectionsStn(:,:,:,:) = MPC_missingValue_R8
1280 uaStations(:) = 'empty'
1281 biasCorrPresentStn(:,:,:) = .false.
1282
1283 stationIndex = 0
1284 stnPrev = "toto"
1285 main_loop: do
1286
1287 read (nulcoeff, *, iostat=ierr) station, sondeType, nProfs
1288 if ( ierr /= 0 ) then
1289 call utl_abort('bcc_readUABcorStn: error reading station line in radiosonde bias correction file ' // biasCorrectionFileName )
1290 end if
1291 if ( trim(station) == 'END' ) exit main_loop
1292 typeIndex = bcc_SondeIndex(sondeType)
1293 if ( typeIndex < 0 ) call utl_abort('bcc_readUABcorStn: Unknown sonde type '//sondeType//' not in rsTypes defined in namelist')
1294
1295 ! If new station, add station to uaStations array.
1296 ! Assumes entries for the same station and different types are consecutive in bcor file
1297 if ( trim(station) /= trim(stnPrev) ) then
1298 stationIndex = stationIndex + 1
1299 uaStations(stationIndex) = station
1300 stnPrev = station
1301 end if
1302
1303 ! Determine which station/sondeType/time-of-day have enough profiles to use the bias corrections
1304 ! and set logical biasCorrPresentStn accordingly.
1305
1306 do group = 1, maxGroups
1307 if ( nProfs >= nProfsMin ) biasCorrPresentStn(stationIndex,typeIndex,group) = .true.
1308 end do
1309
1310 ! Read the bias correction profiles for this station/sondeType.
1311 do group = 1, nGroups
1312 groupIndex = (group-1)*2 + 1
1313 do levelIndex = 1, nMandLevs
1314 read (nulcoeff, *, iostat=ierr) ttBiasNight, ttBiasDay, tdBiasNight, tdBiasDay
1315 if ( ierr /= 0 ) then
1316 call utl_abort('bcc_readUABcorStn: error reading corrections in radiosonde bias correction file ' // biasCorrectionFileName )
1317 end if
1318 if ( ttBiasNight == uaMissingValue ) ttBiasNight = MPC_missingValue_R8
1319 if ( ttBiasDay == uaMissingValue ) ttBiasDay = MPC_missingValue_R8
1320 if ( tdBiasNight == uaMissingValue ) tdBiasNight = MPC_missingValue_R8
1321 if ( tdBiasDay == uaMissingValue ) tdBiasDay = MPC_missingValue_R8
1322 if ( ttBiasNight /= MPC_missingValue_R8 ) ttCorrectionsStn(stationIndex,typeIndex,groupIndex,levelIndex) = -1.0d0*ttBiasNight
1323 if ( ttBiasDay /= MPC_missingValue_R8 ) ttCorrectionsStn(stationIndex,typeIndex,groupIndex+1,levelIndex) = -1.0d0*ttBiasDay
1324 if ( tdBiasNight /= MPC_missingValue_R8 ) tdCorrectionsStn(stationIndex,typeIndex,groupIndex,levelIndex) = -1.0d0*tdBiasNight
1325 if ( tdBiasDay /= MPC_missingValue_R8 ) tdCorrectionsStn(stationIndex,typeIndex,groupIndex+1,levelIndex) = -1.0d0*tdBiasDay
1326 end do
1327 end do
1328
1329 end do main_loop
1330
1331 ierr = fclos(nulcoeff)
1332
1333 end subroutine bcc_readUABcorStn
1334
1335 !-----------------------------------------------------------------------
1336 ! bcc_applyUABcor
1337 !-----------------------------------------------------------------------
1338 subroutine bcc_applyUABcor(obsSpaceData)
1339 !
1340 !:Purpose: To apply bias corrections to radiosonde TT and ES observationa in obsSpaceData
1341 ! This public routine is called by external routines.
1342
1343 ! NOTE: We ADD the correction (negative of the bias) to the raw observation.
1344 ! We first SUBTRACT the old correction (if any) from the observation to remove (reverse) the old correction!
1345 ! Bias correction is put in obsSpaceData OBS_BCOR column.
1346 !
1347 ! NOTE: Bias correction is NOT applied if the sonde type is RS41. Observations from this sonde type are
1348 ! assumed to be unbiased.
1349 !
1350 ! Routine does nothing if uaBiasActive = .false. (just returns to calling routine)
1351 !
1352 implicit none
1353
1354 ! Arguments:
1355 type(struct_obs), intent(inout) :: obsSpaceData
1356
1357 ! Locals:
1358 integer :: headerIndex, bodyIndex, codtyp
1359 integer :: flag, bufrCode, sondeTypeCode, sondeTypeIndex, stnIndex
1360 integer :: date, time, latBand, groupIndex, phase
1361 integer :: countTTCorrections, countESCorrections, countUnknownStype, countMissingStype
1362 integer :: countUnknownStation, countTTCorrByStation, countTTCorrByStype, countTTObs, countRS41
1363 integer :: countCat1, countCat2, countNight, countDay, countDawnDusk
1364 integer :: ttBodyIndex, esBodyIndex
1365 real(8) :: tt, oldCorr, pressure, lat, lon, ttOriginal, es, esOriginal, td
1366 real(8) :: solarElev, corr, p1, p2, p3, p4
1367 real(8) :: timeOfDayX
1368 character(len=9) :: stnid, stnidPrev
1369 character(len=8) :: sondeType
1370 character(len=5) :: sourceCorr
1371 logical :: newStation, debug, stationFound, realRS41
1372
1373 if ( .not. uaBiasActive ) return
1374
1375 write(*,*) "bcc_applyUABcor: start"
1376
1377 debug = .false.
1378
1379 debug = debug .and. ( mmpi_myid == 0 )
1380
1381 ! Read the ascii files containing TT and TD bis profiles by sonde-type and by station
1382 if ( .not.uaRevOnly ) then
1383 call bcc_readUABcorStype(uaBcFileStype, uaNbiasCat)
1384 call bcc_readUABcorStn(uaBcFileStn, uaNprofsMin, uaNbiasCat)
1385 end if
1386
1387 countTTCorrections = 0
1388 countESCorrections = 0
1389 countUnknownStype = 0
1390 countMissingStype = 0
1391 countUnknownStation = 0
1392 countTTCorrByStation = 0
1393 countTTCorrByStype = 0
1394 countTTObs = 0
1395 countRS41 = 0
1396 countCat1 = 0
1397 countCat2 = 0
1398 countNight = 0
1399 countDay = 0
1400 countDawnDusk = 0
1401
1402 call obs_set_current_header_list(obsSpaceData,'UA')
1403
1404 stnidPrev = 'toto'
1405 HEADER: do ! Loop over UA headers (one header per pressure level)
1406 headerIndex = obs_getHeaderIndex(obsSpaceData)
1407 if ( headerIndex < 0 ) exit HEADER
1408
1409 call obs_set_current_body_list(obsSpaceData, headerIndex)
1410
1411 stnid = trim(obs_elem_c(obsSpaceData,'STID',headerIndex))
1412 if ( trim(stnid) /= trim(stnidPrev) ) then
1413 newStation = .true.
1414 stnidPrev = stnid
1415 else
1416 newStation = .false.
1417 end if
1418
1419 ! Get the information needed to apply the bias corrections from the first header for this station
1420 if ( newStation .and. .not.uaRevOnly ) then
1421
1422 ! Station index in list of stations (uaStations) from ua_bcors_stn file
1423 stationFound = .false.
1424 stnIndex = bcc_StationIndex(stnid)
1425 if (debug) write(*,*) 'stnid, index = ', stnid, stnIndex
1426 if ( stnIndex > 0 ) then
1427 stationFound = .true.
1428 else
1429 countUnknownStation = countUnknownStation + 1
1430 if (debug) write(*,*) 'Unknown station (not in ua_bcors_stn file) '//stnid
1431 end if
1432
1433 ! Date and lat,lon
1434 codtyp = obs_headElem_i(obsSpaceData, OBS_ITY, headerIndex) ! CODE TYPE
1435 date = obs_headElem_i(obsSpaceData, OBS_DAT, headerIndex) ! YYYYMMDD
1436 time = obs_headElem_i(obsSpaceData, OBS_ETM, headerIndex) ! HHMM
1437 lat = obs_headElem_r(obsSpaceData, OBS_LAT, headerIndex) ! radians!
1438 lon = obs_headElem_r(obsSpaceData, OBS_LON, headerIndex) ! radians!
1439
1440 ! Sonde phase: 1=asc, 2=desc
1441 phase = bcc_uaPhase(codtyp)
1442
1443 if (uaNbiasCat == 2) then
1444 groupIndex = (phase-1)*2 + 1
1445 else
1446 groupIndex = 1
1447 end if
1448
1449 if ( groupIndex == 1 ) then
1450 countCat1 = countCat1 + 1
1451 else
1452 countCat2 = countCat2 + 1
1453 end if
1454
1455 ! Get the sonde type index in rsTypes read from namelist
1456 sondeTypeCode = obs_headElem_i(obsSpaceData, OBS_RTP, headerIndex) ! sonde type BUFR code (WMO table)
1457 if (debug) then
1458 write(*,*) 'stnid, sondeTypeCode, date, time, lat'
1459 write(*,*) stnid, sondeTypeCode, date, time, lat*MPC_DEGREES_PER_RADIAN_R8
1460 end if
1461 if ( sondeTypeCode == -999 ) then
1462 sondeTypeCode = 0
1463 countMissingStype = countMissingStype + 1
1464 write (*,'(a60)') "bcc_applyUABcor: Missing sonde type at stn "//stnid
1465 end if
1466 call bcc_GetsondeType(sondeTypeCode,sondeType,sondeTypeIndex)
1467 if ( sondeTypeIndex == 0 ) then
1468 countUnknownStype = countUnknownStype + 1
1469 write (*,'(a70,i4)') "bcc_applyUABcor: Unknown sonde-type code at stn "//stnid//". Code = ", sondeTypeCode
1470 end if
1471 if (debug) write(*,*) 'sondeType, sondeTypeIndex = ', sondeType, sondeTypeIndex
1472
1473 ! We assume that a sonde type of "RS41" reported by Chinese UA stations is not correct
1474 realRS41 = .false.
1475 if ( trim(sondeType) == "RS41" ) then
1476 if ( stnid(1:1) == "5" ) then
1477 realRS41 = .false.
1478 else
1479 realRS41 = .true.
1480 end if
1481 end if
1482
1483 ! Time-of-day x-value
1484 call bcc_GetSolarElevation(lat,lon,date,time,solarElev)
1485 call bcc_GetTimeOfDay(solarElev,timeOfDayX)
1486 if ( timeOfDayX == 0.0 ) then
1487 countNight = countNight+1
1488 else if ( timeOfDayX == 1.0 ) then
1489 countDay = countDay+1
1490 else
1491 countDawnDusk = countDawnDusk+1
1492 end if
1493
1494 ! Latitude band
1495 if ( uaNlatBands == 1 ) then
1496 latBand = 1
1497 else
1498 latBand = bcc_LatBand(lat)
1499 if ( latBand == -1 ) then
1500 write(*,*) "uaNlatBands = ", uaNlatBands
1501 call utl_abort("bcc_applyUABcor: uaNlatBands must equal 1 or 5")
1502 end if
1503 end if
1504
1505 if (debug) write(*,*) 'solarElev, timeOfDayX, latBand = ',solarElev, timeOfDayX, latBand
1506
1507 end if
1508
1509 ttBodyIndex = -1
1510 esBodyIndex = -1
1511
1512 BODY: do ! Find the body indices for the TT and ES observations for this header (pressure level)
1513
1514 bodyIndex = obs_getBodyIndex(obsSpaceData)
1515 if ( bodyIndex < 0 ) exit BODY
1516
1517 if ( obs_bodyElem_i(obsSpaceData,OBS_ASS,bodyIndex) /= obs_assimilated ) cycle BODY
1518
1519 bufrCode = obs_bodyElem_i(obsSpaceData, OBS_VNM, bodyIndex )
1520
1521 ! Assumes that normal (non high-precision) codes are used for TT and ES observations in obsSpaceData
1522 if ( bufrCode == BUFR_NETT ) then
1523 ttBodyIndex = bodyIndex
1524 else if ( bufrCode == BUFR_NEES ) then
1525 esBodyIndex = bodyIndex
1526 end if
1527
1528 end do BODY
1529
1530 ! Reverse any old corrections and apply new corrections (if uaRevOnly=.false.)
1531
1532 ! TT bias correction
1533 if ( ttBodyIndex >= 0 ) then
1534 bodyIndex = ttBodyIndex
1535 ttOriginal = obs_bodyElem_r(obsSpaceData, OBS_VAR, bodyIndex )
1536 tt = ttOriginal
1537 pressure = obs_bodyElem_r(obsSpaceData, OBS_PPP, bodyIndex ) * MPC_MBAR_PER_PA_R8
1538 flag = obs_bodyElem_i(obsSpaceData, OBS_FLG, bodyIndex )
1539 oldCorr = obs_bodyElem_r(obsSpaceData, OBS_BCOR, bodyIndex )
1540 corr = MPC_missingValue_R8
1541 if ( debug .and. pressure == 500.0d0 ) write(*,*) 'TTin',ttBodyIndex,pressure,tt
1542 if ( tt /= MPC_missingValue_R8 ) then
1543 if ( btest(flag, 6) .and. oldCorr /= MPC_missingValue_R8 ) then
1544 tt = tt - oldCorr
1545 flag = ibclr(flag, 6)
1546 end if
1547 if ( uaRevOnly ) corr = 0.0d0
1548 if ( .not.uaRevOnly ) then
1549 ! Get the TT correction for this station, sonde type, time-of-day and level
1550 ! stnIndex = -1 if station not found in bcor file
1551 ! sondeTypeIndex = 0 if sonde-type code is not associated with any of the types in
1552 ! namelist (rsTypes) including stypes "Others" and "None" [code=0]
1553 countTTObs = countTTObs + 1
1554 if ( realRS41 ) then
1555 corr = MPC_missingValue_R8
1556 countRS41 = countRS41 + 1
1557 else
1558 call bcc_GetUACorrection('TT',stnIndex,sondeTypeIndex,sondeType,groupIndex,timeOfDayX,latBand,pressure,corr,sourceCorr)
1559 if ( sourceCorr == 'stn' ) then
1560 countTTCorrByStation = countTTCorrByStation + 1
1561 else if ( sourceCorr == 'stype' ) then
1562 countTTCorrByStype = countTTCorrByStype + 1
1563 end if
1564 if ( debug .and. pressure == 500.0d0 ) write(*,*) 'corrTT, source, groupIndex = ',corr, sourceCorr, groupIndex
1565 end if
1566 if ( corr /= MPC_missingValue_R8 ) then
1567 tt = tt + corr
1568 flag = ibset(flag, 6)
1569 countTTCorrections = countTTCorrections + 1
1570 else
1571 if ( .not.realRS41 .and. uaRejUnBcor ) flag = ibset(flag, 11)
1572 end if
1573 end if
1574 end if
1575 if ( debug .and. pressure == 500.0d0 ) write(*,*) 'TTout, corr, flag = ', tt, corr, flag
1576 call obs_bodySet_r( obsSpaceData, OBS_BCOR, bodyIndex, corr )
1577 call obs_bodySet_r( obsSpaceData, OBS_VAR , bodyIndex, tt )
1578 call obs_bodySet_i( obsSpaceData, OBS_FLG , bodyIndex, flag )
1579 end if
1580
1581 ! ES bias correction
1582 if ( esBodyIndex >= 0 .and. ttBodyIndex >= 0 ) then
1583 bodyIndex = esBodyIndex
1584 esOriginal = obs_bodyElem_r(obsSpaceData, OBS_VAR, bodyIndex )
1585 es = esOriginal
1586 pressure = obs_bodyElem_r(obsSpaceData, OBS_PPP, bodyIndex ) * MPC_MBAR_PER_PA_R8
1587 flag = obs_bodyElem_i(obsSpaceData, OBS_FLG, bodyIndex )
1588 oldCorr = obs_bodyElem_r(obsSpaceData, OBS_BCOR, bodyIndex )
1589 corr = MPC_missingValue_R8
1590 if ( debug .and. pressure == 500.0d0 ) write(*,*) 'ESin',esBodyIndex,pressure,es
1591 if ( es /= MPC_missingValue_R8 ) then
1592 if ( btest(flag, 6) .and. oldCorr /= MPC_missingValue_R8 ) then
1593 es = es - oldCorr
1594 flag = ibclr(flag, 6)
1595 end if
1596 if ( uaRevOnly ) corr = 0.0d0
1597 if ( .not.uaRevOnly ) then
1598 if ( realRS41 ) then
1599 corr = MPC_missingValue_R8
1600 else
1601 call bcc_GetUACorrection('TD',stnIndex,sondeTypeIndex,sondeType,groupIndex,timeOfDayX,latBand,pressure,corr,sourceCorr)
1602 if ( debug .and. pressure == 500.0d0 ) write(*,*) 'corrTD, source, groupIndex = ',corr, sourceCorr, groupIndex
1603 end if
1604 if ( corr /= MPC_missingValue_R8 ) then
1605 td = (ttOriginal - es) + corr
1606 es = tt - td
1607 if ( es < 0.0 ) es = 0.0d0
1608 if ( es > 30.0) es = 30.0d0
1609 corr = es - esOriginal
1610 if ( debug .and. pressure == 500.0d0 ) write(*,*) 'ES corr =',corr
1611 flag = ibset(flag, 6)
1612 countESCorrections = countESCorrections + 1
1613 else
1614 if ( .not.realRS41 .and. uaRejUnBcor ) flag = ibset(flag, 11)
1615 end if
1616 end if
1617 end if
1618 if ( debug .and. pressure == 500.0d0 ) write(*,*) 'ESout, corr, flag = ', es, corr, flag
1619 call obs_bodySet_r( obsSpaceData, OBS_BCOR, bodyIndex, corr )
1620 call obs_bodySet_r( obsSpaceData, OBS_VAR , bodyIndex, es )
1621 call obs_bodySet_i( obsSpaceData, OBS_FLG , bodyIndex, flag )
1622 end if
1623
1624 end do HEADER
1625
1626 if ( .not.uaRevOnly ) then
1627 write (*, '(a60, i10)' ) "bcc_applyUABcor: Number of TT obs corrected = ", countTTCorrections
1628 write (*, '(a60, i10)' ) "bcc_applyUABcor: Number of ES obs corrected = ", countESCorrections
1629 write (*, '(a60, i10)' ) "bcc_applyUABcor: Number obs with unknown station = ", countUnknownStation
1630 write (*, '(a60, i10)' ) "bcc_applyUABcor: Number obs with missing sonde type = ", countMissingStype
1631 write (*, '(a60, i10)' ) "bcc_applyUABcor: Number obs with unknown sonde code = ", countUnknownStype
1632 write (*, '(a40)' ) "----------------------------------------"
1633 write (*, '(a60, i10)' ) "bcc_applyUABcor: Total number of TT observations = ", countTTObs
1634 if ( countTTObs > 0 ) then
1635 p1 = (float(countTTCorrByStation)/float(countTTObs))*100.d0
1636 p2 = (float(countTTCorrByStype)/float(countTTObs))*100.d0
1637 p3 = 100.d0 - (p1+p2)
1638 p4 = (float(countRS41)/float(countTTObs))*100.
1639 write (*, '(a60, f7.2)' ) "bcc_applyUABcor: Percentage corrected by STATION = ", p1
1640 write (*, '(a60, f7.2)' ) "bcc_applyUABcor: Percentage corrected by SONDE-TYPE = ", p2
1641 write (*, '(a60, f7.2)' ) "bcc_applyUABcor: Percentage uncorrected = ", p3
1642 write (*, '(a60, f7.2)' ) "bcc_applyUABcor: Percentage RS41 = ", p4
1643 write (*, '(a60, i10)' ) "bcc_applyUABcor: Number of night profiles = ", countNight
1644 write (*, '(a60, i10)' ) "bcc_applyUABcor: Number of day profiles = ", countDay
1645 write (*, '(a60, i10)' ) "bcc_applyUABcor: Number of dawn-dusk profiles = ", countDawnDusk
1646 if (uaNbiasCat == 2) then
1647 write (*, '(a60, i10)' ) "bcc_applyUABcor: Number of ascent profiles = ", countCat1
1648 write (*, '(a60, i10)' ) "bcc_applyUABcor: Number of descent profiles = ", countCat2
1649 end if
1650 end if
1651 end if
1652
1653 if ( allocated(ttCorrections) ) deallocate(ttCorrections)
1654 if ( allocated(tdCorrections) ) deallocate(tdCorrections)
1655 if ( allocated(ttCorrectionsStn) ) deallocate(ttCorrectionsStn)
1656 if ( allocated(tdCorrectionsStn) ) deallocate(tdCorrectionsStn)
1657 if ( allocated(biasCorrPresentStype) ) deallocate(biasCorrPresentStype)
1658 if ( allocated(biasCorrPresentStn) ) deallocate(biasCorrPresentStn)
1659 if ( allocated(uaStations) ) deallocate(uaStations)
1660
1661 if ( allocated(rsTypes) ) deallocate(rsTypes)
1662
1663 write(*,*) "bcc_applyUABcor: end"
1664
1665 end subroutine bcc_applyUABcor
1666
1667
1668 !-----------------------------------------------------------------------
1669 ! bcc_biasActive
1670 !-----------------------------------------------------------------------
1671 function bcc_biasActive(obsFam) result(biasActive)
1672 !
1673 !:Purpose: returns True if bias correction is active for the given conventional observation family
1674 !
1675 implicit none
1676
1677 ! Arguments:
1678 character(len=*),intent(in) :: obsFam
1679 ! Result:
1680 logical :: biasActive
1681
1682 if (.not.initialized) call bcc_readConfig()
1683
1684 select case(trim(obsFam))
1685 case('GP')
1686 biasActive = bcc_gpBiasActive
1687 case('UA')
1688 biasActive = bcc_uaBiasActive
1689 case('AI')
1690 biasActive = bcc_aiBiasActive
1691 case default
1692 biasActive = .false.
1693 end select
1694
1695 end function bcc_biasActive
1696
1697end MODULE biasCorrectionConv_mod
