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_ent.F90 in branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_ent.F90 @ 4332

Last change on this file since 4332 was 4332, checked in by clem, 10 years ago

update LIM3 to fix remaining bugs. Now working in global and regional config.
  • Property svn:keywords set to Id
File size: 27.5 KB
Line 
1MODULE limthd_ent
2   !!======================================================================
3   !!                       ***  MODULE limthd_ent   ***
4   !!                  Redistribution of Enthalpy in the ice
5   !!                        on the new vertical grid
6   !!                       after vertical growth/decay
7   !!======================================================================
8   !! History :  LIM  ! 2003-05 (M. Vancoppenolle) Original code in 1D
9   !!                 ! 2005-07 (M. Vancoppenolle) 3D version
10   !!                 ! 2006-11 (X. Fettweis) Vectorized
11   !!            3.0  ! 2008-03 (M. Vancoppenolle) Energy conservation and clean code
12   !!            4.0  ! 2011-02 (G. Madec) dynamical allocation
13   !!----------------------------------------------------------------------
14#if defined key_lim3
15   !!----------------------------------------------------------------------
16   !!   'key_lim3'                                      LIM3 sea-ice model
17   !!----------------------------------------------------------------------
18   !!   lim_thd_ent   : ice redistribution of enthalpy
19   !!----------------------------------------------------------------------
20   USE par_oce        ! ocean parameters
21   USE dom_oce        ! domain variables
22   USE domain         !
23   USE phycst         ! physical constants
24   USE ice            ! LIM variables
25   USE par_ice        ! LIM parameters
26   USE thd_ice        ! LIM thermodynamics
27   USE limvar         ! LIM variables
28   USE in_out_manager ! I/O manager
29   USE lib_mpp        ! MPP library
30   USE wrk_nemo       ! work arrays
31   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
32
33   IMPLICIT NONE
34   PRIVATE
35
36   PUBLIC   lim_thd_ent         ! called by lim_thd
37
38   REAL(wp) ::   epsi20 = 1.e-20_wp   ! constant values
39   REAL(wp) ::   epsi10 = 1.e-10_wp   !
40   REAL(wp) ::   zzero  = 0._wp      !
41   REAL(wp) ::   zone   = 1._wp      !
42
43   !!----------------------------------------------------------------------
44   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
45   !! $Id$
46   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
47   !!----------------------------------------------------------------------
48CONTAINS
49 
50   SUBROUTINE lim_thd_ent( kideb, kiut, jl )
51      !!-------------------------------------------------------------------
52      !!               ***   ROUTINE lim_thd_ent  ***
53      !!
54      !! ** Purpose :
55      !!           This routine computes new vertical grids
56      !!           in the ice and in the snow, and consistently redistributes
57      !!           temperatures in the snow / ice.
58      !!           Redistribution is made so as to ensure to energy conservation
59      !!
60      !!
61      !! ** Method  : linear conservative remapping
62      !!           
63      !! ** Steps : 1) Grid
64      !!            2) Switches
65      !!            3) Snow redistribution
66      !!            4) Ice enthalpy redistribution
67      !!            5) Ice salinity, recover temperature
68      !!
69      !! References : Bitz & Lipscomb, JGR 99; Vancoppenolle et al., GRL, 2005
70      !!-------------------------------------------------------------------
71      INTEGER , INTENT(in) ::   kideb, kiut   ! Start/End point on which the  the computation is applied
72      INTEGER , INTENT(in) ::   jl            ! Thickness cateogry number
73
74      INTEGER ::   ji,jk   !  dummy loop indices
75      INTEGER ::   ii, ij       ,   &  !  dummy indices
76         ntop0          ,   &  !  old layer top index
77         nbot1          ,   &  !  new layer bottom index
78         ntop1          ,   &  !  new layer top index
79         limsum         ,   &  !  temporary loop index
80         nlayi0,nlays0  ,   &  !  old number of layers
81         maxnbot0       ,   &  !  old layer bottom index
82         layer0, layer1        !  old/new layer indexes
83
84
85      REAL(wp) :: &
86         ztmelts        ,   &  ! ice melting point
87         zqsnic         ,   &  ! enthalpy of snow ice layer
88         zhsnow         ,   &  ! temporary snow thickness variable
89         zswitch        ,   &  ! dummy switch argument
90         zfac1          ,   &  ! dummy factor
91         zfac2          ,   &  ! dummy factor
92         ztform         ,   &  !: bottom formation temperature
93         zaaa           ,   &  !: dummy factor
94         zbbb           ,   &  !: dummy factor
95         zccc           ,   &  !: dummy factor
96         zdiscrim              !: dummy factor
97
98      INTEGER, POINTER, DIMENSION(:) ::   snswi     !  snow switch
99      INTEGER, POINTER, DIMENSION(:) ::   nbot0     !  old layer bottom index
100      INTEGER, POINTER, DIMENSION(:) ::   icsuind   !  ice surface index
101      INTEGER, POINTER, DIMENSION(:) ::   icsuswi   !  ice surface switch
102      INTEGER, POINTER, DIMENSION(:) ::   icboind   !  ice bottom index
103      INTEGER, POINTER, DIMENSION(:) ::   icboswi   !  ice bottom switch
104      INTEGER, POINTER, DIMENSION(:) ::   snicind   !  snow ice index
105      INTEGER, POINTER, DIMENSION(:) ::   snicswi   !  snow ice switch
106      INTEGER, POINTER, DIMENSION(:) ::   snind     !  snow index
107      !
108      REAL(wp), POINTER, DIMENSION(:) ::   zh_i   ! thickness of an ice layer
109      REAL(wp), POINTER, DIMENSION(:) ::   zh_s          ! thickness of a snow layer
110      REAL(wp), POINTER, DIMENSION(:) ::   zqsnow        ! enthalpy of the snow put in snow ice   
111      REAL(wp), POINTER, DIMENSION(:) ::   zdeltah       ! temporary variable
112      REAL(wp), POINTER, DIMENSION(:) ::   zqti_in, zqts_in
113      REAL(wp), POINTER, DIMENSION(:) ::   zqti_fin, zqts_fin
114
115      REAL(wp), POINTER, DIMENSION(:,:) ::   zm0       !  old layer-system vertical cotes
116      REAL(wp), POINTER, DIMENSION(:,:) ::   qm0       !  old layer-system heat content
117      REAL(wp), POINTER, DIMENSION(:,:) ::   z_s       !  new snow system vertical cotes
118      REAL(wp), POINTER, DIMENSION(:,:) ::   z_i       !  new ice system vertical cotes
119      REAL(wp), POINTER, DIMENSION(:,:) ::   zthick0   !  old ice thickness
120      REAL(wp), POINTER, DIMENSION(:,:) ::   zhl0      ! old and new layer thicknesses
121      REAL(wp), POINTER, DIMENSION(:,:) ::   zrl01
122
123      REAL(wp) ::   zinda 
124      !!-------------------------------------------------------------------
125
126      CALL wrk_alloc( jpij, snswi, nbot0, icsuind, icsuswi, icboind, icboswi, snicind, snicswi, snind )   ! integer
127      CALL wrk_alloc( jpij, zh_i, zh_s, zqsnow, zdeltah, zqti_in, zqts_in, zqti_fin, zqts_fin )           ! real
128      CALL wrk_alloc( jpij,jkmax+4, zm0, qm0, z_s, z_i, zthick0, zhl0, kjstart = 0 )
129      CALL wrk_alloc( jkmax+4,jkmax+4, zrl01, kistart = 0, kjstart = 0 )
130
131      zthick0(:,:) = 0._wp
132      zm0    (:,:) = 0._wp
133      qm0    (:,:) = 0._wp
134      zrl01  (:,:) = 0._wp
135      zhl0   (:,:) = 0._wp
136      z_i    (:,:) = 0._wp
137      z_s    (:,:) = 0._wp
138
139      !
140      !------------------------------------------------------------------------------|
141      !  1) Grid                                                                     |
142      !------------------------------------------------------------------------------|
143      nlays0 = nlay_s
144      nlayi0 = nlay_i
145
146      DO ji = kideb, kiut
147         zh_i(ji) = old_ht_i_b(ji) / REAL( nlay_i ) 
148         zh_s(ji) = old_ht_s_b(ji) / REAL( nlay_s )
149      END DO
150
151      !
152      !------------------------------------------------------------------------------|
153      !  2) Switches                                                                 |
154      !------------------------------------------------------------------------------|
155      ! 2.1 snind(ji), snswi(ji)
156      ! snow surface behaviour : computation of snind(ji)-snswi(ji)
157      ! snind(ji) : index which equals
158      !   0 if snow is accumulating
159      !   1 if 1st layer is melting
160      !   2 if 2nd layer is melting ...
161      DO ji = kideb, kiut
162         snind  (ji) = 0
163         zdeltah(ji) = 0._wp
164      ENDDO !ji
165
166      DO jk = 1, nlays0
167         DO ji = kideb, kiut
168            snind(ji)  = jk        *      NINT(MAX(0.0,SIGN(1.0,-dh_s_tot(ji)-zdeltah(ji)))) &
169               + snind(ji) * (1 - NINT(MAX(0.0,SIGN(1.0,-dh_s_tot(ji)-zdeltah(ji)))))
170            zdeltah(ji)= zdeltah(ji) + zh_s(ji)
171         END DO ! ji
172      END DO ! jk
173
174      ! snswi(ji) : switch which value equals 1 if snow melts
175      !              0 if not
176      DO ji = kideb, kiut
177         snswi(ji)     = MAX(0,NINT(-dh_s_tot(ji)/MAX(epsi20,ABS(dh_s_tot(ji)))))
178      END DO ! ji
179
180      ! 2.2 icsuind(ji), icsuswi(ji)
181      ! ice surface behaviour : computation of icsuind(ji)-icsuswi(ji)
182      ! icsuind(ji) : index which equals
183      !     0 if nothing happens at the surface
184      !     1 if first layer is melting
185      !     2 if 2nd layer is reached by melt ...
186      DO ji = kideb, kiut
187         icsuind(ji) = 0
188         zdeltah(ji) = 0._wp
189      END DO !ji
190      DO jk = 1, nlayi0
191         DO ji = kideb, kiut
192            icsuind(ji) = jk          *      NINT(MAX(0.0,SIGN(1.0,-dh_i_surf(ji)-zdeltah(ji)))) &
193               + icsuind(ji) * (1 - NINT(MAX(0.0,SIGN(1.0,-dh_i_surf(ji)-zdeltah(ji)))))
194            zdeltah(ji) = zdeltah(ji) + zh_i(ji)
195         END DO ! ji
196      ENDDO !jk
197
198      ! icsuswi(ji) : switch which equals
199      !     1 if ice melts at the surface
200      !     0 if not
201      DO ji = kideb, kiut
202         icsuswi(ji)  = MAX(0,NINT(-dh_i_surf(ji)/MAX(epsi20 , ABS(dh_i_surf(ji)) ) ) )
203      ENDDO
204
205      ! 2.3 icboind(ji), icboswi(ji)
206      ! ice bottom behaviour : computation of icboind(ji)-icboswi(ji)
207      ! icboind(ji) : index which equals
208      !     0 if accretion is on the way
209      !     1 if last layer has started to melt
210      !     2 if penultiem layer is melting ... and so on
211      !            N+1 if all layers melt and that snow transforms into ice
212      DO ji = kideb, kiut 
213         icboind(ji) = 0
214         zdeltah(ji) = 0._wp
215      END DO
216      DO jk = nlayi0, 1, -1
217         DO ji = kideb, kiut
218            icboind(ji) = (nlayi0+1-jk) *      NINT(MAX(0.0,SIGN(1.0,-dh_i_bott(ji)-zdeltah(ji)))) &
219               &          + icboind(ji) * (1 - NINT(MAX(0.0,SIGN(1.0,-dh_i_bott(ji)-zdeltah(ji))))) 
220            zdeltah(ji) = zdeltah(ji) + zh_i(ji)
221         END DO
222      END DO
223
224      DO ji = kideb, kiut
225         ! case of total ablation with remaining snow
226         IF ( ( ht_i_b(ji) .GT. epsi20 ) .AND. &
227            ( ht_i_b(ji) - dh_snowice(ji) .LT. epsi20 ) ) icboind(ji) = nlay_i + 1
228      END DO
229
230      ! icboswi(ji) : switch which equals
231      !     1 if ice accretion is on the way
232      !     0 if ablation is on the way
233      DO ji = kideb, kiut 
234         icboswi(ji) = MAX(0,NINT(dh_i_bott(ji) / MAX(epsi20,ABS(dh_i_bott(ji)))))
235      END DO
236
237      ! 2.4 snicind(ji), snicswi(ji)
238      ! snow ice formation : calcul de snicind(ji)-snicswi(ji)
239      ! snicind(ji) : index which equals
240      !     0 if no snow-ice forms
241      !     1 if last layer of snow has started to melt
242      !     2 if penultiem layer ...
243      DO ji = kideb, kiut
244         snicind(ji) = 0
245         zdeltah(ji) = 0._wp
246      END DO
247      DO jk = nlays0, 1, -1
248         DO ji = kideb, kiut
249            snicind(ji) = (nlays0+1-jk) &
250               *      NINT(MAX(0.0,SIGN(1.0,dh_snowice(ji)-zdeltah(ji)))) + snicind(ji)   &
251               * (1 - NINT(MAX(0.0,SIGN(1.0,dh_snowice(ji)-zdeltah(ji)))))
252            zdeltah(ji) = zdeltah(ji) + zh_s(ji)
253         END DO
254      END DO
255
256      ! snicswi(ji) : switch which equals
257      !     1 if snow-ice forms
258      !     0 if not
259      DO ji = kideb, kiut
260         snicswi(ji)   = MAX(0,NINT(dh_snowice(ji)/MAX(epsi20,ABS(dh_snowice(ji)))))
261      ENDDO
262
263      !
264      !------------------------------------------------------------------------------|
265      !  3) Snow redistribution                                                      |
266      !------------------------------------------------------------------------------|
267      !
268      !-------------
269      ! Old profile
270      !-------------
271
272      ! by 'old', it is meant that layers coming from accretion are included,
273      ! and that interfacial layers which were partly melted are reduced
274
275      ! indexes of the vectors
276      !------------------------
277      ntop0    =  1
278      maxnbot0 =  0
279
280      DO ji = kideb, kiut
281         nbot0(ji) =  nlays0  + 1 - snind(ji) + ( 1 - snicind(ji) ) * snicswi(ji)
282         ! cotes of the top of the layers
283         zm0(ji,0) =  0._wp
284         maxnbot0 =  MAX ( maxnbot0 , nbot0(ji) )
285      END DO
286      IF( lk_mpp )   CALL mpp_max( maxnbot0, kcom=ncomm_ice )
287
288      DO jk = 1, maxnbot0
289         DO ji = kideb, kiut
290            !change
291            limsum = ( 1 - snswi(ji) ) * ( jk - 1 ) + snswi(ji) * ( jk + snind(ji) - 1 )
292            limsum = MIN( limsum , nlay_s )
293            zm0(ji,jk) =  dh_s_tot(ji) + zh_s(ji) * REAL( limsum )
294         END DO
295      END DO
296
297      DO ji = kideb, kiut
298         zm0(ji,nbot0(ji)) =  dh_s_tot(ji) - REAL( snicswi(ji) ) * dh_snowice(ji) + zh_s(ji) * REAL( nlays0 )
299         zm0(ji,1)         =  dh_s_tot(ji) * REAL( 1 - snswi(ji) ) + REAL( snswi(ji) ) * zm0(ji,1)
300      END DO
301
302      DO jk = ntop0, maxnbot0
303         DO ji = kideb, kiut
304            zthick0(ji,jk)  =  zm0(ji,jk) - zm0(ji,jk-1)            ! layer thickness
305         END DO
306      END DO
307
308      zqts_in(:) = 0._wp
309
310      DO ji = kideb, kiut         ! layer heat content
311         qm0    (ji,1) =  rhosn * (  cpic * ( rtt - REAL( 1 - snswi(ji) ) * tatm_ice_1d(ji)        &
312            &                                         - REAL( snswi(ji) ) * t_s_b      (ji,1)  )   &
313            &                      + lfus  ) * zthick0(ji,1)
314         zqts_in(ji)   =  zqts_in(ji) + qm0(ji,1) 
315      END DO
316
317      DO jk = 2, maxnbot0
318         DO ji = kideb, kiut
319            limsum      = ( 1 - snswi(ji) ) * ( jk - 1 ) + snswi(ji) * ( jk + snind(ji) - 1 )
320            limsum      = MIN( limsum , nlay_s )
321            qm0(ji,jk)  = rhosn * ( cpic * ( rtt - t_s_b(ji,limsum) ) + lfus ) * zthick0(ji,jk)
322            zswitch = 1.0 - MAX (0.0, SIGN ( 1.0, - ht_s_b(ji) ) )
323            zqts_in(ji) = zqts_in(ji) + REAL( 1 - snswi(ji) ) * qm0(ji,jk) * zswitch
324         END DO ! jk
325      END DO ! ji
326
327      !------------------------------------------------
328      ! Energy given by the snow in snow-ice formation
329      !------------------------------------------------
330      ! zqsnow, enthalpy of the flooded snow
331      DO ji = kideb, kiut
332         zqsnow (ji) =  rhosn * lfus
333         zdeltah(ji) =  0._wp
334      END DO
335
336      DO jk =  nlays0, 1, -1
337         DO ji = kideb, kiut
338            zhsnow =  MAX( 0._wp , dh_snowice(ji)-zdeltah(ji) )
339            zqsnow (ji) =  zqsnow (ji) + rhosn*cpic*(rtt-t_s_b(ji,jk))
340            zdeltah(ji) =  zdeltah(ji) + zh_s(ji)
341         END DO
342      END DO
343
344      DO ji = kideb, kiut
345         zqsnow(ji) = zqsnow(ji) * dh_snowice(ji)
346      END DO
347
348      !------------------
349      ! new snow profile
350      !------------------
351
352      !--------------
353      ! Vector index   
354      !--------------
355      ntop1 =  1
356      nbot1 =  nlay_s
357
358      !-------------------
359      ! Layer coordinates
360      !-------------------
361      DO ji = kideb, kiut
362         zh_s(ji)  = ht_s_b(ji) / REAL( nlay_s )
363         z_s(ji,0) =  0._wp
364      ENDDO
365
366      DO jk = 1, nlay_s
367         DO ji = kideb, kiut
368            z_s(ji,jk) =  zh_s(ji) * REAL( jk )
369         END DO
370      END DO
371
372      !-----------------
373      ! Layer thickness
374      !-----------------
375      DO layer0 = ntop0, maxnbot0
376         DO ji = kideb, kiut
377            zhl0(ji,layer0) = zm0(ji,layer0) - zm0(ji,layer0-1)
378         END DO
379      END DO
380
381      DO layer1 = ntop1, nbot1
382         DO ji = kideb, kiut
383            q_s_b(ji,layer1) = 0._wp
384         END DO
385      END DO
386
387      !----------------
388      ! Weight factors
389      !----------------
390      DO layer0 = ntop0, maxnbot0
391         DO layer1 = ntop1, nbot1
392            DO ji = kideb, kiut
393               zinda = MAX( 0._wp, SIGN( 1._wp , zhl0(ji,layer0) - epsi10 ) )
394               zrl01(layer1,layer0) = zinda * MAX(0.0,( MIN(zm0(ji,layer0),z_s(ji,layer1))   &
395                  &                 - MAX(zm0(ji,layer0-1), z_s(ji,layer1-1))) / MAX(zhl0(ji,layer0),epsi10)) 
396               q_s_b(ji,layer1) = q_s_b(ji,layer1) + zrl01(layer1,layer0)*qm0(ji,layer0)   &
397                  &                                * MAX(0.0,SIGN(1.0,REAL(nbot0(ji)-layer0)))
398            END DO
399         END DO
400      END DO
401
402      ! Heat conservation
403      zqts_fin(:) = 0._wp
404      DO jk = 1, nlay_s
405         DO ji = kideb, kiut
406            zqts_fin(ji) = zqts_fin(ji) + q_s_b(ji,jk)
407         END DO
408      END DO
409
410      IF ( con_i .AND. jiindex_1d > 0 ) THEN
411         DO ji = kideb, kiut
412            IF ( ABS ( zqts_in(ji) - zqts_fin(ji) ) * r1_rdtice  >  1.0e-6 ) THEN
413               ii                 = MOD( npb(ji) - 1, jpi ) + 1
414               ij                 = ( npb(ji) - 1 ) / jpi + 1
415               WRITE(numout,*) ' violation of heat conservation : ', ABS ( zqts_in(ji) - zqts_fin(ji) ) * r1_rdtice
416               WRITE(numout,*) ' ji, jj   : ', ii, ij
417               WRITE(numout,*) ' ht_s_b   : ', ht_s_b(ji)
418               WRITE(numout,*) ' zqts_in  : ', zqts_in (ji) * r1_rdtice
419               WRITE(numout,*) ' zqts_fin : ', zqts_fin(ji) * r1_rdtice
420               WRITE(numout,*) ' dh_snowice : ', dh_snowice(ji)
421               WRITE(numout,*) ' dh_s_tot : ', dh_s_tot(ji)
422               WRITE(numout,*) ' snswi    : ', snswi(ji)
423            ENDIF
424         END DO
425      ENDIF
426
427      !---------------------
428      ! Recover heat content
429      !---------------------
430      DO jk = 1, nlay_s
431         DO ji = kideb, kiut
432            zinda = MAX( 0._wp, SIGN( 1._wp , zh_s(ji) - epsi10 ) )       
433            q_s_b(ji,jk) = zinda * q_s_b(ji,jk) / MAX( zh_s(ji) , epsi10 )
434         END DO !ji
435      END DO !jk 
436
437      !---------------------
438      ! Recover temperature
439      !---------------------
440      zfac1 = 1. / ( rhosn * cpic )
441      zfac2 = lfus / cpic 
442      DO jk = 1, nlay_s
443         DO ji = kideb, kiut
444            zswitch = MAX ( 0.0 , SIGN ( 1.0, - ht_s_b(ji) ) )
445            t_s_b(ji,jk) = rtt + ( 1.0 - zswitch ) * ( - zfac1 * q_s_b(ji,jk) + zfac2 )
446         END DO
447      END DO
448      !
449      !------------------------------------------------------------------------------|
450      !  4) Ice redistribution                                                       |
451      !------------------------------------------------------------------------------|
452      !
453      !-------------
454      ! OLD PROFILE
455      !-------------
456
457      !----------------
458      ! Vector indexes
459      !----------------
460      ntop0    =  1
461      maxnbot0 =  0
462
463      DO ji = kideb, kiut
464         ! reference number of the bottommost layer
465         nbot0(ji) =  MAX( 1 ,  MIN( nlayi0 + ( 1 - icboind(ji) ) +        &
466            &                           ( 1 - icsuind(ji) ) * icsuswi(ji) + snicswi(ji) , nlay_i + 2 ) )
467         ! maximum reference number of the bottommost layer over all domain
468         maxnbot0 =  MAX( maxnbot0 , nbot0(ji) )
469      END DO
470
471      !-------------------------
472      ! Cotes of old ice layers
473      !-------------------------
474      zm0(:,0) =  0._wp
475
476      DO jk = 1, maxnbot0
477         DO ji = kideb, kiut
478            ! jk goes from 1 to nbot0
479            ! the ice layer number goes from 1 to nlay_i
480            ! limsum is the real ice layer number corresponding to present jk
481            limsum    =  ( (icsuswi(ji)*(icsuind(ji)+jk-1) + & 
482               (1-icsuswi(ji))*jk))*(1-snicswi(ji)) + (jk-1)*snicswi(ji)
483            zm0(ji,jk)=  REAL(icsuswi(ji))*dh_i_surf(ji) + REAL(snicswi(ji))*dh_snowice(ji) &
484               +  REAL(limsum) * zh_i(ji)
485         END DO
486      END DO
487
488      DO ji = kideb, kiut
489         zm0(ji,nbot0(ji)) =  REAL(icsuswi(ji))*dh_i_surf(ji) + REAL(snicswi(ji))*dh_snowice(ji) + dh_i_bott(ji) &
490            +  zh_i(ji) * REAL(nlayi0)
491         zm0(ji,1)         =  REAL(snicswi(ji))*dh_snowice(ji) + REAL(1-snicswi(ji))*zm0(ji,1)
492      END DO
493
494      !-----------------------------
495      ! Thickness of old ice layers
496      !-----------------------------
497      DO jk = ntop0, maxnbot0
498         DO ji = kideb, kiut
499            zthick0(ji,jk) =  zm0(ji,jk) - zm0(ji,jk-1)
500         END DO
501      END DO
502
503      !---------------------------
504      ! Inner layers heat content
505      !---------------------------
506      qm0(:,:) =  0.0
507      zqti_in(:) = 0.0
508
509      DO jk = ntop0, maxnbot0
510         DO ji = kideb, kiut
511            limsum =  MAX(1,MIN(snicswi(ji)*(jk-1) + icsuswi(ji)*(jk-1+icsuind(ji)) + &
512               (1-icsuswi(ji))*(1-snicswi(ji))*jk,nlay_i))
513            ztmelts = -tmut * s_i_b(ji,limsum) + rtt
514            qm0(ji,jk) = rhoic * ( cpic * (ztmelts-t_i_b(ji,limsum)) + lfus * ( 1.0-(ztmelts-rtt)/ &
515               MIN((t_i_b(ji,limsum)-rtt),-epsi20) ) - rcp*(ztmelts-rtt) ) &
516               * zthick0(ji,jk)
517         END DO
518      END DO
519
520      !----------------------------
521      ! Bottom layers heat content
522      !----------------------------
523      DO ji = kideb, kiut       
524         ztmelts    = REAL( 1 - icboswi(ji) ) * (-tmut * s_i_b  (ji,nlayi0) )   &   ! case of melting ice
525            &       +     REAL( icboswi(ji) ) * (-tmut * s_i_new(ji)        )   &   ! case of forming ice
526            &       + rtt                                                         ! in Kelvin
527
528         ! bottom formation temperature
529         ztform = t_i_b(ji,nlay_i)
530         IF(  num_sal == 2  )   ztform = t_bo_b(ji)
531         qm0(ji,nbot0(ji)) = REAL( 1 - icboswi(ji) )*qm0(ji,nbot0(ji))             &   ! case of melting ice
532            &              + REAL( icboswi(ji) ) * rhoic * ( cpic*(ztmelts-ztform)       &   ! case of forming ice
533            + lfus *( 1.0-(ztmelts-rtt) / MIN ( (ztform-rtt) , - epsi10 ) )      & 
534            - rcp*(ztmelts-rtt) ) * zthick0(ji,nbot0(ji)  )
535      END DO
536
537      !-----------------------------
538      ! Snow ice layer heat content
539      !-----------------------------
540      DO ji = kideb, kiut
541         ! energy of the flooding seawater
542         zqsnic = rau0 * rcp * ( rtt - t_bo_b(ji) ) * dh_snowice(ji) * &
543            (rhoic - rhosn) / rhoic * REAL(snicswi(ji)) ! generally positive
544         ! Heat conservation diagnostic
545         qt_i_in(ji,jl) = qt_i_in(ji,jl) + zqsnic 
546
547         qldif_1d(ji)   = qldif_1d(ji) + zqsnic * a_i_b(ji)
548
549         ! enthalpy of the newly formed snow-ice layer
550         ! = enthalpy of snow + enthalpy of frozen water
551         zqsnic         =  zqsnow(ji) + zqsnic
552         qm0(ji,1)      =  REAL(snicswi(ji)) * zqsnic + REAL( 1 - snicswi(ji) ) * qm0(ji,1)
553
554      END DO ! ji
555
556      DO jk = ntop0, maxnbot0
557         DO ji = kideb, kiut
558            ! Heat conservation
559            zqti_in(ji) = zqti_in(ji) + qm0(ji,jk) * MAX( 0.0 , SIGN(1.0,ht_i_b(ji)-epsi10) ) &
560               &                                   * MAX( 0.0 , SIGN( 1. , REAL(nbot0(ji) - jk) ) )
561         END DO
562      END DO
563
564      !-------------
565      ! NEW PROFILE
566      !-------------
567
568      !---------------
569      ! Vectors index
570      !---------------
571      ntop1 =  1 
572      nbot1 =  nlay_i
573
574      !------------------
575      ! Layers thickness
576      !------------------
577      DO ji = kideb, kiut
578         zh_i(ji) = ht_i_b(ji) / REAL( nlay_i )
579      ENDDO
580
581      !-------------
582      ! Layer cotes     
583      !-------------
584      z_i(:,0) =  0._wp
585      DO jk = 1, nlay_i
586         DO ji = kideb, kiut
587            z_i(ji,jk) =  zh_i(ji) * jk
588         END DO
589      END DO
590
591      !--thicknesses of the layers
592      DO layer0 = ntop0, maxnbot0
593         DO ji = kideb, kiut
594            zhl0(ji,layer0) = zm0(ji,layer0) - zm0(ji,layer0-1)   ! thicknesses of the layers
595         END DO
596      END DO
597
598      !------------------------
599      ! Weights for relayering
600      !------------------------
601      q_i_b(:,:) = 0._wp
602      DO layer0 = ntop0, maxnbot0
603         DO layer1 = ntop1, nbot1
604            DO ji = kideb, kiut
605               zinda = MAX( 0._wp, SIGN( 1._wp , zhl0(ji,layer0) - epsi10 ) )
606               zrl01(layer1,layer0) = zinda * MAX(0.0,( MIN(zm0(ji,layer0),z_i(ji,layer1)) &
607                  - MAX(zm0(ji,layer0-1), z_i(ji,layer1-1)))/MAX(zhl0(ji,layer0),epsi10))
608               q_i_b(ji,layer1) = q_i_b(ji,layer1) & 
609                  + zrl01(layer1,layer0)*qm0(ji,layer0) &
610                  * MAX(0.0,SIGN(1.0,ht_i_b(ji)-epsi10)) &
611                  * MAX(0.0,SIGN(1.0,REAL(nbot0(ji)-layer0)))
612            END DO
613         END DO
614      END DO
615
616      !-------------------------
617      ! Heat conservation check
618      !-------------------------
619      zqti_fin(:) = 0._wp
620      DO jk = 1, nlay_i
621         DO ji = kideb, kiut
622            zqti_fin(ji) = zqti_fin(ji) + q_i_b(ji,jk)
623         END DO
624      END DO
625      !
626      IF ( con_i .AND. jiindex_1d > 0 ) THEN
627         DO ji = kideb, kiut
628            IF ( ABS ( zqti_in(ji) - zqti_fin(ji) ) * r1_rdtice  >  1.0e-6 ) THEN
629               ii                 = MOD( npb(ji) - 1, jpi ) + 1
630               ij                 = ( npb(ji) - 1 ) / jpi + 1
631               WRITE(numout,*) ' violation of heat conservation : ', ABS ( zqti_in(ji) - zqti_fin(ji) ) * r1_rdtice
632               WRITE(numout,*) ' ji, jj   : ', ii, ij
633               WRITE(numout,*) ' ht_i_b   : ', ht_i_b(ji)
634               WRITE(numout,*) ' zqti_in  : ', zqti_in (ji) * r1_rdtice
635               WRITE(numout,*) ' zqti_fin : ', zqti_fin(ji) * r1_rdtice
636               WRITE(numout,*) ' dh_i_bott: ', dh_i_bott(ji)
637               WRITE(numout,*) ' dh_i_surf: ', dh_i_surf(ji)
638               WRITE(numout,*) ' dh_snowice:', dh_snowice(ji)
639               WRITE(numout,*) ' icsuswi  : ', icsuswi(ji)
640               WRITE(numout,*) ' icboswi  : ', icboswi(ji)
641               WRITE(numout,*) ' snicswi  : ', snicswi(ji)
642            ENDIF
643         END DO
644      ENDIF
645
646      !----------------------
647      ! Recover heat content
648      !----------------------
649      DO jk = 1, nlay_i
650         DO ji = kideb, kiut
651            zinda = MAX( 0._wp, SIGN( 1._wp , zh_i(ji) - epsi10 ) )
652            q_i_b(ji,jk) = zinda * q_i_b(ji,jk) / MAX( zh_i(ji) , epsi10 )
653         END DO !ji
654      END DO !jk 
655
656      ! Heat conservation
657      zqti_fin(:) = 0.0
658      DO jk = 1, nlay_i
659         DO ji = kideb, kiut
660            zqti_fin(ji) = zqti_fin(ji) + q_i_b(ji,jk) * zh_i(ji)
661         END DO
662      END DO
663
664      !
665      !------------------------------------------------------------------------------|
666      !  5) Update salinity and recover temperature                                  |
667      !------------------------------------------------------------------------------|
668      !
669      ! Update salinity (basal entrapment, snow ice formation)
670      DO ji = kideb, kiut
671         sm_i_b(ji) = sm_i_b(ji) + dsm_i_se_1d(ji) + dsm_i_si_1d(ji)
672      END DO !ji
673
674      ! Recover temperature
675      DO jk = 1, nlay_i
676         DO ji = kideb, kiut
677            ztmelts    =  -tmut*s_i_b(ji,jk) + rtt
678            !Conversion q(S,T) -> T (second order equation)
679            zaaa         =  cpic
680            zbbb         =  ( rcp - cpic ) * ( ztmelts - rtt ) + q_i_b(ji,jk) / rhoic - lfus
681            zccc         =  lfus * ( ztmelts - rtt )
682            zdiscrim     =  SQRT( MAX(zbbb*zbbb - 4.0*zaaa*zccc,0.0) )
683            t_i_b(ji,jk) =  rtt - ( zbbb + zdiscrim ) / ( 2.0 *zaaa )
684         END DO !ji
685
686      END DO !jk
687      !
688      CALL wrk_dealloc( jpij, snswi, nbot0, icsuind, icsuswi, icboind, icboswi, snicind, snicswi, snind )   ! integer
689      CALL wrk_dealloc( jpij, zh_i, zh_s, zqsnow, zdeltah, zqti_in, zqts_in, zqti_fin, zqts_fin )           ! real
690      CALL wrk_dealloc( jpij,jkmax+4, zm0, qm0, z_s, z_i, zthick0, zhl0, kjstart = 0 )
691      CALL wrk_dealloc( jkmax+4,jkmax+4, zrl01, kistart = 0, kjstart = 0 )
692      !
693   END SUBROUTINE lim_thd_ent
694
695#else
696   !!----------------------------------------------------------------------
697   !!   Default option                               NO  LIM3 sea-ice model
698   !!----------------------------------------------------------------------
699CONTAINS
700   SUBROUTINE lim_thd_ent          ! Empty routine
701   END SUBROUTINE lim_thd_ent
702#endif
703
704   !!======================================================================
705END MODULE limthd_ent
Note: See TracBrowser for help on using the repository browser.