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

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

source: NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/ICE/icethd_zdf_bl99.F90 @ 9939

Last change on this file since 9939 was 9939, checked in by gm, 6 years ago

#1911 (ENHANCE-04): RK3 branche phased with MLF@9937 branche
File size: 42.8 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: icethd_zdf.F90 8420 2017-08-08 12:18:46Z clem $
34   !! Software governed by the CeCILL licence     (./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) ::   ztmelt_i                  ! 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) = 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 * 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) + 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)      = 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!!gm faster coding
225                  ztcond_i(ji,jk) = rcnd_i + zbeta * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) )             &
226                     &                          / MIN( t_i_1d (ji,jk) + t_i_1d (ji,jk+1) - 2._wp * rt0 , -epsi10 )
227!!gm old
228!                  ztcond_i(ji,jk) = rcnd_i + zbeta * 0.5_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) )   &
229!                     &                             / MIN( 0.5_wp * (t_i_1d(ji,jk) + t_i_1d(ji,jk+1)) - rt0 , -epsi10 )
230!!gm
231               END DO
232            END DO
233            !
234         ELSEIF( ln_cndi_P07 ) THEN     !-- Pringle et al formula: k = k0 + beta1.S/T - beta2.T
235            !
236            DO ji = 1, npti
237               ztcond_i(ji,0)      = rcnd_i + 0.09_wp  *  sz_i_1d(ji,1)      / MIN( -epsi10, t_i_1d(ji,1) - rt0 )  &
238                  &                         - 0.011_wp * ( t_i_1d(ji,1) - rt0 )
239               ztcond_i(ji,nlay_i) = rcnd_i + 0.09_wp  *  sz_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji)  - rt0 )  &
240                  &                         - 0.011_wp * ( t_bo_1d(ji) - rt0 )
241            END DO
242            DO jk = 1, nlay_i-1
243               DO ji = 1, npti
244!!gm faster coding
245                  zfac   = t_i_1d (ji,jk) + t_i_1d (ji,jk+1) - 2._wp * rt0
246                  ztcond_i(ji,jk) = rcnd_i + 0.09_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) ) / MIN( -epsi10, zfac )  &
247                     &                     - 0.0055_wp * zfac            !NB:  0.0055 = 1/2 * 0.011
248!!gm old
249!                  ztcond_i(ji,jk) = rcnd_i + 0.09_wp  *   0.5_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) ) /        &
250!                     &                      MIN( -epsi10, 0.5_wp * ( t_i_1d (ji,jk) + t_i_1d (ji,jk+1) ) - rt0 )  &
251!                     &                     - 0.011_wp * ( 0.5_wp * ( t_i_1d (ji,jk) + t_i_1d (ji,jk+1) ) - rt0 )
252!!gm
253               END DO
254            END DO
255            !
256         ENDIF
257         ztcond_i(1:npti,:) = MAX( zkimin, ztcond_i(1:npti,:) )       
258         !
259         !--- G(he) : enhancement of thermal conductivity in mono-category case
260         ! Computation of effective thermal conductivity G(h)
261         ! Used in mono-category case only to simulate an ITD implicitly
262         ! Fichefet and Morales Maqueda, JGR 1997
263         zghe(1:npti) = 1._wp
264         !
265         SELECT CASE ( nn_virtual_itd )
266         !
267         CASE ( 1 , 2 )
268            !
269            zepsilon = 0.1_wp
270            DO ji = 1, npti
271               zcnd_i = SUM( ztcond_i(ji,:) ) / REAL( nlay_i+1, wp )                                ! Mean sea ice thermal conductivity
272               zhe = ( rn_cnd_s * h_i_1d(ji) + zcnd_i * h_s_1d(ji) ) / ( rn_cnd_s + zcnd_i )        ! Effective thickness he (zhe)
273               IF( zhe >=  zepsilon * 0.5_wp * EXP(1._wp) )  &
274                  &   zghe(ji) = MIN( 2._wp, 0.5_wp * ( 1._wp + LOG( 2._wp * zhe / zepsilon ) ) )   ! G(he)
275            END DO
276            !
277         END SELECT
278         !
279         !-----------------
280         ! 4) kappa factors
281         !-----------------
282         !--- Snow
283         DO jk = 0, nlay_s-1
284            DO ji = 1, npti
285               zkappa_s(ji,jk) = zghe(ji) * rn_cnd_s * z1_h_s(ji)
286            END DO
287         END DO
288         DO ji = 1, npti   ! Snow-ice interface
289            zfac = 0.5_wp * ( ztcond_i(ji,0) * zh_s(ji) + rn_cnd_s * zh_i(ji) )
290            IF( zfac > epsi10 ) THEN
291               zkappa_s(ji,nlay_s) = zghe(ji) * rn_cnd_s * ztcond_i(ji,0) / zfac
292            ELSE
293               zkappa_s(ji,nlay_s) = 0._wp
294            ENDIF
295         END DO
296
297         !--- Ice
298         DO jk = 0, nlay_i
299            DO ji = 1, npti
300               zkappa_i(ji,jk) = zghe(ji) * ztcond_i(ji,jk) * z1_h_i(ji)
301            END DO
302         END DO
303         DO ji = 1, npti   ! Snow-ice interface
304            zkappa_i(ji,0) = zkappa_s(ji,nlay_s) * isnow(ji) + zkappa_i(ji,0) * ( 1._wp - isnow(ji) )
305         END DO
306         !
307         !--------------------------------------
308         ! 5) Sea ice specific heat, eta factors
309         !--------------------------------------
310         DO jk = 1, nlay_i
311            DO ji = 1, npti
312               zcpi = rcpi + zgamma * sz_i_1d(ji,jk) / MAX( ( t_i_1d(ji,jk) - rt0 ) * ( ztiold(ji,jk) - rt0 ), epsi10 )
313               zeta_i(ji,jk) = rdt_ice * r1_rhoi * z1_h_i(ji) / MAX( epsi10, zcpi ) 
314            END DO
315         END DO
316
317         DO jk = 1, nlay_s
318            DO ji = 1, npti
319               zeta_s(ji,jk) = rdt_ice * r1_rhos * r1_cpi * z1_h_s(ji)
320            END DO
321         END DO
322         !
323         !----------------------------------------!
324         !                                        !
325         !   JULES COUPLING IS OFF OR EMULATED    !
326         !                                        !
327         !----------------------------------------!
328         !
329         IF( k_jules == np_jules_OFF .OR. k_jules == np_jules_EMULE ) THEN
330            !
331            ! ==> The original BL99 temperature computation is used
332            !       (with qsr_ice, qns_ice and dqns_ice as inputs)
333            !
334            !----------------------------
335            ! 6) surface flux computation
336            !----------------------------
337            ! update of the non solar flux according to the update in T_su
338            DO ji = 1, npti
339               qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_1d(ji) - ztsub(ji) )
340            END DO
341
342            DO ji = 1, npti
343               zfnet(ji) = qsr_ice_1d(ji) - qsr_ice_tr_1d(ji) + qns_ice_1d(ji) ! net heat flux = net - transmitted solar + non solar
344            END DO
345            !
346            !----------------------------
347            ! 7) tridiagonal system terms
348            !----------------------------
349            ! layer denotes the number of the layer in the snow or in the ice
350            ! jm denotes the reference number of the equation in the tridiagonal
351            ! system, terms of tridiagonal system are indexed as following :
352            ! 1 is subdiagonal term, 2 is diagonal and 3 is superdiagonal one
353
354            ! ice interior terms (top equation has the same form as the others)
355            ztrid   (1:npti,:,:) = 0._wp
356            zindterm(1:npti,:)   = 0._wp
357            zindtbis(1:npti,:)   = 0._wp
358            zdiagbis(1:npti,:)   = 0._wp
359
360            DO jm = nlay_s + 2, nlay_s + nlay_i 
361               DO ji = 1, npti
362                  jk = jm - nlay_s - 1
363                  ztrid   (ji,jm,1) =       - zeta_i(ji,jk) *   zkappa_i(ji,jk-1)
364                  ztrid   (ji,jm,2) = 1._wp + zeta_i(ji,jk) * ( zkappa_i(ji,jk-1) + zkappa_i(ji,jk) )
365                  ztrid   (ji,jm,3) =       - zeta_i(ji,jk) *                       zkappa_i(ji,jk)
366                  zindterm(ji,jm)   = ztiold(ji,jk) + zeta_i(ji,jk) * zradab_i(ji,jk)
367               END DO
368            END DO
369
370            jm =  nlay_s + nlay_i + 1
371            DO ji = 1, npti
372               ! ice bottom term
373               ztrid   (ji,jm,1) =       - zeta_i(ji,nlay_i) *   zkappa_i(ji,nlay_i-1)   
374               ztrid   (ji,jm,2) = 1._wp + zeta_i(ji,nlay_i) * ( zkappa_i(ji,nlay_i-1) + zkappa_i(ji,nlay_i) * zg1 )
375               ztrid   (ji,jm,3) = 0._wp
376               zindterm(ji,jm)   = ztiold(ji,nlay_i) + zeta_i(ji,nlay_i) *  &
377                  &              ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i) * zg1 * t_bo_1d(ji) ) 
378            END DO
379
380            DO ji = 1, npti
381               !                               !---------------------!
382               IF( h_s_1d(ji) > 0._wp ) THEN   !  snow-covered cells !
383                  !                            !---------------------!
384                  ! snow interior terms (bottom equation has the same form as the others)
385                  DO jm = 3, nlay_s + 1
386                     jk = jm - 1
387                     ztrid   (ji,jm,1) =       - zeta_s(ji,jk) *   zkappa_s(ji,jk-1)
388                     ztrid   (ji,jm,2) = 1._wp + zeta_s(ji,jk) * ( zkappa_s(ji,jk-1) + zkappa_s(ji,jk) )
389                     ztrid   (ji,jm,3) =       - zeta_s(ji,jk) *                       zkappa_s(ji,jk)
390                     zindterm(ji,jm)   = ztsold(ji,jk) + zeta_s(ji,jk) * zradab_s(ji,jk)
391                  END DO
392                 
393                  ! case of only one layer in the ice (ice equation is altered)
394                  IF( nlay_i == 1 ) THEN
395                     ztrid   (ji,nlay_s+2,3) = 0._wp
396                     zindterm(ji,nlay_s+2)   = zindterm(ji,nlay_s+2) + zeta_i(ji,1) * zkappa_i(ji,1) * t_bo_1d(ji) 
397                  ENDIF
398                 
399                  IF( t_su_1d(ji) < rt0 ) THEN   !--  case 1 : no surface melting
400                     
401                     jm_min(ji) = 1
402                     jm_max(ji) = nlay_i + nlay_s + 1
403                     
404                     ! surface equation
405                     ztrid   (ji,1,1) = 0._wp
406                     ztrid   (ji,1,2) = zdqns_ice_b(ji) - zg1s * zkappa_s(ji,0)
407                     ztrid   (ji,1,3) =                   zg1s * zkappa_s(ji,0)
408                     zindterm(ji,1)   = zdqns_ice_b(ji) * t_su_1d(ji) - zfnet(ji)
409                     
410                     ! first layer of snow equation
411                     ztrid   (ji,2,1) =       - zeta_s(ji,1) *                    zkappa_s(ji,0) * zg1s
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)
415                     
416                  ELSE                            !--  case 2 : surface is melting
417                     !
418                     jm_min(ji) = 2
419                     jm_max(ji) = nlay_i + nlay_s + 1
420                     
421                     ! first layer of snow equation
422                     ztrid   (ji,2,1) = 0._wp
423                     ztrid   (ji,2,2) = 1._wp + zeta_s(ji,1) * ( zkappa_s(ji,1) + zkappa_s(ji,0) * zg1s )
424                     ztrid   (ji,2,3) =       - zeta_s(ji,1) *   zkappa_s(ji,1) 
425                     zindterm(ji,2)   = ztsold(ji,1) + zeta_s(ji,1) * ( zradab_s(ji,1) + zkappa_s(ji,0) * zg1s * t_su_1d(ji) ) 
426                  ENDIF
427                  !                            !---------------------!
428               ELSE                            ! cells without snow  !
429                  !                            !---------------------!
430                  !
431                  IF( t_su_1d(ji) < rt0 ) THEN   !--  case 1 : no surface melting
432                     !
433                     jm_min(ji) = nlay_s + 1
434                     jm_max(ji) = nlay_i + nlay_s + 1
435                     
436                     ! surface equation   
437                     ztrid   (ji,jm_min(ji),1) = 0._wp
438                     ztrid   (ji,jm_min(ji),2) = zdqns_ice_b(ji) - zkappa_i(ji,0) * zg1   
439                     ztrid   (ji,jm_min(ji),3) =                   zkappa_i(ji,0) * zg1
440                     zindterm(ji,jm_min(ji))   = zdqns_ice_b(ji) * t_su_1d(ji) - zfnet(ji)
441                     
442                     ! first layer of ice equation
443                     ztrid   (ji,jm_min(ji)+1,1) =       - zeta_i(ji,1) *                    zkappa_i(ji,0) * zg1
444                     ztrid   (ji,jm_min(ji)+1,2) = 1._wp + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 )
445                     ztrid   (ji,jm_min(ji)+1,3) =       - zeta_i(ji,1) *   zkappa_i(ji,1) 
446                     zindterm(ji,jm_min(ji)+1)   = ztiold(ji,1) + zeta_i(ji,1) * zradab_i(ji,1) 
447                     
448                     ! case of only one layer in the ice (surface & ice equations are altered)
449                     IF( nlay_i == 1 ) THEN
450                        ztrid   (ji,jm_min(ji),1)   = 0._wp
451                        ztrid   (ji,jm_min(ji),2)   = zdqns_ice_b(ji)      -   zkappa_i(ji,0) * 2._wp
452                        ztrid   (ji,jm_min(ji),3)   =                          zkappa_i(ji,0) * 2._wp
453                        ztrid   (ji,jm_min(ji)+1,1) =       - zeta_i(ji,1) *   zkappa_i(ji,0) * 2._wp
454                        ztrid   (ji,jm_min(ji)+1,2) = 1._wp + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2._wp + zkappa_i(ji,1) )
455                        ztrid   (ji,jm_min(ji)+1,3) = 0._wp
456                        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))
457                     ENDIF
458                     
459                  ELSE                            !--  case 2 : surface is melting
460                     
461                     jm_min(ji) = nlay_s + 2
462                     jm_max(ji) = nlay_i + nlay_s + 1
463                     
464                     ! first layer of ice equation
465                     ztrid   (ji,jm_min(ji),1) = 0._wp
466                     ztrid   (ji,jm_min(ji),2) = 1._wp + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 ) 
467                     ztrid   (ji,jm_min(ji),3) =       - zeta_i(ji,1) *   zkappa_i(ji,1)
468                     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)) 
469                     
470                     ! case of only one layer in the ice (surface & ice equations are altered)
471                     IF( nlay_i == 1 ) THEN
472                        ztrid   (ji,jm_min(ji),1) = 0._wp
473                        ztrid   (ji,jm_min(ji),2) = 1._wp + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2._wp + zkappa_i(ji,1) )
474                        ztrid   (ji,jm_min(ji),3) = 0._wp
475                        zindterm(ji,jm_min(ji))   = ztiold(ji,1) + zeta_i(ji,1) * ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) ) &
476                           &                      + t_su_1d(ji) * zeta_i(ji,1) * zkappa_i(ji,0) * 2._wp
477                     ENDIF
478                     
479                  ENDIF
480               ENDIF
481               !
482               zindtbis(ji,jm_min(ji)) = zindterm(ji,jm_min(ji))
483               zdiagbis(ji,jm_min(ji)) = ztrid   (ji,jm_min(ji),2)
484               !
485            END DO
486            !
487            !------------------------------
488            ! 8) tridiagonal system solving
489            !------------------------------
490            ! Solve the tridiagonal system with Gauss elimination method.
491            ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, McGraw-Hill 1984
492            jm_maxt = 0
493            jm_mint = nlay_i+5
494            DO ji = 1, npti
495               jm_mint = MIN(jm_min(ji),jm_mint)
496               jm_maxt = MAX(jm_max(ji),jm_maxt)
497            END DO
498
499            DO jk = jm_mint+1, jm_maxt
500               DO ji = 1, npti
501                  jm = MIN(MAX(jm_min(ji)+1,jk),jm_max(ji))
502                  zdiagbis(ji,jm) = ztrid   (ji,jm,2) - ztrid(ji,jm,1) * ztrid   (ji,jm-1,3) / zdiagbis(ji,jm-1)
503                  zindtbis(ji,jm) = zindterm(ji,jm  ) - ztrid(ji,jm,1) * zindtbis(ji,jm-1  ) / zdiagbis(ji,jm-1)
504               END DO
505            END DO
506
507            ! ice temperatures
508            DO ji = 1, npti
509               t_i_1d(ji,nlay_i) = zindtbis(ji,jm_max(ji)) / zdiagbis(ji,jm_max(ji))
510            END DO
511
512            DO jm = nlay_i + nlay_s, nlay_s + 2, -1
513               DO ji = 1, npti
514                  jk = jm - nlay_s - 1
515                  t_i_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,jm)
516               END DO
517            END DO
518
519            DO ji = 1, npti
520               ! snow temperatures     
521               IF( h_s_1d(ji) > 0._wp ) THEN
522                  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)
523               ENDIF
524               ! surface temperature
525               ztsub(ji) = t_su_1d(ji)
526               IF( t_su_1d(ji) < rt0 ) THEN
527                  t_su_1d(ji) = ( zindtbis(ji,jm_min(ji)) - ztrid(ji,jm_min(ji),3) *  &
528                     &          ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) ) ) / zdiagbis(ji,jm_min(ji))
529               ENDIF
530            END DO
531            !
532            !--------------------------------------------------------------
533            ! 9) Has the scheme converged?, end of the iterative procedure
534            !--------------------------------------------------------------
535            ! check that nowhere it has started to melt
536            ! zdti_max is a measure of error, it has to be under zdti_bnd
537            zdti_max = 0._wp
538            DO ji = 1, npti
539               t_su_1d(ji) = MAX( MIN( t_su_1d(ji) , rt0 ) , rt0 - 100._wp )
540               zdti_max    = MAX( zdti_max, ABS( t_su_1d(ji) - ztsub(ji) ) )     
541            END DO
542
543            DO jk = 1, nlay_s
544               DO ji = 1, npti
545                  t_s_1d(ji,jk) = MAX( MIN( t_s_1d(ji,jk), rt0 ), rt0 - 100._wp )
546                  zdti_max      = MAX( zdti_max, ABS( t_s_1d(ji,jk) - ztsb(ji,jk) ) )
547               END DO
548            END DO
549
550            DO jk = 1, nlay_i
551               DO ji = 1, npti
552                  ztmelt_i      = -tmut * sz_i_1d(ji,jk) + rt0 
553                  t_i_1d(ji,jk) =  MAX( MIN( t_i_1d(ji,jk), ztmelt_i ), rt0 - 100._wp )
554                  zdti_max      =  MAX( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) )
555               END DO
556            END DO
557
558            ! Compute spatial maximum over all errors
559            ! note that this could be optimized substantially by iterating only the non-converging points
560            IF( lk_mpp ) CALL mpp_max( zdti_max, kcom=ncomm_ice )
561         !
562         !----------------------------------------!
563         !                                        !
564         !      JULES COUPLING IS ACTIVE          !
565         !                                        !
566         !----------------------------------------!
567         !
568         ELSEIF( k_jules == np_jules_ACTIVE ) THEN
569            !
570            ! ==> we use a modified BL99 solver with conduction flux (qcn_ice) as forcing term
571            !
572            !----------------------------
573            ! 7) tridiagonal system terms
574            !----------------------------
575            ! layer denotes the number of the layer in the snow or in the ice
576            ! jm denotes the reference number of the equation in the tridiagonal
577            ! system, terms of tridiagonal system are indexed as following :
578            ! 1 is subdiagonal term, 2 is diagonal and 3 is superdiagonal one
579
580            ! ice interior terms (top equation has the same form as the others)
581            ztrid   (1:npti,:,:) = 0._wp
582            zindterm(1:npti,:)   = 0._wp
583            zindtbis(1:npti,:)   = 0._wp
584            zdiagbis(1:npti,:)   = 0._wp
585
586            DO jm = nlay_s + 2, nlay_s + nlay_i 
587               DO ji = 1, npti
588                  jk = jm - nlay_s - 1
589                  ztrid   (ji,jm,1) =       - zeta_i(ji,jk) *   zkappa_i(ji,jk-1)
590                  ztrid   (ji,jm,2) = 1._wp + zeta_i(ji,jk) * ( zkappa_i(ji,jk-1) + zkappa_i(ji,jk) )
591                  ztrid   (ji,jm,3) =       - zeta_i(ji,jk) *                       zkappa_i(ji,jk)
592                  zindterm(ji,jm)   = ztiold(ji,jk) + zeta_i(ji,jk) * zradab_i(ji,jk)
593               END DO
594            ENDDO
595
596            jm =  nlay_s + nlay_i + 1
597            DO ji = 1, npti
598               ! ice bottom term
599               ztrid   (ji,jm,1) =       - zeta_i(ji,nlay_i) *   zkappa_i(ji,nlay_i-1)   
600               ztrid   (ji,jm,2) = 1._wp + zeta_i(ji,nlay_i) * ( zkappa_i(ji,nlay_i-1) + zkappa_i(ji,nlay_i) * zg1 )
601               ztrid   (ji,jm,3) = 0._wp
602               zindterm(ji,jm)   = ztiold(ji,nlay_i) + zeta_i(ji,nlay_i) *  &
603                  &              ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i) * zg1 * t_bo_1d(ji) ) 
604            ENDDO
605
606            DO ji = 1, npti
607               !                               !---------------------!
608               IF( h_s_1d(ji) > 0._wp ) THEN   !  snow-covered cells !
609                  !                            !---------------------!
610                  ! snow interior terms (bottom equation has the same form as the others)
611                  DO jm = 3, nlay_s + 1
612                     jk = jm - 1
613                     ztrid   (ji,jm,1) =       - zeta_s(ji,jk) *   zkappa_s(ji,jk-1)
614                     ztrid   (ji,jm,2) = 1._wp + zeta_s(ji,jk) * ( zkappa_s(ji,jk-1) + zkappa_s(ji,jk) )
615                     ztrid   (ji,jm,3) =       - zeta_s(ji,jk) *                       zkappa_s(ji,jk)
616                     zindterm(ji,jm)   = ztsold(ji,jk) + zeta_s(ji,jk) * zradab_s(ji,jk)
617                  END DO
618                 
619                  ! case of only one layer in the ice (ice equation is altered)
620                  IF ( nlay_i == 1 ) THEN
621                     ztrid   (ji,nlay_s+2,3) = 0._wp
622                     zindterm(ji,nlay_s+2)   = zindterm(ji,nlay_s+2) + zeta_i(ji,1) * zkappa_i(ji,1) * t_bo_1d(ji) 
623                  ENDIF
624                 
625                  jm_min(ji) = 2
626                  jm_max(ji) = nlay_i + nlay_s + 1
627                 
628                  ! first layer of snow equation
629                  ztrid   (ji,2,1) = 0._wp
630                  ztrid   (ji,2,2) = 1._wp + zeta_s(ji,1) * zkappa_s(ji,1)
631                  ztrid   (ji,2,3) =       - zeta_s(ji,1) * zkappa_s(ji,1) 
632                  zindterm(ji,2)   = ztsold(ji,1) + zeta_s(ji,1) * ( zradab_s(ji,1) + qcn_ice_1d(ji) ) 
633                 
634                  !                            !---------------------!
635               ELSE                            ! cells without snow  !
636                  !                            !---------------------!
637                  jm_min(ji) = nlay_s + 2
638                  jm_max(ji) = nlay_i + nlay_s + 1
639                 
640                  ! first layer of ice equation
641                  ztrid   (ji,jm_min(ji),1) = 0._wp
642                  ztrid   (ji,jm_min(ji),2) = 1._wp + zeta_i(ji,1) * zkappa_i(ji,1) 
643                  ztrid   (ji,jm_min(ji),3) =       - zeta_i(ji,1) * zkappa_i(ji,1)
644                  zindterm(ji,jm_min(ji))   = ztiold(ji,1) + zeta_i(ji,1) * ( zradab_i(ji,1) + qcn_ice_1d(ji) )
645                 
646                  ! case of only one layer in the ice (surface & ice equations are altered)
647                  IF( nlay_i == 1 ) THEN
648                     ztrid   (ji,jm_min(ji),1) = 0._wp
649                     ztrid   (ji,jm_min(ji),2) = 1._wp + zeta_i(ji,1) * zkappa_i(ji,1)
650                     ztrid   (ji,jm_min(ji),3) = 0._wp
651                     zindterm(ji,jm_min(ji))   = ztiold(ji,1) + zeta_i(ji,1) *  &
652                        &                                     ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) + qcn_ice_1d(ji) )
653                  ENDIF
654                 
655               ENDIF
656               !
657               zindtbis(ji,jm_min(ji)) = zindterm(ji,jm_min(ji))
658               zdiagbis(ji,jm_min(ji)) = ztrid   (ji,jm_min(ji),2)
659               !
660            END DO
661            !
662            !------------------------------
663            ! 8) tridiagonal system solving
664            !------------------------------
665            ! Solve the tridiagonal system with Gauss elimination method.
666            ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, McGraw-Hill 1984
667            jm_maxt = 0
668            jm_mint = nlay_i+5
669            DO ji = 1, npti
670               jm_mint = MIN(jm_min(ji),jm_mint)
671               jm_maxt = MAX(jm_max(ji),jm_maxt)
672            END DO
673           
674            DO jk = jm_mint+1, jm_maxt
675               DO ji = 1, npti
676                  jm = MIN(MAX(jm_min(ji)+1,jk),jm_max(ji))
677                  zdiagbis(ji,jm) = ztrid   (ji,jm,2) - ztrid(ji,jm,1) * ztrid   (ji,jm-1,3) / zdiagbis(ji,jm-1)
678                  zindtbis(ji,jm) = zindterm(ji,jm)   - ztrid(ji,jm,1) * zindtbis(ji,jm-1)   / zdiagbis(ji,jm-1)
679               END DO
680            END DO
681           
682            ! ice temperatures
683           DO ji = 1, npti
684               t_i_1d(ji,nlay_i) = zindtbis(ji,jm_max(ji)) / zdiagbis(ji,jm_max(ji))
685            END DO
686
687            DO jm = nlay_i + nlay_s, nlay_s + 2, -1
688               DO ji = 1, npti
689                  jk = jm - nlay_s - 1
690                  t_i_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,jm)
691               END DO
692            END DO
693           
694            ! snow temperatures     
695            DO ji = 1, npti
696               IF( h_s_1d(ji) > 0._wp ) THEN
697                  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)
698               ENDIF
699            END DO
700            !
701            !--------------------------------------------------------------
702            ! 9) Has the scheme converged?, end of the iterative procedure
703            !--------------------------------------------------------------
704            ! check that nowhere it has started to melt
705            ! zdti_max is a measure of error, it has to be under zdti_bnd
706            zdti_max = 0._wp
707
708            DO jk = 1, nlay_s
709               DO ji = 1, npti
710                  t_s_1d(ji,jk) = MAX( MIN( t_s_1d(ji,jk), rt0 ), rt0 - 100._wp )
711                  zdti_max = MAX( zdti_max, ABS( t_s_1d(ji,jk) - ztsb(ji,jk) ) )
712               END DO
713            END DO
714           
715            DO jk = 1, nlay_i
716               DO ji = 1, npti
717                  ztmelt_i      = -tmut * sz_i_1d(ji,jk) + rt0 
718                  t_i_1d(ji,jk) =  MAX( MIN( t_i_1d(ji,jk), ztmelt_i ), rt0 - 100._wp )
719                  zdti_max      =  MAX( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) )
720               END DO
721            END DO
722
723            ! Compute spatial maximum over all errors
724            ! note that this could be optimized substantially by iterating only the non-converging points
725            IF( lk_mpp )   CALL mpp_max( zdti_max, kcom=ncomm_ice )
726
727         ENDIF ! k_jules
728         
729      END DO  ! End of the do while iterative procedure
730     
731      IF( ln_icectl .AND. lwp ) THEN
732         WRITE(numout,*) ' zdti_max : ', zdti_max
733         WRITE(numout,*) ' iconv    : ', iconv
734      ENDIF
735     
736      !
737      !-----------------------------
738      ! 10) Fluxes at the interfaces
739      !-----------------------------
740      !
741      ! --- update conduction fluxes
742      !
743      DO ji = 1, npti
744         !                                ! surface ice conduction flux
745         fc_su(ji)   =  -           isnow(ji)   * zkappa_s(ji,0) * zg1s * ( t_s_1d(ji,1) - t_su_1d(ji) )  &
746            &           - ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1  * ( t_i_1d(ji,1) - t_su_1d(ji) )
747         !                                ! bottom ice conduction flux
748         fc_bo_i(ji) =  - zkappa_i(ji,nlay_i) * zg1 * ( t_bo_1d(ji) - t_i_1d(ji,nlay_i) )
749      END DO
750     
751      !
752      ! --- Diagnose the heat loss due to changing non-solar / conduction flux --- !
753      !
754      IF( k_jules == np_jules_OFF .OR. k_jules == np_jules_EMULE ) THEN
755         !
756         DO ji = 1, npti
757            hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) - ( qns_ice_1d(ji) - zqns_ice_b(ji) ) * a_i_1d(ji) 
758         END DO
759         !
760      ELSEIF( k_jules == np_jules_ACTIVE ) THEN
761         !
762         DO ji = 1, npti
763            hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) - ( fc_su(ji)      - qcn_ice_1d(ji) ) * a_i_1d(ji) 
764         END DO
765         !
766      ENDIF
767     
768      !
769      ! --- Diagnose the heat loss due to non-fully converged temperature solution (should not be above 10-4 W-m2) --- !
770      !
771      IF( k_jules == np_jules_OFF .OR. k_jules == np_jules_ACTIVE ) THEN
772         
773         CALL ice_var_enthalpy       
774         
775         ! zhfx_err = correction on the diagnosed heat flux due to non-convergence of the algorithm used to solve heat equation
776         DO ji = 1, npti
777            zdq = - zq_ini(ji) + ( SUM( e_i_1d(ji,1:nlay_i) ) * h_i_1d(ji) * r1_nlay_i +  &
778               &                   SUM( e_s_1d(ji,1:nlay_s) ) * h_s_1d(ji) * r1_nlay_s )
779           
780            IF( k_jules == np_jules_OFF ) THEN
781               
782               IF( t_su_1d(ji) < rt0 ) THEN  ! case T_su < 0degC
783                  zhfx_err = ( qns_ice_1d(ji) + qsr_ice_1d(ji)    - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq * r1_Dt_ice )*a_i_1d(ji)
784               ELSE                          ! case T_su = 0degC
785                  zhfx_err = ( fc_su(ji)      + qsr_ice_tr_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq * r1_Dt_ice )*a_i_1d(ji)
786               ENDIF
787               
788            ELSEIF( k_jules == np_jules_ACTIVE ) THEN
789           
790               zhfx_err = ( fc_su(ji) + qsr_ice_tr_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq * r1_Dt_ice ) * a_i_1d(ji)
791           
792            ENDIF
793            !
794            ! total heat sink to be sent to the ocean
795            hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) + zhfx_err
796            !
797            ! hfx_dif = Heat flux diagnostic of sensible heat used to warm/cool ice in W.m-2   
798            hfx_dif_1d(ji) = hfx_dif_1d(ji) - zdq * r1_Dt_ice * a_i_1d(ji)
799            !
800         END DO
801         !
802         ! --- SIMIP diagnostics
803         !
804         DO ji = 1, npti
805            !--- Conduction fluxes (positive downwards)
806            diag_fc_bo_1d(ji) = diag_fc_bo_1d(ji) + fc_bo_i(ji) * a_i_1d(ji) / at_i_1d(ji)
807            diag_fc_su_1d(ji) = diag_fc_su_1d(ji) + fc_su  (ji) * a_i_1d(ji) / at_i_1d(ji)
808   
809            !--- Snow-ice interfacial temperature (diagnostic SIMIP)
810            zfac = rn_cnd_s * zh_i(ji) + ztcond_i(ji,1) * zh_s(ji)
811            IF( h_s_1d(ji) >= zhs_min ) THEN
812               t_si_1d(ji) = ( rn_cnd_s       * zh_i(ji) * t_s_1d(ji,1) +   &
813                  &            ztcond_i(ji,1) * zh_s(ji) * t_i_1d(ji,1) ) / MAX( epsi10, zfac )
814            ELSE
815               t_si_1d(ji) = t_su_1d(ji)
816            ENDIF
817         END DO
818         !
819      ENDIF
820      !
821      !---------------------------------------------------------------------------------------
822      ! 11) Jules coupling: reset inner snow and ice temperatures, update conduction fluxes
823      !---------------------------------------------------------------------------------------
824      ! effective conductivity and 1st layer temperature (needed by Met Office)
825      DO ji = 1, npti
826         IF( h_s_1d(ji) > 0.1_wp ) THEN
827            cnd_ice_1d(ji) = 2._wp * zkappa_s(ji,0)
828         ELSE
829            IF( h_i_1d(ji) > 0.1_wp ) THEN
830               cnd_ice_1d(ji) = 2._wp * zkappa_i(ji,0)
831            ELSE
832               cnd_ice_1d(ji) = 2._wp * ztcond_i(ji,0) / 0.1_wp
833            ENDIF
834         ENDIF
835         t1_ice_1d(ji) = isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1)
836      END DO
837      !
838      IF( k_jules == np_jules_EMULE ) THEN
839         ! Restore temperatures to their initial values
840         t_s_1d    (1:npti,:) = ztsold(1:npti,:)
841         t_i_1d    (1:npti,:) = ztiold(1:npti,:)
842         qcn_ice_1d(1:npti)   = fc_su (1:npti)
843      ENDIF
844      !
845   END SUBROUTINE ice_thd_zdf_BL99
846
847
848#else
849   !!----------------------------------------------------------------------
850   !!   Default option       Dummy Module             No SI3 sea-ice model
851   !!----------------------------------------------------------------------
852#endif
853
854   !!======================================================================
855END MODULE icethd_zdf_BL99
Note: See TracBrowser for help on using the repository browser.