Opened 5 weeks ago

Last modified 2 hours ago

#3317 accepted help

UM selectively simulates tiles

Reported by: mtodt Owned by: pmcguire
Component: UM Model Keywords: UM output tiles
Cc: pmcguire Platform: ARCHER
UM Version: 11.4

Description

Hi

I've noticed in the output of my UM suite (u-bt231) that Gross Primary Productivity (GPP) on tiles is only realistic for C4 grass (on the order of 10-8 and always positive). GPP for all other PFTs is somewhere between -1*10-24 and 1*10-24. I can see that the regions of those PFTs are correct, i.e. the grid cells that feature non-zero values, it's just that the values are strange and unrealistic.
My suite features 13 tiles in total, 9 of which are PFTs. I've checked output from the original 5-PFT suite I based mine on, and it has the same issue. GPP is only properly simulated for C4 grass.
When I look at other output fields from my simulation that are PFT-specific, I see differences in magnitudes between C4 grass and other PFTs for Potential Plant Respiration and Potential Net Primary Productivity, but values for all PFTs seem similar for canopy evaporation or stomatal conductance.

Now I don't understand what this means and where this issue comes from. Does it mean there's something wrong in how GPP is simulated/calculated? Is there something wrong with my reconfiguration/start dump? Or is this related to how GPP is written and it might actually be simulated correctly?

Patrick, you know more about the code that involves GPP, etc. than me. Does that ring a bell for you?

Thanks a lot for your help in advance!

Cheers
Markus

Change History (29)

comment:1 Changed 4 weeks ago by grenville

Markus

Is it OK to output GPP on domain DIAG?

Grenville

comment:2 Changed 4 weeks ago by mtodt

Hi Grenville

GPP on DIAG features realistic values but only for the grid cells where there's considerable coverage by C4 grasses.

Both my 9-PFT suite and the original 5-PFT suite use the same UM and JULES branches, which we have created, so we're currently checking whether this issue stems from our branches.

Cheers
Markus

comment:3 Changed 4 weeks ago by pmcguire

Hi Markus:
I looked at the output of the 5-PFT acclimation suite (u-bq532q) for the GPP of the different PFTs.
Yes, it only seems right for the C4 grass. Not for the other PFTs.
Did you check the output for non-acclimation (trunk) runs?
Patrick

comment:4 Changed 4 weeks ago by pmcguire

Hi Markus:
I just looked at the output for the 5-PFT non-acclimation suite (u-br310). The GPP for the different PFTs looks OK there.
That data is in /gws/nopw/j04/porcelain_rdg/pmcguire/archer_archive/u-br310 .
The main difference is that photo_model=1(=photo_collatz) for this one instead of photo_model=2(=photo_farquhar).
But you probably have figured this out already.
Patrick

comment:5 Changed 4 weeks ago by pmcguire

Hi Markus:
in this code:

/home/pmcguire/jules/r16016_acclimation_um4/src/science/surface/sf_stom_jls_mod.F90

There is this code snippet:

!-----------------------------------------------------------------------------
! Decide which photosynthesis model will be used for this PFT.
!-----------------------------------------------------------------------------
SELECT CASE ( photo_model )
CASE ( photo_collatz )
  ! C3 and C4 both use a Collatz model.
  pft_photo_model = photo_collatz
CASE ( photo_farquhar )
  ! C3 uses Farquhar, C4 uses Collatz.
  IF ( c3(ft) == 1 ) THEN
    pft_photo_model = photo_farquhar
  ELSE
    pft_photo_model = photo_collatz
  END IF
END SELECT

So since when you are using photo_model = 2 (= photo_farquhar) and it only works
for C4 grasses, this would make sense, all the other PFTs are considered to be C3 and use
pft_photo_model = photo_farquhar. There must be something wrong with pft_photo_model = photo_farquhar.
Patrick

comment:6 Changed 4 weeks ago by mtodt

Hi Patrick

Thanks a lot for having a look at this as well! That's a bit of a bummer. I was hoping the issue was due to one particular parameter or variable that was calculated differently for C3 and C4 PFTs, which would have been easier to find. I suppose the error could be anywhere within the Farquhar part.

Cheers
Markus

comment:7 Changed 4 weeks ago by pmcguire

Hi Markus

I can't figure out right now why the stomatal conductance for the photo_model=photo_collatz=1 version (u-br310) is so high in the summer (July) for trees in the far north of Canada in the summer, where there are no trees and the soil is bare. There must be something wrong.

Patrick

comment:8 Changed 4 weeks ago by pmcguire

Hi Markus

Do you have output from any runs that don't use the acclimation branches of JULES+UM at all?
Can I look at those?

Patrick

comment:9 Changed 4 weeks ago by mtodt

Hi Patrick

No, I do not unfortunately. I just checked but it looks like I didn't get any of my previous test suites to the point that they'd produce output.

Cheers
Markus

comment:10 Changed 4 weeks ago by pmcguire

Hi Markus
Yes. I just did the same. I haven't found any output from the previous test suites I did either. The only output I have found is for these two production suites.
Patrick

comment:11 Changed 4 weeks ago by pmcguire

Hi Markus
I just found that I ran u-bq486 with UM11.0 trunk (no accclimation( for Sept1988, producing 1 month of output.

The output is on ARCHER at:
/work/n02/n02/pmcguire/cylc-run/u-bq486/share/data/History_Data
But the standard suite doesn't have GPP or GPP on PFTs as output, as far as I can see.
Patrick

comment:12 Changed 4 weeks ago by pmcguire

Hi Markus
I see that we're using can_rad_mod=6. I am not sure right now if that is necessary or not. Maybe the Farquhar code doesn't need or require that. And maybe that's why GPP /approx 0.
Patrick

comment:13 Changed 4 weeks ago by pmcguire

Hi Markus:
You might want to look at this part of the code on PUMA:

~pmcguire/jules/r16016_acclimation_um4/src/science/surface/sf_stom_jls_mod.F90

That's where the photosynthesis parameters are defined.
Patrick

comment:14 Changed 4 weeks ago by pmcguire

Markus:
I just found an email that says that the Farquhar model needs to run with can_rad_mod=6.
Patrick

comment:15 Changed 4 weeks ago by pmcguire

Hi Markus:
I just looked at the air temperatures at 1.5 meters height above the ground, for both the Farquhar run (u-bq532) and the Collatz run (u-br310). It looks like for many areas for example, for July 1993, the temperatures in the Farquhar run are 10-15 degrees warmer than in the Collatz run. This might explain the GPP being near 0 in the Farquhar acclimation run for C3 PFTs — the reference temperature is 298K, and if it's 10-15K above that, then maybe the exponential jmax and vcmax temperature factors maybe are exponentially big or small?
Patrick

comment:16 Changed 4 weeks ago by mtodt

Hi Patrick

Sorry for the delay, I've worked on something else today, and thanks for having a look!

The temperature differences are indeed weird, but are they global? There aren't any spatial differences in GPP for the C3 PFTs, they are somewhere between -1*10-24 and 1*10-24 around the globe. Maybe the temperature differences are a consequence of the issues in productivity and respiration?
Anyway, I'll have a look at the code you pointed out.

Cheers
Markus

comment:17 Changed 4 days ago by mtodt

Hi Patrick

I've been trying to add print statements to the acclimation code that give me GPP and growing temperature (and if those are correct other variables and parameters further backward through the calculations) for a couple of locations that I've selected based on the GPP maps. However, I just don't get any useful information. I've calculated what the land_index of those locations should be based on their position in the N216 grids so that I can print them in the JULES routines where calculations are done over land points:

land_index(l) = (j-1)*t_i_length +i

That information I got from

DO l = 1,land_pts
   j = (land_index(l) - 1) / t_i_length + 1
   i = land_index(l) - (j-1) * t_i_length
END DO

where t_i_length is the row length, i is the position within the row, and j is the number of a row.
I've assumed that i refers to longitude and j refers to latitude, but I also checked where these land points would be if i referred to longitude and j referred to latitude. Anyway, what's printed doesn't work out for any of these versions.

I've then counted which land point the locations should refer to, in case my calculation of the land_index is wrong, and printed land_index, longitude, latitude, growing temperature, and GPP for those land points. What I get is:

physiol: land point  22718 , land index  0 , longitude  -3.47306054607043318E-164 , latitude  -3.47306054607043318E-164 , t_growth_gb  0. , GPP  215.3751950264386
physiol: land point  22606 , land index  0 , longitude  -3.47306054607043318E-164 , latitude  -3.47306054607043318E-164 , t_growth_gb  0. , GPP  217.72262550930856
physiol: land point  22660 , land index  0 , longitude  -3.47306054607043318E-164 , latitude  -3.47306054607043318E-164 , t_growth_gb  0. , GPP  216.12691851036638
physiol: land point  20753 , land index  0 , longitude  -3.47306054607043318E-164 , latitude  -3.47306054607043318E-164 , t_growth_gb  0. , GPP  228.74340499029998
physiol: land point  22247 , land index  0 , longitude  -3.47306054607043318E-164 , latitude  -3.47306054607043318E-164 , t_growth_gb  0. , GPP  219.87913983181713
physiol: land point  38129 , land index  0 , longitude  -3.47306054607043318E-164 , latitude  -3.47306054607043318E-164 , t_growth_gb  0. , GPP  1.
physiol: land point  36781 , land index  0 , longitude  -3.47306054607043318E-164 , latitude  -3.47306054607043318E-164 , t_growth_gb  0. , GPP  1.

That's clearly wrong. These specific land points are within the range of land points, which I've checked, so the locations should correspond to some land_index, longitude, and latitude. I just don't know what to make of this. Do you have any idea? My code is under /home/mtodt/branches/jules/r16016_acclimation_um4/, and the print statements are included in src/science/surface/physiol_jls_mod.F90 (as well as src/science/surface/sf_stom_jls_mod.F90).
Thanks a lot for your help!

Cheers
Markus

comment:18 Changed 3 days ago by pmcguire

Hi Markus
Are you still working on this? I see that you commented out the print statements in your code.
The sample output that you have from these print statements is doesn't make sense.
What I would suggest is to print out the entire arrays for one time step, and see if the output then makes sense.
The lats and lons should increment monotonically, for example, and the GPP and and the t_growth_gb should have reasonable values
for all grid boxes.
Patrick

comment:19 Changed 3 days ago by pmcguire

  • Owner changed from um_support to pmcguire
  • Status changed from new to accepted

comment:20 Changed 2 days ago by mtodt

Hi Patrick

Thanks for helping with this! I've tried to print out parts of an array before but got the following error:

lib-4212 : UNRECOVERABLE library error
An internal WRITE tried to write beyond the end of an internal file.
Encountered during a list-directed WRITE to an internal file (character variable)

I assumed that was due to me printing too many values, but when printing only for one time step that could work. How do I do that, though? How do I check for a particular time step?

Cheers
Markus

comment:21 Changed 2 days ago by pmcguire

Just do the run for one time step, and then stop the run.
Or if you want to continue, then check with an if statement whether or not it is the first time step or not.
Patrick

comment:22 Changed 25 hours ago by mtodt

Hi Patrick

What I had done in the meantime was checking minval and maximal of t_growth_gb rather than printing the entire array. I added print statements in physiol_jls and in sf_expl_jls right before physiol_jls gets called. Here are two examples from the first month after starting the run (on 1st September 1960):

Model time:   1960-09-01 00:45:00
t_growth_gb in sf_expl_jls before calling physiol: minimum t_growth_gb over land points 14.901793115654556 , maximum t_growth_gb over land points 14.923597512479047
physiol: minimum t_growth_gb over land points 14. , maximum t_growth_gb over land points 14.

and

Model time:   1960-09-30 22:00:00
t_growth_gb in sf_expl_jls before calling physiol: minimum t_growth_gb over land points -39.661518139185283 , maximum t_growth_gb over land points -29.24153435544163
physiol: minimum t_growth_gb over land points 14. , maximum t_growth_gb over land points 14.

There are several strange things.
1) The range in t_growth_gb over land points in sf_expl_jls is always very small, which suggests it doesn't cover global land points. The minimum and maximum values in sf_expl_jls decrease monotonically without any diurnal variations, though.
2) The range in t_growth_gb over land points in sf_expl_jls decreases to values at the end of September that can only be for (sub-)polar regions. Or there's something wrong with the calculations. (Or my print statement is wrongly placed, of course.)
3) Obviously, something in physiol_jls doesn't quite work. It seems like t_growth_gb is constant spatially and temporally. I've noticed that while in sf_expl_jls it is imported through its call in a higher module, in physiol_jls it is imported by USE fluxes, ONLY: t_growth_gb. Could that be an issue?

Cheers
Markus

comment:23 Changed 24 hours ago by mtodt

Hi Patrick

I've also tried printing the entire GPP array and aborting after the first time step. But atmos_main fails during the first time step as I previously wrote with:

lib-4212 : UNRECOVERABLE library error 
  An internal WRITE tried to write beyond the end of an internal file.
Encountered during a list-directed WRITE to an internal file (character variable)

Cheers
Markus

comment:24 Changed 21 hours ago by pmcguire

Hi Markus:
Maybe for writing out the entire array, you can write out instead every 10th element? What syntax are you using to write out the array?
Maybe you can test the writing of the entire array in a super-small test program?
Patrick

comment:25 Changed 21 hours ago by pmcguire

Hi Markus:
You can see in src/initialisation/shared/allocate_jules_arrays.F90 that we initialise t_growth_gb as a constant of 14.0C on all land points. It will take some time to do the temporal moving average for each grid box to get to the actual value. The filtering time is be default n_day_photo_acclim=30 days. So you might need to wait 60-90 days to get reasonable values.

It would be good to associate lats/lons with the temperatures you're reporting.

Furthermore, if the value right before physiol_jls is 14C everywhere after 1 month, then probably there is something wrong with passing the variables, as you suggest, or that it's getting reset to the initial value of 14C after each time step, or something.
Patrick

comment:26 Changed 20 hours ago by mtodt

Hi Patrick

Thanks for the clarification! I'll run a longer simulation then.

Yes, I've tried printing lon and lat for specific land points, but that didn't work as I've written. I haven't solved that yet.

Regarding the t_growth_gb discrepancy between physiol_jls and sf_expl_jls, might that be due to reading t_growth_gb from fluxes? So maybe it works as follows: since there's no temporally averaged t_growth_gb yet the values from fluxes are still the initial 14C while in sf_expl_jls values are still calculated?
Nonetheless, that the values decrease to ←20C so quickly indicates that something is wrong here, don't you think?

Cheers
Markus

comment:27 Changed 20 hours ago by pmcguire

Hi Markus:
Maybe you can try to get rid of the lines in physiol_jls that say:

USE fluxes, ONLY: t_growth_gb

and replace that with passing t_growth_gb directly as an argument to the SUBROUTINE physiol routine?
You'll probably have to change several files to do this.
Patrick

comment:28 Changed 19 hours ago by pmcguire

Hi Markus:
For writing out arrays with jules_print, maybe you can just write out a smaller array, and it won't have an error message?
Maybe just the first ten elements: array(1:10) ?
Or maybe printing out with a cadence of 10 or 100: array(1:1000:100) or array(1:100:10) . I am not sure if that syntax will work, so maybe you can check that.

That way you can also print out the same thing with lats/lons and maybe photosynthesis parameters, for multiple array elements, more easily.
Patrick

comment:29 Changed 2 hours ago by pmcguire

Hi Markus:
The t_growth_gb array is already written out to the dump file at the beginning of each month, as stash code 206. Maybe you can use xconv to inspect that? You can do this for example on JASMIN in /gws/nopw/j04/porcelain_rdg/pmcguire/archer_archive/u-bq532q.
The geographic maps of this time-averaged temperature (30-day moving average) are not totally unreasonable, at first glance.
Patrick

Note: See TracTickets for help on using tickets.