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

source: branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90 @ 3517

Last change on this file since 3517 was 3517, checked in by gm, 11 years ago

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