Issues with ssh initialization in iceinit with vvl and AGRIF

Context

• ssh initialization to account for ice mass (hence with si3)

Analysis & Fix

• If starting from rest, ice initialization changes the ocean volume according to the ice mass. In the non qco case, this implies modifying scale factors (in the qco case, scale factors explicitly depends on ssh). There's a problem in iceistate.F90 in the following lines:

         DO jk = 1, jpk
            DO_2D( nn_hls, nn_hls, nn_hls, nn_hls)
               IF( snwice_mass(ji,jj) /= 0._wp ) THEN
                  e3t(ji,jj,jk,Kmm) = e3t_0(ji,jj,jk) * ( 1._wp + ssh(ji,jj,Kmm) * r1_ht_0(ji,jj) * tmask(ji,jj,jk) )
                  e3t(ji,jj,jk,Kbb) = e3t_0(ji,jj,jk) * ( 1._wp + ssh(ji,jj,Kbb) * r1_ht_0(ji,jj) * tmask(ji,jj,jk) )
               ENDIF
            END_2D
         END DO

which does not take into account the case where ssh has been uniformly changed by a constant value. One should remove the IF statement and rather have:

         DO jk = 1, jpk
            DO_2D( nn_hls, nn_hls, nn_hls, nn_hls)
               e3t(ji,jj,jk,Kmm) = e3t_0(ji,jj,jk) * ( 1._wp + ssh(ji,jj,Kmm) * r1_ht_0(ji,jj) * tmask(ji,jj,jk) )
               e3t(ji,jj,jk,Kbb) = e3t_0(ji,jj,jk) * ( 1._wp + ssh(ji,jj,Kbb) * r1_ht_0(ji,jj) * tmask(ji,jj,jk) )
            END_2D
         END DO

• During the initialization of AGRIF child grids, the ssh from the Parent grid is used to set the ssh at least over the ghost points of the child grid. The problem is that the ssh has already been updated to account for the ice mass so that ssh is currently updated twice at the boundaries. One should remove the correction prior interpolating the ssh from parent. Hence in agrif_oce_interp.F90:

   SUBROUTINE interpsshn( ptab, i1, i2, j1, j2, before )
      !!----------------------------------------------------------------------
      !!                  ***  ROUTINE interpsshn  ***
      !!----------------------------------------------------------------------  
      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
      LOGICAL                         , INTENT(in   ) ::   before
      !
      !!----------------------------------------------------------------------  
      !
      IF( before) THEN
#if defined key_si3
         IF (l_ini_child.AND.(.NOT.(ln_rstart .OR. nn_iceini_file == 2))) THEN
            IF( ln_ice_embd ) THEN
               ptab(i1:i2,j1:j2) = ssh(i1:i2,j1:j2,Kmm_a)  &
                    &              + snwice_mass(i1:i2,j1:j2) * r1_rho0
            ELSE
               ptab(i1:i2,j1:j2) = ssh(i1:i2,j1:j2,Kmm_a)  &
                    &              + rsshadj * tmask(i1:i2,j1:j2,1)
            ENDIF
         ELSE
            ptab(i1:i2,j1:j2) = ssh(i1:i2,j1:j2,Kmm_a)
         ENDIF
#else

NB:

• This change in ssh in ice_istate is not intuitive and confusing. This is definitively not the right place to do it, or alternatively, ice initialization should occur before domain initialization (if possible).
• In the case of an uniform volume correction (ln_ice_embd=F), the ice mass is computed over the parent grid. Formally, this is not exact since the global ice mass should be done considering the ice mass over all grids (including child grids).

---