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

Last change on this file since 12995 was 12995, checked in by dancopsey, 17 months ago

Do not use isnow in forming zkappa_i but still use it for the coupler

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