Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Changing land masks ... #173

Open
AndyHoggANU opened this issue Oct 29, 2019 · 34 comments
Open

Changing land masks ... #173

AndyHoggANU opened this issue Oct 29, 2019 · 34 comments

Comments

@AndyHoggANU
Copy link
Contributor

Hi All

We are doing some large-scale bulldozing on ACCESS-OM2 ... just to test out potential changes to land ice on glacial changes. This is work with @Acosta-Goncalves

So far, after changing topog.nc and ocean_mask.nc in MOM and kmt.nc in CICE, we are finding divide by zero errors in CICE soon after reading in atmospheric variables.

So, the question is -- does anybody know if YATM/OASIS requires other mask changes, or if regrid weights will need to be computed? @nichannah - hoping you might be able to give us some pointers here.
@Acosta-Goncalves
Copy link

To further clarify, this is the error we're getting:

forrtl: error (73): floating divide by zero

Image              PC                Routine            Line        Source            
cice_auscom_360x3  0000000000946921  Unknown               Unknown  Unknown
cice_auscom_360x3  0000000000944A5B  Unknown               Unknown  Unknown
cice_auscom_360x3  00000000008F2874  Unknown               Unknown  Unknown
cice_auscom_360x3  00000000008F2686  Unknown               Unknown  Unknown
cice_auscom_360x3  0000000000870CB9  Unknown               Unknown  Unknown
cice_auscom_360x3  000000000087B669  Unknown               Unknown  Unknown
libpthread-2.12.s  00002B46BBEF07E0  Unknown               Unknown  Unknown
cice_auscom_360x3  0000000000455BD7  ice_atmo_mp_atmo_         352  ice_atmo.f90
cice_auscom_360x3  00000000006001A8  ice_step_mod_mp_s         395  ice_step_mod.f90
cice_auscom_360x3  000000000040FD4E  cice_runmod_mp_ci         348  CICE_RunMod.f90
cice_auscom_360x3  000000000040C887  MAIN__                     60  CICE.f90
cice_auscom_360x3  000000000040C81E  Unknown               Unknown  Unknown
libc-2.12.so       00002B46BC320D20  __libc_start_main     Unknown  Unknown
cice_auscom_360x3  000000000040C729  Unknown               Unknown  Unknown

The line in question (352) within ice_atmo.F90 is computing:

hol(ij) = vonkar * gravit * zlvl(i,j) &
                   * (tstar(ij)/thva(ij) &
                    + qstar(ij)/(c1/zvir+Qa(i,j))) &
                   / ustar(ij)**2

And the problem, specifically, seems to be in tstar(ij)/thva(ij)

@russfiedler
Copy link

Hi, as I mentioned on the slack channel this implies that the supplied air potential temperature, potT is zero K. It's initialised to 253K, 263K or 273K in init_coupler_flux so it must be getting the zeros from the coupler.

Note that the bug for which I approved the fix yesterday COSIMA/cice5#41 has the potential to cause this.

@Acosta-Goncalves
Copy link

Thanks @russfiedler. I did encounter the bug you mention in cpl_interface.F90 when I compiled the cice in debug mode, but fixing it did not solve the division-by-zero problem.

We've been looking through the code to find where potT could be getting zeros from, but could not find the problem. Only in ice_force.F90is the value of potT ever changed in the code, and it's done via: potT(i,j) = Tair(i,j). Since Tair is also initialized in the init_coupler_fluxalong potT to the same values (253K, ...), we can't find where the problem might be.

Any ideas more specific? Thanks again!

@Acosta-Goncalves
Copy link

I've made the change of only erasing Hawai'i from the map, just one tiny dot, and the problem is exactly the same. It has to be a masking problem, I thought maybe related to the hm mask? But if the kmt.ncfile is correct, I don't know why it's putting 0s where they don't belong...

@aekiss
Copy link
Contributor

aekiss commented Oct 29, 2019

Re. @AndyHoggANU's question, at 1 deg we use distribution_type=cartesian and processor_shape=slenderX1, which allocates pole-to-pole meridional strips to each processor (Craig et al 2015). Land-only blocks are not allocated a processor, so if the modified topog.nc has more ocean blocks there'd need to be some remapping changes. However that won't be an issue in this global model with processor_shape=slenderX1, since there are no land-only meridional strips.

@rmholmes
Copy link
Collaborator

I've run into the same error trying to do a perturbation experiment where I bulldoze the ITF...

@Acosta-Goncalves Acosta-Goncalves mentioned this issue Oct 30, 2019
Merged
@AndyHoggANU
Copy link
Contributor Author

@rmholmes - did you ever find a fix?

@rmholmes
Copy link
Collaborator

@AndyHoggANU nope, I just tried it yesterday briefly. happy to put some time in in the next few days if you can't solve it.

@aidanheerdegen
Copy link
Contributor

aidanheerdegen commented Oct 30, 2019 •
edited
Loading

We just found this old issue where Nic widened the Bering Strait. Alfonso will try using the tools Nic points to here:

#48

Specifically the mask changing tool unmask

https://github.com/COSIMA/topogtools

@russfiedler
Copy link

Do the remapping weights that yatm uses need to be recalculated when the land mask changes? I can't see anywhere in the steps outlined that tell you to do this. If zero weight is assigned because of the mask things will fall apart.

@nichannah
Copy link
Contributor

nichannah commented Oct 30, 2019 •
edited
Loading

Some thoughs:

  • The regrid weights should not need to be recomputed. They put the atmosphere onto the land as well as the ocean - the masking features of the weight generation software are just not used. This is a simplification but I don't think it is a performance problem because of land processor masking used by CICE. Also there is no regridding between ice and ocean so as long as ocean_mask.nc and kmt.nc match that should be fine.

  • I think @aidanheerdegen could be onto something with the unmask.py tool. When cice starts it reads OASIS restart files which are written out with 0's on the land (by OASIS). If there is new land this is OK but any new water points will be uninitialised in the restart. To check if this could be happening open o2i.nc and i2o.nc and see if there is land there. If so then it needs to be removed. unmask.py does this using nearest neighbour.

A longer term fix which hadn't occurred to me until now would be to modify OASIS to do what unmask.py does.

@rmholmes
Copy link
Collaborator

rmholmes commented Oct 30, 2019 •
edited
Loading

There doesn't seem to be any land points or weirdness in ice/RESTART/o2i.nc or ice/RESTART/i20.nc

If it's not in the remap files (in INPUT/ or atmosphere/INPUT), then the only other field I can see that could be a problem is the ocean/INPUT/ssw_atten_depth.nc, which has NaNs on land. However, I tried replacing these land values with something sensible and I still got the same error...

@russfiedler
Copy link

Isn't the attenuation calculated via the GFDL scheme? In that case it's the chl.nc that is used. I fail to see how any of the above relate to the air temperature being zero over the new wet points. This is coming from the atmosphere. Also the line of code where it is failing can only be reached if there is ice in the cell. Clearly, that's not the case for Hawaii or the ITF!

Also, to be perfectly clear, this is definitely happening in the first time step and not the second, correct? If the ocean is supplying a 0K temperature via ocean_sbc.res.nc it could be freezing the water in a strange way and creating some unexpected logic paths.

@Acosta-Goncalves
Copy link

Ok, so plotting the Tair forcing field that goes in to OASIS (tair_ai_matmxx_07.nc) gives values over the whole planet:
Tair IN

Plotting the Tair forcing field that comes out of OASIS (tair_i_cicexx_07.nc) shows all land having been masked (including Hawai'i):
Tair OUT

@aidanheerdegen
Copy link
Contributor

aidanheerdegen commented Oct 31, 2019 •
edited
Loading

So @nichannah this looks like the remap weights do have land masking. Is this because they are conservative?

Looks like the tair field is using the smooth remapping weights, so not conservative.

Is it possible to create remapping weights that don't use destination cell masking?

@aidanheerdegen
Copy link
Contributor

Ok, so I generated some new remapping files with no destination grid mask.

Made a change to the make_remap_weights.py 

$ git diff
diff --git a/tools/make_remap_weights.py b/tools/make_remap_weights.py
index 49dd54c..c6dab97 100755
--- a/tools/make_remap_weights.py
+++ b/tools/make_remap_weights.py
@@ -71,7 +71,7 @@ def convert_to_scrip_output(weights):
 
 def create_weights(src_grid, dest_grid, npes, method,
                    ignore_unmapped=False,
-                   unmasked_src=True, unmasked_dest=False):
+                   unmasked_src=True, unmasked_dest=True):
 
     my_dir = os.path.dirname(os.path.realpath(__file__))

and ran for the MOM1 case

./make_remap_weights.py /short/x77/nah599/access-om2/input/ /g/data1/ua8/JRA55-do/RYF/v1-3/ /short/x77/nah599/access-om2/input/yatm_1deg/ --atm JRA55 --ocean MOM1 --npes 4

Alfonso copied JRA55_MOM1_patch.nc to his input and changed the namcouple file to use this for tair

So. Question: is patch ok? Before it used rmp_jra55_cice_smooth.nc. I guess it is ok to replace all the fields that use the smooth remapping with the patch one we just made?

@aekiss
Copy link
Contributor

aekiss commented Oct 31, 2019

nice work @aidanheerdegen - I'm not sure exactly what smooth is (there's no mention of it here), but I thought we've been using patch for all non-conservative fields (or at least intending to) so that should be fine. It is supposed to be better than bilinear anyway (perhaps that what smooth refers to?)

@russfiedler
Copy link

The OASIS documentation says that you need to provide masks for SCRIPR and CONSERV remapping. If you don't it looks like this is done when you call oasis_enddef. This is done for each grid. It may be that the masking is done no matter what. I haven't delved into the code that far.

The short version is that the remapping weights should be recalculated whenever the land mask is changed. No ifs, no buts, no maybes, no shortcuts...

@PaulSpence
Copy link
Collaborator

PaulSpence commented Oct 31, 2019 via email

@russfiedler
Copy link

It's the same process as setting up the model grids in the first place. It's the changing of old restarts to a new land mask that is the tricky bit and @nichannah has created some nice tools for that.

@PaulSpence
Copy link
Collaborator

PaulSpence commented Oct 31, 2019 via email

@russfiedler
Copy link

Yes, This should be done and put on the wiki.

@nichannah
Copy link
Contributor

'smooth' is another name for patch. Apologies for the confusion there.

So @aidanheerdegen, your change has fixed the problem?

@aidanheerdegen
Copy link
Contributor

@nichannah Seems so yes. @Acosta-Goncalves has run a case with original weights and another with the new patch weights with no destination masking and after a year he can see some small (~4mm) differences in sea level around the equator, so sounds like it is fine. Not quite the same, as I guess there will likely be some small differences with how the interpolation works when it borders land.

We were unsure about making the same change to conserved fields as we weren't sure what would happen in that case.

I'm unsure about what @russfiedler refers to above, about OASIS and the CONSERV directives etc. Do those apply when we're supplying our own remapping weights files?

@aidanheerdegen aidanheerdegen mentioned this issue Oct 31, 2019
Merged
@russfiedler
Copy link

@aidanheerdegen Have a look at sections 4.3 and 5.1 of the OASIS documentation. Reading the documentation a bit more closely it seems that the weights are just read in even when SCRIPR/CONSERVE is specified. I honestly don't understand why SCRIPR and its options are even being specified in the namcouple file as they are overridden and just make things confusing.

Also, the ATM=>ICE fields should definitely be calculated without masking for conserved fields. You don't want to move radiation etc from land to the ocean/ice. There don't appear to be spurious values near the coasts from what I can tell so I guess the actual area is used when calculating the the weights but they get masked after.

@rmholmes
Copy link
Collaborator

rmholmes commented Dec 6, 2019

Hi everyone,

I am resurrecting this thread because I think there are still problems. I am doing an experiment where I widen the ITF by bulldozing some of the Indonesia. Following @Acosta-Goncalves I can get the run working by replacing rmp_jra55_cice_smooth.nc with @aidanheerdegen's JRA55_MOM1_patch.nc in namcouple. However, the short-wave and downward long-wave radiative fluxes are wrong in that region and give me a large surface cooling:
swflx
lw_heat

@Acosta-Goncalves sees the same thing in his remove-Hawaii run. I suspect the same is happening in his Antarctic runs but the radiative fluxes are small there so its not obvious.

So it looks like a new rmp_jra55_cice_conserve.nc is neccessary (this is used for the remapping of swfld and lwfld in namcouple)? There could be problems with rain and snow as well (I haven't output them)?

@aidanheerdegen @russfiedler can you help out here? Thanks!

@aidanheerdegen
Copy link
Contributor

Looks like you're right.

I'm trying to finish off a few things before I take the next week and a bit off, so won't have a chance to look at this.

@nichannah might be able to assist.

@Acosta-Goncalves
Copy link

@aidanheerdegen Didn't you have a script to remove the masking and generate the patch? Or how did you un-mask the land?

@aidanheerdegen
Copy link
Contributor

The previous fix was to a patch (smoothed) regridding, so I just turned off the option that took notice of masking in the destination cells.

@rmholmes said this is an error with a conservative regridding weights file. I don't know the best approach in this case. I guess just redo with altered mask for the destination, but this would require regenerating the weights every time the destination mask is altered.

If you could get away with ignoring the destination mask then that would be easier, but this is conservative regridding, so what happens in that case? Is that what we want to happen. This is not a 20 minute job on a friday afternoon.

@russfiedler
Copy link

Honestly, I'd just go and construct the regridding files, masks etc from scratch rather than go through this updating process. There's clearly something being missed out in the process and the time spent constructing the new remapping file is trivial compared to the running of the model.

@PaulSpence
Copy link
Collaborator

PaulSpence commented Dec 6, 2019 via email

@rmholmes
Copy link
Collaborator

rmholmes commented Dec 6, 2019 •
edited
Loading

I agree with @russfiedler - regenerating the weights everytime the mask changes is probably the way to go. I'll try this monday with make_regrid_weights.

@PaulSpence Alfonso has already documented the steps so far. I'll add to that and then PR the access-om2 wiki.

PS: no real rush with this. If there are problems I'll ping again in January.

@rmholmes
Copy link
Collaborator

rmholmes commented Dec 9, 2019 •
edited
Loading

I've got this working fine by making my own remapping files using my ocean_mask.nc and the make_remap_weights.py script. Using the new conserve and patch remapping weights I no longer have problems with the surface heat fluxes.

I've documented the process in draft form in the last section of the access-om2 wiki: https://github.com/COSIMA/access-om2/wiki/Tutorials
Most of this comes from @Acosta-Goncalves tutorial (thanks!). This should probably be checked by someone who knows more about what they are doing. E.g. @aidanheerdegen when he's back from holidays.

@aekiss
Copy link
Contributor

aekiss commented Oct 28, 2020

For the record, this repo shows how the 0.25deg topog.nc and associated ocean_vgrid.nc, ocean_mask.nc and kmt.nc were generated for /g/data/ik11/inputs/access-om2/input_20201022: https://github.com/COSIMA/make_025deg_topo

ocean_mask.nc and topog.nc were altered by ocean_mask_edits.txt and topog_edits.txt (respectively), which were generated by this gui tool: https://github.com/COSIMA/topogtools/blob/master/editTopo.py

This automated topography generation process makes it easy to alter the topography, but doesn't recalculate the remap weights.

Since the land mask was changed, the remap weights also needed regenerating by https://github.com/COSIMA/access-om2/blob/master/tools/make_remap_weights.sh using Russ' tripole bug fix to ESMF_RegridWeightGen https://github.com/COSIMA/esmf/tree/f536c3e (see #216).

Also see #158 (comment) and /g/data/ik11/inputs/access-om2/input_20201022/README.txt.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees
No one assigned
Labels
None yet
Projects
None yet
Milestone
No milestone
Development

No branches or pull requests
8 participants
@nichannah @aidanheerdegen @PaulSpence @russfiedler @rmholmes @AndyHoggANU @aekiss @Acosta-Goncalves