1module burpFiles_mod
2 ! MODULE burpFiles_mod (prefix='brpf' category='3. Observation input/output')
3 !
4 !:Purpose: To store the filenames of the burp observation files and call
5 ! subroutines in readBurp_mod to read and update burp files.
6 !
7 use codePrecision_mod
8 use mathPhysConstants_mod
9 use utilities_mod
10 use obsSpaceData_mod
11 use burpread_mod
12 use bufr_mod
13 use utilities_mod
14 use obsSubSpaceData_mod
15 use burp_module
16 use obsUtil_mod
17 use obsVariableTransforms_mod
18
19 implicit none
20 save
21 private
22
23 ! public procedures
24 public :: brpf_getDateStamp, brpf_readfile, brpf_updatefile
25 public :: brpf_obsSub_read, brpf_obsSub_update
26
27contains
28
29 !--------------------------------------------------------------------------
30 ! brpf_getDataStamp
31 !--------------------------------------------------------------------------
32 subroutine brpf_getDateStamp(datestamp, burpFileName)
33 implicit none
34
35 ! Arguments:
36 integer, intent(out) :: dateStamp
37 character(len=*), intent(in) :: burpFileName
38
39 ! Locals:
40 integer :: ier, inblks, nulburp, fnom, fclos, numblks
41 logical :: isExist_L
42 integer :: ktime, kdate, kdate_recv, ktime_recv, ihandl, ilong
43 integer :: itime, iflgs, idburp, ilat, ilon, idx, idy
44 integer :: ialt, idelay, idate, irs, irunn, inblk, isup, ixaux
45 integer :: insup, inxaux
46 integer, allocatable :: ibuf(:)
47 integer :: inrecs, mrfcls, mrfopn, mrfopc, mrbhdr, mrfloc, mrfget, mrfmxl
48 integer :: istampobs, inewhh, newdate, nresume, ivals
49 real(8) :: delhh
50 character(len=9) :: clstnid
51
52 !
53 !- Get the date from the burp files
54 !
55
56 ier = mrfopc('MSGLVL','FATAL')
57
58 ivals = 8
59 kdate = -9999
60 ktime = -9999
61 nresume = 0
62 nulburp = 0
63 inquire(file=trim(burpFileName),exist=isExist_L)
64 if ( isExist_L ) then
65 ier = fnom(nulburp,trim(burpFileName),'RND+OLD',0)
66 write(*,*)' Open File : ',trim(burpFileName)
67 if ( ier == 0 ) then
68 inblks = -1
69 inblks = numblks(nulburp)
70 if ( inblks > 0 ) then
71 inrecs= mrfopn(nulburp,'READ')
72 ilong = mrfmxl(nulburp)
73 allocate(ibuf(ilong + 20))
74 ibuf(1) = ilong + 20
75 ihandl = mrfloc(nulburp,0,'>>*******',-1,-1,-1,-1,-1,-1,0)
76 if ( ihandl < 0 ) then
77 ihandl=mrfloc(nulburp,0,'*********',-1,-1,-1,-1,-1,-1,0)
78 else
79 nresume=nresume+1
80 end if
81 if ( ihandl < 0 ) then
82 write(*,*) 'AUCUN ENREGISTREMENT VALIDE DANS LE FICHIER BURP'
83 else
84 if ((kdate < 0.and.ktime < 0).or.nresume == 1) then
85 insup=0
86 inxaux=0
87 ier=mrfget(ihandl,ibuf)
88 ier=mrbhdr(ibuf,itime,iflgs,clstnid,idburp,ilat, &
89 ilon,idx,idy, ialt,idelay,idate,irs,irunn,inblk, &
90 isup,insup,ixaux,inxaux)
91 ktime=itime
92 kdate=idate
93 if (nresume == 1) nresume=2
94 end if
95 end if
96 deallocate(ibuf)
97 ier=mrfcls(nulburp)
98 end if
99 end if
100 ier= fclos(nulburp)
101 end if
102
103 !
104 !- Set reference datestamp
105 !
106 ! Make sure all mpi tasks have a valid date (important for split burp files)
107 call rpn_comm_allreduce(kdate,kdate_recv,1,"MPI_INTEGER","MPI_MAX","GRID",ier)
108 call rpn_comm_allreduce(ktime,ktime_recv,1,"MPI_INTEGER","MPI_MAX","GRID",ier)
109 kdate = kdate_recv
110 ktime = ktime_recv
111 if (nresume >= 1 ) then
112 ier = newdate(datestamp,kdate,ktime*10000,3)
113 else
114 ! Assumes 6-hour windows with reference times being synoptic times.
115 ! Does not require kdate and ktime to be from a resume record.
116 ier = newdate(istampobs,kdate,ktime*10000,3)
117 delhh = 3.0d0
118 call incdatr (datestamp, istampobs, delhh)
119 ier = newdate(datestamp,kdate,inewhh,-3)
120 ktime = ktime/100
121 if (ktime >= 21 .or. ktime < 3) then
122 ktime = 0
123 else if(ktime >= 3 .and. ktime < 9) then
124 ktime = 6
125 else if(ktime >= 9 .and. ktime < 15) then
126 ktime = 12
127 else
128 ktime = 18
129 end if
130 ier = newdate(datestamp,kdate,ktime*1000000,3)
131 ktime = ktime*100
132 end if
133
134 write(*,*)' BURP FILES VALID DATE (YYYYMMDD) : ', kdate
135 write(*,*)' BURP FILES VALID TIME (HHMM) : ', ktime
136 write(*,*)' BURP FILES DATESTAMP : ', datestamp
137
138 end subroutine brpf_getDateStamp
139
140 !--------------------------------------------------------------------------
141 ! brpf_readFile
142 !--------------------------------------------------------------------------
143 subroutine brpf_readFile(obsdat,fileName,familyType,fileIndex)
144 implicit none
145
146 ! Arguments:
147 type (struct_obs), intent(inout) :: obsdat
148 character(len=*), intent(in) :: fileName
149 character(len=*), intent(in) :: familyType
150 integer, intent(in) :: fileIndex
151
152 ! Locals:
153 integer :: bodyIndex, bodyIndexBegin, bodyIndexEnd, headerIndexBegin, headerIndexEnd, headerIndex
154 integer :: numBody, numHeader
155 logical :: burp_chem
156 real(pre_obsReal) :: missingValue
157
158 write(*,*) ' '
159 write(*,*) 'brpf_readFile: Starting'
160 write(*,*) ' '
161 missingValue = real(MPC_missingValue_R8,pre_obsReal)
162
163 bodyIndexBegin = obs_numbody(obsdat) + 1
164 headerIndexBegin = obs_numheader(obsdat) + 1
165 call brpr_readBurp(obsdat, & ! INOUT
166 familyType, fileName, fileIndex) ! IN
167 bodyIndexEnd = obs_numbody(obsdat)
168 headerIndexEnd = obs_numheader(obsdat)
169
170 burp_chem = trim(familyType) == 'CH'
171
172 if ( trim(familyType) /= 'TO' .and. .not.burp_chem ) then
173
174 call ovt_transformObsValues (obsdat, headerIndexBegin, headerIndexEnd )
175 call ovt_adjustHumGZ (obsdat, headerIndexBegin, headerIndexEnd )
176 call obsu_computeVertCoordSurfObs(obsdat, headerIndexBegin, headerIndexEnd )
177
178 end if
179
180 do headerIndex = headerIndexBegin, headerIndexEnd
181
182 call obs_headSet_i(obsdat, OBS_IDF, headerIndex, fileIndex)
183 call obs_setFamily(obsdat, trim(familyType), headerIndex)
184
185 ! For CH family, apply scaling from the element BUFR_SCALE_EXPONENT when present.
186 if (burp_chem) call brpf_setScaleCH(obsdat, headerIndex, forward=.true.)
187
188 end do
189
190 ! initializations
191 do bodyIndex = bodyIndexBegin, bodyIndexEnd
192
193 if ( obs_columnActive_RB(obsdat, OBS_OMA) ) call obs_bodySet_r(obsdat, OBS_OMA , bodyIndex, missingValue )
194 if ( obs_columnActive_RB(obsdat, OBS_OMA0)) call obs_bodySet_r(obsdat, OBS_OMA0, bodyIndex, missingValue )
195 if ( obs_columnActive_RB(obsdat, OBS_OMP) ) call obs_bodySet_r(obsdat, OBS_OMP , bodyIndex, missingValue )
196 if ( obs_columnActive_RB(obsdat, OBS_OMP6)) call obs_bodySet_r(obsdat, OBS_OMP6, bodyIndex, missingValue )
197 if ( obs_columnActive_RB(obsdat, OBS_OER) ) call obs_bodySet_r(obsdat, OBS_OER , bodyIndex, missingValue )
198 if ( obs_columnActive_RB(obsdat, OBS_HPHT)) call obs_bodySet_r(obsdat, OBS_HPHT, bodyIndex, missingValue )
199 if ( obs_columnActive_RB(obsdat, OBS_HAHT)) call obs_bodySet_r(obsdat, OBS_HAHT, bodyIndex, missingValue )
200 if ( obs_columnActive_RB(obsdat, OBS_WORK)) call obs_bodySet_r(obsdat, OBS_WORK, bodyIndex, missingValue )
201 if ( obs_columnActive_RB(obsdat, OBS_SIGI)) call obs_bodySet_r(obsdat, OBS_SIGI, bodyIndex, missingValue )
202 if ( obs_columnActive_RB(obsdat, OBS_SIGO)) call obs_bodySet_r(obsdat, OBS_SIGO, bodyIndex, missingValue )
203 if ( obs_columnActive_RB(obsdat, OBS_ZHA )) call obs_bodySet_r(obsdat, OBS_ZHA , bodyIndex, missingValue )
204
205 end do
206
207 ! For GP family, initialize OBS_OER to element 15032 (ZTD formal error)
208 ! for all ZTD data (element 15031)
209 if ( trim(familyType) == 'GP') then
210 write(*,*)' Initializing OBS_OER for GB-GPS ZTD to formal error (ele 15032)'
211 call obsu_setGbgpsError(obsdat, headerIndexBegin, headerIndexEnd )
212 end if
213
214 numHeader = obs_numHeader(obsdat)
215 numBody = obs_numBody(obsdat)
216 write(*,*) 'brpf_readFile: obs_numheader', numHeader
217 write(*,*) 'brpf_readFile: obs_numbody ', numBody
218
219 end subroutine brpf_readFile
220
221 !--------------------------------------------------------------------------
222 ! brpf_updateFile
223 !--------------------------------------------------------------------------
224 subroutine brpf_updateFile(obsSpaceData,fileName,familyType,fileIndex)
225 implicit none
226
227 ! Arguments:
228 type (struct_obs), intent(inout) :: obsSpaceData
229 character(len=*), intent(in) :: fileName
230 character(len=*), intent(in) :: familyType
231 integer, intent(in) :: fileIndex
232
233 ! Locals:
234 integer :: headerIndex
235
236 call utl_tmg_start(12,'----UpdateBurpFile')
237
238 write(*,*)
239 write(*,*) 'brpf_updateFile: Starting'
240
241 ! CH family: Scaling of the obs related values to be stored in the BURP files
242 if (familytype == 'CH') then
243 call obs_set_current_header_list(obsSpaceData,'CH')
244 HEADER: do
245 headerIndex = obs_getHeaderIndex(obsSpaceData)
246 if (headerIndex < 0) exit HEADER
247 call brpf_setScaleCH(obsSpaceData,headerIndex,forward=.false.)
248 end do HEADER
249 end if
250
251 call brpr_updateBurp(obsSpaceData,familyType,fileName,fileIndex)
252
253 write(*,*) 'brpf_updateFile: Done'
254 write(*,*)
255
256 call utl_tmg_stop(12)
257
258 end subroutine brpf_updateFile
259
260 !--------------------------------------------------------------------------
261 ! brpf_setScaleCH
262 !--------------------------------------------------------------------------
263 subroutine brpf_setScaleCH(obsdat,headerIndex,forward)
264 !
265 !:Purpose: Apply or unapply scaling to CH observations by multiplying
266 ! (or dividing) with 10^{exponent} where the exponent is from
267 ! element BUFR_SCALE_EXPONENT if provided.
268 !
269 !
270 implicit none
271
272 ! Arguments:
273 type(struct_obs), intent(inout) :: obsdat ! struct_obs instance
274 integer, intent(in) :: headerIndex ! header index in obsdat
275 logical, intent(in) :: forward ! applies scaling if .true., unapplies scaling if .false.
276
277 ! Locals:
278 integer :: bodyIndex,rln,nlv
279 real(pre_obsReal) :: obsv
280 real(pre_obsReal) :: vomp, voma, voer, vhpht, scale
281 integer :: nexp,iobs,iexp
282 real(pre_obsReal), allocatable :: expnt(:)
283
284 rln = obs_headElem_i(obsdat,OBS_RLN,headerIndex)
285 nlv = obs_headElem_i(obsdat,OBS_NLV,headerIndex)
286
287 allocate(expnt(nlv))
288
289 ! Count number of power of 10 exponents
290 nexp = 0
291 do bodyIndex = rln, nlv + rln -1
292 if (obs_bodyElem_i(obsdat,OBS_VNM,bodyIndex) == bufr_scale_exponent) then
293 nexp = nexp + 1
294 expnt(nexp) = obs_bodyElem_r(obsdat,OBS_VAR,bodyIndex)
295 end if
296 end do
297
298 if (nexp == 0) then
299 deallocate(expnt)
300 return
301 end if
302
303 if (nexp*2 /= nlv) then
304 ! Skip over obs assuming mantissa was filtered out in brpr_readBurp
305 ! (not inserted in obsSpaceData) due to quality flags.
306 ! Set exponent quality flag to that of a 'Suspicious element'
307
308 do bodyIndex = RLN, NLV + RLN -1
309 call obs_bodySet_r(obsdat,OBS_VAR,bodyIndex, 0.0D0 )
310 call obs_bodySet_i(obsdat,OBS_FLG,bodyIndex, ibset(obs_bodyElem_i(obsdat,OBS_FLG,bodyIndex),02) )
311 call obs_bodySet_i(obsdat,OBS_FLG,bodyIndex, ibset(obs_bodyElem_i(obsdat,OBS_FLG,bodyIndex),04) )
312 call obs_bodySet_i(obsdat,OBS_FLG,bodyIndex, ibset(obs_bodyElem_i(obsdat,OBS_FLG,bodyIndex),09) )
313 end do
314
315 ! write(*,*) 'NLV =',nlv,' Nexp=',nexp
316 ! call utl_abort('brpf_setScaleCH: Inconsistent number of exponents')
317 deallocate(expnt)
318 return
319 end if
320
321 if (forward) then
322
323 ! Apply power of 10 exponents if present
324 iobs = 0
325 do bodyIndex = RLN, NLV + RLN -1
326 if (obs_bodyElem_i(obsdat,OBS_VNM,bodyIndex) /= bufr_scale_exponent) then
327 iobs = iobs + 1
328 obsv = obs_bodyElem_r(obsdat,OBS_VAR,bodyIndex)
329 call obs_bodySet_r(obsdat,OBS_VAR,bodyIndex,obsv*10**(expnt(iobs)) )
330 end if
331 end do
332
333 else
334
335 ! Unapply power of 10 exponents if present
336 iobs=0
337 iexp=0
338 do bodyIndex = RLN, NLV + RLN -1
339 if (obs_bodyElem_i(obsdat,OBS_VNM,bodyIndex) == bufr_scale_exponent) then
340 ! Store scaling exponents
341 iexp = iexp + 1
342 call obs_bodySet_r(obsdat,OBS_OMP,bodyIndex,expnt(iexp))
343 call obs_bodySet_r(obsdat,OBS_OMA,bodyIndex,expnt(iexp))
344 call obs_bodySet_r(obsdat,OBS_OER,bodyIndex,expnt(iexp))
345 call obs_bodySet_r(obsdat,OBS_HPHT,bodyIndex,expnt(iexp))
346 else
347 iobs=iobs+1
348 obsv=obs_bodyElem_r(obsdat,OBS_VAR,bodyIndex)
349 vomp=obs_bodyElem_r(obsdat,OBS_OMP,bodyIndex)
350 voma=obs_bodyElem_r(obsdat,OBS_OMA,bodyIndex)
351 voer=obs_bodyElem_r(obsdat,OBS_OER,bodyIndex)
352 vhpht=obs_bodyElem_r(obsdat,OBS_HPHT,bodyIndex)
353 scale=10**(-expnt(iobs))
354 call obs_bodySet_r(obsdat,OBS_VAR,bodyIndex,obsv*scale )
355 call obs_bodySet_r(obsdat,OBS_OMP,bodyIndex,vomp*scale )
356 call obs_bodySet_r(obsdat,OBS_OMA,bodyIndex,voma*scale )
357 call obs_bodySet_r(obsdat,OBS_OER,bodyIndex,voer*scale )
358 call obs_bodySet_r(obsdat,OBS_HPHT,bodyIndex,vhpht*scale )
359 end if
360 end do
361
362 end if
363
364 end subroutine brpf_setScaleCH
365
366 !--------------------------------------------------------------------------
367 ! brpf_obsSub_read
368 !--------------------------------------------------------------------------
369 function brpf_obsSub_read(filename,stnid,varno,nlev,ndim,block_type,bkstp_opt, &
370 match_nlev_opt,codtyp_opt,numColumns_opt) result(burp_out)
371 !
372 !:Purpose: To retrieve information from observation BURP file. Returns the
373 ! data in a struct_oss_obsdata object. Can retrieve either 1D or 2D
374 ! data from a report.
375 !
376 !:Comments:
377 !
378 ! - BUFR power 10 exponent element (i.e. data with BUFR number
379 ! BUFR_SCALE_EXPONENT) will be applied only to 1D data if present.
380 ! - As burp_out is for a specific input stnid, burp_out%code contains
381 ! only the (lat/long and time coord.) with 22 characters.
382 ! - Exponent BUFR data (i.e. data with BUFR number BUFR_SCALE_EXPONENT)
383 ! will be applied only to 1D data.
384 !
385 !:Arguments:
386 ! :filename: observation family name
387 ! :stnid: station ID of observation
388 ! :varno: BUFR code (if <=0, search through all codes to obtain first
389 ! between 10000 and 16000)
390 ! :nlev: number of levels in the observation (number of rows)
391 ! :ndim: number of dimensions for the retrieved data in each report
392 ! (e.g. ndim=1 for std, ndim=2 for averaging kernels)
393 ! :numColumns_opt: Number of columns (if different from nlev and for ndim=2)
394 ! :block_type: block type indicated by the two rightmost bits of bknat.
395 ! Valid values are 'DATA', 'INFO', '3-D', and 'MRQR'.
396 ! :match_nlev_opt: =.true. (default) causes filtering out of report if
397 ! the report number of levels is different from the input
398 ! argument nlev
399 !
400 implicit none
401
402 ! Arguments:
403 character(len=*), intent(in) :: filename ! BURP file name
404 character(len=9), intent(in) :: stnid ! station ID of observation
405 integer, intent(in) :: varno
406 integer, intent(in) :: nlev ! number of levels in the observation
407 integer, intent(in) :: ndim
408 character(len=4), intent(in) :: block_type
409 integer, optional, intent(in) :: numColumns_opt ! Number of columns (if different from nlev and for ndim=2)
410 integer, optional, intent(in) :: bkstp_opt ! bkstp number of requested block
411 logical, optional, intent(in) :: match_nlev_opt
412 integer, optional, intent(in) :: codtyp_opt(:) ! optional CODTYP list for search
413 ! Result:
414 type(struct_oss_obsdata) :: burp_out ! struct_oss_obsdata object
415
416 ! Locals:
417 character(len=9) :: rep_stnid
418 type(burp_file) :: brp
419 type(burp_rpt) :: rep
420 type(burp_block) :: blk
421 integer :: error,ref_rpt,nrep,ref_blk,varno_ivar
422 integer :: ref_bkstp,nval,ivar,ilev,icount,icodtyp
423 integer :: date,time,ilat,ilon,iele,nele,icol,ivar_previous
424 logical :: match_nlev
425 real(8) :: exponent
426
427 if (present(match_nlev_opt)) then
428 match_nlev = match_nlev_opt
429 else
430 match_nlev = .true.
431 end if
432
433 ! initialize burp file, report, and block
434 call BURP_Init(brp, iostat=error)
435 call BURP_Init(rep, iostat=error)
436 call BURP_Init(blk, iostat=error)
437
438 ! open the burp file
439 call BURP_New(brp, FILENAME=filename, MODE=FILE_ACC_READ, IOSTAT=error)
440 if (error /= 0) call utl_abort('brpf_obsSub_read: Could not find/open BURP file: ' // trim(filename))
441
442 write(*,*) "brpf_obsSub_read: Reading file " // trim(filename)
443 write(*,*) "brpf_obsSub_read: Selecting STNID = ",stnid," BUFR = ",varno," block type = ",block_type
444 write(*,*) "brpf_obsSub_read: nlev = ",nlev," match_nlev = ",match_nlev
445 if (present(bkstp_opt)) write(*,*) "brpf_obsSub_read: bkstp = ",bkstp_opt
446 if (present(codtyp_opt)) write(*,*) "brpf_obsSub_read: codtyp = ",codtyp_opt(:)
447
448 ! get number of reports in file
449 call BURP_Get_Property(brp, NRPTS=nrep)
450
451 ! allocate memory
452 if (ndim == 1) then
453 call oss_obsdata_alloc(burp_out,nrep,dim1=nlev)
454 else
455 if (present(numColumns_opt)) then
456 call oss_obsdata_alloc(burp_out,nrep,dim1=nlev,dim2_opt=numColumns_opt)
457 else
458 call oss_obsdata_alloc(burp_out,nrep,dim1=nlev,dim2_opt=nlev)
459 end if
460 end if
461
462 icount = 0 ! counter of reports with same stnid, number of levels, and varno as input
463 ref_rpt = 0
464
465 ! loop through reports
466 REPORTS: do
467
468 ref_rpt = BURP_Find_Report(brp, REPORT=rep, SEARCH_FROM=ref_rpt, IOSTAT=error)
469
470 if (ref_rpt<0) exit REPORTS
471
472 call BURP_Get_Property(rep, STNID=rep_stnid, DATE=date, TEMPS=time, LATI=ilat, LONG=ilon, IDTYP=icodtyp)
473
474 if (present(codtyp_opt)) then
475 if (.not.any(codtyp_opt(:) == icodtyp)) cycle REPORTS
476 end if
477
478 if (.not.utl_stnid_equal(stnid,rep_stnid)) cycle REPORTS
479
480 ! loop through blocks
481 ref_blk = 0
482 BLOCKS: do
483
484 ref_blk = BURP_Find_Block(rep, BLOCK=blk, SEARCH_FROM=ref_blk, IOSTAT=error)
485 if (ref_blk<0) exit BLOCKS
486
487 call BURP_Get_Property(blk, NELE=nele, NVAL=nval, BKSTP=ref_bkstp, IOSTAT=error)
488
489 if (.not.IS_Burp_Btyp(trim(block_type),BLOCK=blk)) cycle BLOCKS
490 if (match_nlev .and. nval /= nlev) cycle BLOCKS
491 if (present(bkstp_opt)) then
492 if (bkstp_opt /= ref_bkstp) cycle BLOCKS
493 end if
494
495 if (varno > 0) then
496 ivar = BURP_Find_Element(blk, ELEMENT=varno, IOSTAT=error)
497 if (ivar < 0) cycle BLOCKS
498 else
499 ! Search for first data element within elements 10000 and 16000.
500 varno_ivar=-1
501 do ivar=1,nele
502 varno_ivar=BURP_Get_Element(blk, INDEX=ivar, IOSTAT=error)
503 if (varno_ivar >= 10000.and.varno_ivar < 16000) exit
504 end do
505 if (varno_ivar < 10000.or.varno_ivar >= 16000) call utl_abort('brpf_obsSub_read: No valid element found for STNID ' // rep_stnid )
506 end if
507
508 ! required block found if code reaches this point, retrieve data and store in burp_out
509
510 if (nval > nlev) call utl_abort('brpf_obsSub_read: number of levels in the report (' // trim(utl_str(nval)) // &
511 ') exceeds the specified maximum number of levels (' // trim(utl_str(nlev)) // &
512 ') for STNID ' // rep_stnid )
513
514 icount=icount+1
515 burp_out%code(icount) = oss_obsdata_get_header_code(ilon,ilat,date,time,rep_stnid) ! this code is a unique identifier for this report
516
517 if (ndim == 1) then
518 ! retrieve 1D data
519
520 do ilev=1,nval
521 burp_out%data1d(ilev,icount) = BURP_Get_Rval(blk, NELE_IND=ivar, NVAL_IND=ilev, NT_IND=1, IOSTAT=error)
522 end do
523
524 if (ivar < nele) then
525 if (BURP_Get_Element(blk, INDEX=ivar+1, IOSTAT=error) == BUFR_SCALE_EXPONENT) then
526 ! Apply exponent
527 do ilev=1,nval
528 exponent = BURP_Get_Rval(blk, NELE_IND=ivar+1, NVAL_IND=ilev, NT_IND=1, IOSTAT=error)
529 burp_out%data1d(ilev,icount) = burp_out%data1d(ilev,icount) * 10**exponent
530 end do
531 end if
532 end if
533
534 else if (ndim == 2) then
535 ! retrieve 2D data
536
537 icol = 0
538 ivar_previous=0
539 do iele=1,nele
540 ivar = BURP_Get_Element(blk, INDEX=iele, IOSTAT=error)
541 if (ivar == varno) then
542 icol = icol+1
543 ivar_previous=ivar
544 do ilev=1,nval
545 burp_out%data2d(ilev,icol,icount) = BURP_Get_Rval(blk, NELE_IND=iele, NVAL_IND=ilev, NT_IND=1, IOSTAT=error)
546 end do
547 else if (ivar_previous == varno .and. ivar == BUFR_SCALE_EXPONENT) then
548 ivar_previous=0
549 do ilev=1,nval
550 exponent = BURP_Get_Rval(blk, NELE_IND=iele, NVAL_IND=ilev, NT_IND=1, IOSTAT=error)
551 burp_out%data2d(ilev,icol,icount) = burp_out%data2d(ilev,icol,icount) * 10**exponent
552 end do
553 else
554 ivar_previous=0
555 end if
556 end do
557 if (present(numColumns_opt)) then
558 if (icol /= numColumns_opt) call utl_abort('brpf_obsSub_read: number of columns (' // trim(utl_str(icol)) // &
559 ') is not equal to the required number (' // trim(utl_str(numColumns_opt)) // &
560 ') for STNID ' // rep_stnid )
561 else
562 if (icol > nlev ) call utl_abort('brpf_obsSub_read: number of columns (' // trim(utl_str(icol)) // &
563 ') exceeds the maximum number (' // trim(utl_str(nlev)) // &
564 ') for STNID ' // rep_stnid )
565 end if
566 end if
567
568 exit BLOCKS
569
570 end do BLOCKS
571
572 end do REPORTS
573
574 ! resize first dimension of data arrays from length of nrep to icount
575 call utl_resize(burp_out%code,icount)
576 if (ndim == 1) then
577 call utl_resize(burp_out%data1d,nlev,icount)
578 else if (ndim == 2) then
579 call utl_resize(burp_out%data2d,nlev,nlev,icount)
580 end if
581
582 burp_out%nrep = icount
583
584 write(*,*) "brpf_obsSub_read: Reading of file complete. Number of reports found: ",burp_out%nrep
585
586 ! deallocate
587 Call BURP_Free(brp,iostat=error)
588 Call BURP_Free(rep,iostat=error)
589 Call BURP_Free(blk,iostat=error)
590
591 end function brpf_obsSub_read
592
593 !--------------------------------------------------------------------------
594 ! brpf_obsSub_update
595 !--------------------------------------------------------------------------
596 function brpf_obsSub_update(obsdata,filename,varno,block_type,bkstp_opt, &
597 multi_opt) result(nrep_modified)
598 !
599 !:Purpose: To add or modify data in BURP files from data stored in a
600 ! struct_oss_obsdata object. Provided data can be either 1D or 2D.
601 !
602 !:Comments:
603 ! - Currently assumes that all elements of varno(:) are distinct from
604 ! each other.
605 ! - Blocks with new data will overwrite any existing data of the same
606 ! varno. Otherwise the new data will be appended to the block.
607 !
608 !:Arguments:
609 ! :varno: BUFR descriptors. Number of elements must be
610 ! max(1,obsdata%dim2)
611 ! :block_type: block type indicated by the two rightmost bits of bknat.
612 ! Valid values are 'DATA', 'INFO', '3-D', and 'MRQR'.
613 ! :multi_opt: Indicates if intended report are for 'UNI' or 'MULTI'
614 ! level data (description is not accurate)
615 !
616 implicit none
617
618 ! Arguments:
619 type(struct_oss_obsdata), intent(inout) :: obsdata ! Input struct_oss_obsdata object for varno
620 character(len=*), intent(in) :: filename ! BURP file name
621 character(len=4), intent(in) :: block_type
622 integer, intent(in) :: varno(:)
623 integer, optional, intent(in) :: bkstp_opt ! bkstp number of requested block
624 character(len=*), optional, intent(in) :: multi_opt
625 ! Result:
626 integer :: nrep_modified ! Number of modified reports
627
628 ! Locals:
629 integer :: ncount
630 logical :: blk_found
631 integer, parameter :: LNMX=100000, code_len=90
632 character(len=9) :: stnid
633 character(len=code_len) :: code ! Must be at least as large as burp_code_len
634 type(burp_file) :: brp
635 type(burp_rpt) :: rep,rep_new
636 type(burp_block) :: blk
637 integer :: error,ref_rpt,nrep,ref_blk,ndim,dim1,dim2
638 integer :: ref_bkstp,nval,ivar,ilev,istat
639 integer :: date,time,ilat,ilon,iele,nele,k
640 integer, allocatable :: address(:)
641 real(4), allocatable :: new_vals(:,:,:)
642 logical, allocatable :: modify(:)
643
644 ! Check presence of data to update
645 if (obsdata%nrep <= 0) then
646 write(*,*) 'brpf_obsSub_update: Skipped due to absence of data to update.'
647 return
648 end if
649
650 ! Identify dimensions for the input data
651 ndim=obsdata%ndim
652 dim1=obsdata%dim1
653 if (ndim == 1) then
654 dim2=1
655 else
656 dim2=obsdata%dim2
657 end if
658
659 if (size(varno) < dim2) call utl_abort('brpf_obsSub_update: Number of BUFR elements not sufficient. ' // &
660 trim(utl_str(size(varno))) // ' vs ' // trim(utl_str(dim2)))
661
662 if (code_len < oss_obsdata_code_len()) call utl_abort('brpf_obsSub_update: Length of code string' // &
663 ' needs to be increased to ' // &
664 trim(utl_str(oss_obsdata_code_len())))
665
666 ! initialize burp file, report, and block system resources
667 call BURP_Init(brp, iostat=error)
668 call BURP_Init(rep, R2=rep_new, iostat=error)
669 call BURP_Init(blk, iostat=error)
670
671 ! open the burp file in append mode (to replace or add data in a block)
672 call BURP_New(brp, FILENAME=filename, MODE=FILE_ACC_APPEND, IOSTAT=error)
673 if (error /= 0) call utl_abort('brpf_obsSub_update: Could not open BURP file: ' // trim(filename))
674
675 ! get number of reports in file
676 call BURP_Get_Property(brp, NRPTS=nrep)
677
678 allocate(address(nrep),modify(nrep),new_vals(dim1,dim2,nrep))
679 address(:)=0
680 modify(:)=.false.
681 new_vals(:,:,:)=0.
682
683 ! First loop through reports to identify addresses of original file as well as identify if new
684 ! information should be included to that report.
685 ! NOTE: The addresses of all reports have to be saved in their original order to ensure the
686 ! order of the reports in the file is unchanged.
687 ref_rpt=0
688 ncount=0
689 obsdata%irep=1
690 REPORTS1: do
691
692 ref_rpt = BURP_Find_Report(brp, REPORT=rep, SEARCH_FROM=ref_rpt, IOSTAT=error)
693 if (ref_rpt<0) exit REPORTS1
694
695 ncount=ncount+1
696 address(ncount)=ref_rpt
697
698 call BURP_Get_Property(rep, STNID=stnid, DATE=date, TEMPS=time, LATI=ilat, LONG=ilon)
699
700 if (stnid(1:2) == '>>') cycle REPORTS1
701
702 ! Get unique identifier for search from input data
703 code = oss_obsdata_get_header_code(ilon,ilat,date,time,stnid)
704
705 ! Determine if replacement/additional data likely present for this report
706 if (dim1 == 1.and.dim2 == 1) then
707 new_vals(1,1,ncount)=real(oss_obsdata_get_element(obsdata,trim(code),1,stat_opt=istat))
708 else if (dim2 == 1) then
709 new_vals(:,1,ncount)=real(oss_obsdata_get_array1d(obsdata,trim(code),stat_opt=istat))
710 else
711 new_vals(:,:,ncount)=real(oss_obsdata_get_array2d(obsdata,trim(code),stat_opt=istat))
712 end if
713
714 if (istat == 0) then
715 if (present(multi_opt)) then
716 ! loop through blocks to find first data block
717 ref_blk = 0
718 BLOCKS1: do
719 ref_blk = BURP_Find_Block(rep, BLOCK=blk, SEARCH_FROM=ref_blk, IOSTAT=error)
720 if (ref_blk<0) exit BLOCKS1
721 if (IS_Burp_Btyp('DATA',BLOCK=blk)) then
722 if (IS_Burp_Btyp(trim(multi_opt),BLOCK=blk)) modify(ncount) = .true.
723 exit BLOCKS1
724 end if
725 end do BLOCKS1
726 else
727 modify(ncount) = .true.
728 end if
729 end if
730
731 end do REPORTS1
732
733 nrep_modified = count(modify) ! number of reports with same code and, possibly, same number of obs data levels
734
735 ! Generate new report
736 Call BURP_New(rep_new, ALLOC_SPACE=10*LNMX, IOSTAT=error)
737
738 ! second loop through reports to include the new information to the file
739 REPORTS2: do k=1,ncount
740
741 call BURP_Get_Report(brp, REPORT=rep, REF=address(k), IOSTAT=error)
742
743 ! Copy report header
744 Call BURP_Copy_Header(TO=rep_new,FROM=rep)
745 Call BURP_Init_Report_Write(brp,rep_new,IOSTAT=error)
746
747 ! Ignore resume records in BLOCKS do loop
748 call BURP_Get_Property(rep, STNID=stnid)
749
750 if (stnid(1:2) == '>>') then
751 write(*,*) 'brpf_obsSub_update: reading resume record, skip reading blocks'
752 else
753 ! loop through blocks
754 ref_blk = 0
755 BLOCKS: do
756
757 ref_blk = BURP_Find_Block(rep, BLOCK=blk, SEARCH_FROM=ref_blk, IOSTAT=error)
758 if (ref_blk<0) exit BLOCKS
759
760 if (modify(k)) then
761
762 call BURP_Get_Property(blk, NELE=nele, NVAL=nval, BKSTP=ref_bkstp, IOSTAT=error)
763
764 blk_found = IS_Burp_Btyp(trim(block_type),BLOCK=blk) .and. dim1 == nval
765 if (present(bkstp_opt)) blk_found = blk_found .and. bkstp_opt == ref_bkstp
766
767 if (blk_found) then
768 ! Block to be modified has been found, add new data to block.
769 ! If the varno is already in the block, the new data will overwrite the
770 ! existing data, otherwise will append the new data to the block.
771
772 do iele=1,dim2
773 ivar = BURP_Find_Element(blk, ELEMENT=varno(iele), IOSTAT=error)
774 if (ivar < 0) then
775 ivar=nele+1
776 call BURP_Resize_Block(blk,ADD_NELE=1,IOSTAT=error)
777 call BURP_Set_Element(blk,NELE_IND=ivar,ELEMENT=varno(iele),IOSTAT=error)
778 end if
779
780 do ilev=1,nval
781 call BURP_Set_Rval(blk,NELE_IND=ivar,NVAL_IND=ilev,NT_IND=1,RVAL=new_vals(ilev,iele,k),IOSTAT=error)
782 end do
783 end do
784
785 end if
786 end if
787
788 ! The call to BURP_Write_Block has ENCODE_BLOCK and CONVERT_BLOCK set
789 ! to .true. in all cases, including when the block has not been
790 ! modified, due to problems that can occur when writing blocks
791 ! containing negative integers with datyp=4.
792 call BURP_Write_Block(rep_new, BLOCK=blk, ENCODE_BLOCK=.true., CONVERT_BLOCK=.true., IOSTAT=error)
793
794 end do BLOCKS
795 end if
796 call BURP_Delete_Report(brp,rep,IOSTAT=error)
797 call BURP_Write_Report(brp,rep_new,IOSTAT=error)
798
799 end do REPORTS2
800
801 ! deallocate
802 deallocate(address,modify,new_vals)
803 Call BURP_Free(brp,iostat=error)
804 Call BURP_Free(rep,R2=rep_new,iostat=error)
805 Call BURP_Free(blk,iostat=error)
806
807 end function brpf_obsSub_update
808
809end module burpFiles_mod
