PROGRAM tuv !_____________________________________________________________________________! ! Tropospheric Ultraviolet-Visible (TUV) radiation model ! ! Version 4.2 - beta ! ! March 2002 by Madronich et al. ! ! DEK Version 1.0 ! !_____________________________________________________________________________! ! This program is free software; you can redistribute it and/or modify ! ! it under the terms of the GNU General Public License as published by the ! ! Free Software Foundation; either version 2 of the license, or (at your ! ! option) any later version. ! ! The TUV package is distributed in the hope that it will be useful, but ! ! WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTIBI- ! ! LITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public ! ! License for more details. ! ! To obtain a copy of the GNU General Public License, write to: ! ! Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. ! !-----------------------------------------------------------------------------! ! To contact the authors, please mail to: ! ! Sasha Madronich, NCAR/ACD, P.O.Box 3000, Boulder, CO, 80307-3000, USA or ! ! send email to: sasha@ucar.edu or tuv@acd.ucar.edu ! !-----------------------------------------------------------------------------! ! Copyright (C) 1994,95,96,97,98,99,2000,01,02 University Corporation for ! ! Atmospheric Research ! !-----------------------------------------------------------------------------! implicit none ! ...include parameter file include 'params' !============================================================================== ! SECTION 1: VARIABLES AND PARAMETERS !============================================================================== ! ... wavelength grid: integer nw, iw real wl(kw), wc(kw), wu(kw) ! ... altitude grid integer nz, iz real z(kz) ! ... solar zenith angle and azimuth ! ... slant pathlengths in spherical geometry real zen, azim integer nid(0:kz) real dsdh(0:kz,kz) ! ... extra terrestrial solar flux and earth-sun distance ^-2 real etflux(kw), etf(kw) real esfact ! ... ozone absorption cross section real o3xs(kw), s226(kw), s263(kw), s298(kw) ! ... o2 absorption cross section real o2xs(kz,kw), o2xs1(kw) ! ... so2 absorption cross section real so2xs(kw) ! ... no2 absorption cross section real no2xs(kw) ! ... atmospheric optical parameters real tlev(kz), tlay(kz) real airden(kz), dcol(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) ! ... spectral irradiance and actinic flux (scalar irradiance) real edir(kz), edn(kz), eup(kz) real sirrad(kz,kw) real fdir(kz), fdn(kz), fup(kz) real saflux(kz,kw) ! ... spectral weighting functions and weighted radiation integer ns, is real sw(ks,kw), rate(ks,kz), dose(ks) real drdw character*40 label(ks) ! ...photolysis coefficients (j-values) integer nj, ij real sj(kj,kz,kw), valj(kj,kz) real djdw character*40 jlabel(kj) ! ... new sea level pressure, surface dobson, etc. real pmbnew, dobnew, so2new, no2new, visnew, cldnew, albnew ! ... Used in set gridz integer izout real z1, zout, zaird, ztemp ! ... date for etf integer idate ! ... other user-defined variables here: real n2den(kz) ! N2 Concentration, molecules cm-3 real noden(kz) ! NO Concentration, molecules cm-3 real o3den(kz) ! O3 Concentration, molecules cm-3 real o2den(kz) ! O2 Concentration, molecules cm-3 real o2vmr(kz) ! O2 Volume Mixing Ratio, Derived in seto2 real pmlev(kz) ! Pressure (hPa) real dobson ! Dobson Unit variable, returned from seto3 real o3vcol(kz) ! O3 Vertical Column, molecules cm-2 real o3scol(kz) ! O3 Slant Column, molecules cm-2 real o2scol(kz) ! O2 Slant Column, molecules cm-2 real noscol(kz) ! NO Slant Column, molecules cm-2 real etfphot(kw) ! Extra Terr. Flux, units quanta cm-2 sec-1 nm-1 real jno(kz) ! Total JNO, units sec-1, Minsch.+Siskind, 1993 real fnorm(kz) ! fdir(iz) + fdn(iz) + fup(iz) ! ... Look Up Table Arrays real dobson_ref integer numsza, numcolo3, numz, numwc, numalb integer icolo3, isza, ialb parameter (numalb = 6) parameter (numsza = 24) parameter (numcolo3 = 20) parameter (numz = 151) parameter (numwc = 101) real sza(numsza) real alb(numalb) real colo3_ref(numz) real o3col(numz) real colo3_fact(numcolo3) real wlintv(kw) real rad_src_funct(numalb, numcolo3, numsza, numz, numwc) data sza /0.0, 15.0, 30.0, 45.0, 60.0, $ 70.0, 75.0, 80.0, 82.5, 85.0, $ 87.0, 89.0, 90.0, 90.5, 91.0, $ 91.5, 92.0, 92.5, 93.0, 93.5, $ 94.0, 95.0, 96.0, 97.0/ data alb /0.05, 0.1, 0.2, 0.5, 0.75, 1.0/ ! This will cover surface values from 30 to 600 DU data colo3_fact /0.1, 0.2, 0.3, 0.4, 0.5, $ 0.6, 0.7, 0.8, 0.9, 1.0, $ 1.1, 1.2, 1.3, 1.4, 1.5, $ 1.6, 1.7, 1.8, 1.9, 2.0/ real AVO real GASC real PZ character*60 dir character*80 pathname, pathname1 dir = '../SRC/INPUT/' !============================================================================== ! ... SET GRIDS !============================================================================== ! ------------------------------------------------- ! ... Set altitudes (creates altitude grid, ! locates index for selected output, izout) ! ------------------------------------------------- ! ... z1 = surface elevation (km) z(1) = z1 z1 = 0. call gridz(z1,zout,izout,nz,z) ! ------------------------------------------------- ! ... Set wavelengths (creates wavelength ! grid: lower, center, upper of each bin) ! ------------------------------------------------- pathname = TRIM(dir)//'GRIDS/waccm2_ref_101.grid' call gridw(nw,wl,wc,wu,pathname) !============================================================================== ! ... SPECTRAL DATA !============================================================================== ! ------------------------------------------------ ! ... Read extra terrestrial flux data: ! ------------------------------------------------ pathname = TRIM(dir)//'SOLAR_FLUX/' call rdetfl(nw,wl,etflux,pathname) ! ------------------------------------------------ ! ... Read cross section data for ! o2 (will overwrite at lyman-alpha ! and srb wavelengths see subroutine ! la_srb.f) ! ! Read cross sections for SO2 ! Read Cross sections for NO2 n ! ------------------------------------------------ pathname = TRIM(dir)//'XSQY/' call rdo2xs (nw,wl,wc, o2xs1,pathname) call rdso2xs(nw,wl,so2xs,pathname) call rdno2xs(nw,wl,wc, nz, no2xs,pathname) !============================================================================== ! ... SET MODEL ATMOSPHERE !============================================================================== ! ----------------------------------------------- ! ... set temperature profile ! ----------------------------------------------- pathname = TRIM(dir)//'REF_ATMS/' call settmp(nz,z,tlev,tlay,pathname) ! ----------------------------------------------- ! ... set air profile and Rayleigh ! optical depths ! ----------------------------------------------- pmbnew = -999.0 pathname = TRIM(dir)//'REF_ATMS/' call setair(pmbnew,nz,z,nw,wl,airden,dtrl,dcol,pathname) ! ----------------------------------------------- ! ... set N2 conc (molec. cm-3) ! ----------------------------------------------- n2den(:) = airden(:)*0.790 ! ----------------------------------------------- ! ... set NO profile (molec cm-3) ! ----------------------------------------------- pathname = TRIM(dir)//'REF_ATMS/' CALL setno(nz,z,noden,pathname) ! ----------------------------------------------- ! ... set Pressure profile (hPa) ! Pm is input is consistent with ! T and AD ! ----------------------------------------------- ! CALL setpres(nz,z,pmlev,pathname) AVO = 6.022098e+23 GASC = 8.20580e-02 PZ = 1013.26 do iz = 1,nz pmlev(iz) = PZ*airden(iz)*(tlev(iz)*Gasc*1000.)/AVO enddo ! ----------------------------------------------- ! ... Set O2 profile (molec. cm-3) ! and derive dto2 ! ----------------------------------------------- ! ... Set O2 vmr; ! optical depths in lyman-alpha and ! srb will be over-written ! in subroutine la_srb.f ! ----------------------------------------------- pathname = TRIM(dir)//'REF_ATMS/' call seto2(nz,z,nw,wl,airden,dcol, $ o2xs1,dto2,o2den,o2vmr,pathname) ! ----------------------------------------------- ! ... Set O3 reference atmosphere, ! ozone vertical profile and optical ! depths ! ----------------------------------------------- dobnew = 300.0 pathname = TRIM(dir)//'REF_ATMS/' pathname1 = TRIM(dir)//'XSQY/' call seto3(dobnew,nz,z,nw,wl, $ tlay, o3xs, $ dto3,o3den,o3vcol,dobson, $ pathname,pathname1) colo3_ref(:) = o3vcol (:) dobson_ref = dobson ! ----------------------------------------------- ! ... SO2 vertical profile and optical ! depths ! ----------------------------------------------- so2new = 0.0 pathname = TRIM(dir)//'REF_ATMS/' call setso2(so2new,nz,z,nw,wl,so2xs, $ tlay,dtso2,pathname) ! ----------------------------------------------- ! ... NO2 vertical optical depth (also ! has temperature correction) ! ----------------------------------------------- no2new = 0.0 pathname = TRIM(dir)//'REF_ATMS/' call setno2(no2new,nz,z,nw,wl,no2xs, $ tlay,dtno2,pathname) ! ----------------------------------------------- ! ... Cloud vertical profile, optical ! depths, single scattering albedo, ! asymmetry factor ! ----------------------------------------------- cldnew = 0.0 pathname = TRIM(dir)//'REF_ATMS/' call setcld(cldnew,nz,z,nw,wl, $ dtcld,omcld,gcld,pathname) ! ----------------------------------------------- ! ... Aerosol vertical profile, optical ! depths, single scattering albedo, ! asymmetry factor ! ----------------------------------------------- visnew = -999. pathname = TRIM(dir)//'REF_ATMS/' call setaer(visnew,nz,z,nw,wl, $ dtaer,omaer,gaer,pathname) ! ----------------------------------------------- ! ... set albedo ! ----------------------------------------------- albnew = 0.2 call setalb(albnew,nw,wl, $ albedo) !============================================================================== ! ... TIME AND LOCATION !============================================================================== ! ----------------------------------------------- ! ... Specify date and compute ! earth-sun distance correction ! then multiply extra-terrestrial ! spectral irradiance by this factor. ! ----------------------------------------------- idate = 940321 call sundis(idate, esfact) esfact = 1.0 write(kout,*) 'idate = ', idate,' esfact = ', esfact ! ----------------------------------------------- ! ... Convert to quanta s-1 nm-1 cm-2 ! (1.e-4 * (wc*1e-9) / ! (hc = 6.62E-34 * 2.998E8) ) ! ----------------------------------------------- do iw = 1, nw-1 etf(iw) = etflux(iw) * esfact etfphot(iw) = etf(iw)*5.039e11 * wc(iw) enddo !============================================================================== ! ... LOOKUP TABLE LOOPS START HERE... !============================================================================== ! ... rad_src_funct(numalb, numcolo3, numsza, numz, numwc) ! ... rad_src_funct(ialb, icolo3, isza, iz, iw) do ialb = 1,numalb ! numalb = 6 do icolo3 = 1,numcolo3 ! numcolo3 = 20 do isza = 1,numsza ! numsza = 24 write(*,777) "ialb=", ialb, "icolo3=", icolo3, "isza =",isza 777 format(a6,i2,3x,a7,i2,3x,a6,i2,3x,a6,i2) ! ----------------------------------------------- ! ... Set Surface Albedo: ! ----------------------------------------------- do iw = 1, nw-1 albedo(iw) = alb(ialb) enddo ! ----------------------------------------------- ! ... Set solar zenith angle calculation: ! ----------------------------------------------- zen = sza(isza) ! ----------------------------------------------- ! ...Ozone reference atmosphere, ! calculate dobson unit ! ----------------------------------------------- dobnew = dobson_ref * colo3_fact(icolo3) ! ----------------------------------------------- ! ... Set O3 reference atmosphere, ! ozone vertical profile and optical ! depths ! ----------------------------------------------- ! ... dto3 is the only influence on RSF. ! o3vcol will change with dobnew ! o3den is ref atm. This is only used ! for scaling. ! ----------------------------------------------- pathname = TRIM(dir)//'REF_ATMS/' pathname1 = TRIM(dir)//'XSQY/' call seto3(dobnew,nz,z,nw,wl, $ tlay, o3xs, $ dto3,o3den,o3vcol,dobson, $ pathname,pathname1) ! ---------------------------------------------- ! ... Slant path lengths for spherical ! geometry ! ---------------------------------------------- call sphers(nz, z, zen, dsdh, nid) ! ---------------------------------------------- ! ... Derive slant and vertical columns ! for air and O2 ! ---------------------------------------------- call airmas(nz, z, zen, dsdh, nid, dcol, vcol, scol) call o2mas (nz, z, dsdh, nid, o2den, o2scol) write(*,677) "alb", "dobn","zen","o3den1", "o3vcol1", $ "o2scol1", "o2scol20","o2scol40","o2scol60", $ "o2scol80","o2scol100" write(*,678) albedo(1), dobnew, zen, o3den(1), o3vcol(1), & o2scol(1), o2scol(20),o2scol(40),o2scol(60), & o2scol(80),o2scol(100) ! ---------------------------------------------- ! ... Recalculate effective o2 optical ! depth and cross sections for ! lyman-alpha and SRBs, must know ! zenith angle ! ---------------------------------------------- pathname = TRIM(dir)//'XSQY/' call la_srb(nz,z,tlev,nw,wl,o2scol,vcol,scol,o2xs1, $ dto2,o2xs,pathname) !============================================================================== ! ... Radiative Transfer (wavelength loop) !============================================================================== ! ---------------------------------------------- ! ... Main wavelength loop: ! ---------------------------------------------- do iw = 1, nw-1 wlintv(iw) = wu(iw) - wl(iw) ! ---------------------------------------------- ! ... Monochromatic radiative transfer: ! outputs are edir(iz), eup, fdir, ! fdn, fup ! ... dir = direct beam, ! dn = down-welling diffuse, ! up = up-welling diffuse ! ---------------------------------------------- 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) ! ---------------------------------------------- ! ... Calculate Wavelength Dependent ! radiative source function ! ---------------------------------------------- do iz = 1, nz saflux (iz,iw) = fdir(iz) + fdn(iz) + fup(iz) IF (saflux(iz,iw) .LT. 0.0) THEN saflux(iz,iw) = 0.0 ENDIF rad_src_funct(ialb, icolo3, isza, iz, iw) = saflux(iz,iw) end do ! end altitude loop end do !end wavelength loop end do !isza end do !icolo3 end do !ialb !============================================================================== ! ... OUTPUT !============================================================================== call NetCDF_out(wc, wl, wlintv, z, pmlev, alb, sza, colo3_fact, $ colo3_ref, rad_src_funct, etfphot) 678 format(f5.2,2x,f6.1,2x,f7.2,2x,1p,8e13.5) 677 format(1x,a3,4x,a4,3x,a6,2x,8(4xa9)) end program TUV