source: NEMO/trunk/src/ICE/icethd_zdf_bl99.F90 @ 10425

Last change on this file since 10425 was 10425, checked in by smasson, 20 months ago

trunk: merge back dev_r10164_HPC09_ESIWACE_PREP_MERGE@10424 into the trunk

  • Property svn:keywords set to Id
File size: 44.1 KB
Line 
1MODULE icethd_zdf_BL99
2   !!======================================================================
3   !!                       ***  MODULE icethd_zdf_BL99 ***
4   !!   sea-ice: vertical heat diffusion in sea ice (computation of temperatures)
5   !!======================================================================
6   !! History :       !  2003-02  (M. Vancoppenolle) original 1D code
7   !!                 !  2005-06  (M. Vancoppenolle) 3d version
8   !!            4.0  !  2018     (many people)      SI3 [aka Sea Ice cube]
9   !!----------------------------------------------------------------------
10#if defined key_si3
11   !!----------------------------------------------------------------------
12   !!   'key_si3'                                       SI3 sea-ice model
13   !!----------------------------------------------------------------------
14   !!  ice_thd_zdf_BL99 : vertical diffusion computation
15   !!----------------------------------------------------------------------
16   USE dom_oce        ! ocean space and time domain
17   USE phycst         ! physical constants (ocean directory)
18   USE ice            ! sea-ice: variables
19   USE ice1D          ! sea-ice: thermodynamics variables
20   USE icevar         ! sea-ice: operations
21   !
22   USE in_out_manager ! I/O manager
23   USE lib_mpp        ! MPP library
24   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero)
25
26   IMPLICIT NONE
27   PRIVATE
28
29   PUBLIC   ice_thd_zdf_BL99   ! called by icethd_zdf
30
31   !!----------------------------------------------------------------------
32   !! NEMO/ICE 4.0 , NEMO Consortium (2018)
33   !! $Id$
34   !! Software governed by the CeCILL license (see ./LICENSE)
35   !!----------------------------------------------------------------------
36CONTAINS
37
38   SUBROUTINE ice_thd_zdf_BL99( k_jules )
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_jules     ! Jules coupling (0=OFF, 1=EMULATED, 2=ACTIVE)
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_jules == np_jules_OFF .OR. k_jules == np_jules_EMULE ) 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      !                                                          !============================!
209      ! Convergence calculated until all sub-domain grid points have converged
210      ! Calculations keep going for all grid points until sub-domain convergence (vectorisation optimisation)
211      ! but values are not taken into account (results independant of MPI partitioning)
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         SELECT CASE ( nn_virtual_itd )
268         !
269         CASE ( 1 , 2 )
270            !
271            zepsilon = 0.1_wp
272            DO ji = 1, npti
273               zcnd_i = SUM( ztcond_i(ji,:) ) / REAL( nlay_i+1, wp )                                ! Mean sea ice thermal conductivity
274               zhe = ( rn_cnd_s * h_i_1d(ji) + zcnd_i * h_s_1d(ji) ) / ( rn_cnd_s + zcnd_i )        ! Effective thickness he (zhe)
275               IF( zhe >=  zepsilon * 0.5_wp * EXP(1._wp) )  &
276                  &   zghe(ji) = MIN( 2._wp, 0.5_wp * ( 1._wp + LOG( 2._wp * zhe / zepsilon ) ) )   ! G(he)
277            END DO
278            !
279         END SELECT
280         !
281         !-----------------
282         ! 4) kappa factors
283         !-----------------
284         !--- Snow
285         ! Variable used after iterations
286         ! Value must be frozen after convergence for MPP independance reason
287         DO jk = 0, nlay_s-1
288            DO ji = 1, npti
289               IF ( .NOT. l_T_converged(ji) ) &
290                  zkappa_s(ji,jk) = zghe(ji) * rn_cnd_s * z1_h_s(ji)
291            END DO
292         END DO
293         DO ji = 1, npti   ! Snow-ice interface
294            IF ( .NOT. l_T_converged(ji) ) THEN
295               zfac = 0.5_wp * ( ztcond_i(ji,0) * zh_s(ji) + rn_cnd_s * zh_i(ji) )
296               IF( zfac > epsi10 ) THEN
297                  zkappa_s(ji,nlay_s) = zghe(ji) * rn_cnd_s * ztcond_i(ji,0) / zfac
298               ELSE
299                  zkappa_s(ji,nlay_s) = 0._wp
300               ENDIF
301            ENDIF
302         END DO
303
304         !--- Ice
305         ! Variable used after iterations
306         ! Value must be frozen after convergence for MPP independance reason
307         DO jk = 0, nlay_i
308            DO ji = 1, npti
309               IF ( .NOT. l_T_converged(ji) ) &
310                  zkappa_i(ji,jk) = zghe(ji) * ztcond_i(ji,jk) * z1_h_i(ji)
311            END DO
312         END DO
313         DO ji = 1, npti   ! Snow-ice interface
314            IF ( .NOT. l_T_converged(ji) ) &
315               zkappa_i(ji,0) = zkappa_s(ji,nlay_s) * isnow(ji) + zkappa_i(ji,0) * ( 1._wp - isnow(ji) )
316         END DO
317         !
318         !--------------------------------------
319         ! 5) Sea ice specific heat, eta factors
320         !--------------------------------------
321         DO jk = 1, nlay_i
322            DO ji = 1, npti
323               zcpi = rcpi + zgamma * sz_i_1d(ji,jk) / MAX( ( t_i_1d(ji,jk) - rt0 ) * ( ztiold(ji,jk) - rt0 ), epsi10 )
324               zeta_i(ji,jk) = rdt_ice * r1_rhoi * z1_h_i(ji) / MAX( epsi10, zcpi ) 
325            END DO
326         END DO
327
328         DO jk = 1, nlay_s
329            DO ji = 1, npti
330               zeta_s(ji,jk) = rdt_ice * r1_rhos * r1_rcpi * z1_h_s(ji)
331            END DO
332         END DO
333         !
334         !----------------------------------------!
335         !                                        !
336         !   JULES COUPLING IS OFF OR EMULATED    !
337         !                                        !
338         !----------------------------------------!
339         !
340         IF( k_jules == np_jules_OFF .OR. k_jules == np_jules_EMULE ) THEN
341            !
342            ! ==> The original BL99 temperature computation is used
343            !       (with qsr_ice, qns_ice and dqns_ice as inputs)
344            !
345            !----------------------------
346            ! 6) surface flux computation
347            !----------------------------
348            ! update of the non solar flux according to the update in T_su
349            DO ji = 1, npti
350               ! Variable used after iterations
351               ! Value must be frozen after convergence for MPP independance reason
352               IF ( .NOT. l_T_converged(ji) ) &
353                  qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_1d(ji) - ztsub(ji) )
354            END DO
355
356            DO ji = 1, npti
357               zfnet(ji) = qsr_ice_1d(ji) - qtr_ice_top_1d(ji) + qns_ice_1d(ji) ! net heat flux = net - transmitted solar + non solar
358            END DO
359            !
360            !----------------------------
361            ! 7) tridiagonal system terms
362            !----------------------------
363            ! layer denotes the number of the layer in the snow or in the ice
364            ! jm denotes the reference number of the equation in the tridiagonal
365            ! system, terms of tridiagonal system are indexed as following :
366            ! 1 is subdiagonal term, 2 is diagonal and 3 is superdiagonal one
367
368            ! ice interior terms (top equation has the same form as the others)
369            ztrid   (1:npti,:,:) = 0._wp
370            zindterm(1:npti,:)   = 0._wp
371            zindtbis(1:npti,:)   = 0._wp
372            zdiagbis(1:npti,:)   = 0._wp
373
374            DO jm = nlay_s + 2, nlay_s + nlay_i 
375               DO ji = 1, npti
376                  jk = jm - nlay_s - 1
377                  ztrid   (ji,jm,1) =       - zeta_i(ji,jk) *   zkappa_i(ji,jk-1)
378                  ztrid   (ji,jm,2) = 1._wp + zeta_i(ji,jk) * ( zkappa_i(ji,jk-1) + zkappa_i(ji,jk) )
379                  ztrid   (ji,jm,3) =       - zeta_i(ji,jk) *                       zkappa_i(ji,jk)
380                  zindterm(ji,jm)   = ztiold(ji,jk) + zeta_i(ji,jk) * zradab_i(ji,jk)
381               END DO
382            END DO
383
384            jm =  nlay_s + nlay_i + 1
385            DO ji = 1, npti
386               ! ice bottom term
387               ztrid   (ji,jm,1) =       - zeta_i(ji,nlay_i) *   zkappa_i(ji,nlay_i-1)   
388               ztrid   (ji,jm,2) = 1._wp + zeta_i(ji,nlay_i) * ( zkappa_i(ji,nlay_i-1) + zkappa_i(ji,nlay_i) * zg1 )
389               ztrid   (ji,jm,3) = 0._wp
390               zindterm(ji,jm)   = ztiold(ji,nlay_i) + zeta_i(ji,nlay_i) *  &
391                  &              ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i) * zg1 * t_bo_1d(ji) ) 
392            END DO
393
394            DO ji = 1, npti
395               !                               !---------------------!
396               IF( h_s_1d(ji) > 0._wp ) THEN   !  snow-covered cells !
397                  !                            !---------------------!
398                  ! snow interior terms (bottom equation has the same form as the others)
399                  DO jm = 3, nlay_s + 1
400                     jk = jm - 1
401                     ztrid   (ji,jm,1) =       - zeta_s(ji,jk) *   zkappa_s(ji,jk-1)
402                     ztrid   (ji,jm,2) = 1._wp + zeta_s(ji,jk) * ( zkappa_s(ji,jk-1) + zkappa_s(ji,jk) )
403                     ztrid   (ji,jm,3) =       - zeta_s(ji,jk) *                       zkappa_s(ji,jk)
404                     zindterm(ji,jm)   = ztsold(ji,jk) + zeta_s(ji,jk) * zradab_s(ji,jk)
405                  END DO
406                 
407                  ! case of only one layer in the ice (ice equation is altered)
408                  IF( nlay_i == 1 ) THEN
409                     ztrid   (ji,nlay_s+2,3) = 0._wp
410                     zindterm(ji,nlay_s+2)   = zindterm(ji,nlay_s+2) + zeta_i(ji,1) * zkappa_i(ji,1) * t_bo_1d(ji) 
411                  ENDIF
412                 
413                  IF( t_su_1d(ji) < rt0 ) THEN   !--  case 1 : no surface melting
414                     
415                     jm_min(ji) = 1
416                     jm_max(ji) = nlay_i + nlay_s + 1
417                     
418                     ! surface equation
419                     ztrid   (ji,1,1) = 0._wp
420                     ztrid   (ji,1,2) = zdqns_ice_b(ji) - zg1s * zkappa_s(ji,0)
421                     ztrid   (ji,1,3) =                   zg1s * zkappa_s(ji,0)
422                     zindterm(ji,1)   = zdqns_ice_b(ji) * t_su_1d(ji) - zfnet(ji)
423                     
424                     ! first layer of snow equation
425                     ztrid   (ji,2,1) =       - zeta_s(ji,1) *                    zkappa_s(ji,0) * zg1s
426                     ztrid   (ji,2,2) = 1._wp + zeta_s(ji,1) * ( zkappa_s(ji,1) + zkappa_s(ji,0) * zg1s )
427                     ztrid   (ji,2,3) =       - zeta_s(ji,1) *   zkappa_s(ji,1)
428                     zindterm(ji,2)   = ztsold(ji,1) + zeta_s(ji,1) * zradab_s(ji,1)
429                     
430                  ELSE                            !--  case 2 : surface is melting
431                     !
432                     jm_min(ji) = 2
433                     jm_max(ji) = nlay_i + nlay_s + 1
434                     
435                     ! first layer of snow equation
436                     ztrid   (ji,2,1) = 0._wp
437                     ztrid   (ji,2,2) = 1._wp + zeta_s(ji,1) * ( zkappa_s(ji,1) + zkappa_s(ji,0) * zg1s )
438                     ztrid   (ji,2,3) =       - zeta_s(ji,1) *   zkappa_s(ji,1) 
439                     zindterm(ji,2)   = ztsold(ji,1) + zeta_s(ji,1) * ( zradab_s(ji,1) + zkappa_s(ji,0) * zg1s * t_su_1d(ji) ) 
440                  ENDIF
441                  !                            !---------------------!
442               ELSE                            ! cells without snow  !
443                  !                            !---------------------!
444                  !
445                  IF( t_su_1d(ji) < rt0 ) THEN   !--  case 1 : no surface melting
446                     !
447                     jm_min(ji) = nlay_s + 1
448                     jm_max(ji) = nlay_i + nlay_s + 1
449                     
450                     ! surface equation   
451                     ztrid   (ji,jm_min(ji),1) = 0._wp
452                     ztrid   (ji,jm_min(ji),2) = zdqns_ice_b(ji) - zkappa_i(ji,0) * zg1   
453                     ztrid   (ji,jm_min(ji),3) =                   zkappa_i(ji,0) * zg1
454                     zindterm(ji,jm_min(ji))   = zdqns_ice_b(ji) * t_su_1d(ji) - zfnet(ji)
455                     
456                     ! first layer of ice equation
457                     ztrid   (ji,jm_min(ji)+1,1) =       - zeta_i(ji,1) *                    zkappa_i(ji,0) * zg1
458                     ztrid   (ji,jm_min(ji)+1,2) = 1._wp + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 )
459                     ztrid   (ji,jm_min(ji)+1,3) =       - zeta_i(ji,1) *   zkappa_i(ji,1) 
460                     zindterm(ji,jm_min(ji)+1)   = ztiold(ji,1) + zeta_i(ji,1) * zradab_i(ji,1) 
461                     
462                     ! case of only one layer in the ice (surface & ice equations are altered)
463                     IF( nlay_i == 1 ) THEN
464                        ztrid   (ji,jm_min(ji),1)   = 0._wp
465                        ztrid   (ji,jm_min(ji),2)   = zdqns_ice_b(ji)      -   zkappa_i(ji,0) * 2._wp
466                        ztrid   (ji,jm_min(ji),3)   =                          zkappa_i(ji,0) * 2._wp
467                        ztrid   (ji,jm_min(ji)+1,1) =       - zeta_i(ji,1) *   zkappa_i(ji,0) * 2._wp
468                        ztrid   (ji,jm_min(ji)+1,2) = 1._wp + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2._wp + zkappa_i(ji,1) )
469                        ztrid   (ji,jm_min(ji)+1,3) = 0._wp
470                        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))
471                     ENDIF
472                     
473                  ELSE                            !--  case 2 : surface is melting
474                     
475                     jm_min(ji) = nlay_s + 2
476                     jm_max(ji) = nlay_i + nlay_s + 1
477                     
478                     ! first layer of ice equation
479                     ztrid   (ji,jm_min(ji),1) = 0._wp
480                     ztrid   (ji,jm_min(ji),2) = 1._wp + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 ) 
481                     ztrid   (ji,jm_min(ji),3) =       - zeta_i(ji,1) *   zkappa_i(ji,1)
482                     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)) 
483                     
484                     ! case of only one layer in the ice (surface & ice equations are altered)
485                     IF( nlay_i == 1 ) THEN
486                        ztrid   (ji,jm_min(ji),1) = 0._wp
487                        ztrid   (ji,jm_min(ji),2) = 1._wp + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2._wp + zkappa_i(ji,1) )
488                        ztrid   (ji,jm_min(ji),3) = 0._wp
489                        zindterm(ji,jm_min(ji))   = ztiold(ji,1) + zeta_i(ji,1) * ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) ) &
490                           &                      + t_su_1d(ji) * zeta_i(ji,1) * zkappa_i(ji,0) * 2._wp
491                     ENDIF
492                     
493                  ENDIF
494               ENDIF
495               !
496               zindtbis(ji,jm_min(ji)) = zindterm(ji,jm_min(ji))
497               zdiagbis(ji,jm_min(ji)) = ztrid   (ji,jm_min(ji),2)
498               !
499            END DO
500            !
501            !------------------------------
502            ! 8) tridiagonal system solving
503            !------------------------------
504            ! Solve the tridiagonal system with Gauss elimination method.
505            ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, McGraw-Hill 1984
506            jm_maxt = 0
507            jm_mint = nlay_i+5
508            DO ji = 1, npti
509               jm_mint = MIN(jm_min(ji),jm_mint)
510               jm_maxt = MAX(jm_max(ji),jm_maxt)
511            END DO
512
513            DO jk = jm_mint+1, jm_maxt
514               DO ji = 1, npti
515                  jm = MIN(MAX(jm_min(ji)+1,jk),jm_max(ji))
516                  zdiagbis(ji,jm) = ztrid   (ji,jm,2) - ztrid(ji,jm,1) * ztrid   (ji,jm-1,3) / zdiagbis(ji,jm-1)
517                  zindtbis(ji,jm) = zindterm(ji,jm  ) - ztrid(ji,jm,1) * zindtbis(ji,jm-1  ) / zdiagbis(ji,jm-1)
518               END DO
519            END DO
520
521            ! ice temperatures
522            DO ji = 1, npti
523               ! Variable used after iterations
524               ! Value must be frozen after convergence for MPP independance reason
525               IF ( .NOT. l_T_converged(ji) ) &
526                  t_i_1d(ji,nlay_i) = zindtbis(ji,jm_max(ji)) / zdiagbis(ji,jm_max(ji))
527            END DO
528
529            DO jm = nlay_i + nlay_s, nlay_s + 2, -1
530               DO ji = 1, npti
531                  jk = jm - nlay_s - 1
532                  IF ( .NOT. l_T_converged(ji) ) &
533                     t_i_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,jm)
534               END DO
535            END DO
536
537            DO ji = 1, npti
538               ! Variables used after iterations
539               ! Value must be frozen after convergence for MPP independance reason
540               IF ( .NOT. l_T_converged(ji) ) THEN
541                  ! snow temperatures     
542                  IF( h_s_1d(ji) > 0._wp ) THEN
543                     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)
544                  ENDIF
545                  ! surface temperature
546                  ztsub(ji) = t_su_1d(ji)
547                  IF( t_su_1d(ji) < rt0 ) THEN
548                     t_su_1d(ji) = ( zindtbis(ji,jm_min(ji)) - ztrid(ji,jm_min(ji),3) *  &
549                        &          ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) ) ) / zdiagbis(ji,jm_min(ji))
550                  ENDIF
551               ENDIF
552            END DO
553            !
554            !--------------------------------------------------------------
555            ! 9) Has the scheme converged?, end of the iterative procedure
556            !--------------------------------------------------------------
557            ! check that nowhere it has started to melt
558            ! zdti_max is a measure of error, it has to be under zdti_bnd
559
560            DO ji = 1, npti
561
562               zdti_max = 0._wp
563
564               IF ( .NOT. l_T_converged(ji) ) THEN
565                  t_su_1d(ji) = MAX( MIN( t_su_1d(ji) , rt0 ) , rt0 - 100._wp )
566                  zdti_max    = MAX( zdti_max, ABS( t_su_1d(ji) - ztsub(ji) ) )
567
568                  t_s_1d(ji,1:nlay_s) = MAX( MIN( t_s_1d(ji,1:nlay_s), rt0 ), rt0 - 100._wp )
569                  zdti_max = MAX ( zdti_max , MAXVAL( ABS( t_s_1d(ji,1:nlay_s) - ztsb(ji,1:nlay_s) ) ) )
570
571                  DO jk = 1, nlay_i
572                     ztmelts       = -rTmlt * sz_i_1d(ji,jk) + rt0
573                     t_i_1d(ji,jk) =  MAX( MIN( t_i_1d(ji,jk), ztmelts ), rt0 - 100._wp )
574                     zdti_max      =  MAX( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) )
575                  END DO
576
577                  IF ( zdti_max < zdti_bnd ) l_T_converged(ji) = .TRUE.
578
579               ENDIF
580
581            END DO
582
583         !----------------------------------------!
584         !                                        !
585         !      JULES COUPLING IS ACTIVE          !
586         !                                        !
587         !----------------------------------------!
588         !
589         ELSEIF( k_jules == np_jules_ACTIVE ) THEN
590            !
591            ! ==> we use a modified BL99 solver with conduction flux (qcn_ice) as forcing term
592            !
593            !----------------------------
594            ! 7) tridiagonal system terms
595            !----------------------------
596            ! layer denotes the number of the layer in the snow or in the ice
597            ! jm denotes the reference number of the equation in the tridiagonal
598            ! system, terms of tridiagonal system are indexed as following :
599            ! 1 is subdiagonal term, 2 is diagonal and 3 is superdiagonal one
600
601            ! ice interior terms (top equation has the same form as the others)
602            ztrid   (1:npti,:,:) = 0._wp
603            zindterm(1:npti,:)   = 0._wp
604            zindtbis(1:npti,:)   = 0._wp
605            zdiagbis(1:npti,:)   = 0._wp
606
607            DO jm = nlay_s + 2, nlay_s + nlay_i 
608               DO ji = 1, npti
609                  jk = jm - nlay_s - 1
610                  ztrid   (ji,jm,1) =       - zeta_i(ji,jk) *   zkappa_i(ji,jk-1)
611                  ztrid   (ji,jm,2) = 1._wp + zeta_i(ji,jk) * ( zkappa_i(ji,jk-1) + zkappa_i(ji,jk) )
612                  ztrid   (ji,jm,3) =       - zeta_i(ji,jk) *                       zkappa_i(ji,jk)
613                  zindterm(ji,jm)   = ztiold(ji,jk) + zeta_i(ji,jk) * zradab_i(ji,jk)
614               END DO
615            ENDDO
616
617            jm =  nlay_s + nlay_i + 1
618            DO ji = 1, npti
619               ! ice bottom term
620               ztrid   (ji,jm,1) =       - zeta_i(ji,nlay_i) *   zkappa_i(ji,nlay_i-1)   
621               ztrid   (ji,jm,2) = 1._wp + zeta_i(ji,nlay_i) * ( zkappa_i(ji,nlay_i-1) + zkappa_i(ji,nlay_i) * zg1 )
622               ztrid   (ji,jm,3) = 0._wp
623               zindterm(ji,jm)   = ztiold(ji,nlay_i) + zeta_i(ji,nlay_i) *  &
624                  &              ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i) * zg1 * t_bo_1d(ji) ) 
625            ENDDO
626
627            DO ji = 1, npti
628               !                               !---------------------!
629               IF( h_s_1d(ji) > 0._wp ) THEN   !  snow-covered cells !
630                  !                            !---------------------!
631                  ! snow interior terms (bottom equation has the same form as the others)
632                  DO jm = 3, nlay_s + 1
633                     jk = jm - 1
634                     ztrid   (ji,jm,1) =       - zeta_s(ji,jk) *   zkappa_s(ji,jk-1)
635                     ztrid   (ji,jm,2) = 1._wp + zeta_s(ji,jk) * ( zkappa_s(ji,jk-1) + zkappa_s(ji,jk) )
636                     ztrid   (ji,jm,3) =       - zeta_s(ji,jk) *                       zkappa_s(ji,jk)
637                     zindterm(ji,jm)   = ztsold(ji,jk) + zeta_s(ji,jk) * zradab_s(ji,jk)
638                  END DO
639                 
640                  ! case of only one layer in the ice (ice equation is altered)
641                  IF ( nlay_i == 1 ) THEN
642                     ztrid   (ji,nlay_s+2,3) = 0._wp
643                     zindterm(ji,nlay_s+2)   = zindterm(ji,nlay_s+2) + zeta_i(ji,1) * zkappa_i(ji,1) * t_bo_1d(ji) 
644                  ENDIF
645                 
646                  jm_min(ji) = 2
647                  jm_max(ji) = nlay_i + nlay_s + 1
648                 
649                  ! first layer of snow equation
650                  ztrid   (ji,2,1) = 0._wp
651                  ztrid   (ji,2,2) = 1._wp + zeta_s(ji,1) * zkappa_s(ji,1)
652                  ztrid   (ji,2,3) =       - zeta_s(ji,1) * zkappa_s(ji,1) 
653                  zindterm(ji,2)   = ztsold(ji,1) + zeta_s(ji,1) * ( zradab_s(ji,1) + qcn_ice_1d(ji) ) 
654                 
655                  !                            !---------------------!
656               ELSE                            ! cells without snow  !
657                  !                            !---------------------!
658                  jm_min(ji) = nlay_s + 2
659                  jm_max(ji) = nlay_i + nlay_s + 1
660                 
661                  ! first layer of ice equation
662                  ztrid   (ji,jm_min(ji),1) = 0._wp
663                  ztrid   (ji,jm_min(ji),2) = 1._wp + zeta_i(ji,1) * zkappa_i(ji,1) 
664                  ztrid   (ji,jm_min(ji),3) =       - zeta_i(ji,1) * zkappa_i(ji,1)
665                  zindterm(ji,jm_min(ji))   = ztiold(ji,1) + zeta_i(ji,1) * ( zradab_i(ji,1) + qcn_ice_1d(ji) )
666                 
667                  ! case of only one layer in the ice (surface & ice equations are altered)
668                  IF( nlay_i == 1 ) THEN
669                     ztrid   (ji,jm_min(ji),1) = 0._wp
670                     ztrid   (ji,jm_min(ji),2) = 1._wp + zeta_i(ji,1) * zkappa_i(ji,1)
671                     ztrid   (ji,jm_min(ji),3) = 0._wp
672                     zindterm(ji,jm_min(ji))   = ztiold(ji,1) + zeta_i(ji,1) *  &
673                        &                                     ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) + qcn_ice_1d(ji) )
674                  ENDIF
675                 
676               ENDIF
677               !
678               zindtbis(ji,jm_min(ji)) = zindterm(ji,jm_min(ji))
679               zdiagbis(ji,jm_min(ji)) = ztrid   (ji,jm_min(ji),2)
680               !
681            END DO
682            !
683            !------------------------------
684            ! 8) tridiagonal system solving
685            !------------------------------
686            ! Solve the tridiagonal system with Gauss elimination method.
687            ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, McGraw-Hill 1984
688            jm_maxt = 0
689            jm_mint = nlay_i+5
690            DO ji = 1, npti
691               jm_mint = MIN(jm_min(ji),jm_mint)
692               jm_maxt = MAX(jm_max(ji),jm_maxt)
693            END DO
694           
695            DO jk = jm_mint+1, jm_maxt
696               DO ji = 1, npti
697                  jm = MIN(MAX(jm_min(ji)+1,jk),jm_max(ji))
698                  zdiagbis(ji,jm) = ztrid   (ji,jm,2) - ztrid(ji,jm,1) * ztrid   (ji,jm-1,3) / zdiagbis(ji,jm-1)
699                  zindtbis(ji,jm) = zindterm(ji,jm)   - ztrid(ji,jm,1) * zindtbis(ji,jm-1)   / zdiagbis(ji,jm-1)
700               END DO
701            END DO
702           
703            ! ice temperatures
704            DO ji = 1, npti
705               ! Variable used after iterations
706               ! Value must be frozen after convergence for MPP independance reason
707               IF ( .NOT. l_T_converged(ji) ) &
708                  t_i_1d(ji,nlay_i) = zindtbis(ji,jm_max(ji)) / zdiagbis(ji,jm_max(ji))
709            END DO
710
711            DO jm = nlay_i + nlay_s, nlay_s + 2, -1
712               DO ji = 1, npti
713                  IF ( .NOT. l_T_converged(ji) ) THEN
714                     jk = jm - nlay_s - 1
715                     t_i_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,jm)
716                  ENDIF
717               END DO
718            END DO
719           
720            ! snow 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) ) THEN
725                  IF( h_s_1d(ji) > 0._wp ) THEN
726                     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)
727                  ENDIF
728               ENDIF
729            END DO
730            !
731            !--------------------------------------------------------------
732            ! 9) Has the scheme converged?, end of the iterative procedure
733            !--------------------------------------------------------------
734            ! check that nowhere it has started to melt
735            ! zdti_max is a measure of error, it has to be under zdti_bnd
736
737            DO ji = 1, npti
738
739               zdti_max = 0._wp
740
741               IF ( .NOT. l_T_converged(ji) ) THEN
742                  ! t_s
743                  t_s_1d(ji,1:nlay_s) = MAX( MIN( t_s_1d(ji,1:nlay_s), rt0 ), rt0 - 100._wp )
744                  zdti_max = MAX ( zdti_max , MAXVAL( ABS( t_s_1d(ji,1:nlay_s) - ztsb(ji,1:nlay_s) ) ) )
745                  ! t_i
746                  DO jk = 0, nlay_i
747                     ztmelts       = -rTmlt * sz_i_1d(ji,jk) + rt0 
748                     t_i_1d(ji,jk) =  MAX( MIN( t_i_1d(ji,jk), ztmelts ), rt0 - 100._wp )
749                     zdti_max      =  MAX ( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) )
750                  END DO
751
752                  IF ( zdti_max < zdti_bnd ) l_T_converged(ji) = .TRUE.
753
754               ENDIF
755
756            END DO
757
758         ENDIF ! k_jules
759         
760      END DO  ! End of the do while iterative procedure
761     
762      IF( ln_icectl .AND. lwp ) THEN
763         WRITE(numout,*) ' zdti_max : ', zdti_max
764         WRITE(numout,*) ' iconv    : ', iconv
765      ENDIF
766     
767      !
768      !-----------------------------
769      ! 10) Fluxes at the interfaces
770      !-----------------------------
771      !
772      ! --- calculate conduction fluxes (positive downward)
773
774      DO ji = 1, npti
775         !                                ! surface ice conduction flux
776         qcn_ice_top_1d(ji) =  -           isnow(ji)   * zkappa_s(ji,0)      * zg1s * ( t_s_1d(ji,1) - t_su_1d(ji) )  &
777            &                  - ( 1._wp - isnow(ji) ) * zkappa_i(ji,0)      * zg1  * ( t_i_1d(ji,1) - t_su_1d(ji) )
778         !                                ! bottom ice conduction flux
779         qcn_ice_bot_1d(ji) =                          - zkappa_i(ji,nlay_i) * zg1  * ( t_bo_1d(ji ) - t_i_1d (ji,nlay_i) )
780      END DO
781     
782      !
783      ! --- Diagnose the heat loss due to changing non-solar / conduction flux --- !
784      !
785      IF( k_jules == np_jules_OFF .OR. k_jules == np_jules_EMULE ) THEN
786         !
787         DO ji = 1, npti
788            hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) - ( qns_ice_1d(ji) - zqns_ice_b(ji) ) * a_i_1d(ji) 
789         END DO
790         !
791      ELSEIF( k_jules == np_jules_ACTIVE ) THEN
792         !
793         DO ji = 1, npti
794            hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) - ( qcn_ice_top_1d(ji) - qcn_ice_1d(ji) ) * a_i_1d(ji) 
795         END DO
796         !
797      ENDIF
798     
799      !
800      ! --- Diagnose the heat loss due to non-fully converged temperature solution (should not be above 10-4 W-m2) --- !
801      !
802      IF( k_jules == np_jules_OFF .OR. k_jules == np_jules_ACTIVE ) THEN
803         
804         CALL ice_var_enthalpy       
805         
806         ! zhfx_err = correction on the diagnosed heat flux due to non-convergence of the algorithm used to solve heat equation
807         DO ji = 1, npti
808            zdq = - zq_ini(ji) + ( SUM( e_i_1d(ji,1:nlay_i) ) * h_i_1d(ji) * r1_nlay_i +  &
809               &                   SUM( e_s_1d(ji,1:nlay_s) ) * h_s_1d(ji) * r1_nlay_s )
810           
811            IF( k_jules == np_jules_OFF ) THEN
812               
813               IF( t_su_1d(ji) < rt0 ) THEN  ! case T_su < 0degC
814                  zhfx_err = ( qns_ice_1d(ji)     + qsr_ice_1d(ji)     - zradtr_i(ji,nlay_i) - qcn_ice_bot_1d(ji)  &
815                     &       + zdq * r1_rdtice ) * a_i_1d(ji)
816               ELSE                          ! case T_su = 0degC
817                  zhfx_err = ( qcn_ice_top_1d(ji) + qtr_ice_top_1d(ji) - zradtr_i(ji,nlay_i) - qcn_ice_bot_1d(ji)  &
818                     &       + zdq * r1_rdtice ) * a_i_1d(ji)
819               ENDIF
820               
821            ELSEIF( k_jules == np_jules_ACTIVE ) THEN
822           
823               zhfx_err    = ( qcn_ice_top_1d(ji) + qtr_ice_top_1d(ji) - zradtr_i(ji,nlay_i) - qcn_ice_bot_1d(ji)  &
824                  &          + zdq * r1_rdtice ) * a_i_1d(ji)
825           
826            ENDIF
827            !
828            ! total heat sink to be sent to the ocean
829            hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) + zhfx_err
830            !
831            ! hfx_dif = Heat flux diagnostic of sensible heat used to warm/cool ice in W.m-2   
832            hfx_dif_1d(ji) = hfx_dif_1d(ji) - zdq * r1_rdtice * a_i_1d(ji)
833            !
834         END DO
835         !
836      ENDIF
837      !
838      !---------------------------------------------------------------------------------------
839      ! 11) Jules coupling: reset inner snow and ice temperatures, update conduction fluxes
840      !---------------------------------------------------------------------------------------
841      ! effective conductivity and 1st layer temperature (needed by Met Office)
842      DO ji = 1, npti
843         IF( h_s_1d(ji) > 0.1_wp ) THEN
844            cnd_ice_1d(ji) = 2._wp * zkappa_s(ji,0)
845         ELSE
846            IF( h_i_1d(ji) > 0.1_wp ) THEN
847               cnd_ice_1d(ji) = 2._wp * zkappa_i(ji,0)
848            ELSE
849               cnd_ice_1d(ji) = 2._wp * ztcond_i(ji,0) * 10._wp
850            ENDIF
851         ENDIF
852         t1_ice_1d(ji) = isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1)
853      END DO
854      !
855      IF( k_jules == np_jules_EMULE ) THEN
856         ! Restore temperatures to their initial values
857         t_s_1d    (1:npti,:) = ztsold        (1:npti,:)
858         t_i_1d    (1:npti,:) = ztiold        (1:npti,:)
859         qcn_ice_1d(1:npti)   = qcn_ice_top_1d(1:npti)
860      ENDIF
861      !
862      ! --- SIMIP diagnostics
863      !
864      DO ji = 1, npti         
865         !--- Snow-ice interfacial temperature (diagnostic SIMIP)
866         zfac = rn_cnd_s * zh_i(ji) + ztcond_i(ji,1) * zh_s(ji)
867         IF( h_s_1d(ji) >= zhs_min ) THEN
868            t_si_1d(ji) = ( rn_cnd_s       * zh_i(ji) * t_s_1d(ji,1) +   &
869               &            ztcond_i(ji,1) * zh_s(ji) * t_i_1d(ji,1) ) / MAX( epsi10, zfac )
870         ELSE
871            t_si_1d(ji) = t_su_1d(ji)
872         ENDIF
873      END DO
874      !
875   END SUBROUTINE ice_thd_zdf_BL99
876
877#else
878   !!----------------------------------------------------------------------
879   !!   Default option       Dummy Module             No SI3 sea-ice model
880   !!----------------------------------------------------------------------
881#endif
882
883   !!======================================================================
884END MODULE icethd_zdf_BL99
Note: See TracBrowser for help on using the repository browser.