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 trunk/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90 @ 5385

Last change on this file since 5385 was 5385, checked in by cetlod, 9 years ago

merge 2015/dev_r5204_CNRS_PISCES_dcy branch into the trunk, see ticket #1532

  • Property svn:keywords set to Id
File size: 42.5 KB
Line 
1MODULE limthd_dif
2   !!======================================================================
3   !!                       ***  MODULE limthd_dif ***
4   !!                       heat diffusion in sea ice
5   !!                   computation of surface and inner T 
6   !!======================================================================
7   !! History :  LIM  ! 02-2003 (M. Vancoppenolle) original 1D code
8   !!                 ! 06-2005 (M. Vancoppenolle) 3d version
9   !!                 ! 11-2006 (X Fettweis) Vectorization by Xavier
10   !!                 ! 04-2007 (M. Vancoppenolle) Energy conservation
11   !!            4.0  ! 2011-02 (G. Madec) dynamical allocation
12   !!             -   ! 2012-05 (C. Rousset) add penetration solar flux
13   !!----------------------------------------------------------------------
14#if defined key_lim3
15   !!----------------------------------------------------------------------
16   !!   'key_lim3'                                      LIM3 sea-ice model
17   !!----------------------------------------------------------------------
18   USE par_oce        ! ocean parameters
19   USE phycst         ! physical constants (ocean directory)
20   USE ice            ! LIM-3 variables
21   USE thd_ice        ! LIM-3: thermodynamics
22   USE in_out_manager ! I/O manager
23   USE lib_mpp        ! MPP library
24   USE wrk_nemo       ! work arrays
25   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
26   USE sbc_oce, ONLY : lk_cpl
27
28   IMPLICIT NONE
29   PRIVATE
30
31   PUBLIC   lim_thd_dif   ! called by lim_thd
32
33   !!----------------------------------------------------------------------
34   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
35   !! $Id$
36   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
37   !!----------------------------------------------------------------------
38CONTAINS
39
40   SUBROUTINE lim_thd_dif( kideb , kiut )
41      !!------------------------------------------------------------------
42      !!                ***  ROUTINE lim_thd_dif  ***
43      !! ** Purpose :
44      !!           This routine determines the time evolution of snow and sea-ice
45      !!           temperature profiles.
46      !! ** Method  :
47      !!           This is done by solving the heat equation diffusion with
48      !!           a Neumann boundary condition at the surface and a Dirichlet one
49      !!           at the bottom. Solar radiation is partially absorbed into the ice.
50      !!           The specific heat and thermal conductivities depend on ice salinity
51      !!           and temperature to take into account brine pocket melting. The
52      !!           numerical
53      !!           scheme is an iterative Crank-Nicolson on a non-uniform multilayer grid
54      !!           in the ice and snow system.
55      !!
56      !!           The successive steps of this routine are
57      !!           1.  Thermal conductivity at the interfaces of the ice layers
58      !!           2.  Internal absorbed radiation
59      !!           3.  Scale factors due to non-uniform grid
60      !!           4.  Kappa factors
61      !!           Then iterative procedure begins
62      !!           5.  specific heat in the ice
63      !!           6.  eta factors
64      !!           7.  surface flux computation
65      !!           8.  tridiagonal system terms
66      !!           9.  solving the tridiagonal system with Gauss elimination
67      !!           Iterative procedure ends according to a criterion on evolution
68      !!           of temperature
69      !!
70      !! ** Arguments :
71      !!           kideb , kiut : Starting and ending points on which the
72      !!                         the computation is applied
73      !!
74      !! ** Inputs / Ouputs : (global commons)
75      !!           surface temperature : t_su_1d
76      !!           ice/snow temperatures   : t_i_1d, t_s_1d
77      !!           ice salinities          : s_i_1d
78      !!           number of layers in the ice/snow: nlay_i, nlay_s
79      !!           profile of the ice/snow layers : z_i, z_s
80      !!           total ice/snow thickness : ht_i_1d, ht_s_1d
81      !!
82      !! ** External :
83      !!
84      !! ** References :
85      !!
86      !! ** History :
87      !!           (02-2003) Martin Vancoppenolle, Louvain-la-Neuve, Belgium
88      !!           (06-2005) Martin Vancoppenolle, 3d version
89      !!           (11-2006) Vectorized by Xavier Fettweis (UCL-ASTR)
90      !!           (04-2007) Energy conservation tested by M. Vancoppenolle
91      !!------------------------------------------------------------------
92      INTEGER , INTENT(in) ::   kideb, kiut   ! Start/End point on which the  the computation is applied
93
94      !! * Local variables
95      INTEGER ::   ji          ! spatial loop index
96      INTEGER ::   ii, ij      ! temporary dummy loop index
97      INTEGER ::   numeq       ! current reference number of equation
98      INTEGER ::   jk       ! vertical dummy loop index
99      INTEGER ::   nconv       ! number of iterations in iterative procedure
100      INTEGER ::   minnumeqmin, maxnumeqmax
101     
102      INTEGER, POINTER, DIMENSION(:) ::   numeqmin   ! reference number of top equation
103      INTEGER, POINTER, DIMENSION(:) ::   numeqmax   ! reference number of bottom equation
104     
105      REAL(wp) ::   zg1s      =  2._wp        ! for the tridiagonal system
106      REAL(wp) ::   zg1       =  2._wp        !
107      REAL(wp) ::   zgamma    =  18009._wp    ! for specific heat
108      REAL(wp) ::   zbeta     =  0.117_wp     ! for thermal conductivity (could be 0.13)
109      REAL(wp) ::   zraext_s  =  10._wp       ! extinction coefficient of radiation in the snow
110      REAL(wp) ::   zkimin    =  0.10_wp      ! minimum ice thermal conductivity
111      REAL(wp) ::   ztsu_err  =  1.e-5_wp     ! range around which t_su is considered as 0°C
112      REAL(wp) ::   ztmelt_i    ! ice melting temperature
113      REAL(wp) ::   zerritmax   ! current maximal error on temperature
114      REAL(wp) ::   zhsu
115     
116      REAL(wp), POINTER, DIMENSION(:)     ::   isnow       ! switch for presence (1) or absence (0) of snow
117      REAL(wp), POINTER, DIMENSION(:)     ::   ztsub       ! old surface temperature (before the iterative procedure )
118      REAL(wp), POINTER, DIMENSION(:)     ::   ztsubit     ! surface temperature at previous iteration
119      REAL(wp), POINTER, DIMENSION(:)     ::   zh_i        ! ice layer thickness
120      REAL(wp), POINTER, DIMENSION(:)     ::   zh_s        ! snow layer thickness
121      REAL(wp), POINTER, DIMENSION(:)     ::   zfsw        ! solar radiation absorbed at the surface
122      REAL(wp), POINTER, DIMENSION(:)     ::   zqns_ice_b  ! solar radiation absorbed at the surface
123      REAL(wp), POINTER, DIMENSION(:)     ::   zf          ! surface flux function
124      REAL(wp), POINTER, DIMENSION(:)     ::   dzf         ! derivative of the surface flux function
125      REAL(wp), POINTER, DIMENSION(:)     ::   zerrit      ! current error on temperature
126      REAL(wp), POINTER, DIMENSION(:)     ::   zdifcase    ! case of the equation resolution (1->4)
127      REAL(wp), POINTER, DIMENSION(:)     ::   zftrice     ! solar radiation transmitted through the ice
128      REAL(wp), POINTER, DIMENSION(:)     ::   zihic
129     
130      REAL(wp), POINTER, DIMENSION(:,:)   ::   ztcond_i    ! Ice thermal conductivity
131      REAL(wp), POINTER, DIMENSION(:,:)   ::   zradtr_i    ! Radiation transmitted through the ice
132      REAL(wp), POINTER, DIMENSION(:,:)   ::   zradab_i    ! Radiation absorbed in the ice
133      REAL(wp), POINTER, DIMENSION(:,:)   ::   zkappa_i    ! Kappa factor in the ice
134      REAL(wp), POINTER, DIMENSION(:,:)   ::   ztib        ! Old temperature in the ice
135      REAL(wp), POINTER, DIMENSION(:,:)   ::   zeta_i      ! Eta factor in the ice
136      REAL(wp), POINTER, DIMENSION(:,:)   ::   ztitemp     ! Temporary temperature in the ice to check the convergence
137      REAL(wp), POINTER, DIMENSION(:,:)   ::   zspeche_i   ! Ice specific heat
138      REAL(wp), POINTER, DIMENSION(:,:)   ::   z_i         ! Vertical cotes of the layers in the ice
139      REAL(wp), POINTER, DIMENSION(:,:)   ::   zradtr_s    ! Radiation transmited through the snow
140      REAL(wp), POINTER, DIMENSION(:,:)   ::   zradab_s    ! Radiation absorbed in the snow
141      REAL(wp), POINTER, DIMENSION(:,:)   ::   zkappa_s    ! Kappa factor in the snow
142      REAL(wp), POINTER, DIMENSION(:,:)   ::   zeta_s      ! Eta factor in the snow
143      REAL(wp), POINTER, DIMENSION(:,:)   ::   ztstemp     ! Temporary temperature in the snow to check the convergence
144      REAL(wp), POINTER, DIMENSION(:,:)   ::   ztsb        ! Temporary temperature in the snow
145      REAL(wp), POINTER, DIMENSION(:,:)   ::   z_s         ! Vertical cotes of the layers in the snow
146      REAL(wp), POINTER, DIMENSION(:,:)   ::   zindterm    ! 'Ind'ependent term
147      REAL(wp), POINTER, DIMENSION(:,:)   ::   zindtbis    ! Temporary 'ind'ependent term
148      REAL(wp), POINTER, DIMENSION(:,:)   ::   zdiagbis    ! Temporary 'dia'gonal term
149      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ztrid       ! Tridiagonal system terms
150     
151      ! diag errors on heat
152      REAL(wp), POINTER, DIMENSION(:)     :: zdq, zq_ini, zhfx_err
153     
154      ! Mono-category
155      REAL(wp)                            :: zepsilon      ! determines thres. above which computation of G(h) is done
156      REAL(wp)                            :: zratio_s      ! dummy factor
157      REAL(wp)                            :: zratio_i      ! dummy factor
158      REAL(wp)                            :: zh_thres      ! thickness thres. for G(h) computation
159      REAL(wp)                            :: zhe           ! dummy factor
160      REAL(wp)                            :: zkimean       ! mean sea ice thermal conductivity
161      REAL(wp)                            :: zfac          ! dummy factor
162      REAL(wp)                            :: zihe          ! dummy factor
163      REAL(wp)                            :: zheshth       ! dummy factor
164     
165      REAL(wp), POINTER, DIMENSION(:)     :: zghe          ! G(he), th. conduct enhancement factor, mono-cat
166     
167      !!------------------------------------------------------------------     
168      !
169      CALL wrk_alloc( jpij, numeqmin, numeqmax )
170      CALL wrk_alloc( jpij, isnow, ztsub, ztsubit, zh_i, zh_s, zfsw )
171      CALL wrk_alloc( jpij, zf, dzf, zqns_ice_b, zerrit, zdifcase, zftrice, zihic, zghe )
172      CALL wrk_alloc( jpij,nlay_i+1, ztcond_i, zradtr_i, zradab_i, zkappa_i, ztib, zeta_i, ztitemp, z_i, zspeche_i, kjstart=0 )
173      CALL wrk_alloc( jpij,nlay_s+1,           zradtr_s, zradab_s, zkappa_s, ztsb, zeta_s, ztstemp, z_s, kjstart=0 )
174      CALL wrk_alloc( jpij,nlay_i+3, zindterm, zindtbis, zdiagbis  )
175      CALL wrk_alloc( jpij,nlay_i+3,3, ztrid )
176
177      CALL wrk_alloc( jpij, zdq, zq_ini, zhfx_err )
178
179      ! --- diag error on heat diffusion - PART 1 --- !
180      zdq(:) = 0._wp ; zq_ini(:) = 0._wp     
181      DO ji = kideb, kiut
182         zq_ini(ji) = ( SUM( q_i_1d(ji,1:nlay_i) ) * ht_i_1d(ji) * r1_nlay_i +  &
183            &           SUM( q_s_1d(ji,1:nlay_s) ) * ht_s_1d(ji) * r1_nlay_s ) 
184      END DO
185
186      !------------------------------------------------------------------------------!
187      ! 1) Initialization                                                            !
188      !------------------------------------------------------------------------------!
189      DO ji = kideb , kiut
190         isnow(ji)= 1._wp - MAX( 0._wp , SIGN(1._wp, - ht_s_1d(ji) ) )  ! is there snow or not
191         ! layer thickness
192         zh_i(ji) = ht_i_1d(ji) * r1_nlay_i
193         zh_s(ji) = ht_s_1d(ji) * r1_nlay_s
194      END DO
195
196      !--------------------
197      ! Ice / snow layers
198      !--------------------
199
200      z_s(:,0) = 0._wp   ! vert. coord. of the up. lim. of the 1st snow layer
201      z_i(:,0) = 0._wp   ! vert. coord. of the up. lim. of the 1st ice layer
202
203      DO jk = 1, nlay_s            ! vert. coord of the up. lim. of the layer-th snow layer
204         DO ji = kideb , kiut
205            z_s(ji,jk) = z_s(ji,jk-1) + ht_s_1d(ji) * r1_nlay_s
206         END DO
207      END DO
208
209      DO jk = 1, nlay_i            ! vert. coord of the up. lim. of the layer-th ice layer
210         DO ji = kideb , kiut
211            z_i(ji,jk) = z_i(ji,jk-1) + ht_i_1d(ji) * r1_nlay_i
212         END DO
213      END DO
214      !
215      !------------------------------------------------------------------------------|
216      ! 2) Radiation                                                       |
217      !------------------------------------------------------------------------------|
218      !
219      !-------------------
220      ! Computation of i0
221      !-------------------
222      ! i0 describes the fraction of solar radiation which does not contribute
223      ! to the surface energy budget but rather penetrates inside the ice.
224      ! We assume that no radiation is transmitted through the snow
225      ! If there is no no snow
226      ! zfsw    = (1-i0).qsr_ice   is absorbed at the surface
227      ! zftrice = io.qsr_ice       is below the surface
228      ! ftr_ice = io.qsr_ice.exp(-k(h_i)) transmitted below the ice
229      ! fr1_i0_1d = i0 for a thin ice cover, fr1_i0_2d = i0 for a thick ice cover
230      zhsu = 0.1_wp ! threshold for the computation of i0
231      DO ji = kideb , kiut
232         ! switches
233         isnow(ji) = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_s_1d(ji) ) ) 
234         ! hs > 0, isnow = 1
235         zihic(ji) = MAX( 0._wp , 1._wp - ( ht_i_1d(ji) / zhsu ) )     
236
237         i0(ji)    = ( 1._wp - isnow(ji) ) * ( fr1_i0_1d(ji) + zihic(ji) * fr2_i0_1d(ji) )
238      END DO
239
240      !-------------------------------------------------------
241      ! Solar radiation absorbed / transmitted at the surface
242      ! Derivative of the non solar flux
243      !-------------------------------------------------------
244      DO ji = kideb , kiut
245         zfsw   (ji)    =  qsr_ice_1d(ji) * ( 1 - i0(ji) )   ! Shortwave radiation absorbed at surface
246         zftrice(ji)    =  qsr_ice_1d(ji) *       i0(ji)     ! Solar radiation transmitted below the surface layer
247         dzf    (ji)    = dqns_ice_1d(ji)                    ! derivative of incoming nonsolar flux
248         zqns_ice_b(ji) = qns_ice_1d(ji)                     ! store previous qns_ice_1d value
249      END DO
250
251      !---------------------------------------------------------
252      ! Transmission - absorption of solar radiation in the ice
253      !---------------------------------------------------------
254
255      DO ji = kideb, kiut           ! snow initialization
256         zradtr_s(ji,0) = zftrice(ji)     ! radiation penetrating through snow
257      END DO
258
259      DO jk = 1, nlay_s          ! Radiation through snow
260         DO ji = kideb, kiut
261            !                             ! radiation transmitted below the layer-th snow layer
262            zradtr_s(ji,jk) = zradtr_s(ji,0) * EXP( - zraext_s * ( MAX ( 0._wp , z_s(ji,jk) ) ) )
263            !                             ! radiation absorbed by the layer-th snow layer
264            zradab_s(ji,jk) = zradtr_s(ji,jk-1) - zradtr_s(ji,jk)
265         END DO
266      END DO
267
268      DO ji = kideb, kiut           ! ice initialization
269         zradtr_i(ji,0) = zradtr_s(ji,nlay_s) * isnow(ji) + zftrice(ji) * ( 1._wp - isnow(ji) )
270      END DO
271
272      DO jk = 1, nlay_i          ! Radiation through ice
273         DO ji = kideb, kiut
274            !                             ! radiation transmitted below the layer-th ice layer
275            zradtr_i(ji,jk) = zradtr_i(ji,0) * EXP( - rn_kappa_i * ( MAX ( 0._wp , z_i(ji,jk) ) ) )
276            !                             ! radiation absorbed by the layer-th ice layer
277            zradab_i(ji,jk) = zradtr_i(ji,jk-1) - zradtr_i(ji,jk)
278         END DO
279      END DO
280
281      DO ji = kideb, kiut           ! Radiation transmitted below the ice
282         ftr_ice_1d(ji) = zradtr_i(ji,nlay_i) 
283      END DO
284
285      !------------------------------------------------------------------------------|
286      !  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         ! forced mode only : update of latent heat fluxes (sublimation) (always >=0, upward flux)
748         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) ) )
749         !                                ! surface ice conduction flux
750         isnow(ji)       = 1._wp - MAX( 0._wp, SIGN( 1._wp, -ht_s_1d(ji) ) )
751         fc_su(ji)       =  -           isnow(ji)   * zkappa_s(ji,0) * zg1s * (t_s_1d(ji,1) - t_su_1d(ji))   &
752            &               - ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1  * (t_i_1d(ji,1) - t_su_1d(ji))
753         !                                ! bottom ice conduction flux
754         fc_bo_i(ji)     =  - zkappa_i(ji,nlay_i) * ( zg1*(t_bo_1d(ji) - t_i_1d(ji,nlay_i)) )
755      END DO
756
757      ! --- computes sea ice energy of melting compulsory for limthd_dh --- !
758      CALL lim_thd_enmelt( kideb, kiut )
759
760      ! --- diagnose the change in non-solar flux due to surface temperature change --- !
761      IF ( ln_it_qnsice ) hfx_err_dif_1d(:) = hfx_err_dif_1d(:) - ( qns_ice_1d(:)  - zqns_ice_b(:) ) * a_i_1d(:) 
762
763      ! --- diag conservation imbalance on heat diffusion - PART 2 --- !
764      DO ji = kideb, kiut
765         zdq(ji)        = - zq_ini(ji) + ( SUM( q_i_1d(ji,1:nlay_i) ) * ht_i_1d(ji) * r1_nlay_i +  &
766            &                              SUM( q_s_1d(ji,1:nlay_s) ) * ht_s_1d(ji) * r1_nlay_s )
767         IF( t_su_1d(ji) < rt0 ) THEN  ! case T_su < 0degC
768            zhfx_err(ji) = qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq(ji) * r1_rdtice 
769         ELSE                          ! case T_su = 0degC
770            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 
771         ENDIF
772         hfx_err_1d(ji) = hfx_err_1d(ji) + zhfx_err(ji) * a_i_1d(ji)
773
774         ! total heat that is sent to the ocean (i.e. not used in the heat diffusion equation)
775         hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) + zhfx_err(ji) * a_i_1d(ji)
776      END DO 
777
778      !-----------------------------------------
779      ! Heat flux used to warm/cool ice in W.m-2
780      !-----------------------------------------
781      DO ji = kideb, kiut
782         IF( t_su_1d(ji) < rt0 ) THEN  ! case T_su < 0degC
783            hfx_dif_1d(ji) = hfx_dif_1d(ji)  +   &
784               &            ( qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_1d(ji)
785         ELSE                          ! case T_su = 0degC
786            hfx_dif_1d(ji) = hfx_dif_1d(ji) +    &
787               &             ( fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_1d(ji)
788         ENDIF
789         ! correction on the diagnosed heat flux due to non-convergence of the algorithm used to solve heat equation
790         hfx_dif_1d(ji) = hfx_dif_1d(ji) - zhfx_err(ji) * a_i_1d(ji)
791      END DO   
792      !
793      CALL wrk_dealloc( jpij, numeqmin, numeqmax )
794      CALL wrk_dealloc( jpij, isnow, ztsub, ztsubit, zh_i, zh_s, zfsw )
795      CALL wrk_dealloc( jpij, zf, dzf, zerrit, zdifcase, zftrice, zihic, zghe )
796      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 )
797      CALL wrk_dealloc( jpij,nlay_s+1,           zradtr_s, zradab_s, zkappa_s, ztsb, zeta_s, ztstemp, z_s, kjstart = 0 )
798      CALL wrk_dealloc( jpij,nlay_i+3, zindterm, zindtbis, zdiagbis )
799      CALL wrk_dealloc( jpij,nlay_i+3,3, ztrid )
800      CALL wrk_dealloc( jpij, zdq, zq_ini, zhfx_err )
801
802   END SUBROUTINE lim_thd_dif
803
804   SUBROUTINE lim_thd_enmelt( kideb, kiut )
805      !!-----------------------------------------------------------------------
806      !!                   ***  ROUTINE lim_thd_enmelt ***
807      !!                 
808      !! ** Purpose :   Computes sea ice energy of melting q_i (J.m-3) from temperature
809      !!
810      !! ** Method  :   Formula (Bitz and Lipscomb, 1999)
811      !!-------------------------------------------------------------------
812      INTEGER, INTENT(in) ::   kideb, kiut   ! bounds for the spatial loop
813      !
814      INTEGER  ::   ji, jk   ! dummy loop indices
815      REAL(wp) ::   ztmelts  ! local scalar
816      !!-------------------------------------------------------------------
817      !
818      DO jk = 1, nlay_i             ! Sea ice energy of melting
819         DO ji = kideb, kiut
820            ztmelts      = - tmut  * s_i_1d(ji,jk) + rt0
821            t_i_1d(ji,jk) = MIN( t_i_1d(ji,jk), ztmelts ) ! Force t_i_1d to be lower than melting point
822                                                          !   (sometimes dif scheme produces abnormally high temperatures)   
823            q_i_1d(ji,jk) = rhoic * ( cpic * ( ztmelts - t_i_1d(ji,jk) )                           &
824               &                    + lfus * ( 1.0 - ( ztmelts-rt0 ) / ( t_i_1d(ji,jk) - rt0 ) )   &
825               &                    - rcp  *         ( ztmelts-rt0 )  ) 
826         END DO
827      END DO
828      DO jk = 1, nlay_s             ! Snow energy of melting
829         DO ji = kideb, kiut
830            q_s_1d(ji,jk) = rhosn * ( cpic * ( rt0 - t_s_1d(ji,jk) ) + lfus )
831         END DO
832      END DO
833      !
834   END SUBROUTINE lim_thd_enmelt
835
836#else
837   !!----------------------------------------------------------------------
838   !!                   Dummy Module                 No LIM-3 sea-ice model
839   !!----------------------------------------------------------------------
840CONTAINS
841   SUBROUTINE lim_thd_dif          ! Empty routine
842   END SUBROUTINE lim_thd_dif
843#endif
844   !!======================================================================
845END MODULE limthd_dif
Note: See TracBrowser for help on using the repository browser.