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

source: branches/dev_r4028_CNRS_LIM3_MV2014/NEMOGCM/NEMO/LIM_SRC_3/limthd_ent.F90 @ 4506

Last change on this file since 4506 was 4506, checked in by vancop, 10 years ago

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