search

metno / emep-ctm

Open Source EMEP/MSC-W model
GNU General Public License v3.0
29 stars 19 forks source link

Reading in hourly emissions from NetCDF #94

Closed rubenww closed 2 years ago

rubenww commented 2 years ago

Hi, I'm trying to read in hourly emissions from a NetCDF file for a specific region (NL) and a specific component (NH3). These emissions should be used in place of the yearly emissions that I'd want to use for the rest of the components. Following the manual on the new emission format, I have added the following to my namelist file:

! yearly emissions
emis_inputlist(1)%name     = '/data/gridPOLL', ! high resolution NL emissions
emis_inputlist(1)%PollName(1:) = 'sox', 'nox', 'co', 'voc', 'pm25', 'pmco'         ! but exclude NH3
emis_inputlist(2)%name     = '/data/Emis_TNO7grid_EMEP2014.nc',       ! emissions outside NL

emis_inputlist(1)%incl(1:) = 'NL',
emis_inputlist(2)%excl(1:) = 'NL','NOS',

! hourly emissions
Emis_sourceFiles(1)%filename       = '2015_NH3.nc',                                ! high resolution NH3 emissions, hourly
Emis_sourceFiles(1)%projection     = 'lon lat',
Emis_sourceFiles(1)%periodicity    = 'hourly',                                     ! contains hourly fields
Emis_sourceFiles(1)%source(1)%varname     = 'SNAP1',                    ! netcdf variable name
Emis_sourceFiles(1)%source(1)%species     = "NH3         ",                ! 'nh3' was not found
Emis_sourceFiles(1)%source(1)%sector      = 1,
Emis_sourceFiles(1)%source(1)%units       = 'mg/m2/h',
Emis_sourceFiles(1)%source(2)%varname     = 'SNAP2', 
Emis_sourceFiles(1)%source(2)%species     = "NH3         ",
Emis_sourceFiles(1)%source(2)%sector      = 2,

And so on up to SNAP sector 11. The netcdf file '2015_NH3.nc' looks as follows:

netcdf \2015_NH3 {
dimensions:
        lon = 222 ;
        lat = 201 ;
        time = 8760 ;
variables:
        double lon(lon) ;
                lon:units = "degrees_east" ;
                lon:long_name = "longitude" ;
                lon:axis = "X" ;
        double lat(lat) ;
                lat:units = "degrees_north" ;
                lat:long_name = "latitude" ;
                lat:axis = "Y" ;
        double time(time) ;
                time:units = "days since 1900-01-01 00:00:00.0" ;
                time:long_name = "time" ;
                time:calendar = "standard" ;
                time:axis = "T" ;
        float SNAP1(time, lat, lon) ;
                SNAP1:units = "mg" ;
                SNAP1:_FillValue = -999.f ;
                SNAP1:long_name = "Emis_mg_NH3_SNAP1" ;
        float SNAP2(time, lat, lon) ;
                SNAP2:units = "mg" ;
                SNAP2:_FillValue = -999.f ;
                SNAP2:long_name = "Emis_mg_NH3_SNAP2" ;

etc. for the rest of the SNAP sectors. In ncview: image where the pure white indicate the NaN/missing values.

In the log files, I see the following:

 Reading emissions for year        2015
 Emission source number            1 from 
 Emissions:/data/gridPOLL
 Emissions:/data/gridPOLL INPUTLIST-INCL
           1
 including only countries: NL 
including pollutant sox from /data/gridPOLL
including pollutant nox from /data/gridPOLL
including pollutant co from /data/gridPOLL
including pollutant voc from /data/gridPOLL
including pollutant pm25 from /data/gridPOLL
including pollutant pmco from /data/gridPOLL
original emission name /data/gridPOLL
filename redefined as: /data/gridPOLL
 Emission source number            2 from 
 Emissions:/data/Emis_TNO7grid_EMEP2014.nc
 excluding countries: NL NOS 
original emission name /data/Emis_TNO7grid_EMEP2014.nc
filename redefined as: /data/Emis_TNO7grid_EMEP2014.nc
femis: reductions:
  1    1        1.79      1.78      3.00      1.29      1.18      1.15      5.43

and then the list of femis reductions, so it seems the yearly emissions are read in:

#EMTAB Total emissions by countries for /data/gridPOLL (Gg)
EMTAB CC Land                 sox         nox         co          voc         nh3         pm25        pmco
EMTAB  17 NL                 42.72      314.93      578.40      141.23        0.00       15.82       13.98
EMTAB 998        EU          42.72      314.93      578.40      141.23        0.00       15.82       13.98
Emissions: reading emis_inputlist    2/data/Emis_TNO7grid_EMEP2014.nc
#EMTAB Total emissions by countries for /data/Emis_TNO7grid_EMEP2014.nc (Gg)
EMTAB CC Land                 sox         nox         co          voc         nh3         pm25        pmco
EMTAB   3 BE                 33.74      146.39      285.47       88.69       41.23       19.23        7.19
EMTAB   8 FR                  0.14        0.34        1.00        0.10        0.22        0.11        0.06
EMTAB  60 DE                 98.44      237.61      613.23      136.06       68.77       14.04       19.56
EMTAB 998        EU         132.32      384.35      899.69      224.84      110.22       33.39       26.81
EMTAB 999 TOTA         175.04      699.28     1478.09      366.07      110.22       49.21       40.79

For the hourly emissions, I get output of the following type:

 Ini_Em:Initializing Emissions from 2015_NH3.nc
 Will set default resolution (it does not need to be exact) 
   2056.78325701306     
writing config attribute on SNAP1   1   1   1  11
 Ini_Em:country found N/A          68          11
Ini_Em: species found N/A 104 NH3
 Ini_Em:Final emission source parameters 
 SNAP1                                                           
 NH3                                                             
 mg/m2/h                                                         
 N/A                                                             
 hourly                                                                     1
         104          21   1.00000000000000      T T 
 NOTSET                                                          
 NOTSET                                                                    -1
          -1          68           0 F           1           1           1
           1 F           1
including SNAP1 as NH3 sector  1 country N/A

And then

NH3 included in nh3   1   1  68  68   5NH3       1.000
NH3 included in nh3   2   2  68  68   5NH3       1.000
NH3 included in nh3   3   3  68  68   5NH3       1.000
NH3 included in nh3   4   4  68  68   5NH3       1.000
NH3 included in nh3   5   5  68  68   5NH3       1.000
NH3 included in nh3   6   6  68  68   5NH3       1.000
NH3 included in nh3   7   7  68  68   5NH3       1.000
NH3 included in nh3   8   8  68  68   5NH3       1.000
NH3 included in nh3   9   9  68  68   5NH3       1.000
NH3 included in nh3  10  10  68  68   5NH3       1.000
NH3 included in nh3  11  11  68  68   5NH3       1.000

Later, it seems that the emission data is updated from the file:

 Emis: current date has reached past update date            0           0
           0           0           0 
 hourly                                                          
           1  getemis mg/m2/h SNAP1 read as NH3
           1 ENDOFVAL         2015           1           1           1
        1800
 SNAP1 reading new emis from 2015_NH3.nc, record            1
 ReadCDF start: 2015_NH3.nc:SNAP1
 ReadCDF UnDef defined as:   0.000000000000000E+000
 ReadCDF reading 2015_NH3.nc nstart            1
 ReadCDF variable exists: SNAP1
 minlat attribute not found for SNAP1
  FillValue (not counted)  -999.000000000000     
 ReadCDF size variable            1         222
 ReadCDF size variable            2         201
 ReadCDF size variable            3        8760
 data known_projection lon lat
alloc_flag lon lat lon lat    0  222x  201
 ReadCDF using lon namelon
 ReadCDF using lat namelat
 ReadCDF lon bounds   3.26000000000000        7.35000000000000     
 ReadCDF lat bounds   50.5100000000000        54.2100000000000     
 interpol_used: conservative                                                    
 SET Grid resolution:2015_NH3.nc   2056.78325701306     
ReadCDF LL values     3.28    3.26   54.03   50.53   50.51   54.05      3.94    3.25   51.43   50.50
 Grid does not cover 360 degrees    3.26000000000000     
   7.35000000000000     
ReadCDF minmax values     3.25    3.94   50.50   51.43       1      38       1      51
ReadCDF reading chunk        1       1      38
ReadCDF reading chunk        2       1      51
ReadCDF reading chunk        3       1       1
  xtype real            5
 interpol case conservative                                                    
 interpol case conservative projection lon lat

And the same for the rest of the SNAP sectors. However, then finally:

 Emissions per country for 2015_NH3.nc (Gg/year) 
EMTAB CC Land                 sox         nox         co          voc         nh3         pm25        pmco
EMTAB 999 TOTAL               0.00        0.00        0.00        0.00         NaN        0.00        0.00

So it seems there is a problem with NaN's that causes the program to SEGFAULT later on. First, I get a lot of these kind of warnings: RPMARES: not safe to divide...exit? 2 1 0.45 1.22E-03 NaN 1.18E-03 1.47E-01 NaN NaN 7.136E-01 NaN 1.176E-01 6.680E-02 7.571E-03 NaN 4.227E-02 4.530E-01 2.736E+02 And then after current date and time: 2015-01-01 00:05:00, the program crashes:

forrtl: severe (174): SIGSEGV, segmentation fault occurred
Image              PC                Routine            Line        Source             
libifcoremt.so.5   00002B176D7539A9  for__signal_handl     Unknown  Unknown
libpthread-2.17.s  00002B177227C630  Unknown               Unknown  Unknown
emepctm            00000000004A2378  codep_mod_mp_code         171  CoDep_mod.f90
emepctm            0000000000836A2B  rsurface_mod_mp_r         226  Rsurface_mod.f90
emepctm            0000000000574016  drydep_mod_mp_dry         460  DryDep_mod.f90
emepctm            000000000083A5D4  runchem_mod_mp_ru         275  Runchem_mod.f90
emepctm            00000000008D04F9  phychem_mod_mp_ph         237  PhyChem_mod.f90
emepctm            0000000000585E0E  MAIN__                    362  emep_Main.f90
emepctm            0000000000403D82  Unknown               Unknown  Unknown
libc-2.17.so       00002B17726AF555  __libc_start_main     Unknown  Unknown
emepctm            0000000000403C89  Unknown               Unknown  Unknown

Presumably because of the NaNs.

Is this kind of functionality possible? Any help on how to get this to work would be greatly appreciated! All the best, Ruben

gitpeterwind commented 2 years ago

Hi Ruben, I don't know what happened here. It should work. Could you send me the hourly emission file? (I can send you an upload link if you prefer) Also if you can sen me config_emep.nml and the entire log, so that I can find other info, such as projection etc.

rubenww commented 2 years ago

Thank you Peter, the text files are below: config_emep.nml.txt run001_2015_2015_04.out.txt

I have sent a link to download the emissions file to your met.no email address. Thank you so much for taking a look. Best, Ruben

gitpeterwind commented 2 years ago

Thanks for the log and config, but I did not get the emissions. Perhaps the file was too large to be sent by mail? I can send you an upload link if you send me your email address.

rubenww commented 2 years ago

My apologies, the file was still uploading. You should receive an email soon.

gitpeterwind commented 2 years ago

Thanks, I got it now! I can already see a potential problem: In the variable attribute you define SNAP7:_FillValue = -999.f ; , but the values of the variable in the undefined region are different (9.9e36 instead of -999). I am anyway not 100% sure changing FillValue would work anyway, since it is often challenging to test the equality of two real numbers. Could you use zero in the undefined regions instead?

gitpeterwind commented 2 years ago

A more robust alternative where you can keep the emission file unchanged, (and that I will also consider implement in the next version), would be to add a condition <1.E30 in the NetCDF_mod.f90 (lines 3369 and 3446):

 elseif( Rvalues(igjgk)<1E30 .and. (OnlyDefinedValues.or.Rvalues(igjgk)/=FillValue))then
rubenww commented 2 years ago

Thank you for the advice Peter, I have changed the NaNs to zero and the model seems to be running now, I'll update whether it worked OK once it has completed the test run.

I agree that it would be good to add a condition like that to be able to specify missing values. Just to be clear - if I set the missing values to zero now, would they overwrite the values from emis_inputlist(2)%name = '/data/Emis_TNO7grid_EMEP2014.nc' at the same locations or is a sum taken in that case?

gitpeterwind commented 2 years ago

Good! Emissions are always additive: they cannot be overwritten once taken into account.

rubenww commented 2 years ago

After some more thorough testing, I can confirm everything is working as expected. Thanks again for the help!