IODA Bundle
gnssro_bufr2ioda.f90
Go to the documentation of this file.
1 !!!-------- README --------------------------------------------------------------------
2 ! this is a temporary routine to generate netcdf file
3 ! for jedi/ufo/gnssro/ operator test
4 ! Copyright UCAR 2019
5 ! Authors: Hailing Zhang
6 ! Mark Olah
7 
8 !!!--------- to run -----------------------------------------------------------------
9 ! ./gnssro_bufr2ioda yyyymmddhh $bufrfile_input $netcdffile_output
10 !
11 ! Return codes:
12 ! 0 - Success.
13 ! 1 - Unrecoverable system or logical error.
14 ! 2 - All provided observations are invalid. Cannot create NetCDF IODA file.
15 
16 
17 program gnssro_bufr2ioda
18 use netcdf
19 implicit none
20 
21 integer, parameter :: i_short = selected_int_kind(4) !2
22 integer, parameter :: i_kind = selected_int_kind(8)
23 integer, parameter :: r_double = selected_real_kind(15) !8
24 integer, parameter :: r_kind = selected_real_kind(15) !8
25 
26 ! output obs data stucture
27 integer :: ncid
28 integer :: nobs_dimid,nlocs_dimid,nvars_dimid,nrecs_dimid,ndatetime_dimid
29 integer :: varid_lat,varid_lon,varid_time,varid_datetime
30 integer :: varid_said,varid_siid,varid_ptid,varid_sclf,varid_asce,varid_ogce
31 integer :: varid_recn
32 integer :: varid_geoid, varid_rfict
33 integer :: varid_ref,varid_msl,varid_refoe
34 integer :: varid_bnd,varid_bndoe,varid_impp, varid_imph,varid_azim
35 integer :: nlev_dimid
36 integer :: varid_geo_temp,varid_geo_pres,varid_geo_shum,varid_geo_geop, varid_geo_geop_sfc
37 
38 character(len=256) :: infile, outfile
39 character,dimension(8) :: subset
40 character(len=10) :: anatime
41 integer(i_kind) :: ndatetime = 20
42 character(len=20) :: datetime
43 integer(i_kind) :: i,k,m,ireadmg,ireadsb,said,siid,ptid,sclf,asce,ogce
44 integer(i_kind) :: lnbufr = 10
45 integer(i_kind) :: nread,ndata,nvars,nrec, ndata0
46 integer(i_kind) :: anatime_i
47 integer(i_kind) :: idate5(6), idate,iadate5(6)
48 integer(i_kind) :: mincy,minobs
49 logical :: good,outside
50 integer :: refflag, bendflag
51 integer(i_kind),parameter :: mxib=31
52 integer(i_kind) :: ibit(mxib),nib
53 integer(i_kind),parameter :: maxlevs=500
54 integer(i_kind),parameter :: n1ahdr=13
55 integer(i_kind) :: maxobs
56 real(r_kind) :: timeo
57 type gpsro_type
58  integer(i_kind), allocatable, dimension(:) :: said
59  integer(i_kind), allocatable, dimension(:) :: siid
60  integer(i_kind), allocatable, dimension(:) :: sclf
61  integer(i_kind), allocatable, dimension(:) :: ptid
62  integer(i_kind), allocatable, dimension(:) :: recn
63  integer(i_kind), allocatable, dimension(:) :: asce
64  integer(i_kind), allocatable, dimension(:) :: ogce
65  real(r_double), allocatable, dimension(:) :: time
66  character(len=20), allocatable, dimension(:) :: datetime
67  real(r_double), allocatable, dimension(:) :: lat
68  real(r_double), allocatable, dimension(:) :: lon
69  real(r_double), allocatable, dimension(:) :: rfict
70  real(r_double), allocatable, dimension(:) :: azim
71  real(r_double), allocatable, dimension(:) :: geoid
72  real(r_double), allocatable, dimension(:) :: msl_alt
73  real(r_double), allocatable, dimension(:) :: ref
74  real(r_double), allocatable, dimension(:) :: refoe_gsi
75  real(r_double), allocatable, dimension(:) :: bend_ang
76  real(r_double), allocatable, dimension(:) :: impact_para
77  real(r_double), allocatable, dimension(:) :: bndoe_gsi
78 end type gpsro_type
79 
80 type(gpsro_type) :: gpsro_data
81 
82 real(r_double),dimension(n1ahdr) :: bfr1ahdr
83 real(r_double),dimension(50,maxlevs) :: data1b
84 real(r_double),dimension(50,maxlevs) :: data2a
85 real(r_double),dimension(maxlevs) :: nreps_this_roseq2
86 integer(i_kind) :: iret,levs,levsr,nreps_roseq1,nreps_roseq2_int
87 real(r_double) :: pcc,qfro,usage,dlat,dlat_earth,dlon,dlon_earth,freq_chk,freq,azim
88 real(r_double) :: height,rlat,rlon,ref,bend,impact,roc,geoid, bend_error,ref_error,bend_pccf,ref_pccf
89 real(r_double) :: obserr
90 real(r_double), parameter :: missingvalue=-9.9e10
91 logical, parameter :: globalmodel = .true. ! temporary
92 
93 character(10) nemo
94 character(80) hdr1a
95 
96 data hdr1a / 'YEAR MNTH DAYS HOUR MINU PCCF ELRC SAID SIID PTID GEODU SCLF OGCE' /
97 data nemo /'QFRO'/
98 nrec =0
99 ndata=0
100 nvars=2
101 maxobs=0
102 
103 call getarg(1,anatime)
104 call getarg(2,infile)
105 call getarg(3,outfile)
106 
107 read(anatime,'(i10)') anatime_i
108 read(anatime(1:4),'(i4)') iadate5(1)
109 read(anatime(5:6),'(i4)') iadate5(2)
110 read(anatime(7:8),'(i4)') iadate5(3)
111 read(anatime(9:10),'(i4)') iadate5(4)
112 iadate5(5) = 0
113 call w3fs21(iadate5,mincy)
114 
115 open(lnbufr,file=trim(infile),form='unformatted')
116 call openbf(lnbufr,'IN',lnbufr)
117 call datelen(10)
118 call readmg(lnbufr,subset,idate,iret)
119 if (iret/=0) then
120  write(6,*)'READ_GNSSRO: can not open gnssro file!'
121  stop
122 end if
123 
124 ! loop over all message to estimate maxobs
125 do while(ireadmg(lnbufr,subset,idate)==0)
126  read_loop_countmaxobs: do while(ireadsb(lnbufr)==0)
127  call ufbint(lnbufr,bfr1ahdr,n1ahdr,1,iret,hdr1a)
128  call ufbseq(lnbufr,data1b,50,maxlevs,levs,'ROSEQ1')
129  maxobs= maxobs + levs
130  enddo read_loop_countmaxobs
131 end do
132 
133 allocate(gpsro_data%said(maxobs))
134 allocate(gpsro_data%siid(maxobs))
135 allocate(gpsro_data%sclf(maxobs))
136 allocate(gpsro_data%ptid(maxobs))
137 allocate(gpsro_data%recn(maxobs))
138 allocate(gpsro_data%asce(maxobs))
139 allocate(gpsro_data%ogce(maxobs))
140 allocate(gpsro_data%time(maxobs))
141 allocate(gpsro_data%datetime(maxobs))
142 allocate(gpsro_data%lat(maxobs))
143 allocate(gpsro_data%lon(maxobs))
144 allocate(gpsro_data%rfict(maxobs))
145 allocate(gpsro_data%azim(maxobs))
146 allocate(gpsro_data%geoid(maxobs))
147 allocate(gpsro_data%msl_alt(maxobs))
148 allocate(gpsro_data%ref(maxobs))
149 allocate(gpsro_data%refoe_gsi(maxobs))
150 allocate(gpsro_data%bend_ang(maxobs))
151 allocate(gpsro_data%impact_para(maxobs))
152 allocate(gpsro_data%bndoe_gsi(maxobs))
153 
154 !rewind lnbufr
155 call closbf(lnbufr)
156 open(lnbufr,file=trim(infile),form='unformatted')
157 call openbf(lnbufr,'IN',lnbufr)
158 call datelen(10)
159 call readmg(lnbufr,subset,idate,iret)
160 
161 do while(ireadmg(lnbufr,subset,idate)==0)
162  read_loop: do while(ireadsb(lnbufr)==0)
163 ! Read/decode data in subset (profile)
164  call ufbint(lnbufr,bfr1ahdr,n1ahdr,1,iret,hdr1a)
165  call ufbint(lnbufr,qfro,1,1,iret,nemo)
166 ! observation time in minutes
167  idate5(1) = bfr1ahdr(1) ! year
168  idate5(2) = bfr1ahdr(2) ! month
169  idate5(3) = bfr1ahdr(3) ! day
170  idate5(4) = bfr1ahdr(4) ! hour
171  idate5(5) = bfr1ahdr(5) ! minute
172  idate5(6) = 0 ! seconds
173  pcc = bfr1ahdr(6) ! profile per cent confidence
174  roc = bfr1ahdr(7) ! Earth local radius of curvature
175  said = bfr1ahdr(8) ! Satellite identifier
176  siid = bfr1ahdr(9) ! Satellite instrument
177  ptid = bfr1ahdr(10) ! Platform transmitter ID number
178  geoid= bfr1ahdr(11) ! Geoid undulation
179  sclf = bfr1ahdr(12) ! Satellite classification
180  ogce = bfr1ahdr(13) ! Identification of originating/generating centre
181 
182  call w3fs21(idate5,minobs)
183  write(datetime,'(I4,"-",I2.2,"-",I2.2,"T",I2.2,":",I2.2,":",I2.2,"Z")') &
184  idate5(1),idate5(2),idate5(3),idate5(4),idate5(5),idate5(6)
185  timeo=real(minobs-mincy,r_kind)/60.0
186 
187  if( roc>6450000.0_r_kind .or. roc<6250000.0_r_kind .or. &
188  geoid>200_r_kind .or. geoid<-200._r_kind ) then
189  write(6,*)'READ_GNSSRO: profile fails georeality check, skip this report'
190  cycle read_loop
191  endif
192 
193 ! profile check: (1) CDAAC processing - cosmic-1, cosmic-2, sacc, cnofs, kompsat5
194  if ( ((said >= 740).and.(said <=745)).or.((said >= 750).and.(said <= 755)) &
195  .or.(said == 820).or.(said == 786).or.(said == 825) &
196  .or. ogce == 60) then !CDAAC processing
197  if(pcc == 0.0) then
198  write(6,*)'READ_GNSSRO: bad profile 0.0% confidence said=',said,'ptid=',ptid, ' SKIP this report'
199  cycle read_loop
200  endif
201  endif
202 
203  bendflag = 0
204  refflag = 0
205 ! profile check: (2) GRAS SAF processing - metopa-c, oceansat2, megha-tropiques, sacd
206  if ( (said >= 3 .and.said <= 5).or.(said == 421).or.(said == 440).or. (said == 821) ) then
207  call upftbv(lnbufr,nemo,qfro,mxib,ibit,nib)
208  if(nib > 0) then
209  do i = 1, nib
210  if(ibit(i)== 5) then ! bending angle
211  bendflag = 1
212  write(6,*)'READ_GNSSRO: bad profile said=',said,'ptid=',ptid, ' SKIP this report'
213  cycle read_loop
214  endif
215  if(ibit(i)== 6) then ! refractivity
216  refflag = 1
217  exit
218  endif
219  enddo
220  endif
221  endif
222 
223  asce = 0
224  call upftbv(lnbufr,nemo,qfro,mxib,ibit,nib)
225  if ( nib > 0) then
226  do i = 1, nib
227  if(ibit(i)== 3) then
228  asce = 1
229  exit
230  endif
231  enddo
232  end if
233 
234  call ufbint(lnbufr,nreps_this_roseq2,1,maxlevs,nreps_roseq1,'{ROSEQ2}')
235  call ufbseq(lnbufr,data1b,50,maxlevs,levs,'ROSEQ1')
236  call ufbseq(lnbufr,data2a,50,maxlevs,levsr,'ROSEQ3') ! refractivity
237 
238  nrec=nrec+1
239  ndata0 = ndata
240 
241  do k = 1, levs
242  rlat = data1b(1,k) ! earth relative latitude (degrees)
243  rlon = data1b(2,k) ! earth relative longitude (degrees)
244  azim = data1b(3,k)
245  height = data2a(1,k)
246  ref = data2a(2,k)
247  ref_error= data2a(4,k)
248  ref_pccf = data2a(6,k)
249 
250 ! Loop over number of replications of ROSEQ2 nested inside this particular replication of ROSEQ1
251  nreps_roseq2_int = nreps_this_roseq2(k)
252  do i = 1,nreps_roseq2_int
253  m = (6*i)-2
254  freq_chk=data1b(m,k) ! frequency (hertz)
255  if(nint(freq_chk).ne.0) cycle ! do not want non-zero freq., go on to next replication of ROSEQ2
256  freq = data1b(m,k)
257  impact = data1b(m+1,k) ! impact parameter (m)
258  bend = data1b(m+2,k) ! bending angle (rad)
259  bend_error = data1b(m+4,k) ! RMSE in bending angle (rad)
260  enddo
261 
262  bend_pccf=data1b((6*nreps_roseq2_int)+4,k) ! percent confidence for this ROSEQ1 replication
263  good=.true.
264 
265  if ( height<0._r_kind .or. height>60000._r_kind .or. &
266  abs(rlat)>90._r_kind .or. abs(rlon)>360._r_kind ) then
267  good=.false.
268  write(6,*)'READ_GNSSRO: obs fails georeality check, said=',said,'ptid=',ptid
269  endif
270  if ( bend>=1.e+9_r_kind .or. bend<=0._r_kind .or. impact>=1.e+9_r_kind .or. impact<roc .or. bendflag == 1 ) then
271  good=.false.
272  write(6,*)'READ_GNSSRO: obs bend/impact is invalid, said=',said,'ptid=',ptid
273  endif
274 
275  if ( ref>=1.e+9_r_kind .or. ref<=0._r_kind .or. refflag == 1 ) then
276  ref = missingvalue
277  endif
278 
279  if ( abs(azim)>360._r_kind .or. azim<0._r_kind ) then
280  azim = missingvalue
281  endif
282 
283  if(good) then
284  ndata = ndata +1
285  gpsro_data%recn(ndata) = nrec
286  gpsro_data%lat(ndata) = rlat
287  gpsro_data%lon(ndata) = rlon
288  gpsro_data%time(ndata) = timeo
289  gpsro_data%datetime(ndata) = datetime
290  gpsro_data%said(ndata) = said
291  gpsro_data%siid(ndata) = siid
292  gpsro_data%sclf(ndata) = sclf
293  gpsro_data%asce(ndata) = asce
294  gpsro_data%ptid(ndata) = ptid
295  gpsro_data%ogce(ndata) = ogce
296  gpsro_data%ref(ndata) = ref
297  gpsro_data%msl_alt(ndata) = height
298  gpsro_data%bend_ang(ndata) = bend
299  gpsro_data%bndoe_gsi(ndata) = bend_error
300  gpsro_data%impact_para(ndata) = impact!
301  gpsro_data%rfict(ndata) = roc
302  gpsro_data%geoid(ndata) = geoid
303  gpsro_data%azim(ndata) = azim
304  CALL bendingangle_err_gsi(rlat,impact-roc, obserr)
305  gpsro_data%bndoe_gsi(ndata) = obserr
306  if ( ref > missingvalue) then
307  CALL refractivity_err_gsi(rlat,height, globalmodel, obserr)
308  gpsro_data%refoe_gsi(ndata) = obserr
309  else
310  gpsro_data%refoe_gsi(ndata) = missingvalue
311  end if
312  end if
313 
314  end do ! end of k loop
315 
316  if (ndata == ndata0) nrec = nrec -1
317 
318  enddo read_loop
319 end do
320 
321 call closbf(lnbufr)
322 if (nrec==0) then
323  write(6,*) "Error. No valid observations found. Cannot create NetCDF ouput."
324  stop 2
325 endif
326 
327 call check( nf90_create(trim(outfile), nf90_netcdf4, ncid))
328 call check( nf90_def_dim(ncid, 'nlocs', ndata, nlocs_dimid) )
329 call check( nf90_def_dim(ncid, 'ndatetime', ndatetime, ndatetime_dimid) )
330 call check( nf90_put_att(ncid, nf90_global, 'date_time', anatime_i) )
331 call check( nf90_def_var(ncid, "latitude@MetaData", nf90_float, nlocs_dimid, varid_lat) )
332 call check( nf90_def_var(ncid, "longitude@MetaData", nf90_float, nlocs_dimid, varid_lon) )
333 call check( nf90_def_var(ncid, "time@MetaData", nf90_float, nlocs_dimid, varid_time) )
334 call check( nf90_put_att(ncid, varid_time, "longname", "time offset to analysis time" ))
335 call check( nf90_put_att(ncid, varid_time, "units", "hour" ))
336 call check( nf90_def_var(ncid, "datetime@MetaData", nf90_char, (/ ndatetime_dimid, nlocs_dimid /), varid_datetime) )
337 call check( nf90_def_var(ncid, "record_number@MetaData", nf90_int, nlocs_dimid, varid_recn))
338 call check( nf90_put_att(ncid, varid_recn, "longname", "GNSS RO profile identifier" ))
339 call check( nf90_def_var(ncid, "gnss_sat_class@MetaData", nf90_int, nlocs_dimid, varid_sclf))
340 call check( nf90_put_att(ncid, varid_sclf, "longname", "GNSS satellite classification, e.g, 401=GPS, 402=GLONASS" ))
341 call check( nf90_def_var(ncid, "reference_sat_id@MetaData", nf90_int, nlocs_dimid, varid_ptid))
342 call check( nf90_put_att(ncid, varid_ptid, "longname", "GNSS satellite transmitter identifier (1-32)" ))
343 call check( nf90_def_var(ncid, "occulting_sat_id@MetaData", nf90_int, nlocs_dimid, varid_said))
344 call check( nf90_put_att(ncid, varid_said, "longname", "Low Earth Orbit satellite identifier, e.g., COSMIC2=750-755" ))
345 call check( nf90_def_var(ncid, "occulting_sat_is@MetaData", nf90_int, nlocs_dimid, varid_siid))
346 call check( nf90_put_att(ncid, varid_siid, "longname", "satellite instrument"))
347 call check( nf90_def_var(ncid, "ascending_flag@MetaData", nf90_int, nlocs_dimid, varid_asce))
348 call check( nf90_put_att(ncid, varid_asce, "longname", "the original occultation ascending/descending flag" ))
349 call check( nf90_put_att(ncid, varid_asce, "valid_range", "0/descending or 1/ascending" ))
350 call check( nf90_def_var(ncid, "process_center@MetaData",nf90_int,nlocs_dimid, varid_ogce))
351 call check( nf90_put_att(ncid, varid_ogce, "longname", "originally data processing_center, &
352  e.g., 60 for UCAR, 94 for DMI, 78 for GFZ" ))
353 call check( nf90_def_var(ncid, "refractivity@ObsValue", nf90_float, nlocs_dimid, varid_ref) )
354 call check( nf90_put_att(ncid, varid_ref, "longname", "Atmospheric refractivity" ))
355 call check( nf90_put_att(ncid, varid_ref, "_FillValue", real(missingvalue) ))
356 call check( nf90_put_att(ncid, varid_ref, "units", "N" ))
357 call check( nf90_put_att(ncid, varid_ref, "valid_range", "0 - 500 N" ))
358 call check( nf90_def_var(ncid, "refractivity@ObsError", nf90_float, nlocs_dimid, varid_refoe))
359 call check( nf90_put_att(ncid, varid_refoe, "longname", "Input error in atmospheric refractivity" ))
360 call check( nf90_put_att(ncid, varid_refoe, "_FillValue", real(missingvalue) ))
361 call check( nf90_put_att(ncid, varid_refoe, "units", "N" ))
362 call check( nf90_put_att(ncid, varid_refoe, "valid_range", "0 - 10 N" ))
363 call check( nf90_def_var(ncid, "altitude@MetaData", nf90_float, nlocs_dimid, varid_msl) )
364 call check( nf90_put_att(ncid, varid_msl, "longname", "Geometric altitude" ))
365 call check( nf90_put_att(ncid, varid_msl, "units", "Meters" ))
366 call check( nf90_def_var(ncid, "bending_angle@ObsValue", nf90_float, nlocs_dimid, varid_bnd) )
367 call check( nf90_put_att(ncid, varid_bnd, "longname", "Bending Angle" ))
368 call check( nf90_put_att(ncid, varid_bnd, "units", "Radians" ))
369 call check( nf90_put_att(ncid, varid_bnd, "valid_range", "-0.001 - 0.08 Radians" ))
370 call check( nf90_def_var(ncid, "bending_angle@ObsError", nf90_float, nlocs_dimid, varid_bndoe) )
371 call check( nf90_put_att(ncid, varid_bndoe, "longname", "Input error in Bending Angle" ))
372 call check( nf90_put_att(ncid, varid_bndoe, "units", "Radians" ))
373 call check( nf90_put_att(ncid, varid_bndoe, "valid_range", "0 - 0.008 Radians" ))
374 call check( nf90_def_var(ncid, "impact_parameter@MetaData", nf90_float, nlocs_dimid, varid_impp))
375 call check( nf90_put_att(ncid, varid_impp, "longname", "distance from centre of curvature" ))
376 call check( nf90_put_att(ncid, varid_impp, "units", "Meters" ))
377 call check( nf90_put_att(ncid, varid_impp, "valid_range", "6200 - 6600 km" ))
378 call check( nf90_def_var(ncid, "impact_height@MetaData", nf90_float, nlocs_dimid, varid_imph))
379 call check( nf90_put_att(ncid, varid_imph, "longname", "distance from mean sea level" ))
380 call check( nf90_put_att(ncid, varid_imph, "units", "Meters" ))
381 call check( nf90_put_att(ncid, varid_imph, "valid_range", "0 - 200 km" ))
382 call check( nf90_def_var(ncid, "sensor_azimuth_angle@MetaData", nf90_float, nlocs_dimid, varid_azim))
383 call check( nf90_put_att(ncid, varid_azim, "longname", "GNSS->LEO line of sight" ))
384 call check( nf90_put_att(ncid, varid_azim, "_FillValue", real(missingvalue) ))
385 call check( nf90_put_att(ncid, varid_azim, "units", "Degree" ))
386 call check( nf90_put_att(ncid, varid_azim, "valid_range", "0 - 360 degree" ))
387 call check( nf90_def_var(ncid, "geoid_height_above_reference_ellipsoid@MetaData",nf90_float, nlocs_dimid, varid_geoid))
388 call check( nf90_put_att(ncid, varid_geoid, "longname", "Geoid height above WGS-84 ellipsoid" ))
389 call check( nf90_put_att(ncid, varid_geoid, "units", "Meters" ))
390 call check( nf90_put_att(ncid, varid_geoid, "valid_range", "-200 - 200 m" ))
391 call check( nf90_def_var(ncid, "earth_radius_of_curvature@MetaData",nf90_float, nlocs_dimid, varid_rfict))
392 call check( nf90_put_att(ncid, varid_rfict, "longname", ’"Earths local radius of curvature" ))
393 call check( nf90_put_att(ncid, varid_rfict, "units", "Meters" ))
394 call check( nf90_put_att(ncid, varid_rfict, "valid_range", "6200 - 6600 km" ))
395 call check( nf90_enddef(ncid) )
396 
397 call check( nf90_put_var(ncid, varid_lat, gpsro_data%lat(1:ndata)) )
398 call check( nf90_put_var(ncid, varid_lon, gpsro_data%lon(1:ndata)) )
399 call check( nf90_put_var(ncid, varid_time, gpsro_data%time(1:ndata)) )
400 call check( nf90_put_var(ncid, varid_datetime, gpsro_data%datetime(1:ndata)) )
401 call check( nf90_put_var(ncid, varid_recn, gpsro_data%recn(1:ndata)) )
402 call check( nf90_put_var(ncid, varid_said, gpsro_data%said(1:ndata)) )
403 call check( nf90_put_var(ncid, varid_siid, gpsro_data%siid(1:ndata)) )
404 call check( nf90_put_var(ncid, varid_ptid, gpsro_data%ptid(1:ndata)) )
405 call check( nf90_put_var(ncid, varid_sclf, gpsro_data%sclf(1:ndata)) )
406 call check( nf90_put_var(ncid, varid_asce, gpsro_data%asce(1:ndata)) )
407 call check( nf90_put_var(ncid, varid_ogce, gpsro_data%ogce(1:ndata)) )
408 call check( nf90_put_var(ncid, varid_ref, gpsro_data%ref(1:ndata)) )
409 call check( nf90_put_var(ncid, varid_refoe, gpsro_data%refoe_gsi(1:ndata)) )
410 call check( nf90_put_var(ncid, varid_msl, gpsro_data%msl_alt(1:ndata)) )
411 call check( nf90_put_var(ncid, varid_bnd, gpsro_data%bend_ang(1:ndata)) )
412 call check( nf90_put_var(ncid, varid_bndoe, gpsro_data%bndoe_gsi(1:ndata)) )
413 call check( nf90_put_var(ncid, varid_impp, gpsro_data%impact_para(1:ndata)) )
414 call check( nf90_put_var(ncid, varid_imph, gpsro_data%impact_para(1:ndata)-gpsro_data%rfict(1:ndata)-gpsro_data%geoid(1:ndata) ) )
415 call check( nf90_put_var(ncid, varid_azim, gpsro_data%azim(1:ndata)) )
416 call check( nf90_put_var(ncid, varid_geoid, gpsro_data%geoid(1:ndata)) )
417 call check( nf90_put_var(ncid, varid_rfict, gpsro_data%rfict(1:ndata)) )
418 call check( nf90_close(ncid) )
419 
420 deallocate(gpsro_data%said)
421 deallocate(gpsro_data%siid)
422 deallocate(gpsro_data%sclf)
423 deallocate(gpsro_data%ptid)
424 deallocate(gpsro_data%recn)
425 deallocate(gpsro_data%asce)
426 deallocate(gpsro_data%ogce)
427 deallocate(gpsro_data%time)
428 deallocate(gpsro_data%datetime)
429 deallocate(gpsro_data%lat)
430 deallocate(gpsro_data%lon)
431 deallocate(gpsro_data%rfict)
432 deallocate(gpsro_data%azim)
433 deallocate(gpsro_data%geoid)
434 deallocate(gpsro_data%msl_alt)
435 deallocate(gpsro_data%ref)
436 deallocate(gpsro_data%refoe_gsi)
437 deallocate(gpsro_data%bend_ang)
438 deallocate(gpsro_data%impact_para)
439 deallocate(gpsro_data%bndoe_gsi)
440 
441 contains
442  subroutine check(status)
443  integer, intent ( in) :: status
444 
445  if(status /= nf90_noerr) then
446  print *, trim(nf90_strerror(status))
447  stop "Stopped"
448  end if
449 
450  end subroutine check
451 
452 
453 subroutine refractivity_err_gsi(obsLat, obsZ, GlobalModel, obsErr)
454 real(r_double), intent(in) :: obsLat, obsZ
455 real(r_double), intent(out) :: obsErr
456 logical, intent(in) :: GlobalModel
457 real(r_double) :: obsZ_km
458 
459 obsz_km = obsz / 1000.0
460 
461 if( globalmodel ) then ! for global
462 
463  if( obslat>= 20.0 .or.obslat<= -20.0 ) then
464  obserr=-1.321_r_kind+0.341_r_kind*obsz_km-0.005_r_kind*obsz_km**2
465  else
466  if(obsz_km > 10.0) then
467  obserr=2.013_r_kind-0.060_r_kind*obsz_km+0.0045_r_kind*obsz_km**2
468  else
469  obserr=-1.18_r_kind+0.058_r_kind*obsz_km+0.025_r_kind*obsz_km**2
470  endif
471  endif
472  obserr = 1.0_r_kind/abs(exp(obserr))
473 
474 else ! for regional
475  if( obslat >= 20.0 .or.obslat <= -20.0 ) then
476  if (obsz_km > 10.00) then
477  obserr =-1.321_r_kind+0.341_r_kind*obsz_km-0.005_r_kind*obsz_km**2
478  else
479  obserr =-1.2_r_kind+0.065_r_kind*obsz_km+0.021_r_kind*obsz_km**2
480  endif
481  else
482  if(obsz_km > 10.00) then
483  obserr =2.013_r_kind-0.120_r_kind*obsz_km+0.0065_r_kind*obsz_km**2
484  else
485  obserr =-1.19_r_kind+0.03_r_kind*obsz_km+0.023_r_kind*obsz_km**2
486  endif
487  endif
488  obserr = 1.0_r_kind/abs(exp(obserr))
489 endif
490 
491 end subroutine refractivity_err_gsi
492 
493 subroutine bendingangle_err_gsi(obsLat, obsZ, obsErr)
494 real(r_double), intent(in) :: obsLat, obsZ
495 real(r_double), intent(out) :: obsErr
496 real(r_double) :: obsZ_km
497 
498 obsz_km = obsz / 1000.0
499 if((said==41).or.(said==722).or.(said==723).or.(said==4).or.(said==42).or.&
500  (said==3).or.(said==821.or.(said==421)).or.(said==440).or.(said==43)) then
501  if( abs(obslat)>= 40.00 ) then
502 
503  if(obsz_km>12.) then
504  obserr=0.19032_r_kind+0.287535_r_kind*obsz_km-0.00260813_r_kind*obsz_km**2
505  else
506  obserr=-3.20978_r_kind+1.26964_r_kind*obsz_km-0.0622538_r_kind*obsz_km**2
507  endif
508  else
509  if(obsz_km>18.) then
510  obserr=-1.87788_r_kind+0.354718_r_kind*obsz_km-0.00313189_r_kind*obsz_km**2
511  else
512  obserr=-2.41024_r_kind+0.806594_r_kind*obsz_km-0.027257_r_kind*obsz_km**2
513  endif
514  endif
515  obserr = 0.001_r_kind/abs(exp(obserr))
516 
517  else !!!! CDAAC processing
518  if( abs(obslat)>= 40.00 ) then
519  if (obsz_km > 12.00) then
520  obserr=-0.685627_r_kind+0.377174_r_kind*obsz_km-0.00421934_r_kind*obsz_km**2
521  else
522  obserr=-3.27737_r_kind+1.20003_r_kind*obsz_km-0.0558024_r_kind*obsz_km**2
523  endif
524  else
525  if(obsz_km >18.00) then
526  obserr=-2.73867_r_kind+0.447663_r_kind*obsz_km-0.00475603_r_kind*obsz_km**2
527  else
528  obserr=-3.45303_r_kind+0.908216_r_kind*obsz_km-0.0293331_r_kind*obsz_km**2
529  endif
530  endif
531  obserr = 0.001_r_kind/abs(exp(obserr))
532 
533 endif
534 
535 end subroutine bendingangle_err_gsi
536 !!!!!!________________________________________________________
537 
538 !!!!!! SUBROUTINE W3FS21 was copied from GSI/src/libs/w3nco_v2.0.6/w3fs21.f and iw3jdn.f
539 SUBROUTINE w3fs21(IDATE, NMIN)
540  INTEGER IDATE(5)
541  INTEGER NMIN
542  INTEGER IYEAR, NDAYS, IJDN
543  INTEGER JDN78
544  DATA jdn78 / 2443510 /
545  nmin = 0
546 
547  iyear = idate(1)
548 
549  IF (iyear.LE.99) THEN
550  IF (iyear.LT.78) THEN
551  iyear = iyear + 2000
552  ELSE
553  iyear = iyear + 1900
554  ENDIF
555  ENDIF
556 
557 ! COMPUTE JULIAN DAY NUMBER FROM YEAR, MONTH, DAY
558  ijdn = idate(3) - 32075 &
559  + 1461 * (iyear + 4800 + (idate(2) - 14) / 12) / 4 &
560  + 367 * (idate(2)- 2 - (idate(2) -14) / 12 * 12) / 12 &
561  - 3 * ((iyear + 4900 + (idate(2) - 14) / 12) / 100) / 4
562 
563 ! SUBTRACT JULIAN DAY NUMBER OF JAN 1,1978 TO GET THE
564 ! NUMBER OF DAYS BETWEEN DATES
565  ndays = ijdn - jdn78
566  nmin = ndays * 1440 + idate(4) * 60 + idate(5)
567  RETURN
568 END SUBROUTINE w3fs21
569 
570 end program
bendingangle_err_gsi
subroutine bendingangle_err_gsi(obsLat, obsZ, obsErr)
Definition: gnssro_bufr2ioda.f90:494
refractivity_err_gsi
subroutine refractivity_err_gsi(obsLat, obsZ, GlobalModel, obsErr)
Definition: gnssro_bufr2ioda.f90:454
check
subroutine check(status)
Definition: gnssro_bufr2ioda.f90:443
w3fs21
subroutine w3fs21(IDATE, NMIN)
Definition: gnssro_bufr2ioda.f90:540
gnssro_bufr2ioda
Definition: gnssro_mod.f90:8
gpsro_type
Definition: gnssro_bufr2ioda.f90:57