Index: /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/cfgs/SHARED/namelist_ice_ref
===================================================================
--- /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/cfgs/SHARED/namelist_ice_ref	(revision 12893)
+++ /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/cfgs/SHARED/namelist_ice_ref	(revision 12894)
@@ -102,6 +102,6 @@
 &namdyn_adv     !   Ice advection
 !------------------------------------------------------------------------------
-   ln_adv_Pra       = .true.         !  Advection scheme (Prather)
-   ln_adv_UMx       = .false.          !  Advection scheme (Ultimate-Macho)
+   ln_adv_Pra       = .true.          !  Advection scheme (Prather)
+   ln_adv_UMx       = .false.         !  Advection scheme (Ultimate-Macho)
       nn_UMx        =   5             !     order of the scheme for UMx (1-5 ; 20=centered 2nd order)
 /
@@ -123,4 +123,7 @@
    ln_cndflx        = .false.         !  Use conduction flux as surface boundary conditions (i.e. for Jules coupling)
       ln_cndemulate = .false.         !     emulate conduction flux (if not provided in the inputs)
+   nn_qtrice        = 0               !  Solar flux transmitted thru the surface scattering layer:
+                                      !     = 0  Grenfell and Maykut 1977 (depends on cloudiness and is 0 when there is snow) 
+                                      !     = 1  Lebrun 2019 (equals 0.3 anytime with different melting/dry snw conductivities)
 /
 !------------------------------------------------------------------------------
@@ -142,5 +145,8 @@
    rn_cnd_s         =   0.31          !  thermal conductivity of the snow (0.31 W/m/K, Maykut and Untersteiner, 1971)
                                       !     Obs: 0.1-0.5 (Lecomte et al, JAMES 2013)
-   rn_kappa_i       =   1.0           !  radiation attenuation coefficient in sea ice [1/m]
+   rn_kappa_i       =   1.0           !  radiation attenuation coefficient in sea ice                     [1/m]
+   rn_kappa_s       =  10.0           !  nn_qtrice = 0: radiation attenuation coefficient in snow         [1/m]
+   rn_kappa_smlt    =   7.0           !  nn_qtrice = 1: radiation attenuation coefficient in melting snow [1/m]
+   rn_kappa_sdry    =  10.0           !                 radiation attenuation coefficient in dry snow     [1/m]
 /
 !------------------------------------------------------------------------------
Index: /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/doc/namelists/namthd_zdf
===================================================================
--- /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/doc/namelists/namthd_zdf	(revision 12893)
+++ /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/doc/namelists/namthd_zdf	(revision 12894)
@@ -7,4 +7,7 @@
    rn_cnd_s         =   0.31          !  thermal conductivity of the snow (0.31 W/m/K, Maykut and Untersteiner, 1971)
                                       !     Obs: 0.1-0.5 (Lecomte et al, JAMES 2013)
-   rn_kappa_i       =   1.0           !  radiation attenuation coefficient in sea ice [1/m]
+   rn_kappa_i       =   1.0           !  radiation attenuation coefficient in sea ice                     [1/m]
+   rn_kappa_s       =  10.0           !  nn_qtrice = 0: radiation attenuation coefficient in snow         [1/m]
+   rn_kappa_smlt    =   7.0           !  nn_qtrice = 1: radiation attenuation coefficient in melting snow [1/m]
+   rn_kappa_sdry    =  10.0           !                 radiation attenuation coefficient in dry snow     [1/m]
 /
Index: /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/ice.F90
===================================================================
--- /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/ice.F90	(revision 12893)
+++ /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/ice.F90	(revision 12894)
@@ -176,4 +176,5 @@
    !                                      !   = 1  Average N(cat) fluxes then redistribute over the N(cat) ice using T-ice and albedo sensitivity
    !                                      !   = 2  Redistribute a single flux over categories
+                                          ! -- icethd_zdf -- !
    LOGICAL , PUBLIC ::   ln_cndflx        !: use conduction flux as surface boundary condition (instead of qsr and qns) 
    LOGICAL , PUBLIC ::   ln_cndemulate    !: emulate conduction flux (if not provided) 
@@ -182,10 +183,16 @@
    INTEGER, PUBLIC, PARAMETER ::   np_cnd_ON  = 1  !: forcing from conduction flux (SM0L) (compute qcn and qsr_tr via sbcblk.F90 or sbccpl.F90)
    INTEGER, PUBLIC, PARAMETER ::   np_cnd_EMU = 2  !: emulate conduction flux via icethd_zdf.F90 (BL99) (1st round compute qcn and qsr_tr, 2nd round use it)
-
+   INTEGER, PUBLIC ::   nn_qtrice         !: Solar flux transmitted thru the surface scattering layer:
+   !                                      !   = 0  Grenfell and Maykut 1977 (depends on cloudiness and is 0 when there is snow) 
+   !                                      !   = 1  Lebrun 2019 (equals 0.3 anytime with different melting/dry snw conductivities)
+   !
    !                                     !!** ice-vertical diffusion namelist (namthd_zdf) **
    LOGICAL , PUBLIC ::   ln_cndi_U64      !: thermal conductivity: Untersteiner (1964)
    LOGICAL , PUBLIC ::   ln_cndi_P07      !: thermal conductivity: Pringle et al (2007)
-   REAL(wp), PUBLIC ::   rn_kappa_i       !: coef. for the extinction of radiation Grenfell et al. (2006) [1/m]
    REAL(wp), PUBLIC ::   rn_cnd_s         !: thermal conductivity of the snow [W/m/K]   
+   REAL(wp), PUBLIC ::   rn_kappa_i       !: coef. for the extinction of radiation in sea ice, Grenfell et al. (2006) [1/m]
+   REAL(wp), PUBLIC ::   rn_kappa_s       !: coef. for the extinction of radiation in snw (nn_qtrice=0) [1/m]
+   REAL(wp), PUBLIC ::   rn_kappa_smlt    !: coef. for the extinction of radiation in melt snw (nn_qtrice=1) [1/m]
+   REAL(wp), PUBLIC ::   rn_kappa_sdry    !: coef. for the extinction of radiation in dry  snw (nn_qtrice=1) [1/m]
 
    !                                     !!** ice-salinity namelist (namthd_sal) **
Index: /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icesbc.F90
===================================================================
--- /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icesbc.F90	(revision 12893)
+++ /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icesbc.F90	(revision 12894)
@@ -276,5 +276,5 @@
       INTEGER ::   ios, ioptio   ! Local integer
       !!
-      NAMELIST/namsbc/ rn_cio, rn_snwblow, nn_snwfra, nn_flxdist, ln_cndflx, ln_cndemulate
+      NAMELIST/namsbc/ rn_cio, nn_snwfra, rn_snwblow, nn_flxdist, ln_cndflx, ln_cndemulate, nn_qtrice
       !!-------------------------------------------------------------------
       !
@@ -292,10 +292,13 @@
          WRITE(numout,*) '~~~~~~~~~~~~~~~~'
          WRITE(numout,*) '   Namelist namsbc:'
-         WRITE(numout,*) '      drag coefficient for oceanic stress              rn_cio        = ', rn_cio
-         WRITE(numout,*) '      fraction of ice covered by snow (options 0,1,2)  nn_snwfra     = ', nn_snwfra
-         WRITE(numout,*) '      coefficient for ice-lead partition of snowfall   rn_snwblow    = ', rn_snwblow
-         WRITE(numout,*) '      Multicategory heat flux formulation              nn_flxdist    = ', nn_flxdist
-         WRITE(numout,*) '      Use conduction flux as surface condition         ln_cndflx     = ', ln_cndflx
-         WRITE(numout,*) '         emulate conduction flux                       ln_cndemulate = ', ln_cndemulate
+         WRITE(numout,*) '      drag coefficient for oceanic stress                       rn_cio        = ', rn_cio
+         WRITE(numout,*) '      fraction of ice covered by snow (options 0,1,2)           nn_snwfra     = ', nn_snwfra
+         WRITE(numout,*) '      coefficient for ice-lead partition of snowfall            rn_snwblow    = ', rn_snwblow
+         WRITE(numout,*) '      Multicategory heat flux formulation                       nn_flxdist    = ', nn_flxdist
+         WRITE(numout,*) '      Use conduction flux as surface condition                  ln_cndflx     = ', ln_cndflx
+         WRITE(numout,*) '         emulate conduction flux                                ln_cndemulate = ', ln_cndemulate
+         WRITE(numout,*) '      solar flux transmitted thru the surface scattering layer  nn_qtrice     = ', nn_qtrice
+         WRITE(numout,*) '         = 0  Grenfell and Maykut 1977'
+         WRITE(numout,*) '         = 1  Lebrun 2019'
       ENDIF
       !
Index: /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icethd_zdf.F90
===================================================================
--- /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icethd_zdf.F90	(revision 12893)
+++ /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icethd_zdf.F90	(revision 12894)
@@ -85,5 +85,6 @@
       INTEGER  ::   ios, ioptio   ! Local integer
       !!
-      NAMELIST/namthd_zdf/ ln_zdf_BL99, ln_cndi_U64, ln_cndi_P07, rn_cnd_s, rn_kappa_i
+      NAMELIST/namthd_zdf/ ln_zdf_BL99, ln_cndi_U64, ln_cndi_P07, rn_cnd_s, &
+         &                 rn_kappa_i, rn_kappa_s, rn_kappa_smlt, rn_kappa_sdry
       !!-------------------------------------------------------------------
       !
@@ -101,9 +102,12 @@
          WRITE(numout,*) '~~~~~~~~~~~~~~~~'
          WRITE(numout,*) '   Namelist namthd_zdf:'
-         WRITE(numout,*) '      Bitz and Lipscomb (1999) formulation                    ln_zdf_BL99  = ', ln_zdf_BL99
-         WRITE(numout,*) '      thermal conductivity in the ice (Untersteiner 1964)     ln_cndi_U64  = ', ln_cndi_U64
-         WRITE(numout,*) '      thermal conductivity in the ice (Pringle et al 2007)    ln_cndi_P07  = ', ln_cndi_P07
-         WRITE(numout,*) '      thermal conductivity in the snow                        rn_cnd_s     = ', rn_cnd_s
-         WRITE(numout,*) '      extinction radiation parameter in sea ice               rn_kappa_i   = ', rn_kappa_i
+         WRITE(numout,*) '      Bitz and Lipscomb (1999) formulation                      ln_zdf_BL99   = ', ln_zdf_BL99
+         WRITE(numout,*) '      thermal conductivity in the ice (Untersteiner 1964)       ln_cndi_U64   = ', ln_cndi_U64
+         WRITE(numout,*) '      thermal conductivity in the ice (Pringle et al 2007)      ln_cndi_P07   = ', ln_cndi_P07
+         WRITE(numout,*) '      thermal conductivity in the snow                          rn_cnd_s      = ', rn_cnd_s
+         WRITE(numout,*) '      extinction radiation parameter in sea ice                 rn_kappa_i    = ', rn_kappa_i
+         WRITE(numout,*) '      extinction radiation parameter in snw      (nn_qtrice=0)  rn_kappa_s    = ', rn_kappa_s
+         WRITE(numout,*) '      extinction radiation parameter in melt snw (nn_qtrice=1)  rn_kappa_smlt = ', rn_kappa_smlt
+         WRITE(numout,*) '      extinction radiation parameter in dry  snw (nn_qtrice=1)  rn_kappa_sdry = ', rn_kappa_sdry
       ENDIF
       !
Index: /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icethd_zdf_bl99.F90
===================================================================
--- /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icethd_zdf_bl99.F90	(revision 12893)
+++ /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icethd_zdf_bl99.F90	(revision 12894)
@@ -90,5 +90,4 @@
       REAL(wp) ::   zgamma    =  18009._wp    ! for specific heat
       REAL(wp) ::   zbeta     =  0.117_wp     ! for thermal conductivity (could be 0.13)
-      REAL(wp) ::   zraext_s  =  10._wp       ! extinction coefficient of radiation in the snow
       REAL(wp) ::   zkimin    =  0.10_wp      ! minimum ice thermal conductivity
       REAL(wp) ::   ztsu_err  =  1.e-5_wp     ! range around which t_su is considered at 0C 
@@ -101,4 +100,5 @@
       REAL(wp) ::   zhfx_err, zdq             ! diag errors on heat
       !
+      REAL(wp), DIMENSION(jpij) ::   zraext_s     ! extinction coefficient of radiation in the snow
       REAL(wp), DIMENSION(jpij) ::   ztsub        ! surface temperature at previous iteration
       REAL(wp), DIMENSION(jpij) ::   zh_i, z1_h_i ! ice layer thickness
@@ -148,4 +148,15 @@
       ! 1) Initialization
       !------------------
+      !
+      ! extinction radiation in the snow
+      IF    ( nn_qtrice == 0 ) THEN   ! constant 
+         zraext_s(1:npti) = rn_kappa_s
+      ELSEIF( nn_qtrice == 1 ) THEN   ! depends on melting/freezing conditions
+         WHERE( t_su_1d(1:npti) < rt0 )   ;   zraext_s(1:npti) = rn_kappa_sdry   ! no surface melting
+         ELSEWHERE                        ;   zraext_s(1:npti) = rn_kappa_smlt   !    surface melting
+         END WHERE
+      ENDIF
+      !
+      ! layers thicknesses
       DO ji = 1, npti
          zh_s(ji) = MAX( zhs_ssl , h_s_1d(ji) ) * r1_nlay_s ! set a minimum snw thickness for conduction
@@ -183,5 +194,5 @@
          DO ji = 1, npti
             !                             ! radiation transmitted below the layer-th snow layer
-            zradtr_s(ji,jk) = zradtr_s(ji,0) * EXP( - zraext_s * MAX( 0._wp, zh_s(ji) * REAL(jk) - zhs_ssl ) )
+            zradtr_s(ji,jk) = zradtr_s(ji,0) * EXP( - zraext_s(ji) * MAX( 0._wp, zh_s(ji) * REAL(jk) - zhs_ssl ) )
             !                             ! radiation absorbed by the layer-th snow layer
             zradab_s(ji,jk) = zradtr_s(ji,jk-1) - zradtr_s(ji,jk)
@@ -193,4 +204,9 @@
          DO ji = 1, npti
             !                             ! radiation transmitted below the layer-th ice layer
+            zradtr_i(ji,jk) =           za_s_fra(ji)   * zradtr_s(ji,nlay_s)                       &   ! part covered by snow
+               &                                       * EXP( - rn_kappa_i * zh_i(ji) * REAL(jk) ) &
+               &            + ( 1._wp - za_s_fra(ji) ) * qtr_ice_top_1d(ji)                        &   ! part snow free
+               &                                       * EXP( - rn_kappa_i * MAX( 0._wp, zh_i(ji) * REAL(jk) - zhi_ssl ) )
+            
             zradtr_i(ji,jk) = zradtr_i(ji,0) * EXP( - rn_kappa_i * MAX( 0._wp, zh_i(ji) * REAL(jk) - zhi_ssl ) )
             !                             ! radiation absorbed by the layer-th ice layer
Index: /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/OCE/SBC/sbcblk.F90
===================================================================
--- /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/OCE/SBC/sbcblk.F90	(revision 12893)
+++ /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/OCE/SBC/sbcblk.F90	(revision 12894)
@@ -46,5 +46,5 @@
    USE lib_fortran    ! to use key_nosignedzero
 #if defined key_si3
-   USE ice     , ONLY :   u_ice, v_ice, jpl, a_i_b, at_i_b, t_su, rn_cnd_s, hfx_err_dif
+   USE ice     , ONLY :   u_ice, v_ice, jpl, a_i_b, at_i_b, t_su, rn_cnd_s, hfx_err_dif, nn_qtrice
    USE icevar         ! for CALL ice_var_snwblow
 #endif
@@ -917,16 +917,26 @@
       END DO
 
-      ! --- shortwave radiation transmitted below the surface (W/m2, see Grenfell Maykut 77) --- !
-      ztri(:,:) = 0.18 * ( 1.0 - cloud_fra(:,:) ) + 0.35 * cloud_fra(:,:)  ! surface transmission when hi>10cm
-      !
-      DO jl = 1, jpl
-         WHERE    ( phs(:,:,jl) <= 0._wp .AND. phi(:,:,jl) <  0.1_wp )     ! linear decrease from hi=0 to 10cm  
-            qtr_ice_top(:,:,jl) = qsr_ice(:,:,jl) * ( ztri(:,:) + ( 1._wp - ztri(:,:) ) * ( 1._wp - phi(:,:,jl) * 10._wp ) )
-         ELSEWHERE( phs(:,:,jl) <= 0._wp .AND. phi(:,:,jl) >= 0.1_wp )     ! constant (ztri) when hi>10cm
-            qtr_ice_top(:,:,jl) = qsr_ice(:,:,jl) * ztri(:,:)
-         ELSEWHERE                                                         ! zero when hs>0
-            qtr_ice_top(:,:,jl) = 0._wp 
-         END WHERE
-      ENDDO
+      ! --- shortwave radiation transmitted thru the surface scattering layer (W/m2) --- !
+      IF( nn_qtrice == 0 ) THEN
+         ! formulation derived from Grenfell and Maykut (1977), where transmission rate
+         !    1) depends on cloudiness
+         !    2) is 0 when there is any snow
+         !    3) tends to 1 for thin ice
+         ztri(:,:) = 0.18 * ( 1.0 - cloud_fra(:,:) ) + 0.35 * cloud_fra(:,:)  ! surface transmission when hi>10cm
+         DO jl = 1, jpl
+            WHERE    ( phs(:,:,jl) <= 0._wp .AND. phi(:,:,jl) <  0.1_wp )     ! linear decrease from hi=0 to 10cm  
+               qtr_ice_top(:,:,jl) = qsr_ice(:,:,jl) * ( ztri(:,:) + ( 1._wp - ztri(:,:) ) * ( 1._wp - phi(:,:,jl) * 10._wp ) )
+            ELSEWHERE( phs(:,:,jl) <= 0._wp .AND. phi(:,:,jl) >= 0.1_wp )     ! constant (ztri) when hi>10cm
+               qtr_ice_top(:,:,jl) = qsr_ice(:,:,jl) * ztri(:,:)
+            ELSEWHERE                                                         ! zero when hs>0
+               qtr_ice_top(:,:,jl) = 0._wp 
+            END WHERE
+         ENDDO
+      ELSEIF( nn_qtrice == 1 ) THEN
+         ! formulation is derived from the thesis of M. Lebrun (2019).
+         !    It represents the best fit using several sets of observations
+         !    It comes with snow conductivities adapted to freezing/melting conditions (see icethd_zdf_bl99.F90)
+         qtr_ice_top(:,:,:) = 0.3_wp * qsr_ice(:,:,:)
+      ENDIF
       !
 
Index: /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/OCE/SBC/sbccpl.F90
===================================================================
--- /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/OCE/SBC/sbccpl.F90	(revision 12893)
+++ /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/OCE/SBC/sbccpl.F90	(revision 12894)
@@ -2111,22 +2111,32 @@
       IF( .NOT.ln_cndflx ) THEN                              !==  No conduction flux as surface forcing  ==!
          !
-         !                    ! ===> used prescribed cloud fraction representative for polar oceans in summer (0.81)
-         !                    !      should be real cloud fraction instead (as in the bulk) but needs to be read from atm.
-         ztri(:,:) = 0.18 * ( 1.0 - zcloud_fra(:,:) ) + 0.35 * zcloud_fra(:,:)  ! surface transmission when hi>10cm (Grenfell Maykut 77)
-         !
-         DO jl = 1, jpl
-            WHERE    ( phs(:,:,jl) <= 0._wp .AND. phi(:,:,jl) <  0.1_wp )     ! linear decrease from hi=0 to 10cm  
-               zqtr_ice_top(:,:,jl) = qsr_ice(:,:,jl) * ( ztri(:,:) + ( 1._wp - ztri(:,:) ) * ( 1._wp - phi(:,:,jl) * 10._wp ) )
-            ELSEWHERE( phs(:,:,jl) <= 0._wp .AND. phi(:,:,jl) >= 0.1_wp )     ! constant (ztri) when hi>10cm
-               zqtr_ice_top(:,:,jl) = qsr_ice(:,:,jl) * ztri(:,:)
-            ELSEWHERE                                                         ! zero when hs>0
-               zqtr_ice_top(:,:,jl) = 0._wp
-            END WHERE
-         ENDDO
+         IF( nn_qtrice == 0 ) THEN
+            ! formulation derived from Grenfell and Maykut (1977), where transmission rate
+            !    1) depends on cloudiness
+            !       ! ===> used prescribed cloud fraction representative for polar oceans in summer (0.81)
+            !       !      should be real cloud fraction instead (as in the bulk) but needs to be read from atm.
+            !    2) is 0 when there is any snow
+            !    3) tends to 1 for thin ice
+            ztri(:,:) = 0.18 * ( 1.0 - zcloud_fra(:,:) ) + 0.35 * zcloud_fra(:,:)  ! surface transmission when hi>10cm
+            DO jl = 1, jpl
+               WHERE    ( phs(:,:,jl) <= 0._wp .AND. phi(:,:,jl) <  0.1_wp )       ! linear decrease from hi=0 to 10cm  
+                  zqtr_ice_top(:,:,jl) = zqsr_ice(:,:,jl) * ( ztri(:,:) + ( 1._wp - ztri(:,:) ) * ( 1._wp - phi(:,:,jl) * 10._wp ) )
+               ELSEWHERE( phs(:,:,jl) <= 0._wp .AND. phi(:,:,jl) >= 0.1_wp )       ! constant (ztri) when hi>10cm
+                  zqtr_ice_top(:,:,jl) = zqsr_ice(:,:,jl) * ztri(:,:)
+               ELSEWHERE                                                           ! zero when hs>0
+                  zqtr_ice_top(:,:,jl) = 0._wp 
+               END WHERE
+            ENDDO
+         ELSEIF( nn_qtrice == 1 ) THEN
+            ! formulation is derived from the thesis of M. Lebrun (2019).
+            !    It represents the best fit using several sets of observations
+            !    It comes with snow conductivities adapted to freezing/melting conditions (see icethd_zdf_bl99.F90)
+            zqtr_ice_top(:,:,:) = 0.3_wp * zqsr_ice(:,:,:)
+         ENDIF
          !     
       ELSEIF( ln_cndflx .AND. .NOT.ln_cndemulate ) THEN      !==  conduction flux as surface forcing  ==!
          !
-         !                    ! ===> here we must receive the qtr_ice_top array from the coupler
-         !                           for now just assume zero (fully opaque ice)
+         !          ! ===> here we must receive the qtr_ice_top array from the coupler
+         !                 for now just assume zero (fully opaque ice)
          zqtr_ice_top(:,:,:) = 0._wp
          !
