! $Id: ropp_io_read_ncdf_get.f90 6490 2020-09-10 17:36:27Z idculv $ !****is* Reading/ropp_io_read_ncdf_get * ! ! NAME ! ropp_io_read_ncdf_get - Get data from an (already defined) netCDF. ! ! SYNOPSIS ! call ropp_io_read_ncdf_get(data, rec) ! ! DESCRIPTION ! This subroutine gets variables from an already openend / netCDF data ! file. ! Reads core parameters from V1.0 format version and new parameters ! added in V1.1. Additional variables contained in the netCDF file are ! read in and define elements of the Vlist structures in ROprof. ! ! NOTES ! A netCDF file must have been opened before; this subroutine only works ! on the current netcdf file. The netCDF data file is left open. ! ! AUTHOR ! Met Office, Exeter, UK. ! Any comments on this software should be given via the ROM SAF ! Helpdesk at http://www.romsaf.org ! ! COPYRIGHT ! (c) EUMETSAT. All rights reserved. ! For further details please refer to the file COPYRIGHT ! which you should have received as part of this distribution. ! !**** !------------------------------------------------------------------------------- ! 1. Core RO data !------------------------------------------------------------------------------- SUBROUTINE ropp_io_read_ncdf_get_rodata(DATA, rec) ! 1.1 Declarations ! ---------------- USE ropp_utils USE ncdf USE ropp_io, not_this => ropp_io_read_ncdf_get_rodata USE ropp_io_types, ONLY: ROprof, & ThisFmtVer, & PCD_summary, & PCD_bangle, & PCD_refrac IMPLICIT NONE TYPE(ROprof), INTENT(inout) :: DATA INTEGER, OPTIONAL :: rec INTEGER :: n, n1, n2, ierr INTEGER :: irec INTEGER, DIMENSION(8) :: DT8 ! Date/time array CHARACTER(len = 23) :: proc_date CHARACTER(LEN=15) :: cval REAL(dp) :: time, start_time, dtocc_time, now_time !dp defined in DateTimeTypes REAL :: fmtver INTEGER :: status, varid, ndim, TYPE CHARACTER(len = 256) :: routine CHARACTER(len = 256) :: readstr ! 1.2 Error handling ! ------------------ CALL message_get_routine(routine) CALL message_set_routine('ropp_io_read_ncdf_get') ! 1.3 Default parameters ! ---------------------- IF (PRESENT(rec)) THEN irec = rec ELSE irec = 1 ENDIF ! 1.4 (Global) Attributes ! ------------------------ data%FmtVersion = ' ' ; CALL ncdf_getatt('format_version', data%FmtVersion) READ ( data%FmtVersion(11:), fmt=*, iostat=ierr ) fmtver IF ( ierr /= 0 ) data%FmtVersion = ThisFmtVer data%processing_centre = ' ' ; CALL ncdf_getatt('processing_centre', data%processing_centre) ! print*,'data%processing_centre (from NetCDF): ',data%processing_centre ! data%processing_centre = 'NESDIS UMD/ESSIC' ! print*,'data%processing_centre (after modified for BUFR): ',data%processing_centre IF (ncdf_isatt('processing_software')) THEN ! added at V8.0 data%processing_software = ' ' ; CALL ncdf_getatt('processing_software', data%processing_software) ENDIF proc_date = ' ' ; CALL ncdf_getatt('processing_date', proc_date) data%pod_method = ' ' ; CALL ncdf_getatt('pod_method', data%pod_method) data%phase_method = ' ' ; CALL ncdf_getatt('phase_method', data%phase_method) data%bangle_method = ' ' ; CALL ncdf_getatt('bangle_method', data%bangle_method) data%refrac_method = ' ' ; CALL ncdf_getatt('refrac_method', data%refrac_method) data%meteo_method = ' ' ; CALL ncdf_getatt('meteo_method', data%meteo_method) IF(ncdf_isatt('thin_method'))THEN ! added at V1.1 data%thin_method = ' ' ; CALL ncdf_getatt('thin_method', data%thin_method) ENDIF data%software_version = ' ' ; CALL ncdf_getatt('software_version', data%software_version) data%bad = ' ' ; CALL ncdf_getatt('bad', data%bad) IF (proc_date( 1: 4) /= ' ') READ(proc_date( 1: 4), *) data%DTpro%Year IF (proc_date( 6: 7) /= ' ') READ(proc_date( 6: 7), *) data%DTpro%Month IF (proc_date( 9:10) /= ' ') READ(proc_date( 9:10), *) data%DTpro%Day IF (proc_date(12:13) /= ' ') READ(proc_date(12:13), *) data%DTpro%Hour IF (proc_date(15:16) /= ' ') READ(proc_date(15:16), *) data%DTpro%Minute IF (proc_date(18:19) /= ' ') READ(proc_date(18:19), *) data%DTpro%Second IF (proc_date(21:23) /= ' ') READ(proc_date(21:23), *) data%DTpro%Msec ! 1.5 Header variables ! -------------------- CALL ncdf_getvar('occ_id', data%occ_id, rec = irec) CALL ncdf_getvar('gns_id', data%gns_id, rec = irec) CALL ncdf_getvar('leo_id', data%leo_id, rec = irec) CALL ncdf_getvar('stn_id', data%stn_id, rec = irec) CALL ncdf_getvar('occfreq1', data%occfreq1, rec = irec) CALL ncdf_getvar('occfreq2', data%occfreq2, rec = irec) IF ( data%leo_id(1:2) == 'CH' ) THEN data%leo_id = 'CHMP' ELSE IF ( data%leo_id(1:4) == 'GRC1' ) THEN data%leo_id = 'GRAA' ELSE IF ( data%leo_id(1:4) == 'GRC2' ) THEN data%leo_id = 'GRAB' ELSE IF ( data%leo_id(1:2) == 'CO' ) THEN data%leo_id(1:2) = 'C0' ELSE IF ( data%leo_id(1:3) == 'MTP' ) THEN data%leo_id(1:3) = 'MET' ELSE IF ( data%leo_id(1:4) == 'G001' ) THEN ! shouldn't ever be needed data%leo_id = 'FY3C' ELSE IF ( data%leo_id(1:4) == 'G002' ) THEN ! shouldn't ever be needed data%leo_id = 'FY3D' END IF ! 1.6 Date and time ! ----------------- CALL ncdf_getvar('start_time', start_time, rec=irec) CALL ncdf_getvar('year', data%DTocc%Year, & units = data%DTocc%units%Year, & range = data%DTocc%range%Year, & rec = irec) CALL ncdf_getvar('month', data%DTocc%Month, & units = data%DTocc%units%Month, & range = data%DTocc%range%Month, & rec = irec) CALL ncdf_getvar('day', data%DTocc%Day, & units = data%DTocc%units%Day, & range = data%DTocc%range%Day, & rec = irec) CALL ncdf_getvar('hour', data%DTocc%Hour, & units = data%DTocc%units%Hour, & range = data%DTocc%range%Hour, & rec = irec) CALL ncdf_getvar('minute', data%DTocc%Minute, & units = data%DTocc%units%Minute, & range = data%DTocc%range%Minute, & rec = irec) CALL ncdf_getvar('second', data%DTocc%Second, & units = data%DTocc%units%Second, & range = data%DTocc%range%Second, & rec = irec) CALL ncdf_getvar('msec', data%DTocc%Msec, & units = data%DTocc%units%Msec, & range = data%DTocc%range%Msec, & rec = irec) ! 1.6.1 Check consistency: start_time and DTocc (if both are valid) ! should refer to the same epoch within 1ms. Issue a warning if not. CALL Date_and_Time_UTC(Values=DT8) CALL TimeSince(DT8, now_time, 1, Base="JS2000") IF (isroppinrange(data%DTocc)) THEN IF ( isinrange(start_time, (/ 0.001_dp, now_time /)) ) THEN DT8 = (/ data%DTocc%Year, data%DTocc%Month, data%DTocc%Day, 0, & data%DTocc%Hour, data%DTocc%Minute, data%DTocc%Second, & data%DTocc%Msec /) CALL TimeSince(DT8, dtocc_time, 1, Base="JS2000") dtocc_time = dtocc_time + (LeapSeconds(DT8) - 22)*1.D0 ! 22 leapseconds had elapsed by 0Z 01/01/2000 time = ABS(start_time - dtocc_time) IF ( time > 0.0005_dp ) THEN WRITE ( cval, FMT="(F15.3)" ) time IF ( time < 30.0_dp ) THEN ! Probably a result of not accounting for leap-seconds, so reduce severity of message to DIAG CALL message(msg_diag, "'start_time' and yr/mo/dy/hr/mn/sc/ms " // & "timestamps differ by " // & TRIM(ADJUSTL(cval))//" seconds (probably a leap-second issue) - using yr/../ms timestamp") ELSE CALL message(msg_warn, "'start_time' and yr/mo/dy/hr/mn/sc/ms " // & "timestamps differ by " // & TRIM(ADJUSTL(cval))//" seconds - using yr/../ms timestamp") END IF END IF END IF ! If any DTocc element is invalid, substitute converted start_time ! (if that is valid) and issue a warning. ELSE IF ( isinrange(start_time, (/ 0.001_dp, now_time /)) ) THEN CALL message(msg_warn, "One or more of yr/mo/dy/hr/mn/sc/ms times " // & "are invalid - using 'start_time' instead") CALL TimeSince(DT8, start_time, -1, Base="JS2000") start_time = start_time - (LeapSeconds(DT8) - 22)*1.D0 ! 22 leapseconds had elapsed by 0Z 01/01/2000 CALL TimeSince(DT8, start_time, -1, Base="JS2000") data%DTocc%Year = DT8(1) data%DTocc%Month = DT8(2) data%DTocc%Day = DT8(3) data%DTocc%Hour = DT8(5) data%DTocc%Minute = DT8(6) data%DTocc%Second = DT8(7) data%DTocc%Msec = DT8(8) END IF ! NB if neither DTocc nor start_time are valid, do nothing here; missing ! coordinates should be dealt with as part of generic Q/C. END IF ! 1.7 Overall quality ! ------------------- CALL ncdf_getvar('pcd', data%pcd, & units = data%units%pcd, & range = data%range%pcd, & rec = irec) CALL ncdf_getvar('overall_qual', data%overall_qual, & units = data%units%overall_qual, & range = data%range%overall_qual, & rec = irec) ! Yong Chen on 05/24/2024, to reset pcd and overall qual based on ! the postprocess quality flag 'bad' for RFSI and ROPP IF(ncdf_isatt('bad')) THEN readstr = ' ' CALL ncdf_getatt('bad', readstr) IF ( TRIM(readstr(1:1)) /= "0" ) THEN data%pcd = IBSET(data%pcd, PCD_summary) ! PCD bit 1 set: non-nominal quality data%pcd = IBSET(data%pcd, PCD_bangle) ! PCD bit 5 set: non-nominal bangle data%pcd = IBSET(data%pcd, PCD_refrac) ! PCD bit 6 set: non-nominal refrac data%overall_qual = 0.0 ELSE data%overall_qual = 100.0 ENDIF ENDIF ! 1.8 Georeferencing ! ------------------ CALL ncdf_getvar('time', time, rec=irec) CALL ncdf_getvar('lat', data%georef%lat, & units = data%georef%units%lat, & range = data%georef%range%lat, & rec = irec) CALL ncdf_getvar('lon', data%georef%lon, & units = data%georef%units%lon, & range = data%georef%range%lon, & rec = irec) CALL ncdf_getvar('time_offset', data%georef%time_offset, & units = data%georef%units%time_offset, & range = data%georef%range%time_offset, & rec = irec) CALL ncdf_getvar('undulation', data%georef%Undulation, & units = data%georef%units%Undulation, & range = data%georef%range%undulation, & rec = irec) ! Yong Chen on 05/30/2024 set Undulation to zero if not in the range !write(*, *) 'Undulation range is: ', data%georef%range%undulation IF ( data%georef%Undulation < data%georef%range%undulation(1) .or. & data%georef%Undulation > data%georef%range%undulation(2) ) & data%georef%Undulation = 0.0 CALL ncdf_getvar('roc', data%georef%roc, & units = data%georef%units%roc, & range = data%georef%range%roc, & rec = irec) CALL ncdf_getvar('r_coc', data%georef%r_coc, & units = data%georef%units%r_coc, & range = data%georef%range%r_coc, & rec = irec) CALL ncdf_getvar('r_alt', data%georef%r_alt, & units = data%georef%units%r_alt, & range = data%georef%range%r_alt, & rec = irec) CALL ncdf_getvar('viewang', data%georef%viewang, & units = data%georef%units%viewang, & range = data%georef%range%viewang, & rec = irec) CALL ncdf_getvar('azimuth', data%georef%azimuth, & units = data%georef%units%azimuth, & range = data%georef%range%azimuth, & rec = irec) IF ( ncdf_isvar('r_leo_ref') ) & CALL ncdf_getvar('r_leo_ref', data%georef%leo_pod%pos, & units = data%georef%units%r_leo, & range = data%georef%range%r_leo, & rec = irec) IF ( ncdf_isvar('v_leo_ref') ) & CALL ncdf_getvar('v_leo_ref', data%georef%leo_pod%vel, & units = data%georef%units%v_leo, & range = data%georef%range%v_leo, & rec = irec) IF ( ncdf_isvar('r_gns_ref') ) & CALL ncdf_getvar('r_gns_ref', data%georef%gns_pod%pos, & units = data%georef%units%r_gns, & range = data%georef%range%r_gns, & rec = irec) IF ( ncdf_isvar('v_gns_ref') ) & CALL ncdf_getvar('v_gns_ref', data%georef%gns_pod%vel, & units = data%georef%units%v_gns, & range = data%georef%range%v_gns, & rec = irec) ! 1.8.1 Other attributes CALL ncdf_getatt('reference_frame', data%georef%reference_frame%r_coc, varname= 'r_coc') IF ( ncdf_isvar('r_leo_ref') ) & CALL ncdf_getatt('reference_frame', data%georef%leo_pod%pos_frame, varname= 'r_leo_ref') IF ( ncdf_isvar('v_leo_ref') ) & CALL ncdf_getatt('reference_frame', data%georef%leo_pod%vel_frame, varname= 'v_leo_ref') IF ( ncdf_isvar('r_gns_ref') ) & CALL ncdf_getatt('reference_frame', data%georef%gns_pod%pos_frame, varname= 'r_gns_ref') IF ( ncdf_isvar('v_gns_ref') ) & CALL ncdf_getatt('reference_frame', data%georef%gns_pod%vel_frame, varname= 'v_gns_ref') ! 1.9 Background characterisation (if any) ! ---------------------------------------- IF (ncdf_isvar('bg_source')) THEN data%BG%Source = 'TBD' ELSE data%BG%Source = 'NONE' ENDIF IF (data%BG%Source /= 'NONE') THEN CALL ncdf_getvar('bg_source', data%BG%Source, rec = irec) CALL ncdf_getvar('bg_year', data%BG%Year, & units = data%BG%units%Year, & range = data%BG%range%Year, & rec = irec) CALL ncdf_getvar('bg_month', data%BG%Month, & units = data%BG%units%Month, & range = data%BG%range%Month, & rec = irec) CALL ncdf_getvar('bg_day', data%BG%Day, & units = data%BG%units%Day, & range = data%BG%range%Day, & rec = irec) CALL ncdf_getvar('bg_hour', data%BG%Hour, & units = data%BG%units%Hour, & range = data%BG%range%Hour, & rec = irec) CALL ncdf_getvar('bg_minute', data%BG%Minute, & units = data%BG%units%Minute, & range = data%BG%range%Minute, & rec = irec) CALL ncdf_getvar('bg_fcperiod', data%BG%fcPeriod, & units = data%BG%units%fcPeriod, & range = data%BG%range%fcPeriod, & rec = irec) ENDIF ! 1.9.1 QA characterisation (if any) ! ---------------------------------------- CALL ncdf_getvar('snr1avg', data%QA%snr1avg, & units = data%QA%units%snr1avg, & range = data%QA%range%snr1avg, & rec = irec) CALL ncdf_getvar('snr2avg', data%QA%snr2avg, & units = data%QA%units%snr2avg, & range = data%QA%range%snr2avg, & rec = irec) CALL ncdf_getvar('irs', data%QA%irs, & units = data%QA%units%irs, & range = data%QA%range%irs, & rec = irec) CALL ncdf_getvar('smean', data%QA%smean, & units = data%QA%units%smean, & range = data%QA%range%smean, & rec = irec) CALL ncdf_getvar('stdv', data%QA%stdv, & units = data%QA%units%stdv, & range = data%QA%range%stdv, & rec = irec) CALL ncdf_getvar('reldevmax', data%QA%reldevmax, & units = data%QA%units%reldevmax, & range = data%QA%range%reldevmax, & rec = irec) CALL ncdf_getvar('reldevstd', data%QA%reldevstd, & units = data%QA%units%reldevstd, & range = data%QA%range%reldevstd, & rec = irec) CALL ncdf_getvar('reldevmax1', data%QA%reldevmax1, & units = data%QA%units%reldevmax1, & range = data%QA%range%reldevmax1, & rec = irec) CALL ncdf_getvar('reldevstd1', data%QA%reldevstd1, & units = data%QA%units%reldevstd1, & range = data%QA%range%reldevstd1, & rec = irec) CALL ncdf_getvar('difmaxref', data%QA%difmaxref, & units = data%QA%units%difmaxref, & range = data%QA%range%difmaxref, & rec = irec) CALL ncdf_getvar('maxdifphase', data%QA%maxdifphase, & units = data%QA%units%maxdifphase, & range = data%QA%range%maxdifphase, & rec = irec) CALL ncdf_getvar('finalQA', data%QA%finalQA, & units = data%QA%units%finalQA, & range = data%QA%range%finalQA, & rec = irec) ! 1.10 Level1a variables (if any) ! ------------------------------ IF (ncdf_isvar('dtime')) THEN CALL ncdf_getsize('dtime', n, dim = 1) CALL ropp_io_init(data%Lev1a, n) ELSE data%Lev1a%Npoints = 0 ENDIF WRITE(*, *) 'data%Lev1a%Npoints before reading Level1a', data%Lev1a%Npoints IF (data%Lev1a%Npoints > 0) THEN CALL ncdf_getvar('dtime', data%Lev1a%dtime, & units = data%Lev1a%units%dtime, & range = data%Lev1a%range%dtime, & rec = irec) CALL ncdf_getvar('snr_L1ca', data%Lev1a%snr_L1ca, & units = data%Lev1a%units%snr, & range = data%Lev1a%range%snr, & rec = irec) CALL ncdf_getvar('snr_L1p', data%Lev1a%snr_L1p, & units = data%Lev1a%units%snr, & range = data%Lev1a%range%snr, & rec = irec) CALL ncdf_getvar('snr_L2p', data%Lev1a%snr_L2p, & units = data%Lev1a%units%snr, & range = data%Lev1a%range%snr, & rec = irec) CALL ncdf_getvar('phase_L1', data%Lev1a%phase_L1, & units = data%Lev1a%units%phase, & range = data%Lev1a%range%phase, & rec = irec) CALL ncdf_getvar('phase_L2', data%Lev1a%phase_L2, & units = data%Lev1a%units%phase, & range = data%Lev1a%range%phase, & rec = irec) CALL ncdf_getvar('r_gns', data%Lev1a%r_gns, & units = data%Lev1a%units%r_gns, & range = data%Lev1a%range%r_gns, & rec = irec) CALL ncdf_getvar('v_gns', data%Lev1a%v_gns, & units = data%Lev1a%units%v_gns, & range = data%Lev1a%range%v_gns, & rec = irec) CALL ncdf_getvar('r_leo', data%Lev1a%r_leo, & units = data%Lev1a%units%r_leo, & range = data%Lev1a%range%r_leo, & rec = irec) CALL ncdf_getvar('v_leo', data%Lev1a%v_leo, & units = data%Lev1a%units%v_leo, & range = data%Lev1a%range%v_leo, & rec = irec) CALL ncdf_getvar('phase_qual', data%Lev1a%phase_qual, & units = data%Lev1a%units%phase_qual, & range = data%Lev1a%range%phase_qual, & rec = irec) ! 1.10.1 Other attributes CALL ncdf_getatt('reference_frame', data%Lev1a%reference_frame%r_gns, varname = 'r_gns') CALL ncdf_getatt('reference_frame', data%Lev1a%reference_frame%v_gns, varname = 'v_gns') CALL ncdf_getatt('reference_frame', data%Lev1a%reference_frame%r_leo, varname = 'r_leo') CALL ncdf_getatt('reference_frame', data%Lev1a%reference_frame%v_leo, varname = 'v_leo') data%Lev1a%Missing = .FALSE. ENDIF ! 1.11 Level1b variables (if any) ! ------------------------------- IF (ncdf_isvar('lat_tp')) THEN CALL ncdf_getsize('lat_tp', n, dim = 1) CALL ropp_io_init(data%Lev1b, n) ELSE data%Lev1b%Npoints = 0 ENDIF IF (data%Lev1b%Npoints > 0) THEN CALL ncdf_getvar('lat_tp', data%Lev1b%lat_tp, & units = data%Lev1b%units%lat_tp, & range = data%Lev1b%range%lat_tp, & rec = irec) CALL ncdf_getvar('lon_tp', data%Lev1b%lon_tp, & units = data%Lev1b%units%lon_tp, & range = data%Lev1b%range%lon_tp, & rec = irec) CALL ncdf_getvar('azimuth_tp', data%Lev1b%azimuth_tp, & units = data%Lev1b%units%azimuth_tp, & range = data%Lev1b%range%azimuth_tp, & rec = irec) CALL ncdf_getvar('impact_L1', data%Lev1b%impact_L1, & units = data%Lev1b%units%impact, & range = data%Lev1b%range%impact, & rec = irec) CALL ncdf_getvar('impact_L2', data%Lev1b%impact_L2, & units = data%Lev1b%units%impact, & range = data%Lev1b%range%impact, & rec = irec) CALL ncdf_getvar('impact', data%Lev1b%impact, & units = data%Lev1b%units%impact, & range = data%Lev1b%range%impact, & rec = irec) ! Yong Chen on 05/24/2024, replace impact with impact_opt IF (ncdf_isvar('impact_opt')) & ! added at v1.1 CALL ncdf_getvar('impact_opt', data%Lev1b%impact_opt, & units = data%Lev1b%units%impact, & range = data%Lev1b%range%impact, & rec = irec) CALL ncdf_getvar('bangle_L1', data%Lev1b%bangle_L1, & units = data%Lev1b%units%bangle, & range = data%Lev1b%range%bangle, & rec = irec) CALL ncdf_getvar('bangle_L2', data%Lev1b%bangle_L2, & units = data%Lev1b%units%bangle, & range = data%Lev1b%range%bangle, & rec = irec) CALL ncdf_getvar('bangle', data%Lev1b%bangle, & units = data%Lev1b%units%bangle, & range = data%Lev1b%range%bangle, & rec = irec) ! Yong Chen on 05/24/2024, replace bangle with bangle_opt IF (ncdf_isvar('bangle_opt')) & ! added at v1.1 CALL ncdf_getvar('bangle_opt', data%Lev1b%bangle_opt, & units = data%Lev1b%units%bangle, & range = data%Lev1b%range%bangle, & rec = irec) CALL ncdf_getvar('bangle_L1_sigma', data%Lev1b%bangle_L1_sigma, & units = data%Lev1b%units%bangle_sigma, & range = data%Lev1b%range%bangle_sigma, & rec = irec) CALL ncdf_getvar('bangle_L2_sigma', data%Lev1b%bangle_L2_sigma, & units = data%Lev1b%units%bangle_sigma, & range = data%Lev1b%range%bangle_sigma, & rec = irec) CALL ncdf_getvar('bangle_sigma', data%Lev1b%bangle_sigma, & units = data%Lev1b%units%bangle_sigma, & range = data%Lev1b%range%bangle_sigma, & rec = irec) ! Yong Chen on 05/24/2024, replace bangle_sigma with bangle_L1_sigma IF (ncdf_isvar('bangle_opt_sigma')) & ! added at v1.1 CALL ncdf_getvar('bangle_opt_sigma', data%Lev1b%bangle_opt_sigma, & units = data%Lev1b%units%bangle_sigma, & range = data%Lev1b%range%bangle_sigma, & rec = irec) CALL ncdf_getvar('bangle_L1_qual', data%Lev1b%bangle_L1_qual, & units = data%Lev1b%units%bangle_qual, & range = data%Lev1b%range%bangle_qual, & rec = irec) CALL ncdf_getvar('bangle_L2_qual', data%Lev1b%bangle_L2_qual, & units = data%Lev1b%units%bangle_qual, & range = data%Lev1b%range%bangle_qual, & rec = irec) CALL ncdf_getvar('bangle_qual', data%Lev1b%bangle_qual, & units = data%Lev1b%units%bangle_qual, & range = data%Lev1b%range%bangle_qual, & rec = irec) IF (ncdf_isvar('bangle_opt_qual')) & ! added at v1.1 CALL ncdf_getvar('bangle_opt_qual', data%Lev1b%bangle_opt_qual, & units = data%Lev1b%units%bangle_qual, & range = data%Lev1b%range%bangle_qual, & rec = irec) data%Lev1b%Missing = .FALSE. ENDIF ! 1.12 Level2a variables (if any) ! ------------------------------- IF (ncdf_isvar('alt_refrac')) THEN CALL ncdf_getsize('alt_refrac', n, dim = 1) CALL ropp_io_init(data%Lev2a, n) ELSE data%Lev2a%Npoints = 0 ENDIF IF (data%Lev2a%Npoints > 0) THEN CALL ncdf_getvar('alt_refrac', data%Lev2a%alt_refrac, & units = data%Lev2a%units%alt_refrac, & range = data%Lev2a%range%alt_refrac, & rec = irec) CALL ncdf_getvar('geop_refrac', data%Lev2a%geop_refrac, & units = data%Lev2a%units%geop_refrac, & range = data%Lev2a%range%geop_refrac, & rec = irec) CALL ncdf_getvar('refrac', data%Lev2a%refrac, & units = data%Lev2a%units%refrac, & range = data%Lev2a%range%refrac, & rec = irec) CALL ncdf_getvar('refrac_sigma', data%Lev2a%refrac_sigma, & units = data%Lev2a%units%refrac_sigma, & range = data%Lev2a%range%refrac_sigma, & rec = irec) CALL ncdf_getvar('refrac_qual', data%Lev2a%refrac_qual, & units = data%Lev2a%units%refrac_qual, & range = data%Lev2a%range%refrac_qual, & rec = irec) IF (ncdf_isvar('dry_temp')) THEN !For backward compatibility CALL ncdf_getvar('dry_temp', data%Lev2a%dry_temp, & units = data%Lev2a%units%dry_temp, & range = data%Lev2a%range%dry_temp, & rec = irec) CALL ncdf_getvar('dry_temp_sigma', data%Lev2a%dry_temp_sigma, & units = data%Lev2a%units%dry_temp_sigma, & range = data%Lev2a%range%dry_temp_sigma, & rec = irec) CALL ncdf_getvar('dry_temp_qual', data%Lev2a%dry_temp_qual, & units = data%Lev2a%units%dry_temp_qual, & range = data%Lev2a%range%dry_temp_qual, & rec = irec) ENDIF data%Lev2a%Missing = .FALSE. ENDIF ! 1.13 Level2b variables (if any) ! ------------------------------- IF (ncdf_isvar('geop')) THEN CALL ncdf_getsize('geop', n, dim = 1) CALL ropp_io_init(data%Lev2b, n) ELSE data%Lev2b%Npoints = 0 ENDIF IF (data%Lev2b%Npoints > 0) THEN CALL ncdf_getvar('geop', data%Lev2b%geop, & units = data%Lev2b%units%geop, & range = data%Lev2b%range%geop, & rec = irec) CALL ncdf_getvar('geop_sigma', data%Lev2b%geop_sigma, & units = data%Lev2b%units%geop_sigma, & range = data%Lev2b%range%geop_sigma, & rec = irec) CALL ncdf_getvar('press', data%Lev2b%press, & units = data%Lev2b%units%press, & range = data%Lev2b%range%press, & rec = irec) CALL ncdf_getvar('press_sigma', data%Lev2b%press_sigma, & units = data%Lev2b%units%press_sigma, & range = data%Lev2b%range%press_sigma, & rec = irec) CALL ncdf_getvar('temp', data%Lev2b%temp, & units = data%Lev2b%units%temp, & range = data%Lev2b%range%temp, & rec = irec) CALL ncdf_getvar('temp_sigma', data%Lev2b%temp_sigma, & units = data%Lev2b%units%temp_sigma, & range = data%Lev2b%range%temp_sigma, & rec = irec) CALL ncdf_getvar('shum', data%Lev2b%shum, & units = data%Lev2b%units%shum, & range = data%Lev2b%range%shum, & rec = irec) CALL ncdf_getvar('shum_sigma', data%Lev2b%shum_sigma, & units = data%Lev2b%units%shum_sigma, & range = data%Lev2b%range%shum_sigma, & rec = irec) CALL ncdf_getvar('meteo_qual', data%Lev2b%meteo_qual, & units = data%Lev2b%units%meteo_qual, & range = data%Lev2b%range%meteo_qual, & rec = irec) data%Lev2b%Missing = .FALSE. ENDIF ! 1.14 Level2c variables (if any) ! ------------------------------- IF (ncdf_isvar('geop_sfc')) THEN data%Lev2c%Npoints = 1 ELSE data%Lev2c%Npoints = 0 ENDIF IF (data%Lev2c%Npoints > 0) THEN CALL ncdf_getvar('geop_sfc', data%Lev2c%geop_sfc, & units = data%Lev2c%units%geop_sfc, & range = data%Lev2c%range%geop_sfc, & rec = irec) CALL ncdf_getvar('press_sfc', data%Lev2c%press_sfc, & units = data%Lev2c%units%press_sfc, & range = data%Lev2c%range%press_sfc, & rec = irec) CALL ncdf_getvar('press_sfc_sigma', data%Lev2c%press_sfc_sigma, & units = data%Lev2c%units%press_sfc_sigma, & range = data%Lev2c%range%press_sfc_sigma, & rec = irec) CALL ncdf_getvar('press_sfc_qual', data%Lev2c%press_sfc_qual, & units = data%Lev2c%units%press_sfc_qual, & range = data%Lev2c%range%press_sfc_qual, & rec = irec) IF (ncdf_isvar('Ne_max')) THEN !For backward compatibility CALL ncdf_getvar('Ne_max', data%Lev2c%Ne_max, & units = data%Lev2c%units%Ne_max, & range = data%Lev2c%range%Ne_max, & rec = irec) ENDIF IF (ncdf_isvar('Ne_max_sigma')) THEN !For backward compatibility CALL ncdf_getvar('Ne_max_sigma', data%Lev2c%Ne_max_sigma, & units = data%Lev2c%units%Ne_max_sigma, & range = data%Lev2c%range%Ne_max_sigma, & rec = irec) ENDIF IF (ncdf_isvar('H_peak')) THEN !For backward compatibility CALL ncdf_getvar('H_peak', data%Lev2c%H_peak, & units = data%Lev2c%units%H_peak, & range = data%Lev2c%range%H_peak, & rec = irec) ENDIF IF (ncdf_isvar('H_peak_sigma')) THEN !For backward compatibility CALL ncdf_getvar('H_peak_sigma', data%Lev2c%H_peak_sigma, & units = data%Lev2c%units%H_peak_sigma, & range = data%Lev2c%range%H_peak_sigma, & rec = irec) ENDIF IF (ncdf_isvar('H_width')) THEN !For backward compatibility CALL ncdf_getvar('H_width', data%Lev2c%H_width, & units = data%Lev2c%units%H_width, & range = data%Lev2c%range%H_width, & rec = irec) ENDIF IF (ncdf_isvar('H_width_sigma')) THEN !For backward compatibility CALL ncdf_getvar('H_width_sigma', data%Lev2c%H_width_sigma, & units = data%Lev2c%units%H_width_sigma, & range = data%Lev2c%range%H_width_sigma, & rec = irec) ENDIF IF (ncdf_isvar('tph_bangle')) THEN !For backward compatibility CALL ncdf_getvar('tph_bangle', data%Lev2c%tph_bangle, & units = data%Lev2c%units%tph_bangle, & range = data%Lev2c%range%tph_bangle, & rec = irec) ENDIF IF (ncdf_isvar('tpa_bangle')) THEN !For backward compatibility CALL ncdf_getvar('tpa_bangle', data%Lev2c%tpa_bangle, & units = data%Lev2c%units%tpa_bangle, & range = data%Lev2c%range%tpa_bangle, & rec = irec) ENDIF IF (ncdf_isvar('tph_bangle_flag')) THEN !For backward compatibility CALL ncdf_getvar('tph_bangle_flag', data%Lev2c%tph_bangle_flag, & units = data%Lev2c%units%tph_bangle_flag, & range = data%Lev2c%range%tph_bangle_flag, & rec = irec) ENDIF IF (ncdf_isvar('tph_refrac')) THEN !For backward compatibility CALL ncdf_getvar('tph_refrac', data%Lev2c%tph_refrac, & units = data%Lev2c%units%tph_refrac, & range = data%Lev2c%range%tph_refrac, & rec = irec) ENDIF IF (ncdf_isvar('tpn_refrac')) THEN !For backward compatibility CALL ncdf_getvar('tpn_refrac', data%Lev2c%tpn_refrac, & units = data%Lev2c%units%tpn_refrac, & range = data%Lev2c%range%tpn_refrac, & rec = irec) ENDIF IF (ncdf_isvar('tph_refrac_flag')) THEN !For backward compatibility CALL ncdf_getvar('tph_refrac_flag', data%Lev2c%tph_refrac_flag, & units = data%Lev2c%units%tph_refrac_flag, & range = data%Lev2c%range%tph_refrac_flag, & rec = irec) ENDIF IF (ncdf_isvar('tph_tdry_lrt')) THEN !For backward compatibility CALL ncdf_getvar('tph_tdry_lrt', data%Lev2c%tph_tdry_lrt, & units = data%Lev2c%units%tph_tdry_lrt, & range = data%Lev2c%range%tph_tdry_lrt, & rec = irec) ENDIF IF (ncdf_isvar('tpt_tdry_lrt')) THEN !For backward compatibility CALL ncdf_getvar('tpt_tdry_lrt', data%Lev2c%tpt_tdry_lrt, & units = data%Lev2c%units%tpt_tdry_lrt, & range = data%Lev2c%range%tpt_tdry_lrt, & rec = irec) ENDIF IF (ncdf_isvar('tph_tdry_lrt_flag')) THEN !For backward compatibility CALL ncdf_getvar('tph_tdry_lrt_flag', data%Lev2c%tph_tdry_lrt_flag, & units = data%Lev2c%units%tph_tdry_lrt_flag, & range = data%Lev2c%range%tph_tdry_lrt_flag, & rec = irec) ENDIF IF (ncdf_isvar('tph_tdry_cpt')) THEN !For backward compatibility CALL ncdf_getvar('tph_tdry_cpt', data%Lev2c%tph_tdry_cpt, & units = data%Lev2c%units%tph_tdry_cpt, & range = data%Lev2c%range%tph_tdry_cpt, & rec = irec) ENDIF IF (ncdf_isvar('tpt_tdry_cpt')) THEN !For backward compatibility CALL ncdf_getvar('tpt_tdry_cpt', data%Lev2c%tpt_tdry_cpt, & units = data%Lev2c%units%tpt_tdry_cpt, & range = data%Lev2c%range%tpt_tdry_cpt, & rec = irec) ENDIF IF (ncdf_isvar('tph_tdry_cpt_flag')) THEN !For backward compatibility CALL ncdf_getvar('tph_tdry_cpt_flag', data%Lev2c%tph_tdry_cpt_flag, & units = data%Lev2c%units%tph_tdry_cpt_flag, & range = data%Lev2c%range%tph_tdry_cpt_flag, & rec = irec) ENDIF IF (ncdf_isvar('prh_tdry_cpt')) THEN !For backward compatibility CALL ncdf_getvar('prh_tdry_cpt', data%Lev2c%prh_tdry_cpt, & units = data%Lev2c%units%prh_tdry_cpt, & range = data%Lev2c%range%prh_tdry_cpt, & rec = irec) ENDIF IF (ncdf_isvar('prt_tdry_cpt')) THEN !For backward compatibility CALL ncdf_getvar('prt_tdry_cpt', data%Lev2c%prt_tdry_cpt, & units = data%Lev2c%units%prt_tdry_cpt, & range = data%Lev2c%range%prt_tdry_cpt, & rec = irec) ENDIF IF (ncdf_isvar('prh_tdry_cpt_flag')) THEN !For backward compatibility CALL ncdf_getvar('prh_tdry_cpt_flag', data%Lev2c%prh_tdry_cpt_flag, & units = data%Lev2c%units%prh_tdry_cpt_flag, & range = data%Lev2c%range%prh_tdry_cpt_flag, & rec = irec) ENDIF IF (ncdf_isvar('tph_temp_lrt')) THEN !For backward compatibility CALL ncdf_getvar('tph_temp_lrt', data%Lev2c%tph_temp_lrt, & units = data%Lev2c%units%tph_temp_lrt, & range = data%Lev2c%range%tph_temp_lrt, & rec = irec) ENDIF IF (ncdf_isvar('tpt_temp_lrt')) THEN !For backward compatibility CALL ncdf_getvar('tpt_temp_lrt', data%Lev2c%tpt_temp_lrt, & units = data%Lev2c%units%tpt_temp_lrt, & range = data%Lev2c%range%tpt_temp_lrt, & rec = irec) ENDIF IF (ncdf_isvar('tph_temp_lrt_flag')) THEN !For backward compatibility CALL ncdf_getvar('tph_temp_lrt_flag', data%Lev2c%tph_temp_lrt_flag, & units = data%Lev2c%units%tph_temp_lrt_flag, & range = data%Lev2c%range%tph_temp_lrt_flag, & rec = irec) ENDIF IF (ncdf_isvar('tph_temp_cpt')) THEN !For backward compatibility CALL ncdf_getvar('tph_temp_cpt', data%Lev2c%tph_temp_cpt, & units = data%Lev2c%units%tph_temp_cpt, & range = data%Lev2c%range%tph_temp_cpt, & rec = irec) ENDIF IF (ncdf_isvar('tpt_temp_cpt')) THEN !For backward compatibility CALL ncdf_getvar('tpt_temp_cpt', data%Lev2c%tpt_temp_cpt, & units = data%Lev2c%units%tpt_temp_cpt, & range = data%Lev2c%range%tpt_temp_cpt, & rec = irec) ENDIF IF (ncdf_isvar('tph_temp_cpt_flag')) THEN !For backward compatibility CALL ncdf_getvar('tph_temp_cpt_flag', data%Lev2c%tph_temp_cpt_flag, & units = data%Lev2c%units%tph_temp_cpt_flag, & range = data%Lev2c%range%tph_temp_cpt_flag, & rec = irec) ENDIF IF (ncdf_isvar('prh_temp_cpt')) THEN !For backward compatibility CALL ncdf_getvar('prh_temp_cpt', data%Lev2c%prh_temp_cpt, & units = data%Lev2c%units%prh_temp_cpt, & range = data%Lev2c%range%prh_temp_cpt, & rec = irec) ENDIF IF (ncdf_isvar('prt_temp_cpt')) THEN !For backward compatibility CALL ncdf_getvar('prt_temp_cpt', data%Lev2c%prt_temp_cpt, & units = data%Lev2c%units%prt_temp_cpt, & range = data%Lev2c%range%prt_temp_cpt, & rec = irec) ENDIF IF (ncdf_isvar('prh_temp_cpt_flag')) THEN !For backward compatibility CALL ncdf_getvar('prh_temp_cpt_flag', data%Lev2c%prh_temp_cpt_flag, & units = data%Lev2c%units%prh_temp_cpt_flag, & range = data%Lev2c%range%prh_temp_cpt_flag, & rec = irec) ENDIF IF (ncdf_isvar('pblh_bangle')) THEN !For backward compatibility CALL ncdf_getvar('pblh_bangle', data%Lev2c%pblh_bangle, & units = data%Lev2c%units%pblh_bangle, & range = data%Lev2c%range%pblh_bangle, & rec = irec) ENDIF IF (ncdf_isvar('pblh_bangle2')) THEN !For backward compatibility CALL ncdf_getvar('pblh_bangle2', data%Lev2c%pblh_bangle2, & units = data%Lev2c%units%pblh_bangle, & range = data%Lev2c%range%pblh_bangle, & rec = irec) ENDIF IF (ncdf_isvar('pbla_bangle')) THEN !For backward compatibility CALL ncdf_getvar('pbla_bangle', data%Lev2c%pbla_bangle, & units = data%Lev2c%units%pbla_bangle, & range = data%Lev2c%range%pbla_bangle, & rec = irec) ENDIF IF (ncdf_isvar('pbla_bangle2')) THEN !For backward compatibility CALL ncdf_getvar('pbla_bangle2', data%Lev2c%pbla_bangle2, & units = data%Lev2c%units%pbla_bangle, & range = data%Lev2c%range%pbla_bangle, & rec = irec) ENDIF IF (ncdf_isvar('pblh_bangle_flag')) THEN !For backward compatibility CALL ncdf_getvar('pblh_bangle_flag', data%Lev2c%pblh_bangle_flag, & units = data%Lev2c%units%pblh_bangle_flag, & range = data%Lev2c%range%pblh_bangle_flag, & rec = irec) ENDIF IF (ncdf_isvar('pblh_refrac')) THEN !For backward compatibility CALL ncdf_getvar('pblh_refrac', data%Lev2c%pblh_refrac, & units = data%Lev2c%units%pblh_refrac, & range = data%Lev2c%range%pblh_refrac, & rec = irec) ENDIF IF (ncdf_isvar('pblh_refrac2')) THEN !For backward compatibility CALL ncdf_getvar('pblh_refrac2', data%Lev2c%pblh_refrac2, & units = data%Lev2c%units%pblh_refrac, & range = data%Lev2c%range%pblh_refrac, & rec = irec) ENDIF IF (ncdf_isvar('pbln_refrac')) THEN !For backward compatibility CALL ncdf_getvar('pbln_refrac', data%Lev2c%pbln_refrac, & units = data%Lev2c%units%pbln_refrac, & range = data%Lev2c%range%pbln_refrac, & rec = irec) ENDIF IF (ncdf_isvar('pbln_refrac2')) THEN !For backward compatibility CALL ncdf_getvar('pbln_refrac2', data%Lev2c%pbln_refrac2, & units = data%Lev2c%units%pbln_refrac, & range = data%Lev2c%range%pbln_refrac, & rec = irec) ENDIF IF (ncdf_isvar('pblh_refrac_flag')) THEN !For backward compatibility CALL ncdf_getvar('pblh_refrac_flag', data%Lev2c%pblh_refrac_flag, & units = data%Lev2c%units%pblh_refrac_flag, & range = data%Lev2c%range%pblh_refrac_flag, & rec = irec) ENDIF IF (ncdf_isvar('pblh_tdry')) THEN !For backward compatibility CALL ncdf_getvar('pblh_tdry', data%Lev2c%pblh_tdry, & units = data%Lev2c%units%pblh_tdry, & range = data%Lev2c%range%pblh_tdry, & rec = irec) ENDIF IF (ncdf_isvar('pblh_tdry2')) THEN !For backward compatibility CALL ncdf_getvar('pblh_tdry2', data%Lev2c%pblh_tdry2, & units = data%Lev2c%units%pblh_tdry, & range = data%Lev2c%range%pblh_tdry, & rec = irec) ENDIF IF (ncdf_isvar('pblt_tdry')) THEN !For backward compatibility CALL ncdf_getvar('pblt_tdry', data%Lev2c%pblt_tdry, & units = data%Lev2c%units%pblt_tdry, & range = data%Lev2c%range%pblt_tdry, & rec = irec) ENDIF IF (ncdf_isvar('pblt_tdry2')) THEN !For backward compatibility CALL ncdf_getvar('pblt_tdry2', data%Lev2c%pblt_tdry2, & units = data%Lev2c%units%pblt_tdry, & range = data%Lev2c%range%pblt_tdry, & rec = irec) ENDIF IF (ncdf_isvar('pblh_tdry_flag')) THEN !For backward compatibility CALL ncdf_getvar('pblh_tdry_flag', data%Lev2c%pblh_tdry_flag, & units = data%Lev2c%units%pblh_tdry_flag, & range = data%Lev2c%range%pblh_tdry_flag, & rec = irec) ENDIF IF (ncdf_isvar('pblh_temp')) THEN !For backward compatibility CALL ncdf_getvar('pblh_temp', data%Lev2c%pblh_temp, & units = data%Lev2c%units%pblh_temp, & range = data%Lev2c%range%pblh_temp, & rec = irec) ENDIF IF (ncdf_isvar('pblh_temp2')) THEN !For backward compatibility CALL ncdf_getvar('pblh_temp2', data%Lev2c%pblh_temp2, & units = data%Lev2c%units%pblh_temp, & range = data%Lev2c%range%pblh_temp, & rec = irec) ENDIF IF (ncdf_isvar('pblt_temp')) THEN !For backward compatibility CALL ncdf_getvar('pblt_temp', data%Lev2c%pblt_temp, & units = data%Lev2c%units%pblt_temp, & range = data%Lev2c%range%pblt_temp, & rec = irec) ENDIF IF (ncdf_isvar('pblt_temp2')) THEN !For backward compatibility CALL ncdf_getvar('pblt_temp2', data%Lev2c%pblt_temp2, & units = data%Lev2c%units%pblt_temp, & range = data%Lev2c%range%pblt_temp, & rec = irec) ENDIF IF (ncdf_isvar('pblh_temp_flag')) THEN !For backward compatibility CALL ncdf_getvar('pblh_temp_flag', data%Lev2c%pblh_temp_flag, & units = data%Lev2c%units%pblh_temp_flag, & range = data%Lev2c%range%pblh_temp_flag, & rec = irec) ENDIF IF (ncdf_isvar('pblh_shum')) THEN !For backward compatibility CALL ncdf_getvar('pblh_shum', data%Lev2c%pblh_shum, & units = data%Lev2c%units%pblh_shum, & range = data%Lev2c%range%pblh_shum, & rec = irec) ENDIF IF (ncdf_isvar('pblh_shum2')) THEN !For backward compatibility CALL ncdf_getvar('pblh_shum2', data%Lev2c%pblh_shum2, & units = data%Lev2c%units%pblh_shum, & range = data%Lev2c%range%pblh_shum, & rec = irec) ENDIF IF (ncdf_isvar('pblq_shum')) THEN !For backward compatibility CALL ncdf_getvar('pblq_shum', data%Lev2c%pblq_shum, & units = data%Lev2c%units%pblq_shum, & range = data%Lev2c%range%pblq_shum, & rec = irec) ENDIF IF (ncdf_isvar('pblq_shum2')) THEN !For backward compatibility CALL ncdf_getvar('pblq_shum2', data%Lev2c%pblq_shum2, & units = data%Lev2c%units%pblq_shum, & range = data%Lev2c%range%pblq_shum, & rec = irec) ENDIF IF (ncdf_isvar('pblh_shum_flag')) THEN !For backward compatibility CALL ncdf_getvar('pblh_shum_flag', data%Lev2c%pblh_shum_flag, & units = data%Lev2c%units%pblh_shum_flag, & range = data%Lev2c%range%pblh_shum_flag, & rec = irec) ENDIF IF (ncdf_isvar('pblh_rhum')) THEN !For backward compatibility CALL ncdf_getvar('pblh_rhum', data%Lev2c%pblh_rhum, & units = data%Lev2c%units%pblh_rhum, & range = data%Lev2c%range%pblh_rhum, & rec = irec) ENDIF IF (ncdf_isvar('pblh_rhum2')) THEN !For backward compatibility CALL ncdf_getvar('pblh_rhum2', data%Lev2c%pblh_rhum2, & units = data%Lev2c%units%pblh_rhum, & range = data%Lev2c%range%pblh_rhum, & rec = irec) ENDIF IF (ncdf_isvar('pblr_rhum')) THEN !For backward compatibility CALL ncdf_getvar('pblr_rhum', data%Lev2c%pblr_rhum, & units = data%Lev2c%units%pblr_rhum, & range = data%Lev2c%range%pblr_rhum, & rec = irec) ENDIF IF (ncdf_isvar('pblr_rhum2')) THEN !For backward compatibility CALL ncdf_getvar('pblr_rhum2', data%Lev2c%pblr_rhum2, & units = data%Lev2c%units%pblr_rhum, & range = data%Lev2c%range%pblr_rhum, & rec = irec) ENDIF IF (ncdf_isvar('pblh_rhum_flag')) THEN !For backward compatibility CALL ncdf_getvar('pblh_rhum_flag', data%Lev2c%pblh_rhum_flag, & units = data%Lev2c%units%pblh_rhum_flag, & range = data%Lev2c%range%pblh_rhum_flag, & rec = irec) ENDIF data%Lev2c%Missing = .FALSE. ENDIF ! 1.15 Level2d variables (if any) ! ------------------------------- IF (ncdf_isvar('level_coeff_a')) THEN CALL ncdf_getsize('level_coeff_a', n, dim = 1) CALL ropp_io_init(data%Lev2d, n) ELSE data%Lev2d%Npoints = 0 ENDIF IF (data%Lev2d%Npoints > 0) THEN CALL ncdf_getvar('level_type', data%Lev2d%level_type, & rec = irec) CALL ncdf_getvar('level_coeff_a', data%Lev2d%level_coeff_a, & units = data%Lev2d%units%level_coeff_a, & range = data%Lev2d%range%level_coeff_a, & rec = irec) CALL ncdf_getvar('level_coeff_b', data%Lev2d%level_coeff_b, & units = data%Lev2d%units%level_coeff_b, & range = data%Lev2d%range%level_coeff_b, & rec = irec) data%Lev2d%Missing = .FALSE. ENDIF ! 1.16 Level2e variables (if any) ! ------------------------------- n1 = 0 IF (ncdf_isvar('n_e')) THEN CALL ncdf_getsize('n_e', n1, dim = 1) ELSE data%Lev2e%Npoints = 0 ENDIF n2 = 0 IF (ncdf_isvar('ne_peak')) THEN CALL ncdf_getsize('ne_peak', n2, dim = 1) ELSE data%Lev2e%Nlayers = 0 ENDIF IF (ncdf_isvar('n_e') .OR. ncdf_isvar('ne_peak')) THEN CALL ropp_io_init(data%Lev2e, (/n1, n2/)) ENDIF IF (data%Lev2e%Npoints > 0) THEN CALL ncdf_getvar('n_e', data%Lev2e%n_e, & units = data%Lev2e%units%n_e, & range = data%Lev2e%range%n_e, & rec = irec) CALL ncdf_getvar('n_e_sigma', data%Lev2e%n_e_sigma, & units = data%Lev2e%units%n_e_sigma, & range = data%Lev2e%range%n_e_sigma, & rec = irec) CALL ncdf_getvar('r_iono', data%Lev2e%r_iono, & units = data%Lev2e%units%r_iono, & range = data%Lev2e%range%r_iono, & rec = irec) data%Lev2e%Missing = .FALSE. ENDIF IF (data%Lev2e%Nlayers > 0) THEN CALL ncdf_getvar('ne_peak', data%Lev2e%ne_peak, & units = data%Lev2e%units%ne_peak, & range = data%Lev2e%range%ne_peak, & rec = irec) CALL ncdf_getvar('ne_peak_sigma', data%Lev2e%ne_peak_sigma, & units = data%Lev2e%units%ne_peak_sigma, & range = data%Lev2e%range%ne_peak_sigma, & rec = irec) CALL ncdf_getvar('r_peak', data%Lev2e%r_peak, & units = data%Lev2e%units%r_peak, & range = data%Lev2e%range%r_peak, & rec = irec) CALL ncdf_getvar('r_peak_sigma', data%Lev2e%r_peak_sigma, & units = data%Lev2e%units%r_peak_sigma, & range = data%Lev2e%range%r_peak_sigma, & rec = irec) CALL ncdf_getvar('h_zero', data%Lev2e%h_zero, & units = data%Lev2e%units%h_zero, & range = data%Lev2e%range%h_zero, & rec = irec) CALL ncdf_getvar('h_zero_sigma', data%Lev2e%h_zero_sigma, & units = data%Lev2e%units%h_zero_sigma, & range = data%Lev2e%range%h_zero_sigma, & rec = irec) CALL ncdf_getvar('h_grad', data%Lev2e%h_grad, & units = data%Lev2e%units%h_grad, & range = data%Lev2e%range%h_grad, & rec = irec) CALL ncdf_getvar('h_grad_sigma', data%Lev2e%h_grad_sigma, & units = data%Lev2e%units%h_grad_sigma, & range = data%Lev2e%range%h_grad_sigma, & rec = irec) data%Lev2e%Missing = .FALSE. ENDIF ! 1.17 Additional variables (if any) ! ---------------------------------- CALL ropp_io_init(data%vlist) DO varid=1,ncdf_nvars IF(.NOT. ncdf_read(varid))THEN status = nf90_inquire_variable(ncdf_ncid, varid, xtype=TYPE, ndims=ndim) IF (TYPE .NE. NF90_CHAR) THEN ! only read scalar variables IF(ndim == 1)THEN CALL ropp_io_read_ncdf_get_vlistD0d(varid, data%vlist%VlistD0d, irec) ENDIF IF(ndim == 2)THEN CALL ropp_io_read_ncdf_get_vlistD1d(varid, data%vlist%VlistD1d, irec) ENDIF IF(ndim == 3)THEN CALL ropp_io_read_ncdf_get_vlistD2d(varid, data%vlist%VlistD2d, irec) ENDIF ENDIF ENDIF ncdf_read(varid) = .FALSE. ! reset 'read variable' flag ENDDO ! 1.17 Clean up ! ------------- CALL message_set_routine(routine) END SUBROUTINE ropp_io_read_ncdf_get_rodata !------------------------------------------------------------------------------- ! 2. Core RO data (two-dimensional meteorological data) !------------------------------------------------------------------------------- SUBROUTINE ropp_io_read_ncdf_get_rodata_2d(DATA, rec) ! 2.1 Declarations ! ---------------- USE ropp_utils USE ncdf USE ropp_io, not_this => ropp_io_read_ncdf_get_rodata_2d USE ropp_io_types, ONLY: ROprof2d, & ThisFmtVer IMPLICIT NONE TYPE(ROprof2d), INTENT(inout) :: DATA INTEGER, OPTIONAL :: rec INTEGER :: n, ierr INTEGER :: irec INTEGER, DIMENSION(2) :: n2d INTEGER, DIMENSION(8) :: DT8 ! Date/time array CHARACTER(LEN=15) :: cval REAL(dp) :: time, start_time, dtocc_time, now_time !dp defined in DateTimeTypes CHARACTER(len = 23) :: proc_date CHARACTER(len = 256) :: routine REAL :: fmtver INTEGER :: status, varid, ndim, TYPE ! 2.2 Error handling ! ------------------ CALL message_get_routine(routine) CALL message_set_routine('ropp_io_read_ncdf_get') ! 2.3 Default parameters ! ---------------------- IF (PRESENT(rec)) THEN irec = rec ELSE irec = 1 ENDIF ! 2.4 (Global) Attributes ! ------------------------ data%FmtVersion = ' ' ; CALL ncdf_getatt('format_version', data%FmtVersion) READ ( data%FmtVersion(11:), fmt=*, iostat=ierr ) fmtver IF ( ierr /= 0 ) data%FmtVersion = ThisFmtVer data%processing_centre = ' ' ; CALL ncdf_getatt('processing_centre', data%processing_centre) IF (ncdf_isatt('processing_software')) THEN ! added at V8.0 data%processing_software = ' ' ; CALL ncdf_getatt('processing_software', data%processing_software) ENDIF proc_date = ' ' ; CALL ncdf_getatt('processing_date', proc_date) data%pod_method = ' ' ; CALL ncdf_getatt('pod_method', data%pod_method) data%phase_method = ' ' ; CALL ncdf_getatt('phase_method', data%phase_method) data%bangle_method = ' ' ; CALL ncdf_getatt('bangle_method', data%bangle_method) data%refrac_method = ' ' ; CALL ncdf_getatt('refrac_method', data%refrac_method) data%meteo_method = ' ' ; CALL ncdf_getatt('meteo_method', data%meteo_method) IF(ncdf_isatt('thin_method'))THEN ! added at V1.1 data%thin_method = ' ' ; CALL ncdf_getatt('thin_method', data%thin_method) ENDIF data%software_version = ' ' ; CALL ncdf_getatt('software_version', data%software_version) data%bad = ' ' ; CALL ncdf_getatt('bad', data%bad) IF (proc_date( 1: 4) /= ' ') READ(proc_date( 1: 4), *) data%DTpro%Year IF (proc_date( 6: 7) /= ' ') READ(proc_date( 6: 7), *) data%DTpro%Month IF (proc_date( 9:10) /= ' ') READ(proc_date( 9:10), *) data%DTpro%Day IF (proc_date(12:13) /= ' ') READ(proc_date(12:13), *) data%DTpro%Hour IF (proc_date(15:16) /= ' ') READ(proc_date(15:16), *) data%DTpro%Minute IF (proc_date(18:19) /= ' ') READ(proc_date(18:19), *) data%DTpro%Second IF (proc_date(21:23) /= ' ') READ(proc_date(21:23), *) data%DTpro%Msec ! 2.5 Header variables ! -------------------- CALL ncdf_getvar('occ_id', data%occ_id, rec = irec) CALL ncdf_getvar('gns_id', data%gns_id, rec = irec) CALL ncdf_getvar('leo_id', data%leo_id, rec = irec) CALL ncdf_getvar('stn_id', data%stn_id, rec = irec) CALL ncdf_getvar('occfreq1', data%occfreq1, rec = irec) CALL ncdf_getvar('occfreq2', data%occfreq2, rec = irec) IF ( data%leo_id(1:2) == 'CH' ) THEN data%leo_id = 'CHMP' ELSE IF ( data%leo_id(1:4) == 'GRC1' ) THEN data%leo_id = 'GRAA' ELSE IF ( data%leo_id(1:4) == 'GRC2' ) THEN data%leo_id = 'GRAB' ELSE IF ( data%leo_id(1:2) == 'CO' ) THEN data%leo_id(1:2) = 'C0' ELSE IF ( data%leo_id(1:3) == 'MTP' ) THEN data%leo_id(1:3) = 'MET' ELSE IF ( data%leo_id(1:4) == 'G001' ) THEN ! shouldn't ever be needed data%leo_id = 'FY3C' ELSE IF ( data%leo_id(1:4) == 'G002' ) THEN ! shouldn't ever be needed data%leo_id = 'FY3D' END IF ! 2.6 Date and time ! ----------------- CALL ncdf_getvar('start_time', start_time, rec=irec) CALL ncdf_getvar('year', data%DTocc%Year, & units = data%DTocc%units%Year, & range = data%DTocc%range%Year, & rec = irec) CALL ncdf_getvar('month', data%DTocc%Month, & units = data%DTocc%units%Month, & range = data%DTocc%range%Month, & rec = irec) CALL ncdf_getvar('day', data%DTocc%Day, & units = data%DTocc%units%Day, & range = data%DTocc%range%Day, & rec = irec) CALL ncdf_getvar('hour', data%DTocc%Hour, & units = data%DTocc%units%Hour, & range = data%DTocc%range%Hour, & rec = irec) CALL ncdf_getvar('minute', data%DTocc%Minute, & units = data%DTocc%units%Minute, & range = data%DTocc%range%Minute, & rec = irec) CALL ncdf_getvar('second', data%DTocc%Second, & units = data%DTocc%units%Second, & range = data%DTocc%range%Second, & rec = irec) CALL ncdf_getvar('msec', data%DTocc%Msec, & units = data%DTocc%units%Msec, & range = data%DTocc%range%Msec, & rec = irec) ! 2.6.1 Check consistency: start_time and DTocc (if both are valid) ! should refer to the same epoch within 1ms. Issue a warning if not. CALL Date_and_Time_UTC(Values=DT8) CALL TimeSince(DT8, now_time, 1, Base="JS2000") IF (isroppinrange(data%DTocc)) THEN IF ( isinrange(start_time, (/ 0.001_dp, now_time /)) ) THEN DT8 = (/ data%DTocc%Year, data%DTocc%Month, data%DTocc%Day, 0, & data%DTocc%Hour, data%DTocc%Minute, data%DTocc%Second, & data%DTocc%Msec /) CALL TimeSince(DT8, dtocc_time, 1, Base="JS2000") dtocc_time = dtocc_time + (LeapSeconds(DT8) - 22)*1.D0 ! 22 leapseconds had elapsed by 0Z 01/01/2000 time = ABS(start_time - dtocc_time) IF ( time > 0.0005_dp ) THEN WRITE ( cval, FMT="(F15.3)" ) time IF ( time < 30.0_dp ) THEN ! Probably a result of not accounting for leap-seconds, so reduce severity of message to DIAG CALL message(msg_diag, "'start_time' and yr/mo/dy/hr/mn/sc/ms " // & "timestamps differ by " // & TRIM(ADJUSTL(cval))//" seconds (probably a leap-second issue) - using yr/../ms timestamp") ELSE CALL message(msg_warn, "'start_time' and yr/mo/dy/hr/mn/sc/ms " // & "timestamps differ by " // & TRIM(ADJUSTL(cval))//" seconds - using yr/../ms timestamp") END IF END IF END IF ! If any DTocc element is invalid, substitute converted start_time ! (if that is valid) and issue a warning. ELSE IF ( isinrange(start_time, (/ 0.001_dp, now_time /)) ) THEN CALL message(msg_warn, "One or more of yr/mo/dy/hr/mn/sc/ms times " // & "are invalid - using 'start_time' instead") CALL TimeSince(DT8, start_time, -1, Base="JS2000") start_time = start_time - (LeapSeconds(DT8) - 22)*1.D0 ! 22 leapseconds had elapsed by 0Z 01/01/2000 CALL TimeSince(DT8, start_time, -1, Base="JS2000") data%DTocc%Year = DT8(1) data%DTocc%Month = DT8(2) data%DTocc%Day = DT8(3) data%DTocc%Hour = DT8(5) data%DTocc%Minute = DT8(6) data%DTocc%Second = DT8(7) data%DTocc%Msec = DT8(8) END IF ! NB if neither DTocc nor start_time are valid, do nothing here; missing ! coordinates should be dealt with as part of generic Q/C. END IF ! 2.7 Overall quality ! ------------------- CALL ncdf_getvar('pcd', data%pcd, & units = data%units%pcd, & range = data%range%pcd, & rec = irec) CALL ncdf_getvar('overall_qual', data%overall_qual, & units = data%units%overall_qual, & range = data%range%overall_qual, & rec = irec) ! 2.8 Georeferencing ! ------------------ CALL ncdf_getvar('time', time, rec=irec) CALL ncdf_getvar('lat', data%georef%lat, & units = data%georef%units%lat, & range = data%georef%range%lat, & rec = irec) CALL ncdf_getvar('lon', data%georef%lon, & units = data%georef%units%lon, & range = data%georef%range%lon, & rec = irec) CALL ncdf_getvar('time_offset', data%georef%time_offset, & units = data%georef%units%time_offset, & range = data%georef%range%time_offset, & rec = irec) CALL ncdf_getvar('undulation', data%georef%Undulation, & units = data%georef%units%Undulation, & range = data%georef%range%undulation, & rec = irec) ! Yong Chen on 05/30/2024 set Undulation to zero if not in the range !write(*, *) 'Undulation range is: ', data%georef%range%undulation IF ( data%georef%Undulation < data%georef%range%undulation(1) .or. & data%georef%Undulation > data%georef%range%undulation(2) ) & data%georef%Undulation = 0.0 CALL ncdf_getvar('roc', data%georef%roc, & units = data%georef%units%roc, & range = data%georef%range%roc, & rec = irec) CALL ncdf_getvar('r_coc', data%georef%r_coc, & units = data%georef%units%r_coc, & range = data%georef%range%r_coc, & rec = irec) CALL ncdf_getvar('r_alt', data%georef%r_alt, & units = data%georef%units%r_alt, & range = data%georef%range%r_alt, & rec = irec) CALL ncdf_getvar('viewang', data%georef%viewang, & units = data%georef%units%viewang, & range = data%georef%range%viewang, & rec = irec) CALL ncdf_getvar('azimuth', data%georef%azimuth, & units = data%georef%units%azimuth, & range = data%georef%range%azimuth, & rec = irec) IF ( ncdf_isvar('r_leo_ref') ) & CALL ncdf_getvar('r_leo_ref', data%georef%leo_pod%pos, & units = data%georef%units%r_leo, & range = data%georef%range%r_leo, & rec = irec) IF ( ncdf_isvar('v_leo_ref') ) & CALL ncdf_getvar('v_leo_ref', data%georef%leo_pod%vel, & units = data%georef%units%v_leo, & range = data%georef%range%v_leo, & rec = irec) IF ( ncdf_isvar('r_gns_ref') ) & CALL ncdf_getvar('r_gns_ref', data%georef%gns_pod%pos, & units = data%georef%units%r_gns, & range = data%georef%range%r_gns, & rec = irec) IF ( ncdf_isvar('v_gns_ref') ) & CALL ncdf_getvar('v_gns_ref', data%georef%gns_pod%vel, & units = data%georef%units%v_gns, & range = data%georef%range%v_gns, & rec = irec) ! 2.8.1 Other attributes CALL ncdf_getatt('reference_frame', data%georef%reference_frame%r_coc, varname= 'r_coc') IF ( ncdf_isvar('r_leo_ref') ) & CALL ncdf_getatt('reference_frame', data%georef%leo_pod%pos_frame, varname= 'r_leo_ref') IF ( ncdf_isvar('v_leo_ref') ) & CALL ncdf_getatt('reference_frame', data%georef%leo_pod%vel_frame, varname= 'v_leo_ref') IF ( ncdf_isvar('r_gns_ref') ) & CALL ncdf_getatt('reference_frame', data%georef%gns_pod%pos_frame, varname= 'r_gns_ref') IF ( ncdf_isvar('v_gns_ref') ) & CALL ncdf_getatt('reference_frame', data%georef%gns_pod%vel_frame, varname= 'v_gns_ref') ! 2.9 Background characterisation (if any) ! ---------------------------------------- IF (ncdf_isvar('bg_source')) THEN data%BG%Source = 'TBD' ELSE data%BG%Source = 'NONE' ENDIF IF (data%BG%Source /= 'NONE') THEN CALL ncdf_getvar('bg_source', data%BG%Source, rec = irec) CALL ncdf_getvar('bg_year', data%BG%Year, & units = data%BG%units%Year, & range = data%BG%range%Year, & rec = irec) CALL ncdf_getvar('bg_month', data%BG%Month, & units = data%BG%units%Month, & range = data%BG%range%Month, & rec = irec) CALL ncdf_getvar('bg_day', data%BG%Day, & units = data%BG%units%Day, & range = data%BG%range%Day, & rec = irec) CALL ncdf_getvar('bg_hour', data%BG%Hour, & units = data%BG%units%Hour, & range = data%BG%range%Hour, & rec = irec) CALL ncdf_getvar('bg_minute', data%BG%Minute, & units = data%BG%units%Minute, & range = data%BG%range%Minute, & rec = irec) CALL ncdf_getvar('bg_fcperiod', data%BG%fcPeriod, & units = data%BG%units%fcPeriod, & range = data%BG%range%fcPeriod, & rec = irec) ENDIF ! 2.9.1 QA characterisation ! ---------------------------------------- CALL ncdf_getvar('snr1avg', data%QA%snr1avg, & units = data%QA%units%snr1avg, & range = data%QA%range%snr1avg, & rec = irec) CALL ncdf_getvar('snr2avg', data%QA%snr2avg, & units = data%QA%units%snr2avg, & range = data%QA%range%snr2avg, & rec = irec) CALL ncdf_getvar('irs', data%QA%irs, & units = data%QA%units%irs, & range = data%QA%range%irs, & rec = irec) CALL ncdf_getvar('smean', data%QA%smean, & units = data%QA%units%smean, & range = data%QA%range%smean, & rec = irec) CALL ncdf_getvar('stdv', data%QA%stdv, & units = data%QA%units%stdv, & range = data%QA%range%stdv, & rec = irec) CALL ncdf_getvar('reldevmax', data%QA%reldevmax, & units = data%QA%units%reldevmax, & range = data%QA%range%reldevmax, & rec = irec) CALL ncdf_getvar('reldevstd', data%QA%reldevstd, & units = data%QA%units%reldevstd, & range = data%QA%range%reldevstd, & rec = irec) CALL ncdf_getvar('reldevmax1', data%QA%reldevmax1, & units = data%QA%units%reldevmax1, & range = data%QA%range%reldevmax1, & rec = irec) CALL ncdf_getvar('reldevstd1', data%QA%reldevstd1, & units = data%QA%units%reldevstd1, & range = data%QA%range%reldevstd1, & rec = irec) CALL ncdf_getvar('difmaxref', data%QA%difmaxref, & units = data%QA%units%difmaxref, & range = data%QA%range%difmaxref, & rec = irec) CALL ncdf_getvar('maxdifphase', data%QA%maxdifphase, & units = data%QA%units%maxdifphase, & range = data%QA%range%maxdifphase, & rec = irec) CALL ncdf_getvar('finalQA', data%QA%finalQA, & units = data%QA%units%finalQA, & range = data%QA%range%finalQA, & rec = irec) ! 2.10 Level1a variables (if any) ! ------------------------------ IF (ncdf_isvar('dtime')) THEN CALL ncdf_getsize('dtime', n, dim = 1) CALL ropp_io_init(data%Lev1a, n) ELSE data%Lev1a%Npoints = 0 ENDIF IF (data%Lev1a%Npoints > 0) THEN CALL ncdf_getvar('dtime', data%Lev1a%dtime, & units = data%Lev1a%units%dtime, & range = data%Lev1a%range%dtime, & rec = irec) CALL ncdf_getvar('snr_L1ca', data%Lev1a%snr_L1ca, & units = data%Lev1a%units%snr, & range = data%Lev1a%range%snr, & rec = irec) CALL ncdf_getvar('snr_L1p', data%Lev1a%snr_L1p, & units = data%Lev1a%units%snr, & range = data%Lev1a%range%snr, & rec = irec) CALL ncdf_getvar('snr_L2p', data%Lev1a%snr_L2p, & units = data%Lev1a%units%snr, & range = data%Lev1a%range%snr, & rec = irec) CALL ncdf_getvar('phase_L1', data%Lev1a%phase_L1, & units = data%Lev1a%units%phase, & range = data%Lev1a%range%phase, & rec = irec) CALL ncdf_getvar('phase_L2', data%Lev1a%phase_L2, & units = data%Lev1a%units%phase, & range = data%Lev1a%range%phase, & rec = irec) CALL ncdf_getvar('r_gns', data%Lev1a%r_gns, & units = data%Lev1a%units%r_gns, & range = data%Lev1a%range%r_gns, & rec = irec) CALL ncdf_getvar('v_gns', data%Lev1a%v_gns, & units = data%Lev1a%units%v_gns, & range = data%Lev1a%range%v_gns, & rec = irec) CALL ncdf_getvar('r_leo', data%Lev1a%r_leo, & units = data%Lev1a%units%r_leo, & range = data%Lev1a%range%r_leo, & rec = irec) CALL ncdf_getvar('v_leo', data%Lev1a%v_leo, & units = data%Lev1a%units%v_leo, & range = data%Lev1a%range%v_leo, & rec = irec) CALL ncdf_getvar('phase_qual', data%Lev1a%phase_qual, & units = data%Lev1a%units%phase_qual, & range = data%Lev1a%range%phase_qual, & rec = irec) ! 2.10.1 Other attributes CALL ncdf_getatt('reference_frame', data%Lev1a%reference_frame%r_gns, varname = 'r_gns') CALL ncdf_getatt('reference_frame', data%Lev1a%reference_frame%v_gns, varname = 'v_gns') CALL ncdf_getatt('reference_frame', data%Lev1a%reference_frame%r_leo, varname = 'r_leo') CALL ncdf_getatt('reference_frame', data%Lev1a%reference_frame%v_leo, varname = 'v_leo') ENDIF ! 2.11 Level1b variables (if any) ! ------------------------------- IF (ncdf_isvar('lat_tp')) THEN CALL ncdf_getsize('lat_tp', n, dim = 1) CALL ropp_io_init(data%Lev1b, n) ELSE data%Lev1b%Npoints = 0 ENDIF IF (data%Lev1b%Npoints > 0) THEN CALL ncdf_getvar('lat_tp', data%Lev1b%lat_tp, & units = data%Lev1b%units%lat_tp, & range = data%Lev1b%range%lat_tp, & rec = irec) CALL ncdf_getvar('lon_tp', data%Lev1b%lon_tp, & units = data%Lev1b%units%lon_tp, & range = data%Lev1b%range%lon_tp, & rec = irec) CALL ncdf_getvar('azimuth_tp', data%Lev1b%azimuth_tp, & units = data%Lev1b%units%azimuth_tp, & range = data%Lev1b%range%azimuth_tp, & rec = irec) CALL ncdf_getvar('impact_L1', data%Lev1b%impact_L1, & units = data%Lev1b%units%impact, & range = data%Lev1b%range%impact, & rec = irec) CALL ncdf_getvar('impact_L2', data%Lev1b%impact_L2, & units = data%Lev1b%units%impact, & range = data%Lev1b%range%impact, & rec = irec) CALL ncdf_getvar('impact', data%Lev1b%impact, & units = data%Lev1b%units%impact, & range = data%Lev1b%range%impact, & rec = irec) IF (ncdf_isvar('impact_opt')) & ! added at v1.1 CALL ncdf_getvar('impact_opt', data%Lev1b%impact_opt, & units = data%Lev1b%units%impact, & range = data%Lev1b%range%impact, & rec = irec) CALL ncdf_getvar('bangle_L1', data%Lev1b%bangle_L1, & units = data%Lev1b%units%bangle, & range = data%Lev1b%range%bangle, & rec = irec) CALL ncdf_getvar('bangle_L2', data%Lev1b%bangle_L2, & units = data%Lev1b%units%bangle, & range = data%Lev1b%range%bangle, & rec = irec) CALL ncdf_getvar('bangle', data%Lev1b%bangle, & units = data%Lev1b%units%bangle, & range = data%Lev1b%range%bangle, & rec = irec) IF (ncdf_isvar('bangle_opt')) & ! added at v1.1 CALL ncdf_getvar('bangle_opt', data%Lev1b%bangle_opt, & units = data%Lev1b%units%bangle, & range = data%Lev1b%range%bangle, & rec = irec) CALL ncdf_getvar('bangle_L1_sigma', data%Lev1b%bangle_L1_sigma, & units = data%Lev1b%units%bangle_sigma, & range = data%Lev1b%range%bangle_sigma, & rec = irec) CALL ncdf_getvar('bangle_L2_sigma', data%Lev1b%bangle_L2_sigma, & units = data%Lev1b%units%bangle_sigma, & range = data%Lev1b%range%bangle_sigma, & rec = irec) CALL ncdf_getvar('bangle_sigma', data%Lev1b%bangle_sigma, & units = data%Lev1b%units%bangle_sigma, & range = data%Lev1b%range%bangle_sigma, & rec = irec) IF (ncdf_isvar('bangle_opt_sigma')) & ! added at v1.1 CALL ncdf_getvar('bangle_opt_sigma', data%Lev1b%bangle_opt_sigma, & units = data%Lev1b%units%bangle_sigma, & range = data%Lev1b%range%bangle_sigma, & rec = irec) CALL ncdf_getvar('bangle_L1_qual', data%Lev1b%bangle_L1_qual, & units = data%Lev1b%units%bangle_qual, & range = data%Lev1b%range%bangle_qual, & rec = irec) CALL ncdf_getvar('bangle_L2_qual', data%Lev1b%bangle_L2_qual, & units = data%Lev1b%units%bangle_qual, & range = data%Lev1b%range%bangle_qual, & rec = irec) CALL ncdf_getvar('bangle_qual', data%Lev1b%bangle_qual, & units = data%Lev1b%units%bangle_qual, & range = data%Lev1b%range%bangle_qual, & rec = irec) IF (ncdf_isvar('bangle_opt_qual')) & ! added at v1.1 CALL ncdf_getvar('bangle_opt_qual', data%Lev1b%bangle_opt_qual, & units = data%Lev1b%units%bangle_qual, & range = data%Lev1b%range%bangle_qual, & rec = irec) ENDIF ! 2.12 Level2a variables (if any) ! ------------------------------- IF (ncdf_isvar('alt_refrac')) THEN n2d(1) = 0 n2d(2) = 0 CALL ncdf_getsize('alt_refrac', n2d) CALL ropp_io_init(data%Lev2a, n2d) ELSE data%Lev2a%Npoints = 0 data%Lev2a%Nhoriz = 0 ENDIF IF (data%Lev2a%Npoints > 0) THEN CALL ncdf_getvar('alt_refrac', data%Lev2a%alt_refrac, & units = data%Lev2a%units%alt_refrac, & range = data%Lev2a%range%alt_refrac, & rec = irec) CALL ncdf_getvar('geop_refrac', data%Lev2a%geop_refrac, & units = data%Lev2a%units%geop_refrac, & range = data%Lev2a%range%geop_refrac, & rec = irec) CALL ncdf_getvar('refrac', data%Lev2a%refrac, & units = data%Lev2a%units%refrac, & range = data%Lev2a%range%refrac, & rec = irec) CALL ncdf_getvar('refrac_sigma', data%Lev2a%refrac_sigma, & units = data%Lev2a%units%refrac_sigma, & range = data%Lev2a%range%refrac_sigma, & rec = irec) CALL ncdf_getvar('refrac_qual', data%Lev2a%refrac_qual, & units = data%Lev2a%units%refrac_qual, & range = data%Lev2a%range%refrac_qual, & rec = irec) IF (ncdf_isvar('dry_temp')) THEN !For backward compatibility CALL ncdf_getvar('dry_temp', data%Lev2a%dry_temp, & units = data%Lev2a%units%dry_temp, & range = data%Lev2a%range%dry_temp, & rec = irec) CALL ncdf_getvar('dry_temp_sigma', data%Lev2a%dry_temp_sigma, & units = data%Lev2a%units%dry_temp_sigma, & range = data%Lev2a%range%dry_temp_sigma, & rec = irec) CALL ncdf_getvar('dry_temp_qual', data%Lev2a%dry_temp_qual, & units = data%Lev2a%units%dry_temp_qual, & range = data%Lev2a%range%dry_temp_qual, & rec = irec) ENDIF ENDIF ! 2.13 Level2b variables (if any) ! ------------------------------- IF (ncdf_isvar('geop')) THEN n2d(1) = 0 n2d(2) = 0 CALL ncdf_getsize('geop', n2d) CALL ropp_io_init(data%Lev2b, n2d) ELSE data%Lev2b%Npoints = 0 data%Lev2b%Nhoriz = 0 ENDIF IF (data%Lev2b%Npoints > 0) THEN CALL ncdf_getvar('geop', data%Lev2b%geop, & units = data%Lev2b%units%geop, & range = data%Lev2b%range%geop, & rec = irec) CALL ncdf_getvar('geop_sigma', data%Lev2b%geop_sigma, & units = data%Lev2b%units%geop_sigma, & range = data%Lev2b%range%geop_sigma, & rec = irec) CALL ncdf_getvar('press', data%Lev2b%press, & units = data%Lev2b%units%press, & range = data%Lev2b%range%press, & rec = irec) CALL ncdf_getvar('press_sigma', data%Lev2b%press_sigma, & units = data%Lev2b%units%press_sigma, & range = data%Lev2b%range%press_sigma, & rec = irec) CALL ncdf_getvar('temp', data%Lev2b%temp, & units = data%Lev2b%units%temp, & range = data%Lev2b%range%temp, & rec = irec) CALL ncdf_getvar('temp_sigma', data%Lev2b%temp_sigma, & units = data%Lev2b%units%temp_sigma, & range = data%Lev2b%range%temp_sigma, & rec = irec) CALL ncdf_getvar('shum', data%Lev2b%shum, & units = data%Lev2b%units%shum, & range = data%Lev2b%range%shum, & rec = irec) CALL ncdf_getvar('shum_sigma', data%Lev2b%shum_sigma, & units = data%Lev2b%units%shum_sigma, & range = data%Lev2b%range%shum_sigma, & rec = irec) CALL ncdf_getvar('meteo_qual', data%Lev2b%meteo_qual, & units = data%Lev2b%units%meteo_qual, & range = data%Lev2b%range%meteo_qual, & rec = irec) ENDIF ! 2.14 Level2c variables (if any) ! ------------------------------- IF (ncdf_isvar('geop_sfc')) THEN n2d(1) = 1 n2d(2) = 0 CALL ncdf_getsize('geop_sfc', n2d(2), dim = 1) CALL ropp_io_init(data%Lev2c, n2d) ELSE data%Lev2c%Npoints = 0 data%Lev2c%Nhoriz = 0 ENDIF IF (data%Lev2c%Npoints > 0) THEN ! new 2d code CALL ncdf_getvar('dtheta', data%Lev2c%dtheta, & units = data%Lev2c%units%dtheta, & range = data%Lev2c%range%dtheta, & rec = irec) CALL ncdf_getvar('lat_2d', data%Lev2c%lat_2d, & units = data%Lev2c%units%lat_2d, & range = data%Lev2c%range%lat_2d, & rec = irec) CALL ncdf_getvar('lon_2d', data%Lev2c%lon_2d, & units = data%Lev2c%units%lon_2d, & range = data%Lev2c%range%lon_2d, & rec = irec) CALL ncdf_getvar('geop_sfc', data%Lev2c%geop_sfc, & units = data%Lev2c%units%geop_sfc, & range = data%Lev2c%range%geop_sfc, & rec = irec) CALL ncdf_getvar('press_sfc', data%Lev2c%press_sfc, & units = data%Lev2c%units%press_sfc, & range = data%Lev2c%range%press_sfc, & rec = irec) CALL ncdf_getvar('press_sfc_sigma', data%Lev2c%press_sfc_sigma, & units = data%Lev2c%units%press_sfc_sigma, & range = data%Lev2c%range%press_sfc_sigma, & rec = irec) CALL ncdf_getvar('press_sfc_qual', data%Lev2c%press_sfc_qual, & units = data%Lev2c%units%press_sfc_qual, & range = data%Lev2c%range%press_sfc_qual, & rec = irec) ENDIF ! 2.15 Level2d variables (if any) ! ------------------------------- IF (ncdf_isvar('level_coeff_a')) THEN CALL ncdf_getsize('level_coeff_a', n, dim = 1) CALL ropp_io_init(data%Lev2d, n) ELSE data%Lev2d%Npoints = 0 ENDIF IF (data%Lev2d%Npoints > 0) THEN CALL ncdf_getvar('level_type', data%Lev2d%level_type, & rec = irec) CALL ncdf_getvar('level_coeff_a', data%Lev2d%level_coeff_a, & units = data%Lev2d%units%level_coeff_a, & range = data%Lev2d%range%level_coeff_a, & rec = irec) CALL ncdf_getvar('level_coeff_b', data%Lev2d%level_coeff_b, & units = data%Lev2d%units%level_coeff_b, & range = data%Lev2d%range%level_coeff_b, & rec = irec) ENDIF ! 2.16 Additional variables (if any) ! ---------------------------------- CALL ropp_io_init(data%vlist) DO varid=1,ncdf_nvars IF(.NOT. ncdf_read(varid))THEN status = nf90_inquire_variable(ncdf_ncid, varid, xtype=TYPE, ndims=ndim) IF (TYPE .NE. NF90_CHAR) THEN ! only read scalar variables IF(ndim == 1)THEN CALL ropp_io_read_ncdf_get_vlistD0d(varid, data%vlist%VlistD0d, irec) ENDIF IF(ndim == 2)THEN CALL ropp_io_read_ncdf_get_vlistD1d(varid, data%vlist%VlistD1d, irec) ENDIF IF(ndim == 3)THEN CALL ropp_io_read_ncdf_get_vlistD2d(varid, data%vlist%VlistD2d, irec) ENDIF ENDIF ENDIF ncdf_read(varid) = .FALSE. ! reset 'read variable' flag ENDDO ! 2.17 Clean up ! ------------- CALL message_set_routine(routine) END SUBROUTINE ropp_io_read_ncdf_get_rodata_2d !------------------------------------------------------------------------------- ! 3. Error correlation and covariance matrices !------------------------------------------------------------------------------- SUBROUTINE ropp_io_read_ncdf_get_rocorcov(DATA) ! 3.1 Declarations ! ---------------- USE ropp_utils USE ncdf USE ropp_io, not_this => ropp_io_read_ncdf_get_rocorcov USE ropp_io_types, ONLY: ROcorcov IMPLICIT NONE TYPE(ROcorcov), DIMENSION(:), POINTER :: DATA REAL(wp), DIMENSION(:), POINTER :: lat_min => null() REAL(wp), DIMENSION(:), POINTER :: lat_max => null() REAL(wp), DIMENSION(:,:), POINTER :: corr => null() REAL(wp), DIMENSION(:,:), POINTER :: sigma => null() REAL(wp), DIMENSION(:,:), POINTER :: temp_sigma => null() REAL(wp), DIMENSION(:,:), POINTER :: shum_rel_sigma => null() REAL(wp), DIMENSION(:), POINTER :: press_sfc_rel_sigma => null() INTEGER :: i, m, n CHARACTER(len = 256) :: routine ! 3.2 Error handling ! ------------------ CALL message_get_routine(routine) CALL message_set_routine('ropp_io_read_ncdf_get') ! 3.3 Latitude bins ! ----------------- IF (ncdf_isvar('lat_min') .AND. ncdf_isvar('lat_max')) THEN CALL ncdf_getsize('lat_min', m) CALL ncdf_getsize('lat_max', n) ELSE CALL message(msg_fatal, & "NetCDF data file does not seem to contain an error correlation or covariance structure.") ENDIF IF (m /= m) THEN CALL message(msg_fatal, & "Number of latitude bin boundaries in the netCDF data file is inconsistent.") ENDIF CALL ncdf_getvar_alloc('lat_min', lat_min) CALL ncdf_getvar_alloc('lat_max', lat_max) ! 3.4 Error correlation matrices ! ------------------------------ IF (ncdf_isvar('corr')) THEN CALL ncdf_getvar_alloc('corr', corr) ELSE CALL message(msg_fatal, & "NetCDF data file does not seem to contain an error correlation matrix.") ENDIF ! 3.5 Error standard deviations ! ----------------------------- IF (ncdf_isvar('sigma')) THEN CALL ncdf_getvar_alloc('sigma', sigma) ENDIF IF (ncdf_isvar('temp_sigma')) THEN CALL ncdf_getvar_alloc('temp_sigma', temp_sigma) ENDIF IF (ncdf_isvar('shum_rel_sigma')) THEN CALL ncdf_getvar_alloc('shum_rel_sigma', shum_rel_sigma) ENDIF IF (ncdf_isvar('press_sfc_rel_sigma')) THEN CALL ncdf_getvar_alloc('press_sfc_rel_sigma', press_sfc_rel_sigma) ENDIF ! 3.6 Allocate and fill ROPP structure ! ------------------------------------ ALLOCATE(DATA(n)) DO i = 1, n DATA(i)%lat_min = lat_min(i) DATA(i)%lat_max = lat_max(i) ALLOCATE(DATA(i)%corr(SIZE(corr(:,i), 1))) DATA(i)%corr = corr(:,i) IF (ASSOCIATED(sigma)) THEN ALLOCATE(DATA(i)%sigma(SIZE(sigma(:,i), 1))) DATA(i)%sigma = sigma(:, i) ENDIF IF (ASSOCIATED(temp_sigma)) THEN ALLOCATE(DATA(i)%temp_sigma(SIZE(temp_sigma(:,i), 1))) DATA(i)%temp_sigma = temp_sigma(:, i) ENDIF IF (ASSOCIATED(shum_rel_sigma)) THEN ALLOCATE(DATA(i)%shum_rel_sigma(SIZE(shum_rel_sigma(:,i), 1))) DATA(i)%shum_rel_sigma = shum_rel_sigma(:, i) ENDIF IF (ASSOCIATED(press_sfc_rel_sigma)) THEN DATA(i)%press_sfc_rel_sigma = press_sfc_rel_sigma(i) ENDIF ENDDO ! 3.7 Clean up ! ------------ DEALLOCATE(lat_min) DEALLOCATE(lat_max) DEALLOCATE(corr) IF (ASSOCIATED(sigma)) DEALLOCATE(sigma) CALL message_set_routine(routine) END SUBROUTINE ropp_io_read_ncdf_get_rocorcov !------------------------------------------------------------------------------- ! 4. Wrapper for other centres' RO data !------------------------------------------------------------------------------- SUBROUTINE ropp_io_read_ncdf_get_otherdata(DATA, file, centre, rec, & getlevel1b, getextrap, & resolution, getlevel1a, & getdirect, getiono) ! 4.1 Declarations ! ---------------- USE ropp_io, not_this => ropp_io_read_ncdf_get_otherdata USE ropp_io_types, ONLY: ROprof IMPLICIT NONE TYPE(ROprof), INTENT(inout) :: DATA CHARACTER (len = *), INTENT(in) :: file CHARACTER (len=20), INTENT(in) :: centre INTEGER, OPTIONAL :: rec LOGICAL, OPTIONAL :: getlevel1b LOGICAL, OPTIONAL :: getextrap CHARACTER (len=20), OPTIONAL :: resolution CHARACTER (len=20), OPTIONAL :: getlevel1a LOGICAL, OPTIONAL :: getdirect LOGICAL, OPTIONAL :: getiono LOGICAL :: lgetlevel1b = .TRUE. LOGICAL :: lgetextrap = .FALSE. CHARACTER (len=20) :: lresolution = 'thinned' CHARACTER (len=20) :: lgetlevel1a = 'none' LOGICAL :: lgetdirect = .FALSE. LOGICAL :: lgetiono = .FALSE. LOGICAL :: ldummy = .FALSE. ! defaults IF ( PRESENT(getlevel1b) ) lgetlevel1b=getlevel1b IF ( PRESENT(getextrap) ) lgetextrap=getextrap IF ( PRESENT(resolution) ) lresolution=resolution IF ( PRESENT(getlevel1a) ) lgetlevel1a=getlevel1a IF ( PRESENT(getdirect) ) lgetdirect=getdirect IF ( PRESENT(getiono) ) lgetiono=getiono ! call the appropriate data handling function, default is ROPP format SELECT CASE (centre) CASE('UCAR') CALL ropp_io_read_ncdf_get_ucardata(DATA, file) CASE('UCAR1') CALL ropp_io_read_ncdf_get_ucardata(DATA, file, oldform=.TRUE.) CASE('EUM') CALL ropp_io_read_ncdf_get_eumdata(DATA, file, lgetlevel1b, lgetextrap, & lresolution, lgetlevel1a, lgetdirect, & lgetiono,ldummy) write(*, *) 'data%Lev1a%Npoints after ropp_io_read_ncdf_get_eumdata', DATA%Lev1a%Npoints CASE default CALL ropp_io_read_ncdf_get_rodata(DATA, rec) END SELECT END SUBROUTINE ropp_io_read_ncdf_get_otherdata !------------------------------------------------------------------------------- ! 5. UCAR RO data !------------------------------------------------------------------------------- SUBROUTINE ropp_io_read_ncdf_get_ucardata(DATA, file, oldform) ! 5.1 Declarations ! ---------------- USE ncdf USE ropp_utils USE ropp_io_types, ONLY: ROprof IMPLICIT NONE TYPE(ROprof), INTENT(inout) :: DATA CHARACTER(len = *), INTENT(in) :: file LOGICAL, OPTIONAL, INTENT(in) :: oldform CHARACTER(len = 256) :: routine LOGICAL :: xoldform ! 5.2 Error handling ! ------------------ CALL message_get_routine(routine) CALL message_set_routine('ropp_io_read_ucardata') xoldform = .FALSE. IF ( PRESENT(oldform) ) xoldform = oldform ! 5.3 Identify file type ! ---------------------- IF ( ncdf_isvar('Bend_ang') .AND. & (ncdf_isvar('Impact_parm') .OR. ncdf_isvar('Impact_height')) ) THEN CALL ropp_io_read_ucardata_atmPrf(DATA, file, old=xoldform) ELSE IF ( ncdf_isvar('pL1Snr') .AND. ncdf_isvar('pL2Snr') ) THEN CALL ropp_io_read_ucardata_atmPhs(DATA, file) ELSE IF ( ncdf_isvar('MSL_alt') .AND. ncdf_isvar('Pres') ) THEN CALL ropp_io_read_ucardata_atmPrf(DATA, file, old=xoldform) ELSE CALL message(msg_fatal, & "Routine ropp_io_read_ncdf_get_ucardata does not support this" // & "file type. Only atmPrf, ecmPrf, gfsPrf, ncpPrf, sonPrf and atmPhs "// & "files supported. Check input file type.") ENDIF CALL message_set_routine(routine) CONTAINS !------------------------------------------------------------------------------- ! 6. UCAR RO data - atmPrf files !------------------------------------------------------------------------------- SUBROUTINE ropp_io_read_ucardata_atmPrf(DATA, file, old) ! NB this routine only supports UCAR 'Prf' format netCDF files ! (i.e. atmPrf, ecmPrf, gfsPrf, ncpPrf, sonPrf and wetPrf) ! See: http://cdaac-www.cosmic.ucar.edu/cdaac/cgi_bin/fileFormats.cgi?type=atmPrf ! Keyword old ==> use the old definition of impact_height = impact_param - RoC. ! By default from ROPP10.0, the code uses the new definition of ! impact_height = impact_param - RoC - und. ! 6.1 Declarations ! ---------------- USE DateTimeProgs, ONLY: Date_and_Time_UTC USE DateTimeTypes USE ropp_utils USE ncdf USE ropp_io USE ropp_io_types, ONLY: ROprof, & ThisFmtVer, & PCD_summary, & PCD_bangle, & PCD_refrac, & PCD_open_loop, & PCD_rising, & PCD_occultation USE geodesy, ONLY: geometric2geopotential IMPLICIT NONE TYPE(ROprof), INTENT(inout) :: DATA CHARACTER(len = *), INTENT(in) :: file LOGICAL, OPTIONAL, INTENT(in) :: old INTEGER :: n INTEGER :: readint REAL(wp) :: readreal CHARACTER (len = 256) :: readstr REAL(wp), PARAMETER :: freezing_point_in_K = 273.15_wp REAL(wp), PARAMETER :: epsilon_water = 0.621971_wp INTEGER, DIMENSION(8) :: DTnow INTEGER :: nidx INTEGER :: xtype LOGICAL :: lev1b_present, lev2a_present, lev2b_present LOGICAL :: xold ! Local version of old ! 6.2 Interrogate file ! -------------------- xold = .FALSE. IF ( PRESENT(old) ) xold = old lev1b_present = .FALSE. IF( ncdf_isvar('Bend_ang') .AND. & (ncdf_isvar('Impact_parm') .OR. ncdf_isvar('Impact_height')) ) & lev1b_present = .TRUE. lev2a_present = .FALSE. IF ( ncdf_isvar('MSL_alt') .AND. ncdf_isvar('Ref') ) & lev2a_present = .TRUE. lev2b_present = .FALSE. IF ( ncdf_isvar('MSL_alt') .AND. ncdf_isvar('Pres') ) & lev2b_present = .TRUE. ! 6.3 Header variables ! -------------------- ! NB: fileStamp = 'IIII.YYYY.DDD.HH.MM.GGG' readstr = ' ' CALL ncdf_getatt('fileStamp', readstr) data%leo_id = readstr(1:4) IF ( data%leo_id(1:2) == 'CH' ) THEN data%leo_id = 'CHMP' ELSE IF ( data%leo_id(1:4) == 'GRC1' ) THEN data%leo_id = 'GRAA' ELSE IF ( data%leo_id(1:4) == 'GRC2' ) THEN data%leo_id = 'GRAB' ELSE IF ( data%leo_id(1:2) == 'CO' ) THEN data%leo_id(1:2) = 'C0' ELSE IF ( data%leo_id(1:3) == 'MTP' ) THEN data%leo_id(1:3) = 'MET' ELSE IF ( data%leo_id(1:4) == 'G001' ) THEN ! subject to change data%leo_id = 'FY3C' ELSE IF ( data%leo_id(1:4) == 'G002' ) THEN ! subject to change data%leo_id = 'FY3D' END IF READ ( readstr(22:), FMT='(I3.3)' ) readint WRITE ( data%gns_id, FMT='(A1,I3.3)' ) readstr(21:21), readint IF ( ncdf_isatt('fiducial_id') ) THEN readstr = ' ' CALL ncdf_getatt('fiducial_id', readstr) IF(readstr /= " ") data%stn_id = readstr(1:4) END IF readstr = ' ' IF ( ncdf_isatt('occfreq1') ) THEN CALL ncdf_getatt('occfreq1', data%occfreq1) CALL ncdf_getatt('occfreq2', data%occfreq2) ENDIF ! 6.4 PCD flags and overall quality ! --------------------------------- IF ( lev1b_present ) THEN data%PCD = 0 ! All PCD bits unset data%units%overall_qual = "%" CALL ncdf_getatt('bad', readstr) IF ( TRIM(readstr) == "1" ) THEN data%PCD = IBSET(data%PCD, PCD_summary) ! PCD bit 1 set: non-nominal quality data%PCD = IBSET(data%PCD, PCD_bangle) ! PCD bit 5 set: non-nominal bangle data%PCD = IBSET(data%PCD, PCD_refrac) ! PCD bit 6 set: non-nominal refrac data%overall_qual = 0.0 ELSE data%overall_qual = 100.0 END IF IF ( ncdf_isatt('iol') ) THEN CALL ncdf_getatt('iol', readint) IF ( readint == 1 ) THEN data%PCD = IBSET(data%PCD, PCD_open_loop) ! PCD bit 8 set: open loop used END IF END IF IF ( ncdf_isatt('irs') ) THEN CALL ncdf_getatt('irs', readint) IF ( readint == -1 ) THEN data%PCD = IBSET(data%PCD, PCD_rising) ! PCD bit 3 set: rising occultation END IF END IF ELSE data%PCD = 0 ! All PCD bits unset data%PCD = IBSET(data%PCD, PCD_occultation) ! PCD bit 15 set: background data END IF ! 6.5 Date and time ! ----------------- CALL gettime(data%DTocc) ! 6.6 Georeferencing ! ------------------ IF ( ncdf_isatt('lat', xtype=xtype) ) THEN IF ( xtype > NF90_CHAR ) THEN CALL ncdf_getatt('lat', data%georef%lat) ELSE CALL ncdf_getatt('lat', readstr) READ (readstr, *) data%georef%lat ENDIF data%georef%units%lat = "degrees_north" ENDIF IF ( ncdf_isatt('lon', xtype=xtype) ) THEN IF ( xtype > NF90_CHAR ) THEN CALL ncdf_getatt('lon', data%georef%lon) ELSE CALL ncdf_getatt('lon', readstr) READ (readstr, *) data%georef%lon ENDIF data%georef%units%lon = "degrees_east" ENDIF IF ( lev1b_present ) THEN CALL ncdf_getatt('occpt_offset', data%georef%time_offset) data%georef%units%time_offset = "seconds" CALL ncdf_getatt('rgeoid', data%georef%Undulation) data%georef%Undulation = data%georef%Undulation * 1000.0 data%georef%units%Undulation = "meters" CALL ncdf_getatt('rfict', data%georef%roc) data%georef%roc = data%georef%roc * 1000.0 data%georef%units%roc = "meters" CALL ncdf_getatt('curv', data%georef%r_coc) data%georef%r_coc = data%georef%r_coc * 1000.0 data%georef%units%r_coc = "meters" CALL ncdf_getatt('azim', data%georef%azimuth) IF ((data%georef%azimuth < 0.0) .AND. & (data%georef%azimuth > -180.0)) & data%georef%azimuth = data%georef%azimuth + 360.0 data%georef%units%azimuth = "degrees" ENDIF ! 6.7 Level1b variables (if any) ! ------------------------------- IF ( lev1b_present ) THEN CALL ncdf_getsize('Bend_ang', n, dim = 1) CALL ropp_io_init(data%Lev1b, n) ELSE data%Lev1b%Npoints = 0 ENDIF IF (data%Lev1b%Npoints > 0) THEN CALL ncdf_getvar('Lat', data%Lev1b%lat_tp) data%Lev1b%units%lat_tp = "degrees_north" CALL ncdf_getvar('Lon', data%Lev1b%lon_tp) data%Lev1b%units%lon_tp = "degrees_east" CALL ncdf_getvar('Azim', data%Lev1b%azimuth_tp) WHERE ( (data%Lev1b%azimuth_tp < 0.0) .AND. & (data%Lev1b%azimuth_tp > -180.0) ) & data%Lev1b%azimuth_tp = data%Lev1b%azimuth_tp + 360.0 data%Lev1b%units%azimuth_tp = "degrees" IF (ncdf_isvar('Impact_parm')) THEN CALL ncdf_getvar('Impact_parm', data%Lev1b%impact, & units = data%Lev1b%units%impact) ELSE ! If Impact_parm is absent, derive it from Impact_height, if possible IF (ncdf_isvar('Impact_height')) THEN CALL ncdf_getvar('Impact_height', data%Lev1b%impact, & units = data%Lev1b%units%impact) IF ( xold ) THEN data%Lev1b%impact = data%Lev1b%impact + data%georef%roc ! Both are in m at this point ELSE data%Lev1b%impact = data%Lev1b%impact + data%georef%roc + & data%georef%undulation ! All are in m at this point END IF END IF END IF CALL ncdf_getvar('Bend_ang', data%Lev1b%bangle) data%Lev1b%units%bangle = "radians" CALL ncdf_getvar('Bend_ang_stdv', data%Lev1b%bangle_sigma) data%Lev1b%units%bangle_sigma = "radians" ! set the quality for bangle CALL ncdf_getatt('_FillValue', readreal, 'Bend_ang') WHERE (data%Lev1b%bangle > readreal) data%Lev1b%bangle_qual = 100.0 IF ( ncdf_isvar('Bend_ang1') ) THEN ! Leave the unknown sigmas and quals unset CALL ncdf_getvar('Bend_ang1', data%Lev1b%bangle_L1) data%Lev1b%impact_L1 = data%Lev1b%impact END IF IF ( ncdf_isvar('Bend_ang2') ) THEN ! Leave the unknown sigmas and quals unset CALL ncdf_getvar('Bend_ang2', data%Lev1b%bangle_L2) data%Lev1b%impact_L2 = data%Lev1b%impact END IF IF ( ncdf_isvar('Opt_bend_ang') ) THEN CALL ncdf_getvar('Opt_bend_ang', data%Lev1b%bangle_opt) data%lev1b%impact_opt = data%Lev1b%impact data%Lev1b%bangle_opt_sigma = data%Lev1b%bangle_sigma data%Lev1b%bangle_opt_qual = data%Lev1b%bangle_qual END IF ENDIF ! 6.8 Level2a variables (if any) ! ------------------------------- IF ( lev2a_present ) THEN CALL ncdf_getsize('MSL_alt', n, dim = 1) CALL ropp_io_init(data%Lev2a, n) ELSE data%Lev2a%Npoints = 0 ENDIF IF (data%Lev2a%Npoints > 0) THEN CALL ncdf_getvar('MSL_alt', data%Lev2a%alt_refrac, & units = data%Lev2a%units%alt_refrac) ! Geopotential not in UCAR files - generate from altitude IF (data%georef%lat > -999.0) THEN WHERE(data%Lev2a%alt_refrac > -999.0) & data%Lev2a%geop_refrac = geometric2geopotential(data%georef%lat, & data%Lev2a%alt_refrac) ENDIF CALL ncdf_getvar('Ref', data%Lev2a%refrac) data%Lev2a%units%refrac = "1" IF (ncdf_isvar('Ref_stdv')) THEN CALL ncdf_getvar('Ref_stdv', data%Lev2a%refrac_sigma) data%Lev2a%units%refrac_sigma = "1" ! set the quality for ref CALL ncdf_getatt('_FillValue', readreal, 'Ref' ) WHERE( data%Lev2a%refrac > readreal) & data%Lev2a%refrac_qual = 100.0 ENDIF ! Include dry temperature in Level 2a. If the variable 'Bend_ang' exists, ! the UCAR file is an atmPrf file, and 'Temp' is actually dry temperature. IF (ncdf_isvar('Temp') .AND. ncdf_isvar('Bend_ang')) THEN CALL ncdf_getvar('Temp', data%Lev2a%dry_temp) WHERE(data%Lev2a%dry_temp > -999.0) & data%Lev2a%dry_temp = data%Lev2a%dry_temp + freezing_point_in_K data%Lev2a%units%dry_temp = "kelvin" END IF ENDIF ! 6.9 Level2b variables (if any) ! ------------------------------- IF ( lev2b_present ) THEN CALL ncdf_getsize('MSL_alt', n, dim = 1) CALL ropp_io_init(data%Lev2b, n) ELSE data%Lev2b%Npoints = 0 ENDIF IF (data%Lev2b%Npoints > 0) THEN CALL ncdf_getvar('Pres', data%Lev2b%press) data%Lev2b%units%press = "hPa" ! If variable 'Bend_ang' exists, the UCAR file is an atmPrf file, ! and the 'Temp' variable is actually _dry_ temperature, so don't populate ! the level 2b temperature with it. IF ( .NOT. ncdf_isvar('Bend_ang') ) THEN CALL ncdf_getvar('Temp', data%Lev2b%temp) WHERE (data%Lev2b%temp > -999.0) & data%Lev2b%temp = data%Lev2b%temp + freezing_point_in_K data%Lev2b%units%temp = "kelvin" ENDIF IF (ncdf_isvar('Vp')) THEN CALL ncdf_getvar('Vp', data%Lev2b%shum) WHERE (data%Lev2b%shum > -999.0) & data%Lev2b%shum = & (data%Lev2b%shum*epsilon_water) / & (data%Lev2b%press - & (data%Lev2b%shum*(1.0_wp - epsilon_water))) data%Lev2b%units%shum = "kilogram/kilogram" ENDIF ! Get the geopotential from L2a CALL ncdf_getvar('MSL_alt', data%Lev2b%geop, units = 'metres') data%Lev2b%geop = geometric2geopotential(data%georef%lat,data%Lev2b%geop) ! Background time data%bg%Year = data%DTocc%Year data%bg%Month = data%DTocc%Month data%bg%Day = data%DTocc%Day data%bg%Hour = data%DTocc%Hour data%bg%Minute = data%DTocc%Minute nidx = INDEX(file, 'Prf_', .TRUE.) IF ( nidx > 3 ) data%bg%source = file(nidx-3:nidx-1) ! read from filename ENDIF ! 6.10 (Global) Attributes ! ------------------------ data%FmtVersion = ' ' ; data%FmtVersion = ThisFmtVer data%processing_centre = ' ' ; CALL ncdf_getatt('center', data%processing_centre) data%pod_method = ' ' ; data%pod_method = "UNKNOWN" data%phase_method = ' ' ; data%phase_method = "UNKNOWN" data%bangle_method = ' ' ; data%bangle_method = "UNKNOWN" data%refrac_method = ' ' ; data%refrac_method = "UNKNOWN" data%meteo_method = ' ' ; data%meteo_method = "UNKNOWN" data%thin_method = ' ' ; data%thin_method = "UNKNOWN" data%software_version = ' ' ; data%software_version = "UNKNOWN" data%bad = ' ' ; data%bad = "0" nidx = INDEX(file, 'Prf_', .TRUE.) IF ( nidx > 0 ) & ! VNN.nnn from file name data%software_version = "V" // file(nidx+35:nidx+36) // & "." // file(nidx+29:nidx+31) data%processing_software = ' ' ; data%processing_software = "UNKNOWN" ! Set equal to SSSS_VVVV if the filename is ???Prf_IIII.YYYY.DDD.HH.MM.GGG_SSSS.VVVV_nc nidx = INDEX(file, '_nc', .TRUE.) IF ( nidx > 9 ) data%processing_software = file(nidx-9:nidx-1) ! COSMIC has no processing date, set to current utc date/time CALL Date_and_Time_UTC ( Values=DTnow ) data%DTpro%Year = DTnow(1) data%DTpro%Month = DTnow(2) data%DTpro%Day = DTnow(3) data%DTpro%Hour = DTnow(5) data%DTpro%Minute = DTnow(6) data%DTpro%Second = DTnow(7) data%DTpro%Msec = DTnow(8) ! 6.11 Occultation ID ! ------------------- CALL ropp_io_occid(DATA) END SUBROUTINE ropp_io_read_ucardata_atmPrf !------------------------------------------------------------------------------- ! 7. UCAR RO data - atmPhs files !------------------------------------------------------------------------------- SUBROUTINE ropp_io_read_ucardata_atmPhs(DATA, file) ! NB this routine only supports UCAR atmPhs netCDF files. ! See: http://cdaac-www.cosmic.ucar.edu/cdaac/cgi_bin/fileFormats.cgi?type=atmPhs ! 7.1 Declarations ! ---------------- USE DateTimeProgs, ONLY: Date_and_Time_UTC USE DateTimeTypes USE ropp_utils USE ncdf USE ropp_io USE ropp_io_types, ONLY: ROprof, & ThisFmtVer, & PCD_summary, & PCD_offline, & PCD_phase, & PCD_rising IMPLICIT NONE TYPE(ROprof), INTENT(inout) :: DATA CHARACTER(len = *), INTENT(in) :: file INTEGER :: n INTEGER :: nn INTEGER :: readint REAL(wp) :: readreal CHARACTER (len = 256) :: readstr INTEGER, DIMENSION(8) :: DTnow REAL(wp), DIMENSION(:), ALLOCATABLE :: workdata REAL(wp), DIMENSION(:), ALLOCATABLE :: startTime REAL(wp), DIMENSION(:), ALLOCATABLE :: fulltime REAL(wp), DIMENSION(:), ALLOCATABLE :: orbtime REAL(wp), DIMENSION(:), ALLOCATABLE :: txmitLR REAL(wp), DIMENSION(:), ALLOCATABLE :: txmitHR REAL(wp), DIMENSION(:), ALLOCATABLE :: xLR REAL(wp), DIMENSION(:), ALLOCATABLE :: yLR REAL(wp), DIMENSION(:), ALLOCATABLE :: zLR REAL(wp), DIMENSION(:), ALLOCATABLE :: xHR REAL(wp), DIMENSION(:), ALLOCATABLE :: yHR REAL(wp), DIMENSION(:), ALLOCATABLE :: zHR REAL(wp), DIMENSION(:), ALLOCATABLE :: d2txmit ! 2nd derivative ln(N) REAL(wp), DIMENSION(:), ALLOCATABLE :: d2xLR ! 2nd derivative ln(N) REAL(wp), DIMENSION(:), ALLOCATABLE :: d2yLR ! 2nd derivative ln(N) REAL(wp), DIMENSION(:), ALLOCATABLE :: d2zLR ! 2nd derivative ln(N) INTEGER :: ii ! Index INTEGER :: slen INTEGER :: nidx ! 7.2 Header variables ! -------------------- ! fileStamp = 'IIII.YYYY.DDD.HH.MM.GGG' ! or: fileStamp = 'IIIII.YYYY.DDD.HH.MM.GGG' for spire ! or: fileStamp = 'GN05_202604101630_202604101633_G30_A07_c20260410175814' readstr = ' ' CALL ncdf_getatt('fileStamp', readstr) slen = LEN_TRIM(readstr) IF (slen == 23) THEN data%leo_id = readstr(1:4) data%gns_id = readstr(21:23) ELSE IF (slen == 24) THEN data%leo_id = 'S' // readstr(3:5) ! change FM126 to S126 data%gns_id = readstr(22:24) ELSE IF (slen == 54) THEN data%leo_id = readstr(1:4) data%gns_id = readstr(32:34) END IF print*, 'data%leo_id: ',data%leo_id print*, 'data%gns_id: ',data%gns_id readstr = ' ' data%stn_id = ' ' IF ( ncdf_isatt('occfreq1') ) THEN CALL ncdf_getatt('occfreq1', data%occfreq1) CALL ncdf_getatt('occfreq2', data%occfreq2) ENDIF write(*, *) data%occfreq1, data%occfreq2 ! 7.3 Date and time ! ----------------- CALL gettime(data%DTocc) write(*, *) data%DTocc%year ! 7.4 PCD flags and overall quality ! --------------------------------- data%PCD = 0 ! All PCD bits unset data%PCD = IBSET(data%PCD, PCD_offline) ! Assuming that UCAR atmPhs files are non-NRT CALL ncdf_getatt('bad', readstr) IF ( TRIM(readstr) == "1" ) THEN data%PCD = IBSET(data%PCD, PCD_summary) ! PCD bit 1 set: non-nominal quality data%PCD = IBSET(data%PCD, PCD_phase) ! PCD bit 4 set: non-nominal phase data%overall_qual = 0.0 ELSE data%overall_qual = 100.0 END IF data%units%overall_qual = "%" IF ( ncdf_isatt('setting') ) THEN CALL ncdf_getatt('setting', readint) IF ( readint == 0 ) THEN data%PCD = IBSET(data%PCD, PCD_rising) ! PCD bit 3 set: rising occultation ENDIF ENDIF ! 7.5 Level1a variables (if any) ! ------------------------------- IF (ncdf_isvar('time')) THEN CALL ncdf_getsize('time', n, dim = 1) write(*, *) n CALL ropp_io_init(data%Lev1a, n) ELSE data%Lev1a%Npoints = 0 ENDIF write (*, *) 'before level1a data:', data%Lev1a%Npoints IF (data%Lev1a%Npoints > 0) THEN CALL ncdf_getvar('time', data%Lev1a%dtime) CALL ncdf_getvar('caL1Snr', data%Lev1a%snr_L1ca) CALL ncdf_getvar('pL1Snr', data%Lev1a%snr_L1p) CALL ncdf_getvar('pL2Snr', data%Lev1a%snr_L2p) data%Lev1a%units%snr = "volt/volt" CALL ncdf_getvar('exL1', data%Lev1a%phase_L1) CALL ncdf_getvar('exL2', data%Lev1a%phase_L2) data%Lev1a%units%phase = "metres" ! Updated by Yong Chen on 03/20/2023 ! Working for both atmPhs and conPhs format ! for new format IF (ncdf_isvar('orbtime')) THEN !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! added by Xinjia, for new UCAR conPhs format CALL ncdf_getsize('orbtime', nn, dim = 1) ALLOCATE(orbtime(nn)) ALLOCATE(txmitLR(nn)) CALL ncdf_getvar('orbtime', orbtime) CALL ncdf_getvar('txmitLR', txmitLR) ALLOCATE(startTime(1)) CALL ncdf_getvar('startTime', startTime) ALLOCATE(fulltime(n)) fulltime(:) = startTime(1)+data%Lev1a%dtime ALLOCATE(xLR(nn)) ALLOCATE(yLR(nn)) ALLOCATE(zLR(nn)) ALLOCATE(d2xLR(nn)) ALLOCATE(d2yLR(nn)) ALLOCATE(d2zLR(nn)) ALLOCATE(xHR(n)) ALLOCATE(yHR(n)) ALLOCATE(zHR(n)) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! CALL ncdf_getvar('xLeoLR', xLR) CALL ncdf_getvar('yLeoLR', yLR) CALL ncdf_getvar('zLeoLR', zLR) CALL ropp_io_init_spline(orbtime(:), xLR(:), d2xLR) CALL ropp_io_init_spline(orbtime(:), yLR(:), d2yLR) CALL ropp_io_init_spline(orbtime(:), zLR(:), d2zLR) DO ii = 1, n CALL ropp_io_interpol_spline(orbtime, xLR, d2xLR, & fulltime(ii), xHR(ii)) CALL ropp_io_interpol_spline(orbtime, yLR, d2yLR, & fulltime(ii), yHR(ii)) CALL ropp_io_interpol_spline(orbtime, zLR, d2zLR, & fulltime(ii), zHR(ii)) ENDDO data%Lev1a%r_leo(:,1) = xHR(:) data%Lev1a%r_leo(:,2) = yHR(:) data%Lev1a%r_leo(:,3) = zHR(:) data%Lev1a%r_leo(:,:) = data%Lev1a%r_leo(:,:) * 1000.0_wp data%Lev1a%units%r_leo = "metres" data%Lev1a%reference_frame%r_leo = "ECI" DEALLOCATE(xLR) DEALLOCATE(yLR) DEALLOCATE(zLR) DEALLOCATE(d2xLR) DEALLOCATE(d2yLR) DEALLOCATE(d2zLR) DEALLOCATE(xHR) DEALLOCATE(yHR) DEALLOCATE(zHR) ! updating stopped here !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !data%Lev1a%v_leo(:,:) = data%Lev1a%r_leo(:,:) data%Lev1a%units%v_leo = "metres / seconds" !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ALLOCATE(xLR(nn)) ALLOCATE(yLR(nn)) ALLOCATE(zLR(nn)) ALLOCATE(d2xLR(nn)) ALLOCATE(d2yLR(nn)) ALLOCATE(d2zLR(nn)) ALLOCATE(xHR(n)) ALLOCATE(yHR(n)) ALLOCATE(zHR(n)) ALLOCATE(txmitHR(n)) CALL ropp_io_interpol(orbtime, fulltime, txmitLR, txmitHR) CALL ncdf_getvar('xGnssLR', xLR) CALL ncdf_getvar('yGnssLR', yLR) CALL ncdf_getvar('zGnssLR', zLR) CALL ropp_io_init_spline(txmitLR(:), xLR(:), d2xLR) CALL ropp_io_init_spline(txmitLR(:), yLR(:), d2yLR) CALL ropp_io_init_spline(txmitLR(:), zLR(:), d2zLR) DO ii = 1, n CALL ropp_io_interpol_spline(txmitLR, xLR, d2xLR, & txmitHR(ii), xHR(ii)) CALL ropp_io_interpol_spline(txmitLR, yLR, d2yLR, & txmitHR(ii), yHR(ii)) CALL ropp_io_interpol_spline(txmitLR, zLR, d2zLR, & txmitHR(ii), zHR(ii)) ! print *, ii, xHR(ii), yHR(ii), zHR(ii) ENDDO data%Lev1a%r_gns(:,1) = xHR(:) data%Lev1a%r_gns(:,2) = yHR(:) data%Lev1a%r_gns(:,3) = zHR(:) data%Lev1a%r_gns(:,:) = data%Lev1a%r_gns(:,:) * 1000.0_wp data%Lev1a%units%r_gns = "metres" data%Lev1a%reference_frame%r_gns = "ECI" ! print *, 'data%Lev1a%r_gns*****************' ! print *, data%Lev1a%r_gns DEALLOCATE(xLR) DEALLOCATE(yLR) DEALLOCATE(zLR) DEALLOCATE(d2xLR) DEALLOCATE(d2yLR) DEALLOCATE(d2zLR) DEALLOCATE(xHR) DEALLOCATE(yHR) DEALLOCATE(zHR) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !data%Lev1a%v_gns(:,:) = data%Lev1a%r_gns(:,:) data%Lev1a%units%v_gns = "metres / seconds" ELSE ! for old foramt CALL ncdf_getvar('xLeo', data%Lev1a%r_leo(:,1)) CALL ncdf_getvar('yLeo', data%Lev1a%r_leo(:,2)) CALL ncdf_getvar('zLeo', data%Lev1a%r_leo(:,3)) data%Lev1a%r_leo(:,:) = data%Lev1a%r_leo(:,:) * 1000.0_wp data%Lev1a%units%r_leo = "metres" data%Lev1a%reference_frame%r_leo = "ECI" CALL ncdf_getvar('xdLeo', data%Lev1a%v_leo(:,1)) CALL ncdf_getvar('ydLeo', data%Lev1a%v_leo(:,2)) CALL ncdf_getvar('zdLeo', data%Lev1a%v_leo(:,3)) data%Lev1a%v_leo(:,:) = data%Lev1a%v_leo(:,:) * 1000.0_wp data%Lev1a%units%v_leo = "metres / seconds" CALL ncdf_getvar('xGps', data%Lev1a%r_gns(:,1)) CALL ncdf_getvar('yGps', data%Lev1a%r_gns(:,2)) CALL ncdf_getvar('zGps', data%Lev1a%r_gns(:,3)) data%Lev1a%r_gns(:,:) = data%Lev1a%r_gns(:,:) * 1000.0_wp data%Lev1a%units%r_gns = "metres" data%Lev1a%reference_frame%r_gns = "ECI" CALL ncdf_getvar('xdGps', data%Lev1a%v_gns(:,1)) CALL ncdf_getvar('ydGps', data%Lev1a%v_gns(:,2)) CALL ncdf_getvar('zdGps', data%Lev1a%v_gns(:,3)) data%Lev1a%v_gns(:,:) = data%Lev1a%v_gns(:,:) * 1000.0_wp data%Lev1a%units%v_gns = "metres / seconds" ENDIF ! open loop phase model data ALLOCATE(workdata(data%Lev1a%Npoints)) SELECT CASE ( data%leo_id(1:2) ) CASE ('ME','KO','PA','TS','TD','C2','GO','C0','S0','S1','S2', 'GN', 'YM') IF ( ncdf_isvar('xmdldd') ) THEN CALL ncdf_getvar('xmdldd', workdata) CALL ropp_io_addvar_rodataD1d(DATA, & name = "open_loop_lcf", & long_name= "Lost Carrier Flag", & units = "metres", & range = (/-1.0E7_wp, 1.0E7_wp/), & DATA = workdata) ENDIF CASE DEFAULT IF ( ncdf_isvar('xmdl') ) THEN CALL ncdf_getvar('xmdl', workdata) CALL ropp_io_addvar_rodataD1d(DATA, & name = "open_loop_lcf", & long_name= "Lost Carrier Flag", & units = "metres", & range = (/-1.0E7_wp, 1.0E7_wp/), & DATA = workdata) ENDIF END SELECT DEALLOCATE(workdata) ! set the quality for phase CALL ncdf_getatt('_FillValue', readreal, 'exL1') WHERE (data%Lev1a%phase_L1 > readreal) & data%Lev1a%phase_qual = 100.0 ENDIF write(*, *) 'end level1a data' ! 7.6 (Global) Attributes ! ------------------------ data%FmtVersion = ' ' ; data%FmtVersion = ThisFmtVer data%processing_centre = ' ' ; CALL ncdf_getatt('center', data%processing_centre) data%pod_method = ' ' ; data%pod_method = "UNKNOWN" data%phase_method = ' ' ; data%phase_method = "UNKNOWN" data%bangle_method = ' ' ; data%bangle_method = "UNKNOWN" data%refrac_method = ' ' ; data%refrac_method = "UNKNOWN" data%meteo_method = ' ' ; data%meteo_method = "UNKNOWN" data%thin_method = ' ' ; data%thin_method = "UNKNOWN" data%software_version = ' ' ; data%software_version = "UNKNOWN" data%bad = ' ' ; data%bad = "0" nidx = INDEX(file, 'atmPhs_', .TRUE.) IF ( nidx > 0 ) & ! VNN.nnn from file name data%software_version = "V" // file(nidx+38:nidx+39) // & "." // file(nidx+32:nidx+34) data%processing_software = ' ' ; data%processing_software = "UNKNOWN" ! Set equal to SSSS_VVVV if the filename is atmPhs_IIII.YYYY.DDD.HH.MM.GGG_SSSS.VVVV_nc nidx = INDEX(file, '_nc', .TRUE.) IF ( nidx > 9 ) data%processing_software = file(nidx-9:nidx-1) ! COSMIC has no processing date, set to current utc date/time CALL Date_and_Time_UTC ( Values=DTnow ) data%DTpro%Year = DTnow(1) data%DTpro%Month = DTnow(2) data%DTpro%Day = DTnow(3) data%DTpro%Hour = DTnow(5) data%DTpro%Minute = DTnow(6) data%DTpro%Second = DTnow(7) data%DTpro%Msec = DTnow(8) ! 7.7 Occultation ID ! ------------------- CALL ropp_io_occid(DATA) ! 7.8 Clean up ! ------------- CALL message_set_routine(routine) write (*,*) 'end of ropp_io_read_ucardata_atmPhs' END SUBROUTINE ropp_io_read_ucardata_atmPhs !------------------------------------------------------------------------------- ! 15. Splines: essentially a copy of the ropp_io_spline module !------------------------------------------------------------------------------- ! 15.1 Generate second derivative of spline ! ----------------------------------------- SUBROUTINE ropp_io_interpol(x, newx, array, interp) IMPLICIT NONE REAL(wp), DIMENSION(:), INTENT(in) :: x REAL(wp), DIMENSION(:), INTENT(in) :: newx REAL(wp), DIMENSION(:), INTENT(in) :: array REAL(wp), DIMENSION(:), INTENT(out) :: interp INTEGER :: i, j, k DO k = 1, SIZE(newx) j = 2 DO WHILE (j < SIZE(x) .AND. x(j) < newx(k)) j = j + 1 ENDDO i = j - 1 interp(k) = array(j) + & ( (newx(k) - x(j)) / (x(i) - x(j)) * (array(i) - array(j)) ) ENDDO END SUBROUTINE ropp_io_interpol ! 15.2 Generate second derivative of spline ! ----------------------------------------- SUBROUTINE ropp_io_init_spline(x, f, d2) IMPLICIT NONE ! 15.2.1 Declarations REAL(wp), DIMENSION(:), INTENT(in) :: x ! Argument grid (monotonous) REAL(wp), DIMENSION(:), INTENT(in) :: f ! Gridded function REAL(wp), DIMENSION(:), INTENT(out) :: d2 ! 2nd derivative of spline REAL(wp), DIMENSION(:), ALLOCATABLE :: d1 REAL(wp) :: dfl, dfr, df, a, b, c INTEGER :: i, N ! 15.2.2 Initialisation N = SIZE(x) ALLOCATE(d1(N)) d2(:) = 0.0_wp ! 2.3 Drive-through calculation of spline coefficients d1(1) = 0.0_wp d2(N) = 0.0_wp d1(1) = 0.0_wp DO i = 2, N-1 dfl = (f(i) - f(i-1))/(x(i) - x(i-1)) dfr = (f(i+1) - f(i))/(x(i+1) - x(i)) df = (dfr - dfl)/(x(i+1) - x(i-1)) a = (x(i) - x(i-1))/(2*(x(i+1) - x(i-1))) b = (x(i+1) - x(i))/(2*(x(i+1) - x(i-1))) c = 1 + a*d1(i-1) d1(i) = -b/c d2(i) = (3*df - a*d2(i-1))/c ENDDO DO i = N-1, 2, -1 d2(i) = d1(i)*d2(i+1) + d2(i) ENDDO DEALLOCATE(d1) END SUBROUTINE ropp_io_init_spline ! 15.3 Return interpolated function and first two derivatives ! ----------------------------------------------------------- SUBROUTINE ropp_io_interpol_spline(x, f, d2, x_int, f_int, fd_int, fd2_int) IMPLICIT NONE ! 15.3.1 Declarations REAL(wp), DIMENSION(:), INTENT(in) :: x ! Argument grid (monotonous) REAL(wp), DIMENSION(:), INTENT(in) :: f ! Gridded function REAL(wp), DIMENSION(:), INTENT(in) :: d2 ! 2nd derivative of spline REAL(wp), INTENT(in) :: x_int ! Interpolation point REAL(wp), INTENT(out) :: f_int ! Interpolated function value REAL(wp), OPTIONAL, INTENT(out) :: fd_int ! Interpolated 1st derivative REAL(wp), OPTIONAL, INTENT(out) :: fd2_int ! Interpolated 2nd deriv REAL(wp) :: a1, a2, a3 ! Polynomial coefficients REAL(wp) :: dx ! Grid interval REAL(wp) :: dx_t ! Grid-point to interpolation-point distance REAL(wp) :: x_t ! Interpolation point projected to grid extent REAL(wp) :: fd ! Interpolated derivative INTEGER :: i ! Array index INTEGER :: N ! Number of data INTEGER :: i_int ! Interpolation interval index ! 15.3.2 Location of interpolation point inside grid N = SIZE(x) x_t = MIN(MAX(x_int, MIN(x(1),x(N))), MAX(x(1),x(N))) i_int = ropp_io_seek_index(x, x_t) i = MAX(i_int, 1) ! 15.3.3 Calculation of interpolation coefficients dx = x(i+1) - x(i) a2 = d2(i)/2.0_wp a3 = (d2(i+1) - d2(i))/(6.0_wp*dx) a1 = (f(i+1) - f(i))/dx - dx*(a2 + dx*a3) ! 15.3.4 Calculated interpolated value dx_t = x_t - x(i) fd = a1 + dx_t*(2*a2 + dx_t*3*a3) f_int = f(i) + dx_t*(a1 + dx_t*(a2 + dx_t*a3)) + & fd*(x_int - x_t) IF (PRESENT(fd_int)) THEN fd_int = fd ENDIF IF (PRESENT(fd2_int)) THEN fd2_int = 2.0_wp*a2 + dx_t*6.0_wp*a3 ENDIF END SUBROUTINE ropp_io_interpol_spline ! 15.4 Return index of point in grid ! ---------------------------------- FUNCTION ropp_io_seek_index(x, xp) RESULT(ip) IMPLICIT NONE REAL(wp), DIMENSION(:), INTENT(in) :: x ! x-grid (homogeneous) REAL(wp), INTENT(in) :: xp ! point to locate inside grid INTEGER :: ip INTEGER :: N ! Number of grid points INTEGER :: Imin ! Upper estimate of index INTEGER :: Imax ! Lower estimate of index INTEGER :: Dir ! Direction of argument change INTEGER :: It ! Iteration count INTEGER :: di ! Index increment in iterations INTEGER :: is ! Step direction count ! 15.4.1 Grid size and direction calculation N = SIZE(x) Dir = NINT(SIGN(1.0_wp, x(N)-x(1))) ! 15.4.2 Checking if point is inside grid IF ((Dir*xp < Dir*x(1)) .OR. (Dir*xp > Dir*x(N))) THEN ip = 0 RETURN END IF ! 15.4.3 Initial approximation Imin = 1 Imax = N ip = Imin+FLOOR(REAL(Imax - Imin, wp)*(xp - x(Imin))/(x(Imax) - x(Imin))) ip = MAX(1, MIN(ip, N-1)) ! 15.4.4 Iterative index search It = 0 is = 0 Search: DO IF ((Dir*x(ip) <= Dir*xp) .AND. (Dir*xp <= Dir*x(ip+1))) THEN EXIT Search END IF IF (ABS(is) > 1) THEN ip = (Imax + Imin)/2 is = 0 END IF IF (Dir*x(ip+1) < Dir*xp) THEN Imin = ip + 1 di = FLOOR(REAL(Imax - Imin, wp)*(xp - x(Imin))/(x(Imax) - x(Imin))) ip = Imin + di is = is + 1 ELSE IF (Dir*xp < Dir*x(ip)) THEN Imax = ip di = FLOOR(REAL(Imin - Imax, wp)*(xp - x(Imax))/(x(Imin) - x(Imax))) ip = Imax + di is = is - 1 END IF ip = MAX(1, MIN(ip, N-1)) It = It + 1 IF (It > N) THEN ip = -1 EXIT Search END IF END DO Search END FUNCTION ropp_io_seek_index !------------------------------------------------------------------------------- ! 7b. gettime auxiliary routine !------------------------------------------------------------------------------- SUBROUTINE gettime(DTocc) ! Read year, month, day, hour, minute, second from global attributes of ! UCAR atmPhs or atmPrf files, and translate into DTocc structure. USE ropp_io_types, ONLY: DT7type USE ncdf IMPLICIT NONE TYPE(DT7type), INTENT(inout) :: DTocc REAL(wp) :: sec CALL ncdf_getatt('year', DTocc%Year) CALL ncdf_getatt('month', DTocc%Month) CALL ncdf_getatt('day', DTocc%Day) CALL ncdf_getatt('hour', DTocc%Hour) CALL ncdf_getatt('minute', DTocc%Minute) CALL ncdf_getatt('second', sec) DTocc%Second = INT(sec) IF ( sec - DTocc%Second > 0.999_wp ) THEN DTocc%Msec = INT(1000.0_wp*(sec - DTocc%Second)) ELSE DTocc%Msec = NINT(1000.0_wp*(sec - DTocc%Second)) END IF END SUBROUTINE gettime END SUBROUTINE ropp_io_read_ncdf_get_ucardata !------------------------------------------------------------------------------- ! 8. EUM RO data !------------------------------------------------------------------------------- SUBROUTINE ropp_io_read_ncdf_get_eumdata(DATA, file, getlevel1b, getextrap, & resolution, getlevel1a, getdirect, & getiono,ldummy) ! 8.1 Declarations ! ---------------- USE ncdf USE ropp_utils USE ropp_io_types, ONLY: ROprof IMPLICIT NONE TYPE(ROprof), INTENT(INOUT) :: DATA CHARACTER(len = *), INTENT(IN) :: file LOGICAL, INTENT(IN) :: getlevel1b LOGICAL, INTENT(IN) :: getextrap CHARACTER(len = 20), INTENT(IN) :: resolution CHARACTER(len = 20), INTENT(IN) :: getlevel1a LOGICAL, INTENT(IN) :: getdirect LOGICAL, INTENT(IN) :: getiono LOGICAL, INTENT(IN) :: ldummy ! needed to differentiate between ! ropp_io_read_ncdf_get_eumdata and ! ropp_io_read_ncdf_get_otherdata CHARACTER(len = 256) :: routine ! 8.2 Error handling ! ------------------ CALL message_get_routine(routine) CALL message_set_routine('ropp_io_read_eumdata') CALL ropp_io_read_eumdata(DATA, file, getlevel1b, getextrap, resolution, & getlevel1a, getdirect, getiono) write(*, *) 'data%Lev1a%Npoints after ropp_io_read_eumdata', DATA%Lev1a%Npoints CALL message_set_routine(routine) CONTAINS !------------------------------------------------------------------------------- ! 9. EUMETSAT RO data - level 1b and, optionally, level 1a. !------------------------------------------------------------------------------- SUBROUTINE ropp_io_read_eumdata(DATA, file, getlevel1b, getextrap, resolution, & getlevel1a, getdirect, getiono) ! This routine reads the netCDF4 EUMETSAT format into an internal ROPP structure. ! 9.1 Declarations ! ---------------- USE DateTimeTypes USE ropp_utils USE coordinates USE ncdf USE ropp_io USE ropp_io_types, ONLY: ROprof, & ThisFmtVer, & PCD_summary, & PCD_offline, & PCD_phase, & PCD_bangle, & PCD_open_loop, & PCD_rising IMPLICIT NONE TYPE(ROprof), INTENT(INOUT) :: data CHARACTER(len = *), INTENT(IN) :: file LOGICAL, INTENT(IN) :: getlevel1b LOGICAL, INTENT(IN) :: getextrap CHARACTER(len = 20), INTENT(IN) :: resolution CHARACTER(len = 20), INTENT(IN) :: getlevel1a LOGICAL, INTENT(IN) :: getdirect LOGICAL, INTENT(IN) :: getiono TYPE(ROprof) :: data_cl, data_olrs CHARACTER(len = 30) :: sampling_mode CHARACTER(len = 30) :: data_available CHARACTER(len = 30) :: external_navbits_applied CHARACTER(len = 30) :: navbits_internal CHARACTER(len = 30) :: navbits_external CHARACTER(len = 30) :: pmode ! Temporarily needed to read old SE6 data REAL(wp), DIMENSION(:), ALLOCATABLE :: navbit_ext ! OL or RS navbits (external) REAL(wp), DIMENSION(:), ALLOCATABLE :: navbit_int ! OL or RS navbits (internal) REAL(wp), DIMENSION(:), ALLOCATABLE :: temparray ! Temporary array INTEGER, DIMENSION(:), ALLOCATABLE :: olrs_lcf, cl_lcf, lcf ! Lost carrier flag REAL(wp), DIMENSION(:,:), ALLOCATABLE :: r_gns, v_gns ! Temporary pos/velocity arrays REAL(wp), DIMENSION(:,:), ALLOCATABLE :: r_leo, v_leo ! Temporary pos/velocity arrays REAL(wp), DIMENSION(2) :: frequencies ! Temporary frequency array CHARACTER(len = 3), DIMENSION(2) :: codes ! Temporary codes array REAL(wp), DIMENSION(:), ALLOCATABLE :: i_ca_uncorr ! OL or RS I component [V] REAL(wp), DIMENSION(:), ALLOCATABLE :: q_ca_uncorr ! OL or RS Q component [V] REAL(wp), DIMENSION(:), ALLOCATABLE :: exphase_l1_nco ! OL or RS NCO excess phase [m] REAL(wp), DIMENSION(:), ALLOCATABLE :: phase_l1_iq ! I/Q contribution INTEGER, DIMENSION(:), ALLOCATABLE :: tracking_state ! Tracking state ! Variables used for reconstructing OL SNRs until problem fixed at EUMETSAT. REAL(wp), DIMENSION(:), ALLOCATABLE :: dtime ! RS times REAL(wp), DIMENSION(:), ALLOCATABLE :: i_ca ! RS I component [V] REAL(wp), DIMENSION(:), ALLOCATABLE :: q_ca ! RS Q component [V] REAL(wp), DIMENSION(:), ALLOCATABLE :: snr_ca ! RS SNRs [V/V] REAL(wp), DIMENSION(:), ALLOCATABLE :: i_ca_avr ! Averaged I component [V] REAL(wp), DIMENSION(:), ALLOCATABLE :: q_ca_avr ! Averaged Q component [V] INTEGER, DIMENSION(:), ALLOCATABLE :: open_loop_L1, open_loop_L2 ! Yong Chen on 01/06/2026 to handle 1x and 5x different time dimensions INTEGER :: n_t1, n_t5, n_common REAL(wp), PARAMETER :: dt_tol=1.0e-9_wp ! L1 variables REAL(wp),DIMENSION(:), ALLOCATABLE :: dtime_1x(:) REAL(wp),DIMENSION(:), ALLOCATABLE :: exphase_1x(:) REAL(wp),DIMENSION(:), ALLOCATABLE :: snr_1x(:) INTEGER, DIMENSION(:), ALLOCATABLE :: open_loop_1x(:) REAL(wp),DIMENSION(:,:), ALLOCATABLE :: r_transmitter_1x REAL(wp),DIMENSION(:,:), ALLOCATABLE :: v_transmitter_1x REAL(wp),DIMENSION(:,:), ALLOCATABLE :: r_receiver_1x REAL(wp),DIMENSION(:,:), ALLOCATABLE :: v_receiver_1x ! L5 variables REAL(wp),DIMENSION(:), ALLOCATABLE :: dtime_5x(:) REAL(wp),DIMENSION(:), ALLOCATABLE :: exphase_5x(:) REAL(wp),DIMENSION(:), ALLOCATABLE :: snr_5x(:) INTEGER, DIMENSION(:), ALLOCATABLE :: open_loop_5x(:) REAL(wp) :: i_ca_sum ! Temporary sum of I components REAL(wp) :: q_ca_sum ! Temporary sum of Q components REAL(wp) :: dtime_sum ! Temporary sum dtimes REAL(wp) :: angle_iq ! Angle in complex (I,Q) plane REAL(wp) :: open_loop_L1_diff, open_loop_L2_diff INTEGER :: n_rs, j1, j2 ! Number of RS and counters INTEGER :: open_loop_L1_num, open_loop_L2_num CHARACTER(LEN=15) :: cval ! Message string INTEGER :: n, n_cl, n_olrs, first_valid, i, j, k, ierr INTEGER :: readint, readint2, readint3 REAL(EightByteReal) :: readreal, readreal2, ts_cl, ts_olrs, ts1, ts1_max CHARACTER (len = 256) :: readstr, readstr2, fver INTEGER(OneByteInt) :: readbyte1, readbyte2, readbyte3 CHARACTER (len = 256) :: sdir ! directory to science info in netCDF4 CHARACTER (len = 256) :: ddir ! directory to science/data in netCDF4 CHARACTER (len = 256) :: gdir ! directory to science/data in netCDF4 CHARACTER (len = 256) :: tdir ! tmp directory for use ! by Xinjia, 05-05-2025, for METOP-SG CHARACTER (len = 256) :: vname ! actual variable name, stripped of / INTEGER :: groupid, numgrps, status ! group ID (if applicable) LOGICAL :: havegroup ! found group INTEGER, DIMENSION(2) :: ncids REAL(wp), PARAMETER :: f_L1 = 1.57542e9_wp REAL(wp), PARAMETER :: c_light = 299792458.0_wp ! Pi already defined as parameter in coordinates module. Confuses pgf95. REAL(wp), PARAMETER :: pi1 = 3.141592653589793238_wp INTEGER :: have_nb = 0 INTEGER :: DT8(8) ! INTEGER, DIMENSION(:), POINTER :: idx => null() ! Holds 'where' output INTEGER :: nidx ! 9.2 Group path for science, level 1B data ! ----------------------------------------- sdir = '/data/' !ddir = TRIM(sdir) // 'level_1b/' // TRIM(resolution) // '/' ddir = TRIM(sdir) // 'level_1b/high_resolution/' ! <-- changed by Xinjia IF ( getiono ) ddir = TRIM(sdir) // 'ionosphere/' // TRIM(resolution) // '/' ! 9.3 Header variables ! -------------------- ! receiver satellite ID in ROPP format readstr = ' ' CALL ncdf_getatt('spacecraft', readstr) IF ( readstr == 'M02' ) data%leo_id = 'META' IF ( readstr == 'M01' ) data%leo_id = 'METB' IF ( readstr == 'M03' ) data%leo_id = 'METC' IF ( readstr == 'C01' ) data%leo_id = 'C001' IF ( readstr == 'C02' ) data%leo_id = 'C002' IF ( readstr == 'C03' ) data%leo_id = 'C003' IF ( readstr == 'C04' ) data%leo_id = 'C004' IF ( readstr == 'C05' ) data%leo_id = 'C005' IF ( readstr == 'C06' ) data%leo_id = 'C006' IF ( readstr == 'S6A' ) data%leo_id = 'SE6A' IF ( readstr == 'S6B' ) data%leo_id = 'SE6B' IF ( readstr == 'SGA1' ) data%leo_id = 'MESG' ! by Xinjia, 05-05-2025 IF ( readstr(1:2) == 'FM' ) data%leo_id = 'L'//readstr(3:5) ! Unique ID for Lemur/Spire satellites CALL ncdf_getatt('/status/processing/format_version', fver) READ ( fver, fmt=*, iostat=ierr ) ! Clever IF ( ierr /= 0 ) fver = "Unknown" data%FmtVersion = TRIM(fver) CALL message(msg_info, "Format Version " // data%FmtVersion ) CALL ncdf_getatt('/status/processing/processing_mode', pmode) ! in example file 's6x.nc', pmode is 'NTC' ! in ROMEX, pmode is 'NTR' ! in Metop-SG example, pmode is 'NTR' pmode = 'NTC' ! change it to be compliance to s6x.nc ! occulting GNSS satellite readstr = ' ' CALL ncdf_getatt('/data/transmitter/satellite/satellite_prn', readstr) READ ( readstr(2:), FMT='(I3.3)' ) readint WRITE ( data%gns_id, FMT='(A1,I3.3)' ) readstr(1:1), readint ! station ID unused for EUMETSAT data data%stn_id = ' ' ! 9.4 Product Confidence Data ! --------------------------- data%PCD = 0 data%overall_qual = 100.0_wp data%units%overall_qual = '%' IF (getlevel1b) THEN ! Overall nominal / non-nominal tdir = '/quality/' CALL ncdf_getvar(TRIM(tdir)//'overall_quality_ok', readbyte1) IF (readbyte1 == 0) data%PCD = IBSET(data%PCD, PCD_summary) ! BA nominal (set to same value as overall) IF (readbyte1 == 0) data%PCD = IBSET(data%PCD, PCD_bangle) ! set to same value as overall ! quality indicator, currently just set to 100% or 0% data%overall_qual = 100.0_wp * readbyte1 ENDIF ! NRT or offline processing? CALL ncdf_getatt('environment', readstr) IF ( readstr == 'Operational' .OR. readstr == 'Validation' ) THEN data%PCD = IBCLR(data%PCD, PCD_offline) ELSE IF (readstr == 'Offline') data%PCD = IBSET(data%PCD, PCD_offline) END IF ! setting or rising IF (data%FmtVersion < '13.0') THEN IF ( getiono ) THEN CALL ncdf_getatt('/data/ionosphere/occultation/occultation_type', readstr) ELSE CALL ncdf_getatt('/data/occultation/occultation_type', readstr) ENDIF ELSE ! FmtVersion is 14.0 IF ( getiono ) THEN CALL ncdf_getvar('/data/ionosphere/occultation/occultation_type', readstr) ELSE CALL ncdf_getvar('/data/occultation/occultation_type', readstr) ENDIF ENDIF IF (TRIM(readstr) == 'rising') data%PCD = IBSET(data%PCD, PCD_rising) ! 9.5 Date and time ! ----------------- ! start time of occultation tdir = TRIM(sdir) // 'level_1a/' IF ( getiono ) tdir = TRIM(sdir) // 'ionosphere/' CALL ncdf_getvar(TRIM(tdir)//'utc_start_absdate', readint) CALL ncdf_getvar(TRIM(tdir)//'utc_start_abstime', readreal) ! check that we have the right time CALL ncdf_getatt(TRIM(tdir)//'units', readstr, 'utc_start_absdate') ! IF ( readstr == 'days since 2000-01-01 00:00:00.00' ) THEN ! changed by Xinjia, 05-05-2025 IF ( readstr == 'days since 2000-01-01 00:00:00.00' .or. readstr == 'days since 2000-01-01' ) THEN CALL abstimetoDT(readint, readreal, data%DTocc) ELSE CALL message(msg_fatal, & 'Time units found for utc_start_absdate are inconsistent.') ENDIF ! difference between start time and georef time tdir = TRIM(sdir) // 'occultation/' IF ( getiono ) tdir = TRIM(sdir) // 'ionosphere/occultation/' CALL ncdf_getvar(TRIM(tdir)//'utc_georef_absdate', readint2) CALL ncdf_getvar(TRIM(tdir)//'utc_georef_abstime', readreal2) ! check that we have the right time CALL ncdf_getatt(TRIM(tdir)//'units', readstr, 'utc_georef_absdate') IF ( .NOT. (readstr == 'days since 2000-01-01 00:00:00.00') ) & CALL message(msg_fatal, & 'Time units found for utc_georef_absdate are inconsistent.') data%georef%time_offset = 86400.d0*(readint2 - readint) + & ! This is OK as long as there are no leapseconds between the two times (readreal2 - readreal) ! 9.6 Georeferencing ! ------------------ CALL ncdf_getvar(TRIM(tdir)//'latitude', data%georef%lat, & units=data%georef%units%lat) CALL ncdf_getvar(TRIM(tdir)//'longitude', data%georef%lon, & units=data%georef%units%lon) CALL ncdf_getvar(TRIM(tdir)//'undulation', data%georef%Undulation, & units=data%georef%units%Undulation) CALL ncdf_getvar(TRIM(tdir)//'r_curve', data%georef%roc, & units=data%georef%units%roc) CALL ncdf_getvar(TRIM(tdir)//'r_curve_centre_fixed', data%georef%r_coc, & units=data%georef%units%r_coc) CALL ncdf_getvar(TRIM(tdir)//'azimuth_north', data%georef%azimuth, & units=data%georef%units%azimuth) CALL ncdf_getvar(TRIM(tdir)//'position_rec_fixed', data%georef%leo_pod%pos, & units=data%georef%units%r_leo, range = data%georef%range%r_leo) CALL ncdf_getatt(TRIM(tdir)//'long_name', readstr, 'position_rec_fixed') IF ( readstr /= 'Receiver antenna position in Earth fixed coordinates (for SLTA = 0)' .AND. & readstr /= 'Receiver antenna position in Earth fixed coordinates for reference location' ) THEN CALL message(msg_warn, 'Unrecognised long_name ("' // TRIM(readstr) // & '") for position_rec_fixed ... ' // & 'will assume them to be in ECF coordinates') END IF data%georef%leo_pod%pos_frame = 'ECF' IF ( (data%FmtVersion < '13.0') .OR. & (data%leo_id(1:3) == 'SE6' .AND. pmode == 'NRT') ) THEN ! Temporary nonsense for Sentinel-6 with no *_bfr variables CALL ncdf_getvar(TRIM(tdir)//'velocity_rec', data%georef%leo_pod%vel, & units=data%georef%units%v_leo, range = data%georef%range%v_leo) CALL ncdf_getatt(TRIM(tdir)//'long_name', readstr, 'velocity_rec') IF ( readstr /= 'Receiver antenna velocity in Earth centred inertial coordinates (J2000, for SLTA = 0)' .AND. & readstr /= 'Receiver antenna velocity in Earth centred inertial coordinates for reference location' ) THEN CALL message(msg_warn, 'Unrecognised long_name ("' // TRIM(readstr) // & '") for velocity_rec ... ' // & 'will assume them to be in ECI coordinates') END IF ELSE CALL ncdf_getvar(TRIM(tdir)//'velocity_rec_fixed_bfr', data%georef%leo_pod%vel, & units=data%georef%units%v_leo, range = data%georef%range%v_leo) CALL ncdf_getatt(TRIM(tdir)//'long_name', readstr, 'velocity_rec_fixed_bfr') IF ( readstr /= 'Receiver antenna inertial velocity in Earth fixed coordinates for reference location (for WMO BUFR' .AND. & readstr /= 'Receiver antenna inertial velocity in Earth fixed coordinates for reference location (for WMO BUFR)' ) THEN CALL message(msg_warn, 'Unrecognised long_name ("' // TRIM(readstr) // & '") for velocity_rec_fixed_bfr ... ' // & 'will assume them to be in ECI coordinates') END IF END IF data%georef%leo_pod%vel_frame = 'ECI' CALL ncdf_getvar(TRIM(tdir)//'position_gns_fixed', data%georef%gns_pod%pos, & units=data%georef%units%r_gns, range = data%georef%range%r_gns) CALL ncdf_getatt(TRIM(tdir)//'long_name', readstr, 'position_gns_fixed') IF ( readstr /= 'GNSS transmitter position in Earth fixed coordinates (for SLTA = 0)' .AND. & readstr /= 'GNSS transmitter position in Earth fixed coordinates for reference location' ) THEN CALL message(msg_warn, 'Unrecognised long_name ("' // TRIM(readstr) // & '") for position_gns_fixed ... ' // & 'will assume them to be in ECF coordinates') END IF data%georef%gns_pod%pos_frame = 'ECF' IF ( (data%FmtVersion < '13.0') .OR. & (data%leo_id(1:3) == 'SE6' .AND. pmode == 'NRT') ) THEN ! Temporary nonsense for Sentinel-6 with no *_bfr variables CALL ncdf_getvar(TRIM(tdir)//'velocity_gns', data%georef%gns_pod%vel, & units=data%georef%units%v_gns, range = data%georef%range%v_gns) CALL ncdf_getatt(TRIM(tdir)//'long_name', readstr, 'velocity_gns') IF ( readstr /= 'GNSS transmitter velocity in Earth centred inertial coordinates (J2000, for SLTA = 0)' .AND. & readstr /= 'GNSS transmitter velocity in Earth centred inertial coordinates for reference location' ) THEN CALL message(msg_warn, 'Unrecognised long_name ("' // TRIM(readstr) // & '") for velocity_gns ... ' // & 'will assume them to be in ECI coordinates') END IF ELSE CALL ncdf_getvar(TRIM(tdir)//'velocity_gns_fixed_bfr', data%georef%gns_pod%vel, & units=data%georef%units%v_gns, range = data%georef%range%v_gns) CALL ncdf_getatt(TRIM(tdir)//'long_name', readstr, 'velocity_gns_fixed_bfr') IF ( readstr /= 'GNSS transmitter inertial velocity in Earth fixed coordinates for reference location (for WMO BUFR)' ) THEN CALL message(msg_warn, 'Unrecognised long_name ("' // TRIM(readstr) // & '") for velocity_gns_fixed_bfr ... ' // & 'will assume them to be in ECI coordinates') END IF END IF data%georef%gns_pod%vel_frame = 'ECI' ! FIXME: Do we need this? data%georef%reference_frame%r_coc = 'ECF' ! 9.7 Signal characteristics ! -------------------------- !print *, 'before 9.7 Signal characteristics' ! In principle we should read the codes, but this doesn't work because ROPP can ! only read scalar strings from netCDF files, and the codes are a string vector. ! Instead we (temporarily) find the (only two) L1 and L2 codes in a different way. ! The order of checking for the L1 and L2 codes below matters. The least desirable ! codes comes first, but are then overwritten if something more desirable exists. ! All this (and the associated pieces of code below) is not very good, and needs ! to be revised once ROPP has been modified to read a string vector, and when the ! (at the moment ongoing) changes to the format from EUMETSAT's side settles. ! for Metop (not EPS-SG) we know what the codes and frequencies are. In general, the codes are ! used for reading various variables in newer format versions below. For older versions (< 13.0) ! the reading with codes (ca, p1, and p2) appended to variables are hardcoded below as it has always been. IF (data%leo_id(1:3) == 'MET') THEN codes = (/"1c", "2w"/) frequencies = (/1.57542E9_wp, 1.2276E9_wp/) ELSE IF (data%leo_id(1:3) == 'SE6' .AND. pmode == 'NRT') THEN ! older validation version where processing_mode was erroneously set to NRT ! print *, 'before 9.7 /data/level_1a/open_loop/snr_1w' IF ( ncdf_isvar('/data/level_1a/open_loop/snr_1w') ) codes(1) = "1w" ! print *, 'before 9.7 /data/level_1a/open_loop/snr_1p' IF ( ncdf_isvar('/data/level_1a/open_loop/snr_1p') ) codes(1) = "1p" ! print *, 'before 9.7 /data/level_1a/open_loop/snr_1c' IF ( ncdf_isvar('/data/level_1a/open_loop/snr_1c') ) codes(1) = "1c" ! print *, 'before 9.7 /data/level_1a/open_loop/snr_2w' IF ( ncdf_isvar('/data/level_1a/open_loop/snr_2w') ) codes(2) = "2w" ! print *, 'before 9.7 /data/level_1a/open_loop/snr_2p' IF ( ncdf_isvar('/data/level_1a/open_loop/snr_2p') ) codes(2) = "2p" ! print *, 'before 9.7 /data/level_1a/open_loop/snr_2l' IF ( ncdf_isvar('/data/level_1a/open_loop/snr_2l') ) codes(2) = "2l" ! print *, 'before 9.7 /data/level_1a/open_loop/snr_2c' IF ( ncdf_isvar('/data/level_1a/open_loop/snr_2c') ) codes(2) = "2c" ! print *, 'before 9.7 /data/level_1a/open_loop/frequencies' CALL ncdf_getvar('/data/level_1a/open_loop/frequencies', frequencies) ELSE IF ( data%leo_id == 'MESG' ) THEN !print *, 'before 9.7 /data/level_1a/combined for MESG' ! by Xinjia, 05/05/2025. ! ropp_io/ncdf/ncdf_getgroupid_n4.f90 !print *, 'line 3261, test ELSE IF ( data%leo_id == MESG )' ! < --- reach here status = ncdf_getgroupid(ncdf_ncid, '/data/level_1a/combined/signals', & vname, groupid, havegroup) ! a trick to get the E1/E5a or L1/L5 !print *, 'after ncdf_getgroupid,',trim(vname), status status = nf90_inq_grps(groupid, numgrps, ncids) status = NF90_INQ_GRPNAME(ncids(1), codes(1)) status = NF90_INQ_GRPNAME(ncids(2), codes(2)) !print *, 'numgrps=',numgrps !status = ncdf_getgroupid(ncdf_ncid, '/data/level_1a/combined/L1', vname, groupid, havegroup) !print *, 'after ncdf_getgroupid,',vname, havegroup, status !status = ncdf_getgroupid(ncdf_ncid, '/data/level_1a/combined/B1C', vname, groupid, havegroup) !print *, 'after ncdf_getgroupid,',vname, havegroup, status !print *, codes(1) ! E1/L1 !print *, codes(2) ! E5a/L5 !print *, 'before 9.7 /data/level_1a/combined/'//trim(codes(1))//'/frequency' CALL ncdf_getvar('/data/level_1a/combined/'//trim(codes(1))//'/frequency', frequencies(1)) !print *, 'before 9.7 /data/level_1a/combined/'//trim(codes(2))//'/frequency' CALL ncdf_getvar('/data/level_1a/combined/'//trim(codes(2))//'/frequency', frequencies(2)) ELSE ! In newer format version the 'open_loop' group has been renamed to 'combined' and ! codes are appended to the frequency names. !! The out-commented line below doesn't work because ROPP can only read scalar strings from netCDF files ! CALL ncdf_getvar('/data/level_1a/combined/codes', codes) !print *, 'line 3291, test ELSE' ! s6a, confirmed. Xinjia 08192025 ! ropp_io/ncdf/ncdf_isvar.f90 IF ( ncdf_isvar('/data/level_1a/combined/snr_1w') ) codes(1) = "1w" IF ( ncdf_isvar('/data/level_1a/combined/snr_1p') ) codes(1) = "1p" IF ( ncdf_isvar('/data/level_1a/combined/snr_1c') ) codes(1) = "1c" IF ( ncdf_isvar('/data/level_1a/combined/snr_2w') ) codes(2) = "2w" IF ( ncdf_isvar('/data/level_1a/combined/snr_2p') ) codes(2) = "2p" IF ( ncdf_isvar('/data/level_1a/combined/snr_2x') ) codes(2) = "2x" IF ( ncdf_isvar('/data/level_1a/combined/snr_2l') ) codes(2) = "2l" IF ( ncdf_isvar('/data/level_1a/combined/snr_2c') ) codes(2) = "2c" !print *, codes(1) !print *, codes(2) CALL ncdf_getvar('/data/level_1a/combined/frequency_'//codes(1), frequencies(1)) CALL ncdf_getvar('/data/level_1a/combined/frequency_'//codes(2), frequencies(2)) END IF END IF data%occfreq1 = frequencies(1) data%occfreq2 = frequencies(2) !data%signal1%freq = frequencies(1) !data%signal2%freq = frequencies(2) ! 9.7 Level1a variables (if any) ! ------------------------------ ! 9.7.1 Closed Loop Level1a variables (if requested) ! -------------------------------------------------- IF ( TRIM(ADJUSTL(getlevel1a)) /= 'none' ) THEN n_cl = 0 ! will be set differently if CL data is available (or if it is Sentinel-6 or Spire) n_olrs = 0 ! will be set differently if OL or RS data for Metop is available (whichever is requested) SELECT CASE (data%leo_id) CASE ("META", "METB", "METC") IF ( ncdf_isvar('/quality/cl_data_available') ) THEN CALL ncdf_getvar('/quality/cl_data_available', readbyte1) ts_cl = 0.02_wp ! expected (approximate) time in seconds between closed loop samples ELSE readbyte1 = 0 ENDIF CASE ("SE6A", "SE6B") ! by Xinjia readbyte1 = 1 ! Sentinel-6 has only open loop data, but we read it as closed loop because the ! handling of Sentinel-6 data is similar to Metop CL data (e.g., there is both L1 ! and L2, and we just read them directly, we don't care to reconstruct from I's and Q's). data%PCD = IBSET(data%PCD, PCD_open_loop) ts_cl = 0.01_wp ! expected (approximate) time in seconds between samples for Sentinel-6. CASE ("MESG") ! by Xinjia readbyte1 = 1 ! Sentinel-6 has only open loop data, but we read it as closed loop because the ! handling of Sentinel-6 data is similar to Metop CL data (e.g., there is both L1 ! and L2, and we just read them directly, we don't care to reconstruct from I's and Q's). data%PCD = IBSET(data%PCD, PCD_open_loop) ! Yong Chen 12/23/2025 IF( data%gns_id(1:1) =="E") THEN ts_cl = 0.004_wp ! expected (approximate) time in seconds between samples for Galileo (250Hz) ELSE IF (data%gns_id(1:1) =="G" .OR. data%gns_id(1:1) =="C") THEN ts_cl = 0.005_wp ! expected (approximate) time in seconds between samples for GPS and Beidou (200Hz) ELSE ts_cl = 0.005_wp ! default to 200Hz ENDIF CASE DEFAULT IF (data%leo_id(1:1) == "L") THEN ! This is a Lemur/Spire ID. If there are other missions ! with an ID starting with 'L', they should be listed explicitly above. readbyte1 = 1 ! Spire has only combined data, but we read it as closed loop because the ! handling of Spire data is similar to Metop CL data (e.g., there is both L1 ! and L2, and we just read them directly, we don't care to reconstruct from I's and Q's). data%PCD = IBSET(data%PCD, PCD_open_loop) ! assume open loop data. ts_cl = 0.02_wp ! expected (approximate) time in seconds between samples for Spire. ELSE readbyte1 = 0 ENDIF END SELECT IF (readbyte1 == 1) THEN ! CL data exist (or it is Sentinel-6 or Spire) !print *, 'tdir (before 3330): ', tdir IF (data%leo_id(1:3) == 'MET') THEN tdir = TRIM(sdir) // 'level_1a/closed_loop/' ELSE IF (data%leo_id(1:3) == 'SE6' .AND. pmode == 'NRT') THEN ! older validation version of Sentinel-6 data tdir = TRIM(sdir) // 'level_1a/open_loop/' ELSE ! newer format version print *, 'newer format version, level_1a/combined/' ! < --- reach here, msg and s6a tdir = TRIM(sdir) // 'level_1a/combined/' ENDIF ENDIF !print *, 'tdir (after 3330): ', tdir ! by Xinjia, 05/07/2025 ! Yong Chen, 01/06/2026 for common dataset in level_1a/combined/ for 1x and 5x IF ( data%leo_id == 'MESG' ) THEN CALL ncdf_getsize(TRIM(tdir)//trim(codes(1))//'/dtime', n_t1, dim = 1) CALL ncdf_getsize(TRIM(tdir)//trim(codes(2))//'/dtime', n_t5, dim = 1) ALLOCATE(dtime_1x(n_t1), exphase_1x(n_t1), snr_1x(n_t1), & open_loop_1x(n_t1), r_transmitter_1x(3,n_t1), & v_transmitter_1x(3,n_t1), r_receiver_1x(3,n_t1), & v_receiver_1x(3,n_t1)) ALLOCATE(dtime_5x(n_t5), exphase_5x(n_t5), snr_5x(n_t5), & open_loop_5x(n_t5)) ! codes(1)//dtime is the difference with codes(2)//dtime CALL ncdf_getvar(TRIM(tdir)//trim(codes(1))//'/dtime', dtime_1x, & units=data_cl%Lev1a%units%dtime) CALL ncdf_getvar(TRIM(tdir)//trim(codes(2))//'/dtime', dtime_5x, & units=data_cl%Lev1a%units%dtime) ! count common samples (two-pointers) i = 1 j = 1 n_common = 0 do while (i <= n_t1 .and. j <= n_t5) if (abs(dtime_1x(i) - dtime_5x(j)) <= dt_tol) then n_common = n_common + 1 i = i + 1 j = j + 1 else if (dtime_1x(i) < dtime_5x(j)) then i = i + 1 else j = j + 1 end if end do n_cl = n_common print *, 'n_t1, n_t5 and n_common:', n_t1, n_t5, n_common print *, 'dtime:', dtime_1x(1), dtime_1x(n_t1), dtime_5x(1), dtime_5x(n_t5) ELSE CALL ncdf_getsize(TRIM(tdir)//'dtime', n_cl, dim = 1) ENDIF print *, 'n_cl=',n_cl ! <--- the same to the dtime size, msg and s6a IF ( n_cl > 1 ) THEN ! read CL data - needs at least 2 samples to make time differences CALL ropp_io_init(data_cl%Lev1a, n_cl) ! Are closed loop phase measurements okay? IF (data%FmtVersion < '13.0') THEN CALL ncdf_getvar('/quality/cl_snr_ca_ok', readbyte1) CALL ncdf_getvar('/quality/cl_snr_p2_ok', readbyte2) ELSE IF ( data%leo_id == 'MESG' ) THEN !print *, ' ( data%leo_id == MESG ), /quality/snr_l1_ok' ! < -- reach here CALL ncdf_getvar('/quality/snr_l1_ok', readbyte1) CALL ncdf_getvar('/quality/snr_l5_ok', readbyte2) ELSE !print *, 'before snr_l1_ok' !print *, 'line 3392, ELSE, /quality/snr_l1_ok' ! < --- reach here, s6a CALL ncdf_getvar('/quality/snr_l1_ok', readbyte1) CALL ncdf_getvar('/quality/snr_l2_ok', readbyte2) ENDIF IF ( (readbyte1 == 0) .AND. (readbyte2 == 0) ) THEN data%PCD = IBSET(data%PCD, PCD_phase) data%PCD = IBSET(data%PCD, PCD_summary) data%overall_qual = 0.0_wp ENDIF ! Time IF ( data%leo_id == 'MESG' ) THEN CALL ncdf_getvar(TRIM(tdir)//trim(codes(1))//'/r_transmitter_1x', r_transmitter_1x, & units=data_cl%Lev1a%units%r_gns ) CALL ncdf_getvar(TRIM(tdir)//trim(codes(1))//'/v_transmitter_1x', v_transmitter_1x, & units=data_cl%Lev1a%units%v_gns ) CALL ncdf_getvar(TRIM(tdir)//trim(codes(1))//'/r_receiver_1x', r_receiver_1x, & units=data_cl%Lev1a%units%r_leo ) CALL ncdf_getvar(TRIM(tdir)//trim(codes(1))//'/v_receiver_1x', v_receiver_1x, & units=data_cl%Lev1a%units%v_leo ) CALL ncdf_getvar(TRIM(tdir)//trim(codes(1))//'/exphase_1x', exphase_1x, & units=data_cl%Lev1a%units%phase) CALL ncdf_getvar(TRIM(tdir)//trim(codes(1))//'/snr_1x', snr_1x, & units=data_cl%Lev1a%units%snr) CALL ncdf_getvar(TRIM(tdir)//trim(codes(1))//'/open_loop', open_loop_1x) CALL ncdf_getvar(TRIM(tdir)//trim(codes(2))//'/exphase_5x', exphase_5x, & units=data_cl%Lev1a%units%phase) CALL ncdf_getvar(TRIM(tdir)//trim(codes(2))//'/snr_5x', snr_5x, & units=data_cl%Lev1a%units%snr) CALL ncdf_getvar(TRIM(tdir)//trim(codes(2))//'/open_loop', open_loop_5x) ALLOCATE(open_loop_L1(n_cl), open_loop_L2(n_cl)) ALLOCATE(r_gns(3,n_cl), v_gns(3,n_cl), r_leo(3,n_cl), v_leo(3,n_cl)) i = 1 j = 1 n_common = 0 do while (i <= n_t1 .and. j <= n_t5) if (abs(dtime_1x(i) - dtime_5x(j)) <= dt_tol) then n_common = n_common + 1 data_cl%Lev1a%dtime(n_common) = dtime_1x(i) data_cl%Lev1a%phase_L1(n_common) = exphase_1x(i) data_cl%Lev1a%snr_L1ca(n_common) = snr_1x(i) data_cl%Lev1a%snr_L1p(n_common) = snr_1x(i) data_cl%Lev1a%phase_L2(n_common) = exphase_5x(j) data_cl%Lev1a%snr_L2p(n_common) = snr_5x(j) open_loop_L1(n_common) = open_loop_1x(i) open_loop_L2(n_common) = open_loop_5x(j) r_gns(:,n_common) = r_transmitter_1x(:, i) v_gns(:,n_common) = v_transmitter_1x(:, i) r_leo(:,n_common) = r_receiver_1x(:, i) v_leo(:,n_common) = v_receiver_1x(:, i) i = i + 1 j = j + 1 else if (dtime_1x(i) < dtime_5x(j)) then i = i + 1 else j = j + 1 end if end do data_cl%Lev1a%r_gns = TRANSPOSE(r_gns) data_cl%Lev1a%v_gns = TRANSPOSE(v_gns) data_cl%Lev1a%r_leo = TRANSPOSE(r_leo) data_cl%Lev1a%v_leo = TRANSPOSE(v_leo) DEALLOCATE(r_gns, v_gns, r_leo, v_leo) ! --- combine excess phase using the open_loop flag --- open_loop_L1_num = sum(open_loop_L1) open_loop_L2_num = sum(open_loop_L2) print *, 'open_loop number:', open_loop_L1_num, open_loop_L2_num if (open_loop_L1(1) == 1) then ! rising, open_loop first, close_loop after print *, 'rising RO' print *, 'L1 phsse:', data_cl%Lev1a%phase_L1(1), data_cl%Lev1a%phase_L1(n_cl) print *, 'L2 phsse:', data_cl%Lev1a%phase_L2(1), data_cl%Lev1a%phase_L2(n_cl) ! --- phase_L1 and phase_L2 are reversed in raw data. Switch L2. By Xinjia data_cl%Lev1a%phase_L1(:) = data_cl%Lev1a%phase_L1(:) - data_cl%Lev1a%phase_L1(n_cl) data_cl%Lev1a%phase_L2(:) = data_cl%Lev1a%phase_L2(n_cl) - data_cl%Lev1a%phase_L2(:) !print *, 'data_cl%Lev1a%phase_L1 (after):',data_cl%Lev1a%phase_L1(1:3) ! --- align phase_L1 and phase_L2. By Xinjia open_loop_L1_diff = data_cl%Lev1a%phase_L1(open_loop_L1_num) - & data_cl%Lev1a%phase_L1(open_loop_L1_num+1) open_loop_L2_diff = data_cl%Lev1a%phase_L2(open_loop_L2_num) - & data_cl%Lev1a%phase_L2(open_loop_L2_num+1) data_cl%Lev1a%phase_L1(1:open_loop_L1_num) = & data_cl%Lev1a%phase_L1(1:open_loop_L1_num) - open_loop_L1_diff data_cl%Lev1a%phase_L2(1:open_loop_L2_num) = & data_cl%Lev1a%phase_L2(1:open_loop_L2_num) - open_loop_L2_diff print *, 'open_loop diff: ', open_loop_L1_diff, open_loop_L2_diff print *, 'After alignment L1 phsse:', data_cl%Lev1a%phase_L1(1), data_cl%Lev1a%phase_L1(n_cl) print *, 'L2 phsse:', data_cl%Lev1a%phase_L2(1), data_cl%Lev1a%phase_L2(n_cl) else ! setting, close_loop first, open_loop after print *, 'setting RO' print *, 'L1 phsse:', data_cl%Lev1a%phase_L1(1), data_cl%Lev1a%phase_L1(n_cl) print *, 'L2 phsse:', data_cl%Lev1a%phase_L2(1), data_cl%Lev1a%phase_L2(n_cl) ! --- phase_L1 and phase_L2 are reversed in raw data. Switch L2. By Xinjia data_cl%Lev1a%phase_L1(:) = data_cl%Lev1a%phase_L1(:) - data_cl%Lev1a%phase_L1(1) data_cl%Lev1a%phase_L2(:) = data_cl%Lev1a%phase_L2(1) - data_cl%Lev1a%phase_L2(:) !data_cl%Lev1a%phase_L2(:) = data_cl%Lev1a%phase_L2(:) - data_cl%Lev1a%phase_L2(1) !print *, 'data_cl%Lev1a%phase_L1 (after):',data_cl%Lev1a%phase_L1(1:3) ! --- align phase_L1 and phase_L2. By Xinjia open_loop_L1_diff = data_cl%Lev1a%phase_L1(n_cl - open_loop_L1_num +1) - & data_cl%Lev1a%phase_L1(n_cl - open_loop_L1_num) open_loop_L2_diff = data_cl%Lev1a%phase_L2(n_cl - open_loop_L2_num +1) - & data_cl%Lev1a%phase_L2(n_cl - open_loop_L2_num) data_cl%Lev1a%phase_L1(n_cl - open_loop_L1_num+1 : n_cl) = & data_cl%Lev1a%phase_L1(n_cl - open_loop_L1_num+1 : n_cl) - open_loop_L1_diff data_cl%Lev1a%phase_L2(n_cl - open_loop_L2_num+1 : n_cl) = & data_cl%Lev1a%phase_L2(n_cl - open_loop_L2_num+1 : n_cl) - open_loop_L2_diff print *, 'open_loop diff: ', open_loop_L1_diff, open_loop_L2_diff print *, 'After align L1 phsse:', data_cl%Lev1a%phase_L1(1), data_cl%Lev1a%phase_L1(n_cl) print *, 'L2 phsse:', data_cl%Lev1a%phase_L2(1), data_cl%Lev1a%phase_L2(n_cl) endif DEALLOCATE(open_loop_L1, open_loop_L2) DEALLOCATE(dtime_1x, exphase_1x, snr_1x, open_loop_1x, r_transmitter_1x, & v_transmitter_1x, r_receiver_1x, v_receiver_1x) DEALLOCATE(dtime_5x, exphase_5x, snr_5x, open_loop_5x) else CALL ncdf_getvar(TRIM(tdir)//'dtime', & data_cl%Lev1a%dtime, units=data_cl%Lev1a%units%dtime) ALLOCATE(r_gns(3,n_cl), v_gns(3,n_cl), r_leo(3,n_cl), v_leo(3,n_cl)) ! tdir is 'level_1a/combined/' CALL ncdf_getvar(TRIM(tdir)//'r_transmitter', r_gns, & units=data_cl%Lev1a%units%r_gns) CALL ncdf_getvar(TRIM(tdir)//'v_transmitter', v_gns, & units=data_cl%Lev1a%units%v_gns) CALL ncdf_getvar(TRIM(tdir)//'r_receiver', r_leo, & units=data_cl%Lev1a%units%r_leo) CALL ncdf_getvar(TRIM(tdir)//'v_receiver', v_leo, & units=data_cl%Lev1a%units%v_leo) data_cl%Lev1a%r_gns = TRANSPOSE(r_gns) data_cl%Lev1a%v_gns = TRANSPOSE(v_gns) data_cl%Lev1a%r_leo = TRANSPOSE(r_leo) data_cl%Lev1a%v_leo = TRANSPOSE(v_leo) DEALLOCATE(r_gns, v_gns, r_leo, v_leo) endif ! Coordinate transformation to True-of-Date system, not (yet) using SOFA library data_cl%Lev1a%r_gns = eci2eci(data%DTocc%year, data%DTocc%month, & data%DTocc%day, data%DTocc%hour, & data%DTocc%minute,data%DTocc%second, & data_cl%Lev1a%dtime + & data%DTocc%msec/1000.0_wp, & data_cl%Lev1a%r_gns) data_cl%Lev1a%v_gns = eci2eci(data%DTocc%year, data%DTocc%month, & data%DTocc%day, data%DTocc%hour, & data%DTocc%minute,data%DTocc%second, & data_cl%Lev1a%dtime + & data%DTocc%msec/1000.0_wp, & data_cl%Lev1a%v_gns) data_cl%Lev1a%r_leo = eci2eci(data%DTocc%year, data%DTocc%month, & data%DTocc%day, data%DTocc%hour, & data%DTocc%minute,data%DTocc%second, & data_cl%Lev1a%dtime + & data%DTocc%msec/1000.0_wp, & data_cl%Lev1a%r_leo) data_cl%Lev1a%v_leo = eci2eci(data%DTocc%year, data%DTocc%month, & data%DTocc%day, data%DTocc%hour, & data%DTocc%minute,data%DTocc%second, & data_cl%Lev1a%dtime + & data%DTocc%msec/1000.0_wp, & data_cl%Lev1a%v_leo) ! Phase and amplitude IF (data%leo_id /= 'MESG' ) THEN IF (data%FmtVersion < '13.0') THEN CALL ncdf_getvar(TRIM(tdir)//'exphase_ca', data_cl%Lev1a%phase_L1, & units=data_cl%Lev1a%units%phase) CALL ncdf_getvar(TRIM(tdir)//'exphase_p2', data_cl%Lev1a%phase_L2, & units=data_cl%Lev1a%units%phase) CALL ncdf_getvar(TRIM(tdir)//'snr_ca', data_cl%Lev1a%snr_L1ca, & units=data_cl%Lev1a%units%snr) CALL ncdf_getvar(TRIM(tdir)//'snr_p1', data_cl%Lev1a%snr_L1p, & units=data_cl%Lev1a%units%snr) CALL ncdf_getvar(TRIM(tdir)//'snr_p2', data_cl%Lev1a%snr_L2p, & units=data_cl%Lev1a%units%snr) ELSE CALL ncdf_getvar(TRIM(tdir)//'exphase_'//trim(codes(1)), data_cl%Lev1a%phase_L1, & units=data_cl%Lev1a%units%phase) CALL ncdf_getvar(TRIM(tdir)//'exphase_'//trim(codes(2)), data_cl%Lev1a%phase_L2, & units=data_cl%Lev1a%units%phase) CALL ncdf_getvar(TRIM(tdir)//'snr_'//trim(codes(1)), data_cl%Lev1a%snr_L1ca, & units=data_cl%Lev1a%units%snr) CALL ncdf_getvar(TRIM(tdir)//'snr_'//trim(codes(1)), data_cl%Lev1a%snr_L1p, & units=data_cl%Lev1a%units%snr) CALL ncdf_getvar(TRIM(tdir)//'snr_'//trim(codes(2)), data_cl%Lev1a%snr_L2p, & units=data_cl%Lev1a%units%snr) ENDIF ENDIF ! Closed loop lost carrier flag (looking for data gaps) ! by Xinjia, 08/04/2025 ALLOCATE(cl_LCF(n_cl)) cl_LCF(:) = 0 ts1_max = 1.05*ts_cl IF (BTEST(data%PCD, PCD_rising)) THEN ! Rising occultation DO j=n_cl-1,1,-1 ts1 = ABS(data_cl%Lev1a%dtime(j) - data_cl%Lev1a%dtime(j+1)) !print *, 'ts1, rising', ts1, ts1_max IF (ts1 > ts1_max) THEN cl_LCF(j) = IBSET(cl_LCF(j), 3) ! 3 is data gap cl_LCF(j+1) = IBSET(cl_LCF(j+1),3) ! ENDIF ENDDO ELSE ! Setting occultation DO j=2,n_cl ts1 = ABS(data_cl%Lev1a%dtime(j) - data_cl%Lev1a%dtime(j-1)) !print *, 'ts1, setting', ts1, ts1_max IF (ts1 > ts1_max) THEN cl_LCF(j) = IBSET(cl_LCF(j), 3) cl_LCF(j-1) = IBSET(cl_LCF(j-1),3) ENDIF ENDDO ENDIF ! Rising occultation ELSE n_cl = 0 END IF ! n_cl > 1 ELSE CALL message(msg_info, 'Closed loop data requested ' // & 'but are not available in file') ENDIF ! /quality/cl_data_available == 1 ! 9.7.2 Raw Sampling or Open Loop Level1a variables (if requested) ! ---------------------------------------------------------------- IF ( (data%leo_id(1:3) == 'MET') .AND. & (TRIM(ADJUSTL(getlevel1a)) == 'cl+ol') .OR. & (TRIM(ADJUSTL(getlevel1a)) == 'cl+rs') ) THEN IF ( TRIM(ADJUSTL(getlevel1a)) == 'cl+ol' ) THEN sampling_mode = 'open_loop' IF (data%FmtVersion < '13.0') THEN data_available = 'ol_data_available' external_navbits_applied = 'ol_external_navbits_applied' navbits_internal = 'navbits_internal' navbits_external = 'navbits_external' ELSE data_available = 'ol_data_available' external_navbits_applied = 'external_navbits_ok' navbits_internal = 'navbits_1c' navbits_external = 'navbits_1c' ENDIF ts_olrs = 0.02_wp ! expected (approximate) time in seconds between open loop samples ENDIF IF ( TRIM(ADJUSTL(getlevel1a)) == 'cl+rs' ) THEN sampling_mode = 'raw_sampling' IF (data%FmtVersion < '13.0') THEN data_available = 'rs_data_available' external_navbits_applied = 'rs_external_navbits_applied' navbits_internal = 'navbits_internal' navbits_external = 'navbits_external' ELSE data_available = 'ol_data_available' external_navbits_applied = 'external_navbits_ok' navbits_internal = 'navbits_1c' navbits_external = 'navbits_1c' ENDIF ts_olrs = 0.001_wp ! expected (approximate) time in seconds between raw sampling samples ENDIF CALL ncdf_getvar('/quality/'//TRIM(data_available), readbyte1) IF (readbyte1 == 1) THEN ! OL or RS data exist tdir = TRIM(sdir) // 'level_1a/' // TRIM(sampling_mode) // '/' CALL ncdf_getsize(TRIM(tdir)//'dtime', n_olrs, dim = 1) IF ( n_olrs > 1 ) THEN ! read OL or RS data - needs at least 2 samples to make time differences CALL ropp_io_init(data_olrs%Lev1a, n_olrs) ! Time CALL ncdf_getvar(TRIM(tdir)//'dtime', & data_olrs%Lev1a%dtime, units=data_olrs%Lev1a%units%dtime) ! Position and velocity ALLOCATE(r_leo(3,n_olrs), v_leo(3,n_olrs), r_gns(3,n_olrs), v_gns(3,n_olrs)) CALL ncdf_getvar(TRIM(tdir)//'r_transmitter', r_gns, & units=data_olrs%Lev1a%units%r_gns) CALL ncdf_getvar(TRIM(tdir)//'v_transmitter', v_gns, & units=data_olrs%Lev1a%units%v_gns) CALL ncdf_getvar(TRIM(tdir)//'r_receiver', r_leo, & units=data_olrs%Lev1a%units%r_leo) CALL ncdf_getvar(TRIM(tdir)//'v_receiver', v_leo, & units=data_olrs%Lev1a%units%v_leo) data_olrs%Lev1a%r_gns = TRANSPOSE(r_gns) data_olrs%Lev1a%v_gns = TRANSPOSE(v_gns) data_olrs%Lev1a%r_leo = TRANSPOSE(r_leo) data_olrs%Lev1a%v_leo = TRANSPOSE(v_leo) DEALLOCATE(r_gns, v_gns, r_leo, v_leo) ! Coordinate transformation to True-of-Date system, not (yet) using SOFA library data_olrs%Lev1a%r_gns = eci2eci(data%DTocc%year, data%DTocc%month, & data%DTocc%day, data%DTocc%hour, & data%DTocc%minute,data%DTocc%second, & data_olrs%Lev1a%dtime + & data%DTocc%msec/1000.0_wp, & data_olrs%Lev1a%r_gns) data_olrs%Lev1a%v_gns = eci2eci(data%DTocc%year, data%DTocc%month, & data%DTocc%day, data%DTocc%hour, & data%DTocc%minute,data%DTocc%second, & data_olrs%Lev1a%dtime + & data%DTocc%msec/1000.0_wp, & data_olrs%Lev1a%v_gns) data_olrs%Lev1a%r_leo = eci2eci(data%DTocc%year, data%DTocc%month, & data%DTocc%day, data%DTocc%hour, & data%DTocc%minute,data%DTocc%second, & data_olrs%Lev1a%dtime + & data%DTocc%msec/1000.0_wp, & data_olrs%Lev1a%r_leo) data_olrs%Lev1a%v_leo = eci2eci(data%DTocc%year, data%DTocc%month, & data%DTocc%day, data%DTocc%hour, & data%DTocc%minute,data%DTocc%second, & data_olrs%Lev1a%dtime + & data%DTocc%msec/1000.0_wp, & data_olrs%Lev1a%v_leo) ! Phase and amplitude ALLOCATE(navbit_int(n_olrs)) ALLOCATE(navbit_ext(n_olrs)) ALLOCATE(tracking_state(n_olrs)) ! CALL ncdf_getvar(TRIM(tdir)//'navbits_internal', navbit_int) CALL ncdf_getvar('/quality/'//TRIM(external_navbits_applied), have_nb) CALL ncdf_getvar(TRIM(tdir)//'tracking_state', tracking_state) ! IF ( have_nb > 0 ) THEN ! CALL ncdf_getvar(TRIM(tdir)//'navbits_external', navbit_ext) ! ELSE ! navbit_ext(:) = navbit_int(:) ! ENDIF ALLOCATE(phase_l1_iq(n_olrs)) ALLOCATE(i_ca_uncorr(n_olrs)) ALLOCATE(q_ca_uncorr(n_olrs)) ALLOCATE(exphase_l1_nco(n_olrs)) ! FIXME: In newer format versions, the Is and Qs are corrected for navbits; ! the uncorrected versions are no longer available. ! Solution: Use the corrected ones for the newer files, and then set the internal and external ! navbits to zero, which should prevent any later (mis-)corrections being applied. ! Continue to use uncorrected ones with the older files, for backward compatibility. ! We hope and expect that in this case the ropp_pp pre-processing will correctly ! handle the pi phase shifts that result from accumulating the uncorrected Is and Qs. ! If getdirect is in force, phase_L1 will be taken directly from the file anyway. ! In earlier files, open loop snr_L1ca will be overwritten by a direct calculation ! of the SNR derived from navbit-corrected raw sampling Is and Qs. IF (data%FmtVersion < '13.0') THEN CALL ncdf_getvar(TRIM(tdir)//'i_ca_uncorr', i_ca_uncorr) CALL ncdf_getvar(TRIM(tdir)//'q_ca_uncorr', q_ca_uncorr) CALL ncdf_getvar(TRIM(tdir)//'exphase_l1_nco', exphase_l1_nco, & units=data_olrs%Lev1a%units%phase) CALL ncdf_getvar(TRIM(tdir)//'snr_ca', data_olrs%Lev1a%snr_L1ca, & units=data_olrs%Lev1a%units%snr) ! open loop snr_L1ca will be overwritten below ELSE CALL ncdf_getvar(TRIM(tdir)//'i_1c', i_ca_uncorr) ! Despite appearances, these are NOT uncorrected for navigation bits CALL ncdf_getvar(TRIM(tdir)//'q_1c', q_ca_uncorr) ! Despite appearances, these are NOT uncorrected for navigation bits CALL ncdf_getvar(TRIM(tdir)//'exphase_1c_nco', exphase_l1_nco, & units=data_olrs%Lev1a%units%phase) CALL ncdf_getvar(TRIM(tdir)//'snr_1c', data_olrs%Lev1a%snr_L1ca, & units=data_olrs%Lev1a%units%snr) navbit_ext(:) = 0.0_wp navbit_int(:) = 0.0_wp ENDIF DO j=1,n_olrs phase_l1_iq(j) = ATAN2(q_ca_uncorr(j), i_ca_uncorr(j)+TINY(1.0_wp)) ENDDO CALL Accumulate_Phase(phase_l1_iq) data_olrs%lev1a%phase_l1(:) = exphase_l1_nco(:) - & (c_light/(2.0_wp*pi1*data%occfreq1))*phase_l1_iq(:) ! Read the excess phase directly from the EUMETSAT file (if requested) instead ! of deriving from I and Q. This will overwrite phase_L1 and navbits set above. IF (getdirect) THEN IF (data%FmtVersion < '13.0') THEN CALL ncdf_getvar(TRIM(tdir)//'exphase_ca', & data_olrs%Lev1a%phase_L1, units=data_olrs%Lev1a%units%phase) navbit_ext(:) = 0.0_wp navbit_int(:) = 0.0_wp ELSE CALL ncdf_getvar(TRIM(tdir)//'exphase_1c', & data_olrs%Lev1a%phase_L1, units=data_olrs%Lev1a%units%phase) navbit_ext(:) = 0.0_wp navbit_int(:) = 0.0_wp ENDIF ENDIF ! Do a proper average of SNRs (via I and Q). This was needed to fix a temporary ! problem with open loop SNRs at EUMETSAT, which has now been solved. Samples ! are averaged in chunks of 20, unless there is a data gap. When there is a gap, ! only the number of samples up to the gap is averaged and a new average ! is started thereafter. The last average is discarded (going only to n_rs-1). ! This seems to be the way EUMETSAT now constructs open loop samples. ! Note that this calculation uses the navbit-corrected raw sampling Is and Qs, i_ca and q_ca. ! ! EUMETSAT fixed this at PPF version 4.5.1, but this version number is not held in ! the input file. The simplest solution is to use the usual format version ! criterion, especially as this averaging mirrors EUMETSAT's very closely anyway. IF (sampling_mode == 'open_loop' .AND. data%FmtVersion < '13.0') THEN ALLOCATE(i_ca_avr(n_olrs)) ALLOCATE(q_ca_avr(n_olrs)) CALL ncdf_getsize(TRIM(sdir)//'level_1a/raw_sampling/dtime', n_rs, dim = 1) ALLOCATE(dtime(n_rs)) ALLOCATE(i_ca(n_rs)) ALLOCATE(q_ca(n_rs)) ALLOCATE(snr_ca(n_rs)) CALL ncdf_getvar(TRIM(sdir)//'level_1a/raw_sampling/dtime', dtime) CALL ncdf_getvar(TRIM(sdir)//'level_1a/raw_sampling/i_ca', i_ca) CALL ncdf_getvar(TRIM(sdir)//'level_1a/raw_sampling/q_ca', q_ca) CALL ncdf_getvar(TRIM(sdir)//'level_1a/raw_sampling/snr_ca', snr_ca, & units=data_olrs%Lev1a%units%snr) j1=0 j2=0 i_ca_sum=0.0_wp q_ca_sum=0.0_wp dtime_sum=0.0_wp nidx=MAX(1,SUM(MINLOC(dtime, MASK=(dtime > data_olrs%Lev1a%dtime(1))))-10) DO j=nidx,n_rs-1 j2=j2+1 angle_iq=ATAN2(q_ca(j),i_ca(j)+TINY(1.0_wp)) i_ca_sum=i_ca_sum+snr_ca(j)*COS(angle_iq) q_ca_sum=q_ca_sum+snr_ca(j)*SIN(angle_iq) dtime_sum=dtime_sum+dtime(j) IF (j2==20 .OR. dtime(j+1)-dtime(j) > 0.00105_wp) THEN j1=j1+1 i_ca_avr(j1)=i_ca_sum/j2 q_ca_avr(j1)=q_ca_sum/j2 IF (j1==n_olrs) EXIT j2=0 i_ca_sum=0.0_wp q_ca_sum=0.0_wp dtime_sum=0.0_wp ENDIF ENDDO IF (ABS(dtime_sum/j2-data_olrs%Lev1a%dtime(j1)) > 0.00001_wp) THEN WRITE(cval, "(i5,a5,i5)") j1," and ",n_olrs CALL message(msg_warn, "Open loop time does not match - last indices are " // & TRIM(ADJUSTL(cval))) ENDIF DO j=1,n_olrs data_olrs%Lev1a%snr_L1ca(j)=SQRT(i_ca_avr(j)**2+q_ca_avr(j)**2) ENDDO DEALLOCATE(snr_ca) DEALLOCATE(q_ca) DEALLOCATE(i_ca) DEALLOCATE(dtime) DEALLOCATE(q_ca_avr) DEALLOCATE(i_ca_avr) ENDIF ! For setting occultations remove all data with tracking state = 2 at the ! beginning of the record as these are not all valid and may in the processing ! result in cutting off the whole rs/ol section. For rising occultations it is ! better to include them, since the processing will find out where to cut off the ! bad ones at the bottom and thereby keep the good ones with tracking state = 2. IF ( .NOT. BTEST(data%PCD, PCD_rising) ) THEN ! Setting occultation first_valid = -1 ! value to mark that this hasn't been set yet DO j=1,n_olrs IF ((tracking_state(j) /= 2) .AND. (first_valid == -1)) THEN first_valid = j ENDIF ENDDO CALL ropp_io_shrink(data_olrs%Lev1a, first_valid, n_olrs, 1) ALLOCATE(temparray(data_olrs%Lev1a%Npoints)) temparray = navbit_int(first_valid:n_olrs) DEALLOCATE(navbit_int) ALLOCATE(navbit_int(data_olrs%Lev1a%Npoints)) navbit_int = temparray temparray = navbit_ext(first_valid:n_olrs) DEALLOCATE(navbit_ext) ALLOCATE(navbit_ext(data_olrs%Lev1a%Npoints)) navbit_ext = temparray DEALLOCATE(temparray) n_olrs = data_olrs%Lev1a%Npoints ENDIF ! Setting occultation ! data flag ALLOCATE(olrs_LCF(n_olrs)) olrs_LCF(:) = 0 ts1_max = 1.05*ts_olrs IF (BTEST(data%PCD, PCD_rising)) THEN ! Rising occultation DO j=n_olrs-1,1,-1 ts1 = ABS(data_olrs%Lev1a%dtime(j) - data_olrs%Lev1a%dtime(j+1)) IF (ts1 > ts1_max) THEN olrs_LCF(j) = IBSET(olrs_LCF(j), 3) olrs_LCF(j+1) = IBSET(olrs_LCF(j+1),3) ENDIF ENDDO ELSE ! Setting occultation DO j=2,n_olrs ts1 = ABS(data_olrs%Lev1a%dtime(j) - data_olrs%Lev1a%dtime(j-1)) IF (ts1 > ts1_max) THEN olrs_LCF(j) = IBSET(olrs_LCF(j), 3) olrs_LCF(j-1) = IBSET(olrs_LCF(j-1),3) ENDIF ENDDO ENDIF olrs_LCF(:) = IBSET(olrs_LCF(:), 0) ! Open loop mode WHERE (NINT(navbit_ext(:)) == 1) ! External navbit olrs_LCF(:) = IBSET(olrs_LCF(:), 1) ENDWHERE olrs_LCF(:) = IBSET(olrs_LCF(:), 2) ! Navbit quality OK WHERE (NINT(navbit_int(:)) == 1) ! Alternative navbit olrs_LCF(:) = IBSET(olrs_LCF(:), 4) ENDWHERE olrs_LCF(:) = IBSET(olrs_LCF(:), 5) ! Alternative navbit quality ! Setting flag for duplicated CL record WHERE((data_olrs%Lev1a%dtime(1) < data_cl%Lev1a%dtime(:)) .AND. & (data_cl%Lev1a%dtime(:) < data_olrs%Lev1a%dtime(n_olrs))) cl_LCF(:) = IBSET(cl_LCF(:), 6) ENDWHERE ELSE n_olrs = 0 ENDIF ! n_olrs > 1 ELSE CALL message(msg_info, TRIM(sampling_mode) // & ' data requested but are not available in file') ENDIF ! /quality/ == 1 ENDIF ! getlevel1a == 'cl+ol' .OR. getlevel1a == 'cl+rs' ! 9.7.3 Combine CL and RS Level1a variables (if requested) ! --------------------------------------------------------- n = n_cl + n_olrs write(*, *) 'n_cl, n_olrs, n', n_cl, n_olrs, n IF ( n /= 0 ) THEN ! Initialise ro_data structure variables CALL ropp_io_init(data%Lev1a, n) ALLOCATE(lcf(n)) ! Set output variables IF (BTEST(data%PCD, PCD_rising)) THEN ! Rising occultation DO j=1,n_cl data%lev1a%dtime(n_olrs+j) = data_cl%lev1a%dtime(j) data%lev1a%r_gns(n_olrs+j,:) = data_cl%lev1a%r_gns(j,:) data%lev1a%v_gns(n_olrs+j,:) = data_cl%lev1a%v_gns(j,:) data%lev1a%r_leo(n_olrs+j,:) = data_cl%lev1a%r_leo(j,:) data%lev1a%v_leo(n_olrs+j,:) = data_cl%lev1a%v_leo(j,:) data%lev1a%snr_L1ca(n_olrs+j) = data_cl%lev1a%snr_L1ca(j) data%lev1a%snr_L1p(n_olrs+j) = data_cl%lev1a%snr_L1p(j) data%lev1a%snr_L2p(n_olrs+j) = data_cl%lev1a%snr_L2p(j) data%lev1a%phase_L1(n_olrs+j) = data_cl%lev1a%phase_L1(j) data%lev1a%phase_L2(n_olrs+j) = data_cl%lev1a%phase_L2(j) lcf(n_olrs+j) = cl_LCF(j) ENDDO IF (n_olrs /= 0) THEN DO j=1,n_olrs data%lev1a%dtime(j) = data_olrs%lev1a%dtime(j) data%lev1a%r_gns(j,:) = data_olrs%lev1a%r_gns(j,:) data%lev1a%v_gns(j,:) = data_olrs%lev1a%v_gns(j,:) data%lev1a%r_leo(j,:) = data_olrs%lev1a%r_leo(j,:) data%lev1a%v_leo(j,:) = data_olrs%lev1a%v_leo(j,:) data%lev1a%snr_L1ca(j) = data_olrs%lev1a%snr_L1ca(j) data%lev1a%snr_L1p(j) = ropp_MDFV data%lev1a%snr_L2p(j) = ropp_MDFV data%lev1a%phase_L1(j) = data_olrs%lev1a%phase_L1(j) data%lev1a%phase_L2(j) = ropp_MDFV lcf(j) = olrs_LCF(j) ENDDO ! test for gap between cl and (ol or rs) records IF (n_cl > 1) THEN IF (data_olrs%lev1a%dtime(n_olrs) < data_cl%lev1a%dtime(1)-1.5*ts_cl) THEN LCF(n_olrs:n_olrs+1) = IBSET(LCF(n_olrs:n_olrs+1),3) ENDIF ENDIF ! n_cl > 1 ENDIF ! n_olrs /= 0 ELSE ! Setting occultation DO j=1,n_cl data%lev1a%dtime(j) = data_cl%lev1a%dtime(j) data%lev1a%r_gns(j,:) = data_cl%lev1a%r_gns(j,:) data%lev1a%v_gns(j,:) = data_cl%lev1a%v_gns(j,:) data%lev1a%r_leo(j,:) = data_cl%lev1a%r_leo(j,:) data%lev1a%v_leo(j,:) = data_cl%lev1a%v_leo(j,:) data%lev1a%snr_L1ca(j) = data_cl%lev1a%snr_L1ca(j) data%lev1a%snr_L1p(j) = data_cl%lev1a%snr_L1p(j) data%lev1a%snr_L2p(j) = data_cl%lev1a%snr_L2p(j) data%lev1a%phase_L1(j) = data_cl%lev1a%phase_L1(j) data%lev1a%phase_L2(j) = data_cl%lev1a%phase_L2(j) lcf(j) = cl_LCF(j) ENDDO IF (n_olrs /= 0) THEN DO j=1,n_olrs data%lev1a%dtime(n_cl+j) = data_olrs%lev1a%dtime(j) data%lev1a%r_gns(n_cl+j,:) = data_olrs%lev1a%r_gns(j,:) data%lev1a%v_gns(n_cl+j,:) = data_olrs%lev1a%v_gns(j,:) data%lev1a%r_leo(n_cl+j,:) = data_olrs%lev1a%r_leo(j,:) data%lev1a%v_leo(n_cl+j,:) = data_olrs%lev1a%v_leo(j,:) data%lev1a%snr_L1ca(n_cl+j) = data_olrs%lev1a%snr_L1ca(j) data%lev1a%snr_L1p(n_cl+j) = ropp_MDFV data%lev1a%snr_L2p(n_cl+j) = ropp_MDFV data%lev1a%phase_L1(n_cl+j) = data_olrs%lev1a%phase_L1(j) data%lev1a%phase_L2(n_cl+j) = ropp_MDFV lcf(n_cl+j) = olrs_LCF(j) ENDDO IF (n_cl > 1) THEN IF (data_cl%lev1a%dtime(n_cl) < data_olrs%lev1a%dtime(1)-1.5*ts_cl) THEN LCF(n_cl:n_cl+1) = IBSET(LCF(n_cl:n_cl+1),3) ENDIF ENDIF ! n_cl > 1 ENDIF ! n_olrs /= 0 ENDIF ! Rising occultation ! 9.7.4 Missing/invalid data checks ! --------------------------------- WHERE (ropp_io_isnan(data%Lev1a%phase_L1)) data%Lev1a%phase_L1 = ropp_MDFV lcf = IBSET(lcf, 3) ENDWHERE WHERE (ropp_io_isnan(data%Lev1a%phase_L2)) data%Lev1a%phase_L2 = ropp_MDFV ENDWHERE WHERE (ropp_io_isnan(data%Lev1a%snr_L1ca)) data%Lev1a%snr_L1ca = ropp_MDFV lcf = IBSET(lcf, 3) ENDWHERE WHERE (ropp_io_isnan(data%Lev1a%snr_L1p)) data%Lev1a%snr_L1p = ropp_MDFV ENDWHERE WHERE (ropp_io_isnan(data%Lev1a%snr_L2p)) data%Lev1a%snr_L2p = ropp_MDFV ENDWHERE data%Lev1a%reference_frame%r_gns = 'ECI' data%Lev1a%reference_frame%r_leo = 'ECI' data%Lev1a%range%phase = (/ & MIN(MINVAL(data%Lev1a%phase_L1), data%Lev1a%range%phase(1)), & MAX(MAXVAL(data%Lev1a%phase_L1), data%Lev1a%range%phase(2)) /) ! Add lost carrier information to file CALL ropp_io_addvar_rodataD1d(data, & name = 'open_loop_lcf', & long_name= 'Lost Carrier Flag', & units = '', & range = (/-1000000.0_wp, 1000000.0_wp/), & DATA = REAL(lcf,wp) ) DEALLOCATE(lcf) IF ( n_cl /= 0 ) THEN DEALLOCATE(cl_lcf) CALL ropp_io_free(data_cl) ENDIF IF ( n_olrs /= 0 ) THEN DEALLOCATE(olrs_lcf) CALL ropp_io_free(data_olrs) ENDIF ENDIF ! n /= 0 ENDIF ! getlevel1a /= 'none' ! 9.8 Level1b variables (if any and if requested) ! ----------------------------------------------- tdir = '/quality/' IF (data%FmtVersion < '13.0') THEN CALL ncdf_getvar(TRIM(tdir)//'fsi_done', readbyte1) CALL ncdf_getvar(TRIM(tdir)//'go_done', readbyte2) ELSE ! in version 13 and above there are no fsi_done or go_done quality flags ! assume they are set to 1 readbyte1 = 1 readbyte2 = 1 ENDIF IF ( (readbyte1 == 1 .OR. readbyte2 == 1) .AND. (getlevel1b .OR. getiono) ) THEN IF (ncdf_isvar(TRIM(ddir)//'impact')) THEN CALL ncdf_getsize(TRIM(ddir)//'impact', n, dim = 1) CALL ropp_io_init(data%Lev1b, n) ELSE data%Lev1b%Npoints = 0 ENDIF IF (data%Lev1b%Npoints > 0) THEN CALL ncdf_getvar(TRIM(ddir)//'azimuth_tp', data%Lev1b%azimuth_tp, & units=data%Lev1b%units%azimuth_tp) CALL ncdf_getvar(TRIM(ddir)//'impact', data%Lev1b%impact, & units=data%Lev1b%units%impact) CALL ncdf_getvar(TRIM(ddir)//'impact', data%Lev1b%impact_L1, & units=data%Lev1b%units%impact) CALL ncdf_getvar(TRIM(ddir)//'impact', data%Lev1b%impact_L2, & units=data%Lev1b%units%impact) ! CALL ncdf_getvar(TRIM(ddir)//'impact_height', data%Lev1b%impact_height, & ! units=data%Lev1b%units%impact) IF ( getiono ) THEN data%Lev1b%range%bangle(1) = -data%Lev1b%range%bangle(2) ! to allow negative values to pass range-checking CALL ncdf_getvar(TRIM(ddir)//'lat_tp_l1', data%Lev1b%lat_tp, & units=data%Lev1b%units%lat_tp) CALL ncdf_getvar(TRIM(ddir)//'lon_tp_l1', data%Lev1b%lon_tp, & units=data%Lev1b%units%lon_tp) CALL ncdf_getvar(TRIM(ddir)//'bangle_l1', data%Lev1b%bangle_L1, & units=data%Lev1b%units%bangle) CALL ncdf_getvar(TRIM(ddir)//'bangle_l1_sdev', data%Lev1b%bangle_L1_sigma, & units=data%Lev1b%units%bangle) CALL ncdf_getvar(TRIM(ddir)//'bangle_l2', data%Lev1b%bangle_L2, & units=data%Lev1b%units%bangle) CALL ncdf_getvar(TRIM(ddir)//'bangle_l2_sdev', data%Lev1b%bangle_L2_sigma, & units=data%Lev1b%units%bangle) ELSE ! CALL ncdf_getvar(TRIM(ddir)//'dtime_mean', data%Lev1b%dtime_mean, & ! units=data%Lev1b%units%dtime_mean) CALL ncdf_getvar(TRIM(ddir)//'lat_tp', data%Lev1b%lat_tp, & units=data%Lev1b%units%lat_tp) CALL ncdf_getvar(TRIM(ddir)//'lon_tp', data%Lev1b%lon_tp, & units=data%Lev1b%units%lon_tp) CALL ncdf_getvar(TRIM(ddir)//'bangle', data%Lev1b%bangle, & units=data%Lev1b%units%bangle) IF (data%FmtVersion < '13.0') THEN CALL ncdf_getvar(TRIM(ddir)//'bangle_ca', data%Lev1b%bangle_L1, & units=data%Lev1b%units%bangle) ELSE CALL ncdf_getvar(TRIM(ddir)//'bangle_l1', data%Lev1b%bangle_L1, & units=data%Lev1b%units%bangle) ENDIF IF (getextrap) THEN IF (data%FmtVersion < '13.0') THEN CALL ncdf_getvar(TRIM(ddir)//'bangle_ca_p2_diff', data%Lev1b%bangle_L2, & ! placeholder, = L1 minus L2 units=data%Lev1b%units%bangle) ELSE CALL ncdf_getvar(TRIM(ddir)//'bangle_l4', data%Lev1b%bangle_L2, & ! placeholder, = L1 minus L2 units=data%Lev1b%units%bangle) ENDIF data%Lev1b%bangle_L2 = data%Lev1b%bangle_L1 - data%Lev1b%bangle_L2 ! L2 = L1 - (L1 - L2) ELSE IF (data%FmtVersion < '13.0') THEN CALL ncdf_getvar(TRIM(ddir)//'bangle_p2', data%Lev1b%bangle_L2, & units=data%Lev1b%units%bangle) ELSE IF ( data%leo_id == 'MESG' ) THEN CALL ncdf_getvar(TRIM(ddir)//'bangle_l5', data%Lev1b%bangle_L2, & units=data%Lev1b%units%bangle) ELSE CALL ncdf_getvar(TRIM(ddir)//'bangle_l2', data%Lev1b%bangle_L2, & units=data%Lev1b%units%bangle) ENDIF ENDIF ENDIF ENDIF ! remove all NaN from EUM fields ! FIXME: maybe do this in ncdf_getvar? WHERE( .NOT. ropp_io_isnan(data%Lev1b%bangle) ) data%Lev1b%bangle_qual = 100.0_wp WHERE( ropp_io_isnan(data%Lev1b%lat_tp) ) data%Lev1b%lat_tp = ropp_MDFV WHERE( ropp_io_isnan(data%Lev1b%lon_tp) ) data%Lev1b%lon_tp = ropp_MDFV WHERE( ropp_io_isnan(data%Lev1b%azimuth_tp) ) data%Lev1b%azimuth_tp = ropp_MDFV WHERE( ropp_io_isnan(data%Lev1b%bangle) ) data%Lev1b%bangle = ropp_MDFV WHERE( ropp_io_isnan(data%Lev1b%bangle_L1) ) data%Lev1b%bangle_L1 = ropp_MDFV WHERE( ropp_io_isnan(data%Lev1b%bangle_L1_sigma) ) data%Lev1b%bangle_L1_sigma = ropp_MDFV WHERE( ropp_io_isnan(data%Lev1b%bangle_L2) ) data%Lev1b%bangle_L2 = ropp_MDFV WHERE( ropp_io_isnan(data%Lev1b%bangle_L2_sigma) ) data%Lev1b%bangle_L2_sigma = ropp_MDFV ! set the quality for bangle, EUMETSAT data uses NaN for missing data%Lev1b%Missing = .FALSE. ! Why is this here and not in Sec 9.4? Maybe because it makes sense only for level 1b ! Are open loop data used? IF ( ncdf_isvar(TRIM(tdir)//'ol_data_available') ) THEN CALL ncdf_getvar(TRIM(tdir)//'ol_data_available', readbyte1) ELSE SELECT CASE (data%leo_id) CASE ("SE6A", "SE6B", "MESG") readbyte1 = 1 ! Sentinel-6 has no ol_data_available flag CASE DEFAULT IF (data%leo_id(1:1) == "L") THEN ! This is a Lemur/Spire ID. If there are other missions with ! an ID starting with 'L', they should be listed explicitly above. readbyte1 = 1 ! Spire has no ol_data_available flag ELSE readbyte1 = 0 ENDIF END SELECT ENDIF IF (data%FmtVersion < '13.0') THEN CALL ncdf_getvar(TRIM(tdir)//'rs_data_available', readbyte2) ELSE ! in version 13 and above there is no rs_data_available flag ! assume it is the same as ol_data_available readbyte2 = readbyte1 ENDIF IF ( readbyte1 == 1 .OR. readbyte2 == 1) THEN IF ( ncdf_isvar(TRIM(tdir)//'ol_data_used') ) THEN CALL ncdf_getvar(TRIM(tdir)//'ol_data_used', readbyte3) IF ( readbyte3 == 1 ) data%PCD = IBSET(data%PCD, PCD_open_loop) ENDIF ENDIF ENDIF ! data%Lev1b%Npoints > 0 ELSE ! not (fsi_done or go_done) data%Lev1b%Npoints = 0 IF (getlevel1b) THEN CALL message(msg_warn, 'Bending angle data requested ' // & 'but are not available in file') ENDIF ENDIF ! 9.9 (Global) Attributes ! ----------------------- data%processing_centre = ' ' CALL ncdf_getatt('institution', data%processing_centre) data%pod_method = ' ' data%phase_method = ' ' data%bangle_method = ' ' gdir = TRIM(sdir) // 'occultation/' IF ( getiono ) gdir = TRIM(sdir) // 'ionosphere/occultation/' IF (data%FmtVersion < '13.0') THEN CALL ncdf_getatt(TRIM(gdir)//'pod_method', data%pod_method) CALL ncdf_getatt(TRIM(gdir)//'phase_method', data%phase_method) CALL ncdf_getatt(TRIM(gdir)//'retrieval_method', data%bangle_method) ELSE CALL ncdf_getvar(TRIM(gdir)//'pod_method', data%pod_method) CALL ncdf_getvar(TRIM(gdir)//'phase_method', data%phase_method) CALL ncdf_getvar(TRIM(gdir)//'retrieval_method', data%bangle_method) ENDIF data%refrac_method = 'UNKNOWN' data%meteo_method = 'UNKNOWN' IF ( getlevel1b .OR. getiono ) THEN IF (TRIM(resolution) == 'high_resolution') THEN data%thin_method = 'Unthinned data' ELSE data%thin_method = ' ' IF ( getiono ) THEN CALL ncdf_getatt(TRIM(sdir)//'ionosphere/thinned/thinner_method', data%thin_method) ELSE !CALL ncdf_getatt(TRIM(ddir)//'thinner_method', data%thin_method) data%thin_method = 'UNKNOWN' ENDIF ENDIF ELSE data%thin_method = 'UNKNOWN' ENDIF ! Software version is the ROPP software version, in the format 'vnn.mmm'. ! The last digit of the fractional part (mmm) is reserved for indicating ! the bending angle processing method. ! At ROPP8.0, a new variable, 'processing_software', has been introduced, ! which can hold information about other software - in this case, the EUM ! processing code that generated the data in the first place. data%software_version = ' ' ; readstr = ' ' ; readstr2 = ' ' CALL ncdf_getatt('/status/processing/processor_name' , readstr) CALL ncdf_getatt('/status/processing/processor_version' ,readstr2) ! Wipe all characters starting from the first non-digit or non-dot character ! in the version string n = 0 DO k=1,LEN_TRIM(readstr2) IF ( ( n == 1 ) .OR. ( & ( (IACHAR(readstr2(k:k)) < 48) .OR. (IACHAR(readstr2(k:k)) > 57) ) & .AND. & (IACHAR(readstr2(k:k)) /= 46) & ) ) THEN n = 1 readstr2(k:k) = ' ' END IF END DO i = SCAN(readstr2, ".", .FALSE.) j = SCAN(readstr2, ".", .TRUE.) IF (i < 2) THEN CALL message(msg_warn, "Processor version cannot be decoded") ELSE IF (i == j) THEN j = LEN_TRIM(readstr2) + 1 READ(readstr2(1:i-1), FMT=*) readint READ(readstr2(i+1:j-1), FMT=*) readint2 ELSE READ(readstr2(1:i-1), FMT=*) readint READ(readstr2(i+1:j-1), FMT=*) readint2 ENDIF SELECT CASE (TRIM(data%bangle_method)) CASE ("FSI") readint3 = 1 CASE ("GO") readint3 = 0 CASE DEFAULT readint3 = 9 END SELECT WRITE(data%software_version, '("v",i2.2,".",i2.2,i1)') readint, readint2, readint3 ENDIF data%processing_software = TRIM(readstr) // ' ' // & TRIM(readstr2) // ' ' // & TRIM(data%bangle_method) ! Processing time, split up string ! FIXME: Is there a library function for this that does some more format checks? IF (data%FmtVersion < '13.0') THEN CALL ncdf_getatt('/status/processing/creation_time', readstr) SELECT CASE (LEN(TRIM(readstr))) CASE (14) READ(readstr, '(i4,i2,i2,i2,i2,i2)') & data%DTpro%Year, data%DTpro%Month, data%DTpro%Day, & data%DTpro%Hour, data%DTpro%Minute, data%DTpro%Second CASE (17) READ(readstr, '(i4,i2,i2,i2,i2,i2,i3)') & data%DTpro%Year, data%DTpro%Month, data%DTpro%Day, & data%DTpro%Hour, data%DTpro%Minute, data%DTpro%Second, data%DTpro%Msec CASE (19) READ(readstr, '(i4,1x,i2,1x,i2,1x,i2,1x,i2,1x,i2)') & data%DTpro%Year, data%DTpro%Month, data%DTpro%Day, & data%DTpro%Hour, data%DTpro%Minute, data%DTpro%Second CASE (23) READ(readstr, '(i4,1x,i2,1x,i2,1x,i2,1x,i2,1x,i2,1x,i3)') & data%DTpro%Year, data%DTpro%Month, data%DTpro%Day, & data%DTpro%Hour, data%DTpro%Minute, data%DTpro%Second, data%DTpro%Msec CASE DEFAULT CALL message( msg_warn, & 'Invalid /status/processing/creation_time of ' // & TRIM(readstr) // ' ... defaulting DTpro%year etc.' ) END SELECT ELSE ! add here parsing of JS since 2000 CALL ncdf_getvar('/status/processing/creation_time_utc', readreal) ! convert creation time from JS2000 to DT8 CALL TimeSince ( DT8, readreal, 0, Base="JS2000" ) data%DTpro%Year = DT8(1) data%DTpro%Month = DT8(2) data%DTpro%Day = DT8(3) data%DTpro%Hour = DT8(5) data%DTpro%Minute = DT8(6) data%DTpro%Second = DT8(7) data%DTpro%Msec = DT8(8) ENDIF ! 9.10 Occultation ID ! ------------------- CALL ropp_io_occid(DATA) ! 9.11 Clean up ! ------------- ! IF (ASSOCIATED(idx)) DEALLOCATE(idx) END SUBROUTINE ropp_io_read_eumdata !------------------------------------------------------------------------------- ! 10. Accumulate phase !------------------------------------------------------------------------------- SUBROUTINE Accumulate_Phase(Ph, Sign) ! (Array of (accumulated) phase, dir) ! Method: ! Sign = 0 or no Sign: ! Adding +-2*Pi where phase jumps from ! +-Pi to -+Pi, ! Sign > 0: ! Adding +2*Pi where phase jumps from ! - to + ! Sign < 0 ! Adding -2*Pi where phase jumps from ! + to - ! 10.1 Declarations USE typesizes, ONLY: wp => EightByteReal IMPLICIT NONE REAL(wp), DIMENSION(:), INTENT(inout) :: Ph ! Phase --> accumulated phase INTEGER, OPTIONAL, INTENT(in) :: Sign ! Phase change sign ! Pi already defined as parameter in coordinates module. Confuses pgf95. REAL(wp), PARAMETER :: pi1 = 3.141592653589793238_wp INTEGER :: i ! Array index INTEGER :: PSign ! Phase change sign ! 10.2 Determine phase change sign IF (.NOT. PRESENT(Sign)) THEN PSign = 0 ELSE PSign = Sign ENDIF ! 10.3 Accumulate phase IF (PSign == 0) THEN DO i=2,SIZE(Ph) Ph(i) = Ph(i-1) + MODULO(Ph(i)-Ph(i-1)+pi1, 2.0_wp*pi1) - pi1 ENDDO ELSEIF (PSign > 0) THEN DO i=2,SIZE(Ph) Ph(i) = Ph(i-1) + MODULO(Ph(i)-Ph(i-1), 2.0_wp*pi1) ENDDO ELSEIF (PSign < 0) THEN DO i=2,SIZE(Ph) Ph(i) = Ph(i-1) + MODULO(Ph(i)-Ph(i-1)+2.0_wp*pi1, 2.0_wp*pi1) - 2.0_wp*pi1 ENDDO ENDIF END SUBROUTINE Accumulate_Phase !------------------------------------------------------------------------------- ! 11. Handle time format of EUM files !------------------------------------------------------------------------------- SUBROUTINE abstimetoDT(i, r, DTocc) ! Based on 'Practical Ephemeris Calculations' by Oliver Montenbruck ! (Springer-Verlag, 1989). Added handling of hour, seconds as well. ! ! This returns the correct UTC, accounting for leapseconds, because it ! assumes all days have 86400 seconds. If there are leapseconds, some of ! the days will have 86401 secs, but then the corresponding UTC will be ! delayed by the same amount, so the resulting DTocc will be the same. USE ropp_io_types, ONLY: DT7type IMPLICIT NONE INTEGER, INTENT(in) :: i ! number of days since 2000-01-01 00:00:00 REAL(KIND(1.0D0)), INTENT(in) :: r ! seconds since start of day TYPE(DT7type), INTENT(inout) :: DTocc REAL(KIND(1.0D0)) :: h, mi, s, ms INTEGER :: a, b, c, d, e, f, y, m, dd, jd ! calculate JD from input, 2451545 is JD of 2000-01-01 jd = i + 2451545 a = INT(jd+0.5D0) IF (a < 2299161) THEN c = a + 1524 ELSE b = INT( (a - 1867216.25D0) / 36524.25D0 ) c = a + b - INT(b/4D0) + 1525 ENDIF d = INT( ( c-122.1D0 ) / 365.25D0 ) e = INT( 365.25D0*d) f = INT( (c-e) / 30.6001D0 ) dd = c - e - INT( 30.6001*f ) + MOD( ( jd+0.5D0 ), DBLE(a) ) m = f - 1 - 12*INT( f/14D0 ) y = d - 4715 - INT( ( 7+m )/10.D0 ) DTocc%Year = y DTocc%Month = m DTocc%Day = dd h = FLOOR(r/3600.D0) DTocc%Hour = INT(h) mi = FLOOR( (r - h*3600.D0) / 60.D0 ) s = FLOOR( (r - h*3600.D0 - mi*60.D0 ) ) ms = (r - h*3600.D0 - mi*60.D0 - s) * 1000.D0 DTocc%Hour = INT(h) DTocc%Minute = INT(mi) DTocc%Second = INT(s) DTocc%Msec = INT(ms) END SUBROUTINE abstimetoDT !------------------------------------------------------------------------------- ! 12. Find where reals are NaN (used by EUMETSAT to indicate missing data) !------------------------------------------------------------------------------- FUNCTION ropp_io_isnan(x) RESULT(lnan) ! Says where reals are NaNs. USE typesizes, ONLY: wp => EightByteReal IMPLICIT NONE REAL(wp), DIMENSION(:), INTENT(IN) :: x LOGICAL, DIMENSION(SIZE(x)) :: lnan INTEGER :: k lnan(:) = .FALSE. ! g95 doesn't like this. ! WHERE ( & ! (x /= x) .OR. & ! (x + 1.0_wp == x) .OR. & ! ((x > 0) .EQV. (x <= 0)) & ! ) lnan = .TRUE. DO k=1,SIZE(x) IF ( & (x(k) /= x(k)) .OR. & (x(k) + 1.0_wp == x(k)) .OR. & ((x(k) > 0) .EQV. (x(k) <= 0)) & ) lnan(k) = .TRUE. END DO END FUNCTION ropp_io_isnan END SUBROUTINE ropp_io_read_ncdf_get_eumdata