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

source: tags/nemo_v3_2/nemo_v3_2/NEMO/LIM_SRC_3/limitd_th.F90 @ 1878

Last change on this file since 1878 was 1878, checked in by flavoni, 14 years ago

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