! $Id: ropp_fm_bg2ro_2d.f90 4215 2014-06-03 08:40:22Z idculv $

PROGRAM ropp_fm_bg2ro_2d

!****p* Programs/ropp_fm_bg2ro_2d *
!
! NAME
!   ropp_fm_bg2ro_2d
!
! SYNOPSIS
!   Calculate radio occultation pseudo observation from background model
!   data using 2d forward models
!
!   > ropp_fm_bg2ro_2d <infile(s)> -o <outfile> [-comp] [-check_qsat] [-d] [-h] [-v]
!
! ARGUMENTS
!   <infile(s)>   One (or more) input file names.
!
! OPTIONS
!   -o <outfile>  name of output file (default: bg2ro_2d.nc)
!   -comp         include non-ideal gas compressibility
!   -check_qsat   check against saturation
!   -d            output additional diagnostics
!   -h            help
!   -v            version information
!
! DESCRIPTION
!   This program reads 2D model data on model levels from the input data files
!   and calculates vertical profiles of bending angle using the 2d forward 
!   operators. The result is written to an ROPP-formatted 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.
!
! EXAMPLES
!   1) To calculate 2D bending angle from one of the example
!      (single-) files in the data directory, with extra diagnostics:
!
!       > ropp_fm_bg2ro_2d -d ../data/ECMWF_2D_DATA.nc
!
!   2) To calculate 2D bending angle from one of the example
!      multifiles in the data directory, with the non-ideal-gas option:
!
!       > ropp_fm_bg2ro_2d -comp ../data/ECMWF_10_OCCS.nc
!
! SEE ALSO
!   ropp_fm_bangle_2d
!
! 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 ropp_utils
  USE ropp_io
  USE ropp_io_types, ONLY: ROprof2d
  USE ropp_fm
  USE ropp_fm_types
  USE ropp_fm_copy

  IMPLICIT NONE

  TYPE(ROprof2d)                                   :: ro_data
  TYPE(State2dFM)                                  :: state
  TYPE(Obs1dBangle)                                :: obs_bangle
  TYPE(State1dFM)                                  :: state_1d

  INTEGER                                          :: idummy
  INTEGER                                          :: i, iargc, argc, k
  INTEGER                                          :: n_files, n_profiles

  LOGICAL                                          :: give_help
  LOGICAL                                          :: ranchk = .TRUE.
  LOGICAL                                          :: compress   = .FALSE.
  LOGICAL                                          :: check_qsat = .FALSE.
  LOGICAL                                          :: check_qmin = .TRUE.

  CHARACTER(len = 4096), DIMENSION(:), ALLOCATABLE :: ifiles
  CHARACTER(len = 4096)                            :: ofile
  CHARACTER(len =  256)                            :: buffer
  CHARACTER(len =    4)                            :: istr
  CHARACTER(len =    6)                            :: nstr

!-------------------------------------------------------------------------------
! 2. Default settings
!-------------------------------------------------------------------------------

  give_help = .FALSE.
  ofile     = "bg2ro_2d.nc"

  CALL message(msg_noin, '')
  CALL message(msg_noin, &
      '-----------------------------------------------------------------------')
  CALL message(msg_noin, &
      '                  ROPP Forward Model (2D operator)'                     )
  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 ('--no-ranchk')                ! Use no rangecheck on output
           ranchk = .FALSE.
        CASE ('-comp')                      ! Non-ideal gas switch
           compress = .TRUE.
        CASE ('-check_qsat')                ! Check against saturation
           check_qsat = .TRUE.
        CASE ('-nocheck_qmin')              ! Do not check against dryness
           check_qmin = .FALSE.
        CASE('-d')                          ! Additional diagnostic mode
           msg_MODE = VerboseMode
        CASE('-h','--help','?')             ! 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
             n_files = n_files + 1
             ifiles(n_files) = buffer
           END IF
    END SELECT
    i = i + 1
  END DO

  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. Remove pre-existing output file
!-------------------------------------------------------------------------------

  CALL file_delete(ofile, idummy)

!-------------------------------------------------------------------------------
! 5. Set default units
!-------------------------------------------------------------------------------

  CALL ropp_fm_set_units(ro_data)

!-------------------------------------------------------------------------------
! 6. Loop over all input files
!-------------------------------------------------------------------------------

  DO k = 1, n_files

!-------------------------------------------------------------------------------
! 7. Loop over all profiles
!-------------------------------------------------------------------------------

    n_profiles = ropp_io_nrec(ifiles(k))

    DO i = 1, n_profiles

      WRITE(istr, '(i4)') i
      WRITE(nstr, '(i6)') n_profiles
      CALL message(msg_noin, '')
      CALL message(msg_info, "Processing profile " // istr // " of " // nstr )

!-------------------------------------------------------------------------------
! 8. Read data
!-------------------------------------------------------------------------------

      CALL ropp_io_read(ro_data, ifiles(k), rec = i, ranchk = ranchk)
      CALL message(msg_info, "(" // TRIM(ro_data%occ_id) // ") \n")

!-------------------------------------------------------------------------------
! 9. Copy data in RO structure to state and obs vectors
!-------------------------------------------------------------------------------

!  state%check_qsat and state%check_qmin are both used by ropp_fm_roprof2state,
!  and must therefore be set before calling it.

      IF (check_qsat) state%check_qsat = .TRUE. ! Switch on check against saturation

      IF (.NOT. check_qmin) state%check_qmin = .FALSE. ! Switch off check against dryness

      CALL ropp_fm_roprof2state(ro_data, state)

      IF (compress) state%non_ideal = .TRUE.  ! Switch on non-ideal gas

      IF (state%state_ok) THEN

!-------------------------------------------------------------------------------
! 10. Copy data from ROprof2D structure to bending angle obs vector
!-------------------------------------------------------------------------------

        IF (ro_data%lev1b%Npoints > 0) THEN

          CALL ropp_fm_roprof2obs(ro_data, obs_bangle)   !!set to obs levels

        ELSE

          CALL set_obs_levels_bangle(ro_data, obs_bangle) !default

        ENDIF

!-------------------------------------------------------------------------------
! 11. Calculate bending angle (2D operator)
!-------------------------------------------------------------------------------

        CALL ropp_fm_bangle_2d(state, obs_bangle)

!-------------------------------------------------------------------------------
! 12. Copy simulated observations to RO structure
!-------------------------------------------------------------------------------

        CALL ropp_fm_obs2roprof(obs_bangle, ro_data)

!-------------------------------------------------------------------------------
! 13. Perform 1d operation for comparison (saved in bangle_opt)
!-------------------------------------------------------------------------------

        CALL ropp_fm_roprof2state(ro_data, state_1d)

        IF (compress) state%non_ideal = .TRUE.

        CALL ropp_fm_bangle_1d(state_1d, obs_bangle)

        ro_data%lev1b%impact_opt = obs_bangle%impact
        ro_data%lev1b%bangle_opt = obs_bangle%bangle

!-------------------------------------------------------------------------------
! 14. Update RO structure with computed state vector variables
!-------------------------------------------------------------------------------

        CALL ropp_fm_state2roprof(state, ro_data)

      ENDIF

!-------------------------------------------------------------------------------
! 15. Write data
!-------------------------------------------------------------------------------

      CALL ropp_io_write(ro_data, ofile, append = .TRUE., ranchk = ranchk )

!-------------------------------------------------------------------------------
! 16. Clean up
!-------------------------------------------------------------------------------

      CALL ropp_fm_free(state)
      CALL ropp_fm_free(state_1d)
      CALL ropp_fm_free(obs_bangle)
      CALL ropp_io_free(ro_data)

    END DO ! Loop over profiles

  END DO ! Loop over files

  CALL EXIT(msg_exit_status)

CONTAINS

!-------------------------------------------------------------------------------
! 17. Calculate obs levels for bending angle
!-------------------------------------------------------------------------------

  SUBROUTINE set_obs_levels_bangle(ro_data, obs_bangle)

!   17.1 Declarations
!   -----------------

    USE typesizes, ONLY: wp => EightByteReal
    USE geodesy
    USE ropp_io_types, ONLY: ROprof2d
    USE ropp_fm

    IMPLICIT NONE

    TYPE(ROprof2d)    :: ro_data
    TYPE(Obs1dbangle) :: obs_bangle

    INTEGER           :: n,i

!   17.2 Allocate arrays
!   --------------------

    n = 250

    obs_bangle%nobs = n

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

!   17.3 Set scalar arguments of the observation vector
!   ---------------------------------------------------

    IF (ro_data%GEOref%roc > 0.0_wp) THEN
      obs_bangle%r_curve = ro_data%GEOref%roc
      obs_bangle%undulation = ro_data%GEOref%undulation
    ENDIF

!   17.4 Calculate levels to coincide with the geopotential levels
!   --------------------------------------------------------------

    DO i = 1, n

       obs_bangle%impact(i)  =  &
       obs_bangle%r_curve+ obs_bangle%undulation + REAL(i)*200.0_wp + 2500.0_wp

    ENDDO

!   17.5 Fill other arrays
!   ----------------------

    obs_bangle%bangle(:)  = 0.0_wp

  END SUBROUTINE set_obs_levels_bangle

!-------------------------------------------------------------------------------
! 18. Usage information
!-------------------------------------------------------------------------------

  SUBROUTINE usage()
    PRINT *, 'Purpose:'
    PRINT *, '  2D Bending angle forward model.'
    PRINT *, 'Usage:'
    PRINT *, '  > ropp_fm_bg2ro_2d [<options>] <input_file(s)>'
    PRINT *, 'Options:'
    PRINT *, '  -o <output_file>  name of ROPP netCDF output file'
    PRINT *, '  -comp             include non-ideal gas compressibility'
    PRINT *, '  -check_qsat       check against supersaturation'
    PRINT *, '  -nocheck_qmin     do not check for superdryness'
    PRINT *, '  -d                output additional diagnostics'
    PRINT *, '  -h                this help'
    PRINT *, '  -v                version information'
    PRINT *, ''
  END SUBROUTINE usage

!-------------------------------------------------------------------------------
! 19. Version information
!-------------------------------------------------------------------------------

  SUBROUTINE version_info()
    CHARACTER (LEN=40) :: version
    version = ropp_fm_version()
    PRINT *, 'ropp_fm_bg2ro_2d - 2D bending angle forward model.'
    PRINT *, ''
    PRINT *, 'This program is part of ROPP (FM) Release ' // TRIM(version)
    PRINT *, ''
  END SUBROUTINE version_info

END PROGRAM ropp_fm_bg2ro_2d