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

source: branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90 @ 4649

Last change on this file since 4649 was 4649, checked in by clem, 10 years ago

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