Newer
Older
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
MODULE icethd_do
!!======================================================================
!! *** MODULE icethd_do ***
!! sea-ice: sea ice growth in the leads (open water)
!!======================================================================
!! History : ! 2005-12 (M. Vancoppenolle) Original code
!! 4.0 ! 2018 (many people) SI3 [aka Sea Ice cube]
!!----------------------------------------------------------------------
#if defined key_si3
!!----------------------------------------------------------------------
!! 'key_si3' SI3 sea-ice model
!!----------------------------------------------------------------------
!! ice_thd_do : ice growth in open water (=lateral accretion of ice)
!! ice_thd_do_init : initialization
!!----------------------------------------------------------------------
USE dom_oce ! ocean space and time domain
USE phycst ! physical constants
USE sbc_oce , ONLY : sss_m
USE sbc_ice , ONLY : utau_ice, vtau_ice
USE ice1D ! sea-ice: thermodynamics variables
USE ice ! sea-ice: variables
USE icetab ! sea-ice: 2D <==> 1D
USE icectl ! sea-ice: conservation
USE icevar ! sea-ice: operations
USE icethd_sal ! sea-ice: salinity profiles
!
USE in_out_manager ! I/O manager
USE lib_mpp ! MPP library
USE lib_fortran ! fortran utilities (glob_sum + no signed zero)
IMPLICIT NONE
PRIVATE
PUBLIC ice_thd_do ! called by ice_thd
PUBLIC ice_thd_frazil ! called by ice_thd
PUBLIC ice_thd_do_init ! called by ice_stp
!! * Substitutions
# include "do_loop_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/ICE 4.0 , NEMO Consortium (2018)
!! $Id: icethd_do.F90 15388 2021-10-17 11:33:47Z clem $
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
SUBROUTINE ice_thd_do
!!-------------------------------------------------------------------
!! *** ROUTINE ice_thd_do ***
!!
!! ** Purpose : Computation of the evolution of the ice thickness and
!! concentration as a function of the heat balance in the leads
!!
!! ** Method : Ice is formed in the open water when ocean looses heat
!! (heat budget of open water is negative) following
!!
!! (dA/dt)acc = F[ (1-A)/(1-a) ] * [ Bl / (Li*h0) ]
!! where - h0 is the thickness of ice created in the lead
!! - a is a minimum fraction for leads
!! - F is a monotonic non-increasing function defined as:
!! F(X)=( 1 - X**exld )**(1.0/exld)
!! - exld is the exponent closure rate (=2 default val.)
!!
!! ** Action : - Adjustment of snow and ice thicknesses and heat
!! content in brine pockets
!! - Updating ice internal temperature
!! - Computation of variation of ice volume and mass
!! - Computation of a_i after lateral accretion and
!! update h_s_1d, h_i_1d
!!------------------------------------------------------------------------
INTEGER :: ji, jj, jk, jl ! dummy loop indices
!
REAL(wp) :: ztmelts
REAL(wp) :: zdE
REAL(wp) :: zQm ! enthalpy exchanged with the ocean (J/m2, >0 towards ocean)
REAL(wp) :: zEi ! sea ice specific enthalpy (J/kg)
REAL(wp) :: zEw ! seawater specific enthalpy (J/kg)
REAL(wp) :: zfmdt ! mass flux x time step (kg/m2, >0 towards ocean)
!
INTEGER :: jcat ! indexes of categories where new ice grows
REAL(wp) :: zswinew ! switch for new ice or not
!
REAL(wp) :: zv_newice ! volume of accreted ice
REAL(wp) :: za_newice ! fractional area of accreted ice
REAL(wp) :: ze_newice ! heat content of accreted ice
REAL(wp) :: zo_newice ! age of accreted ice
REAL(wp) :: zdv_res ! residual volume in case of excessive heat budget
REAL(wp) :: zda_res ! residual area in case of excessive heat budget
REAL(wp) :: zv_frazb ! accretion of frazil ice at the ice bottom
REAL(wp), DIMENSION(jpl) :: zv_b ! old volume of ice in category jl
REAL(wp), DIMENSION(jpl) :: za_b ! old area of ice in category jl
REAL(wp), DIMENSION(jpij) :: zs_newice ! salinity of accreted ice
REAL(wp), DIMENSION(jpij) :: zh_newice ! thickness of accreted ice
REAL(wp), DIMENSION(jpij) :: zfraz_frac_1d ! relative ice / frazil velocity (1D vector)
!
REAL(wp), DIMENSION(0:nlay_i+1) :: zh_i_old, ze_i_old
!!-----------------------------------------------------------------------!
IF( ln_icediachk ) CALL ice_cons_hsm( 0, 'icethd_do', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft )
IF( ln_icediachk ) CALL ice_cons2D ( 0, 'icethd_do', diag_v, diag_s, diag_t, diag_fv, diag_fs, diag_ft )
!------------------------------------------------------------------------------!
! 1) Compute thickness, salinity, enthalpy, age, area and volume of new ice
!------------------------------------------------------------------------------!
! it occurs if cooling
! Identify grid points where new ice forms
npti = 0 ; nptidx(:) = 0
IF ( qlead(ji,jj) < 0._wp ) THEN
npti = npti + 1
nptidx( npti ) = (jj - 1) * jpi + ji
ENDIF
END_2D
! Move from 2-D to 1-D vectors
IF ( npti > 0 ) THEN
CALL tab_2d_1d( npti, nptidx(1:npti), at_i_1d(1:npti) , at_i )
CALL tab_3d_2d( npti, nptidx(1:npti), a_i_2d (1:npti,1:jpl), a_i (:,:,:) )
CALL tab_3d_2d( npti, nptidx(1:npti), v_i_2d (1:npti,1:jpl), v_i (:,:,:) )
CALL tab_3d_2d( npti, nptidx(1:npti), sv_i_2d(1:npti,1:jpl), sv_i(:,:,:) )
CALL tab_4d_3d( npti, nptidx(1:npti), e_i_2d (1:npti,:,:) , e_i )
CALL tab_2d_1d( npti, nptidx(1:npti), qlead_1d (1:npti) , qlead )
CALL tab_2d_1d( npti, nptidx(1:npti), t_bo_1d (1:npti) , t_bo )
CALL tab_2d_1d( npti, nptidx(1:npti), sfx_opw_1d (1:npti) , sfx_opw )
CALL tab_2d_1d( npti, nptidx(1:npti), wfx_opw_1d (1:npti) , wfx_opw )
CALL tab_2d_1d( npti, nptidx(1:npti), zh_newice (1:npti) , ht_i_new )
CALL tab_2d_1d( npti, nptidx(1:npti), zfraz_frac_1d(1:npti) , fraz_frac )
CALL tab_2d_1d( npti, nptidx(1:npti), hfx_thd_1d(1:npti) , hfx_thd )
CALL tab_2d_1d( npti, nptidx(1:npti), hfx_opw_1d(1:npti) , hfx_opw )
CALL tab_2d_1d( npti, nptidx(1:npti), rn_amax_1d(1:npti) , rn_amax_2d )
CALL tab_2d_1d( npti, nptidx(1:npti), sss_1d (1:npti) , sss_m )
! Convert units for ice internal energy
DO jl = 1, jpl
DO jk = 1, nlay_i
WHERE( v_i_2d(1:npti,jl) > 0._wp )
e_i_2d(1:npti,jk,jl) = e_i_2d(1:npti,jk,jl) / v_i_2d(1:npti,jl) * REAL( nlay_i )
e_i_2d(1:npti,jk,jl) = 0._wp
END WHERE
END DO
END DO
! --- Salinity of new ice --- !
SELECT CASE ( nn_icesal )
CASE ( 1 ) ! Sice = constant
zs_newice(1:npti) = rn_icesal
CASE ( 2 ) ! Sice = F(z,t) [Vancoppenolle et al (2005)]
DO ji = 1, npti
zs_newice(ji) = MIN( 4.606 + 0.91 / zh_newice(ji) , rn_simax , 0.5 * sss_1d(ji) )
END DO
CASE ( 3 ) ! Sice = F(z) [multiyear ice]
zs_newice(1:npti) = 2.3
END SELECT
! ! ==================== !
! ! Start main loop here !
! ! ==================== !
! Keep old ice areas and volume in memory
DO jl = 1, jpl
zv_b(jl) = v_i_2d(ji,jl)
za_b(jl) = a_i_2d(ji,jl)
ENDDO
! --- Heat content of new ice --- !
! We assume that new ice is formed at the seawater freezing point
ztmelts = - rTmlt * zs_newice(ji) ! Melting point (C)
ze_newice = rhoi * ( rcpi * ( ztmelts - ( t_bo_1d(ji) - rt0 ) ) &
& + rLfus * ( 1.0 - ztmelts / MIN( t_bo_1d(ji) - rt0, -epsi10 ) ) &
& - rcp * ztmelts )
! --- Age of new ice --- !
zo_newice = 0._wp
! --- Volume of new ice --- !
zEi = - ze_newice * r1_rhoi ! specific enthalpy of forming ice [J/kg]
zEw = rcp * ( t_bo_1d(ji) - rt0 ) ! specific enthalpy of seawater at t_bo_1d [J/kg]
! clem: we suppose we are already at the freezing point (condition qlead<0 is satisfyied)
zdE = zEi - zEw ! specific enthalpy difference [J/kg]
zfmdt = - qlead_1d(ji) / zdE ! Fm.dt [kg/m2] (<0)
! clem: we use qlead instead of zqld (icethd) because we suppose we are at the freezing point
zQm = zfmdt * zEw ! heat to the ocean >0 associated with mass flux
! Contribution to heat flux to the ocean [W.m-2], >0
hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * zEw * r1_Dt_ice
! Total heat flux used in this process [W.m-2]
hfx_opw_1d(ji) = hfx_opw_1d(ji) - zfmdt * zdE * r1_Dt_ice
! mass flux
wfx_opw_1d(ji) = wfx_opw_1d(ji) - zv_newice * rhoi * r1_Dt_ice
sfx_opw_1d(ji) = sfx_opw_1d(ji) - zv_newice * rhoi * zs_newice(ji) * r1_Dt_ice
! A fraction fraz_frac of frazil ice is accreted at the ice bottom
rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp , - at_i_1d(ji) ) )
zv_frazb = zfraz_frac_1d(ji) * rswitch * zv_newice
zv_newice = ( 1._wp - zfraz_frac_1d(ji) * rswitch ) * zv_newice
! --- Area of new ice --- !
za_newice = zv_newice / zh_newice(ji)
! --- Redistribute new ice area and volume into ice categories --- !
! --- lateral ice growth --- !
! If lateral ice growth gives an ice concentration > amax, then
! we keep the excessive volume in memory and attribute it later to bottom accretion
IF ( za_newice > MAX( 0._wp, rn_amax_1d(ji) - at_i_1d(ji) ) ) THEN ! max is for roundoff error
zda_res = za_newice - MAX( 0._wp, rn_amax_1d(ji) - at_i_1d(ji) )
zdv_res = zda_res * zh_newice(ji)
za_newice = MAX( 0._wp, za_newice - zda_res )
zv_newice = MAX( 0._wp, zv_newice - zdv_res )
! find which category to fill
at_i_1d(ji) = 0._wp
DO jl = 1, jpl
IF( zh_newice(ji) > hi_max(jl-1) .AND. zh_newice(ji) <= hi_max(jl) ) THEN
a_i_2d(ji,jl) = a_i_2d(ji,jl) + za_newice
v_i_2d(ji,jl) = v_i_2d(ji,jl) + zv_newice
jcat = jl
at_i_1d(ji) = at_i_1d(ji) + a_i_2d(ji,jl)
! Heat content
jl = jcat ! categroy in which new ice is put
zswinew = MAX( 0._wp , SIGN( 1._wp , - za_b(jl) ) ) ! 0 if old ice
rswitch = MAX( 0._wp, SIGN( 1._wp , v_i_2d(ji,jl) - epsi20 ) )
e_i_2d(ji,jk,jl) = zswinew * ze_newice + &
& ( 1.0 - zswinew ) * ( ze_newice * zv_newice + e_i_2d(ji,jk,jl) * zv_b(jl) ) &
! --- bottom ice growth + ice enthalpy remapping --- !
DO jl = 1, jpl
! for remapping
zh_i_old(0:nlay_i+1) = 0._wp
ze_i_old(0:nlay_i+1) = 0._wp
zh_i_old(jk) = v_i_2d(ji,jl) * r1_nlay_i
ze_i_old(jk) = e_i_2d(ji,jk,jl) * v_i_2d(ji,jl) * r1_nlay_i
! new volumes including lateral/bottom accretion + residual
rswitch = MAX( 0._wp, SIGN( 1._wp , at_i_1d(ji) - epsi20 ) )
zv_newfra = rswitch * ( zdv_res + zv_frazb ) * a_i_2d(ji,jl) / MAX( at_i_1d(ji) , epsi20 )
a_i_2d(ji,jl) = rswitch * a_i_2d(ji,jl)
v_i_2d(ji,jl) = v_i_2d(ji,jl) + zv_newfra
! for remapping
zh_i_old(nlay_i+1) = zv_newfra
ze_i_old(nlay_i+1) = ze_newice * zv_newfra
sv_i_2d(ji,jl) = sv_i_2d(ji,jl) + zs_newice(ji) * ( v_i_2d(ji,jl) - zv_b(jl) )
e_i_2d(ji,:,jl) = ice_ent2( zh_i_old(:), ze_i_old(:) )
! ! ================== !
! ! End main loop here !
! ! ================== !
!
! Change units for e_i
DO jl = 1, jpl
DO jk = 1, nlay_i
e_i_2d(1:npti,jk,jl) = e_i_2d(1:npti,jk,jl) * v_i_2d(1:npti,jl) * r1_nlay_i
END DO
END DO
! Move 2D vectors to 1D vectors
CALL tab_2d_3d( npti, nptidx(1:npti), a_i_2d (1:npti,1:jpl), a_i (:,:,:) )
CALL tab_2d_3d( npti, nptidx(1:npti), v_i_2d (1:npti,1:jpl), v_i (:,:,:) )
CALL tab_2d_3d( npti, nptidx(1:npti), sv_i_2d(1:npti,1:jpl), sv_i(:,:,:) )
CALL tab_3d_4d( npti, nptidx(1:npti), e_i_2d (1:npti,:,:) , e_i )
CALL tab_1d_2d( npti, nptidx(1:npti), sfx_opw_1d(1:npti), sfx_opw )
CALL tab_1d_2d( npti, nptidx(1:npti), wfx_opw_1d(1:npti), wfx_opw )
CALL tab_1d_2d( npti, nptidx(1:npti), hfx_thd_1d(1:npti), hfx_thd )
CALL tab_1d_2d( npti, nptidx(1:npti), hfx_opw_1d(1:npti), hfx_opw )
!
ENDIF ! npti > 0
!
! the following fields need to be updated on the halos (done in icethd): a_i, v_i, sv_i, e_i
!
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
IF( ln_icediachk ) CALL ice_cons_hsm(1, 'icethd_do', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft)
IF( ln_icediachk ) CALL ice_cons2D (1, 'icethd_do', diag_v, diag_s, diag_t, diag_fv, diag_fs, diag_ft)
!
END SUBROUTINE ice_thd_do
SUBROUTINE ice_thd_frazil
!!-----------------------------------------------------------------------
!! *** ROUTINE ice_thd_frazil ***
!!
!! ** Purpose : frazil ice collection thickness and fraction
!!
!! ** Inputs : u_ice, v_ice, utau_ice, vtau_ice
!! ** Ouputs : ht_i_new, fraz_frac
!!-----------------------------------------------------------------------
INTEGER :: ji, jj ! dummy loop indices
INTEGER :: iter
REAL(wp) :: zvfrx, zvgx, ztaux, zf, ztenagm, zvfry, zvgy, ztauy, zvrel2, zfp, ztwogp
REAL(wp), PARAMETER :: zcai = 1.4e-3_wp ! ice-air drag (clem: should be dependent on coupling/forcing used)
REAL(wp), PARAMETER :: zhicrit = 0.04_wp ! frazil ice thickness
REAL(wp), PARAMETER :: zsqcd = 1.0_wp / SQRT( 1.3_wp * zcai ) ! 1/SQRT(airdensity*drag)
REAL(wp), PARAMETER :: zgamafr = 0.03_wp
!!-----------------------------------------------------------------------
!
!---------------------------------------------------------!
! Collection thickness of ice formed in leads and polynyas
!---------------------------------------------------------!
! ht_i_new is the thickness of new ice formed in open water
! ht_i_new can be either prescribed (ln_frazil=F) or computed (ln_frazil=T)
! Frazil ice forms in open water, is transported by wind, accumulates at the edge of the consolidated ice edge
! where it forms aggregates of a specific thickness called collection thickness.
!
fraz_frac(:,:) = 0._wp
!
! Default new ice thickness
WHERE( qlead(:,:) < 0._wp ) ! cooling
ht_i_new(:,:) = rn_hinew
ELSEWHERE
ht_i_new(:,:) = 0._wp
END WHERE
IF( ln_frazil ) THEN
ztwogp = 2._wp * rho0 / ( grav * 0.3_wp * ( rho0 - rhoi ) ) ! reduced grav
!
DO_2D( 0, 0, 0, 0 )
IF ( qlead(ji,jj) < 0._wp ) THEN ! cooling
! -- Wind stress -- !
ztaux = utau_ice(ji,jj) * smask0(ji,jj)
ztauy = vtau_ice(ji,jj) * smask0(ji,jj)
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
! Square root of wind stress
ztenagm = SQRT( SQRT( ztaux * ztaux + ztauy * ztauy ) )
! -- Frazil ice velocity -- !
rswitch = MAX( 0._wp, SIGN( 1._wp , ztenagm - epsi10 ) )
zvfrx = rswitch * zgamafr * zsqcd * ztaux / MAX( ztenagm, epsi10 )
zvfry = rswitch * zgamafr * zsqcd * ztauy / MAX( ztenagm, epsi10 )
! -- Pack ice velocity -- !
zvgx = ( u_ice(ji-1,jj ) * umask(ji-1,jj ,1) + u_ice(ji,jj) * umask(ji,jj,1) ) * 0.5_wp
zvgy = ( v_ice(ji ,jj-1) * vmask(ji ,jj-1,1) + v_ice(ji,jj) * vmask(ji,jj,1) ) * 0.5_wp
! -- Relative frazil/pack ice velocity -- !
rswitch = MAX( 0._wp, SIGN( 1._wp , at_i(ji,jj) - epsi10 ) )
zvrel2 = MAX( (zvfrx - zvgx)*(zvfrx - zvgx) + (zvfry - zvgy)*(zvfry - zvgy), 0.15_wp*0.15_wp ) * rswitch
! -- fraction of frazil ice -- !
fraz_frac(ji,jj) = rswitch * ( TANH( rn_Cfraz * ( SQRT(zvrel2) - rn_vfraz ) ) + 1._wp ) * 0.5_wp * rn_maxfraz
! -- new ice thickness (iterative loop) -- !
ht_i_new(ji,jj) = zhicrit + ( zhicrit + 0.1_wp ) &
& / ( ( zhicrit + 0.1_wp ) * ( zhicrit + 0.1_wp ) - zhicrit * zhicrit ) * ztwogp * zvrel2
iter = 1
DO WHILE ( iter < 20 )
zf = ( ht_i_new(ji,jj) - zhicrit ) * ( ht_i_new(ji,jj) * ht_i_new(ji,jj) - zhicrit * zhicrit ) - &
& ht_i_new(ji,jj) * zhicrit * ztwogp * zvrel2
zfp = ( ht_i_new(ji,jj) - zhicrit ) * ( 3.0_wp * ht_i_new(ji,jj) + zhicrit ) - zhicrit * ztwogp * zvrel2
ht_i_new(ji,jj) = ht_i_new(ji,jj) - zf / MAX( zfp, epsi20 )
iter = iter + 1
END DO
!
! bound ht_i_new (though I don't see why it should be necessary)
ht_i_new(ji,jj) = MAX( 0.01_wp, MIN( ht_i_new(ji,jj), hi_max(jpl) ) )
!
ELSE
ht_i_new(ji,jj) = 0._wp
ENDIF
!
END_2D
!
ENDIF
END SUBROUTINE ice_thd_frazil
!!-------------------------------------------------------------------
!!
!! ** Purpose :
!! This routine computes new vertical grids in the ice,
!! and consistently redistributes temperatures.
!! Redistribution is made so as to ensure to energy conservation
!!
!!
!! ** Method : linear conservative remapping
!!
!! ** Steps : 1) cumulative integrals of old enthalpies/thicknesses
!! 2) linear remapping on the new layers
!!
!! ------------ cum0(0) ------------- cum1(0)
!! NEW -------------
!! ------------ cum0(1) ==> -------------
!! ... -------------
!! ------------ -------------
!! ------------ cum0(nlay_i+2) ------------- cum1(nlay_i)
!!
!!
!! References : Bitz & Lipscomb, JGR 99; Vancoppenolle et al., GRL, 2005
!!-------------------------------------------------------------------
REAL(wp), DIMENSION(0:nlay_i+1), INTENT(in) :: ph_old, pe_old ! old tickness and enthlapy
REAL(wp), DIMENSION(1:nlay_i) :: ice_ent2 ! new enthlapies (J.m-3, remapped)
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
!
INTEGER :: ji ! dummy loop indices
INTEGER :: jk0, jk1 ! old/new layer indices
REAL(wp) :: zswitch
!
REAL(wp), DIMENSION(0:nlay_i+2) :: zeh_cum0, zh_cum0 ! old cumulative enthlapies and layers interfaces
REAL(wp), DIMENSION(0:nlay_i) :: zeh_cum1, zh_cum1 ! new cumulative enthlapies and layers interfaces
REAL(wp) :: zhnew ! new layers thicknesses
!!-------------------------------------------------------------------
!--------------------------------------------------------------------------
! 1) Cumulative integral of old enthalpy * thickness and layers interfaces
!--------------------------------------------------------------------------
zeh_cum0(0) = 0._wp
zh_cum0 (0) = 0._wp
DO jk0 = 1, nlay_i+2
zeh_cum0(jk0) = zeh_cum0(jk0-1) + pe_old(jk0-1)
zh_cum0 (jk0) = zh_cum0 (jk0-1) + ph_old(jk0-1)
END DO
!------------------------------------
! 2) Interpolation on the new layers
!------------------------------------
! new layer thickesses
zhnew = SUM( ph_old(0:nlay_i+1) ) * r1_nlay_i
! new layers interfaces
zh_cum1(0) = 0._wp
DO jk1 = 1, nlay_i
zh_cum1(jk1) = zh_cum1(jk1-1) + zhnew
END DO
zeh_cum1(0:nlay_i) = 0._wp
! new cumulative q*h => linear interpolation
DO jk0 = 1, nlay_i+2
DO jk1 = 1, nlay_i-1
IF( zh_cum1(jk1) <= zh_cum0(jk0) .AND. zh_cum1(jk1) > zh_cum0(jk0-1) ) THEN
zeh_cum1(jk1) = ( zeh_cum0(jk0-1) * ( zh_cum0(jk0) - zh_cum1(jk1 ) ) + &
& zeh_cum0(jk0 ) * ( zh_cum1(jk1) - zh_cum0(jk0-1) ) ) &
& / ( zh_cum0(jk0) - zh_cum0(jk0-1) )
ENDIF
END DO
END DO
! to ensure that total heat content is strictly conserved, set:
zeh_cum1(nlay_i) = zeh_cum0(nlay_i+2)
! new enthalpies
DO jk1 = 1, nlay_i
zswitch = MAX( 0._wp , SIGN( 1._wp , zhnew - epsi20 ) )
ice_ent2(jk1) = zswitch * MAX( 0._wp, zeh_cum1(jk1) - zeh_cum1(jk1-1) ) / MAX( zhnew, epsi20 ) ! max for roundoff error
END DO
! --- diag error on heat remapping --- !
! comment: if input h_old and eh_old are already multiplied by a_i (as in icethd_do),
! then we should not (* a_i) again but not important since this is just to check that remap error is ~0
! hfx_err_rem_1d(ji) = hfx_err_rem_1d(ji) + a_i_1d(ji) * r1_Dt_ice * &
! & ( SUM( pe_new(ji,1:nlay_i) ) * zhnew(ji) - SUM( eh_old(ji,0:nlay_i+1) ) )
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
SUBROUTINE ice_thd_do_init
!!-----------------------------------------------------------------------
!! *** ROUTINE ice_thd_do_init ***
!!
!! ** Purpose : Physical constants and parameters associated with
!! ice growth in the leads
!!
!! ** Method : Read the namthd_do namelist and check the parameters
!! called at the first timestep (nit000)
!!
!! ** input : Namelist namthd_do
!!-------------------------------------------------------------------
INTEGER :: ios ! Local integer
!!
NAMELIST/namthd_do/ rn_hinew, ln_frazil, rn_maxfraz, rn_vfraz, rn_Cfraz
!!-------------------------------------------------------------------
!
READ ( numnam_ice_ref, namthd_do, IOSTAT = ios, ERR = 901)
901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namthd_do in reference namelist' )
READ ( numnam_ice_cfg, namthd_do, IOSTAT = ios, ERR = 902 )
902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namthd_do in configuration namelist' )
IF(lwm) WRITE( numoni, namthd_do )
!
IF(lwp) THEN ! control print
WRITE(numout,*)
WRITE(numout,*) 'ice_thd_do_init: Ice growth in open water'
WRITE(numout,*) '~~~~~~~~~~~~~~~'
WRITE(numout,*) ' Namelist namthd_do:'
WRITE(numout,*) ' ice thickness for lateral accretion rn_hinew = ', rn_hinew
WRITE(numout,*) ' Frazil ice thickness as a function of wind or not ln_frazil = ', ln_frazil
WRITE(numout,*) ' Maximum proportion of frazil ice collecting at bottom rn_maxfraz = ', rn_maxfraz
WRITE(numout,*) ' Threshold relative drift speed for collection of frazil rn_vfraz = ', rn_vfraz
WRITE(numout,*) ' Squeezing coefficient for collection of frazil rn_Cfraz = ', rn_Cfraz
ENDIF
!
IF ( rn_hinew < rn_himin ) CALL ctl_stop( 'ice_thd_do_init : rn_hinew should be >= rn_himin' )
!
END SUBROUTINE ice_thd_do_init
#else
!!----------------------------------------------------------------------
!! Default option NO SI3 sea-ice model
!!----------------------------------------------------------------------
#endif
!!======================================================================
END MODULE icethd_do