gomspheremain_f.htm

sphere scattering program and statment

c      METHOD :  Tabular interpolation, assuming
c                (1) real index is linear in log(wavelength)
c                    and linear in temperature
c                (2) log(imag. index) is linear in log(wavelength)
c                    and linear in temperature
c     AUTHORS OF subroutine refice( WAVMET, TEMP, ref_ice ):
c             Stephen Warren, Univ. of Washington (1983)
c               ([e-mail])
c             Bo-Cai Gao, JCESS, Univ. of Maryland (1995)
c               ([e-mail])
c             Warren Wiscombe, NASA Goddard (1995)
c               ([e-mail])
c     MODIFICATIONS IN 1995 :
c       Gao, Warren, and (to a small extent) Wiscombe modified the
c       original Warren refice program from 1984 to change values of
c       imaginary refractive index in the 0.161-0.410 and 1.445-2.50
c       micron regions.  The values in 0.161-0.410 were incorrect and
c       the values in 1.445-2.50 were among the most uncertain in 1984.
c       New measurements have made it possible to improve both regions.

c       No changes were made to real refractive indices (by re-doing a
c       Kramers-Kronig analysis), because the values of imaginary
c       index MIM involved are so small (below 0.001) that the
c       resulting changes in real index MRE would be in the third
c       decimal place at most.  (MIM has negligible effect on MRE
c       when MIM << MRE.)

c       The 0.161-0.410 micron region was changed using data provided
c       by Warren, which correct his misinterpretation of Minton's
c       measurements for 0.181-0.185 micron, and incorporate new
c       measurements of Perovich and Govoni (1991) for 0.250-0.400
c       micron.  Warren (1984) correctly represented UV measurements
c       of Seki et al. and visible measurements of Grenfell/Perovich,
c       but he plotted Minton's measurements a factor of 2.3 too low
c       because he misinterpreted base-10 as base-e.  (The UV
c       measurements of Dressler/Schnepp and Shibaguchi et al are also
c       probably expressed as absorption coefficients on base-10;
c       therefore those values also were probably plotted a factor of
c       2.3 too low in Warren's (1984) Figure 2.)

c       The details of how the present imaginary index data for
c       0.161-0.410 micron is obtained are as follows.  Point A in
c       Warren's Figure 2 at 161 nm is joined with a straight line to
c       Minton's corrected point B at 181 nm.  Minton's reported
c       values for 181-185 nm have been smoothed within his stated
c       uncertainty.  Now a smooth curve is drawn to join Minton at
c       185 nm to Perovich/Govoni (PG) at 250 nm.  PG's values from
c       their Table 1 show some unrealistic wiggles that are smaller
c       than their error bars, so a smooth curve was fitted through
c       them and values were taken from the smoothed curve at 10-nm
c       intervals.  PG ends at 400 nm, where Grenfell/Perovich (GP)
c       starts.  At 400 nm we take imaginary index=2.82E-9, the
c       average of PG (2.93E-9) and GP (2.71E-9).

c       The Warren (1984) values of imaginary index in the 1.445-2.50
c       micron region were replaced by those of Kou et al.(1994).  In
c       order to remove the resulting discontinuities near 1.445 and
c       2.5 micron, the Warren values at 1.43 and 1.44 micron were
c       changed to 0.9E-04 and 1.3E-04 respectively, and his values at
c       2.52, 2.55, and 2.565 micron were changed to 8.255E-04,
c       8.578E-04, and 8.739E-04, respectively. The latter change
c       eliminated a small local maximum at 2.5 micron which was not
c       realistic and has never been seen in spectra of snow bracketing
c       that wavelength.


c     REFERENCES :

c       Warren, S., 1984: Optical Constants of Ice from the Ultraviolet
c          to the Microwave, Appl. Opt. 23, 1206-1225

c       Kou, L., D. Labrie, and P. Chylek, 1994: Refractive indices
c          of water and ice in the 0.65- to 2.5-micron spectral range,
c          Appl. Opt. 32, 3531-3540

c       Perovich, D., and J. Govoni, 1991: Absorption Coefficients
c          of Ice from 250 to 400 nm, Geophys. Res. Lett. 18, 1233-1235
c ======================================================================

      IMPLICIT NONE

c     .. Parameters ..

      INTEGER   NWL, NWLT
      PARAMETER ( NWL = 574, NWLT = 62)
c     ..
c     .. Scalar Arguments ..

      REAL*8 TEMP, WAVLEN,WAVMET
      complex*16 ref_ice
c     ..
c     .. Local Scalars ..

      CHARACTER MESSAG*40
      LOGICAL   PASS1
      INTEGER   I, L
      REAL*8      FRAC, MIM, MRE, YHI, YLO
c     ..
c     .. Intrinsic Functions ..

      INTRINSIC LOG, CMPLX, EXP
c     ..
c                                        ** Refractive index table

      REAL*8      TABIM( NWL ), TABIMT( NWLT, 4 ), TABRE( NWL ),
     &          TABRET( NWLT, 4), TEMREF(4), WL( NWL ), WLT( NWLT )

      COMMON / ICEREF / WL, WLT, TABRE, TABRET, TABIM, TABIMT, TEMREF

      DATA  PASS1 / .True. /

       WAVLEN = WAVMET*1.0d+6
       IF( PASS1 )  THEN

         PASS1 = .False.
c                         ** Superficially test if main table messed up

         IF( NWL.LT.100 ) CALL ERRMSG('refice--NWL value bad',.True.)
         IF( WL(1).GT.0.045 ) CALL ERRMSG('refice--WL(1) bad',.True.)
         IF( WL(NWL).LT.166.) CALL ERRMSG('refice--WL(NWL) bad',.True.)

         DO 1  I = 1, NWL

            IF( WL(I).LE.0.0 .OR. TABRE(I).LE.0.5 .OR. TABRE(I).GT.2.0
     &          .OR. TABIM(I).LT.0.0  .OR. TABIM(I).GT.10.0 )  THEN
               WRITE( MESSAG, '(A,I5,A)' )  'refice--table value ', I,
     &                ' out of bounds '
               CALL ERRMSG( MESSAG, .True. )
            END IF

            IF( I.GT.1 .AND. WL(I).LE.WL(I-1) )  THEN
               WRITE( MESSAG, '(A,I5,A)' )  'refice--table WL(', I,
     &                ') not increasing  '
               CALL ERRMSG( MESSAG, .True. )
            END IF

    1    CONTINUE

      END IF


      IF( WAVLEN.LT.WL(1) .OR. WAVLEN.GT.WLT(NWLT) ) THEN

         CALL ERRMSG('refice--wavelength outside table boundaries',
     &                .False.)
         ref_ice = (0.,0.)
         RETURN

      END IF


      IF( WAVLEN.LE.167.) THEN
c                                  ** Wavelength between 0.045 and 167
c                                  ** microns. No temperature dependence
         DO 10 I = 2, NWL
            IF( WAVLEN.LE.WL(I)) GO TO 20
   10    CONTINUE

   20    CONTINUE
         FRAC   = LOG( WAVLEN / WL(I-1) ) / 
     &            LOG( WL(I)  / WL(I-1) )
         MRE    = TABRE(I-1) + FRAC * ( TABRE(I) - TABRE(I-1) )
         MIM    = TABIM(I-1) * ( TABIM(I) / TABIM(I-1) )**FRAC


      ELSE

c               ** Wavelength greater than 167 microns
c               ** (temperature-dependent case)

c        write (*,*)'temp=',TEMP

        IF( TEMP.LT.TEMREF(4) .OR. TEMP.GT.TEMREF(1) ) THEN

            CALL ERRMSG('refice--temperature outside table boundaries',
     &                  .False.)
            ref_ice = (0.,0.)
            RETURN

         END IF
c                         ** Find position in temperature array
         DO 30 L = 2, 4
            IF( TEMP.GE.TEMREF(L) ) GO TO 40
   30    CONTINUE
c                         ** Find position in wavelength array
   40    CONTINUE
         DO 50 I = 2, NWLT
            IF( WAVLEN.LE.WLT(I) ) GO TO 60
   50    CONTINUE

   60    CONTINUE
         FRAC   = LOG( WAVLEN / WLT(I-1) ) /
     &            LOG( WLT(I) / WLT(I-1) )

         YLO    = TABRET(I-1, L) +
     &            FRAC*( TABRET(I, L) - TABRET(I-1, L) )

         YHI    = TABRET( I-1, L-1) +
     &            FRAC*( TABRET(I, L-1) - TABRET(I-1, L-1) )

         MRE    = YLO + ( YHI - YLO) * ( TEMP - TEMREF(L) ) /
     &                                 ( TEMREF(L-1) - TEMREF(L) )

         YLO    = LOG( TABIMT(I-1, L)) +
     &            FRAC*LOG( TABIMT(I, L) / TABIMT(I-1, L) )

         YHI    = LOG( TABIMT(I-1, L-1) ) +
     &            FRAC*LOG( TABIMT(I, L-1) / TABIMT(I-1, L-1) )

         MIM    = EXP( YLO + (YHI - YLO) * (TEMP - TEMREF(L)) /
     &                                     (TEMREF(L-1) - TEMREF(L)) )

      END IF


      ref_ice = CMPLX( MRE, MIM)

      END

      BLOCK DATA ICECON

c        Ice-refractive-index vs. wavelength table for subroutine refice

      IMPLICIT NONE

c     .. Parameters ..

      INTEGER   NWL, NWLT
      PARAMETER ( NWL = 574, NWLT = 62 )
c     ..
c     .. Local Scalars ..

      INTEGER   I
c     ..

      REAL*8      TABIM(NWL), TABIMT( NWLT, 4 ), TABRE(NWL),
     &          TABRET( NWLT, 4 ), TEMREF(4), WL(NWL), WLT(NWLT)

      COMMON / ICEREF / WL, WLT, TABRE, TABRET, TABIM, TABIMT, TEMREF


c       WL, WLT           wavelengths (microns) for temperature-
c                         independent and temperature-dependent
c                         regimes, respectively

c       TABRE, TABRET     real refractive indices for temperature-
c                         independent and temperature-dependent
c                         regimes, respectively

c       TABIM, TABIMT     imaginary refractive indices for temperature-
c                         independent and temperature-dependent
c                         regimes, respectively

c       TEMREF            reference temperatures (-1,-5,-20,-60 deg C)
c                         for TABRET,TABIMT


      DATA  TEMREF / 272.16, 268.16, 253.16, 213.16 /

      DATA  ( WL( I ), TABRE( I ), TABIM( I ), I = 1, 10 ) /
     & 4.43000E-02, 8.34410E-01, 1.64000E-01,
     & 4.51000E-02, 8.36760E-01, 1.73000E-01,
     & 4.59000E-02, 8.37290E-01, 1.83000E-01,
     & 4.68000E-02, 8.37710E-01, 1.95000E-01,
     & 4.77000E-02, 8.38270E-01, 2.08000E-01,
     & 4.86000E-02, 8.40380E-01, 2.23000E-01,
     & 4.96000E-02, 8.47190E-01, 2.40000E-01,
     & 5.06000E-02, 8.55220E-01, 2.50000E-01,
     & 5.17000E-02, 8.60470E-01, 2.59000E-01,
     & 5.28000E-02, 8.62480E-01, 2.68000E-01 /

      DATA  ( WL( I ), TABRE( I ), TABIM( I ) , I = 11, 20 ) /
     & 5.39000E-02, 8.61570E-01, 2.79000E-01,
     & 5.51000E-02, 8.60930E-01, 2.97000E-01,
     & 5.64000E-02, 8.64190E-01, 3.19000E-01,
     & 5.77000E-02, 8.69160E-01, 3.40000E-01,
     & 5.90000E-02, 8.77640E-01, 3.66000E-01,
     & 6.05000E-02, 8.92960E-01, 3.92000E-01,
     & 6.20000E-02, 9.10410E-01, 4.16000E-01,
     & 6.36000E-02, 9.30890E-01, 4.40000E-01,
     & 6.53000E-02, 9.53730E-01, 4.64000E-01,
     & 6.70000E-02, 9.81880E-01, 4.92000E-01 /

      DATA  ( WL( I ), TABRE( I ), TABIM( I ) , I = 21, 30 ) /
     & 6.89000E-02, 1.02334, 5.17000E-01,
     & 7.08000E-02, 1.06735, 5.28000E-01,
     & 7.29000E-02, 1.11197, 5.33000E-01,
     & 7.38000E-02, 1.13134, 5.34000E-01,
     & 7.51000E-02, 1.15747, 5.31000E-01,
     & 7.75000E-02, 1.20045, 5.24000E-01,
     & 8.00000E-02, 1.23840, 5.10000E-01,
     & 8.27000E-02, 1.27325, 5.00000E-01,
     & 8.55000E-02, 1.32157, 4.99000E-01,
     & 8.86000E-02, 1.38958, 4.68000E-01 /

      DATA  ( WL( I ), TABRE( I ), TABIM( I ) , I = 31, 40 ) /
     & 9.18000E-02, 1.41644, 3.80000E-01,
     & 9.30000E-02, 1.40906, 3.60000E-01,
     & 9.54000E-02, 1.40063, 3.39000E-01,
     & 9.92000E-02, 1.40169, 3.18000E-01,
     & 1.03300E-01, 1.40934, 2.91000E-01,
     & 1.07800E-01, 1.40221, 2.51000E-01,
     & 1.10000E-01, 1.39240, 2.44000E-01,
     & 1.12700E-01, 1.38424, 2.39000E-01,
     & 1.14000E-01, 1.38075, 2.39000E-01,
     & 1.18100E-01, 1.38186, 2.44000E-01 /

      DATA  ( WL( I ), TABRE( I ), TABIM( I ) , I = 41, 50 ) /
     & 1.21000E-01, 1.39634, 2.47000E-01,
     & 1.24000E-01, 1.40918, 2.24000E-01,
     & 1.27200E-01, 1.40256, 1.95000E-01,
     & 1.29500E-01, 1.38013, 1.74000E-01,