source: trunk/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90 @ 5202

Last change on this file since 5202 was 5202, checked in by clem, 6 years ago

LIM3: important bug fix to avoid crashes. See ticket #1508

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