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.
limitd_th.F90 in trunk/NEMO/LIM_SRC_3 – NEMO

source: trunk/NEMO/LIM_SRC_3/limitd_th.F90 @ 1264

Last change on this file since 1264 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: 47.8 KB
Line 
1MODULE limitd_th
2#if defined key_lim3
3   !!----------------------------------------------------------------------
4   !!   'key_lim3' :                                   LIM3 sea-ice model
5   !!----------------------------------------------------------------------
6   !!======================================================================
7   !!                       ***  MODULE limitd_th ***
8   !!              Thermodynamics of ice thickness distribution
9   !!                   computation of changes in g(h)     
10   !!======================================================================
11
12   !!----------------------------------------------------------------------
13   !! * Modules used
14   USE dom_ice
15   USE par_oce          ! ocean parameters
16   USE dom_oce
17   USE phycst           ! physical constants (ocean directory)
18   USE ice_oce          ! ice variables
19   USE thd_ice
20   USE limistate
21   USE in_out_manager
22   USE ice
23   USE par_ice
24   USE limthd_lac
25   USE limvar
26   USE iceini
27   USE limcons
28   USE prtctl           ! Print control
29   USE lib_mpp 
30
31   IMPLICIT NONE
32   PRIVATE
33
34   !! * Routine accessibility
35   PUBLIC lim_itd_th        ! called by ice_stp
36   PUBLIC lim_itd_th_rem
37   PUBLIC lim_itd_th_reb
38   PUBLIC lim_itd_fitline
39   PUBLIC lim_itd_shiftice
40
41   !! * Module variables
42   REAL(wp)  ::           &  ! constant values
43      epsi20 = 1e-20   ,  &
44      epsi13 = 1e-13   ,  &
45      zzero  = 0.e0    ,  &
46      zone   = 1.e0
47
48   !!----------------------------------------------------------------------
49   !!   LIM 3.0,  UCL-LOCEAN-IPSL (2008)
50   !! $Id$
51   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
52   !!----------------------------------------------------------------------
53
54
55CONTAINS
56
57   SUBROUTINE lim_itd_th( kt )
58      !!------------------------------------------------------------------
59      !!                ***  ROUTINE lim_itd_th ***
60      !! ** Purpose :
61      !!        This routine computes the thermodynamics of ice thickness
62      !!         distribution
63      !! ** Method  :
64      !!
65      !! ** Arguments :
66      !!           kideb , kiut : Starting and ending points on which the
67      !!                         the computation is applied
68      !!
69      !! ** Inputs / Ouputs : (global commons)
70      !!
71      !! ** External :
72      !!
73      !! ** References :
74      !!
75      !! ** History :
76      !!           (12-2005) Martin Vancoppenolle
77      !!
78      !!------------------------------------------------------------------
79      !! * Arguments
80      INTEGER, INTENT(in) :: kt
81      !! * Local variables
82      INTEGER ::   jl, ja,   &   ! ice category, layers
83         jm,       &   ! ice types    dummy loop index
84         jbnd1,    &
85         jbnd2
86
87      REAL(wp)  ::           &  ! constant values
88         zeps      =  1.0e-10, &
89         epsi10    =  1.0e-10
90
91      IF( kt == nit000 .AND. lwp ) THEN
92         WRITE(numout,*)
93         WRITE(numout,*) 'lim_itd_th  : Thermodynamics of the ice thickness distribution'
94         WRITE(numout,*) '~~~~~~~~~~~'
95      ENDIF
96
97      !------------------------------------------------------------------------------|
98      !  1) Transport of ice between thickness categories.                           |
99      !------------------------------------------------------------------------------|
100      ! Given thermodynamic growth rates, transport ice between
101      ! thickness categories.
102      DO jm = 1, jpm
103         jbnd1 = ice_cat_bounds(jm,1)
104         jbnd2 = ice_cat_bounds(jm,2)
105         IF (ice_ncat_types(jm) .GT. 1 ) CALL lim_itd_th_rem( jbnd1, jbnd2, jm, kt )
106      END DO
107
108      CALL lim_var_glo2eqv ! only for info
109      CALL lim_var_agg(1)
110
111      !------------------------------------------------------------------------------|
112      !  3) Add frazil ice growing in leads.
113      !------------------------------------------------------------------------------|
114
115      CALL lim_thd_lac
116      CALL lim_var_glo2eqv ! only for info
117
118      !----------------------------------------------------------------------------------------
119      !  4) Computation of trend terms and get back to old values     
120      !----------------------------------------------------------------------------------------
121
122      !- Trend terms
123      d_a_i_thd (:,:,:)  = a_i(:,:,:)   - old_a_i(:,:,:) 
124      d_v_s_thd (:,:,:)  = v_s(:,:,:)   - old_v_s(:,:,:)
125      d_v_i_thd (:,:,:)  = v_i(:,:,:)   - old_v_i(:,:,:) 
126      d_e_s_thd(:,:,:,:) = e_s(:,:,:,:) - old_e_s(:,:,:,:) 
127      d_e_i_thd(:,:,:,:) = e_i(:,:,:,:) - old_e_i(:,:,:,:)
128
129      d_smv_i_thd(:,:,:) = 0.0
130      IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) &
131         d_smv_i_thd(:,:,:) = smv_i(:,:,:) - old_smv_i(:,:,:)
132
133      IF(ln_ctl) THEN   ! Control print
134         CALL prt_ctl_info(' ')
135         CALL prt_ctl_info(' - Cell values : ')
136         CALL prt_ctl_info('   ~~~~~~~~~~~~~ ')
137         CALL prt_ctl(tab2d_1=area , clinfo1=' lim_itd_th  : cell area :')
138         CALL prt_ctl(tab2d_1=at_i , clinfo1=' lim_itd_th  : at_i      :')
139         CALL prt_ctl(tab2d_1=vt_i , clinfo1=' lim_itd_th  : vt_i      :')
140         CALL prt_ctl(tab2d_1=vt_s , clinfo1=' lim_itd_th  : vt_s      :')
141         DO jl = 1, jpl
142            CALL prt_ctl_info(' ')
143            CALL prt_ctl_info(' - Category : ', ivar1=jl)
144            CALL prt_ctl_info('   ~~~~~~~~~~')
145            CALL prt_ctl(tab2d_1=a_i   (:,:,jl)   , clinfo1= ' lim_itd_th  : a_i      : ')
146            CALL prt_ctl(tab2d_1=ht_i  (:,:,jl)   , clinfo1= ' lim_itd_th  : ht_i     : ')
147            CALL prt_ctl(tab2d_1=ht_s  (:,:,jl)   , clinfo1= ' lim_itd_th  : ht_s     : ')
148            CALL prt_ctl(tab2d_1=v_i   (:,:,jl)   , clinfo1= ' lim_itd_th  : v_i      : ')
149            CALL prt_ctl(tab2d_1=v_s   (:,:,jl)   , clinfo1= ' lim_itd_th  : v_s      : ')
150            CALL prt_ctl(tab2d_1=e_s   (:,:,1,jl) , clinfo1= ' lim_itd_th  : e_s      : ')
151            CALL prt_ctl(tab2d_1=t_su  (:,:,jl)   , clinfo1= ' lim_itd_th  : t_su     : ')
152            CALL prt_ctl(tab2d_1=t_s   (:,:,1,jl) , clinfo1= ' lim_itd_th  : t_snow   : ')
153            CALL prt_ctl(tab2d_1=sm_i  (:,:,jl)   , clinfo1= ' lim_itd_th  : sm_i     : ')
154            CALL prt_ctl(tab2d_1=smv_i (:,:,jl)   , clinfo1= ' lim_itd_th  : smv_i    : ')
155            DO ja = 1, nlay_i
156               CALL prt_ctl_info(' ')
157               CALL prt_ctl_info(' - Layer : ', ivar1=ja)
158               CALL prt_ctl_info('   ~~~~~~~')
159               CALL prt_ctl(tab2d_1=t_i(:,:,ja,jl) , clinfo1= ' lim_itd_th  : t_i      : ')
160               CALL prt_ctl(tab2d_1=e_i(:,:,ja,jl) , clinfo1= ' lim_itd_th  : e_i      : ')
161            END DO
162         END DO
163      ENDIF
164
165      !- Recover Old values
166      a_i(:,:,:)         = old_a_i (:,:,:)
167      v_s(:,:,:)         = old_v_s (:,:,:)
168      v_i(:,:,:)         = old_v_i (:,:,:)
169      e_s(:,:,:,:)       = old_e_s (:,:,:,:)
170      e_i(:,:,:,:)       = old_e_i (:,:,:,:)
171
172      IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) &
173         smv_i(:,:,:)       = old_smv_i (:,:,:)
174
175   END SUBROUTINE lim_itd_th
176   !
177
178   SUBROUTINE lim_itd_th_rem( klbnd, kubnd, ntyp, kt )
179      !!------------------------------------------------------------------
180      !!                ***  ROUTINE lim_itd_th_rem ***
181      !! ** Purpose :
182      !!        This routine computes the redistribution of ice thickness
183      !!        after thermodynamic growth of ice thickness
184      !!
185      !! ** Method  : Linear remapping
186      !!
187      !! ** Arguments :
188      !!           klbnd, kubnd : Starting and ending category index on which the
189      !!                         the computation is applied
190      !!
191      !! ** Inputs / Ouputs : (global commons)
192      !!
193      !! ** External :
194      !!
195      !! ** References : W.H. Lipscomb, JGR 2001
196      !!
197      !! ** History :
198      !!           largely inspired from CICE (c) W. H. Lipscomb and E.C. Hunke
199      !!
200      !!           (01-2006) Martin Vancoppenolle, UCL-ASTR, translation from
201      !!                     CICE
202      !!           (06-2006) Adaptation to include salt, age and types
203      !!           (04-2007) Mass conservation checked
204      !!------------------------------------------------------------------
205      !! * Arguments
206
207      INTEGER , INTENT (IN) ::  &
208         klbnd ,  &  ! Start thickness category index point
209         kubnd ,  &  ! End point on which the  the computation is applied
210         ntyp  ,  &  ! Number of the type used
211         kt          ! Ocean time step
212
213      !! * Local variables
214      INTEGER ::   ji,       &   ! spatial dummy loop index
215         jj,       &   ! spatial dummy loop index
216         jl,       &   ! ice category dummy loop index
217         zji, zjj, &   ! dummy indices used when changing coordinates
218         nd            ! used for thickness categories
219
220      INTEGER , DIMENSION(jpi,jpj,jpl-1) :: & 
221         zdonor        ! donor category index
222
223      REAL(wp)  ::           &   ! constant values
224         zeps      =  1.0e-10
225
226      REAL(wp)  ::           &  ! constant values for ice enthalpy
227         zindb     ,         &
228         zareamin  ,         &  ! minimum tolerated area in a thickness category
229         zwk1, zwk2,         &  ! all the following are dummy arguments
230         zx1, zx2, zx3,      &  !
231         zetamin   ,         &  ! minimum value of eta
232         zetamax   ,         &  ! maximum value of eta
233         zdh0      ,         &  !
234         zda0      ,         &  !
235         zdamax    ,         &  !
236         zhimin
237
238      REAL(wp), DIMENSION(jpi,jpj,jpl) :: &
239         zdhice           ,  &  ! ice thickness increment
240         g0               ,  &  ! coefficients for fitting the line of the ITD
241         g1               ,  &  ! coefficients for fitting the line of the ITD
242         hL               ,  &  ! left boundary for the ITD for each thickness
243         hR               ,  &  ! left boundary for the ITD for each thickness
244         zht_i_o          ,  &  ! old ice thickness
245         dummy_es
246
247      REAL(wp), DIMENSION(jpi,jpj,jpl-1) :: &
248         zdaice           ,  &  ! local increment of ice area
249         zdvice                 ! local increment of ice volume
250
251      REAL(wp), DIMENSION(jpi,jpj,0:jpl) :: &
252         zhbnew                 ! new boundaries of ice categories
253
254      REAL(wp), DIMENSION(jpi,jpj) :: &
255         zhb0, zhb1             ! category boundaries for thinnes categories
256
257      REAL, DIMENSION(1:(jpi+1)*(jpj+1)) :: &
258         zvetamin, zvetamax     ! maximum values for etas
259
260      INTEGER, DIMENSION(1:(jpi+1)*(jpj+1)) :: &
261         nind_i      ,  &  ! compressed indices for i/j directions
262         nind_j
263
264      INTEGER :: &
265         nbrem             ! number of cells with ice to transfer
266
267      LOGICAL, DIMENSION(jpi,jpj) ::   &  !:
268         zremap_flag             ! compute remapping or not ????
269
270      REAL(wp)  ::           &  ! constant values for ice enthalpy
271         zslope                 ! used to compute local thermodynamic "speeds"
272
273      REAL (wp), DIMENSION(jpi,jpj) :: &  !
274         vt_i_init, vt_i_final,   &  !  ice volume summed over categories
275         vt_s_init, vt_s_final,   &  !  snow volume summed over categories
276         et_i_init, et_i_final,   &  !  ice energy summed over categories
277         et_s_init, et_s_final       !  snow energy summed over categories
278
279      CHARACTER (len = 15) :: fieldid
280
281      !!-- End of declarations
282      !!----------------------------------------------------------------------------------------------
283      zhimin = 0.1      !minimum ice thickness tolerated by the model
284      zareamin = zeps   !minimum area in thickness categories tolerated by the conceptors of the model
285
286      !!----------------------------------------------------------------------------------------------
287      !! 0) Conservation checkand changes in each ice category
288      !!----------------------------------------------------------------------------------------------
289      IF ( con_i ) THEN
290         CALL lim_column_sum (jpl,   v_i, vt_i_init)
291         CALL lim_column_sum (jpl,   v_s, vt_s_init)
292         CALL lim_column_sum_energy (jpl, nlay_i,   e_i, et_i_init)
293         dummy_es(:,:,:) = e_s(:,:,1,:)
294         CALL lim_column_sum (jpl, dummy_es(:,:,:) , et_s_init)
295      ENDIF
296
297      !!----------------------------------------------------------------------------------------------
298      !! 1) Compute thickness and changes in each ice category
299      !!----------------------------------------------------------------------------------------------
300      IF (kt == nit000 .AND. lwp) THEN
301         WRITE(numout,*)
302         WRITE(numout,*) 'lim_itd_th_rem  : Remapping the ice thickness distribution'
303         WRITE(numout,*) '~~~~~~~~~~~~~~~'
304         WRITE(numout,*) ' klbnd :       ', klbnd
305         WRITE(numout,*) ' kubnd :       ', kubnd
306         WRITE(numout,*) ' ntyp  :       ', ntyp 
307      ENDIF
308
309      zdhice(:,:,:) = 0.0
310      DO jl = klbnd, kubnd
311         DO jj = 1, jpj
312            DO ji = 1, jpi
313               zindb             = 1.0-MAX(0.0,SIGN(1.0,-a_i(ji,jj,jl)))     !0 if no ice and 1 if yes
314               ht_i(ji,jj,jl)    = v_i(ji,jj,jl) / MAX(a_i(ji,jj,jl),zeps) * zindb
315               zindb             = 1.0-MAX(0.0,SIGN(1.0,-old_a_i(ji,jj,jl))) !0 if no ice and 1 if yes
316               zht_i_o(ji,jj,jl) = old_v_i(ji,jj,jl) / MAX(old_a_i(ji,jj,jl),zeps) * zindb
317               IF (a_i(ji,jj,jl).gt.1e-6) THEN
318                  zdhice(ji,jj,jl) = ht_i(ji,jj,jl) - zht_i_o(ji,jj,jl) 
319               ENDIF
320            END DO
321         END DO
322      END DO
323
324      !-----------------------------------------------------------------------------------------------
325      !  2) Compute fractional ice area in each grid cell
326      !-----------------------------------------------------------------------------------------------
327      at_i(:,:) = 0.0
328      DO jl = klbnd, kubnd
329         DO jj = 1, jpj
330            DO ji = 1, jpi
331               at_i(ji,jj) = at_i(ji,jj) + a_i(ji,jj,jl)
332            END DO
333         END DO
334      END DO
335
336      !-----------------------------------------------------------------------------------------------
337      !  3) Identify grid cells with ice
338      !-----------------------------------------------------------------------------------------------
339      nbrem = 0
340      DO jj = 1, jpj
341         DO ji = 1, jpi
342            IF ( at_i(ji,jj) .gt. zareamin ) THEN
343               nbrem         = nbrem + 1
344               nind_i(nbrem) = ji
345               nind_j(nbrem) = jj
346               zremap_flag(ji,jj) = .true.
347            ELSE
348               zremap_flag(ji,jj) = .false.
349            ENDIF
350         END DO !ji
351      END DO !jj
352
353      !-----------------------------------------------------------------------------------------------
354      !  4) Compute new category boundaries
355      !-----------------------------------------------------------------------------------------------
356      !- 4.1 Compute category boundaries
357      ! Tricky trick see limitd_me.F90
358      ! will be soon removed, CT
359      ! hi_max(kubnd) = 999.99
360      zhbnew(:,:,:) = 0.0
361
362      DO jl = klbnd, kubnd - 1
363         ! jl
364         DO ji = 1, nbrem
365            ! jl, ji
366            zji = nind_i(ji)
367            zjj = nind_j(ji)
368            !
369            IF ( ( zht_i_o(zji,zjj,jl)  .GT.zeps ) .AND. & 
370               ( zht_i_o(zji,zjj,jl+1).GT.zeps ) ) THEN
371               !interpolate between adjacent category growth rates
372               zslope = ( zdhice(zji,zjj,jl+1)     - zdhice(zji,zjj,jl) ) / &
373                  ( zht_i_o   (zji,zjj,jl+1) - zht_i_o   (zji,zjj,jl) )
374               zhbnew(zji,zjj,jl) = hi_max(jl) + zdhice(zji,zjj,jl) + &
375                  zslope * ( hi_max(jl) - zht_i_o(zji,zjj,jl) )
376            ELSEIF (zht_i_o(zji,zjj,jl).gt.zeps) THEN
377               zhbnew(zji,zjj,jl) = hi_max(jl) + zdhice(zji,zjj,jl)
378            ELSEIF (zht_i_o(zji,zjj,jl+1).gt.zeps) THEN
379               zhbnew(zji,zjj,jl) = hi_max(jl) + zdhice(zji,zjj,jl+1)
380            ELSE
381               zhbnew(zji,zjj,jl) = hi_max(jl)
382            ENDIF
383            ! jl, ji
384         END DO !ji
385         ! jl
386
387         !- 4.2 Check that each zhbnew lies between adjacent values of ice thickness
388         DO ji = 1, nbrem
389            ! jl, ji
390            zji = nind_i(ji)
391            zjj = nind_j(ji)
392            ! jl, ji
393            IF ( ( a_i(zji,zjj,jl) .GT.zeps) .AND. & 
394               ( ht_i(zji,zjj,jl).GE. zhbnew(zji,zjj,jl) ) &
395               ) THEN
396               zremap_flag(zji,zjj) = .false.
397            ELSEIF ( ( a_i(zji,zjj,jl+1) .GT. zeps ) .AND. &
398               ( ht_i(zji,zjj,jl+1).LE. zhbnew(zji,zjj,jl) ) &
399               ) THEN
400               zremap_flag(zji,zjj) = .false.
401            ENDIF
402
403            !- 4.3 Check that each zhbnew does not exceed maximal values hi_max 
404            ! jl, ji
405            IF (zhbnew(zji,zjj,jl).gt.hi_max(jl+1)) THEN
406               zremap_flag(zji,zjj) = .false.
407            ENDIF
408            ! jl, ji
409            IF (zhbnew(zji,zjj,jl).lt.hi_max(jl-1)) THEN
410               zremap_flag(zji,zjj) = .false.
411            ENDIF
412            ! jl, ji
413         END DO !ji
414         ! ji
415      END DO !jl
416
417      !-----------------------------------------------------------------------------------------------
418      !  5) Identify cells where ITD is to be remapped
419      !-----------------------------------------------------------------------------------------------
420      nbrem = 0
421      DO jj = 1, jpj
422         DO ji = 1, jpi
423            IF ( zremap_flag(ji,jj) ) THEN
424               nbrem         = nbrem + 1
425               nind_i(nbrem) = ji
426               nind_j(nbrem) = jj
427            ENDIF
428         END DO !ji
429      END DO !jj
430
431      !-----------------------------------------------------------------------------------------------
432      !  6) Fill arrays with lowermost / uppermost boundaries of 'new' categories
433      !-----------------------------------------------------------------------------------------------
434      DO jj = 1, jpj
435         DO ji = 1, jpi
436            zhb0(ji,jj) = hi_max_typ(0,ntyp) ! 0eme
437            zhb1(ji,jj) = hi_max_typ(1,ntyp) ! 1er
438
439            zhbnew(ji,jj,klbnd-1) = 0.0
440
441            IF ( a_i(ji,jj,kubnd) .GT. zeps ) THEN
442               zhbnew(ji,jj,kubnd) = 3.0*ht_i(ji,jj,kubnd) - 2.0*zhbnew(ji,jj,kubnd-1)
443            ELSE
444               zhbnew(ji,jj,kubnd) = hi_max(kubnd)
445            ENDIF
446
447            IF ( zhbnew(ji,jj,kubnd) .LT. hi_max(kubnd-1) ) &
448               zhbnew(ji,jj,kubnd) = hi_max(kubnd-1)
449
450         END DO !jj
451      END DO !jj
452
453      !-----------------------------------------------------------------------------------------------
454      !  7) Compute g(h)
455      !-----------------------------------------------------------------------------------------------
456      !- 7.1 g(h) for category 1 at start of time step
457      CALL lim_itd_fitline(klbnd, zhb0, zhb1, zht_i_o(:,:,klbnd), &
458         g0(:,:,klbnd), g1(:,:,klbnd), hL(:,:,klbnd), &
459         hR(:,:,klbnd), zremap_flag)
460
461      !- 7.2 Area lost due to melting of thin ice (first category,  klbnd)
462      DO ji = 1, nbrem
463         zji = nind_i(ji) 
464         zjj = nind_j(ji) 
465
466         !ji
467         IF (a_i(zji,zjj,klbnd) .gt. zeps) THEN
468            zdh0 = zdhice(zji,zjj,klbnd) !decrease of ice thickness in the lower category
469            ! ji, a_i > zeps
470            IF (zdh0 .lt. 0.0) THEN !remove area from category 1
471               ! ji, a_i > zeps; zdh0 < 0
472               zdh0 = MIN(-zdh0,hi_max(klbnd))
473
474               !Integrate g(1) from 0 to dh0 to estimate area melted
475               zetamax = MIN(zdh0,hR(zji,zjj,klbnd)) - hL(zji,zjj,klbnd)
476               IF (zetamax.gt.0.0) THEN
477                  zx1  = zetamax
478                  zx2  = 0.5 * zetamax*zetamax 
479                  zda0 = g1(zji,zjj,klbnd) * zx2 + g0(zji,zjj,klbnd) * zx1 !ice area removed
480                  ! Constrain new thickness <= ht_i
481                  zdamax = a_i(zji,zjj,klbnd) * & 
482                     (1.0 - ht_i(zji,zjj,klbnd)/zht_i_o(zji,zjj,klbnd)) ! zdamax > 0
483                  !ice area lost due to melting of thin ice
484                  zda0   = MIN(zda0, zdamax)
485
486                  ! Remove area, conserving volume
487                  ht_i(zji,zjj,klbnd) = ht_i(zji,zjj,klbnd) & 
488                     * a_i(zji,zjj,klbnd) / ( a_i(zji,zjj,klbnd) - zda0 )
489                  a_i(zji,zjj,klbnd)  = a_i(zji,zjj,klbnd) - zda0
490                  v_i(zji,zjj,klbnd)  = a_i(zji,zjj,klbnd)*ht_i(zji,zjj,klbnd)
491               ENDIF     ! zetamax > 0
492               ! ji, a_i > zeps
493
494            ELSE ! if ice accretion
495               ! ji, a_i > zeps; zdh0 > 0
496               IF ( ntyp .EQ. 1 ) zhbnew(zji,zjj,klbnd-1) = MIN(zdh0,hi_max(klbnd)) 
497               ! zhbnew was 0, and is shifted to the right to account for thin ice
498               ! growth in openwater (F0 = f1)
499               IF ( ntyp .NE. 1 ) zhbnew(zji,zjj,0) = 0 
500               ! in other types there is
501               ! no open water growth (F0 = 0)
502            ENDIF ! zdh0
503
504            ! a_i > zeps
505         ENDIF ! a_i > zeps
506
507      END DO ! ji
508
509      !- 7.3 g(h) for each thickness category 
510      DO jl = klbnd, kubnd
511         CALL lim_itd_fitline(jl, zhbnew(:,:,jl-1), zhbnew(:,:,jl), ht_i(:,:,jl), &
512            g0(:,:,jl), g1(:,:,jl), hL(:,:,jl), hR(:,:,jl),     &
513            zremap_flag)
514      END DO
515
516      !-----------------------------------------------------------------------------------------------
517      !  8) Compute area and volume to be shifted across each boundary
518      !-----------------------------------------------------------------------------------------------
519
520      DO jl = klbnd, kubnd - 1
521         DO jj = 1, jpj
522            DO ji = 1, jpi
523               zdonor(ji,jj,jl) = 0
524               zdaice(ji,jj,jl) = 0.0
525               zdvice(ji,jj,jl) = 0.0
526            END DO
527         END DO
528
529         DO ji = 1, nbrem
530            zji = nind_i(ji)
531            zjj = nind_j(ji)
532
533            IF (zhbnew(zji,zjj,jl) .gt. hi_max(jl)) THEN ! transfer from jl to jl+1
534
535               ! left and right integration limits in eta space
536               zvetamin(ji) = MAX(hi_max(jl), hL(zji,zjj,jl)) - hL(zji,zjj,jl)
537               zvetamax(ji) = MIN(zhbnew(zji,zjj,jl), hR(zji,zjj,jl)) - hL(zji,zjj,jl)
538               zdonor(zji,zjj,jl) = jl
539
540            ELSE  ! zhbnew(jl) <= hi_max(jl) ; transfer from jl+1 to jl
541
542               ! left and right integration limits in eta space
543               zvetamin(ji) = 0.0
544               zvetamax(ji) = MIN(hi_max(jl), hR(zji,zjj,jl+1)) - hL(zji,zjj,jl+1)
545               zdonor(zji,zjj,jl) = jl + 1
546
547            ENDIF  ! zhbnew(jl) > hi_max(jl)
548
549            zetamax = MAX(zvetamax(ji), zvetamin(ji)) ! no transfer if etamax < etamin
550            zetamin = zvetamin(ji)
551
552            zx1  = zetamax - zetamin
553            zwk1 = zetamin*zetamin
554            zwk2 = zetamax*zetamax
555            zx2  = 0.5 * (zwk2 - zwk1)
556            zwk1 = zwk1 * zetamin
557            zwk2 = zwk2 * zetamax
558            zx3  = 1.0/3.0 * (zwk2 - zwk1)
559            nd   = zdonor(zji,zjj,jl)
560            zdaice(zji,zjj,jl) = g1(zji,zjj,nd)*zx2 + g0(zji,zjj,nd)*zx1
561            zdvice(zji,zjj,jl) = g1(zji,zjj,nd)*zx3 + g0(zji,zjj,nd)*zx2 + &
562               zdaice(zji,zjj,jl)*hL(zji,zjj,nd)
563
564         END DO ! ji
565      END DO ! jl klbnd -> kubnd - 1
566
567      !!----------------------------------------------------------------------------------------------
568      !! 9) Shift ice between categories
569      !!----------------------------------------------------------------------------------------------
570      CALL lim_itd_shiftice ( klbnd, kubnd, zdonor, zdaice, zdvice )
571
572      !!----------------------------------------------------------------------------------------------
573      !! 10) Make sure ht_i >= minimum ice thickness hi_min
574      !!----------------------------------------------------------------------------------------------
575
576      DO ji = 1, nbrem
577         zji = nind_i(ji)
578         zjj = nind_j(ji)
579         IF ( ( zhimin .GT. 0.0 ) .AND. & 
580            ( ( a_i(zji,zjj,1) .GT. zeps ) .AND. ( ht_i(zji,zjj,1) .LT. zhimin ) ) &
581            ) THEN
582            a_i(zji,zjj,1)  = a_i(zji,zjj,1) * ht_i(zji,zjj,1) / zhimin 
583            ht_i(zji,zjj,1) = zhimin
584            v_i(zji,zjj,1)  = a_i(zji,zjj,1)*ht_i(zji,zjj,1)
585         ENDIF
586      END DO !ji
587
588      !!----------------------------------------------------------------------------------------------
589      !! 11) Conservation check
590      !!----------------------------------------------------------------------------------------------
591      IF ( con_i ) THEN
592         CALL lim_column_sum (jpl,   v_i, vt_i_final)
593         fieldid = ' v_i : limitd_th '
594         CALL lim_cons_check (vt_i_init, vt_i_final, 1.0e-6, fieldid) 
595
596         CALL lim_column_sum_energy (jpl, nlay_i,  e_i, et_i_final)
597         fieldid = ' e_i : limitd_th '
598         CALL lim_cons_check (et_i_init, et_i_final, 1.0e-3, fieldid) 
599
600         CALL lim_column_sum (jpl,   v_s, vt_s_final)
601         fieldid = ' v_s : limitd_th '
602         CALL lim_cons_check (vt_s_init, vt_s_final, 1.0e-6, fieldid) 
603
604         dummy_es(:,:,:) = e_s(:,:,1,:)
605         CALL lim_column_sum (jpl, dummy_es(:,:,:) , et_s_final)
606         fieldid = ' e_s : limitd_th '
607         CALL lim_cons_check (et_s_init, et_s_final, 1.0e-3, fieldid) 
608      ENDIF
609
610   END SUBROUTINE lim_itd_th_rem
611   !
612
613   SUBROUTINE lim_itd_fitline(num_cat, HbL, Hbr, hice, g0, g1, hL, hR, zremap_flag )
614
615      !!------------------------------------------------------------------
616      !!                ***  ROUTINE lim_itd_fitline ***
617      !! ** Purpose :
618      !! fit g(h) with a line using area, volume constraints
619      !!
620      !! ** Method  :
621      !! Fit g(h) with a line, satisfying area and volume constraints.
622      !! To reduce roundoff errors caused by large values of g0 and g1,
623      !! we actually compute g(eta), where eta = h - hL, and hL is the
624      !! left boundary.
625      !!
626      !! ** Arguments :
627      !!
628      !! ** Inputs / Ouputs : (global commons)
629      !!
630      !! ** External :
631      !!
632      !! ** References :
633      !!
634      !! ** History :
635      !! authors: William H. Lipscomb, LANL, Elizabeth C. Hunke, LANL
636      !!          (01-2006) Martin Vancoppenolle
637      !!
638      !!------------------------------------------------------------------
639      !! * Arguments
640
641      INTEGER, INTENT(in) :: num_cat      ! category index
642
643      REAL(wp), DIMENSION(jpi,jpj), INTENT(IN)   ::   &  !:
644         HbL, HbR        ! left and right category boundaries
645
646      REAL(wp), DIMENSION(jpi,jpj), INTENT(IN)   ::   &  !:
647         hice            ! ice thickness
648
649      REAL(wp), DIMENSION(jpi,jpj), INTENT(OUT)  ::   &  !:
650         g0, g1      , & ! coefficients in linear equation for g(eta)
651         hL          , & ! min value of range over which g(h) > 0
652         hR              ! max value of range over which g(h) > 0
653
654      LOGICAL, DIMENSION(jpi,jpj), INTENT(IN)    ::   &  !:
655         zremap_flag
656
657      INTEGER :: &             
658         ji,jj           ! horizontal indices
659
660      REAL(wp) :: &           
661         zh13        , & ! HbL + 1/3 * (HbR - HbL)
662         zh23        , & ! HbL + 2/3 * (HbR - HbL)
663         zdhr        , & ! 1 / (hR - hL)
664         zwk1, zwk2  , & ! temporary variables
665         zacrith         ! critical minimum concentration in an ice category
666
667      REAL(wp)  ::           &  ! constant values
668         zeps      =  1.0e-10
669
670      zacrith       = 1.0e-6
671      !!-- End of declarations
672      !!----------------------------------------------------------------------------------------------
673
674      DO jj = 1, jpj
675         DO ji = 1, jpi
676
677            IF ( zremap_flag(ji,jj) .AND. a_i(ji,jj,num_cat) .gt. zacrith &
678               .AND. hice(ji,jj) .GT. 0.0 ) THEN
679
680               ! Initialize hL and hR
681
682               hL(ji,jj) = HbL(ji,jj)
683               hR(ji,jj) = HbR(ji,jj)
684
685               ! Change hL or hR if hice falls outside central third of range
686
687               zh13 = 1.0/3.0 * (2.0*hL(ji,jj) + hR(ji,jj))
688               zh23 = 1.0/3.0 * (hL(ji,jj) + 2.0*hR(ji,jj))
689
690               IF (hice(ji,jj) < zh13) THEN
691                  hR(ji,jj) = 3.0*hice(ji,jj) - 2.0*hL(ji,jj)
692               ELSEIF (hice(ji,jj) > zh23) THEN
693                  hL(ji,jj) = 3.0*hice(ji,jj) - 2.0*hR(ji,jj)
694               ENDIF
695
696               ! Compute coefficients of g(eta) = g0 + g1*eta
697
698               zdhr = 1.0 / (hR(ji,jj) - hL(ji,jj))
699               zwk1 = 6.0 * a_i(ji,jj,num_cat) * zdhr
700               zwk2 = (hice(ji,jj) - hL(ji,jj)) * zdhr
701               g0(ji,jj) = zwk1 * (2.0/3.0 - zwk2)
702               g1(ji,jj) = 2.0*zdhr * zwk1 * (zwk2 - 0.5)
703
704            ELSE                   ! remap_flag = .false. or a_i < zeps
705
706               hL(ji,jj) = 0.0
707               hR(ji,jj) = 0.0
708               g0(ji,jj) = 0.0
709               g1(ji,jj) = 0.0
710
711            ENDIF                  ! a_i > zeps
712
713         END DO !ji
714      END DO ! jj
715
716   END SUBROUTINE lim_itd_fitline
717   !
718
719   SUBROUTINE lim_itd_shiftice (klbnd, kubnd, zdonor, zdaice, zdvice)
720      !!------------------------------------------------------------------
721      !!                ***  ROUTINE lim_itd_shiftice ***
722      !! ** Purpose : shift ice across category boundaries, conserving everything
723      !!              ( area, volume, energy, age*vol, and mass of salt )
724      !!
725      !! ** Method  :
726      !!
727      !! ** Arguments :
728      !!
729      !! ** Inputs / Ouputs : (global commons)
730      !!
731      !! ** External :
732      !!
733      !! ** References :
734      !!
735      !! ** History :
736      !! authors: William H. Lipscomb, LANL, Elizabeth C. Hunke, LANL
737      !!          (01-2006) Martin Vancoppenolle
738      !!
739      !!------------------------------------------------------------------
740      !! * Arguments
741
742      INTEGER , INTENT (IN) ::  &
743         klbnd ,  &  ! Start thickness category index point
744         kubnd       ! End point on which the  the computation is applied
745
746      INTEGER , DIMENSION(jpi,jpj,jpl-1), INTENT(IN) :: & 
747         zdonor             ! donor category index
748
749      REAL(wp), DIMENSION(jpi,jpj,jpl-1), INTENT(INOUT) :: & 
750         zdaice     ,  &   ! ice area transferred across boundary
751         zdvice            ! ice volume transferred across boundary
752
753      INTEGER :: &
754         ji,jj,jl,      &  ! horizontal indices, thickness category index
755         jl2,           &  ! receiver category
756         jl1,           &  ! donor category
757         jk,            &  ! ice layer index
758         zji, zjj          ! indices when changing from 2D-1D is done
759
760      REAL(wp), DIMENSION(jpi,jpj,jpl) :: &
761         zaTsfn
762
763      REAL(wp), DIMENSION(jpi,jpj) :: &
764         zworka            ! temporary array used here
765
766      REAL(wp) :: &         
767         zdvsnow     ,  &  ! snow volume transferred
768         zdesnow     ,  &  ! snow energy transferred
769         zdeice      ,  &  ! ice energy transferred
770         zdsm_vice      ,  &  ! ice salinity times volume transferred
771         zdo_aice      ,  &  ! ice age times volume transferred
772         zdaTsf      ,  &  ! aicen*Tsfcn transferred
773         zindsn      ,  &  ! snow or not
774         zindb             ! ice or not
775
776      INTEGER, DIMENSION(1:(jpi+1)*(jpj+1)) :: &
777         nind_i      ,  &  ! compressed indices for i/j directions
778         nind_j
779
780      INTEGER :: &
781         nbrem             ! number of cells with ice to transfer
782
783      LOGICAL :: &
784         zdaice_negative       , & ! true if daice < -puny
785         zdvice_negative       , & ! true if dvice < -puny
786         zdaice_greater_aicen  , & ! true if daice > aicen
787         zdvice_greater_vicen      ! true if dvice > vicen
788
789      REAL(wp)  ::           &  ! constant values
790         zeps      =  1.0e-10
791
792      !!-- End of declarations
793
794      !----------------------------------------------------------------------------------------------
795      ! 1) Define a variable equal to a_i*T_su
796      !----------------------------------------------------------------------------------------------
797
798      DO jl = klbnd, kubnd
799         DO jj = 1, jpj
800            DO ji = 1, jpi
801               zaTsfn(ji,jj,jl) = a_i(ji,jj,jl)*t_su(ji,jj,jl)
802            END DO ! ji
803         END DO ! jj
804      END DO ! jl
805
806      !----------------------------------------------------------------------------------------------
807      ! 2) Check for daice or dvice out of range, allowing for roundoff error
808      !----------------------------------------------------------------------------------------------
809      ! Note: zdaice < 0 or zdvice < 0 usually happens when category jl
810      ! has a small area, with h(n) very close to a boundary.  Then
811      ! the coefficients of g(h) are large, and the computed daice and
812      ! dvice can be in error. If this happens, it is best to transfer
813      ! either the entire category or nothing at all, depending on which
814      ! side of the boundary hice(n) lies.
815      !-----------------------------------------------------------------
816      DO jl = klbnd, kubnd-1
817
818         zdaice_negative = .false.
819         zdvice_negative = .false.
820         zdaice_greater_aicen = .false.
821         zdvice_greater_vicen = .false.
822
823         DO jj = 1, jpj
824            DO ji = 1, jpi
825
826               IF (zdonor(ji,jj,jl) .GT. 0) THEN
827                  jl1 = zdonor(ji,jj,jl)
828
829                  IF (zdaice(ji,jj,jl) .LT. 0.0) THEN
830                     IF (zdaice(ji,jj,jl) .GT. -zeps) THEN
831                        IF ( ( jl1.EQ.jl   .AND. ht_i(ji,jj,jl1) .GT. hi_max(jl) )           &
832                           .OR.                                      &
833                           ( jl1.EQ.jl+1 .AND. ht_i(ji,jj,jl1) .LE. hi_max(jl) )           & 
834                           ) THEN                                                             
835                           zdaice(ji,jj,jl) = a_i(ji,jj,jl1)  ! shift entire category
836                           zdvice(ji,jj,jl) = v_i(ji,jj,jl1)
837                        ELSE
838                           zdaice(ji,jj,jl) = 0.0 ! shift no ice
839                           zdvice(ji,jj,jl) = 0.0
840                        ENDIF
841                     ELSE
842                        zdaice_negative = .true.
843                     ENDIF
844                  ENDIF
845
846                  IF (zdvice(ji,jj,jl) .LT. 0.0) THEN
847                     IF (zdvice(ji,jj,jl) .GT. -zeps ) THEN
848                        IF ( ( jl1.EQ.jl .AND. ht_i(ji,jj,jl1).GT.hi_max(jl) )     &
849                           .OR.                                     &
850                           ( jl1.EQ.jl+1 .AND. ht_i(ji,jj,jl1) .LE. hi_max(jl) ) &
851                           ) THEN
852                           zdaice(ji,jj,jl) = a_i(ji,jj,jl1) ! shift entire category
853                           zdvice(ji,jj,jl) = v_i(ji,jj,jl1) 
854                        ELSE
855                           zdaice(ji,jj,jl) = 0.0    ! shift no ice
856                           zdvice(ji,jj,jl) = 0.0
857                        ENDIF
858                     ELSE
859                        zdvice_negative = .true.
860                     ENDIF
861                  ENDIF
862
863                  ! If daice is close to aicen, set daice = aicen.
864                  IF (zdaice(ji,jj,jl) .GT. a_i(ji,jj,jl1) - zeps ) THEN
865                     IF (zdaice(ji,jj,jl) .LT. a_i(ji,jj,jl1)+zeps) THEN
866                        zdaice(ji,jj,jl) = a_i(ji,jj,jl1)
867                        zdvice(ji,jj,jl) = v_i(ji,jj,jl1) 
868                     ELSE
869                        zdaice_greater_aicen = .true.
870                     ENDIF
871                  ENDIF
872
873                  IF (zdvice(ji,jj,jl) .GT. v_i(ji,jj,jl1)-zeps) THEN
874                     IF (zdvice(ji,jj,jl) .LT. v_i(ji,jj,jl1)+zeps) THEN
875                        zdaice(ji,jj,jl) = a_i(ji,jj,jl1)
876                        zdvice(ji,jj,jl) = v_i(ji,jj,jl1) 
877                     ELSE
878                        zdvice_greater_vicen = .true.
879                     ENDIF
880                  ENDIF
881
882               ENDIF               ! donor > 0
883            END DO                   ! i
884         END DO                 ! j
885
886      END DO !jl
887
888      !-------------------------------------------------------------------------------
889      ! 3) Transfer volume and energy between categories
890      !-------------------------------------------------------------------------------
891
892      DO jl = klbnd, kubnd - 1
893         nbrem = 0
894         DO jj = 1, jpj
895            DO ji = 1, jpi
896               IF (zdaice(ji,jj,jl) .GT. 0.0 ) THEN ! daice(n) can be < puny
897                  nbrem = nbrem + 1
898                  nind_i(nbrem) = ji
899                  nind_j(nbrem) = jj
900               ENDIF ! tmask
901            END DO
902         END DO
903
904         DO ji = 1, nbrem 
905            zji = nind_i(ji)
906            zjj = nind_j(ji)
907
908            jl1 = zdonor(zji,zjj,jl)
909            zindb             = MAX( 0.0 , SIGN( 1.0 , v_i(zji,zjj,jl1) - zeps ) )
910            zworka(zji,zjj)   = zdvice(zji,zjj,jl) / MAX(v_i(zji,zjj,jl1),zeps) * zindb
911            IF (jl1 .eq. jl) THEN
912               jl2 = jl1+1
913            ELSE                ! n1 = n+1
914               jl2 = jl 
915            ENDIF
916
917            !--------------
918            ! Ice areas
919            !--------------
920
921            a_i(zji,zjj,jl1) = a_i(zji,zjj,jl1) - zdaice(zji,zjj,jl)
922            a_i(zji,zjj,jl2) = a_i(zji,zjj,jl2) + zdaice(zji,zjj,jl)
923
924            !--------------
925            ! Ice volumes
926            !--------------
927
928            v_i(zji,zjj,jl1) = v_i(zji,zjj,jl1) - zdvice(zji,zjj,jl) 
929            v_i(zji,zjj,jl2) = v_i(zji,zjj,jl2) + zdvice(zji,zjj,jl)
930
931            !--------------
932            ! Snow volumes
933            !--------------
934
935            zdvsnow          = v_s(zji,zjj,jl1) * zworka(zji,zjj)
936            v_s(zji,zjj,jl1) = v_s(zji,zjj,jl1) - zdvsnow
937            v_s(zji,zjj,jl2) = v_s(zji,zjj,jl2) + zdvsnow 
938
939            !--------------------
940            ! Snow heat content 
941            !--------------------
942
943            zdesnow              = e_s(zji,zjj,1,jl1) * zworka(zji,zjj)
944            e_s(zji,zjj,1,jl1)   = e_s(zji,zjj,1,jl1) - zdesnow
945            e_s(zji,zjj,1,jl2)   = e_s(zji,zjj,1,jl2) + zdesnow
946
947            !--------------
948            ! Ice age
949            !--------------
950
951            zdo_aice             = oa_i(zji,zjj,jl1) * zdaice(zji,zjj,jl)
952            oa_i(zji,zjj,jl1)    = oa_i(zji,zjj,jl1) - zdo_aice
953            oa_i(zji,zjj,jl2)    = oa_i(zji,zjj,jl2) + zdo_aice
954
955            !--------------
956            ! Ice salinity
957            !--------------
958
959            zdsm_vice            = smv_i(zji,zjj,jl1) * zworka(zji,zjj)
960            smv_i(zji,zjj,jl1)   = smv_i(zji,zjj,jl1) - zdsm_vice
961            smv_i(zji,zjj,jl2)   = smv_i(zji,zjj,jl2) + zdsm_vice
962
963            !---------------------
964            ! Surface temperature
965            !---------------------
966
967            zdaTsf               = t_su(zji,zjj,jl1) * zdaice(zji,zjj,jl)
968            zaTsfn(zji,zjj,jl1)  = zaTsfn(zji,zjj,jl1) - zdaTsf
969            zaTsfn(zji,zjj,jl2)  = zaTsfn(zji,zjj,jl2) + zdaTsf 
970
971         END DO                 ! ji
972
973         !------------------
974         ! Ice heat content
975         !------------------
976
977         DO jk = 1, nlay_i
978!CDIR NODEP
979            DO ji = 1, nbrem
980               zji = nind_i(ji)
981               zjj = nind_j(ji)
982
983               jl1 = zdonor(zji,zjj,jl)
984               IF (jl1 .EQ. jl) THEN
985                  jl2 = jl+1
986               ELSE             ! n1 = n+1
987                  jl2 = jl 
988               ENDIF
989
990               zdeice = e_i(zji,zjj,jk,jl1) * zworka(zji,zjj)
991               e_i(zji,zjj,jk,jl1) =  e_i(zji,zjj,jk,jl1) - zdeice
992               e_i(zji,zjj,jk,jl2) =  e_i(zji,zjj,jk,jl2) + zdeice 
993            END DO              ! ji
994         END DO                 ! jk
995
996      END DO                   ! boundaries, 1 to ncat-1
997
998      !-----------------------------------------------------------------
999      ! Update ice thickness and temperature
1000      !-----------------------------------------------------------------
1001
1002      DO jl = klbnd, kubnd
1003         DO jj = 1, jpj
1004            DO ji = 1, jpi 
1005               IF ( a_i(ji,jj,jl) .GT. zeps ) THEN
1006                  ht_i(ji,jj,jl)  =  v_i(ji,jj,jl) / a_i(ji,jj,jl) 
1007                  t_su(ji,jj,jl)  =  zaTsfn(ji,jj,jl) / a_i(ji,jj,jl) 
1008                  zindsn          =  1.0 - MAX(0.0,SIGN(1.0,-v_s(ji,jj,jl))) !0 if no ice and 1 if yes
1009               ELSE
1010                  ht_i(ji,jj,jl)  = 0.0
1011                  t_su(ji,jj,jl)  = rtt
1012               ENDIF
1013            END DO                 ! ji
1014         END DO                 ! jj
1015      END DO                    ! jl
1016
1017   END SUBROUTINE lim_itd_shiftice
1018   !
1019
1020   SUBROUTINE lim_itd_th_reb(klbnd, kubnd, ntyp)
1021      !!------------------------------------------------------------------
1022      !!                ***  ROUTINE lim_itd_th_reb ***
1023      !! ** Purpose : rebin - rebins thicknesses into defined categories
1024      !!
1025      !! ** Method  :
1026      !!
1027      !! ** Arguments :
1028      !!
1029      !! ** Inputs / Ouputs : (global commons)
1030      !!
1031      !! ** External :
1032      !!
1033      !! ** References :
1034      !!
1035      !! ** History : (2005) Translation from CICE
1036      !!              (2006) Adaptation to include salt, age and types
1037      !!              (2007) Mass conservation checked
1038      !!
1039      !! authors: William H. Lipscomb, LANL, Elizabeth C. Hunke, LANL
1040      !!          (01-2006) Martin Vancoppenolle (adaptation)
1041      !!
1042      !!------------------------------------------------------------------
1043      !! * Arguments
1044      INTEGER , INTENT (in) ::  &
1045         klbnd ,  &  ! Start thickness category index point
1046         kubnd ,  &  ! End point on which the  the computation is applied
1047         ntyp        ! number of the ice type involved in the rebinning process
1048
1049      INTEGER :: &
1050         ji,jj,          &  ! horizontal indices
1051         jl                 ! category index
1052
1053      INTEGER ::   &  !:
1054         zshiftflag          ! = .true. if ice must be shifted
1055
1056      INTEGER, DIMENSION(jpi,jpj,jpl) :: &
1057         zdonor             ! donor category index
1058
1059      REAL(wp), DIMENSION(jpi, jpj, jpl) :: &
1060         zdaice         , & ! ice area transferred
1061         zdvice             ! ice volume transferred
1062
1063      REAL(wp)  ::           &  ! constant values
1064         zeps      =  1.0e-10, &
1065         epsi10    =  1.0e-10
1066
1067      REAL (wp), DIMENSION(jpi,jpj) :: &  !
1068         vt_i_init, vt_i_final,   &  !  ice volume summed over categories
1069         vt_s_init, vt_s_final       !  snow volume summed over categories
1070
1071      CHARACTER (len = 15) :: fieldid
1072
1073      !!-- End of declarations
1074      !------------------------------------------------------------------------------
1075
1076      !     ! conservation check
1077      IF ( con_i ) THEN
1078         CALL lim_column_sum (jpl,   v_i, vt_i_init)
1079         CALL lim_column_sum (jpl,   v_s, vt_s_init)
1080      ENDIF
1081
1082      !
1083      !------------------------------------------------------------------------------
1084      ! 1) Compute ice thickness.
1085      !------------------------------------------------------------------------------
1086      DO jl = klbnd, kubnd
1087         DO jj = 1, jpj
1088            DO ji = 1, jpi 
1089               IF (a_i(ji,jj,jl) .GT. zeps) THEN
1090                  ht_i(ji,jj,jl) = v_i(ji,jj,jl) / a_i(ji,jj,jl)
1091               ELSE
1092                  ht_i(ji,jj,jl) = 0.0
1093               ENDIF
1094            END DO                 ! i
1095         END DO                 ! j
1096      END DO                    ! n
1097
1098      !------------------------------------------------------------------------------
1099      ! 2) Make sure thickness of cat klbnd is at least hi_max_typ(klbnd)
1100      !------------------------------------------------------------------------------
1101      DO jj = 1, jpj 
1102         DO ji = 1, jpi 
1103
1104            IF (a_i(ji,jj,klbnd) > zeps) THEN
1105               IF (ht_i(ji,jj,klbnd) .LE. hi_max_typ(0,ntyp) .AND. hi_max_typ(0,ntyp) .GT. 0.0 ) THEN
1106                  a_i(ji,jj,klbnd)  = v_i(ji,jj,klbnd) / hi_max_typ(0,ntyp) 
1107                  ht_i(ji,jj,klbnd) = hi_max_typ(0,ntyp)
1108               ENDIF
1109            ENDIF
1110         END DO                    ! i
1111      END DO                    ! j
1112
1113      !------------------------------------------------------------------------------
1114      ! 3) If a category thickness is not in bounds, shift the
1115      ! entire area, volume, and energy to the neighboring category
1116      !------------------------------------------------------------------------------
1117      !-------------------------
1118      ! Initialize shift arrays
1119      !-------------------------
1120
1121      DO jl = klbnd, kubnd
1122         DO jj = 1, jpj 
1123            DO ji = 1, jpi
1124               zdonor(ji,jj,jl) = 0
1125               zdaice(ji,jj,jl) = 0.0
1126               zdvice(ji,jj,jl) = 0.0
1127            END DO
1128         END DO
1129      END DO
1130
1131      !-------------------------
1132      ! Move thin categories up
1133      !-------------------------
1134
1135      DO jl = klbnd, kubnd - 1  ! loop over category boundaries
1136
1137         !---------------------------------------
1138         ! identify thicknesses that are too big
1139         !---------------------------------------
1140         zshiftflag = 0
1141
1142         DO jj = 1, jpj 
1143            DO ji = 1, jpi 
1144               IF (a_i(ji,jj,jl) .GT. zeps .AND. ht_i(ji,jj,jl) .GT. hi_max(jl) ) THEN
1145                  zshiftflag        = 1
1146                  zdonor(ji,jj,jl)  = jl 
1147                  zdaice(ji,jj,jl)  = a_i(ji,jj,jl)
1148                  zdvice(ji,jj,jl)  = v_i(ji,jj,jl)
1149               ENDIF
1150            END DO                 ! ji
1151         END DO                 ! jj
1152         IF( lk_mpp ) CALL mpp_max(zshiftflag)
1153
1154         IF ( zshiftflag == 1 ) THEN
1155
1156            !------------------------------
1157            ! Shift ice between categories
1158            !------------------------------
1159            CALL lim_itd_shiftice (klbnd, kubnd, zdonor, zdaice, zdvice)
1160
1161            !------------------------
1162            ! Reset shift parameters
1163            !------------------------
1164            DO jj = 1, jpj
1165               DO ji = 1, jpi
1166                  zdonor(ji,jj,jl) = 0
1167                  zdaice(ji,jj,jl) = 0.0
1168                  zdvice(ji,jj,jl) = 0.0
1169               END DO
1170            END DO
1171
1172         ENDIF                  ! zshiftflag
1173
1174      END DO                    ! jl
1175
1176      !----------------------------
1177      ! Move thick categories down
1178      !----------------------------
1179
1180      DO jl = kubnd - 1, 1, -1       ! loop over category boundaries
1181
1182         !-----------------------------------------
1183         ! Identify thicknesses that are too small
1184         !-----------------------------------------
1185         zshiftflag = 0
1186
1187         DO jj = 1, jpj
1188            DO ji = 1, jpi
1189               IF (a_i(ji,jj,jl+1) .GT. zeps .AND. &
1190                  ht_i(ji,jj,jl+1) .LE. hi_max(jl)) THEN
1191
1192                  zshiftflag = 1
1193                  zdonor(ji,jj,jl) = jl + 1
1194                  zdaice(ji,jj,jl) = a_i(ji,jj,jl+1) 
1195                  zdvice(ji,jj,jl) = v_i(ji,jj,jl+1)
1196               ENDIF
1197            END DO                 ! ji
1198         END DO                 ! jj
1199
1200         IF(lk_mpp) CALL mpp_max(zshiftflag)
1201         IF (zshiftflag==1) THEN
1202
1203            !------------------------------
1204            ! Shift ice between categories
1205            !------------------------------
1206            CALL lim_itd_shiftice (klbnd, kubnd, zdonor, zdaice, zdvice)
1207
1208            !------------------------
1209            ! Reset shift parameters
1210            !------------------------
1211            DO jj = 1, jpj 
1212               DO ji = 1, jpi 
1213                  zdonor(ji,jj,jl)  = 0
1214                  zdaice(ji,jj,jl)  = 0.0
1215                  zdvice(ji,jj,jl)  = 0.0
1216               END DO
1217            END DO
1218
1219         ENDIF                  ! zshiftflag
1220
1221      END DO                    ! jl
1222
1223      !------------------------------------------------------------------------------
1224      ! 4) Conservation check
1225      !------------------------------------------------------------------------------
1226
1227      IF ( con_i ) THEN
1228         CALL lim_column_sum (jpl,   v_i, vt_i_final)
1229         fieldid = ' v_i : limitd_reb '
1230         CALL lim_cons_check (vt_i_init, vt_i_final, 1.0e-6, fieldid) 
1231
1232         CALL lim_column_sum (jpl,   v_s, vt_s_final)
1233         fieldid = ' v_s : limitd_reb '
1234         CALL lim_cons_check (vt_s_init, vt_s_final, 1.0e-6, fieldid) 
1235      ENDIF
1236
1237   END SUBROUTINE lim_itd_th_reb
1238
1239#else
1240   !!======================================================================
1241   !!                       ***  MODULE limitd_th    ***
1242   !!                              no sea ice model
1243   !!======================================================================
1244CONTAINS
1245   SUBROUTINE lim_itd_th           ! Empty routines
1246   END SUBROUTINE lim_itd_th
1247   SUBROUTINE lim_itd_th_ini
1248   END SUBROUTINE lim_itd_th_ini
1249   SUBROUTINE lim_itd_th_rem
1250   END SUBROUTINE lim_itd_th_rem
1251   SUBROUTINE lim_itd_fitline
1252   END SUBROUTINE lim_itd_fitline
1253   SUBROUTINE lim_itd_shiftice
1254   END SUBROUTINE lim_itd_shiftice
1255   SUBROUTINE lim_itd_th_reb
1256   END SUBROUTINE lim_itd_th_reb
1257#endif
1258END MODULE limitd_th
Note: See TracBrowser for help on using the repository browser.