1module obsErrors_mod
2 ! MODULE obsErrors_mod (prefix='oer' category='2. B and R matrices')
3 !
4 !:Purpose: Subroutines to set up the observation-error standard deviations.
5 !
6 use midasMpi_mod
7 use mathPhysConstants_mod
8 use obsSpaceData_mod
9 use obsSubSpaceData_mod
10 use tovsNL_mod
11 use codtyp_mod
12 use bufr_mod
13 use utilities_mod
14 use earthConstants_mod
15 use gps_mod
16 use horizontalCoord_mod
17 use verticalCoord_mod
18 use columnData_mod
19 use gridStateVector_mod
20 use gridStateVectorFileIO_mod
21 use rmatrix_mod
22 use varNameList_mod
23 use obsfiles_mod
24 use burp_module
25 use rttov_const, only: surftype_sea
26 use statetocolumn_mod
27 implicit none
28 save
29 private
30
31 ! public procedures
32 public :: oer_setObsErrors, oer_SETERRGPSGB, oer_SETERRGPSRO, oer_setErrBackScatAnisIce, oer_sw
33 public :: oer_setInterchanCorr, oer_inflateErrAllsky
34
35 ! public functions
36 public :: oer_getSSTdataParam_char, oer_getSSTdataParam_int, oer_getSSTdataParam_R8
37
38 ! public variables (parameters)
39 public :: oer_ascatAnisOpenWater, oer_ascatAnisIce
40
41 ! Arrays for QC purpose
42 public :: oer_toverrst, oer_cldPredThresh, oer_tovutil
43 public :: oer_errThreshAllsky, oer_useStateDepSigmaObs
44 ! TOVS OBS ERRORS
45 real(8) :: toverrst(tvs_maxChannelNumber,tvs_maxNumberOfSensors)
46 real(8) :: cldPredThresh(tvs_maxChannelNumber,tvs_maxNumberOfSensors,2)
47 real(8) :: errThreshAllsky(tvs_maxChannelNumber,tvs_maxNumberOfSensors,2)
48 integer :: tovutil(tvs_maxChannelNumber,tvs_maxNumberOfSensors)
49 logical :: useStateDepSigmaObs(tvs_maxChannelNumber,tvs_maxNumberOfSensors)
50 real(8) :: oer_toverrst(tvs_maxChannelNumber,tvs_maxNumberOfSensors)
51 real(8) :: oer_cldPredThresh(tvs_maxChannelNumber,tvs_maxNumberOfSensors,2)
52 real(8) :: oer_errThreshAllsky(tvs_maxChannelNumber,tvs_maxNumberOfSensors,2)
53 real(8) :: clearCldPredThresh(tvs_maxChannelNumber,tvs_maxNumberOfSensors)
54 real(8) :: inflateErrAllskyCoeff(tvs_maxNumberOfSensors)
55 integer :: oer_tovutil(tvs_maxChannelNumber,tvs_maxNumberOfSensors)
56 logical :: oer_useStateDepSigmaObs(tvs_maxChannelNumber,tvs_maxNumberOfSensors)
57
58 ! SST data
59 type SSTdataParamsType
60 character(len=20) :: dataType = '' ! type of data: insitu, satellite, pseudo
61 character(len=20) :: instrument = '' ! instrument: drifts, bouys, ships, AVHRR, VIIRS, AMSR2
62 character(len=20) :: sensor = '' ! sensor of satellite data: NOAA19, NOAA20,...
63 character(len=20) :: sensorType = '' ! type of satellite data sensors: infrared, microwave,..
64 integer :: codeType = MPC_missingValue_INT ! data codtype
65 real(8) :: dayError = MPC_missingValue_R8 ! data error for daytime
66 real(8) :: nightError = MPC_missingValue_R8 ! data error for nighttime
67 end type SSTdataParamsType
68 integer, parameter :: maxNumberSSTDatasets = 15
69
70 ! SST namelist variables
71 integer :: numberSSTDatasets = MPC_missingValue_INT ! MUST NOT BE INCLUDED IN NAMELIST!
72 type(SSTdataParamsType) :: SSTdataParams(maxNumberSSTDatasets) ! list of SSTdataParamsType defining SST obs errors
73
74 ! CONVENTIONAL OBS ERRORS
75 real(8) :: xstd_ua_ai_sw(20,11)
76 real(8) :: xstd_sf(9,6)
77 real(8) :: xstd_pr(2)
78 real(8) :: xstd_sc(1)
79 real(8) :: LVLVALUE(9), HGT_ERR(200,9)
80
81 ! CONSTITUENT/CHEMISTRY OBS ERROR STD DEV
82 ! Associated below to routines beginning with the prefix 'chm'
83 type :: struct_chm_std
84 !
85 ! Structure containing information retrieved from auxiliary obs data file
86 ! holding observation std dev information for constituent obs
87 !
88 ! Variable Description
89 ! -------- -----------
90 ! n_stnid Number of sub-families (identified via STNIDs)
91 ! stnids Sub-families (STNIDs; * are wild cards)
92 ! element BUFR element in data block
93 ! source 0: Set entirely from the auxiliary file being read. No
94 ! initial values read from observation files
95 ! 1: Initial values in observation files
96 ! (may be adjusted after input)
97 ! 2: Initial values in observation files for variable number
98 ! of vertical levels (for error std deviations only)
99 ! std_type Index of setup approach (used in combination with source)
100 ! For source value 0 or 1,
101 ! 0: std1 or observation file values (sigma)
102 ! 1: max(std3,std2*ZVAL) if source=0
103 ! max(std3,std2*sigma) otherwise
104 ! 2: sqrt(std3**2+(std2*ZVAL)**2)) if source=0
105 ! sqrt(std3**2+(std2*sigma)**2)) otherwise
106 ! 3: min(std3,max(std2,std1_chm*ZVAL)) if source=0
107 ! min(std3,max(std2,std1_chm*sigma)) otherwise
108 ! 4: sqrt(std2**2+(std1*ZVAL)**2)) if source=0
109 ! sqrt(std2**2+(std1*sigma)**2)) otherwise
110 ! ibegin Position index of start of data for given
111 ! sub-family in the arrays std1,levels,lat
112 ! n_lvl Number of vertical levels (max number when source=2)
113 ! levels Vertical levels (in coordinate of sub-family data)
114 ! n_lat Number of latitudes
115 ! lat Latitudes (degrees; ordered in increasing size)
116 ! std1 See std_type for usage (dependent on vertical level)
117 ! std2 See std_type for usage (independent of vertical level)
118 ! std3 See std_type for usage (independent of vertical level)
119
120 integer :: n_stnid
121 character(len=12), allocatable :: stnids(:)
122 integer, allocatable :: element(:),std_type(:),n_lat(:)
123 integer, allocatable :: source(:),ibegin(:),n_lvl(:)
124 real(8), allocatable :: std1(:),std2(:),std3(:)
125 real(8), allocatable :: levels(:),lat(:)
126
127 ! Array to hold std dev's read from auxiliary obs data/info file
128 type(struct_oss_obsdata), allocatable :: obsStdDev(:)
129
130 end type struct_chm_std
131
132 type(struct_chm_std) :: chm_std
133
134 ! Sea Ice Concentration obs-error standard deviation
135 real(8) :: xstd_sic(9)
136 ! tiepoint standard deviation for ASCAT backscatter anisotropy
137 integer, parameter :: ncells = 21
138 real(8) :: ascatAnisSigmaOpenWater(ncells,12), ascatAnisSigmaIce(ncells,12)
139 real :: oer_ascatAnisOpenWater(ncells,12), oer_ascatAnisIce(ncells,12)
140
141 ! Hydrology
142 real(8) :: xstd_hydro(1)
143
144 integer :: n_sat_type, n_categorie
145 integer :: tbl_m(200), tbl_h(200), tbl_t(200), tbl_g(200)
146 integer :: surfaceObsTypeNumber
147
148 character(len=9) :: SAT_AMV(200,10), SAT_LIST(200), MET_LIST(200)
149 character(len=9) :: HTM_LIST(200), TMG_LIST(200), NSW_LIST(200)
150
151 logical :: useTovsUtil
152 character(len=48) :: obserrorMode
153
154 ! Namelist variables:
155 logical :: new_oer_sw ! use the 'new' method to compute errors for AMV observations
156 logical :: obsfile_oer_sw ! choose to read errors for AMV from the obs files
157 logical :: visAndGustAdded ! choose to read visibility and gust errors in addition to other conv variables
158 logical :: mwAllskyTtInflateByOmp ! choose to inflate all sky TT radiance errors by an amount related to O-P
159 logical :: mwAllskyTtInflateByClwDiff ! choose to inflate all sky TT radiance errors by an amount related to cloud O-P
160 logical :: mwAllskyHuInflateByOmp ! choose to inflate all sky HU radiance errors by an amount related to O-P
161 logical :: mwAllskyHuInflateBySiDiff ! choose to inflate all sky HU radiance errors by an amount related to cloud O-P
162 real(8) :: amsuaClearCldPredThresh(5) ! cloud threshold for considering obs as clear sky
163 real(8) :: amsuaInflateErrAllskyCoeff ! state dependent obs error inflation factor
164 logical :: readOldSymmetricObsErrFile ! choose to read 'old style' obs error file, when only AMSU-A was all sky
165
166contains
167
168 !--------------------------------------------------------------------------
169 ! oer_setInterchanCorr
170 !--------------------------------------------------------------------------
171 subroutine oer_setInterchanCorr()
172 !
173 !:Purpose: Setup of interchannel observation errors correlations
174 !
175 use rmatrix_mod
176 IMPLICIT NONE
177
178 ! Locals:
179 INTEGER :: ISENS
180
181 if (tvs_nsensors == 0) then
182 write(*,*) 'oer_setInterchanCorr: tvs_nsensors is zero, skipping setup'
183 return
184 end if
185
186! Initialization of the correlation matrices
187 call rmat_init(tvs_nsensors,tvs_nobtov)
188 if (rmat_lnondiagr) then
189 do isens = 1, tvs_nsensors
190 if (tvs_isReallyPresent(isens)) call rmat_readCMatrix(tvs_listSensors(:,isens), isens, tvs_ichan(1:tvs_nchan(isens),isens))
191 end do
192 end if
193
194 END SUBROUTINE oer_setInterchanCorr
195
196 !--------------------------------------------------------------------------
197 ! oer_setObsErrors
198 !--------------------------------------------------------------------------
199 subroutine oer_setObsErrors(obsSpaceData, obserrorMode_in, useTovsUtil_opt)
200 !
201 !:Purpose: read and set observation errors (from former sucovo subroutine).
202 !
203 type(struct_obs) :: obsSpaceData
204 character(len=*), intent(in) :: obserrorMode_in
205 logical, optional :: useTovsUtil_opt
206
207 namelist /namoer/ new_oer_sw, obsfile_oer_sw, visAndGustAdded
208 namelist /namoer/ mwAllskyTtInflateByOmp, mwAllskyTtInflateByClwDiff
209 namelist /namoer/ mwAllskyHuInflateByOmp, mwAllskyHuInflateBySiDiff
210 namelist /namoer/ amsuaClearCldPredThresh
211 namelist /namoer/ amsuaInflateErrAllskyCoeff
212 namelist /namoer/ readOldSymmetricObsErrFile
213 integer :: fnom, fclos, nulnam, ierr
214
215 !
216 !- 1. Setup Mode
217 !
218 obserrorMode = obserrorMode_in
219
220 ! Additional key to allow the use of 'util' column in stats_tovs file
221 if (present(useTovsUtil_opt)) then
222 useTovsUtil = useTovsUtil_opt
223 else
224 useTovsUtil = .false.
225 end if
226
227 ! read namelist namoer
228 new_oer_sw = .false.
229 obsfile_oer_sw = .false.
230 visAndGustAdded = .false.
231 mwAllskyTtInflateByOmp = .false.
232 mwAllskyTtInflateByClwDiff = .false.
233 mwAllskyHuInflateByOmp = .false.
234 mwAllskyHuInflateBySiDiff = .false.
235 amsuaClearCldPredThresh(:) = 0.03D0
236 amsuaClearCldPredThresh(1) = 0.05D0
237 amsuaClearCldPredThresh(4) = 0.02D0
238 amsuaInflateErrAllskyCoeff = 13.0D0
239 readOldSymmetricObsErrFile = .true.
240
241 if (utl_isNamelistPresent('namoer','./flnml')) then
242 nulnam = 0
243 ierr = fnom(nulnam,'./flnml','FTN+SEQ+R/O',0)
244 read (nulnam, nml = NAMOER, iostat = ierr)
245 if (ierr /= 0) call utl_abort('oer_setObsErrors: Error reading namelist')
246 if (mmpi_myid == 0) write(*,nml=namoer)
247 ierr = fclos(nulnam)
248 else
249 write(*,*)
250 write(*,*) 'oer_setObsErrors: namoer is missing in the namelist. The default value will be taken.'
251 end if
252
253 !
254 !- 2. Read in the observation std dev errors
255 !
256
257 !- 2.1 Radiance data
258 if (obs_famexist(obsSpaceData,'TO')) then
259 call oer_readObsErrorsTOVS
260 else
261 write(*,*) "oer_setObsErrors: No brightness temperature observations found."
262 end if
263
264 !- 2.2 Conventional data
265 if (obs_famExist(obsSpaceData, 'UA') .or. obs_famExist(obsSpaceData, 'AI') .or. obs_famExist(obsSpaceData, 'SW') .or. &
266 obs_famExist(obsSpaceData, 'SF') .or. obs_famExist(obsSpaceData, 'GP') .or. obs_famExist(obsSpaceData, 'SC') .or. &
267 obs_famExist(obsSpaceData, 'PR')) then
268
269 call oer_readObsErrorsCONV()
270
271 else
272
273 write(*,*) "oer_setObsErrors: No conventional weather observations found."
274
275 end if
276
277 !- 2.3 Constituent data
278 if (obs_famexist(obsSpaceData,'CH')) then
279 call chm_read_obs_err_stddev
280 else
281 write(*,*) "oer_setObsErrors: No CH observations found."
282 end if
283
284 !- 2.4 Sea ice concentration
285 if (obs_famexist(obsSpaceData,'GL')) then
286 call oer_readObsErrorsIce
287 else
288 write(*,*) "oer_setObsErrors: No GL observations found."
289 end if
290
291 !- 2.5 SST
292 if (obs_famexist(obsSpaceData,'TM')) then
293 call oer_readObsErrorsSST
294 else
295 write(*,*) "oer_setObsErrors: No TM observations found."
296 end if
297
298 !- 2.6 Hydrology
299 if (obs_famexist(obsSpaceData,'HY')) then
300 call oer_readObsErrorsHydro
301 else
302 write(*,*) "oer_setObsErrors: No HY observations found."
303 end if
304
305 !
306 !- 3. Set obs error information in obsSpaceData object
307 !
308 call oer_fillObsErrors(obsSpaceData)
309
310 !
311 !- 4. Deallocate temporary storage
312 !
313 if (obs_famExist(obsSpaceData,'CH')) call chm_dealloc_obs_err_stddev
314
315 end subroutine oer_setObsErrors
316
317 !--------------------------------------------------------------------------
318 ! oer_readObsErrorsTOVS
319 !--------------------------------------------------------------------------
320 subroutine oer_readObsErrorsTOVS
321 !
322 !:Purpose: Read the observation error statistics and
323 ! utilization flag for TOVS processing. This information
324 ! resides on an ASCII file and is read using a free format.
325 !
326 implicit none
327
328 ! Locals:
329 integer, parameter :: bgckColumnIndex = 1
330 integer, parameter :: analysisColumnIndex = 2
331 integer,external :: FNOM, FCLOS
332 integer :: IER, ILUTOV, ILUTOV2, JI, obsErrorColumnIndex, JL, JM
333 integer :: INUMSAT, INUMSAT2, ISAT, IPLF
334 integer :: IPLATFORM(tvs_maxNumberOfSensors), ISATID(tvs_maxNumberOfSensors)
335 integer :: IINSTRUMENT(tvs_maxNumberOfSensors), NUMCHN(tvs_maxNumberOfSensors)
336 integer :: NUMCHNIN(tvs_maxNumberOfSensors)
337 integer :: IPLATFORM2, ISATID2, IINSTRUMENT2, NUMCHNIN2
338 integer :: IUTILST(tvs_maxChannelNumber,tvs_maxNumberOfSensors)
339 integer :: ICHN(tvs_maxChannelNumber,tvs_maxNumberOfSensors)
340 integer :: ICHNIN(tvs_maxChannelNumber,tvs_maxNumberOfSensors)
341 integer :: ICHNIN2(tvs_maxChannelNumber)
342 real(8) :: TOVERRIN(tvs_maxChannelNumber,2,tvs_maxNumberOfSensors)
343 real(8) :: cldPredThreshInput(tvs_maxChannelNumber,tvs_maxNumberOfSensors,2)
344 real(8) :: errThreshAllskyInput(tvs_maxChannelNumber,tvs_maxNumberOfSensors,2)
345 real(8) :: tovsObsInflation(tvs_maxChannelNumber,tvs_maxNumberOfSensors)
346 real(8) :: clearCldPredThreshInput(tvs_maxChannelNumber,tvs_maxNumberOfSensors)
347 real(8) :: inflateErrAllskyCoeffInput(tvs_maxNumberOfSensors)
348 integer :: useStateDepSigmaObsInput(tvs_maxChannelNumber,tvs_maxNumberOfSensors)
349 integer :: amsuaChannelOffset, amsuaChannelNum
350 character (len=132) :: CLDUM,CPLATF,CINSTR
351
352 write(*,*) 'oer_readObsErrorsTOVS: reading observation error statistics required for TOVS processing'
353
354 !
355 ! 1. Initialize
356 ! ----------
357 !
358 TOVERRST(:,:) = 0.0D0
359 TOVERRIN(:,:,:) = 0.0D0
360 cldPredThresh(:,:,:) = 0.0d0
361 cldPredThreshInput(:,:,:) = 0.0d0
362 errThreshAllsky(:,:,:) = 0.0d0
363 errThreshAllskyInput(:,:,:) = 0.0d0
364 tovsObsInflation(:,:) = 0.0d0
365 IUTILST(:,:) = 0
366 useStateDepSigmaObs(:,:) = .false.
367 useStateDepSigmaObsInput(:,:) = 0
368 clearCldPredThresh(:,:) = 0.0d0
369 clearCldPredThreshInput(:,:) = 0.0d0
370 inflateErrAllskyCoeff(:) = 0.0d0
371 inflateErrAllskyCoeffInput(:) = 0.0d0
372
373 IPLATFORM(:) = 0
374 NUMCHN(:) = 0
375 NUMCHNIN(:) = 0
376 ICHN(:,:) = 0
377 ICHNIN(:,:) = 0
378 ICHNIN2(:) = 0
379
380 if (tvs_nobtov == 0) return
381
382 !
383 ! 2. Open the file
384 ! -------------
385 !
386 ilutov = 0
387 IER = FNOM(ILUTOV,'stats_tovs','SEQ+FMT',0)
388 if (IER < 0) call utl_abort ('oer_readObsErrorsTOVS: Problem opening file stats_tovs')
389
390 !
391 ! 3. Read number of satellites
392 ! -------------------------
393 !
394 read(ILUTOV,*)
395 read(ILUTOV,*) INUMSAT
396 read(ILUTOV,*)
397
398 if (inumsat > tvs_maxNumberOfSensors) then
399 write(*,'(A)') ' Number of sensors in stats_tovs file is greater than maximum allowed (tvs_maxNumberOfSensors)'
400 call utl_abort('oer_readObsErrorsTOVS')
401 end if
402
403 !
404 ! 4. Read the satellite identification, the number of channels,
405 ! the observation errors and the utilization flags
406 ! ----------------------------------------------------------
407 !
408 write(*,*) 'oer_readObsErrorsTOVS: Reading stats_tovs file.'
409
410 DO JL = 1, INUMSAT
411
412 read(ILUTOV,*)
413 read(ILUTOV,'(A)') CLDUM
414 write(*,'(A)') CLDUM
415 CINSTR=CLDUM
416 call split(CINSTR," ",CPLATF)
417 write(*,*) "CINSTR: ",CINSTR
418 write(*,*) "CPLATF: ",CPLATF
419 read(ILUTOV,*)
420 read(ILUTOV,*) ISATID(JL), NUMCHNIN(JL)
421
422 do JI = 1, 3
423 read(ILUTOV,*)
424 end do
425
426 if (CPLATF == 'FY-3C') THEN
427 CPLATF = 'FY3-3'
428 CINSTR = 'MWHS2'
429 end if
430
431 IPLATFORM(JL) = tvs_getPlatformId(CPLATF)
432
433 if (IPLATFORM(JL) == -1) call utl_abort ('oer_readObsErrorsTOVS: Unknown platform!')
434
435 IINSTRUMENT(JL) = tvs_getInstrumentId(CINSTR)
436
437 if (IINSTRUMENT(JL) == -1) call utl_abort ('oer_readObsErrorsTOVS: Unknown instrument!')
438
439 do JI = 1, NUMCHNIN(JL)
440 read(ILUTOV,*) ICHNIN(JI,JL), TOVERRIN(ICHNIN(JI,JL),1,JL), TOVERRIN(ICHNIN(JI,JL),2,JL), IUTILST(ICHNIN(JI,JL),JL), tovsObsInflation(ICHNIN(JI,JL),JL)
441 end do
442 read(ILUTOV,*)
443
444 end do
445
446 ! read in the parameters to define the user-defined symmetric TOVS errors
447 if (tvs_mwAllskyAssim) then
448 ilutov2 = 10
449 IER = FNOM(ILUTOV2,'stats_tovs_symmetricObsErr','SEQ+FMT',0)
450 if (IER < 0) call utl_abort ('oer_readObsErrorsTOVS: Problem opening symmetricObsErr file.')
451
452 read(ILUTOV2,*)
453 read(ILUTOV2,*) INUMSAT2
454 read(ILUTOV2,*)
455
456 write(*,*) 'oer_readObsErrorsTOVS: Reading symmetricObsErr file.'
457
458 DO JL = 1, INUMSAT2
459
460 read(ILUTOV2,*)
461 read(ILUTOV2,'(A)') CLDUM
462 write(*,'(A)') CLDUM
463 CINSTR = CLDUM
464 call split(CINSTR," ",CPLATF)
465 write(*,*) "CINSTR: ",CINSTR
466 write(*,*) "CPLATF: ",CPLATF
467 read(ILUTOV2,*)
468
469 ! If reading the old style stats_tovs_symmetricObsErr, the the all-sky parameters are available only for AMSUA.
470 if (readOldSymmetricObsErrFile) then
471 read(ILUTOV2,*) ISATID2, NUMCHNIN2
472 if (CINSTR == "AMSUA") inflateErrAllskyCoeffInput(JL) = amsuaInflateErrAllskyCoeff
473 else
474 read(ILUTOV2,*) ISATID2, NUMCHNIN2, inflateErrAllskyCoeffInput(JL)
475 end if
476
477 if (ISATID2 /= ISATID(JL) .or. NUMCHNIN2 /= NUMCHNIN(JL)) &
478 call utl_abort ('oer_readObsErrorsTOVS: problem with ISATID2, NUMCHNIN2 in symmetricObsErr')
479
480 do JI = 1, 3
481 read(ILUTOV2,*)
482 end do
483
484 IPLATFORM2 = tvs_getPlatformId(CPLATF)
485 IINSTRUMENT2 = tvs_getInstrumentId(CINSTR)
486 if (IPLATFORM2 /= IPLATFORM(JL) .or. IINSTRUMENT2 /= IINSTRUMENT(JL)) &
487 call utl_abort ('oer_readObsErrorsTOVS: problem with IPLATFORM2, IINSTRUMENT2 in symmetricObsErr')
488
489 do JI = 1, NUMCHNIN2
490 ! If reading the old style stats_tovs_symmetricObsErr, then the all-sky parameters are available only for AMSUA.
491 if (readOldSymmetricObsErrFile) then
492 read(ILUTOV2,*) ICHNIN2(JI), &
493 cldPredThreshInput(ICHNIN2(JI),JL,1), cldPredThreshInput(ICHNIN2(JI),JL,2), &
494 errThreshAllskyInput(ICHNIN2(JI),JL,1), errThreshAllskyInput(ICHNIN2(JI),JL,2), &
495 useStateDepSigmaObsInput(ICHNIN2(JI),JL)
496 if (CINSTR == "AMSUA") then
497 amsuaChannelOffset = 27
498 amsuaChannelNum = ICHNIN2(JI) - amsuaChannelOffset
499 if (amsuaChannelNum >= 1 .and. amsuaChannelNum <= 5) then
500 clearCldPredThreshInput(ICHNIN2(JI),JL) = amsuaClearCldPredThresh(amsuaChannelNum)
501 end if
502 end if
503 else
504 read(ILUTOV2,*) ICHNIN2(JI), &
505 cldPredThreshInput(ICHNIN2(JI),JL,1), cldPredThreshInput(ICHNIN2(JI),JL,2), &
506 errThreshAllskyInput(ICHNIN2(JI),JL,1), errThreshAllskyInput(ICHNIN2(JI),JL,2), &
507 clearCldPredThreshInput(ICHNIN2(JI),JL), &
508 useStateDepSigmaObsInput(ICHNIN2(JI),JL)
509 end if
510
511 if (ICHNIN2(JI) /= ICHNIN(JI,JL)) &
512 call utl_abort ('oer_readObsErrorsTOVS: problem with ICHNIN2 in symmetricObsErr')
513
514 end do
515 read(ILUTOV2,*)
516
517 end do
518
519 IER = FCLOS(ILUTOV2)
520 if (IER /= 0) call utl_abort ('oer_readObsErrorsTOVS')
521
522 end if
523
524 !
525 ! Select input error to use: if ANAL mode, use ERRANAL (obsErrorColumnIndex=2);
526 ! otherwise use ERRBGCK (obsErrorColumnIndex=1)
527 !
528 if (trim(obserrorMode) == 'analysis' .or. trim(obserrorMode) == 'FSO') THEN
529 obsErrorColumnIndex = analysisColumnIndex
530 ELSE
531 obsErrorColumnIndex = bgckColumnIndex
532 end if
533
534 !
535 ! Fill the observation error array TOVERRST
536 !
537 write(*,*) 'oer_readObsErrorsTOVS: Fill error array TOVERRST.'
538 do JM = 1, INUMSAT
539 do JL = 1, tvs_nsensors
540 if (tvs_platforms (JL) == IPLATFORM(JM) .AND. tvs_satellites(JL) == ISATID(JM)) THEN
541 if (tvs_instruments (JL) == IINSTRUMENT(JM)) THEN
542 NUMCHN(JL)=NUMCHNIN(JM)
543 do JI = 1, tvs_maxChannelNumber
544 TOVERRST(JI,JL) = TOVERRIN(JI,obsErrorColumnIndex,JM)
545 ICHN(JI,JL) = ICHNIN(JI,JM)
546
547 if (tvs_mwAllskyAssim) then
548 cldPredThresh(JI,JL,:) = cldPredThreshInput(JI,JM,:)
549 errThreshAllsky(JI,JL,:) = errThreshAllskyInput(JI,JM,:)
550 clearCldPredThresh(JI,JL) = clearCldPredThreshInput(JI,JM)
551 useStateDepSigmaObs(JI,JL) = (useStateDepSigmaObsInput(JI,JM) == 1)
552
553 if (JI == 1) inflateErrAllskyCoeff(JL) = inflateErrAllskyCoeffInput(JM)
554
555 ! inflate the errThreshAllsky in analysis mode
556 if (obsErrorColumnIndex == analysisColumnIndex) then
557 errThreshAllsky(JI,JL,1) = errThreshAllsky(JI,JL,1) * tovsObsInflation(JI,JM)
558 errThreshAllsky(JI,JL,2) = errThreshAllsky(JI,JL,2) * tovsObsInflation(JI,JM)
559 end if
560 end if
561
562 end do
563 if (trim(obserrorMode) == 'bgck' .or. useTovsUtil) THEN
564 do JI = 1, tvs_maxChannelNumber
565 ! channel Selection using array IUTILST(chan,sat):
566 ! 0 (blacklisted)
567 ! 1 (assmilate)
568 ! 2 (assimilate over open water only)
569 tovutil(JI,JL) = IUTILST(JI,JM)
570 end do
571 end if
572 end if
573 end if
574 end do
575 end do
576
577 !
578 ! Check that oberservation error statistics have been defined for
579 ! all the satellites specified in the namelist.
580 !
581 do JL = 1, tvs_nsensors
582 IPLF = utl_findloc(IPLATFORM(:) , tvs_platforms(JL))
583 ISAT = utl_findloc(ISATID(:), tvs_satellites(JL))
584 if (IPLF == 0 .OR. ISAT == 0) THEN
585 write(*,'(A,I3)') 'oer_readObsErrorsTOVS: Observation errors not defined for sensor #', JL
586 call utl_abort ('oer_readObsErrorsTOVS')
587 end if
588 if (NUMCHN(JL) == 0) THEN
589 write(*,'(A,I3)') 'oer_readObsErrorsTOVS: Problem setting errors for sensor #', JL
590 call utl_abort ('oer_readObsErrorsTOVS')
591 end if
592 end do
593
594 !
595 ! 5. Print out observation errors for each sensor
596 ! --------------------------------------------
597 !
598 if (mmpi_myid == 0) THEN
599 write(*,*) 'Radiance observation errors read from file'
600 write(*,*) '------------------------------------------'
601 do JL = 1, tvs_nsensors
602 write(*,'(A,I2,4(A))') 'SENSOR #', JL, ', Platform: ', tvs_satelliteName(JL), &
603 ', Instrument: ',tvs_instrumentName(JL)
604 if (tvs_mwAllskyAssim .and. any(useStateDepSigmaObs(ICHN(1:NUMCHN(JL),JL),JL))) then
605 write(*,'(A,4(2X,A8),(1X,A9),(2X,A3))') 'Channel','cld1','cld2','sigmaO1','sigmaO2','anlErrInf','use'
606 do JI = 1, NUMCHN(JL)
607 write(*,'(I7,5(2X,F8.4),(2X,L3))') ICHN(JI,JL), &
608 cldPredThresh(ICHN(JI,JL),JL,1), cldPredThresh(ICHN(JI,JL),JL,2), &
609 errThreshAllsky(ICHN(JI,JL),JL,1), errThreshAllsky(ICHN(JI,JL),JL,2), &
610 clearCldPredThresh(ICHN(JI,JL),JL), &
611 useStateDepSigmaObs(ICHN(JI,JL),JL)
612 end do
613 write(*,'(A,(2X,F8.4))') 'inflateErrAllskyCoeff=', &
614 inflateErrAllskyCoeff(JL)
615 else
616 write(*,'(A,2X,A8)') 'Channel','sigmaO'
617 do JI = 1, NUMCHN(JL)
618 write(*,'(I7,1(2X,F8.4))') ICHN(JI,JL),TOVERRST(ICHN(JI,JL),JL)
619 end do
620 end if
621 end do
622 end if
623
624 !
625 ! 6. Close the file
626 ! --------------
627 !
628 IER = FCLOS(ILUTOV)
629 if (IER /= 0) call utl_abort ('oer_readObsErrorsTOVS')
630
631 !
632 ! 7. Fill the publics variables for QC purpose
633 ! --------------
634 oer_toverrst(:,:) = toverrst(:,:)
635 oer_tovutil (:,:) = tovutil(:,:)
636 oer_cldPredThresh(:,:,:) = cldPredThresh(:,:,:)
637 oer_errThreshAllsky(:,:,:) = errThreshAllsky(:,:,:)
638 oer_useStateDepSigmaObs(:,:) = useStateDepSigmaObs(:,:)
639
640 contains
641
642 !--------------------------------------------------------------------------
643 ! compact
644 !--------------------------------------------------------------------------
645 subroutine compact(str)
646 ! Code from Benthien's module: http://www.gbenthien.net/strings/index.html
647 ! Converts multiple spaces and tabs to single spaces; deletes control characters;
648 ! removes initial spaces.
649 implicit none
650
651 ! Arguments:
652 character(len=*), intent(inout) :: str
653
654 ! Locals:
655 character(len=1):: ch
656 character(len=len_trim(str)):: outstr
657 integer isp,k,lenstr,i,ich
658
659 str=adjustl(str)
660 lenstr=len_trim(str)
661 outstr=' '
662 isp=0
663 k=0
664
665 do i=1,lenstr
666 ch=str(i:i)
667 ich=iachar(ch)
668
669 select case(ich)
670 case(9,32) ! space or tab character
671 if(isp==0) then
672 k=k+1
673 outstr(k:k)=' '
674 end if
675 isp=1
676 case(33:) ! not a space, quote, or control character
677 k=k+1
678 outstr(k:k)=ch
679 isp=0
680 end select
681
682 end do
683
684 str=adjustl(outstr)
685
686 end subroutine compact
687
688 !--------------------------------------------------------------------------
689 ! split
690 !--------------------------------------------------------------------------
691 subroutine split(str, delims, before)
692 !
693 !:Comment:
694 ! Code extracted from Benthien's module: http://www.gbenthien.net/strings/index.html
695 ! Routine finds the first instance of a character from 'delims' in the
696 ! the string 'str'. The characters before the found delimiter are
697 ! output in 'before'. The characters after the found delimiter are
698 ! output in 'str'.
699 !
700 implicit none
701
702 ! Arguments:
703 character(len=*), intent(inout) :: str
704 character(len=*), intent(in) :: delims
705 character(len=*), intent(out) :: before
706
707 ! Locals:
708 character :: ch,cha
709 integer lenstr,i,k,ipos,iposa
710 str=adjustl(str)
711 call compact(str)
712 lenstr=len_trim(str)
713
714 if(lenstr == 0) return ! string str is empty
715 k=0
716 before=' '
717 do i=1,lenstr
718 ch=str(i:i)
719
720 ipos=index(delims,ch)
721
722 if(ipos == 0) then ! character is not a delimiter
723 k=k+1
724 before(k:k)=ch
725 cycle
726 end if
727 if(ch /= ' ') then ! character is a delimiter that is not a space
728 str=str(i+1:)
729 exit
730 end if
731
732 cha=str(i+1:i+1) ! character is a space delimiter
733 iposa=index(delims,cha)
734 if(iposa > 0) then ! next character is a delimiter
735 str=str(i+2:)
736 exit
737 else
738 str=str(i+1:)
739 exit
740 end if
741 end do
742 if(i >= lenstr) str=''
743
744 str=adjustl(str) ! remove initial spaces
745
746 end subroutine split
747
748 end subroutine oer_readObsErrorsTOVS
749
750 !--------------------------------------------------------------------------
751 ! oer_readObsErrorsCONV
752 !--------------------------------------------------------------------------
753 subroutine oer_readObsErrorsCONV()
754 !
755 !:Purpose: read observation errors (modification of former readcovo subroutine) of conventional data
756 !
757 implicit none
758
759 ! Locals:
760 integer :: fnom, fclos, ierr, jlev, jelm, jcat, icodtyp, nulstat
761 logical :: LnewExists
762 character (len=128) :: ligne
763
764 if (visAndGustAdded) then
765 surfaceObsTypeNumber = 6
766 else
767 surfaceObsTypeNumber = 4
768 end if
769
770 ! CHECK THE EXISTENCE OF THE NEW FILE WITH STATISTICS
771 inquire(file = 'obserr', exist = LnewExists)
772 if (LnewExists) then
773 write(*,*) '--------------------------------------------------------'
774 write(*,*) 'oer_readObsErrorsCONV: reads observation errors in obserr'
775 write(*,*) '--------------------------------------------------------'
776 else
777 call utl_abort('oer_readObsErrorsCONV: NO OBSERVATION STAT FILE FOUND!!')
778 end if
779
780 ! Read observation errors from file obserr for conventional data
781 nulstat=0
782 ierr=fnom(nulstat, 'obserr', 'SEQ', 0)
783 if (ierr == 0) then
784 write(*,*) 'oer_readObsErrorsCONV: File = ./obserr'
785 write(*,*) ' opened as unit file ',nulstat
786 open(unit=nulstat, file='obserr', status='OLD')
787 else
788 call utl_abort('oer_readObsErrorsCONV: COULD NOT OPEN FILE obserr!!!')
789 end if
790
791 write(*, '(A)') ' '
792
793 do jlev = 1,3
794 read(nulstat, '(A)') ligne
795 write(*, '(A)') ligne
796 end do
797
798 do jlev = 1, 19
799 read(nulstat, *) (xstd_ua_ai_sw(jlev,jelm), jelm=1,11)
800 write(*, '(f6.0,10f6.1)') (xstd_ua_ai_sw(jlev,jelm), jelm=1,11)
801 end do
802
803 do jlev = 1,5
804 read(nulstat, '(A)') ligne
805 write(*, '(A)') ligne
806 end do
807
808 read(nulstat, *) xstd_pr(1),xstd_pr(2)
809 write(*, '(2f6.1)') xstd_pr(1),xstd_pr(2)
810
811 do jlev = 1,5
812 read(nulstat, '(A)') ligne
813 write(*, '(A)') ligne
814 end do
815
816 read(nulstat, *) xstd_sc(1)
817 write(*, '(f8.3)') xstd_sc(1)
818
819 read(nulstat, '(A)') ligne
820 write(*, '(A)') ligne
821
822 do icodtyp = 1,9
823 do jlev = 1,4
824 read(nulstat, '(A)') ligne
825 write(*, '(A)') ligne
826 end do
827 read(nulstat, *) (xstd_sf(icodtyp,jelm), jelm=1,surfaceObsTypeNumber)
828 write(*, '(f6.2,2f6.1,f8.3)') (xstd_sf(icodtyp,jelm), jelm=1,surfaceObsTypeNumber)
829 end do
830
831 !
832 ! Read height assignment errors
833 ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!
834 !
835
836 if(new_oer_sw) then
837
838 do jlev = 1,5
839 read(nulstat, '(A)') ligne
840 write(*, '(A)') ligne
841 end do
842
843 read(nulstat, *) (LVLVALUE(jelm), jelm=1,9)
844 write(*, '(9f6.1)') (LVLVALUE(jelm), jelm=1,9)
845
846 do jlev = 1,2
847 read(nulstat, '(A)') ligne
848 write(*, '(A)') ligne
849 end do
850
851 read(nulstat, '(i10)') n_sat_type
852 write(*,'(i10)') n_sat_type
853 do jlev = 1,n_sat_type
854 read(nulstat, '(10A10)') (SAT_AMV(jlev,jelm), jelm=1,10)
855 write(*, '(10A10)') (SAT_AMV(jlev,jelm), jelm=1,10)
856 end do
857
858 read(nulstat, '(A)') ligne
859 write(*, '(A)') ligne
860
861 tbl_m(:) = 0
862 tbl_h(:) = 0
863 tbl_t(:) = 0
864 tbl_g(:) = 0
865 read(nulstat, '(i10)') n_categorie
866 write(*,'(i10)') n_categorie
867 do jcat = 1,n_categorie
868 read(nulstat, *) SAT_LIST(jcat),MET_LIST(jcat),HTM_LIST(jcat),TMG_LIST(jcat),NSW_LIST(jcat),(HGT_ERR(jcat,jelm), jelm=1,9)
869 write(*, '(5A10,9f6.1)') SAT_LIST(jcat),MET_LIST(jcat),HTM_LIST(jcat),TMG_LIST(jcat),NSW_LIST(jcat),(HGT_ERR(jcat,jelm), jelm=1,9)
870 if(trim(MET_LIST(jcat)) == 'ir') tbl_m(jcat) = 1
871 if(trim(MET_LIST(jcat)) == 'vis') tbl_m(jcat) = 2
872 if(trim(MET_LIST(jcat)) == 'wv') tbl_m(jcat) = 3
873 if(trim(HTM_LIST(jcat)) == 'ebbt') tbl_h(jcat) = 1
874 if(trim(HTM_LIST(jcat)) == 'co2') tbl_h(jcat) = 4
875 if(trim(TMG_LIST(jcat)) == 'sea') tbl_t(jcat) = 1
876 if(trim(TMG_LIST(jcat)) == 'land') tbl_t(jcat) = 2
877 if(trim(TMG_LIST(jcat)) == 'ice') tbl_t(jcat) = 3
878 if(trim(NSW_LIST(jcat)) == 'NH') tbl_g(jcat) = 1
879 if(trim(NSW_LIST(jcat)) == 'SH') tbl_g(jcat) = 2
880 end do
881
882 end if
883
884 write(*, '(A)') ' '
885
886 close(unit = nulstat)
887 ierr = fclos(nulstat)
888
889 end subroutine oer_readObsErrorsCONV
890
891 !--------------------------------------------------------------------------
892 ! oer_readObsErrorsIce
893 !--------------------------------------------------------------------------
894 subroutine oer_readObsErrorsIce
895 !
896 !:Purpose: read observation errors for sea ice concentration analysis
897 !
898 implicit none
899
900 ! Locals:
901 external fnom, fclos
902 integer :: fnom, fclos, ierr, jlev, jelm, nulstat
903 logical :: fileExists
904 character(len=128) :: ligne
905 character(len=15), parameter :: fileName = 'sea_ice_obs-err'
906 ! Variables for the ASCAT backscatter anisotropy
907 integer :: jcell_no, icell_no, imonth
908 real :: tiePoint12, tiePoint13, tiePoint23
909
910 ! CHECK THE EXISTENCE OF THE FILE WITH STATISTICS
911 inquire(file = fileName, exist = fileExists)
912 if (fileExists) then
913 write(*,*) '--------------------------------------------------------'
914 write(*,*) 'oer_readObsErrorsICE: reads observation errors in ',fileName
915 write(*,*) '--------------------------------------------------------'
916 else
917 call utl_abort('oer_readObsErrorsICE: NO OBSERVATION STAT FILE FOUND!!')
918 end if
919
920 ! Read observation errors from file
921 nulstat=0
922 ierr=fnom(nulstat, fileName, 'SEQ', 0)
923 if (ierr == 0) then
924 write(*,*) 'oer_readObsErrorsICE: File = ./',fileName
925 write(*,*) ' opened as unit file ',nulstat
926 open(unit=nulstat, file=fileName, status='OLD')
927 else
928 call utl_abort('oer_readObsErrorsICE: COULD NOT OPEN FILE '//fileName//'!!!')
929 end if
930
931 if (mmpi_myid == 0) write(*, '(A)') ' '
932
933 do jelm = 1, 9
934
935 do jlev = 1, 3
936 read(nulstat, '(A)') ligne
937 if (mmpi_myid == 0) write(*, '(A)') ligne
938 end do
939
940 read(nulstat, *) xstd_sic(jelm)
941 if (mmpi_myid == 0) write(*,*) xstd_sic(jelm)
942
943 do jlev = 1, 2
944 read(nulstat, '(A)') ligne
945 if (mmpi_myid == 0) write(*, '(A)') ligne
946 end do
947
948 end do
949
950 ! Read coefficients for ASCAT backscatter anisotropy
951 do jlev = 1, 5
952 read(nulstat, '(A)') ligne
953 if (mmpi_myid == 0) write(*, '(A)') ligne
954 end do
955
956 do imonth = 1, 12
957 do jlev = 1, 3
958 read(nulstat, '(A)') ligne
959 if (mmpi_myid == 0) write(*, '(A)') ligne
960 end do
961 do jcell_no = 1, ncells
962 read(nulstat, *) icell_no, tiePoint12, tiePoint13, tiePoint23, &
963 oer_ascatAnisIce(jcell_no,imonth), oer_ascatAnisOpenWater(jcell_no,imonth), &
964 ascatAnisSigmaIce(jcell_no,imonth), ascatAnisSigmaOpenWater(jcell_no,imonth)
965 if (mmpi_myid == 0) write(*,*) icell_no, tiePoint12, tiePoint13, tiePoint23, &
966 oer_ascatAnisIce(jcell_no,imonth), oer_ascatAnisOpenWater(jcell_no,imonth), &
967 ascatAnisSigmaIce(jcell_no,imonth), ascatAnisSigmaOpenWater(jcell_no,imonth)
968 end do
969 end do
970
971 if (mmpi_myid == 0) write(*, '(A)') ' '
972
973 close(unit = nulstat)
974 ierr = fclos(nulstat)
975
976 end subroutine oer_readObsErrorsIce
977
978 !--------------------------------------------------------------------------
979 ! oer_readObsErrorsSST
980 !--------------------------------------------------------------------------
981 subroutine oer_readObsErrorsSST
982 !
983 !:Purpose: read observation errors for SST data
984 !
985 implicit none
986
987 ! Locals:
988 integer :: fnom, fclos, nulnam, ierr, indexDataset
989 namelist /namSSTObsErrors/ numberSSTDatasets, SSTdataParams
990
991 if (utl_isNamelistPresent('namSSTObsErrors','./flnml')) then
992 nulnam = 0
993 ierr = fnom(nulnam, './flnml', 'FTN+SEQ+R/O', 0)
994 read (nulnam, nml = namSSTObsErrors, iostat = ierr)
995 if (ierr /= 0) call utl_abort('oer_readObsErrorsSST: Error reading namelist')
996 if (mmpi_myid == 0) write(*,nml=namSSTObsErrors)
997 ierr = fclos(nulnam)
998 if (numberSSTDatasets /= MPC_missingValue_INT) then
999 call utl_abort('oer_readObsErrorsSST: check namSSTObsErrors namelist section: numberSSTDatasets should be removed')
1000 end if
1001 numberSSTDatasets = 0
1002 do indexDataset = 1, maxNumberSSTDatasets
1003 if (trim(SSTdataParams(indexDataset)%dataType) == "" .and. &
1004 trim(SSTdataParams(indexDataset)%instrument) == "" .and. &
1005 trim(SSTdataParams(indexDataset)%sensor) == "" .and. &
1006 trim(SSTdataParams(indexDataset)%sensorType) == "" .and. &
1007 SSTdataParams(indexDataset)%codeType == MPC_missingValue_INT .and. &
1008 SSTdataParams(indexDataset)%dayError == MPC_missingValue_R8 .and. &
1009 SSTdataParams(indexDataset)%nightError == MPC_missingValue_R8 ) exit
1010 numberSSTDatasets = numberSSTDatasets + 1
1011 end do
1012 else
1013 call utl_abort('oer_readObsErrorsSST: namSSTObsErrors is missing in the namelist.')
1014 end if
1015
1016 end subroutine oer_readObsErrorsSST
1017
1018 !--------------------------------------------------------------------------
1019 ! oer_readObsErrorsHydro
1020 !--------------------------------------------------------------------------
1021 subroutine oer_readObsErrorsHydro
1022 implicit none
1023
1024 ! Locals:
1025 external :: fnom, fclos
1026 integer :: fnom, fclos, ierr, nulstat
1027 logical :: fileExists
1028 character(len=15), parameter :: fileName = 'obserr_hydro'
1029 character(len=*) , parameter :: myName = 'oer_readObsErrorsHydro'
1030
1031 inquire(file = fileName, exist = fileExists)
1032 if (fileExists) then
1033 write(*,*) '--------------------------------------------------------'
1034 write(*,*) myName//': reads observation errors in ', fileName
1035 write(*,*) '--------------------------------------------------------'
1036 else
1037 call utl_abort(myName//': NO OBSERVATION STAT FILE FOUND!!')
1038 end if
1039
1040 nulstat=0
1041 ierr=fnom(nulstat, fileName, 'SEQ', 0)
1042 if (ierr == 0) then
1043 write(*,*) myName//': File = ./',fileName
1044 write(*,*) myName//' opened as unit file ',nulstat
1045 open(unit=nulstat, file=fileName, status='OLD')
1046 else
1047 call utl_abort(myName//': COULD NOT OPEN FILE '//fileName//'!!!')
1048 end if
1049
1050 read(nulstat, *) xstd_hydro(1)
1051 write(*, '(f8.3)') xstd_hydro(1)
1052
1053 close(unit = nulstat)
1054 ierr = fclos(nulstat)
1055
1056 end subroutine oer_readObsErrorsHydro
1057
1058 !--------------------------------------------------------------------------
1059 ! oer_fillObsErrors
1060 !--------------------------------------------------------------------------
1061 subroutine oer_fillObsErrors(obsSpaceData)
1062 !
1063 !:Purpose: fill observation errors in obsSpaceData
1064 !
1065 implicit none
1066
1067 ! Arguments:
1068 type(struct_obs), intent(inout) :: obsSpaceData
1069
1070 ! Locals:
1071 integer :: jn, JI, bodyIndex, headerIndex, ityp, iass, idata, idatend, codeType
1072 integer :: sensorIndex
1073 integer :: isat, channelNumber, iplatf, instr, iplatform, instrum
1074 integer :: ilev, nlev, idate, itime
1075 integer :: ielem, icodtyp, header_prev, indexDataset, indexSensor
1076 real(8) :: zlat, zlon, zlev, zval, zwb, zwt, obs_err_stddev, solarZenith
1077 real(8) :: obsValue, obsStdDevError, obsPPP, obsOER
1078 real(8) :: cldPredThresh1, cldPredThresh2, cldPredUsed
1079 real(8) :: errThresh1, errThresh2, sigmaObsErrUsed
1080 real(4), parameter :: minRetrievableClwValue_r4 = 0.0
1081 real(4), parameter :: maxRetrievableClwValue_r4 = 3.0
1082 real(4), parameter :: minRetrievableSiValue_r4 = -10.0
1083 real(4), parameter :: maxRetrievableSiValue_r4 = 30.0
1084 logical :: ifirst, surfTypeIsWater, unsupportedCodeType, unsupportedSensor
1085 character(len=2) :: cfam
1086 character(len=12) :: cstnid
1087
1088 write(*,'(10X, "***********************************")')
1089 write(*,'(10X, "oer_fillObsErrors: starting", /)')
1090 write(*,'(10X, "***********************************")')
1091
1092 header_prev=-1
1093
1094 ! SET STANDARD DEVIATION ERRORS FOR EACH DATA FAMILY
1095 HEADER: do headerIndex = 1, obs_numheader(obsSpaceData)
1096
1097 idata = obs_headElem_i(obsSpaceData, OBS_RLN, headerIndex)
1098 idatend = obs_headElem_i(obsSpaceData, OBS_NLV, headerIndex) + idata - 1
1099 cfam = obs_getFamily (obsSpaceData, headerIndex_opt=headerIndex)
1100 zlat = obs_headElem_r(obsSpaceData, OBS_LAT, headerIndex)
1101 zlon = obs_headElem_r(obsSpaceData, OBS_LON, headerIndex)
1102 codeType = obs_headElem_i(obsSpaceData, OBS_ITY, headerIndex)
1103 iplatf = obs_headElem_i(obsSpaceData, OBS_SAT, headerIndex)
1104 instr = obs_headElem_i(obsSpaceData, OBS_INS, headerIndex)
1105 cstnid = obs_elem_c (obsSpaceData, 'STID' , headerIndex)
1106 idate = obs_headElem_i(obsSpaceData, OBS_DAT, headerIndex)
1107 itime = obs_headElem_i(obsSpaceData, OBS_ETM, headerIndex)
1108
1109 surfTypeIsWater = (tvs_ChangedStypValue(obsSpaceData,headerIndex) == surftype_sea)
1110
1111 nlev = idatend - idata + 1
1112
1113 BODY: do bodyIndex = idata, idatend
1114
1115 ityp = obs_bodyElem_i(obsSpaceData, OBS_VNM, bodyIndex)
1116 iass = obs_bodyElem_i(obsSpaceData, OBS_ASS, bodyIndex)
1117 zval = obs_bodyElem_r(obsSpaceData, OBS_VAR, bodyIndex)
1118
1119 if (iass == obs_assimilated) then
1120
1121 !***********************************************************************
1122 ! TOVS DATA
1123 !***********************************************************************
1124
1125 if (cfam == 'TO') then
1126
1127 if (ityp == BUFR_NBT1 .or. &
1128 ityp == BUFR_NBT2 .or. &
1129 ityp == BUFR_NBT3) then
1130
1131 channelNumber = NINT(obs_bodyElem_r(obsSpaceData, OBS_PPP, bodyIndex))
1132 call tvs_mapSat(iplatf, iplatform, isat)
1133 call tvs_mapInstrum(instr, instrum)
1134
1135 do sensorIndex = 1, tvs_nsensors
1136 if (iplatform == tvs_platforms(sensorIndex) .and. &
1137 isat == tvs_satellites(sensorIndex) .and. &
1138 instrum == tvs_instruments(sensorIndex)) then
1139
1140 ! decide whether or not use the state dependent sigmaObsErrUsed for OBS_OER
1141 if (tvs_mwAllskyAssim .and. &
1142 useStateDepSigmaObs(channelNumber,sensorIndex) .and. &
1143 surfTypeIsWater) then
1144
1145 ! set dummy value for OBS_OER in bgck mode
1146 if (trim(obserrorMode) == 'bgck') then
1147 sigmaObsErrUsed = 1.0d0
1148 else
1149 cldPredThresh1 = cldPredThresh(channelNumber,sensorIndex,1)
1150 cldPredThresh2 = cldPredThresh(channelNumber,sensorIndex,2)
1151 errThresh1 = errThreshAllsky(channelNumber,sensorIndex,1)
1152 errThresh2 = errThreshAllsky(channelNumber,sensorIndex,2)
1153
1154 ! Compute cloudPredictor to use
1155 cldPredUsed = computeCloudPredictor(sensorIndex, headerIndex)
1156
1157 ! Use cloud predictor to compute state-dependent obs error
1158 sigmaObsErrUsed = calcStateDepObsErr(cldPredThresh1, cldPredThresh2, &
1159 errThresh1, errThresh2, cldPredUsed)
1160 end if
1161 else
1162 sigmaObsErrUsed = TOVERRST(channelNumber, sensorIndex)
1163 end if
1164 call obs_bodySet_r(obsSpaceData, OBS_OER, bodyIndex, sigmaObsErrUsed)
1165
1166 ! Utilization flag for AIRS,IASI and CrIS channels (bgck mode only)
1167 if (trim(obserrorMode) == 'bgck' .or. useTovsUtil) then
1168 if (tovutil(channelNumber, sensorIndex) == 0) &
1169 call obs_bodySet_i(obsSpaceData, OBS_FLG, bodyIndex, ibset(obs_bodyElem_i(obsSpaceData, OBS_FLG, bodyIndex), 8))
1170 end if
1171 end if
1172 end do
1173
1174 end if
1175
1176 !***********************************************************************
1177 ! RADIOSONDE DATA
1178 !***********************************************************************
1179
1180 else if (cfam == 'UA') then
1181
1182 zlev = obs_bodyElem_r(obsSpaceData, OBS_PPP, bodyIndex)
1183
1184 if ((ityp == BUFR_NEUS) .or. (ityp == BUFR_NEVS)) then
1185 call obs_bodySet_r(obsSpaceData, OBS_OER, bodyIndex, xstd_sf(1,4))
1186 else if (ityp == BUFR_NETS) then
1187 call obs_bodySet_r(obsSpaceData, OBS_OER, bodyIndex, xstd_sf(1,2))
1188 else if (ityp == BUFR_NESS) then
1189 call obs_bodySet_r(obsSpaceData, OBS_OER, bodyIndex, xstd_sf(1,3))
1190 else if (ityp == BUFR_NEPS) then
1191 call obs_bodySet_r(obsSpaceData, OBS_OER, bodyIndex, xstd_sf(1,1))
1192 else if (ityp == BUFR_NEPN) then
1193 call obs_bodySet_r(obsSpaceData, OBS_OER, bodyIndex, xstd_sf(1,1))
1194 else
1195
1196 if ((ityp == BUFR_NEUU) .or. (ityp == BUFR_NEVV)) then
1197 ielem = 4
1198 else if (ityp == BUFR_NETT) then
1199 ielem = 2
1200 else if (ityp == BUFR_NEES) then
1201 ielem = 3
1202 else if (ityp == BUFR_NEGZ) then
1203 ielem = 5
1204 end if
1205
1206 if ((zlev * MPC_MBAR_PER_PA_R8) >= xstd_ua_ai_sw(1,1)) then
1207
1208 call obs_bodySet_r(obsSpaceData, OBS_OER, bodyIndex, xstd_ua_ai_sw(1, ielem))
1209
1210 else if ((zlev * MPC_MBAR_PER_PA_R8) <= xstd_ua_ai_sw(19, 1)) then
1211
1212 call obs_bodySet_r(obsSpaceData, OBS_OER, bodyIndex, xstd_ua_ai_sw(19, ielem))
1213
1214 else
1215
1216 do jn = 1, 18
1217 if ((zlev * MPC_MBAR_PER_PA_R8) >= xstd_ua_ai_sw(jn + 1, 1)) exit
1218 end do
1219
1220 zwb = log((zlev * MPC_MBAR_PER_PA_R8) / xstd_ua_ai_sw(jn, 1)) / log(xstd_ua_ai_sw(jn + 1, 1) / xstd_ua_ai_sw(jn, 1))
1221 zwt = 1.0D0 - zwb
1222
1223 call obs_bodySet_r(obsSpaceData, OBS_OER, bodyIndex, &
1224 zwt * xstd_ua_ai_sw(jn, ielem) + &
1225 zwb * xstd_ua_ai_sw(jn + 1, ielem))
1226
1227 end if
1228
1229 end if
1230
1231 !***********************************************************************
1232 ! AMV, AIREP, AMDAR DATA
1233 !***********************************************************************
1234
1235 else if (cfam == 'AI' .or. cfam == 'SW') then
1236
1237 zlev=obs_bodyElem_r(obsSpaceData, OBS_PPP, bodyIndex)
1238
1239 if (codeType == codtyp_get_codtyp('windsbufr')) then ! AMV
1240
1241 if ((ityp == BUFR_NEUU) .or. (ityp == BUFR_NEVV)) ielem = 11
1242
1243 else if (codeType == codtyp_get_codtyp('airep')) then ! AIREP
1244 if ((ityp == BUFR_NEUU) .or. (ityp == BUFR_NEVV)) then
1245 ielem = 7
1246 else if (ityp == BUFR_NETT) then
1247 ielem = 6
1248 end if
1249
1250 else if (codeType == codtyp_get_codtyp('amdar') .or. &
1251 codeType == codtyp_get_codtyp('acars') .or. &
1252 codeType == codtyp_get_codtyp('ads')) then
1253 if (ityp == BUFR_NEUU .or. ityp == BUFR_NEVV) then
1254 ielem = 10
1255 else if (ityp == BUFR_NETT) then
1256 ielem = 8
1257 else if (ityp == BUFR_NEES) then
1258 ielem = 9
1259 end if
1260 end if
1261
1262 if ((zlev * MPC_MBAR_PER_PA_R8) >= xstd_ua_ai_sw(1, 1)) then
1263 call obs_bodySet_r(obsSpaceData, OBS_OER, bodyIndex, xstd_ua_ai_sw(1, ielem))
1264 else if ((zlev * MPC_MBAR_PER_PA_R8) <= xstd_ua_ai_sw(19, 1)) then
1265 call obs_bodySet_r(obsSpaceData, OBS_OER, bodyIndex, xstd_ua_ai_sw(19, ielem))
1266 else
1267 do jn = 1, 18
1268 if ((zlev * MPC_MBAR_PER_PA_R8) >= xstd_ua_ai_sw(jn + 1, 1)) exit
1269 end do
1270
1271 zwb = log((zlev * MPC_MBAR_PER_PA_R8) / xstd_ua_ai_sw(jn, 1)) / log(xstd_ua_ai_sw(jn + 1, 1) / xstd_ua_ai_sw(jn, 1))
1272 zwt = 1.0D0 - zwb
1273
1274 call obs_bodySet_r(obsSpaceData, OBS_OER, bodyIndex, &
1275 zwt * xstd_ua_ai_sw(jn, ielem) + &
1276 zwb * xstd_ua_ai_sw(jn + 1, ielem))
1277
1278 end if
1279
1280 !***********************************************************************
1281 ! SURFACE DATA
1282 !***********************************************************************
1283
1284 else if (cfam == 'SF') then
1285 icodtyp = 1 ! Default values
1286 if (codeType == codtyp_get_codtyp('synopnonauto')) icodtyp = 2 ! SYNOP
1287 if (codeType == codtyp_get_codtyp('shipnonauto')) icodtyp = 3 ! SHIP NON-AUTOMATIQUE
1288 if (codeType == codtyp_get_codtyp('synopmobil')) icodtyp = 4 ! DRIBU
1289 if (codeType == codtyp_get_codtyp('drifter')) icodtyp = 5 ! DRifTER
1290 if (codeType == codtyp_get_codtyp('synoppatrol')) icodtyp = 6 ! STATION AUTOMATIQUE
1291 if (codeType == codtyp_get_codtyp('asynopauto')) icodtyp = 7 ! ASYNOP
1292 if (codeType == codtyp_get_codtyp('ashipauto')) icodtyp = 8 ! ASHIP
1293 if (ityp == BUFR_NEUU .or. ityp == BUFR_NEVV .or. &
1294 ityp == BUFR_NEGZ .or. ityp == BUFR_NETT .or. ityp == BUFR_NEES) icodtyp = 9 ! Others
1295 if (ityp == BUFR_NEUS .or. ityp == BUFR_NEVS) then
1296 call obs_bodySet_r(obsSpaceData, OBS_OER, bodyIndex, xstd_sf(icodtyp, 4))
1297 else if (ityp == BUFR_NETS) then
1298 call obs_bodySet_r(obsSpaceData, OBS_OER, bodyIndex, xstd_sf(icodtyp, 2))
1299 else if (ityp == BUFR_NESS) then
1300 call obs_bodySet_r(obsSpaceData, OBS_OER, bodyIndex, xstd_sf(icodtyp, 3))
1301 else if (ityp == BUFR_NEPS) then
1302 call obs_bodySet_r(obsSpaceData, OBS_OER, bodyIndex, xstd_sf(icodtyp, 1))
1303 else if (ityp == BUFR_NEPN) then
1304 if (icodtyp == 2 .or. icodtyp == 7) then
1305 call obs_bodySet_r(obsSpaceData, OBS_OER, bodyIndex, xstd_sf(1, 1))
1306 else
1307 call obs_bodySet_r(obsSpaceData, OBS_OER, bodyIndex, xstd_sf(icodtyp, 1))
1308 end if
1309 else if ((ityp == BUFR_NEUU) .or. (ityp == BUFR_NEVV))then
1310 call obs_bodySet_r(obsSpaceData, OBS_OER, bodyIndex, xstd_sf(icodtyp, 4))
1311 else if (ityp == BUFR_NEGZ) then
1312 call obs_bodySet_r(obsSpaceData, OBS_OER, bodyIndex, xstd_sf(icodtyp, 1))
1313 else if (ityp == BUFR_NETT) then
1314 call obs_bodySet_r(obsSpaceData, OBS_OER, bodyIndex, xstd_sf(icodtyp, 2))
1315 else if (ityp == BUFR_NEES) then
1316 call obs_bodySet_r(obsSpaceData, OBS_OER, bodyIndex, xstd_sf(icodtyp, 3))
1317 else if (ityp == bufr_vis .or. ityp == bufr_logVis) then
1318 if (surfaceObsTypeNumber >= 5) then
1319 call obs_bodySet_r(obsSpaceData, OBS_OER, bodyIndex, xstd_sf(icodtyp, 5))
1320 else
1321 call utl_abort('oer_fillObsErrors: observation error missing for visibility')
1322 end if
1323 else if (ityp == bufr_gust) then
1324 if (surfaceObsTypeNumber >= 6) then
1325 call obs_bodySet_r(obsSpaceData, OBS_OER, bodyIndex, xstd_sf(icodtyp, 6))
1326 else
1327 call utl_abort('oer_fillObsErrors: observation error missing for wind gust')
1328 end if
1329 else if (ityp == BUFR_NEFS) then ! SAR wind speed
1330 call obs_bodySet_r(obsSpaceData, OBS_OER, bodyIndex, 2.0d0)
1331 end if
1332
1333 !***********************************************************************
1334 ! GPS RO DATA
1335 !***********************************************************************
1336
1337 else if (cfam == 'RO') then
1338 ! Process only refractivity data (codtyp 169)
1339 if (codeType == codtyp_get_codtyp('ro')) then
1340 if (ityp == BUFR_NEPS) then
1341 call obs_bodySet_r(obsSpaceData, OBS_OER, bodyIndex, 50.D0)
1342 end if
1343 if (ityp == BUFR_NETT) then
1344 call obs_bodySet_r(obsSpaceData, OBS_OER, bodyIndex, 10.D0)
1345 end if
1346 if (ityp == BUFR_NERF) then
1347 call obs_bodySet_r(obsSpaceData, OBS_OER, bodyIndex, 500.D0)
1348 end if
1349 if (ityp == BUFR_NEBD) then
1350 call obs_bodySet_r(obsSpaceData, OBS_OER, bodyIndex, 0.08D0)
1351 end if
1352 end if
1353
1354 !***********************************************************************
1355 ! GB-GPS SFC MET DATA
1356 !***********************************************************************
1357
1358 ! ERRORS ARE SET TO SYNO SFC OBS ERRORS FROM S/R SUCOVO
1359 ! AND WEIGHTED BY FACTOR YSFERRWGT FOR 3D-VAR FGAT OR 4D-VAR ASSIM.
1360 ! OF TIME-SERIES (YSFERRWGT = 1.0 FOR 3D THINNING)
1361 !
1362 else if (cfam == 'GP') then
1363
1364 if (ityp == BUFR_NEPS) then ! Psfc Error (Pa)
1365 call obs_bodySet_r(obsSpaceData, OBS_OER, bodyIndex, xstd_sf(2, 1))
1366 end if
1367 if (ityp == BUFR_NETS) then ! Tsfc Error (K)
1368 call obs_bodySet_r(obsSpaceData, OBS_OER, bodyIndex, xstd_sf(2, 2))
1369 end if
1370 if (ityp == BUFR_NESS) then ! T-Td Error (K)
1371 call obs_bodySet_r(obsSpaceData, OBS_OER, bodyIndex, xstd_sf(2, 3))
1372 end if
1373 ! ZTD Error (m) (value is formal error, real error set later in s/r seterrgpsgb)
1374 ! If error is missing, set to dummy value (1 m).
1375 if (ityp == BUFR_NEZD) then
1376 if (obs_bodyElem_r(obsSpaceData, OBS_OER, bodyIndex) <= 0.0D0) then
1377 call obs_bodySet_r(obsSpaceData, OBS_OER, bodyIndex, 1.0D0)
1378 end if
1379 end if
1380
1381 !***********************************************************************
1382 ! SCATTEROMETER, WIND PROFILER DATA
1383 !***********************************************************************
1384
1385 else if (cfam == 'SC') then
1386
1387 call obs_bodySet_r(obsSpaceData, OBS_OER, bodyIndex, xstd_sc(1))
1388
1389 else if (cfam == 'PR') then
1390
1391 zlev= obs_bodyElem_r(obsSpaceData, OBS_PPP, bodyIndex) - obs_headElem_r(obsSpaceData, OBS_ALT, headerIndex)
1392
1393 if (zlev >= 6000.) then
1394 call obs_bodySet_r(obsSpaceData, OBS_OER, bodyIndex, xstd_pr(2))
1395 else
1396 call obs_bodySet_r(obsSpaceData, OBS_OER, bodyIndex, xstd_pr(1))
1397 end if
1398
1399 !***********************************************************************
1400 ! ALADIN HORIZONTAL LINE-OF-SIGHT WIND DATA
1401 !***********************************************************************
1402
1403 else if (cfam == 'AL') then
1404
1405 ! TEMPORARILY, hard-code the observation error of AL winds to 2 m/s
1406 call obs_bodySet_r(obsSpaceData, OBS_OER, bodyIndex, 2.0d0)
1407
1408 !***********************************************************************
1409 ! CONSTITUENT DATA (OZONE AND OTHER CHEMICALS)
1410 !***********************************************************************
1411
1412 else if (cfam == 'CH') then
1413
1414 ! Process only retrieved constituent data
1415 ! Also, exclude BUFR_SCALE_EXPONENT element as a data value!
1416
1417 if (codeType == codtyp_get_codtyp('CHEMREMOTE') .or. codeType == codtyp_get_codtyp('CHEMINSITU')) then
1418
1419 ifirst = headerIndex /= header_prev
1420 if (ifirst) then
1421
1422 header_prev = headerIndex
1423
1424 ! Subtract the number of occurences of code BUFR_SCALE_EXPONENT from number of levels
1425 do ji = 0, nlev - 1
1426 if (obs_bodyElem_i(obsSpaceData, OBS_VNM, idata + ji) == BUFR_SCALE_EXPONENT) nlev = nlev-1
1427 end do
1428 end if
1429
1430 zlev = obs_bodyElem_r(obsSpaceData,OBS_PPP,bodyIndex)
1431
1432 ilev = 0
1433 do ji = idata, bodyIndex
1434 if (obs_bodyElem_i(obsSpaceData, OBS_VNM, ji) /= BUFR_SCALE_EXPONENT) ilev = ilev+1
1435 end do
1436
1437 obs_err_stddev = chm_get_obs_err_stddev(cstnid, nlev, ityp, zlat, zlon, idate, itime, zval, zlev, ilev, ifirst)
1438 call obs_bodySet_r(obsSpaceData, OBS_OER, bodyIndex, obs_err_stddev)
1439
1440 end if
1441
1442 !***********************************************************************
1443 ! Sea Surface Temperature
1444 !***********************************************************************
1445
1446 else if (cfam == 'TM') then
1447
1448 if (obs_bodyElem_i(obsSpaceData, OBS_VNM, bodyIndex) /= bufr_sst) cycle BODY
1449
1450 if (codeType == codtyp_get_codtyp('drifter') .or. codeType == codtyp_get_codtyp('shipnonauto') .or. &
1451 codeType == codtyp_get_codtyp('ashipauto') .or. codeType == codtyp_get_codtyp('pseudosfc') ) then
1452
1453 unsupportedCodeType = .true.
1454 dataset_loop: do indexDataset = 1, numberSSTDatasets
1455 if(codeType == SSTdataParams(indexDataset)%codeType) then
1456 unsupportedCodeType = .false.
1457 call obs_bodySet_r(obsSpaceData, OBS_OER, bodyIndex, SSTdataParams(indexDataset)%dayError)
1458 exit dataset_loop
1459 end if
1460 end do dataset_loop
1461
1462 if(unsupportedCodeType) then
1463 write(*,'(a,i5,2a)') 'oer_fillObsErrors: unsupported SST data, codtype ', codeType,', cstnid ', cstnid,' found in dataset_loop!'
1464 call utl_abort('oer_fillObsErrors: unsupported codeType!')
1465 end if
1466
1467 else if (codeType == codtyp_get_codtyp('satob')) then
1468
1469 solarZenith = obs_headElem_r(obsSpaceData, obs_sun, headerIndex)
1470 if (solarZenith == MPC_missingValue_R8) then
1471 write(*,*) 'oer_fillObsErrors: Solar zenith value is missing for satellite SST data ', &
1472 cstnid, ', headerIndex: ', headerIndex,', bodyIndex: ', bodyIndex
1473 call utl_abort('oer_fillObsErrors: Solar zenith value is missing!')
1474 end if
1475
1476 unsupportedSensor = .true.
1477 sensor_loop: do indexSensor = 1, numberSSTDatasets
1478 if (cstnid == trim(SSTdataParams(indexSensor)%sensor)) then
1479 unsupportedSensor = .false.
1480 if (solarZenith <= 90.) then ! day
1481 call obs_bodySet_r(obsSpaceData, OBS_OER, bodyIndex, SSTdataParams(indexSensor)%dayError)
1482 else ! night
1483 call obs_bodySet_r(obsSpaceData, OBS_OER, bodyIndex, SSTdataParams(indexSensor)%nightError)
1484 end if
1485 exit sensor_loop
1486 end if
1487 end do sensor_loop
1488
1489 if(unsupportedSensor) then
1490 write(*,'(3a)') 'oer_fillObsErrors: unsupported satellite SST data sensor ', cstnid,' found in sensor_loop!'
1491 call utl_abort('oer_fillObsErrors: unsupported satellite SST data sensor!')
1492 end if
1493
1494 else
1495
1496 write(*,'(a,i5,2a)') 'oer_fillObsErrors: unsupported SST data, codtype ', codeType,', cstnid ', cstnid,' found!'
1497 call utl_abort('oer_fillObsErrors: unsupported codeType!')
1498
1499 end if
1500
1501 !***********************************************************************
1502 ! Sea Ice Concentration
1503 !***********************************************************************
1504
1505 else if (cfam == 'GL') then
1506
1507 if (index(cstnid,'DMSP') == 1) then
1508 select case(cstnid)
1509 case('DMSP15')
1510 call obs_bodySet_r(obsSpaceData, OBS_OER, bodyIndex, xstd_sic(1))
1511 case('DMSP16','DMSP17','DMSP18')
1512 call obs_bodySet_r(obsSpaceData, OBS_OER, bodyIndex, xstd_sic(2))
1513 case DEFAULT
1514 call utl_abort('oer_fillObsErrors: UNKNOWN station id: '//cstnid)
1515 end select
1516 else if (cstnid == 'GCOM-W1') then
1517 call obs_bodySet_r(obsSpaceData, OBS_OER, bodyIndex, xstd_sic(3))
1518 else if (cstnid(1:6) == 'METOP-') then
1519 if (ityp == BUFR_ICEP) then
1520 call obs_bodySet_r(obsSpaceData, OBS_OER, bodyIndex, xstd_sic(4))
1521 else if (ityp == BUFR_ICES) then
1522 ! This is backscatter anisotropy, obs-error will be set in oer_setErrBackScatAnisIce
1523 ! because the need of background ice concentration.
1524 ! For now just put a large positive value
1525 call obs_bodySet_r(obsSpaceData, OBS_OER, bodyIndex, 1000.0)
1526 else
1527 call utl_abort('oer_fillObsErrors: UNKNOWN varno: '// trim(utl_str(ityp)) // 'station id: '//cstnid)
1528 end if
1529 else if (cstnid == 'noaa-19') then
1530 call obs_bodySet_r(obsSpaceData, OBS_OER, bodyIndex, xstd_sic(5))
1531 else if (cstnid == 'CIS_DAILY') then
1532 call obs_bodySet_r(obsSpaceData, OBS_OER, bodyIndex, xstd_sic(6))
1533 else if (cstnid == 'RS1_IMG') then
1534 call obs_bodySet_r(obsSpaceData, OBS_OER, bodyIndex, xstd_sic(7))
1535 else if (codeType == 178) then ! lake ice
1536 call obs_bodySet_r(obsSpaceData, OBS_OER, bodyIndex, xstd_sic(8))
1537 else if (cstnid == 'CIS_REGIONAL') then
1538 call obs_bodySet_r(obsSpaceData, OBS_OER, bodyIndex, xstd_sic(9))
1539 else if (index(cstnid,'RCM') == 1) then
1540 ! For RCM, obs-error comes from SQLite file
1541 else
1542 call utl_abort('oer_fillObsErrors: UNKNOWN station id: '//cstnid)
1543 end if
1544
1545 !***********************************************************************
1546 ! Hydrology
1547 !***********************************************************************
1548 else if (cfam == 'HY') then
1549
1550 if (obs_bodyElem_i(obsSpaceData, OBS_VNM, bodyIndex) /= bufr_riverFlow) cycle BODY
1551
1552 if (codeType == 12) then ! pseudo-SYNOP
1553 obsValue = obs_bodyElem_r(obsSpaceData, OBS_VAR, bodyIndex)
1554 obsStdDevError = obsValue * xstd_hydro(1)
1555 write(*,*) 'Hydro observation std dev error: ', bodyIndex, obsValue, xstd_hydro(1), obsStdDevError
1556 call obs_bodySet_r(obsSpaceData, OBS_OER, bodyIndex, obsStdDevError)
1557 else
1558 write(*,*) 'oer_fillObsErrors: unsupported codeType for hydro data found in the observations: ', codeType
1559 call utl_abort('oer_fillObsErrors: unsupported codeType')
1560 end if
1561
1562 else if (cfam == 'RA') then
1563
1564 if (ityp == BUFR_radarPrecip .or. ityp == BUFR_logRadarPrecip) then
1565 ! Temporary hardcoded value for log-transformed radar precipitation
1566 call obs_bodySet_r(obsSpaceData, OBS_OER, bodyIndex, 1.0D0)
1567 else if (ityp == bufr_radvel) then
1568 ! Temporary hardcoded value for radar Doppler velocity
1569 call obs_bodySet_r(obsSpaceData, OBS_OER, bodyIndex, 1.0D0)
1570 else
1571 write(*,*) 'varnum = ', ityp
1572 call utl_abort('oer_fillObsErrors: unknown varnum for RA family')
1573 end if
1574
1575 else
1576
1577 write(*,*) 'oer_fillObsErrors: UNKNOWN DATA FAMILY:',cfam
1578
1579 end if
1580
1581 !***********************************************************************
1582 ! Check for case where error should have been set but was
1583 ! not. 3dvar will abort in this case.
1584 !***********************************************************************
1585
1586 obsOER = obs_bodyElem_r(obsSpaceData,OBS_OER,bodyIndex)
1587 if (obsOER <= 0.0D0) then
1588
1589 write(*,*)' PROBLEM OBSERR VARIANCE FAM= ',cfam
1590
1591 write(*,'(1X,"STNID= ",A10,"codeType= ",I5," LAT= ",F10.2," LON = ",F10.2)') &
1592 cstnid, &
1593 codeType, &
1594 zlat*MPC_DEGREES_PER_RADIAN_R8, &
1595 zlon*MPC_DEGREES_PER_RADIAN_R8
1596
1597 obsPPP = obs_bodyElem_r(obsSpaceData,OBS_PPP,bodyIndex)
1598 write(*,'(1X,"ELEMENT= ",I6," LEVEL= ",F10.2," OBSERR = ",E10.2)') &
1599 ityp, obsPPP, obsOER
1600
1601 call utl_abort('oer_fillObsErrors: PROBLEM OBSERR VARIANCE.')
1602
1603 end if
1604
1605 else
1606
1607 call obs_bodySet_r(obsSpaceData, OBS_OER, bodyIndex, MPC_missingValue_R8)
1608
1609 end if ! end of iass == obs_assimilated
1610
1611 end do BODY
1612
1613 end do HEADER
1614
1615 write(*,'(10X,"Fill_obs_errors: Done")')
1616 write(*,'(10X,"---------------",/)')
1617
1618 contains
1619
1620 !--------------------------------------------------------------------------
1621 ! calcStateDepObsErr
1622 !--------------------------------------------------------------------------
1623 function calcStateDepObsErr(cldPredThresh1, cldPredThresh2, &
1624 errThresh1, errThresh2, cldPredUsed) result(sigmaObsErrUsed)
1625 !
1626 !:Purpose: Calculate state-dependent observation error.
1627 !
1628 implicit none
1629
1630 ! Arguments:
1631 real(8), intent(in) :: cldPredThresh1 ! first cloud predictor threshold
1632 real(8), intent(in) :: cldPredThresh2 ! second cloud predictor threshold
1633 real(8), intent(in) :: errThresh1 ! sigmaO corresponding to first cloud predictor threshold
1634 real(8), intent(in) :: errThresh2 ! sigmaO corresponding to second cloud predictor threshold
1635 real(8), intent(in) :: cldPredUsed ! cloud predictor for the obs
1636 ! Result:
1637 real(8) :: sigmaObsErrUsed ! estimated sigmaO for the obs
1638
1639 if (cldPredUsed <= cldPredThresh1) then
1640 sigmaObsErrUsed = errThresh1
1641 else if (cldPredUsed > cldPredThresh1 .and. &
1642 cldPredUsed <= cldPredThresh2) then
1643 sigmaObsErrUsed = errThresh1 + &
1644 (errThresh2 - errThresh1) / &
1645 (cldPredThresh2 - cldPredThresh1) * &
1646 (cldPredUsed - cldPredThresh1)
1647 else
1648 sigmaObsErrUsed = errThresh2
1649 end if
1650
1651 end function calcStateDepObsErr
1652
1653 !--------------------------------------------------------------------------
1654 ! computeCloudPredictor
1655 !--------------------------------------------------------------------------
1656 function computeCloudPredictor(sensorIndex, headerIndex) result(cldPredUsed)
1657 !
1658 !:Purpose: Compute cloud predictor to use for state-dependent observation error.
1659 !
1660 implicit none
1661
1662 ! Arguments:
1663 integer, intent(in) :: sensorIndex ! index of sensor in tvs_nsensors
1664 integer, intent(in) :: headerIndex
1665 ! Result:
1666 real(8) :: cldPredUsed ! cloud predictor. CLW/SI for all-sky temperature/humidity
1667
1668 ! Locals:
1669 integer :: platformId
1670 integer :: satelliteId
1671 integer :: instrumId
1672 real(8) :: clwObs, clwFG
1673 real(8) :: siObs, siFG
1674 real(4) :: cldPredUsed_r4
1675
1676 platformId = tvs_platforms(sensorIndex)
1677 satelliteId = tvs_satellites(sensorIndex)
1678 instrumId = tvs_instruments(sensorIndex)
1679
1680 if (tvs_isInstrumAllskyTtAssim(instrumId)) then
1681 clwObs = obs_headElem_r(obsSpaceData, OBS_CLWO, headerIndex)
1682 clwFG = obs_headElem_r(obsSpaceData, OBS_CLWB, headerIndex)
1683 cldPredUsed = 0.5D0 * (clwObs + clwFG)
1684 cldPredUsed_r4 = real(cldPredUsed)
1685
1686 ! check to ensure CLW is retrieved and properly set
1687 if (cldPredUsed_r4 < minRetrievableClwValue_r4 .or. &
1688 cldPredUsed_r4 > maxRetrievableClwValue_r4) then
1689 write(*,*) 'This observation should have been rejected ', &
1690 'for all-sky temperature in background check!'
1691 write(*,*) 'computeCloudPredictor: platformId=', platformId, &
1692 ', satelliteId=', satelliteId, ', instrumId=', instrumId
1693 write(*,*) 'computeCloudPredictor: clwObs=', clwObs, &
1694 ', clwFG=', clwFG
1695 call utl_abort('computeCloudPredictor: not usable to define obs error with CLW')
1696 end if
1697
1698 else if (tvs_isInstrumAllskyHuAssim(instrumId)) then
1699 siObs = obs_headElem_r(obsSpaceData, OBS_SIO, headerIndex)
1700 siFG = obs_headElem_r(obsSpaceData, OBS_SIB, headerIndex)
1701 cldPredUsed = 0.5D0 * (siObs + siFG)
1702 cldPredUsed_r4 = real(cldPredUsed)
1703
1704 ! check to ensure SI is retrieved and properly set
1705 if (cldPredUsed_r4 < minRetrievableSiValue_r4 .or. &
1706 cldPredUsed_r4 > maxRetrievableSiValue_r4) then
1707 write(*,*) 'This observation should have been rejected ', &
1708 'for all-sky humidity in background check!'
1709 write(*,*) 'computeCloudPredictor: platformId=', platformId, &
1710 ', satelliteId=', satelliteId, ', instrumId=', instrumId
1711 write(*,*) 'computeCloudPredictor: siObs=', siObs, &
1712 ', siFG=', siFG
1713 call utl_abort('computeCloudPredictor: not usable to define obs error with SI')
1714 end if
1715
1716 else
1717 call utl_abort('computeCloudPredictor: instrum is not TT or HU allSky assimilation.')
1718 end if
1719
1720 end function computeCloudPredictor
1721
1722 end subroutine oer_fillObsErrors
1723
1724 !--------------------------------------------------------------------------
1725 ! oer_inflateErrAllsky
1726 !--------------------------------------------------------------------------
1727 subroutine oer_inflateErrAllsky(obsSpaceData, bodyIndex, ompOmaObsColumn, beSilent_opt)
1728 !
1729 !:Purpose: Update OBS_OER with inflated state dependant observation error for all-sky
1730 ! temperature/humidity assimilation.
1731 !
1732 implicit none
1733
1734 ! Arguments:
1735 type(struct_obs), intent(inout) :: obsSpaceData
1736 integer, intent(in) :: bodyIndex
1737 integer, intent(in) :: ompOmaObsColumn ! obsSpaceData OBS_OMP or OBS_OMA column
1738 logical, optional, intent(in) :: beSilent_opt ! prints extra info to listing if .true.
1739
1740 ! Locals:
1741 integer :: headerIndex
1742 integer :: channelNumber_withOffset
1743 integer :: channelNumber, channelIndex
1744 integer :: tovsIndex, sensorIndex
1745 logical :: surfTypeIsWater
1746 real(8) :: clwObs
1747 real(8) :: clwFG
1748 real(8) :: siObs
1749 real(8) :: siFG
1750 real(8) :: deltaE1
1751 real(8) :: deltaE2
1752 real(8) :: ompValue
1753 real(8) :: sigmaObsBeforeInflation
1754 real(8) :: sigmaObsAfterInflation
1755 logical :: chanIsAllskyTt, chanIsAllskyHu
1756 logical :: beSilent
1757
1758 if (present(beSilent_opt)) then
1759 beSilent = beSilent_opt
1760 else
1761 beSilent = .true.
1762 end if
1763
1764 headerIndex = obs_bodyElem_i(obsSpaceData, OBS_HIND, bodyIndex)
1765 tovsIndex = tvs_tovsIndex(headerIndex)
1766 sensorIndex = tvs_lsensor(tovsIndex)
1767
1768 call tvs_getChannelNumIndexFromPPP(obsSpaceData, headerIndex, bodyIndex, &
1769 channelNumber, channelIndex)
1770 channelNumber_withOffset = channelNumber + tvs_channelOffset(sensorIndex)
1771
1772 surfTypeIsWater = (tvs_ChangedStypValue(obsSpaceData,headerIndex) == surftype_sea)
1773
1774 call chanIsAllsky(obsSpaceData, bodyIndex, chanIsAllskyTt, chanIsAllskyHu)
1775
1776 if (chanIsAllskyTt) then
1777 if (.not. surfTypeIsWater .or. &
1778 (.not. mwAllskyTtInflateByOmp .and. .not. mwAllskyTtInflateByClwDiff)) then
1779 return
1780 end if
1781
1782 clwObs = obs_headElem_r(obsSpaceData, OBS_CLWO, headerIndex)
1783 clwFG = obs_headElem_r(obsSpaceData, OBS_CLWB, headerIndex)
1784
1785 sigmaObsBeforeInflation = obs_bodyElem_r(obsSpaceData, OBS_OER, bodyIndex)
1786 ompValue = obs_bodyElem_r(obsSpaceData, ompOmaObsColumn, bodyIndex)
1787
1788 if (.not. beSilent) then
1789 write(*,*) 'oer_inflateErrAllsky: chanIsAllskyTt', &
1790 ', headerIndex=', headerIndex, &
1791 ', sensorIndex=', sensorIndex, &
1792 ', chan_noOff=', channelNumber_withOffset, &
1793 ', chan_no=', channelNumber
1794 write(*,*) 'oer_inflateErrAllsky: clwObs=', clwObs, &
1795 ', clwFG=', clwFG, ', OMP=', ompValue
1796 end if
1797
1798 ! error inflation for cloud placement
1799 deltaE1 = 0.0D0
1800 if (mwAllskyTtInflateByOmp .and. &
1801 ((clwObs - clearCldPredThresh(channelNumber_withOffset,sensorIndex)) * &
1802 (clwFG - clearCldPredThresh(channelNumber_withOffset,sensorIndex)) < 0) .and. &
1803 abs(clwObs - clwFG) >= 0.005) then
1804 deltaE1 = abs(ompValue)
1805 end if
1806
1807 ! error inflation due to cloud liquid water difference
1808 deltaE2 = 0.0D0
1809 if (mwAllskyTtInflateByClwDiff) then
1810 deltaE2 = inflateErrAllskyCoeff(sensorIndex) * abs(clwObs - clwFG) * &
1811 sigmaObsBeforeInflation
1812 end if
1813 deltaE2 = min(deltaE2,3.5D0 * sigmaObsBeforeInflation)
1814
1815 else if (chanIsAllskyHu) then
1816 if (.not. surfTypeIsWater .or. &
1817 (.not. mwAllskyHuInflateByOmp .and. .not. mwAllskyHuInflateBySiDiff)) then
1818 return
1819 end if
1820
1821 siObs = obs_headElem_r(obsSpaceData, OBS_SIO, headerIndex)
1822 siFG = obs_headElem_r(obsSpaceData, OBS_SIB, headerIndex)
1823
1824 sigmaObsBeforeInflation = obs_bodyElem_r(obsSpaceData, OBS_OER, bodyIndex)
1825 ompValue = obs_bodyElem_r(obsSpaceData, ompOmaObsColumn, bodyIndex)
1826
1827 if (.not. beSilent) then
1828 write(*,*) 'oer_inflateErrAllsky: chanIsAllskyHu', &
1829 ', headerIndex=', headerIndex, &
1830 ', sensorIndex=', sensorIndex, &
1831 ', chan_noOff=', channelNumber_withOffset, &
1832 ', chan_no=', channelNumber
1833 write(*,*) 'oer_inflateErrAllsky: siObs=', siObs, &
1834 ', siFG=', siFG, ', OMP=', ompValue
1835 end if
1836
1837 ! error inflation for cloud placement
1838 deltaE1 = 0.0D0
1839 if (mwAllskyHuInflateByOmp .and. &
1840 ((siObs - clearCldPredThresh(channelNumber_withOffset,sensorIndex)) * &
1841 (siFG - clearCldPredThresh(channelNumber_withOffset,sensorIndex)) < 0) .and. &
1842 abs(siObs - siFG) >= 1.0) then
1843 deltaE1 = abs(ompValue)
1844 end if
1845
1846 ! error inflation due to scattering-index difference
1847 deltaE2 = 0.0D0
1848 if (mwAllskyHuInflateBySiDiff) then
1849 deltaE2 = inflateErrAllskyCoeff(sensorIndex) * abs(siObs - siFG) * &
1850 sigmaObsBeforeInflation
1851 end if
1852 deltaE2 = min(deltaE2,3.5D0 * sigmaObsBeforeInflation)
1853
1854 else
1855 return
1856 end if
1857
1858 sigmaObsAfterInflation = sqrt(sigmaObsBeforeInflation ** 2 + &
1859 (deltaE1 + deltaE2) ** 2)
1860
1861 if (.not. beSilent) then
1862 write(*,*) 'oer_inflateErrAllsky: deltaE1=', deltaE1, &
1863 ', deltaE2=', deltaE2
1864 write(*,*) 'oer_inflateErrAllsky: sigmaObs=', &
1865 sigmaObsBeforeInflation, ', sigmaObsInflated=', sigmaObsAfterInflation
1866 end if
1867
1868 call obs_bodySet_r(obsSpaceData, OBS_OER, bodyIndex, sigmaObsAfterInflation)
1869
1870 end subroutine oer_inflateErrAllsky
1871
1872 !--------------------------------------------------------------------------
1873 ! chanIsAllsky
1874 !--------------------------------------------------------------------------
1875 subroutine chanIsAllsky(obsSpaceData, bodyIndex, chanIsAllskyTt, chanIsAllskyHu)
1876 !
1877 !:Purpose: Determine if the tovs instrument/channel combination is all-sky.
1878 !
1879 implicit none
1880
1881 ! Arguments:
1882 type(struct_obs), intent(in) :: obsSpaceData
1883 integer, intent(in) :: bodyIndex
1884 logical, intent(out) :: chanIsAllskyTt ! .true. if channel is all-sky temperature
1885 logical, intent(out) :: chanIsAllskyHu ! .true. if channel is all-sky humidity
1886
1887 ! Locals:
1888 integer :: headerIndex
1889 integer :: channelNumber_withOffset
1890 integer :: channelNumber, channelIndex
1891 integer :: tovsIndex, sensorIndex, instrumId
1892
1893 headerIndex = obs_bodyElem_i(obsSpaceData, OBS_HIND, bodyIndex)
1894 tovsIndex = tvs_tovsIndex(headerIndex)
1895 sensorIndex = tvs_lsensor(tovsIndex)
1896 instrumId = tvs_instruments(sensorIndex)
1897
1898 call tvs_getChannelNumIndexFromPPP(obsSpaceData, headerIndex, bodyIndex, &
1899 channelNumber, channelIndex)
1900 channelNumber_withOffset = channelNumber + tvs_channelOffset(sensorIndex)
1901
1902 chanIsAllskyTt = .false.
1903 chanIsAllskyHu = .false.
1904 if (.not. tvs_mwAllskyAssim .or. &
1905 .not. useStateDepSigmaObs(channelNumber_withOffset,sensorIndex)) return
1906
1907 if (tvs_isInstrumAllskyTtAssim(instrumId)) then
1908 chanIsAllskyTt = .true.
1909 end if
1910 if (tvs_isInstrumAllskyHuAssim(instrumId)) then
1911 chanIsAllskyHu = .true.
1912 end if
1913 if (tvs_isInstrumAllskyTtAssim(instrumId) .and. tvs_isInstrumAllskyHuAssim(instrumId)) then
1914 call utl_abort('chanIsAllsky: all-sky TtHu is not yet implemented')
1915 end if
1916
1917 end subroutine chanIsAllsky
1918
1919 !--------------------------------------------------------------------------
1920 ! readOerFromObsFileForSW
1921 !--------------------------------------------------------------------------
1922 subroutine readOerFromObsFileForSW(obsSpaceData)
1923 !
1924 !:Purpose: Read observation errors for AMVs from the obsFiles so as to use
1925 ! the values computed by a previous MIDAS execution (e.g. GDPS).
1926 !
1927 implicit none
1928
1929 ! Arguments:
1930 type(struct_obs), intent(inout) :: obsSpaceData
1931
1932 ! Locals:
1933 character(len=1060) :: filename
1934 type(BURP_FILE) :: fileIn
1935 type(BURP_BLOCK) :: blkoer
1936 type(BURP_RPT) :: report
1937 character(len=9) :: stnid
1938 integer(kind=int_def) :: error, ref_rpt
1939 integer :: numLevels, numValues, numReports, obsCount
1940 logical :: fileFound
1941 integer :: levelIndex, reportIndex, obsIndex
1942 integer :: uuIndex, vvIndex, headerIndex, bodyIndex, blockIndex, g_btyp_oer
1943 integer :: bodyIndexBeg, bodyIndexEnd
1944 real(8), allocatable :: uu_oer(:), vv_oer(:)
1945
1946 filename = obsf_getFileName('SW',fileFound)
1947 if (fileFound) then
1948 write(*,*) 'oer_readOerFromObsFileForSW: reading OER from the file: ', trim(filename)
1949 else
1950 write(*,*) 'oer_readOerFromObsFileForSW: no obsfile with SW family, returning'
1951 return
1952 end if
1953
1954 g_btyp_oer = 1134
1955
1956 call burp_init(report)
1957 call burp_init(blkoer)
1958 call burp_init(fileIn)
1959 call burp_new(fileIn, FILENAME=filename, MODE=FILE_ACC_READ, IOSTAT=error)
1960
1961 !
1962 ! Scans reports
1963 !
1964 call burp_get_property(fileIn, NRPTS = numReports, FILENAME = filename, IOSTAT = error)
1965 ref_rpt = 0
1966 obsCount = 0
1967
1968 do reportIndex = 1, numReports
1969 ref_rpt = burp_find_report(fileIn, REPORT=report, SEARCH_FROM=ref_rpt, IOSTAT=error)
1970 call burp_get_property(report, ELEV=numLevels, IOSTAT=error)
1971 obsCount = obsCount + numLevels
1972 end do
1973
1974 write(*,*) 'readOerFromObsFileForSW: numReports, obsCount = ',numReports, obsCount
1975
1976 allocate (uu_oer(obsCount))
1977 allocate (vv_oer(obsCount))
1978
1979 !
1980 ! Read observation errors from file
1981 !
1982 ref_rpt = 0
1983 obsIndex = 0
1984
1985 records_in: do
1986
1987 ref_rpt = burp_find_report(fileIn, REPORT=report, SEARCH_FROM=ref_rpt, IOSTAT=error)
1988 if (ref_rpt <0) exit
1989
1990 call burp_get_property(report, STNID=stnid,ELEV=numLevels, IOSTAT=error)
1991
1992 if (stnid(1:2) == ">>") cycle records_in
1993
1994 blockIndex = 0
1995 blockIndex = burp_find_block(report, block=blkoer, search_from=blockIndex, btyp=g_btyp_oer, bfam=10, convert=.false., iostat=error)
1996
1997 call burp_get_property(blkoer, NVAL=numValues, IOSTAT=error)
1998
1999 uuIndex = burp_find_element(blkoer, ELEMENT=BUFR_NEUU, IOSTAT=error)
2000 vvIndex = burp_find_element(blkoer, ELEMENT=BUFR_NEVV, IOSTAT=error)
2001
2002 if (uuIndex == -1) then
2003 write(*,*) 'readOerFromObsFileForSW: WARNING: wind element not found, skipping report'
2004 cycle records_in
2005 end if
2006
2007 do levelIndex = 1, numLevels
2008
2009 obsIndex = obsIndex + 1
2010 if (obsIndex > obsCount) call abort('readOerFromObsFileForSW: Something went wrong')
2011 uu_oer(obsIndex) = (burp_get_tblval(blkoer, NELE_IND=uuIndex, NVAL_IND=numValues, NT_IND=levelIndex, IOSTAT=error) - 4096)/10.0
2012 vv_oer(obsIndex) = (burp_get_tblval(blkoer, NELE_IND=vvIndex, NVAL_IND=numValues, NT_IND=levelIndex, IOSTAT=error) - 4096)/10.0
2013
2014 end do
2015
2016 end do records_in
2017
2018 !
2019 ! Clean-up
2020 !
2021 call burp_free(report)
2022 call burp_free(blkoer)
2023 call burp_free(fileIn)
2024
2025 !
2026 ! Fill obsSpaceData with observation errors from SW file
2027 !
2028 obsIndex = 0
2029 HEADER_UU: do headerIndex = 1, obs_numheader(obsSpaceData)
2030 if (obs_getFamily(obsSpaceData,headerIndex_opt=headerIndex) /= 'SW') cycle HEADER_UU
2031 bodyIndexBeg = obs_headElem_i(obsSpaceData,OBS_RLN,headerIndex)
2032 bodyIndexEnd = obs_headElem_i(obsSpaceData,OBS_NLV,headerIndex) + bodyIndexBeg - 1
2033 do bodyIndex = bodyIndexBeg, bodyIndexEnd
2034 if (obs_bodyElem_i(obsSpaceData,OBS_VNM,bodyIndex) == BUFR_NEUU) then
2035 obsIndex = obsIndex + 1
2036 if (obsIndex > obsCount) call abort('readOerFromObsFileForSW: Something went wrong')
2037 call obs_bodySet_r(obsSpaceData,OBS_OER,bodyIndex,uu_oer(obsIndex))
2038 end if
2039 end do
2040 end do HEADER_UU
2041
2042 obsIndex = 0
2043 HEADER_VV: do headerIndex = 1, obs_numheader(obsSpaceData)
2044 if (obs_getFamily(obsSpaceData,headerIndex_opt=headerIndex) /= 'SW') cycle HEADER_VV
2045 bodyIndexBeg = obs_headElem_i(obsSpaceData,OBS_RLN,headerIndex)
2046 bodyIndexEnd = obs_headElem_i(obsSpaceData,OBS_NLV,headerIndex) + bodyIndexBeg - 1
2047 do bodyIndex = bodyIndexBeg, bodyIndexEnd
2048 if (obs_bodyElem_i(obsSpaceData,OBS_VNM,bodyIndex) == BUFR_NEVV) then
2049 obsIndex = obsIndex + 1
2050 if (obsIndex > obsCount) call abort('readOerFromObsFileForSW: Something went wrong')
2051 call obs_bodySet_r(obsSpaceData,OBS_OER,bodyIndex,vv_oer(obsIndex))
2052 end if
2053 end do
2054 end do HEADER_VV
2055
2056 deallocate (uu_oer)
2057 deallocate (vv_oer)
2058
2059 end subroutine readOerFromObsFileForSW
2060
2061 !--------------------------------------------------------------------------
2062 ! oer_sw
2063 !--------------------------------------------------------------------------
2064 subroutine oer_sw(columnTrlOnTrlLev,obsSpaceData)
2065 !
2066 !:Purpose: Calculate observation errors for AMVs according to the Met-Office
2067 ! situation dependant approach.
2068 !
2069 implicit none
2070
2071 ! Arguments:
2072 type(struct_columnData), intent(in) :: columnTrlOnTrlLev
2073 type(struct_obs), intent(inout) :: obsSpaceData
2074
2075 ! Locals:
2076 integer :: headerIndex,bodyIndex,ilyr,jlev
2077 integer :: iass,ixtr,ivco,ivnm,iqiv,iqiv1,iqiv2,imet,ilsv,igav,ihav,itrn,J_SAT
2078 real(8) :: zvar,zoer
2079 real(8) :: zwb,zwt,ZOTR,ZMOD
2080 real(8) :: zlat,zlon,zlev,zpt,zpb,zpc
2081 real(8) :: SP_WGH,TO_WGH,TO_DSP,E_VHGT,E_DRIFT,E_HEIGHT
2082 character(len=4) :: varName
2083 character(len=4) :: varLevel
2084 character(len=9) :: cstnid
2085 real(8), pointer :: col_ptr_uv(:)
2086 logical :: passe_once, valeurs_defaut, print_debug
2087 logical, save :: firstCall=.true.
2088
2089 ! If requested, just read oer from the burp file (only 1st time)
2090 if(obsfile_oer_sw) then
2091 if (firstCall) then
2092 call readOerFromObsFileForSW(obsSpaceData)
2093 firstCall = .false.
2094 end if
2095 return
2096 end if
2097
2098 if(.not. new_oer_sw) return
2099
2100 valeurs_defaut = .false.
2101 passe_once = .true.
2102 print_debug = .false.
2103
2104 if (firstCall) write(*,*) "Entering subroutine oer_sw"
2105 firstCall = .false.
2106
2107 call obs_set_current_body_list(obsSpaceData, 'SW')
2108 BODY: do
2109 bodyIndex = obs_getBodyIndex(obsSpaceData)
2110 if (bodyIndex < 0) exit BODY
2111
2112 ! Only process pressure level observations flagged to be assimilated
2113 iass=obs_bodyElem_i (obsSpaceData,OBS_ASS,bodyIndex)
2114 ivco=obs_bodyElem_i (obsSpaceData,OBS_VCO,bodyIndex)
2115 if(iass /= obs_assimilated .or. ivco /= 2) cycle BODY
2116
2117 ixtr = obs_bodyElem_i (obsSpaceData,OBS_XTR,bodyIndex)
2118 ivnm = obs_bodyElem_i (obsSpaceData,OBS_VNM,bodyIndex)
2119 zvar = obs_bodyElem_r (obsSpaceData,OBS_VAR,bodyIndex)
2120 zlev = obs_bodyElem_r (obsSpaceData,OBS_PPP,bodyIndex)
2121 zoer = obs_bodyElem_r (obsSpaceData,OBS_OER,bodyIndex)
2122 headerIndex = obs_bodyElem_i (obsSpaceData,OBS_HIND,bodyIndex)
2123
2124 iqiv1 = obs_headElem_i (obsSpaceData,OBS_SWQ1,headerIndex)
2125 iqiv2 = obs_headElem_i (obsSpaceData,OBS_SWQ2,headerIndex)
2126 imet = obs_headElem_i (obsSpaceData,OBS_SWMT,headerIndex)
2127 ilsv = obs_headElem_i (obsSpaceData,OBS_SWLS,headerIndex)
2128 igav = obs_headElem_i (obsSpaceData,OBS_SWGA,headerIndex)
2129 ihav = obs_headElem_i (obsSpaceData,OBS_SWHA,headerIndex)
2130 zlat = obs_headElem_r (obsSpaceData,OBS_LAT,headerIndex)*MPC_DEGREES_PER_RADIAN_R8
2131 zlon = obs_headElem_r (obsSpaceData,OBS_LON,headerIndex)*MPC_DEGREES_PER_RADIAN_R8
2132 cstnid = obs_elem_c (obsSpaceData,'STID' ,headerIndex)
2133
2134 itrn = ilsv
2135 if(igav == 1) itrn = 2
2136 if(imet > 3) imet = 3
2137 if(ihav /= 4) ihav = 1
2138 ! Use the quality score qiv2, but if it is missing then use qiv1
2139 iqiv = iqiv2
2140 if(iqiv < 0) iqiv = iqiv1
2141
2142 if(valeurs_defaut) then
2143 E_DRIFT = 2.5
2144 E_HEIGHT = 8000.0
2145 else
2146 E_DRIFT = 7.5 - 0.05*iqiv
2147 call get_height_error(cstnid,imet,itrn,ihav,zlat,zlev,E_HEIGHT,J_SAT)
2148 end if
2149
2150 varLevel = vnl_varLevelFromVarnum(ivnm)
2151 varName = vnl_varnameFromVarnum(ivnm)
2152
2153 col_ptr_uv=>col_getColumn(columnTrlOnTrlLev,headerIndex,varName)
2154 ilyr = obs_bodyElem_i (obsSpaceData,OBS_LYR,bodyIndex)
2155 ZPT = col_getPressure(columnTrlOnTrlLev,ilyr ,headerIndex,varLevel)
2156 ZPB = col_getPressure(columnTrlOnTrlLev,ilyr+1,headerIndex,varLevel)
2157 ZWB = LOG(zlev/ZPT)/LOG(ZPB/ZPT)
2158 ZWT = 1. - ZWB
2159
2160 ZMOD = ZWB*col_ptr_uv(ilyr+1) + ZWT*col_ptr_uv(ilyr)
2161
2162 TO_WGH = 0.0D0
2163 TO_DSP = 0.0D0
2164 if(zlev > 20000. .and. zlev < 30000. .and. passe_once .and. print_debug) then
2165 write(*,'(3a10,2i10,4f12.3)') 'stlchka',cstnid,varName,iqiv,imet,zlev,ZMOD,ZVAR,E_DRIFT
2166 end if
2167 do jlev = 2, col_getNumLev(columnTrlOnTrlLev,'MM') - 1
2168 ZPT= col_getPressure(columnTrlOnTrlLev,jlev-1,headerIndex,varLevel)
2169 ZPC= col_getPressure(columnTrlOnTrlLev,jlev ,headerIndex,varLevel)
2170 ZPB= col_getPressure(columnTrlOnTrlLev,jlev+1,headerIndex,varLevel)
2171 ZOTR = col_ptr_uv(jlev)
2172 SP_WGH = exp(-0.5*((ZPC - zlev)**2)/(E_HEIGHT**2))*((ZPB - ZPT)/2)
2173 TO_DSP = TO_DSP + SP_WGH*((ZOTR - ZMOD)**2)
2174 TO_WGH = TO_WGH + SP_WGH
2175 if(zlev > 20000. .and. zlev < 30000. .and. passe_once .and. print_debug) then
2176 write(*,'(a10,i10,4f12.3)') 'stlchk',jlev,ZPT,ZPC,ZPB,ZOTR
2177 end if
2178
2179 end do
2180
2181 E_VHGT = sqrt(TO_DSP/TO_WGH)
2182 zoer = sqrt(E_VHGT**2 + E_DRIFT**2)
2183
2184 if(zlev > 20000. .and. zlev < 30000. .and. passe_once .and. print_debug) then
2185 write(*,'(a10,4f10.2)') 'stlchkb',zoer,E_VHGT,E_DRIFT,E_HEIGHT
2186 passe_once = .false.
2187 end if
2188
2189 call obs_bodySet_r(obsSpaceData,OBS_OER,bodyIndex,zoer)
2190
2191 if(print_debug) write(*,'(2a10,6f12.3,4i10)') 'hgterr',cstnid,zlat,zlon,zlev/100.,E_HEIGHT/100.0,E_DRIFT,zoer,imet,itrn,ihav,J_SAT
2192
2193 end do BODY
2194
2195 end subroutine oer_sw
2196
2197 !--------------------------------------------------------------------------
2198 ! get_height_error
2199 !--------------------------------------------------------------------------
2200 subroutine get_height_error(stnid, methode, terrain, htasmet, zlat, zlev, E_HEIGHT, J_SAT)
2201 implicit none
2202
2203 ! Arguments:
2204 character(len=9), intent(in) :: stnid
2205 integer, intent(in) :: methode
2206 integer, intent(in) :: terrain
2207 integer, intent(in) :: htasmet
2208 integer, intent(out) :: J_SAT
2209 real(8), intent(in) :: zlat
2210 real(8), intent(in) :: zlev
2211 real(8), intent(out) :: E_HEIGHT
2212
2213 ! Locals:
2214 integer, parameter :: NLVLVALUE=9
2215 integer :: I_HGT, jelm, jlev, jcat, i_typ, hemisphere
2216 real(8) :: zlev_hpa, ZWB, ZWT
2217 logical :: interpole
2218
2219 if(zlat >= 0) hemisphere = 1
2220 if(zlat < 0) hemisphere = 2
2221
2222
2223 J_SAT = 1 ! default value
2224
2225 i_typ = 0
2226 list1: do jlev = 1,n_sat_type
2227 do jelm = 2, 10
2228 if(trim(stnid(2:9)) == trim(SAT_AMV(jlev,jelm))) then
2229 i_typ = jlev
2230 exit list1
2231 end if
2232 end do
2233 end do list1
2234
2235 if (i_typ /= 0) then
2236 list2: do jcat = 1,n_categorie
2237 if(trim(SAT_AMV(i_typ,1)) == trim(SAT_LIST(jcat))) then
2238 if(tbl_m(jcat) == 0 .or. tbl_m(jcat) == methode) then
2239 if(tbl_h(jcat) == 0 .or. tbl_h(jcat) == htasmet) then
2240 if(tbl_t(jcat) == 0 .or. tbl_t(jcat) == terrain+1) then
2241 if(tbl_g(jcat) == 0 .or. tbl_g(jcat) == hemisphere) then
2242 J_SAT = jcat
2243 exit list2
2244 end if
2245 end if
2246 end if
2247 end if
2248 end if
2249 end do list2
2250 end if
2251
2252 zlev_hpa = zlev/100.
2253 interpole = .true.
2254 do I_HGT=1,NLVLVALUE
2255 if(zlev_hpa >= LVLVALUE(I_HGT)) exit
2256 end do
2257 if(I_HGT == 1) interpole = .false.
2258 if(I_HGT == NLVLVALUE + 1) interpole = .false.
2259 if(I_HGT == NLVLVALUE + 1) I_HGT = NLVLVALUE
2260
2261 if(interpole) then
2262
2263 ZWB = LOG(zlev_hpa/LVLVALUE(I_HGT-1))/LOG(LVLVALUE(I_HGT)/LVLVALUE(I_HGT-1))
2264 ZWT = 1. - ZWB
2265 E_HEIGHT = ZWT*HGT_ERR(J_SAT,I_HGT-1) + ZWB*HGT_ERR(J_SAT,I_HGT)
2266
2267 else
2268
2269 E_HEIGHT = HGT_ERR(J_SAT,I_HGT)
2270
2271 end if
2272
2273!
2274! Retourne l'erreur en (Pa)
2275!
2276 E_HEIGHT = E_HEIGHT*100.0
2277
2278 end subroutine get_height_error
2279
2280 !--------------------------------------------------------------------------
2281 ! oer_SETERRGPSRO
2282 !--------------------------------------------------------------------------
2283 SUBROUTINE oer_SETERRGPSRO(columnTrlOnTrlLev, obsSpaceData, beSilent)
2284 !
2285 !:Purpose: Compute estimated errors for GPSRO observations
2286 !
2287 IMPLICIT NONE
2288
2289 ! Arguments:
2290 type(struct_columnData), intent(in) :: columnTrlOnTrlLev
2291 type(struct_obs), intent(inout) :: obsSpaceData
2292 logical, intent(in) :: beSilent
2293
2294 ! Locals:
2295 integer headerIndex, IDATYP, bodyIndex, iProfile, varNum
2296 REAL*8 zLat, Lat, sLat
2297 REAL*8 zLon, Lon
2298 REAL*8 zAzm
2299 REAL*8, allocatable :: ZPP(:)
2300 REAL*8, allocatable :: ZTT(:)
2301 REAL*8, allocatable :: ZHU(:)
2302 REAL*8, allocatable :: zHeight(:)
2303 REAL*8, allocatable :: ZUU(:)
2304 REAL*8, allocatable :: ZVV(:)
2305 REAL*8 DH,DDH
2306 integer JL, isat, JH, NGPSLEV, NWNDLEV
2307 REAL*8 zMT, Rad, Geo, zP0
2308 REAL*8 HNH1, HJH, SUM0, SUM1, ZMIN, WFGPS, H1, F2, F3, F4
2309 LOGICAL ASSIM, L1, L2, L3
2310 integer NH, NH1
2311 TYPE(GPS_PROFILE) :: PRF
2312 REAL*8 , allocatable :: H (:),AZMV(:)
2313 REAL*8 , allocatable :: ZOBS(:),ZREF(:),ZOFF(:),ZERR(:), ZMHX(:)
2314 TYPE(GPS_DIFF), allocatable :: RSTV(:)
2315
2316 if (.not. beSilent) write(*,*)'ENTER SETERRGPSRO'
2317 !
2318 ! * 1. Initializations
2319 ! * ---------------
2320 !
2321 NGPSLEV=col_getNumLev(columnTrlOnTrlLev,'TH')
2322 NWNDLEV=col_getNumLev(columnTrlOnTrlLev,'MM')
2323 allocate(ZPP (NGPSLEV))
2324 allocate(ZTT (NGPSLEV))
2325 allocate(ZHU (NGPSLEV))
2326 allocate(zHeight (NGPSLEV))
2327 allocate(ZUU (NGPSLEV))
2328 allocate(ZVV (NGPSLEV))
2329 !
2330 allocate(H (gps_RO_MAXPRFSIZE))
2331 allocate(AZMV (gps_RO_MAXPRFSIZE))
2332 allocate(ZOBS (gps_RO_MAXPRFSIZE))
2333 allocate(ZREF (gps_RO_MAXPRFSIZE))
2334 allocate(ZOFF (gps_RO_MAXPRFSIZE))
2335 allocate(ZERR (gps_RO_MAXPRFSIZE))
2336 allocate(RSTV (gps_RO_MAXPRFSIZE))
2337 allocate(ZMHX (gps_RO_MAXPRFSIZE))
2338 !
2339 ! Loop over all header indices of the 'RO' family:
2340 !
2341 call obs_set_current_header_list(obsSpaceData,'RO')
2342 HEADER: do
2343 headerIndex = obs_getHeaderIndex(obsSpaceData)
2344 if (headerIndex < 0) exit HEADER
2345 !
2346 ! * Process only refractivity data (codtyp 169)
2347 !
2348 IDATYP = obs_headElem_i(obsSpaceData,OBS_ITY,headerIndex)
2349 IF (IDATYP == 169) then
2350 !
2351 ! * Scan for requested data values of the profile, and count them
2352 !
2353 ASSIM = .FALSE.
2354 NH = 0
2355 call obs_set_current_body_list(obsSpaceData, headerIndex)
2356 BODY: do
2357 bodyIndex = obs_getBodyIndex(obsSpaceData)
2358 if (bodyIndex < 0) exit BODY
2359 IF (obs_bodyElem_i(obsSpaceData, OBS_ASS, bodyIndex) == obs_assimilated) then
2360 ASSIM = .TRUE.
2361 NH = NH + 1
2362 end if
2363 end do BODY
2364 !
2365 ! * If assimilations are requested, prepare and apply the observation operator
2366 !
2367 IF (ASSIM) then
2368 iProfile=gps_iprofile_from_index(headerIndex)
2369 varNum = gps_vRO_IndexPrf(iProfile, 2)
2370 !
2371 ! * Basic geometric variables of the profile:
2372 !
2373 isat = obs_headElem_i(obsSpaceData,OBS_SAT,headerIndex)
2374 Rad = obs_headElem_r(obsSpaceData,OBS_TRAD,headerIndex)
2375 Geo = obs_headElem_r(obsSpaceData,OBS_GEOI,headerIndex)
2376 zAzm = obs_headElem_r(obsSpaceData,OBS_AZA,headerIndex) / MPC_DEGREES_PER_RADIAN_R8
2377 zMT = col_getHeight(columnTrlOnTrlLev,0,headerIndex,'SF')
2378 !
2379 ! * Profile at the observation location:
2380 !
2381 zLat = obs_headElem_r(obsSpaceData,OBS_LAT,headerIndex)
2382 zLon = obs_headElem_r(obsSpaceData,OBS_LON,headerIndex)
2383 Lat = zLat * MPC_DEGREES_PER_RADIAN_R8
2384 Lon = zLon * MPC_DEGREES_PER_RADIAN_R8
2385 sLat = sin(zLat)
2386 zMT = zMT * ec_rg / gps_gravitysrf(sLat)
2387 zP0 = col_getElem(columnTrlOnTrlLev,1,headerIndex,'P0')
2388 DO JL = 1, NGPSLEV
2389 !
2390 ! * Profile x
2391 !
2392 ZPP(JL) = col_getPressure(columnTrlOnTrlLev,JL,headerIndex,'TH')
2393 ZTT(JL) = col_getElem(columnTrlOnTrlLev,JL,headerIndex,'TT') - MPC_K_C_DEGREE_OFFSET_R8
2394 ZHU(JL) = col_getElem(columnTrlOnTrlLev,JL,headerIndex,'HU')
2395 ZUU(JL) = 0.d0
2396 ZVV(JL) = 0.d0
2397 zHeight(jl) = col_getHeight(columnTrlOnTrlLev,jl,headerIndex,'TH')
2398 end do
2399
2400 if((col_getPressure(columnTrlOnTrlLev,1,headerIndex,'TH') + 1.0d-4) < &
2401 col_getPressure(columnTrlOnTrlLev,1,headerIndex,'MM')) then
2402 ! case with top thermo level above top momentum level (Vcode=5002)
2403 do jl = 1, nwndlev
2404 zuu(jl) = col_getElem(columnTrlOnTrlLev,jl, headerIndex, 'UU')
2405 zvv(jl) = col_getElem(columnTrlOnTrlLev,jl, headerIndex, 'VV')
2406 end do
2407 else
2408 ! case without top thermo above top momentum level or unstaggered (Vcode=5001/4/5)
2409 do jl = 1, nwndlev - 1
2410 zuu(jl) = col_getElem(columnTrlOnTrlLev, jl + 1, headerIndex, 'UU')
2411 zvv(jl) = col_getElem(columnTrlOnTrlLev, jl + 1, headerIndex, 'VV')
2412 end do
2413 zuu(nwndlev) = zuu(nwndlev - 1)
2414 zvv(nwndlev) = zuu(nwndlev - 1)
2415 end if
2416 zuu(ngpslev) = zuu(nwndlev)
2417 zvv(ngpslev) = zuu(nwndlev)
2418 !
2419 ! * GPS profile structure:
2420 !
2421 call gps_struct1sw_v2(ngpslev,zLat,zLon,zAzm,zMT,Rad,geo,zP0,zPP,zTT,zHU,zUU,zVV,zHeight,prf)
2422 !
2423 ! * Prepare the vector of all the observations:
2424 !
2425 NH1 = 0
2426 !
2427 ! * Loop over all body indices for this index_header:
2428 ! * (start at the beginning of the list)
2429 !
2430 call obs_set_current_body_list(obsSpaceData, headerIndex)
2431 BODY_2: do
2432 bodyIndex = obs_getBodyIndex(obsSpaceData)
2433 if (bodyIndex < 0) exit BODY_2
2434 IF (obs_bodyElem_i(obsSpaceData, OBS_ASS, bodyIndex) == obs_assimilated) then
2435 NH1 = NH1 + 1
2436 H(NH1) = obs_bodyElem_r(obsSpaceData, OBS_PPP, bodyIndex)
2437 AZMV(NH1)= zAzm
2438 ZOBS(NH1)= obs_bodyElem_r(obsSpaceData, OBS_VAR, bodyIndex)
2439 ! * Reference value:
2440 IF (varNum == bufr_nebd) then
2441 ZREF(NH1) = 0.025d0*exp(-(H(NH1)-Rad)/6500.d0)
2442 ELSE
2443 ZREF(NH1) = 300.d0*exp(-H(NH1)/6500.d0)
2444 end if
2445 end if
2446 end do BODY_2
2447 !
2448 ! * Apply the observation operator:
2449 !
2450 IF (varNum == bufr_nebd) then
2451 call GPS_BNDOPV1(H, AZMV, NH, PRF, RSTV)
2452 ELSE
2453 call GPS_REFOPV (H, NH, PRF, RSTV)
2454 end if
2455 !
2456 ! * Perform the (H(x)-Y)/R operation:
2457 !
2458 DO NH1 = 1, NH
2459 ZMHX(NH1) = RSTV(NH1)%VAR
2460 !
2461 ! * Normalized offset:
2462 !
2463 IF (.NOT.gps_roBNorm) then
2464 ZOFF(NH1) = (ZOBS(NH1) - ZMHX(NH1)) / ZREF(NH1)
2465 ELSE
2466 ZOFF(NH1) = (ZOBS(NH1) - ZMHX(NH1)) / ZMHX(NH1)
2467 end if
2468 end do
2469 !
2470 ! * The procedure below is well tested to collectively
2471 ! * create error profiles from the data profile, and
2472 ! * intended to be used for these data.
2473 !
2474 DH = 5000.d0
2475 if (varNum == bufr_nebd) then
2476 ZMIN=0.01D0
2477 else
2478 ZMIN=0.002D0
2479 end if
2480
2481 if (varNum == bufr_nerf) then
2482 if (trim(gps_roError) == 'DYNAMIC') then
2483 do NH1 = 1, NH
2484 SUM0=1.d-30
2485 SUM1=0.d0
2486 do JH = 1, NH
2487 if (H(JH) <= gps_HtpMaxEr) then
2488 DDH=H(JH)-H(NH1)
2489 SUM0=SUM0+EXP(-(DDH/DH)**2)
2490 SUM1=SUM1+EXP(-(DDH/DH)**2)*ZOFF(JH)**2
2491 end if
2492 end do
2493 ZERR(NH1)=(SUM1/SUM0)**0.5D0
2494 if (ZERR(NH1) < ZMIN) ZERR(NH1) = ZMIN
2495 end do
2496 else if (trim(gps_roError) == 'STATIC_2018') then
2497 ! this was introduced by Maziar in late 2018 on advice by Josep
2498 do NH1 = 1, NH
2499 HNH1 = H(NH1)
2500 ZERR(NH1) = 0.05d0
2501 L1=(HNH1 <= 10000.d0)
2502 L2=(HNH1 > 10000.d0 .and. HNH1 < 30000.d0)
2503 L3=(HNH1 > 30000.d0)
2504 IF (L1) ZERR(NH1)=0.005d0+0.020d0*(10000.d0-HNH1)/10000.d0
2505 IF (L2) ZERR(NH1)=0.005d0
2506 IF (L3) ZERR(NH1)=0.005d0+0.030d0*(HNH1-30000.d0)/30000.d0
2507 if (ZERR(NH1) < ZMIN) ZERR(NH1) = ZMIN
2508 end do
2509 else if (trim(gps_roError) == 'STATIC_2014') then
2510 ! recipe used in EnKF from Josep by email on February 25 2014
2511 do NH1 = 1, NH
2512 HNH1 = H(NH1)
2513 select case (nint(hnh1))
2514 case(:10000)
2515 ZERR(NH1) = 0.005 + 0.015 * ((10000.0-hnh1)/10000.0)
2516 case (10001:30000)
2517 ZERR(NH1) = 0.005
2518 case (30001:)
2519 ZERR(NH1) = 0.005 + 0.010 * ((hnh1-30000.0)/10000.0)
2520 end select
2521 end do
2522 else
2523 call utl_abort('oer_setErrGPSro: Invalid value for gps_roError')
2524 endif
2525 else
2526 if (trim(gps_roError) == 'DYNAMIC') then
2527 ZMIN = 0.01D0
2528 do NH1 = 1, NH
2529 HNH1=H(NH1)-Rad
2530 SUM0=1.d-30
2531 SUM1=0.d0
2532 DH = 1000.d0 + 0.1d0 * HNH1
2533 do JH = 1, NH
2534 HJH=H(JH)-Rad
2535 if (HJH <= gps_HtpMaxEr) then
2536 DDH=HJH-HNH1
2537 SUM0=SUM0+EXP(-(DDH/DH)**2+(DDH/DH))
2538 SUM1=SUM1+EXP(-(DDH/DH)**2+(DDH/DH))*ZOFF(JH)**2
2539 end if
2540 end do
2541 ZERR(NH1)=(SUM1/SUM0)**0.5D0
2542 if (ZERR(NH1) < ZMIN) ZERR(NH1) = ZMIN
2543 if (H(NH1) < PRF%ast(ngpslev)%Var) ZERR(NH1)=0.08
2544 end do
2545 else if (trim(gps_roError) == 'STATIC_2018') then
2546 do NH1 = 1, NH
2547 ZERR(NH1)=0.05d0
2548 HNH1=H(NH1)-Rad
2549 L1=(HNH1 <= 10000.d0)
2550 L2=(HNH1 > 10000.d0 .and. HNH1 < 30000.d0)
2551 L3=(HNH1 > 30000.d0)
2552 IF (L1) ZERR(NH1)=0.02d0+0.08d0*(10000.d0-HNH1)/10000.d0
2553 IF (L2) ZERR(NH1)=0.02d0
2554 IF (L3) ZERR(NH1)=0.02d0+0.13d0*(HNH1-30000.d0)/30000.d0
2555 IF (isat==3 .or. isat==4 .or. isat==5) ZERR(NH1) = 2*ZERR(NH1)
2556 IF (ZERR(NH1) < ZMIN) ZERR(NH1) = ZMIN
2557 end do
2558 end if
2559 end if
2560
2561 NH1 = 0
2562 !
2563 ! * Loop over all body indices for this index_header:
2564 ! * (start at the beginning of the list)
2565 !
2566 call obs_set_current_body_list(obsSpaceData, headerIndex)
2567 BODY_4: do
2568 bodyIndex = obs_getBodyIndex(obsSpaceData)
2569 if (bodyIndex < 0) exit BODY_4
2570 IF (obs_bodyElem_i(obsSpaceData, OBS_ASS, bodyIndex) == obs_assimilated) then
2571 NH1 = NH1 + 1
2572 H1 = H(NH1)
2573 if (varNum == bufr_nebd) then
2574 H1 = H1 - Rad
2575 F2 = 0.5d0*(erf((H1-22000.d0)/5000.d0)+1.d0)
2576 else
2577 F2 = exp(0.5d0*H1/6500.d0)
2578 endif
2579 F3 = exp(-((H1- 6500.d0)/6500.d0)**2)
2580 F4 = exp(-((H1-14500.d0)/6500.d0)**2)
2581 WFGPS = gps_WGPS(ISAT,1) + gps_WGPS(ISAT,2) * F2 + gps_WGPS(ISAT,3) * F3 + gps_WGPS(ISAT,4) * F4
2582 IF (WFGPS < 1.D0) WFGPS = 1.D0
2583 WFGPS = sqrt(WFGPS)
2584 !
2585 ! * Observation error S
2586 !
2587 if (.NOT.gps_roBNorm) then
2588 call obs_bodySet_r(obsSpaceData, OBS_OER, bodyIndex, ZERR(NH1)*ZREF(NH1)*WFGPS)
2589 else
2590 call obs_bodySet_r(obsSpaceData, OBS_OER, bodyIndex, ZERR(NH1)*ZMHX(NH1)*WFGPS)
2591 end if
2592 end if
2593 end do BODY_4
2594 end if
2595 end if
2596 end do HEADER
2597
2598 deallocate(RSTV)
2599 deallocate(ZERR)
2600 deallocate(ZOFF)
2601 deallocate(ZREF)
2602 deallocate(ZOBS)
2603 deallocate(AZMV)
2604 deallocate(H)
2605 deallocate(ZMHX)
2606
2607 deallocate(zVV)
2608 deallocate(zUU)
2609 deallocate(zHU)
2610 deallocate(zHeight)
2611 deallocate(zTT)
2612 deallocate(zPP)
2613
2614 if (.not.beSilent) write(*,*)'oer_setErrGPSRO: done'
2615
2616 END SUBROUTINE OER_SETERRGPSRO
2617
2618 !--------------------------------------------------------------------------
2619 ! oer_SETERRGPSGB
2620 !--------------------------------------------------------------------------
2621 SUBROUTINE oer_SETERRGPSGB(columnTrlOnTrlLev, obsSpaceData, beSilent, ldata, analysisMode)
2622 !
2623 !:Purpose:
2624 !
2625 ! - Set the observation errors [OBS_OER] and Std(O-P) [OBS_HPHT] for GB-GPS ZTD data.
2626 ! (GPS surface met obs errors are set before in observation_erreurs_mod.ftn90.
2627 ! The ZTD error is also initialized to the "formal error" or to 1.0 if missing.)
2628 ! - Returns ldata=.false. if there are no GPS ZTD data to assimilate
2629 ! and also sets the modgpsztd_mod variable gps_gb_numZTD = 0.
2630 !
2631 IMPLICIT NONE
2632 !!
2633 !! NOTE: gps_gb_YZDERRWGT IN modgpsztd_mod (FROM NML FILE) IS USED FOR ERROR WEIGHTING
2634 !! OF TIME SERIES (FGAT) GPS ZTD OBSERVATIONS TO ACCOUNT FOR TEMPORAL ERROR
2635 !! CORRELATIONS.
2636 !!
2637
2638 ! Arguments:
2639 type(struct_columnData), intent(in) :: columnTrlOnTrlLev
2640 type(struct_obs) , intent(inout) :: obsSpaceData
2641 logical , intent(in) :: beSilent
2642 logical , intent(out) :: ldata
2643 logical , intent(in) :: analysisMode
2644
2645 ! Locals:
2646 integer bodyIndex, headerIndex, ityp, iass, IZTDJ, NBRPDATE, ICOUNT, ICOUNT2
2647 integer nlev_T
2648 LOGICAL LLCZTDE, LLFER, LLFZTDE, LLZTD, LLRZTDE, ASSIM, ERRSET, DEBUG, LESTP
2649 LOGICAL LLZWD
2650 character(len=12) :: cstnid
2651 REAL*8 ZTDERR, ZZTD, ZMINZDE, ZPSFC, ZHD, ZWD, ZTDOER, zlev, zval, ZZWD
2652 REAL*8 ZBTSFC, ZBPSFC, ZBZSFC, ZDZ, ZSTDOMP
2653 !
2654 ! ZZDERMIN = MIN ZTD OER VALUE (M)
2655 ! ZZDERMAX = MAX ZTD OER VALUE (M)
2656 ! ZTDERFAC = MULTIPLICATION FACTOR FOR FORMAL ZTD MEASUREMENT ERROR
2657 ! ZOPEFAC = FRACTION OF REGRESSION EQUATION SD(O-P) TO USE AS ZTD OBSERVATION ERROR
2658 ! ----------------------------------------------------------------------------------
2659 !
2660 REAL*8 ZZDERMIN, ZZDERMAX, ZTDERFAC, ZOPEFAC
2661 DATA ZZDERMIN /0.004D0/
2662 DATA ZZDERMAX /0.030D0/
2663 DATA ZTDERFAC /3.0D0/
2664 DATA ZOPEFAC /1.0D0/
2665 !
2666 ! FOR ESTIMATION OF PSFC (IF MISSING)
2667 ! ZGAMMA = (NEG. OF) TEMPERATURE LAPSE RATE (K/M)
2668 !
2669 REAL*8 ZGAMMA
2670 DATA ZGAMMA /0.0065D0/
2671
2672 ! ----------------------------------------------------------------------------------
2673 ! LINEAR REGRESSION EQUATION CONSTANTS AND COEFFS FOR ZTD ERROR AND STD(O-P):
2674
2675 ! ZRCONST, ZRCOEFF:
2676 ! - From linear regression of Desroziers error estimates binned by observed ZWD.
2677 ! - Gives ZTDerror (mm) as function of ZWD (m):
2678 ! ZTDerror(mm) = ZRCONST + ZRCOEFF*ZWD(m)
2679 ! ZRCONST2, ZRCOEFF2:
2680 ! - From linear regression of Std(O-P) binned by observed ZWD.
2681 ! - Gives Std(O-P) (mm) as function of ZWD (m):
2682 ! Std(O-P)(mm) = ZRCONST2 + ZRCOEFF2*ZWD(m)
2683 ! ----------------------------------------------------------------------------------
2684 REAL*8 ZRCONST, ZRCOEFF, ZRCONST2, ZRCOEFF2
2685 DATA ZRCONST /5.12D0/
2686 DATA ZRCOEFF /26.4D0/
2687 DATA ZRCONST2 /6.67D0/
2688 DATA ZRCOEFF2 /42.6D0/
2689
2690 if (.not.beSilent) write(*,*) 'ENTER SETERRGPSGB'
2691 !
2692 DEBUG = .FALSE.
2693
2694 LLCZTDE = .FALSE.
2695 LLRZTDE = .FALSE.
2696 LLFZTDE = .FALSE.
2697 IF (gps_gb_YZTDERR < 0.0D0) then
2698 LLFZTDE = .TRUE.
2699 else if (gps_gb_YZTDERR > 0.0D0) then
2700 LLCZTDE = .TRUE.
2701 ELSE
2702 LLRZTDE = .TRUE.
2703 end if
2704
2705 nlev_T = col_getNumLev(columnTrlOnTrlLev,'TH')
2706
2707 ldata = .false.
2708 ICOUNT = 0
2709 ICOUNT2 = 0
2710 !
2711 ! Loop over all header indices of the 'GP' family:
2712 !
2713 call obs_set_current_header_list(obsSpaceData,'GP')
2714 HEADER: do
2715 headerIndex = obs_getHeaderIndex(obsSpaceData)
2716 if (headerIndex < 0) exit HEADER
2717 NBRPDATE = obs_headElem_i(obsSpaceData, OBS_DAT, headerIndex)
2718 LLZTD = .FALSE.
2719 LLFER = .FALSE.
2720 LLZWD = .FALSE.
2721 ASSIM = .FALSE.
2722 ERRSET = .FALSE.
2723 ZZTD = -1.0D0
2724 ZZWD = -1.0D0
2725 ZPSFC = -1.0D0
2726 LESTP = .FALSE.
2727 ZSTDOMP = 15.0D0*0.001D0
2728
2729 ! Get Psfc (Pa), Tsfc (K) and model surface height (m) from background profile
2730
2731 ZBPSFC = col_getElem(columnTrlOnTrlLev, 1, headerIndex, 'P0')
2732 ZBTSFC = col_getElem(columnTrlOnTrlLev, nlev_T, headerIndex, 'TT')
2733 ZBZSFC = col_getHeight(columnTrlOnTrlLev, nlev_T, headerIndex, 'TH')
2734 !
2735 ! Loop over all body indices of current report; Set the ZTD error if
2736 ! constant value specified (LLCZTDE=true). Get GPS height and Psfc obs (if any).
2737 ! Get ZTD obs, ZTD formal error and ZWD observation.
2738 !
2739 call obs_set_current_body_list(obsSpaceData, headerIndex)
2740 BODY: do
2741 bodyIndex = obs_getBodyIndex(obsSpaceData)
2742 if (bodyIndex < 0) exit BODY
2743 ityp = obs_bodyElem_i(obsSpaceData, OBS_VNM, bodyIndex)
2744 iass = obs_bodyElem_i(obsSpaceData, OBS_ASS, bodyIndex)
2745 zval = obs_bodyElem_r(obsSpaceData, OBS_VAR, bodyIndex)
2746 ! Store Psfc
2747 IF (ityp == BUFR_NEPS) then
2748 IF (zval > 0.0D0) ZPSFC = zval
2749 end if
2750 ! Set ZTDOER to constant value (if LLCZTDE); get value of ZTD,
2751 ! ZTD formal error (OBS_OER) and antenna height (OBS_PPP).
2752 IF (ityp == BUFR_NEZD) then
2753 IF (LLCZTDE) then
2754 ZTDOER = gps_gb_YZTDERR
2755 ERRSET = .TRUE.
2756 end if
2757 zlev = obs_bodyElem_r(obsSpaceData, OBS_PPP, bodyIndex)
2758 ZTDERR = obs_bodyElem_r(obsSpaceData, OBS_OER, bodyIndex)
2759 IF (ZTDERR /= 1.0D0) LLFER = .TRUE.
2760 IZTDJ = bodyIndex
2761 IF (zval > 0.0D0) then
2762 ZZTD = zval
2763 LLZTD = .TRUE.
2764 end if
2765 IF (iass == obs_assimilated) ASSIM = .TRUE.
2766 end if
2767 IF (ityp == BUFR_NEZW) then
2768 IF (zval > 0.0D0) then
2769 ZZWD = zval
2770 LLZWD = .TRUE.
2771 end if
2772 end if
2773 end do BODY
2774
2775 ! Initialize Std(O-P) to 15 mm for ZTD observation (for bgck mode)
2776 IF (LLZTD .and. .NOT.analysisMode) &
2777 call obs_bodySet_r(obsSpaceData,OBS_HPHT,IZTDJ,ZSTDOMP)
2778
2779 ! Replace formal ZTD error with real error for all ZTD to be assimilated.
2780 ! Set Std(O-P) as function of ZWD for ZTD observation and store in OBS_HPHT.
2781
2782 if (ASSIM) then
2783 if (LLZTD) then
2784 ldata = .true.
2785 ICOUNT = ICOUNT + 1
2786 if (LLZWD) then
2787 ZWD = ZZWD
2788 else
2789 ! If Psfc obs is missing, estimate the pressure from background state
2790 if (ZPSFC < 0.0D0) then
2791 LESTP = .TRUE.
2792 ZDZ = zlev - ZBZSFC
2793 ZPSFC = ZBPSFC * &
2794 (1.0D0-(ZGAMMA/ZBTSFC)*ZDZ)**(ec_rg/(MPC_RGAS_DRY_AIR_R8*ZGAMMA))
2795 ICOUNT2 = ICOUNT2 + 1
2796 end if
2797 ! Compute the hydrostatic delay ZHD (m) from Psfc (Pa)
2798 ZHD = 2.2766D-05 * ZPSFC
2799 ! Compute the wet delay (m) from ZTD and ZHD. Avoid negative ZWD.
2800 if (ZHD > ZZTD) then
2801 ZWD = 0.0D0
2802 else
2803 ZWD = ZZTD - ZHD
2804 end if
2805 end if
2806 ! Std(O-P) for background check. Limit to 30 mm in case ZTD obs is bad (too high).
2807 ZSTDOMP = (ZRCONST2 + ZRCOEFF2*ZWD)*0.001D0
2808 ZSTDOMP = MIN(ZZDERMAX, ZSTDOMP)
2809 ! Compute ZTD error as a function of ZWD using regression coeff (SD(O-P) vs ZWD).
2810 ! Take fraction ZOPEFAC of computed error and convert from mm to m.
2811 ! Ensure error is > ZZDERMIN and < ZZDERMAX
2812 IF (.NOT. ERRSET) then
2813 ZMINZDE = ZRCONST + ZRCOEFF*ZWD
2814 ZMINZDE = ZMINZDE * ZOPEFAC * 0.001D0
2815 IF (LLRZTDE) then
2816 ZTDOER = MAX(ZZDERMIN, ZMINZDE)
2817 ZTDOER = MIN(ZZDERMAX, ZTDOER)
2818 ELSE
2819 IF (LLFER) then
2820 ZTDOER = MAX(ZZDERMIN, ZTDERR*ZTDERFAC)
2821 ELSE
2822 ZTDOER = MAX(ZZDERMIN, ZMINZDE)
2823 ZTDOER = MIN(ZZDERMAX, ZTDOER)
2824 end if
2825 end if
2826 ! Ensure that error is not less than formal error ZTDERR
2827 IF (LLFER) then
2828 IF (DEBUG .and. ICOUNT <= 50) then
2829 write(*,*) cstnid, &
2830 ' FORMAL ERR, OBS ERROR (mm) = ', &
2831 ZTDERR*1000.D0, ZTDOER*1000.D0
2832 end if
2833 ZTDOER = MAX(ZTDOER, ZTDERR)
2834 end if
2835 end if
2836 ! *** APPLY TIME-SERIES WEIGHTING FACTOR TO OBSERVATION ERROR (gps_gb_YZDERRWGT=1 FOR 3D THINNING)
2837 call obs_bodySet_r(obsSpaceData,OBS_OER,IZTDJ, ZTDOER*gps_gb_YZDERRWGT)
2838 IF (.NOT.analysisMode) call obs_bodySet_r(obsSpaceData,OBS_HPHT,IZTDJ, ZSTDOMP)
2839 IF (DEBUG .and. (ICOUNT2 <= 50) .and. LESTP) then
2840 write(*,*) 'TAG SITE ZTD ERROR ELEV PSFC ZWD STDOMP'
2841 write(*,*) 'ERRDEBUG ', cstnid, &
2842 ZZTD*1000.D0, ZTDOER*1000.D0, zlev, ZPSFC/100.D0, ZWD*1000.D0, ZSTDOMP*1000.D0
2843 end if
2844 ELSE
2845 call utl_abort('SETERRGPSGB: ERROR:NEGATIVE ZTD VALUE!')
2846 end if
2847 end if
2848
2849 !
2850 end do HEADER
2851
2852 ! IF (DEBUG) call utl_abort('******DEBUG STOP*******')
2853
2854 IF (.not.ldata) gps_gb_numZTD = 0
2855
2856 IF (ldata .and. .not.beSilent) write(*,*) ' gps_gb_numZTD = ', ICOUNT
2857
2858 if (.not.beSilent) write(*,*) 'EXIT SETERRGPSGB'
2859
2860 END SUBROUTINE OER_SETERRGPSGB
2861
2862 !--------------------------------------------------------------------------
2863 ! oer_setErrBackScatAnisIce
2864 !--------------------------------------------------------------------------
2865 subroutine oer_setErrBackScatAnisIce(obsSpaceData, beSilent, columnTrlOnTrlLev_opt)
2866 !
2867 !:Purpose: Compute estimated errors for ASCAT backscatter anisotropy observations
2868 !
2869 implicit none
2870
2871 ! Arguments:
2872 type(struct_obs), intent(inout) :: obsSpaceData
2873 logical, intent(in) :: beSilent
2874 type(struct_columnData), optional, intent(in) :: columnTrlOnTrlLev_opt
2875
2876 ! Locals:
2877 type(struct_columnData) :: columnTrlOnTrlLev
2878 type(struct_hco), pointer :: hco_trl => null()
2879 type(struct_vco), pointer :: vco_trl => null()
2880 type(struct_gsv) :: stateVector
2881 integer :: headerIndex, bodyIndex
2882 integer :: idate, imonth, varno
2883 integer :: trackCellNum
2884 real(8) :: conc, obsErrStdDev
2885 character(len=*), parameter :: myName = 'oer_setErrBackScatAnisIce'
2886 character(len=8) :: ccyymmdd
2887
2888 if (.not. beSilent) write(*,*) 'ENTER '//myName
2889
2890 if (.not.present(columnTrlOnTrlLev_opt)) then
2891 ! Need background ice concentration
2892 call hco_SetupFromFile( hco_trl, './bgSeaIceConc', '', varName_opt='GL' )
2893 call vco_SetupFromFile( vco_trl, './bgSeaIceConc')
2894 call gsv_allocate( stateVector, 1, hco_trl, vco_trl, dateStamp_opt=-1, &
2895 mpi_local_opt=.true., mpi_distribution_opt='Tiles', &
2896 dataKind_opt=4, hInterpolateDegree_opt='LINEAR', &
2897 varNames_opt=(/'GL'/) )
2898 call gsv_zero( stateVector )
2899 call gio_readFromFile( stateVector, './bgSeaIceConc', '', 'P@', &
2900 unitConversion_opt = .false. )
2901 call col_setVco(columnTrlOnTrlLev, vco_trl)
2902 call col_allocate(columnTrlOnTrlLev, obs_numHeader(obsSpaceData), mpiLocal_opt=.true.)
2903 call s2c_nl( stateVector, obsSpaceData, columnTrlOnTrlLev, hco_trl, &
2904 timeInterpType='NEAREST' )
2905 call gsv_deallocate( stateVector )
2906 end if
2907
2908 !
2909 ! Loop over all header indices of the 'GL' family:
2910 !
2911 call obs_set_current_header_list(obsSpaceData,'GL')
2912 HEADER: do
2913 headerIndex = obs_getHeaderIndex(obsSpaceData)
2914 if (headerIndex < 0) exit HEADER
2915 call obs_set_current_body_list(obsSpaceData, headerIndex)
2916 BODY: do
2917 bodyIndex = obs_getBodyIndex(obsSpaceData)
2918 if (bodyIndex < 0) exit BODY
2919 !
2920 ! * Process only ASCAT backscatter anisotropy observations
2921 !
2922 varno = obs_bodyElem_i(obsSpaceData, OBS_VNM , bodyIndex)
2923 if (varno == BUFR_ICES) then
2924
2925 if (present(columnTrlOnTrlLev_opt)) then
2926 conc = col_getElem(columnTrlOnTrlLev_opt,1,headerIndex,'GL')
2927 else
2928 conc = col_getElem(columnTrlOnTrlLev,1,headerIndex,'GL')
2929 end if
2930 trackCellNum = obs_headElem_i(obsSpaceData, OBS_FOV, headerIndex)
2931 idate = obs_headElem_i(obsSpaceData, OBS_DAT, headerIndex)
2932 write(ccyymmdd, FMT='(i8.8)') idate
2933 read(ccyymmdd(5:6), FMT='(i2)') imonth
2934 obsErrStdDev = SQRT(((1.0-conc)*ascatAnisSigmaOpenWater(trackCellNum,imonth))**2 + &
2935 (conc*ascatAnisSigmaIce(trackCellNum,imonth))**2)
2936
2937 call obs_bodySet_r(obsSpaceData, OBS_OER, bodyIndex, obsErrStdDev)
2938
2939 end if
2940 end do BODY
2941 end do HEADER
2942
2943 if (.not.present(columnTrlOnTrlLev_opt)) then
2944 call col_deallocate(columnTrlOnTrlLev)
2945 end if
2946
2947 if (.not. beSilent) write(*,*) myName//': done'
2948
2949 end subroutine oer_setErrBackScatAnisIce
2950
2951 !--------------------------------------------------------------------------
2952 ! chm_read_obs_err_stddev
2953 !--------------------------------------------------------------------------
2954 subroutine chm_read_obs_err_stddev
2955 !
2956 !:Purpose: To read observation errors from auxiliary file or observation
2957 ! file.
2958 !
2959 implicit none
2960
2961 ! Locals:
2962 integer, parameter :: ndim=1
2963 integer :: istnid
2964
2965 ! read the error std. dev. information from the auxiliary file
2966 call chm_read_obs_err_stddev_file
2967
2968 ! set size of observation sub-space std array
2969 allocate(chm_std%obsStdDev(chm_std%n_stnid))
2970
2971 ! read from the observation file if requested
2972 do istnid=1,chm_std%n_stnid
2973 if (chm_std%source(istnid).ge.1) then
2974
2975 ! retrieve data from stats blocks (with bkstp=14 and block_type='DATA')
2976 chm_std%obsStdDev(istnid) = obsf_obsSub_read('CH',chm_std%stnids(istnid),chm_std%element(istnid), &
2977 chm_std%n_lvl(istnid),ndim,bkstp_opt=14,block_opt='DATA', &
2978 match_nlev_opt=chm_std%source(istnid).eq.1)
2979
2980 end if
2981 end do
2982
2983 end subroutine chm_read_obs_err_stddev
2984
2985 !--------------------------------------------------------------------------
2986 ! chm_read_obs_err_stddev_file
2987 !--------------------------------------------------------------------------
2988 subroutine chm_read_obs_err_stddev_file
2989 !
2990 !:Purpose: To read and store observation error std. dev. as needed for CH
2991 ! family obs.
2992 !
2993 implicit none
2994
2995 ! Locals:
2996 integer :: FNOM, FCLOS
2997 integer :: IERR, JLEV, JELM, nulstat, ios, isize, icount
2998 logical :: LnewExists
2999 character(len=11) :: chemAuxObsDataFile = 'obsinfo_chm'
3000 character (len=128) :: ligne
3001 EXTERNAL FNOM,FCLOS
3002
3003 ! Initialization
3004
3005 chm_std%n_stnid=0
3006
3007 ! CHECK THE EXISTENCE OF THE NEW FILE WITH STATISTICS
3008
3009 INQUIRE(FILE=trim(chemAuxObsDataFile),EXIST=LnewExists)
3010 IF (.not.LnewExists) then
3011 WRITE(*,*) '---------------------------------------------------------------'
3012 WRITE(*,*) 'WARNING! chm_read_obs_err_stddev: auxiliary file ' // trim(chemAuxObsDataFile)
3013 WRITE(*,*) 'WARNING! not available. Default CH family stddev to be applied if needed.'
3014 WRITE(*,*) '---------------------------------------------------------------'
3015 return
3016 ENDIF
3017
3018 ! Read observation error std dev. from auxiliary file for constituent data
3019
3020 NULSTAT=0
3021 IERR=FNOM(NULSTAT,trim(chemAuxObsDataFile),'SEQ',0)
3022 IF (IERR .EQ. 0) THEN
3023 open(unit=nulstat, file=trim(chemAuxObsDataFile), status='OLD')
3024 ELSE
3025 CALL utl_abort('chm_read_obs_err_stddev_file: COULD NOT OPEN AUXILIARY FILE ' // trim(chemAuxObsDataFile))
3026 ENDIF
3027
3028 ! Read error standard deviations for constituents if available.
3029 ! (CH family; ozone and others)
3030
3031 ios=0
3032 read(nulstat,'(A)',iostat=ios,err=10,end=10) ligne
3033 do while (trim(adjustl(ligne(1:12))).ne.'SECTION I:')
3034 read(nulstat,'(A)',iostat=ios,err=10,end=10) ligne
3035 end do
3036
3037 ! Read number of observation set sub-families (STNIDs and ...) and allocate space
3038
3039 read(nulstat,*,iostat=ios,err=10,end=10) chm_std%n_stnid
3040 read(nulstat,*,iostat=ios,err=10,end=10) isize
3041
3042 allocate(chm_std%stnids(chm_std%n_stnid))
3043 allocate(chm_std%std_type(chm_std%n_stnid),chm_std%n_lat(chm_std%n_stnid))
3044 allocate(chm_std%source(chm_std%n_stnid),chm_std%ibegin(chm_std%n_stnid))
3045 allocate(chm_std%element(chm_std%n_stnid),chm_std%n_lvl(chm_std%n_stnid))
3046 allocate(chm_std%std1(isize),chm_std%std2(chm_std%n_stnid),chm_std%std3(chm_std%n_stnid))
3047 allocate(chm_std%levels(isize),chm_std%lat(isize))
3048
3049 chm_std%element(:)=0
3050 chm_std%source(:)=0
3051 chm_std%std_type(:)=0
3052 chm_std%n_lvl(:)=1
3053 chm_std%n_lat(:)=1
3054
3055 ! Begin reading for each sub-family
3056 ! Important: Combination of STNID, BUFR element and number of vertical levels
3057 ! to determine association to the observations.
3058
3059 icount=0
3060 STNIDLOOP: do jelm=1,chm_std%n_stnid
3061 chm_std%ibegin(jelm)=icount+1
3062
3063 ! disregard line of dashes
3064 read(nulstat,'(A)',iostat=ios,err=10,end=10) ligne
3065
3066 ! Read STNID (* as wildcard)
3067 read(nulstat,'(2X,A9)',iostat=ios,err=10,end=10) chm_std%stnids(jelm)
3068
3069 ! Read (1) BUFR element,
3070 ! (2) Flag indication if EOR provided from this auxiliary file or
3071 ! to be read from the observation file,
3072 ! (3) Index specifying OER setup method,
3073 ! (4) Number of vertical levels
3074 ! (5) Number of latitudes
3075 !
3076 ! Important: Combination of STNID, BUFR element and number of vertical levels
3077 ! to determine association to the observations.
3078
3079 read(nulstat,*,iostat=ios,err=10,end=10) chm_std%element(jelm),chm_std%source(jelm), &
3080 chm_std%std_type(jelm), chm_std%n_lvl(jelm), chm_std%n_lat(jelm), &
3081 chm_std%std2(jelm), chm_std%std3(jelm)
3082
3083 if (chm_std%n_lvl(jelm).lt.1) chm_std%n_lvl(jelm)=1
3084 if (chm_std%n_lat(jelm).lt.1) chm_std%n_lat(jelm)=1
3085
3086 if (icount+chm_std%n_lvl(jelm)*chm_std%n_lat(jelm).gt.isize) then
3087 write(*,'(10X,"Max array size exceeded: ",I6)') isize
3088 CALL utl_abort('chm_read_obs_err_stddev_file: PROBLEM READING OBSERR STD DEV.')
3089 end if
3090
3091 ! disregard line of dashes
3092 read(nulstat,'(A)',iostat=ios,err=10,end=10) ligne
3093
3094 ! disregard data section if not needed
3095 if (chm_std%std_type(jelm).eq.1.or.chm_std%std_type(jelm).eq.2.or.(chm_std%source(jelm).ge.1.and.chm_std%std_type(jelm).eq.0)) &
3096 cycle STNIDLOOP
3097
3098 if (chm_std%n_lvl(jelm).eq.1.and.chm_std%n_lat(jelm).eq.1) then
3099
3100 ! Read one value only (independent of level and latitude)
3101
3102 icount=icount+1
3103 read(nulstat,*,iostat=ios,err=10,end=10) chm_std%std1(icount)
3104
3105 else if (chm_std%n_lvl(jelm).eq.1.and.chm_std%n_lat(jelm).gt.1) then
3106
3107 ! Value dependent on latitude only
3108
3109 ! Read reference latitudes (must be in order of increasing size)
3110
3111 read(nulstat,*,iostat=ios,err=10,end=10) &
3112 chm_std%lat(icount+1:icount+chm_std%n_lat(jelm))
3113
3114 ! Read OER-related values
3115
3116 read(nulstat,*,iostat=ios,err=10,end=10) &
3117 chm_std%std1(icount+1:icount+chm_std%n_lat(jelm))
3118
3119 icount=icount+chm_std%n_lat(jelm)
3120
3121 else if (chm_std%n_lvl(jelm).gt.1.and.chm_std%n_lat(jelm).eq.1) then
3122
3123 ! Value dependent on vertical level only
3124
3125 do jlev=1,chm_std%n_lvl(jelm)
3126 icount=icount+1
3127
3128 ! Read vertical level and OER-related value.
3129
3130 read(nulstat,*,iostat=ios,err=10,end=10) &
3131 chm_std%levels(icount),chm_std%std1(icount)
3132
3133 end do
3134
3135 else if (chm_std%n_lvl(jelm).gt.1.and.chm_std%n_lat(jelm).gt.1) then
3136
3137 ! Value dependent on vertical level and latitude
3138
3139 ! Read reference latitudes (must be in order of increasing size)
3140 read(nulstat,*,iostat=ios,err=10,end=10) &
3141 chm_std%lat(icount+1:icount+chm_std%n_lat(jelm))
3142 ! write(*, '(10X,20F9.3)') chm_std%lat(icount+1:icount+chm_std%n_lat(jelm))
3143
3144 do jlev=1,chm_std%n_lvl(jelm)
3145
3146 ! Read vertical level and OER-related lat-dependent values.
3147
3148 read(nulstat,*,iostat=ios,err=10,end=10) &
3149 chm_std%levels(icount+jlev), &
3150 chm_std%std1(icount+(jlev-1)*chm_std%n_lat(jelm)+1:icount+jlev*chm_std%n_lat(jelm))
3151
3152 end do
3153 icount=icount+chm_std%n_lat(jelm)*chm_std%n_lvl(jelm)
3154 end if
3155 end do STNIDLOOP
3156
315710 if (ios.gt.0) then
3158 WRITE(*,*) 'File read error message number: ',ios
3159 CALL utl_abort('chm_read_obs_err_stddev_file: PROBLEM READING OBSERR STD DEV.')
3160 end if
3161
316211 CLOSE(UNIT=NULSTAT)
3163 IERR=FCLOS(NULSTAT)
3164
3165 end subroutine chm_read_obs_err_stddev_file
3166
3167 !--------------------------------------------------------------------------
3168 ! chm_obs_err_stddev_index
3169 !--------------------------------------------------------------------------
3170 subroutine chm_obs_err_stddev_index(CSTNID,NLEV,VARNO,ZLAT,ISTNID,JINT)
3171 !
3172 !:Purpose: To return the station ID and latitude indices corresponding to a
3173 ! measurement.
3174 !
3175 implicit none
3176
3177 ! Arguments:
3178 character(len=*), intent(in) :: CSTNID
3179 integer, intent(in) :: NLEV
3180 integer, intent(in) :: VARNO
3181 real(8), intent(in) :: ZLAT
3182 integer, intent(out) :: ISTNID
3183 integer, intent(out) :: JINT
3184
3185 ! Locals:
3186 integer :: JN,ibegin
3187 real(8) :: lat
3188
3189 ! Important: Combination of STNID, BUFR element and number of vertical levels
3190 ! to determine association to the observations.
3191
3192 ! Find stnid with same number of vertical levels and same BUFR element.
3193 ! Note: * in chm_std%stnids stands for a wildcard
3194
3195 ISTNID=0
3196 DO JN=1,chm_std%n_stnid
3197
3198 ! First compare STNID values allowing for * and blanks in
3199 ! chm_std%stnids(JN) as wildcards
3200
3201 IF (utl_stnid_equal(chm_std%stnids(JN),CSTNID)) THEN
3202 IF ((NLEV.EQ.chm_std%n_lvl(JN) .OR. chm_std%source(JN).eq.2) .AND. VARNO.EQ.chm_std%element(JN)) THEN
3203 ISTNID=JN
3204 exit
3205 END IF
3206 END IF
3207
3208 END DO
3209
3210 IF (ISTNID.EQ.0) THEN
3211 write(*,*) 'chm_obs_err_stddev_index: Error std. dev. is unavailable for STNID ' // trim(CSTNID) // &
3212 ' and NLEV = ' // trim(utl_str(NLEV)) // ' and VARNO = ' // trim(utl_str(VARNO))
3213 write(*,*)
3214 write(*,*) ' Contents of chm_std (n_stnid = ' // trim(utl_str(chm_std%n_stnid)) // '):'
3215 if (chm_std%n_stnid.gt.0) then
3216 write(*,'(A)') ' STNID BUFR NLEVELS'
3217 do jn=1,chm_std%n_stnid
3218 write(*,'(2X,A,2X,I12,2X,I10)') chm_std%stnids(jn),chm_std%element(jn),chm_std%n_lvl(jn)
3219 end do
3220 end if
3221 write(*,*)
3222 call utl_abort('chm_obs_err_stddev_index: Check section I of the auxiliary file.')
3223 ELSE
3224
3225 IF (chm_std%n_lat(ISTNID) .GT. 1) THEN
3226
3227 ! Find latitude index for interpolation.
3228 ! Assuming increasing latitudes in chm_std%lat
3229
3230 lat = zlat / MPC_RADIANS_PER_DEGREE_R8 ! radians to degrees
3231
3232 ibegin=chm_std%ibegin(ISTNID)-1
3233 IF (lat .GE. chm_std%lat(ibegin+chm_std%n_lat(ISTNID))) THEN
3234 JINT=chm_std%n_lat(ISTNID)+1
3235 ELSE
3236 DO JINT=1,chm_std%n_lat(ISTNID)
3237 IF (lat .LE. chm_std%lat(ibegin+JINT)) exit
3238 END DO
3239 END IF
3240
3241 END IF
3242 END IF
3243
3244 end subroutine chm_obs_err_stddev_index
3245
3246 !--------------------------------------------------------------------------
3247 ! chm_get_obs_err_stddev
3248 !--------------------------------------------------------------------------
3249 function chm_get_obs_err_stddev(cstnid,nlev,varno,zlat,zlon,idate,itime,zval,&
3250 zlev,ilev,ifirst) result(obs_err_stddev)
3251 !
3252 !:Purpose: To return the observational error std dev for a CH family
3253 ! measurement
3254 implicit none
3255
3256 ! Arguments:
3257 character(len=*), intent(in) :: CSTNID ! station ID
3258 integer, intent(in) :: NLEV ! number of levels
3259 integer, intent(in) :: VARNO ! BUFR number
3260 real(8), intent(in) :: ZLAT ! latitude (radians)
3261 real(8), intent(in) :: ZLON ! longitude (radians)
3262 integer, intent(in) :: IDATE ! date in YYYYMMDD format
3263 integer, intent(in) :: ITIME ! time in HHMM format
3264 real(8), intent(in) :: ZVAL ! observation values
3265 real(8), intent(in) :: ZLEV ! vertical coordinate value
3266 integer, intent(in) :: ILEV ! observation number in the profile
3267 logical, intent(in) :: IFIRST ! true: first call for a profile
3268 ! Result:
3269 real(8) :: obs_err_stddev
3270
3271 ! Locals:
3272 real(8) :: wgt,zwb,sigma
3273 integer :: ibegin,JLEV,JN,stat
3274 integer, save :: ISTNID,JINT
3275
3276 ! If this call is for the first level for this measurement, get
3277 ! the station ID and latitude indices corresponding to this measurement
3278 if (ifirst) call chm_obs_err_stddev_index(CSTNID,NLEV,VARNO,ZLAT,ISTNID,JINT)
3279
3280 ! Get weighting of error std. dev. if required
3281
3282 if (chm_std%std_type(ISTNID).gt.2 .or. &
3283 (chm_std%source(ISTNID).eq.0 .and. chm_std%std_type(ISTNID).eq.0)) then
3284
3285 IF (chm_std%n_lvl(ISTNID) .GT. 1) THEN
3286
3287 ! Find nearest vertical level (no interpolation)
3288
3289 zwb=1.E10
3290 ibegin=chm_std%ibegin(ISTNID)-1
3291 DO JN=1,chm_std%n_lvl(ISTNID)
3292 IF (zwb .GT. abs(ZLEV-chm_std%levels(ibegin+JN))) THEN
3293 JLEV=JN
3294 zwb=abs(ZLEV-chm_std%levels(ibegin+JN))
3295 END IF
3296 END DO
3297 JLEV=ibegin+(JLEV-1)*chm_std%n_lat(ISTNID)+1
3298 ELSE
3299 JLEV=chm_std%ibegin(ISTNID)
3300 END IF
3301
3302 IF (chm_std%n_lat(ISTNID) .GT. 1) THEN
3303
3304 ! Apply interpolation
3305
3306 JLEV=JLEV+JINT-1
3307 ibegin=chm_std%ibegin(ISTNID)-1
3308 IF (JINT.EQ.1.OR.JINT.GT.chm_std%n_lat(ISTNID)) THEN
3309 wgt=chm_std%std1(JLEV)
3310 ELSE
3311 wgt=(chm_std%std1(JLEV-1)*(chm_std%lat(ibegin+JINT)-ZLAT)+ &
3312 chm_std%std1(JLEV)*(ZLAT-chm_std%lat(ibegin+JINT-1)))/ &
3313 (chm_std%lat(ibegin+JINT)-chm_std%lat(ibegin+JINT-1))
3314 END IF
3315 ELSE
3316 wgt=chm_std%std1(JLEV)
3317 END IF
3318
3319 end if
3320
3321 ! Set the error std. dev.
3322
3323 IF (chm_std%source(ISTNID).EQ.0) THEN
3324
3325 ! Set error standard deviations from scratch using content of
3326 ! previously read content of the auxiliary file.
3327
3328 select case(chm_std%std_type(ISTNID))
3329 case(0)
3330 obs_err_stddev = wgt
3331 case(1)
3332 obs_err_stddev = max(chm_std%std3(ISTNID),chm_std%std2(ISTNID)*ZVAL)
3333 case(2)
3334 obs_err_stddev = sqrt(chm_std%std3(ISTNID)**2+(chm_std%std2(ISTNID)*ZVAL)**2)
3335 case(3)
3336 obs_err_stddev = min(chm_std%std3(ISTNID),max(chm_std%std2(ISTNID),wgt*ZVAL))
3337 case(4)
3338 obs_err_stddev = sqrt(chm_std%std2(ISTNID)**2+(wgt*ZVAL)**2)
3339 case default
3340 call utl_abort('chm_get_obs_err_stddev: std_type = ' // trim(utl_str(chm_std%std_type(ISTNID))) // &
3341 ' for STNID = ' // trim(CSTNID) // ' is not recognized.')
3342 end select
3343
3344 ELSE
3345
3346 ! Adjust error standard deviations read from observation file if requested.
3347
3348 sigma = oss_obsdata_get_element(chm_std%obsStdDev(istnid), oss_obsdata_get_header_code(zlon,zlat,idate,itime,cstnid), ilev, stat_opt=stat)
3349
3350 select case(stat)
3351 case(1)
3352 call utl_abort("chm_get_obs_err_stddev: No reports available for STNID = " // trim(cstnid) // &
3353 ", nlev = " // trim(utl_str(nlev)) // ", varno = " // trim(utl_str(varno)))
3354 case(2)
3355 call utl_abort("chm_get_obs_err_stddev: Report not found for STNID = " // trim(cstnid) // &
3356 ", nlev = " // trim(utl_str(nlev)) // ", varno = " // trim(utl_str(varno)))
3357 end select
3358
3359 select case(chm_std%std_type(ISTNID))
3360 case(0)
3361 obs_err_stddev = sigma
3362 case(1)
3363 obs_err_stddev = max(chm_std%std3(ISTNID),chm_std%std2(ISTNID)*sigma)
3364 case(2)
3365 obs_err_stddev = sqrt(chm_std%std3(ISTNID)**2+(chm_std%std2(ISTNID)*sigma)**2)
3366 case(3)
3367 obs_err_stddev = min(chm_std%std3(ISTNID),max(chm_std%std2(ISTNID),wgt*sigma))
3368 case(4)
3369 obs_err_stddev = sqrt(chm_std%std2(ISTNID)**2+(wgt*sigma)**2)
3370 case default
3371 call utl_abort('chm_get_obs_err_stddev: std_type = ' // trim(utl_str(chm_std%std_type(ISTNID))) // &
3372 ' for STNID = ' // trim(CSTNID) // ' is not recognized.')
3373 end select
3374
3375 END IF
3376
3377 end function chm_get_obs_err_stddev
3378
3379 !--------------------------------------------------------------------------
3380 ! chm_dealloc_obs_err_stddev
3381 !--------------------------------------------------------------------------
3382 subroutine chm_dealloc_obs_err_stddev
3383 !
3384 !:Purpose: To deallocate temporary storage space used for observation errors
3385 ! for the CH family.
3386 !
3387 implicit none
3388
3389 ! Locals:
3390 integer :: istnid
3391
3392 if (chm_std%n_stnid.eq.0) return
3393
3394 if (allocated(chm_std%obsStdDev)) then
3395 do istnid=1,chm_std%n_stnid
3396 if (chm_std%source(istnid).ge.1) call oss_obsdata_dealloc(chm_std%obsStdDev(istnid))
3397 end do
3398 deallocate(chm_std%obsStdDev)
3399 end if
3400
3401 if (allocated(chm_std%stnids)) deallocate(chm_std%stnids)
3402 if (allocated(chm_std%n_lvl)) deallocate(chm_std%n_lvl)
3403 if (allocated(chm_std%std_type)) deallocate(chm_std%std_type)
3404 if (allocated(chm_std%ibegin)) deallocate(chm_std%ibegin)
3405 if (allocated(chm_std%element)) deallocate(chm_std%element)
3406 if (allocated(chm_std%source)) deallocate(chm_std%source)
3407 if (allocated(chm_std%n_lat)) deallocate(chm_std%n_lat)
3408 if (allocated(chm_std%std1)) deallocate(chm_std%std1)
3409 if (allocated(chm_std%std2)) deallocate(chm_std%std2)
3410 if (allocated(chm_std%std3)) deallocate(chm_std%std3)
3411 if (allocated(chm_std%levels)) deallocate(chm_std%levels)
3412 if (allocated(chm_std%lat)) deallocate(chm_std%lat)
3413
3414 end subroutine chm_dealloc_obs_err_stddev
3415
3416 !--------------------------------------------------------------------------
3417 ! oer_getSSTdataParam_char
3418 !--------------------------------------------------------------------------
3419 function oer_getSSTdataParam_char(item, itemIndex) result(value)
3420 !
3421 !:Purpose: get character item value from SSTdataParams derived type
3422 !
3423 implicit none
3424
3425
3426 ! Arguments:
3427 character(len=*), intent(in) :: item
3428 integer , intent(in) :: itemIndex
3429 ! Result:
3430 character(len=20) :: value
3431
3432 select case(trim(item))
3433 case('dataType')
3434 value = SSTdataParams(itemIndex)%dataType
3435 case('instrument')
3436 value = SSTdataParams(itemIndex)%instrument
3437 case('sensor')
3438 value = SSTdataParams(itemIndex)%sensor
3439 case('sensorType')
3440 value = SSTdataParams(itemIndex)%sensorType
3441 case default
3442 call utl_abort('oer_getSSTdataParam_char: invalid item '//(trim(item)))
3443 end select
3444
3445 end function oer_getSSTdataParam_char
3446
3447 !--------------------------------------------------------------------------
3448 ! oer_getSSTdataParam_int
3449 !--------------------------------------------------------------------------
3450 function oer_getSSTdataParam_int(item, itemIndex_opt) result(value)
3451 !
3452 !:Purpose: get integer item value from SSTdataParams derived type
3453 !
3454 implicit none
3455
3456 ! Arguments:
3457 character(len=*), intent(in) :: item
3458 integer, optional, intent(in) :: itemIndex_opt
3459 ! Result:
3460 integer :: value
3461
3462 if (present(itemIndex_opt)) then
3463 select case(trim(item))
3464 case('codeType')
3465 value = SSTdataParams(itemIndex_opt)%codeType
3466 case default
3467 call utl_abort('oer_getSSTdataParam_int: invalid item '//(trim(item)))
3468 end select
3469 else
3470 select case(trim(item))
3471 case('maxNumberSSTDatasets')
3472 value = maxNumberSSTDatasets
3473 case('numberSSTDatasets')
3474 value = numberSSTDatasets
3475 case default
3476 call utl_abort('oer_getSSTdataParam_int: invalid item '//(trim(item)))
3477 end select
3478
3479 end if
3480
3481 end function oer_getSSTdataParam_int
3482
3483 !--------------------------------------------------------------------------
3484 ! oer_getSSTdataParam_R8
3485 !--------------------------------------------------------------------------
3486 function oer_getSSTdataParam_R8(item, itemIndex) result(value)
3487 !
3488 !:Purpose: get real(8) item value from SSTdataParams derived type
3489 !
3490 implicit none
3491
3492 ! Arguments:
3493 character(len=*), intent(in) :: item
3494 integer , intent(in) :: itemIndex
3495 ! Result:
3496 real(8) :: value
3497
3498 select case(trim(item))
3499 case('dayError')
3500 value = SSTdataParams(itemIndex)%dayError
3501 case('nightError')
3502 value = SSTdataParams(itemIndex)%nightError
3503 case default
3504 call utl_abort('oer_getSSTdataParam_R8: invalid item '//(trim(item)))
3505 end select
3506
3507 end function oer_getSSTdataParam_R8
3508
3509end module obsErrors_mod
