Index: /NEMO/branches/2021/dev_r14608_AGRIF_domcfg/cfgs/SHARED/field_def_nemo-oce.xml
===================================================================
--- /NEMO/branches/2021/dev_r14608_AGRIF_domcfg/cfgs/SHARED/field_def_nemo-oce.xml	(revision 14970)
+++ /NEMO/branches/2021/dev_r14608_AGRIF_domcfg/cfgs/SHARED/field_def_nemo-oce.xml	(revision 14971)
@@ -1108,4 +1108,10 @@
     <field id="vtrd_tot"       long_name="j-trend: total momentum trend before atf"        unit="m/s^2"                        />
     <field id="vtrd_atf"       long_name="j-trend: asselin time filter trend"              unit="m/s^2"                        />
+  </field_group>
+
+  <!-- shared variables available with TOP interface -->
+  <field_group id="top_shared" grid_ref="grid_T_3D">
+    <field id="xeps"           long_name="Broadband light attenuation"                     unit="-"                            />
+    <field id="Heup"           long_name="Euphotic layer depth"                            unit="m"     grid_ref="grid_T_2D"   />
   </field_group>
 
Index: /NEMO/branches/2021/dev_r14608_AGRIF_domcfg/cfgs/SHARED/namelist_ref
===================================================================
--- /NEMO/branches/2021/dev_r14608_AGRIF_domcfg/cfgs/SHARED/namelist_ref	(revision 14970)
+++ /NEMO/branches/2021/dev_r14608_AGRIF_domcfg/cfgs/SHARED/namelist_ref	(revision 14971)
@@ -1234,4 +1234,9 @@
    !                       !           = 2 roughness uses rn_hsri and is weighted by 1-fr_i
    !                       !           = 3 roughness uses rn_hsri and is weighted by 1-MIN(1,4*fr_i)
+   nn_mxlice     =     1   !  mixing under sea ice
+                           !     = 0 No scaling under sea-ice
+                           !     = 1 scaling with constant Ice-ocean roughness (rn_hsri)
+                           !     = 2 scaling with mean sea-ice thickness
+                           !     = 3 scaling with max sea-ice thickness
    nn_bc_surf    =     1   !  surface condition (0/1=Dir/Neum)
    nn_bc_bot     =     1   !  bottom condition (0/1=Dir/Neum)
Index: /NEMO/branches/2021/dev_r14608_AGRIF_domcfg/cfgs/SHARED/namelist_top_ref
===================================================================
--- /NEMO/branches/2021/dev_r14608_AGRIF_domcfg/cfgs/SHARED/namelist_top_ref	(revision 14970)
+++ /NEMO/branches/2021/dev_r14608_AGRIF_domcfg/cfgs/SHARED/namelist_top_ref	(revision 14971)
@@ -101,4 +101,15 @@
 /
 !-----------------------------------------------------------------------
+&namtrc_opt      !  light availability in the water column
+!-----------------------------------------------------------------------
+!              !  file name       ! frequency (hours) ! variable  ! time interp. !  clim  ! 'yearly'/ ! weights  ! rotation ! land/sea mask !
+!              !                  !  (if <0  months)  !   name    !   (logical)  !  (T/F) ! 'monthly' ! filename ! pairing  ! filename      !
+   sn_par      = 'par.orca'       ,     24            , 'fr_par'  ,  .true.      , .true. , 'yearly'  , ''       , ''       , ''
+   cn_dir      = './'        ! root directory for the location of the dynamical files
+   ln_varpar   =  .true.     ! Read PAR from file
+   parlux      =  0.43       ! Fraction of shortwave as PAR
+   light_loc   = 'center'    ! Light location in the water cell ('center', 'integral')
+/
+!-----------------------------------------------------------------------
 &namtrc_dmp      !   passive tracer newtonian damping                   (ln_trcdmp=T)
 !-----------------------------------------------------------------------
Index: /NEMO/branches/2021/dev_r14608_AGRIF_domcfg/src/OCE/SBC/sbcwave.F90
===================================================================
--- /NEMO/branches/2021/dev_r14608_AGRIF_domcfg/src/OCE/SBC/sbcwave.F90	(revision 14970)
+++ /NEMO/branches/2021/dev_r14608_AGRIF_domcfg/src/OCE/SBC/sbcwave.F90	(revision 14971)
@@ -24,4 +24,5 @@
    USE bdy_oce        ! open boundary condition variables
    USE domvvl         ! domain: variable volume layers
+   USE zdf_oce,  ONLY : ln_zdfswm ! Qiao wave enhanced mixing 
    !
    USE iom            ! I/O manager library
@@ -310,4 +311,12 @@
             IF( jp_vsd > 0 )   vt0sd(:,:) = sf_sd(jp_vsd)%fnow(:,:,1) * tmask(:,:,1)  ! 2D meridional Stokes Drift at T point
          ENDIF
+
+         ! Read also wave number if needed, so that it is available in
+         ! coupling routines
+         IF( ln_zdfswm .AND. .NOT. cpl_wnum ) THEN     !==wavenumber==!
+            CALL fld_read( kt, nn_fsbc, sf_wn )             ! read wave parameters from external forcing
+            wnum(:,:) = sf_wn(1)%fnow(:,:,1) * tmask(:,:,1)
+         ENDIF
+ 
          !
          IF( jpfld == 4 .OR. ln_wave_test )   &
@@ -506,6 +515,6 @@
             ENDIF
             !
-            ! 3. Wave number (only needed for Qiao parametrisation, ln_zdfqiao=T)
-            IF( .NOT. cpl_wnum ) THEN
+            ! 3. Wave number (only needed for Qiao parametrisation, ln_zdfswm=T)
+            IF( ln_zdfswm .AND. .NOT. cpl_wnum ) THEN
                ALLOCATE( sf_wn(1), STAT=ierror )           !* allocate and fill sf_wave with sn_wnum
                IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_wn structure' )
Index: /NEMO/branches/2021/dev_r14608_AGRIF_domcfg/src/OCE/ZDF/zdfgls.F90
===================================================================
--- /NEMO/branches/2021/dev_r14608_AGRIF_domcfg/src/OCE/ZDF/zdfgls.F90	(revision 14970)
+++ /NEMO/branches/2021/dev_r14608_AGRIF_domcfg/src/OCE/ZDF/zdfgls.F90	(revision 14971)
@@ -26,4 +26,10 @@
    USE zdfmxl         ! mixed layer
    USE sbcwave , ONLY : hsw   ! significant wave height
+#if defined key_si3
+   USE ice, ONLY: hm_i, h_i
+#endif
+#if defined key_cice
+   USE sbc_ice, ONLY: h_i
+#endif
    !
    USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
@@ -51,4 +57,5 @@
    LOGICAL  ::   ln_length_lim     ! use limit on the dissipation rate under stable stratification (Galperin et al. 1988)
    LOGICAL  ::   ln_sigpsi         ! Activate Burchard (2003) modification for k-eps closure & wave breaking mixing
+   INTEGER  ::   nn_mxlice         ! type of scaling under sea-ice (=0/1/2/3)
    INTEGER  ::   nn_bc_surf        ! surface boundary condition (=0/1)
    INTEGER  ::   nn_bc_bot         ! bottom boundary condition (=0/1)
@@ -224,8 +231,26 @@
       !
       ! adapt roughness where there is sea ice
-      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
-         zhsro(ji,jj) = ( (1._wp-zice_fra(ji,jj)) * zhsro(ji,jj) + zice_fra(ji,jj) * rn_hsri )*tmask(ji,jj,1)  + &
-            &           (1._wp - tmask(ji,jj,1))*rn_hsro
-      END_2D
+      SELECT CASE( nn_mxlice )       ! Type of scaling under sea-ice
+      !
+      CASE( 1 )                      ! scaling with constant sea-ice roughness
+         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+            zhsro(ji,jj) = ( (1._wp-zice_fra(ji,jj)) * zhsro(ji,jj) + zice_fra(ji,jj) * rn_hsri )*tmask(ji,jj,1)  + (1._wp - tmask(ji,jj,1))*rn_hsro
+         END_2D
+         !
+      CASE( 2 )                      ! scaling with mean sea-ice thickness
+#if defined key_si3
+         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+            zhsro(ji,jj) = ( (1._wp-zice_fra(ji,jj)) * zhsro(ji,jj) + zice_fra(ji,jj) * hm_i(ji,jj) )*tmask(ji,jj,1)  + (1._wp - tmask(ji,jj,1))*rn_hsro
+         END_2D
+#endif
+         !
+      CASE( 3 )                      ! scaling with max sea-ice thickness
+#if defined key_si3 || defined key_cice
+         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+            zhsro(ji,jj) = ( (1._wp-zice_fra(ji,jj)) * zhsro(ji,jj) + zice_fra(ji,jj) * MAXVAL(h_i(ji,jj,:)) )*tmask(ji,jj,1)  + (1._wp - tmask(ji,jj,1))*rn_hsro
+         END_2D
+#endif
+         !
+      END SELECT
       !
       DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )  !==  Compute dissipation rate  ==!
@@ -329,4 +354,23 @@
          END_2D
          !
+         IF( ln_isfcav) THEN     ! top boundary   (ocean cavity)
+            DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+               IF( mikt(ji,jj) > 1 )THEN
+                  itop   = mikt(ji,jj)       ! k   top w-point
+                  itopp1 = mikt(ji,jj) + 1   ! k+1 1st w-point below the top one
+                  !                                                ! mask at the
+                  !                                                ocean surface
+                  !                                                points
+                  z_en = MAX( rc02r * ustar2_top(ji,jj), rn_emin ) * ( 1._wp - tmask(ji,jj,1) )
+                  !
+                  ! Dirichlet condition applied at:
+                  !     top level (itop)         &      Just below it (itopp1)
+                  zd_lw(ji,jj,itop) = 0._wp   ;   zd_lw(ji,jj,itopp1) = 0._wp
+                  zd_up(ji,jj,itop) = 0._wp   ;   zd_up(ji,jj,itopp1) = 0._wp
+                  zdiag(ji,jj,itop) = 1._wp   ;   zdiag(ji,jj,itopp1) = 1._wp
+                  en   (ji,jj,itop) = z_en    ;   en   (ji,jj,itopp1) = z_en
+               ENDIF
+            END_2D
+         ENDIF
          !
       CASE ( 1 )             ! Neumann boundary condition (set d(e)/dz)
@@ -352,4 +396,24 @@
          END_2D
          !
+         IF( ln_isfcav) THEN     ! top boundary   (ocean cavity)
+            DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+               IF( mikt(ji,jj) > 1 )THEN
+                  itop   = mikt(ji,jj)       ! k   top w-point
+                  itopp1 = mikt(ji,jj) + 1   ! k+1 1st w-point below the top one
+                  !                                                ! mask at the
+                  !                                                ocean surface
+                  !                                                points
+                  z_en = MAX( rc02r * ustar2_top(ji,jj), rn_emin ) * ( 1._wp - tmask(ji,jj,1) )
+                  !
+                  ! Bottom level Dirichlet condition:
+                  !     Bottom level (ibot)      &      Just above it (ibotm1)
+                  !         Dirichlet            !         Neumann
+                  zd_lw(ji,jj,itop) = 0._wp   !   ! Remove zd_up from zdiag
+                  zdiag(ji,jj,itop) = 1._wp   ;   zdiag(ji,jj,itopp1) = zdiag(ji,jj,itopp1) + zd_up(ji,jj,itopp1)
+                  zd_up(ji,jj,itop) = 0._wp   ;   zd_up(ji,jj,itopp1) = 0._wp
+                  en   (ji,jj,itop) = z_en
+               ENDIF
+            END_2D
+         ENDIF
          !
       END SELECT
@@ -607,4 +671,26 @@
          END_2D
          !
+         IF( ln_isfcav) THEN     ! top boundary   (ocean cavity)
+            DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+               IF ( mikt(ji,jj) > 1 ) THEN
+                  itop   = mikt(ji,jj)       ! k   top w-point
+                  itopp1 = mikt(ji,jj) + 1   ! k+1 1st w-point below the top one
+                  !
+                  zdep(ji,jj) = vkarmn * r_z0_top
+                  psi (ji,jj,itop) = rc0**rpp * en(ji,jj,itop)**rmm *zdep(ji,jj)**rnn
+                  zd_lw(ji,jj,itop) = 0._wp
+                  zd_up(ji,jj,itop) = 0._wp
+                  zdiag(ji,jj,itop) = 1._wp
+                  !
+                  ! Just above last level, Dirichlet condition again (GOTM like)
+                  zdep(ji,jj) = vkarmn * ( r_z0_top + e3t(ji,jj,itopp1,Kmm) )
+                  psi (ji,jj,itopp1) = rc0**rpp * en(ji,jj,itop  )**rmm *zdep(ji,jj)**rnn
+                  zd_lw(ji,jj,itopp1) = 0._wp
+                  zd_up(ji,jj,itopp1) = 0._wp
+                  zdiag(ji,jj,itopp1) = 1._wp
+               END IF
+            END_2D
+         END IF
+         !
       CASE ( 1 )             ! Neumman boundary condition
          !
@@ -631,4 +717,33 @@
             psi(ji,jj,ibotm1) = psi(ji,jj,ibotm1) + zflxb / e3w(ji,jj,ibotm1,Kmm)
          END_2D
+         !
+         IF( ln_isfcav) THEN     ! top boundary   (ocean cavity)
+            DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+               IF ( mikt(ji,jj) > 1 ) THEN
+                  itop   = mikt(ji,jj)       ! k   top w-point
+                  itopp1 = mikt(ji,jj) + 1   ! k+1 1st w-point below the top one
+                  !
+                  ! Bottom level Dirichlet condition:
+                  zdep(ji,jj) = vkarmn * r_z0_top
+                  psi (ji,jj,itop) = rc0**rpp * en(ji,jj,itop)**rmm *zdep(ji,jj)**rnn
+                  !
+                  zd_lw(ji,jj,itop) = 0._wp
+                  zd_up(ji,jj,itop) = 0._wp
+                  zdiag(ji,jj,itop) = 1._wp
+                  !
+                  ! Just below cavity level: Neumann condition with flux
+                  ! injection
+                  zdiag(ji,jj,itopp1) = zdiag(ji,jj,itopp1) + zd_up(ji,jj,itopp1) ! Remove zd_up from zdiag
+                  zd_up(ji,jj,itopp1) = 0._wp
+                  !
+                  ! Set psi vertical flux below cavity:
+                  zdep(ji,jj) = r_z0_top + 0.5_wp*e3t(ji,jj,itopp1,Kmm)
+                  zflxb = rsbc_psi2 * ( p_avm(ji,jj,itop) + p_avm(ji,jj,itopp1))   &
+                     &  * (0.5_wp*(en(ji,jj,itop)+en(ji,jj,itopp1)))**rmm * zdep(ji,jj)**(rnn-1._wp)
+                  psi(ji,jj,itopp1) = psi(ji,jj,itopp1) + zflxb / e3w(ji,jj,itopp1,Kmm)
+               END IF
+            END_2D
+         END IF
+
          !
       END SELECT
@@ -799,5 +914,5 @@
       NAMELIST/namzdf_gls/rn_emin, rn_epsmin, ln_length_lim,       &
          &            rn_clim_galp, ln_sigpsi, rn_hsro, rn_hsri,   &
-         &            rn_crban, rn_charn, rn_frac_hs,              &
+         &            nn_mxlice, rn_crban, rn_charn, rn_frac_hs,   &
          &            nn_bc_surf, nn_bc_bot, nn_z0_met, nn_z0_ice, &
          &            nn_stab_func, nn_clos
@@ -839,5 +954,15 @@
          WRITE(numout,*) '      Type of closure                               nn_clos        = ', nn_clos
          WRITE(numout,*) '      Surface roughness (m)                         rn_hsro        = ', rn_hsro
-         WRITE(numout,*) '      Ice-ocean roughness (used if nn_z0_ice/=0)    rn_hsri        = ', rn_hsri
+         WRITE(numout,*) '      type of scaling under sea-ice                 nn_mxlice      = ', nn_mxlice
+         IF( nn_mxlice == 1 ) &
+            WRITE(numout,*) '      Ice-ocean roughness (used if nn_z0_ice/=0) rn_hsri        = ', rn_hsri
+         SELECT CASE( nn_mxlice )             ! Type of scaling under sea-ice
+            CASE( 0 )   ;   WRITE(numout,*) '   ==>>>   No scaling under sea-ice'
+            CASE( 1 )   ;   WRITE(numout,*) '   ==>>>   scaling with constant sea-ice thickness'
+            CASE( 2 )   ;   WRITE(numout,*) '   ==>>>   scaling with mean     sea-ice thickness'
+            CASE( 3 )   ;   WRITE(numout,*) '   ==>>>   scaling with max      sea-ice thickness'
+            CASE DEFAULT
+               CALL ctl_stop( 'zdf_tke_init: wrong value for nn_mxlice, should be 0,1,2,3 ')
+         END SELECT
          WRITE(numout,*)
       ENDIF
Index: /NEMO/branches/2021/dev_r14608_AGRIF_domcfg/src/OCE/ZDF/zdfphy.F90
===================================================================
--- /NEMO/branches/2021/dev_r14608_AGRIF_domcfg/src/OCE/ZDF/zdfphy.F90	(revision 14970)
+++ /NEMO/branches/2021/dev_r14608_AGRIF_domcfg/src/OCE/ZDF/zdfphy.F90	(revision 14971)
@@ -213,5 +213,5 @@
       IF( ioptio /= 1 )    CALL ctl_stop( 'zdf_phy_init: one and only one vertical diffusion option has to be defined ' )
       IF( ln_isfcav ) THEN
-      IF( ln_zdfric .OR. ln_zdfgls )    CALL ctl_stop( 'zdf_phy_init: zdfric and zdfgls never tested with ice shelves cavities ' )
+      IF( ln_zdfric )      CALL ctl_stop( 'zdf_phy_init: zdfric never tested with ice shelves cavities ' )
       ENDIF
       !                                ! shear production term flag
Index: /NEMO/branches/2021/dev_r14608_AGRIF_domcfg/tests/CANAL/EXPREF/file_def_nemo-oce.xml
===================================================================
--- /NEMO/branches/2021/dev_r14608_AGRIF_domcfg/tests/CANAL/EXPREF/file_def_nemo-oce.xml	(revision 14970)
+++ /NEMO/branches/2021/dev_r14608_AGRIF_domcfg/tests/CANAL/EXPREF/file_def_nemo-oce.xml	(revision 14971)
@@ -26,5 +26,4 @@
 	  <field field_ref="utau"  />
 	  <field field_ref="uoce" />
-	  <field_group group_ref="trendU"  />  
 	</file>
 	
@@ -32,5 +31,4 @@
 	  <field field_ref="vtau"  />
 	  <field field_ref="voce" />
-	  <field_group group_ref="trendV"  />  
 	</file>
 	
