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 @ 921

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

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