search

mom-ocean / MOM6

Modular Ocean Model
Other
185 stars 230 forks source link

Trouble with negative depths in land mask #1447

Closed kshedstrom closed 3 years ago

kshedstrom commented 3 years ago

In btstep, the surface height is compared to bathyT and a warning is issued if the water depths go negative. This can be triggered in the land where if one asks for a depth of -5 for tidal wetting and drying applications. SSH eta is 0.0 inside the mask which makes for a dry spot with that depth. It also causes trouble for the barotropic wave speed calculation based on depth.

In ROMS, when wetting and drying is turned on, SSH is always at least some D_crit above the bottom, even inside the land mask. The waves speeds use full water depth, including the contribution from SSH.

Hallberg-NOAA commented 3 years ago

In examining the code, it looks like btstep no longer assumes that is 0 a special geopotential height, and that the total depth is always the sum of bathyT+eta whereever it used in Boussinesq mode. Would it be sufficient to apply a land-mask before issuing a warning on lines 2358 and 2363 of MOM_barotropic.F90 to correct this issue in your case, or did you have something more extensive in mind, @kshedstrom ?

Also, there are several different places where external wave speeds are calculated in the barotropic solver. Could you be more specific about which one is becoming problematic?

kshedstrom commented 3 years ago

It is this one, this time:

#1  0x4df7a8 in __mom_open_boundary_MOD_update_obc_segment_data
    at //import/c1/AKWATERS/kate/ESMG/ESMG-configs/src/MOM6/src/core/MOM_open_boundary.F90:3820
#2  0x1ba4430 in __mom_boundary_update_MOD_update_obc_data
    at //import/c1/AKWATERS/kate/ESMG/ESMG-configs/src/MOM6/src/core/MOM_boundary_update.F90:136
kshedstrom commented 3 years ago

As for your other question, yes, applying a land mask to the checks would be sufficient.

Hallberg-NOAA commented 3 years ago

When I look at what is happening in the relevant code block where this is failing, https://github.com/NOAA-GFDL/MOM6/blob/dev/gfdl/src/core/MOM_open_boundary.F90#L3820-L3825, it strikes me that if we were to replace the existing code:

        segment%Cg(i,J) = sqrt(GV%g_prime(1)*G%bathyT(i,j+jshift))
        segment%Htot(i,J)=0.0
        do k=1,GV%ke
          segment%h(i,J,k) = h(i,j+jshift,k)
          segment%Htot(i,J)=segment%Htot(i,J)+segment%h(i,J,k)
        enddo

with 

    ! segment%Cg(i,J) = sqrt(GV%g_prime(1)*G%bathyT(i,j+jshift))
    segment%Htot(i,J)=0.0
    do k=1,GV%ke
      segment%h(i,J,k) = h(i,j+jshift,k)
      segment%Htot(i,J)=segment%Htot(i,J)+segment%h(i,J,k)
    enddo
    segment%Cg(i,J) = sqrt(GV%g_prime(1)*segment%Htot(i,J)*US%H_to_Z)

it would not only avoid having any problem with negative values of bathyT, it would also be doing the right thing physically. Does anyone foresee any problems with doing it this way?

I am assuming that we would need to wrap this change in a run-time-flag so that people could reproduce the existing answers.

kshedstrom commented 3 years ago

I approve this suggestion.

I for one am not wedded to existing suboptimal answers.

Hallberg-NOAA commented 3 years ago

Are there any remaining aspects to this issue that are unresolved now that PR #1470 (and several other related PRs that avoid the naïve use of G%bathyT) has been merged into dev/gfdl? If not, would you object to our closing this issue, @kshedstrom?

kshedstrom commented 3 years ago

Just a note to self - setting MASKING_DEPTH in SIS_input has to match that in MOM_input or you get trouble.