soilsnow developments from AM3

[JhanSrbinovsky]

CABLE

Description

Developments from AM3 work. Crazy number centre around using an ice density in place of liquid density, or 0.9. * liquid density

iThis is why I left this til the end as there will be differences resulting from these mods

Fixes #494

Type of change

• model improvement

Checklist

• The new content is accessible and located in the appropriate section.
• I have checked that links are valid and point to the intended content.
• I have checked my code/text and corrected any misspellings

Please add a reviewer when ready for review.

---

📚 Documentation preview 📚: https://cable--497.org.readthedocs.build/en/497/

[JhanSrbinovsky]

benchcab WILL reveal changes to output. In the process of downloading/uploading. My 3 site local test on gadi is pretty reasonable. Will post link to ME.org when done

[JhanSrbinovsky]

still sorting out the ME.org plots. My testing shows seasonal diurnal quantities at Tumba (as expected) negligible. One of the sites I test at is Hyttialia (because there is snow): Note seasonal labels appropriate for Sothern Hemisphere

[JhanSrbinovsky]

I think it looks as expected, I'll do a global one and see how it goes (I have to find the code first!)

[JhanSrbinovsky]

This is using GSWP2 forcing. Most of the plots that flashed in front of my eyes as the scripts kept turning over were unremarkable I thought. But the surface SoilTemp caught my eye as it is one that it make sense to do a relative difference as it never gets near zero etc. As with all the others it is pretty bland, except for this bright spot in Northern Russia. On closer inspection this bullseye represents a ~10% diff. But it goes away after a few months. Odd, I wonder what could be causing this? @har917 any ideas? Greenland is slightly cooler in Summer but not nearly as shocking. The Russian 10% difference appears to be an initialisation issue as it doesnt come back towards the end of the year. But still what is so dramatically different here. There was nothing here in AM3 was there?

[har917]

@JhanSrbinovsky could be a feature of the initialisation - could be something to do with the ice fixes.

I'd be cautious about looking at %change in soil temperature - there's little context of a 10% difference in temperature - depending on the units (C vs K) you can identify different regions for concern. e.g. C will be particularly susceptible to divide by near-zero regions, using K will results in even large (+3K) differences only showing up at the 1% level.

[JhanSrbinovsky]

@har917, its only really the ice fixes that are in here. but why should they apply to this bullseye in particular. I just assumed it was Kelvin as almost everything in the UM is K. If it is in C, Im not so surprised. I just assumed this 1 spot was experiencing a ~30degree cooling, so probably C. I'll check

[JhanSrbinovsky]

damn it. its in Kelvin

[JhanSrbinovsky]

@har917 This lazy screenshot shows that almost everything is within -1% to +1% diff. That might be 3K and I can't comment on whether or not that is acceptable (or even expected) BUT the bullseye in Russia is 10 times that. Is this something to worry about, or should we move on?

[har917]

unfortunately I think we need to understand the origin of a ~30K difference in soil temperature - the most likely cause is that something that is now being used is getting through uninitialised (or to a default value of 0).

[JhanSrbinovsky]

I dont have any ideas about what to look at. I might just look at the usual suspects and see what the moisture is doing at the same point. what else?
because it is isolated, and seems to go away by the next year I am hoping it is just a dodgy initialisation. Keep in mind as well that this is a no feedback, forced simulation.

[ccarouge]

Since Ian says there is some digging to do to understand the differences, I've removed the request for review. This way I'm not coming back in here constantly to see if it's ready. @JhanSrbinovsky just ask for reviews again when ready.

[JhanSrbinovsky]

Trouble spot zoomed in (rel. Diff ) T(K)



SoilMoist - Diff


SoilMoistIce - Diff



Iveg



Rather than show you the dozens of field which are unremarkablekable:
snowdepth

Curiously the snow water equivalent does not show this. A clue that I can’t make sense of. I can imagine an initiallized snow depth which is all of a sudden subjected to a different density, but this wouldn't be the only place

[JhanSrbinovsky]

This is interesting, the albedo:

[JhanSrbinovsky]

are we fortuitously hitting a limit?

[har917]

Jhan - can you check snow density? There is a link between snow density and albedo (and a limit on that impact) in surface_albedosn()

This is a region with a lot of (frozen) lakes - remember we changed how met%tk/ssnow%tggsn got used - so that could be a factor as well.

[JhanSrbinovsky]

sounds promising - I'll have to re-run both models - it isnt in my output. I imagine snow density output will be available through an easy switch

[JhanSrbinovsky]

bugger - it isnt!! it is available in a restart but this isnt much good to us. I'll come up with something a bit later.

[JhanSrbinovsky]

ok im back now. Does anyone know if it is possible to change the dump frequency in the offline model? This would be the easiest way I can think of to get a single point out of a global offline run?

[JhanSrbinovsky]

Hey @SeanBryan51, you've recently been looking at the _drivers for offline. Is writing the restart coded after the timestep loop?

[SeanBryan51]

Hey @SeanBryan51, you've recently been looking at the _drivers for offline. Is writing the restart coded after the timestep loop?

Apologies, I haven't looked into the restart IO in the drivers so I'm not too familiar on where the logic is in the code. This has not been touched yet so the restart behavior should be unchanged

[har917]

So canopy%fe should be zero right? So maybe fe is not initialised to zero and remains zero unless a finite canopy response overwrites it?

canopy%fe is the total latent heat flux - so includes contributions from both vegetation %fev (wet and dry leaf components) and soil (%fes) - all of these are initialised in cable_parameters. If I understand the code tracing %fes gets set in define_canopy before it is used (in update_zetar) and definitely should have a value by the time it gets out of cbm.

Basically I don't understand how the code is failing at that point in mpi_master - except is there something weird going on with the MPI code and one or other of the master/worker side exchanges isn't set up properly?

now - seeing your comment - have we really got %tgg=0.0? I understand all/most of the temperatures are in K so it should be in the range 260-300 at initialisation and should be coming from the restart file.

however - since we're enforcing consistency between %isoil and %iveg we may be triggering an entirely new tile which wouldn't/may not have an initialisation in the restart file.

[JhanSrbinovsky]

Sorry I meant NaN (tgg, potev)

[har917]

okay %tgg and %potev being NaNs indicate that something upwind isn't quite correct - and this is consistent with %fe and %fes being NaNs.

Can you check what veg%hc, rough%z0soilsn and rough%z0 are doing? (The main suspect is %rtsoil - and possibly %tss or %qstss - this depends on these 3)

I'm wondering whether we've managed to get it (with the change in surface types) into a state where one part of the code thinks a tile is vegetated and other parts think it's not. I can't see how consistency_ice_veg_soil could get us into that state but worth checking.

[JhanSrbinovsky]

@har917 - I haven't looked at height and roughness etc yet.what would hc have to do with this on ice? Anyway, I suspect tss etc will be screwed up. I was looking at what kills the debug run in cable_parameters. The value of ssat is zero at many points, and at our troublesome point. ssat is always finite in the nml file, but I suspect it is getting this from grid_info. ASSUMING that all the points where it messes up ssat are glacial, clobbering ssat with the ice value from our nml, it crashes at a later point in cable_parameters. Actually I'm finding a bunch of parameters that are zero. I wager that if I used the nml values it would at least get passed this build=debug issue. Another thing I could try is if this is still a problem with the old grid_info. I did consider that how is this working in the UM? I thought perhaps myassumption of these being ice points could be right. Although, I dont know how JULES manages with ice points then. Anyway, quick test seemed to suggest this was indeed the case.

[JhanSrbinovsky]

@har917 I clobbered ssat for ice=9 and it indeed got past the initial crash. The next was one also involving a soil parameter. I decided to just use clobber all the ice points as we do in the UM, which would make sense as thats where @Whyborn 's grid info is coming from. Playing in cable_parameters to see where I should do this I cam across this little gem . Someone has been here before with this problem. As per my inline comment some of the values line up with those in our nml file. Some have a units conversion which I can believe. But some are just whacky. Like for exampl e.g. the density =1455. How does ice end up as dense as clay!?!

[JhanSrbinovsky]

BTW - crash in debug goes away with these clobbering :)

[har917]

what would hc have to do with this on ice?

I am/was wondering whether in enforcing consistency between PFT and soil we've managed to get some weird combination of parameter values at the land point concerned. hc is important as it's critical in determining whether CABLE-roughness thinks we have a vegetation canopy or not - and that is a control on rtsoil, then potev (which has gone to NaN).

However - you've identified that ssat = 0.0 at that point. This will cause NaNs all over the place (there's lots of xx/ssat through the soil_snow code).

Now if I look at consistency_ice_veg_soil again in the sections that it checks whether a land point is soil_ice (so line 2500 and line 2523) the code is only overwriting soil%rhosoil and soil%css - but it's not overwriting the other soil parameters (ssat etc.). I think this is wrong/incomplete - @ccarouge your thoughts?

This perhaps explains why a value of ssat = 0.0 has gotten through the consistency check but I don't understand how it got a value of 0.0 in the first place - since it should have picked up a realistic value when it was set as non-ice soil originally.

[JhanSrbinovsky]

This perhaps explains why a value of ssat = 0.0 has gotten through the consistency check but I don't understand how it got a value of 0.0 in the first place - since it should have picked up a realistic value when it was set as non-ice soil originally.

I think all the zeroes are coming from the grid info file which is from the UM. We clobber ice points, but we dont pass them back so presumably this might be all those zeroes. JULES must be doing something similar and have ice point specific values, otherwise I dont know how they mange to avoid all the potential sources of (fatal( error

[JhanSrbinovsky]

whoaaaaa - it doesnt help us at all with our original problem - it spreads it out - this is day one only - where the original magnitude of the bullseye on day 1 was ~66K

[JhanSrbinovsky]

also keep in mind that this is the ESM1.5 grid info

[JhanSrbinovsky]

woo hoo 👯 using branch main at the branch point of 494 VS 494 and this @Whyborn grid info for both:

This is DAY one! I'll make up some more plots but I am hopeful that even averages over a month will calm this down. I dont think this exceeds what might be expected weather wise from the ice fixes? @har917 I've also got it restarting from a restart using the old grid info so I'll try to clean t his up as well.

[JhanSrbinovsky]

@rml599gh - it has become clear that ice points are misrepresented using the grid info file coming out of ACCESS. Possibly you have encountered this and have made adjustments already? - but this would certainly impact online/offline comparison

[har917]

@JhanSrbinovsky This is looking a lot more reasonable - the differences are mainly in regions where we expect melt and/or rain onto snow to occur - and that was one part of the ice fixes changeset.

I'd expect different (larger?) differences to emerge over time since that will involve soil moisture dynamics with different amounts of water around.

@rml599gh @ccarouge Given that the offline/online comparisons have focused on performance over vegetated tiles I'm not so sure that this is an issue for that work. But we do need the gridinfo that's aligned to ACCESS-ESM to have appropriate values for the soil properties given that we may want to do LSMIP-style runs and use offline-as-ACCESS for spin up purposes in ESM3.

@JhanSrbinovsky I am a little lost with what additional code has had to go in here - in particular around the values for the soil parameters that have been hard-wired. These should be the same as the Greenland values that exist in the default gridinfo file.

[JhanSrbinovsky]

@JhanSrbinovsky I am a little lost with what additional code has had to go in here - in particular around the values for the soil parameters that have been hard-wired. These should be the same as the Greenland values that exist in the default gridinfo file.

Thats what I was expecting. Possibly @rml599gh has some memory of them.

Im getting lost with latest plots showing bizarre things. Have to go pick up Lena now

[rml599gh]

When I was writing my ACCESS-based gridinfo for my offline testing, I did have a special case in the ice grid-cells. I assume I did this based on what was in the UM code for ESM1.5 which, from memory, just assumed one particular soil type under all ice grid-cells and then used the relevant soil parameter values (from or the same as in def_soil_params.txt.
Here's the relevant code from my 'writegridinfo.F90'

           if (iveg(i,j).eq.17) then
                   isoil(i,j) = 9
                   rhosoil(i,j) = 910.
                   css(i,j) = 2100.
                   bch(i,j) = 7.1
                   clay(i,j) = 0.3
                   css(i,j) = 2100.
                   hyds(i,j) = 0.000001
                   sand(i,j) = 0.37
                   sfc(i,j) = 0.301
                   silt(i,j) = 0.33
                   ssat(i,j) = 0.479
                   sucs(i,j) = -0.153
                   swilt(i,j) = 0.216
           else
                   isoil(i,j) = 2
                   rhosoil(i,j) = 1600.
                   css(i,j) = 850.
           endif

[JhanSrbinovsky]

When I was writing my ACCESS-based gridinfo for my offline testing, I did have a special case in the ice grid-cells. I assume I did this based on what was in the UM code for ESM1.5 which, from memory, just assumed one particular soil type under all ice grid-cells and then used the relevant soil parameter values (from or the same as in def_soil_params.txt. Here's the relevant code from my 'writegridinfo.F90'

These are the same values from our nml file but not the ones that someone has added (BLAME doesn't recognise this suggesting it was in the code 3 years ago). In either case they were not necessarily the values in cable_parameters where it was crashing.

I've been wondering where the code came from. It is possible that I did it during HAC/offline comparisson. I dont specifically remember it but there was parameter clobbering at some stages.

[JhanSrbinovsky]

Cant believe that yet again I have to re-write this- grrrr:

The crank values I was seeing was due to a versioning problem, which I am more than happy to share if you’re interested.

I have created a new branch from the main branch as of last night. 494_redux_local. This isolates the AM3 ice fixes. Necessary compensation to use the ESM1.5 gridinfo file is the altered cable_parameters.F90. For your perusal this is committed to 494_redux_local, and present in my local copy of main. Comparison isolates impact of ice fixes.

I thought I could change branch in the PR, Im sure I ve done this before, anyway goto #527 to see Files changed.

[ccarouge]

@JhanSrbinovsky I don't follow your last comment. Are you saying this current PR is now useless, it has been superseeded by #527 and this PR (#497) should be closed? If so, can you close it?

To change the branch at the top click the edit button next to the title and then click the branch name.

[JhanSrbinovsky]

@JhanSrbinovsky I don't follow your last comment. Are you saying this current PR is now useless, it has been superseeded by #527 and this PR (#497) should be closed? If so, can you close it?

To change the branch at the top click the edit button next to the title and then click the branch name.

That was instinctively how I tried to do it but it didnt seem to work. I'll try again and close #527 if successful

[JhanSrbinovsky]

Nope - every time after edit if I click on the branch it t takes me to the branch. The option to change the branch to compare to is there

[JhanSrbinovsky]

The branch is a mess because:

The feature branch 494-merge-soilsnow-.... branched from main ~6 weeks ago. When I went to rebase, it rebased local main which was still at 6 weeks ago. It doesnt want to rebase again because it thinks it has already done it. It was easier to rebranch from main and merge 494-....... which is what rebase does anyway.

[Whyborn]

I've made a new gridinfo with the ice fixes applied as per Rachel's comment. Comparison between soil saturation (ssat) in the old gridinfo vs the new (note the difference in scales):

Previous gridinfo ssat:

Gridinfo ssat with ice fix:

The new gridinfo file is at /g/data/rp23/experiments/2024-03-12_CABLE4-dev/lw5085/CABLE-as-ACCESS/gridinfo/ACCESS-ESM1p5-1p875x1p25-gridinfo-CABLE-withicefix.nc. Something happened to all the NetCDF metadata which I'm working to fix, but I don't think CABLE looks at the gridinfo metadata so it should still work?

Edit: Now has all the metadata attached.

[har917]

@JhanSrbinovsky It's probably worth testing on GSWP2/3 with the @Whyborn's updated grid but without the (re)implemented lines cable_parameters and the removed error catch in cable_mpimaster (i.e. revert the codebase to just include the AM3 soilsnow ice fixes code).

Hopefully getting the parameter layers in the gridinfo correct will solve both problems.

We could look to refine later (especially to remove the error catch) - but it would be good idea to focus the PR.

[JhanSrbinovsky]

sorry just got distracted by a super important phone call - apparently I've won $US1000 in crypto. I just have to l.og into some random website and collect :) - back from lunch - Yeah im confident this will work out of the box and the ice fixes can go in independent of any other code change

[JhanSrbinovsky]

to be clear - when I say ice fixes Im referring to the code changes in soilsnow/ , not the same ice fixes @Whyborn is referring to

[JhanSrbinovsky]

yep. It's the equivalent of what I had this morning.

[JhanSrbinovsky]

moving this to #527