* This file contains the following subroutine, related to setting up
* grids for numerical calculations:
*     gridw
*     gridz
*     buildz
*     gridt
*     gridck
*=============================================================================*

      SUBROUTINE gridw(wstart, wstop, nwint,
     $     nw,wl,wc,wu)

*-----------------------------------------------------------------------------*
*=  PURPOSE:                                                                 =*
*=  Create the wavelength grid for all interpolations and radiative transfer =*
*=  calculations.  Grid may be irregularly spaced.  Wavelengths are in nm.   =*
*=  No gaps are allowed within the wavelength grid.                          =*
*-----------------------------------------------------------------------------*
*=  PARAMETERS:                                                              =*
*=  NW  - INTEGER, number of wavelength grid _points_                     (O)=*
*=  WL  - REAL, vector carrying the lower limit of each wavel. interval   (O)=*
*=  WC  - REAL, vector carrying the center wavel of each wavel. interval  (O)=*
*=              (wc(i) = 0.5*(wl(i)+wu(i), i = 1..NW-1)                      =*
*=  WU  - REAL, vector carrying the upper limit of each wavel. interval   (O)=*
*=
*=  MOPT- INTEGER OPTION for wave-length IF 3 good for JO2                (O)=*
*-----------------------------------------------------------------------------*

      IMPLICIT NONE

      INCLUDE 'params'

* input:

      REAL wstart, wstop
      INTEGER nwint

* output:

      REAL wl(kw), wc(kw), wu(kw)
      INTEGER nw

      integer wn(kw)

* local:

      INTEGER mopt
      REAL wincr
      INTEGER iw, i

      CHARACTER*40 fi
      CHARACTER*20 wlabel

      REAL airout
      INTEGER mrefr

      REAL dum

      LOGICAL ok
*_______________________________________________________________________

**** chose wavelength grid

* some pre-set options
*     mopt = 1    equal spacing
*     mopt = 2    grid defined in data table
*     mopt = 3    user-defined
*     mopt = 4    fast-TUV, troposheric wavelengths only
*     mopt = 5    high resolution grid for O3 isotopologue study
*     mopt = 6    create uniform grid in air-wavelength scale

*     mopt = 10  Landgraf and Crutzen, 1998
*     mopt = 11  fastJ, Wild et al. 2000
*     mopt = 12  fastJ2, Bian and Prather, 2002
*     mopt = 13  UV-b, UV-a, Visible

      mopt = 1

      IF(nwint .EQ. -156) mopt = 2
      IF(nwint .EQ. -7) mopt = 4

      IF(nwint .eq. -10) mopt = 10
      IF(nwint .eq. -11) mopt = 11
      IF(nwint .eq. -12) mopt = 12
      IF(nwint .eq. -13) mopt = 13


      IF (mopt .EQ. 1) GO TO 1
      IF (mopt .EQ. 2) GO TO 2
      IF (mopt .EQ. 3) GO TO 3
      IF (mopt .EQ. 4) GO TO 4
      IF (mopt .EQ. 5) GO TO 5
      IF (mopt .EQ. 6) GO TO 6

      IF (mopt .EQ. 10) GO TO 10
      IF (mopt .EQ. 11) GO TO 11
      IF (mopt .EQ. 12) GO TO 12
      IF (mopt .EQ. 13) GO TO 13

*_______________________________________________________________________

 1    CONTINUE

      wlabel = 'equal spacing'
      nw = nwint + 1
      wincr = (wstop - wstart) / FLOAT (nwint)
      DO iw = 1, nw-1
         wl(iw) = wstart + wincr*FLOAT(iw-1)
         wu(iw) = wl(iw) + wincr
         wc(iw) = ( wl(iw) + wu(iw) )/2.
      ENDDO
      wl(nw) = wu(nw-1)
      GO TO 9

*_______________________________________________________________________

 2    CONTINUE

* Input from table.  In this example:
* Wavelength grid will be read from a file.
* First line of table is:  nw = number of wavelengths (no. of intervals + 1)
* Then, nw wavelengths are read in, and assigned to wl(iw)
* Finally, wu(iw) and wc(iw) are computed from wl(iw)

c      wlabel = 'isaksen.grid'
      wlabel = 'combined.grid'

      fi = 'DATAE1/GRIDS/'//wlabel
      OPEN(unit=kin,file=fi,status='old')
      READ(kin,*) nw
      DO iw = 1, nw
         READ(kin,*) wl(iw)
      ENDDO
      CLOSE(kin)
      DO iw = 1, nw-1
         wu(iw) = wl(iw+1)
         wc(iw) = 0.5*(wl(iw) + wu(iw))
      ENDDO
      GO TO 9

*_______________________________________________________________________

 3    CONTINUE

* user-defined grid.  In this example, a single calculation is used to 
* obtain results for two 1 nm wide intervals centered at 310 and 400 nm:
* interval 1 : 1 nm wide, centered at 310 nm
* interval 3 : 2 nm wide, centered at 400 nm
* (inteval 2 : 310.5 - 399.5 nm, required to connect intervals 1 & 3)

      nw = 4
      wl(1) = 309.5
      wl(2) = 310.5
      wl(3) = 399.5
      wl(4) = 400.5
      DO iw = 1, nw-1
         wu(iw) = wl(iw+1)
         wc(iw) = 0.5*(wl(iw) + wu(iw))
      ENDDO
      GO TO 9

*_______________________________________________________________________

 4    CONTINUE
      wlabel = 'fast-TUV tropospheric grid'
      
      fi = 'DATAE1/GRIDS/fast_tuv.grid'
      OPEN(UNIT=kin,FILE=fi,STATUS='old')
      DO iw = 1, 4
         READ(kin,*)
      ENDDO

* skip wavelength shorter than 289.9 nm

      DO iw = 1, 10
         READ(kin,*)
      ENDDO
      nw = 8
      DO iw = 1, nw-1
         READ(kin,*) dum, wl(iw), dum, dum
      ENDDO
      wl(nw) = dum
      DO iw = 1, nw-1
         wu(iw) = wl(iw+1)
         wc(iw) = 0.5*(wl(iw) + wu(iw))
      ENDDO

      GO TO 9

*_______________________________________________________________________

 5    continue

* use standard grid up to 205.8 nm
* elsewhere, use 10 cm-1 grid to 1000 nm

      wlabel = 'combined.grid'
      fi = 'DATAE1/GRIDS/'//wlabel
      OPEN(unit=kin,file=fi,status='old')
      READ(kin,*) 
      DO iw = 1, 38
         READ(kin,*) wl(iw)
      ENDDO
      CLOSE(kin)

      DO i = 1, 3859
         iw = 3859 - i + 39
         wn(iw) = 10000 + 10*(i-1)
         wl(iw) = 1.E7/float(wn(iw))
      ENDDO

      nw = 3859 + 38
      nwint = nw - 1

      DO iw = 1, nwint
         wu(iw) = wl(iw+1)
         wc(iw) = (wl(iw) + wu(iw))/2.
      ENDDO

      GO TO 9

*_______________________________________________________________________

 6    CONTINUE

***** Correction for air-vacuum wavelength shift:
* The TUV code assumes that all working wavelengths are strictly IN-VACUUM. This is for ALL
* spectral data including extraterrestrial fluxes, ozone (and other) absorption cross sections,
* and various weighting functons (action spectra, photolysis cross sections, instrument spectral
* response functions).  If the original data are specified in-air, conversion to in-vacuum must be
* made when reading those data.

*  Occasionally, users may want their results to be given for wavelengths measured IN-AIR.
*   The shift between IN-VACUUM and IN-AIR wavelengths depends on the index of refraction
*   of air, which in turn depends on the local density of air, which in turn depends on
*   altitude, temperature, etc.
*  Here, we provide users with the option to use a wavelength grid IN-AIR, at the air density
*   corresponding to the selected altitude, airout.
*   The actual radiative transfer calculations will be done strictly with IN-VACUUM values.  

* create grid that will be nicely spaced in air wavelengths.

      wlabel = 'grid in air wavelengths'
      nw = nwint + 1
      wincr = (wstop - wstart) / FLOAT (nwint)
      DO iw = 1, nw-1
         wl(iw) = wstart + wincr*FLOAT(iw-1)
         wu(iw) = wl(iw) + wincr
         wc(iw) = ( wl(iw) + wu(iw) )/2.
      ENDDO
      wl(nw) = wu(nw-1)

* shift by refractive index to vacuum wavelengths, for use in tuv

      airout = 2.45e19
      mrefr = 1
      CALL wshift(mrefr, nw,    wl, airout)
      CALL wshift(mrefr, nwint, wc, airout)
      CALL wshift(mrefr, nwint, wu, airout)

c      do iw = 1, nw-1
c         write(33,333) wl(iw), wc(iw), wu(iw)
c      enddo
c 333  format(3(0pf11.4))

      GO TO 9

*_______________________________________________________________________
* Landgraf and Crutzen 1998
 10   CONTINUE
      nw = 6
      wl(1) = 289.0
      wl(2) = 305.5
      wl(3) = 313.5
      wl(4) = 337.5
      wl(5) = 422.5
      wl(6) = 752.5
      DO iw = 1, nw-1
         wu(iw) = wl(iw+1)
         wc(iw) = 0.5*(wl(iw) + wu(iw))
      ENDDO
      GO TO 9
*_______________________________________________________________________
* Wild 2000
 11   CONTINUE
      nw = 8
      wl(1) = 289.00
      wl(2) = 298.25
      wl(3) = 307.45
      wl(4) = 312.45
      wl(5) = 320.30
      wl(6) = 345.0
      wl(7) = 412.5
      wl(8) = 850.0

      DO iw = 1, nw-1
         wu(iw) = wl(iw+1)
         wc(iw) = 0.5*(wl(iw) + wu(iw))
      ENDDO
      GO TO 9
*_______________________________________________________________________
* Bian and Prather 2002
 12   CONTINUE
      nw = 8
      wl(1) = 291.0
      wl(2) = 298.3
      wl(3) = 307.5
      wl(4) = 312.5
      wl(5) = 320.3
      wl(6) = 345.0
      wl(7) = 412.5
      wl(8) = 850.0

      DO iw = 1, nw-1
         wu(iw) = wl(iw+1)
         wc(iw) = 0.5*(wl(iw) + wu(iw))
      ENDDO
      GO TO 9
*_______________________________________________________________________

*_______________________________________________________________________
* UV-b, UV-A, Vis

 13   CONTINUE
      nw = 4
      wl(1) = 280.0
      wl(2) = 315.0
      wl(3) = 400.0
      wl(4) = 700.0

      DO iw = 1, nw-1
         wu(iw) = wl(iw+1)
         wc(iw) = 0.5*(wl(iw) + wu(iw))
      ENDDO
      GO TO 9
*_______________________________________________________________________

 9    CONTINUE

* check grid for assorted improprieties:

      CALL gridck(kw,nw,wl,ok)

      IF (.NOT. ok) THEN
         WRITE(*,*)'STOP in GRIDW:  The w-grid does not make sense'
         STOP
      ENDIF

*_______________________________________________________________________

      RETURN
      END

*=============================================================================*

      SUBROUTINE gridz(zstart, zstop, nz, z, zout, izout)

*-----------------------------------------------------------------------------*
*=  PURPOSE:                                                                 =*
*=  Create the altitude grid for all interpolations and radiative transfer   =*
*=  calculations.  Grid may be irregularly spaced.  All altitudes are in     =*
*=  kilometers (km).  The altitude at index 1 specifies the surface elevation=*
*-----------------------------------------------------------------------------*
*=  PARAMETERS:                                                              =*
*=  nz  - INTEGER, number of altitude points (levels)                     (O)=*
*=  z   - REAL, vector of altitude levels (in km)                         (O)=*
*-----------------------------------------------------------------------------*

      IMPLICIT NONE

      INCLUDE 'params'

* input surface elevation, altitude for output

      REAL zstart, zstop, zout

* output: altitude working grid, index for output

      REAL z(kz)
      INTEGER nz
      INTEGER izout

* local:

      REAL zincr
      INTEGER i, n, nlev
      LOGICAL ok
*_______________________________________________________________________

* Set vertical grid of the atmosphere: atmospheric level altitudes 
* (in real km), including top-most level.
* User specifies grid (surface at lowest km value), increasing
* upwards:
*     -  nz = total number of user levels
*     -  z(I) = altitude in km for each level.
* z(1) is the elevation of the surface (km asl), and can be specified either
* here or in the main program.
* Non-uniform spacing is possible:
*     - zincr = altitude increment between current and previous level (km)
*     - nlev = number of levels in current equally-spaced section
*     - n = index of top level of equally-spaced section
* Note "levels" are vertical points
*      "layers" are vertical distances between levels


* Grid selection options:
* 1 = standard equally spaced grid, manual
* 2 = standard equally spaced grid, auto generated
* 3 = variable spacing grid, example for snow
* 4 = mirage z-grid for Mexico City
* 5 = arbitrary user-defined grid

      GO TO 1

*-----grid option 2: manual -----------------
* entire grid (nz levels) in increments zincr 

 1    CONTINUE
      WRITE(*,*) 'equally spaced z-grid'
      zincr = (zstop - zstart) / FLOAT(nz - 1)
      z(1) = zstart
      DO i = 2, nz
         z(i) = z(1) + zincr*FLOAT(i-1)
      ENDDO
      GOTO 10

*-----grid option 3: automatic -----------------
* entire grid (nz levels) in increments zincr 

 2    CONTINUE
      WRITE(*,*) 'equally spaced z-grid'
      zincr = (zstop - zstart) / FLOAT(nz - 1)
      nlev = nz-1
      n = 1
      CALL buildz(zincr, nlev, n, z)
      GOTO 10

*-----grid option 4: variable grid example----------------------
*-----copy & edit this section for non-uniform grid----
* the example provided below is high vertical resolution in 
*   snow, with atmosphere above it.

 3    CONTINUE
      WRITE(*,*) 'snow-atmosphere grid'
* 0.-10. cm from ground, in 1 cm increments ( 1 cm = 1e-5 km):
      zincr = 1.e-5
      nlev = 10
      n = 1
      CALL buildz(zincr,nlev,n,z)

* 10-90 cm from ground, in 10 cm increments ( 1 cm = 1e-5 km):
      zincr = 1.e-4
      nlev = 8
      CALL buildz(zincr,nlev,n,z)

* 90-95 cm from ground, in 1x 5 cm increment ( 1 cm = 1e-5 km):
      zincr = 5.e-5
      nlev = 1
      CALL buildz(zincr,nlev,n,z)

* 95-99 cm from ground, in 4x 1 cm increments ( 1 cm = 1e-5 km):
      zincr = 1.e-5
      nlev = 4
      CALL buildz(zincr,nlev,n,z)

* 99-99.5 cm from ground, in 1x 0.5 cm increment ( 1 cm = 1e-5 km):
      zincr = 5.e-6
      nlev = 1
      CALL buildz(zincr,nlev,n,z)

* 99.5 centimeters - 1m, in 0.1 cm increments (1 cm = 1e-5 km):
      zincr = 1.e-6
      nlev = 5
      CALL buildz(zincr,nlev,n,z)

*atmosphere
* 1.-10. m in 1 m increments
      zincr = 1.e-3
      nlev = 9
      CALL buildz(zincr,nlev,n,z)

* 10.-100 m in 10 m increments
      zincr = 1.e-2
      nlev = 9
      CALL buildz(zincr,nlev,n,z)

* 100.- 1000. meters, in 100 m increments
      zincr = 1.e-1
      nlev = 9
      CALL buildz(zincr,nlev,n,z)

* 1.-2. km in 1 km increments
      zincr = 1.
      nlev =  1
      CALL buildz(zincr,nlev,n,z)

* 2.-80. km in 2 km increments
      zincr = 2.
      nlev =  39
      CALL buildz(zincr,nlev,n,z)

      GOTO 10

*-----grid option 4:  grid for Mexico City

 4    CONTINUE
      WRITE(*,*) 'mirage z-grid'

* grid for mirage km: incr(range)i 
* 0.1(0-4)   2-41
* 0.2(4-8)   42-61
*   1(8-30)  62-83
*   2(30-50) 84-93 
*   5(50-80) 94-99

      nz = 99
      z(1) = zstart
      DO i = 2, 41
         z(i) = z(1) + 0.1*FLOAT(i-1)
      ENDDO
      DO i = 42, 61
         z(i) = z(41) + 0.2*FLOAT(i-41)
      ENDDO
      DO i = 62, 83
         z(i) = z(61) + 1.*FLOAT(i-61)
      ENDDO
      DO i = 84, 93
         z(i) = z(83) + 2.*FLOAT(i-83)
      ENDDO
      DO i = 94, 99
         z(i) = z(93) + 5.*FLOAT(i-93)
      ENDDO
      GOTO 10
		
*-----grid option 5: user defined

 5    CONTINUE

* insert your grid values here:
* specify:
*  nz = total number of altutudes
* Table:  z(iz), where iz goes from 1 to nz
* trivial example of 2-layer (3-altitudes) shown below, user should modify
      WRITE(*,*) 'user-defined grid, named...'
      nz = 3
      z(1) = 0.
      z(2) = 10.
      z(3) = 80.

      GOTO 10

*-----end of user options.
*-----grid option 6: high resolution window

6     CONTINUE

* insert your grid values here:
* specify:
*  nz = total number of altutudes
* Table:  z(iz), where iz goes from 1 to nz

      WRITE(*,*) 'user-defined grid, named...'
      goto 10

*-----end of user options.
*------------------------------------------------

 10   CONTINUE      

* Insert additional altitude for selected outputs.

      DO i = 1, nz
         IF(ABS(z(i) - zout) .LT. 0.001) THEN
            izout = i
            GO TO 24
         ENDIF
      ENDDO

* locate index for new altitude

      izout = 0
      DO i = 1, nz
         IF(z(i) .GT. zout) THEN
            izout = i
            GO TO 22
         ENDIF
      ENDDO
 22   CONTINUE
      IF(izout .LE. 1) STOP 'zout not in range - '

* shift overlying levels and insert new point

      nz = nz + 1
      DO i = nz, izout + 1, -1
         z(i) = z(i-1)
      ENDDO
      z(izout) = zout

 24   CONTINUE

* check grid for assorted improprieties:

 99	CONTINUE
      CALL gridck(kz,nz,z,ok)

      IF (.NOT. ok) THEN
         WRITE(*,*)'STOP in GRIDZ:  The z-grid does not make sense'
         STOP
      ENDIF

*_______________________________________________________________________

      RETURN
      END

*=============================================================================*
      SUBROUTINE  buildz(zincr,nlev,n,z)
*-----------------------------------------------------------------------------*
*= Purpose: to construct the altitude grid from parameters in gridz          =*
*-----------------------------------------------------------------------------*

      IMPLICIT NONE
      INCLUDE 'params'

      REAL zincr                   ! i
      INTEGER nlev                 ! i
      INTEGER n                    ! i/o
      REAL z(kz)                   ! i/o
      INTEGER i, j                    ! internal

      j = 0
      DO i = n + 1, n + nlev
        j = j + 1
        z(i) = z(n) + FLOAT(j)*zincr
      ENDDO
      n = n + nlev

      RETURN
      END

*=============================================================================*

      SUBROUTINE gridt(lat, lon, tmzone,
     $     iyear, imonth, iday,
     $     lzenit, tstart, tstop,
     $     nt, t, sza, sznoon, esrm2)

*-----------------------------------------------------------------------------*
*=  Subroutine to create time (or solar zenith angle) grid                   =*
*=  Also computes earth-sun distance (1/R**2) correction.                    =*
*-----------------------------------------------------------------------------*

      IMPLICIT NONE

      INCLUDE 'params'

* INPUTS

      REAL lat, lon, tmzone
      INTEGER iyear, imonth, iday
      LOGICAL lzenit
      INTEGER nt
      REAL tstart, tstop


* OUTPUTS

      REAL t(kt), sza(kt), esrm2(kt), sznoon

* INTERNAL

      INTEGER it
      REAL ut, dt

      INTEGER jday, nday
      LOGICAL oky, okm, okd

      REAL az, el, elnoon, soldia, soldst

*  switch for refraction correction to solar zenith angle. Because
* this is only for the observed sza at the surface, do not use.

      LOGICAL lrefr
      DATA lrefr /.FALSE./

***************

      IF(nt .EQ. 1) THEN
         dt = 0.
      ELSE
         dt = (tstop - tstart) / FLOAT(nt - 1)
      ENDIF

      DO 10 it = 1, nt
         t(it) = tstart + dt * FLOAT(it - 1)

* solar zenith angle calculation:
*  If lzenit = .TRUE., use selected solar zenith angles, also
*  set Earth-Sun distance to 1 AU.

         IF (lzenit) THEN
            sza(it) = t(it)
            esrm2(it) = 1.

*  If lzenit = .FALSE., compute solar zenith angle for specified
* location, date, time of day.  Assume no refraction (lrefr = .FALSE.)
*  Also calculate corresponding
* Earth-Sun correcton factor. 

         ELSE
            CALL calend(iyear, imonth, iday,
     $           jday, nday, oky, okm, okd)
            IF( oky .AND. okm .AND. okd) THEN

               ut = t(it) - tmzone
               CALL sunae(iyear, jday, ut, lat, lon, lrefr,
     &              elnoon, az, el, soldia, soldst )
               sza(it) = 90. - el
               sznoon = 90. - elnoon
               esrm2(it) = 1./(soldst*soldst)
            ELSE
               WRITE(*,*) '**** incorrect date specification'
               STOP ' in gridt '
            ENDIF
         ENDIF
            
 10   CONTINUE
      RETURN
      END

*=============================================================================*

      SUBROUTINE gridck(k,n,x,ok)

*-----------------------------------------------------------------------------*
*=  PURPOSE:                                                                 =*
*=  Check a grid X for various improperties.  The values in X have to comply =*
*=  with the following rules:                                                =*
*=  1) Number of actual points cannot exceed declared length of X            =*
*=  2) Number of actual points has to be greater than or equal to 2          =*
*=  3) X-values must be non-negative                                         =*
*=  4) X-values must be unique                                               =*
*=  5) X-values must be in ascending order                                   =*
*-----------------------------------------------------------------------------*
*=  PARAMETERS:                                                              =*
*=  K  - INTEGER, length of X as declared in the calling program          (I)=*
*=  N  - INTEGER, number of actual points in X                            (I)=*
*=  X  - REAL, vector (grid) to be checked                                (I)=*
*=  OK - LOGICAL, .TRUE. -> X agrees with rules 1)-5)                     (O)=*
*=                .FALSE.-> X violates at least one of 1)-5)                 =*
*-----------------------------------------------------------------------------*

      IMPLICIT NONE

      INCLUDE 'params'

* input:
      INTEGER k, n
      REAL x(k)

* output:
      LOGICAL ok

* local:
      INTEGER i
*_______________________________________________________________________

      ok = .TRUE.

* check if dimension meaningful and within bounds

      IF (n .GT. k) THEN
         ok = .false.
         WRITE(kout,100)
         RETURN
      ENDIF         
  100 FORMAT('Number of data exceeds dimension')

      IF (n .LT. 2) THEN
         ok = .FALSE.
         WRITE(kout,101)
         RETURN
      ENDIF
  101 FORMAT('Too few data, number of data points must be >= 2')

* disallow negative grid values

      IF(x(1) .LT. 0.) THEN
         ok = .FALSE.
         WRITE(kout,105)
         RETURN
      ENDIF
  105 FORMAT('Grid cannot start below zero')

* check sorting

      DO 10, i = 2, n
         IF( x(i) .LE. x(i-1)) THEN
            ok = .FALSE.
            WRITE(kout,110)
            RETURN
         ENDIF
   10 CONTINUE
  110 FORMAT('Grid is not sorted or contains multiple values')
*_______________________________________________________________________

      RETURN
      END
