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_r13648_ASINTER-04_laurent_bulk_ice/src/ICE – NEMO

source: NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/ICE/icethd_zdf_bl99.F90 @ 14021

Last change on this file since 14021 was 14021, checked in by laurent, 3 years ago

Caught up with trunk rev 14020...

  • Property svn:keywords set to Id
File size: 48.2 KB
Line 
1MODULE icethd_zdf_BL99
2   !!======================================================================
3   !!                       ***  MODULE icethd_zdf_BL99 ***
4   !!   sea-ice: vertical heat diffusion in sea ice (computation of temperatures)
5   !!======================================================================
6   !! History :       !  2003-02  (M. Vancoppenolle) original 1D code
7   !!                 !  2005-06  (M. Vancoppenolle) 3d version
8   !!            4.0  !  2018     (many people)      SI3 [aka Sea Ice cube]
9   !!----------------------------------------------------------------------
10#if defined key_si3
11   !!----------------------------------------------------------------------
12   !!   'key_si3'                                       SI3 sea-ice model
13   !!----------------------------------------------------------------------
14   !!  ice_thd_zdf_BL99 : vertical diffusion computation
15   !!----------------------------------------------------------------------
16   USE dom_oce        ! ocean space and time domain
17   USE phycst         ! physical constants (ocean directory)
18   USE ice            ! sea-ice: variables
19   USE ice1D          ! sea-ice: thermodynamics variables
20   USE icevar         ! sea-ice: operations
21   !
22   USE in_out_manager ! I/O manager
23   USE lib_mpp        ! MPP library
24   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero)
25
26   IMPLICIT NONE
27   PRIVATE
28
29   PUBLIC   ice_thd_zdf_BL99   ! called by icethd_zdf
30
31   !!----------------------------------------------------------------------
32   !! NEMO/ICE 4.0 , NEMO Consortium (2018)
33   !! $Id$
34   !! Software governed by the CeCILL license (see ./LICENSE)
35   !!----------------------------------------------------------------------
36CONTAINS
37
38   SUBROUTINE ice_thd_zdf_BL99( k_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)          ::   zq_ini      ! diag errors on heat
128      REAL(wp), DIMENSION(jpij)          ::   zghe        ! G(he), th. conduct enhancement factor, mono-cat
129      REAL(wp), DIMENSION(jpij)          ::   za_s_fra    ! ice fraction covered by snow
130      REAL(wp), DIMENSION(jpij)          ::   isnow       ! snow presence (1) or not (0)
131      REAL(wp), DIMENSION(jpij)          ::   isnow_comb  ! snow presence for met-office
132      REAL(wp), DIMENSION(jpij,nlay_i+nlay_s+1)   ::   zindterm    ! 'Ind'ependent term
133      REAL(wp), DIMENSION(jpij,nlay_i+nlay_s+1)   ::   zindtbis    ! Temporary 'ind'ependent term
134      REAL(wp), DIMENSION(jpij,nlay_i+nlay_s+1)   ::   zdiagbis    ! Temporary 'dia'gonal term
135      REAL(wp), DIMENSION(jpij,nlay_i+nlay_s+1,3) ::   ztrid       ! Tridiagonal system terms
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!!$            !!clem SNWLAY => check why LIM1D does not get this loop. Is nlay_i+5 correct?
542!!$
543!!$            DO jk = jm_mint+1, jm_maxt
544!!$               DO ji = 1, npti
545!!$                  jm = MIN(MAX(jm_min(ji)+1,jk),jm_max(ji))
546!!$                  zdiagbis(ji,jm) = ztrid   (ji,jm,2) - ztrid(ji,jm,1) * ztrid   (ji,jm-1,3) / zdiagbis(ji,jm-1)
547!!$                  zindtbis(ji,jm) = zindterm(ji,jm  ) - ztrid(ji,jm,1) * zindtbis(ji,jm-1  ) / zdiagbis(ji,jm-1)
548!!$               END DO
549!!$            END DO
550            ! clem: maybe one should find a way to reverse this loop for mpi performance
551            DO ji = 1, npti
552               jm_mint = jm_min(ji)
553               jm_maxt = jm_max(ji)
554               DO jm = jm_mint+1, jm_maxt
555                  zdiagbis(ji,jm) = ztrid   (ji,jm,2) - ztrid(ji,jm,1) * ztrid   (ji,jm-1,3) / zdiagbis(ji,jm-1)
556                  zindtbis(ji,jm) = zindterm(ji,jm  ) - ztrid(ji,jm,1) * zindtbis(ji,jm-1  ) / zdiagbis(ji,jm-1)
557               END DO
558            END DO
559
560            ! ice temperatures
561            DO ji = 1, npti
562               ! Variable used after iterations
563               ! Value must be frozen after convergence for MPP independance reason
564               IF ( .NOT. l_T_converged(ji) ) &
565                  t_i_1d(ji,nlay_i) = zindtbis(ji,jm_max(ji)) / zdiagbis(ji,jm_max(ji))
566            END DO
567
568            DO jm = nlay_i + nlay_s, nlay_s + 2, -1
569               DO ji = 1, npti
570                  jk = jm - nlay_s - 1
571                  IF ( .NOT. l_T_converged(ji) ) &
572                     t_i_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,jm)
573               END DO
574            END DO
575
576            ! snow temperatures
577            DO ji = 1, npti
578               ! Variables used after iterations
579               ! Value must be frozen after convergence for MPP independance reason
580               IF ( .NOT. l_T_converged(ji) .AND. h_s_1d(ji) > 0._wp ) &
581                  &   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)
582            END DO
583            !!clem SNWLAY
584            DO jm = nlay_s, 2, -1
585               DO ji = 1, npti
586                  jk = jm - 1
587                  IF ( .NOT. l_T_converged(ji) .AND. h_s_1d(ji) > 0._wp ) &
588                     &   t_s_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_s_1d(ji,jk+1) ) / zdiagbis(ji,jm)
589               END DO
590            END DO
591
592            ! surface temperature
593            DO ji = 1, npti
594               IF( .NOT. l_T_converged(ji) ) THEN
595                  ztsub(ji) = t_su_1d(ji)
596                  IF( t_su_1d(ji) < rt0 ) THEN
597                     t_su_1d(ji) = ( zindtbis(ji,jm_min(ji)) - ztrid(ji,jm_min(ji),3) *  &
598                        &          ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) ) ) / zdiagbis(ji,jm_min(ji))
599                  ENDIF
600               ENDIF
601            END DO
602            !
603            !--------------------------------------------------------------
604            ! 9) Has the scheme converged?, end of the iterative procedure
605            !--------------------------------------------------------------
606            ! check that nowhere it has started to melt
607            ! zdti_max is a measure of error, it has to be under zdti_bnd
608
609            DO ji = 1, npti
610
611               zdti_max = 0._wp
612
613               IF ( .NOT. l_T_converged(ji) ) THEN
614
615                  t_su_1d(ji) = MAX( MIN( t_su_1d(ji) , rt0 ) , rt0 - 100._wp )
616                  zdti_max    = MAX( zdti_max, ABS( t_su_1d(ji) - ztsub(ji) ) )
617
618                  IF( h_s_1d(ji) > 0._wp ) THEN
619                     DO jk = 1, nlay_s
620                        t_s_1d(ji,jk) = MAX( MIN( t_s_1d(ji,jk), rt0 ), rt0 - 100._wp )
621                        zdti_max      = MAX ( zdti_max , ABS( t_s_1d(ji,jk) - ztsb(ji,jk) ) )
622                     END DO
623                  ENDIF
624
625                  DO jk = 1, nlay_i
626                     ztmelts       = -rTmlt * sz_i_1d(ji,jk) + rt0
627                     t_i_1d(ji,jk) =  MAX( MIN( t_i_1d(ji,jk), ztmelts ), rt0 - 100._wp )
628                     zdti_max      =  MAX( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) )
629                  END DO
630
631                  ! convergence test
632                  IF( ln_zdf_chkcvg ) THEN
633                     tice_cvgerr_1d(ji) = zdti_max
634                     tice_cvgstp_1d(ji) = REAL(iconv)
635                  ENDIF
636
637                  IF( zdti_max < zdti_bnd )   l_T_converged(ji) = .TRUE.
638
639               ENDIF
640
641            END DO
642
643         !----------------------------------------!
644         !                                        !
645         !      Conduction flux is on             !
646         !                                        !
647         !----------------------------------------!
648         !
649         ELSEIF( k_cnd == np_cnd_ON ) THEN
650            !
651            ! ==> we use a modified BL99 solver with conduction flux (qcn_ice) as forcing term
652            !
653            !----------------------------
654            ! 7) tridiagonal system terms
655            !----------------------------
656            ! layer denotes the number of the layer in the snow or in the ice
657            ! jm denotes the reference number of the equation in the tridiagonal
658            ! system, terms of tridiagonal system are indexed as following :
659            ! 1 is subdiagonal term, 2 is diagonal and 3 is superdiagonal one
660
661            ! ice interior terms (top equation has the same form as the others)
662            ztrid   (1:npti,:,:) = 0._wp
663            zindterm(1:npti,:)   = 0._wp
664            zindtbis(1:npti,:)   = 0._wp
665            zdiagbis(1:npti,:)   = 0._wp
666
667            DO jm = nlay_s + 2, nlay_s + nlay_i
668               DO ji = 1, npti
669                  jk = jm - nlay_s - 1
670                  ztrid   (ji,jm,1) =       - zeta_i(ji,jk) *   zkappa_i(ji,jk-1)
671                  ztrid   (ji,jm,2) = 1._wp + zeta_i(ji,jk) * ( zkappa_i(ji,jk-1) + zkappa_i(ji,jk) )
672                  ztrid   (ji,jm,3) =       - zeta_i(ji,jk) *                       zkappa_i(ji,jk)
673                  zindterm(ji,jm)   = ztiold(ji,jk) + zeta_i(ji,jk) * zradab_i(ji,jk)
674               END DO
675            ENDDO
676
677            jm =  nlay_s + nlay_i + 1
678            DO ji = 1, npti
679               ! ice bottom term
680               ztrid   (ji,jm,1) =       - zeta_i(ji,nlay_i) *   zkappa_i(ji,nlay_i-1)
681               ztrid   (ji,jm,2) = 1._wp + zeta_i(ji,nlay_i) * ( zkappa_i(ji,nlay_i-1) + zkappa_i(ji,nlay_i) * zg1 )
682               ztrid   (ji,jm,3) = 0._wp
683               zindterm(ji,jm)   = ztiold(ji,nlay_i) + zeta_i(ji,nlay_i) *  &
684                  &              ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i) * zg1 * t_bo_1d(ji) )
685            ENDDO
686
687            DO ji = 1, npti
688               !                               !---------------------!
689               IF( h_s_1d(ji) > 0._wp ) THEN   !  snow-covered cells !
690                  !                            !---------------------!
691                  ! snow interior terms (bottom equation has the same form as the others)
692                  DO jm = 3, nlay_s + 1
693                     jk = jm - 1
694                     ztrid   (ji,jm,1) =       - zeta_s(ji,jk) *   zkappa_s(ji,jk-1)
695                     ztrid   (ji,jm,2) = 1._wp + zeta_s(ji,jk) * ( zkappa_s(ji,jk-1) + zkappa_s(ji,jk) )
696                     ztrid   (ji,jm,3) =       - zeta_s(ji,jk) *                       zkappa_s(ji,jk)
697                     zindterm(ji,jm)   = ztsold(ji,jk) + zeta_s(ji,jk) * zradab_s(ji,jk)
698                  END DO
699
700                  ! case of only one layer in the ice (ice equation is altered)
701                  IF ( nlay_i == 1 ) THEN
702                     ztrid   (ji,nlay_s+2,3) = 0._wp
703                     zindterm(ji,nlay_s+2)   = zindterm(ji,nlay_s+2) + zeta_i(ji,1) * zkappa_i(ji,1) * t_bo_1d(ji)
704                  ENDIF
705
706                  jm_min(ji) = 2
707                  jm_max(ji) = nlay_i + nlay_s + 1
708
709                  ! first layer of snow equation
710                  ztrid   (ji,2,1) = 0._wp
711                  ztrid   (ji,2,2) = 1._wp + zeta_s(ji,1) * zkappa_s(ji,1)
712                  ztrid   (ji,2,3) =       - zeta_s(ji,1) * zkappa_s(ji,1)
713                  zindterm(ji,2)   = ztsold(ji,1) + zeta_s(ji,1) * ( zradab_s(ji,1) + qcn_ice_1d(ji) )
714
715                  !                            !---------------------!
716               ELSE                            ! cells without snow  !
717                  !                            !---------------------!
718                  jm_min(ji) = nlay_s + 2
719                  jm_max(ji) = nlay_i + nlay_s + 1
720
721                  ! first layer of ice equation
722                  ztrid   (ji,jm_min(ji),1) = 0._wp
723                  ztrid   (ji,jm_min(ji),2) = 1._wp + zeta_i(ji,1) * zkappa_i(ji,1)
724                  ztrid   (ji,jm_min(ji),3) =       - zeta_i(ji,1) * zkappa_i(ji,1)
725                  zindterm(ji,jm_min(ji))   = ztiold(ji,1) + zeta_i(ji,1) * ( zradab_i(ji,1) + qcn_ice_1d(ji) )
726
727                  ! case of only one layer in the ice (surface & ice equations are altered)
728                  IF( nlay_i == 1 ) THEN
729                     ztrid   (ji,jm_min(ji),1) = 0._wp
730                     ztrid   (ji,jm_min(ji),2) = 1._wp + zeta_i(ji,1) * zkappa_i(ji,1)
731                     ztrid   (ji,jm_min(ji),3) = 0._wp
732                     zindterm(ji,jm_min(ji))   = ztiold(ji,1) + zeta_i(ji,1) *  &
733                        &                                     ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) + qcn_ice_1d(ji) )
734                  ENDIF
735
736               ENDIF
737               !
738               zindtbis(ji,jm_min(ji)) = zindterm(ji,jm_min(ji))
739               zdiagbis(ji,jm_min(ji)) = ztrid   (ji,jm_min(ji),2)
740               !
741            END DO
742            !
743            !------------------------------
744            ! 8) tridiagonal system solving
745            !------------------------------
746            ! Solve the tridiagonal system with Gauss elimination method.
747            ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, McGraw-Hill 1984
748!!$            jm_maxt = 0
749!!$            jm_mint = nlay_i+5
750!!$            DO ji = 1, npti
751!!$               jm_mint = MIN(jm_min(ji),jm_mint)
752!!$               jm_maxt = MAX(jm_max(ji),jm_maxt)
753!!$            END DO
754!!$
755!!$            DO jk = jm_mint+1, jm_maxt
756!!$               DO ji = 1, npti
757!!$                  jm = MIN(MAX(jm_min(ji)+1,jk),jm_max(ji))
758!!$                  zdiagbis(ji,jm) = ztrid   (ji,jm,2) - ztrid(ji,jm,1) * ztrid   (ji,jm-1,3) / zdiagbis(ji,jm-1)
759!!$                  zindtbis(ji,jm) = zindterm(ji,jm)   - ztrid(ji,jm,1) * zindtbis(ji,jm-1)   / zdiagbis(ji,jm-1)
760!!$               END DO
761!!$            END DO
762            ! clem: maybe one should find a way to reverse this loop for mpi performance
763            DO ji = 1, npti
764               jm_mint = jm_min(ji)
765               jm_maxt = jm_max(ji)
766               DO jm = jm_mint+1, jm_maxt
767                  zdiagbis(ji,jm) = ztrid   (ji,jm,2) - ztrid(ji,jm,1) * ztrid   (ji,jm-1,3) / zdiagbis(ji,jm-1)
768                  zindtbis(ji,jm) = zindterm(ji,jm  ) - ztrid(ji,jm,1) * zindtbis(ji,jm-1  ) / zdiagbis(ji,jm-1)
769               END DO
770            END DO
771
772            ! ice temperatures
773            DO ji = 1, npti
774               ! Variable used after iterations
775               ! Value must be frozen after convergence for MPP independance reason
776               IF ( .NOT. l_T_converged(ji) ) &
777                  t_i_1d(ji,nlay_i) = zindtbis(ji,jm_max(ji)) / zdiagbis(ji,jm_max(ji))
778            END DO
779
780            DO jm = nlay_i + nlay_s, nlay_s + 2, -1
781               DO ji = 1, npti
782                  IF ( .NOT. l_T_converged(ji) ) THEN
783                     jk = jm - nlay_s - 1
784                     t_i_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,jm)
785                  ENDIF
786               END DO
787            END DO
788
789            ! snow temperatures
790            DO ji = 1, npti
791               ! Variables used after iterations
792               ! Value must be frozen after convergence for MPP independance reason
793               IF ( .NOT. l_T_converged(ji) .AND. h_s_1d(ji) > 0._wp ) &
794                  &   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)
795            END DO
796            !!clem SNWLAY
797            DO jm = nlay_s, 2, -1
798               DO ji = 1, npti
799                  jk = jm - 1
800                  IF ( .NOT. l_T_converged(ji) .AND. h_s_1d(ji) > 0._wp ) &
801                     &   t_s_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_s_1d(ji,jk+1) ) / zdiagbis(ji,jm)
802               END DO
803            END DO
804            !
805            !--------------------------------------------------------------
806            ! 9) Has the scheme converged?, end of the iterative procedure
807            !--------------------------------------------------------------
808            ! check that nowhere it has started to melt
809            ! zdti_max is a measure of error, it has to be under zdti_bnd
810
811            DO ji = 1, npti
812
813               zdti_max = 0._wp
814
815               IF ( .NOT. l_T_converged(ji) ) THEN
816
817                  IF( h_s_1d(ji) > 0._wp ) THEN
818                     DO jk = 1, nlay_s
819                        t_s_1d(ji,jk) = MAX( MIN( t_s_1d(ji,jk), rt0 ), rt0 - 100._wp )
820                        zdti_max      = MAX ( zdti_max , ABS( t_s_1d(ji,jk) - ztsb(ji,jk) ) )
821                     END DO
822                  ENDIF
823
824                  DO jk = 1, nlay_i
825                     ztmelts       = -rTmlt * sz_i_1d(ji,jk) + rt0
826                     t_i_1d(ji,jk) =  MAX( MIN( t_i_1d(ji,jk), ztmelts ), rt0 - 100._wp )
827                     zdti_max      =  MAX ( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) )
828                  END DO
829
830                  ! convergence test
831                  IF( ln_zdf_chkcvg ) THEN
832                     tice_cvgerr_1d(ji) = zdti_max
833                     tice_cvgstp_1d(ji) = REAL(iconv)
834                  ENDIF
835
836                  IF( zdti_max < zdti_bnd )   l_T_converged(ji) = .TRUE.
837
838               ENDIF
839
840            END DO
841
842         ENDIF ! k_cnd
843
844      END DO  ! End of the do while iterative procedure
845      !
846      !-----------------------------
847      ! 10) Fluxes at the interfaces
848      !-----------------------------
849      !
850      ! --- calculate conduction fluxes (positive downward)
851      !     bottom ice conduction flux
852      DO ji = 1, npti
853         qcn_ice_bot_1d(ji) = - zkappa_i(ji,nlay_i) * zg1 * ( t_bo_1d(ji ) - t_i_1d (ji,nlay_i) )
854      END DO
855      !     surface ice conduction flux
856      IF( k_cnd == np_cnd_OFF .OR. k_cnd == np_cnd_EMU ) THEN
857         !
858         DO ji = 1, npti
859            qcn_ice_top_1d(ji) = -           isnow(ji)   * zkappa_s(ji,0) * zg1s * ( t_s_1d(ji,1) - t_su_1d(ji) ) &
860               &                 - ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1  * ( t_i_1d(ji,1) - t_su_1d(ji) )
861         END DO
862         !
863      ELSEIF( k_cnd == np_cnd_ON ) THEN
864         !
865         DO ji = 1, npti
866            qcn_ice_top_1d(ji) = qcn_ice_1d(ji)
867         END DO
868         !
869      ENDIF
870      !     surface ice temperature
871      IF( k_cnd == np_cnd_ON .AND. ln_cndemulate ) THEN
872         !
873         DO ji = 1, npti
874            t_su_1d(ji) = ( qcn_ice_top_1d(ji) +          isnow(ji)   * zkappa_s(ji,0) * zg1s * t_s_1d(ji,1) + &
875               &                                ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1  * t_i_1d(ji,1) ) &
876               &          / MAX( epsi10, isnow(ji) * zkappa_s(ji,0) * zg1s + ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1 )
877            t_su_1d(ji) = MAX( MIN( t_su_1d(ji), rt0 ), rt0 - 100._wp )  ! cap t_su
878         END DO
879         !
880      ENDIF
881      !
882      ! --- Diagnose the heat loss due to changing non-solar / conduction flux --- !
883      !
884      IF( k_cnd == np_cnd_OFF .OR. k_cnd == np_cnd_EMU ) THEN
885         !
886         DO ji = 1, npti
887            hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) - ( qns_ice_1d(ji) - zqns_ice_b(ji) ) * a_i_1d(ji)
888         END DO
889         !
890      ENDIF
891      !
892      ! --- Diagnose the heat loss due to non-fully converged temperature solution (should not be above 10-4 W-m2) --- !
893      !
894      IF( k_cnd == np_cnd_OFF .OR. k_cnd == np_cnd_ON ) THEN
895
896         CALL ice_var_enthalpy
897
898         ! zhfx_err = correction on the diagnosed heat flux due to non-convergence of the algorithm used to solve heat equation
899         DO ji = 1, npti
900            zdq = - zq_ini(ji) + ( SUM( e_i_1d(ji,1:nlay_i) ) * h_i_1d(ji) * r1_nlay_i +  &
901               &                   SUM( e_s_1d(ji,1:nlay_s) ) * h_s_1d(ji) * r1_nlay_s )
902
903            IF( k_cnd == np_cnd_OFF ) THEN
904
905               IF( t_su_1d(ji) < rt0 ) THEN  ! case T_su < 0degC
906                  zhfx_err = ( qns_ice_1d(ji)     + qsr_ice_1d(ji)     - zradtr_i(ji,nlay_i) - qcn_ice_bot_1d(ji)  &
907                     &       + zdq * r1_Dt_ice ) * a_i_1d(ji)
908               ELSE                          ! case T_su = 0degC
909                  zhfx_err = ( qcn_ice_top_1d(ji) + qtr_ice_top_1d(ji) - zradtr_i(ji,nlay_i) - qcn_ice_bot_1d(ji)  &
910                     &       + zdq * r1_Dt_ice ) * a_i_1d(ji)
911               ENDIF
912
913            ELSEIF( k_cnd == np_cnd_ON ) THEN
914
915               zhfx_err    = ( qcn_ice_top_1d(ji) + qtr_ice_top_1d(ji) - zradtr_i(ji,nlay_i) - qcn_ice_bot_1d(ji)  &
916                  &          + zdq * r1_Dt_ice ) * a_i_1d(ji)
917
918            ENDIF
919            !
920            ! total heat sink to be sent to the ocean
921            hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) + zhfx_err
922            !
923            ! hfx_dif = Heat flux diagnostic of sensible heat used to warm/cool ice in W.m-2
924            hfx_dif_1d(ji) = hfx_dif_1d(ji) - zdq * r1_Dt_ice * a_i_1d(ji)
925            !
926         END DO
927         !
928      ENDIF
929      !
930      !--------------------------------------------------------------------
931      ! 11) reset inner snow and ice temperatures, update conduction fluxes
932      !--------------------------------------------------------------------
933      ! effective conductivity and 1st layer temperature (needed by Met Office)
934      ! this is a conductivity at mid-layer, hence the factor 2
935      DO ji = 1, npti
936         IF( h_i_1d(ji) >= zhi_ssl ) THEN
937            cnd_ice_1d(ji) = 2._wp * zkappa_comb(ji)
938            !!cnd_ice_1d(ji) = 2._wp * zkappa_i(ji,0)
939         ELSE
940            cnd_ice_1d(ji) = 2._wp * ztcond_i(ji,0) / zhi_ssl ! cnd_ice is capped by: cond_i/zhi_ssl
941         ENDIF
942         t1_ice_1d(ji) = isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1)
943      END DO
944      !
945      IF( k_cnd == np_cnd_EMU ) THEN
946         ! Restore temperatures to their initial values
947         t_s_1d    (1:npti,:) = ztsold        (1:npti,:)
948         t_i_1d    (1:npti,:) = ztiold        (1:npti,:)
949         qcn_ice_1d(1:npti)   = qcn_ice_top_1d(1:npti)
950      ENDIF
951      !
952      ! --- SIMIP diagnostics
953      !
954      DO ji = 1, npti
955         !--- Snow-ice interfacial temperature (diagnostic SIMIP)
956         IF( h_s_1d(ji) >= zhs_ssl ) THEN
957            t_si_1d(ji) = (   rn_cnd_s       * h_i_1d(ji) * r1_nlay_i * t_s_1d(ji,nlay_s)   &
958               &            + ztcond_i(ji,1) * h_s_1d(ji) * r1_nlay_s * t_i_1d(ji,1)      ) &
959               &          / ( rn_cnd_s       * h_i_1d(ji) * r1_nlay_i &
960               &            + ztcond_i(ji,1) * h_s_1d(ji) * r1_nlay_s )
961         ELSE
962            t_si_1d(ji) = t_su_1d(ji)
963         ENDIF
964      END DO
965      !
966   END SUBROUTINE ice_thd_zdf_BL99
967
968#else
969   !!----------------------------------------------------------------------
970   !!   Default option       Dummy Module             No SI3 sea-ice model
971   !!----------------------------------------------------------------------
972#endif
973
974   !!======================================================================
975END MODULE icethd_zdf_BL99
Note: See TracBrowser for help on using the repository browser.