! Name: $Id$ ! ! Purpose: ! Subroutines maia_read_Viirs, maia_Viirs_field_init ! Current Code Owner: NWP SAF ! ! COPYRIGHT ! This software was developed within the context of the EUMETSAT Satellite ! Application Facility on Numerical Weather Prediction (NWP SAF), under the ! Cooperation Agreement dated 29 June 2011, between EUMETSAT and the ! Met Office, UK, by one or more partners within the NWP SAF. The partners ! in the NWP SAF are the Met Office, ECMWF, KNMI and MeteoFrance. ! ! Copyright 2014, EUMETSAT, All Rights Reserved. ! ! History: ! Version Date Comment ! ! 1.0 01/05/2015 Lydie Lavanant ! 1.1 10/01/2017 Pascale Roquet add NWPSAF header and clean code ! 1.2 04/06/2017 Pascale Roquet test if "I" files are present ! 1.3 11/01/2018 Pascale Roquet update for JPSS ! 1.4 25/03/2024 Use attribute Platform_Short_Name instead of Mission_Name Sonia Pere update for CSPP4 subroutine maia_read_Viirs ( idbg, field_M, field_I) ! owner : MF/DP/CMS/R&D ! Authors : lydie lavanant ! date : 10/10/2011 ! Read in HDF5 format the contents of : ! VIIRS_M SDR and Geolocation ! VIIRS_I SDR and Geolocation ! fill the viirs_field structures ! fait quelques suppositions sur la composition des fichiers ! (imageur canaux 4 et 5 dans meme fichier) ! (à verifier et corriger lorsque on recevra des données reelles) #ifdef HAS_LIBHDF5 use hdf5 USE H5LT USE aapp_hdf_mod #endif ! 1.1 01/12/2014 Lydie Lavanant AAPP 7.9 : maia4 version 4.2 USE viirs_sdr_def USE viirs_sdr_mem USE viirs_sdr_io USE viirs_sdr_mod USE MOD_MAIA_CONST USE mod_maia_types #include "aapp_throw.h" IMPLICIT NONE type( debug ), intent(in) :: idbg ! Type (field), intent(inout) :: field_M Type (field), intent(inout) :: field_I #ifdef HAS_LIBHDF5 type(viirs_sdr) :: x integer :: ic, ichan, il, i, is,lig,pix CHARACTER(LEN=8) :: fic_viirsM(16)=(/'virrs_01', 'virrs_02', 'virrs_03', 'virrs_04',& 'virrs_05', 'virrs_06', 'virrs_07', 'virrs_08',& 'virrs_09', 'virrs_10', 'virrs_11', 'virrs_12',& 'virrs_13', 'virrs_14', 'virrs_15', 'virrs_16'/) CHARACTER(LEN=8) :: fic_viirsI(5)=(/'virrsI_1','virrsI_2','virrsI_3','virrsI_4','virrsI_5'/) integer :: dateg(3), iq(3) character(len=22) :: cht integer :: nchannels, npixels, nlines, nscans, ngranules real(kind=8) :: ij0, ij integer :: iscan, ifov, iline, idet integer :: err integer :: hour,minute real :: second logical :: has_I TRY !real, allocatable :: Obs(:,:) !integer, allocatable :: Qual_Obs(:,:), Qual_Scan(:) !----------------------------------------------------------------- ! initialise l interface fortran CALL OPEN_HDF( err ) THROWM(err.ne.0,"Unable to open HDF library") !----------------------------------------------------------------- ! ----------------------------------- ! number of Moderate resolution lines ! ----------------------------------- ! Load only requested channels call viirs_sdr_load('M', x, fic_viirsM(ic_Viirs_M (1:V_maxch_M)), err, & & channels=ic_Viirs_M (1:V_maxch_M),geolocfile='virrsloc', clean=.true.) THROWM(err.ne.0,"VIIRS SDR Load error") call viirs_sdr_info( x, nchannels, npixels, nlines, nscans, Ngranules, err ) THROWM(err.ne.0,"VIIRS SDR Info error") if(idbg%debug) then write(*,*) "nchannels ", nchannels write(*,*) "npixels ", npixels write(*,*) "nlines ", nlines write(*,*) "nscans ", nscans write(*,*) "ngranules ", field_M%nb_granules write(*,*) 'Distributor: ', x%root_att%Distributor write(*,*) 'Mission: ', x% root_att%Mission_Name endif ! allocate the viirs M fields with the number of lines ! -------------------------------------------------- allocate ( field_M % Obs(V_maxch_M,V_maxfovs,NLines), stat = err) ! viirs moderate resolution obs THROW(err.ne.0) allocate ( field_M % Qual_Obs(V_maxch_M,V_maxfovs,NLines), stat = err) THROW(err.ne.0) allocate ( field_M % lon(V_maxfovs,NLines), stat = err) ! longitude (degrees) THROW(err.ne.0) allocate ( field_M % lat(V_maxfovs,NLines), stat = err) ! latitude (degrees) THROW(err.ne.0) allocate ( field_M % satzen(V_maxfovs,NLines), stat = err) ! satellite zenithal angle(deg) THROW(err.ne.0) allocate ( field_M % satazi(V_maxfovs,NLines), stat = err) ! satellite azimutal angle (deg) THROW(err.ne.0) allocate ( field_M % sunzen(V_maxfovs,NLines), stat = err) ! solar zenithal angle(deg) THROW(err.ne.0) allocate ( field_M % sunazi(V_maxfovs,NLines), stat = err) ! solar azimutal angle(deg) THROW(err.ne.0) allocate ( field_M % difazi(V_maxfovs,NLines), stat = err) ! solar azimutal angle(deg) THROW(err.ne.0) allocate ( field_M % ij(NLines) , stat = err) ! time of the line THROW(err.ne.0) allocate ( field_M % locqual(V_maxfovs,NLines), stat=err) ! viirs fov geoloc quality THROW(err.ne.0) allocate ( field_M % scanqual(NLines), stat=err) ! viirs scan quality THROW(err.ne.0) field_M % difazi(:,:) = vmanq call maia_Viirs_field_init ( field_M ) field_M % N_Granule_ID = x%gran0_att(ic_Viirs_M (1))%N_Granule_ID field_M % Distributor = x%root_att%Distributor field_M % Nb_pixels = V_maxfovs ! number of pixels in line for Moderate resolution channels field_M % Nb_chan = V_maxch_M ! number of channels in field field_M % Nb_Lines = NLines ! number of lines in field for Moderate resolution channels field_M % NB_granules = Ngranules if ( x% root_att%Mission_Name .eq. "J01" ) then field_M % mission = "jpss1" elseif ( x% root_att%Mission_Name .eq. "J02" ) then field_M % mission = "jpss2" endif write(*,*) "mission :" ,field_M % mission ! loop on Moderate channels and fill the viirs fields ! --------------------------------------------------- field_M % ScanQual(:)=0 do ic=1, V_maxch_M ichan= ic_Viirs_M (ic) if(ichan .LE. 11) then field_M % Obs(ic,:,:) = x%band(ichan)%Reflectances(:,:)*100 else field_M % Obs(ic,:,:) = x%band(ichan)%BrightnessTemperature(:,:) endif where (field_M % Obs < -99) field_M % Obs = vmanq field_M % Qual_Obs(ic,:,:) = x%band(ichan)%QF1_VIIRSMBANDSDR(:,:) do is=1,nscans do il=(is-1)*MDet_InScan +1, is*MDet_InScan if(x%band(ichan)%QF2_SCAN_SDR(is)/=0) field_M % ScanQual(il) = 1 end do end do if( idbg % debug ) then write (*,*) ">>>>> ichan=",ichan do i=1,10 write (* ,*) " ligne ",i," ",field_M%obs(ic,1:10,i) enddo write (*,*) ">>>>> Qual_Obs" do i=1,10 write (* ,*) " ligne ",i," ",field_M%qual_obs(ic,1:10,i) enddo endif end do ! get the geolocation of the Moderate channels ! -------------------------------------------- dateg(1) = x%geo%StartYear dateg(2) = x%geo%StartMonth dateg(3) = x%geo%StartDay call xdgq(dateg, iq(1:2) ) read(x%agg_att%AggregateBeginningTime,"(i2,i2,f9.6)") hour, minute, second second=second+hour*3600.0 + minute*60.0 if ( idbg % debug ) then write(*,*) "beginning time secondes " , second write(*,*) "x%geo%UTC(1)",x%geo%UTC(1) endif iq(3)=second*1000 if( idbg % debug) write (0,*) "iq start : ",iq call xiqj(iq, field_M % ij_start) call xijc(field_M %ij_start, cht) if( idbg % debug) write(0,*) "UTC start of line 1 ",cht if( idbg % debug) write(0,*) "nscans= ",nscans if( idbg % debug) write (0,*) "NumberOfScans", x%geo%NumberOfScans if( idbg % debug) write (0,*) "x%geo%UTC(nscans)",x%geo%UTC(nscans) dateg(1) = x%geo%EndYear dateg(2) = x%geo%EndMonth dateg(3) = x%geo%EndDay call xdgq(dateg, iq(1:2) ) read(x%agg_att%AggregateEndingTime,"(i2,i2,f9.6)") hour, minute, second second=second+hour*3600.0 + minute*60.0 if( idbg % debug) write(*,*) "ending time " , second iq(3)=second*1000 if( idbg % debug) write (0,*) "iq end : ",iq call xiqj(iq, Field_M % ij_end) call xijc(field_M %ij_end, cht) if( idbg % debug) write(0,*) "UTC start of last line ",cht ! calculations for all lines if( idbg % debug) write(*,*)"MDet_InScan",MDet_InScan do iscan=1,NScans if (x%geo%UTC(iscan) > 0) then if( idbg % debug) write(*,*) "scan x%geo%UTCJulianDay =",iscan,x%geo%UTCJulianDay(iscan) ij = x%geo%UTCJulianDay(iscan) + x%geo%UTC(iscan)/86400.0d03 else ij = vmanq endif do idet=1,MDet_InScan iline = idet + (iscan-1)*MDet_InScan field_M % ij(iline) = ij! start time of the line end do end do do iscan=1,NScans do idet=1,MDet_InScan iline = idet + (iscan-1)*MDet_InScan if(x%geo%QF1_SCAN_VIIRSSDRGEO(iscan)>0 .and. & & x%geo%QF1_SCAN_VIIRSSDRGEO(iscan)<8) field_M % locqual(:,iline) = .false. end do do ifov= 1, npixels if(x%geo%QF2_VIIRSSDRGEO(ifov, iline) >0) field_M % locqual(ifov,iline) =.false. end do end do field_M % Orbit_num_start = x%agg_att%AggregateBeginningOrbitNumber ! Beginning orbit number field_M % Orbit_num_end = x%agg_att%AggregateEndingOrbitNumber ! Endding orbit number field_M % ModScan(:) = x%band(ic_Viirs_M (1))%ModeScan(:) ! Scan operational mode (0: processed) field_M % lat(:,:) = x%geo%Latitude (:,:) ! longitude (degrees) field_M % lon(:,:) = x%geo%Longitude (:,:) ! latitude (degrees) field_M % satzen(:,:) = x%geo%SatelliteZenithAngle (:,:) ! satellite zenithal angle(deg) field_M % satazi(:,:) = x%geo%SatelliteAzimuthAngle (:,:)! satellite azimutal angle (deg) field_M % sunzen(:,:) = x%geo%SolarZenithAngle (:,:) ! solar zenithal angle(deg) field_M % sunazi(:,:) = x%geo%SolarAzimuthAngle (:,:) ! solar azimutal angle(deg) if(idbg%debug) then write(*,'(a5,10f8.2)') 'lat', field_M % lat(1600,1:10) write(*,'(a5,10f8.2)') 'lon', field_M % lon(1600,1:10) endif call viirs_sdr_free(x,err) THROWM(err.ne.0,"Free SDR structure") ! ----------------------------------------------------------- ! number of lines for the Visible Imager resolution channels ! ----------------------------------------------------------- ichan=5 field_I % Nb_Lines = imanq ! number of lines in field for Imager resolution channels inquire(file='virrsloI', exist=has_I) if (has_I) then call viirs_sdr_load('I', x, fic_viirsI, err, geolocfile='virrsloI', clean=.true.) THROWM(err.ne.0,"Error SRD Load I Channels") call viirs_sdr_info( x, nchannels, npixels, nlines, nscans, Ngranules, err ) THROWM(err.ne.0,"Error SRD Info I Channels") else write(*,*) "WARNING : there are no I resolution files" endif !NLines = Max_Ilines_InGran * Nb_Granules if(idbg%debug) then write(*,*) 'NLines Imaging= ', NLines endif ! allocate the viirs I fields with the number of lines ! -------------------------------------------------- allocate ( field_I % Obs(5,V_maxfovs_I,NLines), stat = err) ! viirs moderate resolution obs THROW(err.ne.0) allocate ( field_I % Qual_Obs(5,V_maxfovs_I,NLines), stat = err) THROW(err.ne.0) allocate ( field_I % lon(V_maxfovs_I,NLines), stat = err) ! longitude (degrees) THROW(err.ne.0) allocate ( field_I % lat(V_maxfovs_I,NLines), stat = err) ! latitude (degrees) THROW(err.ne.0) allocate ( field_I % satzen(V_maxfovs_I,NLines), stat = err) ! satellite zenithal angle(deg) THROW(err.ne.0) allocate ( field_I % satazi(V_maxfovs_I,NLines), stat = err) ! satellite azimutal angle (deg) THROW(err.ne.0) allocate ( field_I % sunzen(V_maxfovs_I,NLines), stat = err) ! solar zenithal angle(deg) THROW(err.ne.0) allocate ( field_I % sunazi(V_maxfovs_I,NLines), stat = err) ! solar azimutal angle(deg) THROW(err.ne.0) allocate ( field_I % ij(NLines) , stat = err) ! time of the line THROW(err.ne.0) allocate ( field_I % locqual(V_maxfovs_I,NLines), stat=err) ! viirs fov geoloc quality THROW(err.ne.0) allocate ( field_I % scanqual(NLines), stat=err) ! viirs scan quality THROW(err.ne.0) call maia_Viirs_field_init ( field_I ) if (has_I) then field_I % N_Granule_ID = x%gran0_att(1)%N_Granule_ID field_I % Distributor = x%root_att%Distributor if ( x% root_att%Mission_Name .eq. "J01" ) then field_I % mission = "jpss1" elseif ( x% root_att%Mission_Name .eq. "J02" ) then field_I % mission = "jpss2" endif field_I % NB_granules = Ngranules field_I % Nb_Lines = NLines ! number of lines in field for Imager resolution channels field_I % Nb_pixels = V_maxfovs_I ! number of pixels in line for Imager resolution channels field_I % Nb_chan = 5 ! number of channels in field do ic=1, 5 ichan= ic if(ichan .LT. 4) then field_I % Obs(ic,:,:) = x%band(ichan)%Reflectances(:,:)*100 else field_I % Obs(ic,:,:) = x%band(ichan)%BrightnessTemperature(:,:) endif field_I % Qual_Obs(ic,:,:) = x%band(ichan)%QF1_VIIRSIBANDSDR(:,:) do is=1,nscans do il=(is-1)*IDet_InScan +1, is*IDet_InScan if(x%band(ichan)%QF2_SCAN_SDR(is)/=0) field_I % ScanQual(il) = 1 end do end do end do ! get the geolocation of the Imager channels ! -------------------------------------------- dateg(1) = x%geo%StartYear dateg(2) = x%geo%StartMonth dateg(3) = x%geo%StartDay call xdgq(dateg, iq(1:2) ) read(x%agg_att%AggregateBeginningTime,"(i2,i2,f9.6)") hour, minute, second second=second+hour*3600.0 + minute*60.0 if( idbg % debug) write(*,*) "beginning time " , second iq(3)=second*1000 call xiqj(iq, field_I % ij_start) call xijc(field_I %ij_start, cht) if( idbg % debug ) write(0,*) "UTC start of line 1 ",cht dateg(1) = x%geo%EndYear dateg(2) = x%geo%EndMonth dateg(3) = x%geo%EndDay call xdgq(dateg, iq(1:2) ) read(x%agg_att%AggregateEndingTime,"(i2,i2,f9.6)") hour, minute, second second=second+hour*3600.0 + minute*60.0 if( idbg % debug ) write(*,*) "ending time " , second iq(3)=second*1000 call xiqj(iq, field_I % ij_end) call xijc(field_I %ij_end, cht) if( idbg % debug ) write(0,*) "UTC start of last line ",cht iq(3) = 0 call xiqj(iq, ij0) !ngranules = size(x%band(channel)%ModeGran) !nscans = size(x%band(channel)%ModeScan) ! doit etre fait dans le module geo_proc do iscan=1,NScans if (x%geo%UTC(iscan) > 0) then ij = ij0 + x%geo%UTC(iscan)/86400.0d03 else ij = vmanq endif do idet=1,IDet_InScan iline = idet + (iscan-1)*IDet_InScan field_I % ij(iline) = ij! start time of the line end do end do do iscan=1,NScans do idet=1,IDet_InScan iline = idet + (iscan-1)*IDet_InScan if(x%geo%QF1_SCAN_VIIRSSDRGEO(iscan)>0 .and. & & x%geo%QF1_SCAN_VIIRSSDRGEO(iscan)<8) field_I % locqual(:,iline) = .false. end do do ifov= 1, npixels if(x%geo%QF2_VIIRSSDRGEO(ifov, iline) >0) field_I % locqual(ifov,iline) =.false. end do end do field_I % Orbit_num_start = x%agg_att%AggregateBeginningOrbitNumber ! Beginning orbit number field_I % Orbit_num_end = x%agg_att%AggregateEndingOrbitNumber ! Endding orbit number field_I % ModScan(:) = x%band(1)%ModeScan(:) ! Scan operational mode (0: processed) field_I % lat(:,:) = x%geo%Latitude (:,:) ! longitude (degrees) field_I % lon(:,:) = x%geo%Longitude (:,:) ! latitude (degrees) field_I % satzen(:,:) = x%geo%SatelliteZenithAngle (:,:) ! satellite zenithal angle(deg) field_I % satazi(:,:) = x%geo%SatelliteAzimuthAngle (:,:)! satellite azimutal angle (deg) field_I % sunzen(:,:) = x%geo%SolarZenithAngle (:,:) ! solar zenithal angle(deg) field_I % sunazi(:,:) = x%geo%SolarAzimuthAngle (:,:) ! solar azimutal angle(deg) call viirs_sdr_free(x,err) THROWM(err.ne.0,"Free SDR structure") CALL CLOSE_HDF( err ) THROWM(err.ne.0,"Close HDF library") ! missing data handling for Lydie read routine compatibility do lig=1, field_I % Nb_Lines do pix=1, field_I % Nb_pixels if (field_I % Obs(1,pix,lig) <0 ) field_I % Obs(1,pix,lig)=-999900. if (field_I % Obs(2,pix,lig) <0 ) field_I % Obs(2,pix,lig)=-999900. if (field_I % Obs(3,pix,lig) <0 ) field_I % Obs(3,pix,lig)=-999900. if (field_I % Obs(4,pix,lig) <0 ) field_I % Obs(4,pix,lig)=-9999. if (field_I % Obs(5,pix,lig) <0 ) field_I % Obs(5,pix,lig)=-9999. end do end do endif ! if (has_I) CATCH #endif end subroutine maia_read_Viirs subroutine maia_Viirs_field_init ( viirs) USE MOD_MAIA_CONST USE mod_maia_types IMPLICIT NONE Type (field), intent(inout) :: viirs viirs % mission = 'npp' ! default satellite name viirs % ij_start =vmanq ! Beginning time viirs % ij_end =vmanq ! Endding time viirs % Orbit_num_start =imanq ! Beginning orbit number viirs % Orbit_num_end =imanq ! Endding orbit number viirs % Nb_Granules =imanq ! number of granules in field viirs % Nb_pixels = imanq ! number of pixels in line for Moderate resolution channels viirs % Nb_chan = imanq ! number of channels in field viirs % lon(:,:) =vmanq ! longitude (degrees) viirs % lat(:,:) =vmanq ! latitude (degrees) viirs % satzen(:,:) =vmanq ! satellite zenithal angle(deg) viirs % satazi(:,:) =vmanq ! satellite azimutal angle (deg) viirs % sunzen(:,:) =vmanq ! solar zenithal angle(deg) viirs % sunazi(:,:) =vmanq ! solar azimutal angle(deg) viirs % Obs(:,:,:) = vmanq ! viirs moderate resolution obs viirs % Qual_Obs(:,:,:) =imanq viirs % ij(:) =vmanq ! time of the line viirs % locqual(:,:) =.true. end subroutine maia_Viirs_field_init !******************************************************************************