PROGRAM tuv,25
c*     Tropospheric Ultraviolet-Visible (TUV) radiation model
c*     Version 4.1
c*     July 2000 by Madronich et al.
c*= This program is free software;  you can redistribute it and/or modify     =*
c*= it under the terms of the GNU General Public License as published by the  =*
c*= Free Software Foundation;  either version 2 of the license, or (at your   =*
c*= option) any later version.                                                =*
c*= The TUV package is distributed in the hope that it will be useful, but    =*
c*= WITHOUT ANY WARRANTY;  without even the implied warranty of MERCHANTIBI-  =*
c*= LITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public     =*
c*= License for more details.                                                 =*
c*= To obtain a copy of the GNU General Public License, write to:             =*
c*= Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.   =*
c*= To contact the authors, please mail to:                                   =*
c*= Sasha Madronich, NCAR/ACD, P.O.Box 3000, Boulder, CO, 80307-3000, USA  or =*
c*= send email to:                                            =*
c*= Copyright (C) 1994,95,96,97,98,99  University Corporation for Atmospheric =*
c*= Research                                                                  =*


c* include parameter file

      INCLUDE 'params'

c* ___ SECTION 1: VARIABLES AND PARAMETERS ______________________________

c* wavelength grid:

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

c* altitude grid

      INTEGER nz, iz
      REAL z(kz)

c* solar zenith angle and azimuth
c* slant pathlengths in spherical geometry

      REAL zen, azim
      INTEGER nid(0:kz)
      REAL dsdh(0:kz,kz)

c* extra terrestrial solar flux and earth-Sun distance ^-2

      REAL f(kw), etf(kw)
      REAL esfact

c* ozone absorption cross section and ozone optical depth:

      REAL xso3(kw), s226(kw), s263(kw), s298(kw)

c* O2 absorption cross section

      REAL xso2(kz,kw)

c* SO2 absorption cross section

      REAL xsso2(kw)

c* NO2 absorption cross section

      REAL xsno2(kw)

c* atmospheric optical parameters:

      REAL tlev(kz), tlay(kz)
      REAL airlev(kz), colinc(kz), vcol(kz), scol(kz)
      REAL dtrl(kz,kw)
      REAL dto3(kz,kw), dto2(kz,kw), dtso2(kz,kw), dtno2(kz,kw)
      REAL dtcld(kz,kw), omcld(kz,kw), gcld(kz,kw)
      REAL dtaer(kz,kw), omaer(kz,kw), gaer(kz,kw)
      REAL albedo(kw)

c* spectral irradiance and actinic flux (scalar irradiance):

      REAL edir(kz), edn(kz), eup(kz)
      REAL fdir(kz), fdn(kz), fup(kz)

c* spectral weighting functions and weighted radiation:

      INTEGER ns, is
      REAL sw(ks,kw), rate(ks), dose(ks)
      REAL sirrad, sprate
      CHARACTER*40 label(ks)

c*! j-values:

      INTEGER nj, ij
      REAL sj(kj,kz,kw), valj(kj,kz)
      REAL saflux, deltaj
      CHARACTER*40 jlabel(kj)

c* new sea level pressure, surface dobson, etc.

      REAL pmbnew, dobnew, so2new, no2new

c* Location and time

      REAL alat, along
      INTEGER idate
      REAL dtime, ut, ut0

c* Commonly used looping indices

      integer i, j
      integer  idat, idob, itime, izen

c* Other user-defined variables here:

      integer nzen
      real  sav(100,kz), sza(100)


c* ___ SECTION 2: SET GRIDS _________________________________________________

c* wavelengths

      CALL gridw(nw,wl,wc,wu)

c* for version B only (no sr bands):  limit  to wavelengths longer than 250 nm.

C      IF (wl(iw) .LE. 250.) STOP

c* altitudes (set the surface elevation, z(1), in km asl).

      z(1) = 0.
      CALL gridz(nz,z)

c* ___ SECTION 3: SPECTRAL DATA ____________________________

c* read (and grid) extra terrestrial flux data:

      CALL rdetfl(nw,wl,f)

c* read cross section data for
c*    ozone (temperature-dependent)
c*    SO2
c*    NO2

      CALL rdo3xs(nw,wl,xso3,s226,s263,s298)
      CALL rdso2xs(nw,wl,xsso2)
      CALL rdno2xs(nw,wl,xsno2)

c* ___ SECTION 4: SET MODEL ATMOSPHERE __________________________________

c* temperature profile

      CALL settmp(nz,z,
     $     tlev,tlay)

c*  air profile and Rayleigh optical depths

      pmbnew = -999.
      CALL setair(pmbnew,
     $     nz,z,nw,wl,
     $     airlev,dtrl,colinc)

c* Photo-chemical and photo-biological weigting functions.
c* For pchem, need to know temperature and pressure profiles.
c* Output:
c* from pbiol:  s(ks,kw) - for each weigting function label(ks)
c* from pchem:  sj(kj,kz,kw) - for each reaction jlabel(kj)

      is = 0
      CALL pbiol1(nw,wl,wc,is,sw,label)
      ns = is

      CALL pchem(nw,wl,nz,tlev,airlev,
     $     nj,sj,jlabel)

c* ozone optical depths (must give temperature)

      dobnew = 300.
      CALL setozo(dobnew,
     $     nz,z,nw,wl,
     $     xso3,s226,s263,s298,tlay,
     $     dto3)

c* SO2 optical depth (also has temperature correction)
c* so2new = new column SO2, in Dobson Units

      so2new =  0.
      CALL setso2(so2new,
     $     nz,z,nw,wl,
     $     xsso2, tlay,
     $     dtso2)

c* NO2 optical depth (also has temperature correction)
c* no2new = new column NO2, in Dobson Units

      no2new =  0.
      CALL setno2(no2new,
     $     nz,z,nw,wl,
     $     xsno2, tlay,
     $     dtno2)

c*  cloud and aerosol optical depths:

      CALL setcld(nz,z,nw,wl,
     $     dtcld,omcld,gcld)

      CALL setaer(nz,z,nw,wl,
     $     dtaer,omaer,gaer)

c* surface albedo:

      CALL setalb(nw,wl,albedo)

c* ___ SECTION 5: TIME AND LOCATION _____________________________________

c* specify date and compute earth-sun distance correction

      idate = 940321
      CALL sundis(idate,esfact)
      WRITE(kout,*) 'idate = ', idate,' esfact = ', esfact
      do iw = 1, nw-1
         etf(iw) = f(iw) * esfact

c* specify latitude and longitude

      alat = 0.
      along = 0.
      WRITE(kout,*)'lat = ',alat,' long = ',along

c* below, can  chose between specific time (solar zenith angle is calculated)
c* or can set zenith angle to arbitrary value(s).  Loop DO 20 allows calculation
c* at multiple solar zenith angles  (or multiple times).

c* Set starting time (ut = Universal Time, hrs.), and
c* time increment (dtime, in seconds)

c      ut0 = 0.
c      dtime = 3600./4.

c* initalize time-integrated quantities

cC      call zero1(dose,ks)

c* Loop over time (alternatively, can loop over solar zenith angle)

c      DO 20, itime = 1, 96
C         ut = ut0 + (dtime/3600.) * FLOAT(itime-1)

         ut = 12.

c* solar zenith angle calculation:

       CALL zenith(alat,along,idate,ut,azim,zen)
       WRITE(kout,*) 'ut = ', ut, 'azimuth = ', azim, ' zen = ', zen


c* slant path lengths for spherical geometry

       CALL sphers(nz, z, zen, dsdh, nid)
       CALL airmas(nz, z, zen, dsdh, nid, colinc,
     $      vcol, scol)

c* effective O2 optical depth (SR bands, must know zenith angle!)
c* assign O2 cross section to sj(1,*,*)

       CALL seto2(nz,z,nw,wl,colinc,vcol,scol,dto2,xso2)
       CALL sjo2(nz,nw,xso2,1,sj)

c* ____ SECTION 7: WAVELENGTH LOOP ______________________________________

c* initialize for wavelength integration

       CALL zero1(rate,ks)
       CALL zero2(valj,kj,kz)

c** Main wavelength loop:

       DO 10, iw = 1, nw-1

c** monochromatic radiative transfer:
c*  outputs are  edir(iz), eup, fdir, fdn, fup

         CALL rtlink(nz,z,
     $        iw, albedo(iw), zen,
     $        dsdh,nid,
     $        dtrl,
     $        dto3,
     $        dto2,
     $        dtso2,
     $        dtno2,
     $        dtcld, omcld, gcld,
     $        dtaer,omaer,gaer,
     $        edir, edn, eup, fdir, fdn, fup)

c** surface irradiance and weighted radiation

         iz = 1
         sirrad = etf(iw) * (edir(iz) + edn(iz))

         DO 15, is = 1, ns
            sprate = sirrad * sw(is,iw)
            rate(is) = rate(is) + sprate * (wu(iw) - wl(iw))
 15      CONTINUE

c** spherical irradiance (actinic flux)
c* as a function of altitude
c* convert to quanta s-1 nm-1 cm-2
c* ( 1.e-4 * (wc*1e-9) / (hc = 6.62E-34 * 2.998E8) )

         DO 17 iz = 1, nz
            saflux = etf(iw)* 5.039e11 * wc(iw) *
     $           (fdir(iz) + fdn(iz) + fup(iz))

            DO 16, ij = 1, nj
               deltaj = saflux * sj(ij,iz,iw)
               valj(ij,iz) = valj(ij,iz) + deltaj * (wu(iw) - wl(iw))
 16         CONTINUE

 17      CONTINUE


c*^^^^^^^^^^^^^^^^ end wavelength loop


c** some examples of output:

c* dose rates weighted by specific action spectra:

      DO 35, is = 1, ns
         WRITE(kout,99) is, label(is), rate(is)

c* photolysis rate coefficients (j-values) at surface

      iz = 1
      DO 36, ij = 1, nj
         WRITE(kout,99) ij, jlabel(ij), valj(ij,iz)

 99   FORMAT(I4,1X,A40,1X,1PE10.3)

c*^^^^^^^^^^^^^^^^ end zenith loop