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

source: tags/nemo_v3_2_2/NEMO/LIM_SRC_3/limitd_th.F90 @ 2720

Last change on this file since 2720 was 2477, checked in by cetlod, 13 years ago

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