MCIP v5.1 floating point error in resistcalc.f90

a.kashfi73 · January 22, 2020, 6:54pm

Dear all,

I get following error after running MCIPv5.1.

Program received signal SIGFPE: Floating-point exception - erroneous arithmetic operation.

Backtrace for this error:
#0 0x7f019ce912da in ???
#1 0x7f019ce90503 in ???
#2 0x7f019c2d7f1f in ???
#3 0x55e5854f669b in resistcalc_
at /home/a13/NAS/a13/mcip/src/resistcalc.f90:281
#4 0x55e5854cc2bc in pblsup_
at /home/a13/NAS/a13/mcip/src/pblsup.f90:331
#5 0x55e585477edd in dynflds_
at /home/a13/NAS/a13/mcip/src/dynflds.f90:67
#6 0x55e58543952c in mcip
at /home/a13/NAS/a13/mcip/src/mcip.f90:169
#7 0x55e585439833 in main
at /home/a13/NAS/a13/mcip/src/mcip.f90:79
Floating point exception (core dumped)
Error running mcip

tlspero · January 22, 2020, 7:28pm

@a.kashfi73 –

The traceback shows that MCIP is crashing on line 281 in resistcalc.f90 where there is a division by leaf area index, which I suspect is set to 0.0 when MCIP is crashing. LAI could be 0.0 over water, but the block where LAI is calculated in MCIP should only be activated on land points. Can you check the values of XLAI and XLWMASK for the point where MCIP is crashing?

–Tanya

a.kashfi73 · January 22, 2020, 8:35pm

@tlspero
Dear Tanya,

I thank you for your quick response. I am sorry, but can you please tell me how I can check that?
I just saw LAI and LANDMASK in my WRF output, the values were fine to me(zero at sea, and not zero on land). Can you please tell me what you mean by the point where MCIP is crashing?

Thank you again,

mmallard · January 22, 2020, 8:55pm

Hi,
I think she’s referring the grid cell where mcip is dividing by zero and then crashing. One way to do it would be to add some write statements to resistcalc.f90 outputting i,j and the lai and landmask values just above line 281 (where it crashes) and then recompile and rerun

a.kashfi73 · January 22, 2020, 10:37pm

@mmallard @tlspero

Thank you for your suggestion Megan. Code is now as follows.

print*,xlai(c,r)
print*,xlwmask(c,r)
radf = 1.1 * xrgrnd(c,r) / ( radl * xlai(c,r) ) ! NP89 - EQN34
f1 = (rstmin / rsmax + radf) / (1.0 + radf)

With results as follows.

0.283726513
1.00000000
0.310195982
1.00000000
0.337951839
1.00000000
0.366095901
1.00000000
0.394635946
1.00000000

Program received signal SIGFPE: Floating-point exception - erroneous arithmetic operation.

Backtrace for this error:
#0 0x7f2bcbed22da in ???
#1 0x7f2bcbed1503 in ???
#2 0x7f2bcb318f1f in ???
#3 0x559d436207c3 in resistcalc_
at /home/a13/NAS/a13/mcip/src/resistcalc.f90:282
#4 0x559d435f62bc in pblsup_
at /home/a13/NAS/a13/mcip/src/pblsup.f90:331
#5 0x559d435a1edd in dynflds_
at /home/a13/NAS/a13/mcip/src/dynflds.f90:67
#6 0x559d4356352c in mcip
at /home/a13/NAS/a13/mcip/src/mcip.f90:169
#7 0x559d43563833 in main
at /home/a13/NAS/a13/mcip/src/mcip.f90:79
Floating point exception (core dumped)
Error running mcip

Thank you for your time.

mmallard · January 23, 2020, 2:27pm

That’s puzzling… I expected that LAI would have been set to zero somehow and the error was caused by dividing by zero. But LAI was set to a reasonable value just before the crash.
This type of error could also be caused by trying to do arithmetic operations with really, really small or big numbers. I’d recommend outputting the other variables in the equation (xrgrnd and radl) for completeness to make sure that all terms have reasonable values.

a.kashfi73 · January 23, 2020, 10:44pm

@tlspero @mmallard

Thank you for your time and for your help. I just noticed I had made a mistake in my log file. Sorry for that. There is a LAI zero as follows:

The code is:
IF ( rstmin > 130.0 ) THEN
radl = 30.0 ! [W/m2]
ELSE
radl = 100.0 ! [W/m2]
ENDIF
print*,xlai(c,r)
print*,xlwmask(c,r)
print*,xrgrnd(c,r)
print*,c
print*,r
radf = 1.1 * xrgrnd(c,r) / ( radl * xlai(c,r) ) ! NP89 - EQN34
f1 = (rstmin / rsmax + radf) / (1.0 + radf)

The result is:
45
111
0.599254608
1.00000000
0.00000000
45
112
0.625336707
1.00000000
0.00000000
45
113
0.652132213
1.00000000
0.00000000
45
114
0.00000000
1.00000000
0.00000000
46
40

Program received signal SIGFPE: Floating-point exception - erroneous arithmetic operation.

Backtrace for this error:
#0 0x7f83c0e1f2da in ???
#1 0x7f83c0e1e503 in ???
#2 0x7f83c0265f1f in ???
#3 0x5632d6bed925 in resistcalc_
at /home/a13/NAS/a13/mcip/src/resistcalc.f90:286
#4 0x5632d6bc32bc in pblsup_
at /home/a13/NAS/a13/mcip/src/pblsup.f90:331
#5 0x5632d6b6eedd in dynflds_
at /home/a13/NAS/a13/mcip/src/dynflds.f90:67
#6 0x5632d6b3052c in mcip
at /home/a13/NAS/a13/mcip/src/mcip.f90:169
#7 0x5632d6b30833 in main
at /home/a13/NAS/a13/mcip/src/mcip.f90:79
Floating point exception (core dumped)
Error running mcip

I thank you again.

tlspero · January 23, 2020, 11:09pm

@a.kashfi73 –

What LSM are you using in WRF? Are you using an urban model in WRF?

–Tanya

a.kashfi73 · January 24, 2020, 12:00am

Dear Tanya,

I thank you for your quick response. My WRF sf_surface_physics option is layer thermal diffusion scheme. And I am not using urban model in WRF.

Thank you,

mmallard · January 24, 2020, 5:17pm

Is there any ice at this point? That would one way of having a zero LAI value over a land point.

a.kashfi73 · January 25, 2020, 9:59pm

I am not so sure. How can I check that? I only have SEAICE variable which is zero over my domain.
But what would be the solution for mcip anyways?

Thank you,

mmallard · January 27, 2020, 3:51pm

I was wondering if the land use category was a snow or ice category where LAI would be reduced. Just to clarify, even for cells in such categories, LAI should be set to a minimum value of 0.01. So this could be an odd bug on WRF’s end, or maybe there’s some confusion with the landsea mask. Either way, it’d be good to know the landuse category.

To find A copy of



the LU category for this grid cell, you can check the field LU_INDEX in the wrfout file to get the category number or you could modify the resistcalc code to output it from mcip (“lu” is used just above around line 200).  Then check the wrfout’s attributes to see which LU dataset was used. For example, for runs that used the USGS dataset, the attribute MMINLU = “USGS”.  If you look at WRF’s vegetation parameter table under the “USGS” heading, you’ll find the category numbers on the first column and the category names in the last column.  You’d need to check the corresponding stanzas of the table if you used a different dataset, such as NLCD.
 the vegparm table can be found here:
 href="https://github.com/yyr/wrf/blob/master/run/VEGPARM.TBL" target="_blank">github.com

 href="https://github.com/yyr/wrf/blob/master/run/VEGPARM.TBL" target="_blank">yyr/wrf/blob/master/run/VEGPARM.TBL class="lang-tbl">Vegetation Parameters RGL      HS      SNUP  MAXALB   LAIMIN  LAIMAX  EMISSMIN EMISSMAX ALBEDOMIN ALBEDOMAX  Z0MIN    Z0MAX    ZTOPV    ZBOTV' 200.,   999.,   999.0,   0.04,    46.,    1.00,   1.00,   .880,    .880,     .15,      .15,      .50,     .50,     0.00,    0.00,     'Urban and Built-Up Land' 40.,   100.,   36.25,   0.04,    66.,    1.56,   5.68,   .920,    .985,     .17,      .23,      .05,     .15,     0.50,    0.01,     'Dryland Cropland and Pasture' 40.,   100.,   36.25,   0.04,    66.,    1.56,   5.68,   .930,    .985,     .20,      .25,      .02,     .10,     0.50,    0.01,     'Irrigated Cropland and Pasture' 40.,   100.,   36.25,   0.04,    66.,    1.00,   4.50,   .920,    .985,     .18,      .23,      .05,     .15,     0.50,    0.01,     'Mixed Dryland/Irrigated Cropland and Pasture' 40.,   100.,   36.25,   0.04,    68.,    2.29,   4.29,   .920,    .980,     .18,      .23,      .05,     .14,     0.50,    0.01,     'Cropland/Grassland Mosaic' 70.,    65.,   44.14,   0.04,    60.,    2.00,   4.00,   .930,    .985,     .16,      .20,      .20,     .20,     0.50,    0.01,     'Cropland/Woodland Mosaic' 40.,   100.,   36.35,   0.04,    70.,    0.52,   2.90,   .920,    .960,     .19,      .23,      .10,     .12,     0.50,    0.01,     'Grassland' 300.,   100.,   42.00,   0.03,    60.,    0.50,   3.66,   .930,    .930,     .25,      .30,      .01,     .05,     0.50,    0.10,     'Shrubland' 170.,   100.,   39.18,  0.035,    65.,    0.60,   2.60,   .930,    .950,     .22,      .30,      .01,     .06,     0.50,    0.10,     'Mixed Shrubland/Grassland' 70.,    65.,   54.53,   0.04,    50.,    0.50,   3.66,   .920,    .920,     .20,      .20,      .15,     .15,     5.00,    0.10,     'Savanna' 100.,    30.,   54.53,   0.08,    58.,    1.85,   3.31,   .930,    .930,     .16,      .17,      .50,     .50,     20.0,    11.5,     'Deciduous Broadleaf Forest' 150.,    30.,   47.35,   0.08,    54.,    1.00,   5.16,   .930,    .940,     .14,      .15,      .50,     .50,     14.0,     7.0,     'Deciduous Needleleaf Forest' 150.,    30.,   41.69,   0.08,    35.,    3.08,   6.48,   .950,    .950,     .12,      .12,      .50,     .50,     35.0,     1.0,     'Evergreen Broadleaf Forest' 125.,    30.,   47.35,   0.08,    52.,    5.00,   6.40,   .950,    .950,     .12,      .12,      .50,     .50,     17.0,     8.5,     'Evergreen Needleleaf Forest' 125.,    30.,   51.93,   0.08,    53.,    2.80,   5.50,   .930,    .970,     .17,      .25,      .20,     .50,     18.0,    10.0,     'Mixed Forest' 100.,    30.,   51.75,   0.01,    70.,    0.01,   0.01,   .980,    .980,     .08,      .08,   0.0001,  0.0001,     0.00,    0.00,     'Water Bodies' 40.,   100.,   60.00,   0.01,    68.,    1.50,   5.65,   .950,    .950,     .14,      .14,      .20,     .20,     0.50,    0.01,     'Herbaceous Wetland' has been truncated. show original

a.kashfi73 · January 28, 2020, 1:10am

I thank you for your guidance and for your time. There code now is as follows:

print*,xlai(c,r)
print*,xlwmask(c,r)
print*,xrgrnd(c,r)
print*,c
print*,r
print*,lu
radf = 1.1 * xrgrnd(c,r) / ( radl * xlai(c,r) ) ! NP89 - EQN34

And result is:

0.652132213
1.00000000
0.00000000
45
114
12
0.00000000
1.00000000
0.00000000
46
40
7
MMINLU of WRF output is “MODIFIED_IGBP_MODIS_NOAH”, and for this option 7 is ‘Open Shrublands’.

Thank you again,

mmallard · January 28, 2020, 3:49pm

No problem!
LAI should be constrained to a minimum value of 0.6 for that land use category. By the way, have you verified that there is a zero LAI at this grid cell within the wrfout file that’s being read? I doubt that mcip is overwriting it with a zero, but it’s worth eliminating the possibility.
Also, what WRF version is this? You can tell be looking at the TITLE attribute in your wrfout files.

By the way, if it turns out that this is a single odd value coming in from WRF and you’d like a quick, hard-coded solution to get by it, you could overwrite the zero LAI at that point with a quick if statement. metvars2ctm.f90 (see line 406) plugs in hard-coded LAI values from arrays just in case they are not available to read in from WRF, so you could use those as a reference of what to hard-code it to (see the 7th entry in the “laimod” and “laimnmod” arrays further down).

Of course, this would be a problematic solution if there are several zero’d LAI values to contend with and this is just the first one encountered by the code. So it’d be good to verify where the zero LAI value at that grid cell comes from so we can continue diagnosing the problem.

a.kashfi73 · January 29, 2020, 12:24am

Hello again,

I do appreciate your help. Yes, there are also other points with zero LAI in WRF output. I am using WRFv4.1.3.
And for your suggestion, it also did not work, because the metvars2ctm.f90 does not execute to the line 406. It would execute line 377 first and skip the rest of the script.
ELSE IF ( iflai .AND. ( MAXVAL(lai) < smallnum ) ) THEN
xlai(:,: ) = lai(sc:ec,sr:er)
My fix is to just change the greater sign (>) to the less sign (<) and everything will be fine. Is that the right way do you think?
I just checked the metvars2ctm.f90 in the CMAQv5.2.1 package, the above IF condition is also different.
Please let me know how it goes.

Thank you,

tlspero · January 29, 2020, 7:42pm

@a.kashfi73 –

I would NOT advise you to follow your suggestion.

The MCIP variable “lai” is initialized to 0.0 in init_met.f90. If leaf-area index is available in the WRF output as netCDF variable “LAI” or “LAI_PX”, it will be read and loaded into MCIP variable “lai” or “lai_px” in rdwrfem.f90.

In metvars2ctm.f90, the MCIP internal variable “xlai” is filled on the transfer grid (hence, the “x” in the name of the variable) which is the MCIP output domain plus the boundaries. In the section of code that you identified in metvars2ctm.f90, MCIP is trying to determine how to fill “xlai”, whether it uses leaf-area index from a variable in the WRF output or whether it will use a look-up table based on land use classification.

In the solution that you propose, you are suggesting that leaf-area index might be available in the WRF output file, but you want to bypass it and use a lookup table instead. That is, “MAXVAL(lai)” will still be 0.0 if LAI was not read from the WRF output, but it will be larger than “smallnum” if LAI was in the WRF output and filled with realistic values. We only want to use the lookup table if LAI is not available. However, if LAI is available, we need to ensure that there are reasonable values in the field.

I need a few more pieces of information from you to determine whether MCIP has been properly coded to handle the 5-layer thermal diffusion scheme (which I cannot guarantee that it has).

I need to know whether or not LAI is present in your WRF output file and what kind of values it contains. Therefore, I’d like for you to first type “ncdump -h wrfout | grep LAI $1”, and then paste the result of that in your reply here. That will show whether or not LAI is even in your WRF output file (so we know whether we need to use a lookup table in MCIP or if we are getting LAI directly from WRF).
If LAI is in the file, I’d like for you to visualize the field and tell me the range of values.
I’d also like to see the header to your WRF file. Please type “ncdump -h wrfout” and post the results in your reply. You can either cut-and-paste to the reply or add an attachment.

Thank you. Hopefully we will get this resolved for you soon. I want to determine if this is specific to you or if there is a bug in MCIP for handling the 5-layer thermal diffusion scheme.
–Tanya

a.kashfi73 · January 30, 2020, 3:43pm

Thank you Tanya,

Please see the attachments.

header.txt (38.6 KB) LAI.txt (223 Bytes)

tlspero · January 31, 2020, 1:46pm

@a.kashfi73 –

Thank you for providing the information. I’m still looking into your issue, and I don’t see anything obvious, so I’ll need to do a little more digging. I may need to ask for some additional files, as well.

In the meantime, I’m curious as to why (scientifically) you’ve selected the 5-layer thermal diffusion instead of a more advanced LSM like Noah or P-X. Would you be willing to answer that?

–Tanya

NRiceKYDAQ · February 3, 2020, 1:03pm

Hi Tanya,

I would just like to add that I’m having the exact same error message after running MCIP V5.1 over the weekend. So, this isn’t specific to a.kashfi73. Please let me know if there is any information I can provide with regard to my run/error that might help find the cause of this bug.

Best,

Nathan

tlspero · February 3, 2020, 1:18pm

@NRiceKYDAQ –

I’m sorry that you’re also having the same issue. Can you confirm the setting you used for sf_surface_physics in WRF, as well as the land use classification in your runs?

–Tanya

Topic		Replies	Views
Float point Exception MCIP	3	679	September 5, 2019
MCIP crash error MCIP	13	151	April 6, 2024
Segmentation Fault running MCIP 5.0 MCIP	6	1716	October 9, 2019
A bug found in MCIP v5.1 MCIP	1	336	May 14, 2021
MCIP run: Segmentation fault (core dumped) MCIP	5	434	July 23, 2021

MCIP v5.1 floating point error in resistcalc.f90

Related Topics