source: NEMO/branches/UKMO/NEMO_4.0.1_snow_frac_cpl/src/ICE/icethd_zdf_bl99.F90 @ 13046

Last change on this file since 13046 was 13046, checked in by dancopsey, 5 months ago

Add snow_frac which is fraction of snow for coupling purposes.

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