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
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
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
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 icethd_ent ! sea-ice: thermodynamics, enthalpy
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)
USE lbclnk ! lateral boundary conditions (or mpp links)
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)
!
REAL(wp) :: zv_newfra
!
INTEGER , DIMENSION(jpij) :: jcat ! indexes of categories where new ice grows
REAL(wp), DIMENSION(jpij) :: zswinew ! switch for new ice or not
!
REAL(wp), DIMENSION(jpij) :: zv_newice ! volume of accreted ice
REAL(wp), DIMENSION(jpij) :: za_newice ! fractional area of accreted ice
REAL(wp), DIMENSION(jpij) :: zh_newice ! thickness of accreted ice
REAL(wp), DIMENSION(jpij) :: ze_newice ! heat content of accreted ice
REAL(wp), DIMENSION(jpij) :: zs_newice ! salinity of accreted ice
REAL(wp), DIMENSION(jpij) :: zo_newice ! age of accreted ice
REAL(wp), DIMENSION(jpij) :: zdv_res ! residual volume in case of excessive heat budget
REAL(wp), DIMENSION(jpij) :: zda_res ! residual area in case of excessive heat budget
REAL(wp), DIMENSION(jpij) :: zv_frazb ! accretion of frazil ice at the ice bottom
REAL(wp), DIMENSION(jpij) :: zfraz_frac_1d ! relative ice / frazil velocity (1D vector)
!
REAL(wp), DIMENSION(jpij,jpl) :: zv_b ! old volume of ice in category jl
REAL(wp), DIMENSION(jpij,jpl) :: za_b ! old area of ice in category jl
!
REAL(wp), DIMENSION(jpij,nlay_i,jpl) :: ze_i_2d !: 1-D version of e_i
!
!!-----------------------------------------------------------------------!
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 )
at_i(:,:) = SUM( a_i, dim=3 )
!------------------------------------------------------------------------------!
! 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
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
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(:,:,:) )
DO jl = 1, jpl
DO jk = 1, nlay_i
CALL tab_2d_1d( npti, nptidx(1:npti), ze_i_2d(1:npti,jk,jl), e_i(:,:,jk,jl) )
END DO
END DO
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 )
ze_i_2d(1:npti,jk,jl) = ze_i_2d(1:npti,jk,jl) / v_i_2d(1:npti,jl) * REAL( nlay_i )
ELSEWHERE
ze_i_2d(1:npti,jk,jl) = 0._wp
END WHERE
END DO
END DO
! Keep old ice areas and volume in memory
zv_b(1:npti,:) = v_i_2d(1:npti,:)
za_b(1:npti,:) = a_i_2d(1:npti,:)
! --- 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
! --- Heat content of new ice --- !
! We assume that new ice is formed at the seawater freezing point
DO ji = 1, npti
ztmelts = - rTmlt * zs_newice(ji) ! Melting point (C)
ze_newice(ji) = rhoi * ( rcpi * ( ztmelts - ( t_bo_1d(ji) - rt0 ) ) &
& + rLfus * ( 1.0 - ztmelts / MIN( t_bo_1d(ji) - rt0, -epsi10 ) ) &
& - rcp * ztmelts )
END DO
! --- Age of new ice --- !
zo_newice(1:npti) = 0._wp
! --- Volume of new ice --- !
DO ji = 1, npti
zEi = - ze_newice(ji) * 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
zv_newice(ji) = - zfmdt * r1_rhoi
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(ji) * rhoi * r1_Dt_ice
! salt flux
sfx_opw_1d(ji) = sfx_opw_1d(ji) - zv_newice(ji) * rhoi * zs_newice(ji) * r1_Dt_ice
END DO
! A fraction fraz_frac of frazil ice is accreted at the ice bottom
DO ji = 1, npti
rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp , - at_i_1d(ji) ) )
zv_frazb(ji) = zfraz_frac_1d(ji) * rswitch * zv_newice(ji)
zv_newice(ji) = ( 1._wp - zfraz_frac_1d(ji) * rswitch ) * zv_newice(ji)
END DO
! --- Area of new ice --- !
DO ji = 1, npti
za_newice(ji) = zv_newice(ji) / zh_newice(ji)
END DO
!------------------------------------------------------------------------------!
! 2) 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
DO ji = 1, npti
IF ( za_newice(ji) > MAX( 0._wp, rn_amax_1d(ji) - at_i_1d(ji) ) ) THEN ! max is for roundoff error
zda_res(ji) = za_newice(ji) - MAX( 0._wp, rn_amax_1d(ji) - at_i_1d(ji) )
zdv_res(ji) = zda_res (ji) * zh_newice(ji)
za_newice(ji) = MAX( 0._wp, za_newice(ji) - zda_res (ji) )
zv_newice(ji) = MAX( 0._wp, zv_newice(ji) - zdv_res (ji) )
ELSE
zda_res(ji) = 0._wp
zdv_res(ji) = 0._wp
ENDIF
END DO
! find which category to fill
DO jl = 1, jpl
DO ji = 1, npti
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(ji)
v_i_2d(ji,jl) = v_i_2d(ji,jl) + zv_newice(ji)
jcat(ji) = jl
ENDIF
END DO
END DO
at_i_1d(1:npti) = SUM( a_i_2d(1:npti,:), dim=2 )
! Heat content
DO ji = 1, npti
jl = jcat(ji) ! categroy in which new ice is put
zswinew (ji) = MAX( 0._wp , SIGN( 1._wp , - za_b(ji,jl) ) ) ! 0 if old ice
END DO
DO jk = 1, nlay_i
DO ji = 1, npti
jl = jcat(ji)
rswitch = MAX( 0._wp, SIGN( 1._wp , v_i_2d(ji,jl) - epsi20 ) )
ze_i_2d(ji,jk,jl) = zswinew(ji) * ze_newice(ji) + &
& ( 1.0 - zswinew(ji) ) * ( ze_newice(ji) * zv_newice(ji) + ze_i_2d(ji,jk,jl) * zv_b(ji,jl) ) &
& * rswitch / MAX( v_i_2d(ji,jl), epsi20 )
END DO
END DO
! --- bottom ice growth + ice enthalpy remapping --- !
DO jl = 1, jpl
! for remapping
h_i_old (1:npti,0:nlay_i+1) = 0._wp
eh_i_old(1:npti,0:nlay_i+1) = 0._wp
DO jk = 1, nlay_i
DO ji = 1, npti
h_i_old (ji,jk) = v_i_2d(ji,jl) * r1_nlay_i
eh_i_old(ji,jk) = ze_i_2d(ji,jk,jl) * h_i_old(ji,jk)
END DO
END DO
! new volumes including lateral/bottom accretion + residual
DO ji = 1, npti
rswitch = MAX( 0._wp, SIGN( 1._wp , at_i_1d(ji) - epsi20 ) )
zv_newfra = rswitch * ( zdv_res(ji) + zv_frazb(ji) ) * 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
h_i_old (ji,nlay_i+1) = zv_newfra
eh_i_old(ji,nlay_i+1) = ze_newice(ji) * zv_newfra
END DO
! --- Ice enthalpy remapping --- !
CALL ice_thd_ent( ze_i_2d(1:npti,:,jl) )
END DO
! --- Update salinity --- !
DO jl = 1, jpl
DO ji = 1, npti
sv_i_2d(ji,jl) = sv_i_2d(ji,jl) + zs_newice(ji) * ( v_i_2d(ji,jl) - zv_b(ji,jl) )
END DO
END DO
! Change units for e_i
DO jl = 1, jpl
DO jk = 1, nlay_i
ze_i_2d(1:npti,jk,jl) = ze_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(:,:,:) )
DO jl = 1, jpl
DO jk = 1, nlay_i
CALL tab_1d_2d( npti, nptidx(1:npti), ze_i_2d(1:npti,jk,jl), e_i(:,:,jk,jl) )
END DO
END DO
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
!
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-1,jj ) * umask(ji-1,jj ,1) + utau_ice(ji,jj) * umask(ji,jj,1) ) * 0.5_wp
ztauy = ( vtau_ice(ji ,jj-1) * vmask(ji ,jj-1,1) + vtau_ice(ji,jj) * vmask(ji,jj,1) ) * 0.5_wp
! 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
!
CALL lbc_lnk( 'icethd_frazil', fraz_frac, 'T', 1.0_wp, ht_i_new, 'T', 1.0_wp )
ENDIF
END SUBROUTINE ice_thd_frazil
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