doxygen-example/_o_d_p_s___coordinate_mapping_8f90_source.html

 !
 ! ODPS_CoordinateMapping
 !
 ! Module containing routines to perform data mapping from user pressure space to
 ! internal pressure space for the ODPS algorithm.
 !
 !
 ! CREATION HISTORY:
 !       Written by:     Yong Han & Yong Chen, JCSDA, NOAA/NESDIS 10-NOV-2009


 MODULE odps_coordinatemapping

   ! -----------------
   ! Environment setup
   ! -----------------
   ! Module use
   USE type_kinds,                ONLY: fp
   USE crtm_atmosphere_define,    ONLY: crtm_atmosphere_type, h2o_id
   USE crtm_geometryinfo_define,  ONLY: crtm_geometryinfo_type, &
                                        crtm_geometryinfo_getvalue
   USE odps_define,               ONLY: odps_type
   USE pafv_define,               ONLY: pafv_type, &
                                        pafv_associated
   USE profile_utility_parameters,ONLY: g0, eps, r_dryair
   USE crtm_parameters,           ONLY: zero, one, two, &
                                        earth_radius,   &
                                        degrees_to_radians

   ! Disable implicit typing
   IMPLICIT NONE

   ! ------------
   ! Visibilities
   ! ------------
   ! Everything private by default
   PRIVATE
   ! variables and routines from USE modules
   PUBLIC :: map_input
   PUBLIC :: map_input_tl
   PUBLIC :: map_input_ad
   PUBLIC :: interpolate_profile
   PUBLIC :: interpolate_profile_f1_tl
   PUBLIC :: interpolate_profile_f1_ad
   PUBLIC :: compute_interp_index

   ! Parameters used in the geopotential height calculation routines
   ! a factor used in the virtual temperature Tv calculation
   REAL(fp), PARAMETER :: c = (one/eps - one) / 1000.0_fp
   ! a factor used in the scale height calculation ( H = CC*Tv)
   REAL(fp), PARAMETER :: cc = 0.001_fp*r_dryair/g0
   ! for using in SUBROUTINE LayerAvg
   REAL(fp), PARAMETER :: smalldiff  = 1.0e-20_fp

   ! -----------------
   ! Module parameters
   ! -----------------
   ! RCS Id for the module
   CHARACTER(*), PRIVATE, PARAMETER :: module_rcs_id = &
   '$Id: ODPS_CoordinateMapping.f90 60152 2015-08-13 19:19:13Z paul.vandelst@noaa.gov $'

 CONTAINS


 !################################################################################
 !################################################################################
 !##                                                                            ##
 !##                         ## PUBLIC MODULE ROUTINES ##                       ##
 !##                                                                            ##
 !################################################################################
 !################################################################################

 !------------------------------------------------------------------------------
 !
 ! NAME:
 !       Map_Input
 !
 ! PURPOSE:
 !       Subroutine for data mapping from user pressure space to
 !       internal pressure space for the ODPS algorithm.
 !
 ! CALLING SEQUENCE:
 !       CALL Map_Input( Atm,                   & ! Input
 !                       TC,                    & ! Input
 !                       GeoInfo,               & ! Input
 !                       Temperature,           & ! Output
 !                       Absorber,              & ! Output
 !                       User_Level_LnPressure, & ! Output
 !                       Ref_Level_LnPressure,  & ! Output
 !                       Secant_Zenith,         & ! Output
 !                       PAFV )                   ! In/Output
 !
 ! INPUT ARGUMENTS:
 !
 !        Atm       :     CRTM Atmosphere structure containing the atmospheric
 !                        state data.
 !                        UNITS:      N/A
 !                        TYPE:       CRTM_Atmosphere_type
 !                        DIMENSION:  Scalar
 !                        ATTRIBUTES: INTENT(IN)
 !
 !          TC:           ODPS structure holding tau coefficients
 !                        UNITS:      N/A
 !                        TYPE:       ODPS_type
 !                        DIMENSION:  Scalar
 !                        ATTRIBUTES: INTENT(IN)
 !
 !       GeoInfo     :    CRTM_GeometryInfo structure containing the
 !                        view geometry information.
 !                        UNITS:      N/A
 !                        TYPE:       CRTM_GeometryInfo_type
 !                        DIMENSION:  Scalar
 !                        ATTRIBUTES: INTENT(IN)
 !
 ! OUTPUT ARGUMENTS:
 !
 !       Temperature :    Temperatures on internal pressure grids
 !                        UNITS:      K
 !                        TYPE:       REAL(fp)
 !                        DIMENSION:  Rank 1 array
 !                        ATTRIBUTES: INTENT(IN)
 !
 !         Absorber :     Layer absorber amount on internal pressure grids
 !                        UNITS:      depend on absorber types
 !                        TYPE:       REAL(fp)
 !                        DIMENSION:  Rank 2 array (n_layers x n_abosrbers)
 !                        ATTRIBUTES: INTENT(IN)
 !
 ! User_Level_LnPressure: User level pressure coordinate in Log scale
 !                        UNITS:      N/A  (Pressure in hPa)
 !                        TYPE:       REAL(fp)
 !                        DIMENSION:  Rank 1 array (0:n_userLayers)
 !                        ATTRIBUTES: INTENT(IN)
 !
 ! Ref_Level_LnPressure:  internal level pressure coordinate in Log scale
 !                        UNITS:      N/A  (Pressure in hPa)
 !                        TYPE:       REAL(fp)
 !                        DIMENSION:  Rank 1 array (0:n_refLayers)
 !                        ATTRIBUTES: INTENT(IN)
 !       Secant_Zenith :  Secant zenith angle array on internal pressure grids
 !                        UNITS:      N/A
 !                        TYPE:       REAL(fp)
 !                        DIMENSION:  Rank 1 array (n_refLayers)
 !                        ATTRIBUTES: INTENT(IN)
 !
 ! IN/OUTPUT ARGUMENTS:
 !                 PAFV:  Structure containing FW variables for TL and AD uses
 !                        UNITS:      N/A
 !                        TYPE:       PAFV_type
 !                        DIMENSION:  Scalar
 !                        ATTRIBUTES: INTENT(IN OUT)
 !------------------------------------------------------------------------------
   SUBROUTINE map_input(Atm,                   & ! Input
                        TC,                    & ! Input
                        GeoInfo,               & ! Input
                        Temperature,           & ! Output
                        Absorber,              & ! Output
                        User_Level_LnPressure, & ! Output
                        Ref_Level_LnPressure,  & ! Output
                        Secant_Zenith,         & ! Output
                        H2O_idx,               & ! Output
                        PAFV)                    ! In/Output
     TYPE(crtm_atmosphere_type)   , INTENT(IN)     :: atm
     TYPE(odps_type)              , INTENT(IN)     :: tc
     TYPE(crtm_geometryinfo_type) , INTENT(IN)     :: geoinfo
     REAL(fp)                     , INTENT(OUT)    :: temperature(:)
     REAL(fp)                     , INTENT(OUT)    :: absorber(:,:)
     REAL(fp)                     , INTENT(OUT)    :: user_level_lnpressure(0:)
     REAL(fp)                     , INTENT(OUT)    :: ref_level_lnpressure(0:)
     REAL(fp)                     , INTENT(OUT)    :: secant_zenith(:)
     INTEGER                      , INTENT(OUT)    :: h2o_idx
     TYPE(pafv_type)              , INTENT(INOUT)  :: pafv

     ! Local variables
     REAL(fp) :: tolerance
     REAL(fp) :: sineang
     REAL(fp) :: odps_sfc_fraction, z_offset, s
     REAL(fp) :: z(0:tc%n_layers)  ! Heights of pressure levels
     REAL(fp) :: acc_weighting(atm%n_layers, tc%n_layers)
     INTEGER  :: interp_index(2, tc%n_layers)
     REAL(fp) :: ref_lnpressure(tc%n_layers)
     REAL(fp) :: user_lnpressure(atm%n_layers)
     REAL(fp) :: surface_altitude, sensor_zenith_radian
     ! absorber index mapping from ODPS to user
     INTEGER  :: idx_map(tc%n_absorbers)
     INTEGER  :: j, jj, k, n_odps_layers, n_user_layers, odps_sfc_idx

     ! Define numerical precision tolerance.
     tolerance = epsilon(one)

     n_odps_layers = tc%n_Layers
     n_user_layers = atm%n_layers
     ! Set pressure profiles for interpolations
     ref_lnpressure  = log(tc%Ref_Pressure)
     user_lnpressure = log(atm%Pressure(1:n_user_layers))
     ref_level_lnpressure = log(tc%Ref_Level_Pressure)
     IF(atm%Level_Pressure(0) <= zero)THEN
       ! In this bad case, the top pressure level is set to the half of the next-to-top pressure level
       user_level_lnpressure(0) = log(atm%Level_Pressure(1)/two)
     ELSE
       user_level_lnpressure(0) = log(atm%Level_Pressure(0))
     END IF
     user_level_lnpressure(1:n_user_layers) = log(atm%Level_Pressure(1:n_user_layers))

     ! Find the index at which the ODPS layer contains the user surface pressure level. Then
     ! compute the fraction of the portion between the lower boundary of the ODPS layer and
     ! the user surface pressure level
     odps_sfc_idx = n_odps_layers
     odps_sfc_fraction = zero
     IF(tc%Ref_Level_Pressure(n_odps_layers) > atm%Level_Pressure(n_user_layers))THEN
       DO k = n_odps_layers, 0, -1
         IF(tc%Ref_Level_Pressure(k) < atm%Level_Pressure(n_user_layers))THEN
           odps_sfc_idx = k+1
           odps_sfc_fraction = (ref_level_lnpressure(odps_sfc_idx) - &
                                user_level_lnpressure(n_user_layers)) /&
                               (ref_level_lnpressure(odps_sfc_idx) - &
                                ref_level_lnpressure(k))
           EXIT
         END IF
       END DO
     END IF

     ! Extract the needed GeometryInfo information
     CALL crtm_geometryinfo_getvalue( geoinfo, &
                                      surface_altitude = surface_altitude, &
                                      trans_zenith_radian = sensor_zenith_radian )


     !-----------------------------------------------------------
     ! Interpolate temperautre profile on internal pressure grids
     !-----------------------------------------------------------
     CALL layeravg( ref_lnpressure   , &
                    user_lnpressure  , &
                    acc_weighting    , &
                    interp_index)

     DO k = 1, n_odps_layers
        temperature(k) = sum(acc_weighting(interp_index(1,k):interp_index(2,k), k)  &
                            * atm%Temperature(interp_index(1,k):interp_index(2,k)) )
     END DO

     !-----------------------------------------------------------
     ! Interpolate absorber profiles on internal pressure grids
     !-----------------------------------------------------------
     DO j = 1,tc%n_Absorbers
       idx_map(j) = -1
       DO jj=1, atm%n_Absorbers
        IF( atm%Absorber_ID(jj) == tc%Absorber_ID(j) ) THEN
         idx_map(j) = jj
         EXIT
        END IF
       END DO

       ! save index for water vapor absorption
       h2o_idx = 1  ! default
       IF(tc%Absorber_ID(j) == h2o_id)THEN
         h2o_idx = j
       END IF
       IF(idx_map(j) > 0)THEN

         DO k = 1, n_odps_layers
             absorber(k,j) = sum(acc_weighting(interp_index(1,k):interp_index(2,k), k)  &
                                * atm%Absorber(interp_index(1,k):interp_index(2,k), idx_map(j)) )
         END DO

         DO k=1, n_odps_layers

           IF (absorber(k,j) <= tolerance ) absorber(k,j) = tolerance
           IF (absorber(k,j) < tc%Min_Absorber(k,j) ) absorber(k,j) = tc%Min_Absorber(k,j)
           IF (absorber(k,j) > tc%Max_Absorber(k,j) ) absorber(k,j) = tc%Max_Absorber(k,j)

         END DO
       ELSE ! when the profile is missing, use the referece profile
         absorber(:, j) = tc%Ref_Absorber(:, j)
       END IF

     END DO

     !-----------------------------------------------
     ! Compute height dependent secant zenith angles
     !-----------------------------------------------
     ! Compute geopotential height, which starts from ODPS surface pressure level
     CALL geopotential_height( tc%Ref_Level_Pressure,       &
                               tc%Ref_Temperature,          &
                               tc%Ref_Absorber(:, h2o_idx), &
                               zero,                        &
                               z)
     ! Adjust ODPS surface height for the user surface height. The adjustment includes two parts:
     ! (1) the delta Z from the ODPS surface pressure to the user surface pressure
     ! (2) the user-given surface height
     IF(tc%Ref_Level_Pressure(n_odps_layers) >= atm%Level_Pressure(n_user_layers))THEN
       z_offset = -(z(odps_sfc_idx) + odps_sfc_fraction*(z(odps_sfc_idx-1)-z(odps_sfc_idx))) &
                  + surface_altitude
     ELSE
       ! For the case in which the user surface pressure is larger than the ODPS surface pressure,
       ! the ODPS surface is adjusted for a surface pressure 1013 mb, regardless the user supplied
       ! surface height.
       z_offset = cc*tc%Ref_Temperature(n_odps_layers) &  ! scale height
                  *log(1013.0_fp / tc%Ref_Level_Pressure(n_odps_layers))
     END IF
     z = z + z_offset

     s = (earth_radius + surface_altitude)*sin(sensor_zenith_radian)

     DO k = 1, n_odps_layers
       sineang = s /(earth_radius + z(k))
       secant_zenith(k) = one / sqrt(one - sineang*sineang)
     END DO

     IF ( pafv_associated(pafv) ) THEN
       pafv%Temperature  = temperature
       pafv%Absorber     = absorber
       pafv%interp_index = interp_index
       pafv%Acc_Weighting= acc_weighting
       pafv%idx_map      = idx_map
       pafv%H2O_idx      = h2o_idx
       pafv%Ref_LnPressure = ref_lnpressure
       pafv%User_LnPressure = user_lnpressure
     END IF

   END SUBROUTINE map_input

 !------------------------------------------------------------------------------
 !
 ! NAME:
 !       Map_Input_TL
 !
 ! PURPOSE:
 !       TL Subroutine for data mapping from user pressure space to
 !       internal pressure space for the ODPS algorithm.
 !
 ! CALLING SEQUENCE:
 !       CALL Map_Input_TL( TC,                    & ! Input
 !                          Atm_TL,                & ! Input
 !                          Temperature_TL,        & ! Output
 !                          Absorber_TL,           & ! Output
 !                          PAFV )                   ! In/Output
 !
 ! INPUT ARGUMENTS:
 !          TC:           ODPS structure holding tau coefficients
 !                        UNITS:      N/A
 !                        TYPE:       ODPS_type
 !                        DIMENSION:  Scalar
 !                        ATTRIBUTES: INTENT(IN)
 !
 !        Atm_TL    :     TL CRTM Atmosphere structure containing the atmospheric
 !                        state data.
 !                        UNITS:      N/A
 !                        TYPE:       CRTM_Atmosphere_type
 !                        DIMENSION:  Scalar
 !                        ATTRIBUTES: INTENT(IN)
 !
 !                 PAFV:  Structure containing FW variables for TL and AD uses
 !                        UNITS:      N/A
 !                        TYPE:       PAFV_type
 !                        DIMENSION:  Scalar
 !                        ATTRIBUTES: INTENT(IN)
 ! OUTPUT ARGUMENTS:
 !
 !       Temperature_TL : TL Temperatures on internal pressure grids
 !                        UNITS:      K
 !                        TYPE:       REAL(fp)
 !                        DIMENSION:  Rank 1 array
 !                        ATTRIBUTES: INTENT(OUT)
 !
 !         Absorber_TL :  TL Layer absorber amount on internal pressure grids
 !                        UNITS:      depend on absorber types
 !                        TYPE:       REAL(fp)
 !                        DIMENSION:  Rank 2 array (n_layers x n_abosrbers)
 !                        ATTRIBUTES: INTENT(OUT)
 !
 !
 !------------------------------------------------------------------------------
   SUBROUTINE map_input_tl(TC,              & ! Input
                           Atm_TL,          & ! Input
                           Temperature_TL,  & ! Output
                           Absorber_TL,     & ! Output
                           PAFV)              ! Input
     TYPE(odps_type)              , INTENT(IN)     :: tc
     TYPE(crtm_atmosphere_type)   , INTENT(IN)     :: atm_tl
     REAL(fp)                     , INTENT(OUT)    :: temperature_tl(:)
     REAL(fp)                     , INTENT(OUT)    :: absorber_tl(:,:)
     TYPE(pafv_type)              , INTENT(IN)     :: pafv

     ! Local variables
     REAL(fp) :: tolerance
     ! absorber index mapping from ODPS to user
     INTEGER  :: idx_map(tc%n_absorbers), h2o_idx
     INTEGER  :: j, k, n_odps_layers

     ! Define numerical precision tolerance.
     tolerance = epsilon(one)

     n_odps_layers = tc%n_Layers

     !-----------------------------------------------------------
     ! Interpolate temperautre profile on internal pressure grids
     !-----------------------------------------------------------
     DO k = 1, n_odps_layers
        temperature_tl(k) = &
          sum(pafv%Acc_Weighting(pafv%interp_index(1,k):pafv%interp_index(2,k), k)  &
              * atm_tl%Temperature(pafv%interp_index(1,k):pafv%interp_index(2,k)) )
     END DO

     !-----------------------------------------------------------
     ! Interpolate absorber profiles on internal pressure grids
     !-----------------------------------------------------------
     h2o_idx = pafv%H2O_idx
     idx_map = pafv%Idx_map

     DO j = 1,tc%n_Absorbers
       IF(idx_map(j) > 0)THEN

         DO k = 1, n_odps_layers
            absorber_tl(k,j) = &
             sum(pafv%Acc_Weighting(pafv%interp_index(1,k):pafv%interp_index(2,k), k)  &
                 * atm_tl%Absorber(pafv%interp_index(1,k):pafv%interp_index(2,k), idx_map(j)) )
         END DO

         DO k=1, n_odps_layers
           IF(pafv%Absorber(k,j) <= tolerance ) absorber_tl(k,j) = zero

         END DO
       ELSE ! when the profile is missing, use the referece profile
         absorber_tl(:, j) = zero
       END IF

     END DO

   END SUBROUTINE map_input_tl

 !------------------------------------------------------------------------------
 !
 ! NAME:
 !       Map_Input_AD
 !
 ! PURPOSE:
 !       AD Subroutine for data mapping from user pressure space to
 !       internal pressure space for the ODPS algorithm.
 !
 ! CALLING SEQUENCE:
 !       CALL Map_Input_AD( TC,                    & ! Input
 !                          Temperature_AD,        & ! Input
 !                          Absorber_AD,           & ! Input
 !                          Atm_AD,                & ! In/Output
 !                          PAFV )                   ! Input
 !
 ! INPUT ARGUMENTS:
 !          TC:           ODPS structure holding tau coefficients
 !                        UNITS:      N/A
 !                        TYPE:       ODPS_type
 !                        DIMENSION:  Scalar
 !                        ATTRIBUTES: INTENT(IN)
 !
 !       Temperature_AD : AD Temperatures on internal pressure grids
 !                        UNITS:      K
 !                        TYPE:       REAL(fp)
 !                        DIMENSION:  Rank 1 array
 !                        ATTRIBUTES: INTENT(IN)
 !
 !         Absorber_AD :  AD Layer absorber amount on internal pressure grids
 !                        UNITS:      depend on absorber types
 !                        TYPE:       REAL(fp)
 !                        DIMENSION:  Rank 2 array (n_layers x n_abosrbers)
 !                        ATTRIBUTES: INTENT(IN)
 !
 !                 PAFV:  Structure containing FW variables for TL and AD uses
 !                        UNITS:      N/A
 !                        TYPE:       PAFV_type
 !                        DIMENSION:  Scalar
 !                        ATTRIBUTES: INTENT(IN)
 ! OUTPUT ARGUMENTS:
 !
 !        Atm_AD    :     AD CRTM Atmosphere structure containing the atmospheric
 !                        state data.
 !                        UNITS:      N/A
 !                        TYPE:       CRTM_Atmosphere_type
 !                        DIMENSION:  Scalar
 !                        ATTRIBUTES: INTENT(INOUT)
 !
 !------------------------------------------------------------------------------
   SUBROUTINE map_input_ad(TC,              & ! Input
                           Temperature_AD,  & ! Input
                           Absorber_AD,     & !  Input
                           Atm_AD,          & ! Output
                           PAFV)              ! Input

     TYPE(odps_type)              , INTENT(IN)     :: tc
     REAL(fp)                     , INTENT(INOUT)  :: temperature_ad(:)
     REAL(fp)                     , INTENT(INOUT)  :: absorber_ad(:,:)
     TYPE(crtm_atmosphere_type)   , INTENT(INOUT)  :: atm_ad
     TYPE(pafv_type)              , INTENT(IN)     :: pafv

     ! Local variables
     REAL(fp) :: tolerance
     ! absorber index mapping from ODPS to user
     INTEGER  :: idx_map(tc%n_absorbers), h2o_idx
     INTEGER  :: j, k, n_odps_layers

     ! Define numerical precision tolerance.
     tolerance = epsilon(one)

     !-----------------------------------------------------------
     ! Initialization of forward model part
     !-----------------------------------------------------------

     n_odps_layers = tc%n_Layers

     h2o_idx = pafv%H2O_idx
     idx_map = pafv%Idx_map

     !-----------------------------------------------------------
     ! Interpolate absorber profiles on internal pressure grids
     !-----------------------------------------------------------
     DO j = tc%n_Absorbers, 1, -1
       IF(idx_map(j) > 0)THEN

         DO k=1, n_odps_layers
           IF(pafv%Absorber(k,j) <= tolerance ) absorber_ad(k,j) = zero
         END DO

         DO k = n_odps_layers, 1, -1
            atm_ad%Absorber(pafv%interp_index(1,k):pafv%interp_index(2,k), idx_map(j)) = &
               atm_ad%Absorber(pafv%interp_index(1,k):pafv%interp_index(2,k), idx_map(j)) &
             + pafv%Acc_Weighting(pafv%interp_index(1,k):pafv%interp_index(2,k), k) &
             * absorber_ad(k,j)
         END DO

       ELSE ! when the profile is missing, use the referece profile
         absorber_ad(:, j) = zero
       END IF
     END DO

     !-----------------------------------------------------------
     ! Interpolate temperautre profile on internal pressure grids
     !-----------------------------------------------------------

     DO k = n_odps_layers, 1, -1
       atm_ad%Temperature(pafv%interp_index(1,k):pafv%interp_index(2,k)) = &
          atm_ad%Temperature(pafv%interp_index(1,k):pafv%interp_index(2,k))  &
        + pafv%Acc_Weighting(pafv%interp_index(1,k):pafv%interp_index(2,k), k) &
        * temperature_ad(k)

     END DO

     temperature_ad = zero
     absorber_ad    = zero

   END SUBROUTINE map_input_ad

 !--------------------------------------------------------------------------------
 !S+
 ! NAME:
 !       Geopotential_Height
 !
 ! PURPOSE:
 !       Routine to calculate geopotential height using the hypsometric
 !       equation.
 !
 ! CALLING SEQUENCE:
 !       CALL Geopotential_Height( Level_Pressure,                                    &  ! Input
 !                                 Temperature,                                       &  ! Input
 !                                 Water_Vapor,                                       &  ! Input
 !                                 Surface_Height,                                    &  ! input
 !                                 Level_Height)
 !
 ! INPUT ARGUMENTS:
 !       Level_Pressure:            Level Pressures, which are the boundaries of of the
 !                                  atmospheric layers for the Temperature and Water_vapor arrays
 !                                  UNITS:      hectoPascals, hPa
 !                                  TYPE:       REAL(fp)
 !                                  DIMENSION:  Rank-1 (0:n_layers)
 !                                  ATTRIBUTES: INTENT(IN)
 !
 !       Temperature:               Layer Temperature.
 !                                  UNITS:      Kelvin, K
 !                                  TYPE:       REAL(fp)
 !                                  DIMENSION:  Rank-1 (n_layers)
 !                                  ATTRIBUTES: INTENT(IN)
 !
 !       Surface_Height:            Height of Level_Pressure(n_layer)
 !                                  input arrays.
 !                                  UNITS:      km.
 !                                  TYPE:       REAL(fp)
 !                                  DIMENSION:  Scalar
 !                                  ATTRIBUTES: INTENT(IN)
 !
 !
 !       Water_Vapor:               Layer water vapor mixing radio
 !                                  UNITS:      g/kg
 !                                  TYPE:       REAL(fp)
 !                                  DIMENSION:  Same as Temperature
 !                                  ATTRIBUTES: INTENT(IN)
 !
 ! OUTPUT ARGUMENTS:
 !       Level_Height:              Geopotential Heights of the input Pressure levels.
 !                                  UNITS:      metres, m
 !                                  TYPE:       REAL(fp)
 !                                  DIMENSION:  Same as input Pressure
 !                                  ATTRIBUTES: INTENT(OUT)
 !
 !
 !
 ! CREATION HISTORY:
 !       Written by:     Yong Han, NOAA/NESDIS 5-Feb-2008
 !S-
 !--------------------------------------------------------------------------------

   SUBROUTINE geopotential_height( Level_Pressure,           &  ! Input
                                   Temperature,              &  ! Input
                                   Water_Vapor,              &  ! Input
                                   Surface_Height,           &  ! input
                                   Level_Height)                ! Output

     ! Arguments
     REAL(fp),               INTENT(IN)  :: Level_Pressure(0:)
     REAL(fp),               INTENT(IN)  :: Temperature(:)
     REAL(fp),               INTENT(IN)  :: Water_Vapor(:)
     REAL(fp),               INTENT(IN)  :: Surface_Height
     REAL(fp),               INTENT(OUT) :: Level_Height(0:)

     ! Function result
     REAL(fp) :: Tv, H
     INTEGER  :: k, n_Layers

     n_layers =  SIZE(temperature)
     level_height(n_layers) = surface_height
     DO k = n_layers, 1, -1
       ! virtual temperature computed using an approximation of the exact formula:
       !  Tv = T*(1+w/epsilon)/(1+w), where w is the water vapor mixing ratio
       tv = temperature(k)*(one + c*water_vapor(k))
       h  = cc*tv
       level_height(k-1) = level_height(k) + h*log(level_pressure(k)/level_pressure(k-1))
     END DO

   END SUBROUTINE geopotential_height


 !---------------------------------------------------------------------------------------------
 ! NAME: Interpolate_Profile
 !
 ! PURPOSE:
 !    Given x and u that are ascending arrays, it interpolates y with the abscissa x
 !    on the abscissa u using the following algorithm:
 !       y_int(i) = y(1)  if u(i) < x(1)
 !       y_int(i) = y(nx) if u(i) > x(nx)
 !       y_int(i) = y(ix1) + (y(ix2)-y(ix1))*(u(i) - x(ix1))/(x(ix2)-x(ix1))
 !                        if x(ix1) <= u(i) <= x(ix2)
 !
 !    IThe index array interp_index contains the following content
 !
 !      interp_index(1, i) = 1 and interp_index(2, i) = 1, if u(i) < x(1)
 !      interp_index(1, i) = nx and interp_index(2, i) = nx, if u(i) > x(nx),
 !                                                          where nx = SIZE(x)
 !      x(interp_index(1, i)) <= u(i) <= x(interp_index(2, i)) if x(1) <= u(i) <= x(nx)
 !
 ! CALLING SEQUENCE:
 !            CALL Interpolate_Profile(interp_index, y, x, u, y_int)
 !
 ! INPUT ARGUMENTS:
 !       y:            The data array to be interpolated.
 !                     UNITS:      N/A
 !                     TYPE:       fp
 !                     DIMENSION:  rank-1
 !                     ATTRIBUTES: INTENT(IN)
 !
 !       x:            The abscissa values for y and they must be monotonically
 !                     ascending.
 !                     UNITS:      N/A
 !                     TYPE:       fp
 !                     DIMENSION:  rank-1
 !                     ATTRIBUTES: INTENT(IN)
 !
 !       u:            The abscissa values for the results
 !                     and they must be monotonically ascending
 !                     UNITS:      N/A
 !                     TYPE:       fp
 !                     DIMENSION:  rank-1
 !                     ATTRIBUTES: INTENT(IN)
 !
 !    interp_index:    The index array of dimension (nu x 2), where nu = SIZE(u)
 !                     UNITS:      N/A
 !                     TYPE:       Integer
 !                     DIMENSION:  rank-2
 !                     ATTRIBUTES: INTENT(IN)
 !
 ! OUTPUT ARGUMENTS:
 !       y_int:        The array contains the results
 !                     UNITS:      N/A
 !                     TYPE:       Integer
 !                     DIMENSION:  rank-1
 !                     ATTRIBUTES: INTENT(OUT)
 !
 ! RESTRICTIONS:
 !     To be efficient, this routine does not check that x and u are both
 !     monotonically ascending and the index bounds.
 !
 ! CREATION HISTORY:
 !       Written by:     Yong Han, 10-Dec-2008
 !-----------------------------------------------------------------------------------

   !----------------------------------------------------------------------------
   !  Interpolation routine with interperlation index array already calculated.
   !----------------------------------------------------------------------------
   SUBROUTINE interpolate_profile(interp_index, y, x, u, y_int)
     ! Arguments
     INTEGER,  DIMENSION(:,:), INTENT(IN)  :: interp_index
     REAL(fp), DIMENSION(:),   INTENT(IN)  :: y
     REAL(fp), DIMENSION(:),   INTENT(IN)  :: x
     REAL(fp), DIMENSION(:),   INTENT(IN)  :: u
     REAL(fp), DIMENSION(:),   INTENT(OUT) :: y_int
     ! Local variables
     INTEGER :: i, n, k1, k2

     n = SIZE(interp_index, dim=2)
     DO i = 1, n
       k1 = interp_index(1, i)
       k2 = interp_index(2, i)
       IF( k1 == k2)THEN
         y_int(i) = y(k1)
       ELSE
         CALL interp_linear(y(k1), x(k1), y(k2), x(k2), u(i), y_int(i))
       END IF
     END DO

   END SUBROUTINE interpolate_profile


 !---------------------------------------------------------------------------------------------
 ! NAME: Interpolate_Profile_TL
 !
 ! PURPOSE:
 !     The Tangent_Linear routine of Interpolate_Profile
 ! CALLING SEQUENCE:
 !           CALL Interpolate_Profile_F1_TL(interp_index, y, x, u, y_TL, y_int_TL)
 !      OR   CALL Interpolate_Profile_F2_TL(interp_index, y, x, u, y_TL, u_TL, y_int_TL)
 !      OR   CALL Interpolate_Profile_F3_TL(interp_index, y, x, u, y_TL, x_TL, y_int_TL)
 !
 ! INPUT ARGUMENTS:
 !       y:            The data array to be interpolated.
 !                     UNITS:      N/A
 !                     TYPE:       fp
 !                     DIMENSION:  rank-1
 !                     ATTRIBUTES: INTENT(IN)
 !
 !       x:            The abscissa values for y and they must be monotonically
 !                     ascending.
 !                     UNITS:      N/A
 !                     TYPE:       fp
 !                     DIMENSION:  rank-1
 !                     ATTRIBUTES: INTENT(IN)
 !
 !       u:            The abscissa values for the results
 !                     and they must be monotonically ascending
 !                     UNITS:      N/A
 !                     TYPE:       fp
 !                     DIMENSION:  rank-1
 !                     ATTRIBUTES: INTENT(IN)
 !
 !       y_TL:         The Tangent-linear data array of y
 !                     UNITS:      N/A
 !                     TYPE:       fp
 !                     DIMENSION:  rank-1
 !                     ATTRIBUTES: INTENT(IN)
 !
 !       u_TL:         The Tangent-linear data array of u
 !                     UNITS:      N/A
 !                     TYPE:       fp
 !                     DIMENSION:  rank-1
 !                     ATTRIBUTES: INTENT(IN)
 !
 !       x_TL:         The Tangent-linear data array of x
 !                     UNITS:      N/A
 !                     TYPE:       fp
 !                     DIMENSION:  rank-1
 !                     ATTRIBUTES: INTENT(IN)
 !
 !    interp_index:    The index array of dimension (nu x 2), where nu = SIZE(u)
 !                     UNITS:      N/A
 !                     TYPE:       Integer
 !                     DIMENSION:  rank-2
 !                     ATTRIBUTES: INTENT(IN)
 !
 ! OUTPUT ARGUMENTS:
 !       y_int_TL:     The Tangent-linear array of y_int
 !                     UNITS:      N/A
 !                     TYPE:       Integer
 !                     DIMENSION:  rank-1
 !                     ATTRIBUTES: INTENT(OUT)
 !
 ! RESTRICTIONS:
 !     To be efficient, this routine does not check that x and u are both
 !     monotonically ascending and the index bounds.
 !
 ! CREATION HISTORY:
 !       Written by:     Yong Han, 10-Dec-2008
 !-----------------------------------------------------------------------------------

   SUBROUTINE interpolate_profile_f1_tl(interp_index, x, u, y_TL, y_int_TL)
     ! Arguments
     INTEGER,  DIMENSION(:,:), INTENT(IN)  :: interp_index
     REAL(fp), DIMENSION(:),   INTENT(IN)  :: x
     REAL(fp), DIMENSION(:),   INTENT(IN)  :: u
     REAL(fp), DIMENSION(:),   INTENT(IN)  :: y_tl
     REAL(fp), DIMENSION(:),   INTENT(OUT) :: y_int_tl
     ! Local variables
     INTEGER :: i, n, k1, k2

     n = SIZE(interp_index, dim=2)
     DO i = 1, n
       k1 = interp_index(1, i)
       k2 = interp_index(2, i)
       IF( k1 == k2)THEN
         y_int_tl(i) = y_tl(k1)
       ELSE
         CALL interp_linear_f1_tl(x(k1), x(k2), u(i), y_tl(k1), y_tl(k2), y_int_tl(i))
       END IF
     END DO

   END SUBROUTINE interpolate_profile_f1_tl


 !---------------------------------------------------------------------------------------------
 ! NAME: Interpolate_Profile_AD
 !
 ! PURPOSE:
 !     The Adjoint routine of Interpolate_Profile
 !
 ! CALLING SEQUENCE:
 !            CALL Interpolate_Profile_F1_AD(interp_index, y, x, u, y_int_AD, y_AD)
 !    OR      CALL Interpolate_Profile_F2_AD(interp_index, y, x, u, y_int_AD, y_AD, u_AD)
 !    OR      CALL Interpolate_Profile_F3_AD(interp_index, y, x, u, y_int_AD, y_AD, x_AD)
 !
 ! INPUT ARGUMENTS:
 !       y:            The data array to be interpolated.
 !                     UNITS:      N/A
 !                     TYPE:       fp
 !                     DIMENSION:  rank-1
 !                     ATTRIBUTES: INTENT(IN)
 !
 !       x:            The abscissa values for y and they must be monotonically
 !                     ascending.
 !                     UNITS:      N/A
 !                     TYPE:       fp
 !                     DIMENSION:  rank-1
 !                     ATTRIBUTES: INTENT(IN)
 !
 !       u:            The abscissa values for the results
 !                     and they must be monotonically ascending
 !                     UNITS:      N/A
 !                     TYPE:       fp
 !                     DIMENSION:  rank-1
 !                     ATTRIBUTES: INTENT(IN)
 !
 !       y_int_AD:     The Adjoint array of y_int
 !                     UNITS:      N/A
 !                     TYPE:       Integer
 !                     DIMENSION:  rank-1
 !                     ATTRIBUTES: INTENT(IN)
 !
 !    interp_index:    The index array of dimension (nu x 2), where nu = SIZE(u)
 !                     UNITS:      N/A
 !                     TYPE:       Integer
 !                     DIMENSION:  rank-2
 !                     ATTRIBUTES: INTENT(IN)
 !
 ! OUTPUT ARGUMENTS:
 !       y_AD:         The Adjoint data array of y
 !                     UNITS:      N/A
 !                     TYPE:       fp
 !                     DIMENSION:  rank-1
 !                     ATTRIBUTES: INTENT(IN)
 !
 !       u_AD:         The Adjoint data array of u
 !                     UNITS:      N/A
 !                     TYPE:       fp
 !                     DIMENSION:  rank-1
 !                     ATTRIBUTES: INTENT(IN)
 !
 !       x_AD:         The Adjoint data array of x
 !                     UNITS:      N/A
 !                     TYPE:       fp
 !                     DIMENSION:  rank-1
 !                     ATTRIBUTES: INTENT(IN)
 !
 ! RESTRICTIONS:
 !     To be efficient, this routine does not check that x and u are both
 !     monotonically ascending and the index bounds.
 !
 ! CREATION HISTORY:
 !       Written by:     Yong Han, 10-Dec-2008
 !-----------------------------------------------------------------------------------

   SUBROUTINE interpolate_profile_f1_ad(interp_index, x, u, y_int_AD, y_AD)
     ! Arguments
     INTEGER,  DIMENSION(:,:), INTENT(IN)      :: interp_index
     REAL(fp), DIMENSION(:),   INTENT(IN)      :: x
     REAL(fp), DIMENSION(:),   INTENT(IN)      :: u
     REAL(fp), DIMENSION(:),   INTENT(IN OUT)  :: y_int_ad
     REAL(fp), DIMENSION(:),   INTENT(IN OUT)  :: y_ad
     ! Local variables
     INTEGER :: i, n, k1, k2

     n = SIZE(interp_index, dim=2)
     DO i = n, 1, -1
       k1 = interp_index(1, i)
       k2 = interp_index(2, i)
       IF( k1 == k2)THEN
         y_ad(k1) = y_ad(k1) + y_int_ad(i)
         y_int_ad(i) = zero
       ELSE
         CALL interp_linear_f1_ad(x(k1), x(k2), u(i), y_int_ad(i), y_ad(k1), y_ad(k2))
       END IF
     END DO

   END SUBROUTINE interpolate_profile_f1_ad


 !---------------------------------------------------------------------------------------------
 ! NAME: Compute_Interp_Index
 !
 ! PURPOSE:
 !    Given x and u that are ascending arrays, it computes an index array, interp_index,
 !    such that
 !
 !      interp_index(i, 1) = 1 and interp_index(i, 2) = 1, if u(i) < x(1)
 !      interp_index(i, 1) = nx and interp_index(i, 2) = nx, if u(i) > x(nx),
 !                                                          where nx = SIZE(x)
 !      x(interp_index(i, 1)) <= u(i) <= x(interp_index(i, 2)) if x(1) <= u(i) <= x(nx)
 !
 ! CALLING SEQUENCE:
 !            CALL Compute_Interp_Index(x, u, interp_index)
 !
 ! INPUT ARGUMENTS:
 !       x:            The abscissa values for the data to be interpolated and
 !                     they must be monotonically ascending.
 !                     UNITS:      N/A
 !                     TYPE:       fp
 !                     DIMENSION:  rank-1
 !                     ATTRIBUTES: INTENT(IN)
 !
 !       u:            The abscissa values on which the data are interpolated
 !                     they must be monotonically ascending
 !                     the elements of array x.
 !                     UNITS:      N/A
 !                     TYPE:       fp
 !                     DIMENSION:  rank-1
 !                     ATTRIBUTES: INTENT(IN)
 !
 ! OUTPUT ARGUMENTS:
 !       interp_index: The index array of dimension (nu x 2), where nu = SIZE(u)
 !                     UNITS:      N/A
 !                     TYPE:       Integer
 !                     DIMENSION:  rank-2
 !                     ATTRIBUTES: INTENT(OUT)
 !
 ! RESTRICTIONS:
 !     To be efficient, this routine does not check that x and u are both
 !     monotonically ascending and the index bounds.
 !
 ! CREATION HISTORY:
 !       Written by:     Yong Han, 07-May-2004
 !-----------------------------------------------------------------------------------

   SUBROUTINE compute_interp_index(x, u, interp_index)
     ! Arguments
     REAL(fp), DIMENSION(:),   INTENT(IN)  :: x
     REAL(fp), DIMENSION(:),   INTENT(IN)  :: u
     INTEGER,  DIMENSION(:,:), INTENT(OUT) :: interp_index
     ! Local variables
     INTEGER :: nx, nu, ix, iu, j, k1, k2

     nx = SIZE(x)
     nu = SIZE(u)

     ! Set the indexes to 1 for the elements in u that are smaller than x(1)
     k1 = nu + 1
     lessthan_loop: DO iu = 1, nu
       IF(u(iu) < x(1))THEN
         interp_index(1, iu) = 1
         interp_index(2, iu) = 1
       ELSE
         k1 = iu
         EXIT lessthan_loop
       END IF
     END DO lessthan_loop

     ! Set the indexes to nx for the elements in u that are larger than x(nx)
     k2 = 0
     greaterthan_loop: DO iu = nu, k1, -1
       IF(u(iu) > x(nx))THEN
         interp_index(1, iu) = nx
         interp_index(2, iu) = nx
       ELSE
         k2 = iu
         EXIT greaterthan_loop
       END IF
     END DO greaterthan_loop

     ! Set the indexes for the elements in u that are in the range
     ! between x1(1) and x(nx)
     j = 1
     outer_loop: DO iu = k1, k2
       inner_loop: DO ix = j, nx-1
         IF(u(iu) >= x(ix) .AND. u(iu) <= x(ix+1))THEN
           interp_index(1, iu) = ix
           interp_index(2, iu) = ix+1
           j = ix
           EXIT inner_loop
         ENDIF
       END DO inner_loop
     END DO outer_loop

   END SUBROUTINE compute_interp_index


   !---------------------------------------------
   ! Function for two points linear interpolation
   !---------------------------------------------
   SUBROUTINE interp_linear(y1, x1, y2, x2, x, y)
     REAL(fp), INTENT(IN)  :: y1, x1, y2, x2, x
     REAL(fp), INTENT(OUT) :: y
     y = y1 + (y2-y1)*(x - x1)/(x2 - x1)
   END SUBROUTINE interp_linear


   SUBROUTINE interp_linear_f1_tl(x1, x2, x, y1_TL, y2_TL, y_TL)
     REAL(fp), INTENT(IN)  :: x1, x2, x, y1_TL, y2_TL
     REAL(fp), INTENT(OUT) :: y_TL
     y_tl = y1_tl + (y2_tl-y1_tl)*(x - x1)/(x2 - x1)
   END SUBROUTINE interp_linear_f1_tl


   SUBROUTINE interp_linear_f1_ad(x1, x2, x, y_AD, y1_AD, y2_AD)
     REAL(fp), INTENT(IN)     :: x1, x2, x
     REAL(fp), INTENT(IN OUT) :: y_AD
     REAL(fp), INTENT(IN OUT) :: y1_AD, y2_AD
     ! Local variables
     REAL(fp) :: fac
     fac = (x - x1)/(x2 - x1)
     y1_ad = y1_ad + (one - fac)*y_ad
     y2_ad = y2_ad + fac*y_ad
     y_ad = zero
   END SUBROUTINE interp_linear_f1_ad


 !--------------------------------------------------------------------------------
 !
 ! NAME:
 !       LayerAvg
 !
 ! PURPOSE:
 !    Given px1 (output domain) and px2 (input domain) that are ascending arrays,
 !    it computes the accumulated weighting factors for interpolation from input
 !    to output domainsan, and the index array, interp_index, such that output
 !    variable
 !      to(jo) = sum(pz(interp_index(1,jo):interp_index(2,jo),jo) &
 !                   *(ti(interp_index(1,jo):interp_index(2,jo))))

 !
 ! CALLING SEQUENCE:
 !       CALL LayerAvg(PX1,PX2,PZ,Interp_index)
 !
 ! INPUT ARGUMENTS:
 !       PX1:          The abscissa values for the target data (output domain) and
 !                     they must be monotonically ascending (e.g. lnP; in increasing values).
 !                     UNITS:      N/A
 !                     TYPE:       fp
 !                     DIMENSION:  rank-1 (KN1)
 !                     ATTRIBUTES: INTENT(IN)
 !       PX2:          The abscissa values for the source data (input domain) and
 !                     they must be monotonically ascending (e.g. lnP; in increasing values).
 !                     UNITS:      N/A
 !                     TYPE:       fp
 !                     DIMENSION:  rank-1 (KN2)
 !                     ATTRIBUTES: INTENT(IN)
 !
 ! OUTPUT ARGUMENTS:
 !       PZ:           Resultant accumulated weighting factors for
 !                     interpolation from input to output domains.
 !                     UNITS:      N/A
 !                     TYPE:       fp
 !                     DIMENSION:  rank-2 (KN2,KN1)
 !                     ATTRIBUTES: INTENT(OUT)
 !    interp_index:    The index array of dimension (2 x KN1) for
 !                     Start index for relevant PZ row array segment and
 !                     End index for relevant PZ row array segment, where KN1 = SIZE(PX1)
 !                     UNITS:      N/A
 !                     TYPE:       Integer
 !                     DIMENSION:  rank-2
 !                     ATTRIBUTES: INTENT(IN)
 !
 !     Comments:
 !
 !     1) PX1(i)<PX1(i+1) & PX2(i)<PX2(i+1)
 !
 ! RESTRICTIONS:
 !     To be efficient, this routine does not check that x and u are both
 !     monotonically ascending and the index bounds.
 !
 !     - Journal reference:
 !       Rochon, Y.J., L. Garand, D.S. Turner, and S. Polavarapu.
 !       Jacobian mapping between vertical coordinate systems in data assimilation,
 !       Q. J. of the Royal Met. Soc., 133, 1547-1558 (2007) DOI: 10.1002/qj.117
 !
 ! CREATION HISTORY:
 !       Written by:     Yong Chen, 08-May-2009
 !-----------------------------------------------------------------------------------
 !--------------------------------------------------------------------------------
   SUBROUTINE layeravg( PX1,PX2,PZ,interp_index)
     REAL(fp),     DIMENSION(:),       INTENT(IN)     :: PX1(:)
     REAL(fp),     DIMENSION(:),       INTENT(IN)     :: PX2(:)
     REAL(fp),     DIMENSION(:,:),     INTENT(OUT)    :: PZ(:, :)
     INTEGER,      DIMENSION(:,:),     INTENT(OUT)    :: interp_index

     ! Local variables
     INTEGER  ::  KN1,KN2, ibot,itop,ii,istart,iend, J,IC,ISKIP,KI
     REAL(fp) ::  z1,z2,z3,zw1,zw2,zsum
     REAL(fp) ::  y1,y2,d,w10,w20,dz,dx,dy,dzd,dxd

     kn1 = SIZE(px1)
     kn2 = SIZE(px2)

     istart=1
     iend=kn1
     DO ki=1,kn1
       z2=px1(ki)
 !
       if (ki == 1) then
          z1=2.0*z2-px1(ki+1)
       else
          z1=px1(ki-1)
       endif
 !
       if (ki == kn1) then
          z3=2.0*z2-z1
       else
          z3=px1(ki+1)
       endif
       if (z3 > px2(kn2)) z3=px2(kn2)
 !
       iskip=0
       if (z2 >= px2(kn2)) then
          z3=px2(kn2)
          z2=px2(kn2)
          iskip=1
       endif

 ! --- Determine forward interpolator
 !
       pz(1:kn2,ki)=zero
       ic=0
       do j=istart,kn2-1
         if (px2(j) > z3) go to 1000
 !
         if (px2(j) <= z2 .and. px2(j+1) > z1) then
           itop=0
           ibot=0
           if (z1 < z3) then
              y1=z1
              if (px2(j) > z1) then
                 y1=px2(j)
                 itop=1
              endif
              y2=z2
              if (px2(j+1) < z2) then
                 y2=px2(j+1)
                 ibot=1
              endif
           else
              y1=z2
              if (px2(j) > z2) then
                 y1=px2(j)
                 itop=1
              endif
              y2=z1
              if (px2(j+1) < z1) then
                 y2=px2(j+1)
                 ibot=1
              endif
           endif
 !
 ! ---     Set weights for forward interpolator
 !
           dy=y2-y1
           dz=z1-z2
           if (abs(dz) < smalldiff) then
              write(6,*) 'SUBLAYER: ERROR: dz is <=0. dz = ',dz
              write(6,*) 'z1,z2,z3 = ',z1,z2,z3
              write(6,*) 'px2(j),px2(j+1)    = ',px2(j),px2(j+1)
              return
           else
              dzd=one/dz
           endif
           zw1=(z1-y1)*dzd*dy
           zw2=(z1-y2)*dzd*dy
           w10=zw1
           w20=zw2
           dx=(px2(j+1)-px2(j))
           if (abs(dx) < smalldiff) then
              write(6,*) 'SUBLAYER: ERROR: dx is <=0. dx = ',dx
              write(6,*) 'z1,z2,z3 = ',z1,z2,z3
              write(6,*) 'px2(j),px2(j+1)    = ',px2(j),px2(j+1)
              return
           else
              dxd=one/dx
           endif
 !
           d=(px2(j+1)-z2)*dxd
           if (z1 < z3 .and. ibot == 0) then
              zw1=zw1+zw2*d
              zw2=zw2*(one-d)
           else if (z1 > z3 .and. itop == 0) then
              zw2=zw2+zw1*(one-d)
              zw1=zw1*d
           end if
           pz(j,ki)=pz(j,ki)+zw1
           pz(j+1,ki)=pz(j+1,ki)+zw2
           ic=1
         endif
 !
         if (px2(j) < z3 .and. px2(j+1) >= z2 .and. iskip == 0) then
           itop=0
           ibot=0
           if (z3 < z1) then
              y1=z3
              if (px2(j) > z3) then
                 y1=px2(j)
                 itop=1
              endif
              y2=z2
              if (px2(j+1) < z2) then
                 y2=px2(j+1)
                 ibot=1
              endif
           else
              y1=z2
              if (px2(j) > z2) then
                 y1=px2(j)
                 itop=1
              endif
              y2=z3
              if (px2(j+1) < z3) then
                 y2=px2(j+1)
                 ibot=1
              endif
           endif
 !
 ! ---     Set weights for forward interpolator
 !
           dy=y2-y1
           dz=z3-z2
           if (abs(dz) < smalldiff) then
              write(6,*) 'SUBLAYER: ERROR: dz is <=0. dz = ',dz
              write(6,*) 'z3,z2,z1 = ',z3,z2,z1
              write(6,*) 'px2(j),px2(j+1)    = ',px2(j),px2(j+1)
              return
           else
              dzd=one/dz
           endif
           zw1=(z3-y1)*dzd*dy
           zw2=(z3-y2)*dzd*dy
           w10=zw1
           w20=zw2
           dx=(px2(j+1)-px2(j))
           if (abs(dx) < smalldiff) then
              write(6,*) 'SUBLAYER: ERROR: dx is <=0. dx = ',dx
              write(6,*) 'z3,z2,z1 = ',z3,z2,z1
              write(6,*) 'px2(j),px2(j+1)    = ',px2(j),px2(j+1)
              return
           else
              dxd=one/dx
           endif
 !
           d=(px2(j+1)-z2)*dxd
           if (z3 < z1 .and. ibot == 0) then
              zw1=zw1+zw2*d
              zw2=zw2*(one-d)
           else if (z3 > z1 .and. itop == 0) then
              zw2=zw2+zw1*(one-d)
              zw1=zw1*d
           end if
           pz(j,ki)=pz(j,ki)+zw1
           pz(j+1,ki)=pz(j+1,ki)+zw2
           ic=1
         endif
       enddo
       j=kn2
  1000 continue
       if (ic == 0) pz(j,ki)=one
 !
 !     Normalize sum to unity (instead of calculating and dividing by
 !     weighting denominator)
 !
       do ii=istart,kn2
         if(pz(ii,ki) /= zero)then
           interp_index(1, ki)=ii
           exit
         endif
       enddo
       istart=interp_index(1, ki)

       interp_index(2,ki)=kn2

       do ii=interp_index(1, ki)+1,kn2
         if(pz(ii,ki) == zero)then
           interp_index(2,ki)=ii-1
           exit
         endif
       enddo
       iend=interp_index(2,ki)

       zsum=sum(pz(interp_index(1, ki):interp_index(2,ki),ki))
       pz(interp_index(1,ki):interp_index(2,ki),ki)= &
                  pz(interp_index(1,ki):interp_index(2, ki),ki)/zsum
     ENDDO


   END SUBROUTINE layeravg

 END MODULE odps_coordinatemapping
odps_coordinatemapping::cc
real(fp), parameter cc
Definition: ODPS_CoordinateMapping.f90:51

odps_coordinatemapping::geopotential_height
subroutine geopotential_height(Level_Pressure, Temperature, Water_Vapor, Surface_Height, Level_Height)
Definition: ODPS_CoordinateMapping.f90:614

crtm_parameters::zero
real(fp), parameter, public zero
Definition: CRTM_Parameters.f90:70

type_kinds::fp
integer, parameter, public fp
Definition: Type_Kinds.f90:124

odps_coordinatemapping::module_rcs_id
character(*), parameter, private module_rcs_id
Definition: ODPS_CoordinateMapping.f90:59

profile_utility_parameters::r_dryair
real(fp), parameter, public r_dryair
Definition: Profile_Utility_Parameters.f90:196

crtm_parameters
Definition: CRTM_Parameters.f90:12

crtm_atmosphere_define::h2o_id
integer, parameter, public h2o_id
Definition: CRTM_Atmosphere_Define.f90:250

crtm_geometryinfo_define::crtm_geometryinfo_type
Definition: CRTM_GeometryInfo_Define.f90:104

odps_coordinatemapping::layeravg
subroutine layeravg(PX1, PX2, PZ, interp_index)
Definition: ODPS_CoordinateMapping.f90:1111

crtm_parameters::earth_radius
real(fp), parameter, public earth_radius
Definition: CRTM_Parameters.f90:181

odps_coordinatemapping::c
real(fp), parameter c
Definition: ODPS_CoordinateMapping.f90:49

odps_coordinatemapping::map_input_tl
subroutine, public map_input_tl(TC, Atm_TL, Temperature_TL, Absorber_TL, PAFV)
Definition: ODPS_CoordinateMapping.f90:379

odps_coordinatemapping::interp_linear_f1_tl
subroutine interp_linear_f1_tl(x1, x2, x, y1_TL, y2_TL, y_TL)
Definition: ODPS_CoordinateMapping.f90:1028

odps_coordinatemapping::smalldiff
real(fp), parameter smalldiff
Definition: ODPS_CoordinateMapping.f90:53

crtm_parameters::one
real(fp), parameter, public one
Definition: CRTM_Parameters.f90:71

odps_coordinatemapping::interpolate_profile_f1_tl
subroutine, public interpolate_profile_f1_tl(interp_index, x, u, y_TL, y_int_TL)
Definition: ODPS_CoordinateMapping.f90:800

pafv_define::pafv_type
Definition: PAFV_Define.f90:86

pafv_define
Definition: PAFV_Define.f90:11

odps_coordinatemapping::map_input_ad
subroutine, public map_input_ad(TC, Temperature_AD, Absorber_AD, Atm_AD, PAFV)
Definition: ODPS_CoordinateMapping.f90:487

pafv_define::pafv_associated
elemental logical function, public pafv_associated(self)
Definition: PAFV_Define.f90:223

profile_utility_parameters
Definition: Profile_Utility_Parameters.f90:12

crtm_parameters::two
real(fp), parameter, public two
Definition: CRTM_Parameters.f90:72

odps_coordinatemapping::compute_interp_index
subroutine, public compute_interp_index(x, u, interp_index)
Definition: ODPS_CoordinateMapping.f90:966

odps_coordinatemapping::interp_linear
subroutine interp_linear(y1, x1, y2, x2, x, y)
Definition: ODPS_CoordinateMapping.f90:1021

odps_coordinatemapping
Definition: ODPS_CoordinateMapping.f90:12

odps_define::odps_type
Definition: ODPS_Define.f90:112

odps_coordinatemapping::map_input
subroutine, public map_input(Atm, TC, GeoInfo, Temperature, Absorber, User_Level_LnPressure, Ref_Level_LnPressure, Secant_Zenith, H2O_idx, PAFV)
Definition: ODPS_CoordinateMapping.f90:163

crtm_parameters::degrees_to_radians
real(fp), parameter, public degrees_to_radians
Definition: CRTM_Parameters.f90:88

type_kinds
Definition: Type_Kinds.f90:12

profile_utility_parameters::eps
real(fp), parameter, public eps
Definition: Profile_Utility_Parameters.f90:193

odps_coordinatemapping::interpolate_profile
subroutine, public interpolate_profile(interp_index, y, x, u, y_int)
Definition: ODPS_CoordinateMapping.f90:706

crtm_atmosphere_define::crtm_atmosphere_type
Definition: CRTM_Atmosphere_Define.f90:361

crtm_atmosphere_define
Definition: CRTM_Atmosphere_Define.f90:13

odps_coordinatemapping::interpolate_profile_f1_ad
subroutine, public interpolate_profile_f1_ad(interp_index, x, u, y_int_AD, y_AD)
Definition: ODPS_CoordinateMapping.f90:895

crtm_geometryinfo_define
Definition: CRTM_GeometryInfo_Define.f90:13

odps_coordinatemapping::interp_linear_f1_ad
subroutine interp_linear_f1_ad(x1, x2, x, y_AD, y1_AD, y2_AD)
Definition: ODPS_CoordinateMapping.f90:1035

odps_define
Definition: ODPS_Define.f90:13

crtm_geometryinfo_define::crtm_geometryinfo_getvalue
elemental subroutine, public crtm_geometryinfo_getvalue(gInfo, Geometry, iFOV, Longitude, Latitude, Surface_Altitude, Sensor_Scan_Angle, Sensor_Zenith_Angle, Sensor_Azimuth_Angle, Source_Zenith_Angle, Source_Azimuth_Angle, Flux_Zenith_Angle, Year, Month, Day, Distance_Ratio, Sensor_Scan_Radian, Sensor_Zenith_Radian, Sensor_Azimuth_Radian, Secant_Sensor_Zenith, Cosine_Sensor_Zenith, Source_Zenith_Radian, Source_Azimuth_Radian, Secant_Source_Zenith, Flux_Zenith_Radian, Secant_Flux_Zenith, Trans_Zenith_Radian, Secant_Trans_Zenith, AU_ratio2)
Definition: CRTM_GeometryInfo_Define.f90:421