solar_irr_PC.f

Go to the documentation of this file.

********************************************************************************
*                                          *
*  Name: SOLAR_IRR                                 *
*  Purpose:    obtain solar irradiances above the atmosphere corresponding     *
*              to the measurement time and geographic location.            *
*  Parameters: none                                *
*  Global used:   WAVOBS - wavelengths of input imaging spectrometer data      *
*                 FWHM   - FWHM for each of the channels               *
*                 SOLZNI - solar zenith angle                      *
*                 IDAY   - the day number in the year when the measurement     *
*                          was made                            *
*  Global output: YIRR   - Solar irradiances just above the atmosphere         *
*  Return Codes:  none                                 *
*  Special Considerations:                             *
*       1. The solar irradiance curve is contained in the file 'sun_binary'.   *
*          It is now based on the "thuillier_atlas3.dat" in the 0.2 - 2.4      *
*          micron range (Thuillier, G., et al., Solar irradiance reference     *
*          spectra for two solar active levels, Adv. in Space Research, 34,    *
*          256-261, 2004). Above 2.4 micron, the solar curve contained         *
*          in the data file "SUN01kurucz2005.dat" inside the MODTRAN 5.2 code  *
*          & multiplied by a factor of 1.02 is adopted. Previously, ATREMi     *
*          used a solar curve converted from the file "sun_kur" in MODTRAN3.5. *
*       2. The solar irradiance curve must have enough entries starting at     *
*          about .30 so that the interpolation and smoothing will work         *
*          correctly.                                  *
*       3. This subroutine doesn't contain the statement INCLUDE 'COMMONS_INC'.*
*          As a result, the variables, such as NP_HI, NP_MED, WAVNO_HI,        *
*          TRAN_HI, WAVLN_MED, TRAN_MED, and DWAVNO are all local variables    *
*          used for smoothing the 1 cm-1 resolution solar irradiance data from *
*          MODTRAN3.5 to the desired resolutions of imaging spectrometers.     *
*          These variables have different meanings in most of other            *
*          subroutines.                                    *
*                                          *
********************************************************************************
      SUBROUTINE  solar_irr_pc
 
C  Declarations for common variables
      dimension wavobs(1024),fwhm(1024)
      COMMON /getinput4/ wavobs,fwhm
 
      dimension yirr(1024)
      COMMON /solar_irr1/yirr
 
      COMMON /getinput5/ nobs,hsurf,dlt,dlt2
      COMMON /geometry1/ solzni,solaz,obszni,obsphi,iday
 
 
C     Arrays for wavelength positions and FWHM of measured imaging spectrometer
C         data and for smoothing the medium resolution spectrum (FWHM = 0.2 nm,
C         and point spacing = 0.1 nm) to match coarser resolution spectrum
C         of imaging spectrometer data
 
      parameter(nobs_max = 1024, ninstr_max = 3001)
 
      REAL    FINSTR(NINSTR_MAX)
      INTEGER NCVHF(NOBS_MAX)
 
 
      parameter(np_hi=49933)      !Number of solar irradiance spectral points,
                                   !  covering 51 cm-1 to 49,983 cm-1 at
                                   !  point spacing of 1.00 cm-1
 
      parameter(np_med = 28001)   !Number of medium resolution spectral
                                   !  points from 0.30 um to 3.1 um at 0.1 nm
                                   !  point spacing and 0.2 nm full width at
                                   !  half maximum (FWHM)
 
      REAL WAVNO_HI(NP_HI)         !Wavenumber of solar irradiance data at a
                                   !  point spacing of 1.0 cm-1 
      REAL WAVLN_MED(NP_MED)       !Wavelength of medium resolution data
                                   ! from 0.30 to 3.1 um (0.1 nm point spacing).
      REAL TRAN_HI(NP_HI)          !Solar irradiance data at 1.0 cm-1 point
                                   ! spacing
      REAL TRAN_MED(NP_MED)        !Solar irradiance data at medium resolution 
                                   ! from 0.30 to 3.1 um (0.1 nm point spacing).
 
      REAL VSTART, VEND, DWAVLN    !VSTART: Starting wavelength for calculations
                                   !VEND:   Ending wavelength for calculations
                                   !DWAVLN: Interval between wavelengths
      parameter(vstart = 0.30, vend = 3.1, dwavln = 0.0001)
 
      REAL DWAVNO
      parameter(dwavno = 1.00)  !Point spacing of input solar irradiance
                                 ! data ( 1.0 cm-1).
 
      REAL DLT_MED               ! 0.2-nm medium resolution spectrum.
      parameter(dlt_med = 0.0002)
 
      REAL FACDLT                !Factor to multiply DLT by to determine the
      parameter(facdlt = 2.0)   !  range used in the Gaussian function
 
 
C Programmer's note: CONST2=DLT*SQRT(3.1415926/CONST1)  CONST2=total
C                    area of gaussian curve = DLT*SQRT(pi/(4.0*ln(2)))
      REAL CONST1
      parameter(const1 = 2.7725887)    ! CONST1=4.0*ln(2)=2.7725887
 
      INTEGER INDEX_MED(NP_MED)
      REAL    WAVLN_MED_INDEX(NP_MED), TRAN_MED_INDEX(NP_MED)
C
      REAL    FINSTR_WAVNO(5000), FWHM_WAVNO(NP_MED)
      INTEGER NCVHF_WAVNO(NP_MED)
 
      DATA pi,dtorad /3.1415926,0.0174533/
 
C--- Temp Code ---
        OPEN(31,file='sun_binary_PC',status='old',
     & form='unformatted',access='direct',recl=4*49933)
        read(31,rec=1) (wavno_hi(i), i = 1, 49933)
        read(31,rec=2) (tran_hi(i),  i = 1, 49933)
        close(31)
C--- End of Temp Code
 
      DO i = 1, np_med
         wavln_med(i) = vstart + float(i-1)*dwavln  ! Wavelength of medium
                                                    ! resolution spectrum,
         tran_med(i)  = 1.0                         ! FWHM=.2 nm, .30-3.1 um.
      END DO
 
 
C Note: The grids of WAVNO_HI do not match the grids of 10000./WAVLN_MED.
C       INDEX_MED is a specially designed index for finding close matches
C       between  the two kinds of grids.
C
      DO i = 1, np_med
         index_med(i) = ( (10000./wavln_med(i) - 51.)/dwavno + 1.)
      END DO
C
C Note:     WAVLN_MED_INDEX(I) is very close to WAVLN_MED(I),
C       and WAVLN_MED_INDEX(I) >= WAVLN_MED(I)
C
      DO i = 1, np_med
         wavln_med_index(i) = 10000. /(float(index_med(i)-1)*dwavno
     &                               + 51.)
      END DO
 
C
      DO i = 1, np_med
         fwhm_wavno(i) = 10000.*dlt_med
     &                   /(wavln_med_index(i)*wavln_med_index(i))
         IF(fwhm_wavno(i).LT.1.0) fwhm_wavno(i) = 1.0
      END DO
C
 
      DO i = 1, np_med
         ncvhf_wavno(i) = ( facdlt * fwhm_wavno(i) / dwavno + 1.)
      END DO
 
C Initialize arrays for smoothing medium resolution spectrum (FWHM = 0.2 nm,
C         and point spacing = 0.1 nm) to coarser spectral resolution data
C         from imaging spectrometers.
 
      DO i = 1, nobs
         ncvhf(i) = ( facdlt * fwhm(i) / dwavln + 1.)
      END DO
 
 
C First stage of smoothing - smooth solar irradiance spectrum with 49,933 
C     points at a point spacing of 1 cm-1 to medium resolution spectrum
C     (resolution of 0.2 nm and point spacing of 0.1 nm) with about
C     25,000 points.
C
C     The smoothing is done in wavenumber domain. For a spectrum with a 
C     constant 0.2 nm resolution in wavelength domain, it has variable
C     resolution in wavenumber domain. This effect is properly taken care of
C     in the design of smoothing functions.
C
C     Because the input solar spectrum is in wavenumber units (cm-1), while 
C     the medium resolution spectrum is in wavelength units, the two kinds of 
C     grids do not automatically match. In order to match the grids, arrays 
C     of INDEX_MED and TRAN_MED_INDEX are specially designed. The desired 
C     medium resolution spectrum, TRAN_MED, at constant 0.1 nm point spacing 
C     and 0.2 nm resolution is obtained through linear interpolation of 
C     TRAN_MED_INDEX array. 
C
      DO 466 j = 1, np_med
 
             tran_med_index(j)  = 0.0
             ncvtot_wavno = 2 * ncvhf_wavno(j) - 1
    
             sumins = 0.0
 
          DO 560 i = ncvhf_wavno(j), ncvtot_wavno
             finstr_wavno(i) = 
     &           exp( -const1*(float(i-ncvhf_wavno(j))*dwavno
     &           /fwhm_wavno(j))**2.0)
             sumins = sumins + finstr_wavno(i)
 560      CONTINUE
 
          DO 565 i = 1, ncvhf_wavno(j)-1
             finstr_wavno(i) = finstr_wavno(ncvtot_wavno-i+1)
             sumins = sumins + finstr_wavno(i)
 565      CONTINUE
 
          sumins = sumins * dwavno
 
          DO 570 i = 1, ncvtot_wavno
             finstr_wavno(i) = finstr_wavno(i)*dwavno/sumins
 570      CONTINUE
  
 
          DO 491 k = index_med(j)-(ncvhf_wavno(j)-1), 
     &                            index_med(j)+ncvhf_wavno(j)-1
             tran_med_index(j) = tran_med_index(j) + tran_hi(k)* 
     &                     finstr_wavno(k-index_med(j)+ncvhf_wavno(j))       
 491      CONTINUE
 
 466  CONTINUE
 
C
C Linear interpolation to get TRAN_MED from TRAN_MED_INDEX:
C     (Note that WAVLN_MED_INDEX(J) >= WAVLN_MED(J)    )
C
         tran_med(1)      = tran_med_index(1)  
         tran_med(np_med) = tran_med_index(np_med)  
 
      DO j = 2, np_med-1
           dlt  =  wavln_med_index(j) - wavln_med_index(j-1)
         IF(dlt.LT.1.0e-06) THEN
           tran_med(j) = tran_med_index(j)
         ELSE
           fjm1 = (wavln_med_index(j) - wavln_med(j))        /dlt
           fj   = (wavln_med(j)       - wavln_med_index(j-1))/dlt
           tran_med(j) = fjm1*tran_med_index(j-1) + fj*tran_med_index(j)
         END IF
      END DO
 
 
C Second stage of smoothing - smooth the medium resolution spectrum (resolution 
C     of 0.2 nm and point spacing of 0.1 nm) with about 28,000 points to match
C     the coarser and variable resolution spectrum from imaging spectrometers.
C
C Initialize some index parameters:
      ia = 1000
C
      DO 1466 j =1, nobs
 
             yirr(j) = 0.0
             tran_ia     = 0.0
             tran_iap1   = 0.0
 
             ncvtot = 2 * ncvhf(j) - 1
    
C Calculate instrumental response functions...
             sumins = 0.0
 
          DO 1560 i = ncvhf(j), ncvtot
             finstr(i) = 
     &           exp( -const1*(float(i-ncvhf(j))*dwavln
     &           /fwhm(j))**2.0)
             sumins = sumins + finstr(i)
1560      CONTINUE
 
          DO 1565 i = 1, ncvhf(j)-1
             finstr(i) = finstr(ncvtot-i+1)
             sumins = sumins + finstr(i)
1565      CONTINUE
 
          sumins = sumins * dwavln
 
          DO 1570 i = 1, ncvtot
             finstr(i) = finstr(i)*dwavln/sumins
1570      CONTINUE
  
C  Index searching...
C
          CALL hunt(wavln_med, np_med, wavobs(j), ia)
 
 
C  Smoothing...
C
          DO 1491 k = ia-(ncvhf(j)-1), ia+ncvhf(j)-1
             tran_ia = tran_ia + tran_med(k)* 
     &                     finstr(k-ia+ncvhf(j))       
1491      CONTINUE
 
            ia_p1 = ia + 1
          DO 1492 k = ia_p1-(ncvhf(j)-1), ia_p1+ncvhf(j)-1
             tran_iap1 = tran_iap1 + tran_med(k)* 
     &                     finstr(k-ia_p1+ncvhf(j))       
1492      CONTINUE
C
C Linear interpolation to get YIRR from TRAN_IA and TRAN_IAP1:
C
           dlt_ia  =  wavln_med(ia_p1) - wavln_med(ia)
           fia     = (wavln_med(ia_p1) - wavobs(j)) /dlt_ia
C          FIA_P1  = (WAVOBS(J)     - WAVLN_MED(IA))/DLT_IA
           fia_p1  = 1. - fia
           yirr(j) = fia*tran_ia + fia_p1*tran_iap1
C
1466  CONTINUE
 
 
C--          DO I = 1, NOBS
C--             print*,WAVOBS(I),YIRR(I)
C--          END DO
 
      cossol=cos(solzni)
 
C The Sun-earth distance correction factor, a function of the day in the year
C (see reference: IQBAL, p.3)
 
      gamm=2.0*pi*float(iday-1)/365.
      e0=1.000110+0.034221*cos(gamm)+0.001280*sin(gamm)
     & + 0.000719*cos(2.0*gamm)+0.000077*sin(2.0*gamm)
 
      DO 20 i=1,nobs
        yirr(i)=yirr(i)*cos(solzni)*e0
   20 CONTINUE
 
C      WRITE(*,*)'SOLZNI=',SOLZNI,'COSSOL= ',COS(SOLZNI),'GAMM=',GAMM,
C     &'IDAY=',IDAY,' E0= ',E0
 
C--      WRITE(*,*)'SOLZNI=',SOLZNI,'COSSOL= ',COS(SOLZNI),'GAMM=',GAMM,
C--     &'IDAY=',IDAY,' E0= ',E0
C--      DO 23 I=1,NOBS
C--   23   WRITE(*,*)'I=',I,' WAVOBS(I)=',WAVOBS(I),' YIRR(I)=',YIRR(I),
C--     &YIRR(I)/(COS(SOLZNI)*E0)
 
 
      RETURN
      END