      SUBROUTINE SPHERIC1(HEAT,PMB,T,CMO3,QOF,RMU0,
     1                    XVF,O2VF,AMUVF,CPVF,GVF,EFFH,EFFSRC,TO2,
     2                    KMX,KMXP,LMX,IL1,IL2,
     3                    TO3,ALT,HMTO3,HMTO2,INDX) 
C***********************************************************************
C*                                                                     *
C*              SUBROUTINE SPHERIC1                                    *
C*                                                                     *
C***********************************************************************
C
C to calculate the solar heating due to O2 and O3 taking into account
C effect of sphericity (sun below the horizont).
C Used parameterizations:
C - Strobel (JRG, vol 83, p 6225, 1978) for
C   O3 bands: Hartley (242.5-277.5 nm), Huggins (277.5-360 nm),
C             Herzberg (206-242.5 nm)
C   O2 spectral regions: Herzberg contiuum ((206-242.5 nm),
C                        Schumann-Runge bands (175-205) & continuum (125-175)
C - Shine and Rickaby (in "Ozone in the Atmosphere", edited by R.D. Bojkov and
C                    P. Fabian, 597--600, A. Deepak, Hampton, Virginia, 1989
C   for O3 Chappuis band (407.5-852.5 nm)
C - Mlynczack&Solomon (JGR, vol 98, p 10517, 1993) for
C   heating efficiency in O2 SRC and O3 region shorter than 305.5 nm
C - Rodgers formula for atmospheric mass
C                                               V Fomichev, October, 1999

C Called by MAMRAD
C Calls nothing

      IMPLICIT REAL   (A-H,O-Z),
     +INTEGER (I-N)

C output array (heating in K/s)

      real HEAT(LMX,KMX)

C input arrays: pressure (mb) and temperature (K) at layer's interface,
C               cosine of solar zenith angle, O3 mass mixing ratio at
C               mid-layer and O3 content in a layer (cm-atm).

      REAL   PMB(LMX,KMXP),T(LMX,KMX),CMO3(LMX,KMX),QOF(LMX,KMX)

      REAL   RMU0(LMX)

      REAL   XVF(KMX),O2VF(KMX),AMUVF(KMX),CPVF(KMX),GVF(KMX)
      REAL   EFFH(KMX),EFFSRC(KMX),TO2(KMX)

C     * INTERNAL ARRAYS

      REAL   TO3(LMX,KMX), ALT(LMX,KMX), HMTO3(LMX,KMX), HMTO2(LMX,KMX)
      INTEGER INDX(LMX,KMX)
C------------------------------------------------------------------
c calculate geometric altitude (in meters) from the center of Earth up
c
      do 10 i=IL1,IL2
        hscal=8.31441E3*T(i,1)/(AMuvf(1)*Gvf(1))
        ALT(i,1)=hscal*0.5*log(PMB(i,1)/PMB(i,2))+6.37E6
   10 CONTINUE
c
      do 30 k=2,KMX
        am=0.5*(AMuvf(k)+AMuvf(k-1))
        grav=0.5*(Gvf(k)+Gvf(k-1))
        km=k-1
        kp=k+1
        do 20 i=IL1,IL2
          hscal=8.31441E3*T(i,k)/(am*grav)
          ALT(i,k)=ALT(i,km)+hscal*0.5*log(PMB(i,km)/PMB(i,kp))
   20   continue
   30 continue
c
c calculate O3 column amount (in cm-2) from bottom up
c 
      do 40 i=IL1,IL2
        tO3(i,KMX)=2.687E19*0.5*QOF(i,KMX)
   40 continue
c
      do 60 k=2,KMX
        kr=KMXP-k
        krp=kr+1
        do 50 i=IL1,IL2
          tO3(i,kr)=tO3(i,krp)+2.687E19*0.5*(QOF(i,kr)+QOF(i,krp))
   50   continue
   60 continue
C
C     * CALCULATE NECESSARY ARRAYS TO DETERMINE COLUMN O2,O3 AMOUNTS
C     * ABOVE THE TANGENT HEIGHT.
C     * THIS LOOP WILL NOT VECTORIZE BECAUSE OF COMPUTED GO TO'S.
C     * PERHAPS LOOK AT REVISING THIS FOR NEXT VERSION.
C
      DO 70 K=1,KMX
      DO 70 IL=IL1,IL2
C
C       * TANGENT HEIGHT.
C
        HM=ALT(IL,K)*sqrt(1.-RMU0(IL)*RMU0(IL))

C       * CALCULATIONS ARE TO BE DONE ONLY IF TANGENT HEIGHT ABOVE THE
C       * SURFACE BY 10KM

        IF(HM.GT.6.38E6) THEN
          INDX(IL,K)=1  

C         * TO FIND THE COLUMN O2 AND O3 AMOUNT ABOVE THE TANGENT HEIGHT
C         * LINEAR INTERPOLATION.

          kmin=1
          kmax=K
    1     km=int((kmin+kmax)/2.)
          if(km.eq.kmin) then
            goto 4
          else
            if(hm-ALT(IL,km)) 2,3,3
    2       kmax=km
            goto 1
    3       kmin=km
            goto 1
          end if
    4     continue

c         * 'km' has been found so that ALT(IL,km).le. hm .lt.ALT(IL,km+1).

          hvsp=(hm-ALT(IL,km+1))/(ALT(IL,km+1)-ALT(IL,km))
          hmtO3(il,k)=(tO3(IL,km+1)-tO3(IL,km))*hvsp + tO3(IL,km+1)
          hmtO2(il,k)=(tO2(km+1)-tO2(km))*hvsp + tO2(km+1)
        ELSE
          INDX(IL,K)=0
        ENDIF
   70 CONTINUE

C CALCULATE SOLAR O3 AND O2 HEATING

C FIRST, INITIALIZE FACTOR MULTIPLYING SRB AND SR CONTINUUM HEATINGS
C DEPENDING ON WHERE MODEL TOP IS.

      if(xvf(kmx).gt.8.3) then
         srfact=1.
      else
         srfact=0.
      endif
C
      DO 90 K = 1, KMX

c efficiency for O3 heating at wavelength shorther than 305.5 nm

      effHa=(effH(k)-0.26)/0.74

      DO 80 IL = IL1, IL2

      IF (RMU0(IL) .gt. 0.0) THEN

c **************** SUN ABOVE THE HORIZONT ****************

c O3 and O2 column amount: downward looking paths
        gO3=tO3(IL,K)*35./SQRT(1224.0*RMU0(IL)*RMU0(IL)+1.)
        gO2=tO2(K)*35./SQRT(1224.0*RMU0(IL)*RMU0(IL)+1.)

c O3 bands (CH - after Shine & Rickaby, others - after Strobel):

        CH=4.1847E8*exp(-3.157E-21*gO3)*(0.4967+
     *                                   0.5033)              ! chem.pot
        eHU1=exp(-1.7658E-19*gO3)
        HU=(1.6935E27*(1-eHU1)*(0.75+
     *                          0.25)                         ! chem.pot
     *     +1.1443E27*(eHU1-exp(-4.2249E-18*gO3))*(0.75*effHa+
     *                                         0.25))/gO3     ! chem.pot
        HA=1.4707E10*exp(-8.8E-18*gO3)*(0.7903*effHa+
     *                                  0.2097)               ! chem.pot
        eHZ=exp(-6.6E-24*gO2-4.9E-18*gO3)
        HZ_O3=2.1363E9*eHZ*(0.8079*effHa+
     *                      0.1921)                           ! chem.pot

        heatO3=(CH+HU+HA+HZ_O3)*cmO3(IL,K)

c O2: SRB, and SR & Hz continuums (after Strobel, only net heating for SR):
c NOTE: No SRC and SRB in the stratosphere since conditions for wavelengths < 205nm
c       do not exist!

        HZ_O2=4.7695E3*eHZ*(0.0633+
     *                      0.9367)                           ! chem.pot
        SRC1=2.716E6*exp(-1.E-17*gO2)
        SRC2=(5.902E23*exp(-2.9E-19*gO2)-3.312E23*exp(-1.7E-18*gO2)
     *                               -2.590E23*exp(-1.15E-17*gO2))/gO2
        SRC = (SRC1+SRC2)*effSRC(k)
        if(gO2.gt.1.e18) then
          SRB=1./(1.113E-24*gO2+5.712E-15*sqrt(gO2))
        else
          SRB=1.463E5
        end if

        heatO2=(srfact*(SRC+SRB)+HZ_O2)*O2vf(K)

        HEAT(IL,K)=(heatO3+heatO2)/(AMuvf(K)*Cpvf(K))

      ELSE

c       **************** SUN BELOW THE HORIZONT ****************


c       * calculations are to be done only if tangent height above
c       * the surface by 10km (see where "indx" calculated in earlier
c       * loop)

        IF(INDX(IL,K).EQ.1)                                  THEN


c O3 and O2 amount along the path from atm.lid to a given point

          gO3=35.*(2.*hmtO3(il,k)-tO3(IL,K))
          gO2=35.*(2.*hmtO2(il,k)-tO2(K))

c O3 bands (CH - after Shine & Rickaby, others - after Strobel):

          CH=4.1847E8*exp(-3.157E-21*gO3)*(0.4967+
     *                                     0.5033)            ! chem.pot
          eHU1=exp(-1.7658E-19*gO3)
          HU=(1.6935E27*(1-eHU1)*(0.75+
     *                            0.25)                       ! chem.pot
     *       +1.1443E27*(eHU1-exp(-4.2249E-18*gO3))*(0.75*effHa+
     *                                           0.25))/gO3   ! chem.pot
          HA=1.4707E10*exp(-8.8E-18*gO3)*(0.7903*effHa+
     *                                    0.2097)             ! chem.pot
          eHZ=exp(-6.6E-24*gO2-4.9E-18*gO3)
          HZ_O3=2.1363E9*eHZ*(0.8079*effHa+
     *                        0.1921)                         ! chem.pot

          heatO3=(CH+HU+HA+HZ_O3)*cmO3(IL,K)

c O2: SRB, and SR & Hz continuums (after Strobel, only net heating for SR):
c NOTE: No SRC and SRB in the stratosphere since conditions for wavelengths < 205nm
c       do not exist!

          HZ_O2=4.7695E3*eHZ*(0.0633+
     *                        0.9367)                         ! chem.pot
          SRC1=2.716E6*exp(-1.E-17*gO2)
          SRC2=(5.902E23*exp(-2.9E-19*gO2)-3.312E23*exp(-1.7E-18*gO2)
     *                                 -2.590E23*exp(-1.15E-17*gO2))/gO2
          SRC = (SRC1+SRC2)*effSRC(k)
          if(gO2.gt.1.e18) then
            SRB=1./(1.113E-24*gO2+5.712E-15*sqrt(gO2))
          else
            SRB=1.463E5
          end if

          heatO2=(srfact*(SRC+SRB)+HZ_O2)*O2vf(K)

          HEAT(IL,K)=(heatO3+heatO2)/(AMuvf(K)*Cpvf(K))

        ENDIF                                                 ! indx
      ENDIF                                                   ! RMU0

   80 CONTINUE
   90 CONTINUE

      return
      end
