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/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90 @ 5055

Last change on this file since 5055 was 5055, checked in by clem, 9 years ago

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