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

source: branches/UKMO/dev_r5785_SSS_obsoper/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90 @ 7773

Last change on this file since 7773 was 7773, checked in by mattmartin, 7 years ago

Committing updates after doing the following:

  • merging the branch dev_r4650_general_vert_coord_obsoper@7763 into this branch
  • updating it so that the following OBS changes were implemented correctly on top of the simplification changes:
    • generalised vertical coordinate for profile obs. This was done so that is now the default option.
    • sst bias correction implemented with the new simplified obs code.
    • included the biogeochemical obs types int he new simplified obs code.
    • included the changes to exclude obs in the boundary for limited area models
    • included other changes for the efficiency of the obs operator to remove global arrays.
File size: 42.4 KB
RevLine 
[825]1MODULE limthd_dif
2   !!======================================================================
3   !!                       ***  MODULE limthd_dif ***
4   !!                       heat diffusion in sea ice
5   !!                   computation of surface and inner T 
6   !!======================================================================
[2715]7   !! History :  LIM  ! 02-2003 (M. Vancoppenolle) original 1D code
8   !!                 ! 06-2005 (M. Vancoppenolle) 3d version
9   !!                 ! 11-2006 (X Fettweis) Vectorization by Xavier
10   !!                 ! 04-2007 (M. Vancoppenolle) Energy conservation
11   !!            4.0  ! 2011-02 (G. Madec) dynamical allocation
[4161]12   !!             -   ! 2012-05 (C. Rousset) add penetration solar flux
[825]13   !!----------------------------------------------------------------------
[2528]14#if defined key_lim3
15   !!----------------------------------------------------------------------
16   !!   'key_lim3'                                      LIM3 sea-ice model
17   !!----------------------------------------------------------------------
[3625]18   USE par_oce        ! ocean parameters
19   USE phycst         ! physical constants (ocean directory)
20   USE ice            ! LIM-3 variables
21   USE thd_ice        ! LIM-3: thermodynamics
22   USE in_out_manager ! I/O manager
23   USE lib_mpp        ! MPP library
24   USE wrk_nemo       ! work arrays
25   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
[921]26
[825]27   IMPLICIT NONE
28   PRIVATE
29
[2528]30   PUBLIC   lim_thd_dif   ! called by lim_thd
[825]31
32   !!----------------------------------------------------------------------
[4161]33   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
[1156]34   !! $Id$
[2591]35   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
[825]36   !!----------------------------------------------------------------------
37CONTAINS
38
[4688]39   SUBROUTINE lim_thd_dif( kideb , kiut )
[921]40      !!------------------------------------------------------------------
41      !!                ***  ROUTINE lim_thd_dif  ***
42      !! ** Purpose :
43      !!           This routine determines the time evolution of snow and sea-ice
44      !!           temperature profiles.
45      !! ** Method  :
46      !!           This is done by solving the heat equation diffusion with
47      !!           a Neumann boundary condition at the surface and a Dirichlet one
48      !!           at the bottom. Solar radiation is partially absorbed into the ice.
49      !!           The specific heat and thermal conductivities depend on ice salinity
50      !!           and temperature to take into account brine pocket melting. The
51      !!           numerical
52      !!           scheme is an iterative Crank-Nicolson on a non-uniform multilayer grid
53      !!           in the ice and snow system.
54      !!
55      !!           The successive steps of this routine are
56      !!           1.  Thermal conductivity at the interfaces of the ice layers
57      !!           2.  Internal absorbed radiation
58      !!           3.  Scale factors due to non-uniform grid
59      !!           4.  Kappa factors
60      !!           Then iterative procedure begins
61      !!           5.  specific heat in the ice
62      !!           6.  eta factors
63      !!           7.  surface flux computation
64      !!           8.  tridiagonal system terms
65      !!           9.  solving the tridiagonal system with Gauss elimination
66      !!           Iterative procedure ends according to a criterion on evolution
67      !!           of temperature
68      !!
69      !! ** Arguments :
70      !!           kideb , kiut : Starting and ending points on which the
71      !!                         the computation is applied
72      !!
73      !! ** Inputs / Ouputs : (global commons)
[4872]74      !!           surface temperature : t_su_1d
75      !!           ice/snow temperatures   : t_i_1d, t_s_1d
76      !!           ice salinities          : s_i_1d
[921]77      !!           number of layers in the ice/snow: nlay_i, nlay_s
78      !!           profile of the ice/snow layers : z_i, z_s
[4872]79      !!           total ice/snow thickness : ht_i_1d, ht_s_1d
[921]80      !!
81      !! ** External :
82      !!
83      !! ** References :
84      !!
85      !! ** History :
86      !!           (02-2003) Martin Vancoppenolle, Louvain-la-Neuve, Belgium
87      !!           (06-2005) Martin Vancoppenolle, 3d version
88      !!           (11-2006) Vectorized by Xavier Fettweis (UCL-ASTR)
89      !!           (04-2007) Energy conservation tested by M. Vancoppenolle
90      !!------------------------------------------------------------------
[4688]91      INTEGER , INTENT(in) ::   kideb, kiut   ! Start/End point on which the  the computation is applied
[825]92
[921]93      !! * Local variables
[3294]94      INTEGER ::   ji          ! spatial loop index
95      INTEGER ::   ii, ij      ! temporary dummy loop index
96      INTEGER ::   numeq       ! current reference number of equation
[4870]97      INTEGER ::   jk       ! vertical dummy loop index
[3294]98      INTEGER ::   nconv       ! number of iterations in iterative procedure
99      INTEGER ::   minnumeqmin, maxnumeqmax
[5682]100     
[4688]101      INTEGER, POINTER, DIMENSION(:) ::   numeqmin   ! reference number of top equation
102      INTEGER, POINTER, DIMENSION(:) ::   numeqmax   ! reference number of bottom equation
[5682]103     
[3294]104      REAL(wp) ::   zg1s      =  2._wp        ! for the tridiagonal system
105      REAL(wp) ::   zg1       =  2._wp        !
106      REAL(wp) ::   zgamma    =  18009._wp    ! for specific heat
107      REAL(wp) ::   zbeta     =  0.117_wp     ! for thermal conductivity (could be 0.13)
[4990]108      REAL(wp) ::   zraext_s  =  10._wp       ! extinction coefficient of radiation in the snow
[3294]109      REAL(wp) ::   zkimin    =  0.10_wp      ! minimum ice thermal conductivity
[4688]110      REAL(wp) ::   ztsu_err  =  1.e-5_wp     ! range around which t_su is considered as 0°C
[2715]111      REAL(wp) ::   ztmelt_i    ! ice melting temperature
112      REAL(wp) ::   zerritmax   ! current maximal error on temperature
[5682]113      REAL(wp) ::   zhsu
114     
115      REAL(wp), POINTER, DIMENSION(:)     ::   isnow       ! switch for presence (1) or absence (0) of snow
116      REAL(wp), POINTER, DIMENSION(:)     ::   ztsub       ! old surface temperature (before the iterative procedure )
117      REAL(wp), POINTER, DIMENSION(:)     ::   ztsubit     ! surface temperature at previous iteration
118      REAL(wp), POINTER, DIMENSION(:)     ::   zh_i        ! ice layer thickness
119      REAL(wp), POINTER, DIMENSION(:)     ::   zh_s        ! snow layer thickness
120      REAL(wp), POINTER, DIMENSION(:)     ::   zfsw        ! solar radiation absorbed at the surface
121      REAL(wp), POINTER, DIMENSION(:)     ::   zqns_ice_b  ! solar radiation absorbed at the surface
122      REAL(wp), POINTER, DIMENSION(:)     ::   zf          ! surface flux function
123      REAL(wp), POINTER, DIMENSION(:)     ::   dzf         ! derivative of the surface flux function
124      REAL(wp), POINTER, DIMENSION(:)     ::   zerrit      ! current error on temperature
125      REAL(wp), POINTER, DIMENSION(:)     ::   zdifcase    ! case of the equation resolution (1->4)
126      REAL(wp), POINTER, DIMENSION(:)     ::   zftrice     ! solar radiation transmitted through the ice
127      REAL(wp), POINTER, DIMENSION(:)     ::   zihic
128     
129      REAL(wp), POINTER, DIMENSION(:,:)   ::   ztcond_i    ! Ice thermal conductivity
130      REAL(wp), POINTER, DIMENSION(:,:)   ::   zradtr_i    ! Radiation transmitted through the ice
131      REAL(wp), POINTER, DIMENSION(:,:)   ::   zradab_i    ! Radiation absorbed in the ice
132      REAL(wp), POINTER, DIMENSION(:,:)   ::   zkappa_i    ! Kappa factor in the ice
133      REAL(wp), POINTER, DIMENSION(:,:)   ::   ztib        ! Old temperature in the ice
134      REAL(wp), POINTER, DIMENSION(:,:)   ::   zeta_i      ! Eta factor in the ice
135      REAL(wp), POINTER, DIMENSION(:,:)   ::   ztitemp     ! Temporary temperature in the ice to check the convergence
136      REAL(wp), POINTER, DIMENSION(:,:)   ::   zspeche_i   ! Ice specific heat
137      REAL(wp), POINTER, DIMENSION(:,:)   ::   z_i         ! Vertical cotes of the layers in the ice
138      REAL(wp), POINTER, DIMENSION(:,:)   ::   zradtr_s    ! Radiation transmited through the snow
139      REAL(wp), POINTER, DIMENSION(:,:)   ::   zradab_s    ! Radiation absorbed in the snow
140      REAL(wp), POINTER, DIMENSION(:,:)   ::   zkappa_s    ! Kappa factor in the snow
141      REAL(wp), POINTER, DIMENSION(:,:)   ::   zeta_s      ! Eta factor in the snow
142      REAL(wp), POINTER, DIMENSION(:,:)   ::   ztstemp     ! Temporary temperature in the snow to check the convergence
143      REAL(wp), POINTER, DIMENSION(:,:)   ::   ztsb        ! Temporary temperature in the snow
144      REAL(wp), POINTER, DIMENSION(:,:)   ::   z_s         ! Vertical cotes of the layers in the snow
145      REAL(wp), POINTER, DIMENSION(:,:)   ::   zindterm    ! 'Ind'ependent term
146      REAL(wp), POINTER, DIMENSION(:,:)   ::   zindtbis    ! Temporary 'ind'ependent term
147      REAL(wp), POINTER, DIMENSION(:,:)   ::   zdiagbis    ! Temporary 'dia'gonal term
148      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ztrid       ! Tridiagonal system terms
149     
[4688]150      ! diag errors on heat
[5682]151      REAL(wp), POINTER, DIMENSION(:)     :: zdq, zq_ini, zhfx_err
152     
153      ! Mono-category
154      REAL(wp)                            :: zepsilon      ! determines thres. above which computation of G(h) is done
155      REAL(wp)                            :: zratio_s      ! dummy factor
156      REAL(wp)                            :: zratio_i      ! dummy factor
157      REAL(wp)                            :: zh_thres      ! thickness thres. for G(h) computation
158      REAL(wp)                            :: zhe           ! dummy factor
159      REAL(wp)                            :: zkimean       ! mean sea ice thermal conductivity
160      REAL(wp)                            :: zfac          ! dummy factor
161      REAL(wp)                            :: zihe          ! dummy factor
162      REAL(wp)                            :: zheshth       ! dummy factor
163     
164      REAL(wp), POINTER, DIMENSION(:)     :: zghe          ! G(he), th. conduct enhancement factor, mono-cat
165     
[3625]166      !!------------------------------------------------------------------     
[3610]167      !
[5682]168      CALL wrk_alloc( jpij, numeqmin, numeqmax )
169      CALL wrk_alloc( jpij, isnow, ztsub, ztsubit, zh_i, zh_s, zfsw )
170      CALL wrk_alloc( jpij, zf, dzf, zqns_ice_b, zerrit, zdifcase, zftrice, zihic, zghe )
171      CALL wrk_alloc( jpij,nlay_i+1, ztcond_i, zradtr_i, zradab_i, zkappa_i, ztib, zeta_i, ztitemp, z_i, zspeche_i, kjstart=0 )
172      CALL wrk_alloc( jpij,nlay_s+1,           zradtr_s, zradab_s, zkappa_s, ztsb, zeta_s, ztstemp, z_s, kjstart=0 )
173      CALL wrk_alloc( jpij,nlay_i+3, zindterm, zindtbis, zdiagbis  )
174      CALL wrk_alloc( jpij,nlay_i+3,3, ztrid )
[4688]175
[4990]176      CALL wrk_alloc( jpij, zdq, zq_ini, zhfx_err )
[4688]177
178      ! --- diag error on heat diffusion - PART 1 --- !
179      zdq(:) = 0._wp ; zq_ini(:) = 0._wp     
180      DO ji = kideb, kiut
[5682]181         zq_ini(ji) = ( SUM( q_i_1d(ji,1:nlay_i) ) * ht_i_1d(ji) * r1_nlay_i +  &
182            &           SUM( q_s_1d(ji,1:nlay_s) ) * ht_s_1d(ji) * r1_nlay_s ) 
[4688]183      END DO
184
[921]185      !------------------------------------------------------------------------------!
186      ! 1) Initialization                                                            !
187      !------------------------------------------------------------------------------!
188      DO ji = kideb , kiut
[5682]189         isnow(ji)= 1._wp - MAX( 0._wp , SIGN(1._wp, - ht_s_1d(ji) ) )  ! is there snow or not
[921]190         ! layer thickness
[5682]191         zh_i(ji) = ht_i_1d(ji) * r1_nlay_i
192         zh_s(ji) = ht_s_1d(ji) * r1_nlay_s
[921]193      END DO
[825]194
[921]195      !--------------------
196      ! Ice / snow layers
197      !--------------------
[825]198
[2715]199      z_s(:,0) = 0._wp   ! vert. coord. of the up. lim. of the 1st snow layer
200      z_i(:,0) = 0._wp   ! vert. coord. of the up. lim. of the 1st ice layer
[825]201
[4870]202      DO jk = 1, nlay_s            ! vert. coord of the up. lim. of the layer-th snow layer
[921]203         DO ji = kideb , kiut
[5682]204            z_s(ji,jk) = z_s(ji,jk-1) + ht_s_1d(ji) * r1_nlay_s
[921]205         END DO
206      END DO
[825]207
[4870]208      DO jk = 1, nlay_i            ! vert. coord of the up. lim. of the layer-th ice layer
[921]209         DO ji = kideb , kiut
[5682]210            z_i(ji,jk) = z_i(ji,jk-1) + ht_i_1d(ji) * r1_nlay_i
[921]211         END DO
212      END DO
213      !
214      !------------------------------------------------------------------------------|
[5682]215      ! 2) Radiation                                                       |
[921]216      !------------------------------------------------------------------------------|
217      !
218      !-------------------
219      ! Computation of i0
220      !-------------------
221      ! i0 describes the fraction of solar radiation which does not contribute
222      ! to the surface energy budget but rather penetrates inside the ice.
223      ! We assume that no radiation is transmitted through the snow
224      ! If there is no no snow
225      ! zfsw    = (1-i0).qsr_ice   is absorbed at the surface
226      ! zftrice = io.qsr_ice       is below the surface
[4688]227      ! ftr_ice = io.qsr_ice.exp(-k(h_i)) transmitted below the ice
[5682]228      ! fr1_i0_1d = i0 for a thin ice cover, fr1_i0_2d = i0 for a thick ice cover
229      zhsu = 0.1_wp ! threshold for the computation of i0
[921]230      DO ji = kideb , kiut
231         ! switches
[5682]232         isnow(ji) = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_s_1d(ji) ) ) 
[921]233         ! hs > 0, isnow = 1
[5682]234         zihic(ji) = MAX( 0._wp , 1._wp - ( ht_i_1d(ji) / zhsu ) )     
[825]235
[5682]236         i0(ji)    = ( 1._wp - isnow(ji) ) * ( fr1_i0_1d(ji) + zihic(ji) * fr2_i0_1d(ji) )
[921]237      END DO
[825]238
[921]239      !-------------------------------------------------------
240      ! Solar radiation absorbed / transmitted at the surface
241      ! Derivative of the non solar flux
242      !-------------------------------------------------------
243      DO ji = kideb , kiut
[5682]244         zfsw   (ji)    =  qsr_ice_1d(ji) * ( 1 - i0(ji) )   ! Shortwave radiation absorbed at surface
245         zftrice(ji)    =  qsr_ice_1d(ji) *       i0(ji)     ! Solar radiation transmitted below the surface layer
246         dzf    (ji)    = dqns_ice_1d(ji)                    ! derivative of incoming nonsolar flux
247         zqns_ice_b(ji) = qns_ice_1d(ji)                     ! store previous qns_ice_1d value
[921]248      END DO
[825]249
[921]250      !---------------------------------------------------------
251      ! Transmission - absorption of solar radiation in the ice
252      !---------------------------------------------------------
[825]253
[2715]254      DO ji = kideb, kiut           ! snow initialization
255         zradtr_s(ji,0) = zftrice(ji)     ! radiation penetrating through snow
[921]256      END DO
[825]257
[4870]258      DO jk = 1, nlay_s          ! Radiation through snow
[2715]259         DO ji = kideb, kiut
260            !                             ! radiation transmitted below the layer-th snow layer
[4870]261            zradtr_s(ji,jk) = zradtr_s(ji,0) * EXP( - zraext_s * ( MAX ( 0._wp , z_s(ji,jk) ) ) )
[2715]262            !                             ! radiation absorbed by the layer-th snow layer
[4870]263            zradab_s(ji,jk) = zradtr_s(ji,jk-1) - zradtr_s(ji,jk)
[921]264         END DO
265      END DO
[825]266
[2715]267      DO ji = kideb, kiut           ! ice initialization
[5682]268         zradtr_i(ji,0) = zradtr_s(ji,nlay_s) * isnow(ji) + zftrice(ji) * ( 1._wp - isnow(ji) )
[921]269      END DO
[825]270
[4870]271      DO jk = 1, nlay_i          ! Radiation through ice
[2715]272         DO ji = kideb, kiut
273            !                             ! radiation transmitted below the layer-th ice layer
[5682]274            zradtr_i(ji,jk) = zradtr_i(ji,0) * EXP( - rn_kappa_i * ( MAX ( 0._wp , z_i(ji,jk) ) ) )
[2715]275            !                             ! radiation absorbed by the layer-th ice layer
[4870]276            zradab_i(ji,jk) = zradtr_i(ji,jk-1) - zradtr_i(ji,jk)
[921]277         END DO
278      END DO
[825]279
[2715]280      DO ji = kideb, kiut           ! Radiation transmitted below the ice
[4688]281         ftr_ice_1d(ji) = zradtr_i(ji,nlay_i) 
[921]282      END DO
[834]283
[921]284      !------------------------------------------------------------------------------|
285      !  3) Iterative procedure begins                                               |
286      !------------------------------------------------------------------------------|
287      !
[2715]288      DO ji = kideb, kiut        ! Old surface temperature
[4872]289         ztsub  (ji) =  t_su_1d(ji)                              ! temperature at the beg of iter pr.
290         ztsubit(ji) =  t_su_1d(ji)                              ! temperature at the previous iter
[5682]291         t_su_1d(ji) =  MIN( t_su_1d(ji), rt0 - ztsu_err )       ! necessary
292         zerrit (ji) =  1000._wp                                 ! initial value of error
[921]293      END DO
[825]294
[4870]295      DO jk = 1, nlay_s       ! Old snow temperature
[921]296         DO ji = kideb , kiut
[4872]297            ztsb(ji,jk) =  t_s_1d(ji,jk)
[921]298         END DO
299      END DO
[825]300
[4870]301      DO jk = 1, nlay_i       ! Old ice temperature
[921]302         DO ji = kideb , kiut
[4872]303            ztib(ji,jk) =  t_i_1d(ji,jk)
[921]304         END DO
305      END DO
[825]306
[2715]307      nconv     =  0           ! number of iterations
308      zerritmax =  1000._wp    ! maximal value of error on all points
[825]309
[5682]310      DO WHILE ( zerritmax > rn_terr_dif .AND. nconv < nn_conv_dif )
[921]311         !
[2715]312         nconv = nconv + 1
313         !
[921]314         !------------------------------------------------------------------------------|
315         ! 4) Sea ice thermal conductivity                                              |
316         !------------------------------------------------------------------------------|
317         !
[5682]318         IF( nn_ice_thcon == 0 ) THEN      ! Untersteiner (1964) formula
[921]319            DO ji = kideb , kiut
[5682]320               ztcond_i(ji,0) = rcdic + zbeta * s_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 )
321               ztcond_i(ji,0) = MAX( ztcond_i(ji,0), zkimin )
[921]322            END DO
[4870]323            DO jk = 1, nlay_i-1
[921]324               DO ji = kideb , kiut
[5682]325                  ztcond_i(ji,jk) = rcdic + zbeta * ( s_i_1d(ji,jk) + s_i_1d(ji,jk+1) ) /  &
326                     MIN(-2.0_wp * epsi10, t_i_1d(ji,jk) + t_i_1d(ji,jk+1) - 2.0_wp * rt0)
327                  ztcond_i(ji,jk) = MAX( ztcond_i(ji,jk), zkimin )
[921]328               END DO
329            END DO
330         ENDIF
[825]331
[5682]332         IF( nn_ice_thcon == 1 ) THEN      ! Pringle et al formula included: 2.11 + 0.09 S/T - 0.011.T
[921]333            DO ji = kideb , kiut
[5682]334               ztcond_i(ji,0) = rcdic + 0.090_wp * s_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 )   &
335                  &                   - 0.011_wp * ( t_i_1d(ji,1) - rt0 ) 
[2715]336               ztcond_i(ji,0) = MAX( ztcond_i(ji,0), zkimin )
[921]337            END DO
[4870]338            DO jk = 1, nlay_i-1
[2715]339               DO ji = kideb , kiut
[5682]340                  ztcond_i(ji,jk) = rcdic +                                                                       & 
341                     &                 0.09_wp * ( s_i_1d(ji,jk) + s_i_1d(ji,jk+1) )                              &
342                     &                 / MIN( -2._wp * epsi10, t_i_1d(ji,jk) + t_i_1d(ji,jk+1) - 2.0_wp * rt0 )   &
343                     &               - 0.0055_wp * ( t_i_1d(ji,jk) + t_i_1d(ji,jk+1) - 2.0 * rt0 ) 
[4870]344                  ztcond_i(ji,jk) = MAX( ztcond_i(ji,jk), zkimin )
[2715]345               END DO
346            END DO
347            DO ji = kideb , kiut
[5682]348               ztcond_i(ji,nlay_i) = rcdic + 0.090_wp * s_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji) - rt0 )   &
349                  &                        - 0.011_wp * ( t_bo_1d(ji) - rt0 ) 
[2715]350               ztcond_i(ji,nlay_i) = MAX( ztcond_i(ji,nlay_i), zkimin )
351            END DO
[921]352         ENDIF
[5682]353         
[921]354         !
355         !------------------------------------------------------------------------------|
[5682]356         !  5) G(he) - enhancement of thermal conductivity in mono-category case        |
[921]357         !------------------------------------------------------------------------------|
358         !
[5682]359         ! Computation of effective thermal conductivity G(h)
360         ! Used in mono-category case only to simulate an ITD implicitly
361         ! Fichefet and Morales Maqueda, JGR 1997
362
363         zghe(:) = 1._wp
364
365         SELECT CASE ( nn_monocat )
366
367         CASE (1,3) ! LIM3
368
369            zepsilon = 0.1_wp
370            zh_thres = EXP( 1._wp ) * zepsilon * 0.5_wp
371
372            DO ji = kideb, kiut
373   
374               ! Mean sea ice thermal conductivity
375               zkimean  = SUM( ztcond_i(ji,0:nlay_i) ) / REAL( nlay_i+1, wp )
376
377               ! Effective thickness he (zhe)
378               zfac     = 1._wp / ( rcdsn + zkimean )
379               zratio_s = rcdsn   * zfac
380               zratio_i = zkimean * zfac
381               zhe      = zratio_s * ht_i_1d(ji) + zratio_i * ht_s_1d(ji)
382
383               ! G(he)
384               rswitch  = MAX( 0._wp , SIGN( 1._wp , zhe - zh_thres ) )  ! =0 if zhe < zh_thres, if >
385               zghe(ji) = ( 1._wp - rswitch ) + rswitch * 0.5_wp * ( 1._wp + LOG( 2._wp * zhe / zepsilon ) )
386   
387               ! Impose G(he) < 2.
388               zghe(ji) = MIN( zghe(ji), 2._wp )
389
390            END DO
391
392         END SELECT
393
394         !
395         !------------------------------------------------------------------------------|
396         !  6) kappa factors                                                            |
397         !------------------------------------------------------------------------------|
398         !
399         !--- Snow
[921]400         DO ji = kideb, kiut
[5682]401            zfac                  =  1. / MAX( epsi10 , zh_s(ji) )
402            zkappa_s(ji,0)        = zghe(ji) * rcdsn * zfac
403            zkappa_s(ji,nlay_s)   = zghe(ji) * rcdsn * zfac
[921]404         END DO
[825]405
[4870]406         DO jk = 1, nlay_s-1
[921]407            DO ji = kideb , kiut
[5682]408               zkappa_s(ji,jk)    = zghe(ji) * 2.0 * rcdsn / MAX( epsi10, 2.0 * zh_s(ji) )
[921]409            END DO
410         END DO
[825]411
[5682]412         !--- Ice
[4870]413         DO jk = 1, nlay_i-1
[921]414            DO ji = kideb , kiut
[5682]415               zkappa_i(ji,jk)    = zghe(ji) * 2.0 * ztcond_i(ji,jk) / MAX( epsi10 , 2.0 * zh_i(ji) )
[921]416            END DO
417         END DO
[825]418
[5682]419         !--- Snow-ice interface
[921]420         DO ji = kideb , kiut
[5682]421            zfac                  = 1./ MAX( epsi10 , zh_i(ji) )
422            zkappa_i(ji,0)        = zghe(ji) * ztcond_i(ji,0) * zfac
423            zkappa_i(ji,nlay_i)   = zghe(ji) * ztcond_i(ji,nlay_i) * zfac
424            zkappa_s(ji,nlay_s)   = zghe(ji) * zghe(ji) * 2.0 * rcdsn * ztcond_i(ji,0) / & 
425           &                        MAX( epsi10, ( zghe(ji) * ztcond_i(ji,0) * zh_s(ji) + zghe(ji) * rcdsn * zh_i(ji) ) )
426            zkappa_i(ji,0)        = zkappa_s(ji,nlay_s) * isnow(ji) + zkappa_i(ji,0) * ( 1._wp - isnow(ji) )
[921]427         END DO
[5682]428
[921]429         !
430         !------------------------------------------------------------------------------|
[5682]431         ! 7) Sea ice specific heat, eta factors                                        |
[921]432         !------------------------------------------------------------------------------|
433         !
[4870]434         DO jk = 1, nlay_i
[921]435            DO ji = kideb , kiut
[4872]436               ztitemp(ji,jk)   = t_i_1d(ji,jk)
[5682]437               zspeche_i(ji,jk) = cpic + zgamma * s_i_1d(ji,jk) / MAX( ( t_i_1d(ji,jk) - rt0 ) * ( ztib(ji,jk) - rt0 ), epsi10 )
438               zeta_i(ji,jk)    = rdt_ice / MAX( rhoic * zspeche_i(ji,jk) * zh_i(ji), epsi10 )
[921]439            END DO
440         END DO
[825]441
[4870]442         DO jk = 1, nlay_s
[921]443            DO ji = kideb , kiut
[4872]444               ztstemp(ji,jk) = t_s_1d(ji,jk)
[5682]445               zeta_s(ji,jk)  = rdt_ice / MAX( rhosn * cpic * zh_s(ji), epsi10 )
[921]446            END DO
447         END DO
[5682]448
[921]449         !
450         !------------------------------------------------------------------------------|
[5682]451         ! 8) surface flux computation                                                  |
[921]452         !------------------------------------------------------------------------------|
453         !
[5682]454         IF ( ln_it_qnsice ) THEN
[4990]455            DO ji = kideb , kiut
456               ! update of the non solar flux according to the update in T_su
457               qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_1d(ji) - ztsubit(ji) )
458            END DO
459         ENDIF
460
461         ! Update incoming flux
[921]462         DO ji = kideb , kiut
463            ! update incoming flux
[5682]464            zf(ji)    =          zfsw(ji)              & ! net absorbed solar radiation
465               &         + qns_ice_1d(ji)                ! non solar total flux (LWup, LWdw, SH, LH)
[921]466         END DO
[825]467
[921]468         !
469         !------------------------------------------------------------------------------|
[5682]470         ! 9) tridiagonal system terms                                                  |
[921]471         !------------------------------------------------------------------------------|
472         !
473         !!layer denotes the number of the layer in the snow or in the ice
474         !!numeq denotes the reference number of the equation in the tridiagonal
475         !!system, terms of tridiagonal system are indexed as following :
476         !!1 is subdiagonal term, 2 is diagonal and 3 is superdiagonal one
[825]477
[921]478         !!ice interior terms (top equation has the same form as the others)
479
[4873]480         DO numeq=1,nlay_i+3
[921]481            DO ji = kideb , kiut
482               ztrid(ji,numeq,1) = 0.
483               ztrid(ji,numeq,2) = 0.
484               ztrid(ji,numeq,3) = 0.
[5682]485               zindterm(ji,numeq)= 0.
486               zindtbis(ji,numeq)= 0.
[921]487               zdiagbis(ji,numeq)= 0.
488            ENDDO
489         ENDDO
490
491         DO numeq = nlay_s + 2, nlay_s + nlay_i 
492            DO ji = kideb , kiut
[5682]493               jk                 = numeq - nlay_s - 1
494               ztrid(ji,numeq,1)  =  - zeta_i(ji,jk) * zkappa_i(ji,jk-1)
495               ztrid(ji,numeq,2)  =  1.0 + zeta_i(ji,jk) * ( zkappa_i(ji,jk-1) + zkappa_i(ji,jk) )
496               ztrid(ji,numeq,3)  =  - zeta_i(ji,jk) * zkappa_i(ji,jk)
497               zindterm(ji,numeq) =  ztib(ji,jk) + zeta_i(ji,jk) * zradab_i(ji,jk)
[921]498            END DO
499         ENDDO
500
501         numeq =  nlay_s + nlay_i + 1
502         DO ji = kideb , kiut
[825]503            !!ice bottom term
504            ztrid(ji,numeq,1)  =  - zeta_i(ji,nlay_i)*zkappa_i(ji,nlay_i-1)   
[5682]505            ztrid(ji,numeq,2)  =  1.0 + zeta_i(ji,nlay_i) * ( zkappa_i(ji,nlay_i) * zg1 + zkappa_i(ji,nlay_i-1) )
[825]506            ztrid(ji,numeq,3)  =  0.0
[5682]507            zindterm(ji,numeq) =  ztib(ji,nlay_i) + zeta_i(ji,nlay_i) *  &
508               &                  ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i) * zg1 *  t_bo_1d(ji) ) 
[921]509         ENDDO
[825]510
511
[921]512         DO ji = kideb , kiut
[5682]513            IF ( ht_s_1d(ji) > 0.0 ) THEN
[921]514               !
515               !------------------------------------------------------------------------------|
516               !  snow-covered cells                                                          |
517               !------------------------------------------------------------------------------|
518               !
519               !!snow interior terms (bottom equation has the same form as the others)
520               DO numeq = 3, nlay_s + 1
[5682]521                  jk                  =  numeq - 1
522                  ztrid(ji,numeq,1)   =  - zeta_s(ji,jk) * zkappa_s(ji,jk-1)
523                  ztrid(ji,numeq,2)   =  1.0 + zeta_s(ji,jk) * ( zkappa_s(ji,jk-1) + zkappa_s(ji,jk) )
[4870]524                  ztrid(ji,numeq,3)   =  - zeta_s(ji,jk)*zkappa_s(ji,jk)
[5682]525                  zindterm(ji,numeq)  =  ztsb(ji,jk) + zeta_s(ji,jk) * zradab_s(ji,jk)
[921]526               END DO
[825]527
[921]528               !!case of only one layer in the ice (ice equation is altered)
529               IF ( nlay_i.eq.1 ) THEN
530                  ztrid(ji,nlay_s+2,3)    =  0.0
[5682]531                  zindterm(ji,nlay_s+2)   =  zindterm(ji,nlay_s+2) + zkappa_i(ji,1) * t_bo_1d(ji) 
[921]532               ENDIF
[834]533
[5682]534               IF ( t_su_1d(ji) < rt0 ) THEN
[825]535
[921]536                  !------------------------------------------------------------------------------|
537                  !  case 1 : no surface melting - snow present                                  |
538                  !------------------------------------------------------------------------------|
539                  zdifcase(ji)    =  1.0
540                  numeqmin(ji)    =  1
541                  numeqmax(ji)    =  nlay_i + nlay_s + 1
[825]542
[921]543                  !!surface equation
[5682]544                  ztrid(ji,1,1)  = 0.0
545                  ztrid(ji,1,2)  = dzf(ji) - zg1s * zkappa_s(ji,0)
546                  ztrid(ji,1,3)  = zg1s * zkappa_s(ji,0)
547                  zindterm(ji,1) = dzf(ji) * t_su_1d(ji) - zf(ji)
[825]548
[921]549                  !!first layer of snow equation
[5682]550                  ztrid(ji,2,1)  =  - zkappa_s(ji,0) * zg1s * zeta_s(ji,1)
551                  ztrid(ji,2,2)  =  1.0 + zeta_s(ji,1) * ( zkappa_s(ji,1) + zkappa_s(ji,0) * zg1s )
[921]552                  ztrid(ji,2,3)  =  - zeta_s(ji,1)* zkappa_s(ji,1)
[5682]553                  zindterm(ji,2) =  ztsb(ji,1) + zeta_s(ji,1) * zradab_s(ji,1)
[825]554
[921]555               ELSE 
556                  !
557                  !------------------------------------------------------------------------------|
558                  !  case 2 : surface is melting - snow present                                  |
559                  !------------------------------------------------------------------------------|
560                  !
561                  zdifcase(ji)    =  2.0
562                  numeqmin(ji)    =  2
563                  numeqmax(ji)    =  nlay_i + nlay_s + 1
[825]564
[921]565                  !!first layer of snow equation
566                  ztrid(ji,2,1)  =  0.0
[5682]567                  ztrid(ji,2,2)  =  1.0 + zeta_s(ji,1) * ( zkappa_s(ji,1) + zkappa_s(ji,0) * zg1s )
[921]568                  ztrid(ji,2,3)  =  - zeta_s(ji,1)*zkappa_s(ji,1) 
[5682]569                  zindterm(ji,2) = ztsb(ji,1) + zeta_s(ji,1) *  &
570                     &             ( zradab_s(ji,1) + zkappa_s(ji,0) * zg1s * t_su_1d(ji) ) 
[921]571               ENDIF
572            ELSE
573               !
574               !------------------------------------------------------------------------------|
575               !  cells without snow                                                          |
576               !------------------------------------------------------------------------------|
577               !
[5682]578               IF ( t_su_1d(ji) < rt0 ) THEN
[921]579                  !
580                  !------------------------------------------------------------------------------|
581                  !  case 3 : no surface melting - no snow                                       |
582                  !------------------------------------------------------------------------------|
583                  !
584                  zdifcase(ji)      =  3.0
585                  numeqmin(ji)      =  nlay_s + 1
586                  numeqmax(ji)      =  nlay_i + nlay_s + 1
[825]587
[921]588                  !!surface equation   
589                  ztrid(ji,numeqmin(ji),1)   =  0.0
590                  ztrid(ji,numeqmin(ji),2)   =  dzf(ji) - zkappa_i(ji,0)*zg1   
591                  ztrid(ji,numeqmin(ji),3)   =  zkappa_i(ji,0)*zg1
[5682]592                  zindterm(ji,numeqmin(ji))  =  dzf(ji)*t_su_1d(ji) - zf(ji)
[825]593
[921]594                  !!first layer of ice equation
595                  ztrid(ji,numeqmin(ji)+1,1) =  - zkappa_i(ji,0) * zg1 * zeta_i(ji,1)
[5682]596                  ztrid(ji,numeqmin(ji)+1,2) =  1.0 + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 )
597                  ztrid(ji,numeqmin(ji)+1,3) =  - zeta_i(ji,1) * zkappa_i(ji,1) 
598                  zindterm(ji,numeqmin(ji)+1)=  ztib(ji,1) + zeta_i(ji,1) * zradab_i(ji,1) 
[825]599
[921]600                  !!case of only one layer in the ice (surface & ice equations are altered)
[825]601
[5682]602                  IF ( nlay_i == 1 ) THEN
[921]603                     ztrid(ji,numeqmin(ji),1)    =  0.0
[5682]604                     ztrid(ji,numeqmin(ji),2)    =  dzf(ji) - zkappa_i(ji,0) * 2.0
605                     ztrid(ji,numeqmin(ji),3)    =  zkappa_i(ji,0) * 2.0
606                     ztrid(ji,numeqmin(ji)+1,1)  =  -zkappa_i(ji,0) * 2.0 * zeta_i(ji,1)
607                     ztrid(ji,numeqmin(ji)+1,2)  =  1.0 + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2.0 + zkappa_i(ji,1) )
[921]608                     ztrid(ji,numeqmin(ji)+1,3)  =  0.0
[825]609
[5682]610                     zindterm(ji,numeqmin(ji)+1) =  ztib(ji,1) + zeta_i(ji,1) *  &
611                        &                          ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) )
[921]612                  ENDIF
[825]613
[921]614               ELSE
[825]615
[921]616                  !
617                  !------------------------------------------------------------------------------|
618                  ! case 4 : surface is melting - no snow                                        |
619                  !------------------------------------------------------------------------------|
620                  !
621                  zdifcase(ji)    =  4.0
622                  numeqmin(ji)    =  nlay_s + 2
623                  numeqmax(ji)    =  nlay_i + nlay_s + 1
[825]624
[921]625                  !!first layer of ice equation
626                  ztrid(ji,numeqmin(ji),1)      =  0.0
[5682]627                  ztrid(ji,numeqmin(ji),2)      =  1.0 + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 ) 
[921]628                  ztrid(ji,numeqmin(ji),3)      =  - zeta_i(ji,1) * zkappa_i(ji,1)
[5682]629                  zindterm(ji,numeqmin(ji))     =  ztib(ji,1) + zeta_i(ji,1) *  &
630                     &                             ( zradab_i(ji,1) + zkappa_i(ji,0) * zg1 * t_su_1d(ji) ) 
[825]631
[921]632                  !!case of only one layer in the ice (surface & ice equations are altered)
[5682]633                  IF ( nlay_i == 1 ) THEN
[921]634                     ztrid(ji,numeqmin(ji),1)  =  0.0
[5682]635                     ztrid(ji,numeqmin(ji),2)  =  1.0 + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2.0 + zkappa_i(ji,1) )
[921]636                     ztrid(ji,numeqmin(ji),3)  =  0.0
[5682]637                     zindterm(ji,numeqmin(ji)) =  ztib(ji,1) + zeta_i(ji,1) * ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) )  &
638                        &                       + t_su_1d(ji) * zeta_i(ji,1) * zkappa_i(ji,0) * 2.0
[921]639                  ENDIF
[825]640
[921]641               ENDIF
642            ENDIF
[825]643
[921]644         END DO
[825]645
[921]646         !
647         !------------------------------------------------------------------------------|
[5682]648         ! 10) tridiagonal system solving                                               |
[921]649         !------------------------------------------------------------------------------|
650         !
[825]651
[921]652         ! Solve the tridiagonal system with Gauss elimination method.
653         ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON,
654         ! McGraw-Hill 1984. 
[825]655
[921]656         maxnumeqmax = 0
[4873]657         minnumeqmin = nlay_i+5
[825]658
[921]659         DO ji = kideb , kiut
[5682]660            zindtbis(ji,numeqmin(ji)) =  zindterm(ji,numeqmin(ji))
[921]661            zdiagbis(ji,numeqmin(ji)) =  ztrid(ji,numeqmin(ji),2)
662            minnumeqmin               =  MIN(numeqmin(ji),minnumeqmin)
663            maxnumeqmax               =  MAX(numeqmax(ji),maxnumeqmax)
664         END DO
665
[4870]666         DO jk = minnumeqmin+1, maxnumeqmax
[921]667            DO ji = kideb , kiut
[4870]668               numeq               =  min(max(numeqmin(ji)+1,jk),numeqmax(ji))
[5682]669               zdiagbis(ji,numeq)  =  ztrid(ji,numeq,2)  - ztrid(ji,numeq,1) * ztrid(ji,numeq-1,3)  / zdiagbis(ji,numeq-1)
670               zindtbis(ji,numeq)  =  zindterm(ji,numeq) - ztrid(ji,numeq,1) * zindtbis(ji,numeq-1) / zdiagbis(ji,numeq-1)
[921]671            END DO
672         END DO
673
674         DO ji = kideb , kiut
675            ! ice temperatures
[5682]676            t_i_1d(ji,nlay_i)    =  zindtbis(ji,numeqmax(ji)) / zdiagbis(ji,numeqmax(ji))
[921]677         END DO
678
[5682]679         DO numeq = nlay_i + nlay_s, nlay_s + 2, -1
[921]680            DO ji = kideb , kiut
[4870]681               jk    =  numeq - nlay_s - 1
[5682]682               t_i_1d(ji,jk)  =  ( zindtbis(ji,numeq) - ztrid(ji,numeq,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,numeq)
[921]683            END DO
684         END DO
685
686         DO ji = kideb , kiut
[825]687            ! snow temperatures     
[5682]688            IF (ht_s_1d(ji) > 0._wp) &
689               t_s_1d(ji,nlay_s)     =  ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) )  &
690               &                        / zdiagbis(ji,nlay_s+1) * MAX( 0.0, SIGN( 1.0, ht_s_1d(ji) ) ) 
[825]691
692            ! surface temperature
[5682]693            isnow(ji)   = 1._wp - MAX( 0._wp , SIGN( 1._wp , -ht_s_1d(ji) ) )
[4872]694            ztsubit(ji) = t_su_1d(ji)
[5682]695            IF( t_su_1d(ji) < rt0 ) &
696               t_su_1d(ji) = ( zindtbis(ji,numeqmin(ji)) - ztrid(ji,numeqmin(ji),3) *  &
697               &             ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) ) ) / zdiagbis(ji,numeqmin(ji)) 
[921]698         END DO
699         !
700         !--------------------------------------------------------------------------
[5682]701         !  11) Has the scheme converged ?, end of the iterative procedure         |
[921]702         !--------------------------------------------------------------------------
703         !
704         ! check that nowhere it has started to melt
[5682]705         ! zerrit(ji) is a measure of error, it has to be under terr_dif
[921]706         DO ji = kideb , kiut
[5682]707            t_su_1d(ji) =  MAX(  MIN( t_su_1d(ji) , rt0 ) , 190._wp  )
708            zerrit(ji)  =  ABS( t_su_1d(ji) - ztsubit(ji) )     
[921]709         END DO
[825]710
[4870]711         DO jk  =  1, nlay_s
[921]712            DO ji = kideb , kiut
[5682]713               t_s_1d(ji,jk) = MAX(  MIN( t_s_1d(ji,jk), rt0 ), 190._wp  )
714               zerrit(ji)    = MAX( zerrit(ji), ABS( t_s_1d(ji,jk) - ztstemp(ji,jk) ) )
[921]715            END DO
716         END DO
[825]717
[4870]718         DO jk  =  1, nlay_i
[921]719            DO ji = kideb , kiut
[5682]720               ztmelt_i      = -tmut * s_i_1d(ji,jk) + rt0 
721               t_i_1d(ji,jk) =  MAX( MIN( t_i_1d(ji,jk), ztmelt_i ), 190._wp )
722               zerrit(ji)    =  MAX( zerrit(ji), ABS( t_i_1d(ji,jk) - ztitemp(ji,jk) ) )
[921]723            END DO
724         END DO
[825]725
[921]726         ! Compute spatial maximum over all errors
[2715]727         ! note that this could be optimized substantially by iterating only the non-converging points
728         zerritmax = 0._wp
729         DO ji = kideb, kiut
730            zerritmax = MAX( zerritmax, zerrit(ji) )   
[921]731         END DO
[2715]732         IF( lk_mpp ) CALL mpp_max( zerritmax, kcom=ncomm_ice )
[825]733
734      END DO  ! End of the do while iterative procedure
735
[5682]736      IF( ln_icectl .AND. lwp ) THEN
[1055]737         WRITE(numout,*) ' zerritmax : ', zerritmax
738         WRITE(numout,*) ' nconv     : ', nconv
739      ENDIF
[825]740
[921]741      !
[2715]742      !-------------------------------------------------------------------------!
[5682]743      !   12) Fluxes at the interfaces                                          !
[2715]744      !-------------------------------------------------------------------------!
[921]745      DO ji = kideb, kiut
[2715]746         !                                ! surface ice conduction flux
[5682]747         isnow(ji)       = 1._wp - MAX( 0._wp, SIGN( 1._wp, -ht_s_1d(ji) ) )
748         fc_su(ji)       =  -           isnow(ji)   * zkappa_s(ji,0) * zg1s * (t_s_1d(ji,1) - t_su_1d(ji))   &
749            &               - ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1  * (t_i_1d(ji,1) - t_su_1d(ji))
[2715]750         !                                ! bottom ice conduction flux
[4872]751         fc_bo_i(ji)     =  - zkappa_i(ji,nlay_i) * ( zg1*(t_bo_1d(ji) - t_i_1d(ji,nlay_i)) )
[921]752      END DO
[825]753
[5682]754      ! --- computes sea ice energy of melting compulsory for limthd_dh --- !
755      CALL lim_thd_enmelt( kideb, kiut )
756
757      ! --- diagnose the change in non-solar flux due to surface temperature change --- !
758      IF ( ln_it_qnsice ) THEN
759         DO ji = kideb, kiut
760            hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) - ( qns_ice_1d(ji)  - zqns_ice_b(ji) ) * a_i_1d(ji) 
761         END DO
762      END IF
763
764      ! --- diag conservation imbalance on heat diffusion - PART 2 --- !
765      DO ji = kideb, kiut
766         zdq(ji)        = - zq_ini(ji) + ( SUM( q_i_1d(ji,1:nlay_i) ) * ht_i_1d(ji) * r1_nlay_i +  &
767            &                              SUM( q_s_1d(ji,1:nlay_s) ) * ht_s_1d(ji) * r1_nlay_s )
768         IF( t_su_1d(ji) < rt0 ) THEN  ! case T_su < 0degC
769            zhfx_err(ji) = qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq(ji) * r1_rdtice 
770         ELSE                          ! case T_su = 0degC
771            zhfx_err(ji) = fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq(ji) * r1_rdtice 
772         ENDIF
773         hfx_err_1d(ji) = hfx_err_1d(ji) + zhfx_err(ji) * a_i_1d(ji)
774
775         ! total heat that is sent to the ocean (i.e. not used in the heat diffusion equation)
776         hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) + zhfx_err(ji) * a_i_1d(ji)
777      END DO 
778
[4688]779      !-----------------------------------------
780      ! Heat flux used to warm/cool ice in W.m-2
781      !-----------------------------------------
782      DO ji = kideb, kiut
[5682]783         IF( t_su_1d(ji) < rt0 ) THEN  ! case T_su < 0degC
[4765]784            hfx_dif_1d(ji) = hfx_dif_1d(ji)  +   &
[4872]785               &            ( qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_1d(ji)
[5682]786         ELSE                          ! case T_su = 0degC
[4765]787            hfx_dif_1d(ji) = hfx_dif_1d(ji) +    &
[4872]788               &             ( fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_1d(ji)
[4688]789         ENDIF
[5682]790         ! correction on the diagnosed heat flux due to non-convergence of the algorithm used to solve heat equation
791         hfx_dif_1d(ji) = hfx_dif_1d(ji) - zhfx_err(ji) * a_i_1d(ji)
792      END DO   
[4688]793      !
[5682]794      CALL wrk_dealloc( jpij, numeqmin, numeqmax )
795      CALL wrk_dealloc( jpij, isnow, ztsub, ztsubit, zh_i, zh_s, zfsw )
796      CALL wrk_dealloc( jpij, zf, dzf, zqns_ice_b, zerrit, zdifcase, zftrice, zihic, zghe )
797      CALL wrk_dealloc( jpij,nlay_i+1, ztcond_i, zradtr_i, zradab_i, zkappa_i, ztib, zeta_i, ztitemp, z_i, zspeche_i, kjstart = 0 )
798      CALL wrk_dealloc( jpij,nlay_s+1,           zradtr_s, zradab_s, zkappa_s, ztsb, zeta_s, ztstemp, z_s, kjstart = 0 )
799      CALL wrk_dealloc( jpij,nlay_i+3, zindterm, zindtbis, zdiagbis )
800      CALL wrk_dealloc( jpij,nlay_i+3,3, ztrid )
[4990]801      CALL wrk_dealloc( jpij, zdq, zq_ini, zhfx_err )
[4688]802
803   END SUBROUTINE lim_thd_dif
804
805   SUBROUTINE lim_thd_enmelt( kideb, kiut )
806      !!-----------------------------------------------------------------------
807      !!                   ***  ROUTINE lim_thd_enmelt ***
808      !!                 
809      !! ** Purpose :   Computes sea ice energy of melting q_i (J.m-3) from temperature
810      !!
811      !! ** Method  :   Formula (Bitz and Lipscomb, 1999)
812      !!-------------------------------------------------------------------
813      INTEGER, INTENT(in) ::   kideb, kiut   ! bounds for the spatial loop
814      !
815      INTEGER  ::   ji, jk   ! dummy loop indices
[4990]816      REAL(wp) ::   ztmelts  ! local scalar
[4688]817      !!-------------------------------------------------------------------
818      !
819      DO jk = 1, nlay_i             ! Sea ice energy of melting
[921]820         DO ji = kideb, kiut
[5682]821            ztmelts      = - tmut  * s_i_1d(ji,jk) + rt0
822            t_i_1d(ji,jk) = MIN( t_i_1d(ji,jk), ztmelts ) ! Force t_i_1d to be lower than melting point
823                                                          !   (sometimes dif scheme produces abnormally high temperatures)   
824            q_i_1d(ji,jk) = rhoic * ( cpic * ( ztmelts - t_i_1d(ji,jk) )                           &
825               &                    + lfus * ( 1.0 - ( ztmelts-rt0 ) / ( t_i_1d(ji,jk) - rt0 ) )   &
826               &                    - rcp  *         ( ztmelts-rt0 )  ) 
[921]827         END DO
[4688]828      END DO
829      DO jk = 1, nlay_s             ! Snow energy of melting
830         DO ji = kideb, kiut
[5682]831            q_s_1d(ji,jk) = rhosn * ( cpic * ( rt0 - t_s_1d(ji,jk) ) + lfus )
[921]832         END DO
[4688]833      END DO
[2715]834      !
[4688]835   END SUBROUTINE lim_thd_enmelt
[825]836
837#else
[2715]838   !!----------------------------------------------------------------------
839   !!                   Dummy Module                 No LIM-3 sea-ice model
840   !!----------------------------------------------------------------------
[825]841CONTAINS
842   SUBROUTINE lim_thd_dif          ! Empty routine
843   END SUBROUTINE lim_thd_dif
844#endif
[2528]845   !!======================================================================
[921]846END MODULE limthd_dif
Note: See TracBrowser for help on using the repository browser.