1module bgckOcean_mod
2 ! MODULE bgckOcean_mod (prefix='ocebg' category='1. High-level functionality')
3 !
4 !:Purpose: to perform ocean data background check.
5 !
6 use midasMpi_mod
7 use utilities_mod
8 use obsSpaceData_mod
9 use columnData_mod
10 use codtyp_mod
11 use gridStateVector_mod
12 use gridStateVectorFileIO_mod
13 use horizontalCoord_mod
14 use verticalCoord_mod
15 use statetocolumn_mod
16 use bufr_mod
17 use mathPhysConstants_mod
18 use timeCoord_mod
19 use message_mod
20
21 implicit none
22
23 save
24 private
25
26 ! Public functions/subroutines
27 public :: ocebg_bgCheckSST, ocebg_bgCheckSeaIce
28
29 ! External functions
30 integer, external :: fnom, fclos
31
32 ! mpi topology
33 integer :: myLatBeg, myLatEnd
34 integer :: myLonBeg, myLonEnd
35 integer :: latPerPE, latPerPEmax, lonPerPE, lonPerPEmax
36
37 ! namelist variables with setting default values
38 character(len=20) :: timeInterpType_nl = 'NEAREST' ! 'NEAREST' or 'LINEAR'
39 integer :: numObsBatches = 20 ! number of batches for calling interp setup
40 logical :: checkWinds = .false. ! check winds for last 4 days to amplify error in zone of max wind speed
41 integer :: ndaysWinds = 4 ! number of days in the 'winds' file to detect tropical storm (TS)
42 integer :: timeStepWinds = 6 ! in hours, winds are available every timeStepWinds-hours
43 integer :: windForecastLeadtime = 6 ! in hours, lead time of wind forecast in the input file
44 real(4) :: minLatNH = 10. ! min lat of N. hemisphere lat band where TS is detected
45 real(4) :: maxLatNH = 40. ! max lat of N. hemisphere lat band where TS is detected
46 real(4) :: maxLatExceptionNH = 45. ! max lat of N. hemisphere lat band allows TS to penetrate further North in some months
47 integer :: nmonthsExceptionNH = 0 ! MUST NOT BE INCLUDED IN NAMELIST!
48 character(len=3) :: monthExceptionNH(12) = ' ' ! exceptional months where TS allowed to penetrated further North
49 real(4) :: minLatSH = -35. ! min lat of Southern hemisphere latutude band where TS is detected
50 real(4) :: maxLatSH = -10. ! max lat of Southern hemisphere latutude band where TS is detected
51 real(8) :: smoothLenghtScale = 50000. ! length scale. in m, to smooth the amplification error field
52 real(8) :: globalSelectCriteria(3) = (/5.d0, 25.d0, 30.d0/) ! global selection criteria
53 logical :: separateSelectCriteria = .false. ! apply separate selection criteria: sea/inland waters; insitu/satellite
54 real(8) :: inlandWaterSelectCriteriaSatData(3) = (/5.d0, 25.d0, 30.d0/) ! inland water, satellite selection criteria
55 real(8) :: inlandWaterSelectCriteriaInsitu(3) = (/5.d0, 25.d0, 30.d0/) ! inland water, insitu selection criteria
56 real(8) :: seaWaterSelectCriteriaSatData(3) = (/5.d0, 25.d0, 30.d0/) ! sea water, satellite selection criteria
57 real(8) :: seaWaterSelectCriteriaInsitu(3) = (/5.d0, 25.d0, 30.d0/) ! sea water, insitu selection criteria
58 real(4) :: seaWaterThreshold = 0.1 ! threshold to distinguish inland water from sea water
59 namelist /namOceanBGcheck/ timeInterpType_nl, numObsBatches, checkWinds, ndaysWinds, timeStepWinds, &
60 windForecastLeadtime, minLatNH, maxLatNH, maxLatExceptionNH, nmonthsExceptionNH, &
61 monthExceptionNH, minLatSH, maxLatSH, smoothLenghtScale, &
62 globalSelectCriteria, separateSelectCriteria, inlandWaterSelectCriteriaSatData, &
63 inlandWaterSelectCriteriaInsitu, seaWaterSelectCriteriaSatData, &
64 seaWaterSelectCriteriaInsitu, seaWaterThreshold
65
66 character(len=3), parameter :: months(12) = ['Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec']
67
68 contains
69
70 !----------------------------------------------------------------------------------------
71 ! ocebg_bgCheckSST
72 !----------------------------------------------------------------------------------------
73 subroutine ocebg_bgCheckSST(obsData, dateStamp, columnTrlOnTrlLev, hco)
74 !
75 !:Purpose: to compute SST data background Check
76 !
77
78 implicit none
79
80 ! Arguments:
81 type(struct_obs) , intent(inout) :: obsData ! obsSpaceData object
82 integer , intent(in) :: dateStamp ! date stamp
83 type(struct_columnData), intent(inout) :: columnTrlOnTrlLev ! column data on trl levels
84 type(struct_hco), pointer, intent(in) :: hco ! horizontal trl grid
85
86 ! Locals:
87 type(struct_gsv) :: stateVectorFGE ! state vector containing std B estimation field
88 type(struct_gsv) :: stateVectorAmplFactor ! state vector for error amplification field
89 real(4), pointer :: stateVectorAmplFactor_ptr(:,:,:)
90 type(struct_gsv) :: stateVectorSeaWaterFraction ! statevector for sea water fraction
91 integer :: nulnam, ierr, headerIndex, bodyIndex, obsFlag, obsVarno
92 integer :: numberObs, numberObsRejected
93 integer :: numberObsInsitu, numberObsInsituRejected, codeType
94 real(8) :: OER, OmP, FGE, bgCheck, seaWaterFraction
95 type(struct_columnData) :: columnFGE, columnSeaWaterFraction
96 real(4), pointer :: stateVectorFGE_ptr(:,:,:)
97 logical :: checkMonth, llok
98 integer :: lonIndex, latIndex, monthIndex, exceptMonthIndex
99
100 call msg('ocebg_bgCheckSST', 'performing background check for the SST data...')
101
102 ! get mpi topology
103 call mmpi_setup_lonbands(hco%ni, lonPerPE, lonPerPEmax, myLonBeg, myLonEnd)
104 call mmpi_setup_latbands(hco%nj, latPerPE, latPerPEmax, myLatBeg, myLatEnd)
105
106 ! Read the namelist
107 if (.not. utl_isNamelistPresent('namOceanBGcheck','./flnml')) then
108 if (mmpi_myid == 0) then
109 call msg('ocebg_bgCheckSST', 'namOceanBGcheck is missing in the namelist.')
110 call msg('ocebg_bgCheckSST', 'the default values will be taken.')
111 end if
112 else
113 ! reading namelist variables
114 nulnam = 0
115 ierr = fnom(nulnam, './flnml', 'FTN+SEQ+R/O', 0)
116 read(nulnam, nml = namOceanBGcheck, iostat = ierr)
117 if (ierr /= 0) call utl_abort('ocebg_bgCheckSST: Error reading namelist')
118 ierr = fclos(nulnam)
119 if (nmonthsExceptionNH /=0) then
120 call utl_abort('ocebg_bgCheckSST: check namOceanBGcheck namelist section: nmonthsExceptionNH should be removed')
121 end if
122 do monthIndex = 1, 12
123 if (monthExceptionNH(monthIndex) == ' ') exit
124 nmonthsExceptionNH = nmonthsExceptionNH + 1
125 end do
126 end if
127 call msg('ocebg_bgCheckSST', 'interpolation type: '//timeInterpType_nl)
128 call msg('ocebg_bgCheckSST', 'number obs batches: '//str(numObsBatches))
129 call msg('ocebg_bgCheckSST', 'check winds to detect tropical storms (TS): '//str(checkWinds))
130
131 if (separateSelectCriteria) then
132 call msg('ocebg_bgCheckSST', 'a separate selection criteria will be applied...')
133 write(*,'(a,3f7.2)') 'ocebg_bgCheckSST: inland water selection criteria for INSITU data rejection set to: ', &
134 inlandWaterSelectCriteriaInsitu(:)
135 write(*,'(a,3f7.2)') 'ocebg_bgCheckSST: inland water selection criteria for SATELLITE data rejection set to: ', &
136 inlandWaterSelectCriteriaSatData(:)
137 write(*,'(a,3f7.2)') 'ocebg_bgCheckSST: sea water selection criteria for INSITU data rejection set to: ', &
138 seaWaterSelectCriteriaInsitu(:)
139 write(*,'(a,3f7.2)') 'ocebg_bgCheckSST: sea water selection criteria for SATELLITE data rejection set to: ', &
140 seaWaterSelectCriteriaSatData(:)
141 call msg('ocebg_bgCheckSST', 'sea water fraction threshold is set to: '//str(seaWaterThreshold))
142 ! read sea water fraction
143 call gsv_allocate(stateVectorSeaWaterFraction, 1, hco, columnTrlOnTrlLev%vco, dataKind_opt = 4, &
144 datestamp_opt = -1, mpi_local_opt = .true., &
145 varNames_opt = (/'VF'/), hInterpolateDegree_opt ='LINEAR')
146 call gio_readFromFile(stateVectorSeaWaterFraction, './seaice_analysis', ' ','A', &
147 unitConversion_opt=.false., containsFullField_opt=.true.)
148 ! Convert sea water fraction stateVector to column object
149 call col_setVco(columnSeaWaterFraction, col_getVco(columnTrlOnTrlLev))
150 call col_allocate(columnSeaWaterFraction, col_getNumCol(columnTrlOnTrlLev), varNames_opt = (/'VF'/))
151 call s2c_nl(stateVectorSeaWaterFraction, obsData, columnSeaWaterFraction, hco, varName_opt = 'VF', &
152 timeInterpType = timeInterpType_nl, moveObsAtPole_opt = .true., &
153 numObsBatches_opt = numObsBatches, dealloc_opt = .true.)
154 else
155 write(*,'(a,3f7.2)') 'ocebg_bgCheckSST: global selection criteria for data rejection is set to: ', globalSelectCriteria(:)
156 end if
157
158 ! Read First Guess Error (FGE) and put it into stateVector
159 call gsv_allocate(stateVectorFGE, 1, hco, columnTrlOnTrlLev%vco, dataKind_opt = 4, &
160 hInterpolateDegree_opt = 'NEAREST', &
161 datestamp_opt = -1, mpi_local_opt = .true., varNames_opt = (/'TM'/))
162 call gio_readFromFile(stateVectorFGE, './bgstddev', 'STDDEV', 'X', &
163 unitConversion_opt=.false., containsFullField_opt=.true.)
164 if (checkWinds) then
165 call utl_tmg_start(122, '--checkWindsForSST')
166 call msg('ocebg_bgCheckSST', 'looking for tropical storms...')
167 call msg('ocebg_bgCheckSST', 'number of days with available winds in the input winds file: '//str(ndaysWinds))
168 call msg('ocebg_bgCheckSST', 'winds are provided every: '//str(timeStepWinds)//' hours')
169 call msg('ocebg_bgCheckSST', 'wind forecast lead time: '//str(windForecastLeadtime)//' hours')
170 do exceptMonthIndex = 1, nmonthsExceptionNH
171 checkMonth = .False.
172 loop_month: do monthIndex = 1, 12
173 if (monthExceptionNH(exceptMonthIndex) == months(monthIndex)) then
174 checkMonth = .True.
175 exit loop_month
176 end if
177 end do loop_month
178 if (.not. checkMonth) then
179 write(*,*) 'ocebg_bgCheckSST: month should be one of these: ',months(:)
180 call utl_abort('ocebg_bgCheckSST: unknown month '//monthExceptionNH(exceptMonthIndex))
181 end if
182 end do
183
184 ! amplification error field state vector
185 call gsv_allocate(stateVectorAmplFactor, 1, hco, columnTrlOnTrlLev%vco, dataKind_opt = 4, &
186 hInterpolateDegree_opt = 'LINEAR', datestamp_opt = dateStamp, &
187 mpi_local_opt = .true., varNames_opt = (/'TM'/))
188 call gsv_getField(stateVectorAmplFactor, stateVectorAmplFactor_ptr)
189 stateVectorAmplFactor_ptr(myLonBeg:myLonEnd,myLatBeg:myLatEnd,1) = 1.0d0
190
191 call ocebg_getFGEamplification(stateVectorAmplFactor, dateStamp, hco)
192
193 ! FGE state vector
194 call gsv_getField(stateVectorFGE, stateVectorFGE_ptr)
195
196 ! Apply tropical storm correction to the FGE field:
197 do latIndex = myLatBeg, myLatEnd
198 do lonIndex = myLonBeg, myLonEnd
199 stateVectorFGE_ptr(lonIndex, latIndex, 1) = stateVectorFGE_ptr(lonIndex, latIndex, 1) * &
200 stateVectorAmplFactor_ptr(lonIndex, latIndex, 1)
201 end do
202 end do
203 call gsv_deallocate(stateVectorAmplFactor)
204 call utl_tmg_stop(122)
205 end if
206
207 ! Convert FGE stateVector to column object
208 call col_setVco(columnFGE, col_getVco(columnTrlOnTrlLev))
209 call col_allocate(columnFGE, col_getNumCol(columnTrlOnTrlLev), varNames_opt = (/'TM'/))
210 call s2c_nl(stateVectorFGE, obsData, columnFGE, hco, varName_opt = 'TM', &
211 timeInterpType = timeInterpType_nl, moveObsAtPole_opt = .true., &
212 numObsBatches_opt = numObsBatches, dealloc_opt = .true.)
213
214 numberObs = 0
215 numberObsRejected = 0
216 numberObsInsitu = 0
217 numberObsInsituRejected = 0
218
219 do headerIndex = 1, obs_numheader(obsData)
220
221 bodyIndex = obs_headElem_i(obsData, obs_rln, headerIndex)
222 obsVarno = obs_bodyElem_i(obsData, obs_vnm, bodyIndex)
223 llok = (obs_bodyElem_i(obsData, obs_ass, bodyIndex) == obs_assimilated)
224 if (llok) then
225 if (obsVarno == bufr_sst) then
226
227 FGE = col_getElem(columnFGE, 1, headerIndex, 'TM')
228 OmP = obs_bodyElem_r(obsData, OBS_OMP, bodyIndex)
229 OER = obs_bodyElem_r(obsData, OBS_OER, bodyIndex)
230 codeType = obs_headElem_i(obsData, obs_ity, headerIndex)
231
232 if (FGE /= MPC_missingValue_R8 .and. OmP /= MPC_missingValue_R8) then
233
234 numberObs = numberObs + 1
235 if (codeType /= codtyp_get_codtyp('satob')) numberObsInsitu = numberObsInsitu + 1
236 call obs_bodySet_r(obsData, OBS_HPHT, bodyIndex, FGE)
237 bgCheck = (OmP)**2 / (FGE**2 + OER**2)
238
239 if (separateSelectCriteria) then
240 seaWaterFraction = col_getElem(columnSeaWaterFraction, 1, headerIndex, 'VF')
241 if (seaWaterFraction <= seaWaterThreshold) then
242 if(codeType == codtyp_get_codtyp('satob')) then
243 obsFlag = ocebg_setFlag(obsVarno, bgCheck, inlandWaterSelectCriteriaSatData)
244 else
245 obsFlag = ocebg_setFlag(obsVarno, bgCheck, inlandWaterSelectCriteriaInsitu)
246 end if
247 else
248 if(codeType == codtyp_get_codtyp('satob')) then
249 obsFlag = ocebg_setFlag(obsVarno, bgCheck, seaWaterSelectCriteriaSatData)
250 else
251 obsFlag = ocebg_setFlag(obsVarno, bgCheck, seaWaterSelectCriteriaInsitu)
252 end if
253 end if
254 else
255 obsFlag = ocebg_setFlag(obsVarno, bgCheck, globalSelectCriteria)
256 end if
257
258 if (obsFlag >= 2) then
259 numberObsRejected = numberObsRejected + 1
260 if (codeType /= codtyp_get_codtyp('satob')) numberObsInsituRejected = numberObsInsituRejected + 1
261 write(*,'(i10,a,i5,4f10.4,i5)') numberObsRejected, &
262 ', sensor, codtype, lon, lat, obs.value, OmP, flag: '&
263 //obs_elem_c(obsData, 'STID' , headerIndex)//' data: ', codeType, &
264 obs_headElem_r(obsData, obs_lon, headerIndex) * MPC_DEGREES_PER_RADIAN_R8,&
265 obs_headElem_r(obsData, obs_lat, headerIndex) * MPC_DEGREES_PER_RADIAN_R8,&
266 obs_bodyElem_r(obsData, obs_var, bodyIndex), OmP, obsFlag
267 end if
268
269 ! update background check flags based on bgCheck
270 ! (element flags + global header flags)
271 if (obsFlag == 1) then
272 call obs_bodySet_i(obsData, obs_flg, bodyIndex , ibset(obs_bodyElem_i(obsData, obs_flg, bodyIndex) , 13))
273 else if (obsFlag == 2) then
274 call obs_bodySet_i(obsData, obs_flg, bodyIndex , ibset(obs_bodyElem_i(obsData, obs_flg, bodyIndex) , 14))
275 call obs_bodySet_i(obsData, obs_flg, bodyIndex , ibset(obs_bodyElem_i(obsData, obs_flg, bodyIndex) , 16))
276 call obs_bodySet_i(obsData, obs_flg, bodyIndex , ibset(obs_bodyElem_i(obsData, obs_flg, bodyIndex) , 09))
277 call obs_headSet_i(obsData, obs_st1, headerIndex, ibset(obs_headElem_i(obsData, obs_st1, headerIndex), 06))
278 else if (obsFlag == 3) then
279 call obs_bodySet_i(obsData, obs_flg, bodyIndex , ibset(obs_bodyElem_i(obsData, obs_flg, bodyIndex) , 15))
280 call obs_bodySet_i(obsData, obs_flg, bodyIndex , ibset(obs_bodyElem_i(obsData, obs_flg, bodyIndex) , 16))
281 call obs_bodySet_i(obsData, obs_flg, bodyIndex , ibset(obs_bodyElem_i(obsData, obs_flg, bodyIndex) , 09))
282 call obs_headSet_i(obsData, obs_st1, headerIndex, ibset(obs_headElem_i(obsData, obs_st1, headerIndex), 06))
283 end if
284
285 end if
286 end if
287 end if
288
289 end do
290
291 if (numberObs > 0) then
292 write(*,*)
293 call msg('ocebg_bgCheckSST', 'background check of SST (TM) data is computed')
294 write(*,*) '***************************************************************************************'
295 write(*,'(a, i7,a,i7,a)') 'ocebg_bgCheckSST: total ', numberObsRejected, ' observations out of (ALL) ', numberObs,' rejected'
296 write(*,'(a, i7,a,i7,a)') 'ocebg_bgCheckSST: where ', numberObsInsituRejected, ' insitu observations out of ', &
297 numberObsInsitu,' insitu obs. rejected'
298 write(*,*) '***************************************************************************************'
299 write(*,*)
300 end if
301
302 call gsv_deallocate(stateVectorFGE)
303 call col_deallocate(columnFGE)
304 if (separateSelectCriteria) then
305 call col_deallocate(columnSeaWaterFraction)
306 call gsv_deallocate(stateVectorSeaWaterFraction)
307 end if
308
309 end subroutine ocebg_bgCheckSST
310
311 !----------------------------------------------------------------------------------------
312 ! ocebg_bgCheckSeaIce
313 !----------------------------------------------------------------------------------------
314 subroutine ocebg_bgCheckSeaIce(obsData)
315 !
316 !:Purpose: Compute sea ice data background check.
317 ! The rms of the difference between the observations and
318 ! the background values is calculated over a "swath" of data.
319 ! For satellite observations, this is based on a single passage
320 ! of the satellite over the domain.
321 ! A flag is then set for all observations in the swath
322 ! if the rms value is larger than the threshold
323 ! specified for the observation type.
324 !
325
326 implicit none
327
328 ! Arguments:
329 type(struct_obs), intent(inout) :: obsData ! obsSpaceData object
330
331 ! Locals:
332 integer :: nulnam, ierr
333 integer :: headerIndex, bodyIndex, stationIndex, bodyCount
334 integer :: obsChid, obsDate, obsTime, obsFlag
335 integer :: obsDateStamp
336 integer, allocatable :: numberObs(:), bodyIndexList(:,:)
337 integer, allocatable :: minDateStamp(:), maxDateStamp(:), swathID(:)
338 real, allocatable :: rmsDiff(:)
339 integer :: swathIndex, nSwath
340 real(8) :: OmP, deltaHours
341 integer :: numberObsRejected
342 character(len=12) :: cstnid
343 integer, external :: newdate
344 integer :: imode, prntdate, prnttime
345 integer, parameter :: numStationMax = 14 ! maximum number of 'idStation' values
346
347 ! Namelist variables: (local)
348 integer :: numStation = MPC_missingValue_INT ! MUST NOT BE INCLUDED IN NAMELIST!
349 character(len=12) :: idStation(numStationMax) = 'null' ! list of obsSpaceData 'idStation' values to consider
350 real :: OmpRmsdThresh(numStationMax) = 0.0 ! rejection threshold applied to RMS of O-P for entire swath
351 integer :: maxSwath = 10 ! maximum number of swaths
352 integer :: maxPerSwath = 200000 ! maximum number of data per swath
353 namelist /namIceBGcheck/ numStation, idStation, OmpRmsdThresh, maxSwath, maxPerSwath
354
355 call msg('ocebg_bgCheckSeaIce', 'performing background check for the SeaIce data...')
356
357 ! Read the namelist
358 if (.not. utl_isNamelistPresent('namIceBGcheck','./flnml')) then
359 if (mmpi_myid == 0) then
360 call msg('ocebg_bgCheckSeaIce', 'namIceBGcheck is missing in the namelist.')
361 call msg('ocebg_bgCheckSeaIce', 'the default values will be taken.')
362 end if
363 else
364 ! reading namelist variables
365 nulnam = 0
366 ierr = fnom(nulnam, './flnml', 'FTN+SEQ+R/O', 0)
367 read(nulnam, nml = namIceBGcheck, iostat = ierr)
368 if (ierr /= 0) call utl_abort('ocebg_bgCheckSeaIce: Error reading namelist')
369 ierr = fclos(nulnam)
370 if (numStation /= MPC_missingValue_INT) then
371 call utl_abort('ocebg_bgCheckSeaIce: check namIceBGcheck namelist section: numStation should be removed')
372 end if
373 numStation = 0
374 do stationIndex = 1, numStationMax
375 if (trim(idStation(stationIndex)) == 'null') exit
376 numStation = numStation + 1
377 end do
378 end if
379 write(*, nml = namIceBGcheck)
380
381 allocate(numberObs(maxSwath), bodyIndexList(maxPerSwath,maxSwath))
382 allocate(minDateStamp(maxSwath), maxDateStamp(maxSwath), swathID(maxSwath))
383 allocate(rmsDiff(maxSwath))
384
385 STATION: do stationIndex = 1, numStation
386
387 call msg('ocebg_bgCheckSeaIce', 'doing idStation: '//idStation(stationIndex))
388
389 numberObs(:) = 0
390 rmsDiff(:) = 0.0
391 minDateStamp(:) = tim_getDatestamp()
392 maxDateStamp(:) = tim_getDatestamp()
393 nSwath = 0
394
395 HEADER: do headerIndex = 1, obs_numheader(obsData)
396
397 cstnid = obs_elem_c(obsData, 'STID' , headerIndex )
398 bodyIndex = obs_headElem_i(obsData, obs_rln, headerIndex)
399
400 obsFlag = obs_bodyElem_i(obsData, OBS_FLG, bodyIndex)
401
402 if ( (cstnid /= idStation(stationIndex)) .or. &
403 (obsFlag /= 0 .and. cstnid /= 'CIS_REGIONAL') ) cycle HEADER
404
405 obsChid = obs_headElem_i(obsData, obs_chid, headerIndex)
406 obsDate = obs_headElem_i(obsData, OBS_DAT, headerIndex)
407 obsTime = obs_headElem_i(obsData, OBS_ETM, headerIndex)
408
409 imode = 3 ! printable to stamp
410 ierr = newdate(obsDateStamp, obsDate, obsTime*10000, imode)
411
412 do swathIndex=1, nSwath
413
414 if(obsDateStamp <= maxDateStamp(swathIndex) .and. &
415 obsDateStamp >= minDateStamp(swathIndex) .and. &
416 obsChid == swathID(swathIndex)) exit
417
418 call difdatr(obsDateStamp, minDateStamp(swathIndex), deltaHours)
419
420 if(abs(deltaHours) <= 0.5d0 .and. obsChid == swathID(swathIndex)) exit
421
422 end do
423
424 nSwath = max(nSwath, swathIndex)
425
426 if(swathIndex <= maxSwath) then
427
428 if(numberObs(swathIndex) == 0) then
429 minDateStamp(swathIndex) = obsDateStamp
430 maxDateStamp(swathIndex) = obsDateStamp
431 else
432 call difdatr(obsDateStamp, minDateStamp(swathIndex), deltaHours)
433 if(deltaHours < 0.0d0) minDateStamp(swathIndex) = obsDateStamp
434 call difdatr(obsDateStamp, maxDateStamp(swathIndex), deltaHours)
435 if(deltaHours > 0.0d0) maxDateStamp(swathIndex) = obsDateStamp
436 end if
437
438 swathID(swathIndex) = obsChid
439
440 OmP = obs_bodyElem_r(obsData, OBS_OMP, bodyIndex)
441
442 if (OmP /= MPC_missingValue_R8) then
443
444 numberObs(swathIndex) = numberObs(swathIndex) + 1
445 rmsDiff(swathIndex) = rmsDiff(swathIndex) + (OmP)**2
446 if (numberObs(swathIndex) <= maxPerSwath) then
447 bodyIndexList(numberObs(swathIndex), swathIndex) = bodyIndex
448 else
449 call utl_abort('ocebg_bgCheckSeaIce: Increase maxPerSwath')
450 end if
451
452 end if
453
454 else
455
456 call utl_abort('ocebg_bgCheckSeaIce: Increase maxSwath')
457
458 end if
459
460 end do HEADER
461
462 SWATH: do swathIndex = 1, maxSwath
463
464 if (numberObs(swathIndex) == 0) cycle SWATH
465
466 rmsDiff(swathIndex) = sqrt(rmsDiff(swathIndex)/numberObs(swathIndex))
467
468 call msg('ocebg_bgCheckSeaIce', 'swathIndex = '//str(swathIndex))
469 call msg('ocebg_bgCheckSeaIce', 'Timespan:')
470 imode = -3 ! stamp to printable
471 ierr = newdate(minDateStamp(swathIndex), prntdate, prnttime, imode)
472 call msg('ocebg_bgCheckSeaIce', 'From: '//str(prntdate)//str(prnttime/100))
473 ierr = newdate(maxDateStamp(swathIndex), prntdate, prnttime, imode)
474 call msg('ocebg_bgCheckSeaIce', 'To : '//str(prntdate)//str(prnttime/100))
475
476 write(*,'(a,f10.4)') 'ocebg_bgCheckSeaIce: RMSD = ', rmsDiff(swathIndex)
477 write(*,'(a,i10)') 'ocebg_bgCheckSeaIce: nobs = ', numberObs(swathIndex)
478
479 if (rmsDiff(swathIndex) > OmpRmsdThresh(stationIndex)) then
480
481 numberObsRejected = 0
482 do bodyCount = 1, numberObs(swathIndex)
483 numberObsRejected = numberObsRejected + 1
484 ! update background check flag
485 call obs_bodySet_i(obsData, obs_flg, bodyIndexList(bodyCount,swathIndex), ibset(obs_bodyElem_i(obsData, obs_flg, bodyIndexList(bodyCount,swathIndex)), 10))
486 end do
487
488 write(*,'(a,i7,a)')'ocebg_bgCheckSeaIce: ********** reject: ', numberObsRejected, &
489 ' '//idStation(stationIndex)//' data'
490
491 end if
492
493 end do SWATH
494
495 end do STATION
496
497 deallocate(numberObs, bodyIndexList)
498 deallocate(minDateStamp, maxDateStamp, swathID)
499 deallocate(rmsDiff)
500
501 end subroutine ocebg_bgCheckSeaIce
502
503 !--------------------------------------------------------------------------
504 ! ocebg_setFlag
505 !--------------------------------------------------------------------------
506 function ocebg_setFlag(obsVarno, bgCheck, selectCriteria) result(obsFlag)
507 !
508 !:Purpose: Set background-check flags according to values set in a table.
509 ! Original values in table come from ECMWF.
510 !
511 implicit none
512
513 ! Arguments:
514 integer, intent(in) :: obsVarno ! obsVarno, Universal Field-Identity Numbers defined in bufr_mod
515 real(8), intent(in) :: bgCheck ! normalized background departure
516 real(8), intent(in) :: selectCriteria(:) ! selection criteria for three levels
517 ! Result:
518 integer :: obsFlag ! obs flag
519
520 obsFlag = 0
521
522 if (obsVarno == bufr_sst) then
523 if (bgCheck >= selectCriteria(1) .and. bgCheck < selectCriteria(2)) then
524 obsFlag = 1
525 else if (bgCheck >= selectCriteria(2) .and. bgCheck < selectCriteria(3)) then
526 obsFlag = 2
527 else if (bgCheck >= selectCriteria(3)) then
528 obsFlag = 3
529 end if
530 end if
531
532 end function ocebg_setFlag
533
534 !--------------------------------------------------------------------------
535 ! ocebg_getFGEamplification
536 !--------------------------------------------------------------------------
537 subroutine ocebg_getFGEamplification(stateVectorAmplFactor, dateStamp, hco)
538 !
539 !:Purpose: Read wind speed fields for the last four days.
540 ! In the operations:
541 ! The background error used during the background check is then
542 ! amplified in those regions by a factor that varies from 1,
543 ! where the maximum wind speed is 21m/s or less, to 12,
544 ! where the maximum wind speed is 24m/s or more.
545 ! The factor is then filtered to produce a smoothly varying
546 ! field. This amplified background error is used only to
547 ! perform the background check.
548 !
549 implicit none
550
551 ! Arguments:
552 type(struct_gsv), intent(inout) :: stateVectorAmplFactor ! state vector to save amplification factor
553 integer , intent(in) :: dateStamp ! date stamp
554 type(struct_hco), pointer, intent(in) :: hco ! horizontal trl grid
555
556 ! Locals:
557 type(struct_gsv) :: stateVector ! state vector for surface winds
558 type(struct_vco), pointer :: vco_winds ! vertical grid structure for winds
559 real(4) , pointer :: uu_ptr4d(:,:,:,:) ! pointer to get UU wind component
560 real(4) , pointer :: vv_ptr4d(:,:,:,:) ! pointer to get VV wind component
561 integer :: dataStampList(ndaysWinds * 24 / timeStepWinds) ! datastamp list for wind fields
562 real(4) :: windSpeed
563 integer :: hour, day, monthNumber
564 integer :: yyyy, ndays, timeStepIndex
565 real(8) :: deltaT, lat
566 integer :: lonIndex, latIndex, monthIndex
567 real(4) , pointer :: stateVectorAmplFactor_ptr(:,:,:)
568 real(8) :: amplFactor
569
570 nullify(vco_winds)
571
572 ! looking for the earliest valid time:
573 deltaT = - real(ndaysWinds * 24 - windForecastLeadtime)
574 call incdatr(dataStampList(1), dateStamp, deltaT)
575
576 ! scanning the datastamps up from the earliest to the last valid time
577 deltaT = real(timeStepWinds)
578 do timeStepIndex = 2, ndaysWinds * 24 / timeStepWinds
579 call incdatr(dataStampList(timeStepIndex), dataStampList(timeStepIndex - 1), deltaT)
580 end do
581
582 call vco_SetupFromFile(vco_winds,'./winds')
583 call gsv_allocate(stateVector, ndaysWinds * 24 / timeStepWinds, hco, vco_winds, &
584 dateStampList_opt = dataStampList, dataKind_opt = 4, &
585 varNames_opt=(/'UU','VV'/), mpi_local_opt=.true., &
586 hInterpolateDegree_opt='LINEAR')
587
588 call msg('ocebg_getFGEamplification', 'reading wind speed fields...')
589 do timeStepIndex = 1, ndaysWinds * 24 / timeStepWinds
590 call tim_dateStampToYYYYMMDDHH(dataStampList(timeStepIndex), hour, day, monthNumber, &
591 ndays, yyyy, verbose_opt = .False.)
592 call msg('ocebg_getFGEamplification', ' '//str(timeStepIndex)//str(dataStampList(timeStepIndex))//str(yyyy)//str(monthNumber)//str(day)//str(hour))
593 call gio_readFromFile(stateVector, './winds', ' ', ' ', stepIndex_opt = timeStepIndex, &
594 containsFullField_opt=.true.)
595 end do
596
597 call gsv_getField(stateVector, uu_ptr4d, 'UU')
598 call gsv_getField(stateVector, vv_ptr4d, 'VV')
599 call gsv_getField(stateVectorAmplFactor, stateVectorAmplFactor_ptr)
600
601 write(*,*)
602 call msg('ocebg_getFGEamplification', 'detecting tropical storms (TS)...')
603
604 time_loop: do timeStepIndex = 1, ndaysWinds * 24 / timeStepWinds
605
606 call tim_dateStampToYYYYMMDDHH(dataStampList(timeStepIndex), hour, day, monthNumber, &
607 ndays, yyyy, verbose_opt = .False.)
608 do monthIndex = 1, nmonthsExceptionNH
609 if (months(monthNumber) == monthExceptionNH(monthIndex)) then
610 maxLatNH = maxLatExceptionNH
611 call msg('ocebg_getFGEamplification', 'TS is allowed to penetrate up to '//str(maxLatNH) &
612 //' degrees North; current month: '//months(monthNumber))
613 end if
614 end do
615
616 do latIndex = myLatBeg, myLatEnd
617 do lonIndex = myLonBeg, myLonEnd
618
619 lat = real(hco%lat2d_4(lonIndex, latIndex), 8) * MPC_DEGREES_PER_RADIAN_R8
620 if ((lat > minLatNH .and. lat < maxLatNH) .or. (lat > minLatSH .and. lat < maxLatSH)) then
621 windSpeed = sqrt(uu_ptr4d(lonIndex, latIndex, 1, timeStepIndex) * &
622 uu_ptr4d(lonIndex, latIndex, 1, timeStepIndex) + &
623 vv_ptr4d(lonIndex, latIndex, 1, timeStepIndex) * &
624 vv_ptr4d(lonIndex, latIndex, 1, timeStepIndex))
625 amplFactor = max(1.0d0, min(4.0d0 * max(0.0d0,(windSpeed - 20.75d0)), 12.0d0))
626 ! stateVectorAmplFactor_ptr may already be assigned with a value > 1.0 in a previous state
627 if(amplFactor > stateVectorAmplFactor_ptr(lonIndex, latIndex, 1)) then
628 stateVectorAmplFactor_ptr(lonIndex, latIndex, 1) = amplFactor
629 end if
630 endif
631
632 end do
633 end do
634
635 end do time_loop
636
637 call gsv_deallocate(stateVector)
638
639 call gio_writeToFile(stateVectorAmplFactor, './amplification', 'ORIG')
640 call gsv_smoothHorizontal(stateVectorAmplFactor, smoothLenghtScale)
641 call gio_writeToFile(stateVectorAmplFactor, './amplification', 'SMOOTH')
642
643 end subroutine ocebg_getFGEamplification
644
645end module bgckOcean_mod
