Version 18 (modified by mathiot, 7 years ago) (diff)

Last edited Timestamp?

Author : Mathiot

ticket : #1331

Branch : dev_r4650_UKMO2_ice_shelves


Dynamics for the top layer beneath the ice shelves:

  • top partial step
  • top friction (same option as for the bottom friction)
  • Correction of the divergence locally (as for runoff fresh water flux)

Thermodynamics used to compute the melt rate:

  • ISOMIP formulation (2 equation formulation)
  • 3 equation formulation (Jenkins et al. 1991)
  • Forcing mode (prescribed melt rate) is available.
  • Losh top boundary layer (Losch et al., 2008)

Validation configuration: ISOMIP

ISOMIP is a closed rectangular basin of uniform depth 900 m, spanning 15° of longitude, and latitudes 80° South to 70° South. The whole basin is covered with an ice shelf, with the ice draft rising linearly from 700 m at 80° South to 200 m at 76° South, and remaining constant at 200 m from 76° South to 70° South. The geometry is uniform in the east-west direction. The vertical coordinates is z . All the model setup and geometry is fully described in J.R. Hunter, 2006.


  • Hunter, J.R.: Specification for test models of ice shelf cavities, technical report, Antarctic Climate & Ecosystems Cooperative Research centre, version 7, 2006.
  • Jenkins, A., Hellmer, H.H. and D.H. Holland: The Role of Meltwater Advection in the Formulation of Conservative Boundary Conditions at an Ice–Ocean Interface, Journal of physical oceanography, 31, 2001.
  • Losch, M. Modeling ice shelf cavities in a z coordinate ocean general circulation model Journal of Geophysical Research: Oceans, 2008, 113


  • System reviewer: Andrew Coward
  • Science reviewer: Gurvan Madec


  • Has the test of the compatibility between bathy and iceshelf draft to be done in NEMO or as an offline tool ? From my experience, it is easier to be done offline because it need many hand edit, flood filling algorithm. This is especially true when the data set for bathy, ice shelf draft and grounding line are not from the same source (it means could be not compatible).
  • Has the tracer part of the ice shelf stuff (sbcisf.F90 and part of trasbc.F90) to be call traisf or sbcisf ? I ask this question, because if Losch 2008 parametrisation or BG03 are activated, we apply the fwf and the heat over more than the surface cell. During the merging with the trunk, I saw that the qsr flux is not put in the surface module but in the tra module, so it could be the same for the ice shelf.


Testing could consider (where appropriate) other configurations in addition to NVTK].

NVTK Tested'''YES/NO'''
Other model configurations'''YES/NO'''
Processor configurations testedYES
If adding new functionality please confirm that the
New code doesn't change results when it is switched off
and ''works'' when switched on

(Answering UNSURE is likely to generate further questions from reviewers.)

'Please add further summary details here'

  • Processor configurations tested: 1x8; 8x1; 1x1; 8x4; 4x8; in ISOMIP configuration (12 steps)
  • Volume conservation test with divisf = .true.:
    • configuration ISOMIP vvl after 133 days: bgvolssh = ~1e-2 km3; bgsaline = ~5e-14 psu ; bgtemper = ~2e-6 °C
    • configuration ISOMIP novvl after 133 days: bgvolssh = ~1e-2 km3; bgsaline = ~2e-6 psu ; bgtemper = ~2e-6 °C
  • Restartability tested with ISOMIP: Works only with key_ldfslp activated even if tra/dynldf_iso set to false !!!! if not key_ldfslp, see #1359: Restartability issue without key_ldfslp ticket.
  • option compatibility with ice shelf cavity:
    • STOP if ln_wave (could be a mask issue somewhere in wave module, not tested)
    • STOP if ln_icebergs (could be a mask issue somewhere in iceberg module, not tested)
    • STOP if sbccpl interface used (could be a mask issue somewhere in cpl module, not tested)
    • STOP if not traadv_cen2 or tvd advection scheme (could be a mask and a surface boundary condition issue, not tested)
    • STOP if not zdfcst or zdftke mixing scheme (could be a mask and surface boundary condition issue, not tested)
    • STOP if not vvl_zstar or zps (as it is working with zstar it should work with ztilde … but not tested)
    • STOP if not hpg_sco (comparison in GYRE formulation with a bump in the bathymetry shows similar result between hpg_sco (modified for ice shelf) and hpg_zps)
    • STOP if not dynspg_flt (exp and ts not tested)

Bit Comparability

Does this change preserve answers in your tested standard configurations (to the last bit) ?NO
Does this change bit compare across various processor configurations. (1xM, Nx1 and MxN are recommended)YES in ISOMIP configuration (not yet tested for the other configurations)
Is this change expected to preserve answers in all possible model configurations?NO
Is this change expected to preserve all diagnostics?
,,''Preserving answers in model runs does not necessarily imply preserved diagnostics. ''

If you answered '''NO''' to any of the above, please provide further details:

  • Which routine(s) are causing the difference?
    • domzgr.F90: For ice shelves, I need to define the the e3t/u/v/w=zt/u/v/w(k+1) - zt/u/v/w(k). If not, dep3w(i,j,k) (could be) .NE. dep3w(i,j+1,k) in the full wet cell beneath the ice shelf.
    • zpshde.F90
    • dynhpg.F90

This has been only tested in GYRE configuration with bump at the bottom and ln_zps activated.

  • Why the changes are not protected by a logical switch or new section-version?
  • What is needed to achieve regression with the previous model release (e.g. a regression branch, hand-edits etc). If this is not possible, explain why not.
    • domzgr.F90: Need to replace definition of gdep3w_0, e3uw_0 and e3vw_0, and e3t_1d and e3w_1d by the reference 3.6 version code
    • zpshde.F90: The easiest thing to do is to keep the header and the variable definition.
    • dynhpg.F90: The easiest thing to do is to replace this routine by the reference 3.6 version.
  • What do you expect to see occur in the test harness jobs?
  • Which diagnostics have you altered and why have they changed? All these changes are mainly due to the fact that the surface layer is not always at jk=1 but jk=mikt(ji,jj).
    • sst and sss are now defined as the temperature of the first wet cell. In case of ice shelf it is the cell at the ice/ocean interface.
    • ocean bottom pressure now takes into account ice load (not tested) and total sal and temperature are corrected to works beneath ice shelf (not tested).
    • in diaharm, I replace mask(:,:,1) by mask_i(:,:) but I am really not sure about it.
    • diaptr should work as before.

System Changes

Does your change alter namelists? YES
Does your change require a change in compiler options? NO

If any of these apply, please document the changes required here……. Namelist changes:

 &namsbc        !   Surface Boundary Condition (surface module)
   nn_isf      = 1         !  0=no isf / 1 = presence of ISF / 2 = bg03 parametrisation / 3 = rnf file for isf / 4 = prescribed melt rate 
 &namsbc_isf    !  Top boundary layer (ISF)
 !              ! file name ! frequency (hours) ! variable ! time interpol. !  clim   ! 'yearly'/ ! weights  ! rotation !
 !              !           !  (if <0  months)  !   name   !    (logical)   !  (T/F)  ! 'monthly' ! filename ! pairing  !
 ! nn_isf == 4
   sn_qisf      = 'rnfisf' ,         -12      ,'sohflisf',    .false.      , .true.  , 'yearly'  ,  ''      ,   ''
   sn_fwfisf    = 'rnfisf' ,         -12      ,'sowflisf',    .false.      , .true.  , 'yearly'  ,  ''      ,   ''
 ! nn_isf == 3
   sn_rnfisf    = 'runoffs' ,         -12      ,'sofwfisf',    .false.      , .true.  , 'yearly'  ,  ''      ,   ''
 ! nn_isf == 2 and 3
   sn_depmax_isf = 'runoffs' ,       -12        ,'sozisfmax' ,   .false.  , .true.  , 'yearly'  ,  ''      ,   ''
   sn_depmin_isf = 'runoffs' ,       -12        ,'sozisfmin' ,   .false.  , .true.  , 'yearly'  ,  ''      ,   ''
 ! nn_isf == 2
   sn_Leff_isf = 'rnfisf' ,       0          ,'Leff'         ,   .false.  , .true.  , 'yearly'  ,  ''      ,   ''
 ! for all case
   ln_divisf   = .true.  ! apply isf melting as a mass flux or in the salinity trend. (maybe I should remove this option as for runoff?)
 ! only for nn_isf = 1 or 2
   rn_gammat0  = 1.0e-4   ! gammat coefficient used in blk formula
   rn_gammas0  = 1.0e-4   ! gammas coefficient used in blk formula
 ! only for nn_isf = 1
   nn_isfblk   =  1       ! 1 ISOMIP ; 2 conservative (3 equation formulation, Jenkins et al. 1991 ??)
   rn_hisf_tbl =  30.      ! thickness of the top boundary layer           (Losh et al. 2008)
                          ! 0 => thickness of the tbl = thickness of the first wet cell
   ln_conserve = .true.   ! conservative case (take into account meltwater advection)
   nn_gammablk = 0        ! 0 = cst Gammat (= gammat/s)
                          ! 1 = velocity dependend Gamma (u* * gammat/s)  (Jenkins et al. 2010)
                          !     if you want to keep the cd as in global config, adjust rn_gammat0 to compensate
                          ! 2 = velocity and stability dependent Gamma    Holland et al. 1999
&nambfr        !   bottom/top friction
   rn_tfri1    =    4.e-4  !  top drag coefficient (linear case)
   rn_tfri2    =    2.5e-3 !  top drag coefficient (non linear case). Minimum coeft if ln_loglayer=T
   rn_tfri2_max =   1.e-1  !  max. top drag coefficient (non linear case and ln_loglayer=T)
   rn_tfeb2    =    0.0    !  top turbulent kinetic energy background  (m2/s2)
   rn_tfrz0    =    3.e-3  !  top roughness [m] if ln_loglayer=T
   ln_tfr2d    = .false.   !  horizontal variation of the top friction coef (read a 2D mask file )
   rn_tfrien   =    50.    !  local multiplying factor of tfr (ln_tfr2d=T)


''Please ''summarize'' any changes in runtime or memory use caused by this change……''

IPR issues

Has the code been wholly (100%) produced by NEMO developers staff working exclusively on NEMO?'''YES/ NO '''

If No:

  • Identify the collaboration agreement details
  • Ensure the code routine header is in accordance with the agreement, (Copyright/Redistribution? etc).Add further details here if required……….