source: branches/2017/dev_r7881_no_wrk_alloc/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90 @ 7910

Last change on this file since 7910 was 7910, checked in by timgraham, 3 years ago

All wrk_alloc removed

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