Index: /branches/2017/dev_CNRS_2017/NEMOGCM/CONFIG/TEST_CASES/SAS_BIPER/EXP00/README
===================================================================
--- /branches/2017/dev_CNRS_2017/NEMOGCM/CONFIG/TEST_CASES/SAS_BIPER/EXP00/README	(revision 8961)
+++ /branches/2017/dev_CNRS_2017/NEMOGCM/CONFIG/TEST_CASES/SAS_BIPER/EXP00/README	(revision 8962)
@@ -1,10 +1,24 @@
 -------
-Context
+Purpose
 -------
-The aim of this demonstration case is to study the advection of a patch of ice in a bi-periodic channel with a slab ocean (i.e. 1 ocean layer which does not interact with sea-ice).
+This demonstration case can serve to different purposes:
 
+1) The main one: study of the advection of a patch of ice in a bi-periodic channel with a slab ocean (i.e. 1 ocean layer)
+              ==> set ln_icethd=false in namelist_ice_cfg
+
+2) Secondary: study of ice thermodynamics in the same basin
+              ==> set ln_icedyn=false in namelist_ice_cfg
+              ==> eventually change usrdef_sbc.F90 in MY_SRC to fit your needs
+
+-----------
+Experiments
+-----------
 Two experiments can be configured:
+
 a) A simple channel at 3km horizontal resolution (slightly variable => +- 10%)
+              ==> remove key_agrif in cpp keys
+
 b) Same as a) but with an Agrif zoom in the center 1:3, and 3 ghost cells (defined in par_oce.F90)
+              ==> add key_agrif in cpp keys
 
 ----------
@@ -13,12 +27,12 @@
 Be careful, the model can only run on 1 processor in the N-S dimension
 
-Compile and run the model once to get a mesh_mask.nc file with the following command (note; remove key_agrif in cpp keys for experiment a):
+a) Compile and run the model once to get a mesh_mask.nc file with the following command:
 ../../../makenemo -a TEST_CASES -n SAS_BIPER -m X64_ADA -j 4
 poe ./opa -procs 1
 
-Create the initial condition file for sea-ice (initice.nc) by running this python script: 
+b) Create the initial condition file for sea-ice (initice.nc) by running this python script: 
 python ./make_INITICE.py
 
-Run the model a second time
+c) Run the model a second time
 poe ./opa -procs 1
 
@@ -26,6 +40,6 @@
 What to look at
 ---------------
-One can test
-
+In case of purpose 1, One can test
+---
 a) the advection scheme: Ultimate-Macho (ln_adv_UMx=T) versus Prather (ln_adv_Pra=T)
    for a square (ice concentration) or a gaussian (ice volume)
@@ -36,4 +50,8 @@
 b) the advection through an agrif zoom 1:3
 
+In case of purpose 2, one can test conservation of properties:
+---
+c) ice should not change at all if surface fluxes = 0 and SST = freezing temperature
+
 --------------------------------------
 Interpretation of the results (remarks)
Index: /branches/2017/dev_CNRS_2017/NEMOGCM/CONFIG/TEST_CASES/SAS_BIPER/EXP00/file_def_nemo-lim.xml
===================================================================
--- /branches/2017/dev_CNRS_2017/NEMOGCM/CONFIG/TEST_CASES/SAS_BIPER/EXP00/file_def_nemo-lim.xml	(revision 8961)
+++ /branches/2017/dev_CNRS_2017/NEMOGCM/CONFIG/TEST_CASES/SAS_BIPER/EXP00/file_def_nemo-lim.xml	(revision 8962)
@@ -30,4 +30,5 @@
         <field field_ref="iceshe"           name="sishea" />
         <field field_ref="icestr"           name="sistre" />
+        <field field_ref="sst_m"            name="sst_m" />
         <!-- categories -->
         <field field_ref="snwthic_cat"      name="snthicat"/>
Index: /branches/2017/dev_CNRS_2017/NEMOGCM/CONFIG/TEST_CASES/SAS_BIPER/MY_SRC/usrdef_sbc.F90
===================================================================
--- /branches/2017/dev_CNRS_2017/NEMOGCM/CONFIG/TEST_CASES/SAS_BIPER/MY_SRC/usrdef_sbc.F90	(revision 8961)
+++ /branches/2017/dev_CNRS_2017/NEMOGCM/CONFIG/TEST_CASES/SAS_BIPER/MY_SRC/usrdef_sbc.F90	(revision 8962)
@@ -25,5 +25,4 @@
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
    USE lib_fortran     ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
-   USE wrk_nemo
 
    IMPLICIT NONE
@@ -93,5 +92,5 @@
    END SUBROUTINE usrdef_sbc_ice_tau
 
-   SUBROUTINE usrdef_sbc_ice_flx( kt )
+   SUBROUTINE usrdef_sbc_ice_flx( kt, phs, phi )
       !!---------------------------------------------------------------------
       !!                     ***  ROUTINE usrdef_sbc_ice_flx  ***
@@ -99,8 +98,11 @@
       !! ** Purpose :   provide the surface boundary (flux) condition over sea-ice
       !!---------------------------------------------------------------------
-      REAL(wp), DIMENSION(:,:), POINTER ::   zsnw       ! snw distribution after wind blowing
       INTEGER, INTENT(in) ::   kt   ! ocean time step
+      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   phs    ! snow thickness
+      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   phi    ! ice thickness
+      !!
+      REAL(wp) ::   zfr1, zfr2                 ! local variables
+      REAL(wp), DIMENSION(jpi,jpj) ::   zsnw   ! snw distribution after wind blowing
       !!---------------------------------------------------------------------
-      CALL wrk_alloc( jpi,jpj, zsnw )
       !
       IF( kt==nit000 .AND. lwp)   WRITE(numout,*)' usrdef_sbc_ice : SAS_BIPER case: NO flux forcing'
@@ -133,13 +135,16 @@
       qsr_tot (:,:) = at_i_b(:,:) * qsr_oce(:,:) + SUM( a_i_b(:,:,:) * qsr_ice(:,:,:), dim=3 )
 
-      !--------------------------------------------------------------------
-      ! FRACTIONs of net shortwave radiation which is not absorbed in the
-      ! thin surface layer and penetrates inside the ice cover
-      ! ( Maykut and Untersteiner, 1971 ; Ebert and Curry, 1993 )
-      fr1_i0(:,:) = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice )
-      fr2_i0(:,:) = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice )
-
-      CALL wrk_dealloc( jpi,jpj, zsnw )
-
+      ! --- shortwave radiation transmitted below the surface (W/m2, see Grenfell Maykut 77) --- !
+      zfr1 = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice )            ! transmission when hi>10cm
+      zfr2 = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice )            ! zfr2 such that zfr1 + zfr2 to equal 1
+      !
+      WHERE    ( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) <  0.1_wp )       ! linear decrease from hi=0 to 10cm  
+         qsr_ice_tr(:,:,:) = qsr_ice(:,:,:) * ( zfr1 + zfr2 * ( 1._wp - phi(:,:,:) * 10._wp ) )
+      ELSEWHERE( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) >= 0.1_wp )       ! constant (zfr1) when hi>10cm
+         qsr_ice_tr(:,:,:) = qsr_ice(:,:,:) * zfr1
+      ELSEWHERE                                                         ! zero when hs>0
+         qsr_ice_tr(:,:,:) = 0._wp 
+      END WHERE
+          
    END SUBROUTINE usrdef_sbc_ice_flx
 
Index: /branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/LIM_SRC_3/iceforcing.F90
===================================================================
--- /branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/LIM_SRC_3/iceforcing.F90	(revision 8961)
+++ /branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/LIM_SRC_3/iceforcing.F90	(revision 8962)
@@ -139,5 +139,5 @@
       !
       CASE( jp_usr )              !--- user defined formulation
-                                  CALL usrdef_sbc_ice_flx( kt )
+                                  CALL usrdef_sbc_ice_flx( kt, h_s, h_i )
       CASE( jp_blk )              !--- bulk formulation
                                   CALL blk_ice_flx    ( t_su, h_s, h_i, alb_ice )    ! 
Index: /branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk.F90
===================================================================
--- /branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk.F90	(revision 8961)
+++ /branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk.F90	(revision 8962)
@@ -843,6 +843,5 @@
       REAL(wp) ::   zcoef_dqlw, zcoef_dqla   !   -      -
       REAL(wp) ::   zztmp, z1_lsub           !   -      -
-      REAL(wp) ::   zfrqsr_tr_i              ! sea ice shortwave fraction transmitted below through the ice
-      REAL(wp) ::   zfr1, zfr2, zfac         ! local variables
+      REAL(wp) ::   zfr1, zfr2               ! local variables
       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z_qlw         ! long wave heat flux over ice
       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z_qsb         ! sensible  heat flux over ice
@@ -951,32 +950,15 @@
       END DO
 
-      ! --- absorbed and transmitted shortwave radiation (W/m2) --- !
-      !
-      ! former coding was
-      !     fr1_i0(:,:) = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice )
-      !     fr2_i0(:,:) = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice )
-
-      ! --- surface transmission parameter (i0, Grenfell Maykut 77) --- !
-      zfr1 = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice )            !   standard value
-      zfr2 = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice )            !   zfr2 such that zfr1 + zfr2 to equal 1
-
-      qsr_ice_tr(:,:,:) = 0._wp
-
-      DO jl = 1, jpl
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               !
-               zfac = MAX( 0._wp , 1._wp - phi(ji,jj,jl) * 10._wp )     !   linear weighting factor: =0 for phi=0, =1 for phi = 0.1
-               zfrqsr_tr_i = zfr1 + zfac * zfr2                         !   below 10 cm, linearly increase zfrqsr_tr_i until 1 at zero thickness
-               !
-               IF ( phs(ji,jj,jl) <= 0._wp ) THEN   ;    zfrqsr_tr_i  = zfr1 + zfac * zfr2
-               ELSE                                 ;    zfrqsr_tr_i  = 0._wp                                  !   snow fully opaque
-               ENDIF
-               !
-               qsr_ice_tr(ji,jj,jl) = zfrqsr_tr_i * qsr_ice(ji,jj,jl)   !   transmitted solar radiation 
-               !
-            END DO
-         END DO
-      END DO
+      ! --- shortwave radiation transmitted below the surface (W/m2, see Grenfell Maykut 77) --- !
+      zfr1 = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice )            ! transmission when hi>10cm
+      zfr2 = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice )            ! zfr2 such that zfr1 + zfr2 to equal 1
+      !
+      WHERE    ( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) <  0.1_wp )       ! linear decrease from hi=0 to 10cm  
+         qsr_ice_tr(:,:,:) = qsr_ice(:,:,:) * ( zfr1 + zfr2 * ( 1._wp - phi(:,:,:) * 10._wp ) )
+      ELSEWHERE( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) >= 0.1_wp )       ! constant (zfr1) when hi>10cm
+         qsr_ice_tr(:,:,:) = qsr_ice(:,:,:) * zfr1
+      ELSEWHERE                                                         ! zero when hs>0
+         qsr_ice_tr(:,:,:) = 0._wp 
+      END WHERE
       !
       IF(ln_ctl) THEN
Index: /branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90
===================================================================
--- /branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 8961)
+++ /branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 8962)
@@ -47,5 +47,4 @@
    USE iom            ! NetCDF library
    USE lib_mpp        ! distribued memory computing library
-   USE wrk_nemo       ! work arrays
    USE timing         ! Timing
    USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
@@ -244,5 +243,5 @@
       INTEGER ::   jn          ! dummy loop index
       INTEGER ::   ios, inum   ! Local integer
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zacs, zaos
+      REAL(wp), DIMENSION(jpi,jpj) ::   zacs, zaos
       !!
       NAMELIST/namsbc_cpl/  sn_snd_temp  , sn_snd_alb   , sn_snd_thick, sn_snd_crt   , sn_snd_co2  ,   & 
@@ -258,6 +257,4 @@
       IF( nn_timing == 1 )   CALL timing_start('sbc_cpl_init')
       !
-      CALL wrk_alloc( jpi,jpj,   zacs, zaos )
-
       ! ================================ !
       !      Namelist informations       !
@@ -1028,5 +1025,4 @@
       IF( ln_dm2dc .AND. ln_cpl ) ncpl_qsr_freq = 86400 / ncpl_qsr_freq
 
-      CALL wrk_dealloc( jpi,jpj,   zacs, zaos )
       !
       IF( nn_timing == 1 )   CALL timing_stop('sbc_cpl_init')
Index: /branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/USR/usrdef_sbc.F90
===================================================================
--- /branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/USR/usrdef_sbc.F90	(revision 8961)
+++ /branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/USR/usrdef_sbc.F90	(revision 8962)
@@ -229,6 +229,8 @@
    END SUBROUTINE usrdef_sbc_ice_tau
 
-   SUBROUTINE usrdef_sbc_ice_flx( kt )
+   SUBROUTINE usrdef_sbc_ice_flx( kt, phs, phi )
       INTEGER, INTENT(in) ::   kt   ! ocean time step
+      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   phs    ! snow thickness
+      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   phi    ! ice thickness
    END SUBROUTINE usrdef_sbc_ice_flx
 
