source: NEMO/branches/UKMO/NEMO_4.0.1_remove_0.1m_snow_test/src/ICE/icethd_zdf_bl99.F90 @ 12834

Last change on this file since 12834 was 12834, checked in by dancopsey, 18 months ago

Revert zhs_min back to original value of 0.01

File size: 44.1 KB
Line 
1MODULE icethd_zdf_BL99
2   !!======================================================================
3   !!                       ***  MODULE icethd_zdf_BL99 ***
4   !!   sea-ice: vertical heat diffusion in sea ice (computation of temperatures)
5   !!======================================================================
6   !! History :       !  2003-02  (M. Vancoppenolle) original 1D code
7   !!                 !  2005-06  (M. Vancoppenolle) 3d version
8   !!            4.0  !  2018     (many people)      SI3 [aka Sea Ice cube]
9   !!----------------------------------------------------------------------
10#if defined key_si3
11   !!----------------------------------------------------------------------
12   !!   'key_si3'                                       SI3 sea-ice model
13   !!----------------------------------------------------------------------
14   !!  ice_thd_zdf_BL99 : vertical diffusion computation
15   !!----------------------------------------------------------------------
16   USE dom_oce        ! ocean space and time domain
17   USE phycst         ! physical constants (ocean directory)
18   USE ice            ! sea-ice: variables
19   USE ice1D          ! sea-ice: thermodynamics variables
20   USE icevar         ! sea-ice: operations
21   !
22   USE in_out_manager ! I/O manager
23   USE lib_mpp        ! MPP library
24   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero)
25
26   IMPLICIT NONE
27   PRIVATE
28
29   PUBLIC   ice_thd_zdf_BL99   ! called by icethd_zdf
30
31   !!----------------------------------------------------------------------
32   !! NEMO/ICE 4.0 , NEMO Consortium (2018)
33   !! $Id$
34   !! Software governed by the CeCILL license (see ./LICENSE)
35   !!----------------------------------------------------------------------
36CONTAINS
37
38   SUBROUTINE ice_thd_zdf_BL99( k_cnd )
39      !!-------------------------------------------------------------------
40      !!                ***  ROUTINE ice_thd_zdf_BL99  ***
41      !!
42      !! ** Purpose : computes the time evolution of snow and sea-ice temperature
43      !!              profiles, using the original Bitz and Lipscomb (1999) algorithm
44      !!
45      !! ** Method  : solves the heat equation diffusion with a Neumann boundary
46      !!              condition at the surface and a Dirichlet one at the bottom.
47      !!              Solar radiation is partially absorbed into the ice.
48      !!              The specific heat and thermal conductivities depend on ice
49      !!              salinity and temperature to take into account brine pocket   
50      !!              melting. The numerical scheme is an iterative Crank-Nicolson
51      !!              on a non-uniform multilayer grid in the ice and snow system.
52      !!
53      !!           The successive steps of this routine are
54      !!           1.  initialization of ice-snow layers thicknesses
55      !!           2.  Internal absorbed and transmitted radiation
56      !!           Then iterative procedure begins
57      !!           3.  Thermal conductivity
58      !!           4.  Kappa factors
59      !!           5.  specific heat in the ice
60      !!           6.  eta factors
61      !!           7.  surface flux computation
62      !!           8.  tridiagonal system terms
63      !!           9.  solving the tridiagonal system with Gauss elimination
64      !!           Iterative procedure ends according to a criterion on evolution
65      !!           of temperature
66      !!           10. Fluxes at the interfaces
67      !!
68      !! ** Inputs / Ouputs : (global commons)
69      !!           surface temperature              : t_su_1d
70      !!           ice/snow temperatures            : t_i_1d, t_s_1d
71      !!           ice salinities                   : sz_i_1d
72      !!           number of layers in the ice/snow : nlay_i, nlay_s
73      !!           total ice/snow thickness         : h_i_1d, h_s_1d
74      !!-------------------------------------------------------------------
75      INTEGER, INTENT(in) ::   k_cnd     ! conduction flux (off, on, emulated)
76      !
77      INTEGER ::   ji, jk         ! spatial loop index
78      INTEGER ::   jm             ! current reference number of equation
79      INTEGER ::   jm_mint, jm_maxt
80      INTEGER ::   iconv          ! number of iterations in iterative procedure
81      INTEGER ::   iconv_max = 50 ! max number of iterations in iterative procedure
82      !
83      INTEGER, DIMENSION(jpij) ::   jm_min    ! reference number of top equation
84      INTEGER, DIMENSION(jpij) ::   jm_max    ! reference number of bottom equation
85
86      LOGICAL, DIMENSION(jpij) ::   l_T_converged   ! true when T converges (per grid point)
87!
88      REAL(wp) ::   zg1s      =  2._wp        ! for the tridiagonal system
89      REAL(wp) ::   zg1       =  2._wp        !
90      REAL(wp) ::   zgamma    =  18009._wp    ! for specific heat
91      REAL(wp) ::   zbeta     =  0.117_wp     ! for thermal conductivity (could be 0.13)
92      REAL(wp) ::   zraext_s  =  10._wp       ! extinction coefficient of radiation in the snow
93      REAL(wp) ::   zkimin    =  0.10_wp      ! minimum ice thermal conductivity
94      REAL(wp) ::   ztsu_err  =  1.e-5_wp     ! range around which t_su is considered at 0C
95      REAL(wp) ::   zdti_bnd  =  1.e-4_wp     ! maximal authorized error on temperature
96      REAL(wp) ::   zhs_min   =  0.01_wp      ! minimum snow thickness for conductivity calculation
97      REAL(wp) ::   ztmelts                   ! ice melting temperature
98      REAL(wp) ::   zdti_max                  ! current maximal error on temperature
99      REAL(wp) ::   zcpi                      ! Ice specific heat
100      REAL(wp) ::   zhfx_err, zdq             ! diag errors on heat
101      REAL(wp) ::   zfac                      ! dummy factor
102      !
103      REAL(wp), DIMENSION(jpij) ::   isnow        ! fraction of sea ice that is snow covered (for thermodynamic use only)
104      REAL(wp), DIMENSION(jpij) ::   ztsub        ! surface temperature at previous iteration
105      REAL(wp), DIMENSION(jpij) ::   zh_i, z1_h_i ! ice layer thickness
106      REAL(wp), DIMENSION(jpij) ::   zh_s, z1_h_s ! snow layer thickness
107      REAL(wp), DIMENSION(jpij) ::   zqns_ice_b   ! solar radiation absorbed at the surface
108      REAL(wp), DIMENSION(jpij) ::   zfnet        ! surface flux function
109      REAL(wp), DIMENSION(jpij) ::   zdqns_ice_b  ! derivative of the surface flux function
110      !
111      REAL(wp), DIMENSION(jpij       )     ::   ztsuold     ! Old surface temperature in the ice
112      REAL(wp), DIMENSION(jpij,nlay_i)     ::   ztiold      ! Old temperature in the ice
113      REAL(wp), DIMENSION(jpij,nlay_s)     ::   ztsold      ! Old temperature in the snow
114      REAL(wp), DIMENSION(jpij,nlay_i)     ::   ztib        ! Temporary temperature in the ice to check the convergence
115      REAL(wp), DIMENSION(jpij,nlay_s)     ::   ztsb        ! Temporary temperature in the snow to check the convergence
116      REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   ztcond_i    ! Ice thermal conductivity
117      REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   ztcond_i_cp ! copy
118      REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   zradtr_i    ! Radiation transmitted through the ice
119      REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   zradab_i    ! Radiation absorbed in the ice
120      REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   zkappa_i    ! Kappa factor in the ice
121      REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   zeta_i      ! Eta factor in the ice
122      REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   zradtr_s    ! Radiation transmited through the snow
123      REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   zradab_s    ! Radiation absorbed in the snow
124      REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   zkappa_s    ! Kappa factor in the snow
125      REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   zeta_s      ! Eta factor in the snow
126      REAL(wp), DIMENSION(jpij,nlay_i+3)   ::   zindterm    ! 'Ind'ependent term
127      REAL(wp), DIMENSION(jpij,nlay_i+3)   ::   zindtbis    ! Temporary 'ind'ependent term
128      REAL(wp), DIMENSION(jpij,nlay_i+3)   ::   zdiagbis    ! Temporary 'dia'gonal term
129      REAL(wp), DIMENSION(jpij,nlay_i+3,3) ::   ztrid       ! Tridiagonal system terms
130      REAL(wp), DIMENSION(jpij)            ::   zq_ini      ! diag errors on heat
131      REAL(wp), DIMENSION(jpij)            ::   zghe        ! G(he), th. conduct enhancement factor, mono-cat
132      !
133      ! Mono-category
134      REAL(wp) ::   zepsilon   ! determines thres. above which computation of G(h) is done
135      REAL(wp) ::   zhe        ! dummy factor
136      REAL(wp) ::   zcnd_i     ! mean sea ice thermal conductivity
137      !!------------------------------------------------------------------     
138
139      ! --- diag error on heat diffusion - PART 1 --- !
140      DO ji = 1, npti
141         zq_ini(ji) = ( SUM( e_i_1d(ji,1:nlay_i) ) * h_i_1d(ji) * r1_nlay_i +  &
142            &           SUM( e_s_1d(ji,1:nlay_s) ) * h_s_1d(ji) * r1_nlay_s ) 
143      END DO
144
145      !------------------
146      ! 1) Initialization
147      !------------------
148      DO ji = 1, npti
149
150         ! If the snow thickness drops below zhs_min then reduce the snow fraction instead
151         IF( h_s_1d(ji) < zhs_min ) THEN
152           isnow(ji) = h_s_1d(ji) / zhs_min
153           zh_s(ji) = zhs_min * r1_nlay_s
154         ELSE
155           isnow(ji) = 1.0_wp
156           zh_s(ji) = h_s_1d(ji) * r1_nlay_s
157         END IF
158
159         ! layer thickness
160         zh_i(ji) = h_i_1d(ji) * r1_nlay_i
161
162      END DO
163      !
164      WHERE( zh_i(1:npti) >= epsi10 )   ;   z1_h_i(1:npti) = 1._wp / zh_i(1:npti)
165      ELSEWHERE                         ;   z1_h_i(1:npti) = 0._wp
166      END WHERE
167      !
168      WHERE( zh_s(1:npti) > 0._wp   )   ;   z1_h_s(1:npti) = 1._wp / zh_s(1:npti)
169      ELSEWHERE                         ;   z1_h_s(1:npti) = 0._wp
170      END WHERE
171      !
172      ! Store initial temperatures and non solar heat fluxes
173      IF( k_cnd == np_cnd_OFF .OR. k_cnd == np_cnd_EMU ) THEN
174         !
175         ztsub      (1:npti) = t_su_1d(1:npti)                          ! surface temperature at iteration n-1
176         ztsuold    (1:npti) = t_su_1d(1:npti)                          ! surface temperature initial value
177         t_su_1d    (1:npti) = MIN( t_su_1d(1:npti), rt0 - ztsu_err )   ! required to leave the choice between melting or not
178         zdqns_ice_b(1:npti) = dqns_ice_1d(1:npti)                      ! derivative of incoming nonsolar flux
179         zqns_ice_b (1:npti) = qns_ice_1d(1:npti)                       ! store previous qns_ice_1d value
180         !
181      ENDIF
182      !
183      ztsold (1:npti,:) = t_s_1d(1:npti,:)   ! Old snow temperature
184      ztiold (1:npti,:) = t_i_1d(1:npti,:)   ! Old ice temperature
185
186      !-------------
187      ! 2) Radiation
188      !-------------
189      ! --- Transmission/absorption of solar radiation in the ice --- !
190      zradtr_s(1:npti,0) = qtr_ice_top_1d(1:npti)
191      DO jk = 1, nlay_s
192         DO ji = 1, npti
193            !                             ! radiation transmitted below the layer-th snow layer
194            zradtr_s(ji,jk) = zradtr_s(ji,0) * EXP( - zraext_s * h_s_1d(ji) * r1_nlay_s * REAL(jk) )
195            !                             ! radiation absorbed by the layer-th snow layer
196            zradab_s(ji,jk) = zradtr_s(ji,jk-1) - zradtr_s(ji,jk)
197         END DO
198      END DO
199      !
200      zradtr_i(1:npti,0) = zradtr_s(1:npti,nlay_s) * isnow(1:npti) + qtr_ice_top_1d(1:npti) * ( 1._wp - isnow(1:npti) )
201      DO jk = 1, nlay_i 
202         DO ji = 1, npti
203            !                             ! radiation transmitted below the layer-th ice layer
204            zradtr_i(ji,jk) = zradtr_i(ji,0) * EXP( - rn_kappa_i * zh_i(ji) * REAL(jk) )
205            !                             ! radiation absorbed by the layer-th ice layer
206            zradab_i(ji,jk) = zradtr_i(ji,jk-1) - zradtr_i(ji,jk)
207         END DO
208      END DO
209      !
210      qtr_ice_bot_1d(1:npti) = zradtr_i(1:npti,nlay_i)   ! record radiation transmitted below the ice
211      !
212      iconv    = 0          ! number of iterations
213      !
214      l_T_converged(:) = .FALSE.
215      ! Convergence calculated until all sub-domain grid points have converged
216      ! Calculations keep going for all grid points until sub-domain convergence (vectorisation optimisation)
217      ! but values are not taken into account (results independant of MPI partitioning)
218      !
219      !                                                                            !============================!
220      DO WHILE ( ( .NOT. ALL (l_T_converged(1:npti)) ) .AND. iconv < iconv_max )   ! Iterative procedure begins !
221         !                                                                         !============================!
222         iconv = iconv + 1
223         !
224         ztib(1:npti,:) = t_i_1d(1:npti,:)
225         ztsb(1:npti,:) = t_s_1d(1:npti,:)
226         !
227         !--------------------------------
228         ! 3) Sea ice thermal conductivity
229         !--------------------------------
230         IF( ln_cndi_U64 ) THEN         !-- Untersteiner (1964) formula: k = k0 + beta.S/T
231            !
232            DO ji = 1, npti
233               ztcond_i_cp(ji,0)      = rcnd_i + zbeta * sz_i_1d(ji,1)      / MIN( -epsi10, t_i_1d(ji,1) - rt0 )
234               ztcond_i_cp(ji,nlay_i) = rcnd_i + zbeta * sz_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji)  - rt0 )
235            END DO
236            DO jk = 1, nlay_i-1
237               DO ji = 1, npti
238                  ztcond_i_cp(ji,jk) = rcnd_i + zbeta * 0.5_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) ) /  &
239                     &                         MIN( -epsi10, 0.5_wp * (t_i_1d(ji,jk) + t_i_1d(ji,jk+1)) - rt0 )
240               END DO
241            END DO
242            !
243         ELSEIF( ln_cndi_P07 ) THEN     !-- Pringle et al formula: k = k0 + beta1.S/T - beta2.T
244            !
245            DO ji = 1, npti
246               ztcond_i_cp(ji,0)      = rcnd_i + 0.09_wp  *  sz_i_1d(ji,1)      / MIN( -epsi10, t_i_1d(ji,1) - rt0 )  &
247                  &                           - 0.011_wp * ( t_i_1d(ji,1) - rt0 )
248               ztcond_i_cp(ji,nlay_i) = rcnd_i + 0.09_wp  *  sz_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji)  - rt0 )  &
249                  &                           - 0.011_wp * ( t_bo_1d(ji) - rt0 )
250            END DO
251            DO jk = 1, nlay_i-1
252               DO ji = 1, npti
253                  ztcond_i_cp(ji,jk) = rcnd_i + 0.09_wp  *   0.5_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) ) /        &
254                     &                        MIN( -epsi10, 0.5_wp * ( t_i_1d (ji,jk) + t_i_1d (ji,jk+1) ) - rt0 )  &
255                     &                       - 0.011_wp * ( 0.5_wp * ( t_i_1d (ji,jk) + t_i_1d (ji,jk+1) ) - rt0 )
256               END DO
257            END DO
258            !
259         ENDIF
260
261         ! Variable used after iterations
262         ! Value must be frozen after convergence for MPP independance reason
263         DO ji = 1, npti
264            IF ( .NOT. l_T_converged(ji) ) &
265               ztcond_i(ji,:) = MAX( zkimin, ztcond_i_cp(ji,:) )       
266         END DO
267         !
268         !--- G(he) : enhancement of thermal conductivity in mono-category case
269         ! Computation of effective thermal conductivity G(h)
270         ! Used in mono-category case only to simulate an ITD implicitly
271         ! Fichefet and Morales Maqueda, JGR 1997
272         zghe(1:npti) = 1._wp
273         !
274         IF( ln_virtual_itd ) THEN
275            !
276            zepsilon = 0.1_wp
277            DO ji = 1, npti
278               zcnd_i = SUM( ztcond_i(ji,:) ) / REAL( nlay_i+1, wp )                                ! Mean sea ice thermal conductivity
279               zhe = ( rn_cnd_s * h_i_1d(ji) + zcnd_i * h_s_1d(ji) ) / ( rn_cnd_s + zcnd_i )        ! Effective thickness he (zhe)
280               IF( zhe >=  zepsilon * 0.5_wp * EXP(1._wp) )  &
281                  &   zghe(ji) = MIN( 2._wp, 0.5_wp * ( 1._wp + LOG( 2._wp * zhe / zepsilon ) ) )   ! G(he)
282            END DO
283            !
284         ENDIF
285         !
286         !-----------------
287         ! 4) kappa factors
288         !-----------------
289         !--- Snow
290         ! Variable used after iterations
291         ! Value must be frozen after convergence for MPP independance reason
292         DO jk = 0, nlay_s-1
293            DO ji = 1, npti
294               IF ( .NOT. l_T_converged(ji) ) &
295                  zkappa_s(ji,jk) = zghe(ji) * rn_cnd_s * z1_h_s(ji)
296            END DO
297         END DO
298         DO ji = 1, npti   ! Snow-ice interface
299            IF ( .NOT. l_T_converged(ji) ) THEN
300               zfac = 0.5_wp * ( ztcond_i(ji,0) * zh_s(ji) + rn_cnd_s * zh_i(ji) )
301               IF( zfac > epsi10 ) THEN
302                  zkappa_s(ji,nlay_s) = zghe(ji) * rn_cnd_s * ztcond_i(ji,0) / zfac
303               ELSE
304                  zkappa_s(ji,nlay_s) = 0._wp
305               ENDIF
306            ENDIF
307         END DO
308
309         !--- Ice
310         ! Variable used after iterations
311         ! Value must be frozen after convergence for MPP independance reason
312         DO jk = 0, nlay_i
313            DO ji = 1, npti
314               IF ( .NOT. l_T_converged(ji) ) &
315                  zkappa_i(ji,jk) = zghe(ji) * ztcond_i(ji,jk) * z1_h_i(ji)
316            END DO
317         END DO
318         DO ji = 1, npti   ! Snow-ice interface
319            IF ( .NOT. l_T_converged(ji) ) &
320               zkappa_i(ji,0) = zkappa_s(ji,nlay_s) * isnow(ji) + zkappa_i(ji,0) * ( 1._wp - isnow(ji) )
321         END DO
322         !
323         !--------------------------------------
324         ! 5) Sea ice specific heat, eta factors
325         !--------------------------------------
326         DO jk = 1, nlay_i
327            DO ji = 1, npti
328               zcpi = rcpi + zgamma * sz_i_1d(ji,jk) / MAX( ( t_i_1d(ji,jk) - rt0 ) * ( ztiold(ji,jk) - rt0 ), epsi10 )
329               zeta_i(ji,jk) = rdt_ice * r1_rhoi * z1_h_i(ji) / MAX( epsi10, zcpi ) 
330            END DO
331         END DO
332
333         DO jk = 1, nlay_s
334            DO ji = 1, npti
335               zeta_s(ji,jk) = rdt_ice * r1_rhos * r1_rcpi * z1_h_s(ji)
336            END DO
337         END DO
338         !
339         !----------------------------------------!
340         !                                        !
341         !   Conduction flux is off or emulated   !
342         !                                        !
343         !----------------------------------------!
344         !
345         IF( k_cnd == np_cnd_OFF .OR. k_cnd == np_cnd_EMU ) THEN
346            !
347            ! ==> The original BL99 temperature computation is used
348            !       (with qsr_ice, qns_ice and dqns_ice as inputs)
349            !
350            !----------------------------
351            ! 6) surface flux computation
352            !----------------------------
353            ! update of the non solar flux according to the update in T_su
354            DO ji = 1, npti
355               ! Variable used after iterations
356               ! Value must be frozen after convergence for MPP independance reason
357               IF ( .NOT. l_T_converged(ji) ) &
358                  qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_1d(ji) - ztsub(ji) )
359            END DO
360
361            DO ji = 1, npti
362               zfnet(ji) = qsr_ice_1d(ji) - qtr_ice_top_1d(ji) + qns_ice_1d(ji) ! net heat flux = net - transmitted solar + non solar
363            END DO
364            !
365            !----------------------------
366            ! 7) tridiagonal system terms
367            !----------------------------
368            ! layer denotes the number of the layer in the snow or in the ice
369            ! jm denotes the reference number of the equation in the tridiagonal
370            ! system, terms of tridiagonal system are indexed as following :
371            ! 1 is subdiagonal term, 2 is diagonal and 3 is superdiagonal one
372
373            ! ice interior terms (top equation has the same form as the others)
374            ztrid   (1:npti,:,:) = 0._wp
375            zindterm(1:npti,:)   = 0._wp
376            zindtbis(1:npti,:)   = 0._wp
377            zdiagbis(1:npti,:)   = 0._wp
378
379            DO jm = nlay_s + 2, nlay_s + nlay_i 
380               DO ji = 1, npti
381                  jk = jm - nlay_s - 1
382                  ztrid   (ji,jm,1) =       - zeta_i(ji,jk) *   zkappa_i(ji,jk-1)
383                  ztrid   (ji,jm,2) = 1._wp + zeta_i(ji,jk) * ( zkappa_i(ji,jk-1) + zkappa_i(ji,jk) )
384                  ztrid   (ji,jm,3) =       - zeta_i(ji,jk) *                       zkappa_i(ji,jk)
385                  zindterm(ji,jm)   = ztiold(ji,jk) + zeta_i(ji,jk) * zradab_i(ji,jk)
386               END DO
387            END DO
388
389            jm =  nlay_s + nlay_i + 1
390            DO ji = 1, npti
391               ! ice bottom term
392               ztrid   (ji,jm,1) =       - zeta_i(ji,nlay_i) *   zkappa_i(ji,nlay_i-1)   
393               ztrid   (ji,jm,2) = 1._wp + zeta_i(ji,nlay_i) * ( zkappa_i(ji,nlay_i-1) + zkappa_i(ji,nlay_i) * zg1 )
394               ztrid   (ji,jm,3) = 0._wp
395               zindterm(ji,jm)   = ztiold(ji,nlay_i) + zeta_i(ji,nlay_i) *  &
396                  &              ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i) * zg1 * t_bo_1d(ji) ) 
397            END DO
398
399            DO ji = 1, npti
400               !                               !---------------------!
401               IF( h_s_1d(ji) > 0._wp ) THEN   !  snow-covered cells !
402                  !                            !---------------------!
403                  ! snow interior terms (bottom equation has the same form as the others)
404                  DO jm = 3, nlay_s + 1
405                     jk = jm - 1
406                     ztrid   (ji,jm,1) =       - zeta_s(ji,jk) *   zkappa_s(ji,jk-1)
407                     ztrid   (ji,jm,2) = 1._wp + zeta_s(ji,jk) * ( zkappa_s(ji,jk-1) + zkappa_s(ji,jk) )
408                     ztrid   (ji,jm,3) =       - zeta_s(ji,jk) *                       zkappa_s(ji,jk)
409                     zindterm(ji,jm)   = ztsold(ji,jk) + zeta_s(ji,jk) * zradab_s(ji,jk)
410                  END DO
411                 
412                  ! case of only one layer in the ice (ice equation is altered)
413                  IF( nlay_i == 1 ) THEN
414                     ztrid   (ji,nlay_s+2,3) = 0._wp
415                     zindterm(ji,nlay_s+2)   = zindterm(ji,nlay_s+2) + zeta_i(ji,1) * zkappa_i(ji,1) * t_bo_1d(ji) 
416                  ENDIF
417                 
418                  IF( t_su_1d(ji) < rt0 ) THEN   !--  case 1 : no surface melting
419                     
420                     jm_min(ji) = 1
421                     jm_max(ji) = nlay_i + nlay_s + 1
422                     
423                     ! surface equation
424                     ztrid   (ji,1,1) = 0._wp
425                     ztrid   (ji,1,2) = zdqns_ice_b(ji) - zg1s * zkappa_s(ji,0)
426                     ztrid   (ji,1,3) =                   zg1s * zkappa_s(ji,0)
427                     zindterm(ji,1)   = zdqns_ice_b(ji) * t_su_1d(ji) - zfnet(ji)
428                     
429                     ! first layer of snow equation
430                     ztrid   (ji,2,1) =       - zeta_s(ji,1) *                    zkappa_s(ji,0) * zg1s
431                     ztrid   (ji,2,2) = 1._wp + zeta_s(ji,1) * ( zkappa_s(ji,1) + zkappa_s(ji,0) * zg1s )
432                     ztrid   (ji,2,3) =       - zeta_s(ji,1) *   zkappa_s(ji,1)
433                     zindterm(ji,2)   = ztsold(ji,1) + zeta_s(ji,1) * zradab_s(ji,1)
434                     
435                  ELSE                            !--  case 2 : surface is melting
436                     !
437                     jm_min(ji) = 2
438                     jm_max(ji) = nlay_i + nlay_s + 1
439                     
440                     ! first layer of snow equation
441                     ztrid   (ji,2,1) = 0._wp
442                     ztrid   (ji,2,2) = 1._wp + zeta_s(ji,1) * ( zkappa_s(ji,1) + zkappa_s(ji,0) * zg1s )
443                     ztrid   (ji,2,3) =       - zeta_s(ji,1) *   zkappa_s(ji,1) 
444                     zindterm(ji,2)   = ztsold(ji,1) + zeta_s(ji,1) * ( zradab_s(ji,1) + zkappa_s(ji,0) * zg1s * t_su_1d(ji) ) 
445                  ENDIF
446                  !                            !---------------------!
447               ELSE                            ! cells without snow  !
448                  !                            !---------------------!
449                  !
450                  IF( t_su_1d(ji) < rt0 ) THEN   !--  case 1 : no surface melting
451                     !
452                     jm_min(ji) = nlay_s + 1
453                     jm_max(ji) = nlay_i + nlay_s + 1
454                     
455                     ! surface equation   
456                     ztrid   (ji,jm_min(ji),1) = 0._wp
457                     ztrid   (ji,jm_min(ji),2) = zdqns_ice_b(ji) - zkappa_i(ji,0) * zg1   
458                     ztrid   (ji,jm_min(ji),3) =                   zkappa_i(ji,0) * zg1
459                     zindterm(ji,jm_min(ji))   = zdqns_ice_b(ji) * t_su_1d(ji) - zfnet(ji)
460                     
461                     ! first layer of ice equation
462                     ztrid   (ji,jm_min(ji)+1,1) =       - zeta_i(ji,1) *                    zkappa_i(ji,0) * zg1
463                     ztrid   (ji,jm_min(ji)+1,2) = 1._wp + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 )
464                     ztrid   (ji,jm_min(ji)+1,3) =       - zeta_i(ji,1) *   zkappa_i(ji,1) 
465                     zindterm(ji,jm_min(ji)+1)   = ztiold(ji,1) + zeta_i(ji,1) * zradab_i(ji,1) 
466                     
467                     ! case of only one layer in the ice (surface & ice equations are altered)
468                     IF( nlay_i == 1 ) THEN
469                        ztrid   (ji,jm_min(ji),1)   = 0._wp
470                        ztrid   (ji,jm_min(ji),2)   = zdqns_ice_b(ji)      -   zkappa_i(ji,0) * 2._wp
471                        ztrid   (ji,jm_min(ji),3)   =                          zkappa_i(ji,0) * 2._wp
472                        ztrid   (ji,jm_min(ji)+1,1) =       - zeta_i(ji,1) *   zkappa_i(ji,0) * 2._wp
473                        ztrid   (ji,jm_min(ji)+1,2) = 1._wp + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2._wp + zkappa_i(ji,1) )
474                        ztrid   (ji,jm_min(ji)+1,3) = 0._wp
475                        zindterm(ji,jm_min(ji)+1)   = ztiold(ji,1) + zeta_i(ji,1) * (zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji))
476                     ENDIF
477                     
478                  ELSE                            !--  case 2 : surface is melting
479                     
480                     jm_min(ji) = nlay_s + 2
481                     jm_max(ji) = nlay_i + nlay_s + 1
482                     
483                     ! first layer of ice equation
484                     ztrid   (ji,jm_min(ji),1) = 0._wp
485                     ztrid   (ji,jm_min(ji),2) = 1._wp + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 ) 
486                     ztrid   (ji,jm_min(ji),3) =       - zeta_i(ji,1) *   zkappa_i(ji,1)
487                     zindterm(ji,jm_min(ji))   = ztiold(ji,1) + zeta_i(ji,1) * (zradab_i(ji,1) + zkappa_i(ji,0) * zg1 * t_su_1d(ji)) 
488                     
489                     ! case of only one layer in the ice (surface & ice equations are altered)
490                     IF( nlay_i == 1 ) THEN
491                        ztrid   (ji,jm_min(ji),1) = 0._wp
492                        ztrid   (ji,jm_min(ji),2) = 1._wp + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2._wp + zkappa_i(ji,1) )
493                        ztrid   (ji,jm_min(ji),3) = 0._wp
494                        zindterm(ji,jm_min(ji))   = ztiold(ji,1) + zeta_i(ji,1) * ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) ) &
495                           &                      + t_su_1d(ji) * zeta_i(ji,1) * zkappa_i(ji,0) * 2._wp
496                     ENDIF
497                     
498                  ENDIF
499               ENDIF
500               !
501               zindtbis(ji,jm_min(ji)) = zindterm(ji,jm_min(ji))
502               zdiagbis(ji,jm_min(ji)) = ztrid   (ji,jm_min(ji),2)
503               !
504            END DO
505            !
506            !------------------------------
507            ! 8) tridiagonal system solving
508            !------------------------------
509            ! Solve the tridiagonal system with Gauss elimination method.
510            ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, McGraw-Hill 1984
511            jm_maxt = 0
512            jm_mint = nlay_i+5
513            DO ji = 1, npti
514               jm_mint = MIN(jm_min(ji),jm_mint)
515               jm_maxt = MAX(jm_max(ji),jm_maxt)
516            END DO
517
518            DO jk = jm_mint+1, jm_maxt
519               DO ji = 1, npti
520                  jm = MIN(MAX(jm_min(ji)+1,jk),jm_max(ji))
521                  zdiagbis(ji,jm) = ztrid   (ji,jm,2) - ztrid(ji,jm,1) * ztrid   (ji,jm-1,3) / zdiagbis(ji,jm-1)
522                  zindtbis(ji,jm) = zindterm(ji,jm  ) - ztrid(ji,jm,1) * zindtbis(ji,jm-1  ) / zdiagbis(ji,jm-1)
523               END DO
524            END DO
525
526            ! ice temperatures
527            DO ji = 1, npti
528               ! Variable used after iterations
529               ! Value must be frozen after convergence for MPP independance reason
530               IF ( .NOT. l_T_converged(ji) ) &
531                  t_i_1d(ji,nlay_i) = zindtbis(ji,jm_max(ji)) / zdiagbis(ji,jm_max(ji))
532            END DO
533
534            DO jm = nlay_i + nlay_s, nlay_s + 2, -1
535               DO ji = 1, npti
536                  jk = jm - nlay_s - 1
537                  IF ( .NOT. l_T_converged(ji) ) &
538                     t_i_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,jm)
539               END DO
540            END DO
541
542            DO ji = 1, npti
543               ! Variables used after iterations
544               ! Value must be frozen after convergence for MPP independance reason
545               IF ( .NOT. l_T_converged(ji) ) THEN
546                  ! snow temperatures     
547                  IF( h_s_1d(ji) > 0._wp ) THEN
548                     t_s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) ) / zdiagbis(ji,nlay_s+1)
549                  ENDIF
550                  ! surface temperature
551                  ztsub(ji) = t_su_1d(ji)
552                  IF( t_su_1d(ji) < rt0 ) THEN
553                     t_su_1d(ji) = ( zindtbis(ji,jm_min(ji)) - ztrid(ji,jm_min(ji),3) *  &
554                        &          ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) ) ) / zdiagbis(ji,jm_min(ji))
555                  ENDIF
556               ENDIF
557            END DO
558            !
559            !--------------------------------------------------------------
560            ! 9) Has the scheme converged?, end of the iterative procedure
561            !--------------------------------------------------------------
562            ! check that nowhere it has started to melt
563            ! zdti_max is a measure of error, it has to be under zdti_bnd
564
565            DO ji = 1, npti
566
567               zdti_max = 0._wp
568
569               IF ( .NOT. l_T_converged(ji) ) THEN
570                  t_su_1d(ji) = MAX( MIN( t_su_1d(ji) , rt0 ) , rt0 - 100._wp )
571                  zdti_max    = MAX( zdti_max, ABS( t_su_1d(ji) - ztsub(ji) ) )
572
573                  t_s_1d(ji,1:nlay_s) = MAX( MIN( t_s_1d(ji,1:nlay_s), rt0 ), rt0 - 100._wp )
574                  zdti_max = MAX ( zdti_max , MAXVAL( ABS( t_s_1d(ji,1:nlay_s) - ztsb(ji,1:nlay_s) ) ) )
575
576                  DO jk = 1, nlay_i
577                     ztmelts       = -rTmlt * sz_i_1d(ji,jk) + rt0
578                     t_i_1d(ji,jk) =  MAX( MIN( t_i_1d(ji,jk), ztmelts ), rt0 - 100._wp )
579                     zdti_max      =  MAX( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) )
580                  END DO
581
582                  IF ( zdti_max < zdti_bnd ) l_T_converged(ji) = .TRUE.
583
584               ENDIF
585
586            END DO
587
588         !----------------------------------------!
589         !                                        !
590         !      Conduction flux is on             !
591         !                                        !
592         !----------------------------------------!
593         !
594         ELSEIF( k_cnd == np_cnd_ON ) THEN
595            !
596            ! ==> we use a modified BL99 solver with conduction flux (qcn_ice) as forcing term
597            !
598            !----------------------------
599            ! 7) tridiagonal system terms
600            !----------------------------
601            ! layer denotes the number of the layer in the snow or in the ice
602            ! jm denotes the reference number of the equation in the tridiagonal
603            ! system, terms of tridiagonal system are indexed as following :
604            ! 1 is subdiagonal term, 2 is diagonal and 3 is superdiagonal one
605
606            ! ice interior terms (top equation has the same form as the others)
607            ztrid   (1:npti,:,:) = 0._wp
608            zindterm(1:npti,:)   = 0._wp
609            zindtbis(1:npti,:)   = 0._wp
610            zdiagbis(1:npti,:)   = 0._wp
611
612            DO jm = nlay_s + 2, nlay_s + nlay_i 
613               DO ji = 1, npti
614                  jk = jm - nlay_s - 1
615                  ztrid   (ji,jm,1) =       - zeta_i(ji,jk) *   zkappa_i(ji,jk-1)
616                  ztrid   (ji,jm,2) = 1._wp + zeta_i(ji,jk) * ( zkappa_i(ji,jk-1) + zkappa_i(ji,jk) )
617                  ztrid   (ji,jm,3) =       - zeta_i(ji,jk) *                       zkappa_i(ji,jk)
618                  zindterm(ji,jm)   = ztiold(ji,jk) + zeta_i(ji,jk) * zradab_i(ji,jk)
619               END DO
620            ENDDO
621
622            jm =  nlay_s + nlay_i + 1
623            DO ji = 1, npti
624               ! ice bottom term
625               ztrid   (ji,jm,1) =       - zeta_i(ji,nlay_i) *   zkappa_i(ji,nlay_i-1)   
626               ztrid   (ji,jm,2) = 1._wp + zeta_i(ji,nlay_i) * ( zkappa_i(ji,nlay_i-1) + zkappa_i(ji,nlay_i) * zg1 )
627               ztrid   (ji,jm,3) = 0._wp
628               zindterm(ji,jm)   = ztiold(ji,nlay_i) + zeta_i(ji,nlay_i) *  &
629                  &              ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i) * zg1 * t_bo_1d(ji) ) 
630            ENDDO
631
632            DO ji = 1, npti
633               !                               !---------------------!
634               IF( h_s_1d(ji) > 0._wp ) THEN   !  snow-covered cells !
635                  !                            !---------------------!
636                  ! snow interior terms (bottom equation has the same form as the others)
637                  DO jm = 3, nlay_s + 1
638                     jk = jm - 1
639                     ztrid   (ji,jm,1) =       - zeta_s(ji,jk) *   zkappa_s(ji,jk-1)
640                     ztrid   (ji,jm,2) = 1._wp + zeta_s(ji,jk) * ( zkappa_s(ji,jk-1) + zkappa_s(ji,jk) )
641                     ztrid   (ji,jm,3) =       - zeta_s(ji,jk) *                       zkappa_s(ji,jk)
642                     zindterm(ji,jm)   = ztsold(ji,jk) + zeta_s(ji,jk) * zradab_s(ji,jk)
643                  END DO
644                 
645                  ! case of only one layer in the ice (ice equation is altered)
646                  IF ( nlay_i == 1 ) THEN
647                     ztrid   (ji,nlay_s+2,3) = 0._wp
648                     zindterm(ji,nlay_s+2)   = zindterm(ji,nlay_s+2) + zeta_i(ji,1) * zkappa_i(ji,1) * t_bo_1d(ji) 
649                  ENDIF
650                 
651                  jm_min(ji) = 2
652                  jm_max(ji) = nlay_i + nlay_s + 1
653                 
654                  ! first layer of snow equation
655                  ztrid   (ji,2,1) = 0._wp
656                  ztrid   (ji,2,2) = 1._wp + zeta_s(ji,1) * zkappa_s(ji,1)
657                  ztrid   (ji,2,3) =       - zeta_s(ji,1) * zkappa_s(ji,1) 
658                  zindterm(ji,2)   = ztsold(ji,1) + zeta_s(ji,1) * ( zradab_s(ji,1) + qcn_ice_1d(ji) ) 
659                 
660                  !                            !---------------------!
661               ELSE                            ! cells without snow  !
662                  !                            !---------------------!
663                  jm_min(ji) = nlay_s + 2
664                  jm_max(ji) = nlay_i + nlay_s + 1
665                 
666                  ! first layer of ice equation
667                  ztrid   (ji,jm_min(ji),1) = 0._wp
668                  ztrid   (ji,jm_min(ji),2) = 1._wp + zeta_i(ji,1) * zkappa_i(ji,1) 
669                  ztrid   (ji,jm_min(ji),3) =       - zeta_i(ji,1) * zkappa_i(ji,1)
670                  zindterm(ji,jm_min(ji))   = ztiold(ji,1) + zeta_i(ji,1) * ( zradab_i(ji,1) + qcn_ice_1d(ji) )
671                 
672                  ! case of only one layer in the ice (surface & ice equations are altered)
673                  IF( nlay_i == 1 ) THEN
674                     ztrid   (ji,jm_min(ji),1) = 0._wp
675                     ztrid   (ji,jm_min(ji),2) = 1._wp + zeta_i(ji,1) * zkappa_i(ji,1)
676                     ztrid   (ji,jm_min(ji),3) = 0._wp
677                     zindterm(ji,jm_min(ji))   = ztiold(ji,1) + zeta_i(ji,1) *  &
678                        &                                     ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) + qcn_ice_1d(ji) )
679                  ENDIF
680                 
681               ENDIF
682               !
683               zindtbis(ji,jm_min(ji)) = zindterm(ji,jm_min(ji))
684               zdiagbis(ji,jm_min(ji)) = ztrid   (ji,jm_min(ji),2)
685               !
686            END DO
687            !
688            !------------------------------
689            ! 8) tridiagonal system solving
690            !------------------------------
691            ! Solve the tridiagonal system with Gauss elimination method.
692            ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, McGraw-Hill 1984
693            jm_maxt = 0
694            jm_mint = nlay_i+5
695            DO ji = 1, npti
696               jm_mint = MIN(jm_min(ji),jm_mint)
697               jm_maxt = MAX(jm_max(ji),jm_maxt)
698            END DO
699           
700            DO jk = jm_mint+1, jm_maxt
701               DO ji = 1, npti
702                  jm = MIN(MAX(jm_min(ji)+1,jk),jm_max(ji))
703                  zdiagbis(ji,jm) = ztrid   (ji,jm,2) - ztrid(ji,jm,1) * ztrid   (ji,jm-1,3) / zdiagbis(ji,jm-1)
704                  zindtbis(ji,jm) = zindterm(ji,jm)   - ztrid(ji,jm,1) * zindtbis(ji,jm-1)   / zdiagbis(ji,jm-1)
705               END DO
706            END DO
707           
708            ! ice temperatures
709            DO ji = 1, npti
710               ! Variable used after iterations
711               ! Value must be frozen after convergence for MPP independance reason
712               IF ( .NOT. l_T_converged(ji) ) &
713                  t_i_1d(ji,nlay_i) = zindtbis(ji,jm_max(ji)) / zdiagbis(ji,jm_max(ji))
714            END DO
715
716            DO jm = nlay_i + nlay_s, nlay_s + 2, -1
717               DO ji = 1, npti
718                  IF ( .NOT. l_T_converged(ji) ) THEN
719                     jk = jm - nlay_s - 1
720                     t_i_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,jm)
721                  ENDIF
722               END DO
723            END DO
724           
725            ! snow temperatures     
726            DO ji = 1, npti
727               ! Variable used after iterations
728               ! Value must be frozen after convergence for MPP independance reason
729               IF ( .NOT. l_T_converged(ji) ) THEN
730                  IF( h_s_1d(ji) > 0._wp ) THEN
731                     t_s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) ) / zdiagbis(ji,nlay_s+1)
732                  ENDIF
733               ENDIF
734            END DO
735            !
736            !--------------------------------------------------------------
737            ! 9) Has the scheme converged?, end of the iterative procedure
738            !--------------------------------------------------------------
739            ! check that nowhere it has started to melt
740            ! zdti_max is a measure of error, it has to be under zdti_bnd
741
742            DO ji = 1, npti
743
744               zdti_max = 0._wp
745
746               IF ( .NOT. l_T_converged(ji) ) THEN
747                  ! t_s
748                  t_s_1d(ji,1:nlay_s) = MAX( MIN( t_s_1d(ji,1:nlay_s), rt0 ), rt0 - 100._wp )
749                  zdti_max = MAX ( zdti_max , MAXVAL( ABS( t_s_1d(ji,1:nlay_s) - ztsb(ji,1:nlay_s) ) ) )
750                  ! t_i
751                  DO jk = 1, nlay_i
752                     ztmelts       = -rTmlt * sz_i_1d(ji,jk) + rt0 
753                     t_i_1d(ji,jk) =  MAX( MIN( t_i_1d(ji,jk), ztmelts ), rt0 - 100._wp )
754                     zdti_max      =  MAX ( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) )
755                  END DO
756
757                  IF ( zdti_max < zdti_bnd ) l_T_converged(ji) = .TRUE.
758
759               ENDIF
760
761            END DO
762
763         ENDIF ! k_cnd
764         
765      END DO  ! End of the do while iterative procedure
766     
767      IF( ln_icectl .AND. lwp ) THEN
768         WRITE(numout,*) ' zdti_max : ', zdti_max
769         WRITE(numout,*) ' iconv    : ', iconv
770      ENDIF
771     
772      !
773      !-----------------------------
774      ! 10) Fluxes at the interfaces
775      !-----------------------------
776      !
777      ! --- calculate conduction fluxes (positive downward)
778
779      DO ji = 1, npti
780         !                                ! surface ice conduction flux
781         qcn_ice_top_1d(ji) =  -           isnow(ji)   * zkappa_s(ji,0)      * zg1s * ( t_s_1d(ji,1) - t_su_1d(ji) )  &
782            &                  - ( 1._wp - isnow(ji) ) * zkappa_i(ji,0)      * zg1  * ( t_i_1d(ji,1) - t_su_1d(ji) )
783         !                                ! bottom ice conduction flux
784         qcn_ice_bot_1d(ji) =                          - zkappa_i(ji,nlay_i) * zg1  * ( t_bo_1d(ji ) - t_i_1d (ji,nlay_i) )
785      END DO
786     
787      !
788      ! --- Diagnose the heat loss due to changing non-solar / conduction flux --- !
789      !
790      IF( k_cnd == np_cnd_OFF .OR. k_cnd == np_cnd_EMU ) THEN
791         !
792         DO ji = 1, npti
793            hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) - ( qns_ice_1d(ji) - zqns_ice_b(ji) ) * a_i_1d(ji) 
794         END DO
795         !
796      ELSEIF( k_cnd == np_cnd_ON ) THEN
797         !
798         DO ji = 1, npti
799            hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) - ( qcn_ice_top_1d(ji) - qcn_ice_1d(ji) ) * a_i_1d(ji) 
800         END DO
801         !
802      ENDIF
803     
804      !
805      ! --- Diagnose the heat loss due to non-fully converged temperature solution (should not be above 10-4 W-m2) --- !
806      !
807      IF( k_cnd == np_cnd_OFF .OR. k_cnd == np_cnd_ON ) THEN
808         
809         CALL ice_var_enthalpy       
810         
811         ! zhfx_err = correction on the diagnosed heat flux due to non-convergence of the algorithm used to solve heat equation
812         DO ji = 1, npti
813            zdq = - zq_ini(ji) + ( SUM( e_i_1d(ji,1:nlay_i) ) * h_i_1d(ji) * r1_nlay_i +  &
814               &                   SUM( e_s_1d(ji,1:nlay_s) ) * h_s_1d(ji) * r1_nlay_s )
815           
816            IF( k_cnd == np_cnd_OFF ) THEN
817               
818               IF( t_su_1d(ji) < rt0 ) THEN  ! case T_su < 0degC
819                  zhfx_err = ( qns_ice_1d(ji)     + qsr_ice_1d(ji)     - zradtr_i(ji,nlay_i) - qcn_ice_bot_1d(ji)  &
820                     &       + zdq * r1_rdtice ) * a_i_1d(ji)
821               ELSE                          ! case T_su = 0degC
822                  zhfx_err = ( qcn_ice_top_1d(ji) + qtr_ice_top_1d(ji) - zradtr_i(ji,nlay_i) - qcn_ice_bot_1d(ji)  &
823                     &       + zdq * r1_rdtice ) * a_i_1d(ji)
824               ENDIF
825               
826            ELSEIF( k_cnd == np_cnd_ON ) THEN
827           
828               zhfx_err    = ( qcn_ice_top_1d(ji) + qtr_ice_top_1d(ji) - zradtr_i(ji,nlay_i) - qcn_ice_bot_1d(ji)  &
829                  &          + zdq * r1_rdtice ) * a_i_1d(ji)
830           
831            ENDIF
832            !
833            ! total heat sink to be sent to the ocean
834            hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) + zhfx_err
835            !
836            ! hfx_dif = Heat flux diagnostic of sensible heat used to warm/cool ice in W.m-2   
837            hfx_dif_1d(ji) = hfx_dif_1d(ji) - zdq * r1_rdtice * a_i_1d(ji)
838            !
839         END DO
840         !
841      ENDIF
842      !
843      !--------------------------------------------------------------------
844      ! 11) reset inner snow and ice temperatures, update conduction fluxes
845      !--------------------------------------------------------------------
846      ! effective conductivity and 1st layer temperature (needed by Met Office)
847      DO ji = 1, npti
848         IF( h_i_1d(ji) > 0.1_wp ) THEN
849            cnd_ice_1d(ji) = 2._wp * zkappa_i(ji,0)
850         ELSE
851            cnd_ice_1d(ji) = 2._wp * ztcond_i(ji,0) * 10._wp
852         ENDIF
853         t1_ice_1d(ji) = isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1)
854      END DO
855      !
856      IF( k_cnd == np_cnd_EMU ) THEN
857         ! Restore temperatures to their initial values
858         t_s_1d    (1:npti,:) = ztsold        (1:npti,:)
859         t_i_1d    (1:npti,:) = ztiold        (1:npti,:)
860         qcn_ice_1d(1:npti)   = qcn_ice_top_1d(1:npti)
861
862         !!clem
863         ! remettre t_su_1d, qns_ice_1d et dqns_ice_1d comme avant puisqu'on devrait faire comme si on avant conduction = input
864         !clem
865      ENDIF
866      !
867      ! --- SIMIP diagnostics
868      !
869      DO ji = 1, npti         
870         !--- Snow-ice interfacial temperature (diagnostic SIMIP)
871         zfac = rn_cnd_s * zh_i(ji) + ztcond_i(ji,1) * zh_s(ji)
872         IF( h_s_1d(ji) >= zhs_min ) THEN
873            t_si_1d(ji) = ( rn_cnd_s       * zh_i(ji) * t_s_1d(ji,1) +   &
874               &            ztcond_i(ji,1) * zh_s(ji) * t_i_1d(ji,1) ) / MAX( epsi10, zfac )
875         ELSE
876            t_si_1d(ji) = t_su_1d(ji)
877         ENDIF
878      END DO
879      !
880   END SUBROUTINE ice_thd_zdf_BL99
881
882#else
883   !!----------------------------------------------------------------------
884   !!   Default option       Dummy Module             No SI3 sea-ice model
885   !!----------------------------------------------------------------------
886#endif
887
888   !!======================================================================
889END MODULE icethd_zdf_BL99
Note: See TracBrowser for help on using the repository browser.