FV3 Bundle
gsidiag_rad_bin2nc4.f90
Go to the documentation of this file.
1 program convert_rad_diag
2 
3  use read_diag,only:read_radiag_header, diag_header_fix_list, diag_header_chan_list, diag_data_name_list
4  use read_diag,only:read_radiag_data, diag_data_fix_list, diag_data_extra_list, diag_data_chan_list
5  use ncd_kinds,only: r_quad, r_single
6 
7  use nc_diag_write_mod, only: nc_diag_init, nc_diag_header, nc_diag_metadata, &
8  nc_diag_write, nc_diag_data2d, nc_diag_chaninfo_dim_set, nc_diag_chaninfo
9 
10  implicit none
11 
12  real,parameter:: missing = -9.99e9
13  integer,parameter:: imissing = -999999
14 
15  integer nargs, iargc, n
16  character*256, allocatable :: arg(:)
17 
18  type(diag_header_fix_list ) :: header_fix
19  type(diag_header_chan_list),allocatable :: header_chan(:)
20  type(diag_data_name_list) :: headname
21 
22  type(diag_data_fix_list) :: data_fix
23  type(diag_data_chan_list) ,allocatable :: data_chan(:)
24  type(diag_data_extra_list) ,allocatable :: dataextra(:,:)
25 
26 
27  integer i
28 !!! real(r_quad) :: ret_var
29 !!! real(r_quad) :: ret_stddev
30 
31 ! commandline variables
32  logical :: debug
33  integer :: npred_read
34  logical :: sst_ret
35  integer :: iversion
36  logical :: append_suffix
37 
38  character*256 infn, outfn
39  logical linfile, loutfile
40 
41  integer,parameter :: inlun = 51
42  integer,parameter :: outlun= 52
43  integer,parameter :: nllun = 53
44 
45  integer strlen, iflag
46  integer iuse, ich, nch, ipr
47 
48  logical,dimension(:),allocatable :: luse
49  logical lqcpass
50 
51  integer,parameter :: nvar = 4 ! number of positions in array needed for inline variance calc
52 ! variables for output, all to be allocated as nchan. Variances will be calculated using inline algorithm
53 ! accredited to Welford, according
54 !!! real(r_quad),dimension(:),allocatable :: nobstotal, nobsassim, tbtotal, tbassim, omf_nbc , omf_bc , sigo, jo
55 !!! real(r_quad),dimension(:,:),allocatable :: vomf_nbc, vomf_bc
56 ! total bias and fixed bias terms. When Yanqui's variational angle correction is brought in, this may need to be updated.
57 !!! real(r_quad),dimension(:),allocatable :: totbias , fixbias
58 !!! real(r_quad),dimension(:,:),allocatable :: vtotbias, vfixbias
59 ! variational bias correction variables, which will be allocated as nchan and npred_read
60 !!! real(r_quad),dimension(:,:),allocatable :: biasterms
61 !!! real(r_quad),dimension(:,:,:),allocatable :: vbiasterms
62 ! Definitions for above variables -
63 ! nobstotal - number of observations considered - total
64 ! nobsassim - number of observations that are assimilated
65 ! tbtotal - mean brightness temperature of all observations
66 ! tbassim - mean brightness temperature of assimilated observations
67 ! omf_nbc/vomf_nbc - mean/variance of O-F without bias correction applied
68 ! omf_bc/vomf_bc - mean/variance of O-F with bias correction applied
69 ! sigo - mean observation error of assimilated observations
70 ! jo - mean cost (Jo) of assimilated observations
71 ! totbias/vtotbias - mean/variance of bias correction applied to assimilated observations
72 ! fixbias/vfixbias - mean/variance of fixed scan angle position bias correction (sac.x-derived)
73 ! biasterms
74 ! /vbiasterms - means/variances of the variational terms as defined by npred_read.
75 
76 ! single variables used later for printing purposes
77  integer :: inobstotal, inobsassim
78  real(r_single) :: rvomf_nbc, rvomf_bc, rvtotbias, rvfixbias
79  real(r_single),dimension(:),allocatable :: rvbiasterms
80  character(len=13),dimension(:),allocatable :: chfrwn
81 
82  integer,parameter :: max_npred = 9
83 
84  character(len=10),dimension(max_npred) :: biasnames
85  data biasnames / 'bc_const ', &
86  'bc_satang ', &
87  'bc_tlap ', &
88  'bc_tlap2 ', &
89  'bc_clw ', &
90  'bc_coslat ', &
91  'bc_sinlat ', &
92  'bc_emis ', &
93  'bc_sst ' /
94 
95  real(r_quad) :: cvar, rch
96 
97 
98  nargs = iargc()
99  if( nargs.eq.0 ) then
100  call usage
101  else
102  debug = .false.
103  npred_read = 7
104  sst_ret = .false.
105  iversion = -9999
106  append_suffix = .false.
107 
108  allocate(arg(nargs))
109  do n=1,nargs
110  call getarg(n,arg(n))
111  enddo
112  do n=1,nargs
113  if (trim(arg(n)).eq.'-debug' ) debug=.true.
114  if (trim(arg(n)).eq.'-sst_ret' ) sst_ret=.true.
115  if (trim(arg(n)).eq.'-append_nc4') append_suffix=.true.
116  if (trim(arg(n)).eq.'-npred' ) read ( arg(n+1),* ) npred_read
117  if (trim(arg(n)).eq.'-iversion' ) read ( arg(n+1),* ) iversion
118  enddo
119  endif
120 
121 
122  if (debug) write(*,*)'Debugging on - Verbose Printing'
123 
124  ! get infn from command line
125  call getarg(nargs, infn)
126 
127  strlen = len(trim(infn))
128 
129  write(*,*)'Input bin diag: ',trim(infn)
130  inquire(file=trim(infn), exist=linfile)
131  if (.not. linfile) then
132  write(*,*)trim(infn) // ' does not exist - exiting'
133  call abort
134  endif
135 
136  if (.not. append_suffix) then
137  outfn = infn(1:strlen-3) // 'nc4' ! assumes GMAO diag filename format ending with .bin, and replaces it
138  else
139  outfn = infn(1:strlen) // '.nc4' ! if not GMAO format, use append_suffix = .true. in namelist
140  ! to simply append infile with .nc4 suffix
141  endif
142 
143  write(*,*)'Output NC4 diag: ',trim(outfn)
144  inquire(file=trim(outfn), exist=loutfile)
145  if (loutfile) write(*,*)'WARNING: ' // trim(infn) // ' exists - overwriting'
146 
147  iflag = 0
148 
149  open(inlun,file=infn,form='unformatted',convert='big_endian')
150  call nc_diag_init(outfn)
151 
152  call read_radiag_header( inlun, npred_read, sst_ret, header_fix, header_chan, headname, iflag, debug )
153 
154  call nc_diag_chaninfo_dim_set(header_fix%nchan)
155 
156  call nc_diag_header("Satellite_Sensor", header_fix%isis )
157  call nc_diag_header("Satellite", header_fix%id ) ! sat type
158  call nc_diag_header("Observation_type", header_fix%obstype ) ! observation type
159  call nc_diag_header("Outer_Loop_Iteration", header_fix%jiter )
160  call nc_diag_header("Number_of_channels", header_fix%nchan ) ! number of channels in the sensor
161  call nc_diag_header("Number_of_Predictors", header_fix%npred ) ! number of updating bias correction predictors
162  call nc_diag_header("date_time", header_fix%idate ) ! time (yyyymmddhh)
163  call nc_diag_header("ireal_radiag", header_fix%ireal )
164  call nc_diag_header("ipchan_radiag", header_fix%ipchan )
165  call nc_diag_header("iextra", header_fix%iextra )
166  call nc_diag_header("jextra", header_fix%jextra )
167  call nc_diag_header("idiag", header_fix%idiag )
168  call nc_diag_header("angord", header_fix%angord )
169  call nc_diag_header("iversion_radiag", header_fix%iversion)
170  call nc_diag_header("New_pc4pred", header_fix%inewpc ) ! indicator of newpc4pred (1 on, 0 off)
171  call nc_diag_header("ioff0", header_fix%isens ) ! i think ioff0 = isens
172 
173 
174  nch = header_fix%nchan
175 
176  allocate(luse(nch))
177 
178  if (debug) then
179  write(*,*)'Number of Channels: ',nch
180  write(*,*)'Number of variationalbc predictors: ',npred_read
181  write(*,*)' predictors: ',biasnames(1:npred_read)
182  write(*,*)' iversion=',header_fix%iversion
183  endif
184 
185  if (iversion .gt. 0) then
186  write(*,*)'BE AWARE THAT iversion IS BEING OVERRIDEN!'
187  write(*,*)' iversion diag, override=',header_fix%iversion,iversion
188  write(*,*)' (this was made necessary w/ emis bc...hopefully only temporary)'
189  header_fix%iversion = iversion
190  endif
191 
192  do i=1,nch
193  call nc_diag_chaninfo("chaninfoidx", i )
194  call nc_diag_chaninfo("frequency", header_chan(i)%freq )
195  call nc_diag_chaninfo("polarization", header_chan(i)%polar )
196  call nc_diag_chaninfo("wavenumber", header_chan(i)%wave )
197  call nc_diag_chaninfo("error_variance", header_chan(i)%varch )
198  call nc_diag_chaninfo("mean_lapse_rate", header_chan(i)%tlapmean)
199  call nc_diag_chaninfo("use_flag", header_chan(i)%iuse )
200  call nc_diag_chaninfo("sensor_chan", header_chan(i)%nuchan )
201  call nc_diag_chaninfo("satinfo_chan", header_chan(i)%iochan )
202  end do
203 
204 
205  do while (iflag .ge. 0) ! iflag == 0 means the end of the file
206  call read_radiag_data ( inlun, header_fix, .false., data_fix, data_chan, &
207  dataextra, iflag )
208 
209  if (iflag .lt. 0) cycle
210 
211  do ich=1,nch
212  lqcpass = luse(ich) .and. data_chan(ich)%qcmark .eq. 0
213 
214  call nc_diag_metadata("Channel_Index", i )
215  call nc_diag_metadata("Observation_Class", ' rad' )
216  call nc_diag_metadata("Latitude", data_fix%lat ) ! observation latitude (degrees)
217  call nc_diag_metadata("Longitude", data_fix%lon ) ! observation longitude (degrees)
218 
219  call nc_diag_metadata("Elevation", data_fix%zsges ) ! model (guess) elevation at observation location
220 
221  call nc_diag_metadata("Obs_Time", data_fix%obstime ) ! observation time (hours relative to analysis time)
222 
223  call nc_diag_metadata("Scan_Position", data_fix%senscn_pos ) ! sensor scan position
224  call nc_diag_metadata("Sat_Zenith_Angle", data_fix%satzen_ang ) ! satellite zenith angle (degrees)
225  call nc_diag_metadata("Sat_Azimuth_Angle", data_fix%satazm_ang ) ! satellite azimuth angle (degrees)
226  call nc_diag_metadata("Sol_Zenith_Angle", data_fix%solzen_ang ) ! solar zenith angle (degrees)
227  call nc_diag_metadata("Sol_Azimuth_Angle", data_fix%solazm_ang ) ! solar azimuth angle (degrees)
228  call nc_diag_metadata("Sun_Glint_Angle", data_fix%sungln_ang ) ! sun glint angle (degrees) (sgagl)
229 
230  call nc_diag_metadata("Water_Fraction", data_fix%water_frac ) ! fractional coverage by water
231  call nc_diag_metadata("Land_Fraction", data_fix%land_frac ) ! fractional coverage by land
232  call nc_diag_metadata("Ice_Fraction", data_fix%ice_frac ) ! fractional coverage by ice
233  call nc_diag_metadata("Snow_Fraction", data_fix%snow_frac ) ! fractional coverage by snow
234 
235  if (.not. sst_ret) then
236  call nc_diag_metadata("Water_Temperature", data_fix%water_temp ) ! surface temperature over water (K)
237  call nc_diag_metadata("Land_Temperature", data_fix%land_temp ) ! surface temperature over land (K)
238  call nc_diag_metadata("Ice_Temperature", data_fix%ice_temp ) ! surface temperature over ice (K)
239  call nc_diag_metadata("Snow_Temperature", data_fix%snow_temp ) ! surface temperature over snow (K)
240  call nc_diag_metadata("Soil_Temperature", data_fix%soil_temp ) ! soil temperature (K)
241  call nc_diag_metadata("Soil_Moisture", data_fix%soil_mois ) ! soil moisture
242  call nc_diag_metadata("Land_Type_Index", data_fix%land_type ) ! surface land type
243 
244  call nc_diag_metadata("tsavg5", missing ) ! SST first guess used for SST retrieval
245  call nc_diag_metadata("sstcu", missing ) ! NCEP SST analysis at t
246  call nc_diag_metadata("sstph", missing ) ! Physical SST retrieval
247  call nc_diag_metadata("sstnv", missing ) ! Navy SST retrieval
248  call nc_diag_metadata("dta", missing ) ! d(ta) corresponding to sstph
249  call nc_diag_metadata("dqa", missing ) ! d(qa) corresponding to sstph
250  call nc_diag_metadata("dtp_avh", missing ) ! data type
251  else
252  call nc_diag_metadata("Water_Temperature", missing ) ! surface temperature over water (K)
253  call nc_diag_metadata("Land_Temperature", missing ) ! surface temperature over land (K)
254  call nc_diag_metadata("Ice_Temperature", missing ) ! surface temperature over ice (K)
255  call nc_diag_metadata("Snow_Temperature", missing ) ! surface temperature over snow (K)
256  call nc_diag_metadata("Soil_Temperature", missing ) ! soil temperature (K)
257  call nc_diag_metadata("Soil_Moisture", missing ) ! soil moisture
258  call nc_diag_metadata("Land_Type_Index", missing ) ! surface land type
259 
260  call nc_diag_metadata("tsavg5", data_fix%water_temp ) ! SST first guess used for SST retrieval
261  call nc_diag_metadata("sstcu", data_fix%land_temp ) ! NCEP SST analysis at t
262  call nc_diag_metadata("sstph", data_fix%ice_temp ) ! Physical SST retrieval
263  call nc_diag_metadata("sstnv", data_fix%snow_temp ) ! Navy SST retrieval
264  call nc_diag_metadata("dta", data_fix%soil_temp ) ! d(ta) corresponding to sstph
265  call nc_diag_metadata("dqa", data_fix%soil_mois ) ! d(qa) corresponding to sstph
266  call nc_diag_metadata("dtp_avh", data_fix%land_type ) ! data type
267  endif
268 
269  call nc_diag_metadata("Vegetation_Fraction", data_fix%veg_frac )
270  call nc_diag_metadata("Snow_Depth", data_fix%snow_depth )
271 ! qcdiag1 = slot 25 ; qcdiag2 = slot 26 - simply mapping. not attempting to add logic for missing vals
272  call nc_diag_metadata("tpwc_amsua", missing )
273  call nc_diag_metadata("clw_guess_retrieval", data_fix%qcdiag1 )
274 
275  call nc_diag_metadata("Sfc_Wind_Speed", data_fix%sfc_wndspd )
276  call nc_diag_metadata("Cloud_Frac", data_fix%qcdiag1 )
277  call nc_diag_metadata("CTP", data_fix%qcdiag2 )
278  call nc_diag_metadata("CLW", data_fix%qcdiag1 )
279  call nc_diag_metadata("TPWC", data_fix%qcdiag2 )
280  call nc_diag_metadata("clw_obs", data_fix%qcdiag1 )
281  call nc_diag_metadata("clw_guess", data_fix%qcdiag2 )
282 
283  call nc_diag_metadata("Foundation_Temperature", data_fix%tref ) ! reference temperature (Tr) in NSST
284  call nc_diag_metadata("SST_Warm_layer_dt", data_fix%dtw ) ! dt_warm at zob
285  call nc_diag_metadata("SST_Cool_layer_tdrop", data_fix%dtc ) ! dt_cool at zob
286  call nc_diag_metadata("SST_dTz_dTfound", data_fix%tz_tr ) ! d(Tz)/d(Tr)
287 
288  call nc_diag_metadata("Observation", data_chan(ich)%tbobs ) ! observed brightness temperature (K)
289  call nc_diag_metadata("Obs_Minus_Forecast_adjusted", data_chan(ich)%omgbc ) ! observed - simulated Tb with bias corrrection (K)
290  call nc_diag_metadata("Obs_Minus_Forecast_unadjusted", data_chan(ich)%omgnbc ) ! observed - simulated Tb with no bias correction (K)
291 ! errinv = sqrt(varinv(ich_diag(i)))
292  call nc_diag_metadata("Inverse_Observation_Error", data_chan(ich)%errinv )
293 
294 ! useflag=one
295 ! if (iuse_rad(ich(ich_diag(i))) < 1) useflag=-one
296 
297  call nc_diag_metadata("QC_Flag", data_chan(ich)%qcmark ) ! quality control mark or event indicator
298 
299  call nc_diag_metadata("Emissivity", data_chan(ich)%emiss ) ! surface emissivity
300  call nc_diag_metadata("Weighted_Lapse_Rate", data_chan(ich)%tlap ) ! stability index
301  call nc_diag_metadata("dTb_dTs", data_chan(ich)%tb_tz ) ! d(Tb)/d(Ts)
302 
303  call nc_diag_metadata("BC_Constant", data_chan(ich)%bicons ) ! constant bias correction term
304  call nc_diag_metadata("BC_Scan_Angle", data_chan(ich)%biang ) ! scan angle bias correction term
305  call nc_diag_metadata("BC_Cloud_Liquid_Water", data_chan(ich)%biclw ) ! CLW bias correction term
306  call nc_diag_metadata("BC_Lapse_Rate_Squared", data_chan(ich)%bilap2 ) ! square lapse rate bias correction term
307  call nc_diag_metadata("BC_Lapse_Rate", data_chan(ich)%bilap ) ! lapse rate bias correction term
308  call nc_diag_metadata("BC_Cosine_Latitude_times_Node", data_chan(ich)%bicos ) ! node*cos(lat) bias correction term
309  call nc_diag_metadata("BC_Sine_Latitude", data_chan(ich)%bisin ) ! sin(lat) bias correction term
310  call nc_diag_metadata("BC_Emissivity", data_chan(ich)%biemis ) ! emissivity sensitivity bias correction term
311  if (header_fix%angord .eq. 1) then
312  call nc_diag_metadata("BC_Fixed_Scan_Position", data_chan(ich)%bifix(1) )
313  else if (header_fix%angord .ge. 2) then
314  call nc_diag_data2d('BC_angord ', data_chan(ich)%bifix )
315  endif
316 !
317  enddo
318 
319  enddo
320 
321 ! finalize NCDIAG
322  call nc_diag_write
323 end program convert_rad_diag
324 
325 subroutine usage
326  write(6,100)
327 100 format( "Usage: ",/,/ &
328  " convert_rad_diag.x <options> <filename>",/,/ &
329  "where options:",/ &
330  " -debug : Set debug verbosity",/ &
331  " -sst_ret : SST BC term is included (default: not included)",/ &
332  " -npred INT : Number of preductors (default: 7)",/ &
333  " -iversion INT : Override iversion with INT (default: use internal iversion)",/ &
334  " -append_txt : Append .txt suffix, instead of replace last three",/ &
335  " characters (default: replaced)",/ &
336  " Note: The GMAO diag files end with .bin or .nc4,",/ &
337  " which is where fixed 3-char truncation originates",/,/,/ &
338  " Example:",/ &
339  " convert_rad_diag.x nc_4emily.diag_hirs4_n19_ges.20161202_00z.bin",/ &
340  " Output file:",/ &
341  " nc_4emily.diag_hirs4_n19_ges.20161202_00z.nc4",/ &
342  )
343  stop
344 
345 end subroutine usage
346 
read_diag::diag_header_fix_list
Definition: read_diag.f90:81
read_diag::read_radiag_header
subroutine, public read_radiag_header(ftin, npred_radiag, retrieval, header_fix, header_chan, data_name, iflag, lverbose)
Definition: read_diag.f90:321
crtm_parameters::strlen
integer, parameter, public strlen
Definition: CRTM_Parameters.f90:116
nc_diag_write_mod::nc_diag_init
subroutine nc_diag_init(filename, append)
Definition: nc_diag_write_mod.F90:292
read_diag::diag_data_fix_list
Definition: read_diag.f90:116
read_diag::read_radiag_data
subroutine, public read_radiag_data(ftin, header_fix, retrieval, data_fix, data_chan, data_extra, iflag)
Definition: read_diag.f90:729
usage
subroutine usage
Definition: gsidiag_aod_bin2nc4.f90:198
read_diag::diag_header_chan_list
Definition: read_diag.f90:105
nc_diag_write_mod
Definition: nc_diag_write_mod.F90:19
ncd_kinds::r_quad
integer, parameter, public r_quad
Definition: ncd_kinds.F90:82
read_diag::diag_data_chan_list
Definition: read_diag.f90:150
read_diag
Definition: read_diag.f90:36
convert_rad_diag
program convert_rad_diag
Definition: gsidiag_rad_bin2nc4.f90:1
ncd_kinds
Definition: ncd_kinds.F90:1
read_diag::diag_data_name_list
Definition: read_diag.f90:100
nc_diag_write_mod::nc_diag_write
subroutine nc_diag_write
Definition: nc_diag_write_mod.F90:505
ncd_kinds::r_single
integer, parameter, public r_single
Definition: ncd_kinds.F90:79
read_diag::diag_data_extra_list
Definition: read_diag.f90:171