|
|
(14 intermediate revisions by 2 users not shown) |
Line 35: |
Line 35: |
| === Modification of MAP_A2A for use within GEOS-Chem === | | === Modification of MAP_A2A for use within GEOS-Chem === |
|
| |
|
| '''[mailto:COOPERM2@dal.ca Matt Cooper]''' has replaced the previous GEOS-Chem regridding routines with an algorithm adapted from the [[Regridding in GEOS-Chem#The_MAP_A2A_algorithm|<tt>MAP_A2A</tt> regridding package]] (developed by S.-J. Lin and refined by Bob Yantosca). The previous GEOS-Chem regridding routines sometimes involved two separate regridding processes, passing through the [http://acmg.seas.harvard.edu/geos/doc/man/appendix_2.html#A2.5 GEOS 1° x 1° grid], which could lead to loss of information. The <tt>MAP_A2A</tt> algorithm regrids emissions from any arbitrary horizontal grid to the current model resolution. | | '''[mailto:COOPERM2@dal.ca Matt Cooper]''' has replaced the previous GEOS-Chem regridding routines with an algorithm adapted from the [[Regridding in GEOS-Chem#The_MAP_A2A_algorithm|<tt>MAP_A2A</tt> regridding package]] (developed by S.-J. Lin and refined by Bob Yantosca). The previous GEOS-Chem regridding routines sometimes involved two separate regridding processes, passing through the [[GEOS-Chem_horizontal_grids#Generic_1_x_1_grid|GEOS 1° x 1° grid]], which could lead to loss of information. The <tt>MAP_A2A</tt> algorithm regrids emissions from any arbitrary horizontal grid to the current model resolution. |
|
| |
|
| <tt>MAP_A2A</tt> was introduced into GEOS-Chem v9-01-03 (and tested with 1-month benchmark [[GEOS-Chem v9-01-03 benchmark history#v9-01-03k|'''v9-01-03k''']], which was approved on 27 Apr 2012. | | <tt>MAP_A2A</tt> was introduced into GEOS-Chem v9-01-03 (and tested with 1-month benchmark [[GEOS-Chem v9-01-03 benchmark history#v9-01-03k|'''v9-01-03k''']], which was approved on 27 Apr 2012. |
Line 54: |
Line 54: |
|
| |
|
| --[[User:Bmy|Bob Y.]] 15:50, 28 August 2012 (EDT) | | --[[User:Bmy|Bob Y.]] 15:50, 28 August 2012 (EDT) |
|
| |
| == Previous issues that have are now resolved ==
| |
|
| |
| === Bug fixes for the MAP_A2A regridding algorithm ===
| |
|
| |
| <span style="color:green">'''''This update was tested in the 1-month benchmark simulation [[GEOS-Chem v9-01-03 benchmark history#v9-01-03n|v9-01-03n]] and approved on 08 Jun 2012.'''''</span>
| |
|
| |
| Lee Murray found a typo in the latlon_geos1x1.txt file used for the MAP_A2A regridding algorithm.
| |
|
| |
| '''''[mailto:ltmurray@seas.harvard.edu Lee Murray] wrote:'''''
| |
| :There is definitely an error in the erroneous missing of 0.5°W in the <tt>latlon_geos1x1.txt</tt> file. This led to inherent inconsistencies in how <tt>latlon_geos1x1.txt</tt> and <tt>latlon_generic.txt</tt> were indexed by GEOS-Chem, since <tt>latlon_geos1x1</tt> made up for the missing longitude edge by doubling up both 180.5°W/179.5°E and 179.5°W/180.5°E. <tt>latlon_generic.txt</tt> only repeats 180°W/180°E.
| |
|
| |
| Several bugs were also fixed in the MAP_A2A regridding code to ensure consistency in the use of <tt>IM</tt> and <tt>JM</tt>, the number of longitude and latitude centers on the input grid.
| |
|
| |
| '''''[mailto:psk9@duke.edu Prasad Kasibhatla] wrote:'''''
| |
| :I think we are all on the same page, but just to make sure here is what I suggest we agree on:
| |
|
| |
| :(1) That the calling sequence for <tt>DO_REGRID_A2A</tt> be:
| |
| CALL DO_REGRID_A2A(FILENAME,IM,JM,INGRID,OUTGRID,PERAREA)
| |
| :where
| |
| ::<tt>IM</tt> = number of grid boxes in E-W direction (not number of lon edges)
| |
| ::<tt>JM</tt> = number of grid boxes in N-S direction (not number of lat edges)
| |
| :so that the dimension of <tt>INGRID</tt> is <tt>(IM,JM)</tt> and <tt>IM+1</tt> and <tt>JM+1</tt> are number of lon edges and number of lat edges, respectively.
| |
|
| |
| :(2) That the calling sequence for <tt>MAP_A2A</tt> be
| |
| CALL MAP_A2A( IM, JM, INLON, INSIN, INGRID, &
| |
| IIPAR, JJPAR, LON2, SIN2, OUTGRID, 0, 0)
| |
| :where <tt>INLON</tt>, <tt>INSIN</tt>, and <tt>INGRID</tt> are dimensioned <tt>IM+1</tt>, <tt>JM+1</tt>, and <tt>(IM,JM)</tt>, respectively, and <tt>LON2</tt>, <tt>SIN2</tt>, and <tt>OUTGRID</tt> are dimensioned <tt>IIPAR+1</tt>, <tt>JJPAR+1</tt>, and <tt>(IIPAR, JJPAR)</tt>, respectively.
| |
|
| |
| :So, for example, when the input grid is a generic 1x1 grid, and the output grid is the geos 1x1 grid:
| |
| IM=360
| |
| JM=180
| |
| INLON(1)=-180
| |
| INLON(361)=180
| |
| INSIN(1)=-1
| |
| INSIN(181)=1
| |
| LON2(1)=-180.5
| |
| LON2(361)=179.5
| |
| SIN2(1)=-1
| |
| SIN2(181)=1
| |
|
| |
| One last suggestion - if possible, change variable names <tt>LON2</tt> and <tt>SIN2</tt> to <tt>OUTLON</tt> and <tt>OUTSIN</tt>.
| |
|
| |
| --[[User:Melissa Payer|Melissa Payer]] 11:44, 8 June 2012 (EDT)
| |
|
| |
| === Bug fixes for regional emission masks ===
| |
|
| |
| <span style="color:green">'''''This update was tested in the 1-month benchmark simulation [[GEOS-Chem v9-01-03 benchmark history#v9-01-03q|v9-01-03q]].'''''</span>
| |
|
| |
| <span style="color:red">'''''This code was removed from [[GEOS-Chem v10-01]] and newer versions. EMEP ship emissions are now implemented through the [[HEMCO|HEMCO emissions component]].'''''</span>
| |
|
| |
| '''''[mailto:ltmurray@post.harvard.edu Lee Murray] wrote:'''''
| |
|
| |
| :I believe that there is a bug in the masks for the regional overwrites. I added the following code to the end of <tt>READ_EUROPE_MASK</tt> in <tt>emep_mod.F</tt>:
| |
|
| |
| ! Plot map
| |
| DO J=JJPAR,1,-1
| |
| WRITE(6,'(10F5.1)') EUROPE_MASK(30:39,J)
| |
| ENDDO
| |
|
| |
| :to output a map of EUROPE_MASK values, after it was re-gridded from 1x1 to the GEOS5 4x5 grid.
| |
|
| |
| - READ_EUROPE_MASK: Reading /Data/input/geos/GEOS_1x1/EMEP_200911/EMEP_mask.geos.1x1
| |
| ... etc ...
| |
| 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
| |
| 0.0 7.5 12.5 12.5 12.5 12.5 12.5 12.5 12.5 12.5
| |
| 0.0 12.0 20.0 20.0 20.0 20.0 20.0 20.0 20.0 20.0
| |
| 0.0 12.0 20.0 20.0 20.0 20.0 20.0 20.0 20.0 20.0
| |
| 0.0 12.0 20.0 20.0 20.0 20.0 20.0 20.0 20.0 20.0
| |
| 0.0 12.0 20.0 20.0 20.0 20.0 20.0 20.0 20.0 20.0
| |
| 0.0 12.0 20.0 20.0 20.0 20.0 20.0 20.0 20.0 20.0
| |
| 0.0 12.0 20.0 20.0 20.0 20.0 20.0 20.0 20.0 20.0
| |
| 0.0 10.5 19.5 20.0 20.0 20.0 20.0 20.0 20.0 20.0
| |
| 0.0 0.0 2.0 3.5 9.5 15.5 20.0 20.0 20.0 20.0
| |
| 0.0 0.0 0.0 0.0 0.0 0.0 6.5 17.5 20.0 20.0
| |
| 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 9.5 12.5
| |
| ... etc ...
| |
|
| |
| :You can see that many of the values exceed 1.0, so in many grid boxes, the emission rates for Europe would be multiplied by a factor of 20 when the mask is applied. I have also tested and get similar results for the Canada and US masks. I imagine the same holds for China and Mexico.
| |
|
| |
| :We can get masks in the correct range of 0 to 1 if we use
| |
|
| |
| CALL DO_REGRID_A2A( LLFILENAME, I1x1, J1x1,
| |
| & INGRID, OUTGRID, PERAREA=0 )
| |
| EUROPE_MASK = OUTGRID
| |
|
| |
| :instead of
| |
|
| |
| CALL DO_REGRID_A2A( LLFILENAME, I1x1, J1x1,
| |
| & INGRID, OUTGRID, PERAREA=1 )
| |
| EUROPE_MASK = OUTGRID
| |
|
| |
| :I think that this is because although the masks themselves are unitless, they are actually m2 region / m2 box, they are already in a "per unit area" and don't need to be converted in DO_REGRID_A2A.
| |
|
| |
| :Also, I suspect that once they are re-gridded, all values greater than 0 need to be set to 1 to be consistent with how the masks are applied to the regridded emissions. But someone who is more familiar with the emissions should verify that.
| |
|
| |
| --[[User:Bmy|Bob Y.]] 15:14, 20 December 2012 (EST)
| |
|
| |
| === Updates to improve performance ===
| |
|
| |
| <span style="color:green">'''''This update was tested in the 1-month benchmark simulation [[GEOS-Chem v9-01-03 benchmark history#v9-01-03r|v9-01-03r]] and approved on 28 Aug 2012.'''''</span>
| |
|
| |
| In [[GEOS-Chem v9-01-03]], we have added the following updates to the MAP_A2A regridding scheme:
| |
|
| |
| #<p>The driver routine <tt>DO_REGRID_A2A</tt> now reads the grid input information (i.e. longitude edges and sines of latitude edges) from netCDF files. This replaces the original implementation, which read from ASCII files.</p>
| |
| #<p>Important DO loops in routines <tt>DO_REGRID_A2A</tt>, <tt>XMAP</tt> and <tt>YMAP</tt> (all in <tt>GeosUtil/regrid_a2a_mod.F90</tt>) have now been parallelized with [http://www.openmp.org OpenMP directives].</p>
| |
| #<p>The <tt>QSUM</tt> and <tt>SUM</tt> variables in routines <tt>XMAP</tt> and <tt>YMAP</tt> have been changed from <tt>REAL*4</tt> to <tt>REAL*8</tt>. This eliminates some minor numerical noise past the 6th decimal place.</p>
| |
|
| |
| --[[User:Bmy|Bob Y.]] 15:33, 28 August 2012 (EDT)
| |
|
| |
| === Fix for regridding error in offline CO2 simulation ===
| |
|
| |
| <span style="color:green">'''''We have corrected this issue as a post-release patch on 22 Sep 2012. This issue does not affect the full-chemistry simulations, only the CO2 simulation.'''''</span>
| |
|
| |
| '''''[mailto:RayNassar@ec.gc.ca Ray Nassar] wrote:'''''
| |
|
| |
| :There is a bug in v9-01-03 for the CO2 simulation that causes a crash. It relates to the MAP_A2A regridding and the fact that the net terrestrial exchange file is named with <tt>geos.1x1</tt> (360x181) although it is actually generic 1x1 (360x180) and it is trying to be read as <tt>geos.1x1</tt>.
| |
|
| |
| Please see [[CO2_simulation#Fix_for_regridding_error_in_offline_CO2_simulation|this wiki post on our ''CO2 simulation'' wiki page]] for a full description of the issue.
| |
|
| |
| --[[User:Bmy|Bob Y.]] 15:21, 20 December 2012 (EST)
| |
|
| |
| === Bug in regridding of anthropogenic emissions ===
| |
|
| |
| <span style="color:green">'''''This issue was resolved in [[GEOS-Chem v9-02 benchmark history#v9-02e|GEOS-Chem v9-02e]], which was approved on 09 Jan 2013.'''''</span>
| |
|
| |
| <span style="color:red">'''''NOTE: The code modules listed below were removed in [[GEOS-Chem v10-01]]. All anthropogenic emissions are now handled via the [[HEMCO|HEMCO emissions component]].'''''</span>
| |
|
| |
| It appears that there is a bug in the regridding of the NEI2005 anthropogenic emissions and the David Streets Anthropogenic Emissions. In some instances, the regridding is assuming that the emissions data have units of per unit area (i.e. per m<sup>2</sup> or cm<sup>2</sup>) when they do not (or vice versa). We will rectify this situation shortly.
| |
|
| |
| '''''[mailto:ckeller@seas.harvard.edu Christoph Keller] wrote''''':
| |
|
| |
| :Jenny Fisher asked me why the <tt>PERAREA</tt> flag is switched on in all the <tt>READ_STREETS</tt> calls in <tt>GeosCore/streets_anthro_mod.F</tt> except for NOx and CO? Those flags determine whether the fields are normalized by their area during the regridding or not, and in my opinion the flag should be switched on (I.e. PERAREA=1) for all compounds since all emissions are in kg/yr.
| |
|
| |
| :I'm concerned that we have the wrong regridding for STREETS NOx and CO emissions. E.g. On line 686 we have:
| |
|
| |
| CALL READ_STREETS( FILENAME, 'ANTHSRCE', 1,
| |
| & TAU, TEMP, PERAREA=0 )
| |
|
| |
| :However, these are emissions in kg/yr, so I think we should have:
| |
|
| |
| CALL READ_STREETS( FILENAME, 'ANTHSRCE', 1,
| |
| & TAU, TEMP, PERAREA=1 )
| |
|
| |
| :This only affects CO and NOx, the other species use the PERAREA=1 flag.
| |
|
| |
| :I have made a test run where I wrote out the total Streets emissions before and after regridding, this is what I get:
| |
|
| |
| Original version (PERAREA = 0):
| |
|
| |
| - EMISS_STREETS_ANTHRO: Reading /as/data/geos/GEOS_1x1/Streets_200812/Streets_NOx_ind_2006.generic.1x1
| |
| STREETS sum before regridding (tracer 1): 8483151381.67771
| |
| STREETS sum after regridding (tracer 1): 423913738.256104
| |
|
| |
| With PERAREA = 1:
| |
|
| |
| - EMISS_STREETS_ANTHRO: Reading /as/data/geos/GEOS_1x1/Streets_200812/Streets_NOx_ind_2006.generic.1x1
| |
| STREETS sum before regridding (tracer 1): 8483151381.67771
| |
| STREETS sum after regridding (tracer 1): 8483151381.67771
| |
|
| |
| :I would say the latter makes much more sense! I think this is a bug which has been introduced with the new regridding package in v9-01-03, so this is only very recent.
| |
|
| |
| :The PERAREA flag can be a bit confusing, though, as PERAREA can be interpreted as "is the field I want to regrid in units of per area", whereas in fact it means "do you want to convert to per area before regridding".
| |
|
| |
| '''''[mailto:mpayer@seas.harvard.edu Melissa Payer] wrote:'''''
| |
|
| |
| :After looking carefully through <tt>GeosCore/streets_anthro_mod.F</tt>, we found there was a second error in the code for regridding of the monthly variability factors which essentially negates the bug in NOx emissions.
| |
|
| |
| :For regridding of NOx emissions, we call routine <tt>READ_STREETS</tt> with:
| |
|
| |
| ! File name
| |
| IF ( IS_2006 ) THEN
| |
|
| |
| FILENAME = TRIM( STREETS_DIR ) // 'Streets_NOx_'//
| |
| & SOURCES( NSOURCE ) // '_2006.generic.1x1'
| |
|
| |
| TAU = TAU2006
| |
| ELSE
| |
| FILENAME = TRIM( STREETS_DIR ) //
| |
| & 'Streets_NOx_FF_2000.generic.1x1'
| |
|
| |
| TAU = TAU2000
| |
|
| |
| ENDIF
| |
|
| |
| ! Read data
| |
| CALL READ_STREETS( FILENAME, 'ANTHSRCE', 1,
| |
| & TAU, TEMP, PERAREA=0 )
| |
|
| |
| :As Christoph pointed out, this should be:
| |
|
| |
| ! Read data
| |
| CALL READ_STREETS( FILENAME, 'ANTHSRCE', 1,
| |
| & TAU, TEMP, PERAREA=1 )
| |
|
| |
| :We then apply seasonal variation to both NOx and CO:
| |
|
| |
| ! Seasonal Variation for NOx
| |
| !--------------------------
| |
| ! Monthly Variability for any year. Variability has
| |
| ! been obtained from 2004 1x1 data, for two cases:
| |
| ! FF seasonality for 2000, and TOTAL (BF+FF) for 2006
| |
| ! inventories (phs, 12/2/08)
| |
|
| |
| IF ( IS_2006 ) THEN
| |
| FILENAME = TRIM( STREETS_DIR ) //
| |
| & 'Streets_2004_NOx_MonthFctr_total.generic.1x1'
| |
|
| |
| ELSE
| |
|
| |
| ! redefine the entire path here
| |
| FILENAME = TRIM( DATA_DIR_1x1 ) // 'Streets_200812/'
| |
| & // 'Streets_2004_NOx_MonthFctr_FF.generic.1x1'
| |
|
| |
| ENDIF
| |
|
| |
| CALL READ_STREETS( FILENAME, 'RATIO-2D', 71,
| |
| $ TAUMONTH_2004, TEMP, PERAREA=1 )
| |
|
| |
| NOX = NOX * TEMP
| |
|
| |
| :However, since the monthly scale factors are unitless, the regridding routine should actually be called with
| |
|
| |
| CALL READ_STREETS( FILENAME, 'RATIO-2D', 71,
| |
| $ TAUMONTH_2004, TEMP, PERAREA=0 )
| |
|
| |
| :When we fix both of these bugs, they essentially cancel one another in the NOx totals. This is the reason that we did not see the bug in the 1-month benchmark output. The attached plots show that with both bug fixes applied surface NOx changes by a maximum of +/- 0.0522 ppb over Asia.
| |
|
| |
| :If we only fix the bug in regridding of NOx emissions, we get an increase of up to 300 ppb over China, which is clearly unrealistic. It is therefore obvious that we need to fix both the regridding of NOx emissions and of monthly scale factors.
| |
|
| |
| :These bugs go as far back as [[GEOS-Chem v9-01-03 benchmark history#v9-01-03k|GEOS-Chem v9-01-03k]], when the MAP_A2A regridding algorithm was implemented. However, we did not pick up on it because there was very little change in NOx emissions due to the cancelling effect of the two bugs. The source of the error likely stems from the confusion of the argument <tt>PERAREA</tt>.
| |
|
| |
| :* <tt>PERAREA = 1</tt> is used for units of mass (kg/yr, etc.). The array will be divided by unit area before regridding.
| |
| :* <tt>PERAREA = 0</tt> is used for units of concentration (unitless, v/v, molec/cm2/s, etc).
| |
|
| |
| :We will fix the bugs in <tt>GeosCore/streets_anthro_mod.F</tt> in the next 1-month benchmark ([[GEOS-Chem_v9-02_benchmark_history#v9-02e|GEOS-Chem v9-02e]]). I will also take another thorough look through all of the regridding calls in the code to make sure that we consistently and properly call the regridding routines. This will likely involve renaming the argument <tt>PERAREA</tt> throughout the code to avoid confusion as to what it represents.
| |
|
| |
| --[[User:Bmy|Bob Y.]] 16:11, 19 December 2012 (EST)
| |
|
| |
| === Additional bug fixes for MAP_A2A regridding algorithm ===
| |
|
| |
| <span style="color:green">'''''NOTE: This issue was resolved in [[GEOS-Chem v9-02 benchmark history#v9-02e|GEOS-Chem v9-02e]], which was approved on 09 Jan 2013. This update is included in Adjoint [[GEOS-Chem_Adjoint_v35 | v35j]].'''''</span>
| |
|
| |
| The source code was carefully combed to ensure that <tt>PERAREA</tt> is assigned the proper value when calling the MAP_A2A regridding routines. The following bugs were found and corrected:
| |
|
| |
| #In subroutine <tt>SEASCL_EDGAR_SHIP_SO2 (edgar_mod.F)</tt>, <tt>PERAREA</tt> should be 0 because scale factors are unitless.
| |
| #In subroutine <tt>GET_AEF_05x0666 (megan_mod.F)</tt>, <tt>PERAREA</tt> should be 0 for all calls to <tt>DO_REGRID_A2A</tt> because AEFs are in μg C/m2/hr and do not need to be converted to units per area. This fix ensures that the nested grid code is consistent with the code in <tt>GET_AEF</tt> for global grids.
| |
| #In subroutine <tt>GET_NEI99_WKSCALE_05x0666 (nei2005_anthro_mod.F)</tt>, <tt>PERAREA</tt> should be 0 because scale factors are unitless.
| |
| #In subroutine <tt>READ_BIOFUEL_SO2 (sulfate_mod.F)</tt>, <tt>PERAREA</tt> should be 1 for historical emissions because data are in kg/yr.
| |
|
| |
| To avoid confusion in assigning a value to <tt>PERAREA</tt>, we have renamed the parameter to <tt>IS_MASS</tt> throughout the standard source code, where:
| |
|
| |
| *<tt>IS_MASS = 1</tt> is used for units of mass (kg/yr, etc.). The array will be divided by unit area before regridding.
| |
| *<tt>IS_MASS = 0</tt> is used for units of concentration (unitless, v/v, molec/cm2/s, etc).
| |
|
| |
| --[[User:Melissa Payer|Melissa Payer]] 14:44, 7 January 2013 (EST)
| |
|
| |
| === Bug in grid_mod.F90 ===
| |
|
| |
| <span style="color:green">'''''This update was tested in the 1-month benchmark simulation [[GEOS-Chem_v9-02_benchmark_history#v9-02l|v9-02l]] and approved on 26 Jun 2013. This update is included in Adjoint [[GEOS-Chem_Adjoint_v35 | v35a]].'''''</span>
| |
|
| |
| '''''[mailto:ckeller@seas.harvard.edu Christoph Keller] wrote''''':
| |
|
| |
| :I found a small bug in grid_mod.F: The sine of the last latitude edge (i.e. the north pole) is not calculated, which can cause problems for the regridding.
| |
|
| |
| :I just <span style="color:green">added the following lines</span> to the (last) latitude edge calculation (COMPUTE_GRID in grid_mod.F, L235ff):
| |
|
| |
| !%%%%%% LATITUDE EDGES (the last edge) %%%%%%%%%%%%%%%%%%%%%%%%
| |
|
| |
| ! Test for North Pole
| |
| IF ( J2 == JNP ) THEN
| |
|
| |
| ! North pole case (global grids only)
| |
| DO I = I1, I2
| |
| YEDGE (I,JNP+1,L) = +90d0
| |
| YEDGE_R(I,JNP+1,L) = YEDGE(I,JNP+1,L) * PI_180
| |
|
| |
| <span style="color:green">! Also compute sine of last latitude edge! (ckeller, 02/13/12)
| |
| YSIN(I,J2+1,L) = SIN ( YEDGE_R(I,J2+1,L) )</span>
| |
| ENDDO
| |
|
| |
| ELSE
| |
|
| |
| ! No north pole (nested grids only)
| |
| DO I = I1, I2
| |
| YEDGE (I,J2+1,L) = YEDGE(I,J2,L ) + DLAT(I,J2,L)
| |
| YEDGE_R(I,J2+1,L) = YEDGE(I,J2+1,L) * PI_180
| |
|
| |
| <span style="color:green">'''! Also compute sine of last latitude edge! (ckeller, 02/13/12)
| |
| YSIN(I,J2+1,L) = SIN ( YEDGE_R(I,J2+1,L) )'''</span>
| |
| ENDDO
| |
|
| |
| ENDIF
| |
|
| |
| :Looks like this fixes the problem.
| |
|
| |
| --[[User:Melissa Payer|Melissa Payer]] 16:41, 13 February 2013 (EST)
| |
|
| |
| ==== Update: error observed in v10-01h ====
| |
|
| |
| The fix described here caused an [[#Fix_for_optimization_error_discovered_in_v10-01h|an apparent optimization error]] that ended up inadvertently zeroing the <tt>YSIN</tt> array. This was noticed in [[GEOS-Chem v10-01 benchmark history#v10-01h|GEOS-Chem v10-01h]]. Please see [[#Fix_for_optimization_error_discovered_in_v10-01h|the following section]] for more information about how this issue was resolved.
| |
|
| |
| --[[User:Bmy|Bob Y.]] 17:46, 1 April 2015 (EDT)
| |
|
| |
| === Fix for optimization error discovered in v10-01h ===
| |
|
| |
| <span style="color:green">'''''This bug was fixed for the GEOS-Chem v10-01 public release.'''''</span>
| |
|
| |
| We discovered this apparent optmization error in GEOS-Chem v10-01h. This error was caused the modification labeled [[#Bug in grid_mod.F90|'''''Bug in grid_mod.F90''''']] as described above. This error has now been corrected.
| |
|
| |
| '''''[[User:Bmy|Bob Yantosca]] wrote:'''''
| |
|
| |
| :In our validation of [[GEOS-Chem v10-01 benchmark history#v10-01h|GEOS-Chem v10-01h]], we discovered that the <tt>geos5_4x5_fullchem</tt> simulation was consistently dying with a segmentation fault. Oddly enough, this only occurred when <tt>DEBUG=n</tt> (i.e. when all debugging features were turned off, and when the normal ifort <tt>-O2</tt> optimization was used. This was happening when HEMCO was trying to read in the NEI2011 mask file.
| |
|
| |
| :I traced this to an error in the regridding. I put some debug print statements in the routine <tt>MAP_A2A_R4R4</tt> (in module <tt>GeosUitl/regrid_a2a_mod.F90</tt>) as follows:
| |
|
| |
| print*, '### size sin1: ', size(sin1)
| |
| print*, '### data sin1: ', sin1
| |
| print*, '### size sin2: ', size(sin2)
| |
| print*, '### data sin2: ', sin2
| |
| call flush(6)
| |
|
| |
| ! Otherwise, call YMAP to regrid in the N-S direction
| |
| CALL ymap_r4r4(in, jm, sin1, qtmp(1,1+ig), jn, sin2, q2(1,1+ig), ig, iv)
| |
|
| |
| :Here, the variable <tt>SIN1</tt> holds the sines of the grid box latitude edges on the input grid, and <tt>SIN2</tt> holds the sines of the grid box latitude edges on the output grid. These quantities are needed for the regridding. The debug ouptut revealed that <tt>SIN2</tt> was being set to incorrect values (i.e. mostly zeroes):
| |
|
| |
| HEMCO: Opening /as/data/geos/ExtData/HEMCO/MASKS/v2014-07/USA_LANDMASK_NEI2011_0.1x0.1.nc
| |
| ### size sin1: 401
| |
| ### data sin1: 0.3420201 0.3436597 0.3452982 0.3469357
| |
| 0.3485720 0.3502074 0.3518416 0.3534749 0.3551069
| |
| 0.3567380 0.3583679 0.3599968 0.3616246 0.3632512
| |
| 0.3648767 0.3665012 0.3681245 0.3697468 0.3713678
| |
| ... etc ...
| |
| ### size sin2: 47
| |
| ### data sin2: -1.000000 0.0000000E+00 0.0000000E+00 0.0000000E+00
| |
| 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00
| |
| 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00
| |
| 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00
| |
| 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00
| |
| 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00
| |
| 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00
| |
| 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00
| |
| 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00
| |
| 0.0000000E+00 0.0000000E+00 0.0000000E+00
| |
|
| |
| :I eventually traced the error to these lines of code in <tt>GeosUtil/grid_mod.F90</tt>.
| |
|
| |
| ! Lat edges (radians)
| |
| YDGR(I,J,L) = ( PI_180 * YDG(I,J,L) )
| |
|
| |
| ! mjc - Compute sine of latitude edges (needed for map_a2a regrid
| |
| YSN(I,J,L) = SIN ( YDGR(I,J,L) )
| |
|
| |
| :For some reason, the <tt>YSN(I,J,L)</tt> values were all being set to zero. '''This was not due to a bug in the Fortran code''' (as the lines above are legal Fortran-90 statements). We instead attribute this to an '''optimization error''', perhaps caused by a unforeseen shortcoming of the ifort compiler.
| |
|
| |
| :The fix for this optimization error is simple. I rewrote the code above so that it now looks like this:
| |
|
| |
| <span style="color:red">!------------------------------------------------------------------------------
| |
| ! Prior to 3/26/15:
| |
| ! Fix apparent optimization error. Now pass a scalar to the SIN function.
| |
| ! Also save the result of SIN in a scalar before storing in the YSN array.
| |
| ! Don't know why this happens but this seems to fix it. (bmy, 3/26/15)
| |
| ! YSN(I,J,L) = SIN ( YDGR(I,J,L) )
| |
| !------------------------------------------------------------------------------</span>
| |
| <span style="color:green">YEDGE_VAL = YDGR(I,J,L) ! Lat edge in radians
| |
| YSIN_VAL = SIN( YEDGE_VAL ) ! SIN( lat edge )
| |
| YSN(I,J,L) = YSIN_VAL ! Store in YSN array</span>
| |
|
| |
| :I also had to make this same fix in 3 other locations within <tt>GeosUtil/grid_mod.F90</tt>.
| |
|
| |
| :The fix is this: instead of passing the value of <tt>YDGR(I,J,L)</tt> directly to the Fortran <tt>SIN()</tt> function, I first saved that to a scalar variable, and then passed that to <tt>SIN()</tt>. The output of <tt>SIN()</tt> is then saved to a scalar variable first, before it is used to define the value of <tt>YSN(I,J,L)</tt>. This seems to be more robust.
| |
|
| |
| :In addition, I added some extra printout to the GEOS-Chem log file, so that we can quickly determine if this error is happening again. Good output should look like this:
| |
|
| |
| Grid box latitude edges [degrees]:
| |
| -90.000 -88.000 -84.000 -80.000 -76.000 -72.000 -68.000 -64.000
| |
| -60.000 -56.000 -52.000 -48.000 -44.000 -40.000 -36.000 -32.000
| |
| -28.000 -24.000 -20.000 -16.000 -12.000 -8.000 -4.000 0.000
| |
| 4.000 8.000 12.000 16.000 20.000 24.000 28.000 32.000
| |
| 36.000 40.000 44.000 48.000 52.000 56.000 60.000 64.000
| |
| 68.000 72.000 76.000 80.000 84.000 88.000 90.000
| |
|
| |
| SIN( grid box latitude edges )
| |
| -1.000 -0.999 -0.995 -0.985 -0.970 -0.951 -0.927 -0.899
| |
| -0.866 -0.829 -0.788 -0.743 -0.695 -0.643 -0.588 -0.530
| |
| -0.469 -0.407 -0.342 -0.276 -0.208 -0.139 -0.070 0.000
| |
| 0.070 0.139 0.208 0.276 0.342 0.407 0.469 0.530
| |
| 0.588 0.643 0.695 0.743 0.788 0.829 0.866 0.899
| |
| 0.927 0.951 0.970 0.985 0.995 0.999 1.000
| |
|
| |
| :This error is now fixed in GEOS-Chem v10-01h.
| |
|
| |
| --[[User:Bmy|Bob Y.]] 12:29, 27 March 2015 (EDT)
| |
|
| |
| === Error in regridding regional data files ===
| |
|
| |
| '''''[[User:Bmy|Bob Yantosca]] wrote:'''''
| |
|
| |
| <blockquote>Melissa and I have noticed a problem that was flagged by the [[GEOS-Chem Unit Tester]] in GEOS-Chem versions [[GEOS-Chem v10-01 benchmark history#v10-01h|v10-01h]] and [[GEOS-Chem v10-01 benchmark history#v10-01i|v10-01i]]. It seems that there are some negatives in the new EMEP emissions files (these are regional data) that we have gotten from Aaron van Donkelaar. We’ve been doing a few test simulations to try to isolate the problem.</blockquote>
| |
|
| |
| <blockquote>The output from the HEMCO log shows negatives in the array containing NO from the EMEP emissions file:</blockquote>
| |
|
| |
| Container EMEP_NO
| |
| . . .
| |
| -->Source file : $ROOT/EMEP/v2015-03/EMEP.generic.1x1.nc
| |
| . . .
| |
| -->Source parameter: NO
| |
| . . .
| |
| '''-->Array sum : -1.1036439E+12'''
| |
| '''-->Array min & max : -2.2052954E+10 1.3323866E-10'''
| |
| ... etc ...
| |
|
| |
| <blockquote>Which leads to the following error:</blockquote>
| |
|
| |
| !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
| |
| HEMCO ERROR: Negative emissions in: EMEP_NO. To allow negatives, edit settings
| |
| in the configuration file.
| |
| ERROR LOCATION: HCO_CalcEmis (HCO_CALC_MOD.F90)
| |
| ERROR LOCATION: HCO_RUN (hco_driver_mod.F90)
| |
| !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
| |
|
| |
| '''''[[User:Melissa Payer|Melissa Sulprizio]] wrote:'''''
| |
|
| |
| <blockquote>I did a test, compiling with <tt>DEBUG=yes FPE=yes</tt> and using the prior emissions in the <tt>HEMCO/EMEP/v2014-07/</tt> directory. In this case, we don't get negatives for the EMEP emissions anymore, but we still see negatives from the MIX and NEI1 regional emissions inventories. I think that confirms that the issue has to do with reading and regridding data from regional files as opposed to global files.</blockquote>
| |
|
| |
| '''''[[User:Christoph Keller|Christoph Keller]] wrote:'''''
| |
|
| |
| <blockquote>I found the problem with the MAP_A2A regridding (in module <tt>GeosUtil/regrid_a2a_mod.F90</tt>).
| |
| The variable <tt>qsum</tt> in the ymap routines (i.e. <tt>YMAP_R8R8</tt>, <tt>YMAP_R4R8</tt>, <tt>YMAP_R4R4</tt>, <tt>YMAP_R8R4</tt>), was not being initialized, which resulted in undefined variables being passed to the output array <tt>q2</tt>. This was only a problem when regridding from a nested grid onto a global grid.</blockquote>
| |
|
| |
| <blockquote>I added a bug fix by explicitly initializing <tt>qsum</tt> within the parallel loop in the ymap routines (see the code in '''bold face'''):</blockquote>
| |
|
| |
| !$OMP PARALLEL DO &
| |
| !$OMP DEFAULT( SHARED ) &
| |
| !$OMP PRIVATE( I, J0, J, M, QSUM, MM, DY )
| |
| do 1000 i=1,im
| |
| '''qsum = 0.0d0'''
| |
| ... etc ...
| |
|
| |
| <blockquote>This fixed the problem.</blockquote>
| |
|
| |
| --[[User:Bmy|Bob Y.]] 17:35, 1 April 2015 (EDT)
| |
|
| |
| == Outstanding issues that have not been resolved ==
| |