! $Id: bgrasc2ropp.f90 5714 2019-02-13 15:51:08Z idculv $

PROGRAM bgrasc2ropp

!****x* Programs/bgrasc2ropp *
!
! NAME
!   bgrasc2ropp
!
! SYNOPSIS
!   Convert ASCII file containing background profile data to ROPP netCDF
!
!   > bgrasc2ropp file_in [-o <file_out>] [-n] [-d] [-h] [-v]
!
! ARGUMENTS
!   file_in - ascii file containing Level 2b background data
!             in a Fortran namelist.
!             There is no default for this argument.
!
! OPTIONS
!  -o <file_out> - is the output netCDF file name.
!  -n - switches off range checking before writing ROPP file
!  -d - writes additional diagnostic information to stdout
!  -h - help
!  -v - version information
!
! INPUTS
!   Ascii file holding a Fortran namelist 'bgr_profile' containing those
!   elements of the {bg, GEOref, Lev2b, Lev2c and Lev2d} substructures 
!   of the ROprof structure needed by the ropp forward model. 
!   This namelist could be generated from ECMWF grib data by the sister tool
!   ropp_io/tools/grib2bgrasc.
!
! OUTPUTS
!   ROPP netCDF file, suitable for use by ropp_fm/tools/ropp_fm_bg2ro_1d
!
! MODULES
!   typesizes
!   messages
!   ropp_utils
!   ropp_io_types
!   ropp_io
!   ropp_io_ecmwf
!
! ERRORS
!   Program (shell) return codes:
!     0 = OK
!     1 = At least one Warning occurred
!     2 = At least one Error occurred
!     3 = A Fatal error occurred
!
! CALLS
!   usage
!   ropp_io_init
!   ropp_io_write
!   ropp_io_free
!   ropp_io_version
!   message
!   message_set_routine
!
! DESCRIPTION
!   Conversion of level bg/GEOref/2b/2c/2d background profile data in Fortran
!   namelist format to ROPP netCDF format. The input file is in the same format
!   as that produced by the sister program grib2bgrasc, which extracts a
!   profile from a GRIB2 format model dump.  The data is then written out to a
!   ROPP-standard netCDF file in PWD.
!
! RESTRICTIONS
!   1) Assumes ECMWF-like model level structures.
!   2) Users needing to read in more than 200 levels of data will have to
!      increase Nlevs_max and recompile.
!
! SEE ALSO
!   grib2bgrasc
!
! REFERENCES
!   1. ROPP User Guide - Part I: I/O module.
!      SAF/ROM/METO/UG/ROPP/002
!
! 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.
!
!****

  USE typesizes,     wp => EightByteReal
  USE messages
  USE DateTimeProgs, ONLY: CalToJul
  USE ropp_utils,    ONLY: WHERE, ropp_MDFV, ropp_MDTV, Get_IO_Unit
  USE ropp_io_types, ONLY: ROprof
  USE ropp_io,       ONLY: ropp_io_free,    &
                           ropp_io_init,    &
                           ropp_io_write,   &
                           ropp_io_version
  USE ropp_io_ecmwf

  IMPLICIT NONE

! Fixed parameters
  INTEGER                          :: namelist_lun      ! Namelist logical unit

! Local variables
  TYPE(ROprof)                     :: BGRprof           ! BG profile structure

  CHARACTER (LEN=256)              :: arg               ! Command line argument
  CHARACTER (LEN=256)              :: file_in           ! I/P file name
  CHARACTER (LEN=256)              :: file_out          ! O/P file name
  CHARACTER (LEN=80)               :: outstr            ! Formatted output string
  CHARACTER (LEN=10)               :: number            ! Number as string
  LOGICAL                          :: exists            ! File exists flag
  LOGICAL                          :: iranchk=.TRUE.    ! Range checking flag
  REAL(wp)                         :: swver             ! ROPP version

  INTEGER                          :: j, ih             ! Loop counters / indices
  INTEGER                          :: iarg              ! Command line argument index
  INTEGER                          :: narg              ! No. of command line arguments
  INTEGER                          :: iostatus          ! I/O Status
  REAL(KIND(1.D0))                 :: jul_day           ! Julian day of data and validity time
  INTEGER, DIMENSION(8)            :: CDT               ! Date/time combo for cal2jul

  INTEGER, DIMENSION(:), POINTER   :: idx => NULL()     ! Holds WHERE output (NB this is NOT the F90 intrinsic)
  INTEGER                          :: nidx=0            ! Number of elements if idx

  INTEGER, PARAMETER               :: Nlevs_max=200     ! Maximum number of model levels

  INTEGER                          :: NLevs             ! Number of full model levels
  INTEGER                          :: year, mon, day    ! Background validity date
  INTEGER                          :: hour, min, sec    ! Background validity time
  REAL(wp)                         :: tlead=ropp_MDFV   ! Background forecast range (hr)
  REAL(wp)                         :: lon=ropp_MDFV     ! Tangent point lon (deg)
  REAL(wp)                         :: lat=ropp_MDFV     ! Tangent point lat (deg)
  REAL(wp)                         :: und=ropp_MDFV     ! Tangent point undulation (m)
  REAL(wp)                         :: roc=ropp_MDFV     ! Tangent point radius of curvature (m)
  REAL(wp), DIMENSION(3)           :: coc=(/ropp_MDFV, ropp_MDFV, ropp_MDFV/) ! Tangent point centre of curvature (m)
  REAL(wp)                         :: azi=ropp_MDFV     ! GNSS-->LEO line of sight angle (degT)

  REAL(wp)                         :: Z0=ropp_MDFV      ! Surface geopotential height (m)
  REAL(wp)                         :: P0=ropp_MDFV      ! Surface pressure (Pa)
  REAL(wp), DIMENSION(Nlevs_max)   :: P=ropp_MDFV       ! Pressure (Pa)
  REAL(wp), DIMENSION(Nlevs_max)   :: T=ropp_MDFV       ! Temperature (K)
  REAL(wp), DIMENSION(Nlevs_max)   :: Q=ropp_MDFV       ! Specific humidity (g/kg)
  REAL(wp), DIMENSION(Nlevs_max)   :: Z=ropp_MDFV       ! Geopotential height (m)
  REAL(wp), DIMENSION(Nlevs_max+1) :: Ak=ropp_MDFV      ! Hybrid/Eta level A-coefficient (Pa)
  REAL(wp), DIMENSION(Nlevs_max+1) :: Bk=ropp_MDFV      ! Hybrid/Eta level B-coefficient

  REAL(wp), DIMENSION(Nlevs_max)   :: Tvirt=ropp_MDFV   ! Virtual temperature (K)
  REAL(wp), DIMENSION(Nlevs_max+1) :: Ph=ropp_MDFV      ! Pressure on half levs (hPa)
  REAL(wp), DIMENSION(Nlevs_max)   :: deltaP=ropp_MDFV  ! Change in pressure between half levs (hPa)
  REAL(wp), DIMENSION(Nlevs_max)   :: lnP=ropp_MDFV     ! Log of pressure ratio between half levs
  REAL(wp), DIMENSION(Nlevs_max)   :: alpha=ropp_MDFV   ! Verical interpolation coefficient
  REAL(wp), DIMENSION(Nlevs_max)   :: dZ=ropp_MDFV      ! Function to be interpolated
  REAL(wp), DIMENSION(Nlevs_max+1) :: Zh=ropp_MDFV      ! Geopotential height on half levs (m)

! Constants
  REAL(wp), PARAMETER              :: g_wmo=9.80665_wp  ! Standard gravity (m/s2)
  REAL(wp), PARAMETER              :: R_dry=287.0597_wp ! Dry gas constant (K/kg/K)

! Input to this program
  NAMELIST / bgr_profile / year, mon, day, hour, min, sec, tlead, &
                           lon, lat, und, roc, coc, azi, &
                           Z0, P0, &
                           Nlevs, &
                           P, T, Q, Z, Ak, Bk

!NB: If Z is not in namelist it can be calculated from {Z0, P0, T, Q, Ak, Bk}

! Some compilers may need the following declaration to be commented out
  INTEGER             :: IARGC

!-------------------------------------------------------------
! 1. Initialise
!-------------------------------------------------------------

  CALL message_set_routine ( 'bgrasc2ropp' )

  CALL message(msg_noin, '')
  CALL message(msg_noin, &
       '---------------------------------------------------------------------')
  CALL message(msg_noin, &
       '                     bgrasc to ROPP converter'                        )
  CALL message(msg_noin, &
       '---------------------------------------------------------------------')
  CALL message(msg_noin, '')

  namelist_lun = Get_IO_Unit()

!-------------------------------------------------------------
! 2. Parse command line options
!-------------------------------------------------------------

  narg = IARGC()

  file_in   = ' '       ! no default for i/p file name
  file_out  = ' '       ! assume a default generated from i/p file name

  iarg = 1
  DO WHILE ( iarg <= narg )

    CALL GETARG ( iarg, arg )

    SELECT CASE (arg)
      CASE ('-d', '-D', '--debug')
        msg_MODE = VerboseMode

      CASE ('-h', '-H', '--help', '?')
        narg = 0
        file_in = 'dummy'

      CASE ('-o', '-O', '--output')
        iarg = iarg + 1
        CALL GETARG ( iarg, arg )
        file_out = arg

      CASE ('-n', '-N', '--no-ranchk')
        iranchk = .FALSE.

      CASE ('-v', '-V', '--version')
        CALL version_info()
        CALL EXIT(msg_exit_ok)

      CASE DEFAULT
         IF ( arg(1:1) /= '-' ) THEN
           file_in = arg
         END IF
    END SELECT

    iarg = iarg + 1

  END DO

  IF ( file_in == ' ' ) THEN
    CALL message ( msg_error, 'No input file(s) specified' )
    narg = 0
  END IF

  IF ( narg == 0 ) THEN
    CALL Usage
    CALL EXIT(msg_exit_status)
  END IF

!-------------------------------------------------------------
! 3. Check text input file exists;
!    make output file name if not given on command line
!-------------------------------------------------------------

  INQUIRE ( FILE=file_in, EXIST=exists )
  IF ( .NOT. exists ) &
    CALL message ( msg_fatal, 'Input namelist file ' // TRIM(file_in) // &
                              ' not found' )

  IF ( file_out == ' ' ) THEN
    j = INDEX(file_in, '/', BACK=.TRUE.)
    file_out = TRIM(ADJUSTL(file_in(j+1:))) // '.nc'
  END IF

  IF (msg_MODE  == VerboseMode) THEN
    CALL message ( msg_diag, 'file_in = ' // file_in )
    CALL message ( msg_diag, 'file_out = ' // file_out )
  ENDIF

!-------------------------------------------------------------
! 4. Read namelist
!-------------------------------------------------------------

  CALL message ( msg_info, 'Reading file ' // TRIM(file_in) )

  OPEN ( UNIT=namelist_lun, FILE=file_in, STATUS='OLD', ACTION='READ' )

  READ ( UNIT=namelist_lun, NML=bgr_profile, iostat=iostatus)

  IF ( iostatus > 0 ) THEN
    CALL message ( msg_fatal, 'I/O error while reading text file' )
  ELSE
    IF (msg_MODE == VerboseMode) THEN
      CALL message( msg_diag, 'Contents of bgr_profile namelist:' )
      WRITE (*, NML=bgr_profile)  ! For diagnostic purposes
    ENDIF
  ENDIF

  CLOSE ( UNIT=namelist_lun )

!-------------------------------------------------------------
! 4.1 Check namelist
!-------------------------------------------------------------

  IF (Nlevs > Nlevs_max) &
    CALL message ( msg_fatal, 'Too many levels. Increase Nlevs_max and recompile' )

  IF (Nlevs <= 0) CALL message ( msg_fatal, 'Nlevs missing' )

  IF (P0 < ropp_MDTV) CALL message ( msg_fatal, 'Surface pressure missing' )

  IF (Z0 < ropp_MDTV) CALL message ( msg_fatal, 'Surface GPH missing' )

  IF (ALL(P < ropp_MDTV)) CALL message ( msg_fatal, 'Pressure profile missing' )
  idx => WHERE(P >= ropp_MDTV, nidx)
  IF (nidx /= Nlevs) CALL message ( msg_fatal, 'Wrong number of pressures' )

  IF (ALL(T < ropp_MDTV)) CALL message ( msg_fatal, 'Temperature profile missing' )
  idx => WHERE(T >= ropp_MDTV, nidx)
  IF (nidx /= Nlevs) CALL message ( msg_fatal, 'Wrong number of temperatures' )

  IF (ALL(Q < ropp_MDTV)) CALL message ( msg_fatal, 'Specific humidity profile missing' )
  idx => WHERE(Q >= ropp_MDTV, nidx)
  IF (nidx /= Nlevs) CALL message ( msg_fatal, 'Wrong number of specific humidities' )

! If possible, read Ak from namelist; otherwise use standard values, if possible.

  idx => WHERE(Ak >= ropp_MDTV, nidx)
  IF (nidx /= (Nlevs+1)) THEN

    CALL message ( msg_warn, 'Some Aks missing ... will use ' // &
                             'standard values if possible' )

    SELECT CASE (Nlevs)
      CASE (60)
        Ak(1:Nlevs+1) = Ak60
        Ak(Nlevs+1)   = 1.0e-32_wp
      CASE (91)
        Ak(1:Nlevs+1) = Ak91
        Ak(Nlevs+1)   = 1.0e-32_wp
      CASE (137)
        Ak(1:Nlevs+1) = Ak137
        Ak(Nlevs+1)   = 1.0e-32_wp
      CASE DEFAULT
        CALL message ( msg_warn, 'Non-standard number of levels ... ' // &
                                 'setting all Ak missing' )
        Ak(1:Nlevs+1) = ropp_MDFV
    END SELECT

  ENDIF

! If possible, read Bk from namelist; otherwise use standard values, if possible.

  idx => WHERE(Bk >= ropp_MDTV, nidx)
  IF (nidx /= (Nlevs+1)) THEN

    CALL message ( msg_warn, 'Some Bks missing ... will use ' // &
                             'standard values if possible' )

    SELECT CASE (Nlevs)
      CASE (60)
        Bk(1:Nlevs+1) = Bk60
      CASE (91)
        Bk(1:Nlevs+1) = Bk91
      CASE (137)
        Bk(1:Nlevs+1) = Bk137
      CASE DEFAULT
        CALL message ( msg_warn, 'Non-standard number of levels ... ' // &
                                 'setting all Bk missing' )
        Bk(1:Nlevs+1) = ropp_MDFV
    END SELECT

  ENDIF

! If necessary, recalculate geopotential here, or rangechecking will remove the whole lev1b structure.
! Calculate using Z(p)-Z(p*) = int_from(p*)_to_(p) (-RTv(p')/g0) dlogp'.

  IF (ANY(Z(1:Nlevs) < ropp_MDTV)) THEN 

    !--- Virtual temperatures
    Tvirt  = (1.0_wp + 0.61_wp * Q/1000.0_wp) * T

    !--- Pressure differences
    Ph = Ak + Bk*P0
    deltaP = Ph(1:Nlevs) - Ph(2:Nlevs+1)

    !--- Log of pressure ratio
    lnP = log(Ph(1:Nlevs) / Ph(2:Nlevs+1))

    !--- Interpolation coefficients
    alpha = 1.0_wp - Ph(2:Nlevs+1)/deltaP * lnP
    alpha(Nlevs) = log(2.0_wp)

    !--- Function to be integrated
    dZ = R_dry * Tvirt * lnP / g_wmo

    !--- Calculate geopotential height integral
    Zh(1) = 0.0_wp
    DO ih = 2,Nlevs+1
      Zh(ih) = SUM(dZ(1:ih-1))
    ENDDO
    IF (Z0 > ropp_MDTV) Zh(:) = Zh(:) + Z0

    !--- Interpolate onto full levels
    Z(1:Nlevs) = Zh(1:Nlevs) + alpha(1:Nlevs) * R_dry * Tvirt(1:Nlevs) / g_wmo

  ENDIF

  idx => WHERE(Z >= ropp_MDTV, nidx)
  IF (nidx /= Nlevs) CALL message ( msg_fatal, 'Wrong number of GPHs' )

! Checks on ranges could be put in here, although the whole BGprof structure
! will be range-checked in the call to ropp_io_write below. 

!-------------------------------------------------------------
! 5.0 Initialise ROPP structures & fill in static header data
!     if not in text files
!-------------------------------------------------------------

  CALL ropp_io_init ( BGRprof, 0, 0, 0, NLevs, 1, NLevs+1 )

!-------------------------------------------------------------
! 5.1 Populate ROPP structure with data from namelist
!-------------------------------------------------------------

  BGRprof%bg%source            = 'ECMWF'
  BGRprof%bg%Year              = year
  BGRprof%bg%Month             = mon
  BGRprof%bg%Day               = day
  BGRprof%bg%Hour              = hour
  BGRprof%bg%Minute            = min
!  BGRprof%bg%Second   = sec  ! Not currently part of the ROprof%bg substructure
  BGRprof%bg%Fcperiod          = tlead

  BGRprof%GEOref%lat           = lat
  BGRprof%GEOref%lon           = lon
  BGRprof%GEOref%r_coc         = coc ! depends on azimuth
  BGRprof%GEOref%roc           = roc ! depends on azimuth
  BGRprof%GEOref%azimuth       = azi
  BGRprof%GEOref%undulation    = und

  BGRprof%Lev2b%Npoints        = Nlevs
  BGRprof%Lev2b%press          = P(1:Nlevs)*1.0e-2_wp ! Pa to hPa for ROPP
  BGRprof%Lev2b%temp           = T(1:Nlevs)
  BGRprof%Lev2b%shum           = Q(1:Nlevs)
  BGRprof%Lev2b%geop           = Z(1:Nlevs)
  BGRprof%Lev2b%meteo_qual     = 100.0_wp

  BGRprof%Lev2c%Npoints        = 1
  BGRprof%Lev2c%geop_sfc       = Z0
  BGRprof%Lev2c%press_sfc      = P0*1.0e-2_wp ! Pa to hPa for ROPP
  BGRprof%Lev2c%press_sfc_qual = 100.0_wp

  BGRprof%Lev2d%Npoints        = Nlevs+1
  BGRprof%Lev2d%level_type     = 'ECMWF'
  BGRprof%Lev2d%level_coeff_a  = ak(1:Nlevs+1)*1.0e-2_wp ! Pa to hPa for ROPP
  BGRprof%Lev2d%level_coeff_b  = bk(1:Nlevs+1)

!-------------------------------------------------------------
! 5.2 Convert text parameters to ROPP standard
!-------------------------------------------------------------

  BGRprof%software_version = ropp_io_version()
    IF (INDEX(BGRprof%software_version, ' ') > 1) THEN
      READ (BGRprof%software_version(2:INDEX(BGRprof%software_version, ' ')-1), *) swver
      IF (INT(swver) < 10) THEN
        BGRprof%software_version = 'V0' // BGRprof%software_version(2:INDEX(BGRprof%software_version, ' ')-1)
      ELSE
        BGRprof%software_version = 'V'  // BGRprof%software_version(2:INDEX(BGRprof%software_version, ' ')-1)
      ENDIF
      IF (LEN_TRIM(BGRprof%software_version) < 8) THEN
        DO j=LEN_TRIM(BGRprof%software_version),8
          BGRprof%software_version(j:j) = '0'
        ENDDO
      ENDIF
    ENDIF

  WRITE ( outstr, FMT='(I4,"/",I2,"/",I2," ",I4,"Z")' ) &
    BGRprof%bg%Year, BGRprof%bg%Month, BGRprof%bg%Day, 100*BGRprof%bg%Hour+BGRprof%bg%Minute
  CALL message ( msg_diag, '   Verification time            : '//TRIM(outstr) )

  CDT = (/ BGRprof%bg%Year, BGRprof%bg%Month, BGRprof%bg%Day, 0, BGRprof%bg%Hour, BGRprof%bg%Minute, 0, 0 /)
  CALL CalToJul( CDT,  jul_day, 1 )
  jul_day = jul_day - (BGRprof%bg%Fcperiod/24.d0)
  CALL CalToJul( CDT,  jul_day, -1 )
  WRITE ( outstr, FMT='(I4,"/",I2,"/",I2," ",I4,"Z")' ) &
    CDT(1), CDT(2), CDT(3), 100*CDT(5)+CDT(6)
  CALL message ( msg_diag, '   Data time                    : '//TRIM(outstr) )

  WRITE ( outstr, FMT='(F6.2,",",F7.2)' ) BGRprof%GeoRef%Lat, BGRprof%GeoRef%Lon
  CALL message ( msg_diag, '   Latitude,Longitude           : '//TRIM(outstr) )

  WRITE ( number, FMT='(I6)' ) BGRprof%Lev2a%Npoints
  CALL message ( msg_diag, '   No. of refractivity samples  : '//TRIM(number) )

  WRITE ( number, FMT='(I6)' ) BGRprof%Lev2b%Npoints
  CALL message ( msg_diag, '   No. of geophysical samples   : '//TRIM(number) )

  WRITE ( number, FMT='(I6)' ) BGRprof%Lev2c%Npoints
  CALL message ( msg_diag, '   No. of surface geo. samples  : '//TRIM(number) )

  WRITE ( number, FMT='(I6)' ) BGRprof%Lev2d%Npoints
  CALL message ( msg_diag, '   No. of model coeff. levels   : '//TRIM(number) )

!-------------------------------------------------------------
! 6. Write ROPP netCDF file
!-------------------------------------------------------------

  CALL message ( msg_info, 'Writing ' // TRIM(file_out) )

  CALL ropp_io_write ( BGRprof, file=file_out, ierr=iostatus , ranchk=iranchk)

  IF ( iostatus > 0 ) THEN
    CALL message ( msg_fatal, 'Error while writing output file' )
  END IF

!-------------------------------------------------------------
! 7. Tidy up - deallocate structures & free memory
!-------------------------------------------------------------

  CALL ropp_io_free ( BGRprof )
  CALL message ( msg_noin, ' ' )
  CALL EXIT(msg_exit_status)

CONTAINS

!-------------------------------------------------------------------------------
! 8. Usage (help) information
!-------------------------------------------------------------------------------

  SUBROUTINE Usage()
    PRINT *, 'Purpose:'
    PRINT *, '  Convert an ascii file containing BGR model data, '
    PRINT *, '  in Fortran namelist format to a ROPP netCDF file'
    PRINT *, 'Usage:'
    PRINT *, '  > bgrasc2ropp file_in [-o <file_out>] [-n] [-d] [-h] [-v]'
    PRINT *, '  where:'
    PRINT *, '    file_in is an ascii file holding the Fortran namelist.'
    PRINT *, 'Options:'
    PRINT *, '  -o <file_out> is the output (netCDF) file name'
    PRINT *, '  -n switches off range checking before writing ROPP file'
    PRINT *, '  -d prints out some additional diagnostics to stdout'
    PRINT *, '  -h this help'
    PRINT *, '  -v version information'
    PRINT *, 'Defaults:'
    PRINT *, '  Input  file name : required'
    PRINT *, '  Output file name : $PWD/file_in.nc'
    PRINT *, 'See bgrasc2ropp(1) for details.'
    PRINT *, ''
  END SUBROUTINE Usage

!-------------------------------------------------------------------------------
! 9. Version information
!-------------------------------------------------------------------------------

  SUBROUTINE version_info()
    CHARACTER (LEN=40) :: version
    version = ropp_io_version()
    PRINT *, 'bgrasc2ropp - convert ascii files to ROPP netCDF'
    PRINT *, ''
    PRINT *, 'This program is part of ROPP (IO) Release ' // TRIM(version)
    PRINT *, ''
  END SUBROUTINE version_info

END Program bgrasc2ropp