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

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

Last change on this file since 2528 was 2528, checked in by rblod, 13 years ago

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