From 7848630eab70afb86566c369e271d96a9b7edb07 Mon Sep 17 00:00:00 2001 From: ibanos90 Date: Fri, 28 Oct 2022 12:06:05 -0600 Subject: [PATCH] Remove do_thinning option. If superob is false, skip the number of subsample pixels. If super_ob is true, then no thinning is done. The readme files for abi namelist and ahi converter usage are updated to include more accurate information for users. --- goes_abi/README.md | 6 +- goes_abi/src/goes_abi_converter.f90 | 132 ++++++++++++++-------------- obs2ioda-v2/README.md | 113 ++++++++++++------------ obs2ioda-v2/src/hsd.f90 | 3 + 4 files changed, 129 insertions(+), 125 deletions(-) diff --git a/goes_abi/README.md b/goes_abi/README.md index 44386d3..76fd6ee 100644 --- a/goes_abi/README.md +++ b/goes_abi/README.md @@ -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: @@ -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: diff --git a/goes_abi/src/goes_abi_converter.f90 b/goes_abi/src/goes_abi_converter.f90 index 5ed8ef5..910cfb8 100644 --- a/goes_abi/src/goes_abi_converter.f90 +++ b/goes_abi/src/goes_abi_converter.f90 @@ -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 @@ -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 @@ -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 @@ -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) diff --git a/obs2ioda-v2/README.md b/obs2ioda-v2/README.md index 638b46c..4e3f168 100644 --- a/obs2ioda-v2/README.md +++ b/obs2ioda-v2/README.md @@ -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 @@ -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. diff --git a/obs2ioda-v2/src/hsd.f90 b/obs2ioda-v2/src/hsd.f90 index ce5f145..4f89689 100644 --- a/obs2ioda-v2/src/hsd.f90 +++ b/obs2ioda-v2/src/hsd.f90 @@ -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 @@ -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 @@ -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