      module m_interface_crtm
!??? NOT FINISH THE PART REGARDOING THE AIRS things
     
      use crtm_module

      implicit none
!
      private

      public :: setupf_names
      public :: InitGetbt_crtm
!
!!add these here is not structure clear
      integer, save :: nf_zsfc, nf_u10m, nf_v10m
      integer, save :: nf_temp, nf_sphu, nf_ozon

 contains

!  X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X
!
      subroutine setupf_names (nfields2d,nfields3d,nfields,f_names)
      implicit none
      integer :: nfields2d, nfields3d,nfields
      character(len=4) :: f_names(nfields)


      call find_name (nfields2d,f_names(1:nfields2d),'zsfc', nf_zsfc,.true.)
      call find_name (nfields2d,f_names(1:nfields2d),'u10m', nf_u10m,.true.)
      call find_name (nfields2d,f_names(1:nfields2d),'v10m', nf_v10m,.true.)
      call find_name (nfields3d,f_names(1+nfields2d:nfields),'temp', nf_temp,.true.)
      call find_name (nfields3d,f_names(1+nfields2d:nfields),'sphu', nf_sphu,.true.)
      call find_name (nfields3d,f_names(1+nfields2d:nfields),'ozon', nf_ozon,.true.)


      end subroutine setupf_names
    

!  X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X
!
      subroutine InitGetbt_crtm (prof_ndim4,nlevs,nfields2d,nfields3d,  &
                    prof_info,prof_plevs,prof_2d,prof_3d,lcldtop,pcldtop, &
                    linit,d_type,subtype,btchan,nchan,channelinfo,soza)
     

     
!RY 3/12/08
!NOTE: I have problem to pass the channel information if I use two
!      subroutine.  Tentatively fix: use one subroutine

      use m_kinds

  ! ------------
  ! Module usage
  ! ------------
  USE CRTM_SpcCoeff
  USE SpcCoeff_Define

  ! -- Utility modules
  USE Type_Kinds
  USE Message_Handler

  ! -- CRTM module
  USE CRTM_Module
  USE CRTM_Atmosphere_Define

  ! -- Modules to read in Atmosphere and Surface data
  USE CRTM_Atmosphere_Binary_IO
  USE CRTM_Surface_Binary_IO
  USE CRTM_Surface_Define
  USE CRTM_SensorInfo
  USE CRTM_ChannelInfo
  USE CRTM_ChannelInfo_Define
  USE CRTM_surface_ir_emissivity

      implicit none
      integer,intent(in) :: prof_ndim4, nlevs, nfields2d, nfields3d
      integer,intent(in) :: lcldtop
      real(rkind1),intent(in) :: pcldtop
      real(rkind1),intent(in) :: prof_info(prof_ndim4)
      real(rkind1),intent(in) :: prof_2d(nfields2d)
      real(rkind1),intent(in) :: prof_3d(nlevs,nfields3d)
      real(rkind1),intent(in) :: prof_plevs(nlevs)
      integer,intent(in)      :: nchan
      real(rkind3),intent(inout) :: btchan(nchan)
      real(rkind1)            :: soza ! solar zenith angle
      TYPE(CRTM_ChannelInfo_type ),intent(inout)  :: ChannelInfo

      character(len=*),intent (in) :: d_type
      logical, intent (inout)   :: linit
      integer,intent(in) :: subtype     ! ID # for either satellite or instrument type

! local variable
      real(rkind1) :: lza,panglr,rato
      real(rkind1) :: step,start,pi,deg2rad,rad2deg
  CHARACTER( * ), PARAMETER :: PROGRAM_NAME   = 'Test_setup_crtm'
  character(len=*), parameter :: routine_name='interface_crtm'
  ! ---------
  ! Variables
  ! ---------
  INTEGER :: Error_Status, L, i, j, k, cloud_layer_idx, i_prof
  INTEGER :: Allocate_Status

  CHARACTER( 256 ) :: File_Prefix
  CHARACTER( 256 ) :: SpcCoeff_File
  CHARACTER( 256 ) :: TauCoeff_File
  CHARACTER( 256 ) :: AerosolCoeff_File
  CHARACTER( 256 ) :: CloudCoeff_File
  CHARACTER( 256 ) :: EmisCoeff_File 

  TYPE( CRTM_Atmosphere_type ) :: Atmosphere
  TYPE( CRTM_Surface_type ) :: Surface
  TYPE( CRTM_GeometryInfo_type ) :: GeometryInfo
  TYPE( CRTM_RTSolution_type ), ALLOCATABLE, DIMENSION( : ) :: RTSolution
  TYPE( CRTM_Options_type ) :: Options
  INTEGER, PARAMETER :: UNSET = 0
  INTEGER :: CloudType, AerosolType, stype


   REAL( fp_kind ) :: a, alat, alon, height, scover, u, v
   real :: rlat,rlon
! local variable




!==============================================================
!----- ONLY READ COEFFICIENT FILES ONCE FOR EACH INSTRUMENT
!==============================================================


    if (linit) then

! d_type=airs and subtype=420 is the Aqua satellite
!there are five d_type

         if (d_type == 'AIRS_' .and. subtype ==420) then 
           File_Prefix='airs281SUBSET_aqua'
         elseif (d_type == 'AIRS_' .and. subtype==570) then
           File_Prefix='amsua_aqua'
         elseif (d_type == 'HIRS2' ) then
           File_Prefix='hirs2_n14'
         elseif (d_type == 'HIRS3') then
           if (subtype ==207) File_Prefix='hirs3_n16'
           if (subtype ==208) File_Prefix='hirs3_n17'
         elseif (d_type == 'AMSUA') then
!for AMSUA/B the subtype is the satellite ID
!AMSUA starts from satellite=206
             if (subtype ==206) then
               File_Prefix='amsua_n15'
             elseif (subtype ==207) then
               File_Prefix='amsua_n16'
             elseif (subtype ==208) then
               File_Prefix='amsua_n17'
             elseif (subtype ==209) then
               File_Prefix='amsua_n18'
             endif
          elseif (d_type == 'AMSUB') then
!for AMSUA/B the subtype is the satellite ID
!AMSUA starts from satellite=206
             if (subtype ==206) then
               File_Prefix='amsub_n15'
             elseif (subtype ==207) then
               File_Prefix='amsub_n16'
             elseif (subtype ==208) then
               File_Prefix='amsub_n17'
             endif
         else
          print *,' '
          print ('(a,a4,a)'),'Option d_type=',d_type,' not allowed'
          stop
         endif



     File_Prefix = ADJUSTL( File_Prefix )

     ! --------------------
     ! Create the filenames
     ! --------------------

     SpcCoeff_File = TRIM( File_Prefix )//'.SpcCoeff.bin.Big_Endian'
     TauCoeff_File = TRIM( File_Prefix )//'.TauCoeff.bin.Big_Endian'

     AerosolCoeff_File = 'AerosolCoeff.bin.Big_Endian'
     CloudCoeff_File = 'CloudCoeff.bin.Big_Endian'
     EmisCoeff_File    = 'EmisCoeff.bin.Big_Endian'

     !#----------------------------------------------------------------------------#
     !#                          -- INITIALISE THE CRTM --                         #
     !#----------------------------------------------------------------------------#

     WRITE( *, '( /5x, "Initializing the CRTM..." )' )
     write (6,*) 'SpcCoeff_File=', SpcCoeff_File
     write (6,*) 'TauCoeff_File=', TauCoeff_File

     Error_Status = CRTM_Init(ChannelInfo, &
                              SpcCoeff_File     = SpcCoeff_File, &
                              TauCoeff_File     = TauCoeff_File, &
                              AerosolCoeff_File = AerosolCoeff_File, &
                              EmisCoeff_File    = EmisCoeff_File, &
                              CloudCoeff_File = CloudCoeff_File )

     IF ( Error_Status /= SUCCESS ) THEN 
       CALL Display_Message( PROGRAM_NAME, &
                           'Error initializing CRTM', & 
                            Error_Status)  
       STOP
     END IF
   


! print out channel information
    WRITE( *, '( /5x, "Number of channels indexed: ", i5 )' ) ChannelInfo%n_Channels
    WRITE( *, '(  /2x, "Channel         Sensor             NCEP          WMO           WMO     Channel", &
                 &/2x, " Index        Descriptor         Sensor ID   Satellite ID   Sensor ID   Number", &
                 &/2x, "------------------------------------------------------------------------------" )' )
    DO l = 1, ChannelInfo%n_Channels
      WRITE( *, '( 2x, 2x, i4, 2x, ">", a, "<", 5x, i3, 11x, i3, 11x, i3, 7x, i4 )' ) &
                ChannelInfo%Channel_Index( l ), &
                ChannelInfo%Sensor_Descriptor( l ), &
                ChannelInfo%NCEP_Sensor_ID( l ), &
                ChannelInfo%WMO_Satellite_ID( l ), &
                ChannelInfo%WMO_Sensor_ID( l ), &
                ChannelInfo%Sensor_Channel( l )
    END DO

!set linit=false

     linit=.false.
   endif                                !end of if (init)



    pi      = acos(-one)
    deg2rad = pi/180.0
    rad2deg = one/deg2rad
    if (nchan .ne. ChannelInfo%n_Channels) then
    stop 900
    endif
      

     !#----------------------------------------------------------------------------#
     !#    -- ALLOCATE THE OUTPUT RTSolution ARRAY TO THE NUMBER OF CHANNELS --    #
     !#----------------------------------------------------------------------------#

     Error_Status = CRTM_Allocate_Options( ChannelInfo%n_Channels, Options   )


! 

      Atmosphere%n_Clouds = 0
      Atmosphere%n_Absorbers = 2
       Atmosphere%Climatology = 3
      Atmosphere%n_Aerosols = 0
      Atmosphere%n_Layers = lcldtop

      
      Error_Status = CRTM_Allocate_Atmosphere( Atmosphere%n_Layers, &
                                             Atmosphere%n_Absorbers, &
                                             Atmosphere%n_Clouds, &
                                             Atmosphere%n_Aerosols, &
                                             Atmosphere)
     IF ( Error_Status /= SUCCESS ) THEN
        CALL Display_Message(routine_name, &
                          'Error allocating atmosphere', &
                           Error_Status)
        STOP
     END IF

! allocate optical depth array
     ALLOCATE( RTSolution( ChannelInfo%n_Channels ), STAT = Allocate_Status )
     Error_Status = CRTM_Allocate_RTSolution( Atmosphere%n_Layers, RTSolution)
     IF ( Error_Status /= SUCCESS ) THEN
        CALL Display_Message(routine_name, &
                          'Error allocating RT and Options structure array', &
                           Error_Status)
        STOP
     END IF



      Atmosphere%Absorber_ID(1) = 1
      Atmosphere%Absorber_ID(2) = 3
      Atmosphere%Absorber_Units(1) = 3
      Atmosphere%Absorber_Units(2) = 1
      alat = prof_info(2)                     ! units degrees?
      alon = prof_info(3)                     ! units degrees?
! check whether hirs2/3 have this part.
!for both read_bufrtovs and read_airs: need this part:
      if (alon .ge. 360.0) alon = alon - 360.0
      if (alon .lt. 0.00) alon = alon + 360.0

      height = prof_2d(nf_zsfc)               ! when I create height data I did conversion


!--------------------------------------------------------------------------
! those four angles need to be checked and signed
! LIST and unit!!
!questions about the angle variables:
!thingcheck about these two values
! INFO from Mark Liu
! GeometryInfo%Sensor_Zenith_Angle = local zenith angle,  not sensor scan angle
! GeometryInfo%Source_Zenith_Angle = sun zenith angle
!there are two more angle variables: 
! For test, set:
!need the computation of sensor_zenith_angle, unit...

! compute the local zenith angle based on the formula used in gsi
! NOTE: all angles are in degree when call forward model
!-----------------------------------------------------------
!case 1: bufrtovs:
  if ( d_type(1:4) == 'HIRS' ) then
     step   = 1.80
     start  = -49.5
     rato=1.1363987
  else if (d_type == 'AMSUA' ) then
     step   = 3.+ 1.0/3.0
     start = -48. - 1.0/3.0
  else if ( d_type == 'AMSUB' )  then
     step   = 1.1
     start  = -48.95
  elseif (d_type (1:4) == 'AIRS') then
    if (subtype ==420) then      ! this is for airs
      step   = 1.1 
      start = -48.9
  elseif ( subtype == 570) then    ! this is for aqua_amsua
      step   = 3. + 1./3.
      start = -48. - 1./3.
    endif
  endif


! which one if ifov
!---------------------------------
!!! in bufrtovs:    
! (1) hdr1b ='SAID FOVN YEAR MNTH DAYS HOUR MINU SECO CLAT CLON SAZA SOZA BEARAZ SOLAZI HOLS '
! (2)Sensor_zenith_angle=local zenith angle
! (3) check whether the prof_info(10) is solar zenith angle
!---------------------------------

!INFORMATION FROM RE:
!  prof_info(10)=SAZA
!  prof_info(11)=fovs
! geometryinfo CRTM:
! REFER TO setuprad.f90 that read in the data sets created by
!          read_airs.f90 and read_bufrtovs.f90 and then call CRTM:
! (1) local zenith angle, or local zenith angle (satellite zenith angle)
! Var. name here: lza,    data_s(5) in setuprad.f90 
! (2) look angle (scan angle): 
! Var. name here: panglr, data_s(7) in setuprad.f90
! (3) solar zenith angle: 
! Var. name here: SOZApanglr, data_s(7) in setuprad.f90
look angle (scan angle):    var. name here: panglr here  data_s(7) in setuprad.f90
local zenith angle, or local zenith angle: lza
data_s(5)panglr is the data_s (7)the what
!------------------------
!  The unit of geomet. angles to the crtm are in degree !  those values in bufr data are in degree 

   if (d_type(1:4) == 'HIRS' ) then
       panglr=(start+(prof_info(11)-1.0)*step)*deg2rad
       
! separate hirs2 from hirs3 and (hirs4 also, but we do not use hirs4 here)
         if (d_type == 'HIRS2')  then
           lza = asin(rato*sin(panglr))
           lza = lza*rad2deg
         endif
         if (d_type == 'HIRS3')  then
           lza = prof_info(10)
         endif
!  for all hirs 
       if( prof_info(11) <= 28.) lza=-lza

! lza is local zenith angle
   elseif (d_type == 'AMSUA' .or. d_type == 'AMSUB' ) then
       lza=prof_info(10)
       if((d_type == 'AMSUA' .and. prof_info(11) <= 15.) .or.        &
         (d_type == 'AMSUB' .and. prof_info(11) <= 45.) )      &
       lza=-lza
    elseif (d_type == 'AIRS_' ) then
!
!!!!//// WHERE TO PUTshould have a one for airs_amsua****
      panglr
    lza = start + (prof_info(11)-1.)*step
    if ( prof_info(11)
    end if
!  change lza back to the degree



! FOR NOW, WE DO NOT USE SENSOR_ZENITH_ANGLE
!     GeometryInfo%Sensor_Zenith_Angle=prof_info(10) 
!--------------------------------------------------
!  for all five instruments:
!???????  seems that tovs and airs are different?????

       geometryinfo%sensor_zenith_angle  = lza     
       geometryinfo%source_zenith_angle= soza     
! azimuth angle seems not in use
      GeometryInfo%Sensor_Azimuth_Angle = 0.0   ! 
      GeometryInfo%Source_Azimuth_Angle = 0.0   ! 
      

      u=prof_2d(nf_u10m)
      v=prof_2d(nf_v10m)
!
! copy input profile (above cloud only) and rescale
! change p from Pa to hPa      
! change q from kg/kg to g/kg
! change ozone units
! top pressure is fixed as 0.005 mb
      Atmosphere%Level_Pressure(0) = 0.005
      do i=1,lcldtop
        Atmosphere%Pressure(i) = prof_plevs(i)*0.01  
        Atmosphere%Level_Temperature(i) = prof_3d(i,nf_temp) 
        Atmosphere%Temperature(i) = Atmosphere%Level_Temperature(i)
        Atmosphere%Absorber(i,1) = 1000.*max(prof_3d(i,nf_sphu), 1.e-7)
        Atmosphere%Absorber(i,2) = 604229.*max(prof_3d(i,nf_ozon), 1.e-10)
      enddo
!
! set pressure at interfaces (edges)

      do i=1,lcldtop-1
        Atmosphere%Level_Pressure(i) = 0.5*(Atmosphere%Pressure(i+1) &
                                          + Atmosphere%Pressure(i))
      enddo
      Atmosphere%Level_Pressure(lcldtop)=pcldtop*0.01
!
! Consider the surface as ocean
      Surface%Land_Coverage = 0.0
      Surface%Snow_Coverage = 0.0
      Surface%Ice_Coverage = 0.0
      Surface%Water_Coverage = 1.0
      Surface%Water_Type = 1
      Surface%Water_Temperature = &
        Atmosphere%Level_Temperature(Atmosphere%n_Layers)
      Surface%Wind_Speed = sqrt(u*u+v*v) 
      Surface%Wind_Direction = ACOS(u/Surface%Wind_Speed) * 180.0/3.14159
      Surface%Salinity = 33.0
!
! Call CRTM
      Error_Status = CRTM_Forward( Atmosphere, &
                                     Surface, &
                                GeometryInfo, &
                                 ChannelInfo, &
                                 RTSolution)
!                                 Options = Options )

      if ( Error_Status /= SUCCESS ) then 
        call Display_Message(routine_name, &
                           'Error in CRTM Forward Model',Error_Status)  
        stop
      endif
!
! copy output into obs_values(1:n_channels)
! Runhua: check to see whther this is the full set of channels in the obs 
! report.  Otherwise we must make sure we put the info in the correct place.

      do i = 1, ChannelInfo%n_Channels
        l = ChannelInfo%Channel_Index(i)
       btchan(l)= RTSolution(l)%Brightness_Temperature
      enddo

!  destroy the RTSolution and atmosphere

!
      Error_Status = CRTM_Destroy_RTSolution( RTSolution )
      deallocate ( RTSolution,  STAT = Allocate_Status )
      if ( Error_Status /= SUCCESS ) then
        call Display_Message(PROGRAM_NAME, &
                          'Error destroying RTSolution structure', &
                           Error_Status )
        stop
      endif

      Error_Status = CRTM_Destroy_Atmosphere( Atmosphere )
      if ( Error_Status /= SUCCESS ) then
        call Display_Message(PROGRAM_NAME, &
                          'Error destroying Atmosphere structure', &
                           Error_Status )
        stop
      endif

      end subroutine InitGetbt_crtm

!
!  X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X
!
      subroutine find_name (nfields,f_names,f_name,nf,lstop)
!
! Find the index (nf) corresponding to the first occurance of the name
! of the field (f_name) within a chacter string (f_names).
! The name is a single character.
!
! input
      logical          :: lstop    ! stop if name not found in list
      integer          :: nfields
      character(len=*) :: f_names(nfields)
      character(len=*) :: f_name
!
! output
      integer :: nf
!
! local
      integer :: n
!
      nf=0
      do n=1,nfields
        if ( (nf == 0) .and. (f_names(n) == f_name) ) then
          nf=n
        endif
      enddo
!
      if (lstop .and. (nf==0)) then
        print *,' '
        print *,' requested field name not found:'
        print *,'    requested name = ',f_name
        print *,'    possible names = '
        print ('(10(1x,a4))'),f_names(1:nfields)
        stop
      endif
!
      end subroutine find_name


   end module m_interface_crtm
