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

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

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