!============================================================================== subroutine photolysis_online(nlayer,lswitch, alt, press, temp, $ zmmean, i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h, $ i_h2, i_oh, i_ho2, i_h2o2, i_h2o, $ i_n, i_n2d, i_no, i_no2, i_n2, nesp, rm, $ tau, sza, dist_sol, v_phot) implicit none #include "chimiedata.h" ! input integer :: nesp ! total number of chemical species integer :: nlayer, lswitch integer :: i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h, $ i_h2, i_oh, i_ho2, i_h2o2, i_h2o, $ i_n, i_n2d, i_no, i_no2, i_n2 real, dimension(nlayer) :: press, temp, zmmean ! pressure (hpa)/temperature (k)/mean molecular mass (g.mol-1) real, dimension(nlayer) :: alt ! altitude (km) real, dimension(nlayer,nesp) :: rm ! mixing ratios real :: tau ! integrated aerosol optical depth at the surface real :: sza ! solar zenith angle (degrees) real :: dist_sol ! solar distance (au) ! output real (kind = 8), dimension(nlayer,nb_phot_max) :: v_phot ! photolysis rates (s-1) ! spectral grid ! integer, parameter :: nw = 3789 ! number of spectral intervals (high-res) integer, parameter :: nw = 152 ! number of spectral intervals (low-res) real, dimension(nw) :: wl, wc, wu ! lower, center, upper wavelength for each interval save wl, wc, wu ! solar flux real, dimension(nw) :: f ! solar flux (w.m-2.nm-1) at 1 au real :: factor save f ! cross-sections and quantum yields integer, parameter :: nabs = 10 ! number of absorbing gases real, dimension(nlayer,nw,nd) :: sj ! general cross-section array (cm2) real, dimension(nw) :: xsco2_195, xsco2_295, xsco2_370 ! co2 absorption cross-section at 195-295-370 k (cm2) real, dimension(nw) :: yieldco2 ! co2 photodissociation yield real, dimension(nw) :: xso2_150, xso2_200, xso2_250, xso2_300 ! o2 absorption cross-section at 150-200-250-300 k (cm2) real, dimension(nw) :: yieldo2 ! o2 photodissociation yield real, dimension(nw) :: xso3_218, xso3_298 ! o3 absorption cross-section at 218-298 k (cm2) real, dimension(nw) :: xsh2o ! h2o absorption cross-section (cm2) real, dimension(nw) :: xsh2o2 ! h2o2 absorption cross-section (cm2) real, dimension(nw) :: xsho2 ! ho2 absorption cross-section (cm2) real, dimension(nw) :: xsh2 ! h2 absorption cross-section (cm2) real, dimension(nw) :: yieldh2 ! h2 photodissociation yield real, dimension(nw) :: xsno2, xsno2_220, xsno2_294 ! no2 absorption cross-section at 220-294 k (cm2) real, dimension(nw) :: yldno2_248, yldno2_298 ! no2 quantum yield at 248-298 k real, dimension(nw) :: xsno ! no absorption cross-section (cm2) real, dimension(nw) :: yieldno ! no photodissociation yield real, dimension(nw) :: xsn2 ! n2 absorption cross-section (cm2) real, dimension(nw) :: yieldn2 ! n2 photodissociation yield real, dimension(nw) :: albedo ! surface albedo ! atmosphere real, dimension(nlayer) :: colinc ! air column increment (molecule.cm-2) real, dimension(nlayer,nw) :: dtrl ! rayleigh optical depth real, dimension(nlayer,nw) :: dtaer ! aerosol optical depth real, dimension(nlayer,nw) :: omaer ! aerosol single scattering albedo real, dimension(nlayer,nw) :: gaer ! aerosol asymmetry parameter real, dimension(nlayer,nw) :: dtcld ! cloud optical depth real, dimension(nlayer,nw) :: omcld ! cloud single scattering albedo real, dimension(nlayer,nw) :: gcld ! cloud asymmetry parameter real, dimension(nlayer,nw,nabs) :: dtgas ! optical depth for each gas real, dimension(nlayer,nw) :: dagas ! total gas optical depth real, dimension(nlayer) :: edir, edn, eup ! normalised irradiances real, dimension(nlayer) :: fdir, fdn, fup ! normalised actinic fluxes real, dimension(nlayer) :: saflux ! total actinic flux integer, dimension(0:nlayer) :: nid real, dimension(0:nlayer,nlayer) :: dsdh integer :: j_o2_o, j_o2_o1d, j_co2_o, j_co2_o1d, j_o3_o1d, j_o3_o, $ j_h2o, j_h2o2, j_ho2, j_h2, j_no, j_no2, j_n2 integer :: a_o2, a_co2, a_o3, a_h2o, a_h2o2, a_ho2, a_h2, a_no, $ a_no2, a_n2 integer :: i, ilay, iw, ialt, mopt real :: deltaj logical, save :: firstcall data firstcall / .true. / ! high-res/low-res switch ! ! mopt = 1 high-resolution ! mopt = 2 low-resolution (recommended for gcm use) if (nw >= 3789) then mopt = 1 else mopt = 2 end if ! absorbing gas numbering a_o2 = 1 ! o2 a_co2 = 2 ! co2 a_o3 = 3 ! o3 a_h2o = 4 ! h2o a_h2o2 = 5 ! h2o2 a_ho2 = 6 ! ho2 a_h2 = 7 ! h2 a_no = 8 ! no a_no2 = 9 ! no2 a_n2 = 10 ! no2 ! photodissociation rates numbering j_o2_o = 1 ! o2 + hv -> o + o j_o2_o1d = 2 ! o2 + hv -> o + o(1d) j_co2_o = 3 ! co2 + hv -> co + o j_co2_o1d = 4 ! co2 + hv -> co + o(1d) j_o3_o1d = 5 ! o3 + hv -> o2 + o(1d) j_o3_o = 6 ! o3 + hv -> o2 + o j_h2o = 7 ! h2o + hv -> h + oh j_h2o2 = 8 ! h2o2 + hv -> oh + oh j_ho2 = 9 ! ho2 + hv -> oh + o j_h2 = 10 ! h2 + hv -> h + h j_no = 11 ! no + hv -> n + o j_no2 = 12 ! no2 + hv -> no + o j_n2 = 13 ! n2 + hv -> n + n !===== day/night criterion if (sza <= 120.) then !===== read once absorption cross-sections, quantum yields, and albedo if (firstcall) then ! set wavelength grid call gridw(nw,wl,wc,wu,mopt) ! read and grid solar flux data call rdsolarflux(nw,wl,wc,f) ! read and grid o2 cross-sections call rdxso2(nw,wl,xso2_150,xso2_200,xso2_250,xso2_300,yieldo2) ! read and grid co2 cross-sections call rdxsco2(nw,wl,xsco2_195,xsco2_295,xsco2_370,yieldco2) ! read and grid o3 cross-sections call rdxso3(nw,wl,xso3_218,xso3_298) ! read and grid h2o cross-sections call rdxsh2o(nw,wl,xsh2o) ! read and grid h2o2 cross-sections call rdxsh2o2(nw,wl,xsh2o2) ! read and grid ho2 cross-sections call rdxsho2(nw,wl,xsho2) ! read and grid h2 cross-sections call rdxsh2(nw,wl,wc,xsh2,yieldh2) ! read and grid no cross-sections call rdxsno(nw,wl,xsno,yieldno) ! read and grid no2 cross-sections call rdxsno2(nw,wl,xsno2,xsno2_220,xsno2_294, $ yldno2_248, yldno2_298) ! read and grid n2 cross-sections call rdxsn2(nw,wl,xsn2,yieldn2) ! set surface albedo call setalb(nw,wl,albedo) firstcall = .false. end if !==== air column increments and rayleigh optical depth call setair(nlayer, nw, wl, wc, press, temp, zmmean, colinc, dtrl) !==== set temperature-dependent cross-sections and optical depths ! o2: call seto2(nd, nlayer, nw, wc, mopt, temp, xso2_150, xso2_200, $ xso2_250, xso2_300, yieldo2, j_o2_o, j_o2_o1d, $ colinc, rm(:,i_o2), dtgas(:,:,a_o2), sj) ! co2: call setco2(nd, nlayer, nw, wc, temp, xsco2_195, xsco2_295, $ xsco2_370, yieldco2, j_co2_o, j_co2_o1d, $ colinc, rm(:,i_co2), dtgas(:,:,a_co2), sj) ! o3: call seto3(nd, nlayer, nw, wc, temp, xso3_218, xso3_298, $ j_o3_o, j_o3_o1d, $ colinc, rm(:,i_o3), dtgas(:,:,a_o3), sj) ! h2o2: call seth2o2(nd, nlayer, nw, wc, temp, xsh2o2, j_h2o2, $ colinc, rm(:,i_h2o2), dtgas(:,:,a_h2o2), sj) ! no2: call setno2(nd, nlayer, nw, wc, temp, xsno2, xsno2_220, $ xsno2_294, yldno2_248, yldno2_298, j_no2, $ colinc, rm(:,i_no2), dtgas(:,:,a_no2), sj) !==== temperature independent optical depths and cross-sections ! optical depths do ilay = 1,nlayer do iw = 1,nw-1 dtgas(ilay,iw,a_h2o) = colinc(ilay)*rm(ilay,i_h2o)*xsh2o(iw) dtgas(ilay,iw,a_ho2) = colinc(ilay)*rm(ilay,i_ho2)*xsho2(iw) dtgas(ilay,iw,a_h2) = colinc(ilay)*rm(ilay,i_h2)*xsh2(iw) dtgas(ilay,iw,a_no) = colinc(ilay)*rm(ilay,i_no)*xsno(iw) dtgas(ilay,iw,a_n2) = colinc(ilay)*rm(ilay,i_n2)*xsn2(iw) end do end do ! total gas optical depth dagas(:,:) = 0. do ilay = 1,nlayer do iw = 1,nw-1 do i = 1,nabs dagas(ilay,iw) = dagas(ilay,iw) + dtgas(ilay,iw,i) end do end do end do ! cross-sections do ilay = 1,nlayer do iw = 1,nw-1 sj(ilay,iw,j_h2o) = xsh2o(iw) ! h2o sj(ilay,iw,j_ho2) = xsho2(iw) ! ho2 sj(ilay,iw,j_h2) = xsh2(iw)*yieldh2(iw) ! h2 sj(ilay,iw,j_no) = xsno(iw)*yieldno(iw) ! no sj(ilay,iw,j_n2) = xsn2(iw)*yieldn2(iw) ! n2 end do end do ! ialt = 49 ! print*, 'altitude = ', alt(ialt), ' km' ! do iw = 1,nw ! write(35,'(f8.4,12e13.5)') wc(iw), dtrl(ialt,iw), ! $ (max(dtgas(ialt,iw,i),1.e-30),i = 1,nabs) ! end do !==== set aerosol properties and optical depth call setaer(nlayer,alt,tau,nw,dtaer,omaer,gaer) !==== set cloud properties and optical depth call setcld(nlayer,nw,dtcld,omcld,gcld) !==== slant path lengths in spherical geometry call sphers(nlayer, alt, sza, dsdh, nid) !==== solar flux at mars factor = (1./dist_sol)**2. do iw = 1,nw-1 f(iw) = f(iw)*factor end do !==== initialise photolysis rates v_phot(:,1:nd) = 0. !==== start of wavelength lopp do iw = 1,nw-1 ! monochromatic radiative transfer. outputs are: ! normalized irradiances edir(nlayer), edn(nlayer), eup(nlayer) ! normalized actinic fluxes fdir(nlayer), fdn(nlayer), fup(nlayer) ! where ! dir = direct beam, dn = down-welling diffuse, up = up-welling diffuse call rtlink(nlayer, nw, iw, albedo(iw), sza, dsdh, nid, dtrl, $ dagas, dtcld, omcld, gcld, dtaer, omaer, gaer, $ edir, edn, eup, fdir, fdn, fup) ! spherical actinic flux do ilay = 1,lswitch-1 saflux(ilay) = f(iw)*(fdir(ilay) + fdn(ilay) + fup(ilay)) end do ! photolysis rate integration do i = 1,nd do ilay = 1, lswitch-1 deltaj = saflux(ilay)*sj(ilay,iw,i) v_phot(ilay,i) = v_phot(ilay,i) + deltaj*(wu(iw)-wl(iw)) end do end do ! eliminate small values where (v_phot(:,1:nd) < 1.e-30) v_phot(:,1:nd) = 0. end where end do ! iw else ! night v_phot(:,1:nd) = 0. end if ! sza end subroutine photolysis_online !============================================================================== subroutine gridw(nw,wl,wc,wu,mopt) ! 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. implicit none ! input integer :: nw ! number of wavelength grid points integer :: mopt ! high-res/low-res switch ! output real, dimension(nw) :: wl, wc, wu ! lower, center, upper wavelength for each interval ! local real :: wincr ! wavelength increment integer :: iw, kw ! mopt = 1 high-resolution mode (3789 intervals) ! ! 0-108 nm : 1.0 nm ! 108-124 nm : 0.1 nm ! 124-175 nm : 0.5 nm ! 175-205 nm : 0.01 nm ! 205-365 nm : 0.5 nm ! 365-850 nm : 5.0 nm ! ! mopt = 2 low-resolution mode ! ! 0-60 nm : 6.0 nm ! 60-80 nm : 2.0 nm ! 80-120 nm : 5.0 nm ! 120-123 nm : 0.2 nm ! 123-163 nm : 5.0 nm ! 163-175 nm : 2.0 nm ! 175-205 nm : 0.5 nm ! 205-245 nm : 5.0 nm ! 245-415 nm : 10.0 nm ! 415-815 nm : 50.0 nm if (mopt == 1) then ! high-res ! define wavelength intervals of width 1.0 nm from 0 to 108 nm: kw = 0 wincr = 1.0 do iw = 0, 107 kw = kw + 1 wl(kw) = real(iw) wu(kw) = wl(kw) + wincr wc(kw) = (wl(kw) + wu(kw))/2. end do ! define wavelength intervals of width 0.1 nm from 108 to 124 nm: wincr = 0.1 do iw = 1080, 1239, 1 kw = kw + 1 wl(kw) = real(iw)/10. wu(kw) = wl(kw) + wincr wc(kw) = (wl(kw) + wu(kw))/2. end do ! define wavelength intervals of width 0.5 nm from 124 to 175 nm: wincr = 0.5 do iw = 1240, 1745, 5 kw = kw + 1 wl(kw) = real(iw)/10. wu(kw) = wl(kw) + wincr wc(kw) = (wl(kw) + wu(kw))/2. end do ! define wavelength intervals of width 0.01 nm from 175 to 205 nm: wincr = 0.01 do iw = 17500, 20499, 1 kw = kw + 1 wl(kw) = real(iw)/100. wu(kw) = wl(kw) + wincr wc(kw) = (wl(kw) + wu(kw))/2. end do ! define wavelength intervals of width 0.5 nm from 205 to 365 nm: wincr = 0.5 do iw = 2050, 3645, 5 kw = kw + 1 wl(kw) = real(iw)/10. wu(kw) = wl(kw) + wincr wc(kw) = (wl(kw) + wu(kw))/2. end do ! define wavelength intervals of width 5.0 nm from 365 to 855 nm: wincr = 5.0 do iw = 365, 850, 5 kw = kw + 1 wl(kw) = real(iw) wu(kw) = wl(kw) + wincr wc(kw) = (wl(kw) + wu(kw))/2. end do wl(kw+1) = wu(kw) else if (mopt == 2) then ! low-res * define wavelength intervals of width 6.0 nm from 0 to 60 nm: kw = 0 wincr = 6.0 DO iw = 0, 54, 6 kw = kw + 1 wl(kw) = real(iw) wu(kw) = wl(kw) + wincr wc(kw) = (wl(kw) + wu(kw))/2. END DO * define wavelength intervals of width 2.0 nm from 60 to 80 nm: wincr = 2.0 DO iw = 60, 78, 2 kw = kw + 1 wl(kw) = real(iw) wu(kw) = wl(kw) + wincr wc(kw) = (wl(kw) + wu(kw))/2. END DO * define wavelength intervals of width 5.0 nm from 80 to 120 nm: wincr = 5.0 DO iw = 80, 115, 5 kw = kw + 1 wl(kw) = real(iw) wu(kw) = wl(kw) + wincr wc(kw) = (wl(kw) + wu(kw))/2. END DO * define wavelength intervals of width 0.2 nm from 120 to 123 nm: wincr = 0.2 DO iw = 1200, 1228, 2 kw = kw + 1 wl(kw) = real(iw)/10. wu(kw) = wl(kw) + wincr wc(kw) = (wl(kw) + wu(kw))/2. ENDDO * define wavelength intervals of width 5.0 nm from 123 to 163 nm: wincr = 5.0 DO iw = 123, 158, 5 kw = kw + 1 wl(kw) = real(iw) wu(kw) = wl(kw) + wincr wc(kw) = (wl(kw) + wu(kw))/2. ENDDO * define wavelength intervals of width 2.0 nm from 163 to 175 nm: wincr = 2.0 DO iw = 163, 173, 2 kw = kw + 1 wl(kw) = real(iw) wu(kw) = wl(kw) + wincr wc(kw) = (wl(kw) + wu(kw))/2. ENDDO * define wavelength intervals of width 0.5 nm from 175 to 205 nm: wincr = 0.5 DO iw = 1750, 2045, 5 kw = kw + 1 wl(kw) = real(iw)/10. wu(kw) = wl(kw) + wincr wc(kw) = (wl(kw) + wu(kw))/2. ENDDO * define wavelength intervals of width 5.0 nm from 205 to 245 nm: wincr = 5. DO iw = 205, 240, 5 kw = kw + 1 wl(kw) = real(iw) wu(kw) = wl(kw) + wincr wc(kw) = (wl(kw) + wu(kw))/2. ENDDO * define wavelength intervals of width 10.0 nm from 245 to 415 nm: wincr = 10.0 DO iw = 245, 405, 10 kw = kw + 1 wl(kw) = real(iw) wu(kw) = wl(kw) + wincr wc(kw) = (wl(kw) + wu(kw))/2. ENDDO * define wavelength intervals of width 50.0 nm from 415 to 815 nm: wincr = 50.0 DO iw = 415, 815, 50 kw = kw + 1 wl(kw) = real(iw) wu(kw) = wl(kw) + wincr wc(kw) = (wl(kw) + wu(kw))/2. ENDDO wl(kw+1) = wu(kw) end if ! mopt ! do iw = 1,nw ! write(20,*) iw, wl(iw), wu(iw) ! end do print*, 'number of spectral intervals : ', kw+1 return end !============================================================================== subroutine rdsolarflux(nw,wl,wc,f) ! Read and re-grid solar flux data. use datafile_mod, only: datadir implicit none ! input integer :: nw ! number of wavelength grid points real, dimension(nw) :: wl, wc ! lower and central wavelength for each interval ! output real, dimension(nw) :: f ! solar flux (w.m-2.nm-1) ! local integer, parameter :: kdata = 20000 ! max dimension of input solar flux integer :: msun ! choice of solar flux integer :: iw, nhead, ihead, n, i, ierr, kin real, parameter :: deltax = 1.e-4 real, dimension(kdata) :: x1, y1 ! input solar flux real, dimension(nw) :: yg1 ! gridded solar flux character(len=100) :: fil kin = 10 ! input logical unit ! select desired extra-terrestrial solar irradiance, using msun: ! 18 = atlas3_thuillier_tuv.txt 0-900 nm November 1994 ! Thuillier et al., Adv. Space. Res., 34, 256-261, 2004 msun = 18 if (msun == 18) THEN fil = trim(datadir)//'/solar_fluxes/atlas3_thuillier_tuv.txt' print*, 'solar flux : ', fil OPEN(UNIT=kin,FILE=fil,STATUS='old') nhead = 9 n = 19193 DO ihead = 1, nhead READ(kin,*) ENDDO DO i = 1, n READ(kin,*) x1(i), y1(i) y1(i) = y1(i)*1.e-3 ! mw -> w ENDDO CLOSE (kin) CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.) CALL addpnt(x1,y1,kdata,n, 0.,0.) CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.) CALL addpnt(x1,y1,kdata,n, 1.e+38,0.) CALL inter2(nw,wl,yg1,n,x1,y1,ierr) IF (ierr .NE. 0) THEN WRITE(*,*) ierr, fil STOP ENDIF ! convert to photon.s-1.nm-1.cm-2 ! 5.039e11 = 1.e-4*1e-9/(hc = 6.62e-34*2.998e8) DO iw = 1, nw-1 f(iw) = yg1(iw)*wc(iw)*5.039e11 ! write(25,*) iw, wc(iw), f(iw) ENDDO end if end subroutine rdsolarflux !============================================================================== subroutine addpnt ( x, y, ld, n, xnew, ynew ) *-----------------------------------------------------------------------------* *= PURPOSE: =* *= Add a point to a set of data pairs . x must be in =* *= ascending order =* *-----------------------------------------------------------------------------* *= PARAMETERS: =* *= X - REAL vector of length LD, x-coordinates (IO)=* *= Y - REAL vector of length LD, y-values (IO)=* *= LD - INTEGER, dimension of X, Y exactly as declared in the calling (I)=* *= program =* *= N - INTEGER, number of elements in X, Y. On entry, it must be: (IO)=* *= N < LD. On exit, N is incremented by 1. =* *= XNEW - REAL, x-coordinate at which point is to be added (I)=* *= YNEW - REAL, y-value of point to be added (I)=* *-----------------------------------------------------------------------------* IMPLICIT NONE C calling parameters INTEGER ld, n REAL x(ld), y(ld) REAL xnew, ynew INTEGER ierr C local variables INTEGER insert INTEGER i C----------------------------------------------------------------------- * initialize error flag ierr = 0 * check n>> ERROR (ADDPNT) <<< Cannot expand array ' WRITE(0,*) ' All elements used.' STOP ENDIF insert = 1 i = 2 * check, whether x is already sorted. * also, use this loop to find the point at which xnew needs to be inserted * into vector x, if x is sorted. 10 CONTINUE IF (i .LT. n) THEN IF (x(i) .LT. x(i-1)) THEN print*, x(i-1), x(i) WRITE(0,*) '>>> ERROR (ADDPNT) <<< x-data must be '// > 'in ascending order!' STOP ELSE IF (xnew .GT. x(i)) insert = i + 1 ENDIF i = i+1 GOTO 10 ENDIF * if needs to be appended at the end, just do so, * otherwise, insert at position INSERT IF ( xnew .GT. x(n) ) THEN x(n+1) = xnew y(n+1) = ynew ELSE * shift all existing points one index up DO i = n, insert, -1 x(i+1) = x(i) y(i+1) = y(i) ENDDO * insert new point x(insert) = xnew y(insert) = ynew ENDIF * increase total number of elements in x, y n = n+1 END c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c subroutine inter2(ng,xg,yg,n,x,y,ierr) *-----------------------------------------------------------------------------* *= PURPOSE: =* *= Map input data given on single, discrete points onto a set of target =* *= bins. =* *= The original input data are given on single, discrete points of an =* *= arbitrary grid and are being linearly interpolated onto a specified set =* *= of target bins. In general, this is the case for most of the weighting =* *= functions (action spectra, molecular cross section, and quantum yield =* *= data), which have to be matched onto the specified wavelength intervals. =* *= The average value in each target bin is found by averaging the trapezoi- =* *= dal area underneath the input data curve (constructed by linearly connec-=* *= ting the discrete input values). =* *= Some caution should be used near the endpoints of the grids. If the =* *= input data set does not span the range of the target grid, an error =* *= message is printed and the execution is stopped, as extrapolation of the =* *= data is not permitted. =* *= If the input data does not encompass the target grid, use ADDPNT to =* *= expand the input array. =* *-----------------------------------------------------------------------------* *= PARAMETERS: =* *= NG - INTEGER, number of bins + 1 in the target grid (I)=* *= XG - REAL, target grid (e.g., wavelength grid); bin i is defined (I)=* *= as [XG(i),XG(i+1)] (i = 1..NG-1) =* *= YG - REAL, y-data re-gridded onto XG, YG(i) specifies the value for (O)=* *= bin i (i = 1..NG-1) =* *= N - INTEGER, number of points in input grid (I)=* *= X - REAL, grid on which input data are defined (I)=* *= Y - REAL, input y-data (I)=* *-----------------------------------------------------------------------------* IMPLICIT NONE * input: INTEGER ng, n REAL x(n), y(n), xg(ng) * output: REAL yg(ng) * local: REAL area, xgl, xgu REAL darea, slope REAL a1, a2, b1, b2 INTEGER ngintv INTEGER i, k, jstart INTEGER ierr *_______________________________________________________________________ ierr = 0 * test for correct ordering of data, by increasing value of x DO 10, i = 2, n IF (x(i) .LE. x(i-1)) THEN ierr = 1 WRITE(*,*)'data not sorted' WRITE(*,*) x(i), x(i-1) RETURN ENDIF 10 CONTINUE DO i = 2, ng IF (xg(i) .LE. xg(i-1)) THEN ierr = 2 WRITE(0,*) '>>> ERROR (inter2) <<< xg-grid not sorted!' RETURN ENDIF ENDDO * check for xg-values outside the x-range IF ( (x(1) .GT. xg(1)) .OR. (x(n) .LT. xg(ng)) ) THEN WRITE(0,*) '>>> ERROR (inter2) <<< Data do not span '// > 'grid. ' WRITE(0,*) ' Use ADDPNT to '// > 'expand data and re-run.' STOP ENDIF * find the integral of each grid interval and use this to * calculate the average y value for the interval * xgl and xgu are the lower and upper limits of the grid interval jstart = 1 ngintv = ng - 1 DO 50, i = 1,ngintv * initalize: area = 0.0 xgl = xg(i) xgu = xg(i+1) * discard data before the first grid interval and after the * last grid interval * for internal grid intervals, start calculating area by interpolating * between the last point which lies in the previous interval and the * first point inside the current interval k = jstart IF (k .LE. n-1) THEN * if both points are before the first grid, go to the next point 30 CONTINUE IF (x(k+1) .LE. xgl) THEN jstart = k - 1 k = k+1 IF (k .LE. n-1) GO TO 30 ENDIF * if the last point is beyond the end of the grid, complete and go to the next * grid 40 CONTINUE IF ((k .LE. n-1) .AND. (x(k) .LT. xgu)) THEN jstart = k-1 * compute x-coordinates of increment a1 = MAX(x(k),xgl) a2 = MIN(x(k+1),xgu) * if points coincide, contribution is zero IF (x(k+1).EQ.x(k)) THEN darea = 0.e0 ELSE slope = (y(k+1) - y(k))/(x(k+1) - x(k)) b1 = y(k) + slope*(a1 - x(k)) b2 = y(k) + slope*(a2 - x(k)) darea = (a2 - a1)*(b2 + b1)/2. ENDIF * find the area under the trapezoid from a1 to a2 area = area + darea * go to next point k = k+1 GO TO 40 ENDIF ENDIF * calculate the average y after summing the areas in the interval yg(i) = area/(xgu - xgl) 50 CONTINUE RETURN END !============================================================================== subroutine rdxsco2(nw,wl,xsco2_195,xsco2_295,xsco2_370,yieldco2) *-----------------------------------------------------------------------------* *= PURPOSE: =* *= Read and grid CO2 absorption cross-sections and photodissociation yield =* *-----------------------------------------------------------------------------* *= PARAMETERS: =* *= NW - INTEGER, number of specified intervals + 1 in working (I)=* *= wavelength grid =* *= WL - REAL, vector of lower limits of wavelength intervals in (I)=* *= working wavelength grid =* *= XSCO2 - REAL, molecular absoprtion cross section (cm^2) of CO2 at (O)=* *= each specified wavelength =* *-----------------------------------------------------------------------------* use datafile_mod, only: datadir implicit none ! input integer :: nw ! number of wavelength grid points real, dimension(nw) :: wl ! lower and central wavelength for each interval ! output real, dimension(nw) :: xsco2_195, xsco2_295, xsco2_370 ! co2 cross-sections (cm2) real, dimension(nw) :: yieldco2 ! co2 photodissociation yield ! local integer, parameter :: kdata = 36000 real, parameter :: deltax = 1.e-4 real, dimension(kdata) :: x1, y1, y2, y3, xion, ion real, dimension(nw) :: yg real :: xl, xu integer :: ierr, i, l, n, n1, n2, n3, n4 CHARACTER*100 fil integer :: kin, kout ! input/ouput logical units kin = 10 kout = 30 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c CO2 absorption cross-sections c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c 195K: Chan et al. (1993) + Starck et al. (2006) c + Yoshino et al. (1996) + Parkinson et al. (2003) + extrapolation c 11 lignes d'en-tete et n1 = 34537 c fichier: co2_euv_uv_195k.txt c c 295K: Chan et al. (1993) + Starck et al. (2006) c + Yoshino et al. (1996) + Parkinson et al. (2003) + extrapolation c 11 lignes d'en-tete et n2 = 35360 c fichier: co2_euv_uv_295k.txt c c 370K: Chan et al. (1993) + Starck et al. (2006) c + Lewis and Carver (1983) + extrapolation c 11 lignes d'en-tete et n3 = 3884 c fichier: co2_euv_uv_370k.txt c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c n1 = 34537 n2 = 35360 n3 = 3884 c c 195K: c fil = trim(datadir)//'/cross_sections/co2_euv_uv_195k.txt' print*, 'section efficace CO2 195K: ', fil OPEN(UNIT=kin,FILE=fil,STATUS='old') DO i = 1,11 read(kin,*) END DO DO i = 1, n1 READ(kin,*) x1(i), y1(i) END DO CLOSE (kin) CALL addpnt(x1,y1,kdata,n1,x1(1)*(1.-deltax),0.) CALL addpnt(x1,y1,kdata,n1, 0.,0.) CALL addpnt(x1,y1,kdata,n1,x1(n1)*(1.+deltax),0.) CALL addpnt(x1,y1,kdata,n1, 1.e+38,0.) CALL inter2(nw,wl,yg,n1,x1,y1,ierr) IF (ierr .NE. 0) THEN WRITE(*,*) ierr, fil STOP ENDIF DO l = 1, nw-1 xsco2_195(l) = yg(l) END DO c c 295K: c fil = trim(datadir)//'/cross_sections/co2_euv_uv_295k.txt' print*, 'section efficace CO2 295K: ', fil c OPEN(UNIT=kin,FILE=fil,STATUS='old') DO i = 1,11 read(kin,*) END DO c DO i = 1, n2 READ(kin,*) x1(i), y1(i) END DO CLOSE (kin) CALL addpnt(x1,y1,kdata,n2,x1(1)*(1.-deltax),0.) CALL addpnt(x1,y1,kdata,n2, 0.,0.) CALL addpnt(x1,y1,kdata,n2,x1(n2)*(1.+deltax),0.) CALL addpnt(x1,y1,kdata,n2, 1.e+38,0.) CALL inter2(nw,wl,yg,n2,x1,y1,ierr) IF (ierr .NE. 0) THEN WRITE(*,*) ierr, fil STOP ENDIF c DO l = 1, nw-1 xsco2_295(l) = yg(l) END DO c c 370K: c fil = trim(datadir)//'/cross_sections/co2_euv_uv_370k.txt' print*, 'section efficace CO2 370K: ', fil c OPEN(UNIT=kin,FILE=fil,STATUS='old') DO i = 1,11 read(kin,*) END DO c DO i = 1, n3 READ(kin,*) x1(i), y1(i) END DO CLOSE (kin) c CALL addpnt(x1,y1,kdata,n3,x1(1)*(1.-deltax),0.) CALL addpnt(x1,y1,kdata,n3, 0.,0.) CALL addpnt(x1,y1,kdata,n3,x1(n3)*(1.+deltax),0.) CALL addpnt(x1,y1,kdata,n3, 1.e+38,0.) CALL inter2(nw,wl,yg,n3,x1,y1,ierr) IF (ierr .NE. 0) THEN WRITE(*,*) ierr, fil STOP ENDIF DO l = 1, nw-1 xsco2_370(l) = yg(l) END DO ! photodissociation yield: fil = trim(datadir)// $'/cross_sections/efdis_co2-o2_schunkandnagy2000.txt' print*, 'photodissociation yield CO2: ', fil OPEN(UNIT=kin,FILE=fil,STATUS='old') do i = 1,3 read(kin,*) end do n4 = 17 do i = 1, n4 read(kin,*) xl, xu, ion(i) xion(i) = (xl + xu)/2. ion(i) = max(ion(i), 0.) end do close(kin) CALL addpnt(xion,ion,kdata,n4,xion(1)*(1.-deltax),0.) CALL addpnt(xion,ion,kdata,n4, 0.,0.) CALL addpnt(xion,ion,kdata,n4,xion(n4)*(1.+deltax),1.) CALL addpnt(xion,ion,kdata,n4, 1.e+38,1.) CALL inter2(nw,wl,yieldco2,n4,xion,ion,ierr) IF (ierr .NE. 0) THEN WRITE(*,*) ierr, fil STOP ENDIF ! DO l = 1, nw-1 ! write(kout,*) wl(l), xsco2_195(l), ! $ xsco2_295(l), ! $ xsco2_370(l), ! $ yieldco2(l) ! END DO end subroutine rdxsco2 !============================================================================== subroutine rdxso2(nw,wl,xso2_150,xso2_200,xso2_250,xso2_300, $ yieldo2) *-----------------------------------------------------------------------------* *= PURPOSE: =* *= Read and grid O2 cross-sections and photodissociation yield =* *-----------------------------------------------------------------------------* *= PARAMETERS: =* *= NW - INTEGER, number of specified intervals + 1 in working (I)=* *= wavelength grid =* *= WL - REAL, vector of lower limits of wavelength intervals in (I)=* *= working wavelength grid =* *= XSO2 - REAL, molecular absorption cross section =* *-----------------------------------------------------------------------------* use datafile_mod, only: datadir implicit none ! input integer :: nw ! number of wavelength grid points real, dimension(nw) :: wl ! lower and central wavelength for each interval ! output real, dimension(nw) :: xso2_150, xso2_200, xso2_250, xso2_300 ! o2 cross-sections (cm2) real, dimension(nw) :: yieldo2 ! o2 photodissociation yield ! local integer, parameter :: kdata = 18000 real, parameter :: deltax = 1.e-4 real, dimension(kdata) :: x1, y1, x2, y2, x3, y3, x4, y4 real, dimension(kdata) :: xion, ion real :: factor, xl, xu, dummy integer :: i, ierr, n, n1, n2, n3, n4, nhead integer :: kin, kout ! input/output logical units character*100 fil kin = 10 kout = 30 ! read o2 cross section data nhead = 22 n = 17434 fil = trim(datadir)//'/cross_sections/o2_composite_2018_150K.txt' print*, 'section efficace O2 150K: ', fil OPEN(UNIT=kin,FILE=fil,STATUS='old') DO i = 1,nhead read(kin,*) END DO n1 = n DO i = 1, n1 READ(kin,*) x1(i), y1(i) END DO CLOSE (kin) CALL addpnt(x1,y1,kdata,n1,x1(1)*(1.-deltax),0.) CALL addpnt(x1,y1,kdata,n1, 0.,0.) CALL addpnt(x1,y1,kdata,n1,x1(n1)*(1.+deltax),0.) CALL addpnt(x1,y1,kdata,n1, 1.e+38,0.) CALL inter2(nw,wl,xso2_150,n1,x1,y1,ierr) IF (ierr .NE. 0) THEN WRITE(*,*) ierr, fil STOP ENDIF fil = trim(datadir)//'/cross_sections/o2_composite_2018_200K.txt' print*, 'section efficace O2 200K: ', fil OPEN(UNIT=kin,FILE=fil,STATUS='old') DO i = 1,nhead read(kin,*) END DO n2 = n DO i = 1, n2 READ(kin,*) x2(i), y2(i) END DO CLOSE (kin) CALL addpnt(x2,y2,kdata,n2,x2(1)*(1.-deltax),0.) CALL addpnt(x2,y2,kdata,n2, 0.,0.) CALL addpnt(x2,y2,kdata,n2,x2(n2)*(1.+deltax),0.) CALL addpnt(x2,y2,kdata,n2, 1.e+38,0.) CALL inter2(nw,wl,xso2_200,n2,x2,y2,ierr) IF (ierr .NE. 0) THEN WRITE(*,*) ierr, fil STOP ENDIF fil = trim(datadir)//'/cross_sections/o2_composite_2018_250K.txt' print*, 'section efficace O2 250K: ', fil OPEN(UNIT=kin,FILE=fil,STATUS='old') DO i = 1,nhead read(kin,*) END DO n3 = n DO i = 1, n3 READ(kin,*) x3(i), y3(i) END DO CLOSE (kin) CALL addpnt(x3,y3,kdata,n3,x3(1)*(1.-deltax),0.) CALL addpnt(x3,y3,kdata,n3, 0.,0.) CALL addpnt(x3,y3,kdata,n3,x3(n3)*(1.+deltax),0.) CALL addpnt(x3,y3,kdata,n3, 1.e+38,0.) CALL inter2(nw,wl,xso2_250,n3,x3,y3,ierr) IF (ierr .NE. 0) THEN WRITE(*,*) ierr, fil STOP ENDIF fil = trim(datadir)//'/cross_sections/o2_composite_2018_300K.txt' print*, 'section efficace O2 300K: ', fil OPEN(UNIT=kin,FILE=fil,STATUS='old') DO i = 1,nhead read(kin,*) END DO n4 = n DO i = 1, n4 READ(kin,*) x4(i), y4(i) END DO CLOSE (kin) CALL addpnt(x4,y4,kdata,n4,x4(1)*(1.-deltax),0.) CALL addpnt(x4,y4,kdata,n4, 0.,0.) CALL addpnt(x4,y4,kdata,n4,x4(n4)*(1.+deltax),0.) CALL addpnt(x4,y4,kdata,n4, 1.e+38,0.) CALL inter2(nw,wl,xso2_300,n4,x4,y4,ierr) IF (ierr .NE. 0) THEN WRITE(*,*) ierr, fil STOP ENDIF ! photodissociation yield fil = trim(datadir)// $ '/cross_sections/efdis_co2-o2_schunkandnagy2000.txt' OPEN(UNIT=kin,FILE=fil,STATUS='old') do i = 1,11 read(kin,*) end do n = 9 DO i = 1, n READ(kin,*) xl, xu, dummy, ion(i) xion(i) = (xl + xu)/2. ion(i) = max(ion(i), 0.) END DO CLOSE (kin) CALL addpnt(xion,ion,kdata,n,xion(1)*(1.-deltax),0.) CALL addpnt(xion,ion,kdata,n, 0.,0.) CALL addpnt(xion,ion,kdata,n,xion(n)*(1.+deltax),1.) CALL addpnt(xion,ion,kdata,n, 1.e+38,1.) CALL inter2(nw,wl,yieldo2,n,xion,ion,ierr) IF (ierr .NE. 0) THEN WRITE(*,*) ierr, fil STOP ENDIF end subroutine rdxso2 !============================================================================== subroutine rdxso3(nw,wl,xso3_218,xso3_298) *-----------------------------------------------------------------------------* *= PURPOSE: =* *= Read ozone molecular absorption cross section. Re-grid data to match =* *= specified wavelength working grid. =* *-----------------------------------------------------------------------------* *= PARAMETERS: =* *= NW - INTEGER, number of specified intervals + 1 in working (I)=* *= wavelength grid =* *= WL - REAL, vector of lower limits of wavelength intervals in (I)=* *= working wavelength grid =* *= XSO3_218 REAL, molecular absoprtion cross section (cm^2) of O3 at (O)=* *= each specified wavelength (JPL 2006) 218 K =* *= XSO3_298 REAL, molecular absoprtion cross section (cm^2) of O3 at (O)=* *= each specified wavelength (JPL 2006) 298 K =* *-----------------------------------------------------------------------------* use datafile_mod, only: datadir implicit none ! input integer :: nw ! number of wavelength grid points real, dimension(nw) :: wl ! lower and central wavelength for each interval ! output real, dimension(nw) :: xso3_218, xso3_298 ! o3 cross-sections (cm2) ! local integer, parameter :: kdata = 200 real, parameter :: deltax = 1.e-4 real, dimension(kdata) :: x1, x2, y1, y2 real, dimension(nw) :: yg real :: a1, a2 integer :: i, ierr, iw, n, n1, n2 integer :: kin, kout ! input/output logical units character*100 fil ! JPL 2006 218 K fil = trim(datadir)// $ '/cross_sections/o3_cross-sections_jpl_2006_218K.txt' print*, 'section efficace O3 218K: ', fil OPEN(UNIT=kin,FILE=fil,STATUS='old') n1 = 167 DO i = 1, n1 READ(kin,*) a1, a2, y1(i) x1(i) = (a1+a2)/2. END DO CLOSE (kin) CALL addpnt(x1,y1,kdata,n1,x1(1)*(1.-deltax),0.) CALL addpnt(x1,y1,kdata,n1, 0.,0.) CALL addpnt(x1,y1,kdata,n1,x1(n1)*(1.+deltax),0.) CALL addpnt(x1,y1,kdata,n1, 1.e+38,0.) CALL inter2(nw,wl,yg,n1,x1,y1,ierr) IF (ierr .NE. 0) THEN WRITE(*,*) ierr, fil STOP ENDIF c DO iw = 1, nw-1 xso3_218(iw) = yg(iw) END DO ! JPL 2006 298 K fil = trim(datadir)// $ '/cross_sections/o3_cross-sections_jpl_2006_298K.txt' print*, 'section efficace O3 298K: ', fil OPEN(UNIT=kin,FILE=fil,STATUS='old') n2 = 167 DO i = 1, n2 READ(kin,*) a1, a2, y2(i) x2(i) = (a1+a2)/2. END DO CLOSE (kin) CALL addpnt(x2,y2,kdata,n2,x2(1)*(1.-deltax),0.) CALL addpnt(x2,y2,kdata,n2, 0.,0.) CALL addpnt(x2,y2,kdata,n2,x2(n2)*(1.+deltax),0.) CALL addpnt(x2,y2,kdata,n2, 1.e+38,0.) CALL inter2(nw,wl,yg,n2,x2,y2,ierr) IF (ierr .NE. 0) THEN WRITE(*,*) ierr, fil STOP ENDIF c DO iw = 1, nw-1 xso3_298(iw) = yg(iw) END DO end subroutine rdxso3 !============================================================================== subroutine rdxsh2o(nw, wl, yg) *-----------------------------------------------------------------------------* *= PURPOSE: =* *= Read H2O molecular absorption cross section. Re-grid data to match =* *= specified wavelength working grid. =* *-----------------------------------------------------------------------------* *= PARAMETERS: =* *= NW - INTEGER, number of specified intervals + 1 in working (I)=* *= wavelength grid =* *= WL - REAL, vector of lower limits of wavelength intervals in (I)=* *= working wavelength grid =* *= YG - REAL, molecular absoprtion cross section (cm^2) of H2O at (O)=* *= each specified wavelength =* *-----------------------------------------------------------------------------* use datafile_mod, only: datadir IMPLICIT NONE ! input integer :: nw ! number of wavelength grid points real, dimension(nw) :: wl ! lower and central wavelength for each interval ! output real, dimension(nw) :: yg ! h2o cross-sections (cm2) ! local integer, parameter :: kdata = 500 real, parameter :: deltax = 1.e-4 REAL x1(kdata) REAL y1(kdata) INTEGER ierr INTEGER i, n CHARACTER*100 fil integer :: kin, kout ! input/output logical units kin = 10 fil = trim(datadir)//'/cross_sections/h2o_composite_250K.txt' print*, 'section efficace H2O: ', fil OPEN(UNIT=kin,FILE=fil,STATUS='old') DO i = 1,26 read(kin,*) END DO n = 420 DO i = 1, n READ(kin,*) x1(i), y1(i) END DO CLOSE (kin) CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.) CALL addpnt(x1,y1,kdata,n, 0.,0.) CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.) CALL addpnt(x1,y1,kdata,n, 1.e+38,0.) CALL inter2(nw,wl,yg,n,x1,y1,ierr) IF (ierr .NE. 0) THEN WRITE(*,*) ierr, fil STOP ENDIF end subroutine rdxsh2o !============================================================================== subroutine rdxsh2o2(nw, wl, xsh2o2) *-----------------------------------------------------------------------------* *= PURPOSE: =* *= Read and grid H2O2 cross-sections *= H2O2 + hv -> 2 OH =* *= Cross section: Schuergers and Welge, Z. Naturforsch. 23a (1968) 1508 =* *= from 125 to 185 nm, then JPL97 from 190 to 350 nm. =* *-----------------------------------------------------------------------------* *= PARAMETERS: =* *= NW - INTEGER, number of specified intervals + 1 in working (I)=* *= wavelength grid =* *= WL - REAL, vector of lower limits of wavelength intervals in (I)=* *= working wavelength grid =* *-----------------------------------------------------------------------------* use datafile_mod, only: datadir implicit none ! input integer :: nw ! number of wavelength grid points real, dimension(nw) :: wl ! lower wavelength for each interval ! output real, dimension(nw) :: xsh2o2 ! h2o2 cross-sections (cm2) ! local real, parameter :: deltax = 1.e-4 integer, parameter :: kdata = 100 real, dimension(kdata) :: x1, y1 real, dimension(nw) :: yg integer :: i, ierr, iw, n, idum integer :: kin, kout ! input/output logical units character*100 fil kin = 10 ! read cross-sections fil = trim(datadir)//'/cross_sections/h2o2_composite.txt' print*, 'section efficace H2O2: ', fil OPEN(kin,FILE=fil,STATUS='OLD') READ(kin,*) idum,n DO i = 1, idum-2 READ(kin,*) ENDDO DO i = 1, n READ(kin,*) x1(i), y1(i) ENDDO CLOSE (kin) CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.) CALL addpnt(x1,y1,kdata,n, 0.,0.) CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.) CALL addpnt(x1,y1,kdata,n, 1.e+38,0.) CALL inter2(nw,wl,yg,n,x1,y1,ierr) IF (ierr .NE. 0) THEN WRITE(*,*) ierr, fil STOP ENDIF DO iw = 1, nw - 1 xsh2o2(iw) = yg(iw) END DO end subroutine rdxsh2o2 !============================================================================== subroutine rdxsho2(nw, wl, yg) *-----------------------------------------------------------------------------* *= PURPOSE: =* *= Read ho2 cross-sections =* *= JPL 2006 recommendation =* *-----------------------------------------------------------------------------* *= PARAMETERS: =* *= NW - INTEGER, number of specified intervals + 1 in working (I)=* *= wavelength grid =* *-----------------------------------------------------------------------------* use datafile_mod, only: datadir IMPLICIT NONE ! input integer :: nw ! number of wavelength grid points real, dimension(nw) :: wl ! lower wavelength for each interval ! output real, dimension(nw) :: yg ! ho2 cross-sections (cm2) ! local real, parameter :: deltax = 1.e-4 integer, parameter :: kdata = 100 real, dimension(kdata) :: x1, y1 integer :: i, n, ierr character*100 fil integer :: kin, kout ! input/output logical units kin = 10 **** cross sections from Sander et al. [2003] fil = trim(datadir)//'/cross_sections/ho2_jpl2003.txt' print*, 'section efficace HO2: ', fil OPEN(kin,FILE=fil,STATUS='OLD') READ(kin,*) n DO i = 1, n READ(kin,*) x1(i), y1(i) ENDDO CLOSE(kin) CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.) CALL addpnt(x1,y1,kdata,n, 0.,0.) CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.) CALL addpnt(x1,y1,kdata,n, 1E38,0.) CALL inter2(nw,wl,yg,n,x1,y1,ierr) IF (ierr .NE. 0) THEN WRITE(*,*) ierr, fil STOP ENDIF end subroutine rdxsho2 !============================================================================== subroutine rdxsh2(nw, wl, wc, yg, yieldh2) *-----------------------------------------------------------------------------* *= PURPOSE: =* *= Read h2 cross-sections and photodissociation yield =* *-----------------------------------------------------------------------------* *= PARAMETERS: =* *= NW - INTEGER, number of specified intervals + 1 in working (I)=* *= wavelength grid =* *-----------------------------------------------------------------------------* use datafile_mod, only: datadir implicit none ! input integer :: nw ! number of wavelength grid points real, dimension(nw) :: wl, wc ! lower and central wavelength for each interval ! output real, dimension(nw) :: yg ! h2 cross-sections (cm2) real, dimension(nw) :: yieldh2 ! photodissociation yield ! local integer, parameter :: kdata = 1000 real, parameter :: deltax = 1.e-4 real, dimension(kdata) :: x1, y1, x2, y2 real :: xl, xu integer :: i, iw, n, ierr integer :: kin, kout ! input/output logical units character*100 fil kin = 10 ! h2 cross sections fil = trim(datadir)//'/cross_sections/h2secef.txt' print*, 'section efficace H2: ', fil OPEN(kin,FILE=fil,STATUS='OLD') n = 792 read(kin,*) ! avoid first line with wavelength = 0. DO i = 1, n READ(kin,*) x1(i), y1(i) ENDDO CLOSE(kin) CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.) CALL addpnt(x1,y1,kdata,n, 0.,0.) CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.) CALL addpnt(x1,y1,kdata,n, 1E38,0.) CALL inter2(nw,wl,yg,n,x1,y1,ierr) IF (ierr .NE. 0) THEN WRITE(*,*) ierr, fil STOP ENDIF ! photodissociation yield fil = trim(datadir)//'/cross_sections/h2_ionef_schunknagy2000.txt' OPEN(UNIT=kin,FILE=fil,STATUS='old') n = 19 read(kin,*) DO i = 1, n READ(kin,*) xl, xu, y2(i) x2(i) = (xl + xu)/2. y2(i) = max(1. - y2(i),0.) END DO CLOSE (kin) CALL addpnt(x2,y2,kdata,n,x2(1)*(1.-deltax),0.) CALL addpnt(x2,y2,kdata,n, 0.,0.) CALL addpnt(x2,y2,kdata,n,x2(n)*(1.+deltax),1.) CALL addpnt(x2,y2,kdata,n, 1.e+38,1.) CALL inter2(nw,wl,yieldh2,n,x2,y2,ierr) IF (ierr .NE. 0) THEN WRITE(*,*) ierr, fil STOP ENDIF end subroutine rdxsh2 !============================================================================== subroutine rdxsno2(nw,wl,xsno2,xsno2_220,xsno2_294, $ yldno2_248, yldno2_298) *-----------------------------------------------------------------------------* *= PURPOSE: =* *= read and grid cross section + quantum yield for NO2 =* *= photolysis =* *= Jenouvrier et al., 1996 200-238 nm *= Vandaele et al., 1998 238-666 nm 220K and 294K *= quantum yield from jpl 2006 *-----------------------------------------------------------------------------* *= PARAMETERS: =* *= NW - INTEGER, number of specified intervals + 1 in working (I)=* *= wavelength grid =* *= WL - REAL, vector of lower limits of wavelength intervals in (I)=* *= working wavelength grid =* *= SQ - REAL, cross section x quantum yield (cm^2) for each (O)=* *= photolysis reaction defined, at each defined wavelength and =* *= at each defined altitude level =* *-----------------------------------------------------------------------------* use datafile_mod, only: datadir implicit none ! input integer :: nw ! number of wavelength grid points real, dimension(nw) :: wl ! lower and central wavelength for each interval ! output real, dimension(nw) :: xsno2, xsno2_220, xsno2_294 ! no2 cross-sections (cm2) real, dimension(nw) :: yldno2_248, yldno2_298 ! quantum yields at 248-298 k ! local integer, parameter :: kdata = 28000 real, parameter :: deltax = 1.e-4 real, dimension(kdata) :: x1, x2, x3, x4, x5, y1, y2, y3, y4, y5 real, dimension(nw) :: yg1, yg2, yg3, yg4, yg5 real :: dum, qy integer :: i, iw, n, n1, n2, n3, n4, n5, ierr character*100 fil integer :: kin, kout ! input/output logical units kin = 10 !*************** NO2 photodissociation * Jenouvrier 1996 + Vandaele 1998 (JPL 2006) fil = trim(datadir)//'/cross_sections/no2_xs_jenouvrier.txt' print*, 'section efficace NO2: ', fil OPEN(UNIT=kin,FILE=fil,status='old') DO i = 1, 3 READ(kin,*) ENDDO n1 = 10001 DO i = 1, n1 READ(kin,*) x1(i), y1(i) end do CALL addpnt(x1,y1,kdata,n1,x1(1)*(1.-deltax), 0.) CALL addpnt(x1,y1,kdata,n1, 0., 0.) CALL addpnt(x1,y1,kdata,n1,x1(n1)*(1.+deltax),0.) CALL addpnt(x1,y1,kdata,n1, 1.e+38, 0.) CALL inter2(nw,wl,yg1,n1,x1,y1,ierr) fil = trim(datadir)//'/cross_sections/no2_xs_vandaele_294K.txt' print*, 'section efficace NO2: ', fil OPEN(UNIT=kin,FILE=fil,status='old') DO i = 1, 3 READ(kin,*) ENDDO n2 = 27993 DO i = 1, n2 READ(kin,*) x2(i), y2(i) end do CALL addpnt(x2,y2,kdata,n2,x2(1)*(1.-deltax), 0.) CALL addpnt(x2,y2,kdata,n2, 0., 0.) CALL addpnt(x2,y2,kdata,n2,x2(n2)*(1.+deltax),0.) CALL addpnt(x2,y2,kdata,n2, 1.e+38, 0.) CALL inter2(nw,wl,yg2,n2,x2,y2,ierr) fil = trim(datadir)//'/cross_sections/no2_xs_vandaele_220K.txt' print*, 'section efficace NO2: ', fil OPEN(UNIT=kin,FILE=fil,status='old') DO i = 1, 3 READ(kin,*) ENDDO n3 = 27993 DO i = 1, n3 READ(kin,*) x3(i), y3(i) end do CALL addpnt(x3,y3,kdata,n3,x3(1)*(1.-deltax), 0.) CALL addpnt(x3,y3,kdata,n3, 0., 0.) CALL addpnt(x3,y3,kdata,n3,x3(n3)*(1.+deltax),0.) CALL addpnt(x3,y3,kdata,n3, 1.e+38, 0.) CALL inter2(nw,wl,yg3,n3,x3,y3,ierr) do iw = 1, nw - 1 xsno2(iw) = yg1(iw) xsno2_294(iw) = yg2(iw) xsno2_220(iw) = yg3(iw) end do ! photodissociation efficiency from jpl 2006 fil = trim(datadir)//'/cross_sections/no2_yield_jpl2006.txt' print*, 'quantum yield NO2: ', fil OPEN(UNIT=kin,FILE=fil,STATUS='old') DO i = 1, 5 READ(kin,*) ENDDO n = 25 n4 = n n5 = n DO i = 1, n READ(kin,*) x4(i), y4(i), y5(i) x5(i) = x4(i) ENDDO CLOSE(kin) CALL addpnt(x4,y4,kdata,n4,x4(1)*(1.-deltax),y4(1)) CALL addpnt(x4,y4,kdata,n4, 0.,y4(1)) CALL addpnt(x4,y4,kdata,n4,x4(n4)*(1.+deltax), 0.) CALL addpnt(x4,y4,kdata,n4, 1.e+38, 0.) CALL inter2(nw,wl,yg4,n4,x4,y4,ierr) IF (ierr .NE. 0) THEN WRITE(*,*) ierr, fil STOP ENDIF CALL addpnt(x5,y5,kdata,n5,x5(1)*(1.-deltax),y5(1)) CALL addpnt(x5,y5,kdata,n5, 0.,y5(1)) CALL addpnt(x5,y5,kdata,n5,x5(n5)*(1.+deltax), 0.) CALL addpnt(x5,y5,kdata,n5, 1.e+38, 0.) CALL inter2(nw,wl,yg5,n5,x5,y5,ierr) IF (ierr .NE. 0) THEN WRITE(*,*) ierr, fil STOP ENDIF do iw = 1, nw - 1 yldno2_298(iw) = yg4(iw) yldno2_248(iw) = yg5(iw) end do end subroutine rdxsno2 !============================================================================== subroutine rdxsno(nw, wl, yg, yieldno) *-----------------------------------------------------------------------------* *= PURPOSE: =* *= Read NO cross-sections and photodissociation efficiency =* *= Lida et al 1986 (provided by Francisco Gonzalez-Galindo) =* *-----------------------------------------------------------------------------* *= PARAMETERS: =* *= NW - INTEGER, number of specified intervals + 1 in working (I)=* *= wavelength grid =* *-----------------------------------------------------------------------------* use datafile_mod, only: datadir implicit none ! input integer :: nw ! number of wavelength grid points real, dimension(nw) :: wl ! lower wavelength for each interval ! output real, dimension(nw) :: yg ! no cross-sections (cm2) real, dimension(nw) :: yieldno ! no photodissociation efficiency ! local integer, parameter :: kdata = 110 real, parameter :: deltax = 1.e-4 real, dimension(kdata) :: x1, y1, x2, y2 integer :: i, iw, n, ierr character*100 fil integer :: kin, kout ! input/output logical units kin = 10 ! no cross-sections fil = trim(datadir)//'/cross_sections/no_xs_francisco.txt' print*, 'section efficace NO: ', fil OPEN(kin,FILE=fil,STATUS='OLD') n = 99 DO i = 1, n READ(kin,*) x1(i), y1(i) ENDDO CLOSE(kin) CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.) CALL addpnt(x1,y1,kdata,n, 0.,0.) CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.) CALL addpnt(x1,y1,kdata,n, 1E38,0.) CALL inter2(nw,wl,yg,n,x1,y1,ierr) IF (ierr .NE. 0) THEN WRITE(*,*) ierr, fil STOP ENDIF ! photodissociation yield fil = trim(datadir)//'/cross_sections/noefdis.txt' OPEN(UNIT=kin,FILE=fil,STATUS='old') n = 33 DO i = 1, n READ(kin,*) x2(n-i+1), y2(n-i+1) END DO CLOSE (kin) CALL addpnt(x2,y2,kdata,n,x2(1)*(1.-deltax),0.) CALL addpnt(x2,y2,kdata,n, 0.,0.) CALL addpnt(x2,y2,kdata,n,x2(n)*(1.+deltax),1.) CALL addpnt(x2,y2,kdata,n, 1.e+38,1.) CALL inter2(nw,wl,yieldno,n,x2,y2,ierr) IF (ierr .NE. 0) THEN WRITE(*,*) ierr, fil STOP ENDIF end subroutine rdxsno !============================================================================== subroutine rdxsn2(nw, wl, yg, yieldn2) *-----------------------------------------------------------------------------* *= PURPOSE: =* *= Read n2 cross-sections and photodissociation yield =* *-----------------------------------------------------------------------------* *= PARAMETERS: =* *= NW - INTEGER, number of specified intervals + 1 in working (I)=* *= wavelength grid =* *-----------------------------------------------------------------------------* use datafile_mod, only: datadir implicit none ! input integer :: nw ! number of wavelength grid points real, dimension(nw) :: wl ! lower wavelength for each interval ! output real, dimension(nw) :: yg ! n2 cross-sections (cm2) real, dimension(nw) :: yieldn2 ! n2 photodissociation yield ! local integer, parameter :: kdata = 1100 real, parameter :: deltax = 1.e-4 real, dimension(kdata) :: x1, y1, x2, y2 real :: xl, xu integer :: i, iw, n, ierr integer :: kin, kout ! input/output logical units character*100 fil kin = 10 ! n2 cross sections fil = trim(datadir)//'/cross_sections/n2secef_01nm.txt' print*, 'section efficace N2: ', fil OPEN(kin,FILE=fil,STATUS='OLD') n = 1020 DO i = 1, n READ(kin,*) x1(i), y1(i) ENDDO CLOSE(kin) CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.) CALL addpnt(x1,y1,kdata,n, 0.,0.) CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.) CALL addpnt(x1,y1,kdata,n, 1E38,0.) CALL inter2(nw,wl,yg,n,x1,y1,ierr) IF (ierr .NE. 0) THEN WRITE(*,*) ierr, fil STOP ENDIF ! photodissociation yield fil = trim(datadir)//'/cross_sections/n2_ionef_schunknagy2000.txt' OPEN(UNIT=kin,FILE=fil,STATUS='old') n = 19 read(kin,*) DO i = 1, n READ(kin,*) xl, xu, y2(i) x2(i) = (xl + xu)/2. y2(i) = 1. - y2(i) END DO CLOSE (kin) CALL addpnt(x2,y2,kdata,n,x2(1)*(1.-deltax),0.) CALL addpnt(x2,y2,kdata,n, 0.,0.) CALL addpnt(x2,y2,kdata,n,x2(n)*(1.+deltax),1.) CALL addpnt(x2,y2,kdata,n, 1.e+38,1.) CALL inter2(nw,wl,yieldn2,n,x2,y2,ierr) IF (ierr .NE. 0) THEN WRITE(*,*) ierr, fil STOP ENDIF end subroutine rdxsn2 !============================================================================== subroutine setalb(nw,wl,albedo) *-----------------------------------------------------------------------------* *= PURPOSE: =* *= Set the albedo of the surface. The albedo is assumed to be Lambertian, =* *= i.e., the reflected light is isotropic, and idependt of the direction =* *= of incidence of light. Albedo can be chosen to be wavelength dependent. =* *-----------------------------------------------------------------------------* *= PARAMETERS: =* *= NW - INTEGER, number of specified intervals + 1 in working (I)=* *= wavelength grid =* *= WL - REAL, vector of lower limits of wavelength intervals in (I)=* *= working wavelength grid =* *= ALBEDO - REAL, surface albedo at each specified wavelength (O)=* *-----------------------------------------------------------------------------* implicit none ! input: (wavelength working grid data) INTEGER nw REAL wl(nw) ! output: REAL albedo(nw) ! local: INTEGER iw REAL alb ! 0.015: mean value from clancy et al., icarus, 49-63, 1999. alb = 0.015 do iw = 1, nw - 1 albedo(iw) = alb end do return end !============================================================================== subroutine setair(nlev, nw, wl, wc, press, temp, zmmean, $ colinc, dtrl) *-----------------------------------------------------------------------------* *= PURPOSE: =* *= computes air column increments and rayleigh optical depth =* *-----------------------------------------------------------------------------* implicit none ! input: integer :: nlev, nw real, dimension(nw) :: wl, wc ! lower and central wavelength grid (nm) real, dimension(nlev) :: press, temp, zmmean ! pressure (hpa), temperature (k), molecular mass (g.mol-1) ! output: real, dimension(nlev) :: colinc ! air column increments (molecule.cm-2) real, dimension(nlev,nw) :: dtrl ! rayleigh optical depth ! local: real, parameter :: avo = 6.022e23 real, parameter :: g = 3.72 real :: dp real :: alpha, rho, df, pi real, dimension(nw) :: srayl integer :: ilev, iw c compute column increments do ilev = 1, nlev-1 dp = (press(ilev) - press(ilev+1))*100. colinc(ilev) = avo*0.1*dp/(zmmean(ilev)*g) end do colinc(nlev) = avo*0.1*press(nlev)*100./(zmmean(nlev)*g) ! do ilev = nlev,1,-1 ! print*, ilev, press(ilev), temp(ilev), zmmean(ilev), ! $ colinc(ilev) ! end do do iw = 1, nw - 1 c calcul de section efficace Rayleigh: voir Atreya and Gu, JGR, c 13133-13145, 1994. c c rho : normal depolarization ratio = 0.0774 for CO2 c df : depolarization factor c alpha : average dielectric polarizability = 2.911e-24 cm3 for CO2 pi = acos(-1.0) rho = 0.0774 alpha = 2.911e-24*1.e21 ! nm3 df = (6. + 3.*rho)/(6. - 7.*rho) srayl(iw) = 8.*pi/3.*((2.*pi/wc(iw))**4)*(alpha**2)*df ! nm2 do ilev = 1, nlev dtrl(ilev,iw) = colinc(ilev)*srayl(iw)*1.e-14 ! cm2 end do end do return end !============================================================================== subroutine setco2(nd, nlayer, nw, wc, tlay, xsco2_195, xsco2_295, $ xsco2_370, yieldco2, j_co2_o, j_co2_o1d, $ colinc, rm, dt, sj) *-----------------------------------------------------------------------------* *= PURPOSE: =* *= Set up the CO2 temperature-dependent cross-sections and optical depth =* *-----------------------------------------------------------------------------* implicit none ! input: integer :: nd ! number of photolysis rates integer :: nlayer ! number of vertical layers integer :: nw ! number of wavelength grid points integer :: j_co2_o, j_co2_o1d ! photolysis indexes real, dimension(nw) :: wc ! central wavelength for each interval real, dimension(nw) :: xsco2_195, xsco2_295, xsco2_370 ! co2 cross-sections (cm2) real, dimension(nw) :: yieldco2 ! co2 photodissociation yield real, dimension(nlayer) :: tlay ! temperature (k) real, dimension(nlayer) :: rm ! co2 mixing ratio real, dimension(nlayer) :: colinc ! air column increment (molecule.cm-2) ! output: real, dimension(nlayer,nw) :: dt ! optical depth real, dimension(nlayer,nw,nd) :: sj ! cross-section array (cm2) ! local: integer :: extrapol integer :: i, l real :: temp, sco2 ! extrapol = 0 no extrapolation below 195 k ! extrapol = 1 extrapolation below 195 k extrapol = 0 do i = 1, nlayer if (extrapol == 1) then temp = tlay(i) else temp = max(tlay(i), 195.) end if temp = min(temp, 370.) do l = 1, nw-1 if (temp <= 295.) then if (xsco2_195(l) /= 0. .and. xsco2_295(l) /= 0.) then sco2 = alog(xsco2_195(l)) $ + (alog(xsco2_295(l)) - alog(xsco2_195(l))) $ /(295. - 195.)*(temp - 195.) sco2 = exp(sco2) else sco2 = 0. end if else if (xsco2_295(l) /= 0. .and. xsco2_370(l) /= 0.) then sco2 = alog(xsco2_295(l)) $ + (alog(xsco2_370(l)) - alog(xsco2_295(l))) $ /(370. - 295.)*(temp - 295.) sco2 = exp(sco2) else sco2 = 0. end if end if ! optical depth dt(i,l) = colinc(i)*rm(i)*sco2 ! production of o(1d) for wavelengths shorter than 167 nm if (wc(l) >= 167.) then sj(i,l,j_co2_o) = sco2*yieldco2(l) sj(i,l,j_co2_o1d) = 0. else sj(i,l,j_co2_o) = 0. sj(i,l,j_co2_o1d) = sco2*yieldco2(l) end if end do end do ! do l = 1,nw-1 ! write(35,*) wc(l), sj(1,l,j_co2_o1d), ! $ sj(1,l,j_co2_o), tlay(1) ! end do end subroutine setco2 !============================================================================== subroutine seto2(nd, nlayer, nw, wc, mopt, tlay, xso2_150, $ xso2_200, xso2_250, xso2_300, yieldo2, j_o2_o, $ j_o2_o1d, colinc, rm, dt, sj) *-----------------------------------------------------------------------------* *= PURPOSE: =* *= Set up the O2 temperature-dependent cross-sections and optical depth =* *-----------------------------------------------------------------------------* implicit none ! input: integer :: nd ! number of photolysis rates integer :: nlayer ! number of vertical layers integer :: nw ! number of wavelength grid points integer :: mopt ! high-res/low-res switch integer :: j_o2_o, j_o2_o1d ! photolysis indexes real, dimension(nw) :: wc ! central wavelength for each interval real, dimension(nw) :: xso2_150, xso2_200, xso2_250, ! o2 cross-sections (cm2) $ xso2_300 real, dimension(nw) :: yieldo2 ! o2 photodissociation yield real, dimension(nlayer) :: tlay ! temperature (k) real, dimension(nlayer) :: rm ! o2 mixing ratio real, dimension(nlayer) :: colinc ! air column increment (molecule.cm-2) ! output: real, dimension(nlayer,nw) :: dt ! optical depth real, dimension(nlayer,nw,nd) :: sj ! cross-section array (cm2) ! local: integer :: ilev, iw real :: temp real :: xso2, factor ! correction by factor if low-resolution in schumann-runge bands if (mopt == 1) then factor = 1. else if (mopt == 2) then factor = 0.8 end if ! calculate temperature dependance do ilev = 1,nlayer temp = max(tlay(ilev),150.) temp = min(temp, 300.) do iw = 1, nw-1 if (tlay(ilev) > 250.) then xso2 = xso2_250(iw) + (xso2_300(iw) - xso2_250(iw)) $ /(300. - 250.)*(temp - 250.) else if (tlay(ilev) > 200.) then xso2 = xso2_200(iw) + (xso2_250(iw) - xso2_200(iw)) $ /(250. - 200.)*(temp - 200.) else xso2 = xso2_150(iw) + (xso2_200(iw) - xso2_150(iw)) $ /(200. - 150.)*(temp - 150.) end if if (wc(iw) > 180. .and. wc(iw) < 200.) then xso2 = xso2*factor end if ! optical depth dt(ilev,iw) = colinc(ilev)*rm(ilev)*xso2 ! production of o(1d) for wavelengths shorter than 175 nm if (wc(iw) >= 175.) then sj(ilev,iw,j_o2_o) = xso2*yieldo2(iw) sj(ilev,iw,j_o2_o1d) = 0. else sj(ilev,iw,j_o2_o) = 0. sj(ilev,iw,j_o2_o1d) = xso2*yieldo2(iw) end if end do end do end subroutine seto2 !============================================================================== subroutine seto3(nd, nlayer, nw, wc, tlay, xso3_218, xso3_298, $ j_o3_o, j_o3_o1d, $ colinc, rm, dt, sj) *-----------------------------------------------------------------------------* *= PURPOSE: =* *= Set up the O3 temperature dependent cross-sections and optical depth =* *-----------------------------------------------------------------------------* implicit none ! input: integer :: nd ! number of photolysis rates integer :: nlayer ! number of vertical layers integer :: nw ! number of wavelength grid points integer :: j_o3_o, j_o3_o1d ! photolysis indexes real, dimension(nw) :: wc ! central wavelength for each interval real, dimension(nw) :: xso3_218, xso3_298 ! o3 cross-sections (cm2) real, dimension(nlayer) :: tlay ! temperature (k) real, dimension(nlayer) :: rm ! o3 mixing ratio real, dimension(nlayer) :: colinc ! air column increment (molecule.cm-2) ! output: real, dimension(nlayer,nw) :: dt ! optical depth real, dimension(nlayer,nw,nd) :: sj ! cross-section array (cm2) ! local: ! integer :: ilev, iw real :: temp real, dimension(nw) :: xso3(nw) real, dimension(nw) :: qy1d ! quantum yield for o(1d) production real :: q1, q2, a1, a2, a3 do ilev = 1, nlayer temp = max(tlay(ilev), 218.) temp = min(temp,298.) do iw = 1, nw-1 xso3(iw) = xso3_218(iw) + (xso3_298(iw) - xso3_218(iw)) $ /(298. - 218.) *(temp - 218.) ! optical depth dt(ilev,iw) = colinc(ilev)*rm(ilev)*xso3(iw) end do ! calculate quantum yield for o(1d) production (jpl 2006) temp = max(tlay(ilev),200.) temp = min(temp,320.) do iw = 1, nw-1 if (wc(iw) <= 306.) then qy1d(iw) = 0.90 else if (wc(iw) > 306. .and. wc(iw) < 328.) then q1 = 1. q2 = exp(-825.518/(0.695*temp)) a1 = (304.225 - wc(iw))/5.576 a2 = (314.957 - wc(iw))/6.601 a3 = (310.737 - wc(iw))/2.187 qy1d(iw) = (q1/(q1 + q2))*0.8036*exp(-(a1*a1*a1*a1)) $ + (q2/(q1 + q2))*8.9061*(temp/300.)**2. $ *exp(-(a2*a2)) $ + 0.1192*(temp/300.)**1.5*exp(-(a3*a3)) $ + 0.0765 else if (wc(iw) >= 328. .and. wc(iw) <= 340.) then qy1d(iw) = 0.08 else qy1d(iw) = 0. endif end do do iw = 1, nw-1 sj(ilev,iw,j_o3_o) = xso3(iw)*(1. - qy1d(iw)) sj(ilev,iw,j_o3_o1d) = xso3(iw)*qy1d(iw) end do end do end subroutine seto3 !============================================================================== subroutine seth2o2(nd, nlayer, nw, wc, tlay, xsh2o2, j_h2o2, $ colinc, rm, dt, sj) *-----------------------------------------------------------------------------* *= PURPOSE: =* *= Set up the h2o2 temperature dependent cross-sections and optical depth =* *-----------------------------------------------------------------------------* implicit none ! input: integer :: nd ! number of photolysis rates integer :: nlayer ! number of vertical layers integer :: nw ! number of wavelength grid points integer :: j_h2o2 ! photolysis index real, dimension(nw) :: wc ! central wavelength for each interval real, dimension(nw) :: xsh2o2 ! h2o2 cross-sections (cm2) real, dimension(nlayer) :: tlay ! temperature (k) real, dimension(nlayer) :: rm ! h2o2 mixing ratio real, dimension(nlayer) :: colinc ! air column increment (molecule.cm-2) ! output: real, dimension(nlayer,nw) :: dt ! optical depth real, dimension(nlayer,nw,nd) :: sj ! cross-section array (cm2) ! local: integer :: ilev, iw real :: a0, a1, a2, a3, a4, a5, a6, a7 real :: b0, b1, b2, b3, b4 real :: lambda, suma, sumb, chi, temp, xs A0 = 6.4761E+04 A1 = -9.2170972E+02 A2 = 4.535649 A3 = -4.4589016E-03 A4 = -4.035101E-05 A5 = 1.6878206E-07 A6 = -2.652014E-10 A7 = 1.5534675E-13 B0 = 6.8123E+03 B1 = -5.1351E+01 B2 = 1.1522E-01 B3 = -3.0493E-05 B4 = -1.0924E-07 ! temperature dependance: jpl 2006 do ilev = 1,nlayer temp = min(max(tlay(ilev),200.),400.) chi = 1./(1. + exp(-1265./temp)) do iw = 1, nw-1 if ((wc(iw) >= 260.) .and. (wc(iw) < 350.)) then lambda = wc(iw) sumA = ((((((A7*lambda + A6)*lambda + A5)*lambda + $ A4)*lambda +A3)*lambda + A2)*lambda + $ A1)*lambda + A0 sumB = (((B4*lambda + B3)*lambda + B2)*lambda + $ B1)*lambda + B0 xs = (chi*sumA + (1. - chi)*sumB)*1.e-21 sj(ilev,iw,j_h2o2) = xs else sj(ilev,iw,j_h2o2) = xsh2o2(iw) end if ! optical depth dt(ilev,iw) = colinc(ilev)*rm(ilev)*sj(ilev,iw,j_h2o2) end do end do end subroutine seth2o2 !============================================================================== subroutine setno2(nd, nlayer, nw, wc, tlay, xsno2, xsno2_220, $ xsno2_294, yldno2_248, yldno2_298, j_no2, $ colinc, rm, dt, sj) *-----------------------------------------------------------------------------* *= PURPOSE: =* *= Set up the no2 temperature-dependent cross-sections and optical depth =* *-----------------------------------------------------------------------------* implicit none ! input: integer :: nd ! number of photolysis rates integer :: nlayer ! number of vertical layers integer :: nw ! number of wavelength grid points integer :: j_no2 ! photolysis index real, dimension(nw) :: wc ! central wavelength for each interval real, dimension(nw) :: xsno2, xsno2_220, xsno2_294 ! no2 absorption cross-section at 220-294 k (cm2) real, dimension(nw) :: yldno2_248, yldno2_298 ! no2 quantum yield at 248-298 k real, dimension(nlayer) :: tlay ! temperature (k) real, dimension(nlayer) :: rm ! no2 mixing ratio real, dimension(nlayer) :: colinc ! air column increment (molecule.cm-2) ! output: real, dimension(nlayer,nw) :: dt ! optical depth real, dimension(nlayer,nw,nd) :: sj ! cross-section array (cm2) ! local: integer :: ilev, iw real :: temp, qy ! temperature dependance: jpl 2006 do ilev = 1,nlayer temp = max(220.,min(tlay(ilev),294.)) do iw = 1, nw - 1 if (wc(iw) < 238.) then sj(ilev,iw,j_no2) = xsno2(iw) else sj(ilev,iw,j_no2) = xsno2_220(iw) $ + (xsno2_294(iw) - xsno2_220(iw)) $ /(294. - 220.)*(temp - 220.) end if ! optical depth dt(ilev,iw) = colinc(ilev)*rm(ilev)*sj(ilev,iw,j_no2) end do end do ! quantum yield: jpl 2006 do ilev = 1,nlayer temp = max(248.,min(tlay(ilev),298.)) do iw = 1, nw - 1 qy = yldno2_248(iw) + (yldno2_298(iw) - yldno2_248(iw)) $ /(298. - 248.)*(temp - 248.) sj(ilev,iw,j_no2) = sj(ilev,iw,j_no2)*qy end do end do ! do iw = 1, nw-1 ! write(35,*) wc(iw), sj(1,iw,j_no2), tlay(1) ! end do ! stop end subroutine setno2 !============================================================================== subroutine setaer(nlayer,alt,tau,nw,dtaer,omaer,gaer) *-----------------------------------------------------------------------------* *= PURPOSE: =* *= Set aerosol properties for each specified altitude layer. Properties =* *= may be wavelength dependent. =* *-----------------------------------------------------------------------------* implicit none ! input integer :: nlayer ! number of vertical layers integer :: nw ! number of wavelength grid points real, dimension(nlayer) :: alt ! altitude (km) real :: tau ! integrated aerosol optical depth at the surface ! output real, dimension(nlayer,nw) :: dtaer ! aerosol optical depth real, dimension(nlayer,nw) :: omaer ! aerosol single scattering albedo real, dimension(nlayer,nw) :: gaer ! aerosol asymmetry parameter ! local integer :: ilay, iw real, dimension(nlayer) :: aer ! dust extinction real :: omega, g, scaleh, gamma real :: dz, tautot, q0 omega = 0.622 ! single scattering albedo : wolff et al.(2010) at 258 nm g = 0.88 ! asymmetry factor : mateshvili et al. (2007) at 210 nm scaleh = 10. ! scale height (km) gamma = 0.03 ! conrath parameter dtaer(:,:) = 0. omaer(:,:) = 0. gaer(:,:) = 0. ! optical depth profile: tautot = 0. do ilay = 1, nlayer-1 dz = alt(ilay+1) - alt(ilay) tautot = tautot + exp(gamma*(1. - exp(alt(ilay)/scaleh)))*dz end do ! print*, 'tau = ', tau ! print*, 'q0 = ', q0 q0 = tau/tautot do ilay = 1, nlayer-1 dz = alt(ilay+1) - alt(ilay) dtaer(ilay,:) = q0*exp(gamma*(1. - exp(alt(ilay)/scaleh)))*dz omaer(ilay,:) = omega gaer(ilay,:) = g end do ! do ilay = nlayer,1,-1 ! print*, alt(ilay), q0*exp(gamma*(1. - exp(alt(ilay)/scaleh))), ! $ dtaer(ilay,1) ! end do end subroutine setaer !============================================================================== subroutine setcld(nlayer,nw,dtcld,omcld,gcld) *-----------------------------------------------------------------------------* *= PURPOSE: =* *= Set cloud properties for each specified altitude layer. Properties =* *= may be wavelength dependent. =* *-----------------------------------------------------------------------------* implicit none ! input integer :: nlayer ! number of vertical layers integer :: nw ! number of wavelength grid points ! output real, dimension(nlayer,nw) :: dtcld ! cloud optical depth real, dimension(nlayer,nw) :: omcld ! cloud single scattering albedo real, dimension(nlayer,nw) :: gcld ! cloud asymmetry parameter ! local integer :: ilay, iw c dtcld : optical depth c omcld : single scattering albedo c gcld : asymmetry factor do ilay = 1, nlayer - 1 do iw = 1, nw - 1 dtcld(ilay,iw) = 0. ! no clouds for the moment omcld(ilay,iw) = 0.99 gcld(ilay,iw) = 0.85 end do end do end !============================================================================== subroutine sphers(nlev, z, zen, dsdh, nid) *-----------------------------------------------------------------------------* *= PURPOSE: =* *= Calculate slant path over vertical depth ds/dh in spherical geometry. =* *= Calculation is based on: A.Dahlback, and K.Stamnes, A new spheric model =* *= for computing the radiation field available for photolysis and heating =* *= at twilight, Planet.Space Sci., v39, n5, pp. 671-683, 1991 (Appendix B) =* *-----------------------------------------------------------------------------* *= PARAMETERS: =* *= NZ - INTEGER, number of specified altitude levels in the working (I)=* *= grid =* *= Z - REAL, specified altitude working grid (km) (I)=* *= ZEN - REAL, solar zenith angle (degrees) (I)=* *= DSDH - REAL, slant path of direct beam through each layer crossed (O)=* *= when travelling from the top of the atmosphere to layer i; =* *= DSDH(i,j), i = 0..NZ-1, j = 1..NZ-1 =* *= NID - INTEGER, number of layers crossed by the direct beam when (O)=* *= travelling from the top of the atmosphere to layer i; =* *= NID(i), i = 0..NZ-1 =* *-----------------------------------------------------------------------------* c implicit none * input INTEGER nlev REAL zen, z(nlev) * output INTEGER nid(0:nlev) REAL dsdh(0:nlev,nlev) * more program constants REAL re, ze(nlev) REAL dr real radius parameter (radius = 3393.) * local REAL pi, zenrad, rpsinz, rj, rjp1, dsj, dhj, ga, gb, sm INTEGER i, j, k INTEGER id INTEGER nlay REAL zd(0:nlev-1) *----------------------------------------------------------------------------- pi = acos(-1.0) dr = pi/180. zenrad = zen*dr * number of layers: nlay = nlev - 1 * include the elevation above sea level to the radius of Mars: re = radius + z(1) * correspondingly z changed to the elevation above Mars surface: DO k = 1, nlev ze(k) = z(k) - z(1) END DO * inverse coordinate of z zd(0) = ze(nlev) DO k = 1, nlay zd(k) = ze(nlev - k) END DO * initialize dsdh(i,j), nid(i) DO i = 0, nlev nid(i) = 0 DO j = 1, nlev dsdh(i,j) = 0. END DO END DO * calculate ds/dh of every layer DO 100 i = 0, nlay rpsinz = (re + zd(i)) * SIN(zenrad) IF ( (zen .GT. 90.0) .AND. (rpsinz .LT. re) ) THEN nid(i) = -1 ELSE * * Find index of layer in which the screening height lies * id = i IF( zen .GT. 90.0 ) THEN DO 10 j = 1, nlay IF( (rpsinz .LT. ( zd(j-1) + re ) ) .AND. $ (rpsinz .GE. ( zd(j) + re )) ) id = j 10 CONTINUE END IF DO 20 j = 1, id sm = 1.0 IF(j .EQ. id .AND. id .EQ. i .AND. zen .GT. 90.0) $ sm = -1.0 rj = re + zd(j-1) rjp1 = re + zd(j) dhj = zd(j-1) - zd(j) ga = rj*rj - rpsinz*rpsinz gb = rjp1*rjp1 - rpsinz*rpsinz IF (ga .LT. 0.0) ga = 0.0 IF (gb .LT. 0.0) gb = 0.0 IF(id.GT.i .AND. j.EQ.id) THEN dsj = SQRT( ga ) ELSE dsj = SQRT( ga ) - sm*SQRT( gb ) END IF dsdh(i,j) = dsj / dhj 20 CONTINUE nid(i) = id END IF 100 CONTINUE c RETURN END !============================================================================== SUBROUTINE rtlink(nlev, nw, iw, ag, zen, dsdh, nid, dtrl, $ dagas, dtcld, omcld, gcld, dtaer, omaer, gaer, $ edir, edn, eup, fdir, fdn, fup) implicit none * input integer :: nw ! number of wavelength grid points INTEGER nlev, iw REAL ag REAL zen REAL dsdh(0:nlev,nlev) INTEGER nid(0:nlev) REAL dtrl(nlev,nw) REAL dagas(nlev,nw) REAL dtcld(nlev,nw), omcld(nlev,nw), gcld(nlev,nw) REAL dtaer(nlev,nw), omaer(nlev,nw), gaer(nlev,nw) * output REAL edir(nlev), edn(nlev), eup(nlev) REAL fdir(nlev), fdn(nlev), fup(nlev) * local: REAL dt(nlev), om(nlev), g(nlev) REAL dtabs,dtsct,dscld,dsaer,dacld,daaer INTEGER i, ii real, parameter :: largest = 1.e+36 * specific two ps2str REAL ediri(nlev), edni(nlev), eupi(nlev) REAL fdiri(nlev), fdni(nlev), fupi(nlev) logical, save :: delta = .true. *_______________________________________________________________________ * initialize: do i = 1, nlev fdir(i) = 0. fup(i) = 0. fdn(i) = 0. edir(i) = 0. eup(i) = 0. edn(i) = 0. end do c do i = 1, nlev - 1 dscld = dtcld(i,iw)*omcld(i,iw) dacld = dtcld(i,iw)*(1.-omcld(i,iw)) c dsaer = dtaer(i,iw)*omaer(i,iw) daaer = dtaer(i,iw)*(1.-omaer(i,iw)) c dtsct = dtrl(i,iw) + dscld + dsaer dtabs = dagas(i,iw) + dacld + daaer c dtabs = amax1(dtabs,1./largest) dtsct = amax1(dtsct,1./largest) c c invert z-coordinate: c ii = nlev - i dt(ii) = dtsct + dtabs om(ii) = dtsct/(dtsct + dtabs) IF(dtsct .EQ. 1./largest) om(ii) = 1./largest g(ii) = (gcld(i,iw)*dscld + $ gaer(i,iw)*dsaer)/dtsct end do c c call rt routine: c call ps2str(nlev, zen, ag, dt, om, g, $ dsdh, nid, delta, $ fdiri, fupi, fdni, ediri, eupi, edni) c c output (invert z-coordinate) c do i = 1, nlev ii = nlev - i + 1 fdir(i) = fdiri(ii) fup(i) = fupi(ii) fdn(i) = fdni(ii) edir(i) = ediri(ii) eup(i) = eupi(ii) edn(i) = edni(ii) end do c return end *=============================================================================* SUBROUTINE ps2str(nlev,zen,rsfc,tauu,omu,gu, $ dsdh, nid, delta, $ fdr, fup, fdn, edr, eup, edn) *-----------------------------------------------------------------------------* *= PURPOSE: =* *= Solve two-stream equations for multiple layers. The subroutine is based =* *= on equations from: Toon et al., J.Geophys.Res., v94 (D13), Nov 20, 1989.=* *= It contains 9 two-stream methods to choose from. A pseudo-spherical =* *= correction has also been added. =* *-----------------------------------------------------------------------------* *= PARAMETERS: =* *= NLEVEL - INTEGER, number of specified altitude levels in the working (I)=* *= grid =* *= ZEN - REAL, solar zenith angle (degrees) (I)=* *= RSFC - REAL, surface albedo at current wavelength (I)=* *= TAUU - REAL, unscaled optical depth of each layer (I)=* *= OMU - REAL, unscaled single scattering albedo of each layer (I)=* *= GU - REAL, unscaled asymmetry parameter of each layer (I)=* *= DSDH - REAL, slant path of direct beam through each layer crossed (I)=* *= when travelling from the top of the atmosphere to layer i; =* *= DSDH(i,j), i = 0..NZ-1, j = 1..NZ-1 =* *= NID - INTEGER, number of layers crossed by the direct beam when (I)=* *= travelling from the top of the atmosphere to layer i; =* *= NID(i), i = 0..NZ-1 =* *= DELTA - LOGICAL, switch to use delta-scaling (I)=* *= .TRUE. -> apply delta-scaling =* *= .FALSE.-> do not apply delta-scaling =* *= FDR - REAL, contribution of the direct component to the total (O)=* *= actinic flux at each altitude level =* *= FUP - REAL, contribution of the diffuse upwelling component to (O)=* *= the total actinic flux at each altitude level =* *= FDN - REAL, contribution of the diffuse downwelling component to (O)=* *= the total actinic flux at each altitude level =* *= EDR - REAL, contribution of the direct component to the total (O)=* *= spectral irradiance at each altitude level =* *= EUP - REAL, contribution of the diffuse upwelling component to (O)=* *= the total spectral irradiance at each altitude level =* *= EDN - REAL, contribution of the diffuse downwelling component to (O)=* *= the total spectral irradiance at each altitude level =* *-----------------------------------------------------------------------------* c implicit none c ******* * input: ******* INTEGER nlev REAL zen, rsfc REAL tauu(nlev), omu(nlev), gu(nlev) REAL dsdh(0:nlev,nlev) INTEGER nid(0:nlev) LOGICAL delta ******* * output: ******* REAL fup(nlev),fdn(nlev),fdr(nlev) REAL eup(nlev),edn(nlev),edr(nlev) ******* * local: ******* REAL tausla(0:nlev), tauc(0:nlev) REAL mu2(0:nlev), mu, sum * internal coefficients and matrix REAL lam(nlev),taun(nlev),bgam(nlev) REAL e1(nlev),e2(nlev),e3(nlev),e4(nlev) REAL cup(nlev),cdn(nlev),cuptn(nlev),cdntn(nlev) REAL mu1(nlev) INTEGER row REAL a(2*nlev),b(2*nlev),d(2*nlev),e(2*nlev),y(2*nlev) ******* * other: ******* REAL pifs, fdn0 REAL gi(nlev), omi(nlev), tempg REAL f, g, om REAL gam1, gam2, gam3, gam4 real, parameter :: largest = 1.e+36 real, parameter :: precis = 1.e-7 * For calculations of Associated Legendre Polynomials for GAMA1,2,3,4 * in delta-function, modified quadrature, hemispheric constant, * Hybrid modified Eddington-delta function metods, p633,Table1. * W.E.Meador and W.R.Weaver, GAS,1980,v37,p.630 * W.J.Wiscombe and G.W. Grams, GAS,1976,v33,p2440, * uncomment the following two lines and the appropriate statements further * down. C REAL YLM0, YLM2, YLM4, YLM6, YLM8, YLM10, YLM12, YLMS, BETA0, C > BETA1, BETAn, amu1, subd REAL expon, expon0, expon1, divisr, temp, up, dn REAL ssfc INTEGER nlayer, mrows, lev INTEGER i, j * Some additional program constants: real pi, dr REAL eps PARAMETER (eps = 1.E-3) *_______________________________________________________________________ * MU = cosine of solar zenith angle * RSFC = surface albedo * TAUU = unscaled optical depth of each layer * OMU = unscaled single scattering albedo * GU = unscaled asymmetry factor * KLEV = max dimension of number of layers in atmosphere * NLAYER = number of layers in the atmosphere * NLEVEL = nlayer + 1 = number of levels * initial conditions: pi*solar flux = 1; diffuse incidence = 0 pifs = 1. fdn0 = 0. nlayer = nlev - 1 pi = acos(-1.) dr = pi/180. mu = COS(zen*dr) ************** compute coefficients for each layer: * GAM1 - GAM4 = 2-stream coefficients, different for different approximations * EXPON0 = calculation of e when TAU is zero * EXPON1 = calculation of e when TAU is TAUN * CUP and CDN = calculation when TAU is zero * CUPTN and CDNTN = calc. when TAU is TAUN * DIVISR = prevents division by zero do j = 0, nlev tauc(j) = 0. tausla(j) = 0. mu2(j) = 1./SQRT(largest) end do c IF (.NOT. delta) THEN DO i = 1, nlayer gi(i) = gu(i) omi(i) = omu(i) taun(i) = tauu(i) END DO ELSE * delta-scaling. Have to be done for delta-Eddington approximation, * delta discrete ordinate, Practical Improved Flux Method, delta function, * and Hybrid modified Eddington-delta function methods approximations DO i = 1, nlayer f = gu(i)*gu(i) gi(i) = (gu(i) - f)/(1 - f) omi(i) = (1 - f)*omu(i)/(1 - omu(i)*f) taun(i) = (1 - omu(i)*f)*tauu(i) END DO END IF * * calculate slant optical depth at the top of the atmosphere when zen>90. * in this case, higher altitude of the top layer is recommended which can * be easily changed in gridz.f. * IF (zen .GT. 90.0) THEN IF (nid(0) .LT. 0) THEN tausla(0) = largest ELSE sum = 0.0 DO j = 1, nid(0) sum = sum + 2.*taun(j)*dsdh(0,j) END DO tausla(0) = sum END IF END IF * DO 11, i = 1, nlayer g = gi(i) om = omi(i) tauc(i) = tauc(i-1) + taun(i) * stay away from 1 by precision. For g, also stay away from -1 tempg = AMIN1(abs(g),1. - precis) g = SIGN(tempg,g) om = AMIN1(om,1.-precis) * calculate slant optical depth * IF (nid(i) .LT. 0) THEN tausla(i) = largest ELSE sum = 0.0 DO j = 1, MIN(nid(i),i) sum = sum + taun(j)*dsdh(i,j) END DO DO j = MIN(nid(i),i)+1,nid(i) sum = sum + 2.*taun(j)*dsdh(i,j) END DO tausla(i) = sum IF (tausla(i) .EQ. tausla(i-1)) THEN mu2(i) = SQRT(largest) ELSE mu2(i) = (tauc(i)-tauc(i-1))/(tausla(i)-tausla(i-1)) mu2(i) = SIGN( AMAX1(ABS(mu2(i)),1./SQRT(largest)), $ mu2(i) ) END IF END IF * *** the following gamma equations are from pg 16,289, Table 1 *** save mu1 for each approx. for use in converting irradiance to actinic flux * Eddington approximation(Joseph et al., 1976, JAS, 33, 2452): c gam1 = (7. - om*(4. + 3.*g))/4. c gam2 = -(1. - om*(4. - 3.*g))/4. c gam3 = (2. - 3.*g*mu)/4. c gam4 = 1. - gam3 c mu1(i) = 0.5 * quadrature (Liou, 1973, JAS, 30, 1303-1326; 1974, JAS, 31, 1473-1475): c gam1 = 1.7320508*(2. - om*(1. + g))/2. c gam2 = 1.7320508*om*(1. - g)/2. c gam3 = (1. - 1.7320508*g*mu)/2. c gam4 = 1. - gam3 c mu1(i) = 1./sqrt(3.) * hemispheric mean (Toon et al., 1089, JGR, 94, 16287): gam1 = 2. - om*(1. + g) gam2 = om*(1. - g) gam3 = (2. - g*mu)/4. gam4 = 1. - gam3 mu1(i) = 0.5 * PIFM (Zdunkovski et al.,1980, Conrib.Atmos.Phys., 53, 147-166): c GAM1 = 0.25*(8. - OM*(5. + 3.*G)) c GAM2 = 0.75*OM*(1.-G) c GAM3 = 0.25*(2.-3.*G*MU) c GAM4 = 1. - GAM3 c mu1(i) = 0.5 * delta discrete ordinates (Schaller, 1979, Contrib.Atmos.Phys, 52, 17-26): c GAM1 = 0.5*1.7320508*(2. - OM*(1. + G)) c GAM2 = 0.5*1.7320508*OM*(1.-G) c GAM3 = 0.5*(1.-1.7320508*G*MU) c GAM4 = 1. - GAM3 c mu1(i) = 1./sqrt(3.) * Calculations of Associated Legendre Polynomials for GAMA1,2,3,4 * in delta-function, modified quadrature, hemispheric constant, * Hybrid modified Eddington-delta function metods, p633,Table1. * W.E.Meador and W.R.Weaver, GAS,1980,v37,p.630 * W.J.Wiscombe and G.W. Grams, GAS,1976,v33,p2440 c YLM0 = 2. c YLM2 = -3.*G*MU c YLM4 = 0.875*G**3*MU*(5.*MU**2-3.) c YLM6=-0.171875*G**5*MU*(15.-70.*MU**2+63.*MU**4) c YLM8=+0.073242*G**7*MU*(-35.+315.*MU**2-693.*MU**4 c *+429.*MU**6) c YLM10=-0.008118*G**9*MU*(315.-4620.*MU**2+18018.*MU**4 c *-25740.*MU**6+12155.*MU**8) c YLM12=0.003685*G**11*MU*(-693.+15015.*MU**2-90090.*MU**4 c *+218790.*MU**6-230945.*MU**8+88179.*MU**10) c YLMS=YLM0+YLM2+YLM4+YLM6+YLM8+YLM10+YLM12 c YLMS=0.25*YLMS c BETA0 = YLMS c c amu1=1./1.7320508 c YLM0 = 2. c YLM2 = -3.*G*amu1 c YLM4 = 0.875*G**3*amu1*(5.*amu1**2-3.) c YLM6=-0.171875*G**5*amu1*(15.-70.*amu1**2+63.*amu1**4) c YLM8=+0.073242*G**7*amu1*(-35.+315.*amu1**2-693.*amu1**4 c *+429.*amu1**6) c YLM10=-0.008118*G**9*amu1*(315.-4620.*amu1**2+18018.*amu1**4 c *-25740.*amu1**6+12155.*amu1**8) c YLM12=0.003685*G**11*amu1*(-693.+15015.*amu1**2-90090.*amu1**4 c *+218790.*amu1**6-230945.*amu1**8+88179.*amu1**10) c YLMS=YLM0+YLM2+YLM4+YLM6+YLM8+YLM10+YLM12 c YLMS=0.25*YLMS c BETA1 = YLMS c c BETAn = 0.25*(2. - 1.5*G-0.21875*G**3-0.085938*G**5 c *-0.045776*G**7) * Hybrid modified Eddington-delta function(Meador and Weaver,1980,JAS,37,630): c subd=4.*(1.-G*G*(1.-MU)) c GAM1 = (7.-3.*G*G-OM*(4.+3.*G)+OM*G*G*(4.*BETA0+3.*G))/subd c GAM2 =-(1.-G*G-OM*(4.-3.*G)-OM*G*G*(4.*BETA0+3.*G-4.))/subd c GAM3 = BETA0 c GAM4 = 1. - GAM3 c mu1(i) = (1. - g*g*(1.- mu) )/(2. - g*g) ***** * delta function (Meador, and Weaver, 1980, JAS, 37, 630): c GAM1 = (1. - OM*(1. - beta0))/MU c GAM2 = OM*BETA0/MU c GAM3 = BETA0 c GAM4 = 1. - GAM3 c mu1(i) = mu ***** * modified quadrature (Meador, and Weaver, 1980, JAS, 37, 630): c GAM1 = 1.7320508*(1. - OM*(1. - beta1)) c GAM2 = 1.7320508*OM*beta1 c GAM3 = BETA0 c GAM4 = 1. - GAM3 c mu1(i) = 1./sqrt(3.) * hemispheric constant (Toon et al., 1989, JGR, 94, 16287): c GAM1 = 2.*(1. - OM*(1. - betan)) c GAM2 = 2.*OM*BETAn c GAM3 = BETA0 c GAM4 = 1. - GAM3 c mu1(i) = 0.5 ***** * lambda = pg 16,290 equation 21 * big gamma = pg 16,290 equation 22 * if gam2 = 0., then bgam = 0. lam(i) = sqrt(gam1*gam1 - gam2*gam2) IF (gam2 .NE. 0.) THEN bgam(i) = (gam1 - lam(i))/gam2 ELSE bgam(i) = 0. END IF expon = EXP(-lam(i)*taun(i)) * e1 - e4 = pg 16,292 equation 44 e1(i) = 1. + bgam(i)*expon e2(i) = 1. - bgam(i)*expon e3(i) = bgam(i) + expon e4(i) = bgam(i) - expon * the following sets up for the C equations 23, and 24 * found on page 16,290 * prevent division by zero (if LAMBDA=1/MU, shift 1/MU^2 by EPS = 1.E-3 * which is approx equiv to shifting MU by 0.5*EPS* (MU)**3 expon0 = EXP(-tausla(i-1)) expon1 = EXP(-tausla(i)) divisr = lam(i)*lam(i) - 1./(mu2(i)*mu2(i)) temp = AMAX1(eps,abs(divisr)) divisr = SIGN(temp,divisr) up = om*pifs*((gam1 - 1./mu2(i))*gam3 + gam4*gam2)/divisr dn = om*pifs*((gam1 + 1./mu2(i))*gam4 + gam2*gam3)/divisr * cup and cdn are when tau is equal to zero * cuptn and cdntn are when tau is equal to taun cup(i) = up*expon0 cdn(i) = dn*expon0 cuptn(i) = up*expon1 cdntn(i) = dn*expon1 11 CONTINUE ***************** set up matrix ****** * ssfc = pg 16,292 equation 37 where pi Fs is one (unity). ssfc = rsfc*mu*EXP(-tausla(nlayer))*pifs * MROWS = the number of rows in the matrix mrows = 2*nlayer * the following are from pg 16,292 equations 39 - 43. * set up first row of matrix: i = 1 a(1) = 0. b(1) = e1(i) d(1) = -e2(i) e(1) = fdn0 - cdn(i) row=1 * set up odd rows 3 thru (MROWS - 1): i = 0 DO 20, row = 3, mrows - 1, 2 i = i + 1 a(row) = e2(i)*e3(i) - e4(i)*e1(i) b(row) = e1(i)*e1(i + 1) - e3(i)*e3(i + 1) d(row) = e3(i)*e4(i + 1) - e1(i)*e2(i + 1) e(row) = e3(i)*(cup(i + 1) - cuptn(i)) + $ e1(i)*(cdntn(i) - cdn(i + 1)) 20 CONTINUE * set up even rows 2 thru (MROWS - 2): i = 0 DO 30, row = 2, mrows - 2, 2 i = i + 1 a(row) = e2(i + 1)*e1(i) - e3(i)*e4(i + 1) b(row) = e2(i)*e2(i + 1) - e4(i)*e4(i + 1) d(row) = e1(i + 1)*e4(i + 1) - e2(i + 1)*e3(i + 1) e(row) = (cup(i + 1) - cuptn(i))*e2(i + 1) - $ (cdn(i + 1) - cdntn(i))*e4(i + 1) 30 CONTINUE * set up last row of matrix at MROWS: row = mrows i = nlayer a(row) = e1(i) - rsfc*e3(i) b(row) = e2(i) - rsfc*e4(i) d(row) = 0. e(row) = ssfc - cuptn(i) + rsfc*cdntn(i) * solve tri-diagonal matrix: CALL tridiag(a, b, d, e, y, mrows) **** unfold solution of matrix, compute output fluxes: row = 1 lev = 1 j = 1 * the following equations are from pg 16,291 equations 31 & 32 fdr(lev) = EXP( -tausla(0) ) edr(lev) = mu * fdr(lev) edn(lev) = fdn0 eup(lev) = y(row)*e3(j) - y(row + 1)*e4(j) + cup(j) fdn(lev) = edn(lev)/mu1(lev) fup(lev) = eup(lev)/mu1(lev) DO 60, lev = 2, nlayer + 1 fdr(lev) = EXP(-tausla(lev-1)) edr(lev) = mu *fdr(lev) edn(lev) = y(row)*e3(j) + y(row + 1)*e4(j) + cdntn(j) eup(lev) = y(row)*e1(j) + y(row + 1)*e2(j) + cuptn(j) fdn(lev) = edn(lev)/mu1(j) fup(lev) = eup(lev)/mu1(j) row = row + 2 j = j + 1 60 CONTINUE c RETURN END *=============================================================================* SUBROUTINE tridiag(a,b,c,r,u,n) *_______________________________________________________________________ * solves tridiagonal system. From Numerical Recipies, p. 40 *_______________________________________________________________________ IMPLICIT NONE * input: INTEGER n REAL a, b, c, r DIMENSION a(n),b(n),c(n),r(n) * output: REAL u DIMENSION u(n) * local: INTEGER j REAL bet, gam DIMENSION gam(n) *_______________________________________________________________________ IF (b(1) .EQ. 0.) STOP 1001 bet = b(1) u(1) = r(1)/bet DO 11, j = 2, n gam(j) = c(j - 1)/bet bet = b(j) - a(j)*gam(j) IF (bet .EQ. 0.) STOP 2002 u(j) = (r(j) - a(j)*u(j - 1))/bet 11 CONTINUE DO 12, j = n - 1, 1, -1 u(j) = u(j) - gam(j + 1)*u(j + 1) 12 CONTINUE *_______________________________________________________________________ RETURN END