Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fix bug for thinning option in the goes_abi converter and update README files #5

Open
wants to merge 1 commit into
base: main
Choose a base branch
from
+129 −125
Open

Fix bug for thinning option in the goes_abi converter and update README files #5

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
6 changes: 4 additions & 2 deletions goes_abi/README.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,7 @@ Compile: cd goes_abi/src; make
```
Usage: goes_abi_converter.x
```
Run-time options are specified through ``namelist.goes_abi_converter``.
Run-time options are specified through ``namelist.goes_abi_converter``.
List of GoesReBroadcast files (exclude path) to be processed are specified in a plain text file.

Example namelist.goes_abi_converter:
Expand All @@ -16,7 +16,9 @@ Example namelist.goes_abi_converter:
data_dir = '/data/goes', ! path of the GRB nc files
data_id = 'OR_ABI-L1b-RadF-M6' ! prefix of the downloaded GRB nc files
sat_id = 'G16'
n_subsample = 1
n_subsample = 1 ! skip n_subsample of pixels
do_superob = .false. ! default option is no superobbing (if do_superob = .true., then no thinning is done)
superob_halfwidth = 29
/
```
Example content of flist.txt:
Expand Down
132 changes: 65 additions & 67 deletions goes_abi/src/goes_abi_converter.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -81,10 +81,9 @@ program Goes_ReBroadcast_converter
integer(i_kind) :: n_subsample
logical :: do_superob
integer(i_kind) :: superob_halfwidth
logical :: do_thinning
logical :: write_iodav1

namelist /data_nml/ nc_list_file, data_dir, data_id, sat_id, do_thinning, n_subsample, do_superob, superob_halfwidth
namelist /data_nml/ nc_list_file, data_dir, data_id, sat_id, n_subsample, do_superob, superob_halfwidth

real(r_kind) :: sdtb ! to be done
integer(i_kind) :: istat
Expand Down Expand Up @@ -117,7 +116,6 @@ program Goes_ReBroadcast_converter
data_dir = '.'
data_id = 'OR_ABI-L1b-RadC-M3'
sat_id = 'G16'
do_thinning = .false.
n_subsample = 1
do_superob = .false.
superob_halfwidth = 1
Expand Down Expand Up @@ -842,70 +840,6 @@ subroutine output_iodav1(fname, time_start, nx, ny, nband, got_latlon, lat, lon,

nvars = nband

if ( do_thinning ) then
nlocs = 0
do iline = 1, ny, n_subsample
do isample = 1, nx, n_subsample
if ( .not. got_latlon(isample,iline) ) cycle
if ( sat_zen(isample,iline) > 80.0 ) cycle
! qf (DQF, Data Quality Flag)
! 0:good, 1:conditionally_usable, 2:out_of_range, 3:no_value
! keep only qf=0,1 pixels
if ( all(qf(:,isample,iline) > 1) ) cycle
if ( all(bt(:,isample,iline)<0.0) ) cycle
nlocs = nlocs + 1
end do
end do

write(0,*) 'nlocs = ', nlocs
if ( nlocs <= 0 ) then
return
end if

allocate (name_var_tb(1:nband))
allocate (datetime(nlocs))
allocate (lat_out(nlocs))
allocate (lon_out(nlocs))
allocate (scan_pos_out(nlocs))
allocate (sat_zen_out(nlocs))
allocate (sat_azi_out(nlocs))
allocate (sun_zen_out(nlocs))
allocate (sun_azi_out(nlocs))
allocate (bt_out(nband,nlocs))
allocate (err_out(nband,nlocs))
allocate (qf_out(nband,nlocs))

read(time_start( 1: 4), '(i4)') iyear
read(time_start( 6: 7), '(i2)') imonth
read(time_start( 9:10), '(i2)') iday
read(time_start(12:13), '(i2)') ihour
read(time_start(15:16), '(i2)') imin
read(time_start(18:19), '(i2)') isec

iloc = 0
do iline = 1, ny, n_subsample
do isample = 1, nx, n_subsample
if ( .not. got_latlon(isample,iline) ) cycle
if ( sat_zen(isample,iline) > 80.0 ) cycle
if ( all(qf(:,isample,iline) > 1) ) cycle
if ( all(bt(:,isample,iline)<0.0) ) cycle
iloc = iloc + 1
write(unit=datetime(iloc), fmt='(i4,a,i2.2,a,i2.2,a,i2.2,a,i2.2,a,i2.2,a)') &
iyear, '-', imonth, '-', iday, 'T', ihour, ':', imin, ':', isec, 'Z'
lat_out(iloc) = lat(isample,iline)
lon_out(iloc) = lon(isample,iline)
sat_zen_out(iloc) = sat_zen(isample,iline)
sun_zen_out(iloc) = sun_zen(isample,iline)
bt_out(1:nband,iloc) = bt(1:nband,isample,iline)
qf_out(1:nband,iloc) = qf(1:nband,isample,iline)
scan_pos_out(iloc) = isample
sat_azi_out(iloc) = missing_r
sun_azi_out(iloc) = missing_r
err_out(1:nband,iloc) = 1.0 !missing_r
end do
end do
end if

! Superobbing code is based on WRFDA's AHI superobbing code from https://github.com/wrf-model/WRF/blob/develop/var/da/da_radiance/da_read_obs_netcdf4ahi_jaxa.inc
if ( do_superob ) then
nlocs = 0
Expand Down Expand Up @@ -1025,6 +959,70 @@ subroutine output_iodav1(fname, time_start, nx, ny, nband, got_latlon, lat, lon,

end do fov_loop
end do scan_loop

else

nlocs = 0
do iline = 1, ny, n_subsample
do isample = 1, nx, n_subsample
if ( .not. got_latlon(isample,iline) ) cycle
if ( sat_zen(isample,iline) > 80.0 ) cycle
! qf (DQF, Data Quality Flag)
! 0:good, 1:conditionally_usable, 2:out_of_range, 3:no_value
! keep only qf=0,1 pixels
if ( all(qf(:,isample,iline) > 1) ) cycle
if ( all(bt(:,isample,iline)<0.0) ) cycle
nlocs = nlocs + 1
end do
end do

write(0,*) 'nlocs = ', nlocs
if ( nlocs <= 0 ) then
return
end if

allocate (name_var_tb(1:nband))
allocate (datetime(nlocs))
allocate (lat_out(nlocs))
allocate (lon_out(nlocs))
allocate (scan_pos_out(nlocs))
allocate (sat_zen_out(nlocs))
allocate (sat_azi_out(nlocs))
allocate (sun_zen_out(nlocs))
allocate (sun_azi_out(nlocs))
allocate (bt_out(nband,nlocs))
allocate (err_out(nband,nlocs))
allocate (qf_out(nband,nlocs))

read(time_start( 1: 4), '(i4)') iyear
read(time_start( 6: 7), '(i2)') imonth
read(time_start( 9:10), '(i2)') iday
read(time_start(12:13), '(i2)') ihour
read(time_start(15:16), '(i2)') imin
read(time_start(18:19), '(i2)') isec

iloc = 0
do iline = 1, ny, n_subsample
do isample = 1, nx, n_subsample
if ( .not. got_latlon(isample,iline) ) cycle
if ( sat_zen(isample,iline) > 80.0 ) cycle
if ( all(qf(:,isample,iline) > 1) ) cycle
if ( all(bt(:,isample,iline)<0.0) ) cycle
iloc = iloc + 1
write(unit=datetime(iloc), fmt='(i4,a,i2.2,a,i2.2,a,i2.2,a,i2.2,a,i2.2,a)') &
iyear, '-', imonth, '-', iday, 'T', ihour, ':', imin, ':', isec, 'Z'
lat_out(iloc) = lat(isample,iline)
lon_out(iloc) = lon(isample,iline)
sat_zen_out(iloc) = sat_zen(isample,iline)
sun_zen_out(iloc) = sun_zen(isample,iline)
bt_out(1:nband,iloc) = bt(1:nband,isample,iline)
qf_out(1:nband,iloc) = qf(1:nband,isample,iline)
scan_pos_out(iloc) = isample
sat_azi_out(iloc) = missing_r
sun_azi_out(iloc) = missing_r
err_out(1:nband,iloc) = 1.0 !missing_r
end do
end do
end if

call open_netcdf_for_write(trim(fname),ncfileid)
Expand Down
113 changes: 57 additions & 56 deletions obs2ioda-v2/README.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,101 +3,102 @@
```
Compile: cd obs2ioda-v2/src; make
```
NCEP BUFR library (https://github.com/NOAA-EMC/NCEPLIBS-bufr) is required to compile ``obs2ioda-v2.x``.
NCEP BUFR library (https://github.com/NOAA-EMC/NCEPLIBS-bufr) is required to compile ``obs2ioda-v2.x``.
Edit obs2ioda-v2/src/Makefile to set proper BUFR_LIB before ``make``.

## caveate
NetCDF-Fortran interface does not allow reading/writing NF90_STRING, so ``station_id`` and ``variable_names`` are still written out as
``char station_id(nlocs, nstring)``
``char variable_names(nvars, nstring)``
rather than
``string station_id(nlocs)``
NetCDF-Fortran interface does not allow reading/writing NF90_STRING, so ``station_id`` and ``variable_names`` are still written out as
``char station_id(nlocs, nstring)``
``char variable_names(nvars, nstring)``
rather than
``string station_id(nlocs)``
``string variable_names(nvars)``

## Converting PREPBUFR and BUFR files
```
Usage: obs2ioda-v2.x [-i input_dir] [-o output_dir] [bufr_filename(s)_to_convert] [-split]
```
If [-i input_dir] [-o output_dir] are not specified in the command line, the default is the current working directory.
If [bufr_filename(s)_to_convert] is not specified in the command line, the code looks for file name, **prepbufr.bufr** (also **satwnd.bufr**, **gnssro.bufr**, **amsua.bufr**, **airs.bufr**, **mhs.bufr**, **iasi.bufr**, **cris.bufr**), in the input/working directory. If the file exists, do the conversion, otherwise skip it.
If [-i input_dir] [-o output_dir] are not specified in the command line, the default is the current working directory.
If [bufr_filename(s)_to_convert] is not specified in the command line, the code looks for file name, **prepbufr.bufr** (also **satwnd.bufr**, **gnssro.bufr**, **amsua.bufr**, **airs.bufr**, **mhs.bufr**, **iasi.bufr**, **cris.bufr**), in the input/working directory. If the file exists, do the conversion, otherwise skip it.
If specify ``-split``, the converted file will contain hourly data.

> obs2ioda-v2.x -i input_dir -o output_dir prepbufr.gdas.YYYYMMDD.tHHz.nr

Example output files (date in the output filename is extracted from the input bufr files):
aircraft_obs_YYYYMMDDHH.h5
ascat_obs_YYYYMMDDHH.h5
profiler_obs_YYYYMMDDHH.h5
satwind_obs_YYYYMMDDHH.h5
sfc_obs_YYYYMMDDHH.h5
sondes_obs_YYYYMMDDHH.h5
Example output files (date in the output filename is extracted from the input bufr files):
aircraft_obs_YYYYMMDDHH.h5
ascat_obs_YYYYMMDDHH.h5
profiler_obs_YYYYMMDDHH.h5
satwind_obs_YYYYMMDDHH.h5
sfc_obs_YYYYMMDDHH.h5
sondes_obs_YYYYMMDDHH.h5

> obs2ioda-v2.x -i input_dir -o output_dir gdas.satwnd.tHHz.YYYYMMDD.bufr

Example output files (date in the output filename is extracted from the input bufr files):
satwnd_obs_YYYYMMDDHH.h5 (GOES-16/GOES-17, AVHRR (METOP/NOAA), VIIRS (NPP/NOAA), LEOGEO AMVs)
Example output files (date in the output filename is extracted from the input bufr files):
satwnd_obs_YYYYMMDDHH.h5 (GOES-16/GOES-17, AVHRR (METOP/NOAA), VIIRS (NPP/NOAA), LEOGEO AMVs)

> obs2ioda-v2.x -i input_dir -o output_dir gdas.1bamua.tHHz.YYYYMMDD.bufr

Example output files:
amsua_metop-a_obs_YYYYMMDDHH.h5
amsua_metop-b_obs_YYYYMMDDHH.h5
amsua_n15_obs_YYYYMMDDHH.h5
amsua_n18_obs_YYYYMMDDHH.h5
amsua_n19_obs_YYYYMMDDHH.h5
Example output files:
amsua_metop-a_obs_YYYYMMDDHH.h5
amsua_metop-b_obs_YYYYMMDDHH.h5
amsua_n15_obs_YYYYMMDDHH.h5
amsua_n18_obs_YYYYMMDDHH.h5
amsua_n19_obs_YYYYMMDDHH.h5

> obs2ioda-v2.x -i input_dir -o output_dir gdas.airsev.tHHz.YYYYMMDD.bufr

Example output files:
amsua_aqua_obs_YYYYMMDDHH.h5
Example output files:
amsua_aqua_obs_YYYYMMDDHH.h5

> obs2ioda-v2.x -i input_dir -o output_dir gdas.1bmhs.tHHz.YYYYMMDD.bufr

Example output files:
mhs_metop-a_obs_YYYYMMDDHH.h5
mhs_metop-b_obs_YYYYMMDDHH.h5
mhs_n18_obs_YYYYMMDDHH.h5
mhs_n19_obs_YYYYMMDDHH.h5
Example output files:
mhs_metop-a_obs_YYYYMMDDHH.h5
mhs_metop-b_obs_YYYYMMDDHH.h5
mhs_n18_obs_YYYYMMDDHH.h5
mhs_n19_obs_YYYYMMDDHH.h5

> obs2ioda-v2.x -i input_dir -o output_dir gdas.mtiasi.tHHz.YYYYMMDD.bufr

**the following CRTM SpcCoeff files in little_endian must be present in the working directory for IASI radiance to brightness temperature conversion**
iasi_metop-a.SpcCoeff.bin -> iasi616_metop-a.SpcCoeff.bin
iasi_metop-b.SpcCoeff.bin -> iasi616_metop-b.SpcCoeff.bin
iasi_metop-c.SpcCoeff.bin -> iasi616_metop-c.SpcCoeff.bin
**the following CRTM SpcCoeff files in little_endian must be present in the working directory for IASI radiance to brightness temperature conversion**
iasi_metop-a.SpcCoeff.bin -> iasi616_metop-a.SpcCoeff.bin
iasi_metop-b.SpcCoeff.bin -> iasi616_metop-b.SpcCoeff.bin
iasi_metop-c.SpcCoeff.bin -> iasi616_metop-c.SpcCoeff.bin

Example output files:
iasi_metop-a_obs_YYYYMMDDHH.h5
iasi_metop-b_obs_YYYYMMDDHH.h5
iasi_metop-c_obs_YYYYMMDDHH.h5
Example output files:
iasi_metop-a_obs_YYYYMMDDHH.h5
iasi_metop-b_obs_YYYYMMDDHH.h5
iasi_metop-c_obs_YYYYMMDDHH.h5

> obs2ioda-v2.x -i input_dir -o output_dir gdas.crisf4.tHHz.YYYYMMDD.bufr

**the following CRTM SpcCoeff files in little_endian must be present in the working directory for CrIS radiance to brightness temperature conversion**
_for **cris** bufr file_
cris_npp.SpcCoeff.bin -> cris399_npp.SpcCoeff.bin
cris_n20.SpcCoeff.bin -> cris399_n20.SpcCoeff.bin
_for **crisf4** bufr file_
cris_npp.SpcCoeff.bin -> cris-fsr431_npp.SpcCoeff.bin
cris_n20.SpcCoeff.bin -> cris-fsr431_n20.SpcCoeff.bin
**the following CRTM SpcCoeff files in little_endian must be present in the working directory for CrIS radiance to brightness temperature conversion**
_for **cris** bufr file_
cris_npp.SpcCoeff.bin -> cris399_npp.SpcCoeff.bin
cris_n20.SpcCoeff.bin -> cris399_n20.SpcCoeff.bin
_for **crisf4** bufr file_
cris_npp.SpcCoeff.bin -> cris-fsr431_npp.SpcCoeff.bin
cris_n20.SpcCoeff.bin -> cris-fsr431_n20.SpcCoeff.bin

Example output files:
cris_npp_obs_YYYYMMDDHH.h5
cris_n20_obs_YYYYMMDDHH.h5
Example output files:
cris_npp_obs_YYYYMMDDHH.h5
cris_n20_obs_YYYYMMDDHH.h5

> obs2ioda-v2.x -i input_dir -o output_dir gdas.gpsro.tHHz.YYYYMMDD.bufr

Example output files:
gnssro_obs_YYYYMMDDHH.h5
Example output files:
gnssro_obs_YYYYMMDDHH.h5

## Converting Himawari Standard Data (HSD) FLDK files
```
Usage: obs2ioda-v2.x -i input_dir -ahi -t YYYYMMDDHHNN [-s 3]
Usage: obs2ioda-v2.x -i input_dir -ahi -t YYYYMMDDHHNN [-s 3] [-superob 3]
```

Input files are a list of Himawari Standard Data, e.g. HS_H08_20200815_0000_B14_FLDK_R20_S0210.DAT in the input_dir.
Minute must be specified in the time (-t) option.
Number of pixels to skip is optional.
Input files are a list of Himawari Standard Data, e.g. HS_H08_20200815_0000_B14_FLDK_R20_S0210.DAT in the input_dir.
Minute must be specified in the time (-t) option.
Number of pixels to skip is optional (-s). Default is 1, thus no thinning.
Do superobbing is optional, default is false. If -superob is included in the command line, then, no thinning is done and the parameter superob_halfwidth [-superob superob_halfwidth] should also be specified such that superob_width = 2*superob_halfwidth+1.

## Notes
* The output prefix (before _obs) is defined in define_mod.f90
Expand All @@ -110,9 +111,9 @@ through subroutines set_obtype_conv, set_name_satellite, set_name_sensor.
* The ob errors of AMSU-A/MHS radiances are coded in define_mod.f90. This should be changed in the future to read in from an external error table.
* The ob errors of satwnd-decoded AMVs are from an external error table (obs_errtable).
* Subroutine filter_obs_conv applies some additional QC for conventional observations as in GSI's read_prepbufr.f90 for the global model and can be de-activated through ``-noqc`` command-line option.
100 is added to the @PreQC value when the ob is flagged as not-use by filter_obs_conv.
100 is added to the @PreQC value when the ob is flagged as not-use by filter_obs_conv.
100 is chosen to make the original prepbufr quality marker easily readable.
* Subroutine filter_obs_satwnd applies QC for GOES-16/GOES-17 AMVs as in GSI's read_satwnd.f90.
* Subroutine filter_obs_satwnd applies QC for GOES-16/GOES-17 AMVs as in GSI's read_satwnd.f90.
@PreQC value is set to 15 for rejected obs.


3 changes: 3 additions & 0 deletions obs2ioda-v2/src/hsd.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -543,6 +543,7 @@ subroutine read_HSD(ccyymmddhhnn, inpdir, do_superob, superob_halfwidth)
call get_julian_time(iyear, imonth, iday, ihour, imin, isec, gstime, epochtime)

! do superobbing of ahi_himawari observations
! Superobbing code is based on WRFDA's AHI superobbing code from https://github.com/wrf-model/WRF/blob/develop/var/da/da_radiance/da_read_obs_netcdf4ahi_jaxa.inc
! npixel => x direction => nx
! nline => y direction => ny
if ( do_superob ) then
Expand Down Expand Up @@ -707,6 +708,7 @@ subroutine read_HSD(ccyymmddhhnn, inpdir, do_superob, superob_halfwidth)
end do y_loop

else

nlocs = 0
do jj = 1, nline , subsample
do ii = 1, npixel , subsample
Expand Down Expand Up @@ -747,6 +749,7 @@ subroutine read_HSD(ccyymmddhhnn, inpdir, do_superob, superob_halfwidth)
end do
end if

iloc = 0
! do thinning every subsample pixels
do jj = 1, nline, subsample
do ii = 1, npixel, subsample
Expand Down