! $Id: t_iono.f90 4010 2014-01-10 11:07:40Z idculv $ PROGRAM t_iono !****p* Programs/t_iono * ! ! NAME ! t_iono - tests the compiled fm model ! ! SYNOPSIS ! t_iono ! ! DESCRIPTION ! This program tests the bending angle operator when L1 and L2 are being ! modelled directly. ! ! SEE ALSO ! t_iono_tl.f90 ! t_iono_ad.f90 ! ! AUTHOR ! Met Office, Exeter, UK and ECMWF, Reading 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 ropp_utils USE ropp_io USE ropp_io_types, ONLY: ROprof USE ropp_fm USE ropp_fm_types USE ropp_fm_copy USE ropp_io_test IMPLICIT NONE TYPE(ROprof) :: ro_data TYPE(ROprof) :: ro_data_ref TYPE(State1dFM) :: state TYPE(Obs1dBangle) :: obs_bangle INTEGER :: i, iargc, argc, k, ii, j1, j2 INTEGER :: error INTEGER :: n_profiles REAL(wp) :: diff REAL(wp), PARAMETER :: tol=1.E-4_wp LOGICAL :: give_version LOGICAL :: give_help INTEGER, PARAMETER :: n_files=1 CHARACTER(len = 4096), DIMENSION(n_files) :: ifiles, rfiles CHARACTER(len = 256) :: buffer CHARACTER(len = 64) :: version !------------------------------------------------------------------------------- ! 2. Default settings !------------------------------------------------------------------------------- version = ropp_io_version() give_version = .FALSE. give_help = .FALSE. error = 0 ifiles = (/ '../data/bgr20090401_000329_M02_2030337800_N0007_YYYY.nc' /) rfiles = (/ '../data/bgr20090401_000329_M02_2030337800_N0007_YYYY_iono_ref.nc' /) !------------------------------------------------------------------------------- ! 3. Command line arguments !------------------------------------------------------------------------------- argc = iargc() i = 1 DO WHILE(i <= argc) CALL getarg(i, buffer) SELECT CASE (buffer) CASE('-h', '-help', '--help') ! Give some help give_help = .TRUE. CASE('-V', '-version', '--version') ! Give some version information give_version = .TRUE. CASE DEFAULT END SELECT i = i + 1 END DO IF (give_help) THEN CALL usage() END IF IF (give_version) THEN CALL version_info(version) END IF !------------------------------------------------------------------------------- ! 4. Set default units !------------------------------------------------------------------------------- CALL ropp_fm_set_units(ro_data) !------------------------------------------------------------------------------- ! 5. Loop over all input files !------------------------------------------------------------------------------- DO k = 1, n_files !------------------------------------------------------------------------------- ! 6. Loop over all profiles !------------------------------------------------------------------------------- n_profiles = ropp_io_nrec(ifiles(k)) DO i = 1, n_profiles !------------------------------------------------------------------------------- ! 7. Read data !------------------------------------------------------------------------------- CALL ropp_io_read(ro_data, ifiles(k), rec=i) CALL ropp_io_read(ro_data_ref, rfiles(k), rec=i) !------------------------------------------------------------------------------- ! 8. Copy data in RO structure to state and refrac obs vectors !------------------------------------------------------------------------------- ro_data%lev2c%direct_ion = .TRUE. CALL ropp_fm_roprof2state(ro_data, state) !------------------------------------------------------------------------------- ! 9. Copy data in RO and refrac structure to bending angle obs vector !------------------------------------------------------------------------------- CALL set_obs_levels_bangle(ro_data, obs_bangle) ! Reset to match the values in the reference file j1 = 1 ; j2 = obs_bangle%n_L1 obs_bangle%impact(j1:j2) = ro_data_ref%lev1b%impact j1 = obs_bangle%n_L1 + 1 ; j2 = 2*obs_bangle%n_L1 obs_bangle%impact(j1:j2) = ro_data_ref%lev1b%impact obs_bangle%r_curve = ro_data_ref%GEOref%roc obs_bangle%r_leo = ro_data_ref%GEOref%roc + 8.0E5_wp ! Matches value in other t_ion test files. obs_bangle%undulation = ro_data_ref%GEOref%undulation !------------------------------------------------------------------------------- ! 10. Calculate bending angle !------------------------------------------------------------------------------- CALL ropp_fm_bangle_1d(state, obs_bangle) WRITE(*, '(A)')' ii impact bangle_L1(test) bangle_L1(ref) bangle_L2(test) bangle_L2(ref)' DO ii = 1,obs_bangle%n_L1,10 WRITE (6,'(i6, f12.2, 5e18.8)') & ii,obs_bangle%impact(ii)-obs_bangle%r_curve, & obs_bangle%bangle(ii), & ro_data_ref%Lev1b%bangle_L1(ii), & obs_bangle%bangle(ii+obs_bangle%n_L1), & ro_data_ref%Lev1b%bangle_L2(ii) END DO ! If the new bangles are needed in a netCDF file. ! CALL ropp_fm_obs2roprof(obs_bangle, ro_data) ! CALL ropp_io_write(ro_data, file='out.nc') !------------------------------------------------------------------------------- ! 11. Compare bending angle data with reference !------------------------------------------------------------------------------- PRINT*, 'Comparing against bending angles in ' // TRIM(ADJUSTL(rfiles(k))) j1 = 1 ; j2 = obs_bangle%n_L1 diff = MAXVAL(ABS((ro_data_ref%Lev1b%bangle_L1 - & obs_bangle%bangle(j1:j2)) / & ro_data_ref%Lev1b%bangle_L1)) IF (diff > tol) THEN PRINT*, 'L1 bending angle comparison failed: max difference = ', & REAL(diff*100.0_wp, KIND(1.0)), '%' error = 1 END IF j1 = obs_bangle%n_L1 + 1 ; j2 = 2*obs_bangle%n_L1 diff = MAXVAL(ABS((ro_data_ref%Lev1b%bangle_L2 - & obs_bangle%bangle(j1:j2)) / & ro_data_ref%Lev1b%bangle_L2)) IF (diff > tol) THEN PRINT*, 'L2 bending angle comparison failed: max difference = ', & REAL(diff*100.0_wp, KIND(1.0)), '%' error = 1 END IF ! IF (error == 1) THEN ! PRINT *,'' ! PRINT *,'' ! PRINT *,'******************************' ! PRINT *,'*** ropp_fm (t_iono): FAIL ***' ! PRINT *,'******************************' ! PRINT *,'' ! ELSE ! PRINT *,'' ! PRINT *,'' ! PRINT *,'******************************' ! PRINT *,'*** ropp_fm (t_iono): PASS ***' ! PRINT *,'******************************' ! PRINT *,'' ! END IF CALL ropp_io_success(error/=1, 't_iono_1', 'FM L1 and L2') !------------------------------------------------------------------------------- ! 12. Clean up !------------------------------------------------------------------------------- CALL ropp_io_free(ro_data) CALL ropp_io_free(ro_data_ref) END DO END DO CONTAINS !------------------------------------------------------------------------------- ! 13. Calculate observation levels for bending angle !------------------------------------------------------------------------------- SUBROUTINE set_obs_levels_bangle(ro_data, obs_bangle) ! 13.1 Declarations ! ----------------- USE typesizes, ONLY: wp => EightByteReal USE ropp_io USE ropp_fm USE geodesy IMPLICIT NONE TYPE(ROprof) :: ro_data TYPE(Obs1dbangle) :: obs_bangle REAL(wp) :: tmp INTEGER :: i, n=600 ! 13.2 Allocate arrays ! -------------------- obs_bangle%nobs = n obs_bangle%n_L1 = n / 2 ALLOCATE(obs_bangle%bangle(n)) ALLOCATE(obs_bangle%impact(n)) ALLOCATE(obs_bangle%weights(n)) ! 13.3 Set scalar arguments of the observation vector ! --------------------------------------------------- obs_bangle%g_sfc = gravity(ro_data%GEOref%lat) obs_bangle%r_earth = R_eff(ro_data%GEOref%lat) IF ( ro_data%GEOref%roc > ropp_MDTV ) THEN obs_bangle%r_curve = ro_data%GEOref%roc ELSE obs_bangle%r_curve = 6.371E6_wp ENDIF IF ( ro_data%GEOref%undulation > ropp_MDTV ) THEN obs_bangle%undulation = ro_data%GEOref%undulation ELSE obs_bangle%undulation = 0.0_wp ENDIF ! set height of LEO to be 800 km for testing obs_bangle%r_leo = ro_data%GEOref%roc + 800.0E3_wp WRITE(*, '(A,3(1PE15.5))') 'Ne_max, H_peak, H_width = ', & ro_data%lev2c%ne_max,ro_data%lev2c%h_peak,ro_data%lev2c%h_width WRITE(*, '(A,4(1PE15.5),/)') 'RoC, r_LEO, und, lat = ', & obs_bangle%r_curve, obs_bangle%r_leo, & obs_bangle%undulation, ro_data%GEOref%lat ! 13.4 Calculate levels to coincide with the geopotential levels ! -------------------------------------------------------------- DO i = 1,obs_bangle%n_L1 tmp = 2500.0_wp + REAL(i-1, wp)*250.0_wp + & obs_bangle%r_curve + obs_bangle%undulation obs_bangle%impact(i) = tmp obs_bangle%impact(i+obs_bangle%n_L1) = tmp END DO IF (obs_bangle%undulation > ropp_MDTV) THEN obs_bangle%impact = obs_bangle%impact + obs_bangle%undulation ELSE CALL message(msg_warn, "Undulation missing. Will assume to " // & "be zero when calculating impact parameters.") END IF ! 13.5 Fill other arrays ! ---------------------- obs_bangle%bangle(:) = 0.0_wp obs_bangle%weights(:) = 1.0_wp END SUBROUTINE set_obs_levels_bangle !------------------------------------------------------------------------------- ! 14. Usage information !------------------------------------------------------------------------------- SUBROUTINE usage() PRINT *, 't_iono - Test local compilation of ropp_fm_iono.' PRINT *, 'Valid options are:' PRINT *, ' -h give (this) help.' PRINT *, ' -V give some version information.' PRINT *, '' END SUBROUTINE usage !------------------------------------------------------------------------------- ! 15. Version information !------------------------------------------------------------------------------- SUBROUTINE version_info(version) CHARACTER(len = *) :: version PRINT *, 't_iono - Test local compilation of ropp_fm_iono.' PRINT *, '' PRINT *, 'This program is part of ROPP version ' // TRIM(version) PRINT *, '' END SUBROUTINE version_info END PROGRAM t_iono