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/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90 @ 8420

Last change on this file since 8420 was 8420, checked in by clem, 7 years ago

changing names

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