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

MODULE limistate
   !!======================================================================
   !!                     ***  MODULE  limistate  ***
   !!              Initialisation of diagnostics ice variables
   !!======================================================================
   !! History :  2.0  ! 2004-01 (C. Ethe, G. Madec)  Original code
   !!            4.0  ! 2011-02 (G. Madec) dynamical allocation
   !!             -   ! 2014    (C. Rousset) add N/S initializations
   !!----------------------------------------------------------------------
#if defined key_lim3
   !!----------------------------------------------------------------------
   !!   'key_lim3' :                                    LIM3 sea-ice model
   !!----------------------------------------------------------------------
   !!   lim_istate      :  Initialisation of diagnostics ice variables
   !!   lim_istate_init :  initialization of ice state and namelist read
   !!----------------------------------------------------------------------
   USE phycst           ! physical constant
   USE oce              ! dynamics and tracers variables
   USE dom_oce          ! ocean domain
   USE sbc_oce          ! Surface boundary condition: ocean fields
   USE sbc_ice          ! Surface boundary condition: ice fields
   USE eosbn2           ! equation of state
   USE ice              ! sea-ice variables
   USE par_oce          ! ocean parameters
   USE dom_ice          ! sea-ice domain
   USE limvar           ! lim_var_salprof
   USE in_out_manager   ! I/O manager
   USE lib_mpp          ! MPP library
   USE lib_fortran      ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
   USE wrk_nemo         ! work arrays

   IMPLICIT NONE
   PRIVATE

   PUBLIC   lim_istate      ! routine called by lim_init.F90

   !                          !!** init namelist (namiceini) **
   REAL(wp) ::   rn_thres_sst   ! threshold water temperature for initial sea ice
   REAL(wp) ::   rn_hts_ini_n   ! initial snow thickness in the north
   REAL(wp) ::   rn_hts_ini_s   ! initial snow thickness in the south
   REAL(wp) ::   rn_hti_ini_n   ! initial ice thickness in the north
   REAL(wp) ::   rn_hti_ini_s   ! initial ice thickness in the south
   REAL(wp) ::   rn_ati_ini_n   ! initial leads area in the north
   REAL(wp) ::   rn_ati_ini_s   ! initial leads area in the south
   REAL(wp) ::   rn_smi_ini_n   ! initial salinity 
   REAL(wp) ::   rn_smi_ini_s   ! initial salinity
   REAL(wp) ::   rn_tmi_ini_n   ! initial temperature
   REAL(wp) ::   rn_tmi_ini_s   ! initial temperature

   LOGICAL  ::  ln_iceini    ! initialization or not
   !!----------------------------------------------------------------------
   !!   LIM 3.0,  UCL-LOCEAN-IPSL (2008)
   !! $Id: limistate.F90 8156 2017-06-08 16:13:32Z vancop $
   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
   !!----------------------------------------------------------------------
CONTAINS

   SUBROUTINE lim_istate
      !!-------------------------------------------------------------------
      !!                    ***  ROUTINE lim_istate  ***
      !!
      !! ** Purpose :   defined the sea-ice initial state
      !!
      !! ** Method  :   
      !!                This routine will put some ice where ocean
      !!                is at the freezing point, then fill in ice 
      !!                state variables using prescribed initial 
      !!                values in the namelist            
      !!
      !! ** Steps   :   
      !!                1) Read namelist
      !!                2) Basal temperature; ice and hemisphere masks
      !!                3) Fill in the ice thickness distribution using gaussian
      !!                4) Fill in space-dependent arrays for state variables
      !!                5) Diagnostic arrays
      !!                6) Lateral boundary conditions
      !!
      !! ** Notes   : o_i, t_su, t_s, t_i, s_i must be filled everywhere, even
      !!              where there is no ice (clem: I do not know why, is it mandatory?) 
      !!
      !! History :
      !!   2.0  !  01-04  (C. Ethe, G. Madec)  Original code
      !!   3.0  !  2007   (M. Vancoppenolle)   Rewrite for ice cats
      !!   4.0  !  09-11  (M. Vancoppenolle)   Enhanced version for ice cats
      !!--------------------------------------------------------------------

      !! * Local variables
      INTEGER    :: ji, jj, jk, jl             ! dummy loop indices
      REAL(wp)   :: ztmelts, zdh
      INTEGER    :: i_hemis, i_fill, jl0  
      REAL(wp)   :: ztest_1, ztest_2, ztest_3, ztest_4, ztests, zsigma, zarg, zA, zV, zA_cons, zV_cons, zconv
      REAL(wp), POINTER, DIMENSION(:)     :: zht_i_ini, zat_i_ini, zvt_i_ini, zht_s_ini, zsm_i_ini, ztm_i_ini
      REAL(wp), POINTER, DIMENSION(:,:)   :: zh_i_ini, za_i_ini, zv_i_ini
      REAL(wp), POINTER, DIMENSION(:,:)   :: zswitch    ! ice indicator
      INTEGER,  POINTER, DIMENSION(:,:)   :: zhemis   ! hemispheric index
      !--------------------------------------------------------------------

      CALL wrk_alloc( jpi, jpj, zswitch )
      CALL wrk_alloc( jpi, jpj, zhemis )
      CALL wrk_alloc( jpl,   2, zh_i_ini,  za_i_ini,  zv_i_ini )
      CALL wrk_alloc(   2,      zht_i_ini, zat_i_ini, zvt_i_ini, zht_s_ini, zsm_i_ini, ztm_i_ini )

      IF(lwp) WRITE(numout,*)
      IF(lwp) WRITE(numout,*) 'lim_istate : Ice initialization '
      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ '

      !--------------------------------------------------------------------
      ! 1) Read namelist
      !--------------------------------------------------------------------

      CALL lim_istate_init     !  reading the initials parameters of the ice

      ! surface temperature
      DO jl = 1, jpl ! loop over categories
         t_su  (:,:,jl) = rt0 * tmask(:,:,1)
         tn_ice(:,:,jl) = rt0 * tmask(:,:,1)
      END DO

      ! basal temperature (considered at freezing point)
      CALL eos_fzp( sss_m(:,:), t_bo(:,:) )
      t_bo(:,:) = ( t_bo(:,:) + rt0 ) * tmask(:,:,1) 

      IF( ln_iceini ) THEN

      !--------------------------------------------------------------------
      ! 2) Basal temperature, ice mask and hemispheric index
      !--------------------------------------------------------------------

      DO jj = 1, jpj                                       ! ice if sst <= t-freez + ttest
         DO ji = 1, jpi
            IF( ( sst_m(ji,jj)  - ( t_bo(ji,jj) - rt0 ) ) * tmask(ji,jj,1) >= rn_thres_sst ) THEN 
               zswitch(ji,jj) = 0._wp * tmask(ji,jj,1)    ! no ice
            ELSE                                                                                   
               zswitch(ji,jj) = 1._wp * tmask(ji,jj,1)    !    ice
            ENDIF
         END DO
      END DO


      ! Hemispheric index
      DO jj = 1, jpj
         DO ji = 1, jpi
            IF( fcor(ji,jj) >= 0._wp ) THEN    
               zhemis(ji,jj) = 1 ! Northern hemisphere
            ELSE
               zhemis(ji,jj) = 2 ! Southern hemisphere
            ENDIF
         END DO
      END DO

      !--------------------------------------------------------------------
      ! 3) Initialization of sea ice state variables
      !--------------------------------------------------------------------

      !-----------------------------
      ! 3.1) Hemisphere-dependent arrays
      !-----------------------------
      ! assign initial thickness, concentration, snow depth and salinity to an hemisphere-dependent array
      zht_i_ini(1) = rn_hti_ini_n ; zht_i_ini(2) = rn_hti_ini_s  ! ice thickness
      zht_s_ini(1) = rn_hts_ini_n ; zht_s_ini(2) = rn_hts_ini_s  ! snow depth
      zat_i_ini(1) = rn_ati_ini_n ; zat_i_ini(2) = rn_ati_ini_s  ! ice concentration
      zsm_i_ini(1) = rn_smi_ini_n ; zsm_i_ini(2) = rn_smi_ini_s  ! bulk ice salinity
      ztm_i_ini(1) = rn_tmi_ini_n ; ztm_i_ini(2) = rn_tmi_ini_s  ! temperature (ice and snow)

      zvt_i_ini(:) = zht_i_ini(:) * zat_i_ini(:)   ! ice volume

      !---------------------------------------------------------------------
      ! 3.2) Distribute ice concentration and thickness into the categories
      !---------------------------------------------------------------------
      ! a gaussian distribution for ice concentration is used
      ! then we check whether the distribution fullfills
      ! volume and area conservation, positivity and ice categories bounds
      DO i_hemis = 1, 2

      ztest_1 = 0 ; ztest_2 = 0 ; ztest_3 = 0 ; ztest_4 = 0

      ! note for the great nemo engineers: 
      ! only very few of the WRITE statements are necessary for the reference version
      ! they were one day useful, but now i personally doubt of their
      ! potential for bringing anything useful

      DO i_fill = jpl, 1, -1
         IF ( ( ztest_1 + ztest_2 + ztest_3 + ztest_4 ) .NE. 4 ) THEN
            !----------------------------
            ! fill the i_fill categories
            !----------------------------
            ! *** 1 category to fill
            IF ( i_fill .EQ. 1 ) THEN
               zh_i_ini(1,i_hemis)       = zht_i_ini(i_hemis)
               za_i_ini(1,i_hemis)       = zat_i_ini(i_hemis)
               zh_i_ini(2:jpl,i_hemis)   = 0._wp
               za_i_ini(2:jpl,i_hemis)   = 0._wp
            ELSE

               ! *** >1 categores to fill
               !--- Ice thicknesses in the i_fill - 1 first categories
               DO jl = 1, i_fill - 1
                  zh_i_ini(jl,i_hemis) = hi_mean(jl)
               END DO
               
               !--- jl0: most likely index where cc will be maximum
               DO jl = 1, jpl
                  IF ( ( zht_i_ini(i_hemis) >  hi_max(jl-1) ) .AND. &
                     & ( zht_i_ini(i_hemis) <= hi_max(jl)   ) ) THEN
                     jl0 = jl
                  ENDIF
               END DO
               jl0 = MIN(jl0, i_fill)
               
               !--- Concentrations
               za_i_ini(jl0,i_hemis)      = zat_i_ini(i_hemis) / SQRT(REAL(jpl))
               DO jl = 1, i_fill - 1
                  IF ( jl .NE. jl0 ) THEN
                     zsigma               = 0.5 * zht_i_ini(i_hemis)
                     zarg                 = ( zh_i_ini(jl,i_hemis) - zht_i_ini(i_hemis) ) / zsigma
                     za_i_ini(jl,i_hemis) = za_i_ini(jl0,i_hemis) * EXP(-zarg**2)
                  ENDIF
               END DO
               
               zA = 0. ! sum of the areas in the jpl categories 
               DO jl = 1, i_fill - 1
                 zA = zA + za_i_ini(jl,i_hemis)
               END DO
               za_i_ini(i_fill,i_hemis)   = zat_i_ini(i_hemis) - zA ! ice conc in the last category
               IF ( i_fill .LT. jpl ) za_i_ini(i_fill+1:jpl, i_hemis) = 0._wp
         
               !--- Ice thickness in the last category
               zV = 0. ! sum of the volumes of the N-1 categories
               DO jl = 1, i_fill - 1
                  zV = zV + za_i_ini(jl,i_hemis)*zh_i_ini(jl,i_hemis)
               END DO
               zh_i_ini(i_fill,i_hemis) = ( zvt_i_ini(i_hemis) - zV ) / za_i_ini(i_fill,i_hemis) 
               IF ( i_fill .LT. jpl ) zh_i_ini(i_fill+1:jpl, i_hemis) = 0._wp

               !--- volumes
               zv_i_ini(:,i_hemis) = za_i_ini(:,i_hemis) * zh_i_ini(:,i_hemis)
               IF ( i_fill .LT. jpl ) zv_i_ini(i_fill+1:jpl, i_hemis) = 0._wp

            ENDIF ! i_fill

            !---------------------
            ! Compatibility tests
            !---------------------
            ! Test 1: area conservation
            zA_cons = SUM(za_i_ini(:,i_hemis)) ; zconv = ABS(zat_i_ini(i_hemis) - zA_cons )
            IF ( zconv .LT. 1.0e-6 ) THEN
               ztest_1 = 1
            ELSE 
               ztest_1 = 0
            ENDIF

            ! Test 2: volume conservation
            zV_cons = SUM(zv_i_ini(:,i_hemis))
            zconv = ABS(zvt_i_ini(i_hemis) - zV_cons)

            IF ( zconv .LT. 1.0e-6 ) THEN
               ztest_2 = 1
            ELSE
               ztest_2 = 0
            ENDIF

            ! Test 3: thickness of the last category is in-bounds ?
            IF ( zh_i_ini(i_fill, i_hemis) > hi_max(i_fill-1) ) THEN
               ztest_3 = 1
            ELSE
               ztest_3 = 0
            ENDIF

            ! Test 4: positivity of ice concentrations
            ztest_4 = 1
            DO jl = 1, jpl
               IF ( za_i_ini(jl,i_hemis) .LT. 0._wp ) THEN 
                  ztest_4 = 0
               ENDIF
            END DO

         ENDIF ! ztest_1 + ztest_2 + ztest_3 + ztest_4
 
         ztests = ztest_1 + ztest_2 + ztest_3 + ztest_4

      END DO ! i_fill

      IF(lwp) THEN 
         WRITE(numout,*) ' ztests : ', ztests
         IF ( ztests .NE. 4 ) THEN
            WRITE(numout,*)
            WRITE(numout,*) ' !!!! ALERT                  !!! '
            WRITE(numout,*) ' !!!! Something is wrong in the LIM3 initialization procedure '
            WRITE(numout,*)
            WRITE(numout,*) ' *** ztests is not equal to 4 '
            WRITE(numout,*) ' *** ztest_i (i=1,4) = ', ztest_1, ztest_2, ztest_3, ztest_4
            WRITE(numout,*) ' zat_i_ini : ', zat_i_ini(i_hemis)
            WRITE(numout,*) ' zht_i_ini : ', zht_i_ini(i_hemis)
         ENDIF ! ztests .NE. 4
      ENDIF
      
      END DO ! i_hemis

      !---------------------------------------------------------------------
      ! 3.3) Space-dependent arrays for ice state variables
      !---------------------------------------------------------------------

      ! Ice concentration, thickness and volume, ice salinity, ice age, surface temperature
      DO jl = 1, jpl ! loop over categories
         DO jj = 1, jpj
            DO ji = 1, jpi
               a_i(ji,jj,jl)   = zswitch(ji,jj) * za_i_ini (jl,zhemis(ji,jj))  ! concentration
               ht_i(ji,jj,jl)  = zswitch(ji,jj) * zh_i_ini(jl,zhemis(ji,jj))   ! ice thickness
               ht_s(ji,jj,jl)  = ht_i(ji,jj,jl) * ( zht_s_ini( zhemis(ji,jj) ) / zht_i_ini( zhemis(ji,jj) ) )  ! snow depth
               sm_i(ji,jj,jl)  = zswitch(ji,jj) * zsm_i_ini(zhemis(ji,jj))     ! salinity
               o_i(ji,jj,jl)   = zswitch(ji,jj) * 0._wp                        ! age (0 day)
               t_su(ji,jj,jl)  = zswitch(ji,jj) * ztm_i_ini(zhemis(ji,jj)) + ( 1._wp - zswitch(ji,jj) ) * rt0 ! surf temp

               ! This case below should not be used if (ht_s/ht_i) is ok in namelist
               ! 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 * ht_s(ji,jj,jl) + ( rhoic - rau0 ) * ht_i(ji,jj,jl) ) * r1_rau0 ) 
               ! recompute ht_i, ht_s avoiding out of bounds values
               ht_i(ji,jj,jl) = MIN( hi_max(jl), ht_i(ji,jj,jl) + zdh )
               ht_s(ji,jj,jl) = MAX( 0._wp, ht_s(ji,jj,jl) - zdh * rhoic * r1_rhosn )

               ! ice volume, salt content, age content
               v_i(ji,jj,jl)   = ht_i(ji,jj,jl) * a_i(ji,jj,jl)              ! ice volume
               v_s(ji,jj,jl)   = ht_s(ji,jj,jl) * a_i(ji,jj,jl)              ! snow volume
               smv_i(ji,jj,jl) = MIN( sm_i(ji,jj,jl) , sss_m(ji,jj) ) * v_i(ji,jj,jl) ! salt content
               oa_i(ji,jj,jl)  = o_i(ji,jj,jl) * a_i(ji,jj,jl)               ! age content
            END DO
         END DO
      END DO

      ! for constant salinity in time
      IF( nn_icesal == 1 .OR. nn_icesal == 3 )  THEN
         CALL lim_var_salprof
         smv_i = sm_i * v_i
      ENDIF

      ! Snow temperature and heat content
      DO jk = 1, nlay_s
         DO jl = 1, jpl ! loop over categories
            DO jj = 1, jpj
               DO ji = 1, jpi
                   t_s(ji,jj,jk,jl) = zswitch(ji,jj) * ztm_i_ini(zhemis(ji,jj)) + ( 1._wp - zswitch(ji,jj) ) * rt0
                   ! Snow energy of melting
                   e_s(ji,jj,jk,jl) = zswitch(ji,jj) * rhosn * ( cpic * ( rt0 - t_s(ji,jj,jk,jl) ) + lfus )

                   ! Mutliply by volume, and divide by number of layers to get heat content in J/m2
                   e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * v_s(ji,jj,jl) * r1_nlay_s
               END DO
            END DO
         END DO
      END DO

      ! Ice salinity, temperature and heat content
      DO jk = 1, nlay_i
         DO jl = 1, jpl ! loop over categories
            DO jj = 1, jpj
               DO ji = 1, jpi
                   t_i(ji,jj,jk,jl) = zswitch(ji,jj) * ztm_i_ini(zhemis(ji,jj)) + ( 1._wp - zswitch(ji,jj) ) * rt0 
                   s_i(ji,jj,jk,jl) = zswitch(ji,jj) * zsm_i_ini(zhemis(ji,jj)) !+ ( 1._wp - zswitch(ji,jj) ) * rn_simin
                   ztmelts          = - tmut * s_i(ji,jj,jk,jl) + rt0 !Melting temperature in K

                   ! heat content per unit volume
                   e_i(ji,jj,jk,jl) = zswitch(ji,jj) * rhoic * (   cpic    * ( ztmelts - t_i(ji,jj,jk,jl) ) &
                      +   lfus    * ( 1._wp - (ztmelts-rt0) / MIN((t_i(ji,jj,jk,jl)-rt0),-epsi20) ) &
                      -   rcp     * ( ztmelts - rt0 ) )

                   ! Mutliply by ice volume, and divide by number of layers to get heat content in J/m2
                   e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * v_i(ji,jj,jl) * r1_nlay_i
               END DO
            END DO
         END DO
      END DO

      tn_ice (:,:,:) = t_su (:,:,:)

      ELSE 
         ! if ln_iceini=false
         a_i  (:,:,:) = 0._wp
         v_i  (:,:,:) = 0._wp
         v_s  (:,:,:) = 0._wp
         smv_i(:,:,:) = 0._wp
         oa_i (:,:,:) = 0._wp
         ht_i (:,:,:) = 0._wp
         ht_s (:,:,:) = 0._wp
         sm_i (:,:,:) = 0._wp
         o_i  (:,:,:) = 0._wp

         e_i(:,:,:,:) = 0._wp
         e_s(:,:,:,:) = 0._wp

         DO jl = 1, jpl
            DO jk = 1, nlay_i
               t_i(:,:,jk,jl) = rt0 * tmask(:,:,1)
            END DO
            DO jk = 1, nlay_s
               t_s(:,:,jk,jl) = rt0 * tmask(:,:,1)
            END DO
         END DO
      
      ENDIF ! ln_iceini
      
      at_i (:,:) = 0.0_wp
      DO jl = 1, jpl
         at_i (:,:) = at_i (:,:) + a_i (:,:,jl)
      END DO
      !
      !--------------------------------------------------------------------
      ! 4) Global ice variables for output diagnostics                    | 
      !--------------------------------------------------------------------
      u_ice (:,:)     = 0._wp
      v_ice (:,:)     = 0._wp
      stress1_i(:,:)  = 0._wp
      stress2_i(:,:)  = 0._wp
      stress12_i(:,:) = 0._wp

      !--------------------------------------------------------------------
      ! 5) Moments for advection
      !--------------------------------------------------------------------

      sxopw (:,:) = 0._wp 
      syopw (:,:) = 0._wp 
      sxxopw(:,:) = 0._wp 
      syyopw(:,:) = 0._wp 
      sxyopw(:,:) = 0._wp

      sxice (:,:,:)  = 0._wp   ;   sxsn (:,:,:)  = 0._wp   ;   sxa  (:,:,:)  = 0._wp
      syice (:,:,:)  = 0._wp   ;   sysn (:,:,:)  = 0._wp   ;   sya  (:,:,:)  = 0._wp
      sxxice(:,:,:)  = 0._wp   ;   sxxsn(:,:,:)  = 0._wp   ;   sxxa (:,:,:)  = 0._wp
      syyice(:,:,:)  = 0._wp   ;   syysn(:,:,:)  = 0._wp   ;   syya (:,:,:)  = 0._wp
      sxyice(:,:,:)  = 0._wp   ;   sxysn(:,:,:)  = 0._wp   ;   sxya (:,:,:)  = 0._wp

      sxc0  (:,:,:)  = 0._wp   ;   sxe  (:,:,:,:)= 0._wp   
      syc0  (:,:,:)  = 0._wp   ;   sye  (:,:,:,:)= 0._wp   
      sxxc0 (:,:,:)  = 0._wp   ;   sxxe (:,:,:,:)= 0._wp   
      syyc0 (:,:,:)  = 0._wp   ;   syye (:,:,:,:)= 0._wp   
      sxyc0 (:,:,:)  = 0._wp   ;   sxye (:,:,:,:)= 0._wp   

      sxsal  (:,:,:)  = 0._wp
      sysal  (:,:,:)  = 0._wp
      sxxsal (:,:,:)  = 0._wp
      syysal (:,:,:)  = 0._wp
      sxysal (:,:,:)  = 0._wp

      sxage  (:,:,:)  = 0._wp
      syage  (:,:,:)  = 0._wp
      sxxage (:,:,:)  = 0._wp
      syyage (:,:,:)  = 0._wp
      sxyage (:,:,:)  = 0._wp


      CALL wrk_dealloc( jpi, jpj, zswitch )
      CALL wrk_dealloc( jpi, jpj, zhemis )
      CALL wrk_dealloc( jpl,   2, zh_i_ini,  za_i_ini,  zv_i_ini )
      CALL wrk_dealloc(   2,      zht_i_ini, zat_i_ini, zvt_i_ini, zht_s_ini, zsm_i_ini, ztm_i_ini )

   END SUBROUTINE lim_istate

   SUBROUTINE lim_istate_init
      !!-------------------------------------------------------------------
      !!                   ***  ROUTINE lim_istate_init  ***
      !!        
      !! ** Purpose : Definition of initial state of the ice 
      !!
      !! ** Method : Read the namiceini namelist and check the parameter 
      !!       values called at the first timestep (nit000)
      !!
      !! ** input : 
      !!        Namelist namiceini
      !!
      !! history :
      !!  8.5  ! 03-08 (C. Ethe) original code 
      !!  8.5  ! 07-11 (M. Vancoppenolle) rewritten initialization
      !!-----------------------------------------------------------------------------
      NAMELIST/namiceini/ ln_iceini, rn_thres_sst, rn_hts_ini_n, rn_hts_ini_s, rn_hti_ini_n, rn_hti_ini_s,  &
         &                                      rn_ati_ini_n, rn_ati_ini_s, rn_smi_ini_n, rn_smi_ini_s, rn_tmi_ini_n, rn_tmi_ini_s
      INTEGER :: ios                 ! Local integer output status for namelist read
      !!-----------------------------------------------------------------------------
      !
      REWIND( numnam_ice_ref )              ! Namelist namiceini in reference namelist : Ice initial state
      READ  ( numnam_ice_ref, namiceini, IOSTAT = ios, ERR = 901)
901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceini in reference namelist', lwp )

      REWIND( numnam_ice_cfg )              ! Namelist namiceini in configuration namelist : Ice initial state
      READ  ( numnam_ice_cfg, namiceini, IOSTAT = ios, ERR = 902 )
902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceini in configuration namelist', lwp )
      IF(lwm) WRITE ( numoni, namiceini )

      ! Define the initial parameters
      ! -------------------------

      IF(lwp) THEN
         WRITE(numout,*)
         WRITE(numout,*) 'lim_istate_init : ice parameters inititialisation '
         WRITE(numout,*) '~~~~~~~~~~~~~~~'
         WRITE(numout,*) '   initialization with ice (T) or not (F)       ln_iceini     = ', ln_iceini
         WRITE(numout,*) '   threshold water temp. for initial sea-ice    rn_thres_sst  = ', rn_thres_sst
         WRITE(numout,*) '   initial snow thickness in the north          rn_hts_ini_n  = ', rn_hts_ini_n
         WRITE(numout,*) '   initial snow thickness in the south          rn_hts_ini_s  = ', rn_hts_ini_s 
         WRITE(numout,*) '   initial ice thickness  in the north          rn_hti_ini_n  = ', rn_hti_ini_n
         WRITE(numout,*) '   initial ice thickness  in the south          rn_hti_ini_s  = ', rn_hti_ini_s
         WRITE(numout,*) '   initial ice concentr.  in the north          rn_ati_ini_n  = ', rn_ati_ini_n
         WRITE(numout,*) '   initial ice concentr.  in the north          rn_ati_ini_s  = ', rn_ati_ini_s
         WRITE(numout,*) '   initial  ice salinity  in the north          rn_smi_ini_n  = ', rn_smi_ini_n
         WRITE(numout,*) '   initial  ice salinity  in the south          rn_smi_ini_s  = ', rn_smi_ini_s
         WRITE(numout,*) '   initial  ice/snw temp  in the north          rn_tmi_ini_n  = ', rn_tmi_ini_n
         WRITE(numout,*) '   initial  ice/snw temp  in the south          rn_tmi_ini_s  = ', rn_tmi_ini_s
      ENDIF

   END SUBROUTINE lim_istate_init

#else
   !!----------------------------------------------------------------------
   !!   Default option :         Empty module          NO LIM sea-ice model
   !!----------------------------------------------------------------------
CONTAINS
   SUBROUTINE lim_istate          ! Empty routine
   END SUBROUTINE lim_istate
#endif

   !!======================================================================
END MODULE limistate