1module obsFilter_mod
2 ! MODULE obsFilter_mod (prefix='filt' category='5. Observation operators')
3 !
4 !:Purpose: Various types of filters that are applied to the observations
5 ! mostly to reject them so that they will not be assimilated.
6 !
7 use codePrecision_mod
8 use midasMpi_mod
9 use earthConstants_mod
10 use obsSpaceData_mod
11 use columnData_mod
12 use bufr_mod
13 use tovsNL_mod
14 use gps_mod
15 use utilities_mod
16 use varNameList_mod
17 use physicsFunctions_mod
18 use codtyp_mod
19 use radialVelocity_mod
20 implicit none
21 save
22 private
23
24 ! public variables
25 public :: filt_rlimlvhu
26 ! public procedures
27 public :: filt_setup, filt_topo, filt_suprep
28 public :: filt_surfaceWind, filt_gpsro, filt_backScatAnisIce, filt_iceConcentration, filt_radvel
29 public :: filt_bufrCodeAssimilated, filt_getBufrCodeAssimilated, filt_nBufrCodeAssimilated
30 integer, parameter :: nelemsMax = 30
31 integer, parameter :: nflagsMax = 15
32 integer :: filt_nelems, filt_nflags
33 integer, target :: filt_nlist(nelemsMax)
34 integer :: filt_nlistflg(nflagsMax)
35
36 real(8) :: filt_rlimlvhu
37
38 ! topographic rejection criteria
39 integer, parameter :: numElem = 21
40 real(8) :: altDiffMax(numElem) = & ! default values (in metres)
41 (/ 50.d0, 50.d0, 50.d0, 50.d0, 50.d0, 800.d0, 800.d0, &
42 800.d0, 800.d0, 1000.d0, 50.d0, 50.d0, 50.d0, 50.d0, &
43 50.d0, 50.d0, 50.d0, 50.d0, 50.d0, 50.d0, 1000.d0 /)
44 integer, parameter :: elemList(numElem) = &
45 (/ BUFR_NEDS, BUFR_NEFS, BUFR_NEUS, BUFR_NEVS, BUFR_NESS, BUFR_NETS, BUFR_NEPS, &
46 BUFR_NEPN, BUFR_NEGZ, BUFR_NEZD, BUFR_NEDD, BUFR_NEFF, BUFR_NEUU, BUFR_NEVV, &
47 BUFR_NEES, BUFR_NETT, BUFR_NEAL, bufr_vis , bufr_logVis, bufr_gust, bufr_radvel /)
48
49 integer, parameter :: nTopoFiltFam = 8
50 character(len=2) :: filtTopoList(nTopoFiltFam) = ' '
51
52 character(len=48) :: filterMode
53
54 logical :: initialized = .false.
55
56 ! Namelist variables:
57 logical :: discardlandsfcwind ! choose to reject surface wind obs over land
58 real(8) :: surfaceBufferZone_Pres ! height of buffer zone (in Pa) for rejecting obs near sfc
59 real(8) :: surfaceBufferZone_Height ! height of buffer zone (in m) for rejecting obs near sfc
60 logical :: useEnkfTopoFilt ! choose to use simpler approach (originally in EnKF) for rejecting obs near sfc
61 logical :: rejectGZforAnalysis ! whether to reject geopotential height for analysis update
62
63contains
64
65 !--------------------------------------------------------------------------
66 ! findElemIndex
67 !-------------------------------------------------------------------------
68 function findElemIndex(varNum) result(listIndex)
69 implicit none
70
71 ! Arguments:
72 integer, intent(in) :: varNum
73 ! Result:
74 integer :: listIndex
75
76 ! Locals:
77 integer :: elemIndex
78
79 listIndex = -1
80 do elemIndex=1,numElem
81 if(varNum == elemList(elemIndex)) listIndex = elemIndex
82 end do
83
84 if (listIndex == -1) then
85 write(*,*) 'filterobs_mod-findElemIndex: WARNING: varNum value not found: ',varNum
86 end if
87
88 end function findElemIndex
89
90 !--------------------------------------------------------------------------
91 ! filt_setup
92 !--------------------------------------------------------------------------
93 subroutine filt_setup(filterMode_in)
94 implicit none
95
96 ! Arguments:
97 character(len=*), intent(in) :: filterMode_in
98
99 ! Locals:
100 integer :: nulnam, ierr, elem, jflag, ibit, obsFamilyIndex, elemIndex
101 integer :: fnom, fclos
102 integer :: flagIndex, elementIndex
103
104 character(len=35) :: CREASON(-12:13)
105 data creason/'PASSIVE ', &
106 'SIMULATED ', &
107 'CLOUD EFFECTED RADIANCE ', &
108 'REJECT BY RAW TO BURP ', &
109 'JACOBIAN IMPORTANT ABOVE MODEL TOP', &
110 'ABS OROGRAPH-PHI ', &
111 'MASQUE TERRE-MER ', &
112 'OROGRAPHIE ', &
113 'REJECTED BY QCVAR ', &
114 'REJECTED BY BACKGROUND CHECK ', &
115 'BACKGROUND CHECK LEVEL 3 ', &
116 'BACKGROUND CHECK LEVEL 2 ', &
117 'BACKGROUND GHECK LEVEL 1 ', &
118 'RESERVED ', &
119 'REJECTED BY SELECTION PROCESS ', &
120 'GENERATED BY OI ', &
121 'REJECTION BY OI ', &
122 'ELEMENT ON BLACK LIST ', &
123 'RESERVED ', &
124 'CORRECTED ELEMENT ', &
125 'INTERPOLATED ELEMENT ', &
126 'DOUBTFUL ELEMENT ', &
127 'POSSIBLY ERRONEOUS ELEMENT ', &
128 'ERRONEOUS ELEMENT ', &
129 'ELEMENT EXCEEDS CLIMATE EXTREME ', &
130 'ELEMENT MODIFIED OR GEN BY ADE '/
131
132 ! Namelist variables: (local)
133 integer :: nelems ! MUST NOT BE INCLUDED IN NAMELIST!
134 integer :: nlist(nelemsMax) ! list of bufr element IDs to consider for assimilation
135 integer :: nflags ! MUST NOT BE INCLUDED IN NAMELIST!
136 integer :: nlistflg(nFLagsMax) ! list of flag 'reference numbers' to use for rejecting obs
137 integer :: nelems_altDiffMax ! MUST NOT BE INCLUDED IN NAMELIST!
138 integer :: list_altDiffMax(numElem) ! list of bufr element IDs to apply maximum altitude difference
139 character(len=2) :: list_topoFilt(nTopoFiltFam) ! list of obs family names for applying max altitude
140 real(8) :: value_altDiffMax(numElem) ! value of maximum difference between model sfc and obs altitude
141 real(8) :: rlimlvhu ! pressure level (in hPa) above which humidity (ES) obs are rejected
142
143 namelist /namfilt/nelems, nlist, nflags, nlistflg, rlimlvhu, discardlandsfcwind, &
144 nelems_altDiffMax, list_altDiffMax, value_altDiffMax, surfaceBufferZone_Pres, &
145 surfaceBufferZone_Height, list_topoFilt, useEnkfTopoFilt, rejectGZforAnalysis
146
147 filterMode = filterMode_in
148
149 ! set default values for namelist variables
150 nlist(:) = MPC_missingValue_INT
151 nelems = MPC_missingValue_INT
152 nlistflg(:) = MPC_missingValue_INT
153 nflags = MPC_missingValue_INT
154 list_altDiffMax (:) = MPC_missingValue_INT
155 value_altDiffMax(:) = MPC_missingValue_R8
156 nelems_altDiffMax = MPC_missingValue_INT
157
158 list_topoFilt(:) = '**'
159
160 rlimlvhu = 300.d0
161 discardlandsfcwind = .true.
162
163 surfaceBufferZone_Pres = 5000.0d0 ! default value in Pascals
164 surfaceBufferZone_Height = 400.0d0 ! default value in Metres
165
166 useEnkfTopoFilt = .false.
167 rejectGZforAnalysis = .true.
168
169 nulnam=0
170 ierr=fnom(nulnam,'./flnml','FTN+SEQ+R/O',0)
171 read(nulnam,nml=namfilt,iostat=ierr)
172 if(ierr.ne.0) call utl_abort('filt_setup: Error reading namelist! Hint: did you replace ltopofilt by list_topoFilt?')
173 if(mmpi_myid == 0) write(*,nml=namfilt)
174 ierr=fclos(nulnam)
175
176 filt_rlimlvhu = rlimlvhu
177
178 if (nelems /= MPC_missingValue_INT) then
179 call utl_abort('filt_setup: check NAMFILT namelist section; NELEMS should be removed')
180 end if
181 filt_nelems = 0
182 do elementIndex = 1, nelemsMax
183 if (nlist(elementIndex) /= MPC_missingValue_INT) then
184 filt_nlist(elementIndex) = nlist(elementIndex)
185 filt_nelems = filt_nelems + 1
186 end if
187 end do
188
189 if (nflags /= MPC_missingValue_INT) then
190 call utl_abort('filt_setup: check NAMFILT namelist section; NFLAGS should be removed')
191 end if
192 filt_nflags = 0
193 do flagIndex = 1, nflagsMax
194 if (nlistflg(flagIndex) /= MPC_missingValue_INT) then
195 filt_nlistflg(flagIndex) = nlistflg(flagIndex)
196 filt_nflags = filt_nflags + 1
197 end if
198 end do
199
200 if(mmpi_myid == 0) then
201 write(*,'(1X,"***********************************")')
202 write(*,'(1X," ELEMENTS SELECTED FOR ASSIMILATION:",/)')
203 write(*,'(1X,"***********************************")')
204 do elem=1,filt_nelems
205 write(*,'(15X,I5)') filt_nlist(elem)
206 end do
207 write(*,'(1X,"***********************************")')
208 write(*,*) ' REJECT ELEMENTS WITH REJECT FLAG '
209 write(*,*)' BIT : '
210 do jflag=1,filt_nflags
211 ibit= filt_nlistflg(jflag)
212 write(*,*) ibit,' ',creason(ibit)
213 end do
214 write(*,'(1X,"***********************************")')
215 end if
216
217 !
218 !- Set values for altDiffMax
219 !
220 if (nelems_altDiffMax /= MPC_missingValue_INT) then
221 call utl_abort('filt_setup: check namelist section NAMFILT: nelems_altDiffMax should be removed')
222 end if
223 do elem = 1, numElem
224 if (list_altDiffMax(elem) /= MPC_missingValue_INT .and. value_altDiffMax(elem) /= MPC_missingValue_R8) then
225 elemIndex = findElemIndex(list_altDiffMax(elem))
226 if ( elemIndex >= 1 .and. elemIndex <= numElem ) then
227 altDiffMax(elemIndex) = value_altDiffMax(elem)
228 write(*,*) ' filt_setup: altDiffMax value for ', elemList(elemIndex), ' is set to ', altDiffMax(elemIndex)
229 else
230 call utl_abort('filt_setup: Error in value setting for altDiffMax')
231 end if
232 end if
233 end do
234
235 !
236 !- Set the topographic rejection list
237 !
238 if (all(list_topoFilt(:) == '**')) then
239 ! default list
240 filtTopoList(1) = 'SF'
241 filtTopoList(2) = 'UA'
242 filtTopoList(3) = 'AI'
243 filtTopoList(4) = 'SW'
244 filtTopoList(5) = 'PR'
245 filtTopoList(6) = 'AL'
246 filtTopoList(7) = 'TO'
247 filtTopoList(8) = 'CH'
248 else
249 do obsFamilyIndex = 1, nTopoFiltFam
250 if (list_topoFilt(obsFamilyIndex) /= '**') then
251 filtTopoList(obsFamilyIndex) = list_topoFilt(obsFamilyIndex)
252 end if
253 end do
254 end if
255
256 initialized = .true.
257
258 end subroutine filt_setup
259
260 !--------------------------------------------------------------------------
261 ! filt_suprep
262 !--------------------------------------------------------------------------
263 subroutine filt_suprep(obsSpaceData)
264 !
265 !:Purpose: Select the data in the obsSpaceData which are to be assimilated
266 !
267 implicit none
268
269 ! Arguments:
270 type(struct_obs), intent(inout) :: obsSpaceData
271
272 ! Locals:
273 integer :: bodyIndex, headerIndex
274 integer :: ipres, ivco, ierr, loopIndex
275 integer :: idburp, ivnm, iflg, ibad, iknt, iknt_mpiglobal, ilansea
276 logical :: llok, llrej, llbogus
277
278 call utl_tmg_start(22,'----ObsFiltSuprep')
279
280 if(mmpi_myid == 0) write(*,*) 'starting subroutine filt_suprep'
281
282 iknt = 0
283
284 BODY: do bodyIndex = 1, obs_numbody( obsSpaceData )
285 headerIndex = obs_bodyElem_i( obsSpaceData, OBS_HIND, bodyIndex )
286 ivnm = obs_bodyElem_i( obsSpaceData, OBS_VNM , bodyIndex )
287 idburp = obs_headElem_i( obsSpaceData, OBS_ITY , headerIndex )
288 !
289 ! Unwanted data types via types specified in NLIST
290 !
291 llok = .false.
292 do loopIndex = 1, filt_nelems
293 llok = ( ivnm == filt_nlist( loopIndex ) ) .or. llok
294 end do
295 if (.not.llok) then
296 call obs_bodySet_i(obsSpaceData,OBS_ASS,bodyIndex,obs_notAssimilated)
297 cycle BODY
298 end if
299 !
300 ! Allow gz for bogus data only in analysis case
301 !
302 llbogus = ( idburp == 150 .or. idburp == 151 .or. idburp == 152 .or. idburp == 153 )
303 if ( (filterMode == 'analysis' .or. filterMode == 'FSO') .and. llok .and. ivnm == BUFR_NEGZ .and. .not.llbogus ) then
304 if (rejectGZforAnalysis) llok=.false.
305 end if
306 !
307 ! Ground-based GPS (GP) data (codtyp 189)
308 ! LLOK = .TRUE. DY DEFAULT IF ELEMENT IS IN NLIST
309 ! If gps_gb_LASSMET = .FALSE. don't want to assimilate Ps (BUFR_NEPS),
310 ! Ts (BUFR_NETS), or (T-Td)s (BUFR_NESS)
311 !
312 if ( idburp == 189 ) then
313 if (.not.gps_gb_lassmet .and. ( ivnm == BUFR_NEPS .or. &
314 ivnm == BUFR_NETS .or. &
315 ivnm == BUFR_NESS )) then
316 llok = .false.
317 end if
318 end if
319 !
320 ! Exclude T-Td above level RLIMLVHU (mbs)
321 !
322 ivco = obs_bodyElem_i( obsSpaceData, OBS_VCO, bodyIndex )
323 ipres = nint( obs_bodyElem_r( obsSpaceData, OBS_PPP, bodyIndex))
324 if ( ( ivco == 2) .and. ( ivnm == BUFR_NEES ) .and. &
325 ( ipres < nint( filt_rlimlvhu *100.0d0 )) ) then
326 llok=.false.
327 end if
328 !
329 ! Bad data with quality control flags via bit list specified in NLISTFLG
330 !
331 iflg = obs_bodyElem_i( obsSpaceData, OBS_FLG, bodyIndex )
332 llrej = .false.
333 do loopIndex = 1, filt_nflags
334 ibad = 13 - filt_nlistflg( loopIndex )
335 llrej=( btest(iflg,ibad) ) .or. llrej
336 end do
337 !
338 ! Filter TOVS data: check for invalid land/sea/sea-ice flag
339 !
340 if (ivnm == BUFR_NBT1 .or. ivnm == BUFR_NBT2 .or. ivnm == BUFR_NBT3) then
341 ! Here we have to exclude the ssmis and mwhs2 data since the burp file does not contain the ele 8021 land_sea
342 if ( ( tvs_isIdBurpTovs(idburp) ) .and. &
343 ( idburp /= codtyp_get_codtyp('ssmis') ) .and. &
344 ( idburp /= codtyp_get_codtyp('mwhs2') ) ) then
345 ilansea = tvs_ChangedStypValue( obsSpaceData, headerIndex )
346 if (ilansea < 0 .or. ilansea > 2 ) llok = .false.
347 end if
348 end if
349 !
350 ! SAR winds: assimilates wind speed only for SAR winds
351 !
352 if ( ivnm == BUFR_NEFS ) then
353 if ( idburp .ne. 204 ) then
354 llok = .false.
355 end if
356 end if
357 if ( llok .and. .not. llrej ) then
358 call obs_bodySet_i( obsSpaceData, OBS_ASS, bodyIndex, obs_assimilated )
359 iknt = iknt + 1
360 else
361 call obs_bodySet_i( obsSpaceData, OBS_ASS, bodyIndex, obs_notAssimilated )
362 end if
363
364 end do body
365
366 call rpn_comm_allreduce( iknt, iknt_mpiglobal, 1, "MPI_INTEGER", "MPI_SUM", "GRID", ierr )
367 if(mmpi_myid == 0) write(*,*) ' Number of data to be assimilated: ', iknt_mpiglobal
368
369 if(mmpi_myid == 0) write(*,*) 'end of filt_suprep'
370
371 ! abort if there is no data to be assimilated
372 if (iknt_mpiglobal == 0 ) then
373 call utl_abort('SUPREP. NO DATA TO BE ASSIMILATED')
374 end if
375
376 call utl_tmg_stop(22)
377
378 end subroutine filt_suprep
379
380 !--------------------------------------------------------------------------
381 ! filt_topo
382 !--------------------------------------------------------------------------
383 subroutine filt_topo(columnTrlOnTrlLev, obsSpaceData, beSilent)
384 implicit none
385
386 ! Arguments:
387 type(struct_columnData), intent(in) :: columnTrlOnTrlLev
388 type(struct_obs) , intent(inout) :: obsSpaceData
389 logical, intent(in) :: beSilent
390
391 if (all(filtTopoList(:) == ' ')) then
392
393 if (mmpi_myid == 0 .and. .not.beSilent) then
394 write(*,*)
395 write(*,*)' --------------------------------------------------------------'
396 write(*,*)' - filt_topo: NO topographic filtering '
397 write(*,*)' --------------------------------------------------------------'
398 end if
399
400 else
401
402 if (mmpi_myid == 0 .and. .not.beSilent) then
403 write(*,*)
404 write(*,*)' --------------------------------------------------------------'
405 write(*,*)' - filt_topo: topographic filtering of the following obs family'
406 write(*,*)' --------------------------------------------------------------'
407 end if
408
409 if (any(filtTopoList(:) == 'SF')) call filt_topoSurface (columnTrlOnTrlLev,obsSpaceData,beSilent)
410 if (any(filtTopoList(:) == 'UA')) call filt_topoRadiosonde(columnTrlOnTrlLev,obsSpaceData,beSilent)
411 if (any(filtTopoList(:) == 'AI')) call filt_topoAISW (columnTrlOnTrlLev,obsSpaceData,'AI',beSilent)
412 if (any(filtTopoList(:) == 'SW')) call filt_topoAISW (columnTrlOnTrlLev,obsSpaceData,'SW',beSilent)
413 if (any(filtTopoList(:) == 'PR')) call filt_topoProfiler (columnTrlOnTrlLev,obsSpaceData,beSilent)
414 if (any(filtTopoList(:) == 'AL')) call filt_topoAladin (columnTrlOnTrlLev,obsSpaceData,beSilent)
415 if (any(filtTopoList(:) == 'TO')) call filt_topoTovs (columnTrlOnTrlLev,obsSpaceData,beSilent)
416 if (any(filtTopoList(:) == 'CH')) call filt_topoChemistry (columnTrlOnTrlLev,obsSpaceData,beSilent)
417
418 end if
419
420 end subroutine filt_topo
421
422 !--------------------------------------------------------------------------
423 ! filt_topoSurface
424 !--------------------------------------------------------------------------
425 subroutine filt_topoSurface(columnTrlOnTrlLev, obsSpaceData, beSilent)
426 !
427 !:Purpose: Refuse elements which are too far away from the surface.
428 ! Replace the pressure of elements which are slightly below
429 ! the model surface by the pressure of the trial field.
430 !
431 implicit none
432
433 ! Arguments:s
434 type(struct_columnData), intent(in) :: columnTrlOnTrlLev
435 type(struct_obs) , intent(inout) :: obsSpaceData
436 logical , intent(in) :: beSilent
437
438 ! Locals:
439 real(8) :: altitudeDiff
440 integer :: headerIndex, bodyIndex, familyIndex, elemIndex
441 integer :: ivnm,countAssim
442 integer :: countAcc(numElem),countRej(numElem)
443 integer, parameter :: numFamily = 3
444 character(len=2) :: list_family(numFamily)
445 character(len=4) :: varLevel
446
447 !
448 ! reset gps_gb_dzmax for gb gps ztd to value in namelist file
449 !
450 altDiffMax(findElemIndex(BUFR_NEZD)) = gps_gb_dzmax
451
452 if ( .not.beSilent ) then
453 write(*,*) ' '
454 write(*,*) ' filt_topoSurface: '
455 write(*,*) ' '
456 write(*,*) '*****************************************************'
457 write(*,222) 'ELEMENTS ',(elemList(elemIndex),elemIndex=1,numElem)
458 write(*,223) 'REJECTION BOUNDARY(METRE) ',(altDiffMax(elemIndex),elemIndex=1,numElem)
459 write(*,*) '*****************************************************'
460 write(*,*) ' '
461 end if
462
463 ! Loop over the families of interest
464 list_family(1) = 'SF'
465 list_family(2) = 'UA'
466 list_family(3) = 'GP'
467 FAMILY: do familyIndex = 1, numFamily
468
469 ! set counters to zero
470 countRej(:)=0
471 countAcc(:)=0
472
473 ! loop over all header indices of each family
474 call obs_set_current_header_list(obsSpaceData, &
475 list_family(familyIndex))
476 HEADER: do
477 headerIndex = obs_getHeaderIndex(obsSpaceData)
478 if (headerIndex < 0) exit HEADER
479
480 ! loop over all body indices (still in the same family)
481 call obs_set_current_body_list(obsSpaceData, headerIndex)
482 BODY: do
483 bodyIndex = obs_getBodyIndex(obsSpaceData)
484 if (bodyIndex < 0) exit BODY
485
486 ! skip this obs if it is not on height levels
487 if (obs_bodyElem_i(obsSpaceData,OBS_VCO,bodyIndex).ne.1) cycle BODY
488
489 ! skip this obs if already flagged to not be assimilated
490 if (obs_bodyElem_i(obsSpaceData,OBS_ASS,bodyIndex) == obs_notAssimilated) cycle BODY
491
492 ivnm = obs_bodyElem_i(obsSpaceData,OBS_VNM,bodyIndex)
493 varLevel = vnl_varLevelFromVarnum(ivnm)
494 altitudeDiff = abs( obs_bodyElem_r(obsSpaceData,OBS_PPP,bodyIndex) - &
495 col_getHeight(columnTrlOnTrlLev,col_getNumLev(columnTrlOnTrlLev,varLevel),headerIndex,varLevel) )
496 !
497 ! apply filter to selected elements
498 !
499 elemIndex = findElemIndex(ivnm)
500 if (elemIndex == -1) cycle BODY
501
502 if (altitudeDiff <= altDiffMax(elemIndex)) then
503 ! obs passes the acceptance criteria
504 countAcc(elemIndex) = countAcc(elemIndex)+1
505 else
506 ! obs does not pass the acceptance criteria, reject it
507 call obs_bodySet_i( obsSpaceData,OBS_FLG,bodyIndex, &
508 ibset(obs_bodyElem_i(obsSpaceData,OBS_FLG,bodyIndex),18) )
509 call obs_bodySet_i(obsSpaceData,OBS_ASS,bodyIndex,obs_notAssimilated)
510 countRej(elemIndex) = countRej(elemIndex)+1
511 end if
512 end do BODY
513 end do HEADER
514
515 if ( .not.beSilent ) then
516 write(*,*) ' '
517 write(*,*) '*****************************************************'
518 write(*,*) 'FAMILY = ',list_family(familyIndex)
519 write(*,222) 'ELEMENTS ', (elemList(elemIndex),elemIndex=1,numElem)
520 write(*,222) 'ACCEPTED ',(countAcc(elemIndex),elemIndex=1,numElem)
521 write(*,222) 'REJECTED ',(countRej(elemIndex),elemIndex=1,numElem)
522 write(*,*) '*****************************************************'
523 write(*,*) ' '
524 end if
525222 format(2x,a29,16(2x,i5))
526223 format(2x,a29,16(2x,f5.0))
527
528 end do FAMILY
529
530 countAssim=0
531 do bodyIndex=1,obs_numbody(obsSpaceData)
532 if ( obs_bodyElem_i(obsSpaceData,OBS_ASS,bodyIndex) == obs_assimilated) countAssim=countAssim+1
533 end do
534 if ( .not.beSilent ) write(*,'(1X," NUMBER OF DATA TO BE ASSIMILATED AFTER ADJUSTMENTS: ",i10)') countAssim
535 if ( .not.beSilent ) write(*,* ) ' '
536
537 end subroutine filt_topoSurface
538
539 !--------------------------------------------------------------------------
540 ! filt_topoRadiosonde
541 !--------------------------------------------------------------------------
542 subroutine filt_topoRadiosonde(columnTrlOnTrlLev, obsSpaceData, beSilent)
543 !
544 !:Purpose: Refuse elements which are too far away from the surface of the model
545 ! Refuse elements which are considered in the free atmosphere of
546 ! the RAOB but fall in the surface boundary layer of the model atmosphere.
547 !
548 implicit none
549
550 ! Arguments:
551 type(struct_columnData), intent(in) :: columnTrlOnTrlLev
552 type(struct_obs), intent(inout) :: obsSpaceData
553 logical, intent(in) :: beSilent
554
555 ! Locals:
556 integer :: headerIndex, bodyIndex, listIndex, elemIndex
557 integer :: ivnm, countAssim
558 integer :: itotacc(numElem), itotrej(numElem), isblrej(numElem)
559 integer :: igzacc(numElem), igzrej(numElem), ibndrej(numElem)
560 real(8) :: zval, obsPressure, altitudeDiff
561 real(8) :: obsSfcAltitude, colSfcAltitude, colPressureBelow, colPressureAbove, zdelp
562 logical :: llok
563 real(8) :: geopotential(1), height(1)
564 integer :: nlev_M
565 real(8) :: lat
566
567 if ( .not.beSilent ) then
568 write(*,*) ' '
569 write(*,*) ' filt_topoRadiosonde: '
570 write(*,*) ' '
571 write(*,*) '************************************************'
572 write(*,222) ' ELEMENTS ',(elemList(elemIndex),elemIndex=1,numElem)
573 write(*,223) ' REJECTION BOUNDARY(METRE) ',(altDiffMax(elemIndex),elemIndex=1,numElem)
574 write(*,223) ' REJECTION SBL (PASCAL) ',(surfaceBufferZone_Pres,elemIndex=1,numElem)
575 write(*,*) '************************************************'
576 write(*,*) ' '
577 end if
578
579 ! set counters to zero
580 itotrej(:)=0
581 itotacc(:)=0
582 isblrej(:)=0
583 igzacc(:)=0
584 igzrej(:)=0
585 ibndrej(:)=0
586
587 nlev_M = col_getNumLev(columnTrlOnTrlLev,'MM')
588
589 ! loop over all header indices of the 'UA' family
590 call obs_set_current_header_list(obsSpaceData, 'UA')
591 HEADER: do
592 headerIndex = obs_getHeaderIndex(obsSpaceData)
593 if (headerIndex < 0) exit HEADER
594
595 ! HEIGHT GZ
596
597 ! loop over all body indices (still in the 'UA' family)
598 call obs_set_current_body_list(obsSpaceData, headerIndex)
599 BODY: do
600 bodyIndex = obs_getBodyIndex(obsSpaceData)
601 if (bodyIndex < 0) exit BODY
602
603 ! skip this obs if it is not on pressure level
604 if( obs_bodyElem_i(obsSpaceData,OBS_VCO,bodyIndex).ne.2 ) cycle BODY
605
606 ! skip this obs if already flagged to not be assimilated
607 if( obs_bodyElem_i(obsSpaceData,OBS_ASS,bodyIndex) == obs_notAssimilated ) cycle BODY
608
609 ! skip this obs if it is not GZ
610 ivnm=obs_bodyElem_i(obsSpaceData,OBS_VNM,bodyIndex)
611 listIndex = findElemIndex(ivnm)
612 llok = (ivnm == BUFR_NEGZ .and. listIndex.ne.-1)
613 if (.not. llok ) cycle BODY
614
615 ! convert altitude read from column to geopotential
616 height(1) = col_getHeight(columnTrlOnTrlLev,0,headerIndex,'SF')
617 lat = obs_headElem_r(obsSpaceData,OBS_LAT,headerIndex)
618 call phf_height2geopotential(height,lat,geopotential)
619
620 zval = obs_bodyElem_r(obsSpaceData,OBS_VAR,bodyIndex)
621 altitudeDiff = ( zval - geopotential(1) )/ec_rg
622 ! obs is above surface, so it is ok, lets jump to the next obs
623 if(altitudeDiff >= 0.0d0) cycle BODY
624
625 if(altitudeDiff >= -1.0d0*altDiffMax(listIndex)) then
626 ! obs is an acceptably small distance below the surface
627 itotacc(listIndex) = itotacc(listIndex)+1
628 igzacc(listIndex) = igzacc(listIndex)+1
629 else
630 ! too far below surface, reject
631 call obs_bodySet_i(obsSpaceData,OBS_FLG,bodyIndex, &
632 ibset( obs_bodyElem_i(obsSpaceData,OBS_FLG,bodyIndex), 18 ))
633 call obs_bodySet_i(obsSpaceData,OBS_ASS,bodyIndex,obs_notAssimilated)
634 itotrej(listIndex) = itotrej(listIndex)+1
635 igzrej(listIndex) = igzrej(listIndex)+1
636 end if
637 end do BODY
638 !
639 ! REJECT ELEMENTS OF U,V,T-TD,T BELOW THE MODEL SURFACE
640 ! AND THOSE NON SURFACE ELEMENTS PRESENT IN THE SURFACE
641 ! BOUNDARY LAYER OF THE RAOB OR OF THE MODEL.
642 ! AT THIS POINT WE WANT TO KEEP OBSERVATIONS IN THE FREE
643 ! ATMOSPHERE
644 !
645 !---Special case if station elevation is above model elevation
646 ! we want to define colPressureAbove at a level above the station.
647 ! To approximate that value, we will transform the difference
648 ! between the 2 elevations into a difference in pressure using
649 ! the rule of thumb (1Mb =8 metres)
650 !---Even though TT(element=12001) is not assimmilated
651 ! it is treated as if it were for the evaluation step.
652 ! Otherwise we use observations of TT that are too far
653 ! from the model topography in the verification.
654
655 obsSfcAltitude = obs_headElem_r(obsSpaceData,OBS_ALT,headerIndex)
656 colSfcAltitude = col_getHeight(columnTrlOnTrlLev,0,headerIndex,'SF')
657 altitudeDiff = obsSfcAltitude - colSfcAltitude
658
659 ! Set the body list & start at the beginning of the list
660 call obs_set_current_body_list(obsSpaceData, headerIndex)
661 BODY2: do
662 bodyIndex = obs_getBodyIndex(obsSpaceData)
663 if (bodyIndex < 0) exit BODY2
664
665 if( obs_bodyElem_i(obsSpaceData,OBS_VCO,bodyIndex).ne.2 ) cycle BODY2
666
667 ! skip this obs if already flagged to not be assimilated
668 if( obs_bodyElem_i(obsSpaceData,OBS_ASS,bodyIndex) == obs_notAssimilated ) cycle BODY2
669
670 ivnm = obs_bodyElem_i(obsSpaceData,OBS_VNM,bodyIndex)
671 listIndex = findElemIndex(ivnm)
672 llok = (ivnm.ne.BUFR_NEGZ .and. listIndex.ne.-1)
673 if (.not. llok ) cycle BODY2 ! Proceed with the next bodyIndex
674
675 obsPressure = obs_bodyElem_r(obsSpaceData,OBS_PPP,bodyIndex)
676 colPressureBelow = col_getElem(columnTrlOnTrlLev,1,headerIndex,'P0')
677 colPressureAbove = colPressureBelow - surfaceBufferZone_Pres
678
679 if (useEnkfTopoFilt) then
680 ! Simpler rules used in the EnKF
681 if(obsPressure >= colPressureAbove ) then
682 if(abs(altitudeDiff) >= 50.0D0 .or. obsPressure >= colPressureBelow) then
683 call obs_bodySet_i(obsSpaceData,OBS_FLG,bodyIndex, &
684 ibset(obs_bodyElem_i(obsSpaceData,OBS_FLG,bodyIndex), 18 ))
685 call obs_bodySet_i(obsSpaceData,OBS_ASS,bodyIndex,obs_notAssimilated)
686 itotrej(listIndex) = itotrej(listIndex) + 1
687 ibndrej(listIndex) = ibndrej(listIndex) + 1
688 end if
689 end if
690 else
691 ! Original (and confusing) rules used in Var
692 if (altitudeDiff > 0.0d0) then
693 zdelp = altitudeDiff * 100.d0 / 8.0d0
694 colPressureAbove = colPressureBelow - (zdelp + surfaceBufferZone_Pres)
695 end if
696
697 if(abs(altitudeDiff) <= altDiffMax(listIndex)) then
698 !--Model surface and station altitude are very close
699 ! Accept observation if obsPressure is within the domain
700 ! of the trial field
701 colPressureAbove = col_getPressure(columnTrlOnTrlLev,col_getNumLev(columnTrlOnTrlLev,'MM')-1,headerIndex,'MM')
702 end if
703 if(obsPressure > colPressureBelow ) then
704 call obs_bodySet_i(obsSpaceData,OBS_FLG,bodyIndex, &
705 ibset(obs_bodyElem_i(obsSpaceData,OBS_FLG,bodyIndex), 18 ))
706 call obs_bodySet_i(obsSpaceData,OBS_ASS,bodyIndex,obs_notAssimilated)
707 itotrej(listIndex) = itotrej(listIndex) + 1
708 ibndrej(listIndex) = ibndrej(listIndex) + 1
709 else if(obsPressure <= colPressureBelow .and. obsPressure > colPressureAbove ) then
710 call obs_bodySet_i(obsSpaceData,OBS_FLG,bodyIndex, &
711 ibset( obs_bodyElem_i(obsSpaceData,OBS_FLG,bodyIndex), 18 ))
712 call obs_bodySet_i(obsSpaceData,OBS_ASS,bodyIndex,obs_notAssimilated)
713 itotrej(listIndex) = itotrej(listIndex) + 1
714 isblrej(listIndex) = isblrej(listIndex) + 1
715 end if
716 end if
717
718 end do BODY2
719 end do HEADER
720
721 if ( .not.beSilent ) then
722 write(*,*) ' '
723 write(*,*) '***************************************'
724 write(*,*) 'FAMILY = UA'
725 write(*,222) ' ELEMENTS ', (elemList(elemIndex),elemIndex=1,numElem)
726 write(*,222) ' ACC GZ EXT ',(igzacc(elemIndex),elemIndex=1,numElem)
727 write(*,222) ' ACC TOTAL ',(itotacc(elemIndex),elemIndex=1,numElem)
728 write(*,*) '***************'
729 write(*,222) ' REJ GZ EXT ',(igzrej(elemIndex),elemIndex=1,numElem)
730 write(*,222) ' REJ OUT BND ',(ibndrej(elemIndex),elemIndex=1,numElem)
731 write(*,222) ' REJ SBL ',(isblrej(elemIndex),elemIndex=1,numElem)
732 write(*,222) ' REJ TOTAL ',(itotrej(elemIndex),elemIndex=1,numElem)
733 write(*,*) '***************************************'
734 write(*,*) ' '
735 end if
736222 format(2x,a29,16(2x,i5))
737223 format(2x,a29,16(2x,f6.0))
738
739 countAssim=0
740 do bodyIndex=1,obs_numbody(obsSpaceData)
741 if ( obs_bodyElem_i(obsSpaceData,OBS_ASS,bodyIndex) == obs_assimilated) countAssim=countAssim+1
742 end do
743 if ( .not.beSilent ) write(*,'(1X," NUMBER OF DATA TO BE ASSIMILATED AFTER ADJUSTMENTS:",i10)') countAssim
744 if ( .not.beSilent ) write(*,*) ' '
745
746 end subroutine filt_topoRadiosonde
747
748 !--------------------------------------------------------------------------
749 ! filt_topoAISW
750 !--------------------------------------------------------------------------
751 subroutine filt_topoAISW(columnTrlOnTrlLev, obsSpaceData, obsFamily, beSilent)
752 !
753 !:Purpose: Refuse elements which are too close to the surface.
754 !
755 implicit none
756
757 ! Arguments:
758 type(struct_columnData), intent(in) :: columnTrlOnTrlLev
759 type(struct_obs), intent(inout) :: obsSpaceData
760 logical, intent(in) :: beSilent
761 character(len=2), intent(in) :: obsFamily
762
763 ! Locals:
764 integer :: headerIndex, bodyIndex, elemIndex, listIndex
765 integer :: ivnm, countRej(numElem), countAssim
766 real(8) :: obsPressure, pressureDiff
767
768 if (obsFamily /= 'AI' .and. obsFamily /= 'SW') then
769 call utl_abort('filt_topoAISW: only AI and SW family are handled by this routine. You ask for '//obsFamily)
770 end if
771
772 if ( .not.beSilent ) then
773 write(*,*) ' '
774 write(*,*) ' filt_topoAISW for obsFamily = ', obsFamily
775 write(*,*) ' '
776 write(*,*) '****************************************************'
777 write(*,222) 'ELEMENTS ', (elemList(elemIndex),elemIndex=1,numElem)
778 write(*,223) 'REJECTION BOUNDARY(HPA) ', (surfaceBufferZone_Pres,elemIndex=1,numElem)
779 write(*,*) '****************************************************'
780 write(*,*) ' '
781 end if
782
783 ! set counters to zero
784 countRej(:)=0
785
786 ! loop over all body indices of each family
787 call obs_set_current_body_list(obsSpaceData, obsFamily)
788 BODY: do
789 bodyIndex = obs_getBodyIndex(obsSpaceData)
790 if (bodyIndex < 0) exit BODY
791
792 ! skip this observation if already flagged to not assimilate
793 if(obs_bodyElem_i(obsSpaceData,OBS_ASS,bodyIndex) == obs_notAssimilated) cycle BODY
794
795 !
796 ! reject data too close to the model orography, put to
797 ! model orography, data which is below , but close to the surface.
798 !
799 obsPressure = obs_bodyElem_r(obsSpaceData,OBS_PPP,bodyIndex)
800 headerIndex = obs_bodyElem_i(obsSpaceData,OBS_HIND,bodyIndex)
801 pressureDiff = col_getElem(columnTrlOnTrlLev,1,headerIndex,'P0') - obsPressure
802 if ( pressureDiff < surfaceBufferZone_Pres ) then
803 ivnm=obs_bodyElem_i(obsSpaceData,OBS_VNM,bodyIndex)
804 listIndex = findElemIndex(ivnm)
805 if(listIndex == -1) cycle BODY
806 call obs_bodySet_i(obsSpaceData,OBS_ASS,bodyIndex,obs_notAssimilated)
807 countRej(listIndex)=countRej(listIndex)+1
808 call obs_bodySet_i(obsSpaceData,OBS_FLG,bodyIndex, &
809 ibset(obs_bodyElem_i(obsSpaceData,OBS_FLG,bodyIndex),18 ))
810 end if
811 end do BODY
812
813 if ( .not.beSilent ) then
814 write(*,*) ' '
815 write(*,*) '*****************************************************************'
816 write(*,*) ' FAMILY = ',obsFamily
817 write(*,222) 'ELEMENTS ',(elemList(elemIndex),elemIndex=1,numElem)
818 write(*,222) 'REJECTED ',(countRej(elemIndex),elemIndex=1,numElem)
819 write(*,*) '*****************************************************************'
820 write(*,*) ' '
821 end if
822
823 countAssim=0
824 do bodyIndex=1,obs_numbody(obsSpaceData)
825 if ( obs_bodyElem_i(obsSpaceData,OBS_ASS,bodyIndex) == obs_assimilated) countAssim=countAssim+1
826 end do
827 if ( .not.beSilent ) write(*,'(1X," NUMBER OF DATA TO BE ASSIMILATED AFTER ADJUSTMENTS:",i10)') countAssim
828 if ( .not.beSilent ) write(*,*) ' '
829
830222 format(2x,a29,16(2x,i5))
831223 format(2x,a29,16(2x,f5.0))
832
833end subroutine filt_topoAISW
834
835 !--------------------------------------------------------------------------
836 ! filt_topoProfiler
837 !--------------------------------------------------------------------------
838 subroutine filt_topoProfiler(columnTrlOnTrlLev,obsSpaceData,beSilent)
839 !
840 !:Purpose: Refuse elements which are too far away from the surface of the model
841 ! Refuse elements which are considered in the free atmosphere of
842 ! the RAOB but fall in the surface boundary layer of the model atmosphere.
843 !
844 implicit none
845
846 ! Arguments:
847 type(struct_columnData), intent(in) :: columnTrlOnTrlLev
848 type(struct_obs), intent(inout) :: obsSpaceData
849 logical, intent(in) :: beSilent
850
851 ! Locals:
852 integer :: headerIndex, bodyIndex, listIndex, elemIndex
853 integer :: ivnm, countAssim
854 integer :: itotrej(numElem), isblrej(numElem), ibndrej(numElem)
855 real(8) :: obsAltitude
856 real(8) :: obsSfcAltitude,colSfcAltitude,colAltitudeBelow,colAltitudeAbove
857 logical :: llok, list_is_empty
858
859 if ( .not.beSilent ) then
860 write(*,*) ' '
861 write(*,*) ' filt_topoProfiler: '
862 write(*,*) ' '
863 write(*,*) '************************************************'
864 write(*,222) ' ELEMENTS ',(elemList(elemIndex),elemIndex=1,numElem)
865 write(*,223) ' REJECTION BOUNDARY(METRE) ',(altDiffMax(elemIndex),elemIndex=1,numElem)
866 write(*,223) ' REJECTION SBL (METRE) ',(surfaceBufferZone_Height,elemIndex=1,numElem)
867 write(*,*) '************************************************'
868 write(*,*) ' '
869 end if
870
871 ! set counters to zero
872 itotrej(:)=0
873 isblrej(:)=0
874 ibndrej(:)=0
875
876 ! loop over all header indices of the 'PR' family
877 call obs_set_current_header_list(obsSpaceData, 'PR')
878 HEADER: do
879 headerIndex = obs_getHeaderIndex(obsSpaceData)
880 if (headerIndex < 0) exit HEADER
881
882 ! Set the body list & start at the beginning of the list
883 call obs_set_current_body_list(obsSpaceData, headerIndex,list_is_empty)
884 if (list_is_empty) cycle HEADER ! Proceed to the next HEADER
885
886 !
887 ! REJECT OBS BELOW THE MODEL SURFACE
888 ! AND THOSE NON SURFACE ELEMENTS PRESENT IN THE SURFACE
889 ! BOUNDARY LAYER OF THE RAOB OR OF THE MODEL.
890 ! AT THIS POINT WE WANT TO KEEP OBSERVATIONS IN THE FREE
891 ! ATMOSPHERE
892 !
893 colSfcAltitude = col_getHeight(columnTrlOnTrlLev,0,headerIndex,'SF')
894 obsSfcAltitude = obs_headElem_r(obsSpaceData,OBS_ALT,headerIndex)
895
896 ! loop over all body indices (still in the 'PR' family)
897 BODY: do
898 bodyIndex = obs_getBodyIndex(obsSpaceData)
899 if (bodyIndex < 0) exit BODY
900
901 ! skip this obs if already flagged to not be assimilated
902 if (obs_bodyElem_i(obsSpaceData,OBS_ASS,bodyIndex) == obs_notAssimilated) cycle BODY
903
904 ivnm = obs_bodyElem_i(obsSpaceData,OBS_VNM,bodyIndex)
905 listIndex = findElemIndex(ivnm)
906 llok = (obs_bodyElem_i(obsSpaceData,OBS_VCO,bodyIndex) == 1 &
907 .and. ivnm /= BUFR_NEGZ .and. listIndex /= -1)
908 if (.not. llok ) cycle BODY ! Proceed to the next bodyIndex
909
910 obsAltitude = obs_bodyElem_r(obsSpaceData,OBS_PPP,bodyIndex)
911 colAltitudeBelow = colSfcAltitude
912 if (obsSfcAltitude > colSfcAltitude) then
913 colAltitudeAbove = obsSfcAltitude + surfaceBufferZone_Height
914 else
915 colAltitudeAbove = colSfcAltitude + surfaceBufferZone_Height
916 end if
917 if(abs(obsSfcAltitude-colSfcAltitude) <= altDiffMax(listIndex)) then
918 !----Model surface and station altitude are very close
919 ! Accept observation if obsAltitude is within the domain
920 ! of the trial field
921 colAltitudeBelow = colSfcAltitude
922 colAltitudeAbove = col_getHeight(columnTrlOnTrlLev,col_getNumLev(columnTrlOnTrlLev,'MM')-1,headerIndex,'MM')
923 end if
924 if(obsAltitude < colAltitudeBelow ) then
925 call obs_bodySet_i(obsSpaceData,OBS_FLG,bodyIndex, &
926 ibset( obs_bodyElem_i(obsSpaceData,OBS_FLG,bodyIndex), 18 ))
927 call obs_bodySet_i(obsSpaceData,OBS_ASS,bodyIndex,obs_notAssimilated)
928 itotrej(listIndex)=itotrej(listIndex)+1
929 ibndrej(listIndex)=ibndrej(listIndex)+1
930 else if(obsAltitude >= colAltitudeBelow .and. obsAltitude < colAltitudeAbove ) then
931 call obs_bodySet_i(obsSpaceData,OBS_FLG,bodyIndex, &
932 ibset( obs_bodyElem_i(obsSpaceData,OBS_FLG,bodyIndex), 18 ))
933 call obs_bodySet_i(obsSpaceData,OBS_ASS,bodyIndex,obs_notAssimilated)
934 itotrej(listIndex)=itotrej(listIndex)+1
935 isblrej(listIndex)=isblrej(listIndex)+1
936 end if
937 end do BODY
938 end do HEADER
939
940 if ( .not.beSilent ) then
941 write(*,*) ' '
942 write(*,*) '***************************************'
943 write(*,*) 'FAMILY = PR'
944 write(*,222) ' ELEMENTS ', (elemList(elemIndex),elemIndex=1,numElem)
945 write(*,*) '************'
946 write(*,222) ' REJ OUT BND ',(ibndrej(elemIndex),elemIndex=1,numElem)
947 write(*,222) ' REJ SBL ',(isblrej(elemIndex),elemIndex=1,numElem)
948 write(*,222) ' REJ TOTAL ',(itotrej(elemIndex),elemIndex=1,numElem)
949 write(*,*) '***************************************'
950 write(*,*) ' '
951 end if
952222 format(2x,a29,16(2x,i5))
953223 format(2x,a29,16(2x,f6.0))
954
955 countAssim=0
956 do bodyIndex=1,obs_numbody(obsSpaceData)
957 if (obs_bodyElem_i(obsSpaceData,OBS_ASS,bodyIndex) == obs_assimilated) countAssim=countAssim+1
958 end do
959 if ( .not.beSilent ) write(*,'(1X," NUMBER OF DATA TO BE ASSIMILATED AFTER ADJUSTMENTS:",i10)') countAssim
960 if ( .not.beSilent ) write(*,*) ' '
961
962 end subroutine filt_topoProfiler
963
964 !--------------------------------------------------------------------------
965 ! filt_topoAladin
966 !--------------------------------------------------------------------------
967 subroutine filt_topoAladin(columnTrlOnTrlLev, obsSpaceData, beSilent)
968 !
969 !:Purpose: Refuse elements which are considered to be in the free atmosphere
970 ! of the Aladin instrument but which fall in the surface boundary
971 ! layer of the model atmosphere.
972 !
973 implicit none
974
975 ! Arguments:
976 type(struct_columnData), intent(in) :: columnTrlOnTrlLev
977 type(struct_obs), intent(inout) :: obsSpaceData
978 logical, intent(in) :: beSilent
979
980 ! Locals:
981 integer :: headerIndex, bodyIndex, elemIndex
982 integer :: ivnm, countAssim
983 integer :: countAcc(numElem), countRej(numElem)
984 real(8) :: obsAltitude ! altitide of the observation
985 real(8) :: colSfcAltitude ! altitude of the model's lowest layer
986 real(8) :: colAltitudeAbove ! top of the boundary layer
987 logical :: list_is_empty
988
989 if(.not. beSilent )then
990 write(*,*) ' '
991 write(*,*) ' filt_topoAladin: '
992 write(*,*) ' '
993 write(*,*) '************************************************'
994 write(*,222) ' ELEMENTS ',(elemList(elemIndex),elemIndex=1,numElem)
995 write(*,223) ' REJECTION SBL (METRE) ',(surfaceBufferZone_Height,elemIndex=1,numElem)
996 write(*,*) '************************************************'
997 write(*,*) ' '
998 end if
999
1000 ! set counter to zero
1001 countAcc(:)=0
1002 countRej(:)=0
1003
1004 ! loop over all header indices of the 'AL' family
1005 call obs_set_current_header_list(obsSpaceData, 'AL')
1006 HEADER: do
1007 headerIndex = obs_getHeaderIndex(obsSpaceData)
1008 if (headerIndex < 0) exit HEADER
1009
1010 ! Set the body list & start at the beginning of the list
1011 call obs_set_current_body_list(obsSpaceData, headerIndex,list_is_empty)
1012 if (list_is_empty) cycle HEADER ! Proceed to the next HEADER
1013
1014 !
1015 ! REJECT OBS IN THE SURFACE BOUNDARY LAYER OF THE MODEL.
1016 ! AT THIS POINT WE WANT TO KEEP OBSERVATIONS THAT ARE IN THE FREE
1017 ! ATMOSPHERE
1018 !
1019 colSfcAltitude = col_getHeight(columnTrlOnTrlLev,col_getNumLev(columnTrlOnTrlLev,'MM'), &
1020 headerIndex,'MM')
1021 colAltitudeAbove = colSfcAltitude + surfaceBufferZone_Height
1022
1023 ! Loop over all body indices (still in the 'AL' family)
1024 BODY: do
1025 bodyIndex = obs_getBodyIndex(obsSpaceData)
1026 if (bodyIndex < 0) exit BODY
1027
1028 ! Skip this obs if it is not on height levels
1029 if (obs_bodyElem_i(obsSpaceData,OBS_VCO,bodyIndex) /= 1) cycle BODY
1030
1031 ! Skip this obs if already flagged not to be assimilated
1032 if (obs_bodyElem_i(obsSpaceData,OBS_ASS,bodyIndex) == obs_notAssimilated) cycle BODY
1033
1034 ! Skip this obs if it is not in the list of element types
1035 ivnm = obs_bodyElem_i(obsSpaceData,OBS_VNM,bodyIndex)
1036 elemIndex = findElemIndex(ivnm)
1037 if(elemIndex == -1) cycle BODY
1038
1039
1040 !
1041 ! apply filter to selected elements
1042 !
1043
1044 ! Reject this obs if it is in the boundary layer or below
1045 obsAltitude = obs_bodyElem_r(obsSpaceData,OBS_PPP,bodyIndex)
1046 if(obsAltitude > colAltitudeAbove) then
1047 ! obs passes the acceptance criterion
1048 countAcc(elemIndex) = countAcc(elemIndex)+1
1049
1050 else
1051 ! Flag rejection due to orography
1052 call obs_bodySet_i(obsSpaceData,OBS_FLG,bodyIndex, &
1053 ibset( obs_bodyElem_i(obsSpaceData,OBS_FLG,bodyIndex), 18 ))
1054
1055 ! Do not assimilate the observation
1056 call obs_bodySet_i(obsSpaceData,OBS_ASS,bodyIndex,obs_notAssimilated)
1057 countRej(elemIndex)=countRej(elemIndex)+1
1058 end if
1059 end do BODY
1060 end do HEADER
1061
1062 if(.not. besilent)then
1063 write(*,*) ' '
1064 write(*,*) '***************************************'
1065 write(*,*) 'FAMILY = AL'
1066 write(*,222) ' ELEMENTS ', (elemList(elemIndex),elemIndex=1,numElem)
1067 write(*,222) ' ACCEPTED ', (countAcc(elemIndex),elemIndex=1,numElem)
1068 write(*,222) ' REJECTED ', (countRej(elemIndex),elemIndex=1,numElem)
1069 write(*,*) '***************************************'
1070 write(*,*) ' '
1071 end if
1072222 format(2x,a29,16(2x,i5))
1073223 format(2x,a29,16(2x,f6.0))
1074
1075 countAssim=0
1076 do bodyIndex=1,obs_numbody(obsSpaceData)
1077 if (obs_bodyElem_i(obsSpaceData,OBS_ASS,bodyIndex) == obs_assimilated) countAssim=countAssim+1
1078 end do
1079
1080 if(.not. besilent)then
1081 write(*,'(1X," NUMBER OF DATA TO BE ASSIMILATED AFTER ADJUSTMENTS:",i10)') countAssim
1082 write(*,*) ' '
1083 end if
1084
1085 end subroutine filt_topoAladin
1086
1087 !--------------------------------------------------------------------------
1088 ! filt_topoTovs
1089 !--------------------------------------------------------------------------
1090 subroutine filt_topoTovs(columnTrlOnTrlLev, obsSpaceData, beSilent)
1091 !
1092 !:Purpose: Refuse data which are too close to the surface.
1093 !
1094 implicit none
1095
1096 ! Arguments:
1097 type(struct_columnData), intent(in) :: columnTrlOnTrlLev
1098 type(struct_obs), intent(inout) :: obsSpaceData
1099 logical, intent(in) :: beSilent
1100
1101 ! Locals:
1102 integer :: headerIndex, bodyIndex
1103 integer :: idatyp, countAssim, countRej
1104 real(8), parameter :: minSfcPressure = 80000.d0
1105
1106 if ( .not.beSilent ) then
1107 write(*,* ) ' '
1108 write(*,* ) ' filt_topoTovs: '
1109 write(*,* ) ' '
1110 write(*,* ) '****************************************************'
1111 write(*,222) 'ELEMENTS ', BUFR_NBT3
1112 write(*,223) 'MINIMUM SFC PRESSURE (PA) ', minSfcPressure
1113 write(*,* ) '****************************************************'
1114 write(*,* ) ' '
1115 end if
1116
1117 ! set counters to zero
1118 countRej=0
1119
1120 ! loop over all header indices of the 'TO' family
1121 call obs_set_current_header_list(obsSpaceData, 'TO')
1122 HEADER: do
1123 headerIndex = obs_getHeaderIndex(obsSpaceData)
1124 if (headerIndex < 0) exit HEADER
1125
1126 idatyp = obs_headElem_i(obsSpaceData,OBS_ITY,headerIndex)
1127 if (idatyp .ne. 185) cycle HEADER ! Proceed to the next headerIndex
1128
1129 ! loop over all body indices (still in the 'TO' family)
1130 call obs_set_current_body_list(obsSpaceData, headerIndex)
1131 BODY: do
1132 bodyIndex = obs_getBodyIndex(obsSpaceData)
1133 if (bodyIndex < 0) exit BODY
1134
1135 if (obs_bodyElem_i(obsSpaceData,OBS_VNM,bodyIndex).ne.BUFR_NBT3) cycle BODY
1136
1137 ! reject obs if the model surface pressure is below the minimum specified value
1138 if (col_getElem(columnTrlOnTrlLev,1,headerIndex,'P0') < minSfcPressure) then
1139 call obs_bodySet_i(obsSpaceData,OBS_ASS,bodyIndex,obs_notAssimilated)
1140 countRej=countRej+1
1141 call obs_bodySet_i(obsSpaceData,OBS_FLG,bodyIndex, &
1142 ibset(obs_bodyElem_i(obsSpaceData,OBS_FLG,bodyIndex),9))
1143 call obs_bodySet_i(obsSpaceData,OBS_FLG,bodyIndex, &
1144 ibset(obs_bodyElem_i(obsSpaceData,OBS_FLG,bodyIndex),18))
1145 end if
1146
1147 end do BODY
1148 end do HEADER
1149
1150 if ( .not.beSilent ) then
1151 write(*,*) ' '
1152 write(*,*) '*****************************************************************'
1153 write(*,*) ' FAMILY = TO'
1154 write(*,222) 'ELEMENTS ', BUFR_NBT3
1155 write(*,222) 'REJECTED ',countRej
1156 write(*,*) '*****************************************************************'
1157 write(*,*) ' '
1158 end if
1159222 format(2x,a29,1(4x,i5))
1160223 format(2x,a29,1(2x,f7.0))
1161
1162 countAssim=0
1163 do bodyIndex=1,obs_numbody(obsSpaceData)
1164 if ( obs_bodyElem_i(obsSpaceData,OBS_ASS,bodyIndex) == obs_assimilated) countAssim=countAssim+1
1165 end do
1166 if ( .not.beSilent ) write(*,'(1X," NUMBER OF DATA TO BE ASSIMILATED AFTER ADJUSTMENTS:",i10)') countAssim
1167 if ( .not.beSilent ) write(*,* ) ' '
1168
1169 end subroutine filt_topoTovs
1170
1171 !--------------------------------------------------------------------------
1172 ! filt_surfaceWind
1173 !--------------------------------------------------------------------------
1174 subroutine filt_surfaceWind(obsSpaceData, beSilent)
1175 !
1176 !:Purpose: zap sfc wind components at land stations
1177 !
1178 IMPLICIT NONE
1179
1180 ! Arguments:
1181 type(struct_obs), intent(inout) :: obsSpaceData
1182 logical, intent(in) :: beSilent
1183
1184 ! Locals:
1185 INTEGER, parameter :: JPINEL=2,JPIDLND=9
1186 INTEGER :: J,JID,JDATA
1187 LOGICAL :: LLPRINT
1188 INTEGER :: ITYP,IDBURP
1189 INTEGER :: ILISTEL(JPINEL), IDLND(JPIDLND)
1190 INTEGER :: IKOUNTREJ(JPINEL), IKOUNTT
1191 character(len=2), dimension(2) :: list_family
1192 integer :: index_family, index_header, index_body
1193 character(len=12) :: stnid
1194 real(pre_obsReal) :: obsLAT, obsLON, obsPPP
1195
1196 DATA IDLND / 12, 14, 146, 32, 35, 135, 136, 137, 138 /
1197 !
1198 if ( .not. discardlandsfcwind ) return
1199 !
1200 ILISTEL(1)=BUFR_NEUS
1201 ILISTEL(2)=BUFR_NEVS
1202 if ( .not.beSilent ) then
1203 WRITE(*,* ) ' '
1204 WRITE(*,* ) ' filt_surfaceWind:'
1205 WRITE(*,* ) ' '
1206 WRITE(*,* ) '*****************************************************'
1207 WRITE(*,222)'ELEMENTS REJECTED ',( ILISTEL(J),J=1,jpinel)
1208 WRITE(*,222)'LIST OF IDTYP ',( idlnd(J),J=1,jpidlnd)
1209 WRITE(*,* ) '*****************************************************'
1210 WRITE(*,* ) ' '
1211 end if
1212 LLPRINT = .FALSE.
1213 !cc LLPRINT = .TRUE.
1214 !
1215 ! SET COUNTERS TO ZERO
1216 !
1217 DO J=1,JPINEL
1218 IKOUNTREJ(J)=0
1219 END DO
1220
1221 !
1222 ! Loop over the families of interest
1223 !
1224 list_family(1) = 'SF'
1225 list_family(2) = 'UA'
1226 do index_family = 1,2
1227 if ( .not.beSilent ) WRITE(*,'(2x,A9,2x,A2)')'FAMILY = ',list_family(index_family)
1228
1229 !
1230 ! loop over all header indices of each family
1231 !
1232 ! Set the header list
1233 ! (& start at the beginning of the list)
1234 call obs_set_current_header_list(obsSpaceData, &
1235 list_family(index_family))
1236 HEADER: do
1237 index_header = obs_getHeaderIndex(obsSpaceData)
1238 if (index_header < 0) exit HEADER
1239
1240 !
1241 ! loop over all body indices (still in the same family)
1242 !
1243 ! Set the body list
1244 ! (& start at the beginning of the list)
1245 call obs_set_current_body_list(obsSpaceData, index_header)
1246 BODY: do
1247 index_body = obs_getBodyIndex(obsSpaceData)
1248 if (index_body < 0) exit BODY
1249
1250 ! UNCONDITIONALLY REJECT SURFACE WINDS AT SYNOP/TEMP LAND STATIONS
1251 ITYP=obs_bodyElem_i(obsSpaceData,OBS_VNM,INDEX_BODY)
1252 IDBURP = obs_headElem_i(obsSpaceData,OBS_ITY,INDEX_HEADER)
1253 IF ( ITYP == BUFR_NEUS .OR. ITYP == BUFR_NEVS) THEN
1254 DO JID = 1, JPIDLND
1255 IF(IDBURP == IDLND(JID) .AND. &
1256 obs_bodyElem_i(obsSpaceData,OBS_ASS,INDEX_BODY) == obs_assimilated) THEN
1257 call obs_bodySet_i(obsSpaceData,OBS_FLG,INDEX_BODY, &
1258 ibset( obs_bodyElem_i(obsSpaceData,OBS_FLG,INDEX_BODY), 19))
1259 call obs_bodySet_i(obsSpaceData,OBS_ASS,INDEX_BODY,obs_notAssimilated)
1260 DO J = 1, JPINEL
1261 IF(ITYP ==ILISTEL(J)) THEN
1262 IKOUNTREJ(J)=IKOUNTREJ(J)+1
1263 END IF
1264 END DO
1265 IF(LLPRINT .and. .not.beSilent ) THEN
1266 stnid = obs_elem_c(obsSpaceData,'STID',INDEX_HEADER)
1267 obsLAT = obs_headElem_r(obsSpaceData,OBS_LAT,INDEX_HEADER)
1268 obsLON = obs_headElem_r(obsSpaceData,OBS_LON,INDEX_HEADER)
1269 obsPPP = obs_bodyElem_r(obsSpaceData,OBS_PPP,INDEX_BODY)
1270 WRITE(*,225) 'Rej sfc wind lnd',INDEX_HEADER,ITYP,stnid, &
1271 IDBURP, obsLAT, obsLON, obsPPP
1272 END IF
1273 END IF
1274 END DO
1275 END IF ! BUFR_NEUS or BUFR_NEVS
1276 END DO BODY
1277 END DO HEADER
1278 !
1279 if ( .not.beSilent ) then
1280 WRITE(*,* ) ' '
1281 WRITE(*,* ) '*****************************************************'
1282 WRITE(*,222 )'ELEMENTS ', ( ILISTEL(J),J=1,JPINEL)
1283 WRITE(*,222)'REJECTED ',(IKOUNTREJ(J),J=1,JPINEL)
1284 WRITE(*,* ) '*****************************************************'
1285 WRITE(*,* ) ' '
1286 END IF
1287222 FORMAT(2x,a29,10(2x,i5))
1288223 FORMAT(2x,a29,10(2x,f5.0))
1289224 FORMAT(2x,a17,2x,I6,2X,I5,1x,a9,1x,2(2x,f9.2))
1290225 FORMAT(2x,a13,2x,I6,2X,I5,1x,a9,1x,I6,1x,3(2x,f9.2))
1291 !
1292 END DO ! family
1293 !
1294 IKOUNTT=0
1295 DO JDATA=1,obs_numbody(obsSpaceData)
1296 IF ( obs_bodyElem_i(obsSpaceData,OBS_ASS,JDATA) == obs_assimilated) IKOUNTT=IKOUNTT+1
1297 END DO
1298 if ( .not.beSilent ) WRITE(*, &
1299 '(1X," NUMBER OF DATA ASSIMILATED BY MIDAS AFTER ADJUSTMENTS: ",i10)') &
1300 IKOUNTT
1301 if ( .not.beSilent ) WRITE(*,* ) ' '
1302
1303 end subroutine filt_surfaceWind
1304
1305 !--------------------------------------------------------------------------
1306 ! filt_radvel
1307 !--------------------------------------------------------------------------
1308 subroutine filt_radvel(columnTrlOnTrlLev, obsSpaceData, beSilent)
1309 !
1310 !:Purpose: Filter Radvel observations
1311 ! guarantee that altitude values are
1312 ! within bounds Altitude and check the horizontal distance between levels
1313 ! for further processing
1314 implicit none
1315
1316 ! Arguments:
1317 type(struct_columnData), intent(in) :: columnTrlOnTrlLev
1318 type(struct_obs) , intent(inout) :: obsSpaceData
1319 logical , intent(in) :: beSilent
1320
1321 ! Locals:
1322 integer :: bodyIndex, headerIndex, numLevels, bufrCode, obsflag
1323 integer :: fnom, fclos, nulnam, ierr, levelIndex
1324 real(8) :: obsAltitude, radarAltitude, beamElevation
1325 real(8) :: levelAltLow, levelAltHigh
1326 real(8) :: levelBracketLow, levelBracketHigh
1327 real(8) :: levelRangeNear, levelRangeFar
1328 logical, save :: firstCall = .true.
1329
1330 ! Namelist variables:
1331 real(8), save :: maxRangeInterp ! max allowable horizontal distance between levels (in m) for radar winds
1332
1333 namelist /namradvel/ maxRangeInterp
1334
1335 if (.not.beSilent) then
1336 write(*,*)
1337 write(*,*) 'filt_radvel: begin'
1338 end if
1339
1340 ! reading namelist variables
1341 if (firstCall) then
1342 ! default value
1343 maxRangeInterp = -1.0D0
1344
1345 if ( utl_isNamelistPresent('namradvel', './flnml') ) then
1346 nulnam=0
1347 ierr=fnom(nulnam, './flnml', 'FTN+SEQ+R/O', 0)
1348 read(nulnam, nml=namradvel, iostat=ierr)
1349 if ( ierr /= 0 ) call utl_abort('oop_raDvel_nl: Error reading namelist namradvel')
1350 if ( .not. beSilent ) write(*,nml=namradvel)
1351 ierr = fclos(nulnam)
1352 else if ( .not. beSilent ) then
1353 write(*,*)
1354 write(*,*) 'filt_radvel: namradvel is missing in the namelist. The default value will be taken.'
1355 end if
1356 firstCall = .false.
1357 end if
1358 !
1359 ! Loop over all header indices of the 'RA' family (Doppler Velocity)
1360 !
1361 call obs_set_current_header_list(obsSpaceData, 'RA')
1362 HEADER: do
1363 headerIndex = obs_getHeaderIndex(obsSpaceData)
1364 if ( headerIndex < 0 ) exit HEADER
1365 !
1366 numLevels = col_getNumLev(columnTrlOnTrlLev, 'MM')
1367 ! Elevation beam (PPI)
1368 beamElevation = obs_headElem_r(obsSpaceData, OBS_RELE, headerIndex)
1369 ! Altitude radar
1370 radarAltitude = obs_headElem_r(obsSpaceData, OBS_ALT, headerIndex)
1371 !
1372 ! Loop over all body indices of the 'RA' family (Doppler Velocity)
1373 !
1374 call obs_set_current_body_list(obsSpaceData, headerIndex)
1375 BODY: do
1376 bodyIndex = obs_getBodyIndex(obsSpaceData)
1377 if ( bodyIndex < 0 ) exit BODY
1378 ! Check that this observation has the expected bufr element ID
1379 bufrCode = obs_bodyElem_i(obsSpaceData, OBS_VNM, bodyIndex)
1380 if ( bufrCode /= bufr_radvel ) cycle BODY
1381 ! only process observations flagged to be assimilated
1382 if ( obs_bodyElem_i(obsSpaceData, OBS_ASS, bodyIndex) /= obs_assimilated ) cycle BODY
1383
1384 ! Altitude of observation
1385 obsAltitude = obs_bodyElem_r(obsSpaceData, OBS_PPP, bodyIndex)
1386
1387 ! TODO we need to understand why model height at "numLevels" is not always the lowest
1388 ! Observations are rejected if their altitude is below the height 1 which may not be lowest model level...
1389 levelAltLow = col_getHeight(columnTrlOnTrlLev, numLevels, headerIndex,'MM')
1390 if ( obsAltitude < levelAltLow ) then
1391 call obs_bodySet_i(obsSpaceData, OBS_ASS, bodyindex, obs_notAssimilated)
1392 obsFlag = obs_bodyElem_i(obsSpaceData, OBS_FLG, bodyindex)
1393 call obs_bodySet_i(obsSpaceData, OBS_FLG, bodyindex, IBSET(obsFlag, 11))
1394 cycle BODY
1395 end if
1396
1397 ! Observations are rejected if observation is above the height of first (highest) model level
1398 levelAltHigh = col_getHeight(columnTrlOnTrlLev, 1, headerIndex,'MM')
1399 if ( obsAltitude > levelAltHigh ) then
1400 call obs_bodySet_i(obsSpaceData, OBS_ASS, bodyindex, obs_notAssimilated)
1401 obsFlag = obs_bodyElem_i(obsSpaceData, OBS_FLG, bodyindex)
1402 call obs_bodySet_i(obsSpaceData, OBS_FLG, bodyindex, IBSET(obsFlag, 11))
1403 cycle BODY
1404 end if
1405
1406 ! Levels that bracket the observation from OBS_LYR
1407 ! note to self: like in GEM, level=1 is the highest level
1408 levelIndex = obs_bodyElem_i(obsSpaceData, OBS_LYR, bodyIndex)
1409 levelBracketHigh = col_getHeight(columnTrlOnTrlLev, levelIndex, headerIndex,'MM')
1410 levelBracketLow = col_getHeight(columnTrlOnTrlLev, levelIndex+1, headerIndex,'MM')
1411
1412 ! Observations are rejected if horizontal distance between levels is too large
1413 if ( maxRangeInterp > 0.0 ) then
1414 call rdv_getRangefromH(levelBracketHigh, radarAltitude, beamElevation, levelRangeFar )
1415 call rdv_getRangefromH(levelBracketLow, radarAltitude, beamElevation, levelRangeNear)
1416
1417 if ( abs(levelRangeFar-levelRangeNear) > maxRangeInterp ) then
1418 call obs_bodySet_i(obsSpaceData, OBS_ASS, bodyindex, obs_notAssimilated)
1419 obsFlag = obs_bodyElem_i(obsSpaceData, OBS_FLG, bodyindex)
1420 call obs_bodySet_i(obsSpaceData, OBS_FLG, bodyindex, IBSET(obsFlag, 11))
1421 cycle BODY
1422 end if
1423 end if
1424
1425 end do BODY
1426 end do HEADER
1427 if ( .not. beSilent ) write(*,*) 'filt_radvel: end'
1428 end subroutine filt_radvel
1429
1430 ! filt_gpsro
1431 !--------------------------------------------------------------------------
1432 subroutine FILT_GPSRO(columnTrlOnTrlLev, obsSpaceData, beSilent)
1433 !
1434 !:Purpose: Filter GPSRO observations
1435 ! Guarantee that altitude and observation values are
1436 ! within bounds for further processing
1437 !
1438 !:Note: For noncompliant GPSRO observations:
1439 !
1440 ! - Set assimilable flag to 0
1441 ! - Set bit of cma flag 11 ON
1442 !
1443 use gps_mod
1444 implicit none
1445
1446 ! Arguments:
1447 type(struct_columnData), intent(in) :: columnTrlOnTrlLev
1448 type(struct_obs) , intent(inout) :: obsSpaceData
1449 logical , intent(in) :: beSilent
1450
1451 ! Locals:
1452 INTEGER :: INDEX_HEADER, IDATYP, INDEX_BODY
1453 INTEGER :: JL, ISAT, IQLF, iProfile, IFLG, varNum, IDSC
1454 REAL(8) :: ZMT, Rad, Geo, AZM
1455 REAL(8) :: HNH1, HSF, HTP, HMIN, HMAX, ZOBS, ZREF, ZSAT
1456 LOGICAL :: LLEV, LOBS, LNOM, LSAT, LAZM, LALL, LDSC
1457
1458 if (.not.beSilent) then
1459 write(*,*)
1460 write(*,*) 'filt_gpsro: begin'
1461 end if
1462 !
1463 ! Loop over all header indices of the 'RO' family:
1464 !
1465 call obs_set_current_header_list(obsSpaceData,'RO')
1466 gps_numROProfiles=0
1467 HEADER: do
1468 index_header = obs_getHeaderIndex(obsSpaceData)
1469 if (index_header < 0) exit HEADER
1470 !
1471 ! Process only refractivity data (codtyp 169):
1472 !
1473 IDATYP = obs_headElem_i(obsSpaceData,OBS_ITY,INDEX_HEADER)
1474 if ( IDATYP == 169 ) then
1475 gps_numROProfiles=gps_numROProfiles+1
1476 !
1477 ! Basic variables of the profile:
1478 !
1479 AZM = obs_headElem_r(obsSpaceData,OBS_AZA ,INDEX_HEADER)
1480 ISAT = obs_headElem_i(obsSpaceData,OBS_SAT ,INDEX_HEADER)
1481 IQLF = obs_headElem_i(obsSpaceData,OBS_ROQF,INDEX_HEADER)
1482 Rad = obs_headElem_r(obsSpaceData,OBS_TRAD,INDEX_HEADER)
1483 Geo = obs_headElem_r(obsSpaceData,OBS_GEOI,INDEX_HEADER)
1484 LNOM = .NOT.BTEST(IQLF,16-1)
1485 LAZM = .TRUE.
1486 !
1487 ! Check if the satellite is within the accepted set:
1488 !
1489 ZSAT = ABS(gps_WGPS(ISAT,1))+ABS(gps_WGPS(ISAT,2))+ABS(gps_WGPS(ISAT,3))+ABS(gps_WGPS(ISAT,4))
1490 LSAT = ( ZSAT > 0.d0 )
1491 !
1492 ZMT = col_getHeight(columnTrlOnTrlLev,0,index_header,'SF')
1493 !
1494 ! Acceptable height limits:
1495 !
1496 JL = 1
1497 HTP = col_getHeight(columnTrlOnTrlLev,JL,INDEX_HEADER,'TH')
1498 HSF = ZMT+gps_SurfMin
1499 IF (HSF < gps_HsfMin) HSF=gps_HsfMin
1500 IF (HTP > gps_HtpMax) HTP=gps_HtpMax
1501 HMIN=Geo+HSF
1502 HMAX=Geo+HTP
1503 !
1504 ! Loop over all body indices for this index_header:
1505 ! (start at the beginning of the list)
1506 !
1507 call obs_set_current_body_list(obsSpaceData, INDEX_HEADER)
1508 BODY: do
1509 index_body = obs_getBodyIndex(obsSpaceData)
1510 if (index_body < 0) exit BODY
1511 !
1512 ! Altitude and reference order of magnitude value:
1513 !
1514 HNH1= obs_bodyElem_r(obsSpaceData,OBS_PPP,INDEX_BODY)
1515 varNum = obs_bodyElem_i(obsSpaceData,OBS_VNM,INDEX_BODY)
1516 if (varNum == bufr_nebd) then
1517 HNH1=HNH1-Rad
1518 ZREF = 0.025d0*exp(-HNH1/6500.d0)
1519 LAZM = (-0.1d0 < AZM .AND. AZM < 360.1)
1520 else
1521 ZREF = 300.d0*exp(-HNH1/6500.d0)
1522 LAZM = .TRUE.
1523 end if
1524 !
1525 ! Observation:
1526 !
1527 ZOBS= obs_bodyElem_r(obsSpaceData,OBS_VAR,INDEX_BODY)
1528 !
1529 ! Positively verify that the altitude is within bounds:
1530 !
1531 LLEV= (HNH1 > HMIN) .AND. (HNH1 < HMAX)
1532 !
1533 ! Positively verify that the observable is within bounds:
1534 !
1535 LOBS= (ZOBS > (0.3d0*ZREF)) .AND. (ZOBS < (3.d0*ZREF))
1536 !
1537 ! Mark as not assimilable unless all conditions are satisfied:
1538 !
1539 LALL = LLEV .AND. LOBS .AND. LAZM .AND. LNOM .AND. LSAT
1540 if ( .NOT.LALL ) then
1541 call obs_bodySet_i(obsSpaceData,OBS_ASS,INDEX_BODY, obs_notAssimilated)
1542 IFLG = obs_bodyElem_i(obsSpaceData,OBS_FLG,INDEX_BODY)
1543 call obs_bodySet_i(obsSpaceData,OBS_FLG,INDEX_BODY, IBSET(IFLG,11))
1544 end if
1545 ! Do not assimilate bending in mode gps_Level_RO = gps_Level_RO_Ref:
1546 if (varNum == bufr_nebd .and. gps_Level_RO == gps_Level_RO_Ref) then
1547 call obs_bodySet_i(obsSpaceData,OBS_ASS,INDEX_BODY, obs_notAssimilated)
1548 endif
1549 ! Do not assimilate refractivity in mode gps_Level_RO = gps_Level_RO_Bnd:
1550 if (varNum == bufr_nerf .and. gps_Level_RO == gps_Level_RO_Bnd) then
1551 call obs_bodySet_i(obsSpaceData,OBS_ASS,INDEX_BODY, obs_notAssimilated)
1552 endif
1553 end do BODY
1554 end if
1555 end do HEADER
1556 !
1557 ! List to enumerate and cross-link GPSRO headers 0 <= iProfile < gps_numROProfiles):
1558 !
1559 if (gps_numROProfiles > 0) then
1560 if(.not.allocated(gps_vRO_IndexPrf)) allocate(gps_vRO_IndexPrf(gps_numROProfiles, 10))
1561 iProfile=0
1562 !
1563 ! Loop over all header indices of the 'RO' family:
1564 !
1565 call obs_set_current_header_list(obsSpaceData,'RO')
1566 HEADER2: do
1567 index_header = obs_getHeaderIndex(obsSpaceData)
1568 if (index_header < 0) exit HEADER2
1569 !
1570 ! Process only refractivity data (codtyp 169):
1571 !
1572 IDATYP = obs_headElem_i(obsSpaceData,OBS_ITY,INDEX_HEADER)
1573 if ( IDATYP == 169 ) then
1574 iProfile=iProfile+1
1575 gps_vRO_IndexPrf(iProfile, 1) = INDEX_HEADER
1576 ISAT = obs_headElem_i(obsSpaceData,OBS_SAT ,INDEX_HEADER)
1577 IQLF = obs_headElem_i(obsSpaceData,OBS_ROQF,INDEX_HEADER)
1578 LDSC = .not.btest(IQLF,16-3)
1579 if (LDSC) then
1580 IDSC = 0
1581 else
1582 IDSC = 1
1583 end if
1584 varNum = -1
1585 call obs_set_current_body_list(obsSpaceData, INDEX_HEADER)
1586 !
1587 ! Loop over all body indices
1588 ! For storing varNum of each profile for the 'RO' family
1589 !
1590 BODY2: do
1591 index_body = obs_getBodyIndex(obsSpaceData)
1592 if (index_body < 0) exit BODY2
1593 varNum = obs_bodyElem_i(obsSpaceData,OBS_VNM,INDEX_BODY)
1594 if (varNum > 0) exit BODY2
1595 end do BODY2
1596 gps_vRO_IndexPrf(iProfile, 2) = varNum
1597 gps_vRO_IndexPrf(iProfile, 3) = ISAT
1598 gps_vRO_IndexPrf(iProfile, 4) = IDSC
1599 if (.not.beSilent) write(*,*)'RO Prf', gps_numROProfiles, iProfile, varNum, ISAT, IDSC
1600 end if
1601 end do HEADER2
1602 end if
1603
1604 if (.not.beSilent) write(*,*) 'filt_gpsro: end'
1605 end subroutine FILT_GPSRO
1606
1607 !--------------------------------------------------------------------------
1608 ! filt_backScatAnisIce
1609 !--------------------------------------------------------------------------
1610 subroutine filt_backScatAnisIce(obsSpaceData, beSilent)
1611 !
1612 !:Purpose: Filter scatterometer backscatter anisotropy observations
1613 ! where wind speed is too small
1614 !
1615 !:Note: For noncompliant observations:
1616 !
1617 ! - Set assimilable flag to 0
1618 ! - Set bit of cma flag 13 ON
1619 !
1620 implicit none
1621
1622 ! Arguments:
1623 type(struct_obs), intent(inout) :: obsSpaceData
1624 logical, intent(in) :: beSilent
1625
1626 ! Locals:
1627 integer :: bufrCode, headerIndex, bodyIndex
1628 real(8) :: modelWindSpeed
1629
1630 if (.not. obs_famExist(obsSpaceData,'GL')) return
1631
1632 if (.not. beSilent) then
1633 write(*,*)
1634 write(*,*) ' filt_backScatAnisIce: begin'
1635 end if
1636
1637 ! loop over all body indices
1638 call obs_set_current_body_list( obsSpaceData, 'GL' )
1639
1640 BODY: do
1641
1642 bodyIndex = obs_getBodyIndex( obsSpaceData )
1643 if ( bodyIndex < 0 ) exit BODY
1644
1645 bufrCode = obs_bodyElem_i( obsSpaceData, OBS_VNM, bodyIndex )
1646
1647 if ( bufrCode == BUFR_ICES ) then
1648
1649 headerIndex = obs_bodyElem_i( obsSpaceData, OBS_HIND, bodyIndex )
1650 modelWindSpeed = obs_headElem_r( obsSpaceData, OBS_MWS, headerIndex )
1651
1652 if ( modelWindSpeed < 4.0 ) then
1653 call obs_bodySet_i(obsSpaceData, OBS_ASS, bodyIndex, obs_notAssimilated)
1654 call obs_bodySet_i(obsSpaceData, OBS_FLG, bodyIndex, IBSET(obs_bodyElem_i(obsSpaceData,OBS_FLG,bodyIndex),13))
1655 end if
1656
1657 else
1658
1659 cycle BODY
1660
1661 end if
1662
1663 end do BODY
1664
1665 if (.not. beSilent) write(*,*) ' filt_backScatAnisIce: end'
1666
1667 end subroutine filt_backScatAnisIce
1668
1669 !--------------------------------------------------------------------------
1670 ! filt_iceConcentration
1671 !--------------------------------------------------------------------------
1672 subroutine filt_iceConcentration(obsSpaceData, beSilent)
1673 !
1674 !:Purpose: Filter out observations from satellites
1675 ! not specified in the name list
1676 !
1677 !:Note: For noncompliant observations:
1678 !
1679 ! - Set assimilable flag to 0
1680 !
1681 implicit none
1682
1683 ! Arguments:
1684 type(struct_obs), intent(inout) :: obsSpaceData
1685 logical, intent(in) :: beSilent
1686
1687 ! Locals:
1688 character(len=12) :: cstnid
1689 integer :: nulnam, ierr
1690 integer :: headerIndex, bodyIndex, codeType, platformIndex
1691 integer :: fnom, fclos
1692 logical :: inPlatformList
1693 logical, save :: firstCall = .true.
1694 integer, parameter :: maxPlatformIce = 50
1695
1696 ! Namelist variables:
1697 integer, save :: nPlatformIce ! MUST NOT BE INCLUDED IN NAMELIST!
1698 character(len=12), save :: listPlatformIce(maxPlatformIce) ! list of ice obs 'platforms' (station IDs) to assimilate
1699
1700 namelist /namPlatformIce/ nPlatformIce, listPlatformIce
1701
1702 if (.not. obs_famExist(obsSpaceData,'GL')) return
1703
1704 if (firstCall) then
1705 ! set default values for namelist variables
1706 nPlatformIce = MPC_missingValue_INT
1707 listPlatformIce(:) = '1234567890ab'
1708
1709 if (utl_isNamelistPresent('namPlatformIce','./flnml')) then
1710 nulnam = 0
1711 ierr = fnom(nulnam,'./flnml','FTN+SEQ+R/O',0)
1712 read (nulnam, nml = namPlatformIce, iostat = ierr)
1713 if ( ierr /= 0 ) call utl_abort('filt_iceConcentration: Error reading namelist')
1714 if ( mmpi_myid == 0 ) write(*,nml=namPlatformIce)
1715 ierr = fclos(nulnam)
1716 if (nPlatformIce /= MPC_missingValue_INT) then
1717 call utl_abort('filt_iceConcentration: check namPlatformIce namelist section: nPlatformIce should be removed')
1718 end if
1719 nPlatformIce = 0
1720 do platformIndex = 1, maxPlatformIce
1721 if (listPlatformIce(platformIndex) == '1234567890ab') exit
1722 nPlatformIce = nPlatformIce + 1
1723 end do
1724 else
1725 write(*,*)
1726 write(*,*) 'filt_iceConcentration: namPlatformIce is missing in the namelist. Filtering will be skipped.'
1727 end if
1728 firstCall = .false.
1729 end if
1730
1731 if ( nPlatformIce < 1 ) return
1732
1733 if ( nPlatformIce > maxPlatformIce ) then
1734 call utl_abort('filt_setup: too many elements for listPlatformIce')
1735 end if
1736
1737 if (.not. beSilent) then
1738 write(*,*)
1739 write(*,*) 'filt_iceConcentration: begin'
1740 end if
1741
1742 ! loop over all header indices of the 'GL' family
1743 call obs_set_current_header_list(obsSpaceData, 'GL')
1744
1745 HEADER: do
1746
1747 headerIndex = obs_getHeaderIndex(obsSpaceData)
1748 if (headerIndex < 0) exit HEADER
1749
1750 cstnid = obs_elem_c ( obsSpaceData, 'STID' , headerIndex )
1751 codeType = obs_headElem_i( obsSpaceData, OBS_ITY, headerIndex )
1752
1753 inPlatformList = .false.
1754
1755 PLATFORM: do platformIndex = 1, nPlatformIce
1756
1757 if ( index(cstnid,trim(listPlatformIce(platformIndex))) > 0 .or. &
1758 index(codtyp_get_name(codeType),trim(listPlatformIce(platformIndex))) > 0) then
1759
1760 inPlatformList = .true.
1761 exit PLATFORM
1762
1763 end if
1764
1765 end do PLATFORM
1766
1767 if ( .not. inPlatformList ) then
1768
1769 call obs_set_current_body_list(obsSpaceData, headerIndex)
1770 BODY: do
1771 bodyIndex = obs_getBodyIndex(obsSpaceData)
1772 if (bodyIndex < 0) exit BODY
1773
1774 call obs_bodySet_i(obsSpaceData, OBS_ASS, bodyIndex, obs_notAssimilated)
1775
1776 end do BODY
1777
1778 end if
1779
1780 end do HEADER
1781
1782 if (.not. beSilent) write(*,*) 'filt_iceConcentration: end'
1783
1784 end subroutine filt_iceConcentration
1785
1786 !--------------------------------------------------------------------------
1787 ! filt_topoChemistry
1788 !--------------------------------------------------------------------------
1789 subroutine filt_topoChemistry(columnTrlOnTrlLev, obsSpaceData, beSilent)
1790 !
1791 !:Purpose: Rejects elements which are too far below the model surface
1792 ! or above the model top.
1793 !
1794 !:Comments:
1795 ! Flagging of bit 4 in OBS_FLG done in filt_topoChemistry instead of set_scale_chm
1796 ! since this subroutine is called after chm_setup, allowing use of utl_open_asciifile
1797 !
1798 implicit none
1799
1800 ! Arguments:
1801 type(struct_columnData), intent(in) :: columnTrlOnTrlLev
1802 type(struct_obs), intent(inout) :: obsSpaceData
1803 logical, intent(in) :: beSilent
1804
1805 ! Locals:
1806 integer :: headerIndex, bodyIndex, listIndex, elemIndex, listIndex_stnid
1807 integer :: ivnm, countAssim, jl, icount
1808 real(8) :: obsAltitude, obsPressure, colTopPressure, colSfcPressure
1809 real(8) :: colAltitudeBelow, colAltitudeAbove
1810 logical :: list_is_empty
1811 integer, parameter :: Nmax=100
1812 integer :: Num_stnid_chm,nobslev,Num_chm
1813 character(len=13) :: CstnidList_chm(Nmax)
1814 integer :: countAcc_stnid(Nmax),countRej_stnid(Nmax)
1815 integer :: countRejflg_stnid(Nmax),countRejflg(Nmax)
1816 integer :: countAcc(Nmax),countRej(Nmax),iConstituentList(Nmax)
1817
1818 if (.not.obs_famExist(obsSpaceData,'CH')) return
1819
1820 ! Set counters to zero
1821 countAcc_stnid(:)=0
1822 countRej_stnid(:)=0
1823 countRejflg_stnid(:)=0
1824 Num_stnid_chm=0
1825
1826 countAcc(:)=0
1827 countRej(:)=0
1828 countRejflg(:)=0
1829 Num_chm=0
1830
1831 ! Loop over all header indices of the 'CH' family
1832 call obs_set_current_header_list(obsSpaceData, 'CH')
1833 HEADER: do
1834 headerIndex = obs_getHeaderIndex(obsSpaceData)
1835 if (headerIndex < 0) exit HEADER
1836
1837 ! Set the body list & start at the beginning of the list
1838 call obs_set_current_body_list(obsSpaceData, headerIndex,list_is_empty)
1839
1840 if (list_is_empty) cycle HEADER ! Proceed to next HEADER
1841
1842 ! Set geopotential height and pressure boundaries.
1843
1844 colAltitudeBelow = col_getHeight(columnTrlOnTrlLev,0,headerIndex,'SF')
1845 colAltitudeAbove = col_getHeight(columnTrlOnTrlLev,1,headerIndex,'MM')
1846 colSfcPressure = col_getElem(columnTrlOnTrlLev,1,headerIndex,'P0')
1847 colTopPressure = col_getPressure(columnTrlOnTrlLev,1,headerIndex,'MM')
1848
1849 ! Identify max number of profile points in the profile (exclude BUFR_SCALE_EXPONENT elements)
1850
1851 nobslev = obs_headElem_i(obsSpaceData,OBS_NLV,headerIndex)
1852 bodyIndex =obs_headElem_i(obsSpaceData,OBS_RLN,headerIndex)
1853 do jl=0,obs_headElem_i(obsSpaceData,OBS_NLV,headerIndex)-1
1854 if (obs_bodyElem_i(obsSpaceData,OBS_VNM,bodyIndex+jl) == BUFR_SCALE_EXPONENT) &
1855 nobslev = nobslev-1
1856 end do
1857
1858 ! Identify element index of stnid list for the CH family.
1859
1860 call utl_get_stringId(obs_elem_c(obsSpaceData,'STID',headerIndex),&
1861 nobslev,CstnidList_chm,Num_stnid_chm,Nmax,listIndex_stnid)
1862
1863 ! Loop over all body indices (still in the 'CH' family)
1864 icount=0
1865 BODY: do
1866 bodyIndex = obs_getBodyIndex(obsSpaceData)
1867 if (bodyIndex < 0) exit BODY
1868
1869 ivnm = obs_bodyElem_i(obsSpaceData,OBS_VNM,bodyIndex)
1870 if (ivnm == BUFR_SCALE_EXPONENT) cycle BODY
1871
1872 ! Identify element index of constituent ID list for the CH family.
1873 if (icount == 0) call utl_get_Id(obs_headElem_i(obsSpaceData,OBS_CHM,headerIndex), &
1874 iConstituentList,Num_chm,Nmax,listIndex)
1875 icount=icount+1
1876
1877 ! Check for bit 4 of OBS_FLG, indicating a 'Suspicious element'
1878 if (btest(obs_bodyElem_i(obsSpaceData,OBS_FLG,bodyIndex),4)) &
1879 call obs_bodySet_i(obsSpaceData,OBS_ASS,bodyIndex,obs_notAssimilated)
1880
1881 if (obs_bodyElem_i(obsSpaceData,OBS_ASS,bodyIndex) == obs_notAssimilated) then
1882
1883 ! Already rejected from input marker/flag.
1884
1885 countRej(listIndex)=countRej(listIndex)+1
1886 countRej_stnid(listIndex_stnid)=countRej_stnid(listIndex_stnid)+1
1887 countRejflg(listIndex)=countRejflg(listIndex)+1
1888 countRejflg_stnid(listIndex_stnid)=countRejflg_stnid(listIndex_stnid)+1
1889
1890 else if (obs_bodyElem_i(obsSpaceData,OBS_VCO,bodyIndex) == 1) then
1891
1892 ! Check as a function of altitude.
1893
1894 obsAltitude = obs_bodyElem_r(obsSpaceData,OBS_PPP,bodyIndex)
1895 !write(*,*) 'rejected zzz ',obs_elem_c(obsSpaceData,'STID',headerIndex),obsAltitude,colSfcPressure,colTopPressure
1896 if( obsAltitude < colAltitudeBelow .or. obsAltitude > colAltitudeAbove ) then
1897 call obs_bodySet_i(obsSpaceData,OBS_FLG,bodyIndex, &
1898 ibset( obs_bodyElem_i(obsSpaceData,OBS_FLG,bodyIndex), 18 ))
1899 call obs_bodySet_i(obsSpaceData,OBS_ASS,bodyIndex,obs_notAssimilated)
1900 countRej(listIndex)=countRej(listIndex)+1
1901 countRej_stnid(listIndex_stnid)=countRej_stnid(listIndex_stnid)+1
1902 else
1903 countAcc(listIndex)=countAcc(listIndex)+1
1904 countAcc_stnid(listIndex_stnid)=countAcc_stnid(listIndex_stnid)+1
1905 end if
1906
1907 else if (obs_bodyElem_i(obsSpaceData,OBS_VCO,bodyIndex) == 2) then
1908
1909 ! Check as a function of pressure.
1910
1911 obsPressure = obs_bodyElem_r(obsSpaceData,OBS_PPP,bodyIndex)
1912
1913 if ( obsPressure > colSfcPressure .or. obsPressure < colTopPressure) then
1914 call obs_bodySet_i(obsSpaceData,OBS_ASS,bodyIndex,obs_notAssimilated)
1915 call obs_bodySet_i(obsSpaceData,OBS_FLG,bodyIndex, &
1916 ibset(obs_bodyElem_i(obsSpaceData,OBS_FLG,bodyIndex),18 ))
1917 countRej(listIndex)=countRej(listIndex)+1
1918 countRej_stnid(listIndex_stnid)=countRej_stnid(listIndex_stnid)+1
1919 else
1920 countAcc(listIndex)=countAcc(listIndex)+1
1921 countAcc_stnid(listIndex_stnid)=countAcc_stnid(listIndex_stnid)+1
1922 end if
1923 else
1924 countAcc(listIndex)=countAcc(listIndex)+1
1925 countAcc_stnid(listIndex_stnid)=countAcc_stnid(listIndex_stnid)+1
1926 end if
1927
1928 end do BODY
1929
1930 end do HEADER
1931
1932 if (Num_stnid_chm > 0 .and. .not.beSilent) then
1933 write(*,*) ' '
1934 write(*,*) '*****************************************************************'
1935 write(*,*) ' filt_topoChemistry: '
1936 write(*,*) ' FAMILY = CH'
1937 write(*,222) 'ELEMENTS for CH stnids',(CstnidList_chm(elemIndex),elemIndex=1,Num_stnid_chm)
1938 write(*,223) 'ACCEPTED for CH stnids',(countAcc_stnid(elemIndex),elemIndex=1,Num_stnid_chm)
1939 write(*,223) 'REJECTED for CH stnids',(countRej_stnid(elemIndex),elemIndex=1,Num_stnid_chm)
1940 write(*,223) 'REJECTED due to marker',(countRejflg_stnid(elemIndex),elemIndex=1,Num_stnid_chm)
1941 write(*,*) ' '
1942 write(*,224) 'ELEMENTS for CH ',(iConstituentList(elemIndex),elemIndex=1,Num_chm)
1943 write(*,224) 'ACCEPTED for CH ',(countAcc(elemIndex),elemIndex=1,Num_chm)
1944 write(*,224) 'REJECTED for CH ',(countRej(elemIndex),elemIndex=1,Num_chm)
1945 write(*,223) 'REJECTED due to marker',(countRejflg(elemIndex),elemIndex=1,Num_stnid_chm)
1946 write(*,*) '*****************************************************************'
1947 write(*,*) ' '
1948
1949 countAssim=0
1950 do bodyIndex=1,obs_numbody(obsSpaceData)
1951 if (obs_bodyElem_i(obsSpaceData,OBS_ASS,bodyIndex) == obs_assimilated) countAssim=countAssim+1
1952 end do
1953 write(*,'(1X," NUMBER OF DATA TO BE ASSIMILATED AFTER ADJUSTMENTS (after filter_topoChemistry):",i10)') countAssim
1954 write(*,*) ' '
1955 end if
1956222 format(2x,a29,100(2x,a10))
1957223 format(2x,a29,100(2x,i8,2x))
1958224 format(2x,a29,100(2x,i6))
1959
1960 end subroutine filt_topoChemistry
1961
1962 !--------------------------------------------------------------------------
1963 ! filt_bufrCodeAssimilated
1964 !-------------------------------------------------------------------------
1965 function filt_bufrCodeAssimilated(bufrCode) result(assimilated)
1966 !
1967 !:Purpose: To test if a bufr code part of the assimilated observation list
1968 !
1969 implicit none
1970
1971 ! Arguments:
1972 integer, intent(in) :: bufrCode ! The input bufr code
1973 ! Result:
1974 logical :: assimilated ! Assimilated of not
1975
1976 ! Locals:
1977 integer :: elemIndex
1978
1979 if (.not. initialized) call filt_setup('none')
1980
1981 assimilated = .false.
1982
1983 do elemIndex = 1, filt_nelems
1984 if (filt_nlist(elemIndex) == bufrCode) then
1985 assimilated = .true.
1986 return
1987 end if
1988 end do
1989
1990 end function filt_bufrCodeAssimilated
1991
1992 !--------------------------------------------------------------------------
1993 ! filt_getBufrCodeAssimilated
1994 !-------------------------------------------------------------------------
1995 subroutine filt_getBufrCodeAssimilated(bufrCodeList)
1996 !
1997 !:Purpose: To get the assimilated observation list
1998 !
1999 implicit none
2000
2001 ! Arguments:
2002 integer, intent(out) :: bufrCodeList(filt_nelems) ! The list of assimilated bufr codes
2003
2004 if (.not. initialized) call filt_setup('none')
2005
2006 bufrCodeList(:) = filt_nlist(1:filt_nelems)
2007
2008 end subroutine filt_getBufrCodeAssimilated
2009
2010 !--------------------------------------------------------------------------
2011 ! filt_nBufrCodeAssimilated
2012 !-------------------------------------------------------------------------
2013 function filt_nBufrCodeAssimilated() result(nBufrCode)
2014 !
2015 !:Purpose: To get the number of assimilated observations
2016 !
2017 implicit none
2018
2019 ! Result:
2020 integer :: nBufrCode ! The number of assimilated observations
2021
2022 if (.not. initialized) call filt_setup('none')
2023
2024 nBufrCode = filt_nelems
2025
2026 end function filt_nBufrCodeAssimilated
2027
2028end module obsFilter_mod
