ecearth3/sources/nemo-3.6/CONFIG/ORCA1L75_LIM3/MY_SRC/limvar.F90

MODULE limvar
   !!======================================================================
   !!                       ***  MODULE limvar ***
   !!                 Different sets of ice model variables 
   !!                   how to switch from one to another
   !!
   !!                 There are three sets of variables
   !!                 VGLO : global variables of the model
   !!                        - v_i (jpi,jpj,jpl)
   !!                        - v_s (jpi,jpj,jpl)
   !!                        - a_i (jpi,jpj,jpl)
   !!                        - t_s (jpi,jpj,jpl)
   !!                        - e_i (jpi,jpj,nlay_i,jpl)
   !!                        - smv_i(jpi,jpj,jpl)
   !!                        - oa_i (jpi,jpj,jpl)
   !!                 VEQV : equivalent variables sometimes used in the model
   !!                        - ht_i(jpi,jpj,jpl)
   !!                        - ht_s(jpi,jpj,jpl)
   !!                        - t_i (jpi,jpj,nlay_i,jpl)
   !!                        ...
   !!                 VAGG : aggregate variables, averaged/summed over all
   !!                        thickness categories
   !!                        - vt_i(jpi,jpj)
   !!                        - vt_s(jpi,jpj)
   !!                        - at_i(jpi,jpj)
   !!                        - et_s(jpi,jpj)  !total snow heat content
   !!                        - et_i(jpi,jpj)  !total ice thermal content 
   !!                        - smt_i(jpi,jpj) !mean ice salinity
   !!                        - tm_i (jpi,jpj) !mean ice temperature
   !!======================================================================
   !! History :   -   ! 2006-01 (M. Vancoppenolle) Original code
   !!            3.4  ! 2011-02 (G. Madec) dynamical allocation
   !!----------------------------------------------------------------------
#if defined key_lim3
   !!----------------------------------------------------------------------
   !!   'key_lim3'                                      LIM3 sea-ice model
   !!----------------------------------------------------------------------
   USE par_oce        ! ocean parameters
   USE phycst         ! physical constants (ocean directory) 
   USE sbc_oce        ! Surface boundary condition: ocean fields
   USE ice            ! ice variables
   USE thd_ice        ! ice variables (thermodynamics)
   USE dom_ice        ! ice domain
   USE in_out_manager ! I/O manager
   USE lib_mpp        ! MPP library
   USE wrk_nemo       ! work arrays
   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  

   IMPLICIT NONE
   PRIVATE

   PUBLIC   lim_var_agg          
   PUBLIC   lim_var_glo2eqv      
   PUBLIC   lim_var_eqv2glo      
   PUBLIC   lim_var_salprof      
   PUBLIC   lim_var_bv           
   PUBLIC   lim_var_salprof1d    
   PUBLIC   lim_var_zapsmall
   PUBLIC   lim_var_itd

   !!----------------------------------------------------------------------
   !! NEMO/LIM3 3.5 , UCL - NEMO Consortium (2011)
   !! $Id: limvar.F90 7814 2017-03-20 16:21:42Z clem $
   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
   !!----------------------------------------------------------------------
CONTAINS

   SUBROUTINE lim_var_agg( kn )
      !!------------------------------------------------------------------
      !!                ***  ROUTINE lim_var_agg  ***
      !!
      !! ** Purpose :   aggregates ice-thickness-category variables to all-ice variables
      !!              i.e. it turns VGLO into VAGG
      !! ** Method  :
      !!
      !! ** Arguments : n = 1, at_i vt_i only
      !!                n = 2 everything
      !!
      !! note : you could add an argument when you need only at_i, vt_i
      !!        and when you need everything
      !!------------------------------------------------------------------
      INTEGER, INTENT( in ) ::   kn     ! =1 at_i & vt only ; = what is needed
      !
      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
      !!------------------------------------------------------------------

      !--------------------
      ! Compute variables
      !--------------------
      ! integrated values
      vt_i (:,:) = SUM( v_i, dim=3 )
      vt_s (:,:) = SUM( v_s, dim=3 )
      at_i (:,:) = SUM( a_i, dim=3 )
      et_s(:,:)  = SUM( SUM( e_s(:,:,:,:), dim=4 ), dim=3 )
      et_i(:,:)  = SUM( SUM( e_i(:,:,:,:), dim=4 ), dim=3 )
      !
      DO jj = 1, jpj
         DO ji = 1, jpi
            ato_i(ji,jj) = MAX( 1._wp - at_i(ji,jj), 0._wp )   ! open water fraction
         END DO
      END DO

      IF( kn > 1 ) THEN
         !
         ! mean ice/snow thickness
         DO jj = 1, jpj
            DO ji = 1, jpi
               rswitch      = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) 
               htm_i(ji,jj) = vt_i(ji,jj) / MAX( at_i(ji,jj) , epsi10 ) * rswitch
               htm_s(ji,jj) = vt_s(ji,jj) / MAX( at_i(ji,jj) , epsi10 ) * rswitch
            ENDDO
         ENDDO

         ! mean temperature (K), salinity and age
         smt_i(:,:) = 0._wp
         tm_i(:,:)  = 0._wp
         tm_su(:,:) = 0._wp
         tm_si(:,:) = 0._wp 
         om_i (:,:) = 0._wp
         DO jl = 1, jpl
            
            DO jj = 1, jpj
               DO ji = 1, jpi
                  rswitch      = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) )
                  tm_su(ji,jj) = tm_su(ji,jj) + rswitch * ( t_su(ji,jj,jl) - rt0 ) * a_i(ji,jj,jl) / MAX( at_i(ji,jj) , epsi10 )
                  tm_si(ji,jj) = tm_si(ji,jj) + rswitch * ( t_si(ji,jj,jl) - rt0 ) * a_i(ji,jj,jl) / MAX( at_i(ji,jj) , epsi10 )
                  om_i (ji,jj) = om_i (ji,jj) + rswitch *   oa_i(ji,jj,jl)                         / MAX( at_i(ji,jj) , epsi10 )
               END DO
            END DO
            
            DO jk = 1, nlay_i
               DO jj = 1, jpj
                  DO ji = 1, jpi
                     rswitch = MAX( 0._wp , SIGN( 1._wp , vt_i(ji,jj) - epsi10 ) )
                     tm_i(ji,jj)  = tm_i(ji,jj)  + r1_nlay_i * rswitch * ( t_i(ji,jj,jk,jl) - rt0 ) * v_i(ji,jj,jl)  &
                        &            / MAX( vt_i(ji,jj) , epsi10 )
                     smt_i(ji,jj) = smt_i(ji,jj) + r1_nlay_i * rswitch * s_i(ji,jj,jk,jl) * v_i(ji,jj,jl)  &
                        &            / MAX( vt_i(ji,jj) , epsi10 )
                  END DO
               END DO
            END DO
         END DO
         tm_i  = tm_i  + rt0
         tm_su = tm_su + rt0
         tm_si = tm_si + rt0
         !
      ENDIF
      !
   END SUBROUTINE lim_var_agg


   SUBROUTINE lim_var_glo2eqv
      !!------------------------------------------------------------------
      !!                ***  ROUTINE lim_var_glo2eqv ***
      !!
      !! ** Purpose :   computes equivalent variables as function of global variables 
      !!              i.e. it turns VGLO into VEQV
      !!------------------------------------------------------------------
      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
      REAL(wp) ::   zq_i, zaaa, zbbb, zccc, zdiscrim     ! local scalars
      REAL(wp) ::   ztmelts, zq_s, zfac1, zfac2   !   -      -
      !!------------------------------------------------------------------

      !-------------------------------------------------------
      ! Ice thickness, snow thickness, ice salinity, ice age
      !-------------------------------------------------------
      DO jl = 1, jpl
         DO jj = 1, jpj
            DO ji = 1, jpi
               rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi20 ) )   !0 if no ice and 1 if yes
               ht_i(ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch
            END DO
         END DO
      END DO
      ! Force the upper limit of ht_i to always be < hi_max (99 m).
      DO jj = 1, jpj
         DO ji = 1, jpi
            rswitch = MAX( 0._wp , SIGN( 1._wp, ht_i(ji,jj,jpl) - epsi20 ) )
            ht_i(ji,jj,jpl) = MIN( ht_i(ji,jj,jpl) , hi_max(jpl) )
            a_i (ji,jj,jpl) = v_i(ji,jj,jpl) / MAX( ht_i(ji,jj,jpl) , epsi20 ) * rswitch
         END DO
      END DO

      DO jl = 1, jpl
         DO jj = 1, jpj
            DO ji = 1, jpi
               rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi20 ) )   !0 if no ice and 1 if yes
               ht_s(ji,jj,jl) = v_s (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch
               o_i(ji,jj,jl)  = oa_i(ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch
            END DO
         END DO
      END DO
      
      IF(  nn_icesal == 2  )THEN
         DO jl = 1, jpl
            DO jj = 1, jpj
               DO ji = 1, jpi
                  rswitch = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) - epsi20 ) )   !0 if no ice and 1 if yes
                  sm_i(ji,jj,jl) = smv_i(ji,jj,jl) / MAX( v_i(ji,jj,jl) , epsi20 ) * rswitch
                  !                                      ! bounding salinity
                  sm_i(ji,jj,jl) = MAX( sm_i(ji,jj,jl), rn_simin )
               END DO
            END DO
         END DO
      ENDIF

      CALL lim_var_salprof      ! salinity profile

      !-------------------
      ! Ice temperatures
      !-------------------
      DO jl = 1, jpl
         DO jk = 1, nlay_i
            DO jj = 1, jpj
               DO ji = 1, jpi
                  !                                                              ! Energy of melting q(S,T) [J.m-3]
                  rswitch = MAX( 0.0 , SIGN( 1.0 , v_i(ji,jj,jl) - epsi20 ) )     ! rswitch = 0 if no ice and 1 if yes
                  zq_i    = rswitch * e_i(ji,jj,jk,jl) / MAX( v_i(ji,jj,jl) , epsi20 ) * REAL(nlay_i,wp) 
                  ztmelts = -tmut * s_i(ji,jj,jk,jl) + rt0              ! Ice layer melt temperature
                  !
                  zaaa       =  cpic                  ! Conversion q(S,T) -> T (second order equation)
                  zbbb       =  ( rcp - cpic ) * ( ztmelts - rt0 ) + zq_i * r1_rhoic - lfus
                  zccc       =  lfus * (ztmelts-rt0)
                  zdiscrim   =  SQRT( MAX(zbbb*zbbb - 4._wp*zaaa*zccc , 0._wp) )
                  t_i(ji,jj,jk,jl) = rt0 + rswitch *( - zbbb - zdiscrim ) / ( 2.0 *zaaa )
                  t_i(ji,jj,jk,jl) = MIN( ztmelts, MAX( rt0 - 100._wp, t_i(ji,jj,jk,jl) ) )  ! -100 < t_i < ztmelts
               END DO
            END DO
         END DO
      END DO

      !--------------------
      ! Snow temperatures
      !--------------------
      zfac1 = 1._wp / ( rhosn * cpic )
      zfac2 = lfus / cpic  
      DO jl = 1, jpl
         DO jk = 1, nlay_s
            DO jj = 1, jpj
               DO ji = 1, jpi
                  !Energy of melting q(S,T) [J.m-3]
                  rswitch = MAX( 0._wp , SIGN( 1._wp , v_s(ji,jj,jl) - epsi20 ) )     ! rswitch = 0 if no ice and 1 if yes
                  zq_s  = rswitch * e_s(ji,jj,jk,jl) / MAX( v_s(ji,jj,jl) , epsi20 ) * REAL(nlay_s,wp)
                  !
                  t_s(ji,jj,jk,jl) = rt0 + rswitch * ( - zfac1 * zq_s + zfac2 )
                  t_s(ji,jj,jk,jl) = MIN( rt0, MAX( rt0 - 100._wp , t_s(ji,jj,jk,jl) ) )     ! -100 < t_s < rt0
               END DO
            END DO
         END DO
      END DO

      !-------------------
      ! Mean temperature
      !-------------------
      ! integrated values
      vt_i (:,:) = SUM( v_i, dim=3 )
      vt_s (:,:) = SUM( v_s, dim=3 )
      at_i (:,:) = SUM( a_i, dim=3 )

      tm_i(:,:) = 0._wp
      DO jl = 1, jpl
         DO jk = 1, nlay_i
            DO jj = 1, jpj
               DO ji = 1, jpi
                  rswitch = MAX( 0._wp , SIGN( 1._wp , vt_i(ji,jj) - epsi10 ) )
                  tm_i(ji,jj) = tm_i(ji,jj) + r1_nlay_i * rswitch * ( t_i(ji,jj,jk,jl) - rt0 ) * v_i(ji,jj,jl)  &
                     &            / MAX( vt_i(ji,jj) , epsi10 )
               END DO
            END DO
         END DO
      END DO
      tm_i = tm_i + rt0
      !
   END SUBROUTINE lim_var_glo2eqv


   SUBROUTINE lim_var_eqv2glo
      !!------------------------------------------------------------------
      !!                ***  ROUTINE lim_var_eqv2glo ***
      !!
      !! ** Purpose :   computes global variables as function of equivalent variables
      !!                i.e. it turns VEQV into VGLO
      !! ** Method  :
      !!
      !! ** History :  (01-2006) Martin Vancoppenolle, UCL-ASTR
      !!------------------------------------------------------------------
      !
      v_i(:,:,:)   = ht_i(:,:,:) * a_i(:,:,:)
      v_s(:,:,:)   = ht_s(:,:,:) * a_i(:,:,:)
      smv_i(:,:,:) = sm_i(:,:,:) * v_i(:,:,:)
      !
   END SUBROUTINE lim_var_eqv2glo


   SUBROUTINE lim_var_salprof
      !!------------------------------------------------------------------
      !!                ***  ROUTINE lim_var_salprof ***
      !!
      !! ** Purpose :   computes salinity profile in function of bulk salinity     
      !!
      !! ** Method  : If bulk salinity greater than zsi1, 
      !!              the profile is assumed to be constant (S_inf)
      !!              If bulk salinity lower than zsi0,
      !!              the profile is linear with 0 at the surface (S_zero)
      !!              If it is between zsi0 and zsi1, it is a
      !!              alpha-weighted linear combination of s_inf and s_zero
      !!
      !! ** References : Vancoppenolle et al., 2007
      !!------------------------------------------------------------------
      INTEGER  ::   ji, jj, jk, jl   ! dummy loop index
      REAL(wp) ::   zfac0, zfac1, zsal
      REAL(wp) ::   zswi0, zswi01, zargtemp , zs_zero   
      REAL(wp), POINTER, DIMENSION(:,:,:) ::   z_slope_s, zalpha
      REAL(wp), PARAMETER :: zsi0 = 3.5_wp
      REAL(wp), PARAMETER :: zsi1 = 4.5_wp
      !!------------------------------------------------------------------

      CALL wrk_alloc( jpi, jpj, jpl, z_slope_s, zalpha )

      !---------------------------------------
      ! Vertically constant, constant in time
      !---------------------------------------
      IF(  nn_icesal == 1  )  THEN
         s_i (:,:,:,:) = rn_icesal
         sm_i(:,:,:)   = rn_icesal
      ENDIF

      !-----------------------------------
      ! Salinity profile, varying in time
      !-----------------------------------
      IF(  nn_icesal == 2  ) THEN
         !
         DO jk = 1, nlay_i
            s_i(:,:,jk,:)  = sm_i(:,:,:)
         END DO
         !
         DO jl = 1, jpl                               ! Slope of the linear profile 
            DO jj = 1, jpj
               DO ji = 1, jpi
                  rswitch = MAX( 0._wp , SIGN( 1._wp , ht_i(ji,jj,jl) - epsi20 ) )
                  z_slope_s(ji,jj,jl) = rswitch * 2._wp * sm_i(ji,jj,jl) / MAX( epsi20 , ht_i(ji,jj,jl) )
               END DO
            END DO
         END DO
         !
         zfac0 = 1._wp / ( zsi0 - zsi1 )       ! Weighting factor between zs_zero and zs_inf
         zfac1 = zsi1  / ( zsi1 - zsi0 )
         !
         zalpha(:,:,:) = 0._wp
         DO jl = 1, jpl
            DO jj = 1, jpj
               DO ji = 1, jpi
                  ! zswi0 = 1 if sm_i le zsi0 and 0 otherwise
                  zswi0  = MAX( 0._wp   , SIGN( 1._wp  , zsi0 - sm_i(ji,jj,jl) ) ) 
                  ! zswi01 = 1 if sm_i is between zsi0 and zsi1 and 0 othws 
                  zswi01 = ( 1._wp - zswi0 ) * MAX( 0._wp   , SIGN( 1._wp  , zsi1 - sm_i(ji,jj,jl) ) ) 
                  ! If 2.sm_i GE sss_m then rswitch = 1
                  ! this is to force a constant salinity profile in the Baltic Sea
                  rswitch = MAX( 0._wp , SIGN( 1._wp , 2._wp * sm_i(ji,jj,jl) - sss_m(ji,jj) ) )
                  zalpha(ji,jj,jl) = zswi0  + zswi01 * ( sm_i(ji,jj,jl) * zfac0 + zfac1 )
                  zalpha(ji,jj,jl) = zalpha(ji,jj,jl) * ( 1._wp - rswitch )
               END DO
            END DO
         END DO

         ! Computation of the profile
         DO jl = 1, jpl
            DO jk = 1, nlay_i
               DO jj = 1, jpj
                  DO ji = 1, jpi
                     !                                      ! linear profile with 0 at the surface
                     zs_zero = z_slope_s(ji,jj,jl) * ( REAL(jk,wp) - 0.5_wp ) * ht_i(ji,jj,jl) * r1_nlay_i
                     !                                      ! weighting the profile
                     s_i(ji,jj,jk,jl) = zalpha(ji,jj,jl) * zs_zero + ( 1._wp - zalpha(ji,jj,jl) ) * sm_i(ji,jj,jl)
                     !                                      ! bounding salinity
                     s_i(ji,jj,jk,jl) = MIN( rn_simax, MAX( s_i(ji,jj,jk,jl), rn_simin ) )
                  END DO
               END DO
            END DO
         END DO
         !
      ENDIF ! nn_icesal

      !-------------------------------------------------------
      ! Vertically varying salinity profile, constant in time
      !-------------------------------------------------------

      IF(  nn_icesal == 3  ) THEN      ! Schwarzacher (1959) multiyear salinity profile (mean = 2.30)
         !
         sm_i(:,:,:) = 2.30_wp
         !
         DO jl = 1, jpl
            DO jk = 1, nlay_i
               zargtemp  = ( REAL(jk,wp) - 0.5_wp ) * r1_nlay_i
               zsal =  1.6_wp * (  1._wp - COS( rpi * zargtemp**(0.407_wp/(0.573_wp+zargtemp)) )  )
               s_i(:,:,jk,jl) =  zsal
            END DO
         END DO
         !
      ENDIF ! nn_icesal
      !
      CALL wrk_dealloc( jpi, jpj, jpl, z_slope_s, zalpha )
      !
   END SUBROUTINE lim_var_salprof


   SUBROUTINE lim_var_bv
      !!------------------------------------------------------------------
      !!                ***  ROUTINE lim_var_bv ***
      !!
      !! ** Purpose :   computes mean brine volume (%) in sea ice
      !!
      !! ** Method  : e = - 0.054 * S (ppt) / T (C)
      !!
      !! References : Vancoppenolle et al., JGR, 2007
      !!------------------------------------------------------------------
      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
      !!------------------------------------------------------------------
      !
      bvm_i(:,:)   = 0._wp
      bv_i (:,:,:) = 0._wp
      DO jl = 1, jpl
         DO jk = 1, nlay_i
            DO jj = 1, jpj
               DO ji = 1, jpi
                  rswitch        = ( 1._wp - MAX( 0._wp , SIGN( 1._wp , (t_i(ji,jj,jk,jl) - rt0) + epsi10 ) )  )
                  bv_i(ji,jj,jl) = bv_i(ji,jj,jl) - rswitch * tmut * s_i(ji,jj,jk,jl) * r1_nlay_i  &
                     &                            / MIN( t_i(ji,jj,jk,jl) - rt0, - epsi10 )
               END DO
            END DO
         END DO
         
         DO jj = 1, jpj
            DO ji = 1, jpi
               rswitch      = MAX( 0._wp , SIGN( 1._wp , vt_i(ji,jj) - epsi10 ) )
               bvm_i(ji,jj) = bvm_i(ji,jj) + rswitch * bv_i(ji,jj,jl) * v_i(ji,jj,jl) / MAX( vt_i(ji,jj), epsi10 )
            END DO
         END DO
      END DO
      !
   END SUBROUTINE lim_var_bv


   SUBROUTINE lim_var_salprof1d( kideb, kiut )
      !!-------------------------------------------------------------------
      !!                  ***  ROUTINE lim_thd_salprof1d  ***
      !!
      !! ** Purpose :   1d computation of the sea ice salinity profile
      !!                Works with 1d vectors and is used by thermodynamic modules
      !!-------------------------------------------------------------------
      INTEGER, INTENT(in) ::   kideb, kiut   ! thickness category index
      !
      INTEGER  ::   ji, jk    ! dummy loop indices
      INTEGER  ::   ii, ij    ! local integers
      REAL(wp) ::   zfac0, zfac1, zargtemp, zsal   ! local scalars
      REAL(wp) ::   zalpha, zswi0, zswi01, zs_zero              !   -      -
      !
      REAL(wp), POINTER, DIMENSION(:) ::   z_slope_s
      REAL(wp), PARAMETER :: zsi0 = 3.5_wp
      REAL(wp), PARAMETER :: zsi1 = 4.5_wp
      !!---------------------------------------------------------------------

      CALL wrk_alloc( jpij, z_slope_s )

      !---------------------------------------
      ! Vertically constant, constant in time
      !---------------------------------------
      IF( nn_icesal == 1 )   s_i_1d(:,:) = rn_icesal

      !------------------------------------------------------
      ! Vertically varying salinity profile, varying in time
      !------------------------------------------------------

      IF(  nn_icesal == 2  ) THEN
         !
         DO ji = kideb, kiut          ! Slope of the linear profile zs_zero
            rswitch = MAX( 0._wp , SIGN( 1._wp , ht_i_1d(ji) - epsi20 ) )
            z_slope_s(ji) = rswitch * 2._wp * sm_i_1d(ji) / MAX( epsi20 , ht_i_1d(ji) )
         END DO

         ! Weighting factor between zs_zero and zs_inf
         !---------------------------------------------
         zfac0 = 1._wp / ( zsi0 - zsi1 )
         zfac1 = zsi1 / ( zsi1 - zsi0 )
         DO jk = 1, nlay_i
            DO ji = kideb, kiut
               ii =  MOD( npb(ji) - 1 , jpi ) + 1
               ij =     ( npb(ji) - 1 ) / jpi + 1
               ! zswi0 = 1 if sm_i le zsi0 and 0 otherwise
               zswi0  = MAX( 0._wp , SIGN( 1._wp  , zsi0 - sm_i_1d(ji) ) ) 
               ! zswi01 = 1 if sm_i is between zsi0 and zsi1 and 0 othws 
               zswi01 = ( 1._wp - zswi0 ) * MAX( 0._wp , SIGN( 1._wp , zsi1 - sm_i_1d(ji) ) ) 
               ! if 2.sm_i GE sss_m then rswitch = 1
               ! this is to force a constant salinity profile in the Baltic Sea
               rswitch = MAX( 0._wp , SIGN( 1._wp , 2._wp * sm_i_1d(ji) - sss_m(ii,ij) ) )
               !
               zalpha = (  zswi0 + zswi01 * ( sm_i_1d(ji) * zfac0 + zfac1 )  ) * ( 1._wp - rswitch )
               !
               zs_zero = z_slope_s(ji) * ( REAL(jk,wp) - 0.5_wp ) * ht_i_1d(ji) * r1_nlay_i
               ! weighting the profile
               s_i_1d(ji,jk) = zalpha * zs_zero + ( 1._wp - zalpha ) * sm_i_1d(ji)
               ! bounding salinity
               s_i_1d(ji,jk) = MIN( rn_simax, MAX( s_i_1d(ji,jk), rn_simin ) )
            END DO 
         END DO 

      ENDIF 

      !-------------------------------------------------------
      ! Vertically varying salinity profile, constant in time
      !-------------------------------------------------------

      IF( nn_icesal == 3 ) THEN      ! Schwarzacher (1959) multiyear salinity profile (mean = 2.30)
         !
         sm_i_1d(:) = 2.30_wp
         !
         DO jk = 1, nlay_i
            zargtemp  = ( REAL(jk,wp) - 0.5_wp ) * r1_nlay_i
            zsal =  1.6_wp * ( 1._wp - COS( rpi * zargtemp**( 0.407_wp / ( 0.573_wp + zargtemp ) ) ) )
            DO ji = kideb, kiut
               s_i_1d(ji,jk) = zsal
            END DO
         END DO
         !
      ENDIF
      !
      CALL wrk_dealloc( jpij, z_slope_s )
      !
   END SUBROUTINE lim_var_salprof1d

   SUBROUTINE lim_var_zapsmall
      !!-------------------------------------------------------------------
      !!                   ***  ROUTINE lim_var_zapsmall ***
      !!
      !! ** Purpose :   Remove too small sea ice areas and correct fluxes
      !!
      !! history : LIM3.5 - 01-2014 (C. Rousset) original code
      !!-------------------------------------------------------------------
      INTEGER  ::   ji, jj, jl, jk   ! dummy loop indices
      REAL(wp) ::   zsal, zvi, zvs, zei, zes
      !!-------------------------------------------------------------------
      at_i (:,:) = 0._wp
      DO jl = 1, jpl
         at_i(:,:) = at_i(:,:) + a_i(:,:,jl)
      END DO

      DO jl = 1, jpl

         !-----------------------------------------------------------------
         ! Zap ice energy and use ocean heat to melt ice
         !-----------------------------------------------------------------
         DO jk = 1, nlay_i
            DO jj = 1 , jpj
               DO ji = 1 , jpi
                  rswitch          = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi10 ) )
                  rswitch          = MAX( 0._wp , SIGN( 1._wp, at_i(ji,jj  ) - epsi10 ) ) * rswitch
                  rswitch          = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) - epsi10 ) ) * rswitch
                  rswitch          = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) * rswitch  &
                     &                                       / MAX( a_i(ji,jj,jl), epsi10 ) - epsi10 ) ) * rswitch
                  zei              = e_i(ji,jj,jk,jl)
                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * rswitch
                  t_i(ji,jj,jk,jl) = t_i(ji,jj,jk,jl) * rswitch + rt0 * ( 1._wp - rswitch )
                  ! update exchanges with ocean
                  hfx_res(ji,jj)   = hfx_res(ji,jj) + ( e_i(ji,jj,jk,jl) - zei ) * r1_rdtice ! W.m-2 <0
               END DO
            END DO
         END DO

         DO jj = 1 , jpj
            DO ji = 1 , jpi
               rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi10 ) )
               rswitch = MAX( 0._wp , SIGN( 1._wp, at_i(ji,jj  ) - epsi10 ) ) * rswitch
               rswitch = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) - epsi10 ) ) * rswitch
               rswitch = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) * rswitch  &
                  &                              / MAX( a_i(ji,jj,jl), epsi10 ) - epsi10 ) ) * rswitch
               zsal = smv_i(ji,jj,  jl)
               zvi  = v_i  (ji,jj,  jl)
               zvs  = v_s  (ji,jj,  jl)
               zes  = e_s  (ji,jj,1,jl)
               !-----------------------------------------------------------------
               ! Zap snow energy 
               !-----------------------------------------------------------------
               t_s(ji,jj,1,jl) = t_s(ji,jj,1,jl) * rswitch + rt0 * ( 1._wp - rswitch )
               e_s(ji,jj,1,jl) = e_s(ji,jj,1,jl) * rswitch

               !-----------------------------------------------------------------
               ! zap ice and snow volume, add water and salt to ocean
               !-----------------------------------------------------------------
               ato_i(ji,jj)    = a_i  (ji,jj,jl) * ( 1._wp - rswitch ) + ato_i(ji,jj)
               a_i  (ji,jj,jl) = a_i  (ji,jj,jl) * rswitch
               v_i  (ji,jj,jl) = v_i  (ji,jj,jl) * rswitch
               v_s  (ji,jj,jl) = v_s  (ji,jj,jl) * rswitch
               t_su (ji,jj,jl) = t_su (ji,jj,jl) * rswitch + t_bo(ji,jj) * ( 1._wp - rswitch )
               oa_i (ji,jj,jl) = oa_i (ji,jj,jl) * rswitch
               smv_i(ji,jj,jl) = smv_i(ji,jj,jl) * rswitch

               ! update exchanges with ocean
               sfx_res(ji,jj)  = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsal ) * rhoic * r1_rdtice
               wfx_res(ji,jj)  = wfx_res(ji,jj) - ( v_i(ji,jj,jl)   - zvi  ) * rhoic * r1_rdtice
               wfx_snw(ji,jj)  = wfx_snw(ji,jj) - ( v_s(ji,jj,jl)   - zvs  ) * rhosn * r1_rdtice
               hfx_res(ji,jj)  = hfx_res(ji,jj) + ( e_s(ji,jj,1,jl) - zes ) * r1_rdtice ! W.m-2 <0
            END DO
         END DO
      END DO 

      ! to be sure that at_i is the sum of a_i(jl)
      at_i (:,:) = 0._wp
      DO jl = 1, jpl
         at_i(:,:) = at_i(:,:) + a_i(:,:,jl)
      END DO

      ! open water = 1 if at_i=0
      DO jj = 1, jpj
         DO ji = 1, jpi
            rswitch      = MAX( 0._wp , SIGN( 1._wp, - at_i(ji,jj) ) )
            ato_i(ji,jj) = rswitch + (1._wp - rswitch ) * ato_i(ji,jj)
         END DO
      END DO

      !
   END SUBROUTINE lim_var_zapsmall

   SUBROUTINE lim_var_itd( zhti, zhts, zai, zht_i, zht_s, za_i )
      !!------------------------------------------------------------------
      !!                ***  ROUTINE lim_var_itd   ***
      !!
      !! ** Purpose :  converting 1-cat ice to multiple ice categories
      !!
      !!                  ice thickness distribution follows a gaussian law
      !!               around the concentration of the most likely ice thickness
      !!                           (similar as limistate.F90)
      !!
      !! ** Method:   Iterative procedure
      !!                
      !!               1) Try to fill the jpl ice categories (bounds hi_max(0:jpl)) with a gaussian
      !!
      !!               2) Check whether the distribution conserves area and volume, positivity and
      !!                  category boundaries
      !!              
      !!               3) If not (input ice is too thin), the last category is empty and
      !!                  the number of categories is reduced (jpl-1)
      !!
      !!               4) Iterate until ok (SUM(itest(:) = 4)
      !!
      !! ** Arguments : zhti: 1-cat ice thickness
      !!                zhts: 1-cat snow depth
      !!                zai : 1-cat ice concentration
      !!
      !! ** Output    : jpl-cat 
      !!
      !!  (Example of application: BDY forcings when input are cell averaged)  
      !!
      !!-------------------------------------------------------------------
      !! History : LIM3.5 - 2012    (M. Vancoppenolle)  Original code
      !!                    2014    (C. Rousset)        Rewriting
      !!-------------------------------------------------------------------
      !! Local variables
      INTEGER  :: ji, jk, jl             ! dummy loop indices
      INTEGER  :: ijpij, i_fill, jl0  
      REAL(wp) :: zarg, zV, zconv, zdh, zdv
      REAL(wp), DIMENSION(:),   INTENT(in)    ::   zhti, zhts, zai    ! input ice/snow variables
      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   zht_i, zht_s, za_i ! output ice/snow variables
      INTEGER , POINTER, DIMENSION(:)         ::   itest
 
      CALL wrk_alloc( 4, itest )
      !--------------------------------------------------------------------
      ! initialisation of variables
      !--------------------------------------------------------------------
      ijpij = SIZE(zhti,1)
      zht_i(1:ijpij,1:jpl) = 0._wp
      zht_s(1:ijpij,1:jpl) = 0._wp
      za_i (1:ijpij,1:jpl) = 0._wp

      ! ----------------------------------------
      ! distribution over the jpl ice categories
      ! ----------------------------------------
      DO ji = 1, ijpij
         
         IF( zhti(ji) > 0._wp ) THEN

            ! find which category (jl0) the input ice thickness falls into
            jl0 = jpl
            DO jl = 1, jpl
               IF ( ( zhti(ji) >= hi_max(jl-1) ) .AND. ( zhti(ji) < hi_max(jl) ) ) THEN
                  jl0 = jl
                  CYCLE
               ENDIF
            END DO

            ! initialisation of tests
            itest(:)  = 0
         
            i_fill = jpl + 1                                             !====================================
            DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) )  ! iterative loop on i_fill categories
               ! iteration                                               !====================================
               i_fill = i_fill - 1
               
               ! initialisation of ice variables for each try
               zht_i(ji,1:jpl) = 0._wp
               za_i (ji,1:jpl) = 0._wp
               itest(:)        = 0      
               
               ! *** case very thin ice: fill only category 1
               IF ( i_fill == 1 ) THEN
                  zht_i(ji,1) = zhti(ji)
                  za_i (ji,1) = zai (ji)
                  
               ! *** case ice is thicker: fill categories >1
               ELSE

                  ! Fill ice thicknesses in the (i_fill-1) cat by hmean 
                  DO jl = 1, i_fill - 1
                     zht_i(ji,jl) = hi_mean(jl)
                  END DO
                  
                  ! Concentrations in the (i_fill-1) categories 
                  za_i(ji,jl0) = zai(ji) / SQRT(REAL(jpl))
                  DO jl = 1, i_fill - 1
                     IF ( jl /= jl0 ) THEN
                        zarg        = ( zht_i(ji,jl) - zhti(ji) ) / ( zhti(ji) * 0.5_wp )
                        za_i(ji,jl) =   za_i (ji,jl0) * EXP(-zarg**2)
                     ENDIF
                  END DO
                  
                  ! Concentration in the last (i_fill) category
                  za_i(ji,i_fill) = zai(ji) - SUM( za_i(ji,1:i_fill-1) )
                  
                  ! Ice thickness in the last (i_fill) category
                  zV = SUM( za_i(ji,1:i_fill-1) * zht_i(ji,1:i_fill-1) )
                  zht_i(ji,i_fill) = ( zhti(ji) * zai(ji) - zV ) / MAX( za_i(ji,i_fill), epsi10 ) 
                  
                  ! clem: correction if concentration of upper cat is greater than lower cat
                  !       (it should be a gaussian around jl0 but sometimes it is not)
                  IF ( jl0 /= jpl ) THEN
                     DO jl = jpl, jl0+1, -1
                        IF ( za_i(ji,jl) > za_i(ji,jl-1) ) THEN
                           zdv = zht_i(ji,jl) * za_i(ji,jl)
                           zht_i(ji,jl    ) = 0._wp
                           za_i (ji,jl    ) = 0._wp
                           za_i (ji,1:jl-1) = za_i(ji,1:jl-1) + zdv / MAX( REAL(jl-1) * zhti(ji), epsi10 )
                        END IF
                     ENDDO
                  ENDIF
               
               ENDIF ! case ice is thick or thin
            
               !---------------------
               ! Compatibility tests
               !--------------------- 
               ! Test 1: area conservation
               zconv = ABS( zai(ji) - SUM( za_i(ji,1:jpl) ) )
               IF ( zconv < epsi06 ) itest(1) = 1
            
               ! Test 2: volume conservation
               zconv = ABS( zhti(ji)*zai(ji) - SUM( za_i(ji,1:jpl)*zht_i(ji,1:jpl) ) )
               IF ( zconv < epsi06 ) itest(2) = 1
               
               ! Test 3: thickness of the last category is in-bounds ?
               IF ( zht_i(ji,i_fill) >= hi_max(i_fill-1) ) itest(3) = 1
               
               ! Test 4: positivity of ice concentrations
               itest(4) = 1
               DO jl = 1, i_fill
                  IF ( za_i(ji,jl) < 0._wp ) itest(4) = 0
               END DO
               !                                         !============================
            END DO                                       ! end iteration on categories
               !                                         !============================
         ENDIF ! if zhti > 0
      END DO ! i loop

      ! ------------------------------------------------
      ! Adding Snow in each category where za_i is not 0
      ! ------------------------------------------------ 
      DO jl = 1, jpl
         DO ji = 1, ijpij
            IF( za_i(ji,jl) > 0._wp ) THEN
               zht_s(ji,jl) = zht_i(ji,jl) * ( zhts(ji) / zhti(ji) )
               ! In case snow load is in excess that would lead to transformation from snow to ice
               ! Then, transfer the snow excess into the ice (different from limthd_dh)
               zdh = MAX( 0._wp, ( rhosn * zht_s(ji,jl) + ( rhoic - rau0 ) * zht_i(ji,jl) ) * r1_rau0 ) 
               ! recompute ht_i, ht_s avoiding out of bounds values
               zht_i(ji,jl) = MIN( hi_max(jl), zht_i(ji,jl) + zdh )
               zht_s(ji,jl) = MAX( 0._wp, zht_s(ji,jl) - zdh * rhoic * r1_rhosn )
            ENDIF
         ENDDO
      ENDDO

      CALL wrk_dealloc( 4, itest )
      !
    END SUBROUTINE lim_var_itd


#else
   !!----------------------------------------------------------------------
   !!   Default option         Dummy module          NO  LIM3 sea-ice model
   !!----------------------------------------------------------------------
CONTAINS
   SUBROUTINE lim_var_agg          ! Empty routines
   END SUBROUTINE lim_var_agg
   SUBROUTINE lim_var_glo2eqv      ! Empty routines
   END SUBROUTINE lim_var_glo2eqv
   SUBROUTINE lim_var_eqv2glo      ! Empty routines
   END SUBROUTINE lim_var_eqv2glo
   SUBROUTINE lim_var_salprof      ! Empty routines
   END SUBROUTINE lim_var_salprof
   SUBROUTINE lim_var_bv           ! Emtpy routines
   END SUBROUTINE lim_var_bv
   SUBROUTINE lim_var_salprof1d    ! Emtpy routines
   END SUBROUTINE lim_var_salprof1d
   SUBROUTINE lim_var_zapsmall
   END SUBROUTINE lim_var_zapsmall
   SUBROUTINE lim_var_itd
   END SUBROUTINE lim_var_itd
#endif

   !!======================================================================
END MODULE limvar