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

source: branches/2014/dev_r4743_NOC2_ZTS/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90 @ 4899

Last change on this file since 4899 was 4899, checked in by acc, 9 years ago

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