  !
  ! Description:
  !    Parallel ( openmp ) version of RTTOV_DIRECT, _TL, _AD and _K.   
  !    This file contains the definition of RTTOV_PARALLEL_DIRECT, _TL, _AD and _K. 
  !    RTTOV_PARALLEL_DIRECT is activated when _RTTOV_PARALLEL_DIRECT is defined, etc...
  !    The interfaces of these parallel routines are identical to the original ones,
  !    plus the optional parameter nthreads.
  !
  ! Copyright:
  !
  !    This software was developed within the context of
  !    the EUMETSAT Satellite Application Facility on
  !    Numerical Weather Prediction (NWP SAF), under the
  !    Cooperation Agreement dated 25 November 1998, between
  !    EUMETSAT and the Met Office, UK, by one or more partners
  !    within the NWP SAF. The partners in the NWP SAF are
  !    the Met Office, ECMWF, KNMI and MeteoFrance.
  !
  !    Copyright 2007, EUMETSAT, All Rights Reserved.
  !
  !     *************************************************************
  !
  !
  ! Current Code Owner: SAF NWP
  !
  ! History:
  ! Version   Date     Comment
  ! -------   ----     -------
  !          09/2007    Creation P. Marguinaud & P. Brunel
  !          10/2008    Fix for po sounders P. Marguinaud
  !          11/2009    Add strategy and debug parameters P. Marguinaud
  !


Subroutine &
#ifdef _RTTOV_PARALLEL_DIRECT 
rttov_parallel_direct2  &
#endif
#ifdef _RTTOV_PARALLEL_TL 
rttov_parallel_tl2      &
#endif
#ifdef _RTTOV_PARALLEL_AD 
rttov_parallel_ad2      &
#endif
#ifdef _RTTOV_PARALLEL_K
rttov_parallel_k2       &
#endif
      & ( errorstatus      &
      & , chanprof         &
      & , opts             &
      & , profiles         &
#ifdef _RTTOV_PARALLEL_TL
      & , profiles_tl      &
#endif      
#ifdef _RTTOV_PARALLEL_AD
      & , profiles_ad      &
#endif      
#ifdef _RTTOV_PARALLEL_K
      & , profiles_k       &
#endif      
      & , coefs            &
      & , calcemis         &
      & , emissivity       &
#ifdef _RTTOV_PARALLEL_TL
      & , emissivity_tl    &
      & , cloudemissivity_tl &
#endif      
#ifdef _RTTOV_PARALLEL_AD
      & , emissivity_ad    &
      & , cloudemissivity_ad    &
#endif      
#ifdef _RTTOV_PARALLEL_K
      & , emissivity_k     &
      & , cloudemissivity_k &
#endif      
      & , emissivity_out   &
      & , cloudemissivity  &
#ifdef _RTTOV_PARALLEL_TL
      & , emissivity_out_tl&
#endif      
#ifdef _RTTOV_PARALLEL_AD
      & , emissivity_out_ad&
#endif      
#ifdef _RTTOV_PARALLEL_K
      & , emissivity_out_k &
#endif      
      & , transmission     &
#ifdef _RTTOV_PARALLEL_TL
      & , transmission_tl  &
#endif      
#ifdef _RTTOV_PARALLEL_AD
      & , transmission_ad  &
#endif      
#ifdef _RTTOV_PARALLEL_K
      & , transmission_k   &
#endif      
      & , radiancedata     &
#ifdef _RTTOV_PARALLEL_TL
      & , radiancedata_tl  &
#endif      
#ifdef _RTTOV_PARALLEL_AD
      & , radiancedata_ad  &
#endif      
#ifdef _RTTOV_PARALLEL_K
      & , radiancedata_k   &
#endif      
      & , traj             &
#ifdef _RTTOV_PARALLEL_TL
      & , traj_tl          &
#endif      
#ifdef _RTTOV_PARALLEL_AD
      & , traj_ad          &
#endif      
#ifdef _RTTOV_PARALLEL_K
      & , traj_k           &
#endif  
      &,  pccomp           &
#ifdef _RTTOV_PARALLEL_TL
      &,  pccomp_tl        &
#endif
#ifdef _RTTOV_PARALLEL_AD
      &,  pccomp_ad        &
#endif
#ifdef _RTTOV_PARALLEL_K
      &,  pccomp_k         &
      &,  profiles_k_pc    &
      &,  profiles_k_rec   &
#endif
      &,  channels_rec     &
      &,  nthreads         &    
      &,  strategy         &    
      &,  debug            &    
& )   


  Use rttov_const, Only :  &
       & errorstatus_success,&
       & errorstatus_warning,&
       & errorstatus_fatal,  &  
       & errorstatus_info,   &
       & sensor_id_po
  
! Imported Type Definitions:
  Use rttov_types, Only :            &
        & rttov_options,             &
        & rttov_chanprof,            &
        & rttov_pccomp,              &
        & rttov_coefs,               &
        & profile_Type,              &
        & transmission_type,         &
        & radiance_Type,             &
        & rttov_traj
                      
  
  Use parkind1, Only : jpim,jprb,jplm
  Implicit None



!subroutine arguments:
  Integer(Kind=jpim),            Intent(out)   :: errorstatus ! return flag
  Type(rttov_chanprof),          Intent(in)    :: chanprof(:)
  Type(rttov_options),           Intent(in)    :: opts
  Type(profile_Type),            Intent(in)    :: profiles( : ) ! Atmospheric profiles
  Type(rttov_coefs),             Intent(in)    :: coefs
  Logical(Kind=jplm),            Intent(in)    :: calcemis( : )   ! switch for emmissivity calc.
  Real(Kind=jprb),               Intent(inout) :: emissivity( : ) ! surface emmissivity
  Real(Kind=jprb),               Intent(inout) :: emissivity_out( : ) ! surface emmissivity 
  Real(Kind=jprb),               Intent(inout) :: cloudemissivity( : ) ! cloud emissivity
  Type(transmission_type),       Intent(inout) :: transmission   ! transmittances and layer optical depths
  Type(radiance_Type),           Intent(inout) :: radiancedata   ! radiances (mw/cm-1/ster/sq.m) and degK

#ifdef _RTTOV_PARALLEL_TL
  Type(profile_Type),            Intent(inout) :: profiles_tl( : ) 
  Real(Kind=jprb),               Intent(inout) :: emissivity_tl( : ) 
  Real(Kind=jprb),               Intent(inout) :: emissivity_out_tl(:)
  Real(Kind=jprb),               Intent(inout) :: cloudemissivity_tl( : ) 
  Type(transmission_type),       Intent(inout) :: transmission_tl
  Type(radiance_Type),           Intent(inout) :: radiancedata_tl 
#endif

#ifdef _RTTOV_PARALLEL_AD
  Type(profile_Type),            Intent(inout) :: profiles_ad( : ) 
  Real(Kind=jprb),               Intent(inout) :: emissivity_ad( : ) 
  Real(Kind=jprb),               Intent(inout) :: emissivity_out_ad( : ) 
  Real(Kind=jprb),               Intent(inout) :: cloudemissivity_ad( : )
  Type(transmission_type),       Intent(inout) :: transmission_ad
  Type(radiance_Type),           Intent(inout) :: radiancedata_ad 
#endif

#ifdef _RTTOV_PARALLEL_K
  Type(profile_Type),            Intent(inout) :: profiles_k( : ) 
  Real(Kind=jprb),               Intent(inout) :: emissivity_k( : ) 
  Real(Kind=jprb),               Intent(inout) :: emissivity_out_k( : ) 
  Real(Kind=jprb),               Intent(inout) :: cloudemissivity_k( : ) 
  Type(transmission_type),       Intent(inout) :: transmission_k
  Type(radiance_Type), Optional, Intent(inout) :: radiancedata_k 
#endif
  
  Type(rttov_traj),    Optional, Intent(inout) :: traj
#ifdef _RTTOV_PARALLEL_TL
  Type(rttov_traj),    Optional, Intent(inout) :: traj_tl
#endif
#ifdef _RTTOV_PARALLEL_AD
  Type(rttov_traj),    Optional, Intent(inout) :: traj_ad
#endif
#ifdef _RTTOV_PARALLEL_K
  Type(rttov_traj),    Optional, Intent(inout) :: traj_k
#endif
  Type(rttov_pccomp),  Optional, Intent(inout) :: pccomp
#ifdef _RTTOV_PARALLEL_TL
  Type(rttov_pccomp),  Optional, Intent(inout) :: pccomp_tl
#endif
#ifdef _RTTOV_PARALLEL_AD
  Type(rttov_pccomp),  Optional, Intent(inout) :: pccomp_ad
#endif
#ifdef _RTTOV_PARALLEL_K
  Type(rttov_pccomp),  Optional, Intent(inout) :: pccomp_k
  Type(profile_type),  Optional, Intent(inout) :: profiles_k_pc(:)
  Type(profile_type),  Optional, Intent(inout) :: profiles_k_rec(:)
#endif
  Integer(Kind=jpim),  Optional, Intent(in)    :: channels_rec(:)
  

  Integer(Kind=jpim),  Optional, Intent(in)    :: nthreads
  Integer(Kind=jpim),  Optional, Intent(in)    :: strategy ! 0 = no strategy (RTTOV computes band limits), 1 = profile-wise, 2 = channel-wise
  Logical(Kind=jplm),  Optional, Intent(in)    :: debug

!INTF_END

#ifdef _RTTOV_PARALLEL_DIRECT
#include "rttov_direct2.interface"
#endif
#ifdef _RTTOV_PARALLEL_TL
#include "rttov_tl2.interface"
#endif
#ifdef _RTTOV_PARALLEL_AD
#include "rttov_ad2.interface"
#endif
#ifdef _RTTOV_PARALLEL_K
#include "rttov_k2.interface"
#endif

#include "rttov_alloc_prof.interface"
#include "rttov_alloc_rad.interface"
#include "rttov_errorreport.interface"
#include "rttov_init_prof.interface"
#include "rttov_add_prof.interface"




  Character(len=*), Parameter :: NameOfRoutine = &
#ifdef _RTTOV_PARALLEL_DIRECT
    "rttov_parallel_direct2"
#endif
#ifdef _RTTOV_PARALLEL_TL
    "rttov_parallel_tl2"
#endif
#ifdef _RTTOV_PARALLEL_AD
    "rttov_parallel_ad2"
#endif
#ifdef _RTTOV_PARALLEL_K
    "rttov_parallel_k2"
#endif



#ifdef _RTTOV_PARALLEL_AD
  Integer(Kind=jpim) :: nprofiles12max
  Integer(Kind=jpim) :: k
#endif
  Integer(Kind=jpim) :: iband, nbands, ibandk
  Integer(Kind=jpim) :: ithread, mthreads
  Integer(Kind=jpim) :: nsplits
  Integer(Kind=jpim) :: ichannelk, iprofilek
  Integer(Kind=jpim) :: nlevels
  Integer(Kind=jpim) :: nchannels_po
  Integer(Kind=jpim) :: nchannels
  Integer(Kind=jpim) :: nprofiles
  Integer(Kind=jpim) :: strategy1
  Logical(Kind=jplm) :: debug1
  Integer            :: mstat
  Integer(Kind=jpim),      Allocatable :: errorstatus_x( : )               ! errorstatus/bands
  Integer(Kind=jpim),      Allocatable :: iband1( : ), iband2( : )         ! band limits
  Integer(Kind=jpim),      Allocatable :: ichannel1( : ), ichannel2( : )   ! channel limits
  Integer(Kind=jpim),      Allocatable :: iprofile1( : ), iprofile2( : )   ! profile limits
  Type(rttov_chanprof),    Allocatable :: chanprof_x( :, : )               ! chanprof/bands

  Type(transmission_type), Allocatable :: transmission_x( : )         ! transmission/bands ( pointer assoc )
  Type(radiance_Type),     Allocatable :: radiancedata_x( : )         ! radiancedata/bands ( pointer assoc )

#ifdef _RTTOV_PARALLEL_TL
  Type(transmission_type), Allocatable :: transmission_tl_x( : )  
  Type(radiance_Type),     Allocatable :: radiancedata_tl_x( : )  
#endif

#ifdef _RTTOV_PARALLEL_AD
  Integer(Kind=jpim)                   :: errorstatus_ad_x
  Type(transmission_type), Allocatable :: transmission_ad_x( : )  
  Type(radiance_Type),     Allocatable :: radiancedata_ad_x( : )  
  Type(profile_Type),      Allocatable :: profiles_ad_x( :, : ) 
#endif
  
#ifdef _RTTOV_PARALLEL_K
  Integer(Kind=jpim)                   :: errorstatus_k_x
  Type(transmission_type), Allocatable :: transmission_k_x( : )  
  Type(radiance_Type),     Allocatable :: radiancedata_k_x( : )  
#endif


! Activated if openmp
!  
!$  Integer, External :: omp_get_thread_num
!$  Integer, External :: omp_get_max_threads
!

  nchannels = size(chanprof)
  nprofiles = size(profiles)

  debug1 = .FALSE.
  If( Present( debug ) ) debug1 = debug

  If( opts%addpc ) Then
    Call Rttov_ErrorReport( errorstatus_fatal, "PC not implemented yet in parallel routines !!", NameOfRoutine )
    errorstatus = errorstatus_fatal
    Return
  EndIf

  If( Present( traj ) ) Then
    Call Rttov_ErrorReport( errorstatus_warning, "traj argument will not be used !!", NameOfRoutine )
  EndIf
#ifdef _RTTOV_PARALLEL_TL
  If( Present( traj_tl ) ) Then
    Call Rttov_ErrorReport( errorstatus_warning, "traj_tl argument will not be used !!", NameOfRoutine )
  EndIf
#endif
#ifdef _RTTOV_PARALLEL_AD
  If( Present( traj_ad ) ) Then
    Call Rttov_ErrorReport( errorstatus_warning, "traj_ad argument will not be used !!", NameOfRoutine )
  EndIf
#endif
#ifdef _RTTOV_PARALLEL_K
  If( Present( traj_k ) ) Then
    Call Rttov_ErrorReport( errorstatus_warning, "traj_k argument will not be used !!", NameOfRoutine )
  EndIf
#endif

  nlevels = profiles( 1 ) % nlevels

  strategy1 = 0
  If( Present( strategy ) ) Then
    strategy1 = strategy
  EndIf


  If( Present( nthreads ) ) Then
    mthreads = nthreads
  Else
! Default value for mthreads
    mthreads = 1_jpim
! Activated if openmp
!  
!$  mthreads = omp_get_max_threads()
!
  EndIf

! Check we do not have more bands than channels

  nsplits = 1_jpim
  nbands = mthreads * nsplits


  If( strategy1 .eq. 0 ) Then
    If( coefs%coef%id_sensor .eq. sensor_id_po ) Then
      If( nbands .gt. nprofiles ) nbands = nprofiles
    Else
      If( nbands .gt. nchannels ) nbands = nchannels
    EndIf
  Else If( strategy1 .eq. 1 ) Then
      nbands = nprofiles
  Else If( strategy1 .eq. 2 ) Then
      If( coefs%coef%id_sensor .eq. sensor_id_po ) Then
        errorstatus = errorstatus_fatal
        Call Rttov_ErrorReport( errorstatus_fatal, "Polarized sounders cannot be run channel by channel", NameOfRoutine )
        Return
      EndIf
      nbands = nchannels
  Else
    errorstatus = errorstatus_fatal
    Call Rttov_ErrorReport( errorstatus_fatal, "Unknown strategy", NameOfRoutine )
    Return
  Endif
  
  Allocate( iband1( mthreads ), stat = mstat )
  If( mstat .ne. 0 ) Goto 100
  Allocate( iband2( mthreads ), stat = mstat )
  If( mstat .ne. 0 ) Goto 100
  
  !
  ! we split the input array channels in bands
  !
  ibandk = 1_jpim
  Do ithread = 1, mthreads
    iband1( ithread ) = ibandk
    ibandk = ibandk + nbands / mthreads - 1_jpim
    If( ithread .eq. mthreads ) ibandk = nbands
    iband2( ithread ) = ibandk
    ibandk = ibandk + 1_jpim 
  EndDo


  If( debug1 ) Then
  Write( *, '(A20," = ",I8)' ) "nprofiles", nprofiles
  Write( *, '(A20," = ",I8)' ) "nchannels", nchannels
  Write( *, '(A20," = ",I5)' ) "nbands", nbands
  Endif
  
  Allocate( ichannel1( nbands ), stat = mstat )
  If( mstat .ne. 0 ) Goto 100
  Allocate( ichannel2( nbands ), stat = mstat )
  If( mstat .ne. 0 ) Goto 100
  
  Allocate( iprofile1( nbands ), stat = mstat )
  If( mstat .ne. 0 ) Goto 100
  Allocate( iprofile2( nbands ), stat = mstat )
  If( mstat .ne. 0 ) Goto 100
  
  Allocate( transmission_x( nbands ), stat = mstat )
  If( mstat .ne. 0 ) Goto 100
  Allocate( radiancedata_x( nbands ), stat = mstat )
  If( mstat .ne. 0 ) Goto 100

#ifdef _RTTOV_PARALLEL_TL
  Allocate( transmission_tl_x( nbands ), stat = mstat )
  If( mstat .ne. 0 ) Goto 100
  Allocate( radiancedata_tl_x( nbands ), stat = mstat )
  If( mstat .ne. 0 ) Goto 100
#endif

#ifdef _RTTOV_PARALLEL_AD
  Allocate( transmission_ad_x( nbands ), stat = mstat )
  If( mstat .ne. 0 ) Goto 100
  Allocate( radiancedata_ad_x( nbands ), stat = mstat )
  If( mstat .ne. 0 ) Goto 100
#endif
  
#ifdef _RTTOV_PARALLEL_K
  Allocate( transmission_k_x( nbands ), stat = mstat )
  If( mstat .ne. 0 ) Goto 100
  Allocate( radiancedata_k_x( nbands ), stat = mstat )
  If( mstat .ne. 0 ) Goto 100
#endif
  
  Allocate( chanprof_x( nchannels / nbands + nbands, nbands ), stat = mstat )
  If( mstat .ne. 0 ) Goto 100
  Allocate( errorstatus_x( nbands ), stat = mstat )
  If( mstat .ne. 0 ) Goto 100
  
  
  ichannel1( : ) = 0_jpim
  ichannel2( : ) = 0_jpim
  
  
  iprofile1( : ) = 0_jpim
  iprofile2( : ) = 0_jpim
  
  !
  ! Build the array of lower/upper bounds of bands
  !

  If( ( strategy1 .eq. 0 ) .and. ( coefs%coef%id_sensor .eq. sensor_id_po ) ) Then
    !
    ! The approach is different for PO sensors
    ! because some channels cannot be separated in the
    ! computations; hence we must split the workload
    ! profile-wise
    !
    ichannelk = 1_jpim
    iprofilek = 1_jpim
    Do iband = 1, nbands
      nchannels_po = 0_jpim
      ichannel1( iband ) = ichannelk
      Do While( ( iprofilek .le. nprofiles ) &
          .and. ( ichannelk .lt. nchannels ) )
        Do While( ( iprofilek .le. nprofiles ) &
            .and. ( ichannelk .lt. nchannels ) )
          ichannelk = ichannelk + 1_jpim
          nchannels_po = nchannels_po + 1_jpim
          If( chanprof( ichannelk )%prof .ne. iprofilek ) Then
            iprofilek = iprofilek + 1_jpim
            Exit
          End If
        End Do
        If( nchannels_po .ge. nchannels / nbands ) Exit
      End Do
      ichannel2( iband ) = ichannelk - 1_jpim
    End Do
    ! What remains is given to the last band
    ichannel2( nbands ) = nchannels
  Else If( strategy1 .eq. 1 ) Then
    ichannelk = 1_jpim
    Do iband = 1, nbands
      ichannel1( iband ) = ichannelk
      iprofilek = chanprof(ichannelk)%prof
      Do
        If( ichannelk .gt. nchannels ) Exit
        If( iprofilek .gt. nprofiles ) Exit
        If( iprofilek .ne. chanprof(ichannelk)%prof ) Exit
        ichannelk = ichannelk + 1_jpim
      EndDo
      ichannel2( iband ) = ichannelk-1_jpim
    EndDo
  Else
    !
    ! For IR & MW sensors, we try to make bands
    ! with the same number of channels
    !
    ichannelk = 1_jpim
    Do iband = 1, nbands
      ichannel1( iband ) = ichannelk
      ichannelk = ichannelk + nchannels / nbands - 1_jpim
      If( iband .eq. nbands ) ichannelk = nchannels
      ichannel2( iband ) = ichannelk
      ichannelk = ichannelk + 1_jpim
    End Do
  End If


 
  chanprof_x( :, : )%chan = 0_jpim
  chanprof_x( :, : )%prof = 0_jpim
  
  Do iband = 1, nbands
    iprofile1( iband ) = chanprof( ichannel1( iband ) )%prof
    iprofile2( iband ) = chanprof( ichannel2( iband ) )%prof
    chanprof_x( 1 : ichannel2( iband ) - ichannel1( iband ) + 1, iband )%chan = &
          & chanprof( ichannel1( iband ) : ichannel2( iband ) )%chan
    Do ichannelk = ichannel1( iband ), ichannel2( iband )
      chanprof_x( ichannelk - ichannel1( iband ) + 1, iband )%prof &
        & = chanprof( ichannelk )%prof - iprofile1( iband ) + 1
    EndDo
  EndDo 
  
#ifdef _RTTOV_PARALLEL_AD
  !
  ! for rttov_ad, we have to allocate a temporary array holding
  ! results from every band, and add the contribution of every band
  ! after the parallel loop
  !
  nprofiles12max = maxval( iprofile2( 1 : nbands ) - iprofile1( 1 : nbands ) + 1 )
  Allocate( profiles_ad_x( nprofiles12max, nbands ), stat = mstat )
  If( mstat .ne. 0 ) Goto 100
  Do iband = 1, nbands
    call rttov_alloc_prof( errorstatus_ad_x, nprofiles12max, profiles_ad_x( :, iband ), nlevels, opts, 1_jpim, & 
      init = .true._jplm, coefs=coefs )
    If( errorstatus_ad_x .ne. 0 ) Goto 100
    Call rttov_init_prof(profiles_ad_x( :, iband ))
  EndDo
#endif
  
  
  If( debug1 ) Then

  Do ithread = 1, mthreads
    Write( *, * )
    Write( *, '(A20," = ",I5)' ) "thread", ithread
    Write( *, '(A20," = ",I5)' ) "iband1", iband1( ithread )
    Write( *, '(A20," = ",I5)' ) "iband2", iband2( ithread )
    Write( *, * )
  EndDo

  Do iband = 1, nbands
    Write( *, * )
    Write( *, '(A20," = ",I5)' ) "band", iband
    Write( *, '(A20," = ",100(I5))' ) "channels", chanprof( ichannel1( iband ) : ichannel2( iband ) )%chan
    Write( *, '(A20," = ",100(I5))' ) "lprofiles", chanprof( ichannel1( iband ) : ichannel2( iband ) )%prof
    Write( *, '(A20," = ",100(I5))' ) "chanprof_x%chan", chanprof_x( :, iband )%chan
    Write( *, '(A20," = ",100(I5))' ) "chanprof_x%prof", chanprof_x( :, iband )%prof
    Write( *, * )
  EndDo

  Endif
  
  
  !
  ! we make the correct pointers associations between temporary transmissions
  ! and radiances and arrays passed as input
  !
  Do iband = 1, nbands
  
    Call AssociateTransmission( ichannel1( iband ), ichannel2( iband ), &
     & transmission_x( iband ), transmission )
    Call AssociateRadiance( ichannel1( iband ), ichannel2( iband ), &
     & radiancedata_x( iband ), radiancedata )

#ifdef _RTTOV_PARALLEL_TL
    Call AssociateTransmission( ichannel1( iband ), ichannel2( iband ), &
     & transmission_tl_x( iband ), transmission_tl )
    Call AssociateRadiance( ichannel1( iband ), ichannel2( iband ), &
     & radiancedata_tl_x( iband ), radiancedata_tl )
#endif
    
#ifdef _RTTOV_PARALLEL_AD
    Call AssociateTransmission( ichannel1( iband ), ichannel2( iband ), &
     & transmission_ad_x( iband ), transmission_ad )
    Call AssociateRadiance( ichannel1( iband ), ichannel2( iband ), &
     & radiancedata_ad_x( iband ), radiancedata_ad )
#endif
    
#ifdef _RTTOV_PARALLEL_K
    Call AssociateTransmission( ichannel1( iband ), ichannel2( iband ), &
     & transmission_k_x( iband ), transmission_k )
    If( Present( radiancedata_k ) ) Then
      Call AssociateRadiance( ichannel1( iband ), ichannel2( iband ), &
       & radiancedata_k_x( iband ), radiancedata_k )
    Else
      !
      ! special case for the optional radiancedata_k of rttov_k
      ! we chose to provide the default
      !
      Call rttov_alloc_rad( errorstatus_k_x, Int(ichannel2( iband ) - ichannel1( iband ) + 1,jpim), &
      & radiancedata_k_x( iband ), nlevels-1_jpim, 1_jpim )
      If( errorstatus_k_x .ne. 0 ) Goto 100

      radiancedata_k_x( iband ) % clear( : )       = 0._jprb
      radiancedata_k_x( iband ) % cloudy( : )      = 0._jprb
      radiancedata_k_x( iband ) % bt_clear( : )    = 0._jprb
      radiancedata_k_x( iband ) % upclear( : )     = 0._jprb
      radiancedata_k_x( iband ) % reflclear( : )   = 0._jprb
      radiancedata_k_x( iband ) % overcast( :, : ) = 0._jprb
      radiancedata_k_x( iband ) % bt( : )          = 1._jprb
      radiancedata_k_x( iband ) % total( : )       = 1._jprb
      
    EndIf
#endif
    
    
  End Do

  
  
  
!$OMP PARALLEL DO PRIVATE( iband ) NUM_THREADS( mthreads ) SCHEDULE( DYNAMIC )
  Do iband = 1, nbands

  If( debug1 ) Then
!    Print *, "================================================"
!    Print *, "THREAD = ", OMP_GET_THREAD_NUM(), " BAND = ", iband, &
!  & " ichannel1 = ", ichannel1( iband ), " ichannel2 = ", ichannel2( iband )

    Write( *, * )
#ifdef _RTTOV_PARALLEL_DIRECT
    Write( *, '(A)' ) "========================= rttov_direct2 ========================="
#endif
#ifdef _RTTOV_PARALLEL_TL
    Write( *, '(A)' ) "=========================== rttov_tl2 ==========================="
#endif
#ifdef _RTTOV_PARALLEL_AD
    Write( *, '(A)' ) "=========================== rttov_ad2 ==========================="
#endif
#ifdef _RTTOV_PARALLEL_K
    Write( *, '(A)' ) "=========================== rttov_k2 ============================"
#endif
    Write( *, '(A20)' ) "errorstatus"
    Write( *, '(A20," = ",I5)' ) "nprofiles", iprofile2( iband ) - iprofile1( iband ) + 1
    Write( *, '(A20," = ",I5)' ) "nchannels", ichannel2( iband ) - ichannel1( iband ) + 1
    Write( *, '(A20," = ",100(I5))' ) "channels", chanprof_x( :, iband )%chan
    Write( *, '(A20," = ",100(I5))' ) "lprofiles", chanprof_x( :, iband )%prof
    Write( *, '(A20," = ",L5)' ) "opts", opts
    Write( *, '(A20," = ",A20," ( ",I5," : ",I5," ) ")' ) "profiles",          &
    &"profiles", iprofile1( iband ), iprofile2( iband )
#ifdef _RTTOV_PARALLEL_TL
    Write( *, '(A20," = ",A20," ( ",I5," : ",I5," ) ")' ) "profiles_tl",       &
    &"profiles_tl", iprofile1( iband ), iprofile2( iband )
#endif
#ifdef _RTTOV_PARALLEL_AD
    Write( *, '(A20," = ",A20," ( ",I5," : ",I5," ) ")' ) "profiles_ad",       &
    &"profiles_ad", iprofile1( iband ), iprofile2( iband )
#endif
#ifdef _RTTOV_PARALLEL_K
    Write( *, '(A20," = ",A20," ( ",I5," : ",I5," ) ")' ) "profiles_k",        &
    &"profiles_k", ichannel1( iband ), ichannel2( iband )
#endif
    Write( *, '(A20)' ) "coefs"
    Write( *, '(A20," = ",A20," ( ",I5," : ",I5," ) ")' ) "calcemis",          &
    &"calcemis", ichannel1( iband ), ichannel2( iband )
    Write( *, '(A20," = ",A20," ( ",I5," : ",I5," ) ")' ) "emissivity",        &
    &"emissivity", ichannel1( iband ), ichannel2( iband )
#ifdef _RTTOV_PARALLEL_TL
    Write( *, '(A20," = ",A20," ( ",I5," : ",I5," ) ")' ) "emissivity_tl",     &
    &"emissivity_tl", ichannel1( iband ), ichannel2( iband )
#endif
#ifdef _RTTOV_PARALLEL_AD
    Write( *, '(A20," = ",A20," ( ",I5," : ",I5," ) ")' ) "emissivity_ad",     &
    &"emissivity_ad", ichannel1( iband ), ichannel2( iband )
#endif
#ifdef _RTTOV_PARALLEL_K
    Write( *, '(A20," = ",A20," ( ",I5," : ",I5," ) ")' ) "emissivity_k",      &
    &"emissivity_k", ichannel1( iband ), ichannel2( iband )
#endif
    Write( *, '(A20," = ",A20," ( ",I5," : ",I5," ) ")' ) "emissivity_out",    &
    &"emissivity_out", ichannel1( iband ), ichannel2( iband )
#ifdef _RTTOV_PARALLEL_TL
    Write( *, '(A20," = ",A20," ( ",I5," : ",I5," ) ")' ) "emissivity_out_tl", &
    &"emissivity_out_tl", ichannel1( iband ), ichannel2( iband )
#endif
#ifdef _RTTOV_PARALLEL_AD
    Write( *, '(A20," = ",A20," ( ",I5," : ",I5," ) ")' ) "emissivity_out_ad", &
    &"emissivity_out_ad", ichannel1( iband ), ichannel2( iband )
#endif
#ifdef _RTTOV_PARALLEL_K
    Write( *, '(A20," = ",A20," ( ",I5," : ",I5," ) ")' ) "emissivity_out_k",  &
    &"emissivity_out_k", ichannel1( iband ), ichannel2( iband )
#endif
    Write( *, '(A20)' ) "transmission"
#ifdef _RTTOV_PARALLEL_TL
    Write( *, '(A20)' ) "transmission_tl"
#endif
#ifdef _RTTOV_PARALLEL_AD
    Write( *, '(A20)' ) "transmission_ad"
#endif
#ifdef _RTTOV_PARALLEL_K
    Write( *, '(A20)' ) "transmission_k"
#endif
    Write( *, '(A20)' ) "radiancedata"
#ifdef _RTTOV_PARALLEL_TL
    Write( *, '(A20)' ) "radiancedata_tl"
#endif
#ifdef _RTTOV_PARALLEL_AD
    Write( *, '(A20)' ) "radiancedata_ad"
#endif
#ifdef _RTTOV_PARALLEL_K
    Write( *, '(A20)' ) "radiancedata_k"
#endif
    Write( *, * )
    
  Endif



    Call &
#ifdef _RTTOV_PARALLEL_DIRECT
      & rttov_direct2                                                           &
#endif
#ifdef _RTTOV_PARALLEL_TL
      & rttov_tl2                                                               &
#endif
#ifdef _RTTOV_PARALLEL_AD
      & rttov_ad2                                                               &
#endif
#ifdef _RTTOV_PARALLEL_K
      & rttov_k2                                                                &
#endif
      & ( errorstatus_x( iband )                                               &
      & , chanprof_x( 1 : ichannel2( iband ) - ichannel1( iband ) + 1, iband ) &
      & , opts                                                                 &
      & , profiles( iprofile1( iband ) : iprofile2( iband ) )                  &
#ifdef _RTTOV_PARALLEL_TL
      & , profiles_tl( iprofile1( iband ) : iprofile2( iband ) )               &
#endif
#ifdef _RTTOV_PARALLEL_AD
      & , profiles_ad_x( :, iband )                                            &
#endif
#ifdef _RTTOV_PARALLEL_K
      & , profiles_k( ichannel1( iband ) : ichannel2( iband ) )                &
#endif
      & , coefs                                                                &
      & , calcemis( ichannel1( iband ) : ichannel2( iband ) )                  &
      & , emissivity( ichannel1( iband ) : ichannel2( iband ) )                &
#ifdef _RTTOV_PARALLEL_TL
      & , emissivity_tl( ichannel1( iband ) : ichannel2( iband ) )             &
#endif
#ifdef _RTTOV_PARALLEL_AD
      & , emissivity_ad( ichannel1( iband ) : ichannel2( iband ) )             &
#endif
#ifdef _RTTOV_PARALLEL_K
      & , emissivity_k( ichannel1( iband ) : ichannel2( iband ) )              &
#endif
      & , emissivity_out( ichannel1( iband ) : ichannel2( iband ) )            &
      & , cloudemissivity( ichannel1( iband ) : ichannel2( iband ) )           &
#ifdef _RTTOV_PARALLEL_TL
      & , cloudemissivity_tl( ichannel1( iband ) : ichannel2( iband ) )        &
      & , emissivity_out_tl( ichannel1( iband ) : ichannel2( iband ) )         &
#endif
#ifdef _RTTOV_PARALLEL_AD
      & , emissivity_out_ad( ichannel1( iband ) : ichannel2( iband ) )         &
#endif
#ifdef _RTTOV_PARALLEL_K
      & , emissivity_out_k( ichannel1( iband ) : ichannel2( iband ) )          &
#endif
      & , transmission_x( iband )                                              &
#ifdef _RTTOV_PARALLEL_TL
      & , transmission_tl_x( iband )                                           &
#endif
#ifdef _RTTOV_PARALLEL_AD
      & , transmission_ad_x( iband )                                           &
#endif
#ifdef _RTTOV_PARALLEL_K
      & , transmission_k_x( iband )                                            &
#endif
      & , radiancedata_x( iband )                                              &
#ifdef _RTTOV_PARALLEL_TL
      & , radiancedata_tl_x( iband )                                           &
#endif
#ifdef _RTTOV_PARALLEL_AD
      & , radiancedata_ad_x( iband )                                           &
#endif
#ifdef _RTTOV_PARALLEL_K
      & , radiancedata_k_x( iband )                                            &
#endif
      &,  pccomp                = pccomp                                       &
#ifdef _RTTOV_PARALLEL_TL         
      &,  pccomp_tl             = pccomp_tl                                    &
#endif                                                                         
#ifdef _RTTOV_PARALLEL_AD         
      &,  pccomp_ad             = pccomp_ad                                    &
#endif                                                                         
#ifdef _RTTOV_PARALLEL_K         
      &,  pccomp_k              = pccomp_k                                     &
      &,  profiles_k_pc         = profiles_k_pc                                &
      &,  profiles_k_rec        = profiles_k_rec                               &
#endif                                                                         
      &,  channels_rec          = channels_rec                                 &
      & ) 
   
!     Print *, iband, errorstatus_x(iband )
   
  End Do
!$OMP END PARALLEL DO


#ifdef _RTTOV_PARALLEL_AD
 !
 ! add every contribution in output profile_k array
 !
  Do iband = 1, nbands
    Do iprofilek = iprofile1( iband ), iprofile2( iband )
!      Print *, "profiles_ad( ", iprofilek, " ) = ", &
!      & "profiles_ad( ", iprofilek, " ) + profiles_ad_x( ", &
!      & iprofilek - iprofile1( iband ) + 1, ", ", iband, " )"
      k = iprofilek - iprofile1( iband ) + 1
      Call rttov_add_prof( profiles_ad( iprofilek:iprofilek ), &
                         & profiles_ad( iprofilek:iprofilek ), &
                         & profiles_ad_x( k:k, iband ) )
    EndDo
  EndDo
#endif


#ifdef _RTTOV_PARALLEL_K
  If( .not. Present( radiancedata_k ) ) Then
    Do iband = 1, nbands

      Call rttov_alloc_rad( errorstatus_k_x, Int(ichannel2( iband ) - ichannel1( iband ) + 1,jpim), &
      & radiancedata_k_x( iband ), nlevels, 0_jpim )

    EndDo
  EndIf
#endif


!  Do iband = 1, nbands
!    Call CopyTransmission( ichannel1( iband ), ichannel2( iband ), transmission_x( iband ), transmission )
!  EndDo

  !
  ! we have to construct the global errorstatus array 
  ! from what we got from every thread
  !
  errorstatus = errorstatus_success
  If (Any(errorstatus_x( : ) /= errorstatus_success)) errorstatus = errorstatus_fatal

!  Do iband = 1, nbands
!    Do iprofilek = iprofile1( iband ), iprofile2( iband )
!      Print *, iprofilek, " <= ", iband, iprofilek - iprofile1( iband ) + 1, &
!      & errorstatus_x( iprofilek - iprofile1( iband ) + 1, iband )
!      Call SetMaxErrorStatus( &
!       &  errorstatus( iprofilek ), &
!       &  errorstatus_x( iprofilek - iprofile1( iband ) + 1, iband ) )
!    EndDo
!  EndDo

  
  DeAllocate( chanprof_x )
  DeAllocate( errorstatus_x )

  DeAllocate( transmission_x )
  DeAllocate( radiancedata_x )

#ifdef _RTTOV_PARALLEL_TL
  DeAllocate( transmission_tl_x )
  DeAllocate( radiancedata_tl_x )
#endif
  
#ifdef _RTTOV_PARALLEL_AD
  DeAllocate( transmission_ad_x )
  DeAllocate( radiancedata_ad_x )
#endif

#ifdef _RTTOV_PARALLEL_K
  DeAllocate( transmission_k_x )
  DeAllocate( radiancedata_k_x )
#endif


#ifdef _RTTOV_PARALLEL_AD
  Do iband = 1, nbands
    call rttov_alloc_prof( errorstatus_ad_x, nprofiles12max, profiles_ad_x( :, iband ), nlevels, opts, 0_jpim, coefs = coefs)
  EndDo
  DeAllocate( profiles_ad_x )
#endif
  
  DeAllocate( iprofile2 )
  DeAllocate( iprofile1 )

  DeAllocate( ichannel2 )
  DeAllocate( ichannel1 )

  
  DeAllocate( iband2 )
  DeAllocate( iband1 )

  Return
  
100 Continue
    
    Call Rttov_ErrorReport( errorstatus_fatal, "Memory allocation failed", NameOfRoutine )

    Stop  
  
Contains

  ! ancillary routines

!  Subroutine SetMaxErrorStatus( e1, e2 ) 
!    Integer(Kind=jpim), Intent(inout) :: e1
!    Integer(Kind=jpim), Intent(in) :: e2
!    
!    If( e1 .eq. errorstatus_fatal   ) Return
!    If( e2 .eq. errorstatus_fatal   ) e1 = errorstatus_fatal
!    If( e1 .eq. errorstatus_warning ) Return
!    If( e2 .eq. errorstatus_warning ) e1 = errorstatus_warning
!    If( e1 .eq. errorstatus_info    ) Return
!    If( e2 .eq. errorstatus_info    ) e1 = errorstatus_info
!    If( e1 .eq. errorstatus_success ) Return
!    If( e2 .eq. errorstatus_success ) e1 = errorstatus_success
!    
!  End Subroutine
  
  Subroutine AssociateRadiance( ichannel1, ichannel2, radiancedata_x, radiancedata )
    Integer(Kind=jpim), Intent(in) :: ichannel1, ichannel2
    Type(radiance_Type), Intent(in) :: radiancedata
    Type(radiance_Type), Intent(out) :: radiancedata_x
    
    radiancedata_x % clear &
      => radiancedata % clear( ichannel1 : ichannel2 )

    radiancedata_x % cloudy &
      => radiancedata % cloudy( ichannel1 : ichannel2 )
    
    radiancedata_x % total &
      => radiancedata % total( ichannel1 : ichannel2 )
    
    radiancedata_x % bt &
      => radiancedata % bt( ichannel1 : ichannel2 )
    
    radiancedata_x % bt_clear &
      => radiancedata % bt_clear( ichannel1 : ichannel2 )
    
    radiancedata_x % upclear &
      => radiancedata % upclear( ichannel1 : ichannel2 )
    
    radiancedata_x % dnclear &
      => radiancedata % dnclear( ichannel1 : ichannel2 )
    
    radiancedata_x % reflclear &
      => radiancedata % reflclear( ichannel1 : ichannel2 )
    
    radiancedata_x % overcast &
      => radiancedata % overcast( :, ichannel1 : ichannel2 )
  
    radiancedata_x % up &
      => radiancedata % up( :, ichannel1 : ichannel2 )
  
    radiancedata_x % down &
      => radiancedata % down( :, ichannel1 : ichannel2 )
  
    radiancedata_x % surf &
      => radiancedata % surf( :, ichannel1 : ichannel2 )
  
  End Subroutine

  Subroutine AssociateTransmission( ichannel1, ichannel2, transmission_x, transmission )
    Integer(Kind=jpim), Intent(in) :: ichannel1, ichannel2
    Type(transmission_type), Intent(in) :: transmission
    Type(transmission_type), Intent(out) :: transmission_x
    
    transmission_x % tau_levels & 
      => transmission % tau_levels( :, ichannel1 : ichannel2 )

    transmission_x % tau_total &
      => transmission % tau_total( ichannel1 : ichannel2 )
      
  End Subroutine
  
End Subroutine






