1program midas_obsSelection
2 !
3 !:Purpose: Main program for background check, and thinning of all observation types.
4 !
5 ! ---
6 !
7 !:Algorithm: For background check of each observation type, several quality control tests
8 ! specific to that observation type are performed on each available observation
9 ! to determine if the observation meets the standard to be assimilated later.
10 ! Observation flags are modified for rejected observations as well as for the
11 ! assimilated observations of inferior quality.The thinning is performed afterwards
12 ! to reduce the number of observation data for assimilation. This is to reduce the
13 ! 1) computational cost, and 2) observation error correlation during
14 ! assimilation stage. The background-checked thinned observations are written
15 ! to new files, ready for assimilation.
16 !
17 ! --
18 !
19 ! The computed bias correction values are applied to the observation before the
20 ! background check step for certain observation types (e.g. radiances). One of quality
21 ! control tests during background check is the ``rogue check`` which is the allowed
22 ! distance of the observation from the background state. Innovation vector
23 ! ``y-H(xb)`` is required to measure this distance. The innovations are
24 ! computed at the beginning of the program before the background check starts.
25 !
26 ! --
27 !
28 !============================================== ==============================================================
29 ! Input and Output Files Description of file
30 !============================================== ==============================================================
31 ! ``flnml`` In - Main namelist file with parameters that user may modify
32 ! ``flnml_static`` In - The "static" namelist that should not be modified
33 ! ``trlm_$NN`` (e.g. ``trlm_01``) In - Background state (a.k.a. trial) files for each timestep
34 ! ``analysisgrid`` In - File defining grid for computing the analysis increment
35 ! ``obsfiles_$FAM/obs$FAM_$NNNN_$NNNN`` In - Observation file for each "family" and MPI task
36 ! ``obserr`` In - Observation error statistics
37 ! ``obsfiles_$FAM.updated/obs$FAM_$NNNN_$NNNN`` Out - Updated obs file for each "family" and MPI task
38 ! Remainder are files related to radiance obs:
39 ! ``stats_$SENSOR_assim`` In - Satellite radiance observation errors of difference sensors
40 ! ``stats_tovs`` In - Satellite radiance observation errors
41 ! ``stats_tovs_symmetricObsErr`` In - User-defined symmetric TOVS errors for all sky
42 ! ``Cmat_$PLATFORM_$SENSOR.dat`` In - Inter-channel observation-error correlations
43 ! ``dynbcor.coeffs.$SENSOR.*.coeffs_$SENSOR`` In - Dynamic bias correction file
44 ! ``ceres_global.std`` In - High-res surface type and water fraction for radiance obs
45 ! ``rtcoef_$PLATFORM_$SENSOR.dat`` In - RTTOV coefficient files
46 ! ``rttov_h2o_limits.dat`` In - Min/max humidity limits applied to analysis
47 ! ``ozoneclim98`` In - Ozone climatology
48 !============================================== ==============================================================
49 !
50 ! --
51 !
52 !:Synopsis: Below is a summary of the ``obsSelection`` program calling sequence:
53 !
54 ! - **Initial setups:**
55 !
56 ! - Read the NAMOBSSELECTION namelist and check/modify some values.
57 !
58 ! - Various modules are setup: ``obsFiles_mod``, ``timeCoord_mod``.
59 !
60 ! - Setup horizontal and vertical grid objects for "analysis
61 ! grid" from ``analysisgrid`` file.
62 !
63 ! - Setup ``obsSpaceData`` object and read observations from
64 ! files: ``inn_setupObs``.
65 !
66 ! - Compute and update the stored surface type for some satellite
67 ! radiance instruments.
68 !
69 ! - Setup ``columnData`` module (read list of analysis variables
70 ! from namelist) and allocate column object for storing trial
71 ! on analysis levels.
72 !
73 ! - Setup the observation error statistics in ``obsSpaceData``
74 ! object: ``oer_setObsErrors``
75 !
76 ! - Setup ``gridStateVector`` module to initialize the gridstatevector
77 ! objects.
78 !
79 ! - Applying optional bias corrections to some observation types.
80 !
81 ! - Setup horizontal and vertical grid objects for "trial grid" from
82 ! first trial file: ``trlm_01``.
83 !
84 ! - Allocate a stateVector object on the trial grid and then
85 ! read the trials: ``gio_readTrials``.
86 !
87 ! - **Computation**
88 !
89 ! - Compute ``columnTrlOnTrlLev`` and ``columnTrlOnAnlIncLev`` from
90 ! background state: ``inn_setupColumnsOnTrlLev``,
91 ! ``inn_setupColumnsOnAnlIncLev``.
92 !
93 ! - Compute innovation from background state: ``inn_computeInnovation``.
94 !
95 ! - Do background check for conventional observation:
96 ! ``bgck_bgCheck_conv``.
97 !
98 ! - Update radiance bias correction in ``obsSpaceData`` and apply
99 ! the bias corrections to the observations and innovations for
100 ! radiances: ``bcs_calcBias``, ``bcs_applyBiasCorrection``.
101 !
102 ! - Perform background check for multiple observation types.
103 !
104 ! - If thinning was requested:
105 !
106 ! - add some cloud parameters and set missing observation flags in
107 ! observation files (specific for radiances) and
108 !
109 ! - perform thinning for different observation types.
110 !
111 ! - Write the final background-checked bias corrected results (either
112 ! thinned or not thinned) into the observation file.
113 !
114 ! - If thinning was requested, remove observations which were flagged
115 ! not to be assimilated from the observation file.
116 !
117 ! --
118 !
119 !:Options: `List of namelist blocks <../namelists_in_each_program.html#obsSelection>`_
120 ! that can affect the ``obsSelection`` program.
121 !
122 ! * The use of ``obsSelection`` program is controlled by the namelist block
123 ! ``&NAMOBSSELECTION`` read by the ``obsSelection`` program.
124 !
125 ! * Some of the other relevant namelist blocks used to configure the
126 ! ``obsSelection`` are listed in the following table:
127 !
128 !=========================== ========================= =========================================
129 ! Module Namelist Description of what is controlled
130 !=========================== ========================= =========================================
131 ! ``midas_obsSelection`` ``NAMOBSSELECTION`` whether thinning is performed or not.
132 ! ``biasCorrectionConv_mod`` ``NAMBIASCONV`` variables to perform bias correction
133 ! for conventional observations.
134 ! ``biasCorrectionConv_mod`` ``NAMSONDETYPES`` additional variables to perform bias
135 ! correction for radiosondes conventional
136 ! observations.
137 ! ``backgroundCheck_mod`` ``NAMBGCKCONV`` variables to perform background check
138 ! for conventional observations.
139 ! ``SSTbias_mod`` ``NAMSSTBIASESTIMATE`` variables to perform bias estimation
140 ! and bias correction for satellite SST.
141 ! ``biasCorrectionSat_mod`` ``NAMBIASSAT`` variables to perform bias correction
142 ! for satellite radiances.
143 ! ``multi_ir_bgck_mod`` ``NAMBGCKIR`` Variables to perform background check
144 ! for hyperspectral infrared radiances.
145 ! ``bgckmicrowave_mod`` ``NAMBGCK`` Variables to perform background check
146 ! for microwave radiances.
147 ! ``bgckcsr_mod`` ``NAMCSR`` Variables To perform background check
148 ! for CSR radiances.
149 ! ``bgckssmis_mod`` ``NAMBGCK`` Variables to perform background check
150 ! for SSMIS radiances.
151 ! ``bgckOcean_mod`` ``NAMOCEANBGCHECK`` Variables to perform background check
152 ! for ocean data.
153 ! ``bgckOcean_mod`` ``NAMICEBGCHECK`` Variables to perform background check
154 ! for SST data.
155 ! ``burpread_mod`` ``NAMADDTOBURP`` element IDs to add to the BURP file
156 ! ``thinning_mod`` ``THIN_HYPER`` variables to perform thinning on
157 ! hyperspectral infrared radiances.
158 ! ``thinning_mod`` ``THIN_TOVS`` variables to perform thinning on
159 ! microwave radiances.
160 ! ``thinning_mod`` ``thin_csr`` variables to perform thinning on
161 ! CSR radiances.
162 ! ``thinning_mod`` ``thin_raobs`` variables to perform thinning on
163 ! radiosonde observations.
164 ! ``thinning_mod`` ``thin_scat`` variables to perform thinning on
165 ! scatterometer wind observations.
166 ! ``thinning_mod`` ``thin_aircraft`` variables to perform thinning on
167 ! aircraft observations.
168 ! ``thinning_mod`` ``thin_surface`` variables to perform thinning on
169 ! surface observations.
170 ! ``thinning_mod`` ``thin_gbgps`` variables to perform thinning on
171 ! ground-based GPS observations.
172 ! ``thinning_mod`` ``thin_gpsro`` variables to perform thinning on
173 ! GPS radio-occultation observations.
174 ! ``thinning_mod`` ``thin_aladin`` variables to perform thinning on
175 ! aladin wind observations.
176 ! ``thinning_mod`` ``thin_aladin`` variables to perform thinning on
177 ! aladin wind observations.
178 ! ``timeCoord_mod`` ``NAMTIME`` assimilation time window length,
179 ! temporal resolution of the background
180 ! state.
181 !=========================== ========================= =========================================
182 !
183 use version_mod
184 use codePrecision_mod
185 use ramDisk_mod
186 use midasMpi_mod
187 use utilities_mod
188 use mathPhysConstants_mod
189 use horizontalCoord_mod
190 use verticalCoord_mod
191 use timeCoord_mod
192 use gridStateVector_mod
193 use gridStateVectorFileIO_mod
194 use backgroundCheck_mod
195 use multiIRbgck_mod
196 use innovation_mod
197 use obsSpaceData_mod
198 use columnData_mod
199 use obsFiles_mod
200 use obsErrors_mod
201 use biasCorrectionSat_mod
202 use biasCorrectionConv_mod
203 use thinning_mod
204 use bgckmicrowave_mod
205 use bgckSSMIS_mod
206 use bgckCSR_mod
207 use bgckOcean_mod
208 use sstBias_mod
209
210 implicit none
211
212 integer :: dateStampFromObs, headerIndex, ierr, nulnam
213 type(struct_columnData),target :: columnTrlOnAnlIncLev
214 type(struct_columnData),target :: columnTrlOnTrlLev
215 type(struct_obs), target :: obsSpaceData
216 type(struct_gsv) :: stateVectorTrialHighRes
217 type(struct_hco), pointer :: hco_anl => null()
218 type(struct_vco), pointer :: vco_anl => null()
219 type(struct_hco), pointer :: hco_trl => null()
220 type(struct_vco), pointer :: vco_trl => null()
221 type(struct_hco), pointer :: hco_core => null()
222
223 logical :: allocHeightSfc
224 integer :: get_max_rss, fnom, fclos
225
226 ! Namelist variables
227 logical :: doThinning ! Control whether or not thinning is done
228
229 namelist /namObsSelection/ doThinning
230
231 call ver_printNameAndVersion('obsSelection','Obs Quality Control and Thinning')
232
233 !- 1.0 mpi
234 call mmpi_initialize
235
236 !- 1.1 timings
237 call tmg_init(mmpi_myid, 'TMG_INFO')
238 call utl_tmg_start(0,'Main')
239
240 !- 1.2 Read the namelist for obsSelection program (if it exists)
241 ! set default values for namelist variables
242 doThinning = .false.
243 if (utl_isNamelistPresent('namObsSelection', './flnml')) then
244 nulnam = 0
245 ierr = fnom(nulnam, './flnml', 'FTN+SEQ+R/O', 0)
246 if (ierr /= 0) call utl_abort('midas-obsSelection: Error opening file flnml')
247 read(nulnam, nml = namObsSelection, iostat = ierr)
248 if (ierr /= 0) call utl_abort('midas-obsSelection: Error reading namelist namObsSelection')
249 if (mmpi_myid == 0) write(*,nml = namObsSelection)
250 ierr = fclos(nulnam)
251 else
252 write(*,*)
253 write(*,*) 'midas-obsSelection: Namelist block namObsSelection is missing in the namelist.'
254 write(*,*) ' The default value will be taken.'
255 if (mmpi_myid == 0) write(*, nml = namObsSelection)
256 end if
257
258 !
259 !- Initialize the ram disk
260 !
261 call ram_setup
262
263 !
264 !- Initialize observation file names, but don't use datestamp
265 !
266 call obsf_setup(dateStampFromObs, 'bgck')
267
268 !
269 !- Initialize the Temporal grid and dateStamp from trial file
270 !
271 call tim_setup(fileNameForDate_opt = './trlm_01')
272
273 !
274 !- Initialize constants
275 !
276 if (mmpi_myid == 0) then
277 call mpc_printConstants(6)
278 call pre_printPrecisions
279 end if
280
281 !
282 !- Initialize the Analysis grid
283 !
284 if(mmpi_myid == 0) write(*,*)
285 if(mmpi_myid == 0) write(*,*) 'midas-obsSelection: Set hco parameters for analysis grid'
286 call hco_SetupFromFile(hco_anl, './analysisgrid', 'ANALYSIS', 'Analysis' ) ! IN
287
288 if ( hco_anl % global ) then
289 hco_core => hco_anl
290 else
291 !- Initialize the core (Non-Extended) analysis grid
292 if( mmpi_myid == 0) write(*,*) 'midas-obsSelection: Set hco parameters for core grid'
293 call hco_SetupFromFile( hco_core, './analysisgrid', 'COREGRID', 'AnalysisCore' ) ! IN
294 end if
295
296 !
297 !- Initialisation of the analysis grid vertical coordinate from analysisgrid file
298 !
299 call vco_SetupFromFile(vco_anl, './analysisgrid')
300
301 call col_setVco(columnTrlOnAnlIncLev, vco_anl)
302 write(*,*) 'Memory Used: ', get_max_rss()/1024,'Mb'
303
304 !
305 !- Setup and read observations
306 !
307 call inn_setupObs(obsSpaceData, hco_anl, 'ALL', 'LIKESPLITFILES', 'bgck')
308 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
309
310 ! if ssmis, compute the surface type ele and update obspacedata
311
312 call ssbg_computeSsmisSurfaceType(obsSpaceData)
313
314 ! Only for MWHS2 data and if modLSQ option is set to .true. in nambgck namelist, set values for land qualifier
315 ! indice and terrain type based on calculations
316 if (obs_famExist(obsSpaceData, 'TO')) then
317 call mwbg_computeMwhs2SurfaceType(obsSpaceData)
318 end if
319
320 !
321 !- Basic setup of columnData module
322 !
323 call col_setup
324 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
325
326 !
327 !- Memory allocation for background column data
328 !
329 call col_allocate(columnTrlOnAnlIncLev, obs_numheader(obsSpaceData), mpiLocal_opt = .true.)
330
331 !
332 !- Initialize the observation error covariances
333 !
334 call oer_setObsErrors(obsSpaceData, 'bgck')
335
336 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
337
338 !
339 ! Initialize list of analyzed variables.
340 !
341 call gsv_setup
342
343 ! Apply optional bias corrections
344 if (obs_famExist(obsSpaceData, 'TM')) then
345 call sstb_applySatelliteSSTBiasCorrection(obsSpaceData, hco_anl, &
346 vco_anl, columnTrlOnAnlIncLev)
347 end if
348
349 if (obs_famExist(obsSpaceData, 'AI')) call bcc_applyAIBcor(obsSpaceData)
350 if (obs_famExist(obsSpaceData, 'GP')) call bcc_applyGPBcor(obsSpaceData)
351 if (obs_famExist(obsSpaceData, 'UA')) call bcc_applyUABcor(obsSpaceData)
352
353 ! Reading trials
354 call inn_getHcoVcoFromTrlmFile(hco_trl, vco_trl)
355 allocHeightSfc = (vco_trl%Vcode /= 0)
356
357 call gsv_allocate(stateVectorTrialHighRes, tim_nstepobs, hco_trl, vco_trl, &
358 dateStamp_opt = tim_getDateStamp(), mpi_local_opt = .true., &
359 mpi_distribution_opt = 'Tiles', dataKind_opt = 4, &
360 allocHeightSfc_opt = allocHeightSfc, &
361 hInterpolateDegree_opt = 'LINEAR', beSilent_opt=.false.)
362 call gsv_zero(stateVectorTrialHighRes)
363 call gio_readTrials(stateVectorTrialHighRes)
364 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
365
366 ! Horizontally interpolate trials to trial columns
367 call inn_setupColumnsOnTrlLev(columnTrlOnTrlLev, obsSpaceData, hco_core, &
368 stateVectorTrialHighRes)
369 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
370
371 ! Interpolate trial columns to analysis levels and setup for linearized H
372 call inn_setupColumnsOnAnlIncLev(columnTrlOnTrlLev, columnTrlOnAnlIncLev)
373
374 ! Compute observation innovations and prepare obsSpaceData for minimization
375 call inn_computeInnovation(columnTrlOnTrlLev, obsSpaceData, analysisMode_opt = .false.)
376
377 ! 2.2 Perform the background check
378 ! The routine also calls compute_HBHT and writes to listings & obsSpaceData
379
380 ! Do the conventional data background check
381 call bgck_bgCheck_conv(columnTrlOnAnlIncLev, columnTrlOnTrlLev, hco_anl, obsSpaceData)
382
383 if (obs_famExist(obsSpaceData, 'TO')) then
384
385 ! Satellite radiance bias correction
386 call bcs_calcBias(obsSpaceData, columnTrlOnTrlLev) ! Fill in OBS_BCOR obsSpaceData column
387 ! with computed bias correction
388 call bcs_applyBiasCorrection(obsSpaceData, obs_var, 'TO') ! Apply bias correction to OBS
389 call bcs_applyBiasCorrection(obsSpaceData, obs_omp, 'TO') ! Apply bias correction to O-F
390
391 ! Do the TO background check
392 call irbg_bgCheckIR(columnTrlOnTrlLev, obsSpaceData)
393 call mwbg_bgCheckMW(obsSpaceData)
394 call csrbg_bgCheckCSR(obsSpaceData)
395 call ssbg_bgCheckSsmis(obsSpaceData)
396
397 end if
398
399 ! Do the ocean data background check
400 if (obs_famExist(obsSpaceData, 'TM')) call ocebg_bgCheckSST(obsSpaceData, tim_getDateStamp(), &
401 columnTrlOnTrlLev, hco_trl)
402
403 ! Do the sea ice data gross background check
404 if (obs_famExist(obsSpaceData, 'GL')) call ocebg_bgCheckSeaIce(obsSpaceData)
405
406 if (doThinning) then
407
408 ! Copy original obs files into another directory
409 call obsf_copyObsDirectory('./obsOriginal', direction = 'TO')
410
411 ! 2.3 Write obs files after background check, but before thinning
412 call obsf_writeFiles(obsSpaceData, writeDiagFiles_opt = .false.)
413
414 ! Add cloud parameter data to burp files (AIRS,IASI,CrIS,ATMS,AMSUA,...)
415 if (obs_famExist(obsSpaceData, 'TO')) then
416 call obsf_updateMissingObsFlags(obsSpaceData)
417 call obsf_addCloudParametersAndEmissivity(obsSpaceData)
418 end if
419
420 ! Copy the pre-thinning files into another directory
421 call obsf_copyObsDirectory('./obsBeforeThinning', direction = 'TO')
422
423 ! Copy original obs files back into usual directory
424 call obsf_copyObsDirectory('./obsOriginal', direction = 'FROM')
425
426 ! 2.4 Thinning: Set bit 11 of flag, one observation type at a time
427 call thn_thinHyper(obsSpaceData)
428 call thn_thinTovs(obsSpaceData)
429 call thn_thinCSR(obsSpaceData)
430 call thn_thinRaobs(obsSpaceData)
431 call thn_thinScat(obsSpaceData)
432 call thn_thinSatWinds(obsSpaceData)
433 call thn_thinAircraft(obsSpaceData)
434 call thn_thinSurface(obsSpaceData, 'SF') ! surface data thinning
435 if (obs_famExist(obsSpaceData, 'TM')) then
436 call thn_thinSurface(obsSpaceData, 'TM') ! SST thinning
437 call thn_thinSatSST(obsSpaceData) ! satellite SST thinning
438 end if
439 call thn_thinGbGps(obsSpaceData)
440 call thn_thinGpsRo(obsSpaceData)
441 call thn_thinAladin(obsSpaceData)
442 ! if requested, dump the thinned predictors and coefficients to sqlite
443 call bcs_dumpBiasToSqliteAfterThinning(obsSpaceData)
444
445 end if
446
447
448 ! 3 Write the final results
449
450 ! 3.1 Into the listings
451 write(*,*)
452 write(*,*) '> midas-obsSelection: printing the FIRST header and body'
453 do headerIndex = 1, min(1,obs_numHeader(obsSpaceData))
454 call obs_prnthdr(obsSpaceData, headerIndex)
455 call obs_prntbdy(obsSpaceData, headerIndex)
456 end do
457
458 ! 3.2 Into the observation files
459 write(*,*)
460 write(*,*) '> midas-obsSelection: writing to file'
461 call obsf_writeFiles(obsSpaceData)
462
463 ! Add cloud parameter data to burp files (AIRS,IASI,CrIS,...)
464 if (obs_famExist(obsSpaceData, 'TO')) then
465 call obsf_updateMissingObsFlags(obsSpaceData)
466 call obsf_addCloudParametersAndEmissivity(obsSpaceData)
467 end if
468
469 ! cleaning the observation files
470 if (doThinning) call obsf_cleanObsFiles()
471
472 !
473 ! 4. Ending
474 !
475 write(*,*)
476 write(*,*) '> midas-obsSelection: Ending'
477 call obs_finalize(obsSpaceData) ! deallocate obsSpaceData
478
479 call rpn_comm_finalize(ierr)
480
481 call utl_tmg_stop(0)
482 call tmg_terminate(mmpi_myid, 'TMG_INFO')
483
484end program midas_obsSelection
