source: branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90 @ 5048

Last change on this file since 5048 was 5048, checked in by vancop, 6 years ago

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