ufo_crtm_utils_mod.F90

Go to the documentation of this file.

! (C) Copyright 2018 UCAR
!
! This software is licensed under the terms of the Apache Licence Version 2.0
! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
 
!> Fortran module to provide code shared between nonlinear and tlm/adm radiance calculations
 
MODULE ufo_crtm_utils_mod
 
use fckit_configuration_module, only: fckit_configuration
use iso_c_binding
use kinds
 
use crtm_module
 
use ufo_geovals_mod, only: ufo_geovals, ufo_geoval, ufo_geovals_get_var
use ufo_vars_mod
use obsspace_mod
use ufo_utils_mod, only: cmp_strings
 
implicit none
private
 
public crtm_comm_stat_check
public crtm_conf
public crtm_conf_setup
public crtm_conf_delete
public get_var_name
public load_atm_data
public load_sfc_data
public load_geom_data
public ufo_crtm_skip_profiles
 
PUBLIC load_aerosol_data
public assign_aerosol_names
public calculate_aero_layer_factor
PUBLIC :: qsmith
PUBLIC :: upper2lower
 
PUBLIC :: grav,rv_rd,&
     &aerosol_concentration_minvalue,aerosol_concentration_minvalue_layer
 
REAL(kind_real), PARAMETER :: &
     &rdgas = 2.8704e+2_kind_real,&
     &rvgas = 4.6150e+2_kind_real,&
     &rv_rd = rvgas/rdgas,&
     &esl = 0.621971831,&
     &zvir =  rv_rd - 1_kind_real,&
     &tice = 273.16_kind_real,&
     &grav = 9.81_kind_real,&
     &aerosol_concentration_minvalue=1.e-16_kind_real,&
     &aerosol_concentration_minvalue_layer=tiny(rdgas),& 
     &ozone_default_value=1.e-3_kind_real ! in ppmv in crtm
 
integer, parameter, public :: max_string=800
 
!Type for general config
type crtm_conf
 integer :: n_sensors
 integer :: n_absorbers
 integer :: n_clouds
 integer :: n_aerosols
 integer :: n_surfaces
 character(len=MAXVARLEN), allocatable :: absorbers(:)
 integer, allocatable :: absorber_id(:)
 integer, allocatable :: absorber_units(:)
 character(len=MAXVARLEN), allocatable :: clouds(:,:)
 integer, allocatable :: cloud_id(:)
 character(len=MAXVARLEN), allocatable :: surfaces(:)
 
 character(len=255), allocatable :: sensor_id(:)
 character(len=255) :: endian_type
 character(len=255) :: coefficient_path
 character(len=255) :: &
    irwatercoeff_file, irlandcoeff_file, irsnowcoeff_file, iricecoeff_file, &
    viswatercoeff_file, vislandcoeff_file, vissnowcoeff_file, visicecoeff_file, &
    mwwatercoeff_file
 integer, allocatable :: land_wsi(:)
 real(kind_real) :: cloud_fraction = -1.0_kind_real
 integer :: inspect
 character(len=MAXVARLEN) :: aerosol_option
 character(len=255) :: salinity_option
  character(len=MAXVARLEN) :: sfc_wind_geovars
end type crtm_conf
 
INTERFACE calculate_aero_layer_factor
 
   MODULE PROCEDURE calculate_aero_layer_factor_atm_profile,&
        &calculate_aero_layer_factor_atm
 
END INTERFACE
 
INTERFACE qsmith
 
   MODULE PROCEDURE qsmith_atm,qsmith_profiles
 
END INTERFACE qsmith
 
 
 ! Add more UFO_Absorbers as needed
 ! Note: must have same ordering as CRTM_Absorbers, 
 !       CRTM_Absorber_Id, and CRTM_Absorber_Units
 character(len=MAXVARLEN), parameter :: &
      ufo_absorbers(3) = &
         [ var_mixr, var_co2, var_oz ]
 
 ! copy of ABSORBER_ID_NAME defined in CRTM_Atmosphere_Define
 character(len=MAXVARLEN), parameter :: &
      crtm_absorbers(n_valid_absorber_ids) = &
         absorber_id_name(1:n_valid_absorber_ids)
 integer, parameter :: &
      crtm_absorber_id(n_valid_absorber_ids) = &
         [ h2o_id,  co2_id,   o3_id,  n2o_id, &
            co_id,  ch4_id,   o2_id,   no_id, &
           so2_id,  no2_id,  nh3_id, hno3_id, &
            oh_id,   hf_id,  hcl_id,  hbr_id, &
            hi_id,  clo_id,  ocs_id, h2co_id, &
          hocl_id,   n2_id,  hcn_id, ch3l_id, &
          h2o2_id, c2h2_id, c2h6_id,  ph3_id, &
          cof2_id,  sf6_id,  h2s_id,hcooh_id ]
 integer, parameter :: &
      crtm_absorber_units(3) = [ &
         mass_mixing_ratio_units   & !H2O
       , volume_mixing_ratio_units & !CO2
       , volume_mixing_ratio_units & !O3
        ]
 
 character(len=MAXVARLEN), parameter :: &
      ufo_clouds(n_valid_cloud_categories,2) = &
         reshape( &
            [ var_clw,    var_cli,    var_clr,    var_cls,    var_clg,    var_clh, &
              var_clwefr, var_cliefr, var_clrefr, var_clsefr, var_clgefr, var_clhefr ] &
            , [n_valid_cloud_categories,2] )
 
 ! copy of CLOUD_CATEGORY_NAME defined in CRTM_Cloud_Define
 character(len=MAXVARLEN), parameter :: &
      crtm_clouds(n_valid_cloud_categories) = &
         cloud_category_name(1:n_valid_cloud_categories)
 integer, parameter :: &
      crtm_cloud_id(n_valid_cloud_categories) = &
         [   water_cloud, &
               ice_cloud, &
              rain_cloud, &
              snow_cloud, &
           graupel_cloud, &
              hail_cloud  ]
 
! Surface Variables
 
 character(len=MAXVARLEN), parameter :: &
      ufo_surfaces(4) = &
         [ var_sfc_wtmp, var_sfc_wspeed, var_sfc_wdir, var_sfc_sss]
 
 character(len=MAXVARLEN), parameter :: & 
      crtm_surfaces(4) = &
         [  character(len=maxvarlen):: 'Water_Temperature', 'Wind_Speed', 'Wind_Direction', 'Salinity' ]
 
 character(len=MAXVARLEN), parameter :: &
      validsurfacewindgeovars(2) = [character(len=maxvarlen) :: 'vector', 'uv']
 
contains
 
! ------------------------------------------------------------------------------
 
subroutine crtm_conf_setup(conf, f_confOpts, f_confOper)
 
implicit none
type(crtm_conf),            intent(inout) :: conf
type(fckit_configuration),  intent(in)    :: f_confopts
type(fckit_configuration),  intent(in)    :: f_confoper
 
character(*), PARAMETER :: routine_name = 'crtm_conf_setup'
character(len=255) :: irwatercoeff, viswatercoeff, &
                      irvislandcoeff, irvissnowcoeff, irvisicecoeff, &
                      mwwatercoeff
integer :: jspec, ivar
character(len=max_string) :: message
character(len=:), allocatable :: str
character(len=:), allocatable :: str_array(:)
 
CHARACTER(len=MAXVARLEN), ALLOCATABLE :: var_aerosols(:)
 
 
 !Some config needs to come from user
 !-----------------------------------
 
 !Number of sensors, each call to CRTM will be for a single sensor
 !type (zenith/scan angle will be different)
 conf%n_Sensors = 1
 
 ! absorbers, clouds and aerosols (should match what model will provide)
 
 ! Absorbers
 !----------
 conf%n_Absorbers = 0
 if (f_confoper%has("Absorbers")) &
   conf%n_Absorbers = conf%n_Absorbers + f_confoper%get_size("Absorbers")
 
 allocate( conf%Absorbers     ( conf%n_Absorbers ), &
           conf%Absorber_Id   ( conf%n_Absorbers ), &
           conf%Absorber_Units( conf%n_Absorbers ) )
 
 if (conf%n_Absorbers > 0) then
   call f_confoper%get_or_die("Absorbers",str_array)
   conf%Absorbers(1:conf%n_Absorbers) = str_array
 end if
 
 ! check for duplications
 do jspec = 2, conf%n_Absorbers
   if ( any(conf%Absorbers(jspec-1) == conf%Absorbers(jspec:conf%n_Absorbers)) ) then
     write(message,*) trim(routine_name),' error: ',trim(conf%Absorbers(jspec)),' is duplicated in Absorbers'
     call abor1_ftn(message)
   end if
 end do
 
 ! convert from CRTM names to UFO CF names and define Id and Units
 do jspec = 1, conf%n_Absorbers
   ivar = ufo_vars_getindex(crtm_absorbers, conf%Absorbers(jspec))
   if (ivar < 1 .or. ivar > size(ufo_absorbers)) then
     write(message,*) trim(routine_name),' error: ',trim(conf%Absorbers(jspec)),' not supported by UFO_Absorbers'
     call abor1_ftn(message)
   end if
   conf%Absorbers(jspec) = ufo_absorbers(ivar)
   conf%Absorber_Id(jspec) = crtm_absorber_id(ivar)
   conf%Absorber_Units(jspec) = crtm_absorber_units(ivar)
 end do
 
 
 ! Clouds
 !-------
 conf%n_Clouds = 0
 if (f_confoper%has("Clouds")) &
   conf%n_Clouds = f_confoper%get_size("Clouds")
 allocate( conf%Clouds  ( conf%n_Clouds,2), &
           conf%Cloud_Id( conf%n_Clouds ) )
 if (conf%n_Clouds > 0) then
   call f_confoper%get_or_die("Clouds",str_array) 
   conf%Clouds(1:conf%n_Clouds,1) = str_array
 
   if (f_confoper%has("Cloud_Fraction")) then
     call f_confoper%get_or_die("Cloud_Fraction",conf%Cloud_Fraction)
     if ( conf%Cloud_Fraction < 0.0 .or. &
          conf%Cloud_Fraction > 1.0 ) then
       write(message,*) trim(routine_name),' error: must specify ' // &
                        ' 0.0 <= Cloud_Fraction <= 1.0' // &
                        ' or remove Cloud_Fraction from conf' // &
                        ' and provide as a geoval'
       call abor1_ftn(message)
     end if
   else
     message = trim(routine_name) // &
             ': Cloud_Fraction is not provided in conf.' // &
             ' Will request as a geoval.'
     CALL display_message(routine_name, trim(message), warning )
   end if
 end if
 
 ! check for duplications
 do jspec = 2, conf%n_Clouds
   if ( any(conf%Clouds(jspec-1,1) == conf%Clouds(jspec:conf%n_Clouds,1)) ) then
     write(message,*) trim(routine_name),' error: ',trim(conf%Clouds(jspec,1)), &
                      ' is duplicated in Clouds'
     call abor1_ftn(message)
   end if
 end do
 
 ! convert from CRTM names to UFO CF names and define Id
 do jspec = 1, conf%n_Clouds
   ivar = ufo_vars_getindex(crtm_clouds, conf%Clouds(jspec,1))
   if (ivar < 1 .or. ivar > size(ufo_clouds)) then
     write(message,*) trim(routine_name),' error: ',trim(conf%Clouds(jspec,1)),' not supported by UFO_Clouds'
     call abor1_ftn(message)
   end if
   conf%Clouds(jspec,1:2) = ufo_clouds(ivar,1:2)
   conf%Cloud_Id(jspec)   = crtm_cloud_id(ivar)
 end do
 
 ! Aerosols
 !---------
 IF (f_confopts%has("AerosolOption")) THEN
    call f_confopts%get_or_die("AerosolOption",str)
    conf%aerosol_option = str
    conf%aerosol_option = upper2lower(conf%aerosol_option)
    CALL assign_aerosol_names(conf%aerosol_option,var_aerosols)
    conf%n_Aerosols=SIZE(var_aerosols)
 ELSE
    conf%n_Aerosols  = 0
    conf%aerosol_option = ""
 ENDIF
 
 ! Surface variables
 !----------
 conf%n_Surfaces = 0
 if (f_confoper%has("Surfaces")) &
   conf%n_Surfaces = conf%n_Surfaces + f_confoper%get_size("Surfaces")
 
 allocate( conf%Surfaces    ( conf%n_Surfaces ))
 
 if (conf%n_Surfaces > 0) then
   call f_confoper%get_or_die("Surfaces",str_array)
   conf%Surfaces(1:conf%n_Surfaces) = str_array
 end if
 
 ! check for duplications
 do jspec = 2, conf%n_Surfaces
   if ( any(conf%Surfaces(jspec-1) == conf%Surfaces(jspec:conf%n_Surfaces)) ) then
     write(message,*) 'crtm_conf_setup error: ',trim(conf%Surfaces(jspec)),' is duplicated in Surfaces'
     call abor1_ftn(message)
   end if
 end do
 
 ! convert from CRTM names to UFO CF names and define Id and Units
 do jspec = 1, conf%n_Surfaces
   ivar = ufo_vars_getindex(crtm_surfaces, conf%Surfaces(jspec))
   if (ivar < 1 .or. ivar > size(ufo_surfaces)) then
     write(message,*) 'crtm_conf_setup error: ',trim(conf%Surfaces(jspec)),' not supported by UFO_Surfaces'
     call abor1_ftn(message)
   end if
   conf%Surfaces(jspec) = ufo_surfaces(ivar)
 
 end do
 
 ! select between two surface wind geovals options
 ! valid options: vector [default], uv
 if (f_confoper%get('SurfaceWindGeoVars', str)) then
   conf%sfc_wind_geovars = str
 else
   conf%sfc_wind_geovars = 'vector'
 endif
 if (ufo_vars_getindex(validsurfacewindgeovars, conf%sfc_wind_geovars) < 1) then
   write(message,*) 'crtm_conf_setup error: invalid SurfaceWindGeoVars ',trim(str)
   call abor1_ftn(message)
 end if
 
 ! Sea_Surface_Salinity
 !---------
 IF (f_confopts%get("Salinity",str)) THEN
    conf%salinity_option = str
 ELSE
    conf%salinity_option = 'off'
 END IF
 
 !Allocate SENSOR_ID
 allocate(conf%SENSOR_ID(conf%n_Sensors))
 
 !Get sensor ID from config
 call f_confopts%get_or_die("Sensor_ID",str)
 conf%SENSOR_ID(conf%n_Sensors) = str
 
 !ENDIAN type
 call f_confopts%get_or_die("EndianType",str)
 conf%ENDIAN_TYPE = str
 
 !Path to coefficient files
 call f_confopts%get_or_die("CoefficientPath",str)
 conf%COEFFICIENT_PATH = str
 
 ! Coefficient file prefixes
 irwatercoeff = "Nalli"
 if (f_confopts%has("IRwaterCoeff")) then
    call f_confopts%get_or_die("IRwaterCoeff",str)
    irwatercoeff = str
 end if
 viswatercoeff = "NPOESS"
 if (f_confopts%has("VISwaterCoeff")) then
    call f_confopts%get_or_die("VISwaterCoeff",str)
    viswatercoeff = str
 end if
 irvislandcoeff = "NPOESS"
 if (f_confopts%has("IRVISlandCoeff")) then
    call f_confopts%get_or_die("IRVISlandCoeff",str)
    irvislandcoeff = str
 end if
 irvissnowcoeff = "NPOESS"
 if (f_confopts%has("IRVISsnowCoeff")) then
    call f_confopts%get_or_die("IRVISsnowCoeff",str)
    irvissnowcoeff = str
 end if
 irvisicecoeff = "NPOESS"
 if (f_confopts%has("IRVISiceCoeff")) then
    call f_confopts%get_or_die("IRVISiceCoeff",str)
    irvisicecoeff = str
 end if
 mwwatercoeff = "FASTEM6"
 if (f_confopts%has("MWwaterCoeff")) then
    call f_confopts%get_or_die("MWwaterCoeff",str)
    mwwatercoeff = str
 end if
 
 ! Define water, snow, ice (WSI) categories
 select case (trim(irvislandcoeff))
    case ("USGS")
       allocate(conf%Land_WSI(2))
       conf%Land_WSI(1:2) = (/16,24/)
    case ("IGBP")
       allocate(conf%Land_WSI(2))
       conf%Land_WSI(1:2) = (/15,17/)
    case default
       allocate(conf%Land_WSI(1))
       conf%Land_WSI(1) = -1
 end select
 
 ! IR emissivity coeff files
 conf%IRwaterCoeff_File = trim(irwatercoeff)//".IRwater.EmisCoeff.bin"
 conf%IRlandCoeff_File  = trim(irvislandcoeff)//".IRland.EmisCoeff.bin"
 conf%IRsnowCoeff_File  = trim(irvissnowcoeff)//".IRsnow.EmisCoeff.bin"
 conf%IRiceCoeff_File   = trim(irvisicecoeff)//".IRice.EmisCoeff.bin"
 
 !VIS emissivity coeff files
 conf%VISwaterCoeff_File = trim(viswatercoeff)//".VISwater.EmisCoeff.bin"
 conf%VISlandCoeff_File  = trim(irvislandcoeff)//".VISland.EmisCoeff.bin"
 conf%VISsnowCoeff_File  = trim(irvissnowcoeff)//".VISsnow.EmisCoeff.bin"
 conf%VISiceCoeff_File   = trim(irvisicecoeff)//".VISice.EmisCoeff.bin"
 
 ! MW water emissivity coeff file
 conf%MWwaterCoeff_File = trim(mwwatercoeff)//".MWwater.EmisCoeff.bin"
 
 conf%inspect = 0
 if (f_confopts%has("InspectProfileNumber")) then
   call f_confopts%get_or_die("InspectProfileNumber",conf%inspect)
 endif
 
end subroutine crtm_conf_setup
 
! -----------------------------------------------------------------------------
 
subroutine crtm_conf_delete(conf)
 
implicit none
type(crtm_conf), intent(inout) :: conf
 
 deallocate(conf%SENSOR_ID)
 deallocate(conf%Land_WSI)
 deallocate(conf%Absorbers)
 deallocate(conf%Absorber_Id)
 deallocate(conf%Absorber_Units)
 deallocate(conf%Clouds)
 deallocate(conf%Cloud_Id)
 
end subroutine crtm_conf_delete
 
! ------------------------------------------------------------------------------
 
subroutine crtm_comm_stat_check(stat, PROGRAM_NAME, message, f_comm)
use fckit_mpi_module,   only: fckit_mpi_comm
 
implicit none
 
integer,              intent(in) :: stat
character(*),         intent(in) :: program_name
character(*),         intent(in) :: message
type(fckit_mpi_comm), intent(in) :: f_comm
 
character(max_string) :: rank_message
 
 if ( stat /= success ) THEN
    write(rank_message,*) trim(message)," on rank ",f_comm%rank()
    call display_message( program_name, rank_message, failure )
    call abor1_ftn("Abort from "//program_name)
 end if
 
end subroutine crtm_comm_stat_check
 
! ------------------------------------------------------------------------------
 
subroutine ufo_crtm_skip_profiles(n_Profiles,n_Channels,channels,obss,Skip_Profiles)
! Profiles are skipped when the ObsValue of all channels is missing.
! TODO: Use complete QC information
! It would be more comprehensive to use EffectiveQC or EffectiveError. That
! would require those ObsSpace values to be initialized before calls to
! this subroutine within ufo_radiancecrtm_simobs+ufo_radiancecrtm_tlad_settraj.
use missing_values_mod
 
implicit none
integer,              intent(in)    :: n_profiles, n_channels
type(c_ptr), value,   intent(in)    :: obss
integer(c_int),       intent(in)    :: channels(:)
logical,              intent(inout) :: skip_profiles(:)
 
integer :: jprofile, jchannel
character(len=MAXVARLEN) :: varname
real(kind_real)  :: obsval(n_profiles,n_channels)
!real(kind_real) :: EffObsErr(n_Profiles,n_Channels)
!integer         :: EffQC(n_Profiles,n_Channels)
 
real(c_double) :: missing
 
 ! Set missing value
 missing = missing_value(missing)
 
 obsval = missing
! EffObsErr = missing
! EffQC = 0
 
 do jchannel = 1, n_channels
   call get_var_name(channels(jchannel),varname)
   call obsspace_get_db(obss, "ObsValue", varname, obsval(:,jchannel))
!   call obsspace_get_db(obss, "EffectiveError", varname, EffObsErr(:,jchannel))
!   call obsspace_get_db(obss, "EffectiveQC{iter}", varname, EffQC(:,jchannel))
 enddo
 
 !Loop over all n_Profiles, i.e. number of locations
 do jprofile = 1, n_profiles
   skip_profiles(jprofile) = all(obsval(jprofile,:) == missing)
!                       .OR. all(EffObsErr(jprofile,:) == missing) &
!                       .OR. all(EffQC(jprofile,:) /= 0)
 end do
 
end subroutine ufo_crtm_skip_profiles
 
! ------------------------------------------------------------------------------
 
SUBROUTINE load_atm_data(n_Profiles, n_Layers, geovals, atm, conf)
implicit none
 
integer, intent(in) :: n_profiles, n_layers
type(ufo_geovals), intent(in) :: geovals
type(crtm_atmosphere_type), intent(inout) :: atm(:)
type(crtm_conf) :: conf
 
! Local variables
integer :: k1, jspec
type(ufo_geoval), pointer :: geoval
character(max_string) :: err_msg
 
  ! Populate the atmosphere structures for CRTM
  ! -------------------------------------------
 
  call ufo_geovals_get_var(geovals, var_ts, geoval)
  ! Check model levels is consistent in geovals & crtm
  if (geoval%nval /= n_layers) then
    write(err_msg,*) 'Load_Atm_Data error: layers inconsistent!'
    call abor1_ftn(err_msg)
  endif
 
  do k1 = 1, n_profiles
    atm(k1)%Temperature(1:n_layers) = geoval%vals(:, k1)
  end do
 
  call ufo_geovals_get_var(geovals, var_prs, geoval)
  do k1 = 1, n_profiles
    atm(k1)%Pressure(1:n_layers) = geoval%vals(:, k1) * 0.01  ! to hPa
  end do
 
  call ufo_geovals_get_var(geovals, var_prsi, geoval)
  do k1 = 1, n_profiles
    atm(k1)%Level_Pressure(:) = geoval%vals(:, k1) * 0.01     ! to hPa
    atm(k1)%Climatology = us_standard_atmosphere
  end do
 
  do jspec = 1, conf%n_Absorbers
    ! O3 Absorber has special treatment for Aerosols
    if (cmp_strings(conf%Absorbers(jspec), var_oz) .AND. &
      ufo_vars_getindex(geovals%variables, var_oz) < 0 .AND. &
      (.NOT. cmp_strings(conf%aerosol_option,""))) then
      do k1 = 1, n_profiles
        atm(k1)%Absorber(1:n_layers, jspec) = ozone_default_value
      end do
    else
      call ufo_geovals_get_var(geovals, conf%Absorbers(jspec), geoval)
      do k1 = 1, n_profiles
        atm(k1)%Absorber(1:n_layers, jspec) = geoval%vals(:, k1)
      end do
    end if
    do k1 = 1, n_profiles
      atm(k1)%Absorber_Id(jspec) = conf%Absorber_Id(jspec)
      atm(k1)%Absorber_Units(jspec) = conf%Absorber_Units(jspec)
    end do
  end do
 
  do jspec = 1, conf%n_Clouds
    ! cloud species content
    CALL ufo_geovals_get_var(geovals, conf%Clouds(jspec,1), geoval)
    do k1 = 1, n_profiles
      atm(k1)%Cloud(jspec)%Water_Content = geoval%vals(:, k1)
      atm(k1)%Cloud(jspec)%Type = conf%Cloud_Id(jspec)
    end do
 
    ! effective radius
    CALL ufo_geovals_get_var(geovals, conf%Clouds(jspec,2), geoval)
    do k1 = 1, n_profiles
      atm(k1)%Cloud(jspec)%Effective_Radius = geoval%vals(:, k1)
    end do
  end do
 
  ! When n_Clouds>0, Cloud_Fraction must either be provided as geoval or in conf
  ! If Cloud_Fraction is provided in conf,then use the Cloud_Fraction in conf
  ! If Cloud_Fraction is not provided in conf , then use the Cloud_Fraction in geoval 
  ! and make sure the cloud fraction interpolated from background is in the valid range [0.0,1.0]
  if (conf%n_Clouds > 0) then
    if ( conf%Cloud_Fraction >= 0.0  ) then 
      do k1 = 1, n_profiles
        atm(k1)%Cloud_Fraction(:) = conf%Cloud_Fraction
      end do
    else
      if ( ufo_vars_getindex(geovals%variables, var_cldfrac) > 0 ) then
        CALL ufo_geovals_get_var(geovals, var_cldfrac, geoval)
        do k1 = 1, n_profiles
          where( geoval%vals(:, k1) < 0_kind_real ) geoval%vals(:, k1) = 0_kind_real
          where( geoval%vals(:, k1) > 1_kind_real ) geoval%vals(:, k1) = 1_kind_real
          atm(k1)%Cloud_Fraction(:) =  geoval%vals(:, k1)  
        end do
      end if
    end if
  end if
 
end subroutine load_atm_data
 
! ------------------------------------------------------------------------------
 
subroutine load_sfc_data(n_Profiles, n_Channels, channels, geovals, sfc, chinfo, obss, conf)
 
implicit none
 
integer,                     intent(in)    :: n_profiles, n_channels
type(ufo_geovals),           intent(in)    :: geovals
type(crtm_surface_type),     intent(inout) :: sfc(:)
type(crtm_channelinfo_type), intent(in)    :: chinfo(:)
type(c_ptr), value,          intent(in)    :: obss
integer(c_int),              intent(in)    :: channels(:)
type(crtm_conf),             intent(in)    :: conf
 
type(ufo_geoval), pointer :: geoval, u, v
integer :: k1, n1
integer :: iland
 
! Surface type definitions for default SfcOptics definitions
! for IR and VIS, this is the NPOESS reflectivities.
integer, parameter :: tundra_surface_type         = 10  ! NPOESS Land surface type for IR/VIS Land SfcOptics
integer, parameter :: scrub_surface_type          =  7  ! NPOESS Land surface type for IR/VIS Land SfcOptics
integer, parameter :: coarse_soil_type            =  1  ! Soil type                for MW land SfcOptics
integer, parameter :: groundcover_vegetation_type =  7  ! Vegetation type          for MW Land SfcOptics
integer, parameter :: bare_soil_vegetation_type   = 11  ! Vegetation type          for MW Land SfcOptics
integer, parameter :: sea_water_type              =  1  ! Water type               for all SfcOptics
integer, parameter :: fresh_snow_type             =  2  ! NPOESS Snow type         for IR/VIS SfcOptics
integer, parameter :: fresh_ice_type              =  1  ! NPOESS Ice type          for IR/VIS SfcOptics
 
character(len=MAXVARLEN) :: varname
 
real(kind_real), allocatable :: obstb(:,:)
 
  allocate(obstb(n_profiles, n_channels))
  obstb = 0.0_kind_real
 
  do n1 = 1, n_channels
    call get_var_name(channels(n1), varname)
    call obsspace_get_db(obss, "ObsValue", varname, obstb(:, n1))
  enddo
 
  do k1 = 1, n_profiles
    !Pass sensor information
    sfc(k1)%sensordata%sensor_id        = chinfo(1)%sensor_id
    sfc(k1)%sensordata%wmo_sensor_id    = chinfo(1)%wmo_sensor_id
    sfc(k1)%sensordata%wmo_satellite_id = chinfo(1)%wmo_satellite_id
    sfc(k1)%sensordata%sensor_channel   = channels
 
    !Pass observation value
    do n1 = 1, n_channels
      sfc(k1)%sensordata%tb(n1) = obstb(k1, n1)
    enddo
 
    !Water_type
    sfc(k1)%Water_Type = sea_water_type    !** NOTE: need to check how to determine fresh vs sea water types (salinity???)
  end do
  deallocate(obstb)
 
  if (ufo_vars_getindex(geovals%variables, var_sfc_wspeed) > 0 .and. &
      ufo_vars_getindex(geovals%variables, var_sfc_wdir) > 0) then
    ! Directly use model-provided wind speed and direction
    !Wind_Speed
    call ufo_geovals_get_var(geovals, var_sfc_wspeed, geoval)
    do k1 = 1, n_profiles
      sfc(k1)%Wind_Speed = geoval%vals(1, k1)
    end do
 
    !Wind_Direction
    call ufo_geovals_get_var(geovals, var_sfc_wdir, geoval)
    do k1 = 1, n_profiles
      sfc(k1)%Wind_Direction = geoval%vals(1, k1)
    end do
  else if (ufo_vars_getindex(geovals%variables, var_sfc_u) > 0 .and. &
      ufo_vars_getindex(geovals%variables, var_sfc_v) > 0) then
    ! Convert 2d wind components to speed and direction
    call ufo_geovals_get_var(geovals, var_sfc_u, u)
    call ufo_geovals_get_var(geovals, var_sfc_v, v)
 
    !Wind_Speed
    do k1 = 1, n_profiles
      sfc(k1)%Wind_Speed = sqrt(u%vals(1, k1)**2 + v%vals(1, k1)**2)
    end do
 
    !Wind_Direction
    do k1 = 1, n_profiles
      sfc(k1)%Wind_Direction = uv_to_wdir(u%vals(1, k1), v%vals(1, k1))
    end do
  else
    call abor1_ftn('Load_Sfc_Data error: missing surface wind geovals')
  end if
 
  !Water_Coverage
  call ufo_geovals_get_var(geovals, var_sfc_wfrac, geoval)
  do k1 = 1, n_profiles
    sfc(k1)%Water_Coverage = geoval%vals(1, k1)
  end do
 
  !Water_Temperature
  call ufo_geovals_get_var(geovals, var_sfc_wtmp, geoval)
  do k1 = 1, n_profiles
    sfc(k1)%Water_Temperature = geoval%vals(1, k1)
  end do
 
  !Ice_Coverage
  call ufo_geovals_get_var(geovals, var_sfc_ifrac, geoval)
  do k1 = 1, n_profiles
    sfc(k1)%Ice_Coverage = geoval%vals(1, k1)
  end do
 
  !Ice_Temperature
  call ufo_geovals_get_var(geovals, var_sfc_itmp, geoval)
  do k1 = 1, n_profiles
    sfc(k1)%Ice_Temperature = geoval%vals(1, k1)
  end do
 
  !Snow_Coverage
  call ufo_geovals_get_var(geovals, var_sfc_sfrac, geoval)
  do k1 = 1, n_profiles
    sfc(k1)%Snow_Coverage = geoval%vals(1, k1)
  end do
 
  !Snow_Temperature
  call ufo_geovals_get_var(geovals, var_sfc_stmp, geoval)
  do k1 = 1, n_profiles
    sfc(k1)%Snow_Temperature = geoval%vals(1, k1)
  end do
 
  !Snow_Depth
  call ufo_geovals_get_var(geovals, var_sfc_sdepth, geoval)
  do k1 = 1, n_profiles
    sfc(k1)%Snow_Depth = geoval%vals(1, k1)
  end do
 
  !Land_Coverage
  call ufo_geovals_get_var(geovals, var_sfc_lfrac, geoval)
  do k1 = 1, n_profiles
    sfc(k1)%Land_Coverage = geoval%vals(1, k1)
  end do
 
  !Land_Type
  ! + used to lookup land sfc emiss. for IR and VIS
  ! + land sfc emiss. undefined over water/snow/ice
  call ufo_geovals_get_var(geovals, var_sfc_landtyp, geoval)
  do k1 = 1, n_profiles
    iland = int(geoval%vals(1, k1))
    if (.not.any(iland == conf%Land_WSI)) then
      sfc(k1)%Land_Type = iland
    end if
  end do
 
  !Land_Temperature
  call ufo_geovals_get_var(geovals, var_sfc_ltmp, geoval)
  do k1 = 1, n_profiles
    sfc(k1)%Land_Temperature = geoval%vals(1, k1)
  end do
 
  !Lai
  call ufo_geovals_get_var(geovals, var_sfc_lai, geoval)
  do k1 = 1, n_profiles
    sfc(k1)%Lai = geoval%vals(1, k1)
  end do
 
  !Vegetation_Fraction
  call ufo_geovals_get_var(geovals, var_sfc_vegfrac, geoval)
  do k1 = 1, n_profiles
    sfc(k1)%Vegetation_Fraction = geoval%vals(1, k1)
  end do
 
  !Vegetation_Type
  call ufo_geovals_get_var(geovals, var_sfc_vegtyp, geoval)
  do k1 = 1, n_profiles
    sfc(k1)%Vegetation_Type = int(geoval%vals(1, k1))
  end do
 
  !Soil_Type
  call ufo_geovals_get_var(geovals, var_sfc_soiltyp, geoval)
  do k1 = 1, n_profiles
    sfc(k1)%Soil_Type = int(geoval%vals(1, k1))
  end do
 
  !Soil_Moisture_Content
  call ufo_geovals_get_var(geovals, var_sfc_soilm, geoval)
  do k1 = 1, n_profiles
    sfc(k1)%Soil_Moisture_Content = geoval%vals(1, k1)
  end do
 
  !Soil_Temperature
  call ufo_geovals_get_var(geovals, var_sfc_soilt, geoval)
  do k1 = 1, n_profiles
    sfc(k1)%Soil_Temperature = geoval%vals(1, k1)
  end do
  
  !Sea_Surface_Salinity
  if (cmp_strings(conf%salinity_option, "on")) THEN
    call ufo_geovals_get_var(geovals, var_sfc_sss, geoval)
    do k1 = 1, n_profiles
      sfc(k1)%Salinity = geoval%vals(1, k1)
    end do
  end if
 
  ! Update Ice_Coverage and Land_Coverage for
  ! Soil_Type or Vegetation_Type corresponding to Glacial land ice
  do k1 = 1, n_profiles
    if (sfc(k1)%Land_Coverage > zero .and. &
       (sfc(k1)%Soil_Type == 9 .or. sfc(k1)%Vegetation_Type == 13)) then
      sfc(k1)%Ice_Coverage = min(sfc(k1)%Ice_Coverage + sfc(k1)%Land_Coverage, one)
      sfc(k1)%Land_Coverage = zero
    end if
  end do
 
end subroutine load_sfc_data
 
! ------------------------------------------------------------------------------
 
subroutine load_geom_data(obss,geo,geo_hf)
 
implicit none
type(c_ptr), value,       intent(in)    :: obss
type(crtm_geometry_type), intent(inout) :: geo(:)
type(crtm_geometry_type), intent(inout), optional :: geo_hf(:)
real(kind_real), allocatable :: tmpvar(:)
integer :: nlocs
character(kind=c_char,len=101) :: obsname
 
 call obsspace_obsname(obss, obsname)
 nlocs = obsspace_get_nlocs(obss)
 allocate(tmpvar(nlocs))
 
 call obsspace_get_db(obss, "MetaData", "sensor_zenith_angle", tmpvar)
 geo(:)%Sensor_Zenith_Angle = abs(tmpvar(:)) ! needs to be absolute value
 
 call obsspace_get_db(obss, "MetaData", "solar_zenith_angle", tmpvar)
 geo(:)%Source_Zenith_Angle = tmpvar(:)
 
 call obsspace_get_db(obss, "MetaData", "sensor_azimuth_angle", tmpvar)
 geo(:)%Sensor_Azimuth_Angle = tmpvar(:)
 
 call obsspace_get_db(obss, "MetaData", "solar_azimuth_angle", tmpvar)
 geo(:)%Source_Azimuth_Angle = tmpvar(:)
 
!  For some microwave instruments the solar and sensor azimuth angles can be
!  missing  (given a value of 10^11).  Set these to zero to get past CRTM QC.
 where (geo(:)%Source_Azimuth_Angle < 0.0_kind_real .or. &
        geo(:)%Source_Azimuth_Angle > 360.0_kind_real) &
    geo(:)%Source_Azimuth_Angle = 0.0_kind_real
 where (geo(:)%Sensor_Azimuth_Angle < 0.0_kind_real .or. &
        geo(:)%Sensor_Azimuth_Angle > 360.0_kind_real) &
    geo(:)%Sensor_Azimuth_Angle = 0.0_kind_real
 
 where (abs(geo(:)%Source_Zenith_Angle) > 180.0_kind_real) &
    geo(:)%Source_Zenith_Angle = 100.0_kind_real
 
 call obsspace_get_db(obss, "MetaData", "scan_position", tmpvar)
 geo(:)%Ifov = tmpvar(:)
 
 call obsspace_get_db(obss, "MetaData", "sensor_view_angle", tmpvar) !The Sensor_Scan_Angle is optional
 geo(:)%Sensor_Scan_Angle = tmpvar(:)
 
 where (abs(geo(:)%Sensor_Scan_Angle) > 80.0_kind_real) &
    geo(:)%Sensor_Scan_Angle = 0.0_kind_real
 
!read geophysical values for gmi high frequency channels 10-13.
 if (cmp_strings(trim(obsname),'GMI-GPM') .or. cmp_strings(trim(obsname),'gmi_gpm')) then
    if ( present(geo_hf) ) then
       geo_hf = geo
       if (obsspace_has(obss, "MetaData", "sensor_zenith_angle1")) then
          call obsspace_get_db(obss, "MetaData", "sensor_zenith_angle1", tmpvar)
          geo_hf(:)%Sensor_Zenith_Angle = abs(tmpvar(:)) ! needs to be absolute value
       endif
       if (obsspace_has(obss, "MetaData", "solar_zenith_angle1")) then
          call obsspace_get_db(obss, "MetaData", "solar_zenith_angle1", tmpvar)
          geo_hf(:)%Source_Zenith_Angle = tmpvar(:)
       endif
       if (obsspace_has(obss, "MetaData", "sensor_azimuth_angle1")) then
          call obsspace_get_db(obss, "MetaData", "sensor_azimuth_angle1", tmpvar)
          geo_hf(:)%Sensor_Azimuth_Angle = tmpvar(:)
       endif
       if (obsspace_has(obss, "MetaData", "solar_azimuth_angle1")) then
          call obsspace_get_db(obss, "MetaData", "solar_azimuth_angle1", tmpvar)
          geo_hf(:)%Source_Azimuth_Angle = tmpvar(:)
       endif
       if (obsspace_has(obss, "MetaData", "sensor_view_angle1")) then
          call obsspace_get_db(obss, "MetaData", "sensor_view_angle1", tmpvar)
          geo_hf(:)%Sensor_Scan_Angle = tmpvar(:)
       endif
    endif
 endif
 
 deallocate(tmpvar)
 
end subroutine load_geom_data
 
! ------------------------------------------------------------------------------
 
subroutine get_var_name(n,varname)
 
integer, intent(in) :: n
character(len=*), intent(out) :: varname
 
character(len=6) :: chan
 
 write(chan, '(I0)') n
 varname = 'brightness_temperature_' // trim(chan)
 
end subroutine get_var_name
 
! -----------------------------------------------------------------------------
 
!> \brief Determines the wind direction from U and V components
!!
!! \details **uv_to_wdir** Calculates the wind direction, as measured clockwise
!! from north, similar to an azimuth angle.  Takes the eastward and northward
!! wind component magnitudes, respectively, as arguments.
!! Due to the azimuthal convention used here, the inverse equations are:
!! u = w * cos(wdir * deg2rad)
!! v = w * sin(wdir * deg2rad)
!! where w is the wind speed
function uv_to_wdir(u, v) result(wdir)
 
use ufo_constants_mod, only: zero, one, two, pi, rad2deg
 
implicit none
 
real(kind=kind_real), intent(in) :: u !< eastward_wind
real(kind=kind_real), intent(in) :: v !< northward_wind
real(kind=kind_real) :: wdir
real(kind=kind_real) :: windratio, windangle
integer :: iquadrant
real(kind=kind_real),parameter:: windscale = 999999.0_kind_real
real(kind=kind_real),parameter:: windlimit = 0.0001_kind_real
real(kind=kind_real),parameter:: quadcof(4,2) = &
  reshape((/zero,  one,  one,  two, &
            one,  -one,  one, -one/), (/4,2/))
 
  if (u >= zero .and. v >= zero) iquadrant = 1
  if (u >= zero .and. v <  zero) iquadrant = 2
  if (u <  zero .and. v >= zero) iquadrant = 4
  if (u <  zero .and. v <  zero) iquadrant = 3
 
  if (abs(v) >= windlimit) then
    windratio = u / v
  else
    windratio = zero
    if (abs(u) > windlimit) then
      windratio = windscale * u
    endif
  endif
  windangle = atan(abs(windratio))   ! wind azimuth is in radians
  wdir = ( quadcof(iquadrant, 1) * pi + windangle * quadcof(iquadrant, 2) ) * rad2deg
 
end function uv_to_wdir
 
! -----------------------------------------------------------------------------
 
SUBROUTINE load_aerosol_data(n_profiles,n_layers,geovals,&
     &aerosol_option,atm)
 
    USE crtm_aerosolcoeff, ONLY: aeroc
 
    INTEGER, INTENT(in) :: n_profiles,n_layers
    TYPE(ufo_geovals), INTENT(in) :: geovals
    TYPE(crtm_atmosphere_type), INTENT(inout) :: atm(:)
 
    CHARACTER(*) :: aerosol_option
    CHARACTER(max_string) :: message
    CHARACTER(len=MAXVARLEN) :: varname
 
    TYPE(ufo_geoval), POINTER :: geoval
 
    CHARACTER(*), PARAMETER :: routine_name = 'Load_Aerosol_Data'
 
    REAL(kind_real), DIMENSION(n_layers,n_profiles) :: rh
    INTEGER :: ivar
 
    IF (cmp_strings(aerosol_option, "aerosols_gocart_default")) THEN
       varname=var_rh
       CALL ufo_geovals_get_var(geovals, varname, geoval)
       rh(1:n_layers,1:n_profiles)=geoval%vals(1:n_layers,1:n_profiles)
       WHERE (rh > 1_kind_real) rh=1_kind_real
       CALL assign_gocart_default
    ELSEIF (cmp_strings(aerosol_option, "aerosols_gocart_merra_2")) THEN
       varname=var_rh
       CALL ufo_geovals_get_var(geovals, varname, geoval)
       rh(1:n_layers,1:n_profiles)=geoval%vals(1:n_layers,1:n_profiles)
       WHERE (rh > 1_kind_real) rh=1_kind_real
       CALL assign_gocart_merra_2
    ELSEIF (cmp_strings(aerosol_option, "aerosols_other")) THEN
       CALL assign_other
    ELSE
       message = 'this aerosol not implemented - check later'
       CALL display_message( aerosol_option, message, failure )
       stop
    ENDIF
 
  CONTAINS 
 
    SUBROUTINE assign_gocart_default
 
      INTEGER, PARAMETER :: ndust_bins=5, nseas_bins=4
      REAL(kind_real), DIMENSION(ndust_bins), PARAMETER  :: dust_radii=[&
           &0.55_kind_real,1.4_kind_real,2.4_kind_real,4.5_kind_real,8.0_kind_real]
      
      INTEGER, DIMENSION(nseas_bins), PARAMETER  :: seas_types=[&
           seasalt_ssam_aerosol,seasalt_sscm1_aerosol,seasalt_sscm2_aerosol,seasalt_sscm3_aerosol]
 
      REAL(kind_real), DIMENSION(n_layers) :: layer_factors
      
      INTEGER :: i,k,m
 
      CHARACTER(len=MAXVARLEN) :: varname
 
      DO m=1,n_profiles
 
         CALL calculate_aero_layer_factor(atm(m),layer_factors)
 
         DO i=1,n_aerosols_gocart_default
 
            varname=var_aerosols_gocart_default(i)
            CALL ufo_geovals_get_var(geovals,varname, geoval)
 
            atm(m)%aerosol(i)%Concentration(1:n_layers)=&
                 &max(geoval%vals(:,m)*layer_factors,aerosol_concentration_minvalue_layer)
 
            SELECT CASE (trim(varname))
            CASE (var_sulfate)
               atm(m)%aerosol(i)%type  = sulfate_aerosol
               DO k=1,n_layers
                  atm(m)%aerosol(i)%effective_radius(k)=&
                       &gocart_aerosol_size(atm(m)%aerosol(i)%type, &
                       &rh(k,m))
               ENDDO
 
            CASE (var_bcphobic)
               atm(m)%aerosol(i)%type  = black_carbon_aerosol
               atm(m)%aerosol(i)%effective_radius(:)=&
                    &aeroc%Reff(1,atm(m)%aerosol(i)%type)
            CASE (var_bcphilic)
               atm(m)%aerosol(i)%type  = black_carbon_aerosol
               DO k=1,n_layers
                  atm(m)%aerosol(i)%effective_radius(k)=&
                       &gocart_aerosol_size(atm(m)%aerosol(i)%type, &
                       &rh(k,m))
               ENDDO
 
            CASE (var_ocphobic)
               atm(m)%aerosol(i)%type  = organic_carbon_aerosol
               atm(m)%aerosol(i)%effective_radius(:)=&
                    &aeroc%Reff(1,atm(m)%aerosol(i)%type)
            CASE (var_ocphilic)
               atm(m)%aerosol(i)%type  = organic_carbon_aerosol
               DO k=1,n_layers
                  atm(m)%aerosol(i)%effective_radius(k)=&
                       &gocart_aerosol_size(atm(m)%aerosol(i)%type, &
                       &rh(k,m))
               ENDDO
 
            CASE (var_du001)
               atm(m)%aerosol(i)%type  = dust_aerosol
               atm(m)%aerosol(i)%effective_radius(:)=dust_radii(1)
            CASE (var_du002)
               atm(m)%aerosol(i)%type  = dust_aerosol
               atm(m)%aerosol(i)%effective_radius(:)=dust_radii(2)
            CASE (var_du003)
               atm(m)%aerosol(i)%type  = dust_aerosol
               atm(m)%aerosol(i)%effective_radius(:)=dust_radii(3)
            CASE (var_du004)
               atm(m)%aerosol(i)%type  = dust_aerosol
               atm(m)%aerosol(i)%effective_radius(:)=dust_radii(4)
            CASE (var_du005)
               atm(m)%aerosol(i)%type  = dust_aerosol
               atm(m)%aerosol(i)%effective_radius(:)=dust_radii(5)
 
            CASE (var_ss001)
               atm(m)%aerosol(i)%type  = seas_types(1)
               DO k=1,n_layers
                  atm(m)%aerosol(i)%effective_radius(k)=&
                       &gocart_aerosol_size(atm(m)%aerosol(i)%type, &
                       &rh(k,m))
               ENDDO
            CASE (var_ss002)
               atm(m)%aerosol(i)%type  = seas_types(2)
               DO k=1,n_layers
                  atm(m)%aerosol(i)%effective_radius(k)=&
                       &gocart_aerosol_size(atm(m)%aerosol(i)%type, &
                       &rh(k,m))
               ENDDO
            CASE (var_ss003)
               atm(m)%aerosol(i)%type  = seas_types(3)
               DO k=1,n_layers
                  atm(m)%aerosol(i)%effective_radius(k)=&
                       &gocart_aerosol_size(atm(m)%aerosol(i)%type, &
                       &rh(k,m))
               ENDDO
            CASE (var_ss004)
               atm(m)%aerosol(i)%type  = seas_types(4)
               DO k=1,n_layers
                  atm(m)%aerosol(i)%effective_radius(k)=&
                       &gocart_aerosol_size(atm(m)%aerosol(i)%type, &
                       &rh(k,m))
               ENDDO
            END SELECT
 
         ENDDO
 
      ENDDO
 
 
    END SUBROUTINE assign_gocart_default
 
    SUBROUTINE assign_gocart_merra_2
 
      message = 'this aerosol not implemented in the CRTM - check later'
      CALL display_message( aerosol_option, message, failure )
      stop
 
    END SUBROUTINE assign_gocart_merra_2
 
    SUBROUTINE assign_other
 
      message = 'this aerosol not implemented - check later'
      CALL display_message( aerosol_option, message, failure )
      stop
 
    END SUBROUTINE assign_other
    
  END SUBROUTINE load_aerosol_data
 
  SUBROUTINE assign_aerosol_names(aerosol_option,var_aerosols)
    
    CHARACTER(*), INTENT(in) :: aerosol_option
    CHARACTER(len=MAXVARLEN), ALLOCATABLE, INTENT(out) :: var_aerosols(:)
 
    CHARACTER(max_string) :: err_msg
 
    IF (aerosol_option == "aerosols_gocart_default") THEN
       ALLOCATE(var_aerosols(SIZE(var_aerosols_gocart_default)))
       var_aerosols=var_aerosols_gocart_default
    ELSEIF (aerosol_option == "aerosols_gocart_merra_2") THEN
       ALLOCATE(var_aerosols(SIZE(var_aerosols_gocart_merra_2)))
       var_aerosols=var_aerosols_gocart_merra_2
    ELSEIF (aerosol_option == "var_aerosols_other") THEN
       ALLOCATE(var_aerosols(SIZE(var_aerosols_other)))
       var_aerosols=var_aerosols_other
    ELSE
       WRITE(err_msg,*) 'assign_aerosol_names: aerosol_option not implemented '//trim(aerosol_option)
        call abor1_ftn(err_msg)
     END IF
 
   END SUBROUTINE assign_aerosol_names
 
   SUBROUTINE calculate_aero_layer_factor_atm_profile(atm,layer_factors)
 
     TYPE(crtm_atmosphere_type), INTENT(in) :: atm
     REAL(kind_real), INTENT(out) :: layer_factors(:)
 
     INTEGER :: k
 
     DO k=1,SIZE(layer_factors)
!correct for mixing ratio factor layer_factors 
!being calculated from dry pressure, cotton eq. (2.4)
!p_dry=p_total/(1+1.61*mixing_ratio)
        layer_factors(k)=1e-9_kind_real*(atm%Level_Pressure(k)-&
             &atm%Level_Pressure(k-1))*100_kind_real/grav/&
             &(1_kind_real+rv_rd*atm%Absorber(k,1)*1e-3_kind_real)
     ENDDO
 
   END SUBROUTINE calculate_aero_layer_factor_atm_profile
 
   SUBROUTINE calculate_aero_layer_factor_atm(atm,layer_factors)
 
     TYPE(crtm_atmosphere_type), INTENT(in) :: atm(:)
     REAL(kind_real), INTENT(out) :: layer_factors(:,:)
 
     INTEGER :: k,m
 
     DO k=1,SIZE(layer_factors,1)
        DO m=1,SIZE(layer_factors,2)
!correct for mixing ratio factor layer_factors 
!being calculated from dry pressure, cotton eq. (2.4)
!p_dry=p_total/(1+1.61*mixing_ratio)
           layer_factors(k,m)=1e-9_kind_real*(atm(m)%Level_Pressure(k)-&
                &atm(m)%Level_Pressure(k-1))*100_kind_real/grav/&
                &(1_kind_real+rv_rd*atm(m)%Absorber(k,1)*1.e-3_kind_real)
        ENDDO
     ENDDO
     
   END SUBROUTINE calculate_aero_layer_factor_atm
 
   FUNCTION gocart_aerosol_size( itype, rh ) & ! rh input in 0-1
        &result(r_eff)   ! in micrometer
 
     USE crtm_aerosolcoeff, ONLY: aeroc
     IMPLICIT NONE
 
!
!   modified from a function provided by quanhua liu
!
     INTEGER ,INTENT(in) :: itype
     REAL(kind_real)    ,INTENT(in) :: rh
 
     INTEGER :: j1,j2,m
     REAL(kind_real)    :: h1
     REAL(kind_real)    :: r_eff
 
     j2 = 0
     j1 = 1
     IF ( rh <= aeroc%rh(1) ) THEN
        j1 = 1
     ELSE IF ( rh >= aeroc%rh(aeroc%n_rh) ) THEN
        j1 = aeroc%n_rh
     ELSE
        DO m = 1, aeroc%n_rh-1
           IF ( rh < aeroc%rh(m+1) .AND. rh > aeroc%rh(m) ) THEN
              j1 = m
              j2 = m+1
              h1 = (rh-aeroc%rh(m))/(aeroc%rh(m+1)-aeroc%rh(m))
              EXIT
           ENDIF
        ENDDO
     ENDIF
 
     IF ( j2 == 0 ) THEN
        r_eff = aeroc%reff(j1,itype )
     ELSE
        r_eff = (1_kind_real-h1)*aeroc%reff(j1,itype ) + h1*aeroc%reff(j2,itype )
     ENDIF
 
   END FUNCTION gocart_aerosol_size
 
 
   FUNCTION upper2lower(str) RESULT(string)
 
     IMPLICIT NONE
 
     CHARACTER(*), INTENT(in) :: str
     CHARACTER(LEN(str))      :: string
 
     INTEGER :: ic, i
 
     CHARACTER(26), PARAMETER :: upper = 'ABCDEFGHIJKLMNOPQRSTUVWXYZ'
     CHARACTER(26), PARAMETER :: lower = 'abcdefghijklmnopqrstuvwxyz'
 
!   lowcase each letter if it is lowecase
     string = str
     DO i = 1, len_trim(str)
        ic = index(upper, str(i:i))
        IF (ic > 0) string(i:i) = lower(ic:ic)
     END DO
 
   END FUNCTION upper2lower
 
   INTEGER FUNCTION getindex(names,usrname)
     IMPLICIT NONE
     CHARACTER(len=*),INTENT(in) :: names(:)
     CHARACTER(len=*),INTENT(in) :: usrname
     INTEGER i
     getindex=-1
     DO i=1,SIZE(names)
        IF(cmp_strings(usrname, names(i))) THEN
           getindex=i
           EXIT
        ENDIF
     ENDDO
   END FUNCTION getindex
 
!from fv3
 
   SUBROUTINE qsmith_atm(atm,rh)
 
     TYPE(crtm_atmosphere_type), INTENT(in) :: atm(:)
     REAL(kind_real), INTENT(out),DIMENSION(:,:):: rh
 
     REAL, ALLOCATABLE :: table(:),des(:)
 
     REAL es, qs, q
     REAL ap1, eps10
     REAL Tmin
     INTEGER i, k, it, n_layers, n_profiles
 
     n_layers=SIZE(rh,1)
     n_profiles=SIZE(rh,2)
 
     tmin = tice-160.
     eps10  = 10.*esl
 
     IF( .NOT. ALLOCATED(table) ) CALL qsmith_init(table,des)
 
     DO k=1,n_layers
        DO i=1,n_profiles
           ap1 = 10.*dim(atm(i)%Temperature(k), tmin) + 1.
           ap1 = min(2621., ap1)
           it = ap1
           es = table(it) + (ap1-it)*des(it)
           q=atm(i)%Absorber(k,1)*1.e-3/(1.+atm(i)%Absorber(k,1)*1.e-3)
           qs = esl*es*(1.+zvir*q)/(atm(i)%Pressure(k)*100.)
           rh(k,i) = q/qs
        ENDDO
     ENDDO
 
   END SUBROUTINE qsmith_atm
 
   SUBROUTINE qsmith_profiles(t,sphum,p,rh)
 
     REAL(kind_real), DIMENSION(:,:), INTENT(in) :: t,sphum,p
     REAL(kind_real), DIMENSION(:,:), INTENT(out) :: rh
 
     REAL, ALLOCATABLE :: table(:),des(:)
 
     REAL es, qs, q
     REAL ap1, eps10
     REAL Tmin
     INTEGER i, k, it, n_layers, n_profiles
 
     n_layers=SIZE(t,1)
     n_profiles=SIZE(t,2)
 
     tmin = tice-160.
     eps10  = 10.*esl
 
     IF ( .NOT. ALLOCATED(table) ) CALL qsmith_init(table,des)
 
     DO k=1,n_layers
        DO i=1,n_profiles
           ap1 = 10.*dim(t(k,i), tmin) + 1.
           ap1 = min(2621., ap1)
           it = ap1
           es = table(it) + (ap1-it)*des(it)
           q=sphum(k,i)
           qs = esl*es*(1.+zvir*q)/p(k,i)
           rh(k,i) = q/qs
        ENDDO
     ENDDO
 
   END SUBROUTINE qsmith_profiles
 
   SUBROUTINE qsmith_init(table,des)
     
     REAL, ALLOCATABLE, INTENT(out) :: table(:),des(:)
     INTEGER, PARAMETER:: length=2621 
     INTEGER i
     
     IF( .NOT. ALLOCATED(table) ) THEN
!                            Generate es table (dT = 0.1 deg. C)
        
        ALLOCATE ( table(length) )
        ALLOCATE (  des(length) )
        
        CALL qs_table(length, table)
        
        DO i=1,length-1
           des(i) = table(i+1) - table(i)
        ENDDO
        des(length) = des(length-1)
     ENDIF
     
   END SUBROUTINE qsmith_init
   
   SUBROUTINE qs_table(n,table)
     INTEGER, INTENT(in):: n
     REAL table (n)
     REAL :: dt=0.1
     REAL esbasw, tbasw, tbasi, Tmin, tem, aa, b, c, d, e
     INTEGER i
! Constants
     esbasw = 1013246.0
     tbasw =   373.16
     tbasi =   273.16
     tmin = tbasi - 160.
!  Compute es over water
!  see smithsonian meteorological tables page 350.
     DO  i=1,n
        tem = tmin+dt*real(i-1)
        aa  = -7.90298*(tbasw/tem-1)
        b   =  5.02808*alog10(tbasw/tem)
        c   = -1.3816e-07*(10**((1-tem/tbasw)*11.344)-1)
        d   =  8.1328e-03*(10**((tbasw/tem-1)*(-3.49149))-1)
        e   =  alog10(esbasw)
        table(i)  = 0.1*10**(aa+b+c+d+e)
     ENDDO
     
   END SUBROUTINE qs_table
   
END MODULE ufo_crtm_utils_mod