Possible bug in coupled sea surface slope/tilt field

[anton-seaice]

I compiled mom5 with :

mom5:
      require:
        - '@git.2025.05.000=access-esm1.6'
        - 'build_type="Debug"'

Running access-esm1.6 pre industrial configs give this error during initialisation:

forrtl: severe (408): fort: (33): Shape mismatch: The extent of dimension 2 of array SSLOPE is 21 and the corresponding extent of array <RHS expression> is 23

Image              PC                Routine            Line        Source
mom5_access_cm     0000000000917476  ocean_sbc_mod_mp_        2963  ocean_sbc.F90
mom5_access_cm     0000000000551336  ocean_model_mod_m        2091  ocean_model.F90
mom5_access_cm     0000000000440921  MAIN__                    450  ocean_solo.F90
mom5_access_cm     0000000000411C8D  Unknown               Unknown  Unknown
libc-2.28.so       000014FD54FEC7E5  __libc_start_main     Unknown  Unknown
mom5_access_cm     0000000000411BAE  Unknown               Unknown  Unknown

(There are lots of similar messages, implying it happens on every task).

The failure is in:

#if defined(ACCESS_CM) || defined(ACCESS_OM)
  sslope(:,:,:) = GRAD_BAROTROPIC_P(Thickness%sea_lev(:,:), isc - isd, isc - isd)
  if (isc - isd /= 1) then
    call mpp_error (FATAL, '==>Error from ocean_sbc_mod (sum_ocean_sfc): grad barotropic halos')
  endif 
#endif

sslope is of this size:

  allocate ( sslope(isc:iec, jsc:jec, 2) )

implying some sort of indexing error in the arrays

[dougiesquire]

I compiled mom5 with :

mom5:
      require:
        - '@git.2025.05.000=access-esm1.6'
        - 'build_type="Debug"'

Pandora's Box has been opened

[dougiesquire]

Some notes to self:

I think sslope should be allocated over the data domain, not the computational domain. This will fix the above issue and is consistent with other uses of GRAD_BAROTROPIC_P, e.g. to calculate grad_ps from ps in ocean_barotropic_mod.

Only the computational domain region will be added to Ocean_sfc%gradient in the subsequent loop:

      do j = jsc_bnd,jec_bnd
         do i = isc_bnd,iec_bnd
            ii = i + i_shift
            jj = j + j_shift
            ...
            Ocean_sfc%gradient(i,j,:) = Ocean_sfc%gradient(i,j,:) + sslope(ii,jj,:)
            ...

I'm not sure at this point whether this will change answers.

I do not understand the use of sslope in the avg_ocean_sfc routine. Specifically, why we use sslope(i,j,:) rather than sslope(ii,jj,:) here:

  if(.NOT. avg_sfc_velocity) then 
      if(horz_grid == MOM_BGRID) then 
          do j = jsc_bnd,jec_bnd
             do i = isc_bnd,iec_bnd
                ii = i + i_shift
                jj = j + j_shift
                if(Grd%tmask(ii,jj,1) == 1.0) then
                    ....
                    Ocean_sfc%gradient(i,j,:) = sslope(i,j,:)
                    ....

i_shift and j_shift are needed because arrays in ocean_public_type are allocated using absolute indices, regardless of whether ocean allocations are static. My guess is that i_shift/j_shift are zero in our models (I'll check this) but I think strictly we should be using sslope(ii,jj,:)...