! $Id: ropp_pp_occ_tool.f90 2021 2009-01-16 10:49:04Z frhl $ PROGRAM ropp_pp_occ_tool !****p* Programs/ropp_pp_occ_tool * ! ! NAME ! ropp_pp_occ_tool ! ! SYNOPSIS ! Pre-processing tool to calculate refractivity and dry temperature ! profile from L1 and L2 excess phase radio occultation data using ! geometric optics or wave optics, ionospheric correction, statistical ! optimization, Abel transform and hydrostatic equation. ! ! > ropp_pp_occ_tool [] ! ! ARGUMENTS ! One (or more) input file names. ! ! OPTIONS ! -c name of configuration file ! -o name of ROPP netCDF output file ! -m ionospheric correction method ! [NONE,MSIS,GMSIS,BG,BARO,GBARO], (default GMSIS) ! -mfile model coefficients file ! [MSIS_coeff.nc,BAROCLIM_coeff.nc], (default MSIS_coeff.nc) ! -bfile background model atmospheric profile file path ! (if using BG ionospheric correction method) ! -navfile external navigation bit *_txt file path ! (default internal correction) ! -occ processing method, WO or GO (default WO) ! -filter filtering method, slpoly or optest ! (default slpoly, sliding polynomial) ! -fit apply 2-parameter regression fit to model ! -ellipsoid output height with respect to WGS84 ellipsoid ! (default output wrt EGM96 geoid) ! -full output whole profile ! (default: only output data within the observed range) ! -w process output from ropp_pp_wopt_tool ! -d output additional diagnostics ! -h help ! -v version information ! ! DESCRIPTION ! This program reads RO L1 and L2 excess phase data as a function of time ! from the input data files and calculates vertical profiles L1 and L2 ! bending angles, ionospheric corrected bending angle, refractivity and ! dry temperature. The result is written to a ROPP output file. ! ! NOTES ! If the input file is a multifile, or more than one input files are ! specified, the output file is a multifile. ! ! Already existing output files will be overwritten. ! ! ERRORS ! Program (shell) return codes: ! 0 = OK ! 1 = At least one Warning occurred ! 2 = At least one Error occurred ! 3 = A Fatal error occurred ! ! EXAMPLE ! To calculate bending angle, refractivity and dry temperature from one ! of the example (single-) files in the data directory: ! ! > ropp_pp_occ_tool ../data/input.nc -o example_01.nc ! ! To calculate bending angle, refractivity and dry temperature profiles ! from all singlefiles in the data directory: ! ! > ropp_pp_occ_tool ../data/*.nc -o example_02.nc ! ! Note that the resulting example_02.nc file contains processed data from ! all example profiles. ! ! To calculate bending angle, refractivity and dry temperature profiles ! from all profiles contained in the multifile multi.nc: ! ! > ropp_pp_occ_tool ../data/multi.nc -o example_03.nc ! ! Since the multi_* file was generated by concatenating the other files ! in the data directory, example_02.nc and example_03.nc should be identical ! apart from the file names. ! ! 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. Declarations !------------------------------------------------------------------------------- USE typesizes, ONLY: wp => EightByteReal USE ropp_utils USE ropp_io USE ropp_io_types, ONLY: ROprof, L1atype, L1btype, L2atype USE ropp_pp USE ropp_pp_preproc USE ropp_pp_types, ONLY: PPConfig, PPDiag IMPLICIT NONE TYPE(ROprof) :: ro_data ! Input RO data TYPE(L1atype) :: Lev1a ! Temporary Level1a structure for storage TYPE(L1btype) :: bangle ! Output bending angle profiles TYPE(L1btype) :: out_ba ! Corrected bending angles TYPE(L1btype) :: smt_ba ! Smoothed bending angle TYPE(L1btype) :: mod_ba ! Model bending angle TYPE(L1btype) :: gnos_ba ! For GNOS data TYPE(L2atype) :: out_refrac ! Retrieved refractivity profile TYPE(L2btype) :: dum_meteo ! Level2b structure for dummy meteo variables TYPE(ppConfig) :: config ! Configuration options TYPE(ppDiag) :: diag ! Diagnostic variables INTEGER :: idummy, imin, imax, i1, i2 INTEGER :: i, iargc, argc, k, j INTEGER :: n_files, n_profiles, nbi INTEGER :: ws_go_smooth, ws_go_full, ws_wo, ws_low REAL(wp) :: Pmin, Pmax REAL(wp),parameter :: GLO_KFreq(24)=[ 1,-4,05,06,01,-4,05,06,-2,-7,00, & -1,-2,-7,00,-1, 4,-3,03,02,04,-3,03,02] REAL(wp),parameter :: DFRQ1_GLO=0.56250e6_wp, baseFL1=1602.00e6_wp; REAL(wp),parameter :: DFRQ2_GLO=0.43750e6_wp, baseFL2=1246.00e6_wp; REAL(wp), PARAMETER :: ih1 = 60.0_wp ! impact height limit for reference bangle in km REAL(wp), PARAMETER :: ih2 = 80.0_wp ! impact height limit for reference bangle in km REAL(wp), PARAMETER :: ih3 = 25.0_wp ! impact height limit for reference bangle in km REAL(wp), PARAMETER :: ih4 = 40.0_wp ! impact height limit for reference bangle in km ! Yong Chen 01/25/2025 ! Parameters for quality control tests based on CDAAC recommended values ! see "Standard RO Inversions in the Neutral Atmosphere 2013 - 2020 (Processing Steps and Explanation of Data) ! by Sergey Sokolovskiy, 2021 REAL(wp), PARAMETER :: th_reldevmax = 0.25_wp ! Maximum relative difference between ! bending angle and climate model between 25 and 40 km REAL(wp), PARAMETER :: th_reldevstd = 0.1_wp !3.0e-5_wp ! Standard deviation of the relative difference ! between bending angle and climate model between 25 and 40 km REAL(wp), PARAMETER :: th_difmaxref = 0.5_wp ! Maximum relative difference between refractivity and ! climate model between 10 and 60 km REAL(wp), PARAMETER :: th_snr1avg = 60.0_wp !200.0_wp ! L1 average SNR from 60 to 80 km (v/v) REAL(wp), PARAMETER :: th_maxdifphase = 0.1_wp ! Maximum absolute difference of L1 and L2 excess phase ! finite differences between 20 and 40 km (m/sample) REAL(wp), PARAMETER :: th_stdv = 4.0e-5_wp ! 1.5e-4_wp ! std bangle diference between obs and background from 60 to 80 km (rad) REAL(wp), PARAMETER :: th_smean = 2.0e-5_wp ! 1.0e-4_wp ! mean bangle diference between obs and background from 60 to 80 km (rad) REAL(wp), ALLOCATABLE :: relative_ba(:) INTEGER :: GLO_SLOT_NUM REAL(wp), ALLOCATABLE :: go_l1_ip(:) INTEGER :: ialt REAL(wp), ALLOCATABLE :: alt_tmp(:) REAL(wp), ALLOCATABLE :: refrac_tmp(:) LOGICAL :: give_help LOGICAL :: earth_ellipsoid = .FALSE. LOGICAL :: output_full = .FALSE. LOGICAL :: twofit = .FALSE. LOGICAL :: ranchk = .TRUE. LOGICAL :: wopt_mode = .FALSE. CHARACTER(len = 4096), DIMENSION(:), ALLOCATABLE :: ifiles CHARACTER(len = 4096) :: ofile CHARACTER(len = 4096) :: cfg_file = " " CHARACTER(len = 256) :: buffer CHARACTER(len = 64) :: outstr CHARACTER(len = 256) :: mfile = " " CHARACTER(len = 256) :: bfile = " " CHARACTER(len = 256) :: navfile = " " CHARACTER(len = 10) :: method = " " CHARACTER(len = 10) :: occmethod = " " CHARACTER(len = 10) :: filtmethod = " " CHARACTER(len = 4) :: istr CHARACTER(len = 10) :: nstr CHARACTER(len = 10) :: pstr1 CHARACTER(len = 10) :: pstr2 CHARACTER(len = 10), DIMENSION(3) :: smess !------------------------------------------------------------------------------- ! 2. Default settings !------------------------------------------------------------------------------- give_help = .FALSE. ofile = "ropp_pp_occ.nc" CALL message(msg_noin, '') CALL message(msg_noin, & '----------------------------------------------------------------------') CALL message(msg_noin, & ' ROPP Occultation Pre-processor Tool') CALL message(msg_noin, & '----------------------------------------------------------------------') CALL message(msg_noin, '') !------------------------------------------------------------------------------- ! 3. Command line arguments !------------------------------------------------------------------------------- argc = iargc() i = 1 n_files = 0 ALLOCATE(ifiles(argc)) DO WHILE(i <= argc) CALL getarg(i, buffer) SELECT CASE (buffer) CASE('-o') ! Output file name (netCDF output) CALL getarg(i+1, buffer) ofile = buffer i = i + 1 CASE('-c') ! Configuration file name CALL getarg(i+1, buffer) cfg_file = buffer i = i + 1 CASE('-m') ! Ionospheric correction method CALL getarg(i+1, buffer) method = buffer i = i + 1 CASE('-mfile') ! Model coefficients file CALL getarg(i+1, buffer) mfile = buffer i = i + 1 CASE('-bfile') ! Background atmosphere profile file CALL getarg(i+1, buffer) bfile = buffer i = i + 1 CASE('-navfile') ! External navigation bit file CALL getarg(i+1, buffer) navfile = buffer i = i + 1 CASE('-occ', '--occ') ! Occultation processing method CALL getarg(i+1, buffer) occmethod = buffer i = i + 1 CASE('-filter', '--filter') ! Filtering method CALL getarg(i+1, buffer) filtmethod = buffer i = i + 1 CASE('-fit', '--fit') ! Use two parameter fitting method twofit = .true. CASE('-ellipsoid') ! Output height wrt reference ellipsoid earth_ellipsoid = .TRUE. CASE('-full', '--full') ! Output full height grid (IC plus climatology) output_full = .TRUE. CASE('-d') ! 'Diagnostic' output mode msg_MODE = VerboseMode CASE('-w') ! Read output from WOPT wopt_mode = .TRUE. CASE ('-no-ranchk') ! Use no rangecheck on output ranchk = .FALSE. CASE('-h', '--help', '?') ! Give some help give_help = .TRUE. CASE('-v', '-V', '--version') ! Give some version information CALL version_info() CALL EXIT(msg_exit_ok) CASE default ! Input file name n_files = n_files + 1 ifiles(n_files) = buffer END SELECT i = i + 1 END DO IF ( n_files == 0 .AND. .NOT. give_help ) THEN CALL message ( msg_error, "No input file(s) specified" ) ENDIF IF (argc == 0 .OR. n_files == 0 .OR. give_help) THEN CALL usage() CALL EXIT(msg_exit_status) ENDIF ! 3.1 Read configuration file (if exists), preserving command-line options IF (cfg_file /= " ") THEN CALL message(msg_info, & "Reading configuration file " // TRIM(cfg_file) // ".\n") CALL ropp_pp_read_config(cfg_file, config) ENDIF IF (method /= " ") config%method = method IF (mfile /= " ") config%mfile = mfile IF (bfile /= " ") config%bfile = bfile IF (navfile /= " ") config%navbit_file = navfile IF (occmethod /= " ") config%occ_method = occmethod IF (filtmethod /= " ") config%filter_method = filtmethod IF (twofit) config%nparm_fit = 2 !------------------------------------------------------------------------------- ! 4. Remove pre-existing output file !------------------------------------------------------------------------------- CALL file_delete(ofile, idummy) !------------------------------------------------------------------------------- ! 5. Loop over all input files and profiles !------------------------------------------------------------------------------- DO k = 1, n_files n_profiles = ropp_io_nrec(ifiles(k)) DO i = 1, n_profiles WRITE(istr, '(i4)') i WRITE(nstr, '(i6)') n_profiles CALL message(msg_info, "Processing profile " // istr // " of " // nstr ) !------------------------------------------------------------------------------- ! 6. Read data !------------------------------------------------------------------------------- ! 6.1 Read input file CALL message(msg_info, & "Reading input data file " // TRIM(ifiles(k)) // ".\n") CALL ropp_io_read(ro_data, ifiles(k), rec=i, ranchk=ranchk) CALL message(msg_info, "(" // TRIM(ro_data%occ_id) // ") \n") IF (ro_data%Lev1a%Npoints == 0) THEN config%obs_ok = .FALSE. CALL ropp_io_free(ro_data%vlist) CALL ropp_io_write(ro_data, ofile, append=.TRUE., ranchk=ranchk ) CALL message(msg_fatal, "FAILURE: No Level1a data in file " // & ifiles(k) // "\n") ENDIF ! 6.2 Shrink ro_data to correct size (for multiple profiles) IF (ro_data%lev1a%npoints == 0) THEN CALL message(msg_fatal, "No Level1a data in file " // & TRIM(ifiles(k)) // ". No data to process. \n") CALL ropp_io_free(ro_data%vlist) CALL ropp_io_write(ro_data, ofile, append=.TRUE., ranchk=ranchk ) EXIT ENDIF CALL ropp_io_init(Lev1a, ro_data%lev1a%npoints) Lev1a = ro_data%lev1a ro_data%lev1a = Lev1a CALL ropp_io_free(Lev1a) CALL ropp_io_free(ro_data%Lev1b) ! 6.21 assign frequencies accoording to the input ro_data ! Yong Chen 08/09/2024 IF (ro_data%occfreq1 > ropp_MDTV .and. ro_data%occfreq2 > ropp_MDTV ) THEN f_L1 = ro_data%occfreq1 f_L2 = ro_data%occfreq2 ELSE ! 6.2.2 Change Frequencies acoording to the GNLONASS satellite ! Yong Chen 08/22/2021 ! Here, we need to update the GNSS frequency since the ropp_pp_constants sets ! f_L1 and f_L2 for GPS only, the linear combination using GPS frequencies ! will not yield incorrect results, so update here (only once). print *, ro_data%gns_id IF (INDEX(ro_data%gns_id,"R") ==1 ) THEN read(ro_data%gns_id,'(1X,I03)')GLO_SLOT_NUM f_L1=baseFL1+GLO_KFreq(GLO_SLOT_NUM)*DFRQ1_GLO f_L2=baseFL2+GLO_KFreq(GLO_SLOT_NUM)*DFRQ2_GLO print *, 'GLONASS:F_L1==',f_L1 print *, 'GLONASS:F_L2==',f_L2 ELSE IF (INDEX(ro_data%gns_id,"E") ==1 ) THEN f_L1=1575.420d6 ! Note the Galeo for SPIRE uses E1 as L1, and L7Q (E5b) as L2. f_L2=1207.140d6 ! ! Updated the L1 and L2 frequencies for Spire and PlanetIQ ! Yong Chen 02/06/2024 ELSE IF (INDEX(ro_data%gns_id,"C") ==1 ) THEN IF (ro_data%leo_id(1:2)=='S0' .OR. ro_data%leo_id(1:2)=='S1' & .OR. ro_data%leo_id(1:2)=='S2' ) THEN !f_L1=1561.098d6 ! Note BDS for SPIRE uses B1C_P as L1, and B2b as L2. f_L1=1575.420d6 ! Yong Chen on 08/23/2024 f_L2=1207.140d6 ELSE IF (ro_data%leo_id(1:2)=='GN' ) THEN f_L1=1575.420d6 ! Note BDS for PLANETIQ uses B1C as L1, and B2a as L2. f_L2=1176.450d6 ! ELSE ! default f_L1=1575.420d6 f_L2=1207.140d6 ENDIF ENDIF ! F_L1 and F_L2 based on GNSS ENDIF ! assign f_L1 and F_L2 ! 6.3 Set the bending angle method IF (INDEX(config%occ_method, "WO" ) == 1) THEN ro_data%bangle_method = "CT2 wave optics" ELSE IF (INDEX(config%occ_method, "FS" ) == 1) THEN ro_data%bangle_method = "FSI inversion" ELSE IF (INDEX(config%occ_method, "GO" ) == 1 ) THEN ro_data%bangle_method = "geometrical optics" ELSE CALL message(msg_warn, "Processing method " // config%occ_method // & " not supported") EXIT ENDIF !------------------------------------------------------------------------------- ! 7. Check coordinate frames !------------------------------------------------------------------------------- CALL ropp_pp_set_coordinates(ro_data) !------------------------------------------------------------------------------- ! 8. Determine occultation point georeference information !------------------------------------------------------------------------------- CALL occ_point( ro_data%lev1a%dtime, & ro_data%lev1a%r_leo, ro_data%lev1a%r_gns, & ro_data%georef%lat, ro_data%georef%lon, & ro_data%georef%r_coc, ro_data%georef%roc, & ro_data%georef%r_alt, ro_data%georef%viewang, & ro_data%georef%azimuth, & ro_data%georef%undulation, & ro_data%georef%time_offset, & config%egm96, config%corr_egm96) WRITE(outstr,'(A,1x,i2.2,a1,i2.2,a1,i4.4)') & "Occultation date (dd/mm/yyyy): = ", & ro_data%dtocc%day, "/", & ro_data%dtocc%month, "/", & ro_data%dtocc%year CALL message(msg_diag, outstr) WRITE(outstr,'(A,1x,i2.2,a1,i2.2,a1,i2.2)') & "Occultation time (hh:mm:ss): = ", & ro_data%dtocc%hour, ":", & ro_data%dtocc%minute, ":", & ro_data%dtocc%second CALL message(msg_diag, outstr) WRITE(outstr,'(A,1x,a1,f7.2,a2,f7.2,a2)') & "Occultation point (lat,lon): = ", & "(", ro_data%georef%lat, "N,", & ro_data%georef%lon, "E)" CALL message(msg_diag, outstr) WRITE(outstr,'(A,1x,e15.5)') & "Undulation (m) = ", ro_data%georef%undulation CALL message(msg_diag, outstr) !------------------------------------------------------------------------------- ! 9. Mission-specific pre-processing and input data cut-off !------------------------------------------------------------------------------- ! Do not perform if using output from WOPT IF (.NOT. wopt_mode) THEN CALL ropp_pp_preprocess(ro_data, config, diag) ENDIF !------------------------------------------------------------------------------- ! 10. Initialise bending angle structures !------------------------------------------------------------------------------- IF (config%obs_ok) THEN CALL ropp_io_init(smt_ba, ro_data%Lev1a%Npoints) CALL ropp_io_init(bangle, ro_data%Lev1a%Npoints) IF (ro_data%leo_id(1:2) == 'FY') & CALL ropp_io_init(gnos_ba, ro_data%Lev1a%Npoints) !------------------------------------------------------------------------------- ! 11. Geometric optics processing !------------------------------------------------------------------------------- CALL message(msg_info, & "Retrieving bending angle profiles by GEOMETRIC OPTICS \n") ! 11.1 Calculate smoothed bending angle DO j=1,ro_data%Lev1a%Npoints smt_ba%impact(j) = & impact_parameter(ro_data%Lev1a%r_leo(j,:)-ro_data%georef%r_coc, & ro_data%Lev1a%r_gns(j,:)-ro_data%georef%r_coc ) ENDDO Pmax = MAXVAL(smt_ba%impact) Pmin = MAX(MINVAL(smt_ba%impact), ro_data%georef%roc+2000.0_wp) ws_go_smooth = CEILING( config%fw_go_smooth * & (ro_data%Lev1a%Npoints-1)/ABS(Pmax-Pmin)) WRITE(pstr1, '(F10.3)') Pmin-ro_data%georef%roc WRITE(pstr2, '(F10.3)') Pmax-ro_data%georef%roc CALL message(msg_diag, "Smoothed bending angle. ") CALL message(msg_diag, "Pmin = " // pstr1 // " Pmax = " // pstr2) WRITE(nstr, '(i8)') ws_go_smooth CALL message(msg_diag, "ws_go_smooth " // nstr) IF (Pmax > Pmin) THEN CALL ropp_pp_bending_angle_go( ro_data%Lev1a%dtime, & ro_data%Lev1a%r_leo, & ro_data%Lev1a%r_gns, & ro_data%georef%r_coc, & ro_data%Lev1a%phase_L1, & ro_data%Lev1a%phase_L2, & ws_go_smooth, & config%filter_method, & smt_ba%impact_L1, & smt_ba%bangle_L1, & smt_ba%impact_L2, & smt_ba%bangle_L2 ) CALL ropp_pp_linear_combination( smt_ba%impact_L1, & smt_ba%bangle_L1, & smt_ba%impact_L2, & smt_ba%bangle_L2, & smt_ba%impact, & smt_ba%bangle ) ELSE config%obs_ok = .FALSE. CALL message(msg_warn, "Cannot process profile (smt): Pmin > Pmax") ENDIF ! 11.2 Calculate bending angle bangle%impact = smt_ba%impact_L1 Pmax = MAXVAL(smt_ba%impact_L1) Pmin = MAX(MINVAL(smt_ba%impact_L1), ro_data%georef%roc+2000.0_wp) ws_go_full = CEILING(config%fw_go_full*(ro_data%Lev1a%Npoints - 1)/ & ABS(Pmax - Pmin)) WRITE(pstr1, '(F10.3)') Pmin-ro_data%georef%roc WRITE(pstr2, '(F10.3)') Pmax-ro_data%georef%roc CALL message(msg_diag, "Full resolution bending angle. ") CALL message(msg_diag, "Pmin = " // pstr1 // " Pmax = " // pstr2) WRITE(nstr, '(i8)') ws_go_full CALL message(msg_diag, "ws_go_full " // nstr) IF (Pmax > Pmin) THEN CALL ropp_pp_bending_angle_go(ro_data%Lev1a%dtime, & ro_data%Lev1a%r_leo, & ro_data%Lev1a%r_gns, & ro_data%georef%r_coc, & ro_data%Lev1a%phase_L1, & ro_data%Lev1a%phase_L2, & ws_go_full, & config%filter_method, & bangle%impact_L1, & bangle%bangle_L1, & bangle%impact_L2, & bangle%bangle_L2 ) IF (ro_data%leo_id(1:2) == 'FY') THEN gnos_ba%impact_L1 = bangle%impact_L1 gnos_ba%impact_L2 = bangle%impact_L2 gnos_ba%bangle_L1 = bangle%bangle_L1 gnos_ba%bangle_L2 = bangle%bangle_L2 ENDIF WRITE(nstr, '(i10)') ro_data%Lev1a%Npoints WRITE(pstr2, '(f6.1)') Pmax/1000.0_wp WRITE(pstr1, '(f6.1)') Pmin/1000.0_wp CALL message(msg_info, & TRIM(nstr) // " data points in output between " // & TRIM(pstr1) // "km and " // TRIM(pstr2) // "km \n") ELSE config%obs_ok = .FALSE. CALL message(msg_warn, "Cannot process profile (full): Pmin > Pmax") ENDIF !------------------------------------------------------------------------------- ! 11b. FSI inversion processing !------------------------------------------------------------------------------- IF ( INDEX(config%occ_method, "FS" ) == 1 ) THEN WRITE(pstr1, '(F10.3)') Pmin-ro_data%georef%roc WRITE(pstr2, '(F10.3)') Pmax-ro_data%georef%roc CALL message(msg_diag, "Full resolution bending angle. ") CALL message(msg_diag, "Pmin = " // pstr1 // " Pmax = " // pstr2) ALLOCATE( go_l1_ip( SIZE( bangle%impact_L1 ) ) ) go_l1_ip = bangle%impact_L1 IF (Pmax > Pmin) THEN CALL ropp_pp_bending_angle_fsi(ro_data%Lev1a%dtime, & ro_data%Lev1a%r_leo, & ro_data%Lev1a%r_gns, & ro_data%georef%r_coc, & ro_data%georef%roc, & ro_data%Lev1a%phase_L1, & ro_data%Lev1a%phase_L2, & ro_data%Lev1a%snr_L1ca, & ro_data%Lev1a%snr_L2p, & ro_data%dtocc%month, & config%hmax_wo, & config%opt_DL2, & bangle%impact_L1, & bangle%bangle_L1, & bangle%bangle_L1_sigma, & bangle%impact_L2, & bangle%bangle_L2, & bangle%bangle_L2_sigma, & diag ) CALL message(msg_info, & "Retrieving bending angle profiles by FSI \n") WRITE(nstr, '(i10)') ro_data%Lev1a%Npoints WRITE(pstr2, '(f6.1)') Pmax/1000.0_wp WRITE(pstr1, '(f6.1)') Pmin/1000.0_wp CALL message(msg_info, & TRIM(nstr) // " data points in output between " // & TRIM(pstr1) // "km and " // TRIM(pstr2) // "km \n") IF (ro_data%leo_id(1:2) == 'FY') THEN ! Try to overwrite L2 with extrapolated stuff IF ( diag%L2_min_slta < ropp_MDTV ) THEN CALL message(msg_warn, 'No valid minimum L2 SLTA found ... ' // & 'cannot extrapolate L2') ELSE CALL ropp_pp_bending_angle_gnos(ro_data, config, & gnos_ba, bangle, diag) END IF CALL ropp_io_free(gnos_ba) ENDIF ELSE config%obs_ok = .FALSE. CALL message(msg_warn, "Cannot process profile (wo): Pmin > Pmax") ENDIF ENDIF !------------------------------------------------------------------------------- ! 12. Wave optics processing !------------------------------------------------------------------------------- IF ( INDEX(config%occ_method, "WO" ) == 1 ) THEN CALL message(msg_info, & "Retrieving bending angle profiles by WAVE OPTICS \n") ws_wo = CEILING(config%fw_wo * (bangle%npoints-1.0_wp) / & ABS(Pmax - Pmin)) ws_low = CEILING(config%fw_low * (bangle%npoints-1.0_wp) / & ABS(Pmax - Pmin)) WRITE(nstr,'(I5)') ws_go_smooth CALL message(msg_diag, 'WM = ' // nstr) WRITE(nstr,'(I5)') ws_wo CALL message(msg_diag, 'W = ' // nstr) WRITE(nstr,'(I5)') ws_low CALL message(msg_diag, 'WL = ' // nstr) ALLOCATE( go_l1_ip( SIZE( bangle%impact_L1 ) ) ) go_l1_ip = bangle%impact_L1 IF (Pmax > Pmin) THEN CALL ropp_pp_bending_angle_wo( ro_data%Lev1a%dtime, & ro_data%Lev1a%r_leo, & ro_data%Lev1a%r_gns, & ro_data%georef%r_coc, & ro_data%georef%roc, & ro_data%Lev1a%phase_L1, & ro_data%Lev1a%phase_L2, & ro_data%Lev1a%snr_L1ca, & ro_data%Lev1a%snr_L2p, & ws_go_smooth, & ws_wo, & ws_low, & config%hmax_wo, & config%filter_method, & config%opt_DL2, & config%cff, & config%dsh, & bangle%impact_L1, & bangle%bangle_L1, & bangle%bangle_L1_sigma, & bangle%impact_L2, & bangle%bangle_L2, & bangle%bangle_L2_sigma, & diag ) WRITE(nstr, '(i10)') ro_data%Lev1a%Npoints WRITE(pstr2, '(f6.1)') Pmax/1000.0_wp WRITE(pstr1, '(f6.1)') Pmin/1000.0_wp CALL message(msg_info, & TRIM(nstr) // " data points in output between " // & TRIM(pstr1) // "km and " // TRIM(pstr2) // "km \n") IF (ro_data%leo_id(1:2) == 'FY') THEN ! Try to overwrite L2 with extrapolated stuff IF ( diag%L2_min_slta < ropp_MDTV ) THEN CALL message(msg_warn, 'No valid minimum L2 SLTA found ... ' // & 'cannot extrapolate L2') ELSE CALL ropp_pp_bending_angle_gnos(ro_data, config, & gnos_ba, bangle, diag) END IF CALL ropp_io_free(gnos_ba) ENDIF ELSE config%obs_ok = .FALSE. CALL message(msg_warn, "Cannot process profile (wo): Pmin > Pmax") ENDIF ENDIF !------------------------------------------------------------------------------- ! 13. Check that profiles are monotonously increasing in height ! (index 1 towards surface) and valid bending angles !------------------------------------------------------------------------------- CALL ropp_pp_monotonous(bangle%impact_L1, -1) CALL ropp_pp_monotonous(bangle%impact_L2, -1) CALL ropp_pp_monotonous(smt_ba%impact, -1) IF ( ANY(bangle%bangle_L1 <= ropp_MDTV) .OR. & ANY(bangle%bangle_L2 <= ropp_MDTV) ) THEN config%obs_ok = .FALSE. CALL message(msg_warn, "Cannot process profile: " // & "Invalid bending angles computed \n") ENDIF ENDIF !------------------------------------------------------------------------------- ! 14. Compute un-optimised ionospheric corrected bending angle profile (LC) !------------------------------------------------------------------------------- IF (config%obs_ok) THEN CALL ropp_pp_linear_combination( bangle%impact_L1, bangle%bangle_L1, & bangle%impact_L2, bangle%bangle_L2, & bangle%impact, bangle%bangle ) ! Optionally calculate TP lat-lons to account for bending IF (config%tp_bending) THEN IF ( ANY(ro_data%georef%r_coc < 10.0_wp*ropp_MDTV) ) THEN WRITE(smess(1), '(ES10.3)') ro_data%georef%r_coc(1) WRITE(smess(2), '(ES10.3)') ro_data%georef%r_coc(2) WRITE(smess(3), '(ES10.3)') ro_data%georef%r_coc(3) CALL message(msg_fatal, 'Tangent points can only be computed with a ' // & 'valid centre of curvature. Given (x, y, z)_CoC = (' // & smess(1) // ', ' // smess(2) // ', ' // smess(3) // ')' ) ENDIF CALL tangent_point(ro_data%lev1a%r_leo, ro_data%lev1a%r_gns, & bangle%lat_tp, bangle%lon_tp, bangle%azimuth_tp, & ro_data%georef%r_coc, ip=bangle%impact) ELSE ! Assume straight line ray CALL tangent_point(ro_data%lev1a%r_leo, ro_data%lev1a%r_gns, & bangle%lat_tp, bangle%lon_tp, bangle%azimuth_tp, & ro_data%georef%r_coc) ENDIF !------------------------------------------------------------------------------- ! 15. Define configuration parameters !------------------------------------------------------------------------------- config%r_curve = ro_data%georef%roc config%npoints = bangle%npoints config%Pmax = MAXVAL(bangle%impact_L1(:)) config%Pmin = MINVAL(bangle%impact_L1(:)) !------------------------------------------------------------------------------- ! 16. Interpolate input data to standard grid !------------------------------------------------------------------------------- ! 16.1 Calculate size of standard grid IF ( INDEX(config%method, "NONE" ) == 1 ) THEN Pmax = MAXVAL(bangle%impact_L1(:)) - 5000.0_wp ELSE Pmax = config%ztop_invert + config%r_curve ENDIF Pmin = MINVAL(bangle%impact_L1(:)) nbi = 1 + CEILING((Pmax - Pmin)/config%dpi) WRITE(pstr2, '(f10.3)') Pmax-config%r_curve WRITE(pstr1, '(f10.3)') Pmin-config%r_curve CALL message(msg_diag, "Pmin = " // pstr1 // " Pmax = " // pstr2) WRITE(nstr, '(I6)') nbi CALL message(msg_diag, "Standard grid size nbi = " // nstr) ! 16.2 Initialise standard grid data structures for output CALL ropp_io_init(out_ba, nbi) CALL ropp_io_init(out_refrac, nbi) ! 16.3 Interpolate input data onto standard grid CALL ropp_pp_merge_profile(bangle%impact_L1, bangle%bangle_L1, & bangle%impact_L2, bangle%bangle_L2, & out_ba%impact_L1, out_ba%bangle_L1, & out_ba%impact_L2, out_ba%bangle_L2, & Pmin, Pmax) CALL ropp_pp_merge_profile(bangle%impact_L1,bangle%bangle_L1_sigma, & bangle%impact_L2,bangle%bangle_L2_sigma, & out_ba%impact_L1,out_ba%bangle_L1_sigma, & out_ba%impact_L2,out_ba%bangle_L2_sigma, & Pmin, Pmax) !------------------------------------------------------------------------------- ! 17. Ionospheric correction of bending angle profile by linear combination !------------------------------------------------------------------------------- IF ( INDEX(config%method, "NONE" ) == 1 ) THEN CALL message(msg_info, & "Correcting bending angle profile by LINEAR COMBINATION \n") ro_data%bangle_method = TRIM(ro_data%bangle_method) // "(LC)" CALL ropp_pp_linear_combination( out_ba%impact_L1, & out_ba%bangle_L1, & out_ba%impact_L2, & out_ba%bangle_L2, & out_ba%impact_opt, & out_ba%bangle_opt ) imax = SUM(MINLOC(ABS(out_ba%impact_L1(:)-(config%Pmax-5000.0_wp)))) imin = SUM(MINLOC(ABS(out_ba%impact_L1(:)-config%Pmin))) !------------------------------------------------------------------------------- ! 18. Ionospheric correct bending angle profile by statistical optimization !------------------------------------------------------------------------------- ELSE CALL message(msg_info, & "Correcting bending angle profile by STATISTICAL OPTIMISATION \n") ro_data%bangle_method = TRIM(ro_data%bangle_method) // " (StatOpt)" ! 18.1 Retrieve model bending angles on L1 impact parameter levels CALL ropp_io_init(mod_ba, nbi) IF (INDEX(config%method, "GMSIS" ) == 1 .OR. & INDEX(config%method, "GBARO" ) == 1 ) THEN IF (INDEX(config%sf_method, "serial" ) == 1) THEN CALL ropp_pp_search_model_refraction_new( config%mfile, & ro_data%dtocc%month, & ro_data%georef%lat, & ro_data%georef%lon, & smt_ba%impact, & smt_ba%bangle, & out_ba%impact_L1, & mod_ba%bangle, & config ) ELSE IF (INDEX(config%sf_method, "parallel" ) == 1) THEN CALL ropp_pp_search_model_refraction( config%mfile, & ro_data%dtocc%month, & ro_data%georef%lat, & ro_data%georef%lon, & smt_ba%impact, & smt_ba%bangle, & out_ba%impact_L1, & mod_ba%bangle, & config ) ELSE CALL message(msg_fatal, "Search and fit method " // & config%sf_method // " not supported") EXIT ENDIF ELSE IF (INDEX(config%method, "MSIS" ) == 1 .OR. & INDEX(config%method, "BARO" ) == 1) THEN CALL ropp_pp_model_refraction( config%mfile, & ro_data%dtocc%month, & ro_data%georef%lat, & ro_data%georef%lon, & out_ba%impact_L1, & mod_ba%bangle, & config ) ELSE IF (INDEX(config%method, "BG" ) == 1 ) THEN CALL ropp_pp_bg_refraction( config%bfile, & ro_data%dtocc%month, & ro_data%georef%lat, & ro_data%georef%lon, & out_ba%impact_L1, & mod_ba%bangle, & config ) ELSE config%obs_ok = .FALSE. CALL message(msg_warn, & "Statistical optimisation method " // config%method // & " not supported") EXIT ENDIF ! 18.2 Fit model bending angle with (smoothed) observed bending angles IF (INDEX(config%sf_method, "serial" ) == 1) THEN CALL ropp_pp_fit_model_refraction_new( smt_ba%impact, & smt_ba%bangle, & out_ba%impact_L1, & mod_ba%bangle, & config ) ELSE IF (INDEX(config%sf_method, "parallel" ) == 1) THEN CALL ropp_pp_fit_model_refraction( smt_ba%impact, & smt_ba%bangle, & out_ba%impact_L1, & mod_ba%bangle, & config ) ELSE CALL message(msg_fatal, "Search and fit method " // & config%sf_method // " not supported") EXIT ENDIF ! 18.3 Perform ionospheric correction with statistical optimization imax = SUM(MINLOC(ABS(out_ba%impact_L1(:)-(config%Pmax-5000.0_wp)))) imin = SUM(MINLOC(ABS(out_ba%impact_L1(:)-config%Pmin))) ! Merge observed bending angles with model above observation top WHERE ( out_ba%impact_L1 > out_ba%impact_L1(imax) ) out_ba%bangle_L1(:) = mod_ba%bangle(:) out_ba%bangle_L2(:) = mod_ba%bangle(:) out_ba%bangle_L1_sigma(:) = 0.0_wp out_ba%bangle_L2_sigma(:) = 0.0_wp END WHERE IF (imin < imax) THEN ALLOCATE(diag%ba_model(imax-imin+1)) diag%ba_model = mod_ba%bangle(imin:imax) out_ba%impact_opt(:) = out_ba%impact_L1(:) out_ba%bangle_opt(:) = out_ba%bangle_L1(:) CALL ropp_pp_ionospheric_correction( & out_ba%impact_L1(imin:imax), & ! L1 impact out_ba%bangle_L1(imin:imax), & ! L1 bending out_ba%impact_L2(imin:imax), & ! L2 impact out_ba%bangle_L2(imin:imax), & ! L2 bending out_ba%impact_L1(imin:imax), & ! Model impact mod_ba%bangle(imin:imax), & ! Model bending config, & ! Configuration out_ba%impact_opt(imin:imax), & ! Opt impact out_ba%bangle_opt(imin:imax), & ! Opt bending diag ) ! Diagnostics diag%err_neut = diag%err_neut + & out_ba%bangle_L1_sigma(imin:imax)**2 diag%sq = 100.0_wp * & MAXVAL(SQRT(diag%err_neut(:))/out_ba%bangle_opt(imin:imax)) ! Mean and standard difference between observational and climatology ! bending angles between 60-80 km ! Yong Chen on 08/09/2024 i2 = SUM(MINLOC(ABS(out_ba%impact_opt(:)-config%r_curve - ih2 * 1000.0_wp))) i1 = SUM(MINLOC(ABS(out_ba%impact_opt(:)-config%r_curve - ih1 * 1000.0_wp))) IF (i1 < i2) THEN ro_data%qa%smean = sum(out_ba%bangle_opt(i1:i2) - mod_ba%bangle(i1:i2))/REAL(i2 - i1 + 1, KIND=wp ) ro_data%qa%stdv = sqrt(sum((out_ba%bangle_opt(i1:i2) - mod_ba%bangle(i1:i2) & - ro_data%qa%smean)**2) /REAL(i2 - i1 + 1, KIND=wp ) ) ELSE ro_data%qa%smean = ropp_MDTV ro_data%qa%stdv = ropp_MDTV ENDIF ! Max and standard relative difference between observational and climatology ! bending angles between 25-40 km ! Yong Chen on 01/24/2025 i2 = SUM(MINLOC(ABS(out_ba%impact_opt(:)-config%r_curve - ih4 * 1000.0_wp))) i1 = SUM(MINLOC(ABS(out_ba%impact_opt(:)-config%r_curve - ih3 * 1000.0_wp))) IF (i1 < i2) THEN ALLOCATE(relative_ba(i2-i1+1)) relative_ba(1:i2-i1+1) = (out_ba%bangle_opt(i1:i2) - mod_ba%bangle(i1:i2))/mod_ba%bangle(i1:i2) ro_data%qa%reldevmax = maxval(abs(relative_ba)) ro_data%qa%reldevstd = sqrt(sum( (relative_ba(1:i2-i1+1) - & sum(relative_ba(1:i2-i1+1))/REAL(i2-i1+1, KIND=wp))**2) & /REAL(i2 - i1 + 1, KIND=wp ) ) DEALLOCATE(relative_ba) ELSE ro_data%qa%reldevmax = ropp_MDTV ro_data%qa%reldevstd = ropp_MDTV ENDIF ELSE config%obs_ok = .FALSE. CALL message(msg_warn, & "Cannot perform ionospheric correction: Imin >= Imax \n") ENDIF ENDIF ! 18.4 Interpolate LC bending angles, lat_tp, lon_tp, azimuth_tp to output grid out_ba%impact = out_ba%impact_opt CALL ropp_pp_interpol(bangle%impact, out_ba%impact, & bangle%bangle, out_ba%bangle) ! CALL ropp_pp_interpol(bangle%impact, out_ba%impact, & CALL ropp_pp_interpol(go_l1_ip, out_ba%impact, & bangle%lat_tp, out_ba%lat_tp) ! CALL ropp_pp_interpol(bangle%impact, out_ba%impact, & CALL ropp_pp_interpol(go_l1_ip, out_ba%impact, & bangle%lon_tp, out_ba%lon_tp) ! CALL ropp_pp_interpol(bangle%impact, out_ba%impact, & CALL ropp_pp_interpol(go_l1_ip, out_ba%impact, & bangle%azimuth_tp, out_ba%azimuth_tp) ! Eventually ! CALL ropp_pp_interpolate_latlonaz(go_l1_ip, & ! out_ba%impact, & ! bangle%lat_tp, & ! bangle%lon_tp, & ! out_ba%lat_tp, & ! out_ba%lon_tp, & ! bangle%azimuth_tp, & ! out_ba%azimuth_tp) WRITE(nstr, '(i10)') imax-imin+1 WRITE(pstr2, '(f6.1)') out_ba%impact(imax)/1000.0_wp WRITE(pstr1, '(f6.1)') out_ba%impact(imin)/1000.0_wp CALL message(msg_info, & TRIM(nstr) // " data points in output between " // & TRIM(pstr1) // "km and " // TRIM(pstr2) // "km \n") ENDIF !------------------------------------------------------------------------------- ! 19. Perform inverse Abel transform to compute refractivity !------------------------------------------------------------------------------- IF (config%obs_ok .AND. config%output_lev2a) THEN ro_data%refrac_method = "Abel transform (" // TRIM(config%abel) // ")" ! 19.1 Abel inversion of corrected bending angle profile IF ( INDEX(config%method, "NONE" ) == 1 ) THEN CALL message(msg_info, & "Retrieving refractivity profile by " // TRIM(config%abel) // & " ABEL TRANSFORM \n") IF ( INDEX(config%abel, "LIN" ) == 1 ) THEN CALL ropp_pp_invert_LIN(out_ba%impact_opt, out_ba%bangle_opt, & out_ba%impact_opt, out_refrac%refrac) ELSE IF ( INDEX(config%abel, "EXP" ) == 1 ) THEN CALL ropp_pp_invert_EXP(out_ba%impact_opt, out_ba%bangle_opt, & out_ba%impact_opt, out_refrac%refrac) ELSE config%obs_ok = .FALSE. CALL message(msg_warn, & "Abel integral method " // config%abel // " not supported") ENDIF ! 19.2 Abel inversion of corrected bending angle profile with stat opt ELSE CALL message(msg_info, & "Retrieving refractivity profile by " // TRIM(config%abel) // & " ABEL TRANSFORM with STAT OPT \n") CALL ropp_pp_invert_refraction( config%mfile, & ro_data%dtocc%month, & ro_data%georef%lat, & ro_data%georef%lon, & out_ba%impact_opt, & out_ba%bangle_opt, & out_refrac%geop_refrac, & out_refrac%refrac, & config ) ENDIF !------------------------------------------------------------------------------- ! 20. Compute output height scales !------------------------------------------------------------------------------- ! Override default output height scales wrt geoid if requested IF (earth_ellipsoid) ro_data%georef%undulation = ropp_MDFV IF (ro_data%georef%undulation > ropp_MDTV) THEN CALL message(msg_info, "Writing output altitude scales " // & "with respect to EGM96 GEOID") out_refrac%alt_refrac = & ((out_ba%impact_opt / (1.0_wp + out_refrac%refrac * 1.e-6_wp)) & - (ro_data%GEOref%roc + ro_data%GEOref%undulation)) ELSE IF (.NOT. earth_ellipsoid) & CALL message(msg_warn, "Invalid undulation calculated. " // & "Check for valid EGM geoid coefficient and correction file.") CALL message(msg_info, "Writing output altitude scales " // & "with respect to WGS84 ELLIPSOID.") out_refrac%alt_refrac = & ((out_ba%impact_opt / (1.0_wp + out_refrac%refrac * 1.e-6_wp)) & - ro_data%GEOref%roc) ENDIF ! Yong Chen 03/02/2023, check the lowest altitude, and make sure it is above ZERO ialt = SUM(MINLOC(out_refrac%alt_refrac(:), out_refrac%alt_refrac(:) > 0.0_wp)) IF (ialt > 1 ) THEN ALLOCATE(alt_tmp( SIZE( out_refrac%alt_refrac ) ) ) ALLOCATE(refrac_tmp( SIZE( out_refrac%refrac ) ) ) alt_tmp=out_refrac%alt_refrac alt_tmp(1:ialt-1) = 0.0_wp CALL ropp_pp_interpol(out_refrac%alt_refrac, alt_tmp, & out_refrac%refrac, refrac_tmp) out_refrac%refrac(1:ialt-1) = refrac_tmp(1:ialt-1) out_refrac%alt_refrac(1:ialt-1) = 0.0_wp DEALLOCATE(alt_tmp) DEALLOCATE(refrac_tmp) ENDIF out_refrac%geop_refrac = & geometric2geopotential(ro_data%georef%lat, out_refrac%alt_refrac) !------------------------------------------------------------------------------- ! 21. Compute Tdry and Pdry - though only Tdry is written to ro_data !------------------------------------------------------------------------------- CALL message(msg_info, "Computing dry temperature \n") CALL ropp_io_init(dum_meteo, nbi) ! easy way to get dummy storage dum_meteo%shum = 0.0_wp CALL ropp_pp_tdry(ro_data%georef%lat, out_refrac%alt_refrac, & out_refrac%refrac, dum_meteo%shum, out_refrac%dry_temp, & dum_meteo%press) CALL ropp_io_free(dum_meteo) ENDIF ! (config%obs_ok .AND. config%output_lev2a) !------------------------------------------------------------------------------- ! 22. Copy retrieved profiles to RO structure ! - only output data within observed range !------------------------------------------------------------------------------- IF (config%obs_ok) THEN IF (config%output_lev1b) THEN IF (.NOT. output_full) out_ba%npoints = imax - imin + 1 CALL ropp_io_roprof2roprof(out_ba, ro_data%lev1b) ELSE CALL ropp_io_free(ro_data%lev1b) ENDIF IF (config%output_lev2a) THEN IF (.NOT. output_full) out_refrac%npoints = imax - imin + 1 CALL ropp_io_roprof2roprof(out_refrac, ro_data%lev2a) ELSE CALL ropp_io_free(ro_data%lev2a) ENDIF ENDIF IF (.NOT. config%output_lev1a) CALL ropp_io_free(ro_data%lev1a) CALL ropp_io_free(ro_data%vlist) !! interim, avoid writing lcf data IF (config%output_diag) THEN CALL message(msg_info, "Writing additional diagnostic output \n") CALL ropp_pp_diag2roprof(diag, ro_data) ENDIF !------------------------------------------------------------------------------- ! 23. Write data !------------------------------------------------------------------------------- ! Yong Chen on 01/25/2025 ! setup the final quality flag ro_data%qa%finalQA = 0 IF (ro_data%qa%reldevmax > th_reldevmax .or. & ro_data%qa%reldevstd > th_reldevstd .or. & ro_data%qa%difmaxref > th_difmaxref .or. & ro_data%qa%snr1avg < th_snr1avg .or. & ro_data%qa%maxdifphase > th_maxdifphase .or. & ro_data%qa%stdv > th_stdv .or. & abs(ro_data%qa%smean) > th_smean ) THEN ro_data%qa%finalQA = 1 ENDIF CALL message(msg_info, "Writing to output file " // TRIM(ofile) // "\n") CALL ropp_io_write(ro_data, ofile, append=.TRUE., ranchk=ranchk ) !------------------------------------------------------------------------------- ! 24. Clean up !------------------------------------------------------------------------------- CALL ropp_io_free(ro_data) CALL ropp_io_free(smt_ba) CALL ropp_io_free(bangle) CALL ropp_io_free(out_ba) CALL ropp_io_free(out_refrac) CALL ropp_io_free(mod_ba) CALL ropp_io_free(dum_meteo) IF (ro_data%leo_id(1:2) == 'FY') CALL ropp_io_free(gnos_ba) IF (ASSOCIATED(diag%ba_ion)) DEALLOCATE(diag%ba_ion) IF (ASSOCIATED(diag%err_ion)) DEALLOCATE(diag%err_ion) IF (ASSOCIATED(diag%err_neut)) DEALLOCATE(diag%err_neut) IF (ASSOCIATED(diag%wt_data)) DEALLOCATE(diag%wt_data) IF (ASSOCIATED(diag%ba_model)) DEALLOCATE(diag%ba_model) IF (ASSOCIATED(diag%CTimpact)) DEALLOCATE(diag%CTimpact) IF (ASSOCIATED(diag%CTamplitude)) DEALLOCATE(diag%CTamplitude) IF (ASSOCIATED(diag%CTamplitude_smt)) DEALLOCATE(diag%CTamplitude_smt) IF (ASSOCIATED(diag%CTimpactL2)) DEALLOCATE(diag%CTimpactL2) IF (ASSOCIATED(diag%CTamplitudeL2)) DEALLOCATE(diag%CTamplitudeL2) IF (ASSOCIATED(diag%CTamplitudeL2_smt)) DEALLOCATE(diag%CTamplitudeL2_smt) IF (ALLOCATED(go_l1_ip)) DEALLOCATE(go_l1_ip) END DO END DO CALL EXIT(msg_exit_status) CONTAINS !------------------------------------------------------------------------------- ! 25. Usage information !------------------------------------------------------------------------------- SUBROUTINE usage() PRINT *, 'Purpose:' PRINT *, ' Calculate corrected bending angle, refractivity and ' PRINT *, ' dry temperature from L1 and L2 phase data' PRINT *, 'Usage:' PRINT *, ' > ropp_pp_occ_tool [] ' PRINT *, 'Options' PRINT *, ' -c name of configuration file' PRINT *, ' -o name of ROPP netCDF output file' PRINT *, ' -m ionospheric correction method' PRINT *, ' [NONE,MSIS,GMSIS,BG,BARO,GBARO], (default GMSIS)' PRINT *, ' -mfile model coefficients file' PRINT *, ' [MSIS_coeff.nc,BAROCLIM_coeff.nc], (default MSIS_coeff.nc)' PRINT *, ' -bfile background model atmospheric profile file path' PRINT *, ' (if using BG ionospheric correction method)' PRINT *, ' -navfile external navigation bit *_txt file path' PRINT *, ' (default: use internal correction)' PRINT *, ' -occ processing method, WO or GO (default: WO)' PRINT *, ' -filter filtering method, slpoly or optest' PRINT *, ' (default: slpoly = sliding polynomial)' PRINT *, ' -fit apply 2-parameter regression fit to model' PRINT *, ' -ellipsoid output height with respect to WGS84 ellipsoid' PRINT *, ' (default: output wrt EGM96 geoid)' PRINT *, ' -full output whole profile' PRINT *, ' (default: only output data within the observed range)' PRINT *, ' -w process output from ropp_pp_wopt_tool' PRINT *, ' (default: use other types of input data)' PRINT *, ' -d output additional diagnostics' PRINT *, ' -h this help' PRINT *, ' -v version information' PRINT *, '' END SUBROUTINE usage !------------------------------------------------------------------------------- ! 26. Version information !------------------------------------------------------------------------------- SUBROUTINE version_info() CHARACTER (LEN=40) :: version version = ropp_pp_version() PRINT *, 'ropp_pp_occ_tool - Pre-processor occ tool: ' PRINT *, ' Calculate corrected bending angles, ' PRINT *, ' refractivity and dry temperature ' PRINT *, ' from excess phase' PRINT *, '' PRINT *, 'This program is part of ROPP (PP) Release ' // TRIM(version) PRINT *, '' END SUBROUTINE version_info END PROGRAM ropp_pp_occ_tool