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/2020/dev_r13383_HPC-02_Daley_Tiling/src/ICE – NEMO

source: NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/ICE/icethd_zdf_bl99.F90 @ 13553

Last change on this file since 13553 was 13553, checked in by hadcv, 3 years ago

Merge in trunk up to [13550]

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