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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
icethd_zdf_bl99.F90 in branches/2017/dev_merge_2017/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/2017/dev_merge_2017/NEMOGCM/NEMO/LIM_SRC_3/icethd_zdf_bl99.F90 @ 9019

Last change on this file since 9019 was 9019, checked in by timgraham, 6 years ago

Merge of dev_CNRS_2017 into branch

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 :  LIM  ! 02-2003 (M. Vancoppenolle) original 1D code
7   !!                 ! 06-2005 (M. Vancoppenolle) 3d version
8   !!                 ! 11-2006 (X Fettweis) Vectorization by Xavier
9   !!                 ! 04-2007 (M. Vancoppenolle) Energy conservation
10   !!            4.0  ! 2011-02 (G. Madec) dynamical allocation
11   !!             -   ! 2012-05 (C. Rousset) add penetration solar flux
12   !!----------------------------------------------------------------------
13#if defined key_lim3
14   !!----------------------------------------------------------------------
15   !!   'key_lim3'                                       ESIM sea-ice model
16   !!----------------------------------------------------------------------
17   !!  ice_thd_zdf_BL99 : vertical diffusion computation
18   !!----------------------------------------------------------------------
19   USE dom_oce        ! ocean space and time domain
20   USE phycst         ! physical constants (ocean directory)
21   USE ice            ! sea-ice: variables
22   USE ice1D          ! sea-ice: thermodynamics variables
23   USE icevar         ! sea-ice: operations
24   !
25   USE in_out_manager ! I/O manager
26   USE lib_mpp        ! MPP library
27   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero)
28
29   IMPLICIT NONE
30   PRIVATE
31
32   PUBLIC   ice_thd_zdf_BL99   ! called by icethd_zdf
33
34   !!----------------------------------------------------------------------
35   !! NEMO/ICE 4.0 , NEMO Consortium (2017)
36   !! $Id: icethd_zdf.F90 8420 2017-08-08 12:18:46Z clem $
37   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
38   !!----------------------------------------------------------------------
39CONTAINS
40
41   SUBROUTINE ice_thd_zdf_BL99( k_jules )
42      !!-------------------------------------------------------------------
43      !!                ***  ROUTINE ice_thd_zdf_BL99  ***
44      !!
45      !! ** Purpose : computes the time evolution of snow and sea-ice temperature
46      !!              profiles, using the original Bitz and Lipscomb (1999) algorithm
47      !!
48      !! ** Method  : solves the heat equation diffusion with a Neumann boundary
49      !!              condition at the surface and a Dirichlet one at the bottom.
50      !!              Solar radiation is partially absorbed into the ice.
51      !!              The specific heat and thermal conductivities depend on ice
52      !!              salinity and temperature to take into account brine pocket   
53      !!              melting. The numerical scheme is an iterative Crank-Nicolson
54      !!              on a non-uniform multilayer grid in the ice and snow system.
55      !!
56      !!           The successive steps of this routine are
57      !!           1.  initialization of ice-snow layers thicknesses
58      !!           2.  Internal absorbed and transmitted radiation
59      !!           Then iterative procedure begins
60      !!           3.  Thermal conductivity
61      !!           4.  Kappa factors
62      !!           5.  specific heat in the ice
63      !!           6.  eta factors
64      !!           7.  surface flux computation
65      !!           8.  tridiagonal system terms
66      !!           9.  solving the tridiagonal system with Gauss elimination
67      !!           Iterative procedure ends according to a criterion on evolution
68      !!           of temperature
69      !!           10. Fluxes at the interfaces
70      !!
71      !! ** Inputs / Ouputs : (global commons)
72      !!           surface temperature              : t_su_1d
73      !!           ice/snow temperatures            : t_i_1d, t_s_1d
74      !!           ice salinities                   : sz_i_1d
75      !!           number of layers in the ice/snow : nlay_i, nlay_s
76      !!           total ice/snow thickness         : h_i_1d, h_s_1d
77      !!-------------------------------------------------------------------
78      INTEGER, INTENT(in) ::   k_jules     ! Jules coupling (0=OFF, 1=EMULATED, 2=ACTIVE)
79      !
80      INTEGER ::   ji, jk         ! spatial loop index
81      INTEGER ::   jm             ! current reference number of equation
82      INTEGER ::   jm_mint, jm_maxt
83      INTEGER ::   iconv          ! number of iterations in iterative procedure
84      INTEGER ::   iconv_max = 50 ! max number of iterations in iterative procedure
85      !
86      INTEGER, DIMENSION(jpij) ::   jm_min    ! reference number of top equation
87      INTEGER, DIMENSION(jpij) ::   jm_max    ! reference number of bottom equation
88      !
89      REAL(wp) ::   zg1s      =  2._wp        ! for the tridiagonal system
90      REAL(wp) ::   zg1       =  2._wp        !
91      REAL(wp) ::   zgamma    =  18009._wp    ! for specific heat
92      REAL(wp) ::   zbeta     =  0.117_wp     ! for thermal conductivity (could be 0.13)
93      REAL(wp) ::   zraext_s  =  10._wp       ! extinction coefficient of radiation in the snow
94      REAL(wp) ::   zkimin    =  0.10_wp      ! minimum ice thermal conductivity
95      REAL(wp) ::   ztsu_err  =  1.e-5_wp     ! range around which t_su is considered at 0C
96      REAL(wp) ::   zdti_bnd  =  1.e-4_wp     ! maximal authorized error on temperature
97      REAL(wp) ::   ztmelt_i                  ! 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)   ::   zradtr_i    ! Radiation transmitted through the ice
118      REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   zradab_i    ! Radiation absorbed in the ice
119      REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   zkappa_i    ! Kappa factor in the ice
120      REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   zeta_i      ! Eta factor in the ice
121      REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   zradtr_s    ! Radiation transmited through the snow
122      REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   zradab_s    ! Radiation absorbed in the snow
123      REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   zkappa_s    ! Kappa factor in the snow
124      REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   zeta_s      ! Eta factor in the snow
125      REAL(wp), DIMENSION(jpij,nlay_i+3)   ::   zindterm    ! 'Ind'ependent term
126      REAL(wp), DIMENSION(jpij,nlay_i+3)   ::   zindtbis    ! Temporary 'ind'ependent term
127      REAL(wp), DIMENSION(jpij,nlay_i+3)   ::   zdiagbis    ! Temporary 'dia'gonal term
128      REAL(wp), DIMENSION(jpij,nlay_i+3,3) ::   ztrid       ! Tridiagonal system terms
129      REAL(wp), DIMENSION(jpij)            ::   zq_ini      ! diag errors on heat
130      REAL(wp), DIMENSION(jpij)            ::   zghe        ! G(he), th. conduct enhancement factor, mono-cat
131      !
132      ! Mono-category
133      REAL(wp) ::   zepsilon   ! determines thres. above which computation of G(h) is done
134      REAL(wp) ::   zhe        ! dummy factor
135      REAL(wp) ::   zcnd_i     ! mean sea ice thermal conductivity
136      !!------------------------------------------------------------------     
137
138      ! --- diag error on heat diffusion - PART 1 --- !
139      DO ji = 1, npti
140         zq_ini(ji) = ( SUM( e_i_1d(ji,1:nlay_i) ) * h_i_1d(ji) * r1_nlay_i +  &
141            &           SUM( e_s_1d(ji,1:nlay_s) ) * h_s_1d(ji) * r1_nlay_s ) 
142      END DO
143
144      !------------------
145      ! 1) Initialization
146      !------------------
147      DO ji = 1, npti
148         isnow(ji) = 1._wp - MAX( 0._wp , SIGN(1._wp, - h_s_1d(ji) ) )  ! is there snow or not
149         ! layer thickness
150         zh_i(ji) = h_i_1d(ji) * r1_nlay_i
151         zh_s(ji) = h_s_1d(ji) * r1_nlay_s
152      END DO
153      !
154      WHERE( zh_i(1:npti) >= epsi10 )   ;   z1_h_i(1:npti) = 1._wp / zh_i(1:npti)
155      ELSEWHERE                         ;   z1_h_i(1:npti) = 0._wp
156      END WHERE
157      !
158      WHERE( zh_s(1:npti) >= epsi10 )   ;   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) = qsr_ice_tr_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 * zh_s(ji) * 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) + qsr_ice_tr_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      ftr_ice_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)      = rcdic + zbeta * sz_i_1d(ji,1)      / MIN( -epsi10, t_i_1d(ji,1) - rt0 )
220               ztcond_i(ji,nlay_i) = rcdic + 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) = rcdic + 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)      = rcdic + 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) = rcdic + 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) = rcdic + 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_monocat )
255         !
256         CASE ( 1 , 3 )
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 = cpic + 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_rhoic * 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_rhosn * r1_cpic * 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) - qsr_ice_tr_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) + 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                  ztmelt_i      = -tmut * sz_i_1d(ji,jk) + rt0 
542                  t_i_1d(ji,jk) =  MAX( MIN( t_i_1d(ji,jk), ztmelt_i ), 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) + 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                  ztmelt_i      = -tmut * sz_i_1d(ji,jk) + rt0 
707                  t_i_1d(ji,jk) =  MAX( MIN( t_i_1d(ji,jk), ztmelt_i ), 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      ! --- update conduction fluxes
731      !
732      DO ji = 1, npti
733         !                                ! surface ice conduction flux
734         fc_su(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         fc_bo_i(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) - ( fc_su(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) - fc_bo_i(ji) + zdq * r1_rdtice )*a_i_1d(ji)
773               ELSE                          ! case T_su = 0degC
774                  zhfx_err = ( fc_su(ji)      + qsr_ice_tr_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq * r1_rdtice )*a_i_1d(ji)
775               ENDIF
776               
777            ELSEIF( k_jules == np_jules_ACTIVE ) THEN
778           
779               zhfx_err = ( fc_su(ji) + qsr_ice_tr_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq * r1_rdtice ) * a_i_1d(ji)
780           
781            ENDIF
782            !
783            ! total heat sink to be sent to the ocean
784            hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) + zhfx_err
785            !
786            ! hfx_dif = Heat flux diagnostic of sensible heat used to warm/cool ice in W.m-2   
787            hfx_dif_1d(ji) = hfx_dif_1d(ji) - zdq * r1_rdtice * a_i_1d(ji)
788            !
789         END DO
790         !
791         ! --- SIMIP diagnostics
792         !
793         DO ji = 1, npti
794            !--- Conduction fluxes (positive downwards)
795            diag_fc_bo_1d(ji) = diag_fc_bo_1d(ji) + fc_bo_i(ji) * a_i_1d(ji) / at_i_1d(ji)
796            diag_fc_su_1d(ji) = diag_fc_su_1d(ji) + fc_su  (ji) * a_i_1d(ji) / at_i_1d(ji)
797   
798            !--- Snow-ice interfacial temperature (diagnostic SIMIP)
799            zfac = rn_cnd_s * zh_i(ji) + ztcond_i(ji,1) * zh_s(ji)
800            IF( zh_s(ji) >= 1.e-3 .AND. zfac > epsi10 ) THEN
801               t_si_1d(ji) = ( rn_cnd_s       * zh_i(ji) * t_s_1d(ji,1) +   &
802                  &            ztcond_i(ji,1) * zh_s(ji) * t_i_1d(ji,1) ) / zfac
803            ELSE
804               t_si_1d(ji) = t_su_1d(ji)
805            ENDIF
806         END DO
807         !
808      ENDIF
809      !
810      !---------------------------------------------------------------------------------------
811      ! 11) Jules coupling: reset inner snow and ice temperatures, update conduction fluxes
812      !---------------------------------------------------------------------------------------
813      ! effective conductivity and 1st layer temperature (needed by Met Office)
814      DO ji = 1, npti
815         IF( h_s_1d(ji) > 0.1_wp ) THEN
816            cnd_ice_1d(ji) = 2._wp * zkappa_s(ji,0)
817         ELSE
818            IF( h_i_1d(ji) > 0.1_wp ) THEN
819               cnd_ice_1d(ji) = 2._wp * zkappa_i(ji,0)
820            ELSE
821               cnd_ice_1d(ji) = 2._wp * ztcond_i(ji,0) / 0.1_wp
822            ENDIF
823         ENDIF
824         t1_ice_1d(ji) = isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1)
825      END DO
826      !
827      IF( k_jules == np_jules_EMULE ) THEN
828         ! Restore temperatures to their initial values
829         t_s_1d    (1:npti,:) = ztsold(1:npti,:)
830         t_i_1d    (1:npti,:) = ztiold(1:npti,:)
831         qcn_ice_1d(1:npti)   = fc_su (1:npti)
832      ENDIF
833      !
834   END SUBROUTINE ice_thd_zdf_BL99
835
836
837#else
838   !!----------------------------------------------------------------------
839   !!   Default option       Dummy Module             No ESIM sea-ice model
840   !!----------------------------------------------------------------------
841#endif
842
843   !!======================================================================
844END MODULE icethd_zdf_BL99
Note: See TracBrowser for help on using the repository browser.