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

Last change on this file since 12374 was 12374, checked in by dancopsey, 10 months ago

Get skin temperature to be passed from the atmosphere

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