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

Last change on this file since 2542 was 2528, checked in by rblod, 13 years ago

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