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

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

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

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

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