Opened 6 years ago

Closed 6 years ago

Last modified 4 years ago

#1389 closed Bug (fixed)

icbutl.f90 : Numerical instability caused by the interpolation of ocean properties in ICB module.

Reported by: nacho Owned by: gm
Priority: low Milestone:
Component: OCE Version: release-3.5
Severity: Keywords: ICB
Cc:

Description

  • Setting : ORCA025.L75 with mutlproc and icb


  • The dynamics of icebergs is solved with a 4 steps Runge-Kutta scheme. In the middle of those Runge-Kutta steps, an iceberg can move to a grid point out of the processor.


  • ocean and sea ice properties are interpolated in the new iceberg position in icb_utl_interp() subroutine which uses icb_utl_bilin() function.


  • but icb_utl_bilin() does not handle extra allocation of arrays and can therefore produce out of bounds accesses to arrays (tested with Check Boundary compiler option)


  • This may lead for example to unphysical values for SST, the iceberg melts produceing a big amount of freshwater input, which eventually leads to a numerical instability in the ocean dynamics.


## Suggested bug fix


  • The  proposed solution is to use  icb_utl_bilin_h() instead of icb_utl_bilin() for interpolating properties phi, psst and pcn in icb_utl_interp() function. (indeed icb_utl_bilin_h() works with expanded arrays).


  • those arrays should be defined as expanded arrays in icb_oce.F90 and in icb_utl_copy() as well.


  • I also added a piece of code inside icb_utl_bilin_h(). It forces the iceberg position to be at least one point away from the extended allocation size of the array before computing the interpolation.


  • This is probably a very unlikely situation, but the ICB code could be more robust thanks to that.


  • Here is  the code:
if (ii.lt.mig(1)) then
        ii = 1
      else if (ii.gt.mig(jpi)) then
        ii = jpi
      else
        ii  = mi1( ii  )
      end if

      if (ij.lt.mjg(1)) then
        ij = 1
      else if (ij.gt.mjg(jpj)) then
        ij = jpj
      else
        ij  = mj1( ij  )
      end if

      if (ij.eq.jpj) ij=ij-1
      if (ii.eq.jpi) ii=ii-1  
  • After those modifications, the model keeps stable for more than 15 years of simulation with ICB module.

Commit History (1)

ChangesetAuthorTimeChangeLog
4957acc2014-12-02T15:21:47+01:00

Branch dev_MERGE_2014. Style changes and implementation of ICB bugfix #1389.

Change History (5)

comment:1 Changed 6 years ago by gm

  • Owner changed from NEMO team to gm

comment:2 Changed 6 years ago by acc

  • Resolution set to fixed
  • Status changed from new to closed

Implemented on dev_MERGE_2014 branch @ revision 4957

comment:3 Changed 4 years ago by nicolasmartin

  • Keywords icebergs added; iceberg removed

comment:4 Changed 4 years ago by nicolasmartin

  • Keywords model removed

comment:5 Changed 4 years ago by nemo

  • Keywords ICB added; icebergs removed
Note: See TracTickets for help on using tickets.