source: NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icethd_zdf_bl99.F90 @ 12475

Last change on this file since 12475 was 12475, checked in by dancopsey, 14 months ago
  • Add more print statements.
  • Move away from using snow to ice diagnostics and use a new snow to pond one instead.
File size: 53.9 KB
Line 
1MODULE icethd_zdf_BL99
2   !!======================================================================
3   !!                       ***  MODULE icethd_zdf_BL99 ***
4   !!   sea-ice: vertical heat diffusion in sea ice (computation of temperatures)
5   !!======================================================================
6   !! History :       !  2003-02  (M. Vancoppenolle) original 1D code
7   !!                 !  2005-06  (M. Vancoppenolle) 3d version
8   !!            4.0  !  2018     (many people)      SI3 [aka Sea Ice cube]
9   !!----------------------------------------------------------------------
10#if defined key_si3
11   !!----------------------------------------------------------------------
12   !!   'key_si3'                                       SI3 sea-ice model
13   !!----------------------------------------------------------------------
14   !!  ice_thd_zdf_BL99 : vertical diffusion computation
15   !!----------------------------------------------------------------------
16   USE dom_oce        ! ocean space and time domain
17   USE phycst         ! physical constants (ocean directory)
18   USE ice            ! sea-ice: variables
19   USE ice1D          ! sea-ice: thermodynamics variables
20   USE icevar         ! sea-ice: operations
21   !
22   USE in_out_manager ! I/O manager
23   USE lib_mpp        ! MPP library
24   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero)
25
26   IMPLICIT NONE
27   PRIVATE
28
29   PUBLIC   ice_thd_zdf_BL99   ! called by icethd_zdf
30
31   !!----------------------------------------------------------------------
32   !! NEMO/ICE 4.0 , NEMO Consortium (2018)
33   !! $Id: icethd_zdf_bl99.F90 10926 2019-05-03 12:32:10Z clem $
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.13_wp     ! for thermal conductivity (could be 0.13)
92      REAL(wp) ::   zraext_s  =  10._wp       ! extinction coefficient of radiation in the snow
93      REAL(wp) ::   zkimin    =  0.10_wp      ! minimum ice thermal conductivity
94      REAL(wp) ::   ztsu_err  =  1.e-5_wp     ! range around which t_su is considered at 0C
95      REAL(wp) ::   zdti_bnd  =  1.e-4_wp     ! maximal authorized error on temperature
96      REAL(wp) ::   zhs_min   =  0.1_wp       ! minimum snow thickness for conductivity calculation
97      REAL(wp) ::   ztmelts                   ! ice melting temperature
98      REAL(wp) ::   zdti_max                  ! current maximal error on temperature
99      REAL(wp) ::   zcpi                      ! Ice specific heat
100      REAL(wp) ::   zhfx_err, zdq             ! diag errors on heat
101      REAL(wp) ::   zfac                      ! dummy factor
102      !
103      REAL(wp), DIMENSION(jpij) ::   isnow        ! fraction of sea ice that is snow covered (for thermodynamic use only)
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,nlay_i+3)   ::   zindterm    ! 'Ind'ependent term
127      REAL(wp), DIMENSION(jpij,nlay_i+3)   ::   zindtbis    ! Temporary 'ind'ependent term
128      REAL(wp), DIMENSION(jpij,nlay_i+3)   ::   zdiagbis    ! Temporary 'dia'gonal term
129      REAL(wp), DIMENSION(jpij,nlay_i+3,3) ::   ztrid       ! Tridiagonal system terms
130      REAL(wp), DIMENSION(jpij)            ::   zq_ini      ! diag errors on heat
131      REAL(wp), DIMENSION(jpij)            ::   zghe        ! G(he), th. conduct enhancement factor, mono-cat
132      !
133      ! Mono-category
134      REAL(wp) ::   zepsilon   ! determines thres. above which computation of G(h) is done
135      REAL(wp) ::   zhe        ! dummy factor
136      REAL(wp) ::   zcnd_i     ! mean sea ice thermal conductivity
137      !!------------------------------------------------------------------     
138
139      ! --- diag error on heat diffusion - PART 1 --- !
140      DO ji = 1, npti
141         zq_ini(ji) = ( SUM( e_i_1d(ji,1:nlay_i) ) * h_i_1d(ji) * r1_nlay_i +  &
142            &           SUM( e_s_1d(ji,1:nlay_s) ) * h_s_1d(ji) * r1_nlay_s ) 
143      END DO
144
145         DO ji = 1, npti
146           IF (to_print(ji) == 10) THEN
147             write(numout,*) 'icethd_zdf_bl99 0: t_i_1d(ji,1), e_i_1d(ji,1) = ',t_i_1d(ji,1), ' ', e_i_1d(ji,1)
148           END IF
149         END DO
150
151      !------------------
152      ! 1) Initialization
153      !------------------
154      DO ji = 1, npti
155         
156         ! If the snow thickness drops below zhs_min then reduce the snow fraction instead
157         IF( h_s_1d(ji) < zhs_min ) THEN
158           isnow(ji) = h_s_1d(ji) / zhs_min 
159           zh_s(ji) = zhs_min * r1_nlay_s 
160         ELSE
161           isnow(ji) = 1.0_wp 
162           zh_s(ji) = h_s_1d(ji) * r1_nlay_s 
163         END IF 
164
165         ! layer thickness
166         zh_i(ji) = h_i_1d(ji) * r1_nlay_i
167
168         IF( to_print(ji) == 10 ) THEN
169           write(numout,*)'icethd_zdf_bl99: v_i_1d(ji), a_i_1d(ji), h_i_1d(ji) = ',v_i_1d(ji), ' ', a_i_1d(ji), ' ', h_i_1d(ji)
170         END IF
171      END DO
172      !
173      WHERE( zh_i(1:npti) >= epsi10 )   ;   z1_h_i(1:npti) = 1._wp / zh_i(1:npti)
174      ELSEWHERE                         ;   z1_h_i(1:npti) = 0._wp
175      END WHERE
176      !
177      !
178      WHERE( zh_s(1:npti) > 0._wp   )   ;   z1_h_s(1:npti) = 1._wp / zh_s(1:npti)
179      ELSEWHERE                         ;   z1_h_s(1:npti) = 0._wp
180      END WHERE
181      !
182      ! Store initial temperatures and non solar heat fluxes
183      IF( k_cnd == np_cnd_OFF .OR. k_cnd == np_cnd_EMU ) THEN
184         !
185         ztsub      (1:npti) = t_su_1d(1:npti)                          ! surface temperature at iteration n-1
186         ztsuold    (1:npti) = t_su_1d(1:npti)                          ! surface temperature initial value
187         t_su_1d    (1:npti) = MIN( t_su_1d(1:npti), rt0 - ztsu_err )   ! required to leave the choice between melting or not
188         zdqns_ice_b(1:npti) = dqns_ice_1d(1:npti)                      ! derivative of incoming nonsolar flux
189         zqns_ice_b (1:npti) = qns_ice_1d(1:npti)                       ! store previous qns_ice_1d value
190         !
191      ENDIF
192      !
193      ztsold (1:npti,:) = t_s_1d(1:npti,:)   ! Old snow temperature
194      ztiold (1:npti,:) = t_i_1d(1:npti,:)   ! Old ice temperature
195
196      !-------------
197      ! 2) Radiation
198      !-------------
199      ! --- Transmission/absorption of solar radiation in the ice --- !
200      zradtr_s(1:npti,0) = qtr_ice_top_1d(1:npti)
201      DO jk = 1, nlay_s
202         DO ji = 1, npti
203            !                             ! radiation transmitted below the layer-th snow layer
204            zradtr_s(ji,jk) = zradtr_s(ji,0) * EXP( - zraext_s * h_s_1d(ji) * r1_nlay_s * REAL(jk) )
205            !                             ! radiation absorbed by the layer-th snow layer
206            zradab_s(ji,jk) = zradtr_s(ji,jk-1) - zradtr_s(ji,jk)
207         END DO
208      END DO
209      !
210      zradtr_i(1:npti,0) = zradtr_s(1:npti,nlay_s) * isnow(1:npti) + qtr_ice_top_1d(1:npti) * ( 1._wp - isnow(1:npti) )
211      DO jk = 1, nlay_i 
212         DO ji = 1, npti
213            !                             ! radiation transmitted below the layer-th ice layer
214            zradtr_i(ji,jk) = zradtr_i(ji,0) * EXP( - rn_kappa_i * zh_i(ji) * REAL(jk) )
215            !                             ! radiation absorbed by the layer-th ice layer
216            zradab_i(ji,jk) = zradtr_i(ji,jk-1) - zradtr_i(ji,jk)
217         END DO
218      END DO
219      !
220      qtr_ice_bot_1d(1:npti) = zradtr_i(1:npti,nlay_i)   ! record radiation transmitted below the ice
221      !
222      iconv    = 0          ! number of iterations
223      !
224      l_T_converged(:) = .FALSE.
225      ! Convergence calculated until all sub-domain grid points have converged
226      ! Calculations keep going for all grid points until sub-domain convergence (vectorisation optimisation)
227      ! but values are not taken into account (results independant of MPI partitioning)
228      !
229      !                                                                            !============================!
230      DO WHILE ( ( .NOT. ALL (l_T_converged(1:npti)) ) .AND. iconv < iconv_max )   ! Iterative procedure begins !
231         !                                                                         !============================!
232         iconv = iconv + 1
233         !
234         ztib(1:npti,:) = t_i_1d(1:npti,:)
235         ztsb(1:npti,:) = t_s_1d(1:npti,:)
236         !
237         !--------------------------------
238         ! 3) Sea ice thermal conductivity
239         !--------------------------------
240         IF( ln_cndi_U64 ) THEN         !-- Untersteiner (1964) formula: k = k0 + beta.S/T
241            !
242            DO ji = 1, npti
243               ztcond_i_cp(ji,0)      = rcnd_i + zbeta * sz_i_1d(ji,1)      / MIN( -epsi10, t_i_1d(ji,1) - rt0 )
244               ztcond_i_cp(ji,nlay_i) = rcnd_i + zbeta * sz_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji)  - rt0 )
245            END DO
246            DO jk = 1, nlay_i-1
247               DO ji = 1, npti
248                  ztcond_i_cp(ji,jk) = rcnd_i + zbeta * 0.5_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) ) /  &
249                     &                         MIN( -epsi10, 0.5_wp * (t_i_1d(ji,jk) + t_i_1d(ji,jk+1)) - rt0 )
250               END DO
251            END DO
252            !
253         ELSEIF( ln_cndi_P07 ) THEN     !-- Pringle et al formula: k = k0 + beta1.S/T - beta2.T
254            !
255            DO ji = 1, npti
256               ztcond_i_cp(ji,0)      = rcnd_i + 0.09_wp  *  sz_i_1d(ji,1)      / MIN( -epsi10, t_i_1d(ji,1) - rt0 )  &
257                  &                           - 0.011_wp * ( t_i_1d(ji,1) - rt0 )
258               ztcond_i_cp(ji,nlay_i) = rcnd_i + 0.09_wp  *  sz_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji)  - rt0 )  &
259                  &                           - 0.011_wp * ( 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 + 0.09_wp  *   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                     &                       - 0.011_wp * ( 0.5_wp * ( t_i_1d (ji,jk) + t_i_1d (ji,jk+1) ) - rt0 )
266               END DO
267            END DO
268            !
269         ENDIF
270
271         ! Variable used after iterations
272         ! Value must be frozen after convergence for MPP independance reason
273         DO ji = 1, npti
274            IF ( .NOT. l_T_converged(ji) ) &
275               ztcond_i(ji,:) = MAX( zkimin, ztcond_i_cp(ji,:) )   
276
277            IF (to_print(ji) == 10) THEN
278               write(numout,*)'icethd_zdf_bl99: ln_cndi_U64, ln_cndi_P07, rcnd_i = ',ln_cndi_U64, ln_cndi_P07, rcnd_i
279               write(numout,*)'icethd_zdf_bl99: sz_i_1d(ji,1), epsi10, t_i_1d(ji,1) = ',sz_i_1d(ji,1), epsi10, t_i_1d(ji,1)
280               write(numout,*)'icethd_zdf_bl99: rt0, nlay_i, sz_i_1d(ji,nlay_i), t_bo_1d(ji) = ', rt0, nlay_i, sz_i_1d(ji,nlay_i), t_bo_1d(ji)
281               write(numout,*)'icethd_zdf_bl99: sz_i_1d(ji,:) = ',sz_i_1d(ji,:)
282               write(numout,*)'icethd_zdf_bl99: t_i_1d(ji,:) = ',t_i_1d(ji,:)
283               write(numout,*)'icethd_zdf_bl99: ztcond_i_cp(ji,:) = ',ztcond_i_cp(ji,:)
284               write(numout,*)'icethd_zdf_bl99: zkimin = ',zkimin
285               write(numout,*)'icethd_zdf_bl99: ztcond_i = ',ztcond_i(ji,:)
286            ENDIF
287         END DO
288
289         
290
291         !
292         !--- G(he) : enhancement of thermal conductivity in mono-category case
293         ! Computation of effective thermal conductivity G(h)
294         ! Used in mono-category case only to simulate an ITD implicitly
295         ! Fichefet and Morales Maqueda, JGR 1997
296         zghe(1:npti) = 1._wp
297         !
298         IF( ln_virtual_itd ) THEN
299            !
300            zepsilon = 0.1_wp
301            DO ji = 1, npti
302               zcnd_i = SUM( ztcond_i(ji,:) ) / REAL( nlay_i+1, wp )                                ! Mean sea ice thermal conductivity
303               zhe = ( rn_cnd_s * h_i_1d(ji) + zcnd_i * h_s_1d(ji) ) / ( rn_cnd_s + zcnd_i )        ! Effective thickness he (zhe)
304               IF( zhe >=  zepsilon * 0.5_wp * EXP(1._wp) )  &
305                  &   zghe(ji) = MIN( 2._wp, 0.5_wp * ( 1._wp + LOG( 2._wp * zhe / zepsilon ) ) )   ! G(he)
306            END DO
307            !
308         ENDIF
309         !
310         !-----------------
311         ! 4) kappa factors
312         !-----------------
313         !--- Snow
314         ! Variable used after iterations
315         ! Value must be frozen after convergence for MPP independance reason
316         DO jk = 0, nlay_s-1
317            DO ji = 1, npti
318               IF ( .NOT. l_T_converged(ji) ) &
319                  zkappa_s(ji,jk) = zghe(ji) * rn_cnd_s * z1_h_s(ji)
320            END DO
321         END DO
322         DO ji = 1, npti   ! Snow-ice interface
323            IF ( .NOT. l_T_converged(ji) ) THEN
324               zfac = 0.5_wp * ( ztcond_i(ji,0) * zh_s(ji) + rn_cnd_s * zh_i(ji) )
325               IF( zfac > epsi10 ) THEN
326                  zkappa_s(ji,nlay_s) = zghe(ji) * rn_cnd_s * ztcond_i(ji,0) / zfac
327               ELSE
328                  zkappa_s(ji,nlay_s) = 0._wp
329               ENDIF
330            ENDIF
331         END DO
332
333         !--- Ice
334         ! Variable used after iterations
335         ! Value must be frozen after convergence for MPP independance reason
336         DO jk = 0, nlay_i
337            DO ji = 1, npti
338               IF ( .NOT. l_T_converged(ji) ) &
339                  zkappa_i(ji,jk) = zghe(ji) * ztcond_i(ji,jk) * z1_h_i(ji)
340            END DO
341         END DO
342         DO ji = 1, npti   ! Snow-ice interface
343           
344            IF (to_print(ji) == 10) THEN
345               write(numout,*) 'icethd_zdf_bl99: zkappa_i(ji,0), ztcond_i(ji,0), z1_h_i(ji), isnow(ji) = ',zkappa_i(ji,0), ' ', ztcond_i(ji,0), ' ', z1_h_i(ji), ' ', isnow(ji)
346            END IF
347         END DO
348         !
349         !--------------------------------------
350         ! 5) Sea ice specific heat, eta factors
351         !--------------------------------------
352         DO jk = 1, nlay_i
353            DO ji = 1, npti
354               zcpi = rcpi + zgamma * sz_i_1d(ji,jk) / MAX( ( t_i_1d(ji,jk) - rt0 ) * ( ztiold(ji,jk) - rt0 ), epsi10 )
355               zeta_i(ji,jk) = rdt_ice * r1_rhoi * z1_h_i(ji) / MAX( epsi10, zcpi ) 
356
357               IF (to_print(ji) == 10 .AND. jk == 1 .AND. cat == 1) THEN
358                  write(numout,*) 'zeta_i(ji,jk), z1_h_i(ji), zcpi = ',zeta_i(ji,jk), ' ', z1_h_i(ji), ' ', zcpi
359                  write(numout,*) 'sz_i_1d(ji,jk), t_i_1d(ji,jk), ztiold(ji,jk) = ',sz_i_1d(ji,jk), ' ', t_i_1d(ji,jk), ' ', ztiold(ji,jk)
360                  write(numout,*) 'iconv, l_T_converged(ji) = ',iconv, ' ',l_T_converged(ji)
361               END IF
362            END DO
363         END DO
364
365         DO jk = 1, nlay_s
366            DO ji = 1, npti
367               zeta_s(ji,jk) = rdt_ice * r1_rhos * r1_rcpi * z1_h_s(ji)
368            END DO
369         END DO
370
371         !
372         !----------------------------------------!
373         !                                        !
374         !   Conduction flux is off or emulated   !
375         !                                        !
376         !----------------------------------------!
377         !
378         IF( k_cnd == np_cnd_OFF .OR. k_cnd == np_cnd_EMU ) THEN
379            !
380            ! ==> The original BL99 temperature computation is used
381            !       (with qsr_ice, qns_ice and dqns_ice as inputs)
382            !
383            !----------------------------
384            ! 6) surface flux computation
385            !----------------------------
386            ! update of the non solar flux according to the update in T_su
387            DO ji = 1, npti
388               ! Variable used after iterations
389               ! Value must be frozen after convergence for MPP independance reason
390               IF ( .NOT. l_T_converged(ji) ) &
391                  qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_1d(ji) - ztsub(ji) )
392            END DO
393
394            DO ji = 1, npti
395               zfnet(ji) = qsr_ice_1d(ji) - qtr_ice_top_1d(ji) + qns_ice_1d(ji) ! net heat flux = net - transmitted solar + non solar
396            END DO
397
398         DO ji = 1, npti
399           IF (to_print(ji) == 10) THEN
400             write(numout,*) 'icethd_zdf_bl99 16 iter: t_i_1d(ji,1), e_i_1d(ji,1) = ',t_i_1d(ji,1), ' ', e_i_1d(ji,1)
401           END IF
402         END DO
403            !
404            !----------------------------
405            ! 7) tridiagonal system terms
406            !----------------------------
407            ! layer denotes the number of the layer in the snow or in the ice
408            ! jm denotes the reference number of the equation in the tridiagonal
409            ! system, terms of tridiagonal system are indexed as following :
410            ! 1 is subdiagonal term, 2 is diagonal and 3 is superdiagonal one
411
412            ! ice interior terms (top equation has the same form as the others)
413            ztrid   (1:npti,:,:) = 0._wp
414            zindterm(1:npti,:)   = 0._wp
415            zindtbis(1:npti,:)   = 0._wp
416            zdiagbis(1:npti,:)   = 0._wp
417
418            DO jm = nlay_s + 2, nlay_s + nlay_i 
419               DO ji = 1, npti
420                  jk = jm - nlay_s - 1
421                  ztrid   (ji,jm,1) =       - zeta_i(ji,jk) *   zkappa_i(ji,jk-1)
422                  ztrid   (ji,jm,2) = 1._wp + zeta_i(ji,jk) * ( zkappa_i(ji,jk-1) + zkappa_i(ji,jk) )
423                  ztrid   (ji,jm,3) =       - zeta_i(ji,jk) *                       zkappa_i(ji,jk)
424                  zindterm(ji,jm)   = ztiold(ji,jk) + zeta_i(ji,jk) * zradab_i(ji,jk)
425               END DO
426            END DO
427
428            jm =  nlay_s + nlay_i + 1
429            DO ji = 1, npti
430               ! ice bottom term
431               ztrid   (ji,jm,1) =       - zeta_i(ji,nlay_i) *   zkappa_i(ji,nlay_i-1)   
432               ztrid   (ji,jm,2) = 1._wp + zeta_i(ji,nlay_i) * ( zkappa_i(ji,nlay_i-1) + zkappa_i(ji,nlay_i) * zg1 )
433               ztrid   (ji,jm,3) = 0._wp
434               zindterm(ji,jm)   = ztiold(ji,nlay_i) + zeta_i(ji,nlay_i) *  &
435                  &              ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i) * zg1 * t_bo_1d(ji) ) 
436            END DO
437
438         DO ji = 1, npti
439           IF (to_print(ji) == 10) THEN
440             write(numout,*) 'icethd_zdf_bl99 17.5 iter: t_i_1d(ji,1), e_i_1d(ji,1) = ',t_i_1d(ji,1), ' ', e_i_1d(ji,1)
441           END IF
442         END DO
443
444            DO ji = 1, npti
445               !                               !---------------------!
446               IF( h_s_1d(ji) > 0._wp ) THEN   !  snow-covered cells !
447                  !                            !---------------------!
448                  ! snow interior terms (bottom equation has the same form as the others)
449                  DO jm = 3, nlay_s + 1
450                     jk = jm - 1
451                     ztrid   (ji,jm,1) =       - zeta_s(ji,jk) *   zkappa_s(ji,jk-1)
452                     ztrid   (ji,jm,2) = 1._wp + zeta_s(ji,jk) * ( zkappa_s(ji,jk-1) + zkappa_s(ji,jk) )
453                     ztrid   (ji,jm,3) =       - zeta_s(ji,jk) *                       zkappa_s(ji,jk)
454                     zindterm(ji,jm)   = ztsold(ji,jk) + zeta_s(ji,jk) * zradab_s(ji,jk)
455                  END DO
456                 
457                  ! case of only one layer in the ice (ice equation is altered)
458                  IF( nlay_i == 1 ) THEN
459                     ztrid   (ji,nlay_s+2,3) = 0._wp
460                     zindterm(ji,nlay_s+2)   = zindterm(ji,nlay_s+2) + zeta_i(ji,1) * zkappa_i(ji,1) * t_bo_1d(ji) 
461                  ENDIF
462                 
463                  IF( t_su_1d(ji) < rt0 ) THEN   !--  case 1 : no surface melting
464                     
465                     jm_min(ji) = 1
466                     jm_max(ji) = nlay_i + nlay_s + 1
467                     
468                     ! surface equation
469                     ztrid   (ji,1,1) = 0._wp
470                     ztrid   (ji,1,2) = zdqns_ice_b(ji) - zg1s * zkappa_s(ji,0)
471                     ztrid   (ji,1,3) =                   zg1s * zkappa_s(ji,0)
472                     zindterm(ji,1)   = zdqns_ice_b(ji) * t_su_1d(ji) - zfnet(ji)
473                     
474                     ! first layer of snow equation
475                     ztrid   (ji,2,1) =       - zeta_s(ji,1) *                    zkappa_s(ji,0) * zg1s
476                     ztrid   (ji,2,2) = 1._wp + zeta_s(ji,1) * ( zkappa_s(ji,1) + zkappa_s(ji,0) * zg1s )
477                     ztrid   (ji,2,3) =       - zeta_s(ji,1) *   zkappa_s(ji,1)
478                     zindterm(ji,2)   = ztsold(ji,1) + zeta_s(ji,1) * zradab_s(ji,1)
479                     
480                  ELSE                            !--  case 2 : surface is melting
481                     !
482                     jm_min(ji) = 2
483                     jm_max(ji) = nlay_i + nlay_s + 1
484                     
485                     ! first layer of snow equation
486                     ztrid   (ji,2,1) = 0._wp
487                     ztrid   (ji,2,2) = 1._wp + zeta_s(ji,1) * ( zkappa_s(ji,1) + zkappa_s(ji,0) * zg1s )
488                     ztrid   (ji,2,3) =       - zeta_s(ji,1) *   zkappa_s(ji,1) 
489                     zindterm(ji,2)   = ztsold(ji,1) + zeta_s(ji,1) * ( zradab_s(ji,1) + zkappa_s(ji,0) * zg1s * t_su_1d(ji) ) 
490                  ENDIF
491                  !                            !---------------------!
492               ELSE                            ! cells without snow  !
493                  !                            !---------------------!
494                  !
495                  IF( t_su_1d(ji) < rt0 ) THEN   !--  case 1 : no surface melting
496                     !
497                     jm_min(ji) = nlay_s + 1
498                     jm_max(ji) = nlay_i + nlay_s + 1
499                     
500                     ! surface equation   
501                     ztrid   (ji,jm_min(ji),1) = 0._wp
502                     ztrid   (ji,jm_min(ji),2) = zdqns_ice_b(ji) - zkappa_i(ji,0) * zg1   
503                     ztrid   (ji,jm_min(ji),3) =                   zkappa_i(ji,0) * zg1
504                     zindterm(ji,jm_min(ji))   = zdqns_ice_b(ji) * t_su_1d(ji) - zfnet(ji)
505                     
506                     ! first layer of ice equation
507                     ztrid   (ji,jm_min(ji)+1,1) =       - zeta_i(ji,1) *                    zkappa_i(ji,0) * zg1
508                     ztrid   (ji,jm_min(ji)+1,2) = 1._wp + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 )
509                     ztrid   (ji,jm_min(ji)+1,3) =       - zeta_i(ji,1) *   zkappa_i(ji,1) 
510                     zindterm(ji,jm_min(ji)+1)   = ztiold(ji,1) + zeta_i(ji,1) * zradab_i(ji,1) 
511                     
512                     ! case of only one layer in the ice (surface & ice equations are altered)
513                     IF( nlay_i == 1 ) THEN
514                        ztrid   (ji,jm_min(ji),1)   = 0._wp
515                        ztrid   (ji,jm_min(ji),2)   = zdqns_ice_b(ji)      -   zkappa_i(ji,0) * 2._wp
516                        ztrid   (ji,jm_min(ji),3)   =                          zkappa_i(ji,0) * 2._wp
517                        ztrid   (ji,jm_min(ji)+1,1) =       - zeta_i(ji,1) *   zkappa_i(ji,0) * 2._wp
518                        ztrid   (ji,jm_min(ji)+1,2) = 1._wp + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2._wp + zkappa_i(ji,1) )
519                        ztrid   (ji,jm_min(ji)+1,3) = 0._wp
520                        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))
521                     ENDIF
522                     
523                  ELSE                            !--  case 2 : surface is melting
524                     
525                     jm_min(ji) = nlay_s + 2
526                     jm_max(ji) = nlay_i + nlay_s + 1
527                     
528                     ! first layer of ice equation
529                     ztrid   (ji,jm_min(ji),1) = 0._wp
530                     ztrid   (ji,jm_min(ji),2) = 1._wp + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 ) 
531                     ztrid   (ji,jm_min(ji),3) =       - zeta_i(ji,1) *   zkappa_i(ji,1)
532                     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)) 
533                     
534                     ! case of only one layer in the ice (surface & ice equations are altered)
535                     IF( nlay_i == 1 ) THEN
536                        ztrid   (ji,jm_min(ji),1) = 0._wp
537                        ztrid   (ji,jm_min(ji),2) = 1._wp + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2._wp + zkappa_i(ji,1) )
538                        ztrid   (ji,jm_min(ji),3) = 0._wp
539                        zindterm(ji,jm_min(ji))   = ztiold(ji,1) + zeta_i(ji,1) * ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) ) &
540                           &                      + t_su_1d(ji) * zeta_i(ji,1) * zkappa_i(ji,0) * 2._wp
541                     ENDIF
542                     
543                  ENDIF
544               ENDIF
545               !
546               zindtbis(ji,jm_min(ji)) = zindterm(ji,jm_min(ji))
547               zdiagbis(ji,jm_min(ji)) = ztrid   (ji,jm_min(ji),2)
548               !
549            END DO
550
551         DO ji = 1, npti
552           IF (to_print(ji) == 10) THEN
553             write(numout,*) 'tri diagonal matrix terms zetai = ',zeta_i(ji,:)
554           END IF
555         END DO
556            !
557            !------------------------------
558            ! 8) tridiagonal system solving
559            !------------------------------
560            ! Solve the tridiagonal system with Gauss elimination method.
561            ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, McGraw-Hill 1984
562            jm_maxt = 0
563            jm_mint = nlay_i+5
564            DO ji = 1, npti
565               jm_mint = MIN(jm_min(ji),jm_mint)
566               jm_maxt = MAX(jm_max(ji),jm_maxt)
567            END DO
568
569            DO jk = jm_mint+1, jm_maxt
570               DO ji = 1, npti
571                  jm = MIN(MAX(jm_min(ji)+1,jk),jm_max(ji))
572                  zdiagbis(ji,jm) = ztrid   (ji,jm,2) - ztrid(ji,jm,1) * ztrid   (ji,jm-1,3) / zdiagbis(ji,jm-1)
573                  zindtbis(ji,jm) = zindterm(ji,jm  ) - ztrid(ji,jm,1) * zindtbis(ji,jm-1  ) / zdiagbis(ji,jm-1)
574               END DO
575            END DO
576
577            ! ice temperatures
578            DO ji = 1, npti
579               ! Variable used after iterations
580               ! Value must be frozen after convergence for MPP independance reason
581               IF ( .NOT. l_T_converged(ji) ) &
582                  t_i_1d(ji,nlay_i) = zindtbis(ji,jm_max(ji)) / zdiagbis(ji,jm_max(ji))
583            END DO
584
585            DO jm = nlay_i + nlay_s, nlay_s + 2, -1
586               DO ji = 1, npti
587                  jk = jm - nlay_s - 1
588                  IF ( .NOT. l_T_converged(ji) ) &
589                     t_i_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,jm)
590               END DO
591            END DO
592
593            DO ji = 1, npti
594               ! Variables used after iterations
595               ! Value must be frozen after convergence for MPP independance reason
596               IF ( .NOT. l_T_converged(ji) ) THEN
597                  ! snow temperatures     
598                  IF( h_s_1d(ji) > 0._wp ) THEN
599                     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)
600                  ENDIF
601                  ! surface temperature
602                  ztsub(ji) = t_su_1d(ji)
603                  IF( t_su_1d(ji) < rt0 ) THEN
604                     t_su_1d(ji) = ( zindtbis(ji,jm_min(ji)) - ztrid(ji,jm_min(ji),3) *  &
605                        &          ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) ) ) / zdiagbis(ji,jm_min(ji))
606                  ENDIF
607               ENDIF
608
609            IF (to_print(ji) == 10 ) THEN
610                write(numout,*)'icethd_zdf_bl99 1: t_s_1d(ji,1), e_s_1d(ji,1) = ',t_s_1d(ji,1), ' ', e_s_1d(ji,1)
611                write(numout,*)'icethd_zdf_bl99 1: zindtbis(ji,nlay_s+1), ztrid(ji,nlay_s+1,3), t_i_1d(ji,1) = ', zindtbis(ji,nlay_s+1), ' ', ztrid(ji,nlay_s+1,3), ' ', t_i_1d(ji,1)
612            END IF
613
614            END DO
615         DO ji = 1, npti
616           IF (to_print(ji) == 10) THEN
617             write(numout,*) 'icethd_zdf_bl99 18 iter: t_i_1d(ji,1), e_i_1d(ji,1) = ',t_i_1d(ji,1), ' ', e_i_1d(ji,1)
618           END IF
619         END DO
620            !
621            !--------------------------------------------------------------
622            ! 9) Has the scheme converged?, end of the iterative procedure
623            !--------------------------------------------------------------
624            ! check that nowhere it has started to melt
625            ! zdti_max is a measure of error, it has to be under zdti_bnd
626
627            DO ji = 1, npti
628
629               zdti_max = 0._wp
630
631               IF ( .NOT. l_T_converged(ji) ) THEN
632                  t_su_1d(ji) = MAX( MIN( t_su_1d(ji) , rt0 ) , rt0 - 100._wp )
633                  zdti_max    = MAX( zdti_max, ABS( t_su_1d(ji) - ztsub(ji) ) )
634
635                  t_s_1d(ji,1:nlay_s) = MAX( MIN( t_s_1d(ji,1:nlay_s), rt0 ), rt0 - 100._wp )
636                  zdti_max = MAX ( zdti_max , MAXVAL( ABS( t_s_1d(ji,1:nlay_s) - ztsb(ji,1:nlay_s) ) ) )
637
638                  DO jk = 1, nlay_i
639                     ztmelts       = -rTmlt * sz_i_1d(ji,jk) + rt0
640                     t_i_1d(ji,jk) =  MAX( MIN( t_i_1d(ji,jk), ztmelts ), rt0 - 100._wp )
641                     zdti_max      =  MAX( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) )
642                  END DO
643
644                  IF ( zdti_max < zdti_bnd ) l_T_converged(ji) = .TRUE.
645
646               ENDIF
647
648            END DO
649
650         DO ji = 1, npti
651           IF (to_print(ji) == 10) THEN
652             write(numout,*) 'icethd_zdf_bl99 19 iter: t_i_1d(ji,1), e_i_1d(ji,1) = ',t_i_1d(ji,1), ' ', e_i_1d(ji,1)
653           END IF
654         END DO
655
656         !----------------------------------------!
657         !                                        !
658         !      Conduction flux is on             !
659         !                                        !
660         !----------------------------------------!
661         !
662         ELSEIF( k_cnd == np_cnd_ON ) THEN
663            !
664            ! ==> we use a modified BL99 solver with conduction flux (qcn_ice) as forcing term
665         DO ji = 1, npti
666           IF (to_print(ji) == 10) THEN
667             write(numout,*) 'icethd_zdf_bl99 7.1 iter: t_i_1d(ji,1), e_i_1d(ji,1) = ',t_i_1d(ji,1), ' ', e_i_1d(ji,1)
668           END IF
669         END DO
670            !
671            !----------------------------
672            ! 7) tridiagonal system terms
673            !----------------------------
674            ! layer denotes the number of the layer in the snow or in the ice
675            ! jm denotes the reference number of the equation in the tridiagonal
676            ! system, terms of tridiagonal system are indexed as following :
677            ! 1 is subdiagonal term, 2 is diagonal and 3 is superdiagonal one
678
679            ! ice interior terms (top equation has the same form as the others)
680            ztrid   (1:npti,:,:) = 0._wp
681            zindterm(1:npti,:)   = 0._wp
682            zindtbis(1:npti,:)   = 0._wp
683            zdiagbis(1:npti,:)   = 0._wp
684
685            DO jm = nlay_s + 2, nlay_s + nlay_i 
686               DO ji = 1, npti
687                  jk = jm - nlay_s - 1
688                  ztrid   (ji,jm,1) =       - zeta_i(ji,jk) *   zkappa_i(ji,jk-1)
689                  ztrid   (ji,jm,2) = 1._wp + zeta_i(ji,jk) * ( zkappa_i(ji,jk-1) + zkappa_i(ji,jk) )
690                  ztrid   (ji,jm,3) =       - zeta_i(ji,jk) *                       zkappa_i(ji,jk)
691                  zindterm(ji,jm)   = ztiold(ji,jk) + zeta_i(ji,jk) * zradab_i(ji,jk)
692               END DO
693            ENDDO
694
695            jm =  nlay_s + nlay_i + 1
696            DO ji = 1, npti
697               ! ice bottom term
698               ztrid   (ji,jm,1) =       - zeta_i(ji,nlay_i) *   zkappa_i(ji,nlay_i-1)   
699               ztrid   (ji,jm,2) = 1._wp + zeta_i(ji,nlay_i) * ( zkappa_i(ji,nlay_i-1) + zkappa_i(ji,nlay_i) * zg1 )
700               ztrid   (ji,jm,3) = 0._wp
701               zindterm(ji,jm)   = ztiold(ji,nlay_i) + zeta_i(ji,nlay_i) *  &
702                  &              ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i) * zg1 * t_bo_1d(ji) ) 
703            ENDDO
704
705         DO ji = 1, npti
706           IF (to_print(ji) == 10 .AND. iconv == 1) THEN
707             jm =  nlay_s + nlay_i + 1
708             write(numout,*) 'icethd_zdf_bl99: zindterm(ji,jm), ztiold(ji,nlay_i), zeta_i(ji,nlay_i), zradab_i(ji,nlay_i) = ',zindterm(ji,jm), ' ', ztiold(ji,nlay_i), ' ', zeta_i(ji,nlay_i), ' ', zradab_i(ji,nlay_i)
709             write(numout,*) 'icethd_zdf_bl99: zkappa_i(ji,nlay_i), zg1, t_bo_1d(ji) = ',zkappa_i(ji,nlay_i), ' ', zg1, ' ', t_bo_1d(ji)
710           END IF
711         END DO
712
713            DO ji = 1, npti
714               !                               !---------------------!
715               IF( h_s_1d(ji) > 0._wp ) THEN   !  snow-covered cells !
716                  !                            !---------------------!
717                  ! snow interior terms (bottom equation has the same form as the others)
718                  DO jm = 3, nlay_s + 1
719                     jk = jm - 1
720                     ztrid   (ji,jm,1) =       - zeta_s(ji,jk) *   zkappa_s(ji,jk-1)
721                     ztrid   (ji,jm,2) = 1._wp + zeta_s(ji,jk) * ( zkappa_s(ji,jk-1) + zkappa_s(ji,jk) )
722                     ztrid   (ji,jm,3) =       - zeta_s(ji,jk) *                       zkappa_s(ji,jk)
723                     zindterm(ji,jm)   = ztsold(ji,jk) + zeta_s(ji,jk) * zradab_s(ji,jk)
724                  END DO
725                 
726                  ! case of only one layer in the ice (ice equation is altered)
727                  IF ( nlay_i == 1 ) THEN
728                     ztrid   (ji,nlay_s+2,3) = 0._wp
729                     zindterm(ji,nlay_s+2)   = zindterm(ji,nlay_s+2) + zeta_i(ji,1) * zkappa_i(ji,1) * t_bo_1d(ji) 
730                  ENDIF
731                 
732                  jm_min(ji) = 2
733                  jm_max(ji) = nlay_i + nlay_s + 1
734                 
735                  ! first layer of snow equation
736                  ztrid   (ji,2,1) = 0._wp
737                  ztrid   (ji,2,2) = 1._wp + zeta_s(ji,1) * zkappa_s(ji,1)
738                  ztrid   (ji,2,3) =       - zeta_s(ji,1) * zkappa_s(ji,1) 
739                  zindterm(ji,2)   = ztsold(ji,1) + zeta_s(ji,1) * ( zradab_s(ji,1) + qcn_ice_1d(ji) ) 
740
741                  IF (to_print(ji) == 10) THEN
742                     write(numout,*) 'ztrid for snow = ',ztrid   (ji,2,:)
743                     write(numout,*) 'zindterm for snow = ',zindterm(ji,2)
744                     write(numout,*) 'zindterm using degC = ',zindterm(ji,2) - rt0
745                     write(numout,*) 'nlay_s, nlay_i = ',nlay_s, ' ', nlay_i
746                  ENDIF
747
748                 
749                  !                            !---------------------!
750               ELSE                            ! cells without snow  !
751                  !                            !---------------------!
752                  jm_min(ji) = nlay_s + 2
753                  jm_max(ji) = nlay_i + nlay_s + 1
754                 
755                  ! first layer of ice equation
756                  ztrid   (ji,jm_min(ji),1) = 0._wp
757                  ztrid   (ji,jm_min(ji),2) = 1._wp + zeta_i(ji,1) * zkappa_i(ji,1) 
758                  ztrid   (ji,jm_min(ji),3) =       - zeta_i(ji,1) * zkappa_i(ji,1)
759                  zindterm(ji,jm_min(ji))   = ztiold(ji,1) + zeta_i(ji,1) * ( zradab_i(ji,1) + qcn_ice_1d(ji) )
760
761                  IF (to_print(ji) == 10 .AND. iconv == 1) THEN
762                     write(numout,*) 'zeta_i(ji,1), zkappa_i(ji,1), ztiold(ji,1), qcn_ice_1d(ji) = ',zeta_i(ji,1), ' ', zkappa_i(ji,1), ' ', ztiold(ji,1), ' ', qcn_ice_1d(ji)
763                  ENDIF                 
764                 
765                  ! case of only one layer in the ice (surface & ice equations are altered)
766                  IF( nlay_i == 1 ) THEN
767                     ztrid   (ji,jm_min(ji),1) = 0._wp
768                     ztrid   (ji,jm_min(ji),2) = 1._wp + zeta_i(ji,1) * zkappa_i(ji,1)
769                     ztrid   (ji,jm_min(ji),3) = 0._wp
770                     zindterm(ji,jm_min(ji))   = ztiold(ji,1) + zeta_i(ji,1) *  &
771                        &                                     ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) + qcn_ice_1d(ji) )
772                  ENDIF               
773                 
774               ENDIF
775
776                  IF (to_print(ji) == 10) THEN
777                     write(numout,*)'after minmax: h_s_1d(ji), jm_min(ji), nlay_s = ',h_s_1d(ji), ' ', jm_min(ji), ' ', nlay_s
778                  ENDIF   
779
780               !
781               zindtbis(ji,jm_min(ji)) = zindterm(ji,jm_min(ji))
782               zdiagbis(ji,jm_min(ji)) = ztrid   (ji,jm_min(ji),2)
783               !
784            END DO
785            !
786            !------------------------------
787            ! 8) tridiagonal system solving
788            !------------------------------
789            ! Solve the tridiagonal system with Gauss elimination method.
790            ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, McGraw-Hill 1984
791            jm_maxt = 0
792            jm_mint = nlay_i+5
793            DO ji = 1, npti
794               jm_mint = MIN(jm_min(ji),jm_mint)
795               jm_maxt = MAX(jm_max(ji),jm_maxt)
796            END DO
797           
798            DO jk = jm_mint+1, jm_maxt
799               DO ji = 1, npti
800
801 
802                  jm = MIN(MAX(jm_min(ji)+1,jk),jm_max(ji))
803                  IF (to_print(ji) == 10) THEN
804                     write(numout,*)'before calc zdiagbis: jk, jm, zdiagbis(ji,jm-1) = ',jk, ' ', jm, ' ', zdiagbis(ji,jm-1)
805                     write(numout,*)'before calc zdiagbis:ztrid(ji,jm,2), ztrid(ji,jm,1), ztrid(ji,jm-1,3) = ', ztrid(ji,jm,2), ' ', ztrid(ji,jm,1), ' ', ztrid(ji,jm-1,3)
806                  ENDIF
807
808                  zdiagbis(ji,jm) = ztrid   (ji,jm,2) - ztrid(ji,jm,1) * ztrid   (ji,jm-1,3) / zdiagbis(ji,jm-1)
809                  zindtbis(ji,jm) = zindterm(ji,jm)   - ztrid(ji,jm,1) * zindtbis(ji,jm-1)   / zdiagbis(ji,jm-1)
810               END DO
811            END DO
812           
813            ! ice temperatures
814            DO ji = 1, npti
815               ! Variable used after iterations
816               ! Value must be frozen after convergence for MPP independance reason
817               IF ( .NOT. l_T_converged(ji) ) &
818                  t_i_1d(ji,nlay_i) = zindtbis(ji,jm_max(ji)) / zdiagbis(ji,jm_max(ji))
819            END DO
820
821            DO jm = nlay_i + nlay_s, nlay_s + 2, -1
822               DO ji = 1, npti
823                  IF ( .NOT. l_T_converged(ji) ) THEN
824                     jk = jm - nlay_s - 1
825                     t_i_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,jm)
826                  ENDIF
827               END DO
828            END DO
829           
830            ! snow temperatures     
831            DO ji = 1, npti
832               ! Variable used after iterations
833               ! Value must be frozen after convergence for MPP independance reason
834               IF ( .NOT. l_T_converged(ji) ) THEN
835                  IF( h_s_1d(ji) > 0._wp ) THEN
836                     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)
837                  ENDIF
838               ENDIF
839
840              IF (to_print(ji) == 10) THEN
841                write(numout,*)'icethd_zdf_bl99 2: t_s_1d(ji,1), e_s_1d(ji,1) = ',t_s_1d(ji,1), ' ', e_s_1d(ji,1)
842                write(numout,*)'icethd_zdf_bl99 2: zindtbis(ji,nlay_s+1), ztrid(ji,nlay_s+1,3), t_i_1d(ji,1) = ', zindtbis(ji,nlay_s+1), ' ', ztrid(ji,nlay_s+1,3), ' ', t_i_1d(ji,1)
843                write(numout,*)'icethd_zdf_bl99 2: zdiagbis(ji,nlay_s+1), h_s_1d(ji) = ',zdiagbis(ji,nlay_s+1), ' ', h_s_1d(ji)
844              END IF
845            END DO
846         DO ji = 1, npti
847           IF (to_print(ji) == 10) THEN
848             write(numout,*) 'icethd_zdf_bl99 8 iter: t_i_1d(ji,1), e_i_1d(ji,1) = ',t_i_1d(ji,1), ' ', e_i_1d(ji,1)
849           END IF
850         END DO
851            !
852            !--------------------------------------------------------------
853            ! 9) Has the scheme converged?, end of the iterative procedure
854            !--------------------------------------------------------------
855            ! check that nowhere it has started to melt
856            ! zdti_max is a measure of error, it has to be under zdti_bnd
857
858            DO ji = 1, npti
859
860               zdti_max = 0._wp
861
862               IF ( .NOT. l_T_converged(ji) ) THEN
863                  ! t_s
864                  t_s_1d(ji,1:nlay_s) = MAX( MIN( t_s_1d(ji,1:nlay_s), rt0 ), rt0 - 100._wp )
865                  zdti_max = MAX ( zdti_max , MAXVAL( ABS( t_s_1d(ji,1:nlay_s) - ztsb(ji,1:nlay_s) ) ) )
866                  ! t_i
867                  DO jk = 1, nlay_i
868                     ztmelts       = -rTmlt * sz_i_1d(ji,jk) + rt0 
869                     t_i_1d(ji,jk) =  MAX( MIN( t_i_1d(ji,jk), ztmelts ), rt0 - 100._wp )
870                     zdti_max      =  MAX ( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) )
871                  END DO
872
873                  IF ( zdti_max < zdti_bnd ) l_T_converged(ji) = .TRUE.
874
875               ENDIF
876
877            END DO
878
879         ENDIF ! k_cnd
880
881         DO ji = 1, npti
882           IF (iconv == 1 .AND. to_print(ji) == 10) THEN
883             write(numout,*) 'matrix1 befroe diagonal = ',ztrid(ji,1,1),' ',ztrid(ji,2,1),' ',ztrid(ji,3,1),' ',ztrid(ji,4,1),' ',ztrid(ji,5,1),' ',ztrid(ji,6,1)
884             write(numout,*) 'matrix1 on diagonal = ',ztrid(ji,1,2),' ',ztrid(ji,2,2),' ',ztrid(ji,3,2),' ',ztrid(ji,4,2),' ',ztrid(ji,5,2),' ',ztrid(ji,6,2)
885             write(numout,*) 'matrix1 after diagonal = ',ztrid(ji,1,3),' ',ztrid(ji,2,3),' ',ztrid(ji,3,3),' ',ztrid(ji,4,3),' ',ztrid(ji,5,3),' ',ztrid(ji,6,3)
886             write(numout,*) 'rhs1 for all = ',zindterm(ji,1),' ',zindterm(ji,2),' ',zindterm(ji,3),' ',zindterm(ji,4),' ',zindterm(ji,5),' ',zindterm(ji,6)
887             write(numout,*) 'rhs1 using degC = ',zindterm(ji,1)-rt0,' ',zindterm(ji,2)-rt0,' ',zindterm(ji,3)-rt0,' ',zindterm(ji,4)-rt0,' ',zindterm(ji,5)-rt0,' ',zindterm(ji,6)-rt0
888           END IF
889         END DO
890         
891      END DO  ! End of the do while iterative procedure
892
893      DO ji = 1, npti
894        IF (to_print(ji) == 10) THEN
895          write(numout,*) 'matrix2 befroe diagonal = ',ztrid(ji,1,1),' ',ztrid(ji,2,1),' ',ztrid(ji,3,1),' ',ztrid(ji,4,1),' ',ztrid(ji,5,1),' ',ztrid(ji,6,1)
896          write(numout,*) 'matrix2 on diagonal = ',ztrid(ji,1,2),' ',ztrid(ji,2,2),' ',ztrid(ji,3,2),' ',ztrid(ji,4,2),' ',ztrid(ji,5,2),' ',ztrid(ji,6,2)
897          write(numout,*) 'matrix2 after diagonal = ',ztrid(ji,1,3),' ',ztrid(ji,2,3),' ',ztrid(ji,3,3),' ',ztrid(ji,4,3),' ',ztrid(ji,5,3),' ',ztrid(ji,6,3)
898          write(numout,*) 'rhs2 for all = ',zindterm(ji,1),' ',zindterm(ji,2),' ',zindterm(ji,3),' ',zindterm(ji,4),' ',zindterm(ji,5),' ',zindterm(ji,6)
899          write(numout,*) 'rhs2 using degC = ',zindterm(ji,1)-rt0,' ',zindterm(ji,2)-rt0,' ',zindterm(ji,3)-rt0,' ',zindterm(ji,4)-rt0,' ',zindterm(ji,5)-rt0,' ',zindterm(ji,6)-rt0
900          write(numout,*) 'number of iterations iconv = ',iconv
901        END IF
902      END DO
903
904     
905      IF( ln_icectl .AND. lwp ) THEN
906         WRITE(numout,*) ' zdti_max : ', zdti_max
907         WRITE(numout,*) ' iconv    : ', iconv
908      ENDIF
909     
910      !
911      !-----------------------------
912      ! 10) Fluxes at the interfaces
913      !-----------------------------
914      !
915      ! --- calculate conduction fluxes (positive downward)
916      !     bottom ice conduction flux
917      DO ji = 1, npti
918         qcn_ice_bot_1d(ji) =  - zkappa_i(ji,nlay_i) * zg1  * ( t_bo_1d(ji ) - t_i_1d (ji,nlay_i) )
919
920         IF (to_print(ji) == 1) THEN
921           write(numout,*) 'icethd_zdf_bl99: qcn_ice_bot_1d(ji), zkappa_i(ji,nlay_i), zg1, t_bo_1d(ji ), t_i_1d (ji,nlay_i) = ',qcn_ice_bot_1d(ji), zkappa_i(ji,nlay_i), zg1, t_bo_1d(ji ), t_i_1d (ji,nlay_i)
922         ENDIF
923      END DO
924      !     surface ice conduction flux
925      IF( k_cnd == np_cnd_OFF .OR. k_cnd == np_cnd_EMU ) THEN
926         !
927         DO ji = 1, npti
928            qcn_ice_top_1d(ji) =  -           isnow(ji)   * zkappa_s(ji,0) * zg1s * ( t_s_1d(ji,1) - t_su_1d(ji) )  &
929               &                  - ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1  * ( t_i_1d(ji,1) - t_su_1d(ji) )
930         END DO
931         !
932      ELSEIF( k_cnd == np_cnd_ON ) THEN
933         !
934         DO ji = 1, npti
935            qcn_ice_top_1d(ji) = qcn_ice_1d(ji)
936         END DO
937         !
938      ENDIF
939      !     surface ice temperature
940      IF( k_cnd == np_cnd_ON .AND. ln_cndemulate ) THEN
941         !
942         DO ji = 1, npti
943            t_su_1d(ji) = (  qcn_ice_top_1d(ji) &            ! calculate surface temperature
944               &           +           isnow(ji)   * zkappa_s(ji,0) * zg1s * t_s_1d(ji,1) &
945               &           + ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1  * t_i_1d(ji,1) &
946               &          ) / MAX( epsi10, isnow(ji) * zkappa_s(ji,0) * zg1s + ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1 )
947            t_su_1d(ji) = MAX( MIN( t_su_1d(ji), rt0 ), rt0 - 100._wp )  ! cap t_su
948         END DO
949         !
950      ENDIF
951      !
952      ! --- Diagnose the heat loss due to changing non-solar / conduction flux --- !
953      !
954      IF( k_cnd == np_cnd_OFF .OR. k_cnd == np_cnd_EMU ) THEN
955         !
956         DO ji = 1, npti
957            hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) - ( qns_ice_1d(ji) - zqns_ice_b(ji) ) * a_i_1d(ji) 
958         END DO
959         !
960      ENDIF
961     
962      !
963      ! --- Diagnose the heat loss due to non-fully converged temperature solution (should not be above 10-4 W-m2) --- !
964      !
965      IF( k_cnd == np_cnd_OFF .OR. k_cnd == np_cnd_ON ) THEN
966
967         DO ji = 1, npti
968           IF (to_print(ji) == 10) THEN
969             write(numout,*) 'before ice_var_enthalpy: t_i_1d(ji,1), e_i_1d(ji,1) = ',t_i_1d(ji,1), ' ', e_i_1d(ji,1)
970           END IF
971         END DO
972         
973         CALL ice_var_enthalpy   
974
975         DO ji = 1, npti
976           IF (to_print(ji) == 10) THEN
977             write(numout,*) 'after ice_var_enthalpy: t_i_1d(ji,1), e_i_1d(ji,1) = ',t_i_1d(ji,1), ' ', e_i_1d(ji,1)
978           END IF
979         END DO   
980         
981         ! zhfx_err = correction on the diagnosed heat flux due to non-convergence of the algorithm used to solve heat equation
982         DO ji = 1, npti
983            zdq = - zq_ini(ji) + ( SUM( e_i_1d(ji,1:nlay_i) ) * h_i_1d(ji) * r1_nlay_i +  &
984               &                   SUM( e_s_1d(ji,1:nlay_s) ) * h_s_1d(ji) * r1_nlay_s )
985           
986            IF( k_cnd == np_cnd_OFF ) THEN
987               
988               IF( t_su_1d(ji) < rt0 ) THEN  ! case T_su < 0degC
989                  zhfx_err = ( qns_ice_1d(ji)     + qsr_ice_1d(ji)     - zradtr_i(ji,nlay_i) - qcn_ice_bot_1d(ji)  &
990                     &       + zdq * r1_rdtice ) * a_i_1d(ji)
991               ELSE                          ! case T_su = 0degC
992                  zhfx_err = ( qcn_ice_top_1d(ji) + qtr_ice_top_1d(ji) - zradtr_i(ji,nlay_i) - qcn_ice_bot_1d(ji)  &
993                     &       + zdq * r1_rdtice ) * a_i_1d(ji)
994               ENDIF
995               
996            ELSEIF( k_cnd == np_cnd_ON ) THEN
997           
998               zhfx_err    = ( qcn_ice_top_1d(ji) + qtr_ice_top_1d(ji) - zradtr_i(ji,nlay_i) - qcn_ice_bot_1d(ji)  &
999                  &          + zdq * r1_rdtice ) * a_i_1d(ji)
1000           
1001            ENDIF
1002            !
1003            ! total heat sink to be sent to the ocean
1004            hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) + zhfx_err
1005            !
1006            ! hfx_dif = Heat flux diagnostic of sensible heat used to warm/cool ice in W.m-2   
1007            hfx_dif_1d(ji) = hfx_dif_1d(ji) - zdq * r1_rdtice * a_i_1d(ji)
1008
1009            IF (to_print(ji) == 10) THEN
1010               write(numout,*)'icethd_zdf_bl99: qcn_ice_top_1d(ji) = ',qcn_ice_top_1d(ji)
1011               write(numout,*)'icethd_zdf_bl99: qcn_ice_1d(ji) = ',qcn_ice_1d(ji)
1012               write(numout,*)'icethd_zdf_bl99: hfx_err_dif_1d(ji) = ', hfx_err_dif_1d(ji)
1013               write(numout,*)'icethd_zdf_bl99: hfx_dif_1d(ji) = ',hfx_dif_1d(ji)
1014            ENDIF
1015
1016            !
1017         END DO
1018         !
1019      ENDIF
1020      !
1021      !--------------------------------------------------------------------
1022      ! 11) reset inner snow and ice temperatures, update conduction fluxes
1023      !--------------------------------------------------------------------
1024      ! effective conductivity and 1st layer temperature (needed by Met Office)
1025      DO ji = 1, npti
1026         IF( h_i_1d(ji) > 0.1_wp ) THEN
1027            cnd_ice_1d(ji) = 2._wp * zkappa_i(ji,0)
1028         ELSE
1029            cnd_ice_1d(ji) = 2._wp * ztcond_i(ji,0) * 10._wp
1030         ENDIF
1031         t1_ice_1d(ji) = isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1)
1032
1033         IF (to_print(ji) == 10) THEN
1034            write(numout,*)'icethd_zdf_bl99: cnd_ice_1d(ji), h_s_1d(ji), zkappa_s(ji,0) = ',cnd_ice_1d(ji), ' ', h_s_1d(ji), ' ', zkappa_s(ji,0)
1035            write(numout,*)'icethd_zdf_bl99: ztcond_i(ji,0), h_i_1d(ji), zkappa_i(ji,0) = ',ztcond_i(ji,0), ' ', h_i_1d(ji), ' ', zkappa_i(ji,0)
1036            write(numout,*)'icethd_zdf_bl99: t1_ice_1d(ji), isnow(ji), t_s_1d(ji,1), t_i_1d(ji,1) = ',t1_ice_1d(ji), ' ', isnow(ji), ' ', t_s_1d(ji,1), ' ', t_i_1d(ji,1)
1037         ENDIF
1038      END DO
1039      !
1040      IF( k_cnd == np_cnd_EMU ) THEN
1041         ! Restore temperatures to their initial values
1042         t_s_1d    (1:npti,:) = ztsold        (1:npti,:)
1043         t_i_1d    (1:npti,:) = ztiold        (1:npti,:)
1044         qcn_ice_1d(1:npti)   = qcn_ice_top_1d(1:npti)
1045
1046         !!clem
1047         ! remettre t_su_1d, qns_ice_1d et dqns_ice_1d comme avant puisqu'on devrait faire comme si on avant conduction = input
1048         !clem
1049      ENDIF
1050      !
1051      ! --- SIMIP diagnostics
1052      !
1053      DO ji = 1, npti         
1054         !--- Snow-ice interfacial temperature (diagnostic SIMIP)
1055         zfac = rn_cnd_s * zh_i(ji) + ztcond_i(ji,1) * zh_s(ji)
1056         IF( h_s_1d(ji) >= zhs_min ) THEN
1057            t_si_1d(ji) = ( rn_cnd_s       * zh_i(ji) * t_s_1d(ji,1) +   &
1058               &            ztcond_i(ji,1) * zh_s(ji) * t_i_1d(ji,1) ) / MAX( epsi10, zfac )
1059         ELSE
1060            t_si_1d(ji) = t_su_1d(ji)
1061         ENDIF
1062      END DO
1063      !
1064   END SUBROUTINE ice_thd_zdf_BL99
1065
1066#else
1067   !!----------------------------------------------------------------------
1068   !!   Default option       Dummy Module             No SI3 sea-ice model
1069   !!----------------------------------------------------------------------
1070#endif
1071
1072   !!======================================================================
1073END MODULE icethd_zdf_BL99
Note: See TracBrowser for help on using the repository browser.