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.F90 in branches/UKMO/dev_r8183_ICEMODEL_svn_removed/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/UKMO/dev_r8183_ICEMODEL_svn_removed/NEMOGCM/NEMO/LIM_SRC_3/icethd_zdf.F90 @ 8752

Last change on this file since 8752 was 8752, checked in by dancopsey, 5 years ago

Merged in main ICEMODEL branch (branches/2017/dev_r8183_ICEMODEL) from revision 8587 to 8726.

File size: 47.7 KB
Line 
1MODULE icethd_zdf
2   !!======================================================================
3   !!                       ***  MODULE icethd_zdf ***
4   !!   sea-ice: vertical heat diffusion in sea ice (computation of temperatures)
5   !!======================================================================
6   !! History :  LIM  ! 02-2003 (M. Vancoppenolle) original 1D code
7   !!                 ! 06-2005 (M. Vancoppenolle) 3d version
8   !!                 ! 11-2006 (X Fettweis) Vectorization by Xavier
9   !!                 ! 04-2007 (M. Vancoppenolle) Energy conservation
10   !!            4.0  ! 2011-02 (G. Madec) dynamical allocation
11   !!             -   ! 2012-05 (C. Rousset) add penetration solar flux
12   !!----------------------------------------------------------------------
13#if defined key_lim3
14   !!----------------------------------------------------------------------
15   !!   'key_lim3'                                       ESIM sea-ice model
16   !!----------------------------------------------------------------------
17   USE dom_oce        ! ocean space and time domain
18   USE phycst         ! physical constants (ocean directory)
19   USE ice            ! sea-ice: variables
20   USE ice1D          ! sea-ice: thermodynamics variables
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        ! called by icethd
30   PUBLIC   ice_thd_zdf_init   ! called by icestp
31
32   !!** namelist (namthd_zdf) **
33   LOGICAL  ::   ln_zdf_BL99      ! Heat diffusion follows Bitz and Lipscomb (1999)
34   LOGICAL  ::   ln_cndi_U64      ! thermal conductivity: Untersteiner (1964)
35   LOGICAL  ::   ln_cndi_P07      ! thermal conductivity: Pringle et al (2007)
36   REAL(wp) ::   rn_kappa_i       ! coef. for the extinction of radiation Grenfell et al. (2006) [1/m]
37   REAL(wp), PUBLIC ::   rn_cnd_s ! thermal conductivity of the snow [W/m/K]
38
39   INTEGER  ::   nice_zdf         ! Choice of the type of vertical heat diffusion formulation
40   !                                           ! associated indices:
41   INTEGER, PARAMETER ::   np_BL99    = 1      ! Bitz and Lipscomb (1999)
42
43   INTEGER , PARAMETER ::   np_zdf_jules_OFF    = 0   !   compute all temperatures from qsr and qns
44   INTEGER , PARAMETER ::   np_zdf_jules_SND    = 1   !   compute conductive heat flux and surface temperature from qsr and qns
45   INTEGER , PARAMETER ::   np_zdf_jules_RCV    = 2   !   compute snow and inner ice temperatures from qcnd
46   
47   !!----------------------------------------------------------------------
48   !! NEMO/ICE 4.0 , NEMO Consortium (2017)
49   !! $Id: icethd_zdf.F90 8420 2017-08-08 12:18:46Z clem $
50   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
51   !!----------------------------------------------------------------------
52CONTAINS
53
54      SUBROUTINE ice_thd_zdf
55     
56      !!
57      !!-------------------------------------------------------------------
58      !!                ***  ROUTINE ice_thd_zdf  ***
59      !! ** Purpose :
60      !!           This chooses between the appropriate routine for the
61      !!           computation of vertical diffusion
62      !!
63      !!-------------------------------------------------------------------
64      !!
65   
66      SELECT CASE ( nice_zdf )      ! Choose the vertical heat diffusion solver
67     
68                                       !-------------     
69         CASE( np_BL99 )               ! BL99 solver
70                                       !-------------
71         
72            IF ( nice_jules == np_jules_OFF ) THEN   ! No Jules coupler     
73
74               CALL ice_thd_zdf_BL99 ( np_zdf_jules_OFF )
75                                       
76            ENDIF
77           
78            IF ( nice_jules == np_jules_EMULE ) THEN   ! Jules coupler is emulated
79           
80               CALL ice_thd_zdf_BL99 ( np_zdf_jules_SND )
81               CALL ice_thd_zdf_BL99 ( np_zdf_jules_RCV )
82                                                   
83            ENDIF
84
85            IF ( nice_jules == np_jules_ACTIVE ) THEN   ! Jules coupler is emulated
86           
87               CALL ice_thd_zdf_BL99 ( np_zdf_jules_RCV )
88                                                   
89            ENDIF
90           
91      END SELECT
92     
93   END SUBROUTINE ice_thd_zdf
94   
95   
96   
97   SUBROUTINE ice_thd_zdf_BL99(k_jules)
98      !!-------------------------------------------------------------------
99      !!                ***  ROUTINE ice_thd_zdf_BL99  ***
100      !! ** Purpose :
101      !!           This routine determines the time evolution of snow and sea-ice
102      !!           temperature profiles, using the original Bitz and Lipscomb (1999) algorithm
103      !!
104      !! ** Method  :
105      !!           This is done by solving the heat equation diffusion with
106      !!           a Neumann boundary condition at the surface and a Dirichlet one
107      !!           at the bottom. Solar radiation is partially absorbed into the ice.
108      !!           The specific heat and thermal conductivities depend on ice salinity
109      !!           and temperature to take into account brine pocket melting. The
110      !!           numerical
111      !!           scheme is an iterative Crank-Nicolson on a non-uniform multilayer grid
112      !!           in the ice and snow system.
113      !!
114      !!           The successive steps of this routine are
115      !!           1.  initialization of ice-snow layers thicknesses
116      !!           2.  Internal absorbed and transmitted radiation
117      !!           Then iterative procedure begins
118      !!           3.  Thermal conductivity
119      !!           4.  Kappa factors
120      !!           5.  specific heat in the ice
121      !!           6.  eta factors
122      !!           7.  surface flux computation
123      !!           8.  tridiagonal system terms
124      !!           9.  solving the tridiagonal system with Gauss elimination
125      !!           Iterative procedure ends according to a criterion on evolution
126      !!           of temperature
127      !!           10. Fluxes at the interfaces
128      !!
129      !! ** Inputs / Ouputs : (global commons)
130      !!           surface temperature : t_su_1d
131      !!           ice/snow temperatures   : t_i_1d, t_s_1d
132      !!           ice salinities          : sz_i_1d
133      !!           number of layers in the ice/snow: nlay_i, nlay_s
134      !!           total ice/snow thickness : h_i_1d, h_s_1d
135      !!-------------------------------------------------------------------
136      INTEGER, INTENT(in) ::   k_jules     ! Jules coupling (0=OFF, 1=RECEIVE, 2=SEND)
137     
138      INTEGER ::   ji, jk         ! spatial loop index
139      INTEGER ::   jm             ! current reference number of equation
140      INTEGER ::   jm_mint, jm_maxt
141      INTEGER ::   iconv          ! number of iterations in iterative procedure
142      INTEGER ::   iconv_max = 50 ! max number of iterations in iterative procedure
143     
144      INTEGER, DIMENSION(jpij) ::   jm_min    ! reference number of top equation
145      INTEGER, DIMENSION(jpij) ::   jm_max    ! reference number of bottom equation
146     
147      REAL(wp) ::   zg1s      =  2._wp        ! for the tridiagonal system
148      REAL(wp) ::   zg1       =  2._wp        !
149      REAL(wp) ::   zgamma    =  18009._wp    ! for specific heat
150      REAL(wp) ::   zbeta     =  0.117_wp     ! for thermal conductivity (could be 0.13)
151      REAL(wp) ::   zraext_s  =  10._wp       ! extinction coefficient of radiation in the snow
152      REAL(wp) ::   zkimin    =  0.10_wp      ! minimum ice thermal conductivity
153      REAL(wp) ::   ztsu_err  =  1.e-5_wp     ! range around which t_su is considered at 0C
154      REAL(wp) ::   zdti_bnd  =  1.e-4_wp     ! maximal authorized error on temperature
155      REAL(wp) ::   ztmelt_i                  ! ice melting temperature
156      REAL(wp) ::   zdti_max                  ! current maximal error on temperature
157      REAL(wp) ::   zcpi                      ! Ice specific heat
158      REAL(wp) ::   zhfx_err, zdq             ! diag errors on heat
159      REAL(wp) ::   zfac                      ! dummy factor
160     
161      REAL(wp), DIMENSION(jpij)     ::   isnow       ! switch for presence (1) or absence (0) of snow
162      REAL(wp), DIMENSION(jpij)     ::   ztsub       ! surface temperature at previous iteration
163      REAL(wp), DIMENSION(jpij)     ::   zh_i, z1_h_i ! ice layer thickness
164      REAL(wp), DIMENSION(jpij)     ::   zh_s, z1_h_s ! snow layer thickness
165      REAL(wp), DIMENSION(jpij)     ::   zqns_ice_b  ! solar radiation absorbed at the surface
166      REAL(wp), DIMENSION(jpij)     ::   zfnet       ! surface flux function
167      REAL(wp), DIMENSION(jpij)     ::   zdqns_ice_b ! derivative of the surface flux function
168     
169      REAL(wp), DIMENSION(jpij       )     ::   ztsuold     ! Old surface temperature in the ice
170      REAL(wp), DIMENSION(jpij,nlay_i)     ::   ztiold      ! Old temperature in the ice
171      REAL(wp), DIMENSION(jpij,nlay_s)     ::   ztsold      ! Old temperature in the snow
172      REAL(wp), DIMENSION(jpij,nlay_i)     ::   ztib        ! Temporary temperature in the ice to check the convergence
173      REAL(wp), DIMENSION(jpij,nlay_s)     ::   ztsb        ! Temporary temperature in the snow to check the convergence
174      REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   ztcond_i    ! Ice thermal conductivity
175      REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   zradtr_i    ! Radiation transmitted through the ice
176      REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   zradab_i    ! Radiation absorbed in the ice
177      REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   zkappa_i    ! Kappa factor in the ice
178      REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   zeta_i      ! Eta factor in the ice
179      REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   zradtr_s    ! Radiation transmited through the snow
180      REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   zradab_s    ! Radiation absorbed in the snow
181      REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   zkappa_s    ! Kappa factor in the snow
182      REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   zeta_s      ! Eta factor in the snow
183      REAL(wp), DIMENSION(jpij,nlay_i+3)   ::   zindterm    ! 'Ind'ependent term
184      REAL(wp), DIMENSION(jpij,nlay_i+3)   ::   zindtbis    ! Temporary 'ind'ependent term
185      REAL(wp), DIMENSION(jpij,nlay_i+3)   ::   zdiagbis    ! Temporary 'dia'gonal term
186      REAL(wp), DIMENSION(jpij,nlay_i+3,3) ::   ztrid       ! Tridiagonal system terms
187      REAL(wp), DIMENSION(jpij)            ::   zq_ini      ! diag errors on heat
188      REAL(wp), DIMENSION(jpij)            ::   zghe        ! G(he), th. conduct enhancement factor, mono-cat
189
190      REAL(wp) ::   zfr1, zfr2, zfrqsr_tr_i   ! dummy factor
191     
192      ! Mono-category
193      REAL(wp) :: zepsilon      ! determines thres. above which computation of G(h) is done
194      REAL(wp) :: zhe           ! dummy factor
195      REAL(wp) :: zcnd_i        ! mean sea ice thermal conductivity
196      !!------------------------------------------------------------------     
197
198      ! --- diag error on heat diffusion - PART 1 --- !
199      DO ji = 1, npti
200         zq_ini(ji) = ( SUM( e_i_1d(ji,1:nlay_i) ) * h_i_1d(ji) * r1_nlay_i +  &
201            &           SUM( e_s_1d(ji,1:nlay_s) ) * h_s_1d(ji) * r1_nlay_s ) 
202      END DO
203
204      !------------------
205      ! 1) Initialization
206      !------------------
207      DO ji = 1, npti
208         isnow(ji)= 1._wp - MAX( 0._wp , SIGN(1._wp, - h_s_1d(ji) ) )  ! is there snow or not
209         ! layer thickness
210         zh_i(ji) = h_i_1d(ji) * r1_nlay_i
211         zh_s(ji) = h_s_1d(ji) * r1_nlay_s
212      END DO
213      !
214      WHERE( zh_i(1:npti) >= epsi10 )   ;   z1_h_i(1:npti) = 1._wp / zh_i(1:npti)
215      ELSEWHERE                         ;   z1_h_i(1:npti) = 0._wp
216      END WHERE
217
218      WHERE( zh_s(1:npti) >= epsi10 )   ;   z1_h_s(1:npti) = 1._wp / zh_s(1:npti)
219      ELSEWHERE                         ;   z1_h_s(1:npti) = 0._wp
220      END WHERE
221
222      ! Store initial temperatures and non solar heat fluxes
223      IF ( k_jules == np_zdf_jules_OFF .OR. k_jules == np_zdf_jules_SND  ) THEN ! OFF or SND mode
224
225         ztsub  (1:npti)       =   t_su_1d(1:npti)                          ! surface temperature at iteration n-1
226         ztsuold(1:npti)       =   t_su_1d(1:npti)                          ! surface temperature initial value
227         zdqns_ice_b(1:npti)   =   dqns_ice_1d(1:npti)                      ! derivative of incoming nonsolar flux
228         zqns_ice_b (1:npti)   =   qns_ice_1d(1:npti)                       ! store previous qns_ice_1d value
229
230         t_su_1d(1:npti)       =   MIN( t_su_1d(1:npti), rt0 - ztsu_err )   ! required to leave the choice between melting or not
231
232      ENDIF   
233
234      ztsold (1:npti,:) = t_s_1d(1:npti,:)   ! Old snow temperature
235      ztiold (1:npti,:) = t_i_1d(1:npti,:)   ! Old ice temperature
236
237      !-------------
238      ! 2) Radiation
239      !-------------
240      ! --- Transmission/absorption of solar radiation in the ice --- !
241!     zfr1 = ( 0.18 * ( 1.0 - 0.81 ) + 0.35 * 0.81 )            !   standard value
242!     zfr2 = ( 0.82 * ( 1.0 - 0.81 ) + 0.65 * 0.81 )            !   zfr2 such that zfr1 + zfr2 to equal 1
243
244!     DO ji = 1, npti
245
246!        zfac = MAX( 0._wp , 1._wp - ( h_i_1d(ji) * 10._wp ) )     
247
248!              zfrqsr_tr_i = zfr1 + zfac * zfr2                         !   below 10 cm, linearly increase zfrqsr_tr_i until 1 at zero thickness
249!              IF ( h_s_1d(ji) >= 0.0_wp )   zfrqsr_tr_i = 0._wp     !   snow fully opaque
250
251!              qsr_ice_tr_1d(ji) = zfrqsr_tr_i * qsr_ice_1d(ji)   !   transmitted solar radiation
252
253!              zfsw(ji) = qsr_ice_1d(ji) - qsr_ice_tr_1d(ji)
254!              zftrice(ji) = qsr_ice_tr_1d(ji)
255!              i0(ji) = zfrqsr_tr_i
256
257!     END DO
258
259      zradtr_s(1:npti,0) = qsr_ice_tr_1d(1:npti)
260      DO jk = 1, nlay_s
261         DO ji = 1, npti
262            !                             ! radiation transmitted below the layer-th snow layer
263            zradtr_s(ji,jk) = zradtr_s(ji,0) * EXP( - zraext_s * zh_s(ji) * REAL(jk) )
264            !                             ! radiation absorbed by the layer-th snow layer
265            zradab_s(ji,jk) = zradtr_s(ji,jk-1) - zradtr_s(ji,jk)
266         END DO
267      END DO
268         
269      zradtr_i(1:npti,0) = zradtr_s(1:npti,nlay_s) * isnow(1:npti) + qsr_ice_tr_1d(1:npti) * ( 1._wp - isnow(1:npti) )
270      DO jk = 1, nlay_i 
271         DO ji = 1, npti
272            !                             ! radiation transmitted below the layer-th ice layer
273            zradtr_i(ji,jk) = zradtr_i(ji,0) * EXP( - rn_kappa_i * zh_i(ji) * REAL(jk) )
274            !                             ! radiation absorbed by the layer-th ice layer
275            zradab_i(ji,jk) = zradtr_i(ji,jk-1) - zradtr_i(ji,jk)
276         END DO
277      END DO
278
279      ftr_ice_1d(1:npti) = zradtr_i(1:npti,nlay_i)   ! record radiation transmitted below the ice
280      !
281      iconv    =  0          ! number of iterations
282      zdti_max =  1000._wp   ! maximal value of error on all points
283      !                                                          !============================!
284      DO WHILE ( zdti_max > zdti_bnd .AND. iconv < iconv_max )   ! Iterative procedure begins !
285         !                                                       !============================!
286         iconv = iconv + 1
287         !
288         ztib(1:npti,:) = t_i_1d(1:npti,:)
289         ztsb(1:npti,:) = t_s_1d(1:npti,:)
290         !
291         !--------------------------------
292         ! 3) Sea ice thermal conductivity
293         !--------------------------------
294         IF( ln_cndi_U64 ) THEN         !-- Untersteiner (1964) formula: k = k0 + beta.S/T
295            !
296            DO ji = 1, npti
297               ztcond_i(ji,0)      = rcdic + zbeta * sz_i_1d(ji,1)      / MIN( -epsi10, t_i_1d(ji,1) - rt0 )
298               ztcond_i(ji,nlay_i) = rcdic + zbeta * sz_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji)  - rt0 )
299            END DO
300            DO jk = 1, nlay_i-1
301               DO ji = 1, npti
302                  ztcond_i(ji,jk) = rcdic + zbeta * 0.5_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) ) /  &
303                     &                      MIN( -epsi10, 0.5_wp * (t_i_1d(ji,jk) + t_i_1d(ji,jk+1)) - rt0 )
304               END DO
305            END DO
306            !
307         ELSEIF( ln_cndi_P07 ) THEN     !-- Pringle et al formula: k = k0 + beta1.S/T - beta2.T
308            !
309            DO ji = 1, npti
310               ztcond_i(ji,0)      = rcdic + 0.09_wp * sz_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 )  &
311                  &                        - 0.011_wp * ( t_i_1d(ji,1) - rt0 )
312               ztcond_i(ji,nlay_i) = rcdic + 0.09_wp * sz_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji) - rt0 )  &
313                  &                        - 0.011_wp * ( t_bo_1d(ji) - rt0 )
314            END DO
315            DO jk = 1, nlay_i-1
316               DO ji = 1, npti
317                  ztcond_i(ji,jk) = rcdic + 0.09_wp * 0.5_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) ) /          &
318                     &                     MIN( -epsi10, 0.5_wp * (t_i_1d(ji,jk) + t_i_1d(ji,jk+1)) - rt0 )   &
319                     &                    - 0.011_wp * ( 0.5_wp * (t_i_1d(ji,jk) + t_i_1d(ji,jk+1)) - rt0 )
320               END DO
321            END DO
322            !
323         ENDIF
324         ztcond_i(1:npti,:) = MAX( zkimin, ztcond_i(1:npti,:) )       
325         !
326         !--- G(he) : enhancement of thermal conductivity in mono-category case
327         ! Computation of effective thermal conductivity G(h)
328         ! Used in mono-category case only to simulate an ITD implicitly
329         ! Fichefet and Morales Maqueda, JGR 1997
330         zghe(1:npti) = 1._wp
331         !
332         SELECT CASE ( nn_monocat )
333
334         CASE ( 1 , 3 )
335
336            zepsilon = 0.1_wp
337            DO ji = 1, npti
338               zcnd_i = SUM( ztcond_i(ji,:) ) / REAL( nlay_i+1, wp )                             ! Mean sea ice thermal conductivity
339               zhe = ( rn_cnd_s * h_i_1d(ji) + zcnd_i * h_s_1d(ji) ) / ( rn_cnd_s + zcnd_i )   ! Effective thickness he (zhe)
340               IF( zhe >=  zepsilon * 0.5_wp * EXP(1._wp) ) THEN
341                  zghe(ji) = MIN( 2._wp, 0.5_wp * ( 1._wp + LOG( 2._wp * zhe / zepsilon ) ) )    ! G(he)
342               ENDIF
343            END DO
344
345         END SELECT
346         !
347         !-----------------
348         ! 4) kappa factors
349         !-----------------
350         !--- Snow
351         DO jk = 0, nlay_s-1
352            DO ji = 1, npti
353               zkappa_s(ji,jk) = zghe(ji) * rn_cnd_s * z1_h_s(ji)
354            END DO
355         END DO
356         DO ji = 1, npti   ! Snow-ice interface
357            zfac = 0.5_wp * ( ztcond_i(ji,0) * zh_s(ji) + rn_cnd_s * zh_i(ji) )
358            IF( zfac > epsi10 ) THEN
359               zkappa_s(ji,nlay_s) = zghe(ji) * rn_cnd_s * ztcond_i(ji,0) / zfac
360            ELSE
361               zkappa_s(ji,nlay_s) = 0._wp
362            ENDIF
363         END DO
364
365         !--- Ice
366         DO jk = 0, nlay_i
367            DO ji = 1, npti
368               zkappa_i(ji,jk) = zghe(ji) * ztcond_i(ji,jk) * z1_h_i(ji)
369            END DO
370         END DO
371         DO ji = 1, npti   ! Snow-ice interface
372            zkappa_i(ji,0) = zkappa_s(ji,nlay_s) * isnow(ji) + zkappa_i(ji,0) * ( 1._wp - isnow(ji) )
373         END DO
374         !
375         !--------------------------------------
376         ! 5) Sea ice specific heat, eta factors
377         !--------------------------------------
378         DO jk = 1, nlay_i
379            DO ji = 1, npti
380               zcpi = cpic + zgamma * sz_i_1d(ji,jk) / MAX( ( t_i_1d(ji,jk) - rt0 ) * ( ztiold(ji,jk) - rt0 ), epsi10 )
381               zeta_i(ji,jk) = rdt_ice * r1_rhoic * z1_h_i(ji) / MAX( epsi10, zcpi ) 
382            END DO
383         END DO
384
385         DO jk = 1, nlay_s
386            DO ji = 1, npti
387               zeta_s(ji,jk) = rdt_ice * r1_rhosn * r1_cpic * z1_h_s(ji)
388            END DO
389         END DO
390         
391         !
392         !----------------------------------------!
393         !                                        !
394         !       JULES IF (OFF or SND MODE)       !
395         !                                        !
396         !----------------------------------------!
397         !
398         
399         IF ( k_jules == np_zdf_jules_OFF .OR. k_jules == np_zdf_jules_SND  ) THEN ! OFF or SND mode
400         
401            !
402            ! In OFF mode the original BL99 temperature computation is used
403            ! (with qsr_ice, qns_ice and dqns_ice as inputs)
404            !
405            ! In SND mode, the computation is required to compute the conduction fluxes
406            !
407         
408            !----------------------------
409            ! 6) surface flux computation
410            !----------------------------
411
412            DO ji = 1, npti
413            ! update of the non solar flux according to the update in T_su
414               qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_1d(ji) - ztsub(ji) )
415            END DO
416
417            DO ji = 1, npti
418               zfnet(ji) = qsr_ice_1d(ji) - qsr_ice_tr_1d(ji) + qns_ice_1d(ji) ! net heat flux = net solar - transmitted solar + non solar
419            END DO
420            !
421            !----------------------------
422            ! 7) tridiagonal system terms
423            !----------------------------
424            !!layer denotes the number of the layer in the snow or in the ice
425            !!jm denotes the reference number of the equation in the tridiagonal
426            !!system, terms of tridiagonal system are indexed as following :
427            !!1 is subdiagonal term, 2 is diagonal and 3 is superdiagonal one
428
429            !!ice interior terms (top equation has the same form as the others)
430            ztrid   (1:npti,:,:) = 0._wp
431            zindterm(1:npti,:)   = 0._wp
432            zindtbis(1:npti,:)   = 0._wp
433            zdiagbis(1:npti,:)   = 0._wp
434
435            DO jm = nlay_s + 2, nlay_s + nlay_i 
436            DO ji = 1, npti
437               jk = jm - nlay_s - 1
438               ztrid(ji,jm,1)  = - zeta_i(ji,jk) * zkappa_i(ji,jk-1)
439               ztrid(ji,jm,2)  = 1.0 + zeta_i(ji,jk) * ( zkappa_i(ji,jk-1) + zkappa_i(ji,jk) )
440               ztrid(ji,jm,3)  = - zeta_i(ji,jk) * zkappa_i(ji,jk)
441               zindterm(ji,jm) = ztiold(ji,jk) + zeta_i(ji,jk) * zradab_i(ji,jk)
442            END DO
443            ENDDO
444
445            jm =  nlay_s + nlay_i + 1
446            DO ji = 1, npti
447            !!ice bottom term
448            ztrid(ji,jm,1)  = - zeta_i(ji,nlay_i)*zkappa_i(ji,nlay_i-1)   
449            ztrid(ji,jm,2)  = 1.0 + zeta_i(ji,nlay_i) * ( zkappa_i(ji,nlay_i) * zg1 + zkappa_i(ji,nlay_i-1) )
450            ztrid(ji,jm,3)  = 0.0
451            zindterm(ji,jm) = ztiold(ji,nlay_i) + zeta_i(ji,nlay_i) *  &
452               &            ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i) * zg1 *  t_bo_1d(ji) ) 
453            ENDDO
454
455
456            DO ji = 1, npti
457            !                               !---------------------!
458            IF ( h_s_1d(ji) > 0.0 ) THEN    !  snow-covered cells !
459               !                            !---------------------!
460               ! snow interior terms (bottom equation has the same form as the others)
461               DO jm = 3, nlay_s + 1
462               jk = jm - 1
463               ztrid(ji,jm,1)  = - zeta_s(ji,jk) * zkappa_s(ji,jk-1)
464               ztrid(ji,jm,2)  = 1.0 + zeta_s(ji,jk) * ( zkappa_s(ji,jk-1) + zkappa_s(ji,jk) )
465               ztrid(ji,jm,3)  = - zeta_s(ji,jk)*zkappa_s(ji,jk)
466               zindterm(ji,jm) = ztsold(ji,jk) + zeta_s(ji,jk) * zradab_s(ji,jk)
467               END DO
468
469               ! case of only one layer in the ice (ice equation is altered)
470               IF ( nlay_i == 1 ) THEN
471               ztrid(ji,nlay_s+2,3)  = 0.0
472               zindterm(ji,nlay_s+2) = zindterm(ji,nlay_s+2) + zkappa_i(ji,1) * t_bo_1d(ji) 
473               ENDIF
474
475               IF ( t_su_1d(ji) < rt0 ) THEN   !--  case 1 : no surface melting
476
477               jm_min(ji) = 1
478               jm_max(ji) = nlay_i + nlay_s + 1
479
480               ! surface equation
481               ztrid(ji,1,1)  = 0.0
482               ztrid(ji,1,2)  = zdqns_ice_b(ji) - zg1s * zkappa_s(ji,0)
483               ztrid(ji,1,3)  = zg1s * zkappa_s(ji,0)
484               zindterm(ji,1) = zdqns_ice_b(ji) * t_su_1d(ji) - zfnet(ji)
485
486               ! first layer of snow equation
487               ztrid(ji,2,1)  = - zkappa_s(ji,0) * zg1s * zeta_s(ji,1)
488               ztrid(ji,2,2)  = 1.0 + zeta_s(ji,1) * ( zkappa_s(ji,1) + zkappa_s(ji,0) * zg1s )
489               ztrid(ji,2,3)  = - zeta_s(ji,1)* zkappa_s(ji,1)
490               zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) * zradab_s(ji,1)
491
492               ELSE                            !--  case 2 : surface is melting
493               !
494               jm_min(ji) = 2
495               jm_max(ji) = nlay_i + nlay_s + 1
496
497               ! first layer of snow equation
498               ztrid(ji,2,1)  = 0.0
499               ztrid(ji,2,2)  = 1.0 + zeta_s(ji,1) * ( zkappa_s(ji,1) + zkappa_s(ji,0) * zg1s )
500               ztrid(ji,2,3)  = - zeta_s(ji,1)*zkappa_s(ji,1) 
501               zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) *  &
502               &           ( zradab_s(ji,1) + zkappa_s(ji,0) * zg1s * t_su_1d(ji) ) 
503               ENDIF
504            !                               !---------------------!
505            ELSE                            ! cells without snow  !
506               !                            !---------------------!
507               !
508               IF ( t_su_1d(ji) < rt0 ) THEN   !--  case 1 : no surface melting
509               !
510               jm_min(ji) = nlay_s + 1
511               jm_max(ji) = nlay_i + nlay_s + 1
512
513               ! surface equation   
514               ztrid(ji,jm_min(ji),1)  = 0.0
515               ztrid(ji,jm_min(ji),2)  = zdqns_ice_b(ji) - zkappa_i(ji,0)*zg1   
516               ztrid(ji,jm_min(ji),3)  = zkappa_i(ji,0)*zg1
517               zindterm(ji,jm_min(ji)) = zdqns_ice_b(ji)*t_su_1d(ji) - zfnet(ji)
518
519               ! first layer of ice equation
520               ztrid(ji,jm_min(ji)+1,1)  = - zkappa_i(ji,0) * zg1 * zeta_i(ji,1)
521               ztrid(ji,jm_min(ji)+1,2)  = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 )
522               ztrid(ji,jm_min(ji)+1,3)  = - zeta_i(ji,1) * zkappa_i(ji,1) 
523               zindterm(ji,jm_min(ji)+1) = ztiold(ji,1) + zeta_i(ji,1) * zradab_i(ji,1) 
524
525               ! case of only one layer in the ice (surface & ice equations are altered)
526               IF ( nlay_i == 1 ) THEN
527               ztrid(ji,jm_min(ji),1)    = 0.0
528               ztrid(ji,jm_min(ji),2)    = zdqns_ice_b(ji) - zkappa_i(ji,0) * 2.0
529               ztrid(ji,jm_min(ji),3)    = zkappa_i(ji,0) * 2.0
530               ztrid(ji,jm_min(ji)+1,1)  = -zkappa_i(ji,0) * 2.0 * zeta_i(ji,1)
531               ztrid(ji,jm_min(ji)+1,2)  = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2.0 + zkappa_i(ji,1) )
532               ztrid(ji,jm_min(ji)+1,3)  = 0.0
533               zindterm(ji,jm_min(ji)+1) = ztiold(ji,1) + zeta_i(ji,1) *  &
534               &                      ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) )
535               ENDIF
536
537               ELSE                            !--  case 2 : surface is melting
538
539               jm_min(ji)    =  nlay_s + 2
540               jm_max(ji)    =  nlay_i + nlay_s + 1
541
542               ! first layer of ice equation
543               ztrid(ji,jm_min(ji),1)  = 0.0
544               ztrid(ji,jm_min(ji),2)  = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 ) 
545               ztrid(ji,jm_min(ji),3)  = - zeta_i(ji,1) * zkappa_i(ji,1)
546               zindterm(ji,jm_min(ji)) = ztiold(ji,1) + zeta_i(ji,1) *  &
547               &                    ( zradab_i(ji,1) + zkappa_i(ji,0) * zg1 * t_su_1d(ji) ) 
548
549               ! case of only one layer in the ice (surface & ice equations are altered)
550               IF ( nlay_i == 1 ) THEN
551               ztrid(ji,jm_min(ji),1)  = 0.0
552               ztrid(ji,jm_min(ji),2)  = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2.0 + zkappa_i(ji,1) )
553               ztrid(ji,jm_min(ji),3)  = 0.0
554               zindterm(ji,jm_min(ji)) = ztiold(ji,1) + zeta_i(ji,1) * ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) )  &
555               &                    + t_su_1d(ji) * zeta_i(ji,1) * zkappa_i(ji,0) * 2.0
556               ENDIF
557
558               ENDIF
559            ENDIF
560            !
561            zindtbis(ji,jm_min(ji)) = zindterm(ji,jm_min(ji))
562            zdiagbis(ji,jm_min(ji)) = ztrid(ji,jm_min(ji),2)
563            !
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            DO ji = 1, npti
587            ! ice temperatures
588            t_i_1d(ji,nlay_i) = zindtbis(ji,jm_max(ji)) / zdiagbis(ji,jm_max(ji))
589            END DO
590
591            DO jm = nlay_i + nlay_s, nlay_s + 2, -1
592            DO ji = 1, npti
593               jk    =  jm - nlay_s - 1
594               t_i_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,jm)
595            END DO
596            END DO
597
598            DO ji = 1, npti
599            ! snow temperatures     
600            IF( h_s_1d(ji) > 0._wp ) THEN
601               t_s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) )  &
602               &                / zdiagbis(ji,nlay_s+1)
603            ENDIF
604            ! surface temperature
605            ztsub(ji) = t_su_1d(ji)
606            IF( t_su_1d(ji) < rt0 ) THEN
607               t_su_1d(ji) = ( zindtbis(ji,jm_min(ji)) - ztrid(ji,jm_min(ji),3) *  &
608               &          ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) ) ) / zdiagbis(ji,jm_min(ji))
609            ENDIF
610            END DO
611            !
612            !--------------------------------------------------------------
613            ! 9) Has the scheme converged ?, end of the iterative procedure
614            !--------------------------------------------------------------
615            ! check that nowhere it has started to melt
616            ! zdti_max is a measure of error, it has to be under zdti_bnd
617            zdti_max = 0._wp
618            DO ji = 1, npti
619            t_su_1d(ji) = MAX( MIN( t_su_1d(ji) , rt0 ) , rt0 - 100._wp )
620            zdti_max = MAX( zdti_max, ABS( t_su_1d(ji) - ztsub(ji) ) )     
621            END DO
622
623            DO jk  =  1, nlay_s
624            DO ji = 1, npti
625               t_s_1d(ji,jk) = MAX( MIN( t_s_1d(ji,jk), rt0 ), rt0 - 100._wp )
626               zdti_max = MAX( zdti_max, ABS( t_s_1d(ji,jk) - ztsb(ji,jk) ) )
627            END DO
628            END DO
629
630            DO jk  =  1, nlay_i
631            DO ji = 1, npti
632               ztmelt_i      = -tmut * sz_i_1d(ji,jk) + rt0 
633               t_i_1d(ji,jk) =  MAX( MIN( t_i_1d(ji,jk), ztmelt_i ), rt0 - 100._wp )
634               zdti_max      =  MAX( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) )
635            END DO
636            END DO
637
638            ! Compute spatial maximum over all errors
639            ! note that this could be optimized substantially by iterating only the non-converging points
640            IF( lk_mpp ) CALL mpp_max( zdti_max, kcom=ncomm_ice )
641         !
642         !----------------------------------------!
643         !                                        !
644         !       JULES IF (RCV MODE)              !
645         !                                        !
646         !----------------------------------------!
647         !
648
649         ELSE IF ( k_jules == np_zdf_jules_RCV  ) THEN ! RCV mode
650         
651            !
652            ! In RCV mode, we use a modified BL99 solver
653            ! with conduction flux (qcn_ice) as forcing term
654            !
655            !----------------------------
656            ! 7) tridiagonal system terms
657            !----------------------------
658            !!layer denotes the number of the layer in the snow or in the ice
659            !!jm denotes the reference number of the equation in the tridiagonal
660            !!system, terms of tridiagonal system are indexed as following :
661            !!1 is subdiagonal term, 2 is diagonal and 3 is superdiagonal one
662
663            !!ice interior terms (top equation has the same form as the others)
664            ztrid   (1:npti,:,:) = 0._wp
665            zindterm(1:npti,:)   = 0._wp
666            zindtbis(1:npti,:)   = 0._wp
667            zdiagbis(1:npti,:)   = 0._wp
668
669            DO jm = nlay_s + 2, nlay_s + nlay_i 
670            DO ji = 1, npti
671               jk = jm - nlay_s - 1
672               ztrid(ji,jm,1)  = - zeta_i(ji,jk) * zkappa_i(ji,jk-1)
673               ztrid(ji,jm,2)  = 1.0 + zeta_i(ji,jk) * ( zkappa_i(ji,jk-1) + zkappa_i(ji,jk) )
674               ztrid(ji,jm,3)  = - zeta_i(ji,jk) * zkappa_i(ji,jk)
675               zindterm(ji,jm) = ztiold(ji,jk) + zeta_i(ji,jk) * zradab_i(ji,jk)
676            END DO
677            ENDDO
678
679            jm =  nlay_s + nlay_i + 1
680            DO ji = 1, npti
681            !!ice bottom term
682            ztrid(ji,jm,1)  = - zeta_i(ji,nlay_i)*zkappa_i(ji,nlay_i-1)   
683            ztrid(ji,jm,2)  = 1.0 + zeta_i(ji,nlay_i) * ( zkappa_i(ji,nlay_i) * zg1 + zkappa_i(ji,nlay_i-1) )
684            ztrid(ji,jm,3)  = 0.0
685            zindterm(ji,jm) = ztiold(ji,nlay_i) + zeta_i(ji,nlay_i) *  &
686               &            ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i) * zg1 *  t_bo_1d(ji) ) 
687            ENDDO
688
689
690            DO ji = 1, npti
691            !                              !---------------------!
692            IF ( h_s_1d(ji) > 0.0 ) THEN   !  snow-covered cells !
693               !                           !---------------------!
694               ! snow interior terms (bottom equation has the same form as the others)
695               DO jm = 3, nlay_s + 1
696               jk = jm - 1
697               ztrid(ji,jm,1)  = - zeta_s(ji,jk) * zkappa_s(ji,jk-1)
698               ztrid(ji,jm,2)  = 1.0 + zeta_s(ji,jk) * ( zkappa_s(ji,jk-1) + zkappa_s(ji,jk) )
699               ztrid(ji,jm,3)  = - zeta_s(ji,jk)*zkappa_s(ji,jk)
700               zindterm(ji,jm) = ztsold(ji,jk) + zeta_s(ji,jk) * zradab_s(ji,jk)
701               END DO
702
703               ! case of only one layer in the ice (ice equation is altered)
704               IF ( nlay_i == 1 ) THEN
705               ztrid(ji,nlay_s+2,3)  = 0.0
706               zindterm(ji,nlay_s+2) = zindterm(ji,nlay_s+2) + zkappa_i(ji,1) * t_bo_1d(ji) 
707               ENDIF
708
709               jm_min(ji) = 2
710               jm_max(ji) = nlay_i + nlay_s + 1
711
712               ! first layer of snow equation
713               ztrid(ji,2,1)  = 0.0
714               ztrid(ji,2,2)  = 1.0 + zeta_s(ji,1) * zkappa_s(ji,1)
715               ztrid(ji,2,3)  = - zeta_s(ji,1)*zkappa_s(ji,1) 
716               zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) *  &
717               &           ( zradab_s(ji,1) + qcn_ice_1d(ji) ) 
718
719            !                               !---------------------!
720            ELSE                            ! cells without snow  !
721            !                               !---------------------!
722
723               jm_min(ji)    =  nlay_s + 2
724               jm_max(ji)    =  nlay_i + nlay_s + 1
725
726               ! first layer of ice equation
727               ztrid(ji,jm_min(ji),1)  = 0.0
728               ztrid(ji,jm_min(ji),2)  = 1.0 + zeta_i(ji,1) * zkappa_i(ji,1) 
729               ztrid(ji,jm_min(ji),3)  = - zeta_i(ji,1) * zkappa_i(ji,1)
730               zindterm(ji,jm_min(ji)) = ztiold(ji,1) + zeta_i(ji,1) *  &
731               &                    ( zradab_i(ji,1) + qcn_ice_1d(ji) )
732
733               ! case of only one layer in the ice (surface & ice equations are altered)
734               IF ( nlay_i == 1 ) THEN
735               ztrid(ji,jm_min(ji),1)  = 0.0
736               ztrid(ji,jm_min(ji),2)  = 1.0 + zeta_i(ji,1) * zkappa_i(ji,1)
737               ztrid(ji,jm_min(ji),3)  = 0.0
738               zindterm(ji,jm_min(ji)) = ztiold(ji,1) + zeta_i(ji,1) * ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji)    &
739               &                    + qcn_ice_1d(ji) )
740
741               ENDIF
742
743            ENDIF
744            !
745            zindtbis(ji,jm_min(ji)) = zindterm(ji,jm_min(ji))
746            zdiagbis(ji,jm_min(ji)) = ztrid(ji,jm_min(ji),2)
747            !
748            END DO
749            !
750            !------------------------------
751            ! 8) tridiagonal system solving
752            !------------------------------
753            ! Solve the tridiagonal system with Gauss elimination method.
754            ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, McGraw-Hill 1984
755            jm_maxt = 0
756            jm_mint = nlay_i+5
757            DO ji = 1, npti
758            jm_mint = MIN(jm_min(ji),jm_mint)
759            jm_maxt = MAX(jm_max(ji),jm_maxt)
760            END DO
761
762            DO jk = jm_mint+1, jm_maxt
763            DO ji = 1, npti
764               jm = min(max(jm_min(ji)+1,jk),jm_max(ji))
765               zdiagbis(ji,jm) = ztrid(ji,jm,2)  - ztrid(ji,jm,1) * ztrid(ji,jm-1,3)  / zdiagbis(ji,jm-1)
766               zindtbis(ji,jm) = zindterm(ji,jm) - ztrid(ji,jm,1) * zindtbis(ji,jm-1) / zdiagbis(ji,jm-1)
767            END DO
768            END DO
769
770            DO ji = 1, npti
771            ! ice temperatures
772            t_i_1d(ji,nlay_i) = zindtbis(ji,jm_max(ji)) / zdiagbis(ji,jm_max(ji))
773            END DO
774
775            DO jm = nlay_i + nlay_s, nlay_s + 2, -1
776            DO ji = 1, npti
777               jk    =  jm - nlay_s - 1
778               t_i_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,jm)
779            END DO
780            END DO
781
782            DO ji = 1, npti
783            ! snow temperatures     
784            IF( h_s_1d(ji) > 0._wp ) THEN
785               t_s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) )  &
786               &                / zdiagbis(ji,nlay_s+1)
787            ENDIF
788            END DO
789            !
790            !--------------------------------------------------------------
791            ! 9) Has the scheme converged ?, end of the iterative procedure
792            !--------------------------------------------------------------
793            ! check that nowhere it has started to melt
794            ! zdti_max is a measure of error, it has to be under zdti_bnd
795            zdti_max = 0._wp
796
797            DO jk  =  1, nlay_s
798            DO ji = 1, npti
799               t_s_1d(ji,jk) = MAX( MIN( t_s_1d(ji,jk), rt0 ), rt0 - 100._wp )
800               zdti_max = MAX( zdti_max, ABS( t_s_1d(ji,jk) - ztsb(ji,jk) ) )
801            END DO
802            END DO
803
804            DO jk  =  1, nlay_i
805            DO ji = 1, npti
806               ztmelt_i      = -tmut * sz_i_1d(ji,jk) + rt0 
807               t_i_1d(ji,jk) =  MAX( MIN( t_i_1d(ji,jk), ztmelt_i ), rt0 - 100._wp )
808               zdti_max      =  MAX( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) )
809            END DO
810            END DO
811
812            ! Compute spatial maximum over all errors
813            ! note that this could be optimized substantially by iterating only the non-converging points
814            IF( lk_mpp ) CALL mpp_max( zdti_max, kcom=ncomm_ice )
815
816         ENDIF ! k_jules
817         
818      END DO  ! End of the do while iterative procedure
819     
820      IF( ln_icectl .AND. lwp ) THEN
821         WRITE(numout,*) ' zdti_max : ', zdti_max
822         WRITE(numout,*) ' iconv    : ', iconv
823      ENDIF
824     
825      !
826      !-----------------------------
827      ! 10) Fluxes at the interfaces
828      !-----------------------------
829      !
830      ! --- update conduction fluxes
831      !
832      DO ji = 1, npti
833      !                                ! surface ice conduction flux
834         fc_su(ji)   =  -           isnow(ji)   * zkappa_s(ji,0) * zg1s * (t_s_1d(ji,1) - t_su_1d(ji))   &
835            &           - ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1  * (t_i_1d(ji,1) - t_su_1d(ji))
836      !                                ! bottom ice conduction flux
837         fc_bo_i(ji) =  - zkappa_i(ji,nlay_i) * ( zg1*(t_bo_1d(ji) - t_i_1d(ji,nlay_i)) )
838      END DO
839     
840      !
841      ! --- Diagnose the heat loss due to changing non-solar / conduction flux --- !
842      !
843      DO ji = 1, npti
844            IF ( k_jules == np_zdf_jules_OFF .OR. k_jules == np_zdf_jules_SND  ) THEN 
845                   ! OFF or SND mode
846               hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) - ( qns_ice_1d(ji)  - zqns_ice_b(ji) ) * a_i_1d(ji) 
847            ELSE   ! RCV mode
848               hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) - ( fc_su(ji) - qcn_ice_1d(ji) ) * a_i_1d(ji) 
849            ENDIF
850      END DO
851     
852      !
853      ! --- Diagnose the heat loss due to non-fully converged temperature solution (should not be above 10-4 W-m2) --- !
854      !
855
856      IF ( ( k_jules == np_zdf_jules_OFF ) .OR. ( k_jules == np_zdf_jules_RCV ) ) THEN ! OFF
857     
858         CALL ice_thd_enmelt       
859
860         !     zhfx_err = correction on the diagnosed heat flux due to non-convergence of the algorithm used to solve heat equation
861         DO ji = 1, npti
862            zdq = - zq_ini(ji) + ( SUM( e_i_1d(ji,1:nlay_i) ) * h_i_1d(ji) * r1_nlay_i +  &
863               &                   SUM( e_s_1d(ji,1:nlay_s) ) * h_s_1d(ji) * r1_nlay_s )
864               
865            IF ( ( k_jules == np_zdf_jules_OFF ) ) THEN
866               
867                  IF( t_su_1d(ji) < rt0 ) THEN  ! case T_su < 0degC
868                     zhfx_err = ( qns_ice_1d(ji) + qsr_ice_1d(ji)    - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq * r1_rdtice ) * a_i_1d(ji) 
869                  ELSE                          ! case T_su = 0degC
870                     zhfx_err = ( fc_su(ji)      + qsr_ice_tr_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq * r1_rdtice ) * a_i_1d(ji)
871                  ENDIF
872           
873            ELSE ! RCV CASE
874           
875               zhfx_err = ( fc_su(ji) + qsr_ice_tr_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq * r1_rdtice ) * a_i_1d(ji)
876           
877            ENDIF
878           
879            ! total heat sink to be sent to the ocean
880            hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) + zhfx_err
881
882            ! hfx_dif = Heat flux diagnostic of sensible heat used to warm/cool ice in W.m-2   
883            hfx_dif_1d(ji) = hfx_dif_1d(ji) - zdq * r1_rdtice * a_i_1d(ji)
884   
885         END DO 
886           
887         !
888         ! --- SIMIP diagnostics
889         !
890
891         DO ji = 1, npti
892            !--- Conduction fluxes (positive downwards)
893            diag_fc_bo_1d(ji) = diag_fc_bo_1d(ji) + fc_bo_i(ji) * a_i_1d(ji) / at_i_1d(ji)
894            diag_fc_su_1d(ji) = diag_fc_su_1d(ji) + fc_su(ji)   * a_i_1d(ji) / at_i_1d(ji)
895   
896            !--- Snow-ice interfacial temperature (diagnostic SIMIP)
897            zfac = rn_cnd_s * zh_i(ji) + ztcond_i(ji,1) * zh_s(ji)
898            IF( zh_s(ji) >= 1.e-3 .AND. zfac > epsi10 ) THEN
899               t_si_1d(ji) = ( rn_cnd_s       * zh_i(ji) * t_s_1d(ji,1) +   &
900                  &            ztcond_i(ji,1) * zh_s(ji) * t_i_1d(ji,1) ) / zfac
901            ELSE
902               t_si_1d(ji) = t_su_1d(ji)
903            ENDIF
904         END DO
905     
906      ENDIF
907      !
908      !---------------------------------------------------------------------------------------
909      ! 11) Jules coupling: reset inner snow and ice temperatures, update conduction fluxes
910      !---------------------------------------------------------------------------------------
911      !
912      IF ( k_jules == np_zdf_jules_SND ) THEN   ! --- Jules coupling in "SND" mode
913     
914         ! Restore temperatures to their initial values
915         t_s_1d(1:npti,:)        = ztsold (1:npti,:)
916         t_i_1d(1:npti,:)        = ztiold (1:npti,:)
917         qcn_ice_1d(1:npti)      = fc_su(1:npti)
918         
919      ENDIF
920     
921   END SUBROUTINE ice_thd_zdf_BL99
922
923
924
925   SUBROUTINE ice_thd_enmelt
926      !!-------------------------------------------------------------------
927      !!                   ***  ROUTINE ice_thd_enmelt ***
928      !!                 
929      !! ** Purpose :   Computes sea ice energy of melting q_i (J.m-3) from temperature
930      !!
931      !! ** Method  :   Formula (Bitz and Lipscomb, 1999)
932      !!-------------------------------------------------------------------
933      INTEGER  ::   ji, jk   ! dummy loop indices
934      REAL(wp) ::   ztmelts  ! local scalar
935      !!-------------------------------------------------------------------
936      !
937      DO jk = 1, nlay_i             ! Sea ice energy of melting
938         DO ji = 1, npti
939            ztmelts      = - tmut  * sz_i_1d(ji,jk)
940            t_i_1d(ji,jk) = MIN( t_i_1d(ji,jk), ztmelts + rt0 ) ! Force t_i_1d to be lower than melting point
941                                                                !   (sometimes dif scheme produces abnormally high temperatures)   
942            e_i_1d(ji,jk) = rhoic * ( cpic * ( ztmelts - ( t_i_1d(ji,jk) - rt0 ) )           &
943               &                    + lfus * ( 1._wp - ztmelts / ( t_i_1d(ji,jk) - rt0 ) )   &
944               &                    - rcp  *   ztmelts )
945         END DO
946      END DO
947      DO jk = 1, nlay_s             ! Snow energy of melting
948         DO ji = 1, npti
949            e_s_1d(ji,jk) = rhosn * ( cpic * ( rt0 - t_s_1d(ji,jk) ) + lfus )
950         END DO
951      END DO
952      !
953   END SUBROUTINE ice_thd_enmelt
954
955
956
957   SUBROUTINE ice_thd_zdf_init
958      !!-----------------------------------------------------------------------
959      !!                   ***  ROUTINE ice_thd_zdf_init ***
960      !!                 
961      !! ** Purpose :   Physical constants and parameters associated with
962      !!                ice thermodynamics
963      !!
964      !! ** Method  :   Read the namthd_zdf namelist and check the parameters
965      !!                called at the first timestep (nit000)
966      !!
967      !! ** input   :   Namelist namthd_zdf
968      !!-------------------------------------------------------------------
969      INTEGER  ::   ios, ioptio   ! Local integer output status for namelist read
970      !!
971      NAMELIST/namthd_zdf/ ln_zdf_BL99, ln_cndi_U64, ln_cndi_P07, rn_cnd_s, rn_kappa_i
972      !!-------------------------------------------------------------------
973      !
974      REWIND( numnam_ice_ref )              ! Namelist namthd_zdf in reference namelist : Ice thermodynamics
975      READ  ( numnam_ice_ref, namthd_zdf, IOSTAT = ios, ERR = 901)
976901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namthd_zdf in reference namelist', lwp )
977
978      REWIND( numnam_ice_cfg )              ! Namelist namthd_zdf in configuration namelist : Ice thermodynamics
979      READ  ( numnam_ice_cfg, namthd_zdf, IOSTAT = ios, ERR = 902 )
980902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namthd_zdf in configuration namelist', lwp )
981      IF(lwm) WRITE ( numoni, namthd_zdf )
982      !
983      !
984      IF(lwp) THEN                          ! control print
985         WRITE(numout,*) 'ice_thd_zdf_init: Ice vertical heat diffusion'
986         WRITE(numout,*) '~~~~~~~~~~~~~~~~'
987         WRITE(numout,*) '   Namelist namthd_zdf:'
988         WRITE(numout,*) '      Bitz and Lipscomb (1999) formulation                    ln_zdf_BL99  = ', ln_zdf_BL99
989         WRITE(numout,*) '      thermal conductivity in the ice (Untersteiner 1964)     ln_cndi_U64  = ', ln_cndi_U64
990         WRITE(numout,*) '      thermal conductivity in the ice (Pringle et al 2007)    ln_cndi_P07  = ', ln_cndi_P07
991         WRITE(numout,*) '      thermal conductivity in the snow                        rn_cnd_s     = ', rn_cnd_s
992         WRITE(numout,*) '      extinction radiation parameter in sea ice               rn_kappa_i   = ', rn_kappa_i
993     ENDIF
994
995      !
996      IF ( ( ln_cndi_U64 .AND. ln_cndi_P07 ) .OR. ( .NOT.ln_cndi_U64 .AND. .NOT.ln_cndi_P07 ) ) THEN
997         CALL ctl_stop( 'ice_thd_zdf_init: choose one and only one formulation for thermal conduction (ln_cndi_U64 or ln_cndi_P07)' )
998      ENDIF
999
1000      !                             !== set the choice of ice vertical thermodynamic formulation ==!
1001      ioptio = 0 
1002      !      !--- BL99 thermo dynamics                               (linear liquidus + constant thermal properties)
1003      IF( ln_zdf_BL99 ) THEN   ;   ioptio = ioptio + 1   ;   nice_zdf = np_BL99   ;   ENDIF
1004      IF( ioptio /= 1 )    CALL ctl_stop( 'ice_thd_init: one and only one ice thermo option has to be defined ' )
1005      !
1006   END SUBROUTINE ice_thd_zdf_init
1007
1008#else
1009   !!----------------------------------------------------------------------
1010   !!   Default option       Dummy Module             No ESIM sea-ice model
1011   !!---------------------------------------------------------------------
1012#endif
1013
1014   !!======================================================================
1015END MODULE icethd_zdf
Note: See TracBrowser for help on using the repository browser.