! $Id: acha_cloud_cover_layers_module.f90 1523 2016-02-22 16:30:27Z wstraka $
!--------------------------------------------------------------------------------------
! Clouds from AVHRR Extended (CLAVR-x) 1b PROCESSING SOFTWARE Version 6.0
!
! NAME: cloud_cover_layers_module.f90 (src)
!       CLOUD_COVER_LAYERS (module)
!
! PURPOSE: this module houses routines for computing cloud cover or
!          fraction in each atmospheric layer including the total
!
! DESCRIPTION: 
!
! AUTHORS:
!  Andrew Heidinger, Andrew.Heidinger@noaa.gov
!
! COPYRIGHT
! THIS SOFTWARE AND ITS DOCUMENTATION ARE CONSIDERED TO BE IN THE PUBLIC
! DOMAIN AND THUS ARE AVAILABLE FOR UNRESTRICTED PUBLIC USE. THEY ARE
! FURNISHED "AS IS." THE AUTHORS, THE UNITED STATES GOVERNMENT, ITS
! INSTRUMENTALITIES, OFFICERS, EMPLOYEES, AND AGENTS MAKE NO WARRANTY,
! EXPRESS OR IMPLIED, AS TO THE USEFULNESS OF THE SOFTWARE AND
! DOCUMENTATION FOR ANY PURPOSE. THEY ASSUME NO RESPONSIBILITY (1) FOR
! THE USE OF THE SOFTWARE AND DOCUMENTATION; OR (2) TO PROVIDE TECHNICAL
! SUPPORT TO USERS.
!
! Public routines used in this MODULE:
!
! Cloud Layer Definition
!
! Layered Cloud Layer Meanings
! 00000000 = Clear
! 00000001 = First Layer
! 00000010 = Second Layer
! 00000100 = Third Layer
! 00001000 = Fourth Layer
! 00010000 = Fifth Layer
! 00000101 = Third and First
! 00000110 = Third and Second
! 00000011 = Second and First
! 00000111 = High and Mid and Low
! etc, etc, etc
!
! to extract
! to see if first layer cloud, ibits(layer,0,1)
! to see if second layer cloud, ibits(layer,1,1)
! to see if third layer cloud, ibits(layer,2,1)
!
! Cloud_Cover_Layer_Profile = one bit integer = (0 to 100)
! Supercooled_Profile = one bit integer = (0 to 100)
!
!--------------------------------------------------------------------------------------
module CCL_MODULE

  use CCL_SERVICES_MOD, only : &
           real4, int1, int4, ccl_output_struct,ccl_symbol_struct, &
           ccl_input_struct, ccl_diag_struct

 implicit none
 public:: SETUP_CCL
 public:: DESTROY_CCL
 public:: COMPUTE_CLOUD_COVER_LAYERS
 private:: COMPUTE_BOX_WIDTH

 type(ccl_symbol_struct), private :: symbol

 integer, parameter:: N_Levels_NCEP = 4
 integer, parameter:: N_Levels_ISCCP = 4
 integer, parameter:: N_Levels_NOAT = 6
 real, dimension(N_Levels_NCEP), parameter, private:: CCL_Levels_NCEP =   (/1100.0,631.0,350.0,0.0/)  !hpa
 real, dimension(N_Levels_ISCCP), parameter, private:: CCL_Levels_ISCCP = (/1100.0,680.0,440.0,0.0/)  !hpa
 real, dimension(N_Levels_NOAT), parameter, private:: CCL_Levels_NOAT =   (/0.0,5000.0,10000.0,18000.0,24000.0,99000.0/) !kft
 integer, save, private:: N_Levels
 integer, save, private:: N_Layers
 real, dimension(:), allocatable, save, private:: CCL_Levels
  integer:: Elem_Idx
  integer:: Line_Idx
  integer:: Num_Layer
!  integer:: Layer_Idx, Layer_Top_Idx, Layer_Lower_Idx,Layer_Base_Idx
!  real:: Z, Z_Base, Z_Lower
!ynoh (cira/csu) for ccl mode 3
  integer:: Layer_Idx, Layer_Top_Idx, Layer_Lower_Idx, Layer_Base_Idx, Layer_Lower_Base_Idx                            
  real:: Z, Z_Base, Z_Lower, Z_Lower_Base , Z_Lower_dummy

 include 'ccl_parameters.inc'

 contains
!------------------------------------------------------------------------------
! Setup CCL parameters based on CCL_TYPE
!------------------------------------------------------------------------------
subroutine SETUP_CCL(CCL_Type) 

   integer, intent(in) :: CCL_Type 

   select case (CCL_Type)

    case (2)    !NCEP
      N_Levels = N_Levels_NCEP
      allocate(CCL_Levels(N_Levels))
      CCL_Levels = CCL_Levels_NCEP

    case (1)    !ISCCP
      N_Levels = N_Levels_ISCCP
      allocate(CCL_Levels(N_Levels))
      CCL_Levels = CCL_Levels_ISCCP

    case (0)  !NOAT
      N_Levels = N_Levels_NOAT
      allocate(CCL_Levels(N_Levels))
      CCL_Levels = CCL_Levels_NOAT

  end select 

  N_Layers = N_Levels - 1

end subroutine SETUP_CCL

!------------------------------------------------------------------------------
! Destroy CCL parameters based on CCL_TYPE
!------------------------------------------------------------------------------
subroutine DESTROY_CCL() 

  if (allocated(CCL_Levels)) deallocate(CCL_Levels)

end subroutine DESTROY_CCL

!------------------------------------------------------------------------------
! compute cloud fraction over a nxn array using the Bayesian probability
!------------------------------------------------------------------------------
 subroutine COMPUTE_CLOUD_COVER_LAYERS(Input, Symbol_In, Output, Diag)

  !===============================================================================
  !  Argument Declaration
  !==============================================================================

  type(ccl_input_struct), intent(inout) :: Input
  type(ccl_symbol_struct), intent(inout) :: Symbol_In
  type(ccl_output_struct), intent(inout) :: Output
  type(ccl_diag_struct), intent(inout), optional :: Diag
  integer, save:: Diag_Warning_Flag = 0


  integer:: Num_Elems
  integer:: Num_Lines
  integer :: i,j                    !pixel indices
  integer :: i1,i2,j1,j2            !pixel indices of ccl box
  integer :: i11,i22,j11,j22        !pixels inices of area which ccl result is valid
  integer :: ni, nj                 !size of ccl box
  integer:: N                       !width of ccl box
  integer:: M                       !width of ccl region
  integer :: Num_Good
  integer :: Num_Cloud
  integer:: Num_High, Num_Mid, Num_Low, Num_Clear, Num_All
  real:: Clear_Fraction

  real (kind=4), dimension(:,:), allocatable:: Pixel_Uncertainty
  integer (kind=1), dimension(:,:), allocatable:: pos, len      !ibits arguments

  real (kind=real4), dimension(:,:), allocatable:: Alt
  real (kind=real4), dimension(:,:), allocatable:: Alt_Lower
  real (kind=real4), dimension(:,:), allocatable:: Alt_Base
!ynoh (cira/csu) for ccl mode 3
  real (kind=real4), dimension(:,:), allocatable:: Alt_Lower_Base

!  logical:: USE_BASE, USE_LOWER
  logical:: USE_BASE, USE_LOWER, USE_LOWER_BASE


  USE_BASE = .false.
  USE_LOWER = .false.
  if (Input%CCL_Mode >= 2) USE_BASE = .true.
  if (Input%CCL_Mode == 3) USE_LOWER = .true.
!ynoh (cira/csu) for ccl mode 3  (to use the lower base)
  USE_LOWER_BASE = .true.
  !-------------------------------------------------------------------------------
  ! Total Cloud Fraction and Its Uncertainty
  !-------------------------------------------------------------------------------

  !--- copy input symbol to a module-wide variable
  symbol = Symbol_In

  !--- initialize
  Output%Total_Cloud_Fraction = MISSING_VALUE_REAL4
  Output%Total_Cloud_Fraction_Uncer = MISSING_VALUE_REAL4
  Output%Cloud_Layer = MISSING_VALUE_INTEGER1
  Output%QF = MISSING_VALUE_INTEGER1

  !--- initialize diagnostic output
  if (present(Diag) .and. Diag_Warning_Flag == 0) then
      print *, "CLAVR-x / CCL ===>  Diagnostic Output Turned On"
      Diag_Warning_Flag = 1
  endif
  if (present(Diag)) Diag%Array_1 = Missing_Value_Real4
  if (present(Diag)) Diag%Array_2 = Missing_Value_Real4
  if (present(Diag)) Diag%Array_3 = Missing_Value_Real4

  !--- Determine Box Width 
  call COMPUTE_BOX_WIDTH(Input%Sensor_Resolution_KM,CCL_BOX_WIDTH_KM, N)
  call COMPUTE_BOX_WIDTH(Input%Sensor_Resolution_KM,CCL_SPACING_KM, M)

!==============================================================================
!
! Compute Cloud Cover Layers which is defined here as the fraction of
! high, middle and low-level cloud in a NxN array centered on each pixel
!
! output (through global arrays)
!  Output%Cloud_Fraction = cloud fraction (0-1) for all clouds over array
!  Output%Cloud_Fraction_Uncertainty = cloud fraction uncertainty (0-1)
!  Output%Cloud_Cover_Layer= cloud fraction (0-100) (integer)
!
!  Output%High_Cloud_Fraction = cloud fraction (0-1) for high cloud over array
!  Output%Mid_Cloud_Fraction = cloud fraction (0-1) for mid cloud over array
!  Output%Low_Cloud_Fraction = cloud fraction (0-1) for low cloud over array
!
!  Notes
!  1. cloud fractions have values limited by 1/N^2. If N=1, values are 1/9,2/9...
!  2. Cloud_Fraction is the total cloud amount.  This is computed from the
!     mask only.  If arrays have pixels with failed ACHA, H+M+L /= Total.
!  3. H/M/L boundaries in ccl_parameters.inc
!==============================================================================

  !--------------------------------------------------------------------
  ! compute pixel-level cloud layer flag and H/M/L masks
  !--------------------------------------------------------------------
  Num_Elems = Input%Number_of_Elements
  Num_Lines = Input%Number_of_Lines
  allocate(Pixel_Uncertainty(Num_Elems, Num_Lines))

  !--- make cloud fraction pixel level uncertainty
  Pixel_Uncertainty = MISSING_VALUE_REAL4
  where(Input%Cloud_Probability >= 0.5)
     Pixel_Uncertainty = 1.0 - Input%Cloud_Probability
  endwhere
  where(Input%Cloud_Probability < 0.5)
     Pixel_Uncertainty = Input%Cloud_Probability
  endwhere

!---- Alt should be an input
! call COMPUTE_ALTITUDE_FROM_PRESSURE_CCL(1,Num_Lines,Num_Elems,symbol%YES,Input%Invalid_Data_Mask,Input%Pc,Alt)
! call COMPUTE_ALTITUDE_FROM_PRESSURE_CCL(1,Num_Lines,Num_Elems,symbol%YES,Input%Invalid_Data_Mask,Input%Pc_Base,Alt_Base)
! call COMPUTE_ALTITUDE_FROM_PRESSURE_CCL(1,Num_Lines,Num_Elems,symbol%YES,Input%Invalid_Data_Mask,Input%Pc_Lower,Alt_Lower)

!ynoh (cira/csu) for ccl mode 3
! Previously ACHA%Alt and ACHA%Base_Alt are only computed in process_clavrx.f90

  !--------------------------------------------------------------------
  !Layer_Top_Idx <= N_Layers Compute cloud_layer flag 
  !--------------------------------------------------------------------
  line_loop: do Line_Idx = 1, Num_Lines
    elem_loop: do Elem_Idx = 1, Num_Elems

      !--- select the vertical coordinates based on CCL_Type
      select case (Input%CCL_Type) 

        case (1,2)
          Z = Input%Pc(Elem_Idx,Line_Idx)
          Z_Base = Input%Pc_Base(Elem_Idx,Line_Idx)
          Z_Lower = Input%Pc_Lower(Elem_Idx,Line_Idx)
!ynoh (cira/csu) for ccl mode 3
          Z_Lower_Base = Input%Pc_Lower_Base(Elem_Idx,Line_Idx)
          Z_Lower_dummy = max(Z_Lower, Z_Base)

        case (0)
          Z = Input%Alt(Elem_Idx,Line_Idx)
          Z_Base = Input%Alt_Base(Elem_Idx,Line_Idx)
!ynoh (cira/csu) for ccl mode 3  (three added)
          Z_Lower = Input%Alt_Lower(Elem_Idx,Line_Idx)
          Z_Lower_Base = Input%Alt_Lower_Base(Elem_Idx,Line_Idx)
          Z_Lower_dummy = min(Z_Lower, Z_Base)

      end select 

      if (Z /= MISSING_VALUE_REAL4) Z = max(0.0,Z)
      if (Z_Base /= MISSING_VALUE_REAL4) Z_Base = max(0.0,Z_Base)
      if (Z_Lower /= MISSING_VALUE_REAL4) Z_Lower = max(0.0,Z_Lower)
!ynoh (cira/csu) for ccl mode 3
      if (Z_Lower_Base /= MISSING_VALUE_REAL4) Z_Lower_Base = max(0.0,Z_Lower_Base)

      Layer_Top_Idx = -1
      Layer_Base_Idx = -1
      Layer_Lower_Idx = -1
!ynoh (cira/csu) for ccl mode 3
      Layer_Lower_Base_Idx = -1

      !----- check for valid cloud vertical coordinate
      if (Z == MISSING_VALUE_REAL4) then
          cycle
      endif

      !-- need to intialize pixels with valid data to 0
      Output%Cloud_Layer(Elem_Idx,Line_Idx) = 0

      !---------------------------------------------------------------------
      ! Find Layers
      !---------------------------------------------------------------------

      !--- find top of highest cloud level index
      if (Z /= MISSING_VALUE_REAL4) then
         do Layer_Top_Idx = 1, N_Layers
            if ((Z >= CCL_Levels(Layer_Top_Idx) .and. Z <= CCL_Levels(Layer_Top_Idx + 1)) .or. &
                (Z <= CCL_Levels(Layer_Top_Idx) .and. Z >= CCL_Levels(Layer_Top_Idx + 1))) exit
         enddo
      endif

      !--- find base of highest cloud level index
      if (USE_BASE .and. Z_Base /= MISSING_VALUE_REAL4) then
        do Layer_Base_Idx = 1, N_Layers
         if ((Z_Base >= CCL_Levels(Layer_Base_Idx) .and. Z_Base <= CCL_Levels(Layer_Base_Idx + 1)) .or. &
             (Z_Base <= CCL_Levels(Layer_Base_Idx) .and. Z_Base >= CCL_Levels(Layer_Base_Idx + 1))) exit
        enddo
      endif

      !--- find top of lowest cloud level index
!ynoh (cira/csu) for ccl mode 3
!-- use Z_Lower_dummy instead of Z_Lower to avoid double counting layers in case Z_Lower > Z_Base 
      if (USE_LOWER .and. Z_Lower /= MISSING_VALUE_REAL4 .and. Z_Lower_dummy /= MISSING_VALUE_REAL4 ) then
        do Layer_Lower_Idx = 1, N_Layers
         if ((Z_Lower_dummy >= CCL_Levels(Layer_Lower_Idx) .and. Z_Lower_dummy <= CCL_Levels(Layer_Lower_Idx + 1)) .or. &
             (Z_Lower_dummy <= CCL_Levels(Layer_Lower_Idx) .and. Z_Lower_dummy >= CCL_Levels(Layer_Lower_Idx + 1))) exit
        enddo
      endif

!ynoh (cira/csu) for ccl mode 3
      !--- find base of lowest cloud level index 
      if (USE_LOWER .and. USE_LOWER_BASE .and. &
          Z_Lower_Base /= MISSING_VALUE_REAL4 .and. Z_Lower_Base < Z_Lower_dummy ) then
        do Layer_Lower_Base_Idx = 1, N_Layers
         if ((Z_Lower_Base >= CCL_Levels(Layer_Lower_Base_Idx) .and. Z_Lower_Base <= CCL_Levels(Layer_Lower_Base_Idx + 1)) .or. &
             (Z_Lower_Base <= CCL_Levels(Layer_Lower_Base_Idx) .and. Z_Lower_Base >= CCL_Levels(Layer_Lower_Base_Idx + 1))) exit
        enddo
      endif

      !---------------------------------------------------------------------
      ! Compute Cloud Layer
      !---------------------------------------------------------------------
      if (Layer_Top_Idx > 0 .and. Layer_Top_Idx <= N_Layers) then
          Output%Cloud_Layer(Elem_Idx,Line_Idx) = ibset(Output%Cloud_Layer(Elem_Idx,Line_Idx),Layer_Top_Idx-1)
      endif

      !--- account for base
      do Layer_Idx = 1, N_Layers
        if (USE_BASE .and. Layer_Base_Idx <= Layer_Top_Idx .and. Layer_Idx >= Layer_Base_Idx .and. Layer_Top_Idx >= Layer_Idx) then
           Output%Cloud_Layer(Elem_Idx,Line_Idx) = ibset(Output%Cloud_Layer(Elem_Idx,Line_Idx),Layer_Idx-1)
        endif
      enddo 

      !--- account for lower cloud
!      if (USE_LOWER .and. Layer_Lower_Idx > 0 .and. Layer_Lower_Idx <= N_Layers) then
!ynoh (cira/csu) for ccl mode 3
!-- to avoid double counting layers in case Z_Lower > Z_Base 
      if (USE_LOWER .and. Layer_Lower_Idx > 0 .and. Layer_Lower_Idx <= N_Layers .and. Layer_Lower_Idx < Layer_Base_Idx) then
          Output%Cloud_Layer(Elem_Idx,Line_Idx) = ibset(Output%Cloud_Layer(Elem_Idx,Line_Idx),Layer_Lower_Idx-1)
      endif

!ynoh (cira/csu) for ccl mode 3
      !--- account for lower cloud base
      do Layer_Idx = 1, N_Layers
        if (USE_LOWER .and. USE_LOWER_BASE .and. Layer_Lower_Base_Idx <= Layer_Lower_Idx .and. &
            Layer_Idx >= Layer_Lower_Base_Idx .and. &
            Layer_Lower_Idx >= Layer_Idx .and. Layer_Idx < Layer_Base_Idx) then
          Output%Cloud_Layer(Elem_Idx,Line_Idx) = ibset(Output%Cloud_Layer(Elem_Idx,Line_Idx),Layer_Lower_Idx-1)
        endif
      enddo 

      !print *, "Test ", Z, Layer_Top_Idx, Z_Base, Layer_Base_Idx, Output%Cloud_Layer(Elem_Idx,Line_Idx)

!      write(unit=6,fmt="(I3, I3, I08,B08)") Layer_Top_Idx, Layer_Base_Idx, Output%Cloud_Layer(Elem_Idx,Line_Idx), Output%Cloud_Layer(Elem_Idx,Line_Idx)

    end do elem_loop
  end do line_loop

   where (Input%Cloud_Mask == Symbol%CLEAR .or. Input%Cloud_Mask == Symbol%PROB_CLEAR)
          Output%Cloud_Layer = 0   !0000
   endwhere 

! set QF based on ACHA QF
   where (Input%ACHA_QF == CTH_DQF_GOOD_RETREVIAL ) Output%QF = CCL_DQF_GOOD_RETREVIAL
   where (Input%ACHA_QF == CTH_DQF_MARGINAL_RETREVIAL .or. Input%ACHA_Qf == CTH_DQF_RETREVIAL_ATTEMPTED) &
         Output%QF = CCL_DQF_DEGRADED_RETREVIAL
   where (Input%ACHA_QF == CTH_DQF_BAD_RETREVIAL ) Output%QF = CCL_DQF_BAD_RETREVIAL


 !--------------------------------------------------------------------
 ! compute pixel-level cloud cover for each layer over the box
 ! N = ccl box size
 ! M = ccl result spacing
 !--------------------------------------------------------------------

 line_loop_cover: do j = 1, Num_Lines, 2*M+1

    j1 = max(1,j-N)
    j2 = min(Num_Lines,j+N)
    j11 = max(1,j-M)
    j22 = min(Num_Lines,j+M)

    element_loop_cover: do i = 1, Num_Elems, 2*M+1

      i1 = max(1,i-N)
      i2 = min(Num_Elems,i+N)
      i11 = max(1,i-M)
      i22 = min(Num_Elems,i+M)

      ni = i2 - i1 + 1 
      nj = j2 - j1 + 1 

      if (allocated(pos)) deallocate(pos)
      if (allocated(len)) deallocate(len)

      allocate(pos(ni,nj), len(ni,nj))
      len = 1

      !--- check for a bad pixel pixel
      if (Input%Invalid_Data_Mask(i,j) == Symbol%YES) cycle

      !---- new   layers are inbetween levels
      Num_Clear = count(Input%Cloud_Probability(i1:i2,j1:j2) < 0.5 .and. &
                        Input%Cloud_Probability(i1:i2,j1:j2) /= MISSING_VALUE_REAL4)         
      Num_Cloud = count(Input%Cloud_Probability(i1:i2,j1:j2) >= 0.5)        
      Num_Good = count(Input%Cloud_Probability(i1:i2,j1:j2) /= MISSING_VALUE_REAL4)        
      Num_All = Num_Cloud + Num_Clear

      !--- see if there are any valid mask points, if not skip this pixel
      if (Num_Good < COUNT_MIN_CCL) then
         cycle
      endif

      !--- Total Cloud Fraction
      Output%Total_Cloud_Fraction(i11:i22,j11:j22) = real(Num_Cloud)/real(Num_Good)

      !--- compute the uncertainty of the total cloud fraction
      Output%Total_Cloud_Fraction_Uncer(i11:i22,j11:j22) = sum(Pixel_Uncertainty(i1:i2,j1:j2))/real(Num_Good)

      !--- see if there are any valid CCL points, if not skip this pixel
      if (Num_All == 0 ) then
         cycle
      endif

      do Layer_Idx = 1, N_Layers
          Num_Layer = int(sum(real(ibits(Output%Cloud_Layer(i1:i2,j1:j2),Layer_Idx-1,len))))
          if (Layer_Idx == 1) Output%Fraction_Layer1(i11:i22,j11:j22) = int(100.0*real(Num_Layer)/real(Num_All))
          if (Layer_Idx == 2) Output%Fraction_Layer2(i11:i22,j11:j22) = int(100.0*real(Num_Layer)/real(Num_All))
          if (Layer_Idx == 3) Output%Fraction_Layer3(i11:i22,j11:j22) = int(100.0*real(Num_Layer)/real(Num_All))
          if (Layer_Idx == 4) Output%Fraction_Layer4(i11:i22,j11:j22) = int(100.0*real(Num_Layer)/real(Num_All))
          if (Layer_Idx == 5) Output%Fraction_Layer5(i11:i22,j11:j22) = int(100.0*real(Num_Layer)/real(Num_All))
      enddo 

!write(unit=6,fmt="(I08,B08,5I5)") Output%Cloud_Layer(i,j), Output%Cloud_Layer(i,j), &
!                              Output%Fraction_Layer1(i,j), &
!                              Output%Fraction_Layer2(i,j), &
!                              Output%Fraction_Layer3(i,j), &
!                              Output%Fraction_Layer4(i,j), &
!!                             Output%Fraction_Layer5(i,j)
    end do element_loop_cover
 end do line_loop_cover

 if (allocated(Pixel_Uncertainty)) deallocate(Pixel_Uncertainty)


 end subroutine COMPUTE_CLOUD_COVER_LAYERS
!----------------------------------------------------------------------
!--- determine cirrus box width
!---
!--- Sensor_Resolution_KM = the nominal resolution in kilometers
!--- Box_Width_KM = the width of the desired box in kilometers
!--- Box_Half_Width = the half width of the box in pixel-space
!----------------------------------------------------------------------
subroutine COMPUTE_BOX_WIDTH(Sensor_Resolution_KM,Box_Width_KM, &
                             Box_Half_Width)

   real, intent(in):: Sensor_Resolution_KM
   integer, intent(in):: Box_Width_KM
   integer, intent(out):: Box_Half_Width

   if (Sensor_Resolution_KM <= 0.0) then
       Box_Half_Width = 20
   else
       Box_Half_Width = int((Box_Width_KM / Sensor_Resolution_KM) / 2)
   endif

end subroutine COMPUTE_BOX_WIDTH

!----------------------------------------------------------------------
! routine to interpolate pressure to flight level altitude.
!----------------------------------------------------------------------
subroutine COMPUTE_ALTITUDE_FROM_PRESSURE_CCL(Line_Idx_Min,Num_Lines,Num_Elems,symbol_yes,Bad_Pixel_Mask,Pc_In,Alt_Out)

!--- Based on Sarah Monette calculations for HS3.  Her calculation is based on:
!--- http://psas.pdx.edu/RocketScience/PressureAltitude_Derived.pdf.

  integer, intent(in):: Line_Idx_Min
  integer, intent(in):: Num_Lines
  integer, intent(in):: Num_Elems
  integer(kind=int1), intent(in) :: symbol_yes
  integer(kind=int1), intent(in), dimension(:,:):: Bad_Pixel_Mask
  real, intent(inout), dimension(:,:) :: Alt_Out
  real, intent(in), dimension(:,:) :: Pc_In
  integer:: Elem_Idx
  integer:: Line_Idx
  integer:: Line_Start
  integer:: Line_End
  real:: Pc_Temp
  real:: Alt_Temp

  !--- Constants from Sarah Monette
  real, parameter :: PW1 = 227.9 ! hPa
  real, parameter :: PW2 = 56.89 ! hPa
  real, parameter :: PW3 = 11.01 ! hPa
  real, parameter :: P0 = 1013.25 ! hPa
  real, parameter :: LR_OVER_G = 0.190263
  real, parameter :: Z0 = 145422.16 ! feet
  real, parameter :: LN_1 = -20859.0
  real, parameter :: LN_2 = 149255.0
  real, parameter :: PN_4 = 0.000470034
  real, parameter :: PN_3 = -0.364267
  real, parameter :: PN_2 = 47.5627
  real, parameter :: PN_1 = -2647.45
  real, parameter :: PN_0 = 123842.0

  !--- save elements
  Line_Start = Line_Idx_Min
  Line_End = Line_Start + Num_Lines - 1

  !--- initialize
  Alt_Out = Missing_Value_Real4

  !----------------------------------------------------------
  ! loop through segment
  !----------------------------------------------------------
  Line_Loop_1: do Line_Idx = Line_Start, Line_End
  Element_Loop_1: do Elem_Idx = 1, Num_Elems

     !--- Initialize temporary value each time.
     Pc_Temp = Missing_Value_Real4
     Alt_Temp = Missing_Value_Real4

     !--- skip bad pixels
     if (Bad_Pixel_Mask(Elem_Idx,Line_Idx) == symbol_yes) cycle

     !--- check if cld-top pressure is valid
     if (Pc_In(Elem_Idx,Line_Idx) == Missing_Value_Real4) cycle

     !--- Place valid pressure in temp variable for readability in the
     !--- calculations below.
     Pc_Temp = Pc_In(Elem_Idx,Line_Idx)

     !--- calculated altitude, in feet, from pressure.
     !--- 1st pivot point is directly from the pressure to
     !--- altitude from above reference.
     if (Pc_Temp > PW1) then
       Alt_Temp = (1.0 - (Pc_Temp/P0)**LR_OVER_G) * Z0
     endif

     !--- 2nd pivot point was modeled best with a natural log
     !--- fit.  From Sarah Monette.
     if (Pc_Temp <= PW1 .AND. Pc_Temp >= PW2) then
       Alt_Temp = LN_1 * LOG(Pc_Temp) + LN_2
     endif

     !--- 3rd pivot point. Modeled best with a polynomial
     !--- fit from Sarah Monette.
     if (Pc_Temp < PW2 .AND. Pc_Temp >= PW3) then
       Alt_Temp = (PN_4*Pc_Temp**4) + (PN_3*Pc_Temp**3) + (PN_2*Pc_Temp**2) + &
                  (PN_1*Pc_Temp) + PN_0
     endif

     if (Pc_Temp < PW3 .or. Alt_Temp < 0) then
       Alt_Temp = Missing_Value_Real4
     endif


     !--- Assign final altitude, in feet, to the level2 array.
     Alt_Out(Elem_Idx,Line_Idx) = Alt_Temp

  enddo Element_Loop_1
  enddo Line_Loop_1

end subroutine COMPUTE_ALTITUDE_FROM_PRESSURE_CCL

!-----------------------------------------------------------
! end of MODULE
!-----------------------------------------------------------
end module CCL_MODULE