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/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90 @ 3963

Last change on this file since 3963 was 3963, checked in by clem, 11 years ago

bugs correction + creation of glob_max and glob_min in lib_fortran.F90, see ticket:#1116

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