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

source: tags/nemo_v3_2/nemo_v3_2/NEMO/LIM_SRC_3/limthd_ent.F90 @ 1878

Last change on this file since 1878 was 1878, checked in by flavoni, 14 years ago

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