diff --git a/examples/ts1_tsmlt.json b/examples/ts1_tsmlt.json index 65b3b640..7ec8b7c4 100644 --- a/examples/ts1_tsmlt.json +++ b/examples/ts1_tsmlt.json @@ -176,7 +176,8 @@ "name": "jo2_a", "__reaction": "O2 + hv -> O + O1D", "cross section": { - "netcdf files": [ + "apply O2 bands": true, + "netcdf files": [ { "file path": "data/cross_sections/O2_1.nc", "lower extrapolation": { "type": "boundary" }, diff --git a/src/photolysis_rates.F90 b/src/photolysis_rates.F90 index 682ba14d..c9d3493a 100644 --- a/src/photolysis_rates.F90 +++ b/src/photolysis_rates.F90 @@ -294,7 +294,7 @@ subroutine get( this, la_srb, spherical_geometry, grid_warehouse, & !> Local variables character(len=*), parameter :: Iam = "photolysis rates calculator" integer :: vertNdx, rateNdx, nRates - real(dk), allocatable :: airVcol(:), airScol(:) + real(dk), allocatable :: air_vertical_column(:), air_slant_column(:) real(dk), allocatable :: xsqyWrk(:) real(dk), allocatable :: cross_section(:,:) real(dk), allocatable :: quantum_yield(:,:) @@ -343,13 +343,15 @@ subroutine get( this, la_srb, spherical_geometry, grid_warehouse, & ! O2 photolysis can have special la & srb band handling if( any( this%o2_rate_indices_ == rateNdx ) ) then airProfile => profile_warehouse%get_profile( this%air_profile_ ) - allocate( airVcol( airProfile%ncells_ ), & - airScol( airProfile%ncells_ + 1 ) ) - call spherical_geometry%air_mass( airProfile%exo_layer_dens_, airVcol,& - airScol ) - call la_srb%cross_section( grid_warehouse, profile_warehouse, airVcol,& - airScol, cross_section, spherical_geometry ) - deallocate( airVcol, airScol ) + allocate( air_vertical_column( airProfile%ncells_ ), & + air_slant_column( airProfile%ncells_ + 1 ) ) + call spherical_geometry%air_mass( airProfile%exo_layer_dens_, & + air_vertical_column, & + air_slant_column ) + call la_srb%cross_section( grid_warehouse, profile_warehouse, & + air_vertical_column, air_slant_column, & + cross_section, spherical_geometry ) + deallocate( air_vertical_column, air_slant_column ) deallocate( airProfile ) endif diff --git a/test/data/la_srb_bands.config.json b/test/data/la_srb_bands.config.json index 51bd6d98..a00baaa2 100644 --- a/test/data/la_srb_bands.config.json +++ b/test/data/la_srb_bands.config.json @@ -1,18 +1,24 @@ { "grids" : [ { - "name": "height", - "type": "equal interval", - "units": "km", - "begins at" : 0.0, - "ends at" : 120.0, - "cell delta" : 1.0 + "name": "height", + "type": "equal interval", + "units": "km", + "begins at" : 0.5, + "ends at" : 150.5, + "cell delta" : 1.0 }, { - "name": "wavelength", - "type": "from csv file", - "units": "nm", - "file path": "data/grids/wavelength/combined.grid" + "name": "wavelength", + "type": "from csv file", + "units": "nm", + "file path": "test/data/waccm2_ref_101_mod.grid" + }, + { + "name": "LUT wavelength", + "type": "from csv file", + "units": "nm", + "file path": "test/data/waccm2_ref_101.grid" } ], "cross section parameters file": "data/cross_sections/O2_parameters.txt", @@ -26,9 +32,28 @@ "name": "height", "units": "km" } + }, + { + "name": "air", + "type": "air", + "units": "molecule cm-3", + "file path": "data/profiles/atmosphere/ussa.dens" + }, + { + "name": "O2", + "type": "O2", + "units": "molecule cm-3", + "file path": "data/profiles/atmosphere/ussa.dens" } ], - "O2 estimate" :{ - "scale factor": 0.2095 + "O2 cross section": { + "netcdf files": [ + { + "file path": "data/cross_sections/O2_1.nc", + "lower extrapolation": { "type": "boundary" }, + "interpolator": { "type": "fractional target" } + } + ], + "type": "base" } } diff --git a/test/data/waccm2_ref_101.grid b/test/data/waccm2_ref_101.grid new file mode 100644 index 00000000..3988bcef --- /dev/null +++ b/test/data/waccm2_ref_101.grid @@ -0,0 +1,104 @@ + 102 + 120.0000 + 121.0000 + 122.0000 + 123.5000 + 124.3000 + 125.5000 + 126.3000 + 127.1000 + 130.1000 + 131.1000 + 135.0000 + 140.0000 + 145.0000 + 150.0000 + 155.0000 + 160.0000 + 165.0000 + 168.0000 + 171.0000 + 173.0000 + 174.4000 + 177.0000 + 178.6000 + 180.2000 + 181.8000 + 183.5000 + 185.2000 + 186.9000 + 188.7000 + 190.5000 + 192.3000 + 194.2000 + 196.1000 + 198.0000 + 200.0000 + 202.0000 + 204.1000 + 205.8000 + 208.0000 + 211.0000 + 214.0000 + 217.0000 + 220.0000 + 223.0000 + 226.0000 + 229.0000 + 232.0000 + 235.0000 + 238.0000 + 241.0000 + 244.0000 + 247.0000 + 250.0000 + 253.0000 + 256.0000 + 259.0000 + 263.0000 + 267.0000 + 271.0000 + 275.0000 + 279.0000 + 283.0000 + 287.0000 + 291.0000 + 295.0000 + 298.5000 + 302.5000 + 305.5000 + 308.5000 + 311.5000 + 314.5000 + 317.5000 + 322.5000 + 327.5000 + 332.5000 + 337.5000 + 342.5000 + 347.5000 + 350.0000 + 355.0000 + 360.0000 + 365.0000 + 370.0000 + 375.0000 + 380.0000 + 385.0000 + 390.0000 + 395.0000 + 400.0000 + 405.0000 + 410.0000 + 415.0000 + 420.0000 + 430.0000 + 440.0000 + 450.0000 + 500.0000 + 550.0000 + 600.0000 + 650.0000 + 700.0000 + 750.0000 + diff --git a/test/data/waccm2_ref_101_mod.grid b/test/data/waccm2_ref_101_mod.grid new file mode 100644 index 00000000..2fd1304c --- /dev/null +++ b/test/data/waccm2_ref_101_mod.grid @@ -0,0 +1,104 @@ + 102 + 120.0000 + 121.4000 + 121.9000 + 123.5000 + 124.3000 + 125.5000 + 126.3000 + 127.1000 + 130.1000 + 131.1000 + 135.0000 + 140.0000 + 145.0000 + 150.0000 + 155.0000 + 160.0000 + 165.0000 + 168.0000 + 171.0000 + 173.0000 + 175.4000 + 177.0000 + 178.6000 + 180.2000 + 181.8000 + 183.5000 + 185.2000 + 186.9000 + 188.7000 + 190.5000 + 192.3000 + 194.2000 + 196.1000 + 198.0000 + 200.0000 + 202.0000 + 204.1000 + 206.2000 + 208.0000 + 211.0000 + 214.0000 + 217.0000 + 220.0000 + 223.0000 + 226.0000 + 229.0000 + 232.0000 + 235.0000 + 238.0000 + 241.0000 + 244.0000 + 247.0000 + 250.0000 + 253.0000 + 256.0000 + 259.0000 + 263.0000 + 267.0000 + 271.0000 + 275.0000 + 279.0000 + 283.0000 + 287.0000 + 291.0000 + 295.0000 + 298.5000 + 302.5000 + 305.5000 + 308.5000 + 311.5000 + 314.5000 + 317.5000 + 322.5000 + 327.5000 + 332.5000 + 337.5000 + 342.5000 + 347.5000 + 350.0000 + 355.0000 + 360.0000 + 365.0000 + 370.0000 + 375.0000 + 380.0000 + 385.0000 + 390.0000 + 395.0000 + 400.0000 + 405.0000 + 410.0000 + 415.0000 + 420.0000 + 430.0000 + 440.0000 + 450.0000 + 500.0000 + 550.0000 + 600.0000 + 650.0000 + 700.0000 + 750.0000 + diff --git a/test/unit/la_sr_bands.F90 b/test/unit/la_sr_bands.F90 index e5fb6c12..f1cbfed9 100644 --- a/test/unit/la_sr_bands.F90 +++ b/test/unit/la_sr_bands.F90 @@ -67,22 +67,31 @@ program test_la_sr_bands !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine test_optical_depth( ) - use tuvx_grid, only : grid_t - use tuvx_spherical_geometry, only : spherical_geometry_t + + use tuvx_profile, only : profile_t + use tuvx_grid, only : grid_t + use tuvx_spherical_geometry, only : spherical_geometry_t real(dk), allocatable :: air_vertical_column(:), air_slant_column(:) real(dk), allocatable :: o2_optical_depth(:,:) - class(grid_t), pointer :: height_grid => null( ) ! specified altitude working grid [km] - type(spherical_geometry_t) :: spherical_geometry + class(grid_t), pointer :: height_grid ! specified altitude working grid [km] + class(grid_t), pointer :: wavelength_grid ! [nm] + class(spherical_geometry_t), pointer :: spherical_geometry + class(profile_t), pointer :: air height_grid => grid_warehouse%get_grid( "height", "km" ) - allocate( air_vertical_column( height_grid%ncells_ ), & - air_slant_column( height_grid%ncells_ + 1 ) ) + wavelength_grid => grid_warehouse%get_grid( "wavelength", "nm" ) + air => profile_warehouse%get_profile( "air", "molecule cm-3" ) + + spherical_geometry => spherical_geometry_t( grid_warehouse ) + call spherical_geometry%set_parameters( 45.0_dk, grid_warehouse ) - allocate( o2_optical_depth(120, 38) ) + allocate( air_vertical_column( air%ncells_ ), & + air_slant_column( air%ncells_ + 1 ) ) + call spherical_geometry%air_mass( air%exo_layer_dens_, air_vertical_column,& + air_slant_column ) - air_vertical_column(:) = 1 - air_slant_column(:) = 3 + allocate( o2_optical_depth(height_grid%ncells_, wavelength_grid%ncells_) ) o2_optical_depth(:,:) = 0 ! just checking that it runs. This method apparently requires at least @@ -93,6 +102,9 @@ subroutine test_optical_depth( ) spherical_geometry ) deallocate( height_grid ) + deallocate( wavelength_grid ) + deallocate( air ) + deallocate( spherical_geometry ) deallocate( o2_optical_depth ) deallocate( air_vertical_column ) deallocate( air_slant_column ) diff --git a/test/unit/tuv_doug/CMakeLists.txt b/test/unit/tuv_doug/CMakeLists.txt index b5d422d0..d26e4fc5 100644 --- a/test/unit/tuv_doug/CMakeLists.txt +++ b/test/unit/tuv_doug/CMakeLists.txt @@ -20,6 +20,11 @@ target_sources(tuv_doug inter2.f inter3.f inter4.f + la_srb.f + lymana.f + rdo2xs.f + schum.f + ../test_utils.F90 ) add_subdirectory(JCALC) @@ -29,7 +34,7 @@ target_link_libraries(tuv_doug PUBLIC musica::tuvx musica::musicacore) ################################################################################ # Tests refactored configurations based on Doug's data sets and Fortran code -add_executable(test_data_sets data_sets.F90 ../test_utils.F90) +add_executable(test_data_sets data_sets.F90) target_link_libraries(test_data_sets PUBLIC tuv_doug) if(ENABLE_OPENMP) target_link_libraries(test_data_sets PUBLIC OpenMP::OpenMP_Fortran) @@ -37,3 +42,14 @@ endif() add_tuvx_test(data_sets test_data_sets "" ${CMAKE_BINARY_DIR}) ################################################################################ + +# Tests the Lymann-Alpha and Schumann Runge bands calculations in TUV-x against +# Doug's version used to generate the lookup tables used in CAM +add_executable(test_la_srb_lut test_la_srb.F90) +target_link_libraries(test_la_srb_lut PUBLIC tuv_doug) +if(ENABLE_OPENMP) + target_link_libraries(test_la_srb_lut PUBLIC OpenMP::OpenMP_Fortran) +endif() +add_tuvx_test(la_srb_lut test_la_srb_lut "" ${CMAKE_BINARY_DIR}) + +################################################################################ diff --git a/test/unit/tuv_doug/la_srb.f b/test/unit/tuv_doug/la_srb.f new file mode 100644 index 00000000..eddd17e6 --- /dev/null +++ b/test/unit/tuv_doug/la_srb.f @@ -0,0 +1,196 @@ + SUBROUTINE la_srb(nz,z,tlev,nw,wl,o2col,vcol,scol, + $ o2xs1,dto2,o2xs,pathname) +!---------------------------------------------------------------------------! +! PURPOSE: ! +! Compute equivalent optical depths for O2 absorption, and O2 effective ! +! absorption cross sections, parameterized in the Lyman-alpha and SR bands ! +!---------------------------------------------------------------------------! +! PARAMETERS: ! +! NZ - INTEGER, number of specified altitude levels in the working (I)! +! grid ! +! Z - REAL, specified altitude working grid (km) (I)! +! NW - INTEGER, number of specified intervals + 1 in working (I)! +! wavelength grid ! +! WL - REAL, vector of lxower limits of wavelength intervals in (I)! +! working wavelength grid ! +! CZ - REAL, number of air molecules per cm^2 at each specified (I)! +! altitude layer ! +! ZEN - REAL, solar zenith angle (I)! +! ! +! O2XS1 - REAL, O2 cross section from rdo2xs (I)! +! ! +! DTO2 - REAL, optical depth due to O2 absorption at each specified (O)! +! vertical layer at each specified wavelength ! +! O2XS - REAL, molecular absorption cross section in SR bands at (O)! +! each specified altitude and wavelength. Includes Herzberg ! +! continuum. ! +!---------------------------------------------------------------------------! +! EDIT HISTORY: ! +! 02/02 Major revision only over-write LA and SRB ! +! 02/02 add Koppers and delete Korcarts ! +! 02/98 Included Lyman-alpha parameterization ! +! 03/97 Fix dto2 problem at top level (nz) ! +! 02/97 Changed offset for grid-end interpolation to relative number ! +! (x * (1 +- deltax)) ! +! 08/96 Modified for early exit, no redundant read of data and smaller ! +! internal grid if possible; internal grid uses user grid points ! +! whenever possible ! +! 07/96 Modified to work on internal grid and interpolate final values ! +! onto the user-defined grid ! +!---------------------------------------------------------------------------! + implicit none + include 'params' + + integer nz, nw, iz, iw + real wl(kw), z(kz) + real vcol (kz), scol (kz) + real o2col(kz), o2xs1(kw) + real dto2(kz,kw), o2xs(kz,kw) + real secchi(kz) + real tlev(kz) + character*80 pathname + +!---------------------------------------------------------------------- +! Lyman-alpha variables +! O2 optical depth and equivalent cross section in the +! Lyman-alpha region +!---------------------------------------------------------------------- + integer ila, nla, kla + parameter (kla = 2) + real wlla(kla) + real dto2la(kz, kla-1), o2xsla(kz, kla-1) + save ila + +!---------------------------------------------------------------------- +! Grid on which Koppers' parameterization is defined +! O2 optical depth and equivalent cross section on Koppers' grid +!---------------------------------------------------------------------- + integer isrb, nsrb, ksrb + parameter(ksrb = 18) + real wlsrb(ksrb) + real dto2k(kz, ksrb-1), o2xsk(kz, ksrb-1) + save isrb + + integer i + + logical call1 + data call1/.TRUE./ + save call1 + +!---------------------------------------------------------------------- +! Wavelengths for Lyman alpha and SRB parameterizations: +!---------------------------------------------------------------------- + data nla /1/ + data wlla/ 121.0, 122.0/ + + data nsrb /17/ + data wlsrb/174.4, 177.0, 178.6, 180.2, 181.8, 183.5, 185.2, 186.9, + $ 188.7, 190.5, 192.3, 194.2, 196.1, 198.0, 200.0, 202.0, + $ 204.1, 205.8/ + +!---------------------------------------------------------------------- +! initalize O2 cross sections +!---------------------------------------------------------------------- + DO iz = 1, nz + DO iw =1, nw - 1 + o2xs(iz,iw) = o2xs1(iw) + ENDDO + ENDDO + + IF(wl(1) .GT. wlsrb(nsrb)) RETURN + +!---------------------------------------------------------------------- +! On first call, check that the user wavelength grid, WL(IW), is compatible +! with the wavelengths for the parameterizations of the Lyman-alpha and SRB. +! Also compute and save corresponding grid indices (ILA, ISRB) +!---------------------------------------------------------------------- + IF (call1) THEN + +! locate Lyman-alpha wavelengths on grid + + ila = 0 + DO iw = 1, nw + IF(ABS(wl(iw) - wlla(1)) .LT. 10.*precis) THEN + ila = iw + GO TO 5 + ENDIF + ENDDO + 5 CONTINUE + +! check + IF(ila .EQ. 0) STOP ' Lyman alpha grid mis-match - 1' + DO i = 2, nla + 1 + IF(ABS(wl(ila + i - 1) - wlla(i)) .GT. 10.*precis) THEN + WRITE(*,*) 'Lyman alpha grid mis-match - 2' + STOP + ENDIF + ENDDO + +! locate Schumann-Runge wavelengths on grid + isrb = 0 + DO iw = 1, nw + IF(ABS(wl(iw) - wlsrb(1)) .LT. 10.*precis) THEN + isrb = iw + GO TO 6 + ENDIF + ENDDO + 6 CONTINUE + + +! check + IF(isrb .EQ. 0) STOP ' SRB grid mis-match - 1' + DO i = 2, nsrb + 1 + IF(ABS(wl(isrb + i - 1) - wlsrb(i)) .GT. 10.* precis) THEN + WRITE(*,*) ' SRB grid mismatch - w' + STOP + ENDIF + ENDDO + + IF (call1) call1 = .FALSE. + ENDIF + +!---------------------------------------------------------------------- +! Effective secant of solar zenith angle. +! Use 2.0 if no direct sun (value for isotropic radiation) +! For nz, use value at nz-1 +!---------------------------------------------------------------------- + DO i = 1, nz - 1 + secchi(i) = scol(i)/vcol(i) + IF(scol(i) .GT. largest/10.) secchi(i) = 2. + ENDDO + secchi(nz) = secchi(nz-1) + +!--------------------------------------------------------------------- +! Lyman-Alpha parameterization, output values of O2 optical depth +! and O2 effective (equivalent) cross section +!--------------------------------------------------------------------- + CALL lymana(nz,o2col,secchi,dto2la,o2xsla) + + DO iw = ila, ila + nla - 1 + DO iz = 1, nz + dto2(iz,iw) = dto2la(iz, iw - ila + 1) + o2xs(iz,iw) = o2xsla(iz, iw - ila + 1) + ENDDO + ENDDO + +!---------------------------------------------------------------------- +! Koppers' parameterization of the SR bands, output values of O2 +! optical depth and O2 equivalent cross section +!---------------------------------------------------------------------- + + CALL schum(nz,o2col,tlev,secchi,dto2k,o2xsk,pathname) + DO iw = isrb, isrb + nsrb - 1 + DO iz = 1, nz + dto2(iz,iw) = dto2k(iz, iw - isrb + 1) + o2xs(iz,iw) = o2xsk(iz, iw - isrb + 1) + ENDDO + ENDDO + + + RETURN + END + + + + + diff --git a/test/unit/tuv_doug/lymana.f b/test/unit/tuv_doug/lymana.f new file mode 100644 index 00000000..45f3b268 --- /dev/null +++ b/test/unit/tuv_doug/lymana.f @@ -0,0 +1,120 @@ + SUBROUTINE lymana(nz,o2col,secchi,dto2la,o2xsla) + +*-----------------------------------------------------------------------------* +*= PURPOSE: =* +*= Calculate the effective absorption cross section of O2 in the Lyman-Alpha=* +*= bands and an effective O2 optical depth at all altitudes. Parameterized =* +*= after: Chabrillat, S., and G. Kockarts, Simple parameterization of the =* +*= absorption of the solar Lyman-Alpha line, Geophysical Research Letters, =* +*= Vol.24, No.21, pp 2659-2662, 1997. =* +*-----------------------------------------------------------------------------* +*= PARAMETERS: =* +*= NZ - INTEGER, number of specified altitude levels in the working (I)=* +*= grid =* +*= O2COL - REAL, slant overhead O2 column (molec/cc) at each specified (I)=* +*= altitude =* +*= DTO2LA - REAL, optical depth due to O2 absorption at each specified (O)=* +*= vertical layer =* +*= O2XSLA - REAL, molecular absorption cross section in LA bands (O)=* +*-----------------------------------------------------------------------------* +*= EDIT HISTORY: =* +*= 01/98 Original =* +*-----------------------------------------------------------------------------* +*= This program is free software; you can redistribute it and/or modify =* +*= it under the terms of the GNU General Public License as published by the =* +*= Free Software Foundation; either version 2 of the license, or (at your =* +*= option) any later version. =* +*= The TUV package is distributed in the hope that it will be useful, but =* +*= WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTIBI- =* +*= LITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public =* +*= License for more details. =* +*= To obtain a copy of the GNU General Public License, write to: =* +*= Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. =* +*-----------------------------------------------------------------------------* +*= To contact the authors, please mail to: =* +*= Sasha Madronich, NCAR/ACD, P.O.Box 3000, Boulder, CO, 80307-3000, USA or =* +*= send email to: sasha@ucar.edu =* +*-----------------------------------------------------------------------------* +*= Copyright (C) 1994 - 1998 University Corporation for Atmospheric Research =* +*-----------------------------------------------------------------------------* + + IMPLICIT NONE + +* input: + + INCLUDE 'params' + INTEGER nz + REAL o2col(kz) + REAL secchi(kz) + +* output + + REAL dto2la(kz,*), o2xsla(kz,*) + +* local variables + + DOUBLE PRECISION rm(kz), ro2(kz) + DOUBLE PRECISION b(3), c(3), d(3), e(3) + DATA b/ 6.8431D-01, 2.29841D-01, 8.65412D-02/, + > c/8.22114D-21, 1.77556D-20, 8.22112D-21/, + > d/ 6.0073D-21, 4.28569D-21, 1.28059D-20/, + > e/8.21666D-21, 1.63296D-20, 4.85121D-17/ + + INTEGER iz, i + REAL xsmin +*------------------------------------------------------------------------------* +! sm: set minimum cross section + xsmin = 1.D-20 + + DO iz = 1, nz + rm(iz) = 0.D+00 + ro2(iz) = 0.D+00 + DO i = 1, 3 + rm(iz) = rm(iz) + b(i) * DEXP(-c(i) * DBLE(o2col(iz))) + END DO + ! TUV-x logic difference + ! DO i = 1, 2 + DO i = 1, 3 + ro2(iz) = ro2(iz) + d(i) * DEXP(-e(i) * DBLE(o2col(iz))) + ENDDO + + ENDDO + +* calculate effective O2 optical depths and effective O2 cross sections + DO iz = 1, nz-1 + + IF (rm(iz) .GT. 1.0D-100) THEN + IF (ro2(iz) .GT. 1.D-100) THEN + o2xsla(iz,1) = ro2(iz)/rm(iz) + ELSE + ! TUV-x logic difference + ! o2xsla(iz,1) = 0. + o2xsla(iz,1) = xsmin + ENDIF + + IF (rm(iz+1) .GT. 0.) THEN + + dto2la(iz,1) = LOG(rm(iz+1)) / secchi(iz+1) + $ - LOG(rm(iz)) / secchi(iz) + + ELSE + dto2la(iz,1) = 1000. + ENDIF + ELSE + dto2la(iz,1) = 1000. + o2xsla(iz,1) = xsmin + ENDIF + + ENDDO + +* do top layer separately + + dto2la(nz,1) = 0. + IF(rm(nz) .GT. 1.D-100) THEN + o2xsla(nz,1) = ro2(nz)/rm(nz) + ELSE + o2xsla(nz,1) = xsmin + ENDIF + +*------------------------------------------------------------------------------* + END diff --git a/test/unit/tuv_doug/rdo2xs.f b/test/unit/tuv_doug/rdo2xs.f new file mode 100644 index 00000000..2f874245 --- /dev/null +++ b/test/unit/tuv_doug/rdo2xs.f @@ -0,0 +1,117 @@ + SUBROUTINE rdo2xs(nw,wl,wc,o2xs1,pn) +!---------------------------------------------------------------------------! +! PURPOSE: ! +! Read O2 absorption cross section. Except the SR bands and L-alpha line ! +!---------------------------------------------------------------------------! +! 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 ! +!---------------------------------------------------------------------------! +! EDIT HISTORY: ! +! 02/02 By Xuexi ! +!---------------------------------------------------------------------------! + IMPLICIT NONE + INCLUDE 'params' + +!-----------------------------------------------------------------------------! +! ... input ! +!-----------------------------------------------------------------------------! + real, intent(in) :: wl(kw) + real, intent(in) :: wc(kw) + integer, intent(in) :: nw + +!-----------------------------------------------------------------------------! +! ... output ! +!-----------------------------------------------------------------------------! + real, intent(out) :: o2xs1(kw) + +!... Internal + + integer i, iw, n, kdata, ierr + parameter (kdata = 200) + real x1(kdata), y1(kdata) + real x, y + character*80 pn + +!------------------------------------------------------------------------ +! NOTE: Output O2 xsect, is temporary and will be over-written in +! Lyman-alpha and Schumann-Runge wavelength bands. +!------------------------------------------------------------------------ +! ... data +!------------------------------------------------------------------------ +! Read O2 absorption cross section data: +! 116.65 to 203.05 nm = from Brasseur and Solomon 1986 +! 205 to 240 nm = Yoshino et al. 1988 (same as JPL06) +! +! Note that subroutine seto2.f will over-write values in the +! spectral regions corresponding to: +! Lyman-alpha (LA: 121.4-121.9 nm, Chabrillat and Kockarts +! parameterization +! Schumann-Runge bands (SRB: 174.4-205.8 nm, Koppers +! parameteriaztion) +!----------------------------------------------------------------------- + n = 0 + + OPEN(UNIT=kin,FILE=Trim(pn)//'XS_O2_brasseur.txt') + DO i = 1, 7 + READ(kin,*) + ENDDO + DO i = 1, 78 + READ(kin,*) x, y + IF (x .LE. 204.) THEN + n = n + 1 + x1(n) = x + y1(n) = y +! print*, x1(n), y1(n) + ENDIF + ENDDO + CLOSE(kin) + + OPEN(UNIT=kin, + $ FILE=Trim(pn)//'XS_O2_yoshino.txt',STATUS='old') + DO i = 1, 8 + READ(kin,*) + ENDDO + DO i = 1, 36 + n = n + 1 + READ(kin,*) x, y + y1(n) = y*1.E-24 + x1(n) = x +! print*, x1(n), y1(n) + END DO + CLOSE (kin) + +!----------------------------------------------------------------------------- +! Add termination points and interpolate onto the +! user grid (set in subroutine gridw): +!----------------------------------------------------------------------------- + CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),y1(1)) + CALL addpnt(x1,y1,kdata,n,0. ,y1(1)) + 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,o2xs1, n,x1,y1, ierr) + print*, "* interp4 used in rdo2xs.f" + ierr = 0 + CALL inter4(nw,wl,o2xs1, n+1,x1,y1, ierr) +!--------------------------------------------------------------- +! ... Check routine +! do iw = 1,51 +! print*, iw, wc(iw), o2xs1(iw) +! enddo +! stop +!--------------------------------------------------------------- + IF (ierr .NE. 0) THEN + WRITE(*,*) ierr, 'O2 -> O + O' + STOP + ENDIF + + end subroutine rdo2xs + + + + + + diff --git a/test/unit/tuv_doug/schum.f b/test/unit/tuv_doug/schum.f new file mode 100644 index 00000000..3823c8e0 --- /dev/null +++ b/test/unit/tuv_doug/schum.f @@ -0,0 +1,320 @@ + SUBROUTINE schum(nz, o2col,tlev,secchi,dto2,o2xsk,pathname) + +*-----------------------------------------------------------------------------* +*= PURPOSE: =* +*= Calculate the equivalent absorption cross section of O2 in the SR bands. =* +*= The algorithm is based on parameterization of G.A. Koppers, and =* +*= D.P. Murtagh [ref. Ann.Geophys., 14 68-79, 1996] =* +*= Final values do include effects from the Herzberg continuum. =* +*-----------------------------------------------------------------------------* +*= PARAMETERS: =* +*= NZ - INTEGER, number of specified altitude levels in the working (I)=* +*= grid =* +*= O2COL - REAL, slant overhead O2 column (molec/cc) at each specified (I)=* +*= altitude =* +*= TLEV - tmeperature at each level (I)=* +*= SECCHI - ratio of slant to vertical o2 columns (I)=* +*= DTO2 - REAL, optical depth due to O2 absorption at each specified (O)=* +*= vertical layer at each specified wavelength =* +*= O2XSK - REAL, molecular absorption cross section in SR bands at (O)=* +*= each specified wavelength. Includes Herzberg continuum =* +*-----------------------------------------------------------------------------* + + + IMPLICIT NONE + INCLUDE 'params' + + INTEGER nz + REAL o2col(kz), o2col1(kz) + REAL tlev(kz), secchi(kz) + + REAL dto2(kz,17), o2xsk(kz,17) + + CHARACTER*80 pathname + INTEGER i, k, ktop, ktop1, kbot + + REAL XS(17), X + REAL xslod(17) + LOGICAL firstcall + SAVE firstcall + DATA firstcall /.TRUE./ + + DATA xslod /6.2180730E-21, 5.8473627E-22, 5.6996334E-22, + $ 4.5627094E-22, 1.7668250E-22, 1.1178808E-22, + $ 1.2040544E-22, 4.0994668E-23, 1.8450616E-23, + $ 1.5639540E-23, 8.7961075E-24, 7.6475608E-24, + $ 7.6260556E-24, 7.5565696E-24, 7.6334338E-24, + $ 7.4371992E-24, 7.3642966E-24/ +c------------------------------------------ +C Initialize values +c------------------------------------------ + dto2(:,:) = 0.0 + +c------------------------------------------ +c sm Initialize cross sections to values +c sm at large optical depth +c------------------------------------------ + + DO k = 1, nz + DO i = 1, 17 + o2xsk(k,i) = xslod(i) + ENDDO + ENDDO + +c------------------------------------------ +c Loads Chebyshev polynomial Coeff. +c------------------------------------------ + + if (firstcall) then + call INIT_XS(pathname) + firstcall = .FALSE. + endif + +c------------------------------------------ +c Calculate cross sections +c sm: Set smallest O2col = exp(38.) molec cm-2 +c sm to stay in range of parameterization +c sm given by Koppers et al. at top of atm. +c------------------------------------------ + + ktop = nz + kbot = 0 + +c EXP(38.) = 3.185e16 +c EXP(56.) = 2.091e24 + DO k=1,nz !! loop for alt + o2col1(k) = MAX(o2col(k),EXP(38.)) + + x = ALOG(o2col1(k)) + + IF (x .LT. 38.0) THEN + ktop1 = k-1 + write(*,*) ktop1 + ktop = MIN(ktop1,ktop) + ELSE IF (x .GT. 56.0) THEN + kbot = k + ELSE + CALL effxs( x, tlev(k), xs ) + DO i=1,17 + o2xsk(k,i) = xs(i) + END DO + ENDIF + + END DO !! finish loop for alt + +c------------------------------------------ +c fill in cross section where X is out of range +c by repeating edge table values +c------------------------------------------ + +c sm do not allow kbot = nz to avoid division by zero in +c no light case. + + IF(kbot .EQ. nz) kbot = nz - 1 + + DO k=1,kbot + DO i=1,17 + o2xsk(k,i) = o2xsk(kbot+1,i) + END DO + END DO + + DO k=ktop+1,nz + DO i=1,17 + o2xsk(k,i) = o2xsk(ktop,i) + END DO + END DO + +c------------------------------------------ +c Calculate incremental optical depths +c------------------------------------------ + + DO i=1,17 ! loop over wavelength + + DO k=1,nz-1 ! loop for alt + +c... calculate an optical depth weighted by density +c sm: put in mean value estimate, if in shade + + IF (ABS(1. - o2col1(k+1)/o2col1(k)) .LE. 2.*precis) THEN + + dto2(k,i) = o2xsk(k+1,i)*o2col1(k+1)/(nz-1) + + ELSE + + dto2(k,i) = ABS( + $ ( o2xsk(k+1,i)*o2col1(k+1) - o2xsk(k,i)*o2col1(k) ) + $ / ( 1. + ALOG(o2xsk(k+1,i)/o2xsk(k,i)) + $ / ALOG(o2col1(k+1)/o2col1(k)) ) ) + +c... change to vertical optical depth + + dto2(k,i) = 2. * dto2(k,i) / (secchi(k)+secchi(k+1)) + + ENDIF + + END DO + dto2(nz,i) = 0.0 ! set optical depth to zero at top + + + END DO + + return + end + +C------------------------------------------------------------- + SUBROUTINE EFFXS( X, T, XS ) +C------------------------------------------------------------- +C +C Subroutine for evaluating the effective cross section +C of O2 in the Schumann-Runge bands using parameterization +C of G.A. Koppers, and D.P. Murtagh [ref. Ann.Geophys., 14 +C 68-79, 1996] +C +C method: +C ln(xs) = A(X)[T-220]+B(X) +C X = log of slant column of O2 +C A,B calculated from Chebyshev polynomial coeffs +C AC and BC using NR routine chebev. Assume interval +C is 38 Tests against Doug's LA and SR band calculations + + use musica_assert, only : assert, almost_equal + use musica_constants, only : dk => musica_dk + use musica_config, only : config_t + use tuvx_cross_section, only : cross_section_t + use tuvx_grid_warehouse, only : grid_warehouse_t + use tuvx_profile_warehouse, only : profile_warehouse_t + use tuvx_la_sr_bands, only : la_sr_bands_t + use tuvx_test_utils, only : check_values + + implicit none + + character(len=*), parameter :: my_name = "LUT LA/SRB test" + character(len=*), parameter :: conf_l = 'test/data/la_srb_bands.config.json' + type(config_t) :: grid_config, la_config, profile_config, o2_config + class(grid_warehouse_t), pointer :: grids => null( ) + class(profile_warehouse_t), pointer :: profiles => null( ) + class(la_sr_bands_t), pointer :: la_sr_bands => null( ) + class(cross_section_t), pointer :: o2_cross_section => null( ) + character, allocatable :: buffer(:) + + call la_config%from_file( conf_l ) + call la_config%get( "grids", grid_config, my_name ) + call la_config%get( "profiles", profile_config, my_name ) + call la_config%get( "O2 cross section", o2_config, my_name ) + + grids => grid_warehouse_t( grid_config ) + profiles => profile_warehouse_t( profile_config, grids ) + la_sr_bands => la_sr_bands_t( la_config, grids, profiles ) + o2_cross_section => cross_section_t( o2_config, grids, profiles ) + + call compare_o2_cross_sections( la_sr_bands, grids, profiles ) + + deallocate( grids ) + deallocate( profiles ) + deallocate( la_sr_bands ) + deallocate( o2_cross_section ) + +contains + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + subroutine compare_o2_cross_sections( la_sr_bands, grids, profiles ) + + use tuvx_grid, only : grid_t + use tuvx_profile, only : profile_t + use tuvx_spherical_geometry, only : spherical_geometry_t + + class(la_sr_bands_t), intent(inout) :: la_sr_bands + class(grid_warehouse_t), intent(inout) :: grids + class(profile_warehouse_t), intent(inout) :: profiles + + character(len=80) :: file_path + real(dk), allocatable :: air_vertical_column(:), air_slant_column(:) + real(dk), allocatable :: o2_vertical_column(:), o2_slant_column(:) + real(dk), allocatable :: tuvx_o2_optical_depth(:,:), & + tuvx_o2_cross_section(:,:) + real, dimension(151) :: lut_heights, lut_temperature, & + lut_air_vertical_column, lut_air_slant_column, & + lut_o2_column + real, dimension(700) :: lut_wavelength_edges, lut_wavelength_centers, & + lut_o2_base_cross_section + real, dimension(151,700) :: lut_o2_cross_section, lut_o2_optical_depth + class(grid_t), pointer :: heights, wavelengths, lut_wavelengths + class(spherical_geometry_t), pointer :: geometry + class(profile_t), pointer :: air, o2, temperature + real(dk), allocatable :: solar_zenith_angles(:) + integer :: i_sza, i_height, i_wl, n_heights, n_wavelengths, i_output_height + integer :: output_heights(4) = (/ 1, 50, 100, 150 /) + real(dk) :: rel_tol + + heights => grids%get_grid( "height", "km" ) + wavelengths => grids%get_grid( "wavelength", "nm" ) + lut_wavelengths => grids%get_grid( "LUT wavelength", "nm" ) + air => profiles%get_profile( "air", "molecule cm-3" ) + o2 => profiles%get_profile( "O2", "molecule cm-3" ) + temperature => profiles%get_profile( "temperature", "K" ) + geometry => spherical_geometry_t( grids ) + solar_zenith_angles = (/ 0.0_dk, 13.2_dk, 45.0_dk, 87.3_dk, 90.0_dk /) + allocate( air_vertical_column( air%ncells_ ), & + air_slant_column( air%ncells_ + 1 ) ) + allocate( o2_vertical_column( o2%ncells_ ), & + o2_slant_column( o2%ncells_ + 1) ) + allocate( tuvx_o2_optical_depth( heights%ncells_, wavelengths%ncells_ ) ) + lut_heights(:) = huge(1.0) + lut_temperature(:) = huge(1.0) + lut_air_vertical_column(:) = huge(1.0) + lut_air_slant_column(:) = huge(1.0) + lut_o2_column(:) = huge(1.0) + lut_wavelength_edges(:) = huge(1.0) + lut_wavelength_centers(:) = huge(1.0) + lut_o2_base_cross_section(:) = huge(1.0) + lut_o2_cross_section(:,:) = huge(1.0) + lut_o2_optical_depth(:,:) = huge(1.0) + + do i_sza = 1, size( solar_zenith_angles ) + + ! calculate slant O2 column + call geometry%set_parameters( solar_zenith_angles( i_sza ), grids ) + call geometry%air_mass( air%exo_layer_dens_, & + air_vertical_column, & + air_slant_column ) + call geometry%air_mass( o2%exo_layer_dens_, & + o2_vertical_column, & + o2_slant_column ) + + tuvx_o2_optical_depth(:,:) = 0.0_dk + lut_o2_cross_section(:,:) = 0.0 + lut_o2_optical_depth(:,:) = 0.0 + + ! get TUV-x O2 optical depths and cross sections + call la_sr_bands%optical_depth( grids, profiles, air_vertical_column, & + air_slant_column, tuvx_o2_optical_depth, geometry ) + tuvx_o2_cross_section = o2_cross_section%calculate( grids, profiles ) + call la_sr_bands%cross_section( grids, profiles, air_vertical_column, & + air_slant_column, tuvx_o2_cross_section, geometry ) + + ! get LUT O2 optical depths and cross sections + n_heights = heights%ncells_ + 1 + n_wavelengths = lut_wavelengths%ncells_ + 1 + lut_heights(1:n_heights) = real( heights%edge_(:) ) + lut_temperature(1:n_heights) = real( temperature%edge_val_(:) ) + lut_wavelength_edges(1:n_wavelengths) = real( lut_wavelengths%edge_(:) ) + lut_wavelength_centers(1:n_wavelengths-1) = & + real( lut_wavelengths%mid_(:) ) + lut_o2_column(1:n_heights) = real( o2_slant_column(:) ) + lut_air_vertical_column(1:air%ncells_) = real( air_vertical_column(:) ) + lut_air_slant_column(1:air%ncells_+1) = real( air_slant_column(:) ) + + file_path = "test/unit/tuv_doug/INPUT/XSQY/" + call rdo2xs( n_wavelengths, lut_wavelength_edges, & + lut_wavelength_centers, lut_o2_base_cross_section, & + file_path ) + + call la_srb( n_heights, lut_heights, lut_temperature, & + n_wavelengths, lut_wavelength_edges, lut_o2_column, & + lut_air_vertical_column, lut_air_slant_column, & + lut_o2_base_cross_section, lut_o2_optical_depth, & + lut_o2_cross_section, file_path ) + + do i_height = 1, n_heights - 1 + do i_wl = 1, n_wavelengths - 1 + rel_tol = 1.0e-4 + if ( i_wl == 1 .or. i_wl == 3 ) cycle + if ( i_wl == 20 .or. i_wl == 38 ) rel_tol = 0.5_dk + if ( i_wl == 2 .and. i_height >= 112 ) rel_tol = 0.05_dk + call assert( 624510149, & + almost_equal( tuvx_o2_cross_section( i_height, i_wl ), & + real( lut_o2_cross_section( i_height, i_wl ), kind=dk ),& + relative_tolerance = rel_tol ) ) + call assert( 746904813, & + almost_equal( tuvx_o2_optical_depth( i_height, i_wl ), & + real( lut_o2_optical_depth( i_height, i_wl ), kind=dk ),& + relative_tolerance = rel_tol ) ) + end do + end do + deallocate( tuvx_o2_cross_section ) + end do + + deallocate( heights ) + deallocate( wavelengths ) + deallocate( lut_wavelengths ) + deallocate( air ) + deallocate( o2 ) + deallocate( temperature ) + deallocate( geometry ) + deallocate( tuvx_o2_optical_depth ) + deallocate( air_vertical_column ) + deallocate( air_slant_column ) + deallocate( o2_vertical_column ) + deallocate( o2_slant_column ) + + end subroutine compare_o2_cross_sections + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + +end program test_la_srb \ No newline at end of file