! $Id: t_iono_tl.f90 4010 2014-01-10 11:07:40Z idculv $

PROGRAM t_iono_tl

!****p* Programs/t_iono_tl *
!
! NAME
!    t_iono_tl - tests the compiled fm model and tangent linear functions
!
! SYNOPSIS
!    t_iono_tl
!
! DESCRIPTION
!    This program tests the bending angle operator when L1 and L2 are being
!    modelled directly.
!
! SEE ALSO
!   t_iono.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 typesizes, ONLY: wp => EightByteReal
  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(State1dFM)                                  :: state,state_new,state_tl
  TYPE(Obs1dBangle)                                :: obs_bangle,obs_bangle_new

  INTEGER                                          :: i, iargc, argc, k, ii, imax, imin
  INTEGER                                          :: n_profiles
  INTEGER                                          :: error

  LOGICAL                                          :: give_version
  LOGICAL                                          :: give_help

  REAL(wp)                                         :: cos_alpha, rel_err
  REAL(wp)                                         :: max_alpha, min_err
  REAL(wp), DIMENSION(:), ALLOCATABLE              :: y_tl, delta_y

  INTEGER, PARAMETER                               :: n_files=1
  CHARACTER(len = 4096), DIMENSION(n_files)        :: ifiles
  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' /)

!-------------------------------------------------------------------------------
! 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)

!-------------------------------------------------------------------------------
! 8. Copy data in RO structure to state
!-------------------------------------------------------------------------------

      ro_data%lev2c%direct_ion = .TRUE.

      CALL ropp_fm_roprof2state(ro_data, state)
      CALL ropp_fm_roprof2state(ro_data, state_new)
      CALL ropp_fm_roprof2state(ro_data, state_tl)

      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

!-------------------------------------------------------------------------------
! 9. Copy data in RO and refrac structure to bending angle obs vector
!-------------------------------------------------------------------------------

      CALL set_obs_levels_bangle(ro_data, obs_bangle, obs_bangle_new)

!-------------------------------------------------------------------------------
! 10. Calculate bending angle for original state
!-------------------------------------------------------------------------------

      CALL ropp_fm_bangle_1d(state, obs_bangle)

!-------------------------------------------------------------------------------
! 11. Set the perturbations
!-------------------------------------------------------------------------------

      CALL random_number(state_tl%temp(:))
      CALL random_number(state_tl%pres(:))
      CALL random_number(state_tl%geop(:))
      CALL random_number(state_tl%shum(:))

      state_tl%temp(:) = 0.0_wp
      state_tl%pres(:) = 0.0_wp
      state_tl%geop(:) = 0.0_wp
      state_tl%shum(:) = 0.0_wp

! only perturb ionospheric parameters

      CALL random_number(state_tl%ne_max)
      CALL random_number(state_tl%h_peak)
      CALL random_number(state_tl%h_width)

! scale the parameters

      state_tl%ne_max  = 100.0E11_wp*state_tl%ne_max
      state_tl%h_peak  = 100.0E05_wp*state_tl%h_peak
      state_tl%h_width = 100.0E05_wp*state_tl%h_width

      ALLOCATE(y_tl(SIZE(obs_bangle%bangle)), delta_y(SIZE(obs_bangle%bangle)))

      imax = -10000
      imin = 10000
      max_alpha = -10000.0_wp
      min_err = 10000.0_wp

      WRITE(*, '(A)') '    |dx|/|dx|_init    |dy=Kdx|    (Kdx|H(x+dx)-H(x))/|.||.|' // &
                      '   |H(x+dx)-H(x)-Kdx|/|H(x+dx)-H(x)|'

      DO ii = 1,15

! update perturbed state ...

        state_new%temp = state%temp + state_tl%temp
        state_new%pres = state%pres + state_tl%pres
        state_new%geop = state%geop + state_tl%geop
        state_new%shum = state%shum + state_tl%shum

! ... although only the ionospheric part has changed

        state_new%ne_max  = state%ne_max  + state_tl%ne_max
        state_new%h_peak  = state%h_peak  + state_tl%h_peak
        state_new%h_width = state%h_width + state_tl%h_width

! call operator with perturbed state

        CALL ropp_fm_bangle_1d(state_new, obs_bangle_new)

! call tangent linear with perturbation

        CALL ropp_fm_bangle_1d_tl(state, state_tl, obs_bangle, y_tl)

! save change in bangle

        delta_y = obs_bangle_new%bangle - obs_bangle%bangle

! angle between the vectors

        cos_alpha = DOT_PRODUCT(y_tl, delta_y) / &
                    SQRT(DOT_PRODUCT(y_tl, y_tl)*DOT_PRODUCT(delta_y, delta_y))

! relative error in the tangent linear

        rel_err = SQRT(DOT_PRODUCT(delta_y - y_tl, delta_y - y_tl)) / &
                  SQRT(DOT_PRODUCT(delta_y       , delta_y       ))

        WRITE (*, '(4(1PE20.10))') 10.0_wp**(1-ii), &
                                   SQRT(DOT_PRODUCT(delta_y, delta_y)), &
                                   cos_alpha, &
                                   rel_err

! look for max/minimum values

        IF ( cos_alpha > max_alpha ) THEN
          max_alpha = cos_alpha
          imax = ii
        ENDIF

        IF ( rel_err < min_err ) THEN
          min_err = rel_err
          imin = ii
        END IF

! reduce the size of the perturbation by a factor of 10

        state_tl%temp     = 0.1_wp*state_tl%temp
        state_tl%pres     = 0.1_wp*state_tl%pres
        state_tl%geop     = 0.1_wp*state_tl%geop
        state_tl%shum     = 0.1_wp*state_tl%shum

        state_tl%ne_max   = 0.1_wp*state_tl%ne_max
        state_tl%h_peak   = 0.1_wp*state_tl%h_peak
        state_tl%h_width  = 0.1_wp*state_tl%h_width

      END DO

      DEALLOCATE(y_tl, delta_y)

! check tangent linear

      IF (ABS(imin-imax) > 1) error = 1

!-------------------------------------------------------------------------------
! 12. Clean up
!-------------------------------------------------------------------------------

      CALL ropp_io_free(ro_data)

    END DO

  END DO

!  IF (error == 1) THEN
!     PRINT *,''
!     PRINT *,''
!     PRINT *,'*********************************'
!     PRINT *,'*** ropp_fm (t_iono_tl): FAIL ***'
!     PRINT *,'*********************************'
!     PRINT *,''
!  ELSE
!     PRINT *,''
!     PRINT *,''
!     PRINT *,'*********************************'
!     PRINT *,'*** ropp_fm (t_iono_tl): PASS ***'
!     PRINT *,'*********************************'
!     PRINT *,''
!  END IF

  CALL ropp_io_success(error/=1, 't_iono_tl_1', 'FM_TL L1 and L2')

CONTAINS

!-------------------------------------------------------------------------------
! 15. Calculate observation levels for bending angle
!-------------------------------------------------------------------------------

  SUBROUTINE set_obs_levels_bangle(ro_data, obs_bangle, obs_bangle_new)

!   15.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, obs_bangle_new

    REAL(wp) :: tmp

    INTEGER           :: i, n=400

!   15.2 Allocate arrays
!   --------------------

    obs_bangle%nobs = n
    obs_bangle%n_L1 = n / 2

! original

    ALLOCATE(obs_bangle%bangle(n))
    ALLOCATE(obs_bangle%impact(n))
    ALLOCATE(obs_bangle%weights(n))

! perturbed

    obs_bangle_new%nobs = n
    obs_bangle_new%n_L1 = n / 2

    ALLOCATE(obs_bangle_new%bangle(n))
    ALLOCATE(obs_bangle_new%impact(n))
    ALLOCATE(obs_bangle_new%weights(n))

!   15.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          = obs_bangle%r_curve + 800.0E3_wp

    obs_bangle_new%g_sfc      = obs_bangle%g_sfc
    obs_bangle_new%r_earth    = obs_bangle%r_earth
    obs_bangle_new%r_curve    = obs_bangle%r_curve
    obs_bangle_new%undulation = obs_bangle%undulation
    obs_bangle_new%r_leo      = obs_bangle%r_leo

    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

!   15.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

      obs_bangle_new%impact(i)                  =  tmp
      obs_bangle_new%impact(i+obs_bangle%n_L1)  =  tmp

    END DO

!   15.5 Fill other arrays
!   ----------------------

    obs_bangle%bangle(:)      = 0.0_wp
    obs_bangle%weights(:)     = 1.0_wp

    obs_bangle_new%bangle(:)  = 0.0_wp
    obs_bangle_new%weights(:) = 1.0_wp

  END SUBROUTINE set_obs_levels_bangle

!-------------------------------------------------------------------------------
! 16. Usage information
!-------------------------------------------------------------------------------

  SUBROUTINE usage()

    PRINT *, 't_iono_tl - Test local compilation of ropp_fm_iono_tl.'
    PRINT *, 'Valid options are:'
    PRINT *, '  -h                give (this) help.'
    PRINT *, '  -V                give some version information.'
    PRINT *, ''

  END SUBROUTINE usage

!-------------------------------------------------------------------------------
! 17. Version information
!-------------------------------------------------------------------------------

  SUBROUTINE version_info(version)
    CHARACTER(len = *) :: version
    PRINT *, 't_iono_tl - Test local compilation of ropp_fm_iono_tl.'
    PRINT *, ''
    PRINT *, 'This program is part of ROPP version ' // TRIM(version)
    PRINT *, ''
  END SUBROUTINE version_info

END PROGRAM t_iono_tl