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

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

source: branches/UKMO/dev_r5021_nn_etau_revision/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90 @ 5446

Last change on this file since 5446 was 5446, checked in by davestorkey, 9 years ago

Update UKMO/dev_r5021_nn_etau_revision branch to revision 5442 of the trunk.

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