SUBROUTINE srb_lya_heating (sun_lya,sun_srb,o2,klev,papm1,paphm1,amu0_x,ptte,am)
! 
!   Calculates heating by Schumann Runge bands (175-200 nm) and by Lyman 
!
!   alpha line (121.6 nm) not taken into account by MA-ECHAM5 and 
!   Parameterization based on Zhu(1994) and Nicolet(1985).
!   
!   *srb_lya_heating* is called from *physc*.
!
!   Subroutine arguments:
!
!    real sun_lya,sun_srb,paphm1(klev+1),amu0_x
!    real, dimension (klev) :: o2,papm1,ptte
    INTEGER, PARAMETER      :: dp=4
    INTEGER, INTENT(in)     :: klev
    REAL(dp), INTENT(in)    :: papm1(klev)    ! level pressure (centers) Pa
    REAL(dp), INTENT(IN)    :: paphm1(klev+1) ! half level pressure (edges) Pa
    REAL(dp), INTENT(IN)    :: amu0_x         ! optical mass (Chapman function)
    REAL(dp), INTENT(IN)    :: o2(klev)       ! molecular oxygen number density
    REAL(dp), INTENT(inout) :: ptte(klev)     ! heating rate
    REAL(dp), INTENT(in)    :: am(klev)     ! heating rate
!
!   Local constants and variables:
!
    INTEGER :: jk  
    REAL(dp), PARAMETER :: srb_lya_heating_lowbound = 1100_dp ! Apply parame-
                                                              ! terization for levels
                                                              ! above 1100 Pa (30 km)
    REAL(dp), PARAMETER :: vclo2_above = 2.4E19_dp ! Part from above the model 
                                                   ! top  
    REAL(dp), PARAMETER :: srb_a = 0.143_dp        ! Schumann-Runge constant
    REAL(dp), PARAMETER :: srb_b = 9.64E08_dp      ! Schumann-Runge constant
    REAL(dp), PARAMETER :: sigma_lya = 1.E-20_dp   ! Absorption cross section
                                                   ! of O2 for Lyman alpha line 
                                                   ! [cm2]
    REAL(dp) :: srb_heating, lya_heating               
    REAL(dp) :: vclo2(klev), yy1, yy2
    
    REAL(dp) :: vmr_o2, heat_cap, grav
    
    vmr_o2   = 0.21            ! oxygen VMR
    heat_cap = 1005            ! air heat capacity [J/kg/K]
    grav     = 9.80665         ! gravity accel [m/s**2]
    amo2     = 32.0            ! O2 molecular weight 
    amo3     = 48.0            ! O3 molecular weight 
    amd      = 28.97           ! air molecular weight g/mol
    avonumber= 6.022137E23     ! Avogadro Number 
    api      = 3.141592        ! Pi number
    avoinv=1./avonumber
!
!   Calculations only for positive zenith angles:
!
!       print *,' Klev=',klev,' Optical mass',amu0_x
       IF (amu0_x .LT. 35._dp) THEN         

          vclo2(1) = vclo2_above + o2(1)*(papm1(1)-paphm1(1))/(grav*amd*avoinv*am(1))*amu0_x*1.0e-1_dp          !slant mol/cm**2
!          print *,1,vclo2(1),(papm1(1)-paphm1(1))/(grav*amd*avoinv*am(1))*1.0e-1_dp
          DO jk=2,klev
!
!            Calculations of heating due to oxygen absorption only above ~30 km:
!
             IF (papm1(jk) .LE. srb_lya_heating_lowbound) THEN
                vclo2(jk) = vclo2(jk-1) + o2(jk)*(papm1(jk)-papm1(jk-1))/(grav*amd*avoinv*am(jk))*amu0_x*1.0e-1_dp  !slant mol/cm**2
!                print *,jk,vclo2(jk),(papm1(jk)-papm1(jk-1))/(grav*amd*avoinv*am(jk))*1.0e-1_dp
             ENDIF
          ENDDO
          DO jk=1,klev         
             IF (papm1(jk) .LE. srb_lya_heating_lowbound) THEN
                lya_heating = vmr_o2*sun_lya*(1.725E-18_dp/(vclo2(jk)**0.1175_dp))*1e-1_dp/(avoinv*amd)/heat_cap*exp(-2.115E-18_dp*vclo2(jk)**0.8855_dp)  ! K/s
!
!               Heating by Schumann Runge bands:
!
                yy1 = sqrt(1.0_dp+4.0_dp*2.07e-20_dp*vclo2(jk)/api/0.0152_dp)                                  ! cm2*mol/cm2
                yy2 = (vclo2(1)/vclo2(jk))**0.3_dp*2.07e-24_dp                                                 ! cm2/cm2*m2 = m2
                ro=amd*avoinv*am(jk)*1.e3                                                                      ! kg/m3
                srb_heating = o2(jk)*1.e6*sun_srb*yy2/heat_cap/yy1*exp((-api*0.0152_dp/2._dp)*(yy1-1.0_dp))/ro ! mol/m3*w/m2*m2/J*kg*K/kg*m3 = kg*K/m3*J/s/J = K/s
!                print *,jk,ro,lya_heating*86400.,srb_heating*86400.
!             
!               Total heating [K/s]:
!
                ptte(jk) = srb_heating + lya_heating
             ENDIF ! Layer limitations            
          ENDDO ! klev
!pause
       ENDIF   ! daytime
  END SUBROUTINE srb_lya_heating