! $Id$
!====================================================================
! SUBROUTINE Name: CONVERT_ATMOS_PROF_NWP_RTM
!
! Description:
! This routine interpolate the NWP profiles to profiles with the
! vertical spacing defined by the RTM model.  It operates on profiles
! stored in this module
!
! INPUTS:
!
! Highest_Level_Rtm_Nwp - highest Rtm Level that is below the highest nwp Level
! Lowest_Level_Rtm_Nwp - lowest Rtm Level that is above the lowest nwp Level
! Sfc_Level_Rtm - lowest Rtm Level above the surface
! P_near_Sfc_Nwp - lowest standard nwp Level above surface pressure
!
!====================================================================
 subroutine CONVERT_ATMOS_PROF_NWP_RTM ( )

    integer, intent(in):: Lon_Idx, Lat_Idx
    integer:: k, Lowest_Level_Rtm_Nwp, Highest_Level_Rtm_Nwp,  &
              Sfc_Level_Rtm
    real:: dT_dP_near_Sfc
    real:: dWvmr_dP_near_Sfc
    real:: dZ_dP_near_Sfc
    real:: P_near_Sfc_Nwp
    real:: Wvmr_Sfc
    real:: es
    real:: e
    real:: T_Offset


    !--- initialize indices
    Lowest_Level_Rtm_Nwp = NLevels_Rtm
    Highest_Level_Rtm_Nwp = 1
    Sfc_Level_Rtm = NLevels_Rtm
    P_Near_Sfc_Nwp = P_Std_Nwp(Sfc_Level_Nwp(Lon_Idx,Lat_Idx))
   
    !--- make Wvmr at sfc
    es = vapor(Tmpair_Nwp(Lon_Idx,Lat_Idx))
    e = Rhsfc_Nwp(Lon_Idx,Lat_Idx) * es / 100.0
    Wvmr_Sfc = 1000.0*0.622 * (e / (Psfc_Nwp(Lon_Idx,Lat_Idx) - e))  !(g/kg)

    !--- determine some critical Levels in the Rtm profile
    do k = 1, NLevels_Rtm
       if (P_Std_Rtm(k) > P_Std_Nwp(1)) then
            Highest_Level_Rtm_Nwp = k
            exit
       endif
    enddo

    do k = NLevels_Rtm,1,-1
       if (P_Std_Rtm(k) < Psfc_Nwp(Lon_Idx,Lat_Idx)) then
            Sfc_Level_Rtm = k
            exit
       endif
    enddo

    do k = NLevels_Rtm,1,-1
       if (P_Std_Rtm(k) < P_near_Sfc_Nwp) then
            Lowest_Level_Rtm_Nwp = k
            exit
       endif
    enddo

    !--- compute T and Wvmr lapse rate near the surface
    dT_dP_near_Sfc = 0.0
    dZ_dP_near_Sfc = 0.0
    dWvmr_dP_near_Sfc = 0.0
    if (Psfc_Nwp(Lon_Idx,Lat_Idx) /= P_Std_Nwp(NLevels_Nwp)) then
     dT_dP_near_Sfc =  &
             (T_Nwp(Sfc_Level_Nwp(Lon_Idx,Lat_Idx)) - Tmpair_Nwp(Lon_Idx,Lat_Idx))/ &
             (P_Std_Nwp(Sfc_Level_Nwp(Lon_Idx,Lat_Idx)) - Psfc_Nwp(Lon_Idx,Lat_Idx))
     dWvmr_dP_near_Sfc =  &
             (Wvmr_Nwp(Sfc_Level_Nwp(Lon_Idx,Lat_Idx)) - Wvmr_Sfc)/ &
             (P_Std_Nwp(Sfc_Level_Nwp(Lon_Idx,Lat_Idx)) - Psfc_Nwp(Lon_Idx,Lat_Idx))
     dZ_dP_near_Sfc =  &
             (Z_Nwp(Sfc_Level_Nwp(Lon_Idx,Lat_Idx)) - 0.0)/ &
             (P_Std_Nwp(Sfc_Level_Nwp(Lon_Idx,Lat_Idx)) - Psfc_Nwp(Lon_Idx,Lat_Idx))
    else
     dT_dP_near_Sfc =  &
             (T_Nwp(Sfc_Level_Nwp(Lon_Idx,Lat_Idx)-1) - Tmpair_Nwp(Lon_Idx,Lat_Idx))/ &
             (P_Std_Nwp(Sfc_Level_Nwp(Lon_Idx,Lat_Idx)-1) - Psfc_Nwp(Lon_Idx,Lat_Idx))
     dWvmr_dP_near_Sfc =  &
             (Wvmr_Nwp(Sfc_Level_Nwp(Lon_Idx,Lat_Idx)-1) - Wvmr_Sfc)/ &
             (P_Std_Nwp(Sfc_Level_Nwp(Lon_Idx,Lat_Idx)-1) - Psfc_Nwp(Lon_Idx,Lat_Idx))
     dZ_dP_near_Sfc =  &
             (Z_Nwp(Sfc_Level_Nwp(Lon_Idx,Lat_Idx)-1) - 0.0)/ &
             (P_Std_Nwp(Sfc_Level_Nwp(Lon_Idx,Lat_Idx)-1) - Psfc_Nwp(Lon_Idx,Lat_Idx))
    endif

    !--- compute temperature offset between standard and NWP profiles at top
    !--- this will be added to the standard profiles
    T_Offset = T_Nwp(1) - T_Std_Rtm(Highest_Level_Rtm_Nwp)

    !--- for Rtm Levels above the highest nwp Level, use Rtm standard values
    do k = 1,Highest_Level_Rtm_Nwp-1
          Z_Rtm(k) = Z_Nwp(1)
          T_Prof_Rtm(k) = T_Std_Rtm(k) + T_Offset
          Wvmr_Prof_Rtm(k) = Wvmr_Std_Rtm(k)
          Ozmr_Prof_Rtm(k) = Ozmr_Std_Rtm(k)
    enddo

    !--- Rtm Levels within standard nwp Levels above the surface
    do k = Highest_Level_Rtm_Nwp, Lowest_Level_Rtm_Nwp
          T_Prof_Rtm(k) = T_Nwp(k_Rtm_Nwp(k)) + x_Rtm_Nwp(k) *  &
                     (T_Nwp(k_Rtm_Nwp(k)+1) - T_Nwp(k_Rtm_Nwp(k)))
          Z_Rtm(k) = Z_Nwp(k_Rtm_Nwp(k)) + x_Rtm_Nwp(k) *  &
                     (Z_Nwp(k_Rtm_Nwp(k)+1) - Z_Nwp(k_Rtm_Nwp(k)))
          Wvmr_Prof_Rtm(k) = Wvmr_Nwp(k_Rtm_Nwp(k)) + x_Rtm_Nwp(k) *  &
                     (Wvmr_Nwp(k_Rtm_Nwp(k)+1) - Wvmr_Nwp(k_Rtm_Nwp(k)))
          Ozmr_Prof_Rtm(k) = Ozmr_Nwp(k_Rtm_Nwp(k)) + x_Rtm_Nwp(k) *  &
                     (Ozmr_Nwp(k_Rtm_Nwp(k)+1) - Ozmr_Nwp(k_Rtm_Nwp(k)))
    enddo

   !--- Rtm Levels that are below the lowest nwp Level but above the surface
   do k = Lowest_Level_Rtm_Nwp+1, Sfc_Level_Rtm
           T_Prof_Rtm(k) = Tmpair_Nwp(Lon_Idx,Lat_Idx) + dT_dP_near_Sfc * &
                      (P_Std_Rtm(k) - Psfc_Nwp(Lon_Idx,Lat_Idx))
           Wvmr_Prof_Rtm(k) = Wvmr_Sfc + dWvmr_dP_near_Sfc * &
                      (P_Std_Rtm(k) - Psfc_Nwp(Lon_Idx,Lat_Idx))
           Z_Rtm(k) = dZ_dP_near_Sfc * &
                      (P_Std_Rtm(k) - Psfc_Nwp(Lon_Idx,Lat_Idx))
           Ozmr_Prof_Rtm(k) = Ozmr_Nwp(NLevels_Nwp)
   enddo

   !--- Rtm Levels below the surface
   do k = Sfc_Level_Rtm +1, NLevels_Rtm
         T_Prof_Rtm(k) = Tmpair_Nwp(Lon_Idx,Lat_Idx)
         Z_Rtm(k) = dZ_dP_near_Sfc * &
                   (P_Std_Rtm(k) - Psfc_Nwp(Lon_Idx,Lat_Idx))
         Wvmr_Prof_Rtm(k) = Wvmr_Sfc
         Ozmr_Prof_Rtm(k) = Ozmr_Std_Rtm(k)
   enddo

   !--- if using NCEP reanalysis which has no ozone profile, use default
   if (Nwp_Flag == 2) then
          Ozmr_Prof_Rtm = Ozmr_Std_Rtm
   endif

 end subroutine CONVERT_ATMOS_PROF_NWP_RTM