PROGRAM tuv *-----------------------------------------------------------------------------* *= Tropospheric Ultraviolet-Visible (TUV) radiation model =* *= Version 4.2 =* *= May 2003 =* *-----------------------------------------------------------------------------* *= Developed by Sasha Madronich with important contributions from: =* *= Chris Fischer, Siri Flocke, Julia Lee-Taylor, Bernhard Meyer, =* *= Irina Petropavlovskikh, Xuexi Tie, and Jun Zen. =* *= Special thanks to Knut Stamnes and co-workers for the development of the =* *= Discrete Ordinates code, and to Warren Wiscombe and co-workers for the =* *= development of the solar zenith angle subroutine. Citations for the many =* *= data bases (e.g. extraterrestrial irradiances, molecular spectra) may be =* *= found in the data files headers and/or in the subroutines that read them. =* *= To contact the author, write 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 =* *-----------------------------------------------------------------------------* *= 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. =* *-----------------------------------------------------------------------------* *= Copyright (C) 1994,95,96,97,98,99,2000,01,02,03 University Corporation =* *= for Atmospheric Research =* *-----------------------------------------------------------------------------* IMPLICIT NONE * Include parameter file INCLUDE 'params' * Wavelength grid: INTEGER nw, iw, nwint REAL wl(kw), wc(kw), wu(kw) REAL wstart, wstop * Altitude grid INTEGER nz, iz, izout REAL z(kz), zstart, zstop, zout * Solar zenith angle and azimuth * slant pathlengths in spherical geometry REAL sza(kt), zen INTEGER nid(0:kz) REAL dsdh(0:kz,kz) * Extra terrestrial solar flux and earth-Sun distance ^-2 REAL f(kw), etf(kw) REAL esfact(kt) * Ozone absorption cross section REAL o3xs(kz,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), cair(kz), vcol(kz), scol(kz) REAL dtrl(kz,kw) REAL co3(kz) 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*50 slabel(ks) * Photolysis coefficients (j-values) INTEGER nj, ij REAL sj(kj,kz,kw), valj(kj,kz) REAL djdw CHARACTER*50 jlabel(kj) *cba following added for MCM or saprc99 scheme INTEGER nj4,nj6,nj8 REAL sj4(kj,kz,kw),valj4(kj,kz) REAL sj6(kj,kz,kw), valj6(kj,kz),sj8(kj,kz,kw), valj8(kj,kz) CHARACTER*40 jlabel3(kj), jlabel4(kj), jlabel5(kj) CHARACTER*40 jlabel6(kj), jlabel7(kj),jlabel8(kj), jlabel9(kj) REAL resu4(kgout,kj), resu6(kgout,kj),resu8(kgout,kj) REAL toto1, toto2 **** Re-scaling factors (can be read from input file) * New surface albedo and surface pressure (milli bar) * Total columns of O3, SO2, NO2 (Dobson Units) * Cloud optical depth, altitude of base and top * Aerosol optical depth at 550 nm, single scattering albedo, Angstrom alpha REAL alsurf, psurf REAL o3col, so2col, no2col REAL taucld, zbase, ztop REAL tauaer, ssaaer, alpha * Location: Lat and Lon (deg.), surface elev (km) * Altitude, temperature and pressure for specific outputs REAL lat, lon REAL zaird, ztemp * Time and/or solar zenith angle INTEGER iyear, imonth, iday INTEGER it, nt REAL t(kt), tstart, tstop REAL tmzone, ut LOGICAL lzenit * number of radiation streams INTEGER nstr * input/output control LOGICAL intrct CHARACTER*6 inpfil, outfil INTEGER iout REAL dirsun, difdn, difup CHARACTER*1 again * Other user-defined variables here: * Save arrays for output: LOGICAL lirrad, laflux, lrates, ljvals, lmmech INTEGER isfix, ijfix, itfix, izfix, iwfix INTEGER nms, ims(ks), nmj, imj(kj) REAL svj_zj(kz,kj), svj_tj(kt,kj), svj_zt(kz,kt) REAL svr_zs(kz,ks), svr_ts(kt,ks), svr_zt(kz,kt) REAL svf_zw(kz,kw), svf_tw(kt,kw), svf_zt(kz,kt) REAL svi_zw(kz,kw), svi_tw(kt,kw), svi_zt(kz,kt) * Parameters for shifting wavelengths air <--> vacuum INTEGER mrefr LOGICAL lrefr REAL airout * some usual variables INTEGER i logical llamp * --- END OF DECLARATIONS --------------------------------------------- llamp = .TRUE. !llamp = .FALSE. * re-entry point 1000 CONTINUE * Open log file: c OPEN(UNIT=kout,FILE='tuvlog',STATUS='UNKNOWN') OPEN(UNIT=kout,FILE='tuvlog'//'.txt',STATUS='UNKNOWN') * ___ SECTION 1: SIMPLE INPUT VARIABLES -------------------------------- ******* Read simple input variables from a file: * can read interactively (intrct = .TRUE.) * or in batch mode (intrct = .FALSE.) intrct = .TRUE. c intrct = .FALSE. IF ( .NOT. intrct) inpfil = 'usrinp' CALL rdinp(intrct, $ inpfil, outfil, nstr, lat, lon, tmzone, $ iyear, imonth, iday, zstart, zstop, nz, $ wstart, wstop, nwint, tstart, tstop, nt, $ lzenit, alsurf, psurf, o3col, so2col, no2col, $ taucld, zbase, ztop, tauaer, ssaaer, alpha, $ dirsun, difdn, difup, zout, zaird, ztemp, $ lirrad, laflux, lmmech, lrates, isfix, nms, $ ljvals, ijfix, nmj, iwfix, itfix, izfix, $ ims, slabel, imj, jlabel) IF(outfil .EQ. 'screen') THEN iout = 6 ELSE iout = 30 ENDIF ************* Can overwrite basic inputs here manually: * Input and output files: * inpfil = input file name * outfil = output file name * Radiative transfer scheme: * nstr = number of streams * If nstr < 2, will use 2-stream Delta Eddington * If nstr > 1, will use nstr-stream discrete ordinates * Location (geographic): * lat = LATITUDE (degrees, North = positive) * lon = LONGITUDE (degrees, East = positive) * tmzone = Local time zone difference (hrs) from Universal Time (ut): * ut = timloc - tmzone * Date: * iyear = year (1950 to 2050) * imonth = month (1 to 12) * iday = day of month * Time of day grid: * tstart = starting time, local hours * tstop = stopping time, local hours * nt = number of time steps * lzenith = switch for solar zenith angle (sza) grid rather than time * grid. If lzenith = .TRUE. then * tstart = first sza in deg., * tstop = last sza in deg., * nt = number of sza steps. * esfact = 1. (Earth-sun distance = 1.000 AU) * Vertical grid: * zstart = surface elevation above sea level, km * zstop = top of the atmosphere (exospheric), km * nz = number of vertical levels, equally spaced * (nz will increase by +1 if zout does not match altitude grid) * Wavlength grid: * wstart = starting wavelength, nm * wstop = final wavelength, nm * nwint = number of wavelength intervals, equally spaced * if nwint < 0, the standard atmospheric wavelength grid, not * equally spaced, from 120 to 735 nm, will be used. In this * case, wstart and wstop values are ignored. * Surface condition: * alsurf = surface albedo, wavelength independent * psurf = surface pressure, mbar. Set to negative value to use * US Standard Atmosphere, 1976 (USSA76) * Column amounts of absorbers (in Dobson Units, from surface to space): * Vertical profile for O3 from USSA76. For SO2 and NO2, vertical * concentration profile is 2.69e10 molec cm-3 between 0 and * 1 km above sea level, very small residual (10/largest) above 1 km. * o3col = ozone (O3) * so2col = sulfur dioxide (SO2) * no2col = nitrogen dioxide (NO2) * Cloud, assumed horizontally uniform, total coverage, single scattering * albedo = 0.9999, asymmetry factor = 0.85, indep. of wavelength, * and also uniform vertically between zbase and ztop: * taucld = vertical optical depth, independent of wavelength * zbase = altitude of base, km above sea level * ztop = altitude of top, km above sea level * Aerosols, assumed vertical provile typical of continental regions from * Elterman (1968): * tauaer = aerosol vertical optical depth at 550 nm, from surface to space. * If negative, will default to Elterman's values (ca. 0.235 * at 550 nm). * ssaaer = single scattering albedo of aerosols, wavelength-independent. * alpha = Angstrom coefficient = exponent for wavelength dependence of * tauaer, so that tauaer1/tauaer2 = (w2/w1)**alpha. * Directional components of radiation, weighting factors: * dirsun = direct sun * difdn = down-welling diffuse * difup = up-welling diffuse * e.g. use: * dirsun = difdn = 1.0, difup = 0 for total down-welling irradiance * dirsun = difdn = difup = 1.0 for actinic flux from all directions * dirsun = difdn = 1.0, difup = -1 for net irradiance * Output altitude: * zout = altitude, km, for desired output. * If not within 1 m of altitude grid, an additional * level will be inserted and nz will be increased by +1. * zaird = air density (molec. cm-3) at zout. Set to negative value for * default USSA76 value interpolated to zout. * ztemp = air temperature (K) at zout. Set to negative value for * default USSA76 value interpolated to zout. * Output options, logical switches: * lirrad = output spectral irradiance * laflux = output spectral actinic flux * lmmech = output for NCAR Master Mechanism use * lrates = output dose rates (UVB, UVA, CIE/erythema, etc.) * Output options, integer selections: * isfix: if > 0, output dose rate for action spectrum is=isfix, tabulated * for different times and altitudes. * ijfix: if > 0, output j-values for reaction ij=ijfix, tabulated * for different times and altitudes. * iwfix: if > 0, output spectral irradiance and/or spectral actinic * flux at wavelength iw=iwfix, tabulated for different times * and altitudes. * itfix: if > 0, output spectral irradiance and/or spectral actinic * flux at time it=itfix, tabulated for different altitudes * and wavelengths. * izfix: if > 0, output spectral irradiance and/or spectral actinic * flux at altitude iz=izfix, tabulated for different times * and wavelengths. * nms: number of dose rates that will be reported. Selections must be * made interactively, or by editing input file. * nmj: number of j-values that will be reported. Selections must be * made interactively, or by editing input file. * The following default settings are also found in the input file 'defin1': c inpfil = defin1 c outfil = usrout c nstr = -2 c lat = 0. c lon = 0. c tmzone = 0. c iyear = 2002 c imonth = 3 c iday = 21 c zstart = 0. c zstop = 80. c nz = 81 c wstart = 280. c wstop = 420. c nwint = 140 c tstart = 12. c tstop = 20. c nt = 5 c lzenith = .FALSE. c alsurf = 0.1 c psurf = -999. c o3col = 300. c so2col = 0. c no2col = 0. c tcloud = 0. c zbase = 4. c ztop = 5. c tauaer = 0.235 c ssaaer = 0.99 c alpha = 1. c dirsun = 1. c difdn = 1. c difup = 0. c zout = 0. c zaird = -999. c ztemp = -999. c lirrad = .TRUE. c laflux = .FALSE. c lmmech = .FALSE. c lrates = .TRUE. c isfix = 0 * nms cannot be set here c ljvals = .FALSE. c ijfix = 0 * nmj cannot be set here c iwfix = 0 c itfix = 0 c izfix = 0 IF(nstr .LT. 2) THEN WRITE(kout,*) 'Delta-Eddington 2-stream radiative transfer' ELSE WRITE(kout,*) 'Discrete ordinates ', $ nstr, '-stream radiative transfer' ENDIF WRITE(*,*) 'calculating....' * ___ SECTION 2: SET GRIDS _________________________________________________ * wavelengths (creates wavelength grid: lower, center, upper of each bin) * NOTE: Wavelengths are in vacuum. To use wavelengths in air, see * Section 3 below, where you must set lrefr= .TRUE. CALL gridw(wstart, wstop, nwint, $ nw, wl, wc, wu) * altitudes (creates altitude grid, locates index for selected output, izout) CALL gridz(zstart, zstop, nz, z, zout, izout) if(izfix .gt. 0) izout = izfix * time/zenith (creates time/zenith angle grid, starting at tstart) CALL gridt(lat, lon, tmzone, $ iyear, imonth, iday, $ lzenit, tstart, tstop, $ nt, t, sza, esfact) *cba change the zenith angle : 90=86, 80=>78 DO i=1,nt IF (sza(i).eq.80.) sza(i)=78. IF (sza(i).eq.90.) sza(i)=86. IF (sza(i).eq.100.) sza(i)=90. ENDDO * ___ SECTION 3: SET UP VERTICAL PROFILES OF TEMPERATURE, AIR DENSITY, and OZONE______ ***** Temperature vertical profile, Kelvin * can overwrite temperature at altitude z(izout) CALL vptmp(nz,z, tlev,tlay) IF(ztemp .GT. nzero) tlev(izout) = ztemp ***** Air density (molec cm-3) vertical profile * can overwrite air density at altitude z(izout) CALL vpair(psurf, nz, z, $ airden, cair) IF(zaird .GT. nzero) airden(izout) = zaird ***** Correction for air-vacuum wavelength shift: * The TUV code assumes that all working wavelengths are strictly IN-VACUUM. This is assumed 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). * 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 output altitude, airout = airden(izout), by setting the logical variable * lrefr = .TRUE. (default is lrefr = .FALSE.). The wavelengths specified in gridw.f will be assumed * to be IN-AIR, and will be shifted here to IN-VACUUM values to carry out the calculatons. * The actual radiative transfer calculations will be done strictly with IN-VACUUM values. * If this shift is applied (i.e., if lrefr = .TRUE.), the wavelength grid will be shifted back to air * values just before the output is written. * Note: if this option is used (lref = .TRUE.), the wavelength values will be correct ONLY at the * selected altitude, iz = iout. The wavelength shift will be INCORRECT at all other altitudes. * Note: This option cannot be changed interactively in the input table. It must be changed here. C lrefr = .TRUE. lrefr = .FALSE. IF(lrefr) THEN airout = airden(izout) write(*,*) 'applying air to vacuum wavelength shift', $ izout, airout mrefr = 1 CALL wshift(mrefr, nw, wl, airout) CALL wshift(mrefr, nwint, wc, airout) CALL wshift(mrefr, nwint, wu, airout) ENDIF ***** Ozone vertical profile CALL vpo3(o3col, nz, z, cair, co3) * ___ SECTION 4: READ SPECTRAL DATA ____________________________ * read (and grid) extra terrestrial flux data: * OR read (and grid) lamp data: if(llamp) then CALL rdlamp(nw,wl,etf) else CALL rdetfl(nw,wl,etf) endif * read cross section data for * O2 (will overwrite at Lyman-alpha and SRB wavelengths * see subroutine la_srb.f) * O3 (temperature-dependent) * SO2 * NO2 CALL rdo2xs(nw,wl, o2xs1) CALL rdo3xs(nw,wl,nz,tlay, o3xs) CALL rdso2xs(nw,wl, so2xs) CALL rdno2xs(nw,wl, no2xs) ****** Spectral weighting functions * (Some of these depend on temperature T and pressure P, and therefore * on altitude z. Therefore they are computed only after the T and P profiles * are set above with subroutines settmp and setair.) * Photo-physical set in swphys.f (transmission functions) * Photo-biological set in swbiol.f (action spectra) * Photo-chemical set in swchem.f (cross sections x quantum yields) * Physical and biological weigthing functions are assumed to depend * only on wavelength. * Chemical weighting functions (product of cross-section x quantum yield) * for many photolysis reactions are known to depend on temperature * and/or pressure, and therefore are functions of wavelength and altitude. * Output: * from pphys & pbiol: s(ks,kw) - for each weighting function slabel(ks) * from pchem: sj(kj,kz,kw) - for each reaction jlabel(kj) * For pchem, need to know temperature and pressure profiles. CALL swphys(nw,wl,wc, ns,sw,slabel) CALL swbiol(nw,wl,wc, ns,sw,slabel) CALL swchem(nw,wl,nz,tlev,airden, nj,sj,jlabel) * cba pchem3,4,IN added for jvalues in SGMCM CALL pchem3(nw,wl,nz,tlev,airden, $ nj4,sj4,jlabel4,jlabel5) * sl pchem4 added for jvalues for MCM2001 CALL pchem4(nw,wl,nz,tlev,airden, $ nj6,sj6,jlabel6,jlabel7) * sl pchemIN added for jvalues for MCM2001 CALL pchemIN(nw,wl,nz,tlev,airden, $ nj8,sj8,jlabel8,jlabel9) * jlt pchem08 added for jvalues for difunctional nitrates CALL pchem08(nw,wl,nz,tlev,airden, $ nj6,sj6,jlabel6,jlabel7) c CALL swbiol2(nw,wl,wc, ns,sw,slabel) c CALL swbiol3(nw,wl,wc, ns,sw,slabel) * The following lines serve to print a list of the spectral weighting * functions. If new functions (e.g. action spectra, photo-reactions) * are added, this list should be used to replace the list in the * default input files (defin1, defin2, etc.). The true/false toggle * will be set to F, and should be changed manually to select weighting * functions for output. Note that if many more functions are added, it * may be necessary to increase the parameters ks and kj in the include * file 'params' * The program will stop after writing this list. * Comment out these lines when not generating a new list. ! OPEN(UNIT=50,FILE='spectra.list',STATUS='NEW') ! WRITE(50,500) ! 500 FORMAT(5('='),1X,'Available spectral weighting functions:') ! DO is = 1, ns ! WRITE(50,505) is, slabel(is) ! ENDDO ! WRITE(50,510) ! 510 FORMAT(5('='),1X,'Available photolysis reactions') ! DO ij = 1, nj ! WRITE(50,505) ij, jlabel(ij) ! ENDDO ! 505 FORMAT('F',I3,1X,A50) ! WRITE(50,520) ! 520 FORMAT(66('=')) ! CLOSE (50) ! STOP ****** * ___ SECTION 4: SET ATMOSPHERIC OPTICAL DEPTH INCREMENTS _____________________ * Rayleigh optical depth increments: CALL odrl(nz, z, nw, wl, cair, dtrl) * O2 vertical profile and O2 absorption optical depths * For now, O2 densitiy assumed as 20.95% of air density, can be change * in subroutine. * Optical depths in Lyman-alpha and SRB will be over-written * in subroutine la_srb.f CALL seto2(nz,z,nw,wl,cair,o2xs1, dto2) * Ozone optical depths CALL odo3(nz,z,nw,wl,o3xs,co3, dto3) * SO2 vertical profile and optical depths CALL setso2(so2col,nz,z,nw,wl,so2xs, $ tlay, cair, $ dtso2) * NO2 vertical profile and optical depths CALL setno2(no2col,nz,z,nw,wl,no2xs, $ tlay, cair, $ dtno2) * Cloud vertical profile, optical depths, single scattering albedo, asymmetry factor CALL setcld(taucld,zbase,ztop, $ nz,z,nw,wl, $ dtcld,omcld,gcld) * Aerosol vertical profile, optical depths, single scattering albedo, asymmetry factor CALL setaer(tauaer,ssaaer,alpha, $ nz,z,nw,wl, $ dtaer,omaer,gaer) * Surface albedo CALL setalb(alsurf,nw,wl, $ albedo) * ___ SECTION 5: TIME/SZA LOOP _____________________________________ * Initialize any time-integrated quantities here CALL zero1(dose,ks) * Loop over time or solar zenith angle (zen): DO 20, it = 1, nt zen = sza(it) WRITE(*,200) it, zen, esfact(it) WRITE(kout,200) it, zen, esfact(it) 200 FORMAT('step = ', I4,' sza = ', F9.3, $ ' Earth-sun factor = ', F10.7) * correction for earth-sun distance (SOLAR only, not LAMP) IF(.NOT.llamp) THEN DO iw = 1, nw - 1 etf(iw) = etf(iw) * esfact(it) ENDDO ENDIF * ____ SECTION 6: CALCULATE ZENITH ANGLE-DEPENDENT QUANTITIES __________ * slant path lengths for spherical geometry CALL sphers(nz,z,zen, dsdh,nid) CALL airmas(nz, dsdh,nid, cair,vcol,scol) * Recalculate effective O2 optical depth and cross sections for Lyman-alpha * and Schumann-Runge bands, must know zenith angle * Then assign O2 cross section to sj(1,*,*) CALL la_srb(nz,z,tlev,nw,wl,vcol,scol,o2xs1, $ dto2,o2xs) CALL sjo2(nz,nw,o2xs,1, sj) * ____ SECTION 7: WAVELENGTH LOOP ______________________________________ * initialize for wavelength integration CALL zero2(rate,ks,kz) CALL zero2(valj,kj,kz) * cba following lines added for saprc99 and SGMCM CALL zero2(valj4,kj,kz) CALL zero2(valj6,kj,kz) CALL zero2(valj8,kj,kz) ***** Main wavelength loop: DO 10, iw = 1, nw-1 ** monochromatic radiative transfer. Outputs are: * normalized irradiances edir(iz), edn(iz), eup(iz) * normalized actinic fluxes fdir(iz), fdn(zi), fup(iz) * where * dir = direct beam, dn = down-welling diffuse, up = up-welling diffuse !SM: for lamp, do not use directional components (dirsun, fdir, etc.) IF(.NOT.llamp) CALL rtlink(nstr, nz, $ iw, albedo(iw), zen, $ dsdh,nid, $ dtrl, $ dto3, $ dto2, $ dtso2, $ dtno2, $ dtcld, omcld, gcld, $ dtaer,omaer,gaer, $ edir, edn, eup, fdir, fdn, fup) * Spectral irradiance, W m-2 nm-1, down-welling: DO iz = 1, nz sirrad(iz,iw) = etf(iw) IF(.NOT.llamp)THEN sirrad(iz,iw) = etf(iw) * ( dirsun*edir(iz) + $ difdn*edn(iz) + difup*eup(iz) ) ELSE ENDIF ENDDO ! Spectra Actinic Flux, quanta m-2 s-1 DO iz = 1, nz ! LAMP, SOLAR: direct beam saflux(iz,iw) = etf(iw)* 5.039e11 * wc(iw) ! SOLAR adds directional (diffuse) components IF(.NOT.llamp) THEN saflux(iz,iw) = saflux(iz,iw) * $ (dirsun*fdir(iz) + difdn*fdn(iz) + difup*fup(iz)) ENDIF ENDDO *** Accumulate weighted integrals over wavelength, at all altitudes: DO 18 iz = 1, nz * Weighted irradiances (dose rates) W m-2 DO is = 1, ns drdw = sirrad(iz,iw) * sw(is,iw) rate(is,iz) = rate(is,iz) + drdw * (wu(iw) - wl(iw)) ENDDO * Photolysis rate coefficients (J-values) s-1 DO ij = 1, nj djdw = saflux(iz,iw) * sj(ij,iz,iw) valj(ij,iz) = valj(ij,iz) + djdw * (wu(iw) - wl(iw)) ENDDO * cba - loop for SGMCM DO ij = 1, nj4 djdw = saflux(iz,iw) * sj4(ij,iz,iw) valj4(ij,iz) = valj4(ij,iz) + djdw * (wu(iw) - wl(iw)) ENDDO DO ij = 1, nj6 djdw = saflux(iz,iw) * sj6(ij,iz,iw) valj6(ij,iz) = valj6(ij,iz) + djdw * (wu(iw) - wl(iw)) ENDDO DO ij = 1, nj8 djdw = saflux(iz,iw) * sj8(ij,iz,iw) valj8(ij,iz) = valj8(ij,iz) + djdw * (wu(iw) - wl(iw)) ENDDO 18 CONTINUE * Save irradiances and actinic fluxes for output CALL saver1(it, itfix, iw, iwfix, nz, izout, $ sirrad, saflux, $ svi_zw, svf_zw, svi_zt, svf_zt, svi_tw, svf_tw) 10 CONTINUE *^^^^^^^^^^^^^^^^ end wavelength loop * Save dose rates and j-values for output CALL saver2(it,itfix, nz,izout, ns,isfix,ims, nj,ijfix,imj, $ rate, valj, $ svr_zs, svj_zj, svr_zt, svj_zt, svr_ts, svj_tj) * cba - save the j-value data for izout DO ij=1,nj4 resu4(it,ij)=valj4(ij,izout) ENDDO DO ij=1,nj6 resu6(it,ij)=valj6(ij,izout) ENDDO DO ij=1,nj8 resu8(it,ij)=valj8(ij,izout) ENDDO 20 CONTINUE write(6,*) 'izout=', izout write(6,*) 'z(izout)', z(izout) *^^^^^^^^^^^^^^^^ end time/zenith loop ** reset wavelength scale if needed: IF(lrefr) THEN write(*,*) 'applying vacuum to air wavelength shift' $ , izout, airout mrefr = -mrefr CALL wshift(mrefr, nw, wl, airout) CALL wshift(mrefr, nwint, wc, airout) CALL wshift(mrefr, nwint, wu, airout) ENDIF ** output call outpt1( outfil, iout, $ lirrad, laflux, lrates, ljvals, lmmech, lzenit, $ nms, ims, nmj, imj, $ nz, z, tlev, airden, izout, $ nw, wl, etf, iwfix, $ nt, t, sza, itfix, $ ns, slabel, isfix, nj, jlabel, ijfix, $ svj_zj, svj_tj, svj_zt, $ svr_zs, svr_ts, svr_zt, $ svf_zw, svf_tw, svf_zt, $ svi_zw, svi_tw, svi_zt ) ! jmlt - jvalues @ sza=90 are not always zero. ! toto1 = 90.0 ! toto2 = 0. * cba - write results for MCM2001 * jmlt - and write results for difunctional nitrates DO ij=1,nj8 WRITE(90,'(a2)') '/ ' WRITE(90,'(a2,a40)') '/ ', jlabel8(ij) WRITE(90,'(a40,i2)') jlabel9(ij),nt DO it=1,nt WRITE(90,'(f6.1, 1pe12.4)') sza(it),resu8(it,ij) ENDDO ! WRITE(90,'(f6.1, 1pe12.4)') toto1, toto2 ENDDO DO ij=1,nj4 WRITE(90,'(a2)') '/ ' WRITE(90,'(a2,a40)') '/ ', jlabel4(ij) WRITE(90,'(a40,i2)') jlabel5(ij),nt DO it=1,nt WRITE(90,'(f6.1, 1pe12.4)') sza(it),resu4(it,ij) ENDDO ! WRITE(90,'(f6.1, 1pe12.4)') toto1, toto2 ENDDO DO ij=1,nj6 WRITE(90,'(a2)') '/ ' WRITE(90,'(a2,a40)') '/ ', jlabel6(ij) WRITE(90,'(a40,i2)') jlabel7(ij),nt DO it=1,nt WRITE(90,'(f6.1, 1pe12.4)') sza(it),resu6(it,ij) ENDDO ! WRITE(90,'(f6.1, 1pe12.4)') toto1, toto2 ENDDO WRITE(90,'(a4)') 'END ' *_______________________________________________________________________ IF(intrct) THEN WRITE(*,*) 'do you want to do another calculation?' WRITE(*,*) 'y = yes' WRITE(*,*) 'any other key = no' READ(*,1001) again 1001 FORMAT(A1) IF(again .EQ. 'y' .OR. again .EQ. 'Y') GO TO 1000 ENDIF CLOSE(iout) CLOSE(kout) END