source: branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90 @ 5220

Last change on this file since 5220 was 5220, checked in by smasson, 6 years ago

dev_r5218_CNRS17_coupling: first update

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