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

source: trunk/NEMO/LIM_SRC_3/limthd_ent.F90 @ 1156

Last change on this file since 1156 was 1156, checked in by rblod, 16 years ago

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