source: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90 @ 4333

Last change on this file since 4333 was 4333, checked in by clem, 7 years ago

remove remaining bugs in LIM3, so that it can run in both regional and global config

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