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/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE – NEMO

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

Last change on this file since 12369 was 12369, checked in by dancopsey, 4 years ago

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