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 @ 12919

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

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