New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
limitd_th.F90 in trunk/NEMO/LIM_SRC_3 – NEMO

source: trunk/NEMO/LIM_SRC_3/limitd_th.F90 @ 844

Last change on this file since 844 was 844, checked in by ctlod, 16 years ago

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