1module biasCorrectionSat_mod
2 ! MODULE biasCorrectionSat_mod (prefix='bcs' category='1. High-level functionality')
3 !
4 !:Purpose: Performs the bias correction for satellite radiance
5 ! data. This includes both the traditional approach based
6 ! on regression and the variational bias correction approach
7 ! for estimating the bias. Existing bias correction estimates
8 ! can also be applied to observations.
9 !
10 use utilities_mod
11 use ramDisk_mod
12 use MathPhysConstants_mod
13 use obsSpaceData_mod
14 use controlVector_mod
15 use midasMpi_mod
16 use rttov_const, only : ninst
17 use tovsNL_mod
18 use timeCoord_mod
19 use columnData_mod
20 use codePrecision_mod
21 use localizationFunction_mod
22 use HorizontalCoord_mod
23 use verticalCoord_mod
24 use gridStateVector_mod
25 use gridStateVectorFileIO_mod
26 use stateToColumn_mod
27 use codtyp_mod
28 use timeCoord_mod
29 use clibInterfaces_mod
30 use obserrors_mod
31 use fSQLite
32 use obsfiles_mod
33
34 implicit none
35 save
36 private
37
38 public :: bcs_setup, bcs_calcBias_tl, bcs_calcBias_ad, bcs_writeBias, bcs_finalize
39 public :: bcs_removeBiasCorrection, bcs_refreshBiasCorrection
40 public :: bcs_do_regression, bcs_filterObs, bcs_computeResidualsStatistics, bcs_calcBias
41 public :: bcs_removeOutliers, bcs_applyBiasCorrection
42 public :: bcs_mimicSatbcor
43 public :: bcs_readConfig
44 public :: bcs_dumpBiasToSqliteAfterThinning
45
46 type :: struct_chaninfo
47 integer :: numActivePredictors
48 logical :: isDynamic
49 integer :: channelNum
50 character(len=1) :: bcmode
51 character(len=1) :: bctype
52 integer, allocatable :: predictorIndex(:)
53 real(8), allocatable :: coeff(:)
54 real(8), allocatable :: coeffIncr(:)
55 real(8), allocatable :: coeff_fov(:)
56 real(8), allocatable :: coeff_offset(:)
57 integer :: coeff_nobs
58 real(8), allocatable :: coeffIncr_fov(:)
59 real(8), allocatable :: stddev(:)
60 real(8), allocatable :: coeffCov(:,:)
61 end type struct_chaninfo
62
63 type :: struct_bias
64 type(struct_chaninfo), allocatable :: chans(:)
65 integer :: numscan
66 integer :: numChannels
67 real(8), allocatable :: BHalfScanBias(:,:)
68 real(8), allocatable :: BMinusHalfScanBias(:,:)
69 end type struct_bias
70
71 type(struct_bias), allocatable :: bias(:)
72 type(struct_vco), pointer :: vco_mask => null()
73 type(struct_hco), pointer :: hco_mask => null()
74 type(struct_columnData) :: column_mask
75 logical :: initialized = .false.
76 logical :: bcs_mimicSatbcor
77 logical :: doRegression
78 integer, parameter :: NumPredictors = 7
79 integer, parameter :: NumPredictorsBcif = 6
80 integer, parameter :: maxfov = 120
81 integer, parameter :: maxNumInst = 25
82 integer, parameter :: maxPassiveChannels = 15
83
84 real(8), allocatable :: trialHeight300m1000(:)
85 real(8), allocatable :: trialHeight50m200(:)
86 real(8), allocatable :: trialHeight1m10(:)
87 real(8), allocatable :: trialHeight5m50(:)
88 real(8), allocatable :: RadiosondeWeight(:)
89 real(8), allocatable :: trialTG(:)
90 integer :: nobs
91 integer, external :: fnom, fclos
92 character(len=2), parameter :: predTab(0:7) = [ "SB", "KK","T1", "T2", "T3", "T4", "SV", "TG"]
93 integer :: passiveChannelNumber(maxNumInst)
94 ! Namelist variables
95 character(len=5) :: biasMode ! "varbc" for varbc, "reg" to compute bias correction coefficients by regression, "apply" to compute and apply bias correction
96 logical :: biasActive ! logical variable to activate the module
97 logical :: outstats ! flag to activate output of residual statistics in "reg" mode
98 logical :: mimicSatbcor ! in "reg" mode compute regression coefficients the same way as the original satbcor program
99 logical :: weightedestimate ! flag to activate radiosonde weighting for bias correction computation in "reg" mode
100 logical :: filterObs ! flag to activate additional observation filtering in "reg" mode. If it is .false. only observations selected for assimilation will be used in the linear regression
101 logical :: removeBiasCorrection ! flag to activate removal of an already present bias correction
102 logical :: refreshBiasCorrection !flag to replace an existing bias correction with a new one
103 logical :: centerPredictors ! flag to transparently remove predictor mean in "reg" mode (more stable problem; very little impact on the result)
104 logical :: outCoeffCov ! flag to activate output of coefficients error covariance (useful for EnKF system)
105 logical :: outOmFPredCov ! flag to activate output of O-F/predictors coefficients covariances and correlations
106 real(8) :: scanBiasCorLength ! if positive and .not. mimicSatBcor use error correlation between scan positions with the given correlation length
107 real(8) :: bg_stddev(NumPredictors) ! background error for predictors ("varbc" mode)
108 character(len=7) :: cinst(maxNumInst) ! to read the bcif file for each instrument in cinst
109 character(len=3) :: cglobal(maxNumInst) ! a "global" parameter and
110 integer :: nbscan(maxNumInst) ! the number of scan positions are necessary
111 integer :: passiveChannelList(maxNumInst, maxPassiveChannels)
112 ! To understand the meaning of the following parameters controling filtering,
113 ! please see https://wiki.cmc.ec.gc.ca/images/f/f6/Unified_SatRad_Dyn_bcor_v19.pdf pages 20-22
114 logical :: offlineMode ! flag to select offline mode for bias correction computation
115 logical :: allModeSsmis ! flag to select "ALL" mode for SSMIS
116 logical :: allModeTovs ! flag to select "ALL" mode for TOVS (AMSU-A, AMSU-B, MHS, ATMS, MWHS-2)
117 logical :: allModeCsr ! flag to select "ALL" mode for CSR (GOES, SEVIRI, MVIRI, ABI, etc..)
118 logical :: allModeHyperIr! flag to select "ALL" mode for hyperSpectral Infrared (AIRS, IASI, CrIS)
119 logical :: dumpToSqliteAfterThinning ! option to output all usefull parameters to sqlite files after thinning
120 namelist /nambiassat/ biasActive, biasMode, bg_stddev, removeBiasCorrection, refreshBiasCorrection
121 namelist /nambiassat/ centerPredictors, scanBiasCorLength, mimicSatbcor, weightedEstimate
122 namelist /nambiassat/ cglobal, cinst, nbscan, passiveChannelList, filterObs, outstats, outCoeffCov
123 namelist /nambiassat/ offlineMode, allModeSsmis, allModeTovs, allModeCsr, allModeHyperIr
124 namelist /nambiassat/ dumpToSqliteAfterThinning, outOmFPredCov
125contains
126
127 !-----------------------------------------------------------------------
128 ! bcs_readConfig
129 !-----------------------------------------------------------------------
130 subroutine bcs_readConfig()
131 !
132 !:Purpose: Read nambiassat namelist section
133 !
134 implicit none
135
136 ! Locals:
137 integer :: ierr, nulnam
138 logical, save :: firstCall = .true.
139 integer :: instrumentIndex, channelIndex
140
141 if (.not. firstCall) return
142 firstCall = .false.
143
144 ! set default values for namelist variables
145 biasActive = .false.
146 biasMode = "varbc"
147 bg_stddev(:) = 0.0d0
148 removeBiasCorrection = .false.
149 filterObs = .false.
150 refreshBiasCorrection = .false.
151 centerPredictors = .false.
152 mimicSatbcor = .true.
153 scanBiasCorLength = -1.d0
154 weightedEstimate = .false.
155 outCoeffCov = .false.
156 outOmFPredCov = .false.
157 nbscan(:) = -1
158 cinst(:) = "XXXXXXX"
159 cglobal(:) = "XXX"
160 outstats = .false.
161 offlineMode = .false.
162 allModeSsmis = .true.
163 allModeTovs = .true.
164 allModeCsr = .true.
165 allModeHyperIr = .false.
166 dumpToSqliteAfterThinning = .false.
167 passiveChannelNumber(:) = 0
168 passiveChannelList(:,:) = -1
169 !
170 ! read in the namelist NAMBIASSAT
171 if (utl_isNamelistPresent('nambiassat', './flnml')) then
172 nulnam = 0
173 ierr = fnom(nulnam, './flnml', 'FTN+SEQ+R/O', 0)
174 read(nulnam,nml=nambiassat,iostat=ierr)
175 if (ierr /= 0) call utl_abort('bcs_readConfig: Error reading namelist section nambiassat')
176 if (mmpi_myid == 0) write(*,nml=nambiassat)
177 ierr = fclos(nulnam)
178 else
179 write(*,*)
180 write(*,*) 'bcs_readconfig: nambiassat is missing in the namelist. The default value will be taken.'
181 end if
182
183 bcs_mimicSatbcor = mimicSatbcor
184 doRegression = (trim(biasMode) == "reg")
185
186 do instrumentIndex = 1, maxNumInst
187 do channelIndex = 1, maxPassiveChannels
188 if (passiveChannelList(instrumentIndex,channelIndex) > 0) then
189 passiveChannelNumber(instrumentIndex) = passiveChannelNumber(instrumentIndex) + 1
190 end if
191 end do
192 end do
193
194 end subroutine bcs_readConfig
195
196 !-----------------------------------------------------------------------
197 ! bcs_setup
198 !-----------------------------------------------------------------------
199 subroutine bcs_setup()
200 implicit none
201
202 ! Locals:
203 integer :: cvdim
204 integer :: iSensor, iPredictor, instIndex
205 integer :: iChan
206 integer :: iPred, jPred, kPred, iScan1, iScan2
207 character(len=85) :: bcifFile
208 character(len=10) :: instrName, instrNamecoeff, satNamecoeff
209 logical :: bcifExists
210 ! variables from background coeff file
211 integer :: nfov, exitCode
212 character(len=2) :: predBCIF(tvs_maxchannelnumber,numpredictorsBCIF)
213 integer :: canBCIF(tvs_maxchannelnumber), npredBCIF(tvs_maxchannelnumber), ncanBcif, npredictors
214 character(len=1) :: bcmodeBCIF(tvs_maxchannelnumber), bctypeBCIF(tvs_maxchannelnumber)
215 character(len=3) :: global
216 character(len=128) :: errorMessage
217 real(8), allocatable :: Bmatrix(:,:)
218
219 call bcs_readConfig()
220
221 cvdim = 0
222
223 if (biasActive) then
224
225 if (scanBiasCorLength > 0.d0) call lfn_Setup('FifthOrder')
226
227 allocate(bias(tvs_nSensors))
228
229 do iSensor = 1, tvs_nSensors
230
231 write(*,*) "bcs_setup: iSensor = ", iSensor
232
233 instrName = InstrNametoCoeffFileName(tvs_instrumentName(iSensor))
234 instrNamecoeff = InstrNameinCoeffFile(tvs_instrumentName(iSensor))
235 satNamecoeff = SatNameinCoeffFile(tvs_satelliteName(iSensor))
236
237 bcifFile = 'bcif_' // trim(instrName)
238
239 global = "XXX"
240 nfov = -1
241 do instIndex = 1, size(cinst)
242 if (trim(instrNamecoeff) == trim(cinst(instIndex))) then
243 global = cglobal(instIndex)
244 nfov = nbscan(instIndex)
245 end if
246 end do
247 if (nfov == -1) then
248 write(*,*) "bcs_setup: Problem with instrName ", instrNamecoeff
249 write(*,'(15(A10,1x))') cinst(:)
250 call utl_abort('bcs_setup: check nambiassat namelist')
251 end if
252
253 inquire(file=trim(bcifFile), exist=bcifExists)
254 if (bcifExists) then
255 !read bcif (Bias Correction Information File)
256 call read_bcif(bcifFile, tvs_isNameHyperSpectral(instrName), ncanBcif, &
257 canBCIF, bcmodeBCIF, bctypeBCIF, npredBCIF, predBCIF, global, exitcode)
258 if (exitcode /= 0) then
259 call utl_abort("bcs_setup: Problem in read_bcif while reading " // trim(bcifFile))
260 end if
261
262 bias(iSensor)%numChannels = ncanBcif
263
264 allocate(bias(iSensor)%chans(ncanBcif))
265
266 do ichan=1, ncanBcif
267 bias(iSensor) % chans(ichan) % channelNum = canBCIF(ichan + 1)
268 bias(iSensor) % chans(ichan) % coeff_nobs = 0
269 bias(iSensor) % chans(ichan) % bcmode = bcmodeBCIF(ichan + 1)
270 bias(iSensor) % chans(ichan) % bctype = bctypeBCIF(ichan + 1)
271 bias(iSensor) % chans(ichan) % isDynamic = ( (biasmode == "varbc" .and. bcmodeBCIF(ichan + 1) == "D") .or. biasmode /= "varbc")
272 npredictors = 1 + npredBCIF(ichan + 1)
273 bias(iSensor) % chans(ichan) % numActivePredictors = npredictors
274
275 allocate( bias(iSensor) % chans(ichan) % stddev( npredictors ) )
276 allocate( bias(iSensor) % chans(ichan) % coeffIncr( npredictors ) )
277 allocate( bias(iSensor) % chans(ichan) % coeff( npredictors ) )
278 allocate( bias(iSensor) % chans(ichan) % coeff_offset( npredictors ) )
279 allocate( bias(iSensor) % chans(ichan) % predictorIndex( npredictors ) )
280 bias(iSensor) % chans(ichan) % stddev(:) = 0.d0
281 bias(iSensor) % chans(ichan) % coeffIncr(:) = 0.d0
282 bias(iSensor) % chans(ichan) % coeff(:) = MPC_missingValue_R8
283 bias(iSensor) % chans(ichan) % coeff_offset(:) = 0.d0
284
285 bias(iSensor)%chans(ichan)% predictorIndex(1) = 1 !the constant term is always included
286 jPred = 1
287 do ipred = 1, npredBCIF(ichan + 1)
288 jPred = jPred + 1
289 select case(predBCIF(ichan + 1, ipred))
290 case('T1')
291 kpred = 2
292 case('T2')
293 kpred = 3
294 case('T3')
295 kpred = 4
296 case('T4')
297 kpred = 5
298 case('SV')
299 kpred = 6
300 case('TG')
301 kpred = 7
302 case default
303 write(errorMessage,*) "bcs_setup: Unknown predictor ", predBCIF(ichan+1, ipred), ichan, ipred
304 call utl_abort(errorMessage)
305 end select
306 bias(iSensor)%chans(ichan)%predictorIndex(jPred) = kpred
307 end do
308 end do
309 else
310 call utl_abort("bcs_setup: Error : " // trim(bcifFile) // " not present !")
311 end if
312
313 bias(iSensor)%numscan = nfov
314
315 allocate( bias(iSensor) % BHalfScanBias (nfov,nfov))
316 if (doRegression) allocate( bias(iSensor) % BMinusHalfScanBias (nfov,nfov))
317 allocate( Bmatrix(nfov,nfov))
318 do ichan=1, ncanBcif
319 allocate( bias(iSensor) % chans(ichan) % coeffIncr_fov(nfov))
320 allocate( bias(iSensor) % chans(ichan) % coeff_fov(nfov))
321 bias(iSensor) % chans(ichan) % coeffIncr_fov(:) = 0.d0
322 bias(iSensor) % chans(ichan) % coeff_fov(:) = MPC_missingValue_R8
323 end do
324
325 do ichan = 1, ncanBcif
326 if (bias(iSensor)%chans(ichan)%isDynamic) then
327 do iPredictor = 1, bias(iSensor)%chans(ichan)%numActivePredictors
328 bias(iSensor)%chans(ichan)%stddev(iPredictor) = bg_stddev(bias(iSensor)%chans(ichan)%PredictorIndex(iPredictor))
329 end do
330 end if
331 end do
332
333 if (trim(biasMode) == "varbc") then
334 !change dimension of control vector
335 do ichan = 1, ncanBcif
336 if (bias(iSensor)%chans(ichan)%isDynamic) &
337 cvdim = cvdim + bias(iSensor)%chans(ichan)%numActivePredictors - 1 + bias(iSensor)%numScan
338 end do
339 end if
340
341 if (allocated(Bmatrix)) then
342 if (scanBiasCorLength > 0.d0) then
343 do iScan2 = 1, nfov
344 do iScan1 = 1, nfov
345 Bmatrix(iScan1,iScan2) = bg_stddev(1) * bg_stddev(1) * lfn_Response(1.d0*abs(iScan1-iScan2), scanBiasCorLength)
346 end do
347 end do
348 else
349 Bmatrix(:,:) = 0.d0
350 do iScan1 = 1, nfov
351 Bmatrix(iScan1,iScan1) = bg_stddev(1) * bg_stddev(1)
352 end do
353 end if
354 bias(iSensor)%BHalfScanBias(:,:) = Bmatrix(:,:)
355 call utl_matsqrt(bias(iSensor)%BHalfScanBias, nfov, 1.d0, printInformation_opt=.true.)
356 if (doRegression) then
357 bias(iSensor)%BMinusHalfScanBias(:,:) = Bmatrix(:,:)
358 call utl_matsqrt(bias(iSensor)%BMinusHalfScanBias, nfov, -1.d0, printInformation_opt=.true.)
359 end if
360 deallocate(Bmatrix)
361 end if
362
363 end do
364
365 end if
366
367 if (trim(biasMode) == "varbc" .and. cvdim > 0) then
368 if (mmpi_myid > 0) cvdim = 0 ! for minimization, all coefficients only on task 0
369 call cvm_setupSubVector('BIAS', 'BIAS', cvdim)
370 end if
371
372 call bcs_readCoeffs() ! Read coefficient files in the case of bias correction application (biasMode=="apply")
373
374 end subroutine bcs_setup
375
376 !-----------------------------------------------------------------------
377 ! bcs_readCoeffs
378 !-----------------------------------------------------------------------
379 subroutine bcs_readCoeffs()
380 !
381 !:Purpose: Fill the bias structure with read static and dynamic bias correction coefficient files
382 !
383 implicit none
384
385 ! Locals:
386 integer :: iSensor, iSat, jchannel, jChan
387 integer :: satIndexDynamic, satIndexStatic
388 integer :: chanindexDynamic, chanindexStatic
389 character(len=10) :: instrName, instrNamecoeff, satNamecoeff
390 character(len=64) :: dynamicCoeffFile, staticCoeffFile
391 logical :: corrected
392 integer :: nfov, npredictors
393 character(len=10) :: satsDynamic(tvs_nsensors) ! satellite names
394 integer :: chansDynamic(tvs_nsensors,tvs_maxchannelnumber) ! channel numbers
395 real(8) :: fovbiasDynamic(tvs_nsensors,tvs_maxchannelnumber,maxfov)! bias as F(fov)
396 real(8) :: coeffDynamic(tvs_nsensors,tvs_maxchannelnumber,NumPredictors + 1)
397 integer :: nsatDynamic
398 integer :: nchanDynamic(tvs_nsensors) !number of channels
399 integer :: nfovDynamic
400 integer :: npredDynamic(tvs_nsensors,tvs_maxchannelnumber) !number of predictors
401 character(len=7) :: cinstrumDynamic ! string: instrument (e.g. AMSUB)
402 character(len=2) :: ptypesDynamic(tvs_nsensors,tvs_maxchannelnumber,NumPredictors)
403 integer :: ndataDynamic(tvs_nsensors,tvs_maxchannelnumber) !number of channels
404 character(len=10) :: satsStatic(tvs_nsensors) ! satellite names
405 integer :: chansStatic(tvs_nsensors,tvs_maxchannelnumber) !channel numbers
406 real(8) :: fovbiasStatic(tvs_nsensors,tvs_maxchannelnumber,maxfov)!bias as F(fov)
407 real(8) :: coeffStatic(tvs_nsensors, tvs_maxchannelnumber,NumPredictors + 1)
408 integer :: nsatStatic
409 integer :: nchanStatic(tvs_nsensors) !number of channels
410 integer :: nfovStatic
411 integer :: npredStatic(tvs_nsensors,tvs_maxchannelnumber) ! number of predictors
412 character(len=7) :: cinstrumStatic ! string: instrument (e.g. AMSUB) 9
413 character(len=2) :: ptypesStatic(tvs_nsensors,tvs_maxchannelnumber,NumPredictors)
414 integer :: ndataStatic(tvs_nsensors,tvs_maxchannelnumber) ! number of channels
415
416 if (biasActive .and. biasMode == "apply") then
417
418 ! 1 fichier de coefficient par intrument avec les differentes plateformes
419 ! Cas particulier GEORAD (CSR)
420
421 do iSensor = 1, tvs_nSensors
422
423 write(*,*) "bcs_readCoeffs: iSensor = ", iSensor
424
425 instrName = InstrNametoCoeffFileName(tvs_instrumentName(iSensor))
426 instrNamecoeff = InstrNameinCoeffFile(tvs_instrumentName(iSensor))
427 satNamecoeff = SatNameinCoeffFile(tvs_satelliteName(iSensor))
428
429 dynamicCoeffFile = "coeff_file_" // trim(instrName)
430 staticCoeffFile = "static_coeff_file_" // trim(instrName)
431
432 if ( tvs_isNameGeostationary(instrName)) then
433 dynamicCoeffFile = trim(dynamicCoeffFile) // "." // trim(satNamecoeff)
434 staticCoeffFile = trim(staticCoeffFile) // "." // trim(satNamecoeff)
435 end if
436
437 call read_coeff(satsDynamic, chansDynamic, fovbiasDynamic, coeffDynamic, nsatDynamic, nchanDynamic, nfovDynamic, &
438 npredDynamic, cinstrumDynamic, dynamicCoeffFile, ptypesDynamic, ndataDynamic)
439
440 call read_coeff(satsStatic, chansStatic, fovbiasStatic, coeffStatic, nsatStatic, nchanStatic, nfovStatic, &
441 npredStatic, cinstrumStatic, staticCoeffFile, ptypesStatic, ndataStatic)
442 write(*,*) "bcs_readCoeffs: cinstrumDynamic = ", cinstrumDynamic
443 write(*,*) "bcs_readCoeffs: cinstrumStatic = ", cinstrumStatic
444
445 satIndexDynamic = -1
446 do iSat = 1, nsatDynamic
447 if (trim(satNameCoeff) /= trim(satsDynamic(iSat)) .or. trim(instrNamecoeff) /= trim(cinstrumDynamic)) cycle
448 satIndexDynamic = iSat
449 end do
450
451 satIndexStatic = -1
452 do iSat = 1, nsatStatic
453 if (trim(satNameCoeff) /= trim(satsStatic(iSat)) .or. trim(instrNamecoeff) /= trim(cinstrumStatic)) cycle
454 satIndexStatic = iSat
455 end do
456
457 nfov = bias(iSensor)%numscan
458
459 do jChannel = 1, bias(iSensor)%numChannels
460
461 chanindexDynamic = -1
462 if (satIndexDynamic > 0) then
463 do jChan = 1, nchanDynamic(satIndexDynamic)
464 if (chansDynamic(satIndexDynamic,jChan) == bias(iSensor)%chans(jChannel)%channelNum) then
465 chanindexDynamic = jChan
466 exit
467 end if
468 end do
469 end if
470
471 chanindexStatic = -1
472 if (satIndexStatic > 0) then
473 do jChan = 1, nchanStatic(satIndexStatic)
474 if (chansStatic(satIndexStatic,jChan) == bias(iSensor)%chans(jChannel)%channelNum) then
475 chanindexStatic = jChan
476 exit
477 end if
478 end do
479 end if
480 npredictors = bias(iSensor)%chans(jChannel)%numActivePredictors
481
482 corrected = .true.
483 select case(bias(iSensor)%chans(jChannel)%bcmode)
484 case("D")
485 if (chanindexDynamic > 0) then
486 if (bias(iSensor)%chans(jChannel)%bctype=="C") &
487 bias(iSensor)%chans(jChannel)%coeff(1:npredictors) = coeffDynamic(satIndexDynamic,chanindexDynamic,1:npredictors)
488 if (bias(iSensor)%chans(jChannel)%bctype=="C" .or. bias(iSensor)%chans(jChannel)%bctype=="F") &
489 bias(iSensor)%chans(jChannel)%coeff_fov(1:nfov) = fovbiasDynamic(satIndexDynamic,chanindexDynamic,1:nfov)
490 else if (chanindexStatic > 0) then
491 if (bias(iSensor)%chans(jChannel)%bctype=="C") &
492 bias(iSensor)%chans(jChannel)%coeff(1:npredictors) = coeffStatic(satIndexStatic,chanindexStatic,1:npredictors)
493 if (bias(iSensor)%chans(jChannel)%bctype=="C" .or. bias(iSensor)%chans(jChannel)%bctype=="F") &
494 bias(iSensor)%chans(jChannel)%coeff_fov(1:nfov) = fovbiasStatic(satIndexStatic,chanindexStatic,1:nfov)
495 else
496 corrected = .false.
497 end if
498 case("S")
499 if (chanindexStatic > 0) then
500 if (bias(iSensor)%chans(jChannel)%bctype=="C") &
501 bias(iSensor)%chans(jChannel)%coeff(1:npredictors) = coeffStatic(satIndexStatic,chanindexStatic,1:npredictors)
502 if (bias(iSensor)%chans(jChannel)%bctype=="C" .or. bias(iSensor)%chans(jChannel)%bctype=="F") &
503 bias(iSensor)%chans(jChannel)%coeff_fov(1:nfov) = fovbiasStatic(satIndexStatic,chanindexStatic,1:nfov)
504 else
505 corrected = .false.
506 end if
507 end select
508
509 if (.not. corrected) then
510 write(*,*) "bcs_readCoeffs: Warning: channel ", bias(iSensor)%chans(jChannel)%channelNum, " of ", &
511 trim(instrName), " ", trim(satNamecoeff), " not corrected!"
512 end if
513
514 end do
515
516 end do
517
518 end if
519
520 end subroutine bcs_readCoeffs
521
522 !-----------------------------------------------------------------------
523 ! bcs_computePredictorBiases
524 !-----------------------------------------------------------------------
525 subroutine bcs_computePredictorBiases(obsSpaceData)
526 !
527 !:Purpose: to compute predictor average
528 !
529 implicit none
530
531 ! Arguments:
532 type(struct_obs), intent(inout) :: obsSpaceData
533
534 ! Locals:
535 real(8) :: predictor(NumPredictors)
536 integer :: iobs, nsize, i, j, npred
537 integer :: headerIndex, idatyp, indxtovs
538 integer :: iSensor, iFov, iPredictor, ierr
539 integer :: bodyIndex, jpred, chanIndx
540 real(8), allocatable :: temp_offset(:,:)
541 integer, allocatable :: temp_nobs(:)
542 real(8), allocatable :: temp_offset2(:,:,:)
543 integer, allocatable :: temp_nobs2(:,:)
544
545 if (centerPredictors) then
546 write(*,*) "bcs_computePredictorBiases: start"
547 npred = 0
548 do iSensor = 1, tvs_nSensors
549 npred = max(npred, maxval(bias(iSensor)%chans(:)%numActivePredictors))
550 end do
551 allocate(temp_offset2(tvs_nsensors,maxval(bias(:)%numChannels),2:npred))
552 allocate(temp_nobs2(tvs_nsensors,maxval(bias(:)%numChannels)))
553 temp_offset2(:,:,:) = 0.d0
554 temp_nobs2(:,:) = 0
555
556 call obs_set_current_header_list(obsSpaceData, 'TO')
557 iobs = 0
558 HEADER: do
559 headerIndex = obs_getHeaderIndex(obsSpaceData)
560 if (headerIndex < 0) exit HEADER
561
562 ! process only radiance data to be assimilated?
563 idatyp = obs_headElem_i(obsSpaceData, OBS_ITY, headerIndex)
564 if (.not. tvs_isIdBurpTovs(idatyp)) then
565 write(*,*) 'bcs_computePredictorBiases: warning unknown radiance codtyp present check NAMTOVSINST', idatyp
566 cycle HEADER
567 end if
568 iobs = iobs + 1
569
570 indxTovs = tvs_tovsIndex(headerIndex)
571 if (indxTovs < 0) cycle HEADER
572
573 iSensor = tvs_lsensor(indxTovs)
574
575 call obs_set_current_body_list(obsSpaceData, headerIndex)
576 iFov = obs_headElem_i(obsSpaceData, OBS_FOV, headerIndex)
577
578 BODY: do
579 bodyIndex = obs_getBodyIndex(obsSpaceData)
580 if (bodyIndex < 0) exit BODY
581
582 if (obs_bodyElem_i(obsSpaceData, OBS_ASS, bodyIndex) /= obs_assimilated) cycle BODY
583
584 call bcs_getChannelIndex(obsSpaceData, iSensor, chanIndx, bodyIndex)
585 if (chanindx > 0) then
586 call bcs_getPredictors(predictor, headerIndex, iobs, chanIndx, obsSpaceData)
587 do iPredictor = 2, bias(iSensor)%chans(chanIndx)%NumActivePredictors
588 jPred = bias(iSensor)%chans(chanIndx)%PredictorIndex(iPredictor)
589 temp_offset2(iSensor,chanIndx,iPredictor) = temp_offset2(iSensor,chanIndx,iPredictor) + predictor(jPred)
590 end do
591 temp_nobs2(iSensor,chanIndx) = temp_nobs2(iSensor,chanIndx) + 1
592 end if
593 end do BODY
594 end do HEADER
595
596 allocate(temp_offset(maxval(bias(:)%numChannels),2:npred))
597 allocate(temp_nobs(maxval(bias(:)%numChannels)))
598
599 do iSensor = 1, tvs_nSensors
600 temp_offset(:,:) = 0.0d0
601 temp_offset(:,:) = temp_offset2(iSensor,:,:)
602 call mmpi_allreduce_sumR8_2d( temp_offset, "GRID" )
603
604 do i = 1, bias(iSensor)%numChannels
605 do j = 2, bias(iSensor)%chans(i)%numActivePredictors
606 bias(iSensor)%chans(i)%coeff_offset(j) = temp_offset(i,j)
607 end do
608 end do
609
610 temp_nobs(:) = 0
611 nsize = size(temp_nobs)
612 call rpn_comm_allreduce(temp_nobs2(iSensor,:), temp_nobs, nsize, "mpi_integer", &
613 "mpi_sum", "GRID", ierr)
614 if (ierr /= 0) then
615 call utl_abort("bcs_computePredictorBiases: Erreur de communication MPI 2")
616 end if
617
618 do i = 1, bias(iSensor)%numChannels
619 bias(iSensor)%chans(i)%coeff_nobs = temp_nobs(i)
620 end do
621 do i = 1, bias(iSensor)%numChannels
622 if (bias(iSensor)%chans(i)%coeff_nobs > 0) then
623 bias(iSensor)%chans(i)%coeff_offset(:) = bias(iSensor)%chans(i)%coeff_offset / bias(iSensor)%chans(i)%coeff_nobs
624 end if
625 end do
626 end do
627
628 deallocate(temp_offset)
629 deallocate(temp_nobs)
630 deallocate(temp_offset2)
631 deallocate(temp_nobs2)
632
633 write(*,*) "bcs_computePredictorBiases: end"
634 end if
635
636 end subroutine bcs_computePredictorBiases
637
638 !-----------------------------------------------------------------------
639 ! bcs_calcBias
640 !-----------------------------------------------------------------------
641 subroutine bcs_calcBias(obsSpaceData, columnTrlOnTrlLev)
642 !
643 !:Purpose: to fill OBS_BCOR column of ObsSpaceData body with bias correction computed from read coefficient file
644 !
645 implicit none
646
647 ! Arguments:
648 type(struct_obs), intent(inout) :: obsSpaceData
649 type(struct_columnData), intent(inout) :: columnTrlOnTrlLev
650
651 ! Locals:
652 integer :: headerIndex, bodyIndex, iobs, indxtovs, idatyp
653 integer :: iSensor, iPredictor, chanIndx
654 integer :: iScan, iFov, jPred
655 real(8) :: predictor(NumPredictors)
656 real(8) :: biasCor
657
658 if (.not. biasActive) return
659
660 write(*,*) "bcs_calcBias: start"
661
662 if (.not. allocated(trialHeight300m1000)) then
663 call bcs_getTrialPredictors(obsSpaceData, columnTrlOnTrlLev)
664 end if
665
666 iobs = 0
667 call obs_set_current_header_list(obsSpaceData, 'TO')
668 HEADER: do
669 headerIndex = obs_getHeaderIndex(obsSpaceData)
670 if (headerIndex < 0) exit HEADER
671
672 ! process only radiance data to be assimilated?
673 idatyp = obs_headElem_i(obsSpaceData, OBS_ITY, headerIndex)
674 if (.not. tvs_isIdBurpTovs(idatyp)) then
675 write(*,*) 'bcs_calBias: warning unknown radiance codtyp present check NAMTOVSINST', idatyp
676 cycle HEADER
677 end if
678
679 iobs = iobs + 1
680
681 indxtovs = tvs_tovsIndex(headerIndex)
682 if (indxtovs < 0) cycle HEADER
683
684 iSensor = tvs_lsensor(indxTovs)
685
686 call obs_set_current_body_list(obsSpaceData, headerIndex)
687 iFov = obs_headElem_i(obsSpaceData, OBS_FOV, headerIndex)
688
689 if (bias(iSensor)%numScan > 1) then
690 iScan = iFov
691 else
692 iScan = 1
693 end if
694
695 BODY: do
696 bodyIndex = obs_getBodyIndex(obsSpaceData)
697 if (bodyIndex < 0) exit BODY
698
699 if (obs_bodyElem_i(obsSpaceData, OBS_ASS, bodyIndex) /= obs_assimilated) cycle BODY
700 if (obs_bodyElem_r(obsSpaceData, OBS_VAR, bodyIndex) == MPC_missingValue_R8) cycle BODY
701
702 call bcs_getChannelIndex(obsSpaceData, iSensor, chanIndx, bodyIndex)
703 if (chanindx > 0) then
704 if (bias(iSensor)%chans(chanIndx)%isDynamic .and. bias(iSensor)%numScan >0) then
705 if ( bias(iSensor)%chans(chanIndx)%coeff_fov(iScan)/= MPC_missingValue_R8 .and. &
706 all( bias(iSensor)%chans(chanIndx)%coeff(1:bias(iSensor)%chans(chanIndx)%NumActivePredictors)/= MPC_missingValue_R8) ) then
707 biasCor = bias(iSensor)%chans(chanIndx)%coeff_fov(iScan) + &
708 bias(iSensor)%chans(chanIndx)%coeff(1)
709 call bcs_getPredictors(predictor, headerIndex, iobs, chanIndx, obsSpaceData)
710 do iPredictor = 2, bias(iSensor)%chans(chanIndx)%NumActivePredictors
711 jPred = bias(iSensor)%chans(chanIndx)%PredictorIndex(iPredictor)
712 biasCor = biasCor + predictor(jPred) * bias(iSensor)%chans(chanIndx)%coeff(iPredictor)
713 end do
714 biasCor = -1.d0 * biascor
715 call obs_bodySet_r( obsSpaceData, OBS_BCOR, bodyIndex, biasCor)
716 end if
717 end if
718 end if
719 end do BODY
720 end do HEADER
721
722 write(*,*) "bcs_calcBias: end"
723
724 end subroutine bcs_calcBias
725
726 !-----------------------------------------------------------------------
727 ! bcs_dumpBiasToSqliteAfterThinning
728 !-----------------------------------------------------------------------
729 subroutine bcs_dumpBiasToSqliteAfterThinning(obsSpaceData)
730 !
731 !:Purpose: to dump bias correction coefficients and predictors in dedicated sqlite files
732 !
733 implicit none
734
735 ! Arguments:
736 type(struct_obs), intent(inout) :: obsSpaceData
737
738 ! Locals:
739 integer :: headerIndex, bodyIndex,iobs, indxtovs, idatyp
740 integer :: sensorIndex, iPredictor, chanIndx, codeTypeIndex, fileIndex, searchIndex
741 integer :: iScan, iFov, jPred, burpChanIndex
742 real(8) :: predictor(NumPredictors)
743 real(8) :: biasCor
744 real(8) :: sunzen, sunaz, satzen, sataz
745 type(fSQL_DATABASE), allocatable :: db(:) ! SQLIte file handle
746 type(fSQL_STATEMENT), allocatable :: stmtPreds(:), stmtCoeffs(:) ! precompiled SQLite statements
747 type(fSQL_STATUS) :: stat ! SQLiteerror status
748 character(len=512) :: queryCreate, queryPreds, queryCoeffs, queryTrim
749 logical, allocatable :: first(:)
750 integer, allocatable :: fileIndexes(:), obsOffset(:), dataOffset(:)
751 character(len=*), parameter :: myName = 'bcs_dumpBiasToSqliteAfterThinning::'
752 character(len=30) :: fileNameExtension
753 character(len=4) :: cmyidx, cmyidy
754 integer :: tovsCodeTypeListSize, tovsCodeTypeList(ninst)
755 integer :: tovsFileNameListSize
756 character(len=20) :: tovsFileNameList(30)
757 character(len=256) :: fileName
758 integer :: tovsAllCodeTypeListSize, tovsAllCodeTypeList(ninst)
759
760 if (.not. biasActive) return
761 if (.not. dumpToSqliteAfterThinning) return
762
763 write(*,*) "bcs_dumpBiasToSqliteAfterThinning: start"
764
765 ! get list of all possible tovs codetype values and unique list of corresponding filenames
766 call tvs_getAllIdBurpTovs(tovsAllCodeTypeListSize, tovsAllCodeTypeList)
767 write(*,*) 'bcs_dumpBiasToSqliteAfterThinning: tovsAllCodeTypeListSize = ', tovsAllCodeTypeListSize
768 write(*,*) 'bcs_dumpBiasToSqliteAfterThinning: tovsAllCodeTypeList = ', tovsAllCodeTypeList(1:tovsAllCodeTypeListSize)
769
770 tovsFileNameListSize = 0
771 tovsFileNameList(:) = 'XXXXX'
772 do codeTypeIndex = 1, tovsAllCodeTypeListSize
773 fileName = getObsFileName(tovsAllCodeTypeList(codeTypeIndex))
774 if (all(tovsFileNameList(:) /= fileName)) then
775 tovsFileNameListSize = tovsFileNameListSize + 1
776 tovsFileNameList(tovsFileNameListSize) = fileName
777 end if
778 end do
779 write(*,*) 'bcs_dumpBiasToSqliteAfterThinning: tovsFileNameListSize = ', tovsFileNameListSize
780 write(*,*) 'bcs_dumpBiasToSqliteAfterThinning: tovsFileNameList = ', tovsFileNameList(1:tovsFileNameListSize)
781
782 allocate(db(tovsFileNameListSize))
783 allocate(stmtPreds(tovsFileNameListSize))
784 allocate(stmtCoeffs(tovsFileNameListSize))
785 allocate(first(tovsFileNameListSize))
786 first(:) = .true.
787 allocate(fileIndexes(size(obsf_fileName)))
788 fileIndexes(:) = -1
789 do fileIndex = 1, tovsFileNameListSize
790 do searchIndex = 1, size(obsf_fileName)
791 if (index(trim(obsf_fileName(searchIndex)), trim(tovsFileNameList(fileIndex))) >0) then
792 fileIndexes(searchIndex) = fileIndex
793 end if
794 end do
795 end do
796 write(*,*) 'bcs_dumpBiasToSqliteAfterThinning: fileIndexes', fileIndexes(1:tovsFileNameListSize)
797 allocate(obsOffset(tovsFileNameListSize))
798 allocate(dataOffset(tovsFileNameListSize))
799 do fileIndex = 1, tovsFileNameListSize
800 fileName = tovsFileNameList(fileIndex)
801 write(*,*) 'tovs filename = ', fileName
802 ! get list of codetypes associated with this filename
803 tovsCodeTypeListSize = 0
804 tovsCodeTypeList(:) = MPC_missingValue_INT
805 do codeTypeIndex = 1, tovsAllCodeTypeListSize
806 if (fileName == getObsFileName(tovsAllCodeTypeList(codeTypeIndex))) then
807 tovsCodeTypeListSize = tovsCodeTypeListSize + 1
808 tovsCodeTypeList(tovsCodeTypeListSize) = tovsAllCodeTypeList(codeTypeIndex)
809 end if
810 end do
811
812 write(*,*) 'bcs_dumpBiasToSqliteAfterThinning: tovsCodeTypeListSize = ', tovsCodeTypeListSize
813 write(*,*) 'bcs_dumpBiasToSqliteAfterThinning: tovsCodeTypeList = ', tovsCodeTypeList(1:tovsCodeTypeListSize)
814 call getInitialIdObsData(obsSpaceData, 'TO', obsOffset(fileIndex), dataOffset(fileIndex), &
815 codeTypeList_opt=tovsCodeTypeList(1:tovsCodeTypeListSize))
816 write(*,*) 'bcs_dumpBiasToSqliteAfterThinning: obsOffset(fileIndex), dataOffset(fileIndex)', fileIndex, obsOffset(fileIndex), dataOffset(fileIndex)
817 end do
818
819 iobs = 0
820 call obs_set_current_header_list(obsSpaceData, 'TO')
821 HEADER: do
822 headerIndex = obs_getHeaderIndex(obsSpaceData)
823 if (headerIndex < 0) exit HEADER
824 ! process only radiance data to be assimilated?
825 idatyp = obs_headElem_i(obsSpaceData, OBS_ITY, headerIndex)
826 if (.not. tvs_isIdBurpTovs(idatyp)) then
827 write(*,*) 'bcs_dumpBiasToSqliteAfterThinning: warning unknown radiance codtyp present check NAMTOVSINST', idatyp
828 cycle HEADER
829 end if
830 iobs = iobs + 1
831 fileIndex = fileIndexes(obs_headElem_i(obsSpaceData, OBS_IDF, headerIndex))
832 indxtovs = tvs_tovsIndex(headerIndex)
833 if (indxtovs < 0) cycle HEADER
834 sensorIndex = tvs_lsensor(indxTovs)
835 if (first(fileIndex)) then
836 if (obs_mpiLocal(obsSpaceData)) then
837 write(cmyidy,'(i4.4)') (mmpi_myidy + 1)
838 write(cmyidx,'(i4.4)') (mmpi_myidx + 1)
839 fileNameExtension = trim(cmyidx) // '_' // trim(cmyidy)
840 else
841 fileNameExtension = ' '
842 end if
843
844 fileName = 'obs/bcr' // trim(tovsFileNameList(fileIndex)) &
845 // '_' // trim(filenameExtension)
846
847 call fSQL_open(db(fileIndex), fileName, stat)
848 write(*,*) 'bcs_dumpBiasToSqliteAfterThinning: Open ', trim(fileName)
849 if (fSQL_error(stat) /= FSQL_OK) call handleError(stat, 'fSQL_open: ')
850
851 ! Create the tables
852 queryCreate = 'CREATE TABLE predictors(id_data integer, id_obs integer, predIndex integer, PredictorValue real,' // &
853 'PredictorType varchar(2), bcor real, fov integer, sunzen real, sunaz real, satzen real, sataz real);' // &
854 'CREATE TABLE coeffs2(predIndex integer, coeff real, instrument varchar(10), platform varchar(16), vcoord integer, fov integer);'
855 call fSQL_do_many(db(fileIndex), queryCreate, stat)
856 if (fSQL_error(stat) /= FSQL_OK) call handleError(stat, 'fSQL_do_many: ')
857
858 queryPreds = 'insert into predictors (id_data, id_obs, predIndex, predictorValue,' // &
859 ' predictorType, bcor, fov, sunzen, sunaz, satzen, sataz) values(?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?);'
860 queryCoeffs = 'insert into coeffs2 (predIndex, coeff, instrument, platform, vcoord, fov) values(?, ?, ?, ?, ?, ?); '
861 write(*,*) myName // ' Insert query Predictors = ', trim(queryPreds)
862 write(*,*) myName // ' Insert query Coeffs = ', trim(queryCoeffs)
863 call fSQL_begin(db(fileIndex))
864 call fSQL_prepare(db(fileIndex), queryPreds, stmtPreds(fileIndex), stat)
865 if (fSQL_error(stat) /= FSQL_OK) call handleError(stat, 'fSQL_prepare 1: ')
866 call fSQL_prepare(db(fileIndex), queryCoeffs, stmtCoeffs(fileIndex), stat)
867 if (fSQL_error(stat) /= FSQL_OK) call handleError(stat, 'fSQL_prepare 2: ')
868 first(fileIndex) = .false.
869 end if
870 call obs_set_current_body_list(obsSpaceData, headerIndex)
871 iFov = obs_headElem_i(obsSpaceData, OBS_FOV, headerIndex)
872 sunzen = obs_headElem_r(obsSpaceData, OBS_SUN, headerIndex)
873 sunaz = obs_headElem_r(obsSpaceData, OBS_SAZ, headerIndex)
874 satzen = obs_headElem_r(obsSpaceData, OBS_SZA, headerIndex)
875 sataz = obs_headElem_r(obsSpaceData, OBS_AZA, headerIndex)
876 if (bias(sensorIndex)%numScan > 1) then
877 iScan = iFov
878 else
879 iScan = 1
880 end if
881 BODY: do
882 bodyIndex = obs_getBodyIndex(obsSpaceData)
883 if (bodyIndex < 0) exit BODY
884 if (obs_bodyElem_i(obsSpaceData, OBS_ASS, bodyIndex) /= obs_assimilated) cycle BODY
885 if (obs_bodyElem_r(obsSpaceData, OBS_VAR, bodyIndex) == MPC_missingValue_R8) cycle BODY
886 if (btest(obs_bodyElem_i(obsSpaceData, OBS_FLG, bodyIndex), 11)) cycle BODY
887 call bcs_getChannelIndex(obsSpaceData, sensorIndex, chanIndx, bodyIndex)
888 if (chanindx > 0) then
889 biasCor = 0.0d0
890 if (bias(sensorIndex)%chans(chanIndx)%isDynamic .and. bias(sensorIndex)%numScan > 0) then
891 call bcs_getPredictors(predictor, headerIndex, iobs, chanIndx, obsSpaceData)
892 biasCor = bias(sensorIndex)%chans(chanIndx)%coeff_fov(iScan) + &
893 bias(sensorIndex)%chans(chanIndx)%coeff(1)
894 burpChanIndex = nint(obs_bodyElem_r(obsSpaceData, OBS_PPP, bodyIndex))
895 call fSQL_bind_param(stmtCoeffs(fileIndex), param_index=1, int_var=0)
896 call fSQL_bind_param(stmtCoeffs(fileIndex), param_index=2, real8_var=bias(sensorIndex)%chans(chanIndx)%coeff_fov(iScan))
897 call fSQL_bind_param(stmtCoeffs(fileIndex), param_index=3, char_var=trim(tvs_instrumentName(sensorIndex)))
898 call fSQL_bind_param(stmtCoeffs(fileIndex), param_index=4, char_var=trim(tvs_satelliteName(sensorIndex)))
899 call fSQL_bind_param(stmtCoeffs(fileIndex), param_index=5, int_var=burpChanIndex)
900 call fSQL_bind_param(stmtCoeffs(fileIndex), param_index=6, int_var=iScan)
901 call fSQL_exec_stmt(stmtCoeffs(fileIndex))
902 call fSQL_bind_param(stmtCoeffs(fileIndex), param_index=1, int_var=1)
903 call fSQL_bind_param(stmtCoeffs(fileIndex), param_index=2, real8_var=bias(sensorIndex)%chans(chanIndx)%coeff(1))
904 call fSQL_bind_param(stmtCoeffs(fileIndex), param_index=3, char_var=trim(tvs_instrumentName(sensorIndex)))
905 call fSQL_bind_param(stmtCoeffs(fileIndex), param_index=4, char_var=trim(tvs_satelliteName(sensorIndex)))
906 call fSQL_bind_param(stmtCoeffs(fileIndex), param_index=5, int_var=burpChanIndex)
907 call fSQL_bind_param(stmtCoeffs(fileIndex), param_index=6, int_var=iScan)
908 call fSQL_exec_stmt(stmtCoeffs(fileIndex))
909 do iPredictor = 2, bias(sensorIndex)%chans(chanIndx)%NumActivePredictors
910 jPred = bias(sensorIndex)%chans(chanIndx)%PredictorIndex(iPredictor)
911 biasCor = biasCor + predictor(jPred) * bias(sensorIndex)%chans(chanIndx)%coeff(iPredictor)
912 end do
913 end if
914 biasCor = -1.d0 * biascor
915 do iPredictor = 2, bias(sensorIndex)%chans(chanIndx)%NumActivePredictors
916 jPred = bias(sensorIndex)%chans(chanIndx)%PredictorIndex(iPredictor)
917 call fSQL_bind_param(stmtCoeffs(fileIndex), param_index=1, int_var=jPred)
918 call fSQL_bind_param(stmtCoeffs(fileIndex), param_index=2, real8_var=bias(sensorIndex)%chans(chanIndx)%coeff(iPredictor))
919 call fSQL_bind_param(stmtCoeffs(fileIndex), param_index=3, char_var=trim(tvs_instrumentName(sensorIndex)))
920 call fSQL_bind_param(stmtCoeffs(fileIndex), param_index=4, char_var=trim(tvs_satelliteName(sensorIndex)))
921 call fSQL_bind_param(stmtCoeffs(fileIndex), param_index=5, int_var=burpChanIndex)
922 call fSQL_bind_param(stmtCoeffs(fileIndex), param_index=6, int_var=iScan)
923 call fSQL_exec_stmt (stmtCoeffs(fileIndex))
924 call fSQL_bind_param(stmtPreds(fileIndex), param_index=1, int_var=bodyIndex + dataOffset(fileIndex))
925 call fSQL_bind_param(stmtPreds(fileIndex), param_index=2, int_var=headerIndex + obsOffset(fileIndex))
926 call fSQL_bind_param(stmtPreds(fileIndex), param_index=3, int_var=jPred)
927 call fSQL_bind_param(stmtPreds(fileIndex), param_index=4, real8_var=predictor(jPred))
928 call fSQL_bind_param(stmtPreds(fileIndex), param_index=5, char_var=predtab(jPred))
929 call fSQL_bind_param(stmtPreds(fileIndex), param_index=6, real8_var=biascor)
930 call fSQL_bind_param(stmtPreds(fileIndex), param_index=7, int_var=iScan)
931 call fSQL_bind_param(stmtPreds(fileIndex), param_index=8, real8_var=sunzen)
932 call fSQL_bind_param(stmtPreds(fileIndex), param_index=9, real8_var=sunaz )
933 call fSQL_bind_param(stmtPreds(fileIndex), param_index=10, real8_var=satzen )
934 call fSQL_bind_param(stmtPreds(fileIndex), param_index=11, real8_var=sataz)
935 call fSQL_exec_stmt (stmtPreds(fileIndex))
936 end do
937 end if
938 end do BODY
939 end do HEADER
940 do fileIndex = 1, tovsFileNameListSize
941 if (.not. first(fileIndex)) then
942 call fSQL_finalize(stmtCoeffs(fileIndex))
943 call fSQL_finalize(stmtPreds(fileIndex))
944 queryTrim = 'create table coeffs as select distinct * from coeffs2; drop table coeffs2;'
945 call fSQL_do_many(db(fileIndex), queryTrim, stat)
946 if (fSQL_error(stat) /= FSQL_OK) call handleError(stat, 'fSQL_do_many: ')
947 call fSQL_commit(db(fileIndex), stat)
948 if (fSQL_error(stat) /= FSQL_OK) call handleError(stat, 'fSQL_commit: ')
949 call fSQL_close(db(fileIndex), stat)
950 if (fSQL_error(stat) /= FSQL_OK) call handleError(stat, 'fSQL_close: ')
951 end if
952 end do
953 deallocate(dataOffset)
954 deallocate(obsOffset)
955 deallocate(fileIndexes)
956 deallocate(first)
957 deallocate(stmtCoeffs)
958 deallocate(stmtPreds)
959 deallocate(db)
960 write(*,*) "bcs_dumpBiasToSqliteAfterThinning: end"
961 contains
962
963 subroutine handleError(stat, message)
964 implicit none
965
966 ! Arguments:
967 type(FSQL_STATUS), intent(in) :: stat
968 character(len=*), intent(in) :: message
969
970 write(*,*) message, fSQL_errmsg(stat)
971 call utl_abort(trim(message))
972 end subroutine handleError
973
974 end subroutine bcs_dumpBiasToSqliteAfterThinning
975
976 !---------------------------------------
977 ! bcs_computeResidualsStatistics
978 !----------------------------------------
979 subroutine bcs_computeResidualsStatistics(obsSpaceData, prefix)
980 !
981 !:Purpose: to compute residuals mean and standard deviation by intrument, channel and scan position
982 !
983 implicit none
984
985 ! Arguments:
986 type(struct_obs), intent(inout) :: obsSpaceData
987 character(len=*), intent(in) :: prefix
988
989 ! Locals:
990 real(8), allocatable :: tbias(:,:), tstd(:,:)
991 integer, allocatable :: tcount(:,:)
992 real(8), allocatable :: biasMpiGlobal(:,:), stdMpiGLobal(:,:)
993 integer, allocatable :: countMpiGlobal(:,:)
994 integer :: sensorIndex, headerIndex, bodyIndex
995 integer :: nchans, nscan
996 integer :: iSensor, iScan, chanIndx, iFov
997 real(8) :: OmF, bcor
998 integer :: ierr, nulfile1, nulfile2
999 character(len=10) :: instrName, satNamecoeff
1000 character(len=72) :: errorMessage
1001
1002 if (.not. biasActive) return
1003
1004 if (.not. outstats) return
1005
1006 write(*,*) "bcs_computeResidualsStatistics: start"
1007
1008 SENSORS:do sensorIndex = 1, tvs_nsensors
1009
1010 if (.not. tvs_isReallyPresentMpiGLobal(sensorIndex)) cycle SENSORS
1011
1012 write(*,*) "bcs_computeResidualsStatistics: sensorIndex ", sensorIndex
1013
1014 nchans = bias(sensorIndex)%numChannels
1015 nscan = bias(sensorIndex)%numscan
1016
1017 allocate(tbias(nchans,nscan))
1018 tbias(:,:) = 0.d0
1019 allocate(tstd(nchans,nscan))
1020 tstd(:,:) = 0.d0
1021 allocate(tcount(nchans,nscan))
1022 tcount(:,:) = 0
1023
1024 call obs_set_current_header_list(obsSpaceData, 'TO')
1025
1026 HEADER: do
1027 headerIndex = obs_getHeaderIndex(obsSpaceData)
1028 if (headerIndex < 0) exit HEADER
1029 if (tvs_tovsIndex(headerIndex) < 0) cycle HEADER
1030 iSensor = tvs_lsensor(tvs_tovsIndex(headerIndex))
1031 if (iSensor /= sensorIndex) cycle HEADER
1032
1033 iFov = obs_headElem_i(obsSpaceData, OBS_FOV, headerIndex)
1034 if (nscan > 1) then
1035 iScan = iFov
1036 else
1037 iScan = 1
1038 end if
1039
1040 call obs_set_current_body_list(obsSpaceData, headerIndex)
1041 BODY: do
1042 bodyIndex = obs_getBodyIndex(obsSpaceData)
1043 if (bodyIndex < 0) exit BODY
1044
1045 if (obs_bodyElem_i(obsSpaceData, OBS_ASS, bodyIndex) /= obs_assimilated) cycle BODY
1046 call bcs_getChannelIndex(obsSpaceData, iSensor, chanIndx, bodyIndex)
1047 if (chanindx > 0) then
1048 OmF = obs_bodyElem_r(obsSpaceData, OBS_OMP, bodyIndex)
1049 bcor = obs_bodyElem_r(obsSpaceData, OBS_BCOR, bodyIndex)
1050 if (OmF /= MPC_missingValue_R8 .and. bcor /= MPC_missingValue_R8) then
1051 tbias(chanIndx,iScan) = tbias(chanIndx,iScan) + (OmF + bcor)
1052 tstd(chanIndx,iScan) = tstd(chanIndx,iScan) + (OmF + bcor) ** 2
1053 tcount(chanIndx,iScan) = tcount(chanIndx,iScan) + 1
1054 end if
1055 end if
1056 end do BODY
1057 end do HEADER
1058
1059 allocate( biasMpiGlobal(nchans,nscan) )
1060 allocate( stdMpiGlobal(nchans,nscan) )
1061 allocate( countMpiGlobal(nchans,nscan) )
1062
1063 call mmpi_reduce_sumR8_2d( tbias, biasMpiGlobal, 0, "GRID" )
1064 call mmpi_reduce_sumR8_2d( tstd, stdMpiGlobal, 0, "GRID" )
1065 call rpn_comm_reduce(tcount, countMpiGlobal, size(countMpiGlobal), "mpi_integer", "MPI_SUM", 0, "GRID", ierr)
1066 if (ierr /=0) then
1067 write(errorMessage,*) "bcs_computeResidualsStatistics: MPI communication error 3", ierr
1068 call utl_abort(errorMessage)
1069 end if
1070
1071 if (mmpi_myId == 0) then
1072 where(countMpiGlobal > 0)
1073 biasMpiGlobal = biasMpiGlobal / countMpiGlobal
1074 stdMpiGlobal = sqrt(stdMpiGlobal/ countMpiGlobal - biasMpiGlobal**2)
1075 end where
1076
1077 instrName = InstrNametoCoeffFileName(tvs_instrumentName(sensorIndex))
1078 satNamecoeff = SatNameinCoeffFile(tvs_satelliteName(sensorIndex))
1079
1080 nulfile1 = 0
1081 ierr = fnom(nulfile1, './std_' // trim(instrName) // '_' // trim(satNamecoeff) // trim(prefix) // '.dat', 'FTN+FMT', 0)
1082 nulfile2 = 0
1083 ierr = fnom(nulfile2, './mean_' // trim(instrName) // '_' // trim(satNamecoeff) // trim(prefix) // '.dat', 'FTN+FMT', 0)
1084
1085 do chanIndx = 1, nchans
1086 if (sum(countMpiGlobal(chanIndx, :)) > 0) then
1087 write(nulfile2,'(i4,1x,100e14.6)') chanIndx, biasMpiGlobal(chanindx,:)
1088 write(nulfile1,'(i4,1x,100e14.6)') chanIndx, stdMpiGlobal(chanindx,:)
1089 end if
1090 end do
1091 ierr = fclos(nulfile1)
1092 ierr = fclos(nulfile2)
1093
1094 end if
1095
1096 deallocate(biasMpiGlobal)
1097 deallocate(stdMpiGlobal)
1098 deallocate(countMpiGlobal)
1099 deallocate(tbias)
1100 deallocate(tstd)
1101 deallocate(tcount)
1102
1103 end do SENSORS
1104
1105 write(*,*) "bcs_computeResidualsStatistics: end"
1106
1107 end subroutine bcs_computeResidualsStatistics
1108
1109 !---------------------------------------
1110 ! bcs_removeOutliers
1111 !----------------------------------------
1112 subroutine bcs_removeOutliers(obsSpaceData)
1113 !
1114 !:Purpose: to remove outliers (too large OmF) from linear regression
1115 !
1116 implicit none
1117
1118 ! Arguments:
1119 type(struct_obs), intent(inout) :: obsSpaceData
1120
1121 ! Locals:
1122 real(8), allocatable :: tbias(:,:), tstd(:,:)
1123 integer, allocatable :: tcount(:,:)
1124 real(8), allocatable :: biasMpiGlobal(:,:), stdMpiGLobal(:,:)
1125 integer, allocatable :: countMpiGlobal(:,:)
1126 integer :: sensorIndex, headerIndex, bodyIndex
1127 integer :: nchans, nfiles
1128 integer :: iSensor, chanIndx, timeIndex
1129 real(8) :: OmF
1130 real(8), parameter :: alpha = 5.d0
1131 real(8) :: stepObsIndex
1132 integer :: ierr
1133 character(len=72) :: errorMessage
1134
1135 if (.not. biasActive) return
1136
1137 write(*,*) "bcs_removeOutliers: start"
1138
1139 SENSORS:do sensorIndex = 1, tvs_nsensors
1140
1141 if (.not. tvs_isReallyPresentMpiGLobal(sensorIndex)) cycle SENSORS
1142
1143 write(*,*) "bcs_removeOutliers: sensorIndex ", sensorIndex
1144
1145 nchans = bias(sensorIndex)%numChannels
1146 nfiles = tim_nstepobs
1147
1148 allocate(tbias(nchans,nfiles))
1149 tbias(:,:) = 0.d0
1150 allocate(tstd(nchans,nfiles))
1151 tstd(:,:) = 0.d0
1152 allocate(tcount(nchans,nfiles))
1153 tcount(:,:) = 0
1154
1155 call obs_set_current_header_list(obsSpaceData, 'TO')
1156
1157 HEADER1: do
1158 headerIndex = obs_getHeaderIndex(obsSpaceData)
1159 if (headerIndex < 0) exit HEADER1
1160 if (tvs_tovsIndex(headerIndex) < 0) cycle HEADER1
1161 iSensor = tvs_lsensor(tvs_tovsIndex(headerIndex))
1162 if (iSensor /= sensorIndex) cycle HEADER1
1163
1164 call tim_getStepObsIndex(stepObsIndex, tim_getDatestamp(), &
1165 obs_headElem_i(obsSpaceData, OBS_DAT, headerIndex), &
1166 obs_headElem_i(obsSpaceData, OBS_ETM, headerIndex), tim_nstepobs)
1167
1168 timeIndex = nint(stepObsIndex)
1169 if (timeIndex < 0) cycle HEADER1
1170 call obs_set_current_body_list(obsSpaceData, headerIndex)
1171 BODY1: do
1172 bodyIndex = obs_getBodyIndex(obsSpaceData)
1173 if (bodyIndex < 0) exit BODY1
1174
1175 if (obs_bodyElem_i(obsSpaceData, OBS_ASS, bodyIndex) == obs_assimilated .and. &
1176 .not. btest(obs_bodyElem_i(obsSpaceData, OBS_FLG, bodyIndex), 6)) then
1177 call bcs_getChannelIndex(obsSpaceData, iSensor, chanIndx, bodyIndex)
1178 if (chanindx > 0) then
1179 OmF = obs_bodyElem_r(obsSpaceData, OBS_OMP, bodyIndex)
1180 if (OmF /= MPC_missingValue_R8) then
1181 tbias(chanIndx,timeindex) = tbias(chanIndx,timeindex) + OmF
1182 tstd(chanIndx,timeindex) = tstd(chanIndx,timeindex) + (OmF ** 2)
1183 tcount(chanIndx,timeindex) = tcount(chanIndx,timeindex) + 1
1184 end if
1185 end if
1186 end if
1187 end do BODY1
1188 end do HEADER1
1189
1190 allocate(biasMpiGlobal(nchans,nfiles))
1191 allocate(stdMpiGlobal(nchans,nfiles))
1192 allocate(countMpiGlobal(nchans,nfiles))
1193
1194 call mmpi_reduce_sumR8_2d( tbias, biasMpiGlobal, 0, "GRID" )
1195 call mmpi_reduce_sumR8_2d( tstd, stdMpiGlobal, 0, "GRID" )
1196 call rpn_comm_reduce(tcount, countMpiGlobal, size(countMpiGlobal), "mpi_integer", "MPI_SUM", 0, "GRID", ierr)
1197 if (ierr /=0) then
1198 write(errorMessage,*) "bcs_removeOutliers: MPI communication error 3", ierr
1199 call utl_abort(errorMessage)
1200 end if
1201
1202 if (mmpi_myId == 0) then
1203 where(countMpiGlobal > 0)
1204 biasMpiGlobal = biasMpiGlobal / countMpiGlobal
1205 stdMpiGlobal = sqrt(stdMpiGlobal/ countMpiGlobal - biasMpiGlobal**2)
1206 end where
1207 end if
1208
1209 call rpn_comm_bcast(countMpiGlobal, nchans*nfiles, "mpi_integer", 0, "GRID", ierr)
1210 if (ierr /=0) then
1211 write(errorMessage,*) "bcs_removeOutliers: MPI communication error 4", ierr
1212 call utl_abort(errorMessage)
1213 end if
1214 call rpn_comm_bcast(stdMpiGlobal, nchans*nfiles, "mpi_double_precision", 0, "GRID", ierr)
1215 if (ierr /=0) then
1216 write(errorMessage,*) "bcs_removeOutliers: MPI communication error 5", ierr
1217 call utl_abort(errorMessage)
1218 end if
1219
1220 if (sum(countMpiGlobal) /= 0) then
1221
1222 tcount(:,:) = 0
1223
1224 call obs_set_current_header_list(obsSpaceData, 'TO')
1225
1226 HEADER2: do
1227 headerIndex = obs_getHeaderIndex(obsSpaceData)
1228 if (headerIndex < 0) exit HEADER2
1229 if (tvs_tovsIndex(headerIndex) < 0) cycle HEADER2
1230 iSensor = tvs_lsensor(tvs_tovsIndex(headerIndex))
1231 if (iSensor /= sensorIndex) cycle HEADER2
1232
1233 call tim_getStepObsIndex(stepObsIndex, tim_getDatestamp(), &
1234 obs_headElem_i(obsSpaceData, OBS_DAT, headerIndex), &
1235 obs_headElem_i(obsSpaceData, OBS_ETM, headerIndex), tim_nstepobs)
1236
1237 timeIndex = nint(stepObsIndex)
1238 if (timeIndex < 0) cycle HEADER2
1239
1240 call obs_set_current_body_list(obsSpaceData, headerIndex)
1241 BODY2: do
1242 bodyIndex = obs_getBodyIndex(obsSpaceData)
1243 if (bodyIndex < 0) exit BODY2
1244
1245 if (obs_bodyElem_i(obsSpaceData, OBS_ASS, bodyIndex) == obs_assimilated .and. &
1246 .not. btest(obs_bodyElem_i(obsSpaceData, OBS_FLG, bodyIndex), 6)) then
1247
1248 call bcs_getChannelIndex(obsSpaceData, iSensor, chanIndx, bodyIndex)
1249
1250 if (chanindx > 0) then
1251
1252 OmF = obs_bodyElem_r(obsSpaceData, OBS_OMP, bodyIndex)
1253 if (countMpiGlobal(chanIndx, timeindex) > 2 .and. &
1254 abs(OmF) > alpha * stdMpiGlobal(chanIndx, timeindex)) then
1255 call obs_bodySet_i(obsSpaceData, OBS_ASS, bodyIndex, obs_notAssimilated)
1256 tcount(chanIndx,timeindex) = tcount(chanIndx,timeindex) + 1
1257 end if
1258
1259 end if
1260 end if
1261 end do BODY2
1262 end do HEADER2
1263
1264 do chanIndx = 1, nchans
1265 if (sum(tcount(chanIndx,:)) > 0) then
1266 write(*,'(A,1x,i2,1x,i4,1x,i6,1x,30e14.6)') "bcs_removeOutliers:", sensorindex, chanIndx, sum(tcount(chanIndx,:)), stdMpiGlobal(chanIndx,:)
1267 end if
1268 end do
1269
1270 end if
1271
1272 deallocate(biasMpiGlobal)
1273 deallocate(stdMpiGlobal)
1274 deallocate(countMpiGlobal)
1275 deallocate(tbias)
1276 deallocate(tstd)
1277 deallocate(tcount)
1278
1279 end do SENSORS
1280
1281 write(*,*) "bcs_removeOutliers: end"
1282
1283 end subroutine bcs_removeOutliers
1284
1285 !---------------------------------------
1286 ! bcs_calcBias_tl
1287 !----------------------------------------
1288 subroutine bcs_calcBias_tl(cv_in, obsColumnIndex, obsSpaceData, columnTrlOnTrlLev)
1289 !
1290 !:Purpose: tl of bias computation (for varBC)
1291 !
1292 implicit none
1293
1294 ! Arguments:
1295 real(8), intent(in) :: cv_in(:)
1296 integer, intent(in) :: obsColumnIndex
1297 type(struct_obs), intent(inout) :: obsSpaceData
1298 type(struct_columnData), intent(inout) :: columnTrlOnTrlLev
1299
1300 ! Locals:
1301 integer :: headerIndex, bodyIndex, iobs, indxtovs, idatyp
1302 integer :: iSensor, iPredictor, chanIndx
1303 integer :: iScan, iFov, jPred
1304 real(8) :: predictor(NumPredictors)
1305 real(8), pointer :: cv_bias(:)
1306 real(8), target :: dummy4Pointer(1)
1307 real(8) :: biasCor
1308
1309 if (.not. biasActive) return
1310
1311 if (.not. allocated(trialHeight300m1000)) then
1312 call bcs_getTrialPredictors(obsSpaceData, columnTrlOnTrlLev)
1313 call bcs_computePredictorBiases(obsSpaceData)
1314 call bcs_getRadiosondeWeight(obsSpaceData, lmodify_obserror_opt=.true.)
1315 end if
1316
1317 nullify(cv_bias)
1318 if (mmpi_myid == 0) then
1319 if (cvm_subVectorExists('BIAS')) then
1320 cv_Bias => cvm_getSubVector(cv_in, 'BIAS')
1321 write(*,*) 'bcs_calcBias_tl: maxval(cv_bias) = ', maxval(cv_bias(:))
1322 else
1323 write(*,*) 'bcs_calcBias_tl: control vector does not include bias coefficients'
1324 return
1325 end if
1326 else
1327 cv_bias => dummy4Pointer
1328 end if
1329
1330 ! get bias coefficients
1331 call bcs_cvToCoeff(cv_bias)
1332
1333 ! apply bias increment to specified obs column
1334 iobs = 0
1335 call obs_set_current_header_list(obsSpaceData, 'TO')
1336 HEADER: do
1337 headerIndex = obs_getHeaderIndex(obsSpaceData)
1338 if (headerIndex < 0) exit HEADER
1339
1340 ! process only radiance data to be assimilated?
1341 idatyp = obs_headElem_i(obsSpaceData, OBS_ITY, headerIndex)
1342 if (.not. tvs_isIdBurpTovs(idatyp)) then
1343 write(*,*) 'bcs_calBias_tl: warning unknown radiance codtyp present check NAMTOVSINST', idatyp
1344 cycle HEADER
1345 end if
1346
1347 iobs = iobs + 1
1348
1349 indxtovs = tvs_tovsIndex(headerIndex)
1350 if (indxtovs < 0) cycle HEADER
1351
1352 iSensor = tvs_lsensor(indxTovs)
1353
1354 call obs_set_current_body_list(obsSpaceData, headerIndex)
1355 iFov = obs_headElem_i(obsSpaceData, OBS_FOV, headerIndex)
1356
1357 BODY: do
1358 bodyIndex = obs_getBodyIndex(obsSpaceData)
1359 if (bodyIndex < 0) exit BODY
1360
1361 if (obs_bodyElem_i(obsSpaceData, OBS_ASS, bodyIndex) /= obs_assimilated) cycle BODY
1362
1363 call bcs_getChannelIndex(obsSpaceData, iSensor, chanIndx, bodyIndex)
1364
1365 if (chanindx > 0) then
1366 biasCor = 0.0d0
1367 if (bias(iSensor)%chans(chanIndx)%isDynamic .and. bias(iSensor)%numScan >0) then
1368 call bcs_getPredictors(predictor, headerIndex, iobs, chanindx, obsSpaceData)
1369 do iPredictor = 1, bias(iSensor)%chans(chanIndx)%NumActivePredictors
1370 jPred = bias(iSensor)%chans(chanIndx)%PredictorIndex(iPredictor)
1371 if (iPredictor == 1) then
1372 if (bias(iSensor)%numScan > 1) then
1373 iScan = iFov
1374 else
1375 iScan = 1
1376 end if
1377 biasCor = biasCor + predictor(jPred) * bias(iSensor)%chans(chanIndx)%coeffIncr_fov(iScan)
1378 else
1379 biasCor = biasCor + predictor(jPred) * bias(iSensor)%chans(chanIndx)%coeffIncr(iPredictor)
1380 end if
1381 end do
1382 end if
1383 call obs_bodySet_r(obsSpaceData, obsColumnIndex, bodyIndex, &
1384 obs_bodyElem_r(obsSpaceData, obsColumnIndex, bodyIndex) - biasCor)
1385 end if
1386 end do BODY
1387 end do HEADER
1388
1389
1390 end subroutine bcs_calcBias_tl
1391
1392 !----------------------
1393 ! bcs_getTrialPredictors
1394 !----------------------
1395 subroutine bcs_getTrialPredictors(obsSpaceData, columnTrlOnTrlLev)
1396 !
1397 !:Purpose: get predictors from trial fields
1398 !
1399 implicit none
1400
1401 ! Arguments:
1402 type(struct_columnData), intent(inout) :: columnTrlOnTrlLev
1403 type(struct_obs), intent(inout) :: obsSpaceData
1404
1405 ! Locals:
1406 integer :: headerIndex, idatyp, iobs
1407 real(8) :: height1, height2
1408
1409 if (tvs_nobtov > 0) then
1410 allocate(trialHeight300m1000(tvs_nobtov))
1411 allocate(trialHeight50m200(tvs_nobtov))
1412 allocate(trialHeight5m50(tvs_nobtov))
1413 allocate(trialHeight1m10(tvs_nobtov))
1414 allocate(trialTG(tvs_nobtov))
1415 allocate(RadiosondeWeight(tvs_nobtov))
1416 else
1417 write(*,*) 'bcs_getTrialPredictors: No radiance OBS found'
1418 return
1419 end if
1420
1421 iobs = 0
1422
1423 call obs_set_current_header_list(obsSpaceData, 'TO')
1424
1425 HEADER2: do
1426 headerIndex = obs_getHeaderIndex(obsSpaceData)
1427 if (headerIndex < 0) exit HEADER2
1428 idatyp = obs_headElem_i(obsSpaceData, OBS_ITY, headerIndex)
1429 if (.not. tvs_isIdBurpTovs(idatyp)) then
1430 write(*,*) 'bcs_getTrialPredictors: warning unknown radiance codtyp present check NAMTOVSINST', idatyp
1431 cycle HEADER2
1432 end if
1433 iobs = iobs + 1
1434
1435 height1 = logInterpHeight(columnTrlOnTrlLev, headerIndex, 1000.d0)
1436 height2 = logInterpHeight(columnTrlOnTrlLev, headerIndex, 300.d0)
1437
1438 trialHeight300m1000(iobs) = height2 - height1
1439
1440 height1 = logInterpHeight(columnTrlOnTrlLev, headerIndex, 200.d0)
1441 height2 = logInterpHeight(columnTrlOnTrlLev, headerIndex, 50.d0)
1442
1443 trialHeight50m200(iobs) = height2 - height1
1444
1445 height1 = height2
1446 height2 = logInterpHeight(columnTrlOnTrlLev, headerIndex, 5.d0)
1447
1448 trialHeight5m50(iobs) = height2 - height1
1449
1450 height1 = logInterpHeight(columnTrlOnTrlLev, headerIndex, 10.d0)
1451 height2 = logInterpHeight(columnTrlOnTrlLev, headerIndex, 1.d0)
1452
1453 trialHeight1m10(iobs) = height2 - height1
1454
1455 trialTG(iobs) = col_getElem(columnTrlOnTrlLev, 1, headerIndex, 'TG')
1456
1457 end do HEADER2
1458
1459 if (trialTG(1) > 150.0d0) then
1460 write(*,*) 'bcs_getTrialPredictors: converting TG from Kelvin to deg_C'
1461 trialTG(:) = trialTG(:) - MPC_K_C_DEGREE_OFFSET_R8
1462 end if
1463
1464 trialHeight300m1000(:) = 0.1d0 * trialHeight300m1000(:) ! conversion factor
1465 trialHeight50m200(:) = 0.1d0 * trialHeight50m200(:)
1466 trialHeight5m50(:) = 0.1d0 * trialHeight5m50(:)
1467 trialHeight1m10(:) = 0.1d0 * trialHeight1m10(:)
1468
1469 write(*,*) 'bcs_getTrialPredictors: end'
1470
1471 contains
1472
1473 function logInterpHeight(columnTrlOnTrlLev, headerIndex, p) result(height)
1474 implicit none
1475
1476 ! Arguments:
1477 type(struct_columnData), intent(inout) :: columnTrlOnTrlLev
1478 integer, intent(in) :: headerIndex
1479 real(8), intent(in) :: p
1480 ! Result:
1481 real(8) :: height
1482
1483 ! Locals:
1484 integer :: jk, nlev, ik
1485 real(8) :: zpt, zpb, zwt, zwb
1486 real(8), pointer :: col_ptr(:)
1487
1488 ik = 1
1489 nlev = col_getNumLev(columnTrlOnTrlLev, 'TH')
1490 do jk = 2, nlev - 1
1491 zpt = col_getPressure(columnTrlOnTrlLev, jk, headerIndex, 'TH') * MPC_MBAR_PER_PA_R8
1492 if(p > zpt) ik = jk
1493 end do
1494 zpt = col_getPressure(columnTrlOnTrlLev, ik, headerIndex, 'TH') * MPC_MBAR_PER_PA_R8
1495 zpb = col_getPressure(columnTrlOnTrlLev, ik + 1, headerIndex, 'TH') * MPC_MBAR_PER_PA_R8
1496
1497 zwb = log(p/zpt) / log(zpb/zpt)
1498 zwt = 1.d0 - zwb
1499 col_ptr => col_getColumn(columnTrlOnTrlLev, headerIndex, 'Z_T')
1500
1501 height = zwb * col_ptr(ik+1) + zwt * col_ptr(ik)
1502
1503 end function logInterpHeight
1504
1505 end subroutine bcs_getTrialPredictors
1506
1507 !----------------------
1508 ! bcs_cvToCoeff
1509 !----------------------
1510 subroutine bcs_cvToCoeff(cv_bias)
1511 !
1512 !:Purpose: get coefficient increment from control vector
1513 !
1514 implicit none
1515
1516 ! Arguments:
1517 real(8), intent(in) :: cv_bias(:)
1518
1519 ! Locals:
1520 integer :: index_cv, iSensor, iChannel, iPredictor, iScan
1521 integer :: nsize, ierr
1522
1523 if (mmpi_myid == 0) then
1524 write(*,*) 'bcs_cvToCoeff: start'
1525 index_cv = 0
1526 ! initialize of coeffIncr
1527 do iSensor = 1, tvs_nSensors
1528 if (bias(iSensor)%numScan > 0) then
1529 do iChannel = 1, bias(iSensor)%numChannels
1530 if (bias(iSensor)%chans(iChannel)%isDynamic) then
1531 bias(iSensor)%chans(iChannel)%coeffIncr(:) = 0.0d0
1532 bias(iSensor)%chans(iChannel)%coeffIncr_fov(:) = 0.0d0
1533 end if
1534 end do
1535 end if
1536 end do
1537
1538 do iSensor = 1, tvs_nSensors
1539 if (bias(iSensor)%numScan > 0) then
1540 do iChannel = 1, bias(iSensor)%numChannels
1541 if (bias(iSensor)%chans(iChannel)%isDynamic) then
1542 do iPredictor = 1, bias(iSensor)%chans(iChannel)%numActivePredictors
1543 if (iPredictor == 1) then
1544 do iScan = 1, bias(iSensor)%numScan
1545 index_cv = index_cv + 1
1546 bias(iSensor)%chans(iChannel)%coeffIncr_fov(iScan) = bias(iSensor)%chans(iChannel)%stddev(iPredictor) * cv_bias(index_cv)
1547 end do
1548 else
1549 index_cv = index_cv + 1
1550 bias(iSensor)%chans(iChannel)%coeffIncr(iPredictor) = bias(iSensor)%chans(iChannel)%stddev(iPredictor) * cv_bias(index_cv)
1551 end if
1552 end do !iPredictor
1553 end if ! isDynamic
1554 end do !iChannel
1555 end if
1556 end do !iSensor
1557
1558 end if
1559
1560 ! for constant part
1561 do iSensor = 1, tvs_nSensors
1562 if (bias(iSensor)%numScan > 0) then
1563 nsize = bias(iSensor)%numScan
1564 do iChannel = 1, bias(iSensor)%numChannels
1565 if (bias(iSensor)%chans(iChannel)%isDynamic) &
1566 call rpn_comm_bcast(bias(iSensor)%chans(iChannel)%coeffIncr_fov, nsize, "mpi_double_precision", 0, "GRID", ierr)
1567 end do
1568 end if
1569 end do
1570
1571 ! for predictor part
1572 do iSensor = 1, tvs_nSensors
1573 if (bias(iSensor)%numScan > 0) then
1574 do iChannel = 1, bias(iSensor)%numChannels
1575 if (bias(iSensor)%chans(iChannel)%isDynamic) then
1576 nsize = bias(iSensor)%chans(iChannel)%numActivePredictors
1577 call rpn_comm_bcast(bias(iSensor)%chans(iChannel)%coeffIncr, nsize, "mpi_double_precision", 0, "GRID", ierr)
1578 end if
1579 end do
1580 end if
1581 end do
1582
1583 end subroutine bcs_cvToCoeff
1584
1585 !-----------------------------------
1586 ! bcs_getPredictors
1587 !----------------------------------
1588 subroutine bcs_getPredictors(predictor, headerIndex, obsIndex, chanindx, obsSpaceData)
1589 !
1590 !:Purpose: get predictors
1591 !
1592 implicit none
1593
1594 ! Arguments:
1595 real(8), intent(out) :: predictor(NumPredictors)
1596 integer, intent(in) :: headerIndex
1597 integer, intent(in) :: obsIndex
1598 integer, intent(in) :: chanindx
1599 type(struct_obs), intent(inout) :: obsSpaceData
1600
1601 ! Locals:
1602 integer :: iSensor, iPredictor, jPredictor
1603 real(8) :: zenithAngle
1604
1605 predictor(:) = 0.0d0
1606
1607 do iPredictor = 1, NumPredictors
1608
1609 if (iPredictor == 1) then
1610 ! constant
1611 predictor(iPredictor) = 1.0d0
1612 else if (iPredictor == 2) then
1613 ! Height300-Height1000 (dam) /1000 T1
1614 predictor(iPredictor) = trialHeight300m1000(obsIndex) / 1000.0d0
1615 else if (iPredictor == 3) then
1616 ! Height50-Height200 (dam) /1000 T2
1617 predictor(iPredictor) = trialHeight50m200(obsIndex) / 1000.0d0
1618 else if (iPredictor == 4) then
1619 ! Height5-Height50 (dam) /1000 T3
1620 predictor(iPredictor) = trialHeight5m50(obsIndex) / 1000.0d0
1621 else if (iPredictor == 5) then
1622 ! Height1-Height10 (dam) /1000 T4
1623 predictor(iPredictor) = trialHeight1m10(obsIndex) / 1000.0d0
1624 else if (iPredictor == 6) then
1625 ! SV secant of satellite zenith angle minus one
1626 zenithAngle = obs_headElem_r(obsSpaceData, OBS_SZA, headerIndex)
1627 if (zenithAngle < 75.) predictor(iPredictor) = (1.d0 / cos(zenithAngle * MPC_RADIANS_PER_DEGREE_R8)) - 1.d0
1628 else if (iPredictor == 7) then
1629 ! skin temperature (C) /10
1630 predictor(iPredictor) = trialTG(obsIndex)
1631 end if
1632
1633 end do
1634
1635 iSensor = tvs_lsensor(tvs_tovsIndex(headerIndex))
1636 do iPredictor = 1, bias(iSensor)%chans(chanIndx)%numActivePredictors
1637 jPredictor = bias(iSensor)%chans(chanIndx)%predictorIndex(iPredictor)
1638 predictor(jPredictor) = predictor(jPredictor) - bias(iSensor)%chans(chanindx)%coeff_offset(iPredictor)
1639 end do
1640
1641 end subroutine bcs_getPredictors
1642
1643 !---------------------------------------------
1644 ! bcs_calcBias_ad
1645 !---------------------------------------------
1646 subroutine bcs_calcBias_ad(cv_out, obsColumnIndex, obsSpaceData)
1647 !
1648 !:Purpose: bias computation adjoint (for varBC)
1649 !
1650 implicit none
1651
1652 ! Arguments:
1653 real(8), intent(in) :: cv_out(:)
1654 integer, intent(in) :: obsColumnIndex
1655 type(struct_obs), intent(inout) :: obsSpaceData
1656
1657 ! Locals:
1658 integer :: headerIndex, bodyIndex, iobs, idatyp
1659 integer :: iSensor, iChannel, iPredictor, chanIndx
1660 integer :: iScan, iFOV, jPred
1661 real(8) :: predictor(NumPredictors)
1662 real(8), pointer :: cv_bias(:)
1663 real(8), target :: dummy4Pointer(1)
1664 real(8) :: biasCor
1665
1666 if (.not. biasActive) return
1667
1668 if (mmpi_myid == 0) write(*,*) 'bcs_calcBias_ad: start'
1669
1670 nullify(cv_bias)
1671 if (mmpi_myid == 0) then
1672 if (cvm_subVectorExists('BIAS')) then
1673 cv_bias => cvm_getSubVector(cv_out, 'BIAS')
1674 else
1675 write(*,*) 'bcs_calcBias_ad: control vector does not include bias coefficients'
1676 return
1677 end if
1678 else
1679 cv_bias => dummy4Pointer
1680 end if
1681
1682 ! adjoint of applying bias increment to specified obs column
1683 do iSensor = 1, tvs_nSensors
1684 if (bias(iSensor)%numScan > 0) then
1685 do iChannel = 1, bias(iSensor)%numChannels
1686 if (bias(iSensor)%chans(iChannel)%isDynamic) then
1687 bias(iSensor)%chans(iChannel)%coeffIncr(:) = 0.0d0
1688 bias(iSensor)%chans(iChannel)%coeffIncr_fov(:) = 0.0d0
1689 end if
1690 end do
1691 end if
1692 end do
1693
1694 iobs = 0
1695 call obs_set_current_header_list(obsSpaceData, 'TO')
1696 HEADER: do
1697 headerIndex = obs_getHeaderIndex(obsSpaceData)
1698 if (headerIndex < 0) exit HEADER
1699 idatyp = obs_headElem_i(obsSpaceData, OBS_ITY, headerIndex)
1700 if (.not. tvs_isIdBurpTovs(idatyp)) then
1701 write(*,*) 'bcs_calcBias_ad: warning unknown radiance codtyp present check NAMTOVSINST', idatyp
1702 cycle HEADER
1703 end if
1704 iobs = iobs + 1
1705 if (tvs_tovsIndex(headerIndex) < 0) cycle HEADER
1706
1707 iSensor = tvs_lsensor(tvs_tovsIndex(headerIndex))
1708
1709 call obs_set_current_body_list(obsSpaceData, headerIndex)
1710 iFov = obs_headElem_i(obsSpaceData, OBS_FOV, headerIndex)
1711
1712 BODY: do
1713 bodyIndex = obs_getBodyIndex(obsSpaceData)
1714 if (bodyIndex < 0) exit BODY
1715
1716 if (obs_bodyElem_i(obsSpaceData, OBS_ASS, bodyIndex) /= obs_assimilated) cycle BODY
1717
1718 call bcs_getChannelIndex(obsSpaceData, iSensor, chanIndx, bodyIndex)
1719
1720 if (chanindx > 0) then
1721 if (bias(iSensor)%chans(chanIndx)%isDynamic) then
1722 call bcs_getPredictors(predictor, headerIndex, iobs, chanIndx, obsSpaceData)
1723 biasCor = obs_bodyElem_r(obsSpaceData, obsColumnIndex, bodyIndex)
1724 do iPredictor = 1, bias(iSensor)%chans(chanIndx)%numActivePredictors
1725 jPred = bias(iSensor)%chans(chanIndx)%PredictorIndex(iPredictor)
1726 if (jPred == 1) then
1727 if (bias(iSensor)%numScan > 1) then
1728 iScan = iFov
1729 else
1730 iScan = 1
1731 end if
1732 bias(iSensor)%chans(chanIndx)%coeffIncr_fov(iScan) = bias(iSensor)%chans(chanIndx)%coeffIncr_fov(iScan) &
1733 + predictor(jPred) * biasCor
1734 else
1735 bias(iSensor)%chans(chanIndx)%coeffIncr(iPredictor) = bias(iSensor)%chans(chanIndx)%coeffIncr(iPredictor) &
1736 + predictor(jPred) * biasCor
1737 end if
1738 end do !iPredictor
1739 end if
1740 end if
1741
1742 end do BODY
1743
1744 end do HEADER
1745
1746 ! put the coefficients into the control vector
1747 call bcs_cvToCoeff_ad(cv_bias)
1748
1749 if (mmpi_myid == 0) then
1750 write(*,*) 'bcs_calcBias_ad: maxval(cv_bias) = ', maxval(cv_bias(:))
1751 end if
1752
1753 end subroutine bcs_calcBias_ad
1754
1755 !----------------------------------------------------
1756 ! bcs_cvToCoeff_ad
1757 !----------------------------------------------------
1758 subroutine bcs_cvToCoeff_ad(cv_bias)
1759 !
1760 !:Purpose: adjoint of control vector to coeff transfer (for varBC)
1761 !
1762 implicit none
1763
1764 ! Arguments:
1765 real(8), intent(inout) :: cv_bias(:)
1766
1767 ! Locals:
1768 integer :: index_cv, iSensor, iChannel, iPredictor, iScan
1769 integer :: nChan, nScan
1770 integer :: nsize, iChan
1771 real(8), allocatable :: temp_coeffIncr(:), temp_coeffIncr_fov(:)
1772
1773 write(*,*) "bcs_cvToCoeff_ad: start"
1774
1775 do iSensor = 1, tvs_nSensors
1776 if (bias(iSensor)%numScan > 0) then
1777 nChan = bias(iSensor)%numChannels
1778 do ichan = 1, nChan
1779 if (bias(iSensor)%chans(ichan)%isDynamic) then
1780 nSize = bias(iSensor)%chans(iChan)%numActivePredictors
1781 allocate(temp_coeffIncr(nSize))
1782 temp_coeffIncr(:) = 0.0d0
1783 call mmpi_reduce_sumR8_1d( bias(iSensor)%chans(ichan)%coeffIncr(:), temp_coeffIncr, 0, "GRID" )
1784 bias(iSensor)%chans(ichan)%coeffIncr(:) = temp_coeffIncr(:)
1785 deallocate(temp_coeffIncr)
1786 end if
1787 end do
1788 end if
1789 end do
1790
1791 do iSensor = 1, tvs_nSensors
1792 if (bias(iSensor)%numScan > 0) then
1793 nChan = bias(iSensor)%numChannels
1794 nScan = bias(iSensor)%numScan
1795 nsize = nChan * nScan
1796 if (nsize > 0) then
1797 allocate(temp_coeffIncr_fov(1:nScan))
1798 do ichan = 1, nChan
1799 if (bias(iSensor)%chans(ichan)%isDynamic) then
1800 temp_coeffIncr_fov(:) = 0.0d0
1801 call mmpi_reduce_sumR8_1d( bias(iSensor)%chans(ichan)%coeffIncr_fov, temp_coeffIncr_fov, 0, "GRID" )
1802 bias(iSensor)%chans(iChan)%coeffIncr_fov(:) = temp_coeffIncr_fov(:)
1803 end if
1804 end do
1805 deallocate(temp_coeffIncr_fov)
1806 end if
1807 end if
1808 end do
1809
1810 if (mmpi_myid == 0) then
1811 cv_bias(:) = 0.d0
1812 index_cv = 0
1813 do iSensor = 1, tvs_nSensors
1814 if (bias(iSensor)%numScan > 0) then
1815 do iChannel = 1, bias(iSensor)%numChannels
1816 if (bias(iSensor)%chans(iChannel)%isDynamic) then
1817 do iPredictor = 1, bias(iSensor)%chans(iChannel)%numActivePredictors
1818 if (iPredictor == 1) then
1819 do iScan = 1, bias(iSensor)%numScan
1820 index_cv = index_cv + 1
1821 cv_bias(index_cv) = bias(iSensor)%chans(iChannel)%stddev(iPredictor) * bias(iSensor)%chans(iChannel)%coeffIncr_fov(iScan)
1822 end do
1823 else
1824 index_cv = index_cv + 1
1825 cv_bias(index_cv) = bias(iSensor)%chans(iChannel)%stddev(iPredictor) * bias(iSensor)%chans(iChannel)%coeffIncr(iPredictor)
1826 end if
1827 end do
1828 end if
1829 end do
1830 end if
1831 end do
1832 end if
1833
1834 end subroutine bcs_cvToCoeff_ad
1835
1836 !-----------------------------------------
1837 ! bcs_writeBias
1838 !-----------------------------------------
1839 subroutine bcs_writeBias(cv_in_opt)
1840 !
1841 !:Purpose: to write bias increments and coefficients (varBC)
1842 !
1843 implicit none
1844
1845 ! Arguments:
1846 real(8), optional, intent(in) :: cv_in_opt(:)
1847
1848 ! Locals:
1849 integer :: iSensor, iChannel, iPredictor, iScan
1850 integer :: jSensor, iChannel2
1851 integer :: nulfile_inc, nulfile_fov, ierr
1852 real(8), pointer :: cv_bias(:)
1853 real(8), target :: dummy4Pointer(1)
1854 character(len=80) :: BgFileName
1855 !for background coeff and write out
1856 integer :: iInstr
1857 integer :: numCoefFile, jCoef, kCoef
1858 character(len=10) :: coefInstrName(tvs_nSensors), temp_instrName
1859 character(len=25) :: filecoeff
1860 logical :: coeffExists
1861 ! these variables are not used but need to be present to satisfy bcs_updateCoeff interface
1862 ! some bcs_updateCoeff arguments could be made optional (todo)
1863 integer :: chans(tvs_nSensors, tvs_maxChannelNumber), nsat, nfov
1864 integer :: nchan(tvs_nSensors)
1865 character(len=10) :: sats(tvs_nsensors) ! satellite names
1866 character(len=7) :: cinstrum ! string: instrument (e.g. AMSUB)
1867
1868 if (.not. biasActive) return
1869
1870 if (present(cv_in_opt)) then
1871 nullify(cv_bias)
1872 if (mmpi_myid == 0) then
1873 if (cvm_subVectorExists('BIAS')) then
1874 cv_bias => cvm_getSubVector(cv_in_opt, 'BIAS')
1875 write(*,*) 'bcs_writeBias: maxval(cv_bias) = ', maxval(cv_bias(:))
1876 else
1877 write(*,*) 'bcs_writeBias: control vector does not include bias coefficients'
1878 return
1879 end if
1880 else
1881 cv_bias => dummy4Pointer
1882 end if
1883 call bcs_cvToCoeff(cv_bias)
1884 end if
1885
1886 ! apply transformation to account for predictor offset
1887
1888 do iSensor = 1, tvs_nSensors
1889 if (bias(iSensor)%numScan > 0) then
1890 do iChannel = 1, bias(iSensor)%numChannels
1891 do iPredictor = 2, bias(iSensor)%chans(iChannel)%numActivePredictors
1892 if (mimicSatbcor) then
1893 bias(iSensor)%chans(iChannel)%coeffIncr(1) = bias(iSensor)%chans(iChannel)%coeffIncr(1) - &
1894 bias(iSensor)%chans(iChannel)%coeffIncr(iPredictor) * bias(iSensor)%chans(iChannel)%coeff_offset(iPredictor)
1895 else
1896 do iScan = 1, bias(iSensor)%numScan
1897 bias(iSensor)%chans(iChannel)%coeffIncr_fov(iScan) = bias(iSensor)%chans(iChannel)%coeffIncr_fov(iScan) - &
1898 bias(iSensor)%chans(iChannel)%coeffIncr(iPredictor) * bias(iSensor)%chans(iChannel)%coeff_offset(iPredictor)
1899 end do
1900 end if
1901 end do
1902 if (bias(iSensor)%chans(iChannel)%numActivePredictors > 0 .and. mmpi_myId ==0 .and. (.not. doRegression)) &
1903 write(*,*) "bcs_writeBias: bias(iSensor)%chans(iChannel)%coeffIncr(:) = ", bias(iSensor)%chans(iChannel)%coeffIncr(:)
1904 end do
1905 end if
1906 end do
1907
1908 if (doRegression) then
1909 call bcs_writeCoeff()
1910 return
1911 end if
1912
1913 if (mmpi_myId == 0) then
1914
1915 ! write out bias coefficient increments in ascii file
1916 nulfile_inc = 0
1917 ierr = fnom(nulfile_inc, './satbias_increment.dat', 'FTN+FMT', 0)
1918
1919 do iSensor = 1, tvs_nSensors
1920 write(nulfile_inc,'(/,1X,"Sensor Index=",I3,", Satellite Name=",A15,", Instrument Name=",A15)') &
1921 iSensor, tvs_satelliteName(iSensor), tvs_instrumentName(iSensor)
1922 if (bias(iSensor)%numScan > 0) then
1923 do iChannel = 1, bias(iSensor)%numChannels
1924 if (bias(iSensor)%chans(iChannel)%isDynamic) then
1925 iChannel2 = bias(iSensor)%chans(iChannel)%channelNum
1926 if (sum(bias(iSensor)%chans(iChannel)%coeffIncr(:)) /= 0.0d0) &
1927 write(nulfile_inc,'(3X,"Channel number=",I4)') iChannel2
1928 do iPredictor = 2, bias(iSensor)%chans(iChannel)%numActivePredictors
1929 if (bias(iSensor)%chans(iChannel)%coeffIncr(iPredictor) /= 0.0d0) &
1930 write(nulfile_inc,'(5X,"Predictor number=",I4,", Coefficient=",e12.4)') &
1931 iPredictor, bias(iSensor)%chans(iChannel)%coeffIncr(iPredictor)
1932 end do
1933 end if
1934 end do
1935 end if
1936 end do
1937
1938 ierr = fclos(nulfile_inc)
1939
1940 ! write out fovbias coefficient increments in ascii file
1941 nulfile_fov = 0
1942 ierr = fnom(nulfile_fov, './fovbias_incre.dat', 'FTN+FMT', 0)
1943 do iSensor = 1, tvs_nSensors
1944 write(nulfile_fov,'(/,1X,"Sensor Index=",I3,", Satellite Name=",A15,", Instrument Name=",A15)') &
1945 iSensor, tvs_satelliteName(iSensor), tvs_instrumentName(iSensor)
1946 if (bias(iSensor)%numScan > 0) then
1947 do iChannel = 1, bias(iSensor)%numChannels
1948 if (bias(iSensor)%chans(iChannel)%isDynamic) then
1949 iChannel2 = bias(iSensor)%chans(iChannel)%channelNum
1950 if (sum(bias(iSensor)%chans(iChannel)%coeffIncr_fov(:)) /= 0.0d0) then
1951 write(nulfile_fov,'(3X,"Channel number=",I4)') iChannel2
1952 write(nulfile_fov,*) bias(iSensor)%chans(iChannel)%coeffIncr_fov(:)
1953 end if
1954 end if
1955 end do
1956 end if
1957 end do
1958 ierr = fclos(nulfile_fov)
1959
1960 end if
1961
1962 ! Find the background coeff_file number and name
1963 do iSensor = 1, tvs_nSensors
1964 numCoefFile = 0
1965 jCoef = 0
1966 do jSensor = 1, tvs_nSensors
1967 temp_instrName = InstrNametoCoeffFileName(tvs_instrumentName(jSensor))
1968 filecoeff = 'coeff_file_' // trim(temp_instrName) // ''
1969 inquire(file=trim(filecoeff), exist=coeffExists)
1970
1971 if (coeffExists) then
1972 numCoefFile = numCoefFile + 1
1973 jCoef = jCoef + 1
1974 coefInstrName(jCoef) = temp_instrName
1975 end if
1976 if (jSensor > 1) then
1977 do kCoef = 1, jCoef - 1
1978 if (temp_instrName == coefInstrName(kCoef)) then
1979 numCoefFile = numCoefFile - 1
1980 jCoef = jCoef - 1
1981 end if
1982 end do
1983 end if
1984 end do
1985 end do
1986
1987 ! update coeff_file_instrument and write out
1988 do iInstr=1, numCoefFile
1989 BgFileName = './coeff_file_' // coefInstrName(iInstr)
1990 call bcs_updateCoeff(tvs_nSensors, NumPredictors, BgFileName, sats, chans, nsat, nchan, nfov, cinstrum)
1991 end do
1992
1993 end subroutine bcs_writeBias
1994
1995 !-----------------------------------------
1996 ! bcs_updateCoeff
1997 !-----------------------------------------
1998 subroutine bcs_updateCoeff(maxsat, maxpred, coeff_file, sats, chans, nsat, nchan, nfov, cinstrum, updateCoeff_opt)
1999 !
2000 !:Purpose: to read, and optionaly update and write out, the coeff files (varBC).
2001 !
2002 implicit none
2003
2004 ! Arguments:
2005 integer, intent(in) :: maxsat, maxpred
2006 character(len=*), intent(in) :: coeff_file
2007 logical, optional, intent(in) :: updateCoeff_opt
2008 integer, intent(out) :: chans(maxsat,tvs_maxChannelNumber) ! channel numbers
2009 integer, intent(out) :: nsat
2010 integer, intent(out) :: nfov
2011 integer, intent(out) :: nchan(maxsat) ! number of channels
2012 character(len=10), intent(out) :: sats(maxsat) ! Satellite names
2013 character(len=*), intent(out) :: cinstrum ! instrument (e.g. AMSUB)
2014
2015 ! Locals:
2016 real(8) :: fovbias(maxsat,tvs_maxChannelNumber,maxfov)
2017 real(8) :: coeff(maxsat,tvs_maxChannelNumber,maxpred)
2018 character(len=2) :: ptypes(maxsat,tvs_maxChannelNumber,maxpred)
2019 integer :: npred(maxsat,tvs_maxChannelNumber) ! number of predictors
2020 integer :: ndata(maxsat,tvs_maxChannelNumber)
2021 ! LOCAL for reading background coeff file
2022 integer :: iSat, jChan, kPred, kFov
2023 logical :: verbose
2024 ! update coeff files
2025 real :: fovbias_an(maxsat,tvs_maxChannelNumber,maxfov)
2026 real :: coeff_an(maxsat,tvs_maxChannelNumber,maxpred)
2027 integer :: iSensor, jChannel, iFov, iPred, totPred
2028 character(len=10) :: tmp_SatName, tmp_InstName
2029 ! write out files
2030 integer :: iuncoef2, ierr, numPred
2031 character(len=80):: filename2
2032 logical :: updateCoeff_opt2
2033 ! sats(nsat) = satellite names
2034 ! chans(nsat, nchan(i)) = channel numbers of each channel of each satellite i
2035 ! npred(nsat, nchan(i)) = number of predictors for each channel of each satellite i
2036 ! fovbias(i, j, k) = bias for satellite i, channel j, FOV k k=1,nfov
2037 ! if FOV not considered for instrument, nfov = 1 and fovbias is global bias for channel
2038 ! coeff(i, j, 1) = regression constant
2039 ! coeff(i, j, 2), ..., coeff(i, j, npred(i, j)) = predictor coefficients
2040 ! nsat, nchan, nfov, cinstrum (output) are determined from file
2041 ! if returned nsat = 0, coeff_file was empty
2042 ! maxpred (input) is max number of predictors
2043 ! maxsat (input) is max number of satellites
2044
2045 ! There are three parts in this subroutine, read, update and write out the coeff files
2046 !
2047 !- 1. read in the background coeff files, this program is read_coeff from genbiascorr
2048 !
2049 if (present(updateCoeff_opt)) then
2050 updateCoeff_opt2 = updateCoeff_opt
2051 else
2052 updateCoeff_opt2 = .true.
2053 end if
2054
2055 verbose = .false.
2056
2057 call read_coeff(sats, chans, fovbias, coeff, nsat, nchan, nfov, &
2058 npred, cinstrum, coeff_file, ptypes, ndata)
2059
2060 ! Transfer of coefficients read from file to the bias structure
2061 satloop :do iSat = 1, nsat !backgroud sat
2062 instloop:do iSensor = 1, tvs_nSensors
2063 ! for Satellite Name
2064 tmp_SatName = SatNameinCoeffFile(tvs_satelliteName(iSensor))
2065 ! for Instrument Name
2066 tmp_InstName = InstrNameinCoeffFile(tvs_instrumentName(iSensor))
2067
2068 if (trim(tmp_SatName) /= trim(sats(iSat)) .or. trim(tmp_InstName) /= trim(cinstrum)) cycle instloop
2069 write(*,*) "bcs_updateCoeff: " // tmp_SatName // " " // tmp_InstName
2070
2071 if (.not. allocated(bias(iSensor)%chans)) cycle instloop
2072
2073 chan1loop:do jChan = 1, nchan(iSat)
2074 chan2loop:do jChannel = 1, bias(iSensor)%numChannels
2075
2076 if (chans(iSat, jChan) /= bias(iSensor)%chans(jChannel)%channelNum) cycle chan2loop
2077
2078 ! part 1 for coeffIncr
2079 do iFov = 1, nfov
2080 bias(iSensor)%chans(jchannel)%coeff_fov(iFov) = fovbias(iSat, jChan, iFov)
2081 end do ! iFov
2082
2083 ! part 2 for coeffIncr_fov
2084 totPred = bias(iSensor)%chans(jchannel)%NumActivePredictors
2085 do iPred = 1, totPred
2086 bias(iSensor)%chans(jchannel)%coeff(iPred) = coeff(iSat, jChan, iPred)
2087 end do ! iPred
2088
2089 end do chan2loop ! jChannel
2090 end do chan1loop !jChan
2091
2092 end do instloop
2093 end do satloop
2094
2095 if (.not. updateCoeff_opt2) return
2096
2097 !
2098 !- 2.update coeff and fovbias
2099 !
2100 coeff_an(:,:,:) = coeff(:,:,:)
2101 fovbias_an(:,:,:) = fovbias(:,:,:)
2102
2103 do iSat = 1, nsat !backgroud sat
2104 do iSensor = 1, tvs_nSensors
2105 ! for Satellite Name
2106 tmp_SatName = SatNameinCoeffFile(tvs_satelliteName(iSensor))
2107 ! for Instrument Name
2108 tmp_InstName = InstrNameinCoeffFile(tvs_instrumentName(iSensor))
2109 if (trim(tmp_SatName) /= trim(sats(iSat)) .or. trim(tmp_InstName) /= trim(cinstrum)) cycle
2110 do jChan = 1, nchan(iSat)
2111 do jChannel = 1, bias(iSensor)%numChannels
2112
2113 if (chans(iSat, jChan) /= bias(iSensor)%chans(jChannel)%channelNum) cycle
2114
2115 ! part 1 for coeffIncr
2116 do iFov = 1, nfov
2117 fovbias_an(iSat, jChan, iFov) = fovbias(iSat, jChan, iFov) + bias(iSensor)%chans(jchannel)%coeffIncr_fov(iFov)
2118 end do ! iFov
2119
2120 ! part 2 for coeffIncr_fov
2121 totPred = bias(iSensor)%chans(jchannel)%NumActivePredictors
2122 do iPred = 1, totPred
2123 coeff_an(iSat, jChan, iPred) = coeff(iSat, jChan, iPred) + bias(iSensor)%chans(jchannel)%coeffIncr(iPred)
2124 end do ! iPred
2125
2126 end do ! jChannel
2127 end do !jChan
2128 end do !iSensor
2129 end do ! iSat
2130
2131 !
2132 !- 3. Write out updated_coeff
2133 !
2134
2135 if (mmpi_myId == 0) then
2136
2137 iuncoef2 = 0
2138 filename2 = './anlcoeffs_' // cinstrum
2139 ierr = fnom(iuncoef2, filename2, 'FTN+FMT', 0)
2140
2141 write(*,*) 'bcs_updateCoeff: write in bcs_updateCoeff'
2142
2143 do iSat = 1, nsat
2144 do jChan = 1, nchan(iSat)
2145 numPred = npred(iSat, jChan)
2146 if (sum(abs(coeff_an(iSat, jchan, 1:numpred+1))) /= 0.d0 .and. sum(abs(fovbias_an(iSat, jChan, 1:nfov))) /= 0.d0) then
2147 write(iuncoef2,'(A52,A8,1X,A7,1X,I6,1X,I8,1X,I2,1X,I3)') &
2148 'SATELLITE, INSTRUMENT, CHANNEL, NOBS, NPRED, NSCAN: ', sats(iSat), cinstrum, chans(iSat,jChan), ndata(isat,jchan), numPred, nfov
2149 write(iuncoef2,'(A7,6(1X,A2))') 'PTYPES:', (ptypes(iSat, jChan, kPred), kPred = 1, numPred)
2150 write(iuncoef2,'(120(1x,ES17.10))') (fovbias_an(iSat, jChan, kFov), kFov = 1, nfov)
2151 write(iuncoef2,*) (coeff_an(iSat, jChan, kPred), kPred = 1, numPred + 1)
2152 end if
2153 end do
2154 end do
2155
2156 ierr = fclos(iuncoef2)
2157
2158 write(*,*) 'bcs_updateCoeff: finish writing coeffient file', filename2
2159
2160 end if
2161
2162 end subroutine bcs_updateCoeff
2163
2164 !-----------------------------------------
2165 ! bcs_writeCoeff
2166 !-----------------------------------------
2167 subroutine bcs_writeCoeff()
2168 !
2169 !:Purpose: write out the coeff files (regression case).
2170 !
2171 implicit none
2172
2173 ! Locals:
2174 integer :: iuncoef, numPred, ierr
2175 character(len=80) :: filename
2176 character(len=80) :: instrName, satNamecoeff
2177 integer :: sensorIndex, nchans, nscan, nfov, kpred, kFov, jChan
2178
2179 if (mmpi_myId == 0) then
2180
2181 SENSORS:do sensorIndex = 1, tvs_nsensors
2182
2183 if (.not. tvs_isReallyPresentMpiGLobal(sensorIndex)) cycle SENSORS
2184
2185 write(*,*) "bcs_writeCoeff: sensorIndex ", sensorIndex
2186
2187 nchans = bias(sensorIndex)%numChannels
2188 nscan = bias(sensorIndex)%numscan
2189
2190 instrName = InstrNameinCoeffFile(tvs_instrumentName(sensorIndex))
2191 satNamecoeff = SatNameinCoeffFile(tvs_satelliteName(sensorIndex))
2192
2193 filename = './anlcoeffs_' // trim(instrName)
2194 call utl_open_asciifile(filename, iuncoef)
2195 nfov = bias(sensorIndex)%numScan
2196 do jChan = 1, nchans
2197 if (bias(sensorIndex)%chans(jChan)%coeff_nobs > 0) then
2198 numPred = bias(sensorIndex)%chans(jChan)%numActivePredictors
2199
2200 write(iuncoef,'(A52,A8,1X,A7,1X,I6,1X,I8,1X,I2,1X,I3)') 'SATELLITE, INSTRUMENT, CHANNEL, NOBS, NPRED, NSCAN: ', &
2201 satNameCoeff, instrName, bias(sensorIndex)%chans(jChan)%channelNum, bias(sensorIndex)%chans(jChan)%coeff_nobs, numPred - 1, nfov
2202 write(iuncoef,'(A7,6(1X,A2))') 'PTYPES:', (predtab(bias(sensorIndex)%chans(jChan)%predictorIndex(kPred)), kPred = 2, numPred)
2203 write(iuncoef,'(120(1x,ES17.10))') (bias(sensorIndex)%chans(jChan)%coeff_fov(kFov), kFov = 1, nfov)
2204 write(iuncoef,'(12(1x,ES17.10))') (bias(sensorIndex)%chans(jChan)%coeff(kPred), kPred = 1, numPred)
2205 end if
2206 end do
2207
2208 ierr = fclos(iuncoef)
2209
2210 if (outCoeffCov) then
2211 filename = './anlcoeffsCov_' // trim(instrName)
2212 call utl_open_asciifile(filename, iuncoef)
2213 do jChan = 1, nchans
2214 if (bias(sensorIndex)%chans(jChan)%coeff_nobs > 0) then
2215 numPred = bias(sensorIndex)%chans(jChan)%numActivePredictors
2216
2217 write(iuncoef,'(A38,A8,1X,A7,1X,I6,1X,I2)') 'SATELLITE, INSTRUMENT, CHANNEL, NPRED: ', &
2218 satNameCoeff, instrName, bias(sensorIndex)%chans(jChan)%channelNum, numPred
2219 do kpred =1, numPred
2220 write(iuncoef,'(10e14.6)') bias(sensorIndex)%chans(jChan)%coeffCov(kpred, :)
2221 end do
2222 end if
2223 end do
2224 ierr = fclos(iuncoef)
2225 end if
2226
2227 end do SENSORS
2228
2229 end if
2230
2231 end subroutine bcs_writeCoeff
2232
2233 !-----------------------------------------
2234 ! bcs_removeBiasCorrection
2235 !-----------------------------------------
2236 subroutine bcs_removeBiasCorrection(obsSpaceData, family_opt)
2237 !
2238 !:Purpose: to remove bias correction from OBS_VAR
2239 ! After the call OBS_VAR contains the uncorrected
2240 ! observation and OBS_BCOR is set to zero
2241 !
2242 implicit none
2243
2244 ! Arguments:
2245 type(struct_obs), intent(inout) :: obsSpaceData
2246 character(len=2), optional, intent(in) :: family_opt
2247
2248 ! Locals:
2249 integer :: nbcor
2250 integer :: bodyIndex
2251 real(8) :: biascor, Obs
2252
2253 if (.not. removeBiasCorrection) return
2254
2255 if (mmpi_myid == 0) write(*,*) 'bcs_removeBiasCorrection: start'
2256
2257 nbcor = 0
2258 call obs_set_current_body_list(obsSpaceData, family_opt)
2259
2260 BODY: do
2261 bodyIndex = obs_getBodyIndex(obsSpaceData)
2262 if (bodyIndex < 0) exit BODY
2263 if (obs_bodyElem_i(obsSpaceData, OBS_ASS, bodyIndex) /= obs_assimilated) cycle BODY
2264 biasCor = obs_bodyElem_r(obsSpaceData, OBS_BCOR, bodyIndex)
2265 Obs = obs_bodyElem_r(obsSpaceData, OBS_VAR, bodyIndex)
2266 if (biasCor /= MPC_missingValue_R8 .and. Obs /= MPC_missingValue_R8) then
2267 call obs_bodySet_r(obsSpaceData, OBS_VAR, bodyIndex, real(Obs - biasCor, pre_obsReal))
2268 call obs_bodySet_r(obsSpaceData, OBS_BCOR, bodyIndex, real(0.d0, pre_obsReal))
2269 nbcor = nbcor + 1
2270 end if
2271 end do BODY
2272
2273 if (mmpi_myid == 0) then
2274 write(*,*) 'bcs_removeBiasCorrection: removed bias correction for ', nbcor, ' observations'
2275 write(*,*) 'bcs_removeBiasCorrection exiting'
2276 end if
2277
2278 end subroutine bcs_removeBiasCorrection
2279
2280 !-----------------------------------------
2281 ! bcs_filterObs
2282 !-----------------------------------------
2283 subroutine bcs_filterObs(obsSpaceData)
2284 !
2285 !:Purpose: to filter radiance observations to include into
2286 ! bias correction offline computation
2287 ! (same rules as in bgck.gen_table)
2288 !
2289 implicit none
2290
2291 ! Arguments:
2292 type(struct_obs), intent(inout) :: obsSpaceData
2293
2294 ! Locals:
2295 integer :: bodyIndex, headerIndex, instrumentIndex, sensorIndex
2296 integer, allocatable :: instrumentList(:)
2297 integer :: assim, flag, codtyp, channelNumber
2298 integer :: isatBufr, instBufr, iplatform, isat, inst, idsat, chanIndx
2299 logical :: lHyperIr, lGeo, lSsmis, lTovs
2300 logical :: condition, condition1, condition2, channelIsAllsky
2301 logical :: channelIsPassive
2302 character(len=10) :: instrName
2303
2304 if (.not. filterObs) return
2305
2306 if (mmpi_myid == 0) write(*,*) 'bcs_filterObs: start'
2307
2308 allocate(instrumentList(tvs_nsensors))
2309 instrumentList(:) = -1
2310 do sensorIndex = 1, tvs_nsensors
2311 instrName = InstrNameinCoeffFile(tvs_instrumentName(sensorIndex))
2312 do instrumentIndex = 1, maxNumInst
2313 if (trim(instrName) == trim(cinst(instrumentIndex))) then
2314 instrumentList(sensorIndex) = instrumentIndex
2315 end if
2316 end do
2317 if (instrumentList(sensorIndex) == -1) then
2318 call utl_abort('bcs_filterObs: Instrument ' // trim(instrName) // &
2319 'missing in CINST table fron NAMBIASSAT namelist section')
2320 end if
2321 end do
2322
2323 call obs_set_current_header_list(obsSpaceData, 'TO')
2324 HEADER: do
2325 headerIndex = obs_getHeaderIndex(obsSpaceData)
2326 if (headerIndex < 0) exit HEADER
2327
2328 codtyp = obs_headElem_i(obsSpaceData, OBS_ITY, headerIndex)
2329
2330 lHyperIr = .false.
2331 lGeo = .false.
2332 lSsmis = .false.
2333 lTovs = .false.
2334
2335 select case (codtyp_get_name(codtyp))
2336 case("ssmis")
2337 lSsmis = .true.
2338 case("csr")
2339 lGeo = .true.
2340 case("airs","iasi","cris","crisfsr")
2341 lHyperIr = .true.
2342 case default
2343 lTovs = .true.
2344 end select
2345
2346 isatBufr = obs_headElem_i(obsSpaceData, OBS_SAT, headerIndex) !BUFR element 1007
2347 instBufr = obs_headElem_i(obsSpaceData, OBS_INS, headerIndex) !BUFR element 2019
2348
2349 call tvs_mapSat(isatBufr, iplatform, isat)
2350 call tvs_mapInstrum(instBufr, inst)
2351
2352 idsat = -1
2353 do sensorIndex = 1, tvs_nsensors
2354 if (tvs_platforms(sensorIndex) == iplatform .and. &
2355 tvs_satellites(sensorIndex) == isat .and. &
2356 tvs_instruments(sensorIndex) == inst) then
2357 idsat = sensorIndex
2358 exit
2359 end if
2360 end do
2361 if (idsat == -1) cycle HEADER
2362
2363 call obs_set_current_body_list(obsSpaceData, headerIndex)
2364
2365 BODY: do
2366 bodyIndex = obs_getBodyIndex(obsSpaceData)
2367 if (bodyIndex < 0) exit BODY
2368
2369 ! determine if instrument/channel function in all-sky mode
2370 channelNumber = nint(obs_bodyElem_r(obsSpaceData, OBS_PPP, bodyIndex))
2371 channelIsAllsky = .false.
2372 if ((tvs_isInstrumUsingCLW(tvs_instruments(idsat)) .or. &
2373 tvs_isInstrumUsingHydrometeors(tvs_instruments(idsat))) .and. &
2374 oer_useStateDepSigmaObs(channelNumber, idsat)) then
2375 channelIsAllsky = .true.
2376 end if
2377 channelIsPassive = .false.
2378 if (passiveChannelNumber(instrumentList(idsat)) > 0) then
2379 channelIsPassive = (utl_findloc(passiveChannelList(instrumentList(idsat),1:passiveChannelNumber(instrumentList(idsat))), channelNumber) > 0)
2380 end if
2381 assim = obs_bodyElem_i(obsSpaceData, OBS_ASS, bodyIndex)
2382 if (assim == obs_notAssimilated) then
2383 call bcs_getChannelIndex(obsSpaceData, idsat, chanIndx, bodyIndex)
2384 if (chanIndx > 0) then
2385 flag = obs_bodyElem_i(obsSpaceData, OBS_FLG, bodyIndex)
2386
2387 if (lSsmis) then ! SSM/I and SSMIS
2388 ! Here "good" means data that are not rejected by first QC program satqc_ssmi(s) (bit 7 ON).
2389 ! Bit 11 is ON for data that are unselected by UTIL or for uncorrected data (or both).
2390 ! Data rejected by first QC program satqc_ssmi(s) have bit 7 switched ON only (in addition to bit 9) as
2391 ! rogue/topo checks are skipped. So if bit 16 (rogue) is ON, bit 7 must be off.
2392 if (.not. offlineMode .and. .not. channelIsPassive) then
2393 if (allModeSsmis) then
2394 ! FLAG test: all good data (corrected/selected or not) that have passed all QC (bit 9 OFF)
2395 condition1 = .not. btest(flag, 9) !' AND (FLAG & 512 = 0)'
2396 ! FLAG test: uncorrected good data that failed rogue check only ([bit 9 ON] + bit 6 OFF + bit 16 ON + bit 18 OFF + [bit 7 OFF])
2397 condition2 = .not. btest(flag, 6) .and. btest(flag, 16) .and. .not. btest(flag, 18) !' AND (FLAG & 64 = 0) AND (FLAG & 65536 = 65536) AND (FLAG & 262144 = 0)'
2398 condition = condition1 .or. condition2
2399 else
2400 ! FLAG test: corrected/selected good data that have passed QC (bits 9,11 OFF) --> data to be assimilated
2401 condition = .not. btest(flag, 9) .and. .not. btest(flag, 11) !' AND (FLAG & 512 = 0) AND (FLAG & 2048 = 0)'
2402 end if
2403 else
2404 ! OFFLINE MODE --> want all observations except data rejected for any reason other than rogue innovation check
2405 condition1 = .not. btest(flag, 9) !' AND (FLAG & 512 = 0)'
2406 ! all good data that passed all QC
2407 ! "good" data that failed rogue check [bit 9 ON, bit 7 OFF, bit 18 OFF]
2408 condition2 = btest(flag, 9) .and. .not. btest(flag, 7) .and. .not. btest(flag, 18) !' AND (FLAG & 512 = 512) AND (FLAG & 128 = 0) AND (FLAG & 262144 = 0)'
2409 condition = condition1 .or. condition2
2410 end if
2411 else if(lTovs) then
2412 ! AMSU-A, AMSU-B/MHS, ATMS, MWHS-2
2413 ! In AMSU case, bit 11 is set for data that are not bias corrected or for unselected channels.
2414 ! BUT unlike other instruments, all AMSU data are bias corrected, whether selected or not
2415 ! so bit 11 = unselected channel (like bit 8 for AIRS/IASI)
2416 ! Bit 9 is set for all other rejections including rogue (9+16) and topography (9+18).
2417 ! In addition, bit 7 is set for channels with bad data or data that should not be assimilated.
2418 if (.not. offlineMode .and. .not. channelIsPassive) then
2419 if (allModeTovs) then
2420 ! FLAG test: all data (selected or not) that have passed QC (bit 9 OFF)
2421 condition1 = .not. btest(flag, 9) !' AND (FLAG & 512 = 0)'
2422 ! FLAG test: uncorrected (bit 6 OFF) data that failed rogue check only (bit (9)/16 ON, 18,7 OFF)
2423 ! NOTE: As all AMSU data are normally bias corrected, query2 will return nothing
2424 condition2 = btest(flag, 16) .and. .not. btest(flag, 6) .and. .not. btest(flag, 18) .and. .not. btest(flag, 7)!' AND (FLAG & 64 = 0) AND (FLAG & 65536 = 65536) AND (FLAG & 262144 = 0) AND (FLAG & 128 = 0)'
2425 condition = condition1 .or. condition2
2426 else
2427 ! FLAG test: selected data (bit 11 OFF) that have passed QC (bit 9 OFF)
2428 condition = .not. btest(flag, 9) .and. .not. btest(flag, 11) !' AND (FLAG & 512 = 0) AND (FLAG & 2048 = 0)'
2429 end if
2430 else ! OFFLINE MODE --> want all observations except data rejected for any reason other than rogue check
2431 condition1 = .not. btest(flag, 9) !' AND (FLAG & 512 = 0)'
2432 ! all good data that passed all QC
2433 ! "good" data that failed rogue check [bit 9 ON, bit 7 OFF, bit 18 OFF]
2434 condition2 = btest(flag, 9) .and. .not. btest(flag, 7) .and. .not. btest(flag, 18) !' AND (FLAG & 512 = 512) AND (FLAG & 128 = 0) AND (FLAG & 262144 = 0)'
2435 condition = condition1 .or. condition2
2436 end if
2437 if (channelIsAllsky) condition = condition .and. .not. btest(flag, 23)
2438 else if(lGeo) then ! CSR case
2439 ! No flag check = all data that have passed QC/filtering
2440 ! (FLAG & 2048 = 0) = bit 11 OFF --> corrected/selected data that have passed QC/filtering
2441 if (allModeCsr .or. offlineMode .or. channelIsPassive) then
2442 condition = .true.
2443 else
2444 condition = .not. btest(flag, 18) ! ' AND (FLAG & 2048 = 0)'
2445 endif
2446 else if (lHyperIr) then ! AIRS, IASI and CRIS
2447 ! (FLAG & 2560 = 0) = bits 9, 11 OFF --> data that passed QC (rogue and other)
2448 ! (FLAG & 11010176 = 0) = bits 7,19,21,23 OFF --> "good" data (corrected/selected or not)! (FLAG & 64 = 64) = bit 6 ON --> bias corrected data
2449 ! (FLAG & 256 = 0) = bit 8 OFF --> passed selction (not blacklisted, UTIL=1)
2450 ! (FLAG & 2048 = 2048) = bit 11 ON
2451 ! (FLAG & 65536 = 65536) = bit 16 ON --> rogue check failure
2452 ! (FLAG & 524288 = 0) = bit 19 OFF --> not surface affected [experimental, bit 11 may not be on if data
2453 ! are to be assimilated]
2454 ! (FLAG & 2097152 = 0) = bit 21 OFF --> not rejected due to model top transmittance
2455 ! (FLAG & 8388608 = 0) = bit 23 OFF --> "clear sky" radiance [experimental, bit 11 may not be on if cloudy
2456 ! data are assimilated]! (FLAG & 128 = 0) = bit 7 OFF --> not shortwave channel during day
2457 ! (FLAG & 512 = 0) = bit 9 OFF --> non-erroneous data that passed O-P rogue check
2458 !! AIRS, IASI:! bit 8 ON: blacklisted/unselected channel (UTIL=0)
2459 ! bit 9 ON: erroneous/suspect data (9), data failed O-P check (9+16)
2460 ! bit 11 ON: cloud (11+23), surface (11+19), model top transmittance (11+21), shortwave channel+daytime (11+7)
2461 ! not bias corrected (11) (with bit 6 OFF)
2462 if (.not. offlineMode .and. .not. channelIsPassive) then
2463 if (allModeHyperIr) then
2464 ! good data that have passed all QC (bits 9 and 7,19,21,23 OFF), corrected/selected or not
2465 condition1 = .not. btest(flag, 9) .and. .not. btest(flag, 7) .and. .not. btest(flag, 19) .and. .not. btest(flag, 21) .and. .not. btest(flag, 23) !' AND (FLAG & 512 = 0) AND (FLAG & 11010176 = 0)'
2466 ! uncorrected (6 OFF, [11 ON]) good data (7,19,21,23 OFF) that failed QC rogue check only (bits [9],16 ON), selected or not
2467 condition2 = .not. btest(flag, 6) .and. btest(flag, 11) .and. .not. btest(flag, 17) .and. .not. btest(flag, 19) .and. .not. btest(flag, 21) .and. .not. btest(flag, 23)
2468 !' AND (FLAG & 64 = 0) AND (FLAG & 65536 = 65536) AND (FLAG & 11010176 = 0)'
2469 condition = condition1 .or. condition2
2470 else
2471 ! corrected data that passed all QC and selection excluding cloud/sfc affected obs
2472 condition = .not. btest(flag, 9) .and. .not. btest(flag, 11) .and. .not. btest(flag, 8) .and. .not. btest(flag, 23) .and. .not. btest(flag, 19)
2473 !' AND (FLAG & 2560 = 0) AND (FLAG & 256 = 0) AND (FLAG & 8388608 = 0) AND (FLAG & 524288 = 0)'
2474 end if
2475 else! OFFLINE MODE --> Want all observations except data rejected for any reason other than innovation rogue check
2476 ! Assumes that type S or N correction has been applied to all data/channels (all data "corrected")
2477 ! data that passed all QC
2478 condition1 = .not. btest(flag, 9) .and. .not. btest(flag, 7) .and. .not. btest(flag, 19) .and. .not. btest(flag, 21) .and. .not. btest(flag, 23)
2479 !' AND (FLAG & 512 = 0) AND (FLAG & 11010176 = 0)'
2480 ! good data (7,19,21,23 OFF) that failed QC rogue check only (bits [9],16 ON)
2481 condition2 = btest(flag, 9) .and. btest(flag, 16) .and. .not. btest(flag, 7) .and. .not. btest(flag, 19) .and. .not. btest(flag, 21) .and. .not. btest(flag, 23) !' AND (FLAG & 65536 = 65536) AND (FLAG & 11010176 = 0)'
2482 condition = condition1 .or. condition2
2483 end if
2484 end if
2485
2486 if (condition) assim = obs_Assimilated
2487
2488 call obs_bodySet_i(obsSpaceData, OBS_ASS, bodyIndex, assim)
2489 end if
2490 end if
2491 end do BODY
2492 end do HEADER
2493
2494 deallocate(instrumentList)
2495
2496 end subroutine bcs_filterObs
2497
2498 !-----------------------------------------
2499 ! bcs_applyBiasCorrection
2500 !-----------------------------------------
2501 subroutine bcs_applyBiasCorrection(obsSpaceData, column, family_opt)
2502 !
2503 !:Purpose: to apply bias correction from OBS_BCOR to
2504 ! obsSpaceData column column.
2505 ! After the call obsSpaceData body column contains the corrected
2506 ! observation or O-F and OBS_BCOR is not modified.
2507 implicit none
2508
2509 ! Arguments:
2510 type(struct_obs), intent(inout) :: obsSpaceData
2511 integer, intent(in) :: column !obsSpaceData column
2512 character(len=2), optional, intent(in) :: family_opt
2513
2514 ! Locals:
2515 integer :: nbcor
2516 integer :: bodyIndex
2517 real(8) :: biascor, Obs
2518 integer :: flag
2519
2520 if (.not. biasActive) return
2521 if (trim(biasMode) /= "apply") return
2522
2523 if (mmpi_myid == 0) write(*,*) 'bcs_applyBiasCorrection: start'
2524
2525 nbcor = 0
2526 call obs_set_current_body_list(obsSpaceData, family_opt)
2527
2528 BODY: do
2529 bodyIndex = obs_getBodyIndex(obsSpaceData)
2530 if (bodyIndex < 0) exit BODY
2531 if (obs_bodyElem_i(obsSpaceData, OBS_ASS, bodyIndex) /= obs_assimilated) cycle BODY
2532 biasCor = obs_bodyElem_r(obsSpaceData, OBS_BCOR, bodyIndex)
2533 if (biasCor /= MPC_missingValue_R8) then
2534 Obs = obs_bodyElem_r(obsSpaceData, column, bodyIndex)
2535 if (Obs /= MPC_missingValue_R8) then
2536 call obs_bodySet_r(obsSpaceData, column, bodyIndex, real(Obs + biasCor, pre_obsReal))
2537 flag = obs_bodyElem_i(obsSpaceData, OBS_FLG, bodyIndex)
2538 flag = ibset(flag, 6)
2539 call obs_bodySet_i(obsSpaceData, OBS_FLG, bodyIndex, flag)
2540 nbcor = nbcor + 1
2541 end if
2542 end if
2543 end do BODY
2544
2545 if (mmpi_myid == 0) then
2546 write(*,*) 'bcs_applyBiasCorrection: apply bias correction for ', nbcor, ' observations'
2547 write(*,*) 'bcs_applyBiasCorrection exiting'
2548 end if
2549
2550 end subroutine bcs_applyBiasCorrection
2551
2552 !-----------------------------------------
2553 ! bcs_refreshBiasCorrection
2554 !-----------------------------------------
2555 subroutine bcs_refreshBiasCorrection(obsSpaceData, columnTrlOnTrlLev)
2556 !
2557 !:Purpose: to apply bias correction from read coefficient file to OBS_VAR
2558 ! After the call OBS_VAR contains the corrected observation
2559 ! and OBS_BCOR is set to applied bias correction
2560 !
2561 implicit none
2562
2563 ! Arguments:
2564 type(struct_obs), intent(inout) :: obsSpaceData
2565 type(struct_columnData), intent(inout) :: columnTrlOnTrlLev
2566
2567 if (.not.biasActive) return
2568 if (.not. refreshBiasCorrection) return
2569
2570 if (mmpi_myid == 0) write(*,*) 'bcs_refreshBiasCorrection: start'
2571 call bcs_calcBias(obsSpaceData, columnTrlOnTrlLev)
2572 call bcs_applyBiasCorrection(obsSpaceData, OBS_VAR, "TO")
2573 if (mmpi_myid == 0) write(*,*) 'bcs_refreshBiasCorrection: exit'
2574
2575 end subroutine bcs_refreshBiasCorrection
2576
2577 !-----------------------------------------
2578 ! bcs_getRadiosondeWeight
2579 !-----------------------------------------
2580 subroutine bcs_getRadiosondeWeight(obsSpaceData, lmodify_obserror_opt)
2581 !
2582 !:Purpose: initialize the weights to give more importance to data near radiosonde stations
2583 !
2584 implicit none
2585
2586 ! Arguments:
2587 type(struct_obs), intent(inout) :: obsSpaceData
2588 logical, optional, intent(in) :: lmodify_obserror_opt
2589
2590 ! Locals:
2591 integer :: iobs, headerIndex, idatyp, nobs, bodyIndex, stepIndex
2592 logical :: lmodify_obserror
2593 real(8) :: sigmaObs
2594 type(struct_gsv) :: stateVector_mask, stateVector_mask_4d
2595
2596 lmodify_obserror =.false.
2597
2598 if (present(lmodify_obserror_opt)) lmodify_obserror = lmodify_obserror_opt
2599
2600
2601 if (weightedEstimate) then
2602 call hco_SetupFromFile(hco_mask, './raob_masque.std', 'WEIGHT', GridName_opt='RadiosondeWeight', varName_opt='WT')
2603 call vco_SetupFromFile(vco_mask, './raob_masque.std') ! IN
2604
2605 call gsv_allocate(stateVector_mask, 1, hco_mask, vco_mask, dateStampList_opt=[-1], varNames_opt=["WT"], &
2606 dataKind_opt=4, mpi_local_opt=.true., mpi_distribution_opt="Tiles")
2607 call gsv_allocate(stateVector_mask_4d, tim_nstepobs, hco_mask, vco_mask, dateStampList_opt=[ (-1, stepIndex = 1, tim_nstepobs) ], varNames_opt=["WT"], &
2608 dataKind_opt=4, mpi_local_opt=.true., mpi_distribution_opt="Tiles")
2609
2610 call gio_readFromFile(stateVector_mask, './raob_masque.std', 'WEIGHT', 'O', unitConversion_opt=.false., &
2611 containsFullField_opt=.false.)
2612
2613 do stepIndex = 1, tim_nstepobs
2614 call gsv_copy(stateVector_mask, stateVector_mask_4d, stepIndexOut_opt=stepIndex)
2615 end do
2616 call gsv_deallocate(stateVector_mask)
2617 call col_setVco(column_mask, vco_mask)
2618 nobs = obs_numHeader(obsSpaceData)
2619 call col_allocate(column_mask, nobs, beSilent_opt=.false., varNames_opt=["WT"])
2620
2621 call s2c_nl(stateVector_mask_4d, obsSpaceData, column_mask, hco_mask, &
2622 'NEAREST', varName_opt="WT", moveObsAtPole_opt=.true.)
2623
2624 call obs_set_current_header_list(obsSpaceData, 'TO')
2625 iobs = 0
2626 HEADER: do
2627 headerIndex = obs_getHeaderIndex(obsSpaceData)
2628 if (headerIndex < 0) exit HEADER
2629
2630 ! process only radiance data to be assimilated?
2631 idatyp = obs_headElem_i(obsSpaceData, OBS_ITY, headerIndex)
2632 if (.not. tvs_isIdBurpTovs(idatyp)) then
2633 write(*,*) 'bcs_getRadiosondeWeight: warning unknown radiance codtyp present check NAMTOVSINST', idatyp
2634 cycle HEADER
2635 end if
2636 iobs = iobs + 1
2637
2638 RadiosondeWeight(iobs) = col_getElem(column_mask, 1, headerIndex)
2639
2640 if (lmodify_obserror) then
2641 call obs_set_current_body_list(obsSpaceData, headerIndex)
2642 BODY: do
2643 bodyIndex = obs_getBodyIndex(obsSpaceData)
2644 if (bodyIndex < 0) exit BODY
2645
2646 if (obs_bodyElem_i(obsSpaceData, OBS_ASS, bodyIndex) /= obs_assimilated) cycle BODY
2647
2648 sigmaObs = obs_bodyElem_r(obsSpaceData, OBS_OER, bodyIndex)
2649
2650 sigmaObs = sigmaObs / sqrt(RadiosondeWeight(iobs))
2651 call obs_bodySet_r(obsSpaceData, OBS_OER, bodyIndex, sigmaObs)
2652
2653 end do BODY
2654
2655 end if
2656
2657 end do HEADER
2658 call gsv_deallocate(stateVector_mask_4d)
2659 else
2660 RadiosondeWeight(:) = 1.d0
2661 end if
2662
2663 end subroutine bcs_getRadiosondeWeight
2664
2665 !-----------------------------------------
2666 ! bcs_do_regression
2667 !-----------------------------------------
2668 subroutine bcs_do_regression(columnTrlOnTrlLev, obsSpaceData)
2669 !
2670 !:Purpose: compute the bias correction coefficients by linear regresion
2671 !
2672 implicit none
2673
2674 ! Arguments:
2675 type(struct_obs), intent(inout) :: obsSpaceData
2676 type(struct_columnData), intent(inout) :: columnTrlOnTrlLev
2677
2678 ! Locals:
2679 integer :: iSensor, iChannel, npred, nchans, nscan, ndim, ndimmax
2680 integer :: sensorIndex, iPred1, jPred1, iobs
2681 integer :: headerIndex, idatyp, nPredMax, ierr, iFov, iScan, idim
2682 integer :: indxtovs, bodyIndex, chanIndx, predstart, ntot
2683 real(8) :: OmF, sigmaObs, lambda, norm
2684 real(8) :: predictor(NumPredictors)
2685 real(8), allocatable :: Matrix(:,:,:), Vector(:,:)
2686 real(8), allocatable :: matrixMpiGlobal(:,:,:), vectorMpiGlobal(:,:)
2687 real(8), allocatable :: pIMatrix(:,:), OmFBias(:,:), omfBiasMpiGlobal(:,:)
2688 real(8), allocatable :: BMatrixMinusOne(:,:), LineVec(:,:)
2689 integer, allocatable :: OmFCount(:,:), omfCountMpiGlobal(:,:)
2690 character(len=80) :: errorMessage
2691
2692 write(*,*) "bcs_do_regression: start"
2693 if (.not. allocated(trialHeight300m1000)) then
2694 call bcs_getTrialPredictors(obsSpaceData, columnTrlOnTrlLev)
2695 call bcs_computePredictorBiases(obsSpaceData)
2696 end if
2697
2698 call bcs_getRadiosondeWeight(obsSpaceData)
2699
2700 if (outOmFPredCov) call bcs_outputCvOmPPred(obsSpaceData)
2701
2702 SENSORS:do sensorIndex = 1, tvs_nsensors
2703
2704 if (.not. tvs_isReallyPresentMpiGLobal(sensorIndex)) cycle SENSORS
2705
2706 nchans = bias(sensorIndex)%numChannels
2707 nscan = bias(sensorIndex)%numscan
2708 npredMax = maxval(bias(sensorIndex)%chans(:)%numActivePredictors)
2709
2710 if (mimicSatbcor) then
2711 ndimmax = npredMax
2712 allocate(OmFBias(nchans,nscan))
2713 OmFBias(:,:) = 0.d0
2714 else
2715 ndimmax = npredMax + nscan - 1
2716 end if
2717
2718 allocate(OmFCount(nchans,nscan))
2719 OmFCount(:,:) = 0
2720
2721 allocate(Matrix(nchans,ndimmax,ndimmax))
2722 Matrix(:,:,:) = 0.d0
2723
2724 allocate(Vector(nchans,ndimmax))
2725 Vector(:,:) = 0.d0
2726
2727 allocate(LineVec(1,ndimmax))
2728
2729 allocate(pIMatrix(ndimmax,ndimmax))
2730
2731
2732 ! First pass throught ObsSpaceData to estimate scan biases and count data
2733
2734 call obs_set_current_header_list(obsSpaceData, 'TO')
2735 HEADER1: do
2736 headerIndex = obs_getHeaderIndex(obsSpaceData)
2737 if (headerIndex < 0) exit HEADER1
2738
2739 ! process only radiance data to be assimilated?
2740 idatyp = obs_headElem_i(obsSpaceData, OBS_ITY, headerIndex)
2741 if (.not. tvs_isIdBurpTovs(idatyp)) then
2742 write(*,*) 'bcs_do_regression: warning unknown radiance codtyp present check NAMTOVSINST', idatyp
2743 cycle HEADER1
2744 end if
2745 indxtovs = tvs_tovsIndex(headerIndex)
2746 if (indxtovs < 0) cycle HEADER1
2747
2748 iSensor = tvs_lsensor(indxTovs)
2749 if (iSensor /= sensorIndex) cycle HEADER1
2750 iFov = obs_headElem_i(obsSpaceData, OBS_FOV, headerIndex)
2751 if (nscan > 1) then
2752 iScan = iFov
2753 else
2754 iScan = 1
2755 end if
2756
2757 call obs_set_current_body_list(obsSpaceData, headerIndex)
2758 BODY1: do
2759 bodyIndex = obs_getBodyIndex(obsSpaceData)
2760 if (bodyIndex < 0) exit BODY1
2761
2762 if (obs_bodyElem_i(obsSpaceData, OBS_ASS, bodyIndex) /= obs_assimilated) cycle BODY1
2763 call bcs_getChannelIndex(obsSpaceData, iSensor, chanIndx, bodyIndex)
2764 if (chanindx > 0) then
2765 if (mimicSatBcor) then
2766 OmF = obs_bodyElem_r(obsSpaceData, OBS_OMP, bodyIndex)
2767 OmFBias(chanIndx,iScan) = OmFBias(chanIndx,iScan) + OmF
2768 end if
2769 OmFCount(chanIndx,iScan) = OmFCount(chanIndx,iScan) + 1
2770 end if
2771 end do BODY1
2772 end do HEADER1
2773
2774 if (mimicSatbcor) allocate(omfBiasMpiGlobal(nchans,nscan))
2775
2776 allocate(omfCountMpiGlobal(nchans,nscan))
2777
2778 if (mimicSatbcor) then
2779 call mmpi_reduce_sumR8_2d( OmFBias, omfBiasMpiGlobal, 0, "GRID" )
2780 end if
2781 call rpn_comm_reduce(OmFCount, omfCountMpiGlobal, size(omfCountMpiGlobal), "mpi_integer", "MPI_SUM", 0, "GRID", ierr)
2782
2783 if (ierr /= 0) then
2784 write(errorMessage,*) "bcs_do_regression: MPI communication error 2", ierr
2785 call utl_abort(errorMessage)
2786 end if
2787 if (mimicSatbcor) then
2788 if (mmpi_myId == 0) then
2789 where(omfCountMpiGlobal == 0) omfBiasMpiGlobal = 0.d0
2790 where(omfCountMpiGlobal > 0) omfBiasMpiGlobal = omfBiasMpiGlobal / omfCountMpiGlobal
2791 end if
2792 call rpn_comm_bcast(omfBiasMpiGlobal, size(omfBiasMpiGlobal), "mpi_double_precision", 0, "GRID", ierr)
2793 if (ierr /= 0) then
2794 write(errorMessage,*) "bcs_do_regression: MPI communication error 3", ierr
2795 call utl_abort(errorMessage)
2796 end if
2797 do iChannel = 1, nchans
2798 bias(sensorIndex)%chans(iChannel)%coeff_fov(:) = omfBiasMpiGlobal(iChannel, :)
2799 end do
2800 deallocate(omfBiasMpiGlobal)
2801 end if
2802
2803 ! Second pass to fill matrices and vectors
2804 call obs_set_current_header_list(obsSpaceData, 'TO')
2805 iobs = 0
2806 HEADER2: do
2807 headerIndex = obs_getHeaderIndex(obsSpaceData)
2808 if (headerIndex < 0) exit HEADER2
2809
2810 ! process only radiance data to be assimilated?
2811 idatyp = obs_headElem_i(obsSpaceData, OBS_ITY, headerIndex)
2812 if (.not. tvs_isIdBurpTovs(idatyp)) then
2813 write(*,*) 'bcs_do_regression: warning unknown radiance codtyp present check NAMTOVSINST', idatyp
2814 cycle HEADER2
2815 end if
2816 iobs = iobs + 1
2817 indxtovs = tvs_tovsIndex(headerIndex)
2818 if (indxtovs < 0) cycle HEADER2
2819
2820 iSensor = tvs_lsensor(indxTovs)
2821 if (iSensor /= sensorIndex) cycle HEADER2
2822
2823 call obs_set_current_body_list(obsSpaceData, headerIndex)
2824 iFov = obs_headElem_i(obsSpaceData, OBS_FOV, headerIndex)
2825 if (nscan > 1) then
2826 iScan = iFov
2827 else
2828 iScan = 1
2829 end if
2830 BODY2: do
2831 bodyIndex = obs_getBodyIndex(obsSpaceData)
2832 if (bodyIndex < 0) exit BODY2
2833 if (obs_bodyElem_i(obsSpaceData, OBS_ASS, bodyIndex) /= obs_assimilated) cycle BODY2
2834 call bcs_getChannelIndex(obsSpaceData, iSensor, chanIndx, bodyIndex)
2835 if (chanIndx > 0) then
2836 call bcs_getPredictors(predictor, headerIndex, iobs, chanIndx, obsSpaceData)
2837 OmF = obs_bodyElem_r(obsSpaceData, OBS_OMP, bodyIndex)
2838
2839 if (mimicSatbcor) OmF = OmF - bias(sensorIndex)%chans(chanIndx)%coeff_fov(iScan)
2840
2841 LineVec(:,:) = 0.d0
2842
2843 if (mimicSatbcor) then
2844 idim = 0
2845 predstart = 1
2846 lambda = 1.d0
2847 else
2848 LineVec(1,iScan) = 1.d0
2849 idim = nscan
2850 predstart = 2
2851 sigmaObs = obs_bodyElem_r(obsSpaceData, OBS_OER, bodyIndex)
2852 lambda = 1.d0 / (sigmaObs ** 2)
2853 end if
2854
2855 lambda = lambda * RadiosondeWeight(iobs)
2856
2857 do iPred1 = predstart, bias(iSensor)%chans(chanIndx)%NumActivePredictors
2858 jPred1 = bias(iSensor)%chans(chanIndx)%PredictorIndex(iPred1)
2859 idim = idim + 1
2860 LineVec(1,idim) = predictor(jPred1)
2861 end do
2862 Matrix(chanindx,:,:) = Matrix(chanindx,:,:) + matmul(transpose(LineVec),LineVec) * lambda
2863 Vector(chanIndx,:) = Vector(chanIndx,:) + LineVec(1,:) * OmF * lambda
2864 end if
2865 end do BODY2
2866 end do HEADER2
2867
2868 if (mmpi_myId == 0) then
2869 allocate(matrixMpiGlobal(nchans,ndimmax,ndimmax))
2870 allocate(vectorMpiGlobal(nchans,ndimmax))
2871 else
2872 allocate(matrixMpiGlobal(1,1,1))
2873 allocate(vectorMpiGlobal(1,1))
2874 end if
2875
2876 ! communication MPI pour tout avoir sur tache 0
2877 call mmpi_reduce_sumR8_3d( Matrix, matrixMpiGlobal, 0, "GRID" )
2878 call mmpi_reduce_sumR8_2d( Vector, vectorMpiGlobal, 0, "GRID" )
2879
2880 do iChannel = 1, nchans
2881
2882 if (mmpi_myId == 0) then
2883 ntot = sum(omfCountMpiGlobal(iChannel, :))
2884 bias(sensorIndex)%chans(iChannel)%coeff_nobs = ntot
2885 if (ntot > 0 .and. .not. mimicSatbcor) then
2886 norm = 1.d0 / (ntot)
2887 matrixMpiGlobal(iChannel,:,:) = matrixMpiGlobal(iChannel,:,:) * norm
2888 vectorMpiGlobal(iChannel,:) = vectorMpiGlobal(iChannel,:) * norm
2889 end if
2890
2891 nPred = bias(sensorIndex)%chans(iChannel)%numActivePredictors
2892 if (mimicSatbcor) then
2893 ndim = npred
2894 else
2895 ndim = npred + nscan -1
2896 allocate(BMatrixMinusOne(ndim,ndim))
2897 BMatrixMinusOne(:,:) = 0.d0
2898 BMatrixMinusOne(1:nscan,1:nscan) = matmul(bias(sensorIndex)%BMinusHalfScanBias, bias(sensorIndex)%BMinusHalfScanBias)
2899 do iPred1 = 2, bias(sensorIndex)%chans(iChannel)%numActivePredictors
2900 BMatrixMinusOne(nscan - 1 + iPred1,nscan - 1 + iPred1) = (1.d0 / (bias(sensorIndex)%chans(iChannel)%stddev(iPred1)) ** 2)
2901 end do
2902 matrixMpiGlobal(iChannel,1:ndim,1:ndim) = matrixMpiGlobal(iChannel,1:ndim,1:ndim) + BMatrixMinusOne(:,:)
2903 deallocate(BMatrixMinusOne)
2904 end if
2905
2906 pIMatrix(:,:) = 0.d0
2907 call utl_pseudo_inverse(matrixMpiGlobal(iChannel, 1:ndim, 1:ndim), pIMatrix(1:ndim, 1:ndim))
2908 LineVec(1,1:ndim) = matmul(pIMatrix(1:ndim,1:ndim), vectorMpiGlobal(iChannel,1:ndim))
2909 !call dsymv("L", ndim, 1.d0, pIMatrix, ndim,vectorMpiGlobal(iChannel,:), 1, 0.d0, LineVec(1,1:ndim), 1)
2910 end if
2911
2912 call rpn_comm_bcast(ndim, 1, "mpi_integer", 0, "GRID", ierr)
2913 call rpn_comm_bcast(npred, 1, "mpi_integer", 0, "GRID", ierr)
2914
2915 call rpn_comm_bcast(LineVec(1,1:ndim), ndim, "mpi_double_precision", 0, "GRID", ierr)
2916 if (ierr /= 0) then
2917 write(errorMessage,*) "bcs_do_regression: MPI communication error 6", ierr
2918 call utl_abort(errorMessage)
2919 end if
2920
2921 if (outCoeffCov) then
2922 allocate (bias(sensorIndex)%chans(iChannel)%coeffCov(ndim,ndim))
2923 call rpn_comm_bcast(pIMatrix(1:ndim, 1:ndim), ndim * ndim, "mpi_double_precision", 0, "GRID", ierr)
2924 bias(sensorIndex)%chans(iChannel)%coeffCov(:,:) = pIMatrix(1:ndim,1:ndim)
2925 end if
2926
2927 if (mimicSatbcor) then
2928 bias(sensorIndex)%chans(iChannel)%coeff(:) = LineVec(1,1:npred)
2929 else
2930 bias(sensorIndex)%chans(iChannel)%coeff_fov(:) = LineVec(1,1:nscan)
2931 bias(sensorIndex)%chans(iChannel)%coeff(1) = 0.d0
2932 bias(sensorIndex)%chans(iChannel)%coeff(2:) = LineVec(1,nscan+1:ndim)
2933 end if
2934
2935 end do
2936
2937 deallocate(LineVec)
2938 deallocate(Matrix)
2939 deallocate(Vector )
2940 deallocate(omfCountMpiGlobal)
2941 deallocate(matrixMpiGlobal)
2942 deallocate(vectorMpiGlobal)
2943 deallocate(pIMatrix)
2944 if (allocated(OmFBias)) deallocate(OmFBias)
2945 deallocate(OmFCount)
2946
2947 end do SENSORS
2948
2949 end subroutine bcs_do_regression
2950
2951 !-----------------------------------------
2952 ! bcs_outputCvOmPPred
2953 !-----------------------------------------
2954 subroutine bcs_outputCvOmPPred(obsSpaceData)
2955 !
2956 !:Purpose: compute and output OmF-predictors covariance and correlation matrices
2957 !
2958 implicit none
2959
2960 ! Arguments:
2961 type(struct_obs), intent(inout) :: obsSpaceData
2962
2963 ! Locals:
2964 integer :: sensorIndex, headerIndex, bodyIndex, channelIndex, predictorIndex,predictorIndex2,nchans
2965 integer :: idatyp, indxtovs, iSensor, chanIndx, Iobs, ierr
2966 Real(8):: OmF
2967 real(8), allocatable :: OmFBias(:), Matrix(:,:,:), PredBias(:,:)
2968 integer, allocatable :: Count(:), CountMpiGlobal(:)
2969 real(8), allocatable :: OmFBiasMpiGlobal(:), predBiasMpiGlobal(:,:), MatrixMpiGLobal(:,:,:)
2970 character(len=128) :: errorMessage
2971 real(8) :: vector(1,numPredictors), predictor(numPredictors),correlation(numPredictors,numPredictors)
2972 real(8) :: sigma(numPredictors)
2973 integer :: iuncov, iuncorr
2974
2975 write(*,*) "bcs_outputCvOmPPred: Starting"
2976
2977 SENSORS:do sensorIndex = 1, tvs_nsensors
2978
2979 if (.not. tvs_isReallyPresentMpiGLobal(sensorIndex)) cycle SENSORS
2980 write(*,*) "sensor ", sensorIndex
2981
2982 nchans = bias(sensorIndex)%numChannels ! from bcif
2983
2984 allocate(OmFBias(nchans))
2985 OmFBias(:) = 0.d0
2986
2987 allocate(predBias(nchans,numPredictors))
2988 predBias(:,:) = 0.d0
2989
2990 allocate(Count(nchans))
2991 Count(:) = 0
2992
2993 iobs = 0
2994 ! First pass throught ObsSpaceData to estimate biases and count data
2995 call obs_set_current_header_list(obsSpaceData, 'TO')
2996 HEADER1: do
2997 headerIndex = obs_getHeaderIndex(obsSpaceData)
2998 if (headerIndex < 0) exit HEADER1
2999 ! process only radiance data to be assimilated?
3000 idatyp = obs_headElem_i(obsSpaceData, OBS_ITY, headerIndex)
3001 if (.not. tvs_isIdBurpTovs(idatyp)) then
3002 write(*,*) 'bcs_outputCvOmPPred: warning unknown radiance codtyp present check NAMTOVSINST', idatyp
3003 cycle HEADER1
3004 end if
3005 iobs = iobs + 1
3006 indxtovs = tvs_tovsIndex(headerIndex)
3007 if (indxtovs < 0) cycle HEADER1
3008 iSensor = tvs_lsensor(indxTovs)
3009 if (iSensor /= sensorIndex) cycle HEADER1
3010
3011 call obs_set_current_body_list(obsSpaceData, headerIndex)
3012 BODY1: do
3013 bodyIndex = obs_getBodyIndex(obsSpaceData)
3014 if (bodyIndex < 0) exit BODY1
3015 if (obs_bodyElem_i(obsSpaceData, OBS_ASS, bodyIndex) /= obs_assimilated) cycle BODY1
3016 call bcs_getChannelIndex(obsSpaceData, iSensor, chanIndx, bodyIndex)
3017 if (chanindx > 0) then
3018 OmF = obs_bodyElem_r(obsSpaceData, OBS_OMP, bodyIndex)
3019 OmFBias(chanIndx) = OmFBias(chanIndx) + OmF
3020 count(chanIndx) = count(chanIndx) + 1
3021 call bcs_getPredictors(predictor, headerIndex, iobs, chanIndx, obsSpaceData)
3022 predBias(chanIndx,:) = predBias(chanIndx,:) + predictor(:)
3023 end if
3024 end do BODY1
3025 end do HEADER1
3026
3027 allocate(omfBiasMpiGlobal(nchans))
3028 allocate(countMpiGlobal(nchans))
3029 allocate(predBiasMpiGlobal(nchans,numPredictors))
3030
3031 call mmpi_reduce_sumR8_1d(OmFBias, omfBiasMpiGlobal, 0, "GRID" )
3032 call mmpi_reduce_sumR8_2d(predBias, predBiasMpiGlobal, 0, "GRID" )
3033
3034 call rpn_comm_reduce(count, countMpiGlobal, size(countMpiGlobal), "mpi_integer", "MPI_SUM", 0, "GRID", ierr)
3035 if (ierr /= 0) then
3036 write(errorMessage,*) "bcs_outputCvOmPPred: MPI communication error 1", ierr
3037 call utl_abort(errorMessage)
3038 end if
3039
3040 if (mmpi_myId == 0) then
3041 where(countMpiGlobal == 0) omfBiasMpiGlobal = 0.d0
3042 where(countMpiGlobal > 0) omfBiasMpiGlobal = omfBiasMpiGlobal / countMpiGlobal
3043 do channelIndex =1, nchans
3044 if (countMpiGlobal(channelIndex) == 0) predBiasMpiGlobal(channelIndex,:) = 0.d0
3045 if (countMpiGlobal(channelIndex) > 0) predBiasMpiGlobal(channelIndex,:) = predBiasMpiGlobal(channelIndex,:) / countMpiGlobal(channelIndex)
3046 end do
3047 end if
3048 call rpn_comm_bcast(omfBiasMpiGlobal, size(omfBiasMpiGlobal), "mpi_double_precision", 0, "GRID", ierr)
3049 if (ierr /= 0) then
3050 write(errorMessage,*) "bcs_outputCvOmPPred: MPI communication error 4", ierr
3051 call utl_abort(errorMessage)
3052 end if
3053 call rpn_comm_bcast(predBiasMpiGlobal, size(predBiasMpiGlobal), "mpi_double_precision", 0, "GRID", ierr)
3054 if (ierr /= 0) then
3055 write(errorMessage,*) "bcs_outputCvOmPPred: MPI communication error 3", ierr
3056 call utl_abort(errorMessage)
3057 end if
3058 call rpn_comm_bcast(countMpiGlobal, size(countMpiGlobal), "mpi_integer", 0, "GRID", ierr)
3059 if (ierr /= 0) then
3060 write(errorMessage,*) "bcs_outputCvOmPPred: MPI communication error 4", ierr
3061 call utl_abort(errorMessage)
3062 end if
3063
3064 deallocate(OmFBias)
3065 deallocate(predBias)
3066 deallocate(Count)
3067 allocate(matrix(nchans,numPredictors,numPredictors))
3068 matrix(:,:,:) = 0.d0
3069
3070 ! Second pass to fill covariance Matrix
3071 call obs_set_current_header_list(obsSpaceData, 'TO')
3072 iobs = 0
3073 HEADER2: do
3074 headerIndex = obs_getHeaderIndex(obsSpaceData)
3075 if (headerIndex < 0) exit HEADER2
3076
3077 ! process only radiance data to be assimilated?
3078 idatyp = obs_headElem_i(obsSpaceData, OBS_ITY, headerIndex)
3079 if (.not. tvs_isIdBurpTovs(idatyp)) then
3080 write(*,*) 'bcs_outputCvOmPPred: warning unknown radiance codtyp present check NAMTOVSINST', idatyp
3081 cycle HEADER2
3082 end if
3083 iobs = iobs + 1
3084 indxtovs = tvs_tovsIndex(headerIndex)
3085 if (indxtovs < 0) cycle HEADER2
3086
3087 iSensor = tvs_lsensor(indxTovs)
3088 if (iSensor /= sensorIndex) cycle HEADER2
3089
3090 call obs_set_current_body_list(obsSpaceData, headerIndex)
3091 BODY2: do
3092 bodyIndex = obs_getBodyIndex(obsSpaceData)
3093 if (bodyIndex < 0) exit BODY2
3094 if (obs_bodyElem_i(obsSpaceData, OBS_ASS, bodyIndex) /= obs_assimilated) cycle BODY2
3095 call bcs_getChannelIndex(obsSpaceData, iSensor, chanIndx, bodyIndex)
3096 if (chanIndx > 0) then
3097 call bcs_getPredictors(predictor, headerIndex, iobs, chanIndx, obsSpaceData)
3098 predictor(:) = predictor(:) - predBiasMpiGLobal(chanIndx,:)
3099 OmF = obs_bodyElem_r(obsSpaceData, OBS_OMP, bodyIndex)
3100 OmF = OmF - omfBiasMpiGlobal(chanIndx)
3101 vector(1,1) = OmF
3102 vector(1,2:numPredictors) = predictor(2:numPredictors)
3103 Matrix(chanindx,:,:) = Matrix(chanindx,:,:) + matmul(transpose(vector),vector)
3104 end if
3105 end do BODY2
3106 end do HEADER2
3107
3108 if (mmpi_myId == 0) then
3109 allocate(matrixMpiGlobal(nchans,numPredictors,numPredictors))
3110 else
3111 allocate(matrixMpiGlobal(1,1,1))
3112 end if
3113
3114 ! communication MPI pour tout avoir sur tache 0
3115 call mmpi_reduce_sumR8_3d(matrix, matrixMpiGlobal, 0, "GRID" )
3116 deallocate(matrix)
3117 deallocate(OmFBiasMpiGLobal)
3118 deallocate(predBiasMpiGLobal)
3119 if (mmpi_myId == 0) then
3120 call utl_open_asciifile("covarianceMatrix_" // trim(tvs_instrumentName(sensorIndex)) // &
3121 "_" // trim(tvs_satelliteName(sensorIndex)) // ".dat", iuncov)
3122 call utl_open_asciifile("correlationMatrix_" // trim(tvs_instrumentName(sensorIndex)) // &
3123 "_" // trim(tvs_satelliteName(sensorIndex)) // ".dat", iuncorr)
3124 do channelIndex = 1, nchans
3125 if (countMpiGlobal(channelIndex) > 1) then
3126 matrixMpiGlobal(channelIndex,:,:) = matrixMpiGlobal(channelIndex,:,:) / countMpiGlobal(channelIndex)
3127
3128 write(iuncov,*) "OmF Pred covariance Matrix for channel ", &
3129 bias(sensorIndex)%chans(channelIndex)%channelNum, "instrument ", &
3130 trim(tvs_instrumentName(sensorIndex))," ", &
3131 trim(tvs_satelliteName(sensorIndex))
3132 write(iuncov,'(10x,A6)',advance="no") "OmF"
3133 do predictorIndex = 2, numPredictors
3134 write(iuncov,'(T6,A6,1x)',advance="no") predTab(predictorIndex)
3135 end do
3136 write(iuncov,*)
3137 write(iuncov,'(A6)',advance="no") "Omf"
3138 write(iuncov,'(100f12.6)') matrixMpiGlobal(channelIndex,1,:)
3139 do predictorIndex = 2, numPredictors
3140 write(iuncov,'(A6)',advance="no") predTab(predictorIndex)
3141 write(iuncov,'(100f12.6)') matrixMpiGlobal(channelIndex,predictorIndex,:)
3142 end do
3143 sigma(:) = 0
3144 do predictorIndex = 1,numPredictors
3145 if ( matrixMpiGlobal(channelIndex,predictorIndex,predictorIndex) >0.d0) &
3146 sigma(predictorIndex) = sqrt(matrixMpiGlobal(channelIndex,predictorIndex,predictorIndex))
3147 end do
3148 do predictorIndex = 1, numPredictors
3149 do predictorIndex2 =1, numPredictors
3150 correlation(predictorIndex, predictorIndex2) = &
3151 matrixMpiGlobal(channelIndex,predictorIndex,predictorIndex2) / &
3152 (sigma(predictorIndex) * sigma(predictorIndex2) )
3153 end do
3154 end do
3155 write(iuncorr,*) "OmF Pred correlation Matrix for channel ", &
3156 bias(sensorIndex)%chans(channelIndex)%channelNum,"instrument ", &
3157 trim(tvs_instrumentName(sensorIndex))," ", &
3158 trim(tvs_satelliteName(sensorIndex))
3159 write(iuncorr,'(10x,A6)',advance="no") "OmF"
3160 do predictorIndex = 2, numPredictors
3161 write(iuncorr,'(T6,A6,1x)',advance="no") predTab(predictorIndex)
3162 end do
3163 write(iuncorr,*)
3164 write(iuncorr,'(A6)',advance="no") "Omf"
3165 write(iuncorr,'(100f12.6)') correlation(1,:)
3166 do predictorIndex = 2, numPredictors
3167 write(iuncorr,'(A6)',advance="no") predTab(predictorIndex)
3168 write(iuncorr,'(100f12.6)') correlation(predictorIndex,:)
3169 end do
3170 end if
3171 end do
3172 ierr = fclos(iuncov)
3173 ierr = fclos(iuncorr)
3174 end if
3175
3176 deallocate(countMpiGlobal)
3177 deallocate(matrixMpiGlobal)
3178
3179 end do SENSORS
3180
3181 end subroutine bcs_outputCvOmPPred
3182
3183 !----------------------
3184 ! bcs_Finalize
3185 !----------------------
3186 subroutine bcs_Finalize
3187 !
3188 !:Purpose: release allocated memory for the module
3189 !
3190 implicit none
3191
3192 ! Locals:
3193 integer :: iSensor, iChannel
3194
3195 if (.not. biasActive) return
3196
3197 if (allocated(trialHeight300m1000)) deallocate(trialHeight300m1000)
3198 if (allocated(trialHeight50m200)) deallocate(trialHeight50m200)
3199 if (allocated(trialHeight1m10)) deallocate(trialHeight1m10)
3200 if (allocated(trialHeight5m50)) deallocate(trialHeight5m50)
3201 if (allocated(trialTG)) deallocate(trialTG)
3202 if (allocated(RadiosondeWeight)) deallocate(RadiosondeWeight)
3203
3204 do iSensor = 1, tvs_nSensors
3205 if (allocated(bias(iSensor)%BHalfScanBias)) &
3206 deallocate(bias(iSensor)%BHalfScanBias)
3207 if (allocated(bias(iSensor)%BMinusHalfScanBias)) &
3208 deallocate(bias(iSensor)%BMinusHalfScanBias)
3209 do iChannel =1, bias(iSensor)%numChannels
3210 deallocate(bias(iSensor)%chans(iChannel)%stddev)
3211 deallocate(bias(iSensor)%chans(iChannel)%coeffIncr)
3212 deallocate(bias(iSensor)%chans(iChannel)%predictorIndex)
3213 if (allocated(bias(iSensor)%chans(iChannel)%coeffIncr_fov)) deallocate(bias(iSensor)%chans(iChannel)%coeffIncr_fov)
3214 deallocate(bias(iSensor)%chans(iChannel)%coeff_offset)
3215 if (allocated(bias(iSensor)%chans(iChannel)%coeff)) deallocate(bias(iSensor)%chans(iChannel)%coeff)
3216 if (allocated(bias(iSensor)%chans(iChannel)%coeff_fov)) deallocate(bias(iSensor)%chans(iChannel)%coeff_fov)
3217 if (allocated(bias(iSensor)%chans(iChannel)%coeffCov)) deallocate(bias(iSensor)%chans(iChannel)%coeffCov)
3218 end do
3219 deallocate(bias(iSensor)%chans)
3220 end do
3221
3222 end subroutine bcs_Finalize
3223
3224 !-----------------------------
3225 ! InstrNametoCoeffFileName
3226 !-----------------------------
3227 function InstrNametoCoeffFileName(nameIn) result(nameOut)
3228 implicit none
3229
3230 ! Arguments:
3231 character(len=10), intent(in) :: nameIn
3232 ! Result:
3233 character(len=10) :: nameOut
3234
3235 ! Locals:
3236 character(len=10) :: temp_instrName
3237 integer :: ierr
3238
3239 temp_instrName = nameIn
3240 ierr = clib_tolower(temp_instrName)
3241 if (trim(temp_instrName) == 'mhs') then
3242 nameOut = 'amsub'
3243 else if (trim(temp_instrName) == 'goesimager') then
3244 nameOut = 'cgoes'
3245 else if (trim(temp_instrName) == 'gmsmtsat') then
3246 nameOut = 'mtsat'
3247 else if (trim(temp_instrName) == 'mviri') then
3248 nameOut = 'mets7'
3249 else
3250 nameOut = temp_instrName
3251 end if
3252
3253 end function InstrNametoCoeffFileName
3254
3255 !-----------------------------
3256 ! InstrNameinCoeffFile
3257 !-----------------------------
3258 function InstrNameinCoeffFile(nameIn) result(nameOut)
3259 implicit none
3260
3261 ! Arguments:
3262 character(len=10), intent(in) :: nameIn
3263 ! Result:
3264 character(len=10) :: nameOut
3265
3266 if (trim(nameIn) == 'MHS') then
3267 nameOut = 'AMSUB'
3268 else if (trim(nameIn) == 'GOESIMAGER') then
3269 nameOut = 'CGOES'
3270 else if (trim(nameIn) == 'GMSMTSAT') then
3271 nameOut = 'MTSAT'
3272 else if (trim(nameIn) == 'MVIRI') then
3273 nameOut = 'METS7'
3274 else
3275 nameOut = nameIn
3276 end if
3277
3278 end function InstrNameinCoeffFile
3279
3280 !-----------------------------
3281 ! SatNameinCoeffFile
3282 !-----------------------------
3283 function SatNameinCoeffFile(nameIn) result(nameOut)
3284 implicit none
3285
3286 ! Arguments:
3287 character(len=10), intent(in) :: nameIn
3288 ! Result:
3289 character(len=10) :: nameOut
3290
3291 if (trim(nameIn) == 'MSG1') then
3292 nameOut = 'METSAT8'
3293 else if (trim(nameIn) == 'MSG2') then
3294 nameOut = 'METSAT9'
3295 else if (trim(nameIn) == 'MSG3') then
3296 nameOut = 'METSAT10'
3297 else if (trim(nameIn) == 'MSG4') then
3298 nameOut = 'METSAT11'
3299 else if (trim(nameIn) == 'METEOSAT7') then
3300 nameOut = 'METSAT7'
3301 else
3302 nameOut = nameIn
3303 end if
3304
3305 end function SatNameinCoeffFile
3306
3307 !-----------------------------------------
3308 ! read_bcif
3309 !-----------------------------------------
3310 subroutine read_bcif(bcifFile, hspec, ncan, can, bcmode, bctype, npred, pred, global_opt, exitcode)
3311 !
3312 !:Purpose: to read channel-specific bias correction (BC) information (predictors) for instrument from BCIF.
3313 !
3314 implicit none
3315
3316 ! Arguments:
3317 character(len=*), intent(in) :: bcifFile
3318 logical, intent(in) :: hspec
3319 integer, intent(out) :: exitcode
3320 integer, intent(out) :: ncan
3321 integer, intent(out) :: can(tvs_maxchannelnumber)
3322 integer, intent(out) :: npred(tvs_maxchannelnumber)
3323 character(len=3), intent(in) :: global_opt
3324 character(len=1), intent(out) :: bcmode(tvs_maxchannelnumber)
3325 character(len=1), intent(out) :: bctype(tvs_maxchannelnumber)
3326 character(len=2), intent(out) :: pred(tvs_maxchannelnumber,numpredictorsBCIF)
3327
3328 ! Locals:
3329 character(len=7) :: instrum
3330 integer :: i, j, ier, ii, iun
3331 character(len=64) :: line
3332 integer :: xcan, xnpred, chknp
3333 character(len=1) :: xbcmode, xbctype
3334 character(len=2) :: xpred(numpredictors)
3335
3336 ! Reads channel-specific bias correction (BC) information (predictors) for instrument from BCIF.
3337 ! Channel 0 values are global or default values (optionally applied to all channels).
3338 ! Returns BC information for all channels to calling routine.
3339 !
3340 ! Sample BCIF
3341 ! AMSUA 15 <--- instrument and number of channels
3342 !CHAN MODE TYPE NPRED PRED1 PRED2 PRED3 PRED4 PRED5 PRED6
3343 ! 0 D C 4 T1 T2 T3 T4 XX XX <--- channel "0" global or default values
3344 ! 1 D C 2 T1 T2 XX XX XX XX
3345 ! 2 D C 3 BT T1 T2 XX XX XX
3346 ! 3 S C 2 T3 T4 XX XX XX XX
3347 ! ....
3348 ! ....
3349 ! =================== 24 APRIL 2014 LIST-DIRECTED I/O VERSION ==============================================
3350 ! CALL read_bcif(iunbc, bc_instrum, bc_ncan, bc_can, bc_mode, bc_type, bc_npred, bc_pred, global_opt, exitcode)
3351 !
3352 ! global_opt = NON Read channel-specific data for ALL ncan channels from BCIF (channel 0 ignored)
3353 ! OUI Read channel 0 data and apply to all ncan channels (global values)
3354 ! DEF Read channel 0 data and apply as default values for all ncan channels;
3355 ! then scan the rest of the BCIF for any channel-specific data that will
3356 ! override the default values.
3357 !
3358 ! NOTE: For hyperspectral instruments (e.g. AIRS, IASI, CrIS) the BCIF must always contain records for ALL ncan channels.
3359 ! If global_opt = OUI, only the channel numbers are needed in column 1 (CHAN) to get the list
3360 ! of channel numbers.
3361 ! If global_opt = DEF, other column data (MODE, TYPE, NPRED, PRED1,...) are entered only for
3362 ! those channels for which the default (channel 0) values are to be overridden.
3363 !
3364 ! For standard instruments (AMSU, SSM/I, SSMIS), with consecutive channels 1,2,3,...,ncan:
3365 ! If global_opt = OUI, only the channel 0 record is needed in the BCIF (any other records after
3366 ! channel 0 are ignored (not read).
3367 ! If global_opt = DEF, the channel 0 record (default values) and only records for those channels
3368 ! for which values are different from defaults are needed in the BCIF.
3369
3370 exitcode = -1
3371
3372 iun = 0
3373 ier = fnom(iun, bcifFile, 'FMT', 0)
3374 if (ier /= 0) then
3375 call utl_abort('read_bcif: ERROR - Problem opening the bcif file!' // trim(bcifFile))
3376 end if
3377
3378 read(iun,*) instrum, ncan
3379 read(iun,'(A64)') line
3380
3381 ! For GLOBAL option, read global values from first line (channel 0) and clone to all channels
3382 if (global_opt == 'OUI' .or. global_opt == 'DEF') then
3383 ! Read channel 0 information
3384 read(iun, *, iostat=ier) can(1), bcmode(1), bctype(1), npred(1), (pred(1,j), j = 1, numpredictorsbcif)
3385 if (ier /= 0) then
3386 write(*,*) 'read_BCIF: Error reading channel 0 data!'
3387 exitcode = ier
3388 return
3389 end if
3390 if (can(1) /= 0) then
3391 write(*,*) 'read_BCIF: Channel 0 global values not found!'
3392 exitcode = -1
3393 return
3394 end if
3395 ! Clone channel 0 information to all ncan channels
3396 if (.not. hspec) then
3397 ! For instruments with consecutive channels 1,2,3,...ncan (e.g. AMSU, SSM/I)
3398 ! -- no need to read the channel numbers from the BCIF
3399 do i = 2, ncan+1
3400 can(i) = i - 1
3401 bcmode(i) = bcmode(1)
3402 bctype(i) = bctype(1)
3403 npred(i) = npred(1)
3404 do j = 1, numpredictorsbcif
3405 pred(i,j) = pred(1,j)
3406 end do
3407 end do
3408 else
3409 ! For hyperspectral instruments (channel subsets), read the channel numbers from the BCIF
3410 do i = 2, ncan + 1
3411 read(iun, *, iostat=ier) xcan
3412 if (ier /= 0) then
3413 write(*,*) 'read_BCIF: Error reading channel numbers!'
3414 exitcode = ier
3415 return
3416 end if
3417 can(i) = xcan
3418 bcmode(i) = bcmode(1)
3419 bctype(i) = bctype(1)
3420 npred(i) = npred(1)
3421 do j = 1, numpredictorsbcif
3422 pred(i,j) = pred(1,j)
3423 end do
3424 end do
3425 ! Reposition the file to just after channel 0 record
3426 rewind (iun)
3427 read(iun,*) instrum, ncan
3428 read(iun,'(A64)') line
3429 read(iun,*,iostat=ier) xcan, xbcmode, xbctype, xnpred, (xpred(j), j = 1, numpredictorsbcif)
3430 end if
3431 ! For global_opt == 'DEF' check for channel-specific information and overwrite the default (channel 0) values
3432 ! for the channel with the values from the file
3433 if (global_opt == 'DEF') then
3434 if (.not. hspec) then
3435 do
3436 read(iun,*,iostat=ier) xcan, xbcmode, xbctype, xnpred, (xpred(j), j = 1, numpredictorsbcif)
3437 if (ier < 0) exit
3438 if (ier > 0) then
3439 write(*,*) 'read_BCIF: Error reading file!'
3440 exitcode = ier
3441 return
3442 end if
3443 ii = xcan + 1
3444 if (ii > ncan + 1) then
3445 write(*,*) 'read_BCIF: Channel number in BCIF exceeds number of channels!'
3446 write(*,'(A,1X,I4,1X,I4)') ' Channel, ncan = ', xcan, ncan
3447 exitcode = -1
3448 return
3449 end if
3450 bcmode(ii) = xbcmode
3451 bctype(ii) = xbctype
3452 npred(ii) = xnpred
3453 do j = 1, numpredictorsbcif
3454 pred(ii,j) = xpred(j)
3455 end do
3456 end do
3457 else
3458 ! For hyperspectral instruments
3459 do i = 2, ncan + 1
3460 read(iun, *,iostat=ier) xcan, xbcmode, xbctype, xnpred, (xpred(j), j = 1, numpredictorsbcif)
3461 if (ier /= 0) cycle
3462 bcmode(i) = xbcmode
3463 bctype(i) = xbctype
3464 npred(i) = xnpred
3465 do j = 1, numpredictorsbcif
3466 pred(i,j) = xpred(j)
3467 end do
3468 end do
3469 end if
3470 end if
3471 end if
3472 ! Non-GLOBAL: Read the entire file for channel specific values (all channels)
3473 ! ---------------------------------------------------------------------------------------------------------------------
3474 if (global_opt == 'NON') then
3475 ii = 1
3476 do
3477 read(iun,*,iostat=ier) can(ii), bcmode(ii), bctype(ii), npred(ii), (pred(ii,j), j = 1, numpredictorsbcif)
3478 if (ier < 0) exit
3479 if (ier > 0) then
3480 write(*,*) 'read_BCIF: Error reading file!'
3481 exitcode = ier
3482 return
3483 end if
3484 if (ii == 1) then
3485 if (can(ii) /= 0) then
3486 write(*,*) 'read_BCIF: Channel 0 global/default values not found!'
3487 exitcode = -1
3488 return
3489 end if
3490 end if
3491 ii = ii + 1
3492 end do
3493 if (ii - 2 < ncan) then
3494 write(*,*) 'read_BCIF: Number of channels in file is less than specified value (NCAN). Changing value of NCAN.'
3495 ncan = ii - 2
3496 end if
3497 if (ii > ncan + 2) then
3498 write(*,*) 'read_BCIF: ERROR -- Number of channels in file is greater than specified value (NCAN)!'
3499 exitcode = -1
3500 return
3501 end if
3502 end if
3503 ! ---------------------------------------------------------------------------------------------------------------------
3504 write(*,*) ' '
3505 write(*,*) 'read_BCIF: Bias correction information for each channel (from BCIF):'
3506 write(*,'(1X,A7,1X,I4)') instrum, ncan
3507 write(*,*) line
3508 do i = 1, ncan + 1
3509 chknp = count(pred(i,:) /= 'XX')
3510 if (chknp /= npred(i)) npred(i) = chknp
3511 if (npred(i) == 0 .and. bctype(i) == 'C') bctype(i) = 'F'
3512 write(*,'(I4,2(5X,A1),5X,I2,6(4X,A2))') can(i), bcmode(i), bctype(i), npred(i), (pred(i,j), j = 1, numpredictorsbcif)
3513 end do
3514 write(*,*) 'read_BCIF: exit '
3515
3516 ier = fclos(iun)
3517 exitcode = 0
3518
3519 end subroutine read_bcif
3520
3521 !-----------------------------------------
3522 ! read_coeff
3523 !-----------------------------------------
3524 subroutine read_coeff(sats, chans, fovbias, coeff, nsat, nchan, nfov, npred, cinstrum, coeff_file, ptypes, ndata)
3525 !
3526 !:Purpose: to read radiance bias correction coefficients file
3527 !
3528 implicit none
3529
3530 ! Arguments:
3531 character(len=10), intent(out) :: sats(:) ! dim(maxsat), satellite names 1
3532 integer, intent(out) :: chans(:,:) ! dim(maxsat, maxchan), channel numbers 2
3533 real(8), intent(out) :: fovbias(:,:,:)! dim(maxsat,maxchan,maxfov), bias as F(fov) 3
3534 real(8), intent(out) :: coeff(:,:,:) ! dim(maxsat,maxchan,maxpred+1) 4
3535 integer, intent(out) :: nsat !5
3536 integer, intent(out) :: nchan(:) ! dim(maxsat), number of channels 6
3537 integer, intent(out) :: nfov !7
3538 integer, intent(out) :: npred(:,:) ! dim(maxsat, maxchan), number of predictors !8
3539 character(len=7), intent(out) :: cinstrum ! string: instrument (e.g. AMSUB) 9
3540 character(len=*), intent(in) :: coeff_file ! 10
3541 character(len=2), intent(out) :: ptypes(:,:,:) ! dim(maxsat,maxchan,maxpred) 11
3542 integer, intent(out) :: ndata(:,:) ! dim(maxsat, maxchan), number of channels 12
3543
3544 ! Locals:
3545 character(len=8) :: sat
3546 character(len=120) :: line
3547 integer :: chan
3548 integer :: nbfov, nbpred, i, j, k, ier, istat, ii, nobs
3549 logical :: newsat, fileExists
3550 real :: dummy
3551 integer :: iun
3552 integer :: maxsat
3553 integer :: maxpred
3554
3555 ! sats(nsat) = satellite names
3556 ! chans(nsat,nchan(i)) = channel numbers of each channel of each satellite i
3557 ! npred(nsat,nchan(i)) = number of predictors for each channel of each satellite i
3558 ! fovbias(i,j,k) = bias for satellite i, channel j, FOV k k=1,nfov
3559 ! if FOV not considered for instrument, nfov = 1 and fovbias is global bias for channel
3560 ! coeff(i,j,1) = regression constant
3561 ! coeff(i,j,2), ..., coeff(i,j,npred(i,j)) = predictor coefficients
3562 ! nsat, nchan, nfov, cinstrum (output) are determined from file
3563 ! if returned nsat = 0, coeff_file was empty
3564
3565 inquire(file=trim(coeff_file), exist=fileExists)
3566 if (fileExists) then
3567 iun = 0
3568 ier = fnom(iun, coeff_file, 'FMT', 0)
3569 if (ier /= 0) then
3570 call utl_abort('read_coeff: ERROR - Problem opening the coefficient file! ' // trim(coeff_file))
3571 end if
3572
3573 write(*,*)
3574 write(*,*) 'read_coeff: Bias correction coefficient file open = ', coeff_file
3575
3576 end if
3577
3578 maxsat = size(sats)
3579 maxpred = size(ptypes, dim = 3)
3580
3581 coeff(:,:,:) = 0.0
3582 fovbias(:,:,:) = 0.0
3583 sats(:) = 'XXXXXXXX'
3584 cinstrum = 'XXXXXXX'
3585 chans(:,:) = 0
3586 npred(:,:) = 0
3587 nsat = 0
3588 nchan(:) = 0
3589 nfov = 0
3590 ptypes(:,:,:) = 'XX'
3591
3592 if (.not. fileExists) then
3593 write(*,*) 'read_coeff: Warning- File ' // trim(coeff_file) //'not there.'
3594 return
3595 end if
3596
3597 read(iun,*,iostat=istat)
3598
3599 if (istat /= 0) then
3600 write(*,*) 'read_coeff: Warning- File appears empty.'
3601 ier = fclos(iun)
3602 return
3603 end if
3604
3605 rewind(iun)
3606
3607 ii = 0
3608
3609 ! Loop over the satellites/channels in the file
3610
3611 do
3612 read(iun,'(A)',iostat=istat) line
3613 if (istat < 0) exit
3614 if (line(1:3) == 'SAT') then
3615 newsat = .true.
3616 read(line,'(T53,A8,1X,A7,1X,I6,1X,I8,1X,I2,1X,I3)',iostat=istat) sat, cinstrum, chan, nobs, nbpred, nbfov
3617 if (istat /= 0) then
3618 call utl_abort('read_coeff: ERROR - reading data from SATELLITE line in coeff file!')
3619 end if
3620 do i = 1, maxsat
3621 if (trim(sats(i)) == trim(sat)) then
3622 newsat = .false.
3623 ii = i
3624 end if
3625 end do
3626 if (newsat) then
3627 ii = ii + 1
3628 if (ii > maxsat) then
3629 call utl_abort('read_coeff: ERROR - max number of satellites exceeded in coeff file!')
3630 end if
3631 sats(ii) = sat
3632 if (ii > 1) nchan(ii - 1) = j
3633 j = 1
3634 else
3635 j = j + 1
3636 end if
3637 chans(ii,j) = chan
3638 npred(ii,j) = nbpred
3639 ndata(ii,j) = nobs
3640 if (nbpred > maxpred) then
3641 call utl_abort('read_coeff: ERROR - max number of predictors exceeded in coeff file!')
3642 end if
3643 read(iun,'(A)',iostat=istat) line
3644 if (line(1:3) /= 'PTY') then
3645 call utl_abort('read_coeff: ERROR - list of predictors is missing in coeff file!')
3646 end if
3647 if (nbpred > 0) then
3648 read(line,'(T8,6(1X,A2))', iostat = istat) (ptypes(ii,j,k), k = 1, nbpred)
3649 if (istat /= 0) then
3650 call utl_abort('read_coeff: ERROR - reading predictor types from PTYPES line in coeff file!')
3651 end if
3652 end if
3653 read(iun,*,iostat=istat) (fovbias(ii,j,k), k= 1, nbfov)
3654 if (istat /= 0) then
3655 call utl_abort('read_coeff: ERROR - reading fovbias in coeff file!')
3656 end if
3657 if (nbpred > 0) then
3658 read(iun,*,iostat=istat) (coeff(ii,j,k), k = 1, nbpred + 1)
3659 else
3660 read(iun,*,iostat=istat) dummy
3661 end if
3662 if (istat /= 0) then
3663 call utl_abort('read_coeff: ERROR - reading coeff in coeff file!')
3664 end if
3665 else
3666 exit
3667 end if
3668
3669 end do
3670
3671 if (ii == 0) then
3672 call utl_abort('read_coeff: ERROR - No data read from coeff file!')
3673 end if
3674
3675 nsat = ii
3676 nfov = nbfov
3677 nchan(ii) = j
3678
3679 write(*,*) ' '
3680 write(*,*) 'read_coeff: ------------- BIAS CORRECTION COEFFICIENT FILE ------------------ '
3681 write(*,*) ' '
3682 write(*,*) 'read_coeff: Number of satellites = ', nsat
3683 write(*,*) 'read_coeff: Number of FOV = ', nfov
3684 write(*,*) 'read_coeff: Max number of predictors = ', maxval(npred)
3685 write(*,*) ' '
3686 do i = 1, nsat
3687 write(*,*) 'read_coeff: Satellite = ' // sats(i)
3688 write(*,*) 'read_coeff: Number of channels = ', nchan(i)
3689 write(*,*) 'read_coeff: predictors, fovbias, coeff for each channel: '
3690 do j = 1, nchan(i)
3691 write(*,*) i, chans(i, j)
3692 if (npred(i, j) > 0) then
3693 write(*,'(6(1X,A2))') (ptypes(i,j,k), k = 1, npred(i,j))
3694 else
3695 write(*,'(A)') 'read_coeff: No predictors'
3696 end if
3697 write(*,*) (fovbias(i,j,k), k = 1, nfov)
3698 write(*,*) (coeff(i,j,k), k = 1, npred(i,j) + 1)
3699 end do
3700 end do
3701 write(*,*) 'read_coeff: exit'
3702
3703 ier = fclos(iun)
3704
3705 end subroutine read_coeff
3706
3707 !-----------------------------------------
3708 ! bcs_getChannelIndex
3709 !-----------------------------------------
3710 subroutine bcs_getChannelIndex(obsSpaceData, idsat, chanIndx,indexBody)
3711 !
3712 !:Purpose: to get the channel index (wrt bcif channels)
3713 !
3714 implicit none
3715
3716 ! Arguments:
3717 integer, intent(in) :: idsat
3718 integer, intent(in) :: indexBody
3719 integer, intent(out) :: chanIndx
3720 type(struct_obs), intent(inout) :: obsSpaceData
3721
3722 ! Locals:
3723 logical, save :: first =.true.
3724 integer :: ichan, isensor, indx
3725 integer, allocatable, save :: Index(:,:)
3726
3727 if (first) then
3728 allocate(Index(tvs_nsensors,tvs_maxChannelNumber))
3729 Index(:,:) = -1
3730 do isensor = 1, tvs_nsensors
3731 channels:do ichan = 1, tvs_maxChannelNumber
3732 indexes: do indx =1, bias(isensor)%numChannels
3733 if (ichan == bias(isensor)%chans(indx)%channelNum) then
3734 Index(isensor,ichan) = indx
3735 exit indexes
3736 end if
3737 end do indexes
3738 end do channels
3739 end do
3740 first = .false.
3741 end if
3742
3743 ichan = nint(obs_bodyElem_r(obsSpaceData, OBS_PPP, indexBody))
3744 ichan = max(0, min(ichan, tvs_maxChannelNumber + 1))
3745 ichan = ichan - tvs_channelOffset(idsat)
3746
3747 chanIndx = Index(idsat,ichan)
3748
3749 end subroutine bcs_getChannelIndex
3750
3751 !-----------------------------------------
3752 ! getObsFileName
3753 !-----------------------------------------
3754 function getObsFileName(codetype) result(fileName)
3755 !
3756 !:Purpose: Return the part of the observation file name associated
3757 ! with the type of observation it contains.
3758 !
3759 implicit none
3760
3761 ! Arguments:
3762 integer, intent(in) :: codeType
3763 ! Result:
3764 character(len=20) :: fileName
3765
3766 if (codtyp_get_name(codeType) == 'radianceclear') then
3767 fileName = 'csr'
3768 else if (codtyp_get_name(codeType) == 'mhs' .or. codtyp_get_name(codeType) == 'amsub') then
3769 fileName = 'to_amsub'
3770 else if (codtyp_get_name(codeType) == 'amsua') then
3771 fileName = 'to_amsua'
3772 else if (codtyp_get_name(codeType) == 'ssmi') then
3773 fileName = 'ssmis'
3774 else if (codtyp_get_name(codeType) == 'crisfsr') then
3775 fileName = 'cris'
3776 else
3777 fileName = codtyp_get_name(codeType)
3778 end if
3779
3780 end function getObsFileName
3781
3782 !-----------------------------------------
3783 ! getInitialIdObsData
3784 !-----------------------------------------
3785 subroutine getInitialIdObsData(obsSpaceData, obsFamily, idObs, idData, codeTypeList_opt)
3786 !
3787 !:Purpose: Compute initial value for idObs and idData that will ensure
3788 ! unique values over all mpi tasks
3789 !
3790 implicit none
3791
3792 ! Arguments:
3793 type(struct_obs), intent(inout) :: obsSpaceData
3794 character(len=*), intent(in) :: obsFamily
3795 integer, intent(out) :: idObs
3796 integer, intent(out) :: idData
3797 integer, optional, intent(in) :: codeTypeList_opt(:)
3798
3799 ! Locals:
3800 integer :: headerIndex, numHeader, numBody, codeType, ierr
3801 integer, allocatable :: allNumHeader(:), allNumBody(:)
3802
3803 numHeader = 0
3804 numBody = 0
3805 call obs_set_current_header_list(obsSpaceData, obsFamily)
3806 HEADERCOUNT: do
3807 headerIndex = obs_getHeaderIndex(obsSpaceData)
3808 if (headerIndex < 0) exit HEADERCOUNT
3809 if (present(codeTypeList_opt)) then
3810 codeType = obs_headElem_i(obsSpaceData, OBS_ITY, headerIndex)
3811 if (all(codeTypeList_opt(:) /= codeType)) cycle HEADERCOUNT
3812 end if
3813 numHeader = numHeader + 1
3814 numBody = numBody + obs_headElem_i(obsSpaceData, OBS_NLV, headerIndex)
3815 end do HEADERCOUNT
3816 allocate(allNumHeader(mmpi_nprocs))
3817 allocate(allNumBody(mmpi_nprocs))
3818 call rpn_comm_allgather(numHeader, 1, 'mpi_integer', &
3819 allNumHeader, 1, 'mpi_integer', 'GRID', ierr)
3820 call rpn_comm_allgather(numBody, 1, 'mpi_integer', &
3821 allNumBody, 1, 'mpi_integer', 'GRID', ierr)
3822 if (mmpi_myid > 0) then
3823 idObs = sum(allNumHeader(1:mmpi_myid))
3824 idData = sum(allNumBody(1:mmpi_myid))
3825 else
3826 idObs = 0
3827 idData = 0
3828 end if
3829 deallocate(allNumHeader)
3830 deallocate(allNumBody)
3831
3832 end subroutine getInitialIdObsData
3833
3834end module biasCorrectionSat_mod
