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

source: trunk/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90 @ 3558

Last change on this file since 3558 was 3558, checked in by rblod, 11 years ago

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