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 @ 863

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

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