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

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

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