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.

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

Updates to ocean thermal forcing extrapolation to generalize ocean mask to 3d and diagnosing face melting errors with new mask #123

Open
wants to merge 5 commits into
base: develop
Choose a base branch
from

Conversation

cshafer
Copy link

@cshafer cshafer commented Aug 27, 2024

This PR makes changes to two scripts: (1) the ocean extrapolation script to generalize the validOceanMask and (2) the ice shelf melting script to diagnose issues related to grounded face melting.

  • The ocean extrapolation script was originally written to take in a 2d ocean mask origOceanMaskHoriz that is projected down through the ocean layers to create a 3d validOceanMask. We instead generalize the code to take in as input a 3d mask called orig3dOceanCavityMask that validOceanMask is set to. We also update a typo that makes sure that locations outside the validOceanMask (and not the availOceanMask) are set to the invalid_value_TF value.

  • In the ice shelf melting script, we modify which ocean level is used to populate TFocean. Previously, there was a step that interpolated between ocean layers if the bedrock topography existed between the two but was neither more shallow that the first ocean layer nor deeper than the bottom ocean layer. However, it is possible that the deeper ocean layer below the bedrock topography could contain invalid TF values, and so we modify it to take valid TF data only from above the bedrock topography as opposed to interpolating between the two layers.

  • If the bedrock topography is more shallow than the first ocean layer, then we keep kk = 1 since there is no 'good' ocean layer to switch to. In this case, TFocean could still get invalid TF values.

  • I also add a check to see whether TFocean is being populated with invalid TF values. This will need to be updated to check for invalid values where there should be valid extrapolated TF data in grounded locations as opposed to checking for invalid TF values anywhere.

@cshafer cshafer changed the title Updates to ocean thermal forcing extrapolation and diagnosing Updates to ocean thermal forcing extrapolation to generalize ocean mask to 3d and diagnosing face melting errors with new mask Aug 27, 2024

call mpas_log_write('==HH==: updating halos for the avail/valid ocean masks')
! Update halos
call mpas_timer_start("halo updates")
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@matthewhoffman removed these halo updates in 8b8e18b. Looks like these got entrained during the rebase.

…nd added check for invalid TF values"

This reverts commit 93c7ad7.
This commit makes adjustments to how the seafloor and inland seas
(locations below sea level not connected to the global ocean) are
handled for extraplation.  The adjustments are applied primarily to
availOceanMask, which is the mask of to where ocean data should be
extrapolated.

Changes made:
* Adjust availOceanMask to extend one layer below the seafloor (needed
  so that facemelting has a valid value one level below the seafloor for
  interpolation to the seafloor elevation)
   * to support this, generate error if ismip6shelfMelt_zBndsOcean is
     not populated, because that field is needed for the seafloor
     detection
* Adjust availOceanMask to ignore inland seas - we will not attempt to flood
  fill into areas below sea level that are not connected to the open
  ocean.
   * to support this, create mask of marine locations connected to global open ocean
* Adjustment to where invalid values are assigned to avoid inserting
  them in inland seas.  This will give inland seas either the value from the ocean
  data if it exists, or else TF=0
Eliminates 100s of repetitive log messages per timestep
Copy link

@alexolinhager alexolinhager left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

These changes look good to me! I made a few minor suggestions for edits to the code, but most of my comments are regarding the naming structure of the various masks. I know the naming convention was established in a different PR, but I feel like it's quite challenging to keep track of all of the different mask definitions/names, especially because a lot of them use very similar wording. I added a comment that may help fix this and make it more user friendly.

@@ -1689,6 +1691,10 @@ is the value of that variable from the *previous* time level!
<var name="origOceanMaskHoriz" type="integer" dimensions="nCells" units="none"
description="2D mask for original valid ocean data (e.g., ISMIP6 original TF data)"
/>
<!-- CAS 8/16/2024 -->
<var name="orig3dOceanCavityMask" type="integer" dimensions="nISMIP6OceanLayers nCells Time" units="none"
description="3D mask for original valid ocean data that includes cavities (e.g., SORRM original TF data)"

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I know this naming convention was introduced before this PR, but I still find it somewhat confusing; particularly, the difference between the origOceanMaskHoriz, orig3dOceanCavityMask, validOceanMaskOrig, and the validOceanMask. Can we change these names (and their descriptions) so that each has consistent terminology (e.g., remove the use of "original" or "orig" from descriptions/names of the "valid masks" and vice versa for the "original masks". Right now the use of "valid" and "original" in the descriptions/names of both types of masks makes it challenging to determine the differences between the two

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Farther down "availOceanMask" is described as a "grow mask", even though there is already a different variable called "growOceanMaskHoriz". Same with "validOceanMask" and "seedMask"

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

After looking at this a bit more, could we do the following?:

  1. Rename validOceanMask and validOceanMaskOrig to seedOceanMask, and seedOceanMaskOrig
  2. Remove the variables in (1) from the Registry and add them as allocatable variables, along with seedOceanMaskHorizInit, growOceanMaskHoriz, and seedOceanMaskHoriz (same as how seedOceanMaskHorizOld is treated).
  3. Rename orig3dOceanCavityMask to orig3dOceanMask

That way only the only two variables in the registry are orig3dOceanMask and 'availOceanMask', which can be easily defined and are the input/output variables we care about. I think that would streamline the code and make the distinctions between these variables less critical from a user perspective


! CAS 8/16/2024 Hijacking Holly's code to test using the 3D SORRM cavity TF field. We set 3d valid ocean mask to original 3d ocean cavity mask.
! We don't need to loop through the layers if we're using a 3d mask already
validOceanMask(:,:) = orig3dOceanCavityMask(:,:)

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It looks like you can now remove the validOceanMask(:,:) = 0 at line 187

! validOceanMask(iLayer,iCell) = 1
! endif
! enddo
!enddo

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can these commented lines be fully removed now?

@@ -915,7 +916,8 @@
<var name="ismip6_2dThermalForcingCurrent" packages="ismip6GroundedFaceMelt" />
<var name="forcingTimeStamp" packages="ismip6GroundedFaceMelt" />
<!-- these variables are just for the ocean thermal forcing re-extrapolation scheme -->
<var name="origOceanMaskHoriz" packages="extrapOceanData" />
<var name="origOceanMaskHoriz" packages="extrapOceanData" />

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This variable doesn't seem to be used anymore after the introduction of orig3dOceanCavityMask. Should be able to be safely removed from the Registry and the variable declarations in the main ocean extrap script.

@@ -89,8 +90,10 @@ subroutine li_ocean_extrap_solve(domain, err)
real (kind=RKIND) :: layerTop
real (kind=RKIND), dimension(:,:), pointer :: TFocean, TFoceanOld
real (kind=RKIND), dimension(:,:), pointer :: ismip6shelfMelt_3dThermalForcing, ismip6shelfMelt_zBndsOcean
real (kind=RKIND), dimension(:), pointer :: ismip6shelfMelt_zOcean
real (kind=RKIND), dimension(:), pointer :: thickness, bedTopography
integer, dimension(:), pointer :: origOceanMaskHoriz

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Variable no longer used

call mpas_pool_get_field(scratchPool, 'growMask', growMaskField)
call mpas_allocate_scratch_field(growMaskField, single_block_in = .true.)
growMask => growMaskField % array
growMask(:) = 0

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

growMask and seedMask have pretty similar in names to some of the other mask variables (see above comments). Could we rename these to avoid confusion, even if they are only used here?

! Calculate mask of connected ocean
call mpas_pool_get_subpool(domain % blocklist % structs, 'scratch', scratchPool)
call mpas_pool_get_field(scratchPool, 'seedMask', seedMaskField)
call mpas_allocate_scratch_field(seedMaskField, single_block_in = .true.)

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm not familiar with how these 3 lines work but I trust they are used correctly

realArgs=(/ismip6shelfMelt_zOcean(iLayer)/), intArgs=(/iLayer/))
err = ior(err, 1)
endif
if ((ismip6shelfMelt_zBndsOcean(1,iLayer) > 0.0_RKIND) .or. &

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Should this also be >=?

realArgs=(/ismip6shelfMelt_zOcean(iLayer)/), intArgs=(/iLayer/))
err = ior(err, 1)
endif
if ((ismip6shelfMelt_zBndsOcean(1,iLayer) > 0.0_RKIND) .or. &

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Could you add a comment describing what you are identifying with these next two if-statements? It looks like you are trying to ensure zBndsOcean is below sea level and that it exactly equals zOcean, but that seems to be different from the commit description of determining if zBndsOcean is populated or not

realArgs=(/ismip6shelfMelt_zBndsOcean(2,iLayer)/), intArgs=(/iLayer/))
err = ior(err, 1)
endif
enddo
availOceanMask(:,:) = 0
validOceanMask(:,:) = 0

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You should be able to remove this now

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

4 participants