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

source: trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_ent.F90 @ 2760

Last change on this file since 2760 was 2715, checked in by rblod, 13 years ago

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