Ignore:
Timestamp:
Feb 12, 2020 1:08:46 PM (20 months ago)
Author:
banzhafs
Message:

chemistry model: implemented on-demand emission read mode for LOD 2

File:
1 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/chem_emissions_mod.f90

    r4356 r4403  
    2222! Current revisions:
    2323! ------------------
    24 ! 
     24!
    2525!
    2626! Former revisions:
    2727! -----------------
    2828! $Id$
     29! Implemented on-demand read mode for LOD 2 (ECC)
     30!  - added following module global variables
     31!    - input_file_chem (namesake in netcdf_data_input_mod is local)
     32!    - timestamps
     33!  - added following public subroutines / interfaces
     34!    - chem_emisisons_header_init
     35!    - chem_emisisons_update_on_demand
     36!  - added following local subroutines
     37!    - chem_emisisons_header_init_lod2
     38!    - chem_emisisons_update_on_demand_lod2
     39!  - added following local auxiliary subroutines
     40!    - chem_emissions_init_species ( )
     41!    - chem_emissions_init_timestamps ( )
     42!    - chem_emissions_assign_surface_flux ( )
     43!  - added following local functions
     44!    - chem_emisisons_convert_base_units ( )
     45!    - chem_emissions_mass_2_molar_flux ( )
     46!    - chem_emissions_locate_species ( )
     47!    - chem_emissions_locate_timestep ( )
     48!  - added following error messages
     49!    - CM0468 - LOD mismatch (namelist / chemistry file)
     50!    - CM0469 - Timestamps no in choronological order
     51!  - depreciated unused module variable filename_emis
     52!
     53! 4356 2019-12-20 17:09:33Z suehring
    2954! Minor formatting adjustment
    30 ! 
     55!
    3156! 4242 2019-09-27 12:59:10Z suehring
    3257! Adjust index_hh access to new definition accompanied with new
     
    140165        ONLY:  weight_pres
    141166
     167!
     168!-- 20200203 (ECC)
     169!-- Added new palm_date_time_mod for on-demand emission reading
     170
     171    USE palm_date_time_mod,                                                  &
     172        ONLY:  get_date_time
     173
     174    IMPLICIT NONE
     175
     176!
     177!-- Declare all global variables within the module
     178
     179!
     180!-- 20200203 (ECC) - variable unused
     181!    CHARACTER (LEN=80) ::  filename_emis             !< Variable for the name of the netcdf input file
     182
     183!
     184!-- 20200203 (ECC) new variables for on-demand read mode
     185
     186    CHARACTER(LEN=512), ALLOCATABLE, DIMENSION(:) :: timestamps   !< timestamps in chemistry file
    142187   
    143     IMPLICIT NONE
    144 
    145 !
    146 !-- Declare all global variables within the module
    147    
    148     CHARACTER (LEN=80) ::  filename_emis             !< Variable for the name of the netcdf input file
    149 
    150     INTEGER(iwp) ::  dt_emis                         !< Time Step Emissions
    151     INTEGER(iwp) ::  i                               !< index 1st selected dimension (some dims are not spatial)
    152     INTEGER(iwp) ::  j                               !< index 2nd selected dimension
    153     INTEGER(iwp) ::  i_start                         !< Index to start read variable from netcdf along one dims
    154     INTEGER(iwp) ::  i_end                           !< Index to end read variable from netcdf in one dims
    155     INTEGER(iwp) ::  j_start                         !< Index to start read variable from netcdf in additional dims
    156     INTEGER(iwp) ::  j_end                           !< Index to end read variable from netcdf in additional dims
    157     INTEGER(iwp) ::  len_index                       !< length of index (used for several indices)
    158     INTEGER(iwp) ::  len_index_pm                    !< length of PMs index
    159     INTEGER(iwp) ::  len_index_voc                   !< length of voc index
    160     INTEGER(iwp) ::  z_start                         !< Index to start read variable from netcdf in additional dims
    161     INTEGER(iwp) ::  z_end                           !< Index to end read variable from netcdf in additional dims
     188    CHARACTER(LEN=*),   PARAMETER   :: input_file_chem = 'PIDS_CHEM' !< chemistry file
     189
     190    INTEGER(iwp) ::  dt_emis         !< Time Step Emissions
     191    INTEGER(iwp) ::  i               !< index 1st selected dimension (some dims are not spatial)
     192    INTEGER(iwp) ::  j               !< index 2nd selected dimension
     193    INTEGER(iwp) ::  i_start         !< Index to start read variable from netcdf along one dims
     194    INTEGER(iwp) ::  i_end           !< Index to end read variable from netcdf in one dims
     195    INTEGER(iwp) ::  j_start         !< Index to start read variable from netcdf in additional dims
     196    INTEGER(iwp) ::  j_end           !< Index to end read variable from netcdf in additional dims
     197    INTEGER(iwp) ::  len_index       !< length of index (used for several indices)
     198    INTEGER(iwp) ::  len_index_pm    !< length of PMs index
     199    INTEGER(iwp) ::  len_index_voc   !< length of voc index
     200    INTEGER(iwp) ::  previous_timestamp_index !< index for current timestamp (20200203 ECC)
     201    INTEGER(iwp) ::  z_start         !< Index to start read variable from netcdf in additional dims
     202    INTEGER(iwp) ::  z_end           !< Index to end read variable from netcdf in additional dims
    162203
    163204    REAL(wp) ::  conversion_factor                   !< Units Conversion Factor
     
    185226       MODULE PROCEDURE chem_emissions_setup
    186227    END INTERFACE chem_emissions_setup
    187    
    188     PUBLIC chem_emissions_init, chem_emissions_match, chem_emissions_setup
     228
     229!
     230!-- 20200203 (ECC) new interfaces for on-demand mode
     231
     232!
     233!-- initialization actions for on-demand mode
     234    INTERFACE chem_emissions_header_init
     235       MODULE PROCEDURE chem_emissions_header_init
     236    END INTERFACE chem_emissions_header_init
     237!
     238!-- load emission data for on-demand mode
     239    INTERFACE chem_emissions_update_on_demand
     240       MODULE PROCEDURE chem_emissions_update_on_demand
     241    END INTERFACE chem_emissions_update_on_demand
     242
     243!
     244!-- 20200203 (ECC) update public routines
     245
     246!    PUBLIC chem_emissions_init, chem_emissions_match, chem_emissions_setup
     247
     248    PUBLIC chem_emissions_init, chem_emissions_match, chem_emissions_setup,    &
     249           chem_emissions_header_init, chem_emissions_update_on_demand
    189250!
    190251!-- Public Variables
     
    16891750!-- LSM surfaces
    16901751
     1752
    16911753             DO  m = 1, surf_lsm_h%ns
    16921754
     
    17411803             ENDDO  ! m
    17421804
     1805
     1806
    17431807!
    17441808!-- USM surfaces
     
    17921856
    17931857                ENDIF  ! emis_distribution
    1794  
     1858
    17951859             ENDDO  ! m
    17961860
     
    18081872 END SUBROUTINE chem_emissions_setup
    18091873
     1874
     1875!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     1876!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     1877!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     1878!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     1879!!
     1880!! 20200203 (ECC) - ON DEMAND EMISSION UPDATE MODE
     1881!!
     1882!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     1883!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     1884!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     1885!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     1886
     1887
     1888!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     1889!!
     1890!! WRAPPER / INTERFACE FUNCTIONS
     1891!!
     1892!! NOTE - I find using an explicity wrapper provides much better flow control
     1893!!          over an interface block
     1894!!
     1895!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     1896
     1897!
     1898!-- 20200203 (ECC)
     1899!
     1900!------------------------------------------------------------------------------!
     1901! Description:
     1902! ------------
     1903!>  interface for initiation of emission arrays based on emission LOD
     1904!
     1905!------------------------------------------------------------------------------!
     1906
     1907 SUBROUTINE chem_emissions_header_init
     1908
     1909    IMPLICIT NONE
     1910
     1911    SELECT CASE ( emiss_lod )
     1912       CASE ( 0 )
     1913! do nothing at the moment
     1914       CASE ( 1 )
     1915! do nothing at the moment
     1916       CASE ( 2 )
     1917          CALL chem_emissions_header_init_lod2
     1918    END SELECT
     1919
     1920 END SUBROUTINE chem_emissions_header_init
     1921
     1922
     1923!
     1924!-- 20200203 (ECC)
     1925!
     1926!------------------------------------------------------------------------------!
     1927! Description:
     1928! ------------
     1929!>  interface for initiation of emission arrays based on emission LOD
     1930!
     1931!------------------------------------------------------------------------------!
     1932
     1933 SUBROUTINE chem_emissions_update_on_demand
     1934
     1935    IMPLICIT NONE
     1936
     1937    SELECT CASE ( emiss_lod )
     1938       CASE ( 0 )
     1939! do nothing at the moment
     1940       CASE ( 1 )
     1941! do nothing at the moment
     1942       CASE ( 2 )
     1943          CALL chem_emissions_update_on_demand_lod2
     1944    END SELECT
     1945
     1946 END SUBROUTINE        ! chem_emisisons_update_on_demand
     1947
     1948 
     1949!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     1950!!
     1951!! SUBROUTINES SPECIFIC FOR LOD 2
     1952!!
     1953!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     1954
     1955!
     1956!-- 20200203 (ECC)
     1957!
     1958!------------------------------------------------------------------------------!
     1959! Description:
     1960! ------------
     1961!>  Initiates header for emissions data attributes for LOD 2
     1962!------------------------------------------------------------------------------!
     1963
     1964 SUBROUTINE chem_emissions_header_init_lod2
     1965
     1966    USE control_parameters,                                                 &
     1967        ONLY: coupling_char, message_string
     1968
     1969    USE netcdf_data_input_mod,                                              &
     1970        ONLY: chem_emis_att,                                                &
     1971              open_read_file, close_input_file,                             &
     1972              get_dimension_length, get_variable,                           &
     1973              get_attribute
     1974
     1975    IMPLICIT NONE
     1976
     1977    INTEGER(iwp) ::  ncid      !< chemistry file netCDF handle
     1978    INTEGER(iwp) ::  att_lod   !<  lod attribute in chemistry file
     1979
     1980    IF  ( debug_output )  &
     1981       CALL debug_message( 'chem_emissions_header_init_lod2', 'start' )
     1982
     1983!
     1984!-- opens _chemistry input file and obtain header information
     1985
     1986    CALL open_read_file ( TRIM(input_file_chem) // TRIM(coupling_char), ncid )
     1987!
     1988!-- check if LOD in chemistry file matches LOD in namelist
     1989
     1990    CALL get_attribute ( ncid, 'lod', att_lod, .TRUE. )
     1991   
     1992    IF  ( att_lod /= emiss_lod )  THEN
     1993       message_string = ''   ! to get around unused variable warning / error
     1994       WRITE ( message_string, * )                                           &
     1995             'LOD mismatch between namelist (emiss_lod) and',                &
     1996             CHAR(10), '                    ',                               &
     1997             'chemistry input file (global attributes>lod)'
     1998       CALL message( 'chem_emissions_header_init_lod2', 'CM0468', 1, 2, 0, 6, 0 )
     1999    ENDIF
     2000!
     2001!-- obtain unit conversion factor
     2002
     2003    CALL get_attribute ( ncid, 'units', chem_emis_att%units, .FALSE., "emission_values" )
     2004    conversion_factor = chem_emissions_convert_base_units ( chem_emis_att%units )
     2005!
     2006!-- obtain header attributes
     2007
     2008    CALL chem_emissions_init_species    ( ncid )
     2009    CALL chem_emissions_init_timestamps ( ncid )
     2010!
     2011!-- done reading file
     2012
     2013    CALL close_input_file (ncid)
     2014
     2015!
     2016!-- set previous timestamp index to something different
     2017!-- to trigger a read event later on
     2018
     2019    previous_timestamp_index = -1
     2020 
     2021    IF  ( debug_output )  &
     2022       CALL debug_message( 'chem_emissions_header_init_lod2', 'end' )
     2023
     2024 END SUBROUTINE chem_emissions_header_init_lod2
     2025
     2026!
     2027!-- 20200203 (ECC)
     2028!
     2029!------------------------------------------------------------------------------!
     2030! Description:
     2031! ------------
     2032!> Reads emission data on demand for LOD2
     2033!------------------------------------------------------------------------------!
     2034
     2035 SUBROUTINE chem_emissions_update_on_demand_lod2
     2036
     2037    USE control_parameters,                                                  &
     2038        ONLY: coupling_char,                                                 &
     2039              time_since_reference_point
     2040
     2041    USE netcdf_data_input_mod,                                               &
     2042        ONLY: chem_emis_att,                                                 &
     2043              open_read_file, close_input_file, get_variable
     2044
     2045    USE arrays_3d,                                                           &
     2046        ONLY: pt, hyp
     2047
     2048    USE surface_mod,                                                         &
     2049        ONLY: surf_def_h, surf_lsm_h, surf_usm_h
     2050
     2051
     2052    IMPLICIT NONE
     2053
     2054    CHARACTER(LEN=80) ::  this_timestamp   !< writes out timestamp
     2055
     2056    INTEGER(iwp) ::  i,j,k,m                !< generic counters
     2057    INTEGER(iwp) ::  kmatch                 !< index of matched species
     2058    INTEGER(iwp) ::  ncid                   !< netCDF file handle (chemistry file)
     2059    INTEGER(iwp) ::  time_index_location    !< location of most recent timestamp
     2060
     2061    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:,:) ::  emissions_raw  !< raw emissions data
     2062    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)     ::  emis_distrib   !< surface emissions
     2063    REAL(wp), ALLOCATABLE, DIMENSION(:,:)       ::  mass2mole      !< conversion factor mass 2 molar (ppm) flux
     2064    REAL(wp), ALLOCATABLE, DIMENSION(:,:)       ::  cssws_def_h    !< dummy default surface array
     2065    REAL(wp), ALLOCATABLE, DIMENSION(:,:)       ::  cssws_lsm_h    !< dummy LSM surface array
     2066    REAL(wp), ALLOCATABLE, DIMENSION(:,:)       ::  cssws_usm_h    !< dummy USM surface array
     2067
     2068    IF  ( debug_output )  &
     2069       CALL debug_message ( 'chem_emissions_update_on_demand_lod2', 'start' )
     2070!
     2071!-- obtain current timestamp and locate index
     2072!-- for most recent timestamp element
     2073!-- end subroutine (RETURN) if it is still the same
     2074!-- index as the existing time index
     2075
     2076    this_timestamp = ''   ! string must be initiated before using
     2077    CALL get_date_time( time_since_reference_point, date_time_str=this_timestamp )
     2078
     2079    time_index_location = chem_emissions_locate_timestep                     &
     2080                             ( this_timestamp, timestamps,                   &
     2081                               1, chem_emis_att%dt_emission )
     2082
     2083    IF  ( time_index_location == previous_timestamp_index )  RETURN
     2084
     2085!
     2086!-- begin extract emission data for matched species from netCDF file
     2087
     2088    previous_timestamp_index = time_index_location
     2089
     2090    ALLOCATE ( emis_distrib(n_matched_vars,nys:nyn,nxl:nxr) )
     2091    emis_distrib = 0.0_wp
     2092
     2093!
     2094!-- open netCDF file and allocate temp memory
     2095
     2096    CALL open_read_file( TRIM(input_file_chem) // TRIM(coupling_char), ncid )
     2097    ALLOCATE( emissions_raw(1,1,nys:nyn,nxl:nxr,1) )
     2098
     2099    DO k = 1, n_matched_vars
     2100!
     2101!-- get index for matching species
     2102
     2103       kmatch = chem_emissions_locate_species (                               &
     2104                               spc_names(match_spec_model(k)),                &
     2105                               chem_emis_att%species_name )
     2106!
     2107!-- extract variable as-is
     2108!-- (note C index notations for nx and ny due to MPI and
     2109!-- reversed index dimension order for netCDF Fortran API)
     2110
     2111       emissions_raw = 0.0_wp
     2112
     2113       CALL get_variable ( ncid, 'emission_values', emissions_raw,            &
     2114                           kmatch, nxl+1, nys+1, 1, time_index_location,      &
     2115                           1, nxr-nxl+1, nyn-nys+1, 1, 1, .FALSE. )
     2116!
     2117!-- transfer emission data
     2118
     2119       DO j = nys,nyn
     2120         DO i = nxl,nxr
     2121            emis_distrib(k,j,i) = emissions_raw(1,1,j,i,1) * conversion_factor
     2122         ENDDO
     2123       ENDDO
     2124
     2125    ENDDO  ! k = n_matched_vars
     2126!
     2127!-- netCDF handle and temp memory no longer needed
     2128
     2129    DEALLOCATE( emissions_raw )
     2130    CALL close_input_file( ncid )
     2131!
     2132!-- Set emis_dt since chemistry ODEs can be stiff, the option
     2133!-- to solve them at every RK substep is present to help improve
     2134!-- stability should the need arises
     2135 
     2136    dt_emis = dt_3d
     2137
     2138    IF  ( call_chem_at_all_substeps )                                          &
     2139       dt_emis = dt_emis * weight_pres(intermediate_timestep_count)
     2140!
     2141!-- calculate conversion factor from mass flux to molar flux (mixing ratio)
     2142
     2143    ALLOCATE ( mass2mole(nys:nyn,nxl:nxr) )
     2144    mass2mole = 0.0_wp
     2145
     2146    DO i = nxl, nxr
     2147       DO j = nys, nyn
     2148          mass2mole(j,i) = mass_2_molar_flux ( hyp(nzb), pt(nzb,j,i) )
     2149       ENDDO
     2150    ENDDO
     2151
     2152!
     2153!-- calculate surface fluxes
     2154!-- NOTE - For some reason I cannot pass surf_xxx%cssws as output argument
     2155!--        into subroutine assign_surface_flux ( ).  The contents got mixed up
     2156!--        once the subroutine is finished.  I don't know why and I don't have
     2157!--        time to investigate.  As workaround I declared dummy variables
     2158!--        and reassign them one by one (i.e., in a loop)
     2159!-- ECC 20200206
     2160
     2161!
     2162!-- allocate and initialize dummy surface fluxes
     2163
     2164    ALLOCATE( cssws_def_h(n_matched_vars,surf_def_h(0)%ns) )
     2165    cssws_def_h = 0.0_wp
     2166
     2167    ALLOCATE( cssws_lsm_h(n_matched_vars,surf_lsm_h%ns) )
     2168    cssws_lsm_h = 0.0_wp
     2169
     2170    ALLOCATE( cssws_usm_h(n_matched_vars,surf_usm_h%ns) )
     2171    cssws_usm_h = 0.0_wp
     2172
     2173!
     2174!-- assign and transfer emissions as surface fluxes
     2175
     2176    CALL assign_surface_flux ( cssws_def_h, surf_def_h(0)%ns,                  &
     2177                               surf_def_h(0)%j, surf_def_h(0)%i,               &
     2178                               emis_distrib, mass2mole )
     2179
     2180
     2181    CALL assign_surface_flux ( cssws_lsm_h, surf_lsm_h%ns,                     &
     2182                               surf_lsm_h%j, surf_lsm_h%i,                     &
     2183                               emis_distrib, mass2mole )
     2184
     2185
     2186    CALL assign_surface_flux ( cssws_usm_h, surf_usm_h%ns,                     &
     2187                               surf_usm_h%j, surf_usm_h%i,                     &
     2188                               emis_distrib, mass2mole )
     2189
     2190    DO k = 1, n_matched_vars
     2191
     2192       DO m = 1, surf_def_h(0)%ns
     2193          surf_def_h(0)%cssws(k,m) = cssws_def_h(k,m)
     2194       ENDDO
     2195
     2196       DO m = 1, surf_lsm_h%ns
     2197          surf_lsm_h%cssws(k,m)    = cssws_lsm_h(k,m)
     2198       ENDDO
     2199
     2200       DO m = 1, surf_usm_h%ns
     2201          surf_usm_h%cssws(k,m)    = cssws_usm_h(k,m)
     2202       ENDDO
     2203
     2204    ENDDO
     2205
     2206!
     2207!-- cleaning up
     2208
     2209    DEALLOCATE( cssws_def_h )
     2210    DEALLOCATE( cssws_lsm_h )
     2211    DEALLOCATE( cssws_usm_h )
     2212
     2213    DEALLOCATE ( emis_distrib )
     2214    DEALLOCATE ( mass2mole )
     2215
     2216    IF  ( debug_output )  &
     2217       CALL debug_message ( 'chem_emissions_update_on_demand_lod2', 'end' )
     2218
     2219 END SUBROUTINE        ! chem_emissions_update_on_demand_lod2
     2220
     2221 
     2222!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     2223!!
     2224!! AUXILIARY SUBROUTINES
     2225!!
     2226!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     2227
     2228!
     2229!-- 20200203 (ECC)
     2230!
     2231!------------------------------------------------------------------------------!
     2232! Description:
     2233! ------------
     2234!> look for matched species between emissions attributes and selected
     2235!> chemical mechanisms and determine corresponding molecular weights
     2236!------------------------------------------------------------------------------!
     2237
     2238 SUBROUTINE chem_emissions_init_species ( ncid )
     2239
     2240    USE netcdf_data_input_mod,                                              &
     2241        ONLY: chem_emis_att,                                                &
     2242              open_read_file, close_input_file,                             &
     2243              get_dimension_length, get_variable
     2244
     2245    IMPLICIT NONE
     2246
     2247    INTEGER(iwp) :: ispec  !< generic counter 4 species
     2248
     2249    INTEGER(iwp), INTENT(IN) :: ncid   !< netcdf file ID
     2250
     2251    IF  ( debug_output )  &
     2252       CALL debug_message( 'chem_emissions_init_species', 'start' )
     2253!
     2254!- assign species
     2255
     2256    CALL get_dimension_length ( ncid, chem_emis_att%n_emiss_species, 'nspecies' )
     2257    CALL get_variable ( ncid, 'emission_name', chem_emis_att%species_name,   &
     2258                              chem_emis_att%n_emiss_species )
     2259!
     2260!- backward compatibility for salsa_mod (ECC)
     2261    chem_emis_att%nspec = chem_emis_att%n_emiss_species
     2262!
     2263!-- get a list of matched species between emission_attributes and
     2264!-- selected chemical mechanism
     2265
     2266    emission_output_required = .FALSE.
     2267    CALL  chem_emissions_match( chem_emis_att, n_matched_vars )
     2268
     2269!
     2270!-- if matched species found (at least 1)
     2271!-- allocate memory for emission attributes
     2272!-- assign molecular masses [kg/mol]
     2273!-- see chemistry_model_mod.f90 for reference
     2274
     2275    IF  ( n_matched_vars > 0 )  THEN
     2276
     2277       emission_output_required = .TRUE.
     2278
     2279       ALLOCATE( chem_emis_att%xm(n_matched_vars) )
     2280
     2281       DO  ispec = 1, n_matched_vars
     2282          chem_emis_att%xm(ispec) = 1.0_wp
     2283          SELECT CASE ( TRIM( spc_names(match_spec_model(ispec)) ) )
     2284             CASE ( 'SO2'  );  chem_emis_att%xm(ispec) = xm_S + xm_O * 2
     2285             CASE ( 'SO4'  );  chem_emis_att%xm(ispec) = xm_S + xm_O * 4
     2286             CASE ( 'NO'   );  chem_emis_att%xm(ispec) = xm_N + xm_O
     2287             CASE ( 'NO2'  );  chem_emis_att%xm(ispec) = xm_N + xm_O * 2
     2288             CASE ( 'NH3'  );  chem_emis_att%xm(ispec) = xm_N + xm_H * 3
     2289             CASE ( 'CO'   );  chem_emis_att%xm(ispec) = xm_C + xm_O
     2290             CASE ( 'CO2'  );  chem_emis_att%xm(ispec) = xm_C + xm_O * 2
     2291             CASE ( 'CH4'  );  chem_emis_att%xm(ispec) = xm_C + xm_H * 4
     2292             CASE ( 'HNO3' );  chem_emis_att%xm(ispec) = xm_H + xm_N + xm_O*3
     2293          END SELECT
     2294       ENDDO
     2295   
     2296    ENDIF  ! IF ( n_matched_vars > 0 )
     2297
     2298    IF  ( debug_output )  &
     2299       CALL debug_message( 'chem_emissions_init_species', 'end' )
     2300
     2301 END SUBROUTINE chem_emissions_init_species
     2302
     2303 
     2304!
     2305!-- 20200203 (ECC)
     2306!
     2307!------------------------------------------------------------------------------!
     2308! Description:
     2309! ------------
     2310!> extract timestamps from netCDF input
     2311!------------------------------------------------------------------------------!
     2312
     2313 SUBROUTINE chem_emissions_init_timestamps ( ncid )
     2314
     2315    USE control_parameters,                                                 &
     2316        ONLY: message_string
     2317
     2318    USE netcdf_data_input_mod,                                              &
     2319        ONLY: chem_emis_att,                                                &
     2320              open_read_file, close_input_file,                             &
     2321              get_dimension_length, get_variable
     2322
     2323    IMPLICIT NONE
     2324
     2325    INTEGER(iwp) ::  fld_len   !< string field length
     2326    INTEGER(iwp) ::  itime     !< generic counter (4 species)
     2327
     2328    INTEGER(iwp), INTENT(IN) :: ncid   !< netcdf file handle
     2329
     2330    IF  ( debug_output )  &
     2331       CALL debug_message( 'chem_emissions_read_timestamps', 'start' )
     2332!
     2333!-- import timestamps from netCDF input
     2334
     2335    CALL get_dimension_length ( ncid, chem_emis_att%dt_emission, 'time' )
     2336    CALL get_dimension_length ( ncid, fld_len, 'field_length' )
     2337    CALL get_variable ( ncid, 'timestamp', timestamps, chem_emis_att%dt_emission, fld_len )
     2338!
     2339!-- throw error at first instance of timestamps
     2340!-- not in listed in chronological order
     2341
     2342    DO itime = 2,chem_emis_att%dt_emission
     2343
     2344       IF ( timestamps(itime-1) > timestamps(itime) )  THEN
     2345
     2346           WRITE( message_string, * ) &
     2347                'input timestamps not in chronological order for',             &
     2348                CHAR(10), '                    ',                              &
     2349                'index ', (itime-1), ' : ', TRIM(timestamps(itime-1)), ' and', &
     2350                CHAR(10), '                    ',                              &
     2351                'index ', (itime),   ' : ', TRIM(timestamps(itime))
     2352
     2353          CALL message( 'chem_emissions_read_timestamps', 'CM0469', 1, 2, 0, 6, 0 )
     2354
     2355       ENDIF
     2356
     2357    ENDDO
     2358
     2359    IF  ( debug_output )  &
     2360       CALL debug_message( 'chem_emissions_read_timestamps', 'end' )
     2361
     2362 END SUBROUTINE chem_emissions_init_timestamps
     2363
     2364
     2365!
     2366!-- 20200203 (ECC)
     2367!
     2368!------------------------------------------------------------------------------!
     2369! Description:
     2370! ------------
     2371!> assign emissions as surface fluxes
     2372!
     2373!> NOTE:  For arguments, I originally wanted to use unspecified dimensions,
     2374!>        but I could not get this to work properly, hence the dimensioned
     2375!>        array arguments
     2376!------------------------------------------------------------------------------!
     2377
     2378SUBROUTINE assign_surface_flux ( surf_array, nsurfs, surf_j, surf_i,         &
     2379                                 emis_dist, conv_mole )
     2380
     2381    USE arrays_3d,                                                           &
     2382        ONLY: rho_air
     2383
     2384    USE netcdf_data_input_mod,                                               &
     2385        ONLY: chem_emis_att
     2386
     2387    USE surface_mod   !< for surf_type
     2388
     2389    IMPLICIT NONE
     2390!
     2391!-- input arguments
     2392
     2393    INTEGER(iwp),                    INTENT(IN) ::  nsurfs   !< # surfaces in surf_array
     2394    INTEGER(iwp), DIMENSION(nsurfs), INTENT(IN) ::  surf_i   !< i indices 4 surf. elements
     2395    INTEGER(iwp), DIMENSION(nsurfs), INTENT(IN) ::  surf_j   !< j indices 4 surf. elements
     2396
     2397    REAL(wp), DIMENSION(nys:nyn,nxl:nxr),                INTENT(IN)    ::  conv_mole   !< conv. 2 molar flux
     2398    REAL(wp), DIMENSION(n_matched_vars,nys:nyn,nxl:nxr), INTENT(IN)    ::  emis_dist   !< surf. emissions
     2399
     2400    REAL(wp), DIMENSION(n_matched_vars,nsurfs),          INTENT(INOUT) ::  surf_array  !< surface listing
     2401
     2402!
     2403!-- parameters (magic numbers)
     2404
     2405    CHARACTER(LEN=2),  PARAMETER ::  sp_PM  = 'PM'   !< id string for all PMs
     2406    CHARACTER(LEN=3),  PARAMETER ::  sp_VOC = 'VOC'  !< id string for VOC
     2407
     2408    REAL(wp), PARAMETER ::  mol2ppm = 1.0E+06_wp     !< conversion from mole 2 ppm
     2409!
     2410!-- local variables
     2411
     2412    CHARACTER(LEN=80) ::  this_species_name  !< matched species name
     2413
     2414    INTEGER(iwp) ::  i,j,k,m        !< generic counters
     2415
     2416    REAL(wp) ::  flux_conv_factor   !< conversion factor
     2417
     2418    IF  ( debug_output )  &
     2419       CALL debug_message( 'chem_emissions_header_init_lod2', 'start' )
     2420
     2421    DO k = 1, n_matched_vars
     2422
     2423       this_species_name = spc_names(k)  !< species already matched
     2424
     2425       DO m = 1, nsurfs
     2426
     2427          j = surf_j(m)   ! get surface coordinates
     2428          i = surf_i(m)
     2429
     2430!
     2431!-- calculate conversion factor depending on emission species type
     2432
     2433          flux_conv_factor = rho_air(nzb)
     2434!
     2435!-- account for conversion to different types of emisison species
     2436
     2437          IF       ( TRIM(this_species_name(1:LEN(sp_PM)))  == sp_PM  )  THEN
     2438
     2439                ! do nothing (use mass flux directly)
     2440
     2441          ELSE IF  ( TRIM(this_species_name(1:LEN(sp_VOC))) == sp_VOC )  THEN
     2442
     2443             flux_conv_factor = flux_conv_factor *                           &
     2444                                conv_mole(j,i) * mol2ppm
     2445
     2446          ELSE
     2447
     2448             flux_conv_factor = flux_conv_factor *                           &
     2449                                conv_mole(j,i) * mol2ppm /                   &
     2450                                chem_emis_att%xm(k)
     2451
     2452          ENDIF
     2453!
     2454!-- finally assign surface flux
     2455
     2456          surf_array(k,m) = emis_dist(k,j,i) * flux_conv_factor
     2457
     2458       ENDDO   ! m = 1, nsurfs
     2459
     2460    ENDDO   ! k = 1, n_matched_vars
     2461
     2462
     2463    IF  ( debug_output )  &
     2464       CALL debug_message( 'chem_emissions_header_init_lod2', 'end' )
     2465   
     2466END SUBROUTINE assign_surface_flux
     2467
     2468 
     2469!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     2470!!
     2471!! AUXILIARY FUNCTIONS
     2472!!
     2473!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     2474
     2475!
     2476!-- 20200203 (ECC)
     2477!
     2478!------------------------------------------------------------------------------!
     2479! Description:
     2480! ------------
     2481!> given incoming flux units ( mass / area / time ) provide single-valued
     2482!> conversion factor to ( kg / m2 / s )
     2483!------------------------------------------------------------------------------!
     2484
     2485 FUNCTION chem_emissions_convert_base_units ( units_in ) RESULT ( conv_factor )
     2486
     2487    IMPLICIT NONE
     2488!
     2489!-- function arguments
     2490
     2491    REAL(wp) ::  conv_factor                       !< convertion factor
     2492
     2493    CHARACTER(LEN=*), INTENT(IN) ::  units_in      !< incoming units (ie emt_att%units)
     2494!
     2495!-- parameters (magic numbers)
     2496
     2497    INTEGER(iwp), PARAMETER ::  up2lo = 32         !< convert letter to lower case
     2498!
     2499!-- base unit conversion factors (should be self-explanatory)
     2500
     2501    REAL(wp), PARAMETER ::  g_to_kg       =     1.0E-03_wp
     2502    REAL(wp), PARAMETER ::  miug_to_kg    =     1.0E-09_wp
     2503    REAL(wp), PARAMETER ::  tons_to_kg    =   100.0_wp
     2504
     2505    REAL(wp), PARAMETER ::  hour_per_year =  8760.0_wp
     2506    REAL(wp), PARAMETER ::  s_per_hour    =  3600.0_wp
     2507    REAL(wp), PARAMETER ::  s_per_day     = 86400.0_wp
     2508
     2509    REAL(wp), PARAMETER ::  day_to_s      =     1.0_wp / s_per_day
     2510    REAL(wp), PARAMETER ::  hour_to_s     =     1.0_wp / s_per_hour
     2511    REAL(wp), PARAMETER ::  year_to_s     =     1.0_wp / s_per_hour / hour_per_year
     2512!
     2513!-- local variables
     2514
     2515    CHARACTER(LEN=LEN(units_in)) ::  units_in_lo   !< units in lower case
     2516
     2517    INTEGER(iwp) ::  j,k                           !< generic counters
     2518    INTEGER(iwp) ::  str_len                       !< length of unit string
     2519!
     2520!-- turn units string to lower case
     2521 
     2522    units_in_lo = ''
     2523    str_len = LEN(TRIM(units_in))
     2524
     2525    DO k = 1,str_len
     2526       j = IACHAR( units_in(k:k) )
     2527       units_in_lo(k:k) = ACHAR(j)
     2528       IF  ( (j>=IACHAR("A")) .AND. (j<=IACHAR("Z")) )   &
     2529          units_in_lo(k:k) = ACHAR ( j + up2lo )
     2530    ENDDO
     2531
     2532    conv_factor = 1.0_wp     !< default value (kg/m2/s)
     2533
     2534    SELECT CASE ( TRIM( units_in_lo ) )
     2535       CASE ( 'kg/m2/s'            );   conv_factor = 1.0_wp
     2536       CASE ( 'kg/m2/hour'         );   conv_factor = hour_to_s
     2537       CASE ( 'kg/m2/day'          );   conv_factor = day_to_s
     2538       CASE ( 'kg/m2/year'         );   conv_factor = year_to_s
     2539       CASE ( 'ton/m2/s'           );   conv_factor = tons_to_kg
     2540       CASE ( 'ton/m2/hour'        );   conv_factor = tons_to_kg * hour_to_s
     2541       CASE ( 'ton/m2/day'         );   conv_factor = tons_to_kg * day_to_s
     2542       CASE ( 'ton/m2/year'        );   conv_factor = tons_to_kg * year_to_s
     2543       CASE ( 'g/m2/s'             );   conv_factor = g_to_kg
     2544       CASE ( 'g/m2/hour'          );   conv_factor = g_to_kg * hour_to_s
     2545       CASE ( 'g/m2/day'           );   conv_factor = g_to_kg * day_to_s
     2546       CASE ( 'g/m2/year'          );   conv_factor = g_to_kg * year_to_s
     2547       CASE ( 'micrograms/m2/s'    );   conv_factor = miug_to_kg
     2548       CASE ( 'micrograms/m2/hour' );   conv_factor = miug_to_kg * hour_to_s
     2549       CASE ( 'micrograms/m2/day'  );   conv_factor = miug_to_kg * day_to_s
     2550       CASE ( 'micrograms/m2/year' );   conv_factor = miug_to_kg * year_to_s
     2551       CASE DEFAULT
     2552          message_string = ''   ! to get around unused variable warning / error
     2553          WRITE ( message_string, * ) 'Specified emission units (',          &
     2554                                      TRIM(units_in),                        &
     2555                                      ') not recognized in PALM-4U'
     2556          CALL message ( 'chem_emission_convert_units', 'CM0446', 2, 2, 0, 6, 0 )
     2557    END SELECT
     2558
     2559 END FUNCTION chem_emissions_convert_base_units
     2560
     2561
     2562!
     2563!-- 20200203 (ECC)
     2564!
     2565!------------------------------------------------------------------------------!
     2566! Description:
     2567! ------------
     2568!> calculates conversion factor from mass flux to ppm (molar flux)
     2569!------------------------------------------------------------------------------!
     2570
     2571 FUNCTION mass_2_molar_flux ( rhogh, theta ) RESULT ( conv_factor )
     2572
     2573    USE basic_constants_and_equations_mod,                                   &
     2574        ONLY: p_0, rgas_univ, rd_d_cp
     2575
     2576    IMPLICIT NONE
     2577!
     2578!-- function arguments
     2579
     2580    REAL(wp)             ::  conv_factor   !< conversion factor
     2581    REAL(wp), INTENT(IN) ::  rhogh         !< hydrostatic pressure
     2582    REAL(wp), INTENT(IN) ::  theta         !< potential temperature
     2583
     2584    conv_factor = ( rgas_univ / rhogh ) * theta * ( ( rhogh / p_0 ) ** rd_d_cp )
     2585
     2586 END FUNCTION mass_2_molar_flux
     2587
     2588
     2589!
     2590!-- 20200203 (ECC)
     2591!
     2592!------------------------------------------------------------------------------!
     2593! Description:
     2594! ------------
     2595!> given target sepecies locate index in species array
     2596!> returns 0 if none is found
     2597!------------------------------------------------------------------------------!
     2598
     2599 FUNCTION chem_emissions_locate_species ( this_species, species_array )      &
     2600          RESULT ( species_index )
     2601
     2602    IMPLICIT NONE
     2603!
     2604!-- function arguments
     2605
     2606    INTEGER(iwp) ::  species_index   !> index matching species
     2607
     2608    CHARACTER(LEN=*),  INTENT(IN) ::  this_species       !> target species
     2609    CHARACTER(LEN=25), INTENT(IN) ::  species_array(:)   !> array of species
     2610!
     2611!-- local variables
     2612
     2613    INTEGER(iwp) :: k           !> generic counter
     2614    INTEGER(iwp) :: n_species   !> number of species in species_array
     2615
     2616    n_species = SIZE( species_array, 1 )
     2617
     2618    DO k = 1, n_species
     2619       IF ( TRIM(species_array(k)) == TRIM(this_species) )  EXIT
     2620    ENDDO
     2621
     2622    species_index = 0   !> assume no matching index is found
     2623
     2624    IF ( TRIM(species_array(k)) == TRIM(this_species) )  specieS_index = k
     2625
     2626 END FUNCTION chem_emissions_locate_species
     2627
     2628
     2629!
     2630!-- 20200203 (ECC)
     2631!
     2632!------------------------------------------------------------------------------!
     2633! Description:
     2634! ------------
     2635!> given target timestamp locate most recent timestep in timestamp array
     2636!> using bisection search (since array is sorted)
     2637!------------------------------------------------------------------------------!
     2638
     2639 RECURSIVE FUNCTION chem_emissions_locate_timestep                           &
     2640                       ( this_timestamp, timestamp_array,                    &
     2641                         lower_bound, upper_bound )                          &
     2642                       RESULT ( timestamp_index )
     2643
     2644!
     2645!-- function arguments
     2646
     2647    INTEGER(iwp) ::  timestamp_index   !> index for most recent timestamp in timestamp_array
     2648
     2649    CHARACTER(LEN=*),   INTENT(IN) ::  this_timestamp       !> target timestamp
     2650    CHARACTER(LEN=512), INTENT(IN) ::  timestamp_array(:)   !> array of timestamps
     2651
     2652    INTEGER(iwp), INTENT(IN) ::  lower_bound, upper_bound   !> timestamp_array index bounds
     2653
     2654!
     2655!-- local variables
     2656
     2657    INTEGER(iwp) :: k0,km,k1          !> lower, central, and upper index bounds
     2658!
     2659!-- assign bounds
     2660
     2661    k0 = lower_bound
     2662    k1 = upper_bound
     2663!
     2664!-- make sure k1 is always not smaller than k0
     2665
     2666    IF  ( k0 > k1 )  THEN
     2667       k0 = upper_bound
     2668       k1 = lower_bound
     2669    ENDIF
     2670!
     2671!-- make sure k0 and k1 stay within array bounds by timestamp_array
     2672
     2673    IF  ( k0 < 1                       )  k0 = 1
     2674    IF  ( k1 > SIZE(timestamp_array,1) )  k1 = SIZE(timestamp_array,1)
     2675!
     2676!-- terminate if target is contained within 2 consecutive indices
     2677!-- otherwise calculate central bound (km) and determine new
     2678!-- index bounds for the next iteration
     2679
     2680    IF ( ( k1 - k0 ) > 1 )  THEN
     2681       km = ( k0 + k1 ) / 2
     2682       IF  ( TRIM(this_timestamp) > TRIM(timestamp_array(km)) )  THEN
     2683          k0 = km
     2684       ELSE
     2685          k1 = km
     2686       ENDIF
     2687       timestamp_index = chem_emissions_locate_timestep &
     2688                         ( this_timestamp, timestamp_array, k0, k1 )
     2689    ELSE
     2690       timestamp_index = k0
     2691    ENDIF
     2692
     2693 END FUNCTION chem_emissions_locate_timestep
     2694
     2695
     2696!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     2697!!
     2698!! END OF MODULE
     2699!!
     2700!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     2701
    18102702 END MODULE chem_emissions_mod
Note: See TracChangeset for help on using the changeset viewer.