module radiation_two_stream_app

  !
  ! Kind type parameter
  !
  use dc_types, only: DP, &      ! Double precision.
    &                 STRING, &  ! Strings.
    &                 TOKEN      ! Keywords.

  implicit none

  private

  integer, parameter :: IDSchemeShortWave = 1
  integer, parameter :: IDSchemeLongWave  = 2

  public :: RadiationTwoStreamApp

  character(*), parameter:: module_name = 'radiation_two_stream_app'
                              ! モジュールの名称.
                              ! Module name
  character(*), parameter:: version = &
    & '$Name: dcpam5-20100224 $' // &
    & '$Id: radiation_two_stream_app.f90,v 1.1 2010-01-11 01:28:11 yot Exp $'
                              ! モジュールのバージョン
                              ! Module version

  !--------------------------------------------------------------------------------------

contains

  !--------------------------------------------------------------------------------------

  !******************************************************************************
  !    subroutine twostreamapp
  !     Calculate Radiative Heating/Cooling Rate 
  !     Using Two-Stream Approximation
  !
  !     Caution!
  !     In infrared calculation, one must specify albedo in infrared, which is
  !     generally zero.
  !******************************************************************************
  !     pi
  !     sdecl    : solar declination angle (degree)
  !     day      : one solar day (sec)
  !     mst       : the mean solar time on the planet prime meridian (sec)
  !     im       : number of the vertical layers
  !     grav     : gravitational acceleration
  !     ntemp    : neutral temperature (K)
  !     press    : pressure at midpoint of layer (Pa)
  !     cp       : specific heat at constant pressure
  !     temps    : surface temperature (K)
  !     q        : Heating Rate
  !     gdf      : downward flux at ground
  !     dod067   : Dust optical depth at 0.67 micron meter
  !     qerat    : ratio of extinction efficiency factor 
  !              : between those at a wavelength and at 0.67 micron meter
  !     ssa      : Single Scattering Albedo
  !     af       : Asynmmetry Factor
  !     sf       : Solar Flux at concerning Wavelength Width
  !     albedo   : Surface Albedo
  !     wn1      : Start wavenumber (ir calculation)
  !     wn2      : End wavenumber (ir calculation)
  !     dwn      : wavenumber interval used in wavenumber integration 
  !              : (ir calculation)
  !     sw       : 1 = calculation for visible wavelength
  !              : 2 = calculation for infrared wavelength
  !******************************************************************************

  subroutine RadiationTwoStreamApp( &
    & QeRatio, xyz_ssa, xyz_af, &
    & SolarFluxTOA, &
    & xy_SurfAlbedo, &
    & IDScheme, &
    & xy_InAngle, &
    & xyr_OptDepBase, &
    & xyr_RadFlux &
    & )

    ! USE statements
    !

    ! メッセージ出力
    ! Message output
    !
    use dc_message, only: MessageNotify

    !
    ! Physical constants settings
    !
    use constants, only: Grav,  & ! $ g $ [m s-2].
                                  !
                                  ! Gravitational acceleration
      &                  CpDry, & ! $ C_p $ [J kg-1 K-1].
                                  ! 乾燥大気の定圧比熱.
                                  ! Specific heat of air at constant pressure
      &                  PI,    & ! $ \pi $ .
                                  ! Circular constant
      &                  StB      ! $ \sigma_{SB} $ .
                                  ! ステファンボルツマン定数.
                                  ! Stefan-Boltzmann constant

    ! 格子点設定
    ! Grid points settings
    !
    use gridset, only: imax, & ! 経度格子点数.
                               ! Number of grid points in longitude
      &                jmax, & ! 緯度格子点数.
                               ! Number of grid points in latitude
      &                kmax    ! 鉛直層数.
                               ! Number of vertical level

!!$  use pf_module , only : pfint_gq_array


    real(DP), intent(in ) :: QeRatio
    real(DP), intent(in ) :: xyz_ssa       ( 0:imax-1, 1:jmax, 1:kmax )
    real(DP), intent(in ) :: xyz_af        ( 0:imax-1, 1:jmax, 1:kmax )
    real(DP), intent(in ) :: SolarFluxTOA
    real(DP), intent(in ) :: xy_SurfAlbedo ( 0:imax-1, 1:jmax )
    integer , intent(in ) :: IDScheme
    real(DP), intent(in ) :: xy_InAngle    (0:imax-1, 1:jmax)
                              ! sec (入射角).
                              ! sec (angle of incidence)
    real(DP), intent(in ) :: xyr_OptDepBase( 0:imax-1, 1:jmax, 0:kmax )
    real(DP), intent(out) :: xyr_RadFlux   ( 0:imax-1, 1:jmax, 0:kmax )

!!$    real(DP), intent(in ) :: gt         ( 0:imax-1, 1:jmax, 1:kmax )
!!$    real(DP), intent(in ) :: gts        ( 0:imax-1, 1:jmax )
!!$    real(DP), intent(in ) :: gph        ( 0:imax-1, 1:jmax, 0:kmax )

!!$    real(DP), intent(in ) :: emis  ( 0:imax-1, 1:jmax )
!!$    real(DP), intent(in ) :: wn1, wn2
!!$    integer , intent(in ) :: divnum


!!$  real(DP)    , intent(out) :: &
!!$    gor ( im, jm ), goru ( im, jm ), gord ( im, jm ), &
!!$    gsr ( im, jm ), gsru ( im, jm ), gsrd ( im, jm )

    !
    ! opdep       : Optical depth
    !
    real(DP) :: opdep( 0:imax-1, 1:jmax, 0:kmax )

    !
    ! ssa2    : Delta-Function Adjusted Single-Scattering Albedo
    ! af2     : Delta-Function Adjusted Asymmetry Factor
    ! opdep2  : Delta-Function Adjusted Optical Depth
    !
    real(DP) :: xyz_ssa2  ( 0:imax-1, 1:jmax, 1:kmax )
    real(DP) :: xyz_af2   ( 0:imax-1, 1:jmax, 1:kmax )
    real(DP) :: opdep2    ( 0:imax-1, 1:jmax, 0:kmax )
    real(DP) :: exp_opdep2( 0:imax-1, 1:jmax, 0:kmax )

    !
    ! gam?    : Coefficients of Generalized Radiative Transfer Equation
    !
    real(DP) :: gam1( 0:imax-1, 1:jmax, 1:kmax )
    real(DP) :: gam2( 0:imax-1, 1:jmax, 1:kmax )
    real(DP) :: gam3( 0:imax-1, 1:jmax, 1:kmax )
    real(DP) :: gam4( 0:imax-1, 1:jmax, 1:kmax )

    !
    ! upfl      : Upward Flux
    ! dofl      : Downward Flux
    !
    real(DP) :: xyr_UwFlux( 0:imax-1, 1:jmax, 0:kmax )
    real(DP) :: xyr_DwFlux( 0:imax-1, 1:jmax, 0:kmax )

    !
    ! cosz     : cosine of solar zenith angle
    ! cosz2    : cosz squared
    !
    real(DP) :: xy_cosSZA     ( 0:imax-1, 1:jmax )
    real(DP) :: xy_cosSZAInv  ( 0:imax-1, 1:jmax )
    real(DP) :: xy_cosSZAInvsq( 0:imax-1, 1:jmax )

    !
    ! temporary variables
    !
    real(DP) :: dtau( 0:imax-1, 1:jmax, 1:kmax )

    integer  :: i, j, k, l
!!$    integer  :: m

    real(DP) :: lambda( 0:imax-1, 1:jmax, 1:kmax )
    real(DP) :: lsigma( 0:imax-1, 1:jmax, 1:kmax )
    real(DP) :: smalle( 0:imax-1, 1:jmax, 4, 1:kmax )

    !
    ! cupb      : upward C at bottom of layer
    ! cupt      : upward C at top of layer
    ! cdob      : downward C at bottom of layer
    ! cdot      : downward C at top of layer
    !
    real(DP) :: cupb( 0:imax-1, 1:jmax, 1:kmax )
    real(DP) :: cupt( 0:imax-1, 1:jmax, 1:kmax )
    real(DP) :: cdob( 0:imax-1, 1:jmax, 1:kmax )
    real(DP) :: cdot( 0:imax-1, 1:jmax, 1:kmax )

    real(DP) :: aa( imax*jmax, kmax*2 )
    real(DP) :: ba( imax*jmax, kmax*2 )
    real(DP) :: da( imax*jmax, kmax*2 )
    real(DP) :: ea( imax*jmax, kmax*2 )

    real(DP) :: tmpval( 0:imax-1, 1:jmax )

!!$  real(DP) :: mu1
!!$  real(DP) :: gth( im, jm, km+1 ), pfinth( im, jm, km+1 )
!!$  real(DP) :: b0( ijs:ije, 1, km ), b1( ijs:ije, 1, km )
!!$  real(DP) :: gemis


!!$    real(DP) :: inteup( 0:imax-1, 1:jmax, 0:kmax )
!!$    real(DP) :: intedo( 0:imax-1, 1:jmax, 0:kmax )
!!$    real(DP) :: tmpg, tmph, tmpj, tmpk, alp1, alp2, sig1, sig2
!!$    real(DP) :: tmpcos
!!$!      integer(i4b), parameter :: quadn = 3
!!$    integer(i4b), parameter :: quadn = 5
!!$!      integer(i4b), parameter :: quadn = 8
!!$    real(DP) :: qucos ( quadn ), qudcos( quadn )



!!$      call gauleg( 0.0d0, 1.0d0, quadn, qucos, qudcos )

    do k = 1, kmax
      do j = 1, jmax
        do i = 0, imax-1
          if ( xyz_ssa(i,j,k) >= 1.0d0 ) then
            call MessageNotify( 'E', module_name, 'Single Scattering Albedo = %f', &
              & d = (/ xyz_ssa(i,j,k) /) )
          end if
        end do
      end do
    end do


    if( IDScheme .eq. IDSchemeShortWave ) then
      do j = 1, jmax
        do i = 0, imax-1
          if ( xy_InAngle(i,j) > 0.0d0 ) then
            xy_cosSZA     (i,j) = 1.0d0 / xy_InAngle(i,j)
            xy_cosSZAInv  (i,j) = xy_InAngle(i,j)
            xy_cosSZAInvsq(i,j) = xy_cosSZAInv(i,j)**2
          else
            xy_cosSZA     (i,j) = 0.0d0
            xy_cosSZAInv  (i,j) = 0.0d0
            xy_cosSZAInvsq(i,j) = 0.0d0
          end if
        end do
      end do


!!$      gth(:,:,:) = 1.0d100

!!$    else
!!$
!!$      stop 'sw != 1 is not supported.'
!!$
!!$         do ij = ijs, ije
!!$            cosz ( ij, 1 ) = 1.0d100
!!$            cosz2( ij, 1 ) = 1.0d100
!!$         end do
!!$
!!$         do ij = ijs, ije
!!$            gth( ij, 1, 1 ) = gt( ij, 1, 1 )
!!$         end do
!!$         do k = 2, km+1-1
!!$            do ij = ijs, ije
!!$               gth( ij, 1, k ) = ( gt( ij, 1, k-1 ) + gt( ij, 1, k ) ) * 0.5d0
!!$            end do
!!$         end do
!!$         do ij = ijs, ije
!!$            gth( ij, 1, km+1 ) = gts( ij, 1 )
!!$         end do
!!$!         call pfint_array( wn1, wn2, divnum, im, jm, km+1, gth, pfinth, &
!!$!              ijs, ije )
!!$         call pfint_gq_array( wn1, wn2, divnum, im, jm, km+1, gth, pfinth, &
!!$              ijs, ije )

    else
      stop 'Unexpected sw number in twostreamapp.f90'
    end if


    !
    ! calculation of Dust Optical Depth
    !
    opdep(:,:,:) = xyr_OptDepBase * QeRatio


    if( IDScheme .eq. IDSchemeShortWave ) then

      !
      ! Delta-Function Adjustment
      !
      xyz_af2  = xyz_af / ( 1.0d0 + xyz_af )
      xyz_ssa2 = ( 1.0d0 - xyz_af**2 ) * xyz_ssa / ( 1.0d0 - xyz_ssa * xyz_af**2 )

!!$      opdep2(:,:,:) = ( 1.0d0 - xyz_ssa * xyz_af**2 ) * opdep(:,:,:)
      opdep2(:,:,:) = 0.0d0
      do k = kmax-1, 0, -1
        opdep2(:,:,k) =                                         &
          &   opdep2(:,:,k+1)                                   &
          & + ( 1.0d0 - xyz_ssa(:,:,k+1) * xyz_af(:,:,k+1)**2 ) &
          &   * ( opdep(:,:,k) - opdep(:,:,k+1) )
      end do

      do k = 0, kmax
        exp_opdep2(:,:,k) = exp( -opdep2(:,:,k) * xy_cosSZAInv )
      end do

      !
      ! Eddington approximation
      !
      do k = 1, kmax
        gam1(:,:,k) =  ( 7.0d0 - xyz_ssa2(:,:,k) * ( 4.0d0 + 3.0d0 * xyz_af2(:,:,k) ) ) / 4.0d0
        gam2(:,:,k) = -( 1.0d0 - xyz_ssa2(:,:,k) * ( 4.0d0 - 3.0d0 * xyz_af2(:,:,k) ) ) / 4.0d0
        gam3(:,:,k) =  ( 2.0d0 - 3.0d0 * xyz_af2(:,:,k) * xy_cosSZA )        / 4.0d0
        gam4(:,:,k) = 1.0d0 - gam3(:,:,k)
      end do


!!$    else if( sw .eq. 2 ) then
!!$
!!$         !
!!$         ! In infrared, delta-function adjustment is not done. 
!!$         !
!!$         af2  = af
!!$         ssa2 = ssa
!!$
!!$         do k = 1, km+1
!!$            do ij = ijs, ije
!!$               opdep2( ij, 1, k ) = opdep( ij, 1, k )
!!$            end do
!!$         end do
!!$
!!$         do k = 1, km+1
!!$            do ij = ijs, ije
!!$               exp_opdep2( ij, 1, k ) = 1.0d100
!!$            end do
!!$         end do
!!$
!!$         !
!!$         ! Hemispheric mean approximation
!!$         !
!!$         do k = 1, km
!!$            do ij = ijs, ije
!!$               gam1( ij, 1, k ) = 2.0d0 - ssa2 * ( 1.0d0 + af2 )
!!$               gam2( ij, 1, k ) = ssa2 * ( 1.0d0 - af2 )
!!$               gam3( ij, 1, k ) = 1.0d100
!!$               gam4( ij, 1, k ) = 1.0d100
!!$            end do
!!$         end do

      else
         stop 'Unexpected sw number in twostreamapp.f90'
      end if


      do k = 1, kmax
        dtau(:,:,k) = opdep2(:,:,k-1) - opdep2(:,:,k)
      end do
      !
      ! Avoiding singularity when dtau equal to zero 
      !
      do k = 1, kmax
        do j = 1, jmax
          do i = 0, imax-1
            if( ( IDScheme .eq. IDSchemeLongWave ) .and. &
              & ( dtau(i,j,k) .lt. 1.0d-10 ) ) then
              dtau(i,j,k) = 0.0d0
            end if
          end do
        end do
      end do


      lambda(:,:,:) = sqrt( gam1(:,:,:) * gam1(:,:,:) - gam2(:,:,:) * gam2(:,:,:) )
      lsigma(:,:,:) = gam2(:,:,:) / ( gam1(:,:,:) + lambda(:,:,:) )

      do k = 1, kmax
        tmpval(:,:)     = exp( - lambda(:,:,k) * dtau(:,:,k) )
        smalle(:,:,1,k) = 1.0d0 + lsigma(:,:,k) * tmpval(:,:)
        smalle(:,:,2,k) = 1.0d0 - lsigma(:,:,k) * tmpval(:,:)
        smalle(:,:,3,k) = lsigma(:,:,k) + tmpval(:,:)
        smalle(:,:,4,k) = lsigma(:,:,k) - tmpval(:,:)
      end do


      if( IDScheme .eq. IDSchemeShortWave ) then
        do k = 1, kmax
          cupb(:,:,k) =                                                   &
            &   xyz_ssa2(:,:,k) * SolarFluxTOA * exp_opdep2(:,:,k-1)      &
            &   * ( ( gam1(:,:,k) - xy_cosSZAInv ) * gam3(:,:,k)          &
            &       + gam2(:,:,k) * gam4(:,:,k) )                         &
            &   / ( lambda(:,:,k) * lambda(:,:,k) - xy_cosSZAInvsq )
          cupt(:,:,k) =                                                   &
            &   xyz_ssa2(:,:,k) * SolarFluxTOA * exp_opdep2(:,:,k  )      &
            &   * ( ( gam1(:,:,k) - xy_cosSZAInv ) * gam3(:,:,k)          &
            &       + gam2(:,:,k) * gam4(:,:,k) )                         &
            &   / ( lambda(:,:,k) * lambda(:,:,k) - xy_cosSZAInvsq )
          cdob(:,:,k) =                                                   &
            &   xyz_ssa2(:,:,k) * SolarFluxTOA * exp_opdep2(:,:,k-1)      &
            &   * ( ( gam1(:,:,k) + xy_cosSZAInv ) * gam4(:,:,k)          &
            &       + gam2(:,:,k) * gam3(:,:,k) )                         &
            &   / ( lambda(:,:,k) * lambda(:,:,k) - xy_cosSZAInvsq )
          cdot(:,:,k) =                                                   &
            &   xyz_ssa2(:,:,k) * SolarFluxTOA * exp_opdep2(:,:,k  )      &
            &   * ( ( gam1(:,:,k) + xy_cosSZAInv ) * gam4(:,:,k)          &
            &       + gam2(:,:,k) * gam3(:,:,k) )                         &
            &   / ( lambda(:,:,k) * lambda(:,:,k) - xy_cosSZAInvsq )
        end do

!!$      else if( sw .eq. 2 ) then
!!$
!!$         do k = 1, km
!!$            do ij = ijs, ije
!!$               !
!!$               ! Notice!
!!$               ! Avoiding singularity when dtau equal to zero. 
!!$               ! dtau occationally has much smaller value. 
!!$               ! When this occurs, b1 cannot be calculated correctly. 
!!$               !
!!$               if( dtau( ij, 1, k ) .ne. 0.0d0 ) then
!!$                  b0( ij, 1, k ) = pfinth( ij, 1, k )
!!$                  b1( ij, 1, k ) = ( pfinth( ij, 1, k+1 ) - pfinth( ij, 1, k ) ) / dtau( ij, 1, k )
!!$               else
!!$                  b0( ij, 1, k ) = 0.0d0
!!$                  b1( ij, 1, k ) = 0.0d0
!!$               end if
!!$
!!$            end do
!!$         end do
!!$
!!$         do k = 1, km
!!$            do ij = ijs, ije
!!$               mu1 = ( 1.0d0 - ssa2 ) / ( gam1( ij, 1, k ) - gam2( ij, 1, k ) )
!!$
!!$               cupb( ij, 1, k ) = pix2 * mu1 * ( b0( ij, 1, k ) &
!!$                    + b1( ij, 1, k ) &
!!$                    * ( dtau( ij, 1, k ) + 1.0d0 / ( gam1( ij, 1, k ) + gam2( ij, 1, k ) ) ) )
!!$               cupt( ij, 1, k ) = pix2 * mu1 * ( b0( ij, 1, k ) &
!!$                    + b1( ij, 1, k ) &
!!$                    * ( 0.0d0     + 1.0d0 / ( gam1( ij, 1, k ) + gam2( ij, 1, k ) ) ) )
!!$               cdob( ij, 1, k ) = pix2 * mu1 * ( b0( ij, 1, k ) &
!!$                    + b1( ij, 1, k ) &
!!$                    * ( dtau( ij, 1, k ) - 1.0d0 / ( gam1( ij, 1, k ) + gam2( ij, 1, k ) ) ) )
!!$               cdot( ij, 1, k ) = pix2 * mu1 * ( b0( ij, 1, k ) &
!!$                    + b1( ij, 1, k ) &
!!$                    * ( 0.0d0     - 1.0d0 / ( gam1( ij, 1, k ) + gam2( ij, 1, k ) ) ) )
!!$            end do
!!$         end do
!!$
      else
        stop 'Unexpected sw number in twostreamapp.f90'
      end if


      k = 1
      l = 1
      do j = 1, jmax
        do i = 0, imax-1
          aa( i+imax*(j-1)+1, l ) =                                     &
            & 0.0d0
          ba( i+imax*(j-1)+1, l ) =                                     &
            & smalle(i,j,2,k) - xy_SurfAlbedo(i,j) * smalle(i,j,4,k)
          da( i+imax*(j-1)+1, l ) =                                     &
            & smalle(i,j,1,k) - xy_SurfAlbedo(i,j) * smalle(i,j,3,k)
        end do
      end do
      if( IDScheme .eq. IDSchemeShortWave ) then
        do j = 1, jmax
          do i = 0, imax-1
            ea( i+imax*(j-1)+1, l ) =                                              &
              & - cupb(i,j,k)                                                      &
              & + xy_SurfAlbedo(i,j) * cdob(i,j,k)                                 &
              & + xy_SurfAlbedo(i,j) * SolarFluxTOA * exp_opdep2(i,j,0) * xy_cosSZA(i,j)
          end do
        end do
!!$      else if( sw .eq. 2 ) then
!!$         do ij = ijs, ije
!!$!            gemis = 1.0d0
!!$            gemis = emis( ij, 1 )
!!$            ea( (ij-ijs+1), l ) &
!!$                 = -cupb( ij, 1, km ) + ( 1.0d0 - gemis ) * cdob( ij, 1, km ) &
!!$                 + gemis * pi * pfinth( ij, 1, km+1 )
!!$         end do
      else
        stop 'Unexpected sw number in twostreamapp.f90'
      end if


!!$      do k = 1+1, kmax
!!$
!!$        do j = 1, jmax
!!$          do i = 0, imax-1
!!$
!!$            ! for even l
!!$            !
!!$            l = 2 * (k-1)
!!$            !
!!$            aa( i+imax*(j-1), l ) =                    &
!!$              &   smalle(i,j,3,k) * smalle(i,j,4,k-1)  &
!!$              & - smalle(i,j,1,k) * smalle(i,j,2,k-1)
!!$            ba( i+imax*(j-1), l ) =                    &
!!$              &   smalle(i,j,1,k) * smalle(i,j,1,k-1)  &
!!$              & - smalle(i,j,3,k) * smalle(i,j,3,k-1)
!!$            da( i+imax*(j-1), l ) =                    &
!!$              &   smalle(i,j,2,k) * smalle(i,j,3,k  )  &
!!$              & - smalle(i,j,1,k) * smalle(i,j,4,k  )
!!$            ea( i+imax*(j-1), l ) =                    &
!!$              &   smalle(i,j,3,k) * ( cupt(i,j,k-1) - cupb(i,j,k) ) &
!!$              & - smalle(i,j,1,k) * ( cdot(i,j,k-1) - cdob(i,j,k) )
!!$
!!$            ! for odd l
!!$            !
!!$            l = 2 * (k-1) + 1
!!$
!!$            aa( i+imax*(j-1), l ) = &
!!$              &   smalle(i,j,1,k-1) * smalle(i,j,4,k-1) &
!!$              & - smalle(i,j,2,k-1) * smalle(i,j,3,k-1)
!!$            ba( i+imax*(j-1), l ) = &
!!$              &   smalle(i,j,2,k  ) * smalle(i,j,2,k-1) &
!!$              & - smalle(i,j,4,k  ) * smalle(i,j,4,k-1)
!!$            da( i+imax*(j-1), l ) = &
!!$              &   smalle(i,j,1,k  ) * smalle(i,j,2,k-1) &
!!$              & - smalle(i,j,3,k  ) * smalle(i,j,4,k-1)
!!$            ea( i+imax*(j-1), l ) = &
!!$              &   smalle(i,j,2,k-1) * ( cupt(i,j,k-1) - cupb(i,j,k) ) &
!!$              & - smalle(i,j,4,k-1) * ( cdot(i,j,k-1) - cdob(i,j,k) )
!!$          end do
!!$        end do
!!$      end do


      do k = 1, kmax-1

        do j = 1, jmax
          do i = 0, imax-1

            l = 2 * k     ! equation number
            !
            aa( i+imax*(j-1)+1, l ) =                                     &
              &   smalle(i,j,3,k+1) * smalle(i,j,4,k  )                 &
              & - smalle(i,j,1,k+1) * smalle(i,j,2,k  )
            ba( i+imax*(j-1)+1, l ) =                                     &
              &   smalle(i,j,1,k+1) * smalle(i,j,1,k  )                 &
              & - smalle(i,j,3,k+1) * smalle(i,j,3,k  )
            da( i+imax*(j-1)+1, l ) =                                     &
              &   smalle(i,j,2,k+1) * smalle(i,j,3,k+1)                 &
              & - smalle(i,j,1,k+1) * smalle(i,j,4,k+1)
            ea( i+imax*(j-1)+1, l ) =                                     &
              &   smalle(i,j,3,k+1) * ( cupt(i,j,k  ) - cupb(i,j,k+1) ) &
              & - smalle(i,j,1,k+1) * ( cdot(i,j,k  ) - cdob(i,j,k+1) )

            l = 2 * k + 1 ! equation number
            !
            aa( i+imax*(j-1)+1, l ) =                                     &
              &   smalle(i,j,1,k  ) * smalle(i,j,4,k  )                 &
              & - smalle(i,j,2,k  ) * smalle(i,j,3,k  )
            ba( i+imax*(j-1)+1, l ) =                                     &
              &   smalle(i,j,2,k+1) * smalle(i,j,2,k  )                 &
              & - smalle(i,j,4,k+1) * smalle(i,j,4,k  )
            da( i+imax*(j-1)+1, l ) =                                     &
              &   smalle(i,j,1,k+1) * smalle(i,j,2,k  )                 &
              & - smalle(i,j,3,k+1) * smalle(i,j,4,k  )
            ea( i+imax*(j-1)+1, l ) =                                     &
              &   smalle(i,j,2,k  ) * ( cupt(i,j,k  ) - cupb(i,j,k+1) ) &
              & - smalle(i,j,4,k  ) * ( cdot(i,j,k  ) - cdob(i,j,k+1) )
          end do
        end do
      end do


      k = kmax
      l = 2 * kmax
      do j = 1, jmax
        do i = 0, imax-1
          aa( i+imax*(j-1)+1, l ) = -smalle(i,j,2,k)
          ba( i+imax*(j-1)+1, l ) =  smalle(i,j,1,k)
          da( i+imax*(j-1)+1, l ) = 0.0d0
          ea( i+imax*(j-1)+1, l ) = -cdot(i,j,k) + 0.0d0
        end do
      end do

!!$      call tridiag( imax*jmax, 2*kmax, aa, ba, da, ea )
      call tridiaginv( imax*jmax, 2*kmax, aa, ba, da, ea )

      if( IDScheme .eq. IDSchemeShortWave ) then

        do k = 1, kmax
          do j = 1, jmax
            do i = 0, imax-1
              xyr_UwFlux(i,j,k-1) = &
                &   ea( i+imax*(j-1)+1, 2*k  ) * smalle(i,j,1,k) &
                & + ea( i+imax*(j-1)+1, 2*k-1) * smalle(i,j,2,k) &
                & + cupb(i,j,k)
              xyr_DwFlux(i,j,k-1) = &
                &   ea( i+imax*(j-1)+1, 2*k  ) * smalle(i,j,3,k) &
                & + ea( i+imax*(j-1)+1, 2*k-1) * smalle(i,j,4,k) &
                & + cdob(i,j,k)
            end do
          end do
        end do

        k = kmax
        do j = 1, jmax
          do i = 0, imax-1
            xyr_UwFlux(i,j,k) =                              &
              &   ea( i+imax*(j-1)+1, 2*k  ) * smalle(i,j,3,k) &
              & - ea( i+imax*(j-1)+1, 2*k-1) * smalle(i,j,4,k) &
              & + cupt(i,j,k)
            xyr_DwFlux(i,j,k) =                              &
              &   ea( i+imax*(j-1)+1, 2*k  ) * smalle(i,j,1,k) &
              & - ea( i+imax*(j-1)+1, 2*k-1) * smalle(i,j,2,k) &
              & + cdot(i,j,k)
          end do
        end do

        !
        ! Addition of Direct Solar Insident
        !
        do k = 0, kmax
          xyr_DwFlux(:,:,k) = xyr_DwFlux(:,:,k) &
            & + SolarFluxTOA * exp_opdep2(:,:,k) * xy_cosSZA
        end do

!!$      else if( sw .eq. 2 ) then
!!$
!!$         ! Source function technique described by Toon et al. [1989] 
!!$         ! is used to calculated infrared heating rate. 
!!$         !
!!$         do k = 1, km+1
!!$            do ij = ijs, ije
!!$               upfl( ij, 1, k ) = 0.0d0
!!$               dofl( ij, 1, k ) = 0.0d0
!!$            end do
!!$         end do
!!$
!!$         do l = 1, quadn
!!$            tmpcos = qucos( l )
!!$
!!$            k = 1
!!$            do ij = ijs, ije
!!$               intedo( ij, 1, k ) = 0.0d0
!!$            end do
!!$            do k = 1, km
!!$               do ij = ijs, ije
!!$                  tmpj = ( ea( (ij-ijs+1), k+k-1 ) + ea( (ij-ijs+1), k+k ) ) * lsigma( ij, 1, k ) &
!!$                       * ( lambda( ij, 1, k ) + 2.0d0 )
!!$                  tmpk = ( ea( (ij-ijs+1), k+k-1 ) - ea( (ij-ijs+1), k+k ) ) &
!!$                       * ( 2.0d0 - lambda( ij, 1, k ) )
!!$                  sig1 = pix2 * ( b0( ij, 1, k ) &
!!$                       - b1( ij, 1, k ) &
!!$                       * ( 1.0d0 / ( gam1( ij, 1, k ) + gam2( ij, 1, k ) ) - 0.5d0 ) )
!!$                  sig2 = pix2 * b1( ij, 1, k )
!!$                  intedo( ij, 1, k+1 ) = intedo( ij, 1, k ) &
!!$                       * exp( -dtau( ij, 1, k ) / tmpcos ) &
!!$                       + tmpj / ( lambda( ij, 1, k ) * tmpcos + 1.0d0 ) &
!!$                       * ( 1.0d0 &
!!$                       - exp( -dtau( ij, 1, k ) &
!!$                       * ( lambda( ij, 1, k ) + 1.0d0 / tmpcos ) ) ) &
!!$                       + tmpk / ( lambda( ij, 1, k ) * tmpcos - 1.0d0 ) &
!!$                       * ( exp( -dtau( ij, 1, k ) / tmpcos ) &
!!$                       -exp( -dtau( ij, 1, k ) * lambda( ij, 1, k ) ) ) &
!!$                       + sig1 * ( 1.0d0 - exp( -dtau( ij, 1, k ) / tmpcos ) ) &
!!$                       + sig2 * ( tmpcos * exp( -dtau( ij, 1, k ) / tmpcos ) &
!!$                       + dtau( ij, 1, k ) - tmpcos )
!!$               end do
!!$            end do
!!$
!!$            k = km+1
!!$            do ij = ijs, ije
!!$!               gemis = 1.0d0
!!$               gemis = emis( ij, 1 )
!!$               inteup( ij, 1, k ) = ( 1.0d0 - gemis ) * intedo( ij, 1, km+1 ) &
!!$                    + gemis * pix2 * pfinth( ij, 1, km+1 )
!!$            end do
!!$            do k = km, 1, -1
!!$               do ij = ijs, ije
!!$                  tmpg = ( ea( (ij-ijs+1), k+k-1 ) + ea( (ij-ijs+1), k+k ) ) &
!!$                       * ( 2.0d0 - lambda( ij, 1, k ) )
!!$                  tmph = ( ea( (ij-ijs+1), k+k-1 ) - ea( (ij-ijs+1), k+k ) ) * lsigma( ij, 1, k ) &
!!$                       * ( lambda( ij, 1, k ) + 2.0d0 )
!!$                  alp1 = pix2 * ( b0( ij, 1, k ) &
!!$                       + b1( ij, 1, k ) &
!!$                       * ( 1.0d0 / ( gam1( ij, 1, k ) + gam2( ij, 1, k ) ) - 0.5d0 ) )
!!$                  alp2 = pix2 * b1( ij, 1, k )
!!$                  inteup( ij, 1, k ) = inteup( ij, 1, k+1 ) &
!!$                       * exp( -dtau( ij, 1, k ) / tmpcos ) &
!!$                       + tmpg / ( lambda( ij, 1, k ) * tmpcos - 1.0d0 ) &
!!$                       * ( exp( -dtau( ij, 1, k ) / tmpcos ) &
!!$                       - exp( -dtau( ij, 1, k ) * lambda( ij, 1, k ) ) ) &
!!$                       + tmph / ( lambda( ij, 1, k ) * tmpcos + 1.0d0 ) &
!!$                       * ( 1.0d0 &
!!$                       -exp( -dtau( ij, 1, k ) &
!!$                       * ( lambda( ij, 1, k ) + 1.0d0 / tmpcos ) ) ) &
!!$                       + alp1 * ( 1.0d0 - exp( -dtau( ij, 1, k ) / tmpcos ) ) &
!!$                       + alp2 * ( tmpcos &
!!$                       - ( dtau( ij, 1, k ) + tmpcos ) &
!!$                       * exp( -dtau( ij, 1, k ) / tmpcos ) )
!!$               end do
!!$            end do
!!$
!!$            do k = 1, km+1
!!$               do ij = ijs, ije
!!$                  upfl( ij, 1, k ) = upfl( ij, 1, k ) &
!!$                       + inteup( ij, 1, k ) * qucos( l ) * qudcos( l )
!!$                  dofl( ij, 1, k ) = dofl( ij, 1, k ) &
!!$                       + intedo( ij, 1, k ) * qucos( l ) * qudcos( l )
!!$               end do
!!$            end do
!!$         end do
!!$
      else
        stop 'Unexpected sw number in twostreamapp.f90'
      end if


!!$      do ij = ijs, ije
!!$         gdf( ij, 1 ) = dofl( ij, 1, km+1 )
!!$      end do
!!$      if( sw == 1 ) then
!!$        ! visible
!!$        do ij = ijs, ije
!!$          gdf( ij, 1 ) = upfl( ij, 1, km+1 ) - dofl( ij, 1, km+1 )
!!$        end do
!!$      else
!!$        ! ir
!!$        do ij = ijs, ije
!!$          gdf( ij, 1 ) = dofl( ij, 1, km+1 )
!!$        end do
!!$      end if
!!$      do ij = ijs, ije
!!$        gdf( ij, 1 ) = upfl( ij, 1, km+1 ) - dofl( ij, 1, km+1 )
!!$      end do


!!$      do k = 1, km
!!$         do ij = ijs, ije
!!$            q( ij, 1, k ) &
!!$                 = ( ( upfl( ij, 1, k+1 ) - dofl( ij, 1, k+1 ) ) &
!!$                   - ( upfl( ij, 1, k   ) - dofl( ij, 1, k   ) ) ) &
!!$                   / ( gph ( ij, 1, k+1 ) - gph ( ij, 1, k   ) ) &
!!$                   * grav / cp
!!$         end do
!!$      end do


!!$      if( sw == 1 ) then
!!$        write( 6, * ) 'vis flux=', &
!!$!          -gdf((ijs+ije)/2,1), &
!!$!          -gdf((ijs+ije)/2,1) * ( 1.0d0 - albedo((ijs+ije)/2,1) ), &
!!$          -gdf((ijs+ije)/2,1) * albedo((ijs+ije)/2,1), &
!!$!          upfl((ijs+ije)/2,1,km+1)-dofl((ijs+ije)/2,1,km+1),
!!$          upfl((ijs+ije)/2,1,km+1), &
!!$          upfl((ijs+ije)/2,1,km+1) - gdf((ijs+ije)/2,1) * albedo((ijs+ije)/2,1)
!!$      else
!!$        write( 6, * ) 'ir  flux=', &
!!$          -gdf((ijs+ije)/2,1), &
!!$          upfl((ijs+ije)/2,1,km+1)-dofl((ijs+ije)/2,1,km+1), &
!!$          ' sig * Ts^4 = ', sboltz * gts((ijs+ije)/2,1)**4
!!$      end if

      xyr_RadFlux = xyr_UwFlux - xyr_DwFlux

      if( IDScheme .eq. IDSchemeShortWave ) then
        do j = 1, jmax
          do i = 0, imax-1
            if( xy_cosSZA(i,j) <= 0.0d0 ) then
              do k = 0, kmax
                xyr_RadFlux(i,j,k) = 0.0d0
              end do
            end if
          end do
        end do
      end if


!!$      !
!!$      ! output variables
!!$      !
!!$      do ij = ijs, ije
!!$         goru( ij, 1 ) = upfl( ij, 1, 1 )
!!$         gord( ij, 1 ) = dofl( ij, 1, 1 )
!!$         gsru( ij, 1 ) = upfl( ij, 1, km+1 )
!!$         gsrd( ij, 1 ) = dofl( ij, 1, km+1 )
!!$      end do
!!$      if( sw .eq. 1 ) then
!!$         do ij = ijs, ije
!!$            if( cosz( ij, 1 ) .le. 0.0d0 ) then
!!$               goru( ij, 1 ) = 0.0d0
!!$               gord( ij, 1 ) = 0.0d0
!!$               gsru( ij, 1 ) = 0.0d0
!!$               gsrd( ij, 1 ) = 0.0d0
!!$            end if
!!$         end do
!!$      end if
!!$      do ij = ijs, ije
!!$         gor ( ij, 1 ) = goru( ij, 1 ) - gord( ij, 1 )
!!$         gsr ( ij, 1 ) = gsru( ij, 1 ) - gsrd( ij, 1 )
!!$      end do


    end subroutine RadiationTwoStreamApp


    !------------------------------------------------------------------------------------

    !******************************************************************************
    !      subroutine tridiag
    !      tidiagonal solver
    !******************************************************************************
    !     a(j), b(j), and c(j) are, respectively, the subdiagonal, diagonal, 
    !     and superdiagonal entries in row j. 
    !     a(1) and c(jmx) need not be initialized. 
    !     The output is in f; a, b, and c are unchanged. 
    !******************************************************************************
    !     jmx    : dimensions of all the following arrays
    !     a      : sub (lower) diagonal
    !     b      : center diagonal
    !     c      : super (upper) diagonal
    !     f      : right hand side
    !******************************************************************************

    subroutine tridiaginv( mm, jmx, a, b, c, f )

      integer , intent(in   ) :: mm, jmx
      real(DP), intent(in   ) :: a( mm, jmx ),b( mm, jmx )
      real(DP), intent(in   ) :: c( mm, jmx )
      real(DP), intent(inout) :: f( mm, jmx )


      ! Local variables
      !
      real(DP) :: a2( mm, jmx )
      real(DP) :: b2( mm, jmx )
      real(DP) :: c2( mm, jmx )
      real(DP) :: f2( mm, jmx )
      integer  :: j

      do j = 1, jmx
        a2(:,j) = c(:,jmx-j+1)
        b2(:,j) = b(:,jmx-j+1)
        c2(:,j) = a(:,jmx-j+1)
        f2(:,j) = f(:,jmx-j+1)
      end do

      call tridiag( mm, jmx, a2, b2, c2, f2 )

      do j = 1, jmx
        f(:,j) = f2(:,jmx-j+1)
      end do

    end subroutine tridiaginv

    !******************************************************************************

    subroutine tridiag( mm, jmx, a, b, c, f )

      integer , intent(in   ) :: mm, jmx
      real(DP), intent(in   ) :: a( mm, jmx ),b( mm, jmx )
      real(DP), intent(in   ) :: c( mm, jmx )
      real(DP), intent(inout) :: f( mm, jmx )


      ! Local variables
      !
      real(DP) :: q( mm, jmx ), p
      integer  :: j, m


      ! Forward elimination sweep
      !
      do m = 1, mm
        q( m, 1 ) = - c( m, 1 ) / b( m, 1 )
        f( m, 1 ) =   f( m, 1 ) / b( m, 1 )
      end do

      do j = 2, jmx
        do m = 1, mm
          p         = 1.0d0 / ( b( m, j ) + a( m, j ) * q( m, j-1 ) )
          q( m, j ) = - c( m, j ) * p
          f( m, j ) = ( f( m, j ) - a( m, j ) * f( m, j-1 ) ) * p
        end do
      end do

      ! Backward pass
      !
      do j = jmx - 1, 1, -1
        do m = 1, mm
          f( m, j ) = f( m, j ) + q( m, j ) * f( m, j+1 )
        end do
      end do

    end subroutine tridiag

    !******************************************************************************
    !      subroutine tridiag
    !      tidiagonal solver
    !******************************************************************************
    !     a(j), b(j), and c(j) are, respectively, the subdiagonal, diagonal, 
    !     and superdiagonal entries in row j. 
    !     a(1) and c(jmx) need not be initialized. 
    !     The output is in f; a, b, and c are unchanged. 
    !******************************************************************************
    !     jmx    : dimensions of all the following arrays
    !     a      : sub (lower) diagonal
    !     b      : center diagonal
    !     c      : super (upper) diagonal
    !     f      : right hand side
    !******************************************************************************

    subroutine tridiag1( jmx, a, b, c, f )

      integer , intent(in   ) :: jmx
      real(DP), intent(in   ) :: a(jmx),b(jmx)
      real(DP), intent(in   ) :: c(jmx)
      real(DP), intent(inout) :: f(jmx)


      ! Local variables
      !
      real(DP) :: q(jmx), p
      integer  :: j


      ! Forward elimination sweep
      !
      q( 1 ) = - c( 1 ) / b( 1 )
      f( 1 ) =   f( 1 ) / b( 1 )

      do j = 2, jmx
        p      = 1.0d0 / ( b( j ) + a( j ) * q( j-1 ) )
        q( j ) = - c( j ) * p
        f( j ) = ( f( j ) - a( j ) * f( j-1 ) ) * p
      end do

      ! Backward pass
      !
      do j = jmx - 1, 1, -1
        f( j ) = f( j ) + q( j ) * f( j+1 )
      end do

    end subroutine tridiag1

    !**************************************************************************

!******************************************************************************
!!$
!!$    subroutine twostreamapp_vis( cossza, gph, q, gdf, &
!!$         dod067, qerat, ssa, af, sf, albedo, &
!!$         gor, goru, gord, gsr, gsru, gsrd, &
!!$         ijs, ije )
!!$
!!$      use matype
!!$      use maparam, only : im => imax, jm => jmax, km => kmax
!!$
!!$      real(dp)    , intent(in ) :: cossza( ijs:ije, 1 )
!!$      real(dp)    , intent(in ) :: gph   ( im, jm, km+1 )
!!$      real(dp)    , intent(out) :: q     ( im, jm, km   )
!!$      real(dp)    , intent(out) :: gdf   ( im, jm )
!!$      real(dp)    , intent(in ) :: dod067( im, jm, km+1 )
!!$      real(dp)    , intent(in ) :: qerat, ssa, af
!!$      real(dp)    , intent(in ) :: sf
!!$      real(dp)    , intent(in ) :: albedo( im, jm )
!!$
!!$      real(dp)    , intent(out) :: &
!!$           gor ( im, jm ), goru ( im, jm ), gord ( im, jm ), &
!!$           gsr ( im, jm ), gsru ( im, jm ), gsrd ( im, jm )
!!$
!!$      integer(i4b), intent(in ) :: ijs, ije
!!$
!!$
!!$      !
!!$      ! local variables
!!$      !
!!$      real(dp)     :: gt  ( im, jm, km )
!!$      real(dp)     :: gts ( im, jm )
!!$      real(dp)     :: emis( im, jm )
!!$      real(dp)     :: wn1 = 1.0d100, wn2 = 1.0d100
!!$      integer(i4b) :: divnum = 3
!!$      integer(i4b) :: sw = 1
!!$
!!$      integer(i4b) :: ij, k
!!$
!!$
!!$      do k = 1, km
!!$         do ij = ijs, ije
!!$            gt  ( ij, 1, k ) = 1.0d100
!!$         end do
!!$      end do
!!$      do ij = ijs, ije
!!$         gts ( ij, 1 ) = 1.0d100
!!$         emis( ij, 1 ) = 1.0d100
!!$      end do
!!$
!!$      call twostreamapp( cossza, gt, gts, gph, q, gdf, &
!!$           dod067, qerat, ssa, af, sf, albedo, emis, wn1, wn2, divnum, sw, &
!!$           gor, goru, gord, gsr, gsru, gsrd, &
!!$           ijs, ije )
!!$
!!$
!!$    end subroutine twostreamapp_vis
!!$
!!$!******************************************************************************
!!$
!!$    subroutine twostreamapp_ir( gt, gts, gph, q, gdf, &
!!$         dod067, qerat, ssa, af, emis, wn1, wn2, divnum, &
!!$         gor, goru, gord, gsr, gsru, gsrd, &
!!$         ijs, ije )
!!$
!!$      use matype
!!$      use maparam, only : im => imax, jm => jmax, km => kmax
!!$
!!$      implicit none
!!$
!!$      real(dp)    , intent(in ) :: gt    ( im, jm, km   )
!!$      real(dp)    , intent(in ) :: gts   ( im, jm )
!!$      real(dp)    , intent(in ) :: gph   ( im, jm, km+1 )
!!$      real(dp)    , intent(out) :: q     ( im, jm,   km )
!!$      real(dp)    , intent(out) :: gdf   ( im, jm )
!!$      real(dp)    , intent(in ) :: dod067( im, jm, km+1 )
!!$      real(dp)    , intent(in ) :: qerat, ssa, af
!!$      real(dp)    , intent(in ) :: emis  ( im, jm )
!!$      real(dp)    , intent(in ) :: wn1, wn2
!!$      integer(i4b), intent(in ) :: divnum
!!$
!!$      real(dp)    , intent(out) :: &
!!$           gor ( im, jm ), goru ( im, jm ), gord ( im, jm ), &
!!$           gsr ( im, jm ), gsru ( im, jm ), gsrd ( im, jm )
!!$
!!$      integer(i4b), intent(in ) :: ijs, ije
!!$
!!$
!!$      !
!!$      ! local variables
!!$      !
!!$      real(dp)     :: cossza( ijs:ije, 1 )
!!$      real(dp)     :: sf    = 1.0d100
!!$      real(dp)     :: albedo( im, jm )
!!$      integer(i4b) :: sw = 2
!!$
!!$      integer(i4b) :: ij
!!$
!!$
!!$      do ij = ijs, ije
!!$         albedo( ij, 1 ) = 1.0d0 - emis( ij, 1 )
!!$         cossza( ij, 1 ) = 1.0d100
!!$      end do
!!$
!!$      call twostreamapp( cossza, gt, gts, gph, q, gdf, &
!!$           dod067, qerat, ssa, af, sf, albedo, emis, wn1, wn2, divnum, sw, &
!!$           gor, goru, gord, gsr, gsru, gsrd, &
!!$           ijs, ije )
!!$
!!$
!!$    end subroutine twostreamapp_ir

  end module radiation_two_stream_app
