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

Last change on this file since 4990 was 4990, checked in by timgraham, 6 years ago

Merged branches/2014/dev_MERGE_2014 back onto the trunk as follows:

In the working copy of branch ran:
svn merge svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk@HEAD
1 conflict in LIM_SRC_3/limdiahsb.F90
Resolved by keeping the version from dev_MERGE_2014 branch
and commited at r4989

In working copy run:
svn switch svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk
to switch working copy

Run:
svn merge —reintegrate svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/branches/2014/dev_MERGE_2014
to merge the branch into the trunk - no conflicts at this stage.

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