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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
icethd_zdf_bl99.F90 in NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE – NEMO

source: NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icethd_zdf_bl99.F90 @ 12854

Last change on this file since 12854 was 12854, checked in by clem, 4 years ago

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