source: NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/ICE/icethd_zdf_bl99.F90 @ 12468

Last change on this file since 12468 was 12468, checked in by davestorkey, 9 months ago

2020/KERNEL-03_Storkey_Coward_RK3_stage2:

  1. Alter ABL, ICE and TOP timestepping variables to be consistent with new schema: rdt_abl —> rDt_abl rdt_ice —> rDt_ice r1_rdt_ice —> r1_Dt_ice rdttrc —> rn_Dt (always equal to ocean timestep parameter in namelist) r2dttrc —> rDt_trc (current tracer timestep)
  2. Reinstate rn_scal_load (revert previous change): rn_load —> rn_scal_load

Passes SETTE and bit compares with the trunk@12436.

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