Opened 8 months ago

Closed 6 months ago

Last modified 6 months ago

#3317 closed help (fixed)

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 (75)

comment:1 Changed 8 months ago by grenville

Markus

Is it OK to output GPP on domain DIAG?

Grenville

comment:2 Changed 8 months 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 8 months 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 8 months 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 8 months 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 8 months 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 8 months 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 8 months 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 8 months 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 8 months 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 8 months 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 8 months 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 8 months 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 8 months 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 8 months 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 8 months 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 7 months 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 7 months 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 7 months ago by pmcguire

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

comment:20 Changed 7 months 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 7 months 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 7 months 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 7 months 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 7 months 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 7 months 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 7 months 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 7 months 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 7 months 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 7 months 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

comment:30 Changed 7 months ago by mtodt

Hi Patrick

Thanks fo the suggestions! I hadn't noticed that growing temperature is included in the dump.

I've run for a couple more months, but the issues don't change. Growing temperature in physiol_jls is constantly 14.C even beyond the first month. And the range of growing temperature in sf_expl_jls is completely below -20C for October to December, so that is either wrong or only a subset of land points.

Looking at the reasonable growing temperatures in the dump files (mine included!), it must be me doing something wrong. Maybe at the point I'm including the print statements there's only a subset of land points used that are over Antarctica? That would match the temperature range. I still don't understand why I can't print out individual land points, though. Very frustrating.

Cheers
Markus

comment:31 Changed 7 months ago by pmcguire

Hi Markus
Can you think about trying http://cms.ncas.ac.uk/ticket/3317#comment:27 ?
I just spoke with Pier Luigi about this, and maybe that could at least partially fix the problem?
Patrick

comment:32 Changed 7 months ago by mtodt

Thinking about it, the range of t_growth_gb decreasing so quickly after starting the run probably indicates that the land points are somewhere over Antarctica and just decrease from the initialised values.

Cheers
Markus

comment:33 Changed 7 months ago by mtodt

Hi Patrick

Sorry, comment #27 was also one I hadn't seen before. I don't know why I missed the email notifications for a couple of your comments.

I'll try doing that. I'll probably come back with some questions since you have don the same before for sf_expl_jls, I assume?

Cheers
Markus

comment:34 Changed 7 months ago by mtodt

Hi Patrick

I've changed physiol_jls to read in t_growth_gb as you outlined above, and the ranges of t_growth_gb in physiol_jls and sf_expl_jls are the same now. So we have solved that issue, at least.

Cheers
Markus

PS: Notification emails have ended up in my junk folder over the last two days, for whatever reason. So that's why I didn't realise there were more responses.

comment:35 Changed 7 months ago by pmcguire

Hi Markus:
I am glad that that suggested fix to physio_jls improves things.
Did it affect the GPP on the different PFTs very much?

That's too bad that some of the CMS Helpdesk emails were/are going to your junk folder.
Patrick

comment:36 Changed 7 months ago by mtodt

Hi Patrick

No, unfortunately GPP still looks the same.

Cheers
Markus

comment:37 Changed 7 months ago by pmcguire

Hi Markus
Then maybe you can put some more print statements showing intermediate calculations of the GPP?
Patrick

comment:38 Changed 7 months ago by mtodt

Hi Patrick

Yes, that's what I wanted to do when I tried to print for specific grid cells. Maybe I'll simply go for the maxval and minval approach.

Cheers
Markus

comment:39 Changed 7 months ago by mtodt

Hi Patrick

To keep you updated, I have included various print statements checking gs_type, fsmc_pft, canht_pft, and lai_pft in physiol_jls as well as land_index and surft_index in sf_expl_jls.

Both land_index and surft_index are, at the point where I'm printing them, only 252 entries long, and the entries are 1:252 for land_index and surft_index(:,13). The 13th tile is land ice. For all other tiles surft_index is 0. So either the modules cover for whatever reason only what I assume is the southernmost part of Antarctica, or there's something going on with the printing. I've tracked the order of modules upwards (from where they are called) but I can't find any loop over tiles that would explain why only the last tile is covered by the print statements.

Regarding the PFT-specific parameters in physiol_jls, the printed values for canht_pft are realistic but lai_pft is consistently 10-10 and gs_type and fsmc_pft are 0. I assume that also indicates that only land ice grid cells are covered in physiol_jls, at least by the print statements.

As I mentioned in the opening text of this ticket, fluxes like canopy evaporation or stomatal conductance seem to be correct, so the unsatisfying print statements are probably my fault rather than an actual issue. I'm sorry I still haven't made much progress!

Cheers
Markus

comment:40 Changed 7 months ago by pmcguire

Hi Markus:
I am trying to find your log files where the print statements put their output. But I can't find them on PUMA in ~mtodt/cylc-run/u-bt231 . Is that the right place? Or did you put them somewhere else?
Patrick

comment:41 Changed 7 months ago by mtodt

Hi Patrick

Thanks for having a look! The print statements are written to bt231.fort6.pe000(0), and I save all of the files to /home/n02/n02/mtodt/GPP_Check_Logs/ on ARCHER.

Cheers
Markus

comment:42 Changed 7 months ago by pmcguire

Hi Markus
It looks like the sf_expl_jls is only printing out 14 x 18 grid cell information.
I think the usage is pdims might be wrong. I think tdims might be better to use for the array indexing and the calculations.
You should check this to be sure.
Patrick

comment:43 Changed 7 months ago by pmcguire

Hi Markus:
And maybe the land_pts should not be 14*18=252.
You can trace this back with print statements.
But I think the pdims is wrong, as I suggest in the previous comment.
Patrick

comment:44 Changed 7 months ago by mtodt

Hi Patrick

Thanks! Yeah, that came from the print statements I included just yesterday. But the limited amount of land points is not just an issue at that point where I added the print statement. So far it has been an issue everywhere I added print statements. So it's probably not due to the loop you added using pdims, don't you think?

Cheers
Markus

comment:45 Changed 7 months ago by pmcguire

Hi Markus
I will try to look at the other code.
But the land_pts is the limiter on the loop. This should be much larger.
Patrick

comment:46 Changed 7 months ago by mtodt

Yeah, that's a good point. I'll trace where land_pts stems from.

Cheers
Markus

comment:47 Changed 7 months ago by pmcguire

Hi Markus:
It is possible that it's just using the value of the land_pts for the land points that are being processed by one solitary processor for one rectangular region of the world.

Yes, tracing it back would be good.

Patrick

comment:48 Changed 7 months ago by mtodt

Hi Patrick

That would make sense. The print statements just don't get written out into log files. When I search in my suite directory, the bt231.fort6.pe000(0) file is the only one containing my print statements. Do I have to add anything else to the model code apart from WRITE(jules_message,*) … and CALL jules_print('…',jules_message)?

Cheers
Markus

comment:49 Changed 7 months ago by pmcguire

Hi Markus
If you look at your job.out file, you'll see that there are warnings with the string WARN_JULES_TEMP_FIXES in them.
So you can use this as a model for how to print to the job.out file.

If you do (on PUMA):

grep -r WARN_JULES_TEMP_FIXES ~mtodt/branches/jules/r16016_acclimation_um4/

you'll see that this WARN_JULES_TEMP_FIXES string is defined in:

~mtodt/branches/jules/r16016_acclimation_um4/src/control/shared/jules_science_fixes_mod.F90

If you edit that file:

vi ~mtodt/branches/jules/r16016_acclimation_um4/src/control/shared/jules_science_fixes_mod.F90

you'll see it uses the ereport() routine.

Maybe if you do something similar, you can write to the job.out file?
Patrick

comment:50 Changed 7 months ago by mtodt

Hi Patrick

Thanks, I'll try that! I had seen print statements like that but also a lot that are simply in the WRITE&CALL mould, which is why I tried that.

Cheers
Markus

comment:51 Changed 7 months ago by mtodt

I'm a bit puzzled. I've basically copied the ereport call from other modules and just included my text, but I get the following error during fcm_make_um:

Variable "ERRORSTATUS" must have its data scope explicitly declared because DEFAULT(NONE) was specified.

Have you encountered that one before? I wonder why that occurs now but not for the error reports that are already implemented.

Cheers
Markus

comment:52 Changed 7 months ago by pmcguire

Hi Markus:
This is just a quick guess. Maybe it's because that routine has two 'IMPLICIT NONE' 's?
Maybe both are necessary since one is for the module and one is for the subroutine.
Patrick

comment:53 Changed 7 months ago by mtodt

Hi Patrick

I had included the print statement within an $OMP PARALLEL DO loop, for which DEFAULT(NONE) is set. My run failed because I hadn't included parameter errorstatus in $OMP SHARED(...). The main atmos_main job is currently running.

Cheers
Markus

comment:54 Changed 7 months ago by pmcguire

Hi Markus:
Excellent! I am glad you figured it out.
Did you figure out what the deal is with the limited number of land_pts being printed?
Patrick

comment:55 Changed 7 months ago by mtodt

Hi Patrick

Another note: I initially used errorstatus=0 as this should result in a log_info instead of log_fatal or log_warn. But that didn't get written into any file, so I had to rerun with a number <0 to create a log_warn, which ends up in job.err.

Anyway, I checked the log file and it only writes out land points 1 to 252 as in the previous bt231.fort6.pe0(00) files. But the additional text makes it clear that only print statements for processor 0 are written out, nothing else. So I assume that it is as you said, the entire range of land points are indeed covered but spread out across processors, and for whatever reason only processor 0 creates print statements that end up in the log file.

On one hand I now know that there's probably nothing wrong with the code in that regard, but on the other hand I won't make any progress until I find out how to write print statements for all processors.

Cheers
Markus

comment:56 Changed 7 months ago by mtodt

While I can't figure out how to print variables and parameters over the entirety of land points, I've had a look again at the differences between trunks and acclimation branches for JULES and the UM, in case I've overlooked something or see something differently now. Some notes:

1) If I remember correctly you have said that you had run offline simulations of JULES with the same acclimation branch, and output from those runs looked alright?

2) All differences between the trunk version and acclimation branch of the UM are related to t_growth. But its field in dump files looks realistic, so that can't be the source of the issues with GPP, etc.

If both statements are true, the issue can only stem from (a) some part of the JULES branch where there's an if condition that checks whether JULES is used in standalone mode or as part of the UM, or (b) wrong parameter choices. Can you tell me the number of the suite you used for the offline runs with the acclimation branch so I can compare parameters?

Cheers
Markus

comment:57 Changed 7 months ago by pmcguire

Hi Markus
I didn't do any offline JULES sims with the acclimation code. I just did the coupled runs. I think Becky did the offline JULES runs with the acclimation suite.
Patrick

comment:58 Changed 7 months ago by pmcguire

Hi Markus:
Since it's only printing the output in sf_expl_jls.F90 for you for 252 of the land_pts, and this is only working maybe for one of the processors (which is all ice), maybe you can try using global_land_pts and gather_land_field( ) to print out non-ice GPP PFT information from other processors?

I am actually surprised that it isn't printing out the land_pts for all the processors.

You can also print out the mype processor number (or something like that), accessible from the um_parcore module, to see which processor is being used. Hopefully you can figure out how to get info printed out from all the processors.
Patrick

comment:59 Changed 7 months ago by mtodt

Hi Patrick

I've turned on printing for all processors as outlined by Jeff in this ticket. I don't get past the first 252 land points in any of the processors, though, when printing statements in physiol_jls. So either I'm still doing something wrong or am misunderstanding the code, or the land points actually aren't spread out evenly over the processors as I assumed initially.

Cheers
Markus

comment:60 Changed 7 months ago by pmcguire

Hi Markus:
You can print out land_pts in physio_jls, as well as land_pts in the routines that call physio_jls, hierarchically. I don't know if the land points are spread out evenly over all the processors or not. Maybe the grid cells are spread out evenly, but for a different rectangular subsection of the global grid, there are different numbers of land points in each sub section. You can also print out the processor number for each processor when printing out the t_growth_gb in physio_jls, if you aren't already. And you can also print out the lat/lon of each land point.

But if you're printing out say 32*252 land points (assuming 32 processors), then maybe you already have some land points that are not at the South Pole, and actually have vegetation? So maybe you an start to figure things out now with the GPP?
Patrick

comment:61 Changed 7 months ago by mtodt

Hi Patrick

Yes, I've been printing all through the JULES hierarchy surf_couple_explicit_modsf_expl_jlsphysiol_jlssf_stom_jls, but it's the same everywhere. I've gone into the higher UM modules now to check there. I've also tried printing lat/lon, but that didn't give me any reasonable values, which obviously means I've been doing something wrong. I mean, the model is running and simulating most fields correctly.

Sorry, I've probably meant nodes when talking about processors. After enabling print statements from all of them as outlined by Jeff, I get 576 pe log files. If each of them covered 252 land points that would equate to 145152 land points, which is slightly higher than the amount of land points in the model (139968). So that seems promising, which is why I'm confused that the pe log files repeat the range of land points that they cover.

Anyway, I've got an email from Becky about her suite, and I'm currently comparing her branch and parameter settings to ours.

Cheers
Markus

comment:62 Changed 7 months ago by pmcguire

Hi Markus
I looked on both PUMA and ARCHER for your pe log files. I guess you've already deleted them and started a new run, since I can't find them immediately. Maybe you can save the pe log files somewhere? Or maybe you can start a new copy of the suite, so that when you rerun the model, the old files don't get deleted?

So you're saying that the pe log files print out the same 252 values of t_growth_gb in physio_jls, for each of the 576 processors?
So it's still the ice grid cells for all 576 of the processors?

BTW, 576=32*18, which makes sense, since that's the domain decomposition of the grid.

Maybe for some reason, the coastal grid cells with land_frac<1 are counted differently, so that might be why 145152 != 139968.

Patrick

Last edited 7 months ago by pmcguire (previous) (diff)

comment:63 Changed 7 months ago by mtodt

Hi Patrick

I save them on ARCHER under /home/n02/n02/mtodt/GPP_Check_Logs/, but for runs with output from all nodes I can't save them all due to their file sizes. I have some of them in that directory under bt231.fort6.pe.v17/, they're from the first run with complete output from print statements, but only contain my checks of land_index, lon/lat, t_growth, and GPP for some grid cells I selected.

My notes regarding the last run, in which I tested the [MINVAL, MAXVAL] range of land_index, are:
However, the first dozens of processors all cover land_index 1 to 252, as do most of the processors. Other processors give the print statements "minimum land_index over land points 9223372036854775807 , maximum land_index over land points -9223372036854775808". I don't even know how that could work, those values have to be some kind of filler. Others, only a handful, feature ranges of land_index from 8 to 152, or other such ranges. But apart from those likely fillers, all ranges are within the first 252 land points.

Cheers
Markus

comment:64 Changed 7 months ago by pmcguire

Hi Markus
Do you have the copy of the FORTRAN file that you used for sf_expl_jls.F90 for v17 of the pe logs?
It seems you already made changes to this file, and I can't try to figure out why the pe logs don't make sense
with out looking at the old version of this file.
Patrick

comment:65 Changed 7 months ago by mtodt

Hi Patrick

I'd say sf_expl_jls.F90 was the same as it is now, just that all of my print statements were commented, and in physiol_jls lines 1627 to 1633 and lines 1643 to 1677 were uncommented (and 1679 to 1743 commented or non-existent at that point).

Cheers
Markus

comment:66 Changed 7 months ago by mtodt

Actually, I think I found out what happens to the spreading-out of land points over processors and nodes. I added a print statement to src/control/top_level/atmos_physics2.F90, which yields the first line in the output below. This is from file bt231.fort6.pe327 of the latest run:

dimensions: rows 18 , row length 14 , land points 32 , min land index 27 , max land index 252
...
sf_stom: minimum gpp over land points -1.27054942088145054E-24 , maximum gpp over land points 1.27054942088145054E-24
sf_stom: minimum fsmc over land points 0.15334544354126933 , maximum fsmc over land points 1.
sf_stom: minimum gpp over land points -2.54109884176290107E-24 , maximum gpp over land points 1.27054942088145054E-24
sf_stom: minimum fsmc over land points 0.1760571263858679 , maximum fsmc over land points 1.
sf_stom: minimum gpp over land points -7.62329652528870394E-24 , maximum gpp over land points 2.54109884176290107E-24
sf_stom: minimum fsmc over land points 0.1760571263858679 , maximum fsmc over land points 1.
sf_stom: minimum gpp over land points -2.54109884176290107E-24 , maximum gpp over land points 6.35274710440725268E-25
sf_stom: minimum fsmc over land points 0.18365087967218408 , maximum fsmc over land points 1.
sf_stom: minimum gpp over land points -1.47911419728939717E-34 , maximum gpp over land points 6.35274710440725268E-25
sf_stom: minimum fsmc over land points 0.1760571263858679 , maximum fsmc over land points 1.
sf_stom: minimum gpp over land points -1.27054942088145054E-24 , maximum gpp over land points 3.81164826264435197E-24
sf_stom: minimum fsmc over land points 0.34003772352106887 , maximum fsmc over land points 1.
sf_stom: minimum gpp over land points 1.64349925716335073E-17 , maximum gpp over land points 1.68008603859089821E-7
sf_stom: minimum fsmc over land points 0.34003772352106887 , maximum fsmc over land points 1.
sf_stom: minimum gpp over land points -2.54109884176290107E-24 , maximum gpp over land points 1.27054942088145054E-24
sf_stom: minimum fsmc over land points 0.24214763060194777 , maximum fsmc over land points 0.99999999999999989
sf_stom: minimum gpp over land points -1.27054942088145054E-24 , maximum gpp over land points 1.27054942088145054E-24
sf_stom: minimum fsmc over land points 0.24214763060194777 , maximum fsmc over land points 0.99999999999999989
...
surf_couple_explicit: min t_growth_gb over land points 21.493020441116638 , max t_growth_gb over land points 23.984000568413418

So land_index and land_pts don't continue across nodes, they always starts again at 1 and run to 252 (or less so that the total amount adds up to the amount of land points, which is less than 576x252). Node 327 seems to cover a tropical or sub-tropical area with grid cells on land, which I suppose makes sense since it's just a bit more than halfway through the amount of nodes, and values for fSMC and T_growth are realistic. And GPP values are as we find in the output files, i.e. multiple orders too low for all C3 PFTs and realistic for C4 grass (the 7th line beginning with sf_stom: minimum gpp over…).

So, I can roughly estimate which nodes to look at and printing steps before and after calculating GPP works. Here we go…

Cheers
Markus

Edit: I copied that file to /home/n02/n02/mtodt/GPP_Check_Logs/bt231.fort6.pe327.v18.

Last edited 7 months ago by mtodt (previous) (diff)

comment:67 Changed 7 months ago by mtodt

I think I found a clue. Scrolling through bt231.fort6.pe327.v18 you can see a diurnal cycle in GPP for C4 grass, and the nighttime values feature the same range of -x*10-24 to x*10-24 as all other PFTs do all the time.

comment:68 Changed 7 months ago by pmcguire

Hi Markus
Did you figure anything else out?
Patrick

comment:69 Changed 7 months ago by mtodt

Hi Patrick

Yes, I'm continuing to trace the issue back/upward through the calculations. It is indeed as I said on Friday, GPP has basically nighttime values for all C3 PFTs when Farquhar is on, but GPP looks realistic when I switch to Collatz.

Cheers
Markus

comment:70 Changed 7 months ago by pmcguire

Hi Markus
Since you say that the GPP has nighttime values for all C3 PFTs when Farquhar is on, then
since gpp(l) = cconu * (anetc(l) + rdc(l) * fsmc(l)) in sf_stom_jls_mod.F90, where the multiplier in the 2nd term is rdc (the dark respiration for the canopy, which occurs even at night), I would suspect that it's the first term that is causing problems (anetc, the net photosynthesis of the canopy).
Patrick

comment:71 Changed 7 months ago by mtodt

Hi Patrick

Yes, that is correct. Dark respiration has similar values for Farquhar and Collatz runs, but anetc does not. As I said before, I'm currently tracing the calculations upward/backward: anetcanetlanetl_sun & anetl_shd & fsunwcarb & wlitev_sun & wlitev_shdwlite & je_sun_ratio & je_shd_ratio.
I might send an email to Doug, now that I have some info on what's going on.

Cheers
Markus

comment:72 Changed 6 months ago by mtodt

  • Resolution set to fixed
  • Status changed from accepted to closed

Hi Patrick

As I mentioned in my email, I think I've fixed this issue. Parameters act_jmax_io, act_vcmax_io, alpha_elec, deact_jmax_io, deact_vcmax_io, ds_jmax_io, ds_vcmax_io, g1_stomata_io, jv25_ratio_io, and knl_io were initialised but not set to the values prescribed in namelist jules_pftparm. So they were used by the model but actually were all just 0. I've added these parameters to src/initialisation/um/surf_couple_allocate_mod.F90, and GPP now looks realistic. So I'm going to close this ticket now. Again, thanks a lot for your help over the last weeks!

Cheers
Markus

comment:73 Changed 6 months ago by pmcguire

Hi Markus:
Wow! I am glad you were able to fix this issue. It was probably some mistake I made that caused the problem. My apologies for that.
Patrick

comment:74 Changed 6 months ago by mtodt

Hi Patrick

No worries! I don't know at which stage those parameters should have been added, but it's certainly something that can slip through the cracks. As I said in my email, it wasn't noticeable in offline simulations.

Cheers
Markus

comment:75 Changed 6 months ago by pmcguire

I have started a 20-year run (1988-2018) on ARCHER based upon Markus's updated suite (fixing 2-3 bugs) u-bx260 with the Farquhar photosynthesis and acclimation switches turned on. This has fixed 1988 CO2 vales (343ppm). The suite number is u-bx393. The finalized UM and JULES branches (with Markus's changes) have been checked in to MOSRS.
Patrick

Note: See TracTickets for help on using tickets.