Modified seamount crashes with new version of code

[claireyung]

Hi, I'm just rerunning some experiments with idealised seamount and icemount shapes. I've noticed when I use the latest mom-ocean version (mom-ocean@d68ff1250e3fc8c), some of my runs crash when they used to run fine. The crash is

FATAL: NaN in input field of reproducing_sum(_3d).

Example experiment which is crashing for me: https://github.com/claireyung/IS-PG-MOM6/tree/main/seamount-z-MWIPG-Wright

This is happening to me on both my local mac M1 installation as well as on gadi (Australian supercomputer). Interestingly, my mac only seems to crash sometimes - sometimes it finishes running, and then sometimes it crashes, using the exact same config - I'm not sure if that means something is wrong with my mac installation.

I did a bit of exploration and I think I start getting this issue from this commit onwards: mom-ocean@8b63869

This kind of makes sense because the vertical viscosity algorithm is important to the seamount with vanished layers, and possibly I am using a combination of parameters not covered in other tests?

If it is helpful... I noticed that adding HARMONIC_VISC = True seemed to fix the issue but HARMONIC_BL_SCALE = 1 did not, where I would have expected the latter to be helpful in arresting vanished layer velocities.

[marshallward]

Some masking was changed to accommodate the loop index reordering. For example, kv_bbl and bbl_thick used to be set for all point, but are now masked with do_I. A NaN could have crept in and broken a reproducing sum, even if the answer was unchanged.

This loop (and the v- loop) originally was masked by do_i_shelf() but no longer.

MOM6/src/parameterizations/vertical/MOM_vert_friction.F90

Lines 1761 to 1763 in 83a1048

	do j=js,je ; do i=Isq,Ieq+1
	zcol(i,j) = zcol(i,j) - dz(i,j,k)
	enddo ; enddo

I don't think it could explain your error, but it's a discrepancy.

I did find one possible error, where h_arith and dz_arith are recomputed after the OBC correction was applied. But if you don't have OBCs, then it would just be a redundant calculation. (u-direction example below)

MOM6/src/parameterizations/vertical/MOM_vert_friction.F90

Lines 1781 to 1787 in 83a1048

	topfn = 1. / (1. + 0.09 * z2**6)
	h_arith(I,j) = 0.5 * (h(i+1,j,k) + h(i,j,k))
	dz_arith(I,j) = 0.5 * (dz(i+1,j,k) + dz(i,j,k))
	hvel_shelf(I,j,k) = min(hvel(I,j,k), (1. - topfn) * h_arith(I,j) + topfn * h_harm(I,j))
	dz_vel_shelf(I,j,k) = min(dz_vel(I,j,k), (1. - topfn) * dz_arith(I,j) + topfn * dz_harm(I,j,k))

Have you tried to run with DEBUG = True and compare the checksum output? It may pick up any NaNs before the reproducing_sum() call. I will try to do the same once I have a chance.

[marshallward]

I was able to reproduce this in GPU. It is a random error.

Backtrace:

#0  0x1494479 in __mpp_mod_MOD_mpp_error_basic
	at ../src/mpp/include/mpp_util_mpi.inc:82
#1  0x9406cf in __mom_error_infra_MOD_mom_err
	at /home/marshall.ward/gfdl/mom6/config_src/infra/FMS2/MOM_error_infra.F90:23
#2  0x93ff77 in loc_mom_err
	at /home/marshall.ward/gfdl/mom6/src/framework/MOM_error_handler.F90:191
#3  0x9403a9 in __mom_error_handler_MOD_mom_error
	at /home/marshall.ward/gfdl/mom6/src/framework/MOM_error_handler.F90:178
#4  0x6bb6a3 in __mom_coms_MOD_reproducing_sum_3d
	at /home/marshall.ward/gfdl/mom6/src/framework/MOM_coms.F90:460
#5  0xe39402 in __mom_sum_output_MOD_write_energy
	at /home/marshall.ward/gfdl/mom6/src/diagnostics/MOM_sum_output.F90:749
#6  0x4beec2 in __mom_MOD_step_mom
	at /home/marshall.ward/gfdl/mom6/src/core/MOM.F90:1150
#7  0x8738d0 in mom6
	at /home/marshall.ward/gfdl/mom6/config_src/drivers/solo_driver/MOM_driver.F90:488
#8  0x874e27 in main
	at /home/marshall.ward/gfdl/mom6/config_src/drivers/solo_driver/MOM_driver.F90:27

This is the summation containing the NaN:

 742   ! Calculate the Kinetic Energy integrated over each layer.                    
 743   tmp1(:,:,:) = 0.0                                                             
 744   do k=1,nz ; do j=js,je ; do i=is,ie                                           
 745     tmp1(i,j,k) = (0.25 * GV%H_to_RZ*(areaTm(i,j) * h(i,j,k))) * &              
 746             (((u(I-1,j,k)**2) + (u(I,j,k)**2)) + ((v(i,J-1,k)**2) + (v(i,J,k)**2)))
 747   enddo ; enddo ; enddo                                                         
 748                                                                                 
 749   KE_tot = reproducing_sum(tmp1, isr, ier, jsr, jer, sums=KE, unscale=RZL4_T2_to_J) 

With DEBUG = True this is the first NaN detected:

FATAL: NaN detected: u BT vel update    1 PF[uv]

Checksum diff

    435 u-point: mean=   1.7604644817101703E+02 min=   0.0000000000000000E+00 max=   2.0000000000000000|     435 u-point: mean=   1.7604644817101703E+02 min=   0.0000000000000000E+00 max=   2.00000000000000
    436 u-point: c=     13562 u vertvisc_coef h_[uv]                                                   |     436 u-point: c=     13562 u vertvisc_coef h_[uv]                                    
    437 v-point: mean=   8.8023224088521076E+01 min=   0.0000000000000000E+00 max=   2.0000000000000000|     437 v-point: mean=   9.1137730745482969E+01 min=   0.0000000000000000E+00 max=   2.00000000000000
    438 v-point: c=      7498 v vertvisc_coef h_[uv]                                                   |     438 v-point: c=      6552 v vertvisc_coef h_[uv]                                                 
    439 u-point: mean=   2.0721668635258004E+08 min=   0.0000000000000000E+00 max=   5.9600073405754566|     439 u-point: mean=   2.0721668635258004E+08 min=   0.0000000000000000E+00 max=   5.96000734057545
    440 u-point: c=     52530 u vertvisc_coef a_[uv]                                                   |     440 u-point: c=     52530 u vertvisc_coef a_[uv]                                    
    441 v-point: mean=   2.2855771971750526E+25 min=  -4.4313450191090950E+15 max=   3.8215423911621091|     441 v-point: mean=   9.5842225132413104E+07 min=   0.0000000000000000E+00 max=   6.16089675829023
    442 v-point: c=     26958 v vertvisc_coef a_[uv]                                                   |     442 v-point: c=     27005 v vertvisc_coef a_[uv]                                                 

The u-direction seems OK, the error seems to be in the v-direction.

[marshallward]

There were two serious errors:

diff --git a/src/parameterizations/vertical/MOM_vert_friction.F90 b/src/parameterizations/vertical/MOM_vert_friction.F90
index c5a085298..69f6f927b 100644
--- a/src/parameterizations/vertical/MOM_vert_friction.F90
+++ b/src/parameterizations/vertical/MOM_vert_friction.F90
@@ -1543,7 +1543,7 @@ subroutine vertvisc_coef(u, v, h, dz, forces, visc, tv, dt, G, GV, US, CS, OBC,
     endif
 
     if (OBC%u_W_OBCs_on_PE) then
-      do j=js_W_OBC,je_W_OBC ; do I=Is_W_OBC,Is_W_OBC
+      do j=js_W_OBC,je_W_OBC ; do I=Is_W_OBC,Ie_W_OBC
         if (do_i(I,j) .and. OBC%segnum_u(I,j) < 0) then
           Dmin(I,j) = G%bathyT(i+1,j)
           zi_dir(I,j) = 1
@@ -2145,7 +2145,7 @@ subroutine vertvisc_coef(u, v, h, dz, forces, visc, tv, dt, G, GV, US, CS, OBC,
             zcol(i,j) = zcol(i,j) - dz(i,j,k)
           enddo ; enddo
 
-          do J=Jsq,Jeq ; do i=is,je ; if (do_i_shelf(i,J)) then
+          do J=Jsq,Jeq ; do i=is,ie ; if (do_i_shelf(i,J)) then
             zh(i,J) = zh(i,J) + dz_harm(i,J,k)
 
             hvel_shelf(i,J,k) = hvel(i,J,k)

The i=is,je would be the one breaking your run.

That seems to fix everything for me, can you test it?

[claireyung]

Hi Marshall, thanks so much for looking at this so promptly! I changed the i=is,je you found, and I can confirm that fixed everything for me too.

Thank you!!
Claire

[marshallward]

Thanks for testing it out, and for catching the mistake. I have to do a second rewrite of those functions, so I will run it past you before submitting it to dev/gfdl and main.

[claireyung]

Just adding now that I realised my ISOMIP z coordinate quiet configs are also crashing, I think since that same commit. The indexing fix didn't seem to resolve my issue. The config is https://github.com/claireyung/IS-PG-MOM6/tree/main/z-ISOMIP-quiet. Wondering if you can reproduce that error too or if it's a problem on my end? Sorry that I had not checked these last time, I had checked all the seamount and icemount cases but not the 3D ISOMIP.

Thanks

[marshallward]

Copying post in #973:

I see a blowup in the new version after 1 hour:

WARNING: btstep: eta has dropped below bathyT:  -2.9718959833467682E+02 vs.  -2.9718414221305966E+02 at   4.8500E+02  1.5000E+01     83      8

I see a checksum diff, but it is appearing in the barotropic code, first in BT_force_[uv] and then in DCor_[uv].

It may be worth doing a git bisect to see if this was caused by a different commit.

[claireyung]

Thanks Marshall, might be a different commit - I will explore in more detail.

Thanks

[marshallward]

The previous diffs may have been a result of changes to OBC masking, and not relevant. But I also notice diffs in the barotropic wt_[uv] and gtot_[NESW] (which depends on wt_[uv]`).

I set VISC_REM_BUG = True and VISC_REM_TIMESTEP_BUG, which affect the calculation of wt_[uv], and did not see the eta btstep error.

There was a change to default settings recently, this may have caused issues, masked by the other ice shelf bug which you reported.

This might be time for @Hallberg-NOAA to comment, I don't know much about these flags or the bug that was fixed.

[Hallberg-NOAA]

Although we did change the default settings, all of these should be properly logged, and any parameters that changed in your configurations as a result can be isolated by comparing the MOM_parameter_doc files with the old and new versions of the code.

The issue you are describing, @claireyung, seems to me to be reminiscent of some problems that we were having with other ice-shelf test cases, and which we have subsequently fixed, in part with #962. This commit adds the new runtime parameter RESCALE_STRONG_DRAG which can be set to True, which, along with setting BT_STRONG_DRAG = True, did wonders to clean up large velocities that were occurring in this regions of water beneath ice shelves.

I would also encourage you to take a look at #972, which fixes the problems with the iterations to determine the balanced 3-equation heat and fresh water fluxes under ice sheets. This commit will change answers when ICE_SHELF_BUOYANCY_FLUX_ITT_BUGFIX = true but by properly bracketing the solution we are looking for and choosing appropriately between Newton's method and the false position method it should give much more reliable answers.

[claireyung]

Thanks for pointing me to those new parameters Bob! Confirming with those two strong drag parameters set to True, my quiet zstar ice shelf test case now runs with the latest dev/gfdl code. Still not sure why it used to run and then didn't, but at least I have working configurations now. I will still try to dig further and see if any parameters had changed from before.

And thanks for fixing the 3-eqn parameterisation iteration!!

[OlgaSergienko]

Hi Claire @claireyung, a reason for your crushes is that Alistair has changed
WARNING: btstep: eta has dropped below bathyT: ...
to
FATAL: btstep: eta has dropped below bathyT: ...

Sorry for not letting you know about Bob's BT_STRONG_DRAGfix. Glad to hear that you've got your configuration working now.