Changing land masks ...

AndyHoggANU commented 5 years ago

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 commented 5 years ago

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 _iceatmo.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 commented 5 years ago

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 https://github.com/COSIMA/cice5/pull/41 has the potential to cause this.

Acosta-Goncalves commented 5 years ago

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 commented 5 years ago

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 commented 5 years ago

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 commented 5 years ago

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

AndyHoggANU commented 5 years ago

@rmholmes - did you ever find a fix?

rmholmes commented 5 years ago

@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 commented 5 years ago

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

https://github.com/COSIMA/access-om2/issues/48

Specifically the mask changing tool unmask

https://github.com/COSIMA/topogtools

russfiedler commented 5 years ago

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 commented 5 years ago

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.

There's also a tool to fix ocean restarts after a landmask change. It has some documentation at the beginning of the file. https://github.com/COSIMA/topogtools/blob/master/fix_restarts.py

rmholmes commented 5 years ago

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 commented 5 years ago

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 commented 5 years ago

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

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

aidanheerdegen commented 5 years ago

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 commented 5 years ago

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 commented 5 years ago

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 commented 5 years ago

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 commented 5 years ago

I hope someone is going to make a nice clear document of how to change the land sea mask in these models ....

Thanks Paul

On Thu, Oct 31, 2019 at 5:43 PM russfiedler notifications@github.com wrote:

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...

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/COSIMA/access-om2/issues/173?email_source=notifications&email_token=ABSWJXCYL4XHI5WLCMXJSF3QRJ5ABA5CNFSM4JGCHVFKYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOECWW5FQ#issuecomment-548236950, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABSWJXC3PR7ISI4O2PFNKMLQRJ5ABANCNFSM4JGCHVFA .

-- Paul Spence, PhD Senior Lecturer, Climate Change Research Centre University of New South Wales Sydney, Australia paul.spence@unsw.edu.au http://goog_934148345 www.iampaulspence.com

russfiedler commented 5 years ago

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 commented 5 years ago

Thanks for the quick response Russ. Is there existing documentation on "how to set up the model grids in the first place" for these models? I am thinking of a comparison to the existing MOM-SIS-FMS documentation for creating the model grids.

Thanks Paul

On Thu, Oct 31, 2019 at 5:54 PM russfiedler notifications@github.com wrote:

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 https://github.com/nichannah has created some nice tools for that.

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/COSIMA/access-om2/issues/173?email_source=notifications&email_token=ABSWJXDD3UVUQWFQFT7WIV3QRJ6JFA5CNFSM4JGCHVFKYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOECWXPZA#issuecomment-548239332, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABSWJXA4PVS3ZX3E5BCNTF3QRJ6JFANCNFSM4JGCHVFA .

-- Paul Spence, PhD Senior Lecturer, Climate Change Research Centre University of New South Wales Sydney, Australia paul.spence@unsw.edu.au http://goog_934148345 www.iampaulspence.com

russfiedler commented 5 years ago

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

nichannah commented 5 years ago

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

So @aidanheerdegen, your change has fixed the problem?

aidanheerdegen commented 5 years ago

@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?

russfiedler commented 5 years ago

@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 commented 4 years ago

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 commented 4 years ago

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 commented 4 years ago

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

aidanheerdegen commented 4 years ago

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 commented 4 years ago

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 commented 4 years ago

When this is sorted, I am begging you (who???) to document the correct/complete process.

Thank you, P

Paul Spence, PhD Senior Lecturer Climate Change Research Centre UNSW, Sydney, Aus. www.iampaulspence.com

On Fri, 6 Dec 2019, 4:30 pm russfiedler, notifications@github.com wrote:

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.

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/COSIMA/access-om2/issues/173?email_source=notifications&email_token=ABSWJXBUDSB5Z3UOZ3PMJB3QXHPORA5CNFSM4JGCHVFKYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOEGDCCJY#issuecomment-562438439, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABSWJXFKEDG35QINWUUVB43QXHPORANCNFSM4JGCHVFA .

rmholmes commented 4 years ago

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 commented 4 years ago

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 commented 4 years ago

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 https://github.com/COSIMA/access-om2/issues/216).

Also see https://github.com/COSIMA/access-om2/issues/158#issuecomment-714267777 and /g/data/ik11/inputs/access-om2/input_20201022/README.txt.

COSIMA / access-om2

Changing land masks ... #173