source: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90 @ 4293

Last change on this file since 4293 was 4293, checked in by clem, 7 years ago

small corrections for ice categories hi_max and for shifting ice between categories

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