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

Last change on this file since 1146 was 1146, checked in by rblod, 16 years ago

Add svn Id (first try), see ticket #210

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