source: branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90 @ 5207

Last change on this file since 5207 was 5207, checked in by cetlod, 7 years ago

dev_r5204_CNRS_PISCES_dcy:minor bugs corrections

  • Property svn:keywords set to Id
File size: 43.3 KB
Line 
1MODULE limthd_dif
2   !!======================================================================
3   !!                       ***  MODULE limthd_dif ***
4   !!                       heat diffusion in sea ice
5   !!                   computation of surface and inner T 
6   !!======================================================================
7   !! History :  LIM  ! 02-2003 (M. Vancoppenolle) original 1D code
8   !!                 ! 06-2005 (M. Vancoppenolle) 3d version
9   !!                 ! 11-2006 (X Fettweis) Vectorization by Xavier
10   !!                 ! 04-2007 (M. Vancoppenolle) Energy conservation
11   !!            4.0  ! 2011-02 (G. Madec) dynamical allocation
12   !!             -   ! 2012-05 (C. Rousset) add penetration solar flux
13   !!----------------------------------------------------------------------
14#if defined key_lim3
15   !!----------------------------------------------------------------------
16   !!   'key_lim3'                                      LIM3 sea-ice model
17   !!----------------------------------------------------------------------
18   USE par_oce        ! ocean parameters
19   USE phycst         ! physical constants (ocean directory)
20   USE ice            ! LIM-3 variables
21   USE thd_ice        ! LIM-3: thermodynamics
22   USE in_out_manager ! I/O manager
23   USE lib_mpp        ! MPP library
24   USE wrk_nemo       ! work arrays
25   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
26   USE sbc_oce, ONLY : lk_cpl, l_trcdm2dc
27
28   IMPLICIT NONE
29   PRIVATE
30
31   PUBLIC   lim_thd_dif   ! called by lim_thd
32
33   !!----------------------------------------------------------------------
34   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
35   !! $Id$
36   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
37   !!----------------------------------------------------------------------
38CONTAINS
39
40   SUBROUTINE lim_thd_dif( kideb , kiut )
41      !!------------------------------------------------------------------
42      !!                ***  ROUTINE lim_thd_dif  ***
43      !! ** Purpose :
44      !!           This routine determines the time evolution of snow and sea-ice
45      !!           temperature profiles.
46      !! ** Method  :
47      !!           This is done by solving the heat equation diffusion with
48      !!           a Neumann boundary condition at the surface and a Dirichlet one
49      !!           at the bottom. Solar radiation is partially absorbed into the ice.
50      !!           The specific heat and thermal conductivities depend on ice salinity
51      !!           and temperature to take into account brine pocket melting. The
52      !!           numerical
53      !!           scheme is an iterative Crank-Nicolson on a non-uniform multilayer grid
54      !!           in the ice and snow system.
55      !!
56      !!           The successive steps of this routine are
57      !!           1.  Thermal conductivity at the interfaces of the ice layers
58      !!           2.  Internal absorbed radiation
59      !!           3.  Scale factors due to non-uniform grid
60      !!           4.  Kappa factors
61      !!           Then iterative procedure begins
62      !!           5.  specific heat in the ice
63      !!           6.  eta factors
64      !!           7.  surface flux computation
65      !!           8.  tridiagonal system terms
66      !!           9.  solving the tridiagonal system with Gauss elimination
67      !!           Iterative procedure ends according to a criterion on evolution
68      !!           of temperature
69      !!
70      !! ** Arguments :
71      !!           kideb , kiut : Starting and ending points on which the
72      !!                         the computation is applied
73      !!
74      !! ** Inputs / Ouputs : (global commons)
75      !!           surface temperature : t_su_1d
76      !!           ice/snow temperatures   : t_i_1d, t_s_1d
77      !!           ice salinities          : s_i_1d
78      !!           number of layers in the ice/snow: nlay_i, nlay_s
79      !!           profile of the ice/snow layers : z_i, z_s
80      !!           total ice/snow thickness : ht_i_1d, ht_s_1d
81      !!
82      !! ** External :
83      !!
84      !! ** References :
85      !!
86      !! ** History :
87      !!           (02-2003) Martin Vancoppenolle, Louvain-la-Neuve, Belgium
88      !!           (06-2005) Martin Vancoppenolle, 3d version
89      !!           (11-2006) Vectorized by Xavier Fettweis (UCL-ASTR)
90      !!           (04-2007) Energy conservation tested by M. Vancoppenolle
91      !!------------------------------------------------------------------
92      INTEGER , INTENT(in) ::   kideb, kiut   ! Start/End point on which the  the computation is applied
93
94      !! * Local variables
95      INTEGER ::   ji          ! spatial loop index
96      INTEGER ::   ii, ij      ! temporary dummy loop index
97      INTEGER ::   numeq       ! current reference number of equation
98      INTEGER ::   jk       ! vertical dummy loop index
99      INTEGER ::   nconv       ! number of iterations in iterative procedure
100      INTEGER ::   minnumeqmin, maxnumeqmax
101     
102      INTEGER, POINTER, DIMENSION(:) ::   numeqmin   ! reference number of top equation
103      INTEGER, POINTER, DIMENSION(:) ::   numeqmax   ! reference number of bottom equation
104     
105      REAL(wp) ::   zg1s      =  2._wp        ! for the tridiagonal system
106      REAL(wp) ::   zg1       =  2._wp        !
107      REAL(wp) ::   zgamma    =  18009._wp    ! for specific heat
108      REAL(wp) ::   zbeta     =  0.117_wp     ! for thermal conductivity (could be 0.13)
109      REAL(wp) ::   zraext_s  =  10._wp       ! extinction coefficient of radiation in the snow
110      REAL(wp) ::   zkimin    =  0.10_wp      ! minimum ice thermal conductivity
111      REAL(wp) ::   ztsu_err  =  1.e-5_wp     ! range around which t_su is considered as 0°C
112      REAL(wp) ::   ztmelt_i    ! ice melting temperature
113      REAL(wp) ::   zerritmax   ! current maximal error on temperature
114      REAL(wp) ::   zhsu
115     
116      REAL(wp), POINTER, DIMENSION(:)     ::   isnow       ! switch for presence (1) or absence (0) of snow
117      REAL(wp), POINTER, DIMENSION(:)     ::   ztsub       ! old surface temperature (before the iterative procedure )
118      REAL(wp), POINTER, DIMENSION(:)     ::   ztsubit     ! surface temperature at previous iteration
119      REAL(wp), POINTER, DIMENSION(:)     ::   zh_i        ! ice layer thickness
120      REAL(wp), POINTER, DIMENSION(:)     ::   zh_s        ! snow layer thickness
121      REAL(wp), POINTER, DIMENSION(:)     ::   zfsw        ! solar radiation absorbed at the surface
122      REAL(wp), POINTER, DIMENSION(:)     ::   zqns_ice_b  ! solar radiation absorbed at the surface
123      REAL(wp), POINTER, DIMENSION(:)     ::   zf          ! surface flux function
124      REAL(wp), POINTER, DIMENSION(:)     ::   dzf         ! derivative of the surface flux function
125      REAL(wp), POINTER, DIMENSION(:)     ::   zerrit      ! current error on temperature
126      REAL(wp), POINTER, DIMENSION(:)     ::   zdifcase    ! case of the equation resolution (1->4)
127      REAL(wp), POINTER, DIMENSION(:)     ::   zftrice     ! solar radiation transmitted through the ice
128      REAL(wp), POINTER, DIMENSION(:)     ::   zftrice_mean ! daily mean solar radiation transmitted through the ice
129      REAL(wp), POINTER, DIMENSION(:)     ::   zihic
130     
131      REAL(wp), POINTER, DIMENSION(:,:)   ::   ztcond_i    ! Ice thermal conductivity
132      REAL(wp), POINTER, DIMENSION(:,:)   ::   zradtr_i    ! Radiation transmitted through the ice
133      REAL(wp), POINTER, DIMENSION(:,:)   ::   zradab_i    ! Radiation absorbed in the ice
134      REAL(wp), POINTER, DIMENSION(:,:)   ::   zkappa_i    ! Kappa factor in the ice
135      REAL(wp), POINTER, DIMENSION(:,:)   ::   ztib        ! Old temperature in the ice
136      REAL(wp), POINTER, DIMENSION(:,:)   ::   zeta_i      ! Eta factor in the ice
137      REAL(wp), POINTER, DIMENSION(:,:)   ::   ztitemp     ! Temporary temperature in the ice to check the convergence
138      REAL(wp), POINTER, DIMENSION(:,:)   ::   zspeche_i   ! Ice specific heat
139      REAL(wp), POINTER, DIMENSION(:,:)   ::   z_i         ! Vertical cotes of the layers in the ice
140      REAL(wp), POINTER, DIMENSION(:,:)   ::   zradtr_s    ! Radiation transmited through the snow
141      REAL(wp), POINTER, DIMENSION(:,:)   ::   zradab_s    ! Radiation absorbed in the snow
142      REAL(wp), POINTER, DIMENSION(:,:)   ::   zkappa_s    ! Kappa factor in the snow
143      REAL(wp), POINTER, DIMENSION(:,:)   ::   zeta_s      ! Eta factor in the snow
144      REAL(wp), POINTER, DIMENSION(:,:)   ::   ztstemp     ! Temporary temperature in the snow to check the convergence
145      REAL(wp), POINTER, DIMENSION(:,:)   ::   ztsb        ! Temporary temperature in the snow
146      REAL(wp), POINTER, DIMENSION(:,:)   ::   z_s         ! Vertical cotes of the layers in the snow
147      REAL(wp), POINTER, DIMENSION(:,:)   ::   zindterm    ! 'Ind'ependent term
148      REAL(wp), POINTER, DIMENSION(:,:)   ::   zindtbis    ! Temporary 'ind'ependent term
149      REAL(wp), POINTER, DIMENSION(:,:)   ::   zdiagbis    ! Temporary 'dia'gonal term
150      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ztrid       ! Tridiagonal system terms
151     
152      ! diag errors on heat
153      REAL(wp), POINTER, DIMENSION(:)     :: zdq, zq_ini, zhfx_err
154     
155      ! Mono-category
156      REAL(wp)                            :: zepsilon      ! determines thres. above which computation of G(h) is done
157      REAL(wp)                            :: zratio_s      ! dummy factor
158      REAL(wp)                            :: zratio_i      ! dummy factor
159      REAL(wp)                            :: zh_thres      ! thickness thres. for G(h) computation
160      REAL(wp)                            :: zhe           ! dummy factor
161      REAL(wp)                            :: zkimean       ! mean sea ice thermal conductivity
162      REAL(wp)                            :: zfac          ! dummy factor
163      REAL(wp)                            :: zihe          ! dummy factor
164      REAL(wp)                            :: zheshth       ! dummy factor
165     
166      REAL(wp), POINTER, DIMENSION(:)     :: zghe          ! G(he), th. conduct enhancement factor, mono-cat
167     
168      !!------------------------------------------------------------------     
169      !
170      CALL wrk_alloc( jpij, numeqmin, numeqmax )
171      CALL wrk_alloc( jpij, isnow, ztsub, ztsubit, zh_i, zh_s, zfsw )
172      CALL wrk_alloc( jpij, zf, dzf, zqns_ice_b, zerrit, zdifcase, zftrice, zihic, zghe )
173      CALL wrk_alloc( jpij,nlay_i+1, ztcond_i, zradtr_i, zradab_i, zkappa_i, ztib, zeta_i, ztitemp, z_i, zspeche_i, kjstart=0 )
174      CALL wrk_alloc( jpij,nlay_s+1,           zradtr_s, zradab_s, zkappa_s, ztsb, zeta_s, ztstemp, z_s, kjstart=0 )
175      CALL wrk_alloc( jpij,nlay_i+3, zindterm, zindtbis, zdiagbis  )
176      CALL wrk_alloc( jpij,nlay_i+3,3, ztrid )
177      IF( l_trcdm2dc )  CALL wrk_alloc( jpij, zftrice_mean )
178
179      CALL wrk_alloc( jpij, zdq, zq_ini, zhfx_err )
180
181      ! --- diag error on heat diffusion - PART 1 --- !
182      zdq(:) = 0._wp ; zq_ini(:) = 0._wp     
183      DO ji = kideb, kiut
184         zq_ini(ji) = ( SUM( q_i_1d(ji,1:nlay_i) ) * ht_i_1d(ji) * r1_nlay_i +  &
185            &           SUM( q_s_1d(ji,1:nlay_s) ) * ht_s_1d(ji) * r1_nlay_s ) 
186      END DO
187
188      !------------------------------------------------------------------------------!
189      ! 1) Initialization                                                            !
190      !------------------------------------------------------------------------------!
191      DO ji = kideb , kiut
192         isnow(ji)= 1._wp - MAX( 0._wp , SIGN(1._wp, - ht_s_1d(ji) ) )  ! is there snow or not
193         ! layer thickness
194         zh_i(ji) = ht_i_1d(ji) * r1_nlay_i
195         zh_s(ji) = ht_s_1d(ji) * r1_nlay_s
196      END DO
197
198      !--------------------
199      ! Ice / snow layers
200      !--------------------
201
202      z_s(:,0) = 0._wp   ! vert. coord. of the up. lim. of the 1st snow layer
203      z_i(:,0) = 0._wp   ! vert. coord. of the up. lim. of the 1st ice layer
204
205      DO jk = 1, nlay_s            ! vert. coord of the up. lim. of the layer-th snow layer
206         DO ji = kideb , kiut
207            z_s(ji,jk) = z_s(ji,jk-1) + ht_s_1d(ji) * r1_nlay_s
208         END DO
209      END DO
210
211      DO jk = 1, nlay_i            ! vert. coord of the up. lim. of the layer-th ice layer
212         DO ji = kideb , kiut
213            z_i(ji,jk) = z_i(ji,jk-1) + ht_i_1d(ji) * r1_nlay_i
214         END DO
215      END DO
216      !
217      !------------------------------------------------------------------------------|
218      ! 2) Radiation                                                       |
219      !------------------------------------------------------------------------------|
220      !
221      !-------------------
222      ! Computation of i0
223      !-------------------
224      ! i0 describes the fraction of solar radiation which does not contribute
225      ! to the surface energy budget but rather penetrates inside the ice.
226      ! We assume that no radiation is transmitted through the snow
227      ! If there is no no snow
228      ! zfsw    = (1-i0).qsr_ice   is absorbed at the surface
229      ! zftrice = io.qsr_ice       is below the surface
230      ! ftr_ice = io.qsr_ice.exp(-k(h_i)) transmitted below the ice
231      ! fr1_i0_1d = i0 for a thin ice cover, fr1_i0_2d = i0 for a thick ice cover
232      zhsu = 0.1_wp ! threshold for the computation of i0
233      DO ji = kideb , kiut
234         ! switches
235         isnow(ji) = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_s_1d(ji) ) ) 
236         ! hs > 0, isnow = 1
237         zihic(ji) = MAX( 0._wp , 1._wp - ( ht_i_1d(ji) / zhsu ) )     
238
239         i0(ji)    = ( 1._wp - isnow(ji) ) * ( fr1_i0_1d(ji) + zihic(ji) * fr2_i0_1d(ji) )
240      END DO
241
242      !-------------------------------------------------------
243      ! Solar radiation absorbed / transmitted at the surface
244      ! Derivative of the non solar flux
245      !-------------------------------------------------------
246      DO ji = kideb , kiut
247         zfsw   (ji)    =  qsr_ice_1d(ji) * ( 1 - i0(ji) )   ! Shortwave radiation absorbed at surface
248         zftrice(ji)    =  qsr_ice_1d(ji) *       i0(ji)     ! Solar radiation transmitted below the surface layer
249         dzf    (ji)    = dqns_ice_1d(ji)                    ! derivative of incoming nonsolar flux
250         zqns_ice_b(ji) = qns_ice_1d(ji)                     ! store previous qns_ice_1d value
251      END DO
252
253      IF( l_trcdm2dc ) THEN
254         DO ji = kideb , kiut
255            zftrice_mean(ji) =  qsr_ice_mean_1d(ji) * i0(ji) ! Solar radiation transmitted below the surface layer
256         END DO
257      ENDIF
258
259      !---------------------------------------------------------
260      ! Transmission - absorption of solar radiation in the ice
261      !---------------------------------------------------------
262
263      DO ji = kideb, kiut           ! snow initialization
264         zradtr_s(ji,0) = zftrice(ji)     ! radiation penetrating through snow
265      END DO
266
267      DO jk = 1, nlay_s          ! Radiation through snow
268         DO ji = kideb, kiut
269            !                             ! radiation transmitted below the layer-th snow layer
270            zradtr_s(ji,jk) = zradtr_s(ji,0) * EXP( - zraext_s * ( MAX ( 0._wp , z_s(ji,jk) ) ) )
271            !                             ! radiation absorbed by the layer-th snow layer
272            zradab_s(ji,jk) = zradtr_s(ji,jk-1) - zradtr_s(ji,jk)
273         END DO
274      END DO
275
276      DO ji = kideb, kiut           ! ice initialization
277         zradtr_i(ji,0) = zradtr_s(ji,nlay_s) * isnow(ji) + zftrice(ji) * ( 1._wp - isnow(ji) )
278      END DO
279
280      DO jk = 1, nlay_i          ! Radiation through ice
281         DO ji = kideb, kiut
282            !                             ! radiation transmitted below the layer-th ice layer
283            zradtr_i(ji,jk) = zradtr_i(ji,0) * EXP( - rn_kappa_i * ( MAX ( 0._wp , z_i(ji,jk) ) ) )
284            !                             ! radiation absorbed by the layer-th ice layer
285            zradab_i(ji,jk) = zradtr_i(ji,jk-1) - zradtr_i(ji,jk)
286         END DO
287      END DO
288
289      DO ji = kideb, kiut           ! Radiation transmitted below the ice
290         ftr_ice_1d(ji) = zradtr_i(ji,nlay_i) 
291      END DO
292
293
294      IF( l_trcdm2dc ) THEN
295         DO ji = kideb , kiut
296             ftr_ice_mean_1d(ji) =  ftr_ice_mean_1d(ji)               & 
297                 &                + a_i_1d(ji) * zftrice_mean(ji)     &
298                 &                             * EXP( - rn_kappa_i * ( MAX ( 0._wp , ht_i_1d(ji) ) ) ) &
299                 &                             * EXP( - zraext_s   * ( MAX ( 0._wp , ht_s_1d(ji) ) ) )
300         END DO
301      ENDIF
302
303      !
304      !------------------------------------------------------------------------------|
305      !  3) Iterative procedure begins                                               |
306      !------------------------------------------------------------------------------|
307      !
308      DO ji = kideb, kiut        ! Old surface temperature
309         ztsub  (ji) =  t_su_1d(ji)                              ! temperature at the beg of iter pr.
310         ztsubit(ji) =  t_su_1d(ji)                              ! temperature at the previous iter
311         t_su_1d(ji) =  MIN( t_su_1d(ji), rt0 - ztsu_err )       ! necessary
312         zerrit (ji) =  1000._wp                                 ! initial value of error
313      END DO
314
315      DO jk = 1, nlay_s       ! Old snow temperature
316         DO ji = kideb , kiut
317            ztsb(ji,jk) =  t_s_1d(ji,jk)
318         END DO
319      END DO
320
321      DO jk = 1, nlay_i       ! Old ice temperature
322         DO ji = kideb , kiut
323            ztib(ji,jk) =  t_i_1d(ji,jk)
324         END DO
325      END DO
326
327      nconv     =  0           ! number of iterations
328      zerritmax =  1000._wp    ! maximal value of error on all points
329
330      DO WHILE ( zerritmax > rn_terr_dif .AND. nconv < nn_conv_dif )
331         !
332         nconv = nconv + 1
333         !
334         !------------------------------------------------------------------------------|
335         ! 4) Sea ice thermal conductivity                                              |
336         !------------------------------------------------------------------------------|
337         !
338         IF( nn_ice_thcon == 0 ) THEN      ! Untersteiner (1964) formula
339            DO ji = kideb , kiut
340               ztcond_i(ji,0) = rcdic + zbeta * s_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 )
341               ztcond_i(ji,0) = MAX( ztcond_i(ji,0), zkimin )
342            END DO
343            DO jk = 1, nlay_i-1
344               DO ji = kideb , kiut
345                  ztcond_i(ji,jk) = rcdic + zbeta * ( s_i_1d(ji,jk) + s_i_1d(ji,jk+1) ) /  &
346                     MIN(-2.0_wp * epsi10, t_i_1d(ji,jk) + t_i_1d(ji,jk+1) - 2.0_wp * rt0)
347                  ztcond_i(ji,jk) = MAX( ztcond_i(ji,jk), zkimin )
348               END DO
349            END DO
350         ENDIF
351
352         IF( nn_ice_thcon == 1 ) THEN      ! Pringle et al formula included: 2.11 + 0.09 S/T - 0.011.T
353            DO ji = kideb , kiut
354               ztcond_i(ji,0) = rcdic + 0.090_wp * s_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 )   &
355                  &                   - 0.011_wp * ( t_i_1d(ji,1) - rt0 ) 
356               ztcond_i(ji,0) = MAX( ztcond_i(ji,0), zkimin )
357            END DO
358            DO jk = 1, nlay_i-1
359               DO ji = kideb , kiut
360                  ztcond_i(ji,jk) = rcdic +                                                                       & 
361                     &                 0.09_wp * ( s_i_1d(ji,jk) + s_i_1d(ji,jk+1) )                              &
362                     &                 / MIN( -2._wp * epsi10, t_i_1d(ji,jk) + t_i_1d(ji,jk+1) - 2.0_wp * rt0 )   &
363                     &               - 0.0055_wp * ( t_i_1d(ji,jk) + t_i_1d(ji,jk+1) - 2.0 * rt0 ) 
364                  ztcond_i(ji,jk) = MAX( ztcond_i(ji,jk), zkimin )
365               END DO
366            END DO
367            DO ji = kideb , kiut
368               ztcond_i(ji,nlay_i) = rcdic + 0.090_wp * s_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji) - rt0 )   &
369                  &                        - 0.011_wp * ( t_bo_1d(ji) - rt0 ) 
370               ztcond_i(ji,nlay_i) = MAX( ztcond_i(ji,nlay_i), zkimin )
371            END DO
372         ENDIF
373         
374         !
375         !------------------------------------------------------------------------------|
376         !  5) G(he) - enhancement of thermal conductivity in mono-category case        |
377         !------------------------------------------------------------------------------|
378         !
379         ! Computation of effective thermal conductivity G(h)
380         ! Used in mono-category case only to simulate an ITD implicitly
381         ! Fichefet and Morales Maqueda, JGR 1997
382
383         zghe(:) = 1._wp
384
385         SELECT CASE ( nn_monocat )
386
387         CASE (1,3) ! LIM3
388
389            zepsilon = 0.1_wp
390            zh_thres = EXP( 1._wp ) * zepsilon * 0.5_wp
391
392            DO ji = kideb, kiut
393   
394               ! Mean sea ice thermal conductivity
395               zkimean  = SUM( ztcond_i(ji,0:nlay_i) ) / REAL( nlay_i+1, wp )
396
397               ! Effective thickness he (zhe)
398               zfac     = 1._wp / ( rcdsn + zkimean )
399               zratio_s = rcdsn   * zfac
400               zratio_i = zkimean * zfac
401               zhe      = zratio_s * ht_i_1d(ji) + zratio_i * ht_s_1d(ji)
402
403               ! G(he)
404               rswitch  = MAX( 0._wp , SIGN( 1._wp , zhe - zh_thres ) )  ! =0 if zhe < zh_thres, if >
405               zghe(ji) = ( 1._wp - rswitch ) + rswitch * 0.5_wp * ( 1._wp + LOG( 2._wp * zhe / zepsilon ) )
406   
407               ! Impose G(he) < 2.
408               zghe(ji) = MIN( zghe(ji), 2._wp )
409
410            END DO
411
412         END SELECT
413
414         !
415         !------------------------------------------------------------------------------|
416         !  6) kappa factors                                                            |
417         !------------------------------------------------------------------------------|
418         !
419         !--- Snow
420         DO ji = kideb, kiut
421            zfac                  =  1. / MAX( epsi10 , zh_s(ji) )
422            zkappa_s(ji,0)        = zghe(ji) * rcdsn * zfac
423            zkappa_s(ji,nlay_s)   = zghe(ji) * rcdsn * zfac
424         END DO
425
426         DO jk = 1, nlay_s-1
427            DO ji = kideb , kiut
428               zkappa_s(ji,jk)    = zghe(ji) * 2.0 * rcdsn / MAX( epsi10, 2.0 * zh_s(ji) )
429            END DO
430         END DO
431
432         !--- Ice
433         DO jk = 1, nlay_i-1
434            DO ji = kideb , kiut
435               zkappa_i(ji,jk)    = zghe(ji) * 2.0 * ztcond_i(ji,jk) / MAX( epsi10 , 2.0 * zh_i(ji) )
436            END DO
437         END DO
438
439         !--- Snow-ice interface
440         DO ji = kideb , kiut
441            zfac                  = 1./ MAX( epsi10 , zh_i(ji) )
442            zkappa_i(ji,0)        = zghe(ji) * ztcond_i(ji,0) * zfac
443            zkappa_i(ji,nlay_i)   = zghe(ji) * ztcond_i(ji,nlay_i) * zfac
444            zkappa_s(ji,nlay_s)   = zghe(ji) * zghe(ji) * 2.0 * rcdsn * ztcond_i(ji,0) / & 
445           &                        MAX( epsi10, ( zghe(ji) * ztcond_i(ji,0) * zh_s(ji) + zghe(ji) * rcdsn * zh_i(ji) ) )
446            zkappa_i(ji,0)        = zkappa_s(ji,nlay_s) * isnow(ji) + zkappa_i(ji,0) * ( 1._wp - isnow(ji) )
447         END DO
448
449         !
450         !------------------------------------------------------------------------------|
451         ! 7) Sea ice specific heat, eta factors                                        |
452         !------------------------------------------------------------------------------|
453         !
454         DO jk = 1, nlay_i
455            DO ji = kideb , kiut
456               ztitemp(ji,jk)   = t_i_1d(ji,jk)
457               zspeche_i(ji,jk) = cpic + zgamma * s_i_1d(ji,jk) / MAX( ( t_i_1d(ji,jk) - rt0 ) * ( ztib(ji,jk) - rt0 ), epsi10 )
458               zeta_i(ji,jk)    = rdt_ice / MAX( rhoic * zspeche_i(ji,jk) * zh_i(ji), epsi10 )
459            END DO
460         END DO
461
462         DO jk = 1, nlay_s
463            DO ji = kideb , kiut
464               ztstemp(ji,jk) = t_s_1d(ji,jk)
465               zeta_s(ji,jk)  = rdt_ice / MAX( rhosn * cpic * zh_s(ji), epsi10 )
466            END DO
467         END DO
468
469         !
470         !------------------------------------------------------------------------------|
471         ! 8) surface flux computation                                                  |
472         !------------------------------------------------------------------------------|
473         !
474         IF ( ln_it_qnsice ) THEN
475            DO ji = kideb , kiut
476               ! update of the non solar flux according to the update in T_su
477               qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_1d(ji) - ztsubit(ji) )
478            END DO
479         ENDIF
480
481         ! Update incoming flux
482         DO ji = kideb , kiut
483            ! update incoming flux
484            zf(ji)    =          zfsw(ji)              & ! net absorbed solar radiation
485               &         + qns_ice_1d(ji)                ! non solar total flux (LWup, LWdw, SH, LH)
486         END DO
487
488         !
489         !------------------------------------------------------------------------------|
490         ! 9) tridiagonal system terms                                                  |
491         !------------------------------------------------------------------------------|
492         !
493         !!layer denotes the number of the layer in the snow or in the ice
494         !!numeq denotes the reference number of the equation in the tridiagonal
495         !!system, terms of tridiagonal system are indexed as following :
496         !!1 is subdiagonal term, 2 is diagonal and 3 is superdiagonal one
497
498         !!ice interior terms (top equation has the same form as the others)
499
500         DO numeq=1,nlay_i+3
501            DO ji = kideb , kiut
502               ztrid(ji,numeq,1) = 0.
503               ztrid(ji,numeq,2) = 0.
504               ztrid(ji,numeq,3) = 0.
505               zindterm(ji,numeq)= 0.
506               zindtbis(ji,numeq)= 0.
507               zdiagbis(ji,numeq)= 0.
508            ENDDO
509         ENDDO
510
511         DO numeq = nlay_s + 2, nlay_s + nlay_i 
512            DO ji = kideb , kiut
513               jk                 = numeq - nlay_s - 1
514               ztrid(ji,numeq,1)  =  - zeta_i(ji,jk) * zkappa_i(ji,jk-1)
515               ztrid(ji,numeq,2)  =  1.0 + zeta_i(ji,jk) * ( zkappa_i(ji,jk-1) + zkappa_i(ji,jk) )
516               ztrid(ji,numeq,3)  =  - zeta_i(ji,jk) * zkappa_i(ji,jk)
517               zindterm(ji,numeq) =  ztib(ji,jk) + zeta_i(ji,jk) * zradab_i(ji,jk)
518            END DO
519         ENDDO
520
521         numeq =  nlay_s + nlay_i + 1
522         DO ji = kideb , kiut
523            !!ice bottom term
524            ztrid(ji,numeq,1)  =  - zeta_i(ji,nlay_i)*zkappa_i(ji,nlay_i-1)   
525            ztrid(ji,numeq,2)  =  1.0 + zeta_i(ji,nlay_i) * ( zkappa_i(ji,nlay_i) * zg1 + zkappa_i(ji,nlay_i-1) )
526            ztrid(ji,numeq,3)  =  0.0
527            zindterm(ji,numeq) =  ztib(ji,nlay_i) + zeta_i(ji,nlay_i) *  &
528               &                  ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i) * zg1 *  t_bo_1d(ji) ) 
529         ENDDO
530
531
532         DO ji = kideb , kiut
533            IF ( ht_s_1d(ji) > 0.0 ) THEN
534               !
535               !------------------------------------------------------------------------------|
536               !  snow-covered cells                                                          |
537               !------------------------------------------------------------------------------|
538               !
539               !!snow interior terms (bottom equation has the same form as the others)
540               DO numeq = 3, nlay_s + 1
541                  jk                  =  numeq - 1
542                  ztrid(ji,numeq,1)   =  - zeta_s(ji,jk) * zkappa_s(ji,jk-1)
543                  ztrid(ji,numeq,2)   =  1.0 + zeta_s(ji,jk) * ( zkappa_s(ji,jk-1) + zkappa_s(ji,jk) )
544                  ztrid(ji,numeq,3)   =  - zeta_s(ji,jk)*zkappa_s(ji,jk)
545                  zindterm(ji,numeq)  =  ztsb(ji,jk) + zeta_s(ji,jk) * zradab_s(ji,jk)
546               END DO
547
548               !!case of only one layer in the ice (ice equation is altered)
549               IF ( nlay_i.eq.1 ) THEN
550                  ztrid(ji,nlay_s+2,3)    =  0.0
551                  zindterm(ji,nlay_s+2)   =  zindterm(ji,nlay_s+2) + zkappa_i(ji,1) * t_bo_1d(ji) 
552               ENDIF
553
554               IF ( t_su_1d(ji) < rt0 ) THEN
555
556                  !------------------------------------------------------------------------------|
557                  !  case 1 : no surface melting - snow present                                  |
558                  !------------------------------------------------------------------------------|
559                  zdifcase(ji)    =  1.0
560                  numeqmin(ji)    =  1
561                  numeqmax(ji)    =  nlay_i + nlay_s + 1
562
563                  !!surface equation
564                  ztrid(ji,1,1)  = 0.0
565                  ztrid(ji,1,2)  = dzf(ji) - zg1s * zkappa_s(ji,0)
566                  ztrid(ji,1,3)  = zg1s * zkappa_s(ji,0)
567                  zindterm(ji,1) = dzf(ji) * t_su_1d(ji) - zf(ji)
568
569                  !!first layer of snow equation
570                  ztrid(ji,2,1)  =  - zkappa_s(ji,0) * zg1s * zeta_s(ji,1)
571                  ztrid(ji,2,2)  =  1.0 + zeta_s(ji,1) * ( zkappa_s(ji,1) + zkappa_s(ji,0) * zg1s )
572                  ztrid(ji,2,3)  =  - zeta_s(ji,1)* zkappa_s(ji,1)
573                  zindterm(ji,2) =  ztsb(ji,1) + zeta_s(ji,1) * zradab_s(ji,1)
574
575               ELSE 
576                  !
577                  !------------------------------------------------------------------------------|
578                  !  case 2 : surface is melting - snow present                                  |
579                  !------------------------------------------------------------------------------|
580                  !
581                  zdifcase(ji)    =  2.0
582                  numeqmin(ji)    =  2
583                  numeqmax(ji)    =  nlay_i + nlay_s + 1
584
585                  !!first layer of snow equation
586                  ztrid(ji,2,1)  =  0.0
587                  ztrid(ji,2,2)  =  1.0 + zeta_s(ji,1) * ( zkappa_s(ji,1) + zkappa_s(ji,0) * zg1s )
588                  ztrid(ji,2,3)  =  - zeta_s(ji,1)*zkappa_s(ji,1) 
589                  zindterm(ji,2) = ztsb(ji,1) + zeta_s(ji,1) *  &
590                     &             ( zradab_s(ji,1) + zkappa_s(ji,0) * zg1s * t_su_1d(ji) ) 
591               ENDIF
592            ELSE
593               !
594               !------------------------------------------------------------------------------|
595               !  cells without snow                                                          |
596               !------------------------------------------------------------------------------|
597               !
598               IF ( t_su_1d(ji) < rt0 ) THEN
599                  !
600                  !------------------------------------------------------------------------------|
601                  !  case 3 : no surface melting - no snow                                       |
602                  !------------------------------------------------------------------------------|
603                  !
604                  zdifcase(ji)      =  3.0
605                  numeqmin(ji)      =  nlay_s + 1
606                  numeqmax(ji)      =  nlay_i + nlay_s + 1
607
608                  !!surface equation   
609                  ztrid(ji,numeqmin(ji),1)   =  0.0
610                  ztrid(ji,numeqmin(ji),2)   =  dzf(ji) - zkappa_i(ji,0)*zg1   
611                  ztrid(ji,numeqmin(ji),3)   =  zkappa_i(ji,0)*zg1
612                  zindterm(ji,numeqmin(ji))  =  dzf(ji)*t_su_1d(ji) - zf(ji)
613
614                  !!first layer of ice equation
615                  ztrid(ji,numeqmin(ji)+1,1) =  - zkappa_i(ji,0) * zg1 * zeta_i(ji,1)
616                  ztrid(ji,numeqmin(ji)+1,2) =  1.0 + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 )
617                  ztrid(ji,numeqmin(ji)+1,3) =  - zeta_i(ji,1) * zkappa_i(ji,1) 
618                  zindterm(ji,numeqmin(ji)+1)=  ztib(ji,1) + zeta_i(ji,1) * zradab_i(ji,1) 
619
620                  !!case of only one layer in the ice (surface & ice equations are altered)
621
622                  IF ( nlay_i == 1 ) THEN
623                     ztrid(ji,numeqmin(ji),1)    =  0.0
624                     ztrid(ji,numeqmin(ji),2)    =  dzf(ji) - zkappa_i(ji,0) * 2.0
625                     ztrid(ji,numeqmin(ji),3)    =  zkappa_i(ji,0) * 2.0
626                     ztrid(ji,numeqmin(ji)+1,1)  =  -zkappa_i(ji,0) * 2.0 * zeta_i(ji,1)
627                     ztrid(ji,numeqmin(ji)+1,2)  =  1.0 + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2.0 + zkappa_i(ji,1) )
628                     ztrid(ji,numeqmin(ji)+1,3)  =  0.0
629
630                     zindterm(ji,numeqmin(ji)+1) =  ztib(ji,1) + zeta_i(ji,1) *  &
631                        &                          ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) )
632                  ENDIF
633
634               ELSE
635
636                  !
637                  !------------------------------------------------------------------------------|
638                  ! case 4 : surface is melting - no snow                                        |
639                  !------------------------------------------------------------------------------|
640                  !
641                  zdifcase(ji)    =  4.0
642                  numeqmin(ji)    =  nlay_s + 2
643                  numeqmax(ji)    =  nlay_i + nlay_s + 1
644
645                  !!first layer of ice equation
646                  ztrid(ji,numeqmin(ji),1)      =  0.0
647                  ztrid(ji,numeqmin(ji),2)      =  1.0 + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 ) 
648                  ztrid(ji,numeqmin(ji),3)      =  - zeta_i(ji,1) * zkappa_i(ji,1)
649                  zindterm(ji,numeqmin(ji))     =  ztib(ji,1) + zeta_i(ji,1) *  &
650                     &                             ( zradab_i(ji,1) + zkappa_i(ji,0) * zg1 * t_su_1d(ji) ) 
651
652                  !!case of only one layer in the ice (surface & ice equations are altered)
653                  IF ( nlay_i == 1 ) THEN
654                     ztrid(ji,numeqmin(ji),1)  =  0.0
655                     ztrid(ji,numeqmin(ji),2)  =  1.0 + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2.0 + zkappa_i(ji,1) )
656                     ztrid(ji,numeqmin(ji),3)  =  0.0
657                     zindterm(ji,numeqmin(ji)) =  ztib(ji,1) + zeta_i(ji,1) * ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) )  &
658                        &                       + t_su_1d(ji) * zeta_i(ji,1) * zkappa_i(ji,0) * 2.0
659                  ENDIF
660
661               ENDIF
662            ENDIF
663
664         END DO
665
666         !
667         !------------------------------------------------------------------------------|
668         ! 10) tridiagonal system solving                                               |
669         !------------------------------------------------------------------------------|
670         !
671
672         ! Solve the tridiagonal system with Gauss elimination method.
673         ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON,
674         ! McGraw-Hill 1984. 
675
676         maxnumeqmax = 0
677         minnumeqmin = nlay_i+5
678
679         DO ji = kideb , kiut
680            zindtbis(ji,numeqmin(ji)) =  zindterm(ji,numeqmin(ji))
681            zdiagbis(ji,numeqmin(ji)) =  ztrid(ji,numeqmin(ji),2)
682            minnumeqmin               =  MIN(numeqmin(ji),minnumeqmin)
683            maxnumeqmax               =  MAX(numeqmax(ji),maxnumeqmax)
684         END DO
685
686         DO jk = minnumeqmin+1, maxnumeqmax
687            DO ji = kideb , kiut
688               numeq               =  min(max(numeqmin(ji)+1,jk),numeqmax(ji))
689               zdiagbis(ji,numeq)  =  ztrid(ji,numeq,2)  - ztrid(ji,numeq,1) * ztrid(ji,numeq-1,3)  / zdiagbis(ji,numeq-1)
690               zindtbis(ji,numeq)  =  zindterm(ji,numeq) - ztrid(ji,numeq,1) * zindtbis(ji,numeq-1) / zdiagbis(ji,numeq-1)
691            END DO
692         END DO
693
694         DO ji = kideb , kiut
695            ! ice temperatures
696            t_i_1d(ji,nlay_i)    =  zindtbis(ji,numeqmax(ji)) / zdiagbis(ji,numeqmax(ji))
697         END DO
698
699         DO numeq = nlay_i + nlay_s, nlay_s + 2, -1
700            DO ji = kideb , kiut
701               jk    =  numeq - nlay_s - 1
702               t_i_1d(ji,jk)  =  ( zindtbis(ji,numeq) - ztrid(ji,numeq,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,numeq)
703            END DO
704         END DO
705
706         DO ji = kideb , kiut
707            ! snow temperatures     
708            IF (ht_s_1d(ji) > 0._wp) &
709               t_s_1d(ji,nlay_s)     =  ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) )  &
710               &                        / zdiagbis(ji,nlay_s+1) * MAX( 0.0, SIGN( 1.0, ht_s_1d(ji) ) ) 
711
712            ! surface temperature
713            isnow(ji)   = 1._wp - MAX( 0._wp , SIGN( 1._wp , -ht_s_1d(ji) ) )
714            ztsubit(ji) = t_su_1d(ji)
715            IF( t_su_1d(ji) < rt0 ) &
716               t_su_1d(ji) = ( zindtbis(ji,numeqmin(ji)) - ztrid(ji,numeqmin(ji),3) *  &
717               &             ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) ) ) / zdiagbis(ji,numeqmin(ji)) 
718         END DO
719         !
720         !--------------------------------------------------------------------------
721         !  11) Has the scheme converged ?, end of the iterative procedure         |
722         !--------------------------------------------------------------------------
723         !
724         ! check that nowhere it has started to melt
725         ! zerrit(ji) is a measure of error, it has to be under terr_dif
726         DO ji = kideb , kiut
727            t_su_1d(ji) =  MAX(  MIN( t_su_1d(ji) , rt0 ) , 190._wp  )
728            zerrit(ji)  =  ABS( t_su_1d(ji) - ztsubit(ji) )     
729         END DO
730
731         DO jk  =  1, nlay_s
732            DO ji = kideb , kiut
733               t_s_1d(ji,jk) = MAX(  MIN( t_s_1d(ji,jk), rt0 ), 190._wp  )
734               zerrit(ji)    = MAX( zerrit(ji), ABS( t_s_1d(ji,jk) - ztstemp(ji,jk) ) )
735            END DO
736         END DO
737
738         DO jk  =  1, nlay_i
739            DO ji = kideb , kiut
740               ztmelt_i      = -tmut * s_i_1d(ji,jk) + rt0 
741               t_i_1d(ji,jk) =  MAX( MIN( t_i_1d(ji,jk), ztmelt_i ), 190._wp )
742               zerrit(ji)    =  MAX( zerrit(ji), ABS( t_i_1d(ji,jk) - ztitemp(ji,jk) ) )
743            END DO
744         END DO
745
746         ! Compute spatial maximum over all errors
747         ! note that this could be optimized substantially by iterating only the non-converging points
748         zerritmax = 0._wp
749         DO ji = kideb, kiut
750            zerritmax = MAX( zerritmax, zerrit(ji) )   
751         END DO
752         IF( lk_mpp ) CALL mpp_max( zerritmax, kcom=ncomm_ice )
753
754      END DO  ! End of the do while iterative procedure
755
756      IF( ln_icectl .AND. lwp ) THEN
757         WRITE(numout,*) ' zerritmax : ', zerritmax
758         WRITE(numout,*) ' nconv     : ', nconv
759      ENDIF
760
761      !
762      !-------------------------------------------------------------------------!
763      !   12) Fluxes at the interfaces                                          !
764      !-------------------------------------------------------------------------!
765      DO ji = kideb, kiut
766         ! forced mode only : update of latent heat fluxes (sublimation) (always >=0, upward flux)
767         IF( .NOT. lk_cpl) qla_ice_1d (ji) = MAX( 0._wp, qla_ice_1d (ji) + dqla_ice_1d(ji) * ( t_su_1d(ji) - ztsub(ji) ) )
768         !                                ! surface ice conduction flux
769         isnow(ji)       = 1._wp - MAX( 0._wp, SIGN( 1._wp, -ht_s_1d(ji) ) )
770         fc_su(ji)       =  -           isnow(ji)   * zkappa_s(ji,0) * zg1s * (t_s_1d(ji,1) - t_su_1d(ji))   &
771            &               - ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1  * (t_i_1d(ji,1) - t_su_1d(ji))
772         !                                ! bottom ice conduction flux
773         fc_bo_i(ji)     =  - zkappa_i(ji,nlay_i) * ( zg1*(t_bo_1d(ji) - t_i_1d(ji,nlay_i)) )
774      END DO
775
776      ! --- computes sea ice energy of melting compulsory for limthd_dh --- !
777      CALL lim_thd_enmelt( kideb, kiut )
778
779      ! --- diagnose the change in non-solar flux due to surface temperature change --- !
780      IF ( ln_it_qnsice ) hfx_err_dif_1d(:) = hfx_err_dif_1d(:) - ( qns_ice_1d(:)  - zqns_ice_b(:) ) * a_i_1d(:) 
781
782      ! --- diag conservation imbalance on heat diffusion - PART 2 --- !
783      DO ji = kideb, kiut
784         zdq(ji)        = - zq_ini(ji) + ( SUM( q_i_1d(ji,1:nlay_i) ) * ht_i_1d(ji) * r1_nlay_i +  &
785            &                              SUM( q_s_1d(ji,1:nlay_s) ) * ht_s_1d(ji) * r1_nlay_s )
786         IF( t_su_1d(ji) < rt0 ) THEN  ! case T_su < 0degC
787            zhfx_err(ji) = qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq(ji) * r1_rdtice 
788         ELSE                          ! case T_su = 0degC
789            zhfx_err(ji) = fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq(ji) * r1_rdtice 
790         ENDIF
791         hfx_err_1d(ji) = hfx_err_1d(ji) + zhfx_err(ji) * a_i_1d(ji)
792
793         ! total heat that is sent to the ocean (i.e. not used in the heat diffusion equation)
794         hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) + zhfx_err(ji) * a_i_1d(ji)
795      END DO 
796
797      !-----------------------------------------
798      ! Heat flux used to warm/cool ice in W.m-2
799      !-----------------------------------------
800      DO ji = kideb, kiut
801         IF( t_su_1d(ji) < rt0 ) THEN  ! case T_su < 0degC
802            hfx_dif_1d(ji) = hfx_dif_1d(ji)  +   &
803               &            ( qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_1d(ji)
804         ELSE                          ! case T_su = 0degC
805            hfx_dif_1d(ji) = hfx_dif_1d(ji) +    &
806               &             ( fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_1d(ji)
807         ENDIF
808         ! correction on the diagnosed heat flux due to non-convergence of the algorithm used to solve heat equation
809         hfx_dif_1d(ji) = hfx_dif_1d(ji) - zhfx_err(ji) * a_i_1d(ji)
810      END DO   
811      !
812      CALL wrk_dealloc( jpij, numeqmin, numeqmax )
813      CALL wrk_dealloc( jpij, isnow, ztsub, ztsubit, zh_i, zh_s, zfsw )
814      CALL wrk_dealloc( jpij, zf, dzf, zerrit, zdifcase, zftrice, zihic, zghe )
815      CALL wrk_dealloc( jpij,nlay_i+1, ztcond_i, zradtr_i, zradab_i, zkappa_i, ztib, zeta_i, ztitemp, z_i, zspeche_i, kjstart = 0 )
816      CALL wrk_dealloc( jpij,nlay_s+1,           zradtr_s, zradab_s, zkappa_s, ztsb, zeta_s, ztstemp, z_s, kjstart = 0 )
817      CALL wrk_dealloc( jpij,nlay_i+3, zindterm, zindtbis, zdiagbis )
818      CALL wrk_dealloc( jpij,nlay_i+3,3, ztrid )
819      CALL wrk_dealloc( jpij, zdq, zq_ini, zhfx_err )
820
821   END SUBROUTINE lim_thd_dif
822
823   SUBROUTINE lim_thd_enmelt( kideb, kiut )
824      !!-----------------------------------------------------------------------
825      !!                   ***  ROUTINE lim_thd_enmelt ***
826      !!                 
827      !! ** Purpose :   Computes sea ice energy of melting q_i (J.m-3) from temperature
828      !!
829      !! ** Method  :   Formula (Bitz and Lipscomb, 1999)
830      !!-------------------------------------------------------------------
831      INTEGER, INTENT(in) ::   kideb, kiut   ! bounds for the spatial loop
832      !
833      INTEGER  ::   ji, jk   ! dummy loop indices
834      REAL(wp) ::   ztmelts  ! local scalar
835      !!-------------------------------------------------------------------
836      !
837      DO jk = 1, nlay_i             ! Sea ice energy of melting
838         DO ji = kideb, kiut
839            ztmelts      = - tmut  * s_i_1d(ji,jk) + rt0
840            t_i_1d(ji,jk) = MIN( t_i_1d(ji,jk), ztmelts ) ! Force t_i_1d to be lower than melting point
841                                                          !   (sometimes dif scheme produces abnormally high temperatures)   
842            q_i_1d(ji,jk) = rhoic * ( cpic * ( ztmelts - t_i_1d(ji,jk) )                           &
843               &                    + lfus * ( 1.0 - ( ztmelts-rt0 ) / ( t_i_1d(ji,jk) - rt0 ) )   &
844               &                    - rcp  *         ( ztmelts-rt0 )  ) 
845         END DO
846      END DO
847      DO jk = 1, nlay_s             ! Snow energy of melting
848         DO ji = kideb, kiut
849            q_s_1d(ji,jk) = rhosn * ( cpic * ( rt0 - t_s_1d(ji,jk) ) + lfus )
850         END DO
851      END DO
852      !
853   END SUBROUTINE lim_thd_enmelt
854
855#else
856   !!----------------------------------------------------------------------
857   !!                   Dummy Module                 No LIM-3 sea-ice model
858   !!----------------------------------------------------------------------
859CONTAINS
860   SUBROUTINE lim_thd_dif          ! Empty routine
861   END SUBROUTINE lim_thd_dif
862#endif
863   !!======================================================================
864END MODULE limthd_dif
Note: See TracBrowser for help on using the repository browser.