! $Id: ropp_pp_wopt_2d_tool.f90 4887 2016-05-25 15:48:28Z sti $ PROGRAM ropp_pp_wopt_2D_tool !****p* Programs/ropp_pp_wopt_2d_tool * ! ! NAME ! ropp_pp_wopt_2d_tool ! ! SYNOPSIS ! 2D wave optics propagator tool to simulate excess phases generated by ! a given 2D slice refractivity field. ! ! > ropp_pp_wopt_2d_tool [-o ] ! [-c ] ! [-f] ! [-d] [-h] [-v] ! ! ARGUMENTS ! input file name ! ! OPTIONS ! -o name of output file ! (default: ropp_pp_wopt_2d_tool_out.nc) ! -c name of configuration namelist ! -f estimate bending angles from the calculated excess phases ! using a Full Spectrum Inversion (FSI) method ! -d write out (lots of!) extra diagnostics ! -h help ! -v version information ! ! DESCRIPTION ! This program reads 2D RO refractivity data from a standard ROPP netCDF file, ! and, optionally, the 'configuration parameters' from a namelist. The output ! is an ROPP file including level 1a excess phase and time data. Note that, ! currently, the phases and the amplitudes are the same for both RO ! frequencies. ! ! NOTES ! Assumes a 2D refractivity slice. ! ! 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 level 1a excess phase fields for the refractivity slice ! in the data directory: ! ! > ropp_pp_wopt_2d_tool ../data/CASE_12_ref2d.nc ! ! SEE ALSO ! For example IDL plotting procedures to view the resulting fields see ! ropp_pp/tests/it_pp_wopt_0{1,2,3}.pro ! ! REFERENCES ! Modeling and inverting radio occultation signals in the moist troposphere, ! Sergey V. Sokolovskiy, 2001, Radio Sci 36 (3), pp 441-458, May/June 2001. ! ROM SAF Report 28, SAF/ROM/ECMWF/REP/RSR/028 ! ! AUTHOR ! ECMWF, Reading, UK, and 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 messages USE ropp_utils USE ropp_io USE ropp_io_types, ONLY: ROprof ! USE ropp_pp, ONLY: ropp_pp_version USE ropp_pp, ONLY: ropp_pp_version USE ropp_pp_spline, ONLY: ropp_pp_init_spline USE ropp_pp_constants, ONLY: f_L1, c_light, pi1 => pi USE ropp_pp_wopt IMPLICIT NONE ! Specifications INTEGER, PARAMETER :: nx_default = 401 ! Number of horizontal screens INTEGER :: nx = nx_default INTEGER, PARAMETER :: log2ny_default = 19 ! Number of vertical points in each screen INTEGER :: log2ny = log2ny_default INTEGER :: ny ! =2**log2ny later INTEGER, PARAMETER :: n_leo_default = 20000 ! Number of points (and times) ! for which phases will be calculated INTEGER :: n_leo = n_leo_default INTEGER, PARAMETER :: nsample_default = 100 ! Number of points in each 'sample' of the final screen, ! used when calculating signal at LEO INTEGER :: nsample = nsample_default REAL(wp), PARAMETER :: dx_default = 5.0e3_wp ! Separation of screens REAL(wp) :: dx = dx_default REAL(wp), PARAMETER :: dy_default = 1.0_wp ! Separation of vertical points REAL(wp) :: dy = dy_default REAL(wp), PARAMETER :: ymin_default = -3.0e5_wp ! Lower limit of screen REAL(wp) :: ymin = ymin_default REAL(wp), PARAMETER :: y_apodize_default = 1.2e5_wp ! Height where apodization starts REAL(wp) :: y_apodize = y_apodize_default REAL(wp), PARAMETER :: leo_altitude_default = 8.0e5_wp ! LEO altitude REAL(wp) :: leo_altitude = leo_altitude_default REAL(wp), PARAMETER :: gps_altitude_default = 2.02e7_wp ! GNSS altitude REAL(wp) :: gps_altitude = gps_altitude_default REAL(wp), PARAMETER :: tpt_altitude_default = 8.0e4_wp ! Nominal TPT altitude REAL(wp) :: tpt_altitude = tpt_altitude_default REAL(wp), PARAMETER :: delta_t_default = 0.005_wp ! Time between phase/amp samples REAL(wp) :: delta_t = delta_t_default ! Local COMPLEX(wp), ALLOCATABLE, DIMENSION(:) :: U_init, U_final, U_leo ! Complex amplitudes at first screen, ! last screen and at LEO. COMPLEX(wp), PARAMETER :: ci = (0.0_wp, 1.0_wp) REAL(wp), ALLOCATABLE, DIMENSION(:) :: xpos, ypos, s_geom REAL(wp), ALLOCATABLE, DIMENSION(:) :: time, theta REAL(wp), ALLOCATABLE, DIMENSION(:) :: amp_leo, phase_leo ! Amplitude and accumulated phase at LEO REAL(wp), ALLOCATABLE, DIMENSION(:) :: alpha_fsi, impact_fsi, amp_fsi ! results of FSI calculation REAL(wp), ALLOCATABLE, DIMENSION(:, :) :: zed, log_refrac, d2logndy2 REAL(wp), ALLOCATABLE, DIMENSION(:, :) :: r_leo, r_gps REAL(wp), PARAMETER :: radius_earth = Re ! = Earth radius ! Re = 6378137.0_wp in ropp_utils/coordinates.f90 REAL(wp) :: radius_leo, radius_gps, radius_tpt REAL(wp) :: rad ! Distance from GNSS to point on screen REAL(wp), PARAMETER :: v_leo_nom = 7.4e3_wp ! Nominal LEO velocity REAL(wp) :: theta_min ! LEO angle at (nominal) TPT REAL(wp) :: theta_dot ! LEO angular velocity REAL(wp) :: kval ! 2 pi f / c REAL(wp) :: roc_plus_und ! RoC plus undulation REAL(wp) :: twopi1 ! 2 pi REAL(wp), PARAMETER :: min_impact_height=2.0e3_wp ! Nullify FSI output below this REAL(wp), PARAMETER :: amp_threshold=0.2_wp ! Nullify FSI output below this REAL(wp) :: dtheta ! Angular separation (rads) between input profiles INTEGER :: nlev, nhoriz ! Number of refractivity levels and locations in plane INTEGER :: namelist_lun INTEGER :: iostatus INTEGER :: i, iargc, argc, idummy, n_files INTEGER, DIMENSION(:), POINTER :: idx => NULL() ! array indices INTEGER :: nidx, imid LOGICAL :: give_help = .FALSE. LOGICAL :: cfile_exists = .FALSE. LOGICAL :: run_FSI = .FALSE. CHARACTER(len = 4096) :: ifile='' ! input file CHARACTER(len = 4096) :: cfile='' ! config file CHARACTER(len = 4096) :: ofile='ropp_pp_wopt_2d_tool_out.nc' ! output file CHARACTER(len = 256) :: buffer CHARACTER(len = 8) :: snx CHARACTER(len = 3) :: snlev, snhoriz CHARACTER(len = 10) :: sdtheta CHARACTER(len = 12) :: sroc, sund TYPE(ROprof2d) :: ro_data1 ! Input RO data TYPE(ROprof) :: ro_data ! Output RO data REAL(wp), DIMENSION(:,:), ALLOCATABLE :: xleo ! LEO position by regression REAL(wp), DIMENSION(:,:), ALLOCATABLE :: vleo ! LEO velocity by regression REAL(wp), DIMENSION(:,:), ALLOCATABLE :: xgns ! GPS position by regression REAL(wp), DIMENSION(:,:), ALLOCATABLE :: vgns ! GPS velocity by regression ! Configuration namelist NAMELIST / config_wopt / nx, log2ny, n_leo, nsample, & dx, dy, ymin, y_apodize, & leo_altitude, gps_altitude, tpt_altitude, & delta_t !------------------------------------------------------------------------------- ! 2. Default settings !------------------------------------------------------------------------------- CALL message_set_routine ( 'ropp_pp_wopt_2d_tool' ) CALL message(msg_noin, '') CALL message(msg_noin, & '-----------------------------------------------------------------------') CALL message(msg_noin, & ' ROPP 2D Wave Optics Propagator' ) CALL message(msg_noin, & '-----------------------------------------------------------------------') CALL message(msg_noin, '') !------------------------------------------------------------------------------- ! 3. Command line arguments !------------------------------------------------------------------------------- argc = iargc() i = 1 n_files = 0 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 (namelist) CALL getarg(i+1, buffer) cfile = buffer i = i + 1 CASE('-f', '-F') ! Run FSI on output run_fsi = .TRUE. CASE('-d') ! Additional diagnostic mode msg_MODE = VerboseMode CASE('-h', '--help', '?') ! Give help give_help = .TRUE. CASE('-v', '-V', '--version') ! Version information CALL version_info() CALL EXIT(msg_exit_ok) CASE default ! Input file name IF ( buffer(1:1) /= '-' ) THEN ifile = buffer n_files = n_files + 1 END IF END SELECT i = i + 1 ENDDO IF ( n_files == 0 .AND. .NOT. give_help ) & CALL message ( msg_error, "No input file(s) specified" ) IF (argc == 0 .OR. n_files == 0 .OR. give_help) THEN CALL usage() CALL EXIT(msg_exit_status) ENDIF !------------------------------------------------------------------------------- ! 4. Read configuration file, if given and possible !------------------------------------------------------------------------------- IF ( TRIM(ADJUSTL(cfile)) /= '' ) THEN INQUIRE ( FILE=cfile, EXIST=cfile_exists ) IF ( cfile_exists ) THEN CALL message ( msg_info, 'Reading config namelist from ' // & TRIM(ADJUSTL(cfile)) ) namelist_lun = Get_IO_Unit() ! From ropp_utils OPEN ( UNIT=namelist_lun, FILE=cfile, STATUS='OLD', ACTION='READ', IOSTAT=iostatus ) IF ( iostatus > 0 ) & CALL message ( msg_fatal, 'I/O error while opening config file ' // & TRIM(ADJUSTL(cfile)) ) READ ( UNIT=namelist_lun, NML=config_wopt, IOSTAT=iostatus ) IF ( iostatus > 0 ) & CALL message ( msg_fatal, 'I/O error while reading config namelist ' // & TRIM(ADJUSTL(cfile)) ) IF ( msg_MODE == VerboseMode ) THEN CALL message( msg_diag, 'Contents of config_wopt namelist:' ) WRITE ( *, NML=config_wopt ) ENDIF CLOSE ( UNIT=namelist_lun ) ELSE CALL message( msg_fatal, 'Requested configuration file ' // & TRIM(ADJUSTL(cfile)) // ' does not exist' ) END IF END IF !------------------------------------------------------------------------------- ! 5. Remove pre-existing output file !------------------------------------------------------------------------------- CALL file_delete(ofile, idummy) !------------------------------------------------------------------------------- ! 6. Define the coordinates and velocities !------------------------------------------------------------------------------- ! 6.1 Basic parameters ! -------------------- IF ( (nx/2)*2 == nx ) THEN WRITE ( snx, FMT='(I8)' ) nx CALL message(msg_fatal, 'nx (= ' // snx // ') must be an odd number') END IF ny = 2**log2ny radius_leo = radius_earth + leo_altitude radius_gps = radius_earth + gps_altitude radius_tpt = radius_earth + tpt_altitude theta_dot = v_leo_nom / radius_leo theta_min = ACOS( radius_tpt / radius_leo ) twopi1 = 2.0_wp * pi1 kval = twopi1 * f_l1 / c_light ! L1 wave number ! 6.2 Set ygrid ! ------------- ALLOCATE (ypos(ny)) ypos = ymin + dy * (/ (i, i=0, ny-1) /) ! 6.3 Set xgrid ! ------------- ALLOCATE (xpos(nx)) xpos = dx * (/ (i, i=-(nx-1)/2, (nx-1)/2) /) ! 6.4 Define the (stationary) GNSS position ! ----------------------------------------- ALLOCATE (r_gps(n_leo, 3)) r_gps(:, 1) = radius_tpt r_gps(:, 2) = -SQRT(radius_gps**2 - radius_tpt**2) r_gps(:, 3) = 0.0_wp ! equatorial orbit ! 6.5 Define the (moving) LEO position ! ------------------------------------ ALLOCATE (time(n_leo), theta(n_leo), s_geom(n_leo)) ALLOCATE (r_leo(n_leo, 3)) time = delta_t * (/ (i, i=0, n_leo-1) /) theta = theta_min + theta_dot * time r_leo(:, 1) = radius_leo * COS(theta) r_leo(:, 2) = radius_leo * SIN(theta) r_leo(:, 3) = 0.0_wp ! equatorial orbit ! For later input into CT2. Used in WOP code to calculate RoC and undulation. DO i = 1, n_leo s_geom(i) = SQRT ( SUM ( (r_leo(i, :) - r_gps(i, :))**2 ) ) END DO ! 6.6 Compute radius of curvature and undulation for orbits ! --------------------------------------------------------- CALL ropp_io_init(ro_data%Lev1a, n_leo) CALL message( msg_info, 'Calculating occultation point' ) ! Array allocation ALLOCATE(xleo(n_leo,3)) ALLOCATE(xgns(n_leo,3)) ALLOCATE(vleo(n_leo,3)) ALLOCATE(vgns(n_leo,3)) CALL ropp_pp_satellite_velocities(time, r_leo, r_gps, & xleo, vleo, xgns, vgns) CALL occ_point(time, r_leo, r_gps, vleo, & ! in ro_data%georef%lat, ro_data%georef%lon, & ! out ro_data%georef%r_coc, ro_data%georef%roc, & ! out ro_data%georef%r_alt, ro_data%georef%viewang, & ! out ro_data%georef%azimuth, & ! out ro_data%georef%undulation, & ! out ro_data%georef%time_offset, & ! out ro_data%georef%slta_hgt_max, & ro_data%georef%slta_hgt_min) DEALLOCATE(xleo) DEALLOCATE(xgns) DEALLOCATE(vleo) DEALLOCATE(vgns) IF ( ro_data%georef%roc < ropp_MDTV ) THEN CALL message( msg_fatal, 'Radius of curvature is incalculable ... ' // & 'cannot continue' ) ELSE WRITE (sroc, '(F12.3)') ro_data%georef%roc / 1.0E3_wp CALL message( msg_info, 'Radius of curvature = ' // TRIM(ADJUSTL(sroc)) // ' km.' ) END IF IF ( ro_data%georef%undulation < ropp_MDTV ) THEN CALL message( msg_warn, 'Undulation is incalculable ... ' // & 'continuing with undulation set to zero' ) CALL message( msg_cont, '(Try setting the environment variables ' // & '$GEOPOT_COEF and $GEOPOT_CORR explicitly, e.g. ' ) CALL message( msg_cont, 'GEOPOT_COEF=$ROPP_SRC/ropp_pp/data/egm96.dat and ' // & 'GEOPOT_CORR=$ROPP_SRC/ropp_pp/data/corrcoef.dat).' ) ro_data%georef%undulation = 0.0_wp ELSE WRITE (sund, '(F12.3)') ro_data%georef%undulation CALL message( msg_info, 'Undulation = ' // TRIM(ADJUSTL(sund)) // ' m.' ) END IF ! Radius of curvature plus undulation roc_plus_und = ro_data%georef%roc + ro_data%georef%undulation !------------------------------------------------------------------------------- ! 7. Read the 2D refractivity slice !------------------------------------------------------------------------------- CALL message( msg_info, 'Reading refractivity slice from ' // & TRIM(ADJUSTL(ifile)) ) CALL ropp_io_read(ro_data1, ifile) nlev = ro_data1 % lev2a % npoints nhoriz = ro_data1 % lev2a % nhoriz dtheta = ro_data1 % lev2c % dtheta IF ( nhoriz == 1 ) THEN CALL message( msg_fatal, 'Input refractivity field is one-dimensional ... ' // & 'cannot use with 2D wave optics propagator. \n' // & 'Try again with a 2D field, or use the 1D propagator.' ) ENDIF WRITE (snlev, '(I3)') nlev WRITE (snhoriz, '(I3)') nhoriz WRITE (sdtheta, '(F10.3)') dtheta CALL message(msg_diag, 'Input file parameters: ' // & 'nlev = ' // snlev // '; ' // & 'nhoriz = ' // snhoriz // '; ' // & 'dtheta = ' // sdtheta // ' rad.') ALLOCATE (zed(nlev, nhoriz), log_refrac(nlev, nhoriz)) zed = ro_data1 % lev2a % alt_refrac log_refrac = LOG( MAX(ropp_ZDTV, ro_data1 % Lev2a % refrac) ) ! What about missing data? ! Append any input level 1b data to the output ROprof structure ro_data % lev1b = ro_data1 % lev1b ! Append any input level 2a in the centre profile to the output ROprof structure imid = (nhoriz + 1) / 2 CALL ropp_io_init(ro_data%Lev2a, ro_data1 % lev2a % npoints) ro_data % lev2a % missing = ro_data1 % lev2a % missing ro_data % lev2a % alt_refrac(:) = ro_data1 % lev2a % alt_refrac(:, imid) ro_data % lev2a % geop_refrac(:) = ro_data1 % lev2a % geop_refrac(:, imid) ro_data % lev2a % refrac(:) = ro_data1 % lev2a % refrac(:, imid) ro_data % lev2a % refrac_sigma(:) = ro_data1 % lev2a % refrac_sigma(:, imid) ro_data % lev2a % refrac_qual(:) = ro_data1 % lev2a % refrac_qual(:, imid) ro_data % lev2a % dry_temp(:) = ro_data1 % lev2a % dry_temp(:, imid) ro_data % lev2a % dry_temp_sigma(:) = ro_data1 % lev2a % dry_temp_sigma(:, imid) ro_data % lev2a % dry_temp_qual(:) = ro_data1 % lev2a % dry_temp_qual(:, imid) CALL ropp_io_free(ro_data1) !------------------------------------------------------------------------------- ! 8. Initialise the complex amplitude on the first screen by assuming ! a spherical wave ~ Amp_source*exp(i*k*r)/r coming from the _stationary_ GNSS. !------------------------------------------------------------------------------- CALL message( msg_info, 'Calculating signal on first screen' ) ! 8.1 Calculate U_init (normalised such that |U_init| = 1 where ypos+roc+und=y_GNSS) ! ---------------------------------------------------------------------------------- ALLOCATE (U_init(ny)) DO i = 1, ny rad = SQRT( (r_gps(1, 2) - xpos(1))**2 + & ! distance from LEO to ith position on first screen (r_gps(1, 1) - (ypos(i) + roc_plus_und))**2 ) U_init(i) = ( ABS(r_gps(1, 2) - xpos(1)) / rad ) * EXP( ci * MODULO(kval*rad, twopi1) ) ENDDO ! 8.2 Add U_init to ro_data ! ------------------------- IF ( msg_mode == VerboseMode ) THEN CALL ropp_io_addvar_rodataD1d( ro_data, & name = "Amp_U_init", & long_name = "Amplitude of complex signal on first screen", & units = "", & range = (/ -1.0E7_wp, 1.0E7_wp /), & DATA = ABS(U_init) ) CALL ropp_io_addvar_rodataD1d( ro_data, & name = "Phase_U_init", & long_name = "Phase of complex signal on first screen", & units = "m", & range = (/ -1.0E7_wp, 1.0E7_wp /), & DATA = cumulative_phase(U_init) / kval ) END IF !------------------------------------------------------------------------------- ! 9. Apodize signal on first screen !------------------------------------------------------------------------------- CALL message( msg_info, 'Apodizing signal on first screen' ) CALL ropp_pp_wopt_window(xpos(1), ypos, radius_earth, y_apodize, U_init) !------------------------------------------------------------------------------- ! 10. Initialise spline fitting of log_refrac. d2logndy2 are the spline coefficients. !------------------------------------------------------------------------------- CALL message( msg_info, 'Fitting splines to each refractivity profile' ) ALLOCATE (d2logndy2(nlev, nhoriz)) DO i = 1, nhoriz ! each position in slice CALL ropp_pp_init_spline(zed(:,i), log_refrac(:,i), d2logndy2(:,i)) ENDDO !------------------------------------------------------------------------------- ! 11. Propagate signal through the screens !------------------------------------------------------------------------------- CALL message( msg_info, 'Propagating signal through the screens' ) ALLOCATE (U_final(ny)) CALL ropp_pp_wopt_mps_2d(xpos, & ypos, & roc_plus_und, & ! radius of curvature + plus undulation y_apodize, & kval, & dtheta, & ! angular separation zed, & ! 2D refrac heights log_refrac, & ! 2D log(refrac) d2logndy2, & ! 2D spline coefficients U_init, & ! 2D phase/amp at first screen U_final) ! 2D phase/amp at final screen !------------------------------------------------------------------------------- ! 12. Propagate signal from final screen to LEO orbit !------------------------------------------------------------------------------- CALL message( msg_info, 'Calculating phase and amplitude at LEO' ) ALLOCATE (U_leo(n_leo), phase_leo(n_leo), amp_leo(n_leo)) CALL ropp_pp_wopt_propagate_to_leo(nsample, & U_final, & kval, & ypos, & dy, & xpos(nx), & ! position of final screen r_leo(:, 2), & r_leo(:, 1)-roc_plus_und, & s_geom, & ! |r_GNSS - r_LEO| U_leo, & ! field at LEO phase_leo, & ! accumulated phase at LEO (m) amp_leo) ! amplitude at LEO !------------------------------------------------------------------------------- ! 13. If desired, call FSI to estimate bending angle from the calculated excess phases !------------------------------------------------------------------------------- IF ( run_FSI ) THEN CALL message( msg_info, 'Calculating bending angle with FSI' ) ALLOCATE (impact_fsi(ny), alpha_fsi(ny), amp_fsi(ny)) CALL ropp_pp_wopt_fsi_quick(kval, theta_dot, time, s_geom, & phase_leo, amp_leo, & r_leo(:, 2), r_leo(:, 1)-roc_plus_und+radius_earth, & r_gps(:, 2), r_gps(:, 1), & impact_fsi, alpha_fsi, amp_fsi) ! Add FSI output as extra data to the RO profile ! ---------------------------------------------- impact_fsi = impact_fsi - radius_earth idx => WHERE ( ( (impact_fsi > min_impact_height ) .AND. & (amp_fsi > amp_threshold*MAXVAL(amp_fsi)) ), nidx ) IF ( nidx > 0 ) THEN CALL ropp_io_addvar_rodataD1d( ro_data, & name = "Impact_height_FSI", & long_name = "Impact height from FSI", & units = "m", & range = (/ -1.0E7_wp, 1.0E7_wp /), & DATA = impact_fsi(idx) ) CALL ropp_io_addvar_rodataD1d( ro_data, & name = "Bangle_FSI", & long_name = "Bending angle from FSI", & units = "rad", & range = (/ -1.0E7_wp, 1.0E7_wp /), & DATA = alpha_fsi(idx) ) CALL ropp_io_addvar_rodataD1d( ro_data, & name = "Amp_FSI", & long_name = "Amplitude spectrum from FSI", & units = "", & range = (/ -1.0E7_wp, 1.0E7_wp /), & DATA = amp_fsi(idx) ) ELSE CALL message( msg_error, 'FSI failed to generate any valid bending angles' ) END IF DEALLOCATE (amp_fsi, alpha_fsi, impact_fsi) END IF !------------------------------------------------------------------------------- ! 14. Populate ROprof structure with results and write out !------------------------------------------------------------------------------- IF ( msg_MODE == VerboseMode ) THEN ! 14.1 Add U_final to ro_data ! --------------------------- CALL ropp_io_addvar_rodataD1d( ro_data, & name = "Amp_U_final", & long_name = "Amplitude of complex signal on final screen", & units = "", & range = (/ -1.0E7_wp, 1.0E7_wp /), & DATA = ABS(U_final) ) CALL ropp_io_addvar_rodataD1d( ro_data, & name = "Phase_U_final", & long_name = "Phase of complex signal on final screen", & units = "m", & range = (/ -1.0E7_wp, 1.0E7_wp /), & DATA = cumulative_phase(U_final) / kval ) END IF ! 14.2 Initialise output structure so that it can be read by ropp_pp_occ_tool ! --------------------------------------------------------------------------- ro_data%Lev1a%dtime = time ro_data%Lev1a%phase_L1 = phase_leo ro_data%Lev1a%phase_L2 = phase_leo ! NB: At present, set L2 = L1 ro_data%Lev1a%snr_L1ca = amp_leo ro_data%Lev1a%snr_L1p = amp_leo ! NB: At present, set L1p = L1ca ro_data%Lev1a%snr_L2p = amp_leo ! NB: At present, set L2p = L1ca ro_data%Lev1a%r_leo(:, 1) = r_leo(:, 1) ro_data%Lev1a%r_leo(:, 2) = r_leo(:, 2) ro_data%Lev1a%r_leo(:, 3) = r_leo(:, 3) ro_data%Lev1a%r_gns(:, 1) = r_gps(:, 1) ro_data%Lev1a%r_gns(:, 2) = r_gps(:, 2) ro_data%Lev1a%r_gns(:, 3) = r_gps(:, 3) ro_data%DTocc%year = 2000 ; ro_data%DTpro%year = 2000 ro_data%DTocc%month = 01 ; ro_data%DTpro%month = 01 ro_data%DTocc%day = 01 ; ro_data%DTpro%day = 01 ro_data%DTocc%hour = 00 ; ro_data%DTpro%hour = 00 ro_data%DTocc%minute = 00 ; ro_data%DTpro%minute = 00 ro_data%DTocc%second = 00 ; ro_data%DTpro%second = 00 ro_data%DTocc%msec = 0000 ; ro_data%DTpro%msec = 0000 ro_data%georef%lat = 0 ! Equatorial occultation ro_data%PCD = 0 ! Unset the whole lot ro_data%PCD = IBSET(ro_data%PCD, PCD_offline) ! Not NRT ro_data%PCD = IBCLR(ro_data%PCD, PCD_rising) ! These are always setting occultations ! 14.3 Write out ! -------------- CALL message( msg_info, 'Writing computed excess phase to ' // TRIM(ADJUSTL(ofile)) ) CALL ropp_io_write(ro_data, ofile, append=.FALSE., ranchk=.FALSE.) !------------------------------------------------------------------------------- ! 15. Clean up !------------------------------------------------------------------------------- CALL ropp_io_free(ro_data) DEALLOCATE (U_final, U_init) ! (ny) DEALLOCATE (amp_leo, phase_leo, U_leo) ! (n_leo) DEALLOCATE (d2logndy2, log_refrac, zed) ! (nlev, nhoriz) DEALLOCATE (r_leo, r_gps) ! (n_leo, 3) DEALLOCATE (s_geom, theta, time) ! (n_leo) DEALLOCATE (xpos, ypos) ! (ny) CALL EXIT(msg_exit_status) CONTAINS !------------------------------------------------------------------------------- ! 16. Cumulative phase calculation !------------------------------------------------------------------------------- FUNCTION cumulative_phase(U) RESULT(phase) COMPLEX(wp), INTENT(in) :: U(:) REAL(wp) :: phase(SIZE(U)) REAL(wp), ALLOCATABLE :: phase_temp(:) INTEGER :: k ALLOCATE(phase_temp(SIZE(U))) phase_temp = 0.0_wp WHERE ( ABS(U) > TINY(0.0_wp) ) phase_temp = ATAN2( AIMAG(U), REAL(U) ) phase(1) = 0.0_wp DO k=2,SIZE(U) phase(k) = phase(k-1) + MODULO(phase_temp(k)-phase(k-1)+pi1, twopi) - pi1 END DO DEALLOCATE(phase_temp) END FUNCTION cumulative_phase !------------------------------------------------------------------------------- ! 17. Usage information !------------------------------------------------------------------------------- SUBROUTINE usage() PRINT *, 'Purpose:' PRINT *, ' 2D wave optics propagation tool' PRINT *, 'Usage:' PRINT *, ' > ropp_pp_wopt_2d_tool [] ' PRINT *, 'Options:' PRINT *, ' -o name of ROPP netCDF output file' PRINT *, ' (default: ropp_pp_wopt_2d_tool_out.nc)' PRINT *, ' -c name of file containing the configuration namelist' PRINT *, ' -f run FSI on excess phases to generate bangles' PRINT *, ' -d output additional diagnostics' PRINT *, ' -h this help' PRINT *, ' -v version information' PRINT *, '' END SUBROUTINE usage !------------------------------------------------------------------------------- ! 18. Version information !------------------------------------------------------------------------------- SUBROUTINE version_info() CHARACTER (LEN=40) :: version version = ropp_pp_version() PRINT *, 'ropp_pp_wopt_2d_tool - 2D wave optics propagation tool' PRINT *, '' PRINT *, 'This program is part of ROPP (PP) Release ' // TRIM(version) PRINT *, '' END SUBROUTINE version_info END PROGRAM ropp_pp_wopt_2D_tool