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

source: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limthd_ent.F90 @ 4333

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

remove remaining bugs in LIM3, so that it can run in both regional and global 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.