New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
limthd.F90 in trunk/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90 @ 5134

Last change on this file since 5134 was 5134, checked in by clem, 9 years ago

LIM3: changes to constrain ice thickness after advection. Ice volume and concentration are advected while ice thickness is deduced. This could result in overly high thicknesses associated with very low concentrations (in the 5th category)
  • Property svn:keywords set to Id
File size: 39.7 KB
Line 
1MODULE limthd
2   !!======================================================================
3   !!                  ***  MODULE limthd   ***
4   !!  LIM-3 :   ice thermodynamic
5   !!======================================================================
6   !! History :  LIM  ! 2000-01 (M.A. Morales Maqueda, H. Goosse, T. Fichefet) LIM-1
7   !!            2.0  ! 2002-07 (C. Ethe, G. Madec)  LIM-2 (F90 rewriting)
8   !!            3.0  ! 2005-11 (M. Vancoppenolle)  LIM-3 : Multi-layer thermodynamics + salinity variations
9   !!             -   ! 2007-04 (M. Vancoppenolle) add lim_thd_glohec, lim_thd_con_dh and lim_thd_con_dif
10   !!            3.2  ! 2009-07 (M. Vancoppenolle, Y. Aksenov, G. Madec) bug correction in wfx_snw
11   !!            3.3  ! 2010-11 (G. Madec) corrected snow melting heat (due to factor betas)
12   !!            4.0  ! 2011-02 (G. Madec) dynamical allocation
13   !!             -   ! 2012-05 (C. Rousset) add penetration solar flux
14   !!----------------------------------------------------------------------
15#if defined key_lim3
16   !!----------------------------------------------------------------------
17   !!   'key_lim3'                                      LIM3 sea-ice model
18   !!----------------------------------------------------------------------
19   !!   lim_thd       : thermodynamic of sea ice
20   !!   lim_thd_init  : initialisation of sea-ice thermodynamic
21   !!----------------------------------------------------------------------
22   USE phycst         ! physical constants
23   USE dom_oce        ! ocean space and time domain variables
24   USE oce     , ONLY : fraqsr_1lev 
25   USE ice            ! LIM: sea-ice variables
26   USE sbc_oce        ! Surface boundary condition: ocean fields
27   USE sbc_ice        ! Surface boundary condition: ice fields
28   USE thd_ice        ! LIM thermodynamic sea-ice variables
29   USE dom_ice        ! LIM sea-ice domain
30   USE domvvl         ! domain: variable volume level
31   USE limthd_dif     ! LIM: thermodynamics, vertical diffusion
32   USE limthd_dh      ! LIM: thermodynamics, ice and snow thickness variation
33   USE limthd_sal     ! LIM: thermodynamics, ice salinity
34   USE limthd_ent     ! LIM: thermodynamics, ice enthalpy redistribution
35   USE limthd_lac     ! LIM-3 lateral accretion
36   USE limitd_th      ! remapping thickness distribution
37   USE limtab         ! LIM: 1D <==> 2D transformation
38   USE limvar         ! LIM: sea-ice variables
39   USE lbclnk         ! lateral boundary condition - MPP links
40   USE lib_mpp        ! MPP library
41   USE wrk_nemo       ! work arrays
42   USE in_out_manager ! I/O manager
43   USE prtctl         ! Print control
44   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
45   USE timing         ! Timing
46   USE limcons        ! conservation tests
47   USE limctl
48
49   IMPLICIT NONE
50   PRIVATE
51
52   PUBLIC   lim_thd        ! called by limstp module
53   PUBLIC   lim_thd_init   ! called by sbc_lim_init
54
55   !! * Substitutions
56#  include "domzgr_substitute.h90"
57#  include "vectopt_loop_substitute.h90"
58   !!----------------------------------------------------------------------
59   !! NEMO/LIM3 3.3 , UCL - NEMO Consortium (2010)
60   !! $Id$
61   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
62   !!----------------------------------------------------------------------
63CONTAINS
64
65   SUBROUTINE lim_thd( kt )
66      !!-------------------------------------------------------------------
67      !!                ***  ROUTINE lim_thd  ***       
68      !! 
69      !! ** Purpose : This routine manages ice thermodynamics
70      !!         
71      !! ** Action : - Initialisation of some variables
72      !!             - Some preliminary computation (oceanic heat flux
73      !!               at the ice base, snow acc.,heat budget of the leads)
74      !!             - selection of the icy points and put them in an array
75      !!             - call lim_thd_dif  for vertical heat diffusion
76      !!             - call lim_thd_dh   for vertical ice growth and melt
77      !!             - call lim_thd_ent  for enthalpy remapping
78      !!             - call lim_thd_sal  for ice desalination
79      !!             - call lim_thd_temp to  retrieve temperature from ice enthalpy
80      !!             - back to the geographic grid
81      !!     
82      !! ** References :
83      !!---------------------------------------------------------------------
84      INTEGER, INTENT(in) :: kt    ! number of iteration
85      !!
86      INTEGER  :: ji, jj, jk, jl   ! dummy loop indices
87      INTEGER  :: nbpb             ! nb of icy pts for vertical thermo calculations
88      INTEGER  :: ii, ij           ! temporary dummy loop index
89      REAL(wp) :: zfric_u, zqld, zqfr
90      REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 
91      REAL(wp), PARAMETER :: zfric_umin = 0._wp        ! lower bound for the friction velocity (cice value=5.e-04)
92      REAL(wp), PARAMETER :: zch        = 0.0057_wp    ! heat transfer coefficient
93      !
94      REAL(wp), POINTER, DIMENSION(:,:) ::  zqsr, zqns
95      !!-------------------------------------------------------------------
96      CALL wrk_alloc( jpi, jpj, zqsr, zqns )
97
98      IF( nn_timing == 1 )  CALL timing_start('limthd')
99
100      ! conservation test
101      IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limthd', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
102
103      !------------------------------------------------------------------------!
104      ! 1) Initialization of some variables                                    !
105      !------------------------------------------------------------------------!
106      ftr_ice(:,:,:) = 0._wp  ! part of solar radiation transmitted through the ice
107
108      !--------------------
109      ! 1.2) Heat content   
110      !--------------------
111      ! Change the units of heat content; from J/m2 to J/m3
112      DO jl = 1, jpl
113         DO jk = 1, nlay_i
114            DO jj = 1, jpj
115               DO ji = 1, jpi
116                  !0 if no ice and 1 if yes
117                  rswitch = MAX(  0._wp , SIGN( 1._wp , v_i(ji,jj,jl) - epsi20 )  )
118                  !Energy of melting q(S,T) [J.m-3]
119                  e_i(ji,jj,jk,jl) = rswitch * e_i(ji,jj,jk,jl) / MAX( v_i(ji,jj,jl) , epsi20 ) * REAL( nlay_i )
120               END DO
121            END DO
122         END DO
123         DO jk = 1, nlay_s
124            DO jj = 1, jpj
125               DO ji = 1, jpi
126                  !0 if no ice and 1 if yes
127                  rswitch = MAX(  0._wp , SIGN( 1._wp , v_s(ji,jj,jl) - epsi20 )  )
128                  !Energy of melting q(S,T) [J.m-3]
129                  e_s(ji,jj,jk,jl) = rswitch * e_s(ji,jj,jk,jl) / MAX( v_s(ji,jj,jl) , epsi20 ) * REAL( nlay_s )
130               END DO
131            END DO
132         END DO
133      END DO
134
135      ! 2) Partial computation of forcing for the thermodynamic sea ice model.      !
136      !-----------------------------------------------------------------------------!
137
138      !--- Ocean solar and non solar fluxes to be used in zqld
139      IF ( .NOT. lk_cpl ) THEN   ! --- forced case, fluxes to the lead are the same as over the ocean
140         !
141         zqsr(:,:) = qsr(:,:)      ; zqns(:,:) = qns(:,:)
142         !
143      ELSE                       ! --- coupled case, fluxes to the lead are total - intercepted
144         !
145         zqsr(:,:) = qsr_tot(:,:)  ; zqns(:,:) = qns_tot(:,:)
146         !
147         DO jl = 1, jpl
148            DO jj = 1, jpj
149               DO ji = 1, jpi
150                  zqsr(ji,jj) = zqsr(ji,jj) - qsr_ice(ji,jj,jl) * a_i_b(ji,jj,jl)
151                  zqns(ji,jj) = zqns(ji,jj) - qns_ice(ji,jj,jl) * a_i_b(ji,jj,jl)
152               END DO
153            END DO
154         END DO
155         !
156      ENDIF
157
158      DO jj = 1, jpj
159         DO ji = 1, jpi
160            rswitch  = tmask(ji,jj,1) * MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) ! 0 if no ice
161            !
162            !           !  solar irradiance transmission at the mixed layer bottom and used in the lead heat budget
163            !           !  practically no "direct lateral ablation"
164            !           
165            !           !  net downward heat flux from the ice to the ocean, expressed as a function of ocean
166            !           !  temperature and turbulent mixing (McPhee, 1992)
167            !
168
169            ! --- Energy received in the lead, zqld is defined everywhere (J.m-2) --- !
170            ! REMARK valid at least in forced mode from clem
171            ! precip is included in qns but not in qns_ice
172            IF ( lk_cpl ) THEN
173               zqld =  tmask(ji,jj,1) * rdt_ice *  &
174                  &    (   zqsr(ji,jj) * fraqsr_1lev(ji,jj) + zqns(ji,jj)               &   ! pfrld already included in coupled mode
175                  &    + ( pfrld(ji,jj)**rn_betas - pfrld(ji,jj) ) * sprecip(ji,jj)  *     &   ! heat content of precip
176                  &      ( cpic * ( MIN( tatm_ice(ji,jj), rt0_snow ) - rt0 ) - lfus )   &
177                  &    + ( 1._wp - pfrld(ji,jj) ) * ( tprecip(ji,jj) - sprecip(ji,jj) ) * rcp * ( tatm_ice(ji,jj) - rt0 ) )
178            ELSE
179               zqld =  tmask(ji,jj,1) * rdt_ice *  &
180                  &      ( pfrld(ji,jj) * ( zqsr(ji,jj) * fraqsr_1lev(ji,jj) + zqns(ji,jj) )    &
181                  &    + ( pfrld(ji,jj)**rn_betas - pfrld(ji,jj) ) * sprecip(ji,jj)  *             &  ! heat content of precip
182                  &      ( cpic * ( MIN( tatm_ice(ji,jj), rt0_snow ) - rt0 ) - lfus )           &
183                  &    + ( 1._wp - pfrld(ji,jj) ) * ( tprecip(ji,jj) - sprecip(ji,jj) ) * rcp * ( tatm_ice(ji,jj) - rt0 ) )
184            ENDIF
185
186            ! --- Energy needed to bring ocean surface layer until its freezing (<0, J.m-2) --- !
187            zqfr = tmask(ji,jj,1) * rau0 * rcp * fse3t_m(ji,jj) * ( t_bo(ji,jj) - ( sst_m(ji,jj) + rt0 ) )
188
189            ! --- Energy from the turbulent oceanic heat flux (W/m2) --- !
190            zfric_u      = MAX( SQRT( ust2s(ji,jj) ), zfric_umin ) 
191            fhtur(ji,jj) = MAX( 0._wp, rswitch * rau0 * rcp * zch  * zfric_u * ( ( sst_m(ji,jj) + rt0 ) - t_bo(ji,jj) ) ) ! W.m-2
192            fhtur(ji,jj) = rswitch * MIN( fhtur(ji,jj), - zqfr * r1_rdtice / MAX( at_i(ji,jj), epsi10 ) )
193            ! upper bound for fhtur: the heat retrieved from the ocean must be smaller than the heat necessary to reach
194            !                        the freezing point, so that we do not have SST < T_freeze
195            !                        This implies: - ( fhtur(ji,jj) * at_i(ji,jj) * rtdice ) - zqfr >= 0
196
197            !-- Energy Budget of the leads (J.m-2). Must be < 0 to form ice
198            qlead(ji,jj) = MIN( 0._wp , zqld - ( fhtur(ji,jj) * at_i(ji,jj) * rdt_ice ) - zqfr )
199
200            ! If there is ice and leads are warming, then transfer energy from the lead budget and use it for bottom melting
201            IF( zqld > 0._wp ) THEN
202               fhld (ji,jj) = rswitch * zqld * r1_rdtice / MAX( at_i(ji,jj), epsi10 ) ! divided by at_i since this is (re)multiplied by a_i in limthd_dh.F90
203               qlead(ji,jj) = 0._wp
204            ELSE
205               fhld (ji,jj) = 0._wp
206            ENDIF
207            !
208            ! -----------------------------------------
209            ! Net heat flux on top of ice-ocean [W.m-2]
210            ! -----------------------------------------
211            !     First  step here      : heat flux at the ocean surface + precip
212            !     Second step below     : heat flux at the ice   surface (after limthd_dif)
213            hfx_in(ji,jj) = hfx_in(ji,jj)                                                                                         & 
214               ! heat flux above the ocean
215               &    +             pfrld(ji,jj)   * ( zqns(ji,jj) + zqsr(ji,jj) )                                                  &
216               ! latent heat of precip (note that precip is included in qns but not in qns_ice)
217               &    +   ( 1._wp - pfrld(ji,jj) ) * sprecip(ji,jj) * ( cpic * ( MIN( tatm_ice(ji,jj), rt0_snow ) - rt0 ) - lfus )  &
218               &    +   ( 1._wp - pfrld(ji,jj) ) * ( tprecip(ji,jj) - sprecip(ji,jj) ) * rcp * ( tatm_ice(ji,jj) - rt0 )
219
220            ! -----------------------------------------------------------------------------
221            ! Net heat flux that is retroceded to the ocean or taken from the ocean [W.m-2]
222            ! -----------------------------------------------------------------------------
223            !     First  step here              :  non solar + precip - qlead - qturb
224            !     Second step in limthd_dh      :  heat remaining if total melt (zq_rema)
225            !     Third  step in limsbc         :  heat from ice-ocean mass exchange (zf_mass) + solar
226            hfx_out(ji,jj) = hfx_out(ji,jj)                                                                                       & 
227               ! Non solar heat flux received by the ocean
228               &    +        pfrld(ji,jj) * qns(ji,jj)                                                                            &
229               ! latent heat of precip (note that precip is included in qns but not in qns_ice)
230               &    +      ( pfrld(ji,jj)**rn_betas - pfrld(ji,jj) ) * sprecip(ji,jj)       &
231               &         * ( cpic * ( MIN( tatm_ice(ji,jj), rt0_snow ) - rt0 ) - lfus )  &
232               &    +      ( 1._wp - pfrld(ji,jj) ) * ( tprecip(ji,jj) - sprecip(ji,jj) ) * rcp * ( tatm_ice(ji,jj) - rt0 )       &
233               ! heat flux taken from the ocean where there is open water ice formation
234               &    -      qlead(ji,jj) * r1_rdtice                                                                               &
235               ! heat flux taken from the ocean during bottom growth/melt (fhld should be 0 while bott growth)
236               &    -      at_i(ji,jj) * fhtur(ji,jj)                                                                             &
237               &    -      at_i(ji,jj) *  fhld(ji,jj)
238
239         END DO
240      END DO
241
242      !------------------------------------------------------------------------------!
243      ! 3) Select icy points and fulfill arrays for the vectorial grid.           
244      !------------------------------------------------------------------------------!
245
246      DO jl = 1, jpl      !loop over ice categories
247
248         IF( kt == nit000 .AND. lwp ) THEN
249            WRITE(numout,*) ' lim_thd : transfer to 1D vectors. Category no : ', jl 
250            WRITE(numout,*) ' ~~~~~~~~'
251         ENDIF
252
253         nbpb = 0
254         DO jj = 1, jpj
255            DO ji = 1, jpi
256               IF ( a_i(ji,jj,jl) > epsi10 ) THEN     
257                  nbpb      = nbpb  + 1
258                  npb(nbpb) = (jj - 1) * jpi + ji
259               ENDIF
260            END DO
261         END DO
262
263         ! debug point to follow
264         jiindex_1d = 0
265         IF( ln_icectl ) THEN
266            DO ji = mi0(iiceprt), mi1(iiceprt)
267               DO jj = mj0(jiceprt), mj1(jiceprt)
268                  jiindex_1d = (jj - 1) * jpi + ji
269                  WRITE(numout,*) ' lim_thd : Category no : ', jl 
270               END DO
271            END DO
272         ENDIF
273
274         !------------------------------------------------------------------------------!
275         ! 4) Thermodynamic computation
276         !------------------------------------------------------------------------------!
277
278         IF( lk_mpp )   CALL mpp_ini_ice( nbpb , numout )
279
280         IF( nbpb > 0 ) THEN  ! If there is no ice, do nothing.
281
282            !-------------------------!
283            ! --- Move to 1D arrays ---
284            !-------------------------!
285            CALL lim_thd_1d2d( nbpb, jl, 1 )
286
287            !--------------------------------------!
288            ! --- Ice/Snow Temperature profile --- !
289            !--------------------------------------!
290            CALL lim_thd_dif( 1, nbpb )
291
292            !---------------------------------!
293            ! --- Ice/Snow thickness ---      !
294            !---------------------------------!
295            CALL lim_thd_dh( 1, nbpb )   
296
297            ! --- Ice enthalpy remapping --- !
298            CALL lim_thd_ent( 1, nbpb, q_i_1d(1:nbpb,:) ) 
299                                           
300            !---------------------------------!
301            ! --- Ice salinity ---            !
302            !---------------------------------!
303            CALL lim_thd_sal( 1, nbpb )   
304
305            !---------------------------------!
306            ! --- temperature update ---      !
307            !---------------------------------!
308            CALL lim_thd_temp( 1, nbpb )
309
310            !------------------------------------!
311            ! --- lateral melting if monocat --- !
312            !------------------------------------!
313            IF ( ( ( nn_monocat == 1 ) .OR. ( nn_monocat == 4 ) ) .AND. ( jpl == 1 ) ) THEN
314               CALL lim_thd_lam( 1, nbpb )
315            END IF
316
317            !-------------------------!
318            ! --- Move to 2D arrays ---
319            !-------------------------!
320            CALL lim_thd_1d2d( nbpb, jl, 2 )
321
322            !
323            IF( lk_mpp )   CALL mpp_comm_free( ncomm_ice ) !RB necessary ??
324         ENDIF
325         !
326      END DO
327
328      !------------------------------------------------------------------------------!
329      ! 5) Global variables, diagnostics
330      !------------------------------------------------------------------------------!
331
332      !------------------------
333      ! Ice heat content             
334      !------------------------
335      ! Enthalpies are global variables we have to readjust the units (heat content in J/m2)
336      DO jl = 1, jpl
337         DO jk = 1, nlay_i
338            e_i(:,:,jk,jl) = e_i(:,:,jk,jl) * a_i(:,:,jl) * ht_i(:,:,jl) * r1_nlay_i
339         END DO
340      END DO
341
342      !------------------------
343      ! Snow heat content             
344      !------------------------
345      ! Enthalpies are global variables we have to readjust the units (heat content in J/m2)
346      DO jl = 1, jpl
347         DO jk = 1, nlay_s
348            e_s(:,:,jk,jl) = e_s(:,:,jk,jl) * a_i(:,:,jl) * ht_s(:,:,jl) * r1_nlay_s
349         END DO
350      END DO
351 
352      !------------------------
353      ! Ice natural aging             
354      !------------------------
355      oa_i(:,:,:) = oa_i(:,:,:) + a_i(:,:,:) * rdt_ice /rday
356
357      !----------------------------------
358      ! Change thickness to volume
359      !----------------------------------
360      CALL lim_var_eqv2glo
361
362      CALL lim_var_zapsmall
363      !--------------------------------------------
364      ! Diagnostic thermodynamic growth rates
365      !--------------------------------------------
366      IF( ln_icectl )   CALL lim_prt( kt, iiceprt, jiceprt, 1, ' - ice thermodyn. - ' )   ! control print
367
368      IF(ln_ctl) THEN            ! Control print
369         CALL prt_ctl_info(' ')
370         CALL prt_ctl_info(' - Cell values : ')
371         CALL prt_ctl_info('   ~~~~~~~~~~~~~ ')
372         CALL prt_ctl(tab2d_1=e12t , clinfo1=' lim_thd  : cell area :')
373         CALL prt_ctl(tab2d_1=at_i , clinfo1=' lim_thd  : at_i      :')
374         CALL prt_ctl(tab2d_1=vt_i , clinfo1=' lim_thd  : vt_i      :')
375         CALL prt_ctl(tab2d_1=vt_s , clinfo1=' lim_thd  : vt_s      :')
376         DO jl = 1, jpl
377            CALL prt_ctl_info(' ')
378            CALL prt_ctl_info(' - Category : ', ivar1=jl)
379            CALL prt_ctl_info('   ~~~~~~~~~~')
380            CALL prt_ctl(tab2d_1=a_i   (:,:,jl)   , clinfo1= ' lim_thd  : a_i      : ')
381            CALL prt_ctl(tab2d_1=ht_i  (:,:,jl)   , clinfo1= ' lim_thd  : ht_i     : ')
382            CALL prt_ctl(tab2d_1=ht_s  (:,:,jl)   , clinfo1= ' lim_thd  : ht_s     : ')
383            CALL prt_ctl(tab2d_1=v_i   (:,:,jl)   , clinfo1= ' lim_thd  : v_i      : ')
384            CALL prt_ctl(tab2d_1=v_s   (:,:,jl)   , clinfo1= ' lim_thd  : v_s      : ')
385            CALL prt_ctl(tab2d_1=e_s   (:,:,1,jl) , clinfo1= ' lim_thd  : e_s      : ')
386            CALL prt_ctl(tab2d_1=t_su  (:,:,jl)   , clinfo1= ' lim_thd  : t_su     : ')
387            CALL prt_ctl(tab2d_1=t_s   (:,:,1,jl) , clinfo1= ' lim_thd  : t_snow   : ')
388            CALL prt_ctl(tab2d_1=sm_i  (:,:,jl)   , clinfo1= ' lim_thd  : sm_i     : ')
389            CALL prt_ctl(tab2d_1=smv_i (:,:,jl)   , clinfo1= ' lim_thd  : smv_i    : ')
390            DO jk = 1, nlay_i
391               CALL prt_ctl_info(' ')
392               CALL prt_ctl_info(' - Layer : ', ivar1=jk)
393               CALL prt_ctl_info('   ~~~~~~~')
394               CALL prt_ctl(tab2d_1=t_i(:,:,jk,jl) , clinfo1= ' lim_thd  : t_i      : ')
395               CALL prt_ctl(tab2d_1=e_i(:,:,jk,jl) , clinfo1= ' lim_thd  : e_i      : ')
396            END DO
397         END DO
398      ENDIF
399      !
400      !
401      CALL wrk_dealloc( jpi, jpj, zqsr, zqns )
402
403      IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limthd', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
404      !------------------------------------------------------------------------------|
405      !  6) Transport of ice between thickness categories.                           |
406      !------------------------------------------------------------------------------|
407      IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limitd_th_rem', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
408
409      ! Given thermodynamic growth rates, transport ice between thickness categories.
410      IF( jpl > 1 )   CALL lim_itd_th_rem( 1, jpl, kt )
411      !
412      CALL lim_var_glo2eqv    ! only for info
413      CALL lim_var_agg(1)
414
415      IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limitd_th_rem', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
416      !------------------------------------------------------------------------------|
417      !  7) Add frazil ice growing in leads.
418      !------------------------------------------------------------------------------|
419      IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limthd_lac', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
420      CALL lim_thd_lac
421      CALL lim_var_glo2eqv    ! only for info
422     
423      ! conservation test
424      IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limthd_lac', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
425
426      IF(ln_ctl) THEN   ! Control print
427         CALL prt_ctl_info(' ')
428         CALL prt_ctl_info(' - Cell values : ')
429         CALL prt_ctl_info('   ~~~~~~~~~~~~~ ')
430         CALL prt_ctl(tab2d_1=e12t , clinfo1=' lim_itd_th  : cell area :')
431         CALL prt_ctl(tab2d_1=at_i , clinfo1=' lim_itd_th  : at_i      :')
432         CALL prt_ctl(tab2d_1=vt_i , clinfo1=' lim_itd_th  : vt_i      :')
433         CALL prt_ctl(tab2d_1=vt_s , clinfo1=' lim_itd_th  : vt_s      :')
434         DO jl = 1, jpl
435            CALL prt_ctl_info(' ')
436            CALL prt_ctl_info(' - Category : ', ivar1=jl)
437            CALL prt_ctl_info('   ~~~~~~~~~~')
438            CALL prt_ctl(tab2d_1=a_i   (:,:,jl)   , clinfo1= ' lim_itd_th  : a_i      : ')
439            CALL prt_ctl(tab2d_1=ht_i  (:,:,jl)   , clinfo1= ' lim_itd_th  : ht_i     : ')
440            CALL prt_ctl(tab2d_1=ht_s  (:,:,jl)   , clinfo1= ' lim_itd_th  : ht_s     : ')
441            CALL prt_ctl(tab2d_1=v_i   (:,:,jl)   , clinfo1= ' lim_itd_th  : v_i      : ')
442            CALL prt_ctl(tab2d_1=v_s   (:,:,jl)   , clinfo1= ' lim_itd_th  : v_s      : ')
443            CALL prt_ctl(tab2d_1=e_s   (:,:,1,jl) , clinfo1= ' lim_itd_th  : e_s      : ')
444            CALL prt_ctl(tab2d_1=t_su  (:,:,jl)   , clinfo1= ' lim_itd_th  : t_su     : ')
445            CALL prt_ctl(tab2d_1=t_s   (:,:,1,jl) , clinfo1= ' lim_itd_th  : t_snow   : ')
446            CALL prt_ctl(tab2d_1=sm_i  (:,:,jl)   , clinfo1= ' lim_itd_th  : sm_i     : ')
447            CALL prt_ctl(tab2d_1=smv_i (:,:,jl)   , clinfo1= ' lim_itd_th  : smv_i    : ')
448            DO jk = 1, nlay_i
449               CALL prt_ctl_info(' ')
450               CALL prt_ctl_info(' - Layer : ', ivar1=jk)
451               CALL prt_ctl_info('   ~~~~~~~')
452               CALL prt_ctl(tab2d_1=t_i(:,:,jk,jl) , clinfo1= ' lim_itd_th  : t_i      : ')
453               CALL prt_ctl(tab2d_1=e_i(:,:,jk,jl) , clinfo1= ' lim_itd_th  : e_i      : ')
454            END DO
455         END DO
456      ENDIF
457      !
458      IF( nn_timing == 1 )  CALL timing_stop('limthd')
459
460   END SUBROUTINE lim_thd 
461
462   SUBROUTINE lim_thd_temp( kideb, kiut )
463      !!-----------------------------------------------------------------------
464      !!                   ***  ROUTINE lim_thd_temp ***
465      !!                 
466      !! ** Purpose :   Computes sea ice temperature (Kelvin) from enthalpy
467      !!
468      !! ** Method  :   Formula (Bitz and Lipscomb, 1999)
469      !!-------------------------------------------------------------------
470      INTEGER, INTENT(in) ::   kideb, kiut   ! bounds for the spatial loop
471      !!
472      INTEGER  ::   ji, jk   ! dummy loop indices
473      REAL(wp) ::   ztmelts, zaaa, zbbb, zccc, zdiscrim  ! local scalar
474      !!-------------------------------------------------------------------
475      ! Recover ice temperature
476      DO jk = 1, nlay_i
477         DO ji = kideb, kiut
478            ztmelts       =  -tmut * s_i_1d(ji,jk) + rt0
479            ! Conversion q(S,T) -> T (second order equation)
480            zaaa          =  cpic
481            zbbb          =  ( rcp - cpic ) * ( ztmelts - rt0 ) + q_i_1d(ji,jk) * r1_rhoic - lfus
482            zccc          =  lfus * ( ztmelts - rt0 )
483            zdiscrim      =  SQRT( MAX( zbbb * zbbb - 4._wp * zaaa * zccc, 0._wp ) )
484            t_i_1d(ji,jk) =  rt0 - ( zbbb + zdiscrim ) / ( 2._wp * zaaa )
485           
486            ! mask temperature
487            rswitch       =  1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_i_1d(ji) ) ) 
488            t_i_1d(ji,jk) =  rswitch * t_i_1d(ji,jk) + ( 1._wp - rswitch ) * rt0
489         END DO
490      END DO
491
492   END SUBROUTINE lim_thd_temp
493
494   SUBROUTINE lim_thd_lam( kideb, kiut )
495      !!-----------------------------------------------------------------------
496      !!                   ***  ROUTINE lim_thd_lam ***
497      !!                 
498      !! ** Purpose :   Lateral melting in case monocategory
499      !!                          ( dA = A/2h dh )
500      !!-----------------------------------------------------------------------
501      INTEGER, INTENT(in) ::   kideb, kiut        ! bounds for the spatial loop
502      INTEGER             ::   ji                 ! dummy loop indices
503      REAL(wp)            ::   zhi_bef            ! ice thickness before thermo
504      REAL(wp)            ::   zdh_mel, zda_mel   ! net melting
505      REAL(wp)            ::   zv                 ! ice volume
506
507      DO ji = kideb, kiut
508         zdh_mel = MIN( 0._wp, dh_i_surf(ji) + dh_i_bott(ji) + dh_snowice(ji) )
509         IF( zdh_mel < 0._wp )  THEN
510            zv         = a_i_1d(ji) * ht_i_1d(ji)
511            ! lateral melting = concentration change
512            zhi_bef     = ht_i_1d(ji) - zdh_mel
513            zda_mel     =  a_i_1d(ji) * zdh_mel / ( 2._wp * MAX( zhi_bef, epsi10 ) )
514            a_i_1d(ji)  = MAX( 0._wp, a_i_1d(ji) + zda_mel ) 
515            ! adjust thickness
516            rswitch     = MAX( 0._wp , SIGN( 1._wp , a_i_1d(ji) - epsi20 ) )
517            ht_i_1d(ji) = rswitch * zv / MAX( a_i_1d(ji), epsi20 )
518            ! retrieve total concentration
519            at_i_1d(ji) = a_i_1d(ji)
520         END IF
521      END DO
522     
523   END SUBROUTINE lim_thd_lam
524
525   SUBROUTINE lim_thd_1d2d( nbpb, jl, kn )
526      !!-----------------------------------------------------------------------
527      !!                   ***  ROUTINE lim_thd_1d2d ***
528      !!                 
529      !! ** Purpose :   move arrays from 1d to 2d and the reverse
530      !!-----------------------------------------------------------------------
531      INTEGER, INTENT(in) ::   kn       ! 1= from 2D to 1D
532                                        ! 2= from 1D to 2D
533      INTEGER, INTENT(in) ::   nbpb     ! size of 1D arrays
534      INTEGER, INTENT(in) ::   jl       ! ice cat
535      INTEGER             ::   jk       ! dummy loop indices
536
537      SELECT CASE( kn )
538
539      CASE( 1 )
540
541         CALL tab_2d_1d( nbpb, at_i_1d     (1:nbpb), at_i            , jpi, jpj, npb(1:nbpb) )
542         CALL tab_2d_1d( nbpb, a_i_1d      (1:nbpb), a_i(:,:,jl)     , jpi, jpj, npb(1:nbpb) )
543         CALL tab_2d_1d( nbpb, ht_i_1d     (1:nbpb), ht_i(:,:,jl)    , jpi, jpj, npb(1:nbpb) )
544         CALL tab_2d_1d( nbpb, ht_s_1d     (1:nbpb), ht_s(:,:,jl)    , jpi, jpj, npb(1:nbpb) )
545         
546         CALL tab_2d_1d( nbpb, t_su_1d     (1:nbpb), t_su(:,:,jl)    , jpi, jpj, npb(1:nbpb) )
547         CALL tab_2d_1d( nbpb, sm_i_1d     (1:nbpb), sm_i(:,:,jl)    , jpi, jpj, npb(1:nbpb) )
548         DO jk = 1, nlay_s
549            CALL tab_2d_1d( nbpb, t_s_1d(1:nbpb,jk), t_s(:,:,jk,jl)  , jpi, jpj, npb(1:nbpb) )
550            CALL tab_2d_1d( nbpb, q_s_1d(1:nbpb,jk), e_s(:,:,jk,jl)  , jpi, jpj, npb(1:nbpb) )
551         END DO
552         DO jk = 1, nlay_i
553            CALL tab_2d_1d( nbpb, t_i_1d(1:nbpb,jk), t_i(:,:,jk,jl)  , jpi, jpj, npb(1:nbpb) )
554            CALL tab_2d_1d( nbpb, q_i_1d(1:nbpb,jk), e_i(:,:,jk,jl)  , jpi, jpj, npb(1:nbpb) )
555            CALL tab_2d_1d( nbpb, s_i_1d(1:nbpb,jk), s_i(:,:,jk,jl)  , jpi, jpj, npb(1:nbpb) )
556         END DO
557         
558         CALL tab_2d_1d( nbpb, tatm_ice_1d(1:nbpb), tatm_ice(:,:)   , jpi, jpj, npb(1:nbpb) )
559         CALL tab_2d_1d( nbpb, qsr_ice_1d (1:nbpb), qsr_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) )
560         CALL tab_2d_1d( nbpb, fr1_i0_1d  (1:nbpb), fr1_i0          , jpi, jpj, npb(1:nbpb) )
561         CALL tab_2d_1d( nbpb, fr2_i0_1d  (1:nbpb), fr2_i0          , jpi, jpj, npb(1:nbpb) )
562         CALL tab_2d_1d( nbpb, qns_ice_1d (1:nbpb), qns_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) )
563         CALL tab_2d_1d( nbpb, ftr_ice_1d (1:nbpb), ftr_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) )
564         IF( .NOT. lk_cpl ) THEN
565            CALL tab_2d_1d( nbpb, qla_ice_1d (1:nbpb), qla_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) )
566            CALL tab_2d_1d( nbpb, dqla_ice_1d(1:nbpb), dqla_ice(:,:,jl), jpi, jpj, npb(1:nbpb) )
567         ENDIF
568         CALL tab_2d_1d( nbpb, dqns_ice_1d(1:nbpb), dqns_ice(:,:,jl), jpi, jpj, npb(1:nbpb) )
569         CALL tab_2d_1d( nbpb, t_bo_1d     (1:nbpb), t_bo            , jpi, jpj, npb(1:nbpb) )
570         CALL tab_2d_1d( nbpb, sprecip_1d (1:nbpb), sprecip         , jpi, jpj, npb(1:nbpb) ) 
571         CALL tab_2d_1d( nbpb, fhtur_1d   (1:nbpb), fhtur           , jpi, jpj, npb(1:nbpb) )
572         CALL tab_2d_1d( nbpb, qlead_1d   (1:nbpb), qlead           , jpi, jpj, npb(1:nbpb) )
573         CALL tab_2d_1d( nbpb, fhld_1d    (1:nbpb), fhld            , jpi, jpj, npb(1:nbpb) )
574         
575         CALL tab_2d_1d( nbpb, wfx_snw_1d (1:nbpb), wfx_snw         , jpi, jpj, npb(1:nbpb) )
576         CALL tab_2d_1d( nbpb, wfx_sub_1d (1:nbpb), wfx_sub         , jpi, jpj, npb(1:nbpb) )
577         
578         CALL tab_2d_1d( nbpb, wfx_bog_1d (1:nbpb), wfx_bog         , jpi, jpj, npb(1:nbpb) )
579         CALL tab_2d_1d( nbpb, wfx_bom_1d (1:nbpb), wfx_bom         , jpi, jpj, npb(1:nbpb) )
580         CALL tab_2d_1d( nbpb, wfx_sum_1d (1:nbpb), wfx_sum         , jpi, jpj, npb(1:nbpb) )
581         CALL tab_2d_1d( nbpb, wfx_sni_1d (1:nbpb), wfx_sni         , jpi, jpj, npb(1:nbpb) )
582         CALL tab_2d_1d( nbpb, wfx_res_1d (1:nbpb), wfx_res         , jpi, jpj, npb(1:nbpb) )
583         CALL tab_2d_1d( nbpb, wfx_spr_1d (1:nbpb), wfx_spr         , jpi, jpj, npb(1:nbpb) )
584         
585         CALL tab_2d_1d( nbpb, sfx_bog_1d (1:nbpb), sfx_bog         , jpi, jpj, npb(1:nbpb) )
586         CALL tab_2d_1d( nbpb, sfx_bom_1d (1:nbpb), sfx_bom         , jpi, jpj, npb(1:nbpb) )
587         CALL tab_2d_1d( nbpb, sfx_sum_1d (1:nbpb), sfx_sum         , jpi, jpj, npb(1:nbpb) )
588         CALL tab_2d_1d( nbpb, sfx_sni_1d (1:nbpb), sfx_sni         , jpi, jpj, npb(1:nbpb) )
589         CALL tab_2d_1d( nbpb, sfx_bri_1d (1:nbpb), sfx_bri         , jpi, jpj, npb(1:nbpb) )
590         CALL tab_2d_1d( nbpb, sfx_res_1d (1:nbpb), sfx_res         , jpi, jpj, npb(1:nbpb) )
591         
592         CALL tab_2d_1d( nbpb, hfx_thd_1d (1:nbpb), hfx_thd         , jpi, jpj, npb(1:nbpb) )
593         CALL tab_2d_1d( nbpb, hfx_spr_1d (1:nbpb), hfx_spr         , jpi, jpj, npb(1:nbpb) )
594         CALL tab_2d_1d( nbpb, hfx_sum_1d (1:nbpb), hfx_sum         , jpi, jpj, npb(1:nbpb) )
595         CALL tab_2d_1d( nbpb, hfx_bom_1d (1:nbpb), hfx_bom         , jpi, jpj, npb(1:nbpb) )
596         CALL tab_2d_1d( nbpb, hfx_bog_1d (1:nbpb), hfx_bog         , jpi, jpj, npb(1:nbpb) )
597         CALL tab_2d_1d( nbpb, hfx_dif_1d (1:nbpb), hfx_dif         , jpi, jpj, npb(1:nbpb) )
598         CALL tab_2d_1d( nbpb, hfx_opw_1d (1:nbpb), hfx_opw         , jpi, jpj, npb(1:nbpb) )
599         CALL tab_2d_1d( nbpb, hfx_snw_1d (1:nbpb), hfx_snw         , jpi, jpj, npb(1:nbpb) )
600         CALL tab_2d_1d( nbpb, hfx_sub_1d (1:nbpb), hfx_sub         , jpi, jpj, npb(1:nbpb) )
601         CALL tab_2d_1d( nbpb, hfx_err_1d (1:nbpb), hfx_err         , jpi, jpj, npb(1:nbpb) )
602         CALL tab_2d_1d( nbpb, hfx_res_1d (1:nbpb), hfx_res         , jpi, jpj, npb(1:nbpb) )
603         CALL tab_2d_1d( nbpb, hfx_err_rem_1d (1:nbpb), hfx_err_rem , jpi, jpj, npb(1:nbpb) )
604
605      CASE( 2 )
606
607         CALL tab_1d_2d( nbpb, at_i          , npb, at_i_1d    (1:nbpb)   , jpi, jpj )
608         CALL tab_1d_2d( nbpb, ht_i(:,:,jl)  , npb, ht_i_1d    (1:nbpb)   , jpi, jpj )
609         CALL tab_1d_2d( nbpb, ht_s(:,:,jl)  , npb, ht_s_1d    (1:nbpb)   , jpi, jpj )
610         CALL tab_1d_2d( nbpb, a_i (:,:,jl)  , npb, a_i_1d     (1:nbpb)   , jpi, jpj )
611         CALL tab_1d_2d( nbpb, t_su(:,:,jl)  , npb, t_su_1d    (1:nbpb)   , jpi, jpj )
612         CALL tab_1d_2d( nbpb, sm_i(:,:,jl)  , npb, sm_i_1d    (1:nbpb)   , jpi, jpj )
613         DO jk = 1, nlay_s
614            CALL tab_1d_2d( nbpb, t_s(:,:,jk,jl), npb, t_s_1d     (1:nbpb,jk), jpi, jpj)
615            CALL tab_1d_2d( nbpb, e_s(:,:,jk,jl), npb, q_s_1d     (1:nbpb,jk), jpi, jpj)
616         END DO
617         DO jk = 1, nlay_i
618            CALL tab_1d_2d( nbpb, t_i(:,:,jk,jl), npb, t_i_1d     (1:nbpb,jk), jpi, jpj)
619            CALL tab_1d_2d( nbpb, e_i(:,:,jk,jl), npb, q_i_1d     (1:nbpb,jk), jpi, jpj)
620            CALL tab_1d_2d( nbpb, s_i(:,:,jk,jl), npb, s_i_1d     (1:nbpb,jk), jpi, jpj)
621         END DO
622         CALL tab_1d_2d( nbpb, qlead         , npb, qlead_1d  (1:nbpb)   , jpi, jpj )
623         
624         CALL tab_1d_2d( nbpb, wfx_snw       , npb, wfx_snw_1d(1:nbpb)   , jpi, jpj )
625         CALL tab_1d_2d( nbpb, wfx_sub       , npb, wfx_sub_1d(1:nbpb)   , jpi, jpj )
626         
627         CALL tab_1d_2d( nbpb, wfx_bog       , npb, wfx_bog_1d(1:nbpb)   , jpi, jpj )
628         CALL tab_1d_2d( nbpb, wfx_bom       , npb, wfx_bom_1d(1:nbpb)   , jpi, jpj )
629         CALL tab_1d_2d( nbpb, wfx_sum       , npb, wfx_sum_1d(1:nbpb)   , jpi, jpj )
630         CALL tab_1d_2d( nbpb, wfx_sni       , npb, wfx_sni_1d(1:nbpb)   , jpi, jpj )
631         CALL tab_1d_2d( nbpb, wfx_res       , npb, wfx_res_1d(1:nbpb)   , jpi, jpj )
632         CALL tab_1d_2d( nbpb, wfx_spr       , npb, wfx_spr_1d(1:nbpb)   , jpi, jpj )
633         
634         CALL tab_1d_2d( nbpb, sfx_bog       , npb, sfx_bog_1d(1:nbpb)   , jpi, jpj )
635         CALL tab_1d_2d( nbpb, sfx_bom       , npb, sfx_bom_1d(1:nbpb)   , jpi, jpj )
636         CALL tab_1d_2d( nbpb, sfx_sum       , npb, sfx_sum_1d(1:nbpb)   , jpi, jpj )
637         CALL tab_1d_2d( nbpb, sfx_sni       , npb, sfx_sni_1d(1:nbpb)   , jpi, jpj )
638         CALL tab_1d_2d( nbpb, sfx_res       , npb, sfx_res_1d(1:nbpb)   , jpi, jpj )
639         CALL tab_1d_2d( nbpb, sfx_bri       , npb, sfx_bri_1d(1:nbpb)   , jpi, jpj )
640         
641         CALL tab_1d_2d( nbpb, hfx_thd       , npb, hfx_thd_1d(1:nbpb)   , jpi, jpj )
642         CALL tab_1d_2d( nbpb, hfx_spr       , npb, hfx_spr_1d(1:nbpb)   , jpi, jpj )
643         CALL tab_1d_2d( nbpb, hfx_sum       , npb, hfx_sum_1d(1:nbpb)   , jpi, jpj )
644         CALL tab_1d_2d( nbpb, hfx_bom       , npb, hfx_bom_1d(1:nbpb)   , jpi, jpj )
645         CALL tab_1d_2d( nbpb, hfx_bog       , npb, hfx_bog_1d(1:nbpb)   , jpi, jpj )
646         CALL tab_1d_2d( nbpb, hfx_dif       , npb, hfx_dif_1d(1:nbpb)   , jpi, jpj )
647         CALL tab_1d_2d( nbpb, hfx_opw       , npb, hfx_opw_1d(1:nbpb)   , jpi, jpj )
648         CALL tab_1d_2d( nbpb, hfx_snw       , npb, hfx_snw_1d(1:nbpb)   , jpi, jpj )
649         CALL tab_1d_2d( nbpb, hfx_sub       , npb, hfx_sub_1d(1:nbpb)   , jpi, jpj )
650         CALL tab_1d_2d( nbpb, hfx_err       , npb, hfx_err_1d(1:nbpb)   , jpi, jpj )
651         CALL tab_1d_2d( nbpb, hfx_res       , npb, hfx_res_1d(1:nbpb)   , jpi, jpj )
652         CALL tab_1d_2d( nbpb, hfx_err_rem   , npb, hfx_err_rem_1d(1:nbpb), jpi, jpj )
653         !
654         CALL tab_1d_2d( nbpb, qns_ice(:,:,jl), npb, qns_ice_1d(1:nbpb) , jpi, jpj)
655         CALL tab_1d_2d( nbpb, ftr_ice(:,:,jl), npb, ftr_ice_1d(1:nbpb) , jpi, jpj )
656                 
657      END SELECT
658
659   END SUBROUTINE lim_thd_1d2d
660
661
662   SUBROUTINE lim_thd_init
663      !!-----------------------------------------------------------------------
664      !!                   ***  ROUTINE lim_thd_init ***
665      !!                 
666      !! ** Purpose :   Physical constants and parameters linked to the ice
667      !!              thermodynamics
668      !!
669      !! ** Method  :   Read the namicethd namelist and check the ice-thermo
670      !!              parameter values called at the first timestep (nit000)
671      !!
672      !! ** input   :   Namelist namicether
673      !!-------------------------------------------------------------------
674      INTEGER  ::   ios                 ! Local integer output status for namelist read
675      NAMELIST/namicethd/ rn_hnewice, ln_frazil, rn_maxfrazb, rn_vfrazb, rn_Cfrazb,                       &
676         &                rn_himin, parsub, rn_betas, rn_kappa_i, nn_conv_dif, rn_terr_dif, nn_ice_thcon, &
677         &                nn_monocat
678      !!-------------------------------------------------------------------
679      !
680      IF(lwp) THEN
681         WRITE(numout,*)
682         WRITE(numout,*) 'lim_thd : Ice Thermodynamics'
683         WRITE(numout,*) '~~~~~~~'
684      ENDIF
685      !
686      REWIND( numnam_ice_ref )              ! Namelist namicethd in reference namelist : Ice thermodynamics
687      READ  ( numnam_ice_ref, namicethd, IOSTAT = ios, ERR = 901)
688901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicethd in reference namelist', lwp )
689
690      REWIND( numnam_ice_cfg )              ! Namelist namicethd in configuration namelist : Ice thermodynamics
691      READ  ( numnam_ice_cfg, namicethd, IOSTAT = ios, ERR = 902 )
692902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicethd in configuration namelist', lwp )
693      IF(lwm) WRITE ( numoni, namicethd )
694      !
695      IF ( ( jpl > 1 ) .AND. ( nn_monocat == 1 ) ) THEN
696         nn_monocat = 0
697         IF(lwp) WRITE(numout, *) '   nn_monocat must be 0 in multi-category case '
698      ENDIF
699
700      IF( lk_cpl .AND. parsub /= 0.0 )   CALL ctl_stop( 'In coupled mode, use parsub = 0. or send dqla' )
701      !
702      IF(lwp) THEN                          ! control print
703         WRITE(numout,*)
704         WRITE(numout,*)'   Namelist of ice parameters for ice thermodynamic computation '
705         WRITE(numout,*)'      ice thick. for lateral accretion                        rn_hnewice   = ', rn_hnewice
706         WRITE(numout,*)'      Frazil ice thickness as a function of wind or not       ln_frazil    = ', ln_frazil
707         WRITE(numout,*)'      Maximum proportion of frazil ice collecting at bottom   rn_maxfrazb  = ', rn_maxfrazb
708         WRITE(numout,*)'      Thresold relative drift speed for collection of frazil  rn_vfrazb    = ', rn_vfrazb
709         WRITE(numout,*)'      Squeezing coefficient for collection of frazil          rn_Cfrazb    = ', rn_Cfrazb
710         WRITE(numout,*)'      minimum ice thickness                                   rn_himin     = ', rn_himin 
711         WRITE(numout,*)'      numerical carac. of the scheme for diffusion in ice '
712         WRITE(numout,*)'      switch for snow sublimation  (=1) or not (=0)           parsub       = ', parsub 
713         WRITE(numout,*)'      coefficient for ice-lead partition of snowfall          rn_betas     = ', rn_betas
714         WRITE(numout,*)'      extinction radiation parameter in sea ice               rn_kappa_i   = ', rn_kappa_i
715         WRITE(numout,*)'      maximal n. of iter. for heat diffusion computation      nn_conv_dif  = ', nn_conv_dif
716         WRITE(numout,*)'      maximal err. on T for heat diffusion computation        rn_terr_dif  = ', rn_terr_dif
717         WRITE(numout,*)'      switch for comp. of thermal conductivity in the ice     nn_ice_thcon = ', nn_ice_thcon
718         WRITE(numout,*)'      check heat conservation in the ice/snow                 con_i        = ', con_i
719         WRITE(numout,*)'      virtual ITD mono-category parameterizations (1) or not  nn_monocat   = ', nn_monocat
720      ENDIF
721      !
722   END SUBROUTINE lim_thd_init
723
724#else
725   !!----------------------------------------------------------------------
726   !!   Default option         Dummy module          NO  LIM3 sea-ice model
727   !!----------------------------------------------------------------------
728#endif
729
730   !!======================================================================
731END MODULE limthd
Note: See TracBrowser for help on using the repository browser.