!$Id$ ! ! ! ! ! module viirs_read_mod use hdf5 use readh5dataset implicit none private public :: get_viirs_data integer, parameter, public:: int1 = selected_int_kind(1) integer, parameter, public:: int2 = selected_int_kind(3) integer, parameter, public:: int4 = selected_int_kind(8) integer, parameter, public:: int8 = selected_int_kind(10) integer, parameter, public:: real4 = selected_real_kind(6,37) integer, parameter, public:: real8 = selected_real_kind(15,307) integer, parameter, public:: ipre = real4 ! - bowtie gaps values integer, parameter :: Ny_Pattern = 48 integer, parameter :: Nx_Pattern = 3200 integer(kind=int1), dimension(:,:), allocatable, public:: Gap_Pixel_Mask_pattern integer(kind=int4), dimension(:,:), allocatable, public:: Gap_Line_Idx_pattern type :: syms integer:: NO = 0 integer :: YES = 1 end type syms type(syms) ::sym logical :: gap_pattern_computed =.false. type , public :: viirs_data_config logical :: chan_on_rfl_mband ( 16 ) logical :: chan_on_rad_mband ( 16 ) logical :: chan_on_iband ( 5 ) logical :: chan_on_dnb logical :: viirs_cloud_mask_on logical :: viirs_cloud_type_on character ( len = 256 ) :: gmtco_file integer , dimension( 2 ) :: offset integer , dimension( 2 ) :: count character ( len =255 ) :: dir_1b character ( len =255 ) :: file_gmtco_base character (len = 200 ) :: Ancil_Data_Dir end type viirs_data_config type :: geo_str real , dimension (:,:) , allocatable :: solzen real , dimension (:,:) , allocatable :: satzen real , dimension (:,:) , allocatable :: solaz real , dimension (:,:) , allocatable :: sataz real , dimension (:,:) , allocatable :: relaz real , dimension (:,:) , allocatable :: lunaz real , dimension (:,:) , allocatable :: lunzen real , dimension (:,:) , allocatable :: lunrelaz real , dimension (:,:) , allocatable :: lat real , dimension (:,:) , allocatable :: lon real , dimension (:) , allocatable :: scan_time integer , dimension (:) , allocatable :: ascend real :: Moon_Illum_Frac real :: Moon_Phase_Angle end type geo_str type :: mband_str real,dimension (:,:) , allocatable :: ref real,dimension (:,:) , allocatable :: rad real,dimension (:,:) , allocatable :: bt end type mband_str type :: dnb_on_mband_str real,dimension (:,:) , allocatable :: ref real,dimension (:,:) , allocatable :: rad end type dnb_on_mband_str type :: iband_str real,dimension (:,:) , allocatable :: ref real,dimension (:,:) , allocatable :: bt end type iband_str type :: dnb_str real,dimension (:,:) , allocatable :: ref_lun real,dimension (:,:) , allocatable :: ref_sol real,dimension (:,:) , allocatable :: rad end type dnb_str type :: cloud_products_str integer , dimension(:,:) , allocatable :: cld_mask integer , dimension(:,:) , allocatable :: cld_phase integer , dimension(:,:) , allocatable :: cld_type end type cloud_products_str type :: gap_str integer(kind = int1) , dimension ( : , : ) , allocatable :: mask end type gap_str type, public :: viirs_data_out type (mband_str), dimension(16) :: mband type (iband_str), dimension(5) :: iband type ( geo_str) :: geo type (dnb_str) :: dnb type ( dnb_on_mband_str ) :: dnb_mgrid type ( cloud_products_str) :: prd type ( gap_str ) :: gap real , dimension (:,:) , allocatable :: ref_dnb integer , dimension (:,:) , allocatable :: cld_mask integer , dimension (:,:) , allocatable :: cld_phase integer , dimension (:,:) , allocatable :: cld_type end type viirs_data_out contains ! ! entree point routine ! subroutine get_viirs_data ( config , out ) type ( viirs_data_config ) , intent ( inout ) :: config type ( viirs_data_out ) , intent ( out ) :: out integer :: error_out call check_input (config ) ! -- call call read_viirs_level1b ( config, out , error_out) end subroutine get_viirs_data ! ! check input - should be extended ( testing if file exists etc... ! subroutine check_input ( config ) type ( viirs_data_config ) , intent ( inout ) :: config if ( config % offset (1) <= 0 ) config % offset (1) = 1 if ( config % offset (2) <= 0 ) config % offset (2) = 1 if ( config % count (1) <= 0 ) config % count (1) = 3200 if ( config % count (2) <= 0 ) config % count (2) = 200 end subroutine check_input ! ! ! subroutine read_viirs_level1b (config, out, error_out ) use hdf5 use readh5dataset use file_tools type ( viirs_data_config ) , intent ( in ) :: config type ( viirs_data_out ) , intent ( out ) :: out integer, intent(out) , optional :: error_out character(len=150) :: orbit_identifier integer(kind=int4) :: nx_start , nx_end , ny_start , ny_end integer(kind=int4) :: nx_start_iband , ny_start_iband real, dimension( : , : ) , pointer :: r2d_buffer real, dimension( : , : ) , allocatable :: r2d_data integer ( kind = int8) , dimension(:) , pointer :: i1d_buffer integer, dimension ( : , : ) , pointer :: i2d_buffer integer, dimension ( :, :), allocatable :: cld_type_idps character ( len = 256 ) , pointer , dimension ( :) :: file_arr_dummy character ( len =100 ) , dimension ( 7 ) :: setname_gmtco_list integer :: i_gmtco character ( len =150 ) :: file_gmtco integer :: i_mband character ( len =100) :: setname_mband , setname_mband_fac character ( len =150 ) :: file_mband logical, dimension ( 16) :: is_mband_on integer :: i_iband character ( len =100) :: setname_iband, setname_iband_fac character ( len =150 ) :: file_iband logical, dimension ( 5) :: is_iband_on character ( len =150 ) :: file_dnb_idx character ( len =150 ) :: file_gdnbo, file_svdnb integer :: io_err_stat integer :: i_map integer , dimension (3200) :: d2m_indx character ( len =100) :: setname_iicmo character ( len =150 ) :: file_iicmo integer :: lun integer :: n_files integer :: k ! - IR channel wavenumbers real , parameter :: nu_12 = 2708.3865 real , parameter :: nu_13 = 2460.8193 real , parameter :: nu_14 = 1166.1845 real , parameter :: nu_15 = 935.10476 real , parameter :: nu_16 = 845.79752 real, dimension(16) :: nu_list character ( len = 2 ) :: band_nr_file character ( len = 2 ) :: band_nr_var real , dimension(:) , pointer :: factors logical , dimension(16) :: data_scaled_mband integer, dimension(2) :: dim_buf logical, dimension(:,:), allocatable :: invalid_pixel integer(kind=int4) , dimension(2) :: dim_seg integer(kind=int4) , dimension(2) :: dim_seg_iband ! - this is to delete integer :: i integer , dimension (: ) , allocatable:: time_msec_day integer (kind = int8) , parameter :: microsec_per_day = 86400000000 integer, dimension(2)::shape_buffer real :: missing_value_real4 = -999. integer :: missing_value_int = -999 integer:: scaled_Missing integer :: unscaled_Missing integer , dimension(2) :: offset_mband , offset_iband Scaled_Missing = 65528 error_out = 0 Unscaled_Missing = -999.3 ! - gmtco geo products setname_gmtco_list = [ & 'All_Data/VIIRS-MOD-GEO-TC_All/Latitude' & ! 1 , 'All_Data/VIIRS-MOD-GEO-TC_All/Longitude' & ! 2 , 'All_Data/VIIRS-MOD-GEO-TC_All/MidTime' & ! 3 , 'All_Data/VIIRS-MOD-GEO-TC_All/SatelliteAzimuthAngle' & ! 4 , 'All_Data/VIIRS-MOD-GEO-TC_All/SatelliteZenithAngle' & ! 5 , 'All_Data/VIIRS-MOD-GEO-TC_All/SolarAzimuthAngle' & ! 6 , 'All_Data/VIIRS-MOD-GEO-TC_All/SolarZenithAngle' & ! 7 ] nx_start = config % offset(1) ny_start = config % offset(2) nx_end = nx_start + config % count(1) - 1 ny_end = ny_start + config % count(2) - 1 dim_seg = [ nx_end - nx_start +1, ny_end - ny_start+1 ] file_gmtco = trim(config % dir_1b) // trim(config % file_gmtco_base) offset_mband = [ nx_start -1 , ny_start - 1] do i_gmtco = 1 , 7 if ( i_gmtco == 3 ) then call h5readdataset ( file_gmtco , setname_gmtco_list(i_gmtco) , i1d_buffer ) else call h5readdataset ( file_gmtco , setname_gmtco_list(i_gmtco) ,offset_mband , dim_seg, r2d_buffer ) where ( r2d_buffer < unscaled_missing ) r2d_buffer = missing_value_real4 end if select case (i_gmtco) case(1) if (.not. allocated ( out % geo % lat) ) allocate ( out % geo % lat (dim_seg(1), dim_seg(2)) ) out % geo % lat = r2d_buffer case(2) if (.not. allocated ( out % geo % lon) ) allocate ( out % geo % lon (dim_seg(1), dim_seg(2)) ) out % geo %lon = r2d_buffer case(3) allocate ( time_msec_day ( size ( i1d_buffer))) time_msec_day = ( mod ( i1d_buffer , microsec_per_day ) ) / 1000 if (.not. allocated ( out % geo % scan_time ) ) allocate ( out % geo % scan_time (dim_seg(2)) ) out % geo % scan_time = (/( time_msec_day( ( k -1 ) /16 + 1) , k =ny_start , ny_end )/) deallocate ( time_msec_day ) case(4) if (.not. allocated ( out % geo % sataz) ) allocate ( out % geo % sataz (dim_seg(1), dim_seg(2)) ) out % geo % sataz = r2d_buffer case(5) if (.not. allocated ( out % geo % satzen) ) allocate ( out % geo % satzen (dim_seg(1), dim_seg(2)) ) out % geo % satzen = r2d_buffer case(6) if (.not. allocated ( out % geo % solaz) ) allocate ( out % geo % solaz (dim_seg(1), dim_seg(2)) ) out % geo % solaz = r2d_buffer case(7) if (.not. allocated ( out % geo % solzen) ) allocate ( out % geo % solzen (dim_seg(1), dim_seg(2)) ) out % geo % solzen = r2d_buffer end select if (i_gmtco /= 3) deallocate ( r2d_buffer) if (i_gmtco == 3) deallocate ( i1d_buffer) end do ! - mbands orbit_identifier = trim ( config %file_gmtco_base(11:37))//"*" is_mband_on = config % chan_on_rfl_mband data_scaled_mband = .true. ! - channel 13 is the only without Factors in VIRRS file data_scaled_mband(13) = .false. ! IR channels wavenumbers ! updated on 2013/03/20 AW nu_list = 0 nu_list(12:16) = [nu_12 , nu_13 , nu_14 , nu_15 , nu_16 ] do i_mband = 1 , 16 if ( .not. is_mband_on(i_mband) ) cycle write ( band_nr_file , '(i2.2)' ) i_mband write ( band_nr_var , '(i2)' ) i_mband if ( i_mband < 10) write ( band_nr_var , '(i1)' ) i_mband file_arr_dummy => file_search (trim(config %dir_1b), 'SVM'//trim(band_nr_file)//'_*'//trim(orbit_identifier) , n_files ) if ( n_files == 0 ) cycle file_mband = file_arr_dummy(1) setname_mband = trim('All_Data/VIIRS-M'//trim(band_nr_var)//'-SDR_All/Reflectance' ) if ( i_mband > 11) setname_mband = trim('All_Data/VIIRS-M'//trim(band_nr_var)//'-SDR_All/Radiance' ) setname_mband_fac = trim(setname_mband)//'Factors' if ( data_scaled_mband ( i_mband ) ) then call h5readdataset ( trim(config %dir_1b)//trim(file_mband), setname_mband , offset_mband , dim_seg , i2d_buffer ) dim_buf = shape ( i2d_buffer ) allocate ( invalid_pixel ( dim_buf(1) , dim_buf(2) ) ) invalid_pixel = i2d_buffer >= scaled_missing call h5readdataset ( trim(config %dir_1b)//trim(file_mband), setname_mband_fac , factors ) allocate ( r2d_data ( dim_buf(1) , dim_buf(2) )) r2d_data = i2d_buffer * factors(1) + factors(2) deallocate ( i2d_buffer ) else call h5readdataset ( trim(config %dir_1b)//trim(file_mband), setname_mband , offset_mband , dim_seg , r2d_buffer ) dim_buf = shape ( r2d_buffer ) allocate ( invalid_pixel ( dim_buf(1) , dim_buf(2) ) ) invalid_pixel = r2d_buffer <= unscaled_missing allocate ( r2d_data ( dim_buf(1) , dim_buf(2) )) r2d_data = r2d_buffer deallocate ( r2d_buffer) end if if ( i_mband <= 11) then if (.not. allocated ( out % mband (i_mband) % ref) ) allocate ( out % mband (i_mband) % ref(dim_seg(1), dim_seg(2)) ) out % mband ( i_mband) % ref = 100. * r2d_data where ( invalid_pixel ) out % mband ( i_mband) % ref = missing_value_real4 else call convert_viirs_radiance ( r2d_data , nu_list(i_mband) , missing_value_real4 ) if (.not. allocated ( out % mband (i_mband) % rad) ) allocate ( out % mband (i_mband) % rad(dim_seg(1), dim_seg(2)) ) out % mband ( i_mband) % rad = r2d_data where ( invalid_pixel ) out % mband ( i_mband) % rad = missing_value_real4 end if if ( allocated ( invalid_pixel) ) deallocate ( invalid_pixel ) if ( allocated (r2d_data) ) deallocate ( r2d_data ) end do !- ibands is_iband_on =config % chan_on_iband dim_seg_iband = dim_seg * 2 nx_start_iband = 1 ny_start_iband = ny_start * 2 - 1 offset_iband = [ nx_start_iband -1 , ny_start_iband - 1 ] do i_iband = 1 , 5 if ( .not. is_iband_on(i_iband) ) cycle write ( band_nr_file , '(i2.2)' ) i_iband write ( band_nr_var , '(i1)' ) i_iband file_arr_dummy => file_search (trim(config %dir_1b), 'SVI'//trim(band_nr_file)//'_*'//trim(orbit_identifier) , n_files ) if ( n_files == 0 ) cycle file_iband = file_arr_dummy(1) setname_iband = trim('All_Data/VIIRS-I'//trim(band_nr_var)//'-SDR_All/Reflectance' ) if ( i_iband > 3) setname_iband = trim('All_Data/VIIRS-I'//trim(band_nr_var)//'-SDR_All/BrightnessTemperature' ) setname_iband_fac = trim(setname_iband)//'Factors' call h5readdataset ( trim(config %dir_1b)//file_iband, setname_iband , offset_iband , dim_seg_iband , r2d_buffer ) dim_buf = shape ( r2d_buffer ) allocate ( invalid_pixel ( dim_buf(1) , dim_buf(2) ) ) invalid_pixel = r2d_buffer >= scaled_missing call h5readdataset ( trim(config %dir_1b)//trim(file_iband), setname_iband_fac , factors ) r2d_buffer = r2d_buffer * factors(1) + factors(2) if ( i_iband <= 3) then if (.not. allocated ( out % iband (i_iband) % ref) ) allocate ( out % iband (i_iband) % ref(dim_seg_iband(1), dim_seg_iband(2)) ) out % iband ( i_iband) % ref = 100. * r2d_buffer where ( invalid_pixel ) out % iband ( i_iband) % ref = missing_value_real4 else if (.not. allocated ( out % iband (i_iband) % bt) ) allocate ( out % iband (i_iband) % bt (dim_seg_iband(1), dim_seg_iband(2)) ) out % iband ( i_iband) % bt = r2d_buffer where ( invalid_pixel ) out % iband ( i_iband) % bt = missing_value_real4 end if deallocate ( invalid_pixel ) deallocate ( r2d_buffer ) end do ! - dnb ! - if (config % chan_on_dnb ) then ! - mapping file ( maps from dnb to M-bands resolution) file_dnb_idx = trim(config % Ancil_Data_Dir)//'viirs/dnb2m_indx.txt' lun = getlun() open (unit = lun , file=trim ( file_dnb_idx) , status="old",action="read") read (unit = lun , fmt=* , iostat = io_err_stat) d2m_indx close (unit = lun) !- gdnbo products file_arr_dummy => file_search (trim(config %dir_1b), 'GDNBO*'//trim(orbit_identifier) , n_files ) if ( n_files == 0 ) stop 98 file_gdnbo = file_arr_dummy(1) call h5readdataset ( trim(config %dir_1b)//file_gdnbo,'All_Data/VIIRS-DNB-GEO_All/LunarAzimuthAngle' , r2d_buffer ) where ( r2d_buffer < unscaled_missing ) r2d_buffer = missing_value_real4 do i_map = 1, dim_seg(1) if (.not. allocated ( out % geo % lunaz) ) allocate ( out % geo % lunaz (dim_seg(1), dim_seg(2)) ) out % geo % lunaz ( i_map , 1 : ( ny_end - ny_start + 1 ) ) = r2d_buffer ( d2m_indx (i_map) , ny_start : ny_end ) end do call h5readdataset ( trim(config %dir_1b)//file_gdnbo,'All_Data/VIIRS-DNB-GEO_All/LunarZenithAngle' , r2d_buffer ) where ( r2d_buffer < unscaled_missing ) r2d_buffer = missing_value_real4 do i_map = 1, dim_seg(1) if (.not. allocated ( out % geo % lunzen ) ) allocate ( out % geo % lunzen (dim_seg(1), dim_seg(2)) ) out % geo % lunzen ( i_map , 1 : ( ny_end - ny_start + 1 ) ) = r2d_buffer (d2m_indx ( i_map ) ,ny_start : ny_end ) end do allocate ( out % geo % lunrelaz ( dim_seg(1) , dim_seg(2) ) ) call compute_relative_azimuth ( out % geo % lunaz , out % geo % sataz, out % geo % lunrelaz ) call H5ReadDataset( trim(config %dir_1b)//file_gdnbo , 'All_Data/VIIRS-DNB-GEO_All/MoonIllumFraction', out % geo % Moon_Illum_Frac) call H5ReadDataset( trim(config %dir_1b)//file_gdnbo , 'All_Data/VIIRS-DNB-GEO_All/MoonPhaseAngle', out % geo % Moon_Phase_Angle) file_arr_dummy => file_search (trim(config %dir_1b), 'SVDNB*'//trim(orbit_identifier) , n_files ) if ( n_files == 0 ) stop 97 file_svdnb = file_arr_dummy(1) call h5readdataset ( trim(config %dir_1b)//file_svdnb, 'All_Data/VIIRS-DNB-SDR_All/Radiance' , r2d_buffer ) where ( r2d_buffer < unscaled_missing ) r2d_buffer = missing_value_real4 ! --- remap dnb using indexes if (.not. allocated ( out % dnb_mgrid % rad ) ) allocate ( out % dnb_mgrid % rad (dim_seg(1), dim_seg(2)) ) do i_map = 1, dim_seg(1) out % dnb_mgrid % rad(i_map , 1 : ( ny_end - ny_start + 1 ) ) = r2d_buffer(d2m_indx( i_map ) , ny_start : ny_end ) end do if (.not. allocated ( out % dnb_mgrid % ref ) ) allocate ( out % dnb_mgrid % ref (dim_seg(1), dim_seg(2)) ) call convert_rad_2_sol_ref_dnb ( out % dnb_mgrid % rad & , out % geo % solzen & , missing_value_real4 & , out % dnb_mgrid % ref ) end if ! - cloud products - should be also on a different location! ! - cloud_mask_1b is bad defined, which values can it have sym%YES and NO or 0,1,2? ! - goal is to populate cld_mask_aux , cld_phase_aux and cld_type_aux if ( config % viirs_cloud_mask_on ) then file_arr_dummy => file_search (trim(config %dir_1b), 'IICMO*'//trim(orbit_identifier) , n_files ) if ( n_files == 0 ) stop 97 file_iicmo = file_arr_dummy(1) setname_iicmo='All_Data/VIIRS-CM-IP_All/QF1_VIIRSCMIP' call h5readdataset ( trim(config %dir_1b)//file_iicmo, setname_iicmo , offset_mband , dim_seg , i2d_buffer ) allocate ( out % prd % cld_mask (dim_seg(1), dim_seg(2)) ) out % prd % cld_mask = 0 where ( btest ( i2d_buffer , 2 ) ) out % prd % cld_mask = out % prd % cld_mask + 1 end where where ( btest ( i2d_buffer , 3 ) ) out % prd % cld_mask = out % prd % cld_mask + 2 end where deallocate ( i2d_buffer) ! ! setname_iicmo='All_Data/VIIRS-CM-IP_All/QF6_VIIRSCMIP' call h5readdataset ( trim(config %dir_1b)//file_iicmo, setname_iicmo , offset_mband , dim_seg , i2d_buffer ) shape_buffer = shape (i2d_buffer) allocate ( cld_type_idps ( shape_buffer(1), shape_buffer(2))) cld_type_idps = log2int ( btest( i2d_buffer , 0 ) ) & + 2.* log2int(btest( i2d_buffer , 1 ) ) & + 4. * log2int(btest( i2d_buffer , 2 ) ) deallocate ( i2d_buffer) allocate ( out % prd % cld_phase ( dim_seg(1) , dim_seg(2) ) ) out % prd % cld_phase = 5 where ( cld_type_idps == 1 .or. cld_type_idps == 2 ) out % prd % cld_phase = 0 end where where ( cld_type_idps == 3 ) out % prd % cld_phase = 1 end where where ( cld_type_idps == 4 ) out % prd % cld_phase = 2 end where where ( cld_type_idps == 5 .or. cld_type_idps == 6 .or. cld_type_idps == 7 ) out % prd % cld_phase = 4 end where ! - type allocate ( out % prd % cld_type ( dim_seg(1), dim_seg(2) )) out % prd % cld_type = 10 where ( cld_type_idps == 1 ) out % prd % cld_type = 0 end where where ( cld_type_idps == 2 ) out % prd % cld_type = 1 end where where ( cld_type_idps == 3 ) out % prd % cld_type = 3 end where where ( cld_type_idps == 4 ) out % prd % cld_type = 4 end where where ( cld_type_idps == 5 ) out % prd % cld_type = 6 end where where ( cld_type_idps == 6 ) out % prd % cld_type = 7 end where where ( cld_type_idps == 7 ) out % prd % cld_type = 8 end where deallocate ( cld_type_idps ) where ( out % prd % cld_type == 10 ) out % prd % cld_mask = missing_value_int end where ! ! ! end if ! -horrible global parameters this has to be gone soon! !Num_Scans_Read = Ny_End - Ny_Start + 1 ! Fill in Scan_Number !do i = 1, num_scans_per_segment ! scan_number(i) = ny_start + i - 1 !end do if ( .not. gap_pattern_computed ) then call compute_viirs_bowtie_gap_pattern() end if call fill_viirs_bowtie_gaps ( ny_start , dim_seg(2) , out ) ! xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx ! everything between these xxx have to go to other locations !- ascending (global varaibel ) allocate ( out % geo % ascend ( dim_seg (2) ) ) out % geo % ascend = 0 do i = 1 , ny_end -ny_start - 1 if ( out%geo%lat (dim_seg(1) / 2 , i + 1) <= out%geo%lat(dim_seg(1) / 2 , i ) ) out % geo % ascend ( i ) = 1 end do allocate ( out % geo % relaz ( dim_seg(1) , dim_seg(2) ) ) ! rel azimuths - these are all global variables call compute_relative_azimuth ( out%geo%solaz , out%geo%sataz , out%geo%relaz) ! compute glintzen - these are all global variables !call COMPUTE_GLINTZEN_SCATANGLE(out%geo%Solzen,out%geo%Satzen,out%geo%Relaz, & ! Glintzen,Scatangle) end subroutine read_viirs_level1b !------------------------------------------------------------------------------------- ! subroutine to compute the bowtie gap pattern. ! This pattern repeats every 48 scans ! all VIIRS files should be integer multiples of this pattern ! ! Gap_Pixel_Mask_Pattern = a binary mask that identifies these bowtie gaps ! Gap_Line_Idx = line index for each pixel in pattern including gap pixels ! ! A description of the pattern ! !-------- line type 3 !---- line type 4 ! 12 lines without gaps !---- line type 1 !-------- line type 2 !-------- line type 3 !---- line type 4 ! 12 lines without gaps !---- line type 1 !-------- line type 2 !-------- line type 3 !---- line type 4 ! 12 lines without gaps !---- line type 1 !-------- line type 2 ! ! line types 1 and 4 have 1280 missing pixels ! line types 2 and 3 have 2016 missing pixels ! !------------------------------------------------------------------------------------- subroutine COMPUTE_VIIRS_BOWTIE_GAP_PATTERN() integer (kind=int4), dimension(Ny_Pattern):: Line_Type integer (kind=int4), parameter:: Ngap_1 = 640 !1280 integer (kind=int4), parameter:: Ngap_2 = 1008 !2016 integer (kind=int4), parameter:: Ngap_3 = 1008 !2016 integer (kind=int4), parameter:: Ngap_4 = 640 !1280 integer (kind=int4):: iline integer (kind=int4):: i1 integer (kind=int4):: i2 !--- define the line patterns as described above Line_Type = (/3,4,0,0,0,0,0,0,0,0,0,0,0,0,1,2,3,4,0,0,0,0,0,0,0,0,0,0,0,0,1,2,3,4,0,0,0,0,0,0,0,0,0,0,0,0,1,2/) if (.not. allocated(Gap_Line_Idx_Pattern)) Allocate(Gap_Line_Idx_Pattern(Nx_Pattern,Ny_Pattern)) if (.not. allocated(Gap_Pixel_Mask_Pattern)) Allocate(Gap_Pixel_Mask_Pattern(Nx_Pattern,Ny_Pattern)) do iline = 1 , Ny_Pattern Gap_Line_Idx_Pattern( : , iline ) = -999 Gap_Pixel_Mask_Pattern( : , iline ) = sym%NO if (line_Type(iline) == 1) then i1 = 1 i2 = Ngap_1 Gap_Line_Idx_Pattern( i1 : i2 , iline) = iline - 1 Gap_Pixel_Mask_Pattern( i1 : i2 , iline) = sym%YES i1 = Nx_Pattern - Ngap_1 + 1 i2 = Nx_Pattern Gap_Line_Idx_Pattern( i1 : i2 , iline ) = iline - 1 Gap_Pixel_Mask_Pattern( i1 : i2 , iline ) = sym%YES end if if (Line_Type(iline) == 2) then i1 = 1 i2 = Ngap_1 Gap_Line_Idx_Pattern( i1 : i2 , iline ) = iline - 2 Gap_Pixel_Mask_Pattern( i1 : i2 , iline ) = sym%YES i1 = Ngap_1 + 1 i2 = Ngap_2 Gap_Line_Idx_Pattern( i1 : i2 , iline ) = iline - 1 Gap_Pixel_Mask_Pattern( i1 : i2 , iline ) = sym%YES i1 = Nx_Pattern - Ngap_1 + 1 i2 = Nx_Pattern Gap_Line_Idx_Pattern( i1 : i2 , iline ) = iline - 2 Gap_Pixel_Mask_Pattern( i1 : i2 , iline ) = sym%YES i1 = Nx_Pattern - Ngap_2 + 1 i2 = i1 + (Ngap_2 - Ngap_1) Gap_Line_Idx_Pattern( i1 : i2 , iline ) = iline - 1 Gap_Pixel_Mask_Pattern( i1 : i2 , iline ) = sym%YES end if if (Line_Type(iline) == 3) then i1 = 1 i2 = Ngap_4 Gap_Line_Idx_Pattern(i1:i2,iline) = iline + 2 Gap_Pixel_Mask_Pattern(i1:i2,iline) = sym%YES i1 = Ngap_4 + 1 i2 = Ngap_3 Gap_Line_Idx_Pattern(i1:i2,iline) = iline + 1 Gap_Pixel_Mask_Pattern(i1:i2,iline) = sym%YES i1 = Nx_Pattern - Ngap_4 + 1 i2 = Nx_Pattern Gap_Line_Idx_Pattern(i1:i2,iline) = iline + 2 Gap_Pixel_Mask_Pattern(i1:i2,iline) = sym%YES i1 = Nx_Pattern - Ngap_3 + 1 i2 = i1 + (Ngap_3 - Ngap_4) Gap_Line_Idx_Pattern(i1:i2,iline) = iline + 1 Gap_Pixel_Mask_Pattern(i1:i2,iline) = sym%YES end if if (Line_Type(iline) == 4) then i1 = 1 i2 = Ngap_4 Gap_Line_Idx_Pattern(i1:i2,iline) = iline + 1 Gap_Pixel_Mask_Pattern(i1:i2,iline) = sym%YES i1 = Nx_Pattern - Ngap_4 + 1 i2 = Nx_Pattern Gap_Line_Idx_Pattern( i1 : i2 , iline ) = iline + 1 Gap_Pixel_Mask_Pattern(i1:i2,iline) = sym%YES end if end do end subroutine compute_viirs_bowtie_gap_pattern !------------------------------------------------------------------------------ ! this routine uses the bowtie gap pattern, applies it to an arbitrary ! segment of data and fills in the observations with the closest valid data !------------------------------------------------------------------------------ subroutine fill_viirs_bowtie_gaps ( Line_Start, Number_of_Lines , out ) integer(kind=int4), intent(in):: Line_Start integer(kind=int4), intent(in):: Number_of_Lines type ( viirs_data_out ) , intent (inout) :: out integer(kind=int4):: Line_Offset integer(kind=int4):: Line_in_Pattern integer(kind=int4):: Line_in_Segment integer(kind=int4):: Line_Idx integer(kind=int4):: Elem_Idx integer :: Num_Pix = 3200 integer , dimension(:,: ) , allocatable :: gap_line_idx integer( kind = int1) , dimension(:,: ) , allocatable :: gap_pixel_mask integer :: i_mband integer :: missing_value_int1 = -999 allocate ( gap_pixel_mask (3200, number_of_lines) , gap_line_idx(3200, number_of_lines)) do line_in_segment = 1 , number_of_lines line_in_pattern = mod (line_start -1 + line_in_segment , ny_pattern) if (Line_In_Pattern == 0) Line_In_Pattern = Ny_Pattern line_offset = line_in_segment - line_in_pattern gap_line_idx ( : , line_in_segment ) = Gap_Line_Idx_Pattern(:,Line_in_Pattern) + line_offset where ( gap_line_idx(:,Line_in_Segment) <= 0 ) Gap_Line_Idx(:,Line_in_Segment) = 1 end where where ( gap_line_idx(:,Line_in_Segment) > Number_of_Lines ) Gap_Line_Idx(:,Line_in_Segment) = Number_of_Lines end where ! - write gap mask to output Gap_Pixel_Mask(:,Line_in_Segment) = Gap_Pixel_Mask_Pattern(:,Line_in_Pattern) do Elem_Idx = 1, Num_Pix if (Gap_Pixel_Mask(Elem_Idx,Line_in_Segment)== sym%YES) then Line_Idx = Gap_Line_Idx(Elem_Idx,Line_in_Segment) do i_mband = 1 , 11 if ( .not. allocated( out % mband(i_mband) % ref )) cycle out % mband (i_mband) % ref (elem_idx , line_in_segment) = out%mband(i_mband)%ref(elem_idx,line_idx) end do do i_mband = 12,16 if (.not. allocated (out % mband(i_mband) % rad )) cycle out % mband(i_mband) % rad (elem_idx,line_in_segment) = out % mband(i_mband)%rad(elem_idx,line_idx) end do if ( allocated (out % prd % cld_mask)) then out % prd % cld_type(elem_idx,line_in_segment) = out % prd % cld_type(elem_idx,line_idx) out % prd % cld_phase(elem_idx,line_in_segment) = out % prd % cld_phase(elem_idx,line_idx) out % prd % cld_mask(Elem_Idx,Line_in_Segment) = out % prd % cld_mask(Elem_Idx,Line_Idx) end if end if end do end do if ( allocated ( out % prd % cld_mask ) ) then where( out % prd % cld_mask == Missing_Value_Int1) Gap_Pixel_Mask = sym%YES end where end if deallocate ( gap_line_idx) allocate ( out % gap % mask (3200, number_of_lines)) out % gap % mask = gap_pixel_mask deallocate ( gap_pixel_mask) end subroutine fill_viirs_bowtie_gaps ! this routine should be at a different place subroutine compute_relative_azimuth ( ang1 , ang2, rel_az_out ) real , dimension(:,:), intent(in) :: ang1 real , dimension(:,:), intent(in) :: ang2 real , dimension(:,:) :: rel_az_out rel_az_out = abs (ang1 - ang2 ) where ( rel_az_out > 180 ) rel_az_out = 180.0 - rel_az_out end subroutine compute_relative_azimuth !----------------------------------------------------------------------------------------- ! Extract time information from VIIRS filename - should explore use of header for this ! assumingly called from outside ... !----------------------------------------------------------------------------------------- SUBROUTINE READ_VIIRS_DATE_TIME(Infile,year,doy,start_time,end_time,orbit,orbit_identifier) ! Get the date & time from the file's name CHARACTER(Len=*), INTENT(IN) :: Infile INTEGER, INTENT(OUT) :: year INTEGER, INTENT(OUT) :: doy !day of year INTEGER, INTENT(OUT) :: start_time !millisec INTEGER, INTENT(OUT) :: end_time !millisec INTEGER, INTENT(OUT) :: orbit INTEGER :: month INTEGER :: day INTEGER :: start_hour INTEGER :: start_minute INTEGER :: start_sec INTEGER :: end_hour INTEGER :: end_minute INTEGER :: end_sec character(38), intent(out) :: orbit_identifier ! 1 2 3 4 5 6 7 8 !12345678901234567890123456789012345678901234567890123456789012345678901234567890 !GMODO_npp_d20100906_t2110510_e2112156_b00012_c20110707160532497848_noaa_ops.h5 ! --- Read data from the file name read(Infile(12:15), fmt="(I4)") year read(Infile(16:17), fmt="(I2)") month read(Infile(18:19), fmt="(I2)") day read(Infile(22:23), fmt="(I2)") start_hour read(Infile(24:25), fmt="(I2)") start_minute read(Infile(26:27), fmt="(I2)") start_sec read(Infile(31:32), fmt="(I2)") end_hour read(Infile(33:34), fmt="(I2)") end_minute read(Infile(35:36), fmt="(I2)") end_sec read(Infile(40:44), fmt="(I5)") orbit !---- store orbit number orbit_identifier = infile(7:44) !--- compute day of year call JULIAN(day,month,year,doy) ! --- Calculate start and end time start_time = ((start_hour * 60 + start_minute) * 60 + start_sec) * 1000 end_time = ((end_hour * 60 + end_minute) * 60 + end_sec) * 1000 END SUBROUTINE READ_VIIRS_DATE_TIME ! this is called from outside, why is this needed? !--------------------------------------------------------------------------------- SUBROUTINE GET_NUMBER_OF_SCANS_FROM_VIIRS(Infile,Number_Of_Viirs_Lines,Error_Out) CHARACTER(Len=*), INTENT(IN) :: Infile INTEGER(kind=int4), INTENT(OUT) :: Error_Out INTEGER(KIND=INT4), INTENT(OUT):: Number_of_Viirs_Lines CHARACTER(Len=100) :: Setname integer ,dimension(:), pointer ::test error_out = 0 Setname = 'All_Data/VIIRS-MOD-GEO-TC_All/NumberOfScans' call H5ReadDataset( infile ,setname, test) Number_of_Viirs_Lines = test(1) END SUBROUTINE GET_NUMBER_OF_SCANS_FROM_VIIRS ! ------------------------------------------------------------------ ! - missing comment !!! why is this needed? ! --------------------------------------------------------------------------------------- subroutine CONVERT_VIIRS_RADIANCE(radiance,nu,missing_value) real (kind=real4), dimension(:,:), intent(inout):: radiance real (kind=real4), intent(in):: nu real (kind=real4), intent(in):: missing_value where(radiance /= missing_value) radiance = radiance * (((10000.0 / nu )**2) / 10.0) end where return end subroutine CONVERT_VIIRS_RADIANCE !--------------------------------------------------------------------------------------- subroutine convert_rad_2_sol_ref_dnb ( radiance , solzen , missing_value , reflectance ) real (kind=real4), dimension(:,:), intent(in) :: radiance real (kind=real4), dimension(:,:), intent(in) :: solzen real (kind=real4), intent(in) :: missing_value real (kind=real4), dimension(:,:), intent(out) :: reflectance real (kind=real4), parameter :: fo = 0.044217282 ! dnb solar energy in w/cm2 ( Source? ) real , parameter :: PI = 3.14159265359 real , parameter :: DTOR = PI / 180. real :: sun_earth_distance real :: start_day = 100 sun_earth_distance = 1.0 - 0.016729 * cos ( 0.9856 * ( start_day - 4.0 ) * DTOR ) reflectance = missing_value where(radiance /= missing_value) reflectance = 100*(pi*radiance*sun_earth_distance**2) / (cosd(solzen)*fo) end where end subroutine convert_rad_2_sol_ref_dnb !--------------------------------------------------------------------------------------- ! transforms LOGICAL data type to an integer ! function log2int (b) result (r) logical , dimension(:,:) :: b integer , dimension(:,:),allocatable :: r integer :: dims(2) dims = shape ( b) allocate ( r (dims(1), dims(2) ) ) r = 0 where ( b ) r = 1 end where end function log2int !--------------------------------------------------------------------------------------- !------------------------------------------------- ! subroutine JULIAN(iday,imonth,iyear,jday) ! compute julian day ! input: ! iday - integer day ! imonth - integer month ! iyear - integer year (2 or four digits) ! ! output : jday - julian day !-------------------------------------------------- subroutine JULIAN(iday,imonth,iyear,jday) !-- Computes julian day (1-365/366) integer, intent(in):: iday,imonth,iyear integer, intent(out):: jday integer:: j integer, dimension(12):: jmonth jmonth = reshape ((/31,28,31,30,31,30,31,31,30,31,30,31/),(/12/)) jday = iday if (modulo(iyear,4) == 0) then jmonth(2)=29 endif do j = 1,imonth-1 jday = jday + jmonth(j) end do end subroutine JULIAN end module viirs_read_mod