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

Last change on this file since 10069 was 10069, checked in by nicolasmartin, 2 years ago

Fix mistakes of previous commit on SVN keywords property

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