Calculating change in area of PFTs in CLM5

[James King]

Hi all,

I am trying to calculate the change in forest area in various domains that CLM5 simulates under the various scenarios produced for CMIP6. I am using the file 'landuse.timeseries_0.9x1.25_SSP1-2.6_78pfts_CMIP6_simyr1850-2100_c190214.nc'.

Currently, I am doing this by multiplying grid cell area by PFT type and dividing by the land fraction, i.e.

AREA* PCT_NAT_PFT/100*LANDFRAC_PFT (for a given domain)

then calculating the area-weighted mean of the resulting field using CDO, adding together all the tree PFTs, and plotting. The shape of the curve looks right, but the area units don't and I am not sure why. Is this method correct or is there a better way to do this sort of thing?

Thanks,

James

 

[lawrencepj1]

Hi James,

This is the correct method for calculating areas from the landuse time series file. The area is in km^2 so that may be the problem. Today is a holiday in the US for 4th of July so sorry on the slow response.

Peter

 

[James King]

Hi Peter,

Thanks! IDL and Python scripts that allegedly do the same thing are giving different values for the area. Good to know I'm on the right track though, and apologies for not remembering the date!

James

 

[lawrencepj1]

Hi James

Sorry I just remembered you need to calculate the PCT_NATVEG as part of this calculation. As this is not on the file you need to take this as PCT_NATVEG = 100.0 - PCT_CROP (for the year from the landuse timeseries file) - PCT_GLACIER - PCT_LAKE - PCT_WETLAND (all on the surface data file).

So the area of a PFT should be = LANDFRAC_PFT*AREA* PCT_NATVEG[year] * PCT_NAT_PFT[year]/100

The LANDFRAC_PFT is simply the landfrac calculated for the land model which can be different to the land frac calculated for the ocean model (different resolutions and practicalities of modeling the ocean).

cheers
Peter

 

[lawrencepj1]

So the area of a PFT should be = LANDFRAC_PFT*AREA* PCT_NATVEG[year] / 100 * PCT_NAT_PFT[year]/100

forgot the second / 100 ;-)

 

[James King]

Thanks for this Peter. Is the correct surface data file for this run the one from 1850, 2000, or 2014? I.e.

surfdata_0.9x1.25_hist_78pfts_CMIP6_simyr1850_c190214.nc
surfdata_0.9x1.25_hist_78pfts_CMIP6_simyr2000_c190214.nc
surfdata_0.9x1.25_hist_78pfts_CMIP6_simyr2014_c190516.nc

 

[lawrencepj1]

Hi James

For your analysis the surface data for the glacier, lake and wetland should all be the same. Just make sure you have the same surface data file date as the time series. ie c190214 dates should match as well as the resolution etc.

Peter

 

[James King]

Hi Peter,

Having done this calculation, the results for North America still seem to have dubious units for area - not sure this is km^2...

 

[lawrencepj1]

Hi James

Not sure how you are doing your regional subsetting but here is python code for global area of all trees in km2

print("Reading File: " + surfacefilename)
surfacefile = netcdf4.Dataset(surfacefilename,'r')

inlandfrac=np.asfarray(surfacefile.variables["LANDFRAC_PFT"][:,:],np.float64)
inlandmask=np.asfarray(surfacefile.variables["PFTDATA_MASK"][:,:],int)
inarea=np.asfarray(surfacefile.variables["AREA"][:,:],np.float64)
inlandarea = inlandfrac * inarea
inlonxy = np.asfarray(surfacefile.variables["LONGXY"][:,:],np.float64)
inlon = inlonxy[0,:]
inlatxy = np.asfarray(surfacefile.variables["LATIXY"][:,:],np.float64)
inlat = inlatxy[:,0]
ninlat = inlat.size
ninlon = inlon.size
inpctglacier = np.asfarray(surfacefile.variables["PCT_GLACIER"][:,:],np.float64)
inpctlake = np.asfarray(surfacefile.variables["PCT_LAKE"][:,:],np.float64)
inpctwetland = np.asfarray(surfacefile.variables["PCT_WETLAND"][:,:],np.float64)
inpcturbansub = np.asfarray(surfacefile.variables["PCT_URBAN"][:,:,:],np.float64)
inpcturban = inpcturbansub.sum(axis=0)

surfacefile.close()

print("Reading File: " + landusefilename)
landusefile = netcdf4.Dataset(landusefilename,'r')

# variableunits = inputfile.variables[pftname].units
inpctcrop = np.asfarray(landusefile.variables["PCT_CROP"][startyearindex:endyearindex+1,:,:],np.float64)
variabletype = landusefile.variables["PCT_NAT_PFT"].datatype
# variablelongname = landusefile.variables[pftname].long_name
# variablefillvalue = landusefile.variables["PCT_NAT_PFT"]._FillValue
variablevalues = np.zeros(shape=(numyears,ninlat,ninlon),dtype=float)
for pftidindex in range(numpftids):
pftid = pftids[pftidindex]
pftvalues = np.asfarray(landusefile.variables["PCT_NAT_PFT"][startyearindex:endyearindex+1,pftid,:,:],variabletype)
variablevalues += pftvalues
landusefile.close()

outresults = np.zeros(numyears,dtype=np.float64)
for yearindex in range(numyears):
inpctcropyear = inpctcrop[yearindex,:,:]
inpctnatveg = 100 - inpctglacier - inpctlake - inpctwetland - inpcturban - inpctcropyear
innatvegfrac = inpctnatveg / 100.0
variablevaluesyear = variablevalues[yearindex,:,:]
variablevaluesyeararea = inlandfrac * innatvegfrac * variablevaluesyear / 100.0 * inarea
outresults[yearindex] = variablevaluesyeararea.sum(axis=0).sum(axis=0)
outputfilename = str(sys.argv[5]) + ".csv"
outputstr = "{0:.3f}"

This produces the following plot in km^2 note scaled by 0.000001 to be millions

 

[lawrencepj1]

indexing did not carry over sorry

 

[James King]

Fantastic, thanks!

 

[James King]

Here we are - not quite the same values, but getting there!

 

[Hemraj]

Hi James

Sorry I just remembered you need to calculate the PCT_NATVEG as part of this calculation. As this is not on the file you need to take this as PCT_NATVEG = 100.0 - PCT_CROP (for the year from the landuse timeseries file) - PCT_GLACIER - PCT_LAKE - PCT_WETLAND (all on the surface data file).

So the area of a PFT should be = LANDFRAC_PFT*AREA* PCT_NATVEG[year] * PCT_NAT_PFT[year]/100

The LANDFRAC_PFT is simply the landfrac calculated for the land model which can be different to the land frac calculated for the ocean model (different resolutions and practicalities of modeling the ocean).

cheers
Peter

Hi Peter,
Do we also need to consider PCT_URBAN when we calculate PCT_NATVEG?

 

[James King]

Hi Peter,
Do we also need to consider PCT_URBAN when we calculate PCT_NATVEG?

Hi Hemraj,

I've included it in my calculations following Peter's script above.

 

[HazelMooney]

Hi James

Sorry I just remembered you need to calculate the PCT_NATVEG as part of this calculation. As this is not on the file you need to take this as PCT_NATVEG = 100.0 - PCT_CROP (for the year from the landuse timeseries file) - PCT_GLACIER - PCT_LAKE - PCT_WETLAND (all on the surface data file).

So the area of a PFT should be = LANDFRAC_PFT*AREA* PCT_NATVEG[year] * PCT_NAT_PFT[year]/100

The LANDFRAC_PFT is simply the landfrac calculated for the land model which can be different to the land frac calculated for the ocean model (different resolutions and practicalities of modeling the ocean).

cheers
Peter

Hi Peter,

I'd like to ask about how these calculations correspond to CLM4.5. I have been trying to unravel discrepancies in my area calculations for a little while. The discussion about for CLM5 above makes me wonder if there is a factor missing from my calculations for CLM4.5

In CLM4.5, my understanding is that we do not have PCT_CROP or PCT_NAT_PFT, we have PCT_PFT which has PFTs for crops within it. To calculate the area of a given PFT, i use: LANDFRAC_PFT*AREA*PCT_PFT/100

Am i missing a step?

Thanks

Hazel

 

[slevis]

I think that this makes sense, but I'm not looking at a clm45 fsurdat file right now. To confirm, I recommend looking at how things add up in a default clm45 fsurdat file, rather than in one that you have modified.

A caveat to point out... I'm pretty sure that LANDFRAC_PFT is not identical to the land area where the model runs. The latter should appear in your history output. Depending on what you want to calculate, LANDFRAC_PFT may or may not be the variable to use.