Skip to content
Snippets Groups Projects
iceistate.F90 29.7 KiB
Newer Older
Guillaume Samson's avatar
Guillaume Samson committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438
MODULE iceistate
   !!======================================================================
   !!                     ***  MODULE  iceistate  ***
   !!   sea-ice : Initialization of ice variables
   !!======================================================================
   !! History :  2.0  !  2004-01  (C. Ethe, G. Madec) Original code
   !!            3.0  !  2007     (M. Vancoppenolle)  Rewrite for ice cats
   !!            4.0  !  2018     (many people)       SI3 [aka Sea Ice cube]
   !!----------------------------------------------------------------------
#if defined key_si3
   !!----------------------------------------------------------------------
   !!   'key_si3'                                       SI3 sea-ice model
   !!----------------------------------------------------------------------
   !!   ice_istate       :  initialization of diagnostics ice variables
   !!   ice_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 , ONLY : sst_m, sss_m, ln_ice_embd 
   USE sbc_ice , ONLY : tn_ice, snwice_mass, snwice_mass_b
   USE eosbn2         ! equation of state
# if defined key_qco
   USE domqco         ! Quasi-Eulerian coord.
# elif defined key_linssh
   !                  ! Fix in time coord.
# else
   USE domvvl         ! Variable volume
# endif
   USE ice            ! sea-ice: variables
   USE ice1D          ! sea-ice: thermodynamics variables
   USE icetab         ! sea-ice: 1D <==> 2D transformation
   USE icevar         ! sea-ice: operations
   !
   USE in_out_manager ! I/O manager
   USE iom            ! I/O manager library
   USE lib_mpp        ! MPP library
   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero)
   USE fldread        ! read input fields

# if defined key_agrif
   USE agrif_oce
   USE agrif_ice_interp 
# endif   

   IMPLICIT NONE
   PRIVATE

   PUBLIC   ice_istate        ! called by icestp.F90
   PUBLIC   ice_istate_init   ! called by icestp.F90
   !
   !                             !! ** namelist (namini) **
   LOGICAL, PUBLIC  ::   ln_iceini        !: Ice initialization or not
   INTEGER, PUBLIC  ::   nn_iceini_file   !: Ice initialization:
                                  !        0 = Initialise sea ice based on SSTs
                                  !        1 = Initialise sea ice from single category netcdf file
                                  !        2 = Initialise sea ice from multi category restart file
   REAL(wp) ::   rn_thres_sst
   REAL(wp) ::   rn_hti_ini_n, rn_hts_ini_n, rn_ati_ini_n, rn_smi_ini_n, rn_tmi_ini_n, rn_tsu_ini_n, rn_tms_ini_n
   REAL(wp) ::   rn_hti_ini_s, rn_hts_ini_s, rn_ati_ini_s, rn_smi_ini_s, rn_tmi_ini_s, rn_tsu_ini_s, rn_tms_ini_s
   REAL(wp) ::   rn_apd_ini_n, rn_hpd_ini_n, rn_hld_ini_n
   REAL(wp) ::   rn_apd_ini_s, rn_hpd_ini_s, rn_hld_ini_s
   !
   !                              ! if nn_iceini_file = 1
   INTEGER , PARAMETER ::   jpfldi = 10          ! maximum number of files to read
   INTEGER , PARAMETER ::   jp_hti = 1           ! index of ice thickness    (m)
   INTEGER , PARAMETER ::   jp_hts = 2           ! index of snw thickness    (m)
   INTEGER , PARAMETER ::   jp_ati = 3           ! index of ice fraction     (-)
   INTEGER , PARAMETER ::   jp_smi = 4           ! index of ice salinity     (g/kg)
   INTEGER , PARAMETER ::   jp_tmi = 5           ! index of ice temperature  (K)
   INTEGER , PARAMETER ::   jp_tsu = 6           ! index of ice surface temp (K)
   INTEGER , PARAMETER ::   jp_tms = 7           ! index of snw temperature  (K)
   INTEGER , PARAMETER ::   jp_apd = 8           ! index of pnd fraction     (-)
   INTEGER , PARAMETER ::   jp_hpd = 9           ! index of pnd depth        (m)
   INTEGER , PARAMETER ::   jp_hld = 10          ! index of pnd lid depth    (m)
   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   si  ! structure of input fields (file informations, fields read)
   !
#if defined key_agrif
   REAL(wp), PUBLIC ::   rsshadj   !: initial mean ssh adjustment due to initial ice+snow mass
#endif
   !
   !! * Substitutions
#  include "do_loop_substitute.h90"
   !!----------------------------------------------------------------------
   !! NEMO/ICE 4.0 , NEMO Consortium (2018)
   !! $Id: iceistate.F90 15530 2021-11-23 15:09:32Z clem $
   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
   !!----------------------------------------------------------------------
CONTAINS

   SUBROUTINE ice_istate( kt, Kbb, Kmm, Kaa )
      !!-------------------------------------------------------------------
      !!                    ***  ROUTINE ice_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) Set initial surface and basal temperatures
      !!                2) Recompute or read sea ice state variables
      !!                3) Fill in space-dependent arrays for state variables
      !!                4) snow-ice mass computation
      !!
      !! ** Notes   : o_i, t_su, t_s, t_i, sz_i must be filled everywhere, even
      !!              where there is no ice
      !!--------------------------------------------------------------------
      INTEGER, INTENT(in) :: kt            ! time step 
      INTEGER, INTENT(in) :: Kbb, Kmm, Kaa ! ocean time level indices
      !
      INTEGER  ::   ji, jj, jk, jl         ! dummy loop indices
      REAL(wp) ::   ztmelts, zsshadj, area
      INTEGER , DIMENSION(4)           ::   itest
      REAL(wp), DIMENSION(jpi,jpj)     ::   z2d
      REAL(wp), DIMENSION(jpi,jpj)     ::   zswitch    ! ice indicator
      REAL(wp), DIMENSION(jpi,jpj)     ::   zht_i_ini, zat_i_ini, ztm_s_ini            !data from namelist or nc file
      REAL(wp), DIMENSION(jpi,jpj)     ::   zt_su_ini, zht_s_ini, zsm_i_ini, ztm_i_ini !data from namelist or nc file
      REAL(wp), DIMENSION(jpi,jpj)     ::   zapnd_ini, zhpnd_ini, zhlid_ini            !data from namelist or nc file
      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zti_3d , zts_3d                            !temporary arrays
      !!
      REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   zhi_2d, zhs_2d, zai_2d, zti_2d, zts_2d, ztsu_2d, zsi_2d, zaip_2d, zhip_2d, zhil_2d
      !--------------------------------------------------------------------

      IF(lwp) WRITE(numout,*)
      IF(lwp) WRITE(numout,*) 'ice_istate: sea-ice initialization '
      IF(lwp) WRITE(numout,*) '~~~~~~~~~~'

      !---------------------------
      ! 1) 1st init. of the fields
      !---------------------------
      !
      ! basal temperature (considered at freezing point)   [Kelvin]
      CALL eos_fzp( sss_m(:,:), t_bo(:,:) )
      t_bo(:,:) = ( t_bo(:,:) + rt0 ) * tmask(:,:,1) 
      !
      ! surface temperature and conductivity
      DO jl = 1, jpl
         t_su   (:,:,jl) = rt0 * tmask(:,:,1)  ! temp at the surface
         cnd_ice(:,:,jl) = 0._wp               ! initialisation of the effective conductivity at the top of ice/snow (ln_cndflx=T)
      END DO
      !
      ! ice and snw temperatures
      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
      !
      ! specific temperatures for coupled runs
      tn_ice (:,:,:) = t_i (:,:,1,:)
      t1_ice (:,:,:) = t_i (:,:,1,:)

      ! heat contents
      e_i (:,:,:,:) = 0._wp
      e_s (:,:,:,:) = 0._wp
      
      ! general fields
      a_i (:,:,:) = 0._wp
      v_i (:,:,:) = 0._wp
      v_s (:,:,:) = 0._wp
      sv_i(:,:,:) = 0._wp
      oa_i(:,:,:) = 0._wp
      !
      h_i (:,:,:) = 0._wp
      h_s (:,:,:) = 0._wp
      s_i (:,:,:) = 0._wp
      o_i (:,:,:) = 0._wp
      !
      ! melt ponds
      a_ip     (:,:,:) = 0._wp
      v_ip     (:,:,:) = 0._wp
      v_il     (:,:,:) = 0._wp
      a_ip_eff (:,:,:) = 0._wp
      h_ip     (:,:,:) = 0._wp
      h_il     (:,:,:) = 0._wp
      !
      ! ice velocities
      u_ice (:,:) = 0._wp
      v_ice (:,:) = 0._wp
      !
      !------------------------------------------------------------------------
      ! 2) overwrite some of the fields with namelist parameters or netcdf file
      !------------------------------------------------------------------------
      IF( ln_iceini ) THEN
         !
#if defined key_agrif
         IF ( ( Agrif_Root() ).OR.(.NOT.ln_init_chfrpar ) ) THEN
#endif
            !                             !---------------!
            IF( nn_iceini_file == 1 )THEN ! Read a file   !
               !                          !---------------!
               WHERE( ff_t(:,:) >= 0._wp )   ;   zswitch(:,:) = 1._wp
               ELSEWHERE                     ;   zswitch(:,:) = 0._wp
               END WHERE
               !
               CALL fld_read( kt, 1, si ) ! input fields provided at the current time-step
               !
               ! -- mandatory fields -- !
               zht_i_ini(:,:) = si(jp_hti)%fnow(:,:,1) * tmask(:,:,1)
               zht_s_ini(:,:) = si(jp_hts)%fnow(:,:,1) * tmask(:,:,1)
               zat_i_ini(:,:) = si(jp_ati)%fnow(:,:,1) * tmask(:,:,1)

               ! -- optional fields -- !
               !    if fields do not exist then set them to the values present in the namelist (except for temperatures)
               !
               ! ice salinity
               IF( TRIM(si(jp_smi)%clrootname) == 'NOT USED' ) &
                  &     si(jp_smi)%fnow(:,:,1) = ( rn_smi_ini_n * zswitch + rn_smi_ini_s * (1._wp - zswitch) ) * tmask(:,:,1)
               !
               ! temperatures
               IF    ( TRIM(si(jp_tmi)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tsu)%clrootname) == 'NOT USED' .AND. &
                  &    TRIM(si(jp_tms)%clrootname) == 'NOT USED' ) THEN
                  si(jp_tmi)%fnow(:,:,1) = ( rn_tmi_ini_n * zswitch + rn_tmi_ini_s * (1._wp - zswitch) ) * tmask(:,:,1)
                  si(jp_tsu)%fnow(:,:,1) = ( rn_tsu_ini_n * zswitch + rn_tsu_ini_s * (1._wp - zswitch) ) * tmask(:,:,1)
                  si(jp_tms)%fnow(:,:,1) = ( rn_tms_ini_n * zswitch + rn_tms_ini_s * (1._wp - zswitch) ) * tmask(:,:,1)
               ENDIF
               IF( TRIM(si(jp_tmi)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tms)%clrootname) /= 'NOT USED' ) & ! if T_s is read and not T_i, set T_i = (T_s + T_freeze)/2
                  &     si(jp_tmi)%fnow(:,:,1) = 0.5_wp * ( si(jp_tms)%fnow(:,:,1) + 271.15 )
               IF( TRIM(si(jp_tmi)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tsu)%clrootname) /= 'NOT USED' ) & ! if T_su is read and not T_i, set T_i = (T_su + T_freeze)/2
                  &     si(jp_tmi)%fnow(:,:,1) = 0.5_wp * ( si(jp_tsu)%fnow(:,:,1) + 271.15 )
               IF( TRIM(si(jp_tsu)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tms)%clrootname) /= 'NOT USED' ) & ! if T_s is read and not T_su, set T_su = T_s
                  &     si(jp_tsu)%fnow(:,:,1) = si(jp_tms)%fnow(:,:,1)
               IF( TRIM(si(jp_tsu)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tmi)%clrootname) /= 'NOT USED' ) & ! if T_i is read and not T_su, set T_su = T_i
                  &     si(jp_tsu)%fnow(:,:,1) = si(jp_tmi)%fnow(:,:,1)
               IF( TRIM(si(jp_tms)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tsu)%clrootname) /= 'NOT USED' ) & ! if T_su is read and not T_s, set T_s = T_su
                  &     si(jp_tms)%fnow(:,:,1) = si(jp_tsu)%fnow(:,:,1)
               IF( TRIM(si(jp_tms)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tmi)%clrootname) /= 'NOT USED' ) & ! if T_i is read and not T_s, set T_s = T_i
                  &     si(jp_tms)%fnow(:,:,1) = si(jp_tmi)%fnow(:,:,1)
               !
               ! pond concentration
               IF( TRIM(si(jp_apd)%clrootname) == 'NOT USED' ) &
                  &     si(jp_apd)%fnow(:,:,1) = ( rn_apd_ini_n * zswitch + rn_apd_ini_s * (1._wp - zswitch) ) * tmask(:,:,1) & ! rn_apd = pond fraction => rn_apnd * a_i = pond conc.
                  &                              * si(jp_ati)%fnow(:,:,1) 
               !
               ! pond depth
               IF( TRIM(si(jp_hpd)%clrootname) == 'NOT USED' ) &
                  &     si(jp_hpd)%fnow(:,:,1) = ( rn_hpd_ini_n * zswitch + rn_hpd_ini_s * (1._wp - zswitch) ) * tmask(:,:,1)
               !
               ! pond lid depth
               IF( TRIM(si(jp_hld)%clrootname) == 'NOT USED' ) &
                  &     si(jp_hld)%fnow(:,:,1) = ( rn_hld_ini_n * zswitch + rn_hld_ini_s * (1._wp - zswitch) ) * tmask(:,:,1)
               !
               zsm_i_ini(:,:) = si(jp_smi)%fnow(:,:,1) * tmask(:,:,1)
               ztm_i_ini(:,:) = si(jp_tmi)%fnow(:,:,1) * tmask(:,:,1)
               zt_su_ini(:,:) = si(jp_tsu)%fnow(:,:,1) * tmask(:,:,1)
               ztm_s_ini(:,:) = si(jp_tms)%fnow(:,:,1) * tmask(:,:,1)
               zapnd_ini(:,:) = si(jp_apd)%fnow(:,:,1) * tmask(:,:,1)
               zhpnd_ini(:,:) = si(jp_hpd)%fnow(:,:,1) * tmask(:,:,1)
               zhlid_ini(:,:) = si(jp_hld)%fnow(:,:,1) * tmask(:,:,1)
               !
               ! change the switch for the following
               WHERE( zat_i_ini(:,:) > 0._wp )   ;   zswitch(:,:) = tmask(:,:,1) 
               ELSEWHERE                         ;   zswitch(:,:) = 0._wp
               END WHERE

               !                          !---------------!
            ELSE                          ! Read namelist !
               !                          !---------------!
               ! no ice if (sst - Tfreez) >= thresold
               WHERE( ( sst_m(:,:) - (t_bo(:,:) - rt0) ) * tmask(:,:,1) >= rn_thres_sst )   ;   zswitch(:,:) = 0._wp 
               ELSEWHERE                                                                    ;   zswitch(:,:) = tmask(:,:,1)
               END WHERE
               !
               ! assign initial thickness, concentration, snow depth and salinity to an hemisphere-dependent array
               WHERE( ff_t(:,:) >= 0._wp )
                  zht_i_ini(:,:) = rn_hti_ini_n * zswitch(:,:)
                  zht_s_ini(:,:) = rn_hts_ini_n * zswitch(:,:)
                  zat_i_ini(:,:) = rn_ati_ini_n * zswitch(:,:)
                  zsm_i_ini(:,:) = rn_smi_ini_n * zswitch(:,:)
                  ztm_i_ini(:,:) = rn_tmi_ini_n * zswitch(:,:)
                  zt_su_ini(:,:) = rn_tsu_ini_n * zswitch(:,:)
                  ztm_s_ini(:,:) = rn_tms_ini_n * zswitch(:,:)
                  zapnd_ini(:,:) = rn_apd_ini_n * zswitch(:,:) * zat_i_ini(:,:) ! rn_apd = pond fraction => rn_apd * a_i = pond conc. 
                  zhpnd_ini(:,:) = rn_hpd_ini_n * zswitch(:,:)
                  zhlid_ini(:,:) = rn_hld_ini_n * zswitch(:,:)
               ELSEWHERE
                  zht_i_ini(:,:) = rn_hti_ini_s * zswitch(:,:)
                  zht_s_ini(:,:) = rn_hts_ini_s * zswitch(:,:)
                  zat_i_ini(:,:) = rn_ati_ini_s * zswitch(:,:)
                  zsm_i_ini(:,:) = rn_smi_ini_s * zswitch(:,:)
                  ztm_i_ini(:,:) = rn_tmi_ini_s * zswitch(:,:)
                  zt_su_ini(:,:) = rn_tsu_ini_s * zswitch(:,:)
                  ztm_s_ini(:,:) = rn_tms_ini_s * zswitch(:,:)
                  zapnd_ini(:,:) = rn_apd_ini_s * zswitch(:,:) * zat_i_ini(:,:) ! rn_apd = pond fraction => rn_apd * a_i = pond conc.
                  zhpnd_ini(:,:) = rn_hpd_ini_s * zswitch(:,:)
                  zhlid_ini(:,:) = rn_hld_ini_s * zswitch(:,:)
               END WHERE
               !
            ENDIF



            ! make sure ponds = 0 if no ponds scheme
            IF ( .NOT.ln_pnd ) THEN
               zapnd_ini(:,:) = 0._wp
               zhpnd_ini(:,:) = 0._wp
               zhlid_ini(:,:) = 0._wp
            ENDIF
            
            IF ( .NOT.ln_pnd_lids ) THEN
               zhlid_ini(:,:) = 0._wp
            ENDIF
            
            !----------------!
            ! 3) fill fields !
            !----------------!
            ! select ice covered grid points
            npti = 0 ; nptidx(:) = 0
            DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
               IF ( zht_i_ini(ji,jj) > 0._wp ) THEN
                  npti         = npti  + 1
                  nptidx(npti) = (jj - 1) * jpi + ji
               ENDIF
            END_2D

            ! move to 1D arrays: (jpi,jpj) -> (jpi*jpj)
            CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d (1:npti)  , zht_i_ini )
            CALL tab_2d_1d( npti, nptidx(1:npti), h_s_1d (1:npti)  , zht_s_ini )
            CALL tab_2d_1d( npti, nptidx(1:npti), at_i_1d(1:npti)  , zat_i_ini )
            CALL tab_2d_1d( npti, nptidx(1:npti), t_i_1d (1:npti,1), ztm_i_ini )
            CALL tab_2d_1d( npti, nptidx(1:npti), t_s_1d (1:npti,1), ztm_s_ini )
            CALL tab_2d_1d( npti, nptidx(1:npti), t_su_1d(1:npti)  , zt_su_ini )
            CALL tab_2d_1d( npti, nptidx(1:npti), s_i_1d (1:npti)  , zsm_i_ini )
            CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_1d(1:npti)  , zapnd_ini )
            CALL tab_2d_1d( npti, nptidx(1:npti), h_ip_1d(1:npti)  , zhpnd_ini )
            CALL tab_2d_1d( npti, nptidx(1:npti), h_il_1d(1:npti)  , zhlid_ini )
            
            ! allocate temporary arrays
            ALLOCATE( zhi_2d (npti,jpl), zhs_2d (npti,jpl), zai_2d (npti,jpl), &
               &      zti_2d (npti,jpl), zts_2d (npti,jpl), ztsu_2d(npti,jpl), zsi_2d(npti,jpl), &
               &      zaip_2d(npti,jpl), zhip_2d(npti,jpl), zhil_2d(npti,jpl) )

            ! distribute 1-cat into jpl-cat: (jpi*jpj) -> (jpi*jpj,jpl)
            CALL ice_var_itd( h_i_1d(1:npti)  , h_s_1d(1:npti)  , at_i_1d(1:npti),                  &
               &              zhi_2d          , zhs_2d          , zai_2d         ,                  &
               &              t_i_1d(1:npti,1), t_s_1d(1:npti,1), t_su_1d(1:npti),                  &
               &              s_i_1d(1:npti)  , a_ip_1d(1:npti) , h_ip_1d(1:npti), h_il_1d(1:npti), &
               &              zti_2d          , zts_2d          , ztsu_2d        ,                  &
               &              zsi_2d          , zaip_2d         , zhip_2d        , zhil_2d )

            ! move to 3D arrays: (jpi*jpj,jpl) -> (jpi,jpj,jpl)
            DO jl = 1, jpl
               zti_3d(:,:,jl) = rt0 * tmask(:,:,1)
               zts_3d(:,:,jl) = rt0 * tmask(:,:,1)
            END DO
            CALL tab_2d_3d( npti, nptidx(1:npti), zhi_2d   , h_i    )
            CALL tab_2d_3d( npti, nptidx(1:npti), zhs_2d   , h_s    )
            CALL tab_2d_3d( npti, nptidx(1:npti), zai_2d   , a_i    )
            CALL tab_2d_3d( npti, nptidx(1:npti), zti_2d   , zti_3d )
            CALL tab_2d_3d( npti, nptidx(1:npti), zts_2d   , zts_3d )
            CALL tab_2d_3d( npti, nptidx(1:npti), ztsu_2d  , t_su   )
            CALL tab_2d_3d( npti, nptidx(1:npti), zsi_2d   , s_i    )
            CALL tab_2d_3d( npti, nptidx(1:npti), zaip_2d  , a_ip   )
            CALL tab_2d_3d( npti, nptidx(1:npti), zhip_2d  , h_ip   )
            CALL tab_2d_3d( npti, nptidx(1:npti), zhil_2d  , h_il   )

            ! deallocate temporary arrays
            DEALLOCATE( zhi_2d, zhs_2d, zai_2d , &
               &        zti_2d, zts_2d, ztsu_2d, zsi_2d, zaip_2d, zhip_2d, zhil_2d )

            ! calculate extensive and intensive variables
            CALL ice_var_salprof ! for sz_i
            DO jl = 1, jpl
               DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
                  v_i (ji,jj,jl) = h_i(ji,jj,jl) * a_i(ji,jj,jl)
                  v_s (ji,jj,jl) = h_s(ji,jj,jl) * a_i(ji,jj,jl)
                  sv_i(ji,jj,jl) = MIN( MAX( rn_simin , s_i(ji,jj,jl) ) , rn_simax ) * v_i(ji,jj,jl)
               END_2D
            END DO
            !
            DO jl = 1, jpl
               DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, nlay_s )
                  t_s(ji,jj,jk,jl) = zts_3d(ji,jj,jl)
                  e_s(ji,jj,jk,jl) = zswitch(ji,jj) * v_s(ji,jj,jl) * r1_nlay_s * &
                     &               rhos * ( rcpi * ( rt0 - t_s(ji,jj,jk,jl) ) + rLfus )
               END_3D
            END DO
            !
            DO jl = 1, jpl
               DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, nlay_i )
                  t_i (ji,jj,jk,jl) = zti_3d(ji,jj,jl) 
                  ztmelts          = - rTmlt * sz_i(ji,jj,jk,jl) + rt0 ! melting temperature in K
                  e_i(ji,jj,jk,jl) = zswitch(ji,jj) * v_i(ji,jj,jl) * r1_nlay_i * &
                     &               rhoi * (  rcpi  * ( ztmelts - t_i(ji,jj,jk,jl) ) + &
                     &                         rLfus * ( 1._wp - (ztmelts-rt0) / MIN( (t_i(ji,jj,jk,jl)-rt0), -epsi20 ) ) &
                     &                       - rcp   * ( ztmelts - rt0 ) )
               END_3D
            END DO
            !
#if defined key_agrif
         ELSE
            CALL  agrif_istate_ice
         ENDIF
#endif
         ! Melt ponds
         WHERE( a_i > epsi10 )   ;   a_ip_eff(:,:,:) = a_ip(:,:,:) / a_i(:,:,:)
         ELSEWHERE               ;   a_ip_eff(:,:,:) = 0._wp
         END WHERE
         v_ip(:,:,:) = h_ip(:,:,:) * a_ip(:,:,:)
         v_il(:,:,:) = h_il(:,:,:) * a_ip(:,:,:)
         
         ! specific temperatures for coupled runs
         tn_ice(:,:,:) = t_su(:,:,:)
         t1_ice(:,:,:) = t_i (:,:,1,:)
         !
         ! ice concentration should not exceed amax
         at_i(:,:) = SUM( a_i, dim=3 )
         DO jl = 1, jpl
            WHERE( at_i(:,:) > rn_amax_2d(:,:) )   a_i(:,:,jl) = a_i(:,:,jl) * rn_amax_2d(:,:) / at_i(:,:)
         END DO
        at_i(:,:) = SUM( a_i, dim=3 )
         !
      ENDIF ! ln_iceini
      !
      !----------------------------------------------------------
      ! 4) Adjust ssh and vertical scale factors to snow-ice mass
      !----------------------------------------------------------
      snwice_mass  (:,:) = tmask(:,:,1) * SUM( rhos * v_s + rhoi * v_i + rhow * ( v_ip + v_il ), dim=3  )   ! snow+ice mass
      snwice_mass_b(:,:) = snwice_mass(:,:)
      !
      IF( ln_ice_embd ) THEN            ! embedded sea-ice: deplete the initial ssh below sea-ice area
         !                              ! ----------------
         ssh(:,:,Kmm) = ssh(:,:,Kmm) - snwice_mass(:,:) * r1_rho0
         ssh(:,:,Kbb) = ssh(:,:,Kbb) - snwice_mass(:,:) * r1_rho0
         !
      ELSE                              ! levitating sea-ice: deplete the initial ssh over the whole domain
         !                              ! ------------------
         area    = glob_sum( 'iceistate', e1e2t(:,:) * ssmask(:,:) )
         zsshadj = glob_sum( 'iceistate', snwice_mass(:,:) * r1_rho0 * e1e2t(:,:) ) / area
#if defined key_agrif
         ! Override ssh adjustment in nested domains by the root-domain ssh adjustment;
         ! store the adjustment value in a global module variable to make it retrievable in nested domains
         IF( .NOT.Agrif_Root() ) THEN
Guillaume Samson's avatar
Guillaume Samson committed
            zsshadj = Agrif_Parent(rsshadj)
Guillaume Samson's avatar
Guillaume Samson committed
         ELSE
            rsshadj = zsshadj
         ENDIF
#endif
         IF(lwp) WRITE(numout,'(A23,F10.6,A20)') ' sea level adjusted by ', -zsshadj, ' m to compensate for'
         IF(lwp) WRITE(numout,*) ' the initial snow+ice mass'
         !
         WHERE( ssmask(:,:) == 1._wp )
            ssh(:,:,Kmm) = ssh(:,:,Kmm) - zsshadj
            ssh(:,:,Kbb) = ssh(:,:,Kbb) - zsshadj
         ENDWHERE
         !
      ENDIF
      !
      IF( .NOT.ln_linssh ) THEN
#if defined key_qco
         CALL dom_qco_zgr( Kbb, Kmm )        ! upadte of r3=ssh/h0 ratios
#elif defined key_linssh
         !                                   ! Fix in time : key_linssh case, set through domzgr_substitute.h90
#else
         DO jk = 1, jpk
            DO_2D( nn_hls, nn_hls, nn_hls, nn_hls)
Guillaume Samson's avatar
Guillaume Samson committed
!               IF( snwice_mass(ji,jj) /= 0._wp ) THEN
Guillaume Samson's avatar
Guillaume Samson committed
                  e3t(ji,jj,jk,Kmm) = e3t_0(ji,jj,jk) * ( 1._wp + ssh(ji,jj,Kmm) * r1_ht_0(ji,jj) * tmask(ji,jj,jk) )
                  e3t(ji,jj,jk,Kbb) = e3t_0(ji,jj,jk) * ( 1._wp + ssh(ji,jj,Kbb) * r1_ht_0(ji,jj) * tmask(ji,jj,jk) )
Guillaume Samson's avatar
Guillaume Samson committed
!               ENDIF
Guillaume Samson's avatar
Guillaume Samson committed
            END_2D
         END DO
         !
         CALL dom_vvl_zgr( Kbb, Kmm, Kaa )   ! interpolation of scale factor, depth and water column
#endif
      ENDIF
      !
      !!clem: output of initial state should be written here but it is impossible because
      !!      the ocean and ice are in the same file
      !!      CALL dia_wri_state( 'output.init' )
      !
   END SUBROUTINE ice_istate


   SUBROUTINE ice_istate_init
      !!-------------------------------------------------------------------
      !!                   ***  ROUTINE ice_istate_init  ***
      !!        
      !! ** Purpose :   Definition of initial state of the ice 
      !!
      !! ** Method  :   Read the namini namelist and check the parameter 
      !!              values called at the first timestep (nit000)
      !!
      !! ** input   :  Namelist namini
      !!
      !!-----------------------------------------------------------------------------
      INTEGER ::   ios   ! Local integer output status for namelist read
      INTEGER ::   ifpr, ierror
      !
      CHARACTER(len=256) ::  cn_dir          ! Root directory for location of ice files
      TYPE(FLD_N)                    ::   sn_hti, sn_hts, sn_ati, sn_smi, sn_tmi, sn_tsu, sn_tms, sn_apd, sn_hpd, sn_hld
      TYPE(FLD_N), DIMENSION(jpfldi) ::   slf_i                 ! array of namelist informations on the fields to read
      !
      NAMELIST/namini/ ln_iceini, nn_iceini_file, rn_thres_sst, &
         &             rn_hti_ini_n, rn_hti_ini_s, rn_hts_ini_n, rn_hts_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, rn_tsu_ini_n, rn_tsu_ini_s, rn_tms_ini_n, rn_tms_ini_s, &
         &             rn_apd_ini_n, rn_apd_ini_s, rn_hpd_ini_n, rn_hpd_ini_s, rn_hld_ini_n, rn_hld_ini_s, &
         &             sn_hti, sn_hts, sn_ati, sn_tsu, sn_tmi, sn_smi, sn_tms, sn_apd, sn_hpd, sn_hld, cn_dir
      !!-----------------------------------------------------------------------------
      !
      READ  ( numnam_ice_ref, namini, IOSTAT = ios, ERR = 901)
901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namini in reference namelist' )
      READ  ( numnam_ice_cfg, namini, IOSTAT = ios, ERR = 902 )
902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namini in configuration namelist' )
      IF(lwm) WRITE ( numoni, namini )
      !
      slf_i(jp_hti) = sn_hti  ;  slf_i(jp_hts) = sn_hts
      slf_i(jp_ati) = sn_ati  ;  slf_i(jp_smi) = sn_smi
      slf_i(jp_tmi) = sn_tmi  ;  slf_i(jp_tsu) = sn_tsu   ;   slf_i(jp_tms) = sn_tms
      slf_i(jp_apd) = sn_apd  ;  slf_i(jp_hpd) = sn_hpd   ;   slf_i(jp_hld) = sn_hld
      !
      IF(lwp) THEN                          ! control print
         WRITE(numout,*)
         WRITE(numout,*) 'ice_istate_init: ice parameters inititialisation '
         WRITE(numout,*) '~~~~~~~~~~~~~~~'
         WRITE(numout,*) '   Namelist namini:'
         WRITE(numout,*) '      ice initialization (T) or not (F)                ln_iceini      = ', ln_iceini
         WRITE(numout,*) '      ice initialization from a netcdf file            nn_iceini_file = ', nn_iceini_file
         WRITE(numout,*) '      max ocean temp. above Tfreeze with initial ice   rn_thres_sst   = ', rn_thres_sst
         IF( ln_iceini .AND. nn_iceini_file == 0 ) THEN
            WRITE(numout,*) '      initial snw thickness in the north-south         rn_hts_ini     = ', rn_hts_ini_n,rn_hts_ini_s 
            WRITE(numout,*) '      initial ice thickness in the north-south         rn_hti_ini     = ', rn_hti_ini_n,rn_hti_ini_s
            WRITE(numout,*) '      initial ice concentr  in the north-south         rn_ati_ini     = ', rn_ati_ini_n,rn_ati_ini_s
            WRITE(numout,*) '      initial ice salinity  in the north-south         rn_smi_ini     = ', rn_smi_ini_n,rn_smi_ini_s
            WRITE(numout,*) '      initial surf temperat in the north-south         rn_tsu_ini     = ', rn_tsu_ini_n,rn_tsu_ini_s
            WRITE(numout,*) '      initial ice temperat  in the north-south         rn_tmi_ini     = ', rn_tmi_ini_n,rn_tmi_ini_s
            WRITE(numout,*) '      initial snw temperat  in the north-south         rn_tms_ini     = ', rn_tms_ini_n,rn_tms_ini_s
            WRITE(numout,*) '      initial pnd fraction  in the north-south         rn_apd_ini     = ', rn_apd_ini_n,rn_apd_ini_s
            WRITE(numout,*) '      initial pnd depth     in the north-south         rn_hpd_ini     = ', rn_hpd_ini_n,rn_hpd_ini_s
            WRITE(numout,*) '      initial pnd lid depth in the north-south         rn_hld_ini     = ', rn_hld_ini_n,rn_hld_ini_s
         ENDIF
      ENDIF
      !
      IF( nn_iceini_file == 1 ) THEN                      ! Ice initialization using input file
         !
         ! set si structure
         ALLOCATE( si(jpfldi), STAT=ierror )
         IF( ierror > 0 ) THEN
            CALL ctl_stop( 'ice_istate_ini in iceistate: unable to allocate si structure' )   ;   RETURN
         ENDIF
         !
         DO ifpr = 1, jpfldi
            ALLOCATE( si(ifpr)%fnow(jpi,jpj,1) )
            IF( slf_i(ifpr)%ln_tint )  ALLOCATE( si(ifpr)%fdta(jpi,jpj,1,2) )
         END DO
         !
         ! fill si with slf_i and control print
         CALL fld_fill( si, slf_i, cn_dir, 'ice_istate_ini', 'initialization of sea ice fields', 'numnam_ice' )
         !
      ENDIF
      !
      IF( .NOT.ln_pnd ) THEN
         rn_apd_ini_n = 0. ; rn_apd_ini_s = 0.
         rn_hpd_ini_n = 0. ; rn_hpd_ini_s = 0.
         rn_hld_ini_n = 0. ; rn_hld_ini_s = 0.
         CALL ctl_warn( 'rn_apd_ini & rn_hpd_ini = 0 & rn_hld_ini = 0 when no ponds' )
      ENDIF
      !
      IF( .NOT.ln_pnd_lids ) THEN
         rn_hld_ini_n = 0. ; rn_hld_ini_s = 0.
      ENDIF
      !
   END SUBROUTINE ice_istate_init

#else
   !!----------------------------------------------------------------------
   !!   Default option :         Empty module         NO SI3 sea-ice model
   !!----------------------------------------------------------------------
#endif

   !!======================================================================
END MODULE iceistate