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

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

LIM3: solve ticket #1421

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